# Damage onset modeling in woven composites based on a coupled stress and energy criterion

Aurélien Doitrand, Christian Fagiano, Nicolas Carrere, Vincent Chiaruttini, Martin Hirsekorn

► To cite this version:

Aurélien Doitrand, Christian Fagiano, Nicolas Carrere, Vincent Chiaruttini, Martin Hirsekorn. Damage onset modeling in woven composites based on a coupled stress and energy criterion. Engineering Fracture Mechanics, 2017, 169, pp.189-200. 10.1016/j.engfracmech.2016.11.021. hal-01474372

HAL Id: hal-01474372

https://hal.science/hal-01474372

Submitted on 22 Feb 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Damage onset modeling in woven composites based on a coupled stress and energy criterion

Aurélien Doitrand<sup>a,*</sup>, Christian Fagiano<sup>a</sup>, Nicolas Carrère<sup>b</sup>, Vincent Chiaruttini<sup>a</sup>, Martin Hirsekorn<sup>a</sup>

The crack onset configuration at damage onset in a four-layer plain weave glass fiber/epoxy matrix composite is studied at the mesoscopic scale using a coupled criterion based on both a stress and an energy condition. The possible crack shapes are selected based on optical microscope observations of damage mechanisms on a specimen edge during a tensile test. The crack location, length and orientation, the decohesion length and the strain at damage onset are determined. The damage onset strain is underestimated compared to the experimental value determined by acoustic emission if only a stress criterion is considered. The coupled stress and energy criterion leads to a more reasonable estimate of strain at damage onset.

## 1. Introduction

Textile composites are increasingly used for aeronautical and automotive advanced structural applications. An advantage of woven composites is the large variety of the fiber reinforcements, whose architecture can be varied throughout the structure. Therefore, fewer assembly operations are required compared to classical laminated composite structures, which reduces both production costs and the number of weak points in the structure.

Woven composites offer a high potential for material design, which can be optimized using design tools able to describe the evolution of the mechanical behavior from damage onset to final failure of the materials. Macroscopic phenomenological models have been developed for the prediction of damage evolution and failure in 2D [1,2] and 3D [3–7] textile composites. However, expensive and time consuming experimental identifications of the model parameters are required. Moreover, these parameters have to be re-identified each time the fiber architecture or the constituents change. The number of tests may be reduced by using more predictive models at the mesoscopic scale, which take into account the fiber reinforcement architecture.

Different meso-scale approaches have been used to model the evolution of the material behavior with growing damage [8–12]. An analytical model using a damage mosaic laminate model was proposed by Gao et al. [8]. It allows the calculation of the effective Young's modulus of the damaged composite. However, yarn undulation, which has a significant influence on the damage location [11,13], is not described by such a model. In most published studies, a Finite Element (FE) approach based on continuum damage mechanics (CDM) is used (e.g., [9,11,14]). It consists in detecting damaged elements using a

<footnote><sup>a</sup> ONERA - The French Aerospace Lab, F-92322 Châtillon, France</footnote>

<footnote><sup>b</sup> LBMS, ENSTA Bretagne, Université de Brest, 2 rue F. Verny, F-29806 Brest Cedex, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: aurelien.doitrand@onera.fr (A. Doitrand).</footnote>
## Nomenclature

### Energy criterion

$\Delta W_k$ change in kinetic energy

$\Delta W$ change in potential energy

$W(0)$ potential energy of the undamaged state

$W(d)$ potential energy of the cracked state

$G^c$ critical energy release rate

$\Delta S$ crack surface

$h$ specimen thickness

$E^{eq}$ equivalent elastic modulus

$A(d)$ normalized incremental energy release rate

### Stress criterion

$f$ stress based failure criterion

$f_c$ criterion for transverse cracking

$f_d$ criterion for decohesions

$Y_t$ transverse ply strength

$Z_t$ out-of-plane ply strength

$S_{12}^R$ in-plane shear strength

$S_{13}, S_{23}^R$ out-of-plane shear strengths

$p_{12}, p_{13}, p_{23}$ shape parameters

$Y_c$ transverse compressive ply strength

### Crack shape

$d = (d_1, \dots, d_n)$ parameters describing the crack shape

$(x, y)$ coordinates of the crack center

$(x_c, y_c, z_c)$ coordinates of the crack center localized using the stress criterion

$(x^*, y^*)$ crack location at damage onset determined using the energy criterion

$d_c^*$ crack length

$d_s^*$ crack length at damage onset for a crack localized using the stress criterion

$d_c^t$ crack length at damage onset determined using the energy criterion

$\theta$ crack orientation

$\theta^*$ crack orientation at damage onset

$\mu$ decohesion length

$\mu^*$ decohesion length at damage onset

### Strains

$\epsilon_e$ strain measured experimentally

$\epsilon_e^c$ strain at damage onset determined experimentally

$\epsilon$ applied global strain

$\epsilon_{energy}$ strain calculated using the energy criterion

$\epsilon_{energy}^c$ strain at damage onset calculated using the energy criterion

$\epsilon_s^s$ strain at damage onset for a crack localized using the stress criterion

$\epsilon_s$ strain at damage onset obtained using the stress criterion only

### Material properties

$E_m$ Young's modulus of the matrix

$v_m$ Poisson's ratio of the matrix

$E_f$ Young's modulus of the fibers

$v_f$ Poisson's ratio of the fibers

$E_l$ longitudinal Young's modulus of the yarns

$E_t$ transverse Young's modulus of the yarns

$v_{lt}, v_{lt}$ Poisson's ratios of the yarns

$G_{lt}$ shear modulus of the yarns

