A computational continuum model of poroelastic beds

TL;DR

This paper introduces a validated continuum model for poroelastic beds that accurately predicts fluid-structure interactions at the macroscale, validated against microscopic simulations, and applicable across various engineering fields.

Contribution

The paper develops and validates a non-empirical continuum model for poroelastic beds, capturing microstructural effects and interface behaviors accurately.

Findings

Effective model captures microstructure anisotropy effects

Model predictions agree with microscopic simulations

Applicable to diverse fields like bio-engineering and geophysics

Abstract

Despite the ubiquity of fluid flows interacting with porous and elastic materials, we lack a validated non-empirical macroscale method for characterizing the flow over and through a poroelastic medium. We propose a computational tool to describe such configurations by deriving and validating a continuum model for the poroelastic bed and its interface with the above free fluid. We show that, using stress continuity condition and slip velocity condition at the interface, the effective model captures the effects of small changes in the microstructure anisotropy correctly and predicts the overall behaviour in a physically consistent and controllable manner. Moreover, we show that the performance of the effective model is accurate by validating with fully microscopic resolved simulations. The proposed computational tool can be used in investigations in a wide range of fields, including mechanical engineering, bio-engineering and geophysics.

Equations

$${\rho }_{\mathrm{f}}\left[{\partial }_t\mathrm{u}+\left(\mathrm{u}\cdot \nabla \right)\mathrm{u}\right]$$

$$\nabla \cdot \mathrm{u}$$

$${\rho }_{\mathrm{s}}{\partial }_t^2\mathrm{v}$$

$$\mathrm{u}={\partial }_t\mathrm{v}\text{mbox}and\left\{-p\delta +\mu \left[\nabla \mathrm{u}+{\left(\nabla \mathrm{u}\right)}^T\right]\right\}\cdot \hat{n}=\left\{{\mathrm{C}}^{\mathrm{sk}}:\frac{1}{2}\left[\nabla \mathrm{v}+{\left(\nabla \mathrm{v}\right)}^T\right]\right\}\cdot \hat{n}.$$

$$\left(1-\theta \right){\rho }_{\mathrm{s}}{\partial }_t^2\mathrm{v}$$

$$\mathrm{u}-\theta {\partial }_t\mathrm{v}=-\frac{\mathrm{K}}{\mu }\cdot \nabla p.$$

$$\nabla \cdot \mathrm{u}$$

$$\mathcal{E}{\partial }_{t}p-\nabla \cdot \left(\frac{\mathrm{K}}{\mu }\cdot \nabla p\right)$$

$${\mathrm{v}|}_{\mathrm{rw}}=0\text{mbox}and{\frac{\partial p}{\partial \hat{n}}|}_{\mathrm{rw}}=0.$$

$$p^{-}=p.$$

$$\left[\mathrm{C}:\frac{1}{2}\left(\nabla \mathrm{v}+{\left(\nabla \mathrm{v}\right)}^T\right)-\alpha p\right]\cdot \hat{n}=\left[-p\delta +\mu \left(\nabla \mathrm{u}+{\left(\nabla \mathrm{u}\right)}^T\right)\right]\cdot \hat{n}.$$

$$\mathrm{u}\cdot \hat{\tau }={\partial }_t\mathrm{v}\cdot \hat{\tau }+\left(-\frac{{\mathrm{K}}^{\mathrm{if}}}{\mu }\cdot \nabla p^{-}+\mathrm{L}:\left[\nabla \mathrm{u}+{\left(\nabla \mathrm{u}\right)}^T\right]\right)\cdot \hat{\tau },$$

$$\mathrm{u}\cdot \hat{n}={\partial }_t\mathrm{v}\cdot \hat{n}-\left(\frac{\mathrm{K}}{\mu }\cdot \nabla p^{-}\right)\cdot \hat{n},$$

$$C_{ijkl}^{\text{sk}}=E\left\{\frac{\nu }{\left(1+\nu \right)\left(1-2\nu \right)}{\delta }_{ij}{\delta }_{kl}+\frac{1}{2\left(1+\nu \right)}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)\right\}.$$

$$\mathrm{C}=\left(1-\theta \right){\mathrm{C}}^{\mathrm{sk}}+{\mathrm{C}}^{\mathrm{sk}}:⟨\frac{1}{2}\left[\nabla \chi +{\left(\nabla \chi \right)}^T\right]⟩,$$

$$⟨f⟩=\frac{1}{l^3}{\int }_{V_{\sigma }}fdV,$$

$$\nabla \cdot \left[{\mathrm{C}}^{\mathrm{sk}}:\frac{1}{2}\left\{\nabla \chi +{\left(\nabla \chi \right)}^T\right\}\right]$$

$$\left[{\mathrm{C}}^{\mathrm{sk}}:\frac{1}{2}\left\{\nabla \chi +{\left(\nabla \chi \right)}^T\right\}\right]$$

$$\boldsymbol{\mathrm{{C}}}=\left(\begin{array}[]{cccccc}4.792&0.239&0.239&0&0&0\\ 0.239&4.794&0.239&0&0&0\\ 0.239&0.239&4.792&0&0&0\\ 0&0&0&0.246&0&0\\ 0&0&0&0&0.246&0\\ 0&0&0&0&0&0.246\\ \end{array}\right)\cdot 10^{-2}E.$$

$$\nabla \cdot \left[\frac{{\mathrm{C}}^{\mathrm{sk}}}{E}:\frac{1}{2}\left\{\nabla \eta +{\left(\nabla \eta \right)}^T\right\}\right]$$

$$\left[\frac{{\mathrm{C}}^{\mathrm{sk}}}{E}:\frac{1}{2}\left\{\nabla \eta +{\left(\nabla \eta \right)}^T\right\}\right]$$

$$\alpha =\theta \delta +⟨\frac{{\mathrm{C}}^{\mathrm{sk}}}{E}:\frac{1}{2}\left[\nabla \eta +{\left(\nabla \eta \right)}^T\right]⟩,$$

$$\mathcal{E}=⟨\frac{\nabla \cdot \eta }{E}⟩.$$

$$\alpha =0.9789\delta \text{and}\mathcal{E}=0.1571E^{-1}.$$

$$-\nabla \mathcal{A}+{\nabla }^2\mathcal{K}$$

$$\nabla \cdot \mathcal{K}$$

$$\mathrm{K}=⟨\mathcal{K}⟩=2.32\cdot 10^{-2}{l}^2\delta .$$

$$\boldsymbol{\mathrm{{C}}}=\left(\begin{array}[]{cccccc}4.668&0.283&0.233&0&0.139&0\\ 0.283&4.638&0.209&0&-0.002&0\\ 0.233&0.209&4.031&0&0.084&0\\ 0&0&0&0.168&0&-0.010\\ 0.139&-0.002&0.084&0&0.205&0\\ 0&0&0&-0.010&0&0.258\\ \end{array}\right)\cdot 10^{-2}\,E.$$

$$\boldsymbol{\mathrm{{\alpha}}}=\left(\begin{array}[]{ccc}0.9793&0&-0.0009\\ 0&0.9795&0\\ -0.0009&0&0.9821\\ \end{array}\right),$$

$$\mathcal{E}=0.1465E^{-1}.$$

$$\boldsymbol{\mathrm{{K}}}=\left(\begin{array}[]{ccc}2.51&0&-0.10\\ 0&2.25&0\\ -0.10&0&2.22\end{array}\right)\cdot 10^{-2}\,l^{2}.$$

$$\mathrm{u}=\mathrm{U}+{\mathrm{u}}^{+}.$$

Full text

{fmtext}

1 Introduction

Recent advances in surface micro- and nano-fabrication techniques[1, 2, 3] are providing new technological oppor-tunities of enormous potential. However, we lack high-fidelity models capturing how the underlying small-scale physicochemical processes interact with the large-scale flow and heat- and mass-transport phenomena. The reason is the vast range of scales in both time and space that need be resolved in order to capture the full physical picture, which renders full-scale numerical investigations extremely costly[4]. The development of multi-scale models with reduced complexity is a necessary enabler in an increasing number of applications where the control of events at the small length scales determines the properties of the flow over much larger space and time scales. Examples include, but are not limited to,

the design of novel surfaces to control heat transfer, skin friction or pressure drag, acoustic noise and fluid mixing.

There is thus a clear need to develop techniques that combine computational fluid mechanics with mesoscopic or microscopic models of materials. In this work, we take a step in this direction by modelling the interaction of mesoscopic surface textures with macroscopic flows. By “mesoscopic” we mean structures that can be described by continuum methods, but are still significantly smaller than the large-scale flow phenomena. When these constructions, which can take the form of high-aspect ratio structures (such as fibres) or granular-type structures (such as particles), are saturated with a fluid and can absorb the fluid stress and possibly deform, we have a poroelastic material.

The description of poroelastic medium as a continuum has an extensive and broad history, beginning with the empirical models of porous media by Darcy[5] and followed by the formulation of total stress tensor of deformable porous media by Biot[6, 7]. Today, there exists a significant amount of work[8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29] on poroelastic media adopting top-down approaches based empirical macroscale stress/deformation tests, bottom-up approaches using homogenization and volume averaging techniques, physical approaches using analytical and mechanical models, or other methods.

In this work, we employ the method of homogenization via multi-scale expansion to model flows through and over poroelastic surfaces. Within the general multi-scale analysis platform (MAP) classification proposed by Scheibe et al.[30], the current method falls into the category of methods for which the different scales can be completely decoupled from each other (e.g. “motif B” in [30]). In other words, our work is “formal upscaling”, which means that after the macroscale governing equations have been derived, the parameters of the model are governed by fully decoupled microscale problems (from now on, we use the term “microscale” instead of “mesoscale” to conform with the terminology within two-scale expansion approach). This method limits our investigations to set-ups for which the underlying microscale closure problems are linear. For example, the fluid phase in the poroelastic material must be slow enough such that inertial effects can be neglected or modelled through some kind of linearisation. Essentially all previous upscaling works on deformable porous media has been in the same motif B category. Examples include the method of volume averaging by Whitaker[31] and the method of homogenization, as employed, for example, by Mei & Vernescu[12].

There exists now a number of studies [32, and references within] on development and applications of the motif B multi-scale methods for poroelastic media. It is, however, the authors’ opinion that we still lack a computational framework that lays out the sequential steps needed to be taken in order to obtain the physical parameters describing the poroelastic medium not only in the interior of the material, but also at the interface with freely moving fluid. Ideally, such a framework, when provided a particular microscopic structure in terms of its geometry and properties (skeleton elasticity, connectivity, etc.), will provide the anisotropic macroscopic material properties (permeability, elasticity, etc.) of the effective continuum fields without any fitting parameters obtained from experiments.

To the best of the authors’ knowledge, the theoretical frameworks of Whitaker[16] and Mei & Vernescu[12] have not been validated from microscopic point of view for the flow over a surface, which is coated with porous and elastic material. Validations of these methods most often consider only macroscopic measures. More specifically, while microscale problems in unit cells needed in the upscaling procedures have been presented and solved previously[16, 33, 12], there is no comparison between global macroscale simulations – using properties obtained from those microscale solutions – and corresponding fully resolved simulations. In this work, we do not only compute effective material tensors by numerically solving microscale problems, but also compare the obtained effective continuum description with fully resolved direct numerical simulation (DNS) of the fluid flow inside and above the poroelastic surface. In this way, we can assess the accuracy of the effective model quantitatively both in the interior and near the interface of the medium with a free-flowing fluid.

