A novel finite element method for predicting the tensile properties of 3D braided composites

A novel unit-cell model which considers the structure interference of braiding yarns in the external region is proposed to describe the geometric structure of 3D braided composites. Based on the multi-scale homogenization theory, mesoscopic structure finite element analysis model is established to predict the stiffness matrix of unit-cells. Periodic boundary conditions are applied on the periodic surfaces of unit-cells to reflect the repeating feature. And the node deformation matching equations are applied on the internal connected faces of the surface unit-cell and the corner-cell, so that the displacement continuity conditions and traction continuity conditions are satisfied. In allusion to 3D braided composites with small braiding angle, principal tensile strength forecasting model is developed.


Introduction
3D braided composites are a kind of carbon fiber reinforced composites (CFRC). The 3D four-directional braided preforms which manufactured by four-step 1×1 rectangular braiding process are impregnated and consolidated with matrix materials. Because of the excellent tensile properties, outstanding fatigue resistance, higher damage tolerance and specific energy absorption behavior, CFRC are widely used in industry [1]. Comparing to laminate composites, the 3D braided composites present better structural integrity and interlaminar shear characteristics.
In previous research, representative volume element (RVE) model was provedreliable to study the structural characteristics of 3D braided composites. Chen et al firstly introduced three types of microstructural unit-cell model to analyse the structural geometry of 3D braided composites and derive the mathematical relationships among the structural parameters [2]. The effective elastic properties were predicted based on the three unit-cells model [3,4]. About the same time, Tang et al developed the finite element analysis and simulation of deformation for the prediction of mechanical properties [5,6]. Zhang, Xu and Li et al presented a new multi-unit cell model to describe the skin-core structure of 3D braided composites, however different shape of cross section of fiber yarns were assumed [7][8][9]. Xu and Fang el al. established a multiunit cell model to study the influences of yarn distortion on the stiffness and strength of the braided composites [10,11]. Dong and Huo analysed the effect of pore defect on the mechanical properties of 3D braided composites [12]. Zhang et al established a mesoscale finite element model which considered the interface to numerically predict the stiffness and strength properties of 3D braided composites [13]. Hao and He et al separately proposed damage model to investigate the fatigue behavior of 3D braided composites [14,15]. Wang et al studied the progressive damage process of 3D braided composites [16,17]. This paper is aimed at developing a new analytical method to predict the elastic properties and principal tensile strength of 3D braided composites. A novel parameterization unit-cell model, which redefines the mesoscopic unit-cell structure in the external region, is established to describe the structural characteristics, and multi-scale homogenization method is applied to calculate the stiffness matrix of unit-cells. The finite element model is also used to predict the tensile strength of 3D braided composites. Through the experimental, the simulation results will be verified.

Mesoscale multi-unit cell structural models
As the braided structure of 3D four-directional braided composites is extremely complex, it's difficult to solve the finite element model of entire structures. Fortunately, the mesoscale structure of the braided composites is of characteristic periodic features. Thus, the representative volume element named as unit-cell, which is the smallest periodic structure composing the braided composites, is proposed to be researched and analysed instead of the whole structure.

