2. Mathematical model for thermo-poroelasticity

This section discusses the governing equations for linear homogeneous isotropic thermoporoelasticity and their FE formulation. It skips details for the sake of brevity thus a more detailed treatment can be found in [1–4]. The mathematical formulation considers a bounded domain <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup>n, n <sup>¼</sup> <sup>2</sup>, 3 and its boundary is <sup>Γ</sup> <sup>¼</sup> <sup>∂</sup>Ω, and a time interval of interest �0, <sup>ℑ</sup>½. Let <sup>T</sup> <sup>h</sup> be a non-degenerate, quasi-uniform conforming partition of Ω 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 \phi^\*}{\partial t} + \nabla \cdot \left( -\frac{1}{\mu} \underline{\mathbf{K}} (\nabla p - \rho \mathbf{g} \nabla \mathbf{z}) \right) = q\_\prime \tag{1}$$

where the equation's parameters are f<sup>∗</sup>, a model specific porosity, K represents the absolute permeability tensor. The dynamic viscosity is μ, while r 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 f<sup>∗</sup> is defined by:

$$
\phi^\* = \phi^0 + \alpha \cdot \left(\nabla \cdot \underline{u} - \varepsilon\_v^0\right) + \frac{1}{M}(p - p^0),
\tag{2}
$$

where the additional parameters are accordingly α which is the Biot's constant, u represents the displacement vector, while ε<sup>0</sup> <sup>v</sup> is the initial volumetric strain. Herein M is the Biot's modulus [6], while f<sup>0</sup> and p<sup>0</sup> define for a reference or initial state. The common BCS for the pressure equation imply Neumann or no-flow namely:

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

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 rateindependent 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.

statement, for slightly incompressible single-phase flow and Darcy's law which yields:

μ

where the equation's parameters are f<sup>∗</sup>, a model specific porosity, K represents the absolute permeability tensor. The dynamic viscosity is μ, while r 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

Kð Þ ∇p � rg∇z 

> v <sup>þ</sup>

where the additional parameters are accordingly α which is the Biot's constant, u represents the

1

<sup>v</sup> is the initial volumetric strain. Herein M is the Biot's modulus [6],

¼ q, (1)

<sup>M</sup> <sup>p</sup> � <sup>p</sup><sup>0</sup> , (2)

This section discusses the governing equations for linear homogeneous isotropic thermoporoelasticity and their FE formulation. It skips details for the sake of brevity thus a more detailed treatment can be found in [1–4]. The mathematical formulation considers a bounded domain <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup>n, n <sup>¼</sup> <sup>2</sup>, 3 and its boundary is <sup>Γ</sup> <sup>¼</sup> <sup>∂</sup>Ω, and a time interval of interest �0, <sup>ℑ</sup>½. Let <sup>T</sup> <sup>h</sup> be a non-degenerate, quasi-uniform conforming partition of Ω 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

2. Mathematical model for thermo-poroelasticity

∂f<sup>∗</sup>

porosity f<sup>∗</sup> is defined by:

displacement vector, while ε<sup>0</sup>

<sup>∂</sup><sup>t</sup> <sup>þ</sup> <sup>∇</sup> � � <sup>1</sup>

<sup>f</sup><sup>∗</sup> <sup>¼</sup> <sup>f</sup><sup>0</sup> <sup>þ</sup> <sup>α</sup> � <sup>∇</sup> � <sup>u</sup> � <sup>ε</sup><sup>0</sup>

decomposition and mortar methods for non-matching interfaces is included.

224 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

$$
\nabla p \cdot \hat{\underline{n}} = 0 \text{ on } \Gamma\_{\prime} \tag{3}
$$

one should also consider an initial or reference pressure distribution in the whole domain. Sources and sinks simulate injector and producer wells, respectively. Herein <sup>b</sup><sup>n</sup> 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:

$$-\nabla \cdot \underline{\underline{\sigma}} = \underline{f} \quad \text{in} \ \Omega \ \vdots \ \Gamma = \Gamma\_D^u \cup \Gamma\_N^u$$

