3.1 Formal solution in a unit cell

With the above at our disposal, the analysis can be directed to the homogeneous part ${\mathcal{V}}_{Mh}$ of the entire domain, i.e. excluding the region $I(\unicode[STIX]{x2202}{\mathcal{V}}_{M})$ near the macroscopic boundary of the system. A development similar to that reported by Whitaker (Reference Whitaker1996) may be followed to reach a macroscopic model involving only the average velocity and pressure. However, following the idea used in Barrere, Gipouloux & Whitaker (Reference Barrere, Gipouloux and Whitaker1992, § 2), a shorter alternative procedure, which basically consists of expressing the formal solution of the pore-scale initial boundary value problem over a periodic unit cell in terms of the driving forces, followed by an averaging step, may be adopted to considerably simplify the development. The use of periodicity is more a convenience than a necessity, since the same procedure can be adopted without this restriction. This approach is consistent with the assumption of separation of length scales in the homogeneous region of the porous medium. With this in mind, the pore-scale problem in (2.1) is rewritten in a periodic unit cell in terms of the dependent variables $\boldsymbol{v}$ and $p$ . Unlike $\boldsymbol{v}$ , $p$ is not a periodic field, but the decomposition provided in (3.6) should be used so that $\tilde{p}$ is periodic, yielding

(3.8a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0,\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0, & \displaystyle\end{eqnarray}$$

(3.8b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{v}}{\unicode[STIX]{x2202}t}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}\right)=-\unicode[STIX]{x1D735}\tilde{p}+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}+\underbrace{(-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})}_{source},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0, & \displaystyle\end{eqnarray}$$

(3.8c ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}=\mathbf{0},\quad \text{at }{\mathcal{A}}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}},~t\geqslant 0, & \displaystyle\end{eqnarray}$$

(3.8d ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}=\underbrace{\boldsymbol{v}_{0}}_{source},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},\quad \text{when }t=0, & \displaystyle\end{eqnarray}$$

(3.8e ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \tilde{p}\rangle ^{\unicode[STIX]{x1D6FD}}=0,~t>0, & \displaystyle\end{eqnarray}$$

(3.8f ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}(\boldsymbol{r}+\boldsymbol{l}_{i})=\boldsymbol{v}(\boldsymbol{r}),\quad \tilde{p}(\boldsymbol{r}+\boldsymbol{l}_{i})=\tilde{p}(\boldsymbol{r}),\quad t>0,\quad i=1,2,3, & \displaystyle\end{eqnarray}$$

${\mathcal{V}}_{\unicode[STIX]{x1D6FD}}$

$A_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$

$\unicode[STIX]{x1D6FD}$

$\boldsymbol{l}_{i}$

The nonlinear character of the above problem makes a formal solution very difficult to obtain and, for this reason, a linearization approach is of interest. Indeed, a solution at time $t$ can be sought assuming that the convective velocity exists and is available at a time $t-\unicode[STIX]{x0394}t$ , for any small enough value of $\unicode[STIX]{x0394}t$ . Under these circumstances, the momentum equation (3.8b ) can be approximated using a zeroth-order Taylor expansion in time, leading to the approximation $\boldsymbol{v}|_{t}\approx \boldsymbol{v}|_{t-\unicode[STIX]{x0394}t}\equiv \boldsymbol{v}_{\unicode[STIX]{x1D6E5}}$ . In this way, equation (3.8b ) takes the form

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{v}}{\unicode[STIX]{x2202}t}+\boldsymbol{v}_{\unicode[STIX]{x1D6E5}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}\right)=-\unicode[STIX]{x1D735}\tilde{p}+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}+\underbrace{(-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})}_{source},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0.\end{eqnarray}$$

This type of approximation is a consistent one and is typical in a numerical approach consisting of a linearization of the convective term that makes use of an explicit form of the convective velocity. As will be shown later, the information about the time step and the rest of the terms of the Taylor expansion in time are ultimately not required in the resulting closure problems. With the momentum balance in the form of (3.9), the initial and boundary value problem is a linear one for which a formal solution can be obtained using an integral equation formulation in terms of Green’s functions as shown by Wood & Valdés-Parada (Reference Wood and Valdés-Parada2013). This solution can be written as

(3.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}=\frac{1}{\unicode[STIX]{x1D707}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\,\ast \boldsymbol{\cdot }\,(-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})+\boldsymbol{m}_{0}\right),\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0, & \displaystyle\end{eqnarray}$$

(3.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{p}=\frac{\unicode[STIX]{x2202}\boldsymbol{d}}{\unicode[STIX]{x2202}t}\,\ast \boldsymbol{\cdot }\,(-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})+n_{0},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0. & \displaystyle\end{eqnarray}$$

$\ast \boldsymbol{\cdot }$

$\unicode[STIX]{x1D73F}_{1}$

$\unicode[STIX]{x1D73F}_{2}$

(3.11) $$\begin{eqnarray}\unicode[STIX]{x1D73F}_{1}\,\ast \boldsymbol{\cdot }\,\unicode[STIX]{x1D73F}_{2}=\int _{t_{0}=0}^{t_{0}=t}\unicode[STIX]{x1D73F}_{1}|_{t-t_{0}}\boldsymbol{\cdot }\unicode[STIX]{x1D73F}_{2}|_{t_{0}}\,\text{d}t_{0}=\int _{t_{0}=0}^{t_{0}=t}\unicode[STIX]{x1D73F}_{1}|_{t_{0}}\boldsymbol{\cdot }\unicode[STIX]{x1D73F}_{2}|_{t-t_{0}}\,\text{d}t_{0}.\end{eqnarray}$$

While writing the representations in (3.10), it is meant that $\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}/\unicode[STIX]{x2202}t$ (respectively $\unicode[STIX]{x2202}\boldsymbol{d}/\unicode[STIX]{x2202}t$ ) is the closure variable that maps $(-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})$ onto $\boldsymbol{v}$ (respectively $\tilde{p}$ ) while $\boldsymbol{m}_{0}$ (respectively $n_{0}$ ) is the closure variable that maps $\boldsymbol{v}_{0}$ onto $\boldsymbol{v}$ (respectively $\tilde{p}$ ). At this point, it is worth mentioning that the initial conditions for $\unicode[STIX]{x1D63F}$ and $\boldsymbol{d}$ yield unique solutions that are driven by the sources, as will be provided later in the derivations.

