An Advanced Finite Element Modeling for the Failure of Notched Ceramic Matrix Composite With TFP Patch Reinforcement

An advanced modeling strategy for notched ceramic matrix composite coupons with patch reinforcement was proposed to investigate the failure mechanisms. This model considered the tailored fiber–placed (TFP) yarn details obtained from the design phase and the embedded element concept which was used to successfully overcome the meshing difficulties. Inter-ply “glue” layers were simulated using the surface-based contact cohesive element method, so the delamination due to interfacial material discontinuity and damage can be well reproduced and analyzed. For composite ply, the energy-based composite progressive damage model that is independent of the mesh size was applied. Virtual test campaign was performed with a variety of geometrical and material parameters, and the damage and failure mechanisms based on the stress analysis can be revealed to support the design optimization of patch reinforcement.


INTRODUCTION
Continuous silicon carbide fiber-reinforced silicon carbide ceramic matrix (SiC f /SiC) composites are considered to be candidate accident-tolerant fuel (ATF) cladding materials due to their excellent ablation resistance, corrosion resistance, and radiation tolerance (Snead et al., 2007). However, despite the substantial achievements, there are still some key problems that need to be solved to satisfy the requirement of engineering application. For example, it is often inevitable to manufacture open holes or notches in the engineering composite components. In the assembling process of SiC f / SiC composite flame holders (Berdoyes et al., 2005), some hole flanges are needed for joining ( Figure  1A). In the roller kiln or tunnel kiln system, there are also some holes in the silicon carbide ceramic beam and rod for assembly ( Figures 1B,C). Nevertheless, the SiC f /SiC composites have the properties of high hardness and brittleness, and the stress concentration caused by notches, holes, or attachment points is a threat to structural safety and reliability.
The presence of a hole or notch in a composite laminate is associated with stress concentration and out-of-plane stresses and will produce a change in the failure mechanisms and a failure strength reduction of the laminate compared to a laminate without a notch. In the breakage of an open hole or notch laminate subjected to engineering loads, several failure mechanisms can be observed including matrix cracking, fiber breakage, and delamination. Consequently, local reinforcement around an open hole or notch is often required in the composite structure design. Although the mechanical joining technique with bolted reinforcement widely used in the traditional engineering materials such as steel, aluminum, or plastic provides an easy engineering solution, the drilled holes required for the bolted joints in the composite laminate will further remarkably reduce the strength of the structure by additional fiber breakage and delamination. The detailed damage introduced by the drilling process can be found in the works of Hocheng and Tsao (2006), Tsao and Hocheng (2004), and Srinivasa Rao et al. (2008). Consequently, innovative reinforcement methods for the composite laminate with an open hole or notch are required in engineering structure design.
The tailored fiber placement (TFP) technology that is based on embroidery machines was developed by the Institute of Polymer Research Dresden (Crothers et al.;Gliesche and Feltin;Rothe and Feltin, 1994). According to the designed reinforcement paths, the machine will stitch a yarn onto a base material in almost any desired orientation. Other than isotropic materials whose stress concentration can be only reduced by geometrical optimization, the anisotropic materials, such as FRPs, can be optimized based on the tailored fiber orientations that locally follow the principal stresses. Therefore, the maximum properties can be achieved along the fiber direction. Meanwhile, the stress concentration introduced by notches can be efficiently reduced to the largest extent. Gliesche et al. (2003) applied the TFP technique to a notched textile composite specimen. It was reported that the damage of the reinforced sample started outside of the reinforced area. The strength of the sample was 94% that of the un-notched specimen and 33% higher than that of the one without reinforcement.
The notched sample has been numerically studied since the 1970s. Earlier researchers employed the finite element (FE) method (Hillerborg et al., 1976;Kwon and Berner, 1995;Kwon and Craugh, 2001), in which the damage modes were considered and micromechanical constitutive law was applied. Thereafter, continuum damage modeling (CDM) (Chang and Chang, 1987;Davila et al., 2005;Forghani et al., 2013) and discrete modeling, e.g., cohesive element method (Turon et al., 2010), were developed to study the damage caused by stress concentration.
In contrast to the sample with an open hole or with traditional reinforcements, however, the notched sample reinforced by the TFP patch has rarely been studied. Koricho et al. (2015) investigated the TFP-made specimen in comparison with drilled composite samples, both experimentally and numerically. The FE model of Koricho was based on the 2D Hashin failure criterion. Temmen et al. developed a tool chain to optimize the TFP process. The Puck failure criterion was adopted in order to distinguish the fiber rupture from the inter-fiber fracture of the stitched yarn. Gliesche et al. (2003) carried out research on the reinforcement of the TFP patch on the unidirectional laminate. However, the work of Koricho, Temmen, and Gliesche did not explicitly reproduce the single yarn embedded in the patch that exhibited high inhomogeneity. Meanwhile, the ply-ply and patch-sample interfacial properties, which affect the damage initiation location and propagation, were not accounted for.
Furthermore, the studies on mechanical properties of the reinforced sample with an open hole or notch are mainly focused on polymer matrix composites rather than ceramic matrix composites. This is partly due to the fact that ceramic matrix composite technology is still being developed and the ceramic matrix composite experiments are time-, money-, and energy-consuming. However, the stress concentrations caused by a hole or notch are inevitable with the gradual application of ceramic matrix composite materials. As a result, it is vital to investigate the strength of open-hole ceramic matrix composites and the related reinforcement methods.
In this paper, an advanced three-dimensional modeling strategy for the TFP patch-reinforced ceramic matrix composite sample with an open hole is presented. It possesses the following features: 1) The placed yarn that is consolidated in the patch is fully simulated as a piece-wise orthotropic material using solid elements. It is nested in the matrix mesh using the "embedded element technique." 2) To predict the material damage of the yarn and the sample with an open hole that is strengthened by woven fabrics, an energy-based damage evolution mode proposed by Maimí et al. (2007) is applied. The damage in a ply is triggered by the maximum stress failure criterion.
3) The interface between the TFP patch and the composite sample, as well as those between plies in the composite sample, is simulated using the surface-based cohesive behavioral properties, which can well represent a quasizero-thickness cohesive layer.
A Python script is written to automatically generate the FE model according to user definitions, e.g., coupon size, hole diameter, size of the reinforcement, number of plies, fiber orientations, and material properties. With this FE model, the progressive damage of the TFP yarn, woven composite plies, and interfacial delamination can be analyzed. The model is applied to investigate the damage mechanism of the patch-reinforced sample with a variety of geometrical and material parameters. The results provide a promising opportunity for quick design and optimization.

