Utah FORGE Native State Model

Governing Equations

The native state model considers geomechanics, fluid flow, and heat transfer. Instead of using Actions, such as the PorousFlowUnsaturated Actions in the previous RTES examples, we now build the THM coupled governing equation using Kernals from the scratch. First, we define the primary unknows, which are the pore pressure, temperature, and displacement as

[Variables]
  [./pressure]
  [../]
  [./temperature]
    scaling = 1E-8
  [../]
  [./disp_i]
    scaling = 1E-10
  [../]
  [./disp_j]
    scaling = 1E-10
  [../]
  [./disp_k]
    scaling = 1E-10
  [../]
[]

Note the scaling factors are used to achieve a better convergence rate, see the instruction for details.

The complete governing equations of the THM coupled problem are detialed in the theory page. For this FORGE simulate case, we simplify them as where and denote pore presure and temperature, while and are the stress and identity tensors. is the Biot’s constant, and is the thermal expansion coefficient for rocks, represents the density of porous media including both rock skeleton and filled fluids, and the gravity vector. For the mass convervatiom equation, where is the displacement vector of the solid phase, is the permeability tensor of porous media, is the fluid flow source, and is the Biot’s modulus. For the energy conservation equation, is the thermal conductivity tensor for whole rock media, the specific enthalpy for phase j, is the advective Darcy flux, and is the heat source.

Note the objective of this simulation is to obtain the native state of the FORGE site with a steady state simulation. The pressure () and temperature () time derivative terms, in above equations are omitted. Given the extreme low permeability, the heat conduction term due to fluid flow is also omitted. We also assume the fluid in FORGE site is in the liquid state under the native state. With above reasoning, the Kernals section for the governing equation is listed in the following:

[Kernels]
  [./Darcy_flow]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pressure
    gravity = '0 0 -9.8'
  [../]
  [./heat_conduction]
    type = PorousFlowHeatConduction
    variable = temperature
  [../]
  [./grad_stress_i]
    type = StressDivergenceTensors
    variable = disp_i
    component = 0
  [../]
  [./grad_stress_j]
    type = StressDivergenceTensors
    variable = disp_j
    component = 1
  [../]
  [./grad_stress_k]
    type = StressDivergenceTensors
    variable = disp_k
    component = 2
   # displacements = 'disp_i disp_j disp_k'
  [../]
  [./poro_i]
    type = PoroMechanicsCoupling
    variable = disp_i
    porepressure = pressure
    component = 0
  [../]
  [./poro_j]
    type = PoroMechanicsCoupling
    variable = disp_j
    porepressure = pressure
    component = 1
  [../]
  [./poro_k]
    type = PoroMechanicsCoupling
    variable = disp_k
    porepressure = pressure
    component = 2
  [../]
  [./weight]
    type = Gravity
    variable = disp_k
    value = -9.8
  [../]
[]

In order to obtain the secondary variables, such as the stress tensor, fluid density, viscosity, etc, the AuxKernels section needs to explicitly request these output.

Initial and boundary conditions

For fluid flow field, as shown in Figure 2, we applied

• 101 atmospheric pore pressure on the top surface,

• 34 pore pressure on the bottom surfaced,

• and Zero flux for all side surfaces

Figure 2: Illustration of boundary conditions for coupled three field.

The input lines are

[./pore_pressure_top]
    type = DirichletBC
    variable = pressure
    boundary = 'surface_40m'
    value = 101325 #atmospheric pressure in Pa-assumes water level at ground surface
  [../]
  [./pore_pressure_bottom]
    type = DirichletBC
    variable = pressure
    boundary = 'bottom_40m'
    value = 34000000
[../]

[./Side_Z_Minus_T]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'bottom_40m'
    function = T_on_face_Z_Minus
  [../]
  [./Side_Z_Plus_T]
    type = DirichletBC
    variable = temperature
    boundary = 'surface_40m'
    value= 299
  [../]
[Functions] # temperature on boundary surfaces
  [./T_on_face_Z_Minus]
    type = PiecewiseBilinear0
    data_file = bottom_2021_11_10_plus_50.csv
    xaxis = 1
    yaxis = 0
  [../]
[]

    [./roller_sigma_h_min]
    type = DirichletBC
    preset = true
    variable = disp_i
    value = 0
    boundary = 'bottom_40m_granitoid_40m_front granitoid_40m_surface_40m_front'
  [../]
  [./roller_sigma_h_max]
    type = DirichletBC
    preset = true
    variable = disp_j
    value = 0
    boundary = 'bottom_40m_granitoid_40m_right granitoid_40m_surface_40m_right'
  [../]
  [./bottom_z]
    type = DirichletBC
    preset = true
    variable = disp_k
    value = 0
    boundary = 'bottom_40m'
  [../]
    [./S_h_min_normal]
    type = FunctionNeumannBC
    variable = disp_i
    boundary = 'bottom_40m_granitoid_40m_back granitoid_40m_surface_40m_back'
    function  = '-15709*(1785-z)'
  [../]
  [./Side_Y_Plus_normal]
    type = FunctionNeumannBC
    variable = disp_j
    boundary = 'bottom_40m_granitoid_40m_left granitoid_40m_surface_40m_left'
    function  = '-19640*(1785-z)'
  [../]
  [./Side_Z_Plus_traction_Z]
    type = NeumannBC
    variable = disp_k
    boundary = 'surface_40m'
    value = -101325
  [../]

Initial condition

We only initialized the pore pressure and temperature approximately using the same z depenedent functions as

[ICs]
  [temperature]
    type = FunctionIC
    variable = temperature
    function = '616 - ((2360+z)*7.615e-2)'
  []
  [pressure]
    type = FunctionIC
    variable = pressure
    function = '3.5317e7 - ((2360+z)*8455)'
  []
[]