Torsional and Transversal Stiffness of Orthotropic Sandwich Panels

In the present work, an analytical equation describing the plate torsion test taking into account the transverse shear stiffness in sandwich plates is derived and numerically validated. Transverse shear becomes an important component if the analyzed plate or shell is thick with respect to the in-plane dimensions and/or its core has significantly lower stiffness than the outer faces. The popular example of such a sandwich plate is a corrugated cardboard, widely used in the packaging industry. The flat layers of a corrugated board are usually made of thicker (stronger) material than that used for the corrugated layer, the role of which is rather to keep the outer layers at a certain distance, to ensure high bending stiffness of the plate. However, the soft core of such a plate usually has a low transverse shear stiffness, which is often not considered in the plate analysis. Such simplification may lead to significant calculation errors. The paper presents the generalization of the Reissner’s analytical formula, which describes the torsional stiffness of the plate sample including two transverse shear stiffnesses. The paper also presents the implementation of the numerical model of the plate torsion test including the transverse shear stiffnesses. Both approaches are compared with each other on a wide range of material parameters and different aspect ratios of the specimen. It has been proved that both analytical and numerical formulations lead to an identical result. Finally, the performance of presented formulations is compared with other numerical models using commercial implementation of various Reissner–Mindlin shell elements and other analytical formulas from the literature. The comparison shows good agreement of presented theory and numerical implementation with other existing approaches.


Introduction
Composite materials play an important role in many practical engineering designs. They consist of at least two materials, which are called fiber and matrix material, or are composed of more than one layer. A special class of composite layered materials is sandwich plates; they consist of a soft core and two outer sheets. The internal one supports the faces and provides the overall flexural stiffness of a structure. In composite laminates or sandwich plates and shells, an important role is played by flexural and torsional stiffness. A typical example of such a sandwich plate is the 3-layer corrugated cardboard, in which the core is corrugated and therefore light. The role of the core called "flute" in jargon is to provide shear stiffness while the outer sheets called "liners" ensure the flexural stiffness. Transversal shear effect is especially important in thick plates. resulting from the CLPT. In the analysis of boxes made from corrugated cardboard, quite a different approach is strength prediction based on empirical formulas [23], including the very popular McKee formula [24], which is very important from a practical point of view.
In this paper, a new generalized form of the Reissner analytical formula describing the torsion of the sandwich plate is proposed. The torsional stiffness and both transverse shear stiffness are taken into account in these considerations. The presented formula is compared with the numerical models and other analytical approaches existing in the literature.

Reissner-Mindlin Plate-Governing Equations
A plate of constant thickness t was under consideration herein. Its middle surface lies in the plane x-y, while the "thickness direction" is related to the z-axis. It was assumed that cross sections remain straight during deformation, but they do not necessarily remain normal to the middle surface.
The assumption is typical for shear deformable plates.
The equilibrium equations are expressed in terms of the bending moments m xx and m yy , the twisting moments m xy = m yx and the shear forces q x and q y . An infinitesimal plate element dA = dx · dy subjected to an external transversal load q is depicted in Figure 1, which shows also the stress resultants. The plate equilibrium equations consist of the transversal equilibrium: q x,x + q y,y = q (1) and the rotational equilibrium in x and y: m yy,y + m xy,x − q y = 0, where comma represents the partial derivative of the variable with respect to Cartesian coordinate direction that appears after the comma.
Materials 2020, 13, x FOR PEER REVIEW 3 of 17 cardboard, quite a different approach is strength prediction based on empirical formulas [23], including the very popular McKee formula [24], which is very important from a practical point of view.
In this paper, a new generalized form of the Reissner analytical formula describing the torsion of the sandwich plate is proposed. The torsional stiffness and both transverse shear stiffness are taken into account in these considerations. The presented formula is compared with the numerical models and other analytical approaches existing in the literature.

