Experimental and Numerical Analysis of Progressive Damage of SiC_f/SiC Composite under Three-Point Bending

Abstract

Satin SiCf/SiC composite has a wide range of applications; it is necessary to study its mechanical properties. The progressive failure of five-harness five-layer satin weave SiCf/SiC plate composites was explored experimentally and numerically in this research. The bending properties were derived and elucidated at ambient temperature through a three-point bending experiment. The generation and progression of damage was observed by CCD camera. For quantitative analysis of the strain field evolution, DIC (digital image correlation) was adopted, while the microscopic analyses were performed for the description of the derived failure markings. With the aid of ABAQUS/Explicit, the experiment was subjected to 3D finite element modeling for the reproduction of the material behavioral, where the VUMAT subroutine was used to implement a 3D-altered criterion of the Hashin damage initiation and the progression law of its complementary damage. Intra-deformation interface failure was simulated with a composite interlayer cohesive zone element. The experimentally derived DIC-based strain fields were well-consistent with the numerical outcomes. Deeper investigation was made into the superiority of the 3D modeling, which is ascribed to the predictability for distribution of complex field variables like the free-edge effect and progressive failure accumulation within the critical sample section. The damage mechanism of the satin weave composite was explored in depth and it provides useful guidance for the practical application of the composite.

1. Introduction

Satin weave SiC_f/SiC composites play an important role in the aerospace, naval, and automotive industries. They have many advantages over metals, such as higher specific mechanical properties and a longer life cycle. In contrast, their damage tolerance and energy absorption properties are observed to be better than conventional laminated composite structures. These materials are subjected to different loading conditions during service, including shock, vibration, and fatigue.

Recently, woven composites have attracted much attention because of their excellent mechanical properties in both in-plane and out-of-plane directions. A considerable of work has been done to predict the properties [1,2,3] and to investigate the damage evolution and failure mechanisms [4,5,6].

Woven composites have different weaving forms, which can be characterized by the number of interweavings of radial and weft yarns. The number of times a radial yarn is interwoven with a weft yarn is represented by R, and the number of times a weft yarn is interwoven with a radial yarn is represented by W. In general, they are the same, i.e., R = W = N. As shown in the Figure 1, when N = 2, it is a plain structure, when N = 3, it is a twill structure, when N = 4, it is a four-satin structure, and when N = 5, it is a five-satin structure.

The mechanical properties of woven composites are closely related to N. The larger the N, the fewer the number of yarns interwoven; the lower the curl of yarn, the more it can maintain the mechanical properties of yarn itself; and the smaller the N, the larger the number of yarns interwoven, which helps to improve the stability of the structure, while the frequent interweaving will cause the curl of yarn in space, which will cause some deterioration to the mechanical properties of the yarns [7,8,9,10]. The study of yarn structure can provide a better understanding of the woven composites’ internal mesostructure.

To better simulate the damage response of satin weave composite, it is particularly important to use fidelity damage constitutive theories and numerical models in the modeling process. The damage defects of satin weave composite in a three-bending test exhibit the hybrid failure mode by the coupling action of a multi-failure mechanism in the case of classic numerical analysis. The elements’ failure modes in the finite element (FE) tool consist primarily of complicated processes, such as generation, accumulation, and evolution of damage and deletion of elements. Past decades witnessed the formulation of miscellaneous failure theories for the damage-generation forecasting of varying composite types, which were used for complex failure mechanism identification, including the criteria by Tsai-Wu [11], Hashin [12,13], Puch [14], and Chang-Chang [15]. Application of the foregoing progressive failure theories has been seen for complex failure behavior forecasting in practical application situations. Experimental studies [16,17] explored various woven laminate composites that were subjected to quasi-static tensile loads, thereby achieving their failure mechanism elucidation. Specimen failures were assessed, such as fiber cracks and chief damage modes. The failure mechanism of the composite is complex and there is wide room for research.

Compared to laboratory tests, failure modes of major compound components subjected to complex operational loads are usually more complex, such as matrix cracking and interactions between delaminations. This is because of the usual presence of 3D stresses at the critical failure sites of large engineering components, which result from local stress enrichment such as edge or bolting effects at the cut-outs and open holes. One of the important factors inducing interactions between the intra- and interlaminar failure modes is the free-edge effect. Many studies have observed the free-edge effect in composite laminates [18,19,20,21], which is due to the different material properties between adjacent laminates in different directions, resulting in an out-of-plane load imbalance near the edge of the laminate. Highly concentrated interlaminar stress fields can occur near the free edges, which can trigger interlaminar damage even in loaded composite plates. It shows that the free-edge effect is one of the main reasons for the failure of laminates.

A CCD camera was used to track the onset of damage and its progression in woven SiC_f/PyC/SiC composites plates, and DIC analysis [22] was used to generate dynamic full-field strain contour lines. To determine the failure mode, a detailed microscopic examination of the fracture site was performed. A 3D finite element (FE) model was constructed to represent the behavior of the composite specimens under loading, using the modified 3D Hashin damage model with a suitable damage progress law implemented using a VUMAT subroutine. The delamination between layers was modeled using cohesive zone elements.

