Homogenization theory and nonlinearities in Darcy's law

Abstract

The theory of homogenization of differential equations has become an open field of research in several areas of the exact sciences and has proved to be a powerful tool for understanding the global behavior of heterogeneous materials. Despite knowing that the deduction of Darcy's law through the Navier-Stokes equations has been debated for decades, many questions remain open, mainly regarding more complex boundary conditions, cases involving multiphase flows and the numerical homogenization techniques. It is known today that Darcy's law is presented in the form of a linear relationship only for a range of hydraulic gradient and that this range overlaps the range of laminar flow through soil voids. Therefore, it is proposed in this work to understand the loss of linearity in Darcy's law, based on the theory of homogenization, modifying and exploring the limit results obtained by Allaire in 1991.

Keywords:
Homogenization theory; Navier-Stokes equations; Flow in porous media; Darcy’s law Numerical analysis

1. Introduction

After Darcy, mainly at the beginning of the 20th century, researchers such as Forchheimer and Kochina brought to light the emergence of a non-linear relationship between the hydraulic gradient and the flow velocity, a concept to be formally defined in this work. What should be emphasized here is that Darcy's law loses its linear character from a certain hydraulic gradient, a quantity that is directly proportional to the pressure gradient of the fluid (Bear, 1985). By admitting that the non-linearity is relative to a critical hydraulic gradient and that the pressure gradient is dependent on the domain where it is being evaluated, and knowing that the flow does not occur through the solid particles, but through the pores of the material, the hypothesis of a study on the pore scale is valid. Therefore, a process of homogenization of the differential equations will be presented so that the mathematical objects living in the pore scale are captured on a larger scale, in real geotechnical problems.

For the arguments just presented, a mathematical analysis of the possible mechanisms and triggers of this non-linearity is suggested, which will later be associated with a loss of hydrodynamic stability and a possible transition to the turbulent regime. To mathematically capture these phenomena, a rigorous study of the differential equations that govern the flow problem will be carried out throughout this work. Not only differential equations will be studied, but also the construction of numerical methods capable of exploring the behavior of the fluid in a tortuous and locally discontinuous domain such as the soil. This work will use arguments from functional analysis, general topology, dynamical systems and the modern theory of differential equations to support the simulations carried out by the Finite Element Method.

2. Proof of the convergence of the Stokes problem for unit cells with non-Dirichlet boundary conditions

Let the Stokes problem be studied by Allaire (1991a), being this:

Find $\left(u_{\epsilon },p_{\epsilon }\right)\in {\left[H_0^1\left({\Omega }_{\epsilon }\right)\right]}^N\times L^2\left({\Omega }_{\epsilon }\right)$ solutions of the problem (P1)

$$\nabla p_{\epsilon }-\Delta u_{\epsilon }=f$$, at Ω_ε

$$\begin{array}{c}\nabla \cdot u_{\epsilon }=0\end{array}$$, at Ω_ε (P1)

With $\epsilon \to 0$, $p_{\epsilon }$, $u_{\epsilon }$ and $f$ being the pressure, velocity and the driving forces acting on the fluid. For the problem with non-Dirichlet conditions (slip condition) on part of the boundary and with Dirichlet conditions on the rest, the functional space $H_0^1\left({\Omega }_{\epsilon }\right)$ is no longer adequate. Solutions must be sought in $H^1\left({\Omega }_{\epsilon }\right)\times L^2\left({\Omega }_{\epsilon }\right).$ Let be the variational formulation of the same problem (P1):

$\begin{array}{l}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla p_{\epsilon }\cdot vdV-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\Delta u_{\epsilon }\cdot vdV=\\ \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV\forall v\in H^1{\left({\Omega }_{\epsilon }\right)}^N\end{array}$, at Ω_ε

$\begin{array}{c}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}q\left(\nabla \cdot u_{\epsilon }\right)dV=0\forall q\in L^2\left({\Omega }_{\epsilon }\right)\end{array}$, at Ω_ε (P2)

Take $v=\varphi w_{\epsilon }^k$ and $q=\varphi q_{\epsilon }^k$ where $\varphi \in D\left(\Omega \right)$. The functions $w_{\epsilon }^k\in H^1{\left({\Omega }_{\epsilon }\right)}^N$ and $q\in L^2\left({\Omega }_{\epsilon }\right)$ are sequences whose convergence will be proved from some hypotheses. Allaire (1991a) assumed the problem in ${ℝ}^N$, with no boundary conditions except on the surface of grains (holes) periodically distributed throughout the domain. In this work, we will assume that the problem is restricted to an open set $\Omega$ such that:

$\begin{array}{c}\Omega ={\Omega }_{\epsilon }+\underset{i=1}{\overset{N\left(\epsilon \right)}{{\displaystyle \cup }}}{T}_i^{\epsilon }\end{array}$ (Property 1)

$\begin{array}{c}\partial \Omega =\underset{j=1}{\overset{n}{{\displaystyle \cup }}}{\Gamma }_j,{\Gamma }_i\cap {\Gamma }_j=\varnothing fori\ne j\end{array}$ (Property 2)

$\begin{array}{c}\partial \Omega \cap \partial T_i^{\epsilon }=\varnothing \end{array}$ (Property 3)

$\begin{array}{c}\overline{\Omega }=\Omega +\partial \Omega \end{array}$ (Property 4)

Integrating by parts (P2):

$-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\left(\nabla \cdot v\right)dV+\underset{\partial {\Omega }_{\epsilon }}{{\displaystyle \int }}(p_{\epsilon }v)\cdot ndS+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla u_{\epsilon }\cdot \nabla vdV$

$=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV\forall v\in H^1{\left({\Omega }_{\epsilon }\right)}^N$, at Ω_ε

$\begin{array}{c}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}q\left(\nabla \cdot u_{\epsilon }\right)dV=0\forall q\in L^2\left({\Omega }_{\epsilon }\right)\end{array}$, at Ω_ε (P2)

Note that the boundary condition was omitted in (P2) because the expression refers only to the functional inside the domain and not on its boundary. Replacing $v$ and $q$ by their previously given definitions and expanding (P2):

$\begin{array}{l}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\left(\varphi \nabla \cdot w_{\epsilon }^k+w_{\epsilon }^k\cdot \nabla \varphi \right)dV+\\ \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla u_{\epsilon }\cdot \left(\varphi \nabla w_{\epsilon }^k+w_{\epsilon }^k\nabla \varphi \right)dV\end{array}$

$=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot \varphi w_{\epsilon }^{k}dV\forall v\in H^1{\left({\Omega }_{\epsilon }\right)}^N$, at Ω_ε

$\begin{array}{c}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\varphi q_{\epsilon }^k\left(\nabla \cdot u_{\epsilon }\right)dV=0\forall q\in L^2\left({\Omega }_{\epsilon }\right)\backslash R\end{array}$, at Ω_ε (P2)

Consider the following six hypotheses regarding the functions $w_{\epsilon }^k$ and $q_{\epsilon }^k$, based on the work of Allaire (1991a):

$\begin{array}{c}{w}_{\epsilon }^k\in H^1{\left(\Omega \right)}^N,q_{\epsilon }^k\in L^2\left(\Omega \right)\end{array}$ (Hypothesis 1)

$\begin{array}{c}\nabla \cdot w_{\epsilon }^k=0at\Omega and{w}_{\epsilon }^{k}inside{T}_i^{\epsilon }\end{array}$ (Hypothesis 2)

$\begin{array}{c}\begin{array}{l}{w}_{\epsilon }^k\rightharpoonup e_{k}weaklyin{H}^1\left(\Omega \right)and{q}_{\epsilon }^k\rightharpoonup \\ 0weaklyin{L}^2\left(\Omega \right)\end{array}\end{array}$ (Hypothesis 3)

$\begin{array}{c}\exists {\mu }_k\in W^{-\mathrm{1,}\infty }{\left(\Omega \right)}^N\end{array}$ (Hypothesis 4)

For every sequence $v_{\epsilon }$ and for every $v$ such that:

$\begin{array}{c}\begin{array}{l}{v}_{\epsilon }\rightharpoonup vweaklyin{H}^1{\left(\Omega \right)}^N,\\ v_{\epsilon }=0at{T}_i^{\epsilon }\end{array}\end{array}$ (Hypothesis 5.a)

And for each $\varphi \in D\left(\Omega \right)$, holds:

\[ \[ (\nabla q_{\epsilon}^{k} - \Delta w_{\epsilon}^{k}, \phi v_{\epsilon})_{H^{-1}, H^{1}(\Omega)^{N}} \rightarrow \] \[ (\mu_{k}, \phi v)_{H^{-1}, H^{1}(\Omega)^{N}} \] \] (Hypothesis 5.b)

There exists and extension linear operator $R_{\epsilon }$ such that:

$\begin{array}{c}{R}_{\epsilon }\in L\left(H^1{\left(\Omega \right)}^N;H^1{\left({\Omega }_{\epsilon }\right)}^N\right)\end{array}$ (Hypothesis 6.a)

$\begin{array}{c}u\in H^1{\left({\Omega }_{\epsilon }\right)}^N\Rightarrow R_{\epsilon }\tilde{u}=u\end{array}$, at Ω_ε (Hypothesis 6.b)

$\begin{array}{c}\nabla \cdot u=0at\Omega \Rightarrow \nabla \cdot \left(R_{\epsilon }u\right)=0\text{at}{\Omega }_{\epsilon }\end{array}$ (Hypothesis 6.c)

$\begin{array}{c}{R}_{\epsilon }{u}_{H^1\left({\Omega }_{\epsilon }\right)}\le C{u}_{H^1\left(\Omega \right)},C\text{independent of}\epsilon \end{array}$ (Hypothesis 6.d)

The functions $w_{\epsilon }^k$ and $q_{\epsilon }^k$, $1\le k\le N$ are test functions for the homogenization process by the energy method. In Allaire's work (1991a) all elements of the Hilbert spaces cited in the hypotheses are functions with compact support in $\Omega$. The aim of this section is to prove that convergence also occurs for $H^1$ spaces with functions that guarantee specific boundary conditions that are not of the Dirichlet type. Rewriting (P2) with respect to the six hypotheses presented, we have:

$\begin{array}{l}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }{w}_{\epsilon }^k\cdot \nabla \varphi dV+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla u_{\epsilon }\cdot w_{\epsilon }^k\nabla \varphi dV+\\ \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\varphi \nabla u_{\epsilon }\cdot \nabla w_{\epsilon }^{k}dV\end{array}$

$=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot \varphi w_{\epsilon }^{k}dV,\forall v\in H^1{\left({\Omega }_{\epsilon }\right)}^N$, at Ω_ε

$\begin{array}{c}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\varphi q_{\epsilon }^k\left(\nabla \cdot u_{\epsilon }\right)dV=\mathrm{0,}\forall q\in L^2\left({\Omega }_{\epsilon }\right)\end{array}$, at Ω_ε (P2)

The case studied here has four boundaries where two of them receive Dirichlet conditions (for pressure) and the other two receive slip conditions (free shear stress). Figures 1 to 3 illustrate the problem to be passed to the limit.

Figure 1
Periodically perforated domain.