$$\underline{u} = \underline{0} \ \text{on} \ \Gamma\_D^u \tag{4}$$

$$\underline{t} = \underline{\underline{\sigma}} \cdot \widehat{\underline{u}} \ \text{ on} \ \Gamma\_N^u$$

where σ is the stress tensor, 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., Γ<sup>u</sup> <sup>D</sup>, and Neumann type BCS, i.e., Γ<sup>u</sup> <sup>N</sup>, where the external tractions are known or prescribed. Hooke's law combined with Biot's poroelastic theory defines σ by the following expression:

$$\underline{\underline{\sigma}} \equiv \underline{\underline{\mathcal{C}}} : \underline{\underline{\varepsilon}} - [a(p - p^0) + 3\mathsf{K}\beta(T - T^0)] \underline{\underline{\delta}} ; \underline{\underline{\mathcal{C}}} = \lambda \underline{\underline{\delta}} \otimes \underline{\delta} + 2\mathsf{G} \underline{\underline{\mathcal{C}}} \tag{5}$$

where T ¼ T xð Þ ; t is the temperature, C is the elastic moduli, β corresponds to the coefficient of thermal dilatation while K is the bulk modulus. The Kronecker delta becomes δ while λ, and G, are the Lamé constants, and I represents the fourth-order identity tensor. The strain tensor ε is given by:

$$\underline{\underline{\varepsilon}} = \nabla^s \underline{\underline{u}} = \frac{1}{2} \left[ \nabla \underline{\underline{u}} + \left( \nabla \underline{\underline{u}} \right)^T \right]. \tag{6}$$

One can derive a weak form by substituting Eq. (2) into Eq. (1) and then multiplying by a test function v∈ H<sup>1</sup> <sup>0</sup>ð Þ Ω and integrating over Ω and using the Gauss-divergence theorem, this yields:

$$\begin{split} \int\_{\Omega} \left( \frac{1}{M} \frac{\partial p}{\partial t} v + a v \nabla \cdot \underline{\dot{u}} + \frac{1}{\mu} \underline{K} \cdot \nabla p (\nabla v)^{T} \right) \cdot d\mathbf{x} &= \int\_{\Omega} q \cdot v d\mathbf{x} + \\ \int\_{\Omega} \left( \frac{\rho g}{\mu} \underline{K} \cdot \nabla z (\nabla v)^{T} \right) d\mathbf{x} + \int\_{\partial \Omega\_{N}^{r}} v \frac{1}{\mu} \underline{K} (\nabla p - \rho g \nabla z) \cdot \hat{\underline{n}}^{T} ds. \end{split} \tag{7}$$

A weak form for the equilibrium Eq. (4) can be derived in a similar way, by testing against a given virtual displacement, χ. One arrives at:

$$\int\_{\Omega} \left( \nabla \underline{\chi} \right)^{T} : \underline{\sigma} d\Omega = \int\_{\partial \Omega\_{N}^{\rm u}} \underline{\chi}^{T} \cdot \underline{\mathfrak{t}} ds + \int\_{\Omega} \underline{\chi}^{T} \cdot \underline{f} d\Omega \tag{8}$$

<sup>K</sup> � <sup>u</sup><sup>k</sup>þ<sup>1</sup> <sup>¼</sup> <sup>f</sup> <sup>u</sup> <sup>þ</sup> <sup>Q</sup> pk � <sup>p</sup><sup>0</sup> � �

S<sup>0</sup> ¼ S þ θ � Δt � H

ð

Ω

ð

∂Ω<sup>u</sup> N

<sup>K</sup>∇<sup>Π</sup> � ð Þ <sup>∇</sup><sup>Π</sup> Tdx ; <sup>ω</sup>ð Þ<sup>2</sup> <sup>¼</sup> ð Þ <sup>1</sup>; <sup>1</sup>; <sup>0</sup> <sup>T</sup>

� Πds þ

performed in [2] to demonstrate that CG can compute pressures accurately.

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

