Application Usage link

Physics and Syntax

Application Development link

Framework Development link

MOOSEDocs link

Infrastructure

Questions

Information and Tools

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]
  [./pwater]
    initial_condition = 20E6
  [../]
  [./pgas]
    initial_condition = 20.1E6
  [../]
  [./T]
    initial_condition = 330
    scaling = 1E-5
  [../]
  [./disp_x]
    scaling = 1E-5
  [../]
  [./disp_y]
    scaling = 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]
  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]
  [./mass_water_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./flux_water]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./mass_co2_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./flux_co2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./vol_strain_rate_co2]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./energy_dot]
    type = PorousFlowEnergyTimeDerivative
    use_displaced_mesh = false
    variable = T
  [../]
  [./advection]
    type = PorousFlowHeatAdvection
    use_displaced_mesh = false
    variable = T
  [../]
  [./conduction]
    type = PorousFlowHeatConduction
    use_displaced_mesh = false
    variable = T
  [../]
  [./vol_strain_rate_heat]
    type = PorousFlowHeatVolumetricExpansion
    use_displaced_mesh = false
    variable = T
  [../]
  [./grad_stress_x]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_x
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 0
  [../]
  [./poro_x]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_x
    use_displaced_mesh = false
    component = 0
  [../]
  [./grad_stress_y]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_y
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 1
  [../]
  [./poro_y]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_y
    use_displaced_mesh = false
    component = 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]
  [./disp_z]
  [../]
  [./effective_fluid_pressure]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./mass_frac_phase0_species0]
    initial_condition = 1 # all water in phase=0
  [../]
  [./mass_frac_phase1_species0]
    initial_condition = 0 # no water in phase=1
  [../]
  [./sgas]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./swater]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_rr]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_tt]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_zz]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./porosity]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[AuxKernels]
  [./effective_fluid_pressure]
    type = ParsedAux
    args = 'pwater pgas swater sgas'
    function = 'pwater * swater + pgas * sgas'
    variable = effective_fluid_pressure
  [../]
  [./swater]
    type = PorousFlowPropertyAux
    variable = swater
    property = saturation
    phase = 0
    execute_on = timestep_end
  [../]
  [./sgas]
    type = PorousFlowPropertyAux
    variable = sgas
    property = saturation
    phase = 1
    execute_on = timestep_end
  [../]
  [./stress_rr]
    type = RankTwoScalarAux
    variable = stress_rr
    rank_two_tensor = stress
    scalar_type = RadialStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  [../]
  [./stress_tt]
    type = RankTwoScalarAux
    variable = stress_tt
    rank_two_tensor = stress
    scalar_type = HoopStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  [../]
  [./stress_zz]
    type = RankTwoAux
    variable = stress_zz
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  [../]
  [./porosity]
    type = PorousFlowPropertyAux
    variable = porosity
    property = porosity
    execute_on = 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]
  [./constrain_effective_fluid_pressure]
    type = ParsedFunction
    vars = effective_fluid_pressure_at_wellbore
    vals = effective_fluid_pressure_at_wellbore
    value = '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
  [../]
(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:

[Modules]
  [./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_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_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'
  [../]
  [./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]
    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]
    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
    vars = effective_fluid_pressure_at_wellbore
    vals = effective_fluid_pressure_at_wellbore
    value = '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)

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

# Two-phase borehole injection problem
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [./make3D]
    input = annular
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  [../]
  [./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_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  [../]
  [./rename]
    type = RenameBlockGenerator
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    input = 'injection_area'
  [../]
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater pgas T disp_x disp_y'
    number_fluid_phases = 2
    number_fluid_components = 2
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  [../]
[]

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

[Variables]
  [./pwater]
    initial_condition = 20E6
  [../]
  [./pgas]
    initial_condition = 20.1E6
  [../]
  [./T]
    initial_condition = 330
    scaling = 1E-5
  [../]
  [./disp_x]
    scaling = 1E-5
  [../]
  [./disp_y]
    scaling = 1E-5
  [../]
[]

[Kernels]
  [./mass_water_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./flux_water]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./mass_co2_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./flux_co2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./vol_strain_rate_co2]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./energy_dot]
    type = PorousFlowEnergyTimeDerivative
    use_displaced_mesh = false
    variable = T
  [../]
  [./advection]
    type = PorousFlowHeatAdvection
    use_displaced_mesh = false
    variable = T
  [../]
  [./conduction]
    type = PorousFlowHeatConduction
    use_displaced_mesh = false
    variable = T
  [../]
  [./vol_strain_rate_heat]
    type = PorousFlowHeatVolumetricExpansion
    use_displaced_mesh = false
    variable = T
  [../]
  [./grad_stress_x]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_x
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 0
  [../]
  [./poro_x]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_x
    use_displaced_mesh = false
    component = 0
  [../]
  [./grad_stress_y]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_y
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 1
  [../]
  [./poro_y]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_y
    use_displaced_mesh = false
    component = 1
  [../]
