Third-Order Theory for the Bending Analysis of Laminated Thin and Thick Plates Including the Strain Gradient Effect

The aim of the paper is the development of a third-order theory for laminated composite plates that is able to accurately investigate their bending behavior in terms of displacements and stresses. The starting point is given by the corresponding Reddy’s Third-order Shear Deformation Theory (TSDT). This model is then generalized to consider simultaneously the Classical Laminated Plate Theory (CLPT), as well as the First-order Shear Deformation Theory (FSDT). The constitutive laws are modified according to the principles of the nonlocal strain gradient approach. The fundamental equations are solved analytically by means of the Navier methodology taking into account cross-ply and angle-ply lamination schemes. The numerical applications are presented to highlight the nonlocal effects on static behavior.


Introduction
Higher-order plates theories for laminates have been introduced in the last decades to avoid some issues related to the use of lower-order and simpler approaches, such as the Classical Laminated Plate Theory (CLPT) and First-order Shear Deformation Theory (FSDT) [1,2]. In particular, higher-order equivalent single-layer theories allow to obtain a more accurate interlaminar stress analysis and do not need to introduce the shear correction factor [3][4][5][6]. The displacement field that characterizes the Third-order Shear Deformation Theory (TSDT), for instance, determines a quadratic profile of shear strains and stresses along the thickness [7,8], due to its cubic expansion in the thickness coordinate. Consequently, there is no need for the shear correction factor [9][10][11][12]. The importance of these cubic terms in the analysis of laminates has been recently highlighted in the paper by Petrolo and Carrera [13], in which the best theory diagrams for multilayered structures have been widely discussed.
In general, higher-order approaches have been justified by the use of more and more advanced materials [14][15][16][17] and the need of innovative configurations for the optimal design of structures [18,19]. In particular, their introduction could be essential when these innovative constituents are included in the stacking sequences of multilayered or sandwich structures [20][21][22][23][24][25]. These aspects have been clearly emphasized in the works by Carrera [26][27][28], Carrera and Giunta [29], and Carrera et al. [30]. In these works, accurate and effective higher-order structural models based on a unified formulation have been presented.
The increasing number of applications involving micro-and nanostructures [31][32][33][34] has proven that the size-dependent features of the advanced constituents could have not negligible effects on their mechanical behavior [35][36][37]. These aspects have been highlighted by so-called multiscale analysis [38][39][40][41]. In these circumstances, structural theories based on classical elasticity could turn out to be inadequate to model such innovative mediums [42][43][44][45][46][47]. Nonlocal theories have been developed to overcome these issues, as illustrated in the first papers by Eringen [48,49]. Some other subsequent contributions could be mentioned to this aim. For instance, Luciano and Willis investigated the nonlocal constitutive behavior of an infinite composite laminate [50]. Nanorods made of functionally graded materials were studied by Barretta et al. [51], who included the gradient Eringen model in their framework. The same approach was used for Timoshenko nanobeams in bending [52]. Analogously, the effects of nonlocal elasticity were emphasized in the papers by Apuzzo et al. [53,54] for the investigation of torsional behavior and dynamic response of Bernoulli-Euler nanobeams. More recently, similar remarks and observations have been drawn in [55][56][57]. In general, different nonlocal theoretical frameworks can be developed.
Examples of nonlocal approaches that should be mentioned for completeness purposes are the strain and stress gradient theories [58][59][60], the modified couple stress theory [61][62][63][64], and the ones based on micropolar formulations [65][66][67]. A comprehensive literature review concerning nonlocal elasticity can be found in the paper by Zhao et al. [68] for completeness purposes.
The current paper is developed within the strain gradient theoretical framework [69]. According to this approach, the size effects of the materials are modeled by means of a length-scale parameter, which controls the second-order derivatives of the strain components [70]. Consequently, the stresses are functions not only of the strains in an evaluation point, but also depend on the divergence of the gradient of the strains. These aspects have been emphasized in the papers by Aifantis [71] and by Askes and Aifantis [72], where an overview of gradient elasticity formulations in statics and dynamics has been presented. Therefore, higher-order derivatives of displacements are involved in the threedimensional constitutive laws [73][74][75][76]. These aspects have been also recently emphasized in papers [77,78] where peculiar Finite Element (FE) formulations have been developed to take into account the strain gradient effect. In particular, the use of higher-order Hermite interpolating polynomials for the approximation of both membrane and bending degrees of freedom have been disclosed.
The plate theory that this nonlocal effect is included in is based on the TSDT for laminates [1], in which the various layers that define the stacking sequence have orthotropic features [79,80]. The kinematic model is written to take into account simultaneously lowerorder approaches, such as the CLPT and FSDT, for comparison purposes. After the brief introduction of the main topics of the paper (Section 1), the theory is presented by using a matrix compact notation in Section 2. Here, the strong form of the governing equations is obtained, once the definitions of both strains and stress resultants are carried out. The fundamental system is solved analytically by means of the Navier approach as shown in [73]. The algebraic system of equations is written in Section 3, highlighting the contributions related to the strain gradient effect for cross-ply and angle-ply simply supported laminated composite plates. As remarked in [1], the Navier methodology can be efficiently applied to deal with these configurations. Section 4 is focused on the numerical applications. The proposed approach is validated for both classical and nonlocal elasticity through comparison with the results available in the literature, taking into account the three different structural models. Then, the results are extended to emphasize the influence of the nonlocal parameter on the static behavior, which is expressed in terms of displacements and stress components. The through-the-thickness stress profiles are also provided. Finally, the main achievements are drawn in Section 5. Appendix A collects instead the analytical expressions of the terms involved in the algebraic formulation for cross-ply and angle-ply laminates.