Figure 2
Unit cell.

Figure 3
Homogenized domain.

Let be the following boundary conditions (BC) for ${\Omega }_{\epsilon }$:

$\begin{array}{c}{\Gamma }_1\Rightarrow p_{\epsilon }=p_{\epsilon }^{in}\end{array}$ (BC1)

$\begin{array}{c}{\Gamma }_2\Rightarrow -n\cdot \left(-p_{\epsilon }I+\mu \left(\nabla u_{\epsilon }+\nabla u_{\epsilon }^T\right)\right)\cdot t=0\end{array}$ (BC2)

$\begin{array}{c}{\Gamma }_3\Rightarrow p_{\epsilon }=p_{\epsilon }^{out}\end{array}$ (BC3)

$\begin{array}{c}{\Gamma }_4={\Gamma }_2\end{array}$ (BC4)

$\begin{array}{c}{\Gamma }_5\Rightarrow u_{\epsilon }=0\end{array}$ (BC5)

How then does the functional behave for each boundary? Intuitively, different functionals are expected to act on each type of boundary. It remains to be seen whether they are compatible.

On the boundary ${\Gamma }_1$ one can write the inlet pressure from the Bernouilli relation such that:

$\begin{array}{c}{p}_{\epsilon }{}_0=p_{\epsilon }+\frac{1}{2}\rho {\left|u_{\epsilon }\right|}^2\end{array}$ (1)

Where $p_{\epsilon }$ is the stagnation pressure imposed on the boundary and p is the static pressure on the boundary. It is important to point out at this point that there are no boundary conditions for pressures in the Stokes and Navier-Stokes equations; for this, the Bernoulli’s relation is assumed. Assuming $\rho =1$, the pressure boundary condition becomes a Dirichlet condition for velocities, such that the imposed velocity is:

$\begin{array}{c}|u_{\epsilon }{|}_{{\Gamma }_1}=\sqrt{2\left(p_{\epsilon }{}_0-p_{\epsilon }\right)}\end{array}$ (2. a)

Assuming that this input velocity is purely horizontal, it becomes possible to construct the velocity vector from its magnitude such that:

$\begin{array}{c}{u}_{\epsilon }{}_{{\Gamma }_1}=\left[\sqrt{2\left(p_{\epsilon }{}_0-p_{\epsilon }\right)}\right]e_k,e_k=\left(\begin{array}{c}1\\ 0\end{array}\right)\end{array}$ (2. b)

Note that for the equality in (2.b) to be true, the right-hand side must also be in $H^1$. It is known that $p_{\epsilon }\in L^2\left({\Omega }_{\epsilon }\right)$ and that, being $p_{\epsilon }{}_0\in L^2\left({\Omega }_{\epsilon }\right)$ then, by the linearity of space, $\left(p_{\epsilon }{}_0-p_{\epsilon }\right)=\widehat{p_{\epsilon }}\in L^2\left({\Omega }_{\epsilon }\right),$ therefore, every sequence or subsequence $\widehat{p_{\epsilon }}$ is continuous in ${\Omega }_{\epsilon }$.

On the boundaries ${\Gamma }_2$ and ${\Gamma }_4$, the free sliding condition operates, where the shear stress is considered zero. This condition comes from Navier's hypothesis about the proportionality between sliding velocity and shear stress. One has:

$\begin{array}{c}{n}^T\cdot \left(2\mu \mathbb{D}\left(u\right)-pI\right)\cdot t-\beta u\cdot t=0\end{array}$ (3. a)

where

$\begin{array}{c}D\left(u\right)=\mu \left(\nabla u+\nabla u^T\right)\end{array}$ (3. b)

Where β is the coefficient of friction between the fluid and the boundary, so $\beta =0$ indicates free sliding. It is important to point out that the imposition of a slip boundary condition affects the very formulation of the variational problem, as the condition acts on the stress tensor itself, thus necessitating another variational formulation for this specific boundary. Allaire (1991b), in another work on homogenization with slip conditions at the boundary of the obstacle (hole), proved the coercivity and convergence of the Stokes system to a Brinkman-type law. Here, the free-sliding condition is imposed on the domain boundary and not on the grain surface (which is maintained with Dirichlet), but the result applies equally. It can be proposed that, by $\nabla \cdot u_{\epsilon }=0$, we have:

$\begin{array}{c}\Delta u_{\epsilon }=2\nabla \cdot D\left(u\right)\end{array}$ (4)

The term corresponding to the Laplacian in the variational becomes:

$\begin{array}{c}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}2\nabla \cdot \mathbb{D}\left(u\right)\cdot vdV\end{array}$ (5)

whose integration by parts results in:

$\begin{array}{c}-2\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u\right)\cdot \nabla vdV+2\underset{{\Gamma }_{\mathrm{2,4}}}{{\displaystyle \int }}\mathbb{D}\left(u\right)\cdot v\cdot ndS\end{array}$ (6)

The first term, by symmetry, can be rewritten as:

$\begin{array}{c}\begin{array}{l}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u\right)\cdot \nabla vdV=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u\right)\cdot \frac{\nabla v}{2}dV+\\ \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u\right)\cdot \frac{\nabla v^T}{2}dV=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u\right)\cdot \mathbb{D}\left(v\right)dV\end{array}\end{array}$ (7)

Because, by symmetry, $\mathbb{D}{\left(u\right)}^T=\mathbb{D}\left(u\right)$. Knowing that the Stokes equation can be rewritten with the explicit tensor:

$\begin{array}{c}\nabla \cdot \left(pI-\mu \mathbb{D}\left(u\right)\right)=f\end{array}$ (8)

The variational formulation can then be rewritten as:

$\begin{array}{c}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla \cdot \left(-pI+\mu \mathbb{D}\left(u\right)\right)\cdot vdV=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV\end{array}$ (9. a)

$\begin{array}{c}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\left(-\nabla p_{\epsilon }+\mu \nabla \cdot \mathbb{D}\left(u\right)\right)\cdot vdV=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV\end{array}$ (9. b)

$\begin{array}{c}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla p_{\epsilon }\cdot vdV-\mu \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla \cdot \mathbb{D}\left(u\right)\cdot vdV=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV\end{array}$ (9. c)

$$\begin{array}{c}\begin{array}{l}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\nabla \cdot vdV+\underset{{\Gamma }_{\mathrm{2,4}}}{{\displaystyle \int }}{p}_{\epsilon }v\cdot ndS+2\mu \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u_{\epsilon }\right)\cdot \\ \mathbb{D}\left(v\right)dV+2\mu \underset{{\Gamma }_{\mathrm{2,4}}}{{\displaystyle \int }}\mathbb{D}\left(u_{\epsilon }\right)\cdot v\cdot ndS\end{array}\\ =\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV,at{\Omega }_{\epsilon }\cup {\Gamma }_2\cup {\Gamma }_4\forall v\in H^1\left({\Omega }_{\epsilon }\right)\end{array}$$ (9. d)

The above expression can be further condensed by joining the terms on the boundary, such as:

$\begin{array}{c}\begin{array}{l}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\nabla \cdot vdV+2\mu \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u_{\epsilon }\right)\cdot \mathbb{D}\left(v\right)dV-\\ \underset{{\Gamma }_{\mathrm{2,4}}}{{\displaystyle \int }}\left(-p_{\epsilon }I+2\mu \mathbb{D}\left(u_{\epsilon }\right)\right)\cdot vdS=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV\end{array}\end{array}$ (9. e)

By decomposing $v\in H^1\left({\Omega }_{\epsilon }\right)$ at the boundary into its normal and tangential components:

$\begin{array}{c}v=v\cdot n+{\displaystyle \sum }_{i=1}^{N-1}v\cdot t,v\in H^1{\left({\Omega }_{\epsilon }\right)}^N\end{array}$ (10)

The variational formulation becomes:

$\begin{array}{c}\begin{array}{l}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\nabla \cdot vdV+2\mu \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u_{\epsilon }\right)\cdot \mathbb{D}\left(v\right)dV-\\ \underset{{\Gamma }_{\mathrm{2,4}}}{{\displaystyle \int }}{n}^T\left(-p_{\epsilon }I+2\mu \mathbb{D}\left(u_{\epsilon }\right)\right)nv\cdot ndS\end{array}\\ +\underset{{\Gamma }_{\mathrm{2,4}}}{{\displaystyle \int }}{n}^T\left(-p_{\epsilon }I+2\mu \mathbb{D}\left(u_{\epsilon }\right)\right)tv\cdot tdS=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV\end{array}$ (11)

Where by the impenetrability of the boundary, the variational is reduced to:

$\begin{array}{c}\begin{array}{l}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\nabla \cdot vdV+2\mu \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u_{\epsilon }\right)\cdot \mathbb{D}\left(v\right)dV+\\ \underset{{\Gamma }_{\mathrm{2,4}}}{{\displaystyle \int }}{n}^T\left(-p_{\epsilon }I+2\mu \mathbb{D}\left(u_{\epsilon }\right)\right)tv\cdot tdS\end{array}\\ =\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV\end{array}$ (12)

Being, by condition imposed on the border, $\beta =0$

$\begin{array}{c}\underset{{\Gamma }_{\mathrm{2,4}}}{{\displaystyle \int }}{n}^T\left(-p_{\epsilon }I+2\mu \mathbb{D}\left(u_{\epsilon }\right)\right)tv\cdot tdS=0\end{array}$ (13)

Thus, at ${\Gamma }_2\cup {\Gamma }_4+{\Omega }_{\epsilon }$:

$\begin{array}{c}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\nabla \cdot vdV+2\mu \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u_{\epsilon }\right)\cdot \mathbb{D}\left(v\right)dV=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV\end{array}$ (14)

With that, one has the (apparent) result that such conditions on the boundary of $\Omega$ simulate (within a few possible adjustments) the condition of Allaire (1991a) on ${ℝ}^N$ since the Dirichlet conditions are just inputs of “velocity at the infinity” and the slip conditions recover the Stokes partial differential equation itself. Soon, the problem becomes just a Stokes flow around the cylinder, from where we can henceforth analyze the grain size problem and the emergence of “strange terms”.

With the variational defined in ${\Omega }_{\epsilon }$, $\Omega$ (by extension to be used later) and its boundaries, it is assumed that $v=w_{\epsilon }^k\varphi$ and $q=q_{\epsilon }^k\varphi$ with hypotheses similar to those of Allaire (1991a). It is verified that the variational generated above converge to a Brinkman-type equation, regardless of whether it is on the boundary or inside the domain.

For the variational referring to ${\Omega }_{\epsilon }$:

$\begin{array}{c}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }{w}_{\epsilon }^k\cdot \nabla \varphi dV+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla u_{\epsilon }\cdot \\ w_{\epsilon }^k\nabla \varphi dV+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\varphi \nabla u_{\epsilon }\cdot \nabla w_{\epsilon }^{k}dV\end{array}$

$=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot \varphi w_{\epsilon }^{k}dV,\forall v\in H^1{\left({\Omega }_{\epsilon }\right)}^N$, at Ω_ε

$=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot \varphi w_{\epsilon }^{k}dV,\forall v\in H^1{\left({\Omega }_{\epsilon }\right)}^N$, at Ω_ε (P2)

