A Methodology to Obtain the Accurate RVEs by a Multiscale Numerical Simulation of the 3D Braiding Process

To accurately evaluate the mechanical performance of three-dimensional (3D) braiding composites, it is essential to consider the braiding process and generate realistic representative volume element (RVE) structures. An efficient simulation methodology based on truss elements was used to simulate the 3D four-directional (3D4D) braiding process utilizing the finite element method (FEM) on the macroscale. The goal was to obtain the spatial trajectories of yarns and establish the relationship between the braiding parameters and the preform structure. Based on the initial yarn topology, the yarns were discretized as bundles of virtual sub-yarns. Then, a temperature drop simulation using hybrid elements was implemented to deform the yarn cross-section and obtain the interior, surface, and corner cells on the mesoscale. The simulation results show good agreement with the experiment. A parametric study was deployed to identify the effect of the model input parameters on the computation cost and accuracy. Furthermore, the approach applies to the other braiding processes, such as the cylindrical braiding composite.


Introduction
Composite materials are an effective means of lightweight design [1]. Three-dimensional braiding composites are formed by interweaving yarns in space, and the mechanical performance in each direction can be designed by changing the braiding process [2]. In addition, a rational braiding process can produce near-net-shape products, which eliminates an additional machining process [3]. With these advantages, 3D braiding composites have significant application potential in the aerospace, marine, and transportation industries [4]. Nowadays, 3D braiding equipment is fully automated to produce composite components, which significantly increases the application and production efficiency of braiding composites.
The macro-mechanical properties of 3D braiding composites are primarily determined by the preform's mesostructure and the constituent materials' mechanical properties. Accurate mesostructure is a precondition for calculating the mechanical properties of composite products. The representative volume element (RVE) is the smallest volume over which a measurement can be made to yield a value representative of the whole [5]. Wang et al. [6] established three types of RVEs corresponding to the interior, surface, and corner cells, by observing the interweaving topology of yarns during the braiding process. Chen et al. [7] used an analytical approach combined with experimental observations to study the mathematical relationships between the braiding structure and the braiding process parameters. Zhang [8] analyzed the spatial motion trajectories of yarns for 3D five-directional and 3D full five-directional braiding composite materials. It was assumed that the shapes of the yarn section differed according to the compaction conditions of the yarns. Therefore, three in Section 4. In Section 5, the cylindrical braiding process is simulated to present the method's versatility. Concluding remarks are provided in Section 6.

Three-Dimensional Four-Directional Braiding Process
The 3D braiding fabric is woven by the alternate movements of yarn carriers in row and column directions. Figure 1 is a schematic illustration of the four-step 1 × 1 rectangular braiding process. The arrangement of the carrier pattern (Figure 1a) determines the section shape of the preform. A machine cycle is divided into four movements: in Step 1, the yarn carriers move alternately in one position laterally ( Figure 1c); in Step 2, carrier motions occur in alternate columns ( Figure 1d); Steps 3 and 4 reverse the carrier movement direction in Steps 1 and 2 (Figure 1e,f), as illustrated by the arrows in each step. After these four steps, the carrier pattern is identical to the original pattern ( Figure 1a). One segment of the preform is fabricated and called a braiding pitch, denoted by h. The take-up device moves up one pitch length simultaneously. The carriers regulate the tension force by retracting the yarns into the spool or releasing the yarns from the spool. Therefore, the yarns compact with each other, resulting in a stable and tight preform. The braiding process is fully automated, without a manual jamming process. Figure 1b shows an actual braiding setup, including cylinders in each line, the carriers, and the take-up module [27].   [27]. Reprinted with permission from [27]. Copyright 2020, The Regenerative Engineering Society. (c) Step 1,(d) Step 2, (e) Step 3, and (f) Step 4.  [27]. Reprinted with permission from [27]. Copyright 2020, The Regenerative Engineering Society. (c) Step 1,(d) Step 2, (e) Step 3, and (f) Step 4. The count of carriers in the rectangle pattern is determined by the row and column number of yarns of the main part [m × n], which is given as, where m and n are the row and column yarn numbers in the main part, respectively, and N is the total number of yarns. After 4N/t steps, all carriers return to the original position, where t denotes the greatest common divisor of m and n. These yarns in the [6 × 6] arrangement are divided into six groups according to their trajectories, as shown in Figure 2a. After moving 32 steps, all carriers return to their original positions. A Python script was written to obtain the path of every carrier and generate the boundary condition of the numerical simulation model. (e) (f) Figure 1. Schematic of the 3D four-directional braiding process: (a) the original arrangement, (b) the braiding machine [27]. Reprinted with permission from [27].   Table 1 provides the main parameters of the 3D4D braiding process, as shown in Figure 1b, which significantly affects the braiding preform structure. The preform  Table 1 provides the main parameters of the 3D4D braiding process, as shown in Figure 1b, which significantly affects the braiding preform structure. The preform structure is mainly characterized by the interior braiding angle γ, the surface braiding angle θ, and the pitch length h. In the specimen, the surface braiding angle and the pitch length were easy to observe and measure and thus were used to benchmark with the simulation model, as shown in Figure 2b. The high-strength 12K T700 carbon yarn was used to braid the preform. The primary material properties of the yarn are listed in Table 2.

