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.
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
This basic sink boundary condition is implemented in
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:
A piecewise linear function (for simulating fluid or heat exchange with an external environment via a conductivity term, for instance), implemented in
A half-gaussian (for simulating evapotranspiration, for instance), implemented in
A half-cubic (for simulating evapotranspiration, for example), implemented in
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
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
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.
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.
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_shift, but with
use_mobility = true and
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 [../] 
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 [../] 
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 [../] 
A van Genuchten capillary pressure is used
[./pc] type = PorousFlowCapillaryPressureVG alpha = 1E-6 m = 0.6 [../]
The remainder of the input file is pretty standard, save for the important
[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 [../] 
Below are shown some outputs. Evidently the boundary condition satisfies the requirements.
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
Kernel, while the final term is the flux out through the surface and can be represented by an
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
OutflowBC 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
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:
NeumannBCsimply 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 .
DirichletBCfixes 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.
DiffusionFluxBCwill 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 ().
NeumannBCadds 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).
OutflowBCtogether with a
DiffusionFluxBCremoves fluid through at exactly the rate specified by the advection-diffusion equation, so this is a "free" boundary that is "invisible" to the simulation.