A Phase-Field Description for Pressurized and Non-Isothermal Propagating Fractures

TL;DR

This paper extends phase-field modeling of pressurized fractures to include non-isothermal conditions, incorporating thermodynamical principles, and demonstrates efficient numerical algorithms with applications in 2D and 3D.

Contribution

It introduces a thermodynamically consistent phase-field model for pressurized, non-isothermal fractures and develops robust numerical methods for its simulation.

Findings

Effective semi-smooth Newton solution approach

Successful verification with manufactured solutions

Parallel implementation reduces computational cost

Abstract

In this work, we extend a phase-field approach for pressurized fractures to non-isothermal settings. Specifically, the pressure and the temperature are given quantities and the emphasis is on the correct modeling of the interface laws between a porous medium and the fracture. The resulting model is augmented with thermodynamical arguments and then analyzed from a mechanical perspective. The numerical solution is based on a robust semi-smooth Newton approach in which the linear equation systems are solved with a generalized minimal residual method and algebraic multigrid preconditioning. The proposed modeling and algorithmic developments are substantiated with different examples in two- and three dimensions. We notice that for some of these configurations manufactured solutions can be constructed, allowing for a careful verification of our implementation. Furthermore, crack-oriented predictor-corrector adaptivity and a parallel implementation are used to keep the computational cost reasonable. Snapshots of iteration numbers show an excellent performance of the nonlinear and linear solution algorithms. Lastly, for some tests, a computational analysis of the effects of strain-energy splitting is performed, which has not been undertaken to date for similar phase-field settings involving pressure, fluids or non-isothermal effects.

Tables

Table 1. Table 1: Combined Newton loop at time step t n superscript 𝑡 𝑛 t^{n} .

Input: loading data $(𝒖,\phi ,p_t,p_0,\Theta ,{\Theta }_0,t)$ on ${\Gamma }_D,{\Gamma }_N\subset \partial \Omega$;

solution ${𝑼}_{n-1}:=({𝒖}_{n-1},{\phi }_{n-1})$ from time step $n-1$.

Combined Newton / primal-dual active set iteration $k\ge 1$:

•  Assemble residual $𝑹\left({𝑼}_h^k\right)$,

•  Compute active set ${𝒜}_k=\{i\mid {\left({𝑩}^{-1}\right)}_{ii}{\left({𝑹}_k\right)}_i+c{\left(\delta {𝑼}_h^k\right)}_i>0\}$,

•  Assemble matrix $𝑮={\nabla }^2{E}_{\epsilon }\left({𝑼}_h^k\right)$ and right-hand side $𝑭=-\nabla E_{\epsilon }\left({𝑼}_h^k\right)$,

•  Eliminate rows and columns in ${𝒜}_k$ from $𝑮$ and $𝑭$ to obtain $\widetilde{𝑮}$ and $\widetilde{𝑭}$,

•  Solve linear system $\widetilde{𝑮}\delta {𝑼}_k=\widetilde{𝑭}$, i.e, find $\delta {𝑼}_h^k\in V_h\times W_h$ with,

${\nabla }^2{ℰ}_{\epsilon }\left({𝑼}_h^k\right)(\delta {𝑼}_h^k,𝚿)=-\nabla {ℰ}_{\epsilon }\left({𝑼}_h^k\right)\left(𝚿\right),\forall 𝚿\in V_h\times W_h,$

where ${\nabla }^2{ℰ}_{\epsilon }$ and $\nabla {ℰ}_{\epsilon }$ are defined in Section 4.4 and 4.5.

•  Find a step size using back-tracking line search algorithm to get

${𝑼}_h^{k+1}={𝑼}_h^k+\omega \delta {𝑼}_h^k$ with .

•  if fulfilled, such that

)

set $({𝒖}^k,{\phi }^k)=:({𝒖}_t,{\phi }_t)$ and stop;

•  else $k+1\to k$.

Output: solution $({𝒖}_n,{\phi }_n)=:{𝑼}_n$.

Table 2. Table 2: Case d. Average number of GMRES iterations for one Newton cycle for four different times steps.

Time step (day)	no split	vol./ dev. split
100	2	2
114	6	15
121	6	16
136	16	17

Table 3. Table 3: Case e. Average number of GMRES iterations per Newton cycle for different times steps.

Time step (min)	no split	vol./dev. split
240	2	2
252	10	14
259	4	17
267	13	24

Equations

$$\mathcal{E}(\mathbf{u})$$

$$\mathbf{\sigma }={\mathbf{\sigma }}_0+{\mathbf{\sigma }}_{\epsilon }-{\alpha }_B(p-p_0)\mathbf{I}-3{\alpha }_{\Theta }{K}_d(\Theta -{\Theta }_0)\mathbf{I}.$$

$${\mathbf{\sigma }}_{\mathbf{\epsilon }}=2\mu \mathbf{\epsilon }+\lambda tr(\mathbf{\epsilon })\mathbf{I},\mathbf{\epsilon }(\mathbf{u})=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T),$$

$$\mathcal{E}(\mathbf{u})$$

$$=\frac{1}{2}({\mathbf{\sigma }}_{\mathbf{\epsilon }},\mathbf{\epsilon }(\mathbf{u}){)}_{\Omega }-⟨\mathbf{\tau },\mathbf{u}{⟩}_{{\partial }_N\Omega }-({\alpha }_B(p-p_0),\nabla \cdot \mathbf{u})-(3{\alpha }_{\Theta }{K}_d(\Theta -{\Theta }_0),\nabla \cdot \mathbf{u}).$$

$$G_c{\mathcal{H}}^{d-1}(\mathcal{C}),$$

$$\mathcal{E}(\mathbf{u},\mathcal{C})$$

$$+G_c{\mathcal{H}}^{d-1}(\mathcal{C}).$$

$${\partial }_t\phi \le 0.$$

$$V:=H_0^1(\Omega {)}^d,W:=H^1(\Omega ),W_{in}:=\{\phi \in H^1(\Omega )|0\le \phi \le {\phi }^{old}\}.$$

$$0=t_0<t_1<\dots <t_n<\dots <t_N=T.$$

$${\mathcal{E}}_{\epsilon }(\mathbf{u},\phi )$$

$$-\underset{\text{pressure term}}{\underbrace{({\alpha }_B(p-p_0){\phi }^2,\nabla \cdot \mathbf{u})}}-\underset{\text{thermal term}}{\underbrace{(3{\alpha }_{\Theta }{K}_d(\Theta -{\Theta }_0){\phi }^2,\nabla \cdot \mathbf{u})}}$$

$$\displaystyle\quad+\underbrace{G_{c}\bigg{(}\frac{1}{2\varepsilon}\|1-\varphi\|^{2}+\frac{\varepsilon}{2}\|\nabla\varphi\|^{2}\bigg{)}}_{\text{fracture term}}.$$

$$(((1-\kappa ){\phi }^2+\kappa ){\mathbf{\sigma }}_{\mathbf{\epsilon }},\mathbf{\epsilon }(\mathbf{w}))-⟨\mathbf{\tau },\mathbf{w}⟩$$

$$-({\alpha }_B(p-p_0){\phi }^2,\nabla \cdot \mathbf{w})-(3{\alpha }_{\Theta }{K}_d(\Theta -{\Theta }_0){\phi }^2,\nabla \cdot \mathbf{w})=0\forall \mathbf{w}\in V,$$

$$(1-\kappa )(\phi {\mathbf{\sigma }}_{\mathbf{\epsilon }}(\mathbf{u}):\mathbf{\epsilon }(\mathbf{u}),\psi -\phi )$$

$$-2({\alpha }_B-1)((p-p_0)\nabla \cdot \mathbf{u}\phi ,\psi -\phi )-2(3{\alpha }_{\Theta }{K}_d)((\Theta -{\Theta }_0)\nabla \cdot \mathbf{u}\phi ,\psi -\phi )$$

$$+G_c(\frac{1}{\epsilon }(\phi -1,\psi -\phi )+\epsilon (\nabla \phi ,\nabla (\psi -\phi )))\ge 0\forall \psi \in W\cap L^{\infty }.$$

$$-(p-p_0)\mathbf{I}.$$

$$-C_{\Theta }(\Theta -{\Theta }_0)\mathbf{I},\text{with }{C}_{\Theta }:=C_{\Theta }(x,t).$$

$${\partial }_t\Theta -\nabla \cdot ({\kappa }_{\Theta }\nabla \Theta )=0,$$

$$C_{\Theta }=A_{\Theta }\left(\frac{{\lambda }_{\Theta }}{2{\lambda }_{\Theta }+1}\right),A_{\Theta }=\frac{E_Y\beta }{1-{\nu }_s},$$

$${\lambda }_{\Theta }={\sinh }^{-1}(\frac{{\gamma }_T}{0.5l_0}\sqrt{\pi {\kappa }_{\Theta }t}),$$

$$p_f=p_R\text{and}{\Theta }_f={\Theta }_R\text{on }\mathcal{C},$$

$$\mathbf{\sigma }\mathbf{n}$$

$$=-\underset{\text{Fracture stress}}{\underbrace{[(p-p_0)-C_{\Theta }(\Theta -{\Theta }_0)]}}\mathbf{n}.$$

