Well-posedness and variational numerical scheme for an adaptive model in highly heterogeneous porous media

doi:10.1016/j.jcp.2022.111844

Journal of Computational Physics

Volume 475, 15 February 2023, 111844

https://doi.org/10.1016/j.jcp.2022.111844 Get rights and content

Abstract

Mathematical modeling of fluid flow in a porous medium is usually described by a continuity equation and a chosen constitutive law. The latter, depending on the problem at hand, may be a nonlinear relation between the fluid's pressure gradient and velocity. The actual shape of this relation is normally chosen at the outset of the problem, even though, in practice, the fluid may experience velocities outside of its range of applicability. We propose here an adaptive model, so that the most appropriate law is locally selected depending on the computed velocity. From the analytical point of view, we show well-posedness of the problem when the law is monotone in velocity and show existence in one space dimension otherwise. From the computational point of view, we present a new approach based on regularizing via mollification the underlying dissipation, i.e., the power lost by the fluid to the porous medium through drag. The resulting regularization is shown to converge to the original problem using Γ-convergence on the dissipation in the monotone case. This approach gives rise to a variational numerical scheme which applies to very general problems and which we validate on three test cases.

Introduction

We study the stationary flow of a Newtonian fluid in a fully saturated, highly heterogeneous porous medium. Typically, the heterogeneities come from the lithological and geometrical properties of the medium. Indeed, very different sediments (such as sandstone and carbonates) and fractures with irregular aperture may be involved. These properties impact the permeability of the domain and thus the fluid's velocity.

Our model is based on the following constitutive, or seepage, law, which is in fact a force balance: $λ (u) = - \nabla p + f,$ where u and p are respectively the seepage flux and the fluid's pressure. The term $λ (u)$ is the opposite of the drag force experienced by the fluid and the term f is a vector of external body forces, like gravity. The assumption on the operator λ most recurrent in the literature of porous media is linearity, meaning that (1.1) is Darcy's law [3], [21]. This, however, is known to be valid only for low Reynolds numbers, i.e., low fluid speeds [38], beyond which Darcy's law tends to overestimate velocities. To describe flows at higher speeds more accurately, it is common to add a quadratic term to Darcy's law to penalize high velocities and get the so-called Darcy–Forchheimer law, which is an example of a nonlinear operator λ [1], [14], [20], [25]. Other common laws are obtained by adding a higher-order term or a Laplacian term to Darcy's, yielding Forchheimer's generalized law or Brinkman's law [28], respectively.

These commonly used models are well known for being well posed and providing good predictions in a homogeneous medium. They can also be adapted to give accurate results in a heterogeneous medium by, for instance, taking spatially dependent permeabilities. However, they do not cover the case when the medium's heterogeneities yield an operator λ discontinuous in u. Since a linear law is better adapted to low Reynolds numbers whereas a nonlinear law gives a better description of high-speed regimes, one may expect that allowing λ to be linear under some given speed threshold and to be nonlinear above this threshold should deliver improved results. This consideration motivated us to study discontinuous seepage laws in [17] (which has also been recently explored numerically in [37] in the context of landfill management), and we continue here the work started therein.

