A pure Stokes approach for coupling fluid flow with
porous media flow

Modesar Shakoor*¹ and Chung Hae Park¹

¹IMT Nord Europe, Institut Mines-Télécom, Univ. Lille, Centre for
Materials and Processes, F-59000 Lille, France

January 10, 2024

<figure><img src="image_1.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License.

**Abstract**

Most numerical approaches for coupling fluid flow with porous media flow rely either on Stokes equations in the fluid part of the domain and Darcy's law in the porous part, or on Brinkman's equation. In both cases, difficulties arise at the boundary between the two parts because the equations used in the porous part are not of Stokes type. In this paper, an alternative to Darcy's law is proposed for modeling flows in porous media. This alternative relies on equations of Stokes type where the permeability tensor is replaced by force and stress derivative tensors. Numerical procedures are presented to compute these tensors from simulations at pore scale. Simulations in domains containing both fluid and porous parts are finally conducted simply assuming continuity of velocity and pressure and hence without imposing any condition at the boundary between the two parts. Results

<footnote>*Corresponding author: modesar.shakoor@imt-nord-europe.fr</footnote>
show that the proposed method is accurate and hence a promising alternative to
Darcy's law for problems involving both fluid and porous parts.

Keywords: porous media, Stokes, Darcy, Brinkman, Beavers-Joseph-Saffman

**Article highlights:**

* Alternative to Darcy's law relying on equations of Stokes type for modeling flows in porous media

* Permeability tensor replaced by force and stress derivatives tensors computed from simulations at pore scale

* Flows in domains containing both fluid and porous parts modeled without any boundary conditions between the two parts

# 1 Introduction

Flows in porous media are of great interest for a wide range of applications such as oil engineering, civil engineering and composites engineering. Although Stokes equations can model these flows at the pore scale, they are not convenient for large simulation domains where the pore characteristic length is very small compared to the domain size. Under such conditions, indeed, simulations discretizing and solving simultaneously both the coarse scale of the domain and the fine scale of the pores become intractable from a computational cost point of view. Darcy's law, fortunately, has been proposed to circumvent this difficulty. This law has first been established from experimental measurements (Darcy, 1856), and later on proven from a mathematical point of view in the steady case (Whitaker, 1966; Neuman, 1977; Whitaker, 1986). The main advantage of this law is that an effective property of the porous medium called permeability tensor can be computed beforehand in a pre-processing step, and then be used for coarse scale simulations where pores are not discretized.

For single-phase flows in domains containing both a porous part and a purely fluid part (i.e., without pores), it seems straightforward to rely on Stokes equations in the fluid part and Darcy's law where pores are present. In practice, this is quite challenging due to intrinsic mathematical and physical differences between both models:

* Velocities have different meanings in both parts (Marco Discacciati and Alfio Quar-
Quarteroni, 2009), as the Darcy velocity $v^D$ is related to the fluid velocity $v$ through the relation $v^D = \phi v$, with $\phi$ the porosity (i.e., pore volume fraction).

* Darcy's law involves no differential operator for the velocity while Stokes equations involve a differential operator with second order derivatives of the velocity.

* It is common to eliminate completely the velocity from Darcy's law and solve with only the pressure as unknown, while this is not the case for Stokes equations. More generally, numerical methods for efficiently solving Darcy's law and Stokes equations might be different (Chidyagwai and Rivière, 2010).

These difficulties are usually dealt with using specific boundary conditions at the interface between the fluid and porous parts (Jäger et al., 2001; Layton et al., 2002; Urquiza et al., 2008; Marco Discacciati and Alfio Quarteroni, 2009; Münzenmaier and Starke, 2011; Pacquaut et al., 2012; Eggenweiler and Rybak, 2020):

* Continuity of the normal fluid velocity to ensure mass conservation across the interface.

* Continuity of the fluid normal stress to ensure balance of normal forces across the interface.

* Some condition on the tangential fluid velocity such as the Beavers-Joseph or Beavers-Joseph-Saffman conditions, or even a simple no-slip condition.

Beavers-Joseph and Beavers-Joseph-Saffman conditions have been the subject of many research works in the last decades. As shown in Ref. (Auriault, 2010), they can be used only in restrictive conditions where the pore characteristic length is not too small as compared to the domain size, and they involve an empirical coefficient that is “dedicated to the particular setup at stake”. A numerical proof can be found in Ref. (Eggenweiler and Rybak, 2020). For the cases with a very small pore characteristic length, a boundary condition involving a permeability tensor specific to the interface has been proposed in Ref. (Auriault, 2010). Similarly, the boundary condition proposed in Ref. (Marušić-Paloka and Pažanin, 2022) features a so-called “permeability tensor of the porous wall”. An “interface permeability tensor” also appears in the higher-order boundary conditions
introduced in Ref. (Sudhakar et al., 2021), but it is accompanied by other terms that
can all be pre-computed by solving fine scale problems.

These boundary conditions, with their empirical or pre-computed terms, should be in-
troduced independently from the numerical methods used to solve and couple the fluid and
porous media flows. Unified or monolithic approaches based on least squares (Münzen-
maier and Starke, 2011), the level-set method (Pacquaut et al., 2012) or Brinkman's
equation (Payne and Straughan, 1998; Auriault, 2009) may be used, as well as other
approaches where the two parts are decoupled and discretized separately.

It is important to point out that the interface between the fluid and porous parts is an
artificial mathematical boundary and it can be given any arbitrary thickness as long as
this thickness is small compared to the pore characteristic length (Marciniak-Czochra and
Mikelić, 2012). As an alternative, the so-called one domain approach consists in modeling
this boundary as a new intermediary porous medium whose thickness is not arbitrary but
determined iteratively by the numerical method (Valdés-Parada and Lasseux, 2021a,b).

In a recent work (Blanco et al., 2017), an alternative approach for homogenizing flows
in porous media has been proposed. First, the velocity in the porous medium is averaged
only over the pores, and therefore carries the same physical meaning of fluid velocity
in both fluid and porous parts. Second, the coarse scale model is of Stokes type, and
hence matches the model used in the fluid medium. This model has been extended to the
unsteady case and implemented with two-way coupling at each time increment between
the two scales in Ref. (Shakoor and Park, 2023). To ease the analogy with Darcy's law,
the present work is restricted to the steady case where the pore scale model is given by
Stokes equations, while the coarse scale model is shown to remain of Stokes type.

In Sec. 2, this pure Stokes approach for modeling flows in porous media is described.
Stokes equations for the coarse scale are shown to depend on some effective properties
of the porous medium that can be computed beforehand in a pre-processing step. Two
methodologies for computing those properties in a Finite Element (FE) implementation
are presented in Sec. 3. It is then demonstrated through comparisons with single-scale
simulations (i.e., simulations where pores are directly discretized at the coarse scale) in
Sec. 4 that the same discretization can be used for both fluid and porous parts, and that
no boundary conditions are necessary at the interface.

## 2 Multiscale model

The coarse scale domain is denoted $\Omega^M \subset \mathbb{R}^d$, with Dirichlet boundary conditions $v_D^M$ on $\Gamma_D^M$ and Neumann boundary conditions $t_N^M$ on $\Gamma_N^M$. In a strong form while neglecting body forces, the equations to satisfy at the coarse scale are

$$ egin{cases} 
\boldsymbol{f}^M(\boldsymbol{x}) - \nabla_{\boldsymbol{x}} \cdot \boldsymbol{\sigma}^{M,\mathrm{dev}}(\boldsymbol{x}) + \nabla_{\boldsymbol{x}} p^M(\boldsymbol{x}) = \mathbf{0}, \boldsymbol{x} \in \Omega^M, \\ 
\nabla_{\boldsymbol{x}} \cdot \boldsymbol{v}^M(\boldsymbol{x}) = 0, \boldsymbol{x} \in \Omega^M, \\ 
\boldsymbol{v}^M(\boldsymbol{x}) = \boldsymbol{v}_D^M(\boldsymbol{x}), \boldsymbol{x} \in \Gamma_D^M, \\ 
(\boldsymbol{\sigma}^{M,\mathrm{dev}}(\boldsymbol{x}) - p^M(\boldsymbol{x})\mathbf{I}) \cdot n(\boldsymbol{x}) = \boldsymbol{t}_N^M(\boldsymbol{x}), \boldsymbol{x} \in \Gamma_N^M, 
\end{cases} 
\quad (1) $$

