Porous Flow Tutorial Page 08. The PorousFlowSink and unsaturated flow
The PorousFlowSink has been mentioned in this tutorial already, but now is an appropriate time to study it more carefully. It is a very powerful BC which can handle virtually all physical situations. It also helps enormously in resolving some convergence problems. Therefore, it is well worth studying, so please feel free to take a break from this tutorial and read the info here.
Now that you've finished studying the PorousFlowSink let's apply it in a simple situation. The purely fluid-flow model of Page 01 is extended to include unsaturated flow. Firstly, the PorousFlowBasicTHM Action must be replaced by a PorousFlowUnsaturated Action:
[PorousFlowUnsaturated<<<{"href": "../../syntax/PorousFlowUnsaturated/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
relative_permeability_exponent<<<{"description": "Relative permeability exponent"}>>> = 3
relative_permeability_type<<<{"description": "Type of relative-permeability function. FLAC relperm = (1+m)S^m - mS^(1+m). Corey relperm = S^m. m is the exponent. Here S = (saturation - residual)/(1 - residual)"}>>> = Corey
residual_saturation<<<{"description": "Residual saturation to use in the relative permeability expression"}>>> = 0.1
van_genuchten_alpha<<<{"description": "Van Genuchten alpha parameter used to determine saturation from porepressure"}>>> = 1E-6
van_genuchten_m<<<{"description": "Van Genuchten m parameter used to determine saturation from porepressure"}>>> = 0.6
[]There are evidently some more input parameters:
The
relative_permeabilityparameters. These are documented in Relative permeabilityThe
van_genuchtencapillarity parameters. These are documented in Capillary pressure
The PorousFlowConstantBiotModulus is not needed here (it's only needed by PorousFlowBasicTHM). Because of the multiplication of the fluid mass-balance equation by density, an appropriate nl_abs_tol is in contrast to the calculated on Page 02:
[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-7
[]Finally we get to the PorousFlowSink. In the following, fluid is extracted through injection_area at a constant rate of kg.s.m. However, in the unsaturated region () this rate is modified by the relative permeability because use_relperm = true. This greatly improves convergence and in many cases is physically reasonable as there are limits to pumps and other similar machinery.
[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
[production]
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"}>>> = porepressure
fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy. If not supplied, then this BC can only be a function of temperature"}>>> = 0
flux_function<<<{"description": "The flux. The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE. The functional form is useful for spatially or temporally varying sinks. Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = 1E-2
use_relperm<<<{"description": "If true, then fluxes are multiplied by relative permeability. This can be used in conjunction with other use_*"}>>> = true
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
[]
[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
petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = ' asm lu NONZERO 2'
[]
[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-7
[]
[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
[]An animation of the results is shown in Figure 1.
Figure 1: Porepressure evolution in the production version of the borehole-aquifer-caprock system.
Try to run the simulation with use_relperm = false and you'll see very poor convergence. The total fluid extracted is kg (the product of rate, area and total time), while the model contains kg (the product of porosity, fluid density, and model volume), so the model is not running completely dry, but it requires extremely low porepressures to affect the constant rate of fluid removal ( gets so low that MOOSE starts to suffer from precision loss).
The simulation may be promoted to a full THMC simulation using the approach used in Page 03, Page 04 and Page 06. The arguments made about scaling the variables must be modified to take into account the fluid density appearing in the fluid equation (see Page 06) so the scaling will be approximately for the temperature and for the displacement variables. The porosity may depend on porepressure, temperature, volumetric strain and chemistry.
The simulation described so far uses full upwinding, which is the default in PorousFlow. TVD stabilization (see also numerical diffusion and a worked example of KT stabilization) may be used instead by simply changing the PorousFlowUnsaturated block to:
[PorousFlowUnsaturated<<<{"href": "../../syntax/PorousFlowUnsaturated/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
relative_permeability_exponent<<<{"description": "Relative permeability exponent"}>>> = 3
relative_permeability_type<<<{"description": "Type of relative-permeability function. FLAC relperm = (1+m)S^m - mS^(1+m). Corey relperm = S^m. m is the exponent. Here S = (saturation - residual)/(1 - residual)"}>>> = Corey
residual_saturation<<<{"description": "Residual saturation to use in the relative permeability expression"}>>> = 0.1
van_genuchten_alpha<<<{"description": "Van Genuchten alpha parameter used to determine saturation from porepressure"}>>> = 1E-6
van_genuchten_m<<<{"description": "Van Genuchten m parameter used to determine saturation from porepressure"}>>> = 0.6
stabilization<<<{"description": "Numerical stabilization used. 'Full' means full upwinding. 'KT' means FEM-TVD stabilization of Kuzmin-Turek"}>>> = KT
flux_limiter_type<<<{"description": "Type of flux limiter to use if stabilization=KT. 'None' means that no antidiffusion will be added in the Kuzmin-Turek scheme"}>>> = None
[]