Predicting effective fracture toughness of ZrB2-based ultra-high temperature ceramics by phase-field modeling

Abstract The effective fracture toughness (EFT) of ZrB2-C ceramics with different engineered microarchitectures was numerically evaluated by phase-field modeling. To verify the model, fibrous monoliths (elongated hexagonal ZrB2-rich cells in a continuous C-rich matrix) with different volume fractions of a C-rich phase were considered. Architectures containing 10 and 30 vol% of C-rich phase showed EFT values about 42% more than that of pure ZrB2. Increasing the C-rich phase to 50 vol%, dropped toughness significantly, which is in agreement with the experimental results. Replacing hexagonal cells with cylindrical, triangular, or square cells of the same cross-sectional area changed the toughening mechanism and EFT. The orientation of the interface between the soft and hard phases with respect to the crack orientation also affected the energy required for crack propagation, and in some cases resulted in a higher EFT (even up to 70% of pure ZrB2 fracture toughness) either by suppressing uniform crack propagation or making crack cranking. Results not only show that the model can predict fracture toughness but also provide insight to improve toughness by engineering different microarchitectures.

have applications in extreme environments because of their high melting temperature, excellent strength, a relatively good chemical stability and thermal shock resistance. Examples of UHTCs include borides, nitrides and carbides of metals from groups IV and V in the periodic table [1,2]. Diboride materials especially ZrB 2 and HfB 2 are of the most interest because they have the best oxidation resistance [3,4].
The main weaknesses of UHTCs are their brittle fracture behavior and low damage tolerance, which have limited their applications. Cook and Gordon [5] were the first to introduce the idea that it is possible to control crack propagation in a brittle material by considering particular microstructural features that can change the crack path. Later, Clegg et al. [6] showed that by separating strong phase layers with weak interphases, a brittle composite ceramic can fail in a non-brittle manner by deflecting the crack through different phases. Several attempts have been made to process and characterize mechanical and thermal properties of different engineered architectures for ZrB 2based ceramics [2,[7][8][9][10][11][12][13][14]. For example, Fig. 1 shows load-displacement curves for a single phase ZrB 2 ceramic and a laminated ZrB 2 -C composite; ZrB 2 by itself is brittle, and exhibits no visible inelastic work of fracture after reaching the breaking load. On the other hand, the laminate composite showed non-brittle fracture behavior through crack deflection in a non-elastic form. Similar studies have been done on bioinspired microstructure design of laminates to increase their toughness by providing preferential paths for propagating cracks [15][16][17][18][19]. Chen et al. [20] fabricated Al 2 O 3 /reduced graphene oxide (rGO) fibrous monolithic ceramics with bamboo-like structures that had 475% and 1075% higher fracture toughness and work of fracture than those of the mono-lithic Al 2 O 3 ceramics. In a different study, TaC-based/graphite fibrous monolithic ceramics were investigated by Shahedifar et al. [21] to study the fiber core/shell volume ratio and fiber orientation on the fracture toughness and the work of fracture. Experimental efforts for developing ceramics with engineered microarchitectures have laid out the initial path for improving the fracture toughness of UHTCs; however, these were mostly based on trial and error experiments. Another example is the designed microstructure by Parthasarathy et al. [22] for a nonconventional fiber reinforcement, which is shown in Fig. 2. They showed that a combination of bone-shaped fibers with a surface compression coating, or multilayer compositions, can be used to optimize strength and toughness via an engineered microarchitecture in ceramics.
There are a few methods that analytically predict the fracture behavior of microarchitecture ceramics; however, they can be used only for simple microstructures. Zhang et al. [23] studied fracture of unidirectional nanocomposite structures with parallel staggered platelet reinforcements, and they analytically calculated four dimensionless parameters associated with platelet distributions to show the effect of volume fraction, orientation, and the distribution of platelets on the fracture behavior. In another study, Leguillon et al. [24] used a coupled stress-energy criterion to analytically predict the initiation and propagation of surface cracks in ceramic laminates of Al 2 O 3 /monoclinic ZrO 2 and Al 2 O 3 /tetragonal ZrO 2 under thermo-mechanical loading. Begley et al. [25] used an analytical micromechanical analysis for composites comprising elastic platelets (bricks) bonded together with thin elastic perfectly plastic layers (mortar) under uniaxial loading; they introduced a closed-form solution for the spatial variation of displacements in the bricks as a function of constituent properties, which can be used to calculate the effective properties of the composite, like the work of fracture. Separately, Wang et al. [26] developed a temperature dependent fracture strength model for the laminated ultra-high temperature ZrB 2 -BN ceramic composites; this model was based on the Griffith energy criterion and the concept of energy storage capacity. All these analytical models were developed for laminate microstructures and cannot be used for complex microstructures.
Numerical modeling and simulations of the fracture process, validated by experiments, can be used as powerful tools to predict the failure mechanism and fracture properties of composite ceramics and design microstructure-engineered ceramics with superior fracture properties. Such computational simulations can also reduce the number of trial and error experiments in the design process of microarchitectured composites. There are some computational studies based on classical continuum mechanics for investigating the fracture of composite materials. Delamination in composites is a commonly studied problem. In the majority of cases, delamination is simulated using the cohesive interface model [27][28][29][30][31]. The crack path needs to be predefined in cohesive interface model. Therefore, to study delamination with this model, the crack was restricted just to the interfaces. Extended finite element (XFEM) is another numerical method, based on the Finite Element Method (FEM), which is especially designed for treating discontinuities without predefined their path. The main advantage of XFEM related to FEM is that the finite element mesh does not need to be updated to track the crack path. However, XFEM has difficulties in nonlinear heterogeneous systems. Yan and Park [32] used XFEM [33,34] to model near interface crack growth in a ceramic-metal-ceramic laminate. Although the calculated crack path was similar to the experiments, the crack growth did   2. A schematic designed microstructure for a composite ceramics with enhanced strength and fracture toughness [22]. The matrix and the bone-shaped fibers were made of HfN or HfC 0.67 . not deflect into different layers and remained only at the interfaces between different layers. Also, the fracture toughness of the composite was not determined.
Recently, the phase-field (PF) method has emerged as a powerful and unique meso-scale modeling tool for studying crack formation and propagation in nano and microstructures. Utilizing PF models in simulating crack propagation have some major advantages over the conventional methods like FEM, XFEM, and cohesive zone model. For example, the crack path is not predefined in a PF model, and the model will be truly a predictive tool for crack propagation. PF models can be easily implemented in complex cracking situations including crack merging and branching, even in three-dimensional problems. There are a few PF models that studied crack propagation in microstructure-engineered ceramics. However, all were in dimensionless form [35][36][37][38] or they used a hypothetical set of material properties [39]. For example, Khaderi et al. [37] used a PF model to explore crack propagation and calculate critical energy release rate in a layered mineral/organic composites as a function of elastic modulus mismatch and the thickness of the organic layer. A modified PF model developed by Emdadi et al. [40] predicted the crack path quantitatively for composite ceramics; this model was not applied to study different microstructures, and also the fracture toughness of composites was not determined. To the best of our knowledge, there is no attempt in using a PF model to evaluate the fracture toughness of engineered ceramics with real properties of phases.
This study introduces a quantitative computational framework for evaluating an effective fracture toughness (EFT) for engineered ceramics that can be used not only to study crack propagation in ceramic architectures, but also to design and optimize the microarchitectures and potentially predict the properties of constituent phases that would maximize the fracture toughness and damage tolerance of composite ceramics. This work is focused on demonstrating the capability of the PF computational model in prediction of EFT of composite ceramics, and some experiments are performed to validate the crack path and the trend of the fracture toughness in the fibrous monolith (FM) microarchitectures (FMs with elongated hexagonal ZrB 2 -rich cells in a C-rich matrix) with different volume fractions of a C-rich phase. After verifying the computational model, the effect of different microarchitectures (cylindrical, triangular and square cells) on EFT of ZrB 2 -C composites is investigated.

Experimental procedure
Cellular ZrB 2 -carbon FMs were produced via a thermoplastic forming method [41]. The raw materials and compositions used for the ZrB 2 -carbon FMs are shown in Table 1. To produce the FM's, ZrB 2 , B 4 C, and WC powders were first attrition milled together (Model 01-HD, Union Process, Akron, OH) for 2 h using acetone and tungsten carbide (WC) milling media. After milling, the powder was dried via rotary evaporation (Model Rotavapor R-124, Buchi, Flawil, Germany). The graphite powder was ball milled for 12 h using acetone and WC milling media to break down agglomerates and large flakes of graphite. The graphite slurry was wet sieved through a 270 mesh sieve to remove any remaining large agglomerates or flakes. The graphite slurry was then dried via rotary evaporation. The cell and cell boundary phases were batched according to Table 1 by blending the powders with a thermoplastic polymer and plasticizer using a heated high shear mixer (C.W. Brabender, South Hackensack, NJ). The cell phase was pressed into cylindrical rods and the cell boundary phase was pressed into semicircular shells to be laminated around the cell material. The finished feedrods were~22 mm in diameter. The ratio of cell to cell boundary varied according to the desired composition. The feedrods were then extruded directly to 300 μm diameter filament using a screw-driven extruder and a heated spinneret. The batching and extrusion process are described in more detail elsewhere [8,11,42].
After extrusion, the filament was wound around a cylinder with each subsequent wrap being tight against the one before. The filaments were glued together with a spray adhesive. This produced a uniaxially aligned ribbon of filaments when the wrapped filament was split and removed from the cylinder. These ribbons were then cut into rectangles 45 mm × 35 mm to match the size of the hot-pressing die. Layers of the filament ribbon were loaded into a rectangular graphite die lined with graphite paper that was coated with boron nitride. Enough layers were used to produce a sintered billet~5 mm thick. Once loaded, the layers were laminated together in the graphite die at 130°C to form a solid billet. The entire die was loaded into an atmosphere-controlled furnace and heated at 5°C/h to 600°C under a gas mixture of 90% argon -10% hydrogen to remove the organics. Following binder removal, the die was transferred to a resistively heated graphite hot press (Thermal Technology Inc., Model HP20-3060, Santa Rosa, CA) for sintering. The ramp rate during heating was 25°C/min. Vacuum holds were employed at 1450°C and 1650°C to remove oxide species from the surfaces of powder particles [43]. At 1650°C, the atmosphere was changed to argon and a pressure of 32 MPa was applied. The furnace was then heated to 2000°C and held until densification ceased. Densification was monitored by the travel of the pressing ram in the hot press. The furnace power was then turned off, and the specimen was allowed to cool. After hot pressing, the densified billets were machined into mechanical test bars. Bars were cut from the billet and ground to a final dimension of 3 mm × 4 mm × 45 mm using an automated surface grinder (Chevalier, FSG-3A818, Santa Fe Springs, CA) and diamond abrasives. These bars were cut such that the cells of the specimen architecture were aligned with the long axis of the bars (Fig. 3). Grinding to the final dimension was accomplished using a 400-grit diamond abrasive wheel. The composites examined for this study consisted of cell to cell boundary ratios of 95:5, 90:10, 82.5:17.5, 70:30, and 50:50 by volume. For each composition, 5 bars were tested for flexural and fracture toughness according to ASTM C1161 [44] and ASTM C1421 [45], respectively.

Computational model
A recently modified PF model [40] based on the regularized formulation of Griffith's fracture theory [46] was utilized in this work. In this model, the entire quasi-static process of crack initiation, propagation and branching is governed by the minimization of an energy functional E k (u, ϕ c ). This functional is the variational formulation [47,48] of Griffith's theory where u is the displacement vector, and ϕ c is the scalar crack PF parameter describing a smooth transition between the unbroken (ϕ c = 1) and broken (ϕ c = 0) state of the material. The energy functional is in the form of Eq. (1) in which the crack is considered to be a phase with the evolution equation of Eq. (2).
Eq. (3) shows the mechanical equilibrium equations: In the above equations, F e is the elastic energy density, ε(u) is the strain tensor, G c is the critical energy release rate or crack surface energy in Griffith's theory for admissible crack set Γ ⊂ Ω, t is the external traction applied on the boundary of ∂Ω, and C is the elastic stiffness tensor. k is a positive regularization parameter to regulate the fracture zone, and η k is a small (related to k) residual stiffness to avoid singularity in the first part of the energy in Eq. (1) in fully fractured regions of the domain. M is the mobility of the crack.
A ⁎ in Eq. (1) and B ⁎ in Eq. (3) are the correction parameters in the total free energy functional and the mechanical equilibrium equation which were defined in the modified PF model [40]. These parameters were defined to consider the effect of material strength on crack nucleation and propagation independent of the regularization parameter, and also to ensure that the maximum stress in front of the crack tip is equal to its counterpart predicted by classical linear elastic fracture mechanics (CLEFM). It should be noted that A ⁎ and B ⁎ are only applied in the diffusive crack area where ϕ c b 1, and they are unity in intact regions: σ s is the material strength, and H[ϕ c − 1] is the Heaviside step function. More information on how these two correction parameters, A and B, were defined can be found in Ref. [40]. To prevent any formed crack from healing, irreversibility of the crack PF variable was ensured by introducing a local strain-history field of the maximum strain energy density [49] by: To calculate the EFT, the approach of using a moving boundary condition proposed by Hossein et al. [38] was used. In this method, a steady propagation of a mode-I crack opening displacement was applied as the boundary condition, then the crack was allowed to propagate based on the free energy minimization (Eq. (1)). For a mode-I crack, the displacement field in front of the crack tip is [50]: where r = |r| is the magnitude of the vector r, a vector with its origin at the crack tip that extends to the point of interest in front of the crack, θ is the angle between the orientation of the crack and r (−180°b θ b 180°) with positive values for counter-clockwise rotation, K Ã I is the mode-I stress intensity factor, β = 2(1 − υ) for plane strain and β = 2 (1 + υ) for plane stress conditions, and G is the shear modulus. The time-dependent steady crack opening displacement is shown in Eq. (10), and this displacement is moving with a uniform velocity of v. An effective elastic modulus is considered for the heterogeneous system (E eff ), using the rule of mixtures, Eq. (11). u ¼ u Ã ðx−vt; yÞ on ∂Ω; ð10Þ where E i is the elastic modulus and A i is the total surface area of phase i. The different considered microstructures were kept inside a homogenous elastic system with the elastic modulus of E eff to ensure that the J-integral can be used to evaluate EFT for the heterogeneous system [38]. The energy release rate can be calculated through the pathindependent J-integral [50]: in which Γ is a closed path around the crack. The calculated EFT using Eq. (12) was modified according to the mesh size, h [48], where All simulations in this study were performed using the measured properties and the dimensional form of the governing equations for ZrB 2 -C ceramics with engineered microarchitectures to be able to compare the results to the available experimental results. Since many other PF models for crack propagation have utilized unitless or scaled simulations, the non-dimensional forms of the equations for a two-phase heterogeneous system is provided in Appendix A. The plane strain condition was considered, and the coupled equations, including crack evolution and mechanical equilibrium equations, were solved in a finite element framework with linear Lagrangian elements using the mathematics module COMSOL Multiphysics [51]. An adaptive time step was used to solve nonlinear equations with a maximum time-step size of 0.5 s. All the simulations were performed using a desktop computer with two Xeon Phi processors (E5-2687W-total of 40 CPU cores) and 128 GB RAM. The material parameters for ZrB 2 , and C are summarized in Table 2. It was assumed that h (mesh size) is 0.4k [40], and v = 0.1 μm/s.

Experimental results
The extruded filaments were laid up in a regular array, Fig. 3, and that regularity translated to the densified microstructure, Fig. 4. The compressed appearance of the cellular architecture is due to the lamination and densification processes. The initial consolidation step, lamination, compressed the architecture to fill in the empty space between filaments and then the subsequent densification by hot pressing further compressed the architecture in the same direction. These composites were produced such that the basic unit cell size of the composite, one cell plus its surrounding cell boundary, did not change with composition, only the relative ratio of cell to cell boundary changed.  [41,54]. As compared to the monolithic cell material without the cell boundary, FMs typically exhibit a reduction in strength, but an increase in crack deflection and, therefore, an increase in toughness and work of fracture. The work of fracture of these specimens was calculated based only on the inelastic portion of the stress-displacement curve; using the area under the curve after the first indication of crack initiation. The area under the curve results in units of Joules which can then be divided by twice the cross-sectional area of the bar (twice the 3 mm × 4 mm cross section) to obtain work of fracture in units of J/m 2 (Fig. 5). In this case, the 90:10 composition exhibited the highest work of fracture at 5720 J/m 2 . After initial attempts to measure the fracture toughness (K IC ) were unsuccessful, the work of fracture was selected as the parameter to measure fracture properties. Measurements were attempted according to ASTM C1421 [45], which describes the chevron notch toughness measurement technique. Specimens were notched and tested according to the standard. However, valid tests for K IC require that the crack be driven through the notch in mode-I. For all of the FMs tested, the crack front deflected   out of plane and followed the cell boundary phase as illustrated in Fig. 3. This large-scale crack deflection is indicative of improved toughness and graceful failure, but invalidated the K IC measurements.

Computational results
The initial step was to numerically calculate the critical crack surface energy for a homogenous ZrB 2 rectangular specimen to show how the value of K Ã I affected the calculated critical crack surface energy. Then, the crack path and EFT were studied in ZrB 2 -C ceramics with FM architectures and different C-rich contents. Other cases with different microarchitectures, but the same phase fractions as FM cases, were also investigated to determine the effect of microarchitecture on EFT.

Calculated crack surface energy for homogeneous ZrB 2
A 1 mm × 2 mm rectangular domain of homogeneous ZrB 2 including an initial crack with the length of 2 mm was considered. In this example, the objective was to verify the dependence of the calculated critical crack surface energy of ZrB 2 on the value of K Ã I in the applied mode-I crack opening displacement. The plane-strain condition was assumed in this example with k = 0.01 × 1 mm [55]. Three different values of K Ã I = 1.7, 2.8, and 5.5 MPa ffiffiffiffi ffi m p were considered. Fig. 6 shows the geometry and the propagated crack resulting from a moving boundary condition for a mode-I crack opening displacement.
The calculated J-integrals versus time for the homogeneous ZrB 2 sample were presented in Fig. 7. The value of the J-integral at the beginning of the simulation should be equal to the assumed G Ã IC in the applied displacement. This implies that initiation of a crack in a material with a surface energy of G Ã IC requires the same energy release rate (same Jintegral). However, after the crack initiates, the energy release rate should be equal to the critical crack surface energy of the material for crack propagation. /E). For these calculations, ϕ c = 1 meant that the crack has not initiated yet. Since the calculated J-integrals of 7.3 J/m 2 (for K Ã I = 1.7 MPa) and 15 J/m 2 (for K Ã I = 2.8 MPa) were lower than the critical crack surface energy of ZrB 2 (24 J/m 2 ), the crack did not grow. By moving the position of the crack tip in the applied displacement to the right, the J-integral increased, and, simultaneously, ϕ c decreased from 1 to 0. The crack propagated steadily until the J-integral reached the ZrB 2 critical crack surface energy (24 J/m 2 ) and ϕ c dropped to 0, indicating that the crack had propagated through the material. It is evident from Fig. 7 that steady-state crack propagation started earlier under K Ã , which is expected to occur experimentally, as larger applied displacements make the crack propagates earlier than a smaller displacement. The calculated J-integral dropped quickly for K Ã I = 5.5MPa since it was initially higher than the G c of ZrB 2 , and reached the critical crack surface energy of ZrB 2 . The initial value for the J-integral reflects the considered K Ã I in the applied displacement equation, and its value is equivalent to the critical crack surface energy of the material after steady-state crack propagation. As a result, the calculated crack surface energy after steady-state crack propagation is independent of the considered K Ã I in the applied displacement equation. For all of the multiphase simulations that follow, it was assumed K Ã I ¼ 2:5 MPa ffiffiffiffi ffi m p , and the average of maximum values of the calculated J-integral for the middle 50% of the length of each specimen was considered to be the EFT [38]; this will ensure that the boundary conditions have minimal effects on the value of EFT.

