A micromechanics based elasto-plastic damage model for unidirectional composites under off-axis tensile loads

Nonlinear properties of composite materials are essential for their engineering application. In this work, a three-phase micromechanics bridging model is employed to evaluate the nonlinear behavior of a composite from properties of fiber, matrix and interphase. It is assumed that the matrix elastoplasticity and the interface damage are two major sources of the nonlinearity. The former is described by the J2 flow rule. The latter is approximated by an interphase with stiffness degradation. For an interphase, an equivalent damage stress is introduced to account for the effect of normal and shear stress on the interface damage growth. Further, an explicit empirical equation is developed to relate the equivalent damage stress and the stiffness degradation of an interphase. The present elasto-plastic damage model is validated by comparing with experimental data of a series of composites under off-axis tensile loads.

the nonlinearity of a composite. However, models of Vogler et al. 18 , and Camanho et al. 19 are established at meso-scale. Due to the nature of meso-scale models, no physical mechanism at micro scale can be revealed and the experimental cost is relative higher compared to micromechanics models. Some other researchers analyzed nonlinear behavior of composites at micro-scale, where a kind of plasticity, e.g. J2 flow rule, was applied for matrix elasto-plasticity and a cohesive element method is employed to simulate interface debonding. Examples can be found in the work of Melro et al. 20 , Han et al. 21 , and Pulungan et al. 22 However, such models has to be implemented in a finite element model where no explicit expressions are provided. In addition, the computational cost of a cohesive element method is high as reported 13,23,24 . Hiremath et al. 25 presented an analytical elasto-plastic damage micromechanics model for the stiffness reduction of a composite. Nevertheless, their work mainly focuses on stiffness degradation along fiber direction. It is still dubious whether the model is applicable for a composite under a complex load case. Zhao et al. 26 proposed a plastic damage model for a composite at micro scale, where the damage of matrix is simulated by the variation of micro-crack density, and the plastic deformation is described by frictional sliding along micro-cracks. However, the interface damage is not accounted for in their model. Ju et al. 27 developed a micromechanical elastoplastic damage model for fibrous composite. In their work, an orthotropic fiber is used to approximate a fiber with debonding interface. Ju et al. 's 27 work didn't consider the effect of shear stress on interface damage and is not applicable for a composite with an interphase.
In this work, an analytical micromechanics model is developed for the nonlinear analysis of unidirectional composites. Both the matrix elastoplasticity and interface damage are involved in this model. Firstly, the matrix elasto-plasticity is described by a J2 flow rule. Then, the interface damage is simulated by an interphase with stiffness degradation. A Mohr-Coulomb like criteria is presented to account for the effect of normal and shear stresses on the stiffness degradation of an interphase. Finally, an analytical three-phase bridging model 28 is employed to predict mechanical properties of a composite from properties of fiber, matrix, and interphase. There are three advantages in this work. First, the model is computational efficiency owing to its explicit and analytical expressions. Second, analytical empirical equations are provided for the interphase degradation, where the normalshear stress coupling effect on the interface damage is considered. Last, the present model is not only applicable for a composite with interface damage but also a composite with an interphase, such as fiber coatings. This paper is organized as follows. Section 2 establishes a nonlinear model, in which the bridging model, interphase model, matrix elasto-plasticity, and failure criteria are introduced. An experimental validation is made in section 3. Conclusions are drawn in Section 4.

