Phase-field simulation of thermal cracking with length-scale insensitive degradation functions

Abstract

The present paper applies a length-scale insensitive degradation function to accelerate the simulation of thermal cracking problems. The solid deformation, damage evolution, and heat conduction process are coupled. The solid is treated as a linear elasticity, whose strain is composed of the mechanical strain and the thermal strain induced by thermal expansion. The cracks in solids are represented by the diffusive phase-field variables. The length-scale insensitive degradation function has been employed to decouple the length scale of the phase-field length scale and that of the physical process area, which alleviates the meshing burden of the phase field. The numerical results of our ceramic examples of interest are in good agreement with the experimental results, and this work may have important implications for predicting thermal cracking problems on large scales.

Introduction

The drastic temperature changes can induce thermal shock stress in solids. When the thermal stress value exceeds the strength of the solids, cracks will generate and propagate. ¹ Many efforts have been made to simulate the thermal cracking phenomenon, for example, extended fine element methods ² and remeshing strategies, ³ but these methods are cumbersome in simulating complex fracture networks. ⁴ In comparison, the phase-field method modeled the sharp cracks in a diffusive way. The phase-field method originated in the 1950s and was initially used to address mathematical modeling problems in phase transition processes such as solid-phase transformations and droplet coalescence. The Landau–Ginzburg theory was introduced to describe the phase field in continuous phase transition phenomena. Later, in 1969, Cahn and Hilliard first proposed the phase-field equation, which became an important foundation for phase-field method research. An additional field variable (phase-field variable) varying between 0 and 1 continuously is introduced to represent the cracks, with 0 meaning that the solid is intact and 1 indicating that the solid generates cracks. The advantage of the phase-field method is that it unifies the nucleation, propagation, branching, and coalescence of cracks.

Since its inception, phase-field method has been widely used in computational fracture mechanics, including brittle fracture, dynamic fracture,^5–7 ductile fracture, ⁸ hydraulic fracture,^9–11 fracture in piezoelectric ceramics, ¹² and composite electrode. ³¹ The phase-field fracture model ³² seeks the displacement field and crack set by minimizing the total potential energy of the cracked solid, in order to extend Griffith's theory. Some scholars conducted thermal shock experiments on rectangular ceramic plates using the phase-field cohesive zone model, ¹³ others have implemented a parallel algorithm for the thermoelastic coupled phase-field model and applied it to simulate dynamic brittle fracture and quasistatic brittle fracture of thermoelastic materials. ¹⁴ Thermodynamically consistent phase-field model for thermal fracture is proposed. ¹⁵

Nonetheless, within the phase-field method, the region surrounding a fracture is intricately connected to a minor physical length scale. This scale is associated with the fracture energy and the strength of the material. Achieving precise outcomes necessitates the utilization of an exceedingly fine mesh to accurately capture this length scale. The huge mesh density reduces the efficiency of phase-field method. In order to improve the computational efficiency of the phase-field method, many scholars both domestically and internationally have made significant efforts. For instance, Zhang et al. ³⁴ demonstrated the effectiveness of certain acceleration techniques in dealing with dual-crack systems. Reduced-dimensional models ³⁵ are also an option, which simplifies the problem by reducing the degrees of freedom. Additionally, adaptive mesh refinement methods, such as the one proposed by Goswami et al., ³⁶ can dynamically adjust the mesh density based on the solution characteristics to balance the computational accuracy and efficiency. Yin et al. ³³ propose a novel phase-field fracture method to make the relation between the diffusive length scale and the material parameters adjustable via modifying the energetic degradation functions. Furthermore, some researchers proposed a length-scale insensitive degradation function for phase-field models, and it is widely used to solve the fracture problem of large linear elastic structures. ¹⁶ The novel degradation functions decouple the physical length scale and the phase-field length scale, thus allowing for a phase-field length scale much larger than the physical process zone length scale.

The present paper is intended to apply the length-scale insensitive degradation function to accelerate the simulation of thermal cracking problems, and it is an initial validation of the potential ability of the method for solving large-scale thermal cracking models. The quenching tests of the ceramic plate are considered to validate the feasibility and numerical implementation. The simulated crack morphology and response are compared with the experimental results.

The remainder of the paper is organized as follows. The accelerated phase-field models for thermal cracking are introduced in “Theoretical equation” section. “Finite element discretization” section is the finite element implementation of the governing equations. “Algorithmic examples” section shows some typical examples of thermal shocks involving brittle fracture failure, followed by “Conclusion” section.

Theoretical equation

The theoretical aspect of the thermal fracture phase-field method is discussed, focusing on achieving constant peak stress by modifying the degradation function, regardless of the length scale of the physical process zone surrounding the fracture.

Displacement and damage equations