stress based failure criterion and reducing the local stiffness of these damaged elements with increasing loading. An advantage of this method is the ease of implementation in FE codes. However, erroneous damage propagation directions can be predicted with these models [14,15]. Moreover, regularization methods [16,17] are required in order to avoid damage pattern dependence on FE mesh, which leads to a non-local damage zone. In this case, a very small mesh size compared to the
geometric features of the composite meso-scale structure is required in order to describe the experimentally observed localized cracks, which may be computationally costly.

An approach, which is closer to experimental observations of damage at the meso-scale, is discrete damage modeling [10,18–20], which consists in describing damage by inserting discrete cracks into the FE mesh. In order to determine the strain level at which a crack is generated and the location where it should be inserted in the FE mesh, a crack onset criterion is required. In most published works about damage modeling in textile composites, a stress based criterion is used [9,11,13,14,21]. It seems that such a criterion is satisfying to qualitatively determine the crack location [20,21]. However, it has been shown that a stress criterion is only a necessary condition for damage onset, but not a sufficient one [22]. Indeed, an additional energy condition has to be used for a correct determination of the strain level and of the damage configuration at damage onset. A coupled criterion with stress and energy conditions has been widely used in 2D simulations in order to determine damage onset for different applications [23–32]. An extension to a 3D architecture has been proposed recently by Leguillon [33]. Martin et al. [32] used such a criterion to study the competition between penetration and deflection of a crack at an interface between two different materials. The same author also studied delamination onset in angle-ply laminates [25] and crack onset in the case of a composite plate with an open hole [26]. A coupled criterion has also been employed in order to determine damage onset [27,28] in the case of a single lap joint. Mantič [29,30] used a coupled criterion in order to study fiber-matrix debonding in composites at the microscopic scale. The assumptions formulated in these works for 2D applications of the coupled criterion are also valid for 3D applications, therefore the coupled criterion can be applied to predict damage onset in textile composites at the meso-scale.

In this work, we present a method to determine the strain level at damage onset and the associated crack configuration in a woven composite consisting of four layers of glass fiber plain weave and epoxy matrix. In Section 2, we present the damage mechanisms observed experimentally and a method to insert the observed cracks into a FE mesh of a mesoscopic Representative Unit Cell (RUC) of the studied composite. The stress and the energy conditions of the coupled criterion used for damage onset determination are presented in Section 3. Then, this criterion is applied in Section 4 in order to determine the crack location, length and orientation, the decohesion configuration and the strain at damage onset.

## 2. Discrete damage modeling

### 2.1. Damage mechanisms observed experimentally

The composite under investigation consists of an E-glass fiber plain weave fabric embedded in an Araldite LY564 epoxy matrix. The four layers of the fiber reinforcement were placed in a steel mold and compacted by tightening the screws that keep the mold closed. A rectangular specimen was tested under monotonic tensile loading in order to determine (i) the damage mechanisms through optical microscope observations on the specimen edge and (ii) the strain at damage onset by means of acoustic emission, as in Refs. [12,13].

The observed damage mechanisms are transverse yarn cracks that extend through the entire yarn thickness and decohesions at the interface between cracked yarns and overlapping yarns (Fig. 1a). Some yarn cracks are perpendicular to the loading direction (Fig. 1b) and others are inclined (Fig. 1c). Bilateral (Fig. 1b) and unilateral (Fig. 1c) decohesions are observed around the crack tips. The damage mechanisms in the material result in a spontaneous release of elastic strain energy that is partly dissipated as waves propagating through the material. These acoustic events were detected by sensors attached to the specimen. The damage onset strain level, $\epsilon_e^c = 6.10 \cdot 10^{-3}$, was determined as the strain for which the first acoustic events are detected (Fig. 2).

### 2.2. Discrete damage modeling

Discrete damage modeling at the meso-scale requires the generation of a FE mesh geometrically representative of the composite. It has been shown in several studies that the reinforcement architecture of the composite influences damage onset and propagation [34,35], as well as the energy release rate [18]. Therefore, it is essential to use a geometrical description and a mesh that are as close as possible to the real fiber reinforcement of the material. Such a geometry can be obtained by simulating the compaction of the dry fabric before resin injection. The method described in [20] is used in this work to generate a RUC containing a compacted and nested geometry of four layers of plain weave. Then, a FE mesh consistent with this geometry and conformal at the interfaces between yarns in contact and between the reinforcement and its matrix complement is obtained using the algorithm developed by Grail et al. [36] (Fig. 3a and b). An advantage of this method compared to voxel meshes [20,37,38] is that the discretized model is close to the initial geometry and contains smooth yarn surfaces. Moreover, the FE mesh is adapted for modeling damage observed experimentally such as decohesion between yarns, which would not be possible with a matrix layer between the yarns, as in [11,14,39,40].

Discrete damage modeling consists in remeshing locally an undamaged mesh (Fig. 3c) in order to represent a crack by the remeshed element faces and then doubling the nodes at the crack surface. The remeshing is done using tools initially developed by Chiaruttini et al. [41] for crack propagation simulation in metallic structures. The crack surface is inserted into the RUC mesh (Fig. 3d) using the mesh intersection algorithm detailed in [41]. For the sake of clarity, the mesh element outline has been removed in Fig. 3c and d.
<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) Microscope observations of the edge of the specimen tested under tensile loading in direction of the warp yarns (x direction). (b, c) Zoom on damaged zones in weft yarns with different crack orientations and bilateral (b) and unilateral (c) decohesions.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Strain at damage onset ($ε_e^c = 6.10 \cdot 10^{-3}$) determined by acoustic emission: amplitude and cumulative energy as a function of the strain measured experimentally $ε_e$.</figcaption></figure>

### 3. A coupled criterion for crack onset