It must be mentioned that the formal solution given in (3.10) does not correspond, except when $\boldsymbol{v}_{0}=\mathbf{0}$ , to the one reported by Lions (Reference Lions1981) (see (5.20) therein), under creeping flow conditions, and by Mikelić (Reference Mikelić1994) (see equation (P) therein for the creeping regime solution and the ‘Proof of Theorem 1.4’ in § 2.4 for the inertial case solution in the Laplace domain). The difference lies in the fact that these authors considered the initial condition to be a function only of $\boldsymbol{x}$ as proposed in (5.10) of chap. 2 in Lions (Reference Lions1981). Note that this assumption is physically questionable and it is not retained in the present work just as it was not considered by Allaire (Reference Allaire1992) for the study of unsteady creeping flow in porous media.

The macroscopic mass and momentum conservation equations can now be obtained by applying the superficial averaging operator on (3.8a ) and (3.10a ), respectively. In order to obtain the macroscale mass conservation equation, it is also necessary to make use of the averaging theorem, together with the no-slip boundary condition, leading to

(3.12a ) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\langle \boldsymbol{v}\rangle =0,\quad \text{in }{\mathcal{V}}_{Mh},~t>0.\end{eqnarray}$$

In addition, the macroscale filtration velocity equation is given by

(3.12b ) $$\begin{eqnarray}\langle \boldsymbol{v}\rangle =-\frac{1}{\unicode[STIX]{x1D707}}\left(\frac{\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D63F}\rangle }{\unicode[STIX]{x2202}t}\,\ast \boldsymbol{\cdot }\,(\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}-\unicode[STIX]{x1D70C}\boldsymbol{b})-\langle \boldsymbol{m}_{0}\rangle \right),\quad \text{in }{\mathcal{V}}_{Mh},~t>0.\end{eqnarray}$$

Here, the assumption that $\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}$ can be considered as a constant within ${\mathcal{V}}$ was taken into account. The effective coefficients $\langle \unicode[STIX]{x1D63F}\rangle$ and $\langle \boldsymbol{m}_{0}\rangle$ are determined by solving the closure problems detailed in the following section. Comments about the physics of the upscaled model and its coefficients are also provided below, together with some remarks on the existing related literature.

3.2 Closure problems and macroscopic model

Substitution of the formal solution given in (3.10) into the initial boundary value problem in (3.8) and separating the contributions from the two sources, while maintaining the convolution product with the volume source, leads to the following equations for $\unicode[STIX]{x1D63F}$ and $\boldsymbol{d}$ :

(3.13a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{t_{0}=0}^{t_{0}=t}\left(\unicode[STIX]{x1D735}\cdot \left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right|_{t-t_{0}}\right)\boldsymbol{\cdot }(-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})_{t_{0}}\,\text{d}t_{0}=\mathbf{0},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0,\end{eqnarray}$$

(3.13b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{t_{0}=0}^{t_{0}=t}\left[\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\left(\left.\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t^{2}}\right|_{t-t_{0}}+\boldsymbol{v}_{\unicode[STIX]{x1D6E5}}\cdot \left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right|_{t-t_{0}}\right)+\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D735}\boldsymbol{d}}{\unicode[STIX]{x2202}t}\right|_{t-t_{0}}-\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right|_{t-t_{0}}-\left.\unicode[STIX]{x1D644}\frac{\text{d}{\mathcal{H}}}{\text{d}t}\right|_{t-t_{0}}\right]\nonumber\\ \displaystyle & & \displaystyle \quad \cdot \,(-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})_{t_{0}}\,\text{d}t_{0}=\mathbf{0},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t\geqslant 0,\end{eqnarray}$$

(3.13c ) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{t_{0}=0}^{t_{0}=t}\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right|_{t-t_{0}}\cdot (-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})_{t_{0}}\,\text{d}t_{0}=\mathbf{0},\quad \text{at }{\mathcal{A}}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}},~t\geqslant 0,\end{eqnarray}$$

${\mathcal{H}}$

$\unicode[STIX]{x1D63F}$

$\boldsymbol{d}$

$\langle \unicode[STIX]{x2202}\boldsymbol{d}/\unicode[STIX]{x2202}t\rangle ^{\unicode[STIX]{x1D6FD}}=\mathbf{0}$

In order to satisfy the above equations, valid for any macroscopic forcing $(-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})$ and at any time, $\unicode[STIX]{x1D63F}$ and $\boldsymbol{d}$ must satisfy the following equations, in general (this may be inferred from considering the particular case in which the macroscopic forcing is a constant):

(3.14a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}=\mathbf{0},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0, & \displaystyle\end{eqnarray}$$

(3.14b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\left(\left.\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t^{2}}\right|_{\unicode[STIX]{x1D70F}}+\boldsymbol{v}\boldsymbol{\cdot }\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right|_{\unicode[STIX]{x1D70F}}\right)+\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D735}\boldsymbol{d}}{\unicode[STIX]{x2202}t}\right|_{\unicode[STIX]{x1D70F}}-\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right|_{\unicode[STIX]{x1D70F}}-\unicode[STIX]{x1D644}\left.\frac{\text{d}{\mathcal{H}}}{\text{d}t}\right|_{\unicode[STIX]{x1D70F}}=\mathbf{0}, & \displaystyle \nonumber\\ \displaystyle & \displaystyle \qquad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},\quad 0\leqslant \unicode[STIX]{x1D70F}\leqslant t,~t>0, & \displaystyle\end{eqnarray}$$

(3.14c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}=\mathbf{0},\quad \text{at }{\mathcal{A}}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}},~t>0, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D70F}$

$\unicode[STIX]{x0394}t\rightarrow 0$

