PorousFlowBasicTHM

This action allows simple simulation of fully-saturated, single-phase, single-component fluid flow. The fluid flow may be optionally coupled to mechanics and/or heat flow using the coupling_type flag.

The fluid equation is a simplified form of the full PorousFlow fluid equation (see PorousFlowFullySaturatedMassTimeDerivative and PorousFlowFullySaturatedDarcyBase for the derivation): (1) Note that the fluid-mass time derivative is close to linear, and is perfectly linear if multiply_by_density=false, and this also almost linearises the flow term. Extremely good nonlinear convergence should therefore be expected, but there are some knock-on effects that are documented in PorousFlowFullySaturatedMassTimeDerivative.

In fully-saturated, single-phase simulations upwinding is typically unnecessary. Moreover, the standard PorousFlow Kernels are somewhat inefficient in the fully-saturated case since Material Properties such as relative permeabilities and saturations actually do not need to be computed.

In fully-saturated, single-phase, single-component simulations, the mass lumping is also typically unncessary. More importantly, in many real-life situations, as may be seen in Eq. (1) the time derivative of the fluid mass may be linearised, which leads to improved convergence.

To simulate Eq. (1) the PorousFlowBasicTHM Action employs the following Kernels:

For isothermal simulations, the fluid properties may still depend on temperature, so the temperature input parameter may be set to any real number, or to an AuxVariable if desired.

For anisothermal simulations, the energy equation reads (2) where the final term is only used if coupling with mechanics is also desired. To simulate this DE, PorousFlowBasicTHM uses the following kernels:

For simulations that couple fluid flow to mechanics, the equations are already written in governing equations, and PorousFlowBasicTHM implements these by using the following kernels:

PorousFlowBasicTHM adds many Materials automatically, however, to run a simulation you will need to provide more Materials for each mesh block, depending on your simulation type, viz:

note

Since upwinding and no mass lumping of the fluid mass are used (for simplicity, efficiency and to reduce numerical diffusion), the results may be slightly different to simulations that employ upwinding and mass lumping.

A simple example of using PorousFlowBasicTHM is documented in the PorousFlow tutorial with input file:

# Darcy flow
[Mesh]
  type = AnnularMesh
  dim = 2
  nr = 10
  rmin = 1.0
  rmax = 10
  growth_r = 1.4
  nt = 4
  tmin = 0
  tmax = 1.57079632679
[]

[MeshModifiers]
  [./make3D]
    type = MeshExtruder
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  [../]
  [./shift_down]
    type = Transform
    transform = TRANSLATE
    vector_value = '0 0 -6'
    depends_on = make3D
  [../]
  [./aquifer]
    type = SubdomainBoundingBox
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    depends_on = shift_down
  [../]
  [./injection_area]
    type = ParsedAddSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    depends_on = 'aquifer'
  [../]
  [./rename]
    type = RenameBlock
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    depends_on = 'injection_area'
  [../]
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [./porepressure]
  [../]
[]

[PorousFlowBasicTHM]
  porepressure = porepressure
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = the_simple_fluid
[]

[BCs]
  [./constant_injection_porepressure]
    type = PresetBC
    variable = porepressure
    value = 1E6
    boundary = injection_area
  [../]
[]

