Geometrically nonlinear static and dynamic analysis of composite laminates shells with a triangular finite element

TECHNICAL PAPERS

Geometrically nonlinear static and dynamic analysis of composite laminates shells with a triangular finite element

Liércio André Isoldi^I; Armando Miguel Awruch^II; Paulo Roberto de F. Teixeira^III; Inácio B. Morsch^IV

^[e-mail]. Federal Univ. of Rio Grande do Sul - UFRGS. Graduate Program in Mechanical Engineering. 90035-190 Porto Alegre, RS, Brazil

^IISenior Member, ABCM, [e-mail]. Federal Univ. of Rio Grande do Sul - UFRGS. Graduate Program in Mechanical Engineering/90035-190 Porto Alegre, RS, Brazil

^III [e-mail]

^IV [e-mail]

ABSTRACT

Geometrically nonlinear static and dynamic behaviour of laminate composite shells are analyzed in this work using the Finite Element Method (FEM). Triangular elements with three nodes and six degrees of freedom per node (three displacement and three rotation components) are used. For static analysis the nonlinear equilibrium equations are solved using the Generalized Displacement Control Method (GDCM) while the dynamic solution is performed using the classical Newmark Method with an Updated Lagrangean Formulation (ULF). The system of equations is solved using the Gradient Cojugate Method (GCM) and in nonlinear cases with finite rotations and displacements an iterative-incremental scheme is employed. Numerical examples are presented and compared with results obtained by other authors with different kind of elements and different schemes.

Keywords: geometrically nonlinear analysis, laminate composite, static and dynamic analysis of shells

Introduction

¹It is well known that laminate composite materials are nowadays commonly used in the aeronautical, aerospace, naval and other industries mainly because their attractive properties as compared to isotropic materials, such as higher stiffness/weight, higher strength, higher damping and good properties related to thermal or acoustic isolation, among others.

A triangular finite element (GPL-T9) presented previously by Zhang, Lu and Kuang (1998) and by Teixeira (2001) for isotropic materials was extended to geometrically nonlinear shell analyses.

It was considered the Classical Lamination Theory (CLT) given by Jones (1999), where the complete laminate, having several layers, is analyzed as an equivalent material with only one layer.

To analyze geometrically nonlinear static problems the Generalized Displacement Control Method (Yang and Shieh, 1990) was used because it is suitable for problems with multiple critical points.

For geometrically nonlinear dynamic problems the Newmark’s Method and an Updated Lagrangean Formulation (Bathe, 1996) were employed.

The system of equations was solved using a Pre-conditioned Gradient Conjugate Method and an incremental/iterative scheme (Teixeira, 2001).

Several examples are analyzed and compared with results obtained by other authors, showing that this element, where its mass and stiffness matrices can be implemented analytically, is able to solve structures involving thin plates and shells of composite materials very efficiently.

Nomenclature

A =element area, m² (in²)

a =coefficients of Newmark’s method

B =strain-displacement matrix

C =damping matrix

D =material properties matrix

E =elastic moduli N/m² (lb/in²)

e =linear strain m/m (in/in)

F =nodal vector of internal forces, N (lb)

f =forces, N (lb)

G =shear modulus, N/m² (lb/in²)

h =thickness, m (in)

H =interpolation functions

K =stiffness matrix

L =area coordinates

M =nodal vector of bending internal forces

Mm =mass matrix

MN =nodal vector of bending-membrane coupling internal force

N =nodal vector of membrane internal forces

NM =nodal vector of membrane-bending coupling internal force

Q =elastic constants

 =external forces, N (lb)

R =nodal vector of external loads

r =transformation matrix

S =second Piola-Kirchhoff stress tensor, N/m² (lb/in²)

T =membrane internal forces, N (lb)

t =time, s (s)

u =displacement, m (in)

V =domain (volume), m³ (in³)

w =vertical displacement, m (in)

X =deformation gradient

x =cartesian coordinate

y =cartesian coordinate

z = cartesian coordinate

Greek Symbols

a =coefficient of Newmark’s method

g = angle between x_g and x_l, rad (rad)

D =increment

d =coefficient of Newmark’s method