Effective fracture toughness for ZrB 2 -C fibrous monoliths
An FM structure consists of a major (80-90 vol%) phase that is brittle, which constitutes "strong cells". The cells are surrounded by a thin continuous "cell boundary" phase (10-20 vol%) of a weaker material [13]. In this study, it was assumed that the hexagonal ZrB 2 -rich cells (80 vol% ZrB 2 + 20 vol% C) were covered with a thin C-rich layer (20 vol% ZrB 2 + 80 vol% C). A linear rule of mixtures was used to calculate the effective elastic modulus and effective strength of each phase. According to Table 2, the G c values for C and ZrB 2 are very close, so it was assumed that G c for the C-rich phase, and the ZrB 2 -rich phase were the same as pure C and pure ZrB 2 . The dimension of the ZrB 2 -rich cells and the thickness of the C-rich layer in a 10 vol% C-rich FM ceramics are shown in Fig. 8. The FM ceramics are surrounded by a 0.05 mm thick layer of homogeneous material with the material properties equal to the effective material properties of the heterogeneous system using Eq. (11). The total size of the specimen was 0.7 mm × 2.2 mm. An initial crack was considered in the homogeneous material. k was assumed to be 1% of the specimen height.  In Fig. 9, the crack path and the calculated J-integral for a FM ceramic containing 10 vol% of the C-rich phase are presented. The position of the crack tip is related to its counterpart from the J-integral diagram and is shown in Fig. 9. The value of the J-integral at the beginning of the simulation was associated with the onset of crack propagation in the C-rich phase (~15 J/m 2 ) where the value of ϕ c is starting to decrease from one (in intact C-rich phase). As the value of ϕ c decreased, the Jintegral value increased until reaching the first peak in the J-integral diagram; at the time that ϕ c reached the zero value, the J-integral suddenly fell to the value that is related to crack propagation in the C-rich phase (~22 J/m 2 ). This means as soon as crack nucleation occurs, the crack will propagate in the C-rich phase until encountering the interface between the C-rich cell boundary and the ZrB 2 -rich cell. Then, the Jintegral increases when the crack is deflected at a 60-degree angle and the crack propagates along the edge of the ZrB 2 -rich hexagon. Toughening results from elastic and strength heterogeneities, as well as the microarchitecture. As soon as the crack tip reaches the corner of the hexagon, the crack has a straight path towards the C-rich phase parallel to the x-direction, and the J-integral value decreases to its value in the C-rich phase (~22 J/m 2 ). With the assumption of v = 10μm/s, the applied displacement takes~210 s to move from x = 0 mm (position of the initial crack tip) all the way to the right edge of the specimen. As mentioned before, the average of maximum values of the J-integral for the middle 50% length of the sample (when time is between~50 s and~160 s) was considered to be the EFT. For the FM ceramics with 10 vol% of the C-rich phase, EFT was~31.3 J/m 2 , which is 30% higher than the value of G c in the ZrB 2 -rich phase. Fig. 10(a) shows the effect of different vol% of the C-rich phase on EFT in ZrB 2 -C FM ceramics; for brevity, vol% C is used in the figures.
Increasing the soft C-rich layer content from 10 vol% to 30 vol% does not have a significant effect on the calculated J-integral or the resulting EFT. On the other hand, by increasing the C content to 50 vol%, EFT drops to the fracture toughness of the C-rich phase (~22 J/m 2 ). This behavior is in agreement with the experimental measurements for the flexure stress-deflection response in fourpoint bending tests of ZrB 2 -C FM ceramics as depicted in Fig. 4. The drop after the first peak in the stress-deflection plot was related to the crack initiation. All three microstructures in Fig. 4 showed load retention after the initial crack, which was an indication of nonbrittle fracture behavior for the ceramics. The onset of the initial crack in the FM with 50 vol% C occurred at the lowest flexure stress for the FMs, and therefore it had the lowest non-elastic work of fracture, which represented the lowest fracture toughness. The crack path for each case is presented in Fig. 10(b)-(e). The cracks in these microarchitectured ceramics pass solely through the C-rich phase, which is consistent with the experimental observations ( Fig. 10(f)). As mentioned in Section 4.1, no standard test has been established to measure the EFT of microarchitectured ceramics because of the mixed failure modes. The flexural stress-deflection data from four-point bending tests [44] shows non-brittle behavior of the ceramics and an increased fracture toughness compared to pure-ZrB 2 based on the presence of inelastic work of fracture. The numerical results for ZrB 2 -C FM ceramics confirmed that the model can predict the increased EFT of the microarchitecture ceramics with different vol% of the C-rich phase, as was observed experimentally. After verifying the model, the goal in the next section is to determine how different microarchitectures affect the EFT of the ZrB 2 -C ceramics.

