# Phase field models for thermal fracturing and their variational structures
11footnotemark: 1

S. Alfat<EMAIL_ADDRESS>M. Kimura<EMAIL_ADDRESS>Alifian. M. M<EMAIL_ADDRESS>Physics Education Department, Halu Oleo
University, Kendari, Southeast Sulawesi, Indonesia  Division of Mathematical
and Physical Sciences, Graduate School on Natural Science and Technology,
Kanazawa University, Kakuma, Kanazawa, 920-1192, Japan  Faculty of Mathematics
and Physics, Kanazawa University, Kakuma, Kanazawa, 920-1192, Japan

###### 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 (1956) 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.

###### keywords:

Thermoelasticity , Crack Propagation , Crack Path , Phase Field Model ,
Variational Structure , Energy Equality , Adaptive Finite Element Method

††journal: Journal of LaTeX Templates

[table]capposition=top

##  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
[29, 31, 35], fluid pressure (e.g., in hydraulic fracturing) [26], the
diffusion of hydrogen (or hydrogen embrittlement) [13, 28], chemical reactions
[16], and humidity [31], 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 [5]. Biot’s model is now widely
known as the traditional coupled thermoelasticity model, and it has been
extended to dynamical theory [27] and to various other situations [14, 20, 21,
25, 36, 37]. 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. [6] and Karma et al. [23]. Then, based on the framework of variational
fracture theory [8, 15], the techniques and applications of PFM have been
extensively developed, for example [2, 3, 7, 30, 34]. We refer to [24] 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 [6, 24]. 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 [33].

In this study, we deal with the modeling of thermal fracturing in an isotropic
and homogeneous body by coupling the Biot thermoelasticity model and PFM.
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 [1, 9, 11, 24, 30, 32], but they neglect the strain’s influence on
the heat transfer. To the best of our knowledge, a peridynamics model that
employs the coupled thermoelastic equation was proposed by Gao and Oterkus
[18].

---
Figure 1: A conceptual diagram of energy balance for Biot’s model, F-PFM, and
TF-PFM

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 3.1) and follows the works [24, 34]. In Sections 3.2 and
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.3), we show a partial energy equality for a fixed temperature
(Theorem 3.2). 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 2.4). We prove an energy equality for TF-
PFM2 in Theorem 3.3. 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 Tables 1 and 2.

In Section 4, we show some numerical comparisons between two TF-PFMs using
nondimensionalized 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$. The last section shows
some conclusions and comments on further topics.

Table 1: List of energies Type of Energy |  Definition |  Equation
---|---|---
Elastic |  $E_{el}(u):=\dfrac{1}{2}\displaystyle\int_{\Omega}\sigma[u]:e[u]~{}dx$ |  (2.16)
Thermoelastic |  $E_{el}^{*}(u,\Theta):=\dfrac{1}{2}\displaystyle\int_{\Omega}\sigma^{*}[u,\Theta]:e^{*}[u,\Theta]~{}dx$ |  (2.24)
Thermal |  $E_{th}(\Theta):=\dfrac{\chi}{2\Theta_{0}}\displaystyle\int_{\Omega}\big{|}\Theta(x)-\Theta_{0}\big{|}^{2}~{}dx$ |  (2.28)
Modified elastic |  $\mathcal{E}_{el}(u,z):=\dfrac{1}{2}\displaystyle\int_{\Omega}(1-z)^{2}\sigma[u]:e[u]~{}dx$ |  (3.3)
Modified thermoelastic |  $\mathcal{E}_{el}^{*}(u,\Theta,z):=\dfrac{1}{2}\displaystyle\int_{\Omega}(1-z)^{2}\sigma^{*}[u,\Theta]:e^{*}[u,\Theta]~{}dx$ |  (3.8)
Surface |  $E_{s}(z):=\dfrac{1}{2}\displaystyle\int_{\Omega}\gamma_{*}\Big{(}\epsilon|\nabla z|^{2}+\frac{|z|^{2}}{\epsilon}\Big{)}~{}dx$ |  (3.4)

[b]  Model Strong Form Energy Energy Equality Linear elasticity (2.14)
$E_{el}(u)$ - Biot’s model (2.1a) - (2.1b) $E_{el}(u)+E_{th}(\Theta)$ (2.29)
F-PFM (3.1a) - (3.1b) $\mathcal{E}_{el}(u,z)+E_{s}(z)$ (3.6) TF-PFM1 (3.10a) -
(3.10c) $\mathcal{E}_{el}^{*}(u,\Theta,z)+E_{s}(z)$ (3.12)a TF-PFM2 (3.14a) -
(3.14c) $\mathcal{E}_{el}(u,z)+E_{s}(z)+E_{th}(\Theta)$ (3.15)

Table 2: Different forms of energy equalities

  * a

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

[b]  Symbol Physical meaning [unit] Symbol Physical meaning [unit] $u$
Displacement [$\mathrm{m}$] $\sigma^{*}[u,\Theta]$ Stress tensor with thermal
effect [$\mathrm{Pa}$] $\Theta$ Temperature [$\mathrm{K}$] $e^{*}[u,\Theta]$
Strain tensor with thermal effect [-] $\Theta_{0}$ Reference temperature
[$\mathrm{K}$] $\beta$ Stress thermal modulus
[$\mathrm{Pa}\cdot\mathrm{K}^{-1}$] $z$ Damage variable [-] $\kappa_{0}$
Thermal conductivity [$\mathrm{W}\cdot\mathrm{m}^{-1}\cdot\mathrm{K}^{-1}$]
$\sigma[u]$ Stress tensor [$\mathrm{Pa}$] $\chi$ Volumetric heat capacity
[$\mathrm{J}\cdot\mathrm{K}^{-1}\cdot\mathrm{m}^{-3}$] $e[u]$ Strain tensor
[-] $a_{L}$ Coefficient of linear thermal expansion [${\mathrm{K}}^{-1}$]
$E_{\text{Y}}$ Young’s modulus [$\mathrm{Pa}$] ${\delta}$ Thermoelastic
coupling parameter [-] $\nu_{\text{P}}$ Poisson ratio [-] $\gamma_{*}$
Critical energy release ratea [$\mathrm{Pa}\cdot\mathrm{m}$] $\lambda$, $\mu$
Lamé’s constantsb [$\mathrm{Pa}$] $\epsilon$ 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}\cdot\mathrm{s}$]

Table 3: List of physical properties

  * a

$\gamma_{*}$ is usually denoted by $\text{G}_{\text{c}}$ [4, 24].

  * b

$\lambda$ and $\mu$ are written as
$\lambda=\cfrac{E_{\text{Y}}\nu_{\text{P}}}{(1+\nu_{\text{P}})(1-2\nu_{\text{P}})}$
and $\mu=\cfrac{E_{\text{Y}}}{2(1-\nu_{\text{P}})}$.

To easily understand the relevant notation and symbols in this paper, we
introduce them in this section. Let $\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}.\cdots,x_{d})^{T}\in\mathbb{R}^{d}$, where $~{}^{T}$ denotes the
transposition of a vector or matrix. Let $\nabla$, div, and $\Delta$ be the
gradient, divergence, and Laplacian operators with respect to $x$,
respectively. For simplicity, we write $\dot{u}$, $\dot{\Theta}$, and
$\dot{z}$ as the partial derivatives of $u$, $\Theta$ and $z$ with respect to
$t$, respectively. For simplicity, we often denote $u(t):=u(\cdot,t)$, etc.
The space of the real-valued (symmetric) $d\times d$ matrix is denoted by
$\mathbb{R}^{d\times d}$ ($\mathbb{R}^{d\times d}_{sym}$). 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}(\Omega)$, we refer to the
Lebesgue space on $\Omega$, while $H^{1}(\Omega,\mathbb{R}^{d})$ and
$H^{\frac{1}{2}}(\Gamma_{D}^{u},\mathbb{R}^{d})$ represent the Sobolev space
on $\Omega$ and its trace space on the boundary $\Gamma_{D}^{u}$,
respectively. For more details on Sobolev spaces, we refer to the review in
[19]. In addition, we summarize the physical properties used in this paper in
Table 3.

##  2 Thermoelasticity Model

###  2.1 Formulation of the problem

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

$\displaystyle-\mbox{div}\sigma{[u]}=\beta\nabla\Theta$
$\displaystyle\mbox{in}~{}\Omega\times[0,T],$ (2.1a)
$\displaystyle\chi\dfrac{\partial}{\partial
t}\Theta=\kappa_{0}\Delta\Theta-\Theta_{0}\beta\dfrac{\partial}{\partial
t}(\mbox{div}u)$ $\displaystyle\mbox{in}~{}\Omega\times(0,T],$ (2.1b)

where $\Omega$ is a bounded domain in $\mathbb{R}^{d}$ ($d=2$ or $3$). We
suppose that $\Omega$ is an isotropic elastic body and consider the
thermoelastic coupling between the mechanical deformation and the thermal
expansion in $\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 (2.1a) and (2.1b) are the displacement
$u(x,t)=(u_{1}(x,t)$ $,\cdots,u_{d}(x,t))^{T}\in\mathbb{R}^{d}$ and the
temperature $\Theta(x,t)\in\mathbb{R}$. In addition, the constant
$\Theta_{0}>0$ is a fixed reference temperature. Similarly, strain $e[u]$ and
stress tensors $\sigma[u]$ are defined as

$\displaystyle e[u]:=\frac{1}{2}\left(\nabla u^{T}+(\nabla
u^{T})^{T}\right)\in\mathbb{R}_{sym}^{d\times d},$ (2.2a)
$\displaystyle\sigma[u]:=Ce[u]=\lambda(\mbox{div}u){I}+2\mu
e[u]\in\mathbb{R}_{sym}^{d\times d},$ (2.2b)

where
$C:=(c_{ijkl}),~{}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 (2.2b), (2.1a) is also written in the form

$\displaystyle-\mu\Delta
u-(\lambda+\mu)\nabla(\mbox{div}u)=\beta\nabla\Theta.$

The term $\beta\nabla\Theta$ in (2.1a) and the term
$\Theta_{0}\beta\frac{\partial}{\partial t}(\mbox{div}u)$ in (2.1b) 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$, (2.1a) and (2.1b)
are decoupled.

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

