Boundary conditions

MOOSE's Dirichlet and Neumann boundary conditions enable simulation of simple scenarios. The Porous Flow module includes a very flexible boundary condition that allows many different scenarios to be modelled. The boundary condition is built by adding various options to a basic sink in the following way.

This page may be read in conjuction with the description of some of the tests of the PorousFlow sinks.

Basic sink formulation

The basic sink is (1) where is a MOOSE Function of time and position on the boundary.

If then the boundary condition will act as a sink, while if the boundary condition acts as a source. If applied to a fluid-component equation, the function has units kg.m.s. If applied to the heat equation, the function has units J.m.s. The units of are potentially modified if the extra building blocks enumerated below are used (but the units of the final result, , will always be kg.m.s or J.m.s).

If the fluid flow through the boundary is required, use the save_in command which will save the sink strength to an AuxVariable. This AuxVariable will be the flux (kg.s or J.s) from each node on the boundary, which is the product of and the area attributed to that node. If the total flux (kg.s or J.s) through the boundary is required, integrate the AuxVariable over the boundary using a NodalSum Postprocessor.

This basic sink boundary condition is implemented in PorousFlowSink.

More elaborate options

The basic sink may be multiplied by a MOOSE Function of the pressure of a fluid phase or the temperature: (2) Here the units of are kg.m.s (for fluids) or J.m.s (for heat). and are reference values, that may be AuxVariables to allow spatial and temporal variance. Some commonly use forms have been hard-coded into the Porous Flow module for ease of use in simulations:

In addition, the sink may be multiplied by any or all of the following quantities

  • Fluid relative permeability

  • Fluid mobility (, where is the normal vector to the boundary)

  • Fluid mass fraction

  • Fluid internal energy

  • Thermal conductivity

Fixing fluid porepressures using the PorousFlowSink

Frequently in PorousFlow simulations fixing fluid porepressures using Dirichlet conditions is inappropriate.

  • Physically a Dirichlet condition corresponds to an environment outside a model providing a limitless source (or sink) of fluid and that does not happen in reality.

  • The Variables used may not be porepressure, so a straightforward translation to "fixing porepressure" is not easy.

  • The Dirichlet condition can lead to extremely poor convergence as it forces areas near the boundary close to unphysical or unlikely regions in parameter space.

Physical Intuition Behind PorousFlowSink

It is advantageous to think of what the boundary condition is physically trying to represent before using Dirichlet conditions, and often a PorousFlowSink is more appropriate to use. For instance, at the top of a groundwater-gas model, gas can freely enter the model from the atmosphere, and escape the model to the atmosphere, which may be modelled using a PorousFlowSink. Water can exit the model if the water porepressure increases above atmospheric pressure, but the recharge if porepressure is lowered may be zero. Again this may be modelled using a PorousFlowSink.

When fixing porepressures (or temperatures) using a PorousFlowSink, what should the "conductance" be? Consider an imaginary environment that sits at distance from the model. The imaginary environment provides a limitless source or sink of fluid, but that fluid must travel through the distance between the model and the imaginary environment. The flux of fluid from the model to this environment is (approximately) (3) A similar equation holds for temperature (but in PorousFlow simulations the heat is sometimes supplied by the fluid, in which case the appropriate equation to use for the heat may be the above equation multiplied by the enthalpy). Here is the permeability of the region between the model and the imaginary environment (which can also be thought of as the permeability of the adjacent element), is the relative permeability in this region, and are the fluid density and viscosity. The environment porepressure is and the boundary porepressure is .

If this effectively fixes the porepressure to , since the flux is very large otherwise (physically: fluid is rapidly removed or added to the system by the environment to ensure ). If the boundary flux is almost zero and does nothing.

Numerical Implementation

The PorousFlowSink is implemented in a fairly general way that allows for flexibility in setting combinations of pressure and temperature boundary conditions. Due to its implementation, it is difficult to draw a perfect analogy to the physical flux. Nevertheless, Eq. (3) may be implemented in a number of ways, and one of the most common involves a PorousFlowPiecewiseLinearSink that follows the format of Eq. (2) using and as a piecewise linear function between ordered pairs of pt_vals (on the x-axis) and multiplier (on the y-axis). An example of the function is shown in the figure below. It accepts as an input and returns a value that ends up multiplying to give a flux (as in Eq. (2)). can be thought of as the conductance and is specified with flux_function = C. Its numerical value is discussed below. can be specified using an AuxVariable or set to a constant value using PT_shift = Pe.