e =Green-Lagrange strain tensor, m/m (in/in)

h =nonlinear strain m/m (in/in)

Q =rotations

q =angle between x_g and 1, rad (rad)

L =matrix of cosines

l =cosine of angles between local and global coordinate systems

n =Poisson coefficient

r =specific mass, kg/m³ (lb.s²/in⁴)

t =Cauchy stress tensor, N/m² (lb/in²)

Subscripts

b relative to the bending

bm relative to the bending-membrane coupling

g relative to the global

i,j,m index

k relative to the layer

L linear

l relative to the local

m relative to the membrane

mb relative to the membrane-bending coupling

NL nonlinear

r,s index

x relative to x coordinate

y relative to y coordinate

z relative to z coordinate

g relative to the g angle

Superscripts

B relative to the body

e relative to the element

i iteration

k relative to the layer

n number of layers

S relative to the surface

T transpose

t indicate a configuration at time t

t+D t indicate a configuration at time t+D t

– relative to the global coordinate system

.. indicate a second derivative in time

The Incremental Equilibrium Equation

Starting from the principle of virtual work, the incremental equilibrium equation for geometrically nonlinear static problems, using an Updated Lagrangean Formulation, is given by:

where iDS_ij are the increments of the second Piola-Kirchhoff stress tensor components referred to the configuration at time t, d_tDe_ij are the increments of the virtual Green-Lagrange strain tensor components referred to the configurations at time t (these increments of the virtual strain components are decomposed in the sum of increments of linear strain components d_tDe_ij and nonlinear strain components d_tDh_ij), ^tt _ij are the cartesian components of the Cauchy stress tensor at time t, and ^t+Dt are components of the external forces at time t + Dt given by:

with and being the components of the body force and the surface force, respectively; and du_i the components of virtual displacement vector at time t + Dt. In Eq.(1) and Eq.(2) V and A indicate, respectively, the domain and its boundary.

It is worthwhile to remember that the second Piola-Kirchhoff stress tensor components at time t + Dt referred to their configurations at time t are given by:

where is the determinant of the deformation gradient and

On the other hand, the components of the Green-Lagrange strain tensor in terms of displacement components (u) at time t + Dt, referred to the configuration at time t, are given by:

Equation (5) may be also written as follows:

with the linear and nonlinear terms given, respectively, by:

Linearizing the first term of the left hand side in Eq.(1), one obtains:

where the incremental constitutive relation was used and being ^tD_ijrs the components of the fourth order constitutive tensor containing material properties.

For dynamic problems it is necessary to add an inertia force to the equilibrium equation of motion, which may be written as follows:

where ü_i indicates de second derivative of the displacement components with respect to time and ^t+Dtr is the specific mass at time t + Dt.

The Triangular Finite Element for Thin Plates and Shells

Assuming the Kirchhoff theory, incremental displacement components may be expressed as follows:

where i = 1,2 represent directions corresponding to the local axis x and y, respectively; Du_i are displacement components in the mid plane; w and Dw are the vertical displacement and its respective increment; z is the coordinate value in the perpendicular direction to the mid plane, which is taken as the reference plane.

Substituting Eq.(10) in Eq.(7), it is obtained:

A typical triangular element with membrane and bending degrees of freedom, including drilling, is presented in Fig. 1.

Membrane displacement components field may be interpolated with respect to nodal membrane displacement components as follows:

where the vector of membrane nodal displacement components and the corresponding interpolation functions are given, respectively, by (Zhang, Lu and Kuang, 1998; Teixeira, 2001):

being

with L_i representing area coordinates and denoting by x_i and y_i the nodal coordinates. By cycling i ® j ® m, shape functions for can be obtained. The area coordinates are defined by:

being A the triangular element area, given by:

The displacement component, perpendicular to the element mid plane, may be interpolated with respect to the bending nodal displacement components as follows (Zhang, Lu and Kuang, 1998; Teixeira, 2001):

where the vector of bending nodal displacement components and corresponding interpolation functions are given, respectively, by:

being:

and:

with:

Increments of linear strain components given in Eq.(11) are separated in membrane and bending strain components and they may be expressed in terms of displacement components. Thus, it is obtained:

where the strain-displacement matrices (B) are:

with A being the element area; b_i and c_i are defined in Eq.(15), while L_i is an area coordinate. H_i, H_xi and H_yi are given by Eq.(20), Eq.(21) and Eq.(22), respectively.

The A element described in this section is a conforming element with the compatibility conditions being satisfied in each node and in each element side (see Teixeira, 2001). Using the drilling degree of freedom, numerical accuracy is improved and singularity of the stiffness matrix is avoided for coplanar elements. The total stiffness matrix for each element is obtained by superposition of the membrane and bending matrices.

The Updated Lagrangean Formulation

The incremental equation corresponding to the principle of virtual work in an Updated Lagrangean Formulation (ULF) was presented in Eq.(8). This expression may be written in a compact form as follows:

I₁ is given by

where D_m and D_b are the constitutive matrices for the membrane and bending effects, whereas D_mb = D_bm represents the constitutive matrix which results from membrane-bending coupling effect. These constitutive matrices, referred to the element local coordinate system, are defined by (Jones, 1999; Stegmann and Lind, 2001):

with r_g being the rotation matrix from the global to the local coordinate system, which is defined by:

In Eq.(33), g is the angle formed by the global axis x_g and the local axis x_l, as indicated in Fig. 2, where the fibers reference system is also shown.

In Eq.(30), Eq.(31) and Eq.(32), components of constitutive matrices are given by:

where n is the number of layers; z_k-1 and z_k are the coordinates normal to the lower and upper surfaces of layer k; are elastics constants of the layer k in the global coordinates system (see Fig. 2) defined by the following expressions:

being q_k the angle formed by the global axis x_g and the fibers local axis 1, as its is shown in Fig. 2.

Q_ij are elastic constants in the layer k in the fibers coordinates system and are defined by the following expressions (Jones, 1999):

where and are the elastic moduli of the layer k in the direction of the axis 1 and axis 2, respectively; is the shear modulus in the plane 1-2 of the layer k in the fiber coordinates system; is the Poisson coefficient defined as the relation between the strain in the transversal direction j and the axial strain in the direction i, considering the fiber coordinates system.

The nonlinear term l₂ is defined by the following expressions:

where ^tT contains the membrane internal forces and it is given by:

with

The internal work due to the external forces at time t + Dt, I₃, is given by:

Finally, the virtual work due to internal forces at time t, I₄, is defined by:

where ^tN_m, ^tM_b, (and ^tMN) are nodal vectors of internal forces corresponding to membrane, bending and membrane-bending coupling effects, respectively. These vectors are given by:

The Finite Element Incremental Equilibrium Equations for Static Problems

In static linear problem analysis, only the linear part of the stiffness matrix must be considered, and the finite element discretization at the element level leads to the following incremental matrix equation:

where the linear stiffness matrix, K_L, taking into account membrane (K_m), bending (K_b) and membrane-bending coupling effects (K_bm = ), are given by:

The nodal vectors of external loads referred to membrane and bending effects, respectively, are defined by the following expressions:

Equation (56) is obtained at time t + Dt.

Finally, the equivalent nodal vectors due to internal forces at time t corresponding to membrane and bending effects, respectively, taking into account membrane-bending coupling, are given by:

To calculate ^tN_m and ^tMN, given by Eq.(52) and Eq.(54), respectively, the displacement components due to membrane effects, contained in vector ^tu_m (at time t), are used. These displacement components are obtained performing the difference between the local coordinates at time t and at time 0 (Bathe and Ho, 1981).

To calculate the bending moment components at time t, the following expression is used:

where ^t+DtDu_b is the vector containing increments of displacement components due bending effects at time t + Dt.

The components of ^tNM, defined in Eq.(53), are determined in the same form as ^tM_b.

In a geometrically nonlinear static problem the nonlinear stiffness matrix ^tK_NL must be added to the linear part of the stiffness matrix K_L in Eq.(56). Matrix ^tK_NL is given by:

where: ^tT and G_G are defined by Eq.(48) and Eq.(49), respectively.