ð

Ω

Eq. (5), that renders Eq. (15) unchanged.

1

<sup>M</sup> <sup>Π</sup> � <sup>Π</sup>Tdx ; <sup>Q</sup> <sup>¼</sup>

<sup>B</sup>TC Bdx ; <sup>f</sup> <sup>u</sup> <sup>¼</sup>

K∇p � n � �

Ω

Ω

1 μ

Ω

ð

1 μ

3. Nonlinear heat transfer equation

∂Ω<sup>p</sup> N

S ¼ ð

K ¼ ð

H ¼ ð

f <sup>p</sup> ¼

S<sup>00</sup> ¼ S � ð Þ� 1 � θ Δt � H,

<sup>S</sup><sup>0</sup> � pkþ<sup>1</sup> <sup>¼</sup> <sup>S</sup><sup>00</sup> � pk <sup>þ</sup> <sup>f</sup> <sup>p</sup> � <sup>Δ</sup><sup>t</sup> � <sup>Q</sup><sup>T</sup> <sup>u</sup><sup>k</sup>þ<sup>1</sup> � uk � �

where expressions for the matrices are provided in Eq. (16) and θ is the implicitness factor that lies between 0 and 1, while Δ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

<sup>B</sup><sup>T</sup>αωð Þ <sup>n</sup> � <sup>Π</sup>dx

ð

<sup>Ψ</sup>Tf � dx

; <sup>ω</sup>ð Þ<sup>3</sup> <sup>¼</sup> ð Þ <sup>1</sup>; <sup>1</sup>; <sup>1</sup>; <sup>0</sup>; <sup>0</sup>; <sup>0</sup> <sup>T</sup>

T

Linear Thermo-Poroelasticity and Geomechanics http://dx.doi.org/10.5772/intechopen.71873

K � ∇Πð Þ ∇z

� �dx:

Ω

ð

rg μ

Ω

<sup>t</sup> � <sup>Ψ</sup>Tds <sup>þ</sup>

<sup>Π</sup>Tq � dx <sup>þ</sup>

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

(15)

227

(16)

where <sup>t</sup> <sup>¼</sup> <sup>σ</sup> � <sup>b</sup><sup>n</sup> 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{L}\_k(\mathcal{T}\_h) = \left\{ \boldsymbol{\upsilon} \in L^2(\Omega) \; : \; \forall \mathbf{e} \in \mathcal{T}\_h, \boldsymbol{\upsilon}|\_\boldsymbol{\varepsilon} \in \mathbb{P}\_k(\boldsymbol{\varepsilon}) \right\} \tag{9}$$

where Pkð Þe represents the space of polynomials of total degree less than or equal to k, Ckð Þ T <sup>h</sup> 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(\underline{\mathbf{x}}) = (\underline{\Pi}^e)^T \cdot \underline{p}^e \; ; \underline{u}\_e^h(\underline{\mathbf{x}}) = \underline{\Psi}^e \cdot \underline{\mathbf{u}}^e \tag{10}$$

where Π<sup>e</sup> and Ψ<sup>e</sup> are matrices of shape functions given by:

$$\begin{aligned} \Pi\_i^{\epsilon} &= \psi\_i^{\epsilon}(\underline{\mathbf{x}}) \\ \Psi\_{ij}^{\epsilon} &= \begin{cases} \psi\_k^{\epsilon}(\underline{\mathbf{x}}) & \text{if } j = \mathbf{j} \\ 0 & \text{otherwise} \end{cases} \\ \mathbf{j} &= \mathbf{n} \mathbf{DOF} \cdot (k - 1) + i \; ; k = 1 \dots nm \end{aligned} \tag{11}$$

here nn is the number of nodes in the given element, i ¼ 1…nn j ¼ 1…nn � n and nDOF is the number of degrees of freedom which equals the space dimension, n. Now the engineering strain <sup>b</sup><sup>ε</sup> is defined by:

$$
\hat{\underline{\mathcal{E}}} = \underline{\mathcal{B}} \cdot \underline{\mu}^{\epsilon} \; ; \; \underline{\mathcal{B}} = \underline{\mathcal{D}} \cdot \underline{\Psi}^{\epsilon} \tag{12}
$$

where D ð Þ <sup>n</sup> , <sup>n</sup> <sup>¼</sup> <sup>2</sup>, 3 is defined as:

$$\underline{\boldsymbol{D}}\_{(2)}^{T} = \begin{bmatrix} \boldsymbol{\partial}\_{\boldsymbol{x}} & \boldsymbol{0} & \boldsymbol{\partial}\_{\boldsymbol{y}} \\ \boldsymbol{0} & \boldsymbol{\partial}\_{\boldsymbol{y}} & \boldsymbol{\partial}\_{\boldsymbol{x}} \end{bmatrix}; \ \underline{\boldsymbol{D}}\_{(3)}^{T} = \begin{bmatrix} \boldsymbol{\partial}\_{\boldsymbol{x}} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{\partial}\_{\boldsymbol{y}} & \boldsymbol{\partial}\_{\boldsymbol{z}} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{\partial}\_{\boldsymbol{y}} & \boldsymbol{0} & \boldsymbol{\partial}\_{\boldsymbol{x}} & \boldsymbol{0} & \boldsymbol{\partial}\_{\boldsymbol{z}} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{\partial}\_{\boldsymbol{z}} & \boldsymbol{0} & \boldsymbol{\partial}\_{\boldsymbol{x}} & \boldsymbol{\partial}\_{\boldsymbol{y}} \end{bmatrix}. \tag{13}$$

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:

$$
\begin{bmatrix} 0 & 0 \\ \underline{\mathbf{Q}}^{T} & \underline{\mathbf{S}} \end{bmatrix} \frac{d}{dt} \begin{Bmatrix} \underline{u} \\ \underline{p} \end{Bmatrix} + \begin{bmatrix} \underline{\mathbf{K}} & \underline{\mathbf{-Q}} \\ \underline{0} & \underline{\mathbf{H}} \end{bmatrix} \begin{Bmatrix} \underline{u} \\ \underline{p} \end{Bmatrix} = \begin{Bmatrix} \underline{f\_{u}} \\ \underline{f\_{p}} \end{Bmatrix}.\tag{14}
$$

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{aligned} \underline{\underline{\mathbf{K}}} \cdot \underline{\underline{u}}^{k+1} &= f\_{\underline{u}} + \underline{\underline{\mathbf{Q}}} \left( \underline{p}^{k} - \underline{p}^{0} \right) \\ \underline{\underline{\mathbf{S}'}} \cdot \underline{p}^{k+1} &= \underline{\underline{\mathbf{S}''}} \cdot \underline{p}^{k} + f\_{\underline{p}} \cdot \Delta t - \underline{\underline{\mathbf{Q}}}^{T} (\underline{\mu}^{k+1} - \underline{\mu}^{k}) \\ \underline{\underline{\mathbf{S}'}} &= \underline{\underline{\mathbf{S}}} + \theta \cdot \Delta t \cdot \underline{\underline{\mathbf{H}}} \\ \underline{\underline{\mathbf{S}''}} &= \underline{\underline{\mathbf{S}}} - (1 - \theta) \cdot \Delta t \cdot \underline{\underline{\mathbf{H}}} \end{aligned} \tag{15}$$