Parameterized unit-cell models
Based on the multi-unit cell model developed by Xu [8], the periodic mesostructure of braided composites consists of three types of unit-cells, i.e., the interior unit-cell, the surface unit-cell and the corner unit-cell. Furthermore, the braiding yarns of surface unit-cells and corner unit-cells in different regions of composites take on different configuration. More specifically, the yarn orientation of surface unit-cells in the opposite surfaces of composites takes on antisymmetric permutation, as well as that of corner unit-cells in the catercorner of composites. Moreover, the yarn orientation of surface unit-cells in the adjacent surfaces of composites takes on mirror symmetry permutation, as well as that of corner unit-cells on the same side. Due to these features, both the surface unit-cells and the corner unit-cells are summarized as two kinds of mesostructure.
The periodical distribution of unit-cells is demonstrated in figure 1, W x , T y denotes the width and thickness of composites, respectively. The minimum period of unit-cells along the extension orientation is defined as the pitch height of h.
The topology geometrical models of unit-cells which illustrate the topology realtionships of braiding yarns are illustrated in figure 2. W, T identifies the width and thickness of unit-cells, respectively. While the sub index i, s, c denotes the interior, surface and corner unit-cell, respectively. The angle between the central axis of the interior fiber yarn and the extension orientation is defined as the interior braiding angle of γ.
The establishment of the multi-unit cell model is based on the microstructure configuration of the interior braiding yarns. The SEM micrograph of the cross-section morphology in the preform interior cut longitudinally at a 45°angle with the preform surface was given by Chen [2]. It showed that the interior braiding yarns axes remain straight and the yarn faces contact with each other. In order to embody the distortion as mutual squeezing, the cross-section of interior braiding yarns is assumed as a circumscribed octagon of an ellipse with major radius and minor radius, a and b, respectively, as shown in figure 3.
In the external region of braided composites, the braiding yarns are repositioned by the jamming action, and the path of the yarns are retraced abruptly. So, the braiding yarns axes are distorted and present as complex spatial curves. For simplicity, the yarns axes in the external region is assumed to consist of several similar partial elliptic arcs, as shown in figure 4. Meanwhile, the cross-section of the yarns is irregular and variable because of the distortion of the axes. As shown in figure 3, the cross-section of surface braiding yarns and corner braiding yarns is assumed as a circumscribed regular octagon of a circle with radius, b. That assumption not only embodies the distortion as mutual squeezing, but also minimize the mutual intrusion of the yarns. For simplicity, the change of cross-section from the external region to the interior region is ignored in this model.
According to the hypothesis presented above, the cross-sectional dimension of 3D braided composites and unit-cells can be expressed as: where m, n is the number of yarn carriers used on the machine bed in the rows and columns, respectively. As shown in figure 2(a), the interior braiding angle can be formulated as: Considering the tangent relationship of braiding yarns in the interior unit-cell, the cross-sectional dimension of interior braiding yarns can be obtained: And the cross-sectional area of the interior the yarns is calculated as: The cross-sectional dimension of fiber bundle in the external region can be expressed as: And the cross-sectional area of the yarns is written as: According to the design formulas presented above, all other structure parameters of unit-cells are decided by b and h. With these two basic design parameters, the parameterized structural models of unit-cells will be established in SOLIDWORKS, as displayed in figure 5.

Material properties of constituents
The braiding yarns are composed of thousands of carbon fibers and matrix which fill in the intervals of fibers. The constitutive relationship of fiber is believed to be transversely isotropic and linear elastic, and the constitutive relationship of matrix is assumed isotropic and linear elastic under the assumption of small strain. Consequently, the braiding yarns are modeled as transverse isotropy and unidirectional fiber-reinforced composites in local coordinate system. In the local coordinate of a yarn, local 1-axis follows the yarn path line and local 3-axis is in the upright plane of 1-axis and parallel to the cross-section of composites. Without consideration of the twist of braiding yarns, the elastic constants of braiding yarns can be calculated via the fiber-matrix rule of mixture proposed by Chamis [18], which are given by: where ε is the yarn packing factor. The braiding yarns are regarded as transversely isotropic in elastic.
According to the prediction model presented by Kelly [19], the fibres are assumed to have same strength, and the failure criterion of yarns is decided by the fiber. Based on the assumption that the tensile strain is same between fibres and matrix, the tensile strength of yarns can be computed in the following form: where σ fmax is the maximum tensile stress of fiber, (σ m ) ε fmax is the matrix stress when the matrix stress equals to maximum tensile strain of fiber.

Multi-scale finite element method
As noted earlier, the microstructure of the braided composites is heterogeneous and the finite element (FE) model of entire braided composites is rather difficult to be solved, so the unit-cells are introduced as an alternative. In previous research, the FE models of unit-cells were mutually independent, and the deformation compatibility between two different types of unit-cells was simplified. Additionally, the FE models were not suitable for predicting the tensile strength of the braided composites, because the mesoscopic stress could not be accurately calculated.
In order to make up those inadequate aspects, the multi-scale homogenization method is introduced, which associates the mesoscopic stress with the macroscopic strain. In addition to impose the periodic boundary conditions at the periodic surface, the internal connected faces of surface unit-cells or corner unit-cells must satisfy the displacement continuity condition and traction continuity condition. Moreover, the stiffness matrix of unit-cell is equivalent to the homogenization coefficient, which composed by the average elasticity tensor and the average additional stress that caused by the difference of materials.