where $\mathbf{I}$ is the identity matrix and $n$ is the outgoing normal vector to $\partial\Omega^M$. In an equivalent variational form, the coarse scale problem is to find $\boldsymbol{v}^M \in \mathcal{V}^M$, $p^M \in L^2(\Omega^M)$, such that

$$ \int_{\Omega^M} (\boldsymbol{f}^M(\boldsymbol{x}) \cdot \delta\boldsymbol{v}^M(\boldsymbol{x}) + \boldsymbol{\sigma}^{M,\mathrm{dev}}(\boldsymbol{x}) : \nabla_{\boldsymbol{x}}\delta\boldsymbol{v}^M(\boldsymbol{x})) \, d\boldsymbol{x} + \int_{\Omega^M} (-p^M(\boldsymbol{x})\nabla_{\boldsymbol{x}} \cdot \delta\boldsymbol{v}^M(\boldsymbol{x}) - \delta p^M(\boldsymbol{x})\nabla_{\boldsymbol{x}} \cdot \boldsymbol{v}^M(\boldsymbol{x})) \, d\boldsymbol{x} = \int_{\Gamma_N^M} \boldsymbol{t}_N^M(\boldsymbol{x}) \cdot \delta\boldsymbol{v}^M(\boldsymbol{x}) \, d\boldsymbol{x}, \quad (2) $$

$$ \begin{cases} 
\forall \delta\boldsymbol{v}^M \in \mathcal{V}^M, \\ 
\forall \delta p^M \in L^2(\Omega^M), 
\end{cases} 
\quad (3) $$

with functional spaces

$$ 
\begin{align} 
\mathcal{V}^M &= \{ \boldsymbol{w} \in (H^1(\Omega^M))^d, \boldsymbol{w}(\boldsymbol{x}) = \boldsymbol{v}_D^M(\boldsymbol{x}), \forall (\boldsymbol{x}) \in \Gamma_D^M \}, \nonumber \\ 
H^1(\Omega^M) &= \{ w \in L^2(\Omega^M), \nabla_{\boldsymbol{x}} w \in (L^2(\Omega^M))^d \}, \tag{4} \\ 
L^2(\Omega^M) &= \{ w : \Omega^M \mapsto \mathbb{R}, \int_{\Omega^M} w(\boldsymbol{x})^2 \mathrm{d}\boldsymbol{x} < +\infty \}. \nonumber 
\end{align} 
$$

This is standard and easier to understand if incompressible Newtonian flow is modeled directly at the coarse scale. Stokes equations then lead to:

