Phase Field Study of Ferroelastic Toughening Mechanisms of Polycrystalline t^′-YSZ

1  Introduction

Yttria-stabilized zirconia containing 6–8 wt% Y2O3 (YSZ) serves as the dominant surface coating material for thermal barrier coating systems. This material gains widespread application owing to its superior comprehensive properties, including outstanding low heat conduction performance, high oxygen permeability, and favorable coefficient of thermal expansion [1–5]. Although several candidate oxides exhibit lower thermal conductivity, YSZ coatings demonstrate superior resistance to thermal cycling degradation, which is attributed to their relatively higher fracture toughness [6–8]. The exceptional fracture toughness of YSZ stems from its inherent ferroelastic toughening mechanism, triggered by stress-driven domain reorientation in the vicinity of crack tips. Numerous research endeavors in recent years have focused on developing advanced thermal barrier coating materials with reduced thermal conductivity and elevated temperature tolerance [1–5]. Nevertheless, ferroelastic domain switching remains a critical mechanism to achieve high fracture toughness in newly synthesized oxides [9,10]. A comprehensive understanding of the microstructural effect on ferroelastic toughening is therefore essential for the rational design of advanced top-coat materials.

The as-deposited YSZ top coat exists in the metastable tetragonal (t′) phase. In contrast to the equilibrium tetragonal phase (t phase), the t′-phase exhibits a higher yttria content and reduced tetragonality. Under high-temperature working conditions, the conventional toughening behavior induced by tetragonal-monoclinic (t-m) phase transition becomes thermodynamically inhibited [11,12]. In addition, the t-m transformation is associated with severe shear strain and volumetric expansion, which can induce microcracking within the top coat [13]. Thus, the t-m phase transformation is generally considered to be detrimental to the integrity of TBCs. The exceptional fracture toughness of t′-YSZ coatings comes from its special ferroelastic deformation-induced toughening behavior [14,15]. The tetragonal prime phase of YSZ has three variants with orthogonal spatial orientations. Under mechanical loads, the t′ variants may experience reorientation to accommodate deformation. The domain reorientation (or domain switching) is displacive and is accompanied by eigenstrains. Domain switching (DS) occurring around the crack tip can produce a crack growth shielding effect. DS of t′-YSZ under mechanical loads has been experimentally confirmed by numerous researchers. To name a few, Chan et al. [16] and Srinivasan et al. [17] discovered domain reorientation in compression tests of polydomain YSZ samples. Mercer et al. [14] identified the nucleation of new domains in the wake of indentation cracks in polycrystalline YSZ with mono domain structure. It was inferred that the toughening of the indentation cracks was controlled via new domain formation rather than the migration of original domain boundaries. Bolon and Gentleman [18] observed ferroelastic switching in 18-mol% ceria-stabilized zirconia by using polarized confocal Raman spectroscopy. Domain rearrangement behavior is detected both inside and around indentation cracks. Li et al. [19,20] observed domain reorientation in YSZ thin films with indentation tests. Masuda et al. [21] investigated the orientation-related ferroelastic response and plastic deformation characteristics of t′-phase zirconia pseudosingle crystals with nanodomain microstructures consisting of three different t′ variants. The switching of the domains was monitored with electron microscopy. Luo et al. [22] experimentally identified the ferroelastic domain switching behavior and the associated toughening mechanism in yttrium tantalate-zirconia oxide. Carbogno et al. [23] and Chen et al. [24] validated the ferroelastic structural reorientation mechanism of tetragonal zirconia with density functional theory.

Although the DS behavior of YSZ has been experimentally and theoretically validated, quantitative studies on its toughening effect on crack growth are rarely reported in the literature, especially when the microstructural effects are considered. Actually, the test specimens used in the aforementioned experimental studies possess different domain structures. For example, the tests carried out by Baither et al. [25] and Masuda et al. [21] were based on polydomain single crystals. The polydomain structure manifests as the consequence of the c-t′ phase transformation. By contrast, Mercer et al. [14] conducted indentation tests on single-domain polycrystalline YSZ, where each grain possesses only one YSZ domain. Accordingly, it is of great significance to explore how diverse domain configurations affect the ferroelastic toughening performance.

As a mesoscale approach, the phase field (PF) method is widely recognized for its robustness in characterizing microstructure evolution in solid materials, including phase transformation, ferroelectric domain switching, dislocation slip, and crack propagation [26–29]. It is noted that ferroelectric domain switching and its toughening effect have been intensively studied with the PF approach [30–33]. The PF model for ferroelectric domain switching usually adopts the polarization vectors as the order parameters [34–36]. The free energy functional is constructed based on the Landau-Ginzburg-Devonshire theory [34–36]. However, the cubic, tetragonal (including the t′ phase), and monoclinic phases of YSZ are all non-polar. Thus, the traditional PF model for ferroelectric domain switching is ill-suited for characterizing DS in YSZ.

Recently, Pi et al. [37] and Li et al. [20] developed a PF model for DS in YSZ single crystals. With the above PF model, Pi et al. [38] studied the evolution of the t′ variants around a static crack tip. Bhattacharya and Asle Zaeem [39,40] employed PF modelling to study domain nucleation and the kinetics of domain switching with consideration of back-switching. Zhou et al. [41] proposed a thermodynamically consistent phase-field theory incorporating strain gradient elasticity and interfacial effects to investigate micro- and nano-scale ferroelastic variant switching in t′-YSZ. This PF model is utilized to analyze domain nucleation and evolution behaviors at pre-existing crack tips in single-phase t′-YSZ. Li and Li [42] constructed a phase-field framework to explore domain morphological evolution and mechanical behaviors of polycrystalline t′-YSZ subjected to external loads. The domain evolution around a prefabricated crack tip was also studied. These studies predominantly focused on the nucleation and evolution of domains within stress fields. However, none of the aforementioned studies have explored the toughening effect induced by domain switching on crack propagation.

Actually, the PF method can simulate phase transformation, domain switching, and fracture in a unified framework [33,34]. Zhao et al. [43] integrated the phase-field method for tetragonal-monoclinic phase transition with phase-field fracture theory to probe crack propagation characteristics in tetragonal zirconia. Afterwards, this coupled phase-field framework was further expanded to analyze transformation-triggered microcracking behaviors within polycrystalline tetragonal zirconia materials [13]. Jafarzadeh et al. [44] developed a phase-field method for analyzing nanoscale coupling of fracture and phase transformation. Moshkelgosha and Mamivand [45,46] combined the PF model for tetragonal-monoclinic phase transformation with the variational PF theory for fracture to study the fracture of single crystal and polycrystal tetragonal zirconia. The perturbation of the crack propagation path by the transformed phase was fully discussed. Based on a similar approach, Borzabadi Farahani et al. [47] studied crack nucleation and growth in martensitic phase transformation materials. Lotfolahpour et al. [48] proposed a modified PF model to study the interaction between phase transformation and crack propagation in superelastic shape memory ceramics. The PF model was extended to include the plastic slip mechanism and was applied to study the orientation-dependent deformation and failure of micropillar shape memory ceramics [49]. The above studies focused on the interaction mechanisms between martensitic phase transformation and crack propagation. Sun et al. [50] constructed a phase-field model to investigate ferroelastic toughening of single-domain monocrystalline t′-YSZ. They introduced the energy dissipation rate (EDR) for the first time to quantify the toughening mechanism derived from domain switching. Subsequently, Xiong et al. [51,52] employed this method to explore the toughening characteristics of shape memory alloys.

