Abstract

An efficient C⁰ FE model for the analysis of composites and sandwich laminates with general layup

S.K. Singh^I,* * Author email: sushilbit@yahoo.co.in ; A. Chakrabarti^I; P. Bera^II; J.S.D. Sony^I

^IDepartment of Civil Engineering, Indian Institute of Technology, Roorkee-247 667 – India

^IIDepartment of Mathematics, Indian Institute of Technology, Roorkee-247 667 – India

ABSTRACT

A C⁰ continuous finite element model is developed to model the refined higher order shear deformation theory. The proposed element is an upgraded version of an element based on higher order shear deformation theory. The C⁰ continuity of the present element is compensated in the stiffness matrix calculations. The computational efficiency is achieved by the C⁰ continuous finite element model by satisfying the inter-laminar shear stress continuity at the interfaces and zero transverse shear stress conditions at plate top and bottom. The performance of the upgraded element is illustrated with many numerical examples.

Keywords: Finite Element, interlaminar stress continuity, composites, sandwich.

1 INTRODUCTION

The composite materials are widely used in civil, aerospace and other engineering fields due to their advantage of high stiffness and strength to weight ratio. Laminated composite structures are weak in shear due to their low shear modulus compared to extensional rigidity. Thus the effect of shear deformation is quite significant, making it vulnerable to failure. The Classical Plate Theory [26] under-predicts displacements and over-predicts the natural frequencies and the buckling loads [28]. However this kind of approach is not sufficient for laminated plate due to neglecting the transverse shear deformation in the laminates. In this context a number of plate theories have been developed where the major emphasis is to model the shear deformation in a refined manner.

These plate theories can be divided into two groups and they are (1) single-layer plate theory and (2) layer-wise plate theory. In single-layer theory [19,20,25,26,28,31,32] the deformation of the plate is expressed in terms of unknown parameters of a single plane, which is usually taken as the middle plane of the plate. These are similar to Reissner-Mindlin's plate theory (i.e., the first-order shear deformation theory - FSDT) which requires shear correction factor but there are some improvements, which allow the warping of plate sections to have a higher-order variation of transverse shear stresses/strains along the thickness. In the layer wise theory [18,27,29,30] the deformation of the plate is expressed in terms of unknowns of a number of planes, which are taken at the layer interfaces and also at some intermediate levels in some cases. The mathematical involvement in these plate theories is quite heavy and the solution becomes quite expensive in a multilayered plate, as the unknowns are dependent on number of layers. There is another class of layer wise plate theories [1,6,7,9,14,15,22] where the unknowns of different planes are expressed in terms of those of a particular plane using the condition of shear stress continuity at the layer interfaces. The number of unknowns is reduced in these plate theories considerably [1,6,7,9,14,15,22] .

In this context, the first-order shear deformation theory [32] maybe considered as the simplest option where an arbitrary shear correction factor is used since the transverse shear strain is assumed to have uniform variation over the entire plate thickness. The first order shear deformation theory which assumes a constant transverse shear strain across the thickness direction and a shear correction factor is introduced to correct the discrepancy between the actual transverse shear stress distribution and those assumed in this theory. The performance of first-order shear deformation theory is strongly dependent on shear correction factors [31]. For a better representation of the transverse shear deformations, higher order plate theories (HSDT) are proposed by Lo et al. [19], Manjunatha et al. [20], Reddy [25] and a few others, in which the use of shear correction factor could be eliminated. It gives continuous variation of transverse shear strain across the entire thickness, which leads to discontinuity in the variation of the transverse shear stresses at the layer interfaces. But the actual behavior of laminated plate is the opposite i.e., the transverse shear stress is continuous at the interfaces whereas the strains may be discontinuous. Moreover, the degree of discontinuity in the transverse shear strain is severe especially for sandwich plates due to a wide variation in their material properties.

