Micro-Scale Numerical Simulation for Residual Strength of CFRP After Cyclic Tensile or Out-of-Plane Shear Loadings Fatigue

1. Introduction

Carbon fiber reinforced polymers (CFRPs) are widely used in the aerospace industry due to their excellent mechanical properties. Understanding their durability under long-term use is crucial; however, the current research predominantly relies on experimental approaches, with few studies employing numerical methods to investigate the strength degradation behavior of CFRPs deeply [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18].

Recently, entropy-based fracture criteria have been adopted for matrix resins, defining material failure as the critical value of entropy generation [19,20,21,22,23,24,25,26,27,28,29,30,31]. When a material is unused and undamaged, its molecular arrangement is orderly, resulting in low entropy. Upon loading, the orderly molecular structure partially collapses, leading to an increase in entropy. This collapse is equivalent to the material's inelastic deformation. Thus, with repeated loading, the molecular structure collapses, increasing entropy until it reaches a critical value, causing material failure. The amount of entropy generated at failure is inherent to the material, independent of geometric shape, stress level, and frequency. This allows for the assessment of irreversible cumulative damage, using entropy generation to clarify the degradation of properties and estimate the fatigue life and residual strength of structural materials under periodic loading. Deng et al. performed micro-scale simulations of repeated loading [32]. They used an entropy-based fracture criterion to estimate the fatigue life of CFRPs under multi-amplitude cyclic loading, developing a micro finite element model that separately considered the matrix resin and fibers. However, they focused only on the onset of failure and did not estimate the overall fatigue life of CFRP, failing to link entropy generation to the overall behavior of CFRP. Koyanagi et al. incorporated entropy into Hashin’s criteria, enabling the representation of strength degradation that traditional Hashin’s criteria could not, and conducted mesoscale simulations of cyclic loading [33]. They reproduced the occurrence of transverse cracks, showing stress level and frequency dependence on the number of cracks. However, the coefficient α linking entropy increase to strength degradation remained undetermined, preventing a correct association between entropy generation and strength degradation. If this coefficient α can be determined through micro-scale analysis, it would enable the integration of micro-scale and mesoscale analyses.

Additionally, composite materials like CFRPs have complex microstructures and anisotropic properties. To understand the material properties of the entire product, it is common to perform material tests as mentioned earlier, but the properties obtained from material tests reflect macro-structural characteristics, not micro-structural ones. Therefore, if multi-scale analysis, considering both micro-scale and meso-scale properties, becomes possible, material design can fully leverage the inherent properties of the materials.

In this study, we conducted creep tests on epoxy resin, which is widely used as the matrix resin for CFRP. Using a generalized Maxwell model with 15 elements as the viscoelastic model, we determined the material constants and the necessary constants for numerical simulations. To establish the relationship between entropy generation and the residual strength after cyclic loading, we performed analyses using Abaqus/Standard 2020, with the epoxy resin matrix behavior implemented via the Abaqus/Standard user subroutine UMAT. We carried out micro-scale repeated tensile and out-of-plane shear analyses incorporating an entropy-based damage criterion for the resin part of CFRP. As a result, we established a linear relationship between the strength of CFRP after fatigue loading and the overall entropy of CFRP.

2. Degradable Hashin’s Criteria