[]

[AuxVariables]
  [./disp_z]
  [../]
  [./effective_fluid_pressure]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./mass_frac_phase0_species0]
    initial_condition = 1 # all water in phase=0
  [../]
  [./mass_frac_phase1_species0]
    initial_condition = 0 # no water in phase=1
  [../]
  [./sgas]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./swater]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_rr]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_tt]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_zz]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./porosity]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[AuxKernels]
  [./effective_fluid_pressure]
    type = ParsedAux
    args = 'pwater pgas swater sgas'
    function = 'pwater * swater + pgas * sgas'
    variable = effective_fluid_pressure
  [../]
  [./swater]
    type = PorousFlowPropertyAux
    variable = swater
    property = saturation
    phase = 0
    execute_on = timestep_end
  [../]
  [./sgas]
    type = PorousFlowPropertyAux
    variable = sgas
    property = saturation
    phase = 1
    execute_on = timestep_end
  [../]
  [./stress_rr]
    type = RankTwoScalarAux
    variable = stress_rr
    rank_two_tensor = stress
    scalar_type = RadialStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  [../]
  [./stress_tt]
    type = RankTwoScalarAux
    variable = stress_tt
    rank_two_tensor = stress
    scalar_type = HoopStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  [../]
  [./stress_zz]
    type = RankTwoAux
    variable = stress_zz
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  [../]
  [./porosity]
    type = PorousFlowPropertyAux
    variable = porosity
    property = porosity
    execute_on = timestep_end
  [../]
[]

[BCs]
  [./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
  [../]

  [./cold_co2]
    type = DirichletBC
    boundary = injection_area
    variable = T
    value = 290  # injection temperature
    use_displaced_mesh = false
  [../]
  [./constant_co2_injection]
    type = PorousFlowSink
    boundary = injection_area
    variable = pgas
    fluid_phase = 1
    flux_function = -1E-4
    use_displaced_mesh = false
  [../]

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

[Modules]
  [./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_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_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]
    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]
    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
    vars = effective_fluid_pressure_at_wellbore
    vals = effective_fluid_pressure_at_wellbore
    value = '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)

These are then Joined into std::vectors, ready to be fed into the Kernels (and other Materials, AuxKernels, etc) using the PorousFlowJoiners:

# Two-phase borehole injection problem
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [./make3D]
    input = annular
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  [../]
  [./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_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  [../]
  [./rename]
    type = RenameBlockGenerator
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    input = 'injection_area'
  [../]
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater pgas T disp_x disp_y'
    number_fluid_phases = 2
    number_fluid_components = 2
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  [../]
[]

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

[Variables]
  [./pwater]
    initial_condition = 20E6
  [../]
  [./pgas]
    initial_condition = 20.1E6
  [../]
  [./T]
    initial_condition = 330
    scaling = 1E-5
  [../]
  [./disp_x]
    scaling = 1E-5
  [../]
  [./disp_y]
    scaling = 1E-5
  [../]
[]

[Kernels]
  [./mass_water_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./flux_water]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./mass_co2_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./flux_co2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./vol_strain_rate_co2]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./energy_dot]
    type = PorousFlowEnergyTimeDerivative
    use_displaced_mesh = false
    variable = T
  [../]
  [./advection]
    type = PorousFlowHeatAdvection
    use_displaced_mesh = false
    variable = T
  [../]
  [./conduction]
    type = PorousFlowHeatConduction
    use_displaced_mesh = false
    variable = T
  [../]
  [./vol_strain_rate_heat]
    type = PorousFlowHeatVolumetricExpansion
    use_displaced_mesh = false
    variable = T
  [../]
  [./grad_stress_x]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_x
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 0
  [../]
  [./poro_x]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_x
    use_displaced_mesh = false
    component = 0
  [../]
  [./grad_stress_y]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_y
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 1
  [../]
  [./poro_y]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_y
    use_displaced_mesh = false
    component = 1
  [../]
[]