A particularly important contribution is the insight it provides on how the velocity and the stress interface conditions between the poroelastic medium and the free fluid region model the microscale effects in averaged manner. Gopinath & Mahadevan[29] as well as Minale[34] point out that the effective interface condition for a poroelastic region may in fact be more straightforward to match, compared to the rigid porous material, since there is a natural way to balance the fluid stress from the free fluid with the solid stress of the surface material. However, there exists no validation – where the interface boundary conditions for displacement, fluid velocity and pressure are coupled to the Stokes equations above – for a non-trivial flow, where there is transport across the interface. Most works treating the boundary conditions are empirical[35, 36, 37, 38]. Those contributions which have recently treated the interface problems from first-principles have focused on rigid porous media and one-dimensional problems, such as the laminar channel flow [39, 40] only (for which there is no transfer of mass or momentum between the material and the free fluid), or infiltration flow only[41].

In summary, the objectives of the current work are to (i) present a framework, derived using multi-scale expansion, suitable to model flow through and over poroelastic materials, (ii) validate the framework with respect to the fully resolved direct numerical simulations, and (iii) gain insight on stress transfer near the interface and evaluate the accuracy of interface boundary conditions. Our work is partially the numerical counterpart of the analytical and asymptotic study by Gopinath & Mahadevan[29]; we aim at computing (instead of physically modelling) physical material properties of anisotropic poroelastic materials as well as the interface with the free flow. While Gopinath & Mahadevan[29] consider particularly biologically relevant microstructures of ordered or disordered filaments, we consider connected materials consisting of linked spheres and ellipsoids.

This paper is organized as follows. In section 2, we introduce both the microscale and the macroscale/effective equations governing our problem, which consists of a poroelastic material at the interface with a moving fluid. A method to compute the effective properties of poroelastic materials is presented in section 3. In the same section, we provide effective tensor results for cubic-symmetric and monoclinic-symmetric poroelastic materials and analyse them. A lid-driven cavity problem to investigate poroelastic material response to a steady two-dimensional flow vortex is proposed and solved using the homogenized equations in section 4. In the same section, we report results obtained from resolved direct numerical simulations and explain the shear stress transfer between free fluid and poroelastic medium. In section 5, we discuss the limitations of the presented theory. Finally, in section 6 we conclude this work and outline future directions.

2 Micro- and macroscale equations describing a poroelastic bed

For a dense poroelastic medium exposed to a free flow (Fig. 1), one may define at least two length scales; a microscopic (pore) scale $l$, which characterizes the size of voids in the material and a macroscopic (global) scale $H$, which characterizes the size of the large-scale processes in the medium or nearby. The large scale flow is characterized also by the created pressure difference $\Delta P$. In this section, we present the microscale equations that resolve every scale in the full physical domain, and a set of macroscale equations that model the effective average behavior of the poroelastic bed under a given set of assumptions. Although the effective field equations for describing a fluid-saturated poroelastic material are known [6, 29], we re-derive them from first principles in the supplementary appendices. The main reason is to uncover the detailed microscale problems in unit cells that are required for determining the physical coefficients appearing in the macroscale systems.

2.1 Microscale governing equations

The microscale equations resolving the fluid-structure physics at scale $\sim l$ are the conservation of mass and momentum. For a Newtonian fluid with density $\rho_{\mathrm{f}}$ and viscosity $\mu$, the momentum in the free fluid region and within the porous skeleton is governed by the incompressible Navier-Stokes equations

[TABLE]

where $\boldsymbol{\mathrm{u}}$ and $p$ are velocity and pressure fields, respectively.

The material is defined by the solid skeleton density $\rho_{\mathrm{s}}$ and linear elasticity tensor $\boldsymbol{\mathrm{{C}}}^{\mathrm{sk}}$. Assuming an isotropic material the elasticity tensor is defined by Young’s modulus $E$ and Poisson’s ratio $\nu$. The solid skeleton momentum is governed by a balance between solid inertia and stress, obtained using the linear stress-strain relationship

[TABLE]

where $\boldsymbol{\mathrm{v}}$ is the displacement field of the solid skeleton. In order to couple the fluid and structure problems, the no-slip condition and continuity of stress is prescribed at the boundary between solid skeleton and surrounding fluid, i.e.,

[TABLE]

Here, $\boldsymbol{\mathrm{{\delta}}}$ is the second-rank identity tensor and $\hat{n}$ is unit-normal vector at the boundary. Solving the governing equations everywhere at the pore scale is computationally very expensive due to a globally large domain and due to requirement of fine resolution near the pores. This motivates the development of an alternative continuum description, where the pore fluid and solid are considered as one composite, which given appropriate equations and constitutive relations presents the average behavior of the actual poroelastic bed.

2.2 Effective field equations

We divide the physical domain into two parts; one containing only the free fluid, and the other containing the fluid and solid skeleton (Fig. 2). The free-fluid region is governed by the Navier-Stokes equations (1–2).

The continuum description of the composite (solid and fluid) is based on a separation between the pore-scale and the system scale. Mathematically this can be formulated by a scale separation parameter $\epsilon=l/H\ll 1$.

The effective displacement field $\boldsymbol{\mathrm{v}}$ of the homogenized poroelastic material (see Fig. 2) is governed by,

[TABLE]

Here, $\theta=V_{\mathrm{f}}/V$ is the porosity (which in general is a function of space, but in this work is a constant), where $V_{\mathrm{f}}$ is the fluid phase volume and $V$ is the total volume of the composite in one unit cell (defined later). Moreover, $\boldsymbol{\mathrm{{C}}}$ is the fourth-rank effective elasticity tensor of the material and $\boldsymbol{\mathrm{{\alpha}}}$ is the coefficient for the contribution of the pore pressure ($p$) in the total stress. Whereas the microscale elasticity $\boldsymbol{\mathrm{{C}}}^{\textrm{sk}}$ only depends on $E$ and $\nu$, the effective tensor $\boldsymbol{\mathrm{{C}}}$ in general depends on the porosity, pore geometry (as modelled in the current work) and on the particular type of boundary condition imposed near the interface with the free fluid, or the impermeable walls[29]. The tensor $\boldsymbol{\mathrm{{\alpha}}}$ does not have a correspondence in the microscale and is an effect of solid skeleton deformation due to seepage flow through the pores. We characterize in detail both $\boldsymbol{\mathrm{{C}}}$ and $\boldsymbol{\mathrm{{\alpha}}}$ for two different poroelastic media in section 3.

Moreover, for sufficiently dense poroelastic material the fluid flow between the pores is slow, such that inertial effects are negligible. Therefore the pore pressure is the dominant contribution in the fluid flow and the leading-order equation is the (relative) Darcy’s law,

[TABLE]

This expression relates the gradient of the pore pressure, the solid velocity $\partial_{t}\boldsymbol{\mathrm{v}}$ and the flow field $\boldsymbol{\mathrm{u}}$ in the poroelastic medium to each other. Here $\boldsymbol{\mathrm{{K}}}$ is the interior permeability tensor. The term $\theta\partial_{t}\boldsymbol{\mathrm{v}}$ arises due to viscous friction between the solid structure and the pore fluid; if there is a motion of the solid skeleton, the surrounding pore fluid is dragged along through boundary condition (4).

Finally, the conservation of mass requires that

[TABLE]

where $\boldsymbol{\mathrm{{D}}}=\theta\boldsymbol{\mathrm{{\delta}}}-\boldsymbol{\mathrm{{\alpha}}}$ is a dimensionless second-rank tensor. It determines how the strain of the displacement (caused either by the flow through the pores or by a boundary condition) modifies the solid structure volume within one pore, thus squeezing the pore fluid in or out of the pore. The scalar $\mathcal{E}$ characterizes the change of solid structure volume within one pore with respect to time varying pressure, which, similarly to solid strain, can cause change in pore fluid content and consequently introduce apparent compressibility of the flow field.

The system of equations (5-7) determines the seven unknowns ($\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{v}}$ and $p$), but by combining (6) with (7), the system can be reduced to four unknowns ($\boldsymbol{\mathrm{v}}$ and $p$). The fluid velocity $\boldsymbol{\mathrm{u}}$ can then be computed as a postprocessing step once pressure and displacement fields are known. Following this approach, the equation for the pore pressure is

[TABLE]

Note that the effective system above for describing a poroelastic system can be derived using multi-scale expansion (see section 3 and supplementary appendices), motivated using mixture theory, or physically modelled. The major challenge is to use this developed effective system with appropriate boundary conditions at interfaces with solid walls, free fluids or other structures to describe problems arising from various applications. The aim of this work is to provide a framework from which one can form a fully closed effective system to describe the response and interaction of a poroelastic material with a surrounding free fluid. In the next section, we therefore provide the needed boundary conditions between a poroelastic material and a rigid wall, and also a free fluid.

2.3 Effective interface conditions

2.3.1 Conditions for the poroelastic bed

To solve the governing equations for the poroelastic material, one needs to impose boundary conditions for both the pore-pressure (equation 8) and the displacement field (equation 5). On rigid walls, we impose (similar to [29]) zero displacement and zero transpiration (normal fluid velocity), which through the relative Darcy’s equation (6) leads to

[TABLE]

Here, “$\mathrm{rw}$” means “rigid wall” and $\hat{n}$ is the unit-normal vector at the wall. Physically, the no-slip condition should be satisfied at the wall, but this is not compatible with the leading-order presentation of a poroelastic media based on Darcy’s law. The Darcy’s law only describes the direct proportionality between the pore-pressure gradient and the velocity, and does not include any macroscopic diffusion effects.

At the artificial interface with the free fluid, shown in Fig. 2, we impose a pressure continuity condition

[TABLE]

Here, the pore pressure below the interface is denoted by $p^{-}$ and pressure of the free fluid above the interface with $p$. Note that the choice of condition for pressure depends on the assumptions made about the flow as well as the material geometry. For example, Lācis & Bagheri[42] have shown that, if the interface velocity caused by the shear stress is of the same order as the velocity induced by the pore pressure gradient, the pressure continuity is the leading order boundary condition for any pore geometry. On the other hand, Mikelić et al.[43, 44] have shown that, if velocity contribution from the shear stress at the interface is one order higher than contribution from the pressure gradient, there is a pressure jump for anisotropic pore geometry. However, if the pore geometry is isotropic/cubic-symmetric, the pressure continuity is still applicable. Finally, for the material displacement, we impose continuity of stresses at the interface

[TABLE]

We thus assume, similarly to [29], that the total stress of the free fluid is transferred to the total effective stress of the interior poroelastic medium. In general, the effective elasticity of the composite near the interface could however be different from its value in the interior (it is argued in general to depend on boudary conditions[29]) and, to arrive to a more accurate boundary condition, one could construct an interface-cell with an elasticity problem, similarly as done for the velocity boundary condition[42]. One objective of this paper is to understand how the fluid shear stress is transferred to the solid stress (first across the interface then inside the bed) and if the interface correction is necessary; in section 4, we will show that approximation of the interface effective stress with the interior parameters is able to capture the transfer of stress reasonably well.