Crack onset requires two separate conditions that must be fulfilled. On the one hand, stress must be sufficiently high to damage the material. On the other hand, sufficient energy is required in order to open a crack. Let $d = (d_1, ..., d_n)$ be the parameters describing the crack shape (e.g. its position, its length, its orientation with respect to the loading direction, etc.). The objective is to determine a set of parameters that describes a crack at damage onset that fulfills both the stress and the energy conditions.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Representative unit cell of the composite under investigation: (a) mesh of the RUC, (b) mesh of the yarns only, weft yarn (in purple) under consideration in (c) the undamaged state and (d) the damaged state where a discrete crack (in red) is inserted in the mesh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 3.1. Energy condition

An energy balance between an elastic state before any crack growth and after a crack onset of area $\Delta S$ leads to:

$$ \Delta W + \Delta W_k + G^c \Delta S = 0 \qquad (1) $$

where $\Delta W$ is the change in potential energy, $\Delta W_k$ the change in kinetic energy and $G^c \Delta S$ the fracture energy ($G^c$ is the critical energy release rate of the material). The condition $\Delta W_k \ge 0$ implies

$$ -\Delta W \ge G^c \Delta S \qquad (2) $$

The energy condition compares the difference between the potential energy of the undamaged state, $W(0)$, and that of the cracked state, $W(d)$, with the energy required to open the crack specified by $d$ (Eq. (3)).

$$ -\Delta W = W(0) - W(d) = \Delta S \cdot h \cdot E^{eq} \cdot A(d) \cdot \epsilon^2 \ge \Delta S \cdot G^c \qquad (3) $$

Under the assumption of elasticity, the potential energy of both the undamaged and the damaged material (and hence their difference $\Delta W$) is proportional to the square of the applied global strain $\epsilon$. Through normalization by some material specific parameters $h$ and $E^{eq}$, a dimensionless function $A(d)$ (the normalized incremental energy release rate [24,27]) can be defined, which only depends on the RUC geometry and on the crack parameters $d$. $h$ is a characteristic length of the geometry (for example the specimen thickness). The equivalent elastic modulus $E^{eq}$ can be chosen equal to the yarn transverse modulus (the role of $h$ and $E^{eq}$ is to normalize the energy release rate and their choice does not have any influence on the final result). $A(d)$ can be determined by generating response surfaces of the potential energy balance, $\Delta W$, to the varying parameters ($d_1, \dots, d_n$) by means of FE calculations in the corresponding damaged RUCs. In order to open a crack described by the set of parameters $d$, the released potential energy, $\Delta W$, must be equal or higher than the energy required to generate the crack surface: $\Delta S \cdot G^c$. The strain at damage onset as a function of the crack parameters determined by the energy criterion is thus given by:

$$ \epsilon_{\text{energy}}^{c}(d) = \sqrt{\frac{G^{c}}{A(d) \cdot E^{eq} \cdot h}} \qquad (4) $$

### 3.2. Stress condition

The stress condition states that a failure criterion $f = f(\sigma)$ must be fulfilled (i.e. $f \ge 1$) in a zone covering the whole area of the crack nucleation. The failure criterion used in this work was initially developed for UD plies in a composite laminate [42]. In this model, different damage modes are taken into account: (1) fiber breaking, (2) transverse and (3) out-of-plane matrix cracking. Since, for the material under investigation, damage onset consists in a transverse yarn crack that can be accompanied by decohesions at the crack tip, damage modes 2 and 3 are respectively taken into account for transverse yarn crack and decohesions. The criteria $f_c$ and $f_d$ for modes 2 and 3, respectively, are given in Eqs. (5) and (6). The transverse failure mode (mode 2) has already been used in [12,20] in order to predict damage onset in woven composites, but without considering an energy condition.
<table><caption>**Table 1**

Values of the material parameters used in the coupled criterion.</caption><thead><tr><th>E<sup>eq</sup> (GPa)</th><th>h (mm)</th><th>G<sup>c</sup> [45] (N/mm)</th><th>Y<sub>t</sub> [46] (MPa)</th><th>S<sub>12</sub><sup>R</sup> [46] (MPa)</th><th>S<sub>23</sub><sup>R</sup> [46] (MPa)</th><th>p<sub>12</sub> [42] (MPa<sup>-1</sup>)</th></tr></thead><tbody><tr><td>9.79</td><td>1.68</td><td>0.118</td><td>35</td><td>72</td><td>45</td><td>8.8 &middot; 10<sup>-3</sup></td></tr></tbody></table>

$$f_c = \left( \frac{\langle \sigma_{22} \rangle^+}{Y_t} \right)^2 + \frac{\tau_{12}}{S_{12}^R (1 - p_{12} \sigma_{22})} + \frac{\tau_{23}}{S_{23}^R (1 - p_{23} (\sigma_{22} + \sigma_{33}))} \quad (5)$$

$$f_d = \left( \frac{\langle \sigma_{33} \rangle^+}{Z_t} \right)^2 + \frac{\tau_{13}}{S_{13}^R (1 - p_{13} \sigma_{33})} + \frac{\tau_{23}}{S_{23}^R (1 - p_{23} (\sigma_{22} + \sigma_{33}))} \quad (6)$$

Direction 1 is the direction of the fibers in the yarn local frame and directions 2 and 3 define the transverse isotropy plane, direction 2 being orthogonal to the crack plane. $Y_t$ and $Z_t$ are the transverse and out-of-plane tensile strength. $S_{12}^R$, $S_{13}^R$ and $S_{23}^R$ are the shear strengths. $p_{12}$, $p_{13}$ and $p_{23}$ are shape parameters that take into account coupling between compressive and shear stresses. These parameters can be directly determined from the transverse compressive ply strength $Y_c$ [42]. Since the yarn behavior is transverse isotropic, the transverse and out-of-plane strengths are equal ($Y_t = Z_t$) [43], as well as the shear strengths ($S_{12}^R = S_{13}^R$) and the shape parameters ($p_{12} = p_{13} = p_{23}$). The values of the different parameters used are summarized in Table 1. Damage occurs when any of both criteria becomes greater than 1, the local damage is thus given by the criterion ($f_c$ or $f_d$) that is maximum.

