**3. Progressive Damage Model of Braided Composite**

The fiber bundle of textile composites is generally considered as a transversely isotropic unidirectional lamina in numerical simulation [25–27]. In this part, a progressive damage model of the unidirectional lamina is formulated in terms of damage initiation and damage evolution. This research aims to investigate the internal damage initiation and propagation of a notched specimen, therefore material failure (element deletion) is not introduced in this damage model.

### *3.1. Damage Initiation*

A three-dimensional failure criterion for the fiber bundle was adopted based on Hashin's [28] and Hou's [29] criteria and was incorporated with continuum damage laws. Four distinct failure modes are considered: fiber tensile failure, fiber compression failure, matrix tension failure, and matrix compression failure. The damage initiation criteria are formulated below.

For fiber tension failure (*σ*<sup>11</sup> > 0),

$$f\_{1t} = \left(\frac{\sigma\_{11}}{S\_{1t}}\right)^2 + a\left(\frac{\sigma\_{12} + \sigma\_{31}}{S\_{12}}\right) \ge 1\tag{1}$$

For fiber compression failure (*σ*<sup>11</sup> < 0),

$$f\_{1c} = \left(\frac{\sigma\_{11}}{\mathcal{S}\_{1c}}\right)^2 \ge 1 \tag{2}$$

For matrix tension failure (*σ*<sup>22</sup> > 0),

$$f\_{2t} = \left(\frac{\sigma\_{22}}{S\_{2t}}\right)^2 + \left(\frac{\sigma\_{12}}{S\_{12}}\right)^2 + \left(\frac{\sigma\_{23}}{S\_{23}}\right)^2 \ge 1\tag{3}$$

For matrix compression failure (*σ*<sup>22</sup> < 0),

$$f\_{2\varepsilon} = \frac{1}{4} \left(\frac{-\sigma\_{22}}{S\_{12}}\right)^2 + \frac{\sigma\_{22}}{S\_{2\varepsilon}} \left(\left(\frac{S\_{2\varepsilon}}{2S\_{12}}\right)^2 - 1\right) + \left(\frac{\sigma\_{12}}{S\_{23}}\right)^2 \ge 1\tag{4}$$

where *f* <sup>1</sup>*<sup>t</sup>* , *f* <sup>1</sup>*<sup>c</sup>* , *f* <sup>2</sup>*<sup>t</sup>* , and *f* <sup>2</sup>*<sup>c</sup>* are failure indices corresponding to each damage mode, respectively. The first subscripts, 1, 2, and 3, indicate the fiber axial direction, in-plane transverse direction, and out-of-plane direction, respectively. When the stress state of an element make one of the four failure indices larger than 1, the corresponding damage mode will be initiated in this element and constitutive law entering the stage of damage evolution. *S*1*<sup>t</sup>* , *S*1*<sup>c</sup>* , *S*2*<sup>t</sup>* , *S*2*<sup>c</sup>* , *S*12, and *S*<sup>23</sup> are axial tensile strength, axial compressive strength, transverse tensile strength, transverse compressive strength, longitudinal shear strength, and transverse shear strength of the fiber bundle. *α* is the shear failure coefficient which plays an important role in failure prediction of textile composite. The previous research conducted by Zhang et al. [18] indicated that the coefficient *α* has an obvious impact on the global stress–strain response and mainly on the failure prediction. The value of *α* was determined to be 0.06 for fiber bundles of triaxially-braided composite through correlation with an experimental ultimate strength of straight-sided coupon specimens [18,22]. The same value of parameter *α* is adopted in the present study in consideration that the studied materials are totally the same.

#### *3.2. Damage Evolution*

For the damage evolution behavior, the Murakami–Ohno [30] damage theory is adopted to predict the post-peak softening, and the crack band model developed by Bazant and Oh [30] is also employed to mitigate the mesh size dependency of the proposed mesoscale model in this study. A characteristic element length is introduced into damage evolution expression, aiming to dissipate the constant energy release rate per unit area in the solid element [31,32], and the element dissipated energy can be expressed as

$$\mathbf{G}\_{f,\mathbf{I}} = \frac{1}{2} \sigma\_{\mathbf{eq}}^f \mathbf{e}\_{\mathbf{eq}}^f l\_{\mathbf{c}} \tag{5}$$

where *l<sup>c</sup>* is the characteristic length of element, which calculates by extracting the cubic root of the volume of each element; *Gf*,*<sup>I</sup>* is fracture energy of fiber bundles corresponding to the specific damage mode *I*. The values for the fracture energies of axial and bias fiber bundles used in this study (listed in Table 2) were cited from Li et al. [17]. *σ f eq* and *ε f eq* are the equivalent peak stress and equivalent failure strain, respectively. The evolution of each damage variable is governed by an equivalent displacement expressed by the following equation.

$$d\_{I} = \frac{\delta\_{I,eq}^{f} \left(\delta\_{I,eq} - \delta\_{I,eq}^{0}\right)}{\delta\_{I,eq} \left(\delta\_{I,eq}^{f} - \delta\_{I,eq}^{0}\right)}, \text{ I} = ft, \text{ } f \text{c\\_mt\\_mc} \tag{6}$$

where *δ f <sup>I</sup>*,*eq* is the fully damaged equivalent displacement of the corresponding failure mode and *δ* 0 *I*,*eq* is the equivalent displacement at which the failure criterion is satisfied. For a certain failure mode, the equivalent displacement used in the initiation criteria is expressed in terms of the components corresponding to the effective stress components. Detailed algorithm equations of equivalent displacement and stress for each failure mode can be found in Zhang and coworkers' work [18]. As material parameters of fiber bundles, *δ* 0 *<sup>I</sup>*,*eq* and *δ f <sup>I</sup>*,*eq* be computed by the following equations.

$$\delta\_{\rm I,eq}^f = \frac{2G\_f}{\sigma\_{\rm I,eq}^0} \tag{7}$$

$$
\delta\_{\mathbf{I}, \mathbf{e}\boldsymbol{q}}^{0} = \frac{\delta\_{\mathbf{I}, \mathbf{e}\boldsymbol{q}}}{\sqrt{f\_{\mathbf{I}}}} \tag{8}
$$

**Table 2.** Fracture energies of the fiber bundle.


Here, *σ* 0 *<sup>I</sup>*,*eq* denotes the equivalence stress when each kind of damage criteria is satisfied. Meanwhile, the value of *f<sup>I</sup>* can be obtained from Equations (1)–(4).