In order to overcome the above problem, Srinivas [29], Toledano et al. [30], Li et al. [18], Robbins et al. [27], and some other investigators proposed layer-wise plate theories taking unknowns at each layer interface. These plate theories perform well but they require significant computational involvement in analyzing a multi-layered plate since the number of unknowns increases with the number of layers. A major development in this direction is due to Di Sciuva [14], Murakami [22], Liu et al. [9], and few others. They proposed zigzag plate theory where layer-wise theory is initially used to represent the in-plane displacements having piecewise linear variation across the thickness. The unknowns at the different interfaces are subsequently expressed in terms of those at the reference plane through satisfaction of transverse shear stress continuity at the layer interfaces. A further improvement in this direction is due to Di Sciuva [15], Bhaskar et al. [1], Cho et al. [6,7] and some other investigators who considered the variation of in-plane displacements to be a superposition of a piecewise linearly varying field on an overall higher order variation. Carrera [2] and Demasi [11] considered higher order terms in the displacement field, using Mukarmi's [21] zig-zag function and assumptions for transverse stresses brings about a large number of solution variables. However applying static condensation technique allows to eliminate the unknowns related to the transverse stresses and thus, to derive efficient plate theories [12,13]. Cho et al. [5] have also presented coupled zigzag theory for hybrid plates under thermoelectric load considering global variation for the deflection and a layerwise linear variation of electric potential, which are inadequate to capture the layerwise distribution of deflection due to thermal and potential fields. Kapuria et al. [16,17] have presented zigzag theory for hybrid beams and plates in which number of variable are reduced to FSDT by satisfying interface and boundary conditions, it yield approximately accurate results for cross ply only. Wu et al. [33] proposed C⁰ type higher-order theory for bending analysis of laminated composite and sandwich plates. However in article [33], the analytical formulations and solutions only presented for thermo-mechanical bending analysis of laminated composite and sandwich plates. Wu et al. [34] also proposed C⁰ type finite element based higher-order theory for accurately predicting natural frequencies of sandwich plate with soft core. These theories are usually referred as refined higher order shear deformation theory (RHSDT). However, there are very few C⁰ elements reported in the literature which can model the RHSDT.

Considering all these aspects in view, an attempt has been made in this study to develop an improved FE plate model to accurately predict the deflections and stresses of laminated composites and sandwich plates due to different loadings, boundary and geometric conditions. The plate model has been implemented with a computationally efficient C⁰ finite element based upon refined higher order shear deformation theory. The C⁰ element proposed by Shankara et al. [28] for simple higher order theory is upgraded to model the refined higher order shear deformation theory in the present study. The C⁰ continuity of the present element has been compensated in stiffness matrix calculations. The accuracy of the proposed eight-noded C⁰ element is established by comparing the results with three dimensional elasticity and other finite element solutions.

2 MATHEMATICAL MODEL

The in-plane displacement fields (Figure 1) are typical to those of RHSDT Cho et al. [7] and are as below:

where the subscript α represents the co-ordinate directions [α = 1, 2]

where denotes the in-plane displacements (i.e., u₀ along x-axis and v0 along y-axis) of any point on the mid surface, n_u and n_l are number of upper and lower layers respectively, are the slopes of k-th layer corresponding to upper and lower layers respectively, ξ_α, ϕ_α are the higher order unknown and H(z-z_k) and H(-z+ρ_k) are unit step functions.

The stress-strain relationship of a lamina say k-th lamina having any fiber orientation with respect to structural axes system (x-y-z) may be expressed as

The rigidity matrix [] can be evaluated by material properties and fibre orientation following usual techniques for laminated composites [10].

Now by utilizing the transverse shear stress free boundary condition at the top and bottom of the plate, (σ_3α |_{z = ±h/2} = 0) the components ξ_α and ϕ_α could be expressed as:

Substituting equation (4) and (5) in equation (1), the following expressions may be obtained:

Similarly by imposing the transverse shear stress continuity conditions at the layer interfaces the following expressions for S_α and T_α are obtained as below:

where a_α^k, b_α^k, c_α^k, d_α^k are constants depending on material properties of the individual layers and r is the respective layer interface.

By using equations (2)-(6) the strain field vector can be evaluated by

where {} is the strain field vector and {ε} is the modified strain vector at the reference plane, that is at the mid plane.

where the [H] matrix consists of terms containing z and some term related to material properties.

The nodal unknowns for the present FE model can be defined with respective to the co-ordinate system (x-y-z) as in Figure 1.

The actual displacement fields require C¹ continuity of the transverse displacement for the finite element implementation. In order to avoid the usual difficulties associated with C¹ continuity requirement, is replaced with independent variable w₁ and with w₂. For the sake of convenience to represent all variables as C⁰ continuous, the derivative of w with respect to x and y are expressed as follows

Now the potential energy (V) of the plate/laminate under the action of transverse load of intensity q may be written as

Where,

Using equation (9) the penalty term is expressed as:

where γ is the penalty parameter.

3 FINITE ELEMENT FORMULATION

For the present study, an eight noded quadrilateral C⁰ continuous isoparametric element has been used. The typical node numbering adopted is as shown in figure 2. The element has an arbitrary rectangular geometry in x-y coordinate system. To have a regular rectangular geometry the element is mapped to ξ-η plane. The generalized field variable (u) and the element geometry (x,y) at any point within the element are expressed in terms of nodal variables as follows.

