**1. Introduction**

In recent years, composite structures are more and more widely used as load-bearing structures [1–4]. Three-dimensional (3D) textile composites is a new type of high performance composites developed from the traditional two-dimensional (2D) textile composites in the 1980s. Compared with 2D textile composites, the warp yarns, weft yarns and binder yarns are interlaced with each other not only in the plane, but also in the thickness direction, which can not only improve the specific strength and specific stiffness of composites, but also have other excellent mechanical properties, such as good impact damage resistance, fatigue resistance, etc. [5]. At the same time, it also overcomes the shortcomings of laminated composite that are easy to delaminate after loading. Another significant benefit of 3D textile composites is the ability to manufacture structural component reforms directly from the yarns.

The application of 3D textile composites in load-bearing structures requires a thorough structural analysis. However, precisely predicting the behaviors of textile composites is difficult due to their complex microstructures. Currently, experimental [6], analytical and numerical methods are used to study the elastic behavior of 3D textile composites.

A variety of analytical models have been used to study the mechanical characteristics of 3D braided composites, including fabric geometry model [7], fiber inclination model [8], three cell model [9], mixed volume averaging technique [10], and Mori-Tanaka theories

**Citation:** Xi, S.; Zhong, Y.; Shi, Z.; Yi, Q. Prediction of Bending, Buckling and Free-Vibration Behaviors of 3D Textile Composite Plates by Using VAM-Based Equivalent Model. *Materials* **2022**, *15*, 134. https:// doi.org/10.3390/ma15010134

Academic Editor: Krzysztof Schabowicz

Received: 19 November 2021 Accepted: 21 December 2021 Published: 24 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

combined with stiffness averaging method [11]. Yang et al. [12] introduced the "Fiber Inclination Model" to estimate the elastic properties of 3D textile (woven and braided) composites. The unit cell employed for the study was an assemblage of inclined unidirectional laminae. Based on the classical laminate plate theory (CLPT) and iso-stress/strain assumptions, Ishikawa and Chou [13,14] presented analytical methods to estimate the homogenized response of woven fabric composites. Branch et al. [15] presented a 3D tow inclination model for calculating the elastic constants of three-dimensional braided composites. The global constitutive equation of the composite material was derived by applying an iso-strain method to the unit cell and averaging all tow segments and the matrix inside the unit cell. Yan et al. [16] predicted the properties of 3D braided structures using an analytical model called the Fabric Geometry Model (FGM). Using the stiffness volume average method and Tsai-Wu polynomial failure criterion, Jiang et al. [17] presented a theoretical model based on the helix geometry unit cell for prediction of the effective elastic constants and the failure strength of 3D braided composites under uniaxial load. Analytical models are good at estimating the in-plane properties of textile composites, but they are not so good at predicting the shear and out-of-plane properties.

Although utilizing experimental and analytical models to examine the mechanical characteristics of a laminated composite plate is practical and efficient [18], many studies also employed numerical approaches to investigate the 3D textile composites [19–21]. Tan et al. [22] developed a mesoscale finite element model (MSFEM) in LS-DYNA to simulate impact damage to three-dimensional braided composite plates based on the assumption that the fiber yarn was made up of cylindrical segments. Dong et al. [23] simulated the micro-stress of 3D braided composites by the method of Asymptotic Expansion Homogenization (AEH) combined with finite-element analysis. Tang and Whitcomb [24] used the full multiscale mechanical model to perform progressive failure evaluations of 2 × 2 braided composites. This approach was also used by Potluri and Manan to investigate the mechanical characteristics of braided composite tubes [25], but the stress and strain distribution of the fiber and matrix cannot be determined.

Berdichevsky [26] recently developed the semi-analytical approach-variational asymptotic technique (VAM) to increase the efficiency of numerical methods and the accuracy of analytical methods. It combines the benefits of asymptotic and variational methods, and takes into account all potential deformations [27–29]. The fundamental benefit of adopting VAM for plate analysis is that it can divide the original 3D plate problem into two independent problems using the small parameter of thickness-width ratio, i.e., throughthe-thickness analysis and 2D reference plane analysis [30,31]. Then, VAM was expanded by Zhong and Yu to simulate piezoelectric and piezomagnetic laminates [32], multilayer graded magnetoelectroelastic plates [33] and composite cylindrical shells [34].