Nonlocal Structural Model
The theoretical framework is developed in this Section for a rectangular plate whose planar size is a × b, considering a Cartesian reference system xyz. The plate is made of a sequence of N L orthotropic layers. The generic k-th ply has a thickness h k = z k+1 − z k , with z k+1 , z k being its upper and lower thickness coordinates. The overall plate thickness is given by h = ∑ N L k=1 h k . These geometric features are all shown and specified in Figure 1. The three-dimensional displacements collected in U(x, y, z) = U V W T can be written in terms of the five degrees of freedom, which are three translations u, v, w and two rotations φ x , φ y . The vector u(x, y) = u v w φ x φ y T is conveniently introduced.
The displacement field assumes the following compact aspect: where D (0) is a differential operator given by whereas the matrices I (0) , I (1) , I (3) are defined as The choice of the structural theory defines the values of the constant c 1 and the thickness function F (z). In particular, the TSDT is obtained for c 1 = 4 3h 2 and F = z 3 . By setting c 1 = 0, the FSDT is achieved instead. On the other hand, the CLPT can be defined by using c 1 = 1 and F = z. The fundamental assumptions of each theory, therefore, are different according to this choice [1]. It should be specified that only the FSDT requires the shear correction factor in the definitions of the shear stresses. The value of 5/6 is considered to this aim.
The membrane strain components ε = ε xx ε xy γ xy T and the transverse shear strains γ = γ xz γ yz T can be written as follows: where F = ∂F ∂z . The terms introduced in (4) are discussed below. According to the notation employed by Reddy [1], the membrane strains ε (0) assume the following aspect: On the other hand, the curvatures ε (1) are computed as follows: where the differential operator D (m) is given by Higher-order membrane strains vector ε (3) is defined as in which the following definition is used for the differential operator D (b) : The shear strains are now discussed. The components collected in γ (0) can be written as whereas the differential operator D (0) s is given by Likewise, higher-order shear terms included in γ (2) are defined as The governing equations are derived by means of the principle of virtual displacements [1] δΦ + δL = 0, (14) in which δΦ is the strain energy variation, whereas δL represents the work done by applied external forces. If a laminated composite plate made of N L layers is considered, the variation δΦ can be defined as in which A denotes the plate middle surface. The membrane stresses in the k-th layer are specified by σ (k) = σ and assume the following aspect [77,78]: where is the nonlocal parameter linked to the influence of the micro/macroscale interactions, ∇ 2 = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 is the Laplacian, andQ (k) is the plane-stress-reduced stiffness coefficients matrix given byQ Its terms are computed as a function of the Young's moduli E 1 , E 2 ; Poisson's ratio ν 12 ; and shear modulus G 12 of the orthotropic medium [1]. The constitutive law (16) takes into account the strain gradient effect. Analogously, a similar relation can be written for the in whichQ (k) s is the stiffness coefficients matrix related to the shear shown below whose terms are computed as a function of the shear moduli G 13 , G 23 [1]. It is important to highlight that both membrane and shear stresses are characterized by a classical and a nonlocal part, which is the one multiplied by . The external loads can be collected into the vector q = q x q y q z M x M y T , which includes five load components. Consequently, the work done by external forces δL can be written as The governing equations can be obtained by conveniently introducing the stress resultants as the integrals of the stress components along the thickness of the layer. The following quantities are defined: The system of five differential equations in terms of stress resultants that governs the static behavior of the plates is carried out by performing the proper manipulations starting from the principle of virtual displacements [77]. By using a compact matrix form, it becomes where Boundary conditions in terms of primary and secondary variables are obtained as shown in the book by Reddy [1], since they are not affected by the strain gradient effect.
It is convenient at this point to express the stress resultants in the fundamental Equation (22) as a function of the displacement vector u, recalling the constitutive laws (16) and (18), as well as the definitions of the strains shown in (4). The following relations are carried out, as far as the membrane stress resultants are concerned: where the terms into the constitutive matrices A, B, D, E, F, H are given by for i, j = 1, 2, 6. The following differential operators that appear due to the Laplacian are also introduced: D On the other hand, the stress resultants related to shear forces are defined as follows: in which the terms collected in the constitutive matrices A s , L s , N s , for i, j = 4, 5, are given by where κ s is the shear correction factor. Its value is different from the unity only for the FSDT, in which it is assumed equal to 5/6. The differential operators D

