Least Square Finite Element Model for Analysis of Multilayered Composite Plates under Arbitrary Boundary Conditions

Christian Mathew^1*, Yao Fu²

¹Department of Engineering Mechanics, Virginia Polytechnic Institute and State University, Blacksburg, Virginia, USA.
²Department of Aerospace and Ocean Engineering, Virginia Polytechnic Institute and State University, Blacksburg, Virginia, USA.

DOI: 10.4236/wjet.2024.121003                

Abstract

Laminated composites are widely used in many engineering industries such as aircraft, spacecraft, boat hulls, racing car bodies, and storage tanks. We analyze the 3D deformations of a multilayered, linear elastic, anisotropic rectangular plate subjected to arbitrary boundary conditions on one edge and simply supported on other edge. The rectangular laminate consists of anisotropic and homogeneous laminae of arbitrary thicknesses. This study presents the elastic analysis of laminated composite plates subjected to sinusoidal mechanical loading under arbitrary boundary conditions. Least square finite element solutions for displacements and stresses are investigated using a mathematical model, called a state-space model, which allows us to simultaneously solve for these field variables in the composite structure’s domain and ensure that continuity conditions are satisfied at layer interfaces. The governing equations are derived from this model using a numerical technique called the least-squares finite element method (LSFEM). These LSFEMs seek to minimize the squares of the governing equations and the associated side conditions residuals over the computational domain. The model is comprised of layerwise variables such as displacements, out-of-plane stresses, and in- plane strains, treated as independent variables. Numerical results are presented to demonstrate the response of the laminated composite plates under various arbitrary boundary conditions using LSFEM and compared with the 3D elasticity solution available in the literature.

Keywords

Multilayered Composite and Sandwich Plate, Transverse Stress, Continuity Condition, Arbitrary Boundary Condition, Layerwise Theory, Least-Squares Formulation

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 C⁰ 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 C⁰ 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 $a\times b\times h$ and use rectangular Cartesian coordinates x₁x₂ and x₃ 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 $\left(0,a\right)\times \left(0,b\right)\times \left(-h/2,h/2\right)$ . The laminate is composed of transversely isotropic material with material properties varying smoothly in the transverse directions (x₂ and x₃) [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

${\sigma }_{ij,j}=0$ (1a)

Strain-displacement equations

${ϵ}_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right)$ (1b)

where ${ϵ}_{ij}$ are components of the infinitesimal strain tensor, ${\sigma }_{ij}$ are the components of the Cauchy stress tensor, indices $i,j=1,2,3$ , a comma followed by index j denotes partial differentiation with respect to the position $x_j$ of a material particle, and a repeated index implies summation over the range of the index j, i.e., ${\sigma }_{ij,j}$ is a contracted form of $\frac{\partial {\sigma }_{ij}}{\partial x_j}$ , which if expanded further is written in the following form; $\frac{\partial {\sigma }_{11}}{\partial x_1}+\frac{\partial {\sigma }_{12}}{\partial x_2}+\frac{\partial {\sigma }_{13}}{\partial x_3}=0$ ; $\frac{\partial {\sigma }_{21}}{\partial x_1}+\frac{\partial {\sigma }_{22}}{\partial x_2}+\frac{\partial {\sigma }_{23}}{\partial x_3}=0$ ; $\frac{\partial {\sigma }_{31}}{\partial x_1}+\frac{\partial {\sigma }_{32}}{\partial x_2}+\frac{\partial {\sigma }_{33}}{\partial x_3}=0$ .