2.3.2 Conditions for the free fluid

To solve for the Navier-Stokes equations in the free fluid domain, we need to impose velocity – in principle, one could also have stress condition, but it is already used for the displacements of the poroelastic material – boundary conditions at the interface. In this work, we extend the velocity boundary condition derived by Lācis & Bagheri[42] for a rigid porous bed to include poroelasticity111The discussion can be found in supplementary Appendix B, after equation (B 77).. The condition for the tangential interface velocity is

[TABLE]

where the unit vector $\hat{\tau}$ denotes both tangential directions of the surface. Note that the pressure gradient is the pore-pressure gradient from poroelastic material side of the interface (hence the minus superscript), whereas the flow velocity field $\boldsymbol{\mathrm{u}}$ is on the free-fluid side. The interface velocity has two distinct contributions; i) the no-slip contribution governed by the movement of the solid structure; ii) the slip contribution, which is caused by the porosity of the solid structure and depends both on pore pressure gradient and free-fluid shear. The slip contribution is characterized by the second-rank interface permeability tensor $\boldsymbol{\mathrm{{K}}}^{\mathrm{if}}$ and the third-rank slip length tensor $\boldsymbol{\mathrm{{L}}}$. For a dense material, the first slip term scales as $l^{2}$ and is significantly smaller than the second term, which scales as $l$. The motion of the poroelastic material or the no-slip contribution depends not only on the pore-length $l$, but also on the elasticity of the material and the flow regime. The condition for the normal “penetration” interface velocity component is set by mass conservation, i.e.

[TABLE]

where $\boldsymbol{\mathrm{{K}}}$ is the interior permeability tensor defined in equation (6) of the porous medium and $\nabla p^{-}$ denotes the pressure gradient when approaching the interface from the bed. Similarly as for the tangential component, the velocity has two parts – the no-slip part and the “slip” part. The slip in the normal direction is essentially the fluid mass transport in and out of the poroelastic material, which has to be equal to the relative velocity from the interior. This same condition also arises from the interface cell[42], because the tensor $\boldsymbol{\mathrm{{L}}}$ components corresponding to penetration interface velocity are zero.

In the next section, we will compute the physical parameters which characterize the interior poroelastic medium ($\boldsymbol{\mathrm{{C}}}$, $\boldsymbol{\mathrm{{\alpha}}}$, $\boldsymbol{\mathrm{{K}}}$ and $\mathcal{E}$) as well as the parameters characterizing the interface with the free fluid ($\boldsymbol{\mathrm{{K}}}^{\mathrm{if}},\boldsymbol{\mathrm{{L}}}$) for two specific pore-scale geometries. In section 4, we will investigate in more details interface conditions at boundary between poroelastic material and free fluid by comparing the effective model against the microscopically-resolved system using the lid-driven cavity problem.

3 Properties of poroelastic material exposed to free fluid

We now turn to determining the effective properties of a poroelastic material exposed to the free fluid, by introducing and solving a set of particular microscale problems in an interior unit-cell and in an interface-cell. We focus on a periodic solid structure that can be obtained by duplicating a single pore structure in all directions. The sketch in Fig. 3$a$ shows how the material and its interface with free-fluid is divided into cubic interior cells and elongated rectangular interface-cells. The effective parameters are computed by solving two elasticity problems (to obtain $\boldsymbol{\mathrm{{C}}}$, $\boldsymbol{\mathrm{{\alpha}}}$ and $\mathcal{E}$) and one fluid problem (to obtain $\boldsymbol{\mathrm{{K}}}$) in the interior cell and two fluid problems in the interface-cell (to obtain $\boldsymbol{\mathrm{{K}}}^{\mathrm{if}},\boldsymbol{\mathrm{{L}}}$). These microscale problems are first illustrated for the interior domain using a weakly and a strongly anisotropic microstructure, before we move on to the microscale problems at the interface with the free fluid.

3.1 Interior of the poroelastic material

3.1.1 Microscale problems for elasticity parameters

The effective poroelastic bed obtained from the unit-cell approach can not have fully isotropic elastic properties due to the boundaries between cubic unit cells, as shown in Fig. 3a. The resulting effective elasticity $\boldsymbol{\mathrm{{C}}}$ will at most exhibit a cubic symmetry; a symmetry, which is defined by planes – parallel to the unit-cell sides and diagonal across the unit-cell in all directions – going through the center of the cube.

In order to preserve cubic symmetry, the unit-cell can be filled with any structure that itself is cubic symmetric, such as Wigner-Seitz grains [33, 12] or cube with spherical holes [45]. In this work, we use a sphere at the center of unit-cell with radius $R$, which is connected to neighbouring cells via circular cylinders of radius $r$, as shown in Fig. 3b. The solid skeleton structure is built from isotropic elastic material with elasticity tensor

[TABLE]

The effective elasticity tensor can be determined (see Appendix B for derivation) using a third-rank displacement test tensor $\boldsymbol{\mathrm{{\chi}}}$ as

[TABLE]

where we define transpose of the fourth-rank tensor $\boldsymbol{\nabla}\boldsymbol{\mathrm{{\chi}}}$ acting on the first two indices $\left(\boldsymbol{\nabla}\boldsymbol{\mathrm{{\chi}}}\right)^{T}_{ijkl}=\left(\boldsymbol{\nabla}\boldsymbol{\mathrm{{\chi}}}\right)_{jikl}$. The expression above provides a direct linear link between the micro- and macroscale elasticity tensors through the known porosity $\theta$ and the pore-scale geometry captured by $\boldsymbol{\mathrm{{\chi}}}$ that will be computed by solving equations (16–17). Here, it is assumed that the skeleton elasticity tensor $\boldsymbol{\mathrm{{C^{\mathrm{sk}}}}}$ is constant in space. The brackets denote the volume average over the interior unit-cell volume, which is filled either by solid or by fluid

[TABLE]

where $V_{\sigma}$ is the volume of either the solid or fluid phase. In general, the effective properties should be re-evaluated as the porosity changes. However, we assume small displacements of the effective system so that the porosity is roughly constant and therefore it is enough to compute the averages only once.

The microscale displacement tensor $\boldsymbol{\mathrm{{\chi}}}$ is defined in the solid phase (having volume $V_{s}$) of the unit-cell and is the solution of the following problem,

[TABLE]

Here, the fourth-rank identity tensor appearing in the boundary condition is defined as $\left(\boldsymbol{\mathrm{{\delta^{(4)}}}}\right)_{ijkl}=\delta_{ik}\delta_{jl}$. This equation corresponds to a standard steady linear-elasticity problem generalised to third-rank tensor subjected to different combinations of surface loading (applied on the interface between the solid skeleton and the pore fluid) in order to characterise the response of the structure to all possible surface loading scenarios. To complete the formulation, periodic boundary conditions are applied to solid surfaces, which are in contact with the boundary of the unit-cell. In order to render the solution unique, we impose constraints on average values of displacement[12] using penalty terms in weak formulation[46].

The unit-cell domain has been meshed222We use mesh spacing $\Delta s_{1}=0.025\,l$ at the solid skeleton and $\Delta s_{2}=0.10\,l$ at the boundaries of the unit-cell. We have also carried out the simulation on coarser meshes ($\Delta s_{1}=0.05\,l$ and $\Delta s_{2}=0.20\,l$; $\Delta s_{1}=0.10\,l$ and $\Delta s_{2}=0.40\,l$) and, by observing the convergence of results, we set entries containing only numerical noise to zero. using GMSH software[47] and equations (16 – 17) has been solved using FreeFEM++[48]. The sphere radius is $R=0.3\,l$ and the cylinder radius is $r=0.1\,l$, which results in porosity $\theta=0.85$. Assuming the Poisson’s ratio to be $\nu=0.3$ (due to the linearity of the problem, the solution is valid for any Young’s modulus value $E$), the effective elasticity tensor for the cubic-symmetric poroelastic medium in Voigt notation[49, p. 136] is

[TABLE]

Note that the entries of the effective elasticity tensor are significantly smaller (up to 300 times) compared to these of the skeleton elasticity. This is because most of the solid mass has been removed, leaving only $15\%$ solid volume fraction if compared to a completely filled case, hence the resulting material is much softer. Another observation is that the shear coefficients (diagonal elements in $3\times 3$ bottom right matrix block) are now much smaller compared to pressure-wave coefficients[50, p. 22] (diagonal elements in $3\times 3$ top left matrix block). This effect can be attributed to connecting net of cylinders between spheres; the cylinders are much easier to bend, compared to bulk material, i.e. there is no continuous support from the sides, as there would be in the continuous material case. The effective tensor can be characterized by three independent parameters333The difference of $2\cdot 10^{-5}$ between the first and second second diagonal element is attributed to numerical error that could be caused by, for example, imperfections in the generated mesh. and has the form of a material with cubic symmetry [49, p. 99].

The other two elasticity-related effective parameters $\boldsymbol{\mathrm{{\alpha}}}$ and $\mathcal{E}$ can be computed through expressions444For explicit relationships obtained using multi-scale expansion, see supplementary appendix C. involving the inverse of the fourth-rank elasticity tensor effective elasticity $\boldsymbol{\mathrm{{C}}}$ [29, eqs. 2.4,2.7,2.8]. Alternatively, one can solve for an additional test displacement field $\boldsymbol{\mathrm{\eta}}$ governed by the linear elasticity problem

[TABLE]

Note that this problem is driven by different surface forcing compared to the $\boldsymbol{\mathrm{{\chi}}}$ problem, but is otherwise solved under the same conditions. Given $\boldsymbol{\mathrm{\eta}}$, the pore-pressure contribution tensor is computed from,

[TABLE]

and the coefficient describing the poroelastic material response to the time-variation of the pore pressure is computed from

[TABLE]

By solving (19–20) for the particular microstructure with cubic symmetry shown in Fig. 3$b$, we obtain

[TABLE]

The tensor $\boldsymbol{\mathrm{{\alpha}}}$ is not unity, therefore one can conclude that the solid phase of the poroelastic material is not incompressible. This conclusion agrees with the used elasticity parameter (we use $\nu=0.3$, while incompressible solids have $\nu=0.5$). However, it is non-trivial to use tensor $\boldsymbol{\mathrm{{\alpha}}}$ as a measure of compressibility, because it is also a function of porosity $\theta$.

3.1.2 Microscale problem for permeability tensor

Turning attention to the fluid flow, the poroelastic material is characterized by a permeability [5], which in our framework is obtained by solving the following Stokes problem in the interior unit-cell (Fig. 3b),

[TABLE]

