Localization in Flow of Non-Newtonian Fluids Through Disordered Porous Media

Introduction

Flow through porous media is of great interest in chemical engineering, physics, and biology [1–3]. Previous studies have shown that the disordered characteristics of the pore structure naturally leads to heterogeneous flow patterns [4–8] and preferential channeling [9–11]. Understanding how to control and manipulate these flow patterns can help to optimize catalysts [12, 13] or chromatographic devices [14, 15], and allows to steer chemical reactions inside the porous medium itself [16–18].

In order to understand the physics of important problems like, for example, blood flow through the kidney [19] or oil flow through porous rocks [20, 21], one must also consider the nonlinear constitutive behavior of the fluids involved in these processes. Technological applications which make use of non-Newtonian fluids are ubiquitous nowadays [22–25]. It is, for instance, the case of shear-thinning solvents that are present in dropless paints [26], shear-thickening fluids being used as active dampers [27] and hybrid fluids as components of enhanced body armors [28]. While Newtonian flows in irregular media have been extensively investigated theoretically and confirmed by many experiments, the study of non-Newtonian fluids lack a generalized framework due to their diverse constitutive nature. Non-Newtonian flows through porous media have mainly been studied theoretically [29, 30] and through numerical simulations [31, 32], where the main focus of interest was to find non-Darcian models for the flow of generalized Newtonian fluids [30, 33–37]. In the particular case of power-law fluids, it has been shown that, in spite of the nonlinear nature of the fluid’s rheology and the geometrical complexity of the pore volume, the general behavior of the system can still be quantified in terms of a universal permeability extending over a broad range of Reynolds conditions and power-law exponents [38].

However, quantitative experiments with non-Newtonian materials which go beyond simple bulk measurements [39, 40] are scarce [41] because the design of the experimental pore geometry and the operating conditions need to be adjusted in order to match the nonlinear constitutive regime of the fluid’s rheology. Here we combine the results of microfluidic experiments [42] with fluid dynamics simulations to demonstrate how the nonlinear rheological properties of a fluid can be effectively exploited in order to control the macroscopic transport properties of a flow through the external operational flow conditions. These results have important consequences for the design of chemical reactors and chromatographic systems as well as for the enhancement of oil recovery and transport in porous media in general.

Under steady-state conditions the motion of an incompressible fluid through the interstitial space of a porous medium is described by mass and momentum conservation, respectively,

together with appropriate boundary conditions. The variables ϱ, $\mathbf{u}$ and p are the fluid’s density, velocity and pressure, and $\mathsf{T}$ is the deviatoric stress tensor which depends on the fluid’s rheology. For many fluids, this constitutive relation is well described by a simple linear rheology $\mathsf{T}=2\mu \mathsf{E}$, where ${\mathsf{E}}_{ij}=1/2\left({\partial }_j{u}_i+{\partial }_i{u}_j\right)$ is the shear strain rate tensor [43] and the proportionality constant μ defines the kinematic viscosity. Examples of these so-called Newtonian fluids are water, light oil and most diluted gases. However, many fluids present in industrial products, biology and environmental flows obey much more complex nonlinear constitutive laws [23, 25, 44]. These fluids are called non-Newtonian fluids. The constitutive behavior of most non-Newtonian fluids can be described by a generalization of the Newtonian relation, namely,

Here the apparent viscosity $\mu \left(\dot{\gamma }\right)$ is a nonlinear function of the second principal invariant $\dot{\gamma }=\sqrt{2\mathsf{E}:\mathsf{E}}$ of the shear strain rate tensor $\mathsf{E}$ alone [43]. Examples of fluids—which are often called generalized Newtonian fluids—are colloidal suspensions, protein or polymeric mixtures, heavy petroleum, blood or debris flows, only to mention a few [45–48].

Our analysis is based on experimental results from the setup presented recently in Eberhard et al. [42]. Their main purpose was to map the local viscosity of a non-Newtonian flow in a porous microfluidic channel by means of a high-resolution technique of image velocimetry, namely, Ghost Particle Velocimetry (GVP) [49]. The geometry of the microfluidic chip is shown in Figure 1. As non-Newtonian fluid we used a 0.5 wt% xanthan gum solution, which is a polysaccharide mainly found in food industry [50] and enhanced oil recovery [22, 24]. It has a shear-thinning rheology which changes its apparent viscosity over several orders of magnitude. While polymeric solutions often show viscoelastic behavior [51], the concentration of xanthan gum in our experimental solution was so low that no measurable elastic behavior could be observed during the experiment. The rheology of xanthan gum closely follows a Carreau model,

approaching the viscosity of the solvent (water) ${\mu }_{\infty }=0.001$$\text{Pa}\cdot \text{s}$ in the limit of very high shear [52]. Conversely, for low shear, Eq. 4 reduces to ${\mu }_C\left(\dot{\gamma }\right)={\mu }_0$, corresponding to a constant viscosity ${\mu }_0=24$$\text{Pa}\cdot \text{s}$. At an intermediate range of shear rates, the fluid follows a power-law relation $\mu \sim {\dot{\gamma }}^{n-1}$ with $n=0.3$ in our specific case. The remaining Carreau parameter $\lambda =50\text{s}$ was determined using a nonlinear least-square fit to the experimentally measured values [42]. As shown in Figure 1, the experimental pore structure consisted of a microfluidic device containing a quasi-2D porous medium of size 30 mm × 15 mm and depth of 100 µm. It contains pillars of radius 100 µm that are randomly allocated and can overlap, forming a “swiss-cheese” pore geometry with void fraction approximately equal to 0.8.

