Start | Previous | Next

Porous Flow Tutorial Page 03. Adding heat advection and conduction

This Page adds heat conduction and advection with the fluid. The differential equation governing the temperature evolution is (1)

This equation is nonlinear since there are products of and , as well as the nonlinear function . Most of the nomenclature was used in Page 01, and the additional symbols introduced are:

  • is time (units s)

  • is the density of the rock grains (units kg.m)

  • is the specific heat capacity of the rock grains (units J.kg.K)

  • is the temperature (units K)

  • is the specific heat capacity of the fluid (units J.kg.K)

  • is the thermal conductivity of the rock-fluid system (units J.s.m.K). It is a tensor.

Before attempting to write an input file, a rough estimate of the expected nonlinear residuals must be performed, as discussed in convergence criteria. The residual for the Eq. (1) is approximately

where the parameters , , , , and , have been used in the final expression. In Page 02 the choice Pa.m was made. Choosing K.m yields

Note that this is significantly greater than the for the fluid equation. In the main, MOOSE can handle these types of discrepancies, but it is good practise to scale the variables so that their residuals are of similar magnitude. Therefore, a scaling factor of is applied to the temperature variable.

To model this thermo-hydro system, the PorousFlowBasicTHM action needs to be enhanced to read:

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [porepressure]
  []
  [temperature]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 293
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-8
  []
[]

[PorousFlowBasicTHM<<<{"href": "../../syntax/PorousFlowBasicTHM/index.html"}>>>]
  porepressure<<<{"description": "The name of the porepressure variable"}>>> = porepressure
  temperature<<<{"description": "For isothermal simulations, this is the temperature at which fluid properties (and stress-free strains) are evaluated at.  Otherwise, this is the name of the temperature variable.  Units = Kelvin"}>>> = temperature
  coupling_type<<<{"description": "The type of simulation.  For simulations involving Mechanical deformations, you will need to supply the correct Biot coefficient.  For simulations involving Thermal flows, you will need an associated ConstantThermalExpansionCoefficient Material"}>>> = ThermoHydro
  gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
  fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = the_simple_fluid
[]
(modules/porous_flow/examples/tutorial/03.i)

and some extra properties need to be added to the SimpleFluidProperties:

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

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

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [porepressure]
  []
  [temperature]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 293
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1E-8
  []
[]

[PorousFlowBasicTHM<<<{"href": "../../syntax/PorousFlowBasicTHM/index.html"}>>>]
  porepressure<<<{"description": "The name of the porepressure variable"}>>> = porepressure
  temperature<<<{"description": "For isothermal simulations, this is the temperature at which fluid properties (and stress-free strains) are evaluated at.  Otherwise, this is the name of the temperature variable.  Units = Kelvin"}>>> = temperature
  coupling_type<<<{"description": "The type of simulation.  For simulations involving Mechanical deformations, you will need to supply the correct Biot coefficient.  For simulations involving Thermal flows, you will need an associated ConstantThermalExpansionCoefficient Material"}>>> = ThermoHydro
  gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
  fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = the_simple_fluid
[]

[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
  [constant_injection_porepressure]
    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"}>>> = porepressure
    value<<<{"description": "Value of the BC"}>>> = 1E6
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
  []
  [constant_injection_temperature]
    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"}>>> = temperature
    value<<<{"description": "Value of the BC"}>>> = 313
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = injection_area
  []
[]

[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
  [the_simple_fluid]
    type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../source/fluidproperties/SimpleFluidProperties.html"}>>>
    bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 2E9
    viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 1.0E-3
    density0<<<{"description": "Density at zero pressure and zero temperature"}>>> = 1000.0
    thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0.0002
    cp<<<{"description": "Constant specific heat capacity at constant pressure (J/kg/K)"}>>> = 4194
    cv<<<{"description": "Constant specific heat capacity at constant volume (J/kg/K)"}>>> = 4186
    porepressure_coefficient<<<{"description": "The enthalpy is internal_energy + P / density * porepressure_coefficient.  Physically this should be 1.0, but analytic solutions are simplified when it is zero"}>>> = 0
  []
[]
(modules/porous_flow/examples/tutorial/03.i)

The boundary conditions used are the same as in Page 01 in addition to specifying a constant injection temperature of 313K:

  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 313
    boundary = injection_area
(modules/porous_flow/examples/tutorial/03.i)

Finally, some temperature-related Materials need to be defined

  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    biot_coefficient = 0.8
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caps aquifer'
  []
[]
(modules/porous_flow/examples/tutorial/03.i)

An animation of the results is shown in Figure 1. Readers are encouraged pause and explore the effect of changing parameters such as the rock thermal conductivity.

Figure 1: Temperature evolution in the borehole-aquifer-caprock system.

Start | Previous | Next