where expressions for the matrices are provided in Eq. (16) and θ is the implicitness factor that lies between 0 and 1, while Δ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{aligned} \underline{\mathbf{S}} & \equiv \int\_{\Omega} \underline{\mathbf{I}} \underline{\Pi} \cdot \underline{\Pi}^{T} d\mathbf{x} \; : \; \underline{\mathbf{Q}} = \int\_{\Omega} \underline{\underline{B}}{}^{T} d\underline{\omega}\_{(\mathbf{u})} \cdot \underline{\Pi} d\mathbf{x} \\ \underline{\mathbf{K}} & \equiv \int\_{\Omega} \underline{\underline{B}}{}^{T} \underline{\underline{C}} \; \underline{\underline{B}}{}^{T} \; : \; \underline{\underline{f}}\_{u} = \int\_{\partial \Omega\_{N}^{n}} \underline{\underline{V}}{}^{T} d\mathbf{s} + \int\_{\Omega} \underline{\underline{W}}{}^{T} \underline{\underline{f}}{}^{\cdot} d\mathbf{x} \\ \underline{\underline{H}} & \equiv \int\_{\Omega} \underline{\underline{K}}{}^{T} \underline{\underline{V}}{}^{T} d\mathbf{x} \; : \; \underline{\underline{\omega}}{}^{(2)} \; : \; \underline{\underline{\omega}}{}^{(3)} = (1, 1, 0)^{T} \; : \; \underline{\underline{\omega}}{}^{(3)} = (1, 1, 1, 0, 0, 0)^{T} \end{aligned} \tag{16}$$
 
$$\begin{aligned} \underline{f}\_{p} &= \int\_{\Omega \underline{f}\_{q}} \left( \underline{\underline{K}}{} \underline{\underline{K}} \underline{\nabla p} \cdot \underline{\underline{n}} \right) \cdot \underline{\Pi} ds + \int\_{\Omega} \underline{\Pi}{}^{T} q \cdot \mathbf{dx} + \int\_{\Omega} \left( \frac{\rho \underline{g}}{\mu} \underline{K} \cdot \underline{\nabla \underline{\Pi}}{}(\nabla z)^{T} \right) d\mathbf{x} . \end{aligned}$$

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

ð

∇χ � �<sup>T</sup>

226 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

<sup>C</sup>kð Þ¼ <sup>T</sup> <sup>h</sup> <sup>v</sup><sup>∈</sup> <sup>L</sup><sup>2</sup>

ph

where Π<sup>e</sup> and Ψ<sup>e</sup> are matrices of shape functions given by:

Πi <sup>e</sup> <sup>¼</sup> <sup>ψ</sup><sup>i</sup> e ð Þx

Ψij e ¼ : σdΩ ¼

ð

<sup>χ</sup><sup>T</sup> � tds <sup>þ</sup>

ð

<sup>χ</sup><sup>T</sup> � f d<sup>Ω</sup> (8)

Ω

<sup>e</sup> <sup>∈</sup>Pkð Þ<sup>e</sup> � � (9)

<sup>b</sup><sup>ε</sup> <sup>¼</sup> <sup>B</sup> � ue ; <sup>B</sup> <sup>¼</sup> <sup>D</sup> � <sup>Ψ</sup><sup>e</sup> (12)

3 7

<sup>5</sup>: (13)

: (14)

∂<sup>x</sup> 0 0 ∂<sup>y</sup> ∂<sup>z</sup> 0 0 ∂<sup>y</sup> 0 ∂<sup>x</sup> 0 ∂<sup>z</sup> 0 0 ∂<sup>z</sup> 0 ∂<sup>x</sup> ∂<sup>y</sup>

> p ( )

¼

f u f p ( )

<sup>e</sup> ð Þ¼ <sup>x</sup> <sup>Ψ</sup><sup>e</sup> � ue (10)

(11)

∂Ω<sup>u</sup> N

where <sup>t</sup> <sup>¼</sup> <sup>σ</sup> � <sup>b</sup><sup>n</sup> 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

where Pkð Þe represents the space of polynomials of total degree less than or equal to k, Ckð Þ T <sup>h</sup> 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

> ð Þx if j ¼ j 0 otherwise

j ¼ nDOF � ð Þþ k � 1 i ;k ¼ 1…nn

2 6 4

Finally substituting the generalized Hooke's law Eq. (5) into Eq. (8) and using Eq. (7) leads to

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

K �Q 0 H " # u

here nn is the number of nodes in the given element, i ¼ 1…nn j ¼ 1…nn � n and nDOF is the number of degrees of freedom which equals the space dimension, n. Now the engineering