$$ 
\left\{ 
\begin{array}{ll} 
\boldsymbol{f}^M(\boldsymbol{x}) = \mathbf{0}, & \forall \boldsymbol{x} \in \Omega^M, \\ 
\boldsymbol{\sigma}^{M,\mathrm{dev}}(\boldsymbol{x}) = 2\mu^M(\boldsymbol{x})\nabla_{\boldsymbol{x}}^{S,\mathrm{dev}}\boldsymbol{v}^M(\boldsymbol{x}), & 
\end{array} 
\right. 
\quad (5) $$

where the deviatoric part is

$$ \nabla^{S,\mathrm{dev}} \boldsymbol{v}^M = \mathbb{D} : \nabla \boldsymbol{v}^M, \quad (6) $$
the fourth-order tensor $\mathbb{D}$ being defined by

$$ \mathbb{D} = \left( \frac{\delta_{i,k} \delta_{j,l} + \delta_{i,l} \delta_{j,k}}{2} - \frac{\delta_{i,j} \delta_{k,l}}{d} \right)_{ijkl}, \quad \text{with} \quad \delta_{i,j} = \begin{cases} 1, & i = j, \\ 0, & i \neq j, \end{cases} \qquad (7) $$

and for any matrix $\mathbf{C}$, $\mathbb{D} : \mathbf{C} = \left( \sum_{kl} \mathbb{D}_{ijkl} \mathbf{C}_{kl} \right)_{ij}$ and $\mathbf{C} : \mathbb{D} = \left( \sum_{kl} \mathbf{C}_{kl} \mathbb{D}_{klij} \right)_{ij}$. The problem with this direct approach is that the computational cost is going to blow up if there are pores of a small size compared to $\Omega^M$, as the element or cell size should be smaller than the pore size. The alternative proposed herein is a multiscale approach, in which $\mathbf{f}^M(\mathbf{x})$ and $\boldsymbol{\sigma}^{M,dev}(\mathbf{x})$ are unknown functions of $\mathbf{v}^M(\mathbf{x})$ and $\nabla_{\mathbf{x}}\mathbf{v}^M(\mathbf{x})$, while $p^M$ is still a Lagrange multiplier.

Consequently, a fine scale domain $\Omega^m = \Omega^m(\mathbf{x}) \subset \mathbb{R}^d$ is introduced at each point $\mathbf{x}$ of the coarse scale domain. This domain is exclusively composed of closed and/or open pores, where the fluid can flow. The solid rigid part where the fluid cannot flow is not represented. As shown in Fig. 1, $\partial\Omega^m$ includes the boundaries between the fluid and solid parts: $\Gamma_O^m = \Gamma_O^m(\mathbf{x}) \subset \partial\Omega^m(\mathbf{x})$. Dirichlet boundary conditions $\mathbf{v}^m = \mathbf{0}$ are imposed on $\Gamma_O^m$. The implementation of other types of boundary conditions such as the Navier condition may be considered in the future.

The proposed multiscale approach is kinematic, which means averaging constraints are going to be imposed on the velocity field, while the expressions of $\mathbf{f}^M(\mathbf{x})$ and $\boldsymbol{\sigma}^{M,dev}(\mathbf{x})$ are going to be derived from a so-called principle of multiscale virtual power (Blanco et al., 2016; Shakoor and Park, 2023).

For points $\mathbf{x}$ of the coarse scale domain $\Omega^M$ that are inside the fluid part, there is no solid part at the fine scale. For those points which will correspond to integration points of the coarse scale mesh in the numerical implementation, it is possible to use either the single-scale Stokes model from Eq. (5), or a homogeneous fine scale domain where $\Gamma_O^m = \emptyset$ (Shakoor and Park, 2023).

The details on the principle of multiscale virtual power are given in the Appendix, in particular regarding the derivation of the system of variational equations to solve for
<figure><img src="image_2.png" /><figcaption>Figure 1: Example of square fine scale domain $\Omega^m$ in light gray including three holes of different shapes corresponding to solid parts where the fluid cannot flow. The domain boundary $\partial\Omega^m$ in red includes the boundary $\Gamma_O^m$ where Dirichlet boundary conditions $v^m = 0$ are imposed.</figcaption></figure>

each fine scale domain using the FE method. It can be written in a strong form as

$$ 
\left\{
\begin{array}{l}
-\nabla_y \cdot (2\mu^m(\mathbf{y})\nabla_y^{S,dev}\mathbf{v}^m(\mathbf{y})) + \nabla_y p^m(\mathbf{y}) - \boldsymbol{\alpha} = \mathbf{0}, \forall \mathbf{y} \in \Omega^m, \\
\mathbf{v}^m(\mathbf{y}) = \mathbf{0}, \forall \mathbf{y} \in \Gamma_O^m, \\
(2\mu^m(\mathbf{y})\nabla_y^{S,dev}\mathbf{v}^m(\mathbf{y}) - p^m(\mathbf{y})\mathbf{I}) \cdot \mathbf{n}^m(\mathbf{y}) = \boldsymbol{\beta} \cdot \mathbf{n}^m(\mathbf{y}), \forall \mathbf{y} \in \partial\Omega^m \setminus \Gamma_O^m, \\
\nabla_y \cdot \mathbf{v}^m(\mathbf{y}) = 0, \forall \mathbf{y} \in \Omega^m, \\
\frac{1}{|\Omega^m|} \int_{\Omega^m} \mathbf{v}^m(\mathbf{y})d\mathbf{y} = \mathbf{v}^M(\mathbf{x}), \\
\frac{1}{|\Omega^m|} \int_{\Omega^m} \nabla_y \mathbf{v}^m(\mathbf{y})d\mathbf{y} = \nabla_x \mathbf{v}^M(\mathbf{x}),
\end{array}
\right.
\quad (8) $$

where $\mathbf{n}^m(\mathbf{y})$ is the outgoing normal vector at the fine scale domain boundary. It should be noted that this problem is linear as there is no auto-advection term and Newtonian flow is considered.

This problem involves Lagrange multipliers $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ to impose averaging constraints on the fine scale velocity and its gradient. The force per unit volume and deviatoric stress can be computed using $f^M = \boldsymbol{\alpha}$ and $\boldsymbol{\sigma}^{M,dev} = \mathbb{D} : \boldsymbol{\beta}$ as proven in the Appendix.

To summarize:

* The coarse scale problem in Eq. (2) is similar to the well-known Stokes equations
for unsteady incompressible Newtonian flows, except that it is written in terms of an unknown force per unit volume and an unknown deviatoric stress.

* Small pores are not represented directly within the coarse scale. Instead, they are embedded in fine scale domains, which are placed at each point of the coarse scale domain.

* For each fine scale domain, the boundary value problem in Eq. (8) should be solved. The boundary conditions for this problem come from averaging constraints relating the fine scale velocity field and its gradient to their coarse scale counterparts.

* The force per unit volume and deviatoric stress for the coarse scale are computed from the solution of the fine scale problem.

Because the fine scale problem in Eq. (8) is linear, it is possible to compute the effective properties beforehand and use them for running coarse scale simulations. An FExFE (FE²) strategy with two-way coupling between the scales is hence not necessary.

## 3 Numerical method

As in Ref. (Shakoor and Park, 2023), the Taylor-Hood P2/P1 pair is chosen at both scales. Quadratic interpolation is hence used for the velocity, and linear interpolation for the pressure.

Contrary to the unsteady case for which an FE² scheme was proposed in Ref. (Shakoor and Park, 2023), the steady case only requires the computation of the effective properties in a pre-processing step. These properties are the force derivatives $\frac{\partial f^M}{\partial \mathbf{v}^M}$ and $\frac{\partial f^M}{\partial \nabla_x \mathbf{v}^M}$, and the stress derivatives $\frac{\partial \sigma^{M,dev}}{\partial \mathbf{v}^M}$ and $\frac{\partial \sigma^{M,dev}}{\partial \nabla_x \mathbf{v}^M}$. They are independent of $\mathbf{v}^M$ and $\nabla_x \mathbf{v}^M$ due to the linearity of the fine scale problem. They are necessary to compute the force and stress through the following relations:

$$ 
\begin{aligned}
f^M &= \frac{\partial f^M}{\partial \mathbf{v}^M} \cdot \mathbf{v}^M + \frac{\partial f^M}{\partial \nabla_x \mathbf{v}^M} : \nabla_x \mathbf{v}^M, \\
\sigma^{M,dev} &= \frac{\partial \sigma^{M,dev}}{\partial \mathbf{v}^M} \cdot \mathbf{v}^M + \frac{\partial \sigma^{M,dev}}{\partial \nabla_x \mathbf{v}^M} : \nabla_x \mathbf{v}^M.
\end{aligned} 
\quad (9) 
$$
If the fine scale domain is the same for all points of the porous part, for instance if the pore geometry and distribution do not vary, then these properties can be computed once and for all. In the sequel, the FE solution of the fine scale problem is first detailed, and then two approaches based on the FE method are proposed to compute the force and stress derivatives.

## 3.1 Fine scale problem

The fine scale problem in Eq. (8) can be written as

$$ egin{cases} a_{v^m v^m}(\delta \mathbf{v}^m, \mathbf{v}^m) + a_{v^m p^m}(\delta \mathbf{v}^m, p^m) + a_{v^m \alpha}(\delta \mathbf{v}^m, \mathbf{\alpha}) + a_{v^m \beta}(\delta \mathbf{v}^m, \mathbf{\beta}) = \mathbf{0} \\ a_{p^m v^m}(\delta p^m, \mathbf{v}^m) = \mathbf{0} \\ a_{\alpha v^m}(\delta \mathbf{\alpha}, \mathbf{v}^m) = b_{\alpha}(\delta \mathbf{\alpha}) \\ a_{\beta v^m}(\delta \mathbf{\beta}, \mathbf{v}^m) = b_{\beta}(\delta \mathbf{\beta}) \end{cases} \tag{10} $$

with the bilinear forms

$$ egin{align} a_{v^m v^m}(\delta \mathbf{v}^m, \mathbf{v}^m) &= \frac{1}{|\Omega^m|} \int_{\Omega^m} 2\mu^m(\mathbf{y}) \nabla_y^{S,dev} \mathbf{v}^m(\mathbf{y}) : \nabla_y \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y}, \nonumber \\ a_{v^m p^m}(\delta \mathbf{v}^m, p^m) &= -\frac{1}{|\Omega^m|} \int_{\Omega^m} p^m(\mathbf{y}) \nabla_y \cdot \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y}, \nonumber \\ a_{v^m \alpha}(\delta \mathbf{v}^m, \mathbf{\alpha}) &= -\mathbf{\alpha} \cdot \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \right), \nonumber \\ a_{v^m \beta}(\delta \mathbf{v}^m, \mathbf{\beta}) &= -\mathbf{\beta} : \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \nabla_y \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \right), \tag{11} \\ a_{p^m v^m}(\delta p^m, \mathbf{v}^m) &= -\frac{1}{|\Omega^m|} \int_{\Omega^m} \delta p^m(\mathbf{y}) \nabla_y \cdot \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y}, \nonumber \\ a_{\alpha v^m}(\delta \mathbf{\alpha}, \mathbf{v}^m) &= -\mathbf{\alpha} \cdot \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \right), \nonumber \\ a_{\beta v^m}(\delta \mathbf{\beta}, \mathbf{v}^m) &= -\mathbf{\beta} : \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \nabla_y \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \right) \nonumber \end{align} $$

and the linear forms

$$ egin{align} b_{\alpha}(\delta \mathbf{\alpha}) &= -\delta \mathbf{\alpha} \cdot (\mathbf{v}^M(\mathbf{x})), \tag{12} \\ b_{\beta}(\delta \mathbf{\beta}) &= -\delta \mathbf{\beta} : (\nabla_x \mathbf{v}^M(\mathbf{x})). \nonumber \end{align} $$

After FE discretization, the linear system has the following form:

$$ \left( \begin{array}{cccc} A_{v^m v^m} & A_{v^m p^m} & A_{v^m \alpha} & A_{v^m \beta} \\ A_{v^m p^m}^T & 0 & 0 & 0 \\ A_{v^m \alpha}^T & 0 & 0 & 0 \\ A_{v^m \beta}^T & 0 & 0 & 0 \end{array} \right) \left( \begin{array}{c} U_{v^m} \\ U_{p^m} \\ U_{\alpha} \\ U_{\beta} \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \\ B_{\alpha} \\ B_{\beta} \end{array} \right) \tag{13} $$
where the dimension of the vector $U_{v^m}$ is the number of P2 nodes of the FE mesh times $d$, the dimension of the vector $U_{p^m}$ is the number of P1 nodes of the FE mesh, the dimension of the vector $U_\alpha$ is $d$, and the dimension of the vector $U_\beta$ is $d \times d$.

Rows and columns of the linear system in Eq. (13) that are associated to nodes on $\Gamma_O^m$ are modified to set the velocity to zero at those nodes.

## 3.2 Computation of force and stress derivatives