To simulate the softening process of damage element, a second-order symmetric tensor is used to describe the damage state. The corresponding damaged compliance matrix *S*(*d*) is obtained as Equation (9), and the damaged stiffness matrix *C*(*d*) is the inverse of *S*(*d*).

$$S(d) = \begin{bmatrix} \frac{1}{d\_f E\_{11}} & -\frac{\nu\_{21}}{E\_{22}} & -\frac{\nu\_{31}}{E\_{33}} & \text{zero} \\ -\frac{\nu\_{12}}{E\_{11}} & \frac{1}{d\_m E\_{22}} & -\frac{\nu\_{32}}{E\_{33}} \\ -\frac{\nu\_{13}}{E\_{11}} & -\frac{\nu\_{23}}{E\_{22}} & \frac{1}{E\_{33}} \\ & & & \frac{1}{d\_f d\_m G\_{12}} \\ & & & & \frac{1}{d\_m G\_{23}} \\ & & & & & \frac{1}{d\_f G\_{13}} \end{bmatrix} \tag{9}$$

where *d<sup>f</sup>* and *d<sup>m</sup>* are global damage variables associated with fiber and matrix failure, which are introduced to control the degree of stiffness degeneration of damaged elements, and also satisfy the following equations, respectively.

$$d\_f = 1 - \min\left\{ \max\left(d\_{f\iota}, d\_{f\iota}\right), \gamma\_f \right\} \tag{10}$$

$$d\_m = 1 - \min\{\max(d\_{mt}, d\_{mc}), \gamma\_m\} \tag{11}$$

Here *γ<sup>f</sup>* and *γ<sup>m</sup>* are defined as the damage thresholds of global fiber and matrix damage, respectively. Numerically, the nonzero constants *γ<sup>f</sup>* and *γ<sup>m</sup>* address the singularity issue, and physically represent the effective resultant resistance of the homogeneous damaged elements. The effect of *γ<sup>f</sup>* and *γ<sup>m</sup>* on the global stress–strain responses will be discussed in a later section. Also, the Duvaut

and Lions regularization model [33] is applied to promote the numerical computation and smooth the stiffness degradation process.

### *3.3. Cohesive Element Model for Interface*

Interface is the bridge between the fiber bundle and matrix, which determines how stresses are transferred. The damage status of interface influences significantly the damage initiation and propagation of the composite material [34]. Littell [13] indicates that failure in the transverse tests was a result of edge delamination which occurred quickly and propagated along the bias fibers; and in the axial tensile direction, the subsurface delamination caused the nonlinearities in the global stress–strain response curves. In this study, the tow-to-tow interface is simulated by using the cohesive zone modeling approach, which has been embedded into ABAQUS as an optional element type. The responses of cohesive elements are governed by a typical bilinear traction–separation law, and a quadratic nominal stress criterion is used to describe interfacial damage initiation [32,34,35]. Besides, a power law criterion is adopted, which claims that failure under mixed-mode conditions is governed by a second-order power law interacting of the energies required to cause failure in the individual (normal and two shear) modes. The quadratic nominal stress criterion for damage initiation and a power law criterion for failure are represented in Equations (12) and (13).

$$
\left(\frac{\langle t\_n \rangle}{t\_n^0}\right)^2 + \left(\frac{\langle t\_s \rangle}{t\_s^0}\right)^2 + \left(\frac{\langle t\_l \rangle}{t\_l^0}\right)^2 = 1 \tag{12}
$$

$$\left(\frac{\text{G}\_{\text{fl}}}{\text{G}\_{\text{h}}^{\text{c}}}\right)^{2} + \left(\frac{\text{G}\_{\text{s}}}{\text{G}\_{\text{s}}^{\text{c}}}\right)^{2} + \left(\frac{\text{G}\_{\text{l}}}{\text{G}\_{\text{t}}^{\text{c}}}\right)^{2} = 1\tag{13}$$

where *t<sup>n</sup>* denotes the traction normal stress, and *t<sup>s</sup>* and *t<sup>t</sup>* denote shear stresses. The Macaulay brackets are used to signify that a pure compressive deformation or stress state does not initiate damage. *t* 0 *n* , *t* 0 *s* , and *t* 0 *t* represent the interface strength in normal and two shear directions. Similarly, *Gn*, *G<sup>s</sup>* , and *G<sup>t</sup>* refer to the work done by the traction and its conjugate relative displacement in the normal, first, and second shear directions, respectively; and *G c n* , *G c s* , and *G c t* are critical fracture energies required to cause failure in each of the three directions. Table 3 presents the interface properties, and the values of interface strengths and fracture toughness, which are cited from Zhang et al. [18]. Detailed formulations of the mixed-mode cohesive zone model can be found in the ABAQUS user's manual.

**Table 3.** Strengths and fracture toughness of cohesive elements.


#### **4. Finite Element Model Development**

The mesoscale finite element (FE) model simulates explicitly the fiber bundles of a braided composite structure and defines locally the realistic local fiber volume ratios and bundle orientations of the impregnated bundles. The advantage of the mesoscale model is its ability to analyze the local damage and failure of each component implemented individually through specific failure models for the various constituents and to predict the response of each constituent and their contribution to the global behavior. In this work, a mesoscale finite element model is introduced to study the internal damage and failure mechanism of the notched tensile specimen for 0◦/± 60◦ triaxially-braided composite.

Based on the geometric parameters identified by Zhang et al. [18], a FE model of a single unit cell, which can represent the composite's geometric features, was generated using 8-node solid element. As shown in Figure 2a, the mesh of a unit cell was constructed through TexGEN software by keying the dimensions of the unit cell and fiber bundles. In Figure 2, the axial fiber bundle, +60◦ and −60 ◦