The VAM model is expanded in this article to present an equivalent plate model to replace the original 3D textile composite plate (3D-TCP) for bending, buckling, and freevibration analysis. The influences of structural parameters (binder yarn width, warp or weft yarn width) on the equivalent stiffness of 3D-TCP are investigated. Finally, the effective performance of 2D plain-woven laminate (2D-PWL) and 3D-TCP with the same plate thickness and yarn content are compared. To the best of the authors' knowledge, this technique has never been used to predict the bending, buckling, and free-vibration behaviors of 3D-TCP.

#### **2. Variational Asymptotic Equivalent Model of 3D-TCP**

As illustrated in Figure 1a, the macro-coordinates *xi*(*i* = 1, 2, 3) may be used to describe any point in the 3D-TCP, where *xα*(*α* = 1, 2) are the in-plane coordinates and *x*<sup>3</sup> is the normal coordinate. The 3D finite element model of 3D-TCP may be divided into the 3D unit cell and 2D equivalent plat model (2D-EPM) according to the VAM. It's worth mentioning that the dimensions of the unit cell should be substantially smaller than macro-structure dimensions. As illustrated in Figure 1c, the field variables of 2D-EPM are represented as functions of *x*<sup>1</sup> and *x*2, whereas *x*<sup>3</sup> is disappeared.

**Figure 1.** Decomposition diagram of 3D textile composite plate: (**a**) 3D-FEM; (**b**) Unit cell; (**c**) 2D-EPM.

