- Email: [email protected]
Relations between enhanced strain methods and the HR method
Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677 www.elsevier.com/locate/cma
Relations between enhanced strain methods and the HR method Francesco Marotti de Sciarra * Dipartimento di Scienza delle Costruzioni, Facolt a di Ingegneria, Universit a degli Studi di Napoli Federico II, Piazzale Tecchio 85, I80125 Napoli, Italy Received 21 May 2001; received in revised form 14 December 2001
Abstract The three-field Hu–Washizu functional is the starting point to develop a comparative analysis of approximate linear elastostatic problems. In particular we will analyse the limitation phenomena concerning two mixed methods, respectively denoted the strain gap method and the mixed enhanced strain method, and the two-field Hellinger–Reissner (HR) finite element formulation. Accordingly two original limitation principles are contributed. We provide the theoretical motivation of the numerical observations, provided in computational literature, on the coincidence of the displacement and stress solutions between the enhanced assumed strain method and the HR method if different choices of shape functions are made. It is further shown that the HR method provides the exact stress solution and the exact nodal displacements with one element. The related analytical motivation are thus given. The theoretical analyses presented in the paper are supported by numerical examples on tipical benchmarks. Ó 2002 Elsevier Science B.V. All rights reserved.
Keywords: Hu–Washizu method; Hellinger–Reissner method; Enhanced strain methods; Finite element formulations; Limitation principles
1. Introduction The standard displacement finite element method exhibits severe inaccuracies in several instances of engineering interest. Examples of these problems consist in locking phenomena for nearly incompressible materials and in a poor accuracy when the bending behaviour is dominant and low order elements are adopted. In order to improve the performance of low order finite elements and to overcome locking phenomena, enhanced strain methods have been proposed in the computational literature and have gained considerable
*
Tel.: +39-81-768-3336; fax: +39-81-768-3332. E-mail address: [email protected] (F. Marotti de Sciarra).
0045-7825/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 2 ) 0 0 1 8 9 - 5
2662
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
success in the scientific community as evidenced by the great number of papers on the subject in linear and non-linear elasticity and in elastoplasticity [2,3,7,12,21,22]. In this paper we consider the enhanced assumed strain (EAS) method proposed by Simo and Rifai [22] and the mixed enhanced strain (MES) method proposed by Kasper and Taylor [12]. Enhanced strain methods can be formulated starting from the Hu–Washizu functional in which strains, stresses and displacements are assumed as independent fields. Stress and strain fields are required to be square integrable functions and hence can be interpolated independently in each element of the finite element mesh. The corresponding displacement field must, on the contrary, fulfill the interelement continuity and the external constraint conditions. In the original approach of the EAS method [22] the stress fields are assumed to be orthogonal to the enhanced strains and therefore stresses are dropped out from the variational formulation. As a consequence several a posteriori stress recovery strategies, including the discussed question of a stress recovery based on the stress obtained from the constitutive elastic relation, have been proposed by several authors [2,5,14,22]. A detailed analysis on these issues has been carried out in [21] where a new method, called the strain gap method (SGM), has been contributed. In the MES method the strain field is additively decomposed in two sets. The former contains the stress shape functions and the latter contains the enhanced strain shape functions which are assumed to be orthogonal to the stress shape functions. In [18,21] the basic question of well-posedness of discrete mixed enhanced problems has been addressed, a new enhanced method has been proposed and a critical analysis of both EAS and MES methods [12,22] has been carried out. In the present paper we provide a theoretical and numerical analysis of the relations existing between the finite element solutions obtained from the two-field Hellinger–Reissner (HR) method and from the MES method and the SGM. We show that the stress and the enhanced strain shape functions, usually adopted in the MES method [12], fulfill a limitation principle so that the method does in fact collapse into a two-field HR method. The theoretical treatment is supported by the numerical analysis of a plate subjected to simple bending which is usually assumed as a valuable benchmark for enhanced methods, see e.g. [14,16]. The computational analysis here performed reveales that the two-field HR method provides the exact stress field and the exact nodal displacements with just one element mesh. A detailed explanation of this peculiar result is provided. We contribute a sufficient condition which ensures that the stress solution obtained from the two-field HR method is in fact the exact solution. The proof of the coincidence of the nodal displacements with the exact displacements is based on an application of the virtual work principle. In the second part of the paper, the SGM method is considered. We provide a theoretical foundation to the numerical observation made in [2] concerning the coincidence of the displacement and stress solutions obtained from the EAS method and the two-field HR method under usual choices of stress and enhanced strain shape functions. Finally numerical examples of the limitation phenomenon concerning the SGM, the MES method and the HR method are provided.
2. Mixed variational formulations In this paper we will perform the analysis of mixed methods for a continuous linearly elastic structural model defined on a bounded domain X of an Euclidean space with boundary oX and closure X ¼ X [ oX. The kinematic space V is endowed with a suitable Hilbert topology and its dual F is the space of external forces. To each displacement field v 2 V there corresponds a boundary field Cv 2 oV. The dual
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2663
space oF of oV is the space of boundary tractions t and body forces b belong to the Hilbert space H of square integrable vector fields on the structural model. External forces are then collected in the set ‘ ¼ fb; tg 2 F. Strains e and stresses r belong to the Hilbert space H of square integrable vector fields on X [20], the strain and stress spaces will be also denoted by D S H. The structure is subject to homogeneous boundary conditions so that admissible displacements belong to the subspace L of conforming displacements. The kinematic operator B 2 LinfL; Dg associates the strain fields e with the corresponding displacement fields v and it is a continuous linear map from L into D with closed range [6] and with a finite dimensional kernel [19]. The stress fields belong to the Hilbert space S ¼ fr 2 H: B0o r 2 H g where the differential operator 0 Bo 2 LinfS; H g is the formal adjoint of B and provides the body forces corresponding to the stress field r 2 S. For simplicity, we will not consider in the sequel imposed strains and displacements since they can be carried in without any significant modification of the reasonings presented in this paper. The linear elasticity operator E 2 LinfH; Hg is continuous, symmetric and H-elliptic, that is ððEe; eÞÞH P cE kek2H
8 e 2 D:
The elastic strain energy of the structural model is provided by the convex quadratic functional / : D 7! R [ fþ1g given by /ðeÞ ¼ 12ððEe; eÞÞH . Its conjugate convex functional [10] / : S 7! R [ fþ1g represents the complementary elastic energy and is given by / ðrÞ ¼ sup fððr; eÞÞH /ðeÞg ¼ 12ððCr; rÞÞH ; e2D
1
where C ¼ E is the elastic compliance. The structural problem can be written in an operator form as [17] Bu ¼ e;
B0 r ¼ ‘;
r ¼ d/ðeÞ ¼ Ee;
ð2:1Þ
where B0 2 LinfS; Fg is the dual operator of B and is defined as ððr; BrÞÞH ¼ ððB0 r; rÞÞH for any u 2 V and r 2 S. The solutions of the elastostatic problem can be characterized as stationarity points of the Hu–Washizu functional [24]: W ðe; r; uÞ ¼ /ðeÞ þ ððr; Bu eÞÞH h‘; ui;
ð2:2Þ
with u 2 L, e 2 D, r 2 S,where the linear functional h‘; ui ¼ ðb; uÞH þ ht; Cui yields the virtual work of body forces and boundary tractions. Let us now introduce the strain gap field g 2 D as the difference between the compatible strain Bu 2 D and the strain e 2 D in the form g ¼ Bu e. Accordingly the Hu–Washizu functional can be re-written in the form e ðu; g; rÞ ¼ /ðBu gÞ þ ððr; gÞÞ h‘; ui ¼ 1ððEðBu gÞ; Bu gÞÞ þ ððr; gÞÞ h‘; ui W H H H 2
ð2:3Þ
with u 2 L, g 2 D, r 2 S. The expression of the HR functional Rðr; uÞ can be obtained starting from the HW functional (2.2) by enforcing the constitutive relation in terms of the Fenchel’s equality [10]: /ðeÞ þ / ðrÞ ¼ ððr; eÞÞH
2664
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
to get Rðr; uÞ ¼ / ðrÞ þ ððr; BuÞÞH h‘; ui ¼ 12ððr; CrÞÞH þ ððr; BuÞÞH h‘; ui;
ð2:4Þ
with u 2 L and r 2 S.
3. Discrete formulations The unknown fields fe; g; r; ug 2 D D S V can be interpolated over a finite element Xe by means of the matrices Nee ðxÞ, Neg ðxÞ, Ner ðxÞ and Neu ðxÞ, collecting the strain, strain gap, stress and displacement interpolating shape functions, according to the expressions [11]: e
Deh ¼ feeh : eeh ðxÞ ¼ Nee ðxÞpee and pee 2 Rne ; x 2 Xe g; e
Geh ¼ fgeh : geh ðxÞ ¼ Neg ðxÞpeg and peg 2 Rng ; x 2 Xe g; e
Seh ¼ freh : reh ðxÞ ¼ Ner ðxÞper and per 2 Rnr ; x 2 Xe g; e
Veh ¼ fueh : ueh ðxÞ ¼ Neu ðxÞpeu and peu 2 Rnu ; x 2 Xe g: The domain decomposition induced by the meshing of X depends on a parameter h which goes to zero as the finite mesh is refined. The apex e is referred to quantities pertaining to the element Xe . The dimension of the interpolating subspaces Deh , Geh , Seh and Veh has been denoted respectively by nee , neg , ner and neu . The matrices Nee , Neg , Ner and Neu are assumed to have linearly independent columns and the vectors pee , peg , e pr and peu collect the strain, strain gap, stress and displacement parameters pertaining to the element Xe . Let N be the total number of elements of the finite element discretization, then the interpolating spaces are grouped in the following product spaces: Dh ¼
N Y
Deh ;
e¼1
Gh ¼
N Y
Geh ;
Sh ¼
e¼1
N Y
Seh ;
Vh ¼
N Y
e¼1
Veh :
e¼1
Since strain, strain gap and stress fields can be discontinuous at the element interfaces, the corresponding global interpolated fields are simply the collection of the corresponding field defined over each element Xe : eh ¼ fe1h ; e2h ; . . . ; eN h g 2 Dh ;
gh ¼ fg1h ; g2h ; . . . ; gN h g 2 Gh ;
rh ¼ fr1h ; r2h ; . . . ; rN h g 2 Sh :
As a consequence the virtual work performed by an interpolating stress by an interpolating strain or strain gap field is defined as the sum of the elementwise virtual work, i.e. N N Z N Z X X X ððrh ; eh ÞÞ ¼ hhreh ; eeh iiXe ¼ reh eeh ¼ Ner per Nee pee : e¼1
e¼1
Xe
e¼1
Xe
We consider a conforming finite element approximation so that at the element interface the displacement field is continuous. The conforming interpolated displacement fields uh ¼ fu1h ; u2h ; . . . ; uN h g 2 L h Vh satisfy the homogeneous boundary constraints and the interelement continuity conditions. Here Lh is a proper subspace of the product space Vh . The strain, strain gap and stress global parameters qe , qg and qr are the collection of the element parameters pee , peg and per . On the contrary the displacement element parameters peu are expressed in terms of the global parameters qu by means of the finite element assembly operator Aeu .
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2665
3.1. Discrete variational formulations The interpolated counterpart of the Hu–Washizu functional (2.2) can be obtained by introducing the fields feh ; rh ; uh g 2 Dh Sh Lh in the expression of W to get W ðeh ; rh ; uh Þ ¼ 12ððEeh ; eh ÞÞ þ ððrh ; Buh eh ÞÞ h‘; uh i:
ð3:1Þ
A solution feho ; rho ; uho g 2 Dh Sh Lh of the discretized structural model can be obtained by enforcing the stationarity of the Hu–Washizu functional (3.1) to get 8 8 deh 2 Dh ; < ððEeho ; deh ÞÞ ððrho ; deh ÞÞ ¼ 0 ððdrh ; eho ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0 8 drh 2 Sh ; ð3:2Þ : ððrho ; Bduh ÞÞ ¼ h‘; duh i 8 duh 2 Lh : The matrix form of the discrete structural problem can be obtained starting from the stationarity conditions (3.2) and is given by 3 2 3 2 e e e e qe qe H Q O qe JeT JeT O N e H Je e Q Jr X e e 5 4 JeT QeT Je M qr ¼ q ¼ 4 QT O S 5 qr O JeT r e r S Au r eT eT e q e¼1 O ST O qu qu O Au S J r O u o o N X o ¼ o : ð3:3Þ ¼ eT e fu e¼1 A f u u Here Jee , Jer are the canonical extractors which pick up from the global lists qe and qr the local parameters e pee and per pertaining to the element Xe . The operator AeT u transforms the local forces f u into the corresponding global ones f u and the following local matrices have been defined: Z Z Z e e e eT e eT e e Ne ðxÞE ðxÞNe ðxÞ; Q ¼ Ne ðxÞNr ðxÞ; S ¼ NeT H ¼ r ðxÞB Nu ðxÞ Xe
f eu ¼
Z Xe
NeT u ðxÞbðxÞ þ
Xe
Z
Xe
ðCNeu ÞT ðxÞtðxÞ;
oXe
where B and E denotes the matrix form of the operators B and E. The SGM [21] is obtained starting from the interpolated counterpart of the Hu–Washizu functional e h ðuh ; gh ; rh Þ. In fact by adding up the contributions of each non-assembled element and imposing that the W interpolating displacement uh satisfies the conformity requirement we get e h ðuh ; gh ; rh Þ ¼ 1ððEðBuh gh Þ; Buh gh ÞÞ þ ððrh ; gh ÞÞ h‘; uh i W 2
ð3:4Þ
where fuh ; gh ; rh g 2 Lh Gh Sh . The solutions are obtained by enforcing the stationarity of (3.4) to get 8 < ððEBuho ; Bduh ÞÞH ððEgho ; Bduh ÞÞH ¼ h‘; duh i 8 duh 2 Lh ; ð3:5Þ ððEgho ; dgh ÞÞH þ ððrho ; dgh ÞÞH ððEBuho ; dgh ÞÞH ¼ 0 8 dgh 2 Gh ; : ððdrh ; gho ÞÞH ¼ 0 8 drh 2 Sh : e h ¼ Gh \ S? . The strain gaps belonging The third equation of relation (3.5) amounts to require that gh 2 G h ? e to the subspace G h ¼ Gh \ Sh are referred to as effective strain gaps since they effectively contribute to relax the compatibility condition.
2666
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
The matrix form of the discrete problem is given by 2 3 2 qu X qu AeT Ke Aeu AeT GeT Jeg O N K u u 6 eT e eT e e eT e e 7 e f Jg L Jg Jg P Jr 5 qg ¼ 4 G M qg ¼ 4 Jg G Au q e¼1 q eT e O O JeT O r r r P Jg eT e fu N Au f u X ¼ o ¼ o : o e¼1 o
3 O qu P 5 qg O qr
GT L PT
ð3:6Þ
The component submatrices appearing in (3.6) are defined by Z Z e e e Le ¼ NeT ðxÞE ðxÞN ðxÞ; P ¼ NeT g g g ðxÞNr ðxÞ; Xe Xe Z Z T e e Ge ¼ NeT ðxÞE ðxÞB N ðxÞ; K ¼ ðB Neu Þ ðxÞE ðxÞB Neu ðxÞ: g u Xe
ð3:7Þ
Xe
Remark 1. In [21] the SGM has been proposed on the basis of a suitable variant of the HW functional. It is then performed the analysis of the well-posedness conditions which ensure uniqueness of the discrete solution in terms of displacements, strain gaps and stresses. It is shown that the SGM can be split in the reduced problem, formulated in terms of displacements and strain gaps, and in a stress recovery problem depending on the solution of the reduced problem. The comparison between the EAS method and the SGM shows that the EAS method is partially ill-posed since an orthogonality condition drops out the stress parameters. As a consequence the EAS method coincides only with the reduced problem of the SGM and the stress recovery remains, for the EAS method, an open problem. The expression of the HR functional Rðr; uÞ in terms of interpolating fields is Rh ðrh ; uh Þ ¼ / ðrh Þ þ ððrh ; Buh ÞÞ h‘; uh i
rh 2 Sh ; uh 2 Lh :
ð3:8Þ
An approximate solution frho ; uho g 2 Sh Lh is obtained by enforcing the stationarity of Rh to get: ððdrh ; Crho ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0 8 drh 2 Sh ; ð3:9Þ ððrho ; Bduh ÞÞ h‘; duh i ¼ 0 8 duh 2 Lh : The matrix form of the discrete problem is given by X N qr K JeT ReT Jer JeT Te Aeu qr r r A ¼ ¼ eT eT e qu q TT A T J O u u r e¼1 where Re ¼
Z Xe
1 e NeT r ðxÞE ðxÞNr ðxÞ;
Te ¼
Z
e NeT r ðxÞB Nu ðxÞ:
T O
qr qu
¼
N X AeT f e u
e¼1
o
u
fu ; ¼ o
ð3:10Þ
ð3:11Þ
Xe
From the computational point of view it is more convenient to carry out the assembly operation after the condensation, at the element level, of the strain and stress parameters to put the global discrete problem in terms of the sole displacement parameters qu . The procedure has been illustrated in [18,21]. 4. Mixed enhanced strain method The mixed enhanced strain (MES) method [12] is based on the following decomposition of the interpolating strain space:
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
eh D h ¼ Sh D
2667
ð4:1Þ
e h is the subspace of enhanced strains ~eh given by where D eh ¼ D
N Y
e e; D h
e e ¼ f~ee : ~ee ðxÞ ¼ Ne ðxÞpe and pe 2 Rn~ee ; x 2 Xe g: D ~e ~e ~e h h h
e¼1
e h ¼ fog so that stress and enhanced strain The symbol denotes the direct sum and requires that Sh \ D shape functions must be linearly independent. e h , the interpolated counterpart of the Hu–Washizu funcSetting eh ¼ eh þ ~eh with eh 2 Sh and ~eh 2 D tional W in terms of (eh ; ~eh ; rh ; uh ) is given by Mh ðeh ; ~eh ; rh ; uh Þ ¼ 12ððEeh ; eh ÞÞ þ 12ððE~eh ; ~eh ÞÞ þ ððEeh ; ~eh ÞÞ þ ððrh ; Buh eh ~eh ÞÞ h‘; uh i; The MES method follows from the stationarity conditions of Mh and is given by 8 ððEeho ; deh ÞÞ þ ððE~eho ; deh ÞÞ ððrho ; deh ÞÞ ¼ 0 8 deh 2 Sh ; > > < e h; ððEeho ; d~eh ÞÞ þ ððE~eho ; d~eh ÞÞ ððrho ; d~eh ÞÞ ¼ 0 8 d~eh 2 D > ððdrh ; eho þ ~eho ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0 8 drh 2 Sh ; > : ððrho ; Bduh ÞÞ ¼ h‘; duh i 8 duh 2 Lh :
uh 2 L h :
ð4:2Þ
4.1. The MES method vs the HR principle A sufficient condition for the validity of a limitation phenomenon between the MES method (4.2) and the HR method (3.9) is reported in the next statement. Proposition 4.1. (A limitation principle). The solutions of the approximate problems obtained from the MES and the HR methods coincide in terms of rh and uh if the same interpolating spaces Sh and Lh are assumed and if e h: CSh Sh D
ð4:3Þ
e h Sh Lh be a stationarity point for the MES funcProof. Let the quartet feho ; ~eho ; rho ; uho g 2 Sh D tional. The inclusions (4.3) ensure that eh þ ~eh Crh 2 Dh
e h ; eh ; rh 2 Sh : 8~eh 2 D
Summing up the first and second equation of (4.2) which represents the constitutive relations, and by virtue of the positive definiteness of the elastic stiffness E, we have: Eðeho þ ~eho Þ rho 2 D? h \ EDh ¼ fog; e h . The orthogonality relation ? is since Dh is decomposed in the direct sum of the subspaces Sh and D intended according to the inner product in L2 ðXÞ. Then the interpolating strain eho ¼ eho þ ~eho and the interpolating stress rho fulfil the constitutive relation: eho þ ~eho ¼ Crho which inserted in the second relation of (4.2) yields ððdrh ; Crho ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0 8 drh 2 Sh :
2668
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
The MES method (4.2) can then be expressed by the following reduced set of equations:
ððdrh ; Crho ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0 8 drh 2 Sh ; ððrho ; Bduh ÞÞ ¼ h‘; duh i 8 duh 2 Lh ;
which coincide with the stationarity conditions (3.9) of the HR functional.
It is worth noting that the condition (4.3) is local since interpolating strains and stresses are defined elementwise so that we can restrict ourselves to analyse a single element. In a finite element analysis, it is compelling to verify the condition (4.3) in terms of quantities pertaining to the reference element K. To this end we denote by ve : K 7! Xe the isoparametric map from the reference element K to the generic element Xe of the finite element mesh and we set v ¼ fve ; e ¼ 1; . . . ; Ng. Since e h ¼ vð D e K Þ, to check the inclusion (4.3) it is sufficient to verify Sh ¼ vðSK Þ and D e K: CSK SK D
ð4:4Þ
Let us now show that the strain and stress shape functions adopted for the MES method, see [12,18], meet the above limitation principle for distorted and undistorted elements. In the MES method, the stress subspace for 4-node quadrilateral elements with standard bilinear shape functions (Q1), where K ¼ , are given by [15] 2 3 1 0 0 g 0 S ¼ span4 0 1 0 0 n 5 ð4:5Þ 0 0 1 0 0 and the strain subspace is decomposed, according to (4.1), in the form 2 3 n 0 e where D e ¼ span4 0 g 5: D ¼ S þ D 0 0 For plane stress problems, 2 1 14 m CS ¼ span E 0
the subspace CS becomes 3 m 0 g mn 1 0 mg n 5; 0 2ð1 þ mÞ 0 0
ð4:6Þ
ð4:7Þ
where E and m denote the Young’s modulus and the Poisson’s ratio. A direct inspection shows that the base vectors of the subspace CS , given by the columns of the matrix (4.7), can be obtained as a linear combination of the vectors spanning D , see (4.5) and (4.6). The condition (4.4) is then fulfilled, a limitation phenomenon occurs and no advantage is gained by the use of the MES method in place of the simpler HR method. Additional transformation rules which preserve the point-wise inner product between stress and strain tensors can be envisaged but it seems that they can only be motivated by an a posteriori evaluation of the quality of the numerical results. Examples are provided by the push/pull transformations of differential geometry [13] which provides, for plane problems see e.g. [12,22], the expressions reh ðxÞ ¼ Teo re ½v1 e ðxÞ;
geh ðxÞ ¼
Joe e 1 TeT o g ½ve ðxÞ; e J ½v1 ðxÞ e
ð4:8Þ
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2669
where Teo is the value, at the origin of the reference element, of the matrix field 2 2 3e F11 F122 2F11 F12 5 ðnÞ with n 2 K: Te ðnÞ ¼ 4 F212 F222 2F21 F22 F11 F21 F12 F22 F11 F22 þ F12 F21 Accordingly the inner product in the real space is given by Z Z Z Joe eT e 1 e T reh ðxÞ geh ðxÞ dx ¼ Teo re ½v1 ðxÞ g ½v ðxÞ dx ¼ J rðnÞ gðnÞ dn: o e e o J e ½v1 Xe Xe K e ðxÞ Adopting the interpolations (4.8), the condition (4.3) is not fulfilled in the case of distorted elements so that the MES method and the HR method provide different results. A numerical comparative analysis is provided in Section 6.
5. Limitation phenomena between SGM and HR Let us now show sufficient conditions to be fulfilled in order that the SGM method collapses into the HR method. Proposition 5.1. (Limitation principle between the SGM and HR method). The approximate solution fuh ; rh g provided by the SGM (3.5) and by the two-field HR method (3.9) coincide if the same subspaces Lh and Sh are adopted and if the following condition is met: e h: EBLh Sh G
ð5:1Þ
e h , the condition (5.1) implies that EðBuh gh Þ rh 2 Gh . By Proof. By virtue of the choice Gh ¼ Sh G virtue of the second condition of (3.5) we have EðBuh gh Þ rh 2 Gh [ G? h ¼ fog: Accordingly we have rh ¼ EðBuh gh Þ. The constraint condition provided by the third equation of (3.5) e h , so that we have ensures that the strain gaps are effective, i.e. gh 2 G e h S? gh ¼ Buh Crh 2 G h
()
ððdrh ; Buh ÞÞ ððdrh ; Crh ÞÞ ¼ 0
8 drh 2 Sh
ð5:2Þ
and the first equation of (3.5) yields ððrh ; Bduh ÞÞ ¼ h‘; duh i 8 duh 2 Lh which coincide with the stationarity conditions of the HR method (3.9). The statement is thus proved.
ð5:3Þ
An analogous limitation principle between the EAS method and the HR method has been proved in [5] by adding the further hypothesis that the approximate stress pertaining to the EAS method is obtained by means of the constitutive elastic relation. The above Proposition 5.1 shows that such an assumption is inessential and can thus be avoided. In a finite element analysis, the condition (5.1) deserves a carefull discussion since it involves the conforming displacement subspace Lh and the kinematic operator B. To this end we note that this condition is global due to the presence of Lh which depends on the a priori unknown element assembly operations. This shortcoming can be circumvented by considering the larger non-conforming displacement space Vh Lh so that the global condition (5.1) can be substituted by the e h. sufficient local condition EBVh Sh G
2670
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
The kinematic operator acting on the displacement fields in the reference element K will be denoted by BeK and is defined by BeK uK ðnÞ ¼ Bueh ðxÞ for any x 2 Xe where ueh is the restriction of uh to the element Xe e and n ¼ v1 e ðxÞ is the isoparametric coordinate in the reference element. We then set BK ¼ fBK ; e ¼ 1; . . . ; Ng. Denoting by gK ¼ v1 ðuh Þ the collection of the displacement fields obtained by shifting back to K the e h can then be fields ueh , we have Buh ¼ Bv1 ðuh Þ ¼ BK gK . The sufficient condition above EBVh Sh G written in the reference element K in the form eK: EBK VK SK G
ð5:4Þ
It is apparent that the kinematic operator BK on K depends on the map v and hence the condition (5.4) cannot be checked in terms of quantities pertaining only to the reference element K. These difficulties have not been fully realized in the literature since in previous treatments, see e.g. [2,4,22,25], conditions of this kind has been simply imposed in the reference element K. In fact for a general isoparametric map, the condition (5.4) can be only checked in the form eK EBVK SK G
ð5:5Þ
and confide that the element distorsion does not make it fail on the real elements. On the basis of the previous analysis, let us now provide a theoretical explanation of the numerical observation made in [2] that the EAS method [22] and the HR method yield the same solution in terms of displacements and stresses. e h and Sh with reference to We analyse the different choices proposed in [2] concerning the subspaces G quadrilateral Q1 elements. We consider the Pian and Sumihara stress interpolation and the seven-parameter strain gap interpolation: 2 3 2 3 1 0 0 g 0 n 0 0 0 ng 0 0 e ¼ span4 0 g 0 0 0 ng 0 5: S ¼ span4 0 1 0 0 n 5; G ð5:6Þ 0 0 1 0 0 0 0 n g 0 0 ng Stress and strain gap shape functions in the isoparametric space belong to the four-dimensional space Q1 e is direct so that the subspace S G e generated by the polynomials 1, n, g, and ng [8]. The sum S þ G turns out to be twelve-dimensional. Since the compatible strain subspace BV is given by 2 3 1 0 0 g 0 BV ¼ span4 0 1 0 0 n 5; 0 0 1 n g the condition (5.5) is fulfilled. As a further example we consider the eight-parameter stress interpolation together with a four-parameter enhanced strain (EAS4) interpolation which coincides with the modified incompatible mode approximation of Taylor et al. [23]: 2 3 2 3 1 0 0 g 0 ng 0 0 n 0 0 0 e ¼ span4 0 g 0 0 5: S ¼ span4 0 1 0 0 n 0 ng 0 5; G ð5:7Þ 0 0 1 0 0 0 0 ng 0 0 n g Finally the following twelve-parameter stress interpolation and no enhanced strains were considered in [2]:
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2
1 S ¼ span4 0 0
0 1 0
0 0 1
g 0 0
0 ng n 0 0 0
0 ng 0
0 0 ng
n 0 0 g 0 0
0 0 n
3 0 0 5; g
e ¼ f0g: G
2671
ð5:8Þ
It can be proved, following the arguments shown above, that the choices (5.7) and (5.8) fulfill the condition reported in the Proposition 5.1 so that the SGM collapse into the HR method. It is worth noting that the choice (5.8) of the shape functions fulfils the inclusion EBVK SK which ensures that the displacement solution obtained from the HR method and the standard displacement method coincides [1,9]. Let us now analyse the stress and strain gap shape functions considered in [22]: 2 3 2 3 1 0 0 g 0 n 0 0 0 ng e ¼ span4 0 g 0 0 ð5:9Þ S ¼ span4 0 1 0 0 n 5; G ng 5: 0 0 1 0 0 0 0 n g n2 g2 Following the arguments previously reported it is immediate to prove that the condition (5.5) is fulfilled and the limitation phenomenon occurs.
6. Numerical examples The limitation principle reported in the Proposition 4.1 is tested by means of a simple structural example consisting in a plate of unit thickness subjected to simple bending as shown in Fig. 1. This example is usually adopted as a benchmark for enhanced methods [4,14]. The analytical solution in terms of horizontal and vertical displacements is given by 2p uðx; yÞ ¼ Eb x b2 y ; p vðx; yÞ ¼ Eb ½x2 þ mðy 2 byÞ: The only non-vanishing stress component is rx whose expression is given by 2p b rx ðx; yÞ ¼ y : b 2 Setting a ¼ 10, b ¼ 2, E ¼ 1500, m ¼ 0:25, p ¼ 3000, the analytical solution of the bending problem in terms of displacements and stresses is uðx; yÞ ¼ 2xð1 yÞ, vðx; yÞ ¼ x2 þ 14ðy 2 2y) and rx ðx; yÞ ¼ 3000ð1 yÞ.
Fig. 1. The bending plate problem discretized by means of Q1 elements.
2672
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
A first example is carried out for two meshes of 4-node quadrilateral elements with bi-linear displacement interpolations (Q1): 1 1 rectangular element and 16 8 uniform rectangular elements. The Pian and Sumihara stress shape functions (4.5) and the strain shape functions (4.6) are adopted. The horizontal and vertical displacements and the horizontal stress at the corner node A of the plate are evaluated adopting the MES method, the HR method and the standard displacement formulation and they are compared with the solution of the continuous problem. The results are reported in Table 1. A further example is devoted to analyse how the distortion of the elements affects the discrete solution. A classical benchmark is the plate reported in Fig. 1 when it is modelled with two distorted quadrilateral elements. Figs. 2 and 3 respectively show the vertical displacement of the point A and the horizontal stress rx at the point B against the distortion parameter d. The shape functions (4.5) and (4.6) are adopted for the MES and the HR methods. Figs. 2 and 3 show that the displacement and the stress solutions obtained from the MES method and the HR formulation are identical as predicted by the limitation principle 4.1. The MES-p and the HR-p plots of Figs. 2 and 3 are obtained by including the push-transformation (4.8) in the definition of stress and strain shape functions.
Table 1 The displacement components u and v and the horizontal stress rx at the node A for the following Q1 meshes (panel a) one element mesh; (panel b) 16 8 uniform rectangular mesh MES
HR
Panel a u v rx
DISPL
EXACT
20 100 3000
20 100 3000
2 9 289
20 100 3000
Panel b u v rx
20 100 3000
20 100 3000
19 96 2815
20 100 3000
Fig. 2. Vertical displacement of the node A for a cantilever rectangular plate meshed with two distorted elements.
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2673
Fig. 3. Horizontal stress at the node B for a cantilever rectangular plate meshed with two distorted elements.
The resulting shape functions on the real element do not meet the limitation principle and, in fact, the solutions obtained from these two methods are different. In particular the push-transformation has beneficial effects for the MES method in terms of stresses but it shows a poor performance in terms of displacements for increasing values of d. The displacement and stress solutions obtained from the SGM and the HR method do coincide for undistorted elements if the shape functions (5.9) are used. For distorted elements, the presence of the kinematic operator in the expression (5.1) implies that the interpolating subspaces do not fulfil the inclusion appearing in (5.1) so that the two methods provide different solutions, see Figs. 2 and 3. The examples reveal that, assuming a one element mesh, the interpolating stress rh coincide with the exact stress solution. Moreover, even if the interpolating displacement field uh is quite different from the exact solution, the displacement parameters qu are exact. In the next section we will explain in detail this result.
7. The finite element exact solution Let us now provide a sufficient condition which ensures that the stress solution rho obtained from the HR method coincides with the exact stress solution ro pertaining to the continuous model. In particular we will show that this condition is met by the displacement and stress shape functions adopted for the finite element analysis of the bending problem (Fig. 1). Further we will provide a theoretical motivation to the numerical observation that the nodal displacements qu obtained by the HR method are exact. To this end we prove the next result. Proposition 7.1. (The exact stress solution) Let fro ; uo g 2 S L be the solution of the continuous problem and frho ; uho g 2 Sh Lh be the approximate solution of the HR method.
2674
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
Then if r o 2 Sh
and
Buo 2 BLh S? h
ð7:1Þ
the approximate and the exact stress solutions coincide, i.e. rho ¼ ro . Proof. The exact stress solution ro 2 S fulfils the approximate equilibrium equation (3.9), part 2 ððro ; Bduh ÞÞ h‘; duh i ¼ 0
8 duh 2 Lh
ð7:2Þ
since the inclusion Lh L hold true. It remains to prove that the exact stress solution ro satisfies the approximate elastic compatibility equation (3.9), part 1 which is hereafter rewritten for convenience: ððdrh ; Cro ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0
8 drh 2 Sh :
ð7:3Þ
To this end we note that second condition of (7.1) implies that there exists a unique interpolating conforming displacement uh 2 Lh such that Buo þ Buh 2 S? h
()
ððdrh ; Buo ÞÞ þ ððdrh ; Buh ÞÞ ¼ 0
8 drh 2 Sh :
Further, noting that the equalities Cro ¼ eo ¼ Buo hold true for the exact solution, the orthogonality relation above can be equivalently rewritten in the form ððdrh ; Cro ÞÞ þ ððdrh ; Buh ÞÞ ¼ 0
8 drh 2 Sh :
ð7:4Þ
fro ; uh g
Accordingly, the pair fulfils the approximate equilibrium equation (7.2) and the approximate elastic compatibility (7.4) and hence fro ; uh g is the solution of the HR approximate problem (3.9). By virtue of the uniqueness of the approximate solution we infer that fro ; uh g fro ; uho g so that ro is the stress solution of the approximate problem. The proposition is thus proved. 7.1. The stress solution of the plate bending problem Let us now show that a one element mesh provides the exact stress solution for the linear elastic bending problem presented in the previous section. For a one element Q1 mesh, the set of shape functions spanning the subspace BLh is given by 2 3 1 0 0 g 0 ð7:5Þ BLh ¼ span BNu ¼ span4 0 1 0 0 n 5: 0 0 1 n g The compatible strain Buo 2 D of the continuous problem is given by g 1 Buo ¼ mg E 0
ð7:6Þ
and can be obtained in terms of a linear combination of the columns of the matrix (7.5) as follows: g 0 0 1 1 m Buo ¼ 0 0 g : ð7:7Þ E E E n n 0 The last two vectors appearing in (7.7) are L2 ðÞ––orthogonal to the subspace S generated by the assumed Pian–Sumihara stress interpolation (4.5) so that they belong to S? .
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2675
Noting that the stress solution ro of the continuous problem belong to the interpolating stress subspace S : g ro ¼ 0 2 S ; 0 the requirements of the Proposition 7.1 are fulfilled so that a one element mesh provides the exact stress solution. 7.2. The displacement solution of the plate bending problem Let us now explain why the interpolation fields, assumed for the analysis of the plate bending problem reported in Fig. 1, provide the exact nodal displacements for a one element mesh. Since Crho ¼ Cro ¼ Buo , the elastic compatibility condition (3.9), part 1 can be rewritten as ððdrh ; Buo ÞÞ ¼ ððdrh ; Buho ÞÞ;
8 drh 2 Sh :
ð7:8Þ
The five-dimensional Pian–Sumihara interpolating stress subspaces Sh is spanned by the five columns of the matrix (4.5) which will be denoted by si with i ¼ 1; . . . ; 5. The elastic compatibility relation (7.8) can then be tested in terms of the base vectors fsi g of Sh in the form: ððsi ; Buo ÞÞ ¼ ððsi ; Buho ÞÞ;
i ¼ 1; . . . ; 5:
ð7:9Þ
Since the exact compatible strain solution Buo is given by (7.6), it is apparent that the l.h.s. of (7.9) is vanishing for i ¼ 1, 2, 3, 5. Hence the equality (7.9) yields ððsi ; Buho ÞÞ ¼ 0; i ¼ 1; 2; 3; 5 and
ððsi ; Buho ÞÞ 6¼ 0; i ¼ 4:
ð7:10Þ
The interpolating compatible strain Buho is expressed by a linear combination of the base vectors of BLh , see (7.5), so that the relations (7.10) imply that the deformation modes provided by the columns 1, 2, 3 and 5 of (7.5) cannot be active. In conclusion the hourglass mode in the x-direction provides the only compatible deformation of the approximate model, see Fig. 4. By virtue of the equality (7.9), the virtual work principle allows us to write ððs4 ; Bðuo uho ÞÞÞ ¼ MðDuo Duho Þ ¼ 0
ð7:11Þ
where M denotes the external couple in equilibrium with the stress distribution s4 and Duo , Duho represent the exact and the approximate relative rotations of the beam end sections. The equality (7.11) implies Duo ¼ Duho and, being the left end section fixed, the absolute rotation u of the right end section is exact. It is well known that the curvature of the continuous model is constant so that the displacement of the centroid of the right end section is 12ul (see Fig. 5a). The approximate model provide the same result since
Fig. 4. Hourglass mode in the x-direction.
2676
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
Fig. 5. Deflection of the beam (a) for the continuous model; (b) adopting a one element Q1 analysis.
the rotation to be superimposed to the hourglass mode to restore the external constraints is 12u so that the displacement of the centroid of the right end section is 12ul (see Fig. 5b).
8. Closure The analysis carried out in this paper provide a guide to check if the shape functions to be adopted in the MES method and in the SGM meet a limitation phenomenon so that the same displacement and stress solutions of the HR method are obtained. The SGM has been recently contributed by the authors and includes the EAS method as a special ill-posed case. In particular we have shown that a limitation principle between the MES method and the HR method is fulfilled by the shape functions usually adopted for the interpolating stresses and enhanced strains for undistorted meshes. In this case no advantage is obtained by the adoption of an enhanced method with respect to the two-field HR formulation. Further a limitation principle has been provided between the SGM method and the HR formulation. We have provided a theoretical motivation to the numerical observations made by Andelfinger and Ramm which have shown the coincidence of the displacement and stress solutions between the SGM method and the HR method for several choices of stress and enhanced strain shape functions. The proposed argumentations are based on the properties of the interpolating subspaces and can be easily repeated if different stress and enhanced strain shape functions are considered. Two further results have been pointed out in this paper with reference to the HR method. The former is a sufficient condition, expressed in terms of interpolating subspaces, which ensures that the interpolating stress obtained from the HR method coincides with the exact stress solution. This condition is met by the bending problem of the plate so that the exact stress solution is obtained by means of the HR method with one element. The latter result is the theoretical motivation of the coincidence of the exact displacements with the nodal displacements obtained from the HR method with one element. The theoretical considerations have been supported by the computational analysis of a plate subjected to simple bending, which is usually considered as a benchmark for enhanced methods, and the influence on the solution of the distortion of the elements has also been pointed out. Accordingly, the plate bending problem cannot be considered as an effective benchmark for enhanced methods since the HR method provides the exact solutions with one element. Enhanced methods can then
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2677
share an analogous performance if a limitation phenomenon occur otherwise they behave worse than the two-field method.
Acknowledgements The financial support of the Italian Ministry for Scientific and Technological Research is gratefully acknowledged.
References [1] G. Alfano, F. Marotti de Sciarra, Mixed finite element formulations and related limitation principles: a general treatment, Comp. Meth. Appl. Mech. Engrg. 138 (1996) 105–130. [2] U. Andelfinger, E. Ramm, EAS-elements for two-dimensional, three-dimensional, plate and shell structures and their equivalence to HR-elements, Int. J. Numer. Meth. Engrg. 36 (1993) 1311–1337. [3] U. Andelfinger, E. Ramm, D. Roehl, 2D and 3D enhanced assumed strain elements and their application in plasticity, in: Proceedings Third Conference on Computational Plasticity (COMPLAS III), vol. II, 1992, pp. 1997–2007. [4] K. Arunakirinathar, B.D. Reddy, Further results for enhanced strain methods with isoparametric elements, Comp. Meth. Appl. Mech. Engrg. 127 (1995) 127–143. [5] M. Bischoff, E. Ramm, D. Braess, A class of equivalent enhanced assumed strain and hybrid stress finite elements, Computat. Mech. 22 (1999) 443–449. [6] H. Brezis, Analyse Fonctionnelle, Theorie et Applications, Masson Editeur, Paris, 1983. [7] J. Cesar de Sa, R. Natal Jorge, Enhanced displacement field on low order elements for incompressible problems, in: Proceedings Fifth Conference on Computational Plasticity (COMPLAS V), vol. I, 1997, pp. 565–570. [8] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [9] B. Fraeijs de Veubeke, Displacement and equilibrium models in the finite element method, in: O.C. Zienkiewicz, G.S. Holister (Eds.), Stress Analysis, Wiley, London, 1965, pp. 145–196. [10] J.B. Hiriart-Urruty, C. Lemarechal, Convex Analysis and Minimization Algorithms, vol. 1, Springer, New York, 1993. [11] T.J.R. Hughes, The Finite Element Method—Linear Static and Dynamic Finite Element Analysis, Dover, New York, 2000. [12] E.P. Kasper, R.L. Taylor, A mixed-enhanced strain method problems, University of California at Berkeley, Report UCB/SEMM97/02, 1997. [13] J.E. Marsden, T.J.R. Hughes, Mathematical Foundations of Elasticity, Prentice-Hall, Englewood Clifs, NJ, 1983. [14] U. Perego, Variationally consistent stress recovery for enhanced strain finite elements, Commun. Numer. Math. Engrg. 16 (2000) 151–163. [15] T.H.H. Pian, K. Sumihara, Rational approach for assumed stress finite elements, Int. J. Num. Meth. Engrg. 20 (1985) 1685–1695. [16] B.D. Reddy, J.C. Simo, Stability and convergence of a class of enhanced strain methods, SIAM J. Numer. Anal. 32 (1995) 1705– 1728. [17] G. Romano, L. Rosati, F. Marotti de Sciarra, A variational theory for finite-step elastoplastic problems, Int. J. Solids Struct. 30 (1993) 2317–2334. [18] G. Romano, F. Marotti de Sciarra, M. Diaco, Well-posedness of three field methods with enhanced strains, in: Proceedings of ECCM’ 99, Munich, Germany, 1999. [19] G. Romano, L. Rosati, M. Diaco, Well-posedness of mixed formulations in elasticity, ZAMM 79 (1999) 435–454. [20] G. Romano, Theory of structural models, Part I, Elements of functional analysis, Doctoral Lectures (in Italian), vol. II, University of Napoli Federico, 2001. [21] G. Romano, F. Marotti de Sciarra, M. Diaco, Well-posedness and numerical performances of the strain gap method, Int. J. Numer. Meth. Engrg. 51 (2001) 103–126. [22] J.C. Simo, M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Meth. Engrg. 29 (1990) 1595–1638. [23] R.L. Taylor, P.J. Beresford, E.L. Wilson, A non-conforming element for stress analysis, Int. J. Numer. Meth. Engrg. 22 (1986) 39– 62. [24] K. Washizu, Variational Methods in Elasticity and Plasticity, third ed., Pergamon Press, New York, 1982. [25] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, vol. 1, fourth ed., McGraw-Hill, London, 1989.
Relations between enhanced strain methods and the HR method Francesco Marotti de Sciarra * Dipartimento di Scienza delle Costruzioni, Facolt a di Ingegneria, Universit a degli Studi di Napoli Federico II, Piazzale Tecchio 85, I80125 Napoli, Italy Received 21 May 2001; received in revised form 14 December 2001
Abstract The three-field Hu–Washizu functional is the starting point to develop a comparative analysis of approximate linear elastostatic problems. In particular we will analyse the limitation phenomena concerning two mixed methods, respectively denoted the strain gap method and the mixed enhanced strain method, and the two-field Hellinger–Reissner (HR) finite element formulation. Accordingly two original limitation principles are contributed. We provide the theoretical motivation of the numerical observations, provided in computational literature, on the coincidence of the displacement and stress solutions between the enhanced assumed strain method and the HR method if different choices of shape functions are made. It is further shown that the HR method provides the exact stress solution and the exact nodal displacements with one element. The related analytical motivation are thus given. The theoretical analyses presented in the paper are supported by numerical examples on tipical benchmarks. Ó 2002 Elsevier Science B.V. All rights reserved.
Keywords: Hu–Washizu method; Hellinger–Reissner method; Enhanced strain methods; Finite element formulations; Limitation principles
1. Introduction The standard displacement finite element method exhibits severe inaccuracies in several instances of engineering interest. Examples of these problems consist in locking phenomena for nearly incompressible materials and in a poor accuracy when the bending behaviour is dominant and low order elements are adopted. In order to improve the performance of low order finite elements and to overcome locking phenomena, enhanced strain methods have been proposed in the computational literature and have gained considerable
*
Tel.: +39-81-768-3336; fax: +39-81-768-3332. E-mail address: [email protected] (F. Marotti de Sciarra).
0045-7825/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 2 ) 0 0 1 8 9 - 5
2662
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
success in the scientific community as evidenced by the great number of papers on the subject in linear and non-linear elasticity and in elastoplasticity [2,3,7,12,21,22]. In this paper we consider the enhanced assumed strain (EAS) method proposed by Simo and Rifai [22] and the mixed enhanced strain (MES) method proposed by Kasper and Taylor [12]. Enhanced strain methods can be formulated starting from the Hu–Washizu functional in which strains, stresses and displacements are assumed as independent fields. Stress and strain fields are required to be square integrable functions and hence can be interpolated independently in each element of the finite element mesh. The corresponding displacement field must, on the contrary, fulfill the interelement continuity and the external constraint conditions. In the original approach of the EAS method [22] the stress fields are assumed to be orthogonal to the enhanced strains and therefore stresses are dropped out from the variational formulation. As a consequence several a posteriori stress recovery strategies, including the discussed question of a stress recovery based on the stress obtained from the constitutive elastic relation, have been proposed by several authors [2,5,14,22]. A detailed analysis on these issues has been carried out in [21] where a new method, called the strain gap method (SGM), has been contributed. In the MES method the strain field is additively decomposed in two sets. The former contains the stress shape functions and the latter contains the enhanced strain shape functions which are assumed to be orthogonal to the stress shape functions. In [18,21] the basic question of well-posedness of discrete mixed enhanced problems has been addressed, a new enhanced method has been proposed and a critical analysis of both EAS and MES methods [12,22] has been carried out. In the present paper we provide a theoretical and numerical analysis of the relations existing between the finite element solutions obtained from the two-field Hellinger–Reissner (HR) method and from the MES method and the SGM. We show that the stress and the enhanced strain shape functions, usually adopted in the MES method [12], fulfill a limitation principle so that the method does in fact collapse into a two-field HR method. The theoretical treatment is supported by the numerical analysis of a plate subjected to simple bending which is usually assumed as a valuable benchmark for enhanced methods, see e.g. [14,16]. The computational analysis here performed reveales that the two-field HR method provides the exact stress field and the exact nodal displacements with just one element mesh. A detailed explanation of this peculiar result is provided. We contribute a sufficient condition which ensures that the stress solution obtained from the two-field HR method is in fact the exact solution. The proof of the coincidence of the nodal displacements with the exact displacements is based on an application of the virtual work principle. In the second part of the paper, the SGM method is considered. We provide a theoretical foundation to the numerical observation made in [2] concerning the coincidence of the displacement and stress solutions obtained from the EAS method and the two-field HR method under usual choices of stress and enhanced strain shape functions. Finally numerical examples of the limitation phenomenon concerning the SGM, the MES method and the HR method are provided.
2. Mixed variational formulations In this paper we will perform the analysis of mixed methods for a continuous linearly elastic structural model defined on a bounded domain X of an Euclidean space with boundary oX and closure X ¼ X [ oX. The kinematic space V is endowed with a suitable Hilbert topology and its dual F is the space of external forces. To each displacement field v 2 V there corresponds a boundary field Cv 2 oV. The dual
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2663
space oF of oV is the space of boundary tractions t and body forces b belong to the Hilbert space H of square integrable vector fields on the structural model. External forces are then collected in the set ‘ ¼ fb; tg 2 F. Strains e and stresses r belong to the Hilbert space H of square integrable vector fields on X [20], the strain and stress spaces will be also denoted by D S H. The structure is subject to homogeneous boundary conditions so that admissible displacements belong to the subspace L of conforming displacements. The kinematic operator B 2 LinfL; Dg associates the strain fields e with the corresponding displacement fields v and it is a continuous linear map from L into D with closed range [6] and with a finite dimensional kernel [19]. The stress fields belong to the Hilbert space S ¼ fr 2 H: B0o r 2 H g where the differential operator 0 Bo 2 LinfS; H g is the formal adjoint of B and provides the body forces corresponding to the stress field r 2 S. For simplicity, we will not consider in the sequel imposed strains and displacements since they can be carried in without any significant modification of the reasonings presented in this paper. The linear elasticity operator E 2 LinfH; Hg is continuous, symmetric and H-elliptic, that is ððEe; eÞÞH P cE kek2H
8 e 2 D:
The elastic strain energy of the structural model is provided by the convex quadratic functional / : D 7! R [ fþ1g given by /ðeÞ ¼ 12ððEe; eÞÞH . Its conjugate convex functional [10] / : S 7! R [ fþ1g represents the complementary elastic energy and is given by / ðrÞ ¼ sup fððr; eÞÞH /ðeÞg ¼ 12ððCr; rÞÞH ; e2D
1
where C ¼ E is the elastic compliance. The structural problem can be written in an operator form as [17] Bu ¼ e;
B0 r ¼ ‘;
r ¼ d/ðeÞ ¼ Ee;
ð2:1Þ
where B0 2 LinfS; Fg is the dual operator of B and is defined as ððr; BrÞÞH ¼ ððB0 r; rÞÞH for any u 2 V and r 2 S. The solutions of the elastostatic problem can be characterized as stationarity points of the Hu–Washizu functional [24]: W ðe; r; uÞ ¼ /ðeÞ þ ððr; Bu eÞÞH h‘; ui;
ð2:2Þ
with u 2 L, e 2 D, r 2 S,where the linear functional h‘; ui ¼ ðb; uÞH þ ht; Cui yields the virtual work of body forces and boundary tractions. Let us now introduce the strain gap field g 2 D as the difference between the compatible strain Bu 2 D and the strain e 2 D in the form g ¼ Bu e. Accordingly the Hu–Washizu functional can be re-written in the form e ðu; g; rÞ ¼ /ðBu gÞ þ ððr; gÞÞ h‘; ui ¼ 1ððEðBu gÞ; Bu gÞÞ þ ððr; gÞÞ h‘; ui W H H H 2
ð2:3Þ
with u 2 L, g 2 D, r 2 S. The expression of the HR functional Rðr; uÞ can be obtained starting from the HW functional (2.2) by enforcing the constitutive relation in terms of the Fenchel’s equality [10]: /ðeÞ þ / ðrÞ ¼ ððr; eÞÞH
2664
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
to get Rðr; uÞ ¼ / ðrÞ þ ððr; BuÞÞH h‘; ui ¼ 12ððr; CrÞÞH þ ððr; BuÞÞH h‘; ui;
ð2:4Þ
with u 2 L and r 2 S.
3. Discrete formulations The unknown fields fe; g; r; ug 2 D D S V can be interpolated over a finite element Xe by means of the matrices Nee ðxÞ, Neg ðxÞ, Ner ðxÞ and Neu ðxÞ, collecting the strain, strain gap, stress and displacement interpolating shape functions, according to the expressions [11]: e
Deh ¼ feeh : eeh ðxÞ ¼ Nee ðxÞpee and pee 2 Rne ; x 2 Xe g; e
Geh ¼ fgeh : geh ðxÞ ¼ Neg ðxÞpeg and peg 2 Rng ; x 2 Xe g; e
Seh ¼ freh : reh ðxÞ ¼ Ner ðxÞper and per 2 Rnr ; x 2 Xe g; e
Veh ¼ fueh : ueh ðxÞ ¼ Neu ðxÞpeu and peu 2 Rnu ; x 2 Xe g: The domain decomposition induced by the meshing of X depends on a parameter h which goes to zero as the finite mesh is refined. The apex e is referred to quantities pertaining to the element Xe . The dimension of the interpolating subspaces Deh , Geh , Seh and Veh has been denoted respectively by nee , neg , ner and neu . The matrices Nee , Neg , Ner and Neu are assumed to have linearly independent columns and the vectors pee , peg , e pr and peu collect the strain, strain gap, stress and displacement parameters pertaining to the element Xe . Let N be the total number of elements of the finite element discretization, then the interpolating spaces are grouped in the following product spaces: Dh ¼
N Y
Deh ;
e¼1
Gh ¼
N Y
Geh ;
Sh ¼
e¼1
N Y
Seh ;
Vh ¼
N Y
e¼1
Veh :
e¼1
Since strain, strain gap and stress fields can be discontinuous at the element interfaces, the corresponding global interpolated fields are simply the collection of the corresponding field defined over each element Xe : eh ¼ fe1h ; e2h ; . . . ; eN h g 2 Dh ;
gh ¼ fg1h ; g2h ; . . . ; gN h g 2 Gh ;
rh ¼ fr1h ; r2h ; . . . ; rN h g 2 Sh :
As a consequence the virtual work performed by an interpolating stress by an interpolating strain or strain gap field is defined as the sum of the elementwise virtual work, i.e. N N Z N Z X X X ððrh ; eh ÞÞ ¼ hhreh ; eeh iiXe ¼ reh eeh ¼ Ner per Nee pee : e¼1
e¼1
Xe
e¼1
Xe
We consider a conforming finite element approximation so that at the element interface the displacement field is continuous. The conforming interpolated displacement fields uh ¼ fu1h ; u2h ; . . . ; uN h g 2 L h Vh satisfy the homogeneous boundary constraints and the interelement continuity conditions. Here Lh is a proper subspace of the product space Vh . The strain, strain gap and stress global parameters qe , qg and qr are the collection of the element parameters pee , peg and per . On the contrary the displacement element parameters peu are expressed in terms of the global parameters qu by means of the finite element assembly operator Aeu .
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2665
3.1. Discrete variational formulations The interpolated counterpart of the Hu–Washizu functional (2.2) can be obtained by introducing the fields feh ; rh ; uh g 2 Dh Sh Lh in the expression of W to get W ðeh ; rh ; uh Þ ¼ 12ððEeh ; eh ÞÞ þ ððrh ; Buh eh ÞÞ h‘; uh i:
ð3:1Þ
A solution feho ; rho ; uho g 2 Dh Sh Lh of the discretized structural model can be obtained by enforcing the stationarity of the Hu–Washizu functional (3.1) to get 8 8 deh 2 Dh ; < ððEeho ; deh ÞÞ ððrho ; deh ÞÞ ¼ 0 ððdrh ; eho ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0 8 drh 2 Sh ; ð3:2Þ : ððrho ; Bduh ÞÞ ¼ h‘; duh i 8 duh 2 Lh : The matrix form of the discrete structural problem can be obtained starting from the stationarity conditions (3.2) and is given by 3 2 3 2 e e e e qe qe H Q O qe JeT JeT O N e H Je e Q Jr X e e 5 4 JeT QeT Je M qr ¼ q ¼ 4 QT O S 5 qr O JeT r e r S Au r eT eT e q e¼1 O ST O qu qu O Au S J r O u o o N X o ¼ o : ð3:3Þ ¼ eT e fu e¼1 A f u u Here Jee , Jer are the canonical extractors which pick up from the global lists qe and qr the local parameters e pee and per pertaining to the element Xe . The operator AeT u transforms the local forces f u into the corresponding global ones f u and the following local matrices have been defined: Z Z Z e e e eT e eT e e Ne ðxÞE ðxÞNe ðxÞ; Q ¼ Ne ðxÞNr ðxÞ; S ¼ NeT H ¼ r ðxÞB Nu ðxÞ Xe
f eu ¼
Z Xe
NeT u ðxÞbðxÞ þ
Xe
Z
Xe
ðCNeu ÞT ðxÞtðxÞ;
oXe
where B and E denotes the matrix form of the operators B and E. The SGM [21] is obtained starting from the interpolated counterpart of the Hu–Washizu functional e h ðuh ; gh ; rh Þ. In fact by adding up the contributions of each non-assembled element and imposing that the W interpolating displacement uh satisfies the conformity requirement we get e h ðuh ; gh ; rh Þ ¼ 1ððEðBuh gh Þ; Buh gh ÞÞ þ ððrh ; gh ÞÞ h‘; uh i W 2
ð3:4Þ
where fuh ; gh ; rh g 2 Lh Gh Sh . The solutions are obtained by enforcing the stationarity of (3.4) to get 8 < ððEBuho ; Bduh ÞÞH ððEgho ; Bduh ÞÞH ¼ h‘; duh i 8 duh 2 Lh ; ð3:5Þ ððEgho ; dgh ÞÞH þ ððrho ; dgh ÞÞH ððEBuho ; dgh ÞÞH ¼ 0 8 dgh 2 Gh ; : ððdrh ; gho ÞÞH ¼ 0 8 drh 2 Sh : e h ¼ Gh \ S? . The strain gaps belonging The third equation of relation (3.5) amounts to require that gh 2 G h ? e to the subspace G h ¼ Gh \ Sh are referred to as effective strain gaps since they effectively contribute to relax the compatibility condition.
2666
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
The matrix form of the discrete problem is given by 2 3 2 qu X qu AeT Ke Aeu AeT GeT Jeg O N K u u 6 eT e eT e e eT e e 7 e f Jg L Jg Jg P Jr 5 qg ¼ 4 G M qg ¼ 4 Jg G Au q e¼1 q eT e O O JeT O r r r P Jg eT e fu N Au f u X ¼ o ¼ o : o e¼1 o
3 O qu P 5 qg O qr
GT L PT
ð3:6Þ
The component submatrices appearing in (3.6) are defined by Z Z e e e Le ¼ NeT ðxÞE ðxÞN ðxÞ; P ¼ NeT g g g ðxÞNr ðxÞ; Xe Xe Z Z T e e Ge ¼ NeT ðxÞE ðxÞB N ðxÞ; K ¼ ðB Neu Þ ðxÞE ðxÞB Neu ðxÞ: g u Xe
ð3:7Þ
Xe
Remark 1. In [21] the SGM has been proposed on the basis of a suitable variant of the HW functional. It is then performed the analysis of the well-posedness conditions which ensure uniqueness of the discrete solution in terms of displacements, strain gaps and stresses. It is shown that the SGM can be split in the reduced problem, formulated in terms of displacements and strain gaps, and in a stress recovery problem depending on the solution of the reduced problem. The comparison between the EAS method and the SGM shows that the EAS method is partially ill-posed since an orthogonality condition drops out the stress parameters. As a consequence the EAS method coincides only with the reduced problem of the SGM and the stress recovery remains, for the EAS method, an open problem. The expression of the HR functional Rðr; uÞ in terms of interpolating fields is Rh ðrh ; uh Þ ¼ / ðrh Þ þ ððrh ; Buh ÞÞ h‘; uh i
rh 2 Sh ; uh 2 Lh :
ð3:8Þ
An approximate solution frho ; uho g 2 Sh Lh is obtained by enforcing the stationarity of Rh to get: ððdrh ; Crho ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0 8 drh 2 Sh ; ð3:9Þ ððrho ; Bduh ÞÞ h‘; duh i ¼ 0 8 duh 2 Lh : The matrix form of the discrete problem is given by X N qr K JeT ReT Jer JeT Te Aeu qr r r A ¼ ¼ eT eT e qu q TT A T J O u u r e¼1 where Re ¼
Z Xe
1 e NeT r ðxÞE ðxÞNr ðxÞ;
Te ¼
Z
e NeT r ðxÞB Nu ðxÞ:
T O
qr qu
¼
N X AeT f e u
e¼1
o
u
fu ; ¼ o
ð3:10Þ
ð3:11Þ
Xe
From the computational point of view it is more convenient to carry out the assembly operation after the condensation, at the element level, of the strain and stress parameters to put the global discrete problem in terms of the sole displacement parameters qu . The procedure has been illustrated in [18,21]. 4. Mixed enhanced strain method The mixed enhanced strain (MES) method [12] is based on the following decomposition of the interpolating strain space:
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
eh D h ¼ Sh D
2667
ð4:1Þ
e h is the subspace of enhanced strains ~eh given by where D eh ¼ D
N Y
e e; D h
e e ¼ f~ee : ~ee ðxÞ ¼ Ne ðxÞpe and pe 2 Rn~ee ; x 2 Xe g: D ~e ~e ~e h h h
e¼1
e h ¼ fog so that stress and enhanced strain The symbol denotes the direct sum and requires that Sh \ D shape functions must be linearly independent. e h , the interpolated counterpart of the Hu–Washizu funcSetting eh ¼ eh þ ~eh with eh 2 Sh and ~eh 2 D tional W in terms of (eh ; ~eh ; rh ; uh ) is given by Mh ðeh ; ~eh ; rh ; uh Þ ¼ 12ððEeh ; eh ÞÞ þ 12ððE~eh ; ~eh ÞÞ þ ððEeh ; ~eh ÞÞ þ ððrh ; Buh eh ~eh ÞÞ h‘; uh i; The MES method follows from the stationarity conditions of Mh and is given by 8 ððEeho ; deh ÞÞ þ ððE~eho ; deh ÞÞ ððrho ; deh ÞÞ ¼ 0 8 deh 2 Sh ; > > < e h; ððEeho ; d~eh ÞÞ þ ððE~eho ; d~eh ÞÞ ððrho ; d~eh ÞÞ ¼ 0 8 d~eh 2 D > ððdrh ; eho þ ~eho ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0 8 drh 2 Sh ; > : ððrho ; Bduh ÞÞ ¼ h‘; duh i 8 duh 2 Lh :
uh 2 L h :
ð4:2Þ
4.1. The MES method vs the HR principle A sufficient condition for the validity of a limitation phenomenon between the MES method (4.2) and the HR method (3.9) is reported in the next statement. Proposition 4.1. (A limitation principle). The solutions of the approximate problems obtained from the MES and the HR methods coincide in terms of rh and uh if the same interpolating spaces Sh and Lh are assumed and if e h: CSh Sh D
ð4:3Þ
e h Sh Lh be a stationarity point for the MES funcProof. Let the quartet feho ; ~eho ; rho ; uho g 2 Sh D tional. The inclusions (4.3) ensure that eh þ ~eh Crh 2 Dh
e h ; eh ; rh 2 Sh : 8~eh 2 D
Summing up the first and second equation of (4.2) which represents the constitutive relations, and by virtue of the positive definiteness of the elastic stiffness E, we have: Eðeho þ ~eho Þ rho 2 D? h \ EDh ¼ fog; e h . The orthogonality relation ? is since Dh is decomposed in the direct sum of the subspaces Sh and D intended according to the inner product in L2 ðXÞ. Then the interpolating strain eho ¼ eho þ ~eho and the interpolating stress rho fulfil the constitutive relation: eho þ ~eho ¼ Crho which inserted in the second relation of (4.2) yields ððdrh ; Crho ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0 8 drh 2 Sh :
2668
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
The MES method (4.2) can then be expressed by the following reduced set of equations:
ððdrh ; Crho ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0 8 drh 2 Sh ; ððrho ; Bduh ÞÞ ¼ h‘; duh i 8 duh 2 Lh ;
which coincide with the stationarity conditions (3.9) of the HR functional.
It is worth noting that the condition (4.3) is local since interpolating strains and stresses are defined elementwise so that we can restrict ourselves to analyse a single element. In a finite element analysis, it is compelling to verify the condition (4.3) in terms of quantities pertaining to the reference element K. To this end we denote by ve : K 7! Xe the isoparametric map from the reference element K to the generic element Xe of the finite element mesh and we set v ¼ fve ; e ¼ 1; . . . ; Ng. Since e h ¼ vð D e K Þ, to check the inclusion (4.3) it is sufficient to verify Sh ¼ vðSK Þ and D e K: CSK SK D
ð4:4Þ
Let us now show that the strain and stress shape functions adopted for the MES method, see [12,18], meet the above limitation principle for distorted and undistorted elements. In the MES method, the stress subspace for 4-node quadrilateral elements with standard bilinear shape functions (Q1), where K ¼ , are given by [15] 2 3 1 0 0 g 0 S ¼ span4 0 1 0 0 n 5 ð4:5Þ 0 0 1 0 0 and the strain subspace is decomposed, according to (4.1), in the form 2 3 n 0 e where D e ¼ span4 0 g 5: D ¼ S þ D 0 0 For plane stress problems, 2 1 14 m CS ¼ span E 0
the subspace CS becomes 3 m 0 g mn 1 0 mg n 5; 0 2ð1 þ mÞ 0 0
ð4:6Þ
ð4:7Þ
where E and m denote the Young’s modulus and the Poisson’s ratio. A direct inspection shows that the base vectors of the subspace CS , given by the columns of the matrix (4.7), can be obtained as a linear combination of the vectors spanning D , see (4.5) and (4.6). The condition (4.4) is then fulfilled, a limitation phenomenon occurs and no advantage is gained by the use of the MES method in place of the simpler HR method. Additional transformation rules which preserve the point-wise inner product between stress and strain tensors can be envisaged but it seems that they can only be motivated by an a posteriori evaluation of the quality of the numerical results. Examples are provided by the push/pull transformations of differential geometry [13] which provides, for plane problems see e.g. [12,22], the expressions reh ðxÞ ¼ Teo re ½v1 e ðxÞ;
geh ðxÞ ¼
Joe e 1 TeT o g ½ve ðxÞ; e J ½v1 ðxÞ e
ð4:8Þ
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2669
where Teo is the value, at the origin of the reference element, of the matrix field 2 2 3e F11 F122 2F11 F12 5 ðnÞ with n 2 K: Te ðnÞ ¼ 4 F212 F222 2F21 F22 F11 F21 F12 F22 F11 F22 þ F12 F21 Accordingly the inner product in the real space is given by Z Z Z Joe eT e 1 e T reh ðxÞ geh ðxÞ dx ¼ Teo re ½v1 ðxÞ g ½v ðxÞ dx ¼ J rðnÞ gðnÞ dn: o e e o J e ½v1 Xe Xe K e ðxÞ Adopting the interpolations (4.8), the condition (4.3) is not fulfilled in the case of distorted elements so that the MES method and the HR method provide different results. A numerical comparative analysis is provided in Section 6.
5. Limitation phenomena between SGM and HR Let us now show sufficient conditions to be fulfilled in order that the SGM method collapses into the HR method. Proposition 5.1. (Limitation principle between the SGM and HR method). The approximate solution fuh ; rh g provided by the SGM (3.5) and by the two-field HR method (3.9) coincide if the same subspaces Lh and Sh are adopted and if the following condition is met: e h: EBLh Sh G
ð5:1Þ
e h , the condition (5.1) implies that EðBuh gh Þ rh 2 Gh . By Proof. By virtue of the choice Gh ¼ Sh G virtue of the second condition of (3.5) we have EðBuh gh Þ rh 2 Gh [ G? h ¼ fog: Accordingly we have rh ¼ EðBuh gh Þ. The constraint condition provided by the third equation of (3.5) e h , so that we have ensures that the strain gaps are effective, i.e. gh 2 G e h S? gh ¼ Buh Crh 2 G h
()
ððdrh ; Buh ÞÞ ððdrh ; Crh ÞÞ ¼ 0
8 drh 2 Sh
ð5:2Þ
and the first equation of (3.5) yields ððrh ; Bduh ÞÞ ¼ h‘; duh i 8 duh 2 Lh which coincide with the stationarity conditions of the HR method (3.9). The statement is thus proved.
ð5:3Þ
An analogous limitation principle between the EAS method and the HR method has been proved in [5] by adding the further hypothesis that the approximate stress pertaining to the EAS method is obtained by means of the constitutive elastic relation. The above Proposition 5.1 shows that such an assumption is inessential and can thus be avoided. In a finite element analysis, the condition (5.1) deserves a carefull discussion since it involves the conforming displacement subspace Lh and the kinematic operator B. To this end we note that this condition is global due to the presence of Lh which depends on the a priori unknown element assembly operations. This shortcoming can be circumvented by considering the larger non-conforming displacement space Vh Lh so that the global condition (5.1) can be substituted by the e h. sufficient local condition EBVh Sh G
2670
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
The kinematic operator acting on the displacement fields in the reference element K will be denoted by BeK and is defined by BeK uK ðnÞ ¼ Bueh ðxÞ for any x 2 Xe where ueh is the restriction of uh to the element Xe e and n ¼ v1 e ðxÞ is the isoparametric coordinate in the reference element. We then set BK ¼ fBK ; e ¼ 1; . . . ; Ng. Denoting by gK ¼ v1 ðuh Þ the collection of the displacement fields obtained by shifting back to K the e h can then be fields ueh , we have Buh ¼ Bv1 ðuh Þ ¼ BK gK . The sufficient condition above EBVh Sh G written in the reference element K in the form eK: EBK VK SK G
ð5:4Þ
It is apparent that the kinematic operator BK on K depends on the map v and hence the condition (5.4) cannot be checked in terms of quantities pertaining only to the reference element K. These difficulties have not been fully realized in the literature since in previous treatments, see e.g. [2,4,22,25], conditions of this kind has been simply imposed in the reference element K. In fact for a general isoparametric map, the condition (5.4) can be only checked in the form eK EBVK SK G
ð5:5Þ
and confide that the element distorsion does not make it fail on the real elements. On the basis of the previous analysis, let us now provide a theoretical explanation of the numerical observation made in [2] that the EAS method [22] and the HR method yield the same solution in terms of displacements and stresses. e h and Sh with reference to We analyse the different choices proposed in [2] concerning the subspaces G quadrilateral Q1 elements. We consider the Pian and Sumihara stress interpolation and the seven-parameter strain gap interpolation: 2 3 2 3 1 0 0 g 0 n 0 0 0 ng 0 0 e ¼ span4 0 g 0 0 0 ng 0 5: S ¼ span4 0 1 0 0 n 5; G ð5:6Þ 0 0 1 0 0 0 0 n g 0 0 ng Stress and strain gap shape functions in the isoparametric space belong to the four-dimensional space Q1 e is direct so that the subspace S G e generated by the polynomials 1, n, g, and ng [8]. The sum S þ G turns out to be twelve-dimensional. Since the compatible strain subspace BV is given by 2 3 1 0 0 g 0 BV ¼ span4 0 1 0 0 n 5; 0 0 1 n g the condition (5.5) is fulfilled. As a further example we consider the eight-parameter stress interpolation together with a four-parameter enhanced strain (EAS4) interpolation which coincides with the modified incompatible mode approximation of Taylor et al. [23]: 2 3 2 3 1 0 0 g 0 ng 0 0 n 0 0 0 e ¼ span4 0 g 0 0 5: S ¼ span4 0 1 0 0 n 0 ng 0 5; G ð5:7Þ 0 0 1 0 0 0 0 ng 0 0 n g Finally the following twelve-parameter stress interpolation and no enhanced strains were considered in [2]:
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2
1 S ¼ span4 0 0
0 1 0
0 0 1
g 0 0
0 ng n 0 0 0
0 ng 0
0 0 ng
n 0 0 g 0 0
0 0 n
3 0 0 5; g
e ¼ f0g: G
2671
ð5:8Þ
It can be proved, following the arguments shown above, that the choices (5.7) and (5.8) fulfill the condition reported in the Proposition 5.1 so that the SGM collapse into the HR method. It is worth noting that the choice (5.8) of the shape functions fulfils the inclusion EBVK SK which ensures that the displacement solution obtained from the HR method and the standard displacement method coincides [1,9]. Let us now analyse the stress and strain gap shape functions considered in [22]: 2 3 2 3 1 0 0 g 0 n 0 0 0 ng e ¼ span4 0 g 0 0 ð5:9Þ S ¼ span4 0 1 0 0 n 5; G ng 5: 0 0 1 0 0 0 0 n g n2 g2 Following the arguments previously reported it is immediate to prove that the condition (5.5) is fulfilled and the limitation phenomenon occurs.
6. Numerical examples The limitation principle reported in the Proposition 4.1 is tested by means of a simple structural example consisting in a plate of unit thickness subjected to simple bending as shown in Fig. 1. This example is usually adopted as a benchmark for enhanced methods [4,14]. The analytical solution in terms of horizontal and vertical displacements is given by 2p uðx; yÞ ¼ Eb x b2 y ; p vðx; yÞ ¼ Eb ½x2 þ mðy 2 byÞ: The only non-vanishing stress component is rx whose expression is given by 2p b rx ðx; yÞ ¼ y : b 2 Setting a ¼ 10, b ¼ 2, E ¼ 1500, m ¼ 0:25, p ¼ 3000, the analytical solution of the bending problem in terms of displacements and stresses is uðx; yÞ ¼ 2xð1 yÞ, vðx; yÞ ¼ x2 þ 14ðy 2 2y) and rx ðx; yÞ ¼ 3000ð1 yÞ.
Fig. 1. The bending plate problem discretized by means of Q1 elements.
2672
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
A first example is carried out for two meshes of 4-node quadrilateral elements with bi-linear displacement interpolations (Q1): 1 1 rectangular element and 16 8 uniform rectangular elements. The Pian and Sumihara stress shape functions (4.5) and the strain shape functions (4.6) are adopted. The horizontal and vertical displacements and the horizontal stress at the corner node A of the plate are evaluated adopting the MES method, the HR method and the standard displacement formulation and they are compared with the solution of the continuous problem. The results are reported in Table 1. A further example is devoted to analyse how the distortion of the elements affects the discrete solution. A classical benchmark is the plate reported in Fig. 1 when it is modelled with two distorted quadrilateral elements. Figs. 2 and 3 respectively show the vertical displacement of the point A and the horizontal stress rx at the point B against the distortion parameter d. The shape functions (4.5) and (4.6) are adopted for the MES and the HR methods. Figs. 2 and 3 show that the displacement and the stress solutions obtained from the MES method and the HR formulation are identical as predicted by the limitation principle 4.1. The MES-p and the HR-p plots of Figs. 2 and 3 are obtained by including the push-transformation (4.8) in the definition of stress and strain shape functions.
Table 1 The displacement components u and v and the horizontal stress rx at the node A for the following Q1 meshes (panel a) one element mesh; (panel b) 16 8 uniform rectangular mesh MES
HR
Panel a u v rx
DISPL
EXACT
20 100 3000
20 100 3000
2 9 289
20 100 3000
Panel b u v rx
20 100 3000
20 100 3000
19 96 2815
20 100 3000
Fig. 2. Vertical displacement of the node A for a cantilever rectangular plate meshed with two distorted elements.
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2673
Fig. 3. Horizontal stress at the node B for a cantilever rectangular plate meshed with two distorted elements.
The resulting shape functions on the real element do not meet the limitation principle and, in fact, the solutions obtained from these two methods are different. In particular the push-transformation has beneficial effects for the MES method in terms of stresses but it shows a poor performance in terms of displacements for increasing values of d. The displacement and stress solutions obtained from the SGM and the HR method do coincide for undistorted elements if the shape functions (5.9) are used. For distorted elements, the presence of the kinematic operator in the expression (5.1) implies that the interpolating subspaces do not fulfil the inclusion appearing in (5.1) so that the two methods provide different solutions, see Figs. 2 and 3. The examples reveal that, assuming a one element mesh, the interpolating stress rh coincide with the exact stress solution. Moreover, even if the interpolating displacement field uh is quite different from the exact solution, the displacement parameters qu are exact. In the next section we will explain in detail this result.
7. The finite element exact solution Let us now provide a sufficient condition which ensures that the stress solution rho obtained from the HR method coincides with the exact stress solution ro pertaining to the continuous model. In particular we will show that this condition is met by the displacement and stress shape functions adopted for the finite element analysis of the bending problem (Fig. 1). Further we will provide a theoretical motivation to the numerical observation that the nodal displacements qu obtained by the HR method are exact. To this end we prove the next result. Proposition 7.1. (The exact stress solution) Let fro ; uo g 2 S L be the solution of the continuous problem and frho ; uho g 2 Sh Lh be the approximate solution of the HR method.
2674
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
Then if r o 2 Sh
and
Buo 2 BLh S? h
ð7:1Þ
the approximate and the exact stress solutions coincide, i.e. rho ¼ ro . Proof. The exact stress solution ro 2 S fulfils the approximate equilibrium equation (3.9), part 2 ððro ; Bduh ÞÞ h‘; duh i ¼ 0
8 duh 2 Lh
ð7:2Þ
since the inclusion Lh L hold true. It remains to prove that the exact stress solution ro satisfies the approximate elastic compatibility equation (3.9), part 1 which is hereafter rewritten for convenience: ððdrh ; Cro ÞÞ þ ððdrh ; Buho ÞÞ ¼ 0
8 drh 2 Sh :
ð7:3Þ
To this end we note that second condition of (7.1) implies that there exists a unique interpolating conforming displacement uh 2 Lh such that Buo þ Buh 2 S? h
()
ððdrh ; Buo ÞÞ þ ððdrh ; Buh ÞÞ ¼ 0
8 drh 2 Sh :
Further, noting that the equalities Cro ¼ eo ¼ Buo hold true for the exact solution, the orthogonality relation above can be equivalently rewritten in the form ððdrh ; Cro ÞÞ þ ððdrh ; Buh ÞÞ ¼ 0
8 drh 2 Sh :
ð7:4Þ
fro ; uh g
Accordingly, the pair fulfils the approximate equilibrium equation (7.2) and the approximate elastic compatibility (7.4) and hence fro ; uh g is the solution of the HR approximate problem (3.9). By virtue of the uniqueness of the approximate solution we infer that fro ; uh g fro ; uho g so that ro is the stress solution of the approximate problem. The proposition is thus proved. 7.1. The stress solution of the plate bending problem Let us now show that a one element mesh provides the exact stress solution for the linear elastic bending problem presented in the previous section. For a one element Q1 mesh, the set of shape functions spanning the subspace BLh is given by 2 3 1 0 0 g 0 ð7:5Þ BLh ¼ span BNu ¼ span4 0 1 0 0 n 5: 0 0 1 n g The compatible strain Buo 2 D of the continuous problem is given by g 1 Buo ¼ mg E 0
ð7:6Þ
and can be obtained in terms of a linear combination of the columns of the matrix (7.5) as follows: g 0 0 1 1 m Buo ¼ 0 0 g : ð7:7Þ E E E n n 0 The last two vectors appearing in (7.7) are L2 ðÞ––orthogonal to the subspace S generated by the assumed Pian–Sumihara stress interpolation (4.5) so that they belong to S? .
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2675
Noting that the stress solution ro of the continuous problem belong to the interpolating stress subspace S : g ro ¼ 0 2 S ; 0 the requirements of the Proposition 7.1 are fulfilled so that a one element mesh provides the exact stress solution. 7.2. The displacement solution of the plate bending problem Let us now explain why the interpolation fields, assumed for the analysis of the plate bending problem reported in Fig. 1, provide the exact nodal displacements for a one element mesh. Since Crho ¼ Cro ¼ Buo , the elastic compatibility condition (3.9), part 1 can be rewritten as ððdrh ; Buo ÞÞ ¼ ððdrh ; Buho ÞÞ;
8 drh 2 Sh :
ð7:8Þ
The five-dimensional Pian–Sumihara interpolating stress subspaces Sh is spanned by the five columns of the matrix (4.5) which will be denoted by si with i ¼ 1; . . . ; 5. The elastic compatibility relation (7.8) can then be tested in terms of the base vectors fsi g of Sh in the form: ððsi ; Buo ÞÞ ¼ ððsi ; Buho ÞÞ;
i ¼ 1; . . . ; 5:
ð7:9Þ
Since the exact compatible strain solution Buo is given by (7.6), it is apparent that the l.h.s. of (7.9) is vanishing for i ¼ 1, 2, 3, 5. Hence the equality (7.9) yields ððsi ; Buho ÞÞ ¼ 0; i ¼ 1; 2; 3; 5 and
ððsi ; Buho ÞÞ 6¼ 0; i ¼ 4:
ð7:10Þ
The interpolating compatible strain Buho is expressed by a linear combination of the base vectors of BLh , see (7.5), so that the relations (7.10) imply that the deformation modes provided by the columns 1, 2, 3 and 5 of (7.5) cannot be active. In conclusion the hourglass mode in the x-direction provides the only compatible deformation of the approximate model, see Fig. 4. By virtue of the equality (7.9), the virtual work principle allows us to write ððs4 ; Bðuo uho ÞÞÞ ¼ MðDuo Duho Þ ¼ 0
ð7:11Þ
where M denotes the external couple in equilibrium with the stress distribution s4 and Duo , Duho represent the exact and the approximate relative rotations of the beam end sections. The equality (7.11) implies Duo ¼ Duho and, being the left end section fixed, the absolute rotation u of the right end section is exact. It is well known that the curvature of the continuous model is constant so that the displacement of the centroid of the right end section is 12ul (see Fig. 5a). The approximate model provide the same result since
Fig. 4. Hourglass mode in the x-direction.
2676
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
Fig. 5. Deflection of the beam (a) for the continuous model; (b) adopting a one element Q1 analysis.
the rotation to be superimposed to the hourglass mode to restore the external constraints is 12u so that the displacement of the centroid of the right end section is 12ul (see Fig. 5b).
8. Closure The analysis carried out in this paper provide a guide to check if the shape functions to be adopted in the MES method and in the SGM meet a limitation phenomenon so that the same displacement and stress solutions of the HR method are obtained. The SGM has been recently contributed by the authors and includes the EAS method as a special ill-posed case. In particular we have shown that a limitation principle between the MES method and the HR method is fulfilled by the shape functions usually adopted for the interpolating stresses and enhanced strains for undistorted meshes. In this case no advantage is obtained by the adoption of an enhanced method with respect to the two-field HR formulation. Further a limitation principle has been provided between the SGM method and the HR formulation. We have provided a theoretical motivation to the numerical observations made by Andelfinger and Ramm which have shown the coincidence of the displacement and stress solutions between the SGM method and the HR method for several choices of stress and enhanced strain shape functions. The proposed argumentations are based on the properties of the interpolating subspaces and can be easily repeated if different stress and enhanced strain shape functions are considered. Two further results have been pointed out in this paper with reference to the HR method. The former is a sufficient condition, expressed in terms of interpolating subspaces, which ensures that the interpolating stress obtained from the HR method coincides with the exact stress solution. This condition is met by the bending problem of the plate so that the exact stress solution is obtained by means of the HR method with one element. The latter result is the theoretical motivation of the coincidence of the exact displacements with the nodal displacements obtained from the HR method with one element. The theoretical considerations have been supported by the computational analysis of a plate subjected to simple bending, which is usually considered as a benchmark for enhanced methods, and the influence on the solution of the distortion of the elements has also been pointed out. Accordingly, the plate bending problem cannot be considered as an effective benchmark for enhanced methods since the HR method provides the exact solutions with one element. Enhanced methods can then
F. Marotti de Sciarra / Comput. Methods Appl. Mech. Engrg. 191 (2002) 2661–2677
2677
share an analogous performance if a limitation phenomenon occur otherwise they behave worse than the two-field method.
Acknowledgements The financial support of the Italian Ministry for Scientific and Technological Research is gratefully acknowledged.
References [1] G. Alfano, F. Marotti de Sciarra, Mixed finite element formulations and related limitation principles: a general treatment, Comp. Meth. Appl. Mech. Engrg. 138 (1996) 105–130. [2] U. Andelfinger, E. Ramm, EAS-elements for two-dimensional, three-dimensional, plate and shell structures and their equivalence to HR-elements, Int. J. Numer. Meth. Engrg. 36 (1993) 1311–1337. [3] U. Andelfinger, E. Ramm, D. Roehl, 2D and 3D enhanced assumed strain elements and their application in plasticity, in: Proceedings Third Conference on Computational Plasticity (COMPLAS III), vol. II, 1992, pp. 1997–2007. [4] K. Arunakirinathar, B.D. Reddy, Further results for enhanced strain methods with isoparametric elements, Comp. Meth. Appl. Mech. Engrg. 127 (1995) 127–143. [5] M. Bischoff, E. Ramm, D. Braess, A class of equivalent enhanced assumed strain and hybrid stress finite elements, Computat. Mech. 22 (1999) 443–449. [6] H. Brezis, Analyse Fonctionnelle, Theorie et Applications, Masson Editeur, Paris, 1983. [7] J. Cesar de Sa, R. Natal Jorge, Enhanced displacement field on low order elements for incompressible problems, in: Proceedings Fifth Conference on Computational Plasticity (COMPLAS V), vol. I, 1997, pp. 565–570. [8] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [9] B. Fraeijs de Veubeke, Displacement and equilibrium models in the finite element method, in: O.C. Zienkiewicz, G.S. Holister (Eds.), Stress Analysis, Wiley, London, 1965, pp. 145–196. [10] J.B. Hiriart-Urruty, C. Lemarechal, Convex Analysis and Minimization Algorithms, vol. 1, Springer, New York, 1993. [11] T.J.R. Hughes, The Finite Element Method—Linear Static and Dynamic Finite Element Analysis, Dover, New York, 2000. [12] E.P. Kasper, R.L. Taylor, A mixed-enhanced strain method problems, University of California at Berkeley, Report UCB/SEMM97/02, 1997. [13] J.E. Marsden, T.J.R. Hughes, Mathematical Foundations of Elasticity, Prentice-Hall, Englewood Clifs, NJ, 1983. [14] U. Perego, Variationally consistent stress recovery for enhanced strain finite elements, Commun. Numer. Math. Engrg. 16 (2000) 151–163. [15] T.H.H. Pian, K. Sumihara, Rational approach for assumed stress finite elements, Int. J. Num. Meth. Engrg. 20 (1985) 1685–1695. [16] B.D. Reddy, J.C. Simo, Stability and convergence of a class of enhanced strain methods, SIAM J. Numer. Anal. 32 (1995) 1705– 1728. [17] G. Romano, L. Rosati, F. Marotti de Sciarra, A variational theory for finite-step elastoplastic problems, Int. J. Solids Struct. 30 (1993) 2317–2334. [18] G. Romano, F. Marotti de Sciarra, M. Diaco, Well-posedness of three field methods with enhanced strains, in: Proceedings of ECCM’ 99, Munich, Germany, 1999. [19] G. Romano, L. Rosati, M. Diaco, Well-posedness of mixed formulations in elasticity, ZAMM 79 (1999) 435–454. [20] G. Romano, Theory of structural models, Part I, Elements of functional analysis, Doctoral Lectures (in Italian), vol. II, University of Napoli Federico, 2001. [21] G. Romano, F. Marotti de Sciarra, M. Diaco, Well-posedness and numerical performances of the strain gap method, Int. J. Numer. Meth. Engrg. 51 (2001) 103–126. [22] J.C. Simo, M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Meth. Engrg. 29 (1990) 1595–1638. [23] R.L. Taylor, P.J. Beresford, E.L. Wilson, A non-conforming element for stress analysis, Int. J. Numer. Meth. Engrg. 22 (1986) 39– 62. [24] K. Washizu, Variational Methods in Elasticity and Plasticity, third ed., Pergamon Press, New York, 1982. [25] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, vol. 1, fourth ed., McGraw-Hill, London, 1989.