Modeling Deformation and Fracture of Boron-Based Ceramics with Nonuniform Grain and Phase Boundaries and Thermal-Residual Stress

: A phase ﬁeld framework of elasticity, inelasticity, and fracture mechanics is invoked to study the behavior of ceramic materials. Mechanisms addressed by phase ﬁeld theory include deformation twinning, dislocation slip, amorphization, and anisotropic cleavage fracture. Failure along grain and phase boundaries is resolved explicitly, where Weibull statistics are used to characterize the surface energies of such boundaries. Residual stress incurred by mismatching coefﬁcients of thermal expansion among phases is included. Polycrystalline materials of interest are the ultra-hard ceramics boron carbide (B 4 C) and boron carbide-titanium diboride (B 4 C-TiB 2 ), the latter a dual-phase composite. Recent advancements in processing technology enable the production of these materials via spark-plasma sintering (SPS) at nearly full theoretical density. Numerical simulations invoking biaxial loading (e.g., pure shear) demonstrate how properties and mechanisms at the scale of the microstructure inﬂuence overall strength and ductility. In agreement with experimental inferences, simulations show that plasticity is more prevalent in the TiB 2 phase of the composite and reduces the tendency for transgranular fracture. The composite demonstrates greater overall strength and ductility than monolithic B 4 C in both simulations and experiments. Toughening of the more brittle B 4 C phase from residual stress, in addition to crack mitigation from the stronger and more ductile TiB 2 phase are deemed advantageous attributes of the composite.


Introduction
Polycrystalline ceramic solids typically demonstrate high hardness, a high elastic modulus, and a low mass density relative to engineering metals. These mechanical properties, as well as various thermal, electrical, and optical properties make ceramics attractive for use in numerous industrial applications. However, the intrinsic susceptibility to fracture limits the ductility of ceramic solids, and therefore, methods are continuously sought to inhibit cracking and thereby improve mechanical performance for applications requiring structural integrity.
One possible means of potentially improving mechanical properties involves mixing two or more crystalline phases of differing chemical compositions to produce a polycrystalline ceramic composite. Examples include diamond-silicon carbide [1][2][3], silicon carbide-titanium diboride [4], boron carbide-titanium diboride [5][6][7], boron carbide-silicon carbide [8], boron carbide-zirconium diboride [9], and alumina-zirconium dioxide [10]. Such ceramic blends feature a host of potentially beneficial microstructure features, such as grain and phase boundaries of various strengths, grains with contrasting mass densities and thermoelastic and plastic properties, and residual stress fields. The latter have been noted to improve fracture toughness in some cases [4,5]. Residual stresses often arise during cooling from higher processing temperatures as a result of differing thermal expansion coefficients of constituent phases. Hence, these stresses are referred to as thermal-residual stresses. Such residual stresses may also emerge in single-phase materials from the thermoelastic anisotropy of non-cubic grains [11], though the magnitudes tend to be smaller and are often ignored in the analysis of multi-phase ceramics [4]. The tendency of intergranular versus transgranular fracture is altered by additional phases, and grains of higher strength or ductility may block sharply propagating cracks.
The present work focuses on phase field modeling of the mechanical performance of polycrystalline boron carbide (B 4 C), with or without a second phase of titanium diboride (TiB 2 ). Numerical simulations explicitly resolve individual grains of each phase, as well as grain and phase boundaries. The finite element (FE) method is used, wherein constitutive models invoke phase field theory for elasticity, fracture, and other structural changes. Titanium diboride has a higher elastic modulus, lower Poisson's ratio, higher fracture energy, lower hardness, and higher mass density than B 4 C. (Numerical values of physical properties, with supporting references, are given later in Section 3 of this work). They are also thought to be more ductile as evidenced by dislocation nucleation and motion [12,13], though both ceramics are brittle relative to typical metals at room temperature. The particular material of study is produced by spark-plasma sintering (SPS), where fabrication methods and material properties culled from physical experiments are discussed by Rubink et al. [6]. Polycrystalline ceramics synthesized in this way are nearly fully dense (≥99% theoretical density), with little to no graphitic inclusions and relatively uniform and small grains (e.g., mean sizes on the order of several microns). X-ray diffraction (XRD) and theoretical mechanics analyses show that residual stresses are tensile for TiB 2 grains and compressive for more brittle B 4 C grains [6]. Such residual stresses should prove beneficial for impeding cleavage fractures in the latter phase.
The particular phase field model used in the current research, which includes residual stresses from possibly different physical sources, was developed in several recent works [14,15]. State-of-the-art work implementing phase field fracture models accounts for elastic and fracture anisotropy [16][17][18][19], finite strains [19,20], and dynamic effects from material inertia [21]. Research closely related to the present study incorporates fracture or deformation twinning in a dissipative phase field framework applied to B 4 C single crystals [22]. Over the past decade within the mechanics research community, the phase field method seems to have overtaken in popularity the cohesive zone FE method of modeling fracture in crystalline solids [23][24][25]. The origins of the current variational phase field framework for deformation twinning and fracture are found in prior works: [26,27], respectively. Other recent applications to ceramics focus on the homogenized shock compression response of polycrystalline B 4 C and TiB 2 [14,28]. Generalizations incorporating geometric concepts are reviewed in a recent article [29]. The effects of high confining pressure, which can drastically increase resistance to fracture in ceramics, are not addressed in the current investigation. Advancements in three-dimensional (3D) computational methods at the scale of the microstructure are required to account for such effects, which have been modeled only in a one-dimensional (1D) context for uniaxial strain compression [14,28].
A recent work [15] interrogated the mechanical response of B 4 C and B 4 C-TiB 2 under primarily tensile stress states, which are directly relevant to the flexure experiments reported for these materials [6]. In the context of indentation experiments, tensile principal stress components may arise in biaxial stress fields near the surface outside the indentation region and along the contact perimeter, leading to radial cracking and/or ring cracking depending on the geometry, whereas lateral cracks may arise in the subsurface [30,31]. Surface crack length is used to quantify fracture toughness in the Vickers indentation method, for example [32]. Most simulations by Clayton et al. [15] invoked uniform fracture strengths (i.e., surface energies) for all grain and phase boundary facets, though a few exploratory simulations considered completely random strength distributions [33]. The results showed, for purely tensile loading, that fractures tended to initiate earlier in materials with random versus uniform distributions, leading to a reduced peak strength on average. Other trends regarding the benefits of plasticity and residual stresses, which agreed with experimental inferences [6,34], persisted when either type of distribution was modeled.
The objective of the current research is understanding the effects of composition, microstructure, residual stress, inelastic deformation mechanisms, and imposed biaxial strain magnitude on B 4 C and B 4 C-TiB 2 . The current application of the variational phase field theory with residual stress [14,15] includes the following innovations over previous studies: • Weibull distributions of fracture strengths of grain and phase boundaries, which are known to be more physically realistic for polycrystalline ceramics [35], including B 4 C [36], than uniform or purely random distributions; • Biaxial strain states more pertinent than uniaxial tensile states to surface cracking [37] and also more relevant for failure by shear-driven median cracking in brittle elastic materials for regions directly under the indenter [30,38]; • Revisited comparison of cumulative current and past simulation results with detailed experimental findings [6,34].
The general phase field theory including fracture, inelasticity, and residual stresses is described in Section 2. Material aspects pertinent to the present ceramic systems, including grain and phase boundary descriptions and rendering of microstructures, are detailed in Section 3. Numerical techniques including the imposition of boundary conditions and averaging of field variables for subsequent analysis are explained in Section 4. Model results, including a comparison with past experimental findings, are given in Section 5. Conclusions appear in Section 6.