To characterize the quick change of in-plane material characteristics, the micro-coordinates *yi* = *xi*/*ζ* (*ζ* is a small parameter) are introduced. To obtain the equivalent model of 3D-TCP by using the variational asymptotic method, the 3D displacement field of the original 3D-TCP needs to be represented by 2D plate variables, such as

$$\begin{array}{l} \mu\_{1}(\mathbf{x}\_{a};y\_{i}) = \overline{\eta\_{1}(\mathbf{x}\_{a})} - \zeta y\_{3}\overline{\eta\_{3,1}(\mathbf{x}\_{a})} + \zeta w\_{1}(\mathbf{x}\_{a};y\_{i})\\ \mu\_{2}(\mathbf{x}\_{a};y\_{i}) = \overline{\eta\_{2}(\mathbf{x}\_{a})} - \zeta y\_{3}\overline{\eta\_{3,2}(\mathbf{x}\_{a})} + \zeta w\_{2}(\mathbf{x}\_{a};y\_{i})\\ \mu\_{3}(\mathbf{x}\_{a};y\_{i}) = \overline{\eta\_{3}(\mathbf{x}\_{a})} + \zeta w\_{3}(\mathbf{x}\_{a};y\_{i}) \end{array} \tag{1}$$

where the displacements of the 3D-FEM and 2D-EPM are represented by *ui* and *u*¯*i*, respectively; *wi* are the fluctuation functions to be solved. The underline terms should satisfy the following constraints:

$$h\ddot{u}\_a(\mathbf{x}\_a) = \langle u\_a \rangle + \langle \zeta y\_3 \rangle \ddot{u}\_{3,a}, \quad h\ddot{u}\_3(\mathbf{x}\_a) = \langle u\_3 \rangle \tag{2}$$

where · denotes the volume integral of a unit cell.

The fluctuation functions in Equation (3) are constrained as

$$
\langle \zeta w\_i \rangle = 0 \tag{3}
$$

The strain field can be expressed according to 3D linear elasticity theory, such as

$$
\varepsilon\_{ij} = \frac{1}{2} \left( \frac{\partial u\_i}{\partial x\_j} + \frac{\partial u\_j}{\partial x\_i} \right) \tag{4}
$$

The 3D strain field can be obtained by substituting Equation (1) into Equation (4) and ignoring higher-order terms according to VAM,

$$\begin{aligned} \varepsilon\_{11} &= \gamma\_{11} + \tilde{\zeta}y\_3 \kappa\_{11} + w\_{1,1} \\ 2\varepsilon\_{12} &= 2\gamma\_{12} + 2\tilde{\zeta}y\_3 \kappa\_{12} + w\_{1,2} + w\_{2,1} \\ \varepsilon\_{22} &= \gamma\_{22} + \tilde{\zeta}y\_3 \kappa\_{22} + w\_{2,2} \\ 2\varepsilon\_{13} &= w\_{1,3} + w\_{3,1} \\ 2\varepsilon\_{23} &= w\_{2,3} + w\_{3,2} \\ \varepsilon\_{33} &= w\_{3,3} \end{aligned} \tag{5}$$

where *γαβ* and *καβ* are the in-plane strains and bending curvatures of 2D-EPM, respectively, and may be defined as

$$
\pi\_{a\beta}(\mathbf{x}\_1, \mathbf{x}\_2) = \frac{1}{2} (\bar{u}\_{a,\beta} + \bar{u}\_{\beta,a}), \quad \kappa\_{a\beta}(\mathbf{x}\_1, \mathbf{x}\_2) = -\bar{u}\_{3,a\beta} \tag{6}
$$

We can define the three-dimensional strain field in matrix form to make derivation easier, such as

$$\begin{aligned} \mathcal{E}\_{\mathfrak{k}} &= \begin{bmatrix} \varepsilon\_{11} & \varepsilon\_{22} & 2\varepsilon\_{12} \end{bmatrix}^{\mathrm{T}} = \gamma + \mathfrak{x}\_{3}\mathfrak{k} + I\_{\mathfrak{a}}\mathfrak{w}\_{||,\mathfrak{a}} \\ \mathfrak{2}\mathcal{E}\_{\mathfrak{k}} &= \begin{bmatrix} 2\varepsilon\_{13} & 2\varepsilon\_{23} \end{bmatrix}^{\mathrm{T}} = \mathfrak{w}\_{||,3} + \mathfrak{e}\_{\mathfrak{a}}\mathfrak{w}\_{3,\mathfrak{a}} \\ \mathcal{E}\_{\mathfrak{k}} &= \varepsilon\_{33} = \mathfrak{w}\_{3,3} \end{aligned} \tag{7}$$

where ()|| = [()<sup>1</sup> ()2] T, *<sup>γ</sup>* <sup>=</sup> <sup>0</sup> *<sup>γ</sup>*<sup>11</sup> <sup>2</sup>*γ*<sup>12</sup> *<sup>γ</sup>*<sup>22</sup> <sup>1</sup><sup>T</sup> , *κ* = [*κ*<sup>11</sup> *κ*<sup>12</sup> + *κ*<sup>21</sup> *κ*22] T, and

$$I\_1 = \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ 0 & 0 \end{bmatrix}, I\_2 = \begin{bmatrix} 0 & 0 \\ 1 & 0 \\ 0 & 1 \end{bmatrix}, \mathbf{e}\_1 = \begin{Bmatrix} 1 \\ 0 \end{Bmatrix}, \mathbf{e}\_2 = \begin{Bmatrix} 0 \\ 1 \end{Bmatrix} \tag{8}$$

The strain energy of the 3D-TCP may be expressed briefly as

$$\mathcal{U} = \frac{1}{2} \left\langle \boldsymbol{\mathcal{E}}^{\mathrm{T}} \mathbf{D} \boldsymbol{\mathcal{E}} \right\rangle = \frac{1}{2} \left\langle \begin{Bmatrix} \boldsymbol{\mathcal{E}}\_{\varepsilon} \\ \boldsymbol{2\mathcal{E}}\_{s} \\ \boldsymbol{\mathcal{E}}\_{t} \end{Bmatrix}^{\mathrm{T}} \begin{bmatrix} \mathbf{D}\_{\varepsilon} & \mathbf{D}\_{\varepsilon s} & \mathbf{D}\_{\varepsilon t} \\ \mathbf{D}\_{\varepsilon s}^{\mathrm{T}} & \mathbf{D}\_{s} & \mathbf{D}\_{\varepsilon t} \\ \mathbf{D}\_{\varepsilon t}^{\mathrm{T}} & \mathbf{D}\_{\varepsilon t}^{\mathrm{T}} & \mathbf{D}\_{t} \end{bmatrix} \begin{Bmatrix} \boldsymbol{\mathcal{E}}\_{\varepsilon} \\ \boldsymbol{2\mathcal{E}}\_{s} \\ \boldsymbol{\mathcal{E}}\_{t} \end{Bmatrix} \right\rangle \tag{9}$$

where *<sup>D</sup>e*, *<sup>D</sup>es*, *<sup>D</sup>et*, *<sup>D</sup>s*, *<sup>D</sup>st* and *<sup>D</sup><sup>t</sup>* are the corresponding sub-matrices of a 3D 6 <sup>×</sup> 6 material matrix.

The virtual work done by the applied load may be stated as

$$
\delta\overline{\mathcal{W}}\_{3D} = \delta\overline{\mathcal{W}}\_{2D} + \delta\overline{\mathcal{W}}^\* \tag{10}
$$

where *δW*2*<sup>D</sup>* and *δW*<sup>∗</sup> denote, respectively, the virtual work independent and dependent of the fluctuation function, and

$$
\delta \overline{\mathcal{W}}\_{2D} = \langle p\_i \delta v\_i + q\_d \delta v\_{3,a} \rangle\_\prime \delta \overline{\mathcal{W}}^\ast = \left\langle \langle f\_i \delta w\_i \rangle + \mathfrak{r}\_i \delta w\_i^+ + \beta\_i \delta w\_i^- \right\rangle \tag{11}
$$

where (·)<sup>+</sup> and (·)<sup>−</sup> represent the items acting on the top and bottom of the plate, respectively; *τ<sup>i</sup>* and *β<sup>i</sup>* are the traction forces on the top and bottom surface of the plate, respectively; *fi* are the body forces; *pi* = *fi* + *ai* + *βi*, *qa* = *h*/2(*β<sup>a</sup>* − *aa*) − *x*<sup>3</sup> *fa*.