<sup>e</sup> ð Þ¼ <sup>x</sup> <sup>Π</sup><sup>e</sup> ð Þ<sup>T</sup> � pe ; <sup>u</sup><sup>h</sup>

ψk e

; D<sup>T</sup> ð Þ<sup>3</sup> <sup>¼</sup>

þ

(

ð Þ Ω : ∀e∈ T <sup>h</sup>; vj

Ω

multiplied by shape or interpolation functions [8]:

spaces [7], thus:

strain <sup>b</sup><sup>ε</sup> is defined by:

ð Þ <sup>n</sup> , <sup>n</sup> <sup>¼</sup> <sup>2</sup>, 3 is defined as:

ð Þ<sup>2</sup> <sup>¼</sup> <sup>∂</sup><sup>x</sup> <sup>0</sup> <sup>∂</sup><sup>y</sup> 0 ∂<sup>y</sup> ∂<sup>x</sup> � �

the FE model for linear isotropic poroelasticity, thus:

0 0 Q<sup>T</sup> S " #

d dt

u p ( )

step. Next, one updates the pressures by using the newest displacements:

D<sup>T</sup>

where D

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

$$\rho \mathbf{C}\_p \frac{\partial T}{\partial t} = \nabla \cdot (\kappa \nabla T) + Q\_T \quad \text{on} \quad \Omega \times ]0, \mathfrak{F}[. $$

$$T = \text{g} \quad \text{on} \quad \Gamma\_D^T \times ]0, \mathfrak{F}[. $$

$$m \cdot (\kappa \nabla T) = h \quad \text{on} \quad \Gamma\_N^T \times ]0, \mathfrak{F}[. $$

$$T(x, 0) = T\_0(x) \forall x \in \Omega.$$

<sup>δ</sup><sup>K</sup> <sup>¼</sup> <sup>X</sup> p

Ωe

where described in Section 2. The reader may refer to [11] for a full treatment.

aψ<sup>j</sup> ∇ψ<sup>i</sup>

One often employs the Newton-Raphson algorithm to solve the linearized system of equations in every time step, namely, <sup>J</sup> � <sup>Δ</sup>Tð Þ<sup>ℓ</sup> ¼ �R. One can utilize the same continuous FE space that

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-

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 Γ. After one provides the initial guess on the primary variables on Γ, i.e., γ<sup>k</sup> must be given, then one can solve the problem on the white subdomains (Dirichlet problems), which corre-

2Þ

Let the primary variable be called "displacements" and their gradient "tractions," i.e., normal derivative in the boundary. Then, the tractions on the interface Γ 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 Γ, one can solve for unknown displacements to provide feedback

8 >><

>>:

Lu<sup>2</sup>

u2

∂nu<sup>2</sup>

<sup>k</sup>þ<sup>1</sup> <sup>¼</sup> <sup>f</sup> in <sup>Ω</sup><sup>2</sup>

<sup>k</sup>þ<sup>1</sup> <sup>¼</sup> 0 on <sup>∂</sup>Ω<sup>2</sup> ∩ ∂<sup>Ω</sup>

(25)

<sup>k</sup>þ<sup>1</sup> <sup>¼</sup> <sup>κ</sup><sup>k</sup>þ<sup>1</sup> on <sup>Γ</sup>

∂Kip <sup>∂</sup>Tð Þ<sup>ℓ</sup> j ¼ ð

where the variation term is given by:

4. Domain decomposition methods

Neumann (DN) DDM and Neumann-Neumann.

sponds to step 1 in Eq. (25).

1Þ

8 >>><

>>>:

Lu<sup>1</sup>

u1

u1

<sup>k</sup>þ<sup>1</sup> <sup>¼</sup> <sup>f</sup> in <sup>Ω</sup><sup>1</sup>

<sup>k</sup>þ<sup>1</sup> <sup>¼</sup> <sup>γ</sup><sup>k</sup> on <sup>Γ</sup>

<sup>k</sup>þ<sup>1</sup> <sup>¼</sup> 0 on <sup>∂</sup>Ω<sup>1</sup> ∩ ∂<sup>Ω</sup>

∂Kip <sup>∂</sup>Tð Þ<sup>ℓ</sup> j

� <sup>T</sup>ð Þ<sup>ℓ</sup>

<sup>p</sup> , (23)

Linear Thermo-Poroelasticity and Geomechanics http://dx.doi.org/10.5772/intechopen.71873 229

� �<sup>T</sup> � <sup>∇</sup>ψ<sup>p</sup> � dx: (24)

In (17), Cp is the heat capacity to constant pressure and κ ¼ κð Þ T is the thermal conductivity. QT 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(T, \upsilon) + k(T, \upsilon) - q(Q\_T, \upsilon) - f(h, \upsilon) = 0,\tag{18}$$

where the functions are:

$$\begin{aligned} m(T, v) &= \int\_{\varOmega} v \rho C\_p \mathbf{\hat{o}}\_l T \cdot d\mathbf{x}, \\ k(T, v) &= \int\_{\varOmega} \kappa (\nabla T)^T \cdot \nabla v \cdot d\mathbf{x}, \\ q(Q\_T, v) &= \int\_{\varOmega} v \cdot Q\_T \cdot d\mathbf{x} \; ; \; f(h, v) = \int\_{\varGamma\_h^\*} v \cdot h \cdot ds. \end{aligned} \tag{19}$$

Time discretization renders the local residual for the element e:

$$\begin{aligned} \underline{R} \equiv \underline{M} \cdot \left( \underline{T}^{(\ell)} - \underline{T}^{(m)} \right) + \Delta t \cdot \left( \underline{\underline{K}} \cdot \underline{T} \right)^{(m+\theta)} \\ - \Delta t \cdot \underline{q}^{(m+\theta)} - \Delta t \cdot \underline{f}^{(m+\theta)} = \underline{0} \end{aligned} \tag{20}$$

where the linear operator ð Þ� ð Þ <sup>m</sup>þ<sup>θ</sup> � ð Þ� <sup>1</sup> � <sup>θ</sup> ð Þ <sup>t</sup>¼<sup>t</sup> ð Þ <sup>m</sup> ð Þ <sup>þ</sup> <sup>θ</sup>ð Þ� <sup>t</sup>¼<sup>t</sup> ð Þ<sup>ℓ</sup> ð Þ, <sup>ℓ</sup> <sup>¼</sup> ð Þ <sup>m</sup> <sup>þ</sup> <sup>1</sup> , <sup>M</sup> <sup>K</sup> are the mass and stiffness matrix respectively. Thus the Jacobian is given by:

$$\underline{I} = \frac{\partial \underline{R}}{\partial \underline{L}^{(\ell)}} = \underline{\underline{M}} + \frac{\partial}{\partial \underline{L}^{(\ell)}} \left( \underline{\underline{K}} \cdot \underline{T} \right)^{\binom{\ell \ell}{}} \tag{21}$$

this equation renders once again:

$$\underline{J} = \underline{M} + \Delta t \cdot \left(\underline{\underline{K}} + \delta \underline{\underline{K}}\right) \tag{22}$$

if one assumes that κð Þ¼ T ð Þ a � T þ b ; a, b∈ R, then:

$$\delta \underbrace{\delta \mathcal{K}}\_{\underline{p}} = \sum\_{\underline{p}} \frac{\partial \mathcal{K}\_{\dot{\mathcal{W}}}}{\partial T\_{\underline{j}}^{(\ell)}} \cdot T\_{\underline{p}}^{(\ell)} \,. \tag{23}$$

where the variation term is given by:

rCp ∂T

228 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

no-flux BCS; h ¼ 0 of most domain boundaries.

m Tð Þ¼ ; v

k Tð Þ¼ ; v

q QT ð Þ¼ ; v

where the linear operator ð Þ� ð Þ <sup>m</sup>þ<sup>θ</sup> � ð Þ� <sup>1</sup> � <sup>θ</sup> ð Þ <sup>t</sup>¼<sup>t</sup>

if one assumes that κð Þ¼ T ð Þ a � T þ b ; a, b∈ R, then:

this equation renders once again:

Time discretization renders the local residual for the element e:

ð

vrCp∂tT � dx,

<sup>κ</sup>ð Þ <sup>∇</sup><sup>T</sup> <sup>T</sup> � <sup>∇</sup><sup>v</sup> � dx,

v � QT � dx ; f hð Þ¼ ; v

Ωe

ð

Ωe

ð

Ωe

<sup>R</sup> � <sup>M</sup> � <sup>T</sup>ð Þ<sup>ℓ</sup> � <sup>T</sup>ð Þ <sup>m</sup> � �

mass and stiffness matrix respectively. Thus the Jacobian is given by:

<sup>J</sup> <sup>¼</sup> <sup>∂</sup><sup>R</sup>

� <sup>Δ</sup><sup>t</sup> � <sup>q</sup>ð Þ <sup>m</sup>þ<sup>θ</sup> � <sup>Δ</sup><sup>t</sup> � <sup>f</sup>

<sup>∂</sup>Tð Þ<sup>ℓ</sup> <sup>¼</sup> <sup>M</sup> <sup>þ</sup>

bilinear form:

where the functions are:

<sup>T</sup> <sup>¼</sup> g on <sup>Γ</sup><sup>T</sup>

T xð Þ¼ ; 0 T0ð Þx ∀x∈ Ω:

<sup>n</sup> � ð Þ¼ <sup>κ</sup>∇<sup>T</sup> h on <sup>Γ</sup><sup>T</sup>

<sup>∂</sup><sup>t</sup> <sup>¼</sup> <sup>∇</sup> � ð Þþ <sup>κ</sup>∇<sup>T</sup> QT on <sup>Ω</sup>��0, <sup>ℑ</sup>½,

<sup>D</sup>��0, ℑ½,

(17)

(19)

(20)

(22)

<sup>N</sup>��0, ℑ�,

In (17), Cp is the heat capacity to constant pressure and κ ¼ κð Þ T is the thermal conductivity. QT represents heat sources. Neumann BCS imply heat transfer via Fourier's law: adiabatic or

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

m Tð Þþ ; v k Tð Þ� ; v q QT ð Þ� ; v f hð Þ¼ ; v 0, (18)

ð

v � h � ds:

ð Þ<sup>ℓ</sup> ð Þ, <sup>ℓ</sup> <sup>¼</sup> ð Þ <sup>m</sup> <sup>þ</sup> <sup>1</sup> , <sup>M</sup> <sup>K</sup> are the

ð Þ<sup>ℓ</sup> ð Þ (21)

Γe h

� �ð Þ <sup>m</sup>þ<sup>θ</sup>

þ Δt � K � T

ð Þ <sup>m</sup> ð Þ <sup>þ</sup> <sup>θ</sup>ð Þ� <sup>t</sup>¼<sup>t</sup>

� �

∂ <sup>∂</sup>Tð Þ<sup>ℓ</sup> <sup>K</sup> � <sup>T</sup> � � <sup>t</sup>

J ¼ M þ Δt � K þ δK

ð Þ <sup>m</sup>þ<sup>θ</sup> <sup>¼</sup> <sup>0</sup>,

$$\frac{\partial \mathcal{K}\_{ip}}{\partial T\_j^{(\ell)}} = \int\_{\Omega'} a \boldsymbol{\psi}\_j \left(\boldsymbol{\nabla} \boldsymbol{\psi}\_i\right)^T \cdot \boldsymbol{\nabla} \boldsymbol{\psi}\_p \cdot d\mathbf{x}.\tag{24}$$

One often employs the Newton-Raphson algorithm to solve the linearized system of equations in every time step, namely, <sup>J</sup> � <sup>Δ</sup>Tð Þ<sup>ℓ</sup> ¼ �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.