$\boldsymbol{v}_{\unicode[STIX]{x1D6E5}}\rightarrow \boldsymbol{v}$

$\unicode[STIX]{x1D6E5}$

A subsequent time integration step of (3.14) from $\unicode[STIX]{x1D70F}=0$ to $\unicode[STIX]{x1D70F}=t$ yields

(3.15a ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}\right]_{\unicode[STIX]{x1D70F}=0}^{\unicode[STIX]{x1D70F}=t}=\mathbf{0},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0, & \displaystyle\end{eqnarray}$$

(3.15b ) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\left(\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right|_{\unicode[STIX]{x1D70F}}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}|_{\unicode[STIX]{x1D70F}}\right)+\unicode[STIX]{x1D735}\boldsymbol{d}|_{\unicode[STIX]{x1D70F}}-\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63F}|_{\unicode[STIX]{x1D70F}}-\boldsymbol{I}{\mathcal{H}}|_{\unicode[STIX]{x1D70F}}\right]_{\unicode[STIX]{x1D70F}=0}^{\unicode[STIX]{x1D70F}=t}=\mathbf{0},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0,\qquad \quad & \displaystyle\end{eqnarray}$$

(3.15c ) $$\begin{eqnarray}\displaystyle & [\unicode[STIX]{x1D63F}]_{\unicode[STIX]{x1D70F}=0}^{\unicode[STIX]{x1D70F}=t}=\mathbf{0},\quad \text{at }{\mathcal{A}}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}},~t>0. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D63F}$

$\boldsymbol{d}$

As a consequence, equations (3.15) give rise to the following initial and boundary value problem for the closure variables $\unicode[STIX]{x1D63F}$ and $\boldsymbol{d}$ :

Problem I

(3.16a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}=\mathbf{0},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0, & \displaystyle\end{eqnarray}$$

(3.16b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}\right)=-\unicode[STIX]{x1D735}\boldsymbol{d}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63F}+\unicode[STIX]{x1D644},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0, & \displaystyle\end{eqnarray}$$

(3.16c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}=\mathbf{0},\quad \text{at }{\mathcal{A}}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}},~t>0, & \displaystyle\end{eqnarray}$$

(3.16d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}=\mathbf{0},\quad \text{when }t=0, & \displaystyle\end{eqnarray}$$

(3.16e ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{d}\rangle ^{\unicode[STIX]{x1D6FD}}=\mathbf{0},\quad t>0, & \displaystyle\end{eqnarray}$$

(3.16f ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}(\boldsymbol{r}+\boldsymbol{l}_{i})=\unicode[STIX]{x1D63F}(\boldsymbol{r}),\quad \boldsymbol{d}(\boldsymbol{r}+\boldsymbol{l}_{i})=\boldsymbol{d}(\boldsymbol{r}),\quad t>0,\quad i=1,2,3. & \displaystyle\end{eqnarray}$$

Note that in (3.16b ), $\unicode[STIX]{x1D63F}$ , $\boldsymbol{d}$ and $\boldsymbol{v}$ are all evaluated at the same time $t$ , which improves the exactness of the solution. Practically, this is possible because the pore-scale flow problem in (3.8) (the solution of which provides the field of $\boldsymbol{v}$ at time $t$ ) can be solved independently from the closure problem on $\unicode[STIX]{x1D63F}$ and $\boldsymbol{d}$ . On the basis of this closure problem, it is readily deduced that the contribution from the remaining source (i.e. the initial velocity) leads to the following problem for $\boldsymbol{m}_{0}$ and $n_{0}$ :

Problem II

(3.17a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{m}_{0}=0,\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0, & \displaystyle\end{eqnarray}$$

(3.17b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{m}_{0}}{\unicode[STIX]{x2202}t}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{m}_{0}\right)=-\unicode[STIX]{x1D735}n_{0}+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{m}_{0},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0, & \displaystyle\end{eqnarray}$$

(3.17c ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{m}_{0}=\mathbf{0},\quad \text{at }{\mathcal{A}}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}},~t>0, & \displaystyle\end{eqnarray}$$

(3.17d ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{m}_{0}=\unicode[STIX]{x1D707}\boldsymbol{v}_{0},\quad \text{when }t=0, & \displaystyle\end{eqnarray}$$

(3.17e ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle n_{0}\rangle ^{\unicode[STIX]{x1D6FD}}=0,\quad t>0, & \displaystyle\end{eqnarray}$$

(3.17f ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{m}_{0}(\boldsymbol{r}+\boldsymbol{l}_{i})=\boldsymbol{m}_{0}(\boldsymbol{r}),\quad n_{0}(\boldsymbol{r}+\boldsymbol{l}_{i})=n_{0}(\boldsymbol{r}),\quad t>0,\quad i=1,2,3.\qquad \quad & \displaystyle\end{eqnarray}$$

$\ell _{\unicode[STIX]{x1D6FD}}/L$

(3.18a ) $$\begin{eqnarray}\displaystyle & \hspace{-7.0pt}\displaystyle \frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D707}}\!\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\,\ast \boldsymbol{\cdot }\,(-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})+\boldsymbol{m}_{0}\right)\!\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D63F}\right]=-\unicode[STIX]{x1D735}\boldsymbol{d}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D63F}+\unicode[STIX]{x1D644},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.18b ) $$\begin{eqnarray}\displaystyle & \hspace{-7.0pt}\displaystyle \frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\left[\frac{\unicode[STIX]{x2202}\boldsymbol{m}_{0}}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D707}}\!\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\,\ast \boldsymbol{\cdot }\,(-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b})+\boldsymbol{m}_{0}\right)\!\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{m}_{0}\right]=-\unicode[STIX]{x1D735}n_{0}+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{m}_{0},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t>0. & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

$-\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}+\unicode[STIX]{x1D70C}\boldsymbol{b}$

$\boldsymbol{v}_{0}$

$\boldsymbol{v}$

Letting

(3.19a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D643}_{t}=\langle \unicode[STIX]{x1D63F}\rangle , & \displaystyle\end{eqnarray}$$