[AuxVariables]
  [./disp_z]
  [../]
  [./effective_fluid_pressure]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./mass_frac_phase0_species0]
    initial_condition = 1 # all water in phase=0
  [../]
  [./mass_frac_phase1_species0]
    initial_condition = 0 # no water in phase=1
  [../]
  [./sgas]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./swater]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_rr]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_tt]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_zz]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./porosity]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[AuxKernels]
  [./effective_fluid_pressure]
    type = ParsedAux
    args = 'pwater pgas swater sgas'
    function = 'pwater * swater + pgas * sgas'
    variable = effective_fluid_pressure
  [../]
  [./swater]
    type = PorousFlowPropertyAux
    variable = swater
    property = saturation
    phase = 0
    execute_on = timestep_end
  [../]
  [./sgas]
    type = PorousFlowPropertyAux
    variable = sgas
    property = saturation
    phase = 1
    execute_on = timestep_end
  [../]
  [./stress_rr]
    type = RankTwoScalarAux
    variable = stress_rr
    rank_two_tensor = stress
    scalar_type = RadialStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  [../]
  [./stress_tt]
    type = RankTwoScalarAux
    variable = stress_tt
    rank_two_tensor = stress
    scalar_type = HoopStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  [../]
  [./stress_zz]
    type = RankTwoAux
    variable = stress_zz
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  [../]
  [./porosity]
    type = PorousFlowPropertyAux
    variable = porosity
    property = porosity
    execute_on = timestep_end
  [../]
[]

[BCs]
  [./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
  [../]

  [./cold_co2]
    type = DirichletBC
    boundary = injection_area
    variable = T
    value = 290  # injection temperature
    use_displaced_mesh = false
  [../]
  [./constant_co2_injection]
    type = PorousFlowSink
    boundary = injection_area
    variable = pgas
    fluid_phase = 1
    flux_function = -1E-4
    use_displaced_mesh = false
  [../]

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

[Modules]
  [./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_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_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]
    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]
    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
    vars = effective_fluid_pressure_at_wellbore
    vals = effective_fluid_pressure_at_wellbore
    value = '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)

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. These are then Joined:

# Two-phase borehole injection problem
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [./make3D]
    input = annular
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  [../]
  [./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_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  [../]
  [./rename]
    type = RenameBlockGenerator
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    input = 'injection_area'
  [../]
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater pgas T disp_x disp_y'
    number_fluid_phases = 2
    number_fluid_components = 2
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  [../]
[]

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

[Variables]
  [./pwater]
    initial_condition = 20E6
  [../]
  [./pgas]
    initial_condition = 20.1E6
  [../]
  [./T]
    initial_condition = 330
    scaling = 1E-5
  [../]
  [./disp_x]
    scaling = 1E-5
  [../]
  [./disp_y]
    scaling = 1E-5
  [../]
[]

[Kernels]
  [./mass_water_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./flux_water]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./mass_co2_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./flux_co2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./vol_strain_rate_co2]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./energy_dot]
    type = PorousFlowEnergyTimeDerivative
    use_displaced_mesh = false
    variable = T
  [../]
  [./advection]
    type = PorousFlowHeatAdvection
    use_displaced_mesh = false
    variable = T
  [../]
  [./conduction]
    type = PorousFlowHeatConduction
    use_displaced_mesh = false
    variable = T
  [../]
  [./vol_strain_rate_heat]
    type = PorousFlowHeatVolumetricExpansion
    use_displaced_mesh = false
    variable = T
  [../]
  [./grad_stress_x]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_x
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 0
  [../]
  [./poro_x]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_x
    use_displaced_mesh = false
    component = 0
  [../]
  [./grad_stress_y]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_y
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 1
  [../]
  [./poro_y]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_y
    use_displaced_mesh = false
    component = 1
  [../]
[]

[AuxVariables]
  [./disp_z]
  [../]
  [./effective_fluid_pressure]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./mass_frac_phase0_species0]
    initial_condition = 1 # all water in phase=0
  [../]
  [./mass_frac_phase1_species0]
    initial_condition = 0 # no water in phase=1
  [../]
  [./sgas]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./swater]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_rr]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_tt]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_zz]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./porosity]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[AuxKernels]
  [./effective_fluid_pressure]
    type = ParsedAux
    args = 'pwater pgas swater sgas'
    function = 'pwater * swater + pgas * sgas'
    variable = effective_fluid_pressure
  [../]
  [./swater]
    type = PorousFlowPropertyAux
    variable = swater
    property = saturation
    phase = 0
    execute_on = timestep_end
  [../]
  [./sgas]
    type = PorousFlowPropertyAux
    variable = sgas
    property = saturation
    phase = 1
    execute_on = timestep_end
  [../]
  [./stress_rr]
    type = RankTwoScalarAux
    variable = stress_rr
    rank_two_tensor = stress
    scalar_type = RadialStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  [../]
  [./stress_tt]
    type = RankTwoScalarAux
    variable = stress_tt
    rank_two_tensor = stress
    scalar_type = HoopStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  [../]
  [./stress_zz]
    type = RankTwoAux
    variable = stress_zz
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  [../]
  [./porosity]
    type = PorousFlowPropertyAux
    variable = porosity
    property = porosity
    execute_on = timestep_end
  [../]
[]

