Computational Modeling of Dual-Phase Ceramics with Finsler-Geometric Phase Field Mechanics

A theory invoking concepts from differential geometry of generalized Finsler space in conjunction with diffuse interface modeling is described and implemented in finite element (FE) simulations of dual-phase polycrystalline ceramic microstructures. Order parameters accounting for fracture and other structural transformations, notably partial dislocation slip, twinning, or phase changes, are dimensionless entries of an internal state vector of generalized pseudo-Finsler space. Ceramics investigated in computations are a boron carbide-titanium diboride (B4C-TiB2) composite and a diamond-silicon carbide (C-SiC) composite. Deformation mechanisms-in addition to elasticity and cleavage fracture in grains of any phase-include restricted dislocation glide (TiB2 phase), deformation twinning (B4C and β-SiC phases), and stress-induced amorphization (B4C phase). The metric tensor of generalized Finsler space is scaled conformally according to dilatation induced by cavitation or other fracture modes and densification induced by phase changes. Simulations of pure shear consider various morphologies and lattice orientations. Effects of microstructure on overall strength of each composite are reported. In B4C-TiB2, minor improvements in shear strength and ductility are observed with an increase in the second phase from 10 to 18% by volume, suggesting that residual stresses or larger-scale crack inhibition may be responsible for toughness gains reported experimentally. In diamond-SiC, a composite consisting of diamond crystals encapsulated in a nano-crystalline SiC matrix shows improved strength and ductility relative to a two-phase composite with isolated bulk SiC grains.