Here, $\boldsymbol{\mathrm{{\mathcal{K}}}}$ is a second-rank tensor field; $\mathcal{K}_{ij}$ is the $i$th velocity component of $j$th vectorial Stokes problem associated with the pressure field $\mathcal{A}_{j}$. The three vectorial Stokes problems characterize how the pore flow responds to volume forcing in one spatial direction at a time. This equation system is complemented with no-slip boundary condition at the interface with the solid skeleton, and periodic boundary conditions at the sides of the unit-cell. For uniqueness, we require that the average pressure field is zero. The effective permeability tensor for the cubic symmetric material is then obtained by averaging the field $\boldsymbol{\mathrm{{\mathcal{K}}}}$ over the fluid volume in the unit-cell as

[TABLE]

We observe that the permeability tensor is characterized by a single constant value, which is characteristic for isotropic flow in the poroelastic medium. Hence, the cubic-symmetry of the poroelastic material is not visible in permeability tensor.

3.1.3 Effects of significant anisotropy at the pore-scale

With the aim of providing a framework that can cope with any periodic microstructure, we compute the effective parameters for a strongly anisotropic microstructure, namely, a tilted ellipsoid at the center of the unit-cell, as shown in Fig. 4. Note that the cubic symmetry in the poroelastic material is broken by this structure, since it has less symmetry planes compared to the cubic unit cell.

The ellipsoid semi-axes have lengths $a_{1}$, $a_{2}$ and $a_{3}$ in $x$, $y$ and $z$ directions, respectively. The ellipsoid can be rotated around the $y$ axis by an angle $\phi$ as shown in Fig. 4b. In order to assess how anisotropy at the pore-scale modifies the effective parameters at the macroscale, we choose parameters that break as many symmetries as possible; $a_{1}=0.4\,l$, $a_{2}=0.3\,l$ and $a_{3}=0.2\,l$ and $\phi=30^{\circ}$. This results in a porosity $\theta=0.86$ and structure with only one symmetry plane.

Using the expression (14), we evaluate the effective elasticity tensor for the titled ellipse structure by solving the microscale displacement problem for $\boldsymbol{\mathrm{{\chi}}}$ (16– 17),

[TABLE]

The form of $\boldsymbol{\mathrm{{C}}}$ correponds to a monoclinic material symmetry [49, p. 96], and nearly all the elements of the tensor differ from each other, as expected due to the anisotropic pore-structure. This symmetry property is invariant with respect to transformation of coordinate axes; by using a different coordinate system, one would only be able to reposition the zero entries in the elasticity tensor to different rows or columns. Moreover, the magnitude of the elements are of the same order as for the cubic-symmetric geometry, since the two cases have nearly the same porosity.

By solving for the displacement field $\boldsymbol{\mathrm{\eta}}$ problem (19–20), we obtain via expression (21),

[TABLE]

where we observe that the tensor cannot be characterized by a single constant as for the cubic-symmetric material, but by four distinct coefficients. The tensor $\boldsymbol{\mathrm{{\alpha}}}$ can be understood as a measure of volume force in the poroelastic medium caused by pore-pressure gradient via seepage flow. The diagonal terms are different because of the different lengths of the ellipsoid semi-axes $a_{1}$, $a_{2}$ and $a_{3}$. Due to this difference, the area of the solid structure (see Fig. 4$c$), exposed to flow in $x$, $y$ and $z$ directions, is different. Consequently, it results in a different drag force for the same pore-pressure gradient in different directions. The off-diagonal term is, on the other hand, caused by the tilt of the ellipsoid, which gives raise to non-zero pressure force, when projected in $x$ direction, due to flow in $z$ direction and vice-versa. The elastic response to time variation of pressure is

[TABLE]

The slightly smaller value for the ellipsoid (compared to the sphere) can be attributed to the shape of the pore-structure at the center of the unit-cell. The ellipsoid is thinner in one direction compared to sphere, which leads to larger strain from the same displacement values and consequently a body which is less compressible.

By solving the fluid problem (22–23) at the pore-scale for the monoclinic symmetric geometry we obtain effective permeability tensor

[TABLE]

The second-rank permeability tensor has a similar features as the pore-pressure contribution tensor $\boldsymbol{\mathrm{{\alpha}}}$. It has four independent parameters, which are again set by four geometric parameters – semi-axes lengths $a_{1}$, $a_{2}$ and $a_{3}$ as well as turn angle $\phi$. The off-diagonal term shows that a tilted ellipse generates flow in the $x$ direction if exposed to pressure a gradient along $z$, and vice-versa. This effect is the same as for a tilted plate exposed to an incoming parallel free stream. Due to the tilt, a tangential flow with respect to the incoming free stream appears.

3.2 Microscale problems for interface conditions

In order to finalize the homogenized model, one needs effective tensors for the velocity boundary conditions of the free fluid in contact with the poroelastic material. We can determine the velocity boundary conditions (the interface permeability $\boldsymbol{\mathrm{{K}}}^{\mathrm{if}}$ and the Navier-slip tensor $\boldsymbol{\mathrm{{L}}}$) by solving a set of microscale problems in an interface-cell (Fig. 3$a$) as shown for rigid porous media by Lācis & Bagheri[42].

To derive the microscale problems in the interface unit cell, we decompose the flow above the interface into a fast flow $\boldsymbol{\mathrm{U}}$ and a perturbation $\boldsymbol{\mathrm{u}}^{+}$ flow

[TABLE]

Below the interface there is only slow flow $\boldsymbol{\mathrm{u}}^{-}$. The flow perturbation velocity is the cause of the slip velocity at the interface with porous or poroelastic material[42]. Note that the global pressure difference $\Delta P$ is driving the fast flow $\boldsymbol{\mathrm{U}}$ above the interface as well as the slow flow $\boldsymbol{\mathrm{u}}^{-}$ below the interface, whereas $\boldsymbol{\mathrm{u}}^{+}$ is driven by the processes in the poroelastic material. In order to determine the perturbation velocity, one can arrive to a solution in the interface-cell (Fig. 3$a$) using an assumption of linear superposition

[TABLE]

in which the unknown tensorial fields $\boldsymbol{\mathrm{{\mathcal{K}}}}^{\pm}$ and $\boldsymbol{\mathrm{{\mathcal{L}}}}^{\pm}$ are velocity response fields to test forcing, which corresponds to pressure gradient (volume forcing) and surface stress (surface forcing). Fields $\boldsymbol{\mathrm{{\mathcal{K}}}}^{-}$ and $\boldsymbol{\mathrm{{\mathcal{L}}}}^{-}$ are correspondingly defined below the interface (illustrated using the dash red line), and fields $\boldsymbol{\mathrm{{\mathcal{K}}}}^{+}$ and $\boldsymbol{\mathrm{{\mathcal{L}}}}^{+}$ are defined above the interface, as sketched in Fig. 5. In the next sections, we introduce problems for all of the introduced test fields and solve them for both microscale geometries.

3.2.1 Interface permeability tensor $\boldsymbol{\mathrm{{K}}}^{\mathrm{if}}$

The interface permeability second-rank tensor field below the interface $\boldsymbol{\mathrm{{\mathcal{K}}}}^{-}$ is governed by

[TABLE]

The boundary conditions at the sides of this domain are periodic (same as for interior cell), however, at the bottom of the domain one has to use the interior solution provided by the interior problem (22– 23), see Fig. 5$a$. At the interface, the continuity of permeability fields,

[TABLE]

is employed. Since the permeability field in the upper part of the interface-cell is also unknown, it can be found using an unforced Stokes momentum equation:

[TABLE]

The volume force in Stokes equations above the interface does not exist, because the global driving pressure above the interface acts on the fast flow, not on the perturbation (27). To complete the two-domain formulation of the interface-cell, incompressibility constraints are added, no-slip condition is imposed at the surface of solid skeleton, periodic boundary conditions are applied at the sides, zero stress condition at the top of the cell, and stress continuity at the interface,

[TABLE]

as shown in Fig. 5$a$. One can observe from equations (28–31) that the outlined problem is a combination of Stokes problems exposed to unit volume-forcing in all possible directions below the interface. For example, the first column of the tensor $\boldsymbol{\mathrm{{\mathcal{K}^{\pm}}}}$ field ($\mathcal{K}^{\pm}_{i1}$) corresponds to flow response to unit volume forcing in the $x$ direction, as shown in Fig. 5$a$. Thus, this test problem is used to characterize the pore-pressure gradient (which is a volume force) contribution to velocity near the interface. Finally, after solving the coupled two-domain Stokes problem, one can obtain the effective interface permeability, by employing volume average $\boldsymbol{\mathrm{{K}}}^{\mathrm{if}}=\langle\boldsymbol{\mathrm{{\mathcal{K}}}}^{+}\rangle$ in $l^{3}$ cube above the interface.

Evaluating555Similarly as in interior simulations, we carry out investigations using three different resolutions, $\Delta s_{1}=0.05\,l$ and $\Delta s_{2}=0.20\,l$, $\Delta s_{1}=0.0625\,l$ and $\Delta s_{2}=0.2500\,l$ and $\Delta s_{1}=0.10\,l$ and $\Delta s_{2}=0.40\,l$. The interface-cell is chosen to encompass three structures below the interface and is extending for the same distance in the free fluid. the interface permeability tensor gives

[TABLE]

for cubic-symmetric and monoclinic-symmetric pore-scale geometries, respectively. One can observe that the interface-normal permeability components (for velocity component $u_{z}$, last row of the matrix) are the same as interior ones (set by conservation of mass[42]), but tangential components are different. The interface permeability for the tangential components is larger than the interior one, which can be explained by the fact that there is no friction from the upper part of the interface-cell, while the interior cell would be exposed to friction from neighbouring cells. In addition, the interface permeability matrix for monoclinic symmetric pore geometry is no more symmetric, that is, the effect of anisotropy has vanished for $u_{x}$ velocity component at the interface. The reason for this coefficient vanishing could be the fact that the anisotropy is too weak to overcome dissipation at the cylinder on the top of the last micro-structure and does not contribute for the interface velocity.

3.2.2 Interface Navier-slip tensor $\boldsymbol{\mathrm{{L}}}$

When there is a flow over a porous or poroelastic material, there is a slip velocity induced proportional to free fluid shear stress, as theoretically derived by Mikelić & Jäger[39] for a one-dimensional channel flow and by Lācis & Bagheri[42] for a general three-dimensional set-up.

The slip length problem, unlike the permeability problem, is specific to the interface, i.e. there exists no analogous problem for the interior domain. The slip third-rank tensor fields are also governed by Stokes equations

[TABLE]

below and above the interface, respectively, see Fig. 5$b$. As before, the fields in two domains are connected through continuity condition $\boldsymbol{\mathrm{{\mathcal{L}}}}^{-}=\boldsymbol{\mathrm{{\mathcal{L}}}}^{+}$. Momentum equations in this test problem are not exposed to volume test forcing. The only non-triviality in this coupled two-domain problem is jump in stresses across the interface

[TABLE]

