Pergamon 0045-7949(94)0031 l-4
Con~puler.s & Yrucrur
L. Bernspang et al 28
Fig. 1. Geometry of cantilever beam.
plastic-partly elastic elements is presented. The shape calculated from the derivatives of the displacements. of the plastic zone is continuous and independent In (4) matrix B includes the derivatives of the basis from the element borders. functions and + cr* is the sum of the converged
stress from the previous step (without index i) and the
PROPOSED ALGORITHM stress increment from the previous iteration in the
When solving a plasticity problem, the solution is current step:
usually found by applying the external load in (time-) steos. In each such step, eouilibrium is found by
,I+ $,*f = VC* + ,I+ I&*G 1, (5)
iterations. If the solution is known for step n, the next solution n + 1 is found from the unbalanced force The unbalanced forces will give displacement incre-
vector g in iteration i ments + A II for step 12 + 1 from
+ g = BT + lb *r da _ )I + I& (4) (6)
where K is a stiffness matrix set up in the conven- using the consistent stresses u* (a vector with the tional way of the standard displacement method, stress components) instead of the standard stresses except for the plastic elements where the Hessian can
Fig. 2. Mesh and plastic part with 512 elements for consistent and conventional analysis.
Consistent strain method in finite element plasticity
Table 1. Top deflection for the different meshes with con- sistent and conventional analysis
Elements Consistent analysis Conventional analysis
8 11.52 3.82 32 16.12 8.85
128 16.19 13.35 512 16.22 15.40
be evaluated by making use of all the points where stresses are projected. The new solution is found as
II + I = + U + IA , (7)
From this solution, increments of element strains are calculated Fig. 4. Geometry of thick tube.
These are used in (3) to obtain the consistent strain values from
n= (n+Afi*-+Aci)2dQ-+min. s
With the nodal strain approximation + 'ACHE= @n+h*i this gives
r+ IA;*= M-1 s
@T+lAcidQ (10) R
where M is the standard mass matrix. (Inversion of the mass matrix can be avoided by lumping or by use of a Jacobi iteration that solves the equation on an element level at comparatively low costs.)
The strain increments +A?* are used to evaluate nodal stress increments n+AB*i with the con-
stitutive relation. The evaluation is done in each node instead of in each integration point used in the conventional analysis. In all realistic problems there are more integration points than nodes, so the storage requirement is reduced by this method. The nodal stress increments n+At?*r are used in (5) to evaluate the consistent nodal stresses and in (4) to evaluate the unbalanced force vector so that the local equilibrium condition will be fulfilled with the continuous stress approxima- tion at the end of the iterations. The necessary iterative process (even in the linear elastic case) is legitimized as higher accuracy is obtained, but the most powerful characteristic feature of the method is in the case of plasticity as the continuous stresses evaluated in the proposed way become plastically admissible.
If all nodes of element is elastic,
an element are elastic the entire but if one or more nodes have a
o- 11 1 1 0 20 40 60 80 100 120 140 160 180 200
Elements Fig. 3. Top deflection for four different meshes.
L. BernspAng 61 ul.
Fig. 5. Mesh and plastic part with 135 elements (consistent and conventional analysis)
plastic stress, the element is examined internally. The elastic trial stress of a node is then found from
where D,. is the elastic moduli matrix. The stresses are then checked against the yield surface . If all nodal stresses in an element are elastic, they will be pro- portional to the strains throughout the element, and the approximations of strains and stresses are of equal accuracy. However, if one or more nodes are plastic, the stress will not be proportional to the strain over the element and in order to obtain stresses with the same order of accuracy as the strains the following procedure can be applied. It is possible to evaluate the plastically admissible stress from
Here Ai is the increment of a plastic multiplier and n is the unit normal vector to the yield surface, giving the plastic strain increment I+ At *IJ = Ain, as AC*=A~*-A~*.
These stresses are used to integrate the unbalanced force vector
+ k = s B+ la*&) _ U+ fL-? (13) 0 as well as to evaluate the plastic stiffness matrix of the element. If an element is purely elastic or purely plastic the integration of nodal values and of these
internal values, will give identical results, and the elastic or plastic stiffness matrix can be evaluated without knowledge of the state variables of the element.
For the case of partly elastic-partly plastic el- ements, a separated integration of the elastic part and the plastic part is recommended. The simplest way is to increase the number of integration points even for low order elements. The increase of state variables is limited by the restricted number of partly elas- tic-partly plastic elements. The total number of state variables will remain smaller than in the conventional analysis due to the fact that for the fully plastic part of the domain only the nodal values have to be stored. The plastic zone can now be described inde- pendently of the element borders.
The evaluation of the consistent nodal strains as given in (IO) is a global operation even if it can be done in an iterative process on an element level. It is not possible to use the given relations between con- sistent stresses and strains to obtain an element stiffness matrix that would give convergence in the first iteration of a linear elastic problem. The neces- sity to iterate will increase the cost considerably for purely elastic problems, but for plastic problems iterations cannot be avoided. so there the difference will be less pronounced.
This algorithm is especially suited for low order elements like linear triangles or bilinear quadrilaterals due to the constant derivatives of their displacement functions.
Fig. 6. Mesh and plastic part with 540 elements (consistent and conventional analysis).
Consistent strain method in finite element plasticity 31
In these numerical examples, three-node elements with linear basis functions have been used. For the evaluation of plastically admissible stresses inside partly elastic-partly plastic eIements, 12 integration points have been used. The mesh refinement is ob- tained by dividing each element into four new el- ements of equal size and shape.
Example I. CantiIever beam with moment load
In this example a cantilever beam with dimensions as shown in Fig. 1 is subjected to a moment load at its free end. The material is elasto-plastic with hard- ening, so that a strain level of three times the yield strain corresponds to a stress that is twice the yield stress. Poissons ratio is 0.3.
For this example no exact solution is known. If simple beam theory is applied, the maximum elastic
1.0 1 0.9
0. t , , , , , , , ,
0.0 0.2 0.4 0.6 0.X 1.0 I.2 1.4 I.6
Distance from outer surface Consistent analysis
(b) ,.I -
g 0.8 -
5 0.7 - id ._ n/I - z 3 0.5 -
B 0.4 - P
0.0 I 1 , , , I 1 I 1 0 0 0.2 11.4 A 0 H 1. I.2 1.4 I.6 LX
Distance from wter surface Conventir~nal analysis
Fig. 7. Tangential stress along the vertical boundary with 135 elements.
radius on the
stresses in the clamped end will be the yield stress. One knows, however, that close to the top and bottom sides of the beam, the stresses will be parallel to the side and not to the main axis as assumed in beam theory. From Fig. 2 one finds that the computer solution seems to be very close to such a solution giving plastic zones limited by straight lines from the clamped end.
Four different meshes have been tested, where a new mesh is obtained from the previous mesh by dividing each element into four new elements. The first mesh has eight elements and in Fig. 2 the final mesh of 512 elements is shown. The top deflections for the different meshes are compared in Table 1 and Fig. 3 where results from both consistent and conven- tional analysis are shown.
From Table 1 it seems as if the top deflection tends towards a value just over 16.22, but as stated above, no exact solution is known. Note that the mesh with
(a> 1.1 -
2 a.7 - B * -3 0.6 -
g 0.5 -
1 0.4 -
no d.o.* Distance from outer surface Consistent analysis
I I I 1 , , 1 I I 0.0 (I.2 0.4 f1.6 0.x I.0 I.2 1.4 I.6 1.X
Distance from outer surface Conventiond analysis
Fig. 8. Tangential stress along the vertical radius on the boundary with 540 elements.
32 L. Bernsping ef al.
32 elements and consistent stresses gives a better solution for the top deflection than does the mesh with 512 elements and conventional stresses.
Example 2. Thick-walled tube
In  the theoretical solution of a thick-walled tube (see Fig. 4) in an elastic-perfectly plastic material and Trescas yield criterion is given. The tube is subjected to an internal pressure, causing plastic strains from the inner surface to a radius r = c, indicated with ticks on the two radii in Figs 5 and 6. The figures are organized as in the previous example-used mesh, internal points that are plastic from the new algor- ithm and plastic elements from a conventional analy- sis with constant strain triangles. It is shown that the size of the plastic zone is independent of the element borders and thus the result is less sensitive to different mesh designs.
Distance trom outer surfilce Distance from outer sutike Conventional analysir Conventional analysis
Fig. 9. Radial stress along the vertical radius on the bound- ary with 135 elements.
Fig. 10. Radial stress along the vertical radius on the boundary with 540 elements.
The formulas given in  can be rewritten for the Huber-von Mises criterion giving the radial stresses in the plastic inner part
and the elastic stresses in the outer part
The radius c of the plastic zone can be determined from the condition
at the inner radius r = a.
Consistent strain method in finite element plasticity 33
The corresponding tangential stresses are found from
As the exact solution is known, the exact and the calculated stress distributions can be compared (Figs 7-10). The exact stress distribution is given as a dashed line. With the new algorithm the nodal values 6* are plotted, whereas for the conventional calculation the nodal averages are plotted. Results from the new algorithm are given in parts (a), and those from a conventional analysis in parts (b). Note the improvements in the peak value, as well as the values on the boundaries.
In the algorithm presented here, continuity over element borders is reached by using the nodal values of the stresses and strains in the actual calculations. The local constitutive equation is fulfilled as well as the equilibrium condition in a weak sense at the nodes. The results obtained with the new algorithm are less sensitive to different mesh designs than those obtained from the conventional method. The predic- tion of the plastic zone is already improved on a very coarse mesh and is independent from the element borders even for low order elements. The computer time needed increases compared with that for the
same mesh, but then we find a considerable increase in accuracy. This algorithm is also useful for raising the accuracy without remeshing the problem.
0. C. Zienkiewicz and R. L. Zhu, A simple error estimator and adaptive procedure for practical engin- eering analysis. Inr. J. Numer. Meth. Engng 2.4, 337-357 (1987). N.-E. Wiberg and F. Abdulwahab, An efficient postpro- cessing technique for stress problems based on super- convergent derivatives and equilibrium. In Proc. First European Conference on Numbical Methods in Engin- eering 1992, pp. 25-32. Elsevier Science, Amsterdam (1992). L. Bernspiing, Iterative and adaptive solution tech- niques in computational plasticity. Publ. 91:8, Depart- ment of Structural Mechanics, Chalmers University of Technology, Giiteborg, Sweden (1991). L. Bemsping, Solution techniques for linear and non- linear finite element equations. Publ. 87:3, Department of Structural Mechanics, Chalmers University of Tech- nology, Goteborg, Sweden (1987). G. Cantin, G. Loubignac and G. Touzot, An iterative algorithm to build continuous stress and displacement solutions. Int. J. Numer. Merh. Engng 12, 1493-1506 (1978). J. T. Oden and H. J. Brauchli, On the calculation of consistent stress distributions in finite element approxi- mations. Inr. J. Numer. Mefh. Engng 3, 317-325 (1971). 0. C. Zienkiewicz, L. Kui and S. Nakazawa, Iterative solution of mixed problems and the stress recovery procedure. Commun. ad. Numer. Meth. 1. 3-9 (1985). k. L. Taylor and J:C. Simo, A return mapping algorithm for plane stress elasto-plasticity. Int. J. Nu- mer. Melh. Engng 22, 649670 (1986). R. Hill, The Mathematical Theory of Plasticity. Clarendon Press, Oxford (1950).