Effective fracture toughness for different microarchitectured ceramics
The effect of different microarchitectures on EFT was investigated by replacing hexagonal ZrB 2 -rich cells with triangular, square, or cylindrical cells with the same surface area. Fig. 11(a) shows that the dependency of EFT on volume contents of the C-rich phase in triangle-cell FMs is similar to the hexagonal-cell FM results. After the crack kinks around the triangular ZrB 2 cell, it propagates through the C-rich shell at the edge of triangles (30°angle) (Fig. 11(b)-(d)). Crack trapping at the inverse-k points in the microstructure for FMs containing 10, 15, and 30 vol% of the C-rich phase (for instance, area A marked by a circle in Fig. 11(c)) results in toughening. Then, the crack is deflected in a ZrB 2 cell and follows an inclined path down (instead of a straight path along the x-axis through the ZrB 2 phase) until reaching the C-rich phase again. This means that the total energy was lower when the crack deflected than it propagated straight through the cell. From Fig. 11(a), EFT in microstructure with 30 vol% of C-phase is 34 J/m 2 , which is 41.6% more than G c in pure ZrB 2 . In Fig. 11(e), for the case with 50 vol% of the C-phase, the crack propagates along an almost straight path through the C-rich phase, hence the heterogeneous system EFT is about the same as that of the pure C-rich phase, which is~22 J/m 2 . The numerical results showed that having 10 vol% C-rich phase in a triangular-cell ZrB 2 -C FM ceramic can increase EFT by about 42% compared to that of pure ZrB 2 . Further, increasing the content of the C-rich phase to 30 vol% did not change EFT. Hence, for triangular cells, the EFT does not depend on the amount of the weak phase up to 30 vol%.
We also considered a microarchitecture with square ZrB 2 -rich cells in a C-rich matrix, and the results are presented in Fig. 12(a)-(e). Fig. 12(a) shows the calculated J-integral versus time as the crack propagates through samples containing 10, 30, and 50 vol% of the C-rich phase, as well as a case with 10 vol% of the C-rich phase, but where the initial crack was offset from the center of the domain  and was placed in front of a C-rich layer that is aligned between two ZrB 2 -rich cells, Fig. 12(e). In this case, the crack propagated in a nearly straight line through the entire domain without any noticeable deflection. Nevertheless, the calculated J-integral shows toughening behavior when the crack reaches ZrB 2 -C interfaces. EFT for the sample with 10 vol% of the C-rich phase and an offset initial crack ( Fig. 12(e)) was 38.9 J/m 2 , which is similar to the sample with 10 vol% of the C-rich phase (38.4 J/m 2 ) in Fig. 12(b); the reason is that for both microarchitectures with 10 vol% of the C-rich phase in Fig. 12(b) and (e), the crack path was almost straight, and the toughening came from the increased driving force that was needed for the crack to pass through the vertical interfaces between the hard and soft phases. Some branching of cracks was seen for the case with 30 vol% of the C-rich phase (Fig. 12(c)). EFT for 30 vol% C-rich phase was 40.7 J/m 2 which was 70% higher than G c of the ZrB 2 -rich phase. Even though the crack path in Fig. 12(d) was entirely in the C-rich phase as it was in microstructures with 50 vol% of the C-rich phase for hexagonal and triangular cells, the 90°deflection of the  crack around the ZrB 2 -rich cells affected the crack driving force and resulted in an increased EFT of 33.6 J/m 2 . Having sharp 90°edges in ZrB 2 -rich square cells can result in a higher EFT for the ZrB 2 -C ceramics than having hexagonal or triangular microarchitectures. Fig. 13(a) shows the calculated J-integral for FMs with circular cells having different volume fractions of the C-rich phase ranging from 10 to 50 vol%. In Fig. 13(b), the FM with 10 vol% of C-rich phase, the circular cells are in contact with each other, and the regions between them are filled with C-rich phase. Increasing the driving force for crack propagation is necessary for crack deflection inside circular cells or crack propagation around ZrB 2 -rich circular cells. This raised driving force increased the Jintegral as shown in Fig. 13(a). For the FM with 20 vol% of the C-rich phase (Fig. 13(c)), EFT was 35.5 J/m 2 which is similar to the EFT for the sample with 10 vol% C-rich phase. For the composite with 30 vol% Crich phase (Fig. 13(d)), there is still a 90°crack path rotation around ZrB 2 circular cells and resulted in an EFT of 33.5 J/m 2 . Even the crack path in the composite with 50 vol% C-rich phase (Fig. 13(e)) was not completely straight inside the C-rich phase, the heterogeneous system EFT was about the same as the EFT for pure C-rich phase. The maximum EFT with a circular-cell microstructure was about 61% of the ZrB 2 crack surface energy for C-rich phase fractions less than 20 vol%.

