1. Flow through porous material : the Darcy model

We consider a closed bounded domain \(\Omega\subset\mathbb R ^d\).

Whenever \(\underline u,\underline v\in[L^2(\Omega)]^d\), we denote :

\[ (\underline u,\underline v)_0 = \int_\Omega\underline u\cdot\underline v\]

If \(\varphi,\psi\in L^2(\Omega)\), we define :

\[ (\varphi,\psi)_0=\int_\Omega\varphi\psi\]

If moreover \(A\subset\Omega\) has dimension \(d-1\) or \(d\), then we define the corresponding bilinear form :

\[ \langle\varphi,\psi\rangle_A=\int_A\varphi\psi\]

1.1. Mixed Darcy problem

Find \(\underline u:\Omega\to\mathbb R^d\) and \(\varphi:\Omega\to\mathbb R\) such that :

\[ \left\{\begin{array}[c]{rll} \underline u +\underline{\underline\kappa}\nabla\varphi & =\underline f_1 & \text{in }\Omega\\ \nabla\cdot\underline u & = f_2 & \text{in }\Omega \end{array}\right.\]

where \(\underline f_1\in [L^2(\Omega)]^d\) is the source term and \(f_2\in L^2(\Omega)\) is the incompressibility defect.

Here, \(\underline{\underline\kappa}\) is a positive definite matrix (often assumed symmetric or even diagonal) representing the hydraulic conductivity whose coefficients are in \(m s^{-1}\), and \(\varphi=z+\frac p{\rho g}\) is the fluid charge or piezometric head expressed in \(m\), where \(g\) is the gravity acceleration, \(\rho\) is the fluid density and \(z\) denotes the height from a reference level.


We usually consider the case where \(\underline{\underline\kappa}\) is symmetric, or even diagonal, and \(\underline f_1\equiv0\). If moreover one wishes to consider incompressible flow, it suffices to set \(f_2\equiv0\) as well.

1.2. Primal Darcy Problem

The mixed formulation yields :

\[ \underline u = -\underline{\underline\kappa}\nabla\varphi + \underline f_1\quad\text{ in }\Omega\]

so writing \(\nabla\cdot\underline u = f_2\) gives :

\[ - \nabla\cdot\underline{\underline\kappa}\nabla\varphi + \nabla\cdot\underline f_1 = f_2\quad\text{ in }\Omega\]

In the sequel, we will focus on the case with no source term, i.e. \(\underline f_1\equiv0\). For the sake of simplicity, we will denote \(f=f_2\).

The primal form then becomes :

\[ -\nabla\cdot\underline{\underline\kappa}\nabla\varphi=f\quad\text{ in }\Omega\]

Moreover, we suppose \(\underline{\underline\kappa}\) is diagonal, and we denote \(\kappa_i,i\in\{1,...,d\}\) its diagonal entries. We emphasize the fact that they are all positive constants.

1.3. Boundary conditions

We need to impose boundary conditions on \(\partial\Omega\). As in the Stokes case, we will consider Dirichlet, Neumann and Robin type boundary conditions.

Let us decompose the boundary \(\partial\Omega\) in two relatively open disjoint parts :

\[ \partial\Omega=\overline{\Gamma_{nat}}\cup\overline{\Gamma_{ess}}\]
\[ \Gamma_{nat}\cap\Gamma_{ess}=\emptyset\]

on which we impose the conditions :

\[ \left\{\begin{array}{rll} -(\underline{\underline\kappa}\nabla\varphi)\cdot\underline n & = u_{nat} & \text{ on }\Gamma_{nat}\\ \varphi & = \varphi_{ess} & \text{ on }\Gamma_{ess} \end{array}\right.\]

Thus, the natural condition is a Neumann condition for \(\varphi\) on which depends the variational formulation, and the essential condition is a Dirichlet condition for \(\varphi\) on which depends the test function space.

1.4. Weak formulation

A nice functional space for this primal Darcy problem is the space of all the infinitely differentiable scalar functions on \(\Omega\) vanishing close to the essential boundary : \(V_0 = \left\{\psi\in C^\infty(\Omega)\quad|\quad\left.\psi\right|_{\Gamma_{ess}}=0\right\}\). Nevertheless, it is not complete. Our test function space will be its \(H^1(\Omega)\) completion.


The test space for the primal Darcy problem is the space of function whose trace vanishes on the essential boundary :

\[H_0^1(\Omega)\subset V=\overline{V_0}^{H^1(\Omega)}\subset H^1(\Omega)\]

It inherits the \(H^1(\Omega)\) norm.

Let us multiply the primal form by \(\psi\in V\). The left hand side term becomes :

\[ \begin{array}{rl} \displaystyle\int_\Omega-(\nabla\cdot\underline{\underline\kappa}\nabla\varphi)\psi=&-\displaystyle\int_\Omega\nabla\cdot(\underline{\underline\kappa}(\nabla\varphi)\psi)+\int_\Omega(\underline{\underline\kappa}\nabla\varphi)\cdot\nabla\psi\\ =&\displaystyle-\int_{\Gamma_{nat}}\psi\underline{\underline\kappa}\nabla\varphi\cdot\underline n+\underline{\underline\kappa}\int_\Omega\nabla\varphi\cdot\nabla\psi\\ =&-\langle u_{nat},\psi\rangle_{\Gamma_{nat}}+(\underline{\underline\kappa}\nabla\varphi,\nabla\psi)_0 \end{array}\]

while the right hand side term becomes \((f,\psi)_0\).

To impose the Dirichlet condition \(\varphi=\varphi_{ess}\) on \(\Gamma_{ess}\), one has to involve the trace operator, since the restriction of \(\varphi\) is not properly defined anymore. By surjectivity of the trace operator, \(\varphi_{ess}\) admits an extension (not unique) in \(H^1(\Omega)\) we still denote \(\varphi_{ess}\) for simplicity. The Dirichlet condition becomes :

\[ \varphi-\varphi_{ess}\in V\]

Let \(a:V\times V\to\mathbb R\) be a bilinear form and \(\tilde f:V\to\mathbb R\) a linear form defined by :

\[ a(\varphi,\psi)=(\underline{\underline\kappa}\nabla\varphi,\nabla\psi)_0\quad\text{ ; }\quad\tilde f(\psi)=(f,\psi)_0-\langle u_{nat},\psi\rangle_{\Gamma_{nat}}\]

The weak formulation of the Darcy problem reads as follows :

\[ \left\{\begin{array}{ll} \text{find }\varphi\in H^1(\Omega)\text{ such that }\varphi-\varphi_{ess}\in V\text{ and}&\\ a(\varphi,\psi)=\tilde f(\psi)&\\ \end{array}\right.\]

The Darcy problem in weak form admits an unique solution \(\varphi\).


This comes straightforwardly from Lax-Milgram theorem. It is sufficient to show that \(a\) is bounded and coercive, and that \(\tilde f\) is bounded.

  • Thanks to Cauchy-Schwarz inequality, the boundedness of \(a\) comes easily :

\[ |a(\varphi,\psi)|=|(\underline{\underline\kappa}\nabla\varphi,\nabla\psi)_0|\leqslant\max_{i\in\{1,...,d\}}\{\kappa_i\}|(\nabla\varphi,\nabla\psi)_0|\leqslant\max_{i\in\{1,...,d\}}\{\kappa_i\}\|\varphi\|_V\|\psi\|_V\]
  • The coercivity of \(a\) is also easy to prove :

\[ a(\varphi,\varphi)=(\underline{\underline\kappa}\nabla\varphi,\nabla\varphi)_0\geqslant\min_{i\in\{1,...,d\}}\{\kappa_i\}(\varphi,\varphi)_V=\min_{i\in\{1,...,d\}}\{\kappa_i\}\|\varphi\|_V^2\]
  • Using the Cauchy-Schwarz inequality and the Poincaré inequality, we get the boundedness of \(\tilde f\). We first apply the Riesz representation theorem to find \(\tilde F\in L^2(\Omega)\) representing \(\tilde f\).

\[ |\tilde f(\psi)|=|(\tilde F,\psi)_0|\leqslant\|\tilde F\|_0\|\psi\|_0\leqslant C\|\tilde F\|_0\|\psi\|_V\]

1.5. Numerical simulation