Three-Dimensional Bending Analysis of Multi-Layered Orthotropic Plates by Two-Dimensional Numerical Model

The paper presents effective numerical modelling of multi-layered plates with orthotropic properties. The method called the FEM23 is employed to construct the numerical model. The approach enables a full 3D analysis to be performed while using a 2D finite element mesh. The numerical model for a multi-layered plate is constructed by an assembling procedure, where each layer with orthotropic properties is added to the global numerical model. The paper demonstrates that the FEM23 method is very flexible in defining the multilayered plate, where the thickness of each layer as well as its mechanical orthotropic properties can be defined independently. Several examples of three-layered or nine-layered plates are analyzed in this paper. The results obtained by the FEM23 method coincide with the ones taken from the published papers or calculated with the standard 3D FEM approach. The orthotropic version of the FEM23 can be quite easily applied for other kinds of problems including thermo-mechanics, free vibrations, buckling analysis, or delamination.


Introduction
Multi-layer plates and shells with orthotropic properties in the layers are types of structures widely used in modern aerospace, automotive, and construction industries, due to their lightness, elasticity, durability, and excellent mechanical properties in regard to their thickness and weight. The significant increase in the use of such composite materials has been observed over the past few decades, especially in the aviation industry, where constructive components must show high strength while maintaining a slim shape. The appropriate design of structures made of orthotropic materials requires reliable computational modelling techniques.
The paper proposes a modified method called the FEM23 [1,2], whose main advantage is modelling full three-dimensional (3D) multi-layered structures while using a 2D finite element mesh for calculations. The modification of the FEM23 involves its extension to layers exhibiting orthotropic properties. The paper proves that the FEM23 method can be used effectively to analyze any multi-layered structure combining thick and thin orthotropic layers. The FEM23 is a fully 3D method, similar to the Proper Generalized Decomposition applied in [3,4], which allows for solving a fully 3D model, but with 2D characteristic computational complexity. This approach uses the spatial decomposition of a displacement field combining the in-plane 2D solution with the transverse 1D solution. As a result, another 1D mesh is generated alongside the planar 2D mesh. The method subsequently developed in [5], applies the piecewise fourth-order Lagrange polynomial along the plate's thickness. In previous papers [6,7], the authors have already utilized the FEM23 to model the behavior of multi-layered plates and shells with isotropic layers subjected to mechanical or thermo-mechanical loads, devoting special attention to laminated glass structures-the topic of many numerical or experimental studies [8][9][10][11].
The body of research regarding the analysis of orthotropic multi-layered plates is limited here only to papers published in recent years, to show up-to-date achievements in the area. Rajaneesh et al. in [12] applied a variationally consistent new first-order shear deformation theory to three-layered orthotropic plates using an equivalent-single-layer theory in the finite element. A numerical analysis of elastic waves in the multilayered orthotropic plates was presented in the paper by Liu et al. [13], where the extended Legendre polynomials combined with the anisotropic couple-stress theory were applied to investigate the reflection and transmission of the waves. A layerwise generalized theory of layered orthotropic plates was proposed by Ugrimov and Shupicov [14]. The theory is based on a power series expansion of the displacement vector component in each layer for the transverse direction, where the number of terms retained in the power series is arbitrary and chosen according to the problem being considered. Xu et al. [15] proposed analytical solutions for orthotropic thin plates, where the finite integral transform method was introduced for accurate bending analysis. This type of analysis is limited only to thin plates due to the Kirchhoff-Love plate theory applied to derive the employed method. The vibration response of orthotropic composite plates was analysed by Zhang et al. [16] with a new approach regarding the analytical solution, while 3D semi-analytical vibrations of thick orthotropic plates were investigated by Cui et al. [17]. In this research, the modified Fourier series were applied to obtain admissible displacement functions while boundary conditions were generated by arranging sets of linear springs at the edges. Belyaev et al. [18] proposed a two-dimensional model to analyze deformations of multi-layered orthotropic plates, where a multi-layered plate was replaced by an equivalent plate composed of a monoclinic material with piecewise elastic modules, and an asymptotic expansion was applied to obtain the solution. The buckling analysis of orthotropic laminated plates with asymmetrical boundary conditions was carried out by Schreiber et al. [19]. The framework of Reddy's third-order shear deformation theory was applied for a closed-form solution, and Lévy-type solution was derived for the plate buckling problem.
The multi-layered and orthotropic plates are the subject of both numerical and experimental research. The combined experimental and analytical study of multi-layered laminated glass subjected to low-velocity impacts was presented by Wang et al. [20], while the experimental determination of material properties in an elliptical orthotropic plate was proposed by Marchetti et al. [21]. The method validated subsequently in a honeycomb sandwich panel combined the equivalent thin plate theory with the wave fitting approach. An innovative transparent polymeric laminate consisting of two poly(methyl methacrylate) plies with a thermoplastic polyurethane interlayer was experimentally and numerically investigated by Rühl et al. [22]. The laminate was subjected to low-velocity impact loadings using clamped three-point bending and dart impact tests at different temperatures and velocities. Multi-layered metallic plates were the subject of a number of experiments conducted by Ziya-Shamami et al. [23]. The paper analyzed 15 various testing groups, presenting the data obtained in over a hundred experiments. The plates were subjected to repeated uniform impulsive loadings to study the effects of different charge masses, layering configurations, layering arrangements, layer thicknesses, and stand-off distances on the central deflection of single-, double-, and triple-layered mixed plates made of aluminum alloy and mild steel materials.
In this paper, the FEM23 method is extended in such a way that a multi-layered plate with orthotropic properties can be analyzed. The orthotropy in the numerical model is obtained by the stress-strain constitutive relation defined in the coordinates associated with the orthotropy directions, subsequently transformed to global coordinates and applied to the mathematical model shown in detail in Section 2. A specific form of displacement approximation used for the purpose of this study and the final numerical model are presented in Section 3. The FEM23 approach for orthotropic multi-layered plates is verified in three examples analyzed in Section 4. The paper ends with some conclusions.

