# A phase-field model for thermo-mechanical fracture with an open-source
implementation of it using Gridap in Julia

Ved Prakash  School of Infrastructure
Indian Institute of Technology
Bhubaneswar
Odisha-752050, India
Email<EMAIL_ADDRESS>Akash Kumar Behera  School of Infrastructure
Indian Institute of Technology
Bhubaneswar
Odisha-752050, India
Email<EMAIL_ADDRESS>Mohammad Masiur Rahaman
Corresponding Author Assistant Professor, School of Infrastructure
Indian Institute of Technology
Bhubaneswar
Odisha-752050, India
Email<EMAIL_ADDRESS>
###### Abstract

In this article, we propose a thermodynamically consistent phase-field model
for thermo-mechanical fracture and provide an open-source implementation of
the proposed model using a recently developed finite element toolbox, Gridap
in Julia. Here, we have derived the balance equations for the thermo-
mechanical fracture by invoking the virtual power principle and determined the
constitutive relations for the thermodynamic fluxes based on the satisfaction
of the thermodynamic laws. Our proposed formulation provides an equation of
temperature evolution that can easily accommodate dissipative effects such as
viscous damping. One may consider the proposed model as a non-trivial
extension of a recently developed iso-thermal phase-field model by Dhas et al.
[1] for the non-isothermal case. We provide very compact and user-friendly
open-source codes for implementing the proposed model using Gridap in Julia
that requires very low memory usage and gives a high degree of flexibility to
the users in defining weak forms of the governing partial differential
equations. We have validated the proposed model and its implementation against
such standard results available in the literature as crack propagation in the
cruciform shape material, single edge-notched plate, bi-material beam and a
quenching test.

##  1 Introduction

Analysis of thermally induced stresses and cracks are of great importance in
the design of engineering materials. When materials are exposed to a sudden
change in ambient temperature, it can result in a non-uniform change of volume
that induces thermal stresses and leads to the initiation of cracks. For many
practical problems, an iso-thermal model may not be accurate enough to capture
the experimental observations when the effect due to the difference between
ambient and reference temperature is significant. Specifically, for problems
related to thermal shock, materials with inhomogeneous thermal properties,
materials subjected to cooling during manufacturing, etc., consideration of
thermal effect becomes vital. Some of the practical examples are shrinkage
cracks in the concrete, window materials for aircraft [2], additive
manufacturing [3] etc. Thermally-induced cracks may be complex and may follow
paths very different from the one observed for mechanically induced cracks;
typically, cracks propagate from higher heating area to lower heating area
[4]. Thermo-mechanical coupling can significantly modify the behavior and can
lead to unstable crack growth. Also, understanding the mechanism of thermally
induced fracture is essential in the development of novel engineering
materials. Hence, fracture analysis of materials under thermo-mechanical
loading has become one of the contemporary topics of research interest.

Over the years, several experiments have been carried out to understand the
behavior of thermally induced cracks and their effect on the strength or other
essential parameters. Geyer and Nemat-Nasser [5] considered the crack growth
behavior of a glass plate where cooling from the boundaries resulted in
parallel cracks. They have observed stable crack growth for pre-cooling
cracks. The effect of the width of a plate on speed and nature of crack
propagation under the presence of thermal gradient was studied by Ronsin et
al. [6]. Shao et al. [7] studied how the residual strength of ceramics varies
with the stages of crack propagation. To understand the behavior of infrared
window materials used in hypersonic aircraft Qu et al. [2] conducted numerical
(finite element analysis) and experimental studies on chemical vapor deposited
zinc sulfide (CVD ZnS). Yousef et al. [8] experimentally observed the
microcrack evolution in large-grained polycrystalline alumina along with the
finite element simulation. They also reported variation in macro-structural
properties such as young’s modulus due to reduction in temperature. To study
the stable crack growth in glass Grutzik and Reedy [9] developed a novel bi-
material beam specimen that was tested under temperature cooling conditions.
Jiang et al. [10] studied the crack patterns in ceramic due to thermal shock
both experimentally and numerically using the finite element method.

Modeling of fracture involving realistic conditions is tackled with the help
of numerical approaches. Broadly these approaches can be classified as
discrete or discontinuous and diffusive or continuous. Discrete methods
require modification in the mesh as the fracture process advances, thus making
it computationally expensive. The boundary element method (BEM) [11, 12]
addresses the issue related to mesh by representing the domain by boundaries.
Extended finite element method (XFEM) has also been used for modeling the
thermoelastic fracture problems [13, 14] that uses additional enrichment
functions requiring prior knowledge of the fracture behavior. Meshless
approaches [15] have also been used that removes the difficulties involved
with the discretization of the domain. Ai and Augarde [16] used the meshless
method using adaptive cracking particle method that requires some visibility
criteria to model the discontinuities. Recently Greco et al. [17] used the
moving mesh strategy to model the crack propagation in the FEM framework.
Unlike the discrete methods, the diffusive approaches use the effect of
fracture on the governing laws, hence eliminating the difficulties associated
with modification in the meshing. Gradient damage model, peridynamic and
phase-field model, etc., are some of the popular continuum-based approaches.
Localizing gradient damage model (LGDM) [18] was used for thermoelastic
brittle fracture by Sarkar et al. [19]. Recently using the mathematical
similarity between damaged material and Riemannian manifold a new continuum-
based damage model has been proposed by Das et al. [20]. Along with these
models, some non-local continuum damage [21] studies can also be found in the
literature. In this work, we are using the phase-field approach for modeling
fracture under thermo-mechanical loading.

Several phase-field models (PFM) have been independently developed in the
physics and mechanics community. A brief review can be found in [22]. Our
focus is on the latter in which the sharp discontinuity (crack) is
approximated by a smooth transition of a scalar field variable _s_ (known as
phase-field) between _0_ (completely damaged) to _1_ (no damage). In the
phase-field models, the brittle fracture is modeled based on the criteria of
Griffith [23]. The total energy is considered as a sum of elastic and fracture
parts. Based on the similarity of the energy functional with the work of
Mumford and Shah [24] in image segmentation and approximating it as in [25],
Francfort and Marigo [26] posed the problem as energy minimization that helped
researchers to model fracture problems by a set of governing partial
differential equations. Thus, there is no need to redefine the domain,
additional enrichment function, etc. The work of Miehe et al. [27] was
extended for multiphysics problems, including the effect of finite deformation
thermoelasticity by Miehe et al. [28]. Badnava et al. [29] used a mesh
refinement algorithm along with the PFM for thermoelastic brittle fracture. To
simulate the realistic crack propagation under dynamic loading in a
thermoelastic solid Chu et al. [4] proposed a new elastic energy density and
derived the coupling equations based on the energy balance principle. Mandal
et al. [30] extended the work of Wu [31] and Wu and Nguyen [32] for the
thermoelastic case including a monolithic BFGS
(Broyden–Fletcher–Goldfarb–Shanno) algorithm to solve the coupling equations.
Wang et al. [33] provided an Abaqus implementation of the thermoelastic
problem using explicit time integration, as well as source code that can be
used to adapt it to other multifield problems by researchers. The doctoral
dissertation of Liu [34] shows the capability of the thermo-elastic PFM to
simulate fracture in poly-crystalline materials. Hybrid PFM [22] was used to
model fracture in thick plates by Raghu et al. [35] using the higher-order
shear deformation theory [36]. Along with that, Hybrid PFM has also been used
to numerically simulate the thermoelastic brittle fracture by Tangella et al.
[37].

In this work, we have developed a new phase-field model for thermo-mechanical
brittle fracture. In the proposed model, we have derived the balance equations
by invoking the virtual power principle. Using the first and second law of
thermodynamics, we have derived the constitutive relations for the
thermodynamic fluxes, and the temperature evolution equations, in such a way
that any dissipative effects can be easily incorporated into the formulation.
One may consider the proposed model as an extension to the phase-field model
for brittle fracture, presented by Dhas et al. [1] from iso-thermal to a non-
isothermal case. We have provided the details of an open-source implementation
of the proposed model for various numerical examples using a newly developed
finite element toolbox Gridap [38, 39, 40], in Julia [41, 42]. Julia is a
programming language with advantages of simple code formulation as in dynamic
languages like _MATLAB_ and high performance like static languages, e.g. _C_
or _FORTRAN_. Gridap is written in Julia for simulating numerical problems
that are governed by partial differential equations. One of the prime
advantages of Gridap is, it does not use any back-end support to achieve
performance like other FE libraries as in _FEniCS_ [43]. We have organized the
rest of this article in the following way. In Section 2, we have derived the
coupled governing partial differential equations for modeling fracture under
thermo-mechanical loading. In Section 3, we have provided the corresponding
weak form equations and described the steps involved to implement the proposed
PFM using Gridap in Julia. In Section 4, we have shown the capability of the
proposed model by reproducing a few standard numerical and experimental
results available in the literature. We have finally concluded our article in
Section 5.

##  2 Proposed phase-field model

In this section, we provide the details of our proposed model in a small
deformation set-up. The proposed model can be considered as a non-trivial
extension of a recently developed phase-field model by Dhas et al. [1] from an
iso-thermal to a non-isothermal case.

###  2.1 Kinematics

Consider a material body defined by an open set $\Omega$ with boundary
$\partial\Omega$ in the three-dimensional Euclidean space $\mathbb{E}^{3}$.
Corresponding to the displacement field, let $\partial\Omega_{u}$ and
$\partial\Omega_{t}$ be the Dirichlet and Neumann boundaries, respectively.
Corresponding to the temperature field, let $\partial\Omega_{T}$ and
$\partial\Omega_{q}$ be the Dirichlet and Neumann boundaries, respectively. We
define the displacement field as $\bm{u}:\Omega\rightarrow\mathbb{R}^{3}$ at
any instant of time $t\in\mathrm{R}^{+}$. Within a small deformation set-up,
the strain tensor $\bm{\epsilon}$ is given by

$\bm{\epsilon}(\bm{u})=\frac{1}{2}(\nabla\bm{u}+\nabla\bm{u}^{T}),$ (1)

where $(\cdot)^{T}$ denotes the transpose of a tensor and $\nabla$ is the
gradient operator. For thermo-mechanical problems, $\bm{\epsilon}$ can be
considered to be consisting of an elastic strain $\bm{\epsilon}^{\text{elas}}$
and a thermal strain $\bm{\epsilon}^{\text{th}}$, which can be represented as

$\bm{\epsilon}=\bm{\epsilon}^{\text{elas}}+\bm{\epsilon}^{\text{th}},$ (2)

where

$\bm{\epsilon}^{\text{th}}=\alpha(T-T_{0})\bm{I},$ (3)

for a temperature change from an initial temperature $T_{0}$ to current
temperature $T$, considering $\alpha$ as the thermal coefficient of the
material body with $\bm{I}$ denoting the second-order identity tensor.

