Geomechanics

In PFLOTRAN, a linear elasticity model is assumed as the constitutive model for deformation of the rock. Biot’s model is used to incorporate the effect of flow on geomechanics. In addition, the effect of temperature on geomechanics is considered via the coefficient of thermal expansion. The following governing equations are used:

(1)\[\nabla \cdot [{\boldsymbol{\sigma}}] + \density {\boldsymbol{b}} = 0 \quad \mathrm{in} \; \Omega,\]

(2)\[\begin{split}&{\boldsymbol{\sigma}} = \lambda \text{tr}\left({\boldsymbol{\varepsilon}}\right) + 2\mu {\boldsymbol{\varepsilon}} - \beta p {\boldsymbol{I}} - \alpha T {\boldsymbol{I}}, \\ &{\boldsymbol{\varepsilon}} = \frac{1}{2} \left(\nabla{\boldsymbol{u}}({\boldsymbol{x}}) + [\nabla {\boldsymbol{u}}({\boldsymbol{x}})]^{T} \right), \\ &{\boldsymbol{u}}({\boldsymbol{x}}) = {\boldsymbol{u}}^p({\boldsymbol{x}}) \quad \mathrm{on} \; \Gamma^D,\end{split}\]

(3)\[{\boldsymbol{\sigma}}{\boldsymbol{n}}({\boldsymbol{x}}) = t^p({\boldsymbol{x}}) \quad \mathrm{on} \; \Gamma^N,\]

where \({\boldsymbol{u}}\) is the unknown displacement field, \({\boldsymbol{\sigma}}\) is the Cauchy stress tensor, \(\lambda\), \(\mu\) are Lamé parameters (Young’s modulus and Poisson’s ratio can be related to these two parameters), \({\boldsymbol{b}}\) is the specific body force (which is gravity in most cases), \({\boldsymbol{n}}\) is the outward normal to the boundary \(\Gamma^N\). Also, \({\boldsymbol{u}}^p\) is the prescribed values of \({\boldsymbol{u}}\) on the Dirichlet part of the boundary \(\Gamma^D\), and \({\boldsymbol{t}}^p\) is the prescribed traction on \(\Gamma^N\). Additionally, \(\beta\) is the Biot’s coefficient, \(\alpha\) is the coefficient of thermal expansion, \(p\), \(T\) are the fluid pressure and temperature, obtained by solving subsurface flow problem. Also, \(p_0\) and \(T_0\) are set to initial pressure and temperature in the domain, \({\boldsymbol{\varepsilon}}\) is the strain tensor and \(\text{tr}\) is the trace of a second order tensor, \(\Omega\) is the domain, and \({\boldsymbol{I}}\) is the identity tensor. Note that stress is assumed positive under tension. The above equation also assumes that the resulting stresses and strains are relative to the undeformed configuration. The effect of deformation on the pore structure is accounted for via

(4)\[\porosity = \porosity_0 + \text{tr}({\boldsymbol{\varepsilon}}).\]

Note that the above equations are solved using the finite element method (Galerkin finite element) with the displacements solved for at the vertices. Since, the flow equations are solved via the finite volume method with unknowns such as pressure and temperature solved for at the cell centers, in order to transfer data from the subsurface to geomechanics grid without interpolation, the geomechanics grid is constructed such that the vertices of the geomechanics grid coincide with the cell centers of the subsurface mesh. That is, the dual mesh of the subsurface mesh is used for the geomechanics solve.

Also, the geomechanics grid must be read in as an unstructured grid. Even if one needs to work with a structured grid, the grid must be set up in the unstructured grid format.