Linear Thermo-Poroelasticity and Geomechanics

Abstract

Most engineering applications estimate the deformation induced by loads by using the linear elasticity theory. The discretization process starts with the equilibrium equation and then develops a displacement formulation that employs the Hooke’s law. Problems of practical interest encompass designing of large structures, buildings, subsurface deformation, etc. These applications require determining stresses to compare them with a given failure criteria. One often tackles this way a design or material strength type of problems. For instance, Geomechanics applications in the oil and gas industry assess the induced stresses changes that hydrocarbon production or the injection of fluids, i.e., artificial lift, in a reservoir produce in the surrounding rock mass. These studies often include reservoir compaction and subsidence that pose harmful and costly effects such as in wells casing, cap-rock stability, faults reactivation, and environmental issues as well. Estimating these stress-induced changes and their consequences require accurate elasticity simulations that are usually carried out through finite element (FE) simulations. Geomechanics implies that the flow in porous media simulation must be coupled with mechanics, which causes a substantial increase in CPU time and memory requirements.

1. Introduction

This chapter presents a continuous Galerkin FE formulation for linear isotropic elasticity. It covers in detail how to derive such formulation beginning with the equilibrium equation and the virtual work statement. It also discretizes the continuity equation for slightly compressible single-phase flow to show how to couple different physics with elasticity. It discusses several coupling approaches such as the monolithic and iterative ones, i.e., loosely coupled. This chapter also mentions the affinity of the poroelastic case with the thermoelastic one. It thus also includes thermoelasticity in the treatment herein. It shows concrete numerical examples covering two- and three-dimensional problems of practical interest in thermo-poroelasticity. The sample problems employ triangular, quadrilateral, and hexahedral meshes and include comments about implementing boundary conditions (BCS). An introduction to domain decomposition ideas such as iterative coupling by the BCS, i.e., Dirichlet-Neumann domain decomposition and mortar methods for non-matching interfaces is included.

The treatment herein demonstrates that the continuous Galerkin formulation for linear isotropic elasticity is the foundation to develop codes for mechanics. Indeed, after discretizing linear elasticity is straightforward to extend the implementation to nonlinear mechanics such as rate-independent plasticity. It thus provides some comments about such extension. Applications of practical interest show that industrial size problems will require domain decomposition techniques to handle such simulations in a timely fashion. Unquestionably, domain decomposition techniques can exploit current parallel machines architectures which brings high-performance computing into the picture. For instance, recently the author showed that the Dirichlet-Neumann scheme could handle problems at the reservoir field-level as well as the mortar method decoupled by this last one. Its current results are backed up by papers published in peer-reviewed journals and conferences thus this book chapter summarizes that effort.

2. Mathematical model for thermo-poroelasticity

