COMSOL 6.3 - Two-Phase Flow in Column

The following example analyzes two-phase flow in porous media. Describing how immiscible fluids move through porous media is key to answering many environmental and industrial questions. Unfortunately, multiphase analyses are complicated by the need to solve for multiple dependent variables along with a variety of unknowns. Among them are hydraulic properties that depend on the pressure and saturation levels of each fluid phase.

This problem demonstrates two-phase flow following a U.S. Environmental Protection Agency experimental setup (Ref. 1). This straightforward experiment matches observations for a laboratory column to numerical estimates of two-phase flow. With these column experiments, the researchers evaluated flow for varying fluid pairs (air-water, air-oil, and oil-water) and then match the experimental results to those from computer simulations that employ analytic expressions for retention and permeability. This discussion addresses their work for the Lincoln soil and use formulas from Mualem (Ref. 2) and van Genuchten (Ref. 3) to give hydraulic properties.

Figure 1: Geometry of the two-phase flow column experiments in Hopmans and others (Ref. 1).

This is a multipart example. The first part sets up the two-phase flow model for water and air; the equations solve for pressures. Saturation varies with the solution. An underlying assumption is that at least some residual air and water exist throughout the soil column at all times. The model tracks the gas front as it displaces a wetting fluid by observing saturation rather than assuming a discrete interface. The second part modifies the air-water simulation for air-oil and oil-water mixtures.

Model Definition

In the experimental setup for air and water, air enters from the top surface of a column made of water and sand. The incoming air (the nonwetting phase) forces the water (the wetting phase) toward the outlet at the base of the column. At the inlet, air pressure increases by steps in time, and no water exits through the column top. In moving to the outlet, the water passes through a disc that is impermeable to airflow. Neither the air nor the water can pass through the vertical column walls. The water pressure at the outlet, which changes in time, corresponds to the height of fluid rise in a receiving buret. The column has a total height of 8.34 cm, a 6-cm radius, and the disk is 0.74 cm thick. The experiment covers 170 hours.

Governing Equations and Boundary Conditions

Two-phase flow in porous media follows separate equations for the wetting (w) and nonwetting (nw) fluids:

(1)

(2)

where θs is the total porosity or saturated volume fraction; θr is the residual volume fraction, so the difference θs − θr is the available pore space for phases to move; Se is the effective saturation; t is time; κint is the intrinsic permeability of the porous medium (m2); kr is the relative permeability (a function of saturation for a given fluid); μ is the fluid’s dynamic viscosity (Pa·s); p is pressure (Pa); ρ is the fluid density (kg/m3); g is acceleration of gravity; and D is the coordinate (for example, x, y, or z) of vertical elevation (m).

When considering a continuous fluid field, neither phase ever completely fills the pore space, giving a volume fraction for the wetting phase, θw, and nonwetting phase, θnw. For the wetting phase, θ varies from zero or a small residual value θr to the total porosity, θs. The effective saturation, Se, comes from scaling θ with respect to θs and θr and so varies from 0 to 1. Both θ and Se are functions of the pressures of both fluids in the system.

The capillary pressure pc is commonly defined as the difference between the pressure of the nonwetting and wetting phases

(3)

The available pore space can be completely filled with one fluid at a given time, which relates the effective saturations for each phase

(4)

The specific capacity of the wetting phase

depends on changes in the effective saturation with respect to the capillary pressure as

(5)

in the same way, the specific capacity of the nonwetting phase Cp,nw is defined with the help of Equation 4 as

Since the specific capacity of the two phases is the same but with opposite sign, it is just denoted as Cp.

Using Equation 3, Equation 4, and Equation 5 in Equation 1 and Equation 2 simplifies the numerical model. The governing equations become:

You can solve this system of equations simultaneously for pw and pnw. In this example, the two fluids are incompressible, but that need not be the case. Rearranging terms, and adding the density of each fluid, we obtain two coupled Darcy’s law equations, one for the wetting phase, and another for the nonwetting phase