where Ni is the shape function of the associated node.

Using equation (10) the strain vector in equation (9) can be expressed in term of {δ}_56×1 containing nodal degrees of freedom as

where [B]_17×56 is the strain displacement matrix.

The strain displacement matrix [B]_56×56 is evaluated by using the modified strain vector and the shape functions.

The transverse displacement may now be expressed in terms of the nodal displacement vector as

where the [N] matrix contain shape functions for the eight node element [8].

Also

where [P_x] and [P_y] matrix contains shape functions and their derivatives.

By substituting the expressions equations (15), (16) and (17) in equation (12) and minimizing the expression of the potential energy with respect {δ} equilibrium can be expressed as

where [K] is the element stiffness matrix and {P}the nodal load vector. The integration involved in the above expressions is carried out numerically by following the Guass quadrature method.

4 NUMERICAL RESULTS AND DISCUSSIONS

In order to demonstrate the accuracy and applicability of the present element a number of numerical examples on composites and sandwich laminates are solved by the proposed finite element model based on refined higher order shear deformation theory. The general geometric details of the plate problem considered for the present analysis is shown in Figure 3. The results obtained are presented in the form of different tables and figures. Initially a composite plate problem is solved by using the proposed finite element to model the higher order shear deformation theory (HSDT) with the inclusion of penalty parameters for study of convergence and also to compare the present results with published results. Finally, the proposed FE model is used to generate results based on refined higher order shear deformation theory (RHSDT). A large number of results are compared with the published results.

4.1 Analysis of laminated plates based on higher order shear deformation theory

In this section numerical examples on composite plates are solved by using the present FE model based on HSDT and penalty parameters. The present results are compared with those obtained by Sankara et al. [28].

4.1.1 Cross-ply square laminate subjected to uniformly distributed load

The problem of a three ply (0/90/0) square laminate, simply supported at all the edges and subjected to uniformly distributed load, is studied for different thickness ratio (h/a) ranging from 0.01 to 0.5. The plate is analyzed with different mesh divisions and the deflection obtained at the plate centre is presented with the analytical solution of Reddy [25] and Chakrabarti et al. [3] in Table 1. The result obtained shows the element is capable of predicting results with sufficient accuracy. In this it is observed that the first order shear deformation theory gives lesser deflections as compared to higher order shear deformation theory.

Although the individual layers possess different orientations but they have equal thickness and material property (E1= E2 = 25; G12 = G13 = 0.5E2, G23 = 0.2E2, n12 = 0.25 and n13 = 0.01), which is also applicable to all the subsequent problems.

4.1.2 Cross-ply rectangular laminate subjected to distributed load of sinusoidal variation

The plate is modeled by the proposed element using the mesh arrangement as shown in Figure 3. This mesh arrangement is used in all the examples where the plate is rectangular. A rectangular laminate subjected to a load having a distribution of q(x, y) = q0sin(πx/a)sin(πy/b) is considered in this example. The stacking sequence and boundary conditions are identical to those taken in the previous example. The study is made by taking thickness ratio (h/a) as 0.5, 0.1 and 0.01.

In all the cases the analysis is done with three different mesh divisions and the deflection and stress components obtained at the important locations are presented with the analytical solution of Reddy [25] and numerical solution of Chakrabarti et al. [3] in Table 2. The deflection and stress results of the element are exactly in agreement with the published results.

4.2 Analysis of composites and sandwich laminates based on refined higher order shear deformation theory

In the previous section the results based on higher order shear deformation show a definite improvement over that of first order shear deformation theory, but it requires further improvement specifically for the stresses evaluation to have a comparable value with 3-D elasticity solution Pagano [23]. This improvement is shown in this section by solving some problems of composites and sandwich laminates with the proposed element based on refined higher order shear deformation theory by compensating C⁰ continuity in stiffness matrix calculations.

4.2.1 Simply supported cross-ply square laminate

In order to show the performance of different plate theories a cross-ply (0/90/0) laminate (Figure 3, a = b) subjected to a uniformly distributed load of intensity q(x, y) = q0 sin (πx/a) sin (πy/b) is considered in this example. The full plate is analyzed with different mesh divisions. The results obtained are compared with those obtained by First order Shear Deformation Theory (FSDT), HSDT, RFSDT (Refined FSDT) and 3D elasticity solution. Table 3 shows the present values of central deflections = 100wE₂h³ /(q₀a⁴ ) based on RHSDT are close to the 3D elasticity solution while the FSDT gives the worst results.