$$⟨\mathbf{\sigma }\mathbf{n},\mathbf{u}{⟩}_{{\partial }_N\Omega }=⟨\mathbf{\tau },\mathbf{u}{⟩}_{{\Gamma }_N}+⟨\mathbf{\sigma }\mathbf{n},\mathbf{u}{⟩}_{\partial \mathcal{C}}.$$

$$⟨\mathbf{\sigma }\mathbf{n},\mathbf{u}{⟩}_{\partial \mathcal{C}}$$

$$=-(\nabla \cdot (p-p_0)\mathbf{u})+⟨(p-p_0)\mathbf{n},\mathbf{u}{⟩}_{{\Gamma }_N}-(C_{\Theta }\nabla \cdot (\Theta -{\Theta }_0)\mathbf{u})$$

$$+⟨C_{\Theta }(\Theta -{\Theta }_0)\mathbf{n},\mathbf{u}{⟩}_{{\Gamma }_N}$$

$$=-((p-p_0),\nabla \cdot \mathbf{u})-(\nabla (p-p_0),\mathbf{u})+⟨(p-p_0)\mathbf{n},\mathbf{u}{⟩}_{{\Gamma }_N}$$

$$-(C_{\Theta }(\Theta -{\Theta }_0),\nabla \cdot \mathbf{u})-(C_{\Theta }\nabla (\Theta -{\Theta }_0),\mathbf{u})+⟨C_{\Theta }(\Theta -{\Theta }_0)\mathbf{n},\mathbf{u}{⟩}_{{\Gamma }_N}.$$

$$({\alpha }_B(p-p_0),\nabla \cdot \mathbf{u})\text{and}(3{\alpha }_{\Theta }{K}_d(\Theta -{\Theta }_0),\nabla \cdot \mathbf{u}).$$

$$-⟨\mathbf{\sigma }\mathbf{n},\mathbf{u}{⟩}_{{\partial }_N\Omega }-({\alpha }_B(p-p_0),\nabla \cdot \mathbf{u})-(3{\alpha }_{\Theta }{K}_d(\Theta -{\Theta }_0),\nabla \cdot \mathbf{u})$$

$$=-⟨\mathbf{\tau },\mathbf{u}{⟩}_{{\Gamma }_N}-⟨\mathbf{\sigma }n,\mathbf{u}{⟩}_{\mathcal{C}}-({\alpha }_B(p-p_0),\nabla \cdot \mathbf{u})-(3{\alpha }_{\Theta }{K}_d(\Theta -{\Theta }_0),\nabla \cdot \mathbf{u})$$

$$=-⟨\mathbf{\tau },\mathbf{u}{⟩}_{{\Gamma }_N}-({\alpha }_B(p-p_0),\nabla \cdot \mathbf{u})+((p-p_0),\nabla \cdot \mathbf{u})+(\nabla (p-p_0),\mathbf{u})$$

Full text

1 Introduction

In subsurface hydraulic fractures, it is known that cold injections (e.g., fluids or $CO2$) may open and advance fractures [23, 46, 21, 13, 45, 16, 41, 40, 38]. These effects play a role in oil/gas production, $CO2$ storage, and geothermal energy production. For instance, [21] investigated the effects of cold-water injection on the growth and decay of a pre-existing crack when pressure and temperature are taken into account. They developed an analytical formula to solve their model. Another study incorporating thermal effects was conducted in [13]. Therein, the fracture propagation is based on the displacement discontinuity method. For similar computations the displacement discontinuity method was used in [22]. In [46], the interaction of temperature and fluid stresses with respect to wellbores have been investigated. The authors show that warm fluids avoid (or restrict) fracture development. Recently, [41, 40] studied $CO2$ injections and their influence on the fracture aperture/growth and a thermo-hydro-mechanical approach for fractured geothermal systems.

In the previously mentioned studies, different numerical discretization techniques were adopted. A popular variational approach to fracture was introduced in [18] (with according numerics presented in [8]) and which was complemented with thermodynamical arguments in [4] and [33, 31]. In the last study, the alternative name phase-field fracture was given. The extension of phase-field to hydraulic fractures in poroelasticity was first undertaken in [35] (originally published as ICES preprint [34]). Therein the solution of the fluid equation (of Darcy type) was assumed to be given, yielding a given pressure for the poroelastic mechanics step. On the other hand, phase-field models involving thermal effects in pure elasticity or plasticity have been studied in [11, 32, 30].

The setting presented in [35] constitutes the starting point to formulate the mechanics step of a phase-field fracture model in thermo-poroelasticity with pressure and temperature as given quantities. For the temperature contribution, an elliptic diffusion equation is considered from which a time-dependent decline constant can be derived; see Hagoort [23]. This parameter controls the thermo-poroelastic stress (also called back-stress). Based on analytical formulas derived by Sneddon and Lowengrub [43, 44] for the crack opening displacements (COD, aperture) for pressurized fractures in an elastic medium, Hagoort extended these equations to non-isothermal situations [23]. A very detailed derivation and discussions of the application in thermo-poroelasticity can be found in Tran et al. [45].

In this paper, we follow the ideas outlined in [45] and develop a phase-field model for pressurized and non-isothermal fractures. Here, we use the ideas from [35] and carefully derive interface conditions between the thermo-poroelastic medium and the fracture. The resulting models (i.e., the energy functional and the related Euler-Lagrange equations) are then analyzed in detail from a mechanical perspective. Here, the focus is on thermodynamical arguments and strain-energy splitting, which are also complemented with corresponding numerical tests.

Our proposed model is implemented into the adaptive parallel framework developed in [24, 50, 29] with most recent results on its computational performance provided in [25]. A key purpose in the current work is on a detailed code verification with respect to the chosen discretization, robustness and efficiency. Therefore, some of our settings (in 2D and 3D) are compared with manufactured solutions developed in [45]. Moreover, in some tests, a computational analysis of the effects of strain-energy splitting is performed in order to study the influence in pressurized and non-isothermal fracture propagation. Iteration numbers of the linear solver and the evolution of the nonlinear residual and the primal-dual active set (enforcing the crack irreversibility constraint) method are undertaken to study the performance of the numerical algorithms.

The outline of this paper is as follows: In Section 2 a phase-field model accounting for given pressures and given temperatures is derived. Then, in Section 3 the model is augmented with thermodynamical arguments and then analyzed from a mechanical perspective. The final model and the numerical solution are discussed in Section 4. Several tests demonstrating our developed model are presented in Section 5. Here, comparisons to analytical solutions, grid refinement studies, calculation of quantities of interest such as the aperture, adaptive studies and the performance of the nonlinear and linear solvers are provided. In Section 6 our main findings are summarized.

2 Non-isothermal and pressurized phase-field modeling

In the following, let $\Omega\in\mathbb{R}^{d}$, $d=2,3$ be a smooth open and bounded set. In $\Omega$, a lower dimensional fracture is denoted by $\mathcal{C}\in\mathbb{R}^{d-1}$. We assume either Dirichlet boundaries conditions $\partial\Omega_{D}:=\partial\Omega$ and Neumann condition on $\partial_{N}\Omega:=\Gamma_{N}\cup\partial\mathcal{C}$, where $\Gamma_{N}$ denotes the outer domain boundary and $\partial\mathcal{C}$ the crack boundary. Let $I:=(0,T)$ denote the loading/time interval with $T>0$ being the end time value.

Prototype configurations of the setting are illustrated in Fig. 1. Using a phase-field approach, the surface fracture $\mathcal{C}$ is approximated in $\Omega_{F}\subset\Omega\in\mathbb{R}^{d}$. The intact region, where with no fracture denoted as $\Omega_{R}:=\Omega\backslash\Omega_{F}\subset\Omega\in\mathbb{R}^{d}$ such that $\Omega_{R}\cup\Omega_{F}=\Omega$ and $\Omega_{R}\cap\Omega_{F}=\varnothing$. It has to be noted, that $\Omega_{F}$, i.e. the domain in which the smeared crack phase-field is approximated, and its boundary $\partial\Omega_{F}$ strongly depend on the choice of the phase-field regularization parameter, i.e. $\varepsilon>0$.

We deal with a multi-field problem depending on the displacement field $\bm{u}:\Omega\rightarrow\mathbb{R}^{d}$ and the crack phase-field $\varphi:\Omega\rightarrow[0,1]$ for a given pressure $p:\Omega\rightarrow\mathbb{R}$ and a given temperature $\Theta:\Omega\rightarrow\mathbb{R}$. The limiting values of $\varphi$, namely, $\varphi=1$ and $\varphi=0$ represent the undamaged and fully broken material phases, respectively. Finally, we, denote the $(\cdot,\cdot)$ the standard $L^{2}$ scalar product in the domain $\Omega$ and $\langle\cdot,\cdot\rangle$ the $L^{2}$ scalar product on a boundary part.

2.1 The energy functional in the thermo-poroelastic medium

We start with the following energy functional posed in the thermo-poroelastic medium:

[TABLE]

