A multi-scale method for predicting the properties of 3D braided composite under three-point bending load

A more versatile and efficient multi-scale coupling finite element method for researching the mechanical response of 3D braided composites under three-point bending load is represented in this paper. In the mesoscale, the multiphase representative unit-cell models are established to describe the mesoscopic structure which consists of braiding yarns and matrix. In the macroscale, the unit-cells are regarded as homogeneous material, and the load and constraint conditions are applied on the macroscopic structure model. The multi-scale homogenization theory is introduced to calculate the equivalent stiffness matrixes of mesoscopic unit-cells and build the mathematical relationships between the mesoscopic stress fields and the macroscopic strain fields. According to the element damage criterion, the bending modulus and ultimate load-bearing ability of 3D braided composites are predicted by simulating the progressive damage process of unit-cells Comparing with the experimental result, the predicted result satisfies the required precision for engineering.


Introduction
Comparing with traditional metallic materials and laminated composites, 3D braided composites have exhibited many outstanding advantages, because they comprehensively utilize the characteristics of different components and display better structural integrity. Therefore, 3D braided composites have extensive application in aeronautical, professional racing, sports equipment, medical apparatus and instrument fields, etc. However, it is difficult to accurately calculate the mechanical response and establish the failure criteria for 3D braided composites under various kinds of loads and supports because of the complex internal geometries characteristics. Additionally, the experimental analyses waste lots of time and money due to the variety of braiding parameters which influence the macroscopic and mesoscopic structure of composites. For purpose of taking full advantage of the 3D braided composites reasonably, it is necessary to adequately analyze the physical properties and develop an effective method to investigate the mechanical behavior and failure mechanisms.
It is the precondition for analysing the bending behavior and properties of 3D braided composites that establishing a reasonable mesoscopic structure model and predicting the elastic constants. In the past, Zeng [1] and Jiang [2,3] et al respectively proposed finite element method to predict the effective moduli and the local stress based on the helix geometry unit-cell model. On account of three unit-cell model, Wang [4], Tang [5,6], Chen [7] and Shokrieh [8] and their co-workers presented a series of studies about 3D braided composites which include fabric structure, prediction of mechanical properties, finite element analysis and simulation of deformations. However, these interior unit-cell models were not completely periodic and the surface/corner unit-cell models were not rectangular. Therefore, multi-unit cell models which changed the partition scheme were developed, and the cross section of the fiber bundle was respectively assumed to be an ellipse [9], a hexagon contains an inscribed ellipse [10] or an octagonal contains an inscribed ellipse [11]. But the complicated mesoscopic structure makes it be difficult to forecast the mesoscopic stress field of 3D braided composites. In order to ensure the accuracy, a mesostructure model was established by Wu and co-workers [12,13]  paid on the multi-scale asymptotic expansion homogenization method which studied the mechanical response by coupling macroscale and mesoscale [14][15][16][17]. Most of these methods were based on the finite element simulation of representative volume element (RVE) models. In addition, Sun [18] and Fang [19] et al researched the three-point bending fatigue behavior of 3D braided composites in experimental and finite element analysis approaches.
The present work is aimed at establishing a more versatile and efficient multi-scale finite element method to research the mechanical behaviours of the 3D braided composite beam structure under bending load, in which the multi-scale homogenization method combining with the multiphase unit-cell models are used to calculate the elastic properties of the unit-cells in the mesoscale. This method has taken the difference of the mesostructure between interior region and surface region into account, and the unit-cell models in the mesoscale are from the previous study [20]. Based on the multi-scale homogenization method, the mathematical relationships between the stress fields in the mesoscale and the strain fields in the macroscale are built. Finally, the ultimate load-bearing ability of the 3D braided composite beam structure under bending load is predicted according to the progressive damage of unit-cells in the mesoscale and the complex failure behavior of braided composites in the macroscale. The finite element analysis results are compared with the experimental results for verifying the analytical method. All in all, this method is helpful to the design of the 3D braided composite material.

