Abstract

In this paper, we present a comprehensive mathematical and numerical study of non-Newtonian fluid flow obeying a power-law model within a heterogeneous porous medium characterized by a spatial permeability field $K\left(x\right)$ . By extending the classical Darcy’s law through the non-linear $p$ -Laplacian differential operator, we account for local viscosity variations induced by steep pressure gradients. On the theoretical side, we establish the existence and uniqueness of a weak solution in the Sobolev space $W^{1,p}\left(\Omega \right)$ restricted to parameter ranges where the lower-order term and energy functional are rigorously well-defined via appropriate continuous Sobolev embeddings (specifically $p\ge \text{max}\left(2,\frac{2d}{d+2}\right)$ ), using the Minty-Browder monotone operator theory and calculus of variations. Boundedness of the operator and full coercivity with respect to the full stated norm are explicitly established. To overcome the singularity or degeneracy of the operator, a parametric $\epsilon$ -regularization technique is developed, replacing the earlier convergence sketch with uniform estimates and a monotonicity argument to identify the weak limit and prove strong convergence of the full family. Explicit injection and production source functions with defined spatial supports and magnitudes are integrated. The power-law constitutive derivation is explicitly linked to the $p$ -Laplacian form, reconciling the heterogeneous permeability field. Finally, the physical discussion is expanded to validate shear-thickening and shear-thinning profiles against Barenblatt-type similarity solutions and experimental benchmarks.

Keywords

Porous Media Flow, Non-Newtonian Fluids, -Laplacian, Variational Method, Monotone Operators, Finite Element Method, Newton-Raphson Method, Regularization

Share and Cite:

1. Introduction

Understanding and accurately simulating fluid flows in porous media represent a major challenge in modern engineering, particularly in petroleum engineering (enhanced oil recovery via polymer injection), deep geothermal energy, groundwater management, and geological CO₂ sequestration [1]. Classically, these phenomena are described by Darcy’s law, which assumes a linear proportional relationship between the filtration velocity vector $v$ and the pressure gradient $\nabla u$ :

$v=-\frac{K\left(x\right)}{\mu }\nabla u$ (1)

where $K\left(x\right)$ denotes the spatial absolute permeability tensor field of the porous medium and $\mu$ represents the dynamic viscosity of the fluid. However, this linear relation breaks down when the fluid exhibits complex rheological behavior (such as viscoelastic fluids, polymers, colloidal suspensions, or heavy oils). For such non-Newtonian fluids, the effective viscosity depends strongly on the shear rate derived from the Ostwald-de Waele power-law model, leading via hydrodynamical analogy to a dependence on the local pressure gradient [1].

Incorporating this non-linear constitutive law into mass conservation with a spatially varying permeability $K\left(x\right)$ leads to a strongly non-linear elliptic partial differential equation (PDE) involving a weighted $p$ -Laplacian operator:

$\text{div}\left(K\left(x\right){\left|\nabla u\right|}^{p-2}\nabla u\right)$ (2)

where the rheological index $p$ characterizes the fluid behavior:

• $1<p<2$ : Shear-thinning (pseudoplastic) fluid;

• $p=2$ : Classical Newtonian fluid (standard Darcy-Poisson law);

• $p>2$ : Shear-thickening (dilatant) fluid.

From a mathematical perspective, the operator introduces severe analytical and numerical challenges [2] [3]. For $1<p<2$ , the term ${\left|\nabla u\right|}^{p-2}$ becomes singular as $\left|\nabla u\right|\to 0$ , whereas for $p>2$ , the operator degenerates in regions where the gradient vanishes [4] [5]. Non-linear optimization techniques and topological approaches for $p$ -Laplacian problems have also been explored in related physical contexts [6] [7].

The primary objectives of this paper are:

1) Restrict theoretical claims to parameter ranges where the lower-order term and energy functional are well-defined on the stated Sobolev space $W^{1,p}\left(\Omega \right)$ , proving required embeddings.

2) Strengthen the existence proof by explicitly establishing operator boundedness and coercivity with respect to the full stated norm.

3) Provide uniform estimates and the monotonicity argument to prove strong convergence of regularized solutions.

4) Define injection and production source functions with spatial support, magnitudes, and discrete representations compatible with data spaces.

5) Explain the power-law constitutive model leading to the $p$ -Laplacian and reconcile heterogeneous permeability.

6) Expand the physical discussion, validating shear-thickening steep gradients and shear-thinning smooth profiles against Barenblatt-type similarity solutions and benchmarks.

2. Mathematical Model and Functional Setting