Note that the third term to the left of the equality has the following identity:

$\begin{array}{c}\begin{array}{l}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\varphi \nabla u_{\epsilon }\cdot \nabla w_{\epsilon }^{k}dV+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{u}_{\epsilon }\nabla \varphi \cdot \nabla w_{\epsilon }^{k}dV=\\ -\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\varphi u_{\epsilon }\cdot \Delta w_{\epsilon }^{k}dV\end{array}\end{array}$ (15)

Up to one term on the boundary that disappears because it is a functional on ${\Omega }_{\epsilon }$. Therefore, the same term can be replaced by:

$\begin{array}{c}\begin{array}{l}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\varphi \nabla u_{\epsilon }\cdot \nabla w_{\epsilon }^{k}dV=-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{u}_{\epsilon }\nabla \varphi \cdot \nabla w_{\epsilon }^{k}dV-\\ \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\varphi u_{\epsilon }\cdot \Delta w_{\epsilon }^{k}dV\end{array}\end{array}$ (16)

Integrating by parts the continuity equation (div-free condition):

$\begin{array}{c}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla \left(\varphi q_{\epsilon }^k\right)\cdot u_{\epsilon }dV=\mathrm{0,}\forall q\in L^2\left({\Omega }_{\epsilon }\right)\end{array}$, at Ω_ε(17)

Expanding it as:

$\begin{array}{c}\begin{array}{l}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\varphi \nabla q_{\epsilon }^k\cdot u_{\epsilon }dV+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{q}_{\epsilon }^k\nabla \varphi \cdot u_{\epsilon }=\\ \mathrm{0,}\forall q\in L^2\left({\Omega }_{\epsilon }\right)\backslash R\end{array}\end{array}$, at Ω_ε(18)

Because we are in a strictly real Hilbert space, the order of $\varphi$ in the first integral does not change the result, so:

$\begin{array}{c}\begin{array}{l}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla q_{\epsilon }^k\cdot \varphi u_{\epsilon }dV+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{q}_{\epsilon }^k\nabla \varphi \cdot u_{\epsilon }=\mathrm{0,}\\ \forall q\in L^2\left({\Omega }_{\epsilon }\right)\end{array}\end{array}$, at Ω_ε(19)

The continuity equation is part of the PDE system and therefore, by linearity, the Stokes equation can be added, thus generating the PDE:

$$\begin{array}{l}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }{w}_{\epsilon }^k\cdot \nabla \varphi dV+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla u_{\epsilon }\cdot w_{\epsilon }^k\nabla \varphi dV\\ -\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{u}_{\epsilon }\nabla \varphi \cdot \nabla w_{\epsilon }^{k}dV-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\varphi u_{\epsilon }\cdot \Delta w_{\epsilon }^{k}dV\end{array}$$

$+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla q_{\epsilon }^k\cdot \varphi u_{\epsilon }dV+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{q}_{\epsilon }^k\nabla \varphi \cdot u_{\epsilon }dV$

$\begin{array}{c}=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot \varphi w_{\epsilon }^{k}dV,\forall v\in H^1{\left({\Omega }_{\epsilon }\right)}^N\end{array}$, at Ω_ε(20)

The PDE (20) is only in ${\Omega }_{\epsilon }$. The decision is then made to use the linear extension operator to evaluate the behavior of the PDE in $\Omega$. By the sixth hypothesis, when applying the operator in all terms of EDP (6.20), one recovers:

$\begin{array}{l}-\underset{\Omega }{{\displaystyle \int }}{P}_{\epsilon }\left(p_{\epsilon }\right)w_{\epsilon }^k\cdot \nabla \varphi dV+\underset{\Omega }{{\displaystyle \int }}\nabla \widetilde{u_{\epsilon }}\cdot w_{\epsilon }^k\nabla \varphi dV-\\ \underset{\Omega }{{\displaystyle \int }}\widetilde{u_{\epsilon }}\nabla \varphi \cdot \nabla w_{\epsilon }^{k}dV-\underset{\Omega }{{\displaystyle \int }}\varphi \widetilde{u_{\epsilon }}\cdot \Delta w_{\epsilon }^{k}dV\end{array}$

$+\underset{\Omega }{{\displaystyle \int }}\nabla q_{\epsilon }^k\cdot \varphi \widetilde{u_{\epsilon }}dV+\underset{\Omega }{{\displaystyle \int }}{q}_{\epsilon }^k\nabla \varphi \cdot \widetilde{u_{\epsilon }}dV$

$\begin{array}{c}=\underset{\Omega }{{\displaystyle \int }}f\cdot \varphi w_{\epsilon }^{k}dV,\forall v\in H^1{\left({\Omega }_{\epsilon }\right)}^N\end{array}$, at Ω(21)

Regrouping some terms, specifically the fourth and fifth term to the left of the equality:

$\begin{array}{l}-\underset{\Omega }{{\displaystyle \int }}{P}_{\epsilon }\left(p_{\epsilon }\right)w_{\epsilon }^k\cdot \nabla \varphi dV+\underset{\Omega }{{\displaystyle \int }}\nabla \widetilde{u_{\epsilon }}\cdot w_{\epsilon }^k\nabla \varphi dV\\ -\underset{\Omega }{{\displaystyle \int }}\widetilde{u_{\epsilon }}\nabla \varphi \cdot \nabla w_{\epsilon }^{k}dV\end{array}$

$+\underset{\Omega }{{\displaystyle \int }}\left(\nabla q_{\epsilon }^k-\Delta w_{\epsilon }^k\right)\cdot \varphi \widetilde{u_{\epsilon }}dV+\underset{\Omega }{{\displaystyle \int }}{q}_{\epsilon }^k\nabla \varphi \cdot \widetilde{u_{\epsilon }}dV$

$\begin{array}{c}=\underset{\Omega }{{\displaystyle \int }}f\cdot \varphi w_{\epsilon }^{k}dV,\forall v\in H^1{\left(\Omega \right)}^N\end{array}$, at Ω(22)

Where one immediately recognizes the pair of dualities:

$\begin{array}{c}\begin{array}{l}\underset{\Omega }{{\displaystyle \int }}\left(\nabla q_{\epsilon }^k-\Delta w_{\epsilon }^k\right)\cdot \varphi \widetilde{u_{\epsilon }}dV=\\ (\nabla q_{\epsilon }^k-\Delta w_{\epsilon }^k,\varphi {\widetilde{u_{\epsilon })}}_{H^{-1}\left(\Omega \right),H^1\left(\Omega \right)}\end{array}\end{array}$ (23)

Let the following convergences (C) be:

$\begin{array}{c}\widetilde{u_{\epsilon }}\rightharpoonup u\in H^1{\left(\Omega \right)}^N\text{weakly}\end{array}$ (C1)

$\begin{array}{c}{w}_{\epsilon }^k\rightharpoonup e_k\in H^1{\left(\Omega \right)}^N\text{weakly}\end{array}$ (C2)

$\begin{array}{c}{q}_{\epsilon }^k\rightharpoonup 0\in L^2\left(\Omega \right)\text{weakly}\end{array}$ (C3)

$\begin{array}{c}{P}_{\epsilon }\left(p_{\epsilon }\right)\rightharpoonup p\in L^2\left(\Omega \right)\text{weakly}\end{array}$ (C4)

Finally, the functional on ${\Omega }_{\epsilon }$ converges to a functional on $\Omega$ in the form:

$-\underset{\Omega }{{\displaystyle \int }}p{e}_k\cdot \nabla \varphi dV+\underset{\Omega }{{\displaystyle \int }}\nabla u\cdot e_k\nabla \varphi dV$

$+{\mu }_k,\varphi u_{H^{-1}\left(\Omega \right),H^1\left(\Omega \right)}$

$\begin{array}{c}=\underset{\Omega }{{\displaystyle \int }}f\cdot \varphi e_{k}dV,\forall v\in H^1{\left(\Omega \right)}^N\end{array}$, at Ω(24)

Since $\nabla e_k=0$. Integrating by parts and passing to the limit in $\Omega$:

$\underset{\Omega }{{\displaystyle \int }}\nabla p\cdot \varphi dV-\underset{\Omega }{{\displaystyle \int }}\Delta u\cdot \varphi dV$

$$+{(\varphi {\mu }_k,u)}_{H^{-1}\left(\Omega \right),H^1\left(\Omega \right)}$$

$\begin{array}{c}=\underset{\Omega }{{\displaystyle \int }}f\cdot \varphi e_{k}dV,\forall \varphi \in D\left(\Omega \right)\end{array}$, at Ω(25)

Thus:

$\begin{array}{c}\nabla p-\Delta u+Mu=f\end{array}$, at Ω (PDE. a)

$\begin{array}{c}\nabla \cdot u=0\end{array}$, at Ω (PDE. b)

Where M is independent of u.

Now, let’s focus on the free-sliding boundaries. Let be the variational formulation on the boundary:

$\begin{array}{c}\begin{array}{l}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\nabla \cdot vdV+2\mu \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u_{\epsilon }\right)\cdot \mathbb{D}\left(v\right)dV=\\ \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot vdV\end{array}\end{array}$ (26)

Replacing $v$ with its previously established definition:

$\begin{array}{c}\begin{array}{l}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\nabla \cdot \left(\varphi w_{\epsilon }^k\right)dV+2\mu \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u_{\epsilon }\right)\cdot \mathbb{D}\left(\varphi w_{\epsilon }^k\right)dV=\\ \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot \varphi w_{\epsilon }^{k}dV\end{array}\end{array}$ (27)

The first term to the left of the equality has already been discussed previously for the functional on ${\Omega }_{\epsilon }$. We must now analyze the second term to the left of the equality, referring to the stress tensor, specifically its deviatoric part. At this point, it seems difficult to extract any property from the functional above that might resemble the results obtained so far. However, it is enough to go back a few steps in the development of the variational formulation so that one can take advantage of the incompressibility of the test function field and the velocity field itself.

It is known that:

$\begin{array}{c}2D\left(u_{\epsilon }\right)\cdot D\left(\varphi w_{\epsilon }^k\right)=2D\left(u_{\epsilon }\right)\cdot D\left(v\right)=D\left(u_{\epsilon }\right)\cdot \nabla v\end{array}$ (28)

From this, the following integration by parts follows naturally:

$$\begin{array}{c}\begin{array}{l}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u_{\epsilon }\right)\cdot \nabla vdV=-\frac{1}{2}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla \cdot \left(\nabla u_{\epsilon }+\nabla u_{\epsilon }^T\right)\cdot \\ vdV+\frac{1}{2}\underset{\partial {\Omega }_{\epsilon }={Ã}_{\mathrm{2,4}}}{{\displaystyle \int }}\left(\nabla u_{\epsilon }+\nabla u_{\epsilon }^T\right)\cdot v\cdot ndS\end{array}\end{array}$$ (29)

Where the last term to the right of the equality vanishes because $v\cdot n=0$ in ${\Gamma }_{\mathrm{2,4}}$. Regarding the first term to the right of the equality, it is known that:

$\begin{array}{c}\nabla \cdot \nabla u_{\epsilon }=\Delta u_{\epsilon }\end{array}$ (30. a)