The variation of the total potential energy may be written as

$$\delta \Pi = \delta \mathcal{U} - \delta \mathcal{W}^\* = \frac{1}{2} \delta \left< \mathcal{E}^T \mathcal{D} \mathcal{E} \right> - \left< \left< f\_i \delta w\_i \right> + \pi\_i \delta w\_i^+ + \beta\_i \delta w\_i^- \right> \tag{12}$$

where only the unknown fluctuation function *wi* is changeable.

#### *2.1. Dimensional Reduction Analysis of 3D-TCP*

The zeroth-order fluctuation function may be solved by minimizing the zeroth-order approximation strain energy under the constraint of Equation (9) as

$$
\delta \mathcal{U}\_0 = 0 \tag{13}
$$

where

$$2lI\_0 = \left\langle \begin{array}{c} \left(\gamma + x\_3 \kappa\right)^{\mathrm{T}} \mathbf{D}\_t \left(\gamma + x\_3 \kappa\right) + w\_{||,\beta}^{\mathrm{T}} \mathbf{D}\_s w\_{||,\beta} + w\_{3,3}^{\mathrm{T}} \mathbf{D}\_t w\_{3,3} \\ + 2\left(\gamma + x\_3 \kappa\right)^{\mathrm{T}} \left(\mathbf{D}\_{ct} w\_{||,\beta} + \mathbf{D}\_{ct} w\_{3,3}\right) + 2w\_{||,\beta}^{\mathrm{T}} \mathbf{D}\_{st} w\_{3,3} \end{array} \right\rangle \tag{14}$$

The Lagrange multipliers *λ<sup>i</sup>* are used to impose the constraint on the fluctuation function as

$$
\delta(\Pi + \lambda\_i \langle w\_i \rangle) = 0 \tag{15}
$$

The zeroth-order approximate variational expression can be obtained as

$$
\begin{pmatrix}
\left[\left(\boldsymbol{\gamma} + \mathbf{x}\_{3}\mathbf{x}\right)^{\mathrm{T}}\mathbf{D}\_{\mathrm{cs}} + \boldsymbol{w}\_{||,\mathcal{B}}^{\mathrm{T}}\mathbf{D}\_{\mathrm{s}} + \boldsymbol{w}\_{3,\mathcal{B}}^{\mathrm{T}}\mathbf{D}\_{\mathrm{st}}^{\mathrm{T}}\right] \delta\boldsymbol{w}\_{||,\mathcal{B}} \\
+ \lambda\_{l} \delta\boldsymbol{w}\_{l} + \left[\left(\boldsymbol{\gamma} + \mathbf{x}\_{3}\mathbf{x}\right)^{\mathrm{T}}\mathbf{D}\_{\mathrm{ct}} + \boldsymbol{w}\_{||,\mathcal{B}}^{\mathrm{T}}\mathbf{D}\_{\mathrm{st}} + \boldsymbol{w}\_{3,\mathcal{B}}^{\mathrm{T}}\mathbf{D}\_{\mathrm{l}}\right] \delta\boldsymbol{w}\_{3,\mathcal{B}}
\end{pmatrix} = \boldsymbol{0} \tag{16}
$$

By partly integrating Equation (16), the relevant Euler-Lagrange equation may be obtained as

