(67 intermediate revisions by the same user not shown) | |||
Line 4: | Line 4: | ||
The finite formulation of the 4-node thin shell element based on Co-rotational (CR) and Timoshenko's theory was presented for the analysis of laminated composite structures. Based on the Timoshenko's theory the present thin shell element avoids shear locking behavior, and a bilinear in-plane displacement field is introduced for the coupling of in-plane and bending actions, and such element performs simple formulation and highly efficienct. A highly efficient CR formulation was also established to define the motion of the element with large displacement. The characteristics of the element are more pronounced when shells get thin, each node of the 4-node element has five degrees of freedom. An incremental iterative method based on the Newton Raphson method was used to solve the nonlinear equilibrium equations. A number of numerical examples were given to verify that the formula is computationally efficient and the results showed good agreement. | The finite formulation of the 4-node thin shell element based on Co-rotational (CR) and Timoshenko's theory was presented for the analysis of laminated composite structures. Based on the Timoshenko's theory the present thin shell element avoids shear locking behavior, and a bilinear in-plane displacement field is introduced for the coupling of in-plane and bending actions, and such element performs simple formulation and highly efficienct. A highly efficient CR formulation was also established to define the motion of the element with large displacement. The characteristics of the element are more pronounced when shells get thin, each node of the 4-node element has five degrees of freedom. An incremental iterative method based on the Newton Raphson method was used to solve the nonlinear equilibrium equations. A number of numerical examples were given to verify that the formula is computationally efficient and the results showed good agreement. | ||
− | + | ==1. Introduction== | |
− | + | ||
− | ==1. | + | |
Many nonlinear shell element formulations are based on the total Lagrangian (TL) and update Lagrangian (UL), including complex formulations and increasing calculation interval. However, in order to maximum the computational efficiency, a nonlinear finite element method (FEM) is required. A named “CR” approach was introduced by Wempner [1] and Belytschko [2] firstly. | Many nonlinear shell element formulations are based on the total Lagrangian (TL) and update Lagrangian (UL), including complex formulations and increasing calculation interval. However, in order to maximum the computational efficiency, a nonlinear finite element method (FEM) is required. A named “CR” approach was introduced by Wempner [1] and Belytschko [2] firstly. | ||
Line 18: | Line 16: | ||
The purpose of the paper is to improve the computational efficiency of the thin shell structures modelling by a high efficiency nonlinear 4-node thin shell element. In the local element formulation, a simple and high efficiency laminated composite element was presented et al. [20], the tangential rotation and the shear strain were obtained by Timoshenko laminated composite beam theory, and a bilinear in-plane displacement field was introduced for the coupling of in-plane and bending actions. In order to consider the movement of the element with large translation and rotation, a highly efficienct CR method was obtained. In the current work, the main idea introduced by Pacoste [9] and Battini [11] was employed. At last, a high efficiency tangent stiffness matrix of 4-node thin shell element was obtained in this paper. | The purpose of the paper is to improve the computational efficiency of the thin shell structures modelling by a high efficiency nonlinear 4-node thin shell element. In the local element formulation, a simple and high efficiency laminated composite element was presented et al. [20], the tangential rotation and the shear strain were obtained by Timoshenko laminated composite beam theory, and a bilinear in-plane displacement field was introduced for the coupling of in-plane and bending actions. In order to consider the movement of the element with large translation and rotation, a highly efficienct CR method was obtained. In the current work, the main idea introduced by Pacoste [9] and Battini [11] was employed. At last, a high efficiency tangent stiffness matrix of 4-node thin shell element was obtained in this paper. | ||
− | ==2. | + | ==2. Local deformation== |
In local coordinate, Cen [20] show that tangential rotation and shear strain were obtained by Timoshenko beam theory. Shear strain and rotation fields within the element are then determined by improved rational interpolation, and a bilinear displacement field is introduced for the coupling of in-plane and bending actions. Finally, in order to analyze random laminated composite structures, the 20 degrees of freedom quadrilateral bending element can be used. | In local coordinate, Cen [20] show that tangential rotation and shear strain were obtained by Timoshenko beam theory. Shear strain and rotation fields within the element are then determined by improved rational interpolation, and a bilinear displacement field is introduced for the coupling of in-plane and bending actions. Finally, in order to analyze random laminated composite structures, the 20 degrees of freedom quadrilateral bending element can be used. | ||
− | + | ===2.1 Element shear strain and in-plane strain field of the mid-plane=== | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | ==2.1 | + | |
As shown in Figure1, the local deformational displacements and rotations for random 4-node elements can be written as | As shown in Figure1, the local deformational displacements and rotations for random 4-node elements can be written as | ||
Line 36: | Line 29: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | [ | + | | <math> \bar{a}_{i} =[\begin{array}{ccccc} {u_{i} } & {v_{i} } & {w_{i} } & {\phi _{xi} } & {\phi _{yi} } \end{array}]^{T} \quad \quad \quad i=1,2,3,4 </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (1) | | style="width: 5px;text-align: right;white-space: nowrap;" | (1) | ||
|} | |} | ||
− | based on the local deformation of the element mid-plane, displacement and rotation of random point for the | + | |
+ | <div id='img-1'></div> | ||
+ | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | ||
+ | |- | ||
+ | | [[Image:draft_Aparicio Nogué_939719944-image2.png|534px]] | ||
+ | |- style="text-align: center; font-size: 75%;" | ||
+ | | colspan="1" | '''Figure 1:''' Local displacement and rotation at the mid-plane of a laminated composite element. | ||
+ | |} | ||
+ | |||
+ | based on the local deformation of the element mid-plane, displacement and rotation of random point for the elmenet can be written as | ||
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 47: | Line 49: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>\begin{array}{l} {u(x,y,z)=u_{0} (x,y)-z\phi _{x} (x,y)} \\ {v(x,y,z)=v_{0} (x,y)-z\phi _{y} (x,y)} \\ {w(x,y,z)=w(x,y)} \end{array} |
+ | </math> | ||
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (2) | | style="width: 5px;text-align: right;white-space: nowrap;" | (2) | ||
Line 59: | Line 62: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>\varphi _{s} =\varphi _{si} (1-r)+\varphi _{sj} r+3(1-2\delta )r(1-r)\Gamma</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (3) | | style="width: 5px;text-align: right;white-space: nowrap;" | (3) | ||
Line 69: | Line 72: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>\gamma _{si} =\delta _{i} \Gamma</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (4) | | style="width: 5px;text-align: right;white-space: nowrap;" | (4) | ||
Line 75: | Line 78: | ||
with | with | ||
− | |||
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
|- | |- | ||
Line 89: | Line 91: | ||
|} | |} | ||
− | <div | + | <div id='img-2'></div> |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | |
− | + | |- | |
− | + | | [[Image:draft_Aparicio Nogué_939719944-image7.png|480px]] | |
− | + | |- style="text-align: center; font-size: 75%;" | |
+ | | colspan="1" | '''Figure 2:''' Laminated composite beam element. | ||
+ | |} | ||
Based on the Equation (4), element shear strain of every edge can be written as | Based on the Equation (4), element shear strain of every edge can be written as | ||
Line 119: | Line 123: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math> \delta _{i} =\frac{6\lambda _{i} }{1+12\lambda _{i} } ,\quad \quad \lambda _{i} =\frac{D_{di} }{C_{di} d_{i}^{2} } \quad \quad \quad (i=1,2,3,4)</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (7) | | style="width: 5px;text-align: right;white-space: nowrap;" | (7) | ||
|} | |} | ||
− | where | + | where <math>d_{i} (i=1,2,3,4)</math> is the length of element edge, <math>b_{i} (i=1,2,3,4)</math> is the y orientation distance between two adjacent nodes, <math>c_{i} (i=1,2,3,4)</math> is the <math>x</math> orientation distance between two adjacent nodes. <math>D_{di}</math>, <math>C_{di}</math> are the shear and bending stiffness. |
− | According to the element shear strain of every edge | + | According to the element shear strain of every edge <math>\gamma _{si} (i=1,\; 2,\; 3,\; 4)</math>, <math>\gamma _{s}^{*}</math> can be written as |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 133: | Line 137: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | [ | + | | <math>\gamma _{s}^{*} =\left[\begin{array}{cccc} {d_{1} \gamma _{s1} } & {d_{2} \gamma _{s2} } & {d_{3} \gamma _{s3} } & {d_{4} \gamma _{s4} } \end{array}\right] </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (8) | | style="width: 5px;text-align: right;white-space: nowrap;" | (8) | ||
Line 145: | Line 149: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | [ | + | | <math>\left\{\begin{array}{c} {\gamma _{x1} } \\ {\gamma _{y1} } \end{array}\right\}=\frac{1}{b_{3} c_{4} -b_{4} c_{3} } \left[\begin{array}{cc} {b_{3} } & {-b_{4} } \\ {c_{3} } & {-c_{4} } \end{array}\right]\left\{\begin{array}{c} {\gamma _{s4}^{*} } \\ {\gamma _{s3}^{*} } \end{array}\right\}</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (9) | | style="width: 5px;text-align: right;white-space: nowrap;" | (9) | ||
Line 157: | Line 161: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | [ | + | | <math>\gamma =\left\{\begin{array}{c} {\gamma _{xz} } \\ {\gamma _{yz} } \end{array}\right\}=\left[\begin{array}{c} {\gamma _{x1} N_{1}^{0} +\gamma _{x2} N_{2}^{0} +\gamma _{x3} N_{3}^{0} +\gamma _{x4} N_{4}^{0} } \\ {\gamma _{y1} N_{1}^{0} +\gamma _{y2} N_{2}^{0} +\gamma _{y3} N_{3}^{0} +\gamma _{y4} N_{4}^{0} } \end{array}\right]=B_{s} \bar{a}</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (10) | | style="width: 5px;text-align: right;white-space: nowrap;" | (10) | ||
Line 177: | Line 181: | ||
|} | |} | ||
− | ==2.2 | + | ===2.2 Rotation and in-plane displacement field=== |
The element normal and tangential of the first edge is shown in Figure 1, each edge of the element is assumed to be linearly distributed. | The element normal and tangential of the first edge is shown in Figure 1, each edge of the element is assumed to be linearly distributed. | ||
Line 226: | Line 230: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>\phi _{x} =\sum _{i=1}^{8}N_{i} \phi _{xi} ,\quad \quad \phi _{y} =\sum _{i=1}^{8}N_{i} \phi _{yi}</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (14) | | style="width: 5px;text-align: right;white-space: nowrap;" | (14) | ||
Line 248: | Line 252: | ||
|} | |} | ||
− | The bilinear in-plane displacement | + | The bilinear in-plane displacement <math>u^{0}</math> and <math>v^{0}</math> can be written as |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 255: | Line 259: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>u^{0} =\sum _{i=1}^{4}N_{i}^{0} u_{i} ,\quad \quad v^{0} =\sum _{i=1}^{4}N_{i}^{0} v_{i}</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (16) | | style="width: 5px;text-align: right;white-space: nowrap;" | (16) | ||
|} | |} | ||
− | where | + | where <math>N_{i}^{0} (i=1,2,3,4) </math>are the quadrilateral element shape functions. |
− | ==2.3 | + | ===2.3 Element in-plane strain and curvature field of the mid-plane=== |
− | + | Element in-plane strain can be written as | |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 271: | Line 275: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | |<math> \varepsilon ^{0} =B_{e} \bar{a}</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (17) | | style="width: 5px;text-align: right;white-space: nowrap;" | (17) | ||
|} | |} | ||
− | + | with | |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 283: | Line 287: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>\begin{array}{ | + | | <math>\begin{array}{l} |
B_e=[\begin{array}{cccc} | B_e=[\begin{array}{cccc} | ||
B_{e1} & B_{e2} & B_{e3} & B_{e4} | B_{e1} & B_{e2} & B_{e3} & B_{e4} | ||
Line 291: | Line 295: | ||
0 & \frac{\partial N_i^0}{\partial y} & 0 & 0 & 0\\ | 0 & \frac{\partial N_i^0}{\partial y} & 0 & 0 & 0\\ | ||
\frac{\partial N_i^0}{\partial y} & \frac{\partial N_i^0}{\partial x} & 0 & 0 & 0 | \frac{\partial N_i^0}{\partial y} & \frac{\partial N_i^0}{\partial x} & 0 & 0 & 0 | ||
− | \end{array}\right] | + | \end{array}\right]\\ |
− | a | + | \overline\mathit{\boldsymbol{a}}=\left[\begin{array}{cccc} |
− | + | \overline\mathit{\boldsymbol{a}}_1 & \overline\mathit{\boldsymbol{a}}_2 & \overline\mathit{\boldsymbol{a}}_3 & \overline\mathit{\boldsymbol{a}}_4 | |
\end{array}\right] | \end{array}\right] | ||
− | \end{array}</math> | + | \end{array}\qquad (i=1,2,3,4)</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (18) | | style="width: 5px;text-align: right;white-space: nowrap;" | (18) | ||
Line 307: | Line 311: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>\kappa =B_{b} \bar{a}</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (19) | | style="width: 5px;text-align: right;white-space: nowrap;" | (19) | ||
|} | |} | ||
− | where | + | where <math>B_{b}</math> is bendding strain matrix. |
Based on the constitutive equation of laminated composite plate, the relationship between strain and displacement in mid-plane can be written as | Based on the constitutive equation of laminated composite plate, the relationship between strain and displacement in mid-plane can be written as | ||
Line 321: | Line 325: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>{\ | + | | <math>\mathit{\boldsymbol{\varepsilon }}_p=\left\{\begin{array}{c} |
− | {\ | + | \mathit{\boldsymbol{\varepsilon }}^0\\ |
− | \kappa | + | \mathit{\boldsymbol{\kappa}} |
\end{array}\right\}=\left\{\begin{array}{c} | \end{array}\right\}=\left\{\begin{array}{c} | ||
B_e\\ | B_e\\ | ||
B_b | B_b | ||
− | \end{array}\right\}a | + | \end{array}\right\}\overline\mathit{\boldsymbol{a}}=B_p\overline\mathit{\boldsymbol{a}}</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (20) | | style="width: 5px;text-align: right;white-space: nowrap;" | (20) | ||
|} | |} | ||
− | ==2.4 | + | ===2.4 Element stiffness matrix=== |
The stiffness matrix of the 4-node element is derived from the principle of minimum potential energy, element shear strain and in-plane strain field: | The stiffness matrix of the 4-node element is derived from the principle of minimum potential energy, element shear strain and in-plane strain field: | ||
Line 341: | Line 345: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>K_{l}^{e} =\int _{-1}^{1}\int _{-1}^{1}B_{P}^{T} C_{P} B_{P} \left|J\right|d\xi d\eta + \int _{-1}^{1}\int _{-1}^{1}B_{S}^{T} C_{S} B_{S} \left|J\right|d\xi d\eta</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (21) | | style="width: 5px;text-align: right;white-space: nowrap;" | (21) | ||
|} | |} | ||
− | where shear strain matrix | + | where shear strain matrix <math> B_{S}</math> and in-plane strain matrix <math>B_{P} </math> of the mid-plane can be found in Equations (10) and (20), respectively, <math>C_p</math> is stress-strain matrix of laminated composite plate, <math>C_s</math> is shear stiffness, and <math>\left|J\right| </math> is the Jacobian determinant. The element stiffness matrix was reduced to a <math>20\times 20 </math> matrix, which was calculated by <math>3\times 3 </math> Gauss integral. |
− | + | ==3. Co-rotational formulation== | |
− | + | ||
− | ==3. | + | |
CR approach is more efficient than TL and UL to analyze thin wall structure models. Based on CR description, a lot of plate and shell elements including 3 node and 4 node were established and confirmed by many numerical illustrations and tests. According to CR description, element kinematics from initial undeformed to final deformed attitude can be divided into two sections [7]. The first section is rigid rotation removing deformational part.and translation. The second step is a local deformation and rotation under local coordinate system. Element deformation displacement and rotation were established conveniently in section 2. | CR approach is more efficient than TL and UL to analyze thin wall structure models. Based on CR description, a lot of plate and shell elements including 3 node and 4 node were established and confirmed by many numerical illustrations and tests. According to CR description, element kinematics from initial undeformed to final deformed attitude can be divided into two sections [7]. The first section is rigid rotation removing deformational part.and translation. The second step is a local deformation and rotation under local coordinate system. Element deformation displacement and rotation were established conveniently in section 2. | ||
− | ==3.1 Parameterization of the rigid rotations and coordinate systems== | + | ===3.1 Parameterization of the rigid rotations and coordinate systems=== |
− | In order to represent large three-dimensional rotations, a parameterization approach of the rigid rotations is introduced [9]. Orthogonal matrix <math>R</math> can be obtained by three separate parameter variables. Firstly, the rotational vector can be defined by | + | In order to represent large three-dimensional rotations, a parameterization approach of the rigid rotations is introduced [9]. Orthogonal matrix <math>\boldsymbol R</math> can be obtained by three separate parameter variables. Firstly, the rotational vector can be defined by |
− | + | {| class="formulaSCP" style="width: 100%; text-align: left;" | |
+ | |- | ||
+ | | | ||
+ | {| style="text-align: center; margin:auto;width: 100%;" | ||
+ | |- | ||
+ | | style="text-align: center;" | <math>\boldsymbol{\Psi}= \boldsymbol{e}\psi </math> | ||
+ | | style="width: 5px;text-align: right;white-space: nowrap;" |(22) | ||
+ | |} | ||
+ | |} | ||
− | According to the definition above, as shown in Figure 3, rotation can be obtained by an angle | + | According to the definition above, as shown in Figure 3, rotation can be obtained by an angle <math>\Psi</math> defined by the unit vector <math>e</math>. So the value of <math>\Psi</math> is obtained by |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 372: | Line 382: | ||
|} | |} | ||
− | where <math>{\Psi }_i</math> ( <math>i=1,2,3</math> ) are the components of <math>\Psi </math> . | + | where <math>{\Psi }_i</math> (<math>i=1,2,3</math>) are the components of <math>\Psi </math>. |
− | <div | + | <div id='img-3'></div> |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | |
− | + | |- | |
− | + | | [[Image:draft_Aparicio Nogué_939719944-image49.png|342px]] | |
− | + | |- style="text-align: center; font-size: 75%;" | |
+ | | colspan="1" | '''Figure 3:''' Rotational vector. | ||
+ | |} | ||
In terms of the definition above, the orthogonal matrix <math>R</math> can be written as | In terms of the definition above, the orthogonal matrix <math>R</math> can be written as | ||
Line 387: | Line 399: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>R=I_3+\frac{sin\psi }{\psi }\Psi +\frac{1}{2}{\left(\frac{sin(\psi /2)}{\psi /2}\right)}^2{\Psi }^2</math> | + | | <math>\boldsymbol R =\boldsymbol I_3+\frac{sin\psi }{\psi }\Psi +\frac{1}{2}{\left(\frac{sin(\psi /2)}{\psi /2}\right)}^2{\Psi }^2</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (24) | | style="width: 5px;text-align: right;white-space: nowrap;" | (24) | ||
|} | |} | ||
− | Based on the conception of the rigid rotation, the global rotation matrix <math>R_i^g</math> at point <math>i</math> can be defined by <math>R_r</math> . | + | Based on the conception of the rigid rotation, the global rotation matrix <math>\boldsymbol R_i^g</math> at point <math>i</math> can be defined by <math>\boldsymbol R_r</math> . |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 399: | Line 411: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>R_i^g=I+\frac{sin\psi }{\psi }{\Psi }_i^g+\frac{1}{2}{\left(\frac{sin(\psi /2)}{\psi /2}\right)}^2{\Psi }_i^g{}^2</math> , <math>\psi =\Vert {\Psi }_i^g\Vert </math> | + | | <math>\boldsymbol R_i^g=\boldsymbol I+\frac{sin\psi }{\psi }{\Psi }_i^g+\frac{1}{2}{\left(\frac{sin(\psi /2)}{\psi /2}\right)}^2{\Psi }_i^g{}^2</math> , <math>\psi =\Vert {\Psi }_i^g\Vert </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (25) | | style="width: 5px;text-align: right;white-space: nowrap;" | (25) | ||
|} | |} | ||
− | <div | + | <div id='img-4'></div> |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | |
− | + | |- | |
− | + | | [[Image:draft_Aparicio Nogué_939719944-image57.png|600px]] | |
− | + | |- style="text-align: center; font-size: 75%;" | |
+ | | colspan="1" | '''Figure 4:''' Element kinematics and coordinate systems. | ||
+ | |} | ||
− | As shown in Figure 4, we define <math>e_i</math> as the vectors of the local frame in the current deformed configuration. Then, the matrix <math>R_r</math> [9] defining rotations can be obtained as | + | As shown in Figure 4, we define <math>\boldsymbol e_i</math> as the vectors of the local frame in the current deformed configuration. Then, the matrix <math>\boldsymbol R_r</math> [9] defining rotations can be obtained as |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 417: | Line 431: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>R_r=\left[\begin{array}{ccc} | + | | <math>\boldsymbol R_r=\left[\begin{array}{ccc} |
− | e_1 & e_2 & e_3 | + | \boldsymbol e_1 & \boldsymbol e_2 & \boldsymbol e_3 |
\end{array}\right]</math> | \end{array}\right]</math> | ||
|} | |} | ||
Line 432: | Line 446: | ||
|- | |- | ||
| <math>\begin{array}{c} | | <math>\begin{array}{c} | ||
− | e_1=\frac{\left(r_j^g+u_j^g-r_i^g-u_i^g\right)}{\Vert r_j^g+u_j^g-r_i^g-u_i^g\Vert }\\ | + | \boldsymbol e_1=\frac{\left(\boldsymbol r_j^g+\boldsymbol u_j^g-\boldsymbol r_i^g-\boldsymbol u_i^g\right)}{\Vert \boldsymbol r_j^g+\boldsymbol u_j^g-\boldsymbol r_i^g-\boldsymbol u_i^g\Vert }\\ |
− | e_3=\frac{x_j\times x_k}{\Vert x_j\times x_k\Vert }\\ | + | \boldsymbol e_3=\frac{\boldsymbol x_j\times \boldsymbol x_k}{\Vert \boldsymbol x_j\times \boldsymbol x_k\Vert }\\ |
− | e_2=e_1\times e_3\\ | + | \boldsymbol e_2=\boldsymbol e_1\times \boldsymbol e_3\\ |
− | x_j=r_j^g+u_j^g-r_i^g-u_i^g\\ | + | \boldsymbol x_j=\boldsymbol r_j^g+\boldsymbol u_j^g-\boldsymbol r_i^g-\boldsymbol u_i^g\\ |
− | x_k=r_k^g+u_k^g-r_i^g-u_i^g | + | \boldsymbol x_k=\boldsymbol r_k^g+\boldsymbol u_k^g-\boldsymbol r_i^g-\boldsymbol u_i^g |
\end{array}</math> | \end{array}</math> | ||
|} | |} | ||
Line 442: | Line 456: | ||
|} | |} | ||
− | where <math>u_i^g</math> and <math>u_j^g</math> denote displacement from the initial to current coordinate, respectively. <math>r_i^g</math> and <math>r_j^g</math> denote displacement vector of initial nodes <math>i</math> and <math>j</math> in the global system, respectively. <math>x_j</math> and <math>x_k</math> denote displacement vector along the orientations <math>j</math> and <math>k</math> of the line connecting two nodes, respectively. | + | where <math>\boldsymbol u_i^g</math> and <math>\boldsymbol u_j^g</math> denote displacement from the initial to current coordinate, respectively. <math>\boldsymbol r_i^g</math> and <math>\boldsymbol r_j^g</math> denote displacement vector of initial nodes <math>i</math> and <math>j</math> in the global system, respectively. <math>\boldsymbol x_j</math> and <math>\boldsymbol x_k</math> denote displacement vector along the orientations <math>j</math> and <math>k</math> of the line connecting two nodes, respectively. |
− | + | ||
− | + | ||
+ | The local deformation of point <math>i</math> can be expressed by <math>\boldsymbol {\overline{a}}_i</math> , as shown in section 2. According to the define of rigid frame in the local coordinate, <math>\boldsymbol {\overline{a}}_i</math> can be redefined as | ||
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
|- | |- | ||
Line 451: | Line 464: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>{\overline{a}}_i=R_r^T(r_i^g+u_i^g-r_c^g-u_c^g)-r_c^1</math> | + | | <math>\boldsymbol {\overline{a}}_i=\boldsymbol R_r^T(\boldsymbol r_i^g+\boldsymbol u_i^g-\boldsymbol r_c^g-\boldsymbol u_c^g)-\boldsymbol r_c^1</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (28) | | style="width: 5px;text-align: right;white-space: nowrap;" | (28) | ||
|} | |} | ||
− | As shown in Figure 4, noting that <math>R_r{\overline{R}}_i=R_i^ | + | As shown in Figure 4, noting that <math>\boldsymbol R_r\boldsymbol {\overline{R}}_i=\boldsymbol R_i^g\boldsymbol R_0</math> , the local rotations are defined by the matrices <math>\boldsymbol {\overline{R}}_i</math> <math>\left(i=1,2,3,4\right)</math> |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
as | as | ||
Line 468: | Line 478: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>{\overline{R}}_i=R_r | + | | <math>\boldsymbol {\overline{R}}_i=\boldsymbol R_r^T \boldsymbol R_i^g \boldsymbol R_0</math> , <math>\left(i=1,2,3,4\right).</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (29) | | style="width: 5px;text-align: right;white-space: nowrap;" | (29) | ||
Line 480: | Line 490: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>\delta {\overline{R}}_i=\delta {\tilde{\overline{\varphi }}}_i{\overline{R}}_i,\ | + | | <math>\delta \boldsymbol {\overline{R}}_i=\delta {\tilde{\overline{\varphi }}}_i \boldsymbol {\overline{R}}_i,\delta \boldsymbol R_i^g=</math><math>\delta {\tilde{\varphi }}_i^g \boldsymbol R_i^g,\delta \boldsymbol R_r=</math><math>\delta {\tilde{\varphi }}_r^g \boldsymbol R_r.</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (30) | | style="width: 5px;text-align: right;white-space: nowrap;" | (30) | ||
|} | |} | ||
− | where <math>\delta {\tilde{\overline{\varphi }}}_i</math> , <math>\delta {\tilde{\varphi }}_i^g</math> , <math>\delta {\tilde{\varphi }}_r^g</math> and the corresponding associated vectors <math>\delta {\overline{\varphi }}_i</math> , <math>\delta {\varphi }_i^g</math> , <math>\delta {\varphi }_r^g</math> denote spatial angular variations, which superimposed onto orthogonal matrices <math>{\overline{R}}_i</math> , <math>R_i^g</math> and <math>R_r</math> . | + | where <math>\delta {\tilde{\overline{\varphi }}}_i</math> , <math>\delta {\tilde{\varphi }}_i^g</math> , <math>\delta {\tilde{\varphi }}_r^g</math> and the corresponding associated vectors <math>\delta {\overline{\varphi }}_i</math> , <math>\delta {\varphi }_i^g</math> , <math>\delta {\varphi }_r^g</math> denote spatial angular variations, which superimposed onto orthogonal matrices <math>\boldsymbol {\overline{R}}_i</math> , <math>\boldsymbol R_i^g</math> and <math>\boldsymbol R_r</math> . |
Using Equations (28), (29) and (30), the differentiation of Equations (28) can be given as | Using Equations (28), (29) and (30), the differentiation of Equations (28) can be given as | ||
Line 495: | Line 505: | ||
|- | |- | ||
| <math>\begin{array}{c} | | <math>\begin{array}{c} | ||
− | \delta {\overline{a}}_i=R_r^T\delta u_i^g-R_r^T\delta {\tilde{\varphi }}_r^g\left(r_i^g+u_i^g-r_c^g-u_c^g\right)\\ | + | \delta \boldsymbol {\overline{a}}_i=\boldsymbol R_r^T\delta \boldsymbol u_i^g-\boldsymbol R_r^T\delta {\tilde{\varphi }}_r^g\left(\boldsymbol r_i^g+\boldsymbol u_i^g-\boldsymbol r_c^g-\boldsymbol u_c^g\right)\\ |
− | =\delta u_i^e-\delta {\tilde{\varphi }}_r^e\left({\overline{a}}_i+r_c^1\right)=\delta u_i^e+s_i\delta {\tilde{\varphi }}_r^e | + | =\delta \boldsymbol u_i^e-\delta {\tilde{\varphi }}_r^e\left(\boldsymbol {\overline{a}}_i+\boldsymbol r_c^1\right)=\delta \boldsymbol u_i^e+\boldsymbol s_i\delta {\tilde{\varphi }}_r^e |
\end{array}</math> | \end{array}</math> | ||
|} | |} | ||
Line 509: | Line 519: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>\boldsymbol s_i =\boldsymbol {\overline{a}}_i+\boldsymbol r_c^1</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (32) | | style="width: 5px;text-align: right;white-space: nowrap;" | (32) | ||
Line 521: | Line 531: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math> \begin{array}{ll} \delta\boldsymbol {\overline{R}} |
+ | & = \boldsymbol R_{r}^{T} \delta \boldsymbol R_{i}^{g} \boldsymbol R_{0} +\delta \boldsymbol R_{r}^{T} \boldsymbol R_{i}^{g} \boldsymbol R_{0} \\ | ||
+ | & =\delta \boldsymbol\tilde{\phi }_{i}^{e} \boldsymbol R_{r}^{T} \boldsymbol R_{i}^{g} \boldsymbol R_{0} -\delta \boldsymbol\tilde{\phi }_{r}^{e} \boldsymbol R_{r}^{T} \boldsymbol R_{i}^{g} \boldsymbol R_{0} =\left(\delta \boldsymbol\tilde{\phi }_{i}^{e} -\delta \boldsymbol\tilde{\phi }_{r}^{e} \right)\boldsymbol\bar{R}_{i} \end{array} </math> | ||
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (33) | | style="width: 5px;text-align: right;white-space: nowrap;" | (33) | ||
Line 533: | Line 545: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>\delta \boldsymbol\bar{\phi }_{i} = \delta \boldsymbol{\phi }_{i}^{e}- \boldsymbol{\phi }_{r}^{e} </math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (34) | | style="width: 5px;text-align: right;white-space: nowrap;" | (34) | ||
|} | |} | ||
− | ==3.2 Transformation matrix== | + | ===3.2 Transformation matrix=== |
The transformation matrix is needed to convert the internal force and tangent stiffness matrix. In the global coordinate, the 4-node element displacement is defined as | The transformation matrix is needed to convert the internal force and tangent stiffness matrix. In the global coordinate, the 4-node element displacement is defined as | ||
Line 547: | Line 559: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>\delta a_g={\left[\begin{array}{cccccccc} | + | | <math>\delta \boldsymbol a_g={\left[\begin{array}{cccccccc} |
− | \delta u_1^g{}^T & \delta {\psi }_1^g{}^T & \delta u_2^g{}^T & \delta {\psi }_2^g{}^T & \delta u_3^g{}^T & \delta {\psi }_3^g{}^T & \delta u_4^g{}^{} & \delta {\psi }_4^g{}^T | + | \delta \boldsymbol u_1^g{}^T & \delta {\psi }_1^g{}^T & \delta \boldsymbol u_2^g{}^T & \delta {\psi }_2^g{}^T & \delta \boldsymbol u_3^g{}^T & \delta {\psi }_3^g{}^T & \delta \boldsymbol u_4^g{}^{} & \delta {\psi }_4^g{}^T |
\end{array}\right]}^T</math> | \end{array}\right]}^T</math> | ||
|} | |} | ||
Line 554: | Line 566: | ||
|} | |} | ||
− | where <math>u_i^g\left(i=1, | + | where <math>\boldsymbol u_i^g\left(i=1,2,3,4\right)</math> are the global displacement vectors of the nodes, <math>\delta {\varphi }_i^g\left(i=1,\mbox{ }2,\mbox{ }3,\mbox{ }4\right)</math> denotes the spatial angular variations, defined in Equation (30). |
− | Referring to Equations (23), (25) and (30), the transformation matrix <math>T_m</math> was defined by | + | Referring to Equations (23), (25) and (30), the transformation matrix <math>\boldsymbol T_m</math> was defined by |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 563: | Line 575: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>\delta {\varphi }_i^g=T_m({\psi }_i^g)\delta {\psi }_i^g</math> | + | | <math>\delta {\varphi }_i^g=\boldsymbol T_m({\psi }_i^g)\delta {\psi }_i^g</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (36) | | style="width: 5px;text-align: right;white-space: nowrap;" | (36) | ||
Line 575: | Line 587: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>T_m({\psi }_i^g)=\left(R_i^g+I_3\right)/\sqrt{1-{\psi }_1^2-{\psi }_2^2-{\psi }_3^2}</math> | + | | <math>\boldsymbol T_m({\psi }_i^g)=\left(\boldsymbol R_i^g+\boldsymbol I_3\right)/\sqrt{1-{\psi }_1^2-{\psi }_2^2-{\psi }_3^2}</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (37) | | style="width: 5px;text-align: right;white-space: nowrap;" | (37) | ||
|} | |} | ||
− | where <math>{\psi }_i\ | + | where <math>{\psi }_i\,\left(i=1,2,3\right)</math> are the components of <math>\psi </math> , <math>\boldsymbol I_3</math> is the <math> 3\times 3</math> unit diagonal matrix. |
Using Equations (31), (32) and (34), the matrix can be written as | Using Equations (31), (32) and (34), the matrix can be written as | ||
Line 589: | Line 601: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>\delta \overline{a}=PE^T\delta a_g</math> | + | | <math>\delta \overline{\boldsymbol a}=\boldsymbol{PE}^T\delta \boldsymbol a_g</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (38) | | style="width: 5px;text-align: right;white-space: nowrap;" | (38) | ||
|} | |} | ||
− | The | + | The <math> 5\times 6</math> blocks <math> \boldsymbol P</math> is the matrix removing deformational part: |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 602: | Line 614: | ||
|- | |- | ||
| <math>\begin{array}{c} | | <math>\begin{array}{c} | ||
− | P_{ij}=\left[\begin{array}{cc} | + | \boldsymbol P_{ij}=\left[\begin{array}{cc} |
− | \frac{\partial {\overline{u}}_i^e}{\partial u_j^e} & \frac{\partial {\overline{u}}_i^e}{\partial {\varphi }_j^e}\\ | + | \frac{\partial {\overline{\boldsymbol u}}_i^e}{\partial \boldsymbol u_j^e} & \frac{\partial {\overline{\boldsymbol u}}_i^e}{\partial {\varphi }_j^e}\\ |
− | \frac{\partial {\varphi }_i^e}{\partial u_j^e} & \frac{\partial {\varphi }_i^e}{\partial {\varphi }_j^e} | + | \frac{\partial {\varphi }_i^e}{\partial \boldsymbol u_j^e} & \frac{\partial {\varphi }_i^e}{\partial {\varphi }_j^e} |
− | \end{array}\right]=I_{56}{\delta }_{\mbox{ij}}- | + | \end{array}\right]=\boldsymbol I_{56}{\delta }_{\mbox{ij}}-\boldsymbol A_i\boldsymbol G_j^T\\ |
− | E=diag(R_r,R_r,R_r,R_r,R_r,R_r,R_r,R_r) | + | \boldsymbol E=\boldsymbol{diag}(\boldsymbol R_r,\boldsymbol R_r,\boldsymbol R_r,\boldsymbol R_r,\boldsymbol R_r,\boldsymbol R_r,\boldsymbol R_r,\boldsymbol R_r) |
\end{array}</math> | \end{array}</math> | ||
|} | |} | ||
Line 612: | Line 624: | ||
|} | |} | ||
− | where <math>I_{56}</math> is the | + | where <math>\boldsymbol I_{56}</math> is the <math> 5\times 6</math> unit diagonal matrix. <math>\boldsymbol R_r</math> can be found in section 3.1. <math>{\delta }_{ij}</math> , <math>\boldsymbol A_i</math> , <math>\boldsymbol G_j</math> and more detail process of the element kinematics are defined by local and global displacement and rotation vector [9,11]. |
According to the rotate transformation relation obtained by the change of variables from <math>\delta {\phi }_i^g</math> to <math>\delta {\psi }_i^g</math> . | According to the rotate transformation relation obtained by the change of variables from <math>\delta {\phi }_i^g</math> to <math>\delta {\psi }_i^g</math> . | ||
Line 621: | Line 633: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>\boldsymbol B_{m} =\boldsymbol{diag}(\boldsymbol I_{3} ,\boldsymbol T_{m} (\boldsymbol \psi _{1}^{g} ),\boldsymbol I_{3} ,\boldsymbol T_{m} (\boldsymbol \psi _{2}^{g} ),\boldsymbol I_{3} ,\boldsymbol T_{m} (\boldsymbol \psi _{3}^{g} ),\boldsymbol I_{3} ,\boldsymbol T_{m} (\boldsymbol \psi _{4}^{g} ))</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (40) | | style="width: 5px;text-align: right;white-space: nowrap;" | (40) | ||
Line 633: | Line 645: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>T=PE^ | + | | <math>\boldsymbol T=\boldsymbol {PE}^T\boldsymbol B_m</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (41) | | style="width: 5px;text-align: right;white-space: nowrap;" | (41) | ||
|} | |} | ||
− | ==3.3 Element internal force and tangent stiffness matrix== | + | ===3.3 Element internal force and tangent stiffness matrix=== |
− | According to the transformation matrix <math>T</math> in section 3.2, the element global internal force can be defined as | + | According to the transformation matrix <math>\boldsymbol T</math> in section 3.2, the element global internal force can be defined as |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 647: | Line 659: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | [ | + | | <math> \boldsymbol F =\left[\begin{array}{cccccccc} {\boldsymbol F_{1}^{T} } & {\boldsymbol M_{1}^{T} } & {\boldsymbol F_{2}^{T} } & {\boldsymbol M_{2}^{T} } & {\boldsymbol F_{3}^{T} } & {\boldsymbol M_{3}^{T} } & {\boldsymbol F_{4}^{T} } & {\boldsymbol M_{4}^{T} } \end{array}\right]=\boldsymbol T^{T} \boldsymbol f_{l}</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (42) | | style="width: 5px;text-align: right;white-space: nowrap;" | (42) | ||
|} | |} | ||
− | where <math>F_i</math> and <math>M_i</math> denote the internal force vector of the node <math>i</math> for translation and rotation, respectively. The internal force of local deformation <math>f_l</math> is defined and expressed in [20]. | + | where <math>\boldsymbol F_i</math> and <math>\boldsymbol M_i</math> denote the internal force vector of the node <math>i</math> for translation and rotation, respectively. The internal force of local deformation <math>\boldsymbol f_l</math> is defined and expressed in [20]. |
The element local stiffness matrix and established transformation matrix by CR can be found in sections 2.3 and 3.2, Using Equations (21) and (41), the element global tangent stiffness matrix is defined as | The element local stiffness matrix and established transformation matrix by CR can be found in sections 2.3 and 3.2, Using Equations (21) and (41), the element global tangent stiffness matrix is defined as | ||
Line 661: | Line 673: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | <math>K_T=T^ | + | | <math>\boldsymbol K_T=\boldsymbol T^T \boldsymbol K_l^e \boldsymbol T + \boldsymbol K_{\sigma }</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (43) | | style="width: 5px;text-align: right;white-space: nowrap;" | (43) | ||
Line 674: | Line 686: | ||
|- | |- | ||
| <math>\begin{array}{c} | | <math>\begin{array}{c} | ||
− | K_{\sigma }=diag(0_3,K_{\sigma 1,}0_3,K_{\sigma 2,}0_3,K_{\sigma 3,}0_3,K_{\sigma 4}), | + | \boldsymbol K_{\sigma }=\boldsymbol{diag}(0_3,\boldsymbol K_{\sigma 1,}0_3,\boldsymbol K_{\sigma 2,}0_3,\boldsymbol K_{\sigma 3,}0_3,\boldsymbol K_{\sigma 4}),\\ |
− | K_{\sigma i}=\partial T_m^T(M_i)/\partial {\psi }_i(i=1,2,3,4) | + | \boldsymbol K_{\sigma i}=\partial \boldsymbol T_m^T (\boldsymbol M_i)/\partial {\psi }_i(i=1,2,3,4) |
\end{array}</math> | \end{array}</math> | ||
|} | |} | ||
Line 681: | Line 693: | ||
|} | |} | ||
− | In Equation (43), the element tangent stiffness matrix is reduced from a | + | In Equation (43), the element tangent stiffness matrix is reduced from a <math>24\times 24</math> matrix to a <math>20\times 20</math> matrix. In Equation (44), the <math>0_3</math> denotes <math>3\times 3</math> zero matrix, Equation (44) is explained in detail in [11]. |
− | ==4. | + | ==4. Numerical experiments== |
The 4-node thin shell element is used to carry out linear and geometrically nonlinear analysis. This part is mainly to assess the accuracy and efficiency of the element by numerical experiments. All of the numerical results are compared with references and standard 4-node (QUAD 4) thin shell element, which is established and extensively used based on standard updated Lagrangian formulation (UL). The standard 4-node thin shell element tangent stiffness matrix is still a 24×24 block. Generally speaking, the results from the standard finite element with small grid can be used as the standard of structure design. So the numerical results from the present element can be compared with the results of the QUAD 4, for which the model can be divided into grids small enough. At the same time, the convergence feature of the present method was also verified by analytical solution for simply supported orthotropic symmetry plates. | The 4-node thin shell element is used to carry out linear and geometrically nonlinear analysis. This part is mainly to assess the accuracy and efficiency of the element by numerical experiments. All of the numerical results are compared with references and standard 4-node (QUAD 4) thin shell element, which is established and extensively used based on standard updated Lagrangian formulation (UL). The standard 4-node thin shell element tangent stiffness matrix is still a 24×24 block. Generally speaking, the results from the standard finite element with small grid can be used as the standard of structure design. So the numerical results from the present element can be compared with the results of the QUAD 4, for which the model can be divided into grids small enough. At the same time, the convergence feature of the present method was also verified by analytical solution for simply supported orthotropic symmetry plates. | ||
− | ==4.1 Linear analysis== | + | ===4.1 Linear analysis=== |
− | 4.1.1 Simply supported laminated composite plates | + | ====4.1.1 Simply supported laminated composite plates==== |
− | All edges are provided with a simple supported restraint and a uniform load applied to the surface. The plate with 4 | + | All edges are provided with a simple supported restraint and a uniform load applied to the surface. The plate with <math> 4\times 4</math> meshes is used in the analysis. The geometry and material properties of T700S/5405 can be obtained as |
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 698: | Line 710: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>\begin{array}{l} a=b=100mm\quad ,\quad a/h=100\quad ,\quad E_{1} =131100N/mm^{2} \quad ,\quad \\ E_{2} =8580N/mm^{2}\quad ,\quad G_{12} =G_{13} =G_{23} =5000N/mm^{2} \quad ,\quad \quad v_{12} =0.36 \end{array}</math> |
|} | |} | ||
− | | style="width: 5px;text-align: right;white-space: nowrap;" | | + | | style="width: 5px;text-align: right;white-space: nowrap;" | |
|} | |} | ||
+ | |||
+ | The center deflections are presented in the non-dimensional form using the equation [14] | ||
{| class="formulaSCP" style="width: 100%; text-align: center;" | {| class="formulaSCP" style="width: 100%; text-align: center;" | ||
Line 708: | Line 722: | ||
{| style="text-align: center; margin:auto;" | {| style="text-align: center; margin:auto;" | ||
|- | |- | ||
− | | | + | | <math>\bar{w}=w\left(\frac{E_{2} h^{3} }{qa^{4} } \right)\times 10^{3}</math> |
|} | |} | ||
| style="width: 5px;text-align: right;white-space: nowrap;" | (45) | | style="width: 5px;text-align: right;white-space: nowrap;" | (45) | ||
Line 715: | Line 729: | ||
and are compared with standard element and analytical solution in Table 1. The ply stacking sequence and ply angle influence the stiffness and strength of laminated composite plate in aircraft structural design. Hence, the laminated composite plate is more designable than isotropy plate. | and are compared with standard element and analytical solution in Table 1. The ply stacking sequence and ply angle influence the stiffness and strength of laminated composite plate in aircraft structural design. Hence, the laminated composite plate is more designable than isotropy plate. | ||
− | '''Table 1''' | + | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;font-size:75%;">'''Table 1.''' Maximum displacements in Non-dimensional form.</div> |
− | + | {| style="width: 100%;border-collapse: collapse;font-size:85%;" | |
− | + | ||
− | {| style="width: 100%;border-collapse: collapse;" | + | |
|- | |- | ||
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;"| |
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|Analytical solution |
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|standard element (<math>20\times 20</math>) |
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|Present element (<math>4\times 4</math>) |
|- | |- | ||
− | | style="border-top: 1pt solid black;"| | + | | style="border-top: 1pt solid black;"|[45/-45]<sub>s</sub> |
− | | style="border-top: 1pt solid black; | + | | style="border-top: 1pt solid black;text-align: center;"|<math>7.05\times 10^{-4}</math> |
− | | style="border-top: 1pt solid black; | + | | style="border-top: 1pt solid black;text-align: center;"|<math>7.3\times 10^{-4}</math> |
− | | style="border-top: 1pt solid black; | + | | style="border-top: 1pt solid black;text-align: center;"|<math>7.12\times 10^{-4}</math> |
|- | |- | ||
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;text-align: center;"|[90/0]<sub>s</sub> |
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;text-align: center;"|<math>9.54\times 10^{-4}</math> |
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;text-align: center;"|<math>9.98\times 10^{-4}</math> |
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;text-align: center;"|<math>9.67\times 10^{-4}</math> |
|} | |} | ||
− | |||
− | + | ====4.1.2 Simply-supported laminated composite plates==== | |
− | + | The plate with <math>4\times 4</math> meshes is used in the analysis. Aspect ratio (<math>a/h =100</math>) is used to verify present element. The relevant parameters are as follows: | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | {| class="formulaSCP" style="width: 100%; text-align: left;" | |
+ | |- | ||
+ | | | ||
+ | {| style="text-align: center; margin:auto;width: 100%;" | ||
+ | |- | ||
+ | | style="text-align: center;" | <math>\begin{array}{c} | ||
+ | a=b=100mm\quad ,\quad E_1/E_2=25\quad ,\quad E_2=1.0\times {10}^6N/mm^2\\ | ||
+ | G_{12}=G_{13}=0.5\times {10}^6N/mm^2\quad ,\quad G_{23}=0.2\times {10}^6N/mm^2\quad ,\quad v_{12}=0.25 | ||
+ | \end{array}</math> | ||
+ | | style="width: 5px;text-align: right;white-space: nowrap;" | | ||
+ | |} | ||
+ | |} | ||
− | ''' | + | Maximum displacements are obtained using the equation (45), and are compared with the strandard element and the references by Kim ''et al.'' [13-14] and Pucha ''et al.'' [21] in Table 2. The results show that present element can be well applied to the analysis of the actual engineering structure. |
− | Maximum displacements in Non-dimensional form | + | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;font-size:75%;">'''Table 2'''. Maximum displacements in Non-dimensional form.</div> |
− | {| style="width: 100%;border-collapse: collapse;" | + | {| style="width: 100%;border-collapse: collapse;font-size:85%;" |
|- | |- | ||
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align:left;"|Solution method |
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align:center;"|Mesh |
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align:center;"|5/-5 |
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align:center;"|15/-15 |
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align:center;"|30/-30 |
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align:center;"|45/-45 |
|- | |- | ||
− | | style="border-top: 1pt solid black;"| | + | | style="border-top: 1pt solid black;"|Present |
− | | style="border-top: 1pt solid black; | + | | style="border-top: 1pt solid black;text-align:center;"|<math>4\times 4</math> |
− | | style="border-top: 1pt solid black; | + | | style="border-top: 1pt solid black;text-align:center;"|7.083 |
− | | style="border-top: 1pt solid black; | + | | style="border-top: 1pt solid black;text-align:center;"|9.34 |
− | | style="border-top: 1pt solid black; | + | | style="border-top: 1pt solid black;text-align:center;"|10.77 |
− | | style="border-top: 1pt solid black; | + | | style="border-top: 1pt solid black;text-align:center;"|10.25 |
|- | |- | ||
− | | | + | | Strandard element (UL) |
− | | | + | | style="text-align:center;"|<math>40\times 40</math> |
− | | | + | | style="text-align:center;"|7.145 |
− | | | + | | style="text-align:center;"|9.82 |
− | | | + | | style="text-align:center;"|10.36 |
− | | | + | | style="text-align:center;"|10.17 |
|- | |- | ||
− | | | + | | Kim ''et al.'' [13] |
− | | | + | | style="text-align:center;"|<math>4\times 4</math> |
− | | | + | | style="text-align:center;"|7.015 |
− | | | + | | style="text-align:center;"|9.505 |
− | | | + | | style="text-align:center;"|10.52 |
− | | | + | | style="text-align:center;"|10.12 |
|- | |- | ||
− | | | + | | Kim ''et al.'' (EAS-ANS) [14] |
− | | | + | | style="text-align:center;"|<math>4\times 4</math> |
− | | | + | |style="text-align:center;"| 7.025 |
− | | | + | |style="text-align:center;"| 9.468 |
− | | | + | | style="text-align:center;"|10.56 |
− | | | + | | style="text-align:center;"|10.16 |
|- | |- | ||
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;"|Pucha ''et al.'' [21] |
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;"| |
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;text-align:center;"|7.1298 |
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;text-align:center;"|9.1077 |
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;text-align:center;"|9.1718 |
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;text-align:center;"|9.0793 |
|} | |} | ||
− | |||
− | + | ====4.1.3 Simply supported laminated spherical shell==== | |
− | < | + | A laminated spherical shell shown in Figure 5 is analyzed. The surface of the structures is applied with a uniform load. <math>8\times 8</math> meshes are used to analyze the spherical shell structure with lay-ups of (0/90/0/90/0/90/0/90/0). The relevant parameters are used as follows: |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | {| class="formulaSCP" style="width: 100%; text-align: left;" | |
− | + | |- | |
+ | | | ||
+ | {| style="text-align: center; margin:auto;width: 100%;" | ||
+ | |- | ||
+ | | style="text-align: center;" | <math>\begin{array}{c} | ||
+ | h/a=0.001\quad , \quad E_1/E_2=40\quad , \quad E_2=5.1713\times {10}^9N/mm^2\\ | ||
+ | G_{12}/E_2=0.6\quad , \quad G_{23}/E_2=0.6\quad , \quad v_{12}=0.25\quad , \quad a=b | ||
+ | \end{array}</math> | ||
+ | | style="width: 5px;text-align: right;white-space: nowrap;" | | ||
+ | |} | ||
+ | |} | ||
− | <div | + | <div id='img-5'></div> |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | |
+ | |- | ||
+ | | [[Image:draft_Aparicio Nogué_939719944-image137.png|318px]] | ||
+ | |- style="text-align: center; font-size: 75%;" | ||
+ | | colspan="1" | '''Figure 5:''' Geometry proprieties of laminated spherical shell. | ||
+ | |} | ||
− | The center deflection of the laminated spherical shell structure is compared with the non-dimensional form (as shown in Equation (45)) of the standard element and the references by Kim et al. [13,14] and Noor et al. [22]. The results are tabulated in Table 3. | + | The center deflection of the laminated spherical shell structure is compared with the non-dimensional form (as shown in Equation (45)) of the standard element and the references by Kim ''et al.'' [13,14] and Noor ''et al.'' [22]. The results are tabulated in Table 3. |
− | |||
− | Non-dimensional center deflection of laminated spherical shell | + | <div class="center" style="width: auto; margin-left: auto; margin-right: auto;font-size:75%;">'''Table 3'''. Non-dimensional center deflection of laminated spherical shell.</div> |
− | {| style="width: 100%;border-collapse: collapse;" | + | {| style="width: 100%;border-collapse: collapse;font-size:85%" |
|- | |- | ||
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;"|Solution method |
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align:center;"|Mesh |
− | | style="border-top: | + | | style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align:center;"|Non-dimensional center deflection |
|- | |- | ||
− | | style="border-top: 1pt solid black;"| | + | | style="border-top: 1pt solid black;"|Present |
− | | style="border-top: 1pt solid black; | + | | style="border-top: 1pt solid black;text-align:center;"|<math>8\times 8</math> |
− | | style="border-top: 1pt solid black; | + | | style="border-top: 1pt solid black;text-align:center;"|5.874E-05 |
|- | |- | ||
− | | | + | | Standard element (UL) |
− | | | + | | style="text-align:center;"|<math>40\times 40</math> |
− | | | + | | style="text-align:center;"|5.92E-05 |
|- | |- | ||
− | | | + | | Kim ''et al.'' [13] |
− | | | + | |style="text-align:center;"| <math>8\times 8</math> |
− | | | + | | style="text-align:center;"|5.844E-05 |
|- | |- | ||
− | | | + | | Kim ''et al.'' (EAS-ANS) [14] |
− | | | + | | style="text-align:center;"|<math>8\times 8</math> |
− | | | + | | style="text-align:center;"|5.890E-05 |
|- | |- | ||
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;"|Noor ''et al.'' [22] |
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;text-align:center;"|<math>8\times 8</math> |
− | | style="border-bottom: | + | | style="border-bottom: 1pt solid black;text-align:center;"|5.916E-05 |
|} | |} | ||
− | ==4.2 Nonlinear analysis== | + | ===4.2 Nonlinear analysis=== |
− | 4.2.1 Laminated composite plate | + | ====4.2.1 Laminated composite plate==== |
− | A multi-layer laminated composite plate is used to nonlinear analyze. 8 | + | A multi-layer laminated composite plate is used to nonlinear analyze. <math>8\times8</math> meshes and lay-ups of [45/-45/0/0/45/-45/90/90] are used in this analysis. The length of the plate is <math>a=b=254mm</math>. The total thickness is <math>h=2.114mm</math>. All edges are clamped. The surface of the structures is applied with a uniform load. The material properties are used as follows: |
− | + | {| class="formulaSCP" style="width: 100%; text-align: left;" | |
− | + | |- | |
− | E_1=13.1\times {10}^4N/mm^2,\ | + | | |
− | G_{12}=G_{13}=0.641\times {10}^4N/mm^2,\ | + | {| style="text-align: center; margin:auto;width: 100%;" |
− | \end{array}</math> | + | |- |
+ | | style="text-align: center;" | <math>\begin{array}{c} | ||
+ | E_1=13.1\times {10}^4N/mm^2\quad ,\quad E_2=E_3=1.303\times {10}^4N/mm^2\quad ,\quad v_{12}=v_{23}=v_{13}=0.38\\ | ||
+ | G_{12}=G_{13}=0.641\times {10}^4N/mm^2\quad ,\quad G_{23}=0.4721\times {10}^4N/mm^2 | ||
+ | \end{array}</math> | ||
+ | | style="width: 5px;text-align: right;white-space: nowrap;" | | ||
+ | |} | ||
+ | |} | ||
− | Maximum displacements of the plate is compared with that of the standard element with 40 | + | Maximum displacements of the plate is compared with that of the standard element with <math>40\times 40</math> meshes and the references by Kim ''et al.'' [14] and Lee ''et al.'' [23]. The results are shown in [[#img-6|Figure 6]], normalized deflection: <math>\overline{w}=w/h</math>. |
− | <div | + | <div id='img-6'></div> |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: 100%;" | |
+ | |- | ||
+ | | style="padding:10px;"|[[Image:draft_Aparicio Nogué_939719944-image143.png|416px]] | ||
+ | |- style="text-align: center; font-size: 75%;" | ||
+ | | colspan="1" | '''Figure 6:''' Load-deflection curve of laminated composite plate. | ||
+ | |} | ||
− | + | ====4.2.2 Pinched cantilever cylinder==== | |
− | + | ||
− | + | In this example, nonlinear analysis of pinched cantilever cylinder is carried out. The cylinder is depicted in Figure [[#img-7|7]], which clamped at B face and applied two opposite loads at A face. <math> 16\times 16</math> regular meshes are used to analyze the structure. | |
− | + | <div id='img-7'></div> | |
+ | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: 100%;" | ||
+ | |- | ||
+ | | style="padding:10px;|[[Image:draft_Aparicio Nogué_939719944-image144.png|400px]] | ||
+ | |- style="text-align: center; font-size: 75%;" | ||
+ | | colspan="1" | '''Figure 7:''' Pinched cantilever cylinder. | ||
+ | |} | ||
− | < | + | It is clear that this example reflects the nonextensional deformations for the element, this phenomenon will not occur in actual engineering. By the large deformation analysis, the vertical displacement at A of the cantilever cylinder is compared with that of the standard element (UL) with <math> 50\times 50</math> and the references by Battini [11] and Norachan ''et al.'' [24] as shown in Figure [[#img-8|8]]. The results analyzed by present element showed a satisfied precision and efficiency. |
− | + | ||
− | <div | + | <div id='img-8'></div> |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: 100%;" | |
+ | |- | ||
+ | | style="padding:10px;|[[File:Draft_Samper_831316476_5609_Fig8.png|350px]] | ||
+ | |- style="text-align: center; font-size: 75%;" | ||
+ | | colspan="1" | '''Figura 8:''' Load-displacement curve of cantilever cylinder. | ||
+ | |} | ||
− | + | ====4.2.3 Double-curved thin shell structure==== | |
− | + | Double-curved thin shell structure can be seen everywhere in the actual project. The numerical example is to use present element to carry out the nonlinear analysis of the double-curved structure, which is always used in many engineering structures and shows significant different nonlinear. <math>16\times 10</math> meshes and lay-ups of <math>[45/90/0/-45]_{S}</math> are used in this analysis. The length of <math>b=2000 mm</math> and width of <math>a=1250 mm</math> are shown in Figure 9, and the radius of the two sides in the direction of width are <math>r_1 = 318mm</math>, <math>r_2 = 504mm</math> and <math>r'_1 = 1015mm</math>, <math>r'_2 = 772mm</math>, respectively, the length of the section with radius as <math>r_1</math> in the direction of width is <math>c=645mm</math>, the length of the section with radius as <math>r'_1</math> in the direction of width is <math>d=860mm</math>, the total thickness of the laminates is <math>h=3.2mm</math>. The material properties of T700/QY9511 are used as follows: | |
− | [ | + | |
− | + | {| class="formulaSCP" style="width: 100%; text-align: left;" | |
− | + | |- | |
+ | | | ||
+ | {| style="text-align: center; margin:auto;width: 100%;" | ||
+ | |- | ||
+ | | <math>\begin{array}{l} {E_{1} =1.334\times 10^{5} N/mm,\quad \quad E_{2} =1.008\times 10^{4} N/mm} \\ {G_{12} =G_{13} =G_{23} =4.9\times 10^{3} N/mm,\quad \quad v_{12} =0.31} \end{array}</math> | ||
+ | | style="width: 5px;text-align: right;white-space: nowrap;" | | ||
+ | |} | ||
+ | |} | ||
− | + | <div id='img-9'></div> | |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | |
− | + | |- | |
− | + | | [[Image:draft_Aparicio Nogué_939719944-image157.png|400px]] | |
− | <div | + | |- style="text-align: center; font-size: 75%;" |
− | + | | colspan="1" | '''Figure 9:''' Geometry proprieties of double-curved thin shell structure. | |
− | + | |} | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
− | + | ||
The normal displacement at A of the double-curved thin shell structure is compared with that of the standard 4-node thin shell element with 80×40 meshes in Figure 10, displacements of double-curved thin shell structure are obtained by load of 1.2 KN using standard element and present element, respectively, as shown in Figure 11. The relative error between present element with big grid and standard element with small grid is 4%, and the relative error between big grid and small grid models using standard element is 56.7%. The results showed that the present element is acceptable for the double-curved structure. | The normal displacement at A of the double-curved thin shell structure is compared with that of the standard 4-node thin shell element with 80×40 meshes in Figure 10, displacements of double-curved thin shell structure are obtained by load of 1.2 KN using standard element and present element, respectively, as shown in Figure 11. The relative error between present element with big grid and standard element with small grid is 4%, and the relative error between big grid and small grid models using standard element is 56.7%. The results showed that the present element is acceptable for the double-curved structure. | ||
− | <div | + | <div id='img-10'></div> |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;" | |
+ | |- | ||
+ | | [[Image:draft_Aparicio Nogué_939719944-image158.png|404px]] | ||
+ | |- style="text-align: center; font-size: 75%;" | ||
+ | | colspan="1" | '''Figure 10:''' Load-displacement curve of double-curved thin shell structure. | ||
+ | |} | ||
− | <div | + | <div id='img-11'></div> |
− | + | {| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: 100%;" | |
− | + | ||
− | + | ||
− | + | ||
− | + | ||
|- | |- | ||
− | | | + | | style="padding:10x[[Image:draft_Aparicio Nogué_939719944-image159.png|402px]] |
|- | |- | ||
− | | | + | | (a) using standard element (80×40) |
|- | |- | ||
− | | | + | | [[Image:draft_Aparicio Nogué_939719944-image160.png|396px]] |
|- | |- | ||
− | | | + | | (b) using standard element (16×10) |
|- | |- | ||
− | | | + | | [[Image:draft_Aparicio Nogué_939719944-image161.png|396px]] |
|- | |- | ||
− | | | + | | (c) using present element (16×10) |
+ | |- style="text-align: center; font-size: 75%;" | ||
+ | | colspan="1" | '''Figure 11:''' Displacements of laminated thin cylinder by 1.2 KN. | ||
|} | |} | ||
− | ==5. | + | ==5. Conclusions== |
This paper mainly presented a highly efficient nonlinear computer modelling approach for laminated composite thin shell structures. The analysis method mainly includes 3 points: | This paper mainly presented a highly efficient nonlinear computer modelling approach for laminated composite thin shell structures. The analysis method mainly includes 3 points: | ||
Line 940: | Line 992: | ||
==References== | ==References== | ||
+ | <div class="auto" style="width: auto; margin-left: auto; margin-right: auto;font-size: 85%;"> | ||
[1] G.A. Wempner, Finite elements, finite rotation and small strains of flexible shells, International Journal of Solids and Structures, 5 (1969) 117-153. | [1] G.A. Wempner, Finite elements, finite rotation and small strains of flexible shells, International Journal of Solids and Structures, 5 (1969) 117-153. | ||
Line 988: | Line 1,041: | ||
[24] P. Norachan, S. Suthasupradit, K.D. Kim, A co-rotational 8-node degenerated thin-walled element with assumed natural strain and enhanced assumed strain, Finite Elements in Analysis and Design, 50 (2012) 70-85. | [24] P. Norachan, S. Suthasupradit, K.D. Kim, A co-rotational 8-node degenerated thin-walled element with assumed natural strain and enhanced assumed strain, Finite Elements in Analysis and Design, 50 (2012) 70-85. | ||
+ | </div> |
The finite formulation of the 4-node thin shell element based on Co-rotational (CR) and Timoshenko's theory was presented for the analysis of laminated composite structures. Based on the Timoshenko's theory the present thin shell element avoids shear locking behavior, and a bilinear in-plane displacement field is introduced for the coupling of in-plane and bending actions, and such element performs simple formulation and highly efficienct. A highly efficient CR formulation was also established to define the motion of the element with large displacement. The characteristics of the element are more pronounced when shells get thin, each node of the 4-node element has five degrees of freedom. An incremental iterative method based on the Newton Raphson method was used to solve the nonlinear equilibrium equations. A number of numerical examples were given to verify that the formula is computationally efficient and the results showed good agreement.
Many nonlinear shell element formulations are based on the total Lagrangian (TL) and update Lagrangian (UL), including complex formulations and increasing calculation interval. However, in order to maximum the computational efficiency, a nonlinear finite element method (FEM) is required. A named “CR” approach was introduced by Wempner [1] and Belytschko [2] firstly.
As a new method for improving the computational efficiency of nonlinear element with large displacement problems, the CR has attracted more and more researchers in the last decade. The core of the method is to decompose the analytic process of the element into large displacement and large rotation of the rigid body and local small deformation. When the CR method is used in the process of structural geometric nonlinear analysis, the deformation in the local coordinate system reflects the small strain of the element. However, the translation and rotation of the rigid body reflects the large displacement of the element.
Most of the studies mainly apply the CR method to beam and shell structures. Crisfield et al. [3] based on the CR method to establish a universal model for the analysis of beam, shell and solid structure, the formulas used to build the model are simpler than many of the earlier procedures. Afterwards Hsial et al. [4] and Wen et al. [5] established an axisymmetric thin-walled beam element with second-order accuracy using CR and TL formulation. Tham et al. [6] analyses multi-layered composite plate and shell structures using CR kinematic framework. Cortivo et al. [7] carried out plastic buckling and failure analysis for thin shell structures using CR and ANDES finite element formulations. A triangular shell element formula based on CR method has been obtained by Nour-Omid et al. [8], and then Pacoste [9] and Eriksson et al. [10] developed the instability analysis of the shell structure using the established triangular plate element CR formulations, the difference between them is the parameterization of the rotations.
In order to improve the computational efficiency, Battini [11] introduced three modifications in the CR framework, including local rotations, the number of element local degrees of freedom from 18 to 15 and the parameterization of the global finite rotations. Kim et al. [12] believes that the actual superimposition and the correction factor are not necessary and improve the the transverse shear stiffness, according to it Kim et al. [13] introduced the ANS method and CR formulation with second order accuracy and defines the transverse shear stiffness. To ensure computational efficiency, Kim et al. [14] also improved the 4-node shell element by combining an EAS and ANS. Nowadays, this approach based on beam [15,17,19], rod [18] and shell [16] element is widely used in many fields.
The purpose of the paper is to improve the computational efficiency of the thin shell structures modelling by a high efficiency nonlinear 4-node thin shell element. In the local element formulation, a simple and high efficiency laminated composite element was presented et al. [20], the tangential rotation and the shear strain were obtained by Timoshenko laminated composite beam theory, and a bilinear in-plane displacement field was introduced for the coupling of in-plane and bending actions. In order to consider the movement of the element with large translation and rotation, a highly efficienct CR method was obtained. In the current work, the main idea introduced by Pacoste [9] and Battini [11] was employed. At last, a high efficiency tangent stiffness matrix of 4-node thin shell element was obtained in this paper.
In local coordinate, Cen [20] show that tangential rotation and shear strain were obtained by Timoshenko beam theory. Shear strain and rotation fields within the element are then determined by improved rational interpolation, and a bilinear displacement field is introduced for the coupling of in-plane and bending actions. Finally, in order to analyze random laminated composite structures, the 20 degrees of freedom quadrilateral bending element can be used.
As shown in Figure1, the local deformational displacements and rotations for random 4-node elements can be written as
|
(1) |
Figure 1: Local displacement and rotation at the mid-plane of a laminated composite element. |
based on the local deformation of the element mid-plane, displacement and rotation of random point for the elmenet can be written as
|
(2) |
The element is shown in Figure 2. Based on Timoshenko beam theory, the tangential rotation and the shear strain along each element side can be written as
|
(3) |
|
(4) |
with
|
(5) |
Figure 2: Laminated composite beam element. |
Based on the Equation (4), element shear strain of every edge can be written as
|
(6) |
with
|
(7) |
where is the length of element edge, is the y orientation distance between two adjacent nodes, is the orientation distance between two adjacent nodes. , are the shear and bending stiffness.
According to the element shear strain of every edge , can be written as
|
(8) |
Based on the geometrical relationship of two adjacent edges, shear strain on the crossing point 1 can be written as
|
(9) |
Shear strain of the other nodes can be written similarly. So shear strain field can be written as
|
(10) |
with
|
(11) |
The element normal and tangential of the first edge is shown in Figure 1, each edge of the element is assumed to be linearly distributed.
Based on the Equation (3), the normal and tangential rotations of the first edge can be written as
|
(12) |
Normal and tangential rotations of the other edges can be obtained similarly.
In global coordinate, Normal and tangential rotations of the midpoint on the first edge are represented as
|
(13) |
Element rotation field can be written as
|
(14) |
with
|
(15) |
The bilinear in-plane displacement and can be written as
|
(16) |
where are the quadrilateral element shape functions.
Element in-plane strain can be written as
|
(17) |
with
|
(18) |
According to section 2.2, element curvature field can be given
|
(19) |
where is bendding strain matrix.
Based on the constitutive equation of laminated composite plate, the relationship between strain and displacement in mid-plane can be written as
|
(20) |
The stiffness matrix of the 4-node element is derived from the principle of minimum potential energy, element shear strain and in-plane strain field:
|
(21) |
where shear strain matrix and in-plane strain matrix of the mid-plane can be found in Equations (10) and (20), respectively, is stress-strain matrix of laminated composite plate, is shear stiffness, and is the Jacobian determinant. The element stiffness matrix was reduced to a matrix, which was calculated by Gauss integral.
CR approach is more efficient than TL and UL to analyze thin wall structure models. Based on CR description, a lot of plate and shell elements including 3 node and 4 node were established and confirmed by many numerical illustrations and tests. According to CR description, element kinematics from initial undeformed to final deformed attitude can be divided into two sections [7]. The first section is rigid rotation removing deformational part.and translation. The second step is a local deformation and rotation under local coordinate system. Element deformation displacement and rotation were established conveniently in section 2.
In order to represent large three-dimensional rotations, a parameterization approach of the rigid rotations is introduced [9]. Orthogonal matrix can be obtained by three separate parameter variables. Firstly, the rotational vector can be defined by
|
According to the definition above, as shown in Figure 3, rotation can be obtained by an angle defined by the unit vector . So the value of is obtained by
|
(23) |
where () are the components of .
Figure 3: Rotational vector. |
In terms of the definition above, the orthogonal matrix can be written as
|
(24) |
Based on the conception of the rigid rotation, the global rotation matrix at point can be defined by .
|
(25) |
Figure 4: Element kinematics and coordinate systems. |
As shown in Figure 4, we define as the vectors of the local frame in the current deformed configuration. Then, the matrix [9] defining rotations can be obtained as
|
(26) |
with
|
(27) |
where and denote displacement from the initial to current coordinate, respectively. and denote displacement vector of initial nodes and in the global system, respectively. and denote displacement vector along the orientations and of the line connecting two nodes, respectively.
The local deformation of point can be expressed by , as shown in section 2. According to the define of rigid frame in the local coordinate, can be redefined as
|
(28) |
As shown in Figure 4, noting that , the local rotations are defined by the matrices
as
|
(29) |
The differentiation can be given by
|
(30) |
where , , and the corresponding associated vectors , , denote spatial angular variations, which superimposed onto orthogonal matrices , and .
Using Equations (28), (29) and (30), the differentiation of Equations (28) can be given as
|
(31) |
with
|
(32) |
The differentiation of Equation (30) can be given as
|
(33) |
According to Equations (30) and (33), the following expression is obtained as
|
(34) |
The transformation matrix is needed to convert the internal force and tangent stiffness matrix. In the global coordinate, the 4-node element displacement is defined as
|
(35) |
where are the global displacement vectors of the nodes, denotes the spatial angular variations, defined in Equation (30).
Referring to Equations (23), (25) and (30), the transformation matrix was defined by
|
(36) |
with
|
(37) |
where are the components of , is the unit diagonal matrix.
Using Equations (31), (32) and (34), the matrix can be written as
|
(38) |
The blocks is the matrix removing deformational part:
|
(39) |
where is the unit diagonal matrix. can be found in section 3.1. , , and more detail process of the element kinematics are defined by local and global displacement and rotation vector [9,11].
According to the rotate transformation relation obtained by the change of variables from to .
|
(40) |
Using Equations (38) and (40), the final transformation matrix is obtained as
|
(41) |
According to the transformation matrix in section 3.2, the element global internal force can be defined as
|
(42) |
where and denote the internal force vector of the node for translation and rotation, respectively. The internal force of local deformation is defined and expressed in [20].
The element local stiffness matrix and established transformation matrix by CR can be found in sections 2.3 and 3.2, Using Equations (21) and (41), the element global tangent stiffness matrix is defined as
|
(43) |
with
|
(44) |
In Equation (43), the element tangent stiffness matrix is reduced from a matrix to a matrix. In Equation (44), the denotes zero matrix, Equation (44) is explained in detail in [11].
The 4-node thin shell element is used to carry out linear and geometrically nonlinear analysis. This part is mainly to assess the accuracy and efficiency of the element by numerical experiments. All of the numerical results are compared with references and standard 4-node (QUAD 4) thin shell element, which is established and extensively used based on standard updated Lagrangian formulation (UL). The standard 4-node thin shell element tangent stiffness matrix is still a 24×24 block. Generally speaking, the results from the standard finite element with small grid can be used as the standard of structure design. So the numerical results from the present element can be compared with the results of the QUAD 4, for which the model can be divided into grids small enough. At the same time, the convergence feature of the present method was also verified by analytical solution for simply supported orthotropic symmetry plates.
All edges are provided with a simple supported restraint and a uniform load applied to the surface. The plate with meshes is used in the analysis. The geometry and material properties of T700S/5405 can be obtained as
|
The center deflections are presented in the non-dimensional form using the equation [14]
|
(45) |
and are compared with standard element and analytical solution in Table 1. The ply stacking sequence and ply angle influence the stiffness and strength of laminated composite plate in aircraft structural design. Hence, the laminated composite plate is more designable than isotropy plate.
Analytical solution | standard element () | Present element () | |
[45/-45]s | |||
[90/0]s |
The plate with meshes is used in the analysis. Aspect ratio () is used to verify present element. The relevant parameters are as follows:
|
Maximum displacements are obtained using the equation (45), and are compared with the strandard element and the references by Kim et al. [13-14] and Pucha et al. [21] in Table 2. The results show that present element can be well applied to the analysis of the actual engineering structure.
Solution method | Mesh | 5/-5 | 15/-15 | 30/-30 | 45/-45 |
Present | 7.083 | 9.34 | 10.77 | 10.25 | |
Strandard element (UL) | 7.145 | 9.82 | 10.36 | 10.17 | |
Kim et al. [13] | 7.015 | 9.505 | 10.52 | 10.12 | |
Kim et al. (EAS-ANS) [14] | 7.025 | 9.468 | 10.56 | 10.16 | |
Pucha et al. [21] | 7.1298 | 9.1077 | 9.1718 | 9.0793 |
A laminated spherical shell shown in Figure 5 is analyzed. The surface of the structures is applied with a uniform load. meshes are used to analyze the spherical shell structure with lay-ups of (0/90/0/90/0/90/0/90/0). The relevant parameters are used as follows:
|
Figure 5: Geometry proprieties of laminated spherical shell. |
The center deflection of the laminated spherical shell structure is compared with the non-dimensional form (as shown in Equation (45)) of the standard element and the references by Kim et al. [13,14] and Noor et al. [22]. The results are tabulated in Table 3.
Solution method | Mesh | Non-dimensional center deflection |
Present | 5.874E-05 | |
Standard element (UL) | 5.92E-05 | |
Kim et al. [13] | 5.844E-05 | |
Kim et al. (EAS-ANS) [14] | 5.890E-05 | |
Noor et al. [22] | 5.916E-05 |
A multi-layer laminated composite plate is used to nonlinear analyze. meshes and lay-ups of [45/-45/0/0/45/-45/90/90] are used in this analysis. The length of the plate is . The total thickness is . All edges are clamped. The surface of the structures is applied with a uniform load. The material properties are used as follows:
|
Maximum displacements of the plate is compared with that of the standard element with meshes and the references by Kim et al. [14] and Lee et al. [23]. The results are shown in Figure 6, normalized deflection: .
Figure 6: Load-deflection curve of laminated composite plate. |
In this example, nonlinear analysis of pinched cantilever cylinder is carried out. The cylinder is depicted in Figure 7, which clamped at B face and applied two opposite loads at A face. regular meshes are used to analyze the structure.
Figure 7: Pinched cantilever cylinder. |
It is clear that this example reflects the nonextensional deformations for the element, this phenomenon will not occur in actual engineering. By the large deformation analysis, the vertical displacement at A of the cantilever cylinder is compared with that of the standard element (UL) with and the references by Battini [11] and Norachan et al. [24] as shown in Figure 8. The results analyzed by present element showed a satisfied precision and efficiency.
Figura 8: Load-displacement curve of cantilever cylinder. |
Double-curved thin shell structure can be seen everywhere in the actual project. The numerical example is to use present element to carry out the nonlinear analysis of the double-curved structure, which is always used in many engineering structures and shows significant different nonlinear. meshes and lay-ups of are used in this analysis. The length of and width of are shown in Figure 9, and the radius of the two sides in the direction of width are , and , , respectively, the length of the section with radius as in the direction of width is , the length of the section with radius as in the direction of width is , the total thickness of the laminates is . The material properties of T700/QY9511 are used as follows:
|
Figure 9: Geometry proprieties of double-curved thin shell structure. |
The normal displacement at A of the double-curved thin shell structure is compared with that of the standard 4-node thin shell element with 80×40 meshes in Figure 10, displacements of double-curved thin shell structure are obtained by load of 1.2 KN using standard element and present element, respectively, as shown in Figure 11. The relative error between present element with big grid and standard element with small grid is 4%, and the relative error between big grid and small grid models using standard element is 56.7%. The results showed that the present element is acceptable for the double-curved structure.
Figure 10: Load-displacement curve of double-curved thin shell structure. |
style="padding:10x |
(a) using standard element (80×40) |
(b) using standard element (16×10) |
(c) using present element (16×10) |
Figure 11: Displacements of laminated thin cylinder by 1.2 KN. |
This paper mainly presented a highly efficient nonlinear computer modelling approach for laminated composite thin shell structures. The analysis method mainly includes 3 points:
The test results presented in Section 4 show that all of the numerical results are the same, but the element presented in this paper is more efficient than the standard element. Especially, large deformation analysis of pinched cantilever cylinder (Example 4.2.2) and double-curved thin shell (Example 4.2.3) can be used to detect the results of geometric nonlinear analysis for other engineering structures. It is indicated that the present element is more effective since it can save much computation time, which is more applicable for engineering computation. The present work can be applied in designing and analyzing composite thin shell structure.
The work presented in this article was supported by the National Natural Science Foundation of China (Grant Nos. 51175424 & 51475369), the Basic Research Foundation of Northwestern Polytechnical University (Grant No. JC201238), and the Aeronautical Science Foundation of China (Grant No. 2016ZD53036).
[1] G.A. Wempner, Finite elements, finite rotation and small strains of flexible shells, International Journal of Solids and Structures, 5 (1969) 117-153.
[2] T. Belytschko, L. Schwer, Non-linear transient finite element analysis with convected co-ordinates, International Journal of Numerical Methods and Engineering, 7(9) (1973) 255-271.
[3] M.A. Crisfield, G.F. Moita, A unified co-rotational framework for solids, shells and beams, Solids Structures, 33 (1996) 2969-2992.
[4] K.M. Hsial, Y.L. Wen, A co-rotational formulation for thin-walled beams with monosymmetric open section, Computer Methods in Applied Mechanics and Engineering, 190 (2000) 1163-1185.
[5] Y.L. Wen, K.M. Hsial, Co-rotational formulation for geometric nonlinear analysis of doubly symmetric thin-walled beams, Computer Methods in Applied Mechanics and Engineering, 190 (2001) 6023-6052.
[6] C.L. Tham, Z. Zhang, A. Masud, An elasto-plastic damage model cast in a co-rotational kinematic framework for large deformation analysis of laminated composite shells, Computer Methods in Applied Mechanics and Engineering, 194 (2005) 2641-2660.
[7] N.D. Cortivo, C.A. Felippa, H. Bavestrello, W.T.M. Silva, Plastic buckling and collapse of thin shell structures, using layered plastic modeling and co-rotational ANDES finite elements, Computer Methods in Applied Mechanics and Engineering, 198 (2009) 785-798.
[8] B. Nour-Omid, C.C. Rankin, Finite rotation analysis and consistent linearization using projectors, Computer Methods in Applied Mechanics and Engineering, 93 (1991) 353-384.
[9] C. Pacoste, Co-rotational flat facet triangular elements for shell instability analysis, Computer Methods in Applied Mechanics and Engineering, 156 (1998) 75-110.
[10] A. Eriksson, C. Pacoste, Element formulation and numerical techniques for stability problems in shells, Computer Methods in Applied Mechanics and Engineering, 191 (2002) 3775-3810.
[11] J.M. Battini, A modified corotational framework for triangular shell elements, Computer Methods in Applied Mechanics and Engineering, 196 (2007) 1905-1914.
[12] K,D. Kim, G.R. Lomboy, S.C. Han, A co-rotational 8-node assumed strain shell element for postbuckling analysis of laminated composite plates and shells, Comput. Mech. 30(4) (2003) 330-42.
[13] K.D. Kim, C.S. Lee, S.C. Han, A 4-node co-rotational ANS shell element for laminated composite structures, Composite Structures, 80 (2007) 234-252.
[14] K.D. Kim, S.C. Han, S. Suthasupradit, Geometrically non-linear analysis of laminated composite structures using a 4-node co-rotational shell element with enhanced strains, International Journal of Non-Linear Mechanics, 42 (2007) 864-881.
[15] R. Alsafadie, M. Hjiaj, J.M. Battini, Corotational mixed finite element formulation for thin-walled beams with generic cross-section, Comput. Methods Appl. Mech. Engrg. 199 (2010) 3197-3212.
[16] F.S. Almeida, A.M. Awruch, Corotational nonlinear dynamic analysis of laminated composite shells, Finite Elements in Analysis and Design, 47 (2011) 1131-1145.
[17] T.N. Le, J.M. Battini, M.A. Hjiaj, Consistent 3D corotational beam element for nonlinear dynamic analysis of flexible structures, Comput. Methods Appl. Mech. Engrg. 269 (2014) 538-565.
[18] S. Faroughi, H.H. Khodaparast, M.I. Friswell, Non-linear dynamic analysis of tensegrity structures using a co-rotational method, International Journal of Non-Linear Mechanics, 69 (2015) 55-65.
[19] W. Wang, X.P. Zhu, Z. Zhou, J.B. Duan, A method for nonlinear aeroelasticity trim and stability analysis of very flexible aircraft based on co-rotational theory, Journal of Fluids and Structures, 62 (2016) 209-229.
[20] S. Cen, Y.Q. Long, Z.H. Yao, A new element based on the first-order shear deformation theory for the analysis of laminated composite plates, Engineering Mechanics, 1(19) (2002) 1-8.
[21] N.S. Pucha, J.N. Reddy, A mixed shear flexible finite element for the analysis of laminated plates, Comput. Meth. Appl. Mech. Eng. 44(2) (1984) 213-27.
[22] N A.K. oor, M.D. Mathers, Anisotropy and shear deformation in laminated composite plates, AIAA J, 14 (1975) 282-5.
[23] K.D. Lee, K W. anok-Nukulchai, A nine-node assumed strain finite element for large deformation analysis of laminated shells, Int. J. Num. Meth. Eng. 42 (1998) 777-798.
[24] P. Norachan, S. Suthasupradit, K.D. Kim, A co-rotational 8-node degenerated thin-walled element with assumed natural strain and enhanced assumed strain, Finite Elements in Analysis and Design, 50 (2012) 70-85.
Published on 03/01/18
Accepted on 05/06/17
Submitted on 02/03/17
Volume 34, Issue 1, 2018
DOI: 10.23967/j.rimni.2017.8.001
Licence: CC BY-NC-SA license
Are you one of the authors of this document?