(a)

(b)

Figure 1: A material body is defined by an open set $\Omega$ with Dirichlet
boundaries ($\partial\Omega_{T}$,$\partial\Omega_{u}$) and Neumann boundaries
($\partial\Omega_{q}$,$\partial\Omega_{t}$) subjected to applied temperature,
displacement, heat flux, and traction $\bar{T}$, $\bar{\bm{u}}$, $\bar{q}$,
$\bar{\bm{t}}$, respectively. Sub-figure (a) is the material body with sharp
crack. For phase-field modeling of fracture a diffused representation of crack
is considered in sub-figure (b) by introducing a length scale $l_{s}$ and a
phase-field variable $s$ with $s=0$ denoting a complete damage, $s=1$ no
damage and $s\in(0,1)$ partial damage.

In the PFM, sharp cracks are approximated by employing a diffused
representation of the cracks via a length scale parameter $l_{s}$ and a phase-
field internal variable $s$ (see Fig. 1). Damaged, undamaged, and partially
damaged states of matter can be described using the values of the phase-field
variable as follows: $s=1$ for undamaged state, $s=0$ for fully damaged state
and $s\in(0,1)$ for a partially damaged state. Considering fully damaged
states as the fracture, crack set may be defined as
$\Omega_{s}=\\{\bm{x}\in\Omega|s(\bm{x})=0\\}$. Here, the damage process is
considered to be irreversible that is if $\bm{x}\in\Omega_{s}$ at time $t_{0}$
then $\bm{x}\in\Omega_{s}$ for all $t\geq t_{0}$. Thus, the deformed and
damaged states of the material body may be described by considering phase-
field variable $s$ as an additional kinematic descriptor along with the
displacement vector $\bm{u}$ and the temperature field $T$.

###  2.2 Virtual power principle and the balance equations

To describe the deformation of a material body under external loading, one can
derive the force balances from a virtual power principle by considering a
macro-force and a micro-force system [44]. While stress tensor $\bm{\sigma}$,
traction vector $\bm{t}(\bm{n})$ and body force $\bm{b}$ define the macro-
force system, the micro-force system includes scalar micro-traction
$\chi(\bm{n})$, vector micro-stress $\bm{\xi}$ and scalar micro-stress $\pi$.
It is important to note that the micro-system is also at the same continuum
level of the macro-system and the coinage ‘micro-force’ may be a misnomer. One
can identify $\dot{\bm{u}}$ as the power conjugate of $\bm{t}(\bm{n})$ and
$\bm{b}$, $\dot{s}$ as the power conjugate of $\chi(\bm{n})$ and express the
external mechanical power for any arbitrary sub-domain
$\mathcal{P}\subset\Omega$ as

$P^{\text{ext}}=\int_{\partial\mathcal{P}}\biggl{[}\bm{t}\left(\bm{n}\right)\cdot\dot{\bm{u}}+\chi(\bm{n})\cdot\dot{s}\biggr{]}\,dA+\int_{\mathcal{P}}\bm{b}\cdot\dot{\bm{u}}\,dV,$
(4)

where, $(\dot{\,\,})$ denotes the time derivative of a variable, $\bm{n}$ the
unit normal vector to the boundary $\partial\mathcal{P}$, $dA$ and $dV$ are,
respectively, measures on $\partial\mathcal{P}$ and $\mathcal{P}$. The
internal mechanical power $P^{\text{int}}$ can be defined by the summation of
power expenditure of $\bm{\sigma}$ over $\dot{\bm{\epsilon}}^{\text{elas}}$,
$\bm{\xi}$ over $\nabla\dot{s}$, $\pi$ over $\dot{s}$ and can be given by

$P^{\text{int}}=\int_{\mathcal{P}}\left(\bm{\sigma}:\dot{\bm{\epsilon}}^{\text{elas}}+\bm{\xi}\cdot\nabla\dot{s}+\pi\,\dot{s}\right)dV.$
(5)

We denote the virtual counterparts of $\dot{\bm{u}}$ and $\dot{s}$ by
$\tilde{\bm{u}}$ and $\tilde{s}$, respectively. Introducing the generalized
virtual velocity vector $\mathcal{V}=\left(\tilde{\bm{u}},\tilde{s}\right)$,
one can express the external virtual mechanical power
$\tilde{P}^{\text{ext}}\left(\mathcal{P};\mathcal{V}\right)$ and the internal
virtual mechanical power
$\tilde{P}^{\text{int}}\left(\mathcal{P};\mathcal{V}\right)$ as

$\tilde{P}^{\text{ext}}\left(\mathcal{P};\mathcal{V}\right)=\int_{\partial\mathcal{P}}\biggl{[}\bm{t}\left(\bm{n}\right)\cdot\tilde{\bm{u}}+\chi(\bm{n})\cdot\tilde{s}\biggr{]}dA+\int_{\mathcal{P}}\bm{b}\cdot\tilde{\bm{u}}\,dV,$
(6)

and

$\tilde{P}^{\text{int}}\left(\mathcal{P};\mathcal{V}\right)=\int_{\mathcal{P}}\left(\bm{\sigma}:\tilde{\bm{\epsilon}}^{\text{elas}}+\bm{\xi}\cdot\nabla\tilde{s}+\pi\tilde{s}\right)dV,$
(7)

respectively. At a given temperature from the symmetry of stress tensor one
can show that
$\bm{\sigma}:\tilde{\bm{\epsilon}}^{\text{elas}}=\bm{\sigma}:\nabla{\tilde{\bm{u}}}$
using of which in the virtual power principle, i.e.,
$\tilde{P}^{\text{int}}\left(\mathcal{P};\mathcal{V}\right)=\tilde{P}^{\text{ext}}\left(\mathcal{P};\mathcal{V}\right)$,
one can arrive at

$\begin{split}\int_{\partial\mathcal{P}}\biggl{[}\bm{t}\left(\bm{n}\right)\cdot\tilde{\bm{u}}+\chi(\bm{n})\cdot\tilde{s}]dA+\int_{\mathcal{P}}(\bm{b}\cdot\tilde{\bm{u}})\,dV=\int_{\mathcal{P}}(\bm{\sigma}:\nabla{\tilde{\bm{u}}}+\\\
\bm{\xi}\cdot\nabla\tilde{s}+\pi\tilde{s})\,dV.\end{split}$ (8)

Appropriate choices of $\mathcal{V}=\left(\tilde{\bm{u}},\tilde{s}\right)$ in
Eq. (8) lead to the macro-force and the micro-force balances as described
below.

####  2.2.1 Macro-force balance

Considering $\mathcal{V}=\left(\tilde{\bm{u}},0\right)$, i.e., by substituting
$\tilde{s}=0$ in Eq. (8), and using the divergence theorem the macro-force
balance Eq. may be obtained as

$\int_{\partial\mathcal{P}}\left(\bm{t}\left(\bm{n}\right)-\bm{\sigma}\bm{n}\right).\tilde{\bm{u}}\,dA=-\int_{\mathcal{P}}\left(\nabla\cdot\bm{\sigma}+\bm{b}\right)\cdot\tilde{\bm{u}}\,dV,$
(9)

which holds for all $\tilde{\bm{u}}$ and any arbitrary sub-domain
$\mathcal{P}$. Applying the localization theorem on Eq. (9), one can get that

$\bm{t}\left(\bm{n}\right)=\bm{\sigma}\,\bm{n},$ (10)

which is the macro-traction condition and

$\nabla\cdot\bm{\sigma}+\bm{b}=\bm{0}.$ (11)

called the macro-force balance. Equations (10) and (11) may be identified as
the classical traction condition and the local linear momentum balance
equation, respectively.

####  2.2.2 Micro-force balance

Considering $\mathcal{V}=\left(\bm{0},\tilde{s}\right)$, i.e., by substituting
$\tilde{\bm{u}}=\bm{0}$ in Eq. (8), the equation for micro-force balance can
be derived as

$\int_{\partial\mathcal{P}}\biggl{[}\chi(\bm{n})-\bm{\xi}\cdot\bm{n}\biggr{]}\cdot\tilde{s}\,dA=-\int_{\mathcal{P}}\left(\nabla\cdot\bm{\xi}-\pi\right)\cdot\tilde{s}\,dV.$
(12)

Since Eq. (12) holds for all $\tilde{s}$ and any sub-domain $\mathcal{P}$, Eq.
(12) may be localized as

$\chi(\bm{n})=\bm{\xi}\cdot\bm{n},$ (13)

and

$\nabla\cdot\bm{\xi}-\pi=0.$ (14)

Equations (13) and (14) are called the micro-traction condition and the micro-
force balance, respectively.

###  2.3 Thermodynamics laws and constitutive modeling

To determine constitutive relations for the thermodynamic fluxes, we have
followed the work by Anand and Gurtin [45]. In presence of heat flux $\bm{q}$,
and heat supply $Q_{\text{h}}$, the first law of thermodynamics, i.e., the law
of conservation of energy, can be given by

$\frac{d}{dt}\int_{P}e\,dV=P^{\text{int}}-\int_{\partial
P}\bm{q}\cdot\bm{n}dA+\int_{P}Q_{h}dV,$ (15)

where, $e$ is the internal energy density and the power balance i.e.
$P^{ext}=P^{int}$ is used. From the second law of thermodynamics, we have

$\frac{d}{dt}\int_{P}\eta\,dV\geq-\int_{\partial
P}\frac{\bm{q}\cdot\bm{n}}{T}\,dA+\int_{P}\frac{Q_{\text{h}}}{T}\,dV,$ (16)

where $\eta$ is the entropy per unit volume. The inequality given by Eq. (16)
ensures that the rate of entropy increase is not less than the entropy input
in the system. Using the divergence theorem on Eq. (15) and Eq. (16), and
applying the localization theorem we get

$\dot{e}=-\nabla\cdot\bm{q}+Q_{\text{h}}+\bm{\sigma}:\dot{\bm{\epsilon}}^{\text{elas}}+\bm{\xi}\cdot\nabla\dot{s}+\pi\,\dot{s},$
(17)

and

$\dot{\eta}\geq-\frac{1}{T}\nabla\cdot\bm{q}+\frac{1}{T^{2}}\bm{q}\cdot\nabla
T+\frac{Q_{\text{h}}}{T}.$ (18)

From the relation between Helmholtz free energy and internal energy density
given by $\psi=e-T\eta$, we have

$\dot{\psi}=\dot{e}-\dot{T}\eta-T\dot{\eta}.$ (19)

Employing equations (11), (14), (17), (18) and (19), we get that

$\dot{\psi}-\bm{\sigma}:\dot{\bm{\epsilon}}^{\text{elas}}-\bm{\xi}\cdot\nabla\dot{s}-\pi\,\dot{s}+\eta\,\dot{T}+\frac{1}{T}\bm{q}\cdot\nabla
T\leq 0.$ (20)

