Modeling Non-Newtonian Flow in Porous Media

Non-Newtonian fluids, like ketchup or blood, show strain-rate-dependent behavior, complicating flow modeling through porous media. The complex structure of porous materials — channels, stagnation zones, and isolated pores — causes fluid properties to vary. Without a universal theory, models are often tailored to specific fluid material combinations using measurements or pore-scale simulations. This blog demonstrates how to use these simulations or measurements to develop a homogenized approach for modeling non-Newtonian flow in porous structures.

The Non-Newtonian Fluid Way

From everyday examples, we know that some fluids defy our intuition. A well-known example is ketchup, which will first stubbornly stay in the bottle, and then all of a sudden, the whole plate is covered in red. This phenomenon, known as shear thinning, occurs when viscosity decreases with increasing shear rate. Similarly, most polymers exhibit this behavior. Less common is a shear thickening fluid — where viscosity increases with shear rate. A famous example is Oobleck (a cornstarch suspension), which allows a person to walk on its surface if moving quickly enough.

Shear rate describes how quickly adjacent layers of fluid move relative to each other and depends on the velocity gradient, influenced by the fluid’s velocity \mathbf{u} and the curvature of the flow path. Mathematically, it’s expressed as:

with the strain rate tensor \mathbf{S} = \frac{1}{2}\left(\nabla \mathbf{u} + (\nabla \mathbf{u})^T\right). In porous media, the complex network of diverging and converging pores creates rapidly changing velocities and shear rates. This leads to significant variations in viscosity within the pore space, as illustrated in the following example where a Carreau fluid flows through an irregular pore structure.

Streamlines show the viscosity of a Carreau fluid. Narrow channels with high shear rates exhibit much lower viscosity, while wider regions with lower shear rates maintain higher viscosity, highlighting the influence of pore structure on fluid behavior.

These dynamic viscosity variations pose challenges for modeling non-Newtonian fluid flow. To address these challenges at larger scales, we need effective upscaling methods that translate pore-scale behavior into macroscopic flow characteristics.

Upscaling from Pore-Scale Models

Modeling the flow of non-Newtonian fluids at the pore scale is not feasible for large-scale applications, so upscaling is necessary. One approach involves defining an apparent viscosity \mu_\textrm{app} — essentially the (constant) viscosity that results in the same pressure drop as the non-Newtonian fluid (Ref. 1).

From Darcy’s Law, the apparent viscosity is:

(1)

Here, \kappa is the permeability of the porous medium, v is the velocity, and \nabla p the pressure gradient. These values can be either measured or, as in our case, obtained from pore-scale models of an RVE (representative volume element). The (intrinsic) permeability is a property of the porous matrix and only depends on the pore size, shape, and connectivity. We can compute it as described in this blog post and also shown in this model.

Now, we know the permeability of the structure and let a Carreau fluid flow through it. Its apparent or effective viscosity is defined as

(2)

where the parameters \mu_0 and \mu_\textrm{inf} are the viscosities at zero and infinite shear rate, respectively, \lambda is the relaxation time, and n, the power index. To further justify the need for an upscaling method, we send the same Carreau fluid through a different porous structure that has the same porosity and permeability, and we apply the same pressure drop. As we have learned, the shear rate is influenced by the curvature of the flow path. Since the flow path differs in these different structures, so does the shear rate and hence, what we would measure or compute using (1).

Different porous structures with the same porosity, permeability, and under the same boundary conditions, show varying apparent viscosities, highlighting the influence of both fluid and structure on the shear rate.

So for upscaling of a porous structure, we need an expression that we can use for the shear rate because we cannot directly compute it from the Darcy velocity alone. Let’s call it the apparent shear rate \dot\gamma_\textrm{app} that can then be used in (2). 

There are different formulations that are suitable for different situations. A commonly used formulation is:

(3)

Here, the influence of the pore structure is summarized in the correction factor \alpha. While existing relations for \alpha can be used, more accurate results come from fitting \alpha to measurements or simulation data.

Calculating and Validating the Correction Factor

The correction factor \alpha can be determined using a systematic approach in COMSOL^®. First, compute the apparent shear rate \dot\gamma_\textrm{app} using (2). Also compute the term \frac{|\mathbf{u}|}{\sqrt{\kappa\epsilon_\textrm{p}}}, which we call the normalized velocity. Combine these to compute \alpha from (3). While \alpha can often be assumed constant over a range of pressure drops, it is important to note that it may vary due to changes in flow patterns, such as the onset of turbulence.

In our model, pressure drops ranging from 50 to 5000 Pa/m were simulated. COMSOL^® outputs these results in a table. This table is then used in the Least-Squares Fit function to fit (3) to the data and determine the value of \alpha.

It is also possible to use experimental measurements in the Least-Squares Fit function instead of relying on a pore-scale model.

To validate this approach, we can build a homogenized version of the model based on the Brinkman equations and compare the results. If the outcomes align, this confirms the accuracy and reliability of the approximation.

The Brinkman Equations interface and Darcy’s Law interface now support non-Newtonian flow modeling using the apparent shear rate method including thermal effects. This allows you to choose from six common non-Newtonian relationships. Additionally, you can either input a custom correction factor or select the capillary bundle approach, which offers an alternative description of \alpha suitable for porous media resembling capillary bundles. The figure below shows a screenshot of the COMSOL Multiphysics^® user interface, displaying the settings window for fluid properties.

The results shown below demonstrate that the apparent shear rate method provides accurate and reliable predictions.

Apparent viscosity over normalized velocity for pore-scale model and homogenized approach.

Next Step

Modeling non-Newtonian flow in porous media has unique challenges due to the complex interactions between fluid properties and pore structure. By combining pore-scale simulations or measurements with upscaling techniques, it is possible to derive accurate macroscopic models. The apparent shear rate method, that is now supported in the porous media flow interfaces, provides a reliable approach to account for dynamic viscosity variations and structural influences.

Looking to model the homogenization of non-Newtonian porous media flow yourself? The MPH file and step-by-step instructions are available in the Application Gallery:

In an upcoming blog post, we will explore a real-life example of blood flow through a stenosed pulmonary artery and demonstrate how the apparent shear rate method is applied in this context.

Reference

1. N. Zamani, I. Bondino, R. Kaufmann, A. Skauge, “Computation of polymer in-situ rheology using direct numerical simulation”, Journal of Petroleum Science and Engineering, Volume 159, 2017; https://doi.org/10.1016/j.petrol.2017.09.011.