Multi-scale homogenization method
According to the previous studies, the complete structural model of 3D braided composites is so complex that the calculated quantity of the finite element analysis is huge and the accuracy of calculation cannot be guaranteed. On the other hand, the mesoscopic unit-cell model of 3D braided composites is characterized by periodic structure. Hence, the multi-scale homogenization method is adopted to simplify the finite element analysis process. This method is based on the asymptotic expansion homogenization, and the mesoscopic structure of the composite material is fully taken into account. With using the method, the initial analysis is converted to mesoscopic stress analyses and macroscopic homogenization. Finally, the material properties are evaluated by the analyses of mesoscopic unit-cells.
2.1. The asymptotic expansion homogenization for multi-scale finite element analysis The homogenization theory was firstly developed by Babuška and Lions, and then the theory was introduced into the mechanical analysis of composite materials by other researchers [16,17]. In the asymptotic expansion homogenization method, a dual-scale coordinate system is used to describe the mesoscale and the macroscale. The mesoscopic coordinate y and the macroscopic coordinate x can be correlated by a very small coefficient ε (0<ε=1), and the relational expression is: ε=x/y. If ε equals to zero, the macroscopic structure can be regarded as homogeneous.
Because of the periodicity of mesoscopic structure, the field variable on both macroscale and mesoscale can be expressed as: where Y is the period of the mesoscopic unit-cell, j ε (x) is a scalar, vector or tensor field. When taking the derivative of j ε (x) for x, the relational expression will be formulated as: Displacement fields, strain fields and stress fields are key field variables to analysis the mechanical properties of the unit-cells and the 3D braided composite which is made up of these unit-cells. According to the multi-scale homogenization method, the displacement is assumed to be composed by the smooth global displacement in the macroscale and the local oscillation in the mesoscale, and the mesoscopic displacement is caused by the heterogeneous mesostructure. The displacement field can be written as a series expansion:


In this method, the equilibrium equation of force boundary conditions is expressed as: The equation is still satisfied when ε→0 + , so the following equations hold: According to the past research, equation (7) can be solved and then: That means u i 0 ( ) is independent of the mesoscopic coordinate. Substituting equation (10) into equation (8), the equation will become: In order to simplify the finite element calculation, a characteristic function χ kl (y) is introduced to associate the macroscopic displacement u (0) with the first order mesoscopic displacement u (1) , and the characteristic function can be presented as: Then, the mesoscopic strain and mesoscopic stress will be expressed as: The homogenization effective property coefficients which also mean the stiffness coefficients of the unit-cell are defined as: where |Y| is the volume of unit-cell.

The solution of characteristic correctors
It can be observed that the solution of characteristic correctors χ is the key to solve the mesoscopic strain, the mesoscopic stress and the homogenization coefficients. Substituting equation (12) into equation (11), the equation will be formulated as: Comparing with the equilibrium equation of force boundary conditions, k mn c can be seen as a displacement field which subjected to a body force, and the body force is equal to: The body force f mn is zero inside a single material, and the singularity exists at the boundaries of different materials.
On account of the virtual work principle, equation (16) can be written as a weak form of the integral equation: The first term of the equation can be expanded with applying the integration by parts and Gauss theorem: where n j means the component of the unit outer normal vector for the boundary in direction j.
The second term of equation (19) also can be expanded with applying integration by parts and Gauss theorem: The first term of equation (19) equals to zero because the virtual displacement δv i = 0 on the boundary, and the first term of equation (20) also equals to zero.
Supposing the unit-cell contains p kinds of materials, and they respectively define different regions. The boundaries of these regions are expressed as ∂Y 1 , ∂Y 2 , K, ∂Y p . Summary, equation (18) will be written as: Hence, the solution of χ kl (y) is transformed into a weak form of integrals for linear elasticity problems. For a periodic unit-cell, there are series of internal surface distributed forces acting on the internal boundaries of different material regions, and the component of the interface distributed force in direction i is defined as C ijmn n j .
With the change of mn, the equation will be calculated based on six independent displacement fields. According to the Voigt rules, χis shown as:  F C n C n C n C n C n C n C n C n C n C n C n C n C n C n C n C n C n C n

