Analytical Micromechanics Models for Elastoplastic Behavior of Long Fibrous Composites: A Critical Review and Comparative Study

Elasto-plastic models for composites can be classified into three categories in terms of a length scale, i.e., macro scale, meso scale, and micro scale (micromechanics) models. In general, a so-called multi-scale model is a combination of those at various length scales with a micromechanics one as the foundation. In this paper, a critical review is made for the elastoplastic models at the micro scale, and a comparative study is carried out on most popular analytical micromechanics models for the elastoplastic behavior of long fibrous composites subjected to a static load, meaning that creep and dynamic response are not concerned. Each model has been developed essentially following three steps, i.e., an elastic homogenization, a rule to define the yielding of a constituent phase, and a linearization for the elastoplastic response. The comparison is made for all of the three aspects. Effects of other issues, such as the stress field fluctuation induced by a high contrast heterogeneity, the stress concentration factors in the matrix, and the different approaches to a plastic Eshelby tensor, are addressed as well. Correlation of the predictions by different models with available experimental data is shown.


Introduction
Fibrous composites have been vastly used in engineering due to their high specific moduli and strengths and admirable tailorability. Fibrous composites include long and short fiber reinforcements. Short fibrous composites are superior to the long ones in mass-production benefiting from their excellent formability. But the long fibrous composites have significant advantages over short ones on mechanical properties, based on the continuous long shape and improved alignment. In this work, the long fibrous composites are focused on. Unless specified, a composite in this research work is long fibrous.
The elastoplastic properties of a composite are important for the failure and strength evaluation [1,2], damage evolution [3], and dynamic damping analysis [4]. Mechanics of composites in an elastic range has been well-developed [5,6]. Based on an elastic homogenization and a linearization, an elastoplastic model for a composite can be established. However, much effort is still needed to improve the prediction accuracy and efficiency in a plastic range [7][8][9]. In general, elastoplastic theories for a composite can be classified into three categories by a length scale, i.e., macro scale, meso scale, and micro scale, i.e., micromechanics models.
As shown in Figure 1a, a macro model treats a laminate as a homogeneous anisotropic material so that the mechanical response of a composite structure such as a wind turbine blade can be analyzed. Generally, a macro model is incorporated into a finite element approach [10][11][12][13][14][15], in which elastoplastic properties of a laminate element can be obtained from experiments or a model at a smaller scale. For example, Dano et al. [16] developed a two dimensional numerical model for failure analysis of a fastened joint in a composite laminate, where the effective shear elastoplastic behavior of a laminate plate was described by a nonlinear internal variable obtained experimentally. Cooper and Warrior [17] conducted a finite element crash analysis for a composite structure, where an elastoplastic material element was implemented to describe the nonlinear behavior of the structure.
Materials 2018, 11, x FOR PEER REVIEW 2 of 53 smaller scale. For example, Dano et al. [16] developed a two dimensional numerical model for failure analysis of a fastened joint in a composite laminate, where the effective shear elastoplastic behavior of a laminate plate was described by a nonlinear internal variable obtained experimentally. Cooper and Warrior [17] conducted a finite element crash analysis for a composite structure, where an elastoplastic material element was implemented to describe the nonlinear behavior of the structure. It should be noted that an anisotropic yield rule is necessary to describe the elastoplastic behavior, since a laminate is generally highly anisotropic. For example, Schmidt and Weichert [18] proposed an elastoplastic constitutive model for anisotropic shells using the Hill yield criterion [19]. Brünig [20] developed a numerical algorithm for anisotropic plates with the Tsai-Wu criterion to identify a yielding [21]. Under an assumption of an elastic-perfect plastic behavior, the Tsai-Hill yield criterion [22] was employed by Aykul et al. [23] in an elastoplastic analysis of a steel fiber reinforced aluminum matrix composite. With a macro model, analysis of complicated composite structures can be easily carried out by a numerical method. However, since composite materials are highly anisotropic and a plastic deformation is history dependent, the determination of the critical parameters in an anisotropic yield criterion is time and financial consuming [24,25]. Furthermore, how to establish a proper yield theory for a laminate with arbitrary lay-ups is still an open issue [26][27][28][29].
Meso scale models, also referred as layer-wise models, estimate mechanical properties of a laminate from the information of single layers [30][31][32][33][34][35]. In a meso scale approach (Figure 1b), a single layer is treated as a homogeneous and orthogonally anisotropic media. Its mechanical properties can be directly measured through experiments or obtained from a micromechanics model. One major problem for an elastoplastic model at a meso scale is how to describe the nonlinearity involved. Rotem [36,37] expressed the shear component in the stiffness tensor of a layer as a nonlinear function of a shear strain dependent parameter, while all the other components remained elastic. Pinho et al. [38] and Wolfe et al. [39] utilized a secant and a tangent linearization, respectively, to approximate the elastoplastic behavior of an individual layer. An instantaneous stiffness tensor was obtained by replacing the elastic modulus in it with a secant or tangent one.
To better understand the nonlinearity of a composite, it is necessary to consider an interaction between normal and shear stresses. Puck et al. [40] introduced a concept of stress exposure ratio to account for the interaction between the transverse normal and the in-plane shear stresses. Kress [41] expressed the transverse normal and the in-plane shear stresses in terms of the Hashin's second and fourth invariants [42], respectively. Moreover, to fully address the stress interaction, anisotropic yield criteria for a layer are widely developed and employed [19,21,22,[43][44][45][46][47][48][49][50]. For example, Sun et al. [51] proposed a one parameter plasticity model to describe the elastoplastic behavior of a UD (unidirectional) composite. The Tsai-Wu yield criterion was employed by Pisano et al. [45] to analyze the failure behavior of a pinned-joint composite laminate based on a layer-wise approach. However, It should be noted that an anisotropic yield rule is necessary to describe the elastoplastic behavior, since a laminate is generally highly anisotropic. For example, Schmidt and Weichert [18] proposed an elastoplastic constitutive model for anisotropic shells using the Hill yield criterion [19]. Brünig [20] developed a numerical algorithm for anisotropic plates with the Tsai-Wu criterion to identify a yielding [21]. Under an assumption of an elastic-perfect plastic behavior, the Tsai-Hill yield criterion [22] was employed by Aykul et al. [23] in an elastoplastic analysis of a steel fiber reinforced aluminum matrix composite. With a macro model, analysis of complicated composite structures can be easily carried out by a numerical method. However, since composite materials are highly anisotropic and a plastic deformation is history dependent, the determination of the critical parameters in an anisotropic yield criterion is time and financial consuming [24,25]. Furthermore, how to establish a proper yield theory for a laminate with arbitrary lay-ups is still an open issue [26][27][28][29].
Meso scale models, also referred as layer-wise models, estimate mechanical properties of a laminate from the information of single layers [30][31][32][33][34][35]. In a meso scale approach (Figure 1b), a single layer is treated as a homogeneous and orthogonally anisotropic media. Its mechanical properties can be directly measured through experiments or obtained from a micromechanics model. One major problem for an elastoplastic model at a meso scale is how to describe the nonlinearity involved. Rotem [36,37] expressed the shear component in the stiffness tensor of a layer as a nonlinear function of a shear strain dependent parameter, while all the other components remained elastic. Pinho et al. [38] and Wolfe et al. [39] utilized a secant and a tangent linearization, respectively, to approximate the elastoplastic behavior of an individual layer. An instantaneous stiffness tensor was obtained by replacing the elastic modulus in it with a secant or tangent one.
To better understand the nonlinearity of a composite, it is necessary to consider an interaction between normal and shear stresses. Puck et al. [40] introduced a concept of stress exposure ratio to account for the interaction between the transverse normal and the in-plane shear stresses. Kress [41] expressed the transverse normal and the in-plane shear stresses in terms of the Hashin's second and fourth invariants [42], respectively. Moreover, to fully address the stress interaction, anisotropic yield criteria for a layer are widely developed and employed [19,21,22,[43][44][45][46][47][48][49][50]. For example, Sun et al. [51] proposed a one parameter plasticity model to describe the elastoplastic behavior of a UD (unidirectional) composite. The Tsai-Wu yield criterion was employed by Pisano et al. [45] to analyze resulted in lower computational efficiency [120]. Heinrich et al. [121] claimed that an RVE containing at least 25 fibers could provide a satisfactory prediction accuracy. Determination of an RVE size has also been discussed by some other researchers [122][123][124][125].
RUC based models are applicable to a composite with periodic microstructures. Only one or several fibers are included in an RUC. The boundary condition, a critical point of an RUC model, has to be defined carefully to represent effect of different fiber distribution patterns and loads applied [54,58,89,[126][127][128][129]. In some cases, the uniform boundary conditions are applied to RUCs. For example, Brockenbrough, et al. [130] investigated the effect of fiber distribution patterns on the response of a metal matrix composite with a uniform strain boundary condition. Also, similar treatment can be found in the work of Aboudi [131], Aghdam et al. [132], and Würkner et al. [133]. However, it was reported by Suquet [89] that an RUC model with such a uniform boundary condition could only give an upper or lower bound of the effective properties of a composite. Specifically, the effective stiffness of a composite would be overestimated by an RUC model with a uniform strain while underestimated with a uniform stress assumption. Suquet [89] gave a rigorous definition of the periodic boundary condition with which a better estimation of the effective properties was achieved. Xia, et al. [134,135] proposed a unified periodic boundary condition for an RUC under any combined loads. They indicated that the uniform boundary condition not only over constrained the RUC but might violate the boundary traction periodicity as well. A periodic boundary condition has been widely applied with both RVE [123,[136][137][138] and RUC [139][140][141][142][143]  Elastic RVE and RUC models can be extended to nonlinear cases, providing that the nonlinear constitutive laws for the constituents are available. For example, based on a three-dimensional RVE model, Yuan and Lu [144] conducted a numerical investigation on the elastoplastic behavior of RUC based models are applicable to a composite with periodic microstructures. Only one or several fibers are included in an RUC. The boundary condition, a critical point of an RUC model, has to be defined carefully to represent effect of different fiber distribution patterns and loads applied [54,58,89,[126][127][128][129]. In some cases, the uniform boundary conditions are applied to RUCs. For example, Brockenbrough, et al. [130] investigated the effect of fiber distribution patterns on the response of a metal matrix composite with a uniform strain boundary condition. Also, similar treatment can be found in the work of Aboudi [131], Aghdam et al. [132], and Würkner et al. [133]. However, it was reported by Suquet [89] that an RUC model with such a uniform boundary condition could only give an upper or lower bound of the effective properties of a composite. Specifically, the effective stiffness of a composite would be overestimated by an RUC model with a uniform strain while underestimated with a uniform stress assumption. Suquet [89] gave a rigorous definition of the periodic boundary condition with which a better estimation of the effective properties was achieved. Xia, et al. [134,135] proposed a unified periodic boundary condition for an RUC under any combined loads. They indicated that the uniform boundary condition not only over constrained the RUC but might violate the boundary traction periodicity as well. A periodic boundary condition has been widely applied with both RVE [123,[136][137][138] and RUC [139][140][141][142][143] models.
Elastic RVE and RUC models can be extended to nonlinear cases, providing that the nonlinear constitutive laws for the constituents are available. For example, based on a three-dimensional RVE model, Yuan and Lu [144] conducted a numerical investigation on the elastoplastic behavior of carbon nanotubes (CNTs) reinforced polymer composites. Hoang et al. [124] studied the effect of an RVE size on the prediction for the elastoplastic and elasto-viscoplastic behavior of a two-phase composite. Choosing an RUC, Aghdam et al. [132] analyzed the yield and collapse behavior of a metal matrix composite. With an RUC based multiscale model, Wan et al. [145] studied the compressive behavior of a braided composite after an impact accounting for the elastoplastic deformation of the matrix.

Analytical Micromechanics Models
Analytical micromechanics models have been well developed for elastic problems but are still in progress for elastoplastic ones. Linearization is a practical way to establish an elastoplastic micromechanics model, in which a nonlinear response is discretized and approximated as a series of linear problems. With such a linearization, a micromechanics model can be extended to the analysis of an elastoplastic problem.
The selection of a proper micromechanics model is fundamental. The equivalent inclusion method pioneered by Eshelby [167] can be used to solve an eigenstrain problem with a single inclusion embedded in an unbounded matrix. As long as a fibrous composite is concerned, the interaction between the inclusion and the surrounding fibers is ignored in the Eshelby's method. A self-consistent model (SCM) [55] is advanced from the Eshelby's method by replacing the matrix phase with a medium whose properties equal to the effective properties of the composite. Although better in a prediction accuracy, the implicit formulae resulted from the SCM make it inconvenient for use in engineering. Moreover, the SCM may give a physically nonsense result when the inclusion is rigid or a void [168]. Afterwards, Mori and Tanaka developed a method [56] by keeping the same configuration of the Eshelby model but assigning the homogenized strain field of the matrix as that of the composite. Further, a generalized self-consistent method (GSCM) was proposed [169][170][171], dealing with a configuration that a single fiber surrounded by the matrix is embedded in an infinite homogeneous medium whose properties are the same as those of the composite. In addition, two or three phase concentric cylinder assembly (CCA) models are also widely used to determine the effective properties of the composite [170,[172][173][174][175][176][177][178]. Note that the two-phase and three-phase CCA models are equivalent to the Mori-Tanaka and the GSCM, respectively, when the inclusions are cylindrical [168,179]. Moreover, a double inclusion approach [180][181][182] is proposed to evaluate the mechanical properties of a multiphase composite, of which the Mori-Tanaka model, the SCM, and the GSCM can be considered as a special case.
In addition to the theory of elasticity based models mentioned above, some semi-empirical analytical models are also well-known. With the iso-stress and iso-strain assumptions, the classical rule of mixture [183,184] is resulted. The rule of mixture can give good prediction for the longitudinal modulus and the major Poisson's ratio but deliver poor estimation for the other elastic properties of a composite. Based on experimental results, Chamis et al. [185,186] achieved a model called Chamis' model, which can be considered as a modified version of the rule of mixture. Much better prediction in the transverse and shear moduli is generally seen. Halpin and Tsai [187] presented a simplified version of the SCM, known as Halpin-Tsai's equations.
Huang [188] developed a unified elasto-plastic bridging model in which a bridging tensor is employed to link the homogenized stresses of the matrix with those of the fiber. The bridging tensor elements can be classified as dependent and independent. Whereas the dependent elements are determined from the symmetric condition of the overall compliance tensor of the composite, the independent ones are expressed in the form of Taylor series expansions with respect to the fiber and matrix property parameters. The expansion coefficients, independent of any constituent property, were determined based on existing micromechanics theories in an elastic range. It can be seen that the independent bridging tensor elements thus obtained are simplifications and modifications to the corresponding counterparts of the Mori-Tanaka model [179,189]. As long as a plastic deformation of a constituent e.g., matrix material is concerned, only the corresponding matrix property parameters in the Taylor expansion need to adjust. It has been shown that the bridging model can give better correlation between the predicted and measured effective elastic properties of UD composites than many other famous models [190][191][192][193].
The second step in the establishment of an elastoplastic model is the selection of a linearization. Hill [194] proposed an incremental linearization in which stress and strain increments were linked with an instantaneous stiffness or compliance tensor. Consider a homogeneous material subjected to uniaxial tension. An instantaneous Young's modulus is defined as a tangent to the stress-strain curve of the material at the current stress analysis point. Hence, an incremental linearization scheme is also referred to as a tangent approach. The incremental linearization has been widely used due to its capability on history dependent cases [9,[195][196][197]. Following Hill's work [194], Berveiller and Zaoui [198] and Tandon and Weng [199] proposed a deformation linearization in which the total stress and strain were connected by a secant stiffness or compliance tensor. Under a uniaxial tension, a secant Young's modulus is defined as the secant slope between an objective point and the initial point on the stress-strain curve. Puck and Schurmann [200] pointed out that the secant method offers the advantage of self-correction, meaning that the error induced by the secant method in the current load step will not be transferred to the next step whereas a tangent approach will be done due to its incremental nature. The deformation linearization, referred to as a secant approach, also gains a popularity in application benefiting from its simple formulation and good prediction accuracy [201][202][203]. In comparison, it is seen that the tangent model applied to non-proportional and non-monotonic loads whereas a secant model is generally restricted to proportional and monotonic cases [114,199,202,[204][205][206][207]. In addition to these two linearizations, Dvorak [208] proposed a transformation field analysis (TFA) method. In the TFA, properties of the constituents keep elastic. The plastic strain in a constituent is taken as an eigenstrain whose effect on a composite response is reflected by an influencing function. The TFA has been implemented into an FEA for nonlinear analysis of composites [209].
It has been reported that the nonlinear stress strain curve evaluated by a tangent approach, a secant approach, or the TFA may be too stiff compared with experiments [95,[210][211][212][213]. Existing attempts to address the too stiff problem are mainly in three kinds, i.e., improvements on a linearization, corrections on an equivalent stress, and modifications on the calculation of a plastic Eshelby's tensor.
Regarding the first kind attempt, Molinari [214] suggested that the too stiff response might be amended by a full consideration of an interaction between plastic deformation of the fiber and matrix. He then proposed a non-incremental tangent linearization approach, which was validated by Molinari et al. [215] and Mecier et al. [216]. Using the same linearization [217], Masson et al. [218] and Zaoui et al. [219] proposed an affine formulation for an elastoplastic response of a composite. The affine formulation was further developed and validated by Pierard and Doghri [110], Pierard et al. [220], and Doghri et al. [221]. On the other hand, Wu et al. [112][113][114] believed that better results could be obtained by an incremental-secant scheme they proposed. In their approach, when the load was increased from the current to the next step, the effective properties of a composite were approximately estimated by the secant method. Afterwards, a fictitious elastic unload was introduced, from which residual stresses and strains induced by the plastic deformation were obtained and accumulated [112][113][114]. One feature of this incremental-secant scheme is in its applicability to non-proportional and non-monotonic load problems.
As for the second kind attempt, it was pointed out that the too stiff response might be resulted from an ignorance of a stress field fluctuation in the constituents [205]. Instead of calculating the von-Mises equivalent stress from homogenized stresses of the constituent materials in a composite, Qiu and Weng [205] and Hu [222] derived the von-Mises equivalent stress from an energy approach which can reflect the effect of stress fluctuation on yield behavior. Inspired by their work, Suquet [223] proposed a modified secant model, where the second-moment estimation for stresses was used to capture the stress fluctuation in the constituent phases of a composite. Independently, basing on the work of Talbot and Willis [224], Ponte Castaneda [225,226] proposed a rigorous variational principle to give bounds or estimates for nonlinear composites. It is interesting to find that the modified secant model by Suquet [223] coincides with the variational estimation by Ponte Castaneda [225][226][227] when the complementary energy of the constituent phases is a quadratic form of the stress tensor [228].
Regarding the third attempt, there are four approaches in the literature to determine a plastic Eshelby's tensor. The first one [229][230][231][232] is named as anisotropic Eshelby tensor method, meaning that an instantaneous stiffness tensor of the matrix under an elastoplastic deformation is computed rigorously from a plasticity theory, such as Prandtl-Reuss model, and the corresponding Eshelby tensor is also anisotropic and evaluated through a numerical integration. This approach gives the most rigorous tensor but may result in a too stiff response [65,197,233,234]. The second is referred as isotropicalized matrix method. In this method, the instantaneous stiffness of a matrix can be approximated as isotropic, and the corresponding Eshelby tensor is isotropic as well [65,145,160]. Namely, the Eshelby tensor in a plastic region is similar to that in an elastic one. However, an isotropization has been demonstrated successful only for a proportional load case [212]. The third approach is designated as isotropicalized Eshelby tensor method. In this approach, an instantaneous stiffness tensor of the matrix remains anisotropic whereas the Eshelby tensor is determined according to an isotropic condition. Namely, only the elastic property of the matrix in the Eshelby tensor of an elastic region is replaced with a plastic counterpart of the matrix. This method was first proposed by Doghri and Ouaar [65], and was further applied by many other researchers [111,197,[234][235][236][237][238][239]. However, Huang et al. [240] and Peng et al. [116] pointed out that the isotropization of the Eshelby tensor has no sound physical background. The fourth approach has been recently presented by Peng et al. [116]. At each load step, a reference elastic medium (REM) is introduced whose configuration and properties are identical to the instantaneous ones of the elastoplastic medium (EPM). Based on the REM, Peng et al. obtained a scheme to determine the Eshelby tensor with modified prediction on the elastoplastic behavior of a composite.
Very recently, Huang et al. [103][104][105][107][108][109]241] have shown that the homogenized stresses of the matrix should be converted into true values before an effective property of the composite can be determined in terms of the monolithic fiber and matrix properties. The conversion is achieved by multiplying the homogenized stresses with the respective stress concentration factors (SCFs) of the matrix due to introduction of the fiber, and all of the SCFs in relation to different loads under perfect as well as debonded fiber-matrix interface have been derived [103][104][105][107][108][109]241]. The elastic properties of the composite are independent of the stress values in the constituents, and thus the true stress concept plays no role. However, when the matrix is in a plastic deformation, the true stresses may result in an instantaneous stiffness of the matrix lowered down significantly in comparison with the homogenized counterparts. This is because essentially all of the SCFs are greater, and some are significantly higher, than 1 especially after the interface debonding. In turn, the composite stiffer response resulted from the homogenized stresses in the fiber and matrix can be satisfactorily addressed.
To better understanding the establishment of an elastoplastic model, Figure 3 summarized the process of extending an elastic model to an elastoplastic range.

Micromechanics Models with Imperfect Interface
The afore-mentioned models/theories are mainly built for a composite with a perfect interface. A perfect interface condition means that the stress and displacement fields are continuous at the interface in-between the fiber and matrix. Contrarily, imperfect interface condition implies that the stress or displacement field is discontinuous at the interface. Perfect interface condition is applicable for most engineering applications of composite materials. However, imperfect interface situations do exist in some cases. For example, interface debonding occurs when a composite subjected to a fatigue load that exceeds its elastic limit, e.g., Figure 4a [242]. Cracks or micro-voids are often observed at the interface in a thermo-pressure lamination process of a kind of metal matrix fibrous composite, e.g., Figure 4b [243]. In addition, an interphase is produced due to the chemical reaction between the fiber and matrix (see Figure 4c) [244]. Moreover, coating, also a kind of interphase, is often added to a fiber for the purpose of function design, as shown in Figure 4d [245]. Both the interface debonding/crack and the interphase can be classed to imperfect interface condition. It is reported that the effect of imperfect interface on mechanical properties of a composite is unignorable in some cases [246][247][248].

Micromechanics Models with Imperfect Interface
The afore-mentioned models/theories are mainly built for a composite with a perfect interface. A perfect interface condition means that the stress and displacement fields are continuous at the interface in-between the fiber and matrix. Contrarily, imperfect interface condition implies that the stress or displacement field is discontinuous at the interface. Perfect interface condition is applicable for most engineering applications of composite materials. However, imperfect interface situations do exist in some cases. For example, interface debonding occurs when a composite subjected to a fatigue load that exceeds its elastic limit, e.g., Figure 4a [242]. Cracks or micro-voids are often observed at the interface in a thermo-pressure lamination process of a kind of metal matrix fibrous composite, e.g., Figure 4b [243]. In addition, an interphase is produced due to the chemical reaction between the fiber and matrix (see Figure 4c) [244]. Moreover, coating, also a kind of interphase, is often added to a fiber for the purpose of function design, as shown in Figure 4d [245]. Both the interface debonding/crack and the interphase can be classed to imperfect interface condition. It is reported that the effect of imperfect interface on mechanical properties of a composite is unignorable in some cases [246][247][248].

Micromechanics Models with Imperfect Interface
The afore-mentioned models/theories are mainly built for a composite with a perfect interface. A perfect interface condition means that the stress and displacement fields are continuous at the interface in-between the fiber and matrix. Contrarily, imperfect interface condition implies that the stress or displacement field is discontinuous at the interface. Perfect interface condition is applicable for most engineering applications of composite materials. However, imperfect interface situations do exist in some cases. For example, interface debonding occurs when a composite subjected to a fatigue load that exceeds its elastic limit, e.g., Figure 4a [242]. Cracks or micro-voids are often observed at the interface in a thermo-pressure lamination process of a kind of metal matrix fibrous composite, e.g., Figure 4b [243]. In addition, an interphase is produced due to the chemical reaction between the fiber and matrix (see Figure 4c) [244]. Moreover, coating, also a kind of interphase, is often added to a fiber for the purpose of function design, as shown in Figure 4d [245]. Both the interface debonding/crack and the interphase can be classed to imperfect interface condition. It is reported that the effect of imperfect interface on mechanical properties of a composite is unignorable in some cases [246][247][248].   (c) Interphase produced by chemical reaction in a SiC/Ti-6Al-4V composite [244]; (d) BN coated T300 fiber [245].
In elastic range, micromechanics models with an imperfect interface can be classified into two categories [249], i.e., interface model and interphase model. In an interface model, a crack-like zero-thickness interface is employed to characterize the imperfect interface between fiber and matrix. The stress or displacement fields across the interface are discontinuous. Obviously, an interface model is proper for the cases shown in Figure 4a,b. The linear-spring model [250,251] and interface stress model [252,253] are the two well-known interface models. At the interface, the former assumes that the stress field is continuous while the displacement field is not. The displacement variance is proportional to the stress. The latter assumes a continuity of displacement but discontinuity of stress. This model is generally used to describe interface compression phenomenon. Other models, like dislocation-like model [254], interface sliding model [255,256], and the anti-interpenetration model [257] can be seen as extensions of the linear-spring and interface stress models. An interphase model adds a layer between the fiber and matrix in a two-phase concentric cylindrical assembly (CCA) [174]. The mechanical property of the layer is different from those of the fiber and matrix. Clearly, an interphase model is appropriate for the cases shown in Figure 4c,d. Hashin [174] gave an exact solution for a three-phase CCA model with thin interphase. Further, Benveniste [178] extended Hashin's results to a three-phase CCA model with thick interphase. It is pointed out that an interphase model is equivalent to an interface model when the interphase is thin and soft [249,258,259]. The exact solution of an interphase model is complex. Based on a two-phase bridging model and an equivalent fiber method, Wang [179,189,260] proposed a simplified analytical three-phase model for the analysis of composites with imperfect interface.
In nonlinear range, for a composite with an imperfect interface, it is reasonable to assume that the nonlinearity of a composite mainly comes from the elasto-plastic behavior of the matrix and the progress damage of the interface. Chang et al. [261] developed a progressive damage model for a composite laminate. But their model is in macro-scale, meaning that the stress analysis at micro scale is not available. In addition, the contribution of matrix and interface nonlinearity cannot be distinguished. Ju et al. [262] proposed a micromechanics interface damage model. In their work, the interface damage was approximately described by an inhomogeneity with transverse isotropy. However, their work is only valid for particle reinforced composites. The cohesive element method [3,263] has been widely used in simulation of a composite with an imperfect interface. The cohesive element method has sound physical background and is powerful in the analysis of interface crack propagation. However, the cohesive model has to be implemented into a micro-scale FEM, thereby being of high cost in computation resource. For engineering applications, it is desirable to develop a model with abilities in micro-scale damage analysis, satisfied prediction accuracy, and high computational efficiency.
It is noted that models with imperfect interface have been vastly investigated. But, for most engineering applications, the perfect interface assumption is good enough for mechanical analysis. In this work, a composite with perfect interface is mainly focused on.

Comparison on Elastic Theories
The selection of a proper elastic micromechanics theory is in the fundamental step to establish an elastoplastic model. In this section, summary of different elastic micromechanics models is made. Then, the elastic models are evaluated regarding their capabilities in predicting elastic properties of UD composites. Further, based on a tangent linearization, the elastic models are extended to be valid in an elastoplastic range. A quantitative comparison for them is shown in the next section.

General Framework
Consider an RUC shown in Figure 5. The stress and strain components must be homogenized with respect to the volume of the RUC through Equations (1) and (2).
where σ(x) and ε(x) are, respectively, point-wise stress and strain tensors, σ and ε the homogenized counterparts. Since only the homogenized quantities are dealt with, the over bars are omitted.
where ( ) and ( ) are, respectively, point-wise stress and strain tensors, and the homogenized counterparts. Since only the homogenized quantities are dealt with, the over bars are omitted. For a two-phase composite with fiber and matrix constituents, the stress and strain of a composite are given by Equations (3) and (4).
f and m designate the fiber and matrix, respectively. A quantity with no suffix belongs to a composite. Following Hill [264], there are two fourth-order strain and stress concentration tensors, and , as shown in Equations (5) and (6).
Let and denote compliance and stiffness tensors with Equations (7) and (8).
Equation (11) is also useful: The strain (stress) concentration tensor, ( ), that connects a strain (stress) tensor of a constituent with that of a composite is named as a global strain (stress) concentration tensor. Very often, it is easier to establish a model with a local concentration tensor. Let and represent the local strain and stress concentration tensor, respectively, such as Equations (12) and (13).

= :
= : The connections between the local and the global concentration tensors are found as shown in Equations (14) and (15). For a two-phase composite with fiber and matrix constituents, the stress and strain of a composite are given by Equations (3) and (4).
f and m designate the fiber and matrix, respectively. A quantity with no suffix belongs to a composite. Following Hill [264], there are two fourth-order strain and stress concentration tensors, A r and B r , as shown in Equations (5) and (6).
Let M and L denote compliance and stiffness tensors with Equations (7) and (8).
Then, the effective stiffness and compliance tensor are given by Equations (9) and (10).
Equation (11) is also useful: The strain (stress) concentration tensor, A r (B r ), that connects a strain (stress) tensor of a constituent with that of a composite is named as a global strain (stress) concentration tensor. Very often, it is easier to establish a model with a local concentration tensor. Let T f and P f represent the local strain and stress concentration tensor, respectively, such as Equations (12) and (13).
The connections between the local and the global concentration tensors are found as shown in Equations (14) and (15).
Further, once the stiffness and compliance tensors, L and M, are known, the global concentration tensors are given (see Equations (16) and (17)).
The key to a micromechanics model is the determination of a global or local strain/stress concentration tensor. From Equations (5), (6), (12) and (13), it is found that the determination of concentration tensors requires knowledge of the stress and strain fields in the constituent phases of a composite. However, as shown in Figure 6a, due to the interaction between adjacent fibers, it is arduous to obtain an exact stress/strain field of a constituent in a multi-inclusion model. Therefore, a reference medium is introduced, with which the fiber interaction can be approximated in a single fiber model as shown in Figure 6b. The strain and stiffness tensors of the reference medium are denoted by ε re and L re , respectively. A specific definition of the ε re and L re leads to a specific micromechanics model. = : ( + ) (15) Further, once the stiffness and compliance tensors, and , are known, the global concentration tensors are given (see Equations (16) and (17)).
The key to a micromechanics model is the determination of a global or local strain/stress concentration tensor. From Equations (5), (6), (12) and (13), it is found that the determination of concentration tensors requires knowledge of the stress and strain fields in the constituent phases of a composite. However, as shown in Figure 6a, due to the interaction between adjacent fibers, it is arduous to obtain an exact stress/strain field of a constituent in a multi-inclusion model. Therefore, a reference medium is introduced, with which the fiber interaction can be approximated in a single fiber model as shown in Figure 6b. The strain and stiffness tensors of the reference medium are denoted by and , respectively. A specific definition of the and leads to a specific micromechanics model.

Eshelby Model
It was established on a single fiber model ( Figure 6b). Let it be subjected to a uniform traction . Suppose that [167] = and = where and in Equation (18) stand for the stiffness and strain tensor of the reference medium in Figure 6b, respectively. Since the matrix is infinite, the effect of the fiber on the total strain of the model is neglected, leading to Equation (19).
Making use of the Eshelby equivalent inclusion method, the following Equations (20) and (21)

Eshelby Model
It was established on a single fiber model ( Figure 6b). Let it be subjected to a uniform traction σ 0 . Suppose that [167] L re = L m and ε re = ε m where L re and ε re in Equation (18) stand for the stiffness and strain tensor of the reference medium in Figure 6b, respectively. Since the matrix is infinite, the effect of the fiber on the total strain of the model is neglected, leading to Equation (19).
Making use of the Eshelby equivalent inclusion method, the following Equations (20) and (21) are obtained where ε pt is the perturbed strain tensor due to the presence of the fiber, ε * is an eigenstrain, and S m is an Eshelby tensor. The superscript ES designates the Eshelby method. Comparing Equation (5) with (20), the global and the local strain concentration tensors are the same as shown in Equation (22).
Owing to ignoring the interaction of the inclusion with the surrounding fibers, the Eshelby model is applicable only to a composite with a low fiber volume fraction.

SCM
In a SCM [55], the stiffness and strain of the reference medium equal to those of the composite, shown as Equation (24).
The strain in the fiber is obtained as Equations (25) and (26) where S is the Eshelby tensor from the composite medium. Therefore, Equation (27) is obtained.
Substituting Equations (25)-(27) into (9), the stiffness tensor of the composite is given by Equation (28) is implicit since both sides of it contain the unknown stiffness tensor. Further, the SCM would yield a physical nonsense result for a composite with a rigid or void inhomogeneity [168].

GSCM
A schematic configuration for the GSCM is shown in Figure 7. Hill [264] and Hashin [172] presented exact solutions for four of the five effective elastic constants, E11, , G12, and E22, as Equations (34)-(38) where , , and are the transverse bulk moduli of the composite, fiber, and matrix given, respectively, as Equation (39).
represents a stiffness component of a constituent material.
As for the transverse shear modulus , only the upper and lower bounds were provided by Hill [264] and Hashin [172]. Christensen and Lo [171] derived an explicit solution for when the composite is made of isotropic fiber and matrix. Luo and Weng [265] obtained the displacement fields in the fiber (f), matrix (m), and the composite media of Figure 7. Both the fiber and matrix can be transversely isotropic. Luo and Weng's solutions are shown as Equations (40)-(46). Hill [264] and Hashin [172] presented exact solutions for four of the five effective elastic constants, E 11 , υ 12 , G 12 , and E 22 , as Equations (34)-(38) where k, k f , and k m are the transverse bulk moduli of the composite, fiber, and matrix given, respectively, as Equation (39).
L r ij represents a stiffness component of a constituent material.
As for the transverse shear modulus G 23 , only the upper and lower bounds were provided by Hill [264] and Hashin [172]. Christensen and Lo [171] derived an explicit solution for G 23 when the composite is made of isotropic fiber and matrix. Luo and Weng [265] obtained the displacement fields in the fiber (f ), matrix (m), and the composite media of Figure 7. Both the fiber and matrix can be transversely isotropic. Luo and Weng's solutions are shown as Equations (40)-(46).
where d i , i = 1, 2 . . . 8, are unknown coefficients to be solved using the continuity conditions. Once the homogenized stresses and strains of the three-phases are determined from the displacement fields, the shear modulus is obtained.

Chamis Model
The longitudinal Young's modulus and the major Poisson ratio in Chamis model [271] coincide with the rule of mixture. The other three moduli are given by Equations (52)- (54).

Halpin-Tsai Equations
As stated by Halpin and Kardos [187], the Halpin-Tsai equations were modified from those of GSCM with some engineering based considerations. The expressions for E 11 and υ 12 are consistent with Equations (47) and (48), whereas E 22 is calculated from Equation (37). Table 2 shows the expressions for the other moduli.

12
T , r = f , m are the homogenized stress vectors of the fiber and matrix. The explicit bridging tensor A ij is as Equations (57)-(60).
β and α are the bridging parameters to better correlate the resulting E 22 and G 12 with experiments. If no experiments are available, β = α = 0.3 are mostly recommended. The Equation (61) for a 12 and a 13 are solved from the symmetric condition of the composite compliance, i.e., M ij = M ji .

Quantitative Comparison
A comprehensive quantitative comparison on micromechanics models is desirable both in academic and engineering applications. However, it is almost impossible collect all experiment data in literature for validation. Consider that the WWFE is a worldwide well-known and trustworthy academic activity. Measured elastic properties of the 9 UD composites together with the monolithic fiber and matrix property parameters provided in three world-wide failure exercises (WWFEs) [272][273][274] are used to benchmark the predictions by the different micromechanics models. In addition, three numerical models, i.e., the FEM, FVDAM, and the GMC, are also compared. Please note that all the 9 UD composites are all long fibrous, epoxy matrix, and fiber volume fraction of around 60%. Thus, cautiously speaking, the comparison results are limited to elastic behaviors of long fiber reinforced epoxy matrix with intermediate-high fiber volume fraction.
Since all of the models give essentially the same results for E 11 , only the averaged correlation errors for the other four moduli are shown in Tables 3-6. Table 7 summarizes the overall averaged errors for the five elastic constants. Information of the 9 composites is given in Appendix A Table A1, whereas detailed predictions by all the models are listed in Tables A2-A12.    Table 7 shows that bridging model with α = β = 0.3 gives the best overall prediction accuracy for elastic behaviors of the 9 UD composites among all the homogenization models involved. In addition, the expressions of the bridging model for homogenized stresses of the fiber and matrix are explicit and the simplest, making it convenient in application. Another advantage of the bridging model is in the bridging parameters, α and β, which are semi-empirical to implicitly represent an effect of some uncertain factors such as random fiber arrangement or imperfect interface on a composite response. Considering the adjustability of α, β, it is worthy to expect that the bridging model may also give satisfactory prediction of elastic properties for other kinds of long fiber composites in addition to the 9 UD ones. Bridging model has been programmed into a general-purpose user-subroutine, Bridging model for analysis of composites (BMANC) [275], which can be combined with an FEM software package such as ABAQUS to analyze linear, nonlinear, failure, and strength behaviors of a complex structure with composites involved.
Chamis model is simple and good in prediction accuracy. The FEM and the FVDAM predictions also agree well with the experimental data. Besides, the numerical models can deal with more complicated fiber-induced geometry such as irregular fiber cross-section or random fiber distribution. Once in a while with no experimental data available, the FEM is used to benchmark other kinds of solutions. The GMC is another kind of numerical model but gives inferior accuracy due to the assumption of uniform strains in cells. Better predictions can be attained by a high fidelity generalized method of cell (HFGMC) [276], but its computational efficiency is lower.
SCM, GSCM, Mori-Tanaka, and Halpin-Tsai give comparable predictions. Amongst, GSCM performed the best. On the other hand, Halpin-Tsai equations and Mori-Tanaka model are more widely used due to their simplicity. It should be noted that the prediction accuracy of G 23 by SCM is not satisfactory (61% error), although its results for the other four constants are good. For a composite containing rigid or void inhomogeneity, SCM may lead to non-physical G 23 [168].
The remaining two models, the rule of mixture and Eshelby model, ranked the lowest. The uniform strain/stress assumption used in the rule of mixture violates the continuity condition at fiber/matrix interface, whereas the major error of Eshelby model is resulted from the ignorance of the fiber interaction.

Comparison on Elastoplastic Behavior
In this work, a comparison of elastoplastic theories is restricted to static, monotonic, and proportional loads. A brief introduction is given on whether a model is applicable to non-monotonic and non-proportional load conditions. Three kinds of UD composites are taken for example in the comparative study, i.e., E-glass/Epoxy, IM7/8551-7, and AS4/Peek UD composites. E-glass is a kind of glass fiber and IM7 and AS4 are carbon fiber. Epoxy and 8551-7 are thermoset while Peek is thermoplastic. Thus, the three composites can reasonably represent most common seen fiber reinforced plastic composites. Regarding the load conditions, we choose transverse compression, in-plane shear, and off-axial tension as examples. It is because the longitudinal tensile/compressive and transverse tensile behaviors of a composite are usually linear elastic. Without consideration of interface damage, the nonlinearity of a composite majorly comes from transverse compression and in-plane shear deformation. In addition, off-axial tensile can be seen as a combination of transverse tension and in-plane shear, representing a kind of multi-axial load case. Thus, we choose the three kinds of load cases as benchmark.
A rule is needed to judge the efficiency of different models for the predicted elastoplastic responses. Let us choose three parameters. They are the elastic modulus E or G, the yield point σ Y , and the asymptotic tangent modulus E T asy or G T asy as schematically shown in Figure 8. The latter one, E T asy or G T asy , is defined as the minimum value of the tangents to the tensile or shear stress strain curve of a material [57,95]. A predicted elastic constant affects the prediction accuracy of the stress-strain curve in the elastic part. The evaluation of the yield point determines when a material yields. In fact, the slope variation of the stress strain curve depends on the predicted yield point. In other words, the reduction rate of a predicted modulus is controlled by the yield point. The asymptotic tangent modulus E T asy or G T asy is a new concept introduced in this work. It is defined as the minimum value of the tangents to the tensile or shear stress strain curve of a material. As shown in Equation (62), a tangential modulus is negatively correlated to a plastic strain. Thus, the asymptotic tangent modulus can partly reflect how much plastic strain a material can bear. Equation (62), a tangential modulus is negatively correlated to a plastic strain. Thus, the asymptotic tangent modulus can partly reflect how much plastic strain a material can bear. In addition to the or , , and or , an overall error, as shown in Equation (63), is introduced to characterize the overall prediction performance of a model.
where n is the amount of available experiment curves. As shown in Figure 8, is the enclosed area between the predicted curve and strain-axis, and is the experimental one. For a uniaxial tension/compression or pure shear stress state, the enclosed area corresponds to the accumulated strain energy. For a complex stress state, the area corresponds to a part of the strain energy. Note that can only reflect one valuable side view of a model's capability. When choosing a model, readers should comprehensively evaluate the applicability of a model from their own viewpoint. In addition to the E or G, σ Y , and E T asy or G T asy , an overall error, ER ov as shown in Equation (63), is introduced to characterize the overall prediction performance of a model.

Comparison on Micromechanics Models
where n is the amount of available experiment curves. As shown in Figure 8, R pr is the enclosed area between the predicted curve and strain-axis, and R ex is the experimental one. For a uniaxial tension/compression or pure shear stress state, the enclosed area corresponds to the accumulated strain energy. For a complex stress state, the area corresponds to a part of the strain energy. Note that ER ov can only reflect one valuable side view of a model's capability. When choosing a model, readers should comprehensively evaluate the applicability of a model from their own viewpoint.

Comparison on Micromechanics Models
Among the eight analytical models, Mori-Tanaka method, SCM, bridging model, and Chamis model have been widely used in the analysis of an elastoplastic response with the help of a linearization. For consistency purpose, only the first-moment von-Mises equivalent stress (Section 4.2.1) and the tangent linearization (Section 4.3.2) are incorporated. It should be noted that, for Mori-Tanaka model, bridging model, and self-consistent model, the elastoplastic stiffness of a matrix and the corresponding Eshelby tensor can be anisotropic. However, for Chamis model, the elastoplastic stiffness of a matrix is pre-assumed to be isotropic. Thus, from the viewpoint of consistency, the stiffness of a matrix is isotropicalized (Section 4.4.2) for all the four models. An instantaneous stiffness tensor of the composite by SCM and Mori-Tanaka model are given as Equations (64) and (65) [194,277].
where tan indicates the instantaneous quantity in a tangent form. The subscripts SC and MT denote the quantities from SCM and Mori-Tanaka model, respectively. Regarding bridging model, owing to coupling between normal and shear stresses in a plastic deformation, and elastoplastic bridging tensor is modified to Equations (66)-(69) [278].
The superscript "ep" represents elastoplastic. The off diagonal elements are solved from the condition that the instantaneous compliance tensor of the composite, i.e., M tan ij = M tan ji , is symmetric. M tan ij is given by Equation (70).
As for an elastoplastic Chamis model, Equation (71)-(75) are applied [271]: In general, an instantaneous compliance tensor of the matrix in an elastoplastic Mori-Tanaka model, SCM, and bridging model is anisotropic defined by a classical flow rule. However, an elastoplastic Chamis model is only applicable to cases where the matrix can undergo an isotropic deformation up to failure. Figures 9-11 are the elasto-plastic stress-strain curves of three kinds of UD composites predicted by the four elastoplastic models mentioned above. The experimental data are taken from Kaddour and Hinton [273] and Kawai et al. [279].
As for an elastoplastic Chamis model, Equation (71) In general, an instantaneous compliance tensor of the matrix in an elastoplastic Mori-Tanaka model, SCM, and bridging model is anisotropic defined by a classical flow rule. However, an elastoplastic Chamis model is only applicable to cases where the matrix can undergo an isotropic deformation up to failure. Figures 9-11 are the elasto-plastic stress-strain curves of three kinds of UD composites predicted by the four elastoplastic models mentioned above. The experimental data are taken from Kaddour and Hinton [273] and Kawai et al. [279].    Table 8, it is found that the prediction results of all the four models are too stiff. The yield points and asymptotic tangent moduli provided by the SCM are significantly overestimated for all the cases in Figures 9-11. It is because that, in the configuration of the SCM, a fiber is directly surrounded by a composite medium. Since the composite is much stiffer than the constituent matrix, the plastic deformation of the matrix is underestimated, leading to overestimated stiffness of the composite. In elastic range, the bridging model is recommended. However, regarding yielding, asymptotic modulus and the overall error, all the four model are not satisfactory. More  Table 8, it is found that the prediction results of all the four models are too stiff. The yield points and asymptotic tangent moduli provided by the SCM are significantly overestimated for all the cases in Figures 9-11. It is because that, in the configuration of the SCM, a fiber is directly surrounded by a composite medium. Since the composite is much stiffer than the constituent matrix, the plastic deformation of the matrix is underestimated, leading to overestimated stiffness of the composite. In elastic range, the bridging model is recommended. However, regarding yielding, asymptotic modulus and the overall error, all the four model are not satisfactory. More attempts should be made on modifying the prediction of the yield points and the asymptotic tangent moduli. In addition to a prediction accuracy, the capability on dealing with a tension-shear stress coupling is also critical. From this point, the Chamis model is inferior to the other three models. It is necessary to point out that the comparisons in Figures 9-11 are based on the assumption of an isotropic elastoplastic matrix and the tangent linearization, whereas other factors such as the stress fluctuations and the SCFs are not accounted for. In addition, the averaged error data in Table 8 is useful for readers to evaluate a model but insufficient to tell which model is better. Cautiously speaking, the comparison results are valid only for long fiber reinforced plastic composites with intermediate fiber volume fraction under static, monotonic, and proportional load conditions.

Comparison on Modifications on Yield Stress
As mentioned in Section 4.1, the accuracy in evaluation of a yield point is critical to the prediction capability of an elastoplastic model. A von-Mises equivalent stress is employed to detect yielding. Generally, the equivalent stress is calculated based on the first moment stress (homogenized stress) of a material. However, it is reported that the first moment approach cannot reflect a stress fluctuation, leading to an underestimation of an equivalent stress [205,222,223]. Thus, a second-moment approach is developed [205,222,223]. On the other hand, Huang [103,108,109] pointed out that the homogenized stresses of the matrix in a composite must be converted into true values before an effective property of the composite is evaluated in terms of the fiber and matrix's original property parameters. The conversion is done by multiplying the homogenized quantities with the respective SCFs of the matrix in the composite. Using the true stress concept, the prediction on a yield behavior of a composite can also be significantly improved.
In this section, the tangent linearization is applied to the three approaches, i.e., first moment, second moment, and SCFs. Both the elastoplastic matrix and Eshelby tensor are anisotropic, which are given in Section 4.3.2.

First Moment Approach
For the first moment approach, the J2 flow rule is given as Equations (76)-(78) [114].
where σ 1st eq is the von-Mises equivalent stress based on the first moment approach. σ Y is the yield strength, depending on a work-hardening behavior of a material. The symbol ⊗ denotes a tensor product. I dev is a fourth-order deviatoric tensor. σ and σ are, respectively, the homogenized stresses and the corresponding deviatoric of a material. The key point of the first moment approach is that the von-Mises equivalent stress is calculated from a homogenized stress σ or σ other than the point-wise one σ(x) or σ (x).

Second-Moment Approach
The core concept of the second-moment approach is that a von-Mises equivalent stress is firstly evaluated point-wisely. Then a homogenized von-Mises stress is calculated by volume averaging the obtained point-wise von-Mises stresses (see Equation (79)).
In Equation (79), σ 2nd eq means a von-Mises equivalent stress based on a second-moment approach. The bracket · r represents a volume average manipulation of a phase r⊗ given by Equation (80).
where M r , r = f , m, are the elastic compliance tensors of fiber and matrix, respectively. M is the effective elastic compliance tensor of the composite calculated from a selected homogenization. As addressed by Suquet [223] and Hu [222], the second-moment approach contributes more when the stress/strain fields variation is significant. Theoretically, the field variation in a short fibrous composite is more obvious than that in a long fibrous one. Thus, short fibrous composite is taken as illustration in the comparison between the first and second moment approaches. Doghri et al. [280] conducted a comprehensive study on the second-moment approach based on a kind of aligned short fibrous composite. Figures 12-14 are taken from their work. M-T 1st and 2nd represent Mori-Tanaka model based on the first and second-moment approaches, respectively.
In Equation (79), means a von-Mises equivalent stress based on a second-moment approach. The bracket 〈•〉 represents a volume average manipulation of a phase r⨂ given by Equation (80).
where , = , , are the elastic compliance tensors of fiber and matrix, respectively. is the effective elastic compliance tensor of the composite calculated from a selected homogenization.
As addressed by Suquet [223] and Hu [222], the second-moment approach contributes more when the stress/strain fields variation is significant. Theoretically, the field variation in a short fibrous composite is more obvious than that in a long fibrous one. Thus, short fibrous composite is taken as illustration in the comparison between the first and second moment approaches. Doghri et al. [280] conducted a comprehensive study on the second-moment approach based on a kind of aligned short fibrous composite. Figures 12-14   approach. The bracket 〈•〉 represents a volume average manipulation of a phase r⨂ given by Equation (80).
〈 ( )⨂ ( )〉 = 1 : : where , = , , are the elastic compliance tensors of fiber and matrix, respectively. is the effective elastic compliance tensor of the composite calculated from a selected homogenization.
As addressed by Suquet [223] and Hu [222], the second-moment approach contributes more when the stress/strain fields variation is significant. Theoretically, the field variation in a short fibrous composite is more obvious than that in a long fibrous one. Thus, short fibrous composite is taken as illustration in the comparison between the first and second moment approaches. Doghri et al. [280] conducted a comprehensive study on the second-moment approach based on a kind of aligned short fibrous composite. Figures 12-14 Figure 12, it is found that the second-moment approach can significantly improve the prediction accuracy for a short fiber reinforced polyamide composite under longitudinal tension but makes a low contribution for the case of transverse tension. In order to investigate the reason, FEM plastic strain gradient contours for the composite under the two load cases are plot in Figure 13. It is shown that the strain gradient in the longitudinal tension is much more noticeable than that in the transverse one. The results validate the conclusion of Suquet [223] and Hu [222] that the second-moment approach can make more contributions for cases with significant stress or strain gradients. Please note that the fiber configurations in Figure 13 are not exactly coincident. Theoretically, strain concentration becomes more obviously when the ends of two fibers get closer. But in the second image of Figure 13, the strain distribution is much more homogeneous even at the center-right point where two fiber ends are close to each other. Thus, compared with two coincident images, the images in Figure 13 are more convictive on illustrating the effect of loads on strain distributions. A further comparison is made in Figure 14 to investigate in which situation a secondmoment approach is necessary [280]. Doghri et al. summarized three conditions [280]: (a) the inclusion aspect ratio is larger than 1, (b) the elastic stiffness contrast between a fiber and a matrix is high, (c) work-hardening phenomenon of a matrix is not significant.
When the conditions are satisfied, an application of the second-moment approach can make a significant improvement.   Figure 12, it is found that the second-moment approach can significantly improve the prediction accuracy for a short fiber reinforced polyamide composite under longitudinal tension but makes a low contribution for the case of transverse tension. In order to investigate the reason, FEM plastic strain gradient contours for the composite under the two load cases are plot in Figure 13. It is shown that the strain gradient in the longitudinal tension is much more noticeable than that in the transverse one. The results validate the conclusion of Suquet [223] and Hu [222] that the second-moment approach can make more contributions for cases with significant stress or strain gradients. Please note that the fiber configurations in Figure 13 are not exactly coincident. Theoretically, strain concentration becomes more obviously when the ends of two fibers get closer. But in the second image of Figure 13, the strain distribution is much more homogeneous even at the center-right point where two fiber ends are close to each other. Thus, compared with two coincident images, the images in Figure 13 are more convictive on illustrating the effect of loads on strain distributions. A further comparison is made in Figure 14 to investigate in which situation a second-moment approach is necessary [280]. Doghri et al. summarized three conditions [280]: (a) the inclusion aspect ratio is larger than 1, (b) the elastic stiffness contrast between a fiber and a matrix is high, (c) work-hardening phenomenon of a matrix is not significant.
When the conditions are satisfied, an application of the second-moment approach can make a significant improvement.

SCFs in Matrix
A literature has shown [271] that a predicted matrix-failure controlled strength of a fibrous composite especially under a transverse tension based on the homogenized stresses of the matrix is much lower than the measured counterpart. Liu and Huang [101] pointed out that this phenomenon is resulted from an ignorance of the stress concentrations in the matrix due to introduction of the fiber. They believed that this problem could be tackled by introducing the SCFs to modify the homogenized matrix stresses. The concept of the SCFs in the matrix is originally proposed for a failure prediction based on a bridging model. But it can also be used in an elastoplastic model to estimate the yield behavior of the composite. The SCFs in the matrix of a composite under transverse tension, K t 22 , and compression, K c 22 , are given as Equations (81)-(86) [109].
where σ m u,t and σ m u,c are the tensile and compressive strengths of the matrix, respectively. Owing to transversely isotropic, the SCFs of the matrix in the UD composite along axis 3 are the same as those in the 2nd axis, i.e., K 22 = K 33 , if the transverse stresses of the matrix, σ m 22 and σ m 33 , do not occur simultaneously.
The transverse shear SCF, K 23 , of the matrix is defined according to the Mohr's rule as Equation (87).
where σ m u,s is the shear strength of a matrix. A matrix SCF under in-plane shear is given by Equations (88) and (89), Then, the matrix stress can be modified as Equations (90) and (91).
A comparison is made in Figures 15-17 between the bridging and Mori-Tanaka models with and without the SCFs. The tangent linearization is employed. The involved elastoplastic compliance/stiffness and the corresponding Eshelby's tensors are anisotropic. As mentioned in Section 4, the prediction accuracy of an elastoplastic model depends on three aspects, the elastic modulus, the yield stress, and the asymptotic tangent modulus. From Figures 15-17, it is seen that the prediction error of the bridging model mainly resulted from the overestimated yield stress and the asymptotic tangent modulus. For the cases shown in Figures 15a, 16a and 17, the prediction accuracies are improved significantly by incorporation of the SCFs. However, for the in-plane shear cases, Figures 15b and 16b, the prediction accuracy is still not satisfactory. It is because the introduction of SCFs can modify the prediction on yield stress but cannot reduce errors coming from the overestimated G T asy .
Materials 2018, 11, x FOR PEER REVIEW 26 of 53 modulus, the yield stress, and the asymptotic tangent modulus. From Figures 15-17, it is seen that the prediction error of the bridging model mainly resulted from the overestimated yield stress and the asymptotic tangent modulus. For the cases shown in Figures 15a, 16a and 17, the prediction accuracies are improved significantly by incorporation of the SCFs. However, for the in-plane shear cases, Figures 15b and 16b, the prediction accuracy is still not satisfactory. It is because the introduction of SCFs can modify the prediction on yield stress but cannot reduce errors coming from the overestimated .  modulus, the yield stress, and the asymptotic tangent modulus. From Figures 15-17, it is seen that the prediction error of the bridging model mainly resulted from the overestimated yield stress and the asymptotic tangent modulus. For the cases shown in Figures 15a, 16a and 17, the prediction accuracies are improved significantly by incorporation of the SCFs. However, for the in-plane shear cases, Figures 15b and 16b, the prediction accuracy is still not satisfactory. It is because the introduction of SCFs can modify the prediction on yield stress but cannot reduce errors coming from the overestimated .  The prediction errors of the Mori-Tanaka model come from all the three aspects, the E or G, the , and the or . For most cases, the introduction of SCFs can improve its prediction results. But for some cases, the prediction accuracy becomes worse. It is because the yield stress is underestimated by the Mori-Tanaka model with SCFs. In addition, the errors from the overestimated or cannot be modified by introducing the SCFs, similarly as the bridging model.

Comparison on Linearization
The next step to build an elastoplastic model is in selection of a linearization. For simplicity, only the matrix is treated as elastoplastic while the fiber keeps elastic. Figure 18 shows a schematic of the four most widely used linearization schemes, i.e., the secant linearization, the tangent linearization, the TFA, and the affine formulations. The latest linearization, an incremental-secant scheme, is also introduced later in this section. The Mori-Tanaka model and the first-moment equivalent stress are employed for all the five linearizations. The prediction errors of the Mori-Tanaka model come from all the three aspects, the E or G, the σ Y , and the E T asy or G T asy . For most cases, the introduction of SCFs can improve its prediction results. But for some cases, the prediction accuracy becomes worse. It is because the yield stress is underestimated by the Mori-Tanaka model with SCFs. In addition, the errors from the overestimated E T asy or G T asy cannot be modified by introducing the SCFs, similarly as the bridging model.

Comparison on Linearization
The next step to build an elastoplastic model is in selection of a linearization. For simplicity, only the matrix is treated as elastoplastic while the fiber keeps elastic. Figure 18 shows a schematic of the four most widely used linearization schemes, i.e., the secant linearization, the tangent linearization, the TFA, and the affine formulations. The latest linearization, an incremental-secant scheme, is also introduced later in this section. The Mori-Tanaka model and the first-moment equivalent stress are employed for all the five linearizations. the matrix is treated as elastoplastic while the fiber keeps elastic. Figure 18 shows a schematic of the four most widely used linearization schemes, i.e., the secant linearization, the tangent linearization, the TFA, and the affine formulations. The latest linearization, an incremental-secant scheme, is also introduced later in this section. The Mori-Tanaka model and the first-moment equivalent stress are employed for all the five linearizations.

Secant Linearization
On the secant linearization, it is assumed that the total stresses and strains can be correlated by a "secant" stiffness or compliance tensor through Equations (92) and (93).

Secant Linearization
On the secant linearization, it is assumed that the total stresses and strains can be correlated by a "secant" stiffness or compliance tensor through Equations (92) and (93).
where the superscript sec represents instantaneous quantities in a secant form. The global strain concentration tensor can also be rewritten in a secant form as Equation (94).
The effective secant stiffness tensor is given by Equation (95).
If the Mori-Tanaka model is employed, the B sec f is obtained as Equations (96) and (97).
The non-zero elements of the M sec m are given by Equations (98) where E sec m , G sec m , and υ sec m are the "secant" elastic constants. The secant Young's modulus E sec m is directly obtained from a given uniaxial tension or compression stress-strain curve. Supposing that the bulk modulus keeps elastic, the secant Poisson's ratio and the secant shear modulus are obtained, respectively, as Equations (101) and (102).
Equations (103) It should be noted that Equations (103)-(105) only give the expressions of S sec m for cylindrical fiber reinforced composites. For another shape of the inhomogeneity, the S sec m can also be obtained directly by replacing the elastic properties of the matrix in an elastic Eshelby tensor with the corresponding secant ones.

Tangent Linearization
Equations (106)-(108) are derived from a tangent linearization: where dε, dε r and dσ, dσ r are the strain and stress tensors of the composite, and the constituent phases, respectively. The effective compliance tensor is derived as Equation (109).
Also, if Mori-Tanaka model is employed, the B tan f is obtained as Equations (110) and (111).
The superscript tan represents tangent quantities. Several approaches can be used to determine the M tan m and S tan m . In this section, only the most general one is highlighted. When a material undergoes an elasto-plastic deformation, its mechanical properties become anisotropic. If a J2 flow rule is employed, the tangent instantaneous compliance tensor of a matrix is given as Equations (112) M m is the elastic compliance tensor of a matrix. σ eq is an von-Mises equivalent stress given in Equation (77). σ is a stress deviator. E tan m is a tangent modulus obtained from a uniaxial tension or compression stress-strain curve. The corresponding tangent Poisson's ratio and the tangent shear modulus are given by Equations (115) and (116).
Owing to anisotropy of the matrix in an elasto-plastic region, the corresponding tangent Eshelby tensor has to be calculated by a numerical integration (see Equations (117)-(126)) [230].
cos ω (125) The superscript "m-tan" denotes quantities of the matrix in the tangent form.

Transformation Field Analysis (TFA)
Dvorak [208,281] proposed a transformation field analysis to evaluate elastoplastic behaviors of a composite, in which the plastic strain of a matrix is viewed as an eigenstrain. The stress-strain relations of the composite and the constituent phases are presented in Equations (127)- (130).
where dλ, dλ r , and dε pla , dε where A r and B r are the global strain and stress concentration tensors, respectively. A T r and B T r are corresponding transposed tensors. D rs and F rs are the influence tensors given by Equations (135)-(138), It is reported by Chaboche [57,95] that the prediction results of elastoplastic behavior of composites provided by the traditional TFA are too stiff. Modifications to Equations (131)-(132) were proposed with the correction tensors K r , r = f , m, as shown in Equations (139) and (140) [57,95].
The correction tensors can be obtained by solving Equation (141).
where M tan s , L tan s are the instantaneous tangent compliance and stiffness tensors of the constituents. Since the fiber is linearly elastic, only M tan m needs to be updated. A tan r are the instantaneous strain concentration tensors calculated from M tan s , L tan s , and the related Eshelby tensor. As indicated by Chaboche [57], Equation (141) is indeterminate. But it can be solved by choosing K s = I for the stiffer constituent phase. For a two phase composite with the elastic fibers as reinforcement, we can get Equations (142) and (143).

Affine Formulation
Zaoui and Masson [219] and Masson et al. [218] proposed a new linearization, namely an affine formulation, to improve the prediction accuracy of the viscoplastic behavior of a composite. Chaboche [95] presented a compact version of the affine formulation for an elasto-plastic case. Equations (144)- (147) are the corresponding formulations.
where τ, τ r and η, η r are the pre-stresses and pre-strains of the composite and its constituents. The tangent instantaneous compliance tensor M tan m is given in Equation (112), and L tan m = M tan m −1 .
The instantaneous effective compliance and stiffness tensors M tan and L tan are obtained from a selected homogenization model. The localization equations are given by Equations (148)- (150).
In Equation (148), S tan is given by Equation (117). The instantaneous concentration tensors A tan r and B tan r can be calculated from a selected homogenization model, e.g., Mori-Tanaka model.

Incremental-Secant Scheme
Wu et al. [112] proposed a novel incremental-secant linearization for elasto-plastic responses of a composite. It is claimed that the linearization is applicable to non-monotonic and non-proportional loads [112][113][114]. Figure 19 illustrates the idea of the linearization.
Zaoui and Masson [219] and Masson et al. [218] proposed a new linearization, namely an affine formulation, to improve the prediction accuracy of the viscoplastic behavior of a composite. Chaboche [95] presented a compact version of the affine formulation for an elasto-plastic case. Equations (144)- (147) are the corresponding formulations.
where , and , are the pre-stresses and pre-strains of the composite and its constituents. The tangent instantaneous compliance tensor is given in Equation (112), and = ( ) . The instantaneous effective compliance and stiffness tensors and are obtained from a selected homogenization model. The localization equations are given by Equations (148)- (150).
In Equation (148), is given by Equation (117). The instantaneous concentration tensors and can be calculated from a selected homogenization model, e.g., Mori-Tanaka model.

Incremental-Secant Scheme
Wu et al. [112] proposed a novel incremental-secant linearization for elasto-plastic responses of a composite. It is claimed that the linearization is applicable to non-monotonic and non-proportional loads [112][113][114]. Figure 19 illustrates the idea of the linearization.  Consider a composite at the nth load step with a total stress and strain tensors, σ n and ε n . The corresponding stress and strain increments are ∆σ n+1 , ∆ε n+1 . Then the total stress and strain tensors in the next step are given by Equations (151) and (152).
The subscripts n and (n + 1) denote quantities of the nth and (n + 1)th load steps, respectively. In the incremental-secant scheme, a fictitious elastic unload process is introduced at the nth load step. It should be noted that the residual stresses in constituent phases can be nonzero due to the heterogeneity even no homogenized residual stress exists in the composite. From a fictitious unload state (σ res n , ε res n ), a secant reload procedure is applied from the nth to the (n + 1)th step. Thus, Equations (153)- (155) are obtained. ε n+1 = ε res n + ∆ε rel n+1 (153) σ n+1 = σ res n + ∆σ rel n+1 (154) ∆σ rel n+1 = L sec n+1 ∆ε rel n+1 (155) ε res n and σ res n are the residual strains and stresses after unload. L sec n+1 is the secant stiffness tensor at the (n + 1)th load step, viewing (σ res n , ε res n ) as a starting point. ∆ε rel n+1 and ∆σ rel n+1 are the reload strain and stress increments. It is assumed that at the nth step, all the quantities required are known. Given a total strain ε n+1 for the next load step, the strain increment ∆ε rel n+1 can be obtained from Equation (153). Consequently, the reload stress increment ∆σ rel n+1 is reached by Equation (155), and the total stress σ n+1 is obtained from Equation (154). It should be noted that the homogenized residual stress of the composite is zero, namely σ res n = 0. Thus, in Equation (155), the only unknown argument is the effective secant stiffness tensor L sec n+1 . It can be obtained from the secant properties of the constituents through Equations (156)- (157).
The secant stiffness tensor L i−sec n+1 and the secant global strain concentration tensor A i−sec n+1 depend on the stresses of the phases σ i n+1 . Thus, L i−sec n+1 , A i−sec n+1 , and σ i n+1 can be obtained by solving Equations (156) and (157) iteratively. Then the homogenized stresses of the composite σ n+1 can be determined from Equations (154) and (155), in which the L sec n+1 is evaluated from L i−sec n+1 and A i−sec n+1 by a selected homogenization theory. Alternatively, σ n+1 can also be obtained from Equation (158). (158) To continue the calculation to next load step, it is necessary to get the residual stresses and strains of the composite and the constituents. The residual stresses of the composite are zero, and the corresponding residual strains are given by Equations (159) and (160): where ∆ε unl n+1 are the unload strains shown in Figure 19. M ela is the elastic effective compliance tensor of the composite evaluated by a selected homogenization theory. The residual strains and stresses of the constituents are given by Equations (161) and (162).
The unload strains of the constituents are found as Equation (163).
In this way, the calculation process continues.

Quantitative Comparison on the Linearizations
Before a quantitative study, the applicability of the five linearization models should be briefly introduced. Applicability of a linearization model can be assessed in three aspects, i.e., whether they are applicable to non-monotonic and non-proportional loads, whether a tension-shear coupling can be considered, and whether a numerical integration on an Eshelby tensor can be avoided. Table 9 summarizes the applicability of different linearization models. Essentially, if a linearization is in an incremental form, it is applicable to non-monotonic and non-proportional loading cases. Rewriting an elastoplastic compliance tensor into four blocks, a non-zero off-diagonal block implies that the model can account for a tension-shearing coupling. The third aspect depends on whether Equation (117) is used. A quantitative comparison is made taken experimental results of E-glass/Epoxy, IM7/8551-7, and AS4/Peek UD composites as benchmark. The matrix behavior in the secant and the incremental-secant linearizations is taken as isotropic, whereas that in the other four linearizations can be either isotropic [212] or anisotropic. For a consistency in the comparison, the matrix is assumed to be elasto-plastically isotropic. Each linearization is incorporated with Mori-Tanaka model. Results are shown in Figures 20-22. Again, the comparison results are restricted to static, monotonic, and proportional load cases. In this way, the calculation process continues.

Quantitative Comparison on the Linearizations
Before a quantitative study, the applicability of the five linearization models should be briefly introduced. Applicability of a linearization model can be assessed in three aspects, i.e., whether they are applicable to non-monotonic and non-proportional loads, whether a tension-shear coupling can be considered, and whether a numerical integration on an Eshelby tensor can be avoided. Table 9 summarizes the applicability of different linearization models. Essentially, if a linearization is in an incremental form, it is applicable to non-monotonic and non-proportional loading cases. Rewriting an elastoplastic compliance tensor into four blocks, a non-zero off-diagonal block implies that the model can account for a tension-shearing coupling. The third aspect depends on whether Equation (117) is used. A quantitative comparison is made taken experimental results of E-glass/Epoxy, IM7/8551-7, and AS4/Peek UD composites as benchmark. The matrix behavior in the secant and the incrementalsecant linearizations is taken as isotropic, whereas that in the other four linearizations can be either isotropic [212] or anisotropic. For a consistency in the comparison, the matrix is assumed to be elastoplastically isotropic. Each linearization is incorporated with Mori-Tanaka model. Results are shown in Figures 20-22. Again, the comparison results are restricted to static, monotonic, and proportional load cases.   Since all the five linearizations are combined with Mori-Tanaka model, errors induced by the predicted elastic moduli and yield points cannot be reflected. As for the predicted asymptotic tangent modulus, or , it is seen that the affine formulation delivers a relative soft response, whereas the TFA's predictions are stiffer for the cases in Figure 22a,b. The other four theories present similar results. Overall, the results from all of the five theories are stiff compared with the experiments, especially for the in-plane shear cases.

Comparison on Modifications of a Plastic Eshelby Tensor
Generally speaking, the Eshelby tensor for an anisotropic material should be obtained from a numerical integration. However, some works reported that the precise Eshelby tensor would Since all the five linearizations are combined with Mori-Tanaka model, errors induced by the predicted elastic moduli and yield points cannot be reflected. As for the predicted asymptotic tangent modulus, E T asy or G T asy , it is seen that the affine formulation delivers a relative soft response, whereas the TFA's predictions are stiffer for the cases in Figure 22a,b. The other four theories present similar results. Overall, the results from all of the five theories are stiff compared with the experiments, especially for the in-plane shear cases.

Comparison on Modifications of a Plastic Eshelby Tensor
Generally speaking, the Eshelby tensor for an anisotropic material should be obtained from a numerical integration. However, some works reported that the precise Eshelby tensor would overestimate the instantaneous stiffness of a composite [95,[210][211][212]. Besides, the numerical integration is time consuming compared with an explicit one. Thus, several modifications on a plastic Eshelby tensor are proposed, e.g., an isotropic matrix approach, an isotropic Eshelby tensor approach, and Peng's approach. Including the numerical integration, which is called an anisotropic Eshelby tensor approach, the four methods are compared.
In this section, all the modifications on Eshelby tensor are incorporated with Mori-Tanaka model, first-moment equivalent stress, and tangent linearization.

Anisotropic Eshelby Tensor Approach
In this approach, an elastoplastic compliance or stiffness tensor of the matrix is evaluated by the J2 flow rule shown in Equations (112)- (114). The corresponding Eshelby tensor is obtained from Equation (117), and the effective properties of the composite are derived from Equation (109).

Isotropic Matrix Method
González and LLorca [212] demonstrated that an elastoplastic tangent stiffness tensor of the matrix can be expressed as isotropic if the composite is subjected to an asymmetrically proportional load. In such a case, M tan m and S tan can be directly obtained by replacing the secant moduli E sec m , G sec m , and υ sec m in Equations (98)-(105) with the tangent ones, E tan m , G tan m , and υ tan m .

Isotropic Eshelby Tensor Method
It was suggested by Doghri and Ouaar [65] that a better evaluation of elastoplastic behavior of a composite could be achieved if the elasto-plastic response of the matrix remain anisotropic, but the corresponding Eshelby tensor is defined through an isotropic manner. Specifically, the expressions of M tan m are the same by Equations (112)- (114). The Eshelby tensor is obtained by replacing E sec m , G sec m , and υ sec m in Equations (103)- (105) with E tan m , G tan m , and υ tan m .

Peng's Approach
Peng et al. [116] presented a new method to determine an Eshelby tensor for the elastoplastic behavior of a composite. In Peng's approach, a reference elastic medium is introduced, whose configuration and properties are identical to the elastoplastic matrix in the composite. It is assumed that the elastoplastic behavior of the composite can be characterized by two kinds of eigenstrains. One is induced by the inhomogeneity, and the other by the plastic deformation of the matrix. Based on this, a new determination of the Eshelby tensor is given by Equation (164) [116].
where M m is an elastic compliance tensor of the matrix, and S m is the corresponding elastic Eshelby tensor. L tan m is an instantaneous tangent stiffness tensor of the matrix. If Mori-Tanaka model is employed, the effective compliance tensor of the composite can be obtained by replacing S tan m in Equations (109)- (111) with S * shown in Equation (164).

Quantitative Comparison on Modifications of the Eshelby Tensor
Features of the four approaches to elasto-plastic behavior of a composite are summarized in Table 10. Amongst, the isotropic Eshelby tensor and Peng's approaches can meet all the three aspects mentioned previously. Quantitative comparisons for the four approaches are shown in Figure 23.  From Figure 23a-f, we can see that the anisotropic Eshelby tensor approach presents much stiffer responses than the other three. However, the isotropic Eshelby tensor approach resulted in the − curves in Figure 23a,d-f to be somewhat physically unacceptable. Table 11 lists the stresses of the constituent fiber and matrix predicted by the four approaches under some loads. From Figure 23a-f, we can see that the anisotropic Eshelby tensor approach presents much stiffer responses than the other three. However, the isotropic Eshelby tensor approach resulted in the σ 11 − ε 22 curves in Figure 23a,d-f to be somewhat physically unacceptable. Table 11 lists the stresses of the constituent fiber and matrix predicted by the four approaches under some loads.  Table 11 indicates that the anisotropic Eshelby tensor, the isotropic matrix, and Peng's approaches give similar predictions for the constituent stresses. However, the signs of the majority stress components predicted by the isotropic Eshelby tensor model are changed and the corresponding magnitudes are much larger than those by the three other approaches. Furthermore, a predicted tri-axial tension and tri-axial compression may be attained by the fiber and matrix, respectively, even though a composite is subjected to a longitudinal tension. Due to these, a further development in the isotropic Eshelby tensor model is expected.
Also, the σ 11 − ε 22 curves predicted by Peng's approach in Figure 23b looks quite different from the experiments. Besides, the tendency of the predicted curves by Peng's approach in other cases was consistent with the experiments. Thus, Peng's approach can be seen a reliable method to a large extent. But attention should be paid to the Peng's approach when dealing with a Poisson's effect of a composite under a transverse load.
As mentioned in Section 4, the prediction accuracy can be improved by making use of a second-moment estimation or SCFs of the matrix. Both may be efficient to decrease the too stiff response predicted by the anisotropic Eshelby tensor approach. Considering a reliable applicability, the anisotropic Eshelby tensor approach is still recommended over the three other approaches.

Conclusions
A review and comparative study is made for various elastoplastic micromechanics models of a composite. The comparative study has been carried out regarding four aspects. They are the selection of a homogenization, the modifications on yield stress, the selection of a linearization, and the determination of an Eshelby tensor. Some conclusions can be drawn as follows.
In an elastic range, bridging model performs the best. In an elastoplastic range, Mori-Tanaka model, SCM, bridging model, and Chamis model are often found applicable. Based on an isotropic matrix assumption and a tangent linearization, all of the four models deliver stiffer results for elastoplastic stress-strain curves of some composites. The deviation sources can be attributed to threefold, i.e., the elastic modulus, the yield point, and the asymptotic tangent modulus.
The deviation resulted from an overestimated yield stress can be reduced by using a second stress estimation or the SCFs of the matrix in the composite. The former is only useful when the stress field fluctuation in the composite is significant. The latter can be applied to a more general case. Further, a combination of the elastoplastic bridging model and the SCFs is recommended, owing to its applicability for non-monotonic and non-proportional loads, higher prediction accuracy, and computational efficiency with no evaluation of an Eshelby tensor. Nevertheless, both of the two modifications are less efficient for an overestimated asymptotic tangent modulus.
Different linearizations can give different evaluations on the asymptotic tangent modulus of an elasto-plastic response of the composite, but not significant in general. Thus, attention should be focused on when to apply a linearization. More applications and developments on the incremental-secant linearization are expected because it is efficient in computation, applicable for complicated loads, good in accuracy, and flexible in combination with various homogenizations.
Regarding four kinds of determinations for the Eshelby tensor in an elasto-plastic range, the anisotropic Eshelby tensor approach delivers much stiffer results. The other three approaches present better. But, the isotropic-matrix approach has only been demonstrated applicable to asymmetrically proportional loads. The remaining two approaches need further development on their prediction capability on a Poisson's effect.
Unfortunately, most of the methods reviewed in this work cannot modify an overestimation on the asymptotic tangent modulus to a satisfactory extent, especially for in-plane shear cases. This might be because a perfect fiber/matrix interface bonding has been incorporated with all of the methods. Furthermore, the plastic behavior parameters of a matrix have been determined only from a uniaxial tensile stress-strain curve, defined by a so-called single-parameter plasticity theory. A two-parameter plasticity theory, i.e., both measured uniaxial tensile and in-plane shear stress-strain curves of the matrix are used to define its plastic behavior parameters, would be more pertinent. It is expected that better predictions for elasto-plastic responses of a composite can be achieved from development in both micromechanics models for the composite and plasticity theories for an isotropic matrix.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Appendix A
The following tables show the detail information of the experimental data and simulation results of different models.