To handle mathematically such a discontinuous problem, we make use of a multivalued version of (1.1) in the case when λ involves no space derivatives of the flux (thus excluding Brinkman's law). We show its well-posedness when the drag force is maximal monotone in the flux variable, using classical tools from multivalued operator theory. We also prove existence of solutions when the monotonicity fails and the space dimension d equals one. Consequently, we introduce a regularized, monovalued approximation of the multivalued problem, which can be solved numerically using classical fixed-point and finite-element methods. This regularization is based on the mollification of the dissipation (i.e., the power the fluid loses to the surrounding medium because of drag), and we show that it converges to the original problem using variational results, in particular, the Γ-convergence of the regularized dissipation to the unregularized one when it is convex; when the dissipation is nonconvex, the regularized problem is still shown to have solutions when $d = 1$ , but is not proved to converge in any case. After applying the fixed-point and finite-element methods, we compare the resulting regularized algorithm to that introduced in [17], called the transition-zone tracking algorithm. The latter is based on iteratively locating the zones separating any pair of different speed regimes and solving the appropriate law in every region thus defined; it differs from the algorithm derived in this paper, which, instead of tracking the transition zones sharply, spreads them out smoothly by solving the problem using the same regularized law in the whole medium. The two approaches give very similar results for $d = 1$ and for a combination of two different speed regimes, as we show on a simple test case, but the regularized approach offers the advantage of applying immediately for $d > 1$ and for any number of regimes, as showcased by two other test cases.

This paper is organized as follows. In Section 2, the physical and multivalued framework is introduced and motivated, and in Section 3 the weak formulation is given and the well-posedness results proved in the adequate functional spaces for general constitutive operators including no space derivatives of the flux. Section 4 contains the well-posedness theory specifically formulated for some common examples of constitutive laws. In Section 5, the regularizing approach is introduced and its convergence demonstrated, while in Section 6 we briefly describe the numerical approximation adopted to solve the regularized problem. Section 7 contains the numerical results of the three test cases. Finally, in Section 8, we give conclusions. For the reader's convenience, appendices are provided recalling basic notions on multivalued operators, functionals and mollification.

Section snippets

Physical framework

We denote by $Ω \subset R^{d}$ the porous medium, which we assume to be open, bounded and with Lipschitz boundary ∂Ω; we write n the outward normal unit vector of ∂Ω. The unknowns of the problems discussed throughout the paper are the fluid's pressure $p : Ω \to R$ and the seepage flux $u : Ω \to R^{d}$ defined by $u = ρ ϕ V$ , where ϕ is the medium's porosity and ρ and V are the fluid's density and velocity. This relation between flux and velocity justifies that the terms “flux” and “velocity” may be used interchangeably. We

Mathematical framework

For all $α \in (0, \infty)$ , $β \in [1, \infty)$ and $A \subset R^{d}$ measurable, we denote by $L^{β} (A)$ and $W^{α, β} (A)$ the Lebesgue space of measurable functions on A with integrable βth power and the αth-order Sobolev space associated with $L^{β} (A)$ ; we also write $L^{β} (A)$ for ${(L^{β} (A))}^{d}$ and use ${‖ \cdot ‖}_{β}$ for the canonical norm on $L^{β} (Ω)$ . As usual in these spaces, equality is intended in the almost everywhere sense.

Let $r \in (1, \infty)$ and write $s \in (1, \infty)$ its dual exponent, i.e., $s = r / (r - 1)$ . We fix $q \in L^{r} (Ω)$ , $f \in L^{s} (Ω)$ , $u_{0} \in L^{r} (Σ_{v})$ and $p_{0} \in W^{\frac{1}{r}, s} (Σ_{p})$ , and let $Λ : L^{r} (Ω) ⇉ L^{s}$

Well-posedness for dissipative drag operators

We pick $r \in [2, \infty)$ and let $s \in (1, 2]$ be the dual exponent of r, that is, $s = r / (r - 1)$ . As in Section 3, we fix $q \in L^{r} (Ω)$ , $f \in L^{s} (Ω)$ , $u_{0} \in L^{r} (Σ_{v})$ and $p_{0} \in W^{\frac{1}{r}, s} (Σ_{p})$ , and we let $Λ : L^{r} (Ω) ⇉ L^{s} (Ω)$ be such that, for all $u \in L^{r} (Ω)$ , there holds $Λ (u) \neq \emptyset$ .

Furthermore, we write $M_{sym}^{+}$ the set of $d \times d$ , real, symmetric, positive definite matrices. Given $M \in M_{sym}^{+}$ , we write ${‖ \cdot ‖}_{M}$ the Euclidean norm weighted by $M$ , that is, for all $x \in R^{d}$ , ${‖ x ‖}_{M} = \sqrt{M x \cdot x};$ denoting the spectrum of $M$ by $Sp (M)$ , we have $\sqrt{\min Sp (M)} ‖ x ‖ ⩽ {‖ x ‖}_{M} ⩽ \sqrt{\max Sp (M)} ‖ x ‖ .$ When $M : Ω \to M_{sym}^{+}$

Regularized problem for dissipative operators

We would like to use classical numerical schemes to solve Problem 3.2, such as Picard iterations combined with the Raviart–Thomas or the mixed virtual element methods (cf. Section 6). To this end, we propose first to approximate Problem 3.2 by a monovalued problem obtained from a convolutional regularization of the dissipation. Indeed, we restrict here to $Λ : L^{r} (Ω) ⇉ L^{s} (Ω)$ , $r ⩾ 2$ , being a dissipative operator, i.e., an operator satisfying Assumption 4.1. Also, as done in Sections 3 and 4, we fix $q \in L^{r} ($

Numerical approximation

Let us describe, for all $ε > 0$ , the numerical approximation of Problem (5.3)(ε), which is inherently nonlinear. In fact, even if the law in question is of the jump type discussed in Section 4.3.2 with $m = 0$ (i.e., the law in each speed region is linear), the resulting regularized law $Λ_{ε}$ is nonlinear in u. Inspired by the standard fixed-point algorithm, we propose the following algorithm to solve Problem (5.3)(ε): given $u^{0} \in L^{r} (Ω)$ , find $u^{n}$ such that ${\begin{matrix} div u^{n} = q, \\ ϕ_{ε} ({‖ u^{n - 1} ‖}_{D}^{2}) D u^{n} = - \nabla p^{n} + f, \end{matrix} in Ω,$ for all $n ⩾ 1$

Numerical results

In this section, we propose three test cases to validate and show the capabilities of the proposed model and of the variational numerical scheme of Sections 5 and 6. We focus on drag operators of the jump type discussed in Section 4.3.2 with $D = I$ . First, in Section 7.1, we compare it against the transition-zone tracking algorithm proposed in [17]. The second example, described in Section 7.2, is a problem where three flow regimes may coexist in the domain; we consider both linear and nonlinear

Conclusion

In this work, we have presented a mathematical framework for adaptively choosing the most appropriate constitutive law depending on the developed fluid velocity. The setup is mathematically formulated as a multivalued problem and, under the hypothesis of maximal monotonicity of the drag operator, the problem has been shown to be weakly well posed. If the drag operator fails to be monotone, we have shown existence of weak solutions when $d = 1$ . Moreover, we have derived a monovalued regularization

CRediT authorship contribution statement

All the authors have contributed equally to the work.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References (38)

E. Ahmed et al.
A multiscale flux basis for mortar mixed discretizations of reduced Darcy–Forchheimer fracture models
Comput. Methods Appl. Mech. Eng.
(2019)
N.H. Chieu
Integral of the Clarke subdifferential mapping and a generalized Newton–Leibniz formula
Nonlinear Anal., Theory Methods Appl.
(2010)
A. Fumagalli
Dual virtual element method in presence of an inclusion
Appl. Math. Lett.
(2018)
Y. Huang et al.
An overview on nonlinear porous flow in low permeability porous media
Theor. Appl. Mech. Lett.
(2013)
F.A. Morales et al.
A Darcy–Brinkman model of fractures in porous media
J. Math. Anal. Appl.
(2017)
J.E. Roberts et al.
Mixed and hybrid methods
J.D. Audu et al.
On the well-posedness of generalized Darcy–Forchheimer equation
Bound. Value Probl.
(2018)
J. Bear
Dynamics of Fluids in Porous Media
(1972)
L. Beirão da Veiga et al.
Mixed virtual element methods for general second order elliptic problems on polygonal meshes
ESAIM: M2AN
(2016)
J.M. Borwein
Maximality of sums of two maximal monotone operators in general Banach space
Proc. Am. Math. Soc.
(2007)

A. Braides

Γ-Convergence for Beginners

(2002)

F. Brezzi et al.

Basic principles of mixed virtual element methods

ESAIM: M2AN

(2014)

F.E. Browder

Nonlinear maximal monotone operators in Banach space

Math. Ann.

(1968)

M.A. Christie et al.

SPE-66599-MS, Chapter Tenth SPE Comparative Solution Project: A Comparison of Upscaling Techniques

(2001)

F.H. Clarke

Optimization and Nonsmooth Analysis

(1990)

J. Czarnovska et al.

On the Darboux property of multivalued functions

Demonstr. Math.

(1992)

F. Dassi et al.

The mixed virtual element method on curved edges in two dimensions

Comput. Methods Appl. Mech. Eng.

(2021)

N. Frih et al.

Modeling fractures as interfaces: a model for Forchheimer fractures

Comput. Geosci.

(2008)

A. Fumagalli et al.

Dual virtual element methods for discrete fracture matrix models

Oil Gas Sci. Technol. - Rev. d'IFP Energies nouvelles

(2019)

Cited by (0)

View full text

Well-posedness and variational numerical scheme for an adaptive model in highly heterogeneous porous media

Abstract

Introduction

Section snippets

Physical framework

Mathematical framework

Well-posedness for dissipative drag operators

Regularized problem for dissipative operators

Numerical approximation

Numerical results

Conclusion

CRediT authorship contribution statement

Declaration of Competing Interest

Comput. Methods Appl. Mech. Eng.

Nonlinear Anal., Theory Methods Appl.

Appl. Math. Lett.

Theor. Appl. Mech. Lett.

J. Math. Anal. Appl.

On the well-posedness of generalized Darcy–Forchheimer equation

Bound. Value Probl.

Dynamics of Fluids in Porous Media

Mixed virtual element methods for general second order elliptic problems on polygonal meshes

ESAIM: M2AN

Maximality of sums of two maximal monotone operators in general Banach space

Proc. Am. Math. Soc.

Γ-Convergence for Beginners

Basic principles of mixed virtual element methods

ESAIM: M2AN

Nonlinear maximal monotone operators in Banach space

Math. Ann.

SPE-66599-MS, Chapter Tenth SPE Comparative Solution Project: A Comparison of Upscaling Techniques

Optimization and Nonsmooth Analysis

On the Darboux property of multivalued functions

Demonstr. Math.

The mixed virtual element method on curved edges in two dimensions

Comput. Methods Appl. Mech. Eng.

Modeling fractures as interfaces: a model for Forchheimer fractures

Comput. Geosci.

Dual virtual element methods for discrete fracture matrix models

Oil Gas Sci. Technol. - Rev. d'IFP Energies nouvelles