Accurate finite elements based on shear deformation theory for the analysis of laminated composite plates

 

 

V.E. Verijenko^I; S. Adali^II; E.B. Summers^III

^IAssociate Professor, Department of Mechanical Engineering, University of Natal, King George V Avenue, Durban, 4001 Republic of South Africa
^IIProfessor of Solid Mechanics, Department of Mechanical Engineering, University of Natal, King George V Avenue, Durban, 4001 Republic of South Africa
^IIIPhD candidate, Department of Mechanical Engineering, University of Natal, King George V Avenue, Durban, 4001 Republic of South Africa

 

 

---

ABSTRACT

A finite element formulation for the analysis of laminated composite plates based on a higher-order theory is presented. This formulation leads into a discrete-continuous scheme where the surface of the laminate is discretized with each finite element forming a heterogeneous continuum through the thickness. Rectangular and triangular finite elements are formulated. The degrees of freedom of the nodal points of these elements are independent of the number of layers. Nonlinear laws governing the variation of the components of the displacement vector and of the stress and strain tensors through the thickness of layers are taken into account. The laminate may also exhibit heterogeneous properties in the plane of the plate, where elements with different properties are used as an approximation. The elements are applied to the bending of laminated plates with various loading and boundary conditions and numerical results are obtained. The solutions presented are compared with those obtained using the three-dimensional elasticity theory, and with the closed form solutions of other authors. It is shown that the present approach reduces the number of unknown variables and broadens the field of application of the finite element method.

---

 

 

1 Introduction

Fibre-reinforced composites are now widely used in many engineering applications for their exceptional strength and stiffness properties relative to their weight. Therefore it is important to be able to accurately model the behaviour of structural components such as plates manufactured from these advanced materials. It is well known that the use of classical plate theory, which is based on the Kirchoff assumptions, leads to intolerable errors in the analysis of composite structures. If the phenomenon of transverse shear is neglected, even the deflection of thin laminated plates are underpredicted when the laminae differ significantly in their elastic properties.

 

Nomenclature

U_i, wDisplacement

ϵ_ij, K_ijDeformation of the reference surface

e_ijComponents of the strain tensor

σ_ijComponents of the stress tensor

χShear function

φ_kDistribution functions for shear deformation

ψ_k Distribution functions for tangential displacement

EModulus of elasticity

GShear modulus

VPoisson's ratio

QShear forces

MBending moments

DStiffness characteristics

η, ξ Dimensionless coordinates

KStiffness matrix

v, Degrees of freedom of the nodal points

RVector of nodal forces

Numerous approaches have been suggested which take into account the three-dimensional stress and strain behaviour of multilayered plates based on two-dimensional higher-order theories. Details of different higher-order theories and their finite element modelling may be found in reviews [1-5]. However, as mentioned in [6] and in the review [2], some of the theories exhibit no compatibility between the nonlinear kinematic model, which considers the distortion of the normal, and the system of internal forces and moments equivalent to those obtained using the 'straight line' hypothesis.

A higher-order theory of laminated plates and shells without these drawbacks has been formulated by the authors and is presented in [7-9]. This theory considers plates with transversely isotropic layers of different thicknesses, stiffnesses, and densities, in which the number and sequence of layers are arbitrary. The physical and mechanical characteristics of each layer are variable through the thickness and the layers are assumed to be perfectly bonded. The equations for the tangential components of the displacement vector and the stress and strain tensors consist of similar terms which separately take into account the states of pure bending and transverse shear. This important feature enables the efficient numerical application of this theory using an independent but analogous approximation of the components of the displacement vector.

 

2 Rectangular element of laminated plate

In this section, a rectangular element with transversely isotropic layers of different thicknesses, which takes into account shear deformation, is formulated. The basic equations are expressed for a model which includes transverse shear under the normal loading p₃. Since there is no shear load on the external surfaces of the plate, the tangential displacements of the reference surface are assumed to be negligible. The plate has dimensions of a xb and a total thickness of a₀ + a_n as shown in Figure 1a. The reference surface is positioned such that the maximum transverse shear occurs at this surface.[9] These assumptions are acceptable in many practical engineering applications.[10]

 

 

In the following analysis a subscript after a comma denotes differentiation with respect to the variable following the comma, and k refers to the k-thlayer. The displacements of the plate in the x₁, X₂ and z directions are denoted by u₁, u₂ and w, respectively.

The basic equations for the shear-deformable model are