This section discusses the governing equations for linear homogeneous isotropic thermo-poroelasticity and their FE formulation. It skips details for the sake of brevity thus a more detailed treatment can be found in [1, 2, 3, 4]. The mathematical formulation considers a bounded domain $\mathrm{\Omega }\subset {\mathbb{R}}^n,n=2,3$ and its boundary is $\mathrm{\Gamma }=\partial \mathrm{\Omega }$ , and a time interval of interest $\text{]}0,\Im \text{[}$ . Let ${\mathcal{T}}_h$ be a non-degenerate, quasi-uniform conforming partition of $\mathrm{\Omega }$ composed of triangles or quadrilaterals for two-dimensional problems, and hexahedra or tetrahedra for three-dimensional problems. For instance, Gai [5] thesis showed that deformable porous media, i.e., the reservoir matrix, the single-phase flow model equation derives from the continuity equation, i.e., a mass balance statement, for slightly incompressible single-phase flow and Darcy’s law which yields:

$\frac{\partial {\varphi }^{\ast }}{\partial t}+\nabla \cdot \left(-\frac{1}{\mu }\underset{¯}{\underset{¯}{K}}\left(\nabla p-\mathit{\rho g}\nabla z\right)\right)=q,$ E1

where the equation’s parameters are ${\varphi }^{\ast }$ , a model specific porosity, $\underset{¯}{\underset{¯}{K}}$ represents the absolute permeability tensor. The dynamic viscosity is $\mu$ , while $\rho$ is the fluid density, as well as $g$ , is the gravity acceleration constant, $p$ is the fluid pressure, and $q$ represents sources and sinks. This latter notation is standard in fluid mechanics and reservoir simulation. Finally, the algorithmic porosity ${\varphi }^{\ast }$ is defined by:

${\varphi }^{\ast }={\varphi }^0+\alpha \cdot \left(\nabla \cdot \underset{¯}{u}-{\epsilon }_v^0\right)+\frac{1}{M}\left(p-p^0\right),$ E2

where the additional parameters are accordingly $\alpha$ which is the Biot’s constant, $\underset{¯}{u}$ represents the displacement vector, while ${\epsilon }_v^0$ is the initial volumetric strain. Herein $M$ is the Biot’s modulus [6], while ${\varphi }^0$ and $p^0$ define for a reference or initial state. The common BCS for the pressure equation imply Neumann or no-flow namely:

$\nabla p\cdot \underset{¯}{{\displaystyle \widehat{n}}}=0\mathrm{on}\mathrm{\Gamma },$ E3

one should also consider an initial or reference pressure distribution in the whole domain. Sources and sinks simulate injector and producer wells, respectively. Herein $\underset{¯}{\widehat{n}}$ is the outer unitary normal vector as usual. For the mechanics part, one begins from the equilibrium equation for a quasi-steady process, i.e., Newton second law, which means that one discards the acceleration term:

$\begin{array}{cc}-\nabla \cdot \underset{¯}{\underset{¯}{\sigma }}& =\underset{¯}{f}\text{in}\mathrm{\Omega };\mathrm{\Gamma }={\mathrm{\Gamma }}_D^u\cup {\mathrm{\Gamma }}_N^u\\ \underset{¯}{u}& =\underset{¯}{0}\mathrm{on}{\mathrm{\Gamma }}_D^u\\ \underset{¯}{t}& =\underset{¯}{\underset{¯}{\sigma }}\cdot \underset{¯}{\widehat{n}}\mathrm{on}{\mathrm{\Gamma }}_N^u\end{array}$ E4

where $\underset{¯}{\underset{¯}{\sigma }}$ is the stress tensor, $\underset{¯}{f}$ corresponds to the vector of body forces, such as gravity and electromagnetic effects, for instance. One can decompose BCS in Dirichlet type, i.e., ${\mathrm{\Gamma }}_D^u$ , and Neumann type BCS, i.e., ${\mathrm{\Gamma }}_N^u$ , where the external tractions are known or prescribed. Hooke’s law combined with Biot’s poroelastic theory defines $\underset{¯}{\underset{¯}{\sigma }}$ by the following expression:

$\underset{¯}{\underset{¯}{\sigma }}=\underset{¯}{\underset{¯}{C}}:\underset{¯}{\underset{¯}{\epsilon }}-\left[\alpha \left(p-p^0\right)+3\mathit{K\beta }\left(T-T^0\right)\right]\underset{¯}{\underset{¯}{\delta }};\underset{¯}{\underset{¯}{C}}=\lambda \underset{¯}{\underset{¯}{\delta }}\otimes \underset{¯}{\underset{¯}{\delta }}+2G\underset{¯}{\underset{¯}{I}},$ E5

where $T=T\left(x,t\right)$ is the temperature, $\underset{¯}{\underset{¯}{C}}$ is the elastic moduli, $\beta$ corresponds to the coefficient of thermal dilatation while $K$ is the bulk modulus. The Kronecker delta becomes $\underset{¯}{\underset{¯}{\delta }}$ while $\lambda$ , and $G$ , are the Lamé constants, and $\underset{¯}{\underset{¯}{I}}$ represents the fourth-order identity tensor. The strain tensor $\underset{¯}{\underset{¯}{\epsilon }}$ is given by:

$\underset{¯}{\underset{¯}{\epsilon }}={\nabla }^s\underset{¯}{u}=\frac{1}{2}\left[\nabla \underset{¯}{u}+{\left(\nabla \underset{¯}{u}\right)}^T\right].$ E6

One can derive a weak form by substituting Eq. (2) into Eq. (1) and then multiplying by a test function $v\in H_0^1\left(\mathrm{\Omega }\right)$ and integrating over $\mathrm{\Omega }$ and using the Gauss-divergence theorem, this yields:

$\begin{array}{c}\underset{\mathrm{\Omega }}{\int }\left(\frac{1}{M}\frac{\partial p}{\partial t}v+\mathit{\alpha v}\nabla \cdot \underset{¯}{\dot{u}}+\frac{1}{\mu }\underset{¯}{\underset{¯}{K}}\cdot \nabla p{\left(\nabla v\right)}^T\right)\cdot \mathit{dx}=\underset{\mathrm{\Omega }}{\int }q\cdot \mathit{vdx}+\\ \underset{\mathrm{\Omega }}{\int }\left(\frac{\mathit{\rho g}}{\mu }\underset{¯}{\underset{¯}{K}}\cdot \nabla z{\left(\nabla v\right)}^T\right)\mathit{dx}+\underset{\partial {\mathrm{\Omega }}_N^p}{\int }v\frac{1}{\mu }\underset{¯}{\underset{¯}{K}}\left(\nabla p-\mathit{\rho g}\nabla z\right)\cdot {\underset{¯}{\widehat{n}}}^T\mathit{ds}.\end{array}$ E7

A weak form for the equilibrium Eq. (4) can be derived in a similar way, by testing against a given virtual displacement, $\underset{¯}{\chi }$ . One arrives at:

$\underset{\mathrm{\Omega }}{\int }{\left(\nabla \underset{¯}{\chi }\right)}^T:\underset{¯}{\underset{¯}{\sigma }}\mathit{d}\mathrm{\Omega }=\underset{\partial {\mathrm{\Omega }}_N^u}{\int }{\underset{¯}{\chi }}^T\cdot \underset{¯}{t}\mathit{ds}+\underset{\mathrm{\Omega }}{\int }{\underset{¯}{\chi }}^T\cdot \underset{¯}{f}\mathit{d}\mathrm{\Omega }$ E8

where $\underset{¯}{t}=\underset{¯}{\underset{¯}{\sigma }}\cdot \underset{¯}{\widehat{n}}$ are the tractions applied as Neumann BCS. This is the well-known virtual work statement. The FE space can be taken as a finite-dimensional subspace of the continuous Sobolev spaces [7], thus:

${\mathcal{C}}_k\left({\mathcal{T}}_h\right)=\left\{v\in L^2\left(\mathrm{\Omega }\right):\forall \mathrm{e}\in {\mathcal{T}}_h,{\left.v\right|}_e\in {\mathbb{P}}_k\left(e\right)\right\}$ E9

where ${\mathbb{P}}_k\left(e\right)$ represents the space of polynomials of total degree less than or equal to $k$ , ${\mathcal{C}}_k\left({\mathcal{T}}_h\right)$ is called test functions that are continuous along the given element’s edges. Let one represents the primary variables in the element $e$ , i.e. displacements and pressure, as nodal values multiplied by shape or interpolation functions [8]:

$p_e^h\left(\underset{¯}{x}\right)={\left({\underset{¯}{\mathrm{\Pi }}}^e\right)}^T\cdot {\underset{¯}{p}}^e;{\underset{¯}{u}}_e^h\left(\underset{¯}{x}\right)={\underset{¯}{\underset{¯}{\mathrm{\Psi }}}}^e\cdot {\underset{¯}{u}}^e$ E10

where ${\underset{¯}{\mathrm{\Pi }}}^e$ and ${\underset{¯}{\underset{¯}{\mathrm{\Psi }}}}^e$ are matrices of shape functions given by:

$\begin{array}{cc}{{\mathrm{\Pi }}_i}^e& ={{\psi }_i}^e\left(\underset{¯}{x}\right)\\ {{\mathrm{\Psi }}_{\mathit{ij}}}^e& =\left\{\begin{array}{l}{{\psi }_k}^e\left(\underset{¯}{x}\right)\text{if}j=\underset{¯}{\mathrm{j}}\\ 0\text{otherwise}\end{array}\right)\\ \underset{¯}{\mathrm{j}}& =\text{nDOF}\cdot \left(k-1\right)+i;k=1\dots \mathit{nn}\end{array}$ E11

here $\mathit{nn}$ is the number of nodes in the given element, $i=1\dots \mathit{nn} j=1\dots \mathit{nn}\cdot \mathit{n}$ and nDOF is the number of degrees of freedom which equals the space dimension, $n$ . Now the engineering strain $\underset{¯}{\widehat{\epsilon }}$ is defined by:

$\underset{¯}{\widehat{\epsilon }}=\underset{¯}{\underset{¯}{B}}\cdot {\underset{¯}{u}}^e;\underset{¯}{\underset{¯}{B}}=\underset{¯}{\underset{¯}{D}}\cdot {\underset{¯}{\underset{¯}{\mathrm{\Psi }}}}^e$ E12

where ${\underset{¯}{\underset{¯}{D}}}_{\left(n\right)}$ , $n=2,3$ is defined as:

${\underset{¯}{\underset{¯}{D}}}_{\left(2\right)}^T=\left[\begin{array}{lll}{\partial }_x& 0& {\partial }_y\\ 0& {\partial }_y& {\partial }_x\end{array}\right];{\underset{¯}{\underset{¯}{D}}}_{\left(3\right)}^T=\left[\begin{array}{llllll}{\partial }_x& 0& 0& {\partial }_y& {\partial }_z& 0\\ 0& {\partial }_y& 0& {\partial }_x& 0& {\partial }_z\\ 0& 0& {\partial }_z& 0& {\partial }_x& {\partial }_y\end{array}\right].$ E13

Finally substituting the generalized Hooke’s law Eq. (5) into Eq. (8) and using Eq. (7) leads to the FE model for linear isotropic poroelasticity, thus:

$\left[\begin{array}{ll}0& 0\\ {\underset{¯}{\underset{¯}{\mathrm{Q}}}}^T& \underset{¯}{\underset{¯}{\mathrm{S}}}\end{array}\right]\frac{d}{\mathit{dt}}\left\{\begin{array}{l}\underset{¯}{u}\\ \underset{¯}{p}\end{array}\right\}+\left[\begin{array}{ll}\underset{¯}{\underset{¯}{\mathrm{K}}}& -\underset{¯}{\underset{¯}{\mathrm{Q}}}\\ 0& \underset{¯}{\underset{¯}{\mathrm{H}}}\end{array}\right]\left\{\begin{array}{l}\underset{¯}{u}\\ \underset{¯}{p}\end{array}\right\}=\left\{\begin{array}{l}\underset{¯}{f_u}\\ \underset{¯}{f_p}\end{array}\right\}.$ E14

One can obtain the loose coupling approach in different ways. Eq. (15) shows one possible choice, where one solves the displacements first by taking the pressures from the previous time step. Next, one updates the pressures by using the newest displacements:

$\begin{array}{cc}\underset{¯}{\underset{¯}{\mathrm{K}}}\cdot {\underset{¯}{u}}^{k+1}& =f_u+\underset{¯}{\underset{¯}{\mathrm{Q}}}\left({\underset{¯}{p}}^k-{\underset{¯}{p}}^0\right)\\ \underset{¯}{\underset{¯}{{\mathrm{S}}^{\prime }}}\cdot {\underset{¯}{p}}^{k+1}& =\underset{¯}{\underset{¯}{{{\mathrm{S}}^{″}}^{}}}\cdot {\underset{¯}{p}}^k+\underset{¯}{f_p}\cdot \mathrm{\Delta }t-{\underset{¯}{\underset{¯}{\mathrm{Q}}}}^T\left({\underset{¯}{u}}^{k+1}-{\underset{¯}{u}}^k\right)\\ \underset{¯}{\underset{¯}{{\mathrm{S}}^{\prime }}}& =\underset{¯}{\underset{¯}{\mathrm{S}}}+\theta \cdot \mathrm{\Delta }t\cdot \underset{¯}{\underset{¯}{\mathrm{H}}}\\ \underset{¯}{\underset{¯}{{{\mathrm{S}}^{″}}^{}}}& =\underset{¯}{\underset{¯}{\mathrm{S}}}-\left(1-\theta \right)\cdot \mathrm{\Delta }t\cdot \underset{¯}{\underset{¯}{\mathrm{H}}},\end{array}$ E15

where expressions for the matrices are provided in Eq. (16) and $\theta$ is the implicitness factor that lies between 0 and 1, while $\mathrm{\Delta }t$ represents the time-step size. One can define an iterative coupling scheme in different ways, but they all derive from the loose coupling scheme with incorporating an internal iteration to update lagged quantities. For further details please refer to [4]. Also notice that for thermal stresses, one can derive an equivalent pressure drop, after Eq. (5), that renders Eq. (15) unchanged.

$\begin{array}{cc}\underset{¯}{\underset{¯}{\mathrm{S}}}& =\underset{\mathrm{\Omega }}{\int }\frac{1}{M}\underset{¯}{\mathrm{\Pi }}\cdot {\underset{¯}{\mathrm{\Pi }}}^T\mathit{dx};\underset{¯}{\underset{¯}{\mathrm{Q}}}=\underset{\mathrm{\Omega }}{\int }{\underset{¯}{\underset{¯}{B}}}^T\alpha {\underset{¯}{\omega }}_{\left(n\right)}\cdot \underset{¯}{\mathrm{\Pi }}\mathit{dx}\\ \underset{¯}{\underset{¯}{\mathrm{K}}}& =\underset{\mathrm{\Omega }}{\int }{\underset{¯}{\underset{¯}{B}}}^T\underset{¯}{\underset{¯}{C}}\underset{¯}{\underset{¯}{B}}\mathit{dx};\underset{¯}{f_u}=\underset{\partial {\mathrm{\Omega }}_N^u}{\int }\underset{¯}{t}\cdot {\underset{¯}{\underset{¯}{\mathrm{\Psi }}}}^T\mathit{ds}+\underset{\mathrm{\Omega }}{\int }{\underset{¯}{\underset{¯}{\mathrm{\Psi }}}}^T\underset{¯}{f}\cdot \mathit{dx}\\ \underset{¯}{\underset{¯}{\mathrm{H}}}& =\underset{\mathrm{\Omega }}{\int }\frac{1}{\mu }\underset{¯}{\underset{¯}{K}}\nabla \underset{¯}{\mathrm{\Pi }}\cdot {\left(\nabla \underset{¯}{\mathrm{\Pi }}\right)}^T\mathit{dx};{\underset{¯}{\omega }}_{\left(2\right)}={\left(1,1,0\right)}^T;{\underset{¯}{\omega }}_{\left(3\right)}={\left(1,1,1,0,0,0\right)}^T\\ f_p& =\underset{\partial {\mathrm{\Omega }}_N^p}{\int }\left(\frac{1}{\mu }\underset{¯}{\underset{¯}{K}}\nabla p\cdot \underset{¯}{n}\right)\cdot \underset{¯}{\mathrm{\Pi }}\mathit{ds}+\underset{\mathrm{\Omega }}{\int }{\underset{¯}{\mathrm{\Pi }}}^{T}q\cdot \mathit{dx}+\underset{\mathrm{\Omega }}{\int }\left(\frac{\mathit{\rho g}}{\mu }\underset{¯}{\underset{¯}{K}}\cdot \nabla \underset{¯}{\mathrm{\Pi }}{\left(\nabla z\right)}^T\right)\mathit{dx}.\end{array}$ E16

This section completes with a comment about the Continuous Galerkin (CG) formulation for the pressure (1). It is well-known that the formulation that was presented above for flow it is not locally mass conservative, and thus the resulting fluxes are not continuous across the element edges. But it is also true that accurate flow simulations require the latter, especially for multi-phase flow, though. Nevertheless, one can utilize post-processing techniques to recover locally conservative mass fluxes [2]. This chapter, though, for convenience has restricted its focus to CG methods for flow but has realized that the coupled formulation may be modified to consider mixed FE methods and finite volumes for flow as well as changing CG by post-processing. The author already showed for the simple flow cases reported herein that CG yields to physical pressure fields that can be employed for geomechanics purposes. The precise numerical comparison among CG and Discontinuous Galerkin (DG) solutions was performed in [2] to demonstrate that CG can compute pressures accurately.

3. Nonlinear heat transfer equation

The transient nonlinear heat conduction in a given body is as follows [9, 10, 11]:

$\begin{array}{cc}\rho C_p\frac{\partial T}{\partial t}& =\nabla \cdot \left(\kappa \nabla T\right)+Q_T\mathit{on}\mathrm{\Omega }\times ]0,\Im [,\\ T& =g\mathit{on}{\mathrm{\Gamma }}_D^T\times ]0,\Im [,\\ n\cdot \left(\kappa \nabla T\right)& =h\mathit{on}{\mathrm{\Gamma }}_N^T\times ]0,\Im ],\\ T\left(x,0\right)& =T_0\left(x\right)\forall x\in \mathrm{\Omega }.\end{array}$ E17