4.2.2 Simply supported cross-ply square laminate

The problem of a cross-ply (0/90/90/0) laminated plate (Figure 3, a = b) subjected to a distributed load of intensity of q(x, y) = q0 sin (px/a) sin (py/b) is taken in this example. The variation of normal stress (_x = σ_xh² /(q₀a² )) across the thickness of the plate is presented in Figure 4 with those obtained by 3D elasticity solution [8]. The potential of the element is clearly reflected in the accuracy exhibited in the prediction of stress distribution. The in-plane shear stress, = t_xyh² /(qa² ) at the corner of the plate and the transverse shear stress = τ_xzh/(q₀a) at the centre of the edge of the plate is obtained by RHSDT and is plotted with 3D elasticity solution [23] in Figure 5 and Figure 6 respectively, which shows very good agreement between the results.

4.2.3 Simply supported rectangular laminated sandwich plate under distributed load of sinusoidal variation

A rectangular sandwich plate (0/90/C/0/90) plate (Figure 3) is subjected to a sinusoidal load of intensity q = q0 sin(x/a) sin(y/b). It has a total thickness of h where the thickness of the core is 0. 8h and that of each ply in the top and bottom face sheets is 0. 05h. The plate is analyzed by taking aspect ratio b/a=1. 0 and 2.0 and thickness ratios h/a = 0.1 and 0.2 using different mesh sizes. The material properties used for the core and each laminated face sheets are given in Table 4. The non-dimensional values of central deflection obtained in the present analysis are presented with those obtained from the three-dimensional elasticity solution in Table 5.

As regard to three-dimensional elasticity solution [8], it is relevant to mention that the numerical results based on elasticity solution used here for the comparison are not readily available in the paper of Pagano [3] and hence a software code is used to generate results. The present results are in very good agreement with the elasticity solution [3]. Table 5 also shows that the convergence of the results with mesh refinement is excellent.

4.2.4 Simply supported square sandwich plate under sinusoidal load

A square sandwich (0/C/0) plate, simply supported at all four edges and subjected to distributed load of intensity q = q0 sin(x/a) sin(y/b) is analyzed by taking thickness ratios, h/a = 0. 01, 0.02, 0.05 and 0.25. The thickness of each face is 0. 1h. The material properties are as given in Table 4. The deflection, the normal stresses and the shear stresses at some important points obtained by the present formulation are presented in Table 6 along with the three-dimensional elasticity solution of Pagano [23]. The study has been performed for different mesh divisions to show the convergence characteristics. In general present results are close to the elasticity solution with good convergence.

4.2.5 Sandwich plate having different boundary conditions under distributed load of sinusoidal variation

A sandwich plate (0/90/C/0/90) subjected to a distributed load of intensity q=q0 sin(x/a) sin(y/b) is analyzed in this example. It has a total thickness of h where the thickness of the core is 0.8h and that of each ply in the top and bottom face sheets is 0. 05h. The study has been made for two types of boundary conditions. These are SCSC, i.e. two opposite edges simply supported and other two edges clamped and CCCC, i.e. all edges clamped. The analysis is performed for different thickness ratios (h/a = 0. 01, 0.05, 0.10, 0.20 and 0.50). The material properties used for the core and each laminated face sheets are given in Table 4. The non-dimensional values of central deflection and the stresses calculated at the important points are presented in Table 7.

It may be noted from the results that the percentage increase in the transverse displacement for higher thickness ratios, i.e. from 0.20 to 0.50 is much higher than those for lower thickness ratios, i.e. from 0.01 to 0.05 for both the boundary conditions. The major cause behind it may be that for higher thickness ratios the effect of transverse flexibility of the core and shear deformation effects are more pronounced as compared with that of the lower thickness ratio.

5 CONCLUSIONS

An efficient FE model based on higher order zigzag theory (RHSDT) is presented for the analysis of composites and sandwich plates. The present zigzag theory ensures shear free conditions at the top and bottom of the plate, cubic variation of the in-plane displacements and inter-laminar shear stress continuity at the layer interfaces. An eight node isoparametric element with seven degrees of freedom per node is adopted in this study to model the RHSDT. Convergence of the results is very good indicating reasonably less number of elements required to get the desired results. The results obtained for the deflections and the stresses of a laminated composite and sandwich plate show very good performance of the present formulation and hence it can be recommended for the accurate analysis of laminated composites and sandwich plates.

Received 7 Jan 2011

In revised form 6 Apr 2011

Publication Dates

History