### 3.3. Application of the coupled criterion

The crack configuration at damage onset can be determined as the set of parameters ($d_1, \dots, d_n$) for which the stress and the energy conditions are fulfilled, minimizing the strain at damage onset. In practice, the stress condition is determined by one linear elastic calculation on the undamaged RUC and the energy condition is computed by linear elastic FE calculations in several damaged RUC with different crack configurations (i.e. by varying ($d_1, \dots, d_n$)). It can be noted that the potential energies, $W(0)$ and $W(d)$, are calculated using meshes containing, respectively, an undamaged or a cracked yarn and having exactly the same topology. These meshes only differ because the nodes describing the crack surface are doubled in the cracked mesh. For all the calculations, periodic boundary conditions, as described in [14], are applied in the fabric plane directions of the RUC, whereas the top and bottom surfaces are left free in order to correctly represent the boundary conditions applied to the composite specimen during a tensile test. The matrix behavior is supposed to be linear elastic, and the mechanical properties provided by the manufacturer are: Young's modulus $E_m = 3.2$ GPa and Poisson's ratio $v_m = 0.35$. The yarn behavior is obtained by micro-meso homogenization, as e.g. in [11], using $E_f = 73.6$ GPa [44] and $v_f = 0.3$ for the fibers. The fiber volume fraction in the yarns is chosen in order to obtain an overall fiber volume fraction in the RUC similar to the fiber volume fraction in the real composite. A transverse isotropic elastic behavior is obtained for the yarns, with $E_f = 41.0$ GPa, $E_t = 9.79$ GPa, $v_{tt} = 0.424$, $v_{lt} = 0.320$ and $G_{lt} = 7.21$ GPa, where the index $l$ refers to the fiber direction and the index $t$ to the transverse direction. The local orientation of the yarn material is calculated separately at each integration point through orthogonal projection of its position on the neutral line of the yarn. The tangent to the neutral line at the projected point defines the axis of transverse isotropy of the yarn material (the fiber direction) at the integration point. The material parameters used in the coupled criterion are given in Table 1. The critical energy release rate $G^c$ is supposed to be a material constant, independent of the fracture mode mixity.

## 4. Damage onset prediction in a woven composite

The criterion presented in the previous section is used in this section in order to determine a crack configuration at damage onset. Since the crack shape is not known a priori, hypotheses can be formulated in order to limit the number of possible crack configurations at damage onset and, therefore, the number of calculations required to compute the coupled criterion. The first hypothesis is that the crack plane is parallel to the fiber direction. This is justified by the fact that the fiber strength is much higher than the matrix strength. Therefore, the parameters describing the crack are (i) the position in the yarn, (ii) the crack dimensions, (iii) the angle of the normal to the crack plane with respect to the loading direction and (iv) the decohesion configuration.

### 4.1. Crack dimensions

In order to limit the number of calculations, we first assume that the crack plane is normal to the loading direction. The influence of the crack orientation will be studied in Section 4.3. In published studies concerning woven composite modeling at the meso-scale, damage onset location is usually detected using only a stress criterion. However, fulfilling the stress cri-
terion is only a necessary condition. Crack onset also requires that the energy criterion is fulfilled for a crack covering a zone where the stress criterion is fulfilled everywhere. If the stress criterion is dominant, the crack shape will thus be given by the area where $f_c \ge 1$. The stress criterion on the yarn surface is shown in Fig. 4 for an applied strain of $6.8 \cdot 10^{-3}$ (as is shown below, this is the minimum strain at which the energy criterion is fulfilled for a crack located at the maximum of the stress criterion (Point C in Fig. 4a)). It can be seen from Fig. 4b that the areas where $f_c \ge 1$ have roughly the shape of semi-ellipses in the plane perpendicular to the loading direction located at the point C. We therefore start our analysis with a small semi-elliptical (half axes 0.07 mm and 0.5 mm) crack perpendicular to the loading direction located around point C. The strain required to fulfill the stress criterion on the whole surface of the crack in the undamaged material is $\epsilon = 2.9 \cdot 10^{-3}$. The insertion of such a crack decreases the potential energy by $\Delta W = 3.52 \cdot 10^{-4}$ N mm. However, the energy required to open this crack is: $\Delta S \cdot G^2 = 1.3 \cdot 10^{-2}$. This energy would only be released at a much higher strain: $\epsilon = 17.7 \cdot 10^{-3}$. At this strain level, the stress criterion would be fulfilled in a much larger zone. We therefore repeated this analysis for the strains at which the stress criterion is fulfilled respectively in areas covering half and the whole yarn thickness ($\epsilon = 3.1 \cdot 10^{-3}$ and $\epsilon = 3.5 \cdot 10^{-3}$). Even for the latter case, the energy criterion is only fulfilled at a much larger strain ($\epsilon = 9.4 \cdot 10^{-3}$) than the stress criterion ($\epsilon = 3.5 \cdot 10^{-3}$). The results are summarized in Table 2. The strain required to fulfill the energy criterion decreases with increasing crack size. Therefore, we conclude that the crack at damage onset, for which both the stress and the energy criteria must be fulfilled, must extend through the whole thickness of the yarn, as observed experimentally. García et al. [47,48] showed that theoretically speaking the extension of a crack after the onset is slightly smaller than the yarn extension. However, the difference is very small (about a few fiber diameters) compared to the yarn width.