where $\sigma$ is the stress tensor specified below, $\bm{\varepsilon}({\bm{u}})$ is the linearized strain tensor and $\bm{\tau}$ traction forces on the boundaries $\partial_{N}\Omega:=\Gamma_{N}\cup\partial\mathcal{C}$. For the specific definition of the thermo-poroelastic stress, we apply the constitutive expression derived in [15]:

[TABLE]

Here, $\bm{\sigma}_{0}$ is the initial stress, $\alpha_{B}$ is Biot’s constant, $p$ is the pressure, $p_{0}$ is the initial pressure, $\bm{I}$ is the second order identity tensor, $3\alpha_{\Theta}$ is the volumetric skeleton thermal dilation coefficient (or thermal expansion coefficient). Then, $K_{d}:=\frac{2}{d}\mu+\lambda$ ($d\in\{2,3\}$ still refers to the dimension of problem), $\Theta$ is the temperature, and $\Theta_{0}$ is the initial temperature, and

[TABLE]

where $\mu$ and $\lambda$ are positive material parameters.

Using these definitions and setting $\bm{\sigma}_{0}=0$ (for the convenience of the reader), we can write the energy functional in the thermo-poroelastic domain:

[TABLE]

This functional represents the energy in the thermo-poroelastic domain only. For the total energy, we need to account for the crack energy as well. In the presence of a fracture $\mathcal{C}$, we need to add the fracture energy which is expressed through [18]:

[TABLE]

where $G_{c}$ is the critical elastic energy restitution rate (related to the stress intensity factor through Irvins formula) and $\mathcal{H}^{d-1}$ is the $d-1$-dimensional Hausdorff measure.

Then, we arrive at the total energy functional:

[TABLE]

This functional is composed by contributions of two disjunct domains, namely $\Omega$ and $\mathcal{C}$. For the numerical treatment we regularize (2.1) following [8]. Specifically, the crack energy is approximated through a sequence of elliptic problems, so-called Ambrosio-Tortorelli functionals [2, 3]. Therein, $\mathcal{H}^{d-1}$ is regularized by introducing an additional smoothed indicator variable nowadays called phase-field: $\varphi:\Omega\to[0,1]$. Finally, we account for the crack irreversibility constraint that the crack can only grow:

[TABLE]

For stating the variational formulations, we now introduce three sets:

[TABLE]

As typical in problems with inequality constraints (see e.g., [27, 28]), $W_{in}$ is a nonempty, closed, convex, subset of the linear function space $W$. Due to the inequality constraint (4), $W_{in}$ is not anymore a linear space.

Moreover, we recall that we do not deal with a classical time-dependent problem, but a quasi-static load-incremental evolution, for which the loading interval is discretized using the discrete time (loading) points

[TABLE]

With these preparations, the energy functional (2.1) can be written as:

Formulation 2.1** (Regularized energy functional).**

Let $p_{0},p,\Theta,\Theta_{0}$ to be given with the initial conditions $\bm{u}_{0}=\bm{u}(x,0)$ and $\varphi_{0}=\varphi(x,0)$. For the loading increments $n=1,2,\ldots,N$, find ${\bm{u}}:={\bm{u}}^{n}\in V$ and $\varphi:=\varphi^{n}\in W_{in}$ such that:

[TABLE]

Remark 2.1**.**

Approximation of (2.1) without pressures $p,p_{0}$ and temperatures $\Theta,\Theta_{0}$ has been used in many studies for solid mechanics.

Remark 2.2**.**

In (2.1), $\kappa:=\kappa(h)$ is a positive regularization parameter for the elastic energy and $\varepsilon:=\varepsilon(h)$ is a regularization parameter for the phase-field variable denoting the width of the transition zone in which $\varphi$ changes from [math] to $1$; see Fig. 1.

The stationary points of the energy functional (2.1) are characterized by the first-order necessary conditions, namely the so-called Euler-Lagrange equations, which are obtained by (formal) differentiation in the variables $\bm{u}$ and $\varphi$:

Formulation 2.2** (Euler-Lagrange equations).**

Let $p_{0},p,\Theta,\Theta_{0}$ be given with the initial conditions $\bm{u}_{0}=\bm{u}(x,0)$ and $\varphi_{0}=\varphi(x,0)$. For the loading increments $n=1,2,\ldots,N$, find ${\bm{u}}:={\bm{u}}^{n}\in V$ and $\varphi:=\varphi^{n}\in W_{in}$ such that:

[TABLE]

and

[TABLE]

2.2 Incorporating fracture pressure and temperature

The energy functional (2.1) and the corresponding Euler-Lagrange equations in Formulation 2.2 are still incomplete because only mechanisms in the thermo-poroelastic medium are taken into account so far. In the following, we consider pressure and temperature variations in the fracture as well, acting as additional normal stresses on the interface between the thermo-poroelastic medium and the fracture. Hence, we obtain a modification of the normal stress term $\langle{\bm{\tau}},{\bm{u}}\rangle_{\partial_{N}\Omega}$.

Similar to [35] in which an interface law was derived to account for pressurized fractures, we now derive an expression including temperature effects. For the convenience of understanding, we also recapitulate the steps of the pressure and develop the modified terms for the temperature simultaneously. Assuming that the fracture is a zone of high permeability and its width much smaller than its length, the leading order of the pressure stress in $\mathcal{C}$ is

[TABLE]

Similarly, we assume (according to [45]) that the leading order of the temperature stress in $\mathcal{C}$ is given by

[TABLE]

The constant $C_{\Theta}$ is obtained by following [45] to get a relation for the temperature near the fracture. We briefly sketch the basic idea. Working with the heat conduction equation, while neglecting convection-dominated terms, we obtain:

[TABLE]

with the thermal diffusivity $\kappa_{\Theta}=\frac{K_{r}}{\rho_{r}C_{r}}$, where $\rho_{r}$ is the rock density, $C_{r}$ is the specific heat, $K_{r}$ thermal conductivity. With the help of the heat conduction equation, a parameter measuring the gradient of the temperature at the fracture interface is derived. This law provides an heuristic argument for the temperature evolution without solving the full temperature equation in the whole domain. Specifically,

[TABLE]

where $\beta$ is a linear thermal expansion coefficient, $E_{Y}$ is a Young’s modulus, $\nu_{s}$ is a Poisson’s ratio. The parameter $\lambda_{\Theta}$ is derived by Hagoort [23] based on the heat conduction equation. Then,

[TABLE]

where $1\leq\gamma_{T}\leq\frac{4}{\pi}$.

2.3 Interface laws for pressure and temperature

We pursue in the following by deriving an interface relationship between the thermo-poroelastic medium $\Omega$ and the fracture $\mathcal{C}$. As usually required for coupled problems with interfaces, we need a kinematic and a dynamic coupling condition. The kinematic conditions are to enforce

[TABLE]

where $p_{f}$ is the fracture pressure, $p_{R}$ the thermo-poroelastic pressure, $\Theta_{f}$ the temperature in the fracture and $\Theta_{R}$ the thermo-poroelastic temperature. Since we enforce continuity of these variables on $\mathcal{C}$ we do not distinguish anymore in the rest of the paper and simply use $p$ and $\Theta$.

The dynamic coupling condition represents the continuity of normal stresses on the crack boundary $\partial\mathcal{C}$ and can be written as

[TABLE]

Since the fracture boundary is smeared and not explicitly known, it remains to discuss how the interface condition is imposed. The kinematic conditions are Dirichlet-like conditions and build into the function spaces as usually done. The dynamic conditions are re-written as domain integrals by using Gauss’ divergence theorem. We recall that the entire boundary is composed as

[TABLE]

We now manipulate the last term with the help of Eq. 9 and the divergence theorem, reads

[TABLE]

These terms can be combined with the other terms containing pressures and temperatures in the energy functional (2.1):

[TABLE]

In total we obtain (noticing that all signs change because we deal with $-\langle\bm{\sigma}\bm{n},{\bm{u}}\rangle_{\partial_{N}\Omega}$ in (2.1)):

[TABLE]

To complete, we introduce a phase-field representation of Eq. LABEL:eq_inter1 and weight the domain terms with $\varphi^{2}$ (see Section 2 in [35] for pressurized fractures only) in (2.1), which then yields:

Formulation 2.3** (Energy functional including fracture pressure and temperature).**

Let $p_{0},p,\Theta,\Theta_{0}$ be given with an initial condition $\bm{u}_{0}=\bm{u}(x,0)$ and $\varphi_{0}=\varphi(x,0)$. For the loading increments $n=1,2,\ldots,N$, find $\bm{u}:=\bm{u}^{n}\in V$ and $\varphi:=\varphi^{n}\in W_{in}$ such that:

[TABLE]

where $\tilde{\bm{\tau}}={\bm{\tau}}-(p-p_{0}){\bm{n}}+C_{\Theta}(\Theta-\Theta_{0}){\bm{n}}$.

Remark 2.3**.**

We underline again that $p$ and $\Theta$ are fixed and given in this paper and therefore not solution variables. The extension in which $p$ is also an unknown was first proposed in [36]. A similar development for non-isothermal phase-field fractures in thermo-poroelasticity is ongoing work.

