Analysis of Stress Intensity Factor of a Fibre Embedded in a Matrix

Abstract

The analytical or numerical determination of the stress intensity factor (SIF) in cracked bodies usually assumes the body to be isolated. However, in fibre-reinforced composites, the fibre, which is the main load-carrying component, is embedded in a matrix. To clarify the effect the embedding matrix has on the SIF of the fibre, we propose a 3D computational model of an orthotropic fibre embedded in an isotropic matrix, and compute the SIF using the J-integral method. A parametric analysis based on dimensionless variables explores the effect of the fibre–matrix stiffness ratio as well as the effect of the degree of elastic orthotropy of the fibre. The results show that the SIF is strongly influenced by both factors, and that the matrix reduces the SIF by limiting the crack opening.

1. Introduction

Composites reinforced with technical fibres (carbon, glass, etc.) exhibit the best mechanical performance among structural materials in terms of specific properties (stiffness and strength-to-weight ratio) and are currently used in many industries such as aerospace or wind turbines. However, their low fracture toughness translates into brittle behaviour that often leads to catastrophic failure without prior damage symptoms [1]. Understanding the strength and toughness of composites is the subject of current research. In particular, how the tensile strength of fibres translates into that of the composite still requires clarification [2].

In spite of the small diameter of these technical fibres, typically below 15 microns, their very high strength leads to Irwin’s lengths much smaller than the diameter; that is, they behave like brittle materials. In this scenario, the stress intensity factor, SIF or $K_I$, which depends on the existing surface flaws, geometry, and fibre toughness (referred to as critical SIF or $K_{IC}$) governs the fibre strength [3].

The critical SIF of fibres has been investigated by combining experimental methods with some type of numerical approach to determine the SIF [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]. For instance, Ogihara et al. [4] introduced various types of notches with straight fronts on carbon fibre mono-filaments using a focused ion beam machining system. Then, the authors used the virtual crack closure method to calculate the SIF of the free carbon fibre under tensile load in both isotropic and orthotropic cases to perform the data reduction in the experimental results. Herraez et al. [5] used a different numerical approach based on the J-integral to compute the SIF for different crack lengths. Then, they characterised the strength and toughness of carbon AS4, E-glass, and Kevlar KM2 via tensile tests on notched and un-notched fibres. In both cases, the SIF was computed assuming a free mono-filament, so the effect of the matrix was neglected.

A single fibre in a composite under longitudinal tension can be idealised as a cylinder embedded in a homogeneous material, the matrix, which is loaded axially. Currently, there are numerical and analytical results available in the literature concerning the SIF of isotropic cracked bars, rods, or beams with cylindrical cross-sections, loaded in tension or bending, and either under static or fatigue loading [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]. Nonetheless, to the authors’ best knowledge, no previous work has been carried out that accounts for the effect an embedding material may have on the SIF.

As a result, this current study employs high-fidelity computational micromechanics modelling to compute the SIF along the crack front of a straight-edge crack in a fibre subjected to tension loading while embedded in an isotropic matrix. First, the method is validated by comparing the results of a free isotropic fibre with previously published results. Following that, a parametric study is carried out to investigate the effect of the fibre–matrix stiffness ratio on the SIF. Finally, we investigate how the longitudinal-to-transverse stiffness ratio of the fibre influences the SIF. The findings presented in this paper contribute to a better understanding of the fracture of elastic bodies that extends beyond composite materials. Furthermore, this work provides closed-form equations to describe the SIF for cylindrical-shaped components with a crack surrounded by an elastic and homogeneous medium. These equations are applicable in other practical cases, such as the glass fibre composite rebars reinforcing concrete.

2. Methodology

2.1. Numerical Evaluation of SIF

This work relies on a computational micromechanics simulation performed in Abaqus Standard [21]. The model consisted of a notched fibre embedded in an isotropic matrix subjected to a tensile axial force (Figure 1a), while the lateral surface area was free from stresses. The fibre included a straight-fronted edge crack at $x_1=L/2$, with varying depths (a) as indicated by the grey area in Figure 1b. Matrix and fibres are modelled as elastic materials with a perfect fibre/matrix adhesion using tie constraints at the interface.