(3.19b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D736}=\frac{\langle \boldsymbol{m}_{0}\rangle }{\unicode[STIX]{x1D707}}, & \displaystyle\end{eqnarray}$$

(3.20) $$\begin{eqnarray}\langle \boldsymbol{v}\rangle =-\frac{1}{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D643}_{t}}{\unicode[STIX]{x2202}t}\,\ast \boldsymbol{\cdot }\,(\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}-\unicode[STIX]{x1D70C}\boldsymbol{b})+\unicode[STIX]{x1D736},\quad \text{in }{\mathcal{V}}_{Mh},~t>0\end{eqnarray}$$

or, equivalently,

(3.21) $$\begin{eqnarray}\langle \boldsymbol{v}\rangle =-\frac{1}{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D643}_{t}\,\ast \boldsymbol{\cdot }\,(\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}-\unicode[STIX]{x1D70C}\boldsymbol{b}))+\unicode[STIX]{x1D736},\quad \text{in }{\mathcal{V}}_{Mh},~t>0,\end{eqnarray}$$

and this represents one of the major results of this work. It clearly shows the existence of a memory effect expressed by the convolution product that was anticipated by Auriault (Reference Auriault1980) leading to an unsteady macroscopic momentum equation, which does not resemble the heuristic model given by

(3.22) $$\begin{eqnarray}\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}\langle \boldsymbol{v}\rangle ^{\unicode[STIX]{x1D6FD}}}{\unicode[STIX]{x2202}t}=-(\unicode[STIX]{x1D735}\langle p\rangle ^{\unicode[STIX]{x1D6FD}}-\unicode[STIX]{x1D70C}\boldsymbol{b})-\unicode[STIX]{x1D707}\unicode[STIX]{x1D700}\unicode[STIX]{x1D643}^{-1}\boldsymbol{\cdot }\langle \boldsymbol{v}\rangle ^{\unicode[STIX]{x1D6FD}},\quad \text{in }{\mathcal{V}}_{Mh},~t>0,\end{eqnarray}$$

where $\unicode[STIX]{x1D643}$ is the steady apparent permeability tensor introduced by Whitaker (Reference Whitaker1996) for the average model of steady inertial one-phase flow in homogeneous porous media. The two models only match under steady conditions. This can be proved, for instance, by considering that the macroscopic forcing remains constant after a given time. Under such conditions, $\unicode[STIX]{x1D736}\rightarrow \mathbf{0}$ in the long-time limit. Moreover, in this time limit, closure problem I in (3.16f ) conveniently coincides with the one obtained by Whitaker (Reference Whitaker1996) for steady inertial flow in homogeneous porous media and $\unicode[STIX]{x1D643}_{t}\rightarrow \unicode[STIX]{x1D643}$ . Under these circumstances, the final value theorem applied to (3.20) (or (3.21)) indicates that the average model derived above reduces to the steady form of the macroscopic inertial momentum equation (3.22) also reported by Whitaker (Reference Whitaker1996).

In the macroscopic equation (3.20) (or (3.21)), $\unicode[STIX]{x1D643}_{t}$ is homogeneous to a permeability and will be referred to as the dynamic apparent permeability tensor, the apparent character being inherent to its dependence on inertial effects. It should be noted that, except in the creeping flow regime, the effective coefficient $\unicode[STIX]{x1D643}_{t}$ depends on the initial condition and on the macroscopic forcing through the convective inertial term in closure problem I. In addition, $\unicode[STIX]{x1D736}$ , which has the unit of a velocity, only contains a source due to the initial velocity field, $\boldsymbol{v}_{0}$ , and is zero when $\boldsymbol{v}_{0}=\mathbf{0}$ . However, the values of $\unicode[STIX]{x1D736}$ are also driven by the macroscopic forcing by means of the convective term in problem II. It is important to emphasize this feature and make clear that the effective coefficients are not functions of the average velocity, which would result in a misleading macroscale model given in (3.20). Unfortunately, the dependence of the coefficients on the macroscopic forcing and the initial condition is not trivial and, as a consequence, an explicit functionality of the seepage velocity with them is not easily achievable, in general. Under non-inertial conditions, the dependence of the macroscale velocity on the macroscopic forcing is linear, albeit the dependence on the initial flow is still not trivial. For steady, creeping flow, the well-known linear dependence of the velocity on the macroscopic forcing is recovered in the form of Darcy’s law.

The model reported in (3.20) generalizes to inertial flow the result reported by Allaire (Reference Allaire1992) for the creeping regime, which was also studied by Lions (Reference Lions1981) and later by Mikelić (Reference Mikelić1994). When restricted to this particular type of flow, the macroscale momentum equation (3.20) matches that reported in the latter two references only when the initial flow is zero. The discrepancy observed when $\boldsymbol{v}_{0}\neq \mathbf{0}$ lies in the fact that, as mentioned above (§ 3.1), a particular form of the initial condition was considered by Lions (Reference Lions1981) and Mikelić (Reference Mikelić1994) for the problem in the periodic unit cell. This special form of the boundary condition was, however, not retained by Allaire (Reference Allaire1992) nor in the present work. In the particular case of creeping flow and $\boldsymbol{v}_{0}=\mathbf{0}$ , envisaged by Mei & Vernescu (Reference Mei and Vernescu2010), agreement is also found with the result of this reference. In Mikelić (Reference Mikelić1994), both creeping and inertial flows were considered; unfortunately, no local closure problem was derived for the inertial case. In the work by Allaire (Reference Allaire1992), $\text{A}$ therein identifies with $\unicode[STIX]{x2202}\unicode[STIX]{x1D646}_{t}/\unicode[STIX]{x2202}t$ and $\boldsymbol{a}$ with $\unicode[STIX]{x1D736}$ , $\unicode[STIX]{x1D646}_{t}$ being the dynamic permeability equivalent to $\unicode[STIX]{x1D643}_{t}$ under non-inertial flow conditions. It must be emphasized that $\text{A}$ should not be called ‘permeability’ as it is dimensionally incorrect. This terminology is also improperly used in the works by Mikelić (Reference Mikelić1994) and Mei & Vernescu (Reference Mei and Vernescu2010). Without inertia and when $\boldsymbol{v}_{0}=\mathbf{0}$ , the average model derived above also corresponds to the one presented by Auriault et al. (Reference Auriault, Borne and Chambon1985), Sheng & Zhou (Reference Sheng and Zhou1988) or Zhou & Sheng (Reference Zhou and Sheng1989) (see also Sahimi Reference Sahimi2011) and considered by Johnson et al. (Reference Johnson, Koplik and Dashen1987) that was obtained in Fourier space.