The solid crack region is represented by $\mathrm{\Omega }\subset {\mathrm{R}}^{n_{\dim }}$ ( $n_{\dim }=1,2,3$ ) as shown in Figure 1, while the solid boundary is characterized by $\mathrm{\partial }\mathrm{\Omega }\subset {\mathrm{R}}^{n_{\dim \text{-}1}}$ , the normal outward vector is denoted by $\mathrm{n}$ . The body force per unit is designated by $\mathrm{b}*$ , we prescribe traction boundary condition on $\mathrm{\partial }{\mathrm{\Omega }}_t\subset \mathrm{\partial }\mathrm{\Omega }$ , and displacement boundary condition on $\mathrm{\partial }{\mathrm{\Omega }}_u\subset \mathrm{\partial }\mathrm{\Omega }\mathrm{\setminus }\mathrm{\partial }{\mathrm{\Omega }}_t$ . Sharp cracks in the solid are represented in the domain $\mathrm{S}\subset {\mathrm{R}}^{n_{\dim \text{-}1}}$ , with its boundary $\mathrm{\partial }\mathrm{S}$ and outward normal vector $n_s$ .

Figure 1.

A cracking solid.

In the phase-field method, ¹⁷ sharp fracture $\mathrm{S}$ is in the fracture zone $\mathrm{B}(\mathrm{B}\in \mathrm{\Omega })$ , and the bandwidth of fracture zone $\mathrm{B}$ is measured by the phase-field length scale $b>0$ . If b approaches 0, the sharp fracture $\mathrm{S}$ will recover, and the fracture zone $\mathrm{B}$ is also bounded. The outer boundary of fracture zone $\mathrm{B}$ is denoted as $\mathrm{\partial }\mathrm{B}$ , with the outward normal vector denoted by ${\mathrm{n}}_{\mathrm{B}}$ . The fracture zone $\mathrm{B}$ is continuously updated through the evolution of the fracture phase field $d(x)$ . $d(x)=0$ indicates that the material is intact, whereas $d(x)=1$ indicates the presence of macro-cracks. Once cracks initiated, they will be irreversible, so we specify $d(x)\ge 0$ .

Governing equations

The mechanical control equation and the damage control equation can be written as ²¹ :

