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 singlephase 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 threedimensional problems of practical interest in thermoporoelasticity. 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., DirichletNeumann domain decomposition and mortar methods for nonmatching 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 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 highperformance computing into the picture. For instance, recently the author showed that the DirichletNeumann scheme could handle problems at the reservoir fieldlevel as well as the mortar method decoupled by this last one. Its current results are backed up by papers published in peerreviewed journals and conferences thus this book chapter summarizes that effort.
2. Mathematical model for thermoporoelasticity
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, 2, 3, 4]. The mathematical formulation considers a bounded domain
where the equation’s parameters are
where the additional parameters are accordingly
one should also consider an initial or reference pressure distribution in the whole domain. Sources and sinks simulate injector and producer wells, respectively. Herein
where
where
One can derive a weak form by substituting Eq. (2) into Eq. (1) and then multiplying by a test function
A weak form for the equilibrium Eq. (4) can be derived in a similar way, by testing against a given virtual displacement,
where
where
where
here
where
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:
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:
(15) 
where expressions for the matrices are provided in Eq. (16) and
(16) 
This section completes with a comment about the Continuous Galerkin (CG) formulation for the pressure (1). It is wellknown 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 multiphase flow, though. Nevertheless, one can utilize postprocessing 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 postprocessing. 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]:
In (17),
One can derive a FE formulation for model problem (17) by multiplying by a test function and integrate by parts and applying the Gaussdivergence theorem to arrive at the following bilinear form:
where the functions are:
Time discretization renders the local residual for the element
where the linear operator
this equation renders once again:
if one assumes that
where the variation term is given by:
One often employs the NewtonRaphson algorithm to solve the linearized system of equations in every time step, namely,
4. Domain decomposition methods
Domain Decomposition Methods (DDM) encompass highly efficient algorithms to obtain solutions of largescale discrete problems on parallel supercomputers. 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 DirichletNeumann (DN) DDM and NeumannNeumann.
Let
Let the primary variable be called “displacements” and their gradient “tractions,” i.e., normal derivative in the boundary. Then, the tractions on the interface
It happens that this approach requires at least a twoentry 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
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 nonmatching interfaces. The section first introduces a brief description of nonuniform rational BSplines 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,
where
the parameters in Eq. (28) are as follows:
herein
One can write in a matrix or algebraic form, Eq. (28) as:
The equation above corresponds to the socalled “saddlepoint problem (SPP).” Notice that subdomains are only connected using the Lagrange multiplier
The following line integral defines the projector, for 2D problems, as:
where
where
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,
6.1. Example 1: Twodimensional steady state singlephase 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:
where the domain of interest corresponds to the unitary square and Dirichlet BCS are enforced. The input pressure field is given by:
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 topleft corner of the figure shows the mesh that is employed.
The pressure field is on the righttop corner, and its horizontal derivative is in the bottomleft corner, while the discrepancy between the numerical and exact solutions, i.e., the absolute error, was rendered in the rightbottom 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 BSplines interpolants (NURBS with all weights equal 1) that were constructed by interpolating a sinusoidal wave as the figure shows (see [3] for details). Thirtytwo 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
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., multithreading 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 nonmortar sides. It tackled meshes of size 8, 16, 32, 64, 128 and 256 respectively. Figure 4 displays the resulting convergence rate in a
Finally, Figure 5 shows pressure snapshots that represent four different DirichletNeumann iteration levels evolving from lefttoright and toptobottom. 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
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 socalled 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 payzone 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 payzone such as porosity
Table 2 compiles the mesh dimensions in every direction. The example also contemplates
Case #  Description  Nx  Ny 
 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 
Figure 8 displays the mechanics mesh. The second case on Table 2 corresponds to a layered RS with Young’s modulus
Figure 9 pictures snapshots with the evolution of the vertical displacement
Figure 10 renders pressuredrop 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 toptobottom and lefttoright. 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 redcolored 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
Finally, Table 3 reviews results for the minimum and maximum vertical displacements
Case # 

 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 
The abovecoupled 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 reentrant 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,
Figure 12 shows the domain and the mesh. The BCS are of Dirichlet type on the left (
which is the shorttime linear solution at a time
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.
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
7. Concluding remarks
This chapter introduced how to estimate stressinduced changes using elasticity simulations that are often performed through FE computations. It thus presented a formulation for linear thermoporoelasticity. It covered the nonlinear energy equation as well. It also implemented a comprehensive MFEM on curved interfaces where the classical DNDDM was employed to decouple the global SPP for elasticity, and steady singlephase 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.