The first approach to compute the derivatives for the coarse scale solve is based on the linear system in Eq. (13). It is, indeed, possible to eliminate the velocity from the first line:

$$ \begin{gathered} A_{v^m v^m} \cdot U_{v^m} + A_{v^m p^m} \cdot U_{p^m} + A_{v^m \alpha} \cdot U_{\alpha} + A_{v^m \beta} \cdot U_{\beta} = 0, \tag{14} \\ \Leftrightarrow U_{v^m} = A_{v^m v^m}^{-1} \cdot (-A_{v^m p^m} \cdot U_{p^m} - A_{v^m \alpha} \cdot U_{\alpha} - A_{v^m \beta} \cdot U_{\beta}), \end{gathered} $$

and then inject it in the third line:

$$ \begin{gathered} A_{v^m \alpha}^T \cdot U_{v^m} = B_\alpha, \\ \Leftrightarrow A_{v^m \alpha}^T \cdot (A_{v^m v^m}^{-1} \cdot (-A_{v^m p^m} \cdot U_{p^m} - A_{v^m \alpha} \cdot U_\alpha - A_{v^m \beta} \cdot U_\beta)) = B_\alpha, \\ \Leftrightarrow \left\{ \begin{array}{l} -A_{v^m \alpha}^T \cdot A_{v^m v^m}^{-1} \cdot A_{v^m \alpha} \cdot U_\alpha = B_\alpha - A_{v^m \alpha}^T \cdot A_{v^m v^m}^{-1} \cdot (-A_{v^m p^m} \cdot U_{p^m} - A_{v^m \beta} \cdot U_\beta), \\ -A_{v^m \alpha}^T \cdot A_{v^m v^m}^{-1} \cdot A_{v^m \beta} \cdot U_\beta = B_\alpha - A_{v^m \alpha}^T \cdot A_{v^m v^m}^{-1} \cdot (-A_{v^m p^m} \cdot U_{p^m} - A_{v^m \alpha} \cdot U_\alpha). \end{array} \right. \tag{15} \end{gathered} $$

Since $f^M = \alpha = U_\alpha$ and $v^M = -B_\alpha$, it can be shown that:

$$ \frac{\partial f^M}{\partial v^M} = \frac{\partial \alpha}{\partial v^M} = -\frac{\partial U_\alpha}{\partial B_\alpha} = (A_{v^m \alpha}^T \cdot A_{v^m v^m}^{-1} \cdot A_{v^m \alpha})^{-1} \tag{16} $$

and since $\sigma^{M,dev} = \mathbb{D} : \beta = \mathbb{D} : (\mathbb{P}_\beta \cdot U_\beta)$, with $\mathbb{P}_\beta = (\delta_{(i-1) \times d+j,k})_{ijk}$, it can be shown that:

$$ \begin{aligned} \frac{\partial \sigma^{M,dev}}{\partial v^M} &= \mathbb{D} : \frac{\partial \beta}{\partial v^M} \\ &= -\mathbb{D} : \left( \mathbb{P}_\beta \cdot \frac{\partial U_\beta}{\partial B_\alpha} \right) \tag{17} \\ &= \mathbb{D} : \left( \mathbb{P}_\beta \cdot (A_{v^m \alpha}^T \cdot A_{v^m v^m}^{-1} \cdot A_{v^m \beta})^+ \right), \end{aligned} $$

where, for any matrix $C$, $C^+$ denotes its Moore-Penrose inverse. Similarly, it can be shown that

$$ \frac{\partial f^M}{\partial \nabla_x v^M} = (A_{v^m \beta}^T \cdot A_{v^m v^m}^{-1} \cdot A_{v^m \alpha})^+, \frac{\partial \sigma^{M,dev}}{\partial \nabla_x v^M} = \mathbb{D} : (A_{v^m \beta}^T \cdot A_{v^m v^m}^{-1} \cdot A_{v^m \beta})^{-1}. \tag{18} $$
This approach might not be suitable in practice, as it involves inverting and multiplying several matrices, which might be prone to both human and numerical errors, especially regarding Moore-Penrose inverses. It has not been investigated in this work, in fact, whether these inverses always exist or not. As an alternative, solving the linear system $d+d \times d$ times could yield all derivatives. For instance, in two dimensions (2D), with

$$ \mathbf{v}^M = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \nabla_x \mathbf{v}^M = \begin{pmatrix} 0 & 0 \\ 0 & 0 \end{pmatrix}, \qquad (19) $$

components $\left(\frac{\partial \boldsymbol{\alpha}_i}{\partial \mathbf{v}_k^M}\right)_{i=1...d,k=1}$ of $\frac{\partial \boldsymbol{\alpha}}{\partial \mathbf{v}^M}$ and components $\left(\frac{\partial \boldsymbol{\beta}_{ij}}{\partial \mathbf{v}_k^M}\right)_{i,j=1...d,k=1}$ of $\frac{\partial \boldsymbol{\beta}}{\partial \mathbf{v}^M}$ would be obtained, and with

$$ \mathbf{v}^M = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \nabla_x \mathbf{v}^M = \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}, \qquad (20) $$

components $\left(\frac{\partial \boldsymbol{\alpha}_i}{\partial \nabla_x \mathbf{v}^M_{kl}}\right)_{i=1...d,k=l=1}$ of $\frac{\partial \boldsymbol{\alpha}}{\partial \nabla_x \mathbf{v}^M}$ and components $\left(\frac{\partial \boldsymbol{\beta}_{ij}}{\partial \nabla_x \mathbf{v}^M_{kl}}\right)_{i,j=1...d,k=l=1}$ of $\frac{\partial \boldsymbol{\beta}}{\partial \nabla_x \mathbf{v}^M}$ would be obtained. Consequently, all components of the force and stress derivatives would be obtained from a total of $d+d \times d$ solves.

## 3.3 Coarse scale problem

The weak form of the coarse scale problem in Eq. (2) can be written as

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega^M} \left( \frac{\partial \mathbf{f}^M}{\partial \mathbf{v}^M}(\mathbf{x}) \cdot \mathbf{v}^M(\mathbf{x}) + \frac{\partial \mathbf{f}^M}{\partial \nabla_x \mathbf{v}^M}(\mathbf{x}) : \nabla_x \mathbf{v}^M(\mathbf{x}) \right) \cdot \delta \mathbf{v}^M(\mathbf{x}) \mathrm{d}\mathbf{x} \\ \displaystyle + \int_{\Omega^M} \left( \frac{\partial \boldsymbol{\sigma}^{M,\mathrm{dev}}}{\partial \mathbf{v}^M}(\mathbf{x}) \cdot \mathbf{v}^M(\mathbf{x}) + \frac{\partial \boldsymbol{\sigma}^{M,\mathrm{dev}}}{\partial \nabla_x \mathbf{v}^M}(\mathbf{x}) : \nabla_x \mathbf{v}^M(\mathbf{x}) \right) : \nabla_x \delta \mathbf{v}^M(\mathbf{x}) \mathrm{d}\mathbf{x} \\ \displaystyle - \int_{\Omega^M} p^M(\mathbf{x}) \nabla_x \cdot \delta \mathbf{v}^M(\mathbf{x}) \mathrm{d}\mathbf{x} = \int_{\Gamma_N^M} \mathbf{t}_N^M(\mathbf{x}) \cdot \delta \mathbf{v}^M(\mathbf{x}) \mathrm{d}\mathbf{x}, \\ \displaystyle - \int_{\Omega^M} \delta p^M(\mathbf{x}) \nabla_x \cdot \mathbf{v}^M(\mathbf{x}) \mathrm{d}\mathbf{x} = 0. \end{array} \right. \qquad (21) $$

As pointed out in Ref. (Shakoor and Park, 2023), the linear system to solve after FE discretization is not different from the standard Taylor-Hood saddle point problem:

$$ \left( \begin{array}{cc} \mathbf{A}_{\mathbf{v}^M} \mathbf{v}^M & \mathbf{A}_{\mathbf{v}^M p^M} \\ \mathbf{A}_{\mathbf{v}^M p^M}^T & \mathbf{0} \end{array} \right) \left( \begin{array}{c} \mathbf{U}_{\mathbf{v}^M} \\ \mathbf{U}_{p^M} \end{array} \right) = \left( \begin{array}{c} \mathbf{B}_{\mathbf{v}^M} \\ \mathbf{0} \end{array} \right). \qquad (22) $$
In fact, provided that force and stress derivatives have been computed in pre-processing, there is no specific difficulty in assembling and solving this linear system as compared to standard Taylor-Hood FEs. This is also true for setting up the Dirichlet boundary conditions. As a consequence, although the computational complexity of the coarse scale problem is slightly higher than one-domain or two-domain methods in the sense that solving Stokes is more expensive than solving Darcy, it is still of the same order. The complexity and the number of fine scale problems to solve in order to compute the force and stress derivatives are also of the same order. Since there are no interface conditions and no additional parameters, the proposed method does not involve any empirical coefficients. Contrary to the one-domain approach, moreover, it does not require the introduction of an inter-region between the fluid and porous parts.

# 4 Simulation cases and results

To keep the computational cost reasonable, only 2D simulations are considered in this section. All simulations are conducted using the *FEMS* software (Shakoor, 2022, 2021). Flow through a coarse scale domain containing both a purely fluid part and a porous part is considered. This domain is the square of area $1 \text{ m}^2$ shown in Fig. 2(a). The fluid part is the light part, while the porous part is the dark one. Boundary conditions force the flow to occur along the horizontal direction, which gives the role of obstacle to the dark part. The fluid dynamic viscosity is fixed to $1 \text{ Pa} \cdot \text{s}$ and gravity is neglected.

In the fluid part, the homogeneous fine scale domain shown in Fig. 2(b) is used, which is equivalent to directly using Stokes equations as per Eq. (5) in this part. With this technique, the size of the fine scale domain has no influence on results (Shakoor and Park, 2023), and it is fixed to $l = 1 \text{ m}$. In the porous part, two heterogeneous fine scale domains may be considered: the one with pore area fraction $\phi = 52\%$ presented in Fig. 2(c), and the one with $\phi = 36\%$ presented in Fig. 2(d). In both cases, the fine scale domain is a square of size $l$ with a circular solid at its center. This solid's interior is not discretized. The radius is $R = 0.39l$ for $\phi = 52\%$ and $R = 0.45l$ for $\phi = 52\%$. Both components of the velocity are set to zero at this solid's boundary. Unless otherwise mentioned, the fine scale domain with $\phi = 52\%$ is chosen and its size is set to $l = 0.1 \text{ m}$.
A different and more realistic setup is proposed in Subsec. 4.5.

## 4.1 Single-scale simulations

Single-scale simulations are first conducted in order to establish reference solutions. Two meshes corresponding to two different mesh sizes are used, as shown in Fig. 3. These meshes rely on a fine discretization close to solids, and a coarser discretization in the fluid part. The average edge length $h$ is used to identify mesh size.

The results of pressure and velocity fields are presented in Fig. 4. On the one hand, the pressure changes mainly in the horizontal direction, with a perturbation due to the porous part of the domain. On the other hand, the horizontal velocity component changes mainly at the boundary between fluid and porous parts, as the flow is significantly slowed down in the porous part.

The evolution of the pressure along a horizontal line at $y = 0.25$ m is shown in Fig. 5(a). The two curves are superimposed, which proves that the so-called coarse mesh is actually quite fine enough for this problem. It should be noted that the pressure is arbitrarily represented as zero inside solids. The velocity curve in Fig. 5(b) confirms this conclusion. A mesh size $h = 4$ µm can hence be used close to solids. For the fine scale domains with solids, this means a mesh size $h = 0.04$ l.

## 4.2 Influence of coarse scale mesh size

The proposed multiscale approach relying on pre-computed force and stress derivatives is now considered. The mesh with only two elements shown in Fig. 6(a) is used for the fine scale domain associated to the fluid part, as it is sufficient to model Stokes flow (Shakoor and Park, 2023). The mesh shown in Fig. 6(b) for the fine scale domain associated to the porous part relies on the same mesh size used close to solids in the single-scale mesh shown in Fig. 3(a).

Regarding the coarse scale mesh, two levels of discretization are considered, a coarse mesh with uniform mesh size $h = 0.04$ m, and a fine mesh with $h = 0.02$ m. These meshes are much coarser and easier to construct than those used in single-scale simulations. They rely on a simultaneous discretization of fluid and porous parts, with an explicit
<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Coarse scale domain and boundary conditions with a dark part and a light part, (b) homogeneous fine scale domain used in the light part, (c) porous fine scale domain with $\phi = 52\%$ used in the dark part, (d) porous fine scale domain with $\phi = 36\%$ used in the dark part. Lengths are in meters and velocities in meters per second.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Single-scale meshes: (a) coarse mesh with $h = 4$ µm close to solids, (b) fine mesh with $h = 2$ µm close to solids.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Results using the coarse single-scale mesh: (a) pressure field (in Pa), (b) velocity component $v_x$ (in m·s⁻¹).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Single-scale results: (a) pressure along a horizontal line at $y = 0.25$ m, (b) velocity component $v_x$ along a vertical line at $x = 0.5$ m.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Fine scale meshes: (a) homogeneous fine scale domain, (b) fine scale domain for $\phi = 52\%$ and with $h = 0.04l$.</figcaption></figure>
discretization of the interface between the two parts. It is reminded, however, that no
condition is imposed at this interface.

Results are shown in Fig. 7 regarding pressure and velocity fields. It is clear that
fine perturbations in Fig. 4 due to each solid in single-scale simulations are not captured
with the proposed multiscale approach.

<figure><img src="image_8.png" /><figcaption>Figure 7: Results using the coarse multiscale mesh: (a) coarse scale pressure field (in Pa), (b) coarse scale velocity component $v_x^M$ (in m·s⁻¹).</figcaption></figure>

A comparison between multiscale simulation results using the two levels of discretiza-
tion at the coarse scale and the reference single-scale solution is presented in Fig. 8. Both
pressure and velocity curves prove that the coarse scale mesh using $h = 0.04$ m is already
fine enough.

The curves for the pressure field in Fig. 8(a) demonstrate that the pressure prediction
from the proposed multiscale approach is accurate in the fluid part. The pressure in the
porous part between solids is also accurate but only in an averaged sense. Regarding
velocity curves in Fig. 8(b), there is a slight difference inside both fluid and porous parts.
As shown in Figs. 8(c,d), this seems to be due to a gap between the boundary of the
porous part in the single-scale model and in the multiscale model. In the single-scale
model, indeed, the porous part actually ends at the boundary of the last solid, while in
the multiscale model it is defined by the geometry in Fig. 2(a). This will be investigated
further in the sequel.

<figure><img src="image_9.png" /><figcaption>Figure 8: Multiscale simulation results compared to the reference single-scale solution: (a) pressure along a horizontal line at $y = 0.25$ m, (b) velocity component $v_x^M$ along a vertical line at $x = 0.5$ m, (c) zoom on (a), (d) zoom on (b). The multiscale pressure curve has been offset to ease the comparison.</figcaption></figure>

## 4.3 Influence of fine scale domain size

Simulations with smaller solids are conducted to investigate the effect of fine scale domain
size. The latter is set to $l = 0.025$ m, which results in solids four times smaller and a
number of solids multiplied by 16. The size of the mesh close to solids, consequently, is divided by four in the single-scale mesh. In fine scale meshes, however, reducing $l$ automatically reduces $h$, with no effect on the total number of elements. The number of elements in the coarse scale mesh also remains constant, as the mesh with $h = 0.04\,\text{m}$ is used.

Results are presented in Fig. 9. They are qualitatively similar to those presented in Figs. 4 and 7. In the single-scale result, however, the fine pressure fluctuations are less visible due to the smaller solid size, and it seems that single-scale and multiscale simulation results are closer.

