Two-Phase Flow over a Low Permeable Lens

This example concerns two-phase flow in a porous medium that contains a low permeable lens. The heavier phase infiltrates the porous medium from above, and the low permeable lens is infiltrated only when a critical saturation at the outside of the lens is reached. As the saturation of the heavier phase is discontinuous at the boundary of the lens, this requires the use of the Porous Medium Discontinuity boundary condition.

Model Definition

The porous domain is assumed to be axially symmetric, with a radius of 0.5 m and a height of 0.65 m. The low permeable lens has radius of 0.32 m and a height of 0.12 m. The bottom boundary of the lens is located at a height of 0.35 m. Initially the porous domain, including the lens, is occupied with phase 1. Phase 2 flows into the porous medium at the top boundary through a circle with radius 0.07 m with a uniform and constant mass flux. See Figure 1 below for a graphic representation of the geometry.

Figure 1: Cross section of the axially symmetric geometry.

Initially the porous domain, including the lens, is occupied with phase 1. Phase 2 flows into the domain with a constant mass flux. The properties of the two phases are given in Table 1.

Table 1: Fluid properties.
Quantity	Value	Description
ρ1	1000 kg/m3	Density of phase 1
ρ2	1460 kg/m3	Density of phase 2
μ1	10-3 Pa·s	Dynamics viscosity of phase 1
μ2	0.9·10-3 Pa·s	Dynamics viscosity of phase 2

The properties of the solid matrix and the parameters for the constitutive relations for the relative permeabilities and capillary pressure curves, which are described by the Brooks and Corey model, are given in Table 2.

Table 2: Solid matrix properties and brooks & Corey parameters.
Quantity	Value in lens	value	Description
εp	0.39	0.4	Porosity
κ	3.32·10-11 m2	6.64·10-11 m2	Permeability
sr1	0.12	0.1	Residual saturation of phase 1
sr2	0	0	Residual saturation of phase 2
λp	2	2.7	Pore size distribution inde x
pec	1163.5 Pa	775 Pa	Entry capillary pressure

The initial values for the saturation of phase 1 and the pressure of phase 2 are given in Table 3.

Table 3: Initial values.
quantity	value
s2	0
p	(0.65-z)g_const1000[kg/m^3]

The boundary conditions are given in Table 4. In this table q0,si denotes the normal mass flux of phase i. The number of the boundaries refer to the numbers indicated in Figure 1. The time interval for the simulation is 100 minutes.

Table 4: boundary conditions.
boundary	condition
1,3,5	axial symmetry
2	s2 = 0, q0,s1 = 0
7	q0,s1=0, q0,s2=0.25 kg/(m2·s)
8	q0,s1=0, q0,s2=0
10	s2 = 0, p = (0.65-z)g_const1000[kg/m^3]

Results and Discussion

Due to gravity, the heavier phase 2 infiltrates the porous domain and flows down over the low permeable lens. Since the entry capillary pressure of the lens is higher than the entry capillary pressure of the surrounding material, phase 2 will not enter the lens directly when it reaches the lens. Phase 2 will only enter the lens when a critical saturation is reached. This condition, which applies at boundaries where the porous medium properties, and especially the capillary pressure curves, are discontinuous, is implemented in the model using a Porous Medium Discontinuity boundary condition. This condition allows for a discontinuity in the saturation of phase 2 and determines the critical saturation at which phase 2 enters the low permeable domain. The figures below show that this happens after around 12 minutes. After approximately 60 minutes, phase 2 has reached the bottom of the lens.

This simulation is inspired by a very similar model as discussed in Ref. 1and Ref. 2.

Notes About the COMSOL Implementation

In the present implementation of the model, the dependent variables are the saturation of phase 2, s2, and the pressure of phase 1, p. The equation for the saturation takes as boundary flux the mass flux of phase 2, and the equation for the pressure takes as boundary flux the total mass flux (mass fluxes of phase 1 and 2 added together). The boundary condition at the bottom boundary prescribes the saturation of phase 2 and the mass flux of phase 1. To be able to prescribe the total mass flux in the equation for p, the mass flux of phase 2 is also needed. This mass flux is computed automatically if the saturation condition for phase 2 is implemented as a weak constraint, see the instructions in the Modeling Instructions section.

Figure 2: Isosurfaces of the penetrating phase 2 after 12 minutes. Phase 2 just starts entering the low permeable lens at this instant in time.

