- Email: [email protected]
A highly accurate brick element based on a three-field variational principle for elasto-plastic analysis
Available online at www.sciencedirect.com
Finite Elements in Analysis and Design 39 (2003) 1155 – 1171 www.elsevier.com/locate/!nel
A highly accurate brick element based on a three-!eld variational principle for elasto-plastic analysis Yan P. Caoa , N. Hub;∗ , H. Fukunagac , J. Lud , Zhen H. Yaoa a
b
Department of Engineering Mechanics, Tsinghua University, Beijing 100084, People’s Republic of China Department of Mechanical Engineering, The Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA c Department of Aeronautics and Space Engineering, Tohoku University, 01 Aramaki-Aoba, Aoba-ku, Sendai 980-8579, Japan d LASMIS, Universite de technologie Troyes 12, rue Marie Curie-BP 2060-10010 Troyes Cedex France Received 28 August 2001; received in revised form 24 June 2002; accepted 2 July 2002
Abstract In this paper, a penalty equilibrating 3-D mixed element based on a modi!ed Hu–Washizu variational functional has been proposed to improve the accuracy of the elasto-plastic analysis for irregular meshes. In order to construct this element, a penalty term, which can enforce the stress components to satisfy the equilibrium equations in a weak form, is introduced into the Hu–Washizu three-!eld functional to “penalize” the false stress caused by the distorted meshes. The selection of the penalty factor has also been discussed in detail. Unlike traditional hybrid or mixed elements, the performance of this element is not so sensitive to mesh distortion. Furthermore, compared with hybrid stress elements, the proposed element based on the three-!eld variational principle is more suitable for material non-linear analysis. Numerical examples including elastic and elasto-plastic analyses have demonstrated the improved performance of the proposed element. ? 2002 Elsevier B.V. All rights reserved. Keywords: Three-!eld variational principle; Elasto-plastic analysis; Penalty equilibrium conditions; 3-D mixed element; Mesh distortion
1. Introduction A major goal of the element research in the FEM !eld has been to develop the low-order elements with high accuracy to meet two main requirements, i.e. avoiding locking behavior in bending ∗
Corresponding author. Tel.: +1-410-516-4398; fax: +1-410-516-4316. E-mail address: [email protected] (N. Hu).
0168-874X/03/$ - see front matter ? 2002 Elsevier B.V. All rights reserved. PII: S 0 1 6 8 - 8 7 4 X ( 0 2 ) 0 0 1 6 2 - 2
1156
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
dominated or incompressible problems and reduction of element sensitivity to mesh distortion. Generally, some important initial approaches fall into the so-called reduced integration clari!cation, which may introduce the so-called zero-energy modes into elements. Later on, to overcome the setback of the appearance of spurious zero-energy modes caused by usual reduced integration schemes, Kavanagh and Key [1] put forward the selective reduced integration (SRI), which actually fall within the concept of mixed !nite element methods as shown by Malkus and Hughes [2]. A general method— B-bar method, re!ned from the SRI scheme by Hughes [3], can be derived from the Hu–Washizu variational principle as illustrated by Simo [4]. Flanagan and Belytschko [5] employed a projection operator to illustrate that the usual !nite element approximation can be written as the sum of a de-coupled under-integrated stiOness matrix and a ‘stabilization’ stiOness matrix. Belytschko and Bachrach [6] further presented formulations based upon the Hu–Washizu functional, which lead to the well-known quintessential bending incompressible (QBI) element. However, this element is not frame invariant, i.e. it is dependent upon the coordinate system. Starting from another approach, Wilson et al. [7] introduced a set of incompatible displacements to improve the displacement element’s bending behavior. Taylor et al. [8] modi!ed the incompatible modes so that the distorted element can pass the constant strain/stress patch test. Using the stationary condition for potential energy functional of incompatible discrete system, Wu et al. [9] proposed an energy consistency requirement, which is just the patch test condition (PTC). Based on this condition, a general formulation to produce incompatible functions was presented. Recently Wu et al. [10] proposed a post-treatment approach, i.e., the revise-stiOness method, to construct incompatible elements. Numerical examples have shown that their method can be used to construct various incompatible models successfully. Pian and Sumihara [11] used incompatible displacements to generate an assumed stress !eld for a four-node plane stress element formulated via the Hellinger–Reissner [12] variational principle. In their assumed stress hybrid method, it is not necessary to impose on the equilibrium conditions a priori, but the homogeneous equilibrium equations are enforced to be satis!ed through the Lagrange multiplier terms, i.e., the internal incompatible displacements. This element exhibits excellent characteristics in bending applications. In addition, its sensitivity to mesh distortion from a parallelogram shape appears to be among the smallest of four-node elements. Pian and Tong [13] also proposed a 3-D element using the similar method. Originally, the formulation used in this element required a quadratic perturbation of the element shape [11], which does not possess a clear physical meaning in spite of its elegant mathematical derivations. Subsequently, Pian and Wu [14] presented an alternative approach to obtain the stress interpolation to avoid the perturbation. Wu and Cheung [15] further modi!ed the above 2-D hybrid element by considering the equilibrium conditions through the penalty functions. Rubinstein et al. [16] and Punch and Atluri [17] presented a theoretical framework for the selection of assumed stress !elds based on the Hellinger–Reissner functional. The resulting 2-D element can yield the comparable results to the Pian and Sumihara element [11]. Simo and Rifai [18] put forward four-node elements based upon the Hu–Washizu functional [19], which included ‘enhanced strains’ terms. This kind of elements produced almost identical results to those obtained by Pian and Sumihara [11]. Usually, to avoid locking behavior in the presence of internal constraints, the assumed stress and strain !elds must a priori satisfy these constraints point-wise. Based on the three-!eld Hu–Washizu variational principle, Weissman and Taylor [20] constructed these assumed !elds to possess the aforementioned properties, so that the assumed stress and strain !elds a priori
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1157
satisfy the homogeneous equilibrium equations in a weak sense. From this general theoretical frame, they further developed a mixed !nite element method for in-plane problems [21,22]. The obtained results are identical to those of Pian and Sumihara’s element [11]. Weissman [23] further presented an 8-node three-dimensional brick element using the same theory for elasto-plastic problems. Starting from the Weissman’s theoretical framework of the three-!eld variational, we proposed a 3-D brick element for elastic analysis by introducing a penalty term in the Hu–Washizu functional for enforcing the equilibrium equations to be satis!ed [24]. Recently, Piltner and Taylor [25] proposed a B-Bar element from the Hu–Washizu variational principle. Besides, there are also many other important researches in this !eld, which will not be mentioned for the sake of simplicity. In general, elements derived from the Hellinger–Reissner two-!eld formulation are more computationally eTcient and popular. The Hellinger–Reissner functional is based on the complementary energy function, in which the constitutive laws need to be modi!ed when the elements are used in non-linear materials. Unlike the Hellinger–Ressiner functional, the Hu–Washizu functional is formulated in terms of strain energy function, which is more attractive for non-linear materials. In this paper, based on our previous work [24] for elastic analysis, an 8-node brick element from the Hu– Washizu variational principle is proposed for elasto-plastic analysis. To increase the computational accuracy for irregular meshes, a penalty term is employed in a modi!ed Hu–Washizu functional to enforce the stresses to satisfy the equilibrium equations. This penalty term can eOectively reduce the false stress caused by the distorted meshes. Finally, a mixed element for elasto-plastic analysis is constructed from this functional. In this work, the eOects of the dimension of the penalty factor and a modi!cation scheme have also been discussed. Numerical examples have shown that, compared with the currently existing well-performed elements, the present element predicts better results when meshes are irregular. 2. Theory In this section, the theoretical description of the proposed element is presented, which contains four aspects, i.e., the general weak form of a three-!eld functional, the !nite element implementation from the functional, the iteration strategy for material non-linear analysis and the design of penalty factors. 2.1. Weak form of Hu–Washizu functional for plastic analysis The total free energy is calculated by the following standard three-!eld functional in the Hu– Washizu variational principle, H−W (u; U; ) = [W (U − Up ) + T (∇S u − U)] dV − EXT (u); (1)
where the independent variables are the stress , strain U and displacement u respectively, and ∇S is the symmetric part of the gradient operator. Moreover, Up , the plastic tensor, is given by rate equations, and EXT (u), the potential energy of the external loading, is described by EXT (u) = uT b dV + uT t d: (2)
@
1158
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
Also, the following assumption that U = Ue + Up ;
(3)
where Ue is the elastic strain tensor, is adopted here. Furthermore, W (U − Up ) in Eq. (1) is the stored energy function, which relates the second rank stress tensor, , to the strain tensor, U, as follows: = @ W (U − Up ):
(4)
Generally the strain energy function W is assumed to be quadratic as follows: T
W (Ue ) = 12 Ue C Ue ;
(5)
where the elastic tensor, C , i.e., a rank four tensor, is de!ned by C = @
W:
(6)
To further complete the constitutive law for elasto-plastic analysis, the evolution of the plastic strain and hardening parameters, which are referred to as the ‘Vow rule’ and ‘hardening law’ respectively, are usually needed. Readers can refer to Ref. [23] for further details regarding these issues. Starting from Eq. (1) for the state of iteration n + 1, the modi!cation version by Weissman [23] is !rst adopted here as follows: T H–W (un+1 ; Un+1 ; n+1 ) = [W (U − Up )n+1 + n+1 (∇S un+1 − Un+1 )
+ T @ W (U − Up )n+1 ] dV − EXT (un+1 ); T
(7)
where the appended term in the functional, i.e., @ W (U − Up ) with Lagrange multipliers , denotes the internal constraints, which constrains the enlarged assumed strain distribution, over the minimum required for stability, so that the added terms do not contribute to an energy-like quantity. The obtained element from the modi!ed Hu–Washizu functional in Eq. (7) can produce almost the same results as Pian and Tong’s hybrid element [13] for elastic analysis, which possesses the comparatively stable performance within the currently existing elements for mesh distortion. However, it was still found that this element is sensitive to the distorted meshes, especially the results of stresses. In the plastic analysis, it is known that the accuracy of stresses inVuences the loading history signi!cantly. Then, it is very important to increase the accuracy of stresses in the plastic analysis. Generally, in usual hybrid or mixed elements, to pass the constant strain/stress patch test, a transformation matrix for stress and strain from the natural – co-ordinates to the physical x–y co-ordinates [11,13,23] should be evaluated at the center of an element. However, this constant transform matrix at the element center certainly loses some important information about element geometry when the element shape is irregular. For instance, when using this constant transformation matrix, Pian and Wu [14] found that the assumed stress !eld might lose some typical geometric information when meshes are irregular. The irregular mesh will generate some so-called false stresses, which destroy the satisfaction of equilibrium equations. To address the phenomena found by Pian and Wu [14], Wu and Cheung [15] proposed a penalty equilibrating technique for hybrid stress elements and we extended this technique into mixed 3-D elements from the Hu–Washizu principle
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1159
[24]. For elastic analysis, it was found in our previous research that the accuracy of results can be improved eOectively through this penalty equilibrating technique. Here, our previously proposed technique for 3-D mixed elements in elastic analysis is adopted and developed. Then, the functional in the above equation is further modi!ed as follows by inserting a penalty term from equilibrium conditions, T H–W (un+1 ; Un+1 ; n+1 ) = [W (U − Up )n+1 + n+1 (∇S un+1 − Un+1 )
+ T @ W (U − Up )n+1 ] dV − EXT (un+1 ) 1 − (Dn+1 )T (Dn+1 ) dV; 2 V
(8)
where is the penalty parameter, and D is the partial diOerential operator of equilibrium equations without consideration of body force, which can be expressed as @ @ @ 0 0 0 @x @y @z @ @ @ 0 0 : (9) D= 0 @y @x @z @ @ @ 0 0 0 @z @y @x By considering the possible inaccuracy of the weak form balance of momentum due to irregular meshes and constant transformations, the physical meaning of the added penalty term is that the equilibrium equations in terms of stress is compulsively enforced through the penalty function. If the accuracy of the stresses can be enhanced through the introduced penalty term, it is hopeful that the accuracy of the displacements can also be improved. In Eq. (8), the penalty term without is not of energy unit obviously. To keep the energy unit consistent, should be of the unit: m4 =N , hence this is a very special “penalty spring”. Using the following increment form to replace un+1 , Un+1 and n+1 in the right-hand side of Eq. (8), un+1 = un + Wu;
(10a)
Un+1 = Un + WU;
(10b)
n+1 = n + W
(10c)
and performing the variational operation with respect to variables U, , u and , respectively, the following set of equations can be obtained, @W (U − Up )n+1 n+1 n+1T n+1 T p n+1 n+1 U − dV = 0; (11a) U + @
W (U − U ) U @U V T n+1 (∇S un+1 − Un+1 ) dV − (Dn+1 )T (Dn+1 ) dV = 0; (11b) V
1160
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
T
V
V
n+1 (∇S un+1 ) dV = fext ;
(11c)
T @ W (U − Up )n+1 dV = 0;
(11d)
where Eq. (11a) represents the weak form of the material constitution under the internal constraints, Eq. (11b) denotes the kinematics under the equilibrium constraints, Eq. (11c) is the weak form of balance of momentum and Eq. (11d) is the orthogonal conditions for internal constraints, respectively. 2.2. Finite element implementation From Eqs. (11), the !nite element procedure can be set up using the standard interpolation functions for four kinds of variables U, , u and . Assume the following interpolation functions for strains, displacements, stresses and Lagrange multipliers, WU = WU1 + WU2 = E 1 We1 + E 2 We2 ;
(12a)
Wu = N Wa;
(12b)
W = SWs;
(12c)
= E i z;
(12d)
where the strain contains the components of the compatible strain mode WU1 and the incompatible strain mode WU2 , or enlarged strain terms, respectively. A detailed description of these interpolation functions in Eq. (12) for 8-node brick element can be found from the literature [23,24]. The assumed strain interpolations are constructed so that the matrix E 1 contains minimum distribution required for stability. Also, the following requirement should be satis!ed from the standpoint of De Veubeke [26] to design the assumed strain !elds as follows: WU1 ∩ WU2 = ?; i:e: E 1 ∩ E 2 = ?:
(13)
From Eqs. (11d), (12a) and (12d), the enlarged part of strain, i.e. We2 can be expressed in terms of We1 as follows:
−1 iT alg 2 iT n iT alg 1 (14) E Cn E dV E @ W dV + E Cn E dV We1 ; We2 = − V
V
V
Cnalg
where is the standard algorithmic elasto-plastic tangent by referring to Weissman’s work [23] and other authors. From Eqs. (12a) and (14), the following interpolation for strains can be obtained: WU = E n We1 + E 2 V n ; where n
1
E =E −E Vn = −
V
2
V
E
iT
(15) Cnalg E 2
E iT Cnalg E 2 dV
dV
−1
−1 V
V
E iT Cnalg E 1 dV;
E iT @ W n dV:
(16a) (16b)
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1161
Then, from Eqs. (11a), (11b) and (11c), we can get the following three equations: T E n (@ W n + Cnalg E n We1 + Cnalg E 2 V n − n − SWs) dV = 0;
(17a)
V
V
S T (Ban + BWa − Un − E n We1 − E 2 V n ) dV
−
V
V
P T P dV sn +
V
P T P dV Ws = 0;
(17b)
B T (n + SWs) dV = fext ;
(17c)
where B is the strain displacement matrix, i.e. B = ∇S N and P matrix is equal to DS. Finally, from the above three equations, we can get: Hn We1 − ATn Ws = r1 ;
(18a)
G Wa − An We1 − QWs = r2 ;
(18b)
G T Ws = r3 ;
(18c)
where
Hn =
V
r1 =
V
E
nT
Cnalg E n
dV;
An =
V
T
n
S E dV;
G=
V
T
E n (−@ W n − Cnalg E 2 V n + n ) dV;
S T (−Ban + Un + E 2 V n ) dV + V B T n dV: r3 = fext − r2 =
V
T
S B dV;
Q=
V
P T P dV; (19a) (19b)
P T P dV sn ;
V
(19c) (19d)
From Eqs. (18a)–(18c), the element-level equilibrium equation for the displacement increment vector Wa can be expressed as follows: kne Wa = fne ;
(20)
where kne = G T (An Hn−1 ATn + Q)−1 G ;
(21a)
fne = r3 + G T (An Hn−1 ATn + Q)−1 (An Hn−1 r1 + r2 ):
(21b)
kne
fne ,
and the global system equation can be obtained through a standard assembling Finally using procedure and the nodal displacement increment can be solved during non-equilibrium iterations. A merit of the proposed technique is its simple numerical implementation. However, the computational cost needs to be addressed here. For elastic analysis, which only needs the stiOness matrix, it
1162
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
can be seen from Eq. (21a) that Weissman’s element [23] (without the penalty term Q) is more eT1 T −1 cient by merely obtaining the inverse matrix of A as follows: kne =G T (A− n Hn An )G . Then, one can avoid twice operations of inverse matrix in Eq. (21a). However, for elasto-plastic analysis, the strain !eld has to be recovered by considering that the radial return algorithm is strain-driven. Then, from Eq. (18a), the strain parameter is calculated using the following formulation: We1 =Hn−1 (r1 +ATn Ws). It can be seen, the inverse matrix of Hn has to be obtained regardless of the existence of Q. Therefore, in this case, the computational cost of the present element is similar to that of Weissman’s element. 2.3. Iteration strategy Summarizing the aforementioned algorithm in the previous sections, the sequence of calculations performed in iteration (n + 1) is described as follows: 1. Solve nodal displacement increment Wa from the global equation assembled from the elemental stiOness matrices and load vectors in Eqs. (20) and (21). 2. Update the nodal displacement using the obtained displacement increment as an+1 = an + Wa: 3. Check for the convergence of non-linear computation. The convergence criterion can be chosen for the displacement increment, or others. If the criterion is not satis!ed, go to the next step. Otherwise, jump out of the non-equilibrium iteration and enter into the next loading step. 4. Compute the increment of strain parameter We1 and the increment of stress parameter Ws from Eqs. (18a)–(18c). 5. Compute the increment of strain parameter We2 from Eq. (14). 6. Update the strain and stress parameter as follows: e1n+1 = e1n + We1 ;
e2n+1 = e2n + We2 ;
sn+1 = sn + Ws:
(23)
7. Compute the increment of strain and stress from Eqs. (12a) and (12c) using the interpolation matrices E 1 , E 2 and S. alg 8. Compute @ W (U − UP )n+1 and Cn+1 using a standard strain driven return mapping algorithm (see Weissman’s work [23]). 9. Evaluate elemental stiOness matrices and residual vectors with Eqs. (19)–(21), and assemble the global FEM equation for obtaining a new nodal displacement increment. 10. Set n = n + 1 and go back to Step 1. 2.4. Discussion on penalty factor In our previous work [24], we simply discussed the penalty factor in penalty equilibrating elements. However, there are still some problems unsolved, e.g., the eOects of the dimension of the penalty factor and the proper design of the penalty factor. In this section we will address these problems in detail.
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
E=1.5 ×105MPa, ν=0.25
1m/1000mm
1163
500N
1m/1000mm 2m/2000mm
500N A
y
500N x
z
500N 5m/5000mm
5m/5000mm
Fig. 1. Test of mesh distortion.
Table 1 Displacement of point A in the y-axis direction ( = =E and = 104 ) Displacement (m) 1:03 × 10−6 6:32 × 10−7 1:00 × 10−6
Unit:m Unit:mm Theory
First, we will analyze the dimension of the penalty factor, from Eqs. (19) and (21) it can be found that, in order to keep the energy unit consistent, should be of the unit: m4 =N . In all previous work about penalty equilibrating elements [15,24], was taken as (24) = ; E where = 104 still has the unit of m2 . For the test examples given in Refs. [15,24], the unit of does not cause problems. Therefore, Wu et al. [15] pointed out that the solutions of this penalty-equilibrium element were basically independent of the factor when it is larger than 103 , and then they de!ned = 104 . However, for practical applications, the unit must be considered. In order to show its eOect we discuss an example as shown in Fig. 1. In this example, a pure bending beam modeled with two elements is analyzed with two units: m and mm. Table 1 gives out the results of the penalty element proposed in Ref. [24]. It can be seen from Table 1 that the penalty element produces diOerent results with diOerent units. To overcome this shortfall, we de!ne the penalty factor as follows: =
l20 ; E
(25)
where is a large positive number without dimension, which is taken as =104 , and l0 is a function of the geometry of element, which is de!ned as l0 = 13 (l1 + l2 + l2 ):
(26)
1164
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
2
ζ
η
6 B
3
1
C O
z
5
7
4
A
ξ
x
y
8
Fig. 2. A typical solid element.
Table 2 Displacement of point A in the y-axis direction ( = l20 =E and = 104 ) Displacement (m) 1:04 × 10−6 1:04 × 10−6 1:00 × 10−6
Unit:m Unit:mm Theory
For a typical solid element in Fig. 2, l1 ; l2 ; l3 are de!ned as l1 = |OA| = a21 + b21 + c12 ; l2 = |OB| = a22 + b22 + c22 ; l3 = |OC| = a23 + b23 + c32 ;
(27)
where
a1
a2 a3
b1 b2 b3
c1
1
1 c2 = 8 1 c3 1
−1
−1
1
1
−1
−1
1
−1
−1
1
1
−1
1
1
1
−1
−1
−1
x 1 . −1 .. −1 x8 1
y1 .. . y8
z1 .. . : z8
(28)
The example in Fig. 1 is re-computed by using the penalty factor de!ned in Eq. (25), and the result is presented in Table 2. It can be seen from Table 2 that, using the penalty factor in Eq. (25), the penalty element produces consistent results with diOerent units.
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1165
3. Numerical examples In our previous research, some fundamental requirements for element construction have been checked for the proposed element. For instance, the ability for passing the constant stress/strain patch test was veri!ed. The eigenvalue analysis of the proposed element was carried out to show that the proposed element can overcome the locking at the nearly incompressible limit and does not possess the spurious zero energy modes. Furthermore, various benchmark problems for elastic analysis from the literature have been analyzed to illustrate the eOectiveness of the proposed element, such as Cook’s membrane problem, MacNeal and Harder’s beam problem [27] and a curved beam problem etc. It was found previously that the present element can predict the improved results; especially stresses compared with many other well behaved elements. In the current study, the performance of the proposed element for elasto-plastic analysis (only isotropic hardening model) is evaluated mainly with a number of challenging problems. Some representative elastic problems, which have not been considered in our previous research [24], are also employed. If no other particular statement, the penalty factors de!ned in Eq. (25) is employed. 3.1. Cook’s membrane problem: elasto-plastic analysis To demonstrate the superior performance of the presented element with coarse mesh grids, we !rst employ a tapered panel with unit thickness, which is clamped at one end and loaded by an in-plane bending load at the other end, as shown (plane view) in Fig. 3. This typical structure has been adopted for elastic analysis in our previous work [24] to illustrate the eOectiveness of our F=0.6
16 C
44
y 48 x Fig. 3. Cook’s membrane problem (plane view)—elasto-plastic plane strain simulation.
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171 Deflection at Point C in y-axis Direction
1166
22
20
18 Present Weissman's element (Ref. [23]) 16
14
0
2
4
6
8
10
12
14
16
18
Number of Elements Per Side
Fig. 4. Convergence of Cook’s membrane problem for the proposed element and Weissman’s element [23].
proposed element, and here we extend it into elasto-plastic case. The material constants of the panel are listed as: the Young’s modulus E =1:0, the Poisson’s ratio "= 13 , the uniaxial yield stress #y =0:1 and the tangential modulus in plastic state Ep = 0:01. The dimensions, load and boundary conditions of the panel are elaborated in Fig. 3. For the proposed element and Weissman’s element [23], the vertical displacements along the y-axis direction are shown in Fig. 4 with respect to the number of elements on one side. Fig. 4 illustrates that the reasonable results can be obtained using the present element even by very coarse mesh grids. Although both elements can get the converged results using about 16 elements on one edge, the result of the proposed element is obviously better than that of Weissman’s element [23] under coarse mesh. 3.2. Pinched hemispherical shell with 18◦ hole: elastic analysis The pinched hemispherical shell with an 18◦ hole, proposed by MacNeal and Harder [27], is further used to evaluate the performance of the proposed element for elastic analysis in doubly curved structures. The shell’s dimensions are: radius R = 10:0 and thickness t = 0:04. The material properties are: E = 6:875 × 107 and " = 0:3. The shell is subjected to alternating inward and outward forces 90◦ apart as shown in Fig. 5. According to symmetry boundary conditions, only one quarter of the shell is modelled by uniform meshes as illustrated in Fig. 5. Here, the eOect of the penalty factor has been investigated. In Eq. (25), has been taken as 104 , 105 and 106 , respectively. Fig. 6 presents the convergence of the displacement under the loads (normalized with 0.093, a value quoted in Ref. [28] as obtained by asymptotic expansion) for the proposed element. The results reported by Weissman [23] and those obtained by Hexa8 element (installed in NASTRAN code as an 8-node brick element) by MacNeal and Harder [27] are also plotted for comparison. Once again the present element exhibits obviously superior performance compared with other two outstanding elements. The eOect of the penalty factor mainly concentrates on the cases of coarse mesh. 3.3. ConAned plastic Bow This problem was suggested by Andel!nger et al. [29] as a challenging problem for the performance of element under con!ned plastic Vow. In this problem, a ‘regular brick’ (termed by
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1167
free
symmetry
symmetry
F=1
F=1 free
Normalized diflection at the load point
Fig. 5. Hemispherical shell with 18◦ hole—elastic analysis.
1.0
0.8
0.6 Present element ( κ =10 4 ) Present element ( κ =10 5 )
0.4
Present element ( κ =10 6 ) Weissman's element (Ref. [23]) HEXA8 (Ref. [27])
0.2
0.0 4
6
8
10
12
14
16
Number of elements per side
Fig. 6. Convergence of hemispherical shell problem for the proposed element, Weissman’s element [23] and Hexa8 element by McNeal and Harder [27].
Andel!nger et al. as ‘regular block’) with side length L = 100 and height h = 50 is !xed at its bottom, and a uniform pressure of q = 250 is applied to the center 20 × 20 area of its top surface. Using symmetry conditions, only one quarter is discretized by a uniform 5 ×5×5 mesh, as shown in Fig. 7. The material properties are listed as follows: E = 2:1 × 105 , " = 0:3, #y = 250:0 and Ep = 210:0. Fig. 8 presents the load–deVection curves obtained from the proposed element, Weissman’s element [23] and displacement-based element, respectively. In this case, the shape of the elements is regular, and then the results of the present element and those of Weissman’s element are completely identical. Compared with the results of the present element and Weissman’s element, the load–deVection curve for the displacement-based element exhibits a higher slope. The failure of the displacement-based element to achieve the correct slope is in agreement with the !ndings of Nagtegaal et al. [30] for the double-edged-notched tensile specimen, where this behaviour was observed for very !ne meshes.
1168
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
q=250 symmetry A
symmetry
x
y z
Fig. 7. A regular mesh consisting of 5 × 5 × 5 elements for checking element’s performance under con!ned plastic Vow.
8 7
Load Factor
6 5 4 3 2
Present & Weissman's element (Ref. [23]) Displacement-based element
1 0 0.0
0.5
1.0
1.5
2.0
2.5
Displacement at Center Point
Fig. 8. Load–displacement curves for various elements.
3.4. Clamped thin circular plate under uniform pressure In our previous research [24], a circular plate under a central concentrated load was analyzed elastically, and it was found that the accuracy of the plate’s stress and deVection was improved eTciently by using the proposed element. In this example, the same circular plate is further considered as an elasto-plastic problem. The dimensions and material properties are identical to those stated in Ref. [24] (see Fig. 9). For plastic analysis, the uniaxial yield stress #y is equal to 100.0 and the tangential modulus in plastic state Ep is set to be 0.0. The uniform pressure is 0.12 in this case. Unlike the mesh used in the elastic analysis, the two-layer discretization along the thickness direction is adopted for plastic analysis. Moreover, 3750 8-node incompatible 3-D brick elements in ABAQUS [31] commercial software are
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1169
Fig. 9. Mesh for a clamped circular plate under uniform pressure—elastic analysis and elasto-plastic analysis. Table 3 Central deVection of the plate Number of element
24 54 150 384 768
Central deVection of plate Present element
Weissman’s element (Ref. [23])
ABAQUS
1.111 1.260 1.621 1.860 2.280
1.001 1.230 1.580 1.831 2.270
2.441
employed to get the converged result in the plastic analysis for comparison. The !nal results are shown in Table 3. Inspection of this table reveals that the performance of the proposed element is also better than that of Weissman’s element [23] for elasto-plastic case. 4. Conclusions Based on the Hu–Washizu variational principle, a penalty equilibrating 8-node brick element is proposed in this paper. A technique in this element, by which the equilibrium equations in terms of stress components are enforced to be satis!ed through a penalty term in the functional, can reduce the false stresses caused by the elemental irregular shapes and improve the accuracy of analysis. The penalty factor has been carefully designed. Various numerical examples, such as elastic and elasto-plastic analyses, have shown that, compared with currently existing well-performed elements, the proposed element can predict the improved results when the elements are geometrically distorted.
1170
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
Acknowledgements The !nancial supports for this research from Framatome ANP. in France and Integrated Design Inc. in Japan are greatly appreciated. The project from Framatome is “Joint Research Project Agreement-No. 5” with the title: “Numerical investigation on the eOect of residual stresses on the mechanical behavior of 3-D structures” and the project from Integrated Design is: “Three-dimensional elasto-plastic analysis in civil engineering structures”. References [1] K.T. Kavanagh, S.W. Key, A note on selective and reduced integration techniques in !nite element method, Int. J. Numer. Methods Eng. 4 (1972) 148–150. [2] D.S. Malkus, T.J.R. Hughes, Mixed !nite element methods—reduced and selective integration techniques: a uni!cation of concepts, Comput. Methods Appl. Mech. Eng. 15 (1978) 63–81. [3] T.J.R. Hughes, Generalizing of selective integration procedures to anisotropic and nonlinear media, Int. J. Numer. Methods Eng. 15 (1980) 1413–1418. [4] J.C. Simo, R.L. Taylor, K.S. Pister, Variational and projection methods for the volume constraint in !nite deformation elasto-plasticity, Comput. Methods Appl. Mech. Eng. 51 (1985) 177–208. [5] D.P. Flanagan, T. Belytschko, A uniform strain hexahedron and quadrilaterals with high coarse-mesh accuracy, Int. J. Numer. Methods Eng. 17 (1981) 679–706. [6] T. Belytschko, W.E. Bachrach, ETcient implementation of quadrilaterals with high coarse-mesh accuracy, Comput. Methods Appl. Mech. Eng. 54 (1986) 279–301. [7] E.L. Wilson, R.L. Taylor, W.P. Doherty, J. Ghaboussi, Incompatible displacement modes, in: S.J. Fenves, N. Perrone, A.R. Rpbinson, W.C. Schnobrich (Eds.), Numerical and Computer Models in Structural Mechanics, Academic Press, New York, 1973. [8] R.L. Taylor, P.J. Beresford, E.L. Wilson, A nonconforming element for stress analysis, Int. J. Numer. Methods Eng. 10 (1976) 1211–1219. [9] C.C. Wu, M.G. Huang, T.H.H. Pian, Consistency condition and convergence criteria of incompatible elements: general formulation of incompatible functions and its application, Comput. Structures 27 (1987) 639–644. [10] C.C. Wu, Y.Q. Huang, E. Ramm, A further study on incompatible models: revise-stiOness approach and completeness of trial functions, Comput. Methods Appl. Mech. Eng. 190 (2001) 5923–5934. [11] T.H.H. Pian, K. Sumihara, Rational approach for assumed stress !nite elements, Int. J. Numer. Methods Eng. 20 (1984) 1685–1695. [12] E. Reissner, On a variational theorem in elasticity, J. Math. Phys. 29 (1950) 90–95. [13] T.H.H. Pian, P. Tong, Relations between incompatible displacement model and hybrid stress model, Int. J. Numer. Methods Eng. 22 (1986) 173–181. [14] T.H.H. Pian, C.C. Wu, A rational approach for choosing stress terms for hybrid !nite element formulations, Int. J. Numer. Methods Eng. 26 (1988) 2331–2343. [15] C.C. Wu, Y.K. Cheung, On optimization approaches of hybrid elements, Finite Elements Anal. Des. 21 (1995) 111–128. [16] R. Rubinstein, E.F. Punch, S.N. Atluri, An analysis of, and remedies for, kinematics modes in hybrid-stress elements, Comput. Methods Appl. Mech. Eng. 38 (1983) 63–92. [17] E.F. Punch, S.N. Atluri, Development and testing of stable, invariant isoparametric curvilinear 2- and 3-D hybrid stress elements, Comput. Methods Appl. Mech. Eng. 47 (1984) 331–356. [18] J.C. Simo, M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Methods Eng. 29 (1990) 1595–1638. [19] H.C. Hu, Variational Principle in Elasticity and its Application, Science Press, Beijing, 1981 (in Chinese). [20] S.L. Weissman, R.L. Taylor, Treatment of internal constraints by mixed !nite element methods: uni!cation of concepts, Int. J. Numer. Methods Eng. 33 (1992) 131–141.
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1171
[21] S.L. Weissman, R.L. Taylor, A uni!ed approach to mixed !nite element methods: application to in-plane problems, Comput. Methods Appl. Mech. Eng. 98 (1992) 127–151. [22] S.L. Weissman, A mixed formulation of nonlinear-elastic problems, Comput. Structures 44 (1992) 813–822. [23] S.L. Weissman, High-accuracy low-order three-dimensional brick elements, Int. J. Numer. Methods Eng. 39 (1996) 2337–2361. [24] Y.P. Cao, N. Hu, J. Lu, H. Fukunaga, Z.H. Yao, A 3-D brick element based on Hu–Washizu variational principle for mesh distortion, Int. J. Numer. Methods Eng. 53 (2002) 2529–2548. [25] R. Piltner, R.L. Taylor, A systematic construction of B-Bar functions for linear and Non-linear mixed-enhanced !nite elements for plane elasticity problems, Int. J. Numer. Methods Eng. 44 (1999) 615–639. [26] B. Fraeijs De Veubeke, Displacement and equilibrium models in the !nite element methods, in: O.C. Zienkiewicz, G.S. Holister (Eds.), Stress Analysis, Wiley, London, 1965. [27] R.H. MacNeal, R.L. Harder, A proposed standard set of problems to test !nite element accuracy, Finite Elements Anal. Des. 1 (1985) 3–20. [28] J.C. Simo, D.D. Fox, M.S. Rifai, On a stress resultant geometrically exact shell model, Part II: The linear theory; computational aspects, Comput. Methods Appl. Mech. Eng. 73 (1989) 53–92. [29] U. Andel!nger, E. Ramm, D. Roehl, 2D- and 3D-enhanced assumed strain elements and their application in plasticity, Proceedings of COMPLAST III, Barcelona, Spain, 1992. [30] J.C. Nagtegaal, D.M. Parks, J.C. Rice, On numerically accurate !nite element solutions in the fully plastic range, Comput. Methods Appl. Mech. Eng. 4 (1974) 153–177. [31] HKS ABAQUS/Standard, Version 5.8, User’s Manual, Hibbit, Karlsson and Sorensen Inc., USA, 1998.
Finite Elements in Analysis and Design 39 (2003) 1155 – 1171 www.elsevier.com/locate/!nel
A highly accurate brick element based on a three-!eld variational principle for elasto-plastic analysis Yan P. Caoa , N. Hub;∗ , H. Fukunagac , J. Lud , Zhen H. Yaoa a
b
Department of Engineering Mechanics, Tsinghua University, Beijing 100084, People’s Republic of China Department of Mechanical Engineering, The Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA c Department of Aeronautics and Space Engineering, Tohoku University, 01 Aramaki-Aoba, Aoba-ku, Sendai 980-8579, Japan d LASMIS, Universite de technologie Troyes 12, rue Marie Curie-BP 2060-10010 Troyes Cedex France Received 28 August 2001; received in revised form 24 June 2002; accepted 2 July 2002
Abstract In this paper, a penalty equilibrating 3-D mixed element based on a modi!ed Hu–Washizu variational functional has been proposed to improve the accuracy of the elasto-plastic analysis for irregular meshes. In order to construct this element, a penalty term, which can enforce the stress components to satisfy the equilibrium equations in a weak form, is introduced into the Hu–Washizu three-!eld functional to “penalize” the false stress caused by the distorted meshes. The selection of the penalty factor has also been discussed in detail. Unlike traditional hybrid or mixed elements, the performance of this element is not so sensitive to mesh distortion. Furthermore, compared with hybrid stress elements, the proposed element based on the three-!eld variational principle is more suitable for material non-linear analysis. Numerical examples including elastic and elasto-plastic analyses have demonstrated the improved performance of the proposed element. ? 2002 Elsevier B.V. All rights reserved. Keywords: Three-!eld variational principle; Elasto-plastic analysis; Penalty equilibrium conditions; 3-D mixed element; Mesh distortion
1. Introduction A major goal of the element research in the FEM !eld has been to develop the low-order elements with high accuracy to meet two main requirements, i.e. avoiding locking behavior in bending ∗
Corresponding author. Tel.: +1-410-516-4398; fax: +1-410-516-4316. E-mail address: [email protected] (N. Hu).
0168-874X/03/$ - see front matter ? 2002 Elsevier B.V. All rights reserved. PII: S 0 1 6 8 - 8 7 4 X ( 0 2 ) 0 0 1 6 2 - 2
1156
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
dominated or incompressible problems and reduction of element sensitivity to mesh distortion. Generally, some important initial approaches fall into the so-called reduced integration clari!cation, which may introduce the so-called zero-energy modes into elements. Later on, to overcome the setback of the appearance of spurious zero-energy modes caused by usual reduced integration schemes, Kavanagh and Key [1] put forward the selective reduced integration (SRI), which actually fall within the concept of mixed !nite element methods as shown by Malkus and Hughes [2]. A general method— B-bar method, re!ned from the SRI scheme by Hughes [3], can be derived from the Hu–Washizu variational principle as illustrated by Simo [4]. Flanagan and Belytschko [5] employed a projection operator to illustrate that the usual !nite element approximation can be written as the sum of a de-coupled under-integrated stiOness matrix and a ‘stabilization’ stiOness matrix. Belytschko and Bachrach [6] further presented formulations based upon the Hu–Washizu functional, which lead to the well-known quintessential bending incompressible (QBI) element. However, this element is not frame invariant, i.e. it is dependent upon the coordinate system. Starting from another approach, Wilson et al. [7] introduced a set of incompatible displacements to improve the displacement element’s bending behavior. Taylor et al. [8] modi!ed the incompatible modes so that the distorted element can pass the constant strain/stress patch test. Using the stationary condition for potential energy functional of incompatible discrete system, Wu et al. [9] proposed an energy consistency requirement, which is just the patch test condition (PTC). Based on this condition, a general formulation to produce incompatible functions was presented. Recently Wu et al. [10] proposed a post-treatment approach, i.e., the revise-stiOness method, to construct incompatible elements. Numerical examples have shown that their method can be used to construct various incompatible models successfully. Pian and Sumihara [11] used incompatible displacements to generate an assumed stress !eld for a four-node plane stress element formulated via the Hellinger–Reissner [12] variational principle. In their assumed stress hybrid method, it is not necessary to impose on the equilibrium conditions a priori, but the homogeneous equilibrium equations are enforced to be satis!ed through the Lagrange multiplier terms, i.e., the internal incompatible displacements. This element exhibits excellent characteristics in bending applications. In addition, its sensitivity to mesh distortion from a parallelogram shape appears to be among the smallest of four-node elements. Pian and Tong [13] also proposed a 3-D element using the similar method. Originally, the formulation used in this element required a quadratic perturbation of the element shape [11], which does not possess a clear physical meaning in spite of its elegant mathematical derivations. Subsequently, Pian and Wu [14] presented an alternative approach to obtain the stress interpolation to avoid the perturbation. Wu and Cheung [15] further modi!ed the above 2-D hybrid element by considering the equilibrium conditions through the penalty functions. Rubinstein et al. [16] and Punch and Atluri [17] presented a theoretical framework for the selection of assumed stress !elds based on the Hellinger–Reissner functional. The resulting 2-D element can yield the comparable results to the Pian and Sumihara element [11]. Simo and Rifai [18] put forward four-node elements based upon the Hu–Washizu functional [19], which included ‘enhanced strains’ terms. This kind of elements produced almost identical results to those obtained by Pian and Sumihara [11]. Usually, to avoid locking behavior in the presence of internal constraints, the assumed stress and strain !elds must a priori satisfy these constraints point-wise. Based on the three-!eld Hu–Washizu variational principle, Weissman and Taylor [20] constructed these assumed !elds to possess the aforementioned properties, so that the assumed stress and strain !elds a priori
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1157
satisfy the homogeneous equilibrium equations in a weak sense. From this general theoretical frame, they further developed a mixed !nite element method for in-plane problems [21,22]. The obtained results are identical to those of Pian and Sumihara’s element [11]. Weissman [23] further presented an 8-node three-dimensional brick element using the same theory for elasto-plastic problems. Starting from the Weissman’s theoretical framework of the three-!eld variational, we proposed a 3-D brick element for elastic analysis by introducing a penalty term in the Hu–Washizu functional for enforcing the equilibrium equations to be satis!ed [24]. Recently, Piltner and Taylor [25] proposed a B-Bar element from the Hu–Washizu variational principle. Besides, there are also many other important researches in this !eld, which will not be mentioned for the sake of simplicity. In general, elements derived from the Hellinger–Reissner two-!eld formulation are more computationally eTcient and popular. The Hellinger–Reissner functional is based on the complementary energy function, in which the constitutive laws need to be modi!ed when the elements are used in non-linear materials. Unlike the Hellinger–Ressiner functional, the Hu–Washizu functional is formulated in terms of strain energy function, which is more attractive for non-linear materials. In this paper, based on our previous work [24] for elastic analysis, an 8-node brick element from the Hu– Washizu variational principle is proposed for elasto-plastic analysis. To increase the computational accuracy for irregular meshes, a penalty term is employed in a modi!ed Hu–Washizu functional to enforce the stresses to satisfy the equilibrium equations. This penalty term can eOectively reduce the false stress caused by the distorted meshes. Finally, a mixed element for elasto-plastic analysis is constructed from this functional. In this work, the eOects of the dimension of the penalty factor and a modi!cation scheme have also been discussed. Numerical examples have shown that, compared with the currently existing well-performed elements, the present element predicts better results when meshes are irregular. 2. Theory In this section, the theoretical description of the proposed element is presented, which contains four aspects, i.e., the general weak form of a three-!eld functional, the !nite element implementation from the functional, the iteration strategy for material non-linear analysis and the design of penalty factors. 2.1. Weak form of Hu–Washizu functional for plastic analysis The total free energy is calculated by the following standard three-!eld functional in the Hu– Washizu variational principle, H−W (u; U; ) = [W (U − Up ) + T (∇S u − U)] dV − EXT (u); (1)
where the independent variables are the stress , strain U and displacement u respectively, and ∇S is the symmetric part of the gradient operator. Moreover, Up , the plastic tensor, is given by rate equations, and EXT (u), the potential energy of the external loading, is described by EXT (u) = uT b dV + uT t d: (2)
@
1158
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
Also, the following assumption that U = Ue + Up ;
(3)
where Ue is the elastic strain tensor, is adopted here. Furthermore, W (U − Up ) in Eq. (1) is the stored energy function, which relates the second rank stress tensor, , to the strain tensor, U, as follows: = @ W (U − Up ):
(4)
Generally the strain energy function W is assumed to be quadratic as follows: T
W (Ue ) = 12 Ue C Ue ;
(5)
where the elastic tensor, C , i.e., a rank four tensor, is de!ned by C = @
W:
(6)
To further complete the constitutive law for elasto-plastic analysis, the evolution of the plastic strain and hardening parameters, which are referred to as the ‘Vow rule’ and ‘hardening law’ respectively, are usually needed. Readers can refer to Ref. [23] for further details regarding these issues. Starting from Eq. (1) for the state of iteration n + 1, the modi!cation version by Weissman [23] is !rst adopted here as follows: T H–W (un+1 ; Un+1 ; n+1 ) = [W (U − Up )n+1 + n+1 (∇S un+1 − Un+1 )
+ T @ W (U − Up )n+1 ] dV − EXT (un+1 ); T
(7)
where the appended term in the functional, i.e., @ W (U − Up ) with Lagrange multipliers , denotes the internal constraints, which constrains the enlarged assumed strain distribution, over the minimum required for stability, so that the added terms do not contribute to an energy-like quantity. The obtained element from the modi!ed Hu–Washizu functional in Eq. (7) can produce almost the same results as Pian and Tong’s hybrid element [13] for elastic analysis, which possesses the comparatively stable performance within the currently existing elements for mesh distortion. However, it was still found that this element is sensitive to the distorted meshes, especially the results of stresses. In the plastic analysis, it is known that the accuracy of stresses inVuences the loading history signi!cantly. Then, it is very important to increase the accuracy of stresses in the plastic analysis. Generally, in usual hybrid or mixed elements, to pass the constant strain/stress patch test, a transformation matrix for stress and strain from the natural – co-ordinates to the physical x–y co-ordinates [11,13,23] should be evaluated at the center of an element. However, this constant transform matrix at the element center certainly loses some important information about element geometry when the element shape is irregular. For instance, when using this constant transformation matrix, Pian and Wu [14] found that the assumed stress !eld might lose some typical geometric information when meshes are irregular. The irregular mesh will generate some so-called false stresses, which destroy the satisfaction of equilibrium equations. To address the phenomena found by Pian and Wu [14], Wu and Cheung [15] proposed a penalty equilibrating technique for hybrid stress elements and we extended this technique into mixed 3-D elements from the Hu–Washizu principle
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1159
[24]. For elastic analysis, it was found in our previous research that the accuracy of results can be improved eOectively through this penalty equilibrating technique. Here, our previously proposed technique for 3-D mixed elements in elastic analysis is adopted and developed. Then, the functional in the above equation is further modi!ed as follows by inserting a penalty term from equilibrium conditions, T H–W (un+1 ; Un+1 ; n+1 ) = [W (U − Up )n+1 + n+1 (∇S un+1 − Un+1 )
+ T @ W (U − Up )n+1 ] dV − EXT (un+1 ) 1 − (Dn+1 )T (Dn+1 ) dV; 2 V
(8)
where is the penalty parameter, and D is the partial diOerential operator of equilibrium equations without consideration of body force, which can be expressed as @ @ @ 0 0 0 @x @y @z @ @ @ 0 0 : (9) D= 0 @y @x @z @ @ @ 0 0 0 @z @y @x By considering the possible inaccuracy of the weak form balance of momentum due to irregular meshes and constant transformations, the physical meaning of the added penalty term is that the equilibrium equations in terms of stress is compulsively enforced through the penalty function. If the accuracy of the stresses can be enhanced through the introduced penalty term, it is hopeful that the accuracy of the displacements can also be improved. In Eq. (8), the penalty term without is not of energy unit obviously. To keep the energy unit consistent, should be of the unit: m4 =N , hence this is a very special “penalty spring”. Using the following increment form to replace un+1 , Un+1 and n+1 in the right-hand side of Eq. (8), un+1 = un + Wu;
(10a)
Un+1 = Un + WU;
(10b)
n+1 = n + W
(10c)
and performing the variational operation with respect to variables U, , u and , respectively, the following set of equations can be obtained, @W (U − Up )n+1 n+1 n+1T n+1 T p n+1 n+1 U − dV = 0; (11a) U + @
W (U − U ) U @U V T n+1 (∇S un+1 − Un+1 ) dV − (Dn+1 )T (Dn+1 ) dV = 0; (11b) V
1160
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
T
V
V
n+1 (∇S un+1 ) dV = fext ;
(11c)
T @ W (U − Up )n+1 dV = 0;
(11d)
where Eq. (11a) represents the weak form of the material constitution under the internal constraints, Eq. (11b) denotes the kinematics under the equilibrium constraints, Eq. (11c) is the weak form of balance of momentum and Eq. (11d) is the orthogonal conditions for internal constraints, respectively. 2.2. Finite element implementation From Eqs. (11), the !nite element procedure can be set up using the standard interpolation functions for four kinds of variables U, , u and . Assume the following interpolation functions for strains, displacements, stresses and Lagrange multipliers, WU = WU1 + WU2 = E 1 We1 + E 2 We2 ;
(12a)
Wu = N Wa;
(12b)
W = SWs;
(12c)
= E i z;
(12d)
where the strain contains the components of the compatible strain mode WU1 and the incompatible strain mode WU2 , or enlarged strain terms, respectively. A detailed description of these interpolation functions in Eq. (12) for 8-node brick element can be found from the literature [23,24]. The assumed strain interpolations are constructed so that the matrix E 1 contains minimum distribution required for stability. Also, the following requirement should be satis!ed from the standpoint of De Veubeke [26] to design the assumed strain !elds as follows: WU1 ∩ WU2 = ?; i:e: E 1 ∩ E 2 = ?:
(13)
From Eqs. (11d), (12a) and (12d), the enlarged part of strain, i.e. We2 can be expressed in terms of We1 as follows:
−1 iT alg 2 iT n iT alg 1 (14) E Cn E dV E @ W dV + E Cn E dV We1 ; We2 = − V
V
V
Cnalg
where is the standard algorithmic elasto-plastic tangent by referring to Weissman’s work [23] and other authors. From Eqs. (12a) and (14), the following interpolation for strains can be obtained: WU = E n We1 + E 2 V n ; where n
1
E =E −E Vn = −
V
2
V
E
iT
(15) Cnalg E 2
E iT Cnalg E 2 dV
dV
−1
−1 V
V
E iT Cnalg E 1 dV;
E iT @ W n dV:
(16a) (16b)
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1161
Then, from Eqs. (11a), (11b) and (11c), we can get the following three equations: T E n (@ W n + Cnalg E n We1 + Cnalg E 2 V n − n − SWs) dV = 0;
(17a)
V
V
S T (Ban + BWa − Un − E n We1 − E 2 V n ) dV
−
V
V
P T P dV sn +
V
P T P dV Ws = 0;
(17b)
B T (n + SWs) dV = fext ;
(17c)
where B is the strain displacement matrix, i.e. B = ∇S N and P matrix is equal to DS. Finally, from the above three equations, we can get: Hn We1 − ATn Ws = r1 ;
(18a)
G Wa − An We1 − QWs = r2 ;
(18b)
G T Ws = r3 ;
(18c)
where
Hn =
V
r1 =
V
E
nT
Cnalg E n
dV;
An =
V
T
n
S E dV;
G=
V
T
E n (−@ W n − Cnalg E 2 V n + n ) dV;
S T (−Ban + Un + E 2 V n ) dV + V B T n dV: r3 = fext − r2 =
V
T
S B dV;
Q=
V
P T P dV; (19a) (19b)
P T P dV sn ;
V
(19c) (19d)
From Eqs. (18a)–(18c), the element-level equilibrium equation for the displacement increment vector Wa can be expressed as follows: kne Wa = fne ;
(20)
where kne = G T (An Hn−1 ATn + Q)−1 G ;
(21a)
fne = r3 + G T (An Hn−1 ATn + Q)−1 (An Hn−1 r1 + r2 ):
(21b)
kne
fne ,
and the global system equation can be obtained through a standard assembling Finally using procedure and the nodal displacement increment can be solved during non-equilibrium iterations. A merit of the proposed technique is its simple numerical implementation. However, the computational cost needs to be addressed here. For elastic analysis, which only needs the stiOness matrix, it
1162
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
can be seen from Eq. (21a) that Weissman’s element [23] (without the penalty term Q) is more eT1 T −1 cient by merely obtaining the inverse matrix of A as follows: kne =G T (A− n Hn An )G . Then, one can avoid twice operations of inverse matrix in Eq. (21a). However, for elasto-plastic analysis, the strain !eld has to be recovered by considering that the radial return algorithm is strain-driven. Then, from Eq. (18a), the strain parameter is calculated using the following formulation: We1 =Hn−1 (r1 +ATn Ws). It can be seen, the inverse matrix of Hn has to be obtained regardless of the existence of Q. Therefore, in this case, the computational cost of the present element is similar to that of Weissman’s element. 2.3. Iteration strategy Summarizing the aforementioned algorithm in the previous sections, the sequence of calculations performed in iteration (n + 1) is described as follows: 1. Solve nodal displacement increment Wa from the global equation assembled from the elemental stiOness matrices and load vectors in Eqs. (20) and (21). 2. Update the nodal displacement using the obtained displacement increment as an+1 = an + Wa: 3. Check for the convergence of non-linear computation. The convergence criterion can be chosen for the displacement increment, or others. If the criterion is not satis!ed, go to the next step. Otherwise, jump out of the non-equilibrium iteration and enter into the next loading step. 4. Compute the increment of strain parameter We1 and the increment of stress parameter Ws from Eqs. (18a)–(18c). 5. Compute the increment of strain parameter We2 from Eq. (14). 6. Update the strain and stress parameter as follows: e1n+1 = e1n + We1 ;
e2n+1 = e2n + We2 ;
sn+1 = sn + Ws:
(23)
7. Compute the increment of strain and stress from Eqs. (12a) and (12c) using the interpolation matrices E 1 , E 2 and S. alg 8. Compute @ W (U − UP )n+1 and Cn+1 using a standard strain driven return mapping algorithm (see Weissman’s work [23]). 9. Evaluate elemental stiOness matrices and residual vectors with Eqs. (19)–(21), and assemble the global FEM equation for obtaining a new nodal displacement increment. 10. Set n = n + 1 and go back to Step 1. 2.4. Discussion on penalty factor In our previous work [24], we simply discussed the penalty factor in penalty equilibrating elements. However, there are still some problems unsolved, e.g., the eOects of the dimension of the penalty factor and the proper design of the penalty factor. In this section we will address these problems in detail.
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
E=1.5 ×105MPa, ν=0.25
1m/1000mm
1163
500N
1m/1000mm 2m/2000mm
500N A
y
500N x
z
500N 5m/5000mm
5m/5000mm
Fig. 1. Test of mesh distortion.
Table 1 Displacement of point A in the y-axis direction ( = =E and = 104 ) Displacement (m) 1:03 × 10−6 6:32 × 10−7 1:00 × 10−6
Unit:m Unit:mm Theory
First, we will analyze the dimension of the penalty factor, from Eqs. (19) and (21) it can be found that, in order to keep the energy unit consistent, should be of the unit: m4 =N . In all previous work about penalty equilibrating elements [15,24], was taken as (24) = ; E where = 104 still has the unit of m2 . For the test examples given in Refs. [15,24], the unit of does not cause problems. Therefore, Wu et al. [15] pointed out that the solutions of this penalty-equilibrium element were basically independent of the factor when it is larger than 103 , and then they de!ned = 104 . However, for practical applications, the unit must be considered. In order to show its eOect we discuss an example as shown in Fig. 1. In this example, a pure bending beam modeled with two elements is analyzed with two units: m and mm. Table 1 gives out the results of the penalty element proposed in Ref. [24]. It can be seen from Table 1 that the penalty element produces diOerent results with diOerent units. To overcome this shortfall, we de!ne the penalty factor as follows: =
l20 ; E
(25)
where is a large positive number without dimension, which is taken as =104 , and l0 is a function of the geometry of element, which is de!ned as l0 = 13 (l1 + l2 + l2 ):
(26)
1164
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
2
ζ
η
6 B
3
1
C O
z
5
7
4
A
ξ
x
y
8
Fig. 2. A typical solid element.
Table 2 Displacement of point A in the y-axis direction ( = l20 =E and = 104 ) Displacement (m) 1:04 × 10−6 1:04 × 10−6 1:00 × 10−6
Unit:m Unit:mm Theory
For a typical solid element in Fig. 2, l1 ; l2 ; l3 are de!ned as l1 = |OA| = a21 + b21 + c12 ; l2 = |OB| = a22 + b22 + c22 ; l3 = |OC| = a23 + b23 + c32 ;
(27)
where
a1
a2 a3
b1 b2 b3
c1
1
1 c2 = 8 1 c3 1
−1
−1
1
1
−1
−1
1
−1
−1
1
1
−1
1
1
1
−1
−1
−1
x 1 . −1 .. −1 x8 1
y1 .. . y8
z1 .. . : z8
(28)
The example in Fig. 1 is re-computed by using the penalty factor de!ned in Eq. (25), and the result is presented in Table 2. It can be seen from Table 2 that, using the penalty factor in Eq. (25), the penalty element produces consistent results with diOerent units.
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1165
3. Numerical examples In our previous research, some fundamental requirements for element construction have been checked for the proposed element. For instance, the ability for passing the constant stress/strain patch test was veri!ed. The eigenvalue analysis of the proposed element was carried out to show that the proposed element can overcome the locking at the nearly incompressible limit and does not possess the spurious zero energy modes. Furthermore, various benchmark problems for elastic analysis from the literature have been analyzed to illustrate the eOectiveness of the proposed element, such as Cook’s membrane problem, MacNeal and Harder’s beam problem [27] and a curved beam problem etc. It was found previously that the present element can predict the improved results; especially stresses compared with many other well behaved elements. In the current study, the performance of the proposed element for elasto-plastic analysis (only isotropic hardening model) is evaluated mainly with a number of challenging problems. Some representative elastic problems, which have not been considered in our previous research [24], are also employed. If no other particular statement, the penalty factors de!ned in Eq. (25) is employed. 3.1. Cook’s membrane problem: elasto-plastic analysis To demonstrate the superior performance of the presented element with coarse mesh grids, we !rst employ a tapered panel with unit thickness, which is clamped at one end and loaded by an in-plane bending load at the other end, as shown (plane view) in Fig. 3. This typical structure has been adopted for elastic analysis in our previous work [24] to illustrate the eOectiveness of our F=0.6
16 C
44
y 48 x Fig. 3. Cook’s membrane problem (plane view)—elasto-plastic plane strain simulation.
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171 Deflection at Point C in y-axis Direction
1166
22
20
18 Present Weissman's element (Ref. [23]) 16
14
0
2
4
6
8
10
12
14
16
18
Number of Elements Per Side
Fig. 4. Convergence of Cook’s membrane problem for the proposed element and Weissman’s element [23].
proposed element, and here we extend it into elasto-plastic case. The material constants of the panel are listed as: the Young’s modulus E =1:0, the Poisson’s ratio "= 13 , the uniaxial yield stress #y =0:1 and the tangential modulus in plastic state Ep = 0:01. The dimensions, load and boundary conditions of the panel are elaborated in Fig. 3. For the proposed element and Weissman’s element [23], the vertical displacements along the y-axis direction are shown in Fig. 4 with respect to the number of elements on one side. Fig. 4 illustrates that the reasonable results can be obtained using the present element even by very coarse mesh grids. Although both elements can get the converged results using about 16 elements on one edge, the result of the proposed element is obviously better than that of Weissman’s element [23] under coarse mesh. 3.2. Pinched hemispherical shell with 18◦ hole: elastic analysis The pinched hemispherical shell with an 18◦ hole, proposed by MacNeal and Harder [27], is further used to evaluate the performance of the proposed element for elastic analysis in doubly curved structures. The shell’s dimensions are: radius R = 10:0 and thickness t = 0:04. The material properties are: E = 6:875 × 107 and " = 0:3. The shell is subjected to alternating inward and outward forces 90◦ apart as shown in Fig. 5. According to symmetry boundary conditions, only one quarter of the shell is modelled by uniform meshes as illustrated in Fig. 5. Here, the eOect of the penalty factor has been investigated. In Eq. (25), has been taken as 104 , 105 and 106 , respectively. Fig. 6 presents the convergence of the displacement under the loads (normalized with 0.093, a value quoted in Ref. [28] as obtained by asymptotic expansion) for the proposed element. The results reported by Weissman [23] and those obtained by Hexa8 element (installed in NASTRAN code as an 8-node brick element) by MacNeal and Harder [27] are also plotted for comparison. Once again the present element exhibits obviously superior performance compared with other two outstanding elements. The eOect of the penalty factor mainly concentrates on the cases of coarse mesh. 3.3. ConAned plastic Bow This problem was suggested by Andel!nger et al. [29] as a challenging problem for the performance of element under con!ned plastic Vow. In this problem, a ‘regular brick’ (termed by
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1167
free
symmetry
symmetry
F=1
F=1 free
Normalized diflection at the load point
Fig. 5. Hemispherical shell with 18◦ hole—elastic analysis.
1.0
0.8
0.6 Present element ( κ =10 4 ) Present element ( κ =10 5 )
0.4
Present element ( κ =10 6 ) Weissman's element (Ref. [23]) HEXA8 (Ref. [27])
0.2
0.0 4
6
8
10
12
14
16
Number of elements per side
Fig. 6. Convergence of hemispherical shell problem for the proposed element, Weissman’s element [23] and Hexa8 element by McNeal and Harder [27].
Andel!nger et al. as ‘regular block’) with side length L = 100 and height h = 50 is !xed at its bottom, and a uniform pressure of q = 250 is applied to the center 20 × 20 area of its top surface. Using symmetry conditions, only one quarter is discretized by a uniform 5 ×5×5 mesh, as shown in Fig. 7. The material properties are listed as follows: E = 2:1 × 105 , " = 0:3, #y = 250:0 and Ep = 210:0. Fig. 8 presents the load–deVection curves obtained from the proposed element, Weissman’s element [23] and displacement-based element, respectively. In this case, the shape of the elements is regular, and then the results of the present element and those of Weissman’s element are completely identical. Compared with the results of the present element and Weissman’s element, the load–deVection curve for the displacement-based element exhibits a higher slope. The failure of the displacement-based element to achieve the correct slope is in agreement with the !ndings of Nagtegaal et al. [30] for the double-edged-notched tensile specimen, where this behaviour was observed for very !ne meshes.
1168
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
q=250 symmetry A
symmetry
x
y z
Fig. 7. A regular mesh consisting of 5 × 5 × 5 elements for checking element’s performance under con!ned plastic Vow.
8 7
Load Factor
6 5 4 3 2
Present & Weissman's element (Ref. [23]) Displacement-based element
1 0 0.0
0.5
1.0
1.5
2.0
2.5
Displacement at Center Point
Fig. 8. Load–displacement curves for various elements.
3.4. Clamped thin circular plate under uniform pressure In our previous research [24], a circular plate under a central concentrated load was analyzed elastically, and it was found that the accuracy of the plate’s stress and deVection was improved eTciently by using the proposed element. In this example, the same circular plate is further considered as an elasto-plastic problem. The dimensions and material properties are identical to those stated in Ref. [24] (see Fig. 9). For plastic analysis, the uniaxial yield stress #y is equal to 100.0 and the tangential modulus in plastic state Ep is set to be 0.0. The uniform pressure is 0.12 in this case. Unlike the mesh used in the elastic analysis, the two-layer discretization along the thickness direction is adopted for plastic analysis. Moreover, 3750 8-node incompatible 3-D brick elements in ABAQUS [31] commercial software are
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1169
Fig. 9. Mesh for a clamped circular plate under uniform pressure—elastic analysis and elasto-plastic analysis. Table 3 Central deVection of the plate Number of element
24 54 150 384 768
Central deVection of plate Present element
Weissman’s element (Ref. [23])
ABAQUS
1.111 1.260 1.621 1.860 2.280
1.001 1.230 1.580 1.831 2.270
2.441
employed to get the converged result in the plastic analysis for comparison. The !nal results are shown in Table 3. Inspection of this table reveals that the performance of the proposed element is also better than that of Weissman’s element [23] for elasto-plastic case. 4. Conclusions Based on the Hu–Washizu variational principle, a penalty equilibrating 8-node brick element is proposed in this paper. A technique in this element, by which the equilibrium equations in terms of stress components are enforced to be satis!ed through a penalty term in the functional, can reduce the false stresses caused by the elemental irregular shapes and improve the accuracy of analysis. The penalty factor has been carefully designed. Various numerical examples, such as elastic and elasto-plastic analyses, have shown that, compared with currently existing well-performed elements, the proposed element can predict the improved results when the elements are geometrically distorted.
1170
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
Acknowledgements The !nancial supports for this research from Framatome ANP. in France and Integrated Design Inc. in Japan are greatly appreciated. The project from Framatome is “Joint Research Project Agreement-No. 5” with the title: “Numerical investigation on the eOect of residual stresses on the mechanical behavior of 3-D structures” and the project from Integrated Design is: “Three-dimensional elasto-plastic analysis in civil engineering structures”. References [1] K.T. Kavanagh, S.W. Key, A note on selective and reduced integration techniques in !nite element method, Int. J. Numer. Methods Eng. 4 (1972) 148–150. [2] D.S. Malkus, T.J.R. Hughes, Mixed !nite element methods—reduced and selective integration techniques: a uni!cation of concepts, Comput. Methods Appl. Mech. Eng. 15 (1978) 63–81. [3] T.J.R. Hughes, Generalizing of selective integration procedures to anisotropic and nonlinear media, Int. J. Numer. Methods Eng. 15 (1980) 1413–1418. [4] J.C. Simo, R.L. Taylor, K.S. Pister, Variational and projection methods for the volume constraint in !nite deformation elasto-plasticity, Comput. Methods Appl. Mech. Eng. 51 (1985) 177–208. [5] D.P. Flanagan, T. Belytschko, A uniform strain hexahedron and quadrilaterals with high coarse-mesh accuracy, Int. J. Numer. Methods Eng. 17 (1981) 679–706. [6] T. Belytschko, W.E. Bachrach, ETcient implementation of quadrilaterals with high coarse-mesh accuracy, Comput. Methods Appl. Mech. Eng. 54 (1986) 279–301. [7] E.L. Wilson, R.L. Taylor, W.P. Doherty, J. Ghaboussi, Incompatible displacement modes, in: S.J. Fenves, N. Perrone, A.R. Rpbinson, W.C. Schnobrich (Eds.), Numerical and Computer Models in Structural Mechanics, Academic Press, New York, 1973. [8] R.L. Taylor, P.J. Beresford, E.L. Wilson, A nonconforming element for stress analysis, Int. J. Numer. Methods Eng. 10 (1976) 1211–1219. [9] C.C. Wu, M.G. Huang, T.H.H. Pian, Consistency condition and convergence criteria of incompatible elements: general formulation of incompatible functions and its application, Comput. Structures 27 (1987) 639–644. [10] C.C. Wu, Y.Q. Huang, E. Ramm, A further study on incompatible models: revise-stiOness approach and completeness of trial functions, Comput. Methods Appl. Mech. Eng. 190 (2001) 5923–5934. [11] T.H.H. Pian, K. Sumihara, Rational approach for assumed stress !nite elements, Int. J. Numer. Methods Eng. 20 (1984) 1685–1695. [12] E. Reissner, On a variational theorem in elasticity, J. Math. Phys. 29 (1950) 90–95. [13] T.H.H. Pian, P. Tong, Relations between incompatible displacement model and hybrid stress model, Int. J. Numer. Methods Eng. 22 (1986) 173–181. [14] T.H.H. Pian, C.C. Wu, A rational approach for choosing stress terms for hybrid !nite element formulations, Int. J. Numer. Methods Eng. 26 (1988) 2331–2343. [15] C.C. Wu, Y.K. Cheung, On optimization approaches of hybrid elements, Finite Elements Anal. Des. 21 (1995) 111–128. [16] R. Rubinstein, E.F. Punch, S.N. Atluri, An analysis of, and remedies for, kinematics modes in hybrid-stress elements, Comput. Methods Appl. Mech. Eng. 38 (1983) 63–92. [17] E.F. Punch, S.N. Atluri, Development and testing of stable, invariant isoparametric curvilinear 2- and 3-D hybrid stress elements, Comput. Methods Appl. Mech. Eng. 47 (1984) 331–356. [18] J.C. Simo, M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Methods Eng. 29 (1990) 1595–1638. [19] H.C. Hu, Variational Principle in Elasticity and its Application, Science Press, Beijing, 1981 (in Chinese). [20] S.L. Weissman, R.L. Taylor, Treatment of internal constraints by mixed !nite element methods: uni!cation of concepts, Int. J. Numer. Methods Eng. 33 (1992) 131–141.
Y.P. Cao et al. / Finite Elements in Analysis and Design 39 (2003) 1155 – 1171
1171
[21] S.L. Weissman, R.L. Taylor, A uni!ed approach to mixed !nite element methods: application to in-plane problems, Comput. Methods Appl. Mech. Eng. 98 (1992) 127–151. [22] S.L. Weissman, A mixed formulation of nonlinear-elastic problems, Comput. Structures 44 (1992) 813–822. [23] S.L. Weissman, High-accuracy low-order three-dimensional brick elements, Int. J. Numer. Methods Eng. 39 (1996) 2337–2361. [24] Y.P. Cao, N. Hu, J. Lu, H. Fukunaga, Z.H. Yao, A 3-D brick element based on Hu–Washizu variational principle for mesh distortion, Int. J. Numer. Methods Eng. 53 (2002) 2529–2548. [25] R. Piltner, R.L. Taylor, A systematic construction of B-Bar functions for linear and Non-linear mixed-enhanced !nite elements for plane elasticity problems, Int. J. Numer. Methods Eng. 44 (1999) 615–639. [26] B. Fraeijs De Veubeke, Displacement and equilibrium models in the !nite element methods, in: O.C. Zienkiewicz, G.S. Holister (Eds.), Stress Analysis, Wiley, London, 1965. [27] R.H. MacNeal, R.L. Harder, A proposed standard set of problems to test !nite element accuracy, Finite Elements Anal. Des. 1 (1985) 3–20. [28] J.C. Simo, D.D. Fox, M.S. Rifai, On a stress resultant geometrically exact shell model, Part II: The linear theory; computational aspects, Comput. Methods Appl. Mech. Eng. 73 (1989) 53–92. [29] U. Andel!nger, E. Ramm, D. Roehl, 2D- and 3D-enhanced assumed strain elements and their application in plasticity, Proceedings of COMPLAST III, Barcelona, Spain, 1992. [30] J.C. Nagtegaal, D.M. Parks, J.C. Rice, On numerically accurate !nite element solutions in the fully plastic range, Comput. Methods Appl. Mech. Eng. 4 (1974) 153–177. [31] HKS ABAQUS/Standard, Version 5.8, User’s Manual, Hibbit, Karlsson and Sorensen Inc., USA, 1998.