Phase Field Models for Thermal Fracturing and Their Variational Structures

Abstract

It is often observed that thermal stress enhances crack propagation in materials, and, conversely, crack propagation can contribute to temperature shifts in materials. In this study, we first consider the thermoelasticity model proposed by M. A. Biot and study its energy dissipation property. The Biot thermoelasticity model takes into account the following effects. Thermal expansion and contraction are caused by temperature changes, and, conversely, temperatures decrease in expanding areas but increase in contracting areas. In addition, we examine its thermomechanical properties through several numerical examples and observe that the stress near a singular point is enhanced by the thermoelastic effect. In the second part, we propose two crack propagation models under thermal stress by coupling a phase field model for crack propagation and the Biot thermoelasticity model and show their variational structures. In our numerical experiments, we investigate how thermal coupling affects the crack speed and shape. In particular, we observe that the lowest temperature appears near the crack tip, and the crack propagation is accelerated by the enhanced thermal stress.

1. Introduction

Cracking is a phenomenon that occurs everywhere in our lives, but, if it is allowed to continue, it can cause fatal damage. A crack in a material occurs when the material experiences a continuous overload. However, several other factors, such as thermal expansion and contraction due to temperature changes [1,2,3], fluid pressure (e.g., in hydraulic fracturing) [4], the diffusion of hydrogen (or hydrogen embrittlement) [5,6], chemical reactions [7], and humidity [2], cause cracks in materials. In particular, among these phenomena, cracks due to thermal expansion are interesting to study from the viewpoint of the energy balance between elastic, thermal, and surface energies.

M. A. Biot proposed a theoretical framework for coupled thermoelasticity based on the principle of minimum entropy production [8]. Biot’s model is now widely known as the traditional coupled thermoelasticity model, and it has been extended to dynamical theory [9] and to various other situations [10,11,12,13,14,15]. As shown in Section 2.2, it satisfies an energy balance equality between the elastic and thermal energies.

In fracture mechanics, especially in the modeling and simulation of crack propagation, a phase field approach has been recently recognized as a powerful tool. The phase field model (PFM) for fractures was first proposed by Bourdin et al. [16] and Karma et al. [17]. Then, based on the framework of variational fracture theory [18,19], the techniques and applications of PFM have been extensively developed, for example [20,21,22,23,24,25]. We refer to [26] for further information on the development of PFM for fracture mechanics. PFM for fracture mechanics is derived as a gradient flow of the total energy, which consists of the elastic energy and the surface energy and is known to be consistent with the classical Griffith theory [16,26]. It allows us to handle the complex geometry of multiple, kinked, or branching cracks in both 2D and 3D without a crack path search. Comparisons with the experimental results are investigated in [27].

The aim of our paper is propose an energy-consistent PFM for thermal fracturing by the coupling Biot thermoelasticity model. Naturally, three kinds of energy, i.e., elastic, thermal, and surface energies, appear in our stage, and the exchange and dissipation of those energies are the main interests of our research. An illustration is shown in Figure 1. There are several previous works that address thermal fracturing using PFM [26,28,29,30,31]. In particular, Miehe et al. [23] developed a theory of thermomechanical fracture with a diffusive crack model including various nonlinear effects and demonstrated several interesting numerical examples. However, the strain’s influence on the heat transfer as the Biot model is not involved in those previous works. It means that the previous model can not capture thermal distribution during crack propagation. To the best of our knowledge for the condition, a peridynamics model that employs the coupled thermoelastic equation was proposed by Gao and Oterkus [32]. In this study, for each PDE model, we consider an initial-value and boundary-value problem in a bounded domain ${\mathbb{R}}^d$$(d=2\mathrm{or}3)$ and derive its energy equality which represents the energy dissipation property.

The organization of this paper is as follows: in Section 2, we introduce the linear thermoelasticity model by M.A. Biot and derive its variational principle and energy dissipation property. In addition, we numerically investigate the effect of the thermal coupling term on the elastic and thermoelastic energies in an expanding region.

Section 3 is devoted to PFMs for crack propagation under thermal stress. In Section 3.1, we give a brief review of the irreversible fracturing phase field model (F-PFM) and its energy equality, which guarantees the energy dissipation property (Theorem 2) and follows the works [22,26]. In Section 3.2 and Section 3.3, we propose two types of thermal fracturing phase field models (TF-PFMs). The first model, TF-PFM1, is a straightforward coupling of F-PFM and the Biot thermoelasticity model. Based on the variational principle of the Biot model (Proposition 2), we show a partial energy equality for a fixed temperature (Theorem 3). However, it does not satisfy the energy equality for the total energy, which consists of the elastic, thermal, and surface energies.

The second model, TF-PFM2, presented in Section 3.3 is another natural coupling of F-PFM and the Biot thermoelasticity model based on the energy equality of the Biot model (Theorem 1). We prove an energy equality for TF-PFM2 in Theorem 4. Since we consider several models (Biot’s model, F-PFM, and TF-PFMs) and their energy qualities, for the readers’ convenience, we list the energies and energy equalities for each model in

Table 1. List of energies.

Type of Energy	Definition	Equation
Elastic	${\displaystyle E_{el}\left(u\right):={\displaystyle \frac{1}{2}}{\int }_{\mathsf{\Omega }}\sigma \left[u\right]:e\left[u\right]dx}$	(8)
Thermoelastic	${\displaystyle E_{el}^{*}(u,\mathsf{\Theta }):={\displaystyle \frac{1}{2}}{\int }_{\mathsf{\Omega }}{\sigma }^{*}[u,\mathsf{\Theta }]:e^{*}[u,\mathsf{\Theta }]dx}$	(10)
Thermal	${\displaystyle E_{th}\left(\mathsf{\Theta }\right):={\displaystyle \frac{\chi }{2{\mathsf{\Theta }}_0}}{\int }_{\mathsf{\Omega }}|\mathsf{\Theta }\left(x\right)-{\mathsf{\Theta }}_0{|}^{2}dx}$	(11)
Modified elastic	${\displaystyle {\mathcal{E}}_{el}(u,z):={\displaystyle \frac{1}{2}}{\int }_{\mathsf{\Omega }}{(1-z)}^2\sigma \left[u\right]:e\left[u\right]dx}$	(24)
Modified thermoelastic	${\displaystyle {\mathcal{E}}_{el}^{*}(u,\mathsf{\Theta },z):={\displaystyle \frac{1}{2}}{\int }_{\mathsf{\Omega }}{(1-z)}^2{\sigma }^{*}[u,\mathsf{\Theta }]:e^{*}[u,\mathsf{\Theta }]dx}$	(29)
Surface	${\displaystyle E_s\left(z\right):={\displaystyle \frac{1}{2}}{\int }_{\mathsf{\Omega }}{\gamma }_{*}\left({ϵ|\nabla z|}^2+\frac{{\left|z\right|}^2}{ϵ}\right)dx}$	(25)

and

Table 2. Different forms of energy equalities.

Model	Strong Form	Energy	Energy Equality
Linear elasticity	(6)	$E_{el}\left(u\right)$	-
Biot’s model	(1) and (2)	$E_{el}\left(u\right)+E_{th}\left(\mathsf{\Theta }\right)$	(12)
F-PFM	(22a) and (22b)	${\mathcal{E}}_{el}(u,z)+E_s\left(z\right)$	(27)
TF-PFM1	(31a) and (31c)	${\mathcal{E}}_{el}^{*}(u,\mathsf{\Theta },z)+E_s\left(z\right)$	(33) a
TF-PFM2	(35a) and (35c)	${\mathcal{E}}_{el}(u,z)+E_s\left(z\right)+E_{th}\left(\mathsf{\Theta }\right)$	(36)

^a When a temperature $\mathsf{\Theta }=\mathsf{\Theta }\left(x\right)\in L^2\left(\mathsf{\Omega }\right)$ is given.

.

In Section 4, we show some numerical comparisons between two TF-PFMs using non-dimensionalized equations. We investigate the effects of the thermal coupling in TF-PFM1 and TF-PFM2 on the crack speed and the crack path by changing a dimensionless coupling parameter $\delta$. As noted, although the temperature influences material properties micro-structurally [33], it is not considered in the present study for simplicity. Generally, the effect of micro-structure of material gives the typical crack path, such as: curving, branching, kinking, etc. A clear study of this is addressed by [34]. The last section shows some conclusions and comments on further topics.

To easily understand the relevant notation and symbols in this paper, we introduce them in this section. Let $\mathsf{\Omega }$ be a bounded domain in ${\mathbb{R}}^d$ ($d=2$ or 3). The position in ${\mathbb{R}}^d$ is denoted by $x={(x_1.\dots ,x_d)}^T\in {\mathbb{R}}^d$, where ${ }^T$ denotes the transposition of a vector or matrix. Let ∇, $\mathrm{div}$, and $\mathsf{\Delta }$ be the gradient, divergence, and Laplacian operators with respect to x, respectively. For simplicity, we write $\dot{u}$, $\dot{\mathsf{\Theta }}$, and $\dot{z}$ as the partial derivatives of u, $\mathsf{\Theta }$ and z with respect to t, respectively. For simplicity, we often denote $u\left(t\right):=u(·,t)$, etc. The space of the real-valued (symmetric) $d\times d$ matrix is denoted by ${\mathbb{R}}^{d\times d}$ (${\mathbb{R}}_{sym}^{d\times d}$). The inner product of square matrices $A,B\in {\mathbb{R}}^{d\times d}$ is denoted by $A:B:={\sum }_{i,j=1}^d{A}_{ij}{B}_{ij}$. Using $L^2\left(\mathsf{\Omega }\right)$, we refer to the Lebesgue space on $\mathsf{\Omega }$, while $H^1(\mathsf{\Omega },{\mathbb{R}}^d)$ and $H^{\frac{1}{2}}({\mathsf{\Gamma }}_D^u,{\mathbb{R}}^d)$ represent the Sobolev space on $\mathsf{\Omega }$ and its trace space on the boundary ${\mathsf{\Gamma }}_D^u$, respectively. For more details on Sobolev spaces, we refer to the review in [35]. In addition, we summarize the physical properties used in this paper in

Table 3. List of physical properties.

Symbol	Physical Meaning [Unit]	Symbol	Physical Meaning [Unit]
u	Displacement [$\mathrm{m}$]	${\sigma }^{*}[u,\mathsf{\Theta }]$	Stress tensor with thermal effect [$\mathrm{Pa}$]
$\mathsf{\Theta }$	Temperature [$\mathrm{K}$]	$e^{*}[u,\mathsf{\Theta }]$	Strain tensor with thermal effect [-]
${\mathsf{\Theta }}_0$	Reference temperature [$\mathrm{K}$]	$\beta$	Stress thermal modulus [$\mathrm{Pa}$ · $\mathrm{K}$−1]
z	Damage variable [-]	${\kappa }_0$	Thermal conductivity [$\mathrm{W}$ · $\mathrm{m}$−1 · $\mathrm{K}$−1]
$\sigma \left[u\right]$	Stress tensor [$\mathrm{Pa}$]	$\chi$	Volumetric heat capacity [$\mathrm{J}$ · $\mathrm{K}$−1 · $\mathrm{m}$−3]
$e\left[u\right]$	Strain tensor [-]	$a_L$	Coefficient of linear thermal expansion [$\mathrm{K}$−1]
$E_{\mathrm{Y}}$	Young’s modulus [$\mathrm{Pa}$]	$\delta$	Thermoelastic coupling parameter [-]
${\nu }_{\mathrm{P}}$	Poisson ratio [-]	${\gamma }_{*}$	Critical energy release rate a  [$\mathrm{Pa}·\mathrm{m}$]
$\lambda$, $\mu$	Lamé’s constants b  [$\mathrm{Pa}$]	$ϵ$	Length scale in F-PFM or TF-PFM [$\mathrm{m}$]
t	Time [$\mathrm{s}$]	$\alpha$	Time regularization parameter in F-PFM or TF-PFM [$\mathrm{Pa}·\mathrm{s}$]