2. Materials and Experiments

2.1. Processing of the Five-Harness Satin Weave SiC_f/PyC/SiC Composite Plates

The chemical vapor infiltration (CVI) and precursor infiltration pyrolysis (PIP) procedures were adopted for the fabrication of SiC_f/SiC composites with a pyrolytic carbon (PyC) interface, as shown in Figure 2.

Initially, a 2D paving technology was employed to prepare a 2D preform, whose fiber volume fractions were around 30%. The reinforcements used were the continuous SiC fibers (provided by Xiamen University, China), which were woven in a five-satin weave.

Second, CVI was employed to form a PyC interface coating on the preform fiber surface by utilizing the propylene (C₃H₆) (provided by MACKLIN), where deposition was accomplished at 850 °C.

Finally, the PIP procedure was adopted for fabricating the SiC matrix, where polycarbosilane was used as a precursor and xylene (analytically pure) was used as a diluent. The 50 wt% polycarbosilane (PCS) provided by the National University of Defense Technology Macklin with a xylene solution provided by MACKLIN was used as an impregnating solution for vacuum impregnation.

Afterward, 1 h pyrolysis of the preforms proceeded at 1000 °C. The impregnation–pyrolysis process needs to be repeated until an individual densification cycle had a less than 1% gain of weight.

In the experiment, we processed the sample with a high-speed cutting machine to cut the sample and sanded the edge of the sample to ensure that there was no burr at the edge. We cut the three-bending test specimens from a large SiC_f/PyC/SiC plate composite with aligned warp and weft yarns. The dimensional size was 39.9 mm × 6 mm × 2.9 mm (length × width × thickness). Figure 3 depicts the surface morphology of the SiC_f/PyC/SiC composites.

2.2. Digital Image Correlation

DIC (digital image correlation), a computerized stereo-vision method with the same principles as those possessed by the human eye, was put forward at the end of the 20th century [23,24,25,26] based on DIC [27] and stereo viewing. With this technique, the point correspondence between a sample’s two images, which are acquired from two separate calibrated cameras, is defined by exploiting the DIC algorithm.

For the specimen surface, its 3D-coordinates reconstruction is possible from the point correspondence using triangulation. The principle of DIC is depicted in Figure 4. Derivation of the 3D displacement field is accomplished based on correlation between pre- and post-deformation images (corresponding separately to Figure 4’s initial and displaced geometries). Numerical differentiation is employed to transform the displacement field into a strain field.

The spacing between the two points’ locational coordinates is equivalent to the post-deformation displacement. The correlation factor is equivalent to the equation shown below [28]:

$$C\left(u,v\right)=\frac{{\Sigma }_{i=1}^m{\Sigma }_{j=1}^m\left[f\left(x_i,y_i\right)-\overline{f}\right]\left[g\left(x_i^{*},y_i^{*}\right)-\overline{g}\right]}{\sqrt{{\Sigma }_{i=1}^m{\Sigma }_{j=1}^m{\left[f\left(x_i,y_i\right)-\overline{f}\right]}^2}\sqrt{{\Sigma }_{i=1}^m{\Sigma }_{j=1}^m{\left[g\left(x_i^{*},y_i^{*}\right)-\overline{g}\right]}^2}}$$

(1)

where the correlation factor $C\left(u,v\right)$ is dependent on vertical displacement $v$ and horizontal displacement $u$, while the coordinates $\left(x_i, y_i\right)$and $\left(x_i^{*}, y_i^{*}\right)$ show linkage to the translations between two digital images. The functions $f\left(x_i, y_i\right)$ and $g\left(x_i, y_{\mathrm{i}}\right)$are the grayscale values at point $\left(x,y\right)$ separately in the original and translated images, while $\overline{f}$ and $\overline{g}$ stand for the corresponding means of the intensity matrices.

In the case of a small displacement vertical to the optical axis of camera, Equations (2) and (3) can be used to approximate the correlation of $\left(x_i, y_i\right)$ with $\left(x_i^{*}, y_i^{*}\right)$ through a 2D transformation, where $u_0$ and $v_0$ separately denote the horizontal and vertical translations:

$$x^{*}=x+u_0+\frac{du}{dx}dx+\frac{du}{dy}dy$$

(2)

$$y^{*}=y+v_0+\frac{dv}{dx}dx+\frac{dv}{dy}dy$$

(3)

3. Experimental Setup and Procedure

As depicted in Figure 5a, the three-point bending test of the specimens was performed by the INSTRON 5966 universal electronic tester (INSTRON, Boston, MA, USA) with a ±10 KN force range. Competently, the tester had force-control accuracy within 0.5% and displacement-control accuracy within 0.1%. During the bending tests, the deformation loading rates were set at 0.4 mm per min. Figure 5b shows a schematic diagram of the experimental process.

For the intra-experimental investigation into the strain layout and displacement properties of the SiC_f/PyC/SiC composites, the 2D DIC strategy was employed to assess the strain fields and deformations on the sample surfaces. Spots were sprayed on the sample surfaces prior to the experiment, which was followed by the imaging system arrangement in front of the sample, where the optical axis of the camera was vertical to the spot surface (Figure 5).

