# An eXtended Finite Element Method Implementation in COMSOL Multiphysics:
Thermo-Hydro-Mechanical Modeling of Fluid Flow in Discontinuous Porous Media

Ahmad Jafari  Mohammad Vahab<EMAIL_ADDRESS>Pooyan Broumand  Nasser
Khalili  School of Civil and Environmental Engineering, University of New
South Wales, Sydney 2052, Australia  Department of Civil and Environmental
Engineering, Shiraz University, Shiraz, Iran

###### Abstract

This paper presents the implementation of the eXtended Finite Element Method
(XFEM) in the general-purpose commercial software package COMSOL Multiphysics
for multi-field thermo-hydro-mechanical problems in discontinuous porous
media. To this end, an exclusive enrichment strategy is proposed in compliance
with the COMSOL modeling structure. COMSOL modules and physics interfaces are
adopted to take account of the relevant physical processes involved in thermo-
hydro-mechanical coupling analysis, namely: the mechanical deformation, fluid
flow in porous media and heat transfer. Essential changes are made to the
internal variables of the physics interfaces to ensure consistency in the
evaluation of enriched solution fields. The model preprocessing, level-set
updates, coupling of the relevant physics and postprocessing procedures are
performed adopting a coherent utilization of the COMSOL’s built-in features
along with the COMSOL’s LiveLink for MATLAB functions. The implementation
process, remedies for the treatment of the enriched zones, XFEM framework
setup, multiphysics coupling, numerical integration and numerical solution
strategy are described in detail. The capabilities and performance of the
proposed approach are investigated by examining several multi-field thermo-
hydro-mechanical simulations involving single/multiple discontinuities in
2D/3D porous rock settings.

###### keywords:

XFEM; COMSOL Multiphysics; Thermo-hydro-mechanical coupling; Discontinuous
porous media.

††journal: Elsevier

##  1 Introduction

The thermo-hydro-mechanical (THM) coupled processes due to fluid flow in
deformable porous media, subject to natural discontinuities, has been a
crucial area of interest for the prediction of the physical response in many
engineering problems in geotechnics, mining, petroleum engineering, water
resources, reservoir engineering and environmental engineering, to name a few
[1, 2, 3, 4, 5, 6]. Over the past decades, concurrent experimental [7, 8],
analytical [9] and numerical [10, 11, 12, 13, 14] efforts were undertaken to
develop predictive tools for the analysis of the role of discontinuities on
the THM processes in deformable porous media [15, 2]. With excellent inherent
flexibility in tackling all types of discontinuities, computational
superiority, and other algorithmic advantages (e.g., circumventing the need
for remeshing, data transfer, and mesh refinement for high gradients),
extended finite element method (XFEM) [16, 17, 18, 19] has emerged as one of
the most versatile tools for the study of discontinuities in deformable porous
rock media. In this respect, XFEM has been extensively adopted for the
extension of hydro-mechanical [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
32] and thermo-hydro-mechanical [33, 34, 35, 36, 37, 38] frameworks to the
analysis of fluid flow within partially/fully saturated deformable porous
media in the presence of discontinuities.

Nevertheless, the application of XFEM to real-life engineering problems has
been relatively limited, particularly with respect to the analysis of
multiphysics processes in porous media. This has primarily been due to the
fact that the majority of contributions in this context have been due to in-
house codes [39, 40], with limited availability and capacity for handling
complex geometries, large-scale domain and 3D analysis. To alleviate this
shortcoming, implementations of XFEM in general-purpose FE packages, such as
ABAQUS, has been proposed through development of an elaborate host of user-
defined subroutines with applications in linear elastic fracture mechanics
(LEFM) by Giner et al. [41], cohesive fractures by Xu and Yuan [42] and Haddad
and Sepehrnoori [25], and contact mechanics by Ooi et al. [43]. Such
implementations of XFEM, despite featuring efficient built-in solvers and
advanced meshing tools in 2D/3D settings, have been limited to a limited range
of physics, hindering the generalized applicability of XFEM to the study of
multi-field problems.

COMSOL Multiphysics is considered as one of the pioneering software packages
which enables efficient handling of intricate coupling systems with enhanced
efficiency and reliability [44, 45, 46]. It facilitates both fully-coupled
(monolithic) and segregated (staggered) solution strategies [47], which can be
user-customized depending on the interplaying physics, accessibility of
computational resources, and numerical stability concerns. To date, COMSOL has
been employed as a versatile platform for the phase-field implementation of a
coupled hydro-mechanical analysis of porous media (e.g., see Zhou et al. [48,
49]). However, no such extension of this powerful platform has been made to
include XFEM within a multi-physics platform for multi-field practical
solutions.

In this work, underpinned by our recent developments on the implementation of
XFEM for solid mechanics [50], and by taking advantage of the exceptional
flexibility of this software in dealing with any arbitrary coupling processes,
we present a thermo-hydro-mechanical framework for the XFEM analysis of
discontinuities in deformable porous media. In this context, An exclusive XFEM
framework for the study of THM coupling processes, compatible with the
structure of COMSOL is introduced. The THM coupling process is carried out by
exploiting the generic “Solid Mechanics”, “Darcy’s Law” and “Heat Transfer in
Porous Media” physics interfaces in a conjugate manner to take account of the
standard and enriched terms of the solution field. The discontinuity
interfaces are tracked by means of level set functions introduced via external
MATLAB functions, which accommodate the evolution of cracks prior to (i.e., at
the pre-processing stage) and during the analysis. Moreover, by using internal
functions and variables of the software, stress intensity factor calculation
as well as numerical contact analysis are performed within the Multiphysics
theme. The proposed framework enables robust and efficient tackling of
coupling systems with enhanced reliability in 2D/3D settings. The extended
framework, in turn, facilitates the inclusion of further related physics as
required.

The paper is organized as follows: In section 2, the XFEM formulation of
thermo-hydro-mechanical coupling in discontinuous deformable porous media is
briefly described in conjunction with the weak forms and discretization. In
section 3, the algorithms employed for the implementation of XFEM in COMSOL
are presented, which involves a detailed discussion on the essential steps
related to pre-processing of the interfaces and set-up of modules. Section 4
is devoted to the numerical simulations which explore the performance of the
extended framework in dealing with a range of benchmark problems in the
literature. The concluding remarks are presented in section 5. The models
developed are available from https://github.com/ahmadjafari93/xfem-comsol.git.

##  2 XFEM formulation

In the context of the XFEM, the weak/strong interfaces such as cracks, voids,
and material heterogeneities in the domain are introduced independently from
the background mesh by taking advantage of the partition of unity concept
[51]. Such description is achieved by enriching the standard finite element
discretization with special enrichment functions corresponding to the desired
solution fields [52]. In the case of thermo-hydro-mechanical modeling of fluid
flow in porous media, the fully coupled model is constructed by means of three
independent yet interacting sub-models [33]; a poro-elastic model for the
solid skeleton, a Darcy model for the fluid flow through the pore devoid, and
a heat transfer model for the pore-solid mixture. In the following sections,
the strong forms of the governing equations of each of the processes involved
are presented in detail. Next, the corresponding weak forms and finite element
discretization are obtained in compliance with the structure of the software
[53, 50].

###  2.1 Governing equations

The thermo-hydro-mechanical model is based on Biot’s theory of the pore fluid-
solid mixture [54] in conjunction with the heat transfer equations. The linear
momentum balance equation for the discontinuous porous domain $\Omega$ which
is bounded by $\Gamma$ (Fig. 1) can be expressed as

$\nabla\cdot{\boldsymbol{\sigma}}+\rho{\bf{b}}=0$ (1)

Figure 1: Schematic representation of the porous problem domain, boundary
conditions and discontinuity zone for the corresponding thermo-hydro-
mechanical analysis.

where ${\boldsymbol{\sigma}}$ is the total Cauchy’s stress, ${\bf{b}}$ denotes
the body force vector and $\rho$ is weighted average density of the mixture.
Based on Biot’s theory, the effective stress correlates with the total stress
tensor as ${\boldsymbol{\sigma}^{\prime}}={\boldsymbol{\sigma}}-\alpha
p{\bf{I}}$, where ${\boldsymbol{\sigma}^{\prime}}$ is the effective stress,
$p$ is the pore fluid pressure and ${\bf{I}}$ is the identity tensor. For an
isotropic porous medium, the scalar Biot’s coefficient $\alpha$ is defined as
$\alpha=1-\frac{{{K_{\rm{t}}}}}{{{K_{\rm{s}}}}}$, in which ${{K_{\rm{t}}}}$
and ${{K_{\rm{s}}}}$ are respectively the bulk moduli of the porous medium and
solid particles. Moreover, $\rho=(1-n){\rho_{\rm{s}}}+n{\rho_{\rm{f}}}$, where
${\rho_{\rm{s}}}$ and ${\rho_{\rm{f}}}$ respectively represent the density of
the solid and fluid phases, and $n$ is the porosity of the whole mixture.

The thermally induced strain is expressed as

${{\boldsymbol{\varepsilon}}_{T}}={\boldsymbol{\alpha}_{T}}(T-{T_{{\rm{0}}}})$
(2)

where $T$ and ${T_{{\rm{0}}}}$ are the current and reference temperatures of
the porous medium. $\boldsymbol{\alpha}_{T}$ is the thermal expansivity
tensor, which is defined for an isotropic medium as
$\boldsymbol{\alpha}_{T}={\beta_{\rm{s}}}{\bf{I}}$. Here ${{\beta_{\rm{s}}}}$
is the volumetric thermal expansion coefficient of the porous domain, which is
identical to thermal expansion coefficient of the solid constituent [55]. It
is postulated that the solid-fluid mixture is locally in thermal equilibrium
since the fluid velocity inside the pores is relatively insignificant [56].
Therefore, the constitutive equations for the effective stress can be
expressed as [2]:

${\boldsymbol{\sigma}}^{\prime}={\bf{D}}{\boldsymbol{\varepsilon}}-\left[{\beta_{\rm{s}}}{K_{\rm{t}}}(T-{T_{\text{0}}})\right]{\bf{I}}$
(3)

where ${\bf{D}}$ is the Hooke’s elasticity tensor and
${\boldsymbol{\varepsilon}}=\frac{1}{2}\left({\nabla{\bf{u}}+{\nabla^{\rm{T}}}{\bf{u}}}\right)$
is the total infinitesimal strain tensor, with ${\bf{u}}$ being the
displacement vector. In this study, internal discontinuities are considered as
impermeable boundaries with zero leak-off. In this fashion, the continuity
equation of the fluid phase in the porous matrix can be obtained by employing
the averaging theory in the porous media as follows [1]:

$\alpha\nabla\cdot{\dot{\bf{u}}}+\nabla\cdot{\dot{\bf{w}}_{\rm{f}}}+\frac{1}{{{Q_{\rm{t}}}}}\dot{p}-{\beta_{\rm{t}}}\dot{T}=0$
(4)

where
$\frac{1}{{{Q_{\rm{t}}}}}=\left({\frac{{\alpha-n}}{{{K_{\rm{s}}}}}+\frac{n}{{{K_{\rm{f}}}}}}\right)$
is the compressibility coefficient of the bulk medium. Moreover,
${\beta_{\rm{t}}}={(\alpha-n){\beta_{\rm{s}}}+n{\beta_{\rm{f}}}}$, where
${\beta_{\rm{f}}}$ denotes the fluid coefficients of thermal expansion.
${\dot{\bf{u}}}$ is the solid phase velocity, and ${\dot{\bf{w}}_{\rm{f}}}$ is
the relative velocity of the fluid phase with respect to the solid skeleton
that is expressed by means of Darcy’s law as

${\dot{\bf{w}}_{\rm{f}}}=\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}\left({-\nabla
p+{\rho_{\rm{f}}}{\bf{b}}}\right)$ (5)

in which ${{k_{\rm{f}}}}$ and ${{\mu_{\rm{f}}}}$ are the intrinsic
permeability of the domain and fluid viscosity, respectively. Combining Eqs. 4
and 5, the flow-continuity equation of the pore fluid can be recast into

$\alpha\nabla\cdot{\dot{\bf{u}}}+\nabla\cdot\left({\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}\left({-\nabla
p+{\rho_{\rm{f}}}{\bf{b}}}\right)}\right)+\frac{1}{{{Q_{\rm{t}}}}}\dot{p}-{\beta_{\rm{t}}}\dot{T}=0$
(6)

The heat transfer relation is deduced by combining Fourier’s law and the
conservation of energy equation for the porous media. In resemblance to the
fluid flow formulation, the thermal energy balance equation is derived using
the averaging theory for the energy balance of both solid and fluid phases, as
follows [57]

${\left({\rho
C}\right)_{\text{eff}}}\dot{T}+{\rho_{\rm{f}}}{C_{\rm{f}}}\left({\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}\left({-\nabla
p+{\rho_{\rm{f}}}{\bf{b}}}\right)}\right)\cdot\nabla
T-\nabla\cdot\left({{\lambda_{\text{eff}}}\nabla T}\right)=0$ (7)