Therefore, the following analyses are based on cracks that extend through the whole yarn thickness with straight crack fronts. We first assume that no decohesions are present at the crack tip (this aspect will be addressed in Section 4.4) and that the crack is still perpendicular to the loading direction. In this case, the crack can be described by a single parameter, its length $d_c$.

Several calculations were carried out in a damaged RUC with different crack lengths, and the normalized incremental energy release rate and the strain at damage onset were calculated using the energy criterion (Eq. (4)) for every crack length. Their evolution as a function of the crack length is plotted in Fig. 5. It can be noted that $A(d)$ reaches a maximum for a $d_s^* = 3.4$ mm crack length, corresponding to a damage onset strain: $\epsilon_{energy}^* = 6.8 \cdot 10^{-3}$. Since, at this minimum strain, the stress criterion is fulfilled on the corresponding crack area, the crack length at damage onset is given by the length $d_s$ that maximizes $A(d)$ (respectively minimizes $\epsilon_{energy}$) [31]. In fact, at this minimum strain for damage onset obtained with the energy criterion, the stress criterion is fulfilled everywhere in the whole yarn (minimum $1.17 \cdot 10^{-3}$, see Fig. 4a). The stress criterion thus allows damage onset also at another position in the yarn, not only at its local maximum. In the studied case, the energy criterion determines the strain at damage onset, the location and the crack length being determined by minimizing this strain over the whole yarn, which is typical for tough configurations [48].

## 4.2. Crack length and location

Based on the results presented in the previous section, the crack is supposed to extend to the whole yarn thickness. We first study the case of a crack plane perpendicular to the loading direction. Hence, the crack is described by the parameter set $d = (x, y, d_c)$ where $d_c$ is the crack length and $(x, y)$ the crack location. For several lengths between 0.05 mm and 5 mm, a discrete crack is inserted into the mesh at different positions covering the whole yarn. Then, the strain at damage onset is obtained with the energy criterion (Eq. (4)) for each location and length. Fig. 6a shows, for the studied crack lengths, the minimum strain at damage onset obtained by inserting a crack of the respective length at different positions in the yarn. Fig. 6b shows, for each crack position in the yarn, the minimum strain at damage onset obtained by varying the crack length. It

<figure><img src="image_5.png" /><figcaption>Fig. 4. (a) Stress criterion on the yarn surface for a loading strain of $\epsilon = 6.8 \cdot 10^{-3}$. Its maximum is located at the point C(x_c, y_c, z_c). (b) Loading strain for which the stress criterion locally reaches 1 on a plane perpendicular to the loading direction passing through C.</figcaption></figure>
<caption>**Table 2**</caption>

Dimensions of the three cracks with a semi-elliptical shape and corresponding strain loadings needed to reach the stress criterion on the whole crack surface.

<table><thead><tr><th>Crack</th><th>Small axis (mm)</th><th>Large axis (mm)</th><th>Area (mm²)</th><th>Loading strain</th><th>Strain to fulfill the energy criterion</th></tr></thead><tbody><tr><td>Small</td><td>0.07</td><td>0.5</td><td>0.11</td><td>2.9 · 10<sup>-3</sup></td><td>17.7 · 10<sup>-3</sup></td></tr><tr><td>Half thickness</td><td>0.13</td><td>1.</td><td>0.41</td><td>3.1 · 10<sup>-3</sup></td><td>13.0 · 10<sup>-3</sup></td></tr><tr><td>Whole thickness</td><td>0.3</td><td>1.5</td><td>1.40</td><td>3.5 · 10<sup>-3</sup></td><td>9.4 · 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** (a) Normalized incremental energy release rate and (b) strain at damage onset calculated with the energy criterion as a function of the crack length $d_c$ for cracks located around the local maximum of the stress criterion.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** (a) Minimum of the strain at damage onset over the yarn as a function of the crack length $d_c$. (b) Damage onset strain distribution in the yarn for the crack length minimizing locally the damage onset strain.</figcaption></figure>

should be noted that the damage onset strain reaches a minimum, $\epsilon_{energy}^c = 6.37 \cdot 10^{-3}$, for a $d_c^* = 2.02$ mm crack length and a crack centered at $(x*, y*) = (3.57$ mm, $5.13$ mm). This minimum is smaller than that obtained for cracks centered at the maximum of the stress criterion. In fact, the crack center is located close to, but not exactly at the maximum of the stress criterion ($x_c = 3.31$ mm, $y_c = 6.42$ mm). The energy criterion thus provides the strain at damage onset, as well as the length and the location of the corresponding crack (Fig. 6b). Fig. 7 shows the value of the stress criterion over an area corresponding to the crack located at $(x*, y*)$ and of length $d_c^*$ for a global strain: $\epsilon = \epsilon_{energy}^c = 6.37 \cdot 10^{-3}$. Since it is reached ($f_c \ge 1$) over the whole area, both the stress and the energy conditions are fulfilled at $\epsilon = \epsilon_{energy}^c$, and this strain is the minimum strain for which both conditions are fulfilled.

### 4.3. Crack orientation

As shown in Section 2.1, some transverse yarn crack normals are disoriented with respect to the loading direction. The crack orientation can be estimated by applying the energy criterion to a set of cracks with different orientations. The released potential energy is calculated using the method described in Section 3.3 for several cracks at different orientations from $\theta = -45^\circ$ to $\theta = 45^\circ$. The crack length and position are the same as those identified in Section 4.2. The orientation that minimizes the strain at damage onset is $\theta^* = 0^\circ$. Therefore, the crack is preferentially perpendicular to the loading direction. In
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Failure criterion distribution over the area in the undamaged RUC corresponding to the crack surface determined using the energy criterion.</figcaption></figure>

