Homogenized Finite Element Analysis on Effective Elastoplastic Mechanical Behaviors of Composite with Imperfect Interfaces

A three-dimensional (3D) representative volume element (RVE) model was developed for analyzing effective mechanical behavior of fiber-reinforced ceramic matrix composites with imperfect interfaces. In the model, the fiber is assumed to be perfectly elastic until its tensile strength, and the ceramic material is modeled by an elasto-plastic Drucker-Prager constitutive law. The RVE model is then used to study the elastic properties and the tensile strength of composites with imperfect interfaces and validated through experiments. The imperfect interfaces between the fiber and the matrix are taken into account by introducing some cohesive contact surfaces. The influences of the interface on the elastic constants and the tensile strengths are examined through these interface models.


Introduction
As an important type of ceramic matrix composites, fiber-reinforced ceramics (FRCs) such as carbon-fiber/silicon-carbide (C/SiC) are becoming popular and important due to their unique thermal, OPEN ACCESS mechanical and chemical stability in various environments, high strength and excellent thermal shock resistance of ceramics, and high toughness of carbon fibers at elevated temperature [1]. The assessment of mechanical properties of such composites is, however, much more complex than that of conventional ceramics, as the composites may be partly or highly anisotropic. Usually, the physical and mechanical properties of FRCs depend on properties of their constituents and the corresponding geometry and concentration (e.g., volume fraction of fibers, fiber/matrix interphase structure, fiber weave architecture, and matrix properties). It is noted that when the fibers are embedded into the ceramic matrix to form composites, the matrix bonds fibers together and transfers loads to the fibers through the interfaces between them. Thus, the fiber/matrix interfaces govern to some extent the transverse tensile strength and the fracture behavior of the composite [2]. Experimental studies [3,4] revealed that interfacial properties also play an important role in affecting the macroscopically effective properties of FRCs.
Existing schemes for predicting macroscopically effective properties of composites include Mori-Tanaka method [5,6], self-consistent method [7,8], generalized self-consistent method [9], combination of the Mori-Tanaka method and the iso-stress or iso-strain assumptions [10], Christensen's approach [11], and various mathematical homogenization methods [12,13]. Many works, for example [14][15][16], have been done to study the effects of interfacial properties on effective properties of composites, but the components of composites were usually assumed to be elastic for simplicity in the most existing theoretical models. Ju and Yanase [17] proposed an elasto-plastic damage formulation to predict the overall transverse mechanical behavior of continuous fiber reinforced ductile matrix composites with the framework of micromechanics and homogenization by incorporating the interfacial damage. Alternatively, through investigating interphase effect on elastic and thermal conductivity response of polymer composite materials, Mortazavi et al. [18,19] compared the capability of the Mori-Tanaka method and the three-dimensional (3D) finite element (FE) analysis and concluded that despite complexities for modeling of high volume concentrations and aspect ratios for fillers, FE simulations are more reliable and promising than the other schemes. Based on FE method, Taliercio and Coruzzi [20] estimated in-plane transverse strengths using a representative volume element (RVE), in which the perfect bonding is assumed. Yang and Qin [21,22] investigated effective elastic-plastic properties of fiber-reinforced composites. Caporale et al. [23] implemented an interfacial failure model by connecting the fibers and the matrix at the finite element nodes by normal and tangential brittle-elastic springs, in which the matrix and fibers are considered homogeneous, isotropic and linearly elastic. Rahul-Kumar et al. [24] concluded that the cohesive element can be used to describe the polymer interfacial fracture. These works did not, however, couple the brittle material constitutive law and interfacial debonding in the approaches mentioned above. In addition, those models are not easy to be realized in practical analysis on the effect of interfacial properties on the macroscopically effective elastoplastic properties of composites.
The purpose of this study is to develop a 3D RVE model based on a unidirectional, long-fiber-reinforced ceramic matrix composites, using the computational homogenization FE method which can handle imperfect interface between the fiber and the matrix. Then, the model is incorporated into the commercial FE software ABAQUS through a user subroutine interface. In the RVE, the fiber is assumed to be linear elastic before the stress reaches its tensile strength and the ceramic material is modeled by an elasto-plastic Drucker-Prager constitutive law. The imperfect interfaces between fiber and matrix are taken into account by introducing some cohesive contact surfaces. Making use of the proposed model, comprehensive analyses on the influence of interfacial properties on the macroscopically effective elasto-plastic properties of composites, including the macroscopic stiffness and strength are conducted.