Reissner-Mindlin Plate-Governing Equations
A plate of constant thickness was under consideration herein. Its middle surface lies in the plane -, while the "thickness direction" is related to the -axis. It was assumed that cross sections remain straight during deformation, but they do not necessarily remain normal to the middle surface.
The assumption is typical for shear deformable plates.
The equilibrium equations are expressed in terms of the bending moments and , the twisting moments = and the shear forces and . An infinitesimal plate element = ⋅ subjected to an external transversal load is depicted in Figure 1, which shows also the stress resultants. The plate equilibrium equations consist of the transversal equilibrium: and the rotational equilibrium in and : , where comma represents the partial derivative of the variable with respect to Cartesian coordinate direction that appears after the comma. In the present study, is the transversal displacement and and denote the rotations of the plate cross section. Orientation of the rotations is shown in Figure 1. and denote the shear strains and have the same definitions of the orientations as for the rotations. The kinematic equations take the form: In the present study, w is the transversal displacement and θ x and θ y denote the rotations of the plate cross section. Orientation of the rotations is shown in Figure 1. γ x and γ y denote the shear strains Materials 2020, 13, 5016 4 of 18 and have the same definitions of the orientations as for the rotations. The kinematic equations take the form: The curvatures are defined as: A linear elastic orthotropic material was assumed here. The plate bending stiffness D and the shear stiffness D s are expressed in terms of the material parameters: E 11 , E 22 , ν 12 , ν 21 , G 12 , G 13 , G 23 and the plate thickness t as: where E 11 is the effective stiffness modulus in the x direction, E 22 is the effective stiffness modulus in the y direction, ν 12 is the effective Poisson's ratio in the 1-2 (xy) plane, v 21 = v 12 E 22 /E 11 , G 12 is the effective shear modulus in 1-2 (xy) plane, G 13 is the effective transverse shear modulus in the 1-3 (xz) plane, G 23 is the effective transverse shear modulus in the 2-3 (yz) plane, D 11 and D 22 are bending stiffnesses in two orthogonal directions, D 33 is the twisting stiffness, A 44 is the transverse shear stiffness in the 1-3 (xz) plane and A 55 is the transverse shear stiffness in the 2-3 (yz) plane. The shear correction parameter, α, describes the non-constant distribution of shear stresses through the plate thickness. For a homogenous orthotropic plate with constant thickness, the standard value of α = 5/6 was adopted [17,18]. The relation between bending moments and the curvatures can be expressed as: and the shear forces are related to the shear strains by: where D 11 , D 12 = D 21 , D 22 , D 33 , A 44 and A 55 are defined in Equations (7) and (8).

Torsion of Orthotropic Plates with Transversal Shearing
The foregoing constitutive equations incorporate a limiting-type orthotropy assumption due to sample size considered (a b, see Figure 2a), which causes m yy and q y to be reactive quantities, which enables a simple conversion from plate theory to a beam theory, in a physically justifiable way. This assumption leads to [25,26]: This formula represents the torsional stiffness of the sample including transversal shear stiffness along one of the direction, see Figure 2a, where T is the torque and β is the angle of rotation; a and b are dimensions of the sample. Further research showed that despite the fact that AKK [3] has good agreement with simple, but commonly used numerical solutions, it can be replaced with an exact solution AA, which is represented by Equation (13) or the numerical approach presented in one of the following subsections. The accuracies of the selected methods are checked on various examples in Section 3.

Reissner-Mindlin Composite Laminated Plate
Assuming the coordinate system, the displacement field for composite laminated plates reads: Substituting the torque with a pair of forces, see Figure 2b, and modifying the above equation to incorporate both transversal shear stiffnesses, we received the following analytical form: From hereon, this analytical approach will be abbreviated as AA. Here, R is the reaction force, while w is the vertical displacement. In the inverse form, the expression reads: Determining the transversal shear stiffness of the sample from a static torsion plate test has been recently studied by the authors [3]. The study extended the work of Aviles et al. [2], where the additive form of w/R was used, in which nonphysical constant c was used. In the paper [3], constant c was replaced by the term k 1 k 2 , where k 1 and k 2 are described by material properties and sample geometry (this approach is abbreviated as AKK (approach k 1 k 2 ) in forthcoming sections), so the final form of the expression reads: where Further research showed that despite the fact that AKK [3] has good agreement with simple, but commonly used numerical solutions, it can be replaced with an exact solution AA, which is represented by Equation (13) or the numerical approach presented in one of the following subsections. The accuracies of the selected methods are checked on various examples in Section 3.