A telephoto lens and a CCD camera (12.3M model produced by FLIR company, Portland, OR, USA) were utilized for documenting the images at a frequency of 1 Hz. DIC computational assessment of pre- and post-deformation speckle images was accomplished via Vic-2D 2009 software (Correlated Solutions, Irmo, SC, USA), where the subset size was 29 pixels and the step size was 5 pixels.

4. Numerical Modeling

4.1. Five-Harness Five-Layer Satin Weave SiC_f/SiC Plate Composites Model Hypothesis

To reduce the number of calculations in the simulation, assumptions are made.

They were made when satin weave SiC_f/SiC plate composites model was established: in the process of composite material preparation, the interfacial properties between fibers and matrix are good, and there are no cracks, bubbles, and other defects in the specimens and each layer of the composite material is regarded as anisotropic homogeneous material.

4.2. Numerical Model for Bending Test

A three-dimensional FE model for a three-point bending test of five-harness five-layer satin weave SiC_f/PyC/SiC plate composites was created (Figure 6). The ABAQUS FE-based qusai-static simulations were run under an explicit scheme for time integration, where the simulated workpiece material had a size of 39.9 × 6 × 2.9 mm. For the laminate meshing, eight-node solid elements (C3D8R) with reduced integration and hourglass control were adopted. Since the mechanical behavior of the interface is little influenced by its initial thickness, in this analysis the damage behavior of the composite interface was simulated with the zero-thickness COH3D8 cohesive elements. The optimal size of elements was 0.2 mm in the mesh research after considering a compromise between computational time and outcome accuracy.

Eight-node solid elements with full/reduced integration and hourglass control are utilized.

In Figure 6, the boundary conditions of the finite element model are given. The finite element model is fixed by the way of contact between the indenter with a constant rate and a fixed base, which was meshed by coarse mesh and constrained as a rigid body to improve computational efficiency. In addition, a smaller friction coefficient of 0.05 is adopted to speed up the convergence of the calculation. The Hashin failure criterion [29] and the progressive damage model were used to simulate the damage-mechanical behavior of the composite materials and the cohesive zone modeling was used to simulate the delamination of laminates.

4.3. Damage Initiation Criteria

For the five-harness satin weave SiC_f/PyC/SiC composites, their fracturing under a bending load exhibits pseudoplasticity in contrast to the delamination of the interface and the breakage of yarns. Furthermore, a few failure mechanisms are concurrently integrated regarding the satin weave composites. There are often two components for every failure mechanism: a damage generation criterion and a damage progression law.

Due to the structural properties of five-harness satin weave composites, the in-plane properties (X and Y directions, as shown in Figure 1D and Figure 3) are not significantly different; the properties in the thickness direction (Z-direction) are more different from the in-plane.

The damage generation criterion applied varies by the constituent type: Hashin’s failure criterion is applied for satin weave composites and the traction–separation damage law for the interface of the composite. The corresponding associations of failure modes with their damage generations are formulated as follows:

Tensile failure in X-direction:

$${\varphi }_{Xt}={\left(\frac{{\sigma }_X}{F_X^t}\right)}^2+\alpha {\left(\frac{{\sigma }_{XY}}{F_{XY}^S}\right)}^2+\alpha {\left(\frac{{\sigma }_{XZ}}{F_{XZ}^S}\right)}^2⩾1$$

(4)

Compressive failure in X-direction:

$${\varphi }_{Xc}={\left(\frac{{\sigma }_X}{F_X^c}\right)}^2⩾1$$

(5)

Tensile failure in Y-direction:

$${\varphi }_{Yt}={\left(\frac{{\sigma }_Y}{F_Y^t}\right)}^2+\alpha {\left(\frac{{\sigma }_{XY}}{F_{XY}^S}\right)}^2+\alpha {\left(\frac{{\sigma }_{YZ}}{F_{YZ}^S}\right)}^2⩾1$$

(6)

Yarn compressive failure in Y-direction:

$${\varphi }_{Yc}={\left(\frac{{\sigma }_Y}{F_Y^c}\right)}^2⩾1$$

(7)

Z-directional tensile and shear failures:

$${\varphi }_{Zt}={\left(\frac{{\sigma }_Z}{F_Z^t}\right)}^2+{\left(\frac{{\sigma }_{XZ}}{F_{XZ}^S}\right)}^2+{\left(\frac{{\sigma }_{YZ}}{F_{YZ}^S}\right)}^2⩾1$$

(8)

Z-directional compressive and shear failures of yarn:

$${\varphi }_{Zc}=\left[{\left(\frac{F_Z^c}{2F^s}\right)}^2-1\right]\frac{{\sigma }_Z}{F_Z^c}+{\left(\frac{{\sigma }_Z}{2F^s}\right)}^2+{\left(\frac{{\sigma }_{XZ}}{F_{XZ}^s}\right)}^2+{\left(\frac{{\sigma }_{YZ}}{F_{YZ}^s}\right)}^2⩾1$$

(9)