where ${C_{\text{eff}}}$ and ${\lambda_{\text{eff}}}$ are the effective
specific heat capacity and thermal conductivity, respectively, defined as

$\begin{array}[]{l}{\left({\rho
C}\right)_{\text{eff}}}=(1-n){\rho_{\rm{s}}}{C_{\rm{s}}}+n{\rho_{\rm{f}}}{C_{\rm{f}}}\\\
{\lambda_{\text{eff}}}=(1-n){\lambda_{\rm{s}}}+n{\lambda_{\rm{f}}}\end{array}$
(8)

In Eq. 8, ${C_{\rm{s}}}$, ${C_{\rm{f}}}$, ${\lambda_{\rm{s}}}$ and
${\lambda_{\rm{f}}}$ are the specific heat capacities and thermal heat
conductivities of the solid and fluid phases, respectively. It is worth noting
that heat generation due to phase compression is neglected [1].

The above thermo-hydro-mechanical governing equations are subjected to the
following Dirichlet boundary conditions on the external boundaries:

$\begin{array}[]{*{20}{c}}{{\bf{u=\bar{u}}}}&{{\rm{on}}}&{{\Gamma_{u}}}\\\
{p=\bar{p}}&{{\rm{on}}}&{{\Gamma_{p}}}\\\
{T=\bar{T}}&{{\rm{on}}}&{{\Gamma_{T}}}\end{array}$ (9)

where ${{\bf{\bar{u}}}}$, ${\bar{p}}$ and ${\bar{T}}$ are the prescribed
displacement, pressure and temperature applied on the associated external
boundaries ${{\Gamma_{u}}}$, ${{\Gamma_{p}}}$ and ${{\Gamma_{T}}}$,
respectively. The Neumann boundary conditions on the contrary are described as

$\begin{array}[]{*{20}{c}}{{\boldsymbol{\sigma}}\cdot{{\bf{n}}_{\Gamma}=\bar{\bf
t}}}&{{\rm{on}}}&{{\Gamma_{\rm{t}}}}\\\
{{{\dot{\bf{w}}}}\cdot{{\bf{n}}_{\Gamma}}={{\bar{q}}_{\rm{f}}}}&{{\rm{on}}}&{{\Gamma_{{q_{p}}}}}\\\
{-\left({{\lambda_{\text{eff}}}\nabla
T}\right)\cdot{{\bf{n}}_{\Gamma}}={{\bar{q}}_{T}}}&{{\rm{on}}}&{{\Gamma_{{q_{T}}}}}\end{array}$
(10)

where ${{{\bf{n}}_{\Gamma}}}$ is the normal unit vector of the external
boundary $\Gamma$. In addition, ${{\bar{\bf t}}}$, ${{{\bar{q}}_{\rm{f}}}}$
and ${{{\bar{q}}_{T}}}$ are respectively the applied traction and, fluid and
heat fluxes along the corresponding external boundaries ${{\Gamma_{\rm{t}}}}$,
${{\Gamma_{{q_{p}}}}}$ and ${{\Gamma_{{q_{T}}}}}$. Note the external
boundaries are defined such that
${\Gamma_{u}}\cup{\Gamma_{\rm{t}}}={\Gamma_{p}}\cup{\Gamma_{{q_{p}}}}={\Gamma_{T}}\cup{\Gamma_{{q_{T}}}}=\Gamma$
and
${\Gamma_{u}}\cap{\Gamma_{\rm{t}}}={\Gamma_{p}}\cap{\Gamma_{{q_{p}}}}={\Gamma_{T}}\cap{\Gamma_{q_{T}}}=\emptyset$
hold.

The boundary conditions associated with the impermeable adiabatic internal
discontinuities of the domain subject to fluid flow, thermal flux and contact
constraints are expressed as

$\begin{array}[]{*{20}{c}}{{\boldsymbol{\sigma}}\cdot{{\bf{n}}_{{{\Gamma}_{\rm{d}}}}}{=}{{\bf{t}}^{{\rm{c}}}}}&{{\rm{on}}}&{{\Gamma_{\rm{c}}}\subset{\Gamma_{\rm{d}}}}\\\
{\left[\kern-1.49994pt\left[{{{\dot{\bf{w}}}_{\rm{f}}}}\right]\kern-1.49994pt\right]\cdot{{\bf{n}}_{{{\bf{\Gamma}}_{\rm{d}}}}}=0}&{{\rm{on}}}&{{\Gamma_{\rm{d}}}}\\\
{-\left[\kern-1.49994pt\left[{{\lambda_{{\rm{eff}}}}\nabla
T}\right]\kern-1.49994pt\right]\cdot{{\bf{n}}_{{{\Gamma}_{\rm{d}}}}}={q_{T}^{\rm{c}}}}&{{\rm{on}}}&{{\Gamma_{\rm{c}}}\subset{\Gamma_{\rm{d}}}}\end{array}$
(11)

in which ${\Gamma_{\rm{d}}}$ and ${{\bf{n}}_{{{\Gamma}}_{\rm{d}}}}$ are the
discontinuity interface and its normal unit vector, respectively. The notation
$\left[\kern-1.49994pt\left[\Re\right]\kern-1.49994pt\right]$ represents the
jump in the field
$\Re\in\text{\\{}{\bf{u}},{{{{\bf{w}}}_{\rm{f}}}},T\text{\\}}$ across the
discontinuity faces. Moreover, ${{\bf{t}}^{\rm{c}}}$ and $q_{T}^{\rm{c}}$
represent contact traction vector and heat flux, respectively, that are
imposed on the contact boundary ${\Gamma_{\rm{c}}}$. Given the normal gap
${g_{\rm{N}}}=\left[\kern-1.49994pt\left[{{u_{\rm{N}}}}\right]\kern-1.49994pt\right]=\left[\kern-1.49994pt\left[{{{{\bf{u}}}}}\right]\kern-1.49994pt\right]\cdot{{\bf{n}}_{{{\Gamma}_{\rm{d}}}}}$
(i.e., the normal component of the displacement jump), the Kuhn-Tucker
criterion is employed to elaborate the impenetration condition across the
discontinuity faces as follows

${g_{\rm{N}}}\geq 0,\,\,t_{\rm{N}}^{{\rm{c}}}\leq
0,\,\,{g_{\rm{N}}}t_{\rm{N}}^{{\rm{c}}}=0$ (12)

where
$t_{\rm{N}}^{\rm{c}}={{\bf{t}}^{\rm{c}}}\cdot{\textbf{n}_{{\Gamma_{\rm{d}}}}}$
is the normal contact traction along the discontinuity. Similarly,
$q_{T}^{\rm{c}}={h_{{\rm{cont}}}}\left[\kern-1.49994pt\left[T\right]\kern-1.49994pt\right]$
is the heat flux imposed on the mechanical contact surface, with
${h_{{\rm{cont}}}}$ and
$\left[\kern-1.49994pt\left[T\right]\kern-1.49994pt\right]$ being the thermal
conductivity coefficients and temperature jump across the contact interface,
respectively [35]. Note that since the discontinuity interfaces are stipulated
as impermeable boundaries, no fluid exchange and associated convective heat
transfer may occur between the discontinuity and the surrounding porous
medium. Furthermore, the conductive heat transfer is assumed to be negligible
in the case of discontinuity opening. This is a valid assumption wherever the
thermal conductivity is close to zero, e.g., heat insulation of the void/air
inside the cavity, in the absence of any fluids [58].

###  2.2 Weak form and XFEM discretization

In order to derive the weak form of the governing equations, i.e., Eqs. 1, 6
and 7, three sets of costume tailored tests functions
${\boldsymbol{\eta}}_{u}$, $\eta_{p}$ and $\eta_{T}$ are respectively adopted,
which all satisfy the homogeneous Dirichlet boundary conditions. Multiplying
the governing equations by the test functions and employing the Gauss-
Divergence theorem yields the weak form equations as

$\begin{array}[]{l}\int\limits_{\Omega}{{\nabla^{\rm{s}}}{{\boldsymbol{\eta}}_{u}}:{\boldsymbol{\sigma}}{\rm{d}}\Omega}=\int\limits_{\Omega}{{{\boldsymbol{\eta}}_{u}}\cdot\rho{\bf{b}}{\rm{d}}\Omega}+\int\limits_{{\Gamma_{\rm{t}}}}{{{\boldsymbol{\eta}}_{u}}\cdot{\bf{\bar{t}}}{\rm{d}}\Gamma}\\\
\int\limits_{\Omega}{{\eta_{p}}\alpha\nabla\cdot{\dot{\bf{u}}}{\rm{d}}\Omega}+\int\limits_{\Omega}{\nabla{\eta_{p}}\cdot\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}\nabla
p{\rm{d}}\Omega}+\int\limits_{\Omega}{{\eta_{p}}\frac{1}{{{Q_{\rm{t}}}}}\dot{p}{\rm{d}}\Omega}-\int\limits_{\Omega}{{\eta_{p}}{\beta_{\rm{t}}}\dot{T}{\rm{d}}\Omega}=\int\limits_{{\Gamma_{{q_{p}}}}}{{\eta_{p}}{{\bar{q}}_{\rm{f}}}{\rm{d}}\Gamma}-\int\limits_{\Omega}{{\nabla{\eta_{p}}}\cdot\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}{\rho_{\rm{f}}}{\bf{b}}\rm{d}\Omega}\\\
\int\limits_{\Omega}{{\eta_{T}}{{(\rho
C)}_{{\rm{eff}}}}\dot{T}{\rm{d}}\Omega}+\int\limits_{\Omega}{{\eta_{T}}\left[{{\rho_{\rm{f}}}{C_{\rm{f}}}\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}(-\nabla
p+{\rho_{\rm{f}}}{\bf{b}})}\right]\cdot\nabla
T{\rm{d}}\Omega}+\int\limits_{\Omega}{\nabla{\eta_{T}}\cdot{\lambda_{{\rm{eff}}}}\nabla
T{\rm{d}}\Omega}=\int\limits_{{\Gamma_{{q_{T}}}}}{{\eta_{T}}{{\bar{q}}_{T}}{\rm{d}}\Gamma}\end{array}$
(13)

As depicted in Fig. 1, the displacement, pressure and temperature fields are
discontinuous across ${\Gamma_{\rm{d}}}$. Thus, the corresponding discretized
XFEM fields can be expressed as

$\begin{array}[]{l}{\bf{u}}({\bf{x}})={{\bf{u}}^{{\rm{c}}}}({\bf{x}})+H_{\Gamma_{d}}(\varphi(\mathbf{x})){{\bf{u}}^{{\rm{d}}}}({\bf{x}})\\\
p({\bf{x}})={p^{{\rm{c}}}}({\bf{x}})+H_{\Gamma_{d}}(\varphi(\mathbf{x})){p^{{\rm{d}}}}({\bf{x}})\\\
T({\bf{x}})={T^{{\rm{c}}}}({\bf{x}})+H_{\Gamma_{d}}(\varphi(\mathbf{x})){T^{{\rm{d}}}}({\bf{x}})\end{array}$
(14)

where ${{\bf{u}}^{{\rm{c}}}}$, ${p^{{\rm{c}}}}$, ${T^{{\rm{c}}}}$ and
${{\bf{u}}^{{\rm{d}}}}$, ${p^{{\rm{d}}}}$, ${T^{{\rm{d}}}}$ are associated
with the standard (continuous) and enriched (discontinuous) approximations of
the displacement, pressure and temperature fields, respectively. Moreover,
$H_{\Gamma_{d}}(\varphi(\mathbf{x}))=\mathbb{H}_{\Gamma_{d}}(\varphi(\mathbf{x}))-\mathbb{H}_{\Gamma_{d}}(\varphi(\mathbf{x})^{I})$
is the shifted Heaviside enrichment function that reproduces strong
discontinuity across ${\Gamma_{\rm{d}}}$ [59], in which