Introduction
Polycrystalline ceramics are solid materials that generally are mechanically stiff yet prone to brittle fracture. Heterogeneous composite ceramics consisting of grains of different chemical compositions or phases have been engineered to provide increased strength, ductility, and toughness relative to monolithic ceramics.
The current work focuses on two such ceramic composites, each consisting of two distinct phases. One such dual-phase material contains mostly crystalline boron carbide (B 4 C, trigonal structure) interspersed with secondary grains of crystalline titanium diboride (TiB 2 , hexagonal structure), wherein grain sizes of both phases are comparable. Versions of this material system are discussed in sigl et al. [Sigl and Kleebe (1995); White and Dickey (2011)]. Single-phase boron carbide is prone to brittle failure by cleavage, and at high stresses, may undergo a phase change to an amorphous solid as well as deform by deformation twinning [Chen, McCauley and Hemker (2003); Clayton (2012); An and Goddard (2015)]. Titanium diboride may impart higher fracture resistance by cleavage crack deflection or bridging, and some ductililty may be afforded by dislocation motion that has been suggested from shock experiments [Vanderwalker and Croft (1988)]. Compressive residual stresses induced in B 4 C from TiB 2 by thermal processing may also improve performance.
The second dual-phase material studied is comprised of diamond crystals (C, cubic structure) either interspersed with or encapsulated by silicon carbide micro-or nano-crystals (β-SiC, cubic structure). Though the hardest of natural materials, diamond is prone to cleavage fracture, and incorporation of a softer SiC phase appears to improve overall toughness [Zhao, Qian, Daemen et al. (2004)]. As with the other composite, possible reasons include crack deflection at grain and phase boundaries as well as compressive residual stresses in grains of the stiffer phase. Twinning is also prevalent in β-SiC as a result of its low stacking fault energy [Van Torne (1966); Ning and Ye (1990)].
This paper invokes, in FE simulations, a theory combining aspects of phase field mechanics and differential geometry on generalized Finsler spaces. Perhaps the most comprehensive description of this model, including versions accounting for multiple inelastic deformation mechanisms, finite strain kinematics, and nonlinear elasticity, is given in Clayton et al. [Clayton and Knap (2018a)]. A pseudo-Finsler metric tensor is assigned that depends on a vector of internal state variables, whose entries when rendered dimensionless upon division by a regularization length, are interpreted as order parameters in the context of phase field theory. Dependence of the metric on internal state permits changes in lengths, areas, and volumes of material in conjunction with microstructure evolution. For example, voids may lead to dilatation, or structure collapse in B 4 C may lead to densification [Clayton (2012); An and Goddard (2015)]. Incorporation of this metric tensor in the theory leads to a rescaled energy density and augmented balance laws that enable additional physics (e.g., pressure-shear coupling) without prescription of ad-hoc kinetic equations. The phase field-related components of the model are based on theory in Clayton et al. [Clayton and Knap (2011, 2014, 2016], while the Finsler-geometric concepts entering the theory follow from Clayton [Clayton (2017b[Clayton ( ,c,a, 2018a]. Potential advantages of the present computational implementation over other methods [Tomar, Zhai and Zhou (2004); Clayton (2005Clayton ( , 2009); Bammann and Solanki (2010)] for modeling inelasticity and fracture include the following: intrinsic regularization by gradient-dependent energy, little if any calibration of parameters, and preservation of mesh topology, e.g., no element deletion or separation along element boundaries.
A recent study by Clayton et al. [Clayton, Leavy and Knap (2019)] invoked a simpler version of the present theory towards study of the same two ceramic composites. That work invoked pure phase field mechanics [Clayton and Knap (2016)] without incorporation of the pseudo-Finsler metric. Furthermore, uniaxial compression loading was simulated in Clayton et al. [Clayton, Leavy and Knap (2019)], while pure shear (i.e., biaxial tension-compression) is simulated in the present study. Shear strength is of high interest as it is thought to be an indicator of resistance to ballistic penetration in armor ceramics [Bourne (2008); Clayton (2015Clayton ( , 2016].
Remaining sections of this work are organized as follows. The continuum theory and aspects of its computer implementation are reviewed in §2. Polycrystalline simulations, with results from shear loading, on B 4 C-TiB 2 are given in §3. Parallel treatment of diamond-SiC is contained in §4, followed by conclusions.

Finsler-geometric phase field theory
The continuum theory invoked in subsequent simulations is tersely summarized here. The present approach extends that of Clayton et al. [Clayton, Leavy and Knap (2019)] to include a generalized pseudo-Finsler metric tensor [Clayton (2017c)], rather than the Euclidean (and more specifically, rectangular Cartesian) metric invoked in Clayton et al. [Clayton, Leavy and Knap (2019)]. The present theory implements linearized kinematics [Clayton and Knap (2015)] and quasi-static balance laws following from a variational approach [Clayton (2017b)]. Geometrically nonlinear and dynamic models are discussed elsewhere [Clayton (2017a); Clayton and Knap (2018a)].

Internal state
A vector D of internal state variables coincides with a director vector of Finsler geometry, though here this state vector need not be of unit length. Let X with Cartesian components X K denote the reference position vector of a material particle. In a preferred coordinate system, Entries of D normalized by regularization length l are order parameters ξ and η. Parameter ξ ∈ [0, 1] depicts fracture: ξ(X) = 0∀X ∈ undamaged material, ξ(X) ∈ (0, 1)∀X ∈ partially degraded material, ξ(X) = 1∀X fully failed material. (2) Other structure changes are captured by η ∈ [0, 1]: η(X) = 0∀X ∈ parent elastic crystal, η(X) ∈ (0, 1)∀X ∈ structural boundary zone, η(X) = 1∀X ∈ structurally transformed state. (3) For ceramics considered herein, η > 0 can track shearing by slip or twinning dislocations and monitor localized slip in amorphous bands. It also can represent density variations due to phase changes.
Metric tensor G depends potentially on coordinates X as well as internal state D.
(4) Explicit dependence on X is not necessary for rectangular coordinates. Constants m and k quantify volume changes associated with respective fracture/cavitation and phase changes. The determinant is G = det G = exp(mξ 2 + kη 2 ).
The state-dependent contribution from changes in microstructure is [Clayton and Knap (2018a) (6) The first term on the right accounts for dilatation from cavitation or cracks, the second for compression or extension normal to a material plane from phase changes, the third for inelastic shearing. Vector M denotes the plane of structural transformations, slip, twinning, and/or localized shear, where the orthogonal shear direction is S. Total strain is , and elastic strain is E : Polynomial interpolation functions in (6) are [Clayton, Leavy and Knap (2019) Uniaxial strain normal to the basal plane commensurate with the phase change in B 4 C [Clayton (2012[Clayton ( , 2017c] is quantified by x η = exp( k 2 ) − 1. Twinning on habit plane M is a simple shear whose maximum magnitude is γ 0 . The same functions account for preferred single slip of full or partial dislocations along S with slip plane normal M, where γ 0 is a saturation limit [Clayton (2018b)].

Energy functional and equilibrium equations
A total energy functional for a heterogeneous body Ω with boundary ∂Ω is Strain energy per unit volume is W , phase energy is f , and their sum is ψ = W + f . In a reference state of pseudo-Finsler space, a volume element of the body is scaled as dV =

√
GdV 0 , where dV 0 is the element in Cartesian space. Properties may vary with X, e.g., among grains in a polycrystal. Position X is dropped from the explicit list of arguments in subsequent derivations, where it is understood that properties are homogeneous within each sub-volume of a grain [Clayton and Knap (2016)].
Internal state contributions, including gradient effects, are embedded in function f : (10) Denote the mechanical traction vector by t, the conjugate surface force to η by r, and the conjugate surface force to ξ by s. Stationary points of Ψ = ψdΩ are solutions to the following variational equation at fixed X: Euler-Lagrange equations resulting from this principle and application of Rund's divergence theorem [Rund (1975)] are the static balance of linear momentum and the static equilibrium equations for internal state: Right sides of (13) and (14) vanish identically for a Riemannian metric (m = k = 0) and in such cases coincide with phase field mechanics [Clayton and Knap (2016); Clayton, Leavy and Knap (2019)]. Stress P is symmetric: Natural boundary conditions on ∂Ω with outward unit normal n are t = P · n, r = 2κ κ κ : (∇η ⊗ n), s = 2ω ω ω : (∇ξ ⊗ n).
Strain energy density is linear elastic with damage-dependent moduli λ, µ, K: Accordingly, tangent elastic moduli µ and λ reduce in magnitude as ξ increases. Minimum values are ζµ 0 and ζλ 0 when ξ(X) = 1, corresponding to complete local fracture at X. In the numerical implementation, elements whose nodes contain values of ξ at or approaching unity are considered fractured, but these elements are never removed from the simulation since they contain a very small yet finite stiffness to prevent interpenetration. According to (19), the bulk modulus retains its full initial value K 0 = λ 0 + 2 3 µ 0 for compressive elastic volume changes such that unrealistic local collapse of the material is prohibited. Thus, under compressive pressure, behavior of the material becomes akin to that of a compressible elastic fluid. This approach of addressing failed elements is standard in phase field simulations of fracture [Borden, Verhoosel, Scott et al. (2012); Clayton and Knap (2014)].
Gradient-independent energy functions in (10) are, with H(·) a Heaviside function, The first term in f 0 is a double-well with minima at η = (0, η 0 ), a local maximum at η = η 0 /2, and a cut-off η 0 > 0. Quadratic forms are the remainder of f 0 and g 0 . Constants entering (20) are A, η 0 ,Â, and B. Let Γ denote twin boundary or stacking fault energy, Υ fracture surface energy. Material constants are related by ConstantÂ is an energy barrier for phase changes and can also account for energy of dislocation lines. Valuess for each crystal type follow later in §3 and §4.

Numerical methods
The FE method is used to solve 3D problems in an incremental fashion. Boundary conditions are updated at each load increment. Numerical algorithms apply conjugate gradient energy minimization [Clayton and Knap (2016)]. Solution fields u(X), η(X), ξ(X) minimize Ψ subject to boundary constraints, where the second of (9) is used. In this context, V 0 is the volume of the FE domain in its initial undeformed state in Euclidean space, prior to structural transformations. Force residuals computed numerically from Ψ must be smaller than a prescribed tolerance; (13) and (14) are not evaluated explicitly. Irreversibility is ensured by constraints δξ(X) ≥ 0 for ξ(X) ≥ ξ T and δη(X) ≥ 0 for η(X) ≥ η T , where ξ T = η T = 0.9. According to the present variational model based on energy minimization similar to Bourdin et al. [Bourdin, Francfort and Marigo (2000)], strain rate does not influence the results, which are effectively time-independent. Inertia and stress waves are excluded, as are kinetics of internal state that can be addressed using methods described elsewhere [Clayton and Knap (2018a)].

Boron carbide-titanium diboride polycrystals
Properties and FE renderings of B 4 C-TiB 2 are presented, followed by numerical results for pure shear loading.

Material representation
Material properties and FE models are discussed in brief, noting that detailed justification and sources of property data are available in Clayton et al. [Clayton, Leavy and Knap (2019)]. Parameters entering models of B 4 C and TiB 2 are summarized in Tab. 1.
First consider B 4 C. Twinning occurs on 1010 {0001} and fracture on {0001}. Amorphous shear bands can form inside twins. Twinning shear is γ 0 , and twinning surface energy is Γ.
Lattice collapse normal to (0001) produces x η , and x ξ measures expansion from cavitation [An and Goddard (2015)]. The crystal-to-glass energy barrier isÂ. Fracture energy is Υ and cleavage anisotropy isβ.
Now consider TiB 2 . Fracture occurs on {0001}, with Υ the surface energy [Du and Chen (2018)]. Order parameter η addresses stacking fault energy and line energy of (partial) dislocations via A andÂ, respectively. Dilatation from fractures is represented by x ξ > 0. Denote the Burgers vector of a partial dislocation by b, cumulative dislocation density by ρ D . The lowest energy barrier corresponds to 1120 {0001} slip [Du and Chen (2018)]. Stacking fault energy Γ is inferred from Vanderwalker et al. [Vanderwalker and Croft (1988)]. Letx denote the maximum distance a (partial) dislocation travels prior to arrest. Let ρ S denote the saturation dislocation density. Stored line energy is 1 Clayton (2011)], with full Burgers vector b 0 . The phase field representation of single partial slip in TiB 2 thus obeys Residual stresses may affect toughness [White and Dickey (2011); Scharf and Rubink (2018)]. These may be addressed in future work along with other defects.
Domains consist of 50 crystal polyhedra comprising a cube Ω of edge length L = 10 µm. Average grain size is L/50 1/3 ≈ 2.5 µm. Numerals 1 and 2 pertain to B 4 C and TiB 2 , respectively. Morphologies with volume fractions V 2 of 10% and 18% are simulated. These are shown in Fig. 1. Volume fractions are physically comparable to those of real microstructures [Scharf and Rubink (2018)]. Two random orientation distributions (S, M, N) are assigned to grain lattices of each morphology, such that four simulations are reported in total.  Displacement-controlled boundary conditions on part of ∂Ω correspond to pure shear [Clayton and Knap (2018b)], consisting of tension and compression applied in equal magnitudes along perpendicular Y -and Z-directions. Let¯ denote the load parameter. Essential boundary conditions are applied to four of the six faces of the initially cube-shaped polycrystalline aggregate: Along the remaining two faces of the domain (X = 0, X = L), traction-free conditions allow the body to expand or contract freely in the X-direction. Regarding conjugate forces to internal state, free boundary conditions r = s = 0 on all of ∂Ω. In the current implementation of these boundary conditions, average shear strain¯ is increased incrementally in steps of 5 × 10 −4 . Conjugate gradient energy minimization [Clayton and Knap (2016)] is used to seek an equilibrium configuration of the body at each increment.

Computational results
Microstructures for grains of lattice orientation set 1 are shown in Fig. 2 at an applied shear strain of¯ = 0.01 = 1%. Order parameter η denotes twinning, shear localization, and amorphization in the B 4 C phase and basal slip in the TiB 2 phase. A quadratic scale is used to visualize low transformation magnitudes (0 η 0.25). Elements with ξ > 0.8 are visually removed to show fractures. As explained following (18) and (19), such elements possess very low tangent stiffness, with µ → ζµ 0 as ξ → 1 in shear for example, where herein ζ = 0.01. These elements are never fully removed from the calculation; original mesh topology is retained.
Fracture paths are quite similar for microstructures 1 and 2 in respective Fig. 2(a) and Fig. 2(b). Cleavage cracks tend to traverse the more brittle B 4 C grains. Transformation behaviors (η > 0) are enhanced in microstructure 2 relative to microstructure 1, notably in grains of the TiB 2 phase which have a lower stacking fault energy (0.12 versus 0.54 J/m 2 ) and thus a lower resistance to slip than the twinning threshold of B 4 C. In both of these cases, and in others not shown in Fig. 2 for orientation set 2, inelasticity mechanisms of slip, twinning, and amorphization are relatively scarce (η 1 in most grains), in agreement with experimental observations from samples of comparable sizes [Clayton and Knap (2018a); Scharf and Rubink (2018)].
Volume-averaged maximum shear stress and order parameters are labeledP ,η, andξ, respectively. Definitions of the latter two are obvious, while the stress averaged over the reference domain is Values are shown versus applied shear strain¯ in Fig. 3. Average stresses in Fig. 3(a) demonstrate modest effects of lattice orientation and volume fractions, recalling microstructure 1 vs. 2 in Fig. 1 with 10% and 18% TiB 2 , respectively. Stiffer responses are obtained for the same lattice orientation set when the fraction of TiB 2 is increased.  This results from the higher elastic shear modulus of TiB 2 as well this phase's higher fracture energy, greater ductility due to slip, and lower tendency to cleave (Tab. 1). Average damageξ in Fig. 3(b) evolves similarly for all cases except microstructure 1 with lattice 1. Interestingly, this is the structurally weakest case in terms of average shear stress, yet it demonstrates the lowest accumulation rate of average micro-crack density. The same trend has been observed in simulations of pure B 4 C [Clayton and Knap (2018a,b)]: overall stiffness reduction and failure are most closely correlated with the presence of a few dominant cracks and are inversely correlated with diffuse crack density. Average structural changes are reported in Fig. 3(c). Microstructure 2 demonstrates more inelasticity than microstructure 1 since the former contains more TiB 2 that slips more easily than shear-induced inelastic modes (twinning and localization) in B 4 C.
Tab. 2 lists average maximum shear stress and values of applied shear strain and averaged order parameters at the same peak loading increment. Peak stress P C increases from 1% to 10% when the volume fraction of TiB 2 is increased from 10% to 18%. Transformations η C increase with V 2 for reasons outlined in the context of Fig. 3(c). Trends are less obvious for ductility C and peak damage ξ C . Different trends are evident in these variables for lattice 1 and lattice 2. Lattice orientation affects overall ductility and damage via the action of resolved driving forces for inelasticity and fracture that depend on orientations of basal planes for slip, twinning, and cleavage. Differences due to lattice orientation can be comparable in magnitude to differences due to volume fraction of the second phase.
Experiments suggest fracture toughness of polycrystalline B 4 C may be improved by 10-100% and flexure strength by 10-30% via addition of 13-23% TiB 2 [Scharf and Rubink (2018)]. The current simulations show more modest improvements with increasing TiB 2 . Residual stresses omitted in the simulations may partially explain the discrepancy. Since the FE domain may not be large enough to capture long-crack blunting by TiB 2 , it may be enlarged to include more crystals in the future.

Diamond-silicon carbide polycrystals
Properties and FE representations for diamond (C)-β-SiC are reported, followed by simulation results for pure shear loading.

Material representation
Model features for each crystal type are summarized briefly here; details and supporting references can be found in Clayton et al. [Clayton, Leavy and Knap (2019)]. Properties are given in Tab. 3. First consider diamond. Cleavage occurs on {111}, with surface energy Υ. Dilatation from crack opening is quantified by x ξ > 0, with the same value used for all ceramics of present interest. Silicon carbide can fracture and twin. Cleavage planes are {110}. Twinning in the zinc blende structure occurs on 112 {111} with Van Torne (1966)]. Stacking fault energy Γ is notably low [Ning and Ye (1990)]. Regularization length l is universal among all ceramics addressed. As in §3, residual stresses and other microstructure heterogeneities (initial pores, tertiary phases, etc.) are deferred to future work.
Two very different morphologies are modeled, as shown in Fig. 4. Both domains are cubes of the same dimensions of §3, containing 50 polyhedral grains of average size L/50 1/3 ≈ 2.5 µm. Denote by V 1 and V 2 volume fractions of diamond and SiC phases, respectively. Microstructure 1 features a layer of uniform thickness (200 nm) of nano-crystalline SiC fully surrounding each diamond crystal, contributing V 2 = 20%. Microstructure 2 contains  Boundary conditions for pure shear are prescribed nearly identically to those discussed in §3. For each of the two orientation sets with microstructure 1, one simulation is performed with boundary conditions of (23) and¯ ≥ 0, and a second with the sign of¯ reversed.
In the forthcoming discussion of results of each of these second simulations,¯ refers to the magnitude of (negative) applied shear strain for ease of comparison. Twinning in the nanocrystalline layer is promoted in the former case but precluded in the latter since the driving force, i.e., resolved shear stress for twinning, becomes negative. For microstructure 2,¯ ≥ 0 only. The total number of simulations reported in §4 is six. Fig. 5 are contours of η that show twinning in the SiC phase. Results correspond to lattice orientation set 1, where the habit plane and twinning direction are oriented to promote twinning in the nano-crystalline grain boundary (GB) layer in Fig. 5(a). Quadratic scales enable visualization of low to moderate, i.e., incomplete, twinning (0 η 0.25) that occurs in this GB layer and in a few preferentially oriented bulk grains of SiC in Fig. 5(b). Physically, incomplete twinning corresponds to planar stacking faults typical in β-SiC. Diamond, which comprises the majority of each microstructure, does not undergo twinning or slip, so η = 0 in all diamond crystals. Fracture behaviors are notably different for the two cases shown, where cleavage cracks in the diamond grains are evident at different locations in each microstructure. A dominant cleavage crack cutting diagonally across multiple diamond crystals through entire domain is noteworthy in microstructure 2, while smaller cracks arrested or deflected at phase boundaries are evident in microstructure 1. Intergranular fractures are more common in the latter microstructure.

Shown in
Volume-averaged shear stressP is defined as in (24), and volume-averaged order parametersη andξ are calculated similarly. Values are shown versus applied shear strain¯ in Fig. 6. Stresses for all six simulations of diamond-SiC are reported in Fig. 6(a). Recall that for two of the simulations with microstructure 1, the twin system in the GB phase is oriented with respect to the loading mode such that twinning is precluded. Lattice orientations of diamond grains are maintained identically with their counterparts in the legend immediately above. Initial elastic stiffness is larger in microstructure 2, which has a lower volume fraction than microstructure 1 (V 2 = 0.10 vs. 0.20) of the elastically stiffer diamond phase. However, maximum stresses are often lower for microstructure 2 and always occur at a lower applied strain for microstructure 1. When twinning is suppressed in microstructure 1, effective stiffness increases since local deformation in the GB phase remains elastic until fracture. No trends are discernible for behaviors in the strain softening regime, which differ substantially among microstructures and orientations. The same can be inferred with regard to average damageξ in Fig. 6(b), and correlations with average stress are not all obvious. The rate of damage accumulation does decrease, however, at strains beyond those corresponding to maximum average stress. Transformation behavior quantified in Fig. 6(c) is logically explained. Microstructure 1 demonstrates significant activity due to the preferential orientation of its twinning system in the nano-crystalline GB layers. Minor transformation behavior is reflected byη for microstructure 2 as a result of twinning in a few bulk grains, e.g., see Fig. 5(b). When twinning is suppressed in microstructure 1,η → 0 since the diamond phase is elastic-brittle.
Tabulated results for peak shear stress, corresponding ductility or peak strain, and order parameters are found in Tab. 4. Peak shear stress P C decreases by up to 10% when the GB phase is able to twin in microstructure 1, while ductility is about the same or slightly improved with twinning. Compared to microstructure 1 which contains half as much of the weaker SiC phase by volume (10% versus 20%), microstructure 2 demonstrates lower peak shear stress than three of the four simulations with microstructure 1. Ductility¯ C is    lowered by 15 to 30% in microstructure 2 compared to microstructure 1. Transformation measured by η C is much larger in microstructure 1 (when twinning is enabled) than 2 since it contains a larger percentage of SiC preferentially oriented for twinning. Peak damage ξ C is lower in microstructure 2 than 1, suggesting dominant cracks in the former and diffuse damage in the latter, more durable aggregates. The present findings agree qualitatively with experiments [Zhao, Qian, Daemen et al. (2004)]: fracture toughness of polycrystalline diamond is improved by intergranular nano-crystalline SiC layers.
As was also the case for B 4 C-TiB 2 reported in §3, trends in results for diamond-SiC depend on lattice orientation. In the present material, fracture is strongly affected by orientation of octahedral cleavage planes in the diamond phase, and ductility by orientations of twinning habit planes in the SiC phase. As an exception to other trends implied already, for lattice 2 with slip/twins active, P C is lower in microstructure 1 than microstructure 2. Differences in strength between microstructures due to phase volume fractions and interface morphology depend strongly on enabling of twinning activity or lack thereof. This suggests important coupling exists among mechanisms of twinning, interfacial failure, and cleavage fracture.

Conclusions
A constitutive theory incorporating ideas from Finsler differential geometry and phase field mechanics has been implemented to model deformation and failure of dual-phase ceramics under pure shear loading. Simulations of polycrystalline B 4 C-TiB 2 and diamond-SiC composites demonstrate different twinning, amorphization, slip, and fracture behaviors of these ceramics depending on crystal type, lattice orientation, and phase morphology. Addition of TiB 2 to B 4 C provides an increase in overall shear strength and ductility, where the former phase is more likely to shear inelastically and possesses a higher fracture energy, thereby enabling blunting of cleavage cracks initiated in the more brittle B 4 C phase. The SiC phase, when incorporated as nano-crystalline layers along grain boundaries, apparently obstructs paths of cleavage cracks in brittle diamond crystals. For improved failure resistance, encapsulation of diamond in very fine grains of the softer second phase is recommended over addition of SiC grains of commensurate size of the diamond. For quasi-static deformations of aggregates of 50 bulk crystals with grain sizes on the order of several µm, maximum shear strengths are on the order of 1.5 GPa for B 4 C-TiB 2 and 2.9 GPa for diamond-SiC. Respective applied shear strains at maximum load are on the order of 0.9% and 0.6% for the two material systems.