bias fiber bundles and matrix elements, which fill the space between fiber bundles to form plate, are represented by light blue, dark blue, yellow, and green colors, respectively. Cohesive element layers (colored red), which have the same in-plane size as brick elements but zero thickness, are manually inserted between each two fiber bundles, and fiber bundle and matrix, see Figure 2a. It should be pointed out that the pure matrix is modeled as an elastic perfectly-plastic material. It is assumed that the pure matrix elements will not fail before the fracture of the specimen, due to the much larger failure strain of matrix than that of the fiber bundle. The resultant mesh may have different fiber bundle volume values than the real material, so the fiber volume ratio in each fiber bundle is adjusted to match the real fiber volume content. As identified by Zhang et al. [18], the realistic fiber volume ratios for axial and bias fiber bundles are 77% and 74.5%, respectively. The resultant fiber volume ratio in the present FE model is 86% for the axial fiber bundles and 69% for the bias fiber bundles. Mechanical properties of fiber bundles are listed in Table 4 referring to Zhang's work [18] for the same materials. and − 60 ° bias fiber bundles and matrix elements, which fill the space between fiber bundles to form plate, are represented by light blue, dark blue, yellow, and green colors, respectively. Cohesive element layers (colored red), which have the same in-plane size as brick elements but zero thickness, are manually inserted between each two fiber bundles, and fiber bundle and matrix, see Figure 2a. It should be pointed out that the pure matrix is modeled as an elastic perfectly-plastic material. It is assumed that the pure matrix elements will not fail before the fracture of the specimen, due to the much larger failure strain of matrix than that of the fiber bundle. The resultant mesh may have different fiber bundle volume values than the real material, so the fiber volume ratio in each fiber bundle is adjusted to match the real fiber volume content. As identified by Zhang et al. [18], the realistic fiber volume ratios for axial and bias fiber bundles are 77% and 74.5%, respectively. The resultant fiber volume ratio in the present FE model is 86% for the axial fiber bundles and 69% for the bias fiber bundles. Mechanical properties of fiber bundles are listed in Table 4 referring to Zhang's work [18] for the same materials.

**Figure 2.** (**a**) Composition of a unit cell finite element model. (**b**) Finite element mesh of axial tension notched model. **Figure 2.** (**a**) Composition of a unit cell finite element model. (**b**) Finite element mesh of axial tension notched model.

**Table 4.** Mechanical properties of axial and bias fiber tows.


Figure 2b shows the detailed dimensions of the FE model for the notched specimen. The FE model of the notched tension specimen consists of 208,000 linear brick elements (C3D8R) and 29,412 eight-node three-dimensional cohesive elements (COH3D8), with four unit cells through the width and length direction, respectively. Two whole unit cells (37 mm) are assembled through the width direction of the notched gauge for the axial tension model, which is consistent with the experimental Figure 2b shows the detailed dimensions of the FE model for the notched specimen. The FE model of the notched tension specimen consists of 208,000 linear brick elements (C3D8R) and 29,412 eight-node three-dimensional cohesive elements (COH3D8), with four unit cells through the width and length direction, respectively. Two whole unit cells (37 mm) are assembled through the width direction of the notched gauge for the axial tension model, which is consistent with the experimental specimen. The length of the finite element model is 20.68 mm with four unit cells for axial tension, which are

*S*12 = *S*13 = *S*23 (MPa) 80.60 78.53

specimen. The length of the finite element model is 20.68 mm with four unit cells for axial tension,

shorter than the gauge length of experimental specimens fabricated by Kohlman [14] (Figure 1b). These reasonable simplifications of model size are intended to reduce the computational time, in consideration of the minor effect of the remote sections. The numerical models are solved by ABAQUS Explicit using a 24-core workstation and each job costs ~20 h. By evaluating the various energies generated during the computation process, the accumulated kinetic energy was always less than 1% of the internal energy of the model, which ensures a quasi-static loading status of the problem. time, in consideration of the minor effect of the remote sections. The numerical models are solved by ABAQUS Explicit using a 24-core workstation and each job costs ~20 h. By evaluating the various energies generated during the computation process, the accumulated kinetic energy was always less than 1% of the internal energy of the model, which ensures a quasi-static loading status of the problem. **5. Results and Discussion** 

*Materials* **2019**, *12*, x FOR PEER REVIEW 9 of 19

(Figure 1b). These reasonable simplifications of model size are intended to reduce the computational

#### **5. Results and Discussion** In this section, the damage initiation and propagation behavior of the notched tension specimen

In this section, the damage initiation and propagation behavior of the notched tension specimen modeled using the proposed mesoscale FE scheme is correlated with experimental results. The effect of specimen geometry on the internal damage evolution behavior and the effective strength of this notched specimen are further discussed through numerical parametric studies. These results show the applicability and reliability of the mesoscale FE model for failure study of braided composites. modeled using the proposed mesoscale FE scheme is correlated with experimental results. The effect of specimen geometry on the internal damage evolution behavior and the effective strength of this notched specimen are further discussed through numerical parametric studies. These results show the applicability and reliability of the mesoscale FE model for failure study of braided composites. *5.1. Model Correlation* 

#### *5.1. Model Correlation* In our previous works [15,18,21,22], the failure behavior of straight-side coupon specimens of

In our previous works [15,18,21,22], the failure behavior of straight-side coupon specimens of 2DTBC has been intensively studied. Figure 3 shows the comparison of the experimental and numerical predicted results for the single-layer straight-side coupon specimen under axial tension. The results indicate that the mesoscale FE model can not only predict well the global stress–strain curve, but also the strain distributions which are highly sensitive to the local braided architecture. 2DTBC has been intensively studied. Figure 3 shows the comparison of the experimental and numerical predicted results for the single-layer straight-side coupon specimen under axial tension. The results indicate that the mesoscale FE model can not only predict well the global stress–strain curve, but also the strain distributions which are highly sensitive to the local braided architecture.

**Figure 3.** Comparison of the experimental and numerical predicted results for straight-side coupon specimen. (**a**) Global stress–strain curves. (**b**) Distribution of the axial and in-plane shear strain at global strain level of 2.0%. **Figure 3.** Comparison of the experimental and numerical predicted results for straight-side coupon specimen. (**a**) Global stress–strain curves. (**b**) Distribution of the axial and in-plane shear strain at global strain level of 2.0%.