As previously mentioned, Koyanagi et al. introduced entropy into the Hashin’s failure criteria. Eq. 1 to 4 represent the failure criteria, Eq. 5 to 11 represent strength reduction, and Eq. 12 to 15 represent fracture energy reduction. $X_{\mathrm{T}}$，$X_{\mathrm{C}}$，$Y_{\mathrm{T}}$，$Y_{\mathrm{C}}$，$S_{12}$ (= $S_{13}$), and $S_{23}$ are the tensile strength in the longitudinal direction, compressive strength in the longitudinal direction, tensile strength in the transverse direction, compressive strength in the transverse direction, longitudinal (in-plane) shear strength, and transverse shear strength, respectively. $G_{\mathrm{f}\mathrm{t}},G_{\mathrm{f}\mathrm{c}},G_{\mathrm{m}\mathrm{t}}$, and $G_{\mathrm{m}\mathrm{c}}$ are the fracture toughness in longitudinal tension, longitudinal compression, transverse tension, and transverse compression, respectively. The four damage criteria – fiber tensile mode $\left(e_{\mathrm{f}\mathrm{t}}\right)$, fiber compressive mode $\left(e_{\mathrm{f}\mathrm{c}}\right)$, transverse tensile mode $\left(e_{\mathrm{m}\mathrm{t}}\right)$, and transverse compressive mode $\left(e_{\mathrm{m}\mathrm{c}}\right)$ – are calculated using the degraded strengths. When any of these criteria are met (e≥1), it is judged that damage has occurred. The constants ${\alpha }_{\mathrm{A}\mathrm{T}},{\alpha }_{\mathrm{A}\mathrm{C}}$, and ${\alpha }_o$ used in strength reduction and fracture energy reduction relate entropy $s$ to strength reduction and fracture energy reduction. Currently, these $\alpha$ values are arbitrary constants as they have not yet been determined.

Therefore, We aim to apply the entropy damage criteria to the matrix resin of CFRPs and conduct numerical analyses at the microscale. By correlating the residual strength of CFRPs after fatigue loading with entropy generation and determining $\alpha$ as a function of entropy, I intend to represent mesoscale strength degradation considering the results of microscale numerical simulations. In this study, tensile fatigue analysis $\left(Y_{\mathrm{T}}\right)$ and out-of-plane shear fatigue analysis $\left(S_{12}\right)$ were performed, and ${\alpha }_o\left(s\right)$ was determined.

$$e_{\mathrm{f}\mathrm{t}}={\left({\displaystyle \frac{{\sigma }_{11}}{X_{\mathrm{T}},n_0}}\right)}^2+{\left({\displaystyle \frac{{\tau }_{12}}{S_{12},n_0}}\right)}^2+{\left({\displaystyle \frac{{\tau }_{13}}{S_{13},n_0}}\right)}^2$$

(1)

$$e_{\mathrm{f}\mathrm{c}}={\left({\displaystyle \frac{{\sigma }_{11}}{X_{c,n_0}}}\right)}^2$$

(2)

$$e_{\mathrm{m}\mathrm{t}}={\left({\displaystyle \frac{{\sigma }_{22}+{\sigma }_{33}}{Y_{\mathrm{T},n_0}}}\right)}^2+{\displaystyle \frac{{\tau }_{23}^2-{\sigma }_{22}{\sigma }_{33}}{{\left(S_{23,n_0}\right)}^2}}+{\displaystyle \frac{{\tau }_{12}^2}{{\left(S_{12,n_0}\right)}^2}}+{\displaystyle \frac{{\tau }_{13}^2}{{\left(S_{13,n_0}\right)}^2}}$$

(3)

$$e_{\mathrm{m}\mathrm{c}}=\left[{\left({\displaystyle \frac{Y_{\mathrm{C},n_0}}{2S_{23,n_0}}}\right)}^2-1\right]{\displaystyle \frac{{\sigma }_{22}+{\sigma }_{33}}{Y_{\mathrm{C},n_0}}}+{\left({\displaystyle \frac{{\sigma }_{22}+{\sigma }_{33}}{2S_{23,n_0}}}\right)}^2+{\displaystyle \frac{{\tau }_{23}^2-{\sigma }_{22}{\sigma }_{33}}{{\left(S_{23,n_0}\right)}^2}}+{\displaystyle \frac{{\tau }_{12}^2}{{\left(S_{12,n_0}\right)}^2}}+{\displaystyle \frac{{\tau }_{13}^2}{{\left(S_{13,n_0}\right)}^2}}$$

(4)

$$X_{\mathrm{T},n_0}=\left(1-{\alpha }_{\mathrm{A}\mathrm{T}}{s}_{n_0}\right)X_{\mathrm{T},0}$$

(5)

$$X_{\mathrm{C},n_0}=\left(1-{\alpha }_{\mathrm{A}\mathrm{C}}{s}_{n0}\right)X_{\mathrm{C},0}$$

(6)

$$Y_{\mathrm{T},n_0}=\left(1-{\alpha }_{\mathrm{O}}{s}_{n0}\right)Y_{\mathrm{T},0}$$