Model Validation
The reliability of both the present periodic boundary condition (PBC) and homogeneous boundary conditions (HPC) models is first assessed in estimating the effective elastic constants of FRCs by comparing them with theoretical results. The macroscopic elastic constants of the composites obtained using the present PBCs and HBCs models are depicted in Table 1. For comparison, the overall properties estimated using the Mori-Tanaka method [6,25], the self-consistent method [7,26] and the modified self-consistent method [8], are also calculated here and listed in Table 1. It can be seem from Table 1 that results from the present model show a good agreement with the theoretical results. zero. The damage initiation criterion and evolution law are not defined in this subsection because they influence the macroscopic elastic properties slightly in the initial small elastic deformation stage. Figure 1 plots the change trends of the macroscopic elastic constants with respect to the interfacial stiffness. It can be seen from Figure 1 that the overall elastic constants E1, G12 and G23 decrease gradually as the interfacial stiffness decreases, but they change slightly when the interface is stronger than the average value of the fiber and the matrix. The longitudinal Young's modulus E3 along the fiber direction and v23 are independent of the interfacial stiffness. We list the contour plots of the stress influence functions obtained in the linear perturbation steps including tension along y1, tension along y3, shear along y1y2, and shear along y2y3 directions in case of 5 1 10 k − = ×  in Figure 2, from which it can be found that the bonding of the fiber and the matrix in the composites, subject to the longitudinal tension, will not be affected by the low interfacial stiffness.

Mesh-Sensitive Analysis and Model Validation in Estimating the Ultimate Tensile Strength
Before investigating the ultimate tensile strength, a nonlinear numerical study is first performed to study the sensitivity of the predicted macroscopic responses of the considered FRCs to mesh refinement. Three mesh densities for the unit cell are considered, namely a "coarse", a "medium" and a "fine" mesh, as shown in Figure 3. The three meshes are employed to simulate the macroscopic uniaxial tension tests along y1 axis. The macroscopic elastoplastic responses of FRCs with perfect interfaces are considered here. Figure 4 shows the macroscopic stress-strain curves obtained with the three meshes. The medium and fine meshes predict the same tensile macroscopic 1 σ ε c c x − curves, whereas the coarse mesh overestimates the 1 1 σ ε c c − curve. Thus, no improvement seems to come from the use of a mesh finer than the medium one. The medium mesh will be used in all subsequent simulations. Figure 4 also indicates a good agreement between the FE predictions and the experimental results measured by Heredia et al. [27]. They measured the transverse ultimate tensile strength of C/SiC composites with 22% carbon fibers as 320 ± 30 MPa.

Influence of the Interfacial Properties on the Macroscopic Strength
As mentioned below, the HBCs are less time consuming than PBCs, hence they are more suitable for sufficiently large RVEs or nonlinear analyses. In this subsection, HBCs are chosen to predict the macroscopically ultimate strengths of the composites in order to save the computational time. If we assume the fiber volume fraction to be constant, composite structures may vary their stiffness and strength due to damage accumulation such as matrix cracking and fiber breakage during the loading process of the composite members.
We define the initial damage traction dimensionless parameter The interface fracture energy, Ginterface, as an interfacial property, is selected to define the evolution of debonding in terms of the energy required for failure after the initiation of debonding. As depicted in Section 3.3, the fracture energy Gc is equal to the area under the traction-separation curve, i.e., it must be larger than ( ) ( ) to represent the relative fracture energy. Since the tinterface and Kinterface are varied in the simulations, for simplicity, we define the critical fracture energy dimensionless parameter as