Regarding the ancillary closure problems related to the macroscopic models mentioned above, closure problems I and II, under creeping flow conditions, coincide with those reported by Allaire (Reference Allaire1992). However, under these conditions, closure problem I does not correspond to the ancillary problems reported by Mikelić (Reference Mikelić1994) and by Mei & Vernescu (Reference Mei and Vernescu2010). The closure problems in these two references lack compatibility between the initial and boundary conditions, thus leading to a non-regular solution as admitted by Mikelić (Reference Mikelić1994). This discrepancy is only technical and can be solved by an appropriate change of variables in the Laplace domain, leading to a modification of the initial condition, compatible with the interfacial boundary condition. Consistently, in the creeping regime, the Laplace-transformed version of closure problem I is also identical to that given by Sheng & Zhou (Reference Sheng and Zhou1988), Zhou & Sheng (Reference Zhou and Sheng1989), Mikelić (Reference Mikelić1994) and Mei & Vernescu (Reference Mei and Vernescu2010) for their corresponding so-called ‘permeability’.

It should be noted that in Lions & Masmoudi (Reference Lions and Masmoudi2005) an attempt to use the homogenization method to upscale the unsteady Navier–Stokes equations was presented. Unfortunately, the authors only succeeded to use the two-scale convergence method for the case of perfect flow (i.e. the unsteady Euler equation) and attributed their failure to upscale the unsteady incompressible Navier–Stokes equations to the existence of boundary layers. In fact, in their study no upscaled model or closure problems were provided. The drawback arising from boundary layers was also mentioned by Masmoudi (Reference Masmoudi1998) and was later circumvented by the same author for the case of compressible flow (see Masmoudi Reference Masmoudi2002). Nevertheless, the development in this latter reference yields a semi-stationary macroscopic model. The difficulty encountered by these authors arises from the fact that the solution is sought in a weak sense. It is not present in the approach followed in the current analysis.

Further considering the very particular case of the creeping regime for which the initial flow condition is such that $\boldsymbol{v}_{0}$ obeys a Stokes model, it can be proved that the macroscopic model derived above can be formulated in such a way that closure problem I is the only one that needs to be solved. The proof is provided in the Appendix.

Before illustrating the solution of closure problems I and II with some numerical results, it is of interest to analyse the symmetry and positiveness properties of $\unicode[STIX]{x1D643}_{t}$ , and this is the object of the next section.

3.3 Symmetry properties and positiveness of $\unicode[STIX]{x1D643}_{t}$

The symmetry analysis of $\unicode[STIX]{x1D643}_{t}$ is carried out following the approach developed in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2017) and the reader is referred to this article for details of the derivations, in particular § II.A therein. It must be noted that no special assumption is needed on the pore structure within the periodic unit cell representative of the material on which closure problem I is to be solved.

The analysis starts by redirecting attention to (3.14b ), which is considered at any value of $0\leqslant \unicode[STIX]{x1D70F}\leqslant t$ for a given value of $t$ . Pre-multiplication by $(\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}^{\text{T}}/\unicode[STIX]{x2202}t)|_{t-\unicode[STIX]{x1D70F}}$ , together with a subsequent time integration from $\unicode[STIX]{x1D70F}=0$ to $\unicode[STIX]{x1D70F}=t$ and the application of the superficial averaging operator leads to

(3.23) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\langle \unicode[STIX]{x1D648}^{\text{T}}\,\ast \boldsymbol{\cdot }\,\unicode[STIX]{x1D648}\rangle +\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\langle \unicode[STIX]{x1D648}^{\text{T}}\,\ast \boldsymbol{\cdot }\,\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D648}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =-\langle \unicode[STIX]{x1D648}^{\text{T}}\,\ast \boldsymbol{\cdot }\,\unicode[STIX]{x1D735}\boldsymbol{m}\rangle +\langle \unicode[STIX]{x1D648}^{\text{T}}\,\ast \boldsymbol{\cdot }\,\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D648}\rangle +\frac{\text{d}\unicode[STIX]{x1D643}_{t}^{\text{T}}}{\text{d}t},\quad t>0,\end{eqnarray}$$

where, for simplicity of notation, $\unicode[STIX]{x1D648}=\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}/\unicode[STIX]{x2202}t$ and $\boldsymbol{m}=\unicode[STIX]{x2202}\boldsymbol{d}/\unicode[STIX]{x2202}t$ . Note that $\unicode[STIX]{x1D648}=\mathbf{0}$ at $t=0$ , in accordance with (3.16d ), given the initial condition for $\unicode[STIX]{x1D63F}$ .

The first term on the left-hand side of (3.23) is clearly symmetric. As shown in Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2017), the first term on the right-hand side is zero using the solenoidal character of $\unicode[STIX]{x1D648}$ , the no-slip boundary condition, periodicity and the averaging theorem. On the same basis, it can also be proved that $\langle \unicode[STIX]{x1D648}^{\text{T}}\,\ast \boldsymbol{\cdot }\,\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D648}\rangle =-\langle (\unicode[STIX]{x1D735}\unicode[STIX]{x1D648})^{\text{T}3}\,\ast \boldsymbol{ : }\,\unicode[STIX]{x1D735}\unicode[STIX]{x1D648}\rangle$ , where the superscript T3 stands for the transpose that permutes the first and third indices of a third-order tensor and $\boldsymbol{ : }$ is the double dot product in the sense of the nested convention. This last term can be shown to be symmetric. Finally, the second term on the left-hand side of (3.23) is skew-symmetric. The proof of this can be carried out on the basis of (12) and (13) in § II.A in the work by Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2017), when extended to the case where a convolution product is involved.