where the third-rank interface stress jump tensor is defined as $\left(\boldsymbol{\mathrm{{\mathcal{J}}}}\right)_{ikl}=\delta_{ik}n_{l}$. At the sides of the interface cell we employ periodic boundary conditions, and at the top of the interface cell the zero-stress boundary condition is used, and at the solid structure we have no-slip condition. At the bottom of the interface-cell, however, the no-slip condition is used, because there is no $\boldsymbol{\mathrm{{\mathcal{L}}}}$ field in the interior. By looking at the governing equations of the test problem (33–35) one can conclude that this test problem is exposed to test forcing at the interface. For example, the problem with second and third tensor indices being $j=1$ and $k=3$ ($\mathcal{L}_{i13}$) corresponds to unit forcing in the $x$ direction at the interface, as shown in Fig. 5$b$. Thus this test problem is used to characterize the free-fluid surface stress contribution to the interface velocity. Finally, after solving the coupled two-domain Stokes problem, one can obtain the effective interface slip length, by employing the same volume average $\boldsymbol{\mathrm{{L}}}=\langle\boldsymbol{\mathrm{{\mathcal{L}}}}^{+}\rangle$ above the interface.

Evaluating the slip length tensor at the interface gives

[TABLE]

for cubic-symmetric and monoclinic-symmetric pore-scale geometries, respectively. In the planar interface case, the only meaningful interface problems for the Navier-slip length are those which relate slip length with velocity shear in interface-normal direction[42], therefore we have presented the slip length tensor dot product with unit normal vector $\hat{z}$ normal to the interface with poroelastic material (other entries in this third-rank tensor are zero). The non-zero entries in both matrices corresponds to influence on velocity components $u_{x}$ and $u_{y}$ from velocity strains $\left(\partial_{z}u_{x}+\partial_{x}u_{z}\right)$ and $\left(\partial_{z}u_{y}+\partial_{y}u_{z}\right)$, respectively. For the cubic-symmetric case, both coefficients are similar due to cubic symmetry of the pore-scale geometry. That is, the structure (see Fig. 3$b$) is the same in the $x$ and the $y$ directions. The monoclinic-symmetric structure, on the other hand, is different in the $x$ and the $y$ directions (see Fig. 4$c$). However, the slip-length coefficients are still similar, because (i) the top cylinder acts exactly the same way in both $x$ and $y$ directions and (ii) the tilted ellipsoid slows down the velocity in both $x$ and $y$ directions similarly in integral sense. For $u_{z}$ component, there is no contribution in both cases, because the interface-normal penetration velocity is governed by mass conservation alone.

4 Poroelastic material response to free fluid vortex above it

The purpose of this section is to exemplify the proposed numerical framework by illustrating the response of the cubic-symmetric and monoclinic-symmetric poroelastic materials from section 3 to free fluid vortex above it. In order to create a two-dimensional fluid vortex, we consider a steady low-Reynolds number lid-driven cavity problem, consisting of a free fluid domain $\Omega_{f}$ and a poroelastic domain $\Omega_{p}$. We will validate and characterize the accuracy of the effective continuum description by comparing to a second approach, in which the whole domain is meshed and the flow field as well as the displacement field are resolved at all spatial scales. For the fully resolved numerical studies to be feasible, we do not deform the computational mesh when the microstructure is displaced, which sets an upper limit of displacement we can consider to roughly $\boldsymbol{\mathrm{v}}\lesssim 0.1\,l$. Note that this means we consider a one-way interaction problem, that is, the material elasticity does not influence the free fluid, whereas the free-fluid does induce a displacement of the material. The results obtained through this simplification provides new fundamental insight into the physics at the interface between the free-fluid and porous region. It allows us to study how much of the shear stress from the free fluid results in a stronger pore flow and how much of it is borne by the solid. The theory in sections 2 and 3 and the numerical implementation of effective equations[46] is valid for two-way coupled problems, but the comparison with fully resolved simulations of such systems we leave for future work.

4.1 Effective continuum description

As shown in Fig. 6$a$, the two-dimensional (2D) lid-driven cavity has a depth of $H+d$, a length of $H$ and is infinitely wide. The poroelastic medium is confined to $z\leq 0$ and $-H/2\leq x\leq H/2$. The flow is driven by the top-wall, which moves in the $x$ direction with speed $U_{w}$. In $\Omega_{p}$, we solve the effective equations for a poroelastic medium (5, 8) with the boundary conditions at the wall (9) and at the interface with the free fluid (10,11). For the free fluid in $\Omega_{f}$, we solve equations (1, 2) with velocity boundary condition at the interface with the poroelastic material (12, 13). At the side and top walls of the cavity, we impose

[TABLE]

respectively. The governing equations are discretized using finite-element-method (FEM). The simulation domain along with equations is defined in the FreeFEM++[48] software using uniform mesh spacing666We have investigated results using half of the mesh spacing $\Delta s=0.10\,l$ and observed that the slip velocity changed by $0.4\%$ and the horizontal displacement near the interface changed by $0.8\%$ at the center of the cavity. $\Delta s=0.20\,l$. We implement the weak formulation and solve the time dependent problem in a fully implicit manner, where boundary conditions at the interface are enforced using Lagrange multipliers777Exact weak form can be found in documentation of the open-source software[46].. The effective coefficients, such as elasticity tensor and permeability matrix, are taken from section 3. To render the validation with the resolved model feasible, we select a moderate scale separation parameter $l/H=0.1$.

4.1.1 Effective fluid velocity

We start by presenting the flow field $\boldsymbol{\mathrm{u}}=\left(u_{x},u_{z}\right)$ in the full domain ($\Omega_{p}$ and $\Omega_{f}$) for the cubic-symmetric poroelastic medium with porosity $\theta=0.85$. Streamlines and isocountours of the streamwise fluid velocity are shown in Fig. 7$a$. In the free fluid domain we normalize the results using the top wall velocity $U_{w}$ (in the simulations we have set $U_{w}=100\ U^{d}$), whereas in the poroelastic domain we normalize the results using Darcy’s velocity $U^{d}\equiv l^{2}\Delta P/\left(\mu H\right)$. The moving top wall creates a circulation in the cavity, where the flow for $x>0$ is directed downwards, and upwards for $x<0$. We observe that due to the vortex in the free fluid, there exists a transfer of mass and momentum across the interface. This is characterized in more detail from the slip and infiltration (penetration) velocities very close to the interface ($z=0.01H$) in Figs. 8$a$ and  8$b$, respectively. The streamwise slip velocity (Fig. 8$a$) has a parabolic shape near the poroelastic medium with its minimum (largest magnitude) velocity at the center of the cavity. This velocity component is mainly created by the shear of the free fluid, which is the strongest at $x=0$, before it gradually decays when approaching the sides of the cavity. The penetration velocity (Fig. 8$b$) shows a macroscopic behaviour similar to a sine wave. We observe that for $x>0$, there is a net mass/momentum transport from free fluid region to the poroelastic region, whereas for $x<0$ the net mass flow is in the opposite direction.

Fig. 7$a$ shows that the flow inside the bed circulates in the same direction as in the free-fluid region. Note that the streamwise velocities in the two domains differ by roughly three orders of magnitude. This gives the impression of a clear discontinuity of streamwise velocity component; Figs. 7$b$ and 7$c$ shows the streamwise $u_{x}$ and wall-normal $u_{z}$ velocity profiles for the fixed streamwise position $x=0.1H$. For the streamwise component (see inset of Fig. 7$b$), one can observe that there is a very slow Darcy flow inside the poroelastic material, whereas near the interface there is a jump to the fast slip velocity of the free flow. The leading-order effective equation for flow inside a poroelastic material, which is relative Darcy’s law (5), contains only the pressure contribution from the flow and no viscous fluid stress. As a consequence, the fluid inside the poroelastic layer can only be driven by the normal stress (can not respond to shear stress), which is the reason why the velocity jump arises. In Fig. 7$c$, we see that wall-normal component $u_{z}$ is continuous across the interface, which follows from mass conservation.

We compared the homogenized model results between cubic-symmetric poroelastic material ($\theta=0.85$) and monoclinic-symmetric poroelastic material ($\theta=0.86$) and did not observe any significant difference in flow velocities, despite the different levels of anisotropy in the pore microgeometries. The reason for this outcome is the fact that the introduced anisotropy results only in higher order corrections of permeability tensors (24,26, 32) and Navier-slip tensors (36).

4.1.2 Effective solid displacement

In Figs. 9$a$ and 9$b$ we show displacement $v_{x}$ and $v_{z}$ along $x$-coordinate at a fixed $z=-0.05\,H$ for both cubic-symmetric and monoclinic-symmetric pore-scale geometries. The displacement fields are normalized888Note that due to the linearity of governing equations, the results are valid for any chosen Young’s modulus, as long as the deformation is small enough ($\boldsymbol{\mathrm{v}}\lesssim 0.1\,l$) to neglect the mesh deformation and Poisson’s ratio is $\nu=0.3$. The dimensional normalization is obtained in appendix A. using length scale $H\Delta P/E$. Moreover, the $z$-coordinate is chosen to correspond to the centres of connecting cylinders. In this way, the result can easily be compared to fully micro-resolved simulations (see next section), for which it is not straight-forward to define displacement fields in the fluid part of the pores. For both microstructure geometries, the horizontal displacement $v_{x}$ near the interface (Fig. 9$a$) reveals a similar behaviour as for the slip velocity, i.e., the displacement has parabolic shape with maximum magnitude at the center of the cavity. The horizontal displacement is in the direction of the slip velocity directly above the interface and is caused by the shear stress induced by the overlying flow vortex in the cavity. The vertical displacement $v_{z}$ near the interface (Fig. 9$b$), on the other hand, shows a very similar behaviour to that of the penetration velocity, i.e., the displacement has a sine-like shape. We thus observe that the solid displacement is complying to the tangential and normal fluid velocities near the interface.

In Figs. 9$c$ and 9$d$ we show $v_{x}$ and $v_{z}$ along the $z$-coordinate at a fixed position $x=0.15\,H$. For the horizontal component $v_{x}$ (Fig. 9$c$) two regions can be identified, where the displacement decreases at different rates. The region down to $z\gtrsim-0.2\,H$, for which the decay is very fast, is determined by the shear stress at the interface. Below $z\lesssim-0.2\,H$, it is the slow Darcy flow that induces the small displacement. The vertical displacement (Fig. 9$d$), on the other hand, is entirely governed by the slow penetration velocity inside the medium or the interface-normal stress at the interface, depending on which one is the dominating effect.

We note that the horizontal displacement (Fig. 9$a$) has slightly larger magnitude and is skewed for the monoclinic compared to the cubic material. This difference can be explained by a slightly smaller compression elasticity coefficient $C_{xxxx}$ ($4.668<4.792$) and by the existence of an elasticity coefficient $C_{xzxx}$ for the monoclinic material (see section 33.13.1.3), which relates the strain of $v_{x}$ in the $x$ direction with the stress in the interface-normal direction $z$ on the plane with a fixed $x$ coordinate. The vertical displacement (Fig. 9$b$) is also larger in magnitude and skewed for the monoclinic geometry. The larger amplitude can again be explained by the difference in compression elasticity coefficient $C_{zzzz}$ ($4.031<4.792$), whereas the skewness can be attributed to $C_{xxxz}$, which provides a coupling between the interface-normal stress of the free fluid to strain of the vertical displacement $v_{z}$ in the $x$ direction. We thus conclude that an anisotropy in the microstructure breaks the symmetric behaviour of the solid displacement in the macroscale due to the additional cross terms (such as $C_{xxxz}$) in the effective elasticity tensor.