The incremental equilibrium equation, previously presented, is referred to each element local coordinates, and a transformation to a common global system is necessary in order to perform the assemblage procedure. Thus, considering the cosine of the angles formed by the local coordinates system (x^l, y^l, z^l) and the global coordinates system (x^g, y^g, z^g), they me be given by (see Fig. 3):

where:

being:

and then the transformation matrix from the global to the local coordinates system is given by:

The global equilibrium equations, obtained assemblying all finite elements forming the mesh, are solved using the Generalized Displacement Control Method (GDCM) as proposed by Yang and Shieh (1990), which is suitable for highly nonlinear problems with multiple critical points.

The Finite Element Incremental Equilibrium Equations for Dynamic Problems

The dynamic incremental iterative equilibrium equation, using Newmark’s Method, is given by (Bathe, 1996):

where the superposed index i is referred to the iteration number and

with

In Eq.(67), Eq.(69), Eq.(70) and Eq.(71) the coefficients are given by:

where Dt is the time step, while d and a are two parameters fixed by the user. In Eq.(67) and Eq.(68) Mm, C and F are the mass matrix, the damping matrix and the internal force vector, respectively. It was considered in this work C = 0.

When the convergence of the iterative process is obtained, acceleration, velocity and displacement vectors are updated for a specific time step with the following expressions:

with

To solve the system of Eq.(66), the conjugate gradient method is used, which has less memory requirements with respect to direct methods.

The consistent mass matrix for each element is given by:

where n is the number of layers, h^(k) and r^(k) are the thickness and the specific mass of the layer k, respectively, and H is defined by:

Numerical Applications

Geometrically Nonlinear Static Analysis of a Hinged Cylindrical Laminated Shell With a Concentrate Load at its Central Point

A hinged cylindrical laminate shell with a concentrated load at its central point is shown in Fig. 4. Its geometrical properties are: R(radius) = 2.54m, b(length) = 0.254m, h(total thickness) = 6.35x10^-3m or h = 12.70x10^-3m, j(angle) = 0.1rad. The load is applied at the shell center (point A) and is equal to P = 3000.00 N. Owing to symmetry, one quarter of the shell was modeled. The shell was analyzed for two different configurations: (a) with an isotropic material; (b) with a stacking sequence [90/0/90]. For the isotropic case the material properties are: E(Young modulus) = 3102.75MPa, G(shear modulus) = 1193.37MPa and v(Poisson coefficient) = 0.3. For an orthotropic material the properties are: E₁ = 3300.00MPa, E₂ = 1100.00MPa, G₁₂ = 660.00MPa and v¹² = 0.25. For h = 6.35mm a mesh with 50 triangular elements (generated in 5x5 = 25 rectangular regions) and nodes, and for h = 12.70mm a mesh with 32 triangular elements (generated in 4x4 = 16 rectangular regions) and 25 nodes were used. The boundary conditions are: u_x = Q_y = Q_z = 0.00 on the line , u_x = u_y= w = Q_x = Q_z = 0.00 on the line and u_y = Q_x = Q_z on the line .

Results were compared with those obtained by Sze, Liu and Lo. (2004) using 16x16 S4R elements, with 289 nodes, contained in the software ABAQUS (Hibbit et al., 1998) when h = 6.35mm and 8x8 S4R elements, with 81 nodes, when h = 12.70mm; Sze, Liu and Lo (2004) show that with these meshes very accurate results are obtained. Results for load x displacement for the isotropic and the orthotropic case, respectively, are presented in Fig. 5 and Fig. 6.

This problem, which is an interesting benchmark due to the snapping behavior of the structure, shows an excellent concordance with the results presented by Sze, Liu and Lo (2004), although in the present paper meshes are not so refined as in the reference.

Geometrically Nonlinear Static Analysis of a Pinched Cylindrical Laminated Shell

This application is similar to the example analyzed previously (Fig. 4) but now the stacking sequences are a cross ply [0/90/0] and an angle ply [45/-45]. Two types of boundary conditions are considered: an hinged and a clamped shell. Geometrical and material properties are the same as in the first example, but only the case where h(total thickness) = 12.70x10^-3m is studied. For the hinged shell a load P = 26.70kN is applied at its center point, while for the clamped shell a load P = 89.00kN is applied. A mesh with 72 triangular elements (generated in 6x6 = 36 rectangular regions) and 49 nodes was used.