Figure 1: Depiction of for PorousFlowPiecewiseLinearSink. The function accepts as an input (i.e. the difference between a specified environment pressure and the pressure on the boundary element) and returns a value that multiplies to give the flux out of the domain.

To set a Dirichlet boundary condition , either or the slope of should be very large. It is usually convenient to make the slope of equal to 1 by setting pt_vals = '-1E9 1E9' and multipliers = '-1E9 1E9', and then can be selected appropriately. The range for pt_vals should be sufficiently large that the simulation always occupies the region between the values specified (-1E9 1E9 is a typical choice because porepressures encountered in many simulations are ). This ensures good convergence, and if in doubt set the pt_vals at a higher value than you expect. If falls outside of the range defined in pt_vals, then the slope of . This can be useful to set a boundary condition that will only allow for outflow (e.g. by using pt_vals = '0 1E9', multipliers = '0 1E9').`

The numerical value of the conductance, , is , must be set at an appropriate value for the model. Alternately, a PorousFlowPiecewiseLinearSink may be constructed with the same pt_vals, multipliers and PT_shift, but with use_mobility = true and use_relperm = true, and then . This has three advantages: (1) the MOOSE input file is simpler; (2) MOOSE automatically chooses the correct mobility and relative permeability to use; (3) these parameters are appropriate to the model so it reduces the potential for difficult numerical situations occurring.

Also note, if is too large, the boundary residual will be much larger than residuals within the domain. This results in poor convergence.

So what value should be assigned to ? In the example below, , kg/m, m, , and Pa-s. Therefore m. This value of is small enough to ensure that the Dirichlet boundary condition is satisfied. If is increased to , m, and the simulation has difficulty converging. If , m, and the boundary acts like a source of fluid from a distant reservoir (i.e. it no longer acts like a Dirichlet boundary condition). The value of is simply if use_mobility = true and use_relperm = true.

[BCs]
  [./in_left]
    type = PorousFlowPiecewiseLinearSink
    variable = porepressure
    boundary = 'left'
    pt_vals = '-1e9 1e9' # x coordinates defining g
    multipliers = '-1e9 1e9' # y coordinates defining g
    PT_shift = 2.E6   # BC pressure
    flux_function = 1E-5 # Variable C
    fluid_phase = 0
    save_in = fluxes_out
  [../]
  [./out_right]
    type = PorousFlowPiecewiseLinearSink
    variable = porepressure
    boundary = 'right'
    pt_vals = '-1e9 1e9' # x coordinates defining g
    multipliers = '-1e9 1e9' # y coordinates defining g
    PT_shift = 1.E6   # BC pressure
    flux_function = 1E-6 # Variable C
    fluid_phase = 0
    save_in = fluxes_in
  [../]
[]
(modules/porous_flow/test/tests/sinks/PorousFlowPiecewiseLinearSink_BC_eg1.i)

Injection of fluid at a fixed temperature

A boundary condition corresponding to injection of fluid at a fixed temperature could involve: (1) using a Dirichlet condition for temperature; (2) using without any multiplicative factors. More complicated examples with heat and fluid injection and production are detailed in the test suite documentation.

Fluids that both enter and exit the boundary

Commonly, if fluid or heat is exiting the porous material, multiplication by relative permeability, mobility, mass fraction, enthalpy or thermal conductivity is necessary, while if fluid or heat is entering the porous material this multiplication is not necessary. Two sinks can be constructed in a MOOSE input file: one involving production which is only active for using a Piecewise-linear sink (here is the environmental pressure); and one involving injection, which is only active for using a piecewise-linear sink multiplied by the appropriate factors.

A 2-phase example

To illustrate that the "sink" part of the PorousFlowSink boundary condition works as planned, consider 2-phase flow along a line. The simulation is initially saturated with water, and CO is injected at the left-side. The CO front moves towards the right, but what boundary conditions should be placed on the right so that they don't interfere with the flow?

The answer is that two PorousFlowPiecewiseLinearSink are needed. One pulls out water and the other pulls out CO. They extract the fluids in proportion to their mass fractions, mobilities and relative permeabilities, in the way described above. Let us explore the input file a little.

In this example there are two fluid components:

  • fluid component 0 is water

  • fluid component 1 is CO

and two phases:

  • phase 0 is liquid

  • phase 1 is gas

The water component only exists in the liquid phase and the CO component only exists in the gaseous phase (the fluids are immiscible). The Variables are the water and gas porepressures:

[Variables]
  [./pwater]
    initial_condition = 20E6
  [../]
  [./pgas]
    initial_condition = 20.1E6
  [../]
[]
(modules/porous_flow/test/tests/sinks/injection_production_eg.i)

The pwater Variable is associated with the water component, while the pgas Variable is associated with the CO component:

[Kernels]
  [./mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pwater
  [../]
  [./flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pwater
  [../]
  [./mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = pgas
  [../]
  [./flux1]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    variable = pgas
  [../]
[]
(modules/porous_flow/test/tests/sinks/injection_production_eg.i)

A van Genuchten capillary pressure is used

  [./pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  [../]
(modules/porous_flow/test/tests/sinks/injection_production_eg.i)

The remainder of the input file is pretty standard, save for the important BCs block:

[BCs]
  [./co2_injection]
    type = PorousFlowSink
    boundary = left
    variable = pgas # pgas is associated with the CO2 mass balance (fluid_component = 1 in its Kernels)
    flux_function = -1E-2 # negative means a source, rather than a sink
  [../]

  [./right_water]
    type = PorousFlowPiecewiseLinearSink
    boundary = right

    # a sink of water, since the Kernels given to pwater are for fluid_component = 0 (the water)
    variable = pwater

    # this Sink is a function of liquid porepressure
    # Also, all the mass_fraction, mobility and relperm are referenced to the liquid phase now
    fluid_phase = 0

    # Sink strength = (Pwater - 20E6)
    pt_vals = '0 1E9'
    multipliers = '0 1E9'
    PT_shift = 20E6

    # multiply Sink strength computed above by mass fraction of water at the boundary
    mass_fraction_component = 0

    # also multiply Sink strength by mobility of the liquid
    use_mobility = true

    # also multiply Sink strength by the relperm of the liquid
    use_relperm = true

    # also multiplly Sink strength by 1/L, where L is the distance to the fixed-porepressure external environment
    flux_function = 10 # 1/L
  [../]

  [./right_co2]
    type = PorousFlowPiecewiseLinearSink
    boundary = right
    variable = pgas
    fluid_phase = 1
    pt_vals = '0 1E9'
    multipliers = '0 1E9'
    PT_shift = 20.1E6
    mass_fraction_component = 1
    use_mobility = true
    use_relperm = true
    flux_function = 10 # 1/L
  [../]
[]
(modules/porous_flow/test/tests/sinks/injection_production_eg.i)

Below are shown some outputs. Evidently the boundary condition satisfies the requirements.

Figure 2: Gas saturation at breakthrough (the jaggedness of the line is because gas saturation is an elemental variable).

Figure 3: Porepressures at breakthrough.

Figure 4: Gas saturation at the end of simulation.

Figure 5: Porepressures at the end of simulation.

Figure 6: Evolution of gas saturation (the jaggedness of the line is because gas saturation is an elemental variable).

Boundary terms for the advection-diffusion equation

The theoretical background concerning boundary fluxes is developed in this section. The presentation is general and not specific to Porous Flow. The variable may represent temperature, fluid mass, mass-densty of a species, number-density of a species, etc. The flux quantified below has units depending on the physical meaning of : if is temperature then the flux has units J.s; if is fluid mass then the flux has units kg.s, etc. For sake or argument, is referred to as "fluid" in the discussion below.

The advection equation

The advection equation for is (4) The velocity field advects , meaning that follows the streamlines of . The weak form of this equation reads (5) In a MOOSE input file, the second term is typically represented by a ConservativeAdvection Kernel, while the final term is the flux out through the surface and can be represented by an Outflow or Inflow BC.

The flux out through an area is (6) Here is the outward normal to . If is an arbitrary internal area, this equation is always correct, but if is on the boundary of the simulation, the flux out is dictated by the boundary condition used.

There are various common boundary conditions. For the sake of the arguments below, imagine being an outside boundary of the MOOSE numerical model.

  • If then fluid is moving from the numerical model through the boundary towards the outside of the model. If no BC is used then there is no term in the system of equations that is removing fluid from the boundary. That is, no fluid is allowed out of the model through , and fluid builds up at the boundary. (Even though , does not have to be zero on the boundary.)

  • If then fluid is moving from the numerical model through the boundary towards the outside of the model. An Outflow BC is exactly the final term in Eq. (5). This removes fluid through at exactly the rate specified by the advection equation, so this is a "free" boundary that is "invisible" to the simulation.

  • If then fluid is moving from the numerical model through the boundary towards the outside of the model. If no BC is used then no fluid enters the domain through (and will reduce).

  • If then fluid is moving from the numerical model through the boundary towards the outside of the model. An Inflow BC, which is exactly the last term in Eq. (5) with fixed, may be used to specify an injection rate through this boundary ( has units kg.s.m).

The diffusion equation

The diffusion equation for (representing temperature, fluid mass, mass-density of a species, etc) is (7) with weak form (8) The second is typically represented in MOOSE using the Diffusion Kernel, while the last term is the DiffusionFluxBC BC (assuming ). This last term is the flux out (kg.s, or J.s, or whatever are appropriate units associated with ) through an area .

The flux out through an arbitrary area is (9) Here is the outward normal to . If is an arbitrary internal area, this equation is always correct, but if is on the boundary of the simulation, the flux out is dictated by the boundary condition used.

Inclusion of a BC for in a MOOSE input file will specify something about the flux from the boundary. Some examples are:

  • A NeumannBC simply adds a constant value (kg.s.m) to the Residual and integrates it over the boundary. Adding this constant value corresponds to adding a constant flux of fluid and MOOSE will find the solution that has . The value is the flux into the domain (negative of the out-going flux mentioned above).

  • Using no BC at all on means that there is no boundary term in Eq. (8) so no fluid gets removed from this boundary: it is impermeable. MOOSE will find the solution that has .

  • Using a DirichletBC fixes on the boundary. The physical interpretation is that an external source or sink is providing or removing fluid at the correct rate so that remains fixed.

  • Using a DiffusionFluxBC will remove fluid at exactly the rate specified by the diffusion equation (assuming : otherwise an extension to the DiffusionFluxBC class needs to be encoded into MOOSE). This boundary condition is a "free" boundary that is "invisible" to the simulation.

The advection-diffusion equation

The advection-diffusion equation for (representing temperature, fluid mass, mass-density of a species, etc) is just a combination of the above: (10) with weak form (11)

The flux out (kg.s, or J.s, or whatever are appropriate units associated with ) through an area is (12) Here is the outward normal to .

Inclusion of a BC for in a MOOSE input file will specify something about the flux from the boundary. Some examples are.

  • Using no BCs means the final term Eq. (11) is missing from the system of equations that MOOSE solves. This means that no fluid is entering or exiting the domain from this boundary: the boundary is "impermeable". The fluid will "pile up" at the boundary, if is such that it is attempting to "blow" fluid out of the boundary (). Conversely, the fluid will deplete at the boundary, if is attempting to "blow" fluid from the boundary into the model domain ().

  • A NeumannBC adds a constant value (kg.s.m) to the Residual and integrates it over the boundary. Adding this constant value corresponds to adding a constant flux of fluid and MOOSE will find the solution that has . The value of is the flux into the domain (negative of the out-going flux mentioned above).

  • Using an OutflowBC together with a DiffusionFluxBC removes fluid through at exactly the rate specified by the advection-diffusion equation, so this is a "free" boundary that is "invisible" to the simulation.