$$\begin{array}{c}\{\begin{array}{c}\mathrm{\nabla }\cdot \sigma +b=0\mathrm{in}\mathrm{\Omega }\\ \sigma \cdot n=t*\mathrm{on}\mathrm{\partial }{\mathrm{\Omega }}_t\end{array}\end{array}$$

(1)

$$\begin{array}{c}\{\begin{array}{c}\mathrm{\nabla }\cdot q+Q\le 0\mathrm{in}B\\ q\cdot n_{\mathrm{B}}\ge 0\mathrm{on}\mathrm{\partial }B\end{array}\end{array}$$

(2)

where q represents a flux vector that is work-conjugate to $\mathrm{\nabla }d$ , while Q denotes its corresponding source term. $t*$ represents the force boundary condition applied on $\mathrm{\partial }{\mathrm{\Omega }}_t$ and we will discuss stress tensor $\sigma$ in the following chapters. Besides, b stands for the body force applied, which is expressed by:

$$\begin{array}{c}q={\displaystyle \frac{2b}{c_{\alpha }}{G}_f\mathrm{\nabla }d,Q=-{\omega }^{\mathrm{\prime }}(d)\overline{Y}-{\alpha }^{\mathrm{\prime }}(d){\displaystyle \frac{G_f}{c_{\alpha }b}}}\end{array}$$

(3)

where $G_f$ represents the fracture energy, $\overline{Y}$ denotes the effective crack driving force, $c_{\alpha }$ is a scaling parameter constant, $c_{\alpha }=4\underset{0}{\overset{1}{\int }}\sqrt{\alpha (\beta )}d\beta$ , and these parameters are integral to the characterization of material failure mechanisms within the model. $\alpha (d)$ represents the geometric crack function, which increases with the increase of d, and $\omega (d)$ denotes the energy degradation function, which decreases with the increase of d. In order to inhibit crack healing, the enforcement of a crack irreversibility condition necessitates substituting a variable, denoted as a history variable $\mathcal{H}$ , ¹⁸ representing the maximum value $\overline{Y}$ it has ever attained. Consequently, the source term in equation (3) can be reformulated as:

$$\begin{array}{c}Q=-{\omega }^{\mathrm{\prime }}(d)H-{\alpha }^{\mathrm{\prime }}(d){\displaystyle \frac{G_f}{c_{\alpha }b},H=\max ({\overline{Y}}_0,{\max }_{n\in [0,T]}{\overline{Y}}_n)}\end{array}$$

(4)

where ${\overline{Y}}_n$ denotes the effective crack driving force at a given time $t_n$ . The parameter ${\overline{Y}}_0$ is indicative of a threshold level for $\overline{Y}$ , such that when the material's stress remains beneath this threshold, it maintains its linearly elastic behavior. After the above adjustment, equation (2) becomes equalities.

Constitutive relations

The constitutive equations that define the stress $\sigma$ and the effective crack driving force $\overline{Y}$ ¹⁹ are expressed as follows:

$$\begin{array}{c}\sigma =\omega (d)\overline{\sigma }=\omega (d)E_0:ϵ\end{array}$$

(5)

Considering only elastic strain, using Hooke's law to express the effective stress tensor $\overline{\sigma }$ :

$$\begin{array}{c}\overline{\sigma }=E_0:ϵ\end{array}$$

(6)

where $E_0$ is the fourth-order elasticity tensor,

$$\begin{array}{c}{E}_0={\lambda }_{0}I\otimes I+2{\mu }_{0}I\end{array}$$

(7)

in which ${\lambda }_0:=\frac{v_0{E}_0}{(1-2v_0)(1+v_0)}$ and ${\mu }_0:=\frac{E_0}{2(1+v_0)}$ are the Lame constants. In order to resolve the issue with equation (4), which does not distinguish between material compression and tension, we use the following formulation for the effective crack driving force:

$$\begin{array}{c}\overline{Y}={\displaystyle \frac{1}{2{\overline{E}}_0}{\overline{{\sigma }_1}}^2\to {\overline{Y}}_0={\displaystyle \frac{f_t^2}{2{\overline{E}}_0}}}\end{array}$$

(8)

Here, $\overline{{\sigma }_1}$ represents the maximum principal value within the effective stress tensor $\overline{\sigma }$ . The elongation modulus, denoted by ${\overline{E}}_0$ , characterizes the material's propensity to deform under stress, ${\overline{E}}_0={\lambda }_0+2{\mu }_0$ . Furthermore, $f_t$ signifies the material's failure strength, which corresponds to the peak stress at which the material fails, and they are thermodynamically consistent. ²⁰ Equation (8) indicates that cracks initiate once the $\overline{{\sigma }_1}$ reaches $f_t$ .

Thermal equations

In this section, we introduced the heat conduction process in the solid. ²¹ The external boundary is divided into two parts, one is given the temperature $T*$ , and the other is given the heat flux $h*$ .

Governing equations for the heat transfer process

The following equation describes the heat conduction process in the solid,

$$\begin{array}{c}\{\begin{array}{c}\rho c\dot{T}+\mathrm{\nabla }\cdot h=\mathrm{\Theta }*\mathrm{in}\mathrm{\Omega }\\ h\cdot n=h*\mathrm{on}\mathrm{\partial }{\mathrm{\Omega }}_h\end{array}\end{array}$$

(9)

where $\rho$ represents material density and c stands for the specific heat capacity. $\mathrm{\Theta }*$ represents the volume heat source. The heat flux is calculated by Fourier law:

$$\begin{array}{c}h=-\kappa \cdot \mathrm{\nabla }T=-\omega (d){\kappa }_0\cdot \mathrm{\nabla }T\end{array}$$

(10)

${\kappa }_0={\kappa }_{0}I$ is the second-order isotropic thermal conductivity tensor. Note that ${\kappa }_0$ is a constant here.

Constitutive relations

The strain in solids consists of two parts one is mechanical strain ${ϵ}^{\mathrm{m}}$ , and the other is thermal strain ${ϵ}^{\mathrm{t}}$ caused by thermal expansion:

$$\begin{array}{c}ϵ={ϵ}^{\mathrm{m}}+{ϵ}^{\mathrm{t}}\end{array}$$

(11)

$$\begin{array}{c}{ϵ}^{\mathrm{t}}={\eta }_t(T-T_0)I\end{array}$$

(12)

$T_0$ is the selected reference temperature, and ${\eta }_t$ represents the coefficient of thermal expansion. After the above decomposition, the effective stress tensor $\overline{\sigma }$ needs to be adjusted to the following form:

$$\begin{array}{c}\overline{\sigma }=E_0:{ϵ}^{\mathrm{m}}=E_0:(ϵ-{ϵ}^{\mathrm{t}})\end{array}$$

(13)

We assume that the mechanical properties of materials do not change during thermal shock.^22,23

The length-scale insensitive degradation function

The phase-field formulation for thermal cracking has already been provided. In traditional phase-field modeling, the geometric crack functions and energetic degradation functions typically assume the following forms:

$$\begin{array}{c}\{\begin{array}{c}\alpha (d)=d^2\\ \omega (d)={(1-d)}^2\end{array}\end{array}$$

(14)

The scaling parameters $c_{\alpha }=2$ are specific to AT2 ²⁴ model. For a one-dimensional bar experiencing uniform tension, the application of a conventional degradation function leads to the determination of the peak stress ${\sigma }_c$ , which signifies the failure strength, and ${\sigma }_c$ is expressed as follows:

$$\begin{array}{c}{\sigma }_c=\sqrt{{\displaystyle \frac{27G_f{\overline{E}}_0}{256{\ell }_0}}}\end{array}$$

(15)

${\ell }_0=27G_f{\overline{E}}_0/256{\sigma }_c^2$ is termed as the length scale of the phase field, obtainable through computation. Within the traditional phase-field model, there's no differentiation between the physical length scale and the phase-field length scale, resulting in a numerical approximation of their equality. However, applying this model to large-scale structures leads to computational intensity, rendering it challenging to perform calculations for large engineering structures. To decouple the phase-field length scale from the physical length scale, a new degradation function was introduced by Lo et al. ¹⁶ :

$$\begin{array}{c}\omega (d)=s\left[1-{\left({\displaystyle \frac{s-1}{s}}\right)}^{{(1-d)}^2}\right]\end{array}$$

(16)

where the parameter s is a dimensionless parameter that enables the decoupling of the phase-field length scale and the physical length scale, which means that the dependency between these two scales is separated through s as an intermediary, allowing them to vary independently. The benefit of this approach is that even when the model size increases, the grid density can remain unchanged, thus avoiding the significant increase in computational resources caused by grid refinement.

Regarding the graphical representation of this energy degradation function and the decoupling of phase-field length scale and physical length scale through parameter s, interested readers can refer to Yang et al. ²⁵

Employing this degradation function enables us to redefine the expression for the peak stress in a bar subjected to homogeneous uniaxial tension as follows:

$$\begin{array}{c}{\sigma }_c*={\displaystyle \frac{\sqrt{{\displaystyle \frac{27G_f{\overline{E}}_0}{256{\ell }_0}}}}{\sqrt{(s-1)\ln \left({\displaystyle \frac{s}{s-1}}\right)}}}\end{array}$$

(17)

${\sigma }_c*$ denotes the peak stress obtained through the application of the degradation function that is invariant with respect to the length scale. From this, we can establish the correlation between the physical length scale ${\ell }_p$ and the phase-field length scale ${\ell }_0$ :

$$\begin{array}{c}{\ell }_p=(s-1)\ln \left({\displaystyle \frac{s}{s-1}}\right){\ell }_0\end{array}$$

(18)

By replacing equation (17) with equation (18), ${\sigma }_c^{*}$ can be reformulated as ${\sigma }_c^{*}=\sqrt{(27G_f{\overline{E}}_0)/256{\ell }_p}$ . By assigning various values to the parameter s in the degradation function presented in equation (16), we can retain an equivalent peak stress as seen in the traditional form of the phase-field model while simultaneously augmenting the phase-field length scale. Based on the aforementioned analysis, it is evident that the degradation function, which is insensitive to the length scale, facilitates the separation of the physical and phase-field length scales.

Finite element discretization

This section presents the finite element implementation of the model using the multifield elements. Here, multifield elements consist of displacement degrees of freedom, damage freedom, and temperature degrees of freedom. “Weak form” section describes the weak formal expression of the multifield element under study, and “Finite element discretization” section is the finite element discrete part.

Weak form

Next, we need to rewrite the above control equations into the following weak form: $\mathrm{u}\in {\mathcal{U}}_u$ , $T\in {\mathcal{U}}_T$ and $d\in {\mathcal{U}}_d$ such that:

$$\begin{array}{c}\underset{\mathrm{\Omega }}{\int }\sigma :{\mathrm{\nabla }}^{\mathrm{sym}}\delta \mathbf{u}\mathrm{d}V=\delta {\mathcal{P}}^u\mathrm{\forall }\delta u\in {\mathcal{V}}_u\end{array}$$

(19)

$$\begin{array}{c}\underset{\mathrm{\Omega }}{\int }\rho c\dot{T}\delta T\mathrm{d}V-\underset{\mathrm{\Omega }}{\int }\mathrm{h}\cdot \mathrm{\nabla }(\delta T)\mathrm{d}V=\delta {\mathcal{P}}^T\mathrm{\forall }\delta T\in {\mathcal{V}}_T\end{array}$$

(20)

$$\begin{array}{c}\underset{\mathcal{B}}{\int }(-Q\delta d+\mathrm{q}\cdot \mathrm{\nabla }(\delta d))\mathrm{d}V=0\mathrm{\forall }\delta d\in {\mathcal{V}}_d\end{array}$$

(21)

The virtual power $\delta {\mathcal{P}}^u$ of external forces and the virtual power $\delta {\mathcal{P}}^T$ of heat sources are expressed by

$$\begin{array}{c}\delta {\mathcal{P}}^u=\underset{\mathrm{\Omega }}{\int }b*\cdot \delta \mathrm{u}\mathrm{d}V+\underset{\mathrm{\partial }{\mathrm{\Omega }}_t}{\int }t*\cdot \delta \mathrm{u}\mathrm{d}A\end{array}$$

(22)

$$\begin{array}{c}\delta {\mathcal{P}}^T=\underset{\mathrm{\Omega }}{\int }\mathrm{\Theta }*\delta T\mathrm{d}V-\underset{\mathrm{\partial }{\mathrm{\Omega }}_J}{\int }h*\delta T\mathrm{d}A\end{array}$$

(23)

Next, we introduce test and trial spaces:

$$\begin{array}{c}{\mathcal{U}}_u:=\{\mathrm{u}|\mathrm{u}(\mathrm{x})={\mathrm{u}}^{*}\mathrm{\forall }\mathrm{x}\in \mathrm{\partial }{\mathrm{\Omega }}_u\},{\mathcal{V}}_u:=\{\delta u\delta u(\mathrm{x})=0\mathrm{\forall }\mathrm{x}\in \mathrm{\partial }{\mathrm{\Omega }}_u\}\end{array}$$

(24)

$$\begin{array}{c}{\mathcal{U}}_T:=\{T|T(\mathrm{x})=T^{*}\mathrm{\forall }\mathrm{x}\in \mathrm{\partial }{\mathrm{\Omega }}_T\},{\mathcal{V}}_T:=\{\delta T\delta T(\mathrm{x})=0\mathrm{\forall }\mathrm{x}\in \mathrm{\partial }{\mathrm{\Omega }}_T\}\end{array}$$

(25)

$$\begin{array}{c}{\mathcal{U}}_d:=\{d|d(\mathrm{x})\in [0,1],\dot{d}(\mathrm{x})\ge 0\mathrm{\forall }\mathrm{x}\in \mathrm{B}\},{\mathcal{V}}_d:=\{\delta d\delta d(\mathrm{x})\ge 0\mathrm{\forall }\mathrm{x}\in \mathrm{B}\}\end{array}$$

(26)

Finite element discretization

The displacement field ${\mathrm{u}}^h$ , the temperature field $T^h$ and the damage field $d^h$ are expressed as:

$$\begin{array}{c}{u}^h(\mathrm{x})=\sum _A{N}_A^u(\mathrm{x})a_A=N_{u}a\end{array}$$

(27)

$$\begin{array}{c}{T}^h(\mathrm{x})=\sum _A{N}_A^T(\mathrm{x})T_A=N_T\hat{\mathit{a}}\end{array}$$

(28)

$$\begin{array}{c}{d}^h(\mathrm{x})=\sum _A{N}_A^d(\mathrm{x})d_A=N_d\overline{\mathit{a}}\end{array}$$

(29)

where $N_A^u(\mathrm{x})$ , $N_A^T(\mathrm{x})$ , and ${\mathrm{N}}_A^d(\mathrm{x})$ are the shape functions for the displacement, temperature, and damage field, respectively; and ${\mathrm{N}}_u$ , ${\mathrm{N}}_T$ , and ${\mathrm{N}}_d$ are the corresponding shape function matrices. We use the same shape function to apply it in all fields. The vectors $\mathit{a}$ , $\hat{\mathit{a}}$ , and $\overline{\mathit{a}}$ , respectively, represent the nodal displacement $a_A$ , the nodal temperature $T_A$ and the nodal damage $d_A$ :

$$\begin{array}{c}{\epsilon }^h(\mathrm{x})=\sum _A{B}_A^u(\mathrm{x})a_A=B_{u}a\end{array}$$

(30)

$$\begin{array}{c}\mathrm{\nabla }{T}^h(\mathrm{x})=\sum _A{B}_A^T(\mathrm{x})T_A=B_T\hat{\mathit{a}}\end{array}$$

(31)

$$\begin{array}{c}\mathrm{\nabla }{d}^h(\mathrm{x})=\sum _A{B}_A^d(\mathrm{x})d_A=B_d\overline{\mathit{a}}\end{array}$$

(32)

where the ${\mathrm{N}}_A^u$ , ${\mathrm{B}}_A^u$ , ${\mathrm{B}}_A^T$ , and ${\mathrm{B}}_A^d$ are expressed as follows:

$$\begin{array}{c}\begin{array}{c}{B}_A^u(\mathrm{x})=\left[\begin{array}{c}{\mathrm{\partial }}_x{N}_A^u\\ \begin{array}{c}0\\ {\mathrm{\partial }}_y{N}_A^u\end{array}\end{array}\begin{array}{c}0\\ \begin{array}{c}{\mathrm{\partial }}_y{N}_A^u\\ {\mathrm{\partial }}_x{N}_A^u\end{array}\end{array}\right]\\ {\mathrm{N}}_A^u(\mathrm{x})=\left[\begin{array}{cc}{N}_A^u(\mathrm{x})& 0\\ 0& N_A^u(\mathrm{x})\end{array}\right]\\ B_A^T(\mathrm{x})={\left[\begin{array}{cc}{\mathrm{\partial }}_x{N}_A^T& {\mathrm{\partial }}_y{N}_A^T\end{array}\right]}^{\mathrm{T}}\\ B_A^d(\mathrm{x})={\left[\begin{array}{cc}{\mathrm{\partial }}_x{N}_A^d& {\mathrm{\partial }}_y{N}_A^d\end{array}\right]}^{\mathrm{T}}\end{array}\end{array}$$

(33)

The ${\mathrm{\partial }}_x{N}_A^u$ means the first derivative of $N_A^u(\mathrm{x})$ . The elemental matrices $B_u$ , $B_T$ , and $B_d$ are constructed as:

$$\begin{array}{c}{B}_u=\left[\begin{array}{cccc}{B}_1^u& B_2^u& \dots & B_n^u\end{array}\right]\end{array}$$

(34)

$$\begin{array}{c}{B}_T=\left[\begin{array}{cccc}{B}_1^T& B_2^T& \dots & B_n^T\end{array}\right]\end{array}$$

(35)

$$\begin{array}{c}{B}_d=\left[\begin{array}{cccc}{B}_1^d& B_2^d& \dots & B_n^d\end{array}\right]\end{array}$$

(36)

where n represents the number of nodes per element, the matrix form of the shape functions are as follows:

$$\begin{array}{c}{N}_u=\left[\begin{array}{cccc}{N}_1^u& N_2^u& \dots & N_n^u\end{array}\right]\end{array}$$

(37)

$$\begin{array}{c}{N}_T=\left[\begin{array}{cccc}{N}_1^T& N_2^T& \dots & N_n^T\end{array}\right]\end{array}$$

(38)

$$\begin{array}{c}{N}_d=\left[\begin{array}{cccc}{N}_1^d& N_2^d& \dots & N_n^d\end{array}\right]\end{array}$$

(39)

As a result of the above work, the weakly formal equation of the finite element is semidiscretized to the form described below:

$$\begin{array}{c}{\mathit{r}}_u:={\mathit{f}}_u^{ext}-\underset{\mathrm{\Omega }}{\int }{B}_u^{\mathrm{T}}\mathit{\sigma }\mathrm{d}V=0\end{array}$$

(40)

$$\begin{array}{c}{\mathit{r}}_T:={\mathit{f}}_T^{ext}-\underset{\mathrm{\Omega }}{\int }{B}_T^{\mathrm{T}}\mathrm{h}\mathrm{d}V-C\tilde{\mathit{a}}=0\end{array}$$

(41)

$$\begin{array}{c}{\mathit{r}}_d:=-\underset{\mathrm{B}}{\overset{(}{\int }}-N_d^{\mathrm{T}}Q(\epsilon ,d)+B_d^{\mathrm{T}}\mathrm{q})\mathrm{d}V=0\end{array}$$

(42)

where ${\mathit{r}}_u$ , ${\mathit{r}}_T$ , and ${\mathit{r}}_d$ denote the residual vectors; $\tilde{\mathit{a}}=\mathrm{d}\hat{\mathit{a}}/\mathrm{dt}$ and the external forces are given by

$$\begin{array}{c}{\mathit{f}}_u^{\mathrm{ext}}=\underset{\mathrm{\Omega }}{\int }{N}_u^{\mathrm{T}}b*\mathrm{d}V+\underset{\mathrm{\partial }{\mathrm{\Omega }}_t}{\int }{N}_u^{\mathrm{T}}t*\mathrm{d}A\end{array}$$

(43)

$$\begin{array}{c}{\mathit{f}}_T^{\mathrm{ext}}=\underset{\mathrm{\Omega }}{\int }{N}_T^{\mathrm{T}}\mathrm{\Theta }*\mathrm{d}V-\underset{\mathrm{\partial }{\mathrm{\Omega }}_J}{\int }{N}_T^{\mathrm{T}}\mathrm{h}*\mathrm{d}A\end{array}$$

(44)

The consistent mass matrix $\mathit{M}$ and the heat capacity matrix $\mathit{C}$ are expressed as follows:

$$\begin{array}{c}M=\underset{\mathrm{\Omega }}{\int }\rho N_u^{\mathrm{T}}{N}_u\mathrm{d}V,C=\underset{\mathrm{\Omega }}{\int }\rho c{N}_T^{\mathrm{T}}{N}_T\mathrm{d}V\end{array}$$

(45)

Algorithmic examples

This section primarily presents a pair of representative cases to illustrate the efficacy of the length-scale insensitive degradation function within the phase-field approach to thermal cracking. Plane stress conditions are applied consistently throughout the subsequent examples. The intact material is denoted by the blue region, whereas the red region signifies the presence of cracks. The numerical cases mentioned before employ a staggered procedure to address the nonlinear coupled equations.

Quenching experiment of rectangular ceramic plate

The quenching experiment of rectangular ceramic plates ²⁶ is considered, where a rectangular ceramic plate measuring 50 mm in length and 10 mm in width is initially heated to $T_0$ (300 °C), keeping this temperature unchanged; the plate is then maintained at this temperature for a while before being subjected to free-fall motion, vertically plunging into a water bath at 20 °C. After 10 min, the rectangular plate is removed from the water bath, dried, dyed, and observed. Some scholars have studied this experiment in the literature.^23,27–29 Material parameters are shown in Table 1.

Table 1.

Parameters of quenching experiment of rectangular and circular ceramic plates.

Parameters	Values	Parameters	Values
${\sigma }_c$	180 MPa	${\eta }_t$	7.5 × 10−6 K−1
$G_f$	42.47 N/m	$\rho$	3980 kg/m3
$E_0$	370 GPa	${\kappa }_0$	31 W/(m·K)
$v_0$	0.3	c	880 J/(kg·K)

Figure 2(a) shows the final crack distribution in the ceramic plate quenching experiment. It can be seen that after the thermal shock, the cracks on the surface of the rectangular ceramic plate are almost perpendicular to the edges. Since the specimen is uniform, the cracks exhibit an axisymmetric pattern.

Figure 2.

The quenching experiment of rectangular ceramic plate: (a) The experimentally observed crack pattern; (b) The numerical model depicting the geometry and applied boundary conditions; (c) The finite element mesh diagram in COMSOL.

In the numerical simulation experiment, the thickness of the rectangular plate is set to 1 mm, the initial temperature is 300°C, and the ambient temperature is 20°C. Figure 2(b) shows the two-dimensional numerical simulation model of the ceramic plate. During the initial grid division, a uniform grid is used based on the geometric shape and dimensions of the model, as shown in Figure 2(c). The minimum grid size is set to 0.05 mm, and the complete grid contains 50,000 domain elements and 1200 boundary elements.

Figure 3.

Classical AT2 degradation functions applied to thermal shock, numerical observation of cracking patterns.

The smallest grid size has a significant impact on the simulation results, and it can even directly affect whether the results converge. In this work, the minimum grid scale is set to 0.05, which is a value we determined through experimentation during the numerical simulation process, ensuring convergence across different phase-field characteristic scales. If this minimum grid scale is further increased, the final computational results will not converge under certain phase-field scales, meaning the crack propagation predictions cannot be obtained. On the other hand, reducing the grid size further will lead to an increase in computational time as the number of computational domains grows.

In addition, the initial rectangular plate is intact, so d = 0 acts on all boundaries. The length scale ${\ell }_p$ is not an independent parameter and is set according to equation (15) which equals to $27G_f{\overline{E}}_0/256{\sigma }_c^2$ , where $G_f=42.47\mathrm{N}/\mathrm{m}$ , ${\sigma }_c^{=}180\mathrm{MPa}$ , and ${\overline{E}}_0$ is calculated by:

$$\begin{array}{c}{\overline{E}}_0=E_0{\displaystyle \frac{1-v_0}{(1-2v_0)(1+v_0)}=370\times {\displaystyle \frac{1-0.3}{(1-0.6)(1+0.6)}=498.0769}}\end{array}$$

(46)

So that ${\ell }_p=(27\times 42.47\times 498.0769)/(256\times 180\times 180)=0.068$ .

Subsequently, we incorporate the degradation function that is unaffected by the length scale into our model. The length scale ${\ell }_p$ of the physical process area is always equal to 0.068 mm and remains unchanged, while the phase-field scale parameter ${\ell }_0$ changes with the change of s. Figures 4–7 show the fracture paths under the influence of various degradation functions.

Figure 4.

s = 200, ${\ell }_p=0.068\mathrm{mm}$ , ${\ell }_0={\ell }_p=0.068\mathrm{mm}$ , $\mathrm{h}=0.05\mathrm{mm}$ , ceramic plate of dimensions is 50 mm × 10 mm × 1 mm.

Figure 5.

s = 5, ${\ell }_p=0.068\mathrm{mm}$ , ${\ell }_0=1.12\times {\ell }_p=0.076\mathrm{mm}$ , $\mathrm{h}=0.05\mathrm{mm}$ , ceramic plate of dimensions is 50 mm × 10 mm × 1 mm.

Figure 6.

s = 1.1, ${\ell }_p=0.068\mathrm{mm}$ , ${\ell }_0=4{\ell }_p=0.272\mathrm{mm}$ , $\mathrm{h}=0.05\mathrm{mm}$ , ceramic plate of dimensions 100 mm × 20 mm × 1 mm.

Figure 7.

s = 1.04, ${\ell }_p=0.068\mathrm{mm}$ , ${\ell }_0=8{\ell }_p=0.544\mathrm{mm}$ , $\mathrm{h}=0.05\mathrm{mm}$ , ceramic plate of dimensions 140 mm × 28 mm × 1 mm.

The crack distribution at t = 100 ms is taken as the final crack morphology simulation result. Cracks with lengths exceeding 0.25 mm are considered valid cracks. The distribution of the number of valid cracks under the above parameter settings is shown in the following (Figure 8). The orange and green bar charts represent the number of valid cracks on the long and short edges of the rectangular ceramic, respectively.

Figure 8.

Comparison of the thermal crack number obtained from the present numerical tests and the previous experimental observations. ²⁶

It can be seen that when AT2 and s = 200, the numerical simulation results are closest to the experimental results. Except for these two cases, the number of valid cracks is the closest to the experimental results when s = 1.1. In terms of computation time, the smaller the s value, the less time is consumed in the calculation. In this case, when s = 200, the computation time is 24,194 s, whereas when s = 1.1, it only takes 16,975 s, resulting in a 29.8% improvement in computational efficiency. Therefore, considering both computational efficiency and time, using s = 1.1 is more optimal in this case. In addition, during the experiment, shorter cracks may nucleate at the edge of the plate, and the simulation captures some key features: a series of periodic parallel short cracks nucleate, followed by selective arrest and growth.

It is evident from Figure 9 that quenching process induces a significant temperature gradient in the outer layer while maintaining little temperature variance internally. The propagation of cracks exhibits notable stability, with their tips evolving along with the edges of the central high-temperature region, serving to alleviate stress around the crack tip. Numerical simulations demonstrate the relevance of this degradation process in the thermal fracturing of ceramic plates. The closely resemblance between crack networks generated through stochastic simulations and real crack formations aids in a deeper comprehension of crack development resulting from thermal shock.

Figure 9.

Temperature and maximum principal stress plots of rectangular ceramic plates at t = 0.1 s.

2-D thin circular ceramic specimen in quenching tests

The quenching test of 2-D thin circular ceramic plates ³⁰ is considered; 99% ${\mathrm{Al}}_2{\mathrm{O}}_3$ powder was thermoformed into thin circular plates with dimensions of 13-mm diameter and 1-mm thickness. The bound specimens were heated to the preset temperature $T_0$ and maintained at this temperature for 30 min. The $T_0$ was 400°C. After the model is heated, it is put in free fall into a water bath at 20°C. After 10 minutes, dry the specimen from the water bath, then dye it with blue ink to observe the crack state.

From the experiment result shown in Figure 10, it can be seen that after the specimen is subjected to thermal shock, the crack formed by thermal shock is perpendicular to the circular edge and points to the center. At the same time, we can also observe that the crack pattern on both sides of the specimen is the same, indicating that the crack pattern presented is axisymmetric, which we can easily observe. This issue has been studied in the literature. ²⁹ Material parameters are shown in Table 1. Figures 11–15 show the simulation results of applying different degradation functions acting on circular ceramic plates for thermal shock experiments.

Figure 10.

Thermal shock crack patterns, where the preset temperature $T_0=400{}^{\mathrm{o}}\mathrm{C}$ .

Figure 11.

Classical AT2 degradation function applied to circular ceramic plates, numerically observed crack pattern.

Figure 12.

s = 200, ${\ell }_p=0.068\mathrm{mm}$ , ${\ell }_0={\ell }_p=0.068\mathrm{mm}$ , $\mathrm{h}=0.05\mathrm{mm}$ , ceramic plate of radius remains unchanged.

Figure 13.

s = 5, ${\ell }_p=0.068\mathrm{mm}$ , ${\ell }_0=1.12{\ell }_p=0.076\mathrm{mm}$ , $\mathrm{h}=0.05\mathrm{mm}$ , ceramic plate of radius remains unchanged.

Figure 14.

s = 1.1, ${\ell }_p=0.068\mathrm{mm}$ , ${\ell }_0=4{\ell }_p=0.272\mathrm{mm}$ , $\mathrm{h}=0.05\mathrm{mm}$ , ceramic plate of radius is 13 mm.

Figure 15.

s = 1.04, ${\ell }_p=0.068\mathrm{mm}$ , ${\ell }_0=8{\ell }_p=0.544\mathrm{mm}$ , $\mathrm{h}=0.05\mathrm{mm}$ , ceramic plate of radius is 18.33 mm.

Figures 11–15 illustrate various stages of the fracture process. Observing these images allows for a summarized understanding of the pattern of crack propagation. Initially, thermal shock initiates cracks that propagate uniformly with nearly equal spacing (Figures 11–15(a)). At this stage, crack propagation is rapid, gradually decelerating as thermal stress dissipates, until the number of cracks exceeds the range supported by strain energy. Consequently, only a few cracks continue expanding while others cease propagation, and the distance between cracks increases (Figures 11–15(b)). Subsequently, crack propagation may deviate, attempting to establish equidistant spacing. This process may repeat multiple times until the strain induced by thermal stress can no longer support the expansion of any further cracks. Ultimately, the simulated crack patterns closely resemble the observed patterns (Figures 11–15(c)).

It can be seen from Figure 16 that when the temperature gradient reaches a certain level, cracks initiate from the position of maximum principal stress. Subsequently, stress relaxes at that point, leading to stress redistribution and changes in the location of maximum principal stress, facilitating continual crack initiation. The propagation of cracks arrest when the temperature gradient no longer induces sufficient stress concentration. During thermal shock, stress is uniformly distributed on the specimen's surface, resulting in a certain degree of randomness in the location of maximum principal stress.

Figure 16.

Temperature and maximum principal stress plots of rectangular ceramic plates at t = 0.1 s.

The data in Table 1 are all from.^23,26,28 This work assumes that the material parameters remain constant despite temperature changes. The initial temperatures chosen for analysis are $T_0$  = 300°C, 400°C; see the literature ²¹ for more cases. The specific value of the physical length scale ${\ell }_p$ is set at 0.068 mm, whereas the phase-field length scale ${\ell }_0$ varies depending on the parameter s. For s values of 200, 5, 1.1, and 1.04, the corresponding phase-field length scales are 1, 1.12, 4, and 8, respectively. The simulation time is set to 0.1 ms.

By decreasing the value of s, the phase-field length scale ${\ell }_0$ increase. To obtain an accurate solution, the mesh density is typically set to about 1/5 of the phase-field's length scale. The advantage of selecting a smaller value of s for simulating the same model size is that it allows the corresponding proportion of the phase-field's length scale to be achieved, which significantly reduces the simulation time and enables the correct crack propagation distribution to be determined more quickly. The simulation diagram shows excellent agreement with the observed experimental results. Initially, regardless of whether the plates are rectangular or circular, cracks begin to nucleate at the edges of these plates. As the temperature continues to decrease, some cracks arrested to grow, while others continue to expand. Eventually, the spacing between cracks becomes relatively uniform.

Conclusions

The study effectively integrates a degradation function that is insensitive to the length scale into the phase-field method for modeling thermal cracking. Moreover, the displacement equations, damage equations, and the heat conduction process in the solid are introduced. The research introduces a degradation function that effectively disentangles the physical dimensions from the phase-field dimensions in the context of thermal cracking, leading to a reduced necessity for mesh density to resolve the diffusive zone. To a certain extent, it breaks through the limitation of the original small scale. Through some numerical examples that we are particularly interested in, the crack evolution, nucleation, propagation, and arrest behaviors of cracks are well captured, the applicability of the length-scale insensitive degradation function to the phase-field method of thermal cracking is verified, which opens the door for the extension of phase-field fracture theory to the simulation of large engineering structures and also has a particular guiding role for the fracture simulation of other multiphysical field problems.

In addition, the degree of decoupling between the phase-field length scale and the physical length scale of the material is limited, and further implementation of engineering simulations at larger scales needs to be further investigated. Like brittle fracture, the phase-field method also holds great potential in the field of ductile fracture. In the future, we will also explore the application of length-scale-insensitive degradation functions to ductile fracture models to expand the range of the model's application.