The size of the model was decided based on a parametric study that targeted the best trade-off between computational effort and accuracy. The axial length, L, had to be long enough to guarantee that the remote stresses at the upper and bottom faces were not influenced by the presence of the crack. This led to $L=10D$, where D is the fibre diameter. We selected a large enough thickness of the surrounding matrix, t, to avoid severe deformation in the boundaries of the matrix around the crack region. This criterion prevailed over achieving a realistic fibre volume fraction and led to $t=D$.

To reduce the computational costs, only one-quarter of the model was simulated by taking advantage of the symmetry boundary conditions along the $x_1$$x_2$ and $x_2$$x_3$ planes (Figure 1a). A uniform normal traction stress ($\sigma =1$) was imposed in the upper part of the model. Since the model is linear, SIF and remote stress are proportional. The fibre and matrix were meshed separately. We used a swept and structured mesh around the crack tip to enrich the discretization and capture the high stress gradients therein (Figure 2). The element size, after a convergence study, varied in the range of $D/1500$ (at crack tip) to $D/100$ (in far-field). To account for the stress singularity, quadratic isoparametric elements (C3D20) with the inverse square root singularity were used at the crack tip.

The stress intensity factor evaluation was carried out using the J-integral [22] method. In this approach, the energy release rate is obtained by integration along the contour $\Gamma$ around the crack tip (Figure 3) as follows:

$$J={\int }_{\Gamma }\left(Wd{x}_1-t·\frac{\partial u}{\partial x_2}ds\right)$$

(1)

where $x_1$ and $x_2$ are rectangular coordinates to the crack front, $x_1$ being perpendicular to the crack surface, W is the elastic strain density, t is the traction vector, u is the displacement vector, and $ds$ is an increment of arc length along any contour $\Gamma$. Due to the geometry of the model, the J-integral values are not constant along the crack front. Thus, this work concentrated on the maximum value attained at the crack tip (point A in Figure 1b).

For linear elastic problems of crack opening in mode I, the J-integral equals the energy release rate ($G_I$). The path independence [22] of J-integral was verified and a total of 21 concentric contours surrounding the crack tip were used, as shown in Figure 2 (magnified area). For an isotropic material, the stress intensity factor in mode I, $K_I$, can be calculated from the J-integral as follows:

$$J=\frac{{K_1}^2(1-{\nu }^2)}{E}$$

(2)

where E is the fibre Young’s modulus and $\nu$ the Poisson’s ratio. For an orthotropic material, the J-integral and the $K_I$ are related through the following [23]:

$$J={K_I}^2\sqrt{\frac{a_{11}{a}_{22}}{2}}{\left[\sqrt{\frac{a_{22}}{a_{11}}}+\frac{2a_{12}+a_{66}}{2a_{11}}\right]}^{1/2}$$

(3)

where $a_{ij}$ stands for the compliance tensor of the orthotropic solid under plane strain

$$a_{11}=\frac{1}{E_{11}},a_{22}=\frac{1}{E_{22}},a_{66}=\frac{1}{G_{12}},a_{12}=-\frac{{\nu }_{12}}{E_{11}}$$

(4)

and $E_{11}$ and $E_{22}$ are the fibre Young’s modulus in longitudinal and transverse directions, respectively, and $G_{12}$ is the fibre shear’s modulus. Finally, the non-dimensional stress intensity factor (${\widehat{K}}_I$) is calculated as follows:

$${\widehat{K}}_I=\frac{K_I}{\sigma \sqrt{\pi a}}$$

(5)

where $\sigma$ is the far-field stress.

2.2. Non-Dimensional Parametric Analysis

For an embedded isotropic fibre, the independent elastic constants that have an effect on the SIF are the modulus of elasticity and the Poisson’s ratio of the fibre, $E^f$ and ${\nu }^f$, and of the matrix, $E^m$ and ${\nu }^m$. We performed a non-dimensional analysis based on the definition of the parameter $\alpha$, representing the fibre-to-matrix stiffness ratio:

$$\alpha =\frac{E_1^f-E^m}{E_1^f+E^m}$$

(6)

with $\alpha$ ranging between $\alpha =-1$ (stiffness of the fibre negligible with respect to the matrix, $E^f\ll E^m$) and $\alpha =1$ ($E^f\gg E^m$, which corresponds to the free fibre). For $\alpha =0$, the fibre and matrix have the same modulus of elasticity, so this case corresponds to a homogeneous cylindrical material with an internal crack.

