6.2. Example 2: Coupled flow and mechanics

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

Points Elements Kind of mesh 1814 Triangular 1472 Quadrilateral 7858 Triangular

matched the displacements on the interface, a good accordance is also obtained for the hori-

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 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 log � log plot. The slope of the least-squares straight line is 1.44143, where the coefficient of determination is <sup>R</sup><sup>2</sup> <sup>¼</sup> <sup>84</sup>%.

> 8 256 (1 / h)

. Notice that besides the example only

O(h ^ 3/2) O(h) L^2

also displayed. The discrepancy is of the order of 10�<sup>4</sup>

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

zontal derivative.

<sup>10</sup>-5

<sup>10</sup>-6

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

<sup>10</sup>-4

L^2 error

<sup>10</sup>-3

<sup>10</sup>-2

Table 1. Meshes for example 6.1.

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 f and permeability for the isotropic case Kx ¼ Kz ¼ Ky as shown in Figure 6, whose depiction is three

Figure 5. The numerical L<sup>2</sup> convergence rate for problem 6.1.

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 ct <sup>¼</sup> <sup>1</sup>:<sup>4</sup> � <sup>10</sup>�<sup>5</sup> Psi�<sup>1</sup> (M�<sup>1</sup> <sup>¼</sup> f � ct). This example does not incorporate gravity loading for both flow and mechanics.

Table 2 compiles the mesh dimensions in every direction. The example also contemplates Nz ¼ 10, Nu ¼ 5, No ¼ 7, Nc ¼ 5 (mesh patches on the corners and No and Nu stand for overand under-burden respectively). The table also displays the number of elements, 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.

Figure 8 displays the mechanics mesh. The second case on Table 2 corresponds to a layered RS with Young's modulus Eu <sup>¼</sup> <sup>3</sup> � 104 , Ep <sup>¼</sup> <sup>1</sup> � 104 , Eo <sup>¼</sup> <sup>2</sup> � 104 ½ � Psi , 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

[ ] *n* = 0

= 0, *F* = 0 *<sup>z</sup> <sup>x</sup> u*

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

0 0 = = *z x F u*

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

0 0 = = *z x F u*

Figure 9 pictures snapshots with the evolution of the vertical displacement uz ½ � m 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 compactiondome 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

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

realistic geomechanics model with stiffer surroundings around the RS.

Figure 8. The hexahedral mesh generated for 2nd case in Table 2.

not entirely drain but is still a compelling case for coupled flow and mechanics.

Figure 6. The reservoir's permeability Ky.


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. The hexahedral mesh generated for 2nd case in Table 2.

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 ct <sup>¼</sup> <sup>1</sup>:<sup>4</sup> � <sup>10</sup>�<sup>5</sup> Psi�<sup>1</sup> (M�<sup>1</sup> <sup>¼</sup>

Table 2 compiles the mesh dimensions in every direction. The example also contemplates Nz ¼ 10, Nu ¼ 5, No ¼ 7, Nc ¼ 5 (mesh patches on the corners and No and Nu stand for overand under-burden respectively). The table also displays the number of elements, 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

Figure 8 displays the mechanics mesh. The second case on Table 2 corresponds to a layered RS

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

Case # Description Nx Ny 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

, Eo <sup>¼</sup> <sup>2</sup> � 104 ½ � Psi , while Poisson's ratio,

, Ep <sup>¼</sup> <sup>1</sup> � 104

f � ct). This example does not incorporate gravity loading for both flow and mechanics.

also assumes a zero initial displacement field.

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

with Young's modulus Eu <sup>¼</sup> <sup>3</sup> � 104

Table 2: Mesh sizes and simulations in example 6.2.

Figure 6. The reservoir's permeability Ky.

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 9 pictures snapshots with the evolution of the vertical displacement uz ½ � m 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 compactiondome 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 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

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

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

Case # uz min uz 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

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

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

Figure 12 shows the domain and the mesh. The BCS are of Dirichlet type on the left- (<sup>T</sup> <sup>¼</sup> 103

and right-most (T ¼ 0) sides, and insulation on all other sides: n � ð Þ¼ κ∇T 0. The triangular

kg�<sup>∘</sup> <sup>K</sup> ; <sup>κ</sup> <sup>¼</sup> <sup>1</sup> <sup>þ</sup>

T 1000<sup>∘</sup> K W

<sup>m</sup>�<sup>∘</sup> <sup>K</sup> : (35)

)