The applicability of the mesoscale FE model is further demonstrated through the model validation for the notched specimen. For both the experimental characterization and numerical simulations, the effective strength of the tensile specimen is determined as the ratio of total reaction force at the loading section against the cross-section area at the gauge section (two unit cells wide for both notched and straight-sided coupon specimen). Table 5 compares the numerical predicted and experimental measured effective strength of triaxially-braided composite under axial tension, for both straight-sided coupon and notched specimens. The measured strength values of the straightsided coupon and notched specimen are referred to as Kohlman [14]. It is evident from Table 5 that the mesoscale FE model predicted effective strength values are in good agreement with the experimental values, suggesting the accuracy of this modeling scheme. The effective strength of the notched specimen tends to be lower than the straight-sided coupon specimen, which is due to the presence of stress concentration at the notched section. The applicability of the mesoscale FE model is further demonstrated through the model validation for the notched specimen. For both the experimental characterization and numerical simulations, the effective strength of the tensile specimen is determined as the ratio of total reaction force at the loading section against the cross-section area at the gauge section (two unit cells wide for both notched and straight-sided coupon specimen). Table 5 compares the numerical predicted and experimental measured effective strength of triaxially-braided composite under axial tension, for both straight-sided coupon and notched specimens. The measured strength values of the straight-sided coupon and notched specimen are referred to as Kohlman [14]. It is evident from Table 5 that the mesoscale FE model predicted effective strength values are in good agreement with the experimental values, suggesting the accuracy of this modeling scheme. The effective strength of the notched specimen tends to be lower than the straight-sided coupon specimen, which is due to the presence of stress concentration at the notched section.

acceptable.


**Table 5.** Comparison of simulation predicted and experimental measured the effective tensile strength of the triaxially-braided composite. *Materials* **2019**, *12*, x FOR PEER REVIEW 10 of 19 **Table 5.** Comparison of simulation predicted and experimental measured the effective tensile

> \* Refer to Kohlman et al. [14]. Double-notched 751 765 \* \* Refer to Kohlman et al. [14].

To further demonstrate the reliability of the meso-FE model, the numerically predicted strain distributions are compared with DIC results. As seen in Figure 4, the simulation results match well with the local strain distribution of the notched specimen and additionally capturing the local strain concentration around the notch area and the diagonal distribution along the bias fiber bundles. The axial strain (*εyy*) and shear strain (*εxy*) distributions show higher strain at the pure matrix zone between bias fiber bundles, which coincides with the damage propagation paths of fiber bundle damage. The damage of fiber bundles will induce the decrease of fiber's capability in carrying the load and as a result, more load will be transferred to the matrix material, thereby causing the strain concentration in this area (red and blue in Figure 4). On the other hand, the shear strain shows more obvious strain concentration effect at the notched section. Unlike the traditional isotropic material where strain concentration effect propagates in a circular shape, the shear strain concentration effect of this braided composite propagates along the fiber bundle directions. The mesoscale FE model simulations capture the strain distribution and its magnitude accurately for the axial tension of the notched specimen. The result also reflects the advantages of mesoscale finite element model in studying the internal damage of textile composite that is difficult to achieve from experiments. To further demonstrate the reliability of the meso-FE model, the numerically predicted strain distributions are compared with DIC results. As seen in Figure 4, the simulation results match well with the local strain distribution of the notched specimen and additionally capturing the local strain concentration around the notch area and the diagonal distribution along the bias fiber bundles. The axial strain (*εyy*) and shear strain (*εxy*) distributions show higher strain at the pure matrix zone between bias fiber bundles, which coincides with the damage propagation paths of fiber bundle damage. The damage of fiber bundles will induce the decrease of fiber's capability in carrying the load and as a result, more load will be transferred to the matrix material, thereby causing the strain concentration in this area (red and blue in Figure 4). On the other hand, the shear strain shows more obvious strain concentration effect at the notched section. Unlike the traditional isotropic material where strain concentration effect propagates in a circular shape, the shear strain concentration effect of this braided composite propagates along the fiber bundle directions. The mesoscale FE model simulations capture the strain distribution and its magnitude accurately for the axial tension of the notched specimen. The result also reflects the advantages of mesoscale finite element model in studying the internal damage of textile composite that is difficult to achieve from experiments.

**Figure 4.** Comparison of numerical predicted and experimentally measured strain contours of the notched specimen under axial tension. **Figure 4.** Comparison of numerical predicted and experimentally measured strain contours of the notched specimen under axial tension.

Furthermore, another reason for the localized strain concentration effect is because of the lower local fiber volume ratio at the corresponding regions. The inconspicuous discrepancy between numerically predicted and experimental axial strain contours is very likely results from the idealization of the geometry in the FE model. However, the manufacturing defects like axial fiber bundle undulation are inevitable in actual specimens. The cutting process may also introduce unexpected fiber damage, which could lead to more significant strain concentration area near the notch. Regardless of these factors, the capability and accuracy of this 3D mesoscale model are highly Furthermore, another reason for the localized strain concentration effect is because of the lower local fiber volume ratio at the corresponding regions. The inconspicuous discrepancy between numerically predicted and experimental axial strain contours is very likely results from the idealization of the geometry in the FE model. However, the manufacturing defects like axial fiber bundle undulation are inevitable in actual specimens. The cutting process may also introduce unexpected fiber damage, which could lead to more significant strain concentration area near the notch. Regardless of these factors, the capability and accuracy of this 3D mesoscale model are highly acceptable.

As mentioned before in the introduction section, both Littell [13] and Zhang et al. [15] observed