The inequality given by Eq. (20) needs to be satisfied for a thermodynamic
process to be valid, and one can determine constitutive functions for
$\bm{\sigma}$, $\pi$, $\eta$ and $\bm{\xi}$ on satisfying Eq. (20).

###  2.4 Constitutive response functions

We assume the Helmholtz free energy $\psi$ is a function of
$\bm{\epsilon}^{\text{elas}},\,s,\,\nabla s$, and $T$, i.e.,

$\psi=\hat{\psi}(\bm{\epsilon}^{\text{elas}},s,\nabla s,T).$ (21)

Using the chain rule, the rate of Helmholtz energy can be expressed as

$\dot{\psi}=\partial_{\bm{\epsilon}^{\text{elas}}}\psi:\dot{\bm{\epsilon}}^{\text{elas}}+\partial_{s}\psi\,\dot{s}+\partial_{\nabla
s}\psi\cdot{\nabla\dot{s}}+\partial_{T}\psi\,\dot{T}.$ (22)

Considering $\pi$ as sum of an energetic part $\pi^{\text{en}}$ and a
dissipative part $\pi^{\text{dis}}$, i.e.
$\pi=\pi^{\text{en}}+\pi^{\text{dis}}$, and employing Eq. (22), we can rewrite
Eq. (20) as

$\begin{split}(\bm{\sigma}-\partial_{\bm{\epsilon}^{\text{elas}}}\psi):\dot{\bm{\epsilon}}^{\text{elas}}+(\pi^{\text{en}}-\partial_{s}\psi)\dot{s}+(\bm{\xi}-\partial_{\nabla
s}\psi)\cdot{\nabla\dot{s}}+\\\
\pi^{\text{dis}}\dot{s}-(\eta+\partial_{T}\psi)\dot{T}-\frac{1}{T}\bm{q}\cdot\nabla
T\geq 0.\end{split}$ (23)

Applying Collman and Noll [46] procedure to Eq. (23). we obtain the following
constitutive relations

$\bm{\sigma}=\partial_{\bm{\epsilon}^{\text{elas}}}\psi,$ (24)
$\bm{\xi}=\partial_{\nabla s}\psi,$ (25) $\pi^{\text{en}}=\partial_{s}\psi,$
(26) $\eta=-\partial_{T}\psi,$ (27)

and the inequality given by (23) reduces to

$\pi^{\text{dis}}\dot{s}-\frac{1}{T}\bm{q}\cdot\nabla T\geq 0.$ (28)

Determination of the constitutive relation must be done in such a way that the
inequality constraint given by Eq. (28) always satisfy. From the
irreversiblity condition on damage i.e. $\dot{s}\leq 0$, it can be realised
that a possible choice for $\pi^{\text{dis}}$ could be
$\pi^{\text{dis}}=\mathcal{G}\dot{s}$, where $\mathcal{G}$ is a constitutive
function with $\mathcal{G}\leq 0$. Now, for the satisfaction of the
inequality, the second term also has to be positive. Hence, we define
$\bm{q}=-k\,\nabla T$, where $k$ is the thermal conductivity.

###  2.5 Specialized constitutive relations

We assume the Helmholtz free energy of the system as

$\psi(\bm{\epsilon}^{\text{elas}},s,\nabla
s,T)=\psi_{\text{elas}}(\bm{\epsilon}^{\text{elas}})+\psi_{\text{frac}}(s,\nabla
s)+\psi_{\text{ther}}(T),$ (29)

where $\psi_{\text{elas}}$ is the elastic energy, $\psi_{\text{frac}}$ the
fracture energy and $\psi_{\text{ther}}$ the thermal energy. For imposing the
condition that cracks do not propagate under pure compression, we consider a
decomposition of strain tensor into a deviatoric part
$\bm{\epsilon}_{\text{dev}}$ and a volumetric part
$\bm{\epsilon}_{\text{vol}}$, i.e.,

$\bm{\epsilon}=\bm{\epsilon}_{\text{dev}}+\bm{\epsilon}_{\text{vol}},$ (30)

where

$\bm{\epsilon}_{\text{dev}}=\mathbb{P}_{\text{dev}}\bm{\epsilon};\,\,\,\bm{\epsilon}_{\text{vol}}=\mathbb{P}_{\text{vol}}\bm{\epsilon}.$
(31)

In Eq. (31), $\mathbb{P}_{\text{dev}}$ and $\mathbb{P}_{\text{vol}}$ projects
the strain tensor onto its deviatoric and volumetric parts, respectively. For
the assumption of no damage under pure compression, the elastic energy can be
expressed as

$\psi_{\text{elas}}(\bm{\epsilon}^{\text{elas}},s)=g(s)\,\psi_{\text{elas}}^{+}+\psi_{\text{elas}}^{-}.$
(32)

In Eq. (32), two parts of the elastic energy are given by

$\psi_{\text{elas}}^{+}=\frac{1}{2}g(s)\left(\langle
p\rangle_{+}\bm{I}:\bm{\epsilon}^{\text{elas}}_{\text{vol}}+\mathbb{P}_{\text{dev}}\mathbb{C}\bm{\epsilon}^{\text{elas}}:\bm{\epsilon}^{\text{elas}}_{\text{dev}}\right),$
(33)

and

$\psi_{\text{elas}}^{-}=\frac{1}{2}\langle
p\rangle_{-}\bm{I}:\bm{\epsilon}^{\text{elas}}_{\text{dev}},$ (34)

where $\mathbb{C}$ is the undamaged fourth order elasticity tensor,
$\bm{\epsilon}^{\text{elas}}_{\text{vol}}=\mathbb{P}_{\text{vol}}\,\bm{\epsilon}^{\text{elas}}$,
$\bm{\epsilon}^{\text{elas}}_{\text{dev}}=\mathbb{P}_{\text{dev}}\,\bm{\epsilon}^{\text{elas}}$
and g(s) is a degradation function considered as $s^{2}$ in this work,
$\langle p\rangle_{+}=\frac{1}{2}\left(p+\lvert p\rvert\right)$ and $\langle
p\rangle_{-}=\frac{1}{2}\left(p-\lvert p\rvert\right)$,
$p=\frac{1}{3}\bm{I}:\mathbb{P}_{\text{vol}}\mathbb{C}\bm{\epsilon}^{\text{elas}}$.
The fracture energy $\psi_{\text{frac}}(\nabla s,s)$ may be postulated as

$\psi_{\text{frac}}\left(\nabla
s,s\right)=G_{c}\biggl{(}\frac{\left(1-s\right)^{2}}{2l_{s}}+\frac{l_{s}}{2}\nabla
s\cdot\nabla s\biggr{)},$ (35)

where $G_{c}$ is the critical energy release rate and $l_{s}$ is the phase-
field length scale parameter. Recognizing the fact that the heat flux may get
affected by crack generation we express the thermal conductivity $k$ as

$k=g(s)\,k_{0},$ (36)

where $k_{0}$ is the undamaged thermal conductivity and the expression given
by Eq. (36) can account for reduced thermal conduction due to fracture
considering the adiabatic fracture zone. With the specialized total energy,
one can find out the expressions for thermodynamic fluxes by using the
expressions from equations (24), (25) and (26) as

$\bm{\sigma}=g(s)\left(\langle
p\rangle_{+}\bm{I}+\mathbb{P}_{\text{dev}}\mathbb{C}\bm{\epsilon}^{\text{elas}}\right)+\langle
p\rangle_{-}\bm{I},$ (37) $\bm{\xi}=G_{c}l_{s}\nabla s,$ (38)
$\begin{split}\pi^{en}=g^{\prime}(s)\left(\langle
p\rangle_{+}\bm{I}:\bm{\epsilon}^{\text{elas}}_{\text{vol}}+\mathbb{P}_{\text{dev}}\mathbb{C}\bm{\epsilon}^{\text{elas}}:\bm{\epsilon}^{\text{elas}}_{\text{dev}}\right)-\\\
\frac{G_{c}}{l_{s}}(1-s).\end{split}$ (39)

Using the expression for the above thermodynamic fluxes in Eq. (11) and Eq.
(14), one can arrive at the governing equation for displacement field and
phase-field, respectively.

###  2.6 Temperature evolution

To derive the temperature evolution equation, we are following the approach
similar to the one adopted by Rahaman et al. [47] in a small deformation set-
up with the presence of fracture energy instead of energy contribution due to
plasticity. Employing Eq. (27) and Eq. (29) in the relation $e=\psi+T\eta$, we
get that

$e=\psi_{\text{elas}}+\psi_{\text{frac}}+\psi_{\text{ther}}-T\frac{d\psi_{\text{ther}}}{dT}.$
(40)

For the case of no mechanical strain, differentiation of Eq. (40) leads to

$de=-T\frac{d^{2}\psi_{\text{ther}}}{dT^{2}}dT=\rho\,c\,dT,$ (41)

where $c$ is the specific heat. Employing Eq. (40) and Eq. (41), we get that

$\dot{e}=\dot{\psi}_{\text{elas}}+\dot{\psi}_{\text{frac}}+\rho c\dot{T}.$
(42)

Expressing $\dot{\psi}_{\text{elas}}$ and $\dot{\psi}_{\text{frac}}$ in their
unspecialized form as

$\dot{\psi}_{\text{elas}}=\bm{\sigma}:\dot{\bm{\epsilon}}^{\text{elas}},$ (43)
$\dot{\psi}_{\text{frac}}=\bm{\xi}\cdot\nabla\dot{s}+\pi^{\text{en}}\dot{s},$
(44)

and using Eq. (17) in Eq. (42) we get the temperature evolution equation as
given by

$\rho c\dot{T}=-\nabla\cdot\bm{q}+Q_{\text{h}}+\pi^{\text{dis}}\dot{s}.$ (45)

As can be seen from Eq. (45), our proposed formulation provides a
thermodynamically consistent accommodation of dissipative effects such as
temperature change due to viscous damping through the
$\pi^{\text{dis}}\dot{s}$ in Eq. (45).

###  2.7 Boundary value problem and the strong form

Using the expressions derived in the previous section, one may write the
strong form of the governing PDEs as a boundary value problem and derive the
corresponding weak form of the governing PDEs for the finite element
formulation of the proposed phase-field model. One can express the stress
tensor $\bm{\sigma}$ in terms of kinematic variables as
$\mathbb{C}_{\text{mod}}\,\bm{\epsilon}^{\text{elas}}$ by defining
$\mathbb{C}_{\text{mod}}$ as

$\mathbb{C}_{\text{mod}}=\begin{cases}g(s)\mathbb{C},\,\,\,\,\,\text{for}\,\,p\geq
0&\\\
g(s)\left(\mathbb{P}_{\text{dev}}\mathbb{C}\right)+\mathbb{P}_{\text{vol}}\mathbb{C},\,\,\,\,\text{for}\,\,p<0.\end{cases}$
(46)