(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

${\sigma }_{ij}=C_{ijkl}{ϵ}_{ij}$ (2a)

In matrix form, Equation (2a) is represented as

$\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{23}\\ {\sigma }_{13}\\ {\sigma }_{12}\end{array}\right]=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{12}& C_{22}& C_{23}& 0& 0& 0\\ C_{13}& C_{23}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right]\left[\begin{array}{c}{ϵ}_{11}\\ {ϵ}_{22}\\ {ϵ}_{33}\\ {ϵ}_{23}\\ {ϵ}_{13}\\ {ϵ}_{12}\end{array}\right]$ (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. ${\sigma }_{12}$ and ${\sigma }_{33}$ are switched, so as ${\sigma }_{23}$ and ${\sigma }_{13}$ . The same applies to strain.

${\left[{\sigma }_1,{\sigma }_2,{\sigma }_3,{\sigma }_4,{\sigma }_5,{\sigma }_6\right]}^{\text{T}}={\left[{\sigma }_{11},{\sigma }_{22},{\sigma }_{12},{\sigma }_{13},{\sigma }_{23},{\sigma }_{33}\right]}^{\text{T}}$ (3a)

${\left[{ϵ}_1,{ϵ}_2,{ϵ}_3,{ϵ}_4,{ϵ}_5,{ϵ}_6\right]}^{\text{T}}={\left[{ϵ}_{11},{ϵ}_{22},2{ϵ}_{12},2{ϵ}_{13},2{ϵ}_{23},{ϵ}_{33}\right]}^{\text{T}}$ (3b)

The constitutive law for the three-dimensional stress of the k^th layer of laminae in the principal material coordinate in Voigt matrix notation is given as.

$\left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right]=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& 0& 0& 0& C_{13}\\ C_{12}& C_{22}& 0& 0& 0& C_{23}\\ 0& 0& C_{66}& 0& 0& 0\\ 0& 0& 0& C_{55}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ C_{13}& C_{23}& 0& 0& 0& C_{33}\end{array}\right]\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ {ϵ}_3\\ {ϵ}_4\\ {ϵ}_5\\ {ϵ}_6\end{array}\right]$ (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 k^th layer of laminae in the global coordinate is shown below in Voigt notation.

${\left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right]}^g=\left[\begin{array}{cccccc}{\bar{C}}_{11}& {\bar{C}}_{12}& {\bar{C}}_{16}& 0& 0& {\bar{C}}_{13}\\ {\bar{C}}_{12}& {\bar{C}}_{22}& {\bar{C}}_{26}& 0& 0& {\bar{C}}_{23}\\ {\bar{C}}_{16}& {\bar{C}}_{26}& {\bar{C}}_{66}& 0& 0& {\bar{C}}_{36}\\ 0& 0& 0& {\bar{C}}_{55}& {\bar{C}}_{45}& 0\\ 0& 0& 0& {\bar{C}}_{45}& {\bar{C}}_{44}& 0\\ {\bar{C}}_{13}& {\bar{C}}_{23}& {\bar{C}}_{36}& 0& 0& {\bar{C}}_{33}\end{array}\right]{\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ {ϵ}_3\\ {ϵ}_4\\ {ϵ}_5\\ {ϵ}_6\end{array}\right]}^g$ (5a)

where ${\left({\left[{\sigma }_1,{\sigma }_2,{\sigma }_3,{\sigma }_4,{\sigma }_5,{\sigma }_6\right]}^g\right)}^{\text{T}}={\left[{\sigma }_{xx},{\sigma }_{yy},{\sigma }_{xy},{\sigma }_{xz},{\sigma }_{yz},{\sigma }_{zz}\right]}^{\text{T}}$ and ${\left({\left[{ϵ}_1,{ϵ}_2,{ϵ}_3,{ϵ}_4,{ϵ}_5,{ϵ}_6\right]}^g\right)}^{\text{T}}={\left[{ϵ}_{xx},{ϵ}_{yy},2{ϵ}_{xy},2{ϵ}_{xz},2{ϵ}_{yz},{ϵ}_{zz}\right]}^{\text{T}}$ 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. $\bar{C}$ is a 6 × 6 symmetric matrix of the components of the transformed stiffness matrix and it is related to C by $\bar{C}=TC{T}^{\text{T}}$ , 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, $S^k=\left(u_1^k,u_2^k,u_3^k,{\sigma }_4^k,{\sigma }_5^k,{\sigma }_6^k,{ϵ}_1^k,{ϵ}_2^k,{ϵ}_3^k\right)$ , 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.

${\left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {ϵ}_4\\ {ϵ}_5\\ {ϵ}_6\end{array}\right]}^g=\left[\begin{array}{cccccc}{\hat{C}}_{11}& {\hat{C}}_{12}& {\hat{C}}_{16}& 0& 0& {\hat{C}}_{13}\\ {\hat{C}}_{12}& {\hat{C}}_{22}& {\hat{C}}_{26}& 0& 0& {\hat{C}}_{23}\\ {\hat{C}}_{16}& {\hat{C}}_{26}& {\hat{C}}_{66}& 0& 0& {\hat{C}}_{36}\\ 0& 0& 0& {\hat{C}}_{55}& {\hat{C}}_{45}& 0\\ 0& 0& 0& {\hat{C}}_{45}& {\hat{C}}_{44}& 0\\ -{\hat{C}}_{13}& -{\hat{C}}_{23}& -{\hat{C}}_{36}& 0& 0& {\hat{C}}_{33}\end{array}\right]\left[\begin{array}{c}{ϵ}_1\\ {ϵ}_2\\ {ϵ}_3\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right]$ (5b)

${\hat{C}}_{ij}$ 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

${\hat{C}}_{ij}={\bar{C}}_{ij}-\frac{{\bar{C}}_{i3}}{{\bar{C}}_{33}}{\bar{C}}_{j3}$ for $i,j=1,2,6$

${\hat{C}}_{i3}=\frac{{\bar{C}}_{i3}}{{\bar{C}}_{33}}$ for $i,j=1,2,6$

${\hat{C}}_{i3}=-\frac{{\bar{C}}_{i3}}{{\bar{C}}_{33}}$ for $i,j=1,2,6$ (6)

${\hat{C}}_{ij}=-\frac{{\bar{C}}_{ij}}{{\bar{C}}_{55}{\bar{C}}_{44}-{\bar{C}}_{45}{\bar{C}}_{45}}$ for $i,j=4,5$