[Modules]
  [./FluidProperties]
    [./the_simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 2E9
      viscosity = 1.0E-3
      density0 = 1000.0
    [../]
  [../]
[]

[Materials]
  [./porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  [../]
  [./biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.8
    solid_bulk_compliance = 2E-7
    fluid_bulk_modulus = 1E7
  [../]
  [./permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  [../]
  [./permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  [../]
[]

[Preconditioning]
  active = basic
  [./basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  [../]
  [./preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]

[Outputs]
  exodus = true
[]
(modules/porous_flow/examples/tutorial/01.i)

A TH example that uses PorousFlowBasicTHM is also documented in the PorousFlow tutorial with input file:

# Darcy flow with heat advection and conduction
[Mesh]
  type = AnnularMesh
  dim = 2
  nr = 10
  rmin = 1.0
  rmax = 10
  growth_r = 1.4
  nt = 4
  tmin = 0
  tmax = 1.57079632679
[]

[MeshModifiers]
  [./make3D]
    type = MeshExtruder
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  [../]
  [./shift_down]
    type = Transform
    transform = TRANSLATE
    vector_value = '0 0 -6'
    depends_on = make3D
  [../]
  [./aquifer]
    type = SubdomainBoundingBox
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    depends_on = shift_down
  [../]
  [./injection_area]
    type = ParsedAddSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    depends_on = 'aquifer'
  [../]
  [./rename]
    type = RenameBlock
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    depends_on = 'injection_area'
  [../]
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [./porepressure]
  [../]
  [./temperature]
    initial_condition = 293
    scaling = 1E-8
  [../]
[]

[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydro
  gravity = '0 0 0'
  fp = the_simple_fluid
[]

[BCs]
  [./constant_injection_porepressure]
    type = PresetBC
    variable = porepressure
    value = 1E6
    boundary = injection_area
  [../]
  [./constant_injection_temperature]
    type = PresetBC
    variable = temperature
    value = 313
    boundary = injection_area
  [../]
[]

[Modules]
  [./FluidProperties]
    [./the_simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 2E9
      viscosity = 1.0E-3
      density0 = 1000.0
      thermal_expansion = 0.0002
      cp = 4194
      cv = 4186
      porepressure_coefficient = 0
    [../]
  [../]
[]

[Materials]
  [./porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  [../]
  [./biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.8
    solid_bulk_compliance = 2E-7
    fluid_bulk_modulus = 1E7
  [../]
  [./permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  [../]
  [./permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  [../]

  [./thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    biot_coefficient = 0.8
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  [../]
  [./rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  [../]
  [./thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caps aquifer'
  [../]
[]

[Preconditioning]
  active = basic
  [./basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  [../]
  [./preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]

[Outputs]
  exodus = true
[]
(modules/porous_flow/examples/tutorial/03.i)

A THM example that uses PorousFlowBasicTHM is also documented in the PorousFlow tutorial with input file:

# Darcy flow with heat advection and conduction, and elasticity
[Mesh]
  type = AnnularMesh
  dim = 2
  nr = 10
  rmin = 1.0
  rmax = 10
  growth_r = 1.4
  nt = 4
  tmin = 0
  tmax = 1.57079632679
[]

[MeshModifiers]
  [./make3D]
    type = MeshExtruder
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  [../]
  [./shift_down]
    type = Transform
    transform = TRANSLATE
    vector_value = '0 0 -6'
    depends_on = make3D
  [../]
  [./aquifer]
    type = SubdomainBoundingBox
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    depends_on = shift_down
  [../]
  [./injection_area]
    type = ParsedAddSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    depends_on = 'aquifer'
  [../]
  [./rename]
    type = RenameBlock
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    depends_on = 'injection_area'
  [../]
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  PorousFlowDictator = dictator
  biot_coefficient = 1.0
[]

[Variables]
  [./porepressure]
  [../]
  [./temperature]
    initial_condition = 293
    scaling = 1E-8
  [../]
  [./disp_x]
    scaling = 1E-10
  [../]
  [./disp_y]
    scaling = 1E-10
  [../]
  [./disp_z]
    scaling = 1E-10
  [../]
[]

[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydroMechanical
  gravity = '0 0 0'
  fp = the_simple_fluid
  thermal_eigenstrain_name = thermal_contribution
  use_displaced_mesh = false
[]

[BCs]
  [./constant_injection_porepressure]
    type = PresetBC
    variable = porepressure
    value = 1E6
    boundary = injection_area
  [../]
  [./constant_injection_temperature]
    type = PresetBC
    variable = temperature
    value = 313
    boundary = injection_area
  [../]

  [./roller_tmax]
    type = PresetBC
    variable = disp_x
    value = 0
    boundary = tmax
  [../]
  [./roller_tmin]
    type = PresetBC
    variable = disp_y
    value = 0
    boundary = tmin
  [../]
  [./roller_top_bottom]
    type = PresetBC
    variable = disp_z
    value = 0
    boundary = 'top bottom'
  [../]
  [./cavity_pressure_x]
    type = Pressure
    boundary = injection_area
    variable = disp_x
    component = 0
    factor = 1E6
    use_displaced_mesh = false
  [../]
  [./cavity_pressure_y]
    type = Pressure
    boundary = injection_area
    variable = disp_y
    component = 1
    factor = 1E6
    use_displaced_mesh = false
  [../]
[]

[AuxVariables]
  [./stress_rr]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_pp]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[AuxKernels]
  [./stress_rr]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_rr
    scalar_type = RadialStress
    point1 = '0 0 0'
    point2 = '0 0 1'
  [../]
  [./stress_pp]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = stress_pp
    scalar_type = HoopStress
    point1 = '0 0 0'
    point2 = '0 0 1'
  [../]
[]

[Modules]
  [./FluidProperties]
    [./the_simple_fluid]
      type = SimpleFluidProperties
      bulk_modulus = 2E9
      viscosity = 1.0E-3
      density0 = 1000.0
      thermal_expansion = 0.0002
      cp = 4194
      cv = 4186
      porepressure_coefficient = 0
    [../]
  [../]
[]

[Materials]
  [./porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  [../]
  [./biot_modulus]
    type = PorousFlowConstantBiotModulus
    solid_bulk_compliance = 2E-7
    fluid_bulk_modulus = 1E7
  [../]
  [./permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  [../]
  [./permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  [../]

  [./thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  [../]
  [./rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  [../]
  [./thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caps aquifer'
  [../]

  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5E9
    poissons_ratio = 0.0
  [../]
  [./strain]
    type = ComputeSmallStrain
    eigenstrain_names = thermal_contribution
  [../]
  [./thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = temperature
    thermal_expansion_coeff = 0.001 # this is the linear thermal expansion coefficient
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 293
  [../]
  [./stress]
    type = ComputeLinearElasticStress
  [../]
[]

[Preconditioning]
  active = basic
  [./basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  [../]
  [./preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  [../]
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]

[Outputs]
  exodus = true
[]
(modules/porous_flow/examples/tutorial/04.i)

Input Parameters

  • porepressureThe name of the porepressure variable

    C++ Type:NonlinearVariableName

    Options:

    Description:The name of the porepressure variable

Required Parameters

  • fpuse_brine_materialThe name of the user object for fluid properties. Not required if use_brine is true.

    Default:use_brine_material

    C++ Type:UserObjectName

    Options:

    Description:The name of the user object for fluid properties. Not required if use_brine is true.

  • dictator_namedictatorThe name of the dictator user object that is created by this Action

    Default:dictator

    C++ Type:std::string

    Options:

    Description:The name of the dictator user object that is created by this Action

  • add_darcy_auxTrueAdd AuxVariables that record Darcy velocity

    Default:True

    C++ Type:bool

    Options:

    Description:Add AuxVariables that record Darcy velocity

  • temperature293.0For isothermal simulations, this is the temperature at which fluid properties (and stress-free strains) are evaluated at. Otherwise, this is the name of the temperature variable. Units = Kelvin

    Default:293.0

    C++ Type:std::vector

    Options:

    Description:For isothermal simulations, this is the temperature at which fluid properties (and stress-free strains) are evaluated at. Otherwise, this is the name of the temperature variable. Units = Kelvin

  • simulation_typetransientWhether a transient or steady-state simulation is being performed

    Default:transient

    C++ Type:MooseEnum

    Options:steady transient

    Description:Whether a transient or steady-state simulation is being performed

  • mass_fraction_varsList of variables that represent the mass fractions. With only one fluid component, this may be left empty. With N fluid components, the format is 'f_0 f_1 f_2 ... f_(N-1)'. That is, the N^th component need not be specified because f_N = 1 - (f_0 + f_1 + ... + f_(N-1)). It is best numerically to choose the N-1 mass fraction variables so that they represent the fluid components with small concentrations. This Action will associated the i^th mass fraction variable to the equation for the i^th fluid component, and the pressure variable to the N^th fluid component.

    C++ Type:std::vector

    Options:

    Description:List of variables that represent the mass fractions. With only one fluid component, this may be left empty. With N fluid components, the format is 'f_0 f_1 f_2 ... f_(N-1)'. That is, the N^th component need not be specified because f_N = 1 - (f_0 + f_1 + ... + f_(N-1)). It is best numerically to choose the N-1 mass fraction variables so that they represent the fluid components with small concentrations. This Action will associated the i^th mass fraction variable to the equation for the i^th fluid component, and the pressure variable to the N^th fluid component.

  • use_displaced_meshFalseUse displaced mesh computations in mechanical kernels

    Default:False

    C++ Type:bool

    Options:

    Description:Use displaced mesh computations in mechanical kernels

  • gravity0 0 -10Gravitational acceleration vector downwards (m/s^2)

    Default:0 0 -10

    C++ Type:libMesh::VectorValue

    Options:

    Description:Gravitational acceleration vector downwards (m/s^2)

  • multiply_by_densityFalseIf true, then the Kernels for fluid flow are multiplied by the fluid density. If false, this multiplication is not performed, which means the problem linearises, but that care must be taken when using other PorousFlow objects.

    Default:False

    C++ Type:bool

    Options:

    Description:If true, then the Kernels for fluid flow are multiplied by the fluid density. If false, this multiplication is not performed, which means the problem linearises, but that care must be taken when using other PorousFlow objects.

  • coupling_typeHydroThe type of simulation. For simulations involving Mechanical deformations, you will need to supply the correct Biot coefficient. For simulations involving Thermal flows, you will need an associated ConstantThermalExpansionCoefficient Material

    Default:Hydro

    C++ Type:MooseEnum

    Options:Hydro ThermoHydro HydroMechanical ThermoHydroMechanical

    Description:The type of simulation. For simulations involving Mechanical deformations, you will need to supply the correct Biot coefficient. For simulations involving Thermal flows, you will need an associated ConstantThermalExpansionCoefficient Material

  • inactiveIf specified blocks matching these identifiers will be skipped.

    C++ Type:std::vector

    Options:

    Description:If specified blocks matching these identifiers will be skipped.

  • use_brineFalseUse PorousFlowBrine material for the fluid phase

    Default:False

    C++ Type:bool

    Options:

    Description:Use PorousFlowBrine material for the fluid phase

  • thermal_eigenstrain_namethermal_eigenstrainThe eigenstrain_name used in the ComputeThermalExpansionEigenstrain. Only needed for thermally-coupled simulations with thermal expansion.

    Default:thermal_eigenstrain

    C++ Type:std::string

    Options:

    Description:The eigenstrain_name used in the ComputeThermalExpansionEigenstrain. Only needed for thermally-coupled simulations with thermal expansion.

  • active__all__ If specified only the blocks named will be visited and made active

    Default:__all__

    C++ Type:std::vector

    Options:

    Description:If specified only the blocks named will be visited and made active

  • nacl_index0Index of NaCl variable in mass_fraction_vars, for calculating brine properties. Only required if use_brine is true.

    Default:0

    C++ Type:unsigned int

    Options:

    Description:Index of NaCl variable in mass_fraction_vars, for calculating brine properties. Only required if use_brine is true.

  • number_aqueous_equilibrium0The number of secondary species in the aqueous-equilibrium reaction system. (Leave as zero if the simulation does not involve chemistry)

    Default:0

    C++ Type:unsigned int

    Options:

    Description:The number of secondary species in the aqueous-equilibrium reaction system. (Leave as zero if the simulation does not involve chemistry)

  • number_aqueous_kinetic0The number of secondary species in the aqueous-kinetic reaction system involved in precipitation and dissolution. (Leave as zero if the simulation does not involve chemistry)

    Default:0

    C++ Type:unsigned int

    Options:

    Description:The number of secondary species in the aqueous-kinetic reaction system involved in precipitation and dissolution. (Leave as zero if the simulation does not involve chemistry)

  • biot_coefficient1The Biot coefficient (relevant only for mechanically-coupled simulations)

    Default:1

    C++ Type:double

    Options:

    Description:The Biot coefficient (relevant only for mechanically-coupled simulations)

  • displacementsThe name of the displacement variables (relevant only for mechanically-coupled simulations)

    C++ Type:std::vector

    Options:

    Description:The name of the displacement variables (relevant only for mechanically-coupled simulations)

  • add_stress_auxTrueAdd AuxVariables that record effective stress

    Default:True

    C++ Type:bool

    Options:

    Description:Add AuxVariables that record effective stress

Optional Parameters