Homogenization and Localization of Ratcheting Behavior of Composite Materials and Structures with the Thermal Residual Stress Effect

In this contribution, the ratcheting behavior and local field distribution of unidirectional metal matrix composites are investigated under cyclic loading. To that end, we extended the finite-volume direct averaging micromechanics (FVDAM) theory by incorporating the rule of nonlinear kinematic hardening. The proposed method enables efficient and accurate simulation of the ratcheting behavior of unidirectional composites. The local satisfaction of equilibrium equations of the FVDAM theory facilitates quick and rapid convergence during the plastic iterations. To verify the proposed theory, a finite-element (FE) based unit cell model is constructed with the same mesh discretization. The remarkable correlation of the transverse response and local field distribution generated by the FVDAM and FE techniques demonstrates the effectiveness and accuracy of the proposed models. The stress discontinuities along the fiber/matrix interface that are generic to the finite-element theory are absent in the FVDAM prediction. The effects of thermal residual stresses induced during the consolidation process, as well as fiber orientations, are revealed. The generated results indicate that the FVDAM is well suited for simulating the elastic-plastic ratcheting behavior of metal matrix composites, which will provide the conventional finite-element based technique with an attractive alternative.


Introduction
Metal matrix composites have become and will continue to be one of the most attractive alternatives to conventional materials for many applications in engineering such as turbine engine components, pressure pipelines and vessels. The applications often leave the composites exposed to complicated cyclic loading, which leads to cyclic deformations including ratcheting and endangers safety in the operation. For the safety and full use of the metal matrix composites, a constitutive relation needs to be established to characterize the ratcheting effect of the metal matrix composites. The properties of composites are usually decided by constituent material properties and microstructure parameters. These factors make it very difficult to characterize the mechanical behavior. In addition, during the fabrication of the composites, materials are usually consolidated at a high temperature which induces the residual stresses when being cooled to room temperature. The residual stresses may change the material properties and have a great influence on the thermo-mechanical response. Therefore, the study of the thermal residual effect is critical to the accurate prediction of the ratcheting effect in the metal matrix composites.
Ratcheting is a cyclic accumulation of inelastic deformation, which occurs to materials during cyclic loading with no-zero mean stresses. It has been a critical factor in fatigue life prediction and structure designs. However, simulating ratcheting is a complicated task since the inelastic deformation increases with loading cycle. The ratcheting of monotonic homogeneous metal materials is investigated extensively; one of the most famous cyclic plasticity models is the Armstrong and Frederick model (AF model) [1]. The terms of dynamic recovery of the AF model enable the calculation of the ratcheting effect. However, the dynamic recovery terms in this model are too active, which leads to overprediction of the ratcheting effect. Many researchers have extended the AF model to improve the precision of the prediction, such as Ohno and Wang [2][3][4], Tanaka [5], Chaboche [6][7][8], Abdel-Karim and Ohno [9] and Kang [10]. Details of these works can be seen in the review papers by Ohno [11], Bari and Hassan [12], Chaboche [8], Kang [13], Sajjad [14] and Chen [15]. Meanwhile, experimental characterizations of the cyclic deformations of metal matrix composites have been conducted by researchers, cf., Zhang et al. [16], Kang et al. [17], Han et al. [18] and Kimmig et al. [19]. Compared with monotonic homogeneous metal materials, the properties of the metal matrix composites are decided by off-axis angles, volume fractions and array patterns. Therefore, they must be tested repeatedly for each variation of the composites, which makes the experimental characterization of composite materials a more tedious and expensive task.
Numerical modeling of ratcheting of the composites has been investigated by many researchers during the past decades [20][21][22][23][24][25]. The classical micromechanics models under the framework of the Eshelby theory [26] and the Mori-Tanaka model [20,21] were built to simulate the mechanical response of metal matrix composites. However, classical micromechanics is based on the total deformation theory and cannot simulate the ratcheting effect of composite materials. Therefore, some researchers have built the model based on tangent incremental linearization theory [27] to simulate the elastoplastic response and additive scheme [28] or variational approaches [29] for the viscoplastic response of composites. Conventionally, the finite-element method (FEM) is the most commonly used method to simulate the elastoplastic response of the composites. However, care must be taken since the large modulus contrast near the fiber and matrix interface can produce stress concentrations which requires a more detailed mesh for FEM to avoid stress discontinuity at the interface.
While the FVDAM theory has been studied extensively, little work has been reported concerning the cyclic plasticity. The semi-analytical framework, quick convergence and excellent stability of the FVDAM theory make it particularly suitable for analyzing the cyclic behavior of metal matrix composites.
Therefore, in this research, the FVDAM is further extended to incorporate the ratcheting plasticity model, which allows efficient and accurate simulations of the thermo-mechanical cyclic behavior of the metal matrix composites. The proposed model not only enables the prediction of global ratcheting effects but also provides local information of the microstructures of the metal matrix composites. The validation of the model is provided by FEM simulations. The influence of thermal residual stresses, as well as off-axis angles, is discussed.
The main contributions of this research are listed as follows: (1) The cyclic plastic constitutive model is implemented successfully into the FVDAM theory for the first time, giving the FVDAM the ability to simulate the cyclic behavior of the metal matrix composites. (2) Not only the homogenized cyclic response of the composites but also the local field distributions are obtained using the proposed model, since most of the mean-field and nongeometrical models cannot obtain the local field response.
(3) The effects of thermal residual stresses induced during the consolidation process, as well as fiber orientations, are considered, which makes the calculation of the cyclic behavior closer to the actual working condition. (4) The effectiveness of the new micromechanical model to accurately simulate both the homogenized stress-strain curve and local field distribution is verified by an FEM model with the same mesh discretizations, providing a golden standard for the simulation. The results also show that the proposed model has a better convergence near the stress concentration region.
This paper is presented as follows: Section 2 gives an introduction of the framework of the parametric finite volume direct averaging micromechanics (FVDAM). The process of implementing the Abdel-Karim-Ohno constitutive model into parametric FVDAM theory is demonstrated in Section 3. Section 4 shows the homogenized and local response simulated by the proposed model and provides the validation with FEM theory.