$\mathbb{H}_{\Gamma_{d}}(\varphi(\mathbf{x}))=\left\\{\begin{matrix}1&\varphi(\mathbf{x})\geq
0\\\ -1&\varphi(\mathbf{x})<0\end{matrix}\right.$ (15)

is the Heaviside function and $\varphi(\mathbf{x})$ represents the signed
distance function with respect to discontinuity interface $\Gamma_{d}$ [52].
According to Eq. 14, the infinitesimal strain tensor as well as the pressure
and temperature gradients can be approximated as

$\begin{array}[]{l}{\boldsymbol{\varepsilon}}({\bf{x}})={\nabla^{\rm{s}}}{\bf{u}}({\bf{x}})={\nabla^{\rm{s}}}{{\bf{u}}^{{\rm{c}}}}({\bf{x}})+{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\nabla{{\bf{u}}^{{\rm{d}}}}({\bf{x}})+{\delta_{{\Gamma_{\rm{d}}}}}{({{\bf{u}}^{{\rm{d}}}}({\bf{x}})\otimes{{\bf{n}}_{{\Gamma_{\rm{d}}}}})^{\rm{s}}}\\\
\nabla
p({\bf{x}})=\nabla{p^{{\rm{c}}}}({\bf{x}})+{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\nabla{p^{{\rm{d}}}}({\bf{x}})+{\delta_{{\Gamma_{\rm{d}}}}}p({\bf{x}}){{\bf{n}}_{{\Gamma_{\rm{d}}}}}\\\
\nabla
T({\bf{x}})=\nabla{T^{{\rm{c}}}}({\bf{x}})+{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\nabla{T^{{\rm{d}}}}({\bf{x}})+{\delta_{{\Gamma_{\rm{d}}}}}T({\bf{x}}){{\bf{n}}_{{\Gamma_{\rm{d}}}}}\end{array}$
(16)

where ${\nabla^{\rm{s}}}$ and ${(\cdot)^{\rm{s}}}$ are the symmetric parts of
the spatial gradient operator and tensor, respectively, and
${\delta_{{\Gamma_{\rm{d}}}}}$ is the Dirac’s delta function across
$\Gamma_{\rm{d}}$.

Similar to the trial functions, the test functions can also be decomposed into
continuous and discontinuous parts as follows

$\begin{array}[]{l}{\boldsymbol{\eta}}_{u}={\boldsymbol{\eta}}_{u}^{{\rm{c}}}+{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}})){\boldsymbol{\eta}}_{u}^{{\rm{d}}}\\\
\eta_{p}=\eta_{p}^{{\rm{c}}}+{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\eta_{p}^{{\rm{d}}}\\\
\eta_{T}=\eta_{T}^{{\rm{c}}}+{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\eta_{T}^{{\rm{d}}}\end{array}$
(17)

The presented XFEM formulation is fully compatible with the inherent structure
of COMSOL Multiphysics [50, 53]. By substituting the above decomposition in
Eq. 13, the weak forms of the continuous parts of the governing equations can
be obtained as

$\begin{array}[]{l}\int\limits_{\Omega}{{\nabla^{\rm{s}}}{{\boldsymbol{\eta}}_{u}^{{\rm{c}}}}:{\boldsymbol{\sigma}}{\rm{d}}\Omega}=\int\limits_{\Omega}{{{\boldsymbol{\eta}}_{u}^{{\rm{c}}}}\cdot\rho{\bf{b}}{\rm{d}}\Omega}+\int\limits_{{\Gamma_{\rm{t}}}}{{{\boldsymbol{\eta}}_{u}^{{\rm{c}}}}\cdot{\bf{\bar{t}}}{\rm{d}}\Gamma}\\\
\int\limits_{\Omega}{{\eta_{p}^{{\rm{c}}}}\alpha\nabla\cdot{\dot{\bf{u}}}{\rm{d}}\Omega}+\int\limits_{\Omega}{\nabla{\eta_{p}^{{\rm{c}}}}\cdot\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}\nabla
p{\rm{d}}\Omega}+\int\limits_{\Omega}{{\eta_{p}^{{\rm{c}}}}\frac{1}{{{Q_{\rm{t}}}}}\dot{p}{\rm{d}}\Omega}-\int\limits_{\Omega}{{\eta_{p}^{{\rm{c}}}}{\beta_{\rm{t}}}\dot{T}{\rm{d}}\Omega}=\int\limits_{{\Gamma_{{q_{p}}}}}{{\eta_{p}^{{\rm{c}}}}{{\bar{q}}_{\rm{f}}}{\rm{d}}\Gamma}-\int\limits_{\Omega}{{\nabla}{\eta_{p}^{{\rm{c}}}}\cdot\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}{\rho_{\rm{f}}}{\bf{b}}\rm{d}\Omega}\\\
\int\limits_{\Omega}{{\eta_{T}^{{\rm{c}}}}{{(\rho
C)}_{{\rm{eff}}}}\dot{T}{\rm{d}}\Omega}+\int\limits_{\Omega}{{\eta_{T}^{{\rm{c}}}}\left[{{\rho_{\rm{f}}}{C_{\rm{f}}}\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}(-\nabla
p+{\rho_{\rm{f}}}{\bf{b}})}\right]\cdot\nabla
T{\rm{d}}\Omega}+\int\limits_{\Omega}{\nabla{\eta_{T}^{{\rm{c}}}}\cdot{\lambda_{{\rm{eff}}}}\nabla
T{\rm{d}}\Omega}=\int\limits_{{\Gamma_{{q_{T}}}}}{{\eta_{T}^{{\rm{c}}}}{{\bar{q}}_{T}}{\rm{d}}\Gamma}\end{array}$
(18)

and the discontinuous enriched parts as