In (17), $C_p$ is the heat capacity to constant pressure and $\kappa =\kappa \left(T\right)$ is the thermal conductivity. $Q_T$ represents heat sources. Neumann BCS imply heat transfer via Fourier’s law: adiabatic or no-flux BCS; $h=0$ of most domain boundaries.

One can derive a FE formulation for model problem (17) by multiplying by a test function and integrate by parts and applying the Gauss-divergence theorem to arrive at the following bilinear form:

$m\left(T,v\right)+k\left(T,v\right)-q\left(Q_T,v\right)-f\left(h,v\right)=0,$ E18

where the functions are:

$\begin{array}{cc}m\left(T,v\right)& =\underset{{\mathrm{\Omega }}^e}{\int }\mathit{v\rho }{C}_p{\partial }_{t}T\cdot \mathit{dx},\\ k\left(T,v\right)& =\underset{{\mathrm{\Omega }}^e}{\int }\kappa {\left(\nabla T\right)}^T\cdot \nabla v\cdot \mathit{dx},\\ q\left(Q_T,v\right)& =\underset{{\mathrm{\Omega }}^e}{\int }v\cdot Q_T\cdot \mathit{dx};f\left(h,v\right)=\underset{{\mathrm{\Gamma }}_h^e}{\int }v\cdot h\cdot \mathit{ds}.\end{array}$ E19