As mentioned before in the introduction section, both Littell [13] and Zhang et al. [15] observed out-of-plane warping behavior along the free edges in their tensile tests using the straight-side coupon specimen, which would lead to the premature damage initiation from the free edges. Kolhman et al. [14] also discovered the same phenomena in the notched specimen. Figure 5 compares the out-of-plane displacement (U3) contours from simulation and experiment for the notched tension test. The out-of-plane displacement is formed due to the tension–torsion coupling resulting from the termination of fiber bundles at the free edges. In the simulation results, the specimen shows the antisymmetric distribution of out-of-plane warping at the four corners of the notched specimen, which is related to the antisymmetric mesoscopic structure of this material. In addition, no obvious out-of-plane deformation was observed in the gauge section of the specimen indicates that the failure of the specimen is not affected by the free-edge effect. However, the out-of-plane deformation still consumes partially the energy from external loading, leading to a relatively lower measured modulus of the specimen. Thus, the notched specimen may not be suitable for modulus measurement. This is one of the contents which need to be studied in the future. *Materials* **2019**, *12*, x FOR PEER REVIEW 11 of 19 coupon specimen, which would lead to the premature damage initiation from the free edges. Kolhman et al. [14] also discovered the same phenomena in the notched specimen. Figure 5 compares the out-of-plane displacement (U3) contours from simulation and experiment for the notched tension test. The out-of-plane displacement is formed due to the tension–torsion coupling resulting from the termination of fiber bundles at the free edges. In the simulation results, the specimen shows the antisymmetric distribution of out-of-plane warping at the four corners of the notched specimen, which is related to the antisymmetric mesoscopic structure of this material. In addition, no obvious out-of-plane deformation was observed in the gauge section of the specimen indicates that the failure of the specimen is not affected by the free-edge effect. However, the out-of-plane deformation still consumes partially the energy from external loading, leading to a relatively lower measured modulus of the specimen. Thus, the notched specimen may not be suitable for modulus measurement. This is one of the contents which need to be studied in the future. coupon specimen, which would lead to the premature damage initiation from the free edges. Kolhman et al. [14] also discovered the same phenomena in the notched specimen. Figure 5 compares the out-of-plane displacement (U3) contours from simulation and experiment for the notched tension test. The out-of-plane displacement is formed due to the tension–torsion coupling resulting from the termination of fiber bundles at the free edges. In the simulation results, the specimen shows the antisymmetric distribution of out-of-plane warping at the four corners of the notched specimen, which is related to the antisymmetric mesoscopic structure of this material. In addition, no obvious out-of-plane deformation was observed in the gauge section of the specimen indicates that the failure of the specimen is not affected by the free-edge effect. However, the out-of-plane deformation still consumes partially the energy from external loading, leading to a relatively lower measured modulus of the specimen. Thus, the notched specimen may not be suitable for modulus measurement. This is one of the contents which need to be studied in the future.

*Materials* **2019**, *12*, x FOR PEER REVIEW 11 of 19

**Figure 5.** Comparison of numerical predicted and experimental measured out-of-plane displacement contours of notched tension specimen. **Figure 5.** Comparison of numerical predicted and experimental measured out-of-plane displacement contours of notched tension specimen. **Figure 5.** Comparison of numerical predicted and experimental measured out-of-plane displacement contours of notched tension specimen.

#### *5.2. Progressive Damage Analysis 5.2. Progressive Damage Analysis*

indicates that the element is totally damaged.

*5.2. Progressive Damage Analysis*  The progressive damage process is then studied using the correlated mesoscale FE model, as shown in Figure 6. The status of damage corresponds to the value of the particular damage variable ranging from 0 to 1, where a value of 0 indicates that damage has not occurred yet while a value of 1 The progressive damage process is then studied using the correlated mesoscale FE model, as shown in Figure 6. The status of damage corresponds to the value of the particular damage variable ranging from 0 to 1, where a value of 0 indicates that damage has not occurred yet while a value of 1 indicates that the element is totally damaged. The progressive damage process is then studied using the correlated mesoscale FE model, as shown in Figure 6. The status of damage corresponds to the value of the particular damage variable ranging from 0 to 1, where a value of 0 indicates that damage has not occurred yet while a value of 1 indicates that the element is totally damaged.


**Figure 6.** Comparison of numerical predicted damage development of axial tension. **Figure 6.** Comparison of numerical predicted damage development of axial tension. **Figure 6.** Comparison of numerical predicted damage development of axial tension.

The matrix damage initiated at the global strain level of ~0.5% and spread to the whole central region at the global strain level of 0.80%. The distribution of internal damage inside the bundles

The matrix damage initiated at the global strain level of ~0.5% and spread to the whole central region at the global strain level of 0.80%. The distribution of internal damage inside the bundles reveals that there is no fiber tension damage occurring at these strain levels while matrix tension

The matrix damage initiated at the global strain level of ~0.5% and spread to the whole central region at the global strain level of 0.80%. The distribution of internal damage inside the bundles reveals that there is no fiber tension damage occurring at these strain levels while matrix tension damage (matrix cracking in fiber bundles) is observed to start from the notched area. The matrix tension damage propagates along the notches and the bias fiber bundles; while in the meantime, in the central region matrix tension damage occurs mainly where bias fiber bundles intersect. The interfaces are found to be intact at the first two status of Figure 6 as there is no damage in cohesive elements. This behavior is similar to the observations during the test of straight-sided coupon specimen [13]. *Materials* **2019**, *12*, x FOR PEER REVIEW 12 of 19 damage (matrix cracking in fiber bundles) is observed to start from the notched area. The matrix tension damage propagates along the notches and the bias fiber bundles; while in the meantime, in the central region matrix tension damage occurs mainly where bias fiber bundles intersect. The interfaces are found to be intact at the first two status of Figure 6 as there is no damage in cohesive elements. This behavior is similar to the observations during the test of straight-sided coupon specimen [13].

At the global strain level of 1.8%, unloading is identified followed by fiber tensile damage of axial fiber bundles and almost all elements of bias fiber bundles are dominated by matrix tension damage. The fiber tension damage initiates at the notched area same as the matrix tension damage. It is also found that interface failure appears only near the notches due to localized shear stress concentration. The inert of interface for all other areas indicates the excellent interface properties of this triaxially-braided composite, which is consistent with the previous experimental examinations where delamination is rarely observed during axial tension tests of this material [14]. At the global strain level of 1.8%, unloading is identified followed by fiber tensile damage of axial fiber bundles and almost all elements of bias fiber bundles are dominated by matrix tension damage. The fiber tension damage initiates at the notched area same as the matrix tension damage. It is also found that interface failure appears only near the notches due to localized shear stress concentration. The inert of interface for all other areas indicates the excellent interface properties of this triaxially-braided composite, which is consistent with the previous experimental examinations where delamination is rarely observed during axial tension tests of this material [14].