$\begin{array}[]{l}\int\limits_{\Omega}{\left[{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}})){\nabla^{\rm{s}}}{\boldsymbol{\eta}}_{u}^{{\rm{d}}}}\right]:{\boldsymbol{\sigma}}{\rm{d}}\Omega}=\int\limits_{\Omega}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}})){\boldsymbol{\eta}}_{u}^{{\rm{d}}}\cdot\rho{\bf{b}}{\rm{d}}\Omega}+\int\limits_{{\Gamma_{\rm{t}}}}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}})){\boldsymbol{\eta}}_{u}^{{\rm{d}}}\cdot{\bf{\bar{t}}}{\rm{d}}\Gamma}-\int\limits_{{\Gamma_{\rm{c}}}}{{\boldsymbol{\eta}}_{u}^{{\rm{d}}}\cdot{{\bf{t}}^{{\rm{c}}}}{\rm{d}}\Gamma}\\\
\int\limits_{\Omega}{\left[{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\eta_{p}^{{\rm{d}}}}\right]\alpha\nabla\cdot{\dot{\bf{u}}}{\rm{d}}\Omega}+\int\limits_{\Omega}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\nabla\eta_{p}^{{\rm{d}}}\cdot\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}\nabla
p{\rm{d}}\Omega}+\int\limits_{\Omega}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\eta_{p}^{{\rm{d}}}\frac{1}{{{Q_{\rm{t}}}}}\dot{p}{\rm{d}}\Omega}\\\
-\int\limits_{\Omega}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\eta_{p}^{{\rm{d}}}{\beta_{\rm{t}}}\dot{T}{\rm{d}}\Omega}=\int\limits_{{\Gamma_{{q_{p}}}}}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\eta_{p}^{{\rm{d}}}{{\bar{q}}_{\rm{f}}}{\rm{d}}\Gamma}-\int\limits_{\Omega}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}})){\nabla}\eta_{p}^{{\rm{d}}}\cdot\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}{\rho_{\rm{f}}}{\bf{b}}\rm{d}\Omega}\\\
\int\limits_{\Omega}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\eta_{T}^{{\rm{d}}}{{(\rho
C)}_{{\rm{eff}}}}\dot{T}{\rm{d}}\Omega}+\int\limits_{\Omega}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\eta_{T}^{{\rm{d}}}\left[{{\rho_{\rm{f}}}{C_{\rm{f}}}\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}(-\nabla
p+{\rho_{\rm{f}}}{\bf{b}})}\right]\cdot\nabla T{\rm{d}}\Omega}\\\
+\int\limits_{\Omega}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\nabla\eta_{T}^{{\rm{d}}}\cdot{\lambda_{{\rm{eff}}}}\nabla
T{\rm{d}}\Omega}=\int\limits_{{\Gamma_{{q_{T}}}}}{{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\eta_{T}^{{\rm{d}}}{{\bar{q}}_{T}}{\rm{d}}\Gamma}-\int\limits_{{\Gamma_{\rm{c}}}}{\eta_{T}^{{\rm{d}}}{q_{T}^{\rm{c}}}{\rm{d}}\Gamma}\end{array}$
(19)

In Eq. 19, the integration domain is limited to the support of
$H_{\Gamma_{d}}(\varphi(\mathbf{x}))$ i.e., the enriched zone
${\Omega_{\rm{h}}}$ detected by the signed distance function.

In order to solve the integral Eqs. 18 and 19, the spatial discretization is
employed within the context of XFEM. Considering a Bubnov-Galerkin scheme,
both the trial and test functions are discretized using ${\rm{C}^{0}}$
continuous shape functions ${N_{i}}({\bf{x}})$ for the primary variables
$\aleph=\left\\{{{\bf{u}},p,T}\right\\}$ as

$\left\\{{\begin{array}[]{*{20}{c}}{{\aleph^{{\rm{c}}}}({\bf{x}})=\sum\limits_{i\in{m_{{\rm{std}}}}}{{N_{i}}({\bf{x}}){{\hat{X}}_{i}}}}&{{\rm{in}}\,\Omega}\\\
{{\aleph^{{\rm{d}}}}({\bf{x}})=\sum\limits_{i\in{m_{{\rm{enr}}}}}{{N_{i}}({\bf{x}}){{\tilde{X}}_{i}}}}&{{\rm{in}}\,{\Omega_{\rm{h}}}}\end{array}}\right.$
(20)

where $\hat{X}=\left\\{{{\bf{\hat{u}}},\hat{p},\hat{T}}\right\\}$ and
$\tilde{X}=\left\\{{{\bf{\tilde{u}}},\tilde{p},\tilde{T}}\right\\}$ are the
vectors of nodal displacements, pressure and temperature for the standard and
enriched parts, respectively. Moreover, ${m_{{\rm{std}}}}$ and
${m_{{\rm{enr}}}}$ denote the sets of standard and enriched nodes,
respectively.

##  3 Implementation in COMSOL

###  3.1 Overall framework

COMSOL is a general-purpose multi-field program that renders sophisticated
multiphysics modeling capabilities. In this software, a combination of
available built-in physics interfaces can be incorporated in conjunction with
user-defined physics to investigate formidable multiphysics problems [60, 50].
The thermo-hydro-mechanical coupling analysis of fluid flow in deformable
porous media, considered here, involves three distinct physics that take
account of the mechanical deformation, fluid flow and heat transfer. In
COMSOL, the “Solid Mechanics” physics interface facilitates the most general
toolkit to perform continuum based structural analysis by solving the
equations of motion endowed with suitable constitutive material behaviour. The
fluid flow through the porous domain can be simulated via the “Porous Media
and Subsurface Flow” module. In the case of low-velocity flows, which is
commonly the case in geomechanics, the “Darcy’s Law” physics interface is
employed to incorporate the flow-continuity equations within the porous
medium. The “Heat Transfer” module offers heat transfer analysis in the host
domain. Among the available interfaces in this module, the “Heat Transfer in
Porous Media” physics interface is utilized to include the convection-
diffusion equation through both solid matrix and pore fluid phases. In the
conventional FE modeling approach, the discontinuity in the solution field
needs to be explicitly modeled, while the mesh must conform to the internal
boundaries for all physics interfaces. In contrast, based on the XFEM
approach, for which the theoretical background is presented in the previous
section, one can circumvent the requirement of mesh alignment with the
internal interfaces. Extending the recent work of the authors on the XFEM
implementation of fracturing of solid formations in COMSOL [50], three
distinct COMSOL interfaces are introduced, one for each of the physics
involved. The solution variables associated with these physics interfaces are
recast in the set
$\left\\{{{\aleph^{{\rm{c}}}},{\aleph^{{\rm{d}}}}}\right\\}$, where
${\aleph^{{\rm{c}}}}=\left\\{{{{\bf{u}}^{{\rm{c}}}},{p^{{\rm{c}}}},{T^{{\rm{c}}}}}\right\\}$
constitutes the continuous parts of the solution fields, and
${\aleph^{{\rm{d}}}}=\left\\{{{{\bf{u}}^{{\rm{d}}}},{p^{{\rm{d}}}},{T^{{\rm{d}}}}}\right\\}$
encompasses the discontinuous parts. For the sake of brevity, the allocated
standard parts of the solution domain corresponding to the “Solid Mechanics”,
“Darcy’s Law” and “Heat Transfer in Porous Media” physics interfaces are
denoted by SMstd, DLstd and HTstd, respectively, and their enriched
counterparts are indicated by SMenr, DLenr and HTenr, respectively, in the
remainder of the paper.

Several important features such as integration tools and analytical functions
are inherently incorporated within COMSOL software. However, additional
external functions are imperative to implement in order to successfully
execute the various functionalities required for the XFEM simulation, such as
the level-set analysis in the pre-processing stage. These additional
functions/subroutines are developed in MATLAB and linked to the model through
the COMSOL Live-link for MATLAB software extension [50]. In the following
sections, the pre-processing step, physics interfaces setup (i.e., XFEM
enrichment, thermo-hydro-mechanical coupling, and thermo-mechanical contact),
numerical integration scheme and the numerical solution strategy are
elaborated in detail.

###  3.2 Pre-processing: identification of crack geometry

In XFEM simulations, a pre-processing phase is required prior to the analysis
in order to identify the set of bisected elements by the crack interface. To
this end, a conventional FEM mesh is first generated by neglecting the
geometry of internal interfaces (i.e., cracks) in COMSOL. The mesh elemental
and nodal information is then exported as a “.mphtxt” file to undertake the
level-set analysis in a separate MATLAB script file “preprocess.m”.
Accordingly, a geometrical search algorithm is employed across all elements’
nodes and edges to figure out their configuration with respect to the
discontinuity interfaces. A list consisting of all enriched elements and nodes
is yielded, which in turn, is employed to initialize the XFEM physics
interfaces over the enriched zone ${\Omega_{\rm{h}}}$. In the case of crack
growth, the pre-processing step is repeated in each propagation step to
account for the alterations in the discontinuity interfaces. This is performed
by utilizing two MATLAB functions; “phi.m”, which evaluates the Heaviside
function at any point of interest (e.g., Gauss points and nodes) and modifies
the strain field, pressure and temperature gradients in corresponding physics
interfaces, i.e., SMenr, DLenr and HTenr, respectively. The second function
entitled “interpol.m”, identifies the enriched zone by means of an
interpolation function. This function attributes one, to the enriched nodal
points, and zero, to the remainder of nodal points. During the course of the
analysis, the output variable $\psi$ is interpolated using the MATLAB internal
function scatterInterpolant, which is scripted within the external
function“interpol.m”. In order to provide access to the element
connectivities, nodal coordinates, interpolation values of nodal points and
crack-tip coordinates at any stage of the analysis, all are stored in distinct
MAT files. The above procedure is designated so as to overcome the inherent
restriction of no access to the elements and nodes information encountered
during the course of analysis in COMSOL. In this way, by recalling the MAT
files at any stage of the analysis, the required preprocessing data is
geometrically inferred from the background mesh. It is noteworthy that COMSOL
allows the user to associate a single output with each MATLAB function call.
As a result, the number of user-defined MATLAB functions must match the
variables required during the solution. In addition, all input and output
vectors associated with each of the MATLAB functions must necessarily have
identical array sizes.

Remark 1. In certain circumstances where crack interfaces are aligned with the
element edges (e.g., stationary or mode I propagating cracks), the level set
updating process can be performed by relying on built-in features of COMSOL in
lieu of user-defined MATLAB functions. This can be achieved by taking
advantage of a series of explicit expressions that specify the bisected
elements as well as the Heaviside values at Gauss integration points. This
alleviates the interactions with LiveLink for MATLAB throughout the analysis
and therefore, enhances the computational efficiency considerably.

###  3.3 Module set-up

####  3.3.1 XFEM enrichment

According to the weak formulation presented in section 2.2, two distinct
modules are set up as per physics interfaces (i.e., ”Solid Mechanics”,
”Darcy’s Law” and ”Heat Transfer in Porous Media”) in order to take account of
the standard and enriched terms of the corresponding field variables (i.e.,
$\bf{u}$, $p$ and $T$, respectively). For the continuous parts of the field
variables in standard physics interfaces, no special modification is required
for the definition of the strain field
${{\boldsymbol{\varepsilon}}^{{\rm{std}}}}({\bf{x}})$, pressure
$\nabla{p^{{\rm{std}}}}({\bf{x}})$ and temperature gradients
$\nabla{T^{{\rm{std}}}}({\bf{x}})$, thus

$\begin{array}[]{l}{{\boldsymbol{\varepsilon}}^{{\rm{std}}}}({\bf{x}})={\nabla^{\rm{s}}}{{\bf{u}}^{{\rm{c}}}}({\bf{x}})\\\
\nabla{p^{{\rm{std}}}}({\bf{x}})=\nabla{p^{{\rm{c}}}}({\bf{x}})\\\
\nabla{T^{{\rm{std}}}}({\bf{x}})=\nabla{T^{{\rm{c}}}}({\bf{x}})\end{array}$
(21)

On the other hand, the enriched physics interfaces require modifications to
take into account the discontinuities in the solution domain. This task is
undertaken by adopting the enrichment function provided by the MATLAB function
“phi.m” (i.e., $H_{\Gamma_{d}}(\varphi(\mathbf{x}))$, in here) as

$\begin{array}[]{l}{{\boldsymbol{\varepsilon}}^{{\rm{enr}}}}({\bf{x}})={\nabla^{\rm{s}}}{{\bf{u}}^{{\rm{d}}}}({\bf{x}})\times{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\\\
\nabla{p^{{\rm{enr}}}}({\bf{x}})=\nabla{p^{{\rm{d}}}}({\bf{x}})\times{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\\\
\nabla{T^{{\rm{enr}}}}({\bf{x}})=\nabla{T^{{\rm{d}}}}({\bf{x}})\times{H_{{\Gamma_{\rm{d}}}}}(\varphi({\bf{x}}))\end{array}$
(22)

where ${{\boldsymbol{\varepsilon}}^{{\rm{enr}}}}({\bf{x}})$,
$\nabla{p^{{\rm{enr}}}}({\bf{x}})$ and $\nabla{T^{{\rm{enr}}}}({\bf{x}})$ are
respectively the enriched contributions of the strain field, pressure and
temperature gradients. To apply such modifications, the user needs to enable
“Equation View” option in COMSOL, which facilitates the access to the
definition of all variables related to each physics interface.

In order to ensure the consistency of the solution across all the physics
involved, further modifications need to be introduced in the definition of
stress, velocity and heat flux fields. For this purpose, the overall (i.e.,
contributions from both standard and enriched fields) effective stress
${\boldsymbol{\bar{\sigma}}}$, fluid velocity
$\dot{\bar{\mathbf{w}}}_{\rm{f}}$, and conductive
$\bar{\mathbf{q}}_{\text{cd}}$ and convective $\bar{\mathbf{q}}_{\text{cv}}$
heat transfer vectors, which are respectively associated with Cauchy stress,
Darcy’s velocity field and Conductive and Convective heat transfer vectors,
are modified as

$\begin{array}[]{l}{\boldsymbol{\bar{\sigma}}}={\bf{D}}:({{\boldsymbol{\varepsilon}}^{{\rm{std}}}}+{{\boldsymbol{\varepsilon}}^{{\rm{enr}}}})\\\
\dot{\bar{\mathbf{w}}}_{\rm{f}}=\frac{{{k_{\rm{f}}}}}{{{\mu_{\rm{f}}}}}(-\nabla{p^{{\rm{std}}}}-\nabla{p^{{\rm{enr}}}}+{\rho_{\rm{f}}}({\bf{b}}-\ddot{\bf{u}}))\\\
\bar{\mathbf{q}}_{\text{cd}}={\lambda_{{\rm{eff}}}}\nabla
T={\lambda_{{\rm{eff}}}}(\nabla{T^{{\rm{std}}}}+\nabla{T^{{\rm{enr}}}})\\\
\bar{\mathbf{q}}_{\text{cv}}={\rho_{\rm{f}}}{C_{\rm{f}}}\dot{\bar{\mathbf{w}}}_{\rm{f}}\end{array}$
(23)

It is noteworthy that the above modifications must be identically applied to
both standard and enriched sets of physics interfaces.

The enriched physics interfaces need to be merely introduced and solved over
the enrichment zone ${\Omega_{\rm{h}}}$ and its associated DOFs. However, in
COMSOL the physics interfaces are essentially defined across the whole domain.
In order to overcome such restriction, the field variable $\psi$ (see section
3.2) is applied to prescribe the primary variables of the enriched physics
interfaces such that,

${\aleph^{{\rm{d}}}}=\left\\{{\begin{array}[]{*{20}{c}}{{\aleph^{{\rm{d}}}}}&{{\rm{where}}\,\psi=1}\\\
0&{{\rm{otherwise}}}\end{array}}\right.$ (24)

Using the above definition, the solution domain associated with the enriched
physics interfaces, which is ubiquitously defined over the entire domain, can
be restricted to the enriched zone. This task is performed by using the
“Pointwise Constraints” from the “Domain Constraints” by selecting the same
discretization settings as per their corresponding physics (e.g., shape
function type, element order, etc). This, in turn, results in the elimination
of the redundant DOFs that are located outside the desired enrichment zone.
This strategy significantly improves the computational efficiency of the
model, while problematic singularities in the stiffness matrix emanating from
the incorrect definition of the enriched DOFs over the un-enriched zone are
avoided.

####  3.3.2 Thermo-hydro-mechanical coupling process

COMSOL Multiphysics features excellent capabilities in coupling multiple
physical phenomena. This task can either be performed by selecting predefined
multiphysics interfaces or by administering the coupling process manually. In
the latter, the user controls the multiphysics coupling between the physics
interfaces by manually specifying how each physics affects another [60]. This
technique is particularly advantageous whenever there is no preconfigured
multiphysics interface designated for the physics involved, which is the case
for the XFEM implementation of the thermo-hydro-mechanical coupling theme
proposed here. Each of the physics involved (i.e., geomechanical deformation,
fluid flow and heat transfer processes) essentially interacts with another,
which in turn generates a coupled multiphysics theme as demonstrated in Fig.
2.

Figure 2:  Thermo-hydro-mechanical coupling flowchart.

\- Thermo-mechanical coupling:

Thermally induced strain is accommodated within the Solid Mechanics physics
interface by the use of the “Thermal Expansion” subnode [61]. It is noteworthy
that this term must only be incorporated into the SMstd physics interface, to
ensure the thermal strains are applied consistently over the entire strain
field.

\- Hydro-mechanical coupling:

The poroelastic coupling terms are defined in accordance with Biot’s
poroelasticity theory [62]. In this sense, the effective stress principle is
implemented via the “External Stress” subnodes in both SMstd and SMenr, where
the “stress input” option is set to the “pore pressure”. Moreover, the
variations of the volumetric strain, that directly contribute to the changes
in pore space (i.e., the first term in Eq. 6), are included by means of “Mass
Source” subnodes in the DLstd and DLenr physics interfaces.

\- Hydro-thermal coupling:

The hydro-thermal coupling theme is automatically captured by activating the
“Nonisothermal Flow” within HTstd and HTenr physics interfaces, with no need
to introduce any additional flow and/or heat coupling terms manually [63].

####  3.3.3 Thermo-mechanical contact

COMSOL Multiphysics is endowed with mechanical and thermal contact modeling
algorithms. However, these built-in capabilities can be elaborated only on
pre-determined contact surfaces, which need to be explicitly introduced in the
problem configuration. To facilitate the imposition of XFEM contact
constraints, the $\psi$ field is employed such that the solution domain is
delimited to the enriched zone ${\Omega_{\rm{h}}}$, as explained in section
3.2. Nevertheless, in COMSOL Multiphysics, it is not possible to perform a
numerical integration over a non-predefined geometrical contact surface. This
restriction is circumvented by adding a domain “Weak Contribution” node in
SMenr and HTenr physics interfaces, whereas a penalty contact algorithm is
elaborated to account for contact traction and heat flux in accordance to
Kuhn-Tucker inequalities (i.e., Eq. 12). Moreover, the line integral of the
contact contributions along the contact surface is converted to an equivalent
domain integration over the enriched zone $\Omega_{\rm{h}}$ through the
introduction of a special Dirac delta function $\delta$ over the discontinuity
interface. This process is illustrated schematically in Fig. 3, which pertains
to the case of a horizontal discontinuity interface, supposedly, at a vertical
distance $y_{0}$ from the origin as

$\int\limits_{{\Gamma_{\rm{d}}}}{\Im\cdot{{\bf{n}}_{{{\Gamma}}_{\rm{d}}}}d\Gamma}=\int\limits_{{\Omega_{\rm{h}}}}{\Im\delta(y-{y_{0}})\cdot{{\bf{n}}_{{{\Gamma}}_{\rm{d}}}}d\Omega}$
(25)

where $\Im$ is the normal contact force/heat flux. In this way, the domain
integration over ${{\Omega_{\rm{h}}}}$ (i.e., 2D/3D) is equivalent to the line
integral along the crack dimension ${{\Gamma_{\rm{d}}}}$ (i.e., 1D/2D).

(a)

(b)

Figure 3: a) Numerical line integration of an arbitrary function $\Im$ along
the discontinuity; b) Its equivalent domain integral over the enriched zone by
using Dirac delta function $\delta$.

###  3.4 Numerical integration

The piece-wise continuous polynomials that are employed in the conventional
FEM are accurately integrated by means of lower-order Gauss quadrature rules.
In contrast, due to the presence of various discontinuities and/or
singularities in the primary field variables and, possibly, their
corresponding derivatives, XFEM requires more precise techniques to enable the
accurate evaluation of integrals over the enriched region. In this regard,
high-order Gauss integration schemes, rectangular sub-griding and element
partitioning (e.g., triangular/rectangular) are most frequently employed in
fracture analysis [52]. Nevertheless, only the first approach is viable in
COMSOL and, therefore, it is applied for the numerical simulations.

Remark 2. In order to ensure that an ill-posed and/or singular stiffness
matrix is avoided in the XFEM analysis, a sufficient number of Gauss
integration points must exist at either side of the discontinuity interfaces.
To this end, an additional criterion is introduced so as to examine the
presence of a minimum size of the support domain associated with each node,
for a certain integration order. Accordingly, a node is enriched only if its
support domain i.e., $A^{+}/(A^{+}+A^{-})$ or $A^{-}/(A^{+}+A^{-})$ exceeds a
predefined tolerance $\delta_{\rm{s}}$, where $A^{+}$ and $A^{-}$ are the
bisected areas of each element formed by the discontinuity interface. This
criterion, evidently, avoids the mathematical redundancy of the solution [64].

###  3.5 Solution procedures and analysis strategies

The process of solving a discretized boundary value problem in COMSOL is a
hierarchy that is executed from the “Study” node [65]. From several available
study options in COMSOL, “Stationary” and “Time Dependent” subnodes are among
the frequently used ones for the simulation of porous media, that pertain to
the nature of the numerical analysis (i.e., quasi-static/steady-state and
dynamic/transient, respectively). Whenever studying quasi-static crack
propagation analysis, the prescribed load is applied incrementally. Meanwhile,
the crack interfaces are deemed to propagate once a certain crack growth
criterion is met. This task is performed in COMSOL by using the “Auxiliary
Sweep” section of the “Stationary” subnode, whereas the applied incremental
load is taken as the sweep parameter. In contrast, the “Time Dependent” study
option is inherently history-dependent and thus, the sweep strategy is not
applicable. For further detail on the XFEM implementation of crack growth in
COMSOL, refer to [50].

COMSOL provides two options for solving the nonlinear system of equations
i.e., fully coupled and segregated schemes. Depending on the nonlinearity of
the equations and the degree of coupling, the methods can be employed
interchangeably which, in turn, provides great flexibility for multiphysics
problems in terms of accuracy and efficiency. In this study, the fully coupled
implicit solver is applied in conjunction with Generalized alpha time
integration method, which is a finite difference based scheme, to execute the
transient time-dependent thermo-hydro-mechanical analysis. The step-by-step
implementation procedure adopted for the XFEM thermo-hydro-mechanical modeling
in COMSOL is concisely presented in Algorithm 1, which involves: global/local
variables and geometry definitions, selection of standard and enriched physics
interfaces along with the associated modifications, domain discretization and
mesh generation, and numerical solution technique.

Algorithm 1 Step-by-step implementation of XFEM thermo-hydro-mechanical model
in COMSOL.

1\. Global Definitions

Define all constants (material, load, etc)

Define MATLAB functions (phi.m, interpol.m, etc)

2\. Create Geometry (2D/3D)

3\. Local Variables definition

Define interaction integral equations

For propagating cracks: Assign Global Variable Probe; call MATLAB functions
for crack update,

crack tip and angle (see reference [50])

4\. Select physical model (physics interfaces)

Standard Solid Mechanics (SMstd)

Select material model

Modify stress definitions (Eq. 23)

Add thermo-mechanical coupling using Thermal Expansion subnode

Add hydro-mechanical coupling (effective stress) using External Stress subnode

Enriched Solid Mechanics (SMenr)

Select material model

Modify strain definitions (Eq. 22)

Modify stress definitions (Eq. 23)

Add hydro-mechanical (effective stress) coupling

Apply field variable $\psi$ as a constraint using the domain Pointwise
Constraint option

Standard Darcy’s Law (DLstd)

Select fluid and matrix properties

Modify Darcy’s velocity definitions (Eq. 23)

Add hydro-mechanical coupling (rate of volumetric strain) using Mass Source
subnode

Enriched Darcy’s Law (DLenr)

Select fluid and matrix properties

Modify Darcy’s velocity definitions (Eq. 23)

Modify pressure gradient definitions (Eq. 22)

Add hydro-mechanical coupling (rate of volumetric strain) using Mass Source
subnode

Apply field variable $\psi$ as a constraint using the domain Pointwise
Constraint option

Standard Heat Transfer in Porous Media (HTstd)

Select porous medium and fluid properties

Activate Nonisothermal Flow

Modify Conductive heat flux definitions (Eq. 23)

Modify Convective heat flux definitions (Eq. 23)

Enriched Heat Transfer in Porous Media (HTenr)

Select porous medium and fluid properties

Activate Nonisothermal Flow

Modify Conductive heat flux definitions (Eq. 23)

Modify Convective heat flux definitions (Eq. 23)

Modify heat gradient definitions (Eq. 22)

Apply field variable $\psi$ as a constraint using domain Pointwise Constraint
option

5\. Assign initial and boundary conditions

6\. Domain discretization and mesh generation

7\. Specify Study type

Select Time Dependent (transient) or Stationary (steady-state) analysis

8\. Post-processing and visualization

##  4 Numerical simulations

Several numerical examples are presented in order to illustrate the robustness
and performance of the proposed XFEM implementation procedure. The first
example deals with the validation of the solution strategy in handling
temperature-induced fracture growth in a bi-material composite beam. A thermo-
mechanical analysis is next performed on a cracked domain, for which the
stress intensity factors (SIFs)and temperature distribution are evaluated and
compared to the available analytical solutions in the literature. The
capabilities of the proposed XFEM implementation are further elaborated in the
thermo-mechanical contact analysis conducted on a doubled clamped beam
consisting of a middle crack. Subsequently, an XFEM hydro-mechanical modeling
of impermeable discontinuities in saturated porous media is presented. The
final set of examples is devoted to illustrating the rigor of the proposed
strategy in dealing with challenging thermo-hydro-mechanical simulations of
2D/3D deformable porous rocks in the presence of single/multiple faults.

Unless otherwise specified, linear elastic material behaviour and plane strain
condition are stipulated for all 2D numerical simulations, where bilinear
quadrilateral elements are utilized to discretize the domain. A $20\times 20$
Gauss quadrature is applied over the enriched elements (i.e., enriched physics
interfaces) to achieve satisfactory integration accuracy. In all cases,
$\delta_{s}$ is set to 0.005 to ensure a sufficient number of Gauss
integration points are available at either side of the discontinuities inside
the enriched elements.

###  4.1 Temperature induced fracturing of a notched bi-material beam

The first example studies crack growth in a borosilicate/steel bi-material
beam under thermal loading, which was originally introduced and experimentally
examined in [66]. As shown in Fig. 4, the composite beam consists of a strip
of AISI 304 stainless steel glued to the top edge of a bar of borosilicate
glass, with the respective dimensions of $3.17\text{mm}\times
300\text{mm}\times 6\text{mm}$ and $31.7\text{mm}\times 300\text{mm}\times
6\text{mm}$. An initial notch is placed at a distance of 30 mm from the left
edge of the specimen. The mechanical and thermal material properties of the
domain are presented in Table 1. Although geometrically simple, the specimen
presents a mechanistically complex scenario regarding temperature-induced
crack propagation analysis.

Table 1: Material properties for the bi-material beam problem. |  $E$ (GPa) | $\nu$ |  $\alpha_{T}$ (1/∘C) |  $G_{\rm{f}}$ (N/mm) |  $f_{\rm{t}}$ (MPa)
---|---|---|---|---|---
Borosilicate glass | 64 | 0.2 | $3.25\times 10^{-6}$ | 0.4 | 80
304 Stainless steel | 193 | 0.29 | $17.3\times 10^{-6}$ | - | -

According to the experiment procedure, the specimen is cooled down from the
initial room temperature (i.e., $20^{\circ}$C). The sample undergoes pure
bending induced by the difference in the thermal expansion coefficients of its
constituent materials. Observations show that a crack is initiated at the
notch tip (i.e., within the brittle glass material), which eventually
propagates towards the bi-material interface. Midway through, this trend is
altered by the diversion of the crack propagation trajectory in a direction
parallel to the glass/steel interface. This problem is modeled in COMOSL by
using 14,432 elements under plane stress condition. The average element size
is maintained at $0.5$ mm in the vicinity of the notch zone. For the sake of
simplicity, the material interface is assumed to be fully bonded (no-slip
constraint). The initial crack length and growth increment are set to 11 mm
and $2.5$ mm, respectively. The thermal loading is applied by imposing a
uniform temperature reduction over the entire domain which, in turn, results
in uniform thermal strain ${{\boldsymbol{\varepsilon}}_{\rm{T}}}=\alpha\Delta
T{\bf{I}}$ , where $\Delta T=-1\times i$, with $i=1,..,15$ being the load
increment multiplier. Quasi-static crack growth analysis, using the maximum
hoop stress criterion of $f_{\rm{t}}=80\text{MPa}$, is performed through the
“Auxiliary Sweep” feature in the “Extended Study” section of the “Stationary”
node in COMSOL (for further detail on quasistatic crack growth in COMSOL, see
[50]). Note the mixed-mode SIFs are evaluated following the procedures
explained in A.

Figure 4: The geometry and boundary conditions of the borosilicate/steel
composite beam.

Figs. 5(a) and 5(b) respectively show the simulated and experimental crack
propagation trajectory for the composite specimen. The former is determined
using the field variable $\psi$, that is equal to unity for fully fractured
elements. As can be seen, there is a good agreement between the proposed
numerical solution and the experiment.

(a)

(b)

Figure 5: Crack propagation in the bi-material beam; a) numerical crack
trajectory based on $\psi$, b) experimental crack pattern [66].

###  4.2 Plate with an edge crack subjected to thermal loading

In the second example, a rectangular plate consisting of an adiabatic edge
crack is considered under thermal loading (Fig. 6). This problem was first
introduced by Wilson et al. [67, 68], which presents an analytical study for
an isolated crack in thermo-elastic media. This example was later exploited by
Duflot [69] to develop an XFEM formulation for the modeling of thermo-elastic
fractures. Here, the reference solution is particularly employed to illustrate
the accuracy of the evaluated SIF for the non-isothermal cracked settings
using the proposed XFEM strategy. The problem involves a 2.0 m$\times$0.5 m
plate possessing a half-width sized adiabatic edge crack. The plate is
subjected to prescribed temperatures of $+50^{\rm{o}}{\rm{C}}$ and
$-50^{\rm{o}}{\rm{C}}$ on the right and left edges, respectively. The thermo-
mechanical properties of the plate are presented in Table 2. A structured
uniform mesh of $160\times 40$ elements is supposed, which is used for the
coupled thermo-mechanical analysis with a total duration of $100$ s. Based on
Duflot [69], the normalized stress intensity factor $F_{\rm{I}}$ is evaluated
as

${F_{\rm{I}}}=\frac{{{K_{\rm{I}}}}}{{\frac{E}{{1-\nu}}\alpha{\theta_{0}}\sqrt{\pi
a}}}$ (26)

where $K_{\rm{I}}$ is the mode I SIF, $a$ is the crack length, and
$\theta_{0}$ is the absolute value of the applied temperature on the plate
boundary. The time variation of the computed normalized SIF is compared with
the analytic solution in Fig. 7(a). It is observed that the non-isothermal
solution reaches the steady-state value of 0.491 in 20 s, which indicates less
than $1\%$ relative error with respect to the exact analytical value of 0.495
[68].

Table 2: Material properties for the plate with an edge crack under thermal loading problem. Young’s modulus, $E$ (GPa) | 9
---|---
Poisson’s ratio, $\nu$ | 0.3
Solid density, $\rho_{\rm{s}}$ (kg/m3) | $2\times 10^{3}$
Thermal conductivity, $\lambda_{\rm{s}}$ (W/m ∘C) | $1\times 10^{3}$
Volumetric thermal expansion coefficient, $\beta_{\rm{s}}$ (1/∘C) | $3\times 10^{-7}$

Fig. 7(b) denotes the profile of temperature variation versus the distance
from the crack tip. The simulation results are in an excellent agreement with
the exact solution expressed as $\theta(x)=(2x-w){\theta_{0}}$ [69]. Here,
$\theta(x)$ is the temperature at coordinate $x$ with respect to the crack
tip, and $w$ is the plate’s width. The contours of the temperature, vertical
and horizontal displacements are also depicted in Fig. 8. Evidently, an
intensified contraction of the crack edges can be detected, which is
attributed to the loss of strength along the crack faces.

Figure 6: Geometry and boundary conditions of the rectangular plate with an
edge crack.

(a)

(b)

Figure 7: Thermo-mechanical analysis of edge cracked plate; a) the evolution
of normalized SIF and b) comparison of the COMSOL and analytical temperature
distributions with the distance from crack tip.