(7)

$$Y_{\mathrm{C},n_0}=\left(1-{\alpha }_{\mathrm{O}}{s}_{n0}\right)Y_{\mathrm{C},0}$$

(8)

$$S_{12,n_0}=\left(1-{\alpha }_{\mathrm{O}}{s}_{n0}\right)S_{12,0}$$

(9)

$$S_{13,n_0}=\left(1-{\alpha }_{\mathrm{O}}{s}_{n0}\right)S_{13,0}$$

(10)

$$S_{13,n_0}=\left(1-{\alpha }_{\mathrm{O}}{s}_{n0}\right)S_{13,0}$$

(11)

$$G_{\mathrm{f}\mathrm{t},n_0}={\left(1-{\alpha }_{\mathrm{A}\mathrm{T}}{s}_{n0}\right)}^2{G}_{\mathrm{f}\mathrm{t},0}$$

(12)

$$G_{\mathrm{f}\mathrm{c},n_0}={\left(1-{\alpha }_{\mathrm{A}\mathrm{C}}{s}_{n0}\right)}^2{G}_{\mathrm{f}\mathrm{c},0}$$

(13)

$$G_{\mathrm{m}\mathrm{t},n_0}={\left(1-{\alpha }_{\mathrm{O}}{s}_{n0}\right)}^2{G}_{\mathrm{m}\mathrm{t},0}$$

(14)

$$G_{\mathrm{m}\mathrm{c},n_0}={\left(1-{\alpha }_{\mathrm{O}}{s}_{n0}\right)}^2{G}_{\mathrm{m}\mathrm{c},0}$$

(15)

3. Determination of Material Constants

The most commonly used resin in CFRP is epoxy resin; therefore, this study employed epoxy resin [23]. DGBEA 50 g and 4,4’-DDS 16.3 g were mixed to create test specimens with dimensions of 60 mm in length, 10 mm in width, and 3 mm in thickness. Creep tests were conducted under stresses of 10 MPa, 20 MPa, 30 MPa, and 40 MPa, at temperatures of 40°C, 50°C, 60°C, 70°C, and 80°C for each stress level, with each test lasting 1 hour. The performed creep test data were shifted using the Arrhenius-type temperature-time equivalence principle (Eq. 16). The parameters were set as $R=0.8314\mathrm{k}\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}·\mathrm{K},∆H_a=150\mathrm{k}\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}·\mathrm{K},T_0=300\mathrm{K}$.

$$\log {\alpha }_T=\beta {\displaystyle \frac{∆H_a}{R}}\left({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{T_0}}\right)$$

(16)

Using creep compliance, a time-strain function was created and fitted to the creep test results. The creep compliance and strain function are shown in Eq. 17 and Eq. 18, respectively. The constants were determined as $n=0.2,T_c=8.0\times {10}^7,E_0=3300\left[\mathrm{M}\mathrm{P}\mathrm{a}\right]$. The results above are shown in Figure 1. The creep test results shifted using the temperature-time equivalence principle are indicated by solid lines, and the strain function obtained through fitting is indicated by dashed lines.

$$C\left(t\right)={\displaystyle \frac{1+{\left(\frac{t}{T_c}\right)}^n}{E_0}}$$

(17)

$$\epsilon \left(t\right)=C\left(t\right)\times {\sigma }_0$$

(18)

The stress relaxation curve $E\left(t\right)$, which is necessary for determining the stiffness and viscosity coefficients of the matrix resin, was obtained by taking the inverse of $C\left(t\right)$, resulting in a master curve. Although this conversion requires a Laplace transform, numerical substitution and calculation showed that the values obtained by taking the inverse of $C\left(t\right)$ were almost identical. Therefore, in this study, which deals with numerical calculations, they are considered equivalent. $E\left(t\right)$ is shown in Eq. 19

$$E\left(t\right)={\displaystyle \frac{1}{C\left(t\right)}}={\displaystyle \frac{E_0}{1+{\left(\frac{t}{T_c}\right)}^n}}$$

(19)