${\hat{C}}_{ii}=\frac{C_{jj}}{{\bar{C}}_{55}{\bar{C}}_{44}-{\bar{C}}_{45}{\bar{C}}_{45}}$ for $i,j=4,5$

${\hat{C}}_{33}=\frac{1}{{\bar{C}}_{33}}$

We use a state-space formulation and take $S^k$ as unknown variable at a point in each plate layer. The residual on the k^th layer involves the first-order derivative for the element $S^k$ .

Strain-displacement equations model

for $i,j=1$

${ϵ}_{11}=\frac{1}{2}\left(u_{1,1}+u_{1,1}\right)=u_{1,1}$

$R_1^k=u_{1,1}^k-{ϵ}_1^k$ (7a)

for $i,j=2$

${ϵ}_{22}=\frac{1}{2}\left(u_{2,2}+u_{2,2}\right)=u_{2,2}$

$R_2^k=u_{2,2}^k-{ϵ}_2^k$ (7b)

for $i=1,j=2$

${ϵ}_{12}=\frac{1}{2}\left(u_{1,2}+u_{2,1}\right)$

$2{ϵ}_{12}=u_{1,2}+u_{2,1}$

$R_3^k=u_{1,2}^k+u_{2,1}^k-{ϵ}_3^k$ (7c)

Elastic constitutive equations model

${\sigma }_4={\bar{C}}_{55}{ϵ}_4+{\bar{C}}_{45}{ϵ}_5=2{\bar{C}}_{55}{ϵ}_{13}+2{\bar{C}}_{45}{ϵ}_{23}$

${\sigma }_4={\bar{C}}_{55}{u}_{1,3}+{\bar{C}}_{55}{u}_{3,1}+{\bar{C}}_{45}{u}_{2,3}+{\bar{C}}_{45}{u}_{3,2}$

$R_4^k={\sigma }_4^k-{\bar{C}}_{55}^k{u}_{1,3}^k-{\bar{C}}_{55}^k{u}_{3,1}^k-{\bar{C}}_{45}^k{u}_{2,3}^k-{\bar{C}}_{45}^k{u}_{3,2}^k$ (8a)

${\sigma }_5={\bar{C}}_{45}{ϵ}_4+{\bar{C}}_{44}{ϵ}_5=2{\bar{C}}_{45}{ϵ}_{13}+2{\bar{C}}_{44}{ϵ}_{23}$

${\sigma }_5={\bar{C}}_{45}{u}_{1,3}+{\bar{C}}_{45}{u}_{3,1}+{\bar{C}}_{44}{u}_{2,3}+{\bar{C}}_{44}{u}_{3,2}$

$R_5^k={\sigma }_5^k-{\bar{C}}_{45}^k{u}_{1,3}^k-{\bar{C}}_{45}^k{u}_{3,1}^k-{\bar{C}}_{44}^k{u}_{2,3}^k-{\bar{C}}_{44}^k{u}_{3,2}^k$ (8b)

${\sigma }_6={\bar{C}}_{13}{ϵ}_1+{\bar{C}}_{23}{ϵ}_2+{\bar{C}}_{36}{ϵ}_3+{\bar{C}}_{33}{ϵ}_6$

${\sigma }_6={\bar{C}}_{13}{u}_{1,1}+{\bar{C}}_{23}{u}_{2,2}+{\bar{C}}_{36}\left(u_{1,2}+u_{2,1}\right)+{\bar{C}}_{33}{u}_{3,3}$

$R_6^k={\sigma }_6^k-{\bar{C}}_{13}^k{u}_{1,1}^k+{\bar{C}}_{23}^k{u}_{2,2}^k+{\bar{C}}_{36}^k{u}_{1,2}^k+{\bar{C}}_{36}^k{u}_{2,1}^k+{\bar{C}}_{33}{\bar{C}}_{36}^k{u}_{2,3}^k$ (8c)

Equilibrium equations model

${\sigma }_1={\hat{C}}_{11}{ϵ}_1+{\hat{C}}_{12}{ϵ}_2+{\hat{C}}_{16}{ϵ}_3+{\hat{C}}_{13}{\sigma }_6$

${\sigma }_2={\hat{C}}_{12}{ϵ}_1+{\hat{C}}_{22}{ϵ}_2+{\hat{C}}_{26}{ϵ}_3+{\hat{C}}_{23}{\sigma }_6$

${\sigma }_3={\hat{C}}_{16}{ϵ}_1+{\hat{C}}_{26}{ϵ}_2+{\hat{C}}_{66}{ϵ}_3+{\hat{C}}_{36}{\sigma }_6$

From ${\sigma }_{ij,j}=0$

${\sigma }_{11,1}+{\sigma }_{12,2}+{\sigma }_{13,3}=0$ ; ${\sigma }_{1,1}+{\sigma }_{3,2}+{\sigma }_{4,3}=0$