where is a new unknown function which is called the 'shear function'.[11]

The graphical interpretation of this kinematic model is given in Figure 2. The distribution functions of the tangential displacements Ψ_k and shear deformations φ_k through the thickness of layers are expressed as

 

 

where are the modulus of elasticity, Poisson's ratio and shear modulus, respectively, in the plane of isotropy and are tnese characteristics in transverse direction.

The corresponding static model is given by

The variation of potential energy of the r-th element now may be expressed as

where the generalized displacement = is the r-th element. The forces and the moments are defined as

where the following stiffness characteristics [10] are required

The variation of the external load is

The derivation of the rectangular finite element is now given for multilayered plates. Six degrees of freedom are assigned to each nodal point (see Figure 1b). The first three degrees of freedom correspond to the deformation caused by the bending and are the deflection w and the two angles of rotation αand ß about the axes . The remaining three degrees of freedom correspond to the shear deformation which is analogous to deflection and to which are analogous to angles of rotation. Thus at each nodal point m = 1, 2, 3, 4 of the element r in the local coordinate system, two groups of displacements are defined by

which relate to the deformation of the reference surface.

The approximation of the displacement in the region of the FE may be introduced in the well-known form [12]

which may be written as

where the system of approximation functions

are given by

where are dimensionless coordinates.

The stiffness matrix for element r may be obtained from equations (5) using (12), and may be expressed as

where each block is symmetric. The submatrix K₁₁corresponds to the state of bending and is identical to the stiffness matrix for a homogeneous plate [12]. The sub-matrix K₂₂ corresponds to the state of transverse shear, and the submatrices K₁₂, K₂₁ characterize the interaction of these states. The dimensions of the submatrices in (13) are 12 × 12, and therefore the dimension of the matrix (13) is 24 × 24. The entries of the submatrix K₁₁are given in Table 1 and Appendix A. The submatrix K₁₂is similar to K₁₁where the entries are equal to the corresponding entries of K₁₁when is replaced with . In this case the following stiffnesses are required

 

 

The entries of the submatrix K₂₂can be determined in terms of the entries k_ltas

and α_lt are additional quantities used in the calculation of the coefficients of shear (see Appendix B). It is noted that v₁₁ has been replaced with in the coefficients k_lt of stiffness matrix, where

If the load p₃ = q is uniform over the finite element, then the vector of nodal forces of the element r is given by

In the case of a point load applied in the centre of the rectangular element, the vector of equivalent forces is given by

As seen in equation (8), the load p₃ produces work over the displacements w. Therefore the components of the vector of nodal forces corresponding to the displacements equal zero.

It is noted that the above element may also model the external loads by moments at the nodal points, and this improves the accuracy of the analysis.

The theory governing the finite elements requires two groups of constraints at the boundary of the plate. The first group, referred to as the 'external' boundary conditions, constrains the two-dimensional reference surface of plate and models the general type of support for the plate. The second group, referred to as the 'internal' boundary conditions, models the transverse deformations through the thickness at the boundary of the plate. These boundary conditions are modelled in terms of the degrees of freedom. Table 2 gives the details of the boundary conditions and the corresponding constraints.

 

 

3 Triangular element of laminated plate

This case differs only in the geometry of the finite element (see Figure 3a). At each nodal point m (m = 1,2,3) the group of displacements may be expressed as

 

 

A fourth order polynomial is used to approximate the continuous displacement field.

The dimensionless coordinates

are introduced in order to simplify the formulation of the stiffness matrix (see Figure 3b).

In the new ξ, ηcoordinate system the vectors of the degrees of freedom of the nodal point m are given by

For each of these vectors there are nine approximation functions given by

These functions may be expressed in explicit form as

The system of approximation functions

which are compatible with the degrees of freedom

coordinate system are given as

Then the displacements in the finite element region may be written as

Using equations (5), (23), and (24), the stiffness matrix of a triangular element of a laminated plate may be obtained explicitly. This matrix may also be represented in block form as given in equation (13).

 

4 Numerical results and discussion

The accuracy of the finite elements developed in sections 2 and 3 is now illustrated by means of numerical examples.

