Start | Previous | Next

Porous Flow Tutorial Page 11. Two-phase THM borehole injection

Hold onto your hats, here we're going to build a two-phase THM simulation from scratch. Our starting point is the borehole-reservoir-caprock geometry developed in Page 00.

  • We would like to simulate cold CO injection into an initially water-filled reservoir.

  • There are two phases, and two components. The liquid phase is filled with water only, while the "gas" phases is filled with CO only ("gas" is in quotes because it may be supercritical).

  • There is capillarity between the water and CO.

  • The fluids have certain nontrivial relative permeability curves.

  • Heat is advected with the fluids as well as conducted through the fluid-rock system.

  • The rock can deform elastically.

  • There is no gravity (for simplicity as this effects initial stresses and fluid pressures, which is not difficult to include but the input file becomes even longer).

  • Full coupling between the fluids, the temperature and the displacements is used.

  • The porosity and permeability are functions of the porepressures, saturations, temperature and displacements

  • High precision equations of state are used for the water and CO.

  • At the injection area, CO is injected at a constant rate and at a constant temperature. It pushes mechanically against the borehole wall

  • At the outer boundary fluid is removed slowly using a PorousFlowSink.

  • In this model there are no vertical displacements

Variables

Let's step through the input file. This simulation uses the water and gas porepressures as the independent fluid pressures. A perfectly valid alternative would be to use water saturation and gas porepressure. A scaling factor is needed for the temperature and displacement variables as discussed in the convergence page. The Variables are given initial conditions appropriate to a deep CO sequestration site:

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [pwater]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 20E6
  []
  [pgas]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 20.1E6
  []
  [T]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 330
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-5
  []
  [disp_x]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-5
  []
  [disp_y]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-5
  []
[]
(modules/porous_flow/examples/tutorial/11.i)

Dictator

The PorousFlowDictator is given the variable names as well as the number of fluid phases and components.

  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater pgas T disp_x disp_y'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
(modules/porous_flow/examples/tutorial/11.i)

GlobalParams

The global parameters consist of various parameters that are required by several PorousFlow objects, and using GlobalParams allows the input file to remain organised and a little more succinct:

[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y disp_z'
  gravity = '0 0 0'
  biot_coefficient = 1.0
  PorousFlowDictator = dictator
[]
(modules/porous_flow/examples/tutorial/11.i)

Kernels

The Kernels block is rather big. Please refer to the governing equations page to see what equations are being modelled here. Virtually all the power of PorousFlow is being used. The only things that are missing are the following.

  • Chemical reactions and desorption are not used.

  • Fluid diffusion and dispersion is not active because water only exists in the liquid phase and CO in the gas phase.

  • There is no radioactive decay.

  • There is no plasticity, which also means there is no plastic heating.

  • The solid mechanics is assumed quasi-static.

The Kernels block is

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [mass_water_dot]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pwater
  []
  [flux_water]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pwater
  []
  [vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion<<<{"description": "Component_mass*rate_of_solid_volumetric_expansion.  This Kernel lumps the component mass to the nodes.", "href": "../../source/kernels/PorousFlowMassVolumetricExpansion.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pwater
  []
  [mass_co2_dot]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
  []
  [flux_co2]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 1
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
  []
  [vol_strain_rate_co2]
    type = PorousFlowMassVolumetricExpansion<<<{"description": "Component_mass*rate_of_solid_volumetric_expansion.  This Kernel lumps the component mass to the nodes.", "href": "../../source/kernels/PorousFlowMassVolumetricExpansion.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative<<<{"description": "Derivative of heat-energy-density wrt time", "href": "../../source/kernels/PorousFlowEnergyTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
  []
  [advection]
    type = PorousFlowHeatAdvection<<<{"description": "Fully-upwinded heat flux, advected by the fluid", "href": "../../source/kernels/PorousFlowHeatAdvection.html"}>>>
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
  []
  [conduction]
    type = PorousFlowHeatConduction<<<{"description": "Heat conduction in the Porous Flow module", "href": "../../source/kernels/PorousFlowHeatConduction.html"}>>>
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
  []
  [vol_strain_rate_heat]
    type = PorousFlowHeatVolumetricExpansion<<<{"description": "Energy-density*rate_of_solid_volumetric_expansion.  The energy-density is lumped to the nodes", "href": "../../source/kernels/PorousFlowHeatVolumetricExpansion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
  []
  [grad_stress_x]
    type = StressDivergenceTensors<<<{"description": "Stress divergence kernel for the Cartesian coordinate system", "href": "../../source/kernels/StressDivergenceTensors.html"}>>>
    temperature<<<{"description": "The name of the temperature variable used in the ComputeThermalExpansionEigenstrain.  (Not required for simulations without temperature coupling.)"}>>> = T
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    eigenstrain_names<<<{"description": "List of eigenstrains used in the strain calculation. Used for computing their derivatives for off-diagonal Jacobian terms."}>>> = thermal_contribution
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    component<<<{"description": "An integer corresponding to the direction the variable this kernel acts in. (0 for x, 1 for y, 2 for z)"}>>> = 0
  []
  [poro_x]
    type = PorousFlowEffectiveStressCoupling<<<{"description": "Implements the weak form of the expression biot_coefficient * grad(effective fluid pressure)", "href": "../../source/kernels/PorousFlowEffectiveStressCoupling.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    component<<<{"description": "The component (0 for x, 1 for y and 2 for z) of grad(P)"}>>> = 0
  []
  [grad_stress_y]
    type = StressDivergenceTensors<<<{"description": "Stress divergence kernel for the Cartesian coordinate system", "href": "../../source/kernels/StressDivergenceTensors.html"}>>>
    temperature<<<{"description": "The name of the temperature variable used in the ComputeThermalExpansionEigenstrain.  (Not required for simulations without temperature coupling.)"}>>> = T
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    eigenstrain_names<<<{"description": "List of eigenstrains used in the strain calculation. Used for computing their derivatives for off-diagonal Jacobian terms."}>>> = thermal_contribution
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    component<<<{"description": "An integer corresponding to the direction the variable this kernel acts in. (0 for x, 1 for y, 2 for z)"}>>> = 1
  []
  [poro_y]
    type = PorousFlowEffectiveStressCoupling<<<{"description": "Implements the weak form of the expression biot_coefficient * grad(effective fluid pressure)", "href": "../../source/kernels/PorousFlowEffectiveStressCoupling.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    component<<<{"description": "The component (0 for x, 1 for y and 2 for z) of grad(P)"}>>> = 1
  []
[]
(modules/porous_flow/examples/tutorial/11.i)

AuxVariables and AuxKernels

Some AuxVariables are defined that need further explanation.

  • The solid mechanics needs 3 displacement variables. In the assumptions above, the vertical displacement is always zero. The best way of defining it is as an AuxVariable without an AuxKernel so that it will always be zero, but then it may be coupled into various MOOSE objects (using the displacements = in the GlobalParams)

  • CO is pumped into this model with a constant rate. To achieve this, the CO pressure in the borehole must gradually increase over time. This means the fluid in the borehole will push against the borehole wall with increasing pressure. Hence this pressure must be recorded and coupled into the mechanical BCs. It is recorded in the effective_fluid_pressure AuxVariable.

  • The mass fractions of each component in each phase must be defined, even if they are fixed for always. Since they are unchanging they are most conveniently represented by AuxVariables with certain initial conditions and no AuxKernels.

  • The other AuxVariables are useful for visualising the results.

[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [disp_z]
  []
  [effective_fluid_pressure]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [mass_frac_phase0_species0]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 1 # all water in phase=0
  []
  [mass_frac_phase1_species0]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0 # no water in phase=1
  []
  [sgas]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [swater]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [stress_rr]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [stress_tt]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [stress_zz]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [porosity]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
  [effective_fluid_pressure]
    type = ParsedAux<<<{"description": "Sets a field variable value to the evaluation of a parsed expression.", "href": "../../source/auxkernels/ParsedAux.html"}>>>
    coupled_variables<<<{"description": "Vector of coupled variable names"}>>> = 'pwater pgas swater sgas'
    expression<<<{"description": "Parsed function expression to compute"}>>> = 'pwater * swater + pgas * sgas'
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = effective_fluid_pressure
  []
  [swater]
    type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = swater
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = saturation
    phase<<<{"description": "The index of the phase this auxillary kernel acts on"}>>> = 0
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [sgas]
    type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = sgas
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = saturation
    phase<<<{"description": "The index of the phase this auxillary kernel acts on"}>>> = 1
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [stress_rr]
    type = RankTwoScalarAux<<<{"description": "Compute a scalar property of a RankTwoTensor", "href": "../../source/auxscalarkernels/RankTwoScalarAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_rr
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    scalar_type<<<{"description": "Type of scalar output"}>>> = RadialStress
    point1<<<{"description": "Start point for axis used to calculate some cylindrical material tensor quantities"}>>> = '0 0 0'
    point2<<<{"description": "End point for axis used to calculate some material tensor quantities"}>>> = '0 0 1'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [stress_tt]
    type = RankTwoScalarAux<<<{"description": "Compute a scalar property of a RankTwoTensor", "href": "../../source/auxscalarkernels/RankTwoScalarAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_tt
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    scalar_type<<<{"description": "Type of scalar output"}>>> = HoopStress
    point1<<<{"description": "Start point for axis used to calculate some cylindrical material tensor quantities"}>>> = '0 0 0'
    point2<<<{"description": "End point for axis used to calculate some material tensor quantities"}>>> = '0 0 1'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [stress_zz]
    type = RankTwoAux<<<{"description": "Access a component of a RankTwoTensor", "href": "../../source/auxkernels/RankTwoAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_zz
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    index_i<<<{"description": "The index i of ij for the tensor to output (0, 1, 2)"}>>> = 2
    index_j<<<{"description": "The index j of ij for the tensor to output (0, 1, 2)"}>>> = 2
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [porosity]
    type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = porosity
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = porosity
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]
(modules/porous_flow/examples/tutorial/11.i)

Boundary conditions

The boundary conditions for the displacements are roller on the sides, fixed at the top and bottom (an arbitrary choice made by the creator of this input file) and Pressure boundary conditions on the injection_area:

  [roller_tmax]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = dmax
  []
  [roller_tmin]
    type = DirichletBC
    variable = disp_y
    value = 0
    boundary = dmin
  []
  [pinned_top_bottom_x]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = 'top bottom'
  []
  [pinned_top_bottom_y]
    type = DirichletBC
    variable = disp_y
    value = 0
    boundary = 'top bottom'
  []
  [cavity_pressure_x]
    type = Pressure
    boundary = injection_area
    variable = disp_x
    component = 0
    postprocessor = constrained_effective_fluid_pressure_at_wellbore
    use_displaced_mesh = false
  []
  [cavity_pressure_y]
    type = Pressure
    boundary = injection_area
    variable = disp_y
    component = 1
    postprocessor = constrained_effective_fluid_pressure_at_wellbore
    use_displaced_mesh = false
  []
(modules/porous_flow/examples/tutorial/11.i)

Notice the constrained_effective_fluid_pressure_at_wellbore. This is almost the effective_fluid_pressure AuxVariable defined above, evaluated at the injection_area. But there is a problem at the first timestep, because this uses Material properties that are not properly initialised. So a little bit of deception is used. Firstly, the AuxVariable is evaluated at the injection_area and put into a Postprocessor:

[Postprocessors]
  [effective_fluid_pressure_at_wellbore]
    type = PointValue
    variable = effective_fluid_pressure
    point = '1 0 0'
    execute_on = timestep_begin
    use_displaced_mesh = false
  []
(modules/porous_flow/examples/tutorial/11.i)

Then a Function is made that returns either the value of this Postprocessor or 20E6 (the initial reservoir pressure)

[Functions<<<{"href": "../../syntax/Functions/index.html"}>>>]
  [constrain_effective_fluid_pressure]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = effective_fluid_pressure_at_wellbore
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = effective_fluid_pressure_at_wellbore
    expression<<<{"description": "The user defined function."}>>> = 'max(effective_fluid_pressure_at_wellbore, 20E6)'
  []
[]
(modules/porous_flow/examples/tutorial/11.i)

Finally, the value of this Function is placed into the Postprocessor used in the Pressure BC:

  [constrained_effective_fluid_pressure_at_wellbore]
    type = FunctionValuePostprocessor
    function = constrain_effective_fluid_pressure
    execute_on = timestep_begin
(modules/porous_flow/examples/tutorial/11.i)

The boundary conditions for temperature is a simple preset DirichletBC on the injection_area:

  [cold_co2]
    type = DirichletBC
    boundary = injection_area
    variable = T
    value = 290 # injection temperature
    use_displaced_mesh = false
  []
(modules/porous_flow/examples/tutorial/11.i)

The boundary conditions for the fluids at the injection_area is just a constant injection of CO with rate kg.s.m:

  [constant_co2_injection]
    type = PorousFlowSink
    boundary = injection_area
    variable = pgas
    fluid_phase = 1
    flux_function = -1E-4
    use_displaced_mesh = false
  []
(modules/porous_flow/examples/tutorial/11.i)

At the outer boundary, water and CO are removed if their porepressures rise above their initial values. A PorousFlowPiecewiseLinearSink is used, with an imaginary boundary at fixed porepressure sitting at a distance of m outside the model. The procedure of constructing this sink is described fully in the boundaries documentation. The input-file blocks are

  [outer_water_removal]
    type = PorousFlowPiecewiseLinearSink
    boundary = rmax
    variable = pwater
    fluid_phase = 0
    pt_vals = '0 1E9'
    multipliers = '0 1E8'
    PT_shift = 20E6
    use_mobility = true
    use_relperm = true
    use_displaced_mesh = false
  []
  [outer_co2_removal]
    type = PorousFlowPiecewiseLinearSink
    boundary = rmax
    variable = pgas
    fluid_phase = 1
    pt_vals = '0 1E9'
    multipliers = '0 1E8'
    PT_shift = 20.1E6
    use_mobility = true
    use_relperm = true
    use_displaced_mesh = false
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [tabulated_water]
    type = TabulatedFluidProperties
    fp = true_water
    temperature_min = 275
    pressure_max = 1E8
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = water97_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated_11.csv
  []
  [true_co2]
    type = CO2FluidProperties
  []
  [tabulated_co2]
    type = TabulatedFluidProperties
    fp = true_co2
    temperature_min = 275
    pressure_max = 1E8
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = co2_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = co2_tabulated_11.csv
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = T
  []
  [saturation_calculator]
    type = PorousFlow2PhasePP
    phase0_porepressure = pwater
    phase1_porepressure = pgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'mass_frac_phase0_species0 mass_frac_phase1_species0'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_water
    phase = 0
  []
  [co2]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_co2
    phase = 1
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 4
    s_res = 0.1
    sum_s_res = 0.2
    phase = 0
  []
  [relperm_co2]
    type = PorousFlowRelativePermeabilityBC
    nw_phase = true
    lambda = 2
    s_res = 0.1
    sum_s_res = 0.2
    phase = 1
  []
  [porosity_mat]
    type = PorousFlowPorosity
    fluid = true
    mechanical = true
    thermal = true
    porosity_zero = 0.1
    reference_temperature = 330
    reference_porepressure = 20E6
    thermal_expansion_coeff = 15E-6 # volumetric
    solid_bulk = 8E9 # unimportant since biot = 1
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityKozenyCarman
    block = aquifer
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-12
  []
  [permeability_caps]
    type = PorousFlowPermeabilityKozenyCarman
    block = caps
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-15
    k_anisotropy = '1 0 0  0 1 0  0 0 0.1'
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '2 0 0  0 2 0  0 0 2'
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1100
    density = 2300
  []

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5E9
    poissons_ratio = 0.0
  []
  [strain]
    type = ComputeSmallStrain
    eigenstrain_names = 'thermal_contribution initial_stress'
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = T
    thermal_expansion_coeff = 5E-6 # this is the linear thermal expansion coefficient
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 330
  []
  [initial_strain]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = '20E6 0 0  0 20E6 0  0 0 20E6'
    eigenstrain_name = initial_stress
  []
  [stress]
    type = ComputeLinearElasticStress
  []

  [effective_fluid_pressure_mat]
    type = PorousFlowEffectiveFluidPressure
  []
  [volumetric_strain]
    type = PorousFlowVolumetricStrain
  []
[]

[Postprocessors]
  [effective_fluid_pressure_at_wellbore]
    type = PointValue
    variable = effective_fluid_pressure
    point = '1 0 0'
    execute_on = timestep_begin
    use_displaced_mesh = false
  []
  [constrained_effective_fluid_pressure_at_wellbore]
    type = FunctionValuePostprocessor
    function = constrain_effective_fluid_pressure
    execute_on = timestep_begin
  []
[]

[Functions]
  [constrain_effective_fluid_pressure]
    type = ParsedFunction
    symbol_names = effective_fluid_pressure_at_wellbore
    symbol_values = effective_fluid_pressure_at_wellbore
    expression = 'max(effective_fluid_pressure_at_wellbore, 20E6)'
  []
[]

[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 = 1E3
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1E3
    growth_factor = 1.2
    optimal_iterations = 10
  []
  nl_abs_tol = 1E-7
[]

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

Fluid properties

High-precision equations of state are used for both the water and the CO. Before running the simulation, these are tabulated, and the tabulated versions are used by MOOSE in all computations, which shortens the simulation time:

# Two-phase borehole injection problem
[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  [annular]
    type = AnnularMeshGenerator<<<{"description": "For rmin>0: creates an annular mesh of QUAD4 elements. For rmin=0: creates a disc mesh of QUAD4 and TRI3 elements. Boundary sidesets are created at rmax and rmin, and given these names. If dmin!=0 and dmax!=360, a sector of an annulus or disc is created. In this case boundary sidesets are also created at dmin and dmax, and given these names", "href": "../../source/meshgenerators/AnnularMeshGenerator.html"}>>>
    nr<<<{"description": "Number of elements in the radial direction"}>>> = 10
    rmin<<<{"description": "Inner radius.  If rmin=0 then a disc mesh (with no central hole) will be created."}>>> = 1.0
    rmax<<<{"description": "Outer radius"}>>> = 10
    growth_r<<<{"description": "The ratio of radial sizes of successive rings of elements"}>>> = 1.4
    nt<<<{"description": "Number of elements in the angular direction"}>>> = 4
    dmin<<<{"description": "Minimum degree, measured in degrees anticlockwise from x axis"}>>> = 0
    dmax<<<{"description": "Maximum angle, measured in degrees anticlockwise from x axis. If dmin=0 and dmax=360 an annular mesh is created. Otherwise, only a sector of an annulus is created"}>>> = 90
  []
  [make3D]
    input<<<{"description": "the mesh we want to extrude"}>>> = annular
    type = MeshExtruderGenerator<<<{"description": "Takes a 1D or 2D mesh and extrudes the entire structure along the specified axis increasing the dimensionality of the mesh.", "href": "../../source/meshgenerators/MeshExtruderGenerator.html"}>>>
    extrusion_vector<<<{"description": "The direction and length of the extrusion"}>>> = '0 0 12'
    num_layers<<<{"description": "The number of layers in the extruded mesh"}>>> = 3
    bottom_sideset<<<{"description": "The boundary that will be applied to the bottom of the extruded mesh"}>>> = 'bottom'
    top_sideset<<<{"description": "The boundary that will be to the top of the extruded mesh"}>>> = 'top'
  []
  [shift_down]
    type = TransformGenerator<<<{"description": "Applies a linear transform to the entire mesh.", "href": "../../source/meshgenerators/TransformGenerator.html"}>>>
    transform<<<{"description": "The type of transformation to perform (TRANSLATE, TRANSLATE_CENTER_ORIGIN, TRANSLATE_MIN_ORIGIN, ROTATE, SCALE)"}>>> = TRANSLATE
    vector_value<<<{"description": "The value to use for the transformation. When using TRANSLATE or SCALE, the xyz coordinates are applied in each direction respectively. When using ROTATE, the values are interpreted as the Euler angles phi, theta and psi given in degrees."}>>> = '0 0 -6'
    input<<<{"description": "The mesh we want to modify"}>>> = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator<<<{"description": "Changes the subdomain ID of elements either (XOR) inside or outside the specified box to the specified ID.", "href": "../../source/meshgenerators/SubdomainBoundingBoxGenerator.html"}>>>
    block_id<<<{"description": "Subdomain id to set for inside/outside the bounding box"}>>> = 1
    bottom_left<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '0 0 -2'
    top_right<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '10 10 2'
    input<<<{"description": "The mesh we want to modify"}>>> = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset<<<{"description": "A MeshGenerator that adds element sides to a sideset if the centroid of the side satisfies the `combinatorial_geometry` expression.", "href": "../../source/meshgenerators/ParsedGenerateSideset.html"}>>>
    combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'x*x+y*y<1.01'
    included_subdomains<<<{"description": "A set of subdomain names or ids whose sides will be included in the new sidesets. A side is only added if the subdomain id of the corresponding element is in this set."}>>> = 1
    new_sideset_name<<<{"description": "The name of the new sideset"}>>> = 'injection_area'
    input<<<{"description": "The mesh we want to modify"}>>> = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator<<<{"description": "Changes the block IDs and/or block names for a given set of blocks defined by either block ID or block name. The changes are independent of ordering. The merging of blocks is supported.", "href": "../../source/meshgenerators/RenameBlockGenerator.html"}>>>
    old_block<<<{"description": "Elements with these block ID(s)/name(s) will be given the new block information specified in 'new_block'"}>>> = '0 1'
    new_block<<<{"description": "The new block ID(s)/name(s) to be given by the elements defined in 'old_block'."}>>> = 'caps aquifer'
    input<<<{"description": "The mesh we want to modify"}>>> = 'injection_area'
  []
[]

[UserObjects<<<{"href": "../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pwater pgas T disp_x disp_y'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 2
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1E-6
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.6
  []
[]

[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
  displacements = 'disp_x disp_y disp_z'
  gravity = '0 0 0'
  biot_coefficient = 1.0
  PorousFlowDictator = dictator
[]

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [pwater]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 20E6
  []
  [pgas]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 20.1E6
  []
  [T]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 330
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-5
  []
  [disp_x]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-5
  []
  [disp_y]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-5
  []
[]

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [mass_water_dot]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pwater
  []
  [flux_water]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pwater
  []
  [vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion<<<{"description": "Component_mass*rate_of_solid_volumetric_expansion.  This Kernel lumps the component mass to the nodes.", "href": "../../source/kernels/PorousFlowMassVolumetricExpansion.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pwater
  []
  [mass_co2_dot]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
  []
  [flux_co2]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 1
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
  []
  [vol_strain_rate_co2]
    type = PorousFlowMassVolumetricExpansion<<<{"description": "Component_mass*rate_of_solid_volumetric_expansion.  This Kernel lumps the component mass to the nodes.", "href": "../../source/kernels/PorousFlowMassVolumetricExpansion.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative<<<{"description": "Derivative of heat-energy-density wrt time", "href": "../../source/kernels/PorousFlowEnergyTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
  []
  [advection]
    type = PorousFlowHeatAdvection<<<{"description": "Fully-upwinded heat flux, advected by the fluid", "href": "../../source/kernels/PorousFlowHeatAdvection.html"}>>>
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
  []
  [conduction]
    type = PorousFlowHeatConduction<<<{"description": "Heat conduction in the Porous Flow module", "href": "../../source/kernels/PorousFlowHeatConduction.html"}>>>
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
  []
  [vol_strain_rate_heat]
    type = PorousFlowHeatVolumetricExpansion<<<{"description": "Energy-density*rate_of_solid_volumetric_expansion.  The energy-density is lumped to the nodes", "href": "../../source/kernels/PorousFlowHeatVolumetricExpansion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
  []
  [grad_stress_x]
    type = StressDivergenceTensors<<<{"description": "Stress divergence kernel for the Cartesian coordinate system", "href": "../../source/kernels/StressDivergenceTensors.html"}>>>
    temperature<<<{"description": "The name of the temperature variable used in the ComputeThermalExpansionEigenstrain.  (Not required for simulations without temperature coupling.)"}>>> = T
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    eigenstrain_names<<<{"description": "List of eigenstrains used in the strain calculation. Used for computing their derivatives for off-diagonal Jacobian terms."}>>> = thermal_contribution
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    component<<<{"description": "An integer corresponding to the direction the variable this kernel acts in. (0 for x, 1 for y, 2 for z)"}>>> = 0
  []
  [poro_x]
    type = PorousFlowEffectiveStressCoupling<<<{"description": "Implements the weak form of the expression biot_coefficient * grad(effective fluid pressure)", "href": "../../source/kernels/PorousFlowEffectiveStressCoupling.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    component<<<{"description": "The component (0 for x, 1 for y and 2 for z) of grad(P)"}>>> = 0
  []
  [grad_stress_y]
    type = StressDivergenceTensors<<<{"description": "Stress divergence kernel for the Cartesian coordinate system", "href": "../../source/kernels/StressDivergenceTensors.html"}>>>
    temperature<<<{"description": "The name of the temperature variable used in the ComputeThermalExpansionEigenstrain.  (Not required for simulations without temperature coupling.)"}>>> = T
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    eigenstrain_names<<<{"description": "List of eigenstrains used in the strain calculation. Used for computing their derivatives for off-diagonal Jacobian terms."}>>> = thermal_contribution
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    component<<<{"description": "An integer corresponding to the direction the variable this kernel acts in. (0 for x, 1 for y, 2 for z)"}>>> = 1
  []
  [poro_y]
    type = PorousFlowEffectiveStressCoupling<<<{"description": "Implements the weak form of the expression biot_coefficient * grad(effective fluid pressure)", "href": "../../source/kernels/PorousFlowEffectiveStressCoupling.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    component<<<{"description": "The component (0 for x, 1 for y and 2 for z) of grad(P)"}>>> = 1
  []
[]

[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [disp_z]
  []
  [effective_fluid_pressure]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [mass_frac_phase0_species0]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 1 # all water in phase=0
  []
  [mass_frac_phase1_species0]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0 # no water in phase=1
  []
  [sgas]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [swater]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [stress_rr]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [stress_tt]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [stress_zz]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [porosity]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
  [effective_fluid_pressure]
    type = ParsedAux<<<{"description": "Sets a field variable value to the evaluation of a parsed expression.", "href": "../../source/auxkernels/ParsedAux.html"}>>>
    coupled_variables<<<{"description": "Vector of coupled variable names"}>>> = 'pwater pgas swater sgas'
    expression<<<{"description": "Parsed function expression to compute"}>>> = 'pwater * swater + pgas * sgas'
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = effective_fluid_pressure
  []
  [swater]
    type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = swater
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = saturation
    phase<<<{"description": "The index of the phase this auxillary kernel acts on"}>>> = 0
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [sgas]
    type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = sgas
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = saturation
    phase<<<{"description": "The index of the phase this auxillary kernel acts on"}>>> = 1
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [stress_rr]
    type = RankTwoScalarAux<<<{"description": "Compute a scalar property of a RankTwoTensor", "href": "../../source/auxscalarkernels/RankTwoScalarAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_rr
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    scalar_type<<<{"description": "Type of scalar output"}>>> = RadialStress
    point1<<<{"description": "Start point for axis used to calculate some cylindrical material tensor quantities"}>>> = '0 0 0'
    point2<<<{"description": "End point for axis used to calculate some material tensor quantities"}>>> = '0 0 1'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [stress_tt]
    type = RankTwoScalarAux<<<{"description": "Compute a scalar property of a RankTwoTensor", "href": "../../source/auxscalarkernels/RankTwoScalarAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_tt
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    scalar_type<<<{"description": "Type of scalar output"}>>> = HoopStress
    point1<<<{"description": "Start point for axis used to calculate some cylindrical material tensor quantities"}>>> = '0 0 0'
    point2<<<{"description": "End point for axis used to calculate some material tensor quantities"}>>> = '0 0 1'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [stress_zz]
    type = RankTwoAux<<<{"description": "Access a component of a RankTwoTensor", "href": "../../source/auxkernels/RankTwoAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = stress_zz
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = stress
    index_i<<<{"description": "The index i of ij for the tensor to output (0, 1, 2)"}>>> = 2
    index_j<<<{"description": "The index j of ij for the tensor to output (0, 1, 2)"}>>> = 2
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
  [porosity]
    type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = porosity
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = porosity
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
  [roller_tmax]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    value<<<{"description": "Value of the BC"}>>> = 0
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = dmax
  []
  [roller_tmin]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    value<<<{"description": "Value of the BC"}>>> = 0
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = dmin
  []
  [pinned_top_bottom_x]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    value<<<{"description": "Value of the BC"}>>> = 0
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'top bottom'
  []
  [pinned_top_bottom_y]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    value<<<{"description": "Value of the BC"}>>> = 0
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'top bottom'
  []
  [cavity_pressure_x]
    type = Pressure<<<{"description": "Applies a pressure on a given boundary in a given direction", "href": "../../source/bcs/Pressure.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    component<<<{"description": "The component for the pressure"}>>> = 0
    postprocessor<<<{"description": "Postprocessor that will supply the pressure value"}>>> = constrained_effective_fluid_pressure_at_wellbore
    use_displaced_mesh<<<{"description": "Whether to use the displaced mesh."}>>> = false
  []
  [cavity_pressure_y]
    type = Pressure<<<{"description": "Applies a pressure on a given boundary in a given direction", "href": "../../source/bcs/Pressure.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    component<<<{"description": "The component for the pressure"}>>> = 1
    postprocessor<<<{"description": "Postprocessor that will supply the pressure value"}>>> = constrained_effective_fluid_pressure_at_wellbore
    use_displaced_mesh<<<{"description": "Whether to use the displaced mesh."}>>> = false
  []

  [cold_co2]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../source/bcs/DirichletBC.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    value<<<{"description": "Value of the BC"}>>> = 290 # injection temperature
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation.  Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
  []
  [constant_co2_injection]
    type = PorousFlowSink<<<{"description": "Applies a flux sink to a boundary.", "href": "../../source/bcs/PorousFlowSink.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
    fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy.  If not supplied, then this BC can only be a function of temperature"}>>> = 1
    flux_function<<<{"description": "The flux.  The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE.  The functional form is useful for spatially or temporally varying sinks.  Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = -1E-4
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation.  Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
  []

  [outer_water_removal]
    type = PorousFlowPiecewiseLinearSink<<<{"description": "Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a piecewise linear function.", "href": "../../source/bcs/PorousFlowPiecewiseLinearSink.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = rmax
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pwater
    fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy.  If not supplied, then this BC can only be a function of temperature"}>>> = 0
    pt_vals<<<{"description": "Tuple of pressure values (for the fluid_phase specified).  Must be monotonically increasing.  For heat fluxes that don't involve fluids, these are temperature values"}>>> = '0 1E9'
    multipliers<<<{"description": "Tuple of multiplying values.  The flux values are multiplied by these."}>>> = '0 1E8'
    PT_shift<<<{"description": "Whenever the sink is an explicit function of porepressure (such as a PiecewiseLinear function) the argument of the function is set to P - PT_shift instead of simply P.  Similarly for temperature.  PT_shift does not enter into any use_* calculations."}>>> = 20E6
    use_mobility<<<{"description": "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary.  In this case bare_flux is measured in Pa.m^-1.  This can be used in conjunction with other use_*"}>>> = true
    use_relperm<<<{"description": "If true, then fluxes are multiplied by relative permeability.  This can be used in conjunction with other use_*"}>>> = true
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation.  Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
  []
  [outer_co2_removal]
    type = PorousFlowPiecewiseLinearSink<<<{"description": "Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a piecewise linear function.", "href": "../../source/bcs/PorousFlowPiecewiseLinearSink.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = rmax
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
    fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy.  If not supplied, then this BC can only be a function of temperature"}>>> = 1
    pt_vals<<<{"description": "Tuple of pressure values (for the fluid_phase specified).  Must be monotonically increasing.  For heat fluxes that don't involve fluids, these are temperature values"}>>> = '0 1E9'
    multipliers<<<{"description": "Tuple of multiplying values.  The flux values are multiplied by these."}>>> = '0 1E8'
    PT_shift<<<{"description": "Whenever the sink is an explicit function of porepressure (such as a PiecewiseLinear function) the argument of the function is set to P - PT_shift instead of simply P.  Similarly for temperature.  PT_shift does not enter into any use_* calculations."}>>> = 20.1E6
    use_mobility<<<{"description": "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary.  In this case bare_flux is measured in Pa.m^-1.  This can be used in conjunction with other use_*"}>>> = true
    use_relperm<<<{"description": "If true, then fluxes are multiplied by relative permeability.  This can be used in conjunction with other use_*"}>>> = true
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation.  Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
  []
[]

[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
  [true_water]
    type = Water97FluidProperties<<<{"description": "Fluid properties for water and steam (H2O) using IAPWS-IF97", "href": "../../source/fluidproperties/Water97FluidProperties.html"}>>>
  []
  [tabulated_water]
    type = TabulatedFluidProperties<<<{"description": "Fluid properties using bicubic interpolation on tabulated values provided", "href": "../../source/fluidproperties/TabulatedBicubicFluidProperties.html"}>>>
    fp<<<{"description": "The name of the FluidProperties UserObject"}>>> = true_water
    temperature_min<<<{"description": "Minimum temperature for tabulated data."}>>> = 275
    pressure_max<<<{"description": "Maximum pressure for tabulated data."}>>> = 1E8
    interpolated_properties<<<{"description": "Properties to interpolate if no data file is provided"}>>> = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file<<<{"description": "Name of the CSV file which can be output with the tabulation. This file can then be read as a 'fluid_property_file'"}>>> = water97_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated_11.csv
  []
  [true_co2]
    type = CO2FluidProperties<<<{"description": "Fluid properties for carbon dioxide (CO2) using the Span & Wagner EOS", "href": "../../source/fluidproperties/CO2FluidProperties.html"}>>>
  []
  [tabulated_co2]
    type = TabulatedFluidProperties<<<{"description": "Fluid properties using bicubic interpolation on tabulated values provided", "href": "../../source/fluidproperties/TabulatedBicubicFluidProperties.html"}>>>
    fp<<<{"description": "The name of the FluidProperties UserObject"}>>> = true_co2
    temperature_min<<<{"description": "Minimum temperature for tabulated data."}>>> = 275
    pressure_max<<<{"description": "Maximum pressure for tabulated data."}>>> = 1E8
    interpolated_properties<<<{"description": "Properties to interpolate if no data file is provided"}>>> = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file<<<{"description": "Name of the CSV file which can be output with the tabulation. This file can then be read as a 'fluid_property_file'"}>>> = co2_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = co2_tabulated_11.csv
  []
[]
(modules/porous_flow/examples/tutorial/11.i)

Materials

The capillary pressure relationship is defined by the UserObject:

  [pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
(modules/porous_flow/examples/tutorial/11.i)

As explained on Page 09 and Page 10, there are a set of "fundamental Materials" that compute all porepressures, saturations, temperature and mass fractions as Material properties (as well as their gradients, and the derivatives with respect to the primary Variables, etc). In the case at hand, these fundamental Materials are:

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = T
  []
  [saturation_calculator]
    type = PorousFlow2PhasePP
    phase0_porepressure = pwater
    phase1_porepressure = pgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'mass_frac_phase0_species0 mass_frac_phase1_species0'
  []
(modules/porous_flow/examples/tutorial/11.i)

The water and CO densities, viscosities, enthalpies, and internal energies (as well as derivatives of these) are computed by

  [water]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_water
    phase = 0
  []
  [co2]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_co2
    phase = 1
  []
(modules/porous_flow/examples/tutorial/11.i)

A Corey type of relative permeability is chosen for the liquid phase, and a Brooks-Corey type of relative permeability is chosen for the gas phase:

  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 4
    s_res = 0.1
    sum_s_res = 0.2
    phase = 0
  []
  [relperm_co2]
    type = PorousFlowRelativePermeabilityBC
    nw_phase = true
    lambda = 2
    s_res = 0.1
    sum_s_res = 0.2
    phase = 1
  []
  [porosity_mat]
    type = PorousFlowPorosity
    fluid = true
    mechanical = true
    thermal = true
    porosity_zero = 0.1
    reference_temperature = 330
    reference_porepressure = 20E6
    thermal_expansion_coeff = 15E-6 # volumetric
    solid_bulk = 8E9 # unimportant since biot = 1
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityKozenyCarman
    block = aquifer
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-12
  []
  [permeability_caps]
    type = PorousFlowPermeabilityKozenyCarman
    block = caps
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-15
    k_anisotropy = '1 0 0  0 1 0  0 0 0.1'
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '2 0 0  0 2 0  0 0 2'
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1100
    density = 2300
  []

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5E9
    poissons_ratio = 0.0
  []
  [strain]
    type = ComputeSmallStrain
    eigenstrain_names = 'thermal_contribution initial_stress'
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = T
    thermal_expansion_coeff = 5E-6 # this is the linear thermal expansion coefficient
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 330
  []
  [initial_strain]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = '20E6 0 0  0 20E6 0  0 0 20E6'
    eigenstrain_name = initial_stress
  []
  [stress]
    type = ComputeLinearElasticStress
  []

  [effective_fluid_pressure_mat]
    type = PorousFlowEffectiveFluidPressure
  []
  [volumetric_strain]
    type = PorousFlowVolumetricStrain
  []
[]

[Postprocessors]
  [effective_fluid_pressure_at_wellbore]
    type = PointValue
    variable = effective_fluid_pressure
    point = '1 0 0'
    execute_on = timestep_begin
    use_displaced_mesh = false
  []
  [constrained_effective_fluid_pressure_at_wellbore]
    type = FunctionValuePostprocessor
    function = constrain_effective_fluid_pressure
    execute_on = timestep_begin
  []
[]

[Functions]
  [constrain_effective_fluid_pressure]
    type = ParsedFunction
    symbol_names = effective_fluid_pressure_at_wellbore
    symbol_values = effective_fluid_pressure_at_wellbore
    expression = 'max(effective_fluid_pressure_at_wellbore, 20E6)'
  []
[]

[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 = 1E3
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1E3
    growth_factor = 1.2
    optimal_iterations = 10
  []
  nl_abs_tol = 1E-7
[]

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

Porosity is chosen to depend on porepressures, saturations (actually just the effective porepressure), temperature and mechanical strain using:

  [porosity_mat]
    type = PorousFlowPorosity
    fluid = true
    mechanical = true
    thermal = true
    porosity_zero = 0.1
    reference_temperature = 330
    reference_porepressure = 20E6
    thermal_expansion_coeff = 15E-6 # volumetric
    solid_bulk = 8E9 # unimportant since biot = 1
  []
(modules/porous_flow/examples/tutorial/11.i)

Permeability is chosen to follow the Kozeny-Carman relationship:

  [permeability_aquifer]
    type = PorousFlowPermeabilityKozenyCarman
    block = aquifer
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-12
  []
  [permeability_caps]
    type = PorousFlowPermeabilityKozenyCarman
    block = caps
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-15
    k_anisotropy = '1 0 0  0 1 0  0 0 0.1'
  []
(modules/porous_flow/examples/tutorial/11.i)

The rock thermal conductivity is chosen to be independent of water saturation and to be isotropic (and independent of rock type — reservoir or cap-rock):

  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '2 0 0  0 2 0  0 0 2'
  []
(modules/porous_flow/examples/tutorial/11.i)

while the rock internal energy is also constant:

  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1100
    density = 2300
  []
(modules/porous_flow/examples/tutorial/11.i)

The elasticity tensor of the rock (both reservoir and cap-rock) is assumed isotropic with a Young's modulus of 5GPa:

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5E9
    poissons_ratio = 0.0
  []
(modules/porous_flow/examples/tutorial/11.i)

The strain calculator must take into consideration both the thermal strain (see governing equations) as well as initial effective stress. The initial total stress is assumed to be zero (for simplicity, not because it is physically very likely, but a nonzero value doesn't change the results much), so the initial effective stress is just the initial porepressure

  [strain]
    type = ComputeSmallStrain
    eigenstrain_names = 'thermal_contribution initial_stress'
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = T
    thermal_expansion_coeff = 5E-6 # this is the linear thermal expansion coefficient
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 330
  []
  [initial_strain]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = '20E6 0 0  0 20E6 0  0 0 20E6'
    eigenstrain_name = initial_stress
  []
(modules/porous_flow/examples/tutorial/11.i)

The thermal_contribution eigenstrain name has to be provided to the StressDivergenceTensors Kernels, by the way (see above). Stress is just linear elastic:

  [stress]
    type = ComputeLinearElasticStress
  []
(modules/porous_flow/examples/tutorial/11.i)

Finally, the effective fluid pressure must be computed because it is needed by the Porosity and the solid-fluid coupling, and the volumetric strain feeds into the Porosity:

  [effective_fluid_pressure_mat]
    type = PorousFlowEffectiveFluidPressure
  []
  [volumetric_strain]
    type = PorousFlowVolumetricStrain
  []
[]
(modules/porous_flow/examples/tutorial/11.i)

Executioner

For this model, an IterationAdaptiveDT Timestepper is used. This is because the dynamics at early times, particularly the thermal shock induced by sudden application of cool CO at the injection area, means small timesteps are needed, but after some time larger timesteps can be used.

[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = Newton
  end_time = 1E3
  [TimeStepper<<<{"href": "../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    dt = 1E3
    growth_factor = 1.2
    optimal_iterations = 10
  []
  nl_abs_tol = 1E-7
[]
(modules/porous_flow/examples/tutorial/11.i)

The end

Holey Dooley, you made it to the end, well done!

An animation of the results is shown in Figure 1. The porepressure front moves relatively quickly, followed by the saturation front, and then the temperature front. The solid mechanical deformations are governed mostly by the temperature. By the way, this type of dynamic 2-phase injection problem using PorousFlow has been benchmarked against analytical results with excellent agreement (and will hopefully be written into a PorousFlow Example — we are awaiting permission by funding bodies).

Figure 1: CO saturation, CO porepressure, temperature and hoop stress in the 2-phase CO injection simulation.

Postscript: the same again in 2D

As mentioned on Page 00, this problem is really an axially-symmetric problem, so may be better modelled by MOOSE using its "RZ" coordinate system. The changes to the input file are minimal. Apart from the mesh generation (discussed in Page 00) the changes are itemized below.

There only need by a disp_r Variable in place of disp_x and disp_y:

  [disp_r]
    scaling = 1E-5
(modules/porous_flow/examples/tutorial/11_2D.i)
[GlobalParams]
  displacements = 'disp_r disp_z'
  gravity = '0 0 0'
  biot_coefficient = 1.0
  PorousFlowDictator = dictator
(modules/porous_flow/examples/tutorial/11_2D.i)
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater pgas T disp_r'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
(modules/porous_flow/examples/tutorial/11_2D.i)

There are mechanical Kernels only for disp_r, and the StressDivergenceTensors Kernel is modified:

  [grad_stress_r]
    type = StressDivergenceRZTensors
    temperature = T
    variable = disp_r
    eigenstrain_names = thermal_contribution
    use_displaced_mesh = false
    component = 0
  []
(modules/porous_flow/examples/tutorial/11_2D.i)

The stress AuxKernels are modified. In SolidMechanics with RZ coordinates, the 00 component is , the 11 component is and the 22 component is . Therefore, these AuxKernels read

  [stress_rr_aux]
    type = RankTwoAux
    variable = stress_rr
    rank_two_tensor = stress
    index_i = 0
    index_j = 0
  []
  [stress_tt]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_tt
    index_i = 2
    index_j = 2
  []
  [stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 1
    index_j = 1
  []
(modules/porous_flow/examples/tutorial/11_2D.i)

The boundary conditions for the mechanics become simpler:

  [pinned_top_bottom_r]
    type = DirichletBC
    variable = disp_r
    value = 0
    boundary = 'top bottom'
  []
  [cavity_pressure_r]
    type = Pressure
    boundary = injection_area
    variable = disp_r
    postprocessor = constrained_effective_fluid_pressure_at_wellbore
    use_displaced_mesh = false
  []
(modules/porous_flow/examples/tutorial/11_2D.i)

Finally, the strain calculator needs to be of RZ type:

  [strain]
    type = ComputeAxisymmetricRZSmallStrain
    eigenstrain_names = 'thermal_contribution initial_stress'
  []
(modules/porous_flow/examples/tutorial/11_2D.i)

The reader may check that the 3D and 2D models produce the same answer, although of course the 2D version is much faster due to it having only 176 degrees of freedom compared with the 3D's 1100.

Start | Previous | Next