INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 15—27 (1997) A NEW FORMULATION OF ISOPARAMETRIC FINITE ELEMENTS AND THE RELATIONSHIP BETWEEN HYBRID STRESS ELEMENT AND INCOMPATIBLE ELEMENT ZHAO-PING JIAO Department of Civil Engineering, South China Construction ºniversity, (¼est Campus), 510405 Guangzhou, China THEODORE H. H. PIAN Department of Aeronautics and Astronautics, Massachussets Institute of ¹echnology, Cambridge, MA, º.S.A. SHENG YONG Department of Modern Mechanics, ºniversity of Science and ¹echnology of China, Hei Fei, China SUMMARY A new method of formulating isoparametric finite element is developed, and the element strains are proposed to be resolved into two parts, constant part and higher-order one. The new method indicates two important properties of isoparametric finite element, and the equivalent relationship between hybrid stress elements and incompatible elements. KEY WORDS: formulation; finite elements; isoparametric; hybrid stress; incompatible INTRODUCTION Pian1 pointed out that the incompatible element presented by Wilson et al.2 and the hybrid stress element by Pian3 are equivalent in the special case of rectangular plane element. However, it is difficult to prove the equivalence between these two methods when elements are distorted. The main obstruction of looking for the correlation of the above two methods lies in that the stress model of incompatible elements cannot be separated into two independent sections, constant stresses and higher-order ones, as well as the stress field of hybrid elements. In this paper, a new formulation of quadrilateral plane element Q4 is constructed to be equivalent to the traditional one by means of decomposing the element displacement field into two independent portions, one part contains constant strains of element, another part only includes higher-order strains. By virtue of the new formulation, it becomes very convenient to form the new kinds of incompatible elements model whose behaviour depends on the element characteristic matrix W*. When the matrix W* is designed according to the higher-order stress model of hybrid stress element, the incompatible element which is equivalent to hybrid stress elements can be obtained. On the other hand, for an arbitrary matrix W*, we can always construct a hybrid stress element mode equivalent to the incompatible element, when the higher-order stress model of hybrid stress elements is designed in accordance with the matrix W*. Therefore, the research works mentioned above demonstrate the equivalent correlation between hybrid CCC 0029—5981/97/010015—13 ( 1997 by John Wiley & Sons, Ltd. Received 29 September 1995 Revised 21 May 1996 16 Z.-P. JIAO, T. H. H. PIAN AND S. YONG stress elements and incompatible elements which contain the special case of the compatible element. NEW FORMULATION OF 4-NODE QUADRILATERAL ELEMENT AND ELEMENT SPECIAL PROPERTIES The element displacement u and v of a general quadrilateral element (Figure 1) are interpolated q q as 4 4 u " + NM u , v " + NM v q i i q i i i/1 i/1 (1) NM "1 (1#m m) (1#g g), i"1, 2, 3, 4 i 4 i i (2) where the shape functions are when C D GH a a a 1 4 xN 1 2 3 " + i [f f g g ] i i i i b b b y 4 1 2 3 i/1 i (3) The co-ordinates x and y can be expressed as x"a m#a g#a fg, y"b f#b g#b fg 1 3 2 1 3 2 (4) The determinant of Jacobian is DJD"(a b !a b )#(a b !a b )f#(a b !a b )g"A#Bf#Cg 1 3 3 1 1 2 2 1 2 3 3 2 (5) For example, the displacement u can be written in the form q u "[1 f g fg]a@"[1 x y fg]a q (6) in which equation (4) has been used. Comparing equation (6) with equation (1), we have 4 1 a #(a f#a g#a fg)a #(b f#b g#b fg)a #fga " + (1#f f#g g#f g fg)u i i i i i 1 1 3 2 2 1 3 2 3 4 4 i/1 (7) Figure 1. Quadrilateral element A NEW FORMULATION OF ISOPARAMETRIC FINITE ELEMENTS 17 From equation (7), a (i"1, 2, 3, 4) can be solved and the displacement u can be expressed in i q terms of nodel displacements 4 u "+ N u q i i i/1 (8a) 4 v "+ N v q i i i/1 (8b) N "1#A x#B y, N "C fg i i i) i i# 4 1 1 A " (b f !b g ), B " (a g !a f ) i 4A 3 i 1 i i 4A 1 i 3 i 1 C " (Af g !Bg !Cf ) i i i i i 4A (9) Similarly, where N "N #N i i# i) The element displacements can be expressed again with shape functions GH C DG H u u 4 N 0 4 i i " + N q "Nq"u #u (10) u " q "+ i i q# q) q v 0 N v q i i i/1 i/1 The element displacements in equations (10) are equivalent to that in equation (1), the only differences is that the displacement field here has been decomposed into two independent parts, u and u : qc qh 0 4 N 4 i# u "+ q " + N q "N q (11a) q# i i# i i# 0 N i# i/1 i/1 and C D C D 0 4 4 N i) q " + N q "N q (11b) u "+ i i) i ) q) 0 N i) i/1 i/1 From equation (9), we find that the displacements u are constructed by complete linear q# functions which contain the rigid-body displacements and the constant strains of element, whereas there is only higher-order composition fg in u . q) According to equation (10), it is very easy to calculate strain submatrix B "B #B , i"1, 2, 3, 4 i i# i) (12a) where A 0 i B "LN " 0 B , i# i# i B A i i t 0 1 B "LN "C 0 t "C ( 2 i) i) i i t t 2 1 (12b) 18 Z.-P. JIAO, T. H. H. PIAN AND S. YONG in which C L" L/Lx 0 0 L/Ly D L/Ly T L/Lx and G H G H G t L/Lx 1 b g!b f 1 " 3 1 (fg)" t L/Ly D JD a f!a g 2 1 3 H (13) The new formulas of element displacement lead to two important properties of isoparametric element. Property I: ¹he work done by higher-order stresses due to constant strains is zero, i.e. Pv BT) DB# dv"t PA BT) dA · DB#"0 K (e)"K(e)T" )# #) e (14) e where D is the elasticity matrix and t is the element thickness. From the equation PA G H P P G H PA PA C D e GH t 1 1 b g!b f 0 1 dA" 3 1 df dg" t a f!a g 0 ~1 ~1 1 3 2 we can infer that e BT dA"C i) i e t 0 t 1 2 dA"0 0 t t 2 1 so we have the conclusion of Property I. Property II. ¹o the arbitrary nodel displacements q"q which correspond to constant strains, # the equation below is satisfied: 4 + C q "0 (15) i i# i/1 Since the compatible element Q4 can exactly imitate constant stress conditions, we must have, when q"q , # e"Bq "(B #B )q "e "B q # # ) # # # # and 4 4 B q " + B q "W(f, g) + C q "0 ) # i) i# i i# i/1 i/1 This equation must be satisfied pointwise in the element, so it results in equation (15). It will be discussed below that property II will have important applications in the research of incompatible elements. Based on Property I, the element stiffness submatrix can be written as K(e)"K(e) #K(e) ij ij# ij) (16a) 19 A NEW FORMULATION OF ISOPARAMETRIC FINITE ELEMENTS where PA BTi# DBj# dA"4AtBTi# DBj# K(e) "t ij) PA BTi) DBj) dA"Ci Cj W K(e) "t ij# (16b) e (16c) e in which W"t PA WT DW dA (17) e Because the matrix W is independent of the indices i and j, only a matrix w of dimension 2]2 is needed to be calculated for the numerical integration of strain energy of element Q4 when new formulas are used. NEW METHOD FOR THE CONSTRUCTION OF INCOMPATIBLE ELEMENT MODELS The formulas of incompatible element can be established from the potential energy functional CP P P D 1 n"+ n(e)"+ (Lu)T D(Lu) dv! FT u dv! p6 T u ds (18) q q 2 ve ve spe (e) (e) where F are the body forces, p6 are boundary traction, Lu"e are element strains, D is elasticity matrix, and the displacements u are separated into two components u"u #u "Nq#N k (19a) q j j with u the compatible displacements as shown in equation (10), and u the added incompatible q displacements which are expressed in terms of internal parameters k. The strains e are also divided into two parts correspondingly: e"e #e "Lu #Lu "Bq#B k (19b) q j q j j Based on the stationary condition of total potential energy functional, the internal parameters k can be eliminated on the element level: 4 k"!R~1 G q"! + R~1 G q i i i/1 (20) where Pv BTj D Bj dv, R" e Pv BTi D Bj dv G" i (21) e Substituting equation (20) into (19), we get the modified displacements and strains expressed by means of nodal displacements q: 4 u*"(N!N R~1G)q"N*q" + N* q j i i i/1 4 e*"(B!B R~1G)q"B*q" + B* q j i i i/1 (22a) (22b) 20 Z.-P. JIAO, T. H. H. PIAN AND S. YONG where N*"N #(N !N R~1G )"N #N* (23a) i i# i) j i i# i) B*"B #(B !B R~1G )"B #B* (23b) i i# i) j i i# i) In order to guarantee the incompatible element to pass the constant strain patch test, the following constraint condition must be satisfied: Pv Bj dv"0 (24) e Substituting (24) into (21) gives Pv BTj D(Bi##Bi) ) dv"Pv BTj DBi) dv"Ci Pv BTj DW dv"Ci G0 G" i e e (25) e Thus, the modified higher-order strain submatrix may be written as B* "C (W!B R~1G )"C W* i) i j 0 i Using Property I, the constraint equation (24) may be expressed equally by equation Pv W* dv"0 (26) (27) e Notice that the work done by external forces due to the incompatible displacements has been neglected in equation (18). It is equivalent to Pv F3uj dv#P4p p6 T uj ds"0 e (28) e This convention in formulating the incompatible elements indicates that only the element strain energy (not the external force potential energy) could be changed when u is introduced. j Comparing equation (26) with equation (12), we find that the entire difference between compatible element and incompatible element is only the change of ( to (*. So we describe (* as the characteristic matrix of incompatible element. Because the properties of incompatible elements entirely depend on the matrix (*, it is a natural idea to construct (* directly instead of computing matrices R, R~1 and G . Essentially, 0 the incompatible element formed in this way is a kind of assumed (higher-order) strain element. Choosing a matrix (* arbitrarily based on equation (27), we have K(e) "t ij#) PA BTi# DB*j) dA"tBTi# DCj PA W* dA"0 e (29) e So the element stiffness matrix can be calculated according to equation (16), in which W"t PA (*T D(* dA (30) e Using Property II, we will find, when q"q , # 4 4 B* q " + B* q "W* + C q "W* · 0"0 ) # i) i# i i# i/1 i/1 (31) 21 A NEW FORMULATION OF ISOPARAMETRIC FINITE ELEMENTS Thus, from arbitrary matrix (* satisfying equation. (27), we can correspondingly obtain an incompatible element which can pass the patch test. It will be very convenient, by using the virtual parameter method,5 to modify the matrix (* and make it to satisfy equation (27). On the other hand, eventhough matrix (* does not satisfy equation (27), we can also obtain the incompatible element which can pass the patch test if the coupling term K is neglected in the element stiffness #) matrix. As a result of equation (28), the external force potential energy will not change when incompatible displacements u are added. Therefore, the substance to condense k is only based on j the stationary condition of element higher-order strain energy. This understanding provides a helpful theoretical guidance to construct matrix (* directly. INCOMPATIBLE ELEMENTS EQUIVALENT TO HYBRID STRESS ELEMENTS Hybrid stress element can be derived from the functional !"+ !(e)"+ (e) (e) GP C ve D CP 1 ! rT Sr#rT (Lu ) dv! q 2 ve P4 p6 uq dsDH FT u dv# q T (32) p where S"D~1 is elastic compliance matrix, u are compatible displacements of element shown in q equation (10), r are stresses of element and r"Pb"[I P ]b"b #P b ) # ) ) (33) in which P and I is the higher-order stress model and a 3]3 unit matrix for plane problems, ) respectively. After the condensation of the parameters b, the element stiffness submatrix is K(e)"MT H~1M i j ij (34) in which Pv PT SP dv" S Pv dv S Pv P) dv " H11 HT C 12 sym PT SP dv Pv ) ) e H" e e H 12 H 22 D (35) e Pv PT Bi dv" M" i e Pv Bi dv Pv e e PT B dv ) i C D " M i% M i) (36) To simplify the formulating procedure (not an essential condition), we use the model suggested in Reference 6, wherein the results are identical to that of the original formulation. We modify the higher-order model P by making it satisfy ) Pv P) dv"0 e (37) 22 Z.-P. JIAO, T. H. H. PIAN AND S. YONG This equation leads to H "0 and 12 Pv (Bi##Bi)) dv"Pv Bi# dv#tCi PA W dA"4AtBi# M " PT (B #B ) dv" PT B dv"C t PT W dA"C Q i) i i) P ) i# Pv ) i) i Pv ) v M " i# e e e (38a) e e (38b) e So the element stiffness submatrix may be written as K(e)"[MT MT ] i# i) ij C 0 H~1 11 0 H~1 22 DC D M j# "4At BT DB #C C QT H~1 Q 22 i# j# i j M i) (39) in which we have used 1 H~1"(4At S)~1" D 11 4At Now, let us return to the formulation of incompatible element described in the last section. We know that the added incompatible displacements u ought to be independent of the higher-order j displacements u of the compatible element, otherwise the components of u would vanish q) q) during the condensation of k. So the components related to ( must be deduced from matrix B j when B is constructed, i.e. the matrix B should be formed as j j B "B !W/ (40) j j1 where / is a constant matrix to be determined. Substituting equation (40) into equation (23b), leads to B* "C [W!(B !W/) R~1 G ]"C W* (41) i) i j1 0 i Owing to the factor 1/DJD in matrix W, the form of matrix W* that is composed of W and B will be very complicated, and it is difficult to see the ‘true face’ of matrix W* and to obtain more rational higher-order strain matrix without the factor 1/DJD. In order to gain the matrix W* whose form can be chosen arbitrarily, one approach is to choose matrix W as W#W/R~1G "0 0 (42) Then equation (41) is simplified as B* "!C B R~1G i) i j1 0 From equation (42), we have the following relation (see the appendix): R~1G "!R~1 Q 1 1 0 (43) (44) in which Pv BTj1 DBj1 dv, R " 1 e Pv BTj1 DW dv Q " 1 (45) e Substituting equation (44) into equation (43) yields the modified higher-order strain submatrix which is only dependent on matrix B : j1 * (46) B "C B R~1 Q "C W* 1 i i) i j1 1 23 A NEW FORMULATION OF ISOPARAMETRIC FINITE ELEMENTS By using equation (46), it becomes very simple to construct the desirable higher-order strain matrix. If we choose B "SP j1 ) (47) then we will have Pv BTj1 DBj1 dv"Pv PT) SP) dv"H22 Q " BT DW dv" PT W dv"Q 1 P j1 Pv ) v R " 1 e (48a) e e (48b) e and Pv W* TDW* dv"QT1 R~11 Pv BTj1 DBj1 dv R~11 Q1"QT H~122 Q (49) W" e e Substituting equation (49) into equation (16) and comparing with the element stiffness matrix of hybrid stress model shown in equation (39), we could find immediately that the element stiffness matrix of incompatible element whose higher-order strain matrix is constructed in equations (46) and (47) is entirely equivalent to one kind of hybrid stress element. Numerical solutions show that the incompatible element here and the hybrid stress element are entirely equivalent in different mesh divisions and different higher-order stress models P . It is not difficult to deduce the element ) higher-order displacement field from equation (46): The matrix B j1 4 u*" + C N R~1Q q ) i j1 1 1 i i/1 is expressed by Cartesian co-ordinate x—y in the following form: (50) a x#b y a x#b y 11 11 12 12 B "SP " a x#b y a x#b y j1 ) 21 21 22 22 a x#b y a x#b y 31 31 32 32 (51) We can easily obtain C D b xy#1 a x2#1 (b !a )y2 b xy#1 a x2#1 (b !a )y2 2 12 2 11 2 31 2 32 12 21 22 (52) N " 11 j1 a xy#1 b y2#1 (a !b )x2 a xy#1 b y2#1 (a !b )x2 2 31 2 32 2 21 2 22 11 12 21 22 In order to satisfy equations (27) and (37), we only have to change x to x* and y to y* in the equations above. x*"x!a , y*"y!a 1 2 (53) where 1 a " (a B#a C), 1 3A 1 3 1 a " (b B#b C) 2 3A 1 3 (54) It should be pointed out that the matrix N corresponding to equation (47) cannot be found in j1 general cases when P is expressed by means of isoparametric coordinate f!g. ) 24 Z.-P. JIAO, T. H. H. PIAN AND S. YONG THE HYBRID STRESS ELEMENTS EQUIVALENT TO INCOMPATIBLE ELEMENTS The element displacement field can be expressed as the combination of compatible part u and q incompatible part u : j u*"u #u (55) q j After the condensation of k, element displacements and strains change their forms, respectively, 4 4 4 u*" + N* q " + (N #N* )q " + (N #C N* )q (56a) i i i# i) i i# i 0 i i/1 i/1 i/1 4 4 (56b) e*"¸u*" + B* q " + (B #C W*)q i i i# i i i/1 i/1 To develop hybrid stress element equivalent to incompatible element, we replace the u in (32): q 1 !"+ !(e)"+ ! rTSr#rT(Lu*) dv! FT u dv# p6 T u ds (57) q q 2 ve ve spe (e) (e) under the stationary condition of the above functional, we have GP C d!"+ (e) GP C ve D CP P DH Psp (T!p6 )T duq dsH (Lu*!Sr)T dr!(LTr#F)Tdu ] dv# q P e P #+ [T (a)#T(b)]T du ds#+ rT (Ldu ) dv"0 (58) q j ab (e) ve where T"nr, n is the direction cosine of outward normal on element sides, and (a), (b) represent the neighbouring elements. If the equation P + rT (Ldu ) dv"0 (59) j (e) ve is satisfied the stationary condition (58) will yield all the equilibrium equations except boundary displacement conditions. Considering the limit situation refining finite element meshes, rPr (r # # is constant stress), equation (59) could be changed into stronger constraint condition of single element: Pv Lduj dv"0 (60) e Noticing Pv Luq) dv,0 e we can rewrite equation (60) as 4 or Pv Ldu*) dv"Pv W* dv i/1+ Ci dqi"0 e e Pv W* dv"0 e (61) 25 A NEW FORMULATION OF ISOPARAMETRIC FINITE ELEMENTS Therefore, under the condition (61), we can construct a convergent hybrid stress element model from the functional (57). In general cases, the incompatible elements deduced from (56) may not satisfy equation (61) such as the incompatible element by Wilson. The element stiffness submatrix of these kinds of incompatible elements is K(e)"4AtBT DB #C C W#C ETDB #C BT DE i# j# i j i j# j i# ij (62) in which Pv W* TDW* dv, Pv W* dv W" E" e (63) e Taking the higher-order stress model as P "DW* ) the element stiffness submatrix of hybrid stress element may be written as (64) K(e)"MTH~1M i j ij (65) in which Pv Pv W* dv 4At · S H" e PT SP dv" e sym. Pv W*T DW* dv " C 4At · S ET E W D (66a) e and i# i i# i dv" Pv PTB*i dv"Pv C W*TD(B #C W*) D ETDB #C W D C i# i i# i M" i B #C W* e 4AtB #C E (66b) e Let C A" D A A 11 12 "H~1 AT A 12 22 (67) then we have the following relations. 4AtSA #EAT "I 12 3 11 ETA #WA "I 12 22 2 (68) 4AtSA #EA "0 12 22 ETA #WAT "0 11 12 where I and I are 2]2 and 3]3 unit matrix, respectively. 2 3 Substituting equations (66) and (67) into (65) and using equation (68) to simplify the equation, we have K(e)"MT AM "4AtBT DB #C C W#C ETDB #C BT DE (69) i j i# j# i j i j# j i# ij Comparing (69) with (62), we find at once that the stiffness matrix of incompatible elements is just the same as the one of hybrid stress element. In fact, if and only if the displacement field of incompatible elements may be expressed as shown in equation (56), the equivalent hybrid stress 26 Z.-P. JIAO, T. H. H. PIAN AND S. YONG element can certainly be constructed based on the functional (57) and with the higher-order stress model in equation (64). It is obvious that the hybrid stress element constructed by equation (69) can pass the patch test if the matrix W* in equation (64) satisfies constraint condition (27). If P "DW is chosen, as ) a special example, we will obtain the hybrid stress element that is equivalent to compatible element Q4. DISCUSSION From the above-mentioned research works, we have observed that (i) Because all incompatible elements have their characteristic matrix W* the equivalent hybrid stress elements can always be constructed. (ii) If the higher-order stress model P of hybrid stress element is expressed by means of ) isoparametric co-ordinates, even though the incompatible elements that have the same stiffness matrix can be constructed, usually the adding internal displacements u cannot be j found. The phenomenon above seems to show that the j-type incompatible element is totally contained in the frame of hybrid stress element. CONCLUSION 1. A new formulation of isoparametric element, taking element Q4 for example, has been presented which can reduce the CPU time required for calculating stiffness matrix, and reveal two important properties. 2. Because it is very convenient to construct incompatible element according to the characteristic matrix W* which may be arbitrarily designed, the new method will be a great help to the study of incompatible element. 3. The relationship between hybrid stress element and incompatible element has been established, which provides a helpful theoretical basis for the research of hybrid stress element and incompatible element. 4. All of the discussions in this paper are based on the plane quadrilateral isoparametric element, the further studies are expected to extend these conclusions to the other cases. APPENDIX Equation (42) is equivalent to /R~1G "!I 0 (70) /R~1G /"!/ 0 (71) where I is 2]2 unit matrix. Postmultiplying equation (70) by /, The above equation will hold if the following condition is satisfied: G /"!R 0 (72) A NEW FORMULATION OF ISOPARAMETRIC FINITE ELEMENTS 27 Notice that Pv BTj DBjdv"Pv (BTj1!/TWT)D(Bj1!W/) dv"R1!/T QT1!Q1 /#/TW0 / G " BT DW dv" (BT !/TWT) DW dv"Q !/TW 1 0 0 P j Pv j1 v R" e e e e (73) (74) in which R and Q are expressed in equation (45) and 1 1 Pv WTDW dv W " 0 (75) e Substituting (73) and (74) into (72) gives (76) R "/TQT"RT"Q / 1 1 1 1 When equation (72) is premultiplied by R~1 and equation (76) is premultiplied by R~1, we have 1 R~1 G /"!I (77) 0 and R~1Q /"I 1 Comparing equation (77) with equation (78), we obtain the relation R~1G "!R~1 Q 1 1 0 (78) (79) REFERENCES 1. T. H. H. Pian, ‘On the equivalence of non-conforming element and hybrid stress element’, Appl. Math. Mech. (English edition), 3, 773—776 (1982). 2. E. L. Wilson, R. L. Taylor, W. P. Doherty and J. Ghaboussi, ‘Incompatible displacement models’, in S. J. Fenves et al. (eds.), Numerical and Computer Methods in Structural Mechanics, Academic Press, New York, 1973, pp. 43—57. 3. T. H. H. Pian, ‘Derivation of element stiffness matrices by assumed stress distributions’, AIAA J., 2, 1333—1336 (1964). 4. O. C. Zienkiewicz and R. L. Taylor, ¹he Finite Element Methods, 4th edn. McGraw-Hill, London, 1989. 5. C. C. Wu, M. G. Huang and T. H. H. Pian, ‘Consistency condition and convergence criteria of incompatible elements: General formulation of incompatible functions and its application’, Comput. Struct., 27, 639—664 (1987). 6. Z. P. Jiao, ‘A method for the simplification of the matrix H of hybrid stress element’, Comput. Struct. Mech. Appl., 8, 214—216 (1991) (in chinese). .