In the above equation, $F_X^t, F_Y^t$, and $F_Z^t$ stand for the X-, Y-, and Z-directional tensile strengths, respectively. $F_X^c, F_Y^c$, and $F_Z^c$ stand for the X-, Y-, and Z-directional compressive strengths of yarn, respectively. $F_{XY}^s, F_{YZ}^s$, and $F_{XZ}^s$ signify XY, YZ, and XZ shear strengths, respectively. Due to the in-plane nature of the composite material, it is assumed that $F^s=F_{XZ}^s=F_{YZ}^s$. $\alpha$ is the contributing factor in each mode.

4.4. Damage Evolvement Model

Along with the softening behavior of the material model, the material damage is characterized by localization. Given the declining dissipation of energy with the mesh refinement, the size of the element is an extreme determinant of the numerical solution. Hence, objectivity is lost in the numerical outcome apart from the damage progression equation that is associated with element size, local stress/strain, and other material traits. For the spatial problem, each fracture mode’s fracture-energy density is assumed to correspond with a fixed value, while the failure strain value will fluctuate with the size of the element.

$$\frac{1}{2}{\epsilon }_{I,f}{\sigma }_{I.f}{l}^3=G_I{l}^2$$

(10)

where $l$ refers to the eigenlength of the element, and ${\sigma }_{I,f}, {\epsilon }_{I,f},$and $G_I$ denote the equivalent peak stress, the failure equivalent strain, and the fracture energy density for a failure mode I, respectively.

The equivalent deformation at the failure point is formulated as:

$$X_{\mathrm{eq}}^f={\epsilon }_{I,f }l$$

(11)

Table 1. The equivalence deformations and stresses under various failure modes.

Failure Modes	Equivalence Displacement	Equivalence Stress
$X,{\sigma }_{11}⩾0$	$X_{\mathrm{eq}}^{\mathrm{Xt}}=l\sqrt{\left(\langle {\epsilon }_{11}\rangle {}^2+\alpha {\epsilon }_{12}^2+\alpha {\epsilon }_{31}^2\right)}$	$l\left(\langle {\sigma }_{11}\rangle \langle {\epsilon }_{11}\rangle +\alpha {\sigma }_{12}{\epsilon }_{12}+\alpha {\sigma }_{13}{\epsilon }_{13}\right)/X_{\mathrm{eq}}^{\mathrm{Xt}}$
$X,{\sigma }_{11}<0$	$X_{\mathrm{eq}}^{\mathrm{Xc}}=l\langle -{\epsilon }_{11}\rangle$	$l\langle -{\sigma }_{11}\rangle \langle -{\epsilon }_{11}\rangle /X_{\mathrm{eq}}^{Xc}$
$Y,{\sigma }_{22}⩾0$	$X_{\mathrm{eq}}^{\mathrm{Yt}}=l\sqrt{\left(\langle {\epsilon }_{22}\rangle {}^2+\alpha {\epsilon }_{12}^2+\alpha {\epsilon }_{23}^2\right)}$	$l\left(\langle {\sigma }_{22}\rangle \langle {\epsilon }_{22}\rangle +\alpha {\sigma }_{12}{\epsilon }_{12}+\alpha {\sigma }_{23}{\epsilon }_{23}\right)/X_{\mathrm{eq}}^{\mathrm{Yt}}$
$Y,{\sigma }_{22}<0$	$X_{\mathrm{eq}}^{\mathrm{Yc}}=l\langle -{\epsilon }_{22}\rangle$	$l\langle -{\sigma }_{22}\rangle \langle -{\epsilon }_{22}\rangle /X_{\mathrm{eq}}^{Yc}$
$Z,{\sigma }_{33}⩾0$	$X_{\mathrm{eq}}^{\mathrm{Zt}}=l\sqrt{\left(\langle {\epsilon }_{33}\rangle {}^2+\alpha {\epsilon }_{23}^2+\alpha {\epsilon }_{13}^2\right)}$	$l\left(\langle {\sigma }_{33}\rangle \langle {\epsilon }_{33}\rangle +\alpha {\sigma }_{13}{\epsilon }_{13}+\alpha {\sigma }_{23}{\epsilon }_{23}\right)/X_{\mathrm{eq}}^{\mathrm{Zt}}$
$Z,{\sigma }_{33}<0$	$X_{\mathrm{eq}}^{\mathrm{Zc}}=l\langle -{\epsilon }_{33}\rangle$	$l\langle -{\sigma }_{33}\rangle \langle -{\epsilon }_{33}\rangle /X_{\mathrm{eq}}^{Zc}$

summarizes the equivalent deformations of the initiation damage. Besides, the damage progression equation for every failure mode is defined by the equivalent deformations shown below:

$$d_I=1-\frac{X_{\mathrm{eq}}^{Ii}\left(X_{\mathrm{eq}}^{If}-X_{\mathrm{eq}}^I\right)}{X_{\mathrm{eq}}^I\left(X_{\mathrm{eq}}^{If}-X_{\mathrm{eq}}^{Ii}\right)},\left(I=Xt,Xc,Yt,Yc,Zt,Zc\right)$$

(12)

In the above equation, $X_{eq}^{li}$ and $X_{eq}^{If}$ denote the initiation and full-damage equivalent deformations, respectively, under failure mode I, whose computational formulas are shown below:

$$X_{\mathrm{eq}}^{Ii}=X_{\mathrm{eq}}^I/ \sqrt{{\varphi }_I}$$

(13)

$$X_{\mathrm{eq}}^{If}=2G_I/{\sigma }_{\mathrm{eq}}^{Ii}$$

(14)

where ${\varphi }_I$ stands for the damage generation criterion value, while $G_I$ and ${\sigma }_{eq}^{li}$ stand for the fracture energy density and initiation equivalent stress, respectively, under failure mode I. The computational formula for ${\sigma }_{eq}^{li}$ is displayed below:

$${\sigma }_{\mathrm{eq}}^{Ii}={\sigma }_{\mathrm{eq}}^I/\sqrt{{\varphi }_I}$$

(15)

where ${\sigma }_{eq}^I$ refers to the equivalence stress as formulated in Table 1. The damage progression equation is correlated with the element eigenlength, the local strain, and the fracture energy of the braided composite constituents.

The damage of composites can be elucidated by exploiting the Murakami–Ohno scheme, where the damage state is expressed with three principal parameters, as follows:

$$D=\sum D_i{n}_i\otimes n_i,i=\left(L,T,Z\right)$$

(16)

where $D_i$ and $n_i$ denote the damage tensor’s principal value and principle unit vector, respectively.

To take into account the stress tensor symmetry, the effective stress can be subjected to the following processing:

$${\sigma }^{*}=\frac{1}{2}\left[{(I-D)}^{-1}\sigma +\sigma {(I-D)}^{-1}\right]=M\left(D\right)\sigma$$

(17)

where $\sigma$ and ${\sigma }^{*}$ denote the undamaged and effective stress, respectively. For the damage parameter expression into the stiffness matrix, the Cordebois–Sidoroff principle of energy equivalence is implemented as follows:

$$C\left(D\right)=M^{-1}\left(D\right):C:M^{T,-1}\left(D\right)$$

(18)

where $C$ and $C\left(D\right)$ stand for the undamaged and damaged configurations’ stiffness tensors, respectively.

The matrix $C\left(D\right)$ for damaged stiffness can be formulated with the undamaged stiffness matrix components in conjunction with the damage tensor $D_I$’s principal values. Its more explicit expression is shown below:

$$\begin{array}{c}C\left(D\right)=\left[\begin{array}{cccccc}{b}_X^2{Q}_{11}& b_X{b}_Y{Q}_{12}& b_X{b}_Z{Q}_{13}& 0& 0& 0\\ & b_Y^2{Q}_{22}& b_Y{b}_Z{Q}_{23}& 0& 0& 0\\ & & b_Z^2{Q}_{33}& 0& 0& 0\\ & sym& & b_{XY}{Q}_{44}& 0& 0\\ & & & & b_{XZ}{Q}_{55}& 0\\ & & & & & b_{YZ}{Q}_{66}\end{array}\right]\end{array}$$

(19)

where$b_X=1-D_Y$,$b_Y=1-D_Y$,$b_Z=1-D_Z$, $b_{XY}={\left(\frac{2\left(1-D_X\right)\left(1-D_Y\right)}{2-D_X-D_Y}\right)}^2$,$b_{XZ}={\left(\frac{2\left(1-D_X\right)\left(1-D_Z\right)}{2-D_X-D_Z}\right)}^2$, and$b_{YZ}={\left(\frac{2\left(1-D_Y\right)\left(1-D_Z\right)}{2-D_Y-D_Z}\right)}^2$, and $D_X=max\left(d_{Xt},d_{Xc}\right),{ D}_Y=max\left(d_{Yt},d_{Yc}\right),{ \mathrm{and} D}_Z=max\left(d_{Zt},d_{Zc}\right)$. $Q_{ij}$refers to the undamaged stiffness matrix component. The properties given in

Table 2. Mechanical properties of SiC_f/SiC composite [30].

${\mathit{E}}_{\mathit{X}}\left(GPa\right)$	${\mathit{E}}_{\mathit{Y}}\left(GPa\right)$	${\mathit{E}}_{\mathit{Z}}\left(GPa\right)$	${\mathit{v}}_{\mathit{X}\mathit{Y}}$	${\mathit{v}}_{\mathit{X}\mathit{Z}}={\mathit{v}}_{\mathit{Y}\mathit{Z}}$	${\mathit{G}}_{\mathit{X}\mathit{Y}}\left(GPa\right)$	${\mathit{G}}_{\mathit{X}\mathit{Z}}={\mathit{G}}_{\mathit{Y}\mathit{Z}}\left(GPa\right)$
51.031	50.751	10.53	0.31	0.42	4.36	3.048

and

Table 3. Strength and damage data for SiC_f/SiC composite.

Strength [31,32]	${\mathit{X}}_{\mathit{T}}\left(\mathbf{MPa}\right)$	${\mathit{Y}}_{\mathit{T}}\left(\mathbf{MPa}\right)$	${\mathit{S}}_{\mathit{X}\mathit{Y}}\left(\mathbf{MPa}\right)$	${\mathit{S}}_{\mathit{Y}\mathit{Z}}\left(\mathbf{MPa}\right)$
	$X_C\left(\mathrm{MPa}\right)$	$Y_C\left(\mathrm{MPa}\right)$	$S_{XZ}\left(\mathrm{MPa}\right)$
	754	734	56	56
	618	621	56