Remark 2.4**.**

We notice that the temperature opens the crack if $\Theta<\Theta_{0}$ (see also Remark 3.3), i.e., the current (or injected) temperature is lower than the initial temperature.

Remark 2.5**.**

When normal stresses (traction forces) are prescribed on (parts of) the boundary, the term $\tilde{\bm{\tau}}={\bm{\tau}}-(p-p_{0}){\bm{n}}+C_{\Theta}(\Theta-\Theta_{0}){\bm{n}}$ must be carefully considered; see [42] or [35][Section 5.2].

3 Thermodynamic arguments and mechanical analysis

In this section, we further augment our proposed model with thermodynamic arguments and subsequently analyze our model from a mechanical point of view.

3.1 Extension to a decoupled strain-energy function into volumetric and isochoric response

Since the fracturing material behaves quite differently in bulk and shear parts of the domain, we employ a consistent split for the strain energy density function, i.e.

[TABLE]

Hence, instead of dealing directly with $\bm{\varepsilon}(\bm{u})$, we perform additive decomposition of strain tensor into volume-changing (volumetric part) and volume-preserving (deviatoric part), i.e.

[TABLE]

Here, the volumetric strain is denoted as $\bm{\varepsilon}^{vol}(\bm{u}):=\frac{1}{3}(\bm{\varepsilon}(\bm{u}):\bm{I})\bm{I}$ and the deviatoric strain is denoted as $\bm{\varepsilon}^{dev}(\bm{u}):=\mathbb{P}:\bm{\varepsilon}$. The fourth-order projection tensor $\mathbb{P}:=\mathbb{I}-\frac{1}{3}\bm{I}\otimes\bm{I}$ is introduced to map the full strain tensor to its deviatoric counterpart. Therein, \mathbb{I}_{i,j,k,l}:=\frac{1}{2}\big{(}\delta_{i,k}\delta_{j,l}+\delta_{i,l}\delta_{j,k}\big{)} is the fourth-order symmetric identity tensor. Furthermore, $\mathbb{P}$ possesses the major symmetries, i.e. $\mathbb{P}_{i,j,k,l}=\mathbb{P}_{k,l,i,j}$, and $\mathbb{P}^{n}=\mathbb{P}$ for any given integer $n$. So, a decoupled representation of the strain-energy function into a so-called volumetric and deviatoric contribution are given as follows,

[TABLE]

Therein, the volumetric contribution of the strain energy density function reads,

[TABLE]

where $K_{d}:=\frac{2}{d}\mu+\lambda$ is the bulk modulus and $d\in\{2,3\}$. The deviatoric contribution of the strain energy density function is

[TABLE]

To show that equality holds in Eq.12, the identities $\bm{\varepsilon}^{dev}(\bm{u}):\bm{I}=0$ and $\bm{\varepsilon}^{vol}(\bm{u}):\bm{\varepsilon}^{dev}(\bm{u})=0$ are used. Physically, it is trivial to assume that the degradation induced by the phase field acts only on the tensile and shear counterpart of the elastic strain density function. Hence, it is assumed there is no degradation in compression, which also prevents interpenetration of the crack lips during crack closure [1]. It turns out that the modified strain energy density function for the fracturing material becomes,

[TABLE]

such that a monotonically decreasing quadrature degradation function, i.e.

[TABLE]

describes the degradation of the solid with the evolving crack phase-field parameter $\varphi$, includes $\kappa$ that is a small residual stiffness that is introduced to prevent numerical problems. Additionally, the positive part of the strain energy density function that is the tensile and deviatoric part of full strain energy density function reads

[TABLE]

Therein, $H{{}^{+}}(\nabla.\bm{u})$ is a positive Heaviside function such that if $\nabla.\bm{u}$ is positive return one and otherwise, give zero value. It is noted due to identity $\nabla.\bm{u}=\text{tr}(\bm{\varepsilon})$, the positive Heaviside function indicates the points in a domain where they are in tensile part, i.e. $H{{}^{+}}(\nabla.\bm{u})=1$ and compression part, i.e. $H{{}^{+}}(\nabla.\bm{u})=0$. Negative strain energy density function that is the compression part of the full strain energy density function, is

[TABLE]

The constitutive equation for the modified strain energy density function, whereas the stress response constitutes an additive split of $\bm{\sigma}_{\bm{\varepsilon}}:=\frac{\partial\Psi(\bm{\varepsilon})}{\partial\bm{\varepsilon}}$ to purely tensile contribution, i.e. ${\bm{\sigma}}^{+}_{\bm{\varepsilon}}({\bm{\varepsilon}})$ and a purely compression contribution, i.e. ${\bm{\sigma}}^{-}_{\bm{\varepsilon}}({\bm{\varepsilon}})$, reads

[TABLE]

Therein,

[TABLE]

The decoupled representation of the fourth-order elasticity tensor to relate the work into conjugate pairs of stress and strain tensor by means of additive decomposition of the stress tensor, reads

[TABLE]

with

[TABLE]

where the identity $\frac{\partial\bm{\varepsilon}^{dev}(\bm{u})}{\partial\bm{\varepsilon}(\bm{u})}=\mathbb{P}$ is used.

Formulation 3.1** (Final energy functional with strain-energy density splitting).**

Let $p_{0},p,\Theta_{0},\Theta$ be given with the initial conditions $\bm{u}_{0}=\bm{u}(x,0)$ and $\varphi_{0}=\varphi(x,0)$. For the loading increments $n=1,2,\ldots,N$, find ${\bm{u}}:={\bm{u}}^{n}\in V$ and $\varphi:=\varphi^{n}\in W_{in}$ such that:

[TABLE]

Remark 3.1**.**

In the case of elastic cracks, it can be shown that the phase field unknown satisfies $0\leq\varphi\leq 1$. In order to establish this property for the spatially continuous incremental problem, we need to modify energy functional for the negative values of $\varphi$. Hence, similar to [35][Section 2], we have used $\varphi_{+}$ rather than $\varphi$ in terms where negative $\varphi$ could lead to incorrect physics in the bulk energy, traction, pressure and thermal forces.

Remark 3.2**.**

We briefly mention alternative descriptions for the phase-field approximation and the degradation function. First, using a more general term regarding fracture term is denoted as $\frac{G_{c}}{4c_{w}}(\frac{1}{2\varepsilon}W(\varphi)+\frac{\varepsilon}{2}\|\nabla\varphi\|^{2})$ for the constant $c_{w}$ and $W(\varphi)$, which refers to the local part of the dissipated fracture energy functional. This is the so-called $\operatorname{\texttt{AT-1}}$ fracture model if $c_{w}=\frac{2}{3}$ and $W(\varphi)=\|1-\varphi\|$, see [37], and it is $\operatorname{\texttt{AT-2}}$ for $c_{w}=\frac{1}{2}$ and $W(\varphi)=\|1-\varphi\|^{2}$; see [10]. Second, the quadratic polynomial $g(\varphi)$ can be written in the form of the cubic polynomial as $g(\varphi):=3(1-\kappa)\varphi^{2}-2\varphi^{3}+\kappa$ or quartic polynomial as $g(\varphi):=4(1-\kappa)\varphi^{3}-3\varphi^{4}+\kappa$, see [12].

For $\mathcal{E}_{\varepsilon}({\bm{u}},\varphi)$ given in Formulation 3.1, we derive the characterizing Euler-Lagrange equations by differentiating the energy functional with respect to ${\bm{u}}$ and $\varphi$:

Formulation 3.2** (Final Euler-Lagrange equations).**

Let $p_{0},p,\Theta_{0},\Theta$ be given with the initial conditions $\bm{u}_{0}=\bm{u}(x,0)$ and $\varphi_{0}=\varphi(x,0)$. For $n=1,2,3,\ldots,N$, find $\bm{u}:=\bm{u}^{n}\in V$:

[TABLE]

and find $\varphi:=\varphi^{n}\in W_{in}$:

[TABLE]

3.2 The Euler-Lagrange equations in a strong form

In order to complete our derivations, we derive the strong form of Formulation 3.3 in this section. Using integration by parts, we obtain a quasi-stationary elliptic system for the displacements and the phase-field variable, where the latter one is subject to an inequality constraint in time and therefore needs to be complemented with a complementary condition:

Formulation 3.3** (Strong form of the Euler-Lagrange equations).**

Let a pressure $p:\Omega\rightarrow\mathbb{R}$ and temperature $\Theta:\Omega\rightarrow\mathbb{R}$ and the initial conditions $\bm{u}_{0}=\bm{u}(x,0)$ and $\varphi_{0}=\varphi(x,0)$ be given. For the loading steps $n=1,2,3,\ldots,N$, we solve a displacement equation where we seek $\bm{u}:=\bm{u}^{n}:\Omega\rightarrow\mathbb{R}^{d}$

[TABLE]

[TABLE]

The phase-field system consists of three parts: the PDE, the inequality constraint, and a compatibility condition (in fracture mechanics called Rice condition [39]). Find $\varphi:=\varphi^{n}:\Omega\rightarrow[0,1]$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

3.3 Global balance principle of continuum thermo-mechanics