Next, for the 15 Maxwell models depicted in Figure 4, the stiffness of each element was set to 220 MPa, and the viscosity coefficients were fitted (Figure 2). A strain of 0.1 was assumed. The equation used is shown in Eq. 20. The determined stiffness and viscosity coefficients of each Maxwell element are presented in Table 1.

$${\epsilon }_d={\displaystyle \frac{{\sigma }_{jold}}{{\eta }_j∆t}}, {\sigma }_j=\left(0.1-{\epsilon }_d\right)E_j, \sigma \left(t\right)=\sum _{j=1}^{15}{\sigma }_j$$

(20)

The determination of the non-linear coefficient, as shown in Eq. 22, was conducted. Using the same specimens as in the creep tests, a tensile test was performed at room temperature with a strain rate of 0.5% per second, and the numerical analysis results were fitted. It was determined that $\alpha =10,{\sigma }_0=60\left[\mathrm{M}\mathrm{P}\mathrm{a}\right]$, and $m=22$.

Figure 3. Experimental and Analytical Stress-Strain curve for determining non-linear coefficient.

Figure 3. Experimental and Analytical Stress-Strain curve for determining non-linear coefficient.

4. Numerical Procedures

4.1. Analysis Algorithm for the Matrix Resin

In this study, a generalized Maxwell model shown in Figure 4 was employed as the viscoelastic model, comprising 15 Maxwell elements. The material constants used were determined as previously described. Analyses were conducted using Abaqus/Standard 2020, and the analysis algorithm for the epoxy resin matrix was implemented using the Abaqus/Standard user subroutine UMAT. The UMAT computational flow was shown in Figure 5. First, the strain increment is taken as the viscoelastic strain increment, and a temporally stress increment is calculated. From this temporally stress increment, the temporally spring extension and the temporally stress increment of the Maxwell elements are determined, and the deformation of the dashpot is calculated. Using the temporally spring extension and temporally dashpot deformation, the extension of the spring is computed, and the stress increment for each Maxwell element is determined. By summing the stress increments of each Maxwell element, the total stress increment is obtained. Additionally, the dissipated energy is calculated from the dashpot deformation and stress, and the damage coefficient is determined and carried over to the next increment. The stress of each Maxwell element and the total stress are carried over to the next increment and used in the calculation of the nonlinear coefficients and the dashpot elongation. When the damage coefficient (D) reaches 0.25, the element is considered to have failed, and element deletion is performed. This ensures that the failed elements no longer carry any stress. For the damage coefficient (Eq. 25), following the simple tensile test by Kudo et al. [23], the fracture entropy was set to 20 $\mathrm{k}\mathrm{J}/{\mathrm{m}}^3·\mathrm{K}$.

$$\mathit{E}={\displaystyle \frac{E_0}{\left(1+\nu \right)\left(1-2\nu \right)}}\times \left[\begin{array}{cccccc}\left(1-\nu \right)& \nu & \nu & 0& 0& 0\\ \nu & \left(1-\nu \right)& \nu & 0& 0& 0\\ \nu & \nu & \left(1-\nu \right)& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{1-2\nu }{2}}& 0& 0\\ 0& 0& 0& 0& {\displaystyle \frac{1-2\nu }{2}}& 0\\ 0& 0& 0& 0& 0& {\displaystyle \frac{1-2\nu }{2}}\end{array}\right]$$

(21)

$$g={\displaystyle \frac{1}{1+\alpha {\left(\frac{{\sigma }_{mises}}{{\sigma }_0}\right)}^m}}$$

(22)

$$∆s=\sum _{j=1}^{15}{\displaystyle \frac{\frac{{\sigma }_{temp}\left(i,j\right)+\sigma \left(i,j\right)}{2}·∆b\left(i,j\right)}{T}}$$

(23)

$$s=s_{old}+\sum _{i=1}^6∆Entropy\left(i\right)$$

(24)

$$D=s·{\alpha }_d$$

(25)

4.2. Analysis Model and Conditions