Uniaxial Transverse Tensile Strength along y1 Direction
The uniaxial transverse tension along y1 direction is simulated by using the present HBC models. The relations of the macroscopic stresses 1 σ c with the loading strain for the C/SiC composites with different interfacial stiffness are plotted in Figure 5a, where the critical interfacial damage strength and the critical interfacial fracture energy are assumed to be constant, i.e., 0.1 t =  and G  = 100. The singularity associated with the FE modeling is inevitable at the interface between the fiber and the matrix or the boundaries of the RVE. It must be noted that the Drucker-Prager model is a "smeared crack model", since it does not describe a single crack, but rather associates to any integration point with degraded mechanical properties. So the singularity affects slightly the overall response of the composites. It can be seen from Figure 5a that the interfacial stiffness plays an important role on the uniaxial transverse tensile strength. If we assume that the thickness of the interface is approximately one tenth of the fiber radius, the interface effect can be ignored when the interfacial modulus is stronger than the average of the Young's moduli of the fiber and the matrix. As the interfacial stiffness decreases, the transverse tensile strength decreases significantly. The FE results show a difference in the damage onset in the matrix of the composites with different interfaces. For a perfect interface, the present model predicts the initiation of the damage in four regions near the corner of the RVE, as shown in Figure 5b. For a strong interface, the damage commences near both the corner and the interface, as shown in Figure 5c, while for a weak interface, it commences only near the interface as shown in Figure 5d.
The influence of the interfacial strength on the ultimate transverse tensile strength is shown in Figure 6 for different t  where k  = 0.1 and 100 G =  , It can be seen that the interfacial strength is not the major factor in determining the ultimate transverse tensile strength of the fiber-reinforced ceramic matrix composites. In the calculation, interfaces with different interfacial fracture energy are also considered. Figure 7a shows that the ultimate transverse tensile strength is very sensitive to the interfacial fracture energy when G  is lower than 50. We find that the composite will damage first due to the debonding of the interfaces (seen in Figure 7b) and then begin cracking near the interface in the ceramic matrix (seen in Figure 7c).
where the overall transverse Young's modulus is 1 E which can be found in Table 1, s represents the distance between centre of fibers, r is the radius of fibers, and ( ) max ε T m is the tensile failure strain of matrix. Figure 8 shows that the ultimate transverse tensile strength of C/SiC composites with perfect or imperfect interfaces is insensitive to the fiber volume fraction, and the imperfect interface may reduce the strength enormously. When the fiber volume fraction is very low (1.2% in our simulations), it can be seen from Figure 8 that those three values converge to one certain value (i.e., the tensile strength of the matrix), implying that the interface effect can be ignored only when the fiber volume fraction is very low.

Uniaxial Longitude Tensile Strength along y3 Direction
The uniaxial longitude tension along y3 direction is simulated by using the present HBC model. The relations of the macroscopic stresses 3 σ c with the loading strain for the C/SiC composites, considering imperfect interfaces with different interfacial properties, are plotted in Figure 9. It can be found that the longitudinal tensile strength of the C/SiC composites is almost independent on the interfacial properties. The composites with a weaker interface have a bit higher longitudinal tensile strength because the weaker interface inhibits the interaction of the brittle matrix and the fiber. ( ) As the loading strain e increases until m t e , the composite tensile stress increases and then drops sharply because of the crack of the ceramic matrix when e equals m t e . The corresponding failure strength is  For C/SiC composites having fiber failure strain greater than matrix failure strain, the variation of composite longitudinal tensile strength with fiber volume fraction is governed by: where min f V is the minimum fiber volume fraction below which by adding the fibers to the matrix, the C/SiC will have lower ultimate longitude tensile strength than the matrix. Figure 10 shows that as the fiber volume fraction increases, the ultimate longitude tensile strength of C/SiC composites increases sharply, which satisfies the mixture rule.

Homogenized 3D RVE for FRCs
In the assumption of the two-scale asymptotic homogenization method, unidirectional and long-fiber composites are simplified as composites constructed by periodically and uniformly distributed unit cells, as shown in Figure 11a. An enlarged unit cell (also called RVE) is shown in Figure 11b. Making use of mathematical homogenization, a linear elastic static problem with periodic conditions could be decomposed into uncoupled fine and coarse scale problems. The macroscopically effective mechanical properties of the FRCs could be determined through the estimation of one unit cell at fine scale using the perturbation technology. The computational homogenization approach for linear and nonlinear solid problems, presented in [29] is used in our analysis. where where Θ is the domain of the unit cell; ; and the homogenized constitutive tensor ijmn L is given as: The two-scale algorithm described here can be generalized to account for material and geometric nonlinearities.
It is noted that boundary conditions can significantly affect the macro behavior of the RVEs during the homogenization simulation process. To study the effect, two types of boundary conditions are generally used. If an RVE with 3D periodic boundary conditions (PBCs) is used, the simulation results represent a macro structure consisting of periodically repeated cells. While choosing 3D homogeneous boundary conditions (HBCs), the simulation results would consider the RVE as the macro structure itself with its micro-constituents. HBCs are less time-consuming in computation than periodic boundary conditions, hence they are more suitable for sufficiently large RVEs. In the simulation, both PBCs and HBCs are chosen to estimate macroscopically effective elastic parameters of the FRCs (as shown in Figure 12a), while HBCs are chosen to predict macroscopically ultimate strengths in order to save computational time (as shown in Figure 12b). homogenization boundary conditions (HBCs)-tensile case along y1 direction (only the normal directions are fixed at the boundaries). The blue arrow represents the tension direction, and the red frame represents the configuration after tension.