Solution Procedure
Once the nonlocal governing equations are written in terms of the displacements collected in u, they can be solved analytically by means of the Navier methodology. As illustrated in [1], this approach can be applied only for some peculiar lamination schemes, which are antisymmetric cross-ply and antisymmetric angle-ply, respectively. These two cases are analyzed separately in the following, assuming simply supported boundary conditions for both circumstances. The solution to the current static problem is provided by the algebraic linear system shown below in which ∆ = U mn V mn W mn X mn Y mn T is the vector that collects the unknown coefficients that determine the displacement amplitudes. On the other hand, K and F are the stiffness matrix and the load vector, respectively. The terms included in the symmetric matrix K will be specified for each lamination scheme. On the other hand, the load vector has the following definition, assuming that only transverse surface forces are applied: having introduced the same expansion also for the applied external forces [1].

Cross-Ply Laminates
In order to apply the Navier approach for the cross-ply sequence at issue, the stiffnesses written below are all equal to zero: In this circumstance, the solution can be sought assuming the following expansion: cos αx sin βy sin αx cos βy sin αx sin βy cos αx sin βy sin αx cos βy where • stands for the elementwise product and α = mπ/a, β = nπ/b. The explicit expressions of the coefficients K ij included in the stiffness matrix K, for i, j = 1, . . . , 5, whose size is clearly 5 × 5 are listed in Appendix A, recalling that K ij = K ji .