(a)

(b)

(c)

Figure 8: The final distribution contours of a) temperature
(${}^{\circ}{\rm{C}}$), b) horizontal displacement (mm), and c) vertical
displacement (mm).

###  4.3 Thermo-mechanical analysis of a cracked double-clamped beam

The following example explores the thermo-mechanical contact analysis in a
double clamped beam under simultaneous thermal and mechanical loading. This
problem was initially introduced by Khoei [52] for purely mechanical contact
analysis. As shown in Fig. 9, the domain consists of a 3 m $\times$ 10 m beam,
containing a full depth vertical crack at the middle. The beam is restrained
along its left end, and is subjected to a prescribed uniform horizontal
displacement of $6$ cm at the opposite end. The left edge is also subjected to
a prescribed temperature of 0 ${}^{\rm{o}}{\rm{C}}$, while the right edge is
retained at a constant temperature of 10 ${}^{\rm{o}}{\rm{C}}$. Furthermore, a
vertical loading of $240$ MN is applied along the top edge of the beam, which
is distributed uniformly over a distance of $0.4$ m about its center. The
thermo-mechanical material properties used in this example are summarized in
Table 3. The domain is discretized using a structured mesh consisting of
$155\times 47$ quadrilateral elements. As a result of the presence of the full
depth vertical crack, the beam is divided into two distinct segments.
Considering the symmetry of the problem, no relative vertical displacement
(i.e., slip) may occur between the crack faces. Therefore, only the normal
contact constraint is applied by means of a penalty approach, with no
frictional resistance. A normal penalty parameter of ${10^{9}}$
${\rm{MN/}}{{\rm{m}}^{\rm{3}}}$ is selected for the mechanical phase. To
impose a fully continuous temperature field across the closure zone (i.e.,
active contact zone), a thermal conductivity coefficient of
${h_{{\rm{cont}}}}={10^{5}}$ W/m∘C is considered within the thermal contact
formulation.