Using $\bm{\sigma}=\mathbb{C}_{\text{mod}}\,\bm{\epsilon}^{\text{elas}}$, Eq.
(11) can be re-written as

$\nabla\cdot\biggl{(}\mathbb{C}_{\text{mod}}\,\bm{\epsilon}^{\text{elas}}\biggr{)}+\bm{b}=\bm{0}.$
(47)

Similarly, using Eq. (14), (38) and (39), we can re-write the force balance
corresponding to phase field as

$\begin{split}\nabla\cdot(G_{c}l_{s}\nabla s)-g^{\prime}(s)\left(\langle
p\rangle_{+}\bm{I}:\bm{\epsilon}^{\text{elas}}_{\text{vol}}+\mathbb{P}_{\text{dev}}\mathbb{C}\bm{\epsilon}^{\text{elas}}:\bm{\epsilon}^{\text{elas}}_{\text{dev}}\right)\\\
+\frac{G_{c}}{l_{s}}\left(1-s\right)+\mathcal{G}\dot{s}=0.\end{split}$ (48)

As it can be seen from the above equation, any dissipative effect can easily
be incorporated in to the equation through the term $\mathcal{G}\dot{s}$,
bypassing the empirical route by Miehe et al. [27]. However, in the numerical
examples of the present work, we have not considered any energy dissipation
and implementation of this model for problems with dissipative effect can be
considered as a future scope of the current work. To account for the
irreversiblity condition on damage, a history function
$\mathcal{H}(\mathcal{E})$, where
$\mathcal{H}(f)=\text{max}_{\tau\in[0,t]}f(\tau)$ for any input argument $f$,
of the so-called positive part of elastic energy
$\mathcal{E}=\psi^{+}_{\text{elas}}(\bm{\epsilon}^{\text{elas}})$ is employed
[48]. Using the history function $\mathcal{H}(\mathcal{E})$, Eq. (48) may be
expressed as

$\nabla\cdot\biggl{(}G_{c}l_{s}\nabla
s\biggr{)}-g^{\prime}(s)\mathcal{H}(\mathcal{E})+\frac{G_{c}}{l_{s}}\left(1-s\right)=0,$
(49)

where

$\mathcal{H}(\mathcal{E})=\langle
p\rangle_{+}\bm{I}:\bm{\epsilon}^{\text{elas}}_{\text{vol}}+\mathbb{P}_{\text{dev}}\mathbb{C}\bm{\epsilon}^{\text{elas}}:\bm{\epsilon}^{\text{elas}}_{\text{dev}}.$
(50)

For $\pi^{\text{dis}}=0$ and $\bm{q}=-k\,\nabla T$, the temperature evolution
equation given by Eq. (45) can be re-written as

$\rho c\dot{T}-\nabla\cdot(k\,\nabla T)=Q_{\text{h}}.$ (51)

The governing PDEs (47), (49) and (51) are subjected to boundary conditions,
such as prescribed displacement $\bm{u}=\bar{\bm{u}}$, prescribed temperature
$T=\bar{T}$ and applied traction $\bm{t}(\bm{n})=\bar{\bm{t}}(\bm{n})$ on
$\partial\Omega_{u}$, $\partial\Omega_{T}$, $\partial\Omega_{t}$,
respectively. Equations (47), (49), and (51) together with the boundary
conditions are called the strong form of the governing PDEs. For instance, in
absence of body force i.e. $\bm{b}=\bf{0}$, one possible strong form can be
given as

$\displaystyle\nabla\cdot\bm{\sigma}={\bf{0}}\,\,\,\,\,\left(\text{where}\,\,\bm{\sigma}=\mathbb{C}_{\text{mod}}\,\bm{\epsilon}^{\text{elas}}\right)\,\,\,\,\text{in}\,\,\,\,\Omega,$
(52a)
$\displaystyle\bm{u}=\bar{\bm{u}}\,\,\,\,\text{on}\,\,\,\,\partial\Omega_{u},$
(52b)
$\displaystyle\bm{\sigma}\,\bm{n}=\bm{0}\,\,\,\,\text{on}\,\,\,\,\partial\Omega\backslash\partial\Omega_{u},$
(52c)

corresponding to the macro-force balance,

$\displaystyle\nabla\cdot\biggl{(}G_{c}l_{s}\nabla
s\biggr{)}-2s\,\mathcal{H}(\mathcal{E})+\frac{G_{c}}{l_{s}}\left(1-s\right)=0\,\,\,\,\text{in}\,\,\,\,\Omega,$
(53a) $\displaystyle\nabla
s\cdot\bm{n}=0\,\,\,\,\text{on}\,\,\,\,\partial\Omega\backslash\partial\Omega_{u},$
(53b)

corresponding to the micro-force balance, and

$\displaystyle\rho\,c\,\dot{T}-\nabla\cdot(k\,\nabla
T)=Q_{\text{h}}\,\,\,\,\text{in}\,\,\,\,\Omega,$ (54a)
$\displaystyle{T}=\bar{{T}}\,\,\,\,\text{on}\,\,\,\,\partial\Omega_{T},$ (54b)
$\displaystyle\bm{q}\cdot\bm{n}=0\,\,\,\,\text{on}\,\,\,\,\partial\Omega\backslash\partial\Omega_{T},$
(54c)

corresponding to the temperature evolution. Considering isotropic materials
for the present study, we express the components of the fourth order
elasticity tensor as
$\mathbb{C}_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})$,
where $\lambda$ and $\mu$ are the Lamé parameters and $\delta$ denotes the
Kronecker delta. The Lamé parameters are related to Young’s modulus $E$ and
Poisson’s ratio $\nu$ by $\lambda=E\nu/\left((1+\nu)(1-2\nu)\right)$ and
$\mu=E/\left(2(1+\nu)\right)$.

##  3 The weak form and an outline for finite element implementation in Julia

In this section, we have derived the weak form for the strong form given by
equations (52), (53) and (54) and outlined the finite element formulation
through a numerical example. Let $\bm{v}$, $\phi$ and $w$ be the test
functions corresponding to displacement $\bm{u}$, phase-field $s$, and
temperature field $T$, respectively. The trial spaces with the given
displacement and temperature boundary condition may be given by

$H_{u}=\\{\bm{u}\in
H^{1}(\Omega)^{d};\;\bm{u}=\bar{\bm{u}}\;\text{on}\;\partial\Omega_{u}\\},$
(55) $H_{s}=\\{s\in H^{1}(\Omega)\\},$ (56) $H_{T}=\\{T\in
H^{1}(\Omega);\;T=\bar{T}\;\text{on}\;\partial\Omega_{T}\\},$ (57)

where $\bar{\bm{u}}$ and $\bar{T}$ are prescribed displacement and temperature
on $\partial\Omega_{u}$ and $\partial\Omega_{T}$, respectively, and $d$
corresponds to the dimension of displacement vector. The test spaces may be
defined as

$V_{u}=\\{\bm{v}\in
H^{1}(\Omega)^{d};\;\bm{v}=0\;\text{on}\;\partial\Omega_{u}\\},$ (58)
$V_{s}=\\{\phi\in H^{1}(\Omega)\\},$ (59) $V_{T}=\\{w\in
H^{1}(\Omega);\;w=0\;\text{on}\;\partial\Omega_{T}\\}.$ (60)

One can obtain the following weak form: Find $\bm{u}\in H_{u}$, $s\in H_{s}$
and $T\in H_{T}$ such that for all $\bm{v}\in V_{u}$, $\phi\in V_{s}$ and
$w\in V_{T}$

$\displaystyle
a_{\text{Disp}}((\bm{u},T),(\bm{v},w))=b_{\text{Disp}}(\bm{v},w),$ (61a)
$\displaystyle a_{\text{PF}}(s,\phi)=b_{\text{PF}}(\phi)$ (61b)

where

$\displaystyle\begin{split}a_{\text{Disp}}((\bm{u},T),(\bm{v},w))=\int_{\Omega}(\bm{\epsilon}(\bm{v}):\bm{\sigma}_{T}(\bm{u},T)+w\,\rho\,c\,\dot{T}\\\
+\nabla w\cdot(k\,\nabla T))\,d{\Omega},\end{split}$ (62a) $\displaystyle
b(\bm{v},w)=\int_{\Omega}\left(\bm{\epsilon}(\bm{v}):\bm{\sigma}_{F}+Q_{\text{h}}\right)\,d{\Omega},$
(62b) $\displaystyle
a_{\text{PF}}(s,\phi)=\int_{\Omega}\bigl{(}G_{c}\,l_{s}\,\nabla
s\cdot\nabla\phi+2s\,\phi\,\mathcal{H}+\frac{G_{c}}{l_{s}}\,s\,\phi\bigr{)}\,d\Omega,$
(62c) $\displaystyle b_{\text{PF}}(\phi)=\int_{\Omega}\frac{G_{c}}{l_{s}}\phi
d{\Omega},$ (62d)

with the definitions of $\bm{\sigma}_{T}(\bm{u},T)$ and $\bm{\sigma}_{F}$ as

$\displaystyle\bm{\sigma}_{T}(\bm{u},T)=g(s)\left(\langle
p\rangle_{+}\bm{I}+\mathbb{P}_{\text{dev}}\mathbb{C}(\bm{\epsilon}-\alpha\bm{I}T)\right)+\langle
p\rangle_{-}\bm{I},$ (63a) $\displaystyle\bm{\sigma}_{F}=-g(s)\left(\langle
p\rangle_{+}\bm{I}+\mathbb{P}_{\text{dev}}\mathbb{C}\alpha\bm{I}T_{0}\right)+\langle
p\rangle_{-}\bm{I}.$ (63b)

For numerical implementation of the proposed model, we have adopted the
staggered approach proposed by Miehe et al. [48]. At each step, we have solved
for the temperature and the displacement field from the couple set of
equations (51) and (52) using the MultiFieldFESpace available in gridap. We
have computed the solution for the phase field variable using Eq. (53) based
on the history function given by Eq. (50).

One of the most appealing aspects of Gridap is that it allows one to use the
weak form directly for finite element implementation. We describe all the
steps associated with finite element simulations for the proposed model, such
as creating the mesh file, implementing the weak form, applying boundary
conditions, solving for unknown field variables, and the post-processing
through a Julia code for numerical simulation of brittle fracture in a single
edge notched plate under thermo-mechanical loading (see Section 4.2). One can
easily modify the Julia code for simulating other thermo-mechanical fracture
problems. To replicate the results reported in Section (4), as a first step
one needs to import the Julia packages listed below into a script file
currently created in a jupyter notebook.

⬇

using GridapGmsh

using Gridap

using Gridap.Geometry

using Gridap.TensorValues

using PyPlot

One can define parameters related to applied boundary conditions for a single
edge notched plate under thermo-mechanical loading by the following lines in
Julia.