this case, the crack orientation has been determined separately from the crack length. The validity of this parameter separation will be discussed in Section 5.

### 4.4. Decohesion configuration

As shown in Section 2.1, transverse yarn cracks are often accompanied by decohesions between overlapping yarns. At a $\epsilon = \epsilon_{energy}^c = 6.37 \cdot 10^{-3}$ global strain, for which the onset of a transverse yarn crack without decohesion is predicted by the coupled criterion (see Section 4.2), the stress criterion is not fulfilled at the yarn interface around the crack tip in the undamaged RUC (a maximum of $f_d = 0.43 < 1$ is reached on the yarn surface). Therefore, the imposed loading is not sufficient to generate decohesions, so that, at damage onset, a transverse yarn crack nucleates without decohesion. However, as soon as the crack is nucleated, the stress becomes singular at the crack tip. Therefore, after the crack onset, the stress criterion can be reached at the interface between the yarns close to the crack tip, which allows decohesion nucleation. Martin et al. [32] showed that the decohesion can even nucleate before the crack tip has reached the interface. For simplicity purposes, this configuration will not be studied in this work; therefore, the decohesion is supposed to be generated only when the crack is extended to the whole yarn thickness. The length of the decohesion that appears immediately after the crack onset can be estimated by applying the energy criterion to a set of transverse yarn cracks and decohesions around the crack tip. In this case, the energy balance is established between the state with a transverse crack and a decohesion and the state with a transverse crack only. The released potential energy is calculated using the method described in Section 3.3 for the crack at different orientations studied in Section 4.3, with decohesions of a total length between 0.01 mm and 0.18 mm symmetrical about the crack tip. The strain required to fulfill the energy criterion is calculated using Eq. (4) for each crack configuration. The evolution of this strain as a function of the decohesion length, $\mu$, and the angle, $\theta$, between the crack normal and the applied loading direction, is presented in Fig. 8. It can be seen that for any decohesion length, the strain is minimized for the same crack orientation of $\theta = 0°$ (Fig. 8a). For a given orientation, the energy condition indicates that the smaller the decohesion length, the smaller the onset strain (Fig. 8b). However, the onset strain must be at least equal to $\epsilon_{energy}^c$, otherwise the transverse crack does not nucleate and the stress criterion is not fulfilled at the yarn interface. At this strain level, the released energy is sufficient to generate a decohesion of length: $\mu^* = 0.09$ mm (Fig. 8b). A decohesion of this length is therefore expected to be generated instantaneously after the nucleation of the associated transverse yarn crack. Similar strains fulfilling the energy criterion as a function of decohesion length and crack orientation are obtained in the case of unilateral decohesions. In this case, we can therefore not draw any conclusion about which crack configuration (bilateral or unilateral) is preferred.

## 5. Discussions

### 5.1. Crack localization

In most published studies on woven composite materials, the location of damage onset is usually determined using a stress criterion [11–14,20,21], which seems to provide a correct qualitative estimate of the crack location observed experimentally [20,21]. In the present study, the crack location obtained with the stress criterion ($x_c, y_c$) = (3.31 mm, 6.42 mm) is about 1.3 mm away from the crack location obtained with the coupled criterion. Most of this distance is in the fiber direction,
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Damage onset strain as a function of (a) the crack angle for different decohesion lengths between 0.01 and 0.18 mm (the red curve corresponds to the case of a crack without decohesion) and of (b) the decohesion length for a given crack angle $\theta = 0^{\circ}$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

i.e., parallel to the crack plane. In the studied case, the stress criterion gives therefore a relatively accurate estimate of the location of the crack plane (a distance $\Delta x = 0.26$ mm from the plane predicted by the coupled criterion, the yarn width being 3.6 mm) but the crack is not centered around the maximum of the stress criterion. However, the damage onset strain, obtained using only the stress criterion, is $\epsilon_s = 2.7 \cdot 10^{-3}$, which is about 60% lower than the strain obtained with the coupled criterion $\epsilon_{energy}^c$ and about 55% lower than the experimentally measured damage onset strain ($\epsilon_c^c = 6.1 \cdot 10^{-3}$). The damage onset strain obtained by minimizing $\epsilon_{energy}$ as a function of the crack length at the location determined by the stress criterion is $\epsilon_{energy}^s = 6.8 \cdot 10^{-3}$, which is about 7% higher than $\epsilon_c^c$. The corresponding crack length is $d_s^* = 3.4$ mm, which is much longer than the length of the crack minimizing $\epsilon_{energy}^c$ over the whole yarn ($d_c^* = 2.02$ mm). Hence, the stress criterion without the energy condition yields a strain at damage onset that is much lower than the strain detected experimentally through acoustic emission monitoring. Taking into account both the stress and the energy conditions, a strain much closer to the experimental measurement is obtained.

Nevertheless, we have to point out that the available data does not allow an exact quantitative comparison with experimental data, due to the uncertainties in some of the material parameters used in the model. For example, the fiber Poisson's ratio or the critical energy release rate of the yarn cannot be easily obtained experimentally and are taken from published studies on similar materials. For instance, the value of $G^c$ in a yarn has been determined for a glass/epoxy UD ply with similar constituents [45]. A second point is that in the real material, there may be defects, porosity, or local variations in the reinforcement architecture that may influence damage onset and that are not taken into account in the model. Therefore, for a given set of material parameters, the strain at damage onset is rather overestimated compared to a real material with the same mechanical properties.

## 5.2. Parameter separation