The model of unit-cells and braided composites
In the previous research, a modified multi unit-cell model was established to describe the mesoscopic structural characteristics of 3D braided composites [20]. In this model, the 3D braided composite was regarded as a shellcore structure, and the mesoscopic structure in the surface region was different from that in the interior region. The interior region of the braided composite was composed of periodic structural volume units, and the periodic units were arranged in order, the representative volume unit was named as interior unit-cell. Similarly, the external region was consisted of surface unit-cells and corner unit-cells.
For the braided composite with small aspect ratio, the structural difference of unit-cells in the mesoscale should not be ignored. Hence, not only the interior unit-cell but also the surface unit-cell and the corner unitcell should be studied. The position of braiding yarns inside different unit-cells are illustrated in figures 1, and 2  shows the mesoscopic structure models of different unit-cells. The unit-cell models are built refer to the past paper [20] which proposed the partition scheme of unit-cells and provided the calculation method of structural parameters. These finite element models of mesoscopic unit-cells are used to calculate the characteristic correctors χ and the homogenization effective property matrix C H . The macroscopic structure model of 3D braided composites is demonstrated in figure 3. The unit-cells are defined as homogenized and transversely isotropic material in macroscale, while the elastic constants of unitcells are dependent on the homogenization effective property matrix. The macroscopic model of 3D braided composites is applied to solve the macroscopic strain fields under external loads.

Boundary conditions of representative unit-cells
The periodic boundary conditions are employed on the periodic surfaces of mesoscopic unit-cells. As mentioned in the predecessors' research [21][22][23][24], the periodic boundary conditions can be expressed as: where u i is the component of displacement field for periodic boundaries of mesoscopic unit-cells in direction i; index 'j+' and 'j−' denotes a pair of periodic outer boundaries. In order to fulfil the periodic boundary conditions in the FE model of mesoscopic unit-cells, the nodes on a pair of periodic surfaces need to be coupled each other. Specially, the coupling group of nodes in the intersection lines of periodic surfaces includes four nodes, and the coupling group of nodes in the intersection point of periodic surfaces includes eight nodes. Moreover, any vertex of unit-cell should be fixed because the rigid body motion of the unit-cell should be avoided. The coupling code is written as APDL commands and then applied to the pre-processing of finite element analysis.
Additionally, displacement continuous conditions between contact surfaces of different unit-cells need to be satisfied. Specifically, for the FE model of the mesoscopic surface unit-cell, the displacement fields of nodes on the contact surface with the interior unit-cell are consistent with that on the contact surface of the interior unitcell. For the FE model of the mesoscopic corner unit-cell, the displacement fields of nodes on the contact surface with the surface unit-cell are consistent with that on the contact surface of the surface unit-cell. The displacement continuous conditions are realized by using polynomial interpolation to match the node displacement fields, or making sure that the node positions correspond to each other. The polynomial interpolation algorithm about nodal displacements is written in Matlab, and then applied to the pre-processing of finite element analysis as nodal displacement constraints which expressed in APDL commands.

Strength criterions and damage
The failure strength of composites is decided by the material property of compositions and the mesoscopic structure, while the 3D braided composites are composed of braiding yarns and matrix. Furthermore, the braiding yarn consists of thousands of fibers and matrix which fill in the interspace between fibers. It can be determined that the braiding yarn is a transversely isotropic brittle material, and the matrix is an isotropic elastic material. Hence, Tsai-Wu second-order tensor polynomial is introduced as the failure criterion for braiding yarns, and Von Mises effective stress yield criterion is carried as the failure criterion for matrix.
Tsai-Wu second-order tensor polynomial is expressed as: where F i and F ij denote the strength tensors of braiding yarns, and they can be calculated from the strength parameters of the components (fiber and matrix). Von Mises effective stress yield criterion is expressed as: where σ I , σ II , σ III are the principal stresses, σ m is the uniaxial tensile strength of matrix. These criterions are used to determine whether the elements of the mesoscopic unit-cell are damaged. The mesoscopic stress fields of mesoscopic unit-cells are dependent on the macroscopic strain fields of macroscopic unit-cells. The mesoscopic unit-cell is thought to be damaged when the entire cross section of any fiber bundle is failed.
In the 3-point bending simulation, the load is applied as nodal forces on the upper surface of the 3D braided composite, and the load nodes are on the midline which is in the middle of two supporting constraints. The load and constraint conditions are shown in figure 4, and the dangerous section of the composites is the middle part of the structure. As applied external load increases, the stiffness matrixes of the damaged macroscopic unit-cells are reduced to very small values. When the damages expand to the entire surface and extend into the interior of the dangerous section, the whole structure is supposed to be destroyed completely.