⬇

const T0 = 300

const TAppMax = T0 + 50

const delt = 1e-2

const tMax = 1

const uMax = 1.2e-6

AppVel = uMax/tMax

uMin = 0

uTran = 0.2*uMax

We have imposed the prescribed temperature on Dirchlet boundary
$\partial\Omega_{T}$ in small increments till the temperature reaches at the
maximum applied temperature and kept at constant maximum applied temperature
till the failure of the specimen. For non-linear smooth increment of
temperature, one may use the following function.

⬇

using SymPy

x,$x_{1}$ = symbols("x,$x_{1}$", real = true)

heaviside(x) = 0.5 * (sign(x) + 1)

interval(x, a, b) = heaviside(x-a) - heaviside(x-b)

hS = uMax/10

F(x) = (T0 - TAppMax) * interval(x,-4*hS+uMin,uTran)

$w_{h}$(x,$x_{1}$) = (1/(sqrt(2*pi)*hS))*exp(-(x-$x_{1}$)^2/(2*hS^2))

smoothT = SymPy.integrate(F($x_{1}$)*$w_{h}$(x,$x_{1}$),($x_{1}$,-4*hS,uMax))
+ TAppMax

plot(smoothT,0,uMax)

One may also increase temperature linearly from initial temperature to the
maximum applied temperature based on the convergence requirements using the
below function.

⬇

function Tfun(u)

if u <= uTran

return ((TAppMax - T0)/uTran)*u + T0

else

return TAppMax

end

end

plot(Tfun,0,uMax)

In this work, we have used a non-linear smooth increment of temperature by
choosing option one of the two functions as given below.

⬇

uAppVec = range(0,uMax,length = Int64(floor(tMax/delt)))

AppTOption = 1 ## 1 for smooth and otherwise linear than constant

if AppTOption == 1

TAppVec = smoothT.(uAppVec)

else

TAppVec = Tfun.(uAppVec)

end

Input material parameters related to elastic, fracture and thermal properties
of the notched plate can be defined by the following lines in Julia.

⬇

const E = 340.0e9

const $\nu$ = 0.22

const Gc = 42.47

const lsp = 3.33e-6

const $\eta$ = 1e-8

const $\alpha$ = 8e-6

const c = 0.775

const $\kappa$_mat = 300.0

const $\rho$ = 2450

For the finite element simulations, one needs to have the discretization of
the computational domain that can be generated by writing a mesh file
“SquarePlateWithEdgeNotch.msh” in Julia (see Julia code for meshing of single
edge notched plate) and loading that mesh file to build an instance of
“DiscreteModel” by the following lines.

⬇

model = GmshDiscreteModel("SquarePlateWithEdgeNotch.msh")

writevtk(model,"SquarePlateWithEdgeNotch")

We have defined a function for the constitutive matrix of isotropic material
under plane stress or plane strain condition by the following lines in Julia.

⬇

function ElasFourthOrderConstTensor($E$,$\nu$,PlanarState)

# 1 for Plane Stress and 2 Plane Strain Condition

if PlanarState == 1

C1111 = $E$/(1-$\nu$*$\nu$)

C1122 = ($\nu$*$E$)/(1-$\nu$*$\nu$)

C1112 = 0.0

C2222 = $E$/(1-$\nu$*$\nu$)

C2212 = 0.0

C1212 = $E$/(2*(1+$\nu$))

elseif PlanarState == 2

C1111 = ($E$*(1-$\nu$*$\nu$))/((1+$\nu$)*(1-$\nu$-2*$\nu$*$\nu$))

C1122 = ($\nu$*$E$)/(1-$\nu$-2*$\nu$*$\nu$)

C1112 = 0.0

C2222 = ($E$*(1-$\nu$))/(1-$\nu$-2*$\nu$*$\nu$)

C2212 = 0.0

C1212 = $E$/(2*(1+$\nu$))

end

C_ten =
SymFourthOrderTensorValue(C1111,C1112,C1122,C1112,C1212,C2212,C1122,C2212,C2222)

return C_ten

end

For a plane strain condition, the constitutive tensor can be computed by
calling the above Julia function as given below.

⬇

const C_mat = ElasFourthOrderConstTensor(E,$\nu$,2)

To satisfy the assumption that crack can not propagate under pure compression,
stress and strain tensors are decomposed into a volumetric and a deviatoric
part by introducing the projection operators $\mathbb{P}_{\text{vol}}$ and
$\mathbb{P}_{\text{dev}}$, respectively, which are defined by the following
lines.

⬇

I2 = SymTensorValue{2,Float64}(1.0,0.0,1.0)

I4 = I2 $\otimes$ I2

I4_sym = one(SymFourthOrderTensorValue{2,Float64})

P_vol = (1.0/3)*I4

P_dev = I4_sym - P_vol

To express the stress tensor in terms of kinematics variables, i.e.,
$\bm{\sigma}=\mathbb{C}_{\text{mod}}\,\bm{\epsilon}$, where the expression of
$\mathbb{C}_{\text{mod}}$ is given by Eq. (46), the following function is
defined in Julia.

⬇

$\sigma$_elas($\varepsilon$Elas) = C_mat $\odot$ $\varepsilon$Elas

function $\sigma$_elasMod($\varepsilon$, $\varepsilon$_in, s_in,T,T_in)

$\varepsilon$Elas_in = $\varepsilon$_in - $\alpha$*(T_in-T0)*I2

$\varepsilon$Elas = $\varepsilon$ - $\alpha$*(T-T0)*I2

if tr($\varepsilon$Elas_in) >= 0

$\sigma$ = (s_in^2 + $\eta$)*$\sigma$_elas($\varepsilon$Elas)

elseif tr($\varepsilon$Elas_in) < 0

$\sigma$ = (s_in^2 + $\eta$) *P_dev $\odot$ $\sigma$_elas($\varepsilon$Elas) +
P_vol$\odot$ $\sigma$_elas($\varepsilon$Elas)

end

return $\sigma$

end

The function for $\bm{\sigma}_{T}$ given by Eq. (63) can be defined as given
below.

⬇

function $\sigma$_totMod($\varepsilon$, $\varepsilon$_in,s_in,T,T_in)

$\varepsilon$Elas_in = $\varepsilon$_in -$\alpha$*(T_in-T0)*I2

$\varepsilon$ElasTot = $\varepsilon$ -$\alpha$*T*I2

if tr($\varepsilon$Elas_in) >= 0

$\sigma$T = (s_in^2 + $\eta$)*$\sigma$_elas($\varepsilon$ElasTot)

elseif tr($\varepsilon$Elas_in) < 0

$\sigma$T = (s_in^2 + $\eta$) *P_dev $\odot$
$\sigma$_elas($\varepsilon$ElasTot) + P_vol$\odot$
$\sigma$_elas($\varepsilon$ElasTot)

end

return $\sigma$T

end

Similarly, the function for $\bm{\sigma}_{F}$ given by Eq. (63) can be defined
by the following lines.

⬇

function $\sigma$_thermMod($\varepsilon$_in,s_in,T_in)

$\varepsilon$Elas_in = $\varepsilon$_in - $\alpha$*(T_in-T0)*I2

$\varepsilon$ElasTher = $\alpha$*(T0)*I2

if tr($\varepsilon$Elas_in) >= 0

$\sigma$F = (s_in^2 + $\eta$)*$\sigma$_elas($\varepsilon$ElasTher)

elseif tr($\varepsilon$Elas_in) < 0

$\sigma$F = (s_in^2 + $\eta$)*P_dev $\odot$
$\sigma$_elas($\varepsilon$ElasTher) + P_vol$\odot$
$\sigma$_elas($\varepsilon$ElasTher)

end

return $\sigma$F

end

One can determine the positive part of the elastic free energy
$\psi_{\text{elas}}^{+}(\bm{\epsilon}^{\text{elas}})$ given by Eq. (50) by
writing the following function in Julia.

⬇

function $\Psi$Pos($\varepsilon$_in,T_in)

$\varepsilon$Elas_in = $\varepsilon$_in - $\alpha$*(T_in-T0)*I2

if tr($\varepsilon$Elas_in) >= 0

$\Psi$Plus = 0.5*(($\varepsilon$Elas_in) $\odot$
$\sigma$_elas($\varepsilon$Elas_in))

elseif tr($\varepsilon$Elas_in) < 0

$\Psi$Plus = 0.5*((P_dev $\odot$ $\sigma$_elas($\varepsilon$Elas_in)) $\odot$
(P_dev $\odot$ ($\varepsilon$Elas_in)))

end

return $\Psi$Plus

end

One needs to generate a discrete approximation of the finite element test and
trial spaces of the problem on the discretized computational domain.
Approximation of the finite element spaces associated with the phase field
variable can be done by the following lines in Julia.

⬇

order = 1

reffe_PF = ReferenceFE(lagrangian,Float64,order)

V0_PF = TestFESpace(model,reffe_PF;

conformity=:H1)

U_PF = TrialFESpace(V0_PF)

sh = zero(V0_PF)

One can generate the approximation of the finite element spaces associated
with the displacement field by writing the following lines in Julia.

⬇

reffe_Disp = ReferenceFE(lagrangian,VectorValue{2,Float64},order)

V0_Disp = TestFESpace(model,reffe_Disp;

conformity=:H1,

dirichlet_tags=["BottomEdge","TopEdge"],

dirichlet_masks=[(true,true),(false,true)])

uh = zero(V0_Disp)

Similarly, finite element spaces for the temperature field can be generated as
given below.

⬇

reffe_Temp = ReferenceFE(lagrangian,Float64,order)

V0_Temp = TestFESpace(model,reffe_Temp;

conformity=:H1,

dirichlet_tags=["BottomEdge","TopEdge"])

We have defined a combined finite element space for displacement and
temperature field to develop a multi-field solution algorithm using the
following code.

⬇

V0 = MultiFieldFESpace([V0_Disp,V0_Temp])

To compute the integrals in the weak form given by Eq. (62) numerically, one
needs to define an integration mesh along with an integration rule (for
example, Gauss quadrature) in each of the cells in the triangulation. Using
Gridap, one can easily define the integration mesh and the corresponding
Lebesgue measure by using the built-in functions “Triangulation” and
“Measure”, respectively. For instance, one can use the following lines for
computing the integrals in the weak form given by Eq. (62) defined in the
domain $\Omega$ using a quadrature rule of degree two times the order of
interpolation in the cells of the triangulation.

⬇

degree = 2*order

$\Omega$ = Triangulation(model)

d$\Omega$ = Measure($\Omega$,degree)

One can find normal to a specific load boundary and define a measure on that
boundary to compute a boundary integral using the following built-in functions
available in Gridap.

⬇

labels = get_face_labeling(model)

LoadTagId = get_tag_from_name(labels,"TopEdge")

$\Gamma$_Load = BoundaryTriangulation(model,tags = LoadTagId)