In an orthotropic fibre, the number of elastic constants increases to nine: modulus of elasticity $(E_1^f,E_2^f,E_3^f)$, Poisson’s ratios $({\nu }_{12}^f,{\nu }_{13}^f,{\nu }_{23}^f)$, and shearing moduli $(G_{12}^f,G_{13}^f,G_{23}^f)$. The subscripts 1, 2, and 3 refer to longitudinal and the two transverse directions, respectively. For the sake of simplicity, the fibre was considered transversely isotropic ($E_2^f=E_3^f$, ${\nu }_{12}^f={\nu }_{13}^f$, $G_{12}^f=G_{13}^f$), and some of the elastic constants took a fixed value (${\nu }_{12}^f=0.3$, ${\nu }_{23}^f=0.4$, $G_{12}^f=11.3$, ${\nu }_{12}^f={\nu }^m=0.3$) [24]. Therefore, the elastic constants involved in the parametric study were $E_1^f,E_2^f$, and $E^m$, which were combined in the non-dimensional parameters $\alpha$, defined above in Equation (6), and $\beta$:

$$\beta =\frac{E_1^f-E_2^f}{E_1^f+E_2^f}$$

(7)

where $\beta$ ranges between $\beta =-1$ ($E_1^f\ll E_2^f$) and $\beta =1$ ($E_1^f\gg E_2^f$). The case $\beta =0$ corresponds to the isotropic fibre.

Due to the positive definiteness of the strain energy, the elasticity tensor (compliance matrix) in linear elastic materials must be positive definite [25]. To ensure this, some constraints must be applied to the engineering elastic constants:

$$E>0,G>0,-1<\nu <0.5$$

(8)

for isotropic elasticity, and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E_1,E_2,E_3,& G_{12},G_{13},G_{23}>0\hfill \\ \hfill {\nu }_{12}<& {\left(\frac{E_1}{E_2}\right)}^{1/2}\hfill \\ \hfill {\nu }_{13}<& {\left(\frac{E_1}{E_3}\right)}^{1/2}\hfill \\ \hfill {\nu }_{23}<& {\left(\frac{E_2}{E_3}\right)}^{1/2}\hfill \\ \hfill 1-{\nu }_{12}{\nu }_{12}-{\nu }_{23}{\nu }_{32}& -{\nu }_{31}{\nu }_{13}-2{\nu }_{21}{\nu }_{32}{\nu }_{13}>0\hfill \end{array}\end{array}$$

(9)

for the orthotropic elastic case. These constrains arise due to thermodynamic admissibility; their violation leads to a non-positive strain energy for certain load cases. Based on the linear elastic conditions, stability restricts $\beta$ to $-0.5⩽\beta <1$.

3. Results and Discussion

The numerical approach was verified using two cases of isolated cylindrical bars for which the corresponding SIF has been previously published (free isotropic bar with a surface crack [4,6,26] and a bar with internal crack [27]). This verification is described in the Appendix A.

3.1. Isotropic Fibre Embedded in a Matrix

Figure 4a plots the evolution of the non-dimensional SIF for an isotropic fibre embedded in a matrix for different fibre/matrix stiffness ratios ($\alpha$ equal to −0.95, −0.5, 0.0, 0.5, 0.8, 0.95, and 0.98) and relative crack lengths ($a/D$ equal to 0.07, 0.1, 0.2, 0.3, 0.4, and 0.45). The SIFs of an isotropic free fibre (grey continuous line) were added as a baseline, for comparison. The trend is that an increase in matrix stiffness (decrease in $\alpha$) leads to a significant drop in the SIF of the cracked fibre, especially for longer cracks. Even for a low stiffness matrix ($\alpha$ = 0.98 or $E^m=0.01E^f$), the SIF is reduced by −7% when $a/D$ = 0.2 or by −25% when $a/D$ = 0.4, with respect to the free fibre case. Figure 4b depicts the relative change in SIF for different fibre–matrix stiffness ratios ($\alpha$) compared to the free fibre.

Figure 5 shows the maximum non-dimensional crack mouth opening displacement (CMOD) for the embedded fibre with $\alpha$ being equal to 0.98 and −0.95, as well as that of the free fibre (grey continuous line) for reference. The trends of the curves are very similar to those of the SIF and indicate that the surrounding material restricts the crack opening even in the case of a low stiffness matrix. In fact, the CMOD is negligible when the matrix stiffness is large compared with that of the fibre ($\alpha$ = −0.95, Figure 5). Note that the CMOD shown in Figure 5) corresponds to a unit remote stress.