^a ${\gamma }_{*}$ is usually denoted by ${\mathrm{G}}_{\mathrm{c}}$ [26,36]. ^b $\lambda$ and $\mu$ are written as $\lambda =\frac{E_{\mathrm{Y}}{\nu }_{\mathrm{P}}}{(1+{\nu }_{\mathrm{P}})(1-2{\nu }_{\mathrm{P}})}$ and $\mu =\frac{E_{\mathrm{Y}}}{2(1-{\nu }_{\mathrm{P}})}$.

.

2. Thermoelasticity Model

2.1. Formulation of the Problem

M.A. Biot [8] proposed the following mathematical model for coupled thermoelasticity:

$$\begin{array}{cc}\mathrm{div}\sigma \left[u\right]=\beta \nabla \mathsf{\Theta }\hfill & \hfill \mathrm{in}\mathsf{\Omega }\times [0,T],\end{array}$$

(1)

$$\begin{array}{cc}\chi {\displaystyle \frac{\partial }{\partial t}}\mathsf{\Theta }={\kappa }_0\mathsf{\Delta }\mathsf{\Theta }-{\mathsf{\Theta }}_0\beta {\displaystyle \frac{\partial }{\partial t}}\left(\mathrm{div}u\right)\hfill & \hfill \mathrm{in}\mathsf{\Omega }\times (0,T],\end{array}$$

(2)

where $\mathsf{\Omega }$ is a bounded domain in ${\mathbb{R}}^d$ ($d=2$ or 3). We suppose that $\mathsf{\Omega }$ is an isotropic elastic body and consider the thermoelastic coupling between the mechanical deformation and the thermal expansion in $\mathsf{\Omega }$. The constant $\beta$ is defined by $\beta :=a_L(d\lambda +2\mu )$ with $a_L>0$ as the coefficient of linear thermal expansion and $\mu (>0)$; $\lambda (>-\frac{2\mu }{d})$ are Lamé’s constants.

The unknown functions in (1) and (2) are the displacement $u(x,t)=(u_1(x,t),$$\dots ,u_d{(x,t))}^T\in {\mathbb{R}}^d$ and the temperature $\mathsf{\Theta }(x,t)\in \mathbb{R}$. In addition, the constant ${\mathsf{\Theta }}_0>0$ is a fixed reference temperature. Similarly, strain $e\left[u\right]$ and stress tensors $\sigma \left[u\right]$ are defined as

$$\begin{array}{cc}e\left[u\right]:=\frac{1}{2}\left(\nabla u^T+{(\nabla u^T)}^T\right)\in {\mathbb{R}}_{sym}^{d\times d},\hfill & \hfill \left(3\mathrm{a}\right)\\ \sigma \left[u\right]:=Ce\left[u\right]=\lambda \left(\mathrm{div}u\right)I+2\mu e\left[u\right]\in {\mathbb{R}}_{sym}^{d\times d},\hfill & \hfill \left(3\mathrm{b}\right)\end{array}$$

where $C:=\left(c_{ijkl}\right),c_{ijkl}=\lambda {\delta }_{ij}{\delta }_{kl}+\mu ({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})$ is an isotropic elastic tensor and I is the identity matrix of size d. From (3b), (1) is also written in the form

$$\begin{array}{cc}\hfill & \mu \mathsf{\Delta }u+(\lambda +\mu )\nabla \left(\mathrm{div}u\right)=\beta \nabla \mathsf{\Theta }.\hfill \end{array}$$

The term $\beta \nabla \mathsf{\Theta }$ in (1) and the term ${\mathsf{\Theta }}_0\beta \frac{\partial }{\partial t}\left(\mathrm{div}u\right)$ in (2) represent the body force due to thermal expansion and the heat source due to the volume change rate, respectively. We remark that, when $a_L=0$, (1) and (2) are decoupled.

It is convenient to introduce the following strain and stress tensors, including the thermal effect:

$$\begin{array}{cc}{e}^{*}[u,\mathsf{\Theta }]:=e\left[u\right]-a_L(\mathsf{\Theta }(x,t)-{\mathsf{\Theta }}_0)I\in {\mathbb{R}}_{sym}^{d\times d},\hfill & \hfill \left(4\mathrm{a}\right)\\ {\sigma }^{*}[u,\mathsf{\Theta }]:=C{e}^{*}[u,\mathsf{\Theta }]=\sigma \left[u\right]-\beta (\mathsf{\Theta }(x,t)-{\mathsf{\Theta }}_0)I\in {\mathbb{R}}_{sym}^{d\times d}.\hfill & \hfill \left(4\mathrm{b}\right)\end{array}$$

Using the thermal stress tensor ${\sigma }^{*}[u,\mathsf{\Theta }]$, (1) can be written in the following form:

$$\begin{array}{cc}\hfill & -\mathrm{div}{\sigma }^{*}[u,\mathsf{\Theta }]=0.\hfill \end{array}$$

This means that the force ${\sigma }^{*}[u,\mathsf{\Theta }]$ is in equilibrium in $\mathsf{\Omega }$. In the preceding, Equations (1) and (2) represent the force balance and the thermal diffusion in $\mathsf{\Omega }$, respectively.

The system in (1) and (2) is complemented by the following boundary and initial conditions:

$$\{\begin{array}{ccc}u=u_D(x,t)\hfill & \mathrm{on}{\mathsf{\Gamma }}_D^u\times [0,T],\hfill & \hfill \left(5\mathrm{a}\right)\\ {\sigma }^{*}[u,\mathsf{\Theta }]n=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_N^u\times [0,T],\hfill & \hfill \left(5\mathrm{b}\right)\\ \mathsf{\Theta }={\mathsf{\Theta }}_D(x,t)\hfill & \mathrm{on}{\mathsf{\Gamma }}_D^{\mathsf{\Theta }}\times [0,T],\hfill & \hfill \left(5\mathrm{c}\right)\\ \frac{\partial \mathsf{\Theta }}{\partial n}=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_N^{\mathsf{\Theta }}\times [0,T],\hfill & \hfill \left(5\mathrm{d}\right)\\ \mathsf{\Theta }(x,0)={\mathsf{\Theta }}_{*}\left(x\right)\hfill & \mathrm{in}\mathsf{\Omega },\hfill & \hfill \left(5\mathrm{e}\right)\end{array}$$

where n is the outward unit normal vector along the boundary, $\mathsf{\Gamma }={\mathsf{\Gamma }}_D^u\cup {\mathsf{\Gamma }}_N^u$$(\mathsf{\Gamma }={\mathsf{\Gamma }}_D^{\mathsf{\Theta }}\cup {\mathsf{\Gamma }}_N^{\mathsf{\Theta }})$ with ${\mathsf{\Gamma }}_D^u\cap {\mathsf{\Gamma }}_N^u=\varnothing ({\mathsf{\Gamma }}_D^{\mathsf{\Theta }}\cap {\mathsf{\Gamma }}_N^{\mathsf{\Theta }}=\varnothing )$. The boundaries ${\mathsf{\Gamma }}_D^u$ and ${\mathsf{\Gamma }}_N^u$ (${\mathsf{\Gamma }}_D^{\mathsf{\Theta }}$ and ${\mathsf{\Gamma }}_N^{\mathsf{\Theta }}$) are the Dirichlet and Neumann boundaries for u (for $\mathsf{\Theta }$), respectively. We suppose that the $(d-1)$-dimensional volume of ${\mathsf{\Gamma }}_D^u$ is positive for the solvability of u.

Remark 1.

Instead of boundary conditions (5a) and (5b), we can also consider the following mixed-type condition. When $d=2$, on a part of the boundary (which we denote by ${\mathsf{\Gamma }}_{DN}^u$), $u={(u_1,u_2)}^T$ and

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{u}_1=u_{D1}\hfill & \mathrm{on}{\mathsf{\Gamma }}_{DN}^u,\hfill \\ \left({\sigma }^{*}[u,\mathsf{\Theta }]n\right)·e_2=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_{DN}^u,\hfill \end{array}\right.\end{array}$$

or

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{u}_2=u_{D2}\hfill & \mathrm{on}{\mathsf{\Gamma }}_{DN}^u,\hfill \\ \left({\sigma }^{*}[u,\mathsf{\Theta }]n\right)·e_1=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_{DN}^u,\hfill \end{array}\right.\end{array}$$

where $u_{Di}:={\mathsf{\Gamma }}_{DN}^u\mapsto \mathbb{R}$ is a given horizontal or vertical displacement and $e_1={(1,0)}^T$, $e_2={(0,1)}^T$. These types of mixed boundary conditions are considered in Section 2.3.3 and Section 4.4.1. Even for these mixed-type boundary conditions, we can easily extend the following arguments on weak solutions, variational principles, and energy equalities.

2.2. Variational Principle and Energy Equality

This section aims to show a variational principle and provide an energy equality that implies the energy dissipation property for the system (1) and (2). In linear elasticity theory, a weak form of the boundary value problem for $u_D\in H^{\frac{1}{2}}({\mathsf{\Gamma }}_D^u;{\mathbb{R}}^d)$ is

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}-\mathrm{div}\sigma \left[u\right]=0\hfill & \mathrm{in}\mathsf{\Omega },\hfill \\ u=u_D\hfill & \mathrm{on}{\mathsf{\Gamma }}_D^u,\hfill \\ \sigma \left[u\right]n=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_N^u,\hfill \end{array}\right.\end{array}$$

(6)

which is given by

$$\begin{array}{cc}\hfill & u\in V^u\left(u_D\right),{\int }_{\mathsf{\Omega }}\sigma \left[u\right]:e\left[v\right]dx=0\hspace{1em}\mathrm{for}\mathrm{all}v\in V^u\left(0\right),\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & V^u\left(u_D\right):=\left\{u\in H^1(\mathsf{\Omega };{\mathbb{R}}^d);u{|}_{{\mathsf{\Gamma }}_D^u}=u_D\right\}.\hfill \end{array}$$

(7)

A weak solution uniquely exists and is given by

$$u=\underset{v\in V^u\left(u_D\right)}{argmin}{E}_{el}\left(v\right),$$

where

$$\begin{array}{cc}\hfill & E_{el}\left(v\right):=\frac{1}{2}{\int }_{\mathsf{\Omega }}\sigma \left[v\right]:e\left[v\right]dx\hspace{1em}(v\in H^1\left(\mathsf{\Omega };{\mathbb{R}}^d)\right)\hfill \end{array}$$

(8)

is an elastic energy. This is known as a variational principle [37,38]. For a fixed $\mathsf{\Theta }\left(x\right)$, a weak form for u of (1) and its variational principle are derived as follows:

Proposition 1.

For $u\in H^2(\mathsf{\Omega };{\mathbb{R}}^d)$ and $\mathsf{\Theta }\in H^1\left(\mathsf{\Omega }\right)$,

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}-\mathrm{div}{\sigma }^{*}[u,\mathsf{\Theta }]=0\hfill & \mathrm{in}\mathsf{\Omega },\hfill \\ u=u_D\hfill & \mathrm{on}{\mathsf{\Gamma }}_D^u,\hfill \\ {\sigma }^{*}[u,\mathsf{\Theta }]n=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_N^u,\hfill \end{array}\right.\end{array}$$

is equivalent to the following weak form:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\displaystyle {\int }_{\mathsf{\Omega }}{\sigma }^{*}[u,\mathsf{\Theta }]:e\left[v\right]dx=0}\hfill & \mathrm{for}\mathrm{all} v\in V^u\left(0\right),\hfill \\ u\in V^u\left(u_D\right).\hfill \end{array}\right.\end{array}$$

(9)