Modeling Approach Using the FEA Method
The workflow of the proposed modeling methodology is illustrated in Figure 3, and the procedures are described as follows: Simulation of the 3D4D braiding process on the macroscale: The braiding process was simulated according to Section 2, during which the yarns were discretized as T3D2 truss elements. The parametric model was generated using a Python script and solved by the Abaqus software ( Figure 3a). Yarns were interweaved to achieve the preform. II. Sub-yarn discretization: The yarns in the two pitches of the preform at the stable braiding stage ( Figure 3b) were divided into several bundles of sub-yarns according to the center curve of the yarns achieved in Step I, as shown in Figure 3c. A Python script was written to automate the model generation and the assignment of boundary conditions. III.
Sub-yarn deforming simulation on the mesoscale: The contraction process was accomplished by imposing the temperature drop load with an appropriate thermal expansion coefficient (Figure 3d). The sub-yarns were represented by the hybrid elements consisting of the truss and beam element with coinciding nodes. Then, the yarn section changed from a circle to a realistic shape as the average stress of all sub-yarns reached the tension stress of the carrier retraction. IV.
Mesostructure reconstruction: Since the FEA analyses require yarns with a solid geometry rather than fibers, the chains of sub-yarns must be converted to solid geometry. Another Python script was written to extract the section shapes of different yarns and reconstruct the accurate yarn structure (Figure 3e). Eventually, three RVEs were generated: interior cell, surface cell, and corner cell (Figure 3f-g).

Mesoscale
Step IV: mesostructure Reconstruction Step III : Sub-yarn Deformation

Yarns
Carrier control point

Two pitches
Step I: simulation of the 3D4D braiding process on a macroscale

Simulation of the 3D4D Braiding Process on the Macroscale
The braiding process simulation model takes the machine control and process parameters as inputs and generates the braiding structure as an output. As shown in Figure 3a and Supplementary Video S1 (see Supplementary Materials), the simulation model on the macroscale contains the yarns and nonlinear springs. One end of the yarns

