This benchmark simulates a soil column with fluid injection at the bottom and a production well at the top. It is taken from Kim [1], in detail it coincides with one of his examples (case 2, coupling strength $\tau=1.21$). A brief description of the used staggered scheme follows at the end.
Simulation model with fluid source, sink, observation point and boundary conditions
The fluid enters and leaves only via the source and sink in the domain, there is no flow across the boundaries. The displacements at the bottom are fixed, whereas there is a vertical traction applied on top. Originally the problem is onedimensional, for simulation with OpenGeoSys it is created in two dimensions with corresponding boundary conditions. All parameters are concluded in the following tables.
Property  Value  Unit 

Fluid density  $10^3$  kg/m$^3$ 
Viscosity  $10^{3}$  Pa$\cdot$s 
Fluid compressibility  $27.5\cdot 10^{9}$  Pa$^{1}$ 
Porosity  $0.3$   
Permeability  $493.5\cdot 10^{16}$  m$^2$ 
Young’s modulus (bulk)  $300\cdot 10^6$  Pa 
Poisson ratio (bulk)  $0$   
Biot coefficient  $1.0$   
Solid density  $3\cdot 10^3$  kg/m$^3$ 
Solid compressibility  $0$  Pa$^{1}$ 
Property  Value  Unit 

Height ($y$)  $150$  m 
Width ($x$)  $10$  m 
Finite Elements  $15$ TaylorHood quadrilateral elements  10 m $\times$ 10 m 
Time step  $86.4\cdot10^3$  s 
Hydraulic  Mechanical  