In this section, thermodynamical consistency for the preservation of the principal of balance of energy is shown by considering a sequence of variational substitutions. As a point of departure, considering Formulation 3.1, we define the total free energy functional as

[TABLE]

with the bulk free energy functional

[TABLE]

and the crack free energy functional

[TABLE]

Here, we have introduced the second order crack surface density function per unit volume of the thermo-poroelastic media as,

[TABLE]

We recall that weak forms derived for the displacement and phase-field are given in Formulation 3.2. To derive a global balance of energy, we choose as test functions $\bm{w}=:\dot{\bm{u}}$ and $\psi=:\dot{\varphi}$ and restate Formulation 3.2 for the mechanical part as,

[TABLE]

We define the second order displacement gradient denoted as $\bm{H}:=\nabla{\bm{u}}$, then Eq. (30) reduces to

[TABLE]

that is

[TABLE]

[TABLE]

Herein, the rate of internal mechanical power which describes the response of a domain $\Omega$ done by the stress field is

[TABLE]

Accordingly, through Formulation 3.2 for the phase-field part, we derive

[TABLE]

It is important to note that the inequality $A_{2}(\cdot)(\cdot)$ in Formulation 3.2 becomes an equality in Eq. 36 because we consider the situation in which the inequality constraint (24) is strictly fulfilled; namely $\dot{\varphi}<0$. In this case equality must hold in (23) (and thus in $A_{2}(\cdot)(\cdot)$) since otherwise the compatibility condition (25) is not fulfilled. We can restate Eq. 36 as follows,

[TABLE]

where the left-hand side is considered to be a functional of the rate of the crack phase-field and that is the global crack dissipation functional, i.e. $\dot{\mathcal{E}}_{crack}(\varphi)$, hence

[TABLE]

Through the second term in Eq. 37 and considering Eq. 38, the local form of crack dissipation reads,

[TABLE]

and the third term of Eq. 37 leads to the additional local condition, i.e.

[TABLE]

Herein, $\beta_{\varphi}$ introduced as a crack deriving force is conjugate to the phase-field variable, that is

[TABLE]

In our formulation, in comparison to Miehe et al. [31] we have $\beta_{\varphi}=-f$ (due to the definition of the crack phase-field like the damage variable in the Miehe et al. [31]). Note that $\beta_{\varphi}$ becomes zero in the fracture surface, i.e. $\mathcal{C}$, as $\varphi$ becomes zero.

Remark 3.3**.**

We now look into more detail in the positivity condition for the Eq. LABEL:N20 and discuss term by term. It is noted due to the condition obtained from Eq. 38 which is the positivity of the global crack dissipation functional, and by considering third term in Eq. 37 leads to the positivity of the crack deriving force, i.e. $\beta_{\varphi}\geqslant 0$ and $\dot{\varphi}\leqslant 0$, i.e. Eq. 40. Hence, for the evolving fracture, by means of Eq. LABEL:N20, the first term including $\Psi(\bm{\varepsilon}(\bm{u}))$ is positive due to the normalization condition and for the second and fourth terms we are demanding for pressure $p\geqslant p_{0}$ and thermal $\Theta\leqslant\Theta_{0}$ both conditions hold (i.e. $\beta_{\varphi}\geqslant 0$) to satisfy positiveness of the global crack dissipation functional (it is assumed the gradient term for both pressure and thermal part are neglected and $\alpha_{B}\ll 1$). It has to be noted, due to decoupled strain energy density function we have used (see Section 3.1), results in $\nabla.\bm{u}>0$ in $\Omega_{F}$ (i.e., the transition zone). This can be observed later in Fig. 15 and 18. Additionally, $\dot{\varphi}\leqslant 0$ is imposed to the total energy functional through the primal-dual active set strategy (see Section 4.3). By doing so, the first condition in Eq. 39 is automatically satisfied (due to satisfaction of the third term in Eq. 37. This results in the positivity of the second term in Eq. 37) if and only if both $p\geqslant p_{0}$ and $\Theta\leqslant\Theta_{0}$ hold.

To derive the Karush-Kuhn-Tucker conditions, we consider Eq. 37, and the functional derivative with respect to the $\varphi$ defined as:

[TABLE]

that is

[TABLE]

Due to $\frac{\partial\gamma_{\varepsilon}(\varphi,\nabla\varphi)}{\partial\nabla\varphi}.\bm{n}=\varepsilon\nabla\varphi.\bm{n}=0$, that is Euler-Lagrange crack phase-field equation on the boundary of the given domain. The right hand side of Eq. 37, i.e. $-\int_{\Omega}\beta_{\varphi}\dot{\varphi}\,\mathrm{d}{\textbf{x}}$, can also be obtained by means of Eq. 34 and hence Eq. 37 is restated as,

[TABLE]

that is a global form of the balance of energy for the coupled two-field problem describing the evolution of internal energy and crack dissipation energy, i.e. $\dot{\mathcal{E}}_{crack}(\varphi)+P^{int}$, in a system due to the external loads, i.e. $P^{ext}$. Additionally, by means of the second and the third terms in Eq. 37, it turns out that

[TABLE]

which is the balance law for the evolution of the crack phase-field which ensuring the principal of maximum dissipation during the crack phase-field evolution (see e.g. [31]) and so-called as a compatibility condition; we refer the reader also to Section 3.2. Observe that one may satisfy this global irreversibility constraint of crack evolution, i.e. Eq. 45 by ensuring locally a positive variational derivative of the crack surface function and a positive evolution of the crack phase field, i.e.

[TABLE]

The former condition is ensured in the subsequent treatment by a constitutive assumption that relates the functional derivative to a positive energetic driving force. The latter constraint is a natural assumption that relates the fracture phase field for the non-reversible evolution of crack phase-field.

Remark 3.4**.**

It is noted within loading state, i.e. $\dot{\varphi}<0$, due to the compatibility condition, i.e. Eq. 45, along with Eq. 46, one may observe $\beta_{\varphi}+G_{c}\delta_{\varphi}\gamma_{\varepsilon}(\varphi,\nabla\varphi)=0$ and in the unloading state, i.e. $\dot{\varphi}=0$, we have $\beta_{\varphi}+G_{c}\delta_{\varphi}\gamma_{\varepsilon}(\varphi,\nabla\varphi)<0$. Equation 46 along with Eq. 45 refer to the Karush-Kuhn-Tucker conditions for the phase-field fracturing problem.

Remark 3.5**.**

It turns out that, if the positivity of the crack deriving force is satisfied (see Remark 3.3), i.e. $\beta_{\varphi}\geqslant 0$, and according to Eq. 46 we have $G_{c}\delta_{\varphi}\gamma_{\varepsilon}(\varphi,\nabla\varphi)\leqslant-\beta_{\varphi}$ and hence results to $\delta_{\varphi}\gamma_{\varepsilon}(\varphi,\nabla\varphi)\leqslant 0$, that is the first inequality condition shown in Eq. 39.

4 Numerical solution and the final discrete model

In this section, we briefly describe spatial discretization first. The solution algorithm is then based on a quasi-monolithic approach for which a Newton solver is employed as described in [24]. The crack irreversibility condition in Eq. 4 is treated with a primal-dual active set method. Both techniques can be gathered into one single combined Newton solver [24, 29, 25].

4.1 Spatial discretization

The computational domain is subdivided into quadrilateral or hexahedral element domains. Both subproblems are discretized with a Galerkin finite element method using $H^{1}$-conforming bilinear (2D) or trilinear (3D) elements, i.e., the ansatz and test space uses $Q_{1}^{c}$-finite elements, e.g., for details, we refer readers to the [14]. Consequently, the discrete spaces have the property $V_{h}\subset V$ and $W_{h}\subset W$.

4.2 Problem statement of the compact minimization problem

The starting point for the discrete solution is Formulation 3.1, written in compact form as:

[TABLE]

For the following, we set $U=(\bm{u},\varphi)\in V\times W$. Discretizing

[TABLE]

with the time step size $\delta t:=t^{n+1}-t^{n}$, the incremental problem can be rewritten as

[TABLE]

where $\Phi=0\times W$, so that the constraint acts on the phase-field variable only, and $\bm{U}^{old}$ is the solution from the last time-step (or the initial condition). The minimization of (47) is numerically challenging, due to the following reasons:

the energy functional $\mathcal{E}_{\varepsilon}({\bm{u}},\varphi)$ may admit several local minimizers. Thus finding the global minimum is in general non-feasible, e.g., [9] for discussions on the pure elasticity case. The existence of a minimizer for pressurized fractures was established in [35]. Since the temperature is a given quantity, as the pressure, the proof for existence of the current $\mathcal{E}_{\varepsilon}({\bm{u}},\varphi)$ goes along the same lines as in [9], 2. 2.

the irreversibility of crack phase-field, i.e. $|\mathcal{C}_{t-1}|\leq|\mathcal{C}_{t}|$, is required to provide a thermodynamically consistent minimization problem by having a positive crack dissipation inequality and enforcing on the temporal derivative of the phase-field function, see e.g. [33], 3. 3.