$$\begin{aligned} \left[ \left( \gamma + \mathbf{x}\_3 \kappa \right)^{\mathrm{T}} \mathbf{D}\_{\mathrm{cs}} + \left( \boldsymbol{w}\_{\parallel, 3} \right)^{\mathrm{T}} \mathbf{D}\_{\mathrm{s}} + \boldsymbol{w}\_{3, 3} \mathbf{D}\_{\mathrm{sf}} \right]\_{\mathrm{sf}} &= \lambda\_{\parallel} \\ \left[ \left( \gamma + \mathbf{x}\_3 \kappa \right)^{\mathrm{T}} \mathbf{D}\_{\mathrm{cf}} + \left( \boldsymbol{w}\_{\parallel, 3} \right)^{\mathrm{T}} \mathbf{D}\_{\mathrm{sf}} + \boldsymbol{w}\_{3, 3} \mathbf{D}\_{\mathrm{f}} \right]\_{\mathrm{,3}}^{\mathrm{3}} &= \lambda\_3 \end{aligned} \tag{17}$$

where *<sup>λ</sup>*|| = [*λ*<sup>1</sup> *<sup>λ</sup>*2] T.

According to the free surface condition, the expressions in square brackets of Equation (17) should be zero at the top and bottom of the plate, such as

$$\begin{aligned} \left[ \left( \gamma + \mathbf{x}\_3 \kappa \right)^{\mathrm{T}} \mathbf{D}\_{\mathrm{c}s} + w\_{||,3}^{\mathrm{T}} \mathbf{D}\_{\mathrm{s}} + w\_{3,3} \mathbf{D}\_{\mathrm{st}}^{\mathrm{T}} \right]^{+/-} &= \mathbf{0} \\ \left[ \left( \gamma + \mathbf{x}\_3 \kappa \right)^{\mathrm{T}} \mathbf{D}\_{\mathrm{c}t} + w\_{||,3}^{\mathrm{T}} \mathbf{D}\_{\mathrm{st}} + w\_{3,3} \mathbf{D}\_{\mathrm{t}} \right]^{+/-} &= \mathbf{0} \end{aligned} \tag{18}$$

where the superscript "+ / −" denotes the items at the top and bottom of the plate. *<sup>w</sup>*|| and *<sup>w</sup>*<sup>3</sup> can be solved by putting Equation (18) into Equation (17), such as

$$\boldsymbol{w}\_{\parallel} = \left\langle - (\boldsymbol{\gamma} + \boldsymbol{x}\_{3} \boldsymbol{\kappa}) \overline{\boldsymbol{D}}\_{\rm cs} (\boldsymbol{D}\_{\rm s})^{-1} \right\rangle^{\rm T}, \boldsymbol{w}\_{\rm 3} = \left\langle - (\boldsymbol{\gamma} + \boldsymbol{x}\_{3} \boldsymbol{\kappa}) \overline{\boldsymbol{D}}\_{\rm t} (\boldsymbol{D}\_{\rm t})^{-1} \right\rangle \tag{19}$$

where

$$\begin{cases} \overline{\mathbf{D}}\_{\rm cs} = \mathbf{D}\_{\rm cs} - \overline{\mathbf{D}}\_{\rm ct} (\mathbf{D}\_{\rm st})^{\rm T} \left(\overline{\mathbf{D}}\_{\rm t}\right)^{-1}, \quad \overline{\mathbf{D}}\_{\rm ct} = \mathbf{D}\_{\rm ct} - \mathbf{D}\_{\rm cs} (\mathbf{D}\_{\rm s})^{-1} \mathbf{D}\_{\rm st}, \\\ \overline{\mathbf{D}}\_{\rm t} = \mathbf{D}\_{\rm t} - \left(\mathbf{D}\_{\rm st}\right)^{\rm T} \left(\mathbf{D}\_{\rm s}\right)^{-1} \mathbf{D}\_{\rm st} \end{cases} \tag{20}$$

Substituting Equation (20) into Equation (14), we obtain