MODEL CONSTRUCTION TFP Yarn Mesh
Following the designed pattern in CAD software (Figures 2A, C), the yarn is placed and stitched onto the base material to produce the reinforcement ( Figures 2B, D). With the chemical vapor infiltration (CVI) technique, the TFP reinforcement could be consolidated within the SiC matrix to produce the composite reinforcement patches. Thereafter, the patch and the sample could be joined with the Ti 3 SiC 2 tapes and the in situ reaction TiC gradient layer via FAST (Zhou et al., 2018), while the open holes of the patch and base sample are aligned together to form the TFP patch reinforcement.
The TFP pattern is defined in CAD software in the format of "STandard for the Exchange of Product model data" (STP), from which the points' coordinates can be extracted. Based on these point coordinates, the yarn's 3D mesh is generated ( Figure 3A). Each element is defined with a local coordinate that follows the local fiber orientation ( Figure 3B). A typical feature is demonstrated in Figure 2D that the stitched yarn has a transitional cross-sectional area and inter-yarn space. The yarn width and space gradually reduce to the minimum, respectively, when the angle ϕ decreases from 90°to 0°. This feature can be well reproduced in the given mesh ( Figures 3A,C). During the change of the width, the cross-sectional area of the yarn is maintained constantly. Thus, the yarn thickness of the mesh varies inversely proportional to the yarn width.
The instant width of the TFP yarn can be calculated as where the angle θ ∈ (−π, π ] is the angle between the horizontal axis OA and the tangent line OB ( Figure 3A). O is an arbitrary point on the TFP path. ηis the width reduction factor from measurement. ω 0 is the maximum width of the yarn on the major axis, while η · ω 0 is the minimum width on the minor axis. Apparently, ω reaches ω 0 when θ equals zero or π and has the minimum when θ is ± 0.5π. To keep the cross-sectional area constant, the instant thickness at point O is where T 0 is the yarn thickness on the major axis.