Figure 3: Isosurfaces of the penetrating phase 2 after 60 minutes. Phase 2 has now reached the bottom of the low permeable lens.

References

1. R. Helmig, Multiphase Flow and Transport Processes in the Subsurface – A Contribution to the Modeling of Hydrosystems, Springer Verlag, 1997.

2. P. Bastian, Numerical Computation of Multiphase Flows in Porous Media, Habilitationsschrift Universität Kiel, 1999.

Porous_Media_Flow_Module/Fluid_Flow/low_permeable_lens

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.5.

In the text field, type 0.65.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.32.

In the text field, type 0.12.

Locate the section. In the text field, type 0.35.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.07.

In the text field, type 0.65.

Click

Phase Transport in Porous Media (phtr)

In the window, under click .

In the window for , locate the section.

Select the check box.

Phase and Porous Media Transport Properties 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the pec text field, type 1163.5.

Locate the section. From the ρs1 list, choose . From the μs1 list, choose . In the srs1 text field, type 0.12.

Locate the section. From the ρs2 list, choose . In the associated text field, type 1460[kg/m^3].

From the μs2 list, choose . In the associated text field, type 0.0009[Pa*s].

Phase and Porous Media Transport Properties 2

In the toolbar, click

and choose .

Select Domain 1 only.

In the window for , locate the section.

From the list, choose .

In the pec text field, type 755.

In the λp text field, type 2.7.

Locate the section. From the ρs1 list, choose . From the μs1 list, choose . In the srs1 text field, type 0.1.

Locate the section. From the ρs2 list, choose . In the associated text field, type 1460[kg/m^3].

From the μs2 list, choose . In the associated text field, type 0.0009[Pa*s].

Mass Flux 1

In the toolbar, click

and choose .

Select Boundary 7 only.

In the window for , locate the section.

Select the check box.

In the q0,s2 text field, type 0.25.

Volume Fraction 1

In the toolbar, click

and choose .

Select Boundaries 2 and 10 only.

In the window for , locate the section.

Select the check box.

Click the

button in the toolbar.

In the dialog box, in the tree, select the check box for the node .

Click .

In the window for , click to expand the section.

Select the check box.

Porous Medium Discontinuity 1

In the toolbar, click

and choose .

Select Boundaries 4, 6, and 9 only.

Darcy’s Law (dl)

In the window, under click .

In the window for , locate the section.

Select the check box.

Gravity 1

In the Gravity node specify a reference position to ensure that the pressure is zero at the upper boundary of the model.

In the window, under click .

In the window for , locate the section.

Select the check box.

Specify the rref vector as


r	r
0.65	z

Porous Matrix 1

In the window, under click .

In the window for , locate the section.

From the εp list, choose . In the associated text field, type 0.39.

From the κ list, choose . In the associated text field, type 3.32e-11[m^2].

Porous Medium 2

In the toolbar, click

and choose .

Select Domain 1 only.

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the εp list, choose . In the associated text field, type 0.4.

From the κ list, choose . In the associated text field, type 6.64e-11[m^2].

Initial Values 1

In the window, under click .

In the window for , locate the section.

Click the button. This way the initial pressure field is forced to equal the hydraulic pressure.

Inlet 1

In the toolbar, click

and choose .

Select Boundary 7 only.

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the N0 text field, type 0.25.

Pressure 1

In the toolbar, click

and choose .

Select Boundary 10 only.

In the window for , locate the section.

In the p0 text field, type (0.65-z)*g_const*1000[kg/m^3] to compensate for hydrostatic pressure.

Mass Flux 1

In the toolbar, click

and choose .

Select Boundary 2 only.

In the window for , locate the section.

In the N0 text field, type s2_lm.

Mesh 1

Free Triangular 1

In the toolbar, click

Size

In the window, click .

In the window for , locate the section.

Click the button.

Locate the section. In the text field, type 0.01.

Study 1

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

In the text field, type range(0,60,6000).

Solution 1 (sol1)

In the toolbar, click

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose . This setting ensures that the variable scaling is based on the supplied initial hydrostatic pressure profile, which in this case gives a better scaling than the automatic setting. Scaling is important to obtain well weighted error estimates and avoid ill-conditioned matrices which may hamper or slow down the solution procedure.

In the toolbar, click