Multi-scale homogenization method
The homogenization theory was firstly developed by Babuška and Lions, and then the theory was introduced into the modelling and analysis of composites. According to the homogenization theory, the analysis for properties of composites is transformed into mesoscopic structure analysis and macroscopic homogenization process. In other words, the unit-cells are regarded as homogeneous materials on the macroscopic scale, while the elastic properties of unit-cells are decided by the mesoscale model of unit-cells.
A dual-scale coordinates, which includes macroscopic coordinate x and mesoscopic coordinate y, is introduced in the method. And the ratio of macroscopic unit length to mesoscopic unit length is defined as ε (0<ε«1), which can be expressed as: y+nY=x/ε (n=0, 1, 2, K). Because of the heterogeneity of composites, the displacements, strains and stresses in mesoscale are different from those in macroscale. And these parametric functions are simultaneously depended on the macroscopic coordinates and the mesoscopic coordinates. Considering the periodicity of mesoscale structure, these functions are Y-periodic on the mesoscopic coordinates, where Y is the period of unit-cells. Summary, the parametric functions are written as: Then taking the derivative of j ε (x) with respect to x, the relationship is formulated as: The displacement function can be expanded incrementally with ε: According to the derivation rule given by equation (12), the strain tensor is calculated as: A characteristic function χ kl (y) is introduced to associate the macroscopic displacement u (0) with the 1st order mesoscopic displacement u (1) , which is expressed as: In summary, the mesoscopic stress will be given by: where |Y| is the volume of unit-cell.

Boundary conditions
As metioned earlier in the literature, the 3D braided composites are composed of periodic unit-cells. The unitcells are considered to have same deformation pattern, and there is no separation or overlap between adjacent unit-cells. In order to satisfy the displacement continuity condition and traction continuity condition, the periodic boundary conditions must be imposed to the periodic surfacesof unit-cells. For interior unit-cell, all three pairs of opposite surfaces are periodic. However, that the surface unit-cell has periodicty only in two pairs of opposite surfaces without the outside surface and it's opposite, and the opposite surfaces of corner unit-cell is periodical just in direction of extension.
According to the previous work [20][21][22], the difference of displacement field between a pair of periodic surfaces is: where u i is the displacement of periodic surface; index 'j+' and 'j-' means along the positive X j direction and the negative X j direction, respectively; x k is the Cartesian coordinate of a unit-cell point; e ik is the average strain of the unit-cell.
As !x j k is a constant for each pair of periodic surfaces, once e ik is specified, the equation can be expressed as: By setting the displacement linear constraint equations to the corresponding nodes on the periodic surfaces of the unit-cell, the equation (26) will be carried out, and the displacement continuity conditions are satisfied. It was proved that the periodic boundary conditions also ensured the traction continuity condition due to the minimum total potential energy principle in the FE model of the unit-cell [23].
Besides the periodic boundary conditions applied on the periodic surfaces of the unit-cell, the boundary conditions must be defined on the external surfaces and the internal connected faces of surface unit-cells and corner unit-cells. For the surface unit-cell, the internal surface has contact with the interior unit-cell. And the internal surfaces of the corner unit-cell have contact with the surface unit-cell. The node displacement matching equations are applied to the internal connected faces, and that ensure the displacement continuity condition and traction continuity condition satisfied between two different types of unit-cells. Additionally, the external surfaces of the surface unit-cell and corner unit-cell are fixed without the direction of strain applied in the FE analysis of the unit-cell.

Finite element model
The parameterized unit-cell models are imported into the HyperMesh so that the finite element models are established, and then the models are solved in the ANSYS/Workbench. Obviously, the structures of yarns and matrix in mesoscopic unit-cell models are irregular, so 3D tetrahedral element (SOLID187) is adopted for the discretization of unit-cells. On the macro-scale, the unit-cells are regarded as homogenized material, and the 3D braided composites are composed of these homogeneous structures. 3D hexahedral element (SOLID186) is adopted for the discretization of 3D braided composites. In order to guarantee the accuracy of analysis, the mesh size should be as small as possible. However, the solution of finite element model will spend more time and increase the computation cost at the same time. According to a compromise consideration, a suitable mesh size is selected for the discretization of model.
Before the analysis, the boundary conditions referred above must be applied to the FE models. For the periodic surfaces, the periodic boundary conditions can be applied by setting constraint equations to the opposite surfaces. For the internal connected faces, the node dispalcements are applied by inserting APDL commands, and the node dispalcements are calculated by the mathcing algorithm written in MATLAB. For the external surfaces, the boundary conditions are performed by defining the displacements in different directions., which will be set to 0 or free.
According to equation (24), the characteristic function χ kl (y) should be solved before calculating the stiffness matrix of the unit-cell. As shows in equation (22), χ kl (y) is decided by mesoscopic displacements and macroscopic strains. So, it can be solved by obtaining the displacements of elements under six types of strain loading which are mutually independent. Herein, the strain loading is realized by applying external displacement to the relevant boundary surfaces.
In order to reduce the influence of stress concentration, the length of braided composites model is 3 h, and only the simulation result of middle part is adopted. Based on the imitation results, the unit-cells can be regarded as transverse isotropy, and the elastic constants can be confirmed through the flexibility matrix of unit-cells. As illustrated in equation (23), the mesoscopic stresses of elements rely on the macroscopic strains, which are simulated by setting corresponding external displacement loadings to the 3D braided composites model.