4.2 Validation to fully microscale-resolved simulations

In order to validate the effective model, a three-dimensional (3D) lid-driven cavity problem is considered (same depth of $H+d$ and length $H$ as before), where the flow is periodic in $y$ direction over length $l$, which is set by the pore-sctructure. Therefore, it is sufficient to consider the domain shown in Fig. 6$b$ (spanning over one microstructure in $y$ direction) for direct numerical simulations (DNS).

The domain shown in Fig. 6$b$ is defined using GMSH software [47] and meshed using spacing $\Delta s_{1}=0.04\,l$ at the solid skeleton and $\Delta s_{2}=0.16\,l$ at the boundaries of the cavity999We have carried out same simulations using mesh spacings $\Delta s_{1}=0.05\,l$ and $\Delta s_{2}=0.20\,l$. The change in fluid slip velocity was everywhere below $0.3\%$. The change in horizontal solid displacement was everywhere below $5.5\%$, with maximum at the center of the interface.. We import the generated mesh into FreeFEM++ [48], in which we define and solve the governing equations (1–3) with boundary conditions (4, 37). At the walls, we use no-slip velocity for fluid and zero displacement for solid structure. To simplify the numerical task, we assume a steady flow and neglect inertial effects, which simplifies the Navier-Stokes equations (1–2) to the linear Stokes equations. We use Taylor-Hood (P2+P1) elements for the Stokes system, and quadratic elements (P2) for the solid skeleton linear elasticity system (3). The resulting linear algebraic system is solved using a GMRES iterative linear solver up to a relative tolerance $e=10^{-10}$. For the solid skeleton, we choose the isotropic material that was used as the starting point for the computation of effective properties in section 3 (Poisson’s ratio $\nu=0.3$ and Young’s modulus $E$ unspecified). We assume that deformations are small enough such that the computational mesh can be kept static; this limits the range of length scales $\bar{l}=H\Delta P/E$ – based on Young’s modulus $E$, large scale length $H$ and pressure difference $\Delta P$ – that can be considered.

4.2.1 Fully-resolved fluid flow

In Fig. 10$a$, we show the resolved stream-wise velocity (slip velocity) near the interface at a distance $z=0.01\,H$ and coordinates $y=0$ and $y=0.49\,l$. The poroelastic bed in this case is built using cubic-symmetric geometry. We observe that the slip velocity is slightly slower at the center of the cavity slice ($y=0$) compared to the $xz$-plane near the periodic boundaries ($y=0.49\,l$). This is because the bulk of the solid material (sphere) is located at the center of each volume element, therefore the surrounding fluid is slowed down more than the fluid near the cylindrical obstructions close to the periodic boundaries. The stream-wise velocity distribution in the $y$ direction at three $x$ and $z$ coordinates is shown in Fig. 10$b$. Moving away from the poroelastic material (increasing $z$ coordinate) leads to a rapid dissipation of velocity variation; that is, the flow velocity approaches a constant value with respect to the $y$-coordinate. We continue by averaging the DNS quantities in the $y$-direction as,

[TABLE]

where $f$ denotes any of $u_{x}$, $u_{y}$, $u_{z}$, $p$, $v_{x}$, $v_{y}$ and $v_{z}$. The function $I_{d}$ is an indicator function, which is $I_{d}=1$ in the domain, where field $f$ is defined, and $I_{d}=0$ elsewhere. This average is also known as intrinsic average[31], because it is normalized by fluid or solid volumes separately. For convenience, we omit the “bar” notation further on.

The average slip velocity and penetration velocity variations over the $x$ coordinate sampled at $z=0.01\,H$ are shown in Figs. 8$a$ and  8$b$, respectively, together with the effective model curves; one can conclude that the macroscale model is accurate. Note that the effective model is fully non-empirical without any fitting parameters. The penetration velocity (Fig. 8$b$), shows micro-scale oscillations over one pore-scale structure in the resolved simulations, which are by construction not captured by the effective equations. The penetration velocity is somewhat under-predicted by the effective model, which may be an indication of a pressure jump [43, 44], not modelled in the current work.

4.2.2 Fully-resolved solid displacement

Fig. 11$a$ shows the deformation of the cubic-symmetric structure due to the free fluid vortex above. The figure serves only as an illustration, where the solid structure is displaced after the computation, since our implementation is restricted to static meshes. The solid displacement obtained from DNS is compared to the effective model over a horizontal slice at $z=-0.05\,H=-0.5\,l$ in Figs. 11$b$ and 11$c$. Similarly as for the fluid velocity, we observe that the microscale features of the displacement field are not captured by the leading-order effective model, but that the macroscale behavior is in good agreement.

Fig. 12$a$ compares the horizontal displacement along the $z$-coordinate. This shows that the effective model overestimates the horizontal displacement close to the interface, but is very accurate below $z\lesssim 0.05\,H$. This indicates that the inaccuracy is caused by the interface stress boundary condition (11). The vertical displacement (Fig. 12$b$) is overestimated over the entire depth of the cavity, which suggests that the vertical displacement is mostly governed by the stress at the interface between the free fluid and the poroelastic medium. Despite the inaccuracy at the interface, the model captures the essential qualitative features of the displacement field, including the two different regions of decay of the horizontal displacement. The seemingly good agreement of the horizontal displacement in Fig. 11$b$ arises due to the micro-scale variations of the DNS solution, as seen in Fig. 12$a$

Our results show that the free-fluid shear stress transfer to the poroelastic material shear stress is a macroscopic phenomenon. The influence of the interface stress is ranging over a number of pore-structures; the shear region in Fig. 12$a$ extends over distance of $z\approx 0.3\,H$ containing $2$ to $3$ unit-cell structures. However, the free-fluid shear stress transfer to pore fluid shear stress is a microscopic phenomenon. The influence of the stress is ranging over less than one pore-structure; the region over which velocity decays to Darcy’s velocity in Fig. 7$b$ extends over distance of z\approx 0.05\,H\ {\color[rgb]{0,0,0}=0.5\,l} containing only half of one unit-cell structure. Over this short distance, all the fluid shear stress is transferred over to the solid skeleton, thus contributing to large shear region in Fig. 12$a$.

As a final comment, the effective properties ($\boldsymbol{\mathrm{{C}}},\boldsymbol{\mathrm{{\alpha}}}$ and $\mathcal{E}$) of the interior can be used at the interface for the stress boundary condition (11). The reason is that the shear stress from the free fluid is eventually borne only by solid skeleton. As illustrated in Fig. 13, left, from the microscopic – or “reality” – point of view, the free fluid shear stress is transferred to both solid skeleton and pore fluid. However, roughly a half of a microstructure below the interface, the shear stress of the pore fluid is transferred to the surrounding solid skeleton via viscous friction. Therefore the elasticity tensor from the interior is a reasonable estimate for the elasticity tensor at the interface. In other words, the model (Fig. 13, right), in which free fluid shear stress is transferred directly to solid skeleton only, is a reasonable approximation. This is, however, not the case for the velocity boundary conditions (12) and (13), for which the micro-scale viscous dissipation has to be modelled by Navier-slip term ($\boldsymbol{\mathrm{{L}}}$); otherwise, as shown by Lācis & Bagheri[42, tab. 2], the predicted interface velocity would be around $1/\epsilon$ times smaller compared to DNS. This is due to the fact that the pore fluid undergoes a rapid acceleration caused by the free fluid shear at the interface, which renders the interior velocity a very inaccurate estimate for the interface velocity.

4.3 Capturing anisotropic effects with the effective model

The displacement fields of the two geometries (cubic-symmetric and monoclinic-symmetric) obtained from DNS are compared in Fig. 14. For the monoclinic-symmetric geometry the horizontal displacement is roughly $6\%$ larger, and the vertical displacement is roughly $10\%$ larger (compared to cubic-symmetric geometry). This difference cannot be explained by the change of porosity only (around $1\%$), therefore the reason must be the introduced anisotropy of the pore geometry. Moreover, in the horizontal and the vertical displacement fields one can observe a small asymmetry between left and right half of the cavity, which is a sign of symmetry breaking due to anisotropic effects.

In section 44.14.1.2, a difference was also observed between the two microstructures in the effective model, which was explained by comparing the entries of the effective elasticity tensor. Fig. 14, obtained using fully resolved simulations, can be compared with the corresponding figure obtained from the effective model (Fig. 9). We observe that the relative change between those two geometries is accurately captured by the effective model. This indicates that using the interior effective elasticity, when imposing stress boundary conditions at the interface, provides an effective model that is practically useful in capturing differences in the macroscale, caused by anisotropy within the microstructure.

5 Discussion

The numerical framework presented here is based on first solving a set of microscale problems from section 3 to compute the effective tensors ($\boldsymbol{\mathrm{{C}}}$, $\boldsymbol{\mathrm{{K}}}$, $\boldsymbol{\mathrm{{\alpha}}}$, $\boldsymbol{\mathrm{{K}}}^{\mathrm{if}}$ and $\boldsymbol{\mathrm{{L}}}$), followed by solving the macroscopic effective equations coupled through interface conditions to the free fluid (section 2). The equations formed in both these steps stem from multi-scale homogenization approach and rest on a number of assumptions, similarly as explained in literature[51, 52] for the method of volume averaging. These assumptions restrict the range of physical parameters that we are able to model using effective continuum theory. It is therefore essential to understand if the admissible range of parameters are relevant for solving physical problems that arise in nature and engineering. In this section we discuss the necessary assumptions used in the development of the current framework and how they relate to physical constraints.

The first and foremost prerequisite is scale separation. Based on the numerical tests in the literature[42] and section 4, we have determined the practical limit on the scale separation parameter

[TABLE]

which is less restrictive than the asymptotic limit $l/H\rightarrow 0$. In other words, one does not need to be in the asymptotic limit in order to apply the current framework. We compared the effective model results with predictions from fully resolved simulations in section 4 and observed good agreement with respect to the flow and the displacement fields for the scale separation parameter $l/H=0.1$. It is interesting to point out that Auriault[53] has shown the next order corrections to the Darcy’s law to be zero for macroscopically homogeneous porous media, therefore the Darcy’s law holds well also in the case of poor scale separation. This would explain the good agreement between the model and DNS for the flow velocity. It is possible that similar conclusion could be drawn also from the correctors of the elasticity problem.

We assume that the pore Reynolds number is smaller than one, i.e.

[TABLE]

Here, recall that the $U^{d}\equiv l^{2}\Delta P/\left(\mu H\right)$ is the definition of seepage velocity. The above assumption is a good one for poroelastic materials that are densely packed, resulting in a slow flow through the pores, which can be described using steady linear Stokes equations – such as those formulated for $\boldsymbol{\mathrm{{K}}}$, $\boldsymbol{\mathrm{{K}}}^{\mathrm{if}}$ and $\boldsymbol{\mathrm{{L}}}$ in equations (22–23), (28–31) and (33–35). For larger pore Reynolds numbers the microscopic problems will become non-linear, which renders the current multi-scale approach, based on linear decomposition, unfeasible. A possible workaround is to use some kind of linearisation, similarly as done by Zampogna & Bottaro[38].