Analysis sequence
Because the finite element analysis softwares do not contain the multi-scale honogenization method developed above, secondary development programs based on ANSYS and Matlab are applied to establish the finite element analysis models. For ANSYS/Workbench, APDL Command is internally installed and can be used to insert algorithms. Besides, the simulation results are exported to files and the files will be used in the post-processing algorithms which written by Matlab. The detail analysis process of the multi-scale homogenization finite element method for the 3D braided composite beam is provided in figure 5.
The detailed analysis sequence is as follows: (1) The finite element analysis of the representative unit-cell models in the mesoscale. Firstly, the mesoscopic unit-cell models are established based on the partition scheme proposed in the past research, and the 3D 10node tetrahedral elements are employed to mesh the unit-cells in HyperMesh. Secondly, the corresponding nodes on the periodic boundaries of the unit-cells are coupled, and the nodal displacements of the contact surface for the surface unit-cell and the corner unit-cell are defined to match the adjacent interior unit-cell and the adjacent surface unit-cells, respectively. Thirdly, the interface distributed forces are applied in six independent load steps. After the finite element solutions, the equivalent displacement matrix χ and the homogenization effective property coefficients matrix C H can be obtained.
(2) The establishment of the finite element analysis model of the 3D braided composite beam structure under bending load. Firstly, the unit-cells in the macroscale are replaced by the homogeneous cuboid, and the engineering elastic constants of the macroscopic unit-cell are from the equivalent stiffness matrix which calculated in the mesoscopic finite element analysis process. Secondly, the 3D 20-node rectangular isoparametric elements are employed to mesh the macrostructure, and the grid density at constraint and load location is higher than others. Thirdly, the macrostructure is constrained as a simply supported beam, and the external load is applied in the middle line of upper surface. The constraint conditions and load are shown in figure 4. In this work, the macroscopic displacement field and strain field under bending load are calculated.
(3) The damage criterion of the mesoscopic unit-cells. The mesoscopic stress field of the unit-cell can be calculated according to the multi-scale homogenization method, and the strength criterions of the braiding yarns and matrix are provided above. The whole unit-cell is thought to be damaged when the entire cross section of any fiber bundle is failed.
(4) The failure criteria of the macroscopic structure. Firstly, the stiffness reduction is applied to the damaged macroscopic unit-cells, while the reduction coefficient is 0.99. Secondly, the mechanical response of the macroscopic structure is updated, then the damage of the mesoscopic unit-cells is judged again. This procedure is repeated unit no unit-cell is damaged. If the dangerous section of the macroscopic structure has not destroyed completely, the external load will be increased. When a segment of the macroscopic structure elements is failed, the 3D braided composite beam is broken, and the load is the ultimate load.

Experimental result
The 3D four-directional braided composites with four-step 1×1 rectangular braiding process were investigated in this work, and the experimental data were obtained from Zhen's research [25]. The braided preform was fabricated by T300-12 K carbon fiber tows, then it was impregnated and consolidated with TDE-85 epoxy resin by using Resin Transfer Moulding. The braiding array was 23×4, the braiding angle was 20°, and the presupposed fiber volume fraction was 45%. Three specimens with same parameters were produced and they were cut to 250 mm×25 mm×5 mm. The three-point bending test were performed on a 100 kN test machine (MTS 810) with a constant velocity of 2 mm min −1 . The specimen was placed on two supporting rollers with spacing of 160 mm. The elastic properties of the carbon fiber and epoxy resin are given in table 1.
The bending modulus and bending strength can be calculated respectively as: Table 1. Elastic properties of carbon fiber and matrix.
where, b and h is the width and thickness of the specimen; L is the spacing of two supporting rollers; F is the ultimate bending load; ΔF and Δf is the increment of the load and maximum deflection.
The experimental results are given in table 2.