FIGURE 1. Sketch of the experimental setup showing the whole pore geometry together with the simulated flow field in the mid plane for the non-Newtonian case at qin 5 µL/min. The porous region is L 30 mm long, B 15 mm wide and has a height of h 100 µm. The radius of the circular pillars is also 100 mµ. The gray shaded region marks the part of the porous device where the experimental flow velocity measurement has been performed.

Figure 2 compares the velocimetry measurements obtained from ref. 42 in a section of the mid plane of the microfluidic chip with those obtained from numerical simulations calculated with exactly the same pore geometry, fluid properties and flow conditions. More precisely, the flow rates for the presented cases are q_in = 0.05 µL/min (Figure 2A) and q_in = 5 µL/min (Figure 2B) for the xanthan case, and q_in = 5 µL/min for the measurement with water (Figure 2C). The color scale has been normalized to the 95% quantile of the velocity distribution to facilitate the comparison of the flow fields at different flow rates. Although differences between the experimentally measured and simulated velocity fields can be visually detected, they are mostly local and could be explained by the natural difficulties of exactly reproducing in the mathematical model the detailed features of the flow, the fluid rheology, and the flow operational conditions. To perform numerical simulations, the computational mesh was generated by capturing the two-dimensional technical drawing of the device geometry into Ansys’ meshing module [53]. In the horizontal plane, we first created an unstructured quadrilateral mesh with an average cell size of $\approx 2{\text{μm}}^2$. The three-dimensional structure was obtained by extruding ten vertical layers, resulting in a computational mesh composed of approximately 30 million unstructured hexahedral cells (corresponding to roughly 20 µm³/cell), which was then imported into Ansys Fluent^TM [53]. Non-slip boundary conditions were applied on all solid walls of the microfluidic chip, which is a reasonable assumption on surfaces without hyper-hydrophilic coatings and at scales much larger than the polymer coil size [54]. The fluid was injected via a constant velocity inlet corresponding to the inflow rate reported in the experiment. The density of the non-Newtonian fluid used in the computational simulations matches exactly the experimental value of the xanthan solution, namely, ρ_xan = 1.0 g/cm³. As for the rheology of the fluid, we used Eq. 4 to interpolate the local viscosity values obtained from the independent rheometer experiments [42]. Finally, the steady-state flow solution in terms of the velocity and pressure fields was calculated using a second-order integration scheme and convergence was achieved if the residuals reached a threshold of ${10}^{-6}$.

FIGURE 2. Comparison of the experimentally measured [42] and simulated velocity fields in the mid plane of the porous device for the two non-Newtonian cases at (A) qin 0.05 µL/min, (B) qin 5 µL/min, and for the Newtonian case (C). Note the different velocity ranges between (A–C). For comparison, the color scales has been normalized to the 95% quantile of the velocity distribution.

In order to check the variability of our results with respect to the disorder level of the pore space, numerical simulations have also been performed with three additional realizations of the swiss-cheese geometry, but keeping the same physico-chemical properties of the fluid, operational parameters of the flow, and boundary conditions. Considering the independence of the rheometry and velocimetry measurements, the excellent agreement between results from the numerical model and experiments (Figures 2A–C) clearly demonstrates the global consistency of our methodological approach. In order to highlight the differences between the Newtonian and non-Newtonian flows and the tendency for stronger localization in the non-Newtonian flow, we show in Figure 3A the contour plot of the ratio between the local velocity magnitudes measured with the xanthan solution and water normalized by their respective mean velocities. In both cases, the applied flow rate was set to q_in = 5 µL/min.

FIGURE 3. (A) Contour plot of the ratio between the local velocity magnitudes measured with the non-Newtonian (xanthan) and Newtonian (water) fluids, normalized by the corresponding mean velocities in the observation mid planes, both obtained at qin 5 µL/min. (B) Participation ratio as a function of Reynolds number. Gray markers label the participation values obtained from the simulations of four different realizations of the pore structure. The realization which corresponds to the experimental device is marked in darker gray. The black solid line was obtained by averaging the four realizations of the “swiss-cheese” pore geometry. The participation ratio calculated from the experimentally measured velocity fields of the xanthan gum solution at q_in 0.05 µL/min and q_in 5 µL/min are labeled with yellow and blue stars, respectively.

Analysis

The channeling effect present in the flow fields shown in Figure 2 can be statistically quantified in terms of their spatial distributions of kinetic energy $e\propto {\left|\mathbf{u}\right|}^2$. This is performed here in terms of a measure utilized in previous studies on localization of vibrational modes in harmonic chains [55], namely, the participation number $\text{Π}$ defined as,