$\begin{array}{c}\nabla \cdot \nabla u_{\epsilon }^T=\nabla \left(\nabla \cdot u_{\epsilon }\right)=0\end{array}$ (30. b)

Then:

$\begin{array}{c}\begin{array}{l}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mathbb{D}\left(u_{\epsilon }\right)\cdot \nabla vdV=-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\Delta u_{\epsilon }\cdot vdV=\\ -\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\Delta u_{\epsilon }\cdot \varphi w_{\epsilon }^{k}dV\end{array}\end{array}$ (31)

Then, the variational formulation is rewritten such that:

$\begin{array}{c}\begin{array}{l}-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\nabla \cdot \left(\varphi w_{\epsilon }^k\right)dV-\mu \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\Delta u_{\epsilon }\cdot \varphi w_{\epsilon }^{k}dV=\\ \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot \varphi w_{\epsilon }^{k}dV\end{array}\end{array}$ (32)

Where Stokes is recovered for every test function $v=\varphi w_{\epsilon }^k$. Once you have the Stokes system in hand, for such test functions, ones arrive at the same homogenized equation:

$\begin{array}{c}\nabla p-\Delta u+Mu=f,\mu =1\end{array}$ (33)

This time, developed on the boundary and with the same matrix M. This result is the core of the question about the compatibility of variationals on the boundaries and in the domain, for non-Dirichlet boundary conditions.

The results obtained refer to the open $\Omega$ problem composed of ${\Omega }_{\epsilon }$ and $T_{\epsilon }^i$. According to the procedure established by Allaire (1991a), the cell problem must be a microscale representation of the global problem, that is, the Stokes problem must be solved in a domain whose boundary conditions reflect (or simulate) the conditions in the macroscopic problem. For a bounded open domain, how should the cell problem be formulated? This issue will be addressed in the next section and, with it, results on the limits of M in relation to grain size.

2.1. The cell problem

Paraphrasing Allaire (1991a), “the functions $w_{\epsilon }^k\in H^1\left({\Omega }_{\epsilon }\right)$ and $q_{\epsilon }^k\in L^2\left({\Omega }_{\epsilon }\right)$seem mysterious at first glance. But these have physical meaning.” In fact, these need to make physical sense and this condition will be dealt with in this session. Unlike a generic variational problem, the test function in an energy homogenization problem is a specific function, carefully chosen to satisfy both mathematical (regularity) and physical (coherence with the model) criteria.

Despite the numerical nature of this work, the analysis of the convergence of the Stokes problem for one of its three possible limits – Stokes, Brinkman or Darcy – must be done from the analytical solution. This subsection will then be dedicated to presenting the hypotheses that will allow the resolution of the Stokes system in the unit cell, inspired by the hypotheses presented by Allaire (1991a) and Cionarescu & Murat (1982). Let the following differential problem be:

Find $\left(w_{\epsilon }^k,q_{\epsilon }^k\right)$ solutions of the system:

$\begin{array}{c}\nabla q_{\epsilon }^k-\Delta w_{\epsilon }^k=0,at{\Omega }_{\epsilon }^{unit}\end{array}$ (Cell - problem. a)

$\begin{array}{c}\nabla \cdot w_{\epsilon }^k,at{\Omega }_{\epsilon }^{unit}\end{array}$ (Cell - problem. b)

With the following conditions (Cond.):

$\begin{array}{c}\nabla q_{\epsilon }^k-\Delta w_{\epsilon }^k=0,at{B}_{\epsilon }^k\end{array}$ (Cond.1)

$\begin{array}{c}{w}_{\epsilon }^k=0,at{T}_{\epsilon }^k\end{array}$ (Cond.2)

$\begin{array}{c}{w}_{\epsilon }^k=U_{\infty },at{D}_{\epsilon }^k\end{array}$ (Cond.3)

$\begin{array}{c}{w}_{\epsilon }^k=U_{\infty },at\partial {\Omega }_{\epsilon }^{unit}\end{array}$ (Cond.4)

But here, consider the problem in polar (radial) coordinates and that the pressure field has a previously known “behavior”. A periodic pressure field will be adopted such that:

$\begin{array}{c}{q}_{\epsilon }^k={\displaystyle \sum }_{k=1}^{\infty }{q}_{\epsilon }{e}^{\frac{ir}{\epsilon }k}=f\left(r\right)\end{array}$ (34)

Symmetric (does not depend on φ or θ) and k is the spatial frequency of the solution. Hence, the problem to be solved becomes a Poisson problem on the ring, such that the solution satisfies:

$\begin{array}{c}\Delta w_{\epsilon }^k=f\left(r\right),at{B}_{\epsilon }^k\end{array}$ (Cond.5)

$\begin{array}{c}{w}_{\epsilon }^k=0,at{T}_{\epsilon }^k\end{array}$ (Cond.6)

$\begin{array}{c}{w}_{\epsilon }^k=U_{\infty },at{D}_{\epsilon }^k\end{array}$ (Cond.7)

$\begin{array}{c}{w}_{\epsilon }^k=U_{\infty },at\partial {\Omega }_{\epsilon }^{unit}\end{array}$ (Cond.8)

The problem can be further rewritten as:

$\begin{array}{c}\frac{{\partial }^2{w}_{\epsilon }^k}{\partial r^2}+\frac{1}{r}\frac{\partial w_{\epsilon }^k}{\partial r}=f\left(r\right),at{B}_{\epsilon }^k\end{array}$ (35. a)

Alternatively

$\begin{array}{c}\frac{{\partial }^2{w}_{\epsilon }^k}{\partial r^2}+\frac{1}{r}\frac{\partial w_{\epsilon }^k}{\partial r}={\displaystyle \sum }_{k=1}^{\infty }{q}_{\epsilon }{e}^{\frac{ir}{\epsilon }k},at{B}_{\epsilon }^k\end{array}$ (35. b)

With same conditions. The right-hand side of the PDE can be rewritten according to the Voss-Weyl formula such that:

$\begin{array}{c}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial w_{\epsilon }^k}{\partial r}\right)={\displaystyle \sum }_{k=1}^{\infty }{q}_{\epsilon }{e}^{\frac{ir}{\epsilon }k},at{B}_{\epsilon }^k\end{array}$ (36)

The general solution becomes easier to find from this approach, as follows:

$\begin{array}{c}\begin{array}{l}{w}_{\epsilon }^k\left(r,\theta \right)=-\frac{\epsilon i}{k}{E}_i\left(\frac{ikr}{\epsilon }\right)+\frac{{\epsilon }^2}{k^2}\\ E_i\left(\frac{ikr}{\epsilon }\right)+C_1\log \left(r\right)+C_2\end{array}\end{array}$ (37. a)

Where $E_i\left(r\right)$ is the exponential-integral function that can be expressed by the convergent Ramanujan series (for small r):

$\begin{array}{c}\begin{array}{l}{E}_i\left(r\right)=\gamma +\ln \left(r\right)+e^{\frac{r}{2}}{\displaystyle \sum }_{n=1}^{\infty }\frac{{\left(-1\right)}^{n-1}{r}^n}{n!2^{n-1}}\\ {\displaystyle \sum }_{k=0}^{\frac{n-1}{2}}\frac{1}{2k+1}\end{array}\end{array}$ (37. b)

with $\gamma$ being the Euler-Mascheroni constant.

Note that applying the boundary conditions on $r=a_{\epsilon }$ and $R=\epsilon$, considering that:

$\begin{array}{c}{e}^{\frac{r}{2}}{\displaystyle \sum }_{n=1}^{\infty }\frac{{\left(-1\right)}^{n-1}{r}^n}{n!2^{n-1}}{\displaystyle \sum }_{k=0}^{\frac{n-1}{2}}\frac{1}{2k+1}\to \alpha r^p{e}^{\frac{r}{2}}\end{array}$ (38)

For large n, k, a result of the form is expected:

$\begin{array}{c}{w}_{\epsilon }^k\left(r,\theta \right)={\displaystyle \sum }_{k=1}^{\infty }\begin{array}{c}\begin{array}{l}-\frac{\epsilon i}{k}\left(\gamma +\ln \left(\frac{ikr}{\epsilon }\right)+\alpha ik{r}^p{e}^{\frac{ikr}{2\epsilon }}\right)\\ +\frac{{\epsilon }^2}{k^2}\left(\gamma +\ln \left(\frac{ikr}{\epsilon }\right)+\alpha ik{r}^p{e}^{\frac{ikr}{2\epsilon }}\right)\end{array}\\ +C_1\log \left(r\right)+C_2\end{array}\end{array}$ (39)

Simplifying:

$\begin{array}{c}\begin{array}{l}{w}_{\epsilon }^k\left(r,\theta \right)={\displaystyle \sum }_{k=1}^{\infty }\left(-\frac{\epsilon i}{k}+\frac{{\epsilon }^2}{k^2}\right)\\ \left(\gamma +\ln \left(\frac{ikr}{\epsilon }\right)+\alpha ik{r}^p{e}^{\frac{ikr}{2\epsilon }}\right)+C_1\log \left(r\right)+C_2\end{array}\end{array}$ (40)

generates:

$$\begin{array}{c}{w}_{\epsilon }^k\left(r,\theta \right)={\displaystyle \sum }_{k=1}^{\infty }\begin{array}{c}\left(-\frac{\epsilon i}{k}+\frac{{\epsilon }^2}{k^2}\right)\left(\begin{array}{l}\gamma +\ln \left(\frac{ik\beta a_{\epsilon }}{{\epsilon }^2}\right)+\\ \alpha ik{\left(\frac{\beta a_{\epsilon }}{\epsilon }\right)}^p{e}^{\frac{ik\beta a_{\epsilon }}{2{\epsilon }^2}}\end{array}\right)\\ +C_1\log \left(\frac{\beta a_{\epsilon }}{\epsilon }\right)+C_2\end{array}\end{array}$$ (41)

For $\beta \in ℝ$. Where it is noted that the bounding of the solution will explicitly depend on the relationship between the radii (grain and unit cell) and that this ratio will be controlled by the product of r by the logarithm of the ratio. This result is very similar to that found by Allaire (1990, 1991a). Here, precisely, one can estimate the convergence for the limits 0, +∞ or $0<C_0<+\infty$:

$\begin{array}{c}\epsilon {\left|\ln \left(\frac{a_{\epsilon }}{\epsilon }\right)\right|}^{\frac{1}{p}},p\in {ℝ}^{+}\end{array}$ (42)

3. The particle size distribution function and the grain size decay

In geotechnical engineering, it is common to classify soils based on their granulometric distribution, that is, to analyze the composition of the soil according to the size (diameter) of the grains that compose it. Figure 4 is an example of this type of test.

Figure 4
Particle size distribution of soils

Some parameters can be extracted from this type of analysis, such as the coefficient of uniformity ($C_u$) and the coefficient of curvature ($C_c$), calculated as:

$\begin{array}{c}{C}_u=\frac{D_{60}}{D_{10}}\end{array}$ (43. a)

$\begin{array}{c}{C}_c=\frac{D_{30}^2}{D_{60}\times D_{10}}\end{array}$ (43. b)