TFP Patch Reinforcement
The geometric outer contour of the TFP patch reinforcement is a standard circle or a standard ellipse with a hole. In the traditional continuous model, the matrix mesh is obtained by Boolean subtraction of the TFP patch and the yarn mesh. However, the extremely tiny inter-yarn space between yarns makes the matrix mesh too small, which leads to poor mesh quality and high computational cost. An alternative is the "mesh superposition technique" (McGlockton et al., 2003;Yang and Cox, 2010), with which the yarn mesh is embedded into the mesh of the matrix material ( Figure 4). Actually, the "mesh superposition technique" is a method of multi-point constraints in Abaqus. Geometrically, the matrix mesh includes the yarn mesh. The matrix mesh is the host mesh, which is considered an independent model from the point of view of translational degrees of freedom (DOFs). The kinematic constraints are created between the nodes of the embedded and the host elements. The translational degrees of freedom of the embedded node are interpolated as the values of the DOFs of the surrounding host nodes.
The mesh superposition technique is an effective solution to the problem of composite meshing, as the volume mesh of the separate parts is much simpler and more efficient than that in a full model.

Sample With Open Hole
As shown in Figure 5A, the coupon is partitioned into specific sections (1) and (2), to apply the local refined structural mesh. The locally well-structured mesh serves the purpose of dealing with high stress concentration near the hole.
In Figure 5B, each ply of the composite is modeled. There are two ways to create the interface between neighboring plies or between the TFP patch and the base laminate sample: "cohesive element," and "surface-based cohesive behavior model" (SBCBM). Compared to the cohesive element, the SBCBM, which correlates the inter-ply traction and separation between plies, is physically realistic and numerically efficient. Usually, the inter-ply interface is so thin that the thickness, which is normally needed by cohesive elements, is negligibly small and hard to obtain. Moreover, the SBCBM does not require the contacting surface pair possessing a coincident mesh by sharing common nodes. Bonded by the SBCBM, an FE model of the single-side TFP-patched composite laminate sample is assembled together and demonstrated in Figure 6. This model contains 28,376 elements, which are linear solid elements C3D8 containing eight integration points, and 40,876 nodes in total. The model will be solved using the Abaqus explicit solver.

Traction-Separation Failure Model
As aforementioned, the ply interfaces are modeled using the SBCBM. The corresponding linear damage evolution is demonstrated in Figure 7. Instead of using strain and stress, this damage evolution model is a correlation between the traction force and the separation displacement. In the elastic zone (δ < δ 0 ), no material damage occurs. Once the traction reaches T f (A in Figure 7), the interfacial material starts to degrade by irreversibly losing its properties (A to B). After B, the interface loses all the properties and delamination takes place. The corresponding separation displacements to points A and B are δ 0 and δ f , respectively. Under a three-dimensional stress state, the separation δ stands for δ n ,δ s , and δ t , while the traction Trepresents t n , t s , and t t . The subscripts n, s, and t mean the normal direction, first shear direction, and second shear direction, respectively. The critical fracture energy G (G I , G II , and G III for Mode I, II, and III fracture energies, respectively) and cohesive strength (t 0 n , t 0 s , and t 0 t ) should be provided as inputs. The initial stiffness K is automatically selected by the system as contact penalty.
When the damage evolution is assumed to be linear as shown in Figure 7, the fracture energy G is represented as Frontiers in Materials | www.frontiersin.org August 2021 | Volume 8 | Article 701193  The fracture energies, G I , G II , andG III can be obtained through the double cantilever bending (DCB) test, end notched flexure (ENF) test, and edge crack torsion (ECT) test, respectively. In contrast, the cohesive strength T f is difficult to retain directly. They can be, nevertheless, reversely calculated by fitting the simulation results to the experimental data (Diehl, 2008).

