Prediction of Bending, Buckling and Free-Vibration Behaviors of 3D Textile Composite Plates by Using VAM-Based Equivalent Model

To solve the microstructure-related complexity of a three-dimensional textile composite, a novel equivalent model was established based on the variational asymptotic method. The constitutive modeling of 3D unit cell within the plate was performed to obtain the equivalent stiffness, which can be inputted into the 2D equivalent model (2D-EPM) to perform the bending, free-vibration and buckling analysis. The correctness and effectiveness of the 2D-EPM was validated by comparing with the results from 3D FE model (3D-FEM) under various conditions. The influence of yarn width and spacing on the equivalent stiffness was also discussed. Finally, the effective performances of 3D textile composite plate and 2D plain-woven laminate with the same thickness and yarn content were compared. The results revealed that the bending, buckling and free-vibration behaviors predicted by 2D-EPM were in good agreement with 3D-FEM, and the local field distributions within the unit cell of 3D textile composite plate were well captured. Compared with the 2D plain-woven laminate, the displacement of 3D textile composite plate was relatively larger under the uniform load, which may due to the fact that the through-the-thickness constrains of the former are only dependent on the binder yarns, while the warp yarns and weft yarns of the latter are intertwined closely.


Introduction
In recent years, composite structures are more and more widely used as load-bearing structures [1][2][3][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 To characterize the quick change of in-plane material characteristics, the micro-coordinates y i = x i /ζ (ζ 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 u 1 (x α ; y i ) =ū 1 (x α ) − ζy 3ū3,1 (x α ) + ζw 1 (x α ; y i ) u 2 (x α ; y i ) =ū 2 (x α ) − ζy 3ū3,2 (x α ) + ζw 2 (x α ; y i ) u 3 (x α ; y i ) =ū 3 (x α ) + ζw 3 (x α ; y i ) (1) where the displacements of the 3D-FEM and 2D-EPM are represented by u i andū i , respectively; w i are the fluctuation functions to be solved. The underline terms should satisfy the following constraints: hū α (x α ) = u α + ζy 3 ū 3,α , hū 3 (x α ) = u 3 (2) where · denotes the volume integral of a unit cell. The fluctuation functions in Equation (3) are constrained as The strain field can be expressed according to 3D linear elasticity theory, such as The 3D strain field can be obtained by substituting Equation (1) into Equation (4) and ignoring higher-order terms according to VAM, ε 11 = γ 11 + ζy 3 κ 11 + w 1,1 2ε 12 = 2γ 12 + 2ζy 3 κ 12 + w 1,2 + w 2,1 where γ αβ and κ αβ are the in-plane strains and bending curvatures of 2D-EPM, respectively, and may be defined as We can define the three-dimensional strain field in matrix form to make derivation easier, such as where () || = [() 1 () 2 ] T , γ = γ 11 2γ 12 γ 22 T , κ = [κ 11 κ 12 + κ 21 κ 22 ] T , and The strain energy of the 3D-TCP may be expressed briefly as where D e , D es , D et , D s , D st and D t are the corresponding sub-matrices of a 3D 6 × 6 material matrix. The virtual work done by the applied load may be stated as where δW 2D and δW * denote, respectively, the virtual work independent and dependent of the fluctuation function, and where (·) + and (·) − represent the items acting on the top and bottom of the plate, respectively; τ i and β i are the traction forces on the top and bottom surface of the plate, respectively; f i are the body forces; The variation of the total potential energy may be written as where only the unknown fluctuation function w i is changeable.

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 where The Lagrange multipliers λ i are used to impose the constraint on the fluctuation function as The zeroth-order approximate variational expression can be obtained as By partly integrating Equation (16), the relevant Euler-Lagrange equation may be obtained as 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 where the superscript "+ / −" denotes the items at the top and bottom of the plate. w || and w 3 can be solved by putting Equation (18) into Equation (17), such as where Substituting Equation (20) into Equation (14), we obtain where A, D and B are tensile, bending, and coupling stiffness sub-matrices, respectively, and can be expressed as 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 1 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.

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 The local strain field can be recovered as The Hooke's law may be used to recover the local stress field as

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 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 1 between the weft yarns along the Y direction; (3) the interval length l 2 between the warp yarns along the X direction, where l 2 = nl 1 (n < 1); (4) the warp or weft yarn width b 1 and the thickness h 1 ; and (5) the binder yarn width b 2 = nb 1 and the thickness h 2 = 0.5h 1 . The X, Y and Z direction of the unit cell are the same as y 1 , y 2 and y 3 in the theoretical derivation, respectively.  In this section, these parameters are set as: b 1 = 0.8, h 1 = 0.2 mm, l 1 = 1.2 mm, b 2 = 0.4 mm, h 2 = 0.1 mm, l 2 = 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.  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.

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.   Table 4 shows that the displacement clouds of U 3 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 3 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 3 direction of 2D-EPM under bending load, which can not smoothly simulate the continuous deformation along the thickness direction.
2D-EPM under bending load, which can not smoothly simulate the continuous deformation along the thickness direction.

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. 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

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. 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.

Stress Selected Unit Cell Recovered Unit Cell Max. Error
Von Mises 0.70%    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.

Global Buckling Analysis
In this section, the global buckling of 3D-TCP under different conditions shown in Figure 12 Table 8 lists the first six buckling modes and loads of 3D-TCP predicted by the two models under the conditions in Case 6. The first six buckling modes predicted by 3D-FEM and 2D-EPM are consistent, and the maximum error of buckling critical load in each buckling mode is only 2%. The calculation time of 2D-EPM in buckling analysis is about 18 times faster than 3D-FEM, verifying the effectiveness of 2D-EPM in global buckling analysis of 3D-TCP.   Table 9 lists the first buckling modes and critical loads predicted by the two models under the conditions in Case 4, Case 5 and Case 7. The first buckling modes predicted by the two models are nearly identical, and the maximum error of the buckling load is only 2.69%, which verifies the accuracy of 2D-EPM in buckling analysis of 3D-TCP under different conditions. Table 9. Comparison of the first buckling modes and critical loads (N) of 3D-TCP in different cases predicted by two models.  Tables 10 and 11 show the first three vibration modes and natural frequencies of 3D-TCP under the boundary conditions in Cases 4 and 7. It is clear that the vibration modes predicted by 2D-EPM agree with those predicted by 3D-FEM. For example, the first, second buckling modes, respectively, have one and two half-waves along the x 1 axis, and the third buckling mode has two half-waves along the x 2 axis under the boundary condition in Case 4. The maximum natural frequency error is 8.67%, indicating 2D-EPM has high accuracy in free-vibration analysis of 3D-TCP.

Influence of Structural Parameters on Equivalent Stiffness
The structure of 3D-TCP is complex and has many parameters (see Figure 3). In Section 3.2, we can see that the local stress and strain in binder yarns are relatively greater, indicating the binder yarns plays a very important role in preventing interlayer separation in 3D-TCP. Therefore, it is very important to study the influence of binder yarn width on the equivalent stiffness, which can also provide guidance for the design of 3D-TCP. Secondly, the influence of warp (weft) yarn width on the equivalent stiffness are also investigated. Table 12 lists the structural parameters of 3D-TCP used in the parameter analysis.  Figures 13 and 14 show that the binder width and warp (weft) width have great influence on the equivalent tensile stiffness A 11 and bending stiffness D 11 . The values of A 11 and D 11 increase with the increasing warp (weft) width, but decrease with the increasing binder width. The main reason is that the increase of binder width will lead to the decrease of yarn content, further resulting in the decrease of A 11 and D 11 . However, the smaller the binder yarn width is, the smaller the constraint in the thickness direction will be, resulting in the easy delamination. Therefore, the binder yarn width should be adjusted to ensure enough stiffness and integrity in the design of 3D-TCP.

Comparison of Effective Performance between 2D-PWL and 3D-TCP with the Same Thickness
The 3D-TCP is developed from the 2D plain-woven laminate (2D-PWL), and has better mechanical properties. To compare the effective performance of the two textile composite plate, we establish the 2D-EPM of 2D-PWL and 3D-TCP with the same plate thickness and yarn content.
The structure parameters of unit cell within the 6-layered 2D-PWL as shown in Figure 15

Comparison of Bending Behaviors
The boundary conditions used in bending analysis are shown in Figure 16, and the predicted bending displacements are shown in Table 15. It can be observed that the bending displacement of 3D-TCP is greater than that of 2D-PWL under uniform load (Case 8 and Case 9). This may due to the fact that the warp yarns and weft yarns of 3D-TCP are not interlaced, and are only constrained by the binding yarns in the thickness direction. While the displacement of 3D-TCP is smaller than that of 2D-PWL under the concentrated force and bending moment (Case 10 and Case 11), indicating that the torsion resistance of 3D-TCP are better than 2D-PWL. The above results are consistent with the obtained equivalent stiffness in Tables 13 and 14. That is, the bending stiffness D 11 and D 22 of 2D-PWL are greater than those of 3D-TCP, while the torsional stiffness D 66 of 2D-PWL is smaller than that of 3D-TCP.   Table 16 shows the first four buckling modes predicted by 2D-PWL and 3D-TCP under the conditions in Case 7. It can be observed that the buckling modes predicted by 2D-PWL and 3D-TCP are almost the same. That is, the first and third buckling modes, respectively, have one and two half-waves along the Y direction, and the second buckling mode has two half-waves along the X direction, the fourth buckling mode has two asymmetric half-waves along the X and Y directions. It is worth noting that the first four critical loads of 3D-TCP are greater than 2D-PWL, indicating that 3D-TCP has better stability under lateral load.

Comparison of Free-Vibration Characteristics
The first vibration mode and natural frequency of 2D-PWL and 3D-TCP under different boundary conditions in Figure 12 are compared as shown in Table 17. The vibration modes of 3D-TCP are consistent with those of 2D-PWL, while the natural frequency of 3D-TCP is smaller than that of 2D-PWL, which is consistent with the obtained equivalent stiffness.

Conclusions
The VAM-based equivalent model (2D-EPM) of 3D textile composite plate is established for bending, buckling and free-vibration analysis. The following conclusions can be obtained: (1) The maximum errors of bending displacement and buckling load between 3D-FEM and 2D-EPM are within the range of engineering accuracy, and the displacement distributions along the analysis path predicted by two models have the same trend with small differences. Furthermore, the local stress, strain and displacement distributions within the recovered unit cell are well captured. In addition, the computational efficiency of 2D-EPM is greatly improved by reducing the number of nodes and elements.
(2) The width of binder/warp yarn mainly affect the stiffness components A 11 and D 11 . That is, A 11 and D 11 increase with the increasing width of warp or weft yarn, while decrease with increasing width of binder yarn, which may due to the fact that the increasing width of binder yarn will decrease the yarn content.
(3) Compared with 2D plain-woven laminate with the same thickness and yarn content, the 3D textile composite plate has smaller equivalent bending stiffness and larger torsional stiffness, resulting in large displacement and small natural frequency. This may be because the 3D textile composite plate is only constrained by the binding yarns in the thickness direction, while the warp yarns and weft yarns are intertwined closely in the 2D-PWL.  Data Availability Statement: Data available on request due to restrictions, e.g., privacy or ethical. The data presented in this study are available on request from the corresponding author. The data are not publicly available due to subsequent analyzes and publications.

Conflicts of Interest:
The authors declare no conflict of interest.