Figure 7 shows the numerically predicted stress distributions of axial and biaxial tows before the specimen failure. The stresses on both sides of the notches are not symmetrical because of the asymmetry geometry feature. As seen from the stress contours, most of the applied loads are afforded by the axial bundles across the notched sections. Similar to the strain distribution contours in Figure 4, the shear stress concentration at the notched area is more significant than the normal stress. The shear stress concentration will then result in shear tension-dominated failure of axial fiber bundles at the notched area and tension dominated failure of axial fiber bundles at the central area, corresponding to the inclined (along bias fiber bundle direction) failure pattern of the two axial fiber bundles near the notches and horizontal failure pattern of fiber bundles at other region (fiber tension damage at global strain level 1.8% as shown in Figure 6). Similarly, the stress distribution of bias fiber bundles also explains the initiation and propagation behavior of matrix tension damage behavior at global strain levels 0.5% and 0.8% (Figure 6). Figure 7 shows the numerically predicted stress distributions of axial and biaxial tows before the specimen failure. The stresses on both sides of the notches are not symmetrical because of the asymmetry geometry feature. As seen from the stress contours, most of the applied loads are afforded by the axial bundles across the notched sections. Similar to the strain distribution contours in Figure 4, the shear stress concentration at the notched area is more significant than the normal stress. The shear stress concentration will then result in shear tension-dominated failure of axial fiber bundles at the notched area and tension dominated failure of axial fiber bundles at the central area, corresponding to the inclined (along bias fiber bundle direction) failure pattern of the two axial fiber bundles near the notches and horizontal failure pattern of fiber bundles at other region (fiber tension damage at global strain level 1.8% as shown in Figure 6). Similarly, the stress distribution of bias fiber bundles also explains the initiation and propagation behavior of matrix tension damage behavior at global strain levels 0.5% and 0.8% (Figure 6).

**Figure 7.** Numerical predicted stress contours before the instant of failure under axial tension. **Figure 7.** Numerical predicted stress contours before the instant of failure under axial tension.

Overall, matrix tension damage appears firstly around the notched areas and propagates along the bias fiber bundles, resulting in a slight decrease of effective stiffness. With the increase of external loads, the fiber tension damage and delamination will initiate around the notches due to shear stress and strain concentration in this region. The notch-induced stress/strain concentration effect disturbs Overall, matrix tension damage appears firstly around the notched areas and propagates along the bias fiber bundles, resulting in a slight decrease of effective stiffness. With the increase of externalloads, the fiber tension damage and delamination will initiate around the notches due to shear stressand strain concentration in this region. The notch-induced stress/strain concentration effect disturbsonly a limited region of the specimen and the failure of the specimen is mainly due to normal tension

only a limited region of the specimen and the failure of the specimen is mainly due to normal tension stress-induced fiber tension damage in the axial fiber bundles. Thus, the measured tensile strength of

The damage areas at strain level 1.8% of Figure 6 are in good correlation with the strain

a notched specimen can be representative and can be a lower bound of the realistic strength.

stress-induced fiber tension damage in the axial fiber bundles. Thus, the measured tensile strength of a notched specimen can be representative and can be a lower bound of the realistic strength.

The damage areas at strain level 1.8% of Figure 6 are in good correlation with the strain concentration region of the experimental results reported by Kolhman et al. [14]. As observed by Kolhman, the highest strain appears near the notch tips, which mainly attribute to the fiber tension failure. The relatively high strain concentrated along bias fiber is a result of matrix damage within the fiber bundles. The results of this section demonstrate the capability of a mesoscale model in predicting the internal damage initiation and propagation behavior of textile composites using various specimen shapes.

#### *5.3. Numerical Parametric Study*

Numerical studies were carried out to further investigate the features of the proposed mesoscale FE model and how the specific material parameters contribute to the global response of the model. This is a difficult task mainly because of the tedious modeling process and the considerable computational quantity. Due to nonuniformity of the strain distribution for the entire specimen, smooth stress–strain curves can hardly be produced in the notched tension tests, which were used mainly for strength measurement in Kohlman's work [14]. For numerical comparison, in the meso-FE simulations, stress–strain curves are generated based on the method of "digital strain gauge" proposed by Littell [13]. The macroscopic effective stress for the analyzed section is computed by determining the summation of the reaction forces on the face of the loaded cross-section divided by the area of the cross-section between two notches. The effective strain is calculated by dividing the relative displacement along the loading direction between two nodes by the initial distance between these two nodes before loading, where the relative displacement is calculated as the midpoints on each end (through load direction) section of the gauge region.

Figure 8a shows the variation of global stress–strain responses with different damage thresholds (*γ<sup>f</sup>* and *γm*). The dashed line in the figure corresponds to the experimental measured effective strength of this notched specimen by Kohlman [14]. As mentioned previously, *γ<sup>f</sup>* and *γ<sup>m</sup>* are the damage thresholds of *d<sup>f</sup>* and *dm*, which are the global damage variables associated with fiber-dominated and matrix-dominated failure, respectively. *d<sup>f</sup>* and *d<sup>m</sup>* control the descent of element stiffness along with damage accumulation, while *γ<sup>f</sup>* and *γ<sup>m</sup>* determine the extent of stiffness degradation. As we can see in Figure 8a, there is a sudden force drop with damage accumulation when *γ<sup>f</sup>* is lower than 0.3 and *γ<sup>m</sup>* is lower than 0.5, which is not physically consistent with the experimental stress–strain response of this material where the slope of the stress–strain curve tends to degrade gradually. Numerically, the stress–strain curves become smooth when *γ<sup>f</sup>* and *γ<sup>m</sup>* are larger than 0.3 and 0.5, respectively. Combined with the progressive damage behavior of this specimen discussed in the previous sections, the effect of *γ<sup>f</sup>* and *γ<sup>m</sup>* is concluded as follows; a smaller value of *γ<sup>f</sup>* leads to premature unloading of the stress–strain curve and *γ<sup>m</sup>* impacts the extent of nonlinearity of the curve. This is because the external load is mainly carried by axial fiber bundles under axial tension and the fiber tension damage mainly presences in the axial bundles. Then a rapid degradation of stiffness caused by fiber damage will result in a loss of load-bearing capacity instantaneously in the local area and correspond to a sharp drop of the stress–strain curve. On the other hand, *γ<sup>m</sup>* may not affect much the global responses at the initial stage due to the anisotropic feature of fiber bundles. However, as the matrix damage accumulates and propagates, *γ<sup>m</sup>* becomes more significant and could result in unexpected slope change followed by the initiation of fiber damage (see Figure 8a *γ<sup>f</sup>* = 0.3 and *γ<sup>m</sup>* = 0.5). Also, the predicted effective strength is found to be sensitive to the parameters *γ<sup>f</sup>* and *γm*. In this work, the numerical predicted effective strength correlates with experimental results when *γ<sup>f</sup>* = 0.3 and *γ<sup>m</sup>* = 0.5.

**Figure 8.** Numerical predicted stress-stain responses of notched tension specimen with: (**a**) different damage thresholds and (**b**) different shear failure coefficients. **Figure 8.** Numerical predicted stress-stain responses of notched tension specimen with: (**a**) different damage thresholds and (**b**) different shear failure coefficients.

