An Equivalent Homogenization Theoretical Method for Composite Sandwich Cylinders Subjected to Pure Bending

Abstract

An equivalent theoretical homogenization method was proposed for composite sandwich cylinders subjected to pure bending. Firstly, based on a homogeneous orthotropic layer hypothesis, the trapezoidal corrugated sandwich core was found to be equivalent in a homogenization orthotropic layer with the nine equivalent mechanical properties. Then, Lekhnitskii’s theory, based on a unified connection parameter method, was introduced and applied in the equivalent composite sandwich cylinder. The method developed by Lekhnitskii is suitable for arbitrary combinations of winding layers with different winding angles and materials. Additionally, the bending stiffness of the equivalent sandwich cylinder could be calculated. By developing user subroutine of UMAT, the numerical calculation results were in a good agreement with the results of the proposed method. Further, according to the Hill–Tsai strength criterion and the maximum strain criterion, parametric study was done for specified bending stiffness and specified bending strength. The results show that the influence of core parameters on the specified bending stiffness and strength are lower than that of the skin parameters. Additionally, larger skin thickness and smaller winding angles could improve the specified bending stiffness and specified bending strength of the composite corrugated sandwich cylinders.

1. Introduction

Composite sandwich cylinders made of internal and outer skins and a core material with certain topological configurations are widely used in aerospace, aircraft, and building industries due to their high specified stiffness, specified strength, damage tolerance, and blast resistance [1,2,3]. In the past decades, many different lattice core structures have been proposed. It has been proved that corrugated and periodic lattice cores, including Lattice block [4], Tetrahedral truss [5], Pyramidal truss [6], Kagome lattice [7], Woven lattice [8], and Isogrid [9,10], have the potential to enhance numerous applications in aerospace structures. Among these cores, corrugated core is the most widely used because of its simple forming process, low cost, and relatively good mechanical properties. However, most published reports [4,7,11,12,13] deal with axial-compressive behaviors of the composite sandwich cylinders, but rare investigations are involved in their bending behaviors, which is a significant performance index.

A composite thick wall pipe is generally anisotropic. After deformation caused by various types of loading, such as internal or external pressure, axial pressure, or torsion, the axis of the composite thick wall pipe is still straight. However, under bending condition, deflection occurs on the axis of the thick wall pipe, which makes the deformation prediction more complex. Earlier, Lekhnitskii [14] used the flexibility method of a double stress function to obtain two partial differential equations for the case of cylindrical anisotropy. Later, Jolicoeur and Cardou [15] developed and expanded Lekhnitskii’s theory. The developed theory is only suitable for spiral wound layers and not for special winding layers with orientation n of 0° and 90°. Xia et al. [16] treated cross-winding layers to 90° winding layers. A single stress function was used to establish differential equations because out-of-plane shear stresses are zero, which limits the scope of application of their method. In 2014, Zhang and Hoa [17,18] proposed a unified connection parameter method. This method overcomes the problem of singular parameters in the presence of a special winding layer and is suitable for composite tubes of general stacking sequence including 0°, 90°, and [$\phi$/−$\phi$].

The periodic topological core structure in composite sandwich cylinders is another complex factor. In the past few decades, the research on cores’ equivalent theory has mainly focused on the in-plane stiffness of sandwich panels. For corrugated cores, Libove and Hubka et al. [19] proposed that shear effects were an important phenomenon that could not be neglected. Based on the Classical Plate Theory (CPT), Xia et al. [20] derived the in-plane equivalent mechanical properties for different corrugation cores. Bartolozzi et al. [21] derived the in-plane equivalent mechanical properties of a sinusoidal-shaped corrugated core using an energy method. Then, Bartolozzi et al. further extended the energy method to make it suitable for arbitrary shapes of corrugated cores [22]. Based on the classical Euler–Bernoulli and line hypotheses, Magnucka-Blandzi et al. [23,24] analyzed the bending and buckling behavior and the effect of transverse shearing deformation on sandwich beams with sinusoidal corrugated cores. They also studied the effect of transverse shearing deformation for short and long beams. For corrugated cores, equivalent out-of-plane properties are also important and cannot be ignored [25]. Therefore, Zhang and Lei et al. [26] studied the out-of-plane compressive performance and energy absorption of multi-layer graded sinusoidal corrugated sandwich panels. Ghahfarokhi et al. [27] use a new approach to reveal the damage mode of global buckling and the equivalent stiffness of composite sandwich cylindrical shells with lattice cores under uniaxial compression. In the context of shear deformation theory, Yadav et al. [28] studied circular cylindrical sandwich shells with open/closed cellular core the face sheets. Additionally, their effective mechanical properties were determined using the Gibson and Ashby model and the rule of mixture. Guo et al. [29] proposed a unified semi-analytical method for dynamic homogenization and vibration of lattice-truss-core sandwich beams. Based on the mechanics of materials approach and the classical plate theory, Pydah et al. [30] analytically investigate infinitesimal cylindrical bending deformations of two-layered triangular corrugated and webcore linearly elastic sandwich panels. Kolpakov et al. [31] derived the equivalent stiffnesses for corrugated plates by using a two-step approximation method. For corrugated core sandwich panels, Buannic et al. [32] applied a specific study to derive the transverse shear stiffness and determined an equivalent Reissner–Mindlin homogeneous plate.

This paper proposes an equivalent homogenization theoretical method for composite sandwich cylinder under pure bending. Firstly, based on a homogeneous orthotropic layer hypothesis, the trapezoidal-shaped corrugated sandwich core was approximately equivalent to an orthotropic homogenization layer. This layer was described with nine equivalent mechanical properties. Then, Lekhnitskii’s theory, based on a unified connection parameter method, was introduced and applied in the equivalent composite sandwich cylinder. By developing a user subroutine of UMAT, the results of the equivalent homogenization theoretical method were compared with the numerical results. Further, according to the characteristics of the sandwich cylinder, two different strength criteria for skin and core are discussed. The influence of variable parameters on specified bending stiffness and strength is discussed with a parametric study.

2. Materials and Methods

2.1. Analytical Methods for Composite Cylinders Subjected to Pure Bending

After years of development, the range of application of the elasticity theory for composite cylinders subjected to pure bending has been continuously expanded. This elasticity theory is an important basis for the sandwich cylinder subjected to pure bending. Additionally, the theory was introduced in this section.

2.1.1. Analytical 3D Elasticity Method by Lekhnitskii

A composite tube with multiple layers under a pure bending load $M_x$ and $M_y$ is shown in Figure 1 [18], where $b_n$, $b_{n+1}$ and ${\phi }_n(n=1,2,\cdot \cdot \cdot ,N)$ are the inner and outer radius of the nth layer and its filament winding angle, respectively. Additionally, ${\sigma }_z$, is the main bending stress component along the axis in the composite pipe. The in-plane stress components of ${\sigma }_r$, ${\sigma }_{\theta }$, and ${\tau }_{r\theta }$ are parallel to the cross-section are the secondary stress components, and the stress components of ${\tau }_{\theta z}$ and ${\tau }_{rz}$ are the out-of-plane stress perpendicular to the cross-section.

Relationship between strain and stress of single-layer composite winding layer in material principal axis coordinate system can be expressed as

$$\left\{\begin{array}{l}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\\ {\gamma }_{23}\\ {\gamma }_{13}\\ {\gamma }_{12}\end{array}\right\}=S\left\{\begin{array}{l}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\tau }_{23}\\ {\tau }_{13}\\ {\tau }_{12}\end{array}\right\}=\left[\begin{array}{llllll}{S}_{11}& S_{12}& S_{13}& 0& 0& 0\\ & S_{22}& S_{23}& 0& 0& 0\\ & & S_{33}& 0& 0& 0\\ & & & S_{44}& 0& 0\\ & & \mathrm{sym}& & S_{55}& 0\\ & & & & & S_{66}\end{array}\right] \left\{\begin{array}{l}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\tau }_{23}\\ {\tau }_{13}\\ {\tau }_{12}\end{array}\right\}$$