Table 3: Material properties for the double-clamped beam. Young’s modulus, $E$ (GPa) | 10
---|---
Poisson’s ratio, $\nu$ | 0.3
Density, $\rho_{\rm{s}}$ (kg/m3) | $10^{3}$
Thermal conductivity, $\lambda_{\rm{s}}$ (W/m ∘C) | $150$
Specific heat capacity, Cs (J/kg ∘C) | 100
Volumetric thermal expansion coefficient, $\beta_{\rm{s}}$ (1/∘C) | $23.86\times 10^{-6}$
Figure 9: Geometry and boundary conditions of the double-clamped beam under
thermo-mechanical loading.

In this example, the active contact zone is not known a priori; alternatively,
it is determined following an iterative solution procedure which mitigates the
non-positive values of the normal gap (i.e.,
$\left[\kern-1.49994pt\left[{{u_{N}}}\right]\kern-1.49994pt\right]$=0). As
explained in section 3.3.3, this task is executed by using the enriched
portion of the displacement field, which is explicitly derived from SMenr.
Both mechanical and thermal contact contributions are applied via exploiting
the domain Weak Contribution subnodes that are integrated using a $20\times
20$ Gauss integration quadrature. The use of higher order Gauss quadrature
ensures that the integration of the contact force/heat flux along crack faces
is accurately performed by using their corresponding equivalent domain
integrals. Fig. 10 depicts the distribution of temperature, horizontal heat
flux $q_{x}$, the horizontal displacement $u_{x}$, and normal stress in $x$
direction $\sigma_{xx}$, at various stages of the time-dependent analysis.
Fig. 10(b) illustrates a continuous and increasing heat flux across the active
contact zone with a significant flux concentration in the vicinity of the
crack-tip. The results also demonstrate the formation of a continuous
temperature field over the contact zone, particularly as the simulation
reaches the steady-state condition. This can be observed likewise in Fig. 11,
where the temperature profiles are plotted along both crack faces. As
predicted, in opening mode, the Heaviside enrichment function has induced a
jump in the temperature field.

(a)

(b)

(c)

(d)

Figure 10: The contours of a) temperature (${}^{\circ}{\rm{C}}$), b)
horizontal heat flux (${\rm{N/s}}{\rm{.m}}$), c) horizontal displacement (m),
and d) normal stress in $x$ direction ${\sigma_{x}}$ (Pa) at $10^{3}$s,
$5\times 10^{3}$s, $10^{4}$s, $2\times 10^{4}$s and the steady-state condition
(from top to bottom).  Figure 11: Temperature profiles along the both sides of
the discontinuity at the end of simulation.

###  4.4 Hydro-mechanical modeling of impermeable discontinuities in an
embankment dam

This example investigates the applicability of the proposed framework in
dealing with the hydro-mechanical simulation of impermeable discontinuities in
porous media. Sheet piles, frequently used in dam engineering practice, are
structural barriers inserted beneath embankment dams to increase the pore
water pressure dissipation. This, in turn, reduces the uplift pressure applied
to the dam base as well as the occurrence of sand boiling phenomenon (a.k.a.
quicksand) downstream of the dam [70]. In this example, using the XFEM, an
impermeable discontinuity is introduced to emulate the presence of sheet piles
on seepage flow in a porous medium. As illustrated in Fig. 12, a 13.5 m
$\times$ 36 m silty sand formation is considered which is subjected to the
vertical uniform surcharge of $12$ ton/m from the embankment dam as well as a
$100$ kPa pressure from the upstream water reservoir. Two vertical sheet piles
with lengths of $4.6$ m and $2.6$ m are inserted beneath the dam at a lateral
spacing of $15.5$ m apart from each other. The downstream water level is
supposed to be at the ground level (i.e., $p=0$). To investigate the effect of
the mesh size on the solution accuracy, the domain is discretized by three
different mesh configurations using the average element sizes of $0.5$ m,
$0.25$ m and $0.17$ m. The time-dependent analysis is performed for the period
of $110$ s. The mechanical and hydraulic properties used for this study are
provided in Table 4.

Table 4: Material properties for the embankment dam problem containing impermeable discontinuities. Young’s modulus, $E$ (GPa) | 9
---|---
Poisson’s ratio, $\nu$ | 0.4
Solid density, $\rho_{\rm{s}}$ (kg/m3) | $2\times 10^{3}$
Fluid density, $\rho_{\rm{f}}$ (kg/m3) | $10^{3}$
Porosity, n | 0.3
Bulk modulus of solid, $K_{\rm{s}}$ (MPa) | $1\times 10^{14}$
Bulk modulus of fluid, $K_{\rm{f}}$ (MPa) | $2\times 10^{3}$
Fluid viscosity, $\mu_{\rm{f}}$ (Pa.s) | $2\times 10^{-3}$
Permeability, $k_{\rm{f}}$ (m2) | $1\times 10^{-12}$
Figure 12: Embankment dam problem with impermeable discontinuities: geometry
and boundary conditions.

Fig. 13 illustrates the evolution of pressure at three instants of the
analysis i.e., $5$ s, $15$ s and $110$ s. In the latter case (Fig. 13c), the
flow streamlines are also depicted to highlight the effects of the sheet piles
on the deviation of the flow path, particularly along the dam-domain interface
zone. To further elaborate on this issue, the time variation of vertical
displacement and pore water pressure are plotted in Fig. 14 for five critical
points A, B, C, D and E, which are positioned on either side of the sheet pile
as well as the midpoint beneath the dam body (indicated in Fig. 9). The
results obtained by the classic finite element analysis are also shown in Fig.
14 for the sake of comparison. Notably, the FE analysis is performed by
stipulating that the sheet piles are narrow 2D devoid with impermeable
boundaries. According to Fig. 14, an excellent agreement is observed between
the proposed XFEM strategy and the reference FEM solution for both
displacement and pressure fields. Meanwhile, the convergence of the simulation
results corresponding to the three different mesh configurations is deemed
satisfactory.

Figure 13: The pressure distribution contours at 5 s, 15 s and 110 s (values
are in kPa).

(a)

(b)

Figure 14: The evolution of a) vertical displacement, and b) pressure at
various critical points in the porous domain; comparison of proposed XFEM with
the classic FEM strategy

###  4.5 Thermo-hydro-mechanical simulation of faults in 2D/3D porous media

Thermo-hydro-mechanical modeling of fractures and faults manifest in a wide
range of geomechanics engineering applications, including geothermal energy
extraction, oil and gas recovery and deep geological disposal of nuclear
wastes [71]. The following numerical simulations are presented to demonstrate
the remarkable applicability of the proposed framework in dealing with the
2D/3D thermo-hydro-mechanical analysis of impermeable faults in porous media.
The thermo-hydro-mechanical material properties elaborated in this study are
presented in Table 5.

Table 5: Material properties for the thermo-hydro-mechanical analysis of faults in the porous medium. Young’s modulus, $E$ (GPa) | 1.6
---|---
Poisson’s ratio, $\nu$ | 0.33
Solid density, $\rho_{\rm{s}}$ (kg/m3) | $2\times 10^{3}$
Fluid density, $\rho_{\rm{f}}$ (kg/m3) | $10^{3}$
Porosity, n | 0.3
Bulk modulus of solid, $K_{\rm{s}}$ (MPa) | $1\times 10^{14}$
Bulk modulus of fluid, $K_{\rm{f}}$ (MPa) | $2\times 10^{3}$
Fluid viscosity, $\mu_{\rm{f}}$ (Pa.s) | $2\times 10^{3}$
Permeability, $k_{\rm{f}}$ (m2) | $1\times 10^{-12}$
Thermal conductivity of solid, $\lambda_{\rm{s}}$ (W/m ∘C) | 2.88
Thermal conductivity of fluid, $\lambda_{\rm{f}}$ (W/m ∘C) | 0.6
Solid specific heat capacity, Cs (J/kg ∘C) | $1.17\times 10^{3}$
Fluid specific heat capacity, Cf (J/kg ∘C) | $4.2\times 10^{3}$
Volumetric thermal expansion coefficient, $\beta_{\rm{s}}$ (1/∘C) | $6.6\times 10^{-6}$

####  4.5.1 Inclined fault in 2D porous media

The first problem investigates the effects of the presence of an inclined
impermeable fault, with the length of 0.3 m and inclination angle of $\beta$
with respect to the horizon, on the thermo-hydro-mechanical response of a 1 m
$\times$ 1 m square porous block (Fig. 15). The bottom edge of the block is
subjected to a prescribed temperature of 50 ${}^{\rm{o}}{\rm{C}}$ and normal
fluid flux of $10^{-4}$ m/s. The top edge is assumed to be in a drained
condition, where both pressure and temperature vanish (i.e., $p=0$; $T=0$).
The domain is discretized by using a structured mesh of $91\times 91$
elements. Various fault angles i.e., $\beta$=0,
$30^{\rm{o}},45^{\rm{o}},60^{\rm{o}}$ and $90^{\rm{o}}$, are used to study the
significance of strong discontinuity onto both pressure and temperature
fields, over a total period of $10^{4}$ s.

Figure 15: Thermo-hydro-mechanical modeling of an inclined fault: geometry and
boundary conditions.

The distribution contours including the vertical displacement, pressure and
temperature fields for $\beta=45^{\rm{o}}$ are shown in Fig. 16. As observed,
the discontinuity associated with the fault zone is evident in all three
fields, where the Heaviside enrichment is applied. Consistently, Fig 16(c)
shows the flow streamlines, initially aligned vertically, are severely
diverted due to the presence of the impermeable fault, before retrieving their
initial orientation again. In Fig. 17, the evolution of pore pressure at two
neighboring mid-points, positioned at either side of the fault (indicated in
Fig. 15), is plotted for all fault angles. As Fig. 17(a) shows, for smaller
inclination angles $\beta$, the built-up pressure is greater at the bottom
edge of the fault, where the fluid injection occurs. A reversed pattern is
observed for pressure development on the top edge of the discontinuity (Fig.
17(b)). This mechanism is attributed to the fact that the more inclined a
fault is –with respect to the flow streamlines –the more likely it is to
hinder the fluid flow and therefore, the pressure is built up at a higher rate
(Fig. 17(c)).

(a)

(b)

(c)

Figure 16: The distribution contours of a) vertical displacement (mm), b)
pressure (MPa) and c) temperature (${}^{\circ}{\rm{C}}$) at the end of the
analysis.

(a)

(b)

(c)

Figure 17: The pressure evolution at a) bottom and b) top of the
discontinuity, and c) associated pressure difference.

Fig. 18 elaborates the time-history of the temperature on either side of the
fault, as well as the corresponding temperature difference. While the fault
angle is observed to have a negligible effect on the temperature field along
the lower face of the fault (Fig. 18(a)), the upper face undergoes significant
changes in the temperature evolution (Fig. 18(b)). This can be attributed to
the fact that the convective heat transfer to this zone is delayed as a result
of the presence of an impermeable discontinuity that acts as a thermal
barrier. In resemblance to the pressure field, as the fault inclination angle
is reduced the temperature jump due to the presence of the discontinuity is
decreased. As shown in Fig. 18(c), the greatest temperature jump is related to
the case of horizontal fault. In all cases, the temperature jump eventually
decreases to its steady-state condition, when the convective heat diffusion is
fully developed throughout the entire domain.