Homogenization results
The mesoscopic finite element analysis is the basis for preditcing the mechanical response of the macroscopic finite element model under three-point bending load. The equivalent displacement fields of unit-cells under different load-step are simulated in the finite element analysis processes. Figure 6 shows the distribution of χ 11 under interface load F b11 . It is obvious that the boundaries of the unit-cell are periodic and continuity. According to the simulate results under six step of finite element processes and the computational formula of homogenization effective property coefficients which given by equation (15), the equivalent stiffness matrixes of the unit-cells are calculated, and the engineering elasticity coefficients calculated from the compliance matrix are shown in table 3.

Results comparison and discussion
The mechanical response of the 3D braided composite under three-point bending load is simulated on the macroscopic structure model, and the structural deformation is shown in figure 7. With the load increases, the deflection of the 3D braided composite beam increases in linear. The load-deflection curves of simulation and experiment results are shown in figure 8. According to the computational formula given in equation (27), the bending modulus of the 3D braided composite without damaged was calculated, and the simulated result was 82.742 GPa. Comparing with the experimental result, the simulation results of the bending modulus was 7.04% higher. The deviation is because the stiffness matrixes simulated by homogenization method are usually higher than the actual values. The bending strength of the 3D braided composite was predicted according to the analysis sequence shown in figure 5, and the dangerous section of the composite was the middle part of the structure. When the bending   load was 4200 N, one of the fiber bundle in the mesoscopic corner unit-cell was completely fracture, so the mesoscopic corner unit-cell was deemed to have failed. At this monent, 16.71% elements of all fiber bundles were damaged according to the strength criterions, but the matrix elements were almost no damage because the fiber bundles were carrying most of the load for 3D braided composite with small braiding angle. Then the stiffness matrix of the macroscopic corner unit-cell is reduced by 99% and the macroscopic strain fileds are recalculated. According to the new macroscopic simulation result, other mesoscopic unit-cells were not failure. When the bending load reached 4300 N, the mesoscopic surface unit-cell near the damaged corner unit-cell failed because a completely broken fiber bundle appeared. After several recalculations, allsurface unit-cells of the middle part failed. When all the surface unit-cells of the dangerous section damaged, the 3D braided composite beam was considered to be completely ineffective. While the ultimate bending load was 4300 N, and figure 9 shows the change in the proportion of the failure elements. According to the computational formula given in equation (28), the bending strength of the 3D braided composite was 949 MPa, and the simulated result was 2.15% less than the experimental result. It can be found that the 3D braided composite under ultimate bending load is brittle rupture, that is because the braiding angle is small. So, the fiber bundles bear most of the load, and the main reason to lead the failure of the 3D braided composite is the fiber brittle fracture.

Conclusion
In this paper, the multi-scale finite element analysis method was developed to predict the bending modulus and bending strength of 3D braided composites. In this method, the multi-scale homogenization method and multi unit-cells model were used to calculate the equivalent stiffness matrix of unit-cells in the mesoscale. Additionally, the multi-scale homogenization method was applied to establish the mathematical relationships between the stress field of mesoscopic heterogeneous unit-cells and the strain field of macroscopic homogeneous unit-cells. Based on the Tsai-Wu second-order tensor polynomial strength criterion, the damage of fiber bundles in the mesoscopic unit-cells was judged. Then the progressive damage process of unit-cells as load increases was simulated. The bending modulus and bending strength were predicted according to the simulated mechanical response of the 3D braided composite under three-point bending load. The experimental results demonstrated that the prediction model was receivable.