Browse files

Merge pull request #12819 from permcody/birdsell

Clarify how PorousFlowPiecewiseLinearSink can specify BCs
  • Loading branch information...
lindsayad committed Feb 4, 2019
2 parents 1bc46a1 + f664240 commit 377deb68257fb0c6d70fe91cfc49bd3d81b5684a
@@ -38,6 +38,7 @@ This basic sink boundary condition is implemented in [`PorousFlowSink`](PorousFl
The basic sink may be multiplied by a MOOSE Function of the pressure
of a fluid phase *or* the temperature:
s = f(t, x) \times g(P^{\beta} - P_{\mathrm{e}}) \ \ \ \textrm{or}\ \ \ s = f(t, x)
\times g(T - T_{\mathrm{e}}) \ .
@@ -56,7 +57,7 @@ In addition, the sink may be multiplied by any or all of the following

- Fluid relative permeability
- Fluid mobility ($k_{ij}n_{i}n_{j}k_{r} \rho / \mu$, where $n$ is the normal vector to the boundary)
- Fluid mobility ($k_{ij}n_{i}n_{j} \rho / \mu$, where $n$ is the normal vector to the boundary)
- Fluid mass fraction
- Fluid internal energy
- Thermal conductivity
@@ -75,6 +76,8 @@ Dirichlet conditions is inappropriate.
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.
@@ -94,50 +97,65 @@ environment. The flux of fluid from the model to this environment is
f = \frac{\rho k_{nn}
\frac{\rho k_{nn}
k_{\mathrm{r}}}{\mu}\frac{P-P_{\mathrm{e}}}{L} \ .
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 $k_{nn}$ is the permeability of the
region between the model and the imaginary environment,
region between the model and the imaginary environment
(which can also be thought of as the permeability of the adjacent element),
$k_{\mathrm{r}}$ is the relative permeability in this region, $\rho$
and $\mu$ are the fluid density and viscosity. The environment
porepressure is $P_{\mathrm{e}}$ and the boundary porepressure is

If $L\sim 0$ this effectively fixes the porepressure to $P\sim
P_{\mathrm{e}}$, since $f$ is very large otherwise (physically:
P_{\mathrm{e}}$, since the flux is very large otherwise (physically:
fluid is rapidly removed or added to the system by the environment to
ensure $P=P_{\mathrm{e}}$). If $L\sim\infty$ the boundary flux is
almost zero and does nothing.

[eq:fix_pp_bc] may be implemented in a number of ways, the 2 most
common being the following.
### Numerical Implementation

A [`PorousFlowPiecewiseLinearSink`](
may be constructed that models
f = C (P -
P_{\mathrm{e}}) \ ,
that is, with `pt_vals = '-1E9 1E9'`, `multipliers = '-1E9 1E9'`,
`PT_shift = Pe` and `flux_function = C`. Here the `1E9` is just an
example: the point is that it should be much greater than any expected
porepressure) and $P_{\mathrm{e}}$ provided by an AuxVariable (or set
to a constant value). The numerical value of the conductance, $C$, is
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:fix_pp_bc] may be implemented in a number of ways, and one of the most
common involves a [`PorousFlowPiecewiseLinearSink`](
that follows the format of [eq:s_f_g] using $f(x,t)=C$ and $g(P-P_{\mathrm{e}})$ 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 $g$ is shown in the figure below. It accepts $P-P_{\mathrm{e}}$ as an input and returns
a value that ends up multiplying $C$ to give a flux (as in [eq:s_f_g]). $C$ can be thought of as the conductance and is specified with `flux_function = C`. Its numerical value is discussed below. $P_{\mathrm{e}}$ can be specified using an AuxVariable or set to a constant value using `PT_shift = Pe`.

!media piecewiselinear_g_function.png style=width:50%;margin-left:10px caption=Depiction of $g(P-P_{\mathrm{e}})$ for PorousFlowPiecewiseLinearSink. The function accepts $P-P_{\mathrm{e}}$ 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 $C$ to give the flux out of the domain. id=PiecewiseLinear_g_Function

To set a Dirichlet boundary condition $P = P_{\mathrm{e}}$, either $C$ or the slope of $g$ should be very large.
It is usually convenient to make the slope of $g$ equal to 1 by setting `pt_vals = '-1E9 1E9'` and `multipliers = '-1E9 1E9'`, and then $C$ 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 $O(10^6)$).
This ensures good convergence, and if in doubt set the `pt_vals` at a higher value than you expect.
If $P - P_{\mathrm{e}}$ falls outside of the range defined in `pt_vals`, then the slope of $g = 0$. 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, $C$, is
$\rho k_{nn}k_{\mathrm{r}}/\mu/L$, must be set at an appropriate value
for the model.

Alternately, a
may be constructed with the same `pt_vals`, `multipliers` and `PT_shift`, but with `use_mobility = true` and `use_relperm =
true`, and then $C = 1/L$. 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 occuring.
numerical situations occurring.

Also note, if $C \times g$ 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 $C$? In the example below, $C = 10^{-5}$, $\rho \sim 10^3$ kg/m$^3$, $k = 10^{-15}$ m$^2$, $k_r = 1$, and $\mu \sim 10^{-3}$ Pa-s. Therefore $L \sim 10^{-4}$ m. This value of $L$ is small enough to ensure that the Dirichlet boundary condition is satisfied. If $C$ is increased to $10^{-2}$, $L \sim 10^{-7}$ m, and the simulation has difficulty converging. If $C = 10^{-11}$, $L\sim10^2$ 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 $C$ is simply $1/L$ if `use_mobility = true` and `use_relperm = true`.

!listing modules/porous_flow/test/tests/sinks/PorousFlowPiecewiseLinearSink_BC_eg1.i start=[BCs] end=[Postprocessors]

## Injection of fluid at a fixed temperature

@@ -148,7 +166,7 @@ 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 permeaility,
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 $P>P_{\mathrm{e}}$ using a
@@ -0,0 +1,144 @@
## This is an example input file showing how to set a Type I (Dirichlet) BC with PorousFlowPiecewiseLinearSink
## Problem setup:
## - The boundaries are set to P(x = 0) = 2e6 Pa, P(x = 1) = 1e6 and run to steady state.
## - The 2d domain is 1 m x 1 m
## - The permeability is set to 1E-15 m2, fluid viscosity = 1E-3 Pa-s
## - The steady state flux is calculated q = -k/mu*grad(P) = 1e-6 m/s
## Problem verification (in csv output):
## - The flux in and out of the domain are 1e-6 m/s (matching steady state solution)
## - The pressure at the left and right boundaries are set to 2e6 and 1e6 Pa, respectively

type = GeneratedMesh
dim = 2
nx = 5
xmin = 0
xmax = 1
ny = 2
ymin = 0
ymax = 1

PorousFlowDictator = dictator

initial_condition = 1.5e6 # initial pressure in domain

porepressure = porepressure
coupling_type = Hydro
gravity = '0 0 0'
fp = the_simple_fluid


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

type = NodalSum
boundary = 'left'
variable = fluxes_out
execute_on = 'timestep_end'
type = NodalSum
boundary = 'right'
variable = fluxes_in
execute_on = 'timestep_end'
type = SideAverageValue
boundary = 'left'
variable = porepressure
execute_on = 'timestep_end'
type = SideAverageValue
boundary = 'right'
variable = porepressure
execute_on = 'timestep_end'

type = SimpleFluidProperties
thermal_expansion = 0
viscosity = 1.0E-3
density0 = 1000.0

type = PorousFlowPorosity
porosity_zero = 0.1
type = PorousFlowConstantBiotModulus
biot_coefficient = 0.8
solid_bulk_compliance = 2E-7
fluid_bulk_modulus = 1E7
type = PorousFlowPermeabilityConst
permeability = '1E-15 0 0 0 1E-15 0 0 0 1E-15'
#### The following Material give porepressure at nodes, which is required for PorousFlowPiecewiseLienarSink
type = PorousFlow1PhaseFullySaturated
at_nodes = true
porepressure = porepressure

type = Transient
solve_type = Newton
end_time = 1E6
dt = 1E5
nl_abs_tol = 1E-10

csv = true
@@ -0,0 +1,12 @@
Oops, something went wrong.

0 comments on commit 377deb6

Please sign in to comment.