A quantitative comparison is presented in Fig. 10. It should be noted that these curves show velocity and pressure evolution along lines that do not cross solids, as compared to Fig. 8. As a consequence, the multiscale model appears to be much more accurate inside the porous part. This is also a consequence of reducing $l$.

It is more interesting to look at the gap between the two curves, especially in Fig. 10(d). It is reduced with a smaller fine scale domain size, which leads to a more accurate multiscale prediction also inside the fluid part.

More importantly, the proposed multiscale approach, without any boundary condition between fluid and porous parts, predicts very accurately both pressure and velocity values at the boundary of the porous part. This is particularly remarkable for this simulation where the pore characteristic length is very small compared to the simulation domain size. It is quite advantageous as the computational cost of the proposed approach is actually comparable to the Darcy approach, since it involves the same number of velocity-pressure degrees of freedom. Moreover, it becomes more interesting than a single-scale approach from a computational cost point of view when the pore characteristic length decreases. This was already demonstrated for the unsteady case in Ref. (Shakoor and Park, 2023), but it is more obvious in the present work where force and stress derivatives for coarse scale simulations are pre-computed from a few fine scale simulations. For the results shown in Fig. 9, indeed, the single-scale simulation involves 455974 velocity degrees of freedom while the multiscale simulation involves only 6566 coarse scale velocity degrees of freedom. Force and stress derivatives, moreover, are pre-computed using the fine scale mesh shown in Fig. 6(b) through $2+2 \times 2 = 6$ solutions involving 1008 fine scale velocity
<figure><img src="image_10.png" /><figcaption>Figure 9: Results for a fine scale domain size $l = 0.025$ m: (a) pressure field (in Pa) from the single-scale simulation, (b) velocity component $v_x$ (in m·s⁻¹) from the single-scale simulation, (c) coarse scale pressure field (in Pa) from the multiscale simulation, (d) coarse scale velocity component $v_x^M$ (in m·s⁻¹) from the multiscale simulation.</figcaption></figure>
degrees of freedom.

<figure><img src="image_11.png" /><figcaption>Figure 10: Multiscale and single-scale results for a fine scale domain size $l = 0.025$ m: (a) pressure along a horizontal line at $y = 0.25$ m, (b) velocity component $v_x^M$ along a vertical line at $x = 0.5$ m, (c) zoom on (a), (d) zoom on (b). The multiscale pressure curve has been offset to ease the comparison.</figcaption></figure>

## 4.4 Influence of porosity

Investigating the influence of porosity is also interesting as decreasing it leads to larger solids and hence decreases the gap between the boundary of the porous part in single-scale and multiscale models. The porous fine scale domain shown in Fig. 2(d) is used to compute force and stress derivatives. Mesh size is kept at $h = 0.04l$, with $l = 0.1$ m. The
single-scale mesh is similar to that shown in Fig 3(a), except that it has larger solids.

Pressure and velocity curves are shown in Fig. 11. Reducing the porosity improves the performance of the multiscale model. This is clearly demonstrated by comparing Figs. 11(b) and 8(b). More precisely, the gap between the two curves is reduced in Figs. 11(c,d) as compared to Figs. 8(c,d). The proposed multiscale approach is therefore more accurate with a smaller porosity. This could be expected as the gap between the boundary of the porous part in single-scale and multiscale models is decreased with a smaller porosity.

Another possible explanation could be a capability of the multiscale model to accurately approximate a completely solid part with no porosity.

## 4.5 Application to a woven textile composite

Finally, simulations are conducted for a more realistic case inspired from the woven textile composite material characterized in Ref. (Syerko et al., 2023). As shown in Fig. 12(a), the rectangular domain of area $1 \times 0.1\text{m}^2$ contains a porous part which is an elliptical yarn of major radius 0.25 m and minor radius 0.025 m. This leads to an aspect ratio of 10, which is similar to what is observed in the experimental images (Syerko et al., 2023). Boundary conditions and fluid properties are identical to previous simulations. The fine scale domain is the one with pore area fraction $\phi = 52\%$ presented in Fig. 2(c), with three sizes $l = 10 \text{ mm}, 5 \text{ mm}, 2.5 \text{ mm}$. In this application, the holes in Fig. 2(c) actually represent fibers, and a yarn typically contains thousands of fibers whose diameter depends on the material (carbon, glass, etc.). Composite manufacturing processes involve liquid resin flow around yarns as well as within yarns around the fibers.

Although it has been refined close to the porous part boundary to improve the accuracy, the coarse scale mesh shown in Fig. 12(a) contains only 14190 elements. The domains used for single-scale simulations are presented in Fig. 12(b-d). They are discretized with meshes that contain up to 4118396 elements for $l = 2.5 \text{ mm}$.

Results are presented in Fig. 13, where the velocity component $v_x^M$ is shown along a diagonal line. For this example with more realistic porous medium boundary topology, the accuracy of multiscale predictions is again increased when the pore characteristic length
<figure><img src="image_12.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d)</figcaption></figure>

Figure 11: Multiscale and single-scale results for a porosity $\phi = 36\%$: (a) pressure along a horizontal line at $y = 0.25$ m, (b) velocity component $v_x^M$ along a vertical line at $x = 0.5$ m, (c) zoom on (a), (d) zoom on (b). The multiscale pressure curve has been offset to ease the comparison.
<figure><img src="image_16.png" /><figcaption>Figure 12: (a) Coarse scale mesh with an elliptical porous part for the composite application, (b) single-scale domain for $l = 10$ mm, (c) single-scale domain for $l = 5$ mm, (d) single-scale domain for $l = 2.5$ mm, (e) zoom on (d).</figcaption></figure>
becomes very small compared to the simulation domain size. This is particulary relevant
for this application where the most representative setup is the one with $l = 2.5$ mm and
it features the greatest number of fibers.

<figure><img src="image_17.png" /><figcaption>Figure 13: Multiscale and single-scale results for the composite application: (a) velocity component $v_x^M$ along the diagonal line joining the lower left corner to the upper right corner of the coarse scale domain, (b) zoom on (a). The vertical line on each graph marks the porous medium boundary.</figcaption></figure>

# 5 Conclusions

An alternative approach relying on equations of Stokes type instead of Darcy's law has been proposed in this paper to model flow in porous media. These equations depend on force and stress derivatives tensors that can be computed from pore scale simulations. The latter are based on a fictitious fine scale domain where pores are discretized. The equations to solve for the fine scale domain in order to compute the force and stress derivatives tensors have been presented, as well as two Finite Element (FE) procedures for the numerical implementation.

Similarly to the permeability tensor in Darcy's law, the derivatives tensors can be used for simulations at the coarse scale with discretizations that are no longer required to deal with very small pores. The main advantage of the proposed approach is that it eases greatly the task of modeling flows in domains containing both fluid and porous parts.
While with Darcy's law, specific boundary conditions including the well-known Beavers-Joseph-Saffman condition should be imposed at the boundary between fluid and porous parts, the proposed approach circumvents this difficulty. This advantage both in terms of modeling and numerical implementation is obtained thanks to the use of equations of Stokes type instead of Darcy's law, which means the equations used in both fluid and porous parts are of the same type in the proposed approach.

Numerical simulations considering a domain including a fluid and a porous part have been conducted. Particular attention has been given to pressure and velocity evolution across the boundary between the two parts. It has been shown by comparisons with single-scale simulations that the proposed approach is accurate, and that its accuracy is improved with a lower pore characteristic length or volume fraction.

The proposed pure Stokes approach for coupling fluid flow with porous media flow opens quite a range of interesting possibilities. Among them, the implementation of numerical methods relying on an implicit representation of the interface between fluid and porous parts (e.g., the level-set method) should be facilitated.

# 6 Appendix

This appendix details the derivation of the proposed multiscale model and in particular the governing equations for the fine scale problem.

## 6.1 Kinematic averaging

The average of the fine scale velocity field is constrained to be equal to the coarse scale velocity at the corresponding point $x$ of the coarse scale domain:

$$v^M(x) = \frac{1}{|\Omega^m|} \int_{\Omega^m} v^m(y) dy. \qquad (A.1)$$

Similarly, the average of the fine scale velocity gradient field is constrained to be equal to the coarse scale velocity gradient at the corresponding point $x$ of the coarse scale domain:

$$\nabla_x v^M(x) = \frac{1}{|\Omega^m|} \int_{\Omega^m} \nabla_y v^m(y) dy. \qquad (A.2)$$

These constraints are additional to the boundary condition $v^m = 0$ at $\Gamma_O^m$.
## 6.2 Kinematic admissibility

Kinematic admissibility requires the averaging constraints for the velocity field to be satisfied. These constraints will all be added using Lagrange multipliers except for the Dirichlet boundary conditions on $\Gamma_O^m$, thus the functional spaces for the fine scale unknowns are

$$ v^m \in \mathcal{V}^m, \qquad (A.3) \\ \mathcal{V}^m = \{ w \in (H^1(\Omega^m))^d, w(y) = 0, \forall y \in \Gamma_O^m \}. $$

Given that additional Lagrange multipliers are set, these functional spaces can also be used for virtual variations.

## 6.3 Virtual powers

The coarse scale virtual power $P^M$ is based on the coarse scale problem stated in Eq. (2):

$$ P^M(\delta V^M, \delta G^M) = f^M \cdot \delta V^M + \sigma^{M,dev} : \delta G^M - p^M \mathrm{tr}(\delta G^M), \qquad (A.4) \\ \forall (\delta V^M, \delta G^M) \in \mathbb{R}^d \times \mathbb{R}^{d \times d}. $$

The averaged fine scale virtual power $P^m$ is based on Stokes equations for incompressible Newtonian flow written at the fine scale:

$$ P^m(\delta V^M, \delta G^M, \delta v^m, \delta p^m, \delta \alpha, \delta \beta) = \left\{ \begin{array}{l} \displaystyle \frac{1}{|\Omega^m|} \int_{\Omega^m} 2\mu^m(y) \nabla_y^{S,dev} v^m(y) : \nabla_y \delta v^m(y) \mathrm{d}y \\ \displaystyle -\frac{1}{|\Omega^m|} \int_{\Omega^m} p^m(y) \nabla_y \cdot \delta v^m(y) \mathrm{d}y \\ \displaystyle -\frac{1}{|\Omega^m|} \int_{\Omega^m} \delta p^m(y) \nabla_y \cdot v^m(y) \mathrm{d}y \\ \displaystyle -\delta \alpha \cdot \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} v^m(y) \mathrm{d}y - v^M(x) \right) \\ \displaystyle -\delta \beta : \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \nabla_y v^m(y) \mathrm{d}y - \nabla_x v^M(x) \right) \\ \displaystyle -\alpha \cdot \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \delta v^m(y) \mathrm{d}y - \delta V^M \right) \\ \displaystyle -\beta : \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \nabla_y \delta v^m(y) \mathrm{d}y - \delta G^M \right) \\ \\ \forall (\delta V^M, \delta G^M) \in \mathbb{R}^d \times \mathbb{R}^{d \times d}, \\ \\ \forall (\delta v^m, \delta p^m) \in \mathcal{V}^m \times L^2(\Omega^m), \\ \\ \forall (\delta \alpha, \delta \beta) \in \mathbb{R}^d \times \mathbb{R}^{d \times d}. \end{array} \right. \qquad (A.5) $$
The divergence constraint on the fine scale velocity gradient has been added through La-
grange multiplier $p^m \in L^2(\Omega^m)$ and the kinematic averaging constraints on the fine scale
velocity field and its gradient have been added through Lagrange multipliers $(\alpha, \beta) \in$
$\mathbb{R}^d \times \mathbb{R}^{d \times d}$. Similarly, Lagrange multipliers have been added for the variations of the fine
scale unknowns.

6.4 Principle of multiscale virtual power

The principle of multiscale virtual power establishes the balance between the two scales.
Therefore, the problem is to find $(\mathbf{v}^m, p^m, \boldsymbol{\alpha}, \boldsymbol{\beta}) \in \mathcal{V}^m \times L^2(\Omega^m) \times \mathbb{R}^d \times \mathbb{R}^{d \times d}$, with
$(\mathbf{v}^M, p^M) \in \mathcal{V}^M \times L^2(\Omega^M)$ such that

$$
P^M(\delta \mathbf{V}^M, \delta \mathbf{G}^M) = P^m(\delta \mathbf{V}^M, \delta \mathbf{G}^M, \delta \mathbf{v}^m, \delta p^m, \delta \boldsymbol{\alpha}, \delta \boldsymbol{\beta}), \quad (A.6)
$$

which is equivalent to

$$
\left.
\begin{array}{l}
\displaystyle \mathbf{f}^M \cdot \delta \mathbf{V}^M \\
\displaystyle + \boldsymbol{\sigma}^{M, \mathrm{dev}} : \delta \mathbf{G}^M \\
\displaystyle - p^M \mathrm{tr}(\delta \mathbf{G}^M)
\end{array}
\right\}
=
\left\{
\begin{array}{l}
\displaystyle \frac{1}{|\Omega^m|} \int_{\Omega^m} 2\mu^m(\mathbf{y}) \nabla_{\mathbf{y}}^{S, \mathrm{dev}} \mathbf{v}^m(\mathbf{y}) : \nabla_{\mathbf{y}} \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \\
\displaystyle - \frac{1}{|\Omega^m|} \int_{\Omega^m} p^m(\mathbf{y}) \nabla_{\mathbf{y}} \cdot \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \\
\displaystyle - \frac{1}{|\Omega^m|} \int_{\Omega^m} \delta p^m(\mathbf{y}) \nabla_{\mathbf{y}} \cdot \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \\
\displaystyle - \delta \boldsymbol{\alpha} \cdot \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} - \mathbf{v}^M(\mathbf{x}) \right) \\
\displaystyle - \delta \boldsymbol{\beta} : \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \nabla_{\mathbf{y}} \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} - \nabla_{\mathbf{x}} \mathbf{v}^M(\mathbf{x}) \right) \\
\displaystyle - \boldsymbol{\alpha} \cdot \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} - \delta \mathbf{V}^M \right) \\
\displaystyle - \boldsymbol{\beta} : \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \nabla_{\mathbf{y}} \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} - \delta \mathbf{G}^M \right)
\end{array}
\right\} \quad (A.7)
\\[1.5em]
\forall (\delta \mathbf{V}^M, \delta \mathbf{G}^M) \in \mathbb{R}^d \times \mathbb{R}^{d \times d}, \\
\forall (\delta \mathbf{v}^m, \delta p^m) \in \mathcal{V}^m \times L^2(\Omega^m), \\
\forall (\delta \boldsymbol{\alpha}, \delta \boldsymbol{\beta}) \in \mathbb{R}^d \times \mathbb{R}^{d \times d}.
$$

By taking arbitrary variations from Eq. (A.7), it is possible to derive the problem to
solve at the fine scale with all necessary boundary conditions, and also to derive the
expressions of the missing terms for the coarse scale problem in Eq. (2).
## 6.5 Fine scale problem

As the spaces for the variations are very large, it is possible to zero out $\delta V^M$, $\delta G^M$. This leads to the problem of finding $(v^m, p^m, \alpha, \beta) \in V^m \times L^2(\Omega^m) \times \mathbb{R}^d \times \mathbb{R}^{d \times d}$, with $(v^M, p^M) \in V^M \times L^2(\Omega^M)$ such that