Damage Model
The damage initiation is triggered by the quadratic nominal stress criterion, as shown in the following equation: The Macaulay brackets 〈 · 〉 imply that compressive stress will not cause delamination. The mixed mode damage is governed by Benzeggagh-Kenane mixed mode behavior (Benzeggagh and Kenane, 1996): where ξ is the material constant determined by experiments. G is the crucial fracture energy obtained from the experiment, too. It is very convenient when the fracture energy along the first and second shear directions is the same (G s G t ). Then, Eq. 5 is rewritten as Considering the interfacial properties for epoxy matrix composites and SiC matrix composites are different, the different input data are listed in Table 1.

Composite Ply Model
Continuum Damage Model Muñoz et al. (2015) proposed a continuum damage model to describe the progressive damage of the composite ply. Five damage functions that represent fiber loading directions are defined in Table 2.
In Table 2, the symbols φ, E, S, ε, and τ stand for the damage activation function, Young's modulus, the strength, the strain, and the shear stress, respectively. The superscripts t and c mean tension and compression, respectively. In the post-damage phase, the compliance matrix was defined as follows (Maimí et al., 2007): where C is the compliance matrix; d N (N 1, 2, 6) are the damage variables; and μ is Poisson's ratio.
In-plane shear failure φ 6 |τ12 | S12 corresponding to loadings applied to the material in the local frame. N 1+, 1− represents tensile failure and compressive failure in direction 1. N 2+, 2− represents tensile failure and compressive failure in direction 2. N 6 represents in-plane shear failure. The damage activation functions have the following expression: where r N (N 1+, 1−, 2+, 2−, 6) define the elastic boundary in the post-damage phase, before further damage takes place. Maimí et al. (2007) related r N to the damage variables d N (N 1+, 1−, 2+, 2−, 6) by exponential law: where A N (N 1+, 1−, 2+, 2−, 6) are material-dependent coefficients, defined as where E N , X N , and G N (N 1+, 1−, 2+, 2−, 6) are the elastic moduli, strengths, and fracture energies. l p stands for the characteristic length of a finite element. The parameter A N , according to Bazant's model (Bažant and Oh, 1983), guarantees that the damage propagation and energy release are independent of the mesh size (Maimí et al., 2007). The fracture energies G N control element deletion.

Input of Consolidated TFP Yarn and Woven Laminates
The consolidated TFP yarn can be treated as a homogeneous and transverse isotropic material. The strength components, stiffness components, and Poisson's ratios of the impregnated yarns are calculated ( Table 3) Table 4).

Strength Prediction
A comprehensive study on strength and fracture criteria in fabric composite samples with a circular hole was carried out by virtual test. The detailed elastic properties of the woven SiC f /SiC composite samples that contain [0, 45, −45, 0] plies are listed in Table 5.
To run a calculation, the necessary material properties as input for the model are already provided in Material Models. The dimensions of the notched sample are described in Figure 8.
The radius of the hole varies from 2 to 10 mm, in order to investigate the dependence of failure on hole size.
The calculated strengths of the drilled sample with different hole radii (2-10 mm) are listed in Table 6. The residual strength of the laminate with an open-hole radius of 10 mm is only 36.32% of the strength of the un-notched laminate. In this case, it is necessary to adopt reinforcement methods for the composite laminate with an open hole or notch.