Finite-Volume Direct Averaging Micromechanics (FVDAM)
The displacement equations of the parametric FVDAM are presented by a multiscale expansion. The index of the quadrilateral subvolumes are noted as (q). In this case, the study of unidirectional composites can be seen as a generalized plane strain problem, where fiber direction y 1 cannot affect stress components in other directions. Therefore, a two-dimensional unit cell can be used to analyze this problem (see Figure 1). The number of vertices (y ] for the face F m are: The qth quadrilateral subvolume is generated by the transformation from the reference coordinates η − ξ to the actual discretized microstructure (see Figure 2): where  The displacement equation in each qth subvolume is generated by using the two-scale expansion in the slow scale x and fast scale y: which produces the local strain-displacement relations in the actual coordinates y 2 − y 3 with the form: ij are the fluctuating strain components. The second order Legendre expansion is applied to approximate the u (q) where W (q) i(mn) are unknown coefficients which can be expressed by displacement and tractions in the qth subvolume.
The interfacial displacements and tractions on the pth face are given below: For face (p = 1, 3): where t , from Hook's law.
In the FVDAM theory, the inverse of the volume-averaged JacobianĴ is applied to mapping the surface-averaged interfacial displacement in the actual coordinates y 2 − y 3 onto the reference coordinates η − ξ. The surface-averaged interfacial displacement in η − ξ can be expressed by the first and second order coefficients in representation of the displacement field; these coefficients can be written as W i = [W i(10), W i(01), W i(20), W i(02) ] T . Therefore, we have: where By using the surface-averaged subvolume equilibrium equations: where L p represents the length of the pth face, the expressions of W i in the qth subvolume can be expressed byû (q) . By using the relation between traction and fluctuating displacements in a surface-average sense, the following equation can be derived: show the orientation of each of the four faces of the qth subvolume, and Γ (q) and G (q) contain vectors concerning thermal and plastic contributions. Each element inT (q) andÛ (q) has three components, ij is a 3 × 3 matrix derived by the geometric and mechanical properties in each subvolume.
The global stiffness matrix can then be assembled using the continuity equations: Then the global system of equations can be formed after applying periodic boundary conditions and eliminating the singularity of the global stiffness matrices. The resulting systems are expressed below: where ∆C contains the differences of adjacent local stiffness matrices, Γ and G represent the thermal effect and plastic information. The fluctuating surface-averaged interfacial displacementsÛ can be determined by solving Equation (13). The coefficients W i can be calculated subsequently, which enables the calculation of local fields in every subvolume. The results ofÛ can then be used to derive the Hill's strain concentration matrix A (q) by using the relation below: where the components of A (q) are derived by applying only one component of macroscopic strain at a time. D (q) is a vector concerning thermal and plastic contributions, which can be presented by macroscopic strainsε at each loading step. The homogenized Hook's law can be calculated by taking the volume average of the stress field as below:σ where v (q) = V (q) /V represents the volume fraction.
The macroscopic constitutive equations for the composite material can be derived using the localization relations and stress-strain relations in the volume-averaged sensē represents the homogenized stiffness matrix; the volume-averaged thermal and plastic stresses can be derived below The vector G in Equation (13) denotes plastic contributions, which are an integration of plastic strains on the surface of each subvolume. The total plastic strain at each material point consists of the contribution of the previous step and the contribution of current load increment In order to give the FVDAM model the capability of simulating the ratcheting effect, a cyclic plastic constitutive relation should be implemented in this model when calculating the plastic strain increment dε pl(q) (η, ξ) in Equation (18). The detail of the implementation is described in Section 3.