Example 1. As a first example, a three-layered square plate is considered with the dimensions 2a x 2a and under a uniform load p₃ = q = 10⁵Nm^-2. The 'external' boundary conditions are taken as simply supported and the 'internal' boundary conditions are taken as flexible out of plane of the edge, but rigid in this plane. The properties of the external (bearing) layers 1 and 3 are h₁ = h₃ = 1 × 10^-3m, E₁ = E₃ = 6.8 ×10⁴MPa, G₁ = G₃ = 26 150MPa, and v_x = v₃ = 0.3. The properties of the, filler layer k = 2 are h₂ = 15 × 10^-3 m, E₂ = 4 800MPa, G₂ = 380 MPa, and v₂ = 0.3. The plate is divided into four rectangular finite elements. Noting the symmetry, the constraints of the nodal displacements may be expressed as

The general vector of the nodal displacements is

Taking into account equation (17), the general vector of corresponding nodal forces is expressed as

The entries of the stiffness matrix [K] may be determined using Table 1, Appendix A, Appendix B, and equation (14) where D₀₀ = 11.06Nm, c_2l = 17.44× 10^-1m², c₂₂ = 17.48 × 10^-4m², and v₁₁ = v₁₂ = v₂₂ = v = 0.3. The system of algebraic equations may be written as

Using the block principle [8] the solution is divided into two simpler systems

The solution gives the deflection at the centre of the plate as w₄ = 3.56 × 10^-3m. The analytical solution for this case gives w_max = 3.05 × 10^-3m. When the number of finite elements along each side equals 4, 6, and 8, the corresponding maximum deflections of the plate are obtained as w_max = (3.19, 3.11,3.08) 10^-3m, respectively. The solution using Reissner's theory gives w_max = 3.16 × 10^-3m whereas the classical theory gives w_max = 2.16 × 10"³m which is obviously inaccurate.

Example 2. Consider the bending of a square simply supported homogeneous plate of dimensions a× a subjected to sinusoidal loading q sin(πx₁/a) sin(πx₁/a). The finite element results are compared with the exact solution given in [13]. The application of the finite elements developed in the previous sections produced accurate results even for a plate with a thickness ratio of h/a = l/5. In this case the number of unknown is 3.5 times less than that of the case where 3-D elements are used (Table 3). Table 3 also gives the results for a three-layered symmetrical plate subjected to a sinusoidal loading. The thickness of the bearing layers is taken as h₁ = h₃ = h/6 and the thickness of the filler layer as h₂ = (4/6) h. The layers are isotropic with the elastic properties E₁ = E₃ = 10³E₂, v₁ = v₂ = 0.3. The boundary conditions are the same as in the first example. The thickness ratio h/a equals 1/10. The discrete-continuous scheme (DCS) solution is compared with the analytical solution [14] and the 3-D finite element solution. Accurate results are obtained with the number of unknowns being approximately 10 times less than that of the case of 3-D elements as indicated in Table 3.

 

 

Example 3. Consider a three-layered square plate supported at each of the four corners with the following characteristics: dimensions a = 25h = 1 m; total thickness of the laminate h = 40 ×10^-3 m where the thickness of the bearing layers h₁ = h₃ = 2×l0^-3m, and the thickness of the filler layer is h₂ = 36 × 10^-3 m; elastic moduli E₁ = E₃ = 7 × 10⁴ MPa, E₂ = 70 MPa; Poisson's ratios v₁ = v₂ = v₃= 0.3.

The maximum deflections and stresses in the bearing layers are given in Table 4 for the various types of 'internal' boundary conditions and loading (see Table 2). The types of internal boundary conditions listed are: (1) rigid in and out of the plane of the edge; (2) flexible out of the plane of the edge but rigid in the plane: (3) flexible in the plane but rigid out of the plane; (4) rigid at the corners with no constraints along the sides; and (5) no internal constraints.

 

 

The results for the deflections indicate three different types of boundary behaviour: firstly, for the first and second type of constraints, when the shear in the plane of the edges vanishes and the deflections are almost identical; secondly, for the third and fourth type of constraints, shear is allowed in the plane of the edge and the deflections increase by a factor of 1.5 to 2.2 in comparison with the first case for point and uniform loads; and thirdly, when there are no internal constraints to resist the shear, the deflections increase by a factor of 1.7 to 2.8 in comparison with the first case for the same type of loading. The results obtained for normal stresses at the centre of the plate indicate that the type of internal constraints have only a minor influence on these stresses in comparison with the type of support. The stresses on the edge depend more on the internal constraints than do the stresses away from the edge, and are 1.2 times greater for the cases 3-5 than those for the case 1.