(1)

where

$$\begin{array}{l}{S}_{11} =\frac{1}{E_{11}},S_{22}=\frac{1}{E_{22}},S_{33}=\frac{1}{E_{33}},S_{44}=\frac{1}{G_{23}},S_{55}=\frac{1}{G_{13}},S_{66}=\frac{1}{G_{12}}\\ S_{12}=\frac{-{\nu }_{21}}{E_{22}}=\frac{-{\nu }_{12}}{E_{11}},S_{13}=\frac{-{\nu }_{31}}{E_{33}}=\frac{-{\nu }_{13}}{E_{11}},S_{23}=\frac{-{\nu }_{32}}{E_{33}}=\frac{-{\nu }_{23}}{E_{22}}\end{array}$$

(2)

Based on the winding angle of ${\phi }_n(n=1,2,\cdot \cdot \cdot ,N)$, the constitutive equation of each layer transformed from the principal material coordinate system to the cylindrical system can be obtained as follows [18]:

$$\left\{\begin{array}{l}{\epsilon }_r\\ {\epsilon }_{\theta }\\ {\epsilon }_z\\ {\gamma }_{\theta z}\\ {\gamma }_{rz}\\ {\gamma }_{r\theta }\end{array}\right\}=C\left\{\begin{array}{l}{\sigma }_r\\ {\sigma }_{\theta }\\ {\sigma }_z\\ {\tau }_{\theta z}\\ {\tau }_{rz}\\ {\tau }_{r\theta }\end{array}\right\}=\left[\begin{array}{llllll}{C}_{11}& C_{12}& C_{13}& C_{14}& 0& 0\\ & C_{22}& C_{23}& C_{24}& 0& 0\\ & & C_{33}& C_{34}& 0& 0\\ & & & C_{44}& 0& 0\\ & & \mathrm{sym}& & C_{55}& C_{56}\\ & & & & & C_{66}\end{array}\right]\left\{\begin{array}{l}{\sigma }_r\\ {\sigma }_{\theta }\\ {\sigma }_z\\ {\tau }_{\theta z}\\ {\tau }_{rz}\\ {\tau }_{r\theta }\end{array}\right\}$$

(3)

where

$$C={\mathrm{T}}_{}^{T}S{\mathrm{T}}_{}^{}$$

(4)

in which

$$\mathrm{T}=\left[\begin{array}{cccccc}0& {\sin }^2\phi & {\cos }^2\phi & -2\sin \phi \cos \phi & 0& 0\\ 0& {\cos }^2\phi & {\sin }^2\phi & 2\sin \phi \cos \phi & 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& -\sin \phi & -\cos \phi \\ 0& 0& 0& 0& \cos \phi & -\sin \phi \\ 0& \sin \phi \cos \phi & -\sin \phi \cos \phi & {\sin }^2\phi -{\cos }^2\phi & 0& 0\end{array}\right]$$

(5)

To simplify the expression, the following reduced elastic constants ${\beta }_{ij}$ were defined by Lekhnitskii [14]

$${\beta }_{ij}=C_{ij}-\frac{C_{i3}{C}_{3j}}{C_{33}}$$

(6)

For the pure bending problem of composite pipes, Lekhnitskii [14] first proposed the flexibility method containing two stress functions, $F\left(r,\theta \right)$ and $\Psi \left(r,\theta \right)$. The in-plane stress component is expressed as the following form of stress function $F\left(r,\theta \right)$.

$${\sigma }_r=\frac{1}{r}\frac{\partial F}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}F}{\partial {\theta }^2}, {\sigma }_{\theta }=\frac{{\partial }^{2}F}{\partial r^2}, {\tau }_{r\theta }=-\frac{{\partial }^2}{\partial r\partial \theta }\left(\frac{F}{r}\right)$$

(7)

Meanwhile, the out-of-plane stress component is expressed as follows.

$${\tau }_{rz}=\frac{1}{r}\frac{\partial \Psi }{\partial \theta }, {\tau }_{\theta z}=-\frac{\partial \Psi }{\partial r}$$

(8)

Additionally, the three-dimensional elastic mechanical problem of pure bending of a composite pipe is a generalized plane strain problem, and its bending stress can be expressed as follows:

$${\sigma }_z=\frac{1}{C_{33}}\left(k_{x}r\sin \theta +k_{y}r\cos \theta -C_{13}{\sigma }_{\mathrm{r}}-C_{23}{\sigma }_{\theta }-C_{34}{\tau }_{\theta z}\right)$$

(9)

A system of ordinary differential equations can be obtained with separate variables, and the solution is sought as follows.

$$F=f\left(r\right)\left(k_x\sin \theta -k_y\cos \theta \right)$$

(10)

$$\Psi =\psi \left(r\right)\left(k_x\sin \theta -k_y\cos \theta \right)$$

(11)

where $k_x$ and $k_y$ are the curvature in the plane perpendicular to directions x and y, respectively.

According to the compatibility equations [14], the following ordinary differential equation about the stress function can be obtained.