$$ 
\left. 
\begin{aligned}
& \frac{1}{|\Omega^m|} \int_{\Omega^m} 2\mu^m(\mathbf{y}) \nabla_{\mathbf{y}}^{S,dev} \mathbf{v}^m(\mathbf{y}) : \nabla_{\mathbf{y}} \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \\
& - \frac{1}{|\Omega^m|} \int_{\Omega^m} p^m(\mathbf{y}) \nabla_{\mathbf{y}} \cdot \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \\
& - \frac{1}{|\Omega^m|} \int_{\Omega^m} \delta p^m(\mathbf{y}) \nabla_{\mathbf{y}} \cdot \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \\
& - \delta \alpha \cdot \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} - \mathbf{v}^M(\mathbf{x}) \right) \\
& - \delta \beta \cdot \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \nabla_{\mathbf{y}} \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} - \nabla_{\mathbf{x}} \mathbf{v}^M(\mathbf{x}) \right) \\
& - \alpha \cdot \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \right) \\
& - \beta \cdot \left( \frac{1}{|\Omega^m|} \int_{\Omega^m} \nabla_{\mathbf{y}} \delta \mathbf{v}^m(\mathbf{y}) \mathrm{d}\mathbf{y} \right) \\
& \quad \forall (\delta \mathbf{v}^m, \delta p^m) \in V^m \times L^2(\Omega^m), \\
& \quad \forall (\delta \alpha, \delta \beta) \in \mathbb{R}^d \times \mathbb{R}^{d \times d}.
\end{aligned}
\right\} = 0 \qquad (A.8)
$$

## 6.6 Coarse scale force per unit volume and stress

By zeroing out all variations except $\delta V^M$ in Eq. (A.7), one obtains

$$ f^M \cdot \delta V^M = -\alpha \cdot (-\delta V^M) = \alpha \cdot \delta V^M, \forall \delta V^M \in \mathbb{R}^d \Leftrightarrow f^M = \alpha. \qquad (A.9) $$

Similarly, it is possible to obtain

$$ 
\begin{aligned}
\sigma^{M,dev} : \delta G^M - p^M \mathrm{tr}(\delta G^M) &= -\beta : (-\delta G^M) = \beta : \delta G^M, \forall \delta G^M \in \mathbb{R}^{d \times d}, \\
\Leftrightarrow \sigma^{M,dev} &= \beta + p^M \mathbf{I} = \mathbb{D} : \beta.
\end{aligned}
\qquad (A.10)
$$

## References

Auriault, J. L. (2009). On the domain of validity of Brinkman's equation. Transport in Porous Media, 79(2):215-223.

Auriault, J.-L. (2010). About the Beavers and Joseph Boundary Condition. Transport in Porous Media, 83(2):257-266.
Blanco, P. J., Clausse, A., and Feijóo, R. A. (2017). Homogenization of the Navier-Stokes equations by means of the Multi-scale Virtual Power Principle. *Computer Methods in Applied Mechanics and Engineering*, 315:760–779.

Blanco, P. J., Sánchez, P. J., de Souza Neto, E. A., and Feijóo, R. A. (2016). Variational Foundations and Generalized Unified Theory of RVE-Based Multiscale Models. *Archives of Computational Methods in Engineering*, 23(2):191–253.

Chidyagwai, P. and Rivière, B. (2010). Numerical modelling of coupled surface and subsurface flow systems. *Advances in Water Resources*, 33(1):92–105.

Darcy, H. (1856). *Les fontaines publiques de la ville de Dijon*. Dalmont, Paris.

Eggenweiler, E. and Rybak, I. (2020). Unsuitability of the Beavers–Joseph interface condition for filtration problems. *Journal of Fluid Mechanics*, 892:A10.

Jäger, W., Mikelic, A., and Neuss, N. (2001). Asymptotic Analysis of the Laminar Viscous Flow Over a Porous Bed. *SIAM Journal on Scientific Computing*, 22(6):2006–2028.

Layton, W. J., Schieweck, F., and Yotov, I. (2002). Coupling fluid flow with porous media flow. *SIAM Journal on Numerical Analysis*, 40(6):2195–2218.

Marciniak-Czochra, A. and Mikelić, A. (2012). Effective Pressure Interface Law for Transport Phenomena between an Unconfined Fluid and a Porous Medium Using Homogenization. *Multiscale Modeling & Simulation*, 10(2):285–305.

Marco Discacciati and Alfio Quarteroni (2009). Navier-Stokes/Darcy Coupling: Modeling, Analysis, and Numerical Approximation. *Revista Matematica Complutense*, 22(2):315–426.

Marušić-Paloka, E. and Pažanin, I. (2022). The effective boundary condition on a porous wall. *International Journal of Engineering Science*, 173:103638.

Münzenmaier, S. and Starke, G. (2011). First-Order System Least Squares for Coupled Stokes–Darcy Flow. *SIAM Journal on Numerical Analysis*, 49(1):387–404.
Neuman, S. P. (1977). Theoretical derivation of Darcy's law. *Acta Mechanica*, 25(3-4):153-170.

Pacquaut, G., Bruchon, J., Moulin, N., and Drapier, S. (2012). Combining a level-set method and a mixed stabilized P1/P1 formulation for coupling Stokes-Darcy flow. *International Journal for Numerical Methods in Fluids*, 69:459-480.

Payne, L. E. and Straughan, B. (1998). Analysis of the boundary condition at the interface between a viscous fluid and a porous medium and related modelling questions. *Journal des Mathematiques Pures et Appliquees*, 77(4):317-354.

Shakoor, M. (2021). FEMS – A Mechanics-oriented Finite Element Modeling Software. *Computer Physics Communications*, 260:107729.

Shakoor, M. (2022). FEMS — Finite Element Modeling Software. https://hal.science/hal-03781711.

Shakoor, M. and Park, C. H. (2023). Computational homogenization of unsteady flows with obstacles. *International Journal for Numerical Methods in Fluids*, 95(4):499-527.

Sudhakar, Y., Lācis, U., Pasche, S., and Bagheri, S. (2021). Higher-Order Homogenized Boundary Conditions for Flows Over Rough and Porous Surfaces. *Transport in Porous Media*, 136:1-42.

Syerko, E., Schmidt, T., May, D., Binetruy, C., Advani, S., Lomov, S., Silva, L., Abaimov, S., Aissa, N., Akhatov, I., Ali, M., Asiaban, N., Broggi, G., Bruchon, J., Caglar, B., Digonnet, H., Dittmann, J., Drapier, S., Endruweit, A., Guilloux, A., Kandinskii, R., Leygue, A., Mahato, B., Martínez-Lera, P., Matveev, M., Michaud, V., Middendorf, P., Moulin, N., Orgéas, L., Park, C., Rief, S., Rouhi, M., Sergeichev, I., Shakoor, M., Shishkina, O., Swolfs, Y., Tahani, M., Umer, R., Vanclooster, K., and Vorobeyev, R. (2023). Benchmark exercise on image-based permeability determination of engineering textiles: Microscale predictions. *Composites Part A: Applied Science and Manufacturing*, 167:107397.

Urquiza, J., N'Dri, D., Garon, A., and Delfour, M. (2008). Coupling Stokes and Darcy equations. *Applied Numerical Mathematics*, 58(5):525-538.
Valdés-Parada, F. J. and Lasseux, D. (2021a). A novel one-domain approach for modeling flow in a fluid-porous system including inertia and slip effects. *Physics of Fluids*, 33(2).

Valdés-Parada, F. J. and Lasseux, D. (2021b). Flow near porous media boundaries including inertia and slip: A one-domain approach. *Physics of Fluids*, 33(7):1–15.

Whitaker, S. (1966). The equations of motion in porous media. *Chemical Engineering Science*, 21(3):291–300.

Whitaker, S. (1986). Flow in porous media I: A theoretical derivation of Darcy’s law. *Transport in Porous Media*, 1(1):3–25.

## Funding

This work was supported by the French Agence Nationale de la Recherche (ANR) through the MISSA projet (ANR-22-CE46-0002).

## Competing Interests

The authors have no relevant financial or non-financial interests to disclose.

## Author Contributions

**Modesar Shakoor and Chung Hae Park:** Conceptualization, Validation, Formal analysis, Resources, Writing - Review & Editing and Supervision; **Modesar Shakoor:** Methodology, Software, Investigation, Data Curation, Writing - Original Draft, Visualization, Project administration and Funding acquisition.