Pore-Scale Flow

This example uses Creeping Flow (Stokes Flow) to solve the flow in the interstices of a porous medium. The example comes from pore-scale flow experiments conducted by Arturo Keller, Maria Auset, and Sanya Sirivithayapakorn of the University of California, Santa Barbara. To produce the example geometry they used electron microscope images. This type of nonconventional pore-scale modeling with COMSOL Multiphysics sheds new light on the movement of large particulates and colloids moving through variable-pore geometries in the subsurface. Several of these researchers have published results from their COMSOL Multiphysics modeling in the publication Water Resources Research (Ref. 1 and Ref. 2).

Keller, Auset, and Sirivithayapakorn designed their lab experiments on the basis of scanning electron microscope (SEM) images of thinly sliced rock sections (Figure 1). They etched the geometric patterns from the images onto a solid with an elaborate process similar to the etching of silicon wafers. They then transferred these images to DXF files, which they finally imported into COMSOL Multiphysics.

Figure 1: Scanning electron microscope image of the repeat pattern in the silicon wafer. The scale at bottom indicates that pore throat and body dimensions are on the order of 1 μm–100 μm (Ref. 1).

It is typical to represent fluid flow in the subsurface as a continuum process using average or “continuous” properties for the bulk rather than detailing the shape and orientation of each solid particle within a porous medium. Inserting the bulk properties into an equation, such as Darcy’s law or Brinkman equations, gives an average flow rate for the total volume. While bulk approximations typically produce excellent estimates sufficient for considering flow over large areas, they might miss the between-grain nuances that a close-up Stokes flow analysis would give.

This exercise is divided in two models: The first model takes one of the SEM images of Keller, Auset, and Sirivithayapakorn and solves for the flow velocity and pressure drop in the pore throats using the Creeping Flow interface. The geometry is imported as a binary file and only the pore space is meshed but not the solid regions. The second model is devoted to the modeling of the whole slice, by importing the SEM image and deriving porous medium properties, such as porosity and permeability for further use in the Brinkman Equations interface.

Model Definition

The entire model covers 640 μm by 320 μm. Water moves from right to left across the geometry. The flow in the pores does not penetrate the solid grains. The inlet and outlet fluid pressures are known. Assume no flow at the top and bottom boundaries.

Figure 2: A 640 μm by 320 μm geometry and boundary conditions.

Since the channels are at most 0.1 mm in width and the maximum velocity is lower than 10−4 m/s, the maximum Reynolds number is less than 0.01. Since the Reynolds number is far less than one, the example uses the Creeping Flow interface to solve Stokes equations, instead of solving the Navier–Stokes equations with the Laminar Flow interface. The fluid is considered isothermal and with constant density. Owing to the problem’s small scale, the example does not consider gravity.

The Creeping Flow interface solves Stokes equations in the channels. The incompressible assumption together with the stationary condition reads

here, p is the pressure (SI unit: Pa), u is the velocity field (SI unit: m/s), μ is the dynamic viscosity of the fluid (SI unit: Pa·s), and ρ is the density of the fluid (SI unit: kg/m3).

At the physical boundaries, the inlet pressure and the outlet pressure are known. Velocities are zero at the grain boundaries, which implies a no-slip condition. The flow is symmetric about the top and bottom boundaries. Table 1 summarizes the boundary conditions.

Table 1: Boundary Conditions.
Boundary Type	Boundary Condition	Value
Inlet	Pressure	p = p0
Outlet	Pressure	p = 0
Grain walls	Wall	no slip
Symmetry sides	Symmetry	—

Here p0 is a specified pressure drop. Table 2 collects the relevant model data.

Table 2: Model Data.
Quantity	Value	Description
ρ0	1000 kg/m3	Fluid density
μ0	0.001 kg/(m·s)	Fluid dynamic viscosity
p0	0.715 Pa	Pressure drop

Results and Discussion

Figure 3 shows the COMSOL Multiphysics solution predicted with the creeping flow analysis for the fluid velocity field in the pore spaces of a microscale porous slice. The velocity magnitude is higher in the narrowest pores than at the inlet, tending to decrease in stretches where the channels’ cross-sectional area increases.

Figure 3: Surface and arrow plots of the velocity field calculated by the Creeping Flow interface.

Model Definition

The second model takes a completely different approach than the first model. Here, the scanning electron microscope image is imported and physical properties are derived from the color scale. As opposed to consumer cameras, SEM images are grayscale, but in this example, the color code is binary; see Figure 4.

Instead of solving for the creeping flow in the channels, the incompressible, stationary Brinkman equations with the Stokes–Brinkman assumption are used:

Here, p is the pressure, u is Darcy’s velocity field, μ is the dynamic viscosity of the fluid, εp is the porosity, and κ is the permeability of the medium.