Damage [33]	$G_{Xt}=G_{Yt}\left(\mathrm{N}/\mathrm{mm}\right)$	$G_{Xc}=G_{Yc}\left(\mathrm{N}/\mathrm{mm}\right)$	$G_{Zt}\left(\mathrm{N}/\mathrm{mm}\right)$	$G_{Zc}\left(\mathrm{N}/\mathrm{mm}\right)$
	8	1.5	0.632	0.306

were used to achieve behavioral simulation modeling of the studied five-harness satin weave composite. In Table 3, $X_T$ and $X_c$ represent the tensile and compressive strength of the composite in the X direction, respectively, and so on; $G_{Xt}$ and $G_{Xc}$ represent the tensile and compression fracture energy density of the composite in the X direction, respectively, and so on.

4.5. Traction-Separation Damage Law

Modeling of inter-ply delamination is accomplished with the aid of ABAQUS/Explicit using the CEs (cohesive elements). Camaho et al.’s [34] bilinear CEs model is employed in the present work.

The secondary stress criterion for damage initiation of the cohesive element is as follows:

$${\left(\frac{\langle {\sigma }_n\rangle }{N_{max}}\right)}^2+{\left(\frac{{\sigma }_s}{S_{max }}\right)}^2+{\left(\frac{{\sigma }_t}{T_{max }}\right)}^2=1$$

(20)

In the following equations, ${\delta }_i$ stands for the hybrid-mode deformation,${\delta }_m$ denotes the equivalent deformation, ${\delta }_m^{max}$ stands for the maximum hybrid-mode deformation, ${\delta }_m^0$ represents the effective deformation at the damage generation, ${\delta }_m^f$ represents the hybrid-mode deformation at complete damage (BK Criteria), and$D_{mix}^s$ refers to the hybrid-mode damage parameters for the CEs. Their respective expressions are displayed below:

$$\begin{array}{c}{\delta }_m=\sqrt{{\delta }_n^2+{\delta }_s^2+{\delta }_t^2} , {\delta }_m^{max}=\max \left({\delta }_m^{max},{\delta }_m\right)\end{array}$$

(21)