${\hat{C}}_{11}{ϵ}_{1,1}+{\hat{C}}_{12}{ϵ}_{2,1}+{\hat{C}}_{16}{ϵ}_{3,1}+{\hat{C}}_{13}{\sigma }_{6,1}+{\hat{C}}_{16}{ϵ}_{1,2}+{\hat{C}}_{26}{ϵ}_{2,2}+{\hat{C}}_{66}{ϵ}_{3,2}+{\bar{C}}_{36}{\sigma }_{6,2}+{\sigma }_{4,3}=0$

$R_7^k={\hat{C}}_{11}^k{ϵ}_{1,1}^k+{\hat{C}}_{12}^k{ϵ}_{2,1}^k+{\hat{C}}_{16}^k{ϵ}_{3,1}^k+{\hat{C}}_{13}^k{\sigma }_{6,1}+{\hat{C}}_{16}^k{ϵ}_{1,2}^k+{\hat{C}}_{26}^k{ϵ}_{2,2}^k+{\hat{C}}_{66}^k{ϵ}_{3,2}^k+{\hat{C}}_{36}^k{\sigma }_{6,2}+{\sigma }_{4,3}$ (9a)

${\sigma }_{21,1}+{\sigma }_{22,2}+{\sigma }_{23,3}=0$ ; ${\sigma }_{3,1}+{\sigma }_{2,2}+{\sigma }_{5,3}=0$

${\hat{C}}_{16}{ϵ}_{1,1}+{\hat{C}}_{26}{ϵ}_{2,1}+{\hat{C}}_{66}{ϵ}_{3,1}+{\hat{C}}_{36}{\sigma }_{6,1}+{\hat{C}}_{12}{ϵ}_{1,2}+{\hat{C}}_{22}{ϵ}_{2,2}+{\hat{C}}_{26}{ϵ}_{3,2}+{\hat{C}}_{23}{\sigma }_{6,2}+{\sigma }_{5,3}=0$

$R_8^k={\hat{C}}_{16}^k{ϵ}_{1,1}^k+{\hat{C}}_{26}^k{ϵ}_{2,1}^k+{\hat{C}}_{66}^k{ϵ}_{3,1}^k+{\hat{C}}_{36}^k{\sigma }_{6,1}+{\hat{C}}_{12}^k{ϵ}_{1,2}^k+{\hat{C}}_{22}^k{ϵ}_{2,2}^k+{\hat{C}}_{26}^k{ϵ}_{3,2}^k+{\hat{C}}_{23}^k{\sigma }_{6,2}+{\sigma }_{5,3}$ (9b)

${\sigma }_{31,1}+{\sigma }_{32,2}+{\sigma }_{33,3}=0$ ; ${\sigma }_{4,1}+{\sigma }_{5,2}+{\sigma }_{6,3}=0$

$R_9^k={\sigma }_{4,1}^k+{\sigma }_{5,2}^k+{\sigma }_{6,3}^k$ (9c)

Altogether, there are nine residuals, $R_{\alpha }^k=0$ , where $\alpha =1,2,3,\cdots ,9$ 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 $S^k$ 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.

$R^k={\hat{C}}_{11}\left({ϵ}_1\right)+{\hat{C}}_{12}\left({ϵ}_2\right)+{\hat{C}}_{16}\left({ϵ}_3\right)+{\hat{C}}_{13}{\sigma }_6$ (10)

Three residual boundary conditions are obtained for each bounding face of a cubic composite, resulting in 18 residuals denoted by $R_{10}^k$ through $R_{27}^k$ .

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, $\sigma \left(x_1,x_2,\pm \frac{h}{2}\right)$ , 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 ${\sigma }_{ij}$ . Nonzero normal and zero tangential tractions are typically prescribed on these surfaces. The traction component prescribed is ${\sigma }_{i3}$ sinusoidal surface traction components. ${\sigma }_{i3}\left(x_1,x_2,\pm \frac{h}{2}\right)=\left[0,0,\pm {\sigma }_0\sin \left(\frac{m\pi x_1}{a}\right)\sin \left(\frac{n\pi x_2}{b}\right)\right]$ are prescribed on $\pm \frac{h}{2}$ ,

Where ${\sigma }_0$ 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 $u_i$ , ${ϵ}_{i3}$ and ${\sigma }_{i3}$ are equal through the interface.

Table 1. Boundary and loading condition on two surfaces.

${\sigma }_{i3t}^k={\sigma }_{i3b}^{k+1}$

$u_{it}^k=u_{ib}^{k+1}$ (11)

$k=1,2$ and $i=1,2,3$

${ϵ}_{11t}^k={ϵ}_{11b}^{k+1}$ , ${ϵ}_{22t}^k={ϵ}_{22b}^{k+1}$ , ${ϵ}_{12t}^k={ϵ}_{12b}^{k+1}$