Implementation of the Cyclic Plasticity Model
After the discretization in Section 2, the aim of the model becomes the calculation of the plastic strain incremental of each and every material point in Equation (18). In this paper, we introduce the Abdel-Karim-Ohno model [9] into the FVDAM theory, since the Abdel-Karim-Ohno model is a low-order nonlinear constitutive model which makes it easy to implement in the algorithm [46]. We begin with a brief summary of the Abdel-Karim-Ohno model. For the elastoplastic problems, the small strain ε is assumed to be a combination of elastic deformation ε e and plastic deformation ε pl and back stresses α are divided into N parts as shown below: Then, by using Hook's law and the flow rule related to Mises-type yield surface F = 0, the function of ε e and ε pl in Equation (19) can be written as follows: where (:) means inner product calculation, C is the stiffness matrix,λ can be calculated usingḞ = 0 and F = 0, s and a are deviatoric parts of σ and α. Therefore, in the deviatoric space, a is the center of the yield surface and Y represents the radius of the yield surface. Normally, a and Y will move and expand during the inelastic deformation. In this paper, we consider kinematic hardening, where the radius of the yield surface Y is constant and the center of the yield surface a moves. Since a is the deviatoric part of α, a can also be decomposed into N different parts: The evolution of a is described as below: where ζ i is a material parameter, µ i is the weighing coefficient ranging from 0 to 1, µ i = 1 corresponds to the AF model, µ i = 0 corresponds to the first version of the Ohno-Wang model, r i can be determined by the unidirectional tensile test and denotes Macaulay's brackets (i.e., For numerical implementation, the constitutive equations (Equations (19)- (25)) should be discretized over time t. The problem can then be described as follows: Knowing the variables of every point in the qth subvolume in the reference coordinates (i.e., σ (q) ) at t n and giving the strain increment dε (q) (η, ξ), try to find the plastic strain increment dε pl(q) (η, ξ) at t n+1 , which can be solved by applying the backward Euler method and the return mapping method [46]. After the plastic strain increment dε pl(q) (η, ξ) is calculated, Equation (18) can be solved and vector G in Equation (13) can be derived, which allows the updating of D (q) in Equation (14) and the follow-up calculation of the homogenized and local response of the composites under cyclic loading.
The whole algorithm is illustrated in Figure 3. The outer loop controls the macroscopic loading path. We start with a thermal cooling process (∆T = −550 • C), which has been discreted into n parts [∆t, ..., n∆t = ∆T]. Then, we apply the macroscopic mechanical loading, which has been discreted into m parts [∆ε, ..., m∆ε].
The inner loop demonstrates the process of obtaining the local and homogenized response of the composite material under thermal and cyclic mechanical loading. The effective thermal expansion constants are used to generate thermal residual stressesσ th of the composites to simulate the cooling process. After thermal loading, cyclic mechanical loading is applied and the cyclic plasticity constitutive equation is implemented to calculate cyclic plastic strain ε pl in the qth subvolume. When the plastic strain in all subvolumes meets the requirement of convergence, the local response can be obtained. Then the homogenized response can be updated using the volume-averaging method.
During the calculation of the plastic strain in the qth subvolume, we applied the Newton-Raphson method, the convergence of which is quadratic [47]. The convergence criterion is written in terms of the accumulated plastic strain increment ∆p (q) i+1 as follows: where k denotes the number of local iterations and i represents the loading step.
The numerical successive method we use is similar to the one proposed by Kaboyashi and Ohno [46]. The detailed demonstration of convergence can also be found in [46].

Unidirectional Transverse Loading
In order to verify the effectiveness of the FVDAM theory to accurately predict both the homogenized stress-strain curve and local field distribution, a finite-element unit cell model is built with the same mesh discretizations, providing a golden standard for the FVDAM simulation.
In this paper, we simulated the strain-controlled unidirecional loading experiment, which means we control the macroscopicε 22 . However, other components ofε are not zero. As in the unidirectional experiment, we only load on direction 2 and other directions are not constrained on the material. Therefore, knowingε 22 and the global stiffness matrix, other components ofε can be calculated using Hook's law.
To simulate the response of the composite materials with thermal residual stresses, a square unit cell with a fiber was built. The volume fraction of fiber is 0.3. The material properties are listed in Tables 1 and 2. The large modulus contrast and the strong fiber-fiber interaction at this volume fraction produce highly localized stress fields, which are very demanding for the accuracy of the method. We first compare the local stress distributions within the composite microstructure after a cool down temperature of ∆T = −550 • C simulated by the FVDAM and FEM analyses (see Figure 4). The thermal expansion constants α * = C * −1 ∑ N q q=1 ν (q) (A (q) ) T Γ (q) are obtained using the results of Levin [48]. Because the composites are not constrained during the consolidation process, the macroscopic stresses in Equation (16) are taken as zeros. To generate the smooth thermal stress fields, the unit cells were discretized into 96 × 96 subvolumes/elements in both approaches. We note that the thermal residual stresses of the present material systems generated by both approaches are in good agreement, illustrating the excellent predictive capability of the FVDAM theory.  To illustrate the effect of thermal residual stresses, Figure 5 compares the transverse thermo-mechanical tensile response of the composite with different values of µ = 0, 0.15, 1. The temperature changes are prescribed as ∆T = −550 • C. The stress-strain response generated without thermal residual stress effects (∆T = 0 • C) are included in the figures for comparison. For verification purposes, the corresponding FEM results under the same loading conditions are enclosed in the figures as well. As observed, the correlations between the two methods are in good agreement. The thermal residual stresses provide a slight enhancement in homogenized stresses in the elastic-plastic region while the elastic response almost remains the same. A close look at the figures also indicates that the initiation of plasticity in the macroscopic response occurs earlier (or at a lower stress state) due to the significant tensile stresses induced by the thermal cooling process prior to the mechanical loading. The local field distributions considering thermal effects for µ = 0.15 at the applied strain of ε 22 = 1% are shown in Figure 6. The results of FVDAM correlates well with the FEM prediction and no apparent difference is observed. The FEM's prediction, however, shows small discontinuities in the transverse normal stress field in the localized region near the fiber/matrix interface (see Figure 6b), which is not observed in the FVDAM results. This is because the large modulus contrast produces the large stress concentration and deformation gradient at the fiber/matrix interface. The FEM analysis requires a more detailed mesh discretization in the affected region to avoid stress discontinuity. The FVDAM theory, on the contrary, satisfies the equilibrium equations in each subvolume at large, which allows the excellent convergency in the region of large stress concentration.
To make the small discontinuities in Figure 6b more visible, we took the transverse normal stress σ 22 field and enlarged the region of the large stress concentration area (see Figure 7a). Small discontinuities appear in the FEM's prediction, which are not observed in the FVDAM results.
Then, In order to get a clearer picture of the problem, we take 2 points from the transverse normal stress σ 22  It is worth noting that the FVDAM and FEM are two different methods that share certain similarities. The formulation of local/global stiffness matrices and the procedure of discretizing unit cells are similar. However, the local stiffness matrices in FVDAM satisfy the equilibrium equations in each subvolume at large. The closed form expressions between tractions and displacements are readily available, hence the continuity conditions of both displacements and tractions are applied directly. On the contrary, the FEM method is based on the minimization of global potential and satisfies the continuity conditions of the node displacements solely.

Cyclic Loading
We further investigated the cyclic behavior of the unidirectional composite materials under unidirectional strain loading. Figure 8 shows the applied transverse strainε 22 which varies between 1% and 2% as a function of the incremental number. To demonstrate the thermal effect, the temperature changes of ∆T = − 550 • C and ∆T = 0 • C were applied to the composite material system prior to the mechanical loading.  Figure 9 illustrates the validity of the FVDAM for predicting the cyclic response of the unidirectional composite materials. As shown in the figure, the FVDAM predictions are in good agreement with the FEM results. The ratcheting effect relies heavily on the plastic parameter µ. In the case of µ = 0 where the constitutive model degrades to the first version of the Ohno-Wang model, the ratcheting behavior is negligible. The maximum and the minimum stressesσ 22 at theε 22 = 1% and ε 22 = 2% strain in the cyclic curve stay on two horizontally parallel lines, indicating no ratcheting effect in this case. In the case of µ = 0.15, the ratcheting behavior becomes obvious. The peak and trough stresses decrease gradually with the increase of the cyclic number. Increasing the plastic parameter will further accentuate the differences in the composite transverse response due to the ratcheting effect. In the case of µ = 1 where the constitutive model degrades to the AF model, the ratcheting effect becomes the most significant.
The influence of the thermal residual stresses on the cyclic response of the present material system is negligible under the transverse loading case. This is consistent with expectations because a similar response was found in [49] and Pindera and Bansal [50] for the transverse response of a unidirectional metal matrix composite. The effect of residual stresses, which is pronounced under the off-axis loading, will be discussed in the next section.

Off-Axis Loading
We end this paper by illustrating the cyclic response under off-axis loading. In order to establish the stress-strain relations in the global coordinate systems, we rotate the principal coordinates (1, 2, 3) around the '3' axis by θ degrees. As shown in Figure 10a, axis '1' is parallel to the fiber, axis '2' is perpendicular to the fiber in-plane direction and axis '3' is perpendicular to the fiber out-of-plane direction. According to the coordinate transformation relations in the theory of elasticity, the stresses σ θ = [σ xx ,σ yy ,σ zz ,σ yz ,σ xz ,σ xy ] T in the global coordinates (x, y, z) are related to the stressesσ = [σ 11 ,σ 22 ,σ 33 ,σ 23 ,σ 13 ,σ 12 ] T in the principal coordinates (1, 2, 3) (see Figure 10b) through: Similarly, the strainε θ = [ε xx ,ε yy ,ε zz , 2ε yz , 2ε xz , 2ε xy ] T in the global coordinates (x, y, z) can be transformed to stressε = [ε 11 ,ε 22 ,ε 33 , 2ε 23 , 2ε 13 , 2ε 12 ] T in the principal coordinates (1,2,3): Therefore, the stress-strain relation in the global coordinates can be written as: whereC * = T −1 1 C * T 2 is the stiffness matrix in the global coordinates and (σ pl +σ th ) θ = T −1 1 (σ pl + σ th ) are the thermo-plastic stresses in the global coordinates. Once the applied strainε θ is obtained, the strain componentsε in the principal coordinate system can be calculated usingε = T 2εθ . The stress-strain relations in the global coordinates can be determined after solving Equations (13), (16) and (31). The cyclic stress-strain response of the composite materials with different off-axis angles (θ = 0 • , 15 • , 30 • , 45 • , 60 • and 75 • ) are illustrated in Figure 11. It is revealed that the cyclic response of the unidirectional composites is sensitive to the off-axis angle. The composite response is the stiffest in the case of 0 • since the loading direction is along the fiber direction. With the increase of the off-axis angle, the overall stresses in the global coordinate system decrease first then increase. The thermal residual stresses have a more pronounced effect under the small off-axis loading angle cases but to a lesser extent under the moderate and large off-axis loading angle cases. The ratcheting effects are more obvious under the off-axis loading cases due to the significant out-of-plane shear-coupling. Figure 11. The cyclic stress-strain response with consideration of the thermal effect for the composite materials with different fiber orientations.

Discussions
The existing micromechanics models for simulating the elastic-plastic response of the metal matrix composite materials are often under the framework of the Eshelby and the Mori-Tanaka theory. Most of these approaches, however, are not capable of simulating the distributions and evolutions of locals in the composites. Conventionally, the finite-element micromechanics model is the prevailing way and a golden standard for simulating both local fields and homogenized stress-strain responses of composite materials. Arbitrary microstructures can be easily addressed by the finite-element technique through the readily available commercial packages. The cyclic behavior, as well as other physical effects, can be simulated using the user programmable feature. A major shortcoming of the finite element theory is that the construction of local/global system of equations is under the framework of the variational principle. Hence, the simulation of composite materials often requires a detailed mesh near the interface due to the large modulus contrast between constituent materials. Compared with the FEM, the FVDAM is a semi-analytical method which provides an accuracy comparable to the FEM but with better efficiency. What is more, for the fatigue life prediction of the composite materials subjected to cyclic loading, since the proposed model can successfully calculate the local fields of the composites, the simulation of damage evolution and ultimately failure could be a worthwhile pursuit.

Summaries and Conclusions
The FVDAM model is extended by incorporating the cyclic plasticity model with the nonlinear kinematic hardening rule originally proposed by Abdel-Karim and Ohno, which enables efficient and accurate simulation of elastoplastic ratcheting behavior and local field distributions of the metal matrix composite materials. Comparisons of the results generated by the FVDAM and the finite-element method not only demonstrate the accuracy of the FVDAM theory but the effectiveness and efficiency. The generated results also indicate that the FVDAM has better convergence and stability than the FEM theory since the FVDAM eliminates the stress discontinuities in the vicinity of the interface which are observed in the FE simulations.
The FVDAM is used to investigate the influence of the thermal residual stress and fiber-orientations on the ratcheting behavior of metal matrix composites. The numerical simulations show that the thermal residual stresses provide a slight enhancement in homogenized stresses. However, the influence of the thermal residual stresses on the cyclic response of the present material system is negligible under the transverse loading case. Ratcheting effects are more obvious under the off-axis loading cases due to the significant out-of-plane shear-coupling. Last but not least, the thermal residual stresses have a more pronounced effect under the lower off-axis loading cases.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: FVDAM Finite-volume direct averaging micromechanics FE Finite element FEM Finite element method