Variational Formulation of the Problem
The starting point of the presented FEM23 method for the multi-layered orthotropic plates is a single-layered orthotropic plate model subsequently extended to a multi-layered one. Two Cartesian coordinates, orthotropic (ξ, η, θ) and global (x, y, z), are distinguished in the orthotropic plate. The θ and z directions are the same and oriented transversely, while the (ξ, η) pair is associated with in-plane orthotropy in a layer, and the (x, y) pair is correlated with the in-plane geometry of the plate. Although the (x, y) pair is common for the whole plate, the (ξ, η) directions are oriented differently in each layer. The orthotropic coordinates are used to derive constitutive relations in the layers, while the global ones are used to describe the model for the whole multi-layered plate.
The constitutive (Hooke's) relation for orthotropic materials is expressed with the orthotropic coordinates as follows:σ whereũ is the displacement vector, whileσ andε are the stress and Cauchy's strain tensors, all defined with the use of the orthotropic coordinates.D is a symmetric constitutive tensor of the fourth order in the orthotropic coordinates, demonstrating the following non-zero values [24]:D E i represents the Young moduli of elasticity, G ij is shear moduli, while ν ij are Poisson ratios, all in orthotropic directions. Using the symmetric properties of theD tensor, the following relation is derived: The stress and strain tensors from Equation (1) satisfy the equilibrium equation, which can be written in the variational form for the whole considered domain Ω: which has to be satisfied for allṽ whereṽ = 0 on Γ u . (3) refers to the tensor inner product. For the sake of further analysis, it can be shown that, due to the symmetric properties of theD tensor, the inner product can be expressed with full displacement gradients:

The second integral in Equation
where∇ is the gradient operator defined with the orthotropic coordinates.
Assuming that Q is the transformation matrix from the orthotropic to global coordinates, the relation between the displacement gradients expressed with the orthotropic and global coordinates is as follows:∇ũ where u is the displacement vector in the global coordinates, and ∇ is the gradient operator in the global coordinates. The displacement gradient is defined in this paper as When the relation in Equation (5) is subsequently applied to the inner product in Equation (4), the following relations are derived: where D is the constitutive tensor written in the global coordinates with its components expressed as In this paper, the orthotropic and global plate coordinates are always rotated around the z axis. Consequently, the transformation matrix Q takes the following form: where c = cos(α), s = sin(α) and α constitute a rotation angle between the first axes of the global and orthotropic coordinate systems around the z axis. The relation between those two coordinates is illustrated in Figure 1. Taking into account the form of the transformation matrix shown in Equation (9) the Hook tensor defined with the global coordinates in Equation (8), takes the following form: The integral Equation (3) expressed with the orthotropic coordinates can now be rewritten, so that all the components are described with the global coordinates Using the method described in the authors' previous papers, for instance, [1,2,6], the inner product in Equation (7) is then expressed with the derivatives in the transverse and in-plane directions, consequently changing into where∇ u = ∂u x ∂x ∂u y ∂x ∂u z ∂x ∂u x ∂y ∂u y ∂y ∂u z ∂y (12) The D i tensors constitute adequate parts of the Hooke's D tensor, defined as follows: The relations in Equation (13) provide definitions of D i tensors in the relation D one. However, the explicit definitions of the non-zero components of the tensors are as follows Equation (11) is now applied to Equation (10). Moreover, the volume integrals are separated into the in-plane and transverse ones: where Γ m is the mid-surface of the plate. The 2D integration over the Γ m surface can be performed using a 2D finite element mesh, while the transverse one-dimensional integration can be completed using standard 1D numerical integration schemes where z i are the integration points in the transverse direction. Equation (19) is written for a single-layered orthotropic plate. Assuming that the plate comprises M orthotropic layers, each with its own orthogonal orientation, the variational equation for the whole plate is as follows: where D j k are material tensors defined for the jth layer of the plate and z j i are the transverse integration points for the jth layer.

Approximations
The details concerning spatial approximation have already been provided by Jaśkowiec et al. [1]. This paper presents only an outline of the approximation scheme which is derived first for a single layer and then extended to a multi-layered case.
The spatial approximation of the displacement in the jth layer is a combination of the in-plane approximation and the transverse one. Using the 1D transverse approximation along the plate's thickness, the displacement vector u (j) (x, y, z) is expressed in the following way: where S j is the number of approximation surfaces specified along the thickness of the jth layer dependent on the transverse approximation order, N i is the i-th Lagrange polynomial of the appropriate S j − 1 order, while u (j)i (x, y) is the displacement vector on the i-th approximation surface of the jth layer. The same kind of the displacement decomposition has been previously proposed in [25,26] for hierarchical shell models. The approximation order depends on the layer's thickness. The first order of the transverse approximation may be applied when the layer is thin. However, higher orders of the transverse approximation have to be utilized when the thickness of the layer is significant in relation to the outer dimensions of the plate. On the other hand, every vector u i (x, y) is, as in previous papers such as [1,[25][26][27][28], approximated using the same approximation field constructed on a single 2D mesh: where Φ(x, y) is the approximation matrix defined on a 2D in-plane mesh, andǔ (j)i is the vector of nodal displacements in the i-th approximation surface of the j layer. Using the Equations (21) and (22), the displacement vector in the entire jth layer is approximated as follows: where the global vector of degrees of freedomǔ (j) comprises the vectors associated with the subsequent surfaces of the j layerǔ Two consecutive, adjacent layers share a common surface and, consequently, the degrees of freedom associated with it. The common degrees of freedom of the two adjacent layers satisfy the relationǔ The variational Equation (20) requires the in-plane gradient and the transverse derivative of the displacement field. Such derivatives for the displacement approximation in Equation (23) are expressed as follows: After substituting the approximation shown in Equation (23) and derivatives from Equation (26) to Equation (20), the following discrete system of the equation is obtained: where particular matrices and the right-hand-side vector are generated with the assembling procedure across the plate's layerš where A (j) is the assembling operator for the jth layer. The definitions of other components in Equation (28) are as follows: where z σ is the surface on the top layer, where tractiont is applied.

Examples
This section presents three examples demonstrating the accuracy and efficiency of the FEM23 method for orthotropic multi-layered plates. In the examples, the results yielded by the FEM23 approach are compared with the ones from other published study. The first example, Section 4.1, focuses on a simply supported three-layered square plate whose layers demonstrate the same mechanical properties, but their orthotropic directions are rotated. The results yielded with the layerwise theory and 3D FEM are compared. Additionally, various types of simply supported boundary conditions and the case study of plates with varying thickness are analyzed. The second example, Section 4.2, investigates a plate with layers of varying values of thickness. In the first case, the plate consists of three layers, and in the second one, of nine layers constructed with three different materials. This example shows how the FEM23 method copes with multi-layered structures made of different orthotropic materials with various thickness. The FEM23 results are verified with the zigzag theory and 3D FEM. In the third example, the plate under consideration is a three-layered one, in which the orthotropic directions are rotated in relation to each other. In this example, the FEM23 results are again compared with the ones obtained using the zigzag theory and the standard 3D FEM.

Example 1
In this example, a simply supported three-layered square plate loaded by the pressure on the upper outer surface is analyzed. The side length of the plate is a, while the thickness of the plate is h with h 3 thickness for each layer. The geometric and material properties of the plate are shown in Figure 2. In this example, the results obtained by the FEM23 are compared with the ones presented in [29,30], where layerwise plate theories were applied.
In the FEM23 a full 3D analysis is performed, using standard 3D finite elements in the ABAQUS application. The fact that the plate is simply supported can be easily applied in the case of plate theory. However, in a 3D analysis, such boundary conditions are ambiguous. In a simply supported 3D plate, the outer points at the bottom of the plate may be blocked and in another approach-the outer points at the mid-surface of the plate. Moreover, the plate can be hard or simply supported. Figure 3 presents graphically different versions of the simply supported boundary conditions in a 3D analysis. In this example, the results obtained by the FEM23 method are compared with the ones taken from the papers by Vuksanović et al. [30] and Carrera and Ciuffreda [29], where multilayer plate theories were applied, as well as with the results yielded with the standard 3D FEM method using the Abaqus package. In the FEM approach, each layer of the plate is discretized by a single layer of hexahedral elements. Table 1 shows all the results in the form of a dimensionless deflectionū z and stress componentσ x , both evaluated at two central points of the plate:ū The results yielded with the FEM23 and FEM presented in Table 1 coincide for all sorts of applied supports and for every plate thickness, which proves that the FEM23 provides an accurate 3D analysis. The results are also in good correlation with the ones obtained with plate theories, but only when boundary conditions are applied to the midsurface nodes in the FEM23 and a thin plate is considered. When the relation a/h increases, both types of results become inconsistent, which is natural for plate models based on the Kirchhoff-Love theory efficient only for thin plates.
The results are illustrated with the plots of the dimensionless deflections and the stressσ x for the plates with a/h = 4 and a/h = 100 along the plate's thickness shown in the middle of the plate (see Figures 4 and 5). It can be noted that the application of the boundary conditions plays a significant part in the case of the 3D model. Different results are yielded when the support is applied at the bottom edge or in the middle of the plate's thickness. The bottom edge support case is very sensitive to both hard and soft application of the boundary condition. Such sensitivity is not observed when the middle edge support is used. The following Figures show the comparison of the obtained results with the ones based on the plate theory taken from [29]. In the case of thin plate, the results of the [mh] type boundary conditions almost overlap (see Figure 5). However, when a thick plate is considered (see Figure 4), these two types of results differ, which is typical for the plate theory approaches.   The FEM23 method uses full 3D modelling with the calculations performed on a 2D domain. In this example, the FEM23 results are compared with the ones obtained by the standard 3D FEM method with Abaqus application. As can be seen in Table 1, both approaches yield very similar results for all types of boundary conditions. Both methods are also compared in Figure 6, where the presented 3D deflection and stress maps are very similar to each other. It can be observed that the relative differences of the maximum deflections and stresses are around 7 · 10 −3 and 2 · 10 −1 , respectively. Figure 6. Example 1: The maps of deflection u z (a,b) and the stress component σ xx (c,d) presented on the half of the deformed plate obtained by 3D FEM and FEM23 for the bs boundary type.

Example 2
The following example, excerpted from [31], applying the Hellinger-Reissner refined zigzag theory (HR-RZT), analyzes a square multi-layered plate of a = 2 m total thickness h =0.1 m], with a constant load q 0 = 0.1 N/mm 2 . The plate is constructed of three types of materials called A, B, and C, whose properties are provided in Table 2. The material A is an orthotropic carbon-fiber reinforced plastic, the material B represents an orthotropic titanium honeycomb structure, and the material C corresponds to an isotropic foam made of polyvinyl chloride. The analysis involves, as in considered paper [31], two types of multi-layered plates: the ones consisting of three and nine layers. The structures of the two plates are presented in Table 3, and the nine-layered plate is shown in Figure 7. It is assumed that the angle of the orthotropic direction α is equal to zero for each layer. The calculations performed in the FEM23 method utilize a quadrilateral mesh structure with 42 × 42 nine-node elements. In this example, the coordinates' origin is located in the middle of the plate. The boundary conditions are as follows: u z = 0 for the points (±a/2, y, 0) and (x, ±a/2, 0); u x = u y = 0 for the point (−a/2, −a/2, 0); u y = 0 for the point (a/2, −a/2, 0) (31)  The results presented in the form of diagrams of deflection and selected quantities along the plate's thickness are shown in Figures 8 and 9 for the L1 and L2 plates, respectively. The diagrams provide a comparison of the FEM23 results with the ones yielded by the HR-RZT and full 3D FEM approaches provided in the paper [31]. In the case of the threelayered plate L1, all three types of results coincide. When the case of the nine-layered L2 plate is considered, the results of the 3D FEM and FEM23 are almost the same and differ only slightly when compared with the HR-RZT results. All the above observations confirm the effectiveness of the FEM23 approach for a numerical analysis of multi-layered phototropic plates.

Example 3
The following example investigates a three-layered 50 × 50 mm composite square plate of 1 mm thickness subjected to a uniform pressure of q z = 1 N/mm 2 . Each layer is made of the same orthotropic material with the 0 • /90 • /0 • orientations (for 0 • , axis 1 corresponds to the global x axis). The material properties are presented in Table 4. In the given example, the FEM23 calculations are performed with the use of a 2D 70 × 70 structural mesh consisting of nine-node quadrilateral elements. The coordinates' origin is located in the middle of the plate, and the boundary conditions remain the same as presented in Equation (31). The FEM23 results are compared with the ones yielded by the HR-RZT [31] and 3D FEM [32] methods. Figure 10 shows the diagrams of two shear stress components σ xz and σ yz along the plate's thickness at two selected points of the plate generated with the FEM23, 3D FEM, and HR-RZT, respectively. In this case, the FEM23 results almost coincide with the results of the two other approaches, with only a slight difference observed between them. (d) Figure 9. Example 2. Results for plate L2: deflectionū z (x, 0, 0) (a), displacementū x (3/7l, 0, z) (b), normal stress σ x (3/7l, 0, z) (c), shear stress σ xz (3/7l, 0, z) (d).

Conclusions
This paper has presented the FEM23 method adapted for a numerical analysis of orthotropic multi-layered plates. The novelty of the paper lies in using the FEM23 approach to derive a numerical model for a single orthotropic layer and extend its application to multi-layered plates. The FEM23 approach utilizes a 2D finite element mesh to perform a full 3D analysis of multi-layered structures. The method, in contrast to the standard 3D FEM, can efficiently deal with structures comprising a combination of thick and thin layers made from different materials. The numerical model generated with the FEM23 is based on the 2D-1D decomposition in the physical model, spatial integration, and approximation. The physical properties of the layers are defined for each one independently. The layers' orthotropic properties are introduced with the fourth order Hooke's tensor for orthotropic materials expressed with the local, orthotropic coordinates. Due to the symmetry of Hooke's tensor, the strain tensor is substituted by the displacement gradient, and the transformation to the global coordinates of the tensor is derived. Then, the numerical models of the individual layers are combined to assemble the global model based on a 2D finite element mesh. After postprocessing, full 3D results are obtained.
The paper is illustrated with examples analyzing multi-layered plates subjected to external loads. In each case, the FEM23 results are compared with the ones taken from several published papers, based on plate theories, or obtained with the standard 3D FEM method. The FEM23 and the 3D FEM results generally coincide. In the case of thin plates, the FEM23 results are comparable with the ones yielded by plate theories. However, when thicker plates are considered, the results begin to diverge, which is common for plate theories. Modelling multi-layered plates is relatively easy in the FEM23 method, where the layers, thick or thin, demonstrating homogenous or orthotropic properties, are combined to assemble the final numerical model.
The full potential of the FEM23 method is yet to be discovered, especially with its Matlab-based software still evolving to be more efficient and able to model a spectrum of mechanical issues. In the near future, the applications of the presented method in the context of multi-layered anisotropic structures may include: buckling and post-buckling analysis, dynamics, thermo-mechanical analysis, or delamination processes using large deformation descriptions.