*U*2*<sup>D</sup>* = <sup>1</sup> 2 (*γ* + *x*3*κ*) <sup>T</sup>*De*(*γ* + *x*3*κ*) = <sup>1</sup> 2 *γ κ* <sup>T</sup> *A B B*T *D γ κ* = <sup>1</sup> 2 ⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ *γ*<sup>11</sup> *γ*<sup>22</sup> 2*γ*<sup>12</sup> *κ*<sup>11</sup> *κ*<sup>22</sup> 2*κ*<sup>12</sup> ⎫ ⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎭ T⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *A*<sup>11</sup> *A*<sup>12</sup> *A*<sup>16</sup> *B*<sup>11</sup> *B*<sup>12</sup> *B*<sup>16</sup> *A*<sup>12</sup> *A*<sup>22</sup> *A*<sup>26</sup> *B*<sup>12</sup> *B*<sup>22</sup> *B*<sup>26</sup> *A*<sup>16</sup> *A*<sup>26</sup> *A*<sup>66</sup> *B*<sup>16</sup> *B*<sup>26</sup> *B*<sup>66</sup> *B*<sup>11</sup> *B*<sup>12</sup> *B*<sup>16</sup> *D*<sup>11</sup> *D*<sup>12</sup> *D*<sup>16</sup> *B*<sup>12</sup> *B*<sup>22</sup> *B*<sup>26</sup> *D*<sup>12</sup> *D*<sup>22</sup> *D*<sup>26</sup> *B*<sup>16</sup> *B*<sup>26</sup> *B*<sup>66</sup> *D*<sup>16</sup> *D*<sup>26</sup> *D*<sup>66</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ *γ*<sup>11</sup> *γ*<sup>22</sup> 2*γ*<sup>12</sup> *κ*<sup>11</sup> *κ*<sup>22</sup> 2*κ*<sup>12</sup> ⎫ ⎪⎪⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎪⎪⎭ (21)

where *A*, *D* and *B* are tensile, bending, and coupling stiffness sub-matrices, respectively, and can be expressed as

$$\begin{array}{l} \mathbf{A} = \langle \langle \overline{\mathcal{D}}\_{\varepsilon} \rangle \rangle\_{\prime}, \mathbf{B} = \langle \langle \mathbf{x}\_{3} \overline{\mathcal{D}}\_{\varepsilon} \rangle \rangle\_{\prime}, \mathbf{D} = \langle \langle \mathbf{x}\_{3}^{2} \overline{\mathcal{D}}\_{\varepsilon} \rangle \rangle\_{\prime} \\\ \overline{\mathcal{D}}\_{\varepsilon} = \mathbf{D}\_{\varepsilon} - \overline{\mathcal{D}}\_{\varepsilon \varepsilon} \mathbf{D}\_{s}^{-1} \mathbf{D}\_{\varepsilon \varepsilon}^{\mathrm{T}} - \overline{\mathcal{D}}\_{\varepsilon t} \mathbf{D}\_{ct}^{\mathrm{T}} / \overline{\mathcal{D}}\_{t} \end{array} \tag{22}$$

The stiffness matrix provides the essential information of 3D-TCP and may be easily utilized in the shell elements in a finite element software to perform macroscopic plate analysis. Because it only concerns the 2D field variables in terms of the macro-coordinates *x*<sup>1</sup> and *x*2, the macroscopic behavior of the plate is governed by the the strain energy in Equation (21). As a result, the 2D-EPM may be used to represent the original 3D-TCP in the global analysis, and can be solved using the linear analysis solver in a finite element soft package like ABAQUS/Standard.

#### *2.2. Local Field Analysis*

A well-established equivalent model may be applied not only to global analysis, but also to local field analysis. To improve the equivalent model, the local field recovery relations should be given.

Equation (3) may be used to recover the local 3D displacement field, such as

$$u\_i = \vec{u}\_i + \begin{bmatrix} \vec{u}\_{1,1} & \vec{u}\_{1,2} & \vec{u}\_{1,3} \\ \vec{u}\_{2,1} & \vec{u}\_{2,2} & \vec{u}\_{2,3} \\ \vec{u}\_{3,1} & \vec{u}\_{3,2} & \vec{u}\_{3,3} \end{bmatrix} \begin{Bmatrix} y\_1 \\ y\_2 \\ y\_3 \end{Bmatrix} + \zeta w\_i \tag{23}$$

The local strain field can be recovered as

$$\mathcal{E}\_c^0 = \gamma + \mathbf{x}\_3 \kappa, \quad \mathfrak{L}\_s^0 = -\mathfrak{w}\_{||.3^\prime} \quad \mathcal{E}\_t^0 = \mathfrak{w}\_{3,3} \tag{24}$$

The Hooke's law may be used to recover the local stress field as

$$
\sigma = \overline{D}\mathcal{E} \tag{25}
$$

#### **3. Validation Example**

In this part, numerical examples of bending, buckling, and free vibration of 3D-TCP under various conditions are utilized to validate the accuracy and efficiency of 2D-EPM. The comparative analysis is depicted in Figure 2. The relative error between 2D-EPM and 3D-FEM is calculated as