The non-dimensional SIF for an isotropic fibre embedded in a matrix as a function of the $\alpha$ parameter can be approximated by the following equation, which was obtained by fitting the results of this investigation as follows:

$$\begin{array}{cc}\hfill {\widehat{K}}_I& =A+B\left(\frac{a}{D}\right)+C{\left(\frac{a}{D}\right)}^2\hfill \\ \hfill A& =A_1+\frac{A_2}{{\left(1-\alpha \right)}^{A_3}}\hfill \\ \hfill B& =B_1+\frac{B_2\alpha }{{\left(1-\alpha \right)}^{B_3}}\hfill \\ \hfill C& =C_1+\frac{C_2}{{\left(1-\alpha \right)}^{C_3}}\hfill \end{array}$$

(10)

The use of this equation produces data that deviate by less than $2.8\%$ from any of the FE data points we have examined. The values of the coefficients of Equation (10) are listed in

Table 1. Coefficients for the fitting equation given in Equation (10).

Coefficients of Equation (10)
$A_1$	1.187	$B_1$	−0.4268	$C_1$	−0.4647
$A_2$	−0.4349	$B_2$	−0.1893	$C_2$	0.7246
$A_3$	−0.3255	$B_3$	−0.07555	$C_3$	0.4295

. The usefulness of this equation extends to elastic cylindrical cracked bars surrounded by a homogeneous and isotropic material, provided that the elastic constraints specified in the previous section are respected.

3.2. Orthotropic Free Fibre

Figure 6a shows the non-dimensional SIF for an orthotropic free fibre for different longitudinal/transverse stiffness ratios ($\beta$ equal to −0.5, 0.0, 0.5, 0.8, 0.9, and 0.99) and relative crack lengths ($a/D$ equal to 0.07, 0.1, 0.2, 0.3, 0.4, and 0.45). The case $\beta$ = 0 corresponds to the free fibre with isotropic modulus (the fibre is not completely isotropic as ${\nu }_{23}=0.4$ and ${\nu }_{12}={\nu }_{13}=0.3$). The explored range of $\beta$ involves transverse modulus, $E_2^f$, ranging from three times the longitudinal modulus ($\beta =-0.5$) to being negligible ($\beta =1$). The curves indicate that the degree of orthotropy has a strong influence on the SIF, and that the trend is opposite depending on whether the transverse modulus is smaller ($\beta >0$) or larger ($\beta >0$) than the axial modulus. For example, for $\beta =-0.5$ ($E_2^f=3E_1^f$), the SIF is increased by +30% over the isotropic case; while $\beta =0.5$ ($E_2^f=0.33E_1^f$) involves a decrease of −25%. The relative change in the SIF of orthotripc fibres with respect to the isotropic case does not depend significantly on the crack length (Figure 6b).

The finite element results of the non-dimensional SIF for an orthotropic free fibre were fitted to the following equation:

$$\begin{array}{cc}\hfill {\widehat{K}}_I& =0.9257+A+B\left(\frac{a}{D}\right)+C{\left(\frac{a}{D}\right)}^2+D{\left(\frac{a}{D}\right)}^3\hfill \\ \hfill A& =A_1\beta +A_2{\beta }^2+A_3{\beta }^3+A_4{\beta }^4\hfill \\ \hfill B& =B_1+B_2{\beta }^1+B_3{\beta }^2+B_4{\beta }^3\hfill \\ \hfill C& =C_1+C_2\beta +C_3{\beta }^2\hfill \\ \hfill D& =D_1+D_2\beta \hfill \end{array}$$

(11)

This equation produces data that deviate by less than 2.0% from any of the FE data points we have examined. The values of the coefficients of Equation (11) are listed in

Table 2. Coefficients for the fitting of equation given in Equation (11).

Coefficients of Equation (11)
$A_1$	−0.9022	$B_1$	0.932	$C_1$	−4.647	$D_1$	17.67
$A_2$	0.3677	$B_2$	−0.7516	$C_2$	4.415	$D_2$	−10.7
$A_3$	1.833	$B_3$	1.81	$C_3$	−1.928
$A_4$	−1.993	$B_4$	−1.727

. The application of this equation can be extended to elastic cylindrical cracked bars with orthotropic properties provided that the elastic constraints specified in Equation (9) are met.

3.3. Orthotropic Fibre Embedded in a Matrix