Time discretization renders the local residual for the element $e$ :

$\begin{array}{l}\underset{¯}{R}\equiv \underset{¯}{\underset{¯}{M}}\cdot \left({\underset{¯}{T}}^{\left(\ell \right)}-{\underset{¯}{T}}^{\left(m\right)}\right)+\mathrm{\Delta }t\cdot {\left(\underset{¯}{\underset{¯}{K}}\cdot \underset{¯}{T}\right)}^{\left(m+\theta \right)}\\ -\mathrm{\Delta }t\cdot {\underset{¯}{q}}^{\left(m+\theta \right)}-\mathrm{\Delta }t\cdot {\underset{¯}{f}}^{\left(m+\theta \right)}=\underset{¯}{0},\end{array}$ E20

where the linear operator ${\left(\cdot \right)}^{\left(m+\theta \right)}\equiv \left(1-\theta \right){\left(\cdot \right)}^{\left(t=t^{\left(m\right)}\right)}+\theta {\left(\cdot \right)}^{\left(t=t^{\left(\ell \right)}\right)},$ $\ell =\left(m+1\right), \underset{¯}{\underset{¯}{M}} \underset{¯}{\underset{¯}{K}}$ are the mass and stiffness matrix respectively. Thus the Jacobian is given by:

$\underset{¯}{\underset{¯}{J}}=\frac{\partial \underset{¯}{R}}{\partial {\underset{¯}{T}}^{\left(\ell \right)}}=\underset{¯}{\underset{¯}{M}}+\frac{\partial }{\partial {\underset{¯}{T}}^{\left(\ell \right)}}{\left(\underset{¯}{\underset{¯}{K}}\cdot \underset{¯}{T}\right)}^{\left(t^{\left(\ell \right)}\right)}$ E21

this equation renders once again:

$\underset{¯}{\underset{¯}{J}}=\underset{¯}{\underset{¯}{M}}+\mathrm{\Delta }t\cdot \left(\underset{¯}{\underset{¯}{K}}+\delta \underset{¯}{\underset{¯}{K}}\right)$ E22

if one assumes that $\kappa \left(T\right)=\left(a\cdot T+b\right)$ ; $a,b\in \mathbb{R}$ , then:

$\delta \underset{¯}{\underset{¯}{K}}=\sum _p\frac{\partial K_{\mathit{ip}}}{\partial T_j^{\left(\ell \right)}}\cdot T_p^{\left(\ell \right)},$ E23

where the variation term is given by:

$\frac{\partial K_{\mathit{ip}}}{\partial T_j^{\left(\ell \right)}}=\underset{{\mathrm{\Omega }}^e}{\int }a{\psi }_j{\left(\nabla {\psi }_i\right)}^T\cdot \nabla {\psi }_p\cdot \mathit{dx}.$ E24

One often employs the Newton-Raphson algorithm to solve the linearized system of equations in every time step, namely, $\underset{¯}{\underset{¯}{J}}\cdot {\underset{¯}{\text{Δ}\mathit{T}}}^{\left(\ell \right)}=-\underset{¯}{R}$ . One can utilize the same continuous FE space that where described in Section 2. The reader may refer to [11] for a full treatment.

4. Domain decomposition methods

Domain Decomposition Methods (DDM) encompass highly efficient algorithms to obtain solutions of large-scale discrete problems on parallel super-computers. They mainly consist of partitioning the domain into various subdomains and then getting the global solution through the resolution of the subdomain problems [12, 13] often in an iterative fashion. These methods can be seen as an iterative coupling by the internal and thus unknown BCS. There is a broad literature covering these approaches, and that is why this chapter, therefore, presents a short introduction for the sake of completeness. The recommended references include Bjorstad and Widlund [14], Bramble et al. and Marini and Quarteroni [15], who considered the Dirichlet-Neumann (DN) DDM and Neumann-Neumann.

