Model concepts for two-phase flow

2.4 Pore-network modeling 35

For the square cross-section, only Eq. (2.49, Blunt A) deviates clearly from the single-phase value at r_am/r_ij = 0, while it exactly recovers the single-phase value at r_am/r_ij = 1, as does Eq. (2.51, Bakke). The observations are similar for the triangu-lar throat cross-section. Here, Eq. (2.51, Bakke) overestimates the single-phase value at r_am/r_ij = 0. While Eq. (2.52, Valvatne) and Eq. (2.53, Blunt B) exactly repro-duce the single-phase conductance for both geometries at r_am/r_ij = 0, the values for r_am/r_ij →1, corresponding to the single-phase conductance of a circular tube, are over-and underestimated, respectively.

0.00 0.25 0.50 0.75 1.00 ram/rij [-]

0.6 0.7 0.8 0.9 1.0

normalizedgn,ij[-]

0.00 0.25 0.50 0.75 1.00 ram/rij [-]

0.4 0.6 0.8 1.0 Eq. (2.49, Blunt A): _8µ^πreff4

nl_ij

Eq. (2.52, Valvatne): ^kA_µ^totAnG

nl_ij

Eq. (2.51, Bakke): ^r_8µ^2effAn

nl_ij

Eq. (2.53, Blunt B): ^kA_µ ^2nG

nl_ij

Figure 2.7 –Non-wetting phase conductances. Comparison of normalized non-wetting phase conductances for a square (left) and equilateral triangle (right) throat cross-section as given by Eqs. (2.49) to (2.53).

described in the literature. Such pore-network models can be generally divided into quasi-static and dynamic models.

Quasi-static PNM The first type of pore-network models is based on the assumption that the fluids within the pore space are at rest or, in other words, that there is no pres-sure gradient within the individual fluid phases. This implies that there is a constant capillary pressure throughout the entire fluid-fluid interface which corresponds to an externally applied global capillary pressure. A drainage or imbibition process can be modeled by incrementally increasing or decreasing this global pressure difference which will result in discrete saturation changes, corresponding to equilibrium states. Since fluid flow is not considered, only the critical threshold capillary pressures described in Section 2.4.2 are relevant for quasi-static models. Fluid configurations between the equilibrium states can not be modeled.

Drainage corresponds to an invasion percolation process where the bonds (throats) of the network are filled subsequently based on the ascending order of the entry capillary pressures, under the constraint that only throats adjacent to already filled throats or the inlet of the domain may be invaded. Invasion percolation also holds for imbibition with piston-type advance but here, the sites (pore bodies) are the limiting elements which means that smaller pores and their neighboring throats will be filled first. In contrast to invasion percolation, percolation allows a pore filling throughout the entire domain, regardless whether the newly filled element is connected to an already filled one or not.

This corresponds to imbibition with layer-swelling which assumes a globally connected wetting phase via the pore throat corners [Blunt, 2017]. Quasi-static models have been used successfully for the prediction of rock properties such as capillary pressure-saturation curves [e.g., Oren et al., 1998, Valvatne and Blunt, 2004], but they are limited to situations where equilibrium conditions can be safely assumed or capillary forces dominate the system. This is given for low capillary numbers,

Ca= µv

γ , (2.54)

where µ and v are the dynamic viscosity and the specific bulk flow rate, i.e., Darcy velocity (Section 2.3), of the invading phase and γ is the interfacial tension.

2.4 Pore-network modeling 37

Dynamic PNM Dynamic pore-network models account for time-dependent phase dis-placement processes and are capable of modeling non-equilibrium capillary pressure states under the consideration of the fluid phases’ viscosities. Dynamic models solve for pressure fields and phase fluxes based on the conservation of mass or volume, in close analogy to conventional Darcy-type models. Given their highly non-linear be-havior and the requirement to capture both very slow and fast displacement processes, dynamic models are computationally more complex and resource intensive than quasi-static models [Blunt, 2017]. Joekar-Niasar and Hassanizadeh [2012a] extensively review the different types of solution strategies for dynamic two-phase pore-network models.

They discern two general types of dynamic models: the first one assigns a single pres-sure to each pore body [e.g., Koplik and Lasseter, 1985, Mogensen and Stenby, 1998, Al-Gharbi and Blunt, 2005], assuming either the exclusive presence of a single phase or the concept of an equivalent pressure that accounts for both phases. While this ap-proach decreases the computational complexity of the problem, Al-Gharbi and Blunt [2005] observed inconsistencies with respect to equivalent quasi-static simulation re-sults for networks with angular cross sections. The second type of algorithm solves for an individual pressure field for each phase. This type of two-pressure model was first introduced by Thompson [2002], following a sequential solution strategy based on a decoupling between pressure and saturation and thus a faster solution of the linearized problems. This is very similar to the IMPES (implicit pressure, explicit saturation) method widely used in reservoir models [e.g., Chen et al., 2004]. Having solved for the pressures, the saturation field is updated explicitly, which requires very small time steps in order to maintain numerical stability. However, no close match to quasi-static results could be achieved and the method was unsuited for very low Ca. This was improved by Joekar-Niasar et al. [2010a] who introduced a semi-implicit saturation up-date in order to account for the highly non-linear nature of the processes. The model has found application in a number of other works [e.g., Qin, 2015, Khayrat and Jenny, 2017]. However, the efficiency of the model is still limited by the maximal admissible time step size which is governed by the filling dynamics of the pores at the saturation front [Bierbaum, 2019].