; Cp <sup>¼</sup> <sup>1</sup>:<sup>0</sup> <sup>W</sup> � <sup>s</sup>

subsidence simulations [18].

Table 3. Simulations performed in example 6.2.

linear isotropic thermal conductivity,

<sup>r</sup> <sup>¼</sup> <sup>1</sup>:0 kg=m3

6.3. Example 3: Nonlinear heat transfer: arch problem

Figure 11. Subsidence profiles after 40 years of evolution.

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

the thickness distribution in the zdirection is not uniform but instead graded toward the edge) are not evenly-spaced which tells why these inactive spots appear.

Finally, Table 3 reviews results for the minimum and maximum vertical displacements uz for all cases considered above. Notice that the differences between them are less than 3% for uz min and 8% for uz 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


Table 3. Simulations performed in example 6.2.

Figure 11. Subsidence profiles after 40 years of evolution.

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.

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 thickness distribution in the zdirection is not uniform but instead graded toward the

Figure 9. Snapshots at 10 and 20 years of evolution showing the vertical-displacement field uz, the pay-zone displays

pressure.

Finally, Table 3 reviews results for the minimum and maximum vertical displacements uz for all cases considered above. Notice that the differences between them are less than 3% for uz min and 8% for uz 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

edge) are not evenly-spaced which tells why these inactive spots appear.

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

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

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/m^3}; \ \mathrm{C}\_{\mathbb{P}} = 1.0 \frac{\mathrm{W} \cdot \mathrm{s}}{\mathrm{kg} \cdot \mathrm{\text{\textdegree K}}}; \ \mathrm{\textdegree} = \left(1 + \frac{T}{1000 \,\mathrm{\textdegree K}}\right) \frac{\mathrm{W}}{\mathrm{m} \cdot \mathrm{\textdegree K}}. \tag{35}$$

Figure 12 shows the domain and the mesh. The BCS are of Dirichlet type on the left- (<sup>T</sup> <sup>¼</sup> 103 ) and right-most (T ¼ 0) sides, and insulation on all other sides: n � ð Þ¼ κ∇T 0. The triangular mesh consists of 7985 points and 15,539 elements. The domain lengths are 1 m�0.5 m. The initial temperature distribution was taken to be [9]:

$$T(x, y, t^\*) = 10^3 \text{erfc}\left(\frac{x}{2\sqrt{\kappa t^\*}}\right) \text{K}\_\prime \tag{36}$$

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

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

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

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-

[1] Lewis R, Schrefler B. The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media. 2nd ed. New York: John Wiley and Sons; 1998 [2] Florez H. Domain decomposition methods for geomechanics [Ph.D. thesis]. The Univer-

[3] Florez H, Wheeler M. A mortar method based on NURBS for curved interfaces. Computer Methods in Applied Mechanics and Engineering. 2016;310:535-566. ISSN 0045-

0027 and also acknowledgments Dr. Belsay Borges for proofreading the manuscript.

refer to [10] for further details about this thermo-elasticity example.

7. Concluding remarks

coupled with geomechanics.

Acknowledgements

Author details

Address all correspondence to: florezg@gmail.com

sity of Texas at Austin; 2012

7825. DOI: 10.1016/j.cma.2016.07.030

The University of Texas at El Paso, Adelphi, Maryland, USA

Horacio Florez

References

which is the short-time linear solution at a time t <sup>∗</sup> for a plane semi-infinite medium. In the analysis, it is assumed κ ¼ 1 and t <sup>∗</sup> <sup>¼</sup> <sup>0</sup>:0005 s in the calculation of the initial conditions.

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 E ¼ 30 Ksi and ν ¼ 0:3 and the coefficient of thermal dilatation <sup>β</sup> <sup>¼</sup> <sup>1</sup> � <sup>10</sup>�<sup>5</sup>K�<sup>1</sup> 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

Figure 12. The mesh for the arch-problem.

Figure 13. Temperature, Th, (left) and mean stress, Sm, snapshots (right).

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.