Due to interfacial continuity on x₃, 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 $y=0$ and b in all cases. i.e., the laminate and sandwich plates are simply supported on $y=0$ and b; ${\sigma }_{22}=0$ , $u_1=0$ , $u_3=0$ in all cases. The boundary conditions at the edges on $x=0$ 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 C₀ 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 k^th 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

$J=\frac{1}{2}\left({\displaystyle \underset{k=1}{\overset{p}{\sum }}{\displaystyle {\int }_{\Omega }{\displaystyle {\int }_{h^k}{\left(R_a^k\right)}^2{d}_{x_3}{d}_{x_1}{d}_{x_2}}}}+{\displaystyle \underset{k=1}{\overset{p}{\sum }}{\displaystyle {\int }_{\Gamma }{\displaystyle {\int }_{h^k}{\left(R_b^k\right)}^2{d}_{x_3}{d}_{x_1}{d}_{x_2}}}}\right)$ (12)

where $R_a^k$ 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 B₁B₁ denotes clamped edge an ${\sigma }_{13}$ and ${\sigma }_{23}$ .

for the bounding surfaces of a composite and $R_b^k$ are residuals for the boundary conditions for the k^th layer in the bounding surface and b goes from 10 to 27.

$s^k$ is expressed in terms of the basis functions that are a product of C^o continuous Lagrange polynomials of degree N₁, N₂, N₃ in x₁, x₂, x₃ directions.

$s_a^k={\displaystyle \underset{i=1}{\overset{N_1}{\sum }}{\displaystyle \underset{j=1}{\overset{N_2}{\sum }}{\displaystyle \underset{m=1}{\overset{N_3}{\sum }}{S}_{a,ijm}^k{\psi }_i\left(x_1\right){\psi }_j\left(x_2\right){\psi }_m\left(x_3\right)}}}$ (13)

N_i is the number of nodes in each axis, altogether there are $9\times N_1\times N_2\times N_3$ unknowns for each layer.

The basis function is defined as

${\psi }_i\left(\xi \right)=\frac{\left(\xi -1\right)\left(1+\xi \right){L^{\prime }}_{p-1}\left(\xi \right)}{p\left(p-1\right)L_{p-1}\left({\xi }_i\right)\left(\xi -{\xi }_i\right)}$ (14)

the basis function ${\psi }_i$ is expressed in terms of its natural coordinates and P-th order, $L_{p-1}\left(\xi \right)$ is the Lagrange polynomial and ${L^{\prime }}_{p-1}\left(\xi \right)$ , its derivative, ${\xi }_i$ is the root of the Legendre polynomial $P_n\left(\xi \right)=0$ of n-order. Substituting Equation (14) into Equation (13) leads to

