INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 4313—4339 (1997) CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS OF HYBRID FINITE ELEMENTS W. FENG1, S. V. HOA1* AND Q. HUANG2 1Department of Mechanical Engineering, Concordia ºniversity, 1455 de Maisonneuve Blvd. ¼., H549, Montreal, Que., Canada H3G 1M8 2 Shanghai ºniversity, Shanghai, China ABSTRACT A classification method is presented to classify stress modes in assumed stress fields of hybrid finite element based on the eigenvalue examination and the concept of natural deformation modes. It is assumed that there only exist m ("n!r) natural deformation modes in a hybrid finite element which has n degrees of freedom and r rigid-body modes. For a hybrid element, stress modes in various assumed stress fields proposed by different researchers can be classified into m stress mode groups corresponding to m natural deformation modes and a zero-energy stress mode group corresponding to rigid-body modes by the m natural deformation modes. It is proved that if the flexibility matrix [H] is a diagonal matrix, the classification of stress modes is unique. Each stress mode group, except the zero-energy stress mode group, contains many stress modes that are interchangeable in an assumed stress field and do not cause any kinematic deformation modes in the element. A necessary and sufficient condition for avoiding kinematic deformation modes in a hybrid element is also presented. By means of the m classified stress mode groups and the necessary and sufficient condition, assumed stress fields with the minimum number of stress modes can be constructed and the resulting elements are free from kinematic deformation modes. Moreover, an assumed stress field can be constructed according to the problem to be solved. As examples, 2-D, 4-node plane element and 3-D, 8-node solid element are discussed. ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) No. of Figures: 0. No. of Tables: 2. No. of References: 23. KEY WORDS: finite element; stress modes; classification 1. DEFINITION In order to clarify the presentation, it is useful to give the definitions of the stress field, stress mode, stress matrix, natural deformation mode, rigid-body mode, zero-energy stress mode, kinematic deformation mode and stress mode group firstly. Stress field: In the hybrid formulation, a stress field is assumed independently from beginning, MpN"[P]MbN (1) * Correspondence to: S. V. Hoa, Department of Mechanical Engineering, Concordia University, 1455 de Maisonneuve Blvd. W., H 549 Montreal, Quebec H3G 1M8, Canada CCC 0029—5981/97/234313—27$17.50 ( 1997 John Wiley & Sons, Ltd. Received 27 December 1995 Revised 8 April 1997 4314 W. FENG, S. V. HOA AND Q. HUANG For example, for a 2-D, 4-node plane element, one of the assumed stress fields is b 1 1 0 0 [x 1 p x b p " 1 [1 0 [y 0 2 y F p 0 0 1 x y xy b 5 (2a) and one of the assumed stress fields for a 3-D, 8-node solid element is p x p y p z " p xy p yz p zx 1 1 [1 0 0 0 z 0 y z 1 [1 [1 0 0 0 z x 0 [z 0 y 0 0 0 yz 0 0 x 0 0 0 0 0 zx 0 [x [y 0 0 0 0 0 xy z z 0 0 0 1 0 2 0 0 0 0 x y 0 0 0 0 1 0 0 0 0 0 0 0 0 z 0 0 0 0 1 0 0 0 0 0 0 0 x [x [x 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 y [2y 0 0 0 ] b 1 b 2 F b 18 0 (2b) They can be expressed in the form MpN"[Mp N Mp N . . . Mp N] 1 2 m "[P ]MbN b 1 b 2 F b m (2c) In which, the parameters b are the corresponding stress parameters and [P] is the stress matrix. i Stress modes: Stress modes are vectors which are functions of the co-ordinates. For example, vectors Mp N in equation (2c) are stress modes. i Stress matrix: An assumed stress matrix consists of several stress modes. In equation (2c), the matrix [P] is a stress matrix. Natural deformation modes: Deformation modes in an element that are independent from each other. In this work, the eigenvectors of an element stiffness matrix are regarded as the natural deformation modes of the element. Rigid-body modes: Rigid-body modes are displacement modes in an element that do not produce deformation energy. Zero-energy stress modes: Stress modes in the stress matrix of a hybrid element that do not produce deformation energy. The eigenvalues of the element stiffness matrix corresponding to these stress modes equal zero. These zero-energy stress modes correspond to rigid-body modes. Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4315 Kinematic deformation modes: Deformation modes in an element corresponding to spurious zero stiffness may be caused by unsuitable numerical integration technique or unsuitable assumed stress fields. In this work, it is assumed that suitable numerical integration technique is used. Even though zero-energy stress modes may not appear in the stress matrix, spurious zero stiffness can occur. This happens when more than one stress mode is picked from one eligible stress mode group. This zero stiffness mode is called kinematic deformation mode in this work. The exception for this case is when the stress modes are interchangeable between different stress mode groups. This exception is discussed in Section 6. Stress mode group: A stress mode group contains many stress modes that are interchangeable in the stress matrix [P] and do not cause kinematic deformation modes. The m stress mode groups correspond to m natural deformation modes and the zero-energy stress mode group corresponds to rigid-body modes. If an element has n degrees of freedom, its nodal displacement vector must consist of n components. The displacement distribution in the element can be represented by m ("n!r) natural deformation modes and r rigid-body modes. In this work, it will be shown that the m natural deformation modes can classify stress modes in various assumed stress fields for a hybrid element. All existing stress modes are classified into m ("n!r) stress mode groups corresponding to m natural deformation modes and a zero-energy stress mode group corresponding to rigid-body modes in the element. For a hybrid element, the displacement field is assumed, MuN"[N]MdN (3) where [N] is the shape functions and MdN is a nodal displacement vector. The hybrid element stiffness matrix [K] in nodal displacement space can be obtained by using Hellinger—Reissner variational principle, [K]"[G]T[H]~1[G] (4) [H]MbN"[G]MdN (5) PV[P]T[S][P] d» [G]" [P]T[B] d» P (6) and where [H]" V in which, [S] is the material properties matrix relating stresses to strains, [B] is the geometric matrix relating strains to displacements, [H] is the flexibility matrix, and [G] is the leverage matrix. 2. INTRODUCTION Since the time the hybrid finite element was developed by Pian1 in 1964, various hybrid finite elements have been proposed.2,3 In the development of hybrid finite elements, a major problem has obstructed the application of this method. That is the lack of a rational way for deriving ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4316 W. FENG, S. V. HOA AND Q. HUANG the optimal assumed stress modes. Recently, several partial hybrid finite elements have been proposed4~6 and the problem concerning the method for deriving assumed stress field stands out. If an assumed stress field does not contain enough stress modes, the rank of the element stiffness matrix will be less than the total degrees of deformation freedom and the numerical solution of the finite element model will not be stable. It is possible to suppress kinematic deformation modes by adding stress modes of higher-order term, but this cannot guarantee that all kinematic deformation modes are suppressed. Moreover, each extra term will add more stiffness7 and overuse of stress modes will cost more computational time because the calculation of element stiffness matrix requires inversion of the flexibility matrix. Therefore, m ("n!r) least-order stress modes are considered to be best and are optimal with respect to computer resources.7,8 In this work, the method to derive this kind of assumed stress field will be given using the classified stress mode groups and the necessary and sufficient condition for avoiding kinematic deformation modes in a hybrid element. Some mathematical basis for the stability of the numerical solution of the finite element model has been established and a number of approaches for obtaining the optimal stress modes have been proposed. Fraeijs de Veubeke,9 Pian and Tong10 presented a necessary condition to avoid kinematic deformation modes, m@*n!r (7) in which m is the total number of stress modes in an assumed stress field, n is the total number of nodal displacements, and r is the number of rigid-body modes in an element. Brezzi,11 Babuska et al.12 presented necessary and sufficient conditions for stability and convergence of a hybrid element. However, it is difficult to use these conditions because of the abstract concept and the complex analysis involved. Ahmad and Irons13 suggested use of an eigenvalue technique to assess a hybrid element and to determine the kinematic modes. Huang14 introduced the concept of natural deformation mode and natural stress mode and developed a modal analysis technique to find natural stress modes for hybrid elements. He gave a set of uncoupled stress modes yet without zero-energy stress modes by means of modal analysis of 3-D, 8-node hybrid element. Using group theory, Punch and Atluri8,15 tried to solve the zero-energy mode problem. They gave two assumed stress fields for 2-D, 4-node hybrid element, eight assumed stress fields for 3-D, 8-node hybrid element and 384 assumed stress fields for 3-D, 20-node hybrid element. Pian and Chen16 used the product MpNTMeN, the deformation energy due to the assumed stresses and displacements, to determine the necessary assumed stress modes. The assumed stress fields for 2-D, 4-node plane element and 3-D, 8-node solid element were provided. Pian and Tong17 introduced the internal displacement parameters to relax the stress equilibrium condition and used isoparametric interpolation to construct hybrid element. Pian and Wu18 introduced incompatible displacements to maintain completeness. The initial choice of stress terms are unconstrained and complete polynomials. The additional displacements are used as Lagrange multipliers to enforce the stress equilibrium constraint. Chen et al.19,20 constrained the stress by setting the inner product of the non-constant stress modes with the deformation derived from the incompatible displacement to zero. Sze21 used orthogonal lower- and higher-order stress modes to construct hybrid element. It allows the partition of the element stiffness matrix into a lower- and a higher-order stiffness matrix. When the lower-order stiffness turns out to be identical to the sub-integrated element, the higher-order stiffness matrix plays the role of stabilization matrix. Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4317 For the hybrid element, an assumed stress field can be constructed by various approaches.22 One type of hybrid element usually has many different assumed stress fields. However, the relationship among them has not been investigated although different approaches to assume stress field of hybrid elements have been proposed. In this work, a classification method is presented to classify stress modes in various assumed stress fields corresponding to the 2-D, 4-node plane element and the 3-D, 8-node solid element. It will reveal the relationship between different assumed stress fields and answer the question why there exist kinematic deformation modes in a hybrid element even when the criterion (m@'n!r) is satisfied. Also, a systematic procedure is given for constructing assumed stress fields of hybrid elements that contain minimum number of stress modes and are free from kinematic deformation modes. 3. CLASSIFICATION METHOD The hybrid formulation based on the Hellinger—Reissner principle relaxes the stress equilibrium condition. The stress field would satisfy the equilibrium equations only in a variational sense. Therefore, the stress field may be described in the isoparametric co-ordinate system of the element, which would make the element less sensitive to mesh distortion. The use of isoparametric co-ordinates also preserves the invariance of the resulting stiffness matrix. There are many assumed stress fields for 2-D, 4-node plane element and 3-D, 8-node solid element. For example, Pian16 proposed an assumed stress field for 2-D, 4-node plane element and another for 3-D, 8-node solid element. Punch and Atluri8,15 gave two assumed stress fields for 2-D, 4-node plane element, and eight assumed stress fields for 3-D, 8-node solid element. Huang14 presented an assumed stress field for 3-D, 8-node solid element. Although each of these assumed stress fields may contain the same number of stress modes, the stress modes in these fields are different. In order to study the relationship between different stress modes, the concept of natural deformation modes is used. A finite element has a finite number of degrees of freedom. For instance, a 2-D, 4-node displacement element has (n") 8 degrees of freedom, and a 3-D, 8-node displacement element has (n") 24 degrees of freedom. In an element, there exist (n!r) natural deformation modes and r rigid-body modes. The displacement distribution in the element can be represented by them. The equilibrium equation of a displacement element is [K]MdN"MFN (8) If the nodal force vector is proportional to the nodal displacement vector, the equilibrium equation becomes a eigenvalue equation. It can be expressed as follows: ( [K]!j[I])MdN"0 (9) [K] is an n]n element stiffness matrix. This equation will give (n!r) non-zero eigenvalues and r zero eigenvalues, and (n!r) eigenvectors corresponding to the (n!r) non-zero eigenvalues. The (n!r) eigenvectors MdN (i"1, 2, 3, . . . , m) depend only on the geometry and material i properties of the element, and they are unique. If vectors Md N (i"1, 2, . . . , m) are the eigenveci tors of the stiffness matrix [K], they must satisfy the condition: Md NTMd N"0, iOj i j (10) Md NTMd N"1, i"j i j ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4318 W. FENG, S. V. HOA AND Q. HUANG In this work, these eigenvectors are considered to be the natural deformation modes of the element. They can represent the deformation in the element. The strain modes derived from the assumed displacement directly also can represent the deformation in the element, but they cannot classify stress modes in hybrid element. Moreover, the strain modes cannot provide a general way to derive necessary stress modes for all types of hybrid elements. For example, they can be used to derive a set of stress modes for conventional hybrid elements, but they cannot be used to derive stress modes for partial hybrid element. In the hybrid element, the eigenvectors and eigenvalues of the stiffness matrix will be sensitive to the assumed stress modes. The eigenvalue examination will give r zero eigenvalues corresponding to rigid-body mode and m ("n!r) non-zero eigenvalues corresponding to natural deformation modes if the assumed stress field is suitable. For a hybrid element, if various stress modes can be classified, at least m stress mode groups must exist because the stiffness matrix of hybrid element must have m non-zero eigenvalues, except zero-energy stress mode group. Otherwise, the hybrid element will contain kinematic deformation modes. On the other hand, no matter how many stress modes there are in a stress matrix [P], the maximum number of non-zero eigenvalues of an element stiffness matrix is always equal to or less than m. Therefore, the number of stress mode groups is also equal to or less than m. Thus, it can be considered that there exist and only exist m stress mode groups except zero-energy modes. All stress modes in various assumed stress matrices can be classified into the m stress mode groups corresponding to m natural deformation modes and the zero-energy stress mode group corresponding to rigid-body modes. Postulate: There exist and only exist m ("n!r) natural deformation modes in a hybrid element. All stress modes in assumed stress matrices can be classified into m stress mode groups corresponding to m natural deformation modes and a zero energy mode group corresponding to rigid-body modes of the element which has n degrees of freedom and r rigid-body modes. Based on this postulate, it can be considered that an assumed stress field can be represented by stress modes in the m stress mode groups related to m natural deformation modes, except zero-energy stress modes. This can be expressed as follows: b1 b n1 n2 nm m MpN"[P]MbN" + Mp N, + Mp N, . . . , + Mp N 2 " + [P ]Mb N l1 l2 lm i i F l1/1 l2/1 lm/1 i/1 b C D (11) m where n is the number of stress modes in the ith stress mode group, and [P ] and Mb N i i i (i"1, 2, . . . , m) are the stress matrices and stress-coefficient vectors related to the ith stress mode group which corresponds to the ith natural deformation mode. They are C ni [P ]" M0N . . . M0N, + Mp N, M0N . . . M0N li i li/1 D (12) and Mb N"[0 . . . 0 b 0 . . . 0]T (13) i i The summation sign in the matrix [P] is not necessary to be carried out. One can have only one stress mode or many. Many stress modes will increase the size of stress matrix [P], and this is not desirable. The stress mode which belongs to the ith stress mode group can be expressed in the Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4319 form Mp N"[P]Mb N i i (14) Therefore, the vector Mb N represents the ith stress mode group which corresponds to the natural i deformation mode Md N (i"1, 2, . . . , m). Using equation (5), we have i [H]Mb N"[G]Md N i i (15) If the stress matrix [P] does not contain any stress mode which belongs to the ith stress mode group, the value of b in the vector Mb N should be equal to zero. Then, we can add a new stress i i mode into the stress matrix [P]. The new stress mode will be classified by m natural deformation modes. Corresponding to the ith natural deformation mode Md N, the condition to check whether i the new stress mode belongs to the ith stress mode group can be expressed in the form, b "0 if new stress mode does not belong to ith stress mode group i b O0 if new stress mode belongs to ith stress mode group i Using equations (9) (10) and (15), the eigenvalues are obtained as follows: j "Md NT [K]Md N"Mb NT[H]Mb N i i i i i (16) Because all of the diagonal elements in the flexibility matrix [H] may not be equal to zero, the classification condition above becomes j "0 if new stress mode does not belong to ith stress mode group i j O0 if new stress mode belongs to ith stress mode group i Before classifying stress modes, one can find a number of initial stress modes since there are many ways to derive assumed stress matrices for a hybrid element. For example, Pian and Chen16 used the product MpNTMeN to determine the necessary assumed stress modes, and gave an assumed stress matrix for 2-D, 4-node plane element and another for 3-D, 8-node solid element. Punch and Atluri8,15 used group theory to obtain two assumed stress matrices for 2-D, 4-node hybrid element, eight assumed stress matrices for 3-D, 8-node hybrid element, and 384 assumed stress matrices for 3-D, 20-node hybrid element. One can also derive an assumed stress matrix using iso-function method.23 In the stress matrix derived using iso-function method, the number of the stress modes is larger than m ("n!r). In order to present a systematic procedure for classifying stress modes and constructing assumed stress fields, the iso-function method23 is used to derive initial stress modes to be classified in this work. This is because the hybrid element constructed by the iso-function stress matrix has the same eigenvalues and eigenvectors as conventional displacement element. Also, the method using iso-function is straightforward and can be followed easily. After obtaining initial stress modes, one can use eigenvalue examination to find m representative stress modes that represent m stress mode groups corresponding to m natural deformation modes. The stress matrix consisted of the m representative stress modes is a optimal stress matrix. Then, all existing stress modes can be classified into m#1 stress mode groups. Its detail is presented as follows. ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4320 W. FENG, S. V. HOA AND Q. HUANG 3.1. Determination of optimal stress matrix from the stress matrix derived by iso-function method Step 1. Derive an initial stress matrix [P] by iso-function method. The number of initial *40 stress modes in the matrix is always larger than m ("n!r). In order to select m necessary stress modes, these initial stress modes have to be classified into (n!r) stress mode groups. Step 2. Select stress modes in the order from low-order term to high-order term. Now select a stress mode from the existing stress matrix [P] and form an assumed stress matrix [P ]. The *40 1 element stiffness matrix [K] corresponding to stress matrix [P ] can be obtained by using 1 equations (4) and (6). If the eigenvalue examination gives a non-zero eigenvalue, the stress mode is a non-zero-energy stress mode; otherwise, it is a zero-energy stress mode. Repeating the eigenvalue examination to check whether a stress mode is a zero- energy stress mode for all stress modes in the existing stress matrix [P] . *40 Take all zero-energy stress modes out and keep non-zero-energy stress modes in the matrix [P] . All zero-energy stress modes form a zero-energy stress mode group. *40 Step 3. Take a non-zero-energy stress mode from the existing stress matrix [P] and form an *40 assumed stress matrix [P ]. The stress mode Mp N is the representative stress mode which 1 1 represents group 1 of stress modes. Step 4. Add another stress mode selected in the existing stress matrix [P] into the assumed *40 stress matrix [P ] and form a new stress matrix [P ], 1 2 [P ]"[Mp N Mp N] 2 1 2 (17) Step 5. The eigenvalue examination gives the eigenvalues of the stiffness matrix. If there is only one non-zero eigenvalue, continue to step 6. If there are two non-zero eigenvalues, go to step 7. Step 6. In this case, the added stress mode belongs to group 1 of stress modes. Take the second stress mode out and put it in group 1 of stress modes. Then, go back to step 4. Step 7. The two stress modes belong to two different groups of stress modes. The second stress mode Mp N is the representative stress mode which represents group 2 of stress modes. 2 Step 8. Add another stress mode selected from the matrix [P] into the assumed stress matrix *40 [P ] and form a new stress matrix [P ], 2 3 [P ] " [Mp N Mp N Mp N] 3 1 2 3 (18) Step 9. The element stiffness matrix [K] and its eigenvalues are calculated. If there are only two non-zero eigenvalues, continue to step 10. If there are three non-zero eigenvalues, go to step 11. Step 10. In this case, the new stress mode Mp N belongs to one of the two stress mode groups. 3 Construct the matrices [P@ ] and [P@@ ] as follows: 2 2 [P@ ]"[Mp N Mp N] or [PA ]"[Mp N Mp N] 2 2 3 2 1 3 (19) If the stiffness matrix corresponding to the stress matrix [P@ ] has two non-zero eigenvalues, the 2 stress mode Mp N belongs to the group 2 of stress modes. Otherwise, the stress mode Mp N belongs 3 3 to the group 1 of stress modes. Put the stress mode Mp N into the corresponding stress mode 3 group, and go back to step 8. Step 11. In this case, the three stress modes belong to three different stress mode groups. The added stress mode Mp N is the representative stress mode which represents group 3 of stress 3 modes. Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4321 Step 12. Add one more stress mode selected from the matrix [P] into the matrix [P ] and *40 3 form a new stress matrix [P ], 4 [P ]"[Mp N Mp N Mp N Mp N] 4 1 2 3 4 (20) and so on. Repeating the same process until m representative stress modes that represent m stress mode groups are obtained. The m("n!r) representative stress modes correspond to m natural deformation modes and form a optimal stress matrix [P] from the existing stress matrix [P] . 015 *40 Classification of other stress modes Step 13. After m representative stress modes are obtained, other initial stress modes that remain in the existing stress matrix [P] can be classified into the m stress mode groups. Many *40 other stress modes derived by different methods also can be classified into the m stress mode groups corresponding to m natural deformation modes and the zero-energy stress mode group corresponding to rigid-body modes. Based on the optimal stress matrix [P] , any remaining stress mode in [P] can be classified 015 *40 by using it to replace each and every stress mode in the matrix [P] in order. Once the 015 eigenvalue examination results in m non-zero eigenvalues, the representative stress mode which is replaced and the remaining stress mode which replaces the stress mode in [P] belong to the 015 same stress mode group. Put the remaining stress mode into the corresponding stress mode group and recover the optimal stress matrix [P] . Then, classify another remaining stress mode. 015 Step 14. Repeating the same process until all remaining stress modes are classified. Thus, all existing stress modes are classified into m#1 different mode groups. Every stress mode group contains many interchangeable stress modes. For a stress mode derived by other method, if eigenvalue examination always gives m!1 non-zero eigenvalues when this stress mode replaces each and every stress mode in the matrix [P] , this stress mode is a zero- energy 015 stress mode. 4. ILLUSTRATION FOR THE CLASSIFICATION OF STRESS MODES As an illustration for the above procedure, the stress modes presented in References 8, 14—16 and those derived by iso-function method are classified. 4.1. 2-D, 4-node plane hybrid element 4.1.1. Determination of optimal stress matrix from the existing stress matrix derived by isofunction method. The 2-D, 4-node plane element has (n") 8 degrees of freedom and (r") 3 rigid-body modes. So it has (m"n!r") 5 natural deformation modes. Firstly, an assumed stress matrix can be derived from the assumed displacement field of the element by the isofunction method.23 1 0 0 x y 0 [P ]" 0 1 0 0 I 0 0 1 0 ( 1997 John Wiley & Sons, Ltd. 0 0 0 0 x y 0 0 0 0 (21) 0 x y Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4322 W. FENG, S. V. HOA AND Q. HUANG The number of stress modes in the stress matrix is larger than m ("5). The eigenvalue examination indicates that the eigenvalues and eigenvectors of the hybrid element stiffness matrix constructed by the assumed stress matrix [P ] are the same as that of displacement element I stiffness matrix. Here, we take the stress modes in the stress matrix as initial stress modes to be classified. There are nine stress modes in the matrix [P ], I 1 0 0 Mp N" 0 , Mp N" 1 , 1 2 0 0 Mp N" 0 3 1 x Mp N" 0 4 0 , Mp N" 5 y 0 , Mp N" 6 0 0 Mp N" y 7 0 , 0 Mp N" 0 8 x , Mp N" 9 0 x 0 (22) 0 0 y The stress matrix derived by iso-function method contains a few unnecessary stress modes. It will save computation time for calculating element stiffness matrix if the number of the stress modes can be reduced to m ("n!r). In order to do it, the initial stress modes in the existing stress matrix have to be classified into m stress mode groups. First of all, one must find m representative stress modes corresponding to m natural deformation modes. Following steps 2—12 in the procedure of the classification method given in the above section, one can obtain five representative stress modes Mp p p p p N corresponding to (m") of 5 natural deformation 1 2 3 5 6 modes and the zero-energy stress modes Mp N and Mp N corresponding to rigid-body modes. The 4 7 eigenvalues of the stiffness matrix related to Mp p p p p N are not equal to zero, and the 1 2 3 5 6 eigenvalue of stiffness matrix related to Mp N or Mp N is equal to zero. The 5 representative stress 4 7 modes form an optimal stress matrix [P ] from the existing stress matrix [P ], II I 1 0 0 y 0 [P ]"[p p p p p ]" 0 1 0 0 x II 1 2 3 5 6 0 0 1 0 0 (23) The stress matrix is the same as that given by Pian.16 4.1.2. Classification of other stress modes. After obtaining the optimal stress matrix, one can classify stress modes in the existing stress matrix [P ] into (m#1") 6 stress mode groups by I following steps 13 and 14 in the procedure, Tension mode (Group 1): Mp N 1 Tension mode (Group 2): Mp N 2 Shear mode (Group 3): Mp N 3 Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4323 Bending mode (Group 4): Mp N, Mp N 5 8 Bending mode (Group 5): Mp N, Mp N 6 9 Zero-energy stress mode (Group 6): Mp N, Mp N 4 7 The first five stress mode groups correspond to five natural deformation modes and the zeroenergy stress mode group corresponds to rigid-body modes. There are many methods to derive initial stress modes. For example, in the two assumed stress matrices derived by means of group theory8,15 for the same finite element, there are four stress modes that are different from stress modes Mp N—Mp N: 1 9 Mp N" 10 1 1 0 , 1 Mp N" [1 11 0 , 0 Mp N" [y 12 x , !x Mp N" 0 13 y (24) Moreover, one may want to introduce some stress modes of high-order terms into the assumed stress matrix [P] in order to describe special stress distribution in a local region of a structure to be solved. For example, x2 Mp N" 0 14 0 , 0 Mp N" x2 15 0 , 0 Mp N" 0 16 x2 y2 Mp N" 0 17 0 , 0 Mp N" y2 18 0 , 0 Mp N" 0 19 y2 xy Mp N" 0 20 0 , 0 Mp N" xy 21 0 , Mp N" 22 (25) 0 xy 0 According to the steps 13 and 14 in the procedure of classification method, these new stress modes Mp N—Mp N can also be classified into the 6 stress mode groups above, 10 22 Tension mode (Group 1): Mp N, Mp N, Mp N, Mp N 1 10 14 17 Tension mode (Group 2): Mp N, Mp N, Mp N, Mp N 2 11 15 18 Shear mode (Group 3): Mp N, Mp N, Mp N 3 16 19 Bending mode (Group 4): Mp N, Mp N, Mp N 5 8 12 Bending mode (Group 5): Mp N, Mp N, Mp N 6 9 13 Zero-energy stress mode (Group 6): Mp N, Mp N, Mp N, Mp N, Mp N 4 7 20 21 22 More high-order stress modes can be classified into the 6 stress mode groups above by using the classification method. If the flexibility matrix [H] is a diagonal matrix, the classification of the stress modes is unique (see Section 6). ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4324 W. FENG, S. V. HOA AND Q. HUANG 4.2. 3-D, 8-node solid hybrid element 4.2.1. Determination of optimal stress matrix from the existing stress matrix derived by isofunction method. The 3-D, 8-node solid element has (n") 24 degrees of freedom and (r") 6 rigid-body modes. So it has (m"n!r") 18 natural deformation modes. By means of iso-function method,23 an initial stress matrix [P] can be derived from the assumed displacement field of the element as follows: z 1 0 0 0 0 0 x 0 0 0 0 0 y 0 0 0 0 0 0 1 0 0 0 0 0 x 0 0 0 0 0 y 0 0 0 0 0 0 0 1 0 0 0 [P] " ISO 0 0 0 1 0 0 0 0 x 0 0 0 0 0 y 0 0 0 0 0 0 0 0 x 0 0 0 0 0 y 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 x 0 0 0 0 0 y 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 x 0 0 0 0 0 y 0 0 0 0 0 z 0 0 0 0 z 0 0 0 z 0 0 z 0 0 0 0 0 0 z xy 0 0 0 0 yz 0 0 0 0 zx 0 0 0 0 0 xy 0 0 0 0 yz 0 0 0 0 zx 0 0 0 0 0 xy 0 0 0 0 yz 0 0 0 0 zx 0 0 0 0 0 0 0 0 0 0 yz 0 0 0 0 zx 0 0 0 0 xy 0 0 0 0 0 0 0 0 0 0 zx 0 0 0 0 xy 0 0 0 0 yz 0 0 0 0 0 (26) There are 39 stress modes to be classified in the stress matrix. The number of stress modes is larger than m ("18). The eigenvalue examination shows that the eigenvalues and eigenvectors of the hybrid element stiffness matrix corresponding to the assumed stress matrix [P] are the ISO same as that of displacement element stiffness matrix. We take the stress modes in the matrix [P] as initial stress modes to be classified. The 39 stress modes in the matrix [P] are ISO ISO numbered as follows: 1 0 0 Mp p p p p p N" 1 2 3 4 5 6 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 x 0 0 x 0 0 Mp p p p p p N" 7 8 9 10 11 12 0 0 0 0 0 0 Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) 0 0 0 0 0 x 0 0 0 0 x 0 0 0 0 x 0 0 0 0 x 0 0 (27) 0 (28) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4325 y 0 0 Mp p p p p p N" 13 14 15 16 17 18 0 0 0 0 0 0 0 0 y 0 0 0 0 0 y 0 0 0 0 0 y 0 0 0 0 0 y 0 0 0 0 0 y (29) z 0 0 Mp p p p p p N" 19 20 21 22 23 24 0 0 0 0 0 0 0 0 z 0 0 0 0 0 z 0 0 0 0 0 z 0 0 0 0 0 z 0 0 0 0 0 z (30) xy 0 0 xy 0 0 Mp p p p p N" 25 26 27 28 29 0 0 0 0 0 0 0 0 0 0 xy 0 0 0 0 0 0 xy 0 0 0 xy 0 0 (31) yz 0 0 0 0 0 yz 0 0 0 0 0 yz 0 0 Mp p p p p N" 30 31 32 33 34 0 0 0 yz 0 0 0 0 0 0 0 0 0 0 yz (32) zx 0 0 zx 0 0 Mp p p p p N" 35 36 37 38 39 0 0 0 0 0 0 (33) 0 0 0 0 zx 0 0 0 zx 0 0 0 zx 0 0 0 0 0 In order to reduce the number of stress modes in the assumed stress matrix [P] , these stress ISO modes have to be classified one by one in the order from low-order term to high-order term. Following steps 2—12 in the procedure of the classification, one can obtain (m") 18 representative ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4326 W. FENG, S. V. HOA AND Q. HUANG stress modes Mp p p p p p p p p p p p p p p p p p N corresponding to 18 1 2 3 4 5 6 8 9 11 13 15 18 19 20 22 27 30 36 natural deformation modes. These representative stress modes form a optimal stress matrix [P ] 1 from the existing stress matrix [P] as follows: *40 [P ]"[p p p p p p p p p p p p p p p p p p ] 1 1 2 3 4 5 6 13 20 9 19 8 15 22 11 18 30 36 27 1 0 0 0 0 0 y 0 0 1 0 0 0 0 0 z 0 z 0 0 0 0 yz 0 0 0 0 x 0 0 0 0 0 zx 0 0 0 0 1 0 0 0 0 0 x 0 0 y 0 0 0 0 0 xy 0 0 0 1 0 0 0 0 0 0 0 0 z 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 x 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 y 0 0 0 " 0 (34) This stress matrix is the same as that proposed by Pian.16 4.2.2. Classification of other stress modes. Following steps 13 and 14 in the procedure, other stress modes that remain in the existing stress matrix [P] can be classified into m#1 ("19) *40 stress mode groups as follows: Tension and compressive modes (3 groups): [Mp N , Mp N , Mp N ] 1 G1 2 G2 3 G3 Pure shear modes (3 groups): [Mp N , Mp N , Mp N ] 4 G4 5 G5 6 G6 Bending modes (6 groups): [Mp p N , Mp p N , Mp p N , 8 16 G7 9 24 G8 13 10 G9 Mp p N , Mp p N , Mp p N ] 15 23 G10 19 12 G11 20 17 G12 Torsion modes (3 groups): [Mp N , Mp N , Mp N ] 11 G13 18 G14 22 G15 Saddle modes (3 groups): [Mp p p N , Mp p p N , Mp p p N ] 29 30 38 G16 28 33 36 G17 27 34 39 G18 Zero-energy stress modes (1 group): [Mp , p , p , p , p , p , p , p , p N ] 7 14 21 25 26 31 32 35 37 G19 The first 18 stress mode groups correspond to (m"n!r") 18 natural deformation modes and the last group corresponds to rigid-body modes. There are many other ways to derive initial stress modes. For example, in the assumed stress matrix presented in Reference 14, there are 12 stress modes that are different from stress modes Mp N!Mp N. These stress modes can be 1 39 expressed as follows: Tension and compressive modes: 1 1 [1 1 [1 [1 1 0 2 Mp N" , Mp N" , Mp N" 40 41 42 0 0 0 0 0 0 0 0 0 Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) (35) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4327 Symmetric bending modes: z z 0 Mp N" , 43 0 0 0 0 y x 0 x y Mp N" , Mp N" 44 45 0 0 0 0 0 0 (36) 0 x [x Mp N" , 47 0 0 0 (37) Anti-symmetric bending modes: z [z 0 Mp N" , 46 0 0 0 y 0 !y Mp N" 48 0 0 0 Torsion modes: 0 0 0 0 0 0 0 0 0 Mp N" , Mp N" , Mp N" 49 50 51 z z z x [x x y 0 !2y (38) In the stress matrices derived by means of the symmetric group theory,8,15 there are 11 stress modes that are different from stress modes Mp N—Mp N. They can be expressed in the form. 1 51 Tension and compressive mode: 0 1 [1 Mp N" 52 0 0 0 ( 1997 John Wiley & Sons, Ltd. (39) Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4328 W. FENG, S. V. HOA AND Q. HUANG Torsion modes: 0 0 0 Mp N" 53 0 [x y (40) Bending modes: 2x 0 0 2y 0 0 0 0 2z Mp p p N" 54 55 56 0 [y [x 0 [z [y 0 [x [z 0 0 0 0 0 0 0 0 0 and Mp p p N" 57 58 59 x 0 y 0 [z [y x [z 0 (41) Saddle modes: 0 0 0 0 0 0 0 0 0 Mp p p N" 60 61 62 [2yz x2#y2 [2xz y2#z2 [2xy [2xz [2xy x2#z2 [2yz (42) Other stress modes may be also needed in an assumed stress matrix in order to describe special stress distribution in a local region of a structure to be analysed. For instance, Bending modes: z z z Mp N" , 63 0 0 0 Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) x y x y x y Mp N" , Mp N" 64 65 0 0 0 0 0 0 (43) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4329 Saddle modes: yz yz yz Mp N" , 66 0 0 0 zx xy zx xy zx xy Mp N" , Mp N" 67 68 0 0 0 0 0 0 (44) z2 z2 0 Mp N" 69 0 0 0 (45) Tension and compressive mode: According to steps 13 and 14 in the proposed procedure of classification, the 30 new stress modes Mp N—Mp N can be classified into different stress mode groups above as follows: 40 69 Tension and compressive modes (3 groups): [Mp , p , p N , Mp , p N , Mp , p , p N ] 1 40 69 G1 2 41 G2 3 42 52 G3 Pure shear modes (3 groups): [Mp N , Mp N , Mp N ] 4 G4 5 G5 6 G6 Bending modes (6 groups): [Mp , p , p , p , p N , Mp , p , p , p N , 8 16 44 54 64 G7 9 24 47 57 G8 Mp , p , p , p , p N , Mp , p , p , p N , 13 10 45 55 65 G9 15 23 48 58 G10 Mp , p , p , p , p N , Mp , p , p , p N ] 19 12 43 56 63 G11 20 17 46 59 G12 Torsion modes (3 groups): [Mp , p N , Mp , p N , Mp , p , p N ] 11 49 G13 18 50 G14 22 51 53 G15 Saddle modes (3 groups): [Mp , p , p , p , p N , Mp , p , p , p , p N , 29 30 38 66 60 G16 28 33 36 67 61 G17 Mp , p , p , p , p N ] 27 34 39 68 62 G18 Zero-energy stress modes (1 group): [Mp , p , p , p , p , p , p , p , p N ] 7 14 21 25 26 31 32 35 37 G19 More stress modes can be classified into the stress mode groups above. If the flexibility matrix [H] is a diagonal matrix, the stress modes are uncoupled and the classification of the stress modes is unique (see Section 6). Otherwise, some stress modes may appear in more than one group. 5. CONSTRUCTION OF ASSUMED STRESS MATRICES As shown above, by means of the proposed procedure for the classification of stress mode, stress modes can be classified into m ("n!r) stress mode groups corresponding to m natural deformation modes and a zero-energy mode group corresponding to rigid-body modes. Each natural deformation mode is related to a stress mode group except zero-energy mode group, and each stress mode group may contain many different stress modes that are interchangeable in the stress matrix [P]. ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4330 W. FENG, S. V. HOA AND Q. HUANG The classification of stress modes reveals the relationship among the different stress modes that are used in the different stress matrices for any type of a hybrid element proposed by different researchers. In order to avoid kinematic deformation mode, the stress matrix [P] must contain m stress modes at least. No matter how many stress modes there are in the stress matrix [P], the order of the stiffness matrix is equal to or less than m. Therefore, m stress modes are necessary and sufficient to form a stress matrix for avoiding kinematic deformation modes in the hybrid element. Moreover, in view of the classification of stress modes, the m stress modes in the stress matrix [P] must come from m different stress mode groups except zero-energy mode group. Thus, for a hybrid element to be free from kinematic deformation mode, we have ¹he necessary and sufficient condition: The number of stress modes in an assumed stress matrix must be equal to or more than m ("n!r) and at least m stress modes in the stress matrix [P] must be chosen from m different stress mode groups corresponding to m natural deformation modes of an element which has n degrees of freedom and r rigid-body modes. In this statement, the necessary condition is that the number of stress modes for a hybrid element must be equal to or more than m ("n!r). It was presented by Veubeke9 and Pian.10 The sufficient condition is that the stress matrix [P] must contain m stress modes chosen from m different stress mode groups corresponding to m natural deformation modes. This condition explains why in some examples there exist kinematic deformation modes even when the necessary condition (m@'n!r) is satisfied. In these examples, the stress modes in the stress matrix [P] do not come from m different stress mode groups except the zero-energy mode group. For a hybrid element, overuse of stress modes will result in over-rigid element,7 and will cost more computational time because the calculation of element stiffness matrix requires an inversion of the flexibility matrix [H]. Therefore, an assumed stress field, its stress matrix contains m ("n!r) least-order stress modes and its resulting finite element is free from kinematic deformation modes, is considered to be best and is optimal with respect to computer resources.7,8 By means of the m classified stress mode groups and the necessary and sufficient condition, this kind of stress matrices can be constructed. Furthermore, it is convenient to construct an assumed stress matrix according to the problem to be solved because there are many stress modes in every stress mode group for choice. The procedure of constructing stress matrix is presented as follows: 1. Using the iso-function method,23 one can derive a number of initial stress modes to be classified. 2. One may put the initial stress modes one by one into stress matrix [P] in the order from low-order term to high-order term. By means of the classification method, one can obtain m representative stress modes corresponding to m natural deformation modes. These representative stress modes form a optimal stress matrix [P] from the existing stress 015 matrix [P] . *40 3. One may obtain other initial stress modes derived by different methods. Following the steps 13 and 14 in the procedure of the classification, one can classify all initial stress modes into m#1 different stress mode groups. 4. By means of the m#1 classified stress mode groups and the necessary and sufficient condition above, many stress matrices [P] can be constructed according to the problem to be solved. It is necessary to choose one stress mode at least from each group except the zero-energy mode group in order to avoid kinematic deformation modes. Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) ( 1997 John Wiley & Sons, Ltd. 4331 CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS The necessary steps have been illustrated in the section above. The following gives some examples to illustrate the procedure for constructing a stress matrix [P] which has minimum number of stress modes. 5.1. 2-D, 4-node plane element By means of the m#1 stress mode groups classified above and the necessary and sufficient condition for avoiding kinematic deformation modes, one can choose one stress mode from each stress mode group except zero-energy mode group to form a stress matrix. For example, 1 1 0 [x 0 [P ]"[p p p p p ]" 1 [1 0 [y III 10 11 3 12 13 0 0 1 x 0 (46) y and 1 0 y 1 0 [P ]"[p p p p p ]" 1 [1 0 0 x IV 10 11 3 7 5 0 0 1 0 0 (47) Five stress modes in each stress matrix come from five different stress mode groups corresponding to five natural deformation modes. The two stress matrices are the same as that proposed by Atluri.8,15 More stress matrices can be constructed on purpose. For example, 1 0 y [x 1 [P ]"[p p p p p ]" 1 [1 0 0 V 10 11 3 7 13 0 0 1 0 0 (48) y and x2 [P ]"[p p p p p ]" 0 VI 14 11 3 12 13 0 1 0 0 [1 0 [y 0 1 x [x 0 (49) y The eigenvalue examination shows that the hybrid element constructed by [P ]—[P ] are free I VI from kinematic deformation modes as shown in Table I. In the last column of the table, the Table I. Eigenvalues of element stiffness matrices (2-D, 4-node plane element, v"0·3) [P ] I 0·4945 0·4945 0·7692 0·7692 1·4290 [P ] II [P ] III [P ] IV [P ] V [P ] VI Disp. 0·3333 0·3333 0·7692 0·7692 1·4290 0·09259 0·09259 0·7692 0·7692 1·4290 0·3333 0·3333 0·7692 0·7692 1·4290 0·09259 0·3333 0·7692 0·7692 1·4290 0·09259 0·09259 0·7692 0·7692 1·4290 0·4945 0·4945 0·7692 0·7692 1·4290 ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4332 W. FENG, S. V. HOA AND Q. HUANG eigenvalues of displacement element stiffness matrix are given. If stress modes in a stress matrix [P] come from m ((m) stress mode groups, the hybrid element will have kinematic deformation 1 modes even if the number of stress modes is larger than m. This is why a hybrid element contains kinematic deformation modes when the necessary condition (m@'n!r) is satisfied. A stress matrix [P] must have m stress modes corresponding to m natural deformation modes of an element. 5.2. 3-D, 8-node solid element Using the same way as 2-D case, one can choose m stress modes from m classified stress mode groups except zero-energy mode group above to form the eight stress matrices [P ]—[P ] 2 9 proposed by Atluri et al.8,15 as follows: [P ]"[p p p p p p p p p p p p p p p p p p ] 2 40 41 42 4 5 6 49 50 53 54 55 56 57 58 59 60 61 62 1 1 0 0 0 0 0 0 0 2x 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 2y 0 0 0 0 0 0 0 [1 0 0 0 0 0 0 0 0 2z 0 0 0 0 0 0 z 0 0 y x 0 [2xz [2yz x2#y2 [2xy [2xz x2#z2 [2yz 1 [1 1 0 0 0 0 1 0 0 z 0 0 0 0 1 0 x [x [x 0 0 0 0 0 1 y " 0 y [y [x [z [y 0 [z [z [y y2#z2 0 [x [z 0 x 0 [2xy (50) [P ]"[p p p p p p p p p p p p p p p p p p ] 3 40 41 42 4 5 6 49 50 53 54 55 56 57 58 59 27 36 30 1 0 0 0 0 0 0 0 2x 0 0 0 0 0 0 0 yz 1 0 0 0 0 0 0 0 2y 0 0 0 0 0 xz 0 [1 0 0 0 0 0 0 0 0 2z 0 0 0 xy 0 0 z 0 0 y x 0 0 0 0 0 0 0 0 1 1 [1 1 0 0 0 0 1 0 0 z 0 0 0 0 1 0 x [x [x 0 0 0 0 0 1 y " !y [x y 0 0 [z [y [z 0 p 53 p 53 p 53 p 53 p 53 p 53 p 54 p 54 p 45 p 45 p 45 p 45 0 [x [z [z [y 0 x 0 0 p p 46 p 46 p 59 p 59 p 46 p 46 p 60 p 27 p 60 p 27 p 60 p 27 p 61 p 36 p 61 p 36 p 61 p 36 (51) and [P ]"[p 4 40 [P ]"[p 5 40 [P ]"[p 6 40 [P ]"[p 7 40 [P ]"[p 8 40 [P ]"[p 9 40 p p p p p p 41 41 41 41 41 41 p 42 p 42 p 42 p 42 p 42 p 42 p p 4 5 p p 4 5 p p 4 5 p p 4 5 p p 4 5 p p 4 5 p 6 p 6 p 6 p 6 p 6 p 6 p 49 p 49 p 49 p 49 p 49 p 49 p p p p p p Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) 50 50 50 50 50 50 p 55 p 55 p 44 p 44 p 44 p 44 p 56 p 56 p 43 p 43 p 43 p 43 p 48 p 48 p 57 p 57 p 48 p 48 p p p p p 47 47 58 58 47 47 p ] 62 p ] 30 p ] 62 (52) p ] 30 p ] 62 p ] 30 ( 1997 John Wiley & Sons, Ltd. 4333 CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS The assumed stress matrix given by Huang14 also can be formed by means of the same way, [P ]"[p p p p p p p p p p p p p p p p p p ] 10 40 41 42 4 5 6 43 44 45 46 47 48 49 50 51 27 36 30 0 y 0 0 0 yz 0 0 1 1 [1 0 0 0 z 0 y z 1 [1 [1 0 0 0 z x 0 [z x 0 0 0 0 zx 0 [x !y 0 0 0 0 0 xy z z 0 0 0 0 1 0 2 0 0 0 0 x y 0 0 0 0 1 0 0 0 0 0 0 0 0 z 0 0 0 0 1 0 0 0 0 0 0 0 x [x [x 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 y [2y 0 0 0 " 0 (53) Moreover, many stress matrices [P] can be constructed on purpose. Three new stress matrices are given as follows: [P* ]"[p p p p p p p p p p p p p p p p p p ] 1 40 41 42 4 5 6 63 64 65 46 47 48 49 50 51 66 67 68 0 y 0 0 0 yz xz xy 1 1 [1 0 0 0 z x y z 1 [1 [1 0 0 0 z x y [z x 0 0 0 yz xz xy [x [y 0 0 0 yz xz xy z z 0 0 0 0 1 0 2 0 0 0 z x y 0 0 0 0 1 0 0 0 0 0 0 0 0 z 0 0 0 0 1 0 0 0 0 0 0 0 x [x [x 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 y [2y 0 0 0 " 0 (54) [P*]"[p p p p p p p p p p p p p p p p p p ] 2 40 41 42 4 5 6 43 44 45 57 58 59 22 11 18 30 36 27 1 0 0 0 0 0 z 0 y 0 0 0 0 0 0 yz 0 0 0 1 0 0 0 0 z x 0 0 0 0 0 0 0 0 xz 0 0 0 1 0 0 0 0 x y 0 0 0 0 0 0 0 0 xy 0 0 0 1 0 0 0 0 0 y x 0 z 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 [z [y 0 x 0 0 0 0 0 0 0 (55) " 0 0 0 0 0 1 0 0 0 [z x 0 0 0 y [P*]"[p p p p p p p p p p p p p p p p p p ] 3 69 41 42 4 5 6 63 64 65 46 47 48 49 50 51 66 67 68 z2 1 [1 0 0 0 z x y z z2 [1 [1 0 0 0 z x y [z 0 y 0 0 0 yz xz xy x 0 0 0 0 yz zx xy [x [y 0 0 0 yz xz xy z z 0 0 0 0 0 2 0 0 0 z x y 0 0 0 0 1 0 0 0 0 0 0 0 0 z 0 0 0 0 1 0 0 0 0 0 0 0 x [x [x 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 y [2y 0 0 0 " ( 1997 John Wiley & Sons, Ltd. 0 (56) Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4334 W. FENG, S. V. HOA AND Q. HUANG Table II. Eigenvalues of element stiffness matrices (3-D, 8-node solid element, v"0·3) [P ] 2 0·07123 0·07123 0·07123 0·1282 0·1282 0·1282 0·1282 0·1282 0·07246 0·07246 0·07246 0·5128 0·7692 0·7692 0·7692 0·7692 0·7692 2·5000 [P ] 4 [P ] [P ] [P ] 1 9 10 [P ], [P* ] 7 2 [P* ] 1 [P* ] 3 0·07123 0·07123 0·07123 0·2564 0·2564 0·2564 0·1282 0·1282 0·07264 0·07264 0·07264 0·5128 0·7692 0·7692 0·7692 0·7692 0·7692 2·5000 0·1111 0·1111 0·1111 0·2564 0·2564 0·2564 0·1282 0·1282 0·4762 0·4762 0·4762 0·5128 0·7692 0·7692 0·7692 0·7692 0·7692 2·5000 0·1111 0·1111 0·1111 0·1282 0·1282 0·1282 0·1282 0·1282 0·4762 0·4762 0·4762 0·5128 0·7692 0·7692 0·7692 0·7692 0·7692 2·5000 0·09259 0·09259 0·09259 0·2564 0·2564 0·2564 0·1282 0·1282 0·5556 0·5556 0·5556 0·5128 0·7692 0·7692 0·7692 0·7692 0·7692 2·5000 0·09259 0·09259 0·09259 0·2564 0·2564 0·2564 0·1282 0·1282 0·5556 0·5556 0·5556 0·5128 0·7692 0·7692 0·7692 0·7692 0·7692 0·8065 The results of eigenvalue examination are given in Table II. It shows that each of the stiffness matrices constructed by the assumed stress matrices [P ]—[P ], [P* ], [P* ] and [P* ] has 1 10 1 2 3 m non-zero eigenvalues. The resulting hybrid elements do not have any kinematic deformation modes. More assumed stress matrices can also be constructed by means of this method. If one stress mode group is missed except the zero-energy mode group in the process of choosing stress modes, the hybrid element will contain kinematic deformation modes. In the previous work, it is proposed to suppress kinematic deformation modes by adding stress modes of high-order term. Actually, it cannot guarantee that all kinematic deformation modes are suppressed. If the high-order stress modes do not belong to the stress mode groups which are missed in the construction of the assumed stress matrix except the zero-energy mode group, adding stress modes of high-order term cannot improve the hybrid element any more. Moreover, overuse of stress modes will result in over-rigid elements.7 Therefore, an ideal situation is to choose m("n!r) least-order stress modes, but with the suppression of all kinematic deformation modes. Thus, an assumed stress matrix [P] can be constructed by choosing m stress modes from m stress mode groups that correspond to m natural deformation modes. 6. UNIQUENESS OF STRESS MODE GROUPS When stress modes are classified, it is observed that if the flexibility matrix [H] is a diagonal matrix, the classification of stress modes is unique; otherwise, some stress modes may be interchangeable between two stress mode groups. For example, the stress modes Mp N and Mp N 1 2 for 2-D, 4-node plane hybrid element may be interchanged between group 1 and group 2. This makes the flexibility matrix [H] not diagonal when the stress matrix [P] consists of (Mp N, Mp N, 1 2 Mp N, Mp N, Mp N). Therefore, the first two stress mode groups for 4-node plane element may 3 5 6 Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4335 become Tension mode (Group 1): Mp N (or Mp N), Mp N, Mp N, Mp N 1 2 10 14 17 Tension mode (Group 2): Mp N (or Mp N), Mp N, Mp N, Mp N 2 1 11 15 18 However, the stress modes Mp N and Mp N for the plane element cannot be interchanged between 10 11 the two groups because the matrix [H] is diagonal when the stress matrix [P] consists of (Mp N, 10 Mp N, Mp N, Mp N, Mp N). In fact, if the flexibility matrix [H] is a diagonal matrix, the stress modes 11 3 5 6 that form the stress matrix [P] are a set of uncoupled stress modes. It has been proved by Huang14 that if the matrix [H] is a diagonal matrix, the stiffness matrix satisfies the superposition principle: m [K]" + [K ] i i/1 (57) where [K ]"[G ]T[H ]~1[G ] i i i i Pv [Pi ]T[S][Pi ] d» [G ]" [P ]T[B] d» i Pv i [H ]" i (58) and [K]"[G]T[H]~1[G] Pv [P]T[S][P] d» [G]" [P]T[B] d» Pv [H]" (59) in which [P]"[Mp N Mp N Mp N . . . Mp N] 1 2 3 m [P ]"[M0N M0N . . . Mp N . . . M0N] i i (60) G1 G [G]" 2 ... G m Therefore, the elastic energy of the element is decomposable if the flexibility matrix [H] is a diagonal matrix. ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4336 W. FENG, S. V. HOA AND Q. HUANG 6.1. Classification condition of stress modes A hybrid element stiffness matrix [K] can be formulated using equations (59) and (60). Its eigenvalues and eigenvectors are calculated from equation (9). The eigenvectors Md N (i"1, i 2, . . . , m) satisfy the condition (10). Using equations (9) and (10), the eigenvalue equation is changed to j "Md NT[K] Md N i i i (61) For any stress mode Mp N among m stress modes Mp , p , . . . , p N, the stiffness matrix [K ] can j 1 2 m j be derived using equations (58) and (60). Corresponding to the ith natural deformation mode, we have j "Md NT[K ]Md N i i j i (62) According to the classification condition of stress modes, if the stress mode Mp N belongs to the j ith stress mode group corresponding to the natural deformation mode Md N, the eigenvalue j is i i a non-zero value; otherwise, the eigenvalue j equals zero. This condition can be expressed in the i form Md NT[K ]Md N"0, i j i Md NT[K ]Md N"j , i j i i iOj (63) i"j If the stress mode Mp N is a zero-energy stress mode, all eigenvalue j (i"1, 2, . . . , m) equal zero. j i ¹heorem 1. If and only if the flexibility matrix [H] is a diagonal matrix, the eigenvalues obtained from separate mode equations ( [K ]!j[I])MdN"0, i"1, 2, 3, . . . , m i (64) are equal to the eigenvalues obtained from the total equation ( [K]!j[I])MdN"0 (65) In which, the matrices [K ] and [K] are defined in equations (58) and (59). ¹his was stated as i a postulate in Reference.14 Proof. From the equations (61), (59) and (15), we have j"Md NT[K]Md N"Mb NT[H]Mb N i i i i (66) Because the matrix [H] is a diagonal matrix, we have m [H]" + [H ] j j/1 (67) Thus, using equations (12), (13) and (58), the eigenvalue of the matrix [K] is m j" + Mb NT[H ]Mb N"Mb NT[H ]Mb N i j i i i i j/1 Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) (68) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4337 Furthermore, using equations (15) and (63), we obtain j"Mb NT[H ]Mb N"Md NT[K ]Md N"j K (69) i i i i i i i ¹heorem 2. If the flexibility matrix [H] is a diagonal matrix, the classification of m stress modes is unique. Proof. If a stress mode among m stress modes that form the stress matrix [P] appears in more than one stress mode group, one of the m stress mode groups must contain two stress modes. Assume that the stress modes Mp N and Mp N belong to the ith stress mode group Mb N correspondi j i ing to the natural deformation modes Md N. Thus, we have i j "Md NT[K ]Md N and j "Md NT [K ]Md N (70) ii i i i ij i j i Corresponding to the natural deformation mode Md N, we can obtain the eigenvalue of the i stiffness matrix [K] formulated by m stress modes as follows: j"Md NT [K]Md N (71) i i Because the flexibility matrix [H] is diagonal, the stiffness matrix satisfies the superposition principle. From equations (57) and (70), we obtain m j"Md NT [K]Md N" + Md NT [K ]Md N"j #j i i i k i ii ij k/1 using Theorem 1, we have j"j ii From the equations (72) and (73), we obtain (72) (73) j "0 (74) ij According to the condition of classification, the stress modes Mp N does not belong to the ith j stress mode group. The stress mode group Mb N only contains Mp N. Therefore, the stress modes i i Mp N cannot appear in two stress mode groups Mb N and Mb N. Thus, if the matrix [H] is diagonal, j i j the classification of m stress modes is unique. K 7. CONCLUSION A new method for classifying stress modes in assumed stress matrices is presented. It is assumed that there are m("n!r) natural deformation modes of an element which has n degrees of freedom and r rigid-body modes. For any type of hybrid element, all stress modes in various stress matrices derived by different methods can be classified into m stress mode groups corresponding to m natural deformation modes and a zero-energy mode group corresponding to rigid-body modes. If the flexibility matrix [H] is diagonal, the deformation energy of the element is decomposable and the classification of stress modes is unique. The necessary and sufficient condition for avoiding kinematic deformation modes is that an assumed stress matrix [P] must contain m stress modes chosen from m different stress mode groups, except zero-energy ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997) 4338 W. FENG, S. V. HOA AND Q. HUANG mode group. The reason of the existence of kinematic deformation modes when the criterion (m@'n!r) is satisfied is that the stress modes in the assumed stress matrix [P] are not chosen from m different stress mode groups corresponding to m natural deformation modes. The classification method can be applied to any type of hybrid elements and be used for two purposes: 1. to determine the optimal stress matrix from the existing stress matrix [P] or any other *40 stress matrix [P] derived using other method, and classify stress modes into m different stress mode groups; 2. to construct many new assumed stress matrices by using minimum number of stress modes according to the problems to be analysed. These stress matrices are without zero-energy stress modes, and the resulting element stiffness matrices are free from kinematic deformation modes. The classification of stress modes reveals the relationship among the different assumed stress fields for any type of hybrid element proposed by different researchers. An assumed stress matrix [P], which consists of m("n!r) least-order stress modes and results in the element stiffness matrix without kinematic deformation modes, is considered to be best and is optimal with respect to computer resources because overuse of stress modes will result in over-rigid element and cost more computational time. REFERENCES 1. T. H. H. Pian, ‘Derivation of element stiffness matrices by assumed stress distributions’, AIAA J., 2, 1333—1336 (1964). 2. T. H. H. Pian, ‘State-of-the-art development of hybrid/mixed finite element method’, Finite Elements Anal. Des., 21, 5—20 (1995). 3. T. H. H. Pian and P. Tong, ‘Mixed and hybrid finite element methods’, in H. Kardestuncer and D. H. Norrie (ed.), Finite Element Handbook, McGraw-Hill, New York, 1987. 4. W. Feng and S. V. Hoa, ‘A partial hybrid degenerated plate/shell element for the analysis of laminated composites’, Int. J. Numer. Meth. Engng., 39, 3625—3639 (1996). 5. J. Han and S. V. Hoa, ‘A three-dimensional multilayer composite finite element for stress analysis of composite laminates’, Int. J. Numer. Meth. Engng., 36, 3903—3914 (1993). 6. S. V. Hoa and W. Feng, ‘Finite elements for analysis of composite structures’, in S. V. Hoa (ed.), Computer-aided Design of Polymer-Matrix Composite Structures, Marcel Dekker, New York, 1995. 7. R. D. Henshell, ‘On hybrid finite elements’, Proc. Brunel ºniversity Conf. of the Institute of Math and Its Applications, April, 1972. 8. E. F. Punch and S. N. Atluri, ‘Development and testing of stable, invariant, isoparametric curvilinear 2- and 3-D hybrid-stress elements’, Comput. Meth. Appl. Mech. Engng., 47, 331—356 (1984). 9. B. M. Fraeijs de Veubake, ‘Bending and stretching of plates—special models for upper and lower bounds’, Proc. Conf. on Matrix Methods in Structure Mechanics, AFFDL-TR-M66-80, 863—886 (1965). 10. T. H. H. Pian and P. Tong, ‘Basis of finite element methods for solid continua’, Int. J. Numer. Meth. Engng., 1, 3—28 (1969). 11. F. Brezzi, ‘On the existence, uniqueness and approximation of saddle point problems arising from Lagrange multipliers’, RAIRO 8(R2), 129—151 (1974). 12. I. Babuska, J. T. Oden and J. K. Lee, ‘Mixed-hybrid finite element approximations of second-order elliptic boundary-value problems’, Comput. Meth. Appl. Mech. Engng., 11, 175—206 (1977). 13. S. Ahmad and B. M. Irons, ‘An assumed stress approach to refined isoparametric finite elements in three dimensions’, Finite Element Method in Engineering, University of New South Wales, 1974, pp. 85—100. 14. Q. Huang, ‘Three dimensional composite finite element for stress analysis of anisotropic laminate structures’, Ph.D. ¹hesis, Concordia University, Montreal, Canada, 1989. 15. R. Rubinstein, E. F. Punch and S. N. Atluri, ‘An analysis of, and remedies for, kinematic models in hybrid-stress finite elements: Selection of stable, invariant stress fields’, Comput. Meth. Appl. Mech. Engng., 38, 63—92 (1983). 16. T. H. H. Pian and D. P. Chen, ‘On the suppression of zero energy deformation modes’, Int. J. Numer. Meth. Engng., 19, 1741—1752 (1983). Int. J. Numer. Meth. Engng., 40, 4313—4339 (1997) ( 1997 John Wiley & Sons, Ltd. CLASSIFICATION OF STRESS MODES IN ASSUMED STRESS FIELDS 4339 17. T. H. H. Pian and P. Tong, ‘Relations between incompatible displacement model and hybrid stress model’, Int. J. Numer. Meth. Engng., 22, 173—181 (1986). 18. T. H. H. Pian and C. C. Wu, ‘A rational approach for choosing stress terms for hybrid finite element formulations’, Int. J. Numer. Meth. Engng., 26, 2331—2343 (1988). 19. K. Y. Sze, C. L. Chow and W. J. Chen, ‘A rational formulation of iso-parametric hybrid stress element for three dimensional stress analysis’, Finite Elements Anal. Des., 7, 61—72 (1990). 20. W. J. Chen and Y. K. Cheung, ‘Three-dimensional 8-node and 20-node refined hybrid isoparametric elements’, Int. J. Numer. Meth. Engng., 35, 1871—1889 (1992). 21. K. Y. Sze, ‘Control of spurious mechanisms for 20-node and transition sub-integrated hexahedral elements’, Int. J. Numer. Meth. Engng., 37, 2235—2250 (1994). 22. C. C. Wu and Y. K. Cheung, ‘On optimization approaches of hybrid stress elements’, Finite Elements Anal. Des., 21, 111—128 (1995). 23. J. Han and S. V. Hoa, ‘Three dimensional composite elements for finite element analysis of anisotropic laminated structures’, Proc. 2nd Int. Conf. on Computer Aided Design in Composite Material ¹echnology, Computational Mechanics Publications, 1990, pp. 189—200. . ( 1997 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 4313— 4339 (1997)

1/--страниц