Conclusion
In this paper, we demonstrated the ability of a modified FP model in evaluating the EFT of composite ceramics, and we applied this computational framework to determine the EFT of different engineered microarchitectures of ZrB 2 -C. In this model, nondimensional forms of the equations for a two-phase heterogeneous system are presented. Many other FP models for crack propagation  have utilized the dimensionless simulations, and the length scale used to non-dimensionalize the equations ( L 0 ¼ ffiffiffiffiffiffiffiffiffiffi ffi G c =E p ) can be used only in homogeneous systems, because in heterogeneous systems, different phases have different ffiffiffiffiffiffiffiffiffiffi ffi G c =E p and k = k/L 0 which can cause serious numerical issues.
To validate the model, a cellular ZrB 2 -C FMs with different vol% of C-phase were produced via a thermoplastic forming method. The specimens were tested to measure the work of fracture as an indicator of the fracture property. Then, EFT was numerically calculated by applying a moving mode-I crack opening displacement boundary condition for different percentages of the C-rich phase in ZrB 2 -C FM ceramics. Hexagonal-cell FM with 10 and 30 vol% of the C-rich phase had similar EFT values, and by increasing the C-rich phase content to 50 vol% significantly the EFT dropped significantly; these results were consistent with the experimental results for fracture behavior. To the best of our knowledge, this is the first time that a PF model is used in real dimensions and with the actual material properties to calculate the EFT of ZrB 2 -C engineered microarchitecture ceramics.
After verifying the model, the effect of different microarchitectures on the EFT of ZrB 2 -C UHTC was studied. Triangular-cell ZrB 2 -C FMs ceramics behave similar to FMs with the hexagonal cells. EFT was increased to 42% compared to the G c of pure ZrB 2 by having 10-30 vol% C-rich phase. Increasing the vol% of C-rich phase to 50% eliminated any crack deflection and resulted in a nearly straight crack path, dropping the EFT to the G c of the C-rich phase (~22 J/m 2 ). On the other hand, FM with square cells with 50 vol% of the C-rich phase had a much higher EFT than the G c of the pure C-rich phase (~45% higher), which was a result of crack deflection around the hard-phase. The EFT in circular-cell ZrB 2 -C FM ceramics can be increased up to 70% of G c in pure ZrB 2 by having 30 vol% C-rich phase. Replacing hexagonal ZrB 2rich cells with cylindrical, triangular, or square cells with the same phase fractions showed that the orientation of the interface between the soft and hard phase with respect to the crack growth direction affects the energy required for crack propagation. The relative orientations of crack and the soft-hard interface can suppress crack propagation or promote crack deflection, either of which results in a higher value of EFT.
In this study, different ratios of the same hard and soft phases were considered. Other factors, such as different compositions for the soft and hard phases that change the properties of each phase and more complex distributions of phases should be studied in the future to find potential routes in further improving the EFT of composite ceramics with engineered microstructures. Reliable predictions of the EFT of engineered architectures can be used to guide the design of composites with enhanced damage tolerance.

Date availability
Data and results will be made available upon reasonable requests.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.