Stress Analysis
In Figure 9, the tension and shear contours immediately before the damage onset for a [0, 45, −45, 0] laminate, which has a hole radius of R 10 mm, are demonstrated. The stress concentration factor, which is defined as the ratio of the maximum tensile stress value on the edge of the hole (Figures 9A-D) to the value of external tensile stress (90 MPa), is 4.4 and 2.2 for 0 o and 45 o plies, respectively. The external tensile stress means loading force divided by the material cross-section without the hole. The stress concentration factor indicates that the fiber rupture will take place firstly in 0 o plies and then extend to ±45 o plies ( Figure 10). Under the on-axis loading, the ±45 o plies ( Figures 9F, G) bear the shear loadings effectively, three times those of 0 o plies (legends in Figures 9E, H).
When element deletion is activated, the aforementioned damage initiation and propagation can be clearly demonstrated, as shown in Figure 10. The cracks (due to fiber rupture, damage variable d 1 0.99) start in 0 o plies when the sample is applied with external stress that equals 122 MPa ( Figures 10A, B). At this moment, except the shear damage ( Figure 10C), there is no fiber rupture occurring in ±45 o plies ( Figure 10B). The cracks propagate in 0 o plies till the applied stress reaches 146 MPa, and suddenly, the ±45 o plies are torn while applied stress snaps back to 118 MPa ( Figure 10D).

Delamination
It is obvious that stress concentration "amplifies" the loading at the hole circumference. Furthermore, due to high material discontinuity at the interfaces between the 0 o ply and the 45 o ply, the in-plane shear stress ( Figure 11A) and out-of-plane shear stress ( Figures 11B, C) help to form the mild interfacial damage ( Figure 11D) in an early stage when the external stress is only 45 MPa. Remarkable delamination presents when the 0 o laminate   (Abdi et al., 2016).     Figure 11F). The delamination between ±45 o plies, where the in-plane shear and out-of-plane shear are neglectable, is milder than that between the 0 o ply and the 45 o ply ( Figures  11D, E).

Numerical Tests on Notched Sample With Patch Reinforcement
With the help of the FE model, virtual test campaign is carried out for notched samples reinforced by circular annular (CA) and elliptical annular (EA) patches, based on the patching parameters listed in Table 7.
The dimensions of the virtual test samples are delineated in Figure 12A, where the patch matrix material is set to invisible. The geometrical parameters of the samples are as follows: L 150 mm, W 60 mm, T 2 mm, r 10 mm, D 47 mm, a 41 mm, and b 80 mm, where D, a, and b are the outer diameter of the CA patch and outer major axis and minor axis of the EA patch, respectively. The single-side reinforced (SSR) and doubleside reinforced (DSR) samples are sketched in Figure 12B.

Study on Patching Parameters
In Figures 13A, B, the ultimate applied stress and apparent moduli of the CA patch-reinforced sample, as well as those of the undrilled sample and drilled sample, are delineated. The abbreviations OP, NP, WS, WD, YS, and YD represent the original sample without hole, sample with a hole, and drilled sample with woven SSR, with woven DSR, with TFP SSR, and with TFP DSR, respectively. From Figures 13A, B, two conclusions can be drawn: 1) The TFP patch exhibits overwhelming advantages in terms of enhancing the strength. In contrast to the woven fabric, the TFP yarn seems to be more efficient in resisting the maximum principal stress. 2) DSR samples show higher strengths and moduli than SSR samples. SSR ( Figure 13C) cannot omit the stress concentration on the unpatched side ( Figure 13E), where the damage and crack occur earlier. In contrast, the DSR samples ( Figure 13D), where the same total amount of reinforcement (evenly separated on double sides) as in SSR is used, can ideally balance the load, eliminate the stress concentration, and delay the damage onset ( Figure 13F).
In order to analyze the stress distribution around the hole, S22 contours are sampled clockwise from the front and back sides of the sample along the hole edges ( Figure 14). For the SSR sample, the back edge has patch reinforcement, while the front has not (Figures 13C, E). It is worth noting that the stress (solidsquare-symboled curve in Figure 14) obtained from the front edge in Figures 13C, E shows larger oscillation and higher peak values of 400 MPa, while the edge with patch reinforcements reads a value of about 350 MPa. This higher local stress leads to earlier damage onset and well explains the aforementioned conclusion ii).