The shear failure coefficient *α* controls the contribution of shear stress on fiber tension damage. The results exhibited in Figure 8b are similar to the conclusion for the same material reported in Zhang [18]. It is found that the effect of *α* on effective strength is not apparent when the value of *α* is lower than 0.5. As expected, the notch induced shear stress/strain concentration produces slightly higher sensitivity of the global stress–strain responses to the shear coefficient. The results of parameter analysis also declare that the present constitutive model of fiber bundles needs further improvement to enhance accuracy and applicability. The shear failure coefficient *α* controls the contribution of shear stress on fiber tension damage. The results exhibited in Figure 8b are similar to the conclusion for the same material reported in Zhang [18]. It is found that the effect of *α* on effective strength is not apparent when the value of *α* is lower than 0.5. As expected, the notch induced shear stress/strain concentration produces slightly higher sensitivity of the global stress–strain responses to the shear coefficient. The results of parameter analysis also declare that the present constitutive model of fiber bundles needs further improvement to enhance accuracy and applicability.

#### *5.4. Effect of Specimen Shape 5.4. Effect of Specimen Shape*

The mesoscale FE model proposed in this study can also be used to assess the rationality of the designed specimen. The effect of geometrical characteristics of the notched specimen on the test results was not considered by Kohlmen [14], due to the enormous consumption of time and The mesoscale FE model proposed in this study can also be used to assess the rationality of the designed specimen. The effect of geometrical characteristics of the notched specimen on the test results was not considered by Kohlmen [14], due to the enormous consumption of time and availability of specimens. In this work, the correlated mesoscale FE model is employed to address this concern.

availability of specimens. In this work, the correlated mesoscale FE model is employed to address this concern. Referring to the 0°/± 60° braided architecture shown in Figure 2a, a unit cell of this braided composite can be discretized into four adjacent subcell regions depending on the presence of axial and braided tows or lack thereof. A subcell-based modeling approach has been used by many researchers [12,36] to investigate the static and impact behavior of this triaxially-braided composite. Subcells A and C contain both axial and bias bundles while subcells B and D contain only bias bundles, so the distinction of fiber volume in each subcell leads to the different local mechanical properties in a unit cell. As from the top view and front-side view of the unit cell (Figure 2a), subcells C and D are antisymmetric against subcells A and B in the thickness direction. The section studies the relationship of different subcell distributions in the gauge region with damage evolution behavior and the ultimate strength property of the notched specimen. Figure 9 compares the geometry and damage contours of different notched specimens, where the subcells adjacent to the notches of the specimens are inequitable. The label above each specimen represents the two subcells which are closest to the two notches. For instance, "A-D" indicates that subcell A and D are located in the ends Referring to the 0◦/± 60◦ braided architecture shown in Figure 2a, a unit cell of this braided composite can be discretized into four adjacent subcell regions depending on the presence of axial and braided tows or lack thereof. A subcell-based modeling approach has been used by many researchers [12,36] to investigate the static and impact behavior of this triaxially-braided composite. Subcells A and C contain both axial and bias bundles while subcells B and D contain only bias bundles, so the distinction of fiber volume in each subcell leads to the different local mechanical properties in a unit cell. As from the top view and front-side view of the unit cell (Figure 2a), subcells C and D are antisymmetric against subcells A and B in the thickness direction. The section studies the relationship of different subcell distributions in the gauge region with damage evolution behavior and the ultimate strength property of the notched specimen. Figure 9 compares the geometry and damage contours of different notched specimens, where the subcells adjacent to the notches of the specimens are inequitable. The label above each specimen represents the two subcells which are closest to the two notches. For instance, "A-D" indicates that subcell A and D are located in the ends of gauge region while "A/2-A/2" expresses that half of subcell A connects with the notches.

of gauge region while "A/2-A/2" expresses that half of subcell A connects with the notches.

**Figure 9.** Numerical predictions of matrix damage in fiber bundles for notched tension specimen with **Figure 9.** Numerical predictions of matrix damage in fiber bundles for notched tension specimen with different architectures across the gauge section.

different architectures across the gauge section.

As from Figure 9, the damages of all specimens initiate and concentrate at the notched area at a global strain level of 0.5%. For specimens "A-D" and "B-A", the damage initiates and concentrates at one of the two notched zones, and the damage initiation locations of the two specimens are asymmetrically suggesting the asymmetric fabric architecture of the two specimens. This damage behavior is related to the fabrics subcells layout (Figure 9), where subcells A and C consist of both axial and bias bundles and are supposed to be stronger in the axial direction than subcells B and D. Thus, when subcells A/C and B/D are subjected to the same extent of load, subcells B and D are likely to damage earlier. Similarly, specimens "B/2-B/2" and "A/2-A/2" show symmetric damage progression behavior, as they both have the same fabric subcells layout at the notched area. Besides, comparing the specimens "B/2-B/2" and "A/2-A/2", it is found that the presence of axial fiber bundles along the notches restrict the amount of matrix damage and its propagation in fiber bundles, especially obvious at the global strain level of 0.75% where the matrix damage area is much smaller than that in the other three specimens. However, matrix damage in fiber bundles has little influence on the axial stiffness of fiber bundles and shows the negligible impact to the ultimate strength of the notched specimen. Table 6 lists the predicted effective strength of the four different notched specimens shown in Figure 9. While all specimens have more or less the same effective strength, specimen "A/2-A/2" shows relatively lower strength attributing to the incomplete axial fiber bundles in the gauge region. In the numerical simulations, the crossing of the fiber damage through the width (*x-*direction) of a single fiber bundle corresponds to the fracture of the specimen. Thus, the presence of incomplete axial fiber bundles at both ends of the gauge region in specimen "A/2-A/2" results in an earlier failure of the specimen. It should be noted that the fiber volume ratio in the gauge region of these four kinds of specimens is the same to ensure the comparability. As from Figure 9, the damages of all specimens initiate and concentrate at the notched area at a global strain level of 0.5%. For specimens "A-D" and "B-A", the damage initiates and concentrates at one of the two notched zones, and the damage initiation locations of the two specimens are asymmetrically suggesting the asymmetric fabric architecture of the two specimens. This damage behavior is related to the fabrics subcells layout (Figure 9), where subcells A and C consist of both axial and bias bundles and are supposed to be stronger in the axial direction than subcells B and D. Thus, when subcells A/C and B/D are subjected to the same extent of load, subcells B and D are likely to damage earlier. Similarly, specimens "B/2-B/2" and "A/2-A/2" show symmetric damage progression behavior, as they both have the same fabric subcells layout at the notched area. Besides, comparing the specimens "B/2-B/2" and "A/2-A/2", it is found that the presence of axial fiber bundles along the notches restrict the amount of matrix damage and its propagation in fiber bundles, especially obvious at the global strain level of 0.75% where the matrix damage area is much smaller than that in the other three specimens. However, matrix damage in fiber bundles has little influence on the axial stiffness of fiber bundles and shows the negligible impact to the ultimate strength of the notched specimen. Table 6 lists the predicted effective strength of the four different notched specimens shown in Figure 9. While all specimens have more or less the same effective strength, specimen "A/2-A/2" shows relatively lower strength attributing to the incomplete axial fiber bundles in the gauge region. In the numerical simulations, the crossing of the fiber damage through the width (*x*-direction) of a single fiber bundle corresponds to the fracture of the specimen. Thus, the presence of incomplete axial fiber bundles at both ends of the gauge region in specimen "A/2-A/2" results in an earlier failure of the specimen. It should be noted that the fiber volume ratio in the gauge region of these four kinds of specimens is the same to ensure the comparability.