Simulation of the 3D4D Braiding Process on the Macroscale
The braiding process simulation model takes the machine control and process parameters as inputs and generates the braiding structure as an output. As shown in Figure 3a and Supplementary Video S1 (see Supplementary Materials), the simulation model on the macroscale contains the yarns and nonlinear springs. One end of the yarns is coupled to the take-up point, and the other end is connected to a spring element with constant tension force. Figure 2a displays the movement paths of the carriers applied to these control points. The take-up point moves upward at a constant speed corresponding to the carriers' speed. The yarn-to-yarn friction coefficients for parallel-and perpendicular-to-yarn orientations were experimentally measured to be 0.26 and 0.46 [23]. However, the inter-yarn friction coefficient cannot be changed according to the yarn's cross angle in the braiding simulation. Therefore, the interactions between yarns are defined using the Coulomb friction model with a friction coefficient of 0.36, which is the average of two directions. The braiding process is a quasi-static process, and a stable fabric structure requires multiple braiding cycles. It is essential to speed up the model calculation without compromising the computational accuracy. Reducing the modulus decreases the wave velocity in the material, facilitating a faster analysis speed. Comparative calculations show that setting the yarn's modulus as 20,000 MPa is adequate to avoid stretching the yarn and the contact penetration.
Since a yarn consists of a substantial number of fibers, the bending modulus of yarns is much lower than the value calculated using Young's modulus in the fiber direction. The primary purpose of the braiding process simulation is to achieve the topology structure of the preform. Using truss elements (T3D2) that neglect the bending stiffness approximates the actual state of yarns in the braiding process [22,28]. In Section 3.3, the yarn section deformation due to squeezing will be computed. The comparative analysis results show that the element size setting at half of the yarn radius can achieve minor penetration. Meanwhile, the contact control parameters were set to avoid a reduction in the contact thickness owing to the small element-size-to-truss-radius ratio [29]. The primary process parameters of the 3D4D braiding process in Table 1 and the trajectory of each yarn carrier were parameterized to facilitate modeling. Then, the corresponding Abaqus input file was automatically generated through the user-written Python script. Figure 4a shows that the simulation result of the braiding preform is stable and the length of each pitch is also equal. The yarns in the same group have the same path, as shown in Figure 4b,c. Figure 4e shows the test specimen. In addition, another Python script was written to achieve the pitch length and braiding angle along the braiding direction. The center curve of every yarn was drawn through Abaqus's application programming interface method. These curves were cut by sequential planes perpendicular to the braiding direction with a small resolution distance. The matching degree of two pair intersections was computed using the CDIST function in SciPy [30]. As shown in Figure 4d, the pitch length corresponding to one machine cycle was found when the error in the matching degree between two sections reached the minimum value. The interior braiding angle γ calculated from the yarn path in group 1 (Figure 5b) is the angle between the oblique line and the Z axis (the braiding direction) in the 45 • plane. The surface braiding angle θ calculated from the yarn path in group 5 is the angle between the oblique line and the Z axis in the XOZ plane.
It is well known that the braiding process parameters must match each other to obtain a consistent fabric structure. For example, the fixed height of the interweaving point and the distance between two yarn carriers determine the braiding angle. The take-up speed and the time for each carrier to move one step decided the interweaving point of the braiding process. Therefore, the model with the take-up velocity as zero was implemented to investigate the relationship between braiding angle and interweaving point height. Despite the fluctuations in the braiding angle and height measurements of the experimental specimens, the simulation results are in good agreement with the experimental results, as shown in Figure 5. It is evident that the relationship between braiding angle and interweaving height is nonlinear, and the angle decreases as the height increases ( Figure 5a). In addition, pitch length is positively correlated with braiding interweaving height (Figure 5b). The take-up machine moves up one pitch length in one braiding cycle to obtain a consistent preform. Meanwhile, it is necessary to match the pitch length and the interweaving point height by controlling the speed of the take-up device.
script was written to achieve the pitch length and braiding angle along the braiding direction. The center curve of every yarn was drawn through Abaqus's application programming interface method. These curves were cut by sequential planes perpendicular to the braiding direction with a small resolution distance. The matching degree of two pair intersections was computed using the CDIST function in SciPy [30]. As shown in Figure 4d, the pitch length corresponding to one machine cycle was found when the error in the matching degree between two sections reached the minimum value. The interior braiding angle calculated from the yarn path in group 1 (Figure 5b) is the angle between the oblique line and the Z axis (the braiding direction) in the 45° plane. The surface braiding angle calculated from the yarn path in group 5 is the angle between the oblique line and the Z axis in the XOZ plane. It is well known that the braiding process parameters must match each other to obtain a consistent fabric structure. For example, the fixed height of the interweaving point and the distance between two yarn carriers determine the braiding angle. The takeup speed and the time for each carrier to move one step decided the interweaving point of the braiding process. Therefore, the model with the take-up velocity as zero was implemented to investigate the relationship between braiding angle and interweaving point height. Despite the fluctuations in the braiding angle and height measurements of the experimental specimens, the simulation results are in good agreement with the experimental results, as shown in Figure 5. It is evident that the relationship between braiding angle and interweaving height is nonlinear, and the angle decreases as the height increases ( Figure 5a). In addition, pitch length is positively correlated with braiding interweaving height ( Figure 5b). The take-up machine moves up one pitch length in one braiding cycle to obtain a consistent preform. Meanwhile, it is necessary to match the pitch length and the interweaving point height by controlling the speed of the take-up device.

Sub-Yarn Discretization
The yarn center curves are obtained from the simulation method described in Section 3.1. However, the yarn section is circular due to the non-deformability of the truss element section, which is different from the actual yarn. To simulate the yarn squeezing deformation, the single yarn was discretized as several virtual sub-yarns whose fiber volume fraction is equal to the actual yarn, which was 65.3% in T700 12K [31]. The diameter of the sub-yarns denoted as is determined by the equation as follows: where and are the real yarn area and the fiber volume fraction, respectively, and and are the number and diameter of the virtual sub-yarns, respectively. The sub-yarn discretization is implemented in a custom Python script. Figure 6 shows the schematic of the sub-yarn discretization process.
… is the centerline of the yarn, the coordinates of which are extracted from the simulation result of the 3D braiding process in Section 3.1. , , … are the sections vertical to the centerline of the yarns.
is a point at the circumference of section , and vector ⃗ is achieved by rotating vector ⃗ counterclockwise 90° with ⃗ as the axis. Then, a line parallel to vector ⃗ is drawn from point , which intersects with at a point denoted as . intersects with section at . If the yarn is twisted, the corresponding rotating angle can be calculated by the formula below.