Example 4. Consider the analysis of the plate in Example 1 using triangular elements. Results are given for various numbers of elements in Table 5, where the deflection is given in the form and the stresses in the external layers are given in the form . Deflections are also determined from the analytical solutions of various theories. Reissner's theory gives ; the classical theory gives . Table 5 shows that the triangular elements give accurate results.

 

 

5 Summary and conclusions

The finite element formulation presented above leads to a discrete-continuous scheme for the analysis of laminated composite plates where each finite element forms a heterogeneous continuum through the thickness. On the basis of this scheme, rectangular and triangular finite elements are developed which take into account the deformation of transverse shear. Moreover, the degrees of freedom of the nodal points of these elements are independent of the number of layers.

The approximations and degrees of freedom related to the different types of stress and strain states of the plate yield similar coefficients in the stiffness matrix and in the blocks comprising the stiffness matrix. A significant number of coefficients are the same as those of the classical theory, which simplifies the calculation of the stiffness matrix and allows the experience gained in using the finite element method in similar applications on the basis of classical theory to be extended to a non-classical approach based on the present higher-order theory.

The accuracy of the elements is illustrated by comparing the results with exact, analytical and finite element solutions of other authors. The results predicted by the approach presented here are found to be in excellent agreement with three-dimensional elasticity solutions.

The present investigation indicates that the finite elements proposed in this study are highly efficient and accurate, and may easily be incorporated into existing finite element packages.

 

Acknowledgement

The authors gratefully acknowledge the support of the Foundation for Research Development through a Core Programmes research grant.

 

References

[1] Dudchenko AA, Lur'e SA & Obraztsov IF. Anisotropic multilayer plates and shells. Summary of Science and Engineering. Mechanica Deformiruemogo Tela, 1983, 15, pp.3-68 (in Russian).         [ Links ]

[2] Noor AK & Burton WS. Assessment of computational models for multilayered composite shells'. Appl. Mech. Rev., 1990, 43(4), pp.67-97.         [ Links ]

[3] Reddy JN. A review of refined theories of laminated composite plates. Shock Vib. Dig., 1990, 22(7), pp. 3-17.         [ Links ]

[4] Reddy JN. On refined theories of composite laminates. Meccamca, 1990, 25(4), pp.230-238.         [ Links ]

[5] Librescu L & Reddy JN. A few remarks concerning several refined theories of anisotropic composite laminated plates. Int. J. Eng. Set., 1989, 27(5), pp.515-527.         [ Links ]

[6] Piskunov VG, Verijenko VE &Prisyazhnyuk VK. Calculation of inhomogeneous shells and plates using finite element methods. 1987, Veisha Shkola, Kiev, 200pp. (in Russian).

[7] Verijenko VE. Nonclassical theory of elasto-plastic strain of laminated transversely isotropic tapered shells. Strengths of Materials, 1987, 19(4), pp.547-553.         [ Links ]

[8] Verijenko VE, Adali S & Summers EB. Discrete-continuous scheme of FEM for analysing laminated composite shells and plates. Proceedings of FEMSA 92, 1992, pp.589-600.         [ Links ]

[9] Piskunov VG, Verijenko VE, Adali S &Summers EB. A higher order theory for the analysis of laminated plates and shells with shear and normal deformation. Int. J. Eng. Sci., 1993, 31(6), pp. 967-988.         [ Links ]

[10] Piskunov VG & Verijenko VE. Linear and non-linear problems in the analysis of laminated structures. 1986, Budivelnik, Kiev, 176pp. (in Russian).

[11] Verijenko VE & Prisyazhnyuk VK. Linear and non-linear models of contact interaction of laminated constructions with elastic basement. International Congress IKM-X Weimar Germany, 1984, 6, pp.175-182.         [ Links ]

[12] Clough RW & Penzien J. Dynamics of structures. 1975, McGraw-Hill, New York.

[13] Vlasov BF. On the bending of a rectangular thick plate. Vest. Moscow University, 1957, 2, pp.22-34 (in Russian).         [ Links ]

[14] Brukker, LE. Some variants of the equations for three-layered plates. Journal of Analysis of Elements of Aviation Structures, 1965, 3, pp.74-99 (in Russian).         [ Links ]

 

 

First received July 1992
Final version May 1993

 

 

Appendix A

Elements of submatrix K₁₁ for rectangular FE

 

Appendix B

Additional quantities for computation of the sub-matrix K₂₂ for rectangular FE