- PorousFlowDictatorThe UserObject that holds the list of PorousFlow variable names.
C++ Type:UserObjectName
Controllable:No
Description:The UserObject that holds the list of PorousFlow variable names.
PorousFlowHeatEnergy
Calculates the sum of heat energy of fluid phase(s) and/or the porous skeleton in a region
This Postprocessor calculates the heat energy of fluid phase(s) using where all variables are defined in nomenclature.
The phases that the heat energy is summed over can be entered in the phase input parameter. Multiple indices can be entered.
By default, the additional heat energy due to the porous material is added to the heat energy of the fluid phase(s). This contribution can be ignored by setting include_porous_skeleton = false.
The flag use_displaced_mesh = false is set internally by this Postprocessor, and the parameter cannot be altered by the user, even for simulations with solid-mechanical deformation. The reason is that this postprocessor uses the strain calculated by SolidMechanics to automatically compensate for deformed meshes. Further information may be found here. Therefore:
For mechanically-coupled simulations, you must provide a
base_namethat is identical to that used by the SolidMechanics strain calculator, so that this Postprocessor can retrieve the correct strain. The most common use-case is that you provide nobase_nameto the SolidMechanics strain calculator and hence nobase_nameto this Postprocessor.For non-mechanically-coupled simulations, you must not provide a
base_namethat is used in any SolidMechanics strain calculators. The most common use-case is that you have no SolidMechanics strain calculators, so you needn't worry about providing anybase_nameto this postprocessor. However, there is a possibility that you have a SolidMechanics strain calculator but you don't want to couple mechanics to PorousFlow. In that case, supplybase_name = non_existent, or similar, so that this Postprocessor doesn't retrieve any strain.
Input Parameters
- base_nameFor non-mechanically-coupled systems with no TensorMechanics strain calculators, base_name need not be set. For mechanically-coupled systems, base_name should be the same base_name as given to the TensorMechanics object that computes strain, so that this Postprocessor can correctly account for changes in mesh volume. For non-mechanically-coupled systems, base_name should not be the base_name of any TensorMechanics strain calculators.
C++ Type:std::string
Controllable:No
Description:For non-mechanically-coupled systems with no TensorMechanics strain calculators, base_name need not be set. For mechanically-coupled systems, base_name should be the same base_name as given to the TensorMechanics object that computes strain, so that this Postprocessor can correctly account for changes in mesh volume. For non-mechanically-coupled systems, base_name should not be the base_name of any TensorMechanics strain calculators.
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- include_porous_skeletonTrueInclude the heat energy of the porous skeleton
Default:True
C++ Type:bool
Controllable:No
Description:Include the heat energy of the porous skeleton
- kernel_variable_number0The PorousFlow variable number (according to the dictatory) of the heat-energy kernel. This is required only in the unusual situation where a variety of different finite-element interpolation schemes are employed in the simulation
Default:0
C++ Type:unsigned int
Controllable:No
Description:The PorousFlow variable number (according to the dictatory) of the heat-energy kernel. This is required only in the unusual situation where a variety of different finite-element interpolation schemes are employed in the simulation
- phaseThe index(es) of the fluid phase that this Postprocessor is restricted to. Multiple indices can be entered.
C++ Type:std::vector<unsigned int>
Controllable:No
Description:The index(es) of the fluid phase that this Postprocessor is restricted to. Multiple indices can be entered.
Optional Parameters
- allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
Default:False
C++ Type:bool
Controllable:No
Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
- execute_onTIMESTEP_ENDThe 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.
Default:TIMESTEP_END
C++ Type:ExecFlagEnum
Options:XFEM_MARK, FORWARD, ADJOINT, HOMOGENEOUS_FORWARD, ADJOINT_TIMESTEP_BEGIN, ADJOINT_TIMESTEP_END, NONE, INITIAL, LINEAR, LINEAR_CONVERGENCE, NONLINEAR, NONLINEAR_CONVERGENCE, POSTCHECK, TIMESTEP_END, TIMESTEP_BEGIN, MULTIAPP_FIXED_POINT_END, MULTIAPP_FIXED_POINT_BEGIN, MULTIAPP_FIXED_POINT_CONVERGENCE, FINAL, CUSTOM, TRANSFER
Controllable:No
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.
- execution_order_group0Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
Default:0
C++ Type:int
Controllable:No
Description:Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
- force_postauxFalseForces the UserObject to be executed in POSTAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in POSTAUX
- force_preauxFalseForces the UserObject to be executed in PREAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREAUX
- force_preicFalseForces the UserObject to be executed in PREIC during initial setup
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREIC during initial setup
Execution Scheduling Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- outputsVector of output names where you would like to restrict the output of variables(s) associated with this object
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
Advanced Parameters
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Material Property Retrieval Parameters
Input Files
- (modules/porous_flow/test/tests/sinks/s12.i)
- (modules/porous_flow/test/tests/energy_conservation/except02.i)
- (modules/porous_flow/test/tests/energy_conservation/heat04_fullysat_action.i)
- (modules/porous_flow/test/tests/energy_conservation/except01.i)
- (modules/porous_flow/test/tests/energy_conservation/heat04.i)
- (modules/porous_flow/test/tests/energy_conservation/except03.i)
- (modules/porous_flow/test/tests/energy_conservation/heat02.i)
- (modules/porous_flow/test/tests/energy_conservation/heat03.i)
- (modules/porous_flow/test/tests/energy_conservation/heat04_rz.i)
- (modules/porous_flow/test/tests/energy_conservation/heat04_action_KT.i)
- (modules/porous_flow/test/tests/energy_conservation/heat04_action.i)
- (modules/porous_flow/test/tests/energy_conservation/heat01.i)
- (modules/porous_flow/test/tests/fluidstate/water_vapor_phasechange.i)
- (modules/porous_flow/test/tests/energy_conservation/heat03_rz.i)
- (modules/porous_flow/test/tests/dirackernels/hfrompps.i)
(modules/porous_flow/test/tests/sinks/s12.i)
# The PorousFlowEnthalpy sink adds heat energy corresponding to injecting a 1kg/s/m^2 (flux_function = -1)
# of fluid at pressure 0.5 (given by the AuxVariable p_aux) and the input temperature is 300 (given by the T_in parameter).
# SimpleFluidProperties are used, with density0 = 10, bulk_modulus = 1, thermal_expansion = 0, and cv = 1E-4
# density = 10 * exp(0.5 / 1 + 0) = 16.4872
# internal energy = 1E-4 * 300 = 0.03
# enthalpy = 0.03 + 0.5/16.3872 = 0.0603265
# This is applied over an area of 100, so the total energy flux is 6.03265 J/s.
# This the the rate of change of the heat energy reported by the PorousFlowHeatEnergy Postprocessor
[Mesh]
type = GeneratedMesh
dim = 3
nx = 2
ny = 2
nz = 2
xmin = 0
xmax = 10
ymin = 0
ymax = 10
zmin = 0
zmax = 10
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pp temp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureConst
[]
[]
[AuxVariables]
[p_aux]
initial_condition = 0.5
[]
[]
[Variables]
[pp]
initial_condition = 1
[]
[temp]
initial_condition = 2
[]
[]
[Kernels]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[heat_conduction]
type = PorousFlowEnergyTimeDerivative
variable = temp
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1
density0 = 10
thermal_expansion = 0
cv = 1E-4
[]
[]
[Materials]
[ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = pp
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2
density = 3
[]
[]
[BCs]
[left_p]
type = PorousFlowSink
variable = pp
boundary = left
flux_function = -1
[]
[left_T]
type = PorousFlowEnthalpySink
variable = temp
boundary = left
T_in = 300
fp = simple_fluid
flux_function = -1
porepressure_var = p_aux
[]
[]
[Postprocessors]
[total_heat_energy]
type = PorousFlowHeatEnergy
phase = 0
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 0.25
num_steps = 2
[]
[Outputs]
file_base = s12
[csv]
type = CSV
[]
[]
(modules/porous_flow/test/tests/energy_conservation/except02.i)
# checking that the heat-energy postprocessor throws the correct error if the kernel_variable_number is illegal
[Mesh]
type = GeneratedMesh
dim = 1
nx = 3
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[temp]
[]
[]
[ICs]
[tinit]
type = FunctionIC
function = '100*x'
variable = temp
[]
[pinit]
type = FunctionIC
function = x
variable = pp
[]
[]
[Kernels]
[dummyt]
type = TimeDerivative
variable = temp
[]
[dummyp]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temp pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1
density0 = 1
viscosity = 0.001
thermal_expansion = 0
cv = 1.3
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2.2
density = 0.5
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[]
[Postprocessors]
[total_heat]
type = PorousFlowHeatEnergy
kernel_variable_number = 2
[]
[rock_heat]
type = PorousFlowHeatEnergy
[]
[fluid_heat]
type = PorousFlowHeatEnergy
include_porous_skeleton = false
phase = 0
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1 1 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
[Outputs]
execute_on = 'timestep_end'
file_base = except01
csv = true
[]
(modules/porous_flow/test/tests/energy_conservation/heat04_fullysat_action.i)
# heat04, but using an action
#
# The sample is a single unit element, with fixed displacements on
# all sides. A heat source of strength S (J/m^3/s) is applied into
# the element. There is no fluid flow or heat flow. The rise
# in temperature, porepressure and stress, and the change in porosity is
# matched with theory.
#
# In this case, fluid mass must be conserved, and there is no
# volumetric strain, so
# porosity * fluid_density = constant
# Also, the energy-density in the rock-fluid system increases with S:
# d/dt [(1 - porosity) * rock_density * rock_heat_cap * T + porosity * fluid_density * fluid_heat_cap * T] = S
# Also, the porosity evolves according to THM as
# porosity = biot + (porosity0 - biot) * exp( (biot - 1) * P / fluid_bulk + rock_thermal_exp * T)
# Finally, the effective stress must be exactly zero (as there is
# no strain).
#
# Let us assume that
# fluid_density = dens0 * exp(P / fluid_bulk - fluid_thermal_exp * T)
# Then the conservation of fluid mass means
# porosity = por0 * exp(- P / fluid_bulk + fluid_thermal_exp * T)
# where dens0 * por0 = the initial fluid mass.
# The last expression for porosity, combined with the THM one,
# and assuming that biot = 1 for simplicity, gives
# porosity = 1 + (porosity0 - 1) * exp(rock_thermal_exp * T) = por0 * exp(- P / fluid_bulk + fluid_thermal_exp * T) .... (A)
#
# This stuff may be substituted into the heat energy-density equation:
# S = d/dt [(1 - porosity0) * exp(rock_thermal_exp * T) * rock_density * rock_heat_cap * T + porosity * fluid_density * fluid_heat_cap * T]
#
# If S is constant then
# S * t = (1 - porosity0) * exp(rock_thermal_exp * T) * rock_density * rock_heat_cap * T + porosity * fluid_density * fluid_heat_cap * T
# with T(t=0) = 0 then Eqn(A) implies that por0 = porosity0 and
# P / fluid_bulk = fluid_thermal_exp * T - log(1 + (por0 - 1) * exp(rock_thermal_exp * T)) + log(por0)
#
# Parameters:
# A = 2
# fluid_bulk = 2.0
# dens0 = 3.0
# fluid_thermal_exp = 0.5
# fluid_heat_cap = 2
# por0 = 0.5
# rock_thermal_exp = 0.25
# rock_density = 5
# rock_heat_capacity = 0.2
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.5
cv = 2
cp = 2
bulk_modulus = 2.0
density0 = 3.0
[]
[]
[PorousFlowFullySaturated]
coupling_type = ThermoHydroMechanical
displacements = 'disp_x disp_y disp_z'
porepressure = pp
temperature = temp
dictator_name = Sir
biot_coefficient = 1.0
gravity = '0 0 0'
fp = the_simple_fluid
stabilization = none
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = Sir
block = 0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[pp]
[]
[temp]
[]
[]
[BCs]
[confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[]
[confinez]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[]
[]
[Kernels]
[heat_source]
type = BodyForce
function = 1
variable = temp
[]
[]
[Functions]
[err_T_fcn]
type = ParsedFunction
symbol_names = 'por0 rte temp rd rhc m0 fhc source'
symbol_values = '0.5 0.25 t0 5 0.2 1.5 2 1'
expression = '((1-por0)*exp(rte*temp)*rd*rhc*temp+m0*fhc*temp-source*t)/(source*t)'
[]
[err_pp_fcn]
type = ParsedFunction
symbol_names = 'por0 rte temp rd rhc m0 fhc source bulk pp fte'
symbol_values = '0.5 0.25 t0 5 0.2 1.5 2 1 2 p0 0.5'
expression = '(bulk*(fte*temp-log(1+(por0-1)*exp(rte*temp))+log(por0))-pp)/pp'
[]
[]
[AuxVariables]
[porosity]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[porosity]
type = PorousFlowPropertyAux
property = porosity
variable = porosity
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosity
thermal = true
fluid = true
mechanical = true
ensure_positive = false
biot_coefficient = 1.0
porosity_zero = 0.5
thermal_expansion_coeff = 0.25
solid_bulk = 2
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 0.2
density = 5.0
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '0 0 0 0 0 0 0 0 0'
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '0 0 0 0 0 0 0 0 0'
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = pp
[]
[t0]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = temp
[]
[porosity]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = porosity
[]
[stress_xx]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_xx
[]
[stress_yy]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_yy
[]
[stress_zz]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_zz
[]
[fluid_mass]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'timestep_end'
outputs = 'console csv'
[]
[total_heat]
type = PorousFlowHeatEnergy
phase = 0
execute_on = 'timestep_end'
outputs = 'console csv'
[]
[err_T]
type = FunctionValuePostprocessor
function = err_T_fcn
[]
[err_P]
type = FunctionValuePostprocessor
function = err_pp_fcn
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-12 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 5
[]
[Outputs]
execute_on = 'initial timestep_end'
file_base = heat04_fullysat_action
csv = true
[]
(modules/porous_flow/test/tests/energy_conservation/except01.i)
# checking that the heat-energy postprocessor throws the correct error if the phase number is entered incorrectly
[Mesh]
type = GeneratedMesh
dim = 1
nx = 3
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[temp]
[]
[]
[ICs]
[tinit]
type = FunctionIC
function = '100*x'
variable = temp
[]
[pinit]
type = FunctionIC
function = x
variable = pp
[]
[]
[Kernels]
[dummyt]
type = TimeDerivative
variable = temp
[]
[dummyp]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temp pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1
density0 = 1
viscosity = 0.001
thermal_expansion = 0
cv = 1.3
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2.2
density = 0.5
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[]
[Postprocessors]
[total_heat]
type = PorousFlowHeatEnergy
phase = 1
[]
[rock_heat]
type = PorousFlowHeatEnergy
[]
[fluid_heat]
type = PorousFlowHeatEnergy
include_porous_skeleton = false
phase = 0
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1 1 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
[Outputs]
execute_on = 'timestep_end'
file_base = except01
csv = true
[]
(modules/porous_flow/test/tests/energy_conservation/heat04.i)
# The sample is a single unit element, with fixed displacements on
# all sides. A heat source of strength S (J/m^3/s) is applied into
# the element. There is no fluid flow or heat flow. The rise
# in temperature, porepressure and stress, and the change in porosity is
# matched with theory.
#
# In this case, fluid mass must be conserved, and there is no
# volumetric strain, so
# porosity * fluid_density = constant
# Also, the energy-density in the rock-fluid system increases with S:
# d/dt [(1 - porosity) * rock_density * rock_heat_cap * T + porosity * fluid_density * fluid_heat_cap * T] = S
# Also, the porosity evolves according to THM as
# porosity = biot + (porosity0 - biot) * exp( (biot - 1) * P / fluid_bulk + rock_thermal_exp * T)
# Finally, the effective stress must be exactly zero (as there is
# no strain).
#
# Let us assume that
# fluid_density = dens0 * exp(P / fluid_bulk - fluid_thermal_exp * T)
# Then the conservation of fluid mass means
# porosity = por0 * exp(- P / fluid_bulk + fluid_thermal_exp * T)
# where dens0 * por0 = the initial fluid mass.
# The last expression for porosity, combined with the THM one,
# and assuming that biot = 1 for simplicity, gives
# porosity = 1 + (porosity0 - 1) * exp(rock_thermal_exp * T) = por0 * exp(- P / fluid_bulk + fluid_thermal_exp * T) .... (A)
#
# This stuff may be substituted into the heat energy-density equation:
# S = d/dt [(1 - porosity0) * exp(rock_thermal_exp * T) * rock_density * rock_heat_cap * T + porosity * fluid_density * fluid_heat_cap * T]
#
# If S is constant then
# S * t = (1 - porosity0) * exp(rock_thermal_exp * T) * rock_density * rock_heat_cap * T + porosity * fluid_density * fluid_heat_cap * T
# with T(t=0) = 0 then Eqn(A) implies that por0 = porosity0 and
# P / fluid_bulk = fluid_thermal_exp * T - log(1 + (por0 - 1) * exp(rock_thermal_exp * T)) + log(por0)
#
# Parameters:
# A = 2
# fluid_bulk = 2.0
# dens0 = 3.0
# fluid_thermal_exp = 0.5
# fluid_heat_cap = 2
# por0 = 0.5
# rock_thermal_exp = 0.25
# rock_density = 5
# rock_heat_capacity = 0.2
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.5
cv = 2
cp = 2
bulk_modulus = 2.0
density0 = 3.0
[]
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[pp]
[]
[temp]
[]
[]
[BCs]
[confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[]
[confinez]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[]
[]
[Kernels]
[grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[]
[grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[]
[grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[]
[poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 1.0
variable = disp_x
component = 0
[]
[poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 1.0
variable = disp_y
component = 1
[]
[poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 1.0
component = 2
variable = disp_z
[]
[poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
variable = pp
fluid_component = 0
[]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[temp]
type = PorousFlowEnergyTimeDerivative
variable = temp
[]
[poro_vol_exp_temp]
type = PorousFlowHeatVolumetricExpansion
variable = temp
[]
[heat_source]
type = BodyForce
function = 1
variable = temp
[]
[]
[Functions]
[err_T_fcn]
type = ParsedFunction
symbol_names = 'por0 rte temp rd rhc m0 fhc source'
symbol_values = '0.5 0.25 t0 5 0.2 1.5 2 1'
expression = '((1-por0)*exp(rte*temp)*rd*rhc*temp+m0*fhc*temp-source*t)/(source*t)'
[]
[err_pp_fcn]
type = ParsedFunction
symbol_names = 'por0 rte temp rd rhc m0 fhc source bulk pp fte'
symbol_values = '0.5 0.25 t0 5 0.2 1.5 2 1 2 p0 0.5'
expression = '(bulk*(fte*temp-log(1+(por0-1)*exp(rte*temp))+log(por0))-pp)/pp'
[]
[]
[AuxVariables]
[stress_xx]
order = CONSTANT
family = MONOMIAL
[]
[stress_xy]
order = CONSTANT
family = MONOMIAL
[]
[stress_xz]
order = CONSTANT
family = MONOMIAL
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[stress_yz]
order = CONSTANT
family = MONOMIAL
[]
[stress_zz]
order = CONSTANT
family = MONOMIAL
[]
[porosity]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
[]
[stress_xy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xy
index_i = 0
index_j = 1
[]
[stress_xz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xz
index_i = 0
index_j = 2
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[]
[stress_yz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yz
index_i = 1
index_j = 2
[]
[stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[]
[porosity]
type = PorousFlowPropertyAux
property = porosity
variable = porosity
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temp pp disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[porosity]
type = PorousFlowPorosity
thermal = true
fluid = true
mechanical = true
ensure_positive = false
biot_coefficient = 1.0
porosity_zero = 0.5
thermal_expansion_coeff = 0.25
solid_bulk = 2
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 0.2
density = 5.0
[]
[ppss]
type = PorousFlow1PhaseFullySaturated
porepressure = pp
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
temperature_unit = Kelvin
fp = the_simple_fluid
phase = 0
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = pp
[]
[t0]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = temp
[]
[porosity]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = porosity
[]
[stress_xx]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_xx
[]
[stress_yy]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_yy
[]
[stress_zz]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_zz
[]
[fluid_mass]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'timestep_end'
outputs = 'console csv'
[]
[total_heat]
type = PorousFlowHeatEnergy
phase = 0
execute_on = 'timestep_end'
outputs = 'console csv'
[]
[err_T]
type = FunctionValuePostprocessor
function = err_T_fcn
[]
[err_P]
type = FunctionValuePostprocessor
function = err_pp_fcn
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-12 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 5
[]
[Outputs]
execute_on = 'initial timestep_end'
file_base = heat04
[csv]
type = CSV
[]
[]
(modules/porous_flow/test/tests/energy_conservation/except03.i)
# Checking that the heat energy postprocessor correctly throws a paramError when an incorrect
# strain base_name is given
[Mesh]
type = GeneratedMesh
dim = 1
nx = 3
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pp]
[]
[temp]
[]
[]
[ICs]
[tinit]
type = FunctionIC
function = '100*x'
variable = temp
[]
[pinit]
type = FunctionIC
function = x
variable = pp
[]
[]
[Kernels]
[dummyt]
type = TimeDerivative
variable = temp
[]
[dummyp]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temp pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1
density0 = 1
viscosity = 0.001
thermal_expansion = 0
cv = 1.3
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2.2
density = 0.5
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[]
[Postprocessors]
[heat]
type = PorousFlowHeatEnergy
base_name = incorrect_base_name
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
[Outputs]
execute_on = 'timestep_end'
[]
(modules/porous_flow/test/tests/energy_conservation/heat02.i)
# checking that the heat-energy postprocessor correctly calculates the energy
# 1phase, constant porosity
[Mesh]
type = GeneratedMesh
dim = 1
nx = 3
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[temp]
[]
[pp]
[]
[]
[ICs]
[tinit]
type = FunctionIC
function = '100*x'
variable = temp
[]
[pinit]
type = FunctionIC
function = 'x'
variable = pp
[]
[]
[Kernels]
[dummyt]
type = TimeDerivative
variable = temp
[]
[dummyp]
type = TimeDerivative
variable = pp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temp pp'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 1
density0 = 1
viscosity = 0.001
thermal_expansion = 0
cv = 1.3
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2.2
density = 0.5
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[]
[Postprocessors]
[total_heat]
type = PorousFlowHeatEnergy
phase = 0
[]
[rock_heat]
type = PorousFlowHeatEnergy
[]
[fluid_heat]
type = PorousFlowHeatEnergy
include_porous_skeleton = false
phase = 0
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1 .999 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
[Outputs]
execute_on = 'timestep_end'
file_base = heat02
csv = true
[]
(modules/porous_flow/test/tests/energy_conservation/heat03.i)
# The sample is a single unit element, with roller BCs on the sides
# and bottom. A constant displacement is applied to the top: disp_z = -0.01*t.
# There is no fluid flow or heat flow.
# Heat energy conservation is checked.
#
# Under these conditions (here L is the height of the sample: L=1 in this case):
# porepressure = porepressure(t=0) - (Fluid bulk modulus)*log(1 - 0.01*t)
# stress_xx = (bulk - 2*shear/3)*disp_z/L (remember this is effective stress)
# stress_zz = (bulk + 4*shear/3)*disp_z/L (remember this is effective stress)
# Also, the total heat energy must be conserved: this is
# fluid_mass * fluid_heat_cap * temperature + (1 - porosity) * rock_density * rock_heat_cap * temperature * volume
# Since fluid_mass is conserved, and volume = (1 - 0.01*t), this can be solved for temperature:
# temperature = initial_heat_energy / (fluid_mass * fluid_heat_cap + (1 - porosity) * rock_density * rock_heat_cap * (1 - 0.01*t))
#
# Parameters:
# Bulk modulus = 2
# Shear modulus = 1.5
# fluid bulk modulus = 0.5
# initial porepressure = 0.1
# initial temperature = 10
#
# Desired output:
# zdisp = -0.01*t
# p0 = 0.1 - 0.5*log(1-0.01*t)
# stress_xx = stress_yy = -0.01*t
# stress_zz = -0.04*t
# t0 = 11.5 / (0.159 + 0.99 * (1 - 0.01*t))
#
# Regarding the "log" - it comes from preserving fluid mass
#
# Note that the PorousFlowMassVolumetricExpansion and PorousFlowHeatVolumetricExpansion Kernels are used
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = dictator
block = 0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[pp]
initial_condition = 0.1
[]
[temp]
initial_condition = 10
[]
[]
[BCs]
[confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[]
[basefixed]
type = DirichletBC
variable = disp_z
value = 0
boundary = back
[]
[top_velocity]
type = FunctionDirichletBC
variable = disp_z
function = -0.01*t
boundary = front
[]
[]
[Kernels]
[grad_stress_x]
type = StressDivergenceTensors
variable = disp_x
component = 0
[]
[grad_stress_y]
type = StressDivergenceTensors
variable = disp_y
component = 1
[]
[grad_stress_z]
type = StressDivergenceTensors
variable = disp_z
component = 2
[]
[poro_x]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.3
variable = disp_x
component = 0
[]
[poro_y]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.3
variable = disp_y
component = 1
[]
[poro_z]
type = PorousFlowEffectiveStressCoupling
biot_coefficient = 0.3
component = 2
variable = disp_z
[]
[poro_vol_exp]
type = PorousFlowMassVolumetricExpansion
variable = pp
fluid_component = 0
[]
[mass0]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pp
[]
[temp]
type = PorousFlowEnergyTimeDerivative
variable = temp
[]
[poro_vol_exp_temp]
type = PorousFlowHeatVolumetricExpansion
variable = temp
[]
[]
[AuxVariables]
[stress_xx]
order = CONSTANT
family = MONOMIAL
[]
[stress_xy]
order = CONSTANT
family = MONOMIAL
[]
[stress_xz]
order = CONSTANT
family = MONOMIAL
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[stress_yz]
order = CONSTANT
family = MONOMIAL
[]
[stress_zz]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_xx]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xx
index_i = 0
index_j = 0
[]
[stress_xy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xy
index_i = 0
index_j = 1
[]
[stress_xz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_xz
index_i = 0
index_j = 2
[]
[stress_yy]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yy
index_i = 1
index_j = 1
[]
[stress_yz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_yz
index_i = 1
index_j = 2
[]
[stress_zz]
type = RankTwoAux
rank_two_tensor = stress
variable = stress_zz
index_i = 2
index_j = 2
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temp pp disp_x disp_y disp_z'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureVG
m = 0.5
alpha = 1
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 0.5
density0 = 1
viscosity = 1
thermal_expansion = 0
cv = 1.3
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[vol_strain]
type = PorousFlowVolumetricStrain
[]
[eff_fluid_pressure]
type = PorousFlowEffectiveFluidPressure
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2.2
density = 0.5
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pp
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
[]
[simple_fluid]
type = PorousFlowSingleComponentFluid
fp = simple_fluid
phase = 0
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '0.5 0 0 0 0.5 0 0 0 0.5'
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = 'console csv'
execute_on = 'initial timestep_end'
point = '0 0 0'
variable = pp
[]
[t0]
type = PointValue
outputs = 'console csv'
execute_on = 'initial timestep_end'
point = '0 0 0'
variable = temp
[]
[zdisp]
type = PointValue
outputs = csv
point = '0 0 0.5'
use_displaced_mesh = false
variable = disp_z
[]
[stress_xx]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_xx
[]
[stress_yy]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_yy
[]
[stress_zz]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_zz
[]
[fluid_mass]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[total_heat]
type = PorousFlowHeatEnergy
phase = 0
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[rock_heat]
type = PorousFlowHeatEnergy
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[fluid_heat]
type = PorousFlowHeatEnergy
include_porous_skeleton = false
phase = 0
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-14 1E-8 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 2
end_time = 10
[]
[Outputs]
execute_on = 'initial timestep_end'
file_base = heat03
[csv]
type = CSV
[]
[]
(modules/porous_flow/test/tests/energy_conservation/heat04_rz.i)
# The sample is a single unit element in RZ coordinates
# A constant velocity is applied to the outer boundary is free to move as a source injects heat and fluid into the system
# There is no fluid flow or heat flow.
# Heat energy conservation is checked.
# Mass conservation is checked
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 1
xmin = 1
xmax = 2
ymin = -0.5
ymax = 0.5
coord_type = RZ
[]
[GlobalParams]
displacements = 'disp_r disp_z'
PorousFlowDictator = dictator
block = 0
biot_coefficient = 0.3
[]
[Variables]
[disp_r]
[]
[disp_z]
[]
[pp]
initial_condition = 0.1
[]
[temp]
initial_condition = 10
[]
[]
[BCs]
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'bottom top'
[]
[rmin_fixed]
type = DirichletBC
variable = disp_r
value = 0
boundary = left
[]
[contract]
type = FunctionDirichletBC
variable = disp_r
function = -0.01*t
boundary = right
[]
[]
[PorousFlowFullySaturated]
coupling_type = ThermoHydroMechanical
porepressure = pp
temperature = temp
fp = simple_fluid
[]
[DiracKernels]
[heat_source]
type = PorousFlowPointSourceFromPostprocessor
point = '1.5 0 0'
variable = temp
mass_flux = 10
[]
[fluid_source]
type = PorousFlowPointSourceFromPostprocessor
point = '1.5 0 0'
variable = pp
mass_flux = 1
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 0.5
density0 = 1
viscosity = 1
thermal_expansion = 0
cv = 1.3
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeAxisymmetricRZSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2.2
density = 0.5
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '0.5 0 0 0 0.5 0 0 0 0.5'
[]
[thermal_cond]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '1 0 0 0 1 0 0 0 1'
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = 'console csv'
execute_on = 'initial timestep_end'
point = '1 0 0'
variable = pp
[]
[t0]
type = PointValue
outputs = 'console csv'
execute_on = 'initial timestep_end'
point = '1 0 0'
variable = temp
[]
[rdisp]
type = PointValue
outputs = 'csv console'
point = '2 0 0'
use_displaced_mesh = false
variable = disp_r
[]
[fluid_mass]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[total_heat]
type = PorousFlowHeatEnergy
phase = 0
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[rock_heat]
type = PorousFlowHeatEnergy
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[fluid_heat]
type = PorousFlowHeatEnergy
include_porous_skeleton = false
phase = 0
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 2
end_time = 10
[]
[Outputs]
execute_on = 'initial timestep_end'
[csv]
type = CSV
[]
[]
(modules/porous_flow/test/tests/energy_conservation/heat04_action_KT.i)
# heat04, but using an action with KT stabilization.
# See heat04.i for a full discussion of the results.
# The KT stabilization should have no impact as there is no flow, but this input file checks that MOOSE runs.
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.5
cv = 2
cp = 2
bulk_modulus = 2.0
density0 = 3.0
[]
[]
[PorousFlowUnsaturated]
coupling_type = ThermoHydroMechanical
displacements = 'disp_x disp_y disp_z'
porepressure = pp
temperature = temp
dictator_name = Sir
biot_coefficient = 1.0
gravity = '0 0 0'
fp = the_simple_fluid
van_genuchten_alpha = 1.0E-12
van_genuchten_m = 0.5
relative_permeability_type = Corey
relative_permeability_exponent = 0.0
stabilization = KT
flux_limiter_type = superbee
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = Sir
block = 0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[pp]
[]
[temp]
[]
[]
[BCs]
[confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[]
[confinez]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[]
[]
[Kernels]
[heat_source]
type = BodyForce
function = 1
variable = temp
[]
[]
[Functions]
[err_T_fcn]
type = ParsedFunction
symbol_names = 'por0 rte temp rd rhc m0 fhc source'
symbol_values = '0.5 0.25 t0 5 0.2 1.5 2 1'
expression = '((1-por0)*exp(rte*temp)*rd*rhc*temp+m0*fhc*temp-source*t)/(source*t)'
[]
[err_pp_fcn]
type = ParsedFunction
symbol_names = 'por0 rte temp rd rhc m0 fhc source bulk pp fte'
symbol_values = '0.5 0.25 t0 5 0.2 1.5 2 1 2 p0 0.5'
expression = '(bulk*(fte*temp-log(1+(por0-1)*exp(rte*temp))+log(por0))-pp)/pp'
[]
[]
[AuxVariables]
[porosity]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[porosity]
type = PorousFlowPropertyAux
property = porosity
variable = porosity
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosity
thermal = true
fluid = true
mechanical = true
ensure_positive = false
biot_coefficient = 1.0
porosity_zero = 0.5
thermal_expansion_coeff = 0.25
solid_bulk = 2
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 0.2
density = 5.0
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '0 0 0 0 0 0 0 0 0'
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '0 0 0 0 0 0 0 0 0'
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = pp
[]
[t0]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = temp
[]
[porosity]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = porosity
[]
[stress_xx]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_xx
[]
[stress_yy]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_yy
[]
[stress_zz]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_zz
[]
[fluid_mass]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'timestep_end'
outputs = 'console csv'
[]
[total_heat]
type = PorousFlowHeatEnergy
phase = 0
execute_on = 'timestep_end'
outputs = 'console csv'
[]
[err_T]
type = FunctionValuePostprocessor
function = err_T_fcn
[]
[err_P]
type = FunctionValuePostprocessor
function = err_pp_fcn
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-12 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 5
[]
[Outputs]
execute_on = 'initial timestep_end'
file_base = heat04_action
csv = true
[]
(modules/porous_flow/test/tests/energy_conservation/heat04_action.i)
# heat04, but using an action
#
# The sample is a single unit element, with fixed displacements on
# all sides. A heat source of strength S (J/m^3/s) is applied into
# the element. There is no fluid flow or heat flow. The rise
# in temperature, porepressure and stress, and the change in porosity is
# matched with theory.
#
# In this case, fluid mass must be conserved, and there is no
# volumetric strain, so
# porosity * fluid_density = constant
# Also, the energy-density in the rock-fluid system increases with S:
# d/dt [(1 - porosity) * rock_density * rock_heat_cap * T + porosity * fluid_density * fluid_heat_cap * T] = S
# Also, the porosity evolves according to THM as
# porosity = biot + (porosity0 - biot) * exp( (biot - 1) * P / fluid_bulk + rock_thermal_exp * T)
# Finally, the effective stress must be exactly zero (as there is
# no strain).
#
# Let us assume that
# fluid_density = dens0 * exp(P / fluid_bulk - fluid_thermal_exp * T)
# Then the conservation of fluid mass means
# porosity = por0 * exp(- P / fluid_bulk + fluid_thermal_exp * T)
# where dens0 * por0 = the initial fluid mass.
# The last expression for porosity, combined with the THM one,
# and assuming that biot = 1 for simplicity, gives
# porosity = 1 + (porosity0 - 1) * exp(rock_thermal_exp * T) = por0 * exp(- P / fluid_bulk + fluid_thermal_exp * T) .... (A)
#
# This stuff may be substituted into the heat energy-density equation:
# S = d/dt [(1 - porosity0) * exp(rock_thermal_exp * T) * rock_density * rock_heat_cap * T + porosity * fluid_density * fluid_heat_cap * T]
#
# If S is constant then
# S * t = (1 - porosity0) * exp(rock_thermal_exp * T) * rock_density * rock_heat_cap * T + porosity * fluid_density * fluid_heat_cap * T
# with T(t=0) = 0 then Eqn(A) implies that por0 = porosity0 and
# P / fluid_bulk = fluid_thermal_exp * T - log(1 + (por0 - 1) * exp(rock_thermal_exp * T)) + log(por0)
#
# Parameters:
# A = 2
# fluid_bulk = 2.0
# dens0 = 3.0
# fluid_thermal_exp = 0.5
# fluid_heat_cap = 2
# por0 = 0.5
# rock_thermal_exp = 0.25
# rock_density = 5
# rock_heat_capacity = 0.2
[Mesh]
type = GeneratedMesh
dim = 3
nx = 1
ny = 1
nz = 1
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
[]
[FluidProperties]
[the_simple_fluid]
type = SimpleFluidProperties
thermal_expansion = 0.5
cv = 2
cp = 2
bulk_modulus = 2.0
density0 = 3.0
[]
[]
[PorousFlowUnsaturated]
coupling_type = ThermoHydroMechanical
displacements = 'disp_x disp_y disp_z'
porepressure = pp
temperature = temp
dictator_name = Sir
biot_coefficient = 1.0
gravity = '0 0 0'
fp = the_simple_fluid
van_genuchten_alpha = 1.0E-12
van_genuchten_m = 0.5
relative_permeability_type = Corey
relative_permeability_exponent = 0.0
[]
[GlobalParams]
displacements = 'disp_x disp_y disp_z'
PorousFlowDictator = Sir
block = 0
[]
[Variables]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[pp]
[]
[temp]
[]
[]
[BCs]
[confinex]
type = DirichletBC
variable = disp_x
value = 0
boundary = 'left right'
[]
[confiney]
type = DirichletBC
variable = disp_y
value = 0
boundary = 'bottom top'
[]
[confinez]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'back front'
[]
[]
[Kernels]
[heat_source]
type = BodyForce
function = 1
variable = temp
[]
[]
[Functions]
[err_T_fcn]
type = ParsedFunction
symbol_names = 'por0 rte temp rd rhc m0 fhc source'
symbol_values = '0.5 0.25 t0 5 0.2 1.5 2 1'
expression = '((1-por0)*exp(rte*temp)*rd*rhc*temp+m0*fhc*temp-source*t)/(source*t)'
[]
[err_pp_fcn]
type = ParsedFunction
symbol_names = 'por0 rte temp rd rhc m0 fhc source bulk pp fte'
symbol_values = '0.5 0.25 t0 5 0.2 1.5 2 1 2 p0 0.5'
expression = '(bulk*(fte*temp-log(1+(por0-1)*exp(rte*temp))+log(por0))-pp)/pp'
[]
[]
[AuxVariables]
[porosity]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[porosity]
type = PorousFlowPropertyAux
property = porosity
variable = porosity
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosity
thermal = true
fluid = true
mechanical = true
ensure_positive = false
biot_coefficient = 1.0
porosity_zero = 0.5
thermal_expansion_coeff = 0.25
solid_bulk = 2
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 0.2
density = 5.0
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '0 0 0 0 0 0 0 0 0'
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '0 0 0 0 0 0 0 0 0'
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = pp
[]
[t0]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = temp
[]
[porosity]
type = PointValue
outputs = 'console csv'
execute_on = 'timestep_end'
point = '0 0 0'
variable = porosity
[]
[stress_xx]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_xx
[]
[stress_yy]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_yy
[]
[stress_zz]
type = PointValue
outputs = csv
point = '0 0 0'
variable = stress_zz
[]
[fluid_mass]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'timestep_end'
outputs = 'console csv'
[]
[total_heat]
type = PorousFlowHeatEnergy
phase = 0
execute_on = 'timestep_end'
outputs = 'console csv'
[]
[err_T]
type = FunctionValuePostprocessor
function = err_T_fcn
[]
[err_P]
type = FunctionValuePostprocessor
function = err_pp_fcn
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1E-12 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 5
[]
[Outputs]
execute_on = 'initial timestep_end'
file_base = heat04_action
csv = true
[]
(modules/porous_flow/test/tests/energy_conservation/heat01.i)
# checking that the heat-energy postprocessor correctly calculates the energy
# 0phase, constant porosity
[Mesh]
type = GeneratedMesh
dim = 1
nx = 3
xmin = 0
xmax = 1
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[temp]
[]
[]
[ICs]
[tinit]
type = FunctionIC
function = '100*x'
variable = temp
[]
[]
[Kernels]
[dummy]
type = TimeDerivative
variable = temp
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'temp'
number_fluid_phases = 0
number_fluid_components = 0
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temp
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2.2
density = 0.5
[]
[]
[Postprocessors]
[total_heat]
type = PorousFlowHeatEnergy
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
petsc_options_iname = '-ksp_type -pc_type -snes_atol -snes_rtol -snes_max_it'
petsc_options_value = 'bcgs bjacobi 1 .999 10000'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 1
end_time = 1
[]
[Outputs]
execute_on = 'timestep_end'
file_base = heat01
csv = true
[]
(modules/porous_flow/test/tests/fluidstate/water_vapor_phasechange.i)
# Tests correct calculation of properties in PorousFlowWaterVapor as a phase change
# from liquid to a two-phase model occurs due to a pressure drop.
# A single 10 m^3 element is used, with constant mass and heat production using
# a Dirac kernel. Initial conditions correspond to just outside the two-phase region in
# the liquid state.
#
# An identical problem can be run using TOUGH2, with the following outputs after 1,000s
# Pressure: 8.58 Mpa
# Temperature: 299.92 K
# Vapor saturation: 0.00637
[Mesh]
type = GeneratedMesh
dim = 3
xmax = 10
ymax = 10
zmax = 10
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pliq]
initial_condition = 9e6
[]
[h]
scaling = 1e-3
[]
[]
[ICs]
[hic]
type = PorousFlowFluidPropertyIC
variable = h
porepressure = pliq
property = enthalpy
temperature = 300
temperature_unit = Celsius
fp = water
[]
[]
[DiracKernels]
[mass]
type = ConstantPointSource
point = '5 5 5'
variable = pliq
value = -1
[]
[heat]
type = ConstantPointSource
point = '5 5 5'
variable = h
value = -1.344269e6
[]
[]
[AuxVariables]
[pressure_gas]
order = CONSTANT
family = MONOMIAL
[]
[pressure_water]
order = CONSTANT
family = MONOMIAL
[]
[enthalpy_gas]
order = CONSTANT
family = MONOMIAL
[]
[enthalpy_water]
order = CONSTANT
family = MONOMIAL
[]
[saturation_gas]
order = CONSTANT
family = MONOMIAL
[]
[saturation_water]
order = CONSTANT
family = MONOMIAL
[]
[density_water]
order = CONSTANT
family = MONOMIAL
[]
[density_gas]
order = CONSTANT
family = MONOMIAL
[]
[viscosity_water]
order = CONSTANT
family = MONOMIAL
[]
[viscosity_gas]
order = CONSTANT
family = MONOMIAL
[]
[temperature]
order = CONSTANT
family = MONOMIAL
[]
[e_gas]
order = CONSTANT
family = MONOMIAL
[]
[e_water]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[enthalpy_water]
type = PorousFlowPropertyAux
variable = enthalpy_water
property = enthalpy
phase = 0
execute_on = 'initial timestep_end'
[]
[enthalpy_gas]
type = PorousFlowPropertyAux
variable = enthalpy_gas
property = enthalpy
phase = 1
execute_on = 'initial timestep_end'
[]
[pressure_water]
type = PorousFlowPropertyAux
variable = pressure_water
property = pressure
phase = 0
execute_on = 'initial timestep_end'
[]
[pressure_gas]
type = PorousFlowPropertyAux
variable = pressure_gas
property = pressure
phase = 1
execute_on = 'initial timestep_end'
[]
[saturation_water]
type = PorousFlowPropertyAux
variable = saturation_water
property = saturation
phase = 0
execute_on = 'initial timestep_end'
[]
[saturation_gas]
type = PorousFlowPropertyAux
variable = saturation_gas
property = saturation
phase = 1
execute_on = 'initial timestep_end'
[]
[density_water]
type = PorousFlowPropertyAux
variable = density_water
property = density
phase = 0
execute_on = 'initial timestep_end'
[]
[density_gas]
type = PorousFlowPropertyAux
variable = density_gas
property = density
phase = 1
execute_on = 'initial timestep_end'
[]
[viscosity_water]
type = PorousFlowPropertyAux
variable = viscosity_water
property = viscosity
phase = 0
execute_on = 'initial timestep_end'
[]
[viscosity_gas]
type = PorousFlowPropertyAux
variable = viscosity_gas
property = viscosity
phase = 1
execute_on = 'initial timestep_end'
[]
[temperature]
type = PorousFlowPropertyAux
variable = temperature
property = temperature
execute_on = 'initial timestep_end'
[]
[e_water]
type = PorousFlowPropertyAux
variable = e_water
property = internal_energy
phase = 0
execute_on = 'initial timestep_end'
[]
[egas]
type = PorousFlowPropertyAux
variable = e_gas
property = internal_energy
phase = 1
execute_on = 'initial timestep_end'
[]
[]
[Kernels]
[mass]
type = PorousFlowMassTimeDerivative
variable = pliq
[]
[heat]
type = PorousFlowEnergyTimeDerivative
variable = h
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pliq h'
number_fluid_phases = 2
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureBC
pe = 1e5
lambda = 2
pc_max = 1e6
[]
[fs]
type = PorousFlowWaterVapor
water_fp = water
capillary_pressure = pc
[]
[]
[FluidProperties]
[water]
type = Water97FluidProperties
[]
[]
[Materials]
[watervapor]
type = PorousFlowFluidStateSingleComponent
porepressure = pliq
enthalpy = h
temperature_unit = Celsius
capillary_pressure = pc
fluid_state = fs
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1e-14 0 0 0 1e-14 0 0 0 1e-14'
[]
[relperm0]
type = PorousFlowRelativePermeabilityCorey
n = 2
phase = 0
[]
[relperm1]
type = PorousFlowRelativePermeabilityCorey
n = 3
phase = 1
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.2
[]
[internal_energy]
type = PorousFlowMatrixInternalEnergy
density = 2650
specific_heat_capacity = 1000
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
end_time = 1e3
nl_abs_tol = 1e-12
[TimeStepper]
type = IterationAdaptiveDT
dt = 10
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Postprocessors]
[density_water]
type = ElementAverageValue
variable = density_water
execute_on = 'initial timestep_end'
[]
[density_gas]
type = ElementAverageValue
variable = density_gas
execute_on = 'initial timestep_end'
[]
[viscosity_water]
type = ElementAverageValue
variable = viscosity_water
execute_on = 'initial timestep_end'
[]
[viscosity_gas]
type = ElementAverageValue
variable = viscosity_gas
execute_on = 'initial timestep_end'
[]
[enthalpy_water]
type = ElementAverageValue
variable = enthalpy_water
execute_on = 'initial timestep_end'
[]
[enthalpy_gas]
type = ElementAverageValue
variable = enthalpy_gas
execute_on = 'initial timestep_end'
[]
[sg]
type = ElementAverageValue
variable = saturation_gas
execute_on = 'initial timestep_end'
[]
[sw]
type = ElementAverageValue
variable = saturation_water
execute_on = 'initial timestep_end'
[]
[pwater]
type = ElementAverageValue
variable = pressure_water
execute_on = 'initial timestep_end'
[]
[pgas]
type = ElementAverageValue
variable = pressure_gas
execute_on = 'initial timestep_end'
[]
[temperature]
type = ElementAverageValue
variable = temperature
execute_on = 'initial timestep_end'
[]
[enthalpy]
type = ElementAverageValue
variable = h
execute_on = 'initial timestep_end'
[]
[pliq]
type = ElementAverageValue
variable = pliq
execute_on = 'initial timestep_end'
[]
[liquid_mass]
type = PorousFlowFluidMass
phase = 0
execute_on = 'initial timestep_end'
[]
[vapor_mass]
type = PorousFlowFluidMass
phase = 1
execute_on = 'initial timestep_end'
[]
[liquid_heat]
type = PorousFlowHeatEnergy
phase = 0
execute_on = 'initial timestep_end'
[]
[vapor_heat]
type = PorousFlowHeatEnergy
phase = 1
execute_on = 'initial timestep_end'
[]
[e_water]
type = ElementAverageValue
variable = e_water
execute_on = 'initial timestep_end'
[]
[e_gas]
type = ElementAverageValue
variable = e_gas
execute_on = 'initial timestep_end'
[]
[]
[Outputs]
csv = true
perf_graph = false
[]
(modules/porous_flow/test/tests/energy_conservation/heat03_rz.i)
# The sample is a single unit element in RZ coordinates
# A constant velocity is applied to the outer boundary: disp_r = -0.01*t.
# There is no fluid flow or heat flow.
# Heat energy conservation is checked.
# Mass conservation is checked
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 1
xmin = 1
xmax = 2
ymin = -0.5
ymax = 0.5
coord_type = RZ
[]
[GlobalParams]
displacements = 'disp_r disp_z'
PorousFlowDictator = dictator
block = 0
biot_coefficient = 0.3
[]
[Variables]
[disp_r]
[]
[disp_z]
[]
[pp]
initial_condition = 0.1
[]
[temp]
initial_condition = 10
[]
[]
[BCs]
[plane_strain]
type = DirichletBC
variable = disp_z
value = 0
boundary = 'bottom top'
[]
[rmin_fixed]
type = DirichletBC
variable = disp_r
value = 0
boundary = left
[]
[contract]
type = FunctionDirichletBC
variable = disp_r
function = -0.01*t
boundary = right
[]
[]
[PorousFlowFullySaturated]
coupling_type = ThermoHydroMechanical
porepressure = pp
temperature = temp
fp = simple_fluid
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 0.5
density0 = 1
viscosity = 1
thermal_expansion = 0
cv = 1.3
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeElasticityTensor
C_ijkl = '1 1.5'
# bulk modulus is lambda + 2*mu/3 = 1 + 2*1.5/3 = 2
fill_method = symmetric_isotropic
[]
[strain]
type = ComputeAxisymmetricRZSmallStrain
[]
[stress]
type = ComputeLinearElasticStress
[]
[porosity]
type = PorousFlowPorosity
porosity_zero = 0.1
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 2.2
density = 0.5
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '0.5 0 0 0 0.5 0 0 0 0.5'
[]
[thermal_cond]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '1 0 0 0 1 0 0 0 1'
[]
[]
[Postprocessors]
[p0]
type = PointValue
outputs = 'console csv'
execute_on = 'initial timestep_end'
point = '1 0 0'
variable = pp
[]
[t0]
type = PointValue
outputs = 'console csv'
execute_on = 'initial timestep_end'
point = '1 0 0'
variable = temp
[]
[rdisp]
type = PointValue
outputs = 'csv console'
point = '2 0 0'
use_displaced_mesh = false
variable = disp_r
[]
[fluid_mass]
type = PorousFlowFluidMass
fluid_component = 0
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[total_heat]
type = PorousFlowHeatEnergy
phase = 0
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[rock_heat]
type = PorousFlowHeatEnergy
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[fluid_heat]
type = PorousFlowHeatEnergy
include_porous_skeleton = false
phase = 0
execute_on = 'initial timestep_end'
outputs = 'console csv'
[]
[]
[Preconditioning]
[andy]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = Newton
dt = 2
end_time = 10
[]
[Outputs]
execute_on = 'initial timestep_end'
[csv]
type = CSV
[]
[]
(modules/porous_flow/test/tests/dirackernels/hfrompps.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 3
ny = 3
[]
[GlobalParams]
PorousFlowDictator = dictator
[]
[Variables]
[pressure]
[]
[temperature]
scaling = 1E-6
[]
[]
[ICs]
[pressure_ic]
type = ConstantIC
variable = pressure
value = 1e6
[]
[temperature_ic]
type = ConstantIC
variable = temperature
value = 400
[]
[]
[Kernels]
[P_time_deriv]
type = PorousFlowMassTimeDerivative
fluid_component = 0
variable = pressure
[]
[P_flux]
type = PorousFlowAdvectiveFlux
fluid_component = 0
variable = pressure
gravity = '0 -9.8 0'
[]
[energy_dot]
type = PorousFlowEnergyTimeDerivative
variable = temperature
[]
[heat_conduction]
type = PorousFlowHeatConduction
variable = temperature
[]
[heat_advection]
type = PorousFlowHeatAdvection
variable = temperature
gravity = '0 -9.8 0'
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pressure temperature'
number_fluid_phases = 1
number_fluid_components = 1
[]
[pc]
type = PorousFlowCapillaryPressureConst
[]
[]
[Functions]
[mass_flux_in_fn]
type = PiecewiseConstant
direction = left
xy_data = '
0 0
100 0.1
300 0
600 0.1
1400 0
1500 0.2'
[]
[T_in_fn]
type = PiecewiseLinear
xy_data = '
0 400
600 450'
[]
[]
[FluidProperties]
[simple_fluid]
type = SimpleFluidProperties
bulk_modulus = 2e9
density0 = 1000
thermal_expansion = 0
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temperature
[]
[ppss]
type = PorousFlow1PhaseP
porepressure = pressure
capillary_pressure = pc
[]
[massfrac]
type = PorousFlowMassFraction
at_nodes = true
[]
[fluid_props]
type = PorousFlowSingleComponentFluid
phase = 0
fp = simple_fluid
[]
[relperm]
type = PorousFlowRelativePermeabilityCorey
n = 1
phase = 0
[]
[fp_mat]
type = FluidPropertiesMaterialPT
pressure = pressure
temperature = temperature
fp = simple_fluid
[]
[rock_heat]
type = PorousFlowMatrixInternalEnergy
specific_heat_capacity = 830.0
density = 2750
[]
[thermal_conductivity]
type = PorousFlowThermalConductivityIdeal
dry_thermal_conductivity = '2.5 0 0 0 2.5 0 0 0 2.5'
[]
[permeability]
type = PorousFlowPermeabilityConst
permeability = '1.0E-15 0 0 0 1.0E-15 0 0 0 1.0E-14'
[]
[porosity]
type = PorousFlowPorosityConst
porosity = 0.1
[]
[]
[DiracKernels]
[source]
type = PorousFlowPointSourceFromPostprocessor
variable = pressure
mass_flux = mass_flux_in
point = '0.5 0.5 0'
[]
[source_h]
type = PorousFlowPointEnthalpySourceFromPostprocessor
variable = temperature
mass_flux = mass_flux_in
point = '0.5 0.5 0'
T_in = T_in
pressure = pressure
fp = simple_fluid
[]
[]
[Preconditioning]
[preferred]
type = SMP
full = true
petsc_options_iname = '-pc_type'
petsc_options_value = ' lu '
[]
[]
[Postprocessors]
[total_mass]
type = PorousFlowFluidMass
execute_on = 'initial timestep_end'
[]
[total_heat]
type = PorousFlowHeatEnergy
[]
[mass_flux_in]
type = FunctionValuePostprocessor
function = mass_flux_in_fn
execute_on = 'initial timestep_end'
[]
[avg_temp]
type = ElementAverageValue
variable = temperature
execute_on = 'initial timestep_end'
[]
[T_in]
type = FunctionValuePostprocessor
function = T_in_fn
execute_on = 'initial timestep_end'
[]
[]
[Executioner]
type = Transient
solve_type = Newton
nl_abs_tol = 1e-14
dt = 100
end_time = 2000
[]
[Outputs]
csv = true
execute_on = 'initial timestep_end'
file_base = hfrompps
[]