the minimization problem is characterized by localization of the crack phase-field in bands of width of order $\varepsilon$. From a practical and numerical analysis point of view, $\varepsilon>h$ must hold (at least one element has to be existed to cover regularized phased-field). In more detail, $h=o(\varepsilon)$. The regularization parameter is typically a very small dimensionless value and for the accurate fracture response (i.e. converge toward the sharp crack profile) which should tend to [math] in the limit $h\rightarrow 0$ to resolve the bands, see e.g. [20], 4. 4.

the linear system of equations arises from Hessian matrix of the $\mathcal{E}_{\varepsilon}({\bm{u}},\varphi)$ are typically badly conditioned due to the presence of crack phase-field localizations band where the elastic stiffness varies rapidly from the intact value to zero, see e.g. [17].

4.3 A combined Newton method: treating crack irreversibility and

solving the nonlinear problem

In following, we will address how to resolve the issues mentioned above. To this end, we first describe Newton’s method for solving the unconstrained minimization problem $\min\mathcal{E}_{\varepsilon}(\bm{U})$ in Eq. 47 for the total energy functional $\mathcal{E}_{\varepsilon}(\bm{U})$ given in formulation 3.1. We construct a sequence ${\bm{U}}^{0},{\bm{U}}^{1},\dots,{\bm{U}}^{N}$ with

[TABLE]

where the update $\delta{\bm{U}}^{k}$ is computed as the solution of the linear system (details in Section 4.4):

[TABLE]

If we assume the constraints on the phase-field on Eq. 47 hold for the initial guess ${\bm{U}}^{0}$ (we will start with the solution from the last time step, which satisfies the constraint), the condition

[TABLE]

implies that ${\bm{U}}^{k+1}={\bm{U}}^{k}+\delta{\bm{U}}^{k}\leq{\bm{U}}^{k}\leq\dots\leq{\bm{U}}^{0}\leq{\bm{U}}^{old}\textup{ on }\Phi$. In a variational formulation, the previous Newton method reads:

[TABLE]

where ${\bm{\Psi}}:=[\bm{w},\psi]$ is denoted as the total test function. Crack irreversibility is taken care of by removing the corresponding rows and columns in which the constraint is is active. This yields then a reduced system. In our implementation, we combine two Newton methods (active set and the nonlinear iteration for the PDE solution; see Section 4.4) into a single update loop with variable $\delta{\bm{U}}^{k}$. This Newton loop contains a back-tracking line search to improve the convergence radius. This yields Algorithm 1.

Remark 4.1** (Stopping criteria).**

We remark that the above algorithm has two stopping criteria that must be achieved simultaneously:

[TABLE]

Remark 4.2**.**

It is important to distinguish between the full residual ${\bm{R}}({\bm{U}}^{k}_{h})$ and $\widetilde{\bm{R}}({\bm{U}}_{h}^{k})$. The latter is the residual on the inactive set, which can be computed by eliminating the active set constraints from the former.

4.4 On the Jacobian and the residual inside Newton’s method

For solving the PDE problem at each Newton step, we focus on a monolithic scheme in which all equations are solved simultaneously resulting in one semi-linear form. However, it is well known that the energy functional (3) is not convex simultaneously in both solution variables $\bm{u}$ and $\varphi$; but separately in each variable while keeping the other fixed. Consequently, solving the Euler-Lagrange equations in a straightforward way is not possible and influences the robustness and efficiency of the solution scheme because of an (possibly) indefinite Hessian matrix $G$ [19, 48, 49].

The critical terms arise in the Hessian matrix $G$, are the cross terms, includes for mechanical term i.e. \big{(}((1-\kappa){\varphi_{+}^{2}}+\kappa)\bm{\sigma}^{+}_{\bm{\varepsilon}}(\bm{u}),{\bm{\varepsilon}}(\bm{w})\big{)}, for pressure term i.e. $(\alpha_{B}-1)((p-p_{0}){\varphi_{+}^{2}},\nabla\cdot\bm{w})$, $(\nabla(p-p_{0})\varphi_{+}^{2},\bm{w})$ and for thermal term i.e. $(3\alpha_{\Theta}K+C_{\Theta})((\Theta-\Theta_{0}){\varphi_{+}^{2}},\nabla\cdot\bm{w})$ and $(C_{\Theta}\nabla(\Theta-\Theta_{0}){\varphi}^{2},\bm{w})$. To this end, according to [24], for having a convex energy functional, $A_{1}$ is linearized in the direction of the $\bm{u}$ and $\varphi$ by linear extrapolation and time-lagging of the phase-field, i.e. $\varphi\approx\tilde{\varphi}:=\varphi(\varphi^{n-2},\varphi^{n-1})$ in order to obtain a convex energy functional, that is

[TABLE]

Here, $\varphi^{n-2},\varphi^{n-1}$ denote the solutions to previous time steps, denoted as $t^{n-2}$ and $t^{n-1}$, see Fig. 2.

Hence, linearization of $A_{1}(\bm{u})(\bm{w})$ (see Formulation 3.2) in the direction of $\varphi$, i.e. $\delta\varphi$, within tangent stiffness matrix (the Hessian matrix) is neglected (due to the given $\tilde{\varphi}$).

In the following, we state monolithic formulations for the displacement-phase-field system. The presentation is similar to [24]. Specifically, the phase-field variable is time-lagged in the first term of the displacement equation. For fully monolithic solution algorithms and their performance we refer the reader to [48, 49, 19]. In the following, we deal with equalities since we removed the inequality constraint via the primal-dual active set strategy. The residual reads:

[TABLE]

The corresponding Jacobian is built by computing the directional derivative $A^{\prime}({\bm{U}})(\delta{\bm{U}},{\Psi})$. Then:

[TABLE]

4.5 On the linear equation system at each Newton step

4.5.1 Spatial discretization and block structure

In this section, we consider the structure and solution of the linear discrete system (50) arising in each Newton step. For spatial discretization, we use the previously introduced spaces $V_{h}\times W_{h}$ with vector-valued basis

[TABLE]

where the basis functions are primitive (they are only non-zero in one component), so we can separate them into displacement and phase-field basis functions and sort them accordingly:

[TABLE]

where $N_{\bm{u}}+N_{\varphi}=N_{s}$. This is now used to transform (50) into a system of the form

[TABLE]

where ${\bm{M}}$ is a block matrix (the Jacobian) and $F$ the right-hand side consisting of the residuals. The block structure is

[TABLE]

with entries coming from (54):

[TABLE]

Remark 4.3**.**

Since we replaced $\varphi^{2}$ by $\tilde{\varphi}^{2}$ in the displacement equation, the block $M^{{\bm{u}}\varphi}_{i,j}$ is zero and the Hessian matrix $\bm{G}$ has triangular structure. In the other case, all blocks would be nonzero; see [48, 49].

The right-hand side consists of the corresponding residuals (see semi-linear form (53)).

In particular, we have

[TABLE]

The criterion for convergence of the Newton method is based on the relative residual norm that is, $\texttt{Residual}:\|\bm{F}(\bm{x}_{k+1})\|\leq\texttt{Tol}_{\texttt{N-R}}\|\bm{F}(\bm{x}_{k})\|$ for the user-prescribed $\texttt{Tol}_{\texttt{N-R}}$.

In the matrix, the degrees of freedom that belong to Dirichlet conditions (here only displacements since we assume Neumann conditions for the phase-field) are strongly enforced by replacing the corresponding rows and columns as usual in a finite element code. In a similar fashion, the rows and columns that belong to nodes of the active set are removed from the matrix. Corresponding right-hand side values in the vector $F$ are set to zero; see as well Table 1.

4.5.2 Solution of the linear equation systems

The linear system arising at each Newton step is solved iteratively using a generalized minimal residual (GMRES) scheme with a block diagonal preconditioner, i.e. ${\bm{P}}^{-1}$, as follows

[TABLE]

where we approximate the blocks using a single V-cycle of algebraic multigrid (by Trilinos ML, [26]), and due to the in-dependency to the mesh element size, the whole solver scheme is nearly optimal as recently demonstrated for pressurized phase-field fractures in [25]. The main reason why this solver works nicely is that the block-diagonal terms are of elliptic type.

5 Numerical studies

The performance of the proposed model is investigated by means of a series of representative numerical examples in two- and three spatial dimensions. The proposed formulation is considered to be a canonically consistent and robust scheme for treating pressurized fractures and non-isothermal setting in thermo-poroelastic media. The emphasis is on verifications of the programming code against analytical solutions obtained with Sneddon-Lowengrub’s formula [44] (extended to non-isothermal configurations in [45]) and using realistic material properties. Mesh refinement studies are provided to study the accuracy of our approach. Investigations of the effects strain-energy splitting and solver outcomes are studied in addition. The implementation is based on deal.II [5, 6] and specifically on the programming code from [25].

5.1 Geometries and parameters

5.1.1 Two-dimensional setup