Angle-Ply Laminates
The Navier approach can be applied for angle-ply laminates if the following stiffnesses are all equal to zero: In this case, the solution is obtained assuming the following expansion: sin αx cos βy cos αx sin βy sin αx sin βy cos αx sin βy sin αx cos βy The expressions of the coefficients K ij , for i, j = 1, . . . , 5, of the stiffness matrix K, are listed in Appendix A, presenting only those terms that are different from the ones valid for cross-ply sequences.
The solution in terms of ∆ can be easily obtained from Equation (32), for both laminates under consideration. Once the amplitudes in ∆ are computed, definitions (34) and (36) allow us to obtain the displacements within the reference domain. Consequently, the strains can be also deduced. The membrane stress components σ (k) can be evaluated as well, following the procedure illustrated in the book by Reddy [1] through the constitutive relations previously presented. On the other hand, the shear stress components τ (k) are determined from the three-dimensional equilibrium elasticity equations [1,73]. The complete procedure is omitted for conciseness purposes.

Numerical Results
The current Section aims to present the results of the static analyses. Due to the general features of the theoretical approach, the solutions are presented for different nonlocal theories, which are CLPT, FSDT and TSDT, setting properly the values of c 1 and F . As far as the mechanical features are concerned, the ratio between the longitudinal and transverse Young's moduli E 1 /E 2 is specified in each application, whereas the other quantities are taken as G 12 = G 13 = 0.5E 2 , G 23 = 0.2E 2 , ν 12 = 0.25 for Material 1, or G 12 = G 13 = 0.6E 2 , G 23 = 0.5E 2 , ν 12 = 0.25 for Material 2. The lamination schemes, instead, are denoted by (θ (1) / . . . /θ (k) / . . . /θ (N L ) ), where θ (k) stands for the orientation of the k-th layer. The results are presented in dimensionless form. In particular, the central deflectionw is given byw On the other hand, the stress components are evaluated as follows, unless differently specified: It should be specified that the values of the stresses presented in this Section are all related to the classical component, which can be deducted from definitions (16) and (18) following the approach used in [73,77]. The analyses are carried out for increasing values of the dimensionless nonlocal parameter ( /a) 2 , in order to show the effect of the strain gradient on the static solutions.
The first application aims to investigate the central deflectionw as a function of side-to-thickness ratio a/h of a square plate for different lamination schemes: crossply (0/90/90/0) and angle-ply (45/−45). The results are shown in Figure 2, assuming E 1 /E 2 = 25 for the cross-ply (Material 1) and E 1 /E 2 = 40 for the angle-ply (Material 2). The graphs include the structural response obtained by means of the three theories. Each model is related to a different line style: solid line for TSDT, dotted for FSDT, and dash-dotted for CLPT. The color, instead, is linked to the value of the nonlocal parameter. The same choice is also kept in the following figures. It can be observed that the greater is the nonlocal effect and the lower is the vertical deflection, independently from the theory. In other words, the central deflection is reduced by increasing the value of the nonlocal parameter ( /a) 2 . The FSDT and TSDT, moreover, are characterized by a comparable behavior and are highly affected by plate thickness. By increasing the ratio a/h, their displacements tend to the results of the CLPT, which do not depend on that geometric ratio. The corresponding curves, in fact, are described by rectilinear functions. Similar behaviors are obtained in both lamination schemes.
In the next test, a (0/90/0) cross-ply square plate is considered. Material 1 is taken into account in this circumstance, assuming E 1 /E 2 = 25. The results are presented in Table 1 for different values of a/h, varying the structural theory. Where available, the analytical solutions by Reddy [1] are provided, in terms of both displacement and stress components. The reference results are clearly presented only for classical elasticity, assuming ( /a) 2 = 0. The comparison proves a very good agreement between the current approach and the reference one.  Table 2, instead, aims to extend these results to the nonlocal elasticity framework, varying the ratio ( /a) 2 . The same reduction of the displacement values observed in Figure 2 can also be seen numerically in this circumstance. The three theories provide really close results for thin plates even if the strain gradient effect is taken into account. It should be noted from these tables that noticeable differences can be seen especially in terms of membrane stresses in thick configurations (defined by a/h = 4), if the results related to the different theories are compared. As it will be highlighted in the following paragraphs, this is due to the fact that the TSDT is characterized by nonlinear stress profiles. This nonlinearity is particularly emphasized for thicker plates, whereas it is reduced for lower values of a/h. A (0/90/90/0) cross-ply lamination scheme is considered in the next application. Even in this case, Material 1 is taken into account to describe the orthotropic features of the layers, with E 1 /E 2 = 25. The results are shown in Table 3 for the classical elasticity. In the same Table, the solutions shown in [1] are presented for comparison purposes. As in the previous test, the present configuration is analyzed also in the nonlocal framework. Results for increasing values of ( /a) 2 are included in Table 4 for the sake of completeness. The reference solutions are all related to classical elasticity. In the following application, the current methodology is also verified with regard to nonlocal elasticity, taking into account the solutions presented in [77]. To this aim, only CLPT is considered in accordance with the reference paper, due to the availability of the literature. Several cross-ply lamination schemes are considered, assuming E 1 /E 2 = 25 and Material 1 as orthotropic features. The results are presented in Table 5 for various stacking sequences, varying the value of ( /a) 2 , in terms of central deflections and membrane stress components. In this circumstance, the following dimensionless quantities are considered for the stresses: whereas the same expression introduced before is used forσ xx . The results shown in Table 5 prove that the present solutions are in very good agreement with the reference ones, also in the framework of nonlocal elasticity. The next application is focused on the bending analysis of antisymmetric angle ply laminates, characterized by (θ/−θ/ . . . ) as the lamination scheme and θ being the arbitrary orientation of each layer. Firstly, the analyses are presented in terms of dimensionless central deflectionw, considering square plates made of Material 2 and E 1 /E 2 = 40 as orthotropic ratio. The results are presented for several side-to-thickness values a/h, orthotropic angles θ, and different structural theories in Table 6, as far as classical elasticity is concerned, which means ( /a) 2 = 0. The values are in good agreement with the ones taken as references [1].
As in the previous cases, the analyses are extended to nonlocal elasticity by increasing the value of ( /a) 2 but keeping the same geometric ratios and mechanical features. The results are presented in Table 7. Even in these circumstances, the differences between TSDT and FSDT decrease for lower values of thickness, and the displacements tend to the ones of the CLPT.
A (−45/45) 4 laminated plate square plate is considered in the next application. Each orthotropic layer is made of Material 1, with E 1 /E 2 = 25. This test aims to evaluate the stress components in angle-ply configurations. With respect to the dimensionless values shown in (38), only the transverse shear stressσ xz is specified in a different thickness coordinate as specified belowσ It should be recalled that in the following applicationσ xx =σ yy andσ xz =σ yz , therefore, their values are not repeated twice. The results related to the classical elasticity are shown in Table 8, varying a/h and the structural theory. Even if the displacements become closer independently from the considered approach, the stresses assume different values especially if the thickness is greater. In this circumstance, the reference solutions are available only for the FSDT [1]. Table 6. Dimensionless central displacement 10 2w of a angle-ply square plate with (θ/−θ/ . . . ) as lamination scheme, for ( /a) 2 = 0. The number of plies is denoted by N L .     On the other hand, the nonlocal counterpart of the current application is shown in Table 9. The increasing values of ( /a) 2 emphasize the differences in terms of stresses, especially if the plates are thicker. The last tests aim to present the stress analysis in terms of the through-the-thickness distributions of the various components, highlighting the differences that could arise by varying the structural approach (TSDT, FSDT, and CLPT) for different geometric ratios. The effect of ( /a) 2 is also investigated. Firstly, a (0/90/90/0) cross-ply square plate is analyzed, for a/h equal to 4 and 10, respectively. The orthotropic features of the layers are obtained by setting E 1 /E 2 = 25 and selecting Material 1 as constituent. The membrane and shear dimensionless stresses are given bȳ wherez = 2z/h stands for the dimensionless thickness coordinate. The through-thethickness stress distributions in Figure 3 are related to the case a/h = 4, which is representative for thick plates. It can be observed that in both classical and nonlocal elasticity, the TSDT is characterized by significantly different profiles due to the higher-order features of the displacement field. This aspect is more evident in the membrane stress components, which are obtained by means of the application of constitutive laws. In fact, it should be recalled that the TSDT allows a cubic representation of the stress profiles, whereas a linear variation is associated with the FSDT and CLPT. This feature gives rise to noticeable differences if thick plates are investigated, as it can be seen from the plot ofσ xx in Figure 3. The variation of the transverse stresses, instead, is always characterized by curved and continuous profiles since they are equilibrium-derived, according to the procedure shown in [1]. A lower value of thickness is considered in the graphs shown in Figure 4, assuming a/h = 10. Due to this choice, the stress distributions tend to the same value, for both ( /a) 2 = 0 and ( /a) 2 = 0.10, independently from the theory. Therefore, the differences among the theories is practically negligible starting from a/h = 10, which is typically the geometric ratio that characterizes moderately thick plates.  A similar analysis is carried out for a (−45/45) angle-ply square plate, for a/h equal to 4 and 10, respectively. In this case, the orthotropic features of the layers are given by E 1 /E 2 = 40. Material 2 is taken into account as constituent. With respect to the dimensionless expressions introduced in (41), a different value ofσ yz is considered, which is defined below:σ The graphical plots are shown in Figures 5 and 6 for a/h = 4 and a/h = 10, respectively. It can be observed that in these circumstances, the FSDT and CLPT are overlapped. The nonlinear behavior of TSDT is clearly more evident for thicker plates, for both classical and nonlocal elasticity.

