Porous Flow Tutorial Page 05. Using a realistic equation of state for the fluid

It is trivial to add a realistic equation of state to any PorousFlow simulation. For instance, the high-precision Water97 equation of state may be used:

# Darcy flow with heat advection and conduction, using Water97 properties
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
    input = annular
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [porepressure]
    initial_condition = 1E6
  []
  [temperature]
    initial_condition = 313
    scaling = 1E-8
  []
[]

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

[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 2E6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 333
    boundary = injection_area
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = Water97FluidProperties
  []
[]

(The name "the_simple_fluid" could also be changed to something more appropriate.)

In this case, the initial and boundary conditions must also be changed because the Water97 equation of state is not defined for (an absolute vacuum). For example:

[Variables]
  [porepressure]
    initial_condition = 1E6
  []
  [temperature]
    initial_condition = 313
    scaling = 1E-8
  []
[]

and

[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 2E6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 333
    boundary = injection_area
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = Water97FluidProperties
  []
[]

[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
[]

Using realistic high-precision equations of state can cause PorousFlow to run quite slowly, because the equations of state are so complicated to evaluate. It is always recommended to use TabulatedFluidProperties in the following way:

# Darcy flow with heat advection and conduction, using Water97 properties
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
    input = annular
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [porepressure]
    initial_condition = 1E6
  []
  [temperature]
    initial_condition = 313
    scaling = 1E-8
  []
[]

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

[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 2E6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 333
    boundary = injection_area
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [tabulated_water]
    type = TabulatedFluidProperties
    fp = true_water
    temperature_min = 275
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = water97_tabulated.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated.csv
  []
[]

and

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

Another advantage of using TabulatedfluidProperties is that the bounds on pressure and temperature imposed by the original "true" equation of state can be extrapolated past, which can help to remove problems of MOOSE trying to sample outside the original's region of validity (and thus causing the timestep to be cut).

Finally, in the case at hand, it is worth reminding the reader that PorousFlowBasicTHM is based on the assumptions of a constant, large fluid bulk modulus, and a constant fluid thermal expansion coefficient, which are incorrect for some fluids. User beware!

Start | Previous | Next