- Email: [email protected]
A mixed finite element model for plane strain elastic—plastic analysis Part II. Application to the 4-node bilinear element
Computer methods in applied mechanics and englneerlng
ElSEVlER
Comput.
Methods
Appl. Mech. Engrg.
141 (1997) 81-93
A mixed finite element model for plane strain elastic-plastic analysis Part II. Application to the 4-node bilinear element A. Capsoni, L. Corradi” Department
oj Structural
Engineering,
Politecnico
Received
di Milano,
3 1 March
Piazza Leonardo
1995; revised
da Vinci 32, Milano,
I-20133,
Italy
18 April 1996
Abstract In this paper, the results of part I of this study are applied to the 4-node bilinear element. The trial functions for different fields consistently with the a-priori conditions established, the rate problem is generalized to cover finite loading steps and the solution presented. Computations referred to some discriminating examples from the literature demonstrate that the proposed model capable of predicting the perfectly plastic collapse load for pathological plane strain problems, exhibits excellent performances meshes and is rather insensitive to mesh distortion.
are defined strategy is is actually with coarse
1. Introduction The first part of this study [l] was devoted to a mixed finite element formulation of the plane strain problem for elastically isotrolpic solids with a plastic behavior governed by von Mises condition. It was shown that if the trial functions for different fields comply with some conditions, the resulting discrete scheme is able to prevent locking when plasticity spreads, while preserving its stress redistribution capability. The main purpose of this contribution is to demonstrate, with reference to the 4-node bilinear element, how such a model can be actually established and used for computations. General results are specialized to the element considered .and suitable approximations for all the variables involved are presented; the rate problem is then generalized to cover finite loading steps, as required by computations; the solution strategy is next presented and, finally, the efficiency of the formulation is .demonstrated by solving some problems which are known as pathological. None of the above steps contains contributions that are actually new. The models for different fields are essentially ones that have been already proposed [2,3], with minor modifications introduced to meet the a-priori requirements. The transition from rates to finite increments rests on a standard backward difference integration of the rate law, and ,aclassical step-wise holonomic problem [4] is obtained as a result. The solution algorithm is a strain-driven return mapping one, adapting existing procedures to the problem at hand [5,6]. The model proves capable of correctly predicting perfectly plastic collapse, is highly accurate in coarse meshes and rather insensitive to shape distortion, but its performances are shared by other formulations. Nevertheless, it is felt that the proposal deserves some attention for at least two reasons. First, the problem is formulated in terms of generalized variables in Prager’s sense, which permits the use of
*Corresponding
author.
00457825197/$17.00 0 1997 Elsevier Science S.A. All rights reserved PII SOO45-7825(96)01099-7
82
A. Capsoni, L. Corradi 1 Comput. Methods Appl. Mech. Engrg. 141 (1997) 81-93
results available in this context. In particular, some devices resting on Mathematical Programming theory [ [I], Ref. 221 can be exploited at different stages of the solution process to overcome numerical difficulties. Secondly, the constitutive law is fully reduced to the plane; except that, for the expressions of matrices governing the material behavior, the elastoplastic plane strain problem is driven to a form identical to a plane stress one and, as in the elastic case, both problems can be dealt with in the same manner, with significant conceptual and computational simplifications.
2. Definition of field approximations The isoparametric 4-node element is considered. Actual and intrinsic coordinates and x = {x y}‘, respectively. Nodal coordinates are collected in the vectors X, = 1X, X, X, X4)‘; and the following
are denoted by X = {X Y}’
(1)
Y, = {Y, Y, Y3 YqK
definitions
are introduced
(2a-d)
On this basis, the Jacobian
matrix of the transformation,
of components
J, = 8Xj lax;, is written as
(3) and its determinant
can be expressed
in the following
form
j(x) = det(J) = Ykm X, = - X:m Y, = j"+ x j,+ y j,, where m is the skew-symmetric m =m’+xm,
(4)
matrix
+ym,,
W
1 ino= ~(y,xi -x,yb);
1 m, = +c
xi - x,ct);
1 my = -$y,c’
-c yi>.
(5b-d)
and j”= YLrn’X, (and so on). Here and in the sequel, the symbol ()” is used to indicate the values assumed by the relevant quantities at the element centroid (e.g., j” =j(x =O)). On the parent element, the displacement model z&)=iV(x)U is defined by the usual bilinear shape functions and the compatible strain distribution Vr&)=B(x)U is derived from it through customary isoparametric procedures. The explicit expression of matrix B(x) reads
(6)
The element rigid body motions r are separated from generalized modes, by introducing the non-singular transformation [7]
u = [a, Al{;};
{;}=.,A,-tl=[qu,
where, by the very definition for A, is [8]
(or natural)
strains q, controlling
the straining
U&b)
of rigid body motions, it must be B(x),4 r = 0. It can be verified that the expression
83
A. Capsoni, L. Corradi I Comput. Methods Appl. Mech. Engrg. 141 (1997) 81-93
(8) In vector q, constant and linear (hourglass) modes are separated by writing qf ={qolq,}t; following expression for matrix A is assumed
A,
A = L%,A,l;
rn;;I; =[:jJ”)-‘.
=[:
A,
consistently, the
(9a-c)
From Eq. (6) and (9) one obtains
(100) where To is the trareformation matrix 8, T=
Jt2
[ J,,J,,
2J,, J*, J& 2J,,J22 JnJm J,,J,, + JnJx J;,
(11)
1 ’
computed at the element centroid. Eqs. (lOa,b) define the matrix b(x) expressing compatible strains as function of q through [ [l], Equation (23a)l. Matrix [A, A] can be inverted in closed form by introducing four vectors defined by the condition [a p, pY 71’ =: b -q y, cl-‘.
(12)
Such vectors exist if the isoparametric mapping is non degenerate [8] and, on their basis, matrices C, and C in Eq. (7b) can be computed. They read t
(13a,b) Eq. (13b) expresses the element compatibility matrix, The very fact that it is possible to define it, states that the 5X8 matrix C is of full rank. The models for static fields and for strains [ [l], Eqs. (26), (31)] are assumed to be as follows 0
3y-j,ljO s(x) = +j[~,lT’%,(x)l;
i,(x) =
0
3x-j,/jO
0
~(T”)-f&2(x)
1
9
( 14a,b)
0
1 ,
(154
(15b,c) Finally, plastic multipliers i are bilinearly interpolated over the Gauss points of a 2X2 grid. The row matrix L(x) in [ [l], Eq. (35)] then reads L(x) = IL, (x)r,,(x)L,(x)L,(x)l;
q$x) = $1 + 3x,x)(1 +
where xx, yR (5 1/-hi> are the Gauss point coordinates.
3YgY),
( l&b)
84
A. Capsoni, L. Corradi
I Comput. Methods Appl. Mech. Engrg.
141 (1997)
81-93
A comment on the above models is in order. The expression Eq. (9b) of A, is demanded by the condition BA, =I, while Eq. (SC) is, in principle, only one among several possible choices. A, is the product of two matrices, the first governing the hourglass modes in the undistorted (parent) element; to preserve invariance under the isoparametric transformation, this matrix is multiplied by J-‘, for simplicity evaluated at the element centroid [2]. In Eq. (lo), this role is played by matrix To which, for the same reason, also affects modes other than constants in Eqs. (14), (15). The stress model Eqs. (14) is essentially the one of Pian and Sumihara 133, in which the origin of linear modes was shifted so as to comply with the biorthogonality condition [ [l], Eq. (27)]. The strain model is selected so that shearing strains are constant in the undistorted element. It can be verified that the requirements [[ 11, Eq. (33)] are met and that condition [ [l], Eq. (34)], ruling out locking at the incompressibility limit, is also fulfilled. Moreover, stress and strain models do not include those derived by the displacement approximation and limitation principles do not apply.
3. Finite loading steps In computations rates have to be replaced by finite increments and the problem is now formulated in this context. The state of the system is known in configuration 2 and a transition 2 + 2 = x + AZ, in which loads undergo increments AF from their values P in 2, is considered. The compatibility and equilibrium conditions and the elastic part of the constitutive law are linear and hold independently of the step amplitude. The definitions x=t3VKlt3P, p=aV, /a~ of static internal variables are preserved, while the normality law and Prager’s consistency rule are replaced by the relations AP ={g}xAk AAsO;
ZAA=AA. &SO;
&AA = 0,
( 18a-c)
where a superposed cap denotes the values assumed at the end of the step. In Eqs. (17) normals are evaluated at the final position, consistently with a backward difference integration of the rate relationships. Eq. (18) establish the step-wise reversible, or holonomic, nature of the law: the two sign restrictions and the complementarity condition in fact imply AA=0 if & < 0, showing that the irreversible nature of plastic behavior is no longer accounted for within each step, even if unloading from a plastic state is possible at the beginning of the step itself. The variational statement generalizing [ [l], Eq. (21)] to finite steps was established in [ [ 11, Ref. 161. The relevant functional reads *=I
R ;(Ae [
- A.P)‘D(A& - AZ’) + Aa’V(Au) - A&(AE - Al’) + v,(P) - i%f’
- (i% - ?)‘fiAh
1
dx - W(Au),
+ V,(i)
(19)
where n(a,x) = (a@/&) is the normal to the plane strain yield function and 3 =n( &, i). Note that the terms involving AA are now integrated over the entire volume, since yield functions can be activated even if it is di < 0 when the step starts. For simplicity, linear hardening only is considered. In that case, the expressions for the two stored energy contributions read [ [l], Ref. 181 V,(P) = +P’H P;
V,(A)=ooA+;A2,
(2Oa,b)
where H=(E/(1-2~)*)~+3cM-‘, r and M being the matrices defined by [ [l], Eqs. (16b) and (lSc)]. In the above relations, a0 is the initial yield limit and c and a are non-negative quantities, possibly varying with x but staying unaltered during the elastic-plastic evolution. The discrete functional is obtained by introducing the trial functions. Their structures provide some simplifications with respect to the general case. In particular, the strain model Eqs. (15) makes submatrix W h in
A. Capsoni, L. Corradi I Comput. Methods Appl. Mech. Engrg. 141 (1997) 81-93
85
the expression [[ 11, IS+ (39)] of W* non-singular even in the perfectly plastic case. Therefore, this matrix can be condensed as for E* and the parameters governing enhanced fields directly eliminated (obviously, for perfect plasticity the expression [ [l], Eq. (6?b)] of W, is recovered). When account is taken of this circumstance, the following problem is arrived at min AU,
max
Aq. Am. AA AL!. AX
Ic,subject to AA 3 0,
(214
where, $ = ;(Aq
- Arr)‘E(Aq - A?r) + AQ’CAV - AQ’(Aq - Am)
+ +Wii
- f’An
+ R;AA + $&A
- (Q - f)‘fiAA
- AF’AV,
(21b)
where,
R, = R L’(x)cr,,(x) dx, I and the remaining conditions read
matrices
(22) and vectors are defined as, in Ref. [ 11. It is readily verified that the saddle point
C’AQ = AF;
Aq = CAV;
AX = WAn;
Arr = fiAA,
AQ = E(Aq - A@,
(23a-c) (24a,b)
@=A’(&-&=wAA-~~o;
AA>O;
@AA
= 0,
(24c-e)
where E = R, + WA is the generalized yield limit in configuration 2. It is worth mentioning that the model Bq. (16) for A implies that the plastic portion of the constitutive law for the element simply enforces the material law at the Gauss points [ [ 11, Ref. 61. However, as Bq. (24~) indicates, vector !P depends on the parameters governing the static fields and the elastic range for each element is the intersection of four domains in the generalized stress space, bounded by a non-regular yield surface; hence, the problem Eqs. (24) must be handled at the element level. By performing integrations numerically, account taken of Eq. (16), one obtains
w,@kg1
(25)
The Gauss weights on a 2 X 2 grid being unitary, it is W,=j(xK); for non-degenerate quantities, which ensures that the current elastic domain ?PsO is convex.
elements,
these are positive
4. Solution strategy Eqs. (23a-c), governing AU = K-‘(AI;
the elastic part of the problem,
+ C’EArr);
AQ = E(CAV - Am)
are easily solved to give (26ab)
where K=C’EC is the elastic stiffness matrix for the discrete scheme. The non-linear problem Eqs. (24), enforcing compliance with the plasticity condition, is now considered. An elastic predictor/plastic corrector strategy, similar to that proposed in Ref. [51, is employed for its solution. However, to account for perfect plasticity, when W is singular and cannot be inverted, the role of driving variables is assigned to plastic strains rather than to static internal variables. The procedure is initiated by an elastic trial, defining the stress increment predictor as
A. Capsoni, L. Corradi
86
AQ’E C K-‘AF;
Q”=&
I Comput. Methods Appl. Mech. Engrg.
141 (1997) 81-93
+AQ”.
(27)
On this basis, the current values No of the outward normal matrix and w”= [N”]‘(Qo-8)-R of the generalized yield functions are computed. If it is Y” ~0 componentwise, the process is elastic and no corrections are required. In general, however, a subvector of 9’ will be positive, namely !Pu”,[No,]‘(Q+ AQ” -s)
- 8, > 0.
(28)
It is tentatively assumed that only, and all, the components of this vector are activated in the process (VW=0, AA,>O). Then, Eqs. (24) reduce to the problem rP = Am - N,AA, = 0,
(29a)
r,=E-‘(AQ-AQO)+N,AAa=O,
(29b)
rw=NL(Q+AQ-_%-WAn)-Ww,AA,-Z?,=O,
(29~)
where N, =(alu,laQ) is a function of the difference Q-X. Eqs. (29a), (29b), (29~) are solved by subsequent linearization. The symbol SO’+’ is used for increments in each iteration, to be added to previously computed values A()’ (i=O refers to the elastic trial). Note that [ [I], JZq. (15)] define CDas the difference between the positively homogeneous, of degree one, function F=F(a -x) and the yield limit /? and that EZq. (25) establishes that the components of vector W have the same structure. This implies 8[Nb@ + AQ -8- WAtr)] = Nh(8Q - W&r) and WN,A4)
= NJ4
AAk] (SQ -8X) = N,S/1, + M,(8Q -W&r),
+ [z, (2)’
(30)
where the matrix, (31) is symmetric and positive semidefinite as long as AA,aO. Each iteration is brought to the solution of the problem t:’
=a&+’
r;’
= E-’
,.;I
=
-Nb,s/i;’ SQ’”
+
N’,
_M’$jQi+‘_W&‘+‘)+r;=O, SAz’
+
M’,(SQ’+’
_
W
,&+‘)
(324 +
r;
=
0,
-WSmTi+‘)-w,6A~‘+r~=0.
[N’,]‘(aQ’+’
(32b) (32~)
Let Eqs. (32a) be written as [I + M’,W] Sa’+’ = N’, SAZ’ + ML SQ”’
- rb.
(33)
The following properties can be established (proofs are given in Appendix A): (a) matrix [Z+ M, W] is non singular and one can write z=[Z+M,W]-‘,
(34)
(b) the following matrices are symmetric and at least positive semidefinite .n, = wg
a, = EMm;
(35a,b)
(36a-c)
(cl P=z-M,w~=z-=M~w=~-M,=:"'w
Property (a) states that Eq. (33) can be solved for a&+‘; if the result is introduced into Eqs. (32b), (32c), and account taken of Eqs. (36), one obtains [iywi+l
+
;pAb”
+
7
‘e =
o,
(374
A. Capsoni, L. Corradi [$J’se”‘+
;i$Ab+l
I Comput. Methods Appl. Mech. Engrg.
141 (1997)
+;‘,=o,
87
81-93
Wb)
where, =E-’
[ii]-’
_k @. ,,
+E’N;;
;b, = W, + [N’,]‘@N’,,
;;&+&J+;;
(38a,b)
,i,=,,+N;+;.
(39a-c)
The two matrices defned by Eqs. (38) are symmetric &nd positive definite and semidefinite respectively, since and w are so and because of property (b) above. If N’, is of full column rank, Eqs. (37a), (37b) yield
E
(404
w’”
=
6~21
Ei$i$‘,
+
FL>.
The current estimate of the step solution Am’+’ =A&+&+‘;
(Job) is now obtained
by writing
A,$+A&+a/l;‘;
AQ”
= AQi + @+I.
(41a-c)
The matrices in Eqs. (32a), (32b), (32~) are updated on this basis and the process is repeated until residuals become sufficiently small. Two points in the above procedure require a comment. In some instances (when w =0 and the columns of matrix &b, turn out to be linearly dependent) the matrix to be inverted in Eqs. (40a) is singular even prior to collapse. This is a consequence of the non-regular nature of the element generalized yield function and simply means that the solutmn of the linearized problem is no longer unique in terms of S/i,, but it is still so for 6m and 8Q. To account for this occurrence (which is rare in plane strain situations due to the stabilizing effect of fictitious hardening) _.a pivotal strategy for matrix inversion is employed, which allows for the discarding of the excess columns of N. The procedure, borrowed from Mathematical Programming algorithms, was successfully tested for piecewise linear plasticity [ [l], Ref. 221 and has proved effective in this context as well. More relevant is that Eqs. (24c-e) were handled by assuming that the yield functions active in the step are those indicated by the elastic trial. It may happen that Eqs. (40a) predicts a negative value for some components of AA,, or that some components pk, not originally in !&, become positive. When this occurs, subvectors and submatrices are redefined. For such a procedure convergence cannot be proved, but inconveniences were never experienced. At this point the constitutive law, Eqs. (24), are complied with, but equilibrium is to be restored. Conceptually, this can be done by introducing the computed value of i% in Eqs. (26): an iterative scheme is thus produced, enforcing the elastic-plastic constitutive law inside equilibrium restoring loops. This would amount to a modified Newton-Raphson method, which exhibits slow convergence, and more efficient strategies are available. The one employed is essentially that proposed in Ref. [5]; the n-th equilibrium iteration solves the linear problem J&AU,
= R, _ , ,
(42)
where,
Kep= C'i,,C: E,,=E-E&w, is the tangent
(elastic-plastic)
R=&C'@
algorithmic
+AQ)=AF-C’AQ,
++g+jir$, stiffness
(430)
matrix, and (43c)
is the unbalanced nodal force vector. In Eqs. (43), all solution dependent quantities have final values at the end of the interior iterations that enforce the plasticity conditions. The procedure is standard and does not exhibit any peculiarities worth discussing.
88
A. Capsoni, L. Corradi
I Comput. Methods Appl. Mech. Engrg.
141 (1997)
81-93
5. Numerical examples The performances of the proposed model are demonstrated by solving some problems which are known as pathological. The first is the perfectly plastic double-edge-notched tensile specimen, for which a displacement approach is unable to predict the limit load [9]. Dimensions are indicated in Fig. 1 and the following values are assumed for the material constants: Young’s modulus E =70,Poisson ratio v=O.3 and uniaxial yield limit a, = 0.243. The accepted value of the limit load is FL = b(2 + a-)~,, / \15= 1.44. One quarter of the specimen is subdivided into 15 X5 square elements. 35 displacement controlled steps are performed in the interval 0~(EA)l(a,W)~20, where A is the relative displacement between the specimen ends. The computed load-displacement curve, illustrated in Fig. 2, approaches asymptotically the value FL. A second example is a ring of mean radius R = 10.5 and thickness h = 1.O, subject to equal and opposite forces (Fig. 3). This is a classical elastic benchmark case [lo], also studied as a plastic plane stress problem [6], and it is now analyzed in plane strain conditions. It can be considered a severe plastic test, since the structure exhibits an essentially flexural behavior and plastic strains are confined to limited regions about the two symmetry axes. By neglecting shear effects, the von Mises yield function is brought to the following form [ 1 l]
where N is the axial force and M the bending moment. The ring being statically determinate in terms of N, a mechanism involving plastic hinges in the four symmetry sections predicts the kinematic load
FL=~N'R(I-dl-
8M;l(N,R)*),
0
which can be regarded as a good estimate of the collapse value. For an uniaxial yield limit a,=40. FL=4.41 The remaining material constants are assumed as E = lo6 and v=O.25. L=30
Fig.
1. Double-edge-notched
specimen.
iii;,:;
0
Fig. 2. Dimensionless
4
load-displacement
a
12
curve for the perfectly
16
plastic double-edge-notched
specimen.
one gets
A. Capsoni, L. Corradi I Comput. Methods Appl. Mech. Engrg. 141 (1997) 81-93
Fig. 3. Plane strain circular
ring problem
89
and mesh adopted.
One quarter of the ring is subdivided into 16 X 2 equal elements (Fig. 3) and the displacement w at the load point is increased by performing 50 steps of amplitude Aw =0.1(~0R21Eh). As Fig. 4 shows, the load approaches the value: FL = 4.54, exceeding the comparison value by 3% owing to the distance of the Gauss points from the plastic hinge locations. The figure also depicts the load-displacement curves produced by the displacement approach, still unable to predict collapse, and by the so-called B-bar method [12], which overestimates the limit load by more than 20%. The third example, an axisymmetric plane strain cylinder of internal and external radii r, = 3.0 and re = 9.0, subjected to internal pressure p, was also suggested as a test for element evaluation in the incompressible range [13]. Here, the material is assumed as perfectly plastic, with E = 1.0, v=O.3, a, =0.2, and its response up to collapse is computed. In a sense, the situation is opposite to the previous one, in that plasticity progressively spreads and eventually extends over the entire solid. The analytical limit pressure is [14] pL = $gJn(,-,
/ri) = 1.269v0 = 0.254,
The mesh proposed in Ref. [13] was assumed, with five unequal elements covering a sector of 10 degrees, and results show that the computed limit load exceeds by only 0.3% the analytical value (Fig. 5). The displacement approach still fails: a pressure p = 1.07 pL is predicted for the internal displacement u = 20 ri(cro/E) and the P-U curve keeps growing afterwards, with a constant even if small slope. The problem known as Cook’s problem, a tapered solid clamped at one end and loaded at the other by a
F 12.00 ,
000 000
1.00
200
Fig. 4. Load-displacement
300
4M)
curve for the circular
500
ring.
wEh sR
A. Capsoni, L. Corradi I Comput. Methods
90
Fig. 5. Axisymmetric
Appl. Mech. Engrg. 141 (1997) 81-93
plane strain cylinder-Mesh
Fig. 6. Cook’s problem4X4
and load-displacement
and 8X8
curve.
meshes.
parabolically distributed vertical force, is considered as a final example (Fig. 6, where the geometry is defined). Linear isotropic and kinematic hardening are now included by assuming, as in Refs. [2,6], E = 70.,
u = l/3,
u0 = 0.243,
c = 0.135,
a = 0.015,
c and a being the constants appearing in Eqs. (20). 18 loading steps of amplitude AF =O. 1 are applied. Three meshes are examined, with 4, 8 and 16 elements on each side, and results are depicted in Fig. 7. The intermediate mesh appears adequate and very little is gained by its refinement. The displacement for F = 1.8 is very nearly the same as predicted in Ref. [2]. Comparison of the load-displacement curve obtained with those produced, with the same mesh, by the displacement and B-bar methods demonstrates significant improvements (Fig. 8).
6. Concluding
remarks
A mixed finite element model for plane strain elastic-plastic analysis was proposed and its effectiveness was demonstrated by solving a number of meaningful examples. One of the main goals of the formulation was to
A. Capsoni, L. Corradi
I Comput. Methods Appl. Mech. Engrg.
141 (1997)
91
81-93
F 2.00
t 0400
2 112s
1
0.00 -f 000
Fig. 7. Load-displ,xement
I
I
I
I
,
0.50
100
1.50
2.00
2.50
curve for Cook’s problem-Effects
ow_r 000
Fig. 8. Load-displacement
I
I
I
I
loo
150
200
1
050
of mesh refinement
curves for Cook’s problem-Comparison
with alternative
u
(displacements
formulations
,
for F= 1.8 are indicated)
u
2.50
(displacements
for F=
I .8 are indicated).
avoid locking when plasticity spreads and strains become nearly deviatoric, while retaining the same stress redistribution capability of a displacement model. These features were assessed in Ref. [l] and numerical tests indicate that they are actually incorporated. In fact, for all the perfectly plastic cases the collapse load was predicted with excellent accuracy, but slightly overestimated, indicating that the mixed nature of the formulation does not entail an excess of flexibility. To avoid locking, the well known device of supplementing the strain approximation with enhanced fields was used, but the proposed formulation provides some simplifications. The procedure followed is conceptually analogous to that in Ref. [6] where, however, the current elastic-plastic stiffness matrix is condensed at each stage of the solution process. Here, the operation is performed on the elastic stiffness and hardening matrices separately; for perfect plasticity (or linear hardening) they are both constant and can be condensed once and for all, with significant saving. The stress redistribution capability was preserved by expressing the stress and strain approximations in terms of the generalized variables dictated by the displacement model. With respect to existing proposals, they introduce only mino:r modifications. However, a formulation in terms of generalized variables permits the use of numerical devices borrowed from Mathematical Programming concepts. Due to space limitations, this aspect was only partially discussed; the way of dealing with the non-regular element yield surface, however, is an example of the effectiveness of these tools. To this purpose, it is mentioned that, even if displacement control is preferable for perfect plasticity (and was systematically employed in the examples), in principle load control could also be used. In this case, the collapse load is detected by a. test on matrix A,, defined by [ [ 11, Eq. (67)], again exploiting Mathematical Programming theory. Obviously, the load step amplitude must become very small when approaching the collapse limit.
92
A. Capsoni, L. Corradi
I Comput. Methods Appl. Mech. Engrg.
141 (1997)
81-93
Examples also indicate that the discrete scheme performs well with coarse meshes and is rather insensitive to mesh distortion, an expected occurrence, since the trial functions used are essentially existing ones exhibiting this feature. Finally, it is underlined that the formulation rests on a constitutive law which was fully reduced to the plane by exploiting the deviatoric nature of plastic flow associated to the von Mises condition. This makes the finite element model specifically oriented to a particular, even if relevant, class of problems, but the significant simplifications which become possible more than compensate for the limited range of applicability.
Acknowledgments This research was supported by the Italian Ministry (MURST).
Appendix
of University
and Scientific
and Technological
Research
A
Proof of properties (a), (b) and (c) in Section 4. The solution strategy crucially relies on the fact that the matrix o=z+fW,w,
(A.1)
is non-singular and this assertion is now proved. It is recalled that both M, and W are symmetric and at least positive semidefinite block-diagonal matrices, each 5 X5 block referring to a finite element. Therefore, only the element level needs to be considered. The property is easily proved when material kinematic hardening is present and, hence, W is strictly positive definite. By writing @= [W-’ +M,]W it is immediately recognized that 0, as the product of two symmetric and positive definite matrices, has a positive determinant. For perfectly plastic or purely isotropically hardening materials, matrix W is given by [ [l], Eq. (67b)]; namely 11000 11000
EL! >o.
0 0
0 0
If 0 were singular, condition becomes
0 0
0 0
(A.2)
(1 - 2V)2
0 0
a vector x #O would exist such that Ox = 0. When account is taken of Eqs. (A. l), (A.2), the
=O.
For this equation solution 1 + a(M,, a&f,,
(A.3)
to be satisfied by a vector x #O, the following
homogeneous
system must have a non trivial
+M,,) (A.4)
+ M,,)
where the components
of M,
are denoted
as M,j.The matrix
being
symmetric
and positive
semidefinite,
M,,=M,, and (ASa-c)
M ,,20, M,, a09 M,,M,,-M;, 20. The condition
that the determinant
of the coefficient
matrix in Eqs. (A.4) be equal to zero yields
93
A. Capsoni, L. Corradi I Comput. Methods Appl. Mech. Engrg. 14I (1997) 81-93
M,, +q2
(A.6)
+2h4,, = +o,
Eqs. (A.Sa,b) and (A.6) imply 2M,,
+M*,+2M,,)(M,,
taking account
I7=(M,,
and, hence, M,,+M,,
- 2M,, ~0, thus establishing
the inequality (A.7a)
+&-2M,,)
of Eq. (A.~c), one can also write
+ M,,)* -4M;,
3
(M,,+ M,2)2- 4M,, M,,=(M,,- M,,)'>O.
(A.7b)
The contradiction between Eqs. (A.7a), (A.7b) proves that the equation OX = 0 does not have solutions other than x =0 and, hence, that matrix 0 is non-singular. Property (a) having now been assessed, one can write 0 =Z +M,W= E ’ and the condition ss-’ = u-l L c=z E ==Z-EMu\c Taking
S=Z-MM,Ws,
account of the symmetry
M,S’
of both Ma and W, from Eqs. (A.8) one obtains E’W = W - &‘WM,W,
= M, - M,WM,=‘;
hence M, = [Ma + M, WM,]$, 0, =Mu = =[M,
64.8)
W= P’[W f WM, W] and
- Ma WM,s:“;
d2, = WEs’= E”[W-
WM,W]E
64.9)
The above relations show that &?, and az are symmetric, positive semidefinite matrices, since W and M, are individually so, and assesses property (b). This result also implies sMo =M,Et, WE= s:“‘W and from Eqs. (A.8) one obtains =?=I-M,z;“‘W, Together
with Eqs. l(A.8) the above relation
proves property
(c).
References [I] A. Capsoni and L. Corradi, A mixed finite element model for plane strain elastic-plastic analysis: Part I-Formulation and assessment of the overall behavior, Comp. Methods Appl. Mech. Engrg. 141 (1997) 67. [2] J.C. Simo and MS. Rifai, A class of assumed strain methods and the method of incompatible modes, Int. J. Numer. Methods Engrg. 29 (1990) 1595. [3] T.H.H. Pian and H. Sumihara, Rational approach for assumed stress finite elements, Int. .J. Numer. Methods Engrg. 20 (1984) 1682. [4] J.B. Martin, Plasticity: Theory and Applications (M.I.T. Press, Cambridge, MA, 1975). [5] J.C. Simo, J.G. Kennedy and S. Govindjee, Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms, Int. J. Numer. Methods Engrg. 26 (1988) 2161. [6] S.L. Weissman and M. Jamjian, Two-dimensional elastoplasticity-Approximation by mixed finite elements, Int. I. Numer. Methods Engrg. 36 (1993) 3703. [7] J.H. Argyris, Contitma and Discontinua, in: Proc. 1st Conf. Matrix Methods in Structural Mechanics, AFFDL, TR 66-88, Dayton, OH (1966) 11-92. [8] T. Belytschko, J.S. Ong, W.K. Liu and J.M. Kennedy, Hourglass control in linear and nonlinear problems, Comp. Methods Appl. Mech. Engrg. 43 (1984) 251. ]9] J.C. Nagtegaal, D.M. Parks and J.R. Rice, On numerically accurate finite element solutions in the fully plastic range, Comp. Methods Appl. Mech. Engrg. 4 (1974) 153. [IO] T. Belytschko and E. Bachrach, Efficient implementation of quadrilaterals with high coarse-mesh accuracy, Comp. Methods Appl. Mech. Engrg. 54 (1986) 279. [l I] W.F. Chen and D.J Han, Plasticity for structural engineers (Springer Verlag, 1988). [12] T.J.R. Hughes, Generalization of selective integration procedures to anisotropic and nonlinear media, Int. J. Numer. Methods Engrg. 15 (1980) 1413. [ 131 R.H. MacNeal and !R.L. Harder, A proposed standard set of problems to test finite element accuracy, Finite Element Analysis Design 1 (1985, 3. [14] W. Prager and Ph.G.Jr. Hodge, Theory of perfectly plastic solids (Wiley, New York, 195 I).
ElSEVlER
Comput.
Methods
Appl. Mech. Engrg.
141 (1997) 81-93
A mixed finite element model for plane strain elastic-plastic analysis Part II. Application to the 4-node bilinear element A. Capsoni, L. Corradi” Department
oj Structural
Engineering,
Politecnico
Received
di Milano,
3 1 March
Piazza Leonardo
1995; revised
da Vinci 32, Milano,
I-20133,
Italy
18 April 1996
Abstract In this paper, the results of part I of this study are applied to the 4-node bilinear element. The trial functions for different fields consistently with the a-priori conditions established, the rate problem is generalized to cover finite loading steps and the solution presented. Computations referred to some discriminating examples from the literature demonstrate that the proposed model capable of predicting the perfectly plastic collapse load for pathological plane strain problems, exhibits excellent performances meshes and is rather insensitive to mesh distortion.
are defined strategy is is actually with coarse
1. Introduction The first part of this study [l] was devoted to a mixed finite element formulation of the plane strain problem for elastically isotrolpic solids with a plastic behavior governed by von Mises condition. It was shown that if the trial functions for different fields comply with some conditions, the resulting discrete scheme is able to prevent locking when plasticity spreads, while preserving its stress redistribution capability. The main purpose of this contribution is to demonstrate, with reference to the 4-node bilinear element, how such a model can be actually established and used for computations. General results are specialized to the element considered .and suitable approximations for all the variables involved are presented; the rate problem is then generalized to cover finite loading steps, as required by computations; the solution strategy is next presented and, finally, the efficiency of the formulation is .demonstrated by solving some problems which are known as pathological. None of the above steps contains contributions that are actually new. The models for different fields are essentially ones that have been already proposed [2,3], with minor modifications introduced to meet the a-priori requirements. The transition from rates to finite increments rests on a standard backward difference integration of the rate law, and ,aclassical step-wise holonomic problem [4] is obtained as a result. The solution algorithm is a strain-driven return mapping one, adapting existing procedures to the problem at hand [5,6]. The model proves capable of correctly predicting perfectly plastic collapse, is highly accurate in coarse meshes and rather insensitive to shape distortion, but its performances are shared by other formulations. Nevertheless, it is felt that the proposal deserves some attention for at least two reasons. First, the problem is formulated in terms of generalized variables in Prager’s sense, which permits the use of
*Corresponding
author.
00457825197/$17.00 0 1997 Elsevier Science S.A. All rights reserved PII SOO45-7825(96)01099-7
82
A. Capsoni, L. Corradi 1 Comput. Methods Appl. Mech. Engrg. 141 (1997) 81-93
results available in this context. In particular, some devices resting on Mathematical Programming theory [ [I], Ref. 221 can be exploited at different stages of the solution process to overcome numerical difficulties. Secondly, the constitutive law is fully reduced to the plane; except that, for the expressions of matrices governing the material behavior, the elastoplastic plane strain problem is driven to a form identical to a plane stress one and, as in the elastic case, both problems can be dealt with in the same manner, with significant conceptual and computational simplifications.
2. Definition of field approximations The isoparametric 4-node element is considered. Actual and intrinsic coordinates and x = {x y}‘, respectively. Nodal coordinates are collected in the vectors X, = 1X, X, X, X4)‘; and the following
are denoted by X = {X Y}’
(1)
Y, = {Y, Y, Y3 YqK
definitions
are introduced
(2a-d)
On this basis, the Jacobian
matrix of the transformation,
of components
J, = 8Xj lax;, is written as
(3) and its determinant
can be expressed
in the following
form
j(x) = det(J) = Ykm X, = - X:m Y, = j"+ x j,+ y j,, where m is the skew-symmetric m =m’+xm,
(4)
matrix
+ym,,
W
1 ino= ~(y,xi -x,yb);
1 m, = +c
xi - x,ct);
1 my = -$y,c’
-c yi>.
(5b-d)
and j”= YLrn’X, (and so on). Here and in the sequel, the symbol ()” is used to indicate the values assumed by the relevant quantities at the element centroid (e.g., j” =j(x =O)). On the parent element, the displacement model z&)=iV(x)U is defined by the usual bilinear shape functions and the compatible strain distribution Vr&)=B(x)U is derived from it through customary isoparametric procedures. The explicit expression of matrix B(x) reads
(6)
The element rigid body motions r are separated from generalized modes, by introducing the non-singular transformation [7]
u = [a, Al{;};
{;}=.,A,-tl=[qu,
where, by the very definition for A, is [8]
(or natural)
strains q, controlling
the straining
U&b)
of rigid body motions, it must be B(x),4 r = 0. It can be verified that the expression
83
A. Capsoni, L. Corradi I Comput. Methods Appl. Mech. Engrg. 141 (1997) 81-93
(8) In vector q, constant and linear (hourglass) modes are separated by writing qf ={qolq,}t; following expression for matrix A is assumed
A,
A = L%,A,l;
rn;;I; =[:jJ”)-‘.
=[:
A,
consistently, the
(9a-c)
From Eq. (6) and (9) one obtains
(100) where To is the trareformation matrix 8, T=
Jt2
[ J,,J,,
2J,, J*, J& 2J,,J22 JnJm J,,J,, + JnJx J;,
(11)
1 ’
computed at the element centroid. Eqs. (lOa,b) define the matrix b(x) expressing compatible strains as function of q through [ [l], Equation (23a)l. Matrix [A, A] can be inverted in closed form by introducing four vectors defined by the condition [a p, pY 71’ =: b -q y, cl-‘.
(12)
Such vectors exist if the isoparametric mapping is non degenerate [8] and, on their basis, matrices C, and C in Eq. (7b) can be computed. They read t
(13a,b) Eq. (13b) expresses the element compatibility matrix, The very fact that it is possible to define it, states that the 5X8 matrix C is of full rank. The models for static fields and for strains [ [l], Eqs. (26), (31)] are assumed to be as follows 0
3y-j,ljO s(x) = +j[~,lT’%,(x)l;
i,(x) =
0
3x-j,/jO
0
~(T”)-f&2(x)
1
9
( 14a,b)
0
1 ,
(154
(15b,c) Finally, plastic multipliers i are bilinearly interpolated over the Gauss points of a 2X2 grid. The row matrix L(x) in [ [l], Eq. (35)] then reads L(x) = IL, (x)r,,(x)L,(x)L,(x)l;
q$x) = $1 + 3x,x)(1 +
where xx, yR (5 1/-hi> are the Gauss point coordinates.
3YgY),
( l&b)
84
A. Capsoni, L. Corradi
I Comput. Methods Appl. Mech. Engrg.
141 (1997)
81-93
A comment on the above models is in order. The expression Eq. (9b) of A, is demanded by the condition BA, =I, while Eq. (SC) is, in principle, only one among several possible choices. A, is the product of two matrices, the first governing the hourglass modes in the undistorted (parent) element; to preserve invariance under the isoparametric transformation, this matrix is multiplied by J-‘, for simplicity evaluated at the element centroid [2]. In Eq. (lo), this role is played by matrix To which, for the same reason, also affects modes other than constants in Eqs. (14), (15). The stress model Eqs. (14) is essentially the one of Pian and Sumihara 133, in which the origin of linear modes was shifted so as to comply with the biorthogonality condition [ [l], Eq. (27)]. The strain model is selected so that shearing strains are constant in the undistorted element. It can be verified that the requirements [[ 11, Eq. (33)] are met and that condition [ [l], Eq. (34)], ruling out locking at the incompressibility limit, is also fulfilled. Moreover, stress and strain models do not include those derived by the displacement approximation and limitation principles do not apply.
3. Finite loading steps In computations rates have to be replaced by finite increments and the problem is now formulated in this context. The state of the system is known in configuration 2 and a transition 2 + 2 = x + AZ, in which loads undergo increments AF from their values P in 2, is considered. The compatibility and equilibrium conditions and the elastic part of the constitutive law are linear and hold independently of the step amplitude. The definitions x=t3VKlt3P, p=aV, /a~ of static internal variables are preserved, while the normality law and Prager’s consistency rule are replaced by the relations AP ={g}xAk AAsO;
ZAA=AA. &SO;
&AA = 0,
( 18a-c)
where a superposed cap denotes the values assumed at the end of the step. In Eqs. (17) normals are evaluated at the final position, consistently with a backward difference integration of the rate relationships. Eq. (18) establish the step-wise reversible, or holonomic, nature of the law: the two sign restrictions and the complementarity condition in fact imply AA=0 if & < 0, showing that the irreversible nature of plastic behavior is no longer accounted for within each step, even if unloading from a plastic state is possible at the beginning of the step itself. The variational statement generalizing [ [l], Eq. (21)] to finite steps was established in [ [ 11, Ref. 161. The relevant functional reads *=I
R ;(Ae [
- A.P)‘D(A& - AZ’) + Aa’V(Au) - A&(AE - Al’) + v,(P) - i%f’
- (i% - ?)‘fiAh
1
dx - W(Au),
+ V,(i)
(19)
where n(a,x) = (a@/&) is the normal to the plane strain yield function and 3 =n( &, i). Note that the terms involving AA are now integrated over the entire volume, since yield functions can be activated even if it is di < 0 when the step starts. For simplicity, linear hardening only is considered. In that case, the expressions for the two stored energy contributions read [ [l], Ref. 181 V,(P) = +P’H P;
V,(A)=ooA+;A2,
(2Oa,b)
where H=(E/(1-2~)*)~+3cM-‘, r and M being the matrices defined by [ [l], Eqs. (16b) and (lSc)]. In the above relations, a0 is the initial yield limit and c and a are non-negative quantities, possibly varying with x but staying unaltered during the elastic-plastic evolution. The discrete functional is obtained by introducing the trial functions. Their structures provide some simplifications with respect to the general case. In particular, the strain model Eqs. (15) makes submatrix W h in
A. Capsoni, L. Corradi I Comput. Methods Appl. Mech. Engrg. 141 (1997) 81-93
85
the expression [[ 11, IS+ (39)] of W* non-singular even in the perfectly plastic case. Therefore, this matrix can be condensed as for E* and the parameters governing enhanced fields directly eliminated (obviously, for perfect plasticity the expression [ [l], Eq. (6?b)] of W, is recovered). When account is taken of this circumstance, the following problem is arrived at min AU,
max
Aq. Am. AA AL!. AX
Ic,subject to AA 3 0,
(214
where, $ = ;(Aq
- Arr)‘E(Aq - A?r) + AQ’CAV - AQ’(Aq - Am)
+ +Wii
- f’An
+ R;AA + $&A
- (Q - f)‘fiAA
- AF’AV,
(21b)
where,
R, = R L’(x)cr,,(x) dx, I and the remaining conditions read
matrices
(22) and vectors are defined as, in Ref. [ 11. It is readily verified that the saddle point
C’AQ = AF;
Aq = CAV;
AX = WAn;
Arr = fiAA,
AQ = E(Aq - A@,
(23a-c) (24a,b)
@=A’(&-&=wAA-~~o;
AA>O;
@AA
= 0,
(24c-e)
where E = R, + WA is the generalized yield limit in configuration 2. It is worth mentioning that the model Bq. (16) for A implies that the plastic portion of the constitutive law for the element simply enforces the material law at the Gauss points [ [ 11, Ref. 61. However, as Bq. (24~) indicates, vector !P depends on the parameters governing the static fields and the elastic range for each element is the intersection of four domains in the generalized stress space, bounded by a non-regular yield surface; hence, the problem Eqs. (24) must be handled at the element level. By performing integrations numerically, account taken of Eq. (16), one obtains
w,@kg1
(25)
The Gauss weights on a 2 X 2 grid being unitary, it is W,=j(xK); for non-degenerate quantities, which ensures that the current elastic domain ?PsO is convex.
elements,
these are positive
4. Solution strategy Eqs. (23a-c), governing AU = K-‘(AI;
the elastic part of the problem,
+ C’EArr);
AQ = E(CAV - Am)
are easily solved to give (26ab)
where K=C’EC is the elastic stiffness matrix for the discrete scheme. The non-linear problem Eqs. (24), enforcing compliance with the plasticity condition, is now considered. An elastic predictor/plastic corrector strategy, similar to that proposed in Ref. [51, is employed for its solution. However, to account for perfect plasticity, when W is singular and cannot be inverted, the role of driving variables is assigned to plastic strains rather than to static internal variables. The procedure is initiated by an elastic trial, defining the stress increment predictor as
A. Capsoni, L. Corradi
86
AQ’E C K-‘AF;
Q”=&
I Comput. Methods Appl. Mech. Engrg.
141 (1997) 81-93
+AQ”.
(27)
On this basis, the current values No of the outward normal matrix and w”= [N”]‘(Qo-8)-R of the generalized yield functions are computed. If it is Y” ~0 componentwise, the process is elastic and no corrections are required. In general, however, a subvector of 9’ will be positive, namely !Pu”,[No,]‘(Q+ AQ” -s)
- 8, > 0.
(28)
It is tentatively assumed that only, and all, the components of this vector are activated in the process (VW=0, AA,>O). Then, Eqs. (24) reduce to the problem rP = Am - N,AA, = 0,
(29a)
r,=E-‘(AQ-AQO)+N,AAa=O,
(29b)
rw=NL(Q+AQ-_%-WAn)-Ww,AA,-Z?,=O,
(29~)
where N, =(alu,laQ) is a function of the difference Q-X. Eqs. (29a), (29b), (29~) are solved by subsequent linearization. The symbol SO’+’ is used for increments in each iteration, to be added to previously computed values A()’ (i=O refers to the elastic trial). Note that [ [I], JZq. (15)] define CDas the difference between the positively homogeneous, of degree one, function F=F(a -x) and the yield limit /? and that EZq. (25) establishes that the components of vector W have the same structure. This implies 8[Nb@ + AQ -8- WAtr)] = Nh(8Q - W&r) and WN,A4)
= NJ4
AAk] (SQ -8X) = N,S/1, + M,(8Q -W&r),
+ [z, (2)’
(30)
where the matrix, (31) is symmetric and positive semidefinite as long as AA,aO. Each iteration is brought to the solution of the problem t:’
=a&+’
r;’
= E-’
,.;I
=
-Nb,s/i;’ SQ’”
+
N’,
_M’$jQi+‘_W&‘+‘)+r;=O, SAz’
+
M’,(SQ’+’
_
W
,&+‘)
(324 +
r;
=
0,
-WSmTi+‘)-w,6A~‘+r~=0.
[N’,]‘(aQ’+’
(32b) (32~)
Let Eqs. (32a) be written as [I + M’,W] Sa’+’ = N’, SAZ’ + ML SQ”’
- rb.
(33)
The following properties can be established (proofs are given in Appendix A): (a) matrix [Z+ M, W] is non singular and one can write z=[Z+M,W]-‘,
(34)
(b) the following matrices are symmetric and at least positive semidefinite .n, = wg
a, = EMm;
(35a,b)
(36a-c)
(cl P=z-M,w~=z-=M~w=~-M,=:"'w
Property (a) states that Eq. (33) can be solved for a&+‘; if the result is introduced into Eqs. (32b), (32c), and account taken of Eqs. (36), one obtains [iywi+l
+
;pAb”
+
7
‘e =
o,
(374
A. Capsoni, L. Corradi [$J’se”‘+
;i$Ab+l
I Comput. Methods Appl. Mech. Engrg.
141 (1997)
+;‘,=o,
87
81-93
Wb)
where, =E-’
[ii]-’
_k @. ,,
+E’N;;
;b, = W, + [N’,]‘@N’,,
;;&+&J+;;
(38a,b)
,i,=,,+N;+;.
(39a-c)
The two matrices defned by Eqs. (38) are symmetric &nd positive definite and semidefinite respectively, since and w are so and because of property (b) above. If N’, is of full column rank, Eqs. (37a), (37b) yield
E
(404
w’”
=
6~21
Ei$i$‘,
+
FL>.
The current estimate of the step solution Am’+’ =A&+&+‘;
(Job) is now obtained
by writing
A,$+A&+a/l;‘;
AQ”
= AQi + @+I.
(41a-c)
The matrices in Eqs. (32a), (32b), (32~) are updated on this basis and the process is repeated until residuals become sufficiently small. Two points in the above procedure require a comment. In some instances (when w =0 and the columns of matrix &b, turn out to be linearly dependent) the matrix to be inverted in Eqs. (40a) is singular even prior to collapse. This is a consequence of the non-regular nature of the element generalized yield function and simply means that the solutmn of the linearized problem is no longer unique in terms of S/i,, but it is still so for 6m and 8Q. To account for this occurrence (which is rare in plane strain situations due to the stabilizing effect of fictitious hardening) _.a pivotal strategy for matrix inversion is employed, which allows for the discarding of the excess columns of N. The procedure, borrowed from Mathematical Programming algorithms, was successfully tested for piecewise linear plasticity [ [l], Ref. 221 and has proved effective in this context as well. More relevant is that Eqs. (24c-e) were handled by assuming that the yield functions active in the step are those indicated by the elastic trial. It may happen that Eqs. (40a) predicts a negative value for some components of AA,, or that some components pk, not originally in !&, become positive. When this occurs, subvectors and submatrices are redefined. For such a procedure convergence cannot be proved, but inconveniences were never experienced. At this point the constitutive law, Eqs. (24), are complied with, but equilibrium is to be restored. Conceptually, this can be done by introducing the computed value of i% in Eqs. (26): an iterative scheme is thus produced, enforcing the elastic-plastic constitutive law inside equilibrium restoring loops. This would amount to a modified Newton-Raphson method, which exhibits slow convergence, and more efficient strategies are available. The one employed is essentially that proposed in Ref. [5]; the n-th equilibrium iteration solves the linear problem J&AU,
= R, _ , ,
(42)
where,
Kep= C'i,,C: E,,=E-E&w, is the tangent
(elastic-plastic)
R=&C'@
algorithmic
+AQ)=AF-C’AQ,
++g+jir$, stiffness
(430)
matrix, and (43c)
is the unbalanced nodal force vector. In Eqs. (43), all solution dependent quantities have final values at the end of the interior iterations that enforce the plasticity conditions. The procedure is standard and does not exhibit any peculiarities worth discussing.
88
A. Capsoni, L. Corradi
I Comput. Methods Appl. Mech. Engrg.
141 (1997)
81-93
5. Numerical examples The performances of the proposed model are demonstrated by solving some problems which are known as pathological. The first is the perfectly plastic double-edge-notched tensile specimen, for which a displacement approach is unable to predict the limit load [9]. Dimensions are indicated in Fig. 1 and the following values are assumed for the material constants: Young’s modulus E =70,Poisson ratio v=O.3 and uniaxial yield limit a, = 0.243. The accepted value of the limit load is FL = b(2 + a-)~,, / \15= 1.44. One quarter of the specimen is subdivided into 15 X5 square elements. 35 displacement controlled steps are performed in the interval 0~(EA)l(a,W)~20, where A is the relative displacement between the specimen ends. The computed load-displacement curve, illustrated in Fig. 2, approaches asymptotically the value FL. A second example is a ring of mean radius R = 10.5 and thickness h = 1.O, subject to equal and opposite forces (Fig. 3). This is a classical elastic benchmark case [lo], also studied as a plastic plane stress problem [6], and it is now analyzed in plane strain conditions. It can be considered a severe plastic test, since the structure exhibits an essentially flexural behavior and plastic strains are confined to limited regions about the two symmetry axes. By neglecting shear effects, the von Mises yield function is brought to the following form [ 1 l]
where N is the axial force and M the bending moment. The ring being statically determinate in terms of N, a mechanism involving plastic hinges in the four symmetry sections predicts the kinematic load
FL=~N'R(I-dl-
8M;l(N,R)*),
0
which can be regarded as a good estimate of the collapse value. For an uniaxial yield limit a,=40. FL=4.41 The remaining material constants are assumed as E = lo6 and v=O.25. L=30
Fig.
1. Double-edge-notched
specimen.
iii;,:;
0
Fig. 2. Dimensionless
4
load-displacement
a
12
curve for the perfectly
16
plastic double-edge-notched
specimen.
one gets
A. Capsoni, L. Corradi I Comput. Methods Appl. Mech. Engrg. 141 (1997) 81-93
Fig. 3. Plane strain circular
ring problem
89
and mesh adopted.
One quarter of the ring is subdivided into 16 X 2 equal elements (Fig. 3) and the displacement w at the load point is increased by performing 50 steps of amplitude Aw =0.1(~0R21Eh). As Fig. 4 shows, the load approaches the value: FL = 4.54, exceeding the comparison value by 3% owing to the distance of the Gauss points from the plastic hinge locations. The figure also depicts the load-displacement curves produced by the displacement approach, still unable to predict collapse, and by the so-called B-bar method [12], which overestimates the limit load by more than 20%. The third example, an axisymmetric plane strain cylinder of internal and external radii r, = 3.0 and re = 9.0, subjected to internal pressure p, was also suggested as a test for element evaluation in the incompressible range [13]. Here, the material is assumed as perfectly plastic, with E = 1.0, v=O.3, a, =0.2, and its response up to collapse is computed. In a sense, the situation is opposite to the previous one, in that plasticity progressively spreads and eventually extends over the entire solid. The analytical limit pressure is [14] pL = $gJn(,-,
/ri) = 1.269v0 = 0.254,
The mesh proposed in Ref. [13] was assumed, with five unequal elements covering a sector of 10 degrees, and results show that the computed limit load exceeds by only 0.3% the analytical value (Fig. 5). The displacement approach still fails: a pressure p = 1.07 pL is predicted for the internal displacement u = 20 ri(cro/E) and the P-U curve keeps growing afterwards, with a constant even if small slope. The problem known as Cook’s problem, a tapered solid clamped at one end and loaded at the other by a
F 12.00 ,
000 000
1.00
200
Fig. 4. Load-displacement
300
4M)
curve for the circular
500
ring.
wEh sR
A. Capsoni, L. Corradi I Comput. Methods
90
Fig. 5. Axisymmetric
Appl. Mech. Engrg. 141 (1997) 81-93
plane strain cylinder-Mesh
Fig. 6. Cook’s problem4X4
and load-displacement
and 8X8
curve.
meshes.
parabolically distributed vertical force, is considered as a final example (Fig. 6, where the geometry is defined). Linear isotropic and kinematic hardening are now included by assuming, as in Refs. [2,6], E = 70.,
u = l/3,
u0 = 0.243,
c = 0.135,
a = 0.015,
c and a being the constants appearing in Eqs. (20). 18 loading steps of amplitude AF =O. 1 are applied. Three meshes are examined, with 4, 8 and 16 elements on each side, and results are depicted in Fig. 7. The intermediate mesh appears adequate and very little is gained by its refinement. The displacement for F = 1.8 is very nearly the same as predicted in Ref. [2]. Comparison of the load-displacement curve obtained with those produced, with the same mesh, by the displacement and B-bar methods demonstrates significant improvements (Fig. 8).
6. Concluding
remarks
A mixed finite element model for plane strain elastic-plastic analysis was proposed and its effectiveness was demonstrated by solving a number of meaningful examples. One of the main goals of the formulation was to
A. Capsoni, L. Corradi
I Comput. Methods Appl. Mech. Engrg.
141 (1997)
91
81-93
F 2.00
t 0400
2 112s
1
0.00 -f 000
Fig. 7. Load-displ,xement
I
I
I
I
,
0.50
100
1.50
2.00
2.50
curve for Cook’s problem-Effects
ow_r 000
Fig. 8. Load-displacement
I
I
I
I
loo
150
200
1
050
of mesh refinement
curves for Cook’s problem-Comparison
with alternative
u
(displacements
formulations
,
for F= 1.8 are indicated)
u
2.50
(displacements
for F=
I .8 are indicated).
avoid locking when plasticity spreads and strains become nearly deviatoric, while retaining the same stress redistribution capability of a displacement model. These features were assessed in Ref. [l] and numerical tests indicate that they are actually incorporated. In fact, for all the perfectly plastic cases the collapse load was predicted with excellent accuracy, but slightly overestimated, indicating that the mixed nature of the formulation does not entail an excess of flexibility. To avoid locking, the well known device of supplementing the strain approximation with enhanced fields was used, but the proposed formulation provides some simplifications. The procedure followed is conceptually analogous to that in Ref. [6] where, however, the current elastic-plastic stiffness matrix is condensed at each stage of the solution process. Here, the operation is performed on the elastic stiffness and hardening matrices separately; for perfect plasticity (or linear hardening) they are both constant and can be condensed once and for all, with significant saving. The stress redistribution capability was preserved by expressing the stress and strain approximations in terms of the generalized variables dictated by the displacement model. With respect to existing proposals, they introduce only mino:r modifications. However, a formulation in terms of generalized variables permits the use of numerical devices borrowed from Mathematical Programming concepts. Due to space limitations, this aspect was only partially discussed; the way of dealing with the non-regular element yield surface, however, is an example of the effectiveness of these tools. To this purpose, it is mentioned that, even if displacement control is preferable for perfect plasticity (and was systematically employed in the examples), in principle load control could also be used. In this case, the collapse load is detected by a. test on matrix A,, defined by [ [ 11, Eq. (67)], again exploiting Mathematical Programming theory. Obviously, the load step amplitude must become very small when approaching the collapse limit.
92
A. Capsoni, L. Corradi
I Comput. Methods Appl. Mech. Engrg.
141 (1997)
81-93
Examples also indicate that the discrete scheme performs well with coarse meshes and is rather insensitive to mesh distortion, an expected occurrence, since the trial functions used are essentially existing ones exhibiting this feature. Finally, it is underlined that the formulation rests on a constitutive law which was fully reduced to the plane by exploiting the deviatoric nature of plastic flow associated to the von Mises condition. This makes the finite element model specifically oriented to a particular, even if relevant, class of problems, but the significant simplifications which become possible more than compensate for the limited range of applicability.
Acknowledgments This research was supported by the Italian Ministry (MURST).
Appendix
of University
and Scientific
and Technological
Research
A
Proof of properties (a), (b) and (c) in Section 4. The solution strategy crucially relies on the fact that the matrix o=z+fW,w,
(A.1)
is non-singular and this assertion is now proved. It is recalled that both M, and W are symmetric and at least positive semidefinite block-diagonal matrices, each 5 X5 block referring to a finite element. Therefore, only the element level needs to be considered. The property is easily proved when material kinematic hardening is present and, hence, W is strictly positive definite. By writing @= [W-’ +M,]W it is immediately recognized that 0, as the product of two symmetric and positive definite matrices, has a positive determinant. For perfectly plastic or purely isotropically hardening materials, matrix W is given by [ [l], Eq. (67b)]; namely 11000 11000
EL! >o.
0 0
0 0
If 0 were singular, condition becomes
0 0
0 0
(A.2)
(1 - 2V)2
0 0
a vector x #O would exist such that Ox = 0. When account is taken of Eqs. (A. l), (A.2), the
=O.
For this equation solution 1 + a(M,, a&f,,
(A.3)
to be satisfied by a vector x #O, the following
homogeneous
system must have a non trivial
+M,,) (A.4)
+ M,,)
where the components
of M,
are denoted
as M,j.The matrix
being
symmetric
and positive
semidefinite,
M,,=M,, and (ASa-c)
M ,,20, M,, a09 M,,M,,-M;, 20. The condition
that the determinant
of the coefficient
matrix in Eqs. (A.4) be equal to zero yields
93
A. Capsoni, L. Corradi I Comput. Methods Appl. Mech. Engrg. 14I (1997) 81-93
M,, +q2
(A.6)
+2h4,, = +o,
Eqs. (A.Sa,b) and (A.6) imply 2M,,
+M*,+2M,,)(M,,
taking account
I7=(M,,
and, hence, M,,+M,,
- 2M,, ~0, thus establishing
the inequality (A.7a)
+&-2M,,)
of Eq. (A.~c), one can also write
+ M,,)* -4M;,
3
(M,,+ M,2)2- 4M,, M,,=(M,,- M,,)'>O.
(A.7b)
The contradiction between Eqs. (A.7a), (A.7b) proves that the equation OX = 0 does not have solutions other than x =0 and, hence, that matrix 0 is non-singular. Property (a) having now been assessed, one can write 0 =Z +M,W= E ’ and the condition ss-’ = u-l L c=z E ==Z-EMu\c Taking
S=Z-MM,Ws,
account of the symmetry
M,S’
of both Ma and W, from Eqs. (A.8) one obtains E’W = W - &‘WM,W,
= M, - M,WM,=‘;
hence M, = [Ma + M, WM,]$, 0, =Mu = =[M,
64.8)
W= P’[W f WM, W] and
- Ma WM,s:“;
d2, = WEs’= E”[W-
WM,W]E
64.9)
The above relations show that &?, and az are symmetric, positive semidefinite matrices, since W and M, are individually so, and assesses property (b). This result also implies sMo =M,Et, WE= s:“‘W and from Eqs. (A.8) one obtains =?=I-M,z;“‘W, Together
with Eqs. l(A.8) the above relation
proves property
(c).
References [I] A. Capsoni and L. Corradi, A mixed finite element model for plane strain elastic-plastic analysis: Part I-Formulation and assessment of the overall behavior, Comp. Methods Appl. Mech. Engrg. 141 (1997) 67. [2] J.C. Simo and MS. Rifai, A class of assumed strain methods and the method of incompatible modes, Int. J. Numer. Methods Engrg. 29 (1990) 1595. [3] T.H.H. Pian and H. Sumihara, Rational approach for assumed stress finite elements, Int. .J. Numer. Methods Engrg. 20 (1984) 1682. [4] J.B. Martin, Plasticity: Theory and Applications (M.I.T. Press, Cambridge, MA, 1975). [5] J.C. Simo, J.G. Kennedy and S. Govindjee, Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms, Int. J. Numer. Methods Engrg. 26 (1988) 2161. [6] S.L. Weissman and M. Jamjian, Two-dimensional elastoplasticity-Approximation by mixed finite elements, Int. I. Numer. Methods Engrg. 36 (1993) 3703. [7] J.H. Argyris, Contitma and Discontinua, in: Proc. 1st Conf. Matrix Methods in Structural Mechanics, AFFDL, TR 66-88, Dayton, OH (1966) 11-92. [8] T. Belytschko, J.S. Ong, W.K. Liu and J.M. Kennedy, Hourglass control in linear and nonlinear problems, Comp. Methods Appl. Mech. Engrg. 43 (1984) 251. ]9] J.C. Nagtegaal, D.M. Parks and J.R. Rice, On numerically accurate finite element solutions in the fully plastic range, Comp. Methods Appl. Mech. Engrg. 4 (1974) 153. [IO] T. Belytschko and E. Bachrach, Efficient implementation of quadrilaterals with high coarse-mesh accuracy, Comp. Methods Appl. Mech. Engrg. 54 (1986) 279. [l I] W.F. Chen and D.J Han, Plasticity for structural engineers (Springer Verlag, 1988). [12] T.J.R. Hughes, Generalization of selective integration procedures to anisotropic and nonlinear media, Int. J. Numer. Methods Engrg. 15 (1980) 1413. [ 131 R.H. MacNeal and !R.L. Harder, A proposed standard set of problems to test finite element accuracy, Finite Element Analysis Design 1 (1985, 3. [14] W. Prager and Ph.G.Jr. Hodge, Theory of perfectly plastic solids (Wiley, New York, 195 I).