Next, we make an assumption on how the macroscopic free fluid time scale $\Delta\tau$ that force the poroelastic bed is related to the microscopic time scale $l/U^{d}$ inside the bed. Specifically, the frequency $1/\Delta\tau$, at which the free fluid interacts with the poroelastic bed, has to satisfy,

[TABLE]

This essentially states that the changes in the free fluid must be slower compared to the time a fluid parcel needs to travel a pore length $l$. Thus, from the microscale viewpoint, the external macroscopic forcing is slow, which in turn has the consequence that the corresponding microscale problems (22–23), (28–31) and (33–35) are steady. If external forcing frequency is higher than $U_{d}/l$, the Stokes problems (22–23), (28–31) and (33–35) will become time-dependent, which in turn would require solving convolution integrals in order to take into account the time history, see Mei & Vernescu[12, eq. 6.6.11].

Another restriction of the current method is on the relative size of the characteristic normal stresses of the flow and the solid skeleton. That is, we assume that the macroscale global $\Delta P$ relates to the characteristic Young modulus of the solid skeleton as follows,

[TABLE]

This assumption holds in many engineering configurations, where elastic properties of materials are often of order MPa or GPa, while the pressure difference generated by moving fluids is commonly of order kPa (in wind tunnel experiments, for example). It also holds for many biological systems; the elasticity moduli are of order MPa or kPa[54, 55, 56] where biological materials are often exposed to much slower fluid flows[57] and consequently smaller fluid forces.

The final requirement is that inertial effects of the solid skeleton at the microscopic level are small. We expect that the pore-flow is sufficiently slow such that inertial effects of the solid skeleton in the microscale are not excited. This implies a relation between the characteristic macroscopic pressure, the solid density and characteristic time scale:

[TABLE]

This results in linear, time-independent microscale solid test problems (16–17) and (19–20). If this restriction is not obeyed, the test problems should be complemented using pore-scale inertial effects.

Note that the above assumptions are not unique and in order to understand mathematically why these particular choices are made, we refer the reader to appendix A. At the end, we want to understand whether the resulting effective equations governing the poroelastic bed can be used to describe the length and time scales that are physically relevant. In order to do so, we compare restrictions, involving the external macroscopic forcing time scale $\Delta\tau$ (41,43), with the intrinsic poroelastic time scale $\Delta\tau_{\mathrm{p}}$ and the time scale $\Delta\tau_{\mathrm{i}}$ related to the speed of waves in the effective bed. The former time-scale $\Delta\tau_{\mathrm{p}}$ characterizes the time for the pressure field to equilibrate via fluid transport in the medium and therefore determines how fast the poroelastic bed can respond to external forcing [58]. The latter time scale $\Delta\tau_{\mathrm{i}}$ characterizes the time for wave propagation, and therefore also the time it takes for information to propagate through the effective bed. Following Skotheim & Mahadevan[58], the time scales are given by

[TABLE]

where $k$ is the characteristic permeability and $C_{\mathrm{eff}}$ is the characteristic effective elasticity coefficient. These time scales can be compared to the assumptions discussed above using estimates $k\sim l^{2}$ and $C_{\mathrm{eff}}\sim\theta_{s}E$, where $\theta_{s}=1-\theta$ is the solid volume fraction. Using assumptions (41) and (42) from the current work, we then arrive to inequality $\Delta\tau\geq\theta_{s}\Delta\tau_{\mathrm{p}}$. This shows that the macroscopic driving force time scale can be of the same order as the poroelastic material time scale, therefore the model derived in the current work allows for description of poroelastic effects in the medium. Furthermore, using assumptions (42) and (43), we arrive to inequality $\Delta\tau\geq\sqrt{\theta_{s}}\Delta\tau_{\mathrm{i}}$. This indicates that the model explained in this work allows for description of travelling waves through poroelastic medium. To sum up, the current multi-scale approach should allow for the description of problems where $\Delta\tau$ is of the order of $\Delta\tau_{p}$ (to capture strong fluid-elasticity interaction [58]) and also of the order of $\Delta\tau_{i}$ (to capture travelling waves or elastic instabilities).

6 Conclusions

We have considered the problem of a free-fluid interacting with a poroelastic bed, by deriving and validating an effective continuum model for the bed and its interface with the above free fluid. Although the effective equations of the interior of the bed are well-established, their coupling to a non-trivial vortical free fluid through a set of interface conditions have not been considered and validated before. The imposed interface conditions are (i) the velocity boundary condition for the free fluid; (ii) the pressure continuity boundary condition for the pore-pressure; and (iii) the stress continuity boundary condition for the displacements of poroelastic media. The first two conditions have been derived from first-principles using certain assumptions, while the stress boundary condition is postulated a posteriori. In particular, for condition (iii) the interior effective parameters ($\boldsymbol{\mathrm{{C}}}$ and $\boldsymbol{\mathrm{{\alpha}}}$) were used, whereas for condition (i) the interface permeability and the Navier-slip tensor were computed. This asymmetry in the treatment of the boundary condition is motivated by the fact that for the velocity boundary condition – which is needed to solve for the free-fluid – the transfer of shear-stress to pore flow requires a new Navier-slip tensor $\boldsymbol{\mathrm{{L}}}$, because Darcy’s law cannot accommodate any shear. Therefore, the interfacial velocity condition requires a special treatment. For the stress boundary condition – which is needed to solve for the displacement in the bed – both normal and tangential stresses induced by the free flow above can be matched by the corresponding stresses of the solid. We have shown that using condition (iii), based on the interior coefficients, provides a satisfactory model; on the one hand, it does result in a – small, but not insignificant – discrepancy with the fully-microscopic simulations that we use as a validation; on the other hand, it captures the effects of small changes in the microstructure anisotropy correctly and predicts the overall behaviour in a physically consistent and controllable manner. We thus believe that this approach for modelling the interaction of poroelastic beds with freely moving fluids is a viable framework for engineers. The practical limits of the derived model have been discussed, where we show that the proposed model can be employed for any physically relevant poroelastic material. The corresponding codes of the numerical implementation used in the present work have been released as an open-source software[46].

In the future, we will further improve the model by treating the stress boundary condition using appropriate interface-cells. We also want to validate the proposed model for unsteady flows, where inertia is not negligible. Finally, the next step for the model development is to extend it to significant surface deformations, that would require moving the interface and porosity variation in space.

\dataccess

The core of the codes used for producing the plots in the current paper are available as open source software in Github repository[46]. The same repository contains the flow and displacement fields of DNS simulations presented in this paper.

\aucontribute

The idea of the present work was conceived through discussions between all authors. U.L. has derived the homogenized equations with feedback from all other authors. The numerical implementation of the homogenized and microscale equations was carried out by U.L. The reported numerical simulations were performed by U.L. with feedback from S.B. The results from microscale solvers were independently verified by G.A.Z. using OpenFOAM for cubic-symmetric material. Results were analysed and paper was written by U.L. and S.B. with feedback from G.A.Z.

\competing

Authors declare no competing interests.

\funding

The work by U.L and S.B. on this topic has been funded by the Swedish Research Council (VR-2014-5680) and the Göran Gustafsson foundation.

\ack

G.A.Z. acknowledges several insightful interactions with Prof. A. Bottaro.

\disclaimer

The views and opinions expressed in this article are those of the authors.

Appendix A Governing micro-scale equations and assumptions

In the current appendix, we explain the non-dimensionalization of the pore-scale governing equations, as well as required assumptions and scale estimates for the derivation of the macroscale governing equations and accompanying pore-scale test problems. The governing equations, presented in main paper (1–4), are summarized here for convenience:

[TABLE]

where we have used $\Sigma_{ij}$ and $\sigma_{ij}$ to denote fluid and solid stress tensors, respectively. We have also introduced strain rate tensor for the fluid $\varepsilon_{ij}\left(u\right)=0.5\left(u_{i,j}+u_{j,i}\right)$ and strain tensor for the solid $\varepsilon_{ij}\left(v\right)=0.5\left(v_{i,j}+v_{j,i}\right)$. Here $\Omega_{f}$ is the fluid domain (both free fluid and fluid between solid skeleton – pore fluid), $\Omega_{s}$ is the solid skeleton domain and $\Gamma_{s}$ is boundary between solid and fluid. For convenience, we use the index notation, where summation is carried out over repeating indices and comma indicates a derivative. Now we assume that there is scale separation between two length scales, $l$ and $H$, with $\epsilon=l/H\ll 1$. Other parameters in this problem are $\rho_{\mathrm{f}}$, $\rho_{\mathrm{s}}$, $\Delta P$, $\mu$, $\Delta\tau$ and $\boldsymbol{\mathrm{{C}}}^{\mathrm{sk}}$, where $\Delta P$ is a characteristic pressure difference, $\Delta\tau$ is a characteristic time scale of the processes being studied and $\boldsymbol{\mathrm{{C}}}^{\mathrm{sk}}$ is linear solid skeleton elasticity tensor. Based on these parameters, we can define a velocity scale $U^{d}\equiv l^{2}\Delta P/\left(\mu H\right)$ arising from momentum balance at the microscale[42, eq. A 1]. We estimate the scales of flow, pressure, displacement fields and derivatives in the interior as

[TABLE]

Based on these estimates, we choose to render the equations dimensionless (we temporarily use “tilde” to denote dimensional quantities), using the relationships

[TABLE]

Note that non-dimensionalization using physical parameters $\Delta P$, $\mu$, $l$, $\Delta\tau$ and $H$ is one particular choice and different options are possible. With this choice we have the dimensionless order of all terms close to unity

[TABLE]

These orders are required later on, when the multi-scale expansion is carried out. Under the chosen normalization (50), the governing equations become

[TABLE]

This system of non-dimensional equations is valid everywhere and so far there are no additional assumptions employed. The total list of physical parameters in the problem is $\rho_{\mathrm{f}}$, $\rho_{\mathrm{s}}$, $\Delta P$, $\mu$, $l$, $H$, $\Delta\tau$, and $\boldsymbol{\mathrm{{C}}}^{\mathrm{sk}}$. This adds up to $7+n$ parameters, where $n$ depends on the properties of elastic skeleton, and for the current work $n=2$ (we assume that it is isotropic and elasticity is characterized using $2$ parameters). According to Buckingham $\pi$ theorem, we then have to have 6 dimensionless parameters, which are

[TABLE]

Here, we have separated a Young modulus out of the skeleton elasticity and the dimensionless tensor $C^{\mathrm{sk}}_{ijkl}$ is characterized by only one scalar – Poisson’s ratio $\nu$. If the skeleton is built from an anisotropic material, then the dimensionless tensor $C^{\mathrm{sk}}_{ijkl}$ would imply a larger number of non-dimensional parameters. For the isotropic case, this form is useful to describe behaviour of isotropic materials with similar Poisson’s ratio but different Young’s modulus. In addition, since the resulting equation system in the Stokes limit is linear, in the main paper we plot displacements normalized with dimensionless elasticity parameter $v_{i}\cdot\bar{\bar{E}}$, which in dimensional setting becomes $\tilde{v}_{i}\cdot E/\left(H\Delta P\right)$. In other words, the values presented in the displacement plots would correspond to simulation with $\bar{\bar{E}}=1$.