**Table 6.** The effective strength of notched specimen with different architecture across the gauge **Table 6.** The effective strength of notched specimen with different architecture across the gauge section.


To investigate the effect of notch geometry on the effective strength of the specimen numerically, specimens with three different notch sizes (dimension in *y*-direction as shown in Figure 10a) and five different widths of the gauge section (dimension in *x*-direction as shown in Figure 10b) were studied. To investigate the effect of notch geometry on the effective strength of the specimen numerically, specimens with three different notch sizes (dimension in *y-*direction as shown in Figure 10a) and five different widths of the gauge section (dimension in *x*-direction as shown in Figure 10b) were studied.

*Materials* **2019**, *12*, x FOR PEER REVIEW 16 of 19

**Figure 10.** Mesoscale models of notched specimens with different geometrical characteristics. (**a**) Different notch sizes. (**b**) Different widths of the gauge section (with complete axial tows between notches). **Figure 10.** Mesoscale models of notched specimens with different geometrical characteristics. (**a**) Different notch sizes. (**b**) Different widths of the gauge section (with complete axial tows between notches).

The numerically predicted effective strength results are also listed in Table 7, where the effective strength values are found to be insensitive to the height (*y*-direction value) of the notch size for the cases studied in this work. However, the axial ultimate strength for the specimen changes with the changing of widths for the notch, as shown in Table 7. The predicted values of specimens denoted by "1.5UC" and "2.5UC" are slightly higher than the other three models, which may attribute to the integrity of axial tows in the gauge region. On the other hand, axial fibers among gauge regions of the other three specimens are all incomplete resulting in earlier failure of the specimens than that of specimens with complete fiber bundles. The numerically predicted effective strength results are also listed in Table 7, where the effective strength values are found to be insensitive to the height (*y*-direction value) of the notch size for the cases studied in this work. However, the axial ultimate strength for the specimen changes with the changing of widths for the notch, as shown in Table 7. The predicted values of specimens denoted by "1.5UC" and "2.5UC" are slightly higher than the other three models, which may attribute to the integrity of axial tows in the gauge region. On the other hand, axial fibers among gauge regions of the other three specimens are all incomplete resulting in earlier failure of the specimens than that of specimens with complete fiber bundles.


**Table 7.** The predicted strength of notched specimens with different geometrical characteristics. **Table 7.** The predicted strength of notched specimens with different geometrical characteristics.

 3 UC 760 \* With complete axial tows between notches.

\* With complete axial tows between notches.

The fiber damage behavior of two different specimens ("1.5UC" and "2UC") is shown in Figure 11. As we can see, fiber damage starts in the specimen "2UC" at the global strain level of 1.4%. The damage will spread quickly across the axial fiber bundles of half-width and propagate into the central area. By contrast, fiber damage is not observed in the "1.5UC" specimens at the global strain level of 1.4%. Fiber damage at a global strain level of 1.8% was found to be less significant compared to that in the specimen "2UC". The results further confirm the necessity of including complete axial fiber The fiber damage behavior of two different specimens ("1.5UC" and "2UC") is shown in Figure 11. As we can see, fiber damage starts in the specimen "2UC" at the global strain level of 1.4%. The damage will spread quickly across the axial fiber bundles of half-width and propagate into the central area. By contrast, fiber damage is not observed in the "1.5UC" specimens at the global strain level of 1.4%. Fiber damage at a global strain level of 1.8% was found to be less significant compared to that in the specimen "2UC". The results further confirm the necessity of including complete axial fiber bundles in

bundles in the gauge section when preparing the notched specimens, which could avoid possible variation from cutting and provides more stable and reliable measured properties of the material.

the gauge section when preparing the notched specimens, which could avoid possible variation from cutting and provides more stable and reliable measured properties of the material. *Materials* **2019**, *12*, x FOR PEER REVIEW 17 of 19

**Figure 11.** Comparison of predicted damage in in-complete (**a**) and complete (**b**) axial fiber bundles in gauge region. **Figure 11.** Comparison of predicted damage in in-complete (**a**) and complete (**b**) axial fiber bundles in gauge region.

The appropriate design of the specimen is the prerequisite for determining accurately the mechanical properties of braided composites, which is a big challenge due to the lack of advanced test standards. It is an onerous and costly process to optimize the specimen through an experiment The appropriate design of the specimen is the prerequisite for determining accurately the mechanical properties of braided composites, which is a big challenge due to the lack of advanced test standards. It is an onerous and costly process to optimize the specimen through an experiment only, and the presented numerical method is an alternative and efficient tool for specimen design.

only, and the presented numerical method is an alternative and efficient tool for specimen design. Overall, the fabrics architecture across the gauge section has little influence on the effective strength of the notched sample. However, the simulation capability of the mesomechanical model in predicting the impacts of mesoscale geometric features on the global effective properties of the samples is valuable and advantageous. Overall, the fabrics architecture across the gauge section has little influence on the effective strength of the notched sample. However, the simulation capability of the mesomechanical model in predicting the impacts of mesoscale geometric features on the global effective properties of the samples is valuable and advantageous.
