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

Abstract

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.

1. 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, 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.

2. Results and Discussion

2.1. 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.

Table 1. Comparison of the present PBC and HBC models with some other theoretical solutions.

Models	E1 (GPa)	E3 (GPa)	G12 (GPa)	G23 (GPa)	v23
Present PBC model	392.0	391.0	164.9	165.6	0.179
Present HBC model	391.1	393.5	167.6	173.5	0.174
Mori-Tanaka’s method [6,25]	391.7	391.0	165.6	165.1	0.179
Self-consistent method [7,26]	391.6	391.0	165.5	165.2	0.180
Modified self-consistent method [8]	386.6	389.0	161.6	165.6	0.179

Table 1. Comparison of the present PBC and HBC models with some other theoretical solutions.

Models	E1 (GPa)	E3 (GPa)	G12 (GPa)	G23 (GPa)	v23
Present PBC model	392.0	391.0	164.9	165.6	0.179
Present HBC model	391.1	393.5	167.6	173.5	0.174
Mori-Tanaka’s method [6,25]	391.7	391.0	165.6	165.1	0.179
Self-consistent method [7,26]	391.6	391.0	165.5	165.2	0.180
Modified self-consistent method [8]	386.6	389.0	161.6	165.6	0.179

2.2. Influence of the Interfacial Properties on the Overall Elastic Properties

In the cohesive model, the interface penalty stiffness K_interface is defined as a function of the interface thickness, h_interface, and the elastic modulus of the interface, E_interface, i.e., $K_{\text{interface}}=E_{\text{interface}}/h_{\text{interface}}$ We introduce an interfacial stiffness dimensionless parameter $\tilde{k}=K_{\text{interface}}/\left[(E_m+E_f)/2h_{\text{interface}}\right]$ to represent relative modulus compared with the average of elastic moduli of matrix and fiber. In all simulations, we assume the interfacial thickness h_interface as one tenth of the carbon-fiber radius. $\tilde{k}=1$ when E_interface equals to $(E_m+E_f)/2$. K_nn = K_tt = K_ss = K_interface, and the other K_ij ($i\ne j$) are specified as 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 E₁, G₁₂ and G₂₃ 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 E₃ along the fiber direction and v₂₃ 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 y₁, tension along y₃, shear along y₁y₂, and shear along y₂y₃ directions in case of $\tilde{k}=1\times {10}^{-5}$ 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.

Figure 1. Influence of interface stiffness on the effective elastic constants (the arrows in red and in lightblue represent prefect interface and $\tilde{k}$ = 1, respectively).

Figure 1. Influence of interface stiffness on the effective elastic constants (the arrows in red and in lightblue represent prefect interface and $\tilde{k}$ = 1, respectively).

Figure 2. Stress influence functions obtained in the (a) tension along the y₁; (b) tension along the y₃; (c) shear along y₁y₂; and (d) shear along y₂y₃ linear perturbation steps in case of $\tilde{k}=1\times {10}^{-5}$.

Figure 2. Stress influence functions obtained in the (a) tension along the y₁; (b) tension along the y₃; (c) shear along y₁y₂; and (d) shear along y₂y₃ linear perturbation steps in case of $\tilde{k}=1\times {10}^{-5}$.

2.3. 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.

Figure 3. Mesh-sensitivity analysis in the nonlinear analyses: (a) coarse; (b) medium; and (c) fine mesh densities.

Figure 3. Mesh-sensitivity analysis in the nonlinear analyses: (a) coarse; (b) medium; and (c) fine mesh densities.

The three meshes are employed to simulate the macroscopic uniaxial tension tests along y₁ 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 ${\mathrm{\sigma }}_1^c-{\mathrm{\epsilon }}_x^c$ curves, whereas the coarse mesh overestimates the ${\mathrm{\sigma }}_1^c-{\mathrm{\epsilon }}_1^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.

Figure 4. Influence of the mesh size on the macroscopic response of C/SiC composites subject to uniaxial tension along y₁ direction.

Figure 4. Influence of the mesh size on the macroscopic response of C/SiC composites subject to uniaxial tension along y₁ direction.