The application of the coupled criterion requires a large number of calculations in order to compute the energy criterion as a function of the parameters describing a crack. The number of calculations grows exponentially with the number of

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** (a) 3D and (b) 2D visualisations of damage onset strain ($\epsilon_{energy}$) as a function of the crack length $d_c$ and the crack angle $\theta$.</figcaption></figure>
parameters describing the crack. Therefore, uncoupling some parameters in the minimization of the damage onset strain would greatly reduce the computational cost. Although the crack length cannot be determined separately from the crack location in the yarn, the crack orientation has been decoupled from the remaining parameters in the analysis presented in Section 4. In order to confirm this hypothesis, the damage onset strain has been calculated as a function of the crack length and of the crack orientation for cracks centered at (x*, y*). The results are presented in Fig. 9, which indicates that the evolution of the damage onset strain with the crack orientation is independent of the crack length, and vice versa. As explained in Section 4.4, for the material under consideration, decohesions are only generated once a transverse crack is already nucleated. Therefore, the decohesion configuration can be determined separately from the transverse crack configuration.

## 6. Conclusion

The classical approach for damage onset prediction in woven composites with a stress criterion leads to an underestimate of the strain at damage onset. A more reasonable estimate is obtained by coupling stress and energy conditions, which also allows the determination of the crack configuration (location, length, orientation and decohesion length). The number of calculations required to determine the configuration of a crack at damage onset grows, in theory, exponentially with the number of parameters describing the crack, which can be computationally expensive. However, the number of calculations may be reduced by uncoupling some parameters. To the authors' knowledge, this is the first time that a coupled criterion with a stress and an energy condition has been applied to a complex 3D case of a woven composite. Future work will cover the determination of damage mechanism kinetics in woven composites using the presented approach for multiple cracks, as well as a study of the influence of the material parameter variability on the prediction of damage onset.

## References

[1] Hochard C, Aubourg PA, Charles JP. Modelling of the mechanical behaviour of woven-fabric CFRP laminates up to failure. Compos Sci Technol 2001;61:221-30.

[2] Barbero EJ, Lonetti P, Sikkil KK. Finite element continuum damage modeling of plain weave reinforced composites. Compos Part B 2005;37(23):137-47.

[3] Maire JF, Chaboche JL. A new formulation of continuum damage mechanics (CDM) for composite materials. Aerospace Sci Technol 1997;1(4):247-57.

[4] Marcin L, Carrère N, Maire JF. A macroscopic visco-elastic-damage model for three dimensional woven fabric composites. In: Proceedings of ECCM13 - 13th European conference on composite materials, Stockholm, Sweden.

[5] Rakotoarisoa C, Laurin F, Hirsekorn M, Maire JF, Olivier L. Development of a fatigue model for 3D woven polymer matrix composites based on a damage model. In: Proceedings of ECCM15 - 15th European conference on composite materials, Venice, Italy [paper 101].

[6] Elias A, Kaminski M, Laurin F, Gornet L. Study of low velocity impact defects in organic interlock woven composite. In: Proceedings of ECCM 16 - 16th European conference on composite materials, Seville, Spain.

[7] Hurmane A, Irisarri FX, Laurin F, Leclercq S, Benzeggagh M. Strength analysis of woven interlock composites subjected to compressive loading: experiments and simulations. In: Proceedings of ECCM 16 - 16th European conference on composite materials, Seville, Spain.

[8] Gao XL, Li K, Mall S. A mechanics-of-materials model for predicting Young's modulus of damaged woven fabric composites involving three damage modes. Int J Solids Struct 2003;40:981-99.

[9] Zako M, Uetsuji Y, Kurashiki T. Finite element analysis of damaged woven fabric composite materials. Compos Sci Technol 2003;63:507-16.

[10] Couegnat G. Approche multiéchelle du comportement mécanique de matériaux composites à renfort tissé [Ph.D. thesis]. Université de Bordeaux 1; 2008. <https://tel.archives-ouvertes.fr/tel-00403885/>.

[11] Melro P, Camanho P, Andrade Pires FM, Pinho ST. Numerical simulation of the non-linear deformation of 5-harness satin weaves. Comput Mater Sci 2012;61(0):116-26.

[12] Doitrand A, Fagiano C, Chiaruttini V, Leroy FH, Mavel A, Hirsekorn M. Experimental characterization and numerical modeling of damage at the mesoscopic scale of woven polymer matrix composites under quasi-static tensile loading. Compos Sci Technol 2015;119:1-11.

[13] Daggumati S, Van Paepegem W, Degrieck J, Xu J, Lomov SV, Verpoest I. Local damage in a 5-harness satin weave composite under static tension: Part II Meso-FE modelling. Compos Sci Technol 2010;70:1934-41.

[14] Lomov SV, Ivanov DS, Verpoest I, Zako M, Kurashiki T, Nakai H, et al. Meso-FE modelling of textile composites: road map, data flow and algorithms. Compos Sci Technol 2007;67(9):1870-91.

[15] Gorbatikh L, Ivanov D, Lomov SV, Verpoest I. On modelling of damage evolution in textile composites on meso-level via property degradation approach. Compos Part A 2007;38:2433-42.

[16] Allix O, Feissel P, Thevenet P. A delay damage mesomodel of laminates under dynamic loading: basic aspects and identification issues. Comput Struct 2003;81:1177-91.

[17] Maimí P, Camanho PP, Mayugo JA, Davila CG. A continuum damage model for composite laminates: Part II Computational implementation and validation. Mech Mater 2007;39:909-19.

[18] Le Page BH, Guild FJ, Ogin SL, Smith PA. Finite element simulation of woven fabric composites. Compos Part A 2004;35:861-72.

[19] Obert E, Daghia F, Ladevèze P, Ballere L. Micro and meso modeling of woven composites: transverse cracking kinetics and homogenization. Compos Struct 2014;117:212-21.