(a)

(b)

(c)

Figure 18: The temperature evolution at a) bottom and b) top of the
discontinuity, and c) associated temperature difference.

####  4.5.2 Randomly distributed faults in 2D porous media

Natural geological formations are heterogeneous and often contain numerous
impermeable faults and/or inclusions (e.g., hard rocks), which could be
scattered throughout their bulk. [72]. The following example illustrates the
immense capability of the proposed implementation strategy in dealing with the
simulation of multiple faults in porous media (Fig. 19). Similar boundary
conditions and mesh configuration are adopted in a replica of the previous
example. The domain contains 11 equally sized faults with a length of 0.2 m
each, which are studied for a total simulation time of 8000 s. Fig. 20 shows
the pressure and temperature contours at the end of the simulation. The
promising performance of the strong discontinuity enrichment can be easily
recognized across the faults in both fields. The distribution of the heat flux
accompanied by the flow streamlines is depicted in Fig. 21 at the final stage
of the analysis. Evidently, the fluid flow and the heat flux are both
prevented across the contained adiabatic impermeable discontinuities,
resulting in heat flux concentration at the crack tip zones.

Figure 19: Randomly distributed faults in a porous medium: problem geometry
and boundary conditions.

(a)

(b)

Figure 20: The distribution contours of a) pressure (MPa) and b) temperature
(${}^{\circ}{\rm{C}}$) at the end of thermo-hydro-mechanical analysis of
multiple faults problem. Figure 21: The final distribution of heat flux
contour and its associated fluid flow streamlines (values are in
${\rm{N/s}}{\rm{.m}}$).

####  4.5.3 Faults in three dimensional porous media

The final example is presented to show the general nature of the proposed
thermo-hydro-mechanical XFEM framework to the study of three dimensional
problems in geomechanics. Suppose a cubic domain with a side length of 50 m
that encompasses a penny-shaped impermeable discontinuity of diameter
$2{a_{0}}=20$ m located at its center (Fig. 22). All material properties are
the same as the previous examples. The bottom surface is subjected to a
prescribed temperature of 50 ${{}^{\circ}}{\rm{C}}$ and a constant inflow rate
of $10^{-4}$ m/s, while both the pressure and temperature are assumed to
vanish over the top surface (i.e, $p=0$; $T=0$). The remainder surfaces are
considered to be undrained, with no fluid flow/heat flux. All faces of the
domain are constrained in their corresponding normal directions, except the
top surface, so as to emulate the in-situ boundary conditions. The porous
domain is discretized using 2,312 brick elements, clustered in the vicinity of
the internal discontinuity with an average element size of 1.5 m, in
conjunction with 45,415 tetrahedral elements, elsewhere, that is simulated for
the total duration of $3\times{10^{5}}$ s. Fig. 23 illustrates the contours of
vertical displacement $u_{z}$ as well as pressure and temperature fields, over
two cross-shaped planes perpendicular to the discontinuity. Clearly, the
discontinuity induced by the penny-shaped inclusion can be observed in all
three distribution contours. This is further elaborated by noting the flow and
heat flux streamlines that are depicted in Figs. 23(b) and 23(c), where a
diversion from the far-field vertical alignment can be observed adjacent to
the discontinuity region. The promising results presented here showcase the
flexibility of the implemented technique in dealing with intricate scenarios
in the 2D/3D thermo-hydro-mechanical analysis of porous media with
single/multiple faults.

Figure 22: Three dimensional thermo-hydro-mechanical simulation of a fault in
porous media; problem geometry and boundary conditions.

(a)

(b)

(c)

Figure 23: 3D thermo-hydro-mechanical simulation results on two perpendicular
planes; a) distribution of vertical displacement $u_{z}$ (mm), b) pressure
distribution (MPa) together with fluid flow streamlines, and c) temperature
distribution (${}^{\circ}{\rm{C}}$) with heat flux streamlines.

##  5 Conclusions

In this study, an XFEM implementation for modeling thermo-hydro-mechanical
problems in COMSOL Multiphysics commercial software is presented. The
framework is applied for multi-field fracture analysis in 2D and 3D settings.
The enrichment strategy is implemented by exploiting specific weak forms of
the governing equations, which are consistent with the structure of the
software. In order to take account of the standard and enriched solution
fields, two sets of distinct ”Solid Mechanics”, ”Darcy’s Law” and ”Heat
Transfer in Porous Media” physics interfaces associated with the mechanical
deformation of the bulk, fluid flow within the pores and heat transfer through
the matrix are adopted, respectively. COMSOL’s built-in features and external
MATLAB functions are employed in the modeling procedures of preprocessing and
level-set updating, various coupling effects between the physics and
postprocessing.

Several numerical examples are presented to investigate the robustness and
accuracy of the proposed implementation in geomechanical applications. In the
first example, the crack growth in a bi-material composite beam subjected to
thermal loading is studied and validated against the experimental
observations. The second example investigates the thermo-mechanical analysis
of a cracked plate, in which a comprehensive verification against an exact
analytical solution is carried out. The coupled thermo-mechanical analysis is
further elaborated in the third example to explore the thermo-mechanical
contact within a doubled clamped beam. The next example is devoted to the
hydro-mechanical modeling of discontinuities in porous media. Finally, a set
of examples investigating the thermo-hydro-mechanical XFEM modeling of
single/multiple discontinuities in 2D/3D deformable porous media are presented
in the last example. In this fashion, the robustness and superiority of the
proposed XFEM framework in dealing with intricate real-world multi-field
problems is demonstrated. The simplicity of the proposed algorithms, the
feasibility of inclusion of further physics, and the ease of incorporation of
alternative enrichment strategies facilitate the use of state-of-the-art XFEM
simulation of multiphysics, to the extent that the need for in-house codes
could be alleviated. Showcasing the practicality of the inclusion of the
enumerated scenarios in COMSOL Multiphysics is left to future studies.

## Appendix A

For a non-isothermal cracked domain, the stress intensity factor (SIF) can be
obtained based on the J-integral concept as [73, 69]

$J=\int\limits_{{\Gamma_{\rm{J}}}}{\left[{{w_{{\rm{mec}}}}\cdot{n_{x^{\prime}}}-({\boldsymbol{\sigma}}\cdot{\nabla_{x^{\prime}}}{\bf{u}})\cdot{\bf{n}}}\right]}d\Gamma$
(27)

where
${w_{{\rm{mec}}}}=\frac{1}{2}{\boldsymbol{\sigma}}:({{\boldsymbol{\varepsilon}}}-{{\boldsymbol{\varepsilon}}_{{\rm{T}}}}){\bf{I}}$
is the mechanical strain energy density. In addition, ${\bf{n}}$ and
${{n_{x^{\prime}}}}$ are the outward normal vector and its horizontal
component associated with the integral path ${{\Gamma_{\rm{J}}}}$, that
encloses the crack tip (in the local crack coordinates $x’-y’$). For mixed-
mode loading in two dimensions, the relation between the stress intensity
factors and J-integral can be expressed as,

$J=\frac{K_{\rm{I}}^{2}+K_{\rm{II}}^{2}}{{E}^{\prime}}$ (28)

where ${E}^{\prime}$ is defined as $E/(1-\upsilon^{2})$ and $E$ for plane
strain and plane stress problems, respectively.

The stress intensity factors can be calculated by means of the interaction
integral method with satisfactory accuracy [73]. Accordingly, an auxiliary
displacement field is employed in addition to the actual displacement field of
the problem, denoted by states (2) and (1), respectively, such that

$I^{(1+2)}=\int_{\Gamma_{\rm{J}}}\left[W^{(1,2)}\cdot
n_{{x}^{\prime}}-(\boldsymbol{\sigma}^{{(1)}}\nabla_{{x}^{\prime}}\mathbf{u}^{(2)}+\boldsymbol{\sigma}^{{(2)}}\nabla_{{x}^{\prime}}\mathbf{u}^{(1)})\right]d\Gamma$
(29)

where
$W^{(1,2)}=\boldsymbol{\sigma}^{(1)}\cdot\boldsymbol{\varepsilon}^{(2)}=\boldsymbol{\sigma}^{(2)}\cdot\boldsymbol{\varepsilon}^{(1)}$
denotes the interaction strain energy. Combining Eqs. 27, 28 and 29, the
interaction integral is related to the SIFs corresponding to the actual and
auxiliary states as

$I^{(1+2)}=\frac{2}{{E}^{\prime}}(K_{\text{I}}^{(1)}K_{\text{I}}^{(2)}+K_{\text{II}}^{(1)}K_{\text{II}}^{(2)})$
(30)

The stress intensity factors of a mixed-mode crack domain can be therefore
calculated by taking advantage of the pure mode I (i.e., $K_{\rm{I}}^{(2)}=1$,
$K_{\rm{II}}^{(2)}=0$) or pure mode II (i.e., $K_{\rm{I}}^{(2)}=0$,
$K_{\rm{II}}^{(2)}=1$) asymptotic solutions. Moreover, the alternative domain
form of Eq. 29 can also be employed by the use of special weighting functions
[33].

In this study, the internal variables and built-in mathematical functions of
COMSOL are employed to carry out the stress intensity factor calculations. To
this end, a list of all required variables is defined, which includes: normal
vector components, actual and auxiliary stress and strain fields and associate
derivatives, the displacement field, the crack tip coordinates, the
interaction strain energy density, etc. In this study, the numerical
integration of Eq. 29 is carried out by exploiting the built-in circular path
integration operator circint in COMSOL. Alternatively, the equivalent domain
form of interaction integral can be applied using the diskint operator for a
circular area encompassing the crack tip.