Sub-Yarn Discretization
The yarn center curves are obtained from the simulation method described in Section 3.1. However, the yarn section is circular due to the non-deformability of the truss element section, which is different from the actual yarn. To simulate the yarn squeezing deformation, the single yarn was discretized as several virtual sub-yarns whose fiber volume fraction is equal to the actual yarn, which was 65.3% in T700 12K [31]. The diameter of the sub-yarns denoted as D v f is determined by the equation as follows: where A real and V f v f are the real yarn area and the fiber volume fraction, respectively, and n v f and D v f are the number and diameter of the virtual sub-yarns, respectively. The sub-yarn discretization is implemented in a custom Python script. Figure 6 shows the schematic of the sub-yarn discretization process. C 1 C 2 . . . C n is the centerline of the yarn, the coordinates of which are extracted from the simulation result of the 3D braiding process in Section 3.1. S 0 , S 1 , . . . S n are the sections vertical to the centerline of the yarns. U 1 is a point at the circumference of section S 1 , and vector −−→ C 1 V 1 is achieved by rotating vector Then, a line parallel to vector −−→ C 1 C 2 is drawn from point U 1 , which intersects with S 2 at a point denoted as N 2 . C 2 N 2 intersects with section S 2 at Q 2 . If the yarn is twisted, the corresponding rotating angle θ i can be calculated by the formula below.
where φ is the twist of the yarn. The points C i , U i , and V i can be determined in sequence by repeating the earlier process. Then, the coordinate transform relationship M i from C 0 U 0 V 0 to C i U i V i is as follows: Polymers 2022, 14, 4210 9 of 20 where is the twist of the yarn. The points , , and can be determined in sequence by repeating the earlier process. Then, the coordinate transform relationship from to is as follows: However, in some exceptional cases, the coordinate matrix does not have an inverse matrix, making it impossible to calculate the transformation matrix . A new method to handle all situations was proposed, including one translational transformation and two rotational transformations. The method includes the following steps: The translational transformation was conducted to ensure that points and coincide, as shown in Figure 7a. The coordinates of points and are denoted as , , and , , , respectively. The vector ⃗ is denoted as ( , , ) . The translation matrix is as follows: II. The rotational transformation was deployed to ensure that points and coincide. The rotation axis ⃗ is determined by the cross product of vectors ⃗ and ⃗ , as shown in Figure 7b. The unit vector ⃗ of the rotation direction and the rotation angle are calculated as follows: The rotation matrix is presented in Equation (8), where = ( ) and = ( ). The blue lines (Figure 7b) are sub-yarns after rotation. However, in some exceptional cases, the coordinate matrix does not have an inverse matrix, making it impossible to calculate the transformation matrix M i . A new method to handle all situations was proposed, including one translational transformation and two rotational transformations. The method includes the following steps: The translational transformation was conducted to ensure that points C 0 and C i coincide, as shown in Figure 7a. The coordinates of points C 0 and C i are denoted The rotational transformation was deployed to ensure that points U i and U i coincide. The rotation axis Figure 7b. The unit vector → n of the rotation direction and the rotation angle ϕ i1 are calculated as follows: The rotation matrix M i2 is presented in Equation (8), where c = cos(ϕ i1 ) and s = sin(ϕ i1 ). The blue lines (Figure 7b) are sub-yarns after rotation. III.
Taking the vector −−→ C i U i as the axis, the angle ϕ i2 between −−→ C i V i and −−→ C i V i is rotated so that point V i coincides with point V i and the corresponding rotation matrix M i3 is obtained, as illustrated in Figure 7c. After the three transformations, points C 0 , U 0 , and V 0 coincide at points C i , U i , and V i , respectively. All nodes in the original sub-yarn layout can be transferred to the target location by the following expression: where refers to the section number and refers to the number of sub-yarns. The subscripts and represent the section number and the sub-yarn number, respectively. The coordinate of the jth point on the ith section is ( , , , , , ). Figure 7d shows the discretization results of the yarns with 61 chains of sub-yarns. There is no interpenetration among the yarn sections. It is also able to achieve twisted yarn through this method.

Sub-Yarn Deformation on the Mesoscale
The extrusion deformation among yarns is the dominant form in sub-yarn shrinkage All nodes in the original sub-yarn layout can be transferred to the target location by the following expression: where m refers to the section number and n refers to the number of sub-yarns. The subscripts i and j represent the section number and the sub-yarn number, respectively. The coordinate of the jth point on the ith section is x p i,j , y p i,j , z p i,j . Figure 7d shows the discretization results of the yarns with 61 chains of sub-yarns. There is no interpenetration among the yarn sections. It is also able to achieve twisted yarn through this method.