Reissner-Mindlin Composite Laminated Plate
Assuming the xyz coordinate system, the displacement field for composite laminated plates reads: where u, v and w are three unknown mid-surface displacements of the plate, while θ x and θ y are two rotations of the normal on the plane xz and yz from Equation (19): noting that γ x = −φ x and γ y = −φ y . In both the xz and yz vertical planes, the normal rotation was obtained as the sum of two rotations: (i) the corresponding slope of the middle plane of the plate and (ii) the additional rotation φ, which results from the lack of orthogonality of the normal to the middle plane after deformation, see Figure 3. Therefore, the rotations θ x and θ y cannot be calculated from the deflection only and become independent variables. This is the fundamental difference between Reissner-Mindlin and Kirchhoff-Love plate theories.
Materials 2020, 13, x FOR PEER REVIEW 6 of 17 where , and are three unknown mid-surface displacements of the plate, while and are two rotations of the normal on the plane and from Equation (19): noting that = − and = − .
In both the and vertical planes, the normal rotation was obtained as the sum of two rotations: (i) the corresponding slope of the middle plane of the plate and (ii) the additional rotation , which results from the lack of orthogonality of the normal to the middle plane after deformation, see Figure 3. Therefore, the rotations and cannot be calculated from the deflection only and become independent variables. This is the fundamental difference between Reissner-Mindlin and Kirchhoff-Love plate theories. Relationships between strains (membrane-, bending-and shear-) and displacements read: The stress-strain relations in local coordinates are given by: Relationships between strains (membrane-ε, bending-κ and shear-γ) and displacements read: Materials 2020, 13, 5016 The stress-strain relations in local coordinates are given by: or in the compact form The strain energy, U, reads where and A, B, D and D s are stiffnesses of the plates given by: For an n-layer laminate with homogeneous orthotropic material within each of the k-th layers ( Figure 4), Equations (32) and (33) can be rewritten where C k and C sk . are the in-plane constitutive matrix, defined in Equations (25) and (26), for the k-th layer.  The work, , done by the in plane and transverse load is given by: The energy functional, Π, of the plate was finally obtained as follows:

FEM Formulation of the Laminate Plate Element
For decades a finite element (FE) analysis has been a popular method for modeling advanced engineering problems. Due to its popularity and universality, the developers of FE software more and more often provide their users new functions, extending its capabilities. Following this trend, the software often allows one to include users' material or element subroutines tailored for particular needs. Therefore, the FEM approach in comparison to the analytical approach is easier to be applied in the modern engineering tools. In the laminated plate element, the field variables, , were approximated according to the associated node values as follows: where represents number of nodes in the element; ( , ) is a shape function associated with node ; = is the displacement vector of the node degrees of freedom. The membrane, bending and shear strains associated to the displacement in Equation (40) can be therefore obtained as follows:  The matrix B becomes 0 if the laminate consists of homogeneous material, or the properties of a material are symmetrical with respect to the middle plane (z = 0). The membrane and bending effects were then uncoupled and the neutral plane coincided with the plane xy. This means that the bending moments did not cause any membrane strains and the normal forces did not cause any curvature.
The work, V, done by the in plane and transverse load is given by: The energy functional, Π, of the plate was finally obtained as follows:

FEM Formulation of the Laminate Plate Element
For decades a finite element (FE) analysis has been a popular method for modeling advanced engineering problems. Due to its popularity and universality, the developers of FE software more and more often provide their users new functions, extending its capabilities. Following this trend, the software often allows one to include users' material or element subroutines tailored for particular needs. Therefore, the FEM approach in comparison to the analytical approach is easier to be applied in the modern engineering tools. In the laminated plate element, the field variables, d, were approximated according to the associated node values as follows: where N n represents number of nodes in the element; N j (x, y) is a shape function associated with node j; is the displacement vector of the node degrees of freedom. The membrane, bending and shear strains associated to the displacement in Equation (40) can be therefore obtained as follows: Materials 2020, 13, 5016 Substituting Equations (40) and (42)-(44) into Equation (39) leads to: By using Lagrange's equations for the energy expression in Equation (48), the characteristic equation of the system was obtained as follows: where It is noted that the shear locking phenomenon can appear in Equations (50) and (51) as the plate thickness decreased. To overcome this adverse the Reissner-Mindlin plate quadrilateral element with assumed transverse shear strain fields was adopted here. We used the plate element with 4 nodes with linear shear field, which initially was developed by Bathe and Dvorkin [27,28]. Its formulation bases on auxiliary transverse shear modes proposed by MacNeal [29,30] and Hughes et al. [31]. Later Donea and Lamain [32] and Onate et al. [33,34] derived the element using assumed strain concepts.
A standard 4-noded Q4 element [33,34] is characterized by a bilinear interpolation of deflections and rotations. The assumed transverse shear strain field is defined here in the natural system ξ, η as The α i parameters can be found by taking the natural transverse shear strains γ ξ at the all four middle-edge points shown in Figure 5, with where β i is the angle between the ξ i and ξ axis. This leads to: where The strains γ i ξ are related to γ i ξ , γ i η by  The relationship between the Cartesian transverse shear strains at the middle-edge points and the natural transverse shear strains is: The Cartesian transverse shear strains at the four sampling points are then related to the nodal displacements by: The substitute transverse shear strain matrix was obtained by = . (61) Here, the element is described as QLLL (for quadrilateral, bilinear deflection and rotations and The relationship between the Cartesian transverse shear strains at the middle-edge points and the natural transverse shear strains is: (59) Materials 2020, 13, 5016

of 18
The Cartesian transverse shear strains at the four sampling points are then related to the nodal displacements by: The substitute transverse shear strain matrix was obtained bȳ Here, the element is described as QLLL (for quadrilateral, bilinear deflection and rotations and linear transverse shear strain fields). The element satisfies conditions [35][36][37] n θ + n w ≥ n γ , where n θ , n w and n γ are the number of variables included in the interpolation of the rotations, the deflection and the transverse shear strains, respectively. In order to preserve the element from spurious mechanisms the full 2 ×2 quadrature for all terms were used in computation of the stiffness matrix. Since the shear forces and bending moments are constant along each natural direction, the fine meshes are required for certain applications.
The product AP −1 T in Equation (61) is The assumed transverse shear strain field can be also written in the direct form:

Results
The efficiency of the analytical approach AA presented in the paper was verified in reference to several models using different methods, both numerical and analytical. Results from the numerical methods were obtained from the commercial FE software, namely Abaqus FEA, which enables a few types of FEs to model plates. Results from the AKK method [3] was used in the comparison as the example of the quasi-analytical method, see Equation (14).
To verify the method in a broad application, key geometrical and material parameters of the samples were changed in a wide range. The following parameters were analyzed: the aspect ratio of the sample and the shear stiffnesses, i.e., D 33 , A 44 and A 55 . On the basis of preliminary tests, it was concluded that two sizes of samples should be analyzed for corrugated materials, in which a × b equals, e.g., 75 mm × 75 mm and 125 mm × 25 mm; the aspect ratios equal 1:1 and 5:1, respectively. Regarding torsion stiffness, D 33 had two values, namely, D 33 = 450 Nmm and D 33 = 900 Nmm. The transversal shear stiffness, A 44 and A 55 , varied from 5 to 250 N/mm, with a few selected values in between.
These parameters were used as an input for different methods, which are marked in the study as STRI65, S32, S34, S4/S4R, QLLL, AA and AKK. In all of those methods, the force was applied in the corners, like in Figure 2b. STRI65 refers to the FE model, in which two 6-node triangular thin shell FEs, using five degrees of freedom per node, were used, see Figure 6a. S32 refers to the FE model, in which two 3-node triangular general-purpose shell FEs were used, see Figure 6b. S34 refers to the FE model, in which four 3-node triangular general-purpose shell FEs were used, see Figure 6c. S4/S4R refers to the FE models, in which one 4-node general-purpose shell FE was used with full and reduced integration scheme, respectively, see Figure 6d. QLLL refers to the FE model, in which 4-node shell FE presented in [33,34] was used, see Equations (52)-(66) and Figure 5. AA refers to the analytical approach presented in this paper, see Equation (12), while AKK refers to the analytical approach presented in the paper of Garbowski, Gajewski and Grabski [3]. The results obtained from application all these methods are shown in Figures 7-12. All data used for generation the plots are included in Supplementary Materials. equals, e.g., 75 mm × 75 mm and 125 mm × 25 mm; the aspect ratios equal 1:1 and 5:1, respectively. Regarding torsion stiffness, had two values, namely, = 450 Nmm and = 900 Nmm. The transversal shear stiffness, and , varied from 5 to 250 N/mm, with a few selected values in between.
These parameters were used as an input for different methods, which are marked in the study as STRI65, S32, S34, S4/S4R, QLLL, AA and AKK. In all of those methods, the force was applied in the corners, like in Figure 2b. STRI65 refers to the FE model, in which two 6-node triangular thin shell FEs, using five degrees of freedom per node, were used, see Figure 6a. S32 refers to the FE model, in which two 3-node triangular general-purpose shell FEs were used, see Figure 6b. S34 refers to the FE model, in which four 3-node triangular general-purpose shell FEs were used, see Figure 6c. S4/S4R refers to the FE models, in which one 4-node general-purpose shell FE was used with full and reduced integration scheme, respectively, see Figure 6d. QLLL refers to the FE model, in which 4-node shell FE presented in [33,34] was used, see Equations (52)-(66) and Figure 5. AA refers to the analytical approach presented in this paper, see Equation (12), while AKK refers to the analytical approach presented in the paper of Garbowski, Gajewski and Grabski [3]. The results obtained from application all these methods are shown in Figures 7-12. All data used for generation the plots are included in Supplementary Materials.  In the first row of each graph (Figures 7-12), the reaction forces, , are presented for particular parameters of × , , and , while in the second row these reaction forces in reference to the analytical solutions of the AA model, see Equation (12)

Discussion
The results demonstrated in Section 3 applied to a variety of cardboard structures, which is reflected in a wide shear stiffness represented by the output reaction forces. In this study, the reaction forces were in the range from 0.5 to 2.5 N. Such comprehensive verification allows several useful conclusions to be drawn for a more robust modeling of the corrugated cardboard structures. First, it should be noted that the analysis with the QLLL element produced results that were perfectly consistent with the analytical solution AA, see  This fact was observed in each of the 196 analyzed cases. S34 method was usually the second-best FE method under consideration, comparing its accuracy with respect to the reference analytical solution, for instance see Figure 7a or Figure 12b. This was not always obvious, e.g., see Figure 11a, or true, e.g., see Figure 9b or Figure 10a. In the first row of each graph (Figures 7-12), the reaction forces, R, are presented for particular parameters of a × b, D 33 , A 44 and A 55 , while in the second row these reaction forces in reference to the analytical solutions of the AA model, see Equation (12)

Discussion
The results demonstrated in Section 3 applied to a variety of cardboard structures, which is reflected in a wide shear stiffness represented by the output reaction forces. In this study, the reaction forces were in the range from 0.5 to 2.5 N. Such comprehensive verification allows several useful conclusions to be drawn for a more robust modeling of the corrugated cardboard structures. First, it should be noted that the analysis with the QLLL element produced results that were perfectly consistent with the analytical solution AA, see  This fact was observed in each of the 196 analyzed cases. S34 method was usually the second-best FE method under consideration, comparing its accuracy with respect to the reference analytical solution, for instance see Figure 7a or Figure 12b. This was not always obvious, e.g., see Figure 11a, or true, e.g., see Figure 9b or Figure 10a.
In the results, there were crucial differences while comparing rectangular (125 mm × 25 mm) and square (75 mm × 75 mm) samples. The biggest divergence from the analytical solution were observed in STRI65 and S4/S4R methods for rectangular sample, especially in cases, when the sample had a large transversal shear stiffness along short dimension and small transversal shear stiffness along long dimension, see Figure 11b (A 44 = 5 N/mm), and Figure 12a (A 55 > 50 N/mm). In these cases, the difference in relation to the reference solution was even four times greater. Moreover, for these cases, the reaction force was relatively low, i.e., about 0.5 N. The STRI65 and S4/S4R approached greatly overestimated transversal shear stiffness of the samples.
On contrary, also in rectangular samples, if the material had large transversal shear stiffness along longer dimension of the sample, the differences between all methods and the reference solution were very low, see Figures 10b and 12b. In all these methods, apart from the S34 case, the reactions changed their relation to the reference solution by crossing the value of 1.0. In other words, the analyzed methods for some values of A 55 overestimated the analytical solution and for other values of A 55 , underestimated. Going further, while considering AKK in the rectangular samples; the method in 50% of cases may be considered as equally good comparing to the S34 method, but both methods were worse than S32. In the other 50% of cases, AKK was worse than S32 and S34, which were in these cases equally good methods. In all cases, AKK was more accurate than STRI65 and S4/S4R.
In the results for square samples, over an entire range of transversal shear stiffness considered, AKK was competitive to the S34 method, but usually slightly worse. In these cases, the S32 method was usually much worse than AKK or S34. This was shown, for example, in Figure 8a.
It is worth noting that for the square sample, due to the symmetry of the material (A 44 vs. A 55 ), and thus the symmetry of the results, these plots were not presented to avoid repetition.
One of the limitations of the study was considering the finite ranges of the material parameters of D 33 , A 44 and A 55 . The ranges adopted here, even if they did not appear to be physically reasonable for the typical properties of corrugated board, may be of interest when analyzing sandwich panels made of other materials. Another limitation was fixing the two sets of in-plane dimensions and only one thickness of the samples. The selection of the samples dimensions was followed by a practical laboratory tests of corrugated board [38], while the selected thickness is one of the most commonly used in the corrugated cardboards packaging industry.

Conclusions
In the paper, the classic governing equations of the Reissner-Mindlin plate were presented with a particular application to composite laminates. The formula for torsion of the orthotropic plate including transversal shear stiffness was defined and compared with the recent considerations from the literature. Transverse shear is an important component if the analyzed plate is thick with respect to the in-plane dimensions and/or its core has significantly lower stiffness than the outer faces. Different methods for modeling such a plate, both numerical and analytical, were selected to compute their accuracies comparing to the exact, analytical solution. Numerical methods used varied in the type of finite element used and number of elements taken into consideration. To cover a wide range of material properties, the computations were performed for selected in-plane and transversal shear moduli. Samples with two different aspect ratios, rectangular and square, were tested.
It was concluded that the 4-noded Reissner-Mindlin plate quadrilateral with assumed linear transverse shear strain field gave the best performance. Results obtained via this method were exact with the analytical solution. The results from other methods, with 6-node triangular elements or 4-node quadrilateral elements (both, with full integration or reduced), with commonly used finite elements were essentially worse. In the worst cases, the overestimation of the force was even up to four times. In this comparison, the 3-node triangular element performance was surprisingly moderately good. Moreover, the differences between rectangular and square samples were significant. It should be noted that in rectangular samples, if the material had a large transversal shear stiffness along a longer dimension of the sample, the differences of all methods to the reference solution were very low.
Supplementary Materials: The following are available online at http://www.mdpi.com/1996-1944/13/21/5016/s1. All computational results are presented in CSV files, including the data, based on which Figures 7-12 were generated. In these files, 1st and 2nd columns are a and b dimensions of the box, respectively, 3rd column is the thickness of the wall, the materials parameters: G 12 , G 13 , G 23 , D 33 , A 44 and A 55 . are given in the columns 4th to 9th, respectively. The columns 10th-16th include results of the calculations using: STRI65, S32, S34, S4/S4R, QLL, AKK and AA, respectively.