Constitutive Model for the Components
Unidirectional long fiber-reinforced ceramic matrix composites are investigated in the present study. The fibers are supposed to be linear elastic when the stress level is below its tensile strength. The brittle behavior of the ceramics is described using the Drucker-Prager yield criterion which has already been implemented into the commercial FE software ABAQUS [30]. This model was originally developed for plain concrete subjected to multiaxial stresses and has been successfully used to estimate the transverse strengths for the ceramic matrix composites [20].
The Drucker-Prager failure surface is given by: where φ is the material's angle of friction and d is its cohesion (see Figure 13a). The equivalent compressive stress p is expressed as a function of the principal stresses 1 2 σ , σ and 3 σ : Here we denote σ as the Mises equivalent stress: We can determine the material's angle of friction φ and its cohesion d through the uniaxial tensile strength t f and uniaxial compressive strength c f : The cohesion d is equal to yield stress in the case of fc = ft, i.e., no difference between compressive and tensile strengths. For SiC ceramics, the maximum tensile strength ft and maximum compressive strength fc are 310 and 3900 MPa, respectively. Making use of Equations (14) and (15), the calculated material's angle of friction φ = 68.6°, dilation angle θ = φ = 68.6° which satisfies the associate flow. A sharp post-peak drop in strength is defined for approaching the behavior of a perfectly brittle ceramic material, as shown in Figure 13b. The post-peak strain softening behavior of the ceramics is inputted in the Drucker-Prager model, and can be simulated by means of the modified Riks method in ABAQUS. For carbon fiber, perfectly elastic behavior and a tensile strength of 1390 MPa is assumed. If the stress level exceeds the maximum strength for matrix and/or fiber, the Young's modulus E is degraded to 1% from its initial value at a particular integration point, while the shear modulus G is reduced to 20% of the initial value under the assumption that some shear stiffness remains due to the friction still present on the failure plane [31], which is realized through the user subroutine USDFLD in ABAQUS. The material properties used in the simulations are listed in Table 2. The superscripts m and f appearing in Table 1 and afterwards represent the variables associated with the matrix and fiber, respectively.

Cohesive Interfacial Model for Fiber/Matrix Interfaces
The cohesive elements employ failure criteria that combine aspects of strength-based analysis to predict the onset of the softening process at the interface and fracture mechanics to predict debonding propagation. If the interface thickness is negligibly small, it can be straightforward to define the surface-based cohesive response of the cohesive layer directly in terms of traction versus separation (see Figure 14a), which will spend less computational time compared with the cohesive element in ABAQUS. The available traction-separation model in ABAQUS assumes initially linear elastic behavior followed by the initiation and evolution of damage.
The nominal traction stress vector, t, consists of three-components: tn, ts and tt, which represent the normal (along the local 3-direction) and the two shear tractions (along the local 1-and 2-directions), respectively. The corresponding separations are denoted by δn, δs and δt. Denoting T0 as the original thickness of the cohesive element, the nominal strains can be defined as: damage evolution can be defined based on the energy that is dissipated as a result of the damage process, also called the fracture energy. The fracture energy c G is equal to the area under the traction-separation curve shown in Figure 14b.

Conclusions
Based on the homogenization method with periodic or homogenization boundary conditions, a 3D RVE model is developed. The proposed model has been validated with the theoretical results for the composites with perfect bonding between the fiber and the matrix. From the study, we found that for composites with imperfect interfaces, as the interface stiffness decreases, the E1, G12 and G23 decrease, but E3 and v23 remain almost as constants.
The obtained numerical results for the transverse tensile strength of the composites with perfect bonding agree well with those from experiments.
The imperfect interfaces between the fiber and the matrix are taken into account by introducing cohesive contact surfaces. The influences of the interface on the elastic constants and the tensile strengths are examined using the interfacial model. It is found that the imperfect interface can induce different damage onset in the matrix of the composites with different interfacial properties. In contrast the interfacial strength, the interfacial stiffness and fracture energy can significantly influence the transverse tensile strength of the composites, while the longitudinal tensile strength of the composites is almost independent of the interfacial properties, and it increases sharply with an increase in the fiber volume fraction, satisfying the mixture rule. The results indicate that the proposed approach is simple to use and efficient for performing realistic numerical analyses on complex 3-D fiber-reinforced composites.