Phase Field Mechanics
The general phase field theory of Clayton [14] accounts for finite strains, temperature changes, rate effects, and initial stresses of arbitrary origin. This theory was recently extended [15] with a detailed, thermodynamically consistent treatment of residual stresses, e.g., arising from thermal processing history. Here, a geometrically linearized theory is implemented, which is sufficient for the loading protocols of the current study. The applied strain rates are low enough that a quasi-static, variational approach suffices, wherein viscous dissipation is negligible. The theory of Clayton et al. [14,15] is presented next in Section 2 in abbreviated form, specialized thereafter in Sections 3-5 to low-rate loading of B 4 C and B 4 C-TiB 2 .

Coordinates and Order Parameters
Denote by X, with covariant components X K , the position of a material particle. Rectangular Cartesian coordinates are used exclusively, with {E K } the fixed coordinate basis and E K = ∂X/∂X K = ∂ K X a basis vector (K = 1, 2, 3 in ambient Euclidean threespace). A solid body occupies a reference domain Ω 0 with boundary ∂Ω 0 . The outward unit normal to ∂Ω 0 is N(X). For a quasi-static theory, the explicit time (t) dependence of fields is omitted.
Two order parameter fields are η(X) and ξ(X). Field η measures structural changes such as slip and deformation twinning, depending on the crystal type. Field ξ measures fracture, i.e., damage in the solid involving the formation of free surfaces. In the phase field method, slip bands, twin boundaries, and fracture surfaces are diffuse or regularized, rather than resolved discretely. The values of order parameters η ∈ [0, 1], ξ ∈ [0, 1] physically represent the following: (1)
From thermodynamic arguments requiring non-negative dissipation [14,39] with rate independence restrictions, Governing equations can alternatively be derived with variational methods [26,27]. Assume Ψ is of the functional form in (10) and that the following variational principle holds: with V a region enclosed by surface S upon which natural boundary conditions are applied. Assuming field variables are sufficiently differentiable, (13) produces [26,40] ∇ · ∂Ψ ∂h δΨ δη Equations (14)- (16) are consistent with (7)-(9) and (12). These hold for quasi-statics, isothermal states, and inviscid kinetics. The variational derivative with respect to generic argument x is δ(·)/δx [41]. The dependence of Ψ on h is now replaced with E of (4): A change of variables from h to E with (6) leads to Recalling the present theory is geometrically linear, Cauchy stress obeys, from (12), (17), and (18), The internal microforces in (12) become, with (19), (20), and the symmetry of P, where thermoelastic driving forces are the rightmost terms after each equality involving the structure tensors of (6). Denote by W Ψ the strain energy, f Ψ η a phase energy with prominent dependence on η and its gradient, and f Ψ ξ on ξ and its gradient. Then,