Let $L$ be an abstract linear differential operator, such as the Laplace operator, for instance. The DN-DDM scheme implies solving a series of problems in the proper sequence that requires a coloring tool (see Figure 1). Let the Dirichlet subdomains be colored in white while the Neumann subdomains are in black. Notice that the interface between subdomains is denoted by $\mathrm{\Gamma }$ . After one provides the initial guess on the primary variables on $\mathrm{\Gamma }$ , i.e., ${\gamma }^k$ must be given, then one can solve the problem on the white subdomains (Dirichlet problems), which corresponds to step 1 in Eq. (25).

$1\left)\left\{\begin{array}{l}{{\mathit{Lu}}_1}^{k+1}=f\text{in}{\mathrm{\Omega }}_1\\ {u_1}^{k+1}=0\mathrm{on}\partial {\mathrm{\Omega }}_1\cap \partial \mathrm{\Omega }\\ {u_1}^{k+1}={\gamma }^k\mathrm{on}\mathrm{\Gamma }\end{array}\right)2\right)\left\{\begin{array}{l}{{\mathit{Lu}}_2}^{k+1}=f\text{in}{\mathrm{\Omega }}_2\\ {u_2}^{k+1}=0\mathrm{on}\partial {\mathrm{\Omega }}_2\cap \partial \mathrm{\Omega }\\ {\partial }_n{u_2}^{k+1}={\kappa }^{k+1}\mathrm{on}\mathrm{\Gamma }\end{array}\right)$ E25

Figure 1.

It depicts the DNDDM.

Let the primary variable be called “displacements” and their gradient “tractions,” i.e., normal derivative in the boundary. Then, the tractions on the interface $\mathrm{\Gamma }$ must be computed after first solving step 1 on the white subdomains. They are then passed through communication to solve the second step on the black subdomains, i.e., Neumann subdomains. On this latter, since the tractions are known on $\mathrm{\Gamma }$ , one can solve for unknown displacements to provide feedback on the next iteration level. Both displacements and tractions are often over-relaxed to improve the convergence rate. The given relaxation parameters, referred in Eq. (26) as ${\theta }^D$ and ${\theta }^N$ , must lie between 0 and 1:

$\begin{array}{cc}{\kappa }^{\left(k+1\right)}& =\left(-{\theta }^N\cdot {\partial }_n{u_1}^{k+1}+\left(1-{\theta }^N\right)\cdot {\partial }_n{u_2}^k\right)\mathrm{on}\mathrm{\Gamma }\\ {\gamma }^{\left(k+1\right)}& =\left({\theta }^D\cdot {u_2}^{k+1}+\left(1-{\theta }^D\right)\cdot {u_1}^k\right)\mathrm{on}\mathrm{\Gamma }.\end{array}$ E26

It happens that this approach requires at least a two-entry coloring tool or even more, i.e., there may be subdomains with mixed interfaces, colored as gray [12]. There is a lack of parallelism in the sense that black subdomains must wait for the white ones to communicate their tractions. An initial guess for tractions should be prescribed to mitigate this issue, but this latter is not feasible in most cases. A straightforward way to obtain an initial estimate for the multiplier ${\gamma }^k$ is by computing the so-called coarse-run that implies solving the same problem in a coarser mesh and interpolating over $\mathrm{\Gamma }$ by using the smaller’s problem FE space. The reader may refer to the literature [16, 17] for further reading and proof of convergence and also revise [2] for a more detailed description that includes implementation details, which this chapter omits for the sake of brevity.

5. The mortar FE method (MFEM)

The primary goal here is to extend MFEM to glue curved interfaces such as the one shown in Figure 2 where MFEM treats non-matching interfaces. The section first introduces a brief description of non-uniform rational B-Splines curves and surfaces (NURBS) in [2, 3, 18]. The reader is referred to those references that cover the topics of computational geometry, in particular how to build these NURBS entities. Let MFEM be described for linear isotropic elasticity regarding bilinear forms, $a$ and $\mathrm{\Upsilon }$ defined in Eq. (27) below [2, 3],

$\begin{array}{cc}a\left(\underset{¯}{u},\underset{¯}{v}\right)& =\underset{\mathrm{\Omega }}{\int }\underset{¯}{\epsilon }{\left(\underset{¯}{v}\right)}^T\cdot \underset{¯}{\underset{¯}{C}}\cdot \underset{¯}{\epsilon }\left(\underset{¯}{u}\right)\mathit{dx};l\left(\underset{¯}{v}\right)=\underset{\partial {\mathrm{\Omega }}^N}{\int }{\underset{¯}{t}}^T\cdot \underset{¯}{v}\mathit{ds}+\underset{\mathrm{\Omega }}{\int }{\underset{¯}{f}}^T\cdot \underset{¯}{v}\mathit{dx}\\ \mathrm{\Upsilon }\left(\underset{¯}{u},\underset{¯}{\mathrm{\Phi }}\right)& =\underset{\mathrm{\Gamma }}{\int }{\left[\underset{¯}{u}\right]}^T\cdot \underset{¯}{\mathrm{\Phi }}\mathit{ds};\left[\underset{¯}{u}\right]=\left({\underset{¯}{u}}^{\left(1\right)}-{\underset{¯}{u}}^{\left(2\right)}\right)\end{array}$ E27

where $\mathrm{\Upsilon }$ stands for the gluing condition among subdomain interfaces and the jump $\left[\underset{¯}{u}\right]$ on the displacements is required to vanish in an integral or “weak” sense, thus:

$\left\{\begin{array}{l}a\left({\underset{¯}{u}}_h,{\underset{¯}{v}}_h\right)+\mathrm{\Upsilon }\left({\underset{¯}{v}}_h,{\underset{¯}{\mathrm{\Lambda }}}_h\right)=l\left({\underset{¯}{v}}_h\right)\\ \mathrm{\Upsilon }\left({\underset{¯}{u}}_h,{\underset{¯}{\mathrm{\Phi }}}_h\right)=0\end{array}\right)$ E28

the parameters in Eq. (28) are as follows: ${\underset{¯}{\mathrm{\Phi }}}_h$ represents the mortar space while ${\underset{¯}{v}}_h$ corresponds to the weighting space and ${\underset{¯}{\mathrm{\Lambda }}}_h$ is the Lagrange multiplier space, i.e., the linear combination of mortar functions, often polynomial functions, and Lagrange multiplier degrees of freedom. Let ${{\overline{\mathcal{T}}}_h}^M$ be a conforming partition of the so-called parametric space, $\overline{\mathrm{\Omega }}$ , whose image serves as the mortar’s geometrical entity, i.e., curve or surface, composed of line-segments ( $n=2)$ or quadrilaterals ( $n=3$ ). One takes the mortar space as a finite-dimensional subspace of the continuous Sobolev spaces, that is:

${\mathcal{C}}_k\left({{\overline{\mathcal{T}}}_h}^M\right)=\left\{\mathrm{\Phi }\in L^2\left(\overline{\mathrm{\Omega }}\right):\forall {\mathrm{e}}^{\mathrm{M}}\in {{\overline{\mathcal{T}}}_h}^M,{\left.\mathrm{\Phi }\right|}_{e^M}\in {\mathbb{P}}_k\left(e^M\right)\right\}$ E29

herein ${\mathbb{P}}_k\left(e^M\right)$ stands for the space of polynomials of total degree less than or equal to $k$ while ${\mathcal{C}}_k\left({{\overline{\mathcal{T}}}_h}^M\right)$ represents test functions that are continuous along the edges of $e^M$ .

Figure 2.

Ω₁ is in the top, Ω₂ is in the bottom, and the interface Γ is the bold curve.

One can write in a matrix or algebraic form, Eq. (28) as:

$\left[\begin{array}{lll}\left[k^{\left(1\right)}\right]& \left[0\right]& {\left[{\mathrm{\Upsilon }}^{\left(1\right)}\right]}^T\\ \left[0\right]& \left[k^{\left(2\right)}\right]& -{\left[{\mathrm{\Upsilon }}^{\left(2\right)}\right]}^T\\ \left[{\mathrm{\Upsilon }}^{\left(1\right)}\right]& -\left[{\mathrm{\Upsilon }}^{\left(2\right)}\right]& \left[0\right]\end{array}\right]\cdot \left[\begin{array}{l}{\underset{¯}{u}}^{\left(1\right)}\\ {\underset{¯}{u}}^{\left(2\right)}\\ \underset{¯}{\mathrm{\Lambda }}\end{array}\right]=\left[\begin{array}{l}{\underset{¯}{l}}^{\left(1\right)}\\ {\underset{¯}{l}}^{\left(2\right)}\\ \underset{¯}{0}\end{array}\right]\cdot$ E30

The equation above corresponds to the so-called “saddle-point problem (SPP).” Notice that subdomains are only connected using the Lagrange multiplier $\underset{¯}{\mathrm{\Lambda }}$ if they happen to be known (it is well-known that for elasticity, the multipliers are the unknown tractions on the interface), then one can decouple the system in Eq. (30) and then one just needs to perform subdomain solves. For the SPP (30), one may match displacements or tractions in the interface. The Dirichlet-Neumann scheme that the section presents is only a particular case of the most general Robin-Robin domain decomposition scheme [2, 3]. The rectangular matrices $\left[{\mathrm{\Upsilon }}^{\left(i\right)}\right],i=1\dots 2$ , are denoted as projectors since they permit to map to and from the given mortar space [2, 3].

The following line integral defines the projector, for 2-D problems, as:

${\mathrm{\Upsilon }}_{\mathit{ij}}^{\left(k\right)}=\underset{\overline{\mathrm{\Omega }}}{\int }{\phi }_j^{\left(k\right)}\left(\xi \right){\mathrm{\Phi }}_i\left(\xi \right)\cdot ‖d_{\xi }\underset{¯}{\mathcal{C}}\left(\xi \right)‖\mathit{d\xi }$ E31

where ${\phi }_j^{\left(k\right)}$ represents the global non-mortar side interpolation functions and ${\mathrm{\Phi }}_i$ are the mortar space basis functions, while $‖d_{\xi }\underset{¯}{\mathcal{C}}‖$ is the length of the tangent vector associated to the B-Spline or NURBS curve. Similarly, 3-D problems imply:

${\mathrm{\Upsilon }}_{\mathit{ij}}^{\left(k\right)}=\underset{\overline{\mathrm{\Omega }}}{\iint }{\phi }_j^{\left(k\right)}\left(\underset{¯}{\xi }\right){\mathrm{\Phi }}_i\left(\underset{¯}{\xi }\right)\cdot ‖{\partial }_{\xi }\underset{¯}{\mathcal{S}}\left(\underset{¯}{\xi }\right)\times {\partial }_{\eta }\underset{¯}{\mathcal{S}}\left(\underset{¯}{\xi }\right)‖\mathit{d\xi }\mathit{d\eta }$ E32

where $‖{\partial }_{\xi }\underset{¯}{\mathcal{S}}\times {\partial }_{\eta }\underset{¯}{\mathcal{S}}‖$ is the norm of the surface’s normal vector. Particular quadrature rules to compute these integrals must be developed. The reader should refer to [3] for a detailed explanation including the proper algorithm in pseudo-code.

6. Numerical examples

The author implemented these FE models in the Integrated Parallel Finite Element Analysis program (IPFA) that is a C++ application whose main characteristics are described in [2, 12]. IPFA employs standard continuous Lagrange polynomials as shape functions for the space discretization in each subdomain, ${\mathbb{P}}_k\left(e\right)$ , as well as mortar space ${\mathbb{P}}_k\left(e^M\right)$ . It also utilizes piecewise linear polynomials for the space discretizations in all examples herein that were run on a MacBook Pro laptop equipped with an Intel(R) Quad-Core(TM) i7-4870HQ CPU @ 2.5 GHz and 16 GB of RAM. The author chose this laptop for the sake of convenience, in particular, the availability of debugging tools free of charge, such as the Microsoft Visual Studio Community. Aside, one can achieve some level of parallelism due to the multi-core technology.

6.1. Example 1: Two-dimensional steady state single-phase flow

The example is a manufactured problem where the solution is a priori chosen. Then, one substitutes the given pressure field in the governing partial differential equation to obtain the source term, i.e., loading, that reproduces the input field. The problem in strong form looks like:

$-\nabla \cdot \left(\underset{¯}{\underset{¯}{K}}\nabla p\right)=f\text{in}\mathrm{\Omega };p=p_0\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{D}}=\mathrm{\Gamma },$ E33

where the domain of interest corresponds to the unitary square and Dirichlet BCS are enforced. The input pressure field is given by:

$p\left(x,y\right)=\mathit{xy}\cdot \left(x-1\right)\cdot \left(y-1\right)\cdot \text{exp}\left[-\left(x^2+y^2\right)\right];\underset{¯}{\underset{¯}{K}}=\underset{¯}{\underset{¯}{\mathrm{I}}}.$ E34