Theory
Before introduce the bridging model in detail, proper nouns involved should be specified. The two-phase composite denotes a composite consisted of a fiber and a matrix. The three-phase composite refers to a composite containing an interphase in addition to a fiber and a matrix, as shown in Fig. 1(a). The inner two-phase composite represents a local structure made of a fiber and an interphase, as shown in Fig. 1(b). The global two-phase composite means a composite made of a matrix and an equivalent fiber, as shown in Fig. 1(c). The inner two-phase composite and the equivalent fiber are identical essentially but have different names for the convenience of presentation in different circumstances. Additionally, the noun "matrix" without an attributive adjective means the constituent material in a composite. In other cases, the noun "matrix" with a specified adjective, such as bridging matrix, stiffness matrix, means a mathematical rectangular array. Moreover, parameters with superscripts f, c, ef, and m denote quantities of a fiber, an interphase, an equivalent fiber, and a matrix. The notation description works throughout the paper.
Two-phase bridging model. The original two-phase bridging model 29 is established with a perfect interface assumption. Based on the bridging model, the effective compliance matrix of a unidirectional composite is given by Eq. (1). In a unidirectional fibrous composite laminate, direction "1" is parallel to a fiber. Direction "2" is in-plane and vertical to a fiber. Direction "3" is out-of-plane and through thickness. In eq. (2), E 11 and E 22 are longitudinal and transverse Young's moduli of a composite, respectively. G 12 and G 23 are in-plane and out-of-plane shear moduli, respectively. υ 12 and υ 23 are in-plane and out-of-plane major Poison's ratios. The compliance matrixes of a fiber and a matrix can be obtained by replacing parameters of a composite in Eq. (2) with their counterparts of a fiber and a matrix.
Wang et al. 30 have proved that a rigorous bridging matrix can be derived from stress fields of a concentrated cylindrical assembly (CCA) model subjected to six different loads. The rigorous format of a bridging matrix is given in Eq. (8).
The elements with superscript "r" denotes parameters derived rigorously. Wang In other words, the effective compliance matrix remains symmetric, although the bridging matrix is non-symmetric. Further, Wang et al. 30 made a comparison between the rigorous bridging matrix and the Huang's version. The conclusion is that the Huang's version can be seen as a simplification of the rigorous one, when a fiber is much stiffer than a matrix in a composite, which is the common case of fiber reinforced epoxy composites.
Please note that Eqs. (3)(4)(5)(6)(7) are valid in elastic range. Based on the tangent linearization theory, mechanical properties of a material can be seen as linear in an infinitesimal load step. Thus, the bridging model can be extended to plastic range. However, the normal and shear stress may be coupled in plastic deformation. Strictly, the instantaneous bridging matrix in a load step should be calculated from an anisotropic Eshelby tensor 31 . However, such operation is too complex to be applied in engineering. Huang 29 reported that the instantaneous bridging matrix can be rewritten as Eq. (9). The diagonal elements in Eq. (9) are given by substituting elastic parameters in Eqs. (3)(4)(5)(6) with tangent instantaneous ones, as shown in Eqs. (10)(11)(12). .
Parameters with superscript or subscript t denote tangent instantaneous quantities. In each load step, the effective compliance matrix of a composite must be symmetric. The off-diagonal elements in Eq. (9) can be obtained by satisfying the symmetric condition shown in Eq. (13).
Formulation (13) provides fifteen equations, which is enough to solve all the fifteen off-diagonal elements in Eq. (9). Equations (3) and (9) are the key of Huang's elastoplastic micromechanics model, which has been validated by many researchers [32][33][34] . In this work, the Huang's two-phase bridging model is employed to establish the three-phase bridging model. three-phase bridging model. The two-phase bridging model is only valid for a composite with perfect interface. For a composite with interface damage, an imperfect interface condition must be considered. It is reported that an interphase model 35,36 can approximate an imperfect interface. Thus, a three-phase CCA model as shown in Fig. 1(a) is employed for a composite with imperfect interface. Wang et al. 28 gave an exact solution for a three-phase bridging matrix. Further, using an equivalent fiber method, a concise three-phase bridging model is available 37 . Figure 1 illustrates the establishment process of a concise three-phase bridging model. Firstly, Take the fiber and interphase as an inner two-phase composite ( Fig. 1(b)), where a two-phase bridging model is applicable. Regard the inner composite as an equivalent fiber. Then, a global two-phase composite ( Fig. 1(c)) consisted of a matrix and an equivalent fiber is generated, where the two-phase bridging model can be applied once again.
Based on the bridging model, an effective compliance matrix of an inner two-phase composite is given by Eq. (14). ij Again, the equivalent fiber and the inner two-phase composite are identical essentially. M ij . For a global two-phase composite, the effective compliance matrix is given as Eq. (20) are the compliance matrixes of a matrix and an equivalent fiber, respectively. A ij II       is a bridging matrix for a global two-phase composite, whose expressions are given in Eqs. (21)(22)(23)(24)(25). II   II  II  II   II   II   II   II   II   11  12  13   22   33 44   (27)(28)(29)(30).
12 is a contract incremental stress tensor of the external load subjected to a composite.
Elasto-plasticity of matrix. The elastoplastic deformation of a matrix is one main source of the total nonlinearity of a composite. J2 flow rule is a well-known plastic flow rule that applicable for both metal and nonmetal materials. When an isotropic work hardening law and the J2 flow rule is used, the instantaneous compliance matrix of a matrix material is given by Eq. (31).
respectively, elastic and plastic parts of the compliance matrix, whose expressions are shown in Eqs. (32)(33)(34).    is the tensile yielding strength at the n th load step. eq m σ is the von-Mises stress of a matrix, whose expression is shown in Eq. (37). , can be obtained by a uniaxial test of a matrix. Alternatively, it can also be obtained from a tangential shear modulus, as shown in Eq. (38).
In plastic deformation, the Poisson's ratio t m υ can be set as 0.5, according to the non-compressive law. The tangential shear modulus, G t m , is given by Eq. (39).
, and G T m n , are, respectively, the shear yielding strength and shear tangential modulus at the n th load step. eq m τ is the octahedral shear stress of a matrix, whose expression is given in Eq. (40). Since the shear nonlinear behavior is more significant than tension, it is better to use a shear stress-strain curve when available.
interphase model. Nonlinear behavior induced by interface damage is more significant than longitudinal or transverse loads for a composite under off-axis tension. Thus, this work focuses on a UD composite subjected to an off-axis tension load. The case of off-axis compression would be investigated in future work. Many researchers 35,36 pointed out that the interface damage can be approximated by an interphase with stiffness degradation. Following their suggestion, a three-phase CCA model is utilized to simulate a composite with interface damage in this work. In general, degraded stiffness of an interphase must be determined experimentally. However, more parameters lead to higher experimental cost. One goal of our model is to reduce the parameter amount in an interphase model. Theoretically, in each infinitesimal load step, an off-axis load subjected to a UD composite can be seen as a combination of longitudinal normal stress, transverse normal stress, and longitudinal shear stress. It is well known that a UD composite is nearly linear elastic under a longitudinal tensile/compressive or transverse tensile load. For transverse compressive load, a UD composite is nonlinear to some extent. However, it majorly results from the matrix plasticity instead of interface damage. Some reports 4,38,39 show that interface damage plays a key role in shear nonlinear deformation. Thus, for off-axis tensile load considered in this work, only the shear modulus of an interphase G c is assumed to depend on the damage states while the Young's modulus E c remains constant.
Following the continuum damage theory, the compliance matrix of an interphase is given as Eq. (38).
where E c and υ c are, respectively, elastic Young's modulus and Poisson ratio of an interphase. G D c is a shear damage modulus whose values are sensitive to the interphase damage state.
Ideally, the shear damage modulus G D c is majorly controlled by shear stress in an interphase. However, experimental results 40 show that transverse normal stress can also affect the evolution of G D c . Approximately, a damaged interface can be simulated by an interphase with micro-cracks, as shown in Fig. 2. The micro-cracks expand or contract when a composite subjected to transverse tension or compression. In other words, transverse tension accelerates while transverse compression retards the micro-crack propagation, thereby affecting the stiffness degradation of an interphase. To consider the interaction between normal and shear stress, the Mohr-Coulomb failure theory is proposed for soil material 41 . Further, it is widely used in analysis of matrix plasticity 42 , interface damage 43 , and composite failure 40 . In this work, a Mohr-Coulomb-like damage criterion is developed as shown in Eqs. (42) and (43) to determine G D c , with which the normal-shear stress coupling effect on the damage evolution of an interphase can be considered. σ will increase the value of f G c , thereby accelerating the damage process. When σ c 22 is negative, the f G c will be reduced, meaning that the damage process is delayed. Thus, Eq. (42) is able to consider the effect of a transverse normal stress on the shear damage process. Equation (43) is the evolution law for the shear modulus of an interphase, where the damage driving stress is employed as a criterion.
It is very difficult to measure the mechanical properties of an interphase in a composite. Alternatively, engineering reverse calculation is a proper way. Figure 3 illustrates the reverse calculation process. Two parameters, G D c and b, need to be reverse calculated in each load step, thereby requiring two experimental stress-strain curves. Thus, step one it to select two experimental stress-strain curves of a UD composite under different off-axis loads.
Step two, determine bridging parameters, α and β. Strictly, the two parameters should be determined by fitting experimental data. If no experimental result is available, α = β = .
0 3 is recommended. According to Huang's work 29 , the two bridging parameters keeps unchanged throughout the loading process.
Step three, calculate the  Step seven, obtain the average value of b, and slightly adjust its value to minimize the error between the predicted and experimental results. Finally, a relation between f G c and G D c is built to describe the mechanical behavior of an interphase.
Please note that, according to Huang's work 29 , α and β keep unchanged throughout a load case. When experimental data is available, the value of α and β should be adjusted by fitting the elastic properties of a composite instead of setting α = β = 0.3. If α and β are not well-determined, the error resulting from elastic range will be transferred to nonlinear part, making the curve fitted interphase properties may deviate from reasonable values to some extent. Besides, the model shown in Eqs. (42,43) is built based on a basic assumption that the shear damage modulus G D c n , and the damage driving stress f G c is one-to-one mapping. In other words, the present model is only applicable for a situation that the damage evolution is monotonous while not applicable for load-unload process. failure criteria. Failure modes. For a UD composite in planar stress state, there are five failure modes, i.e., longitudinal tension, longitudinal compression, transverse tension, transverse compression, and in-plane shear 44 . Similarly, a real stress state can be seen a combination of five stress components, i.e., longitudinal tensile and compressive normal stresses, transverse tensile and compressive normal stresses, and in-plane shear stress. A UD composite is failure when any one of the five failure modes occurs. Matrix failure. Huang 29 presented a generalized criterion for matrix tensile failure, which is given as Eq. (45).   I  m  II  m  ut  m  III  m   I  m  II  m  III  m  ut  m  III  m   3  3 1/3   3  3  3 1/3 where σ σ σ , ,   where σ uc m is matrix compressive strength.