where the total volume of the fluid in a domain is given by $V={\displaystyle {\int }_{\text{Ω}}1\text{d}\omega }$. Thus the participation ratio varies between $\text{Π}=1$, corresponding to a limiting state of equal partition of kinetic energy $\left(e\left(\mathbf{x}\right)=\text{const}\cdot \forall \mathbf{x}\in \text{Ω}\right)$ and the value $\text{Π}\approx 0$ for a sufficiently large system $\left(V\to \infty \right)$, indicating strong localization, namely, the presence of intense channeling effects in the flow field [56]. For reference, the Reynolds number is defined here as,

where v is the mean velocity at the entrance of the pore zone and $d_p$ is the diameter of the solid obstacles. Figure 3B shows how the participation ratio obtained with our computational model varies as a function of Reynolds number. The thick black curve is obtained by averaging the results from the pore geometry used in the microfluidic experiment (dark gray) with four additional realizations of the “swiss-cheese” pore geometry (light gray markers) having the same porosity. The spread of the markers therefore gives a good indication for the statistical variability of the participation ratio for different pore geometries. The values of the participation ratios obtained from the two experimentally measured velocity fields are marked with yellow and blue stars, respectively, and are in good agreement with the numerical calculations.

At low Reynolds numbers, the participation ratio is practically constant at $\text{Π}\approx 0.6$, similar to a Stokes flow with fixed viscosity $\mu ={\mu }_0=24$$\text{Pa}\cdot \text{s}$. By increasing the flow rate, the local shear within the interstitial pore space also generally increases to eventually reach a point where its range of variability overlaps with the range in which the rheology of the xanthan gum solution follows a power-law behavior. At this point, the local viscosity spans over a wide range of values, leading to a drop in the participation ratio by almost 20% at a Reynolds numbers around $\text{Re}=3\times {10}^{-4}$.

Interestingly, while the absolute value of the participation ratio varies slightly from realization to realization, the location of the participation ratio minimum is determined by the fluid’s rheology and does not seem to be influenced by the details of the pore geometry. More precisely, we find for the two experimentally measured velocity fields a participation ratio of $\text{Π}=0.556$ for q_in = 0.05 µL/min and $\text{Π}=0.520$ for q_in = 5 µL/min. For the simulations, the mean participation ratios averaged over four realizations yield $\text{Π}=0.582\pm 0.002$ for q_in = 0.05 µL/min and $\text{Π}=0.525\pm 0.002$ for q_in = 5 µL/min. The participation ratios of the two simulations which share the same pore geometry as the experimental setup are $\text{Π}=0.582$ (q_in = 0.05 µL/min) and $\text{Π}=0.528$ (q_in = 5 µL/min). After reaching the minimum, where the flow is most heterogeneously distributed in the porous medium, the participation ratio increases again as the flow rate pushes the fluid’s rheology beyond the power-law regime, ultimately approaching the flow pattern of water. In this limiting case, the Newtonian behavior is recovered, but inertial effects on the flow should prevent the participation ratio to reach the same value obtained for very low Reynolds numbers, namely, $\text{Π}\approx 0.6$.

A question that naturally arises is how rheology influences the flow’s heterogeneity in space. At low and high Reynolds number, the Carreau fluid has an almost constant viscosity equal to the low and high shear limits ${\mu }_0=24$$\text{Pa}\cdot \text{s}$ and ${\mu }_{\infty }=0.001$$\text{Pa}\cdot \text{s}$, respectively. In the intermediate regime, however, the local viscosity covers a broad spectrum of values [42]. In this case, both experimental and simulation results reveal that the interplay between the disordered geometry of the pore space and the fluid rheology leads to a larger flow heterogeneity and therefore to a stronger localization pattern.

Conclusion

It is well accepted that flow and transport processes in porous media are fundamentally controlled by the complex interplay between the fluid and the pore space structure. Here we showed that these processes can be tailored by tuning the rheology of the fluid. Precisely, the heterogeneity of the pore scale structure of the medium causes a high variability of shear rates that, when matched with the nonlinear viscosity window of the non-Newtonian fluid, can substantially enhance macroscopic properties of the system like the participation ratio. These effects may be exploited to improve filters and catalysts or to enhance chemical reactions by spreading the transported chemicals more uniformly throughout the pore space. In particular, localization should have a deleterious influence on the effectiveness of catalysts subjected to flow, for example, in a packed bed chemical reaction. Precisely, the preferential channeling at the minimum of the participation number should be avoided to maximize the activity of the surface area available for reaction in the system.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

HS and JA designed the research. UE carried out the experiments with input from ES, MH, and JJ-M All authors discussed the results. HS and JA wrote the paper with input from all other authors.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We acknowledge Prof. R. Stocker and E. Burmeister from the Department of Civil, Environmental and Geomatic Engineering at ETH Zurich for kindly providing access to laboratory equipment and material. ES acknowledges SNSF PRIMA Grant 179834. JA acknowledges financial support from the Brazilian agencies CNPq, CAPES and FUNCAP. HS and JA acknowledge financial support from Petrobras (“Física do Petróleo em Meios Porosos”, Project Number: F0185. RA acknowledges financial support from the Brazilian agency CNPq and the INCT-SC Project.

References