The p-version and hp-version of the finite element method have now been widely accepted as efficient, accurate and flexible tools for analysing many linear problems in computational mechanics. It is yet still unclear if and how this behaviour can also be observed for non-linear problems. Only few publications have addressed problems of elasto-plasticity until now and theoretical results on the expected rate of convergence or an optimal choice of the element size and p-degree are even more rare. Results of the p-version for two different non-linear material models for small displacements and small strains will be presented in the following, demonstrating that the p-version allows for an efficient treatment of elasto-plastic problems.
 

The deformation theory of plasticity

The deformation theory of plasticity, first proposed by Hencky, is valid for an isotropic material under radial loading up to strains lower than the damage threshold. The main difference between the deformation theory and the flow theory of Prandlt-Reuss is the computation of the plastic strains. In flow theory the plastic strains are given by a rate equation. Due to the postulated proportional loading in the case of the the deformation theory the flow rule, describing the change of plastic strains, can be formally integrated, leading to an algebraic equation for determining the plastic strains. Due to the lack of an internal time the deformation theory is not able to memorize the loading history, which plays a fundamental role in plasticity.

To get an impression of the accuracy and efficiency of high order elements for the deformation theory of plasticity we consider the problem being defined by Stein as benchmark of the German research project Adaptive finite element methods in applied mechanics Figure 1 shows a quarter of a square plate (h=b=100) with central hole (R=10) and unit thickness (t=1), loaded by uniform traction of magnitude 450. The material is assumed to be elastic-perfectly-plastic and plane-strain conditions are considered. A von Mises yield criterion is applied and the following material parameters are assumed: shear modulus=80193.8, bulk modulus=164206.0 and yield stress=450.0. Results of interest are the displacements at point 4 and point 5.
 

Three different meshes with 2, 4 and 10 quadrilateral elements taking advantage of the blending function method are chosen to mesh the plate. A series of computations with polynomial degrees up to p=17 was carried out.

In order to draw a comparison to an adaptive h-version we refer to the results of Barthold, Schmidt and Stein. The computations there were performed with the Q1-P0 element. A mesh consisting of 64 Q1-P0 elements was refined in 10 steps using the equilibrium criterion by BABUSKA and MILLER, yielding 875 elements. The results of the p-version are also to be compared to a sequence of graded meshes with Q1-P0 elements obtained by Barthold, Schmidt and Stein.

Comparing the results of the uniform p-version with those of the h-version based on a sequence of grades meshes, we observe that the efficiency of the p-version is superior (see Figures 4,5). The mesh with 10 elements supplies with a polynomial degree of p=5 and 540 degrees of freedom an accuracy which is higher than the one of the adaptive h-version. The corresponding computational time on a workstation with an alpha processor with 500 MHz amounts to only 1 second. 

In Figure 6 the von Mises stress is plotted for three different discretizations: the mesh with 2 elements and p=17, the mesh with 4 elements and p=9 and the mesh with 10 elements and p=9. Obviously, the different approximations show only small deviations from each other. 

The flow theory of plasticity

A physically more realistic model problem for elasto-plastic behaviour is the so-called flow theory proposed by Prandtl and Reuss. Here we consider a von Mises yield criterion accounting for non-linear isotropic hardening. It is assumed that strains and displacements are small. To test the p-version for this material model we consider the three-dimensional problem being defined by Stein as benchmark of the German research project Adaptive finite element methods in applied mechanics . 

A thick-walled plate with circular hole under monotonous load is discretized with 48 hexahedral elements. Due to symmetry only an eighth of the plate has to be discretized. The curved boundary of the hole was taken care of with the blending function method. The load was raised monotonously within 61 load steps up to a factor of 4.15.
 
 
 

The first result to be considered is the displacement component ux at point 2. As a reference solution we use the results of WIENERS having been obtained with a fine mesh consisting of 1048576 Q1-P0 hexahedral elements resulting in an equation system of 3368499 unkowns. Comparing the results for p>=4 with the reference solution, only very small deviations are visible.

In Figure 10 the stress component syy at point 2 is plotted. The maximum of the relative error for the approximation based on p=8 is 1.1% being very accurate for a three-dimensional non-linear computation.

The significance of a three-dimensional computation is proven in Figure 11 where the von Mises stress and the Gaussian points where yielding occurs are plotted for the last load step, showing the variation of the plastic region over the plate thickness.

A complex example

Finally, a last example is considered in order to demonstrate that the high order finite element approach may also be applied to more complex problems. Figure 12 depicts a structural joint which is part of a cupola. 

Figure 12: Structural joint discretized with 162 hexahedral elements

The structural joint made of steel is composed of rods which are merged with a cylindrical shell. In order to increase the stiffness of the structural joint a sheet plate is integrated concentrically. The material is described by the J2 flow theory, including nonlinear isotropic hardening. To reduce the numerical effort, only one eighth of the structural joint is discretized and symmetry boundary conditions are introduced to excluded the rigid body motion mode. The system is subjected to an axial force caused by a uniform tension acting on the free end of the rod. Furthermore, it is assumed that-- due to some inaccuracies during the installation of the cupola - a small bending moment is applied to the structural joint. This bending moment is modelled by a uniform pressure of acting normal to the axis of the rod. A classical finite element approach would call for special elements in order to model the transition from shell to solid elements. However, this approach would entail all the difficulties which are due to the coupling of elements of different dimensions. In particular, the transition from the shell to solid part of the structural joint, where stress concentrations are expected, would be a source of significant errors when dimensionally reduced models are applied. The structural joint is discretized with 162 high order hexahedral elements, as depicted in Figure 12. To resolve the singularities, the mesh is refined at reentrant corners and edges.



                                

Figure 13: Convergence of displacement

The trunk space with  p = p_x = p_h = p_z = 1,...,7, is applied to perform a series of computations where the whole load history is recomputed  for each polynomial degree. To get an impression of the accuracy of the p-extension the results for p=1,2,3,...,7 are plotted for points A,B,C,D (see Figure 13). From this it is obvious that the finite element solution with p = 7 yields a reliable approximation.  The evolution of plastic zone is sketched in Figure 14 for three different load steps with l = 2.5, l = 3.45, and l = 3.6. In Figure 15 the convergence of the von Mises stress at point E is plotted, demonstrating the accuracy of the p-extension.
 

Figure 14: Evolution of the plastic zone

Figure 15: Convergence of the von Mises stress at point E

Related movies