From the above, the time rate of change of the dynamic apparent permeability tensor can hence be expressed as

(3.24) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D643}_{t}}{\unicode[STIX]{x2202}t} & = & \displaystyle \underbrace{\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}^{\text{T}}}{\unicode[STIX]{x2202}t}\,\ast \boldsymbol{\cdot }\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right\rangle +\left\langle \left(\unicode[STIX]{x1D735}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right)^{\text{T}3}\,\ast \boldsymbol{ : }\,\left(\unicode[STIX]{x1D735}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right)\right\rangle }_{symmetric\,part}\nonumber\\ \displaystyle & & \displaystyle -\underbrace{\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\left\langle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}^{\text{T}}}{\unicode[STIX]{x2202}t}\,\ast \boldsymbol{\cdot }\,\left(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\boldsymbol{v}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}{\unicode[STIX]{x2202}t}\right)\right)\right\rangle }_{skew\text{-}symmetric\,part},\quad t>0,\end{eqnarray}$$

which provides the decomposition of $\unicode[STIX]{x1D643}_{t}$ into its irreducible parts because the operation of time integration does not alter the symmetry properties of a tensor. It should be noted that the skew-symmetric part is only due to the existence of inertial transport, thus extending the result given in § II.A of the work by Lasseux & Valdés-Parada (Reference Lasseux and Valdés-Parada2017) to unsteady conditions. As a corollary, it can be concluded that, in the creeping regime, $\unicode[STIX]{x1D646}_{t}$ is a symmetric tensor at all times (and therefore the intrinsic permeability is also a symmetric tensor).

At this point, the focus should be laid upon the positiveness property of $\unicode[STIX]{x1D643}_{t}$ and, for convenience, the analysis is carried out in the Laplace domain. A starting point of the derivation is the following identity, which holds for any constant but arbitrary vector  $\unicode[STIX]{x1D740}$ :

(3.25) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D740}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D643}}_{t}\boldsymbol{\cdot }\unicode[STIX]{x1D740}=\unicode[STIX]{x1D740}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D643}}_{t}^{\text{T}}\boldsymbol{\cdot }\unicode[STIX]{x1D740}, & & \displaystyle\end{eqnarray}$$

where overlined variables denote the Laplace transform of their temporal counterparts. This implies that

(3.26) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D740}\boldsymbol{\cdot }(\overline{\unicode[STIX]{x1D643}}_{t}+\overline{\unicode[STIX]{x1D643}}_{t}^{\text{T}})\boldsymbol{\cdot }\unicode[STIX]{x1D740}=2\unicode[STIX]{x1D740}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D643}}_{t}\boldsymbol{\cdot }\unicode[STIX]{x1D740}. & & \displaystyle\end{eqnarray}$$

Applying the Laplace transform to (3.24), adding the result to its transpose and pre- and post-multiplying the ensuing expression by $\unicode[STIX]{x1D740}$ while making use of (3.26) yields

(3.27) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D740}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D643}}_{t}\boldsymbol{\cdot }\unicode[STIX]{x1D740}=\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}s^{2}\langle \unicode[STIX]{x1D740}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D63F}}^{\text{T}}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D63F}}\boldsymbol{\cdot }\unicode[STIX]{x1D740}\rangle +s\langle \unicode[STIX]{x1D740}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\overline{\unicode[STIX]{x1D63F}})^{\text{T}3}\,:\,\unicode[STIX]{x1D735}\overline{\unicode[STIX]{x1D63F}}\boldsymbol{\cdot }\unicode[STIX]{x1D740}\rangle . & & \displaystyle\end{eqnarray}$$

Here $s$ denotes the symbolic Laplace variable. The first term on the right-hand side of this expression can be written as

(3.28) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D740}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D63F}}^{\text{T}}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D63F}}\boldsymbol{\cdot }\unicode[STIX]{x1D740}\rangle =\langle (\overline{\unicode[STIX]{x1D63F}}\boldsymbol{\cdot }\unicode[STIX]{x1D740})\boldsymbol{\cdot }(\overline{\unicode[STIX]{x1D63F}}\boldsymbol{\cdot }\unicode[STIX]{x1D740})\rangle =\langle (\overline{\unicode[STIX]{x1D63F}}\boldsymbol{\cdot }\unicode[STIX]{x1D740})^{2}\rangle , & & \displaystyle\end{eqnarray}$$

which is a positive definite quantity. Turning attention to the second average term on the right-hand side of (3.27) and making use of the Gibbs and Einstein notations, one can write

(3.29) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D740}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\overline{\unicode[STIX]{x1D63F}})^{\text{T}3}\,:\,\unicode[STIX]{x1D735}\overline{\unicode[STIX]{x1D63F}}\boldsymbol{\cdot }\unicode[STIX]{x1D740}\rangle & = & \displaystyle \langle \unicode[STIX]{x1D740}_{k}(\unicode[STIX]{x1D735}\overline{\unicode[STIX]{x1D63F}})_{kij}^{\text{T}3}(\unicode[STIX]{x1D735}\overline{\unicode[STIX]{x1D63F}})_{jil}\unicode[STIX]{x1D740}_{l}\rangle \nonumber\\ \displaystyle & = & \displaystyle \langle \unicode[STIX]{x1D740}_{k}(\unicode[STIX]{x1D735}\overline{\unicode[STIX]{x1D63F}})_{jik}(\unicode[STIX]{x1D735}\overline{\unicode[STIX]{x1D63F}})_{jil}\unicode[STIX]{x1D740}_{l}\rangle \nonumber\\ \displaystyle & = & \displaystyle \langle ((\unicode[STIX]{x1D735}\overline{\unicode[STIX]{x1D63F}})_{jik}\unicode[STIX]{x1D740}_{k})^{2}\rangle ,\end{eqnarray}$$