$\displaystyle
e^{*}[u,\Theta]:=e[u]-a_{L}(\Theta(x,t)-\Theta_{0})I\in\mathbb{R}^{d\times
d}_{sym},$
$\displaystyle\sigma^{*}[u,\Theta]:=Ce^{*}[u,\Theta]=\sigma[u]-\beta(\Theta(x,t)-\Theta_{0})I\in\mathbb{R}^{d\times
d}_{sym}.$

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

$\displaystyle-\mbox{div}\sigma^{*}[u,\Theta]=0.$

This means that the force $\sigma^{*}[u,\Theta]$ is in equilibrium in
$\Omega$. In the preceding equation, (2.1) represents the force balance and
the thermal diffusion in $\Omega$, respectively.

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

$\displaystyle u=u_{D}(x,t)$
$\displaystyle\mbox{on}~{}\Gamma_{D}^{u}\times[0,T],$ (2.4a)
$\displaystyle\sigma^{*}[u,\Theta]n=0$
$\displaystyle\mbox{on}~{}\Gamma_{N}^{u}\times[0,T],$ (2.4b)
$\displaystyle\Theta=\Theta_{D}(x,t)$
$\displaystyle\mbox{on}~{}{\Gamma}_{D}^{\Theta}\times[0,T],$ (2.4c)
$\displaystyle\frac{\partial\Theta}{\partial n}=0$
$\displaystyle\mbox{on}~{}{\Gamma}_{N}^{\Theta}\times[0,T],$ (2.4d)
$\displaystyle\Theta(x,0)=\Theta_{*}(x)$
$\displaystyle\mbox{in}~{}\Omega,~{}\qquad\qquad$ (2.4e)

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

###### Remark 2.1

Instead of boundary conditions (2.4a) and (2.4b), we can also consider the
following mixed-type condition. When $d=2$, on a part of the boundary (which
we denote by $\Gamma_{DN}^{u}$), $u=(u_{1},u_{2})^{T}$ and