injection over area $10$m$\times 10$m  $+1.16\cdot 10^{4}$ kg/(m$^3$s)   
production over area $10$m$\times 10$m  $1.16\cdot 10^{4}$ kg/(m$^3$s)   
top  no flow  $\sigma_{yy}=2.125\cdot 10^6$ Pa 
left  no flow  $u_x=0$ 
right  no flow  $u_x=0$ 
bottom  no flow  $u_y=0$ 
initial state  $p(x,y)=2.125\cdot 10^6$ Pa  $u_x(x,y)=u_y(x,y)=0$ 
The gravity related terms are neglected in both: the Darcy velocity and the momentum balance equation.
Note that 100 time steps were used for the following results, whereas the provided input file is set to 1 time step (1 day = 86400 s). Kim plots his results over non dimensional time, referring to the time at which the produced fluid volume equals the pore volume of the domain (450 days).
Pressure at observation point (marked by circle) versus time (t=0…100 days) and spatial pressure distribution at t=100 days
For each time step run alternating simulations of the hydraulic (H) problem and the mechanical (M) problem until a convergence criteria is met. The fixedstress split starts with the mass balance (H) followed by the momentum balance (M). These coupling iterations (H,M,H,M,…) add another iteration level compared to the monolithic formulation (HM). However, due to splitting into smaller problems this may result in a speedup.
The weak form of mass (scalar test function $v$) and momentum balance (vectorial test function $\mathbf{v}$) after discretization in time (Eulerbackward) reads at the $n^\mathrm{th}$ time step for the field equations (boundary terms omitted) $$ \eqalign{ \int_V \left( \alpha \frac{ \varepsilon_V(\mathbf{u}^{n}) \varepsilon_V(\mathbf{u}^{n1})}{\Delta t}v +S\frac{p^np^{n1}}{\Delta t}v +\frac{k}{\mu}\nabla p^n \cdot \nabla v \right) \mathrm{d}V &=&\int_V \left(s v\frac{k}{\mu}\varrho_\mathrm{f} \nabla \cdot \mathbf{g} v \right) \mathrm{d}V, \cr \int_V \left( \boldsymbol{\sigma}^\mathrm{eff}(\mathbf{u}^n):\boldsymbol{\varepsilon}(\mathbf{v}) \alpha p^n \varepsilon_V(\mathbf{v}) \right) \mathrm{d}V &=& \int_V \ \varrho_\mathrm{b} \mathbf{g}\cdot \mathbf{v} \ \mathrm{d}V, } $$ where $\alpha$ denotes Biot coefficient, $S$ storage coefficient, $k$ permeability, $\mu$ viscosity, $\varrho_\mathrm{f}$ fluid density and $\varrho_\mathrm{b}$ bulk density. The fixedstress split makes the assumption of constant volumetric stress to eliminate the displacement from the mass balance. At the $i^\mathrm{th}$ coupling iteration this assumption leads to $$ K_\mathrm{b} \varepsilon_V(\mathbf{u}^{n,i})\alpha p^{n,i} = K_\mathrm{b} \varepsilon_V(\mathbf{u}^{n,i1})\alpha p^{n,i1}, $$ where $K_\mathrm{b}$ denotes the drained bulk modulus. Thus the weak form of the mass balance (H) depends only on pressure $p^{n,i}$ as unknown $$ \int_V \left( \beta_\mathrm{FS}\frac{p^{n,i}p^{n,i1}}{\Delta t}v +S\frac{p^{n,i}p^{n1}}{\Delta t}v +\frac{k}{\mu}\nabla p^{n,i} \cdot \nabla v \right) \mathrm{d}V =\int_V \left( s v\frac{k}{\mu}\varrho_\mathrm{f} \nabla \cdot \mathbf{g} v \alpha\frac{ \varepsilon_V(\mathbf{u}^{n,i1}) \varepsilon_V(\mathbf{u}^{n1})}{\Delta t} v \right) \mathrm{d}V, $$ where $\beta_\mathrm{FS}=\beta_\mathrm{FS}^\mathrm{ph}=\frac{\alpha^2}{K_\mathrm{b}}$ would correspond to the physically motivated volumetric stress. However, choosing another artificial volumetric stress, i.e. using a different bulk modulus than $K_\mathrm{b}$, may accelerate the convergence of the coupling iterations. The analysis by Mikelic and Wheeler [2] revealed that $\beta_\mathrm{FS}^\mathrm{MW}=\frac{\alpha^2}{2K_\mathrm{b}}$ is generally close to the apriori unknown optimal value $\beta_\mathrm{FS}^\mathrm{opt}$, for more information see user guide  conventions.
The obtained pressure $p^{n,i}$ is then inserted into the momentum equation (M) with $\mathbf{u}^{n,i}$ as unknown $$ \int_V \ \boldsymbol{\sigma}^\mathrm{eff}(\mathbf{u}^{n,i}):\boldsymbol{\varepsilon}(\mathbf{v}) \ \mathrm{d}V = \int_V \left( \varrho_\mathrm{b} \mathbf{g}\cdot \mathbf{v} + \alpha p^{n,i} \varepsilon_V(\mathbf{v}) \right) \mathrm{d}V. $$ These two steps (H,M) are repeated until convergence is reached and thus the solution of the $n^\mathrm{th}$ time step is found.
[1] Kim, J. and Tchelepi, H.A. and Juanes, R. (2009): Stability, Accuracy and Efficiency of Sequential Methods for Coupled Flow and Geomechanics. SPE International, vol. 16, p. 249262, DOI:https://doi.org/10.2118/119084PA https://onepetro.org/SJ/articleabstract/16/02/249/204235
[2] Mikelic, A. and Wheeler, M.F. (2013): Convergence of iterative coupling for coupled flow and geomechanics. Computational Geosciences, vol. 17, p. 455461, DOI:https://doi.org/10.1007/s105960129318y https://link.springer.com/content/pdf/10.1007/s105960129318y.pdf
This article was written by Dominik Kern and Wenqing Wang. If you are missing something or you find an error please let us know.
Generated with Hugo 0.101.0
in CI job 282302

Last revision: November 22, 2022
Commit: [PL/THM] Implement freezing for temperature eq. 68ebbec
 Edit this page on