[BCs]
  [./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
  [../]

  [./cold_co2]
    type = DirichletBC
    boundary = injection_area
    variable = T
    value = 290  # injection temperature
    use_displaced_mesh = false
  [../]
  [./constant_co2_injection]
    type = PorousFlowSink
    boundary = injection_area
    variable = pgas
    fluid_phase = 1
    flux_function = -1E-4
    use_displaced_mesh = false
  [../]

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

[Modules]
  [./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_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_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]
    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]
    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
    vars = effective_fluid_pressure_at_wellbore
    vals = effective_fluid_pressure_at_wellbore
    value = '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:

# Two-phase borehole injection problem
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [./make3D]
    input = annular
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  [../]
  [./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_subdomain_ids = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  [../]
  [./rename]
    type = RenameBlockGenerator
    old_block_id = '0 1'
    new_block_name = 'caps aquifer'
    input = 'injection_area'
  [../]
[]

[UserObjects]
  [./dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater pgas T disp_x disp_y'
    number_fluid_phases = 2
    number_fluid_components = 2
  [../]
  [./pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  [../]
[]

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

[Variables]
  [./pwater]
    initial_condition = 20E6
  [../]
  [./pgas]
    initial_condition = 20.1E6
  [../]
  [./T]
    initial_condition = 330
    scaling = 1E-5
  [../]
  [./disp_x]
    scaling = 1E-5
  [../]
  [./disp_y]
    scaling = 1E-5
  [../]
[]

[Kernels]
  [./mass_water_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./flux_water]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  [../]
  [./mass_co2_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./flux_co2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./vol_strain_rate_co2]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  [../]
  [./energy_dot]
    type = PorousFlowEnergyTimeDerivative
    use_displaced_mesh = false
    variable = T
  [../]
  [./advection]
    type = PorousFlowHeatAdvection
    use_displaced_mesh = false
    variable = T
  [../]
  [./conduction]
    type = PorousFlowHeatConduction
    use_displaced_mesh = false
    variable = T
  [../]
  [./vol_strain_rate_heat]
    type = PorousFlowHeatVolumetricExpansion
    use_displaced_mesh = false
    variable = T
  [../]
  [./grad_stress_x]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_x
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 0
  [../]
  [./poro_x]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_x
    use_displaced_mesh = false
    component = 0
  [../]
  [./grad_stress_y]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_y
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 1
  [../]
  [./poro_y]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_y
    use_displaced_mesh = false
    component = 1
  [../]
[]

[AuxVariables]
  [./disp_z]
  [../]
  [./effective_fluid_pressure]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./mass_frac_phase0_species0]
    initial_condition = 1 # all water in phase=0
  [../]
  [./mass_frac_phase1_species0]
    initial_condition = 0 # no water in phase=1
  [../]
  [./sgas]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./swater]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_rr]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_tt]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./stress_zz]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./porosity]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[AuxKernels]
  [./effective_fluid_pressure]
    type = ParsedAux
    args = 'pwater pgas swater sgas'
    function = 'pwater * swater + pgas * sgas'
    variable = effective_fluid_pressure
  [../]
  [./swater]
    type = PorousFlowPropertyAux
    variable = swater
    property = saturation
    phase = 0
    execute_on = timestep_end
  [../]
  [./sgas]
    type = PorousFlowPropertyAux
    variable = sgas
    property = saturation
    phase = 1
    execute_on = timestep_end
  [../]
  [./stress_rr]
    type = RankTwoScalarAux
    variable = stress_rr
    rank_two_tensor = stress
    scalar_type = RadialStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  [../]
  [./stress_tt]
    type = RankTwoScalarAux
    variable = stress_tt
    rank_two_tensor = stress
    scalar_type = HoopStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  [../]
  [./stress_zz]
    type = RankTwoAux
    variable = stress_zz
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  [../]
  [./porosity]
    type = PorousFlowPropertyAux
    variable = porosity
    property = porosity
    execute_on = timestep_end
  [../]
[]

[BCs]
  [./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
  [../]

  [./cold_co2]
    type = DirichletBC
    boundary = injection_area
    variable = T
    value = 290  # injection temperature
    use_displaced_mesh = false
  [../]
  [./constant_co2_injection]
    type = PorousFlowSink
    boundary = injection_area
    variable = pgas
    fluid_phase = 1
    flux_function = -1E-4
    use_displaced_mesh = false
  [../]

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

[Modules]
  [./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_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_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]
    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
  [../]

  [./effective_fluid_pressure]
    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
    vars = effective_fluid_pressure_at_wellbore
    vals = effective_fluid_pressure_at_wellbore
    value = '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)

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:

# Two-phase borehole injection problem
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
(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]
  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
[]
(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
    thermal_eigenstrain_name = thermal_contribution
    use_displaced_mesh = false
    component = 0
  [../]
(modules/porous_flow/examples/tutorial/11_2D.i)

The stress AuxKernels are modified. In TensorMechanics 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
    component = 0
    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