In the first five examples, we present two-dimensional settings. The geometrical parameter denoted as $a$ in Fig. 3 is set to 100, and, hence a reservoir of size is $(0,200)^{2}\;m$. The Tran et al. problem is considered in an infinite domain with $\mathcal{C}$ with a pre-defined crack of length $2l_{0}$ in the $y=a$ plane and is restricted in $a-l_{0}\leq|\mathcal{C}|\leq a+l_{0}$. We set the half crack length as a $l_{0}=10\;m$. We set the initial values for displacement and phase-field as $\bm{u}_{0}:=0\in\Omega$ and $\varphi_{0}:=1\in\Omega_{R}$ and $\varphi_{0}:=0\in\Omega_{F}$. The (finite) computational domain is subdivided into quadrilateral element domains. It is noted, that all numerical examples are computed by parallel computing on $8$ processors; see Fig. 3[b]. The different sub-domains are associated with different processors. Depending on mesh refinement, the workload for each processor is adjusted dynamically at each time step. Typically, we set $2\times 10^{4}$ quadrilateral elements for each processor.

5.1.2 Three-dimensional setup

In the last two examples, we present three-dimensional test cases. The configuration setup is displayed in Fig. 4 and the values of the material parameters are the same as in the two-dimensional test cases. The geometrical parameter denoted as $a$ in Fig. 4 is set to 50, hence, the reservoir of size is a $(0,100)^{3}\;m$ cube. We consider an infinite domain with $\mathcal{C}$ as a pre-defined crack with a given radius $l_{0}$ in the $y=a$ plane and is restricted in $\mathcal{C}(x,y):=\{(x,y)\in\Omega:(x-a)^{2}+(z-a)^{2}\leq l^{2}_{0}\}$. We set the crack radius length as a $l_{0}=10\;m$. The (finite) computational domain is subdivided into hexahedral element domains.

5.1.3 Material parameters

The material parameters are taken from [45] and related to practical field problems. The geo-mechanical parameters are Young’s modulus $E_{Y}=1.5\times 10^{10}Pa$, Poisson’s ratio $\nu=0.15$, and the linear thermal-expansion coefficient $\beta=10^{-5}1/C$ and finally the critical energy release rate $G_{c}=10^{10}\;N/m$.

For the flow and temperature, we have the thermal diffusivity $\kappa_{\Theta}=10^{-6}m^{2}/s$, the initial stress $p_{0}=12130\>$KPa and in the stationary-in-length test cases (only evolution of the aperture), we use the constant pressure $\bar{p}=15834\>$KPa. Additionally, for the non-isothermal setting, the initial temperature $\Theta_{0}=100\;C$, the apparent temperature is set to $\Theta=70\>C$.

For all examples, we set $\alpha_{B}=0$ and $\alpha_{\Theta}=0$. All material properties are fixed for the following numerical examples, unless indicated otherwise.

5.1.4 Model parameters

The phase-field parameters are chosen as $\kappa=10^{-10}$, and $\varepsilon=10\times\sqrt{h}$ (respecting the condition $h<\varepsilon$) in 2D and $\varepsilon=2h$ in the 3D test cases. Using the adaptive mesh refinement, we are not dealing with a uniform mesh and hence the domain is divided into coarser and finer mesh elements, i.e. $\Omega=\Omega_{c}\cup\Omega_{f}$. Let $h^{max}_{c}:=max(h_{c})$ in $\Omega_{c}$ and $h^{min}_{f}:=min(h_{f})$ in $\Omega_{f}$. On the one hand $\varepsilon$, enters in the constitutive modeling (that is in our PDE). On the other hand, it implicitly depends on the discretization of a domain, i.e. $h^{max}_{c}=o(\varepsilon_{c})$ and $h^{min}_{f}=o(\varepsilon_{f})$. Thus we define, $\varepsilon:=max(\varepsilon_{c},\varepsilon_{f})$ resulting in $\varepsilon\geqslant h_{c}>h_{f}$ and thus $\varepsilon>h$ holds in every point of the domain.

The stopping criterion of the Newton method, i.e. the relative residual norm that is $\texttt{Residual}:\|\bm{F}(\bm{x}_{k+1})\|\leq\texttt{Tol}_{\texttt{N-R}}\|\bm{F}(\bm{x}_{k})\|$, is set to $\texttt{Tol}_{\texttt{N-R}}=10^{-10}$.

The threshold value for the local predictor-corrector mesh refinement is $\texttt{TOL}_{\varphi}=0.9$. That is, we refine when

[TABLE]

5.2 Test scenarios

We propose the following numerical tests:

• •

Case a. Two-dimensional problem with constant pressure and without thermal effects.

• •

Case b. Two-dimensional problem with constant pressure and fixed Hagoort’s decline constant value. That is time-independent problem and considered to validate proposed formulation for fixed thermal effect.

• •

Case c. Two-dimensional problem with constant pressure and evolving Hagoort’s decline constant through time. That is a time-dependent problem and used to validate proposed formulation for a temperature variation in time.

• •

Case d. Observing crack propagation in the two-dimensional problem, higher temperature difference is considered. That is time-dependent problem with a constant pressure and evolving Hagoort’s decline constant through the time.

• •

Case e. Observing crack propagation in the two-dimensional problem, evolving in time for pressure and Hagoort’s decline constant are considered. This is a time-dependent problem.

• •

Case f. Three-dimensional problem with constant pressure and a fixed Hagoort’s decline constant value.

• •

Case g. Three-dimensional problem with increasing pressure and Hagoort’s decline constant are considered.

5.3 Quantities of interest

In our numerical tests, we study the following aspects:

• •

Crack opening displacements (COD, also known as aperture) for the first three cases:

[TABLE]

where ${\bm{n}}:={\bm{n}}^{+}=-{\bm{n}}^{-}$ is the normal vector (see Fig. 5 b) and $x_{0}$ the $x$-coordinate of the integration line. We note that the integration is perpendicular to the crack direction. Here, the crack is aligned with the $x$-axis and therefore integration into the normal direction coincides with the $y$-direction. The analytical solution for the COD derived by Tran et al. [45] reads:

[TABLE]

where $l_{0}$ is the half-length of the crack and $0<\rho<1$, for $\rho=x/l_{0}$ with $x$ being the distance to the origin of the crack. In three dimensions, the formula reads:

[TABLE]

For pressurized fracture without any thermal effect, i.e., $\Theta=\Theta_{0}=0$, we obtain the formulas derived in Sneddon and Lowengrub [44][Section 2.4 and Section 3.3].

• •

Fracture length and path for propagating fractures.

• •

Both stopping criteria of the nonlinear solver (Section 4.3), i.e. the relative residual norm and the active set constraint.

• •

For some fixed time steps, the average number of GMRES iterations within per Newton cycle.

• •

The evolution of the mechanical strain energy functional, i.e. mechanical term in Formulation 3.1 and dissipated fracture energy functional, i.e. the fracture term in Formulation 3.1.

• •

Comparisons regarding no split and vol./dev. split of strain density function.

5.4 Case a. Sneddon-Lowengrub’s 2D setting with pressure

In the first example, we consider a stationary setting with a pressurized fracture. The geometrical setup is the same as illustrated in Fig. 3 and the material parameters are listed in Section 5.1.3. Related results were presented in [35, 47, 7].

We set the time step size $\delta t=1\>s$ for a given constant pressure $\bar{p}=10^{-3}\>KPa$. This test case is computed in a quasi-stationary manner: that is, we solve several pseudo-time steps. The goals of this test are to observe the crack opening displacement, and in particular to show the crack tip approximation. The initial mesh is five times uniformly refined, then, local mesh adaptivity is applied for five further levels.

The COD findings for both our computational model and the analytical calculation are shown in Fig. 6. It is observed the crack tips at $90m$ and $110m$ show excellent agreement with the analytical solution and convergence towards Sneddon-Lowengrub’s manufactured solution.

5.5 Case b. 2D setting with pressure and constant temperature

In this second example, we turn our attention now to a pressurized, non-isothermal configuration. Specifically, we follow Tran’s et al. [45] 2D problem with a constant pressure and a fixed Hagoort’s decline constant. We keep the geometry, all parameters, and tolerances as in the first example. Additionally, we set $\lambda_{\Theta}=10^{-4}$ constant and we work still in a full stationary setting where the time step loop terminates after one step. The crack opening displacement both obtained through phase-field fracturing modeling and the analytical calculation are shown in Fig. 7 for the same refinement levels as in the previous example.

The crack phase-field pattern for different time steps on the locally adaptive refined meshes are illustrated in Fig. 8.

5.6 Case c. 2D setting with pressure, temperature and decline constant

While Hagoort’s decline constants were constant in Case b, we now account for a temperature variation in time. We keep the geometry, all parameters, and tolerances as in the first example. The only two changes are that we set the time step size to one day $\delta t=86400s$ and observe $T=365$ days. This test case has a slight time-dependence since we use $\lambda_{\Theta}:=\lambda_{\Theta}(t)$ being time-dependent and defined through Eq. 8. This scenario allows us to study the fracture behavior for a constant pressure like in the production processes.

Our numerical findings regarding the maximum COD is shown in Fig. 9. We consider different levels for uniformly refined meshes. Findings for level 5 and level 9 are shown in Fig. 9 a and b, respectively. In between we computed the levels 6-8 as well, but are not shown to keep the figures legibly. We observe excellent agreement of our numerical model with the analytical solution.

