1D radial heat and mass transport

Description

An analytical solution to the problem of 1D radial coupled heat and mass transport was initially developed by Avdonin (1964), and later by Ross et al. (1982) (in a similar fashion as for the 1D Cartesian model).

The problem consists of a 1D radial model where cold water is injected into a warm semi-infinite reservoir at a constant rate. The top and bottom surfaces of the reservoir are bounded by caprock which is neglected in the modelling to simplify the problem. Instead, these boundaries are treated as no-flow and adiabatic boundary conditions.

For the simple case of a 1D radial model bounded on the upper and lower surfaces by no-flow and adiabatic boundaries, a simplified solution for the temperature profile $var element = document.getElementById("moose-equation-c033693f-2d47-4142-b176-54a6042d9cac");katex.render("T(r, t)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ can be obtained (Updegraff, 1989)

$(1)var element = document.getElementById("moose-equation-8445117f-b598-439a-9b74-13cc8b734ece");katex.render("\\frac{T(r, t) - T(r, 0)}{T(0, t) - T(r, 0)} = \\frac{\\Gamma(\\nu, \\omega^2/(4 \\tau))}{\\Gamma(\\nu)}, ", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where $var element = document.getElementById("moose-equation-c6e7db08-5ec5-4ce5-bf45-f9d9b6519153");katex.render("T(r, 0)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the initial temperature in the reservoir, $var element = document.getElementById("moose-equation-020d17d7-d2a3-4936-9bc5-bb9208b43709");katex.render("T(0, t)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the temperature of the injected water, $var element = document.getElementById("moose-equation-65fef12f-323e-4f2f-8f69-980e0b9d8312");katex.render("\\Gamma(x)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the gamma function, $var element = document.getElementById("moose-equation-55e28115-3166-4b69-977e-58c5506b2aa0");katex.render("\\Gamma(x, a)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the lower incomplete gamma function,

$var element = document.getElementById("moose-equation-de892a60-d14e-419f-bd0a-971d5c5026d6");katex.render("\\nu = \\frac{Q\\rho_w cp_w}{4 \\pi h \\lambda_m},", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

$var element = document.getElementById("moose-equation-fee76f9c-0c96-4dc0-92fb-03f0258a6eff");katex.render("\\omega = \\frac{2 r}{h},", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

and

$var element = document.getElementById("moose-equation-36489531-19bf-48fc-bf57-ba74c7695216");katex.render("\\tau = \\frac{4 \\lambda_m t}{\\rho_m cp_m h^2},", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where $var element = document.getElementById("moose-equation-0a213b0e-a25b-41c1-94cd-08bfb98db000");katex.render("\\rho_w", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the density of water, $var element = document.getElementById("moose-equation-3f256ad3-7cc2-482e-a297-52716024c3b4");katex.render("cp_w", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the specific heat capacity of water, $var element = document.getElementById("moose-equation-fec72bfe-510e-4892-be91-a413e2a6e078");katex.render("\\rho_m", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the density of the fully saturated medium ( $var element = document.getElementById("moose-equation-7e2b725a-869f-4090-b531-76850a887747");katex.render("\\rho_m = \\phi \\rho_w + (1 - \\phi \\rho_r)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ where $var element = document.getElementById("moose-equation-a239fa85-6b61-452b-a9c5-75f087454af8");katex.render("\\phi", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is porosity and $var element = document.getElementById("moose-equation-88352625-c7ce-45b9-aa0a-c5b302cb8a1f");katex.render("\\rho_r", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the density of the dry rock), $var element = document.getElementById("moose-equation-3db0ffcd-c2c8-4386-a67d-db8324066d21");katex.render("cp_m", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the specific heat capacity of the fully saturated porous medium, $var element = document.getElementById("moose-equation-1d38b1db-b2f3-4ffa-bbe0-56402ddaacab");katex.render("\\lambda_m", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the thermal conductivity of the fully saturated reservoir, $var element = document.getElementById("moose-equation-9a48a9f5-7e6b-4281-8474-8095900c10f7");katex.render("Q", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the volumetric flow rate, and $var element = document.getElementById("moose-equation-cabcaa2d-d328-4e5b-8c9d-2e9ac0e9574f");katex.render("h", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the height of the reservoir.

Model

This problem was considered in a code comparison by Updegraff (1989), so we use identical parameters in this verification problem, see Table 1.

Table 1: Model properties

Property	Value
Length	1,000 m
Pressure	5 MPa
Temperature	170 $var element = document.getElementById("moose-equation-592d5ba1-908c-49b3-a1b4-56720729e094");katex.render("^{\\circ}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ C
Permeability	$var element = document.getElementById("moose-equation-eeb52e65-4b18-4950-9cf1-663abc5f9650");katex.render("1.8 \\times 10^{-11}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ m $var element = document.getElementById("moose-equation-21574483-5e82-4d91-b579-c35e4b4a09fb");katex.render("^2", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$
Porosity	0.2
Saturate density	2,500 kg m $var element = document.getElementById("moose-equation-1a0874a9-4af8-4b2d-9934-26a0dea9674d");katex.render("^{-3}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$
Saturated thermal conductivity	25 W m $var element = document.getElementById("moose-equation-305e4574-2d81-4914-afd0-fd52e57715d6");katex.render("^{-1}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ K
Saturated specific heat capacity	1,000 J kg $var element = document.getElementById("moose-equation-55cad6c3-4bd9-4361-bb27-b911fc7275d5");katex.render("^{-1}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ K
Mass flux rate	0.1 kg s $var element = document.getElementById("moose-equation-df8a808d-67cd-48a4-ac7e-efaf9e5278e8");katex.render("{-1}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

Input file

The input file used to run this problem is

# Cold water injection into 1D radial hot reservoir (Avdonin, 1964)
#
# To generate results presented in documentation for this problem,
# set xmax = 1000 and nx = 200 in the Mesh block, and dtmax = 1e4
# and end_time = 1e6 in the Executioner block.

[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 1
  nx = 50
  xmin = 0.1
  xmax = 5
  bias_x = 1.05
  rz_coord_axis = Y
  coord_type = RZ
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [temperature]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
[]

[AuxKernels<<<{"href": "../../../../syntax/AuxKernels/index.html"}>>>]
  [temperature]
    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"}>>> = temperature
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = temperature
    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."}>>> = 'initial timestep_end'
  []
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pliquid]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 5e6
  []
  [h]
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1e-6
  []
[]

[ICs<<<{"href": "../../../../syntax/ICs/index.html"}>>>]
  [hic]
    type = PorousFlowFluidPropertyIC<<<{"description": "An initial condition to calculate one fluid property (such as enthalpy) from pressure and temperature", "href": "../../../../source/ics/PorousFlowFluidPropertyIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = h
    porepressure<<<{"description": "Fluid porepressure"}>>> = pliquid
    property<<<{"description": "The fluid property that this initial condition is to calculate"}>>> = enthalpy
    temperature<<<{"description": "Fluid temperature"}>>> = 170
    temperature_unit<<<{"description": "The unit of the temperature variable"}>>> = Celsius
    fp<<<{"description": "The name of the user object for the fluid"}>>> = water
  []
[]

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

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [source]
    type = PorousFlowSink<<<{"description": "Applies a flux sink to a boundary.", "href": "../../../../source/bcs/PorousFlowSink.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pliquid
    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)"}>>> = injection_rate
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
  []
  [pright]
    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"}>>> = pliquid
    value<<<{"description": "Value of the BC"}>>> = 5e6
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
  []
  [hleft]
    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"}>>> = h
    value<<<{"description": "Value of the BC"}>>> = 678.52e3
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
  []
  [hright]
    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"}>>> = h
    value<<<{"description": "Value of the BC"}>>> = 721.4e3
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
  []
[]

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [mass]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pliquid
  []
  [massflux]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pliquid
  []
  [heat]
    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"}>>> = h
  []
  [heatflux]
    type = PorousFlowHeatAdvection<<<{"description": "Fully-upwinded heat flux, advected by the fluid", "href": "../../../../source/kernels/PorousFlowHeatAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = h
  []
  [heatcond]
    type = PorousFlowHeatConduction<<<{"description": "Heat conduction in the Porous Flow module", "href": "../../../../source/kernels/PorousFlowHeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = h
  []
[]

[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."}>>> = 'pliquid h'
    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"}>>> = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG<<<{"description": "van Genuchten capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureVG.html"}>>>
    pc_max<<<{"description": "Maximum capillary pressure (Pa). Must be > 0. Default is 1e9"}>>> = 1e6
    sat_lr<<<{"description": "Liquid residual saturation.  Must be between 0 and 1. Default is 0"}>>> = 0.1
    m<<<{"description": "van Genuchten exponent m. Must be between 0 and 1, and optimally should be set to >0.5"}>>> = 0.5
    alpha<<<{"description": "van Genuchten parameter alpha. Must be positive"}>>> = 1e-5
  []
  [fs]
    type = PorousFlowWaterVapor<<<{"description": "Fluid state class for water and vapor", "href": "../../../../source/userobjects/PorousFlowWaterVapor.html"}>>>
    water_fp<<<{"description": "The name of the user object for water"}>>> = water
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [water]
    type = Water97FluidProperties<<<{"description": "Fluid properties for water and steam (H2O) using IAPWS-IF97", "href": "../../../../source/fluidproperties/Water97FluidProperties.html"}>>>
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [watervapor]
    type = PorousFlowFluidStateSingleComponent<<<{"description": "Class for single component multiphase fluid state calculations using pressure and enthalpy", "href": "../../../../source/materials/PorousFlowFluidStateSingleComponent.html"}>>>
    porepressure<<<{"description": "Variable that is the porepressure of the liquid phase"}>>> = pliquid
    enthalpy<<<{"description": "Enthalpy of the fluid"}>>> = h
    temperature_unit<<<{"description": "The unit of the temperature variable"}>>> = Celsius
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
    fluid_state<<<{"description": "Name of the FluidState UserObject"}>>> = fs
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1.8e-11 0 0 0 1.8e-11 0 0 0 1.8e-11'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
    s_res<<<{"description": "The residual saturation of the phase j. Must be between 0 and 1"}>>> = 0.1
    sum_s_res<<<{"description": "Sum of residual saturations over all phases.  Must be between 0 and 1"}>>> = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 1
    sum_s_res<<<{"description": "Sum of residual saturations over all phases.  Must be between 0 and 1"}>>> = 0.1
  []
  [internal_energy]
    type = PorousFlowMatrixInternalEnergy<<<{"description": "This Material calculates the internal energy of solid rock grains, which is specific_heat_capacity * density * temperature.  Kernels multiply this by (1 - porosity) to find the energy density of the porous rock in a rock-fluid system", "href": "../../../../source/materials/PorousFlowMatrixInternalEnergy.html"}>>>
    density<<<{"description": "Density of the rock grains"}>>> = 2900
    specific_heat_capacity<<<{"description": "Specific heat capacity of the rock grains (J/kg/K)."}>>> = 740
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal<<<{"description": "This Material calculates rock-fluid combined thermal conductivity by using a weighted sum.  Thermal conductivity = dry_thermal_conductivity + S^exponent * (wet_thermal_conductivity - dry_thermal_conductivity), where S is the aqueous saturation", "href": "../../../../source/materials/PorousFlowThermalConductivityIdeal.html"}>>>
    dry_thermal_conductivity<<<{"description": "The thermal conductivity of the rock matrix when the aqueous saturation is zero"}>>> = '20 0 0  0 20 0  0 0 20'
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [smp]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = NEWTON
  end_time = 1e3
  nl_abs_tol = 1e-8
  [TimeStepper<<<{"href": "../../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    dt = 100
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [injection_area]
    type = AreaPostprocessor<<<{"description": "Computes the \"area\" or dimension - 1 \"volume\" of a given boundary or boundaries in your mesh.", "href": "../../../../source/postprocessors/AreaPostprocessor.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    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."}>>> = initial
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [line]
    type = ElementValueSampler<<<{"description": "Samples values of variables on elements.", "href": "../../../../source/vectorpostprocessors/ElementValueSampler.html"}>>>
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = temperature
    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."}>>> = 'initial timestep_end'
  []
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = true
  [csv]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    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."}>>> = final
  []
[]

Note that the test file is a reduced version of this problem. To recreate these results, follow the instructions at the top of the input file.

Results

The results for the temperature profile after 13,000 seconds are shown in Figure 1. Excellent agreement between the analytical solution and the MOOSE results are observed.

Figure 1: Comparison between Avdonin (1964) result and MOOSE at $var element = document.getElementById("moose-equation-eb316516-8ea2-4b7b-be9b-6c433355c1fb");katex.render("t = 10^6", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ s (left); and $var element = document.getElementById("moose-equation-467211cc-02ec-4baf-9b0a-113968619d6d");katex.render("t = 10^9", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ s (right).

This model also admits a similarity solution $var element = document.getElementById("moose-equation-bc88557f-69c8-40be-b175-d326c987ae82");katex.render("\\zeta = r^2/t", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ (Moridis and Pruess, 1992). Again, excellent agreement between the analytical solution and the MOOSE results are observed, see Figure 2

Figure 2: Similarity solution for 1D radial problem.

References

NA Avdonin. Some formulas for calculating the temperature field of a stratum subject to thermal injection. Neft'i Gaz, 3:37–41, 1964.

@article{avdonin1964,
    author = "Avdonin, NA",
    title = "Some formulas for calculating the temperature field of a stratum subject to thermal injection",
    journal = "Neft'i Gaz",
    volume = "3",
    pages = "37--41",
    year = "1964"
}

GJ Moridis and K Pruess. Tough simulations of updegraff's set of fluid and heat flow problems. LBL-32611, Lawrence Berkeley Lab, 1992.

@article{moridis1992,
    author = "Moridis, GJ and Pruess, K",
    title = "TOUGH simulations of Updegraff's set of fluid and heat flow problems",
    journal = "LBL-32611, Lawrence Berkeley Lab",
    year = "1992"
}

B Ross, JW Mercer, SD Thomas, and BH Lester. Benchmark problems for repository siting models. Technical Report, GeoTrans, 1982.

@techreport{ross1982,
    author = "Ross, B and Mercer, JW and Thomas, SD and Lester, BH",
    title = "Benchmark problems for repository siting models",
    year = "1982",
    institution = "GeoTrans"
}

C David Updegraff. Comparison of strongly heat-driven flow codes for unsaturated media. Technical Report, Nuclear Regulatory Commission, 1989.

@techreport{updegraff1989,
    author = "Updegraff, C David",
    title = "Comparison of strongly heat-driven flow codes for unsaturated media",
    year = "1989",
    institution = "Nuclear Regulatory Commission"
}

(modules/porous_flow/test/tests/fluidstate/coldwater_injection_radial.i)

# Cold water injection into 1D radial hot reservoir (Avdonin, 1964)
#
# To generate results presented in documentation for this problem,
# set xmax = 1000 and nx = 200 in the Mesh block, and dtmax = 1e4
# and end_time = 1e6 in the Executioner block.

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 50
  xmin = 0.1
  xmax = 5
  bias_x = 1.05
  rz_coord_axis = Y
  coord_type = RZ
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [temperature]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [temperature]
    type = PorousFlowPropertyAux
    variable = temperature
    property = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Variables]
  [pliquid]
    initial_condition = 5e6
  []
  [h]
    scaling = 1e-6
  []
[]

[ICs]
  [hic]
    type = PorousFlowFluidPropertyIC
    variable = h
    porepressure = pliquid
    property = enthalpy
    temperature = 170
    temperature_unit = Celsius
    fp = water
  []
[]

[Functions]
  [injection_rate]
    type = ParsedFunction
    symbol_values = injection_area
    symbol_names = area
    expression = '-0.1/area'
  []
[]

[BCs]
  [source]
    type = PorousFlowSink
    variable = pliquid
    flux_function = injection_rate
    boundary = left
  []
  [pright]
    type = DirichletBC
    variable = pliquid
    value = 5e6
    boundary = right
  []
  [hleft]
    type = DirichletBC
    variable = h
    value = 678.52e3
    boundary = left
  []
  [hright]
    type = DirichletBC
    variable = h
    value = 721.4e3
    boundary = right
  []
[]

[Kernels]
  [mass]
    type = PorousFlowMassTimeDerivative
    variable = pliquid
  []
  [massflux]
    type = PorousFlowAdvectiveFlux
    variable = pliquid
  []
  [heat]
    type = PorousFlowEnergyTimeDerivative
    variable = h
  []
  [heatflux]
    type = PorousFlowHeatAdvection
    variable = h
  []
  [heatcond]
    type = PorousFlowHeatConduction
    variable = h
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pliquid h'
    number_fluid_phases = 2
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    pc_max = 1e6
    sat_lr = 0.1
    m = 0.5
    alpha = 1e-5
  []
  [fs]
    type = PorousFlowWaterVapor
    water_fp = water
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [watervapor]
    type = PorousFlowFluidStateSingleComponent
    porepressure = pliquid
    enthalpy = h
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.8e-11 0 0 0 1.8e-11 0 0 0 1.8e-11'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
    sum_s_res = 0.1
  []
  [internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2900
    specific_heat_capacity = 740
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '20 0 0  0 20 0  0 0 20'
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 1e3
  nl_abs_tol = 1e-8
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 100
  []
[]

[Postprocessors]
  [injection_area]
    type = AreaPostprocessor
    boundary = left
    execute_on = initial
  []
[]

[VectorPostprocessors]
  [line]
    type = ElementValueSampler
    sort_by = x
    variable = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  perf_graph = true
  [csv]
    type = CSV
    execute_on = final
  []
[]

Description
Model
Input file
Results
References