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 anAuxKernel
so that it will always be zero, but then it may be coupled into various MOOSE objects (using thedisplacements =
in theGlobalParams
)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 theeffective_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 noAuxKernels
.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.