Another perception is to observe again the COD over the fracture (similar to Case a and Case b). This can be shown by fixing the time steps (i.e. 50, 150 and 365 days). Our findings are illustrated in Fig. 10 and show again excellent agreement.

5.7 Case d. 2D setting for the crack growth due to temperature variation

Having studied slight time dependencies with variations in the COD, we now consider a first test case with a propagating fracture. The propagation will be caused by temperature variations. For temperature variations, we use $\lambda_{\Theta}:=\lambda_{\Theta}(t)$. Hagoort’s decline constants for this example are shown in Fig. 11 a. The second novelty in this numerical test are comparisons of strain-energy split and the non-split version. The third goal are very detailed studies on the performance of the linear and nonlinear solvers.

Again, we keep most parameters as in the first example except a higher temperature difference to $220$ degree Celsius (i.e. $\Theta_{0}=300\>C\;\text{and}\;\Theta=80\>C$) and the initial stress $p_{0}=12130\>$KPa and constant pressure as mentioned earlier is set to $\bar{p}=15834\>$KPa.

In order to facilitate fracture propagation, the critical energy release rate is drastically reduced to $G_{c}=5.5\times 10^{5}\;N/m$. The time step size is $\delta t=86400s$. As it is shown in Fig. 11 a, we are now interested in four specific time steps, namely $100,114,121$ and $136$ days.

5.7.1 Crack propagation and locally refined meshes

Figure 12 displays a sequence of crack patterns during different time steps without considering any strain-energy split, i.e. Formulation 2.3 and considering vol./dev. split, i.e. Formulation 3.1, of the strain energy density function. In fact, until time step $100$, we observe almost no growth. Then, we have crack propagation towards the boundaries.

The functionality of predictor-corrector mesh refinement is shown in Fig. 13 to present the evolution of the locally refined mesh when the crack is growing. As mentioned earlier in Section 5.1.4, $\mathrm{Tol}_{\varphi}$ is set to 0.9. Each cell that has at least one support point with value $\varphi(x_{i})<\texttt{Tol}_{\varphi}$ for $x_{i}\in\Omega_{e}$ will be refined unless we are already at the maximum desired refinement level is reached.

5.7.2 Analysis of the strain-energy splitting

We now turn to the second goal that is novel in this form in the published literature regarding pressurized/non-isothermal fracture settings with phase-field modeling.

A qualitative representation of the tensile and compression counterparts are shown in Fig. 14. We define $\theta_{\bm{u}}:=H{{}^{+}}(\nabla.\bm{u})=H{{}^{+}}(tr(\bm{\varepsilon}))$ that is an indicator to detect the compression region, i.e. $\theta_{\bm{u}}<0$, and, reversily, to detect the tensile region, i.e. $\theta_{\bm{u}}>0$. If the vol./dev. split of the strain energy density function is used, we are in the crack region when the phase-field energy exceeds its critical value $G_{c}$ and additionally $\theta_{\bm{u}}>0$. This gives us an additional constraint to our physical model. Hence, only the tensile stress is degraded and compression stress is free from the effect of the phase-field variable; see Eq. 19.

Figure 15 compares the incremental elastic strain energy, i.e. the mechanical term in the Formulation 3.1, and the crack dissipated energy, i.e. fracture term in Formulation 3.1. For comparison, the strain energy density function without split and vol./dev. split are shown. As observed in Fig. 15, it turns out that dissipated fracture energy resulting from the crack phase-field evolution evolve with a sharp ascending behavior. This clearly shows that the sharp transition between the intact state and the fully cracked state of the solid body arised. Thereafter, dissipated fracture energy evolves towards a steady-state response. While the solid body completely fractured, no change in the dissipated fracture energy is observed anymore.

5.7.3 Performance of nonlinear and linear solvers

A detail analysis of the convergence of the primal-dual active set framework and the Newton solver are shown in Fig. 16 for the finest mesh level.

Finally, the average number of GMRES iterations of linear iterations per Newton cycle for different times steps is shown in Table 2. Specifically, the preconditioner (Section 4.5.2) works extremely well.

5.8 Case e. 2D setting for pressurized crack propagation in non-isothermal settings

In this example, we extend the previous Case d to study now crack propagation due to temperature effects , i.e. $\lambda_{\Theta}:=\lambda_{\Theta}(t)$ and increasing pressure, i.e. $p:=p(t)$. From the application viewpoint this example is closer to a fracture process with a characteristic time scale much smaller than in Case d. Hence, we set the time step size to one day $\delta t=60s$ and simulate $270$ minutes. The increasing pressure is prescribed as:

[TABLE]

We set the constant value $\alpha:=\frac{1}{2}$ and $t$ is the time and given by $t=n\delta t$ with $n$ being the time step number. We keep all remaining parameters as in the first example and we set $\Theta_{0}=100\;C$ and $\Theta=70\;C$ degree Celsius (recall $\Theta<\Theta_{0}$, see Remark 3.3).

5.8.1 Crack propagation and locally refined meshes

As shown in Fig. 17 b, we are now interested in four time steps, namely 240, 252, 259 and 267 minutes to observe the crack path. There again (as in Case d) only slight differences at the same points, suggesting that strain-splitting has nearly no effect in these configurations.

5.8.2 Analysis of the strain-energy splitting

A qualitative representation of the tensile region, i.e. $\theta_{\bm{u}}>0$, and compression counterpart, i.e. $\theta_{\bm{u}}<0$, are displayed in Fig. 18.

Based on the results, shown in Fig. 18 (and also Fig. 14), an alternative scalar variable for adaptive mesh refinement could be $\theta_{\bm{u}}:=H{{}^{+}}(\nabla.\bm{u})$ rather than using the phase-field variable. If $\theta_{\bm{u}}=1$, the material tends to tensile stress and hence crack propagation might occur and if $\theta_{\bm{u}}=0$ the material is under compression stress and we are away from the fracture zone. The main advantages of using $\theta_{\bm{u}}$ rather than $\varphi\leq\operatorname{\texttt{TOL}_{\varphi}}$ are two-fold. First, the above suggestion is based on the (physical) displacement field (and not on phase-field) and secondly, there is no requirement about choosing the correct threshold value for $\operatorname{\texttt{TOL}_{\varphi}}$.

5.8.3 Performance of nonlinear and linear solvers

As in Case d, we study the solver performances. The evolution of the primal-dual active set and the Newton residuals are shown in Fig. 19.

The average number of GMRES iterations of linear iterations per Newton cycle for different times steps is shown in Table 3. Specifically, the preconditioner (Section 4.5.2) works again extremely well.

5.9 Case f. 3D setting with pressure, temperature and decline constant

In this test, we extend Case c to a three-dimensional setting. The geometry and the material parameters are described in the Sections 5.1 and 5.6, respectively. The maximum width evolution is displayed in Figure 20 and shows a good fit with Tran’s et al. [45] manufactured 3D solution for computing the maximal COD.

5.10 Case g. 3D setting for the Pressurized crack propagation in isothermal and non-isothermal settings

In this final example, the 2D Case e is extended to three dimensions. The geometry and material parameters are listed in the Sections 5.1 and 5.8. As in the $2D$ test case, the time step size is $\delta t=60s$. We compute $27$ loading (time) steps. The fracture radii are displayed in Figure 21. Graphical illustrations of the phase-field variable, i.e., the crack path, are displayed in Figure 22.

6 Conclusions

In this work, we developed a phase-field fracture model for pressurized and non-isothermal cracks in thermo-poroelasticity. We concentrated on the so-called mechanics step in which pressure and temperature are given quantities. Using the thermo-poroelastic stress tensor and interface laws, we derived an energy functional and the corresponding Euler-Lagrange equations. We then proceeded with a careful analysis from the mechanical point of view. As computational framework, we used a programming code that allows for adaptive mesh refinement and is fully parallelized. Therein, the nonlinear problem is formulated within a monolithic framework and is solved with a semi-smooth Newton method accounting as well for the crack irreversibility constraint. Based on these computational capacities, we ran in total seven test scenarios. Therein, we compared to analytical (manufactured) solutions taken from the literature with an excellent agreement of our numerical solutions. Moreover, we carried out mesh refinement studies and gave further insight into the performances of the linear and nonlinear solvers. In these studies, we also analyzed the effects on crack growth and solver behavior when strain-energy splitting is used or not. The outcome of our numerical simulations let us come to the conclusion that we have developed a robust and efficient computational framework for pressurized and non-isothermal configurations when the pressure and temperatures are given. In future studies, one idea is to treat the temperature as unknown to be determined by the governing partial differential equation. Here, we could proceed analogously to fluid-filled phase-field fracture models proposed in [36]. Another task will be more computations of practical problems with multiple interacting fractures.

Acknowledgments

This work is supported by the German Research Foundation, Priority Program 1748 (DFG SPP 1748) named Reliable Simulation Techniques in Solid Mechanics. Development of Non-standard Discretization Methods, Mechanical and Mathematical Analysis in the sub-project Structure Preserving Adaptive Enriched Galerkin Methods for Pressure-Driven 3D Fracture Phase-Field Models (WI 4367/2-1).