Where D₆₀, D₁₀ and D₃₀ indicate that 60%, 10% and 30% of the grains have a diameter smaller than X mm. Both parameters $C_u$ and $C_c$ characterize how uniform a soil is, that is, how well distributed its granulometry is. The ABNT-NBR-6502/22 (ABNT, 2022) proposes that a soil is considered uniform if $C_u<5$. ASTM-D2487-17 (ASTM, 2017) proposes that a soil is uniform if $C_c>3$ or $C_c<1$. With these definitions known, we propose the following ansatz that will be justified later:

Ansatz (grain size decay): the function a(ε) used in homogenization theory can be represented as:

$\begin{array}{c}a\left(\epsilon \right)=e^{-\frac{C_u}{{\epsilon }^{C_c}}}\end{array}$ (44)

Where ε must be understood as the size of the representative elementary volume.

Applying Equation 44 directly to the critical grain size Equation 42, using p = 2 as suggested by Allaire (1991a), taking ε→0, we have:

$\begin{array}{c}{\sigma }^2\to 0,C_c<2(\text{Darcy)}\end{array}$ (45. a)

$\begin{array}{c}{\sigma }^2\to +\infty ,C_c>2\text{(Navier-Stokes)}\end{array}$ (45. b)

$\begin{array}{c}{\sigma }^2\to 0<C_0<+\infty ,C_c=2\text{(Brinkman)}\end{array}$ (45. c)

It is now up to prove that these limits make some physical sense in this new proposed approach. Before proving, the following geotechnical analogue can be made:

• Situation 1: $C_c<2$ implies a well-graded soil, making it possible to observe grains for any sample size. Therefore, the voids of the larger grains are occupied by smaller grains, increasing the total specific surface, consequently increasing the friction of the fluid in the solids, and preventing the development of nonlinearities at high speeds. Darcy's law is then formed.

• Situation 2: $C_c>2$ implies a uniform soil, poorly graded, meaning that from a certain size of ε small enough, no more grains are observed, but free fluid. The pores become larger, and the total specific surface does not reach high values. The fluid percolates almost freely, similar to flow governed directly by Navier-Stokes. The term “similar” here becomes necessary as it will later be proved that the presence of grains will generate an apparent viscosity for the flow.

• Situation 3: $C_c=2$ implies a medium uniform soil. It is the critical distribution for it to be neither a free flow nor a diffusion flow as envisaged by Darcy. There is a flow of the Navier-Stokes type, but with the presence of a drag term, multiplying the velocity, as a form of damping generated by the presence of some grains. However, nonlinear terms remain relevant.

It is relevant to comment that, physically, a sample with size ε→0 does not exist. This problem will also appear when the problem is treated numerically, since there is a computational limit that does not allow taking the grain size to zero.

So let Equation 42 be rewritten as:

$\begin{array}{c}\sigma \left(\epsilon \right)=\epsilon {\left|\ln \left(\frac{a_{\epsilon }}{\epsilon }\right)\right|}^{\frac{1}{2}}\end{array}$ (46)

The different limits of (48) are represented in Figure 5.

Figure 5
Limits for a(ε) and σ(ε).

4. Numerical convergence of the homogenization problem of Navier Stokes equations in two dimensions

It is proposed in this brief introduction that, for reasons of numerical limitation, Allaire's theory of homogenization (1991a) had to be modified in fundamental aspects, affecting the hypotheses about the test functions that will no longer be able to adopt certain limits.

4.1 The first steps and first struggles

Let the following problem be numerically homogenized by the conventional finite element method:

Find $\left(u_{\epsilon },p_{\epsilon }\right)\in {\left[H_0^1\left({\Omega }_{\epsilon }\right)\right]}^N\times L^2\left({\Omega }_{\epsilon }\right)$ solutions of the problem (P1-Num)

$\begin{array}{c}\rho u_{\epsilon }\left(u_{\epsilon }\cdot \nabla \right)+\nabla p_{\epsilon }-\Delta u_{\epsilon }=f,at{\Omega }_{\epsilon }\end{array}$ (P1 - Num.1)

$\begin{array}{c}\nabla \cdot u_{\epsilon }=0,at{\Omega }_{\epsilon }\end{array}$ (P1 - Num.2)

$\begin{array}{c}{p}_{\epsilon }=p_{in},at{\Gamma }_1\end{array}$ (P1 - Num.3)

$\begin{array}{c}{p}_{\epsilon }=p_{out},at{\Gamma }_3\end{array}$ (P1 - Num.4)

$\begin{array}{c}{u}_{\epsilon }=\mathrm{0,}at{\Gamma }_{2}e{\Gamma }_4\end{array}$ (P1 - Num.5)

So that when ε→0, $\left(u_{\epsilon },p_{\epsilon }\right)\rightharpoonup \left(u,p\right)$ weakly, in some sense. In the theory of analysis in abstract spaces, it was possible to set the size of the grain, or the hole, to zero. However, when dealing with the numerical problem, the grain size cannot be taken to zero.

It is known that a fluid loses energy when percolating through a porous medium by friction with the solid walls. Therefore, it is known that the greater the specific surface of the grain, the greater the friction with the fluid and the greater the dissipation of energy during percolation. A greater loss of energy, induces in geotechnical theory, a lower permeability coefficient, if Darcy’s law is assumed as true. The numerical consequence is relevant: since it is not possible to disappear with the grains in the limit ε→0, each simulation with smaller grains increases considerably the specific surface no matter how much the porosity is maintained.

In Allaire's theory (1991a), a cell problem with a single obstacle is proposed, in order to represent in a unitary cell, the same porosity of the total domain, since this is just an infinite and periodic repetition of these cells. In the numerical treatment carried out here, it was necessary to find a way to introduce the specific surface of the grains and the porosity in the unit cell. How to do?

4.2 The multigrain cell problem

With the mathematical impasse raised in the previous paragraph, it became necessary to propose a new way of representing the unit cell, contemplating the porosity and the specific surface of the grains. The proposal is to use more than one grain in the unit cell periodically and symmetrically distributed. To do this, simply solve the following nonlinear system of two algebraic equations:

$\begin{array}{c}1-n\pi r^2=\varphi \end{array}$ (47. a)

$\begin{array}{c}\frac{2n}{r}=M_s^{unit}\end{array}$ (47. b)

Where (49.a) is the porosity for a square domain of unitary sides and (47.b) represents the unit specific surface of the grain. This second equation needs a more detailed explanation.

Suppose you want to generate a unit cell with only one grain that meets the requirements of porosity and specific surface. Also suppose that this problem has a solution for a certain radius r of the grain. Therefore, the specific surface of this grain is calculated, in two dimensions, as:

$\begin{array}{c}{M}_s^{unit}=\frac{A_s}{V_s}=\frac{2\pi r}{\pi r^2}=\frac{2}{r}\end{array}$ (48)

But it was said earlier that it is unlikely, unless a coincidence, that a radius r fulfills the requirements of porosity and specific surface. Therefore, it becomes natural to question whether it is possible to divide this grain into grains of smaller radius so that a greater number of grains fulfills the two aforementioned requirements.

A cell with more than one grain is generated, as represented by Figure 6.

Figure 6
Representation of the unit cell with more than one grain.

In the unit cell, the following problem is solved:

Find $\left(w,q\right)$solutions of the problem (P1-Num)

$\begin{array}{c}\rho w\left(w\cdot \nabla \right)+\nabla q-\Delta w=f,at{\Omega }_{unit}\end{array}$ (P_Cel - Num.1)

$\begin{array}{c}\nabla \cdot w=0,at{\Omega }_{unit}\end{array}$ (P_Cel - Num.2)

$\begin{array}{c}w=w_0,at{\Gamma }_1\end{array}$ (P_Cel - Num.3)

$\begin{array}{c}w=w_0,at{\Gamma }_3\end{array}$ (P_Cel - Num.4)

$\begin{array}{c}{n}^T\cdot \left(2\mu \mathbb{D}\left(w\right)-qI\right)\cdot t-\beta w\cdot t=0,at{\Gamma }_{2}e{\Gamma }_4\end{array}$ (P_Cel - Num.5)

It is notable that $\left(w,q\right)$ are analogous to the velocity and pressure fields mentioned on the global scale ${\Omega }_{\epsilon }$. The domain in which one wants to solve the problem ($P_{Cel}-Num$) is graphically represented in Figure 7.

Figure 7
$P_{Cel}-Num$ domain.

4.3 Homogenization via direct numerical simulation (DNS)

The objective of this subsection is to investigate the convergence, in whatever sense, of a homogenization process by direct numerical simulation, that is, gradually decreasing the size of the grains in relation to the domain, increasing their quantity to maintain the porosity and respecting a specific grain size decay function.

Let the problem (P1-Num) be the one we want to solve using the conventional finite element method. The momentum conservation equations can be written in a variational form such that:

$\begin{array}{c}\begin{array}{l}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\rho u_{\epsilon }\left(u_{\epsilon }\cdot \nabla \right)\cdot \text{φ}d{\Omega }_{\epsilon }+\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\nabla p_{\epsilon }\cdot \phi d{\Omega }_{\epsilon }-\\ \underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mu \Delta u_{\epsilon }\cdot \phi d{\Omega }_{\epsilon }=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot \phi d{\Omega }_{\epsilon }\end{array}\end{array}$ (49. a)

where, integrating by parts:

$$\begin{array}{c}\begin{array}{l}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\rho u_{\epsilon }\left(u_{\epsilon }\cdot \nabla \right)\cdot \phi d{\Omega }_{\epsilon }-\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\nabla \cdot \phi d{\Omega }_{\epsilon }\\ +\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\mu \nabla u_{\epsilon }\cdot \nabla \phi d{\Omega }_{\epsilon }\end{array}\\ +\underset{\partial {\Omega }_{\epsilon }}{{\displaystyle \int }}{p}_{\epsilon }\phi \cdot ndS+\underset{\partial {\Omega }_{\epsilon }}{{\displaystyle \int }}\mu \nabla u_{\epsilon }\cdot \phi \cdot ndS=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}f\cdot \phi d{\Omega }_{\epsilon }\end{array}$$ (49. b)

One then has $u_{\epsilon }\in H^1\left({\Omega }_{\epsilon }\right)$ and $p_{\epsilon }\in L^2\left({\Omega }_{\epsilon }\right)$, at least. Let $V\subset H^1\left({\Omega }_{\epsilon }\right)$ be a linear subspace. Having V finite dimension, taking advantage of the properties of orthogonal projection in Hilbert spaces, an element of V is proposed as:

$\begin{array}{c}{u}_{\epsilon }^{\text{'}}={\displaystyle \sum }_{n=0}^N{U}_n{\Phi }_n,N<+\infty \end{array}$ (50)

With $U_n\in {ℝ}^N$ and ${\Phi }_n\in V\subset H^1$. Similarly, there is an element of V:

$\begin{array}{c}{p}_{\epsilon }^{\text{'}}={\displaystyle \sum }_{n=0}^M{P}_n{\Phi }_n,M<N\end{array}$ (51)

With $P_n\in {ℝ}^M$. By the orthogonality of V in relation to $H^1\left({\Omega }_{\epsilon }\right)$, it is known that:

$\begin{array}{c}{u}_{\epsilon }^{\text{'}}-u_{\epsilon }{}_2+p_{\epsilon }^{\text{'}}-p_{\epsilon }\le z-u_{\epsilon }{}_2+w-p_{\epsilon }\end{array}$ (52)

$\forall z,w\in V$

That is, $u_{\epsilon }^{\text{'}}$ and $p_{\epsilon }^{\text{'}}$ minimize the solution approximation error.

The continuity equation can be formulated in its variational form as:

$\begin{array}{c}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\left(\nabla \cdot u_{\epsilon }\right)\Psi d{\Omega }_{\epsilon }=0\end{array}$ (53)

Where $\Psi \in E$, being $E\subset L^2\left({\Omega }_{\epsilon }\right)$.

For the direct numerical simulation, a first-order polynomial function (Lagrange polynomial) was used for ${\Phi }_n$ and a constant, merely continuous function for ${\Psi }_n$. In terms of finite element theory, triangular elements of type T2 were chosen, that is, triangular elements with three nodes per side through which a first-order polynomial for the velocity and a zeroth order polynomial for pressure pass.

The formulation assumed that the characteristic size of the elements was 1% of the domain size, yielding about ${10}^4$ elements in a square domain of unit sides. Close to the holes, the mesh was refined to 1% of the hole perimeter (Figure 8). The non-linearity of the convective term was treated with the Newton-Raphson iterative method. In other words, the residual of each component of the variational formulation was iterated with changes in the nodal value until a minimum was reached. The element used in the approach is shown in Figure 9.

Figure 8
Discretization of the finite element space.

Figure 9
T2-type element.

4.4 Error estimation and the influence of convective nonlinearity: a compact perturbation

To estimate the error due to approximation of the solution, a theorem is needed (Brenner & Scott, 2008):

Theorem (error estimate): Suppose the relation in (8.6) is true. Let $u\in H^k\left(\Omega \right)$ and also $u\in H^s\left(\Omega \right)$. Also let $m<s<k$. Then $\forall u\in H^s\left(\Omega \right)$:

$\begin{array}{c}u-u_h{}_{H^m\left(\Omega \right)}\le C{h}^{s-m}{u}_{H^s\left(\Omega \right)}\end{array}$ (Theorem)

Where $u_h\in M\subset H^m\left(\Omega \right).$

The proof of this theorem can be found in (Brenner & Scott, 2008). The theorem also applies to Sobolev spaces and allows one to obtain:

$\begin{array}{c}u-u_h{}_{W_p^1\left(\Omega \right)}=o\left(h^{s-1}\right)\end{array}$ (54)

Applying it to the problem to be solved in this work, we obtain:

$\begin{array}{c}{u}_{\epsilon }-u_{\epsilon }^{\text{'}}{}_{H^1\left({\Omega }_{\text{ε}}\right)}=o\left(h^{2-1}\right)\end{array}$ (55. a)

Since $u_{\epsilon }\in H^2\left({\Omega }_{\epsilon }\right)$. Assuming a square domain of unit sides and a mesh with a characteristic size of 1% of the domain size:

$\begin{array}{c}{u}_{\epsilon }-u_{\epsilon }^{\text{'}}{}_{H^1\left({\Omega }_{\text{ε}}\right)}=o\left({\mathrm{0,01}}^1\right)\end{array}$ (55. b)

That is, the calculated velocities have an error bounded above by 1% of the actual value.

The convective term is classified as a compact perturbation. In other words, the operator:

$\begin{array}{c}Tu{\text{'}}_{\epsilon }=\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}u{\text{'}}_{\epsilon }\left(u{\text{'}}_{\epsilon }\cdot \nabla \right)\cdot \Phi d{\Omega }_{\epsilon }\end{array}$ (56)

It is said to be compact because it belongs to the space of linear operators of $V\subset H_1\left({\Omega }_{\epsilon }\right)$. The space V is compact and therefore every linear operator in its dual is also compact.

4.5 Simulation script

First, a square domain of unit dimensions of 1m x 1m is constructed. Next, the grain size decay function is chosen. Adopt, for example, the equation given by ansatz (44), such that (42) tends to zero for ε → 0, with p = 2. Take an initial value of ε < 1 (less than the whole domain) and calculate the grain size. Then, choose a porosity ϕ and using relation (47.a), determine how many grains will be needed to fulfill the desired porosity.

A mesh is generated with a characteristic size of 1% of the domain size, as well as the hole perimeter is discretized with 1% of its size to refine the mesh in the vicinity of the hole. The parameters ρ = 1000 kg/m³ and μ = 0,001 Pa⋅s are determined. Using Figure 1 as a reference, the following boundary conditions are imposed:

$\begin{array}{c}{u}_{\epsilon }=0,\text{at}{\Gamma }_2\text{and}{\Gamma }_4\end{array}$ (BC6)

$\begin{array}{c}{p}_{\epsilon }=0Pa,\text{at}{\Gamma }_3\end{array}$ (BC7)

$\begin{array}{c}{p}_{\epsilon }=p_0{}_{\epsilon },\text{at}{\Gamma }_1\end{array}$ (BC8)

Where the value of $p_0{}_{\epsilon }$ is gradually increased to investigate the evolution of the velocity field with the variation of the hydraulic gradient. As the speed varies pointwise, it was determined that the average velocity in the domain would be used as a reference, that is:

$\begin{array}{c}\overline{u_{\epsilon }}=\frac{1}{\left|\Omega \right|}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}{u}_{\epsilon }d{\Omega }_{\epsilon }\end{array}$ (57)

Assuming a flow fundamentally horizontal (gravitational forces are neglected), the hydraulic gradient is written in simplified form as:

$\begin{array}{c}i=\nabla p_{\epsilon }\cong \frac{1}{L}{p}_{\epsilon }\cong \frac{p_0{}_{\epsilon }}{1}=p_0{}_{\epsilon }\end{array}$ (58)

So, one is looking for something like:

$\begin{array}{c}\overline{u_{\epsilon }}=Ki\end{array}$ (59)

The simulation is performed, (57) and (58) are calculated, then the value of K is determined from (59). The value of $p_0{}_{\epsilon }$ is increased and the process is repeated until a curve of the relationship between the average velocity and the hydraulic gradient is constructed. The results are presented in the figures Figure 10 and Figure 11. Note that the behavior of the curve corresponding to a porosity of 0.7 differs from the rest. It is believed to be the beginning of an oscillatory behavior that could not be captured by the numerical scheme used. Note also that in Figure 10, only one of the tests starts with $\epsilon$ close to 1. This is due to the required computational cost for such situation. For such porosity values (0.3), the amount of elements needed for a low enough epsilon exceeds the limits of the machine used for the simulations. Therefore, it was preferred to evaluate the convergence trend for higher epsilon values.

Figure 10
DNS convergence for Darcy.

Figure 11
DNS convergence for Navier-Stokes.

4.6 Homogenization using the multigrain cell problem: the Darcy case

Assume a grain size decay function a(ε) that makes σ(ε)→0 to ε→0. For the sake of simplicity, assume the case ϕ = 0,6 and σ(ε)→0.

First, the number of grains and their radius must be calculated so that the multigrain unit cell is created. Assuming that the direct numerical simulations were able to go up to ε = 0.0833m and that the function that defines the grain size is a(ε)=0.2√ε and $\sigma \left(\epsilon \right)=\epsilon {\left|\ln \left(\frac{a\left(\epsilon \right)}{\epsilon }\right)\right|}^{\frac{1}{2}}$, one has $a\left(\epsilon \right)=a\left(0.0833\right)=0.0577m$.

Knowing that ϕ = 0.6, then the nonlinear system to be solved becomes:

$\begin{array}{c}1-n\pi r^2=0.6\end{array}$ (60. a)

$\begin{array}{c}\frac{2n}{r}=\frac{2}{0.0577}\end{array}$ (60. b)

Resulting in n = 3.369 grains and r = 0.1944 m. As it is not possible to obtain a non-integer number of grains, the value of the nearest integer is assumed, therefore, n = 3. A unitary cell is then built with three grains of radius 0.1944 m, where the problem specified by the equations of $\left(P_{Cel}-Num\right)$ and its due boundary conditions are solved. Once the velocity field solution is in hand, M₀ is calculated, given by:

$\begin{array}{c}{M}_0=\underset{{\Omega }_{unit}}{{\displaystyle \int }}\nabla w:\nabla wd{\Omega }_{unit}\end{array}$ (61)

Where the tensor product of the integrand, in two dimensions, is written in extended form as:

$\begin{array}{c}{M}_0=\underset{{\Omega }_{unit}}{{\displaystyle \int }}{\left({\partial }_{x}u\right)}^2+{\left({\partial }_{y}u\right)}^2+{\left({\partial }_{x}v\right)}^2+{\left({\partial }_{y}v\right)}^{2}d{\Omega }_{unit}\end{array}$ (62)

Where $M_0\in ℝ$. Remember here that the tensor product A:B for dim(A) = m and dim(B) = n has dimension dim(A : B) = m + n - 4.

The value found was $M_0=46.40$. It is known that:

$\begin{array}{c}K={\left(\frac{\mu M_0}{{\sigma }^2}\right)}^{-1}\end{array}$ (63)

Therefore, $\mu M_0=0.001\times 46.40=0.04640$ is already known. It is then enough to find the value of σ(ε) taken directly from the curve:

$\begin{array}{c}\sigma \left(\epsilon \right)=\sigma \left(0.0833\right)=0.0833{\left|\ln \left(\frac{a\left(0.0833\right)}{0.0833}\right)\right|}^{\frac{1}{2}}=0.05045\end{array}$ (64)

thus

$\begin{array}{c}{\sigma }^2={0.05045}^2=0.00254\end{array}$ (65)

Then:

$\begin{array}{c}K={\left(\frac{0.04640}{0.00254}\right)}^{-1}=0.0549\end{array}$ (66)

The mean value found by the direct numerical simulation was 0.0604. Figure 12 shows the convergence of the DNS to the predicted limit.

Figure 12
DNS convergence for Darcy with calculated limit.

5. Homogenization for different Darcy limits: Navier-Stokes apparent viscosity and the Brinkman problem

It has been mentioned that the limits found by Allaire (1991a) cannot be obtained numerically due to limitations of the machines used. For equivalent limits to be reached, modifications in the theory will be proposed here. These modifications act at the core of Allaire's deduction, making terms that once disappeared in his work, gain non-negligible importance in the solution of the partial differential equations presented here.

Let the coefficient $M=\mu M_0$ where $M_0$ is:

$\begin{array}{c}{M}_0=\underset{{\Omega }_{unit}}{{\displaystyle \int }}\nabla w:\nabla wd{\Omega }_{unit}\end{array}$ (67)

Where $w\in H^1\left({\Omega }_{unit}\right)$ and it is known that $v=w\varphi \in H^1\left({\Omega }_{\epsilon }\right)\forall \varphi \in D\left(\Omega \right)$. By the Riesz representation theory, $\exists !z\in H^1\left({\Omega }_{\epsilon }\right)$ such that:

$\begin{array}{c}\underset{{\Omega }_{\epsilon }}{{\displaystyle \int }}\beta z\cdot u_{\epsilon }d{\Omega }_{\epsilon }=\Delta u_{\epsilon }\end{array}$ (68)

For $\beta \in ℝ$. The above result is also verified for the extension in $\Omega$ by the linearity of Hilbert spaces. Realize now that:

$\begin{array}{c}\frac{\mu }{{\sigma }^2}\underset{{\Omega }_{unit}}{{\displaystyle \int }}\nabla w:\nabla wd{\Omega }_{unit}=z\in H^1\left({\Omega }_{unit}\right)\end{array}$ (69)

${\Omega }_{unit}$ being a homeomorphism of ${\Omega }_{\epsilon }$, and the PDE in the unit cell being identical to the PDE in the perforated domain, except for the viscosity and density constants which are obviously limited. With that, it is deduced that:

$\begin{array}{c}\frac{\mu M_0}{{\sigma }^2}u=\beta \Delta u\end{array}$ (70)

Making the general homogenized PDE (PDE):

$\begin{array}{c}\rho u\cdot \left(u\cdot \nabla \right)+\nabla p-\left(\mu +\beta \right)\Delta u=f\end{array}$ (71)

Note that for the limit in Navier Stokes found by Allaire (1991a), β→0 for ${\sigma }^2\to +\infty$. Equation (71) can be interpreted as a Navier-Stokes system with an apparent viscosity, which is greater than the actual viscosity of the fluid:

$\begin{array}{c}\rho u\cdot \left(u\cdot \nabla \right)+\nabla p-{\mu }^{\text{*}}\Delta u=f\end{array}$ (72)

This apparent viscosity only exists in the numerical problem because the grains do not reach the limit of zero and therefore $M_s$ starts to have a drag and friction effect large enough in the solution to be neglected.

To exemplify the procedure and make it less abstract, consider the case of a soil with porosity $\varphi =0.5$ and $a\left(\epsilon \right)=e^{-\frac{0.8}{{\epsilon }^{2.1}}}$, that is, with a limit at Navier-Stokes. The direct numerical simulation managed to reach the value of $\epsilon =0.5m$ . Thus, the value of $a\left(\epsilon \right)=a\left(0.5\right)=0.033$.

The unit cell must be constructed for $\varphi =0.5$ and $M_s^{unit}=\frac{2}{0.033}=60.60$. The solution of the system of non-linear algebraic equations to define the unit cell configuration resulted in n = 5.267 and r = 0.174. The value of $M_0$ calculated in the unit cell was 11931.06. Furthermore, $\sigma =0.99$ was found, so:

$\begin{array}{c}{\mu }^{\text{*}}=\frac{\mu M_0}{{\sigma }^2}=\frac{\mathrm{0,001}\times 11931.06}{{0.99}^2}=11.93\end{array}$ (73)

Using 11.93 as apparent viscosity and calculating the average velocity within the domain, with a hydraulic gradient ranging from $1\times {10}^{-15}Pa/m²$ to $1\times {10}^{-2}Pa/m²$, a value of average permeability coefficient of 0.007892. Direct numerical simulations (DNS) showed a permeability coefficient tending to a limit at 0.01. Convergence is graphically represented in Figure 13.

Figure 13
DNS convergence for Navier-Stokes with calculated limit.

5.1 Brinkman's limit: a gray area between Darcy and Navier Stokes

So far only numerical convergence for Darcy and Navier-Stokes has been presented. These two limits are the two extremes to which the Navier Stokes equations can converge when treated by the homogenization process by the well-known energy method. However, as initially raised by Allaire (1991a) and studied more closely by Cionarescu & Murat (1982), there is a limit for Brinkman when $\sigma \to 0<C_0<+\infty$. How then to identify this limit in the numerical approach proposed here?

Sticking directly to Allaire's theory (1991a), the only way to achieve convergence for Brinkman would be with a function of the kind:

$\begin{array}{c}a\left(\epsilon \right)=e^{-\frac{\alpha }{{\epsilon }^2}}\end{array}$ (74)

Where α can be any positive real, but the exponent of ε must necessarily be 2. This is due to the fact that σ(ε) is written in the form of a natural logarithm raised to the 1/2 power. When trying to bring this concept to a more applied perspective, only soils whose grain size variation varied exactly with expression (74) would be treated with Brinkman's formulation.

The hypothesis that arises is that for some soils, depending on the size of the representative volume, a function a(ε) equal to (74) or with exponents close to 2, would generate approximately the same value of ${\sigma }^2$. This would suggest that soils with theoretical limits in Navier Stokes could be represented by Brinkman. What is suggested is that the Brinkman limit, in the numerical approach, is just a transitional regime between Navier Stokes and Darcy, acting as an overlap between these two limits.

6. Comparison with laboratory results

It should not be forgotten that this work is a paper in geotechnical engineering. Theoretical results are good, but results that match reality are even better. In this sense, this section will be devoted to the construction of a simple mathematical model for predicting hydraulic conductivity, based on the results of Allaire's theory of homogenization, using parameters known from geotechnical engineering.

6.1 A hydraulic conductivity model based on D₅₀, C_u and C_c

Let the definitions of $C_c$, $C_u$ and σ(ε) be given in (45.a), (45.b) and (46). The model that is proposed from here on will depend only on these parameters and on a characteristic size of the soil particles, represented by the D₅₀.

It is intended to predict the hydraulic conductivity of a homogeneous and isotropic soil from the following sequence of steps:

• From the granulometric curve of the soil in question, identify the values of D₁₀, D₃₀, D₅₀ and D₆₀ so that the values of C_u and C_c can be calculated.

• The typical grain size decay curve as a function of the size of the elemental representative element will be given by (46), that is:

  $\begin{array}{c}a\left(\epsilon \right)=e^{-\frac{C_u}{{\epsilon }^{C_c}}}\end{array}$(75)

• The value of ε that represents the chosen sample will be arbitrated as six times the value of D₅₀, that is, ${\epsilon }_{ref}=6D_{50}$. This approach is not ad-hoc. The work of Hattamleh et. al (2009) points out that this relationship is feasible for granular soils.

• The value of σ(ε) will be adjusted by a function A that depends only on the ratio of $C_u/C_c$. Otherwise, we have:

  $\begin{array}{c}\sigma {\left(\epsilon \right)}^{\text{*}}=A\left(\frac{C_u}{C_c}\right)\sigma \left(\epsilon \right)\end{array}$(76)

• Knowing the porosity of the soil and calculating the unitary specific surface as:

  $\begin{array}{c}{M}_s^{unit}=\frac{2\pi \text{r}}{\pi {\text{r}}^2}=\frac{2}{r}=\frac{4}{D_{50}}\end{array}$(77)

• The nonlinear system of equations proposed in (60.a) and (60.b) is solved to find the radius and the number of grains used in the multigrain unit cell.

• In the multigrain unit cell, the problem $\left(P_{Cel}-Num\right)$ is solved with $w_0=1$ and the value of $M_0$ given by (67) is calculated.

• The value of the average hydraulic conductivity of the soil, directly in cm/s, is given by:

  $\begin{array}{c}K=\frac{\sigma {\left(\epsilon \right)}^{\text{*}}{}^2}{\mu M_0}\end{array}$(78)

Being μ the viscosity of the fluid that saturates the soil, in Pa⋅s. Remember here that the relationship between permeability and hydraulic conductivity is given by $K=k\rho g/\mu$ where K, k, μ, ρ and g are respectively the hydraulic conductivity, intrinsic permeability, fluid viscosity, fluid density and the gravitational constant.

Note that this procedure provides an average hydraulic conductivity value, i.e. independent of the limiting equation. If one wants to know the point of non-linearity, one must simulate the problem (PDE) assuming ${\mu }^{\text{*}}=\frac{\mu M_0}{{\sigma }^{\text{*}}{}^2}$ as increased viscosity for Navier Stokes or $M=\frac{\mu M_0}{{\sigma }^{\text{*}}{}^2}$ for the drag (damping) term in the Brinkman limit.

The equation $A\left(\frac{C_u}{C_c}\right)$ was adjusted for 36 different permeability tests (as shown in Table 1), with soils of different porosities and granulometric curves. The adjustment was made by checking for which value of A, the value of K corresponded to that found in the laboratory. Then, graphs of A were drawn in relation to the specific surface of the grains, porosity and average grain size. None of these relationships offered a fit result with $R^2$ close enough to 1 to establish a correlation. However, when comparing A with $C_u/C_c$, a power law curve with $R^2$ above 0.46 was obtained, even with tests whose values were drastically different from those expected for the material tested. The function found for A was:

$\begin{array}{c}A\left(\frac{C_u}{C_c}\right)=0.0428{\left(\frac{C_u}{C_c}\right)}^{-1.754}\end{array}$ (79)

The relationship between parameter A and the value of the ratio between C_u and C_c is graphically represented in Figure 14. Note that the value of A tends to a vertical asymptote for C_u/C_c →3/2. For the C_u/C_c range between 2 and 4 with permeability values lower than expected, the tested soils had a fines content that may have affected water transport. Soils with high fines content do not have a predictable hydraulic behavior by the model.

Table 1
Laboratory tests used for calibration of the model.

Reference	e	ф	D50 (mm)	Cc	Cu	A	K-Lab (cm/s)
Yang et al. (2016)	0.525	0.34426	0.5	0.75	3	0.0055	1.84e-3
	0.525	0.34426	1	1.03	5.6	0.002	3.83e-3
	0.525	0.34426	2	1.45	11	0.0015	3.68e-3
	0.559	0.35856	2	0.75	5.6	0.0017	1.23e-3
	0.559	0.35856	0.5	1.03	11	0.001	3.35e-3
	0.559	0.35856	1	1.45	3	0.0011	8.99e-3
	0.572	0.36387	1	0.75	11	0.002	1.36e-3
	0.572	0.36387	2	1.03	3	0.0006	3.18e-2
	0.572	0.36387	0.5	1.45	5.6	0.0002	9.44e-4
Akbulut et al. (2014)	0.664	0.39904	1.27	2.14	7.41	0.001	1.89e-1
	0.796	0.44321	1.63	0.89	1.35	0.075	1.26
	0.892	0.47146	0.84	0.93	1.41	0.05	4.18e-1
	0.641	0.39062	1.27	2.14	7.41	0.007	1.73e-1
	0.788	0.44072	1.63	0.89	1.35	0.073	1.23
	0.843	0.45741	0.84	0.93	1.41	0.05	3.75e-1
	0.495	0.3311	1.27	2.14	7.41	0.001	1.39e-1
	0.645	0.3921	1.63	0.89	1.35	0.07	1.07
	0.707	0.41418	0.84	0.93	1.41	0.05	3.96e-1
	0.966	0.49135	1.27	2.14	7.41	0.001	2.69e-1
	1,059	0.51433	1.63	0.89	1.35	0.09	1.46
	1,153	0.53553	0.84	0.93	1.41	0.07	5.97e-1
Calabar & Akbulut (2016)	0.86	0.46237	3.1	0.93	1.5	0.03	1.61
	0.82	0.45055	1.63	0.89	1.35	0.036	0.35
	0.77	0.43503	0.84	0.93	1.41	0.034	0.13
	0.72	0.4186	0.5	1	1.19	0.02	0.05
	0.6	0.375	0.35	0.98	1.16	0.018	0.03
	0.87	0.46524	1	0.73	4.39	0.011	0.03
	0.82	0.45055	0.69	1.01	3.63	0.006	0.02
	0.79	0.44134	0.5	1.44	3.47	0.011	0.02
	0.7	0.41176	0.39	1.74	3.27	0.0004	0.02
	0.63	0.3865	0.32	1.29	3.09	0.0015	0.01
	0.56	0.35897	0.16	0.92	2.07	0.012	0.01
	0.74	0.42529	2.2	0.98	2.6	0.015	0.26
	0.67	0.4012	1.18	0.88	1.93	0.022	0.12
	0.52	0.34211	2	1.22	3.28	0.0083	0.14
	0.7	0.41176	0.99	0.8	2.41	0.025	0.06
	0.55	0.35484	0.7	0.93	1.71	0.016	0.04

Figure 14
Relation between A and C_u/C_c.

Figure 15 shows two different data sets. The blue dots are the hydraulic conductivity values measured in the laboratory for 36 different tests. These same data were used to feed the model. The orange dots represent the prediction of these permeabilities by the model itself. The power law curve fit for both cases (prediction and measured) highlights the similar behavior between the laboratory data and those predicted by the proposed model. Note also that for high values of C_u/C_c we have the lowest model errors. This means that for soils with a tendency to converge to Darcy, the model is more accurate.

Figure 15
Comparison of laboratory results and model prediction for K.

Figure 16 shows the results of direct numerical simulations of the Navier-Stokes equations with varying viscosity, thus simulating the various scenarios of increased viscosity. It is noticed that non-linearities appear quickly for low viscosity values and that, for soils with increased viscosity above 100, the relationship between the hydraulic gradient and the flow velocity is perfectly linear for gradients up to 1000 Pa/m.

Figure 16
Loss of linearity in function of the hydraulic gradient for different augmented viscosities.

It was identified in Figure 17 for which gradient values the non-linearity begins to be observed. Relating these values to their respective C_u/C_c values identifies a power law behavior with $R^2=0.9816$. The same relationship occurs between the critical gradient and the hydraulic conductivity of the soil. This relationship is represented in Figure 18.

Figure 17
Relation between the ratio C_u/C_c and the critical hydraulic gradient.

Figure 18
Relation between the estimated hydraulic conductivity and the critical hydraulic gradient.

Using the formulation presented in Bear (1985) for the Reynolds number in porous media, it was calculated for which Reynolds values the linearity of Darcy's law ceases to apply. Take:

$\begin{array}{c}{R}_e=\frac{u{D}_{50}}{\frac{\mu }{\rho }}=\frac{u{D}_{50}}{\nu }\end{array}$ (80)

Reynolds values 1 ≤ R_e ≤ 10 indicate transition from linear to non-linear regime, i.e., Darcy's law would no longer apply. Comparing this limit interval range with those found in the numerical simulations and listed in Tables 2 and 3, it can be seen that these Reynolds values are those calculated for the last limit in Navier-Stokes and the limit in Brinkman. In the numerical experiment, soils that lost the linearity of the relationship between hydraulic gradient and flow velocity had a relationship 0.26 < C_u/C_c < 0.52.

Table 2
Relation between augmented viscosity and the PDE limit.

u*	M0	Cc	a(ε) (mm)	ε (mm)	A	σ2	σ2 - fit	Limit model
0.001	12848.4	27.3	1	6	14.14	912.12	908.6	NSI
0.01	12848.4	14.2	1	6	4.49	289.82	285.95	NSI
0.1	12848.4	7.3	0.99	6	1.39	90.21	91.86	NSI
1	12848.4	3.8	0.99	6	0.44	28.72	28.87	NSI
10	12848.4	1.95	0.97	6	0.14	9.05	9.3	Brinkman
100	12848.4	1	0.84	6	0.04	3.01	3	Darcy
1000	12848.4	0.5	0.66	6	0.012	1	1.01	Darcy
10000	12848.4	0.25	0.53	6	0.003	0.33	0.3415	Darcy

Table 3
Critical Reynolds number in function of the ratio C_u/C_c.

Cc	Cu	Cu/Cc	i-critical	μ*	K (cm/s)	u-critical (m/s)	Re	Limit model
27.3	1	0.03663	1.00E-05	0.001	70.71	7.00E-07	7.07E-04	Navier-Stokes
14.2	1	0.07042	1.00E-04	0.01	22.25	2.23E-06	2.23E-03	Navier-Stokes
7.3	1	0.13699	1.00E-03	0.1	7.149	7.15E-06	7.15E-03	Navier-Stokes
3.8	1	0.26316	1	1	2.247	2.25E-03	2.23E+00	Navier-Stokes
1.95	1	0.51282	1.00E+02	10	0.724	7.24E-02	7.24E+01	Brinkman
1	1	1	1.00E+04	100	0.233	2.34E+00	2.34E+03	Darcy
0.5	1	2	1.00E+06	1000	0.078	7.88E+01	7.88E+04	Darcy

7. Conclusion

A numerical homogenization process was proposed, based on the finite element method, in its standard formulation (Bubnov-Galerkin). Homogenization in a limited domain had its existence driven by the need to impose conditions on the boundaries of a finite element space. This process consists of the direct numerical simulation (DNS) of a multigrain unit cell, which takes into account two main variables that affect soil permeability: porosity and the specific surface of the grains. From the resolution of a system of algebraic equations, it is determined how many grains must be placed symmetrically in the unit cell and what is the radius of these grains. In the unit cell, a horizontal flow is simulated with specific boundary conditions justified throughout the thesis. With the calculated velocity field, the tensor M is determined.

In Allaire's theory, the M tensor is divided by a factor ${\sigma }^2$ which is a way of measuring the ratio between grain size and domain size. The value of σ depends on the grain decay curve and can have its limit going to +∞ ,0 or 0 < C₀ < +∞. In the first case, the PDE would converge to Navier-Stokes, in the second to Darcy and in the third to Brinkman. Since the grains can become very small, but not zero, the specific surface in numerical simulations assumes large, non-negligible values. What was proposed in this work was:

• Modify the measure σ(ε) by a pre-factor $A\left(\frac{C_u}{C_C}\right)$.

• Define a size for the representative elementary volume from the average grain size, more specifically $\epsilon =6\times D_{50}$.

• Calculate σ(ε) and divide the value of M found by the value of ${\sigma }^2$.

Depending on the value of C_c, define whether the value of $M/{\sigma }^2$ will be used as increased viscosity for Navier-Stokes, Brinkman drag coefficient or the inverse of the hydraulic conductivity in Darcy.

List of symbols

$a\left(\epsilon \right)$ Grain size decay function

$i$ Hydraulic gradient

$k$ Intrinsic permeability

$r$ Grain radius

$w_{\epsilon }^k$ Test function depending on $\epsilon$

$A$ Correction function depending on the grain size distribution

$C_c$ Curvature coefficient

$C_u$ Uniformity coefficient

$D\left(\Omega \right)$ Space of infinite differentiable functions with compact support

$H_0^k$ Hilbert space of functions k-times differentiable and compact support

$H^k$ Hilbert space of functions k-times differentiable

$K$ Hydraulic conductivity

$L^p$ Lebesgue space with p Hölder coefficient

$M_0$ Measure of the tortuosity based on the velocity gradient

$M_s^{unit}$ Unitary specific surface area

$R_e$ Reynold’s number

$W^{p,q}$ Sobolev space with conjugated p-q Hölder coefficients

$\mu$ Fluid viscosity

${\mu }^{\text{*}}$ Augmented viscosity

$\Delta$ Laplacian operator

$\sigma$ Measure of the ratio between the size of the grains and the size of the domain

$\varphi$ Test function or porosity

$\Omega$ Domain

${\Omega }_{\epsilon }$$\epsilon$-dependent domain

$\partial \Omega$ Boundary of the domain

${ℝ}^N$ N-dimensional space of the reals

${\cdot }_X$ Norm of a X-space

$\nabla$ Nabla operator

$\epsilon$ Size of the representative element volume.

Data availability

The datasets generated analyzed in the course of the current study are available from the corresponding author upon request.

• Discussion open until February 28, 2025

References

• ABNT NBR 6502. (2022). Solos e Rochas – terminologia ABNT – Associação Brasileira de Normas Técnicas, Rio de Janeiro, RJ (in portuguese).
• Akbulut, N., Wisniewski, M., & Calabar, A.F. (2014). Influences of grain shape and size distribution on permeability. In Proceedings of the TC105 ISSMGE International Symposium on Geomechanics from Micro to Macro (pp. 1451-1454). Cambridge.
• Allaire, G. (1990). Homogénéisation des équations de Navier-Stokes [Doctoral thesis]. Université Pierre et Marie Curie (in french).
• Allaire, G. (1991a). Homogenization of the Navier Stokes equations in open sets perforated with tiny holes I. Abstract framework, a volume distribution of holes. Archives of Rational Mechanics, 113, 209-259. https://doi.org/10.1007/BF00375065.
• Allaire, G. (1991b). Homogenization of the Navier Stokes equations with a slip boundary condition. Communications on Pure and Applied Mathematics, 44, 605-641. http://doi.org/10.1002/cpa.3160440602.
• ASTM D2487-17. (2017). Standard Practice for Classification of Soils for Engineering Purposes (Unified Soil Classification System) ASTM International, West Conshohocken, PA. https://doi.org/10.1520/D2487-17.
• Bear, J. (1985). Dynamics of fluids in porous media Dover Publications.
• Brenner, S., Scott, R. (2008). The mathematical theory of finite element methods. In A. Bloch, C. L. Epstein, A. Goriely & L. Greengard (Eds.),. Texts in applied mathematics (pp. 280-281). Springer. https://doi.org/10.1007/978-0-387-75934-0.
• Calabar, A. F., Akbulut, N. (2016). Evaluation of actual and estimated hydraulic conductivity of sands with different gradation and shape. SpringerPlus, 5, 820. https://doi.org/10.1186/s40064-016-2472-2.
• Cionarescu, D., & Murat, F. (1982). Un terme étrange venu d’ailleurs. Nonlinear Partial Differential Equations and Their Applications, II and III. In Collége de France Seminar (pp. 98-138). Paris, France.
• Hattamleh, O., Razavi, M., & Muhunthan, B. (2009). Experimental determination of representative elementary volume of sands using X-ray computed tomography. WIT Transactions on Engineering Sciences, 64, 145-153. http://doi.org/10.2495/MC090141.
• Yang, B., Liu, Y., Wan, F., Yang, T., Feng, J., Zhao, X., & Zheng, D. (2016). Experimental study on influence of particle-size distribution on permeability coefficient of sand. Journal of Southwest Jiaotong University, 51(5), 855-861.

Publication Dates

• Publication in this collection
  08 July 2024
• Date of issue
  2024

History

• Received
  30 Nov 2022
• Accepted
  01 Mar 2024