Sub-Yarn Deformation on the Mesoscale
The extrusion deformation among yarns is the dominant form in sub-yarn shrinkage simulations. The section shape of the yarn is mainly influenced by fiber bending stiffness that affects the realignment resistances, so the bending stiffness of virtual sub-yarns cannot be neglected. A new mesh technique was used to consider the appropriate bending stiffness of yarns by overlaying the truss element (T3D2) with beam elements (B31) with the coincident nodes [32], as shown in Figure 8a. The hybrid element, where the truss element determines the tensile stiffness (fiber orientation), while the beam element represents the bending property of the sub-yarns with a lower Young's modulus, can simulate the real mechanical property. The beam radius is the same as the truss's radius. The bending stiffness of the beam elements is calculated as follows: where I = 0.25πr 4 is the second moment of inertia and r is the radius of fiber; n v f and n r f denote the number of sub-yarns and actual fiber, respectively; E beam and E r f denote the Young's modulus of the beam element and the actual fiber, respectively.
Polymers 2022, 14, 4210 11 of 20 simulate the real mechanical property. The beam radius is the same as the truss's radius. The bending stiffness of the beam elements is calculated as follows: where = 0.25 is the second moment of inertia and is the radius of fiber; and denote the number of sub-yarns and actual fiber, respectively; and denote the Young's modulus of the beam element and the actual fiber, respectively.
The 3D4D braiding composite is periodic in the braiding direction. Taking into account the edge effects and the computing cost, two pitches in a stable braiding state were chosen for sub-yarn deformation compression analysis. The end nodes of every bundle of sub-yarns were coupled to a control point using a uniform coupling constraint. The periodic boundary conditions (PBC) were applied to the two corresponding control points (Figure 8b). The linear constraint equations (*Equation) in Abaqus were adopted to realize the displacement constraint of the control points, and the formula was as follows: where refers to the displacement of the freedom , subscripts and refer to the control point numbers and refers to the quantity of yarn. To achieve a realistic braiding fabric structure, the temperature dropping load was imposed on all elements to ensure that the yarns squeezed with each other, which was equivalent to exerting tension force during the braiding process. The temperature drop stopped until the average stress of all elements was equal to the tension stress produced by the tension force of the yarn carrier. The following equation determines the average stress: where is the average stress; , and , are the tension stress of the beam and truss elements, respectively; and denotes the total element number of beam elements. The virtual yarn structure and the experiment results were compared, as shown in  The 3D4D braiding composite is periodic in the braiding direction. Taking into account the edge effects and the computing cost, two pitches in a stable braiding state were chosen for sub-yarn deformation compression analysis. The end nodes of every bundle of sub-yarns were coupled to a control point using a uniform coupling constraint. The periodic boundary conditions (PBC) were applied to the two corresponding control points (Figure 8b). The linear constraint equations (*Equation) in Abaqus were adopted to realize the displacement constraint of the control points, and the formula was as follows: where u i refers to the displacement of the freedom i, subscripts a and b refer to the control point numbers and N refers to the quantity of yarn. To achieve a realistic braiding fabric structure, the temperature dropping load was imposed on all elements to ensure that the yarns squeezed with each other, which was equivalent to exerting tension force during the braiding process. The temperature drop stopped until the average stress of all elements was equal to the tension stress produced by the tension force of the yarn carrier. The following equation determines the average stress: where S avg is the average stress; s beam 11,i and S truss 11,i are the tension stress of the beam and truss elements, respectively; and N denotes the total element number of beam elements.
The virtual yarn structure and the experiment results were compared, as shown in Figure 9. The localized distortions are observed, and the fibers are no longer parallel to the yarn centerline due to the lateral compression of adjacent yarns. The axes of yarns are spatial curves, and the section deformations are complicated. The shape of the interior yarn in simulation (Figure 9a) [13]. Reprinted with permission from [13], Copyrigh 2017, Springer Science Business Media B.V., (c) the simulation surface yarn, (d) the micro-CT surfac yarn [13]. Reprinted with permission from [13], Copyright 2017, Springer Science Business Medi B.V, (e) the surface in the simulation, (f) the surface in the experiment [33]. Reprinted wit permission from [33], Copyright 2012 Elsevier Ltd., (g) the 45° cross-section in simulation, and (h the 45° cross-section in the experiment [34]. Reprinted with permission from [34], Copyright 201 Elsevier Ltd.

Mesostructure Reconstruction
Most braiding composite FEA analyses are based on the yarn level. Fiber-leve geometry needs to be transferred to yarn-level geometry. The center fiber (Figure 10a determines the sequential cross-sections' geometric centers and normals. The detaile Figure 9. Comparison of the virtual structure and the experimental results: (a) the simulation interior yarn, (b) the micro-CT interior yarn [13]. Reprinted with permission from [13], Copyright 2017, Springer Science Business Media B.V., (c) the simulation surface yarn, (d) the micro-CT surface yarn [13]. Reprinted with permission from [13], Copyright 2017, Springer Science Business Media B.V, (e) the surface in the simulation, (f) the surface in the experiment [33]. Reprinted with permission from [33], Copyright 2012 Elsevier Ltd., (g) the 45 • cross-section in simulation, and (h) the 45 • cross-section in the experiment [34]. Reprinted with permission from [34], Copyright 2013 Elsevier Ltd.

Mesostructure Reconstruction
Most braiding composite FEA analyses are based on the yarn level. Fiber-level geometry needs to be transferred to yarn-level geometry. The center fiber (Figure 10a) determines the sequential cross-sections' geometric centers and normals. The detailed extraction process of the yarn geometry involves three steps: The evenly distributed rays are drawn along the circumferential direction, with the center of the section as the starting point, as shown in Figure 10c. The center of the big circle moves outward from the center in steps of 0.001 mm on the ray until it does not intersect all sub-yarn sections. The big circle's radius is set as 3-5 times the fiber radius to achieve a smooth and precise profile. II.
The centers of all big circles are connected to form a curve. Then, the curve is offset by the radius of the big circle inward to obtain a curve that surrounds all fibers, that is, the cross-sectional profile of the yarn. III.
All the cross-sectional profiles are joined to a multi-section surface, as shown in Figure 10b. These procedures are implemented with user-written Python scripts to obtain the coordinate sub-yarn nodes from the FEA result and generate realistic yarn geometry in Catia software. center of the section as the starting point, as shown in Figure 10c. The center of the big circle moves outward from the center in steps of 0.001 mm on the ray until it does not intersect all sub-yarn sections. The big circle's radius is set as 3-5 times the fiber radius to achieve a smooth and precise profile. II.
The centers of all big circles are connected to form a curve. Then, the curve is offset by the radius of the big circle inward to obtain a curve that surrounds all fibers, that is, the cross-sectional profile of the yarn. III.
All the cross-sectional profiles are joined to a multi-section surface, as shown in Figure 10b. These procedures are implemented with user-written Python scripts to obtain the coordinate sub-yarn nodes from the FEA result and generate realistic yarn geometry in Catia software.  Figure 11a shows the extracted structure in one pitch. All yarns are anti-symmetric with respect to the center of the preform. In the 6 × 6 braiding preform, there are four interior cells, eight surface cells, and four corner cells, as presented in Figure 11b. Furthermore, the most detailed description of the mesostructure can be procured by comparing the changes in the sequential cross-sections of the preform along the braiding direction [9]. As shown in Figure 11b-g, four sections in one pitch are cut along the braiding direction, with equal distance between each section, as ℎ 5 ⁄ . It is obvious that the yarn layouts of Figure 11b,f are almost identical, indicating that the mesostructure of 3D4D braiding composites is repeated for every pitch.  Figure 11a shows the extracted structure in one pitch. All yarns are anti-symmetric with respect to the center of the preform. In the [6 × 6] braiding preform, there are four interior cells, eight surface cells, and four corner cells, as presented in Figure 11b. Furthermore, the most detailed description of the mesostructure can be procured by comparing the changes in the sequential cross-sections of the preform along the braiding direction [9]. As shown in Figure 11b-g, four sections in one pitch are cut along the braiding direction, with equal distance between each section, as h/5. It is obvious that the yarn layouts of Figure 11b,f are almost identical, indicating that the mesostructure of 3D4D braiding composites is repeated for every pitch.
The essential purpose of our work is to obtain accurate RVEs for subsequent mechanical calculations. The solid structure models of three types of RVEs, including the interior, surface, and corner cells, are periodic, as shown in Figure 4f-h. Thus, it is conducive to generate the periodic mesh and impose the periodic conditions in subsequent mechanical calculations. The essential purpose of our work is to obtain accurate RVEs for subsequent mechanical calculations. The solid structure models of three types of RVEs, including the interior, surface, and corner cells, are periodic, as shown in Figure 4f-h. Thus, it is conducive to generate the periodic mesh and impose the periodic conditions in subsequent mechanical calculations.

Parametric Study
The parameter study is performed to verify the accuracy and convergence performance of the finite element model. How to set the number of virtual sub-yarns to balance the calculation accuracy and efficiency, and the effects of sub-yarn stiffness, tension force, and yarn twist on the structure need further study.

Geometrical Convergence
Owing to the enormous computational effort required, it is impossible to model each fiber as a single chain of beams. A parametric study of the sub-yarn number must be performed to balance the discrete precision of yarn cross-section deformation and the computation cost. Four models, discretized as the 19, 61, 92, and 133 chains of sub-yarns per yarn with the same yarn bending stiffness, were simulated to study the convergence of the model to reproduce an accurate braid structure. Figure 12a shows the result after the temperature drop process, and the cross-section in the middle was chosen to observe the deformation. Figure 12b illustrates the average stress representing the tension force of yarns versus the temperature drop value during the deformation. It is noted that each curve could be divided into three parts: linear, non-linear, and linear. The curves in the first part are flat, with low tension stress, and a kinematic response is evident where the virtual fibers rearrange and fill the gaps. The second parts show non-linear behavior, where the fibers reorient with increasing contact. In the third part, the average stress increases rapidly when fibers are sufficiently rearranged, resulting in a higher fiber volume fraction.
A higher number of sub-yarns leads to greater freedom of deformation of the yarn, so the average stress rises more slowly. The curve of 92 sub-yarns almost coincides with the curve of 133 sub-yarns, demonstrating that the model converges using 92 sub-yarns.

Parametric Study
The parameter study is performed to verify the accuracy and convergence performance of the finite element model. How to set the number of virtual sub-yarns to balance the calculation accuracy and efficiency, and the effects of sub-yarn stiffness, tension force, and yarn twist on the structure need further study.

Geometrical Convergence
Owing to the enormous computational effort required, it is impossible to model each fiber as a single chain of beams. A parametric study of the sub-yarn number must be performed to balance the discrete precision of yarn cross-section deformation and the computation cost. Four models, discretized as the 19, 61, 92, and 133 chains of sub-yarns per yarn with the same yarn bending stiffness, were simulated to study the convergence of the model to reproduce an accurate braid structure. Figure 12a shows the result after the temperature drop process, and the cross-section in the middle was chosen to observe the deformation. Figure 12b illustrates the average stress representing the tension force of yarns versus the temperature drop value during the deformation. It is noted that each curve could be divided into three parts: linear, non-linear, and linear. The curves in the first part are flat, with low tension stress, and a kinematic response is evident where the virtual fibers rearrange and fill the gaps. The second parts show non-linear behavior, where the fibers reorient with increasing contact. In the third part, the average stress increases rapidly when fibers are sufficiently rearranged, resulting in a higher fiber volume fraction.
A higher number of sub-yarns leads to greater freedom of deformation of the yarn, so the average stress rises more slowly. The curve of 92 sub-yarns almost coincides with the curve of 133 sub-yarns, demonstrating that the model converges using 92 sub-yarns. Meanwhile, the ratio of total width change (σ tw = tw/tw 0 ), representing the overall degree of deformation when the average stresses are equal to 1 MPa, is plotted in Figure 12c, showing that the number of sub-yarns has almost no effect on the overall deformation. Furthermore, three types of yarns, surface, corner, and interior, are compared in Figure 12d. The surface and corner yarns change from ellipse to pentagon, while the interior yarns change from ellipse to hexagonal. The cross-sectional shape of 92 sub-yarns is the same as that of 133 sub-yarns, and there is no penetration between the sub-yarns. Hereafter, 92 sub-yarns were used to model the yarn. degree of deformation when the average stresses are equal to 1 MPa, is plotted in Figure  12c, showing that the number of sub-yarns has almost no effect on the overall deformation. Furthermore, three types of yarns, surface, corner, and interior, are compared in Figure 12d. The surface and corner yarns change from ellipse to pentagon, while the interior yarns change from ellipse to hexagonal. The cross-sectional shape of 92 sub-yarns is the same as that of 133 sub-yarns, and there is no penetration between the sub-yarns. Hereafter, 92 sub-yarns were used to model the yarn.

The Bending Stiffness
The changes in fiber bending stiffness by imposing beam elements of different bending stiffness affect yarn spreading. Figure 13a

The Bending Stiffness
The changes in fiber bending stiffness by imposing beam elements of different bending stiffness affect yarn spreading. Figure 13a illustrates the average stress versus temperature drop values, which are plotted for three different moduli of the beam elements, of 133.6 MPa, 1336 MPa (the baseline, (EI) yarn = 0.3), and 13,360 MPa, with the same friction coefficient and tensile stiffness of the truss elements. The results show that the average stress rises faster as the bending modulus increases, indicating that the bending modulus cannot be neglected. Moreover, the degree of cross-sectional deformation is characterized by the ratio (σ = l/w) of the length (l) to the width (w) at a specific tension force (see Figure 12a). A smaller bending modulus of yarns, whether surface yarn, corner yarn, or interior yarn, leads to lower resistance of the yarn spreading to transverse deformation (Figure 13b), resulting in a larger ratio σ. The total width variation (σ tw ) increases as the bending stiffness decreases at the same yarn tension force, resulting in a higher yarn volume fraction. specific tension force (see Figure 12a). A smaller bending modulus of yarns, whether surface yarn, corner yarn, or interior yarn, leads to lower resistance of the yarn spreading to transverse deformation (Figure 13b), resulting in a larger ratio . The total width variation ( ) increases as the bending stiffness decreases at the same yarn tension force, resulting in a higher yarn volume fraction.

The Tension Force of Carriers and the Twisting of Yarns
The 3D braiding process is fully automated and does not require manual jamming. The yarn carriers can control the yarn tension force at a constant value by retracting or releasing the yarns. Three simulations were performed with an increasing tension force of 0.5 N, 1 N, and 3 N, respectively, corresponding to the average stress of 0.98 MPa, 1.96 MPa, and 2.94 MPa. As shown in Figure 14a, the degree of yarn squeezing increases with the tensioning force and vice versa for the total width. The twisted yarn has a higher resistance to cross-section deformation than the straight yarn [35]. Using the same tension force, three models, containing 55 twist/m, 110 twist/m, and no twist, were generated, the results of which are compared in Figure 14b. The ratio of section deformation degree in the twist condition is smaller than that in the no-twist situation, while the total width exhibits the reverse relationship. The twisting

The Tension Force of Carriers and the Twisting of Yarns
The 3D braiding process is fully automated and does not require manual jamming. The yarn carriers can control the yarn tension force at a constant value by retracting or releasing the yarns. Three simulations were performed with an increasing tension force of 0.5 N, 1 N, and 3 N, respectively, corresponding to the average stress of 0.98 MPa, 1.96 MPa, and 2.94 MPa. As shown in Figure 14a, the degree of yarn squeezing increases with the tensioning force and vice versa for the total width. specific tension force (see Figure 12a). A smaller bending modulus of yarns, whether surface yarn, corner yarn, or interior yarn, leads to lower resistance of the yarn spreading to transverse deformation (Figure 13b), resulting in a larger ratio . The total width variation ( ) increases as the bending stiffness decreases at the same yarn tension force, resulting in a higher yarn volume fraction.

The Tension Force of Carriers and the Twisting of Yarns
The 3D braiding process is fully automated and does not require manual jamming. The yarn carriers can control the yarn tension force at a constant value by retracting or releasing the yarns. Three simulations were performed with an increasing tension force of 0.5 N, 1 N, and 3 N, respectively, corresponding to the average stress of 0.98 MPa, 1.96 MPa, and 2.94 MPa. As shown in Figure 14a, the degree of yarn squeezing increases with the tensioning force and vice versa for the total width. The twisted yarn has a higher resistance to cross-section deformation than the straight yarn [35]. Using the same tension force, three models, containing 55 twist/m, 110 twist/m, and no twist, were generated, the results of which are compared in Figure 14b. The ratio of section deformation degree in the twist condition is smaller than that in the no-twist situation, while the total width exhibits the reverse relationship. The twisting The twisted yarn has a higher resistance to cross-section deformation than the straight yarn [35]. Using the same tension force, three models, containing 55 twist/m, 110 twist/m, and no twist, were generated, the results of which are compared in Figure 14b. The ratio of section deformation degree in the twist condition is smaller than that in the no-twist situation, while the total width exhibits the reverse relationship. The twisting process increases the lateral stiffness of the yarn, and a suitable twist can reduce fiber pilling during the braiding process. However, too much twist increases the shear force among the fiber filaments, which is not conducive to the overall tensile strength. Therefore, it is necessary to consider yarn twisting in simulation.

Application
The 3D braiding composite shafts are widely used in aerospace, wind power, and automotive transmission systems. These shafts can be integrally braided using the cylindrical braiding machine shown in Figure 15a [36,37]. The yarn arrangement consists of three layers in the radial direction, 18 columns in the circumferential direction, and the side yarns. The yarn carriers move alternately in one position in radial and circumferential directions. After four steps, the take-up device lifts one pitch height upward correspondingly. Figure 15b shows the simulation result of the braiding process, and Figure 15c shows the result of yarn deformation on the mesoscale. The yarns are arranged periodically in the circumferential direction. Figure 15d-f shows the RVE, and the yarn cross-section changes from a circle to a complex shape. Therefore, the application in the 3D cylindrical braiding process demonstrates the versatility of the simulation methodology.

Application
The 3D braiding composite shafts are widely used in aerospace, wind power, and automotive transmission systems. These shafts can be integrally braided using the cylindrical braiding machine shown in Figure 15a [36,37]. The yarn arrangement consists of three layers in the radial direction, 18 columns in the circumferential direction, and the side yarns. The yarn carriers move alternately in one position in radial and circumferential directions. After four steps, the take-up device lifts one pitch height upward correspondingly. Figure 15b shows the simulation result of the braiding process, and Figure 15c shows the result of yarn deformation on the mesoscale. The yarns are arranged periodically in the circumferential direction. Figure 15d-f shows the RVE, and the yarn cross-section changes from a circle to a complex shape. Therefore, the application in the 3D cylindrical braiding process demonstrates the versatility of the simulation methodology.

Conclusions
The present research developed a multiscale simulation methodology based on the virtual sub-yarns to establish accurate RVEs. The simulation results were compared with the micro-CT and SEM structures in the literature. A parametric study was conducted to validate the model. The multiscale method yields accurate results with a relatively small computational cost. Furthermore, the method could give some predictive guidance on the process design in the early stage. The following conclusions can be drawn: (1) The braiding angle decreases as the height of the braiding point increases and vice

Conclusions
The present research developed a multiscale simulation methodology based on the virtual sub-yarns to establish accurate RVEs. The simulation results were compared with the micro-CT and SEM structures in the literature. A parametric study was conducted to validate the model. The multiscale method yields accurate results with a relatively small computational cost. Furthermore, the method could give some predictive guidance on the process design in the early stage. The following conclusions can be drawn: (1) The braiding angle decreases as the height of the braiding point increases and vice versa for pitch length. (2) To achieve a realistic interior cell, surface cell, and corner cell, a temperature drop simulation with hybrid elements was conducted to make the yarn cross-section deformable. The simulation result shows good agreement with the experiment. (3) The parametric study shows that the use of 92 sub-yarns represents the cross-sectional deformation well. A smaller bending stiffness leads to a larger total width variation, resulting in a higher yarn volume fraction at the same tension force. The yarn deformation degree is positively related to the tension force of carriers. Finally, the twisted yarn has a higher resistance to cross-sectional deformation than the straight yarn. (4) The simulation methodology presented here is a versatile tool that can study the relationship between the different braiding processes and the mesostructure of the preform.
This paper proposes a numerical method to obtain the accurate RVEs. Future works will integrate the source codes (see Supplementary Materials) into Abaqus Plug-ins to further accelerate engineers' research and development work.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/polym14194210/s1. Video S1: the simulation results of the 3D4D braiding process, Video S2: the sub-yarn deformation, and the source codes of our workflow.