Conclusions
A theoretical framework able to simultaneously and accurately study thick and thin laminated composite plates has been presented in the paper. In particular, the chosen kinematic model is able to describe several theories if properly set, which are the CLPT, FSDT, and TSDT. The theories have been modified to include the strain gradient effect, in order to take into account nonlocal contributions in the evaluation of stresses. The proposed approach is general, and the fundamental static equations have been written for arbitrary configurations by using a compact matrix notation. In order to provide an analytical solution, the Navier methodology has been applied to deal with cross-ply and angle-ply lamination schemes. The explicit definitions needed in the solution procedure have been presented, highlighting also the contributions linked to the strain gradient effect. The approach has been validated numerically by means of comparison with the results accessible in the literature, if available, for both classical and nonlocal elasticity. The results have been presented in terms of displacements and stresses, varying the geometric features, mechanical properties, lamination schemes, as well as the influence of the nonlocal effect. The differences that have been observed among the theories have been emphasized. In particular, the results have proven that the stress components are noticeably different by changing the plate theory, especially for higher values of thickness. Finally, it should be specified that the proposed solutions could be used for further advancements of this topic in nonlocal elasticity and could be taken as benchmarks for future comparisons, especially if numerical methods are developed.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Algebraic Form of the Stiffness Matrix
As specified in (32), the static problems at issue are governed by the algebraic linear system K∆ = F. The terms included in the symmetric matrix K are presented in this Appendix, for the two cases investigated in the paper.
Appendix A.1. Cross-Ply Laminates As far as cross-ply lamination schemes are concerned, the terms of the first row of the matrix K are given by K 13 = −c 1 E 11 α 3 + (E 12 + 2E 66 )αβ 2 + 2 E 11 α 5 + α 3 β 2 + (E 12 + 2E 66 ) α 3 β 2 + αβ 4 , (A3) The terms of the second row are defined as As far as the terms of the third row are concerned, one gets Instead, the terms of the fourth row are given by (A14) Finally, the terms included in the last row assume the following aspect: (A15)