This subsection presents the effect a surrounding material has on the SIF of an orthotropic cracked fibre. This study involves both parameters $\alpha$ and $\beta$ ($\beta$ equal to −0.5, 0.0, 0.5, 0.8, and 0.9). We investigated two cases of fibre–matrix stiffness ratios: $\alpha$ equal to −0.95 and 0.95.

The outcome of this study (Figure 7a,b) reinforces the trends presented in the previous cases: (i) the presence of the matrix reduces the SIF and this effect scales with the stiffness of the matrix, and (ii) the degree of orthotropy of the fibre has a strong influence on the SIF, increasing or decreasing it depending on whether the transverse stiffness is higher or lower, respectively, than the axial stiffness. The strong effect of surrounding material on the SIF should be attributed to the constraints imposed by the matrix on the crack opening (Figure 5).

The complete set of results presented is directly applicable to the understanding of failure in unidirectional fiber-reinforced composites. For instance, a widely-used structural carbon composite combines an AS4 carbon fibre with an epoxy matrix. For that case, $\alpha =0.95$ and $\beta =0.89$ (carbon fibres are non-isotropic [1]). An epoxy reinforced with glass fibres (i.e., S2/MTM44-1) corresponds to the case $\alpha =0.91$ and $\beta =0.0$ (glass fibres are elastically isotropic).

However, the computational model relies on a significant simplification of the physics of the problem. Fibre (albeit not carbon fibres) and matrix are considered homogeneous elastic materials. The fibre is considered to be surrounded by a homogeneous material with the mechanical properties of the matrix. In reality, the fibre is surrounded by a heterogeneous material including nearby fibres which will lead to a higher equivalent stiffness of the embedding material, accentuating the crack opening constraint we have described and, thus, further reducing the SIF. The matrix is assumed to surround the fibre while most probably penetrating the small surface notches or defects in the fibres, thus further constraining the crack opening and reducing the SIF.

In spite of these simplifications, analytical or numerical models of fibre failure and/or fragmentation in a composite, relying on the stress intensity factor, will benefit from the information presented. For example, current fragmentation models rely on the experimental characterization of the fibre strength using a Weibull distribution. This characterization is performed on dry isolated fibres. The present work demonstrates that this characterization cannot be transferred to the fragmentation model while ignoring the effect of the matrix.

In addition, the situation described in this investigation, a cylindrical surface-notched specimen surrounded by a material of different stiffness, goes beyond the field of fibre-reinforced composites.

The computational FEM model developed here can be extended to include material non-linearities. In particular, the model can be refined to allow for more accurate modelling of micro-mechanisms by introducing fibre–matrix interface debonding and/or matrix plasticity. The impact of these phenomena will be explored in a future paper.

4. Conclusions

This work presented a computational micromechanics analysis of the stress intensity factor of a single fibre with a straight-fronted edge crack. The fibre was considered elastic and either isotropic or orthotropic, and free or embedded in an elastic material of different stiffness. The study was performed using non-dimensional variables for the fibre–matrix stiffness ratio and fibre longitudinal/transverse stiffness ratio. Only the modulus of elasticity of fibre and matrix were explored, the rest of the elastic parameters were kept constant throughout the study. The stress intensity factor evaluation was carried out using the J-integral.

Embedding the fibre in an elastic matrix reduced the SIF over the free-fibre case. The higher the matrix stiffness was, the more substantial the decrease in SIF. This effect is due to the constraint of the matrix over the crack opening. For instance, for the model with a very low stiffness matrix ($E^m=0.01E^f$), the SIF is reduced by −7% at $a/D$ = 0.2 or by −25% at $a/D$ = 0.4, with respect to the free-fibre case. We proposed a general non-dimensional equation for the non-dimensional SIF as a function of the non-dimensional crack length and the ratio between the stiffness of the cracked cylinder and that of the surrounding material.

The elastic orthotropy of the fibre greatly affected the SIF in that it increased for the case of a transverse stiffness higher than that of the axial, and decreased for the opposite case. For example, for a fibre with $E_2^f=3E_1^f$, the SIF is increased by +30% over the isotropic case ($E_2^f=E_1^f$); while a fibre with $E_2^f=0.33E_1^f$ involves a decrease of −25%. An equation obtained by fitting was also provided for this case.

Although the study presented here was motivated by an interest in improving the understanding of fibre failure in fibre-reinforced composite materials, the findings and the equations also provide a contribution to the broader field of SIF in elastic cracked bodies.