2.4. 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 $\tilde{t}=t_{\text{interface}}/f_t$, where t_nn = t_ss = t_tt = t_interface. The interface fracture energy, G_interface, 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 G_c is equal to the area under the traction-separation curve, i.e., it must be larger than ${\left(t_{\text{interface}}\right)}^2/\left(2K_{\text{interface}}\right)$. So we introduce an interfacial critical fracture energy dimensionless parameter $\tilde{G}=G_{\text{interface}}/\left({t_{\text{interface}}}^2/K_{\text{interface}}\right)$ to represent the relative fracture energy. Since the t_interface and K_interface are varied in the simulations, for simplicity, we define the critical fracture energy dimensionless parameter as $\tilde{G}=G_{\text{interface}}/\left[{\left(0.1f_t\right)}^2/0.1\frac{(E_m+E_f)}{2h_{\text{interface}}}\right]$.

2.4.1. Uniaxial Transverse Tensile Strength along y₁ Direction

The uniaxial transverse tension along y₁ direction is simulated by using the present HBC models. The relations of the macroscopic stresses ${\mathrm{\sigma }}_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., $\tilde{t}=0.1$ and $\tilde{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 $\tilde{t}$ where $\tilde{k}$ = 0.1 and $\tilde{G}=100$, 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 $\tilde{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).

Figure 5. (a) Simulated macroscopic stresses with respect to the loading strain for the C/SiC composites with different interfacial stiffness subject to transverse tension along y_{1 s}direction; inset (b) cracking near the interface for a weak interfacial stiffness; (c) cracking near the corners for a perfectly bonded interface between the fiber and the ceramics; and (d) cracking near the interface and the corners almost at the same time. The black arrows point different loading strain, and the red arrow represents the descent direction of $\tilde{k}$.

Figure 5. (a) Simulated macroscopic stresses with respect to the loading strain for the C/SiC composites with different interfacial stiffness subject to transverse tension along y_{1 s}direction; inset (b) cracking near the interface for a weak interfacial stiffness; (c) cracking near the corners for a perfectly bonded interface between the fiber and the ceramics; and (d) cracking near the interface and the corners almost at the same time. The black arrows point different loading strain, and the red arrow represents the descent direction of $\tilde{k}$.

Figure 6. Simulated macroscopic stresses with respect to the loading strain for the C/SiC composites with different interfacial strengths subject to transverse tension along the y₁ direction. Both the black and red arrows represent the descent direction of $\tilde{t}$.

Figure 6. Simulated macroscopic stresses with respect to the loading strain for the C/SiC composites with different interfacial strengths subject to transverse tension along the y₁ direction. Both the black and red arrows represent the descent direction of $\tilde{t}$.

Figure 7. (a) Simulated macroscopic stresses with respect to the loading strain for the C/SiC composites with different interfacial fracture energy subject to transverse tension along y₁ direction; inset (b) debonding on the interface; and (c) cracking near the interface in the matrix. The black arrows point different loading strain, and the red arrow represents the descent direction of $\tilde{G}$.

Figure 7. (a) Simulated macroscopic stresses with respect to the loading strain for the C/SiC composites with different interfacial fracture energy subject to transverse tension along y₁ direction; inset (b) debonding on the interface; and (c) cracking near the interface in the matrix. The black arrows point different loading strain, and the red arrow represents the descent direction of $\tilde{G}$.

The ultimate transverse tensile strength is analytically given as [28]:

$${({{\mathrm{\sigma }}^T}_1)}_{\text{max}}=E_1{\left({\mathrm{\epsilon }}_1^T\right)}_{\text{max}}$$

(1)

The ultimate tensile strain of composites ${\left({\mathrm{\epsilon }}_1^T\right)}_{\text{max}}$ can be expressed in terms of ${\left({\mathrm{\epsilon }}_m^T\right)}_{\text{max}}$ as,

$${\left({\mathrm{\epsilon }}_1^T\right)}_{\text{max}}=\left[\frac{2r}{s}\frac{E_m}{E_f}+\left(1-\frac{2r}{s}\right)\right]{\left({\mathrm{\epsilon }}_m^T\right)}_{\text{max}}$$

(2)

where the overall transverse Young’s modulus is $E_1$ which can be found in Table 1, s represents the distance between centre of fibers, r is the radius of fibers, and ${\left({\mathrm{\epsilon }}_m^T\right)}_{\text{max}}$ 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.

Figure 8. Ultimate transverse tensile strength with respect to the fiber volume fraction.

Figure 8. Ultimate transverse tensile strength with respect to the fiber volume fraction.

2.4.2. Uniaxial Longitude Tensile Strength along y₃ Direction

The uniaxial longitude tension along y₃ direction is simulated by using the present HBC model. The relations of the macroscopic stresses ${\mathrm{\sigma }}_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.

The fiber tensile strength is expressed as $f_t^f=E_f{e}_t^f$, and the matrix tensile strength is expressed as $f_t^m=E_m{e}_t^m$. The composite tensile stress ${\mathrm{\sigma }}_3^c$ can be written as the function of the loading strain e as:

$${\mathrm{\sigma }}_3^c=\{\begin{array}{c}\left[E_f{V}_f+E_m\left(1-V_f\right)\right]ee\le e_t^m\\ V_f{E}_{f}e{e}_t^m<e\le e_t^f\end{array}$$

(3)

As the loading strain e increases until $e_t^m$, the composite tensile stress increases and then drops sharply because of the crack of the ceramic matrix when e equals $e_t^m$. The corresponding failure strength is $f_t^{cm}=\left[E_f{V}_f+E_m\left(1-V_f\right)\right]e_t^m$), as shown in Figure 9a,b. After the matrix fails, only the fibers of the composite are subjected to the loading, so the composite’s ultimate longitudinal tensile strength is $f_t^{cf}=V_f{E}_f{e}_t^f$, as shown in Figure 9a. From our theoretical and simulated curves shown in Figure 9 it can be seen that the matrix fractures first and then the carbon fiber fractures, because carbon fiber’s failure strain $e_t^f$ is greater than SiC failure strain $e_t^m$. The comparison shows a great agreement between the simulated results and the theoretical curve.

Figure 9. (a) Simulated macroscopic stresses in respect to the loading strain for the C/SiC composites with different interfacial properties subject to longitudinal tension along the y₃ direction; and (b) the enlarged view in the range of loading strain 0.0 to 0.002.

Figure 9. (a) Simulated macroscopic stresses in respect to the loading strain for the C/SiC composites with different interfacial properties subject to longitudinal tension along the y₃ direction; and (b) the enlarged view in the range of loading strain 0.0 to 0.002.

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:

$$f_t^c=\{\begin{array}{c}{E}_f\frac{f_t^m}{E_m}{V}_f+f_t^m\left(1-V_f\right)0\le V_f\le V_{f \text{min}}\\ V_f{f}_t^f{V}_{f \text{min}} <V_f\le 1\end{array}$$

(4)

where $V_{f \text{min}}$ 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.

Figure 10. Ultimate longitudinal tensile strengths with respect to the fiber volume fraction. The arrow represents the ultimate longitudinal tensile strength with respect to the fiber volume fraction calculated using Equation (4).

Figure 10. Ultimate longitudinal tensile strengths with respect to the fiber volume fraction. The arrow represents the ultimate longitudinal tensile strength with respect to the fiber volume fraction calculated using Equation (4).

3. Experimental Section

3.1. 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.

Figure 11. (a) Sketch of fiber-reinforced ceramic matrix composites; and (b) an enlarged 3D unit cell.

Figure 11. (a) Sketch of fiber-reinforced ceramic matrix composites; and (b) an enlarged 3D unit cell.

For coarse scale problem we have:

$${\overline{L}}_{ijmn}{\mathrm{\epsilon }}_{mn,x_j}^c+{\overline{b}}_i=0 \text{on} \Omega$$

(5)

$$u_i^c\left(x\right)={\overline{u}}_i \text{on}{\mathrm{\Gamma }}_u;{\overline{\mathrm{\sigma }}}_{ij}{n}_j={\overline{t}}_i \text{on} {\mathrm{\Gamma }}_t$$

(6)

where $\mathrm{\Omega }$ is the domain of the coarse scale problem; ${\mathrm{\Gamma }}_u$ and ${\mathrm{\Gamma }}_t$ the displacement and traction boundaries, respectively; x is the coarse scale position vectors, respectively; $\mathrm{\zeta }$ satisfies $0<\mathrm{\zeta }\le 1$; $u_i^c$ is the coarse scale displacement as a function of x; ${\mathrm{\epsilon }}_{mn}^c=\frac{1}{2}\left(\frac{\partial u_m^c}{\partial x_n}+\frac{\partial u_n^c}{\partial x_m}\right)$ is the coarse scale strain, and ${\overline{b}}_i$ is the average unit cell body force; ${\overline{u}}_i$, ${\overline{b}}_i$ and ${\overline{\mathrm{\sigma }}}_{ij}$ represent displacement, traction and stress boundary conditions, respectively. Summation convention is employed for repeated indices.

For unit cell problem,

$${\left[L_{ijkl}\left({\mathrm{\chi }}_{\left(k,y_l\right)mn}+I_{klmn}\right)\right]}_{,y}=0 \text{on} \mathrm{\Theta }$$

(7)

$${\mathrm{\chi }}_{imn}\left(y\right)={\mathrm{\chi }}_{imn}\left(y+Y\right)\text{on} \partial \mathrm{\Theta }; {\mathrm{\chi }}_{imn}\left(y\right)=0 \text{on} \partial {\mathrm{\Theta }}^{\text{vert}}$$

(8)

where $\mathrm{\Theta }$ is the domain of the unit cell; $\partial {\mathrm{\Theta }}^{\text{vert}}$ the vertices of the unit cell; $y=x/\mathrm{\zeta }$ the fine scale position vector; Y the period of the associated function; $I_{klmn}=\left({\mathrm{\delta }}_{mk}{\mathrm{\delta }}_{nl}+{\mathrm{\delta }}_{nk}{\mathrm{\delta }}_{ml}\right)$; ${\mathrm{\chi }}_{\left(k,y_l\right)mn}=\frac{1}{2}\left(\frac{\partial {\mathrm{\chi }}_{kmn}}{\partial y_l}+\frac{\partial {\mathrm{\chi }}_{lmn}}{\partial y_k}\right)$; and the homogenized constitutive tensor ${\overline{L}}_{ijmn}$ is given as:

$${\overline{L}}_{ijmn}=\frac{1}{\left|\mathrm{\Theta }\right|}{\displaystyle \underset{\mathrm{\Theta }}{\mathrm{\int }}{\mathrm{\sigma }}_{ij}^{mn}\left(y\right)}d\mathrm{\Theta }$$

(9)

where ${\mathrm{\sigma }}_{ij}^{mn}\left(y\right)$ is the stress influence function defined as:

$${\mathrm{\sigma }}_{ij}^{mn}\left(y\right)=L_{ijkl}\left({\mathrm{\chi }}_{\left(k,y_l\right)mn}+I_{klmn}\right)$$

(10)

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).

Figure 12. (a) An 3D RVE with periodic boundary conditions (PBCs); and (b) homogenization boundary conditions (HBCs)―tensile case along y₁ 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.

Figure 12. (a) An 3D RVE with periodic boundary conditions (PBCs); and (b) homogenization boundary conditions (HBCs)―tensile case along y₁ 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.

3.2. 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:

Fs = σ^––p tan φ–d = 0

(11)

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 ${\mathrm{\sigma }}_1,{\mathrm{\sigma }}_2$ and ${\mathrm{\sigma }}_3$:

$$p=-\frac{1}{3}\left({\mathrm{\sigma }}_1+{\mathrm{\sigma }}_2+{\mathrm{\sigma }}_3\right)$$

(12)

Here we denote $\overline{\mathrm{\sigma }}$ as the Mises equivalent stress:

$$\overline{\mathrm{\sigma }}=\sqrt{\frac{1}{2}\left[{\left({\mathrm{\sigma }}_1-{\mathrm{\sigma }}_2\right)}^2+{\left({\mathrm{\sigma }}_2-{\mathrm{\sigma }}_3\right)}^2+{\left({\mathrm{\sigma }}_3-{\mathrm{\sigma }}_1\right)}^2\right]}$$

(13)

We can determine the material’s angle of friction φ and its cohesion d through the uniaxial tensile strength $f_t$ and uniaxial compressive strength $f_c$:

$$\text{tan} \mathrm{\phi }=\frac{3\left(f_c-f_t\right)}{f_c+f_d}$$

(14)

$$d=\frac{2f_t{f}_c}{f_t+f_c}$$

(15)

The cohesion d is equal to yield stress in the case of f_c = f_t, i.e., no difference between compressive and tensile strengths. For SiC ceramics, the maximum tensile strength f_t and maximum compressive strength f_c 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.

Figure 13. (a) Yield surfaces in the p-$\overline{\text{σ}}$ plane in the Drucker-Prager model; and (b) uniaxial stress-strain curve.

Figure 13. (a) Yield surfaces in the p-$\overline{\text{σ}}$ plane in the Drucker-Prager model; and (b) uniaxial stress-strain curve.

Table 2. Properties of the component materials of C/SiC composites.

Material Properties	SiC	Carbon-Fiber
Young’s modulus (GPa)	Em 400	Ef 350
Poisson’s ratio	vm 0.14	vf 0.3
Tensile strength (MPa)	$f_t^m$ 310	$f_t^f$ 1380
Compressive strength (MPa)	$f_c^m$ 3900	-
Volume fraction	Vm 0.78	$V_f$ 0.22

Table 2. Properties of the component materials of C/SiC composites.

Material Properties	SiC	Carbon-Fiber
Young’s modulus (GPa)	Em 400	Ef 350
Poisson’s ratio	vm 0.14	vf 0.3
Tensile strength (MPa)	$f_t^m$ 310	$f_t^f$ 1380
Compressive strength (MPa)	$f_c^m$ 3900	-
Volume fraction	Vm 0.78	$V_f$ 0.22

3.3. 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: t_n, t_s and t_t, 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 T₀ as the original thickness of the cohesive element, the nominal strains can be defined as:

$${\mathrm{\epsilon }}_n=\frac{{\mathrm{\delta }}_n}{T_0},{\mathrm{\epsilon }}_s=\frac{{\mathrm{\delta }}_s}{T_0},{\mathrm{\epsilon }}_t=\frac{{\mathrm{\delta }}_t}{T_0}$$

(16)

The elastic behavior can then be written as:

$$t=\left\{\begin{array}{l}{t}_n\\ t_s\\ t_t\end{array}\right\}=\left[\begin{array}{ccc}{K}_{nn}& K_{ns}& K_{nt}\\ K_{sn}& K_{ss}& K_{st}\\ K_{tn}& K_{ts}& K_{tt}\end{array}\right]\left\{\begin{array}{l}{\mathrm{\delta }}_n\\ {\mathrm{\delta }}_s\\ {\mathrm{\delta }}_t\end{array}\right\}$$

(17)

Damage of the traction-separation response for cohesive surface is defined within the same general framework used for conventional materials [32]. A quadratic stress damage initiation criterion and an energy damage evolution law are defined for modeling the debonding of the interfaces, as shown in Figure 14b. The quadratic stress criterion suggests that damage initiates when a quadratic interaction function involving the contact stress ratios reaches the value of one. This criterion can be represented as:

$${\left\{\frac{\langle t_n\rangle }{t_n^0}\right\}}^2+{\left\{\frac{t_s}{t_s^0}\right\}}^2+{\left\{\frac{t_t}{t_t^0}\right\}}^2=1$$

(18)

where t_n, t_s, and t_t represent the contact stress normal to the interface, along the first and the second shear directions, respectively $t_n^0$, $t_s^0$ and $t_t^0$ represent the peak values of the contact stress when the separation is either purely normal to the interface or purely in the first or the second shear direction, respectively. The symbol $\langle \rangle$ represents the Macaulay bracket with the usual interpretation, indicating that a purely compressive stress state does not initiate damage. In the damage evolution criterion, 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 $G_c$ is equal to the area under the traction-separation curve shown in Figure 14b.

Figure 14. (a) Imperfect interface defined by surface-based cohesive behavior; and (b) typical traction-separation response.

Figure 14. (a) Imperfect interface defined by surface-based cohesive behavior; and (b) typical traction-separation response.

4. 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 E₁, G₁₂ and G₂₃ decrease, but E₃ and v₂₃ 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.