Prediction of elastic properties and strength
The stiffness matrix of composites can be calculated as: Where C i , C s1 , C s2 , C c1 , C c2 is the stiffness matrix of interior, surface, corner unit-cell, respectively; V i , V s1 , V s2 , V c1 , V c2 is the volume of interior, surface, corner unit-cells, respectively; V is the volume of composites.
For braided composites with small braiding angle, the main failure mode is the brittle fracture of braiding yarns, and the relationship between stress and strain is almost linear. It can be proved that the resin matrix is still linear elastic before the braiding yarns broken. Hence, the displacement increment is applied on the end surfaces of composites. Then, the equivalent principal stresses of braiding yarns are calculated according to the simulation results. It is assumed that the composites threaten to break when very few elements of yarns are damaged. Due to the interior unit-cells account for a significant proportion of the braided composites, the fracture of braiding yarns in interior unit-cells is identified as failure criterion of composites. The percentage of breaking elongation are determined by the corresponding displacement load. Then, the ultimate tensile strength of composites can be calculated based on the constitutive equations.

Experimental results
A piece of 20×6 three-dimensional four-directional braided composites was manufactured to verify the effectiveness of the multi-scale finite element method. T300-3k was selected as the carbon fiber yarns, and the material of matrix was LT-5080 epoxy resin. And the elastic properties of the component materials are listed in table 1. The sectional dimension of sample is 11.059 mm×3.568 mm, and the pitch length equals to 3 mm. The tension test was performed on an MTS 810 test systems at a cross-head speed of 5 mm min −1 . Both the ends of specimens were reinforced by aluminium flakes, and it can avoid the broken of specimens which generated due to the clamping force of testing machine.
According to the test results, the principal tensile modulus of composites was 61 GPa, and the principal tensile strength was 486 MPa, and the corresponding breaking elongation was 0.787%. Figure 6 shows the stressstrain curves of longitudinal tensile property.

Simulation results
According to the finite element analysis model, the engineering elasticity coefficients of unit-cells and composites are predicted, and the results are shown in table 2. Comparing with the experimental result, the simulation results of principal tensile moduli is 4.72% less.
For the finite element model of composites, when the displacement load reaches to 0.072 mm, which means the uniaxial tensile strain equals to 0.8%, the failure elements of fibres appear in the interior unit-cell. It can be calculated that the principal tensile strength of composites is 464.97 MPa, which is 4.33% less than the experimental result.

Comparison and discussion
It can be found that the simulation results are smaller than the experimental results. For the predict results of principal tensile moduli, that is because the surface unit-cell model and corner unit-cell model have not shown the distortion of braiding yarns perfectly. In fact, the structural stiffness is higher in the external region because of the extrusion and deformation of braiding yarns. On the other hand, the boundary conditions of internal contract faces of the surface unit-cell and the corner-cell are specially setted, which effectively reduce the error of simulation results. For the prediction result of the principal tensile strength of the composites, the error comparing with experimental result is due to the failure criterion selected. Although the fiber damage has happened in the interior unit-cell, the proportion of damage elements in the surface unit-cell and the corner unit-cell is considerably less. Because of the friction between the fibres, the structure will not break immediately when minor fibres damaged.

Conclusion
In this paper, a unique unit-cell structure model was adopted to research the 3D braided composites. The braiding pitch height h and the interior braiding angle γ were chosen to be the key design parameters of composites. The assumed yarn path in surface and corner unit-cell greatly decreased the structure interference of braiding yarns under the hypothesis that the cross section of braiding yarns was a regular octagon. Then, multi-scale homogenization method was applied to established the prediction model of elastic properties and principal tensile strength of composites. Based on the finite element analysis model, the mesoscale strain or stress of any element of composites could be calculated. The experimental results demonstrated that the prediction model was receivable.

Acknowledgments
The work was supported by National Key R&D Program of China (No. 2016YFC0600805).