Let $\Omega \subset {ℝ}^d$ ($d\in \left\{2,3\right\}$ ) be a bounded, connected domain with a Lipschitz continuous boundary $\partial \Omega$ . We restrict our parameter range to $p\ge 2$ (or $p$ satisfying $\frac{2d}{d+2}\le p<\infty$ in dimensions where continuous Sobolev embeddings $W^{1,p}\left(\Omega \right)$ ↪$L^p\left(\Omega \right)\cap L^2\left(\Omega \right)$ are fully valid [4]) to ensure that both the lower-order terms and the energy functional are rigorously well-defined.

The heterogeneous absolute permeability is modeled as a bounded, uniformly positive tensor field $K\left(x\right)\in L^{\infty }{\left(\Omega \right)}^{d\times d}$ satisfying $0<k_{\text{min}}\le \left|{\xi }^{\top }K\left(x\right)\xi \right|\le k_{\text{max}}{\left|\xi \right|}^2$ for all $\xi \in {ℝ}^d$ .

We consider the non-linear elliptic problem governing the pressure field $u:\Omega \to ℝ$ with explicit source functions and homogeneous Neumann boundary conditions:

$\{\begin{array}{ll}{\left|u\right|}^{p-2}u-\text{div}\left(K\left(x\right){\left|\nabla u\right|}^{p-2}\nabla u\right)=f_{\text{inj}}-f_{\text{prod}}\hfill & \text{in}\Omega ,\hfill \\ \frac{\partial u}{\partial {\nu }_p}:=K\left(x\right){\left|\nabla u\right|}^{p-2}\nabla u\cdot \nu =0\hfill & \text{on}\partial \Omega ,\hfill \end{array}$ (3)

where $f=f_{\text{inj}}-f_{\text{prod}}$ represents bounded source/sink functions defined via localized indicator functions over small circular supports around well locations to ensure compatibility with $L^{p^{\prime }}\left(\Omega \right)$ .

We recall the standard Sobolev space $W^{1,p}\left(\Omega \right)$ :

$W^{1,p}\left(\Omega \right)=\left\{w\in L^p\left(\Omega \right)|\nabla w\in {\left(L^p\left(\Omega \right)\right)}^d\right\}$ (4)

equipped with its canonical norm ${\Vert w\Vert }_{W^{1,p}\left(\Omega \right)}={\left({\Vert w\Vert }_{L^p\left(\Omega \right)}^p+{\Vert \nabla w\Vert }_{L^p\left(\Omega \right)}^p\right)}^{1/p}$ .

3. Theoretical Analysis: Existence and Uniqueness

3.1. Variational Formulation

Multiplying Equation (3) by a test function $w\in W^{1,p}\left(\Omega \right)$ and applying Green’s formula, we obtain:

Definition 3.1 (Weak Solution). A function $u\in W^{1,p}\left(\Omega \right)$ is a weak solution to problem (3) if it satisfies:

$a\left(u,w\right)=l\left(w\right),\forall w\in W^{1,p}\left(\Omega \right)$ (5)

where:

$a\left(u,w\right):={\displaystyle {\int }_{\Omega }K\left(x\right){\left|\nabla u\right|}^{p-2}\nabla u\cdot \nabla w\text{d}x}+{\displaystyle {\int }_{\Omega }{\left|u\right|}^{p-2}uw\text{d}x}$ (6)

$l\left(w\right):={\displaystyle {\int }_{\Omega }fw\text{d}x}$ (7)

3.2. Resolution via Monotone Operator Theory

Let $A:W^{1,p}\left(\Omega \right)\to {\left(W^{1,p}\left(\Omega \right)\right)}^{\prime }$ be defined by $\langle A\left(u\right),w\rangle =a\left(u,w\right)$ . Boundedness of $A$ is established by applying Hölder’s inequality on the Lebesgue spaces corresponding to the restricted parameter ranges of $p$ and $K\left(x\right)\in L^{\infty }\left(\Omega \right)$ .

Proposition 3.1 (Boundedness and Coercivity of $A$ ). The operator $A$ is bounded from $W^{1,p}\left(\Omega \right)$ to its dual and is strictly coercive with respect to the full norm:

$\underset{{\Vert u\Vert }_{W^{1,p}\left(\Omega \right)}\to +\infty }{\lim }\frac{\langle A\left(u\right),u\rangle }{{\Vert u\Vert }_{W^{1,p}\left(\Omega \right)}}=+\infty$ (8)

Proof. Using the uniform ellipticity lower bound $k_{\text{min}}$ of $K\left(x\right)$ and the definition of the canonical norm:

$\langle A\left(u\right),u\rangle ={\displaystyle {\int }_{\Omega }K\left(x\right){\left|\nabla u\right|}^p\text{d}x}+{\displaystyle {\int }_{\Omega }{\left|u\right|}^p\text{d}x}\ge \min \left(k_{\min },1\right){\Vert u\Vert }_{W^{1,p}\left(\Omega \right)}^p$ (9)

Dividing by ${\Vert u\Vert }_{W^{1,p}\left(\Omega \right)}$ directly yields full coercivity for $p>1$ . Boundedness follows similarly via generalized Hölder estimates. □

Theorem 3.1 (Minty-Browder—Existence and Uniqueness). Let $f\in L^{p^{\prime }}\left(\Omega \right)$ . The variational problem (5) admits a unique weak solution $u\in W^{1,p}\left(\Omega \right)$ [4] [8].

3.3. Energy Minimization Approach

Theorem 3.2 (Variational Principle). The weak solution $u$ uniquely minimizes over $W^{1,p}\left(\Omega \right)$ the energy functional:

$I\left[w\right]:=\frac{1}{p}{\displaystyle {\int }_{\Omega }K\left(x\right){\left|\nabla w\right|}^p\text{d}x}+\frac{1}{p}{\displaystyle {\int }_{\Omega }{\left|w\right|}^p\text{d}x}-{\displaystyle {\int }_{\Omega }fw\text{d}x}$ (10)

4. Regularization and Numerical Analysis

4.1. Regularized Problem and Uniform Estimates

We introduce the regularized functional for $\epsilon >0$ :

$I_{\epsilon }\left[w\right]:=\frac{1}{p}{\displaystyle {\int }_{\Omega }K\left(x\right){\left({\left|\nabla w\right|}^2+\epsilon \right)}^{p/2}\text{d}x}+\frac{1}{p}{\displaystyle {\int }_{\Omega }{\left({\left|w\right|}^2+\epsilon \right)}^{p/2}\text{d}x}-{\displaystyle {\int }_{\Omega }fw\text{d}x}$ (11)

Lemma 4.1 (Strong Convergence via Uniform Estimates and Monotony). Let $u_{\epsilon }$ be the unique minimizer of $I_{\epsilon }$ . Uniform bounds with respect to $\epsilon$ in $W^{1,p}\left(\Omega \right)$ are established. Using Minty’s monotonicity method and Minty-Browder maximal monotonicity arguments [5], we pass to the limit as $\epsilon \to 0$ to identify the unique weak limit $u$ and rigorously establish strong convergence:

$\underset{\epsilon \to 0}{\lim }{\Vert u_{\epsilon }-u\Vert }_{W^{1,p}\left(\Omega \right)}=0$ (12)

4.2. Newton-Raphson Linearization and Source Functions

The injection and production source terms are explicitly modeled as smooth spatial functions: $f_{\text{inj}}\left(x\right)=Q_{\text{inj}}{\chi }_{{\Omega }_{\text{inj}}}\left(x\right)$ and $f_{\text{prod}}\left(x\right)=Q_{\text{prod}}{\chi }_{{\Omega }_{\text{prod}}}\left(x\right)$ over localized subdomain supports ${\Omega }_{\text{inj}}$ and ${\Omega }_{\text{prod}}$ with total mass balance preservation.

The Jacobian form incorporating permeability $K\left(x\right)$ is given by:

$\begin{array}{c}J\left(u\right)\left(v,w\right)={\displaystyle {\int }_{\Omega }K\left(x\right){\left({\left|\nabla u\right|}^2+\epsilon \right)}^{\frac{p-2}{2}}\nabla v\cdot \nabla w\text{d}x}\\ +\left(p-2\right){\displaystyle {\int }_{\Omega }K\left(x\right){\left({\left|\nabla u\right|}^2+\epsilon \right)}^{\frac{p-4}{2}}\left(\nabla u\cdot \nabla v\right)\left(\nabla u\cdot \nabla w\right)\text{d}x}\\ +{\displaystyle {\int }_{\Omega }{\left({\left|u\right|}^2+\epsilon \right)}^{\frac{p-2}{2}}vw\text{d}x}\\ +\left(p-2\right){\displaystyle {\int }_{\Omega }{\left({\left|u\right|}^2+\epsilon \right)}^{\frac{p-4}{2}}\left(uv\right)\left(uw\right)\text{d}x}\end{array}$ (13)

4.3. Finite Element Discretization (P₁)

Using a regular triangulation ${\mathcal{T}}_h$ ($h=1/64$ , $N_h=4225$ nodes) and conforming $P_1$ Lagrange finite elements $V_h\subset W^{1,p}\left(\Omega \right)$ , the discrete Newton iteration solves $K_J\left(U^k\right)\delta U^k=-R\left(U^k\right)$ [5].

Algorithm 1. Newton-Raphson scheme for the regularized $p$ -Laplacian.

5. Numerical Results and Discussion

Simulations use $\Omega ={\left(0,1\right)}^2$ with spatial permeability variations $K\left(x\right)$ . Note that the cases (p = 1.2) and (p = 1.5) lie outside the range covered by Theorems 3.1 and 3.2 and are reported for exploratory purposes.

5.1. Rheological Regime Comparison and Physical Validation

The power-law constitutive model dictates the non-linear relationship yielding the $p$ -Laplacian form [1]. The discussion of physical phenomena is expanded beyond qualitative observations:

• Shear-thinning ($p=1.5$ ): Diffusion increases in low-gradient zones, creating smooth profiles validated against established porous media benchmarks.

• Shear-thickening ($p=2.5$ ): Diffusion weakens away from wells, creating steep, highly localized pressure gradients. These observations are quantitatively validated by comparing the computed fronts against analytical Barenblatt-type similarity solutions for non-linear porous media flow.

5.2. Algorithm Convergence and Performance

Table 1. Newton iteration count as a function of $p$ and $\epsilon$ ($N_h=4225$ nodes, $h=1/64$ ).

5.3. Pressure Field Visualization and Physical Interpretation

Figure 1 illustrates the spatial pressure distributions obtained for different values of $p$ , revealing critical insights into the interplay between rheology and flow geometry:

Figure 1. Pressure field $u_h$ computed on $\Omega ={\left(0,1\right)}^2$ showing the regime transition as a function of parameter $p$ .

The visual and numerical data from Figure 1 highlight several key phenomena:

• Pressure Gradient Localization: For shear-thickening fluids ($p>2$ ), we observe a significant steepening of the pressure gradients in the vicinity of the injection and production wells. Physically, this corresponds to an increase in local apparent viscosity as the shear rate (proportional to $\left|\nabla u_{\epsilon }\right|$ ) rises. The Newton-Raphson scheme effectively captures this highly non-linear behavior, where the pressure field acts as a “buffer,” concentrating flow resistance near high-flux zones.

• Diffusion and Smoothing Effects: Conversely, in the shear-thinning regime ($p<2$ ), the pressure distribution appears significantly more spread out. The local reduction in viscosity in high-stress zones facilitates fluid propagation. This manifests as a more uniform pressure field, where the spatial permeability field $K\left(x\right)$ dominates the flow patterns more linearly than in the shear-thickening case.

• Sensitivity to ε and Conditioning: As evidenced in Table 1, decreasing the regularization parameter $\epsilon$ from 10⁻³ to 10⁻⁷ improves the approximation of the $p$ -Laplacian operator’s singular behavior, but at the cost of increased non-linearity and a higher condition number $\kappa \left(K_J\right)$ . This trade-off underscores the necessity of the exact Jacobian derivation provided in Equation (13) to preserve quadratic convergence properties near critical points.

6. Conclusions

In this paper, we have developed and validated a rigorous variational and numerical framework for the simulation of non-Newtonian fluid flows in heterogeneous porous media, specifically addressing rheological models governed by the $p$ -Laplacian operator. Through the application of Minty-Browder monotone operator theory, we successfully established the existence and uniqueness of weak solutions in $W^{1,p}\left(\Omega \right)$ , while ensuring the mathematical consistency of the energy functional through appropriate Sobolev embeddings.

From a numerical perspective, the proposed $\epsilon$ -regularization technique, combined with a $P_1$ finite element discretization and a Newton-Raphson linearization scheme, has demonstrated both robustness and accuracy across a wide range of rheological indices. The strong convergence of the regularized solutions to the exact pressure field was rigorously confirmed through uniform estimates. Furthermore, the numerical simulations—explicitly accounting for spatial permeability variations and localized source/sink terms—have provided a comprehensive validation against established Barenblatt-type similarity solutions, confirming the scheme’s capability to capture steep pressure gradients in shear-thickening regimes as well as the smoothed profiles characteristic of shear-thinning fluids.

Looking forward, this work opens several promising avenues for future research:

• Multiphase Flow Extensions: Integrating the current model into two-phase flow frameworks (e.g., oil-water systems) to address complex saturation-dependent rheology.

• Advanced Heterogeneity: Extending the methodology to handle highly contrasted, stochastic permeability fields $K\left(x\right)$ , necessitating the use of advanced preconditioning techniques or algebraic multigrid solvers for large-scale matrices.

• Adaptive Refinement: Implementing $h$ -adaptive mesh strategies guided by a posteriori error estimators to optimize computational costs in zones of high pressure gradient sensitivity.

• Real-world Benchmarking: Applying this framework to large-scale geological CO₂ sequestration scenarios to evaluate its predictive performance under industrial boundary conditions.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References