- Email: [email protected]
A non-linear corotational 4-node plane element
Mechanics Research Communications 35 (2008) 408–413
Contents lists available at ScienceDirect
Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom
A non-linear corotational 4-node plane element Jean-Marc Battini Department of Mechanics, KTH, Royal Institute of Technology, SE-10044 Stockholm, Sweden
a r t i c l e
i n f o
Article history: Received 14 September 2007 Received in revised form 28 February 2008 Available online 16 March 2008
a b s t r a c t A non-linear corotational 4-node plane element is presented. The main interested of the approach is that accurate and efficient linear elements can be transformed to non-linear elements with little work and without considering the assumptions used to derive these linear elements. Two numerical examples are presented. Ó 2008 Elsevier Ltd. All rights reserved.
Keywords: Corotational formulation Plane elements Non-linear analyses
1. Introduction Due to their simplicity and their efficiency in large-scale analyses, low-order finite elements are often used in structural mechanics. For this reason, a lot of research has been done in order to develop accurate and efficient 4-node plane elements, see e.g. the paper of Fredriksson and Saabye Ottosen (2004) for a complete review. However, many of the proposed solutions concern only linear formulations. The reason is that these linear elements are by themselves complicated and their transformation to non-linear formulations is not always straightforward (see e.g. Simo and Armero, 1992; Simo et al., 1993). This is particularly the case if the template formulation proposed by Felippa (2006) is used. One simple way, at least for large displacements and small strains problems, to derive non-linear elements is to used the corotational approach. Corotational beams and shells elements have been widely studied, see e.g. the review article of Felippa and Haugen (2005). However, to the author’s knowledge, no complete work and no numerical analyses have been done with plane or solid elements. Besides, this method seems particularly interesting in the present context since it transforms linear elements to non-linear ones with very little work, and without considering the theory used to derive these linear elements. The main idea of the corotational approach is to decompose the motion of the element into rigid-body and pure deformational parts, through the use of a reference system which continuously rotates and translates with the element. The deformational response is captured at the level of the local reference frame, whereas the geometric non-linearity induced by the large rigid-body motion, is incorporated in the transformation matrices relating local and global internal force vectors and tangent stiffness matrices. Assuming the pure deformation part to be small, a geometrical linear theory can be used in the local system. The main interest is that for elements with the same number of nodes and degrees of freedom, the transformation matrices are the same. Using this property, different existing high performing geometrical linear elements can be used as local formulations. As mentioned above, this is particularly interesting for 4-node plane elements for which several advanced linear formulations are available in the literature.
E-mail address: [email protected] 0093-6413/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2008.03.002
J.-M. Battini / Mechanics Research Communications 35 (2008) 408–413
409
The corotational framework for 4-node plane elements is given in Section 2. Using for the local formulation two accurate linear elements found in the literature, two numerical problems presenting large displacements and small strains are presented in Section 3. The results obtained with the corotational approach are compared to the ones obtained with ANSYS. 2. Corotational formulation 2.1. Element kinematics The element kinematics is defined in Fig. 1. (O, X, Y) is the global coordinate system. ri ¼ ðX i ; Y i Þ and ui ¼ ðU i ; V i Þ denote the global coordinates and global displacements of the node i. The origin of the local coordinate system (C, x, y) is calculated by rc ¼ ðX c ; Y c Þ
Xc ¼
1 ðX 1 þ X 2 þ X 3 þ X 4 Þ 4
Yc ¼
1 ðY 1 þ Y 2 þ Y 3 þ Y 4 Þ 4
ð1Þ
In the initial configuration, the axes of the local and global coordinate systems are parallel. roi ¼ ðxi ; yi Þ defines the local coordinates of the node i. The motion of the element from the initial to the deformed configuration is split in two steps. The first one is a rigid translation and rotation of the undeformed element. The rigid translation is defined by uc , the global displacement of the point C, which is calculated by uc ¼ ðU c ; V c Þ
Uc ¼
1 ðU 1 þ U 2 þ U 3 þ U 4 Þ 4
Vc ¼
1 ðV 1 þ V 2 þ V 3 þ V 4 Þ 4
ð2Þ
The rigid rotation, specified by the angle h, defines the orientation of the local coordinate system in the deformed configuration. The second step is a local deformation with respect to the local frame. This step is defined by the local deformational dis i ¼ ðu i ; v i Þ; i ¼ 1; . . . ; 4 at the nodes. These displacements are calculated as placements u i xi u cos h sin h X i þ U i X c U c ¼ ð3Þ i yi v sin h cos h Y i þ V i Y c V c The rigid rotation h is taken such that the square of the Euclidean norm of the nodal local deformational displacements X4 2i 2 þ v ð4Þ u i¼1 i is minimized. Introducing Eq. (3) in Eq. (4), taking the variation with respect to h and setting the resulting expression equal to zero, gives P4 ½xi ðY i þ V i Y c V c Þ yi ðX i þ U i X c U c Þ tan h ¼ Pi¼1 ð5Þ 4 i¼1 ½xi ðX i þ U i X c U c Þ þ yi ðY i þ V i Y c V c Þ This equation gives two values for the rotation h, (h and p þ h). One value corresponds to the minimization of expression (4), the other value to its maximization. One simple way to choose between these values is to calculate the expression (4) with both angles and to take the angle which gives the smallest result. Finally, the local and global displacement vectors are 1 pl ¼ ½ u
1 v
2 u
2 v
3 u
3 v
4 u
4 T v
pg ¼ ½ U 1
V1
U2
V2
Fig. 1. Element kinematics.
U3
V3
U4
V 4 T
ð6Þ
410
J.-M. Battini / Mechanics Research Communications 35 (2008) 408–413
2.2. Transformation matrices The corotational framework is defined by a change of variables between local and global quantities. Several existing linear formulations can be taken for the local stiffness matrix Kl and internal force vector f l consistent with pl . The global tangent stiffness matrix Kg and internal force vector f g consistent with pg can be calculated by equating the virtual work in the two systems, as dpTg f g ¼ dpTl f l
ð7Þ
which, introducing the relations dpl ¼ Bdpg
df g ¼ Kg dpg
f l ¼ Kl pl
ð8Þ
gives f g ¼ BT f l
Kg ¼ BT Kl B þ Kh
oðBT f l Þ opg
Kh ¼
ð9Þ fl
Differentiation of Eq. (3) gives i ydi du cos h sin h dU i dU c ¼ þ dh i dv xdi sin h cos h dV i dV c
i þ xi xdi ¼ u i þ yi ydi ¼ v
The expression of dh can be obtained from Eq. (5) as 4 X cos h sin h dU i dU c 1 ½ x y dh ¼ i i 4 P sin h cos h dV i dV c i¼1 ðxi xdi þ yi ydi Þ
ð10Þ
ð11Þ
i¼1
From Eqs. (10) and (11), and omitting dU see Eq. (8), is expressed as B ¼ PET
c
and dV c since finite elements are invariant under pure translation, the matrix B,
P ¼ I AG
ð12Þ
with A ¼ ½ yd1 G¼
4 P
xd1 1
yd2
xd2
½ y1
yd3
x1
xd3
yd4
x2
y3
y2
xd4 x3
T
y4
ð13Þ x4
ð14Þ
ðxi xdi þ yi ydi Þ
i¼1
E ¼ diagðR; R; R; RÞ
R¼
cos h sin h
sin h cos h
ð15Þ
Using Eq. (12), the matrix Kh , see Eq. (9), is defined by Kh ¼ dEPT f l þ EdPT f l
ð16Þ T
After some work, and noting that A G
T
¼ I, it is obtained
Kh ¼ E½FT G GT FPET
ð17Þ
with PT f l ¼ ½ n1 F ¼ ½ n2
n2
n1
n3 n4
n4 n3
n5 n6
n6
n7
n5
n8 T n8
n7
ð18Þ ð19Þ
3. Numerical examples Two linear elements, Qm6 and Qnew, have been used for the local formulation. The Qm6 element has been introduced by Taylor et al. (1976), whereas the Qnew element has been developed by Fredriksson and Saabye Ottosen (2004). Two problems, presenting large displacements and small strains have been analysed. The results obtained with the present corotational formulation have been compared to the results obtained with ANSYS. The analyses in ANSYS have been performed using the element PLANE182 with KEYOPT(1) = 2. This element is a non-linear formulation of the Qm6 element based on the updated Lagrangian formulation. This element is implemented for large strains, but since the strains in the two test problems remain small, it can be used for the proposed comparison. For each structure, two meshes have been tested: a fine one which gives the converged solutions and a coarse one in order to test the ability of the present corotational formulation in obtaining good results with coarse meshes.
411
J.-M. Battini / Mechanics Research Communications 35 (2008) 408–413
3.1. Example 1: angle frame This example is described in Fig. 2. The structure is clamped at one end (all displacements zero) and loaded by a uniformly distributed horizontal force F at the other end. The fine mesh consists of 304 square elements whereas the coarse one consists of 19 distorted elements. The results obtained using 20 steps and the tolerance parameter 105 for the nodal forces
Fig. 2. Example 1: angle frame.
4.5 4 Qnew coarse mesh 3.5 Qm6 coarse mesh
Load F / 10000
3 ANSYS coarse mesh 2.5 fine mesh 2 1.5 1 0.5 0 0
1
2
3
4
5
6
7
Horizontal displacement of point A Fig. 3. Example 1: load–displacement diagram.
Table 1 Example 1 – total number of Newton–Raphson iterations
Coarse mesh Fine mesh
ANSYS
Qm6
Qnew
66 72
60 60
60 60
412
J.-M. Battini / Mechanics Research Communications 35 (2008) 408–413
Fig. 4. Example 2: thin ring.
11 10 9 8
Load F / 1000
7 6 Qnew coarse mesh 5 Qm6 coarse mesh 4 ANSYS coarse mesh 3 fine mesh 2 1 0 0
5
10
15
20
25
30
Vertical displacement of point A Fig. 5. Example 2: load–displacement diagram.
(command CNVTOL in ANSYS) are shown in Fig. 3. With the fine mesh, the three elements give exactly the same results. With the coarse mesh, ANSYS and Qm6 give almost the same results whereas Qnew gives slightly better results. The total number of Newton–Raphson iterations are given in Table 1. It can be observed that the convergence rate of the corotational approach is slightly better than the one in ANSYS. 3.2. Example 2: thin ring A half ring, clamped at both ends is loaded at its top by a concentrated force 2F. Due to the symmetry, only half of the structure is modelled, with 80 4 ¼ 320 elements and 20 1 ¼ 20 elements, see Fig. 4. The results are shown in Fig. 5. The three elements give the same results, both with the coarse and fine meshes. It can also be noted that the results obtained with the coarse mesh are very closed to the ones obtained with the fine mesh. 4. Conclusion This paper has described a corotational framework for 4-node plane elements. The linear elements Qm6 and Qnew, found in the literature, has been used as local formulations. Two numerical applications have shown that the proposed formulation
J.-M. Battini / Mechanics Research Communications 35 (2008) 408–413
413
can be used to analyse plane problems presenting large displacements and small strains. Very good numerical results have been obtained, both with fine and coarse meshes, and also with distorted elements. The main interest of the proposed approach is that, using the described corotational framework, efficient linear elements are automatically transformed to nonlinear formulations. References Felippa, C., 2006. Supernatural QUAD4: A template formulation. Comput. Meth. Appl. Mech. Eng. 195, 5316–5342. Felippa, C., Haugen, B., 2005. A unified formulation of small-strain corotational finite elements: I. Theory. Comput. Meth. Appl. Mech. Eng. 194, 2285–2335. Fredriksson, M., Saabye Ottosen, N., 2004. Fast and accurate 4-node quadrilateral. Int. J. Numer. Meth. Eng. 61, 1809–1834. Simo, J., Armero, F., 1992. Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes. Int. J. Numer. Meth. Eng. 33, 1413–1449. Simo, J., Armero, F., Taylor, R., 1993. Improved versions of assumed enhanced strain tri-linear elements for 3d finite deformation problems. Comput. Meth. Appl. Mech. Eng. 110, 359–386. Taylor, R., Beresford, P., Wilson, E., 1976. A non-conforming element for stress analysis. Int. J. Numer. Meth. Eng. 10, 1211–1219.
Contents lists available at ScienceDirect
Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom
A non-linear corotational 4-node plane element Jean-Marc Battini Department of Mechanics, KTH, Royal Institute of Technology, SE-10044 Stockholm, Sweden
a r t i c l e
i n f o
Article history: Received 14 September 2007 Received in revised form 28 February 2008 Available online 16 March 2008
a b s t r a c t A non-linear corotational 4-node plane element is presented. The main interested of the approach is that accurate and efficient linear elements can be transformed to non-linear elements with little work and without considering the assumptions used to derive these linear elements. Two numerical examples are presented. Ó 2008 Elsevier Ltd. All rights reserved.
Keywords: Corotational formulation Plane elements Non-linear analyses
1. Introduction Due to their simplicity and their efficiency in large-scale analyses, low-order finite elements are often used in structural mechanics. For this reason, a lot of research has been done in order to develop accurate and efficient 4-node plane elements, see e.g. the paper of Fredriksson and Saabye Ottosen (2004) for a complete review. However, many of the proposed solutions concern only linear formulations. The reason is that these linear elements are by themselves complicated and their transformation to non-linear formulations is not always straightforward (see e.g. Simo and Armero, 1992; Simo et al., 1993). This is particularly the case if the template formulation proposed by Felippa (2006) is used. One simple way, at least for large displacements and small strains problems, to derive non-linear elements is to used the corotational approach. Corotational beams and shells elements have been widely studied, see e.g. the review article of Felippa and Haugen (2005). However, to the author’s knowledge, no complete work and no numerical analyses have been done with plane or solid elements. Besides, this method seems particularly interesting in the present context since it transforms linear elements to non-linear ones with very little work, and without considering the theory used to derive these linear elements. The main idea of the corotational approach is to decompose the motion of the element into rigid-body and pure deformational parts, through the use of a reference system which continuously rotates and translates with the element. The deformational response is captured at the level of the local reference frame, whereas the geometric non-linearity induced by the large rigid-body motion, is incorporated in the transformation matrices relating local and global internal force vectors and tangent stiffness matrices. Assuming the pure deformation part to be small, a geometrical linear theory can be used in the local system. The main interest is that for elements with the same number of nodes and degrees of freedom, the transformation matrices are the same. Using this property, different existing high performing geometrical linear elements can be used as local formulations. As mentioned above, this is particularly interesting for 4-node plane elements for which several advanced linear formulations are available in the literature.
E-mail address: [email protected] 0093-6413/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2008.03.002
J.-M. Battini / Mechanics Research Communications 35 (2008) 408–413
409
The corotational framework for 4-node plane elements is given in Section 2. Using for the local formulation two accurate linear elements found in the literature, two numerical problems presenting large displacements and small strains are presented in Section 3. The results obtained with the corotational approach are compared to the ones obtained with ANSYS. 2. Corotational formulation 2.1. Element kinematics The element kinematics is defined in Fig. 1. (O, X, Y) is the global coordinate system. ri ¼ ðX i ; Y i Þ and ui ¼ ðU i ; V i Þ denote the global coordinates and global displacements of the node i. The origin of the local coordinate system (C, x, y) is calculated by rc ¼ ðX c ; Y c Þ
Xc ¼
1 ðX 1 þ X 2 þ X 3 þ X 4 Þ 4
Yc ¼
1 ðY 1 þ Y 2 þ Y 3 þ Y 4 Þ 4
ð1Þ
In the initial configuration, the axes of the local and global coordinate systems are parallel. roi ¼ ðxi ; yi Þ defines the local coordinates of the node i. The motion of the element from the initial to the deformed configuration is split in two steps. The first one is a rigid translation and rotation of the undeformed element. The rigid translation is defined by uc , the global displacement of the point C, which is calculated by uc ¼ ðU c ; V c Þ
Uc ¼
1 ðU 1 þ U 2 þ U 3 þ U 4 Þ 4
Vc ¼
1 ðV 1 þ V 2 þ V 3 þ V 4 Þ 4
ð2Þ
The rigid rotation, specified by the angle h, defines the orientation of the local coordinate system in the deformed configuration. The second step is a local deformation with respect to the local frame. This step is defined by the local deformational dis i ¼ ðu i ; v i Þ; i ¼ 1; . . . ; 4 at the nodes. These displacements are calculated as placements u i xi u cos h sin h X i þ U i X c U c ¼ ð3Þ i yi v sin h cos h Y i þ V i Y c V c The rigid rotation h is taken such that the square of the Euclidean norm of the nodal local deformational displacements X4 2i 2 þ v ð4Þ u i¼1 i is minimized. Introducing Eq. (3) in Eq. (4), taking the variation with respect to h and setting the resulting expression equal to zero, gives P4 ½xi ðY i þ V i Y c V c Þ yi ðX i þ U i X c U c Þ tan h ¼ Pi¼1 ð5Þ 4 i¼1 ½xi ðX i þ U i X c U c Þ þ yi ðY i þ V i Y c V c Þ This equation gives two values for the rotation h, (h and p þ h). One value corresponds to the minimization of expression (4), the other value to its maximization. One simple way to choose between these values is to calculate the expression (4) with both angles and to take the angle which gives the smallest result. Finally, the local and global displacement vectors are 1 pl ¼ ½ u
1 v
2 u
2 v
3 u
3 v
4 u
4 T v
pg ¼ ½ U 1
V1
U2
V2
Fig. 1. Element kinematics.
U3
V3
U4
V 4 T
ð6Þ
410
J.-M. Battini / Mechanics Research Communications 35 (2008) 408–413
2.2. Transformation matrices The corotational framework is defined by a change of variables between local and global quantities. Several existing linear formulations can be taken for the local stiffness matrix Kl and internal force vector f l consistent with pl . The global tangent stiffness matrix Kg and internal force vector f g consistent with pg can be calculated by equating the virtual work in the two systems, as dpTg f g ¼ dpTl f l
ð7Þ
which, introducing the relations dpl ¼ Bdpg
df g ¼ Kg dpg
f l ¼ Kl pl
ð8Þ
gives f g ¼ BT f l
Kg ¼ BT Kl B þ Kh
oðBT f l Þ opg
Kh ¼
ð9Þ fl
Differentiation of Eq. (3) gives i ydi du cos h sin h dU i dU c ¼ þ dh i dv xdi sin h cos h dV i dV c
i þ xi xdi ¼ u i þ yi ydi ¼ v
The expression of dh can be obtained from Eq. (5) as 4 X cos h sin h dU i dU c 1 ½ x y dh ¼ i i 4 P sin h cos h dV i dV c i¼1 ðxi xdi þ yi ydi Þ
ð10Þ
ð11Þ
i¼1
From Eqs. (10) and (11), and omitting dU see Eq. (8), is expressed as B ¼ PET
c
and dV c since finite elements are invariant under pure translation, the matrix B,
P ¼ I AG
ð12Þ
with A ¼ ½ yd1 G¼
4 P
xd1 1
yd2
xd2
½ y1
yd3
x1
xd3
yd4
x2
y3
y2
xd4 x3
T
y4
ð13Þ x4
ð14Þ
ðxi xdi þ yi ydi Þ
i¼1
E ¼ diagðR; R; R; RÞ
R¼
cos h sin h
sin h cos h
ð15Þ
Using Eq. (12), the matrix Kh , see Eq. (9), is defined by Kh ¼ dEPT f l þ EdPT f l
ð16Þ T
After some work, and noting that A G
T
¼ I, it is obtained
Kh ¼ E½FT G GT FPET
ð17Þ
with PT f l ¼ ½ n1 F ¼ ½ n2
n2
n1
n3 n4
n4 n3
n5 n6
n6
n7
n5
n8 T n8
n7
ð18Þ ð19Þ
3. Numerical examples Two linear elements, Qm6 and Qnew, have been used for the local formulation. The Qm6 element has been introduced by Taylor et al. (1976), whereas the Qnew element has been developed by Fredriksson and Saabye Ottosen (2004). Two problems, presenting large displacements and small strains have been analysed. The results obtained with the present corotational formulation have been compared to the results obtained with ANSYS. The analyses in ANSYS have been performed using the element PLANE182 with KEYOPT(1) = 2. This element is a non-linear formulation of the Qm6 element based on the updated Lagrangian formulation. This element is implemented for large strains, but since the strains in the two test problems remain small, it can be used for the proposed comparison. For each structure, two meshes have been tested: a fine one which gives the converged solutions and a coarse one in order to test the ability of the present corotational formulation in obtaining good results with coarse meshes.
411
J.-M. Battini / Mechanics Research Communications 35 (2008) 408–413
3.1. Example 1: angle frame This example is described in Fig. 2. The structure is clamped at one end (all displacements zero) and loaded by a uniformly distributed horizontal force F at the other end. The fine mesh consists of 304 square elements whereas the coarse one consists of 19 distorted elements. The results obtained using 20 steps and the tolerance parameter 105 for the nodal forces
Fig. 2. Example 1: angle frame.
4.5 4 Qnew coarse mesh 3.5 Qm6 coarse mesh
Load F / 10000
3 ANSYS coarse mesh 2.5 fine mesh 2 1.5 1 0.5 0 0
1
2
3
4
5
6
7
Horizontal displacement of point A Fig. 3. Example 1: load–displacement diagram.
Table 1 Example 1 – total number of Newton–Raphson iterations
Coarse mesh Fine mesh
ANSYS
Qm6
Qnew
66 72
60 60
60 60
412
J.-M. Battini / Mechanics Research Communications 35 (2008) 408–413
Fig. 4. Example 2: thin ring.
11 10 9 8
Load F / 1000
7 6 Qnew coarse mesh 5 Qm6 coarse mesh 4 ANSYS coarse mesh 3 fine mesh 2 1 0 0
5
10
15
20
25
30
Vertical displacement of point A Fig. 5. Example 2: load–displacement diagram.
(command CNVTOL in ANSYS) are shown in Fig. 3. With the fine mesh, the three elements give exactly the same results. With the coarse mesh, ANSYS and Qm6 give almost the same results whereas Qnew gives slightly better results. The total number of Newton–Raphson iterations are given in Table 1. It can be observed that the convergence rate of the corotational approach is slightly better than the one in ANSYS. 3.2. Example 2: thin ring A half ring, clamped at both ends is loaded at its top by a concentrated force 2F. Due to the symmetry, only half of the structure is modelled, with 80 4 ¼ 320 elements and 20 1 ¼ 20 elements, see Fig. 4. The results are shown in Fig. 5. The three elements give the same results, both with the coarse and fine meshes. It can also be noted that the results obtained with the coarse mesh are very closed to the ones obtained with the fine mesh. 4. Conclusion This paper has described a corotational framework for 4-node plane elements. The linear elements Qm6 and Qnew, found in the literature, has been used as local formulations. Two numerical applications have shown that the proposed formulation
J.-M. Battini / Mechanics Research Communications 35 (2008) 408–413
413
can be used to analyse plane problems presenting large displacements and small strains. Very good numerical results have been obtained, both with fine and coarse meshes, and also with distorted elements. The main interest of the proposed approach is that, using the described corotational framework, efficient linear elements are automatically transformed to nonlinear formulations. References Felippa, C., 2006. Supernatural QUAD4: A template formulation. Comput. Meth. Appl. Mech. Eng. 195, 5316–5342. Felippa, C., Haugen, B., 2005. A unified formulation of small-strain corotational finite elements: I. Theory. Comput. Meth. Appl. Mech. Eng. 194, 2285–2335. Fredriksson, M., Saabye Ottosen, N., 2004. Fast and accurate 4-node quadrilateral. Int. J. Numer. Meth. Eng. 61, 1809–1834. Simo, J., Armero, F., 1992. Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes. Int. J. Numer. Meth. Eng. 33, 1413–1449. Simo, J., Armero, F., Taylor, R., 1993. Improved versions of assumed enhanced strain tri-linear elements for 3d finite deformation problems. Comput. Meth. Appl. Mech. Eng. 110, 359–386. Taylor, R., Beresford, P., Wilson, E., 1976. A non-conforming element for stress analysis. Int. J. Numer. Meth. Eng. 10, 1211–1219.