$$\{\begin{array}{l}(\frac{{\beta }_{22}}{r}\left(r\frac{d^4}{d{r}^4}+2\frac{d^3}{d{r}^3}\right)-\frac{{\beta }_{11}+2{\beta }_{12}+{\beta }_{66}}{r^4}\left(r^2\frac{d^2}{d{r}^2}-r\frac{d}{dr}+1\right))f\left(r\right)\\ +(\frac{{\beta }_{14}}{r^2}\left(\frac{d}{dr}+r\frac{d^2}{d{r}^2}\right)-\frac{{\beta }_{24}}{r}\left(r\frac{d^3}{d{r}^3}+2\frac{d^2}{d{r}^2}\right)+\frac{{\beta }_{56}}{r^2}\frac{d}{dr})\psi \left(r\right)\\ =\frac{2}{r}\frac{C_{13}-C_{23}}{C_{33}}\\ (\frac{{\beta }_{14}}{r^3}\left(r^2\frac{d^2}{d{r}^2}-r\frac{d}{dr}-1\right)+\frac{{\beta }_{24}}{r}\left(-r\frac{d^3}{d{r}^3}+\frac{d^2}{d{r}^2}\right)+\frac{{\beta }_{56}}{r^3}\left(-r\frac{d}{dr}-1\right))f\left(r\right)\\ +\left(\frac{{\beta }_{44}}{r}\left(r\frac{d^2}{d{r}^2}+\frac{d}{dr}\right)-\frac{{\beta }_{55}}{r^2}\right)\psi \left(r\right)=\frac{2C_{34}}{C_{33}}\end{array}$$

(12)

2.1.2. Solutions for Composite Tube with Multi-Layer

Based on Lekhnitskii [14], Jolicoeu and Cardou [15] considered the fourth-order ordinary-differential-compatibility equation of pure bending problem of composite tube with multi-layer and arbitrary winding angle. System (12) is of the Cauchy–Euler type. The form of the solution of $f\left(r\right)=K{r}^{m+1}$ and $\psi \left(r\right)=Kg{r}^m$ were assumed. The complete solutions of the pure bending system for the stress function are

$$F=\left(k_x\sin \theta -k_y\cos \theta \right)\left({\displaystyle \sum _{i=1}^4\frac{K_i}{m_i}{r}^{m_i+1}+\frac{{\mu }_1}{2}{r}^3}\right)$$

(13)

$$\Psi =\left(k_x\sin \theta -k_y\cos \theta \right)\left({\displaystyle \sum _{i=1}^4{K}_i{g}_i{r}^{m_i}+{\mu }_2{r}^2}\right)$$

(14)

where $m_i$ are four characteristic roots given by

$$\{\begin{array}{l}{m}_1\\ m_2\end{array}=\pm \sqrt{\frac{-b+\sqrt{b^2-4ac}}{2a}}$$

(15)

$$\{\begin{array}{l}{m}_3\\ m_4\end{array}=\pm \sqrt{\frac{-b-\sqrt{b^2-4ac}}{2a}}$$

(16)

with

$$\begin{array}{ll}a=& {\beta }_{22}{\beta }_{44}-{\beta }_{24}^2\\ b=& {\beta }_{24}\left(2{\beta }_{14}+{\beta }_{24}+2{\beta }_{56}\right)\\ & -{\beta }_{44}\left({\beta }_{11}+2{\beta }_{12}+{\beta }_{22}+{\beta }_{66}\right)-{\beta }_{22}{\beta }_{55}+{\beta }_{14}^2\\ c=& {\beta }_{55}\left({\beta }_{11}+2{\beta }_{12}+{\beta }_{22}+{\beta }_{66}\right)-{\beta }_{56}^2\end{array}$$

(17)

${\mu }_1$ and ${\mu }_2$ are the two unknown constants. The solution is written in matrix form as

$$\left\{\begin{array}{l}{\mu }_1\\ {\mu }_2\end{array}\right\}={\left[\begin{array}{cc}-2{\beta }_{14}-6{\beta }_{24}+{\beta }_{56}& 4{\beta }_{44}-{\beta }_{55}\\ -{\beta }_{11}-2{\beta }_{12}+3{\beta }_{22}-{\beta }_{66}& 2{\beta }_{14}-2{\beta }_{24}+{\beta }_{56}\end{array}\right]}^{-1}\frac{1}{C_{33}}\left\{\begin{array}{c}2C_{34}\\ C_{13}-C_{23}\end{array}\right\}$$

(18)

$g_i$ are four connection parameters obtained from

$$g_i=\frac{{\beta }_{24}{m}_i^2+\left({\beta }_{14}+{\beta }_{24}\right)m_i-{\beta }_{56}}{{\beta }_{44}{m}_i^2-{\beta }_{55}}; i=1 to 4$$

(19)

Additionally, $K_i$ are four combination coefficients that are calculated with the continuity boundary condition.

2.1.3. The Unified Connection Parameter Method

According to the above method proposed by Jolicoeur and Cardou [15], Zhang et al. [18] found the connection parameters of $g_3$ and $g_4$ are infinite in $\underset{\phi \to 0,90}{\left|\lim g_3\right|}=\infty$ and $\underset{\phi \to 0,90}{\left|\lim g_4\right|}=\infty$. In order to overcome this issue, the unified connection parameter method was proposed by

$$K_j^{\ast }=\{\begin{array}{l}{K}_j{g}_j, \phi \ne 90 or 0\\ K_j{g}_j, \phi =90 or 0\end{array} j=3,4$$

(20)

$$g_j^{\ast }=\{\begin{array}{l}{g}_j^{-1}, \phi \ne 90 or 0\\ 0, \phi =90 or 0\end{array} j=3,4$$

(21)

Then, the four combination coefficients of $K_1,K_2,K_3,K_4$ are transformed into $K_1,K_2,K_3^{\ast },K_4^{\ast }$. Additionally, the six stresses are as follows.

$$\{\begin{array}{l}{\sigma }_r=\left(k_x\sin \theta -k_y\cos \theta \right)\left({\displaystyle \sum _{i=1}^2{K}_i}{r}^{m_i-1}+{\displaystyle \sum _{j=3}^4{K}_j^{\ast }}{g}_j^{\ast }{r}^{m_j-1}+{\mu }_{1}r\right)\\ {\sigma }_{\theta }=\left(k_x\sin \theta -k_y\cos \theta \right)\left({\displaystyle \sum _{i=1}^2{K}_i\left(m_i+1\right)}{r}^{m_i-1}+{\displaystyle \sum _{j=3}^4{K}_j^{\ast }}{g}_j^{\ast }\left(m_j+1\right)r^{m_j-1}+3{\mu }_{1}r\right)\\ {\tau }_{r\theta }=\left(k_x\cos \theta +k_y\sin \theta \right)\left(-{\displaystyle \sum _{i=1}^2{K}_i}{r}^{m_i-1}-{\displaystyle \sum _{j=3}^4{K}_j^{\ast }}{g}_j^{\ast }{r}^{m_j-1}-{\mu }_{1}r\right)\\ {\tau }_{rz}=\left(k_x\cos \theta +k_y\sin \theta \right)\left({\displaystyle \sum _{i=1}^2{K}_i}{g}_i{r}^{m_i-1}+{\displaystyle \sum _{j=3}^4{K}_j^{\ast }}{r}^{m_j-1}+{\mu }_{2}r\right)\\ {\tau }_{\theta z}=\left(k_x\sin \theta -k_y\cos \theta \right)\left(-{\displaystyle \sum _{i=1}^2{K}_i}{g}_i{m}_i{r}^{m_i-1}-{\displaystyle \sum _{j=3}^4{K}_j^{\ast }}{m}_j{r}^{m_j-1}-2{\mu }_{2}r\right)\\ {\sigma }_z=\frac{1}{C_{33}}\left(k_{x}r\sin \theta -k_{y}r\cos \theta -C_{13}{\sigma }_r-C_{23}{\sigma }_{\theta }-C_{34}{\tau }_{\theta z}\right)\end{array}$$

(22)

The displacements are expressed as follows

$$\{\begin{array}{l}{u}_r=\left(k_x\sin \theta -k_y\cos \theta \right)\left(-\frac{z^2}{2}+{\displaystyle \sum _{i=1}^2{K}_i}{U}_i^{\prime }{r}^{m_i}+{\displaystyle \sum _{j=3}^4{K}_j^{\ast }}{U}_j^{\ast }{r}^{m_j}+U_5^{\prime }{r}^2+\upsilon \right)\\ u_{\theta }=\left(k_x\cos \theta +k_y\sin \theta \right)\left(-\frac{z^2}{2}+{\displaystyle \sum _{i=1}^2{K}_i}{V}_i^{\prime }{r}^{m_i}+{\displaystyle \sum _{j=3}^4{K}_j^{\ast }}{V}_j^{\ast }{r}^{m_j}+V_5^{\prime }{r}^2+\upsilon \right)\\ w=zr\left(k_x\sin \theta -k_y\cos \theta \right)+\left(k_x\cos \theta +k_y\sin \theta \right)\left({\displaystyle \sum _{i=1}^2{K}_i}{W}_i^{\prime }{r}^{m_i}+{\displaystyle \sum _{j=3}^4{K}_j^{\ast }}{W}_j^{\ast }{r}^{m_j}+W_5^{\prime }{r}^2\right)\end{array}$$

(23)

where functions $U$, $V$, and $W$ are finally written as

$$\{\begin{array}{l}{U}_i^{\prime }=\frac{1}{m_i}\left({\beta }_{11}+{\beta }_{12}\left(m_i+1\right)-{\beta }_{14}{g}_i{m}_i\right)\\ V_i^{\prime }=\frac{1}{m_i}\left({\beta }_{11}+{\beta }_{12}-{\beta }_{22}{m}_i\left(m_i+1\right)-g_i{m}_i\left({\beta }_{14}-{\beta }_{24}{m}_i\right)\right)\\ W_i^{\prime }=\frac{1}{m_i}\left({\beta }_{55}{g}_i-{\beta }_{56}\right)\end{array}$$

(24)

$$\{\begin{array}{l}{U}_j^{\ast }=\frac{1}{m_j}\left(\left({\beta }_{11}+{\beta }_{12}\left(m_i+1\right)\right)g_j^{\ast }-{\beta }_{14}{m}_i\right)\\ V_j^{\ast }=\frac{1}{m_j}\left(\left({\beta }_{11}+{\beta }_{12}-{\beta }_{22}{m}_i\left(m_i+1\right)\right)g_j^{\ast }-m_i\left({\beta }_{14}-{\beta }_{24}{m}_i\right)\right)\\ W_j^{\ast }=\frac{1}{m_j}\left({\beta }_{55}-{\beta }_{56}{g}_j^{\ast }\right)\end{array}$$

(25)

$$\{\begin{array}{l}{U}_5^{\prime }=\frac{1}{2}\left({\mu }_1\left({\beta }_{11}+3{\beta }_{12}\right)-2{\beta }_{14}{\mu }_2+\frac{C_{13}}{C_{33}}\right)\\ V_5^{\prime }=\frac{1}{2}\left({\mu }_1\left({\beta }_{11}+{\beta }_{12}-6{\beta }_{22}\right)-2{\mu }_2\left({\beta }_{14}-2{\beta }_{24}\right)+\frac{C_{13}-2C_{23}}{C_{33}}\right)\\ W_5^{\prime }=\frac{1}{2}\left({\beta }_{55}{\mu }_2-{\beta }_{56}{\mu }_1\right)\end{array}$$

(26)

where $i=1,2$ and $j=3,4$.

In the context of no relative slip between layers in the composite tubecase, the stresses ${\sigma }_r$, ${\tau }_{r\theta }$, and ${\tau }_{rz}$ and displacements $u_r$, $u_{\theta }$, and $w$ must meet the continuity condition. For a composite tube with N layers, there are two free surfaces yielding two equations each, plus (N − 1) interfaces giving five equations each. The value ${\nu }_1$ is set to zero [15]. This gives 5N equations for the same number of arbitrary constants ($K_1,K_2,K_3^{\ast },K_4^{\ast }$, and $\nu$ for each cylinder). Then, the 5N arbitrary constants can be calculated. The end conditions of pure bending are given by

$$\{\begin{array}{l}{M}_x={\displaystyle \sum _{n=0}^N{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{b_n}^{b_{n+1}}{\sigma }_z{r}^2\sin \theta }}}d\theta dr\\ M_y=-{\displaystyle \sum _{n=0}^N{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{b_n}^{b_{n+1}}{\sigma }_z{r}^2\cos \theta }}}d\theta dr\end{array}$$

(27)

By using Equations (9) and (27), the bending stiffness $\left(EI\right)$ is written as [18]

$$\{\begin{array}{l}{M}_x=\left(EI\right)k_x\\ M_y=\left(EI\right)k_{\mathrm{y}}\end{array}$$

(28)

$$\begin{array}{ll}\left(EI\right)=& {\displaystyle \sum _{n=1}^N\frac{\pi }{{\left(C_{33}\right)}_n}\{{\displaystyle \sum _{i=1}^2{\left(K_i\right)}_n[{\left(C_{13}\right)}_n+{\left(C_{23}\right)}_n\left({\left(m_i\right)}_n+1\right)}}\\ & -{\left(C_{34}\right)}_n{\left(m_i\right)}_n{\left(g_i\right)}_n]\frac{b_n^{{\left(m_i\right)}_n+2}-b_{n+1}^{{\left(m_i\right)}_n+2}}{{\left(m_i\right)}_n+2}\\ & +{\displaystyle \sum _{j=3}^4{\left(K_j^{\ast }\right)}_n[{\left(C_{13}\right)}_n{\left(g_j^{\ast }\right)}_n+{\left(C_{23}\right)}_n\left({\left(m_j\right)}_n+1\right)}{\left(g_j^{\ast }\right)}_n\\ & -{\left(C_{34}\right)}_n{\left(m_j\right)}_n]\frac{b_n^{{\left(m_j\right)}_n+2}-b_{n+1}^{{\left(m_j\right)}_n+2}}{{\left(m_j\right)}_n+2}\\ & +\left[{\left({\mu }_1\right)}_n\left({\left(C_{13}\right)}_n+3{\left(C_{23}\right)}_n\right)-2{\left({\mu }_2\right)}_n{\left(C_{34}\right)}_n-1\right]\frac{b_n^4-b_{n+1}^4}{4}\}\end{array}$$

(29)

2.2. Equivalent Properties of Trapezoidal Corrugated Cores

In recent years [12,33,34], the sandwich structures with trapezoidal corrugated-cores, as shown in Figure 2, have been of great concern due to their simple forming process, low cost, and sufficient bonding area. Therefore, the precise equivalent parameters of corrugated core have an important influence on the strength and stiffness design of sandwich structure.

Unlike triangular-shaped cores and sinusoidal cores, the trapezoidal cores are composed of two parts: the horizontal segment and the corrugated segment. Based on Mindlin–Reissner theory, H. Mohammadi [35] and M. Shaban [25] established a quarter trapezoidal cell model and the boundary conditions (BCs) of a dummy moment $M_0$. A vertical force $V$ and a horizontal force $H$ were applied at the end of the horizontal part. Although their calculated results of equivalent mechanical properties are consistent with the simulation under the same BCs, their methods are only suitable for the trapezoidal corrugated cores without skins. The actual BCs are different between trapezoidal-shaped corrugated sandwich plates with two skins and without skins. In order to accurately calculate the equivalent parameters, the correct boundary condition should be to apply forces and constraints at the conjunction points of the core with skins [22].

As shown in Figure 3, a half-trapezoidal cell model represented by its centerline was proposed by Giorgio Bartolozzi [21,22], where p₁, p₂, and p₃ are the lengths of three parts (Part O–A, Part A–B, and Part B–C, respectively) along the x-direction. The corrugation cores can be considered periodic with period 2 × p, dip angle β, and height h.

Generally, the corrugated cores are made of lightweight metal by a stamping and folding process [36] or carbon fiber woven composite by a mold-pressing method [12]. Additionally, the carbon fiber woven composite with symmetric orthogonal ply has an in-plane quasi-isotropic characteristic, which could be treated as a metal material.

2.2.1. Transverse Shear Modulus G_xz

The corrugation part (Part A–B in Figure 3) of a trapezoidal core is the main deformation element that determines the mechanical characteristics of the whole core. Indeed, it is necessary to calculate the deformation of Part A–B.

As shown in Figure 4, the analytical model of Part A–B has thickness t_c and unit width b = 1 in the y-direction. The origin of the XYZ coordinate system is set at point A, which is the conjunction points of the core with skins, and there the structure is clamped. In order to determine the equivalent shear modulus G_xz, a horizontal force $H$, a dummy moment $M_0$, and a vertical force $V$ are applied at point B. Then, the Part A-B will only produce the horizontal displacement ${\delta }_H$, which is also called pure shear deformation (see Figure 5).

The inner forces of arbitrary point in Part A–B are as follows.

$$\begin{array}{l}M=H\left(h-x\times \tan \beta \right)+V\left(p_2-x\right)-M_0\\ N=H\cos \beta -V\sin \beta \\ T=H\sin \beta +V\cos \beta \end{array}$$

(30)

Based on the Castilian second theorem [11], the rotation displacement, vertical displacement, and horizontal displacement of the free end of point B can be derived as:

$$\{\begin{array}{l}{\delta }_H={\displaystyle {\int }_0^{p_2}\left(\frac{M}{EI}\frac{\partial M}{\partial H}+\frac{N}{EA}\frac{\partial N}{\partial H}+\frac{T}{G{A}^{\prime }}\frac{\partial T}{\partial H}\right)}\frac{dx}{\cos \phi }\\ {\delta }_V={\displaystyle {\int }_0^{p_2}\left(\frac{M}{EI}\frac{\partial M}{\partial V}+\frac{N}{EA}\frac{\partial N}{\partial V}+\frac{T}{G{A}^{\prime }}\frac{\partial T}{\partial V}\right)}\frac{dx}{\cos \phi }\\ {\delta }_{M_0}={\displaystyle {\int }_0^{p_2}\left(\frac{M}{EI}\frac{\partial M}{\partial M_0}+\frac{N}{EA}\frac{\partial N}{\partial M_0}+\frac{T}{G{A}^{\prime }}\frac{\partial T}{\partial M_0}\right)}\frac{dx}{\cos \phi }\end{array}$$

(31)

According to the boundary conditions shown in Equation (32), the horizontal displacement ${\delta }_H$ can be calculated from the equations shown in Equation (33).

$${\delta }_V=0, {\delta }_{M_0}=0, H=1$$

(32)

$$\left[\begin{array}{l}{\delta }_H\\ {\delta }_V\\ {\delta }_{M_0}\end{array}\right]=\left[\begin{array}{l}{C}_{1,1}\hspace{1em}{C}_{1,2}\hspace{1em}{C}_{1,3}\\ \hspace{1em}\hspace{1em} C_{2,2}\hspace{1em}{C}_{2,3}\\ \mathrm{sym}\hspace{1em}\hspace{1em}\hspace{1em} C_{3,3}\end{array}\right]\left[\begin{array}{l}H\\ V\\ M_0\end{array}\right]$$

(33)

where

$$\{\begin{array}{l}{C}_{1,1}=\frac{h^2{p}_2-h\tan \beta p_2{}^2+\frac{{\tan }^2\beta }{3}{p}_2{}^3}{EI\cos \beta }+\frac{\cos \beta \cos \beta p_2}{\cos \beta EA}+\frac{\sin \beta \sin \beta p_2}{\cos \beta G{A}^{\prime }}\\ C_{1,2}=\frac{h{p}_2{}^2-\frac{p_2\tan \beta +h}{2}{p}_2{}^2+\frac{\tan \beta }{3}{p}_2{}^3}{EI\cos \beta }-\frac{\sin \beta \cos \beta p_2}{EA\cos \beta }+\frac{\cos \beta \sin \beta p_2}{G{A}^{\prime }\cos \beta }\\ C_{1,3}=\frac{-h{p}_2+\frac{\tan \beta }{2}{p}_2{}^2}{EI\cos \beta }\\ C_{2,2}=\frac{\frac{p_2{}^3}{3}}{EI\cos \beta }+\frac{\sin \beta \sin \beta p_2}{EA\cos \beta }+\frac{\cos \beta \cos \beta p_2}{G{A}^{\prime }\cos \beta }\\ C_{2,3}=\frac{-\frac{p_2{}^2}{2}}{EI\cos \beta }\\ C_{3,3}=\frac{p_2}{EI\cos \beta }\end{array}$$

(34)

and $EI=EA\times \frac{t_c{}^2}{12}$, $G{A}^{\prime }=EA\times \frac{\kappa }{2\left(1+\nu \right)}$, and $\kappa$ is the shear correction factor of the beam ($\kappa =0.85$). The horizontal displacement can be found with

$${\delta }_H=\frac{\left|\begin{array}{l}{C}_{1,1}\hspace{1em}{C}_{1,2}\hspace{1em}{C}_{1,3}\\ \hspace{1em}\hspace{1em}{C}_{2,2}\hspace{1em}{C}_{2,3}\\ \mathrm{sym}\hspace{1em}\hspace{1em}\hspace{1em} C_{3,3}\end{array}\right|}{\left|\begin{array}{l}{C}_{2,2} C_{2,3}\\ C_{2,3} C_{3,3}\end{array}\right|}H$$

(35)

Therefore, the equivalent transverse shear modulus of Part A–B can be obtained as follows:

$$G_{xz2}=\frac{{\tau }_{xz}}{{\gamma }_{xz}}=\frac{F_x}{A_{xy}}/\frac{{\delta }_x}{l_z}=\frac{H}{p_2}/\frac{{\delta }_H}{2h}$$

(36)

In the case of containing upper and lower skins, since the horizontal sections (O–A and B–C) of the trapezoidal corrugated core are constrained by the skin and have little effect on the transverse shear deformation of the core under the action of the horizontal force in the x direction, the equivalent transverse shear modulus of the whole corrugated core can be derived with

$$G_{xz}=G_{xz1}\frac{p_1}{p}+G_{xz2}\frac{p_2}{p}+G_{xz3}\frac{p_3}{p}\approx G_{xz2}\frac{p_2}{p}$$

(37)

2.2.2. Elastic Modulus in the X-Direction $E_x$ and Poisson’s Ratio ${\nu }_{xz}$

To determine the equivalent elastic modulus in the x-direction $E_x$, it is necessary to determine the horizontal displacement ${\delta }_H$ of the upper end due to a horizontal force $H$. The calculation approach is like $G_{xz}$, but there is no dummy vertical force $V$. Therefore, according to Equation (33) and considering the conditions of ${\delta }_{M_0}=0$, $H=1$ and $V=0$ at point B are finally imposed. The arising system of three equations with three unknowns can be solved to determine ${\delta }_H$.

$${\delta }_H=C_{1,1}-\frac{C^2{}_{1,3}}{C_{3,3}}$$

(38)

$${\delta }_V=C_{2,1}-\frac{C_{2,3}{C}_{3,1}}{C_{3,3}}$$

(39)

Additionally, the equivalent Young’s modulus $E_{x2}$ of Part A–B is derived as follows:

$$E_{x2}=\frac{{\sigma }_x}{{\epsilon }_x}=\frac{F_x}{A_{yz}}/\frac{{\delta }_x}{l_x}=\frac{1}{h}/\frac{{\delta }_H}{p_2}$$

(40)

Hence, the equivalent elastic modulus in the x-direction $E_x$ of the whole trapezoidal-shaped corrugation is

$$E_x=p/\left(\frac{p_1}{E_{x1}}+\frac{p_2}{E_{x2}}+\frac{p_3}{E_{x3}}\right)$$

(41)

where $E_{x1}$ and $E_{x3}$ represent the equivalent elastic modulus of Part O–A and Part B–C, respectively, and

$$E_{x1}=E_{x2}=E\times t_c/h$$

(42)

Furthermore, Poisson’s ratio ${\nu }_{xz}$ can be derived as

$${\nu }_{xz}=-\frac{{\epsilon }_z}{{\epsilon }_x}=\frac{{\delta }_V}{h}/\frac{({\delta }_H+\frac{p_1+p_3}{t_{c}E})}{p}$$

(43)

2.2.3. Elastic Modulus in the Y-Direction $E_y$

The equivalent elastic modulus $E_y$ in the y direction can be calculated by Equation (44), which scales the Young’s modulus of the core $E$ to the ratio between the sectional area of corrugated core in X–Z plane, $A_y$, and the section area of the equivalent volume, $A_{eqy}$.

$$E_y=E\frac{A_y}{A_{eqy}}=E\frac{t_{c}l}{hp}$$

(44)

$$l=p_1+\frac{p_2}{\cos \beta }+p_3$$

(45)

2.2.4. Transverse Shear Modulus in Y–Z Plane $G_{yz}$

The transverse shear modulus $G_{yz}$ can be obtained by simplifying the shape of trapezoidal corrugated core. In the sample domain of half a period, the core is stretched to a flattened panel. This corresponds to creating a local curved system of reference along the center line of the core sheet and integrating along that local coordinate. Then, the panel can be considered as a straight beam with a rectangular section. Based on Timoshenko beam theory, the displacement in the y direction due to an applied force $H$ in the same direction is derived as

$${\delta }_{zy}^{re}=\frac{H{l}^{\prime }}{\frac{5}{6}\times G\times t_{c}b}$$

(46)

where $l^{\prime }=\frac{P_2}{\cos \beta }$ and b is the width of the core in the y direction. The shear deformation of the equivalent sample section, ${\delta }_{zy}^{eq}$, under the same conditions, is

$${\delta }_{zy}^{eq}=\frac{Hh}{\frac{5}{6}\times G_{yz2}{p}_{2}b}$$

(47)

By imposing that ${\delta }_{zy}^{re}$ and ${\delta }_{zy}^{eq}$ are equal, the equivalent shear modulus $G_{yz2}$ of Part A–B is obtained as

$$G_{yz2}=G\frac{h{t}_c}{p_2{l}^{\prime }}$$

(48)

Then, the transverse shear modulus in the Y–Z plane of the trapezoidal corrugation is approximated as:

$$G_{yz}=G_{yz1}\frac{p_1}{p}+G_{yz2}\frac{p_2}{p}+G_{yz3}\frac{p_3}{p}\approx G_{yz2}\frac{p_2}{p}$$

(49)

2.2.5. In-Plane Shear Modulus $G_{xy}$

The derivation method of in-plane shear modulus $G_{xy}$ is similar to the derivation method of transverse shear modulus in the Y–Z plane $G_{yz}$.The only difference is that now the force $H^{\prime }$ is acting in the x direction. Therefore, the shear deformations computed for the approximated section ${\delta }_{xy}^{re}$ and for the equivalent sample section ${\delta }_{xy}^{eq}$ are:

$${\delta }_{xy}^{re}=\frac{H^{\prime }b}{\frac{5}{6}\times G{t}_{c}l}$$

(50)

$${\delta }_{xy}^{eq}=\frac{H^{\prime }b}{\frac{5}{6}\times G_{xy}hp}$$

(51)

By imposing that ${\delta }_{xy}^{re}$ and ${\delta }_{xy}^{eq}$ are equal, the equivalent shear modulus $G_{xy}$ is obtained as

$$G_{xy}=G_{al}\frac{t_{c}l}{hp}$$

(52)

2.2.6. In-Plane Poisson’s Ratio ${\nu }_{xy}$

Poisson’s ratio ${\nu }_{xy}$ of typical orthotropic materials can be obtained from the following relation:

$${\nu }_{xy}{E}_y={\nu }_{yx}{E}_x$$

(53)

The Poisson’s ratio in the X–Y plane ${\nu }_{yx}$ can be assumed to be equal to the constituent material Poisson’s ratio [22].

2.2.7. Out-of-Plane Elastic Modulus in Z-Direction $E_z$

As represented in Figure 2, the derivation of elastic modulus in the z direction is similar to that of the elastic modulus in the x direction. By imposing a vertical force $V$ in the z direction to the point B, positive vertical displacement, ${\delta }_V$, and negative horizontal displacement in the x-direction, ${\delta }_H$, occur. Additionally, in order to avoid rotation in the X–Z plane, a dummy moment $M_0$ will be produced at point B. Therefore, according to Equation (33) and considering the conditions of ${\delta }_{M_0}=0$, $H=0$, and $V=1$, the arising system of three equations with three unknowns can be solved to determine ${\delta }_V$ and ${\delta }_H$:

$${\delta }_V=C_{2,2}-\frac{C^2{}_{2,3}}{C_{3,3}}$$

(54)

$${\delta }_H=C_{1,2}-\frac{C_{1,3}{C}_{2,3}}{C_{3,3}}$$

(55)

Additionally, the equivalent Young’s modulus $E_z$ is obtained as:

$$E_z=\frac{{\sigma }_z}{{\epsilon }_z}=\frac{F_z}{A_{xy}}/\frac{{\delta }_V}{l_z}=\frac{1}{p}/\frac{{\delta }_V}{h}$$

(56)

2.2.8. Poisson’s Ratio ${\nu }_{zy}$

The Poisson’s ratio, ${\nu }_{zy}$, is obtained from the following relation:

$${\nu }_{zy}{E}_y={\nu }_{yz}{E}_z$$

(57)

Similar to the calculation of ${\nu }_{yx}$, the Poisson’s ratio of ${\nu }_{yz}$ can also be assumed to be equal to the constituent material of Poisson’s ratio. Additionally, then, ${\nu }_{zy}$ can be calculated by Equation (57).

2.2.9. Relative Density ${\rho }_{\mathrm{core}}$

The relative density of the equivalent material is computed by scaling the density of the constituent material $\rho$ proportionally to the occupied equivalent volume, thus

$${\rho }_{\mathrm{core}}=\rho \frac{t_{c}l}{hp}$$

(58)

Generally, the present theoretical results are almost consistent with the calculated results of the FEA shell elements models with the same boundary conditions [21,22,23,24,25]. By comparing the results of the FEA model and Xia’s method [20], Bartolozzi et al. [22] proved that the proposed method can still be used with good accuracy for trapezoidal corrugated cores with relatively short horizontal segments. Furthermore, compared with Xia’s method, the proposed method can further calculate the out-of-plane equivalent properties.

2.3. The Theoretical Equivalent Homogenization Method for a Sandwich Cylinder

Composite sandwich tubes are widely used in many fields such as civil, aerospace, and military. In the civil field, high-pressure long-tube trailers are an important way to transport hydrogen, compressed natural gas, and other fuels in a short distance. Compared with the traditional steel container, the composite cylinder has higher efficiency of storage and transportation, but, as a result of increases in composite cylinder volume, a longer bottle body would lower the stiffness of the vessel. Therefore, large-volume full-winding hydrogen bottles should not only meet the requirements of their internal pressure strength but should also still meet the tube body stiffness requirements in service conditions. In the field of aeronautics and astronautics, a lightweight composite allio sandwich tube structure can meet the requirements of the bearing performance of large engine shells under complex working conditions, and the dialectical unity of its bending stiffness, axial compression performance, and lightweight is continuously pursued in the design process. In the military field, the improvement of the bending stiffness of the lightweight composite rocket launcher can effectively reduce the lateral vibration during launch, reduce the velocity deviation of the projectile leaving the launcher, and improve the hit accuracy. Therefore, the bending stiffness of composite sandwich tubes is a significant performance index. It is urgent to establish an effective analysis method to design the bending stiffness of this kind of structure.

An equivalent homogenization theoretical method is proposed to analyze the pure bending of a sandwich cylinder. As shown in Figure 6a, the sandwich cylinder is composed of the inner skin, lightweight core, and outer skin. The geometric parameters of the trapezoidal corrugated core are shown in Figure 6b. According to the equivalent theory model of trapezoidal corrugated core in Section 3, the core could be approximated to an orthotropic layer with equivalent mechanical properties. As shown in Figure 6c, the core could be simplified into an orthotropic homogeneous layer with the same height as h_c. By using the developed Lekhnitskii’s theory, the equivalent bending stiffness of the corrugated-core sandwich cylinder can be calculated.

3. Results

To compare with the proposed equivalent homogenization theoretical method, a 3D finite element model of the sandwich cylinder was established and calculated by ABAQUS 6.13. As shown in Figure 6b, the geometric parameters of inner skin, outer skin, and the trapezoidal corrugated core are listed in

Table 1. The geometric parameters of the sandwich cylinder.

Parameters	Dimensions
Inner radius, $r_o$	95 mm
Inner skin thickness, tfi	1.6 mm
Winding angles of inner skin, $\phi$	[90/0]
Outer skin thickness, tfo	1.6 mm
Winding angles of outer skin, $\phi$	[90/0]
Cell number	50
Height of cores, hc	6.8 mm
Angle of corrugation, θ	53.57°
Corrugated rib thickness, tc	1 mm
Length of upper horizontal part, 2 × p1	2 mm
Length of below horizontal part, 2 × p3	2 mm
Length of trapezoidal core, 2 × p	12.56 mm

. The inner and outer skin have two layers with the same thickness of 0.8 mm. The material of the core is standard aluminum with $E$ = 71,000 MPa and $\nu$ = 0.33. The material of the inner and outer skin is carbon fiber/epoxy resin with anisotropic elastic properties listed in

Table 2. Material Properties.

Material Properties	Carbon Fiber/Epoxy Resin	Glass Fiber/Epoxy Resin
E11/MPa	1.55 × 105	6.0 × 104
E22 = E33/MPa	1.21 × 104	1.1 × 104
G23/MPa	3.2 × 103	3.2 × 103
G13 = G12/MPa	4.4 × 103	7.6 × 103
v12 = v13	0.248	0.26
v23	0.458	0.458

.

By developing the user subroutine of UMAT, the orthotropic linear elastic constitutive model is applied in the ABAQUS Implicit Solver. Considering that many meshes require a considerable amount of computing resources, the semi-model of the sandwich cylinder was established. As shown in Figure 7, there are 1000 regular elements in the circumferential direction, 2 regular elements of each layer in the thickness direction, and 200 regular elements in the axial direction. The linearly distributed load with the equivalent moment of $M_x$ = 6 × 10⁶ N·mm and $M_y$ = 0 N·mm was applied on the semi-model. Additionally, the Y–Z section of the model is symmetrically constrained.

As shown in Figure 8, under the equivalent moment, $M_x$ = 6 × 10⁶ N·mm, the axial stress component at the position of θ = 90° is greater than that at other positions. In order to measure the bending curvature of the cylinder, we chose a length of the cylinder to extract the deformed coordinates of the nodes at the position of θ = 90° and 270°. Additionally, the coordinates of the deformed central axis can be calculated, as shown in Figure 8. Because the deformation of the cylinder is symmetrical to the x-axis, only the coordinates of the central axis in the Y–Z plane need to be considered. It is known that every three adjacent nodes in Y–Z plane can determine a circle with a radius, R. The cylinder’s bending radius can be determined by averaging the corresponding radius of all nodes in the central axis.

According to the proposed equivalent homogenization theoretical method, the trapezoidal corrugated core with geometric parameters listed in Table 1 needs to be equivalent to an orthotropic homogenization layer. The mechanical properties are shown in

Table 3. The equivalent parameters for trapezoidal profile geometries.

Ex (MPa)	Ey (MPa)	vyx	Gxy (MPa)	Ez (MPa)	Gxz (MPa)	Gyz (MPa)	vyz	vxz
295.74	17,948.8	0.33	6747.68	453.30	3314.44	3420.03	0.33	0.771

.

Then, the curvature of the equivalent cylinder can be calculated by Equation (28). The numerical and theoretical results are shown in

Table 4. The results of the equivalent homogenization theoretical model and the FE model.

Results	${\mathit{k}}_{\mathit{x}}/{\mathbf{mm}}^{-1}$	${\mathit{\sigma }}_{\mathit{z}}$/MPa $\left(\mathit{\theta }=90°,\mathit{r}={\mathit{r}}_0\right)$	${\mathit{\sigma }}_{\mathit{z}}$/MPa $\left(\mathit{\theta }={90}^{\mathbf{o}},\mathit{r}={\mathit{r}}_0+{\mathit{t}}_{\mathbf{fi}}+{\mathit{h}}_{\mathit{c}}+{\mathit{t}}_{\mathbf{fo}}\right)$
FE model	9.356 × 10−6	11.647	166.4
Theoretical method	9.672 × 10−6	11.118	158.0
Relative errors	−3.34%	4.54%	5.05%

, which shows that the results of the equivalent method were in good agreement with numerical results.

4. Discussion

4.1. The Bending Strength of Sandwich Cylinder

The advantage of the theoretical method is that the core layer is simplified as a homogenized layer, making the sandwich cylinder’s three-dimensional elastic analysis possible.

Intuitively, the bending stiffness of the cylinder can be improved by reducing the winding angle of the inner skin and outer skin. The calculation of the bending strength can reflect the utilization rate of the materials in each layer of the sandwich cylinder. Therefore, the bending strength of sandwich cylinder is also important.

4.1.1. Strength Criterion for the Cylinder’s Skin

The inner and outer skins of the cylinder are composed of winding layers with different materials, angles, and thicknesses. The failure of skins belongs to the progressive failure of each layer. Therefore, it is of great significance to find out the earliest failure location and the failure load. In recent years, many strength criteria [37] have been developed for fiber-reinforced resin matrix composites. In this paper, the Hill–Tsai strength criterion is used to evaluate the strength failure of the winding layers of skin. Additionally, the failure determination formula is

$$\begin{array}{ll}{f}_{\mathrm{s}}=& \frac{{\sigma }_1^2}{X_1^2}+\frac{{\sigma }_2^2}{Y_1^2}+\frac{{\sigma }_3^2}{Z_1^2}-\left(\frac{1}{X_2^2}+\frac{1}{Y_2^2}-\frac{1}{Z_2^2}\right){\sigma }_1{\sigma }_2\\ & -\left(\frac{1}{X_3^2}+\frac{1}{Y_3^2}-\frac{1}{Z_3^2}\right){\sigma }_1{\sigma }_3\\ & -\left(\frac{1}{X_4^2}+\frac{1}{Y_4^2}-\frac{1}{Z_4^2}\right){\sigma }_2{\sigma }_3+\frac{{\tau }_{12}^2}{S_{12}^2}+\frac{{\tau }_{23}^2}{S_{23}^2}+\frac{{\tau }_{31}^2}{S_{31}^2}\end{array}$$

(59)

where $f_{\mathrm{s}}$ is the failure factor; ${\sigma }_1$, ${\sigma }_2$, ${\sigma }_3$, ${\tau }_{12}$, ${\tau }_{13}$, and ${\tau }_{23}$ are the stress components in the material principal coordinate system of 1-2-3 directions; and $S_{12}$, $S_{13}$, and $S_{23}$ are the shear strength in the corresponding principal axis directions.

If ${\sigma }_1>0$, then $X_1=X_T$; otherwise $X_1=X_C$.

If ${\sigma }_2>0$, then $X_2=X_T$, $Y_1=Y_2=Y_T$ and $Z_2=Z_T$; otherwise $X_2=X_C$, $Y_1=Y_2=Y_C$, $Z_2=Z_C$.

If ${\sigma }_3>0$, then $X_3=X_T$, $Y_3=Y_T$ and $Z_1=Z_3=Z_T$; otherwise $X_3=X_C$, $Y_3=Y_C$ and $Z_1=Z_3=Z_C$.

If ${\sigma }_2<0$ or ${\sigma }_3<0$, then $X_4=X_C$, $Y_4=Y_C$ and $Z_4=Z_C$; otherwise $X_4=X_T$, $Y_4=Y_T$ and $Z_4=Z_T$.

where $X_{\mathrm{c}}$, $Y_{\mathrm{c}}$, and $Z_{\mathrm{c}}$ represent compressive strength and $X_{\mathrm{t}}$, $Y_{\mathrm{t}}$ and $Z_{\mathrm{t}}$ represent tensile strength in the principal directions. When the stress components of a material point make $f_{\mathrm{s}}>1$, the material point is considered to be damaged. According to the explicit stress expression of the developed Lekhnitskii’s theory mentioned in Section 2.3, it is possible to predict the failure location. As shown in Equation (44), the stress components need to be transformed into corresponding stress components in material principal coordinate system.

$$\left\{\begin{array}{l}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\tau }_{23}\\ {\tau }_{13}\\ {\tau }_{12}\end{array}\right\}=\mathrm{T}\left\{\begin{array}{l}{\sigma }_r\\ {\sigma }_{\theta }\\ {\sigma }_z\\ {\tau }_{\theta z}\\ {\tau }_{rz}\\ {\tau }_{r\theta }\end{array}\right\}$$

(60)

4.1.2. Strength Criterion for the Core

For a sandwich cylinder, the neutral layer is coplanar with the axis, which will cause a complex stress distribution along with the circumferential and radial directions. Therefore, the actual stress state of the core with multi-corrugation shapes is more complicated. Although the stress of each layer is not necessarily continuous, the strain along the circumferential and radial direction must be continuous. For a trapezoidal corrugated core, the axial stiffness is much greater than that in other directions. Additionally, the main strain is in the axial direction of the cylinder, which causes that the ${\sigma }_z$ to be much greater than other stress components at each layer, as shown in Figure 2. When the core’s axial displacement reaches the maximum uniaxial strain of the core material, failure is considered to have occurred in the core. Therefore, the maximum strain criterion was applied in this paper to evaluate the failure of cores, which is expressed as

$$\{\begin{array}{l}{f}_{\mathrm{c}}={\left(\frac{{\epsilon }_{z,\max }^{core}}{{\epsilon }_{\mathrm{f}}^{T\left(C\right)}}\right)}^2\\ {\epsilon }_{\mathrm{f}}^{T\left(C\right)}=\frac{X_{T\left(C\right)}}{E_{\mathrm{core}}}\end{array}$$

(61)

where ${\epsilon }_{z,\max }^{\mathrm{core}}$ is the maximum strain component in the z-direction of the core layer, ${\epsilon }_{\mathrm{f}}^{T\left(C\right)}$ is the critical failure strain of the core material, and $X_{T\left(C\right)}$ is the strength of the core material. Generally, the corrugated cores are made of lightweight metal or woven carbon fiber composite with symmetric orthogonal ply. For metal material, $X_{T\left(C\right)}$ is the yield strength and $E_{\mathrm{core}}$ is the elastic modulus.

When the core is made of woven carbon fiber composite with ply order of [0/90]_ns, $X_T$, $X_C$, and $E_{\mathrm{core}}$ are the equivalent tensile strength, compressive strength, and elastic modulus of the [0/90]_ns laminated plate, respectively. According to classical laminated plate theory, the relationship between the internal force and the deformation in the middle plane of a laminate is established as

$$\left[\begin{array}{l}{\epsilon }_x\\ {\epsilon }_y\\ {\epsilon }_{xy}\\ k_x\\ k_y\\ k_z\end{array}\right]=\left[\begin{array}{l}{A}^{\prime } B^{\prime }\\ B^{\prime } D^{\prime }\end{array}\right]\left[\begin{array}{l}{N}_x\\ N_y\\ N_{xy}\\ M_x\\ M_y\\ M_{xy}\end{array}\right]$$

(62)

where matrix elements $A^{\prime }$, $B^{\prime }$, and $D^{\prime }$ are the tensile, coupling, and bending flexibility, respectively. For the laminated plate with symmetric ply, the coupling flexibility matrix is a null matrix. The relationship above can be simplified into

$$\left[\begin{array}{l}{\epsilon }_x\\ {\epsilon }_y\\ {\epsilon }_{xy}\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& 0\\ a_{12}& a_{22}& 0\\ 0& 0& a_{66}\end{array}\right]\left[\begin{array}{l}{N}_x\\ N_y\\ N_{xy}\end{array}\right]=h\left[\begin{array}{ccc}{a}_{11}& a_{12}& 0\\ a_{12}& a_{22}& 0\\ 0& 0& a_{66}\end{array}\right]\left[\begin{array}{l}\overline{{\sigma }_x}\\ \overline{{\sigma }_y}\\ \overline{{\sigma }_{xy}}\end{array}\right]$$

(63)

where h is the thickness of the laminated. The effective in-plane engineering constant can be presented as

$$\{\begin{array}{l}{E}_x=\frac{1}{h{a}_{11}},{\nu }_{xy}=\frac{-a_{12}}{a_{11}}\\ E_y=\frac{1}{h{a}_{22}},G_{yx}=\frac{1}{h{a}_{66}}\end{array}$$

(64)

When the layup configuration of the core is orthogonal and symmetrical, the laminate has an in-plane quasi-isotropic characteristic. Therefore, the equivalent theory model of trapezoidal corrugated cores is also applicable to the above composite cores. The in-plane elastic modulus of the core can be determined as $E_{\mathrm{core}}=E_x=E_y$.

4.2. Parametric Study for Specified Bending Stiffness and Specified Bending Strength

There are many parameters that affect the mechanical properties of sandwich cylinders. Based on the geometric parameters listed in Table 4, the specified bending stiffness (the ratio of bending stiffness to structure density) and specified bending strength (the ratio of bending strength to structure density) were studied under variable parameters of skin thickness, winding angles of skin, corrugated rib thickness, and core materials. The material properties of winding layers and cores are listed in

Table 5. Material property for skins and cores.

Winding Layers (Skins)	E1 = 158 GPa, E2 = E3 = 14.4 GPa; u12 = u13 = 0.413; u23 = 0.5; G12 = G13 = 4.4 GPa; G23 = 3.7 GPa; XT = 1850 Mpa; XC = 1470 MPa; YT = 29 MPa; YC = 140 MPa; S12 = S13 = S23 = 65 Mpa; ρ = 1.7 g/cm3
Cores	Standard aluminum alloy	E = 71 GPa, $\mathit{\nu }$ = 0.33, X = ${\sigma }_s$ = 240 MPa, ρ = 2.8 g/cm3
	Carbon fiber woven composite with [0/90]ns	E = 54.5 GPa, $\nu$ = 0.06, G = 4400 MPa XT = XC = 473 MPa, ρ = 1.7 g/cm3

. When the material of the core is woven carbon fiber composite, the effective in-plane engineering constant can be calculated according to classical laminated plate theory. However, it is more convenient to test the material property of a [0/90]_ns laminated plate [11] directly. The thickness of each winding layer is 0.2 mm. In the process of the parametric study, the intermediate diameter of the core was kept unchanged.

According to the density of each layer listed in Table 5 and the core density, ${\rho }_{\mathrm{core}}$, calculated by Equation (58), the equivalent density of the sandwich cylinder can be obtained. As shown in Figure 9, the influence of core material on standard aluminum alloy (SAA) and carbon fiber woven composite (CFWC) with [0/90]_ns was studied and contrasted. The winding angles and thickness of inner and outer skin remain the same. With the increase in winding angle, the specified bending stiffness of the sandwich cylinder drops rapidly at winding angle interval of 0° to ±45°. It keeps at an almost constant level at the interval of ±45° to ±90°. This is because the trapezoidal corrugated core is the main contributor to the specified bending stiffness at winding angle interval of ±45° to ±90°. The cores made of woven carbon fiber composite with lower in-plane elastic modulus show a higher specified bending stiffness due to its lightweight materials characteristic. Further, the specified bending strength with the same parameters is shown in Figure 10. For aluminum alloy cores, the unique failure mode is its yield failure. When the core becomes the main contributor to the specified bending stiffness at winding angle intervals of ±45° to ±90°, the specified bending strength no longer changes significantly. The sandwich cylinder with the composite core has different failure modes. When the skin winding angels are ±15°, ±30°, and ±45°, the failure mode is skin failure, which is mainly caused by the resin matrix tension at the positions of $\theta ={270}^{\mathrm{o}}$ and $r=r_0$. Therefore, due to the transition of the failure modes, the bending strength reaches its minimum at the winding angle of ±30°. The other failure mode is core failure at winding angels of 0°, ±60°, ±75°, and 90°. The phenomenon above proves that the multi-winding angles are necessary to improve the skin stiffness of other directions and prevent the matrix failure from tensile load in advance.

In order to reveal the influence of core parameters on specified bending stiffness and strength further, the intermediate diameter and the material of cores remain the same. The thickness and height of the corrugated cores were considered. As shown in Figure 11 and Figure 12, the smaller thickness of corrugated ribs and larger height of corrugated cores can improve the specified bending stiffness and strength of the sandwich cylinder. For variable core parameters, the failure modes are all skin failure, which is mainly caused by the matrix tension. Due to the effect of mass density, when the skin thickness is 1.5 mm, there is an inflection point in the specific bending strength curve shown in Figure 11. However, the influence of core parameters on the specified bending stiffness and strength are lower than that of the skin parameters. This is because the mass density of cores is greatly affected by the core parameters, which could balance the variation of stiffness and strength of the corrugated core sandwich cylinder.

5. Conclusions

A theoretically equivalent homogenization method for composite sandwich cylinders under pure bending was proposed to overcome the bending stiffness and strength calculation, which has good accuracy with the FE model.

Based on a homogeneous orthotropic layer hypothesis, the trapezoidal corrugated sandwich core was equivalent to a material with nine equivalent mechanical properties. Lekhnitskii’s theory based on a unified connection parameter method was introduced and applied in the equivalent composite sandwich cylinder. A case study shows that the developed Lekhnitskii theory is suitable for an arbitrary combination of winding layers with different angles or materials. The bending stiffness of the equivalent sandwich cylinder could be calculated.

By developing the user subroutine of UMAT, the numerical calculation results were in good agreement with the results of the proposed method. Further, according to the Hill–Tsai strength criterion and the maximum strain criterion, parametric study for specified bending stiffness and specified bending strength was done. The results show that the influence of core parameters on the specified bending stiffness and strength is lower than that of the skin parameters. Larger skin thickness and smaller winding angles could improve the specified bending stiffness and strength of the composite corrugated sandwich cylinders.