Results are presented in Fig 7. and Fig 8. and compared with those obtained by Yeom and Lee (1989) and a good agreement was obtained. These authors used 36 degenerated solid elements with nine nodes and five degrees of freedom per node. In this case the total number of nodes is 85. In the shell with hinged edges softening-stiffening responses with inflexion points are observed for both stacking sequences, but this behavior is modified when hinged edges are substituted by clamped edges. It may be observed that results obtained in this paper are very similar to those of the reference with a mesh where a bi-quadratic degenerated quadrilateral element using Lagrangean polynomials as shape functions were substituted by two triangles with three nodes.

Geometrically Nonlinear Dynamic Analysis of a Spherical Laminated Shell Under an External Uniform Pressure

An hinged spherical laminate shell is presented in Fig. 9, where R(radius) = 10.00m, h(total thickness) = 10.00x10^-3m, b(length of the projected side) = 1.00m. The stacking sequence is an angle ply [45/-45] and the material properties are: E₁ = 250.00GPa, E₂ = 10.00GPa, G₁₂ = 5.00GPa, v₁₂ = 0.25 and r(specific mass)=1.00x10⁸kg/m³.Only one quarter of the shell was modeled using a mesh with 32 triangular elements (generated in 4x4 = 16 rectangular regions) with 25 nodes, taking into account symmetry conditions. The following boundary conditions were prescribed: u_y = Q_x = Q_z = 0.00 on the line , on the line , u_x = w = Q_y = Q_z = 0.00 on the line and u_x = Q_y = Q_z on the line . The adopted time step and uniform pressure acting on the external surface were, respectively, Dt = 3.00x10^-2s and q = 2000.00Pa.

The dynamic response at the central point A is shown in Fig. 10 and it is compared with the result presented by To and Wang (1998), and a good concordance is observed. In the present work the same response were obtained with the consistent and the lumped mass matrix.

To and Wang (1998) used the linear hybrid laminated composite triangular shell elements (HLCTS), with eighteen degrees of freedom and a first-order transverse shear deformation theory. These authors used a mesh with 36 triangular elements (generated dividing by four a region formed by 3x3 = 9 rectangles) with 25 nodes.

Teixeira (2001) used a modified consistent mass matrix taking into account only displacement degree of freedom, i. e. H in Eq.(78) was substituted by:

In both cases, in this specific example, results are close to those obtained by To and Wang (1998).

Geometrically Nonlinear Dynamic Analysis of a Cylindrical Laminated Shell Subjected to an Uniform Internal Pressure

A cylindrical shell subjected to an uniform internal pressure is presented in Fig. 11, were the geometric properties are: R(radius) = 2.54m, a(arc length) = 0.508m, b(length) = 0.508m, h(total thickness) = 1.27x10^-3m and j(angle) = 0.10rad. Material properties are: E₁ = 137.90GPa, E₂ = 9.86GPa, G₁₂ = 5.24GPa, v₁₂ = 0.30 and r = 1562.20kg/m³. Two stacking sequences were considered: (a) a symmetric stacking sequence with eight layers [0/-45/90/45]s and (b) a cross ply [0/90]. In both cases layers with the same thickness were adopted. It is considered that the shell boundary is clamped (u_x= u_y = w = Q_x = Q_y= Q_z = 0.00). The adopted time step and uniform internal pressure were, respectively, Dt = 5.00x10^-2ms and q = 6895.00Pa. A mesh with 128 elements (generated in 8x8 = 64 rectangular regions) was used.

Results for the displacement component w (in the direction of the axis z_g) are presented in Fig. 12 and Fig. 13 and comparison with those presented by To and Wang (1999) and Wu, Yang and Saigal (1987) are also shown. To and Wang (1999) used the HLCTS element while Wu, Yang and Saigal (1987) used a curved high-order quadrilateral shell element.