$\displaystyle\left\\{\begin{array}[]{ll}u_{1}=u_{D1}&\qquad\mbox{on}~{}\Gamma_{DN}^{u},\\\
(\sigma^{*}[u,\Theta]n)\cdot
e_{2}=0&\qquad\mbox{on}~{}\Gamma_{DN}^{u},\end{array}\right.$ (2.7)

or

$\displaystyle\left\\{\begin{array}[]{ll}u_{2}=u_{D2}&\qquad\mbox{on}~{}\Gamma_{DN}^{u},\\\
(\sigma^{*}[u,\Theta]n)\cdot
e_{1}=0&\qquad\mbox{on}~{}\Gamma_{DN}^{u},\end{array}\right.$ (2.10)

where $u_{Di}:=\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 (2.1). In
linear elasticity theory, a weak form of the boundary value problem for
$u_{D}\in H^{\frac{1}{2}}(\Gamma_{D}^{u};\mathbb{R}^{d})$ is

$\displaystyle\left\\{\begin{array}[]{ll}-\mbox{div}\sigma[u]=0&\qquad\mbox{in}~{}\Omega,\\\
u=u_{D}&\qquad\mbox{on}~{}\Gamma_{D}^{u},\\\
\sigma[u]n=0&\qquad\mbox{on}~{}\Gamma_{N}^{u},\end{array}\right.$ (2.14)

which is given by

$\displaystyle u\in
V^{u}(u_{D}),~{}\int_{\Omega}\sigma[u]:e[v]~{}dx=0\quad\text{for all}~{}v\in
V^{u}(0),$

where

$\displaystyle V^{u}(u_{D}):=\left\\{u\in
H^{1}(\Omega;\mathbb{R}^{d});~{}u\big{|}_{\Gamma_{D}^{u}}=u_{D}\right\\}.$
(2.15)

A weak solution uniquely exists and is given by

$u=\operatorname*{argmin}_{v\in V^{u}(u_{D})}{E_{el}(v)},$

where

$\displaystyle
E_{el}(v):=\frac{1}{2}\int_{\Omega}\sigma[v]:e[v]~{}dx\quad(v\in
H^{1}\left(\Omega;\mathbb{R}^{d})\right)$ (2.16)

is an elastic energy. This is known as a variational principle [10, 12]. For a
fixed $\Theta(x)$, a weak form for $u$ of (2.1a) and its variational principle
are derived as follows.

###### Proposition 2.2

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

$\displaystyle\left\\{\begin{array}[]{ll}-{\mbox{div}}\sigma^{*}[u,\Theta]=0&\qquad\mbox{in}~{}\Omega,\\\
u=u_{D}&\qquad\mbox{on}~{}\Gamma_{D}^{u},\\\
\sigma^{*}[u,\Theta]n=0&\qquad\mbox{on}~{}\Gamma_{N}^{u},\end{array}\right.$
(2.20)

is equivalent to the following weak form:

$\displaystyle\left\\{\begin{array}[]{ll}\displaystyle\int_{\Omega}\sigma^{*}[u,\Theta]:e[v]~{}dx=0&\quad\mbox{for
all}~{}v\in V^{u}(0),\\\\[8.0pt] u\in V^{u}(u_{D}).&\end{array}\right.$ (2.23)

###### Proof

For $v\in V^{u}(0)$, we have

$\displaystyle\int_{\Omega}\left(-\mbox{div}\sigma^{*}[u,\Theta]\right)\cdot
v~{}dx$ $\displaystyle=$
$\displaystyle\int_{\Omega}\sigma^{*}[u,\Theta]:e[v]~{}dx-\int_{\Gamma_{N}^{u}}(\sigma^{*}[u,\Theta]n)\cdot
v~{}ds.$

The equivalency immediately follows from this equation.

######  Proposition 2.3 (Variational principle)

For a given $\Theta\in L^{2}(\Omega)$, $u_{D}\in
H^{\frac{1}{2}}(\Gamma_{D}^{u};\mathbb{R}^{d})$, there exists a unique weak
solution $u\in H^{1}(\Omega;\mathbb{R}^{d})$ that satisfies (2.23).
Furthermore, the solution $u$ is a unique minimizer of the variational
problem:

$u=\operatorname*{argmin}_{v\in V^{u}(u_{D})}{E_{el}^{*}(v,\Theta)},$

where

$\displaystyle
E_{el}^{*}(v,\Theta)=\frac{1}{2}\int_{\Omega}\sigma^{*}[v,\Theta]:e^{*}[v,\Theta]~{}dx.$
(2.24)

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

###### Proof

The unique existence of a weak solution for $u$ is shown by the Lax-Milgram
theorem [10] since (2.23) is written as

$\displaystyle\left\\{\begin{array}[]{ll}\displaystyle\int_{\Omega}\sigma[u]:e[v]~{}dx=\displaystyle\int_{\Omega}\beta(\Theta-\Theta_{0})\mbox{div}v~{}dx,\\\\[3.0pt]
u\in V^{u}(u_{D})\qquad(\mbox{for all}~{}v\in V^{u}(0)).\end{array}\right.$
(2.27)

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

${}^{\exists}a_{0}>0~{}\text{such that}~{}\int_{\Omega}\sigma[v]:e[v]~{}dx\geq
a_{0}\|v\|^{2}_{H^{1}(\Omega;\mathbb{R}^{2})},\quad\mbox{for all}~{}v\in
V^{u}(0).$

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

$\displaystyle\sigma^{*}[u+v,\Theta]=\sigma^{*}[u,\Theta]+\sigma[v],$
$\displaystyle e^{*}[u+v,\Theta]=e^{*}[u,\Theta]+e[v],$
$\displaystyle\sigma^{*}[u,\Theta]:e[v]=e^{*}[u,\Theta]:\sigma[v],$

we have

$\displaystyle E_{el}^{*}(u+v,\Theta)-E_{el}^{*}(u,\Theta)$
$\displaystyle\qquad\quad=\frac{1}{2}\int_{\Omega}\sigma^{*}[u+v,\Theta]:e^{*}[u+v,\Theta]~{}dx-\frac{1}{2}\int_{\Omega}\sigma^{*}[u,\Theta]:e^{*}[u,\Theta]~{}dx\qquad$
$\displaystyle\qquad\quad=\int_{\Omega}\sigma^{*}[u,\Theta]:e[v]~{}dx+\frac{1}{2}\int_{\Omega}\sigma[v]:e[v]~{}dx$
$\displaystyle\qquad\quad=\frac{1}{2}\int_{\Omega}\sigma[v]:e[v]~{}dx\geq 0.$

This shows that $u$ is a minimizer of $E^{*}_{el}(u,\Theta)$ among
$V^{u}(u_{D})$.

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}(0)$, we have

$\displaystyle
0=\frac{d}{ds}E^{*}_{el}(u+sv,\Theta)\big{|}_{s=0}=\int_{\Omega}\sigma^{*}[u,\Theta]:e[v]~{}dx.$

Hence, $u$ is a weak solution. Summarizing the above, there exists a unique
weak solution to (2.16), and $u$ is a weak solution if and only if it is a
minimizer of $E^{*}_{el}$ among $V^{u}(u_{D})$.

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

$\displaystyle
E_{th}(\Theta):=\frac{\chi}{2\Theta_{0}}\int_{\Omega}\left|\Theta(x)-\Theta_{0}\right|^{2}~{}dx.$
(2.28)

######  Theorem 2.4 (Energy equality for Biot’s model)

Let $(u(x,t),\Theta(x,t))$ be a sufficiently smooth solution to (2.1) and
(2.4). In addition, we suppose that $u_{D}$ does not depend on $t$ and
$\Theta_{D}=\Theta_{0}$. Then

$\displaystyle\frac{d}{dt}\Big{(}E_{el}(u(t))+E_{th}(\Theta(t))\Big{)}$
$\displaystyle=-\frac{\kappa_{0}}{\Theta_{0}}\int_{\Omega}\left|\nabla\Theta(t)\right|^{2}~{}dx\leq
0.$ (2.29)

###### Proof

Since

$\displaystyle\frac{d}{dt}\left(\frac{1}{2}\sigma[u]:e[u]\right)$
$\displaystyle=\sigma[u]:e[\dot{u}]$
$\displaystyle=(\sigma^{*}[u,\Theta]+\beta(\Theta-\Theta_{0})I):e[\dot{u}]$
$\displaystyle=\sigma^{*}[u,\Theta]-\beta(\Theta-\Theta_{0})\mbox{div}\dot{u}$
(2.30)

we obtain

$\displaystyle\frac{d}{dt}E_{el}(u(t))$
$\displaystyle=\frac{1}{2}\int_{\Omega}\frac{d}{dt}\left(\sigma[u]:e[u]\right)~{}dx$
$\displaystyle=\int_{\Omega}\sigma^{*}[u,\Theta]:e[\dot{u}]~{}dx+\int_{\Omega}\beta(\Theta-\Theta_{0})(\mbox{div}\dot{u})~{}dx$
$\displaystyle=\int_{\Omega}\beta(\Theta-\Theta_{0})(\mbox{div}\dot{u})~{}dx.$
(2.31)

Substituting (2.1b) into (2.31) and using the boundary conditions (2.4c) and
(2.4d) for $\Theta$, we obtain

$\displaystyle\frac{d}{dt}E_{el}(u(t))$
$\displaystyle=\int_{\Omega}\Big{(}\frac{1}{\Theta_{0}}(\Theta-\Theta_{0})\Big{\\{}\kappa_{0}\Delta\Theta-\chi\frac{\partial\Theta}{\partial
t}\Big{\\}}\Big{)}dx$
$\displaystyle=\frac{\kappa_{0}}{\Theta_{0}}\int_{\Gamma}(\Theta-\Theta_{0})\frac{\partial\Theta}{\partial
n}ds-\frac{\kappa_{0}}{\Theta_{0}}\int_{\Omega}\big{|}\nabla\Theta\big{|}^{2}dx-\frac{d}{dt}\Big{(}\frac{\chi}{2\Theta_{0}}\int_{\Omega}\big{|}\Theta-\Theta_{0}\big{|}^{2}dx\Big{)}$
$\displaystyle=-\frac{\kappa_{0}}{\Theta_{0}}\int_{\Omega}\big{|}\nabla\Theta\big{|}^{2}dx-\frac{d}{dt}E_{th}(\Theta(t)).$

This gives the energy equality for (3.2).

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

$\displaystyle{W}(u):=\sigma[u]:e[u],$ (2.32)
$\displaystyle{W}^{*}(u,\Theta):=\sigma^{*}[u,\Theta]:e^{*}[u,\Theta],$ (2.33)

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

###  2.3 Numerical Experiment

####  2.3.1 Nondimensional setting

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

$\displaystyle\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{\Theta}=\frac{\Theta-\Theta_{0}}{c_{\Theta}},~{}\tilde{a}_{L}=\frac{c_{x}c_{\Theta}}{c_{u}}a_{L},~{}\tilde{\beta}=1,$
(2.34)

where $c_{x}$, $c_{t}$, $c_{u}$, $c_{e}$, and $c_{\Theta}>0$ are the scaling
parameters. Let $c_{x}$ [$\mathrm{m}$], $c_{e}$ [$\mathrm{Pa}$], and
$c_{\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

$\displaystyle
c_{t}:=\frac{c_{x}^{2}{\chi}}{\kappa_{0}}~{}[$\mathrm{s}$],\qquad
c_{u}:=\frac{c_{\Theta}c_{x}\beta}{c_{e}}~{}[$\mathrm{m}$],$ (2.35)

where $\chi$ [$\mathrm{Pa}\cdot\mathrm{K}^{-1}$], $\kappa_{0}$
[$\mathrm{Pa}\cdot\mathrm{m}^{2}\cdot\mathrm{s}^{-1}\cdot\mathrm{K}^{-1}$] and
$\beta=a_{L}(d\lambda+2\mu)$ [$\mathrm{Pa}\cdot\mathrm{K}^{-1}$]. Then, (2.1)
is written in the following nondimensional form:

$\displaystyle-\widetilde{\mbox{div}}\tilde{\sigma}{[\tilde{u}]}=\tilde{\nabla}\tilde{\Theta}$
$\displaystyle\mbox{in}~{}\tilde{\Omega}\times[0,\tilde{T}],$ (2.36a)
$\displaystyle\dfrac{\partial}{\partial\tilde{t}}\tilde{\Theta}=\tilde{\Delta}\tilde{\Theta}-\delta\dfrac{\partial}{\partial\tilde{t}}(\widetilde{\mbox{div}}\tilde{u})$
$\displaystyle\mbox{in}~{}\tilde{\Omega}\times(0,\tilde{T}].$ (2.36b)

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

$\displaystyle{\delta}$ $\displaystyle:=$
$\displaystyle\frac{\Theta_{0}\beta^{2}}{c_{e}\chi}~{}[-],$

and $\delta>0$. If we choose $\delta=0$, (2.36b) is decoupled from (2.36a),
and the temperature field $\tilde{\Theta}$ in (2.36a) 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:

$\displaystyle\tilde{e}[\tilde{u}]:=\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[u],$
(2.37a)
$\displaystyle\tilde{\sigma}[\tilde{u}]:=\tilde{C}\tilde{e}[\tilde{u}]=\frac{c_{x}}{c_{u}c_{e}}\sigma[u],$
(2.37b)
$\displaystyle\widetilde{{W}}(\tilde{u}):=\tilde{\sigma}[\tilde{u}]:\tilde{e}[\tilde{u}]=\frac{c_{e}}{(\beta
c_{\Theta})^{2}}{W}[u],$ (2.37c)
$\displaystyle\tilde{\sigma}^{*}[\tilde{u},\tilde{\Theta}]:=\tilde{\sigma}[\tilde{u}]-\tilde{\Theta}I=\frac{1}{\beta
c_{\Theta}}\sigma^{*}[u,\Theta],$ (2.37d)
$\displaystyle\tilde{e}^{*}[\tilde{u},\tilde{\Theta}]:=\tilde{e}[\tilde{u}]-\tilde{a}_{L}\tilde{\Theta}I=\frac{c_{x}}{c_{u}}\sigma^{*}[u,\Theta],$
(2.37e)
$\displaystyle\widetilde{{W}}^{*}(\tilde{u},\tilde{\Theta}):=\tilde{\sigma}^{*}[\tilde{u},\tilde{\Theta}]:\tilde{e}^{*}[\tilde{u},\tilde{\Theta}]=\frac{c_{e}}{(\beta
c_{\Theta})^{2}}{W}^{*}[u,\Theta].$ (2.37f)

In the following section, we apply these nondimensional forms and omit $\sim$
for simplicity.

####  2.3.2 Numerical setup and time discretization

In the following examples, we set Young’s modulus $E_{\text{Y}}=1$, Poisson’s
ratio $\nu_{\text{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 nondimensional form of (2.36). We consider
two numerical examples for (2.36), an L-shaped cantilever domain and a square
domain with a crack (more precisely, a very sharp notch), as illustrated in
Figure 2.

|
---|---
Figure 2: An L-shaped cantilever (left) and a cracked domain (right) with the
subdomain $\mathcal{A}$ as an observation area.

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

$\displaystyle\left\\{\begin{array}[]{ll}-\displaystyle{{\text{div}}{\sigma}^{*}[{u}^{k},\Theta^{k-1}]=0}&\displaystyle\mbox{in}~{}\Omega,\\\\[5.0pt]
\displaystyle\frac{{\Theta}^{k}-{\Theta}^{k-1}}{\Delta{t}}-{\Delta}{\Theta}^{k}+{\delta}{\text{div}}\left(\frac{{u}^{k}-{u}^{k-1}}{\Delta{t}}\right)=0&\displaystyle\mbox{in}~{}\Omega,\end{array}\right.$
(2.40)

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

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

$\displaystyle\mathcal{W}(\mathcal{A})$
$\displaystyle:=\frac{1}{|\mathcal{A}|}\int_{\mathcal{A}}{W}(u)~{}dx,$
$\displaystyle\mathcal{W}^{*}(\mathcal{A})$
$\displaystyle:=\frac{1}{|\mathcal{A}|}\int_{\mathcal{A}}{W}^{*}(u,\Theta)~{}dx,$

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

$\displaystyle\Delta\mathcal{W}(\mathcal{A}):=\mathcal{W}(\mathcal{A})\big{|}_{\delta}-\mathcal{W}(\mathcal{A})\big{|}_{\delta=0}~{},$
$\displaystyle\Delta\mathcal{W}^{*}(\mathcal{A}):=\mathcal{W}^{*}(\mathcal{A})\big{|}_{\delta}-\mathcal{W}^{*}(\mathcal{A})\big{|}_{\delta=0}~{}.$

In the following examples, we use the software FreeFEM [17] with P2 elements
and unstructured meshes. For the time interval and time step, we use $0\leq
t\leq 0.1$ and $\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 $\Gamma_{D}^{u}$
and $\Gamma_{DN}^{u}$, respectively, and define
$\Gamma_{N}^{u}:=\Gamma\setminus(\Gamma_{D}^{u}\cup\Gamma_{DN}^{u})$. The
boundary conditions for $u$ are

$\displaystyle
u=0,~{}\mbox{on}~{}\Gamma_{D}^{u},\qquad\left\\{\begin{array}[]{ll}\sigma_{11}^{*}[u,\Theta]n=0,\\\
u_{2}=-0.1t\end{array}\right.~{}\mbox{on}~{}\Gamma_{DN}^{u},\qquad\sigma^{*}[u,\Theta]n=0~{}\mbox{on}~{}\Gamma_{N}^{u}.$
(2.44)

For $\Theta$, we suppose $\frac{\partial\Theta}{\partial n}=0$ on $\Gamma$ and
the initial temperature $\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 (2.40). The total number of triangular
meshes $=18215$ and the number of nodes (the vertices of the triangles)
$=9301$.

|  |
---|---|---
|  |
${{t}=0}$ | ${{t}=0.05}$ | ${{t}=0.1}$
Figure 3: Snapshots of $\mbox{div}u$ (upper) and the temperature (lower) of
the L-shape cantilever for $t=0,~{}0.05,~{}0.1$ using $\delta=0.1$. Near the
re-entrant corner, the domain is expanded ($\mbox{div}u>0$), and the
temperature decreases. On the other hand, near the bottom boundary, the domain
is compressed ($\mbox{div}u<0$), and the temperature increases.

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}(u)$
and $\mathcal{W}^{*}(u)$, as shown in Figure 4 (a)-(b). Here, a larger
$\delta$ value implies larger $\mathcal{W}(\mathcal{A})$ and
$\mathcal{W}^{*}(\mathcal{A})$ values (Figure 4 (d)-(e)). In addition, we also
observe that $\mathcal{W}^{*}(\mathcal{A})$ is larger than
$\mathcal{W}(\mathcal{A})$ for each $\delta>0$ (Figure 4 (c)).

figure[3] []

[]

[]

[2] []

[]

Figure 4: Profiles of (a) $\mathcal{W}(\mathcal{A})$, (b)
$\mathcal{W}^{*}(\mathcal{A})$, (c)
$\mathcal{W}^{*}(\mathcal{A})-\mathcal{W}(\mathcal{A})$, (d)
$\Delta\mathcal{W}(\mathcal{A})$ and (e) $\Delta\mathcal{W}^{*}(\mathcal{A})$
in an L-shaped cantilever during the loading process.

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 [30], 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

$\displaystyle\left\\{\begin{array}[]{ll}u_{1}=0,\\\ u_{2}=\pm
t\end{array}\right.~{}\mbox{on}~{}\Gamma_{\pm
D}^{u},\qquad\sigma^{*}[u,\Theta]n=0~{}\mbox{on}~{}\Gamma_{N}^{u},$ (2.47)

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

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

|
---|---
Figure 5: Snapshot of $\mbox{div}u$ on the subdomain $\mathcal{A}$ (left) and temperature $\Theta$ in $\Omega$ (right) using $\delta=0.1$ at $t=0.1$. |
---|---
Figure 6: Profile of $\mbox{div}u$ (left) and temperature $\Theta$ (right)
using $\delta=0.1$ along the $x_{1}$ axis, i.e., $x_{2}=0$, $0.5\leq x_{1}\leq
1$, during the loading process.

From Figure 5 (left), we conclude that the area that expands the most (i.e.,
$\mbox{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
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 $\mbox{div}u$ is plotted along the $x_{1}$ axis and $\mbox{div}u$
is increasing over time; i.e., the heat source term $\mbox{div}\dot{u}$ in
(2.1b) is positive.

Similar to Section 2.3.3, for each $\delta>0$, we obtain variations of
$\mathcal{W}(\mathcal{A})$ and $\mathcal{W}^{*}(\mathcal{A})$ 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}^{*}(\mathcal{A})$ is larger than $\mathcal{W}(\mathcal{A})$. This
suggests that the thermoelastic energy density ${W}^{*}(u,\Theta)$ has a
higher value than the elastic energy density ${W}(u)$. These observations are
confirmed by the comparison of our thermal fracturing phase field models.

|
---|---
Figure 7: Profiles of $\mathcal{W}(\mathcal{A})$ (left) and
$\mathcal{W}^{*}(\mathcal{A})$ (right) in subdomain $\mathcal{A}$ during the
loading process.

##  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 [34, 24], we introduce fracturing PFM (we call it
F-PFM) in this section. Let $\Omega$ be a bounded (uncracked) domain in
$\mathbb{R}^{d}$ and
$\Gamma:=\partial\Omega=\Gamma_{D}^{u}\cup\Gamma_{N}^{u}$, similar to Section
2. In F-PFM, a crack in $\Omega$ at time $t$ is described by a damage variable
$z(x,t)\in[0,1]$ for $x\in\overline{\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:

$\displaystyle-{\mbox{div}}\left((1-z)^{2}{\sigma}[{u}]\right)=0$
$\displaystyle\mbox{in}~{}{\Omega}\times[0,{T}],$ (3.1a)
$\displaystyle{\alpha}\frac{\partial
z}{\partial{t}}=\left({\epsilon}~{}{\mbox{div}}\left({\gamma_{*}}{\nabla}z\right)-\frac{{\gamma_{*}}}{{\epsilon}}z+(1-z){W}(u)\right)_{+}$
$\displaystyle\mbox{in}~{}{\Omega}\times[0,{T}],$ (3.1b)

with the following boundary and initial conditions:

$\displaystyle u=u_{D}(x,t)$
$\displaystyle\mbox{on}~{}\Gamma_{D}^{u}\times[0,{T}],$ (3.2a)
$\displaystyle\sigma[u]n=0$
$\displaystyle\mbox{on}~{}\Gamma_{N}^{u}\times[0,{T}],$ (3.2b)
$\displaystyle\frac{\partial z}{\partial n}=0$
$\displaystyle\mbox{on}~{}\Gamma\times[0,{T}],$ (3.2c) $\displaystyle
z({x,0})=z_{*}(x)$ $\displaystyle\mbox{in}~{}\Omega,$ (3.2d)

where the displacement $u:\overline{\Omega}\times[0,T]\mapsto\mathbb{R}^{d}$
and the damage variable $z:\overline{\Omega}\times[0,T]\mapsto[0,1]$ are
unknowns. The parameters $\alpha>0$ and $\epsilon>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
$\text{G}_{\text{c}}$), and the elastic energy density is defined by
${W}={W}(u):=\sigma[u]:e[u]$. In (3.1b), the term ${W}$ works as a driving
force for $z$.

|
---|---

Figure 8: Illustration of the phase field approximation of the cracked surface
in an elastic body.

The symbol $(~{})_{+}$ on the right-hand side in (3.1b) denoted the positive
part $(s)_{+}:=\text{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}(z)$, where

$\displaystyle\mathcal{E}_{el}(u,z)$
$\displaystyle:=\frac{1}{2}\int_{\Omega}(1-z)^{2}\sigma[u]:e[u]~{}dx,$ (3.3)
$\displaystyle E_{s}(z)$
$\displaystyle:=\frac{1}{2}\int_{\Omega}\gamma_{*}\left(\epsilon\big{|}\nabla
z\big{|}^{2}+\frac{|z|^{2}}{\epsilon}\right)dx.$ (3.4)

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

$u(t)=\operatorname*{argmin}_{u\in V(u_{D}(t))}{\mathcal{E}_{el}(u,z(t))},$
(3.5)

and (3.1b) becomes a gradient flow of the energy
$\min\limits_{u}\mathcal{E}_{el}(u,z)+E_{s}(z)$.

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}_{\text{Y}}=(1-z)^{2}E_{\text{Y}}$. The energy $E_{s}(z)$ is
regularized surface energy, which approximates the crack area ($d=3$) or
length ($d=2$) as $\epsilon\mapsto 0$. Please see [24] for more details. The
following energy equality for F-PFM is shown in [24] ([34] for the antiplane
setting).

######  Theorem 3.1 (Energy equality for F-PFM)

Let $(u(x,t),z(x,t))$ be a sufficiently smooth solution to (3.1) and (3.2). If
$u_{D}$ is independent of $t$, then we have

$\displaystyle\frac{d}{dt}\left(\mathcal{E}_{el}(u(t),z(t))+E_{s}(z(t))\right)=-\alpha\int_{\Omega}\left|\dot{z}\right|^{2}dx\leq
0.$ (3.6)

###### Proof

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

$\displaystyle\frac{d}{dt}\left(\mathcal{E}_{el}(u(t),z(t))+E_{s}(z(t))\right)$
$\displaystyle=\int_{\Omega}(1-z)^{2}\sigma[u]:{e}[\dot{u}]~{}dx+\int_{\Omega}\left(\gamma_{*}\epsilon\nabla
z\cdot\nabla\dot{z}+\big{(}\frac{\gamma_{*}}{\epsilon}z-(1-z){W}(u)\big{)}\dot{z}\right)~{}dx$
$\displaystyle=\int_{\Gamma}(1-z)^{2}\underbrace{\left(\sigma[u]n\right)}_{0}\cdot\dot{u}~{}ds-\int_{\Omega}\underbrace{\mbox{div}\left((1-z)^{2}\sigma[u]\right)}_{0}\cdot\dot{u}~{}dx$
$\displaystyle\quad+\int_{\Gamma}\gamma_{*}\epsilon\underbrace{\frac{\partial
z}{\partial n}}_{0}\dot{z}~{}ds-\int_{\Omega}\mathcal{H}\dot{z}~{}dx,$ (3.7)

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

$\displaystyle\frac{d}{dt}\left(\mathcal{E}_{el}(u(t),z(t))+E_{s}(z(t))\right)$
$\displaystyle\qquad=-\int_{\Omega}\mathcal{H}\dot{z}~{}dx=-\int_{\Omega}\mathcal{H}\frac{(\mathcal{H})_{+}}{\alpha}~{}dx=-\int_{\Omega}\frac{(\mathcal{H})_{+}^{2}}{\alpha}~{}dx=-\int_{\Omega}\alpha\left|\dot{z}\right|^{2}dx.\qquad\text{\qed}$

###  3.2 Thermal fracturing phase field model 1 (TF-PFM1)

To combine the Biot model in (2.1) and F-PFM in (3.1), their variational
principles for $u$, Proposition 2.3 and (3.5), suggest that we consider the
following modified thermoelastic energy:

$\displaystyle\mathcal{E}_{el}^{*}(u,\Theta,z)$ $\displaystyle:=$
$\displaystyle\frac{1}{2}\int_{\Omega}(1-z)^{2}{\sigma}^{*}[u,\Theta]:{e}^{*}[u,\Theta]~{}dx,$
(3.8)

and a variational principle:

$u(t)=\operatorname*{argmin}_{u\in
V(u_{D}(t))}{\mathcal{E}_{el}^{*}(u,\Theta(t),z(t))}.$ (3.9)

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

For heat equation (2.1b), 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(z)>0$ in this section, and we set
it as $\kappa(z)=(1-z)^{2}\kappa_{0}$ in Section 4.

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

$\displaystyle-{\mbox{div}}\left((1-z)^{2}{\sigma}^{*}[{u},\Theta]\right)=0$
$\displaystyle\mbox{in}~{}{\Omega}\times[0,{T}],$ (3.10a)
$\displaystyle{\alpha}\frac{\partial
z}{\partial{t}}=\left({\epsilon}~{}{\mbox{div}}({\gamma_{*}}{\nabla}z)-\frac{{\gamma_{*}}}{{\epsilon}}z+(1-z){W}^{*}(u,\Theta)\right)_{+}$
$\displaystyle\mbox{in}~{}{\Omega}\times[0,{T}],$ (3.10b)
$\displaystyle\chi\frac{\partial{\Theta}}{\partial{t}}=~{}{\mbox{div}}\left({\kappa(z)}{\nabla}{\Theta}\right)-\Theta_{0}(1-z)^{2}\beta\frac{\partial}{\partial{t}}({\mbox{div}}{u})$
$\displaystyle\mbox{in}~{}{\Omega}\times(0,{T}],$ (3.10c)

Similar to (2.1) and (3.1), the boundary and the initial conditions to solve
(3.10) are presented as follows:

$\displaystyle u=u_{D}(x,t)$
$\displaystyle\mbox{on}~{}\Gamma_{D}^{u}\times[0,{T}],$ (3.11a)
$\displaystyle\sigma^{*}[u,\Theta]n=0$
$\displaystyle\mbox{on}~{}\Gamma_{N}^{u}\times[0,{T}],$ (3.11b)
$\displaystyle\Theta={\Theta}_{D}(x,t)$
$\displaystyle\mbox{on}~{}{\Gamma}_{D}^{\Theta}\times[0,{T}],$ (3.11c)
$\displaystyle\frac{\partial\Theta}{\partial n}=0$
$\displaystyle\mbox{on}~{}{\Gamma}_{N}^{\Theta}\times[0,{T}],$ (3.11d)
$\displaystyle\frac{\partial z}{\partial n}=0$
$\displaystyle\mbox{on}~{}\Gamma\times[0,{T}],$ (3.11e) $\displaystyle
z(x,0)=z_{*}(x)$ $\displaystyle\mbox{in}~{}\Omega,$ (3.11f)
$\displaystyle\Theta(x,0)=\Theta_{*}(x)$ $\displaystyle\mbox{in}~{}\Omega.$
(3.11g)

In the following, for simplicity, we define

$\displaystyle\sigma_{z}^{*}[u,\Theta]:=(1-z)^{2}\sigma^{*}[u,\Theta].$

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

######  Theorem 3.2 (Energy equality for TF-PFM1)

We suppose that $u_{D}\in H^{\frac{1}{2}}(\Gamma_{D}^{u};\mathbb{R}^{2})$ and
$\Theta\in L^{2}(\Theta)$ are given and do not depend on $t$. If $u(x,t)$ and
$z(u,t)$ are sufficiency smooth and satisfy (3.10a), (3.10b), (3.11a),
(3.11b), (3.11e), and (3.11f), the following energy equality holds:

$\displaystyle\frac{d}{dt}\left(\mathcal{E}_{el}^{*}(u(t),\Theta,z(t))+E_{s}(z(t))\right)=-\alpha\int_{\Omega}\left|{\dot{z}}\right|^{2}dx\leq
0.$ (3.12)

###### Proof

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

$\displaystyle\frac{d}{dt}\left(\mathcal{E}_{el}^{*}(u(t),\Theta,z(t))+E_{s}(z(t))\right)$
$\displaystyle=\frac{1}{2}\frac{d}{dt}\int_{\Omega}\left(\sigma^{*}_{z}[u,\Theta]:e^{*}[u,\Theta]\right)dx+\frac{1}{2}\frac{d}{dt}\int_{\Omega}\gamma_{*}\left(\epsilon\big{|}\nabla
z\big{|}^{2}+\frac{|z|^{2}}{\epsilon}\right)~{}dx$
$\displaystyle=\int_{\Omega}\sigma_{z}^{*}[u,\Theta]:e[\dot{u}]~{}dx+\int_{\Omega}\left(\gamma_{*}\epsilon\nabla
z\cdot\nabla\dot{z}+\big{(}\frac{\gamma_{*}}{\epsilon}z-(1-z){W}^{*}(u,\Theta)\big{)}\dot{z}\right)dx$
$\displaystyle=\int_{\Gamma}\underbrace{\sigma_{z}^{*}[u,\Theta]n}_{0}\cdot
e[\dot{u}]~{}ds-\int_{\Omega}\underbrace{\mbox{div}\sigma_{z}^{*}[u,\Theta]}_{0}\cdot
e[\dot{u}]~{}dx+\gamma_{*}\epsilon\int_{\Gamma}\underbrace{\frac{\partial
z}{\partial n}}_{0}\dot{z}~{}ds$
$\displaystyle\qquad\qquad-\int_{\Omega}\mathcal{H}^{*}\dot{z}~{}dx,$ (3.13)

where we also define $\mathcal{H}^{*}:=\epsilon\mbox{div}(\gamma_{*}\nabla
z)-\frac{\gamma_{*}}{\epsilon}z+(1-z){W}^{*}(u,\Theta)$. Since (3.10b) is
changed to $\alpha\dot{z}=(\mathcal{H}^{*})_{+}$, similar to that in Section
3.1, we conclude that

$\displaystyle\frac{d}{dt}\left(\mathcal{E}_{el}^{*}(u(t),\Theta,z(t))+E_{s}(z(t))\right)=-\alpha\int_{\Omega}\left|\dot{z}\right|^{2}~{}dx\leq
0,$

which is equivalent to (3.12).

###  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,\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,\Theta)$ and $E_{el}(u)$. The variational
principle holds for $E_{el}^{*}(u,\Theta)$ (Proposition 2.3), and the energy
equality holds for $E_{el}(u)$ (Theorem 2.4). This motivates us to consider
another type of thermal fracturing PFM based on elastic energy $E_{el}(u)$. We
call the following thermal fracturing model TF-PFM2:

$\displaystyle-{\mbox{div}}\left((1-z)^{2}{\sigma}^{*}[{u},{\Theta}]\right)=0$
$\displaystyle\mbox{in}~{}{\Omega}\times[0,{T}],$ (3.14a)
$\displaystyle{\alpha}\frac{\partial
z}{\partial{t}}=\left({\epsilon}~{}{\mbox{div}}({\gamma_{*}}{\nabla}z)-\frac{{\gamma_{*}}}{{\epsilon}}z+(1-z){W}(u)\right)_{+}$
$\displaystyle\mbox{in}~{}{\Omega}\times[0,{T}],$ (3.14b)
$\displaystyle\chi\frac{\partial{\Theta}}{\partial{t}}=~{}{\mbox{div}}\left({\kappa(z)}{\nabla}{\Theta}\right)-\Theta_{0}(1-z)^{2}\beta\frac{\partial}{\partial{t}}({\mbox{div}}{u})$
$\displaystyle\mbox{in}~{}{\Omega}\times(0,{T}].$ (3.14c)

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

######  Theorem 3.3 (Energy equality for TF-PFM2)

We suppose that $(u(x,t),$ $\Theta(x,t),z(x,t))$ is a sufficiently smooth
solution for (3.14) and (3.11). If $u_{D}$ is independent of $t$ and
$\Theta_{D}=\Theta_{0}$, then the following energy equality holds:

$\displaystyle\frac{d}{dt}\left(\mathcal{E}_{el}(u(t),z(t))+E_{s}(z(t))+E_{th}(\Theta(t))\right)$
$\displaystyle\hskip
56.9055pt=-\frac{1}{\Theta_{0}}\int_{\Omega}\kappa(z)\left|\nabla\Theta\right|^{2}~{}dx-\alpha\int_{\Omega}\left|\dot{z}\right|^{2}dx\leq
0.$ (3.15)

###### Proof

Since the relation in (2.30) is written as

$\displaystyle\frac{d}{dt}\left(\frac{1}{2}W(u)\right)=\sigma^{*}[u,\Theta]-\beta(\Theta-\Theta_{0})\mbox{div}\dot{u},$

we obtain

$\displaystyle\frac{d}{dt}\left(\frac{1}{2}(1-z)^{2}W(u)\right)=\sigma_{z}^{*}[u,\Theta]:e[\dot{u}]+\beta(1-z)^{2}(\Theta-\Theta_{0})\mbox{div}\dot{u}-(1-z)\dot{z}W(u).$

Hence, we have

$\displaystyle\frac{d}{dt}\mathcal{E}_{el}(u(t),z(t))+\frac{d}{dt}E_{s}(z(t))$
$\displaystyle\quad\quad=\int_{\Omega}\frac{d}{dt}\left(\frac{1}{2}(1-z)^{2}W(u)\right)~{}dx+\int_{\Omega}\left(\epsilon\mbox{div}(\gamma_{*}\nabla
z)-\frac{\gamma_{*}}{\epsilon}z\right)\dot{z}~{}dx$
$\displaystyle\quad\quad=\underbrace{\int_{\Omega}\sigma_{z}^{*}[u,\Theta]:e[\dot{u}]~{}dx}_{0}+\int_{\Omega}\beta(1-z)^{2}(\Theta-\Theta_{0})\mbox{div}\dot{u}~{}dx-\int_{\Omega}\mathcal{H}\dot{z}~{}dx$
$\displaystyle\quad\quad=\int_{\Omega}\beta(1-z)^{2}(\Theta-\Theta_{0})\mbox{div}\dot{u}~{}dx-\int_{\Omega}\alpha|\dot{z}|^{2}~{}dx,$
(3.16)

where $\mathcal{H}=\epsilon\mbox{div}(\gamma_{*}\nabla
z)-\frac{\gamma_{*}}{\epsilon}z+(1-z)W(u)$.

On the other hand,

$\displaystyle\frac{d}{dt}E_{th}(\Theta(t))=\frac{\chi}{\Theta_{0}}\int_{\Omega}(\Theta-\Theta_{0})\dot{\Theta}~{}dx$
$\displaystyle\quad\quad\quad=\frac{1}{\Theta_{0}}\int_{\Omega}(\Theta-\Theta_{0})\left\\{\mbox{div}(\kappa(z)\nabla\Theta)-\Theta_{0}\beta(1-z)^{2}\mbox{div}\dot{u}\right\\}~{}dx$
$\displaystyle\quad\quad\quad=-\frac{1}{\Theta_{0}}\int_{\Omega}\kappa(z)|\nabla\Theta|^{2}~{}dx-\int_{\Omega}\beta(1-z)^{2}(\Theta-\Theta_{0})\mbox{div}\dot{u}~{}dx.$
(3.17)

Taking a sum of these equalities (3.16)-(3.17), we obtain the energy equality
(3.15).

##  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 Nondimensional setting

In the following numerical examples, we suppose
$\kappa(z)=(1-z)^{2}\kappa_{0}$. For convenience, we consider the
nondimensional form with (2.34), (2.35), (2.37), and

$\displaystyle\tilde{\epsilon}=\frac{\epsilon}{c_{x}},~{}\tilde{\gamma}_{*}=\frac{{c_{e}\gamma_{*}}}{c_{x}(\beta
c_{\Theta})^{2}},~{}\tilde{\alpha}=\frac{c_{e}\alpha}{c_{t}(\beta
c_{\Theta})^{2}},~{}\tilde{a}_{L}=\frac{c_{x}c_{\Theta}}{c_{u}}{a}_{L},~{}\tilde{\beta}=1.$

Then, TF-PFM1 in (3.10) is expressed in the following nondimensional form:

$\displaystyle-{\mbox{div}}\left((1-z)^{2}{\sigma}[{u}]\right)=(1-z)^{2}\nabla\Theta$
$\displaystyle\mbox{in}~{}{\Omega}\times[0,{T}],$ (4.1a)
$\displaystyle{\alpha}\frac{\partial
z}{\partial{t}}=\left({\epsilon}~{}{\mbox{div}}({\gamma_{*}}{\nabla}z)-\frac{{\gamma_{*}}}{{\epsilon}}z+(1-z){W}^{*}(u,\Theta)\right)_{+}$
$\displaystyle\mbox{in}~{}{\Omega}\times[0,{T}],$ (4.1b)
$\displaystyle\frac{\partial{\Theta}}{\partial{t}}=~{}{\mbox{div}}\left((1-z)^{2}{\nabla}{\Theta}\right)-(1-z)^{2}\delta\frac{\partial}{\partial{t}}({\mbox{div}}{u})$
$\displaystyle\mbox{in}~{}{\Omega}\times(0,{T}].$ (4.1c)

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

$\displaystyle{\alpha}\frac{\partial
z}{\partial{t}}=\left({\epsilon}~{}{\mbox{div}}({\gamma_{*}}{\nabla}z)-\frac{{\gamma_{*}}}{{\epsilon}}z+(1-z){W}(u)\right)_{+}$
$\displaystyle\mbox{in}~{}{\Omega}\times[0,{T}].$ (4.2)

###  4.2 Time discretization

To solve problem (4.1), we adopt the following semi-implicit time
discretization scheme [24, 34].

$\displaystyle-{\mbox{div}}\left((1-z^{k-1})^{2}{\sigma}^{*}[{u}^{k},{\Theta}^{k-1}]\right)=0,$
(4.3a)
$\displaystyle{\alpha}\frac{\tilde{z}^{k}-z^{k-1}}{\Delta{t}}={\epsilon}~{}{\mbox{div}}\left({\gamma_{*}}{\nabla}\tilde{z}^{k}\right)-\frac{{\gamma_{*}}}{{\epsilon}}\tilde{z}^{k}+\left(1-\tilde{z}^{k}\right){W}^{*}({u}^{k-1},{\Theta}^{k-1}),$
(4.3b) $\displaystyle z^{k}:=\mbox{max}\left(\tilde{z}^{k},z^{k-1}\right),$
(4.3c)
$\displaystyle\dfrac{{\Theta}^{k}-{\Theta}^{k-1}}{\Delta{t}}=~{}{\mbox{div}}\left((1-z^{k-1}){\nabla}{\Theta}^{k}\right)-(1-z^{k-1})\delta{\mbox{div}}\Big{(}\dfrac{{u}^{k}-{u}^{k-1}}{\Delta{t}}\Big{)}.$
(4.3d)

For TF-PFM2, (4.3b) is replaced by

$\displaystyle{\alpha}\frac{\tilde{z}^{k}-z^{k-1}}{\Delta{t}}={\epsilon}~{}{\mbox{div}}\left({\gamma_{*}}{\nabla}\tilde{z}^{k}\right)-\frac{{\gamma_{*}}}{{\epsilon}}\tilde{z}^{k}+\left(1-\tilde{z}^{k}\right){W}({u}^{k-1}),$
(4.4)

where $u^{k}$, $z^{k}$, and $\Theta^{k}$ are the approximations of $u$, $z$,
$\Theta$, respectively, at time ${t}_{k}:=k\Delta t(k=1,2,3,\cdots)$. Since
the adaptive mesh technique in the FEM is often effective and accurate in
numerical experiments with phase field models, problems (4.3) - (4.4) 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 [17] and
executed on a desktop with an Intel(R) Core i7$-$7820X<EMAIL_ADDRESS>GHz, 16 core
processor, and 64 GB RAM.

|
---|---
Figure 9: Domain for Section 4.3 with $z_{*}(x)$ as the initial crack (left)
and the adaptive mesh for the initial crack (right).

###  4.3 Thermoelastic effect on the crack speed

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

Table 4: List of the nondimensional parameters for Sections 4.3 and 4.4 Parameter |  $E_{Y}$ |  $\nu_{P}$ |  $a_{L}$ |  $\alpha$ |  $\epsilon$ |  $\gamma_{*}$ |  $\Theta_{*}$
---|---|---|---|---|---|---|---
Value |  1 |  0.3 |  0.7 |  0.001 |  0.01 |  5.08 |  0

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

---
Figure 10: Snapshots of crack propagation with F-PFM, TF-PFM1, and TF-PFM2 in
$(-1,1)\times(-0.35,0.35)$ at $t=0.4,~{}0.6,~{}0.8,~{}1$ (left to right). For
TF-PFM1 and TF-PFM2, we use the thermoelasticity coupling parameter
$\delta=0.5$, and the color represents the value of $z$.

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 Figures 10 and 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-PFM.

The temperature distributions obtained by TF-PFM1 and TF-PFM2 are shown in
Figure 12. In the equation for $\Theta$, the heat resource is given by
$-(1-z)^{2}\delta\frac{d}{dt}(\mbox{div}u)$. During crack propagation
$(0.4\leq t\leq 0.8)$, the areas near the crack tip, the upper-right corner,
and lower-right corner are continuously expanding when $\mbox{div}u>0$ and
$\frac{\partial}{\partial t}(\mbox{div}u)>0$. Therefore, due to the negative
source $-\frac{\partial}{\partial t}(\mbox{div}u)$, 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.

|
---|---
|
Figure 11: Comparison of the profiles of $z$ obtained by F-PFM, TF-PFM1, and TF-PFM2 along the line $x_{2}=0$ at (a) $t=0.4$, (b) $t=0.6$, (c) $t=0.8$, and (d) $t=1$. |  |  |
---|---|---|---
|  |  |
Figure 12: Snapshots of the temperatures obtained by TF-PFM1 (upper) and TF-
PFM2 (lower) at $t=0.4,~{}0.6,~{}0.8,~{}1$ (left to right); the color
represents the value of $\Theta$.

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.

|
---|---
|
Figure 13: Profile of the elastic (left) and surface energy (right) under thermal expansion during crack propagation using TF-PFM1 (top) and TF-PFM2 (bottom). |
---|---

Figure 14: Mode I (left) and Mode I+II (right) for the study of the crack path
under thermal expansion and the loading process. Here, the initial damage
$z_{*}(x)$ is illustrated by the red initial crack in the figures.

###  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:

$\displaystyle C_{\pm}:=\left(\begin{array}[]{l}-\frac{1}{2}\\\
\pm\frac{5}{8}\end{array}\right)\in\mathbb{R}^{2},$ (4.7) $\displaystyle
H_{\pm}:=\Big{\\{}x\in\mathbb{R}^{2};~{}\left|x-C_{\pm}\right|\leq\frac{3}{20}\Big{\\}},$
$\displaystyle\Omega:=(-1,1)^{2}\setminus(H_{+}\cup H_{-}),$

and we define

$\displaystyle\Gamma_{DN1}^{u}:=\Gamma\cap\\{x_{1}=1\\},~{}\Gamma_{DN2}^{u}:=\partial
H_{+}\cup\partial
H_{-},~{}\Gamma_{N}^{u}:=\Gamma\setminus(\Gamma_{DN1}^{u}\cup\Gamma_{DN2}^{u}),$
$\displaystyle\Gamma_{\pm D}^{\Theta}:=\Gamma\cap\\{x_{2}=\pm
1\\},~{}\Gamma_{N}^{\Theta}:=\Gamma\setminus(\Gamma_{+D}^{\Theta}\cup\Gamma_{-D}^{\Theta}).$

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

$\displaystyle\left\\{\begin{array}[]{l}u_{1}=0\\\
\sigma_{12}^{*}=0\end{array}\right.~{}\mbox{on}~{}\Gamma_{DN1}^{u},\quad\left\\{\begin{array}[]{l}(\sigma^{*}n)\cdot
e_{1}=0\\\ u_{2}=\pm 8t\end{array}\right.~{}\mbox{on}~{}{\partial
H_{\pm}},\quad\sigma^{*}[u,\Theta]n=0~{}\mbox{on}~{}\Gamma_{N}^{u},$ (4.12)
$\displaystyle\Theta=\Theta_{D}~{}\mbox{on}~{}\Gamma_{+D}^{\Theta},\quad\Theta=0~{}\mbox{on}~{}\Gamma_{-D}^{\Theta},\quad\frac{\partial\Theta}{\partial
n}=0~{}\mbox{on}~{}\Gamma_{N}^{\Theta}.$

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

For $z$, similar to the previous example (Section 4.3) , we set
$\frac{\partial z}{\partial n}=0$ on $\Gamma$ and choose the initial value as
$z_{*}(x):=\exp{(-(x_{2}/\eta)^{2})}/(1+\exp{((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
$\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 [22].

|  |  |
---|---|---|---
|  |  |
|  |  |
Figure 15: Snapshots of the crack paths. F-PFM (upper), TF-PFM1 (middle), and TF-PFM2 (lower) at $t=0.4,~{}0.6,~{}0.8,~{}1$ (left to right). For TF-PFM1 and TF-PFM2, we set $\Theta_{D}=10$ and $\delta=0.5$. Here, the color represents the value of $z$. |  |  |
---|---|---|---
|  |  |
Figure 16: Snapshots of the temperature gradient during thermal expansion and
crack growth under the given temperature $\Theta_{D}=10$. TF-PFM1 (top) and
TF-PFM2 (bottom) at $t=0.4,~{}0.6,~{}0.8,~{}1$ (left to right); the color
represents the value of $\Theta$.

Figure 17 shows the crack paths for different temperature gradients
$\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.

|
---|---
Figure 17: Comparison of the crack paths using TF-PFM1 (left) and TF-PFM2
(right) with the given temperature variations under Mode I at the final
computational time $t=1$.

The temperature distributions during crack growth are shown in Figure 16.
There exists a temperature discontinuity along the crack path, which is caused
by $\kappa(z)=(1-z)^{2}\kappa_{0}$. It approximately represents a thermal
insulation condition across the crack.

####  4.4.2 Mode I+II

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

$\displaystyle\Gamma_{\pm D}^{u}:=\Gamma\cap\\{x_{2}=\pm
1\\},~{}\Gamma_{N}^{u}:=\Gamma\setminus(\Gamma_{+D}^{u}\cup\Gamma_{-D}^{u}),$
$\displaystyle\Gamma_{\pm D}^{\Theta}:=\Gamma\cap\\{x_{2}=\pm
1\\},,~{}\Gamma_{N}^{\Theta}:=\Gamma\setminus(\Gamma_{+D}^{\Theta}\cup\Gamma_{-D}^{\Theta}).$

The boundary conditions for $u$ are given as follows:

$\displaystyle\left\\{\begin{array}[]{ll}u_{1}=\pm 3\sin(\pi/3)t,\\\ u_{2}=\pm
3\cos(\pi/3)t\end{array}\right.~{}\mbox{on}~{}\Gamma_{\pm
D}^{u},\qquad\sigma^{*}[u,\Theta]n=0~{}\mbox{on}~{}\Gamma_{N}^{u}.$ (4.15)

The boundary conditions for $\Theta$ and $z$ are the same as those in Section
4.4.1. The initial crack profile is given as
$z_{*}(x):=\exp{(-(x_{2}/\eta)^{2})}/(1+\exp{((x_{1}-0.5)/\eta)})-\exp{(-(x_{2}/\eta)^{2})}/(1+\exp{((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
$\Theta_{D}=0,~{}2,~{}3,~{}5,~{}6$.

|
---|---
Figure 18: Comparison of the crack paths using TF-PFM1 (left) and TF-PFM2
(right) with the given temperature variations under Mode I+II at the final
computational time.

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
$\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.

|
---|---
Figure 19: Comparison of the crack paths using TF-PFM1 and TF-PFM2 when
$\Theta=5$ (left) and $\Theta=6$ (right) at the final computational time.

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.

##  5 Summary and Future Works

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

For the Biot model, we studied a variational principle (Proposition 2.3) and
energy equality (Theorem 2.4), which were related to different energies
$E^{*}_{el}(u,\Theta)$ and $E_{el}(u)+E_{th}(\Theta)$, respectively (see
Tables 1 and 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}(z)$ and admits energy equality
(Theorem 3.1).

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,\Theta)$ in
TF-PFM1 (3.10b) and $W(u)$ in TF-PFM2 (3.14b).

Consequently, we established partial energy equality for TF-PFM1 (Theorem 3.2)
and energy equality for TF-PFM2 (Theorem 3.3). 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.

  1. 1.

The thermoelastic coupling parameter $\delta$ in TF-PFM1 and TF-PFM2 enhances
crack propagation (Figure 10).

  2. 2.

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 17 \- 19).

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

Table 5: Numerical comparison of TF-PFM1 and TF-PFM2 Models |  Driving force |  Energy consistency |  Straight crack speed |  Crack path
---|---|---|---|---
TF-PFM1 |  $W^{*}(u,\Theta)=\sigma^{*}[u,\Theta]:e^{*}[u,\Theta]$ |  Partially satisfied |  Faster |  Less curved
TF-PFM2 |  $W(u)=\sigma[u]:e[u]$ |  Fully satisfied |  Slower |  More curved
Remarks |  $W^{*}(u,\Theta)>W(u)$ (Figure 4) |  Theorems 3.2, & 3.3 |  Figure 13 |  Figures 17 & 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 [3, 24] should be introduced in our
PFM.

##  Appendix A Weak forms for (2.40)

The implicitly time-discretized problem (2.40) is solved with the following
boundary conditions. We set the initial temperature $\Theta^{0}:=\Theta_{*}$
and set $\Theta^{-1}=\Theta_{*}$, which is a temperature of $t=-\Delta t$. For
a given $\Theta^{k-1}$, the boundary value problem of $u^{k}$ is given as
follows:

$\displaystyle\left\\{\begin{array}[]{ll}-{{\text{div}}{\sigma}^{*}[{u}^{k},\Theta^{k-1}]=0}&\mbox{in}~{}\Omega,\\\
u^{k}=u_{D}(\cdot,t_{k})&\text{on}~{}\Gamma_{D}^{u},\\\
\sigma^{*}[u^{k},\Theta^{k-1}]n=0&\text{on}~{}\Gamma_{N}^{u},\end{array}\right.\qquad(k=0,1,2,\cdots).$
(1.4)

We define a weak form for (1.4) as

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle u^{k}\in
V^{u}(u_{D}(\cdot,t_{k})),\\\
\displaystyle\int_{\Omega}\sigma^{*}[u^{k},\Theta^{k-1}]:e[v]~{}dx=0,\qquad(\mbox{for
all}~{}v\in V^{u}(0)),\end{array}\right.$ (1.7)

where $V^{u}(\cdot)$ is defined by (2.15). The second equation of (1.7) is
equivalent to

$\displaystyle\int_{\Omega}\sigma[u^{k}]:e[v]~{}dx=\int_{\Omega}\Theta^{k-1}\mbox{div}v~{}dx.$
(1.8)

Similarly, for given $u^{k-1}$ and $u^{k}$, the boundary value problem of
$\Theta^{k}$ is given as follows:

$\displaystyle{\left\\{\begin{array}[]{ll}\displaystyle\frac{\Theta^{k}-\Theta^{k-1}}{\Delta
t}-\Delta\Theta^{k}+\delta\mbox{div}\left(\frac{u^{k}-u^{k-1}}{\Delta
t}\right)=0&\mbox{in}~{}\Omega,\\\
\displaystyle\Theta^{k}=0&\text{on}~{}\Gamma_{D}^{\Theta},\\\
\displaystyle\frac{\partial\Theta^{k}}{\partial
n}=0&\text{on}~{}\Gamma_{N}^{\Theta},\end{array}\right.~{}(k=1,2,\cdots)}.$
(1.12)

We define a weak form for 1.12 as

$\displaystyle\left\\{\begin{array}[]{ll}\displaystyle\Theta^{k}\in
V^{\Theta},\\\
\displaystyle\int_{\Omega}\left(\frac{\Theta^{k}-\Theta^{k-1}}{\Delta
t}\right)\psi~{}dx+\int_{\Omega}\nabla\Theta^{k}\cdot\nabla\psi~{}dx\\\
\displaystyle\hskip
50.0pt+\delta\int_{\Omega}\mbox{div}\left(\frac{u^{k}-u^{k-1}}{\Delta
t}\right)\psi~{}dx=0~{}(\text{for all}~{}\psi\in
V^{\Theta}(0)),\end{array}\right.$ (1.16)

where $V^{\Theta}:=\\{\psi\in
H^{1}(\Omega);~{}\psi|_{\Gamma_{D}^{\Theta}}=0\\}$.

###### Proposition 1

We suppose that the $(d-1)$-dimensional volume of $\Gamma_{D}^{u}$ is
positive. If $\Theta^{0}=\Theta_{*}\in L^{2}(\Omega)$ and
$u_{D}(\cdot,t_{k})\in H^{\frac{1}{2}}(\Gamma_{D}^{u})$ $(k=0,1,2,\cdots)$,
then weak solutions $u^{k}$ $(k=0,1,2,\cdots)$ for (1.7) and $\Theta^{k}$
$(k=1,2,\cdots)$ for (1.16) uniquely exist.

###### Proof

At each time step, the unique solvabilities of (1.7) and (1.16) follow from
the Lax-Milgram theorem [19, 12]. More precisely, first we solve $u^{0}$ by
(1.7). Then, for $k=1,2,\cdots$, we can obtain $u^{k}$ by (1.7) and
$\Theta^{k}$ by (1.16), sequentially.

##  Appendix B Divergence of $u$ around the crack tip

We want to observe the contracting and expanding areas around the crack tip
area. Here, we show an analytical solution for $\mbox{div}u$ around the crack
tip. We consider Mode I as the type of loading; then, we analytically obtain
the following crack tip displacement field:

$\displaystyle
u_{1}=\frac{K_{I}}{2\mu}\sqrt{\frac{r}{2\pi}}\cos\left(\frac{\theta}{2}\right)\left[\xi-1+2\sin^{2}\left(\frac{\theta}{2}\right)\right],$
(2.1) $\displaystyle
u_{2}=\frac{K_{I}}{2\mu}\sqrt{\frac{r}{2\pi}}\sin\left(\frac{\theta}{2}\right)\left[\xi+1-2\cos^{2}\left(\frac{\theta}{2}\right)\right],$
(2.2)

where $K_{I}$, $\mu$, $\xi=3-4\nu_{P}$, and ($r$, $\theta$) are the Mode I
stress intensity factor, Lamé’s constant, plane strain and polar coordinates
for the crack tip, respectively.

Now, we can calculate $\mbox{div}u$ as follows:

$\displaystyle\mbox{div}u$ $\displaystyle=$ $\displaystyle\left(\frac{\partial
r}{\partial x_{1}}\frac{\partial}{\partial r}+\frac{\partial\theta}{\partial
x_{1}}\frac{\partial}{\partial\theta}\right)u_{1}+\left(\frac{\partial
r}{\partial x_{2}}\frac{\partial}{\partial r}+\frac{\partial\theta}{\partial
x_{2}}\frac{\partial}{\partial\theta}\right)u_{2}$ (2.3) $\displaystyle=$
$\displaystyle\frac{K_{I}}{2\mu\sqrt{2\pi
r}}\Big{(}\xi\cos\Big{(}\frac{\theta}{2}\Big{)}-\cos\Big{(}\frac{\theta}{2}\Big{)}\Big{)}=\frac{K_{I}(\xi-1)}{2\mu\sqrt{2\pi
r}}\cos\left(\frac{\theta}{2}\right).$

Assume a crack is growing as

$\displaystyle\Sigma(t)=\left\\{(x_{1},0)^{T}\big{|}~{}-\infty<x_{1}\leq
v_{0}t\right\\}.$

Then, we obtain the following displacement at time $t$

$\displaystyle\tilde{u}(x,t)\approx
u(x-v_{0}te_{1}),~{}\text{where}~{}e_{1}:=(1,0)^{T},$

and we also obtain $\mbox{div}u$ at time $t$

$\displaystyle\mbox{div}\tilde{u}(x,t)$
$\displaystyle=\mbox{div}u(x-v_{0}te_{1})$
$\displaystyle\frac{\partial}{\partial
t}\mbox{div}\tilde{u}(x,t)\big{|}_{t=0}$
$\displaystyle=-v_{0}\frac{\partial}{\partial
x_{1}}\mbox{div}u=\frac{v_{0}K_{I}(\xi-1)}{4\mu\sqrt{2\pi
r^{3}}}\cos\left(\frac{3\theta}{2}\right).$ (2.4)

Now, we set $E_{Y}=1$, $\nu_{P}=0.3$, $K_{I}=5$, and $v_{0}=0.05$, and then we
obtain the displacement, $\mbox{div}u$ and the $\frac{\partial}{\partial
x_{1}}\mbox{div}u$ profiles through (2.1), (2.2), (2.3), and (2.4),
respectively. From Figure 20, a compressing area exists at the crack tip.

Figure 20: Profile of displacement $[u_{1},u_{2}]$ (left), $\mbox{div}u$
(right), and $\frac{\partial}{\partial x_{1}}(\mbox{div}{u})$ (bottom) around
the crack under Mode I.

## Acknowledgments

This work was supported by the MEXT (the Ministry of Education, Culture,
Sports, Science, and Technology) scholarship in Japan. This work was also
partially supported by JSPS KAKENHI, grant numbers JP20H01812 and JP20KK0058.

## References

  * [1] Ai, W., and Augarde, C. E., 2019. Thermoelastic fracture modelling in 2D by an adaptive cracking particle method without enrichment functions. International Journal of Mechanical Sciences. 160, 343-357. https://doi.org/10.1016/j.ijmecsci.2019.06.033.
  * [2] Alfat, S., Kimura, M., Firihu, M. Z., and Rahmat, 2018. Numerical investigation of shape domain effect to its elasticity and surface energy using adaptive finite element method, AIP Conference Proceedings. 1964(1), 020011. https://doi.org/10.1063/1.5038293.
  * [3] Amor, H., Marigo, J.-J., and Maurini, C., 2009. Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments. J. Mech. Phys. Solids. 57, 1209-1229. https://doi.org/10.1016/j.jmps.2009.04.011.
  * [4] Anderson, T. L., 2017. Fracture Mechanics: Fundamentals and Applications, fourth ed. CRC press.
  * [5] Biot, M. A., 1956. Thermoelasticity and irreversible thermodynamics. Journal of applied physics. 27(3), 240-253. https://doi.org/10.1063/1.1722351.
  * [6] Bourdin, B., Francfort, G. A., and Marigo, J.-J., 2000. Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids. 48(4), 797-826. https://doi.org/10.1016/S0022-5096(99)00028-9.
  * [7] Bourdin, B., 2007. Numerical implementation of the variational formulation of brittle fracture. Interfaces Free Boundaries. 9, 411-430. doi: 10.4171/IFB/171.
  * [8] Bourdin, B., Francfort, G. A., and Marigo, J.-J., 2008. The Variational Approach to Fracture. Journal of Elasticity. 91, 5–148. https://doi.org/10.1007/s10659-007-9107-3.
  * [9] Bourdin, B., Marigo, J.-J., Maurini, C., and Sicsic, P., 2014. Morphogenesis and propagation of complex cracks induced by thermal shocks. Physical review letters. 112(1), 014301. https://doi.org/10.1103/PhysRevLett.112.014301.
  * [10] Ciarlet, P.G., 2002. The Finite Element Method for Elliptic Problems. Society for Industrial and Applied Mathematics.
  * [11] Duflot, M., 2008. The extended finite element method in thermoelastic fracture mechanics. International Journal for Numerical Methods in Engineering. 74(5), 827-847. https://doi.org/10.1002/nme.2197.
  * [12] Duvaut, G., and Lions, J. L., 1976. Inequalities in Mechanics and Physics. Springer, Berlin, Heidelberg.
  * [13] Dwivedi, S. K., and Vishwakarma, M., 2018. Hydrogen embrittlement in different materials: a review. International Journal of Hydrogen Energy. 43(46), 21603-21616. https://doi.org/10.1016/j.ijhydene.2018.09.201.
  * [14] Entezari, A., Filippi, M., Carrera, E., and Kouchakzadeh, M. A., 2018. 3D dynamic coupled thermoelastic solution for constant thickness disks using refined 1D finite element models. Applied Mathematical Modelling. 60, 273-285. https://doi.org/10.1016/j.apm.2018.03.015.
  * [15] Francfort, G. A., and Marigo, J.-J., 1998. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids. 46(8), 1319-1342. https://doi.org/10.1016/S0022-5096(98)00034-9.
  * [16] Freiman, S. W., 1984. Effects of chemical environments on slow crack growth in glasses and ceramics. Journal of Geophysical Research: Solid Earth. 89(B6), 4072-4076. https://doi.org/10.1029/JB089iB06p04072.
  * [17] Hecht, F., 2012. New development in FreeFem++. Journal of numerical mathematics. 20(3-4), 251-266. https://doi.org/10.1515/jnum-2012-0013.
  * [18] Gao, Y., and Oterkus, S., 2019. Ordinary state-based peridynamics modelling for fully coupled thermoelastic problems. Continuum Mechanics and Thermodynamics. 31(4), 907-937. https://doi.org/10.1007/s00161-018-0691-1.
  * [19] Girault, V., and Raviart, P. A., 1979. Finite element approximation of the Navier-Stokes equations. Lecture Notes in Mathematics, Berlin Springer Verlag, 749.
  * [20] Green, A. E., and Lindsay, K. A., 1972. Thermoelasticity. Journal of Elasticity. 2, 1-7. https://doi.org/10.1007/BF00045689.
  * [21] Green, A. E., and Naghdi, P. M., 1991. A re-examination of the basic postulates of thermomechanics. Proceedings of the Royal Society of London, Series A: Mathematical and Physical Sciences. 432(1885), 171-194. https://doi.org/10.1098/rspa.1991.0012.
  * [22] Jaskowiec, J., 2017. A model for heat transfer in cohesive cracks. Computers & Structures, 180, 89-103. https://doi.org/10.1016/j.compstruc.2016.01.009.
  * [23] Karma, A., Kessler, D. A., and Levine, H., 2001. Phase-field model of mode III dynamic fracture. Physical Review Letters. 87(4), 045501. https://doi.org/10.1103/PhysRevLett.87.045501.
  * [24] Kimura, M., Takaishi, T., Alfat, S., Nakano, T., and Tanaka, Y., Irreversible phase field models for crack growth in industrial applications: thermal stress, viscoelasticity, hydrogen embrittlement. SN Applied Sciences. 3(781). https://doi.org/10.1007/s42452-021-04593-6.
  * [25] Kouchakzadeh, M. A., and Entezari, A., 2015. Analytical solution of classic coupled thermoelasticity problem in a rotating disk. Journal of Thermal Stresses. 38(11), 1267-1289. https://doi.org/10.1080/01495739.2015.1073529.
  * [26] Kou, M., Liu, X., Tang, S., and Wang, Y., 2019. 3-D X-ray computed tomography on failure characteristics of rock-like materials under coupled hydro-mechanical loading. Theoretical and Applied Fracture Mechanics. 104, 102396. https://doi.org/10.1016/j.tafmec.2019.102396.
  * [27] Lord, H. W., and Shulman, Y., 1967. A generalized dynamical theory of thermoelasticity. Journal of the Mechanics and Physics of Solids, 15(5), 299-309. https://doi.org/10.1016/0022-5096(67)90024-5.
  * [28] Louthan Jr., M. R., Caskey Jr., G. R., Donovan, J. A., and Rawl Jr., D. E., 1972. Hydrogen embrittlement of metals. Materials Science and Engineering. 10, 357-368. https://doi.org/10.1016/0025-5416(72)90109-7.
  * [29] Mackin, T. J., et al., 2002. Thermal cracking in disc brakes. Engineering Failure Analysis. 9(1), 63-76. https://doi.org/10.1016/S1350-6307(00)00037-6.
  * [30] Miehe, C., Schaenzel, L. M., and Ulmer, H., 2015. Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids. Computer Methods in Applied Mechanics and Engineering. 294, 449-485. https://doi.org/10.1016/j.cma.2014.11.016.
  * [31] Nara, Y., Morimoto, K., Yoneda, T., Hiroyoshi, N., and Kaneko, K., 2011. Effects of humidity and temperature on subcritical crack growth in sandstone. International Journal of Solids and Structures. 48(7-8), 1130-1140. https://doi.org/10.1016/j.ijsolstr.2010.12.019.
  * [32] Nguyen, M. N., Bui, T. Q., Nguyen, N. T., and Truong, T. T., 2017. Simulation of dynamic and static thermoelastic fracture problems by extended nodal gradient finite elements. International Journal of Mechanical Sciences. 134, 370-386. https://doi.org/10.1016/j.ijmecsci.2017.10.022.
  * [33] Nguyen, T.T., Yvonnet, J., Bornert, M., et al., 2016. On the choice of parameters in the phase field method for simulating crack initiation with experimental validation. Int J Fract 197. 213–226. https://doi.org/10.1007/s10704-016-0082-1.
  * [34] Takaishi, T., and Kimura, M., 2009. Phase field model for mode III crack growth in two dimensional elasticity. Kybernetika. 45(4), 605-614.
  * [35] Vivekanandan, A., and Ramesh, K., 2020. Study of crack interaction effects under thermal loading by digital photoelasticity and finite elements. Experimental Mechanics, 60(3), 295-316. https://doi.org/10.1007/s11340-019-00561-9.
  * [36] Zheng, B. J., Gao, X. W., Yang, K., Zhang, C. Z., 2015. A novel meshless local Petrov–Galerkin method for dynamic coupled thermoelasticity analysis under thermal and mechanical shock loading. Engineering Analysis with Boundary Elements. 60, 154-161. https://doi.org/10.1016/j.enganabound.2014.12.001.
  * [37] Zhou, F. X., Li, S. R., and Lai, Y. M., 2011. Three-dimensional analysis for transient coupled thermoelastic response of a functionally graded rectangular plate. Journal of Sound and Vibration, 330(16), 3990-4001. https://doi.org/10.1016/j.jsv.2011.03.015.