Study on Patch Shape
In Figures 15A, B, the ultimate applied stress and apparent moduli for CA and EA patch-reinforced samples, as well as those      Figure 15A). Comparing the stress curves from the patched side, the curve (solid triangles) of the EA-patched sample is much flatter than the CA counterpart (solid circles), which suggests the EA patch has more efficiency in reducing the local stress concentration.
The tensile stress along the longitudinal axis of the notched coupons and patches (vertical red lines sketched on the sample and patches) on the contacted surfaces of the EA and CA patch-reinforced samples is extracted along with and without the nominal elastic energy along the axes (Figure 17). Stress curves are extracted from the patch side of the SSR sample. When applied with the same external strain, the former stress curve shows higher values in zones I and II (outside of the patching area). The local stress reaches the maximum of 200 MPa on the boundary between zones I and III, as well as that between zones II and III. This implies that the damage will take place on the drilled sample next to the edge of the patch when the patch is strong enough and perfectly bonded.
The total elastic energy stored in the sample can be represented as where G e is the total elastic energy and σ, ε, V, E, and E are the local stress, the external strain, the total volume of the sample, Young's modulus, and the apparent Young's modulus, respectively. Using the above equation, the apparent tensile modulus can be evaluated. To explain the reason why the elliptical TFP single-side reinforced sample possesses a higher apparent modulus than the circular TFP single-side reinforced sample (229.4 vs. 211.5 GPa) ( Figure 15B), a concept of nominal elastic energy can be introduced as where L and w are the length of the sample and the width of the elements lying on the loading sampling path. ne G and ne E are the nominal elastic energy stored in the elements along the sampling path shown in Figure 17 and the nominal elastic modulus, respectively. The tensile stress σ is the function of the nodal coordinate on the loading path, as demonstrated in Figure 17.
Hence, the ratio of the nominal moduli, R n , for elliptical and circular TFP single-side reinforced samples reads where subscripts and superscripts "Ellip," "Circ," "S," and "P" stand for elliptical and circular patches and the drilled sample and patch, respectively. The integrals in the above equation are calculated and presented in Figure 17, too. It is apparent that  R n is larger than 1.0. The factor R n reveals the reason why the apparent modulus of the elliptical TFP single-side reinforced sample (229.4 GPa) is larger than that of the circular one (211.5 GPa) ( Figure 15B). By applying the same deformation, more elastic energy is stored in the elliptical patch and zones I and II of the drilled coupon for the elliptical TFP single-side reinforced sample than for the counterpart.

CONCLUSION
In this work, an advanced progressive damage FE model for the notched ceramic matrix composite laminate reinforced by the special TFP composite patch has been established. An energybased composite progressive damage model that is independent of the mesh size is applied. This model also accounts for the TFP yarn details that are inherited from TFP path design. The embedded element concept is used to successfully avoid the meshing difficulties. The inter-ply "glue" layer is simulated using the surface-based contact method, and the delamination due to interfacial material discontinuity and damage can be well reproduced and analyzed. A virtual test campaign was carried out based on a variety of geometrical and material parameters. Based on the stress analysis, the related damage mechanism was revealed. From the FE analysis, the following conclusions can be drawn: 1) The optimized TFP patch is more efficient in minimizing the stress concentration around the hole edges than the woven fabric.
2) The double-side reinforcement, compared to the single-side counterpart, reduces the stress concentration introduced by the geometrical asymmetry on the unpatched side and thus highly increases the modulus and strength.
3) The elliptical patch has shown better mechanical performance. The apparent strength reaches 79% of the un-notched sample, while the modulus exceeds that of the original sample.
In summary, the proposed FE model well simulated the damage behavior of the patch-reinforced ceramic matrix composite sample. Through the FE analysis, the double-side elliptical sample reinforced by the TFP yarn is highly recommended. In the future work, nevertheless, to correctly predict the orientation of crack propagation, especially that in 45 o plies, the enhanced FE method or the orientated mesh should be utilized.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, and further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
GZ carried out investigation, ran the software, contributed to testing code components, and wrote the first draft of the manuscript. JT performed the methodology. GZ, JT, JW, YC, and YF curated the data. JW, YC, and YF were responsible for modeling and simulation. YC, SX, XJ, SZ, and JX contributed to supervision and funding acquisition. SL revised the first draft of the manuscript. SZ and JX contributed to project administration. GZ and JX contributed to supporting algorithms and reviewed and edited the paper.