Curves for model validation
Interphase failure. The interphase failure is subordinate for longitudinal tension/compression and transverse compression. However, it is critical for a UD composite under an off-axis load which is the concern of this work. As shown in Eq. (42), a Mohr-Coulomb like damage driving stress of an interphase is developed to consider the effect of transverse stress on shear damage evolution. Generally, for a UD composite under small off-axial angle load, the interphase failure is shear stress dominant. Thus, it is straightforward to take the damage driving stress as a failure criterion. Specifically, an interphase fails when its damage driving stress exceeds a critical value. However, when the off-axial angle increases, the failure mode transforms from in-plane shear to transverse failure. Thus, a maximum normal stress criterion is employed as a complementary criterion. Overall, the interphase failure criterion is given as Eq. (47).

Composites
Empirical equations (MPa) Table 1 lists the elastic properties of the three composites and their constituent materials. Table 2 gives the elastoplastic properties of the three epoxies. Figure 4 shows the off-axis tensile test results of the three UD composites. interphase properties. Before the interphase properties is involved, the value of the bridging parameters α, β must be determined. In general, the α, β are set as 0.3 for most polymeric fibrous composites. Nevertheless, it is better to calibrate the α, β according to experimental data when available. In this work, they can be adjusted by fitting the predicted elastic properties of a composite to measured ones.
In addition, the volume fraction of an interphase should also be set beforehand. As reported by many researchers 35,36 , the interface damage can be simulated by an interphase stiffness degradation. In other words, the interphase is an approximation of an interface but not means a real interphase material. The interface thickness is zero while the interphase thickness is non-zero. If a linear-spring interface model is employed, there exists relations among the interface properties, interphase properties, and the interphase thickness, as shown in the equation below.
n c c c t c c k n and k t are normal and tangent elastic constants of an interface model, and t c , λ c , and G c are thickness and Lame's constants of an interphase model. From the equation, it can be found that, the interphase properties vary with the interphase thickness for a target interface model. However, the values of the interphase properties won't affect the interface model, as long as the relation between the interphase thickness and interphase properties remain unchanged. Thus, for a composite, one can define the interphase thickness arbitrarily in a reasonable range, and then calculate the interphase properties from the interface parameters. According to Mikata et al. 49 and Benveniste et al. 50 , a general interphase volume fraction ranges from 1% to 10%. Further, it is easier to reversely calculate interphase properties with a larger interphase volume fraction. Thus, for the three composites involved in this work, the authors reasonably set all the interphase volume fractions as 10%, although their interface parameters may be different.
Since there are two independent variables in the interphase mode, the shear modulus G D c and the coupling parameter b, two off-axis stress-strain curves are necessary to determine the f G G c D c − relation of an interphase. In order to obtain G D c n , , a shear dominant small off-axis angle test is needed. To obtain b, a 45° off-axis stress-strain curve is advisable.
Following the reverse calculation process shown in Fig. 3, parameters involved in the interphase model can be obtained. Table 3 lists the pre-determined parameters and experimental curves required in the reverse calculation and model validation process.  Table 4 lists the empirical interphase models of the three composites obtained by the reverse calculation process. Experimental validation and discussion. As shown in Table 3, two off-axis curves of a composite are used for reverse calculation of interphase properties. Other off-axial curves are used as validation. Figures 6-8 illustrate the comparison between theoretical and experimental stress-strain curves. Tables 5-7 are strength prediction results of three composites.     www.nature.com/scientificreports www.nature.com/scientificreports/ Overall speaking, it is found from Figs. 6-8 that the predicted results agree well with experimental data. The predicted nonlinear behavior of all the cases shown in Figs. 6-8 are coincident to experimental data, which proves the model's capability on evaluating nonlinear deformation of UD composite under off-axial tension. When focusing on the off-axial tensile strength, most of the strength prediction errors are less than 10%. Only the strength error of AS4/PEEK UD composite under 60° is slightly larger than 10%, which is also good enough. Thus, it is evident that the model performs good in strength prediction. When it comes to the ultimate strain, the prediction errors become larger, but most cases are still satisfied (less than 15%).
For the case of 30° off-axial tension, the prediction result of IM7/8552 UD composite is in high accuracy compared with experimental data. Only the tangential modulus slightly deviates from the experimental value around the failure point, which is not critical for engineering application. However, there are some deviations between the predicted and experimental curves of AS4/PEEK and T300/7901 UD composites under 30° off-axial tension. Such deviations may result from the angle error induced by manual lay-up process. Thus, an off-axial angle of 33° instead of 30° is used for the prediction of the two cases. Then, good correlations are achieved for the nonlinear deformation of the two composites. Nevertheless, the ultimate error in Fig. 8(b) is still not ignorable. It may because the stress based failure criteria is employed. When a material approaches failure, a slight variation of stress may introduce to a significant strain fluctuation. Overall, it should be noted that the ultimate strength and strain are acceptable whether 30° or 33° is used, considering potential source of experimental error.
For the cases of Fig. 6(d), the strength prediction accuracy is acceptable. However, the error of the ultimate strain is relative large. The error may result from three aspects. Firstly, the stress-based failure criterion may be not able to capture the strain variation when a material is approach failure. Secondly, the developed interphase model assumes that the transverse Young's modulus of an interphase, E c , keeps elastic throughout a load case. However, in fact, E c may also degrade with the damage evolution of an interphase. The degradation of E c of AS4/PEEK may be more significant than the other two composites, since the PEEK epoxy is more ductile than the other two epoxies. Lastly, when the off-axial angle becomes large, the failure mode gradually transferred from in-plane failure dominant to tension failure. For T300/7901 UD composite, the coupling parameter b is 0.85, indicating that the effect of transverse stress on failure can be fully considered. However, the engineering reversed coupling parameter b for AS4/PEEK is 0.2, meaning that the effect of transverse stress may be underestimated. Nevertheless, the prediction error is still reasonable (13% for strength and 22 for ultimate strain), considering the fact that the experimental scatter of a composite is generally much larger than a homogeneous material.
The model developed in this work incorporates the contribution of interface damage evolution and matrix elastoplasticity on the nonlinear behavior of a composite. A Mohr-Coulomb like failure criterion is proposed for the stiffness degradation of an interphase, with which the effect of transverse normal stress on the shear property can be considered. Theoretically, the current model is applicable for a UD composite under off-axial tension and compression. However, due to the lack of experimental data, the model is only assessed for UD composites under off-axial tension. For the nonlinear behavior and ultimate strength, the current model provides satisfied prediction results. For the ultimate strain, the prediction error becomes larger but still in a reasonable range. Thus, overall speaking, it is believed that the model is applicable for UD composites under off-axial tension. It also has potential applicability on off-axial compression cases, once validated by experiments. In addition, further investigation on the degradation of Young's modulus and other failure criteria will be made in future works.
conclusions An analytical micromechanics model with matrix elastoplasticity and interface damage is developed for the simulation of nonlinear behavior of a composite. In the model, an interphase with stiffness degradation is employed to simulate the progress interface damage in a composite. For the convenience of engineering application, an empirical model is present to describe the relationship between the shear damage driving stress and the damage property of an interphase. The model is validated by a series of off axial tension tests of three kinds of UD composites. In addition, owing to the explicit and analytical form, the model is of high computational efficiency and convenient to be applied in engineering.  Table 7. Strength prediction results of T300/7901 UD composite.