Proof. 

For $v\in V^u\left(0\right)$, we have

$$\begin{array}{ccc}\hfill {\int }_{\mathsf{\Omega }}\left(-\mathrm{div}{\sigma }^{*}[u,\mathsf{\Theta }]\right)·v dx& =& {\int }_{\mathsf{\Omega }}{\sigma }^{*}[u,\mathsf{\Theta }]:e\left[v\right]dx-{\int }_{{\mathsf{\Gamma }}_N^u}\left({\sigma }^{*}[u,\mathsf{\Theta }]n\right)·v ds.\hfill \end{array}$$

The equivalency immediately follows from this equation: □

Proposition 2

(Variational principle). For a given $\mathsf{\Theta }\in L^2\left(\mathsf{\Omega }\right)$, $u_D\in H^{\frac{1}{2}}({\mathsf{\Gamma }}_D^u;{\mathbb{R}}^d)$, there exists a unique weak solution $u\in H^1(\mathsf{\Omega };{\mathbb{R}}^d)$ that satisfies (9). Furthermore, the solution u is a unique minimizer of the variational problem:

$$u=\underset{v\in V^u\left(u_D\right)}{argmin}{E}_{el}^{*}(v,\mathsf{\Theta }),$$

where

$$\begin{array}{cc}\hfill & E_{el}^{*}(v,\mathsf{\Theta })=\frac{1}{2}{\int }_{\mathsf{\Omega }}{\sigma }^{*}[v,\mathsf{\Theta }]:e^{*}[v,\mathsf{\Theta }]dx.\hfill \end{array}$$

(10)

We remark that $E_{el}^{*}(v,\mathsf{\Theta })$ represents thermoelastic energy.

Proof. 

The unique existence of a weak solution for u is shown by the Lax–Milgram theorem [38] since (9) is written as

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle {\int }_{\mathsf{\Omega }}\sigma \left[u\right]:e\left[v\right]dx={\int }_{\mathsf{\Omega }}\beta (\mathsf{\Theta }-{\mathsf{\Theta }}_0)\mathrm{div}v dx,}\hfill \\ u\in V^u\left(u_D\right)\hspace{1em}(\mathrm{for}\mathrm{all} v\in V^u\left(0\right)).\hfill \end{array}\right.\end{array}$$

The coercivity of the above weak form is known as Korn’s second inequality [38]:

$$\begin{array}{cc}\hfill & {}^{\exists }{a}_0>0\mathrm{such}\mathrm{that}{\int }_{\mathsf{\Omega }}\sigma \left[v\right]:e\left[v\right]dx\ge a_0{\parallel v\parallel }_{H^1(\mathsf{\Omega };{\mathbb{R}}^2)}^2,\mathrm{for}\mathrm{all} v\in V^u\left(0\right).\hfill \end{array}$$

For a weak solution u and any $v\in V^u\left(0\right)$, using the equalities

$$\begin{array}{cc}\hfill & {\sigma }^{*}[u+v,\mathsf{\Theta }]={\sigma }^{*}[u,\mathsf{\Theta }]+\sigma \left[v\right],\hfill \\ \hfill & e^{*}[u+v,\mathsf{\Theta }]=e^{*}[u,\mathsf{\Theta }]+e\left[v\right],\hfill \\ \hfill & {\sigma }^{*}[u,\mathsf{\Theta }]:e\left[v\right]=e^{*}[u,\mathsf{\Theta }]:\sigma \left[v\right],\hfill \end{array}$$

we have

$$\begin{array}{cc}\hfill E_{el}^{*}(u+v,\mathsf{\Theta })-E_{el}^{*}(u,\mathsf{\Theta })& =\frac{1}{2}{\int }_{\mathsf{\Omega }}{\sigma }^{*}[u+v,\mathsf{\Theta }]:e^{*}[u+v,\mathsf{\Theta }]dx-\frac{1}{2}{\int }_{\mathsf{\Omega }}{\sigma }^{*}[u,\mathsf{\Theta }]:e^{*}[u,\mathsf{\Theta }]dx\hfill \\ \hfill & ={\int }_{\mathsf{\Omega }}{\sigma }^{*}[u,\mathsf{\Theta }]:e\left[v\right]dx+\frac{1}{2}{\int }_{\mathsf{\Omega }}\sigma \left[v\right]:e\left[v\right]dx\hfill \\ \hfill & =\frac{1}{2}{\int }_{\mathsf{\Omega }}\sigma \left[v\right]:e\left[v\right]dx\ge 0.\hfill \end{array}$$

This shows that u is a minimizer of $E_{el}^{*}(u,\mathsf{\Theta })$ among $V^u\left(u_D\right)$.

On the other hand, if u is a minimizer, the first variation of $E_{el}^{*}$ vanishes at u; i.e., for all $v\in V^u\left(0\right)$, we have

$$\begin{array}{cc}\hfill & 0=\frac{d}{ds}{E}_{el}^{*}(u+sv,\mathsf{\Theta }){|}_{s=0}={\int }_{\mathsf{\Omega }}{\sigma }^{*}[u,\mathsf{\Theta }]:e\left[v\right]dx.\hfill \end{array}$$

Hence, u is a weak solution. Summarizing the above, there exists a unique weak solution to (8), and u is a weak solution if and only if it is a minimizer of $E_{el}^{*}$ among $V^u\left(u_D\right)$. □

The next theorem represents a dissipation of the sum of the elastic and thermal energies during the thermomechanical process. We define thermal energy as

$$\begin{array}{c}\hfill E_{th}\left(\mathsf{\Theta }\right):=\frac{\chi }{2{\mathsf{\Theta }}_0}{\int }_{\mathsf{\Omega }}{\left|\mathsf{\Theta }\left(x\right)-{\mathsf{\Theta }}_0\right|}^{2}dx.\end{array}$$

(11)

Theorem 1

(Energy equality for Biot’s model). Let $\left(u\right(x,t),\mathsf{\Theta }(x,t\left)\right)$ be a sufficiently smooth solution to (1), (2) and (5). In addition, we suppose that $u_D$ does not depend on t and ${\mathsf{\Theta }}_D={\mathsf{\Theta }}_0$. Then,

$$\begin{array}{cc}\hfill \frac{d}{dt}\left(E_{el}\left(u\left(t\right)\right)+E_{th}\left(\mathsf{\Theta }\left(t\right)\right)\right)& =-\frac{{\kappa }_0}{{\mathsf{\Theta }}_0}{\int }_{\mathsf{\Omega }}{\left|\nabla \mathsf{\Theta }\left(t\right)\right|}^{2}dx\le 0.\hfill \end{array}$$

(12)

Proof. 

Since

$$\begin{array}{cc}\hfill \frac{d}{dt}\left(\frac{1}{2}\sigma \left[u\right]:e\left[u\right]\right)& =\sigma \left[u\right]:e\left[\dot{u}\right]\hfill \\ \hfill & =({\sigma }^{*}[u,\mathsf{\Theta }]+\beta (\mathsf{\Theta }-{\mathsf{\Theta }}_0)I):e\left[\dot{u}\right]\hfill \\ \hfill & ={\sigma }^{*}[u,\mathsf{\Theta }]:e\left[\dot{u}\right]+\beta (\mathsf{\Theta }-{\mathsf{\Theta }}_0)\mathrm{div}\dot{u}\hfill \end{array}$$

(13)

we obtain

$$\begin{array}{cc}\hfill \frac{d}{dt}{E}_{el}\left(u\left(t\right)\right)& =\frac{1}{2}{\int }_{\mathsf{\Omega }}\frac{d}{dt}\left(\sigma \left[u\right]:e\left[u\right]\right)dx\hfill \\ \hfill & ={\int }_{\mathsf{\Omega }}{\sigma }^{*}[u,\mathsf{\Theta }]:e\left[\dot{u}\right]dx+{\int }_{\mathsf{\Omega }}\beta (\mathsf{\Theta }-{\mathsf{\Theta }}_0)\left(\mathrm{div}\dot{u}\right)dx\hfill \\ \hfill & ={\int }_{\mathsf{\Omega }}\beta (\mathsf{\Theta }-{\mathsf{\Theta }}_0)\left(\mathrm{div}\dot{u}\right)dx.\hfill \end{array}$$

(14)

Substituting (2) into (14) and using the boundary conditions (5c) and (5d) for $\mathsf{\Theta }$, we obtain

$$\begin{array}{cc}\hfill \frac{d}{dt}{E}_{el}\left(u\left(t\right)\right)& ={\int }_{\mathsf{\Omega }}\left(\frac{1}{{\mathsf{\Theta }}_0}(\mathsf{\Theta }-{\mathsf{\Theta }}_0)\left\{{\kappa }_0\mathsf{\Delta }\mathsf{\Theta }-\chi \frac{\partial \mathsf{\Theta }}{\partial t}\right\}\right)dx\hfill \\ \hfill & =\frac{{\kappa }_0}{{\mathsf{\Theta }}_0}{\int }_{\mathsf{\Gamma }}(\mathsf{\Theta }-{\mathsf{\Theta }}_0)\frac{\partial \mathsf{\Theta }}{\partial n}ds-\frac{{\kappa }_0}{{\mathsf{\Theta }}_0}{\int }_{\mathsf{\Omega }}|\nabla \mathsf{\Theta }{|}^{2}dx-\frac{d}{dt}\left(\frac{\chi }{2{\mathsf{\Theta }}_0}{\int }_{\mathsf{\Omega }}|\mathsf{\Theta }-{\mathsf{\Theta }}_0{|}^{2}dx\right)\hfill \\ \hfill & =-\frac{{\kappa }_0}{{\mathsf{\Theta }}_0}{\int }_{\mathsf{\Omega }}|\nabla \mathsf{\Theta }{|}^{2}dx-\frac{d}{dt}{E}_{th}\left(\mathsf{\Theta }\left(t\right)\right).\hfill \end{array}$$

This gives the energy equality for (5). □

As shown in Proposition 2 and Theorem 1, Biot’s thermoelasticity model is related to both energies $E_{el}\left(u\right)$ and $E_{el}^{*}(u,\mathsf{\Theta })$. We denote their energy densities as follows:

$$\begin{array}{cc}\hfill & W\left(u\right):=\sigma \left[u\right]:e\left[u\right],\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & W^{*}(u,\mathsf{\Theta }):={\sigma }^{*}[u,\mathsf{\Theta }]:e^{*}[u,\mathsf{\Theta }],\hfill \end{array}$$

(16)

where $W\left(u\right)$ and $W^{*}(u,\mathsf{\Theta })$ are the elastic and thermoelastic energy densities, respectively.

2.3. Numerical Experiment

2.3.1. Non-Dimensional Setting

In the following numerical examples, we introduce a non-dimensional form of Biot’s model. We consider the following scaling for x, t, u, C (or $\lambda$, $\mu$), and $\mathsf{\Theta }$:

$$\begin{array}{c}\hfill \tilde{x}=\frac{x}{c_x},\tilde{t}=\frac{t}{c_t},\tilde{u}=\frac{u}{c_u},\tilde{C}=\frac{C}{c_e},\tilde{\mathsf{\Theta }}=\frac{\mathsf{\Theta }-{\mathsf{\Theta }}_0}{c_{\mathsf{\Theta }}},{\tilde{a}}_L=\frac{c_x{c}_{\mathsf{\Theta }}}{c_u}{a}_L,\tilde{\beta }=1,\end{array}$$

(17)

where $c_x$, $c_t$, $c_u$, $c_e$, and $c_{\mathsf{\Theta }}>0$ are the scaling parameters. Let $c_x$ [$\mathrm{m}$], $c_e$ [$\mathrm{Pa}$], and $c_{\mathsf{\Theta }}$ [$\mathrm{K}$] be characteristic scales for the length of the domain, the size of the elastic tensor, and the temperature, respectively. The parameters $c_t$ and $c_u$ are defined as