which proves that this term is also positive definite. Since $s$ is positive, it can be concluded that $\overline{\unicode[STIX]{x1D643}}_{t}$ is a positive definite tensor. A corollary of the above is that $\overline{\unicode[STIX]{x1D646}}_{t}$ is a symmetric positive definite tensor that hence admits an inverse. However, the proof provided here does not allow us to conclude that the same applies to $\overline{\unicode[STIX]{x1D643}}_{t}$ .

3.4 Closure problem solution

From the above derivations, it follows that closure problems I and II can be solved using standard unsteady Navier–Stokes solvers. The mathematical structure of these problems indicates that closure variable $\unicode[STIX]{x1D63F}$ is a time-increasing field because the initial condition is homogeneous and the source term in the momentum-like equation is a positive constant. Consequently, the effective coefficient $\unicode[STIX]{x1D643}_{t}$ should be expected to exhibit similar dynamics. On the contrary, in problem II, the only source is the initial condition and therefore both the $\boldsymbol{m}_{0}$ field and coefficient $\unicode[STIX]{x1D736}$ can be expected to be time-decaying.

Before proceeding to the validation of the average model, it is convenient to direct attention to the solution of the local closure problems I and II so that the dynamics of the effective coefficients can be examined. To this end, it is worth expressing the closure problems and the effective coefficients in the following dimensionless forms:

Problem $\text{I}^{\ast }$

(3.30a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D63F}^{\ast }=\mathbf{0},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t^{\ast }>0, & \displaystyle\end{eqnarray}$$

(3.30b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}^{\ast }}{\unicode[STIX]{x2202}t^{\ast }}+Re\,\boldsymbol{v}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D6FB}^{\ast }\unicode[STIX]{x1D63F}^{\ast }=-\unicode[STIX]{x1D6FB}^{\ast }\boldsymbol{d}^{\ast }+\unicode[STIX]{x1D6FB}^{\ast 2}\unicode[STIX]{x1D63F}^{\ast }+\unicode[STIX]{x1D644},\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t^{\ast }>0, & \displaystyle\end{eqnarray}$$

(3.30c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}^{\ast }=\mathbf{0},\quad \text{at }{\mathcal{A}}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}},~t^{\ast }>0, & \displaystyle\end{eqnarray}$$

(3.30d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}^{\ast }=\mathbf{0},\quad \text{when }t^{\ast }=0, & \displaystyle\end{eqnarray}$$

(3.30e ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \boldsymbol{d}^{\ast }\rangle ^{\unicode[STIX]{x1D6FD}}=\mathbf{0},\quad t^{\ast }>0, & \displaystyle\end{eqnarray}$$

(3.30f ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D63F}^{\ast }(\boldsymbol{r}^{\ast }+\boldsymbol{l}_{i}^{\ast })=\unicode[STIX]{x1D63F}^{\ast }(\boldsymbol{r}^{\ast }),\quad \boldsymbol{d}^{\ast }(\boldsymbol{r}^{\ast }+\boldsymbol{l}_{i}^{\ast })=\boldsymbol{d}^{\ast }(\boldsymbol{r}^{\ast }),\quad t^{\ast }>0,\quad i=1,2,3, & \displaystyle\end{eqnarray}$$

(3.31) $$\begin{eqnarray}\unicode[STIX]{x1D643}_{t}^{\ast }=\langle \unicode[STIX]{x1D63F}^{\ast }\rangle .\end{eqnarray}$$

Problem II $^{\ast }$

(3.32a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{\ast }\boldsymbol{\cdot }\boldsymbol{m}_{0}^{\ast }=0,\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t^{\ast }>0, & \displaystyle\end{eqnarray}$$

(3.32b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{m}_{0}^{\ast }}{\unicode[STIX]{x2202}t^{\ast }}+Re\,\boldsymbol{v}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D6FB}^{\ast }\boldsymbol{m}_{0}^{\ast }=-\unicode[STIX]{x1D735}n_{0}^{\ast }+\unicode[STIX]{x1D6FB}^{\ast 2}\boldsymbol{m}_{0}^{\ast },\quad \text{in }{\mathcal{V}}_{\unicode[STIX]{x1D6FD}},~t^{\ast }>0, & \displaystyle\end{eqnarray}$$

(3.32c ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{m}_{0}^{\ast }=\mathbf{0},\quad \text{at }{\mathcal{A}}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}},~t^{\ast }>0, & \displaystyle\end{eqnarray}$$