In order to define physical properties from the image color code, the following relations have been implemented for the porosity and permeability:

(1)

(2)

Here, im1 is an image function derived from the SEM image, which in this example ranges from 0 to 1 as a function of position. Other expressions can be implemented when importing RGB or grayscale images.

Figure 4: SEM image. The color code is blue for 0 and red for 1. COMSOL Multiphysics can handle grayscale and RGB images.

After solving the Brinkman equations with material parameters derived from the SEM image, it is possible to observe a very similar pressure and velocity profiles in the porous slice as obtained with the creeping flow model. Compare Figure 5 with the results in Figure 3.

Figure 5: Surface and arrow plots of the velocity field calculated by the Brinkman Equations interface. The porosity and permeability are taken from the SEM image, as written in Equation 1 and Equation 2.

Notes About the COMSOL Implementation

In this model, an external image is used to infer material properties, such as porosity and permeability. The same technique can be used to derive other material properties, like density and thermal or electric conductivity.

To find out more about how to import images, see the chapter about the Definitions node in the COMSOL Multiphysics Reference Manual.

References

1. M. Auset and A.A. Keller, “Pore-scale Processes that Control Dispersion of Colloids in Saturated Porous Media,” Water Resour. Res., vol. 40, no. 3, p. W03503, 2004.

2. S. Sirivithayapakorn and A.A. Keller, “Transport of Colloids in Saturated Porous Media: A Pore-scale Observation of the Size Exclusion Effect and Colloid Acceleration,” Water Resour. Res., vol. 39, no. 4, 2003.

Subsurface_Flow_Module/Fluid_Flow/pore_scale_flow

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

In the window, under click .

In the window for , locate the section.

From the list, choose .

Import 1 (imp1)

In the toolbar, click

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file pore_scale_flow.mphbin.

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
rho0	1000[kg/m^3]	1000 kg/m³	Fluid density
eta0	0.001[kg/(m*s)]	0.001 kg/(m·s)	Dynamic viscosity
p0	0.715[Pa]	0.715 Pa	Pressure drop

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rho0	kg/m³	Basic
Dynamic viscosity	mu	eta0	Pa·s	Basic

Creeping Flow (spf)

Inlet 1

In the window, under right-click and choose .

Select Boundaries 2231 and 2232 only (the straight edges on the far right).

In the window for , locate the section.

From the list, choose .

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

Outlet 1

In the toolbar, click

and choose .

Select Boundaries 1, 4, 7, 10, 13, and 16 only (the straight edges on the far left).

Symmetry 1

In the toolbar, click

and choose .

Select Boundaries 31, 59, 247, 342, 517, 601, 720, 825, 878, 1149, 1300, 1421, 1488, 1748, 1756, 1823, 1912, and 2025 only.

Use the control in the toolbar to select all straight upper boundaries at once. Repeat to select all lower boundaries at once. Alternatively use the button next to the list and enter the numbers listed above.

Mesh 1

In the window, under right-click and choose .

Study 1

In the toolbar, click

Results

Velocity (spf)

The first default plot shows the magnitude of the velocity field. Add a streamline plot to visualize the velocity field with the following steps.

Streamline 1

In the toolbar, click

Select Boundaries 2231 and 2232 only.

In the window for , locate the section.

Find the subsection. From the list, choose .

From the list, choose .

In the toolbar, click

Root

Next, set up the second model using the Brinkman equations.

Add Component

In the window, right-click the root node and choose .

Add Physics

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Find the subsection. In the table, clear the check box for .

Click in the window toolbar.

In the toolbar, click

to close the window.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the table, clear the check box for .

Find the subsection. In the tree, select .

Click in the window toolbar.

In the window, click the root node.

In the toolbar, click

to close the window.

Global Definitions

Parameters 1

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file pore_scale_flow_parameters.txt.

Geometry 2

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type L.

In the text field, type H.

Click

Global Definitions

Define a function that will be used when setting up the material properties.

Image 1 (im1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type L.

In the text field, type H.

Locate the section. Click

Browse to the model’s Application Libraries folder and double-click the file pore_scale_flow_structure.png.

Click

Results

SEM image

In the window for , type SEM image in the text field.

View 2D 3

In the window, expand the node.

Axis

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Click

, and compare with Figure 4.

Materials

The functions for porosity and permeability are defined according to Equation 1 and Equation 2.

Material 2 (mat2)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rho0	kg/m³	Basic
Dynamic viscosity	mu	eta0	Pa·s	Basic
Porosity	epsilon	1-0.99*im1(x,y)	1	Basic
Permeability	kappa_iso ; kappaii = kappa_iso, kappaij = 0	k0/(100*im1(x,y)+0.1)	m²	Basic