Initially, the water and air in the column follow hydrostatic distributions. The boundary conditions allow the water to exit only from the base of the soil column. For the wetting phase, the boundary conditions are

where n is the unit vector normal to the boundary.

Because air enters at the column top but never exits, the boundary conditions for the nonwetting phase are

Retention and Permeability Relationships

You can set up this two-phase flow analysis using interpolation from experimental data, arbitrary mathematical formulas, or results from other equations in the model to define how θ, C, Se, kr, and pc vary simultaneously. The existing model uses retention and permeability relationships from Ref. 2 and Ref. 3 that express changes in θ, C, Se, and kr as a function of the capillary pressure pc. Because pc is large and because changes in θ, C, Se, and kr are small, these expressions transform capillary pressure to the equivalent height of water or capillary pressure head as in Hc = pc/(ρwater g). The hydraulic properties relative to the wetting fluid in van Genuchten retention model are

where α, n, m, and L denote soil characteristics. Note that with two-phase flow, the van Genuchten–Mualem formulas hinge on the value of Hc.

For the nonwetting fluid, the properties

arise from the definitions for the wetting phase.

Different Fluid Pairs

When switching between air-water, air-oil, and oil-water experiments, the authors used different scaling with interfacial tensions according to Leverett (Ref. 4). The Leverett scaling adjusts the nonwetting phase pressure at the column top to produce the same volume of wetting fluid outflow at the column bottom regardless of the fluid pair. With Leverett scaling, switching between fluid pairs requires using the correct fluid properties ρ and μ for the fluid pair and adjusting the boundary and initial pressures according to

In these equations, σ represents the interfacial tension between the different fluids, and the subscripts denote the fluid pair. These values appear in a table at the end of this section. For example, σao/σaw equals 0.373, and σow/σaw equals 0.534; further, the first nonwetting phase pressure head (in meters of water) is 0.4 m for the air-water system, 0.1 m for the air-oil system, and 0.2 m for the water-oil system.

Because relative permeability and retention properties for a porous medium depend on the fluid moving through it, switching fluid pairs also requires switching the retention and permeability properties in the model. This requirement can mean inserting new experimental data or adjusting mathematical formulas. In this model, the authors assessed the permeability and retention parameters were assessed by curve fitting to analytic formulas. They adjusted the parameters α, n, κs, and θr to get the best fit for each fluid. A review of the data tables that follow reveals that the ratios in the α values for the different fluid pairs roughly equals the σ ratios just given.

Implementation: Step Change on a Boundary

The following step-by-step instructions define the timing of the stepped nonwetting phase pressures at the inlet by using an interpolation function. Interpolating in COMSOL Multiphysics is straightforward. Add an function under , set up the table with the times and corresponding pressure heads, assign a name to the interpolation function, and use it where the function is needed. To activate the functions created, simply enter the function name (for example, Hp_nw_t) along with the argument — that is, the time t in parentheses. For example, you define the nonwetting pressure as

p_nw = Hp_nw_t(t)*rho_water*g_const

The density of water appears in the equation because Ref. 1 defines the boundary pressure as a height of water.

Figure 2: Interpolation function used for Hp_nw_t(t).

Data

The data used in this model correspond to the air-water experiments for the Lincoln sand as reported in Ref. 1:


Variable	Expression	Description
ρw	1000 kg/m3	Fluid density, water
μw	1·10-3 Pa·s	Dynamic viscosity, water
ρa	1.28 kg/m3	Fluid density, air
μa	1.81·10-5 Pa·s	Dynamic viscosity, air
κint	2480 mD	Intrinsic permeability, column
κs	13.57 mD	Permeability, disc
θs	0.32	Saturated volume fraction, columnTab
θs,u	0.5	Saturated volume fraction, disc
pnw,top	0.2 m water	Initial pressure head at inlet, nonwetting phase

The van Genuchten parameters for the different fluid pairs are the following:


Variable	Description	Unit	Air-Water	Air-Oil	Oil-Water
qr,w	Residual volume fraction		0.021	0.00001	0.0072
α	alpha parameter	m-1	1.89	5.29	3.58
n	n parameter, column		2.811	3.002	3.1365
L	L parameter, column		0.5	0.5	0.5
κs	Permeability, disc	m-2	2.48·10-12	1.09·10-12	0.94·10-12

Pressure head at the air inlet increments in time according to the table below, see Figure 2.


Pressure head (m water)	Start time (hours)
0.4	0
0.6	21.25
0.8	45.25
1.0	69
1.5	93
2	122.5
4	155

At the water outlet, the fluid level in the receiving buret increases linearly in time from 0 m to 0.1 m.

Results and Discussion

Figure 3 shows an early-time snapshot from the COMSOL Multiphysics solution for two-phase flow in a laboratory column. The shading depicts the effective saturation of the nonwetting phase (air), while the arrows give the wetting phase (water) velocities.

Figure 3: Solution to two-phase flow model at 0.1 hours: nonwetting phase saturation (surface plot). Results correspond to air-water experiment on Lincoln soil from the US EPA.

The image illustrates the nonwetting fluid entering the soil column and displacing the wetting fluid. The nonwetting phase enters because it is being forced into the inlet with a multistep pressure change.

Figure 4 shows the stepped pressure head used at the inlet boundary along with the capillary pressure in the column at various elevations. Specifying the point locations during postprocessing circumvents the need to plan observation sites during input. The solution to the two-phase flow problem provided is an excellent match to the results of Ref. 1.

That the capillary pressure head and the air inlet pressure in Figure 4 track together is what made the laboratory setup successful. To get high resolution on the permeability and retention behaviors, the authors in Ref. 1 set the pressure steps large enough that the impact is instantaneous in the soil column. As shown in Figure 5, the permeability changes instantaneously throughout the column.

Figure 4: Inlet air pressure head (solid line) and capillary pressure head (dashed lines) for air-water flow in Lincoln soil.

Figure 5: Relative permeability functions for water (solid lines) and air (dashed lines) for Lincoln soil at three different points.

Solutions for two-phase flow for the air-oil and oil-water systems appear in Figure 6 and Figure 7, respectively.

Figure 6: Inlet-air pressure head (solid line) and capillary pressure head (dashed lines) for air-oil flow in Lincoln soil.

Figure 7: Inlet-air pressure head (solid line) and capillary pressure head (dashed lines) for oil-water flow in Lincoln soil.

The COMSOL Multiphysics results for the air-oil and oil-water two-phase flow problems prove to be in agreement to the results shown in Ref. 1. Through Leverett scaling you set the inlet pressure so that the air-oil and oil-water systems would produce the volume outflow rate from the air-water experiment. As with the air-water system, the capillary pressure head and air-inlet pressure for the air-oil experiment track instantaneously. For the water-oil system, however, there is a lag between the nonwetting and wetting phase pressures.

References

1. J.W. Hopmans, M.E. Grismer, J. Chen, and Y.P. Liu, Parameter Estimation of Two-fluid Capillary Pressure Saturation and Permeability Functions, U.S. Environmental Protection Agency EPA/600/R-98/046, Cincinnati, Ohio, 1998.

2. Y. Mualem, “A new model for predicting the hydraulic permeability of unsaturated porous media”, Water Res. Research, vol. 12, pp. 513–522, 1976.

3. M.Th. van Genuchten, “A closed-form equation for predicting the hydraulic of conductivity of unsaturated soils”, Soil Sci. Soc. Am. J., vol. 44, pp. 892–898, 1980.

4. M.C. Leverett, “Capillary behavior in porous solids”, Trans. AIME, vol. 142, pp. 152–169, 1941.

Subsurface_Flow_Module/Fluid_Flow/twophase_flow_column

Modeling Instructions

From the menu, choose .

New

In the window, click