A good relatively agreement is observed, especially with respect to the response obtained by To and Wang (1999). These authors used a mesh with 36 triangular elements (originated dividing in four a region formed by 9 rectangles) with 25 nodes. Some differences may be observed with respect to the results (especially for amplitudes) presented by Wu, Yang and Saigal (1987). In this work a very refined mesh was used in order to confirm results obtained by To and Wang (1999) with respect to those presented by Wu, Yang and Saigal (1987).

Final Remarks

An efficient triangular finite element for the geometrically nonlinear static and dynamic analysis of laminate shells was presented in this work. The mass and the non linear stiffness matrix may be explicitly implemented. Results of all the examples show a very good agreement with those presented by others authors using different types of elements and formulations. In most of the examples the meshes employed in this work are less refined with respect to those used in the corresponding references.

By definition, stress components s_zz, t_xz and t_yz are considered with a null value in the Classical Lamination Theory. Actually these stress components are not null and they may assume important values in the layer interfaces causing delamination. Physically, one of the reasons of this failure mode is the fact that s_zz, t_xz and t_yz change significantly from one layer to the other due to the abrupt differences of elastic properties. Thus, this theory is not suitable to study delamination.

In future works topics such as shells optimization analysis and shells with smart materials will be studied.

Acknowledgments

The authors wish to thank CNPq (Brazilian National Research Council) for its financial support.

Paper accepted February, 2008.

Technical Editor: Nestor A. Zouain Pereira.

• Bathe, K-J., 1996, "Finite Element Procedures", Ed. Prentice-Hall, New Jersey, United States of America, 1037 p.
• Bathe, K-J. and Ho, L-W., 1981, "A Simple and Effective Element for Analysis of General Shell Structures", Computers and Structures, Vol. 13, pp. 673-681.
• Hibbit, Karlsson and Sorensen Inc., 1998, ABAQUS user's Manual, version 5.8.
• Jones, R. M., 1999, "Mechanics of Composite Materials" Ed. Taylor & Francis, Philadelphia, United States of America, 519 p.
• Stegmann, J. and Lund, E., 2001, "Notes on Structural Analysis of Composite Shell Structures", Aalborg University, Aalborg, 90 p.
• Sze, K. Y., Liu, X. H. and Lo, S. H., 2004, "Popular Benchmark Problems for Geometric Nonlinear Analysis of Shells", Finite Elements in Analysis and Design , Vol. 40, pp. 1551-1569.
• Teixeira, P. R. F. 2001, "Numerical Simulation of the Interaction of Three Dimensional Compressible and Incompressible Flows with Elastic Structures using the Finite Element Method" (In Portuguese), Ph. D. Thesis, Federal University of Rio Grande do Sul, R.S., Brazil, 237 p.
• To, C. W. S. and Wang, B., 1998, "Transient Responses of Geometrically Nonlinear Laminated Composite Shell Structures", Finite Elements in Analysis and Design, Vol. 31, pp. 117-134.
• To, C. W. S. and Wang, B., 1999, "Hybrid Strain Based Geometrically Nonlinear Laminated Composite Triangular Shell Finite Elements", Finite Elements in Analysis and Design, Vol. 33, pp. 83-124.
• Wu, C. Y., Yang, T. Y., Saigal, S., 1987, "Free and Forced Nonlinear Dynamics of Composite Shell Structures", Journal of Composite Materials, Vol. 21, pp. 898-909.
• Yang, Y-B. and Shieh, M-S., 1990, "Solution Method for Nonlinear Problems with Multiple Critical Points, AIAA Journal, Vol. 28, pp. 2110-2116.
• Yeom, C. H., and Lee, S. W., 1989, "An Assumed Strain Finite Element Model for Large Deflection Composite Shells", International Journal for Numerical Methods in Engineering, Vol. 28, pp. 1749-1768.
• Zhang, Q., Lu, M. and Kuang, W., 1998, "Geometric Non-linear Analysis of Space Shell Structures Using Generalized Conforming Flat Shell Elements – for Space Shell Structures", Communications in Numerical Methods in Engineering, Vol. 14, pp. 941-957.

Publication Dates

• Publication in this collection
  25 Apr 2008
• Date of issue
  Mar 2008

History

• Accepted
  Feb 2008
• Received
  0000