Least Square Finite Element Model for Analysis of Multilayered Composite Plates under Arbitrary Boundary Conditions ()
1. Introduction
The combination of one or more engineering materials on a macroscopic level produces a composite which has better engineering properties than parent materials. Multilayered composite plates are extensively used in aerospace, automotive, shipbuilding industry, and other structural components due to their high strength than the parent material. Their elastic modulus and ultimate strength are tailored to meet various design requirements. They are also applied in sandwich panels, in which two stiff outer panels, known as skin or face sheets, are bound to a softer inner material called core. Laminates are made of numerous laminae also known as ply or layer with each laminae having unidirectional fibers. The volume fiber-to-matrix ratio largely depends on the specific purpose of the design. It has proven very challenging to accurately predict the response characteristics of multilayered structures due to their intrinsic anisotropy, heterogeneity, and the low ratio of transverse shear modulus to the in-plane Young’s modulus [1] .
Literature has reported many higher-order plate theories for the analysis of laminated plates. These studies are categorized into experimental, analytical, or numerical studies of a composite’s response to static, transient loading. Reddy [2] [3] [4] studied the theoretical difference in equivalent single-layer or layerwise variable descriptions as well as in the chosen unknown variables, either displacement, stress, or mixed formulations. In equivalent single-layer variable descriptions, the variables are introduced for the whole plate or shell, whereas in layerwise variable descriptions, each layer is seen as an independent plate or shell, so the number of independent variables is dependent on the number of plies. Pagano [5] [6] demonstrated the three-dimensional elasticity solutions, which showed that multilayered composite structures may exhibit complicating effects introduced by anisotropic behavior, like high transverse deformability, zig-zag effects, and interlaminar continuity. Multilayered composite structures can show high in-plane anisotropy due to different mechanical-physical properties in different in-plane directions [7] . Tauchert [8] provided exact elasticity solutions to the plane-strain deformation of orthotropic simply supported laminates using the method of displacement potentials. Reddy [3] has developed well known simple third-order shear deformation theory for the mechanical analysis of laminated composite plates. A major challenge encountered in solving problems for sandwich and laminate plates is accurately finding transverse stresses that can cause delamination between adjacent plies while satisfying the continuity conditions [9] . The linear elasticity equations have been solved by expressing the three displacements as double Fourier series in the in-plane coordinates for simply supported edges, and static mechanical loading. Ordinary differential equations are deduced for them in the thickness direction [6] [10] . Vel and Batra [1] [11] utilized the Eshelby-Stroh formalism to satisfy boundary conditions at the edges in the sense of Fourier series that rely on St. Venant’s principle, for arbitrary boundary conditions at the edges. Reissner [12] [13] proposed a mixed formulation as a tool to variationally derive governing equilibrium and constitutive equations in terms of independent variables such as displacements and transverse stresses. Equivalent single-layer models (ESLM) made use of this mixed formulation, but least-squares finite element method (LSFEM) is required to accurately describe the local response of multilayered composite structures. Specifically, Carrera [14] [15] [16] [17] developed LSFEMs based on standard finite element method where the Reissner’s functional was minimized to completely and a priori fulfill C0 continuous functions in the thickness z-direction requirements with very successful results. There after the layerwise mixed model by Carrera was extended by Garcia Lage et al. [18] [19] to apply to piezoelectric/magnetic plates to achieve good results as well.
Finite element models are normally based on weak formulations, with displacement-based or mixed formulations. One of the difficult issues encountered in finite element modeling is shear-locking problems, where computational difficulties arise when modeling thin plates or shells. Finite element models, either layerwise or equivalent single layer, displacement-based or mixed, are often based on weak formulations. The alternative is weighted residual formulations, among which the least-squares formulation is quite unique in its basic idea of minimizing the error introduced in the approximation of the governing equations.
Numerical integration techniques take care of shear-locking during computation and sometimes the shear-locking can be handled by using higher-order elements that experience less locking but at the expense of a slower convergence, [7] . An additional issue in mixed weak formulations is that the finite element approximation spaces must satisfy a condition called the inf-sup condition [20] [21] . Weighted residual formulations can serve as an alternative to weak formulations. The basic idea of least-squares formulations is that it minimizes the error introduced in the approximation of the governing equations. In other words, the combination of the least-squares variational principle and a mixed formulation leads to an unconstrained variational minimization problem, where the finite element approximating spaces can be independently chosen. Hence, stability requirements for the inf-sup condition are not needed. The matrices obtained from mixed least-squares formulations are symmetric and positive-definite whereas the ones formed from weak formulations are generally symmetric but not positive definite in general. Initial studies on mixed least-squares formulations show interesting theoretical and computational merits. Pontaza and Reddy [22] were the first to develop a mixed least-squares model for bending of single-layer isotropic laminates, using the classical plate theory (CPT) and first-order shear deformation theory (FSDT). In addition, Pontaza [23] went further to demonstrate the advantage of least-squares finite element models applied to both solid and fluid mechanics. Moleiro et al. [24] [25] and Dervin [9] extended Pontaza and Reddy’s mixed least-squares FSDT finite element model to static and free vibration analysis of multilayered composite and sandwich plates. It is imperative to note that these works indicate that mixed least-squares models are not sensitive to shear-locking when high-order C0 basis or shape functions with full integration in the two-dimensional approximations. Pontaza and Reddy [22] [23] also demonstrated that there is an exponential decay of the least squares functional for increasing the expansion of the polynomial order of the two-dimensional approximations.
The present work provides a least square finite element model also known as layerwise mixed model for the static analysis of multilayered plates using LSFEM under arbitrary boundary conditions, and mechanical loads. The present model is able to fully account for various types of boundary conditions. We utilize a particular type of FEM, called the least-squares finite element method (LSFEM). The benefit of least-squares formulation is that it leads to a variational unconstrained minimization problem, where the finite element approximating spaces can be chosen independently. The purpose of this research work is to determine the efficacy of the proposed methodology in predicting stresses in laminated composite structures when subjected to static mechanical loading under arbitrary boundary conditions.
1.1. Problem Formulation
We consider a rectangular laminate plate of dimension
and use rectangular Cartesian coordinates x1x2 and x3 to describe the position of a material particle in the unstressed reference configuration. Figure 1 describes the static deformation of an N-layer laminate plate occupying the region
. The laminate is composed of transversely isotropic material with material properties varying smoothly in the transverse directions (x2 and x3) [26] . Fiber direction aligned with longitudinal direction is the axis of transverse isotropy. In the absence of conservative body force and internal heat sources, the static analysis of infinitesimal deformation of a linearly elastic lamina of a multilayered plate is governed by the mechanical equilibrium equations shown in Equation (1a) and the geometrical relation between the components of the strain tensor and displacement vector shown in Equation (1b).
Mechanical static equilibrium equations
(1a)
Strain-displacement equations
(1b)
where
are components of the infinitesimal strain tensor,
are the components of the Cauchy stress tensor, indices
, a comma followed by index j denotes partial differentiation with respect to the position
of a material particle, and a repeated index implies summation over the range of the index j, i.e.,
is a contracted form of
, which if expanded further is written in the following form;
;
;
.
![]()
(a) (b)
Figure 1. (a) Geometry and global coordinates of a composite plate; (b) Notation of a multilayered plate.
1.2. Constitutive Equation
The constitutive equations for an orthotropic elastic material in the principal material coordinate is defined as
(2a)
In matrix form, Equation (2a) is represented as
(2b)
Voigt notation is applied as given by Equation (3) where the in-plane stresses are separated from the out-of-plane stresses i.e.
and
are switched, so as
and
. The same applies to strain.
(3a)
(3b)
The constitutive law for the three-dimensional stress of the kth layer of laminae in the principal material coordinate in Voigt matrix notation is given as.
(4)
C is a 6 × 6 symmetric matrix of elasticities of the layer material in the stress-free configuration [27] .
The constitutive relation for a three-dimensional stress state of the kth layer of laminae in the global coordinate is shown below in Voigt notation.
(5a)
where
and
in the global coordinates.
Here, subscript x, y, z denotes the global coordinates, but we will maintain 1, 2, ∙∙∙, 6 in Voigt notation in the global coordinates for convenience.
is a 6 × 6 symmetric matrix of the components of the transformed stiffness matrix and it is related to C by
, and T is a 6 × 6 transformation matrix available in the literature [3] .
In other to calculate the residuals needed to solve for the unknown field variables namely, displacements, transverse stresses, and in-plane strains,
, we compute the mixed form of modulus of elasticity also known as mixed Hook’s law. The elastic constitutive equations as highlighted above are used and a layerwise variable description is used for displacements, transverse stresses, and in-plane strains, taken as independent field variables. Therefore, the in-plane stresses, as well as the transverse strains are written as a function of in-plane strains, and transverse stresses respectively. Thus, there are a total of 9unknown independent field variables for each layer.
(5b)
is computed by expressing the in-plane stresses and transverse strains in terms of in-plane strains and transverse stresses (layerwise continuous variables and out-of-plane quantities) [27] [28] given by
for
for
for
(6)
for
for
We use a state-space formulation and take
as unknown variable at a point in each plate layer. The residual on the kth layer involves the first-order derivative for the element
.
Strain-displacement equations model
for
(7a)
for
(7b)
for
(7c)
Elastic constitutive equations model
(8a)
(8b)
(8c)
Equilibrium equations model
From
;
(9a)
;
(9b)
;
(9c)
Altogether, there are nine residuals,
, where
which leads to 9 equations with 9 unknowns. The displacements boundary conditions are therefore directly impose because they are independent variables of the present formulation. Hence we write the residual boundary conditions in terms of the variable
since the in-plane stresses are not computed directly in LSFEM. The least-square finite element model allows the introduction of additional residuals in the functional in the least square sense, corresponding to boundary conditions imposed in a least-squares sense. In this case, the constitutive relation described in Equation (2) is used to write the in-plane stresses in terms of the model-independent field variables, so that additional residuals can be included in the least-squares functional.
(10)
Three residual boundary conditions are obtained for each bounding face of a cubic composite, resulting in 18 residuals denoted by
through
.
1.3. Boundary Conditions
The mechanical static analysis focuses on multilayered plates under transverse mechanical loads applied on the plate top surface and/or bottom surface. Therefore, a transverse mechanical load is considered by transverse normal stress,
, which is applied on the plate’s top and/or bottom surfaces. The mechanical boundary conditions prescribed on the top and/or bottom surfaces are surface traction vector
. Nonzero normal and zero tangential tractions are typically prescribed on these surfaces. The traction component prescribed is
sinusoidal surface traction components.
are prescribed on
,
Where
is the load intensity, m, and n are any possible integers (Table 1).
To satisfy the continuity condition, we ensure that transverse normal stresses, in-plain strains, and displacements through the layer interface are equal. In other words, layers interfaces are perfectly bonded together. The values of
,
and
are equal through the interface.
![]()
Table 1. Boundary and loading condition on two surfaces.
(11)
and
,
,
Due to interfacial continuity on x3, it shows that the jump in values of transverse stresses, in-plane strains, and displacement across the interface is zero. “t” denotes the top and “b” the bottom surfaces of a layer, respectively. As pointed earlier, the boundary conditions are simply supported on
and b in all cases. i.e., the laminate and sandwich plates are simply supported on
and b;
,
,
in all cases. The boundary conditions at the edges on
and a can be clamped, simply supported, free edge, slippery surface, and other four arbitrary boundary conditions [1] , that are summarized in Table 2.
1.4. Least Square Finite Element Model
The Least square finite model is developed to allow the use of basis functions C0 in the two-dimensional approximations, to reduce the higher regularity requirements common to all weighted residual formulations. Thus, the system of governing equations is preferably transformed into an equivalent first-order system, which then requires additional independent variables. The least-squares functional is developed by measuring the norms of the residuals of all the governing equations. Generally, the least-squares functional is derived by the sum of the squared residuals of each governing equation for the kth layer, considering of each lamina.
The functional J is defined according to [9] in terms of the residuals with contributions from the material layer as
(12)
where
is residual boundary conditions, and the repeated index a goes from 1 to 9, Ω denotes the in-plane domain of the plate. The second surface integral is
![]()
Table 2. Boundary conditions at the edges.
So B1B1 denotes clamped edge an
and
.
for the bounding surfaces of a composite and
are residuals for the boundary conditions for the kth layer in the bounding surface and b goes from 10 to 27.
is expressed in terms of the basis functions that are a product of Co continuous Lagrange polynomials of degree N1, N2, N3 in x1, x2, x3 directions.
(13)
Ni is the number of nodes in each axis, altogether there are
unknowns for each layer.
The basis function is defined as
(14)
the basis function
is expressed in terms of its natural coordinates and P-th order,
is the Lagrange polynomial and
, its derivative,
is the root of the Legendre polynomial
of n-order. Substituting Equation (14) into Equation (13) leads to
(15)
Substituting the resulting Equation (15) into Equation (12), then the integral is evaluated numerically by using the Gauss-Lobatto quadrature rules in each of x1x2 and x3. Strain energy J is minimized to obtain an algebraic equation by taking the derivative of J with respect to
.
(16)
The layerwise mixed least-squares model for static analysis of multilayered composite plates ultimately yields the following symmetric positive definite system of linear equations KA = F
(17)
2. Results
We study the problem analytically analyzed by Vel et al. [1] and Moleiro [29] to verify our algorithm and to establish the accuracy of computed results for composite and sandwich plate respectively. Each lamina consists of a unidirectional fiber-reinforced material that is modeled as orthotropic and is of equal thickness layers, square laminate of side length a, with the sinusoidal loads applied only on the top surface. In the sandwich plate the face sheet is modeled as orthotropic whereas the core is modeled as isotropic. We express the results of the composite plate in the following nondimensionalized form according to Pagano [5] .
,
(18a)
,
Similarly, we express the results of the sandwich plate in the following nondimensionalized form according to Moleiro [29]
,
(18b)
,
where
,
.
The two load distributions given in Table 1 are considered. The following three cases of lamination are considered and compared with the analytical solutions obtained in Vel and Batra [1] :
Case 1; A two-ply [0˚/90˚] laminate with the fibers parallel to the x1 and the x2 directions in the bottom and the top layers, respectively.
Case 2; A three-ply laminate [0˚/90˚/0˚] with the fibers parallel to the x1, x2, and x1 directions in the bottom, middle, and top layers, respectively.
Case 3; A three-layer sandwich plate made of (isotropic) polyvinyl chloride (PVC) foam core with thickness 2h/3, along with (orthotropic) carbon fiber reinforced polymer (CFRP)composite skins of 0° and thickness h/6 each skin for load case a.
Our results of the selected quantities are verified with those in Vel [1] for plate aspect ratio a/h = 5, and 10 in the following tables below.
The mechanical properties of the test material [6] , the (orthotropic) carbon fiber reinforced polymer CFRP and isotropic polyvinyl chloride (PVC) foam PVC [29] [30] [31] are given in Table 3 below.
Case 1; A two-ply [0˚/90˚] composite laminate
We present the nondimensionalized results for this Case 1 composite laminate [0˚/90˚]. The result of nondimensionalized quantities for this Case 1 composite laminate is made of test material composite layers, under mechanical load are shown in detail in the following Tables 4-6, each for a different aspect (side to thickness) ratio, i.e., a/h = 5, 10 respectively. We also provide a clearer description of this Case 1 composite laminate behaviour under mechanical loading, Figures 2-5 describe through-the-thickness distributions of a number of mechanical quantities, each for a different aspect (side to thickness) ratio, i.e., a/h = 5, 10 respectively.
Case 2; A three-ply laminate [0˚/90˚/0˚]
We present the nondimensionalized results for Case 2 composite laminate [0˚/90˚/0˚] made of test material composite layers. The results of nondimensionalized quantities are shown in detail in the following Tables 7-9, each for a different aspect (side to thickness) ratio, i.e., a/h = 5, and 10 respectively. Figures 6-9 present the through-the-thickness distributions of a number of mechanical quantities, each for a different aspect ratio, i.e., a/h = 5, and 10 respectively.
Case 3; A three-layer sandwich plate [0˚/core/0˚]
We consider a sandwich plate made of PVC foam core and CFRP composite skin, under mechanical load and arbitrary boundary conditions. The same intensity of the mechanical load is also considered for the sandwich plate as required for the nondimensionalized form in Equation (18b). the numerical results are displayed in Table 10 for a different aspect (side to thickness) ratio, i.e., a/h = 4, 10, and 10. Figures 10-13 present the through-the-thickness distributions of a number of mechanical quantities, each for a different aspect ratio, i.e., a/h = 4, 10, and 10 respectively.
3. Discussion
The condition number of the matrix K is improved by using non-dimensional variables throughout the formulations and subsequently improve and reduce the error. In LSFEM, the basis functions are not required to satisfy the essential boundary conditions unlike in the Ritz method where the essential boundary conditions must be satisfied. The Ritz method allows one to apply penalty method in order to enforce the essential boundary condition whereas in finite element method, each ply or laminae is discretized into disjointed domains. The accuracy of the numerical solution can be improved by either increasing the order of polynomials in the basis functions or reducing the size of element or both.
![]()
Table 3. Mechanical properties of the materials used.
where E, G (in Gpa), and v are Young’s modulus, shear modulus, and Poisson’s ratio respectively.
![]()
Table 4. Comparison of present solution with the 3D exact solution of Vel [1] for [0˚/90˚] laminate subjected mechanical load case a, and a/h = 5; the superscript of stress and displacement quantities indicate the in-plane and through-thickness location: 1 = (a/2, a/2, −h/2), 2 = (a/2, a/2, h/2), 3 = (a/8, 0, −h/2), 4 = (a/8, a/2, 0), 5 = (a/2, 0, 0), 6 = (a/2, a/2, 0), 7 = (a/4, a/2, h/2), 8 = (a/2, a/4, −h/2). This applies to the quantities in Tables 5-7.
![]()
Table 5. Comparison of present solutions with the 3D exact solution of Vel [1] for [0˚/90˚] laminate subjected mechanical load case b, and a/h = 5.
![]()
Table 6. Comparison of present solutions with the 3D exact solution of Vel [1] for [0˚/90˚] laminate subjected under mechanical load case a, and a/h = 10.
![]()
(a) (b)
Figure 2. Comparison of present results and Vel [1] for
, and
, both on B4, load case a, and a/h = 5, 10, case 1.
![]()
(a) (b)
Figure 3. Comparison of present results and Vel [1] for
, on B3 and
, on B6, load case a, and a/h = 5, 10, case 1.
![]()
(a) (b)
Figure 4. Comparison of present results and Vel [1] for
, on B8, and
, on B7, load case a, and a/h = 5, 10, case 1.
![]()
(a) (b)
Figure 5. Comparison of present results and Vel [1] for
, on B1 and
, on B2, load case a, and a/h = 5, 10, case 1.
![]()
Table 7. Comparison of present solutions with 3D exact solution of Vel [1] for [0˚/90˚/0˚] laminate subjected mechanical load case a; the superscript of stress and displacement quantities indicate the in-plane and through-thickness location: 1 = (a/2, a/2, −h/2), 2 = (a/2, a/2, h/2), 3 = (a/2, a/2, −2h/3), 4 = (a/2, a/2, 2h/3), 5 = (a/2, a/2, 0), 6 = (a/2, 0, 0). This applies to the quantities in Tables 7-9.
![]()
Table 8. Comparison of present solutions with 3D exact solution of Vel [1] for [0˚/90˚/0˚] laminate subjected mechanical load case a.
![]()
Table 9. Comparison of present solutions with 3D exact solution of Vel [1] for [0˚/90˚/0˚] laminate subjected mechanical load case a.
![]()
(a) (b)
Figure 6. Comparison of present results and Vel [1] for
,
, both on B4 load case a, and a/h = 5, 10, case 2.
In the LSFEM model, case 1, for instance, we used mesh size = 1, the number of degrees of freedom DOF, generated (for N1 = N2 = 6, and N3 = 4) = 5733. The continuity conditions hold for the fact that the values of
,
and
at the top of the layer are the same as
,
and
at the bottom of layer respectively. After refining the mesh size by increasing the number of polynomials, the corresponding number of degrees of freedom (for N1 = N2 = 8, and
![]()
(a) (b)
Figure 7. Comparison of present results and Vel [1] for
, on B3,
, on B6 load case a, and a/h = 5, 10, case 2.
![]()
(a) (b)
Figure 8. Comparison of present results and Vel [1] for
, on B8 and
, on B7 load case a, and a/h = 5, 10, case 2.
![]()
(a) (b)
Figure 9. Comparison of present results and Vel [1] for
, on B1 and
, on B2, load case a, and a/h = 5, 10, case 2.
![]()
Table 10. Comparison of present solutions of [0˚/core/0˚] sandwich laminate subjected to mechanical load case (a) and arbitrary boundary conditions and compared with 3D exact solution of Moleiro [29] ; the superscript of stress and displacement quantities indicate the in-plane and through-thickness location: 1 = (a/2, a/2, −h/2), 2 = (a/2, a/2, h/2), 3 = (0, 0, h/2), 4 = (0, a/2, 0), 5 = (a/2, 0, 0), 6 = (a/2, a/2, 0), 7 = (0, a/2, −h/2).
![]()
(a) (b)
Figure 10. Comparison of present results and Moleiro [29] for
,
, both on B5 load case a, and a/h = 4, 10, 50, case 3.
![]()
(a) (b)
Figure 11. Comparison of present results and Moleiro [29] for
,
, both on B5, load case a, and a/h = 4, 10, 5, case 3.
![]()
(a) (b)
Figure 12. Comparison of present results and Moleiro [29] for
,
, both on B5, load case a, and a/h = 4, 10, 50, case 3.
![]()
(a) (b)
Figure 13. Comparison of present results and Moleiro [29] for
,
, both on B5, load case a, and a/h = 4, 10, 50, case 3.
N3 = 4) = 9477. The error for in-plane-stress error reduced from 1.19% to 0.15% and out-of-plane normal stress reduced marginally from 1% to 0.3% but the computation error for the shear stress
reduced from 19.19% to 15.2%. The-through-the-thickness variations remain the same for the two DOFs. To verify the efficiency of the present model and its predictive capabilities. Figures 2-7 show that the predicted distributions through the plate thickness for the variables
,
,
,
,
,
,
for a/h = 5, 10, alongside the exact three-dimensional solutions. The plots and tables provide detailed descriptions of the composite laminate [0˚/90˚] and [0˚/90˚/0˚] static analysis under arbitrary boundary condition and corresponding exact solutions. In all the figures, both the exact three-dimensional solutions and the present model solutions are represented using the corresponding symbolic toolbox.
4. Conclusion
We present a three-dimensional linear elasticity numerical solution of a least square finite element model formulation for static analysis of multilayered composite and sandwich plates under arbitrary boundary conditions and compare our results with exact solution available in literature. The model is developed in a least squares sense by taking the three transverse stresses, three strain-displacement relations in the xy-plane and three displacement components as independent variables and using the least-squares method to minimize the residuals in the expressions for these nine field variables. The numerical examples show the accuracy and capabilities of our least square finite element model. The present model results judgement is based on a comprehensive comparison with the exact three-dimensional elasticity solutions by Vel and Batra [1] and Moleiro for sandwich plate [29] . The numerical examples focus on the static analysis of the simply supported square multilayered composite and sandwich plate on y = 0, b, and arbitrary BC support on x = 0, a, under a sinusoidal transverse load either at top or bottom surface of the laminate. The ranges of aspect ratios a/h considered are 5, 10 for composite plate and 4, 10, 50 for sandwich plate, i.e., from thin plate to thick plate. Even though we show that the present least-squares finite element models give very accurate results and in excellent agreement with the exact three-dimensional solutions, the accuracy of this model depends on the thickness of plate. For a thin plate of a/h = 10 upwards, the model produces highly accurate results, that is, as the thickness of the plate decreases, the accuracy of the model increases, and computation error diminishes. To increase the accuracy of the model for a very thick plate, such as a/h = 1, 2, 4, we increase the z-polynomial order and refine the mesh, however this increases the computation cost and time. In the future we will introduce a combination of thermal and mechanical load. In fact, in addition to mechanical loading of composite and sandwich structures, thermal effects can be fully coupled into the model to analyze the static or dynamic thermo-mechanical deformations. This can be accomplished by adding the energy equation and an appropriate constitutive law for the heat fluxes into the mechanical and thermal equilibrium governing equations and incorporating the temperature field and transverse heat fluxes into the list of field unknown variables to solve for.