(3.32d ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{m}_{0}^{\ast }=\boldsymbol{v}_{0}^{\ast },\quad \text{when }t^{\ast }=0, & \displaystyle\end{eqnarray}$$

(3.32e ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle n_{0}^{\ast }\rangle ^{\unicode[STIX]{x1D6FD}}=0,\quad t^{\ast }>0, & \displaystyle\end{eqnarray}$$

(3.32f ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{m}_{0}^{\ast }(\boldsymbol{r}^{\ast }+\boldsymbol{l}_{i}^{\ast })=\boldsymbol{m}_{0}^{\ast }(\boldsymbol{r}^{\ast }),\quad n_{0}^{\ast }(\boldsymbol{r}^{\ast }+\boldsymbol{l}_{i}^{\ast })=n_{0}(\boldsymbol{r}^{\ast }),\quad t^{\ast }>0,\quad i=1,2,3. & \displaystyle\end{eqnarray}$$

(3.33) $$\begin{eqnarray}\unicode[STIX]{x1D736}^{\ast }=\langle \boldsymbol{m}_{0}^{\ast }\rangle .\end{eqnarray}$$

The above problems are written in terms of the following definitions:

(3.34) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}t^{\ast }={\displaystyle \frac{\unicode[STIX]{x1D707}t}{\unicode[STIX]{x1D70C}\ell ^{2}}},\quad \boldsymbol{x}^{\ast }={\displaystyle \frac{\boldsymbol{x}}{\ell }},\quad \boldsymbol{v}^{\ast }={\displaystyle \frac{\boldsymbol{v}}{v_{ref}}},\quad \unicode[STIX]{x1D63F}^{\ast }={\displaystyle \frac{\unicode[STIX]{x1D63F}}{\ell ^{2}}},\quad \boldsymbol{d}^{\ast }={\displaystyle \frac{\boldsymbol{d}}{\ell }},\quad \unicode[STIX]{x1D643}_{t}^{\ast }={\displaystyle \frac{\unicode[STIX]{x1D643}_{t}}{\ell ^{2}}},\\ \boldsymbol{m}_{0}^{\ast }={\displaystyle \frac{\boldsymbol{m}_{0}}{\unicode[STIX]{x1D707}v_{ref}}},\quad n_{0}^{\ast }={\displaystyle \frac{n_{0}\ell }{\unicode[STIX]{x1D707}v_{ref}}},\quad \unicode[STIX]{x1D736}^{\ast }={\displaystyle \frac{\unicode[STIX]{x1D736}}{v_{ref}}},\quad Re={\displaystyle \frac{\unicode[STIX]{x1D70C}v_{ref}\ell }{\unicode[STIX]{x1D707}}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\ell$ is the size of the geometrical periodic unit cell (see figure 2). As expected, the closure problems are dependent on $\boldsymbol{v}^{\ast }$ . Consequently, the values of the effective coefficients are sensitive to the different flow conditions considered in the pore-scale model. Here and in the rest of this work, it is assumed that the macroscopic driving force is in the horizontal $x$ -axis direction. In the following section, several case studies are considered for validation with DNS, for which the predictions of the effective coefficients corresponding only to case study I (see § 4.1) are shown here for the purpose of illustration. All the simulations presented in this work were performed using the commercial finite element software COMSOL Multiphysics 5.2. The direct MUMPS solver included in the program was used and standard (spatial and temporal) meshing convergence analyses were carried out in order to ensure that results are independent of these numerical degrees of freedom. In accordance with the case study of § 4.1, the dimensionless macroscopic pressure gradient experiences a step change at $t^{\ast }=0$ from 0.1 to 1.0 as expressed in (4.5), while $Re$ , ranging from $10^{3}$ up to $10^{6}$ , is maintained the same before and after the change of the macroscopic pressure gradient.

Figure 3.

Dynamics of the $xx$ -component of the apparent permeability tensor (a–c) and of the $x$ -component of vector $\unicode[STIX]{x1D736}$  (d) for the flow problem described in § 4.1. The predictions result from solving the closure problems I and II in the unit cell depicted in figure 2 taking: (a)  $\unicode[STIX]{x1D700}=0.4$ and (b)  $\unicode[STIX]{x1D700}=0.6$ . In (c), results are presented after normalization with the steady value of the apparent permeability. In (d), results are normalized with the initial value of $\unicode[STIX]{x1D6FC}_{x}$ .

In figure 3, predictions of the dimensionless $xx$ -component of $\unicode[STIX]{x1D643}_{t}$ (i.e. $H_{txx}^{\ast }$ ) and of the $x$ -component of $\unicode[STIX]{x1D736}$ (i.e. $\unicode[STIX]{x1D6FC}_{x}$ ), resulting from the closure problem solution in the unit cell depicted in figure 2, are reported for two porosity values ( $\unicode[STIX]{x1D700}=0.4$ and $0.6$ ). Values of the Reynolds number were kept smaller than the critical one characteristic of the first Hopf bifurcation, which, for the structure under consideration, is identified to be ${\sim}10^{6}$ for $\unicode[STIX]{x1D700}=0.4$ and ${\sim}10^{5}$ for $\unicode[STIX]{x1D700}=0.6$ (see Agnaou et al. Reference Agnaou, Lasseux and Ahmadi2016, figure 10). Consequently, the solution remains time-independent after steady state is reached. Results on $H_{txx}^{\ast }$ , shown in figure 3(a,b), indicate that the influence of inertial transport is experienced not only at early times (i.e. $t^{\ast }<10^{-4}$ ), but also during intermediate times (i.e. $t^{\ast }>10^{-3}$ ) and until steady state is reached (i.e. $H_{txx}\rightarrow H_{xx}$ ). By comparing the dynamics of the apparent permeability shown in these two panels, it is clear that porosity plays quite a significant role, because the time at which $H_{txx}$ reaches $H_{xx}$ , for $\unicode[STIX]{x1D700}=0.6$ , is significantly larger than the one for $\unicode[STIX]{x1D700}=0.4$ . In this latter case, it is worth noticing that the curves of $H_{txx}$ remain almost unchanged over the whole range of time for $Re\leqslant 10^{4}$ , whereas in figure 3(b) the curves for $Re=10^{3}$ and $Re=10^{4}$ can be clearly differentiated. These observations are consistent with those reported by Lasseux, Abbasian-Arani & Ahmadi (Reference Lasseux, Abbasian-Arani and Ahmadi2011) under steady conditions.

The shape of the curves shown in figure 3(a,b) suggests a subsequent normalization by the steady-state value, $H_{xx}$ , of $H_{txx}$ , so that the dependence on the Reynolds number is no longer present. This is indeed the case as reported in figure 3(c), which follows from the ideas proposed by Sheng & Zhou (Reference Sheng and Zhou1988). The normalized dynamic apparent permeability may then be represented by a linear combination of exponential linear functions of time. The same type of normalization can be applied to $\unicode[STIX]{x1D6FC}_{x}^{\ast }$ , but in this case with respect to the initial value, since it is the maximum one. Results are presented in figure 3(d), showing that this normalization yields a master curve for the dynamics of $\unicode[STIX]{x1D736}$ for Reynolds numbers as high as $10^{6}$ . The dynamics of $\unicode[STIX]{x1D736}$ exhibits a temporal dependence, which can be represented by a Boltzmann-type function. This contrasts with the purely exponential decay suggested by Mikelić (Reference Mikelić1994) for periodic structures.

The master curves for the effective coefficients are distinct for different porosities, and, for the simple geometry considered here, it appears that the steady state is reached faster, as the solid phase occupies a larger fraction of the unit cell, which was also the case for $H_{txx}$ . These results evidence that the effective coefficients are bound functions of time and are sensitive to the topology of the unit cell and to the flow conditions, in general.