$$\begin{array}{c}\hfill c_t:=\frac{c_x^2\chi }{{\kappa }_0}\left[\mathrm{s}\right],c_u:=\frac{c_{\mathsf{\Theta }}{c}_x\beta }{c_e}\left[\mathrm{m}\right],\end{array}$$

(18)

where $\chi$ [$\mathrm{Pa}·{\mathrm{K}}^{-1}$], ${\kappa }_0$ [$\mathrm{Pa}·{\mathrm{m}}^2·{\mathrm{s}}^{-1}·{\mathrm{K}}^{-1}$] and $\beta =a_L(d\lambda +2\mu )$ [$\mathrm{Pa}·{\mathrm{K}}^{-1}$]. Then, (1) and (2) are written in the following non-dimensional form:

$$\{\begin{array}{ccc}\widetilde{\mathrm{div}}\tilde{\sigma }\left[\tilde{u}\right]=\tilde{\nabla }\tilde{\mathsf{\Theta }}\hfill & \hfill \mathrm{in}\tilde{\mathsf{\Omega }}\times [0,\tilde{T}],& \hfill \left(19\mathrm{a}\right)\\ {\displaystyle \frac{\partial }{\partial \tilde{t}}}\tilde{\mathsf{\Theta }}=\tilde{\mathsf{\Delta }}\tilde{\mathsf{\Theta }}-\delta {\displaystyle \frac{\partial }{\partial \tilde{t}}}(\widetilde{\mathrm{div}}\tilde{u})\hfill & \hfill \mathrm{in}\tilde{\mathsf{\Omega }}\times (0,\tilde{T}].& \hfill \left(19\mathrm{b}\right)\end{array}$$

The system (19) has only three parameters, $\tilde{\lambda }$, $\tilde{\mu }$, and $\delta$. The parameter $\delta$ is a non-dimensional thermoelastic coupling parameter defined by

$$\begin{array}{ccc}\hfill \delta & :=& \frac{{\mathsf{\Theta }}_0{\beta }^2}{c_e\chi }[-],\hfill \end{array}$$

and $\delta >0$. If we choose $\delta =0$, (19b) is decoupled from (19a), and the temperature field $\tilde{\mathsf{\Theta }}$ in (19a) is essentially a given function. In the following example, the case $\delta =0$ is referred to as the uncoupled case.

Under the above scaling, we denote the (thermo)elastic strain, stress tensors, and (thermo)elastic energy densities as follows:

$$\begin{array}{cc}\tilde{e}\left[\tilde{u}\right]:=\frac{1}{2}\left(\frac{\partial {\tilde{u}}_i}{\partial {\tilde{x}}_j}+\frac{\partial {\tilde{u}}_j}{\partial {\tilde{x}}_i}\right)=\frac{c_x}{c_u}e\left[u\right],\hfill & \hfill \left(20\mathrm{a}\right)\\ \tilde{\sigma }\left[\tilde{u}\right]:=\tilde{C}\tilde{e}\left[\tilde{u}\right]=\frac{c_x}{c_u{c}_e}\sigma \left[u\right],\hfill & \hfill \left(20\mathrm{b}\right)\\ \tilde{W}\left(\tilde{u}\right):=\tilde{\sigma }\left[\tilde{u}\right]:\tilde{e}\left[\tilde{u}\right]=\frac{c_e}{{\left(\beta c_{\mathsf{\Theta }}\right)}^2}W\left[u\right],\hfill & \hfill \left(20\mathrm{c}\right)\\ {\tilde{\sigma }}^{*}[\tilde{u},\tilde{\mathsf{\Theta }}]:=\tilde{\sigma }\left[\tilde{u}\right]-\tilde{\mathsf{\Theta }}I=\frac{1}{\beta c_{\mathsf{\Theta }}}{\sigma }^{*}[u,\mathsf{\Theta }],\hfill & \hfill \left(20\mathrm{d}\right)\\ {\tilde{e}}^{*}[\tilde{u},\tilde{\mathsf{\Theta }}]:=\tilde{e}\left[\tilde{u}\right]-{\tilde{a}}_L\tilde{\mathsf{\Theta }}I=\frac{c_x}{c_u}{\sigma }^{*}[u,\mathsf{\Theta }],\hfill & \hfill \left(20\mathrm{e}\right)\\ {\tilde{W}}^{*}(\tilde{u},\tilde{\mathsf{\Theta }}):={\tilde{\sigma }}^{*}[\tilde{u},\tilde{\mathsf{\Theta }}]:{\tilde{e}}^{*}[\tilde{u},\tilde{\mathsf{\Theta }}]=\frac{c_e}{{\left(\beta c_{\mathsf{\Theta }}\right)}^2}{W}^{*}[u,\mathsf{\Theta }].\hfill & \hfill \left(20\mathrm{f}\right)\end{array}$$

In the following section, we apply these non-dimensional forms and omit ∼ for simplicity.

2.3.2. Numerical Setup and Time Discretization

In the following examples, we set Young’s modulus $E_{\mathrm{Y}}=1$, Poisson’s ratio ${\nu }_{\mathrm{P}}=0.32$, the coefficient of linear thermal expansion $a_L=0.475$ and the thermoelasticity coupling parameter $\delta =0.0,0.1,0.5$ in the non-dimensional form of (19). We consider two numerical examples for (19), an L-shaped cantilever domain and a square domain with a crack (more precisely, a very sharp notch), as illustrated in Figure 2.

We apply the following implicit time discretization for (19):

$$\begin{array}{cc}\hfill & \left\{\begin{array}{cc}{\displaystyle -\mathrm{div}{\sigma }^{*}[u^k,{\mathsf{\Theta }}^{k-1}]=0}\hfill & {\displaystyle \mathrm{in}\mathsf{\Omega },}\hfill \\ {\displaystyle \frac{{\mathsf{\Theta }}^k-{\mathsf{\Theta }}^{k-1}}{\mathsf{\Delta }t}-\mathsf{\Delta }{\mathsf{\Theta }}^k+\delta \mathrm{div}\left(\frac{u^k-u^{k-1}}{\mathsf{\Delta }t}\right)=0}\hfill & {\displaystyle \mathrm{in}\mathsf{\Omega },}\hfill \end{array}\right.\hfill \end{array}$$

(21)

where $u^k$ and ${\mathsf{\Theta }}^k$ are approximations to u and $\mathsf{\Theta }$ at $t=k\mathsf{\Delta }t(k=0,1,2,\dots )$. At each time step $k=1,2,\dots$, we solve (21) with given boundary and initial conditions (5) using the finite element method. The details of the weak forms for (21) and their unique solvability are described in Appendix A.

In observation area $\mathcal{A}$ illustrated in Figure 2, we define the average of (thermo)elastic energy densities in $\mathcal{A}$ as follows:

$$\begin{array}{cc}\hfill \mathcal{W}\left(\mathcal{A}\right)& :=\frac{1}{\left|\mathcal{A}\right|}{\int }_{\mathcal{A}}W\left(u\right)dx,\hfill \\ \hfill {\mathcal{W}}^{*}\left(\mathcal{A}\right)& :=\frac{1}{\left|\mathcal{A}\right|}{\int }_{\mathcal{A}}{W}^{*}(u,\mathsf{\Theta })dx,\hfill \end{array}$$

and the differences between $\mathcal{W}\left(\mathcal{A}\right)$ and ${\mathcal{W}}^{*}\left(\mathcal{A}\right)$ for each $\delta >0$ and for $\delta =0$ are defined by

$$\begin{array}{c}\hfill \mathsf{\Delta }\mathcal{W}\left(\mathcal{A}\right):=\mathcal{W}\left(\mathcal{A}\right){|}_{\delta }-\mathcal{W}\left(\mathcal{A}\right){|}_{\delta =0},\\ \mathsf{\Delta }{\mathcal{W}}^{*}\left(\mathcal{A}\right):={\mathcal{W}}^{*}\left(\mathcal{A}\right){|}_{\delta }-{\mathcal{W}}^{*}\left(\mathcal{A}\right){|}_{\delta =0}.\hfill \end{array}$$

In the following examples, we use the software FreeFEM [39] with P2 elements and unstructured meshes. For the time interval and time step, we use $0\le t\le 0.1$ and $\mathsf{\Delta }t=1\times {10}^{-4}$, respectively.

2.3.3. L-Shape Cantilever

Here, we consider the L-shaped cantilever whose left side is fixed, and the vertical displacement $u_2$ is given on the right side, as illustrated in Figure 2 left. We denote the left and right boundaries by ${\mathsf{\Gamma }}_D^u$ and ${\mathsf{\Gamma }}_{DN}^u$, respectively, and define ${\mathsf{\Gamma }}_N^u:=\mathsf{\Gamma }\backslash ({\mathsf{\Gamma }}_D^u\cup {\mathsf{\Gamma }}_{DN}^u)$. The boundary conditions for u are

$$\begin{array}{cc}\hfill & u=0,\mathrm{on}{\mathsf{\Gamma }}_D^u,\hspace{1em}\left\{\begin{array}{c}{\sigma }_{11}^{*}[u,\mathsf{\Theta }]n=0,\hfill \\ u_2=-0.1t\hfill \end{array}\right.\mathrm{on}{\mathsf{\Gamma }}_{DN}^u,\hspace{1em}{\sigma }^{*}[u,\mathsf{\Theta }]n=0\mathrm{on}{\mathsf{\Gamma }}_N^u.\hfill \end{array}$$

For $\mathsf{\Theta }$, we suppose $\frac{\partial \mathsf{\Theta }}{\partial n}=0$ on $\mathsf{\Gamma }$ and the initial temperature ${\mathsf{\Theta }}_{*}=0$. Although we adopt the above slightly modified boundary conditions in this example, the previous arguments are valid with small modifications, and we omit their details.

We apply the finite element method to (21). The total number of triangular meshes = 18,215 and the number of nodes (the vertices of the triangles) $=9301$.

As shown in the lower part of Figure 3, we observe that the highest temperature is in the contracting area and the lowest is in the expanding area. Furthermore, there exists a contribution $\delta$ for each $\delta >0$ during the loading process. Although the disparity is small, the thermoelastic coupling parameter $\delta$ contributes to the variations in $\mathcal{W}\left(u\right)$ and ${\mathcal{W}}^{*}\left(u\right)$, as shown in Figure 4a,b. Here, a larger $\delta$ value implies larger $\mathcal{W}\left(\mathcal{A}\right)$ and ${\mathcal{W}}^{*}\left(\mathcal{A}\right)$ values (Figure 4d,e). In addition, we also observe that ${\mathcal{W}}^{*}\left(\mathcal{A}\right)$ is larger than $\mathcal{W}\left(\mathcal{A}\right)$ for each $\delta >0$ (Figure 4c).

In the L-shape cantilever case for each $\delta >0$, we conclude that the thermal coupling parameter enhances the singularity of (thermo)elastic energy in the expanding area. The (thermo)elastic energy plays a role in the driving force in the phase field model [23], which means that the parameter $\delta$ can accelerate crack growth in the expanding area.

2.3.4. Cracked Domain

Here, we consider a cracked domain with vertical displacements on the top and bottom sides, and the other sides are free traction, as shown in Figure 2 right. The boundary conditions for u are

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{u}_1=0,\hfill \\ u_2=\pm t\hfill \end{array}\right.\mathrm{on}{\mathsf{\Gamma }}_{\pm D}^u,\hspace{1em}{\sigma }^{*}[u,\mathsf{\Theta }]n=0\mathrm{on}{\mathsf{\Gamma }}_N^u,\end{array}$$

where ${\mathsf{\Gamma }}_{+D}^u$ and ${\mathsf{\Gamma }}_{-D}^u$ denote the top and bottom boundaries of $\mathsf{\Omega }$, respectively, and ${\mathsf{\Gamma }}_N^u:=\mathsf{\Gamma }\backslash ({\mathsf{\Gamma }}_{+D}^u\cup {\mathsf{\Gamma }}_{-D}^u)$. For $\mathsf{\Theta }$, we suppose $\frac{\partial \mathsf{\Theta }}{\partial n}=0$ on ${\mathsf{\Gamma }}_N^{\mathsf{\Theta }}=\mathsf{\Gamma }$ and the initial temperature ${\mathsf{\Theta }}_{*}=0$.

We use the finite element method to solve (21). Therefore, the total number of triangular meshes and the number of nodes (the vertices of the triangles) are 11,176 and 5722, respectively.

From Figure 5 left, we conclude that the area that expands the most (i.e., $\mathrm{div}u$ is largest) appears near the crack tip. This can be compared with the analytical solution for the linear elasticity in a cracked domain in Appendix B. We also observe that the region with the lowest temperature appears to the right of the crack tip in Figure 5 right. From the temporal change in the temperature along the $x_1$ axis plotted in Figure 6 right, we also observe that the lowest temperature region appears in $0.5<x_1<0.6$ and that the temperature decreases over time. This is shown in Figure 6 left, where the value of $\mathrm{div}u$ is plotted along the $x_1$ axis and $\mathrm{div}u$ is increasing over time; i.e., the heat source term $\mathrm{div}\dot{u}$ in (2) is positive. Experimentally, the lowest temperature around the crack tips does not match with the studies of of Zehnder et al. [40], Rusinek et al. [41], and Wang et al. [42]. They record that the highest temperature occurs around the crack tips, which result from a plastic zone around the crack tips. In the present study, we do not consider the plastic zone. However, it would be more possible to use a thermo-viscous-elasticity condition. We will consider and study it using the thermo-viscous-elasticity equation in the future work [43].

Similar to Section 2.3.3, for each $\delta >0$, we obtain variations of $\mathcal{W}\left(\mathcal{A}\right)$ and ${\mathcal{W}}^{*}\left(\mathcal{A}\right)$ in subdomain $\mathcal{A}$ (Figure 7), where the subdomain $\mathcal{A}$ corresponds to the area that expands the most. From Figure 7, it is observed that ${\mathcal{W}}^{*}\left(\mathcal{A}\right)$ is larger than $\mathcal{W}\left(\mathcal{A}\right)$. This suggests that the thermoelastic energy density $W^{*}(u,\mathsf{\Theta })$ has a higher value than the elastic energy density $W\left(u\right)$. These observations are confirmed by the comparison of our thermal fracturing phase field models.

3. Crack Propagation under Thermal Stress

This section is devoted to the phase field models for thermal fracturing, which are the main purpose of this paper.

3.1. Fracturing Phase Field Model (F-PFM)

According to the works [22,26], we introduce fracturing PFM (we call it F-PFM) in this section. Let $\mathsf{\Omega }$ be a bounded (uncracked) domain in ${\mathbb{R}}^d$ and $\mathsf{\Gamma }:=\partial \mathsf{\Omega }={\mathsf{\Gamma }}_D^u\cup {\mathsf{\Gamma }}_N^u$, similar to Section 2. In F-PFM, a crack in $\mathsf{\Omega }$ at time t is described by a damage variable $z(x,t)\in [0,1]$ for $x\in \overline{\mathsf{\Omega }}$ with space regularization. The cracked and uncracked regions are represented by $z\approx 1$ and $z\approx 0$, respectively, and $z\in (0,1)$ indicates slight damage. A typical example of a straight crack in a square domain is illustrated in Figure 8.

The F-PFM is described as:

$$\{\begin{array}{ccc}-\mathrm{div}\left({(1-z)}^2\sigma \left[u\right]\right)=0\hfill & \mathrm{in}\mathsf{\Omega }\times [0,T],\hfill & \hfill \left(22\mathrm{a}\right)\\ \alpha \frac{\partial z}{\partial t}={\left(ϵ\mathrm{div}\left({\gamma }_{*}\nabla z\right)-\frac{{\gamma }_{*}}{ϵ}z+(1-z)W\left(u\right)\right)}_{+}\hfill & \hfill \mathrm{in}\mathsf{\Omega }\times [0,T],& \hfill \left(22\mathrm{b}\right)\end{array}$$

with the following boundary and initial conditions:

$$\{\begin{array}{ccc}u=u_D(x,t)\hfill & \mathrm{on}{\mathsf{\Gamma }}_D^u\times [0,T],\hfill & \hfill \left(23\mathrm{a}\right)\\ \sigma \left[u\right]n=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_N^u\times [0,T],\hfill & \hfill \left(23\mathrm{b}\right)\\ \frac{\partial z}{\partial n}=0\hfill & \mathrm{on}\mathsf{\Gamma }\times [0,T],\hfill & \hfill \left(23\mathrm{c}\right)\\ z\left(x,0\right)=z_{*}\left(x\right)\hfill & \mathrm{in}\mathsf{\Omega },\hfill & \hfill \left(23\mathrm{d}\right)\end{array}$$

where the displacement $u:\overline{\mathsf{\Omega }}\times [0,T]\mapsto {\mathbb{R}}^d$ and the damage variable $z:\overline{\mathsf{\Omega }}\times [0,T]\mapsto [0,1]$ are unknowns. The parameters $\alpha >0$ and $ϵ>0$ are small numbers related to regularization in time and space, respectively. The critical energy release rate is denoted by ${\gamma }_{*}$ (which is often denoted by ${\mathrm{G}}_{\mathrm{c}}$), and the elastic energy density is defined by $W=W\left(u\right):=\sigma \left[u\right]:e\left[u\right]$. In (22b), the term W works as a driving force for z.

The symbol ${( )}_{+}$ on the right-hand side in (22b) denoted the positive part ${\left(s\right)}_{+}:=\max (s,0)$, and it represents the irreversible property of crack growth.

F-PFM is derived as a unidirectional gradient flow of the total energy ${\mathcal{E}}_{el}(u,z)+E_s\left(z\right)$, where

$$\begin{array}{cc}\hfill {\mathcal{E}}_{el}(u,z)& :=\frac{1}{2}{\int }_{\mathsf{\Omega }}{(1-z)}^2\sigma \left[u\right]:e\left[u\right]dx,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill E_s\left(z\right)& :=\frac{1}{2}{\int }_{\mathsf{\Omega }}{\gamma }_{*}\left(ϵ|\nabla z{|}^2+\frac{{|z|}^2}{ϵ}\right)dx.\hfill \end{array}$$

(25)

More precisely, $u\left(t\right)$ obeys the following variational principle:

$$u\left(t\right)=\underset{u\in V\left(u_D\left(t\right)\right)}{argmin}{\mathcal{E}}_{el}(u,z\left(t\right)),$$

(26)

and (22b) becomes a gradient flow of the energy $\underset{u}{min}{\mathcal{E}}_{el}(u,z)+E_s\left(z\right)$.

We remark that ${\mathcal{E}}_{el}(u,z)$ is a modified elastic energy, which corresponds to the elastic energy with a damaged Young’s modulus ${\tilde{E}}_{\mathrm{Y}}={(1-z)}^2{E}_{\mathrm{Y}}$. The energy $E_s\left(z\right)$ is regularized surface energy, which approximates the crack area ($d=3$) or length ($d=2$) as $ϵ\mapsto 0$. Please see [26] for more details. The following energy equality for F-PFM is shown in [26] ([22] for the antiplane setting).

Theorem 2

(Energy equality for F-PFM). Let $\left(u\right(x,t),z(x,t\left)\right)$ be a sufficiently smooth solution to (22) and (23). If $u_D$ is independent of t, then we have

$$\begin{array}{c}\hfill \frac{d}{dt}\left({\mathcal{E}}_{el}(u\left(t\right),z\left(t\right))+E_s\left(z\left(t\right)\right)\right)=-\alpha {\int }_{\mathsf{\Omega }}{\left|\dot{z}\right|}^{2}dx\le 0.\end{array}$$

(27)

Proof. 

Differentiating the total energy in t and applying integration by parts, we obtain

$$\begin{array}{cc}\frac{d}{dt}\hfill & \left({\mathcal{E}}_{el}(u\left(t\right),z\left(t\right))+E_s\left(z\left(t\right)\right)\right)\hfill \\ & ={\int }_{\mathsf{\Omega }}{(1-z)}^2\sigma \left[u\right]:e\left[\dot{u}\right]dx+{\int }_{\mathsf{\Omega }}\left({\gamma }_{*}ϵ\nabla z·\nabla \dot{z}+\left(\frac{{\gamma }_{*}}{ϵ}z-(1-z)W\left(u\right)\right)\dot{z}\right)dx\hfill \\ & ={\int }_{\mathsf{\Gamma }}{(1-z)}^2\underset{0}{\underbrace{\left(\sigma \left[u\right]n\right)}}·\dot{u} ds-{\int }_{\mathsf{\Omega }}\underset{0}{\underbrace{\mathrm{div}\left({(1-z)}^2\sigma \left[u\right]\right)}}·\dot{u} dx+{\int }_{\mathsf{\Gamma }}{\gamma }_{*}ϵ\underset{0}{\underbrace{\frac{\partial z}{\partial n}}}\dot{z} ds-{\int }_{\mathsf{\Omega }}\mathcal{H}\dot{z} dx,\hfill \end{array}$$

(28)

where we define $\mathcal{H}:=ϵ\mathrm{div}\left({\gamma }_{*}\nabla z\right)-\frac{{\gamma }_{*}}{ϵ}z+(1-z)W\left(u\right)$. Since (22b) is written as $\alpha \dot{z}={\left(\mathcal{H}\right)}_{+}$, using the equality $\mathcal{H}{\left(\mathcal{H}\right)}_{+}={\left(\mathcal{H}\right)}_{+}^2$, we conclude that

$$\begin{array}{cc}\hfill & \frac{d}{dt}\left({\mathcal{E}}_{el}(u\left(t\right),z\left(t\right))+E_s\left(z\left(t\right)\right)\right)\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=-{\int }_{\mathsf{\Omega }}\mathcal{H}\dot{z} dx=-{\int }_{\mathsf{\Omega }}\mathcal{H}\frac{{\left(\mathcal{H}\right)}_{+}}{\alpha } dx=-{\int }_{\mathsf{\Omega }}\frac{{\left(\mathcal{H}\right)}_{+}^2}{\alpha } dx=-{\int }_{\mathsf{\Omega }}\alpha {\left|\dot{z}\right|}^{2}dx.\hfill \end{array}$$

 □

3.2. Thermal Fracturing Phase Field Model 1 (TF-PFM1)

To combine the Biot model in (1) and (2) and F-PFM in (22), their variational principles for u, Proposition 2, and (26) suggest that we consider the following modified thermoelastic energy:

$$\begin{array}{ccc}\hfill {\mathcal{E}}_{el}^{*}(u,\mathsf{\Theta },z)& :=& \frac{1}{2}{\int }_{\mathsf{\Omega }}{(1-z)}^2{\sigma }^{*}[u,\mathsf{\Theta }]:e^{*}[u,\mathsf{\Theta }]dx,\hfill \end{array}$$

(29)

and a variational principle:

$$u\left(t\right)=\underset{u\in V\left(u_D\left(t\right)\right)}{argmin}{\mathcal{E}}_{el}^{*}(u,\mathsf{\Theta }\left(t\right),z\left(t\right)).$$

(30)

From the definition of the modified thermoelastic energy (29), it is natural to replace the driving force term $W\left(u\right)=\sigma \left[u\right]:e\left[u\right]$ in (22b) by the thermoelastic energy density $W^{*}(u,\mathsf{\Theta }):={\sigma }^{*}[u,\mathsf{\Theta }]:e^{*}[u,\mathsf{\Theta }]$.

For heat Equation (2), since $\beta =a_L(d\lambda +2\mu )$ and Lamè’s constants ($\lambda ,\mu$) are replaced by damaged constants (${(1-z)}^2\lambda$, ${(1-z)}^2\mu$), $\beta$ should also be replaced by damaged constant ${(1-z)}^2\beta$. The thermal conductivity ${\kappa }_0$ is also considered to be modified by z because the heat is usually insulated across the crack. We suppose $\kappa =\kappa \left(z\right)>0$ in this section, and we set it as $\kappa \left(z\right)={(1-z)}^2{\kappa }_0$ in Section 4.

Summarizing the above statements, we obtain the following thermal fracturing model, PFM 1 (TF-PFM1):

$$\{\begin{array}{ccc}-\mathrm{div}\left({(1-z)}^2{\sigma }^{*}[u,\mathsf{\Theta }]\right)=0\hfill & \mathrm{in}\mathsf{\Omega }\times [0,T],\hfill & \hfill \left(31\mathrm{a}\right)\\ \alpha \frac{\partial z}{\partial t}={\left(ϵ\mathrm{div}({\gamma }_{*}\nabla z)-\frac{{\gamma }_{*}}{ϵ}z+(1-z)W^{*}(u,\mathsf{\Theta })\right)}_{+}\hfill & \mathrm{in}\mathsf{\Omega }\times [0,T],\hfill & \hfill \left(31\mathrm{b}\right)\\ \chi \frac{\partial \mathsf{\Theta }}{\partial t}=\mathrm{div}\left(\kappa \left(z\right)\nabla \mathsf{\Theta }\right)-{\mathsf{\Theta }}_0{(1-z)}^2\beta \frac{\partial }{\partial t}\left(\mathrm{div}u\right)\hfill & \mathrm{in}\mathsf{\Omega }\times (0,T],\hfill & \hfill \left(31\mathrm{c}\right)\end{array}$$

Similar to (1), (2) and (22), the boundary and the initial conditions to solve (31) are presented as follows:

$$\{\begin{array}{ccc}u=u_D(x,t)\hfill & \mathrm{on}{\mathsf{\Gamma }}_D^u\times [0,T],\hfill & \hfill \left(32\mathrm{a}\right)\\ {\sigma }^{*}[u,\mathsf{\Theta }]n=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_N^u\times [0,T],\hfill & \hfill \left(32\mathrm{b}\right)\\ \mathsf{\Theta }={\mathsf{\Theta }}_D(x,t)\hfill & \mathrm{on}{\mathsf{\Gamma }}_D^{\mathsf{\Theta }}\times [0,T],\hfill & \hfill \left(32\mathrm{c}\right)\\ \frac{\partial \mathsf{\Theta }}{\partial n}=0\hfill & \mathrm{on}{\mathsf{\Gamma }}_N^{\mathsf{\Theta }}\times [0,T],\hfill & \hfill \left(32\mathrm{d}\right)\\ \frac{\partial z}{\partial n}=0\hfill & \mathrm{on}\mathsf{\Gamma }\times [0,T],\hfill & \hfill \left(32\mathrm{e}\right)\\ z(x,0)=z_{*}\left(x\right)\hfill & \mathrm{in}\mathsf{\Omega },\hfill & \hfill \left(32\mathrm{f}\right)\\ \mathsf{\Theta }(x,0)={\mathsf{\Theta }}_{*}\left(x\right)\hfill & \mathrm{in}\mathsf{\Omega }.\hfill & \hfill \left(32\mathrm{g}\right)\end{array}$$

In the following, for simplicity, we define

$$\begin{array}{c}\hfill {\sigma }_z^{*}[u,\mathsf{\Theta }]:={(1-z)}^2{\sigma }^{*}[u,\mathsf{\Theta }].\end{array}$$

As a natural extension of Proposition 2 and Theorem 1, we obtain the following “partial” energy equality for TF-PFM1.

Theorem 3

(Energy equality for TF-PFM1). We suppose that $u_D\in H^{\frac{1}{2}}({\mathsf{\Gamma }}_D^u;{\mathbb{R}}^2)$ and $\mathsf{\Theta }\in L^2\left(\mathsf{\Theta }\right)$ are given and do not depend on t. If $u(x,t)$ and $z(u,t)$ are sufficiency smooth and satisfy (31a), (31b), (32a), (32b), (32e), and (32f), the following energy equality holds:

$$\begin{array}{c}\hfill \frac{d}{dt}\left({\mathcal{E}}_{el}^{*}(u\left(t\right),\mathsf{\Theta },z\left(t\right))+E_s\left(z\left(t\right)\right)\right)=-\alpha {\int }_{\mathsf{\Omega }}{\left|\dot{z}\right|}^{2}dx\le 0.\end{array}$$

(33)

Proof. 

Under this condition, let us derive ${\mathcal{E}}_{el}^{*}(u\left(t\right),\mathsf{\Theta },z\left(t\right))$ and $E_s\left(z\left(t\right)\right)$ with respect to t.

$$\begin{array}{cc}\frac{d}{dt}\hfill & \left({\mathcal{E}}_{el}^{*}(u\left(t\right),\mathsf{\Theta },z\left(t\right))+E_s\left(z\left(t\right)\right)\right)\hfill \\ & =\frac{1}{2}\frac{d}{dt}{\int }_{\mathsf{\Omega }}\left({\sigma }_z^{*}[u,\mathsf{\Theta }]:e^{*}[u,\mathsf{\Theta }]\right)dx+\frac{1}{2}\frac{d}{dt}{\int }_{\mathsf{\Omega }}{\gamma }_{*}\left(ϵ|\nabla z{|}^2+\frac{{\left|z\right|}^2}{ϵ}\right)dx\hfill \\ & ={\int }_{\mathsf{\Omega }}{\sigma }_z^{*}[u,\mathsf{\Theta }]:e\left[\dot{u}\right]dx+{\int }_{\mathsf{\Omega }}\left({\gamma }_{*}ϵ\nabla z·\nabla \dot{z}+\left(\frac{{\gamma }_{*}}{ϵ}z-(1-z)W^{*}(u,\mathsf{\Theta })\right)\dot{z}\right)dx\hfill \\ & ={\int }_{\mathsf{\Gamma }}\underset{0}{\underbrace{{\sigma }_z^{*}[u,\mathsf{\Theta }]n}}·e\left[\dot{u}\right]ds-{\int }_{\mathsf{\Omega }}\underset{0}{\underbrace{\mathrm{div}{\sigma }_z^{*}[u,\mathsf{\Theta }]}}·e\left[\dot{u}\right]dx+{\gamma }_{*}ϵ{\int }_{\mathsf{\Gamma }}\underset{0}{\underbrace{\frac{\partial z}{\partial n}}}\dot{z} ds-{\int }_{\mathsf{\Omega }}{\mathcal{H}}^{*}\dot{z} dx,\hfill \end{array}$$

(34)

where we also define ${\mathcal{H}}^{*}:=ϵ\mathrm{div}({\gamma }_{*}\nabla z)-\frac{{\gamma }_{*}}{ϵ}z+(1-z)W^{*}(u,\mathsf{\Theta })$. Since (31b) is changed to $\alpha \dot{z}={\left({\mathcal{H}}^{*}\right)}_{+}$, similar to that in Section 3.1, we conclude that

$$\begin{array}{cc}\hfill & \frac{d}{dt}\left({\mathcal{E}}_{el}^{*}(u\left(t\right),\mathsf{\Theta },z\left(t\right))+E_s\left(z\left(t\right)\right)\right)=-\alpha {\int }_{\mathsf{\Omega }}{\left|\dot{z}\right|}^{2}dx\le 0,\hfill \end{array}$$

which is equivalent to (33). □

3.3. Thermal Fracturing Phase Field Model 2 (TF-PFM2)

In the previous section, we proposed TF-PFM1 based on the thermoelastic energy $E_{el}^{*}(u,\mathsf{\Theta })$. We proved a variational principle but proved only partial energy equality. As shown in Section 2.2, the Biot model is related to both energies $E_{el}^{*}(u,\mathsf{\Theta })$ and $E_{el}\left(u\right)$. The variational principle holds for $E_{el}^{*}(u,\mathsf{\Theta })$ (Proposition 2), and the energy equality holds for $E_{el}\left(u\right)$ (Theorem 1). This motivates us to consider another type of thermal fracturing PFM based on elastic energy $E_{el}\left(u\right)$. We call the following thermal fracturing model TF-PFM2:

$$\{\begin{array}{ccc}-\mathrm{div}\left({(1-z)}^2{\sigma }^{*}[u,\mathsf{\Theta }]\right)=0\hfill & \mathrm{in}\mathsf{\Omega }\times [0,T],\hfill & \hfill \left(35\mathrm{a}\right)\\ \alpha \frac{\partial z}{\partial t}={\left(ϵ\mathrm{div}({\gamma }_{*}\nabla z)-\frac{{\gamma }_{*}}{ϵ}z+(1-z)W\left(u\right)\right)}_{+}\hfill & \mathrm{in}\mathsf{\Omega }\times [0,T],\hfill & \hfill \left(35\mathrm{b}\right)\\ \chi \frac{\partial \mathsf{\Theta }}{\partial t}=\mathrm{div}\left(\kappa \left(z\right)\nabla \mathsf{\Theta }\right)-{\mathsf{\Theta }}_0{(1-z)}^2\beta \frac{\partial }{\partial t}\left(\mathrm{div}u\right)\hfill & \mathrm{in}\mathsf{\Omega }\times (0,T].\hfill & \hfill \left(35\mathrm{c}\right)\end{array}$$

The associated boundary and initial conditions are given by (32). For this model, we can show the following energy equality.

Theorem 4

(Energy equality for TF-PFM2). We suppose that $\left(u\right(x,t),$$\mathsf{\Theta }(x,t),z(x,t))$ is a sufficiently smooth solution for (35) and (32). If $u_D$ is independent of t and ${\mathsf{\Theta }}_D={\mathsf{\Theta }}_0$, then the following energy equality holds:

$$\begin{array}{cc}\hfill & \frac{d}{dt}\left({\mathcal{E}}_{el}(u\left(t\right),z\left(t\right))+E_s\left(z\left(t\right)\right)+E_{th}\left(\mathsf{\Theta }\left(t\right)\right)\right)=-\frac{1}{{\mathsf{\Theta }}_0}{\int }_{\mathsf{\Omega }}\kappa \left(z\right){\left|\nabla \mathsf{\Theta }\right|}^{2}dx-\alpha {\int }_{\mathsf{\Omega }}{\left|\dot{z}\right|}^{2}dx\le 0.\hfill \end{array}$$

(36)

Proof. 

Since the relation in (13) is written as

$$\begin{array}{cc}\hfill & \frac{d}{dt}\left(\frac{1}{2}W\left(u\right)\right)={\sigma }^{*}[u,\mathsf{\Theta }]:e\left[\dot{u}\right]+\beta (\mathsf{\Theta }-{\mathsf{\Theta }}_0)\mathrm{div}\dot{u},\hfill \end{array}$$

we obtain

$$\begin{array}{cc}\hfill & \frac{d}{dt}\left(\frac{1}{2}{(1-z)}^{2}W\left(u\right)\right)={\sigma }_z^{*}[u,\mathsf{\Theta }]:e\left[\dot{u}\right]+\beta {(1-z)}^2(\mathsf{\Theta }-{\mathsf{\Theta }}_0)\mathrm{div}\dot{u}-(1-z)\dot{z}W\left(u\right).\hfill \end{array}$$

Hence, we have

$$\begin{array}{cc}\frac{d}{dt}{\mathcal{E}}_{el}(u\left(t\right),z\left(t\right))\hfill & +\frac{d}{dt}{E}_s\left(z\left(t\right)\right)\hfill \\ & ={\int }_{\mathsf{\Omega }}\frac{d}{dt}\left(\frac{1}{2}{(1-z)}^{2}W\left(u\right)\right)dx+{\int }_{\mathsf{\Omega }}\left(ϵ\mathrm{div}({\gamma }_{*}\nabla z)-\frac{{\gamma }_{*}}{ϵ}z\right)\dot{z} dx\hfill \\ & =\underset{0}{\underbrace{{\int }_{\mathsf{\Omega }}{\sigma }_z^{*}[u,\mathsf{\Theta }]:e\left[\dot{u}\right]dx}}+{\int }_{\mathsf{\Omega }}\beta {(1-z)}^2(\mathsf{\Theta }-{\mathsf{\Theta }}_0)\mathrm{div}\dot{u} dx-{\int }_{\mathsf{\Omega }}\mathcal{H}\dot{z} dx\hfill \\ & ={\int }_{\mathsf{\Omega }}\beta {(1-z)}^2(\mathsf{\Theta }-{\mathsf{\Theta }}_0)\mathrm{div}\dot{u} dx-{\int }_{\mathsf{\Omega }}\alpha {\left|\dot{z}\right|}^2 dx,\hfill \end{array}$$

(37)

where $\mathcal{H}=ϵ\mathrm{div}({\gamma }_{*}\nabla z)-\frac{{\gamma }_{*}}{ϵ}z+(1-z)W\left(u\right)$.

On the other hand,

$$\begin{array}{cc}\hfill \frac{d}{dt}{E}_{th}\left(\mathsf{\Theta }\left(t\right)\right)& =\frac{\chi }{{\mathsf{\Theta }}_0}{\int }_{\mathsf{\Omega }}(\mathsf{\Theta }-{\mathsf{\Theta }}_0)\dot{\mathsf{\Theta }} dx\hfill \\ \hfill & =\frac{1}{{\mathsf{\Theta }}_0}{\int }_{\mathsf{\Omega }}(\mathsf{\Theta }-{\mathsf{\Theta }}_0)\left\{\mathrm{div}(\kappa \left(z\right)\nabla \mathsf{\Theta })-{\mathsf{\Theta }}_0\beta {(1-z)}^2\mathrm{div}\dot{u}\right\} dx\hfill \\ \hfill & =-\frac{1}{{\mathsf{\Theta }}_0}{\int }_{\mathsf{\Omega }}\kappa \left(z\right){|\nabla \mathsf{\Theta }|}^2 dx-{\int }_{\mathsf{\Omega }}\beta {(1-z)}^2(\mathsf{\Theta }-{\mathsf{\Theta }}_0)\mathrm{div}\dot{u} dx.\hfill \end{array}$$

(38)

Taking a sum of these equalities (37) and (38), then we obtain the energy equality (36). □

4. Numerical Experiments

In this section, we conduct numerical experiments to test F-PFM, TF-PFM1, and TF-PFM2, which were derived in Section 3, and report the numerical results. Through the numerical experiments, we observe the effect of thermal coupling on the crack speed and the crack path during its growth process.

4.1. Non-Dimensional Setting

In the following numerical examples, we suppose $\kappa \left(z\right)={(1-z)}^2{\kappa }_0$. For convenience, we consider the non-dimensional form with (17), (18), (20) and

$$\widetilde{ϵ}=\frac{ϵ}{c_x},\hspace{1em}{\tilde{\gamma }}_{*}=\frac{c_e{\gamma }_{*}}{c_x{\left(\beta c_{\mathsf{\Theta }}\right)}^2},\hspace{1em}\tilde{\alpha }=\frac{c_e\alpha }{c_t{\left(\beta c_{\mathsf{\Theta }}\right)}^2},\hspace{1em}{\tilde{a}}_L=\frac{c_x{c}_{\mathsf{\Theta }}}{c_u}{a}_L,\hspace{1em}\tilde{\beta }=1.$$

Then, TF-PFM1 in (31) is expressed in the following non-dimensional form:

$$\{\begin{array}{ccc}\mathrm{div}({(1-z)}^2\sigma \left[u\right])={(1-z)}^2\nabla \mathsf{\Theta }\hfill & \mathrm{in}\mathsf{\Omega }\times [0,T],\hfill & \hfill \left(39\mathrm{a}\right)\\ \alpha \frac{\partial z}{\partial t}={\left(ϵ\mathrm{div}({\gamma }_{*}\nabla z)-\frac{{\gamma }_{*}}{ϵ}z+(1-z)W^{*}(u,\mathsf{\Theta })\right)}_{+}\hfill & \mathrm{in}\mathsf{\Omega }\times [0,T],\hfill & \hfill \left(39\mathrm{b}\right)\\ \frac{\partial \mathsf{\Theta }}{\partial t}=\mathrm{div}\left({(1-z)}^2\nabla \mathsf{\Theta }\right)-{(1-z)}^2\delta \frac{\partial }{\partial t}\left(\mathrm{div}u\right)\hfill & \mathrm{in}\mathsf{\Omega }\times (0,T].\hfill & \hfill \left(39\mathrm{c}\right)\end{array}$$

For TF-PFM2, we change (39b) to:

$$\alpha \frac{\partial z}{\partial t}={\left(ϵ\mathrm{div}({\gamma }_{*}\nabla z)-\frac{{\gamma }_{*}}{ϵ}z+(1-z)W\left(u\right)\right)}_{+}\hspace{1em}\hspace{1em}\mathrm{in}\mathsf{\Omega }\times [0,T].$$

(40)

4.2. Time Discretization

To solve problem (39), we adopt the following semi-implicit time discretization scheme [22,26]:

$$\{\begin{array}{cc}-\mathrm{div}\left({(1-z^{k-1})}^2{\sigma }^{*}[u^k,{\mathsf{\Theta }}^{k-1}]\right)=0,\hfill & \hfill \left(41\mathrm{a}\right)\\ \alpha \frac{{\tilde{z}}^k-z^{k-1}}{\mathsf{\Delta }t}=ϵ\mathrm{div}\left({\gamma }_{*}\nabla {\tilde{z}}^k\right)-\frac{{\gamma }_{*}}{ϵ}{\tilde{z}}^k+\left(1-{\tilde{z}}^k\right)W^{*}(u^{k-1},{\mathsf{\Theta }}^{k-1}),\hfill & \hfill \left(41\mathrm{b}\right)\\ z^k:=\max \left({\tilde{z}}^k,z^{k-1}\right),\hfill & \hfill \left(41\mathrm{c}\right)\\ {\displaystyle \frac{{\mathsf{\Theta }}^k-{\mathsf{\Theta }}^{k-1}}{\mathsf{\Delta }t}}=\mathrm{div}\left((1-z^{k-1})\nabla {\mathsf{\Theta }}^k\right)-(1-z^{k-1})\delta \mathrm{div}\left({\displaystyle \frac{u^k-u^{k-1}}{\mathsf{\Delta }t}}\right).\hfill & \hfill \left(41\mathrm{d}\right)\end{array}$$

For TF-PFM2, Ref. (41b) is replaced by

$$\begin{array}{cc}\hfill & \alpha \frac{{\tilde{z}}^k-z^{k-1}}{\mathsf{\Delta }t}=ϵ\mathrm{div}\left({\gamma }_{*}\nabla {\tilde{z}}^k\right)-\frac{{\gamma }_{*}}{ϵ}{\tilde{z}}^k+\left(1-{\tilde{z}}^k\right)W\left(u^{k-1}\right),\hfill \end{array}$$

(42)

where $u^k$, $z^k$, and ${\mathsf{\Theta }}^k$ are the approximations of u, z, $\mathsf{\Theta }$, respectively, at time $t_k:=k\mathsf{\Delta }t(k=1,2,3,\dots )$. Since the adaptive mesh technique in the FEM is often effective and accurate in numerical experiments with phase field models, problems (41) and (42) are calculated using adaptive finite elements with P2 elements with a minimum mesh size of $h_{min}=2\times {10}^{-3}$ and a maximum mesh size of $h_{max}=0.1$. The adaptive mesh control at each time step is performed by the adaptmesh() command in FreeFEM based on the variable z. An example of the adaptive mesh is illustrated in Figure 9 right. In addition, the code for the following numerical experiments in the current study is written on FreeFEM [39] and executed on a desktop with an Intel(R) Core i7−7820X [email protected] GHz, 16 core processor, and 64 GB RAM.

4.3. Thermoelastic Effect on the Crack Speed

We set a square domain $\mathsf{\Omega }:={(-1,1)}^2\subset {\mathbb{R}}^2$ with the initial crack $z_{*}\left(x\right):=e^{(-{(x_2/\eta )}^2)}$$/(1+e^{(x_1/\eta )})$ and $\eta =1.5\times {10}^{-2}$. The initial mesh is adapted to $z_{*}\left(x\right)$, as illustrated in Figure 9 right. The material constants for the following examples in the non-dimensional form are listed in

Table 4. List of the non-dimensional parameters for Section 4.3 and Section 4.4.

Parameter	$E_Y$	${\nu }_P$	${\kappa }_0$	$a_L$	$\alpha$	$ϵ$	${\gamma }_{*}$	${\mathsf{\Theta }}_{*}$
Value	1	0.3	1.	0.7	0.001	0.01	5.08	0

.

The boundary conditions for u and $\mathsf{\Theta }$ are illustrated in Figure 9 left. For z, we set $\frac{\partial z}{\partial n}=0$ on $\mathsf{\Gamma }$.

In Figure 10, the numerical results obtained by F-PFM, TF-PFM1, and TF-PFM2 are shown in the upper, middle, and bottom parts, respectively, where we set $\delta =0.5$ for TF-PFM1 and TF-PFM2. In addition, the profile of z on line $x_2=0$ is shown in Figure 11. From Figure 10 and Figure 11, we observe that the crack propagation rate obtained by F-PFM is slower than that obtained by the others, and that the crack propagation rate obtained by TF-PFM1 is slightly faster than that obtained by TF-PFM2.

The temperature distributions obtained by TF-PFM1 and TF-PFM2 are shown in Figure 12. In the equation for $\mathsf{\Theta }$, the heat resource is given by $-{(1-z)}^2\delta \frac{d}{dt}\left(\mathrm{div}u\right)$. During crack propagation $(0.4\le t\le 0.8)$, the areas near the crack tip, the upper-right corner, and lower-right corner are continuously expanding when $\mathrm{div}u>0$ and $\frac{\partial }{\partial t}\left(\mathrm{div}u\right)>0$. Therefore, due to the negative source $-\frac{\partial }{\partial t}\left(\mathrm{div}u\right)$, lower temperatures are observed in those areas. On the other hand, at $t=1$, due to the sudden compression caused by the total fracture, positive heat is generated, and a higher temperature is observed, especially near the upper-right and lower-right corners. In this condition, it does not allow for temperature discontinuities along the crack even if we set $\kappa \left(z\right)={(1-z)}^2{\kappa }_0$.

To see how the thermoelastic coupling parameters contribute to enhanced crack propagation, we consider $\delta =0,0.1,0.2,0.5$ for TF-PFM1 and TF-PFM2, and their elastic and surface energies are plotted in Figure 13. From Figure 13, we observe that faster crack propagation occurs with a larger coupling parameter. The figure also shows that crack propagation using TF-PFM1 is faster than that using TF-PFM2.

4.4. Thermoelastic Effect on the Crack Path

In this section, we investigate the effect of the thermoelastic coupling parameter on crack path selection using our proposed models. Under a given temperature gradient, we consider crack propagation of an opening mode (Mode I) and a mixed mode (Mode I + II). In the following numerical examples, we also use the parameters in Table 4.

4.4.1. Mode I

We use an edge-cracked square domain, which is shown in Figure 14 left. We set the domain as follows:

$$\begin{array}{cc}\hfill & C_{\pm }:=\left(\begin{array}{c}-\frac{1}{2}\hfill \\ \pm \frac{5}{8}\hfill \end{array}\right)\in {\mathbb{R}}^2,\hfill \\ \hfill & H_{\pm }:=\left\{x\in {\mathbb{R}}^2;\left|x-C_{\pm }\right|\le \frac{3}{20}\right\},\hfill \\ \hfill & \mathsf{\Omega }:={(-1,1)}^2\backslash (H_{+}\cup H_{-}),\hfill \end{array}$$

and we define

$$\begin{array}{cc}\hfill & {\mathsf{\Gamma }}_{DN1}^u:=\mathsf{\Gamma }\cap \{x_1=1\},{\mathsf{\Gamma }}_{DN2}^u:=\partial H_{+}\cup \partial H_{-},{\mathsf{\Gamma }}_N^u:=\mathsf{\Gamma }\backslash ({\mathsf{\Gamma }}_{DN1}^u\cup {\mathsf{\Gamma }}_{DN2}^u),\hfill \\ \hfill & {\mathsf{\Gamma }}_{\pm D}^{\mathsf{\Theta }}:=\mathsf{\Gamma }\cap \{x_2=\pm 1\},{\mathsf{\Gamma }}_N^{\mathsf{\Theta }}:=\mathsf{\Gamma }\backslash ({\mathsf{\Gamma }}_{+D}^{\mathsf{\Theta }}\cup {\mathsf{\Gamma }}_{-D}^{\mathsf{\Theta }}).\hfill \end{array}$$

The boundary conditions for u and $\mathsf{\Theta }$ are given as follows:

$$\begin{array}{c}\left\{\begin{array}{c}{u}_1=0\hfill \\ {\sigma }_{12}^{*}=0\hfill \end{array}\right.\mathrm{on}{\mathsf{\Gamma }}_{DN1}^u,\hspace{1em}\left\{\begin{array}{c}\left({\sigma }^{*}n\right)·e_1=0\hfill \\ u_2=\pm 8t\hfill \end{array}\right.\mathrm{on}\partial H_{\pm },\hspace{1em}{\sigma }^{*}[u,\mathsf{\Theta }]n=0\mathrm{on}{\mathsf{\Gamma }}_N^u,\hfill \\ \mathsf{\Theta }={\mathsf{\Theta }}_D\mathrm{on}{\mathsf{\Gamma }}_{+D}^{\mathsf{\Theta }},\hspace{1em}\mathsf{\Theta }=0\mathrm{on}{\mathsf{\Gamma }}_{-D}^{\mathsf{\Theta }},\hspace{1em}\frac{\partial \mathsf{\Theta }}{\partial n}=0\mathrm{on}{\mathsf{\Gamma }}_N^{\mathsf{\Theta }}.\hfill \end{array}$$

The initial condition for $\mathsf{\Theta }$ is given as ${\mathsf{\Theta }}_{*}=0$.

For z, similar to the previous example (Section 4.3), we set $\frac{\partial z}{\partial n}=0$ on $\mathsf{\Gamma }$ and choose the initial value as $z_{*}\left(x\right):=e^{(-{(x_2/\eta )}^2)}/(1+e^{((x_1+0.2)/\eta )})$ with $\eta =1.5\times {10}^{-2}$. In this numerical experiment, we apply the thermoelastic coupling parameter $\delta =0.5$.

Figure 15 shows the different crack paths obtained by the three models when ${\mathsf{\Theta }}_D=10$. Straight cracks occur in the F-PFM path since the thermal effect is ignored there. On the other hand, crack curves occur in the TF-PFM1 and TF-PFM2 paths. Here, the crack path is more curved in the TF-PFM2 path than in the TF-PFM1 path. These results show good qualitative agreement with the results reported in [44].

Figure 16 shows the crack paths for different temperature gradients ${\mathsf{\Theta }}_D=0,3,5,7,10$ obtained by TF-PFM1 (left) and TF-PFM2 (right). A larger temperature gradient generates a more curved crack path, and TF-PFM2 obtains a more curved crack path than TF-PFM1. Both have significant differences in the magnitude of angle deviation but have the same crack path directions. Therefore, it is clear that thermal expansion changes the crack path.

The temperature distributions during crack growth are shown in Figure 17. There exists a temperature discontinuity along the crack path, which is caused by $\kappa \left(z\right)={(1-z)}^2{\kappa }_0$. It approximately represents a thermal insulation condition across the crack. Different from the previous condition in Section 4.3, although we involve $\frac{\partial }{\partial t}\left(\mathrm{div}u\right)$, its contribution is small.

4.4.2. Mode I + II

According to the numerical experiment in [26], we consider the following setting for mixed mode crack propagation under a thermal gradient. Let $\mathsf{\Omega }:={(-1,1)}^2\in {\mathbb{R}}^2$, as shown in Figure 14 right, and $\mathsf{\Gamma }:=\partial \mathsf{\Omega }$. We set

$$\begin{array}{cc}\hfill & {\mathsf{\Gamma }}_{\pm D}^u:=\mathsf{\Gamma }\cap \{x_2=\pm 1\},{\mathsf{\Gamma }}_N^u:=\mathsf{\Gamma }\backslash ({\mathsf{\Gamma }}_{+D}^u\cup {\mathsf{\Gamma }}_{-D}^u),\hfill \\ \hfill & {\mathsf{\Gamma }}_{\pm D}^{\mathsf{\Theta }}:=\mathsf{\Gamma }\cap \{x_2=\pm 1\},{\mathsf{\Gamma }}_N^{\mathsf{\Theta }}:=\mathsf{\Gamma }\backslash ({\mathsf{\Gamma }}_{+D}^{\mathsf{\Theta }}\cup {\mathsf{\Gamma }}_{-D}^{\mathsf{\Theta }}).\hfill \end{array}$$

The boundary conditions for u are given as follows:

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{u}_1=\pm 3sin(\pi /3)t,\hfill \\ u_2=\pm 3cos(\pi /3)t\hfill \end{array}\right.\mathrm{on}{\mathsf{\Gamma }}_{\pm D}^u,\hspace{1em}{\sigma }^{*}[u,\mathsf{\Theta }]n=0\mathrm{on}{\mathsf{\Gamma }}_N^u.\hfill \end{array}$$

The boundary conditions for $\mathsf{\Theta }$ and z are the same as those in Section 4.4.1. The initial crack profile is given as $z_{*}\left(x\right):=e^{(-{(x_2/\eta )}^2)}/(1+e^{((x_1-0.5)/\eta )})-e^{(-{(x_2/\eta )}^2)}/(1+e^{((x_1+0.5)/\eta )})$ with $\eta =1.5\times {10}^{-2}$. We fix the thermoelastic coupling parameter $\delta =0.15$ and change the temperature gradient to ${\mathsf{\Theta }}_D=0,2,3,5,6$.

Figure 18 shows the crack paths obtained by TF-PFM1 and TF-PFM2. The cracks are kinked, and the kink angle becomes larger when the thermal gradient ${\mathsf{\Theta }}_D$ increases. The two models provide similar results, but the kink angle in the TF-PFM2 crack is larger than that in the TF-PFM1 crack, as shown in Figure 19. Therefore, we conclude that thermal expansion changes the crack path.

Here, we do not show the temperature distribution during thermal expansion. We observe that the temperature distribution is quite similar to that of Mode I in Section 4.4.1, and a temperature discontinuity exists along the crack path during temperature injection. As mentioned, since it is relatively difficult to find the available experimental result for the thermal fracturing under mode I + II, we do not compare our result with the experimental result.

At the end of this section, we give a remark on the extendability of our TF-PFM to anisotropic material. When the material has a strong anisotropy, we have to take into account the anisotropies on the elasticity tensor C, a coefficient of linear thermal expansion $a_L$, and the critical energy release rate ${\gamma }_{*}$, especially among many material properties. For C and $a_L$, we can easily include an anisotropic effect by using an anisotropic tensor in (3b) and replace the matrix $a_{L}I$ in (4a) by anisotropic one. On the other hand, for ${\gamma }_{*}$, it has not well succeeded to include the anisotropy, which means that the dependency of the crack surface direction, even in the standard PFM, as far as the authors’ knowledge.

5. Conclusions and Future Works

We proposed two thermal fracturing phase field models, TF-PFM1 and TF-PFM2, by coupling the Biot thermoelasticity model [8] and the fracturing phase field model (F-PFM) by Takaishi–Kimura [22,26].

For the Biot model, we studied a variational principle (Proposition 2) and energy equality (Theorem 1), which were related to different energies $E_{el}^{*}(u,\mathsf{\Theta })$ and $E_{el}\left(u\right)+E_{th}\left(\mathsf{\Theta }\right)$, respectively (see Table 1 and Table 2).

On the other hand, F-PFM has a gradient flow structure with respect to the total energy ${\mathcal{E}}_{el}(u,z)+E_s\left(z\right)$ and admits energy equality (Theorem 2).

As the first model, TF-PFM1 was derived based on the variational principle of the Biot model and the gradient flow structure of F-PFM, while TF-PFM2 is based on the energy equalities of the Biot model and F-PFM. The difference between them is the driving force term for the crack: $W^{*}(u,\mathsf{\Theta })$ in TF-PFM1 (31b) and $W\left(u\right)$ in TF-PFM2 (35b).

Consequently, we established partial energy equality for TF-PFM1 (Theorem 3) and energy equality for TF-PFM2 (Theorem 4). From the viewpoint of energy consistency, both models are satisfactory, but TF-PFM2 is more energetically consistent than TF-PFM1.

Based on the obtained numerical experiments, the following conclusions can be drawn.

• The thermoelastic coupling parameter $\delta$ in TF-PFM1 and TF-PFM2 enhances crack propagation (Figure 10).
• TF-PFM1 accelerates the crack speed more than TF-PFM2 (Figure 11). On the other hand, the effect of the temperature gradient on the crack path in TF-PFM2 is larger than that in TF-PFM1 (Figure 16, Figure 17, Figure 18 and Figure 19).

The analytical and numerical comparisons between the two models are briefly summarized in

Table 5. Numerical comparison of TF-PFM1 and TF-PFM2.

Models	Driving Force	Energy Consistency	Straight Crack Speed	Crack Path
TF-PFM1	$W^{*}(u,\mathsf{\Theta })={\sigma }^{*}[u,\mathsf{\Theta }]:e^{*}[u,\mathsf{\Theta }]$	Partially satisfied	Faster	Less curved
TF-PFM2	$W\left(u\right)=\sigma \left[u\right]:e\left[u\right]$	Fully satisfied	Slower	More curved
Remarks	$W^{*}(u,\mathsf{\Theta })>W\left(u\right)$ (Figure 4)	Theorems 3 & 4	Figure 13	Figure 16 and Figure 19

.

In this study, we did not consider the unilateral contact condition along the crack for the sake of simplicity. To further improve TF-PFM, the ideal unilateral condition for fracturing PFM [21,26] should be introduced in our PFM. In addition, there are many other effects that should be included in the model. For example, although we assumed that the critical energy release rate ${\gamma }_{*}\left(x\right)$ is a priory given, it may depend on the temperature in the real material. A possible extension of our model is to suppose that ${\gamma }_{*}$ depends on $\mathsf{\Theta }$ linearly as

$$\begin{array}{c}\hfill {\gamma }_{*}(x,\mathsf{\Theta })=\overline{\gamma }(1-{\alpha }_0(\mathsf{\Theta }-{\mathsf{\Theta }}_0))\end{array}$$

for some $\overline{\gamma }>0$ and $\alpha >0$ [45].

Such relatively easy extendability is one of the advantages of PFM. However, we should remark that the energy equalities which we derived in this paper may not be valid for all extended models.