A Recursive Methodology to Determine the Mechanical Response of Thin Laminated Plates in Bending

AbstrAct: The paper’s objective is to present the development of a recursive methodology which is based on Adomian Decompostion Method in order to evaluate the mechanical response of thin laminated plates in linear bending. By the equivalent layer concept, the linear relation between the equivalent stresses and the strains, namely ABD matrix, is established. When viewed by the Adomian Decompostion Method perspective, it generates an interesting idea: each layer influence on the plate’s response can be recursively inserted into a base solution by a rearrangement of the plate’s properties. This base solution is previously obtained and it is, in this paper, an isotropic plate response for the same loading and boundary conditions. This approach can significantly increase optimization and delamination studies, given the simplicity on the layers modification, both in fiber orientation and constitutive properties, as these are considered on the recursive procedure. The pb-2 Rayleigh-Ritz Method is used to approximate the solution space and to generate analytic response surfaces. The methodology is applied to symmetrical and unsymmetrical stacking cases for different boundary conditions sets and loading types and the obtained responses are compared to those found on the literature. A study of case complements the methodology analysis: a simplified landing gear door is modeled considering a set of loading conditions as well as different stacking configurations. Good correspondence was found in all studied cases.


INTRODUCTION
Nowadays and mainly in the aerospace, automotive, and marine industries, sandwich and laminated plates are applied in several types of structures due to their mechanical properties (Gay et al. 2003;Mangino et al. 2007).It is well-known that sandwich plates are designed either to protect the inner layer from the external environment/types of loading or to enhance the inertia of the structure, while laminated plates are designed to maximize the strength/weight ratio as well as to fit properly the structural element to the problem's loading and boundary conditions.Despite the knowledge on sandwich and laminated plate mechanical behavior (Reddy 1997;Ferreira et al. 2011;Carrera 1996;Wang et al. 2002;Vel and Batra 2000;Altenbach 1998), there is great interest in new methodologies which can not only solve the problem in a better fashion but also introduce new information about what is to be analyzed.
Several authors have worked on the determination of sandwich and laminated plates' mechanical behavior.Many mathematical models have been discussed and various physical phenomena have been added to them.Reddy (1997) has described many principles related to this topic, from the Classical Laminated Plates Theory (CLPT), herein used, to non-linear analysis of laminated shells.Ferreira et al. (2011) have analyzed thin and thick isotropic and cross-ply plates using a Unified Formulation, which has been developed by Carrera (1996), presented numerical results for plate static deflection and free vibration.Wang et al. (2002) have studied thick rectangular composite plates using the reproducing kernel particle method and shown good results for linear bending, free vibration, and buckling loads for several boundary conditions sets.Vel and Batra (2000) have applied the Eshelby-Stroh formalism in order to analyze the generalized plain strain in anisotropic linear elastic laminated plates and presented normal and transverse stress distributions for clamped-clamped, cantilever (clampedfree), and clamped-simply supported boundary conditions.Altenbach (1998) has presented a complete review in the state-of-art of theories for laminated and sandwich plates, prior to 1998.In spite of the years, his study is still relevant: not only the collection of several aspects and assumptions of analyzing inhomogeneous plates has been presented but also the advantages and disadvantages of using such structures in engineering have been widely discussed.
The Adomian Decomposition Method (ADM) (Adomian 1994) is a recursive methodology which has been used in several scientific areas to solve linear and non-linear ordinary and partial differential equations.It is well-known that ADM is a special case of homotopy methods (Öziş and Yildirim 2008;Li 2009), which is basically the determination of the solution by applying perturbations to the response of a similar problem.Essentially, in ADM, the differential operator is decomposed into 2 or 3 terms -namely linear, remainder, and non-linear -and its solution is expanded in infinite series.Recursively, each solution term is determined by a correlation among the operators parts and previous results.Rao (2010) has determined the Ricatti generalized differential equation solution by applying the method.Karimpour and Ganji (2008) have employed the ADM and obtained the exact solution for lateral motion of thin plates subjected to lateral and in-plane loadings.Given the orthogonal nature of the basis functions, some of the obtained solutions are analytic.Biazar et al. (2004) have demonstrated that this decomposition can be applied to ordinary non-linear equations and they obtained good results.Cheniguel and Ayadi (2011) have used the method in order to solve the heat equation with non-local boundary conditions.For some problems, the authors have achieved the exact solution of the problem.Olivares (2003) has applied the method to partial differential equations and obtained excellent solutions, when compared to those found in the literature.As one can notice, the method was utilized in different types of problems, differing in their physical areas, and for all of them good results were obtained when compared to other methodologies.
The pb-2 Rayleigh-Ritz Method (pb-2 RRM) is widely used in the solution of the thin and thick plate's differential equation.Basically, a weighted kinematically admissible interpolation basis is constructed to approximate the problem's degrees-of-freedom.By the minimization of the Total Potential Energy, the weight of each function is determined.The pb-2 modification eases the interpolation functions necessity of being kinematically admissible by multiplying the entire basis by specific monomials (Bhat's beam functions; Bhat 1985).These functions insert zeros into the interpolation functions, forcing them to satisfy the boundary conditions set.Bhat (1985) and Liew (1992) used the method for obtaining the transverse displacement response of thin isotropic plates.Improved by Liew and Wang (1993), the methodology has been also considered to solve plate free vibration and to analyze buckling loads.Liew and Lam (1991) employed it in order to determine isotropic and anisotropic free vibration responses.The main objective of using pb-2 RRM is to ease the employment of boundary conditions sets and to obtain global semi-analytical solutions.
As aforementioned, one step of ADM is the decomposition of the differential operator into 2 (or 3) terms.Herein, this decomposition follows the idea of an additive decomposition of the constitutive tensor (Browaeys and Chevrot 2004) which is based on a constitutive hierarchy (Forte and Vianello 1996;Chadwick et al. 2001;Ting 2003).The concept of constitutive symmetry relies on the relationship between the set of planes of 2 different materials, for instance.If the union of such sets is complete, these materials belong to the same symmetry.Following a similar consideration, the idea of constitutive hierarchy is related to the set of planes of 2 different symmetries so that if one is a subgroup of another, a hierarchic relation (with lower and higher symmetries) can be devised.Thus, due to the fact that there is a finite number of constitutive symmetries, one can assembly a constitutive hierarchic tree (Chadwick et al. 2001;Ting 2003).Being the constitutive relationship reduced, in the case of classical plate's theory, a new hierarchic tree is required and it is presented herein.
The additive decomposition extracts a lower symmetry (with more symmetric planes) from a higher one (with fewer symmetric planes) by applying projections into the reduced constitutive tensor.Basically, these projections insert symmetric planes into the set of the lower symmetry resulting into a new constitutive tensor which corresponds to a higher symmetry.The original symmetry then is reconstructed by adding the remainder and the projected tensor.Independently from the symmetry of the original tensor, an isotropic tensor can be always extracted.The reason is simple: the isotropic A Recursive Methodology to Determine the Mechanical Response of Thin Laminated Plates in Bending symmetry is irreducible to any anisotropic symmetry by being the highest symmetry in terms of constitutive hierarchy (Forte and Vianello 1996;Chadwick et al. 2001;Ting 2003).Consequently, one can extract an isotropic portion from all laminae and the remainder tensor, for each laminated, concerns to anisotropic terms.
This methodology can greatly accelerate optimization processes in laminated plates due to the fact that the resultant of the anisotropic tensor is not inversed in the recursive procedure, easing the development of the derivatives with respect either to the plate's stacking or to the fiber orientation.The additive decomposition (Browaeys and Chevrot 2004) is developed to 3-D constitutive tensor.In order to apply this procedure to the reduced constitutive tensor, a new projector is presented.
The paper's aim is to analyze a new methodology which solves the governing equation of laminated thin plate in bending by inserting the material's heterogeneity in a recursive fashion.The thin laminated plate's differential equation is derived and in the process, the well-known ABD matrix is developed.This matrix relates equivalent stresses and the plate's strains.Further, the internal energy as well as the external potential are determined to construct the structure total potential energy functional, necessary to pb-2 RRM.A base solution -an isotropic one -is determined and any modification to the fiber orientation or/and to the lamina stacking is viewed as a modification on the plate's stiffness (ABD matrix) and it is recursively inserted into the isotropic solution by ADM.
This paper is divided in 5 sections.The second one describes the process to derive the diff erential equation for sandwich and laminated thin plates.Th e implementation of pb-2 Rayleigh-Ritz to the CLPT is presented in the third section.In the fourth one, Adomian Decomposition Method is described and applied to the pb-2 RRM, resulting in the recursive procedure.Aft erwards, the numeric results of 5 examples are shown and discussed in the fi ft h section.Finally, the conclusions are presented.In this paper, both matrix and index notation are used.Th e summation with repeated indices are considered, otherwise written, as well as the ranges of Greek and small-caps Roman letters are from 1 to 2 and 1 to 3, respectively.In case of high-caps Roman letters, the range is given at the fi rst opportunity.Vectors and matrices are written respectively in small-caps and high-caps bold letters and all vectors are column-vectors.

Reference surface
where x i represent the Cartesian coordinates; u i denote the plate's displacement in i direction and x 0 α correspond to the membrane displacement.
In Fig. 1, h denotes the plate's thickness, n the outwards normal vector of the plate's boundary дΛ.Th e plate's volume and mid-surface are represented respectively by Ω and Λ, while the lateral loading corresponds to q 3 .Moreover, the plate's stacking is shown, presenting the inhomogeneity in x 3 direction.By the equivalent layer concept hypothesis (Altenbach 1998), the displacement field of a laminated plate is the same as a homogeneous one.Th us, the strain vector, ε, can be derived as: (1) (2) where: noting that R is a 3 × 3 transformation matrix and θ (l)  is the directive angle between the global and local Cartesian coordinate systems of each lamina.
Th e stresses resultants, namely moments -M and normals -N, are defi ned as: ε 0 and κ correspond to the membrane strains and the plate's curvatures, respectively.Th e stresses fi eld, σ, can be introduced as: where l is a positive integer which identifies the laminae (Fig. 1); Q (l) is the reduced constitutive tensor transformed to the global coordinates.Given the plate's inhomogeneous nature, each lamina develops its own stress field.When described in the lamina's principal coordinate system (local coordinates), the reduced constitutive tensor is written as: which represents an orthorhombic material.Being the analysis based on the plate's coordinate axes, the constitutive tensor in Eq. 4 is transformed by the following rule (Reddy 1997): where L is the number of laminae; z l denotes the Euclidian distance between l th lamina's base and the mid-surface.Equations 6 can be described in a compact and well-known form, which is: Th e laminated stiff ness matrix, P, or just ABD matrix, has some intrinsic properties: A and D are respectively the membrane and transverse stiff ness along with B, which is the extensional-bending matrix.
With the stress resultants defi ned, the diff erential equations of classical laminated thin plates can be derived.Translating the equilibrium equations in the stress resultants, M and N, one can fi nd (Reddy 1997): with the differential equations and its boundary conditions, one can choose the most appropriate method to solve the problem.The objective herein is to apply the ADM so as to improve the analysis of sandwich and laminated plates.
The pb-2 RRM uses the Principle of Minimal Potential Energy (PMPE) in order to derive the structure's equilibrium position.The Potential Energy, Π, is determined by the sum of the internal energy, U, and the external potential, V e .The first is described as the inner-product between the stress and strain vectors.As a result, one can write: and by Eqs. 2 and 7: where q m = {q 1 q 2 } T .
In considering the coupled problem, the full diff erential equations system is determined as: where Th e external potential is determined as:

PB-2 RAYLEIGH-RITZ METHOD
By the pb-2 RRM, the plate displacement is interpolated using as base independent kinematic admissible functions.The modification, named pb-2, simplifies the boundary conditions employment by the multiplication of special monomials over the entire base.So as to avoid the manipulation of different family functions and to use the Gaussian quadrature, polynomials are chosen to populate the interpolation basis.Given that, one can write the interpolated displacement field as: The boundary conditions with respect to the displacement and the stress resultants are exactly the same as in CPT (Reddy 1997), due to the equivalent layer concept.Consequently, where ω 1 , ω 2 and ω 3 group the interpolation functions, already adjusted to satisfy their respective boundary conditions; λ 1 , λ 2 and λ 3 are their respective weighting vectors.Moreover: logy, which might not be the best approach in this case, the remainder term can be used to modify either the plate's stacking or the laminae fi ber orientation leading into improvements in optimization routines.Nevertheless, ADM is proposed to solve non-linear diff erential equations.By understanding its properties in the linear analysis as well as its requirements for convergence, one can propose further and better explanations over non-linear analysis of any type.
Th ree main characteristics of ADM play signifi cant role in the methodology.Given the linear diff erential L(д)u = f, these characteristics can be described as: 1. Th e decomposition of the diff erential operator: The vectors' size, n, is related to the basis enrichment and it is set in order to achieve the desired convergence degree.Furthermore, n is derived from the highest polynomial degree, n p , in the interpolation functions by n = (n p + 1) (n p + 2)/2.Concerning the pb-2 RRM, corresponds to the variation on the functional.By applying Eq. 14 into the internal energy, Eq. 12, and the external potential, Eq. 13, one can obtain: where B = ∂W Via the PMPE, one can determine the equilibrium position in which the potential energy is minimal.Th is functional and its fi rst variation is determined by: 2. Th e expansion of the solution terms in an infi nite series: 3. Th e determination of each solution's term in a recursive manner: (15) Inserting Eq. 15 into Eq.16, one can obtain:

ADOMIAN DECOMPOSITION METHOD
Many approaches are possible in the solution of Eq. 17.One of them is the ADM.Despite being a recursive methodo-Depending on how the decomposition of the differential operator is performed, the recursive system may have some physical meaning.The proposition is to decompose the constitutive tensor, in the lamina's level, into 2 components: 1 isotropic and 1 remainder.Then, considering it in the plate's level, an isotropic and a remainder part is considered, which corresponds exactly to the ADM mathematical form.Thus, one can rewrite the reduced constitutive tensor as: where (1) and ( 2) are the decomposed constitutive tensor.
As on can note, Eq. 19 is related to the diff erential operator decomposition of ADM, since the internal energy is linear to the constitutive properties.In Eq. 19, P (1) must be positive-defi nite in order to allow the matrix inversion.Th is is not required for P (2) .Actually, none of the properties are required on P (2)  for the procedure.Th e transverse displacement can be expanded as an infi nite series (proposition 2) as: which, in other words, results in the pb-2 RRM constant vector expansion.
Equation 20 is valid only if the same interpolation basis is used for all expanded terms.Nonetheless, this is not a simplifi cation: it is a requirement for all functions to obey the imposed boundary conditions.Lastly, the system is reassembled as: Th e recursive procedure stops when the convergence is reached.It is easy to notice that the convergence criterion depends on the eigenvalues of [N (1) ] -1 N (2) due to the fact that the system in Eq. 22 can be regarded as a pure stretching/ contraction of λ (m) : by the singular value decomposition (SVD), the eigenvectors transform spatially the constant vectors λ (m) and λ (m-1) , turning Eq. 22 into an independent ordinary system of variables.As a result, the convergence can be reached if and only if the absolute value of each eigenvalue [N (1) ] -1 N (2) of is smaller than the unity.On other hand, if just one eigenvalue is larger than the unity, this value will tend to infi nity, leading the series into a divergent one.Furthermore, the closest the eigenvalue is to the unity, slower the convergence will be.Th e resultant matrix of [N (1) ] -1 N (2) is not defective, given the properties of the interpolation basis and the positivity of N (1) .Th is property is required so as to decompose the matrix in its singular values.
The plate's stiffness decomposition is straightforward, as shown in Eq. 19.Noticing the necessity of P (1) being positive-definite, which results in N (1) positive-definite, and the restriction to the maximum absolute value of [N (1) ] -1 N (2) eigenvalues, 3 conditions are essential to the recursive assignment convergence: Th e fi rst and the second are directly related to the necessity of all eigenvalues of [N (1) ] -1 N (2) to be lesser than the unity while the last comes from the requirement of [N (1) ] -1 having an inverse.Th e fi rst requirement is readily satisfi ed due to the constitutive tensor positivity, despite the plate's inhomogeneity.In order to ensure the second and third inequalities, an isotropic material is simulated with laminae elastic constants.Considering a slightly modifi ed version of the projectors (Browaeys and Chevrot 2004), M, one can write it as: where each term of the expanded solution is recursively obtained (proposition 3) as: so as to apply them to the reduced constitutive tensor in a vectorized form, described by: as: where POS denotes the spatial position of assessment.
The material properties used in all numerical results are: E ( where (H) and (L) denote higher and lower with respect to the symmetric hierarchy (Forte and Vianello 1996;Chadwick et al. 2001;Ting 2003); it is possible to extract an isotropic part from each lamina.
As aforementioned, the constitutive hierarchy is defi ned in the 3-D elasticity.For the reduced constitutive tensor, the hierarchy is simplifi ed as follows.
The projections of Eq. 24 extract an isotropic tensor of any anisotropic one (monoclinic, orthorhombic and cubic as in Fig. 2).Being isotropic, the last condition (Eq.23) is fulfilled.The second condition must be always tested.However, given the large difference between the 2 main directions in laminated plates, this requirement is almost readily satisfied.

NUMERIC RESULTS
In order to determine the accuracy and convergence of the proposed methodology and, therefore, if the ADM could be used to simulate laminated plates, 3 cases were compared to analytic/numeric solutions available in the literature.Th e results are made non-dimensional so as to ease the results comparison.Th ey are described by: Th e recursive methodology is applied initially for a plate with all the fi bers aligned in the same direction with a uniform loading and boundary conditions as SSSS (Case 1), where the fi rst letter concerns the plate's south edge, turning to righthand direction (Liew and Wang 1993;Liew 1992;Liew and Lam 1991).Th e results obtained are presented in Table 1.It is possible to notice the good agreement of the solution via ADM and the one from (Reddy 1997) using the CLPT.
It is also possible to analyze a second case (Case 2) having 3 layers stacked with with the same boundary conditions as the previous case.The stacking is symmetric what results in decoupling between transverse and membrane displacement.Once again the results correspond to each other, which indicates the validity of the method used.
Figures 3a and 3b show the convergence curves for the proposed methodology, noticing that the vertical axis is non-dimensional and is set as the ratio between the i-th iteration and the converged value.The stopping criterion was defined as the difference between 2 subsequent steps, in an integral fashion (L 2 distance between 2 subsequent solutions).The tolerance for all numeric results was set as 10 -9 .Despite appearing small, since the method is not computationally expensive this criterion is used to assure good results.For the first case the convergence is reached after 366 iterations.The points are plotted every 5 iterations in both cases.Figure 3b shows that the method converges on the iteration 257.As one can note, a good convergence was reached with approximately 75 iterations, in both cases.Another case of study was a plate with 9 layers so as to analyze the methodology for a larger number of laminae.Th e boundary conditions were set as SSSS with uniform loading.Stresses diagrams as well as the convergence are plotted in Fig. 3. Th e stresses, in this case, have dimension in MPa.Th ere are, in Figs.4a and 4b, 2 main diagonals represented in the dashed line.Th is happens because the layer is rotated , which means that every 2 layers are in the same direction.Th e highest stress is, as expected, in the last plate perpendicular to the mid layer (σ 1 and σ 2 ).Convergence is presented in Fig. 4d, which has been reached aft er 146 iterations.
In order to show that the methodology can be applied in non-symmetric stacking, aleatory fi ber orientations were  Lisboa TV, Geiger FP, Marczak RJ Element Model (FEM), according to (Ansys® 2016), which is constructed by linear shell elements.As aforementioned, in these cases, there is coupling between the membrane and bending effects.The results are expressed in Table 2.For all cases, the boundary conditions are CCCC while the loading is uniform and unitary.

Design results:
Door designed with intentional warp provides stress against latches, to prevent rattling.Door can restrain landing gear that is inadvertently deployed yet can be opened if door hydraylics fail.
Five of six stringers cobonded with door skin for maximum weight savings height equal to 2m, 0.9m and 0.9m, respectively, spatially positioned at 45 o from the main part.

•
The curved part of the landing gear door is inserted into the problem as structural forces applied to the main part of the door.In considering these simplifications, the geometry of the component is then reduced to a plate with dimensions 2.2 m × 0.9m × 0.9m, respectively, height × length × thickness.The boundary conditions considered are FFCF as if the top edge of the door were rigid-connected to the plane's fuselage while the curved part is connected to the bottom edge of the plate.Five cases of different stacking and fiber orientation are analyzed, as described in Table 3.These cases are chosen in order to embrace several types of stacking in laminated plates.The forces applied to the plate (Fig. 4) are: • Aerodynamic in the plate, q a , which is uniform and applied in the out-plane direction.The last example, in which the methodology is used, concerns a landing gear door.The idea of this example is to develop some analytical solution for pre-analysis of this component.The model was based on the main landing gear door of an Airbus A350-900 (AIRBUS® 2015; Fig. 5a).
The geometry is simplified (Fig. 4b) by the following assumptions:

•
The reinforcements are ignored.

•
The curved part is viewed as a plane trapezoidal geometry, with the larger base, smaller base and the (a) (b) Table 5. Transverse displacement using the presented methodology and FEM.

•
Dragging force, q d , which depends on the door's lateral area and the drag coefficient and it is applied in the in-plane direction.
• Consideration of the curved part equivalent forces: in-plane, out-plane and applied bending moments are considered and all depend on the plate's area.In total, there are different loadings: 6 come from the curved part of the landing gear door (a resultant of the dragging loading, 2 bending and 1 twisting moments as well as the 2 resultant forces applied to the plate) and 2 forces are equivalent to the out-plane loading in that part.The analysis was developed both utilizing the presented methodology and FEM.Two points were assessed (K 1 and K 2 ) and they correspond to the vertices of the rectangle, in the free edge opposite to the clamped one (same edge which the curved part is connected).For the numeric results, and are considered 1 N/m 2 and 0.01 N/m 2 , respectively.These values were taken to combine both effects (aerodynamic and dragging) with same order.
Tables 4 and 5 show the results for the membrane and transverse displacement respectively.Slightly differences can be observed due to the small differences in the way the forces are applied in each model.One example is the torsion resultant to the dragging loading in the curved part.In the presented methodology, there is no degree-offreedom which interpolates such behavior.Consequently, this effect was considered by an equivalent loading (in-plane forces).On the other hand, in the FEM model, this degreeof-freedom could be directly considered.Despite the small differences, one can observe that the methodology can solve complex laminated problems, both in unsymmetrical stacking and complicated types of stacking and combination   of loadings.Moreover, the solutions obtained have an analytical base, giving the possibility of analytic derivatives, which can definitely ease optimization processes as well as delamination analysis.with the stacking modification.In the analyzed plates, u 0 2 was the mostly sensible to the stacking modification.

CONCLUSIONS
A new methodology to determine the behavior of laminated plates was proposed.The equations were derived according to the pb-2 Rayleigh-Ritz Method and the Adomian Decomposition Method in order to obtain the stresses and displacements of such structures.The results were then compared, for 2 cases, with analytic solutions available in the literature, showing good agreement.Another case was tested and 2 main diagonals appeared in the stress results, since there are 2 main directions on the laminated plate.The methodology was applied to unsymmetrical stacking and excellent correspondence with both analytic and numeric solutions was found.Furthermore, a case of study was presented and the landing gear door of A350-900 was

REFERENCES
Adomian G (1994) Solving Frontier problems of physics: the decomposition method.Dordrecht: Kluwer Academic Publishers.analyzed.Although the model was simplified in a geometric sense, the results herein obtained could have served as guide to designers so as to give information on the behavior of the structure facing such loading conditions.The analytical procedures as well as the numeric results have shown that the recursive procedure is a powerful tool that can be used to solve many types of engineering problems in a very straightforward fashion.Moreover, together with pb-2 Rayleigh-Ritz Method, one can derive analytic derivatives of any problem variable, which can ease optimization processes.

Figure 2 .
Figure 2. Constitutive hierarchy for the reduced constitutive tensor.

Figure 3 .Figure 4 .
Figure 3. Convergence curves for (a) Case 1 (orientation angle and (b) Case 2 (orientat ion angle) with the displacement normalized by the fi nal value.
considered and compared to analytic solutions (Bhaskar and Dhaoya 2009) as well as to numeric results obtained by the Finite

Figure 5 .
Figure 5. (a) Main landing door of an Airbus A350-900 (Sloan and Reque 2014); (b) Mathematical model with the restriction and the aerodynamic and dragging forces.

Table 1 .
Comparison between the CLPT and the ADM for displacements and stresses for 2 different angles of staking.

Table 3 .
Cases analyzed considering different types of staking together with laminae orientation.

Table 2 .
Non-dimensional center transverse displacement in plates with non-symmetric stacking.

Table 4 .
Membrane displacements using the presented methodology and FEM.