Twinning, Plastic Shear, and Dilatation
Distortion h η arises from twinning shear or localized slip, described here on a single plane with unit normal m. The shearing direction is the unit vector s, and s · m = 0. The saturation magnitude of shear is γ 0 ; for deformation twinning, this is often known from the crystal structure [42]. Isochoric inelastic distortion, inelastic strain, and the second structure tensor of (6) are, respectively, The interpolation function φ is defined, for argument x, as follows, with the corresponding derivative: In ductile solids, residual dilatation may arise from nucleated voids. In brittle solids, dilatation may be caused by shear-induced bulking, especially in high-pressure regimes [43] or granular flow regimes upon impact and dynamic comminution [28,44]. In the present applications, the average pressure and strain rates are small, fragmentation is minimal, and bulking is assumed negligible. Thus,

Isotropic Elastic Formulation
In a heterogeneous polycrystal with grains of different crystal types or phases, each crystal is assigned to a different sub-body Ω 0 with potentially different physical properties. Elasticity is isotropic, but fracture is not whenα > 0, whereα is defined briefly. Slip and twinning are, by the kinematic definition (26), always anisotropic.
Elastic strain energy includes isotropic elasticity and isotropic residual stress of thermodynamic origin [15]: The deviatoric thermoelastic strain is E = E − 1 3 (tr E )1, where 1 is the unit tensor: (1) I J = δ I J . The tangent bulk modulus is B, and the tangent shear modulus is G. A nonlinear pressure-volume response is enabled via β 0 ≈ 2B 0 − 1 [14], where B 0 is the ambient pressure derivative of the bulk modulus. Thermoelastic coefficients are usually constants in a homogeneous single crystal, but generally can depend on the order parameters: Superscripts (·) 0|0 and (·) 0|1 correspond to undamaged states before and after structural transformation, with (ξ = 0, η = 0) and (ξ = 0, η = 1), respectively. The bulk modulus is preserved under compressive pressure, but the shear modulus is not [20,27,45]. The indicator function for volumetric elastic compression versus extension (corresponding to positive versus negative pressure) isω; the interpolator φ is defined in (27); the degradation function is χ(ξ) ∈ [ζ 0 , 1], with 0 < ζ 0 1 [40]: The right-continuous unit Heaviside function is {H : For each anisotropic domain, cleavage fracture may localize on a critical crystallographic plane with unit normal vector M(η(X)). Denote the cleavage anisotropy factor [45] byα ≥ 0, where for isotropic f Ψ ξ ,α = 0. The fracture regularization length is l ξ > 0; the fracture surface energy is Υ. Then, Energy functions for η-equilibrium contain surface energy Γ, bulk energyÂ, and phase field regularization length l η . The damage-transformation coupling function [40] isχ(ξ), typically assumed asχ = χ orχ = 1 (i.e., full or no coupling). Assigning the first of (24), the first of (9) follows: Total stress P is computed from (29), where the residual contributionσ 0 is spherical for isotropy: Elastic strain E is measured with respect to a datum state ( E = 0) at which P =σ 0 1. In a cubic or an isotropic (poly)crystalline solid, lattice spacing at the datum state could be determined from the inversion of a pressure-volume equation of state at ambient temperature and datum pressureσ 0 . A few remarks regarding terms entering the free energy function of (23) are in order. The surface energy of cracks is contained in the term f Ψ ξ , which follows the usual sum of quadratic forms in ξ and ∇ξ of phase field fracture models [46], further extended for damage anisotropy [17,45]. The surface energy of twin boundaries and stacking faults, as well as the possible energy of phase transformations and dislocation lines, is contained in the term f Ψ η , which again follows the usual sum of quadratic forms [41]. The energy f Ψ η also degrades with local damage through functionχ(ξ) [40], regardless of the local stress/strain state. Elastic strain energy density W Ψ consists of three terms in (29): the volumetric strain energy scaled by the tangent bulk modulus B, the deviatoric strain energy scaled by the tangent shear modulus G, and the residual strain energy scaled by the factor Λ defined later in Section 2.6. According to (30) and (31), the tangent bulk modulus is interpolated from the bulk moduli of coexisting phases, where the interpolated value is degraded by local damage ξ under tensile elastic strain states (tensile pressure), but not for compressive elastic strain states (compressive pressure) [27,47]. The tangent shear modulus is interpolated from the shear moduli of coexisting phases, degraded by local damage ξ regardless of local stress/strain state. The residual strain energy attributed to thermal-residual stress, which is hydrostatic for isotropic materials, is degraded similarly, but not identically, to the bulk strain energy according to (39) and (40), i.e., a reduction for tensile states, but not compressive states. Accordingly, fracture will not occur if the loading is of purely hydrostatic compression, while Mode I cracks can nucleate and propagate for tensile states, regardless of shear. For combined compressive pressure and shear stress, the fracture tendency is reduced, but not eliminated, since the reduction of the shear modulus acts as a driving force for Mode II-/III-type cracks. Under such compressive-shear states, a brittle material can become pulverized (e.g., severely reduced tangent shear modulus), yet fluid-like, since the bulk modulus is not degraded. Further details are available in related work [14,15], noting that alternative schemes for tensile versus compressive stress states have also been implemented, e.g., [46].

Initial Stress
The model of Clayton et al. [15] is invoked for residual stresses in the ambient state. Define function y = tr E / C , with C > 0 a constant. The integral of the first of (27) is, with y the argument, Function Λ(tr E , ξ) controls the degradation of initial stress due to increasing ξ, for tensile states: Then the total isotropic stress contribution fromσ 0 and its contribution to the thermodynamic driving force for fracture are, respectively, The critical tensile pressure P C and critical volume strain C are defined as [15] The initial elastic modulus is E 0 = (2 + 2ν 0 )G 0 , and Poisson's ratio is ν 0 . In (42), these elastic constants are all evaluated at the initial (parent) state with η = 0 if nominally phase-dependent.

Material Properties and Polycrystalline Microstructures
The materials of study are B 4 C-23 vol. % TiB 2 and, for comparison, polycrystalline B 4 C. Differences from [15] implemented here are the Weibull statistics for grain boundary fracture properties [35] and a focus on the effects of different physical mechanisms (i.e., residual stresses, plasticity) under biaxial strain boundary conditions. Fracture and plasticity incurred deep in the material directly under indentation are affected by the confining pressure of the surrounding material. In the stress field for Vickers indentation of a homogeneous isotropic elastic body, all principal stresses are compressive, so fractures are shear-induced [37,38]. However, biaxial residual stress fields with a tensile component (typically smaller than the compressive component) arise in elastic-plastic or fractured materials [38]. Pure shear (biaxial) loading is thus thought more qualitatively representative than uniaxial tensile loading for stress fields incurred during indentation, though surface stress distributions for pyramidal indenters are non-uniform and not equi-biaxial at most locations.
Tensile loading conditions, with and without superposed simple shear, were thoroughly investigated [15] with the intent of comparison to flexure experiments. Only a few demonstrative simulations with biaxial and simple shear were executed in prior work [15], without the isolation of the individual effects of the full range of physical mechanisms available in the phase field model.

Boron Carbide Phase
The properties and parameters for bulk grains of B 4 C are listed in Table 1. Isotropic elasticity is assumed [34,48], and nonzero β 0 is invoked in the compressive regime. To prevent tensile numerical instabilities [15], 1 3 (29), where (·) = 1 2 [(·) + |(·)|]. Twinning occurs on 1010 {0001} and fracture on basal {0001} planes, giving M. Amorphous shear bands are enabled to form within twinned domains [49][50][51] as η → 1. Twinning shear is γ 0 = 2 3 a/c, with a and c lattice parameters. The surface energy resisting twinning is Γ. The energy barrierÂ is obtained from the difference between the free energies of crystalline and amorphous phases. The preferred twinning system in a randomly oriented B 4 C grain is the one of among three possible basal twinning systems with the maximum Schmid factor, producing s and m [52]. Cleavage anisotropyα is prescribed similarly based on prior works [27,45,52], andχ = χ in (35) [40].
The choices of the regularization lengths l ξ and l η were analyzed in detail in a prior work on ceramic composites [52]; that discussion is now summarized. Here, the minimum acceptable length permitted by the mesh resolution and dimensions of the simulation domain [45,48] is assigned, with l ξ = l η for convenience. If explicitly resolving atomicscale phenomena, regularization lengths can correlate physically with perturbed atomic displacements in the vicinity of crack faces [53] or twin boundaries [26,54]. In such nanoscale simulations, parameters correlate with molecular mechanics and take values of nm or smaller. Parametric studies with different regularization lengths for twinning, phase transformations, and fractures are reported in several prior works [26,45,52,55]. Generally, increasing the regularization length while holding the surface energy constant leads to a reduction in nucleation stress and peak strength, though the relative importance depends on the particular problem. In one prior study [52], varying the ratio l η /l ξ from 1 2 to 2 produced similar mechanisms in the simulation results, but the predicted widths of regularized cracks and twin boundaries increased and the homogenized peak stresses decreased, with increasing regularization lengths.
Partial twinning is reversible upon unloading, as observed in other investigations [56,57]. The irreversibility of amorphization is imposed for η ≥ η T , with η T = 0.9, and likewise for fracture, with ξ T = 0.9. Constraints enforcing irreversibility in numerical simulations are described in prior work [40,51]. Different choices of ξ T were further investigated in the preparation of prior works on phase field fracture [27] and fracture with twinning [40]. Too small a value of the threshold for the irreversibility of fracture locks in permanent, diffuse damage contours, especially under relatively homogeneous loading protocols, wherein sharp cracks do not nucleate immediately. On the other hand, omitting the threshold for fracture irreversibility entirely leads to healing of fully developed cracks upon elastic unloading as the system seeks to minimize total elastic and surface energy. Similar remarks apply for amorphized, twinned, or slipped regions. The aforementioned prior investigations found that, for the present constitutive theory and variational numerical implementation, the presently prescribed value of 0.9 for ξ T and η T provides a suitable compromise, producing physically realistic results for relatively brittle materials of study.

Titanium Diboride Phase
The properties and parameters for grains of TiB 2 are likewise listed in Table 1. Cleavage occurs on {0001}, giving M, where Υ is obtained from DFT [34]. Slip occurs on families of systems with the lowest stacking fault (SF) energy barrier Γ from DFT [34]: basal slip 1120 {0001} with (full) Burgers vector b 0 = a. In the phase field model for restricted slip in TiB 2 [34,48], η is related to the magnitude b of the Burgers vector of a partial dislocation and cumulative dislocation density ρ D . It is is used to interpolate from null to saturation shear γ 0 on a single system in (26). Regularization length, cleavage anisotropy, and irreversibility are addressed similarly to B 4 C; values are justified in prior work [15]. The active slip system was chosen among six possible basal systems (three pairs with parallel s of opposing signs) by the maximum applied Schmid factor [52], giving s and m. The maximum distance a typical dislocation travels prior to arrest isx = 50 nm [12]. Denote by ρ S = 10 15 m −2 the saturation density of dislocations. The stored elastic energy density of dislocation lines is [42]. Then [34,48],

Grain and Phase Boundaries
As in some prior work [3,15,52], larger crystals of each phase are encapsulated in an isotropic binder with effective elastic and fracture properties. This binding matrix accounts for the effects of phase and grain boundaries and very small grains that cannot be resolved discretely in FE simulations. Fracture properties can vary among boundary regions (e.g., different grain and phase boundary facets), each defined by initial position vectorX. The matrix has the same volume fractions of constituents as the entire composite. Thus, isotropic elastic properties of the composite [14] are used for the matrix. Superscripts (·) (0) , (·) (1) , (·) (2) correspond to (B 4 C, TiB 2 , boundary) phases. Let (υ (0) , υ (1) ) = (0.77, 0.23) be the initial volume fractions of the (B 4 C, TiB 2 ) phases of the composite, where υ (0) = 1 − υ (1) . In the absence of experimental or atomic simulation data, the fracture surface energy for boundary regions follows a mixture rule with possible heterogeneity, the latter scaled according to global Weibull statistics via location-dependent functionr: When homogeneous boundary energies are imposed, Υ (2) = Υ 0 is independent ofX. When heterogeneous energies are imposed, Υ (2) contains dimensionless multiplierr, where a different instantiation ofr is applied at each boundary regionX. Previous works [15,33] considered completely random (uniform) fracture energy distributions on the approximate intervalr ∈ (0, 2). Weibull distributions on the fracture toughness and strengths of grain boundaries were investigated in cohesive FE simulations of polycrystalline alumina [35] with Weibull modulus m = 10, roughly characteristic of macroscopic strength statistics [36]. The macroscopic Weibull modulus has not been measured for B 4 C-TiB 2 , so m = 10 was used for the present phase field simulations, with the shape factor given by tensile failure stress P 0 ∝ √ Υ 0 [14,66]. The Weibull probability density function f for failure stress P ≥ 0 is then In a given simulation, fieldf was initialized over all grain boundary facets such that P follows (45) in discretized form, squared in (44) to account for Υ ∝ P 2 in the phase field approach (f = f ). The median value of P from (45) is P 0 · (ln 2) 1/m ≈ 0.964P 0 , so the median energy is 0.929Υ 0 , which is slightly lower than that of the homogeneous condition (m → ∞). The anisotropic mechanisms of twinning and dislocation slip were omitted in isotropic boundary regions, andα = 0 accordingly. For simulations of pure B 4 C, a distinction was neither necessary nor imposed between the properties of bulk crystals and grain boundary regions. Since pure B 4 C tends to fracture transgranularly, rather than intergranularly [6], the resolution of the distributions of its boundary energies was deemed uncritical here, though the distributions could be considered given sufficient microscopic data. Macroscopic failure statistics emerged regardless, due to grain-to-grain variations of cleavage plane orientations.

Residual Stresses
The B 4 C-TiB 2 composite is processed at temperature T P T 0 , where T 0 is room temperature. Cooling from T P leads to unequal contraction among B 4 C and TiB 2 grains due to different coefficients of thermal expansion A 0 of each phase. Values ofσ 0 were calculated for each bulk phase using the analytical solution of Taya et al. [4]. In the simulations, each crystal was assigned a uniform value ofσ 0 during a seeding step. Initial residual stresses were redistributed as an outcome of the solution, to maintain stress and order parameter equilibrium, before (see Figure 1b) and after displacement boundary conditions were applied. Boundary regions haveσ 0 = 0 since they are nominally self-equilibrated; however, after equilibration, boundary regions can be initially stressed depending on the morphology and composition of the surrounding microstructure. As in the original deriva- Residual stresses are self-equilibrated: υ (0)σ (0) 0 + υ (1)σ (1) 0 = 0. The fracture toughness of a given phase with compressive residual stress was increased by K R [4]: The mean grain diameter of the less-abundant phase is d. Fracture energy is updated via [15] Υ (i) = (K (i) with K C the fracture toughness corresponding to Υ 0 via the usual linear elastic fracture mechanics relation [14,15]. The outcomes of this model wereσ 0 < 0 for B 4 C (residual compression) andσ 0 > 0 for TiB 2 (residual tension). Fracture energy Υ (2) increased by K (48) inserted into (44), since negative (i.e., compressive) residual stress toughens subscale B 4 C constituents. The corresponding K (2) R was extracted from (48) applied to the binder. The values in Table 1 were obtained from α (0) = 4.8 × 10 −6 /K [58], α (1) = 8 × 10 −6 /K [59], d = 3 µm, and T P − T 0 = 1873 K [6]. Parameter C was found from (42).

Numerical Methods
The FE method with an implicit formulation is invoked for solution of boundary value problems of deforming polycrystals that implement the variational phase field theory and constitutive models and parameters of the respective Sections 2 and 3.

Geometric Rendering and Investigated Parameters
Each polycrystal contains on the order of 50 polyhedral bulk grains. These belong to the B 4 C or TiB 2 phase with the associated parameters of Table 1. As discussed in Section 3.3, a very thin layer of binder material, just wide enough to resolve cracks in the regularized phase field approach, covers each bulk grain. This layer is partitioned into many domains. Each domain corresponds to a boundary layer facet between neighboring crystals, and each domain has a centroidal coordinate vectorX. Heterogeneous fracture properties were implemented from (44).
A representative FE mesh of a B 4 C-TiB 2 polycrystal is shown in Figure 1a. The volume of the domain, Ω 0 , is 10 3 µm 3 . The three phases (B 4 C, TiB 2 , and boundary) are differentiated by color in Figure 1a. Volume fractions of (B 4 C, TiB 2 ), denoted (υ (0) , υ (1) ) = (0.77, 0.23), were exactly reproduced. The average grain diameter was 3 µm; too few grains were included to attempt matching experimentally determined size and orientation distributions. Random lattice orientations were assigned. Some simulations modeled pure B 4 C rather than the composite, in which case, all domains were logically assigned properties of B 4 C. Meshes consisted of approximately 1.5 M hexahedral elements with full integration. As discussed in prior 3D phase field modeling works on boron carbide-and silicon carbide-based ceramic composites with comparable mesh densities [15,52,67], FE meshes are sufficiently refined to resolve regularization lengths for fracture, slip, and twinning and, thereby, mitigate the potential mesh size dependence of solutions.
The simulation results of Section 5 investigated the effects of the following for biaxial loading conditions corresponding to pure shear, rather than uniaxial extension prioritized in a prior study [15]:

Boundary Conditions and Homogenization
Conjugate gradient energy minimization was applied in implicit FE simulations at each imposed displacement increment δ i . Fields u(X), η(X), ξ(X) were sought that minimized Ψ subject to boundary constraints, in accordance with (13). A load parameter¯ is defined as where the edge length of domain Ω 0 is L 0 = 10 µm;¯ quantifies the average strain. Displacement boundary conditions imposed on ∂Ω 0 correspond to pure shear (i.e., biaxial tensioncompression). For example, the macroscopic strain tensor¯ for tension-compression in the Principal directions E J and E K were permutated among J, K = 1, 2, 3, where J = K. (In contour figure legends, The aggregate is traction-free in the direction orthogonal to the plane of shear. An equilibrated stress distribution at¯ = 0 for a residually stressed B 4 C-TiB 2 aggregate is shown in Figure 1b. Note that B 4 C grains were initially compressed, TiB 2 initially extended, and boundaries may be in tension or compression depending on location. Free natural boundary conditions t ξ = 0 and t η = 0 [40,66] are assigned in (16) over all faces of the cube comprising ∂Ω 0 upon which tensile displacement boundary conditions are not applied. On faces of ∂Ω 0 where essential tensile increments δ i are enforced, conditions ξ = 0 and η = 0 are invoked to eliminate premature fracture or yielding of the material on such boundaries [15] (i.e., failure at specimen grips is inhibited). The average stress, fracture parameter, and twinned or slipped fraction areP Local pressure is p(X) = − 1 3 trP(X), and the local von Mises stress is Σ(X) = { 1 2 [(P 11 − P 22 ) 2 + (P 22 − P 33 ) 2 + (P 33 − P 11 ) 2 + 6(P 2 12 + P 2 23 + P 2 31 )]} 1/2 ≥ 0, where P I J = P I J (X). When found from (51), p = − 1 3 trP,Σ = 1 2 [(P 11 −P 22 ) 2 + (P 22 −P 33 ) 2 + (P 33 −P 11 ) 2 + 6(P 2 12 +P 2 23 +P 2 31 )] 1/2 (52) are the average pressure and average effective shear stress.

Model Results
Eighteen simulations (Sims), all differing from those in prior work [15], were undertaken with the pure shear boundary conditions of (50). The features are listed in Table 2. Polycrystals were assigned one of two material categories, either B 4 C-23 vol. % TiB 2 (Sims 1-12) or all B 4 C (Sims 13-18), one of three lattice orientation distributions, and one of three permutations of loading directions J, K corresponding to applied strain¯ JK . Only three out of a maximum of six permutations, when sign reversals were considered, were simulated. Twinning and slip were disabled in Sims 4-6 and 16-18. Residual stresses were omitted in Sims 7-9 and 13-18, where the latter six contain all B 4 C. Weibull distributions of grain and phase boundary fracture stress, (44) withr differing from unity, were used in most representations of the composite (Sims 1-9); otherwise, Υ (2) is spatially constant. Shown in Figure 2 are the results for Sim 2, which corresponds to the B 4 C-23 vol. % TiB 2 composite with slip and twinning enabled, residual stresses enabled, and a Weibull distribution of effective grain and phase boundary strengths. Stress normal to the extension direction is visualized in Figure 2a, where local relaxation is evident in many fractured regions. The fracture parameter contour for ξ of Figure 2b demonstrates a preponderance of fractures along grain and phase boundaries, with very sparse fractures within large grains of either phase. Inelasticity order parameter η in Figure 2c shows the largest values in TiB 2 grains (0.5 η ≤ 1.0) preferentially oriented for basal plane slip. Low values in B 4 C crystals (0.05 η 0.25) correspond to diffuse twinning or stacking fault formation prior to shear localization. However, a few localized zones where η approaches unity do arise in B 4 C, which correspond to amorphous shear band formation. (c) inelasticity Fracture contours for six simulations (Sims 1, 4, 7, 10, 14, 17), i.e., one of each grouping with similar characteristics in Table 2, are presented in Figure 3. Comparing Figure 3a with Figure 3b, transgranular fractures, particularly in the TiB 2 phase, are more prevalent when η representing basal slip is suppressed. Similarly, in Figure 3c, the tendency for transgranular fracture is increased in the absence of residual stress. When uniform (i.e., constant) grain and phase boundary strengths were imposed rather than Weibull statistics, fewer distinct cracks emerged and at different locations, as is evident in Figure 3d.
With or without twinning and shear localization (i.e., amorphous banding) enabled, the pure B 4 C material is dominated by transgranular fractures, along the same plane of maximum shear stress in Figure 3e,f. The fractured zone is physically wider in the latter case at the same¯ , but η does not appreciably alter the crack morphology. It was found, however, that the shear strain at peak load was reduced by an average of 0.08% when inelastic shear was incorporated in the simulations of pure B 4 C. Thus, ductility is compromised by the transition from twinning to amorphous banding, leading to earlier catastrophic fracture in the monolithic material.   The effective average stressΣ of (52) is reported versus the applied strain¯ of (50) for 15 simulations in Figure 4, where results from the full composite model (Sims 1, 2, 3) are contrasted with those with other model features enabled/disabled in Sims 4 through 15 in each subfigure. For every simulation, the maximum value ofΣ attained is listed in Table 2, along with the following average field variables at the applied strain corresponding to peak stress: average pressurep, average inelasticity order parameterη, and average fracture order parameterξ. The following trends are noted from Figure 4,   When slip and twinning mechanisms are potentially active, peakΣ appears higher, as noted in related work under different boundary conditions [15]. However, further analysis revealed that the applicability of this trend is sensitive to modeling assumptions and parameters for the phase field representation of inelasticity, as well as the microstructure (e.g., lattice orientation) and other details entering the stress state homogenization procedures. For example, peak strength was found to decrease, but ductility was found to increase, when twinning was more profuse in a diamond-SiC ceramic composite [52]. In another study [34], the tendency for strength increase or decrease was found to depend on lattice orientation affecting the activity, or lack thereof, of inelastic deformation systems. Above a small threshold, or with the availability of more systems, plasticity may reduce hardness and peak load capacity as the tendency to flow inelastically increases [52]. As noted in the present findings and experiments [68], a transition from inelastic shear (e.g., twinning) to amorphization may inhibit ductility in B 4 C. Regardless, a consistent trend among current and prior work is a decrease in average fracture activityξ with increasing average slip or deformation twinningη. Along with local crack tip blunting, this should contribute to a higher toughness, even if ultimate strength is reduced.

Comparison with Experiments and Prior Modeling
Mechanical experiments reported by Rubink et al. [6] include static and dynamic flexure via three-point bending, Chevron notched bending for fracture toughness, and Vickers indentation for hardness and fracture toughness. Scanning electron microscopy (SEM) was used to investigate post-mortem fracture paths in samples recovered from these tests. Transmission electron microscopy (TEM) was used to further investigate fracture and dislocation plasticity mechanisms at finer resolution, where a focused ion beam (FIB) was used to extract small-scale samples in the vicinity of cracks induced by Vickers indentation. Synchrotron X-ray diffraction (XRD) was used to measure lattice spacing and infer residual stress states of grains of each phase.
Analysis of SEM images showed crack deflection and crack bridging by second-phase TiB 2 particles. The high surface roughness of the cracked material suggested intergranular fracture in the composite. In contrast, pure B 4 C failed predominantly by transgranular fracture, similar to observations in other studies [68]. The same physical mechanisms were observed in both static and dynamic flexure. The analysis of TEM images corroborated the fracture tendencies from SEM, i.e., more crack branching in the composite than the pure B 4 C. The TEM images further showed stacking faults and dislocations in the TiB 2 grains, as observed in earlier work [12,13].
The analysis of XRD peaks in the composite showed definitive compressive lattice strains for B 4 C and milder tensile lattice strains for TiB 2 . The addition of TiB 2 to B 4 C yielded higher flexure strength, a higher elastic modulus, and higher fracture toughness [6]. The lower toughness of pure B 4 C was quantified by longer cracks produced by Vickers indentation; however, the higher hardness of pure B 4 C was quantified by smaller residual indents. Similar improvements in the mechanical properties were measured in earlier work on a previous generation of SPS ceramics [34]. Improvements in tensile strength and toughness were likewise measured in hot-pressed B 4 C-TiB 2 composites of lower TiB 2 content (10%) and with Si/B doping [68,69].
The following differences in the experimental findings are summarized for the dualphase composite with 23% by volume of the second phase, i.e., B 4 C-23 vol. % TiB 2 versus pure B 4 C, both produced by SPS [ The first six items listed above are beneficial for the composite, the last two detrimental. The apparently beneficial properties of dislocation plasticity for mitigating fracture in the TiB 2 phase have been noted for other ceramics such as aluminum nitride [70]. Twins and amorphous shear bands have been observed in the B 4 C phase. Twins have been suggested to improve the mechanical properties of B 4 C [71,72], while solid-state amorphization is thought to degrade strength, since it tends to be accompanied by fracture and fragmentation [73]. Doping with silicon (Si) and boron (B) has been shown to reduce the tendency for amorphization in B 4 C [68,69]; the addition of TiB 2 to such doped materials produces increased flexure strength and toughness through the mechanisms listed above [6], as well as grain refinement.
The present model results corroborate the above experimental findings. Under biaxial loading, when the Weibull distributions of boundary strengths and thermal-residual stresses (where compressive stresses intrinsically toughen the more brittle B 4 C phase) are implemented, the phase field predictions of improved peak strength, reduced fracture tendency (i.e., increased overall toughness), and increased plasticity and ductility with the addition of TiB 2 all concur with experimental data [6,34]. For supporting quantitative evidence, see Figure 5, which compares numerical averages, over three simulation cases, for the composite with and without residual stress and pure B 4 C. Average strength, ductility, and plastic deformation are all higher in the composite with residual stress than in the composite when residual stress is omitted and are also all higher in the composite relative to pure B 4 C. In both the model and experiment, intergranular fracture arises more often than transgranular fracture in the composite, as opposed to strongly preferential cleavage fracture in the monolithic B 4 C, which is clear in the model results upon comparison of Figure 3a with Figure 3e. Prior recent work [15] for tensile loading conditions and uniform boundary strengths produced similar trends; furthermore, the toughness, albeit at a lower scale of resolution than the Vickers experiments, was quantified for pre-cracked samples, providing a respectable validation of the model and parameters. The reduction in peak strength with increased variability of the properties, attributed to more weak links for fracture in a probabilistic sense, is shared among tensile and shear loading conditions and among random and Weibull boundary strength distributions.
Overall, including TiB 2 is recommended to improve strength and ductility. However, hardness is reduced and mass density is increased with the addition of TiB 2 . Thus, applications requiring simultaneously high hardness, strength, and ductility, as well as low weight will require moderation of the ideal fraction of TiB 2 , i.e., a trade-off among physical properties. Furthermore, a maximum improvement in strength and toughness exists at a fraction of TiB 2 under 100% since pure TiB 2 tends to demonstrate lower flexure strength and lower fracture toughness than the B 4 C-23 vol. %TiB 2 composite [34].   (c) plasticity measured by averaged slip/twinning order parameterη at peak stress.

Conclusions
A phase field model of fracture and inelasticity accounting for thermal-residual stresses was implemented in numerical simulations of boron-carbide-based ceramics and ceramic composites. The new computations included the effects of the Weibull distributions of grain and phase boundary energies for biaxial loading of polycrystalline B 4 C and B 4 C-TiB 2 aggregates. The results confirm the experimental observations and trends observed in prior numerical work that focused on tensile loading and uniform boundary strengths. Improved mechanical properties manifest from the toughening effects of residual stresses, plasticity, and cleavage crack blockage by stronger TiB 2 grains and more tortuous intergranular crack paths exhibited by the composite relative to the dominant transgranular failure exhibited by the pure B 4 C.