d$\Gamma$_Load = Measure($\Gamma$_Load,degree)

n_$\Gamma$_Load = get_normal_vector($\Gamma$_Load)

To find the values of variables that are defined by using the built-in
function “CellState” at the Gauss points, one need to use the “project”
function as defined below.

⬇

function project($q$,model,d$\Omega$,order)

reffe = ReferenceFE(lagrangian,Float64,order)

V = FESpace(model,reffe,conformity=:L2)

a($u$,$v$) = $\int$($u$*$v$)*d$\Omega$

b($v$) = $\int$($v$*$q$)*d$\Omega$

op = AffineFEOperator(a,b,V,V)

qh = Gridap.solve(op)

return qh

end

In the present study, a staggered scheme is used to update the solution from
the pseudo time $t_{n}$ to $t_{n+1}$ [48]. Given the displacement and
temperature field and the history function at the time $t_{n}$, one can update
the phase-field at the time $t_{n+1}$ by using the following function.

⬇

function stepPhaseField(uh_in,$\Psi$PlusPrev_in)

a_PF(s,$\Phi$) = $\int$( Gc*lsp*$\nabla$($\Phi$)$\cdot$ $\nabla$(s) +
2*$\Psi$PlusPrev_in*s*$\Phi$ + (Gc/lsp)*s*$\Phi$ )*d$\Omega$

b_PF($\Phi$) = $\int$( (Gc/lsp)*$\Phi$ )*d$\Omega$

op_PF = AffineFEOperator(a_PF,b_PF,U_PF,V0_PF)

sh_out = Gridap.solve(op_PF)

return sh_out

end

Following function has been used for defining the product between the
conductivity tensor and gradient of temperature in the weak form corresponding
to temperature field (see Eq. (61)).

⬇

$\kappa$GradTemp($\nabla$,s_in) = (s_in^2 + $\eta$)*$\kappa$_mat*$\nabla$

Approximating the time derivative of temperature by using backward difference
method as

$\dot{T}=\frac{T_{n+1}-T_{n}}{t_{n+1}-t_{n}}=\frac{T_{n+1}-T_{n}}{\Delta t},$
(64)

one can update the displacement vector and temperature field at time $t_{n+1}$
by solving coupled set of equations using the values of temperature vector and
the history function at time $t_{n}$, and the computed value of phase-field at
time $t_{n+1}$. Using a built-in function MultiFieldFESpace in Gridap, we have
written a function for updating the displacement and temperature as given
below.

⬇

function stepDispTemp(uh_in,sh_in,T_in,vApp,TApp,delt)

uApp1(x) = VectorValue(0.0,0.0)

uApp2(x) = VectorValue(0.0,vApp)

U_Disp = TrialFESpace(V0_Disp,[uApp1,uApp2])

Tapp1(x) = T0

Tapp2(x) = TApp

Tg = TrialFESpace(V0_Temp,[Tapp1, Tapp2])

U = MultiFieldFESpace([U_Disp,Tg])

a((u,T),(v,w)) = $\int$( ($\varepsilon$(v) $\odot$
($\sigma$_totMod$\circ$($\varepsilon$(u),$\varepsilon$(uh_in),sh_in,T,T_in)))
+ $\nabla$(w)$\cdot$($\kappa$GradTemp$\circ$ ($\nabla$(T),sh_in)) +
(($\rho$*c*T*w)/delt))*d$\Omega$

b((v,w)) = $\int$((($\rho$*c*T_in*w)/delt) - ($\varepsilon$(v) $\odot$
($\sigma$_thermMod$\circ$($\varepsilon$(uh_in),sh_in,T_in))))*d$\Omega$

op = AffineFEOperator(a,b,U,V0)

uhTh = Gridap.solve(op)

uh_out,Th_out = uhTh

return uh_out,Th_out

end

Once both displacement and phase-field are determined at time $t_{n+1}$, one
can update the energy history function at time $t_{n+1}$ by defining it as

$\mathcal{H}_{n+1}=\begin{cases}\mathcal{E}_{n+1}\,\,\text{for}\,\,\mathcal{E}_{n+1}>\mathcal{E}_{n}&\\\
\mathcal{E}_{n}\,\,\text{otherwise},\end{cases}$ (65)

where the history function $\mathcal{H}(\mathcal{E})$ of elastic free energy
$\mathcal{E}=\psi_{\text{elas}}^{+}(\bm{\epsilon}^{\text{elas}})$ need to be
defined in Julia (see Section 2.7) that can be achieved by writing the
following lines.

⬇

function new_EnergyState($\psi$PlusPrev_in,$\psi$hPos_in)

$\psi$Plus_in = $\psi$hPos_in

if $\psi$Plus_in >= $\psi$PlusPrev_in

$\psi$Plus_out = $\psi$Plus_in

else

$\psi$Plus_out = $\psi$PlusPrev_in

end

true,$\psi$Plus_out

end

Finally, one can simulate the thermo-mechanical brittle fracture in the
notched plate by applying a monotonic displacement control loading and solve
for the unknown field variables at each loading step by writing the main
routine in Julia that uses the above-defined functions as listed below.

⬇

t = 0.0

innerMax = 10

count = 0

tol = 1e-8

Load = Float64[]

Displacement = Float64[]

push!(Load, 0.0)

push!(Displacement, 0.0)

$\Psi$PlusPrev = CellState(0.0,d$\Omega$)

sPrev = CellState(1.0,d$\Omega$)

sh = project(sPrev,model,d$\Omega$,order)

ThPrev = CellState(T0,d$\Omega$)

Th = project(ThPrev,model,d$\Omega$,order)

while t .< tMax

count = count .+ 1

t = t + delt

vApp = AppVel*t

TApp = TAppVec[count]

for inner = 1:innerMax

$\Psi$hPlusPrev = project($\Psi$PlusPrev,model,d$\Omega$,order)

RelErr = abs(sum($\int$( Gc*lsp*$\nabla$(sh)$\cdot$ $\nabla$(sh) +
2*$\Psi$hPlusPrev*sh*sh + (Gc/lsp)*sh*sh)*d$\Omega$ - $\int$(
(Gc/lsp)*sh)*d$\Omega$))/abs(sum($\int$( (Gc/lsp)*sh)*d$\Omega$))

sh = stepPhaseField(uh,$\Psi$hPlusPrev)

uh,Th = stepDispTemp(uh,sh,Th,vApp,TApp,delt)

$\Psi$hPos_in = $\Psi$Pos$\circ$($\varepsilon$(uh),Th)

update_state!(new_EnergyState,$\Psi$PlusPrev,$\Psi$hPos_in)

if RelErr < tol

break

end

end

Node_Force = sum($\int$( n_$\Gamma$_Load $\cdot$
($\sigma$_elasMod$\circ$($\varepsilon$(uh),$\varepsilon$(uh),sh,Th,Th)) )
*d$\Gamma$_Load)

push!(Load, Node_Force[2])

push!(Displacement, vApp)

end

One can create output files at each loading step to view the results in
ParaView. For instance, one can write a “.vtu” file to save data for the
solution of the displacement field, phase-field, and temperature field at each
loading step by using the following lines in Julia.

⬇

writevtk($\Omega$,"results_PhaseFieldThermoElastic",cellfields=

["uh"=>uh,"s"=>sh,"T"=>Th])

One may generate the load-displacement curve by using the plot command as
given below.

⬇

plot(Displacement,Load)

Figure 2: Geometry and boundary conditions for a cruciform shape notched plate
($\bar{u}_{1}$ and $\bar{u}_{2}$ corresponds to the specified displacement
along the direction $x_{1}$ and $x_{2}$ , respectively, $\bar{T}$ corresponds
to the specified temperature, $\partial\Omega_{1}$, $\partial\Omega_{2}$,
$\partial\Omega_{3}$, $\partial\Omega_{4}$ denotes left, right, top, and
bottom boundaries, respectively).

##  4 Numerical simulations

In this section, we have considered a few standard problems available in the
literature for validating our proposed phase-field model and its Julia
implementation. The validation problems include cruciform shaped material
subjected to different boundary conditions, single edge notched specimen under
tension, bi-material beam subjected to reduction of temperature, and a ceramic
beam quenched to lower temperature. We have used constant strain triangle
(CST) elements for finite element solutions of the problems. The simulation
results show the capability of the proposed model to reproduce the results
available in the literature.

###  4.1 Cruciform shape plate with notch

We have considered a plate of cruciform shape with a notch as shown in Fig. 2.
The problem has been introduced by Prasad et al. [49], and since then, it has
been considered as a benchmark problem for numerical analysis of thermo-
mechanical fracture. We have taken $L$ as $50$ mm and chosen material
properties from [49] and [30], which are $\emph{E}=218400$ $N/m^{2}$,
$\nu=0.2$, $\emph{G}_{c}=0.0002\,N/m$, $\rho=0.0$, $\alpha=6\times
10^{-4}\,\tccentigrade$, $k_{0}=1.0\,W/m\tccentigrade$,
$c=1.0J/kg\,\tccentigrade$ . For numerical simulations, we have assumed plane
stress condition similar to Mandal et al. [30], and we have discretized the
domain using CST elements. The Julia code for creating the mesh file for a
particular case is given in Julia code for meshing of cruciform shape notched
plate.

We have considered the following three loading conditions, including the
boundary conditions at the edges for all three cases shown in Fig. 2.

  1. a.

Application of quasi-static load in the form of incremental displacement of
$\bar{u}_{2}=3.5\times 10^{-7}$m at each time step on $\partial\Omega_{3}$,
with $\bar{T}=0$ on $\partial\Omega_{3}$ and $\partial\Omega_{4}$ for the
consideration of pure mechanical loading case.

  2. b.

Increasing the temperature on $\partial\Omega_{3}$ with simultaneous reduction
on $\partial\Omega_{4}$ without any restriction on $\bar{u}_{2}$ on
$\partial\Omega_{3}$ to consider a pure thermal loading case. At each time
step, temperature increases by $0.1\tccentigrade$ on $\partial\Omega_{3}$
while reducing the temperature on $\partial\Omega_{4}$ by $-0.1\tccentigrade$
and keeping the temperature on boundaries $\partial\Omega_{1}$ and
$\partial\Omega_{2}$ at the reference temperature.

  3. c.

The combined effect of mechanical and thermal loading cases. With increment in
displacement as $\bar{u}_{2}=1.92\times 10^{-7}$m on $\partial\Omega_{3}$ and
temperature increment of $0.06\tccentigrade$ and $-0.06\tccentigrade$ on
$\partial\Omega_{3}$ and $\partial\Omega_{4}$ , respectively at each time
step.

(a)

(b)

(c)