[20] Doitrand A, Fagiano C, Irisarri FX, Hirsekorn M. Comparison between voxel and consistent meso-scale models of woven composites. Compos Part A 2015;73:143-54.

[21] Faes JC, Rezaei A, Van Paepegem W, Degrieck J. Accuracy of 2D FE models for prediction of crack initiation in nested textile composites with inhomogeneous intra-yarn fiber volume fractions. Compos Struct 2015. http://dx.doi.org/10.1016/j.compstruct.2015.12.024. Accepted Manuscript [in press].

[22] Leguillon D. Strength or Toughness? A criterion for crack onset at a notch. Eur J Mech A/Sol 2002;21:61-72.

[23] Leguillon D, Martin E, Sevecek O, Bermejo R. Application of the coupled stress-energy criterion to predict the fracture behaviour of layered ceramics designed with internal compressive stresses. Eur J Mech A/Sol 2015;54:94-104.

[24] Martin E, Leguillon D, Lacroix C. A revisited criterion for crack deflection at an interface in a brittle material. Compos Sci Technol 2001;61:1671-9.

[25] Martin E, Leguillon D, Carrère N. A twofold strength and toughness criterion for the onset of free-edge shear delamination in angle-ply laminates. Int J Sol Struct 2010;47:1297-305.

[26] Martin E, Leguillon D, Carrère N. A coupled strength and toughness criterion for the prediction of the open hole tensile strength of a composite plate. Int J Sol Struct 2012;49:3915-22.

[27] Carrère N, Martin R, Leguillon D. Comparison between models based on a coupled criterion for the prediction of the failure of adhesively bonded joints. Eng Fract Mech 2015;138:185-201.
[28] Moradi A, Carrère N, Leguillon D, Martin E, Cognard JY. Strength prediction of bonded assemblies using a coupled criterion under elastic assumptions: effect of material and geometrical parameters. Int J Adhes Adhes 2013;47:73-82.

[29] Mantič V. Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion. Int J Sol Struct 2009;46:1287-304.

[30] Mantič V, García IG. Crack onset and growth at the fibre-matrix interface under a remote biaxial transverse load. Application of a coupled stress and energy criterion. Int J Sol Struct 2012;49:2273-90.

[31] Martin E, Leguillon D. Energetic conditions for interfacial failure in the vicinity of a matrix crack in brittle matrix composites. Int J Sol Struct 2004;41:6937-48.

[32] Martin E, Poitou B, Leguillon D, Gatt JM. Competition between deflection and penetration at an interface in the vicinity of a main crack. Int J Fract 2008;151:247-68.

[33] Leguillon D. An attempt to extend the 2D coupled criterion for crack nucleation in brittle materials to the 3D case. Theroret App Fract Mech 2014;74:7-17.

[34] De Carvalho NV, Pinho ST, Robinson P. An experimental study of failure initiation and propagation in 2D woven composites under compression. Compos Sci Technol 2011;71,1316-25.

[35] Tang X, Whitcomb J. Progressive failure behaviors of 2D woven composites. J Compos Mater 2003;37:1239-59.
[36] Grail G, Hirsekorn M, Wendling A, Hivet G, Hambli R. Consistent finite element mesh generation for meso-scale modeling of textile composites with preformed and compacted reinforcements. Compos Part A 2013;55:143-51.

[37] De Carvalho NV, Pinho ST, Robinson P. Reducing the domain in the mechanical analysis of periodic structures, with application to woven composites. Compos Sci Technol 2011;71(7):969-79.

[38] Potter E, Pinho ST, Robinson P, Iannucci L, McMillan AJ. Mesh generation and geometrical modeling of 3D woven composites with variable tow cross-sections. Comput Mater Sci 2012;51(1):103-11.

[39] De Carvalho NV, Pinho ST, Robinson P. Numerical modelling of woven composites: biaxial loading. Compos Part A 2012;43(8):1326-37.

[40] Stig F, Hallström S. A modelling framework for composites containing 3D reinforcement. Compos Struct 2012;94(9):2895-901.

[41] Chiaruttini V, Riolo V, Feyel F. Advanced remeshing techniques for complex 3D crack propagation. In: 13th International conference on fracture, Beijing, China. p. 547-55.

[42] Charrier JS, Carrère N, Laurin F, Bretheau T, Goncalves-Novo E, Mahdi S. Proposition of 3D progressive failure approach and validation on tests cases. In: Proceedings of ECCM14 - 14th European conference on composite materials, Budapest, Hungary.

[43] Alfonso L, Uguen A, Badulescu C, Cognard JY, Bonnemains T, Lolive E, et al. Determination of the 3D failure envelope of a composite based on a modified Arcan test device. Compos Struct 2015;131:585-93.

[44] Rocher JE, Allaoui S, Hivet G, Gillibert J, Blond E. Experimental characterization and modeling of GF/PP commingled yarns tensile behavior. J Compos Mater 2015;49(21):2609-24.

[45] Benzeggagh M, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos Sci Technol 1996;56(4):439-49.

[46] Soden PD, Hinton MJ, Kaddour AS. Lamina properties, lay-up configurations and loading conditions for a range of fibre-reinforced composite laminates. Compos Sci Technol 1998;58:1011-22.

[47] García IG, Mantič V, Blázquez París F. Transverse crack onset and growth in cross-ply [0/90]s laminates under tension. Application of a coupled stress and energy criterion. Int J Sol Struct 2014;51(23-24):3844-56.

[48] García IG, Carter BJ, Ingraffea AR, Mantič V. A numerical study of transverse cracking in cross ply laminates by 3D finite fracture mechanics. Compos Part B 2016;95:475-87.