## References

  * Khalili and Loret [2001] N. Khalili, B. Loret,  An elasto-plastic model for non-isothermal analysis of flow and deformation in unsaturated porous media: formulation,  International Journal of Solids and Structures 38 (2001) 8305–8330.
  * Khalili and Selvadurai [2003] N. Khalili, A. Selvadurai,  A fully coupled constitutive model for thermo-hydro-mechanical analysis in elastic media with double porosity,  Geophysical Research Letters 30 (2003).
  * Gelet et al. [2012] R. Gelet, B. Loret, N. Khalili,  A thermo-hydro-mechanical coupled model in local thermal non-equilibrium for fractured hdr reservoir with double porosity,  Journal of Geophysical Research: Solid Earth 117 (2012).
  * Bower and Zyvoloski [1997] K. Bower, G. Zyvoloski,  A numerical model for thermo-hydro-mechanical coupling in fractured rock,  International Journal of Rock Mechanics and Mining Sciences 34 (1997) 1201–1211.
  * Watanabe et al. [2010] N. Watanabe, W. Wang, C. I. McDermott, T. Taniguchi, O. Kolditz,  Uncertainty analysis of thermo-hydro-mechanical coupled processes in heterogeneous porous media,  Computational Mechanics 45 (2010) 263–280.
  * Sarris and Papanastasiou [2012] E. Sarris, P. Papanastasiou,  Modeling of hydraulic fracturing in a poroelastic cohesive formation,  International Journal of Geomechanics 12 (2012) 160–167.
  * Zhang et al. [2020] Y. Zhang, G.-F. Zhao, Q. Li,  Acoustic emission uncovers thermal damage evolution of rock,  International Journal of Rock Mechanics and Mining Sciences 132 (2020) 104388\.
  * Zhou et al. [2011] Y. Zhou, K.-w. Xia, X. Li, H. Li, G. Ma, J. Zhao, Z. Zhou, F. Dai,  Suggested methods for determining the dynamic strength parameters and mode-i fracture toughness of rock materials,  in: The ISRM Suggested Methods for Rock Characterization, Testing and Monitoring: 2007-2014, Springer, 2011, pp. 35–44.
  * Terzaghi et al. [1996] K. Terzaghi, R. B. Peck, G. Mesri,  Soil mechanics,  New York: John Wiley & Sons (1996).
  * Zienkiewicz et al. [1999] O. C. Zienkiewicz, A. Chan, M. Pastor, B. Schrefler, T. Shiomi, Computational geomechanics, volume 613, Citeseer, 1999\.
  * Coussy [2004] O. Coussy, Poromechanics, John Wiley & Sons, 2004.
  * Sarris and Papanastasiou [2011] E. Sarris, P. Papanastasiou,  The influence of the cohesive process zone in hydraulic fracturing modelling,  International Journal of Fracture 167 (2011) 33–45.
  * Tamizdoust and Ghasemi-Fare [2020] M. M. Tamizdoust, O. Ghasemi-Fare,  A fully coupled thermo-poro-mechanical finite element analysis to predict the thermal pressurization and thermally induced pore fluid flow in soil media,  Computers and Geotechnics 117 (2020) 103250.
  * Iranmanesh and Pak [2018] M. A. Iranmanesh, A. Pak,  Extrinsically enriched element free galerkin method for heat and fluid flow in deformable porous media involving weak and strong discontinuities,  Computers and Geotechnics 103 (2018) 179–192.
  * Gawin et al. [2006] D. Gawin, F. Pesavento, B. A. Schrefler,  Hygro-thermo-chemo-mechanical modelling of concrete at early ages and beyond. part i: hydration and hygro-thermal phenomena,  International journal for numerical methods in engineering 67 (2006) 299–331.
  * Belytschko and Black [1999] T. Belytschko, T. Black,  Elastic crack growth in finite elements with minimal remeshing,  International Journal for Numerical Methods in Engineering 45 (1999) 601–620.
  * Moës et al. [1999] N. Moës, J. Dolbow, T. Belytschko,  A finite element method for crack growth without remeshing,  International journal for numerical methods in engineering 46 (1999) 131–150.
  * Dolbow et al. [2000] J. Dolbow, N. Moës, T. Belytschko,  Discontinuous enrichment in finite elements with a partition of unity method,  Finite elements in analysis and design 36 (2000) 235–260.
  * Sukumar et al. [2000] N. Sukumar, N. Moës, B. Moran, T. Belytschko,  Extended finite element method for three-dimensional crack modeling,  International journal for numerical methods in engineering 48 (2000) 1549–1570.
  * Réthoré et al. [2007] J. Réthoré, R. d. Borst, M.-A. Abellan,  A two-scale approach for fluid flow in fractured porous media,  International Journal for Numerical Methods in Engineering 71 (2007) 780–800.
  * Réthoré et al. [2008] J. Réthoré, R. De Borst, M.-A. Abellan,  A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks,  Computational Mechanics 42 (2008) 227–238.
  * Khoei and Mohammadnejad [2011] A. Khoei, T. Mohammadnejad,  Numerical modeling of multiphase fluid flow in deforming porous media: A comparison between two-and three-phase models for seismic analysis of earth and rockfill dams,  Computers and Geotechnics 38 (2011) 142–166.
  * Khoei et al. [2014] A. R. Khoei, M. Vahab, E. Haghighat, S. Moallemi,  A mesh-independent finite element formulation for modeling crack growth in saturated porous media based on an enriched-fem technique,  International Journal of Fracture 188 (2014) 79–108.
  * Salimzadeh and Khalili [2015] S. Salimzadeh, N. Khalili,  A three-phase xfem model for hydraulic fracturing with cohesive crack propagation,  Computers and Geotechnics 69 (2015) 82–92.
  * Haddad and Sepehrnoori [2016] M. Haddad, K. Sepehrnoori,  Xfem-based czm for the simulation of 3d multiple-cluster hydraulic fracturing in quasi-brittle shale formations,  Rock Mechanics and Rock Engineering 49 (2016) 4731–4748.
  * Khoei et al. [2018] A. Khoei, M. Vahab, M. Hirmand,  An enriched–fem technique for numerical simulation of interacting discontinuities in naturally fractured porous media,  Computer Methods in Applied Mechanics and Engineering 331 (2018) 197–231.
  * Vahab et al. [2019] M. Vahab, A. Khoei, N. Khalili,  An x-fem technique in modeling hydro-fracture interaction with naturally-cemented faults,  Engineering Fracture Mechanics 212 (2019) 269–290.
  * Cruz et al. [2019] F. Cruz, D. Roehl, E. do Amaral Vargas Jr,  An xfem implementation in abaqus to model intersections between fractures in porous rocks,  Computers and Geotechnics 112 (2019) 135–146.
  * Parchei-Esfahani et al. [2020] M. Parchei-Esfahani, B. Gee, R. Gracie,  Dynamic hydraulic stimulation and fracturing from a wellbore using pressure pulsing,  Engineering Fracture Mechanics 235 (2020) 107152.
  * Mikaeili and Schrefler [2018] E. Mikaeili, B. Schrefler,  Xfem, strong discontinuities and second-order work in shear band modeling of saturated porous media,  Acta Geotechnica 13 (2018) 1249–1264.
  * Jafari et al. [2021] A. Jafari, M. Vahab, N. Khalili,  Fully coupled xfem formulation for hydraulic fracturing simulation based on a generalized fluid leak-off model,  Computer Methods in Applied Mechanics and Engineering 373 (2021) 113447\.
  * Jin and Arson [2020] W. Jin, C. Arson,  Fluid-driven transition from damage to fracture in anisotropic porous media: a multi-scale xfem approach,  Acta Geotechnica 15 (2020) 113–144.
  * Khoei et al. [2012] A. Khoei, S. Moallemi, E. Haghighat,  Thermo-hydro-mechanical modeling of impermeable discontinuity in saturated porous media with x-fem technique,  Engineering Fracture Mechanics 96 (2012) 701–723.
  * Li et al. [2016] S. Li, X. Li, D. Zhang,  A fully coupled thermo-hydro-mechanical, three-dimensional model for hydraulic stimulation treatments,  Journal of Natural Gas Science and Engineering 34 (2016) 64–84.
  * Khoei and Bahmani [2018] A. Khoei, B. Bahmani,  Application of an enriched fem technique in thermo-mechanical contact problems,  Computational Mechanics 62 (2018) 1127–1154.
  * Nguyen-Tuan et al. [2019] L. Nguyen-Tuan, S. S. Nanthakumar, T. Lahmer,  Identification of multiple flaws in dams using inverse analysis based on hydro-mechanical xfem and level sets,  Computers and Geotechnics 110 (2019) 211–221.
  * Zeng et al. [2020] Q. Zeng, J. Yao, J. Shao,  An extended finite element solution for hydraulic fracturing with thermo-hydro-elastic–plastic coupling,  Computer Methods in Applied Mechanics and Engineering 364 (2020) 112967\.
  * Khoei and Mortazavi [2020] A. Khoei, S. Mortazavi,  Thermo-hydro-mechanical modeling of fracturing porous media with two-phase fluid flow using x-fem technique,  International Journal for Numerical and Analytical Methods in Geomechanics 44 (2020) 2430–2472.
  * Flemisch et al. [2018] B. Flemisch, I. Berre, W. Boon, A. Fumagalli, N. Schwenck, A. Scotti, I. Stefansson, A. Tatomir,  Benchmarks for single-phase flow in fractured porous media,  Advances in Water Resources 111 (2018) 239–258.
  * Keilegavlen et al. [2021] E. Keilegavlen, R. Berge, A. Fumagalli, M. Starnoni, I. Stefansson, J. Varela, I. Berre,  Porepy: An open-source software for simulation of multiphysics processes in fractured porous media,  Computational Geosciences 25 (2021) 243–265.
  * Giner et al. [2009] E. Giner, N. Sukumar, J. Tarancón, F. Fuenmayor,  An abaqus implementation of the extended finite element method,  Engineering fracture mechanics 76 (2009) 347–368.
  * Xu and Yuan [2009] Y. Xu, H. Yuan,  On damage accumulations in the cyclic cohesive zone model for xfem analysis of mixed-mode fatigue crack growth,  Computational Materials Science 46 (2009) 579–585.
  * Ooi et al. [2018] G. T. C. Ooi, S. Roy, S. Sundararajan,  Investigating the effect of retained austenite and residual stress on rolling contact fatigue of carburized steel with xfem and experimental approaches,  Materials Science and Engineering: A 732 (2018) 311–319.
  * Pepper and Heinrich [2017] D. W. Pepper, J. C. Heinrich, The finite element method: basic concepts and applications with MATLAB, MAPLE, and COMSOL, CRC press, 2017\.
  * Hokr et al. [2018] M. Hokr, P. Rálek, K. Sosna, J. Franěk,  In situ measurements of thermally loaded rock and evaluation of an experiment with a 3d numerical model,  International Journal of Rock Mechanics and Minings Sciences (2018).
  * Michalec et al. [2021] Z. Michalec, R. Blaheta, M. Hasal, T. Ligurskỳ,  Fully coupled thermo-hydro-mechanical model with oversaturation and its validation to experimental data from febex experiment,  International Journal of Rock Mechanics and Mining Sciences 139 (2021) 104567\.
  * Segura and Carol [2008] J. M. Segura, I. Carol,  Coupled hm analysis using zero-thickness interface elements with double nodes. part i: Theoretical model,  International journal for numerical and analytical methods in geomechanics 32 (2008) 2083–2101.
  * Zhou et al. [2018] S. Zhou, T. Rabczuk, X. Zhuang,  Phase field modeling of quasi-static and dynamic crack propagation: Comsol implementation and case studies,  Advances in Engineering Software 122 (2018) 31–49.
  * Zhou et al. [2019] S. Zhou, X. Zhuang, T. Rabczuk,  Phase-field modeling of fluid-driven dynamic cracking in porous media,  Computer Methods in Applied Mechanics and Engineering 350 (2019) 169–198.
  * Jafari et al. [2021] A. Jafari, P. Broumand, M. Vahab, N. Khalili,  An extended finite element method implementation in comsol multiphysics: Solid mechanics,  arXiv preprint arXiv:2109.03153 (2021).
  * Babuska and Melenk [1997] I. Babuska, J. M. Melenk,  The partition of unity method,  International Journal for Numerical Methods in Engineering 40 (1997) 727–758.
  * Khoei [2014] A. R. Khoei, Extended finite element method: theory and applications, John Wiley & Sons, 2014\.
  * Borja [2008] R. I. Borja,  Assumed enhanced strain and the extended finite element methods: A unification of concepts,  Computer Methods in Applied Mechanics and Engineering 197 (2008) 2789–2803.
  * Biot [1956] M. A. Biot,  General solutions of the equations of elasticity and consolidation for a porous material,  Journal of Applied Mechanics 23 (1956) 91–96.
  * Khalili et al. [2010] N. Khalili, A. Uchaipichat, A. Javadi,  Skeletal thermal expansion coefficient and thermo-hydro-mechanical constitutive relations for saturated homogeneous porous media,  Mechanics of materials 42 (2010) 593–598.
  * Salimzadeh et al. [2018] S. Salimzadeh, A. Paluszny, H. M. Nick, R. W. Zimmerman,  A three-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal systems,  Geothermics 71 (2018) 212–224.
  * Lewis et al. [1998] R. W. Lewis, R. W. Lewis, B. Schrefler, The finite element method in the static and dynamic deformation and consolidation of porous media, John Wiley & Sons, 1998.
  * Madhusudana and Fletcher [1986] C. Madhusudana, L. S. Fletcher,  Contact heat transfer-the last decade,  AIAA journal 24 (1986) 510–523.
  * Liu and Borja [2008] F. Liu, R. I. Borja,  A contact algorithm for frictional crack propagation with the extended finite element method,  International Journal for Numerical methods in engineering 76 (2008) 1489–1512.
  * Mul [2019] COMSOL Multiphysics 5.5, Introduction to COMSOL Multiphysics (2019).
  * mul [2019] COMSOL Multiphysics 5.5, Structural Mechanics Module User’s Guide (2019).
  * mul [2018a] COMSOL Multiphysics 5.4, Subsurface Flow Module User’s Guide (2018a).
  * mul [2018b] COMSOL Multiphysics 5.4, Heat Transfer Module User’s Guide (2018b).
  * Mohammadi [2008] S. Mohammadi, Extended finite element method: for fracture analysis of structures, John Wiley & Sons, 2008.
  * mul [2019] COMSOL Multiphysics 5.5, COMSOL Multiphysics Reference Manual (2019).
  * Grutzik and Reedy [2018] S. J. Grutzik, E. D. Reedy,  Crack path selection in thermally loaded borosilicate/steel bibeam specimen,  Experimental Mechanics 58 (2018) 1–10.
  * Wilson and Yu [1979] W. Wilson, I.-W. Yu,  The use of the j-integral in thermal stress crack problems,  International Journal of Fracture 15 (1979) 377–387.
  * Shih et al. [1986] C. Shih, B. Moran, T. Nakamura,  Energy release rate along a three-dimensional crack front in a thermally stressed body,  International Journal of fracture 30 (1986) 79–102.
  * Duflot [2008] M. Duflot,  The extended finite element method in thermoelastic fracture mechanics,  International Journal for Numerical Methods in Engineering 74 (2008) 827–847.
  * Das [2019] B. M. Das, Advanced soil mechanics, CRC press, 2019.
  * Najari and Selvadurai [2014] M. Najari, A. Selvadurai,  Thermo-hydro-mechanical response of granite to temperature changes,  Environmental earth sciences 72 (2014) 189–198.
  * Vahab et al. [2021] M. Vahab, M. Hirmand, A. Jafari, N. Khalili,  Numerical analysis of multiple hydro-fracture growth in layered media based on a non-differentiable energy minimization approach,  Engineering Fracture Mechanics 241 (2021) 107361.
  * Anderson [2017] T. L. Anderson, Fracture mechanics: fundamentals and applications, CRC press, 2017\.