$$Error = \frac{|\text{ 2D - EPM results} - \text{3D -FEM results}|}{\text{3D -FEM results}} \times 100\%. \tag{26}$$

**Figure 2.** Comparative analysis of 2D-EPM and 3D-FEM.

The microstructure of 3D textile composite plate as shown in Figure 3 is very complex, including three layers with two weft yarns in each layer, two layers with two warp yarns in each layer and two binder yarns. The geometry of unit cell is determined by several parameters as shown in Figure 4: (1) the interval *t* between layers along the *Z* direction; (2) the interval length *l*<sup>1</sup> between the weft yarns along the *Y* direction; (3) the interval length *l*<sup>2</sup> between the warp yarns along the *X* direction, where *l*<sup>2</sup> = *nl*<sup>1</sup> (*n* < 1); (4) the warp or weft yarn width *b*<sup>1</sup> and the thickness *h*1; and (5) the binder yarn width *b*<sup>2</sup> = *nb*<sup>1</sup> and the thickness *h*<sup>2</sup> = 0.5*h*1. The *X*, *Y* and *Z* direction of the unit cell are the same as *y*1, *y*<sup>2</sup> and *y*<sup>3</sup> in the theoretical derivation, respectively.

**Figure 3.** The unit cell of 3D textile composite plate.

**Figure 4.** Structural parameters of unit cell within the 3D-TCP. (**a**) *Y* direction; (**b**) *X* direction.

In this section, these parameters are set as: *b*<sup>1</sup> = 0.8, *h*<sup>1</sup> = 0.2 mm, *l*<sup>1</sup> = 1.2 mm, *b*<sup>2</sup> = 0.4 mm, *h*<sup>2</sup> = 0.1 mm, *l*<sup>2</sup> = 1.0 mm, and *t* = 0 mm. The dimensions of unit cell are 3.6 mm × 2.4 mm × 1.32 mm. The constituents in the composite material are carbon

fiber (T-300) and epoxy resin 3601, and their properties are obtained from ref. [35], as shown in Table 1. The yarn has elliptical cross-section as shown in Figure 5b, and usually modeled as unidirectional composites with hexagonal pack as shown in Figure 5a. The variational asymptotic homogenization method is used to obtain the equivalent engineering constants of the yarn with 64% fiber volume fraction, as shown in Table 2. The geometry model of unit cell within the 3D-TCP is generated by open source TexGen software. The yarns and matrix, as well as the yarns themselves, are linked by common nodes, indicating that different parts are perfectly connected. The equivalent stiffness matrix of 3D-TCP obtained by the present model is provided in Table 3.

**Figure 5.** Unit cell model of the yarn.

**Table 1.** Constituent properties of matrix and fiber used in the yarn [35].


**Table 2.** Equivalent engineering constants of the yarn with a fiber volume fraction of 64%.


**Table 3.** Equivalent stiffness matrix of 3D-TCP (unit: SI).


The unit cell is repeated 12 times in the *X* direction and 8 times along the *Y* direction to construct the 3D-FEM of a 3D textile composite plate. The dimensions of this plate are 28.8 mm long, 28.8 mm wide, and 1.32 mm thick, respectively. After the mesh convergence study, a total of 10,609 shell elements (S4R) and 288,000 solid elements (C3D20) are used in 2D-EPM and 3D-FEM, respectively.

#### *3.1. Bending Analysis*

Three combinations of boundary and load conditions in Figure 6 are used in bending analysis, in which C represents clamped boundary, F for free boundary, and Path represents the comparative analysis path. Table 4 shows the bending behaviors predicted by 3D-FEM and 2D-EPM under various conditions.

**Figure 6.** Combinations of boundary and load conditions for bending analysis. (**a**) Case 1; (**b**) Case 2; and (**c**) Case 3.

**Table 4.** Comparison of displacement *U*<sup>3</sup> predicted by two models under various boundary and load conditions (unit: mm).

Table 4 shows that the displacement clouds of *U*<sup>3</sup> predicted by the two models are consistent, and the maximum error is 7.70% in Case 1, which may due to the fact that the shear strain is negligible in Case 2 and Case 3, but relatively large in Case 1. To more clearly illustrate the details of the displacement distributions (especially out-of-plane distributions), the displacement of *U*<sup>3</sup> along the analysis path are compared in Figure 7. It can be observed that the displacement error between 3D-FEM and 2D-EPM is relatively very small even for 3D-TCP with complex microstructures, and the displacement curve of 3D-FEM is smoother than that of 2D-EPM. The main reason is that there is only one node in *x*<sup>3</sup> direction of