Figure 3 shows the pressure field that corresponds to the problem 6.1 whose discretization encompasses three subdomains: two of them (the top and bottom ones) consist of triangular meshes while the one in the middle was discretized by a regular Cartesian quadrilateral mesh. The top-left corner of the figure shows the mesh that is employed.

Figure 3.

The MFEM solution to problem 6.1.

The pressure field is on the right-top corner, and its horizontal derivative is in the bottom-left corner, while the discrepancy between the numerical and exact solutions, i.e., the absolute error, was rendered in the right-bottom corner. Table 1 represents the number of elements and points of each mesh from top to bottom. The mortars as geometrical entities correspond to two B-Splines interpolants (NURBS with all weights equal 1) that were constructed by interpolating a sinusoidal wave as the figure shows (see [3] for details). Thirty-two quadratic mortar elements per curve were utilized to glue these three subdomains. A direct frontal solver was used to solve the global SPP in Eq. (30) [3]. The results that are summarized on Figure 3 are in good agreement with the analytical solution. The absolute error against the correct answer is also displayed. The discrepancy is of the order of ${10}^{-4}$ . Notice that besides the example only matched the displacements on the interface, a good accordance is also obtained for the horizontal derivative.

Points	Elements	Kind of mesh
980	1814	Triangular
1560	1472	Quadrilateral
4090	7858	Triangular

Table 1.

Meshes for example 6.1.

Whether or not one utilizes the SPP approach, the local problems are completely disconnected. This fact can be exploited to reduce the computational time significantly. Indeed, these sub problems can be handled in separate threads using a shared memory approach, i.e., multi-threading assembling. A convergence analysis was also performed, by successively running refined meshes [3] and by keeping a refinement ratio of 2:1 between subdomains. The exercise used a piecewise quadratic mortar space where the number of mortar elements equals the number of coarse edges in the non-mortar sides. It tackled meshes of size 8, 16, 32, 64, 128 and 256 respectively. Figure 4 displays the resulting convergence rate in a $\mathit{log}-\mathit{log}$ plot. The slope of the least-squares straight line is 1.44143, where the coefficient of determination is $R^2=84\%$ . This slope agrees with the theory that predicts a rate of $\mathcal{O}\left(h^{3/2}\right)$ [2, 3]. However, the resulting slope is slightly lower because of numerical errors, such as quadrature and linear system solving errors.

Figure 4.

Snapshots showing the evolution of the DN-DDM applied to problem 6.1.

Finally, Figure 5 shows pressure snapshots that represent four different Dirichlet-Neumann iteration levels evolving from left-to-right and top-to-bottom. The fact that no initial guess for pressure was provided explains the mismatch in the first snapshot. That is why one needs to eliminate discrepancies by running the process to match up those subdomains, i.e., the traction residual in the interface must vanish, which for this case occurs in just a handful of iterations. The stop criterion precisely involves the residual in the tractions in the interface that is required to fall below the given tolerance. For this particular problem, the iterative process spent six iterations to achieve a residual lower than ${10}^{-6}$ [3].

Figure 5.

The numerical $L^2$ convergence rate for problem 6.1.

6.2. Example 2: Coupled flow and mechanics

This example analyzes a coupled flow and mechanics simulation in a reconstructed reservoir (RS) model with different meshes for the flow and mechanics physics [18]. The author proposed such a reconstruction workflow in [18] which permits this latter feature by computing a projection operator to mapping pressures from the original flow mesh into the so-called reference mechanics mesh. Toward that end, the example employs the slightly compressible flow formulation loosely combined with the mechanics model as shown in Eq. (15). The objective is to show a realistic field level RS compaction and subsidence coupled computation. The goal is thus working three different cases for the mechanics part in which one only changes the resolution of the reconstructed mechanics mesh in the pay-zone while preserving the mechanical properties constant as well as the geometry, BCS, and the depletion scenario. The exercise admits the actual static properties as being in the pay-zone such as porosity $\varphi$ and permeability for the isotropic case $K_x=K_z=K_y$ as shown in Figure 6, whose depiction is three times exaggerated in the $z-$ direction. The numerical values are assumed to be as follows: the fluid viscosity is 0.01325 cp and the total compressibility is $c_t=1.4\times {10}^{-5}{\mathrm{Psi}}^{-1}$ ( $M^{-1}=\varphi \cdot c_t$ ). This example does not incorporate gravity loading for both flow and mechanics.

Figure 6.

The reservoir’s permeability $K_y$ .

Table 2 compiles the mesh dimensions in every direction. The example also contemplates $\mathrm{Nz}=10,\mathrm{Nu}=5,\mathrm{No}=7,\mathrm{Nc}=5$ (mesh patches on the corners and No and Nu stand for over- and under-burden respectively). The table also displays the number of elements, $\mathit{ne}$ , degrees of freedom (DOF) and timing data for all three cases. The example considers 30 vertical producer wells as revealed in Figure 6. The initial condition encompasses a constant pressure of 10,000 Psi in the whole pay-zone while the pressure in the producer wells is set at 5000 Psi. This assumption resembles a depletion scenario. BCS correspond to no-flow on all RS faces for the pressure equation, while Figure 7 depicts BCS for mechanics that are the typical traction free surface on the top and far-field on all remainder planes. Notice that the far-field BCS implies that the displacement in the perpendicular direction to the given plane is zero. The example also assumes a zero initial displacement field.

Case #	Description	Nx	Ny	$\mathit{ne}$	DOF	Assembling time
One	1/4 of RS	35	13	15,960	51,830	0 min, 19 s 75 ms
Two	1/2 of RS	70	26	48,279	159,120	0 min, 59 s 89 ms
Three	1/1 of RS	140	49	156,408	506,160	3 min, 14 s 89 ms

Table 2:

Mesh sizes and simulations in example 6.2.

Figure 7.

The BCS for the mechanics problem in the $x-z$ plane (the pay-zone is highlighted in red).

Figure 8 displays the mechanics mesh. The second case on Table 2 corresponds to a layered RS with Young’s modulus $E_u=3\times {10}^4,E_p=1\times {10}^4,E_o=2\times {10}^4\left[\mathrm{Psi}\right]$ , while Poisson’s ratio, $v=0.25$ , is constant in the whole domain. In Figures 8 through 10 the graphs are 6 times exaggerated in the $z-$ direction for better visualization. The subscript letters symbolize the under-burden, pay-zone and over-burden levels respectively. The goal is representing a more realistic geomechanics model with stiffer surroundings around the RS.

Figure 8.

The hexahedral mesh generated for $2\mathrm{nd}$ case in Table 2.