In this paper, the ferroelastic toughening behavior of polycrystal t′-YSZ is investigated by employing the PF approach developed in our previous study [53]. The impact of the initial domain structure on the toughening effect is discussed through a comparative analysis of two distinct polycrystals. In polycrystal I, each grain is characterized by a single t′-phase domain. In contrast, polycrystal II exhibits a polydomain microstructure, which forms as the product of the cubic to tetragonal prime phase transformation. The evolution of the domain structure in the two polycrystals is simulated using a unified phase field model. Subsequently, the phase-field model describing domain switching is coupled with phase-field fracture theories [54,55] to analyze crack propagation features in polycrystalline t′-YSZ. The ferroelastic toughening effect exerted by domain switching on local crack expansion is quantitatively evaluated via energy dissipation rate calculation. Moreover, comprehensive discussions are carried out to clarify how initial domain configurations affect crack growth characteristics and the ferroelastic toughening effect.

2  Theoretical Foundation

This part elaborates on the theoretical foundation of the phase-field model for domain switching (DS). Afterwards, the coupling scheme between this domain-switching phase-field model and phase-field fracture theory is elaborated in detail.

2.1 PF Model of Domain Switching for Polycrystalline t′-YSZ

To simplify the analysis, this study employs the two-dimensional (2D) plane strain assumption. For the 2D case, the tetragonal prime phase possesses two variants, which are oriented perpendicularly to each other. The two variants are represented by η=±1 respectively, where the a-axis of variant I (η=−1) is along the horizontal direction and that of variant II (η=1) is along the vertical direction. The unstable cubic phase corresponds to η=0 in the PF model. In order to investigate the microstructure effect on the ferroelastic toughening, two polycrystal microstructures are considered.

In polycrystal I, every grain constitutes an individual domain of the t′ phase. In polycrystal II, the c-t′ phase transformation gives rise to a polydomain structure, where the orthogonal variants form twins to relieve the deformation energy. It should be mentioned that twinned t′ variants have been identified in air plasma-sprayed (APS) coatings [15]. However, the formation of twinned microstructure is not prevalent in APS coatings. Thus, the two polycrystal microstructures mentioned above can co-exist in APS coatings in certain circumstances.

The domain switching behavior in two polycrystals is described with a unified model. The system energy consists of the elastic strain energy, the chemical energy, and the gradient energy [56,57]:

F=Fel+Fch+Fgr.(1)

The chemical energy density is written as [50,53]:

fch(η)=Gb(η4−2η2+1).(2)

The variation of chemical energy density vs. order parameter (OP) is shown in Fig. 1. The two minimum points at η=±1 are corresponding to the two variants of the t′ phase, respectively. The maximum point at η=0 is corresponding to the cubic phase. In polycrystal I, at the onset, each grain is exclusively populated by one of the two t′ variants. Domain switching can take place if the driving force overcomes the energy barrier Gb. Polycrystal II commences in the cubic phase. Under rapid cooling, the a phase transformation occurs, accompanied by the simultaneous formation of t′ twins. DS may occur under external loading, manifested as the migration of domain boundaries. The critical stress for DS is correlated with the value of Gb. As elaborated in our earlier research [53], the predicted coercive stress for domain switching is about 400 MPa when Gb is 3×106J/m3, which is close to the experimental result [25]. Gb is set at 3×106J/m3 in the following studies.

Figure 1: The chemical energy density [53].

The gradient energy density and complete elastic energy density are expressed as [57]:

fgr(η)=β2|∇η|2,(3)

fel(η,u)=12Cijklεijelεklel,(4)

with

εijel=εij−εij∗=12(ui,j+uj,i)−εij∗,(5)

where β is the gradient energy coefficient; εij is the total strain; εijel represents the elastic strain; εij∗ denotes the eigenstrain for DS, which is dependent on the crystal orientation. For polycrystal t′-YSZ, the eigenstrain can be written as [50,53]:

εij∗=Rik(θ)Ril(θ)εkl∗0,(6)

where θ is the orientation angle of the a-axis; Rij(θ) denotes the rotation matrix [57]:

[Rij(θ)]=[cos⁡θsin⁡θ−sin⁡θcos⁡θ].(7)

εkl∗0 represents the eigenstrain when the a-axis is along the horizontal direction. For the mono domain t′ phase, the eigenstrain for domain reorientation from variant I to variant II is written as [53]:

[εkl∗0]=[εs00−εs].(8)

Here, εs=(c−a)/a, where a and c denote the lattice parameters of the tetragonal prime phase. According to Sun et al. [50], εs takes 0.0232 in this study. As variants I and II are represented by η=−1 and η=1 respectively, Eq. (8) can be modified as:

[εkl∗0]=[(η+1)εs/200−(η+1)εs/2].(9)

For t′-YSZ with polydomain structure, the eigenstrain is calculated by designating the cubic phase as the reference state. εkl∗0 can be written as:

[εij∗0]=[ηεs/200−ηεs/2].(10)

Combing Eqs. (6), (9) and (10), we have:

{εijel=12(ui,j+uj,i)−(η+1)εijs(single domain)εijel=12(ui,j+uj,i)−ηεijs(multi domain),(11)

with

[εijs]=[R(θ)][εs/200−εs/2][RT(θ)].(12)

where [R(θ)] is the rotation matrix.

In the PF modelling, the OP η usually experiences numeric overflow, i.e., its value may exceed the expected range [−1, 1]. Fang et al. [53] showed that this problem can be resolved by using the truncation function pw(⋅) in COMSOL. It is noted that Rajendran et al. [58] used a similar approach to simulate martensitic transformation. Eq. (11) is thus written as [53]:

{εijel=12(ui,j+uj,i)−[pw(η)+1]εijs(single domain)εijel=12(ui,j+uj,i)−pw(η)εijs(multi domain)(13)

with

pw(x)={−1(x≤−1)x(−1<x<1)1(x≥1)(14)

The truncation function is smoothed below to enhance numerical stability [53]:

d[pw(x)]dx={0(x≤−1)1(−1<x<1)0(x≥1)(15)

2.2 PF Model Coupling Fracture and Domain Switching

In t′-YSZ, the high stress around crack tips triggers domain switching, further affecting fracture characteristics. To simulate this coupled physical behavior, a phase-field coupling model is established by integrating the domain switching phase-field method and variational fracture theory. The total energy of this model is defined as [53–55]:

F=Fch+Fgr+∫Ωψε(η,ε,ϕ)dΩ+∫ΩGc(ϕ22l+l2|∇ϕ|2)dΩ−∫Ωb⋅udΩ−∫∂Ωtt⋅udS,(16)

where ϕ is the crack phase field; l denotes the length scale parameter of the crack phase; Gc represents the intrinsic fracture toughness. The third component stands for degraded elastic strain energy, while the fourth corresponds to fracture energy. The final two terms denote work induced by body force and surface force separately.

The degraded elastic energy density is expressed as [53]:

ψε(ε,η,ϕ)=12g(ϕ)Cijklεijelεklel=12[(1−ξ)(1−ϕ)2+ξ]Cijklεijelεklel,(17)

where g(ϕ) is the degradation function; Cijkl denotes the fourth-order elastic tensor; ξ is a dimensionless parameter with small values, which is introduced to improve the numerical stability. We assume ξ=10−6 in the numerical studies.

The governing equations of the phase fields can be derived with the Ginzburg-Landau theory [56,57]:

η˙=−LδFδη(r,t)(18)

φ˙+M[Gclφ−GclΔφ+g′(φ)ψε]=0,(19)

where L and M are the mobility parameters; Δ is the Laplacian operator. Eq. (19) contains a viscous term, which is helpful to improve the numerical stability [59,60].

The elastic strain energy should be decomposed to properly characterize the crack propagation path [61–63]. Here, we adopt the energy decomposition method proposed by Amor et al. [62]. Only the tensile and deviatoric parts of the strain energy are effective for the evolution of the crack phase field:

ψε(ε,ϕ)=Kn2⟨tr(εel)⟩−2+g(ϕ)[Kn2⟨tr(εel)⟩+2+μ(εDel:εDel)](20a)

ψε+=Kn2⟨tr(εel)⟩+2+μ(εDel:εDel).(20b)

In Eq. (20), εDel denotes the deviatoric part of the elastic strain; μ represents the shear modulus; Kn stands for the n-dimensional bulk modulus and ⟨x⟩±=(x±|x|)/2.

The irreversible condition for the evolution of the crack phase field is enforced by introducing an auxiliary history field [59,60]:

H(x,t)=maxs∈[0,t]ψε+.(21)

The history field H(x,t) satisfies the Kuhn-Tucker conditions for loading and unloading [59]:

ψε+−H≤0,H˙≥0,H˙(ψε+−H)=0.(22)

Replacing ψε with H(x,t) in Eq. (19), the evolution equation of the crack phase field is written as:

φ˙+M[Gclφ−GclΔφ+h′(φ)H]=0.(23)

The equilibrium equation is solved simultaneously to obtain the displacement and stress fields:

∇⋅σ=0,(24)

where the Cauchy stress is given by:

σ=∂ψε∂ε=Kn⟨tr(εel)⟩−I+g(ϕ)[Kn⟨tr(εel)⟩+I+2μεDel].(25)

Eqs. (18), (23), and (24) are implemented in Comsol Multiphysics and solved iteratively. The staggered solution strategy is adopted to solve the multi-physical governing equations.

3  Results and Discussions

As shown in Fig. 2, a square plate with dimension 10µm×10µm containing an edge crack of length 5.4µm is considered. The crack tip is surrounded with a polycrystal region with dimension 5µm×6µm. Outside of the polycrystal region, the material is assumed to be linearly elastic and nontransformable. The polycrystal structure is generated by a Voronoi script, which consists of 49 grains with an average grain size of 1 μm. The grain boundaries (GBs) are 20 nm in thickness and are inert with respect to domain evolution. The inert grain boundary model is adopted to characterize the inhibitory effect of grain boundaries on domain switching. Nevertheless, grain boundaries are still allowed to undergo deformation and fracture. The elastic modulus of the grain boundary zone is assumed to be slightly lower than that of the bulk grains. This assumption is based on the fact that atoms at grain boundaries typically exist in a disordered state with high energy, which inherently weakens the mechanical response of the grain boundary region. The grain orientations are randomly distributed in the range [−45∘,45∘], as shown in Fig. 2b.

Figure 2: The numerical model: (a) The geometry of the model and the boundary conditions; (b) The grain orientation of each grain.

The growth of the edge crack and the accompanying domain structure evolution are simulated with the PF model introduced in the previous section. We assume that the crack propagates forward in a quasi-static manner. The mobility parameters are set as L=2000m3/J/s and M=0.5m/J/s, respectively. A tearing displacement load u∗ is exerted at the left side of the model to drive the crack growth. The right side of the model is hinged. The model parameters and material constants are listed in Table 1. The PF model presented in Section 2 is solved with the finite element method. The region through which the crack propagates is meshed with a refined finite element grid, where the smallest resolved mesh size is 5nm.

3.1 Crack Growth in Polycrystal I

In this section, the crack growth behavior in polycrystal I is studied. For polycrystal I, each grain is composed of a mono domain of the t′ phase. We assume that the material is initially composed of the domain represented by η=−1. Under external loading, the initial domain can switch to the orthogonal tetragonal domain represented by η=1. Firstly, the crack growth behavior under a monotonically increasing displacement load is studied. The displacement load is initially applied at a rate of 100nm/µs until u∗=50nm. After that, the displacement load increases at a lower rate (5nm/µs) to improve the numerical stability. The crack growth and the domain structure evolution are shown in Fig. 3.

Figure 3: Crack growth and domain switching in polycrystal I: (a) t=0µs an u∗=0nm; (b) t=0.3µs and u∗=30nm; (c) t=0.35µs and u∗=35nm; (d) t=0.5µs and u∗=50nm; (e) t=6µs and u∗=77.5nm; (f) t=7µs and u∗=82.5nm; (g) t=8µs and u∗=87.5nm; (h) t=9.6µs and u∗=95.5nm.

In Fig. 3 and the subsequent colored figures, the color bar denotes the value of the order parameter η, while the curved white lines represent the crack profile. As depicted in Fig. 3a,b, as the load increases, the new t′ domain represented by η=1 starts to nucleate at the crack tip. The new domain grows to the adjacent grains with the increase of the load and the advancement of the crack tip. Besides the crack tip, the new domain can also nucleate at the grain boundary or at the triple junction of grains, as shown in Fig. 3e–h. The overall domain pattern is influenced by the geometry and orientations of the grains near the crack tip. Fig. 3e,f show that a portion of the nucleated new domain switches back to the original variant with the advancement of the crack tip. The back-switching of the ferroelastic domain was recently studied by Bhattacharya and Asle Zaeem [40]. However, most of the nucleated new domains are still retained in the material and left in the wake of the growing crack. With the increase in the load and the forward movement of the crack tip, the proportion of the new domain keeps growing. Meanwhile, the nucleated new domains deflect the crack growth path. The crack path thus becomes tortuous with the extension of the crack tip.

Fig. 4 presents the impact of DS on the crack growth path. Fig. 4a shows the tortuous crack path disturbed by the nucleated new domains. If we intentionally ignore DS, the crack growth path is depicted in Fig. 4b. Fig. 4b shows that the crack propagates nearly straight forward when DS is switched off. The minor fluctuations in the growth path are attributed to the difference in the elastic moduli between the grain boundaries and the grain interiors. As a tortuous crack growth path facilitates more energy dissipation, crack deflection is usually regarded as a fracture toughening mechanism [64,65]. Thus, Fig. 4a suggests that DS in t′-YSZ can induce a fracture toughening effect by deflecting the crack growth path.

Figure 4: Comparison of the crack growth paths in Polycrystal I with consideration DS (a) or ignoring DS (b).

The domain reorientation activated around the crack tip is associated with eigenstrains, which yield a toughening effect on crack propagation [35,36]. Domain switching-induced toughening can be quantitatively evaluated by computing energy dissipation rate (EDR) [50,53,66,67], which is defined as the energy consumed per each crack growth length. In this work, the toughening effect characterized by the EDR is termed as the local toughening effect, which is differentiated from the toughening effect associated with crack deflection.

The EDR is calculated with the following formula [50,53]:

D=limδs→0(−δU−δWδs)≈−ΔU−ΔWΔs,(26)

where Δs denotes the increment of the crack length; ΔW and ΔU denote the increment of the external work and the elastic energy at each crack growth step, respectively. The elastic energy and the external work are calculated with the following formula [50,53]:

U=∫Ωψ(ε,ϕ)dΩ,(27)

W=∫Ωb⋅udΩ+∫∂Ωt⋅udS.(28)

An incremental tensile displacement load is applied to the numerical model to facilitate the calculation of the EDR. The displacement load increases monotonically until u∗=u0. After that, the displacement load increases incrementally. At each increment step, the load is held constant for a period of time so that the domain structure can reach the equilibrium state. At each step, the new crack tip is identified by setting ϕ=0.5. The crack growth length is then calculated with respect to the new crack tip. A small load increment is adopted in the simulation process so that the crack only propagates forward for a very short distance at each increment step. The increments of the elastic deformation energy and the external work are calculated at each load step. The EDR is then calculated based on Eq. (26). For crack growth in polycrystal I, the initial load u0 and the load increment Δu are 75 and 0.2 nm respectively. The initial load is applied at a rate of 75mm/s and is hold constant for 2µs. In the following incremental steps, the incremental load is applied at a rate of 0.4mm/s and is hold constant for 1.5µs. The evolution process of the domain pattern is similar to that presented in Fig. 3. Due to the strong nonlinearity of the multi-physical coupling problem and the small load increment, the computational cost is very high. The EDR is thus only calculated for a limited crack growth length to ensure the computational cost remains within tolerable limits. EDR vs. crack growth length is illustrated in Fig. 5, while corresponding domain morphologies are presented in Fig. 6.

Figure 5: The EDR for crack growth in polycrystal I with consideration of DS.

Figure 6: The domain patterns corresponding to different crack growth stages: (a) cl=864.4nm; (b) cl=932.5nm; (c) cl=1080.0nm; (d) cl=1131.7nm; (e) cl=1434.3nm; (f) cl=2060.1nm. (a–f) corresponds to points A1–A6 marked in Fig. 5, respectively.

For the PF fracture approach, the effective intrinsic fracture toughness is influenced by the mesh size and the regulation parameter [55]. According to the discussions by Bourdin et al. [55], the effective fracture toughness is estimated as Gc(1+h/2l), where h is the mesh size and l is the regulation parameter. In this study, l is 20 nm and h is about 5 nm along the crack propagation path. Thus the effective intrinsic fracture toughness is about 7.6N/m. In Fig. 5, the effective intrinsic fracture toughness is denoted by the horizontal line for comparison.

In Fig. 5, the predicted EDRs exhibit numerical oscillations as they are calculated at discrete load steps. Fig. 5 indicates that the average EDR exceeds the inherent fracture toughness at the initial crack propagation stage, proving domain reorientation produces a notable toughening impact. This effect is closely associated with new domain nucleation near crack tips (Fig. 3a–e). However, starting from point A1, the EDR decreases and even falls below the intrinsic fracture toughness between points A2 and A3. The domain corresponding patterns are shown in Fig. 6a–c, respectively. A comparison of these figures clearly shows that back-switching of the nucleated new domain occurs, which accounts for the decline in EDR from point A1 to point A3. Comparing with Fig. 6c, a new domain reorientation nucleus can be observed below the crack in Fig. 6d, which can reasonably explain the increase of the EDR from point A3 to point A4 shown in Fig. 5. Similarly, the back-switching of the new domain observed in Fig. 6e and the expansion of the domain reorientation zone in Fig. 6f account for the decline and increase of the EDR from point A4 to point A6 in Fig. 5, respectively. Fig. 5 also shows that the crack growth rate undergoes abrupt changes at distinct load steps, which is attributed to the complex interaction between crack growth and domain evolution. Especially, from point A5 to point A6, the load increment is only 0.2nm, but the crack length and domain pattern exhibit a great difference, as shown in Fig. 6e,f. The results depicted in Figs. 5 and 6 clearly indicate that the domain structure evolution distinctly influences the growth path, growth rate and local growth resistance of the crack. Due to the strong nonlinearity of the multi-physical coupling problem and the small load increment, the computational cost is extremely high when calculating the EDR. Consequently, only a limited crack growth length is simulated. Fig. 5 demonstrates that DS induces an overall toughening effect on crack growth.

It should be noted that the predicted EDR reflects the local toughening effect induced by domain switching, which may be lower than the experimentally measured macroscopic fracture toughness of YSZ. The actual microstructures of t′-YSZ coatings are far more complex, containing inherent internal defects such as interfaces, micropores, and microcracks. The evolution of these microdefects dissipates additional energy and acts as an extra toughening mechanism.

Furthermore, energy dissipation mechanisms, including plastic slip, are not considered in the current work. In addition, only two tetragonal variants, rather than three, are adopted in the present model. Introducing the above factors helps align calculated EDR with experimental toughness values, which will be explored in subsequent work.

For the purpose of comparison, the crack growth is simulated by ignoring DS in Fig. 7, where the initial load u0 is adjusted to 65 nm and the load increment remains unchanged. As shown in Fig. 7a, the predicted EDR is nearly identical to the effective intrinsic toughness. This is reasonable as the crack propagates in a linearly elastic material when DS is ignored. For the crack growth stage considered, the crack grows nearly straight forward. Only slight oscillations can be observed in the crack growth path, which is consistent with the phenomenon observed in the monotonic loading case (Fig. 4b).

Figure 7: The EDR and crack growth path for crack growth in polycrystal I when ignoring DS: (a) The EDR; (b) The crack growth path.

Comparing Fig. 7 with Fig. 6, it is observed that the crack propagates a greater distance under the same displacement load when domain evolution is ignored. This means the material has become more compliant by considering domain reorientation. This phenomenon represents another manifestation of the toughening effect induced by DS.

3.2 Crack Growth in Polycrystal II

In this section, the crack growth in polydomain t′-YSZ (polycrystal II) is studied. The initial polycrystal structure and the orientations of the grains are the same as those described in Fig. 2. For polydomain t′-YSZ, the initial polydomain structure results from the c-t′ phase transformation. The c-t′ transformation is modeled via the aforementioned phase-field approach, starting from the cubic phase (η=0). When simulating the c-t′ phase transformation, the evolution of the crack phase field is turned off by setting M=0 and the energy barrier is temporarily increased to 6.8×106J/m3 to facilitate the completion of the transformation. In atmospheric plasma-sprayed (APS) coatings, the c-t′ transformation occurs during rapid cooling [1]. Energy barrier temporary growth corresponds to molten-phase deposition. Random perturbation serves as phase transformation trigger. The remaining parameters are presented in Table 1.

The c-t′ phase transformation process is presented in Fig. 8. It can be observed that the material starts from the perturbed cubic phase. The cubic phase transforms to the t′ phase. The phase transformation leads to a twinned tetragonal microstructure, which is beneficial to accommodate the shear transformation strain. Recently, Petrich et al. [68] proposed a stochastic model to describe the herringbone microstructures in tetragonal zirconia, which form when an initial single cubic crystal cools from high temperatures. The influence of crystallographic domain morphology on the fracture behavior will be explored in our follow-up research.

Figure 8: The c-t′ phase transformation process: (a) t=0µs; (b) t=40µs; (c) t=60µs; (d) t=80µs; (e) t=100µs; (f) t=200µs.

After the c-t′ phase transformation, the value of Gb is adjusted back to 3.0×106J/m3 to maintain consistency with that of the mono domain case. The mobility parameter is set to M=0.5m/J/s to activate the evolution of the crack phase field. The mobility parameter for the order parameter η keeps unchanged. Firstly, the crack growth under a monotonically increasing displacement load is considered. The displacement load increases with the rate 87.5nm/µs until u∗=70nm. After that, u∗ increases with a slower rate (5nm/µs) to improve the numerical stability. The crack growth path and the evolution of the domain structure are shown in Fig. 9.

Figure 9: Crack growth in polycrystal II under monotonic loading: (a) t=0µs,u∗=0nm; (b) t=2.8µs,u∗=80nm; (c) t=4.8µs,u∗=90nm; (d) t=6.8µs,u∗=100nm; (e) t=8.8µs,u∗=110nm; (f) t=9.8µs,u∗=115nm.

Fig. 9 shows that as the load increases, the domain boundaries start to migrate, resulting in domain merging. As the crack propagates forward, most of the domain boundaries diminish in the crack tip region. Back-switching is only observed in some local regions in the wake of the crack. The back-switching is due to the change of the stress state with the advancement of the crack tip. However, back-switching is not observed in most of the grains left behind the crack tip. In Fig. 9f, the grains around the crack have been predominantly occupied by mono domains. A comparison of Fig. 9a,f clearly shows that massive DS has occurred in the polycrystal. In the meantime, DS perturbs the stress field around the crack tip, leading to crack deflection and resulting in a tortuous crack growth path. Fig. 10 shows a complex interaction mechanism between the crack growth and DS.

Figure 10: Crack growth in polycrystal II under monotonic loading with DS ignored: (a) t=0µs,u∗=0nm; (b) t=3µs,u∗=53nm; (c) t=8µs,u∗=78nm; (d) t=12.9µs,u∗=102.5nm.

A comparative study on crack growth in polycrystal II with ignoring DS is presented in Fig. 10. Since the domain evolution has been switched off, the initial twinned domain morphology stays unchanged with crack advances. Nevertheless, internal stress within twin domains can locally alter crack propagation trajectory, as shown in Fig. 10b,d. A comparison of Figs. 9b and 10b reveals that the initiation displacement load of crack growth is significantly larger when considering DS. A similar trend is observed during the crack growth stage. This implies that the material is more compliant with considering DS. Notably, analogous conclusions have been drawn for crack growth in polycrystal I. Fig. 11 compares crack propagation trajectories. Tortuous paths appear regardless of domain switching consideration, indicating that domain switching and twin domain internal stress both induce crack deflection. However, the crack growth paths for these two cases are obviously different.

Figure 11: A comparison of the crack growth paths in Polycrystal II: (a) Considering DS; (b) Ignoring DS.

Next, crack growth under incremental loading is simulated for the purpose of calculating the EDR. Again, a comparative study is conducted by ignoring DS. For the case with consideration of DS, the initial load u0 and the incremental load Δu are 102 and 0.5 nm, respectively. The initial load is applied at a rate of 34nm/µs and then held constant for 2µs. After that, the incremental load is applied sequentially at a rate of 1nm/µs and then held constant for 2µs at each step. The crack growth and domain structure evolution are similar to those shown in Fig. 9. For the case with ignorance of DS, the initial load is applied at a rate of 61.5nm/µs until u0=61.5nm and then held constant for 1µs. The incremental load Δu=0.5nm is applied at a rate of 1nm/µs and held constant for 0.5µs at each load step. The crack growth and domain structure evolution are similar to those shown in Fig. 10. To reduce computational costs arising from the high expense of the incremental loading case, only a limited crack growth length has been simulated for both scenarios. The EDRs are calculated using Eq. (26) and plotted against the crack growth length in Fig. 12.

Figure 12: The EDR for crack growth in polycrystal II.

Fig. 12 indicates that with DS taken into account, the predicted EDR is markedly greater than the effective inherent fracture toughness, demonstrating that DS produces local toughening during crack propagation. By comparison, when DS is excluded, EDR stays marginally higher than the inherent fracture toughness throughout most crack growth processes. Given that the EDR values are obtained via numerical calculation at separate load steps and are susceptible to numerical deviations, the internal stress field of the original twin domain structure exerts barely noticeable impacts on local crack growth resistance for the studied crack growth stage.

Fig. 12 also shows that the crack growth rate changes abruptly at specific points. For the load step starting from point B1, the crack growth length is distinctly larger, and the EDR increases dramatically. For the load steps between point B2 and point B3, the crack growth lengths are very small, and the EDR fluctuates dramatically. The crack topology and the domain pattern for the marked points are shown in Fig. 13. As illustrated in Fig. 13a,b, for the load step between point B1 and point B2, the crack exhibits unstable growth. The crack growth length is distinctly larger compared with the other load steps, although the load increment is tiny. The domain pattern corresponding to B2 also differs greatly from that corresponding to B1. For the crack growth stage between B2 (corresponding to Fig. 13b) and B3 (corresponding to Fig. 13c), the crack tip is approaching the grain boundary. The orientations of the two neighboring grains (−2° and 33°) differ substantially. The obstruction of the grain boundary to the crack growth may account for the low crack growth rate and the resulting fluctuations in the EDR. The increase of EDR at point C1 is attributed to a similar mechanism, i.e., the complex stress state around the grain boundary has retarded crack growth (Fig. 13d).

Figure 13: Domain patterns at typical load steps: (a) cl=961.5nm; (b) cl=1367.8nm; (c) cl=1375.9nm; (d) cl=1385.0nm. (a–d) Correspond to points B1, B2, B3, and C1 marked in Fig. 12, respectively.

The comparative studies presented in Sections 3.1 and 3.2 yield critical insights. First, it can be observed that the intricate initial multi-domain structure in Polycrystal II exerts merely trivial influences on local crack propagation resistance when its evolution is restrained. Nevertheless, as shown in Fig. 10b,d, the internal stress field of the multidomain structure can still trigger local deflection of the crack propagation path. It should be mentioned that the anisotropy of the elastic constants and the fracture toughness was not considered in this study. Incorporating these factors is anticipated to enhance the distinguishability of the initial twinned poly-domain pattern’s influence on crack growth paths. Secondly, the domain evolution during crack growth exhibits distinct characteristics in Polycrystal I and Polycrystal II. In mono-domain t′-YSZ, intense stress concentration at the crack tip induces domain reorientation in the tip vicinity. The concurrent evolution of domain structure and crack propagation forms a tortuous crack path. Due to the complex interplay between crack propagation and domain structure, the crack growth rate fluctuates sharply. Domain reorientation induces a local toughening effect on crack growth, which can be effectively characterized using the concept of energy dissipation rate. In poly-domain t′-YSZ, domain evolution primarily occurs via domain boundary migration. As the crack tip advances, domain merging ensues. Though the internal stress field of the original multi-domain structure imposes little effect on the local crack propagation resistance, domain boundary migration yields remarkable local toughening for crack extension. Thirdly, local toughening closely relates to domain evolution near the crack tip. Domain reverse switching and the presence of grain boundaries markedly affect the crack growth rate and fracture resistance. For both polycrystal structures investigated in this study, domain switching improves material compliance.

The findings presented in this study provide critical insights into the fracture and ferroelastic toughening mechanisms of polycrystalline thermal barrier coatings (TBCs) with diverse initial microstructures. Moreover, the phase-field model developed in this work and the proposed methodology for characterizing crack growth resistance through the energy dissipation rate concept establish a robust framework for evaluating the fracture properties of TBCs at the mesoscale. In our future work, the actual microstructural morphology of TBCs will be reconstructed based on electron backscatter diffraction(EBSD) characterization. The proposed PF model will be further extended to evaluate the fracture toughness of real TBC specimens. In addition, plastic slip is not incorporated in the current simulations. Under high stress concentration, local plastic deformation tends to occur. The interaction between plastic slip and grain boundaries can change crack propagation paths and induce additional toughening effects. Such mechanisms will be integrated into our subsequent research.

4  Conclusions

This work systematically investigates the fracture behaviors of single-domain and multi-domain t′-YSZ polycrystals via the phase-field approach. A unified phase-field model is established to describe the coupled evolution of domain structure and crack propagation. The numerical results show that the coupled evolution of crack phase and domain structure leads to tortuous crack growth paths in both polycrystals. Besides inducing crack deflection, crack-tip stress-driven domain reorientation enhances crack growth resistance in both materials, despite distinct domain evolution characteristics. Comparative analyses with domain switching neglected case clearly demonstrate that the toughening effect originates from domain reorientation, instead of the internal stress field of the initial domain structure. As the crack tip advances, back-switching of the domain structure emerges locally in the crack wake zone. The back-switching of the domain structure and the grain boundaries greatly affect the crack growth rate and the fracture resistance. These findings provide critical insights into the fracture toughening mechanisms of polycrystalline t′-YSZ. It is worth noting that this study only focuses on 2D simulations. In future research, the proposed PF model will be extended to investigate three-dimensional crack propagation, and more realistic microstructures of thermal barrier coatings will be taken into account.

Acknowledgement: The numerical computations in this work were performed on the High-Performance Computing (HPC) Platform of Huazhong University of Science and Technology. The authors would like to thank Ruobing Li for her kind assistance in formatting the manuscript.

Funding Statement: This research was funded by the National Natural Science Foundation of China, grant number 11872023.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Jun Luo; methodology, Zhou Fang, Jiaqi Zhong; formal analysis, Zhou Fang, Jiaqi Zhong, Jun Luo; investigation, Zhou Fang, Jiaqi Zhong, Jun Luo, Yuanzun Sun; data curation, Zhou Fang, Jiaqi Zhong; writing—original draft preparation, Zhou Fang; writing—review and editing, Jun Luo; project administration, Jun Luo; funding acquisition, Jun Luo. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data available on request from the authors.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Darolia R. Thermal barrier coatings technology: critical review, progress update, remaining challenges and prospects. Int Mater Rev. 2013;58(6):315–48. doi:10.1179/1743280413Y.0000000019. [Google Scholar] [CrossRef]

2. Vaßen R, Jarligo MO, Steinke T, Mack DE, Stöver D. Overview on advanced thermal barrier coatings. Surf Coat Technol. 2010;205(4):938–42. doi:10.1016/j.surfcoat.2010.08.151. [Google Scholar] [CrossRef]

3. Padture NP. Advanced structural ceramics in aerospace propulsion. Nat Mater. 2016;15(8):804–9. doi:10.1038/nmat4687. [Google Scholar] [PubMed] [CrossRef]

4. Evans AG, Clarke DR, Levi CG. The influence of oxides on the performance of advanced gas turbines. J Eur Ceram Soc. 2008;28(7):1405–19. doi:10.1016/j.jeurceramsoc.2007.12.023. [Google Scholar] [CrossRef]

5. Smialek JL, Miller RA. Revisiting the birth of 7YSZ thermal barrier coatings: stephan stecura. Coatings. 2018;8(7):255. doi:10.3390/coatings8070255. [Google Scholar] [CrossRef]

6. Kumar V, Balasubramanian K. Progress update on failure mechanisms of advanced thermal barrier coatings: a review. Prog Org Coat. 2016;90(93):54–82. doi:10.1016/j.porgcoat.2015.09.019. [Google Scholar] [CrossRef]

7. Evans AG, Mumm DR, Hutchinson JW, Meier GH, Pettit FS. Mechanisms controlling the durability of thermal barrier coatings. Prog Mater Sci. 2001;46(5):505–53. doi:10.1016/S0079-6425(00)00020-7. [Google Scholar] [CrossRef]

8. Dwivedi G, Viswanathan V, Sampath S, Shyam A, Lara-Curzio E. Fracture toughness of plasma-sprayed thermal barrier ceramics: influence of processing, microstructure, and thermal aging. J Am Ceram Soc. 2014;97(9):2736–44. doi:10.1111/jace.13021. [Google Scholar] [CrossRef]

9. Hu Z, Lan H, Sun X, Wu Y, Sun Y, Zhang H, et al. Thermo-physical properties and ferroelastic toughening of non-transformable t′-yttria-stabilized Hafnia-zirconia solid solution. Ceram Int. 2024;50(9):15292–302. doi:10.1016/j.ceramint.2024.02.005. [Google Scholar] [CrossRef]

10. Luo L, Chen Y, Zhou M, Shan X, Lu J, Zhao X. Progress update on extending the durability of air plasma sprayed thermal barrier coatings. Ceram Int. 2022;48(13):18021–34. doi:10.1016/j.ceramint.2022.04.044. [Google Scholar] [CrossRef]

11. Loganathan A, Gandhi AS. Effect of phase transformations on the fracture toughness of t′ yttria stabilized zirconia. Mater Sci Eng A. 2012;556:927–35. doi:10.1016/j.msea.2012.07.095. [Google Scholar] [CrossRef]

12. Ren X, Pan W. Mechanical properties of high-temperature-degraded yttria-stabilized zirconia. Acta Mater. 2014;69:397–406. doi:10.1016/j.actamat.2014.01.017. [Google Scholar] [CrossRef]

13. Zhu J, Luo J. Study of transformation induced intergranular microcracking in tetragonal zirconia polycrystals with the phase field method. Mater Sci Eng A. 2017;701(5):69–84. doi:10.1016/j.msea.2017.06.060. [Google Scholar] [CrossRef]

14. Mercer C, Williams JR, Clarke DR, Evans AG. On a ferroelastic mechanism governing the toughness of metastable tetragonal-prime (t′) yttria-stabilized zirconia. Proc Math Phys Eng Sci. 2007;463(2081):1393–408. [Google Scholar]

15. Krogstad JA, Leckie RM, Krämer S, Cairney JM, Lipkin DM, Johnson CA, et al. Phase evolution upon aging of air plasma sprayed t′-zirconia coatings: II-microstructure evolution. J Am Ceram Soc. 2013;96(1):299–307. doi:10.1111/j.1551-2916.2012.05460.x. [Google Scholar] [CrossRef]

16. Chan CJ, Lange FF, Rühle M, Jue JF, Virkar AV. Ferroelastic domain switching in tetragonal zirconia single crystals—microstructural aspects. J Am Ceram Soc. 1991;74(4):807–13. doi:10.1111/j.1151-2916.1991.tb06929.x. [Google Scholar] [CrossRef]

17. Srinivasan GV, Jue JF, Kuo SY, Virkar AV. Ferroelastic domain switching in polydomain tetragonal zirconia single crystals. J Am Ceram Soc. 1989;72(11):2098–103. doi:10.1111/j.1151-2916.1989.tb06038.x. [Google Scholar] [CrossRef]

18. Bolon AM, Gentleman MM. Raman spectroscopic observations of ferroelastic switching in ceria-stabilized zirconia. J Am Ceram Soc. 2011;94(12):4478–82. doi:10.1111/j.1551-2916.2011.04737.x. [Google Scholar] [CrossRef]

19. Li J, Yang L, Zhou Y, Li Z, Liu Z, Tsai M. Preparationofepitaxialyttria-stabilizedzirconiawithnon-quilibrium t′ phase and ferroelastic domain by pulsed laser deposition. Appl Surf Sci. 2021;537:147790. doi:10.1016/j.apsusc.2020.147790. [Google Scholar] [CrossRef]

20. Li J, Zhou Q, Yang L, Zhou Y, Zhao J, Huang J, et al. Ferroelastic deformation mechanism and mechanical properties of [001]-oriented YSZ film by indentation. J Alloys Compd. 2021;889(7):161557. doi:10.1016/j.jallcom.2021.161557. [Google Scholar] [CrossRef]

21. Masuda H, Morita K, Watanabe M, Hara T, Yoshida H, Ohmura T. Ferroelastic and plastic behaviors in pseudo-single crystal micropillars of nontransformable tetragonal zirconia. Acta Mater. 2021;203:116471. doi:10.1016/j.actamat.2020.11.013. [Google Scholar] [CrossRef]

22. Luo C, Li C, Cao K, Li J, Luo J, Zhang Q, et al. Ferroelastic domain identification and toughening mechanism for yttrium tantalate-zirconium oxide. J Mater Sci Technol. 2022;127:78–88. doi:10.1016/j.jmst.2022.04.007. [Google Scholar] [CrossRef]

23. Carbogno C, Levi CG, Van de Walle CG, Scheffler M. Ferroelastic switching of doped zirconia: modeling and understanding from first principles. Phys Rev B. 2014;90(14):144109. doi:10.1103/physrevb.90.144109. [Google Scholar] [CrossRef]

24. Chen MH, Thomas JC, Natarajan AR, Van der Ven A. Effects of strain on the stability of tetragonal ZrO2. Phys Rev B. 2016;94(5):054108. doi:10.1103/physrevb.94.054108. [Google Scholar] [CrossRef]

25. Baither D, Bartsch M, Baufeld B, Tikhonovsky A, Foitzik A, Rühle M, et al. Ferroelastic and plastic deformation oft′-zirconia single crystals. J Am Ceram Soc. 2001;84(8):1755–62. doi:10.1111/j.1151-2916.2001.tb00911.x. [Google Scholar] [CrossRef]

26. Chen LQ. Phase-field models for microstructure evolution. Annu Rev Mater Res. 2002;32(1):113–40. doi:10.1146/annurev.matsci.32.112001.132041. [Google Scholar] [CrossRef]

27. Mamivand M, Asle Zaeem M, El Kadiri H. Phase field modeling of stress-induced tetragonal-to-monoclinic transformation in zirconia and its effect on transformation toughening. Acta Mater. 2014;64:208–19. doi:10.1016/j.actamat.2013.10.031. [Google Scholar] [CrossRef]

28. Zhu J, Luo J, Sun Y. Study of the fracture behavior of tetragonal zirconia polycrystal with a modified phase field model. Materials. 2020;13(19):4430. doi:10.3390/ma13194430. [Google Scholar] [PubMed] [CrossRef]

29. Fang B, Luo J. New insights into the thermal aging mechanism of yttria stabilized zirconia: a phase field study. Ceram Int. 2022;48(16):23365–80. doi:10.1016/j.ceramint.2022.04.328. [Google Scholar] [CrossRef]

30. Fang D, Li F, Liu B, Zhang Y, Hong J, Guo X. Advances in developing electromechanically coupled computational methods for piezoelectrics/ferroelectrics at multiscale. Appl Mech Rev. 2013;65(6):060802. doi:10.1115/1.4025633. [Google Scholar] [CrossRef]

31. Wang J, Zhang TY. Phase field simulations of polarization switching-induced toughening in ferroelectric ceramics. Acta Mater. 2007;55(7):2465–77. doi:10.1016/j.actamat.2006.11.041. [Google Scholar] [CrossRef]

32. Song YC, Soh AK, Ni Y. Phase field simulation of crack tip domain switching in ferroelectrics. J Phys D Appl Phys. 2007;40(4):1175. doi:10.1088/0022-3727/40/4/040. [Google Scholar] [CrossRef]

33. Abdollahi A, Arias I. Phase-field modeling of the coupled microstructure and fracture evolution in ferroelectric single crystals. Acta Mater. 2011;59(12):4733–46. doi:10.1016/j.actamat.2011.03.030. [Google Scholar] [CrossRef]

34. Abdollahi A, Arias I. Numerical simulation of intergranular and transgranular crack propagation in ferroelectric polycrystals. Int J Fract. 2012;174(1):3–15. doi:10.1007/s10704-011-9664-0. [Google Scholar] [CrossRef]

35. Xu BX, Schrade D, Gross D, Mueller R. Fracture simulation of ferroelectrics based on the phase field continuum and a damage variable. Int J Fract. 2010;166(1):163–72. doi:10.1007/s10704-010-9520-7. [Google Scholar] [CrossRef]

36. Sluka T, Webber KG, Colla E, Damjanovic D. Phase field simulations of ferroelastic toughening: the influence of phase boundaries and domain structures. Acta Mater. 2012;60(13–14):5172–81. doi:10.1016/j.actamat.2012.06.023. [Google Scholar] [CrossRef]

37. Pi ZP, Zhang F, Chen JB, Zhu W, Yang L, Zhou YC. Multiphase field theory for ferroelastic domain switching with an application to tetragonal zirconia. Comput Mater Sci. 2019;170:109165. doi:10.1016/j.commatsci.2019.109165. [Google Scholar] [CrossRef]

38. Pi ZP, Wang KL, Yang L, Zhou YC. On the theoretical and phase field modeling of the stress state associated with ferroelastic twin nucleation and propagation near crack tip. Eng Fract Mech. 2020;235:107200. doi:10.1016/j.engfracmech.2020.107200. [Google Scholar] [CrossRef]

39. Bhattacharya A, Asle Zaeem M. Mechanism of nucleation in ferroelastic domain switching. Scr Mater. 2024;252(1):116273. doi:10.1016/j.scriptamat.2024.116273. [Google Scholar] [CrossRef]

40. Bhattacharya A, Asle Zaeem M. Kinetics of ferroelastic domain switching with and without back-switching events: a phase-field study. Acta Mater. 2025;286(1):120702. doi:10.1016/j.actamat.2024.120702. [Google Scholar] [CrossRef]

41. Zhou Q, Wei Y, Zhou Y, Yang L. Phase field modeling of ferroelastic variant switching in yttria-stabilized t′zirconia with strain gradient elasticity and interface tension. Sci China Technol Sci. 2024;67(5):1443–57. doi:10.1007/s11431-022-2486-x. [Google Scholar] [CrossRef]

42. Li C, Li J. Phase field study of ferroelastic domain switching behavior related to grain size in polycrystalline tetragonal zirconia. J Am Ceram Soc. 2025;108(7):e20495. doi:10.1111/jace.20495. [Google Scholar] [CrossRef]

43. Zhao T, Zhu J, Luo J. Study of crack propagation behavior in single crystalline tetragonal zirconia with the phase field method. Eng Fract Mech. 2016;159(9):155–73. doi:10.1016/j.engfracmech.2016.03.035. [Google Scholar] [CrossRef]

44. Jafarzadeh H, Levitas VI, Farrahi GH, Javanbakht M. Phase field approach for nanoscale interactions between crack propagation and phase transformation. Nanoscale. 2019;11(46):22243–7. doi:10.1039/c9nr05960a. [Google Scholar] [PubMed] [CrossRef]

45. Moshkelgosha E, Mamivand M. Phase field modeling of crack propagation in shape memory ceramics—application to zirconia. Comput Mater Sci. 2020;174:109509. doi:10.1016/j.commatsci.2019.109509. [Google Scholar] [CrossRef]

46. Moshkelgosha E, Mamivand M. Concurrent modeling of martensitic transformation and crack growth in polycrystalline shape memory ceramics. Eng Fract Mech. 2021;241:107403. doi:10.1016/j.engfracmech.2020.107403. [Google Scholar] [CrossRef]

47. Borzabadi Farahani E, Sobhani Aragh B, Voges J, Juhre D. On the crack onset and growth in martensitic micro-structures; a phase-field approach. Int J Mech Sci. 2021;194(13):106187. doi:10.1016/j.ijmecsci.2020.106187. [Google Scholar] [CrossRef]

48. Lotfolahpour A, Huber W, Asle Zaeem M. A phase-field model for interactive evolution of phase transformation and cracking in superelastic shape memory ceramics. Comput Mater Sci. 2023;216(15):111844. doi:10.1016/j.commatsci.2022.111844. [Google Scholar] [CrossRef]

49. Lotfolahpour A, Asle Zaeem M. Orientation-dependent deformation and failure of micropillar shape memory ceramics: a 3D phase-field study. Extreme Mech Lett. 2024;73:102245. doi:10.1016/j.eml.2024.102245. [Google Scholar] [CrossRef]

50. Sun Y, Luo J, Zhu J. Ferroelastic toughening of single crystalline yttria-stabilized t′ zirconia: a phase field study. Eng Fract Mech. 2020;233:107077. doi:10.1016/j.engfracmech.2020.107077. [Google Scholar] [CrossRef]

51. Xiong J, Xu B, Kang G. Phase field simulations on the rate- and grain-size-dependent crack propagation of polycrystalline NiTi shape memory alloy. Fatigue Fract Eng Mater Struct. 2024;47(6):2174–94. doi:10.1111/ffe.14296. [Google Scholar] [CrossRef]

52. Xiong J, Xu B, Hu J, Kang G. Phase-field simulation on grain-size dependent fracture of cyclically loaded NiTi-SMA. Int J Mech Sci. 2025;289:110041. doi:10.1016/j.ijmecsci.2025.110041. [Google Scholar] [CrossRef]

53. Fang Z, Luo J, Sun Y, Huang H, Zhong J. Phase field study of crack growth in t′ yttria stabilized zirconia with initial domain structures. Mater Today Commun. 2023;37(6):107489. doi:10.1016/j.mtcomm.2023.107489. [Google Scholar] [CrossRef]

54. Bourdin B, Francfort GA, Marigo JJ. Numerical experiments in revisited brittle fracture. J Mech Phys Solids. 2000;48(4):797–826. doi:10.1016/S0022-5096(99)00028-9. [Google Scholar] [CrossRef]

55. Bourdin B, Francfort GA, Marigo JJ. The variational approach to fracture. J Elast. 2008;91(1):5–148. doi:10.1007/s10659-007-9107-3. [Google Scholar] [CrossRef]

56. Ter Haar D. Collected papers of L. D. Landau. Oxford, UK: Pergamon; 1965. [Google Scholar]

57. Sun Y, Luo J, Zhu J, Zhou K. A non-isothermal phase field study of the shape memory effect and pseudoelasticity of polycrystalline shape memory alloys. Comput Mater Sci. 2019;167:65–76. doi:10.1016/j.commatsci.2019.05.036. [Google Scholar] [CrossRef]

58. Rajendran MK, Kuna M, Budnitzki M. Undercooling versus stress induced martensitic phase transformation: the case of MgO-partially stabilized zirconia. Comput Mater Sci. 2020;174(3):109460. doi:10.1016/j.commatsci.2019.109460. [Google Scholar] [CrossRef]

59. Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits. Comput Meth Appl Mech Eng. 2010;199(45–48):2765–78. doi:10.1016/j.cma.2010.04.011. [Google Scholar] [CrossRef]

60. Ambati M, Gerasimov T, De Lorenzis L. A review on phase-field models of brittle fracture and a new fast hybrid formulation. Comput Mech. 2015;55(2):383–405. doi:10.1007/s00466-014-1109-y. [Google Scholar] [CrossRef]

61. Wu JY, Nguyen VP. A length scale insensitive phase-field damage model for brittle fracture. J Mech Phys Solids. 2018;119(8):20–42. doi:10.1016/j.jmps.2018.06.006. [Google Scholar] [CrossRef]

62. Amor H, Marigo JJ, Maurini C. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. J Mech Phys Solids. 2009;57(8):1209–29. doi:10.1016/j.jmps.2009.04.011. [Google Scholar] [CrossRef]

63. Zhou S, Zhuang X, Rabczuk T. Phase field modeling of brittle compressive-shear fractures in rock-like materials: a new driving force and a hybrid formulation. Comput Meth Appl Mech Eng. 2019;355(2):729–52. doi:10.1016/j.cma.2019.06.021. [Google Scholar] [CrossRef]

64. Awaji H, Ebisudani M, Choi SM, Ohashi T. Crack deflection toughening mechanism in brittle materials. In: Fracture resistance testing of monolithic and composite brittle materials. West Conshohocken, PA, USA: ASTM; 2002. p. 137–51. doi:10.1520/stp10476s. [Google Scholar] [CrossRef]

65. Zhao H, Li Z, Gao H, Lu L. Fracture and toughening mechanisms in nanotwinned and nanolayered materials. MRS Bull. 2022;47(8):839–47. doi:10.1557/s43577-022-00376-5. [Google Scholar] [CrossRef]

66. Nairn JA. Simulation of crack growth in ductile materials. Eng Fract Mech. 2005;72(6):961–79. doi:10.1016/j.engfracmech.2004.08.006. [Google Scholar] [CrossRef]

67. Xu W, Ren Y, Xiao S, Liu B. A finite crack growth energy release rate for elastic-plastic fracture. J Mech Phys Solids. 2023;181(3):105447. doi:10.1016/j.jmps.2023.105447. [Google Scholar] [CrossRef]

68. Petrich L, Derrien K, Schmidt V, Guinebretière R, Moulinec H, Castelnau O. Thermo-elastic micromechanical modeling of tetragonal ZrO2 with a herringbone microstructure inherited from the cubic phase. Mater Des. 2025;257:114425. doi:10.1016/j.matdes.2025.114425. [Google Scholar] [CrossRef]