$s_a^k={\displaystyle \underset{i=1}{\overset{N_1}{\sum }}{\displaystyle \underset{j=1}{\overset{N_2}{\sum }}{\displaystyle \underset{m=1}{\overset{N_3}{\sum }}{S}_{a,ijm}^k\frac{\left(\xi -1\right)\left(1+\xi \right){L^{\prime }}_{p-1}\left(\xi \right)}{p\left(p-1\right)L_{p-1}\left({\xi }_i\right)\left(\xi -{\xi }_i\right)}{\psi }_j\left(x_2\right){\psi }_m\left(x_3\right)}}}$ (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 x₁x₂ and x₃. Strain energy J is minimized to obtain an algebraic equation by taking the derivative of J with respect to $S_{a,ijm}^k$ .

$\frac{\partial J}{\partial S_{a,ijm}}=0$ (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

$\left[\begin{array}{c}{F}^{u_1}\\ F^{u_2}\\ F^{u_3}\\ F^{{\sigma }_{13}}\\ ⋮\\ F^{{ϵ}_{12}}\end{array}\right]=\left[\begin{array}{cccccc}{K}^{u_1{u}_1}& K^{u_1{u}_2}& K^{u_1{u}_3}& K^{u_1{\sigma }_{13}}& \cdots & K^{u_1{ϵ}_{12}}\\ K^{u_2{u}_1}& K^{u_2{u}_2}& K^{u_2{u}_3}& K^{u_2{\sigma }_{12}}& \cdots & K^{u_2{ϵ}_{12}}\\ K^{u_3{u}_1}& K^{u_3{u}_2}& K^{u_3{u}_3}& K^{u_3{\sigma }_{12}}& \cdots & K^{u_3{ϵ}_{12}}\\ K^{{\sigma }_{13}{u}_1}& K^{{\sigma }_{13}{u}_2}& K^{{\sigma }_{13}{u}_3}& K^{{\sigma }_{13}{\sigma }_{12}}& \cdots & K^{{\sigma }_{13}{ϵ}_{12}}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ K^{{ϵ}_{12}{u}_1}& K^{{ϵ}_{12}{u}_2}& K^{{ϵ}_{12}{u}_3}& K^{{ϵ}_{12}{\sigma }_{12}}& \cdots & K^{{ϵ}_{12}{ϵ}_{12}}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ {\sigma }_{13}\\ ⋮\\ {ϵ}_{12}\end{array}\right]$ (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] .

$\left({\bar{\sigma }}_{xx},{\bar{\sigma }}_{yy},{\bar{\sigma }}_{xy}\right)=\frac{1}{{\sigma }_o{\left(a/2\right)}^2}\left({\sigma }_{xx},{\sigma }_{yy},{\sigma }_{xy}\right)$ , ${\bar{\sigma }}_{zz}=\frac{1}{{\sigma }_o}{\sigma }_{zz}$

$\left({\bar{\sigma }}_{xz},{\bar{\sigma }}_{yz}\right)=\frac{1}{{\sigma }_o\left(a/2\right)}\left({\sigma }_{xz},{\sigma }_{yz}\right)$ (18a)

$\left({\bar{u}}_1,{\bar{u}}_2\right)=\frac{E_2}{q_{o}h{\left(a/2\right)}^3}\left(u_1,u_2\right)$ , ${\bar{u}}_3=\frac{E_2}{{\sigma }_{o}h{\left(a/2\right)}^3}$

Similarly, we express the results of the sandwich plate in the following nondimensionalized form according to Moleiro [29]

$\left({\bar{\sigma }}_{xx},{\bar{\sigma }}_{yy},{\bar{\sigma }}_{xy}\right)=\frac{10h^2}{{\sigma }_o{a}^2}\left({\sigma }_{xx},{\sigma }_{yy},{\sigma }_{xy}\right)$ , ${\bar{\sigma }}_{zz}=\frac{1}{{\sigma }_o}{\sigma }_{zz}$

$\left({\bar{\sigma }}_{xz},{\bar{\sigma }}_{yz}\right)=\frac{1}{{\sigma }_{o}a}\left({\sigma }_{xz},{\sigma }_{yz}\right)$ (18b)

$\left(\bar{u},\bar{v}\right)=\frac{10E{E}^{\ast }{h}^2}{{\sigma }_o{a}^3}\left(u,v\right)$ , $\bar{w}=\frac{10E^{\ast }{h}^3}{{\sigma }_o{a}^4}w$

where $E^{\ast }=100\text{GPa}$ , ${\sigma }_o=1\text{MPa}$ .

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 x₁ and the x₂ directions in the bottom and the top layers, respectively.

Case 2; A three-ply laminate [0˚/90˚/0˚] with the fibers parallel to the x₁, x₂, and x₁ 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 ${\bar{\sigma }}_{xx}\left(\frac{a}{2},\frac{a}{2},z\right)$ , and ${\bar{\sigma }}_{yy}\left(\frac{a}{2},\frac{a}{2},z\right)$ , both on B₄, load case a, and a/h = 5, 10, case 1.

(a) (b)

Figure 3. Comparison of present results and Vel [1] for ${\bar{\sigma }}_{xy}\left(\frac{a}{8},0,z\right)$ , on B₃ and ${\bar{\sigma }}_{xz}\left(\frac{a}{8},\frac{a}{2},z\right)$ , on B₆, load case a, and a/h = 5, 10, case 1.

(a) (b)

Figure 4. Comparison of present results and Vel [1] for ${\bar{\sigma }}_{yz}\left(\frac{a}{8},\frac{a}{2},z\right)$ , on B₈, and ${\bar{\sigma }}_{zz}\left(\frac{a}{2},\frac{a}{2},z\right)$ , on B₇, load case a, and a/h = 5, 10, case 1.

(a) (b)

Figure 5. Comparison of present results and Vel [1] for $\bar{U}\left(\frac{a}{2},\frac{a}{2},z\right)$ , on B₁ and $\bar{W}\left(\frac{a}{2},\frac{a}{2},z\right)$ , on B₂, 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 ${\bar{\sigma }}_{xx}\left(\frac{a}{2},\frac{a}{2},z\right)$ , ${\bar{\sigma }}_{yy}\left(\frac{a}{2},\frac{a}{2},z\right)$ , both on B₄ 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 N₁ = N₂ = 6, and N₃ = 4) = 5733. The continuity conditions hold for the fact that the values of ${\sigma }_{13}$ , ${\sigma }_{23}$ and ${\sigma }_{33}$ at the top of the layer are the same as ${\sigma }_{13}$ , ${\sigma }_{23}$ and ${\sigma }_{33}$ 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 N₁ = N₂ = 8, and

(a) (b)

Figure 7. Comparison of present results and Vel [1] for ${\bar{\sigma }}_{xy}\left(\frac{a}{8},0,z\right)$ , on B₃, ${\bar{\sigma }}_{xz}\left(\frac{a}{8},\frac{a}{2},z\right)$ , on B₆ load case a, and a/h = 5, 10, case 2.

(a) (b)

Figure 8. Comparison of present results and Vel [1] for ${\bar{\sigma }}_{yz}\left(\frac{a}{2},0,z\right)$ , on B₈ and ${\bar{\sigma }}_{zz}\left(\frac{a}{2},\frac{a}{2},z\right)$ , on B₇ load case a, and a/h = 5, 10, case 2.

(a) (b)

Figure 9. Comparison of present results and Vel [1] for $\bar{U}\left(\frac{a}{2},\frac{a}{2},z\right)$ , on B₁ and $\bar{W}\left(\frac{a}{2},\frac{a}{2},z\right)$ , on B₂, 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 ${\bar{\sigma }}_{xx}\left(\frac{a}{2},\frac{a}{2},z\right)$ , ${\bar{\sigma }}_{yy}\left(\frac{a}{2},\frac{a}{2},z\right)$ , both on B₅ load case a, and a/h = 4, 10, 50, case 3.

(a) (b)

Figure 11. Comparison of present results and Moleiro [29] for ${\bar{\sigma }}_{xy}\left(0,0,z\right)$ , ${\bar{\sigma }}_{xz}\left(0,\frac{a}{2},z\right)$ , both on B₅, load case a, and a/h = 4, 10, 5, case 3.

(a) (b)

Figure 12. Comparison of present results and Moleiro [29] for ${\bar{\sigma }}_{yz}\left(0,\frac{a}{2},z\right)$ , ${\bar{\sigma }}_{zz}\left(\frac{a}{2},\frac{a}{2},z\right)$ , both on B₅, load case a, and a/h = 4, 10, 50, case 3.

(a) (b)

Figure 13. Comparison of present results and Moleiro [29] for $\bar{U}\left(0,a/2,z\right)$ , $\bar{W}\left(a/2,a/2,z\right)$ , both on B₅, load case a, and a/h = 4, 10, 50, case 3.

N₃ = 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 ${\sigma }_{13}$ 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 ${\bar{\sigma }}_{xx}$ , ${\bar{\sigma }}_{yy}$ , ${\bar{\sigma }}_{xy}$ , ${\bar{\sigma }}_{xz}$ , ${\bar{\sigma }}_{yz}$ , ${\bar{\sigma }}_{zz}$ , ${\bar{U}}_3$ 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.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Vel, S.S. and Batra, R. (1999) Analytical Solution for Rectangular Thick Laminated Plates Subjected to Arbitrary Boundary Conditions. AIAA Journal, 37, 1464-1473. https://doi.org/10.2514/2.624
[2]	Reddy, J.N. (1984) A Simple Higher-Order Theory for Laminated Composite Plates. Journal of Applied Mechanics, 51, 745-752. https://doi.org/10.1115/1.3167719
[3]	Reddy, J.N. (2003) Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. CRC Press, Boca Raton. https://doi.org/10.1201/b12409
[4]	Reddy, J.N. (2019) Introduction to the Finite Element Method. McGraw-Hill Education, New York.
[5]	Pagano, N. and Hatfield, S.J. (1972) Elastic Behavior of Multilayered Bidirectional Composites. AIAA Journal, 10, 931-933. https://doi.org/10.2514/3.50249
[6]	Pagano, N.J. (1970) Exact Solutions for Rectangular Bidirectional Composites and Sandwich Plates. Journal of Composite Materials, 4, 20-34. https://doi.org/10.1177/002199837000400102
[7]	Moleiro, F., Mota Soares, C.M., Mota Soares, C.A. and Reddy, J.N. (2011) A Layerwise Mixed Least-Squares Finite Element Model for Static Analysis of Multilayered Composite Plates. Computers & Structures, 89, 1730-1742. https://doi.org/10.1016/j.compstruc.2010.10.008
[8]	Tauchert, T. and Howard, H. (2012) Least-Squares Residual Solution for Displacement Control of a Composite Piezothermoelastic Cylinder. Journal of Thermal Stresses, 35, 61-74. https://doi.org/10.1080/01495739.2012.637742
[9]	Burns, D. and Batra, R.C. (2021) Static Deformations of Fiber-Reinforced Composite Laminates by the Least-Squares Method. In: Challamel, N., Kaplunov, J. and Takewaki, I., Eds., Modern Trends in Structural and Solid Mechanics 1: Statics and Stability, Publisher, city, 1-16. https://doi.org/10.1002/9781119831891.ch1
[10]	Srinivas, S. and Rao, A. (1970) Bending, Vibration and Buckling of Simply Supported Thick Orthotropic Rectangular Plates and Laminates. International Journal of Solids and Structures, 6, 1463-1481. https://doi.org/10.1016/0020-7683(70)90076-4
[11]	Vel, S.S. and Batra, R. (2000) The Generalized Plane Strain Deformations of Thick Anisotropic Composite Laminated Plates. International Journal of Solids and Structures, 37, 715-733. https://doi.org/10.1016/S0020-7683(99)00040-2
[12]	Reissner, E. (1984) On a Certain Mixed Variational Theorem and a Proposed Application. International Journal for Numerical Methods in Engineering, 20, 1366-1368. https://doi.org/10.1002/nme.1620200714
[13]	Reissner, E. (1986) On a Mixed Variational Theorem and on Shear Deformable Plate Theory. International Journal for Numerical Methods in Engineering, 23, 193-198. https://doi.org/10.1002/nme.1620230203
[14]	Carrera, E. (1998) Evaluation of Layerwise Mixed Theories for Laminated Plates Analysis. AIAA Journal, 36, 830-839. https://doi.org/10.2514/2.444
[15]	Carrera, E., Brischetto, S. and Nali, P. (2011) Plates and Shells for Smart Structures: Classical and Advanced Theories for Modeling and Analysis. John Wiley & Sons, New York. https://doi.org/10.1002/9781119950004
[16]	Carrera, E. (2000) A Priori vs. a Posteriori Evaluation of Transverse Stresses in Multilayered Orthotropic Plates. Composite Structures, 48, 245-260. https://doi.org/10.1016/S0263-8223(99)00112-9
[17]	Carrera, E. (2002) Theories and Finite Elements for Multilayered, Anisotropic, Composite Plates and Shells. Archives of Computational Methods in Engineering, 9, 87-140. https://doi.org/10.1007/BF02736649
[18]	Lage, R.G., Soares, C.M., Soares, C.M. and Reddy, J. (2004) Layerwise Partial Mixed Finite Element Analysis of Magneto-Electro-Elastic Plates. Computers & Structures, 82, 1293-1301. https://doi.org/10.1016/j.compstruc.2004.03.026
[19]	Lage, R.G., Soares, C.M., Soares, C.M. and Reddy, J. (2004) Modelling of Piezolaminated Plates Using Layerwise Mixed Finite Elements. Computers & Structures, 82, 1849-1863. https://doi.org/10.1016/j.compstruc.2004.03.068
[20]	Bathe, K. (1996) Finite Element Procedures. Prentice Hall, Englewood Cliffs.
[21]	Bathe, K.J. (2001) The Inf-Sup Condition and Its Evaluation for Mixed Finite Element Methods. Computers & Structures, 79, 243-252. https://doi.org/10.1016/S0045-7949(00)00123-1
[22]	Pontaza, J. and Reddy, J. (2004) Mixed Plate Bending Elements Based on Least- Squares Formulation. International Journal for Numerical Methods in Engineering, 60, 891-922. https://doi.org/10.1002/nme.991
[23]	Pontaza, J.P. (2005) Least-Squares Variational Principles and the Finite Element Method: Theory, Formulations, and Models for Solid and Fluid Mechanics. Finite Elements in Analysis and Design, 41, 703-728. https://doi.org/10.1016/j.finel.2004.09.002
[24]	Moleiro, F., Mota Soares, C.M., Mota Soares, C.A. and Reddy, J.N. (2010) Layerwise Mixed Least-Squares Finite Element Models for Static and Free Vibration Analysis of Multilayered Composite Plates. Composite Structures, 92, 2328-2338. https://doi.org/10.1016/j.compstruct.2009.07.005
[25]	Moleiro, F., Soares, C.M., Soares, C.M. and Reddy, J. (2008) Mixed Least-Squares Finite Element Model for the Static Analysis of Laminated Composite Plates. Computers & Structures, 86, 826-838. https://doi.org/10.1016/j.compstruc.2007.04.012
[26]	Vel, S.S. and Batra, R. (2002) Exact Solution for Thermoelastic Deformations of Functionally Graded Thick Rectangular Plates. AIAA Journal, 40, 1421-1433. https://doi.org/10.2514/2.1805
[27]	Demasi, L. (2012) Partially Zig-Zag Advanced Higher Order Shear Deformation Theories Based on the Generalized Unified Formulation. Composite Structures, 94, 363-375. https://doi.org/10.1016/j.compstruct.2011.07.022
[28]	Demasi, L. (2007) Three-Dimensional Closed Form Solutions and Exact Thin Plate Theories for Isotropic Plates. Composite Structures, 80, 183-195. https://doi.org/10.1016/j.compstruct.2006.04.073
[29]	Moleiro, F., Mota Soares, C.M. and Carrera, E. (2019) Three-Dimensional Exact Hygro-Thermo-Elastic Solutions for Multilayered Plates: Composite Laminates, Fibre Metal Laminates and Sandwich Plates. Composite Structures, 216, 260-278. https://doi.org/10.1016/j.compstruct.2019.02.071
[30]	Moleiro, F., Carrera, E., Li, G., Cinefra, M. and Reddy, J.N. (2019) Hygro-Thermo- Mechanical Modelling of Multilayered Plates: Hybrid Composite Laminates, Fibre Metal Laminates and Sandwich Plates. Composites Part B: Engineering, 177, Article ID: 107388. https://doi.org/10.1016/j.compositesb.2019.107388
[31]	Mathew, C. and Ejiofor, O. (2023) Mechanics and Computational Homogenization of Effective Material Properties of Functionally Graded (Composite) Material Plate FGM. International Journal of Scientific and Research Publications, 13, 128-150. https://doi.org/10.29322/IJSRP.13.09.2023.p14120