For the CFRP model, a three-dimensional unit cell analysis model of CFRP, as shown in Figure 6, was utilized. The dimensions of the 3D unit cell were 39 μm × 39 μm × 0.3 μm. The carbon fibers had an elastic modulus of 14 GPa, a Poisson's ratio of 0.3, and a diameter of 6 μm, with 30 fibers included at a volume fraction of 56%. This study did not consider the fiber/matrix interface. The element type used was an 8-node brick element. The epoxy resin matrix was assumed to be a nonlinear viscoelastic material, following a constitutive equation that considers entropy damage. The analysis model applied periodic boundary conditions to all edges using the key degree of freedom method [3,8,12,34,35,36,37,38,39].

For the cyclic tensile analysis, as shown in Figure 7, Dummy-Y was subjected to repeated tensile loading in the Y-axis direction. Displacement control was used to apply strains of 0.1 %, 0.2 %, 0.3 %, 0.4 %, 0.5 %, 0.7 % and 0.8% for 5, 10, 20, 40, and 80 cycles. The average entropy increase of all elements at the end of each cycle was calculated, and after completing the cycles, tensile fracture was induced at a strain rate of 0.1 %/s to determine the residual strength as the maximum stress at fracture.

For the cyclic out-of-plane shear analysis, as shown in Figure 8, the central node was pinned, and a Dummy-X was subjected to repeated shear loading in the Z-direction. Displacement control was used to apply strains of 0.1 %, 0.2 %, 0.3 %, 0.4 %, 0.5 %, 0.6 %, 0.7% and 0.8% for 5, 10, 20, 40, and 80 cycles. The average entropy increase of all elements at the end of each cycle was calculated, and after completing the cycles, tensile fracture was induced at a strain rate of 0.1 %/s to determine the residual strength as the maximum stress at fracture.

5. Results

As an example of the analysis results, Figure 9 shows the damage state after applying repeated tensile loading at a strain of 0.8% for 10 cycles, 20 cycles, 40 cycles, and 80 cycles. The maximum damage value is set to D = 0.25, indicating that the damage coefficient increases from blue to red. Similarly, Figure 10 shows the damage state after applying repeated out-of-plane shear loading at a strain of 0.8% for 10 cycles, 20 cycles, 40 cycles, and 80 cycles. Again, the maximum damage value is set to D = 0.06, indicating that the damage coefficient increases from blue to red. In both analyses, as the number of load cycles increases, the damage also increases.

The relationship between residual strength and entropy derived from cyclic tensile and cyclic out-of-plane shear analyses was fitted with ${\alpha }_o\left(s\right)$ as a linear function and is shown in Figure 11 for cyclic tensile analysis and Figure 12 for cyclic out-of-plane shear analysis. In both the cyclic tensile analysis and the cyclic out-of-plane shear analysis, it was found that an increase in entropy leads to a decrease in residual strength. The determined ${\alpha }_o\left(s\right)$ for each analysis are presented in Eq. 26 and Eq. 27. The residual strength was standardized.

$${\alpha }_o\left(s\right)=-0.14s+1$$

(26)

$${\alpha }_o\left(s\right)=-0.39s+1$$

(27)

6. Conclusions

In this study, we conducted micro-scale numerical simulations to evaluate the residual strength of CFRP subjected to cyclic transverse and out-of-plane shear loading fatigue using a finite element method. By incorporating an entropy-based damage criterion for the matrix resin, we established a linear relationship, which is ${\alpha }_o\left(s\right)$ in Hashin’s’s criterion, between entropy accumulation and strength degradation. The simulations demonstrated that increased entropy, resulting from cyclic loading, directly correlates with decreased residual strength in CFRP. Material constants for the epoxy resin matrix were determined through creep and tensile tests, and a generalized Maxwell model was used to represent viscoelastic behavior. The Abaqus/Standard 2020 software facilitated these analyses with the epoxy resin matrix behavior implemented via a UMAT subroutine. The linear relationship between entropy and residual strength was validated for both cyclic tensile and shear loading conditions. This suggests that using entropy generation allows for a comprehensive assessment of damage under various loading histories. This approach enhances experimental insights with numerical predictions, offering a comprehensive understanding of CFRP strength degradation under fatigue loading. The findings underscore the potential of micro-scale numerical simulations to accurately predict meso-scale strength degradation, thus contributing significantly to the material design and durability assessment of CFRP in practical applications.