Supported physical systems

JutulDarcy supports a number of different systems. These are JutulSystem instances that describe a particular type of physics for porous media flow. We describe these in roughly the order of complexity that they can model.

Summary

The general form of the flow systems we will discuss is a conservation law for $N$ components on residual form:

$$R=\frac{\partial }{\partial t}{M}_i+\mathrm{\nabla }\cdot {\overrightarrow{V}}_i-Q_i,\mathrm{\forall }i\in \{1,\dots ,N\}$$

Here, $M_i$ is the conserved quantity (usually masses) for component $i$ and ${\overrightarrow{V}}_i$ the velocity of the conserved quantity. $Q_i$ represents source terms that come from direct sources SourceTerm, boundary conditions (FlowBoundaryCondition) or from wells (setup_well, setup_vertical_well).

The following table gives an overview of the available features that are described in more detail below:

System	Number of phases	Number of components	$M$	$V$
SinglePhaseSystem	1	1	$\rho \varphi$	$\rho \overrightarrow{v}$
ImmiscibleSystem	Any	(Any)	$S_{\alpha }{\rho }_{\alpha }\varphi$	${\rho }_{\alpha }{\overrightarrow{v}}_{\alpha }$
StandardBlackOilSystem	2-3	(2-3)	${\rho }_o^s(b_g{S}_g+R_s{b}_o{S}_o)$	$b_g{\overrightarrow{v}}_g+R_s{b}_g{\overrightarrow{v}}_o$
MultiPhaseCompositionalSystemLV	2-3	Any	${\rho }_l{X}_i{S}_l+{\rho }_v{Y}_i{S}_v$	${\rho }_l{X}_i{\overrightarrow{v}}_l+{\rho }_v{Y}_i{\overrightarrow{v}}_v$

Phases

Phases are defined using specific types. Some constructors take a list of phases present in the model. Phases do not contain any data themselves and the distinction between different phases applies primarily for well controls.

JutulDarcy.LiquidPhase Type

LiquidPhase()

AbstractPhase subtype for liquid-like phases.

source

JutulDarcy.AqueousPhase Type

AqueousPhase()

AbstractPhase subtype for water-like phases.

source

JutulDarcy.VaporPhase Type

VaporPhase()

AbstractPhase subtype for vapor or gaseous phases.

source

Implementation details

In the above the discrete version of $M_i$ is implemented in the update function for JutulDarcy.TotalMasses that should by convention be named JutulDarcy.update_total_masses!. The discrete component fluxes are implemented by JutulDarcy.component_mass_fluxes!.

JutulDarcy.TotalMasses Type

TotalMasses()

Variable that defines total component masses in each cell of the domain.

source

JutulDarcy.update_total_masses! Function

update_total_masses!(totmass, tv, model, arg..., ix)

Update total masses for a given system. Number of input arguments varies based on physical system under consideration.

source

JutulDarcy.component_mass_fluxes! Function

component_mass_fluxes!(q, face, state, model, flux_type, kgrad, upw)

Implementation of component fluxes for a given system for a given face. Should return a StaticVector with one entry per component.

source

The source terms are implemented by Jutul.apply_forces_to_equation! for boundary conditions and sources, and Jutul.update_cross_term_in_entity! for wells. We use Julia's multiple dispatch to pair the right implementation with the right physics system.

Jutul.apply_forces_to_equation! Function

apply_forces_to_equation!(diag_part, storage, model, eq, eq_s, force, time)

Update an equation with the effect of a force. The default behavior for any force we do not know about is to assume that the force does not impact this particular equation.

source

Missing docstring.

Missing docstring for Jutul.update_cross_term_in_entity!. Check Documenter's build log for details.

Single-phase flow

JutulDarcy.SinglePhaseSystem Type

SinglePhaseSystem(phase = LiquidPhase(); reference_density = 1.0)

A single-phase system that only solves for pressure.

source

The simplest form of porous media flow is the single-phase system.

$$r(p)=\frac{\partial }{\partial t}(\rho \varphi )+\mathrm{\nabla }\cdot (\rho \overrightarrow{v})-\rho q$$

$\rho$ is the phase mass density and $\varphi$ the apparent porosity of the medium, i.e. the void space in the rock available to flow. Where the velocity $\overrightarrow{v}$ is given by Darcy's law that relates the the pressure gradient $\mathrm{\nabla }p$ and hydrostatic head to the velocity field:

$$\overrightarrow{v}=-\frac{\mathbf{K}}{\mu }(\mathrm{\nabla }p+\rho g\mathrm{\nabla }z)$$

Here, $\mathbf{K}$ is a positive-definite permeability tensor, $\mu$ the fluid viscosity, $g$ the magnitude of gravity oriented down and $z$ the depth.

Single-phase implementation

The SinglePhaseSystem is a dedicated single phase system. This is mathematically equivalent to an ImmiscibleSystem when set up with a single phase. For single phase flow, the fluid Pressure is the primary variable in each cell. The equation supports two types of compressibility: That of the fluid where density is a function $\rho (p)$ of pressure and that of the pores where the porosity $\varphi (p)$ changes with pressure.

Tip

JutulDarcy uses the notion of depth rather than coordinate when defining buoyancy forces. This is consistent with the convention in the literature on subsurface flow.

Multi-phase, immiscible flow

JutulDarcy.ImmiscibleSystem Type

ImmiscibleSystem(phases; reference_densities = ones(length(phases)))
ImmiscibleSystem((LiquidPhase(), VaporPhase()), reference_densities = (1000.0, 700.0))

Immiscible flow system: Each component exists only in a single phase, and the number of components equal the number of phases.

Set up an immiscible system for the given phases with optional reference densitites. This system is easy to specify with Pressure and Saturations as the default primary variables. Immiscible system assume that there is no mass transfer between phases and that a phase is uniform in composition.

source

The flow systems immediately become more interesting if we add more phases. We can extend the above single-phase system by introducing the phase saturation of phase with label $\alpha$ as $S_{\alpha }$. The phase saturation represents the volumetric fraction of the rock void space occupied by the phase. If we consider a pair of phases $\{n,w\}$ non-wetting and wetting we can write the system as

$$r_{\alpha }=\frac{\partial }{\partial t}(S_{\alpha }{\rho }_{\alpha }\varphi )+\mathrm{\nabla }\cdot ({\rho }_{\alpha }{\overrightarrow{v}}_{\alpha })-{\rho }_{\alpha }{q}_{\alpha }=0,\alpha \in \{n,w\}$$

This requires an additional closure such that the amount of saturation of all phases exactly fills the available fluid volume:

$$S_w+S_n=1,1\ge S_{\alpha }\ge 0\alpha \in \{n,w\}$$

This equation is local and linear in the saturations and can be eliminated to produce the classical two-equation system for two-phase flow,

$$\begin{gathered} r_n=\frac{\partial }{\partial t}((1-S_w){\rho }_n\varphi )+\mathrm{\nabla }\cdot ({\rho }_n{\overrightarrow{v}}_n)-{\rho }_n{q}_n=0, \\ {r}_w=\frac{\partial }{\partial t}(S_w{\rho }_w\varphi )+\mathrm{\nabla }\cdot ({\rho }_w{\overrightarrow{v}}_w)-{\rho }_w{q}_w=0. \end{gathered}$$

To complete this description we also need expressions for the phase fluxes. We use the standard multiphase extension of Darcy's law,

$${\overrightarrow{v}}_{\alpha }=-\mathbf{K}\frac{k_{r\alpha }}{{\mu }_{\alpha }}(\mathrm{\nabla }p\alpha +{\rho }_{\alpha }g\mathrm{\nabla }z)$$

Here, we have introduced the relative permeability of the phase $k_{r\alpha }$, an empirical relationship between the saturation and the flow rate. Relative permeability is a complex topic with many different relationships and functional forms, but we limit the discussion to monotone, non-negative functions of their respective saturations, for example a simple Brooks-Corey type of $k_{r\alpha }(S_{\alpha })=S_{\alpha }^2$. We have also introduced separate phase pressures $p_{\alpha }$ that account for capillary pressure, e.g. $p_w=p_n+p_c(S_w)$.

Primary variables

Immiscible implementation

The ImmiscibleSystem implements this system for any number of phases. The primary variables for this system is a single reference Pressure and phase Saturations. As we do not solve for the volume closure equation, there is one less degree of freedom associated with the saturations than there are number of phases.

Black-oil: Multi-phase, pseudo-compositional flow

JutulDarcy.StandardBlackOilSystem Type

StandardBlackOilSystem(; rs_max = nothing,
                         rv_max = nothing,
                         phases = (AqueousPhase(), LiquidPhase(), VaporPhase()),
                         reference_densities = [786.507, 1037.84, 0.969758])

Set up a standard black-oil system. Keyword arguments rs_max and rv_max can either be nothing or callable objects / functions for the maximum Rs and Rv as a function of pressure. phases can be specified together with reference_densities for each phase.

NOTE: For the black-oil model, the reference densities significantly impact many aspects of the PVT behavior. These should generally be set consistently with the other properties.

source

The black-oil equations is an extension of the immiscible description to handle limited miscibility between the phases. Originally developed for certain types of oil and gas simulation, these equations are useful when the number of components is low and tabulated values for dissolution and vaporization are available.

The assumptions of the black-oil model is that the "oil" and "gas" pseudo-components have uniform composition throughout the domain. JutulDarcy supports two- and three-phase black oil flow. The difference between two and three phases amounts to an additional immiscible aqueous phase that is identical to that of the previous section. For that reason, we focus on the miscible pseudo-components, oil:

$$r_o={\rho }_o^s(\frac{\partial }{\partial t}((b_o{S}_o+R_v{b}_g(1-S_o))\varphi )+\mathrm{\nabla }\cdot (b_o{\overrightarrow{v}}_o+R_v{b}_o{\overrightarrow{v}}_g)-q_o^s),$$

and gas,

$$r_g={\rho }_g^s(\frac{\partial }{\partial t}((b_g{S}_g+R_s{b}_o{S}_o)\varphi )+\mathrm{\nabla }\cdot (b_g{\overrightarrow{v}}_g+R_s{b}_g{\overrightarrow{v}}_o)-q_g^s)$$

The model uses the notion of surface (or reference densities) ${\rho }_o^s,{\rho }_g^s$ to define the densities of the component at specific pressure and temperature conditions where it is assumed that all "gas" has moved to the vapor phase and the defined "oil" is only found in the liquid phase. Keeping this definition in mind, the above equations can be divided by the surface densities to produce a surface volume balance equation where we have defined b_o and b_g as the dimensionless reciprocal formation volume factors that relate a volume at reservoir conditions to surface volumes and R_s for the dissolved volume of gas in the oil phase when brought to surface conditions. R_v is the same definition, but for oil vaporized into the gas phase.

Blackoil implementation

The StandardBlackOilSystem implements the black-oil equations. It is possible to run cases with and without water, with and without $R_s$ and with and without $R_v$. The primary variables for the most general case is the reference Pressure, an ImmiscibleSaturation for the aqueous phase and the special BlackOilUnknown that will represent either $S_o$, $R_s$ or $R_v$ on a cell-by-cell basis depending on what phases are present and saturated.

A full description of the black-oil equations is outside the scope of this documentation. Please see [1] for more details.

Compositional: Multi-phase, multi-component flow

JutulDarcy.MultiPhaseCompositionalSystemLV Type

MultiPhaseCompositionalSystemLV(equation_of_state)
MultiPhaseCompositionalSystemLV(equation_of_state, phases = (LiquidPhase(), VaporPhase()); reference_densities = ones(length(phases)), other_name = "Water")

Set up a compositional system for a given equation_of_state from MultiComponentFlash with two or three phases. If three phases are provided, the phase that is not a liquid or a vapor phase will be treated as immiscible in subsequent simulations and given the name from other_name when listed as a component.

source

The more general case of multi-component flow is often referred to as a compositional model. The typical version of this model describes the fluid as a system of $N$ components where the phases present and fluid properties are determined by an equation-of-state. This can be highly accurate if the equation-of-state is tuned for the mixtures that are encountered, but comes at a significant computational cost as the equation-of-state must be evaluated many times.

JutulDarcy implements a standard compositional model that assumes local instantaneous equilibrium and that the components are present in up to two phases with an optional immiscible phase added. This is sometimes referred to as a "simple water" or "dead water" description. By default the solvers use MultiComponentFlash.jl to solve thermodynamic equilibrium. This package implements the generalized cubic approach and defaults to Peng-Robinson.

Assume that we have two phases liquid and vapor referred to as $l$ and $v$ with the Darcy flux given as in the preceeding sections. We can then write the residual equation for each of the $M$ components by the liquid and vapor mole fractions $X_i,Y_i$`of that component as:

$$r_i=\frac{\partial }{\partial t}(({\rho }_l{X}_i{S}_l+{\rho }_v{Y}_i{S}_v)\varphi )+\mathrm{\nabla }\cdot ({\rho }_l{X}_i{\overrightarrow{v}}_l+{\rho }_v{Y}_i{\overrightarrow{v}}_v)-Q_i,M\in \{1,\dots ,M\}$$

For additional details, please see Chapter 8 - Compositional Simulation with the AD-OO Framework in [2].

Compositional implementation

The MultiPhaseCompositionalSystemLV implements the compositional model. The primary variables for the most general case is the reference Pressure, an ImmiscibleSaturation for the optional immiscible phase and $M-1$ OverallMoleFractions.

Multi-phase thermal flow

Thermal effects are modelled as an additional residual equation that comes in addition to the flow equations.

$$r_i=\frac{\partial }{\partial t}({\rho }_r{U}_r(1-\varphi )+\sum _{\alpha }{\rho }_{\alpha }{S}_{\alpha }{U}_{\alpha }\varphi )+\mathrm{\nabla }\cdot (\sum _{\alpha }(H_{\alpha }{\rho }_{\alpha }{v}_{\alpha }-S_{\alpha }{\lambda }_{\alpha }\mathrm{\nabla }T)-{\lambda }_r\mathrm{\nabla }T)-Q_e$$

JutulDarcy.add_thermal_to_model! Function

add_thermal_to_model!(model::MultiModel)

Add energy conservation equation and thermal primary variable together with standard set of parameters to existing flow model. Note that more complex models require additional customization after this function call to get correct results.

source