2D-EPM under bending load, which can not smoothly simulate the continuous deformation along the thickness direction.

**Figure 7.** Comparison of displacements along analysis paths under various boundary and load conditions. (**a**) Case 1; (**b**) Case 2; and (**c**) Case 3.

#### *3.2. Local Field Recovery*

The internal structure of 3D-TCP is complex, and the warp yarns, weft yarns and binder yarns are intertwined with each other. The study of the local stress, strain and displacement distributions is of significance to the failure analysis of textile structures. In this section, the local fields within the unit cell at the center of the plate are recovered under the conditions in Case 2. Two paths (as shown in Figure 8) are selected to analyze the local stress, strain and displacement distribution.

**Figure 8.** Path selection in local field recovery analysis. (**a**) Path 1 along the *Z* direction; and (**b**) Path 2 along the *Y* direction.

The local displacement distributions within the recovered unit cell from 2D-EPM are similar with those in the selected unit cell from 3D-FEM, as shown in Table 5, with a maximum error of 1.18%. Figure 9a shows that the maximum displacement of *U* is located in the middle of Path 1, where the displacements of weft yarns are greater than those of matrix. Figure 9b shows that the curve of *U* along Path 2 is divided into two segments with a length of 1.2 mm (the interval length between the weft yarns). The larger value of *U* is located in the suspended area between the weft yarns, while the smaller value of *U* is located at 0.5 mm and 1.7 mm of Path 2, which belongs to the superimposed area of the warp yarns and weft yarns. The local stress distributions within the recovered unit cell from 2D-EPM are similar with those in the selected unit cell from 3D-FEM, as shown in Table 6, and the maximum error is is 5.17% in *σ*22, indicating that the recovered local stress fields from 2D-EPM are accurate.

**Table 5.** Comparison of local displacement field within the unit cell under the conditions in Case 2 (unit: mm).

**Figure 9.** Comparison of local displacement curves along Path 1 and Path 2 of the unit cell predicted by two models under the conditions in Case 2. (**a**) *U* along Path 1; (**b**) *U* along Path 2.

**Table 6.** Comparison of local stress field within the unit cell under the conditions in Case 2 (unit: MPa).

**Table 6.** *Cont.*

Figure 10a,c show that the curves of von Mises stress and *σ*<sup>11</sup> along Path 1 are divided into five segments, representing two layers of warp yarns and three layers of matrix, respectively. The local stresses are distributed unevenly within the unit cell of 3D-TCP. The local stress in the matrix is relatively small, while the local stress in the warp yarns fluctuates greatly, indicating that the yarns along the thickness direction are main bearing components. Figure 10b,d show that the curves of von Mises stress and *σ*<sup>22</sup> along Path 2 are also divided into five sections. The stress at the superimposed area of warp yarns and weft yarns (0.25–0.95 mm, 1.45–2.15 mm) is relatively small, while the stress in the suspended area of the weft yarns fluctuates greatly, which indicates that it is easy to be damaged under the loading.

**Figure 10.** Comparison of local stress distributions along Path 1 and Path 2 of the unit cell predicted by two models under the conditions in Case 2. (**a**) Von mises stress along Path 1; (**b**) Von mises stress along Path 2; (**c**) *σ*<sup>11</sup> along Path 1; and (**d**) *σ*<sup>22</sup> along Path 2.

Table 7 shows that the local strain distributions are consistent with those of local stress distributions, and the error of local strain between recovered unit cell and selected unit cell is less than 1%. Figure 11 shows the local strain distributions along Path 1 and Path 2 of the unit cell predicted by two models under the conditions in Case 2. It can be observed that the local strain distributions within the unit cell of 3D-TCP are non-homogeneous. The strain curves are divided into several segments due to the interpenetration of warp yarn, weft yarn and matrix.

**Figure 11.** Comparison of local strain curves along Path 1 and Path 2 of the unit cell predicted by two models under the conditions in Case 2. (**a**) *ε*<sup>11</sup> along Path 1; (**b**) *ε*<sup>11</sup> along Path 2; (**c**) *ε*<sup>13</sup> along Path 1; (**d**) *ε*<sup>13</sup> along Path 2; (**e**) *ε*<sup>33</sup> along Path 1; and (**f**) *ε*<sup>33</sup> along Path 2.


**Table 7.** Comparison of local strain field within the unit cell under the conditions in Case 2.