Figure 9 pictures snapshots with the evolution of the vertical displacement $u_z\left[\mathrm{m}\right]$ and the RS pressure. A compaction dome naturally grows just above the area where the most significant pressure-drop happens. The pattern of deformation is the typical scenario where a compaction-dome rests on the top (blue color) while a build-up occurs in the bottom of the RS (rendered in red color). The deformation caused by the pressure-drop is localized because this reservoir does not entirely drain but is still a compelling case for coupled flow and mechanics.

Figure 9.

Snapshots at 10 and 20 years of evolution showing the vertical-displacement field $u_z$ , the pay-zone displays pressure.

Figure 10 renders pressure-drop snapshots at 10 years of production. Each picture draws the original RS mesh and the reference mechanic’s mesh for all cases that Table 2 covered, from top-to-bottom and left-to-right. Notice that the action of the projection operator improves with the refinement of the reference mechanics mesh as one should expect. The monotone pressure behavior, which does not drastically change across neighboring elements in the original RS mesh, may explain this improvement. Though, some items remain red-colored because they are inactive. That happens due to the interpolation error that tends to smooth out the RS topology. Perhaps it is not clear in the picture, but the reference mechanics mesh’s layers (since the thickness distribution in the $z-$ direction is not uniform but instead graded toward the edge) are not evenly-spaced which tells why these inactive spots appear.

Figure 10.

Snapshots showing pressure-drop [Psi] evolution at 10 years.

Finally, Table 3 reviews results for the minimum and maximum vertical displacements $u_z$ for all cases considered above. Notice that the differences between them are less than 3 $\%$ for $u_z$ min and 8 $\%$ for $u_z$ max, which proves the consistency of the projection operator. The shape of the compaction dome and the subsidence profile are alike as well. Notice that this is the case for linear isotropic elasticity. For non-linear elasticity or rate-independent plasticity probably one may expect more significant differences, though. The table also displays timing data, which reveals how the computational burden grows with the mesh refinement (see also the time spent to assemble the stiffness matrix in the last column of Table 2). Figure 11 zooms out the snapshot corresponding at 40 years to reveal the subsidence in the surface. The plot is exaggerated several times. It also exposes the subsidence profile on the surface in the centerline of the mechanics mesh in the most extended direction. The differences between the three cases are minimal; it seems that the profile does converge toward a mesh independent solution, which is not far from the last row on the table.

Case #	$u_z$ min	$u_z$ max	Runtime
One	−6.652	2.693	4 min, 34 s 23 ms
Two	−6.511	2.961	7 min, 53 s 84 ms
Three	−6.469	2.752	23 min, 42 s 18 ms

Table 3.

Simulations performed in example 6.2.

Figure 11.

Subsidence profiles after 40 years of evolution.

The above-coupled flow and geomechanics computation, which used the reconstructed model, confirmed that the procedure is quite useful to tackle realistic reservoir compaction and subsidence simulations [18].

6.3. Example 3: Nonlinear heat transfer: arch problem

The example addresses the interesting problem that has been investigated by several researchers [9, 10]. Its distinctive features are the two re-entrant corners. Near sharp corners, there may be singularities in the solution, which cause the spatial derivatives of the solution to become unbounded. The material properties are constant density and specific heat, and a linear isotropic thermal conductivity,

$\rho =1.0\mathrm{kg}/{\mathrm{m}}^3;C_p=1.0\frac{\mathrm{W}\cdot \mathrm{s}}{\mathrm{kg}{\cdot }^{\circ }\mathrm{K}};\kappa =\left(1+\frac{T}{{1000}^{\circ }\mathrm{K}}\right)\frac{\mathrm{W}}{\mathrm{m}{\cdot }^{\circ }\mathrm{K}}.$ E35

Figure 12 shows the domain and the mesh. The BCS are of Dirichlet type on the left- ( $T={10}^3$ ) and right-most ( $T=0$ ) sides, and insulation on all other sides: $n\cdot \left(\kappa \nabla T\right)=0$ . The triangular mesh consists of 7985 points and 15,539 elements. The domain lengths are 1 m $\times$ 0.5 m. The initial temperature distribution was taken to be [9]:

$T\left(x,y,t^{\ast }\right)={10}^3\text{erfc}{\left(\frac{x}{2\sqrt{\kappa t^{\ast }}}\right)}^{˚}\mathrm{K},$ E36

which is the short-time linear solution at a time $t^{\ast }$ for a plane semi-infinite medium. In the analysis, it is assumed $\kappa =1$ and $t^{\ast }=0.0005\mathrm{s}$ in the calculation of the initial conditions.

Figure 12.

The mesh for the arch-problem.

Figure 13 shows temperature field snapshots for different times increasing from top to bottom. The example simulates 0.1 s with a fully implicit approach. It is observed that a heating front quickly travels from left to right as expected due to the temperature gradient. The temperature scale in the color maps is from 0 to 1000°K. As a qualitative benchmark, the temperature profile reported by Winget and Hughes [9] accords very well with the results herein.

Figure 13.

Temperature, $T_h$ , (left) and mean stress, $S_m$ , snapshots (right).

The example finalizes with a simple loosely coupled thermal and mechanics computation. It takes the temperature variation that the arch problem experiences as driving force for the mechanical problem. It assumes linear isotropic elasticity with $E=30$ Ksi and $\nu =0.3$ and the coefficient of thermal dilatation $\beta =1\cdot {10}^{-5}{K}^{-1}$ and the bulk modulus. The bottommost edges are clamped while the remainders are traction free. The right column in Figure 13 includes three snapshots that depict the mean stress. Dilatation grows from the upper-right corner while compression appears from the upper-left corner, which are clearly observed in the results. The figure depicts the magnitude of the induced thermal stresses. The reader should refer to [10] for further details about this thermo-elasticity example.

7. Concluding remarks

This chapter introduced how to estimate stress-induced changes using elasticity simulations that are often performed through FE computations. It thus presented a formulation for linear thermo-poroelasticity. It covered the nonlinear energy equation as well. It also implemented a comprehensive MFEM on curved interfaces where the classical DN-DDM was employed to decouple the global SPP for elasticity, and steady single-phase flow. The coupled flow and geomechanics computation that utilizes the reconstructed model showed that this workflow is valuable to tackle realistic reservoir compaction and subsidence simulations. The research presented herein unfolds new prospects to further parallel codes for reservoir simulation coupled with geomechanics.

Acknowledgments

The author recognizes the financial support of the project “Reduced-Order Parameter Estimation for Underbody Blasts” financed by the Army Research Laboratory, through the Army High-Performance Computing Research Center under Cooperative Agreement W911NF-07-2-0027 and also acknowledgments Dr. Belsay Borges for proofreading the manuscript.