$${\delta }_m^0=\left\{\begin{array}{cc}{\delta }_t^0{\delta }_n^0\sqrt{\frac{1+{\beta }^2}{{\left({\delta }_t^0\right)}^2+{\left(\beta {\delta }_t^0\right)}^2}}\hfill & {\delta }_t>0\hfill \\ \sqrt{{\left({\delta }_n^0\right)}^2+{\left({\delta }_s^2\right)}^2}\hfill & {\delta }_t⩽0\hfill \end{array}\right.$$

(22)

$${\delta }_m^f=\frac{2}{K_i{\delta }_m^0}\left[G_n+\left(G_s-G_{\mathrm{n}}\right)\left.{\left(\frac{G_s+G_t}{G_n+G_s+G_t}\right)}^{\eta }\right]\right.$$

(23)

$$D_{mix}^s=\frac{{\delta }_m^f\left({\delta }_m^{max}-{\delta }_m^0\right)}{{\delta }_m^{max}\left({\delta }_m^f-{\delta }_m^0\right)}$$

(24)

In the equations, $\langle · \rangle$ denotes the MacAuley bracket, and $G_n$,$G_s$, and $G_t$ stand for the inter-laminar fracture energies under modes N, S, and T, respectively, which represent cracks upon the separate impositions of opposite normal/tangential force-deformations on the upper and lower element surfaces, as illustrated in Figure 7. Furthermore, ${\delta }_n , {\delta }_s , \mathrm{and} {\delta }_t$ refer to the deformations in the normal, first shear, and second shear directions, respectively, and η, with a value of 1.45, represents the variable defined in the Benzeggagh–Kenane fracture criterion [35]. The mode mixing ratio β, on the other hand, can be formulated as:

$$\beta =\sqrt{\frac{\left({\delta }_n^2+{\delta }_s^2\right)}{{\delta }_t^2}}, {\delta }_t>0$$

(25)

where $T_m$ denotes the nominal stress under the hybrid mode:

$$\left\{\begin{array}{c}{T}_m=K_i{\delta }_m, {\delta }_m^{max}⩽{\delta }_m^0 \hfill \\ T_m=\left(1-D_{mix}^s\right)K_i, {\delta }_m^0⩽{\delta }_m^{max}<{\delta }_m^f \hfill \\ T_m=0, {\delta }_m^{max}⩾{\delta }_m^f\hfill \end{array}\right.$$

(26)

where $K_i$(i = n,s,t standing for the model N, S, and T delaminations, respectively) represents the constitutive stiffness for CEs.

Table 4. Interface properties of five-harness satin weave SiC_f/PyC/SiC composite [37].

$N_{max}\left(\mathrm{MPa}\right)$	45	$S_{max}=T_{max}\left(\mathrm{MPa}\right)$	35
$G_n\left(\mathrm{N}/\mathrm{mm}\right)$	0.3	$G_s=G_t\left(\mathrm{N}/\mathrm{mm}\right)$	1.0

details the interface layer-related constants. The stiffness for resin-rich interpile interface was expressed as $K=E_{33}/t$ by Daudeville et al. [36], where $E_{33}$ (=10.5 GPa) denoted the thickness-directional stiffness of the material and t (=35 μm) represented the thickness. Accordingly, the interface stiffness during the simulations was set at about 3 × 10⁵ N/mm³.

The creation of jobs was accomplished based on the foregoing settings, and the errors were examined by outputting the corresponding input file. For the computational efficiency enhancement while preserving accuracy, the mass scaling factor was assigned as 10^3 with reference to [38].

4.6. Validation of Mesh Convergence

Mesh quality is a crucial factor for the forecasting accuracy. The mesh selection was accomplished through refinement analysis, where the models’ normalized Bending moduli of elasticity (divided by the effective model moduli with a 0.8 mm edge length of elements). The edge lengths of elements were set at 0.8, 0.3, and 0.15 mm herein. Limited variations (below 2%) are noted between the normalized moduli of varying meshes, as depicted in Figure 8. It is also distinctly visible from the local diagrams that the meshed models show variations with the size of element. Insignificant changes in the trends and outcomes from the foregoing analyses are expected with further refinements of the meshes, and the outcomes are well-consistent with each other. However, the computational duration prolongs significantly with the diminishing mesh size. The element size is finalized as 0.3 mm after trading off between the simulation accuracy and the computing expenditure.

5. Results and Discussions

In this section, the damage mechanisms in satin weave SiC_f/PyC/SiC composite plates is discussed. Experimental and numerically obtained load-displacement curves, damage variable fields of the model, electron microscope views of material fractures, and the strain fields on the sample surface were used to analyze the material’s behavior.

According to the experimental standard [39], a total of five samples were tested in this experiment, as shown in Figure 9, and the experimental results have good consistency. The force-displacement curve of sample #3 (bule one) is good, which is compared with the simulated curve, and its bending modulus is 48 GPa and bending strength is 246 MPa, respectively.

Figure 10 displays the force-deformation graphs derived experimentally and computationally during the three-point bending test of SiC_f/PyC/SiC composite. Given the explanation of nonlinearity effects by the constitutive law, the numerical model-based forecasts consistently satisfy the experimental findings. The outstanding mechanical performance of the woven composite endows the specimen with competent load endurance at larger deformations. There are three phases regarding its response. Phase 1, the linear elastic stage, is noted in both the experimental and simulated curves, with displacement values approximately ranging from 0 mm to 0.1 mm. A relevant feature can be an elastic modulus $E_{bending}$. In phase 2 (about 0.1 mm < displacement < 0.11 mm in the experiment and 0.1 mm < displacement < 0.15 mm in the simulation), the force-deformation graph is no longer linear. Relative warp and weft reorientation describes the nonlinear property of the woven composite in this phase and the PyC interface layer plays a certain role in hindering the formation of main cracks, where the microcracking commences with the appearance in matrix. Eventually, in phase 3 (about 0.11 mm < displacement in the experiment and 0.15 mm < displacement in the simulation), fiber fracturing happened in the studied composite. By the load action, such damage progressed persistently to result in delamination. To sum up, the behavior of the force-deformation graphs is thus ‘ductile nonlinear’ [31].

The morphology of SiC_f/PyC/SiC composites on the microscopic scale is shown in Figure 11. During the pyrolysis process, due to the volume contraction of the precursor of the composite, a series of micropores are formed during the densification process, as shown in Figure 11a. As shown in Figure 11b, there is a distinct PyC interface between the fibers and the matrix with a thickness of about 4 μm, and no significant gap is found between the fibers and the matrix. During the load-bearing process of the composite material, microcracks initiated at the existing micropores in the composite material, then formed large cracks. Due to the laminar structure of the PyC interface, it has a certain deflection effect on the cracks, which expand in the interface, dissipate certain strain energy, and delay the fiber fracture, and the final result graph is shown in Figure 11c.

The full-field strain distributions of the specimen surface were extracted with the recorded images during the experiments. Figure 12 shows the full-field ${\epsilon }_{xx}$, ${\epsilon }_{yy}$, and ${\epsilon }_{xy}$ strain contours under the force of 180 N, 210 N, 165 N, 180 N, 145 N, 130 N, and 30 N, respectively, where each corresponds to the letters a ~ g as shown in Figure 10. It can be seen from Figure 12 that the strain distributions on the surface of the satin woven composite exhibit regular variation. Since the value of ${\epsilon }_{xx}$ is greater than the values of ${\epsilon }_{yy}$ and ${\epsilon }_{xy}$ on strain distribution obtained by DIC, the primary damage mode of the specimen is X-directional tensile damage. When the composite is in the damage extension mode, the X-directional tensile strain at the bottom of the specimen expands rapidly until the material complete failure.

Since the matrix strain is larger and strength lower, generation and propagation of microcracks are noted on the sample surface along the matrix (Figure 11). The matrix cracks are suggestive of yarn microfailure in the composite as well. Due to the declined load capacity of failed yarns and shifting of more load to the adjoining matrix, enhanced stress and strain of the matrix result, as well as the commencement of cracking. At an applied load of approximately 180 N, the micro-cracks are generated in various matrix sites, which subsequently propagate to form a primary crack, followed by the sample fracturing. With the intensification of the applied load, a primary crack is formed by small propagated cracks. The fracturing and pull-out of fibers are noted when the stress on the cross-section of cracks approximates the tensile strength. At stresses below 170 N, there are no apparent cracks on the surface. Since the existence of the PyC interface allows the fiber damage to proceed in several stages, not as a sudden fracture, giving the composite pseudoplastic characteristic.

Figure 13 shows the specimen surface distribution of ${\epsilon }_{xx}$ upon the specimen damage generation (SDV20 represents ${\epsilon }_{xx}$), which is derived based on the FE simulations. For the specimen, its lower side is under tensile strain and the upper side is under compressive strain and contributes to the three-point bending applied to the specimen, the lower and upper sides of which are separately under tensile and compressive stresses primarily. Benefited from the mappings, the deformation of large areas of the material surface can be determined, which is consistent with the strain distribution measured by DIC method in Figure 12. The tensile and compressive strains on the surface of the specimen are focused in a limited range because of the PyC interface layer, indicating that a PyC interface with suitable stiffness can concentrate the material’s deformation in a short range. It is possible to get a satisfactory correlation between experimentally determined and anticipated failure location.

Figure 14 shows the damage distribution for the fiber tension and matrix compression on the sample (SDV1 represents tensile damage and SDV2 represents compression damage). Compared with Figure 12, it is found that Figure 12 and Figure 14 have a good consistency, the main damage mode of the composite on the lower side is tensile damage, and the upper side is mainly compressive damage attributed to the tensile strain mainly applied on the lower side of the specimen and the compressive strain mainly applied on the upper side, which can also be well-verified in the strain contours of each stage in Figure 12. As seen in the local graphs, both tensile and compression damage are significantly reduced when passing through different layups. This is because the PyC interface layer plays a certain hindering effect on the damage extension. The main mechanism is that the proper interfacial strength deflects the tiny cracks that will cause damage to the fibers, thus pushing back the process of forming main cracks and delaying the failure process of the material to some extent, and finally make the specimen exhibit pseudoplasticity.

As shown in Figure 15 (SDEG represents interface damage), the interface damage initiates at the free edge of the specimen and attributed to the free edge effect, then the interfacial damage continues to expand until the interface delaminates. In the three-dimensional stress state, the junction of the composite laminate will show a certain free-edge effect, which is mainly manifested as the stress concentration phenomenon, making the delamination damage occur when the material is far below its service strength, leading to the early failure of the material, so the interface damage value at the free edge is the biggest. However, due to the existence of the PyC interface and the appropriate interface stiffness and strength of the interface, the tiny cracks will deflected, delaying the time of the formation of main cracks to the fibers, and increasing the overall damage tolerance of the specimens, so the satin weave SiC_f/PyC/SiC composite plates shows a certain pseudoplastic characteristics. It can be noted that the delamination is an effective toughening mechanism of the laminated composite material whereby the fracture energy dissipates at the expenses of a translaminar cracking conducting to the failure.

The difference between experiment and simulation might be related to the following:

• The FE model neglects some details of the specimen, which would greatly increase the computational efficiency and convergence difficulty but have a slight influence on the simulation results.
• Some technical defects, such as geometric tolerance and voids in the matrix, are inevitable, so that might lead to the material properties having a deviation from actual material properties.
• There might be some tiny gaps between the testing machine, testing tools, and the joint at the beginning of the tests.

6. Conclusions

In the present paper, the failure behavior for five-harness satin weave SiC_f/PyC/SiC composite was investigated experimentally and numerically through three-point bending experimentation. The failure generation and progression in the studied composite were elucidated by employing the DIC procedure. The following conclusions were reached:

• The five-harness five-layer satin weave SiC_f/PyC/SiC composite was subjected to a three-point bending experimentation for subsequent usage in the created FE model, where the material behavior modeling was accomplished following the altered 3D Hashin failure criterion that was implemented using VUMAT subroutine. The experimental findings were well-consistent with the simulations regarding both the strain fields and the load-deformation curves.
• The existence of a PyC interface increases the damage accumulation time before the complete failure of the composite, increases the warning time of material failure, improves the safety of the application of the composite, and enables the composite to be applied in wider applications.
• In a complex state of stress, even if the mechanical properties of different composite layers are similar, there will be a certain free-edge effect when there is an interface layer, and a stress concentration phenomenon will occur. The phenomenon of free edge will cause delamination damage of the composite, degrade the bearing potential of the composite, and make the composite fail in advance.
• It was found that fiber fracture is the dominating failure mechanism for the studied specimens.
• Each layer of the composite material is regarded as anisotropic homogeneous material, so the influence of the weaving mode of warp and weft on the material properties cannot be obtained. It can be further studied by establishing a mesoscopic model.