To be able to use the multi-scale expansion, one has to estimate the relative orders of all the terms. We assume the following:

[TABLE]

Now the relative magnitude between different terms is completely illustrated by the $\epsilon$ pre-factors in the governing non-dimensional equations (52–56) in the interior of the poroelastic material. This, however, does not hold in the free fluid. In other words, in the free fluid region, it is not possible to use the Darcy’s law as a governing equation.

In the main paper section 5 introduced restrictions of this theory has been obtained by relaxing the equality conditions from (58). By “relaxing” we mean that the theory should not only be applicable to one value of each of the dimensionless parameter, as set by (58), but for a range of dimensionless parameters. The generalization from the equality assumptions (58) to the inequality restrictions is done as fallows:

The equality assumption $\bar{\bar{E}}=\mathcal{O}\left(1\right)$ is relaxed to the inequality restriction based on numerical tests in the main paper section 4, in which one could see that the theory works also for very stiff materials $\bar{\bar{E}}\gg\mathcal{O}\left(1\right)$. Therefore we conclude that the theory works for a range of dimensionless parameters $\bar{\bar{E}}\geq\mathcal{O}\left(1\right)$, which includes the theoretical assumption $\bar{\bar{E}}=\mathcal{O}\left(1\right)$ and numerical validation at $\bar{\bar{E}}\gg\mathcal{O}\left(1\right)$. 2. 2.

The equality assumption $\bar{\bar{\rho}}=\mathcal{O}\left(1\right)$ is relaxed to inequality restriction $\bar{\bar{\rho}}\leq\mathcal{O}\left(1\right)$ based on physical understanding that the developed method works also in the non-inertial regime; that is, when the solid density is sufficiently small that inertial effects can be neglected, the parameter $\bar{\bar{\rho}}$ can be set to zero. 3. 3.

Finally, the equality assumption $\bar{\bar{f}}=\mathcal{O}\left(1\right)$ is relaxed to inequality restriction $\bar{\bar{f}}\leq\mathcal{O}\left(1\right)$ based on numerical tests in the main paper section 4. There we have validated the model equations for an “infinitely slow” or steady test problem, for which $\bar{\bar{f}}=0$. Uniting the two parameter limits, for which the theory should work, we get $\bar{\bar{f}}\leq\mathcal{O}\left(1\right)$.

Although the equality assumptions on the dimensionless parameters (58) are strict and in principle one could set them to the specific value in the equations, we retain these coefficients in the equations to facilitate tests away from the outlined assumptions.

Appendix B Derivation of homogenized effective equations

We are now ready to derive the homogenized governing equations for the poroelastic material in macroscale. Additionally, we obtain the microscale test problems for determining the effective material properties. We introduce the macroscale and microscale coordinates

[TABLE]

respectively. These coordinates are appropriate to describe the macroscopic and microscopic variations and are related to each other by $X_{i}=\epsilon x_{i}$. In the new coordinates, the spatial derivative is given by

[TABLE]

where $\left(\right)_{,i_{0}}$ denotes the derivative with respect to $X_{i}$ and $\left(\right)_{,i_{1}}$ with respect to $x_{i}$. Then we introduce a multi-scale expansions as

[TABLE]

which we insert into the main equations (52–56). The first two orders of fluid momentum equation (52) after the expansion with corresponding stress tensors are

[TABLE]

where we have introduced the strain rate tensor in the micro-scale $\varepsilon^{1}_{ij}\left(u\right)=0.5\left(u_{i,j_{1}}+u_{j,i_{1}}\right)$. The first two orders of the fluid continuity equation (53) after the expansion are

[TABLE]

The first three orders of the solid displacement equation (56) after the expansion with corresponding stress tensors are

[TABLE]

Here we have used double-dot to denote second derivative in time $\ddot{v}_{i}=\partial_{t}^{2}v_{i}$. Note that the orders of the solid equations are shifted to one order lower value to match the orders of the stress boundary condition between solid structure and pore geometry (55). The boundary conditions (54–55) for all orders $n\in\left(0,\infty\right)$ take the form

[TABLE]

where the time derivative is denoted using dot notation $\dot{v}_{i}=\partial_{t}v_{i}$ and there are no lower order values from the fluid stresses $\sigma_{ij}^{(-1)}\hat{n}_{j}=\Sigma_{ij}^{(-1)}\hat{n}_{j}=0$.

Solving $\mathcal{O}\left(1\right)$-problem for fluid and $\mathcal{O}\left(\epsilon^{-1}\right)$-problem for elasticity, gives us

[TABLE]

that is, the leading order pressure and displacement depend only on the macroscale. The fluid $\mathcal{O}\left(\epsilon\right)$-problem is

[TABLE]

This problem can be solved using the ansatz

[TABLE]

where the second term in the velocity ansatz is added to satisfy the velocity boundary condition at the surface with the solid skeleton. In principle, one could introduce a second proportionality tensor, but the resulting problem would always be trivial to solve, yielding factor $\delta_{ij}$ before the velocity of the solid skeleton. Exactly the same conclusion holds also for the derivation of the velocity boundary condition [42], therefore the time derivative of the displacement appears in the boundary condition in the main paper (12) and (13). The solid $\mathcal{O}\left(1\right)$-problem is

[TABLE]

This problem can be solved using the ansatz

[TABLE]

where division by $\bar{\bar{E}}$ is done in order to arrive with test problem for $\eta_{i}$, which is independent from the dimensionless parameter $\bar{\bar{E}}$.

Inserting these ansatzes into the governing equations, one can group different terms together and form a particular solution by setting individual group contributions to zero. For the fluid problem, only one Stokes system has to be solved as a closure problem, which is

[TABLE]

where $\Gamma$ is the boundary with solid skeleton. This problem is solved in one unit-cell using periodic boundary conditions at all sides for both $K_{ij}$ and $A_{j}$ fields. This test problem is explained in the main paper near equations (22–23) in dimensional setting. For the elasticity, two quantities have to be solved for. The first one is for the displacement tensor $\chi_{ijk}$, which is governed by

[TABLE]

where elasticity $\bar{\bar{E}}$ has been cancelled out as a common factor. In this case, the displacement field is also exposed to periodic boundary conditions at all unit-cell sides. In order to ensure uniqueness of the solution, one can constrain the average values of the displacement tensors or enforce point constraints. Finally, there is a problem for proportionality vector before pressure

[TABLE]

where $\bar{\bar{E}}$ has been cancelled in the product of the elasticity tensor with the pre-factor $\eta/\bar{\bar{E}}$ in the displacement ansatz. Also this problem needs to be fixed using periodic boundary conditions and constraints on average values. The elasticity test problems are explained in the main paper near equations (16–16) and (19–20) in dimensional setting.

To derive the governing equations, the volume average operator, as defined in the main paper, equation (15), is used. We employ the volume averaging on the ansatz of the fluid velocity and arrive to

[TABLE]

where $\theta=V_{f}/l^{3}$ is the volume fraction of the fluid (porosity). Note that $\dot{v}_{i}^{(0)}$ is independent of the microscale, therefore averaging only results in $\theta$ pre-factor. In order to continue derivations, we express the volume average of the first order velocity as

[TABLE]

where the closed surface integral has been divided into a part over the boundary between solid and fluid $\Gamma$ and a part over the periodic boundary $\Pi$. Due to periodicity, the second term is zero, and the first one can be rewritten using boundary condition (46) at the solid skeleton as

[TABLE]

where we have added an integral of zero value, taken over the periodic boundaries of the solid. We have also transferred between adjacent domains with $\hat{n}^{s}_{i}=-\hat{n}^{f}_{i}$. Hence in this example we clearly see that the application of Gauss theorem in the microstructure relies on the assumption of periodicity. Nevertheless, the resulting model could also be used if the periodicity assumption is satisfied only approximately. Using the linear expression for the first order displacement (B 80), we arrive with the continuity equations for the macroscale

[TABLE]

where we have used the assumption that the porosity is uniform, and thus the derivative and volume averaging commute. To obtain the final equation in the macroscale, we introduce a stress tensor, which is defined over the whole volume as

[TABLE]

which can be differentiated with respect to the macroscale, and the result is

[TABLE]

Now we average the divergence of the total stress tensor

[TABLE]

where we have used equations (B 64) and (B 69) to rewrite the divergence of zeroth order stress tensors. In the same way as in equations (B 89 – B 90), we can use Gauss theorem, periodicity and the boundary conditions, to write:

[TABLE]

On the other hand, the zeroth order stress tensors can be averaged using explicit expressions, for the fluid stress we get

[TABLE]

whereas for the solid stress we find

[TABLE]

To sum up, the final equation in the macroscale is

[TABLE]

where the effective elasticity tensor is

[TABLE]

and the tensor in front of the pore pressure is

[TABLE]

The equations to be solved for the homogenized macroscale model expressed using relative Darcy’s law (B 88) in the conservation of mass law (B 91) are

[TABLE]

where all variables and coordinates are now dimensional. These are the equations explained in the main paper (5) and (8). We note that if the solid density is very small, as allowed by the restriction (43) in the main paper, then the inertial term in solid displacement field can be neglected. Here we have also used the equality

[TABLE]

which has been derived by Mei & Vernescu[12, p. 258–259] and confirmed within the current work numerically. Now we can relate the parameters given by dimensionless unit-cell problems to the parameters used in main paper (denoted using tilde) as

[TABLE]

The pore pressure contribution tensor $\alpha_{ij}$ is dimensionless in both settings. Note that one can formulate the unit-cell problems also in dimensional setting, as written in the main paper. The microscale problems in the associated open-source software repository[46] are dimensionless, and one should use the relationships above in order to obtain coefficients for the dimensional equations.

Appendix C Obtaining all model parameters using elasticity tensors

In the current section, we demonstrate how to arrive to all necessary elastic properties using only elasticity of solid skeleton and effective elasticity of solid skeleton, similarly as employed by Gopinath & Mahadevan[29]. First, we define an inverse of fourth-rank tensor $\boldsymbol{\mathrm{{A}}}$ in index notation as

[TABLE]

In order to obtain pore pressure contribution tensors using elasticity matrices, we start by multiplying effective elasticity tensor (B 99) by the inverse of the skeleton elasticity and obtain

[TABLE]

which we then multiply by the identity matrix from left to get

[TABLE]

Now, in order to match the expression for the tensor $\alpha_{ij}$ (B 103), we subtract the obtained result from the identity matrix

[TABLE]

which due to symmetry in elasticity tensors $C_{ijkl}=C_{klij}$ can be rewritten as

[TABLE]

which is the final expression for the pore-pressure contribution tensor. In tensorial notation, this expression can be written using double contraction operator as

[TABLE]

This expression is exactly the same one as reported by Gopinath & Mahadevan[29, eq. 2.4]. Carrying out similar derivation for the coefficient $\mathcal{E}$, one can obtain

[TABLE]

Therefore if one is successful in finding the effective elasticity tensor by other means, the pore-pressure contribution tensor $\alpha$ and coefficient $\mathcal{E}$ can be recovered using expressions (C 110–C 111).