Figure 3: Simulated fracture results for cruciform specimen for (a) pure
mechanical loading ($3.4\times 10^{-5}$m) (b) pure thermal loading at
($8\tccentigrade$) and (c) combined mechanical and thermal loading
($1.53\times 10^{-5}m\,\,\text{and}\,\,4.8\tccentigrade$)

The length of the notch was taken to be $0.2L$ with an inclination of
$135^{\circ}$ counterclockwise from horizontal. The length scale parameter
$l_{s}$ was taken to be $5\times 10^{-4}$m. Furthermore, a refined mesh was
used in the expected portion of the domain where the crack may propagate with
the edge length of $l_{s}/4$. The obtained crack pattern for all the three
cases is demonstrated in Fig. 3. As can be seen from Fig. 3, the crack
propagation for pure mechanical loading case is completely horizontal, where
for the thermal loading case the crack initiates vertically from the notch and
eventually starts bending right. For the combined case the crack propagates in
between the purely mechanical and pure thermal case i.e. the crack propagates
diagonally. The obtained crack propagation directions shows similar trend to
that obtained by Duflot [13], Prasad et al. [49], and Mandal et al. [30].

###  4.2 Single edge notched plate under tension

We have considered a notched plate as shown in Fig. 4 and chosen the material
properties as $\emph{E}=32\times 10^{9}$ $N/m^{2}$, $\nu=0.2$,
$\emph{G}_{c}=3\,N/m$, $\rho=2450kg/m^{3}$, $\alpha=8\times
10^{-6}\,\tccentigrade$, $k_{0}=300\,W/m\tccentigrade$,
$c=0.775J/kg\,\tccentigrade$ , taken from [37]. We have used CST elements for
the finite element discretization of the notched plate and the Julia code for
creating the mesh file is given in Julia code for meshing of single edge
notched plate.

Figure 4: Geometry and boundary conditions for a single edge notched plate
under tension ($\bar{T}_{t}$ and $\bar{T}_{b}$ are the specified temperature
at the top and bottom boundary, respectively, $\bar{\bm{u}}_{2}$ is the
specified displacement along the direction $x_{2}$, dimensions are in mm).

(a)

(b)

(c)

Figure 5: Simulated fracture results for single edge notched plate under
tension at an applied displacement $\bar{\bm{u}}_{2}=8.4\times 10^{-4}mm.$ and
applied temperature increment (a) $\Delta T=-75\tccentigrade$ (b) $\Delta
T=0\tccentigrade$ and (c) $\Delta T=+75\tccentigrade$.

We have reproduced the results for the following cases assuming plane strain
condition [37]:

  1. a.

Application of quasi-static load in the form of incremental displacement at
the top boundary, with restrained displacement boundary conditions (see Fig.
4) while maintaining the temperature of the domain at $300\,\tccentigrade$ for
the consideration of pure mechanical loading case.

  2. b.

Incrementing $\bar{T}_{t}$ upto $(T_{0}+50)\tccentigrade$ and keeping
$\bar{T}_{b}$ at $T_{0}$ for the consideration of thermo-mechanical coupling
along with incremental displacement at top. A detailed description of the
implementation of this case has been discussed in section 3.

  3. c.

Decrementing $\bar{T}_{t}$ upto $(T_{0}-50)\tccentigrade$ and keeping
$\bar{T}_{b}$ at $T_{0}$ for the consideration of thermo-mechnical coupling
under opposite gradient along with incremental displacement at top.

  4. d.

Incrementing $\bar{T}_{t}$ up to $(T_{0}+75)\tccentigrade$ and keeping
$\bar{T}_{b}$ at $T_{0}$ for the consideration of thermo-mechanical coupling
with higher temperature gradient than the previous (see b) along with
incremental displacement at top.

  5. e.

Decrementing $\bar{T}_{t}$ upto $(T_{0}-75)\tccentigrade$ and keeping
$\bar{T}_{b}$ at $T_{0}$ for the consideration of thermo-mechnical coupling
under opposite higher temperature gradient than previous (see c) gradient
along with incremental displacement at top.

Figure 6: Load-displacement curve for single edge notched plate under tension
for $\Delta T$ = $-75\tccentigrade$, $-50\tccentigrade$, $0\tccentigrade$,
$+50\tccentigrade$, $+75\tccentigrade$.

The crack propagation for load case (a), (d) and (e) has been shown in Fig. 5.
The delay in crack propagation was observed with the increasing temperature at
the top boundary, which can also be seen in the load-displacement plots with
the delay in peak (see Fig. 6). Case (b) shows a similar behavior as case (d)
by showing delay in crack propagation in comparison to pure mechanical case
but apparent early crack propagation in comparison to case (d) characterized
by the load peak in-between case (a) and case (d). And since in case (c) the
temperature is in between $0\tccentigrade$ and $-75\tccentigrade$ the position
of load peak shows that the obtained crack position will also be in between
the one obtained for case (a) and case (e). We have reproduced the results
reported in Tangella et al. [37].

###  4.3 Notched bi-material beam

We have chosen a bi-material beam as reported in [9] to simulate the behavior
of fracture due to the dissimilar thermal coefficient of the materials and the
advantage of using gridap in case of multi-material problems by defining the
multi-material constitutive laws based on the assigned material tags. The
geometry of the considered model has been shown in Fig. 7(a). The crack
position, height and material properties is chosen from Mandal et al. [30],
which are $\emph{E}=64\times 10^{3}$ MPa, $\nu=0.2$, $\emph{G}_{c}=0.4\,N/mm$,
$\alpha=3.25\times 10^{-6}\,\tccentigrade$ corresponding to borosilicate
glass, and $\emph{E}=193\times 10^{3}$ MPa, $\nu=0.29$, $\alpha=193\times
10^{-6}\,\tccentigrade$ corresponding to $304$ stainless steel. We have
considered the length of the specimen as $150$ mm. A notch of width $1$ mm and
height $15$ mm has been assumed at a distance of $30$ mm from the left edge.
The length scale parameter $l_{s}$ has been chosen to be $0.3$ mm; the mesh is
refined in the area where the crack propagation was expected with the elements
of size $l_{s}/4$. Plane stress condition has been assumed for the simulation.
Julia code for the finite element meshing of the specimen is given in Julia
code for meshing of bi-material beam. The problem consists of cooling the
material from a reference temperature. The reference temperature has been
taken as 300 $\tccentigrade$. The temperature of the domain is slowly reduced
at each time step. Since we are reducing the temperature uniformly at each
time step over the whole domain, we need not consider the temperature
evolution equation. The resulting displacement can be directly obtained by
using the governing equation corresponding to the displacement field. This is
similar to the one adopted by Mandal et al. [30] wherein the load has been
applied in the form of temperature strain based on the temperature change. The
crack propagates because of thermal stresses due to the different coefficients
of thermal expansion. Fig. 7(b) shows the crack path obtained in the domain.
Also, Fig. 8 shows the experimental crack pattern reported in [9] and the
crack path obtained by this study near the crack tip. One can observe that the
predicted crack path by the proposed model is in good agreement with the
experimental results. Prediction of crack path using the proposed model is
similar to the prediction by Mandal et al. [30] that also does not show
experimentally observed little drop in the crack path before going straight .

(a)

(b)

Figure 7: (a) Geometry of bi-material beam (dimensions are in mm) (b) crack
propagation in the bi-material beam due to gradual reduction of temperature
inducing due to the difference in temperature coefficient of the materials.

(a)

(b)

(c)

Figure 8: Crack path in a bi-material beam: (a) experimentally observed [9]
(b) prediction by the PF-CZM model [30] and (c) prediction using the proposed
model.

###  4.4 Quenching test

The experimental results generated by Shao et al. [7] for quenching test is
considered as a benchmark problem for numerical modeling of thermo-mechanical
problems. For the present case, we have considered a specimen with an initial
temperature of $680$ K and $880$ K quenched into an ambient temperature of
$300$ K. For the model development, we have considered only a quarter of the
actual specimen due to the symmetry of the problem, as shown in Fig. 9 with a
height of $4.9$ mm and length of $50$ mm. The horizontal and vertical
displacements have been fixed on the right and top edges, respectively. The
bottom edge and right edge have been subjected to quenching through conduction
by applying a fixed $300$ K temperature. The length scale parameter has been
selected as (length of the specimen)/250, and the element size has been
selected such that $l_{s}$ is at least two times the mesh size. The Julia code
used for creating the finite element mesh for quenching test specimen is given
in Julia code for meshing of quenching test specimen. For the simulation, we
have considered the same material properties as mentioned in Section 4.2 under
plane stress condition similar to Mandal et al. [30]. We haven’t considered
the degradation of the $k$ considering the complex flow of water through the
developed cracks for this problem. The distribution of temperature field
across the domain is given in Fig. 10 and phase-field variable across the
domain is demonstrated in Fig. 11 along with the experimentally found results
by [7]. The results show that, as the initial temperature increases, i.e., for
higher temperature difference, the crack density is higher, and the depth of
crack propagation is also higher. The obtained results are found to be in good
agreement with the experimental results of Shao et al. [7], and numerical
results by Wang et al. [33], and Tangella et al. [37].

(a)

(b)

Figure 9: Geometry and boundary condition for quenching test: (a) complete
domain (b) quarter of the domain with symmetric boundary conditions analyzed
in the present simulation (dimensions are in mm).

(a)

(b)

Figure 10: Variation of temperature field across the domain for (a) model with
initial temperature 680 K and (b) model with initial temperature 880 K.

(a)

(b)

(c)

(d)

Figure 11: Crack patterns in quenching tests: (a) experimentally observed
crack pattern for a test with initial temperature 680 K, (b) predicted crack
pattern by the proposed phase-field for a test with initial temperature 680 K,
(c) experimentally observed crack pattern for a test with initial temperature
880 K, and (d) predicted crack pattern by the proposed phase-field model for a
test with initial temperature 880 K.

##  5 Conclusions

In this work, we have developed a thermodynamically consistent phase-field
model for brittle fracture under thermo-mechanical loading. We have provided
an open-source implementation of the developed model using the finite element
toolbox, Gridap in Julia. Our proposed model provides an equation for
temperature evolution that can account for the effects of dissipative forces
such as viscous damping whenever needed. We have validated the proposed model
and its implementation against a few representative examples available in the
literature. We have outlined all the steps to implement our proposed model in
Julia for simulating fracture in a single edge notched plate that one can
easily modify to generate results for other thermo-mechanical fracture
problems. Numerical results for a cruciform-shaped plate show the capability
of our proposed model to reproduce the behavior under steady-state conditions.
We have reproduced the effect of thermo-mechanical coupling on the global
response (load-displacement curve) for a single edge notched plate under
tension. The proposed model shows excellent agreement with the experimentally
observed crack paths for a bi-material beam and the quenching test.

Our proposed implementation is available with an open-source license. Hence,
it eliminates the technical barrier for practitioners and researchers
interested in exploring the phase-field model for solving a wide range of
thermo-mechanical brittle fracture problems. Moreover, our proposed
implementation will expose the users to many built-in packages of Julia that
may be useful for researchers who want to extend the proposed implementation
for the case of ductile fracture or other applications. The proposed
multifield solution scheme can be implemented very easily for other coupled
problems with the help of Gridap. Most importantly, the availability of an
efficient and compact open-source code will provide better flexibility to the
researchers and will establish a high standard for the research work.

## Data accessibility

The present work does not generate any experimental data. Julia scripts used
as the source codes to generate the numerical results presented in the
manuscript can be downloaded as jupyter notebook files from Julia code for
phase-field modeling of thermo-mechanical fracture.

## Declaration of competing interest

Authors of this article declare that they have no conflict of interest.

## Acknowledgments

Authors of this article gratefully acknowledge support from the Indian
Institute of Technology Bhubaneswar under the grant SP107.

## References

  * [1] Dhas, B., Rahaman, M., Akella, K., Roy, D., Reddy, J., et al., 2018.  “A phase-field damage model for orthotropic materials and delamination in composites”.  Journal of Applied Mechanics, 85(1).
  * [2] Qu, Z., Cheng, X., He, R., Pei, Y., Zhang, R., and Fang, D., 2016.  “Rapid heating thermal shock behavior study of cvd zns infrared window material: Numerical and experimental study”.  Journal of Alloys and Compounds, 682, pp. 565–570.
  * [3] Gouge, M., and Michaleris, P., 2017.  Thermo-mechanical modeling of additive manufacturing.  Butterworth-Heinemann.
  * [4] Chu, D., Li, X., and Liu, Z., 2017.  “Study the dynamic crack path in brittle material under thermal shock loading by phase field modeling”.  International Journal of Fracture, 208(1), pp. 115–130.
  * [5] Geyer, J. F., and Nemat-Nasser, S., 1982.  “Experimental investigation of thermally induced interacting cracks in brittle solids”.  International Journal of Solids and Structures, 18(4), pp. 349–356.
  * [6] Ronsin, O., Heslot, F., and Perrin, B., 1995.  “Experimental study of quasistatic brittle crack propagation”.  Physical review letters, 75(12), p. 2352.
  * [7] Shao, Y., Zhang, Y., Xu, X., Zhou, Z., Li, W., and Liu, B., 2011.  “Effect of crack pattern on the residual strength of ceramics after quenching”.  Journal of the American Ceramic Society, 94(9), pp. 2804–2807.
  * [8] Yousef, S. G., Rödel, J., Fuller Jr, E., Zimmermann, A., and El-Dasher, B., 2005\.  “Microcrack evolution in alumina ceramics: experiment and simulation”.  Journal of the American Ceramic Society, 88(10), pp. 2809–2816.
  * [9] Grutzik, S. J., and Reedy, E. D., 2018.  “Crack path selection in thermally loaded borosilicate/steel bibeam specimen”.  Experimental Mechanics, 58(1), pp. 1–10.
  * [10] Jiang, C., Wu, X., Li, J., Song, F., Shao, Y., Xu, X., and Yan, P., 2012.  “A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock”.  Acta Materialia, 60(11), pp. 4540–4550.
  * [11] Prasad, N., Aliabadi, M., and Rooke, D., 1994.  “The dual boundary element method for thermoelastic crack problems”.  International Journal of Fracture, 66(3), pp. 255–272.
  * [12] Prasad, N., Aliabadi, M., and Rooke, D., 1996.  “The dual boundary element method for transient thermoelastic crack problems”.  International journal of solids and structures, 33(19), pp. 2695–2718.
  * [13] Duflot, M., 2008.  “The extended finite element method in thermoelastic fracture mechanics”.  International Journal for Numerical Methods in Engineering, 74(5), pp. 827–847.
  * [14] Zamani, A., and Eslami, M. R., 2010.  “Implementation of the extended finite element method for dynamic thermoelastic fracture initiation”.  International Journal of Solids and Structures, 47(10), pp. 1392–1404.
  * [15] Pant, M., Singh, I., and Mishra, B., 2010.  “Numerical simulation of thermo-elastic fracture problems using element free galerkin method”.  International Journal of Mechanical Sciences, 52(12), pp. 1745–1755.
  * [16] 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, pp. 343–357.
  * [17] Greco, F., Ammendolea, D., Lonetti, P., and Pascuzzo, A., 2021.  “Crack propagation under thermo-mechanical loadings based on moving mesh strategy”.  Theoretical and Applied Fracture Mechanics, p. 103033.
  * [18] Poh, L. H., and Sun, G., 2017.  “Localizing gradient damage model with decreasing interactions”.  International Journal for Numerical Methods in Engineering, 110(6), pp. 503–522.
  * [19] Sarkar, S., Singh, I., and Mishra, B., 2020.  “A thermo-mechanical gradient enhanced damage method for fracture”.  Computational Mechanics, 66(6), pp. 1399–1426.
  * [20] Das, S., Sharma, S., Ramaswamy, A., Roy, D., and Reddy, J., 2021.  “A geometrically inspired model for brittle damage in compressible elastomers”.  Journal of Applied Mechanics, 88(8), p. 081002.
  * [21] Bazˇant, Z. P., and Pijaudier-Cabot, G., 1988.  “Nonlocal Continuum Damage, Localization Instability and Convergence”.  Journal of Applied Mechanics, 55(2), 06, pp. 287–293.
  * [22] Ambati, M., Gerasimov, T., and De Lorenzis, L., 2015.  “A review on phase-field models of brittle fracture and a new fast hybrid formulation”.  Computational Mechanics, 55(2), pp. 383–405.
  * [23] Griffith, A. A., 1921.  “Vi. the phenomena of rupture and flow in solids”.  Philosophical transactions of the royal society of london. Series A, containing papers of a mathematical or physical character, 221(582-593), pp. 163–198.
  * [24] Mumford, D. B., and Shah, J., 1989.  “Optimal approximations by piecewise smooth functions and associated variational problems”.  Communications on pure and applied mathematics.
  * [25] Ambrosio, L., and Tortorelli, V. M., 1990.  “Approximation of functional depending on jumps by elliptic functional via t-convergence”.  Communications on Pure and Applied Mathematics, 43(8), pp. 999–1036.
  * [26] 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), pp. 1319–1342.
  * [27] Miehe, C., Welschinger, F., and Hofacker, M., 2010.  “Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations”.  International journal for numerical methods in engineering, 83(10), pp. 1273–1311.
  * [28] 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, pp. 449–485.
  * [29] Badnava, H., Msekh, M. A., Etemadi, E., and Rabczuk, T., 2018.  “An h-adaptive thermo-mechanical phase field model for fracture”.  Finite Elements in Analysis and Design, 138, pp. 31–47.
  * [30] Mandal, T. K., Nguyen, V. P., Wu, J.-Y., Nguyen-Thanh, C., and de Vaucorbeil, A., 2021.  “Fracture of thermo-elastic solids: Phase-field modeling and new results with an efficient monolithic solver”.  Computer Methods in Applied Mechanics and Engineering, 376, p. 113648.
  * [31] Wu, J.-Y., 2017.  “A unified phase-field theory for the mechanics of damage and quasi-brittle failure”.  Journal of the Mechanics and Physics of Solids, 103, pp. 72–99.
  * [32] Wu, J.-Y., and Nguyen, V. P., 2018.  “A length scale insensitive phase-field damage model for brittle fracture”.  Journal of the Mechanics and Physics of Solids, 119, pp. 20–42.
  * [33] Wang, T., Ye, X., Liu, Z., Liu, X., Chu, D., and Zhuang, Z., 2020.  “A phase-field model of thermo-elastic coupled brittle fracture with explicit time integration”.  Computational Mechanics, 65(5), pp. 1305–1321.
  * [34] Liu, Z., 2021.  “Phase-field modelling of crack initiation and propagation in polycrystalline materials under thermomechanical loadings”.
  * [35] Raghu, P., Rajagopal, A., and Reddy, J., 2020.  “Thermodynamically consistent variational approach for modeling brittle fracture in thick plates by a hybrid phase field model”.  Journal of Applied Mechanics, 87(2), p. 021002.
  * [36] Reddy, J. N., 1984.  “A Simple Higher-Order Theory for Laminated Composite Plates”.  Journal of Applied Mechanics, 51(4), 12, pp. 745–752.
  * [37] Tangella, R. G., Kumbhar, P., and Annabattula, R. K., 2021.  “Hybrid phase-field modeling of thermo-elastic crack propagation”.  International Journal for Computational Methods in Engineering Science and Mechanics, pp. 1–16.
  * [38] Verdugo, F., and Badia, S., 2019.  “A user-guide to gridap–grid-based approximation of partial differential equations in julia”.  arXiv preprint arXiv:1910.01412.
  * [39] Badia, S., and Verdugo, F., 2020.  “Gridap: An extensible finite element toolbox in julia”.  Journal of Open Source Software, 5(52), p. 2520.
  * [40] Rahaman, M. M., 2021.  “An open-source implementation of a phase-field model for brittle fracture using gridap in julia”.  arXiv preprint arXiv:2109.05927.
  * [41] Bezanson, J., Karpinski, S., Shah, V. B., and Edelman, A., 2012.  “Julia: A fast dynamic language for technical computing”.  arXiv preprint arXiv:1209.5145.
  * [42] Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B., 2017.  “Julia: A fresh approach to numerical computing”.  SIAM Review, 59(1), pp. 65–98.
  * [43] Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M., and Wells, G., 2015.  “Archive of numerical software: The fenics project version 1.5”.  University Library Heidelberg.
  * [44] Gurtin, M. E., 1996.  “Generalized ginzburg-landau and cahn-hilliard equations based on a microforce balance”.  Physica D: Nonlinear Phenomena, 92(3-4), pp. 178–192.
  * [45] Anand, L., and Gurtin, M. E., 2003.  “Thermal effects in the superelasticity of crystalline shape-memory materials”.  Journal of the Mechanics and Physics of Solids, 51(6), pp. 1015–1058.
  * [46] Collman, B., and Noll, W., 1963.  “The thermodynamics of elastic materials with heat conduction”.  Arch. Ration. Mech. Anal, 13, p. 167.
  * [47] Rahaman, M. M., Dhas, B., Roy, D., and Reddy, J., 2017.  “A dynamic flow rule for viscoplasticity in polycrystalline solids under high strain rates”.  International Journal of Non-Linear Mechanics, 95, pp. 10–18.
  * [48] Miehe, C., Hofacker, M., and Welschinger, F., 2010.  “A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits”.  Computer Methods in Applied Mechanics and Engineering, 199(45-48), pp. 2765–2778.
  * [49] Prasad, N., Aliabadi, M., and Rooke, D., 1994.  “Incremental crack growth in thermoelastic problems”.  International journal of fracture, 66(3), pp. R45–R50.