Start | Previous | Next

Porous Flow Tutorial Page 01. A single fluid

This tutorial page describes how PorousFlow can be used to solve a very simple fluid-flow problem. The physics are described by the equation (1) This is the simplest type of equation that PorousFlow solves. It is often used in thermo-poro-elasticity, but PorousFlow is often employed to solve much more complex equations. See governing equations for further details, and PorousFlowFullySaturatedMassTimeDerivative and PorousFlowFullySaturatedDarcyBase for the derivation of Eq. (1) from the general governing equations. In this equation:

  • is the fluid porepressure (units of Pa)

  • an over-dot represents a time derivative

  • is the Biot Modulus (units of Pa)

  • is the Biot coefficient (dimensionless)

  • is the rate of volumetric strain of the solid rock (units s)

  • is the effective volumetric thermal expansion coefficient (units K)

  • is the temperature (units K)

  • represents a spatial derivative

  • is the permeability tensor (units m)

  • is the fluid viscosity (units Pa.s)

  • is the fluid density (units kg.m)

  • is the gravitational acceleration (units m.s)

The units suggested here are not mandatory: any consistent unit system may be used in MOOSE. For instance, in reservoir physics it is often convenient to express everything in MPa and years rather than Pa and seconds.

The Biot modulus is and the effective volumetric thermal expansion coefficient is In these equations

  • is the porosity (dimensionless)

  • is the bulk modulus of the fluid (units Pa)

  • is the bulk modulus of the drained porous skeleton (units Pa)

  • is the volumetric thermal expansion coefficient of the drained porous skeleton (units K)

  • is the volumetric thermal expansion coefficient of the fluid (units K)

The derivation of Eq. (1) from the full PorousFlow equations assumes that and are constant.

In this tutorial page we will be solving fluid flow only, so the and in Eq. Eq. (1) are ignored (set to zero).

All PorousFlow input files must contain a PorousFlowDictator, and almost all PorousFlow objects (Kernels, Materials, etc) require the PorousFlowDictator to be provided. This enables PorousFlow to make consistency checks on the number of fluid phases, components, chemical reactants, etc. Therefore, most input files specify its name in the Globals block:

[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
[]
(modules/porous_flow/examples/tutorial/01.i)

Most PorousFlow simulations require fluid properties to be supplied. In this instance, the SimpleFluidProperties are used, which assume a constant fluid bulk modulus and viscosity:

# Darcy flow
[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]
  []
[]

[PorousFlowBasicTHM<<<{"href": "../../syntax/PorousFlowBasicTHM/index.html"}>>>]
  porepressure<<<{"description": "The name of the porepressure variable"}>>> = porepressure
  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"}>>> = Hydro
  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
  []
[]

[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
  []
[]
(modules/porous_flow/examples/tutorial/01.i)

The DE of Eq. (1) is implemented in the following way

[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
  [porepressure]
  []
[]

[PorousFlowBasicTHM<<<{"href": "../../syntax/PorousFlowBasicTHM/index.html"}>>>]
  porepressure<<<{"description": "The name of the porepressure variable"}>>> = porepressure
  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"}>>> = Hydro
  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/01.i)

There is just one variable — the porepressure — and there is no coupling with heat or mechanics. Gravity is set to zero. The PorousFlowBasicTHM has other optional inputs that you are encouraged to explore, including setting the temperature to a non-default value, or to the value of an AuxVariable (your fluid properties may depend on temperature, even in an isothermal situation).

In this tutorial page, the only boundary condition is to fix the porepressure to 1MPa at the injection area (the other boundaries default to zero flux):

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

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

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [porosity]
    type = PorousFlowPorosity<<<{"description": "This Material calculates the porosity PorousFlow simulations", "href": "../../source/materials/PorousFlowPorosity.html"}>>>
    porosity_zero<<<{"description": "The porosity at zero volumetric strain and reference temperature and reference effective porepressure and reference chemistry.  This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)"}>>> = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus<<<{"description": "Computes the Biot Modulus, which is assumed to be constant for all time.  Sometimes 1 / BiotModulus is called storativity", "href": "../../source/materials/PorousFlowConstantBiotModulus.html"}>>>
    biot_coefficient<<<{"description": "Biot coefficient"}>>> = 0.8
    solid_bulk_compliance<<<{"description": "Reciprocal of the drained bulk modulus of the porous skeleton.  If strain = C * stress, then solid_bulk_compliance = de_ij de_kl C_ijkl.  If the grain bulk modulus is Kg then 1/Kg = (1 - biot_coefficient) * solid_bulk_compliance."}>>> = 2E-7
    fluid_bulk_modulus<<<{"description": "Fluid bulk modulus"}>>> = 1E7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquifer
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = caps
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []
[]

[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
  active<<<{"description": "If specified only the blocks named will be visited and made active"}>>> = basic
  [basic]
    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
  []
  [preferred_but_might_not_be_installed]
    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
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = ' lu       mumps'
  []
[]

[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-13
[]

[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
[]
(modules/porous_flow/examples/tutorial/01.i)

The porosity, Biot modulus and permeability are defined through the Materials block:

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [porosity]
    type = PorousFlowPorosity<<<{"description": "This Material calculates the porosity PorousFlow simulations", "href": "../../source/materials/PorousFlowPorosity.html"}>>>
    porosity_zero<<<{"description": "The porosity at zero volumetric strain and reference temperature and reference effective porepressure and reference chemistry.  This must be a real number or a constant monomial variable (not a linear lagrange or other type of variable)"}>>> = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus<<<{"description": "Computes the Biot Modulus, which is assumed to be constant for all time.  Sometimes 1 / BiotModulus is called storativity", "href": "../../source/materials/PorousFlowConstantBiotModulus.html"}>>>
    biot_coefficient<<<{"description": "Biot coefficient"}>>> = 0.8
    solid_bulk_compliance<<<{"description": "Reciprocal of the drained bulk modulus of the porous skeleton.  If strain = C * stress, then solid_bulk_compliance = de_ij de_kl C_ijkl.  If the grain bulk modulus is Kg then 1/Kg = (1 - biot_coefficient) * solid_bulk_compliance."}>>> = 2E-7
    fluid_bulk_modulus<<<{"description": "Fluid bulk modulus"}>>> = 1E7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquifer
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = caps
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []
[]
(modules/porous_flow/examples/tutorial/01.i)

The result is shown in Figure 1

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

Start | Previous | Next