3 Conceptual and numerical models for the individual sub-domains

This chapter explains the mathematical models and the numerical discretization schemes used for the individual sub-domains. It is followed by a description of the interface con-ditions required for the coupling of the sub-models in Chapter 4 and an elaboration of the coupled model’s implementation in the open-source simulator DuMu^x (see Chap-ter 5). As shown in Fig. 3.1, the coupled model comprises up to three sub-domains: the free-flow domain Ω^FF, the interface region of the porous medium Ω^PNMand an optional bulk porous medium domain Ω^REV with low process activity.

Figure 3.1 – Conceptual coupled model. Sketch of the coupled model comprising the three sub-domains Ω^FF, Ω^PNM and Ω^REV. The latter domain is optional. The domains are coupled at the interfaces Γ^FF and Γ^REV. δ and Lare length scales.

Ω^PNMrepresents a “thin” interface region separating Ω^FF and Ω^REV, henceδL. The thickness δ strongly depends on the coupled system under consideration and should

be determined experimentally or by means of a high-resolution simulation. It should probably at least cover the extent of an REV to meet the requirements of a “non-simple interface” [Hassanizadeh and Gray, 1989]. In future work, the possibility of estimating δbased on characteristics of the free flow (such as the boundary layer thickness) and the bulk porous medium (such as permeability or porosity) could be assessed. We assume that δ is in the order of several pore diameters.

The bulk porous medium domain Ω^REV is optional. As explained later in Section 4.2, formulating physically consistent coupling conditions between Ω^PNM and Ω^REV is chal-lenging and not the focus of this work where we rather concentrate on the coupling between Ω^FF and Ω^PNM. We nevertheless present a first, somewhat tentative set of coupling conditions for the two porous-media models in the hope of providing a foun-dation on which future work can build.

The conceptual models and numbers of fluid phases considered in each sub-domain are:

ˆ Ω^FF: (Navier-)Stokes, single-phase flow (1p)

ˆ Ω^PNM: pore-network model, single-phase or two-phase flow (1p or 2p)

ˆ Ω^REV: Darcy, single-phase or two-phase flow (1p or 2p)

Component transport and non-isothermal flow can be considered for all sub-domains.

The fluid phases’ physical properties, such as density and viscosity, can be solution-dependent. The models are implemented in a way such that the influence of gravity can be incorporated or neglected. For all numerical results presented in this work, gravity is neglected.

Furthermore, we always consider laminar flow regimes in all sub-domains with creeping flow in the porous medium (Ω^FF and Ω^REV) while Recan be greater than one in Ω^FF. A monolithic coupling approach (see Chapter 5) is followed in this work and the fully implicit backward Euler method is used for discretizing the coupled model in time.

We therefore do not address temporal discretization for the individual sub-models.

We employ a Newton-Raphson scheme to solve the potentially non-linear system of equations.

3.1 (Navier-)Stokes model (free flow) 41

Treatment of solution-dependent parameters (secondary variables) All balance equations in the following sections are formulated in a general form such that fluid properties like density%, viscosityµ, specific internal energyu, specific internal enthalpy h, heat conductivity λ or diffusion coefficients D do not have to be constant but may depend on the system’s state, i.e., the current values of the respective primary variables.

We will specify for each numerical example shown in Chapters 7 and 8 whether constant or solution-dependent quantities are used.

In order to avoid repetition, we briefly describe here how these solution-dependent pa-rameters (secondary variables) are treated. We employ the fluidsystem framework of DuMu^x [Lauser, 2012] (see Chapter 5) which provides all necessary constitutive re-lations (see Chapter 2) for a wide range of different substances and mixtures. These may be either equations of state such as theIdeal Gas Law (Eq. (2.7)) or other supple-mentary empirical equations as well as tabulated values such as the industry standard IAPWS [Wagner and Pruß, 2002] for water.

3.1 (Navier-)Stokes model (free flow)

As explained above, only one single phase may be present in Ω^FF. If, for instance, air flow over a partly water-saturated soil is considered, water may only enter the free-flow region in gaseous form. For brevity, we neglect the phase subscript α in this section.