# Modelling wave dispersion in fluid saturating periodic scaffolds

Eduard Rohan<EMAIL_ADDRESS>Robert Cimrman<EMAIL_ADDRESS>European
Centre of Excellence, NTIS – New Technologies for Information Society, Faculty
of Applied Sciences, University of West Bohemia, Univerzitní 8, 30100 Pilsen,
Czech Republic

###### Abstract

Acoustic waves in a slightly compressible fluid saturating porous periodic
structure are studied using two complementary approaches: 1) the periodic
homogenization (PH) method provides effective model equations for a general
dynamic problem imposed in a bounded medium, 2) harmonic acoustic waves are
studied in an infinite medium using the Floquet-Bloch (FB) wave decomposition.
In contrast with usual simplifications, the advection phenomenon of the
Navier-Stokes equations is accounted for. For this, an acoustic approximation
is applied to linearize the advection term. The homogenization results are
based the periodic unfolding method combined with the asymptotic expansion
technique providing a straight upscaling procedure which leads to the
macroscopic model defined in terms of the effective model parameters. These
are computed using the characteristic responses of the porous microstructure.
Using the FB theory, we derive dispersion equations for the scaffolds
saturated by the inviscid, or the viscous, barotropic fluids, whereby the
advection due to a permanent flow in the porous structures is respected. A
computational study is performed for the numerical models obtained using the
finite element discretization. For the FB methods-based dispersion analysis,
quadratic eigenvalue problems must be solved. The numerical examples show
influences of the microstructure size and of the advection generating an
anisotropy of the acoustic waves dispersion.

###### keywords:

homogenization , Navier-Stokes equations, porous media , acoustic waves ,
Floquet-Bloch wave decomposition , wave dispersion ,

##  1 Introduction

Modelling of acoustic waves in fluid saturated porous media has been treated
mostly using the homogenization theory [27, 12], or using the phenomenological
models based on the theory of porous media [3, 4]. Although a large body of
literature devoted to this topic exists, some acoustic phenomena related to
the nonlinearity arising from the flow model deserve a further attention. New
challenges for the modelling of these phenomena emerge due to applications in
biomedicine [26], chemistry and smart material design [31]. As a step forward,
we aim to explore models of the acoustic wave propagation in viscous and
inviscid fluids saturating periodic rigid porous structures (scaffolds). Thus,
we consider waves propagating only in the fluid flowing through channels with
a periodic structure. To give an example of such structures, these can be
created by sintered ceramic fibres, see [16], where elastic waves in the solid
were studied without effects of the fluid. On the contrary, here we consider
waves propagating in the fluid only, while neglecting compliance of the solid.

It is well known that the standard Darcy law describing slow viscous flows in
porous materials can be derived rigorously by the asymptotic homogenization of
the steady Stokes flow [27], cf. [32] where the unfolding method of
homogenization was employed. Homogenization of the non-stationary
incompressible Stokes flow was treated using the two-scale convergence method
in [1]. The obtained Darcy permeability serves as the time-convolution kernel
which, upon being transformed in the frequency domain, yields the so-called
dynamic permeability. In this way, the viscodynamic operator [23] can be
introduced which comprises the tortuosity effects. Flows in periodic channels
governed by the non-stationary Navier-Stokes (N-S) equations including the
advection term were treated in [20], where the critical case of scaling the
viscosity and the velocity with respect to the pressure was discovered. Using
this critical scaling, the formal asymptotic expansion technique of
homogenization was employed in [33] to derive the Forchheimer law from the N-S
equations with the inertia terms. Upscaling of compressible, or incompressible
flows governed by the N-S equations was reported also e.g. in [21, 19, 6, 15],
computational aspects were studied in a number of papers [17, 24, 22].

In this paper, we consider acoustic wave propagation in periodic scaffolds
saturated by Newtonian, or inviscid, slightly compressible barotropic fluids.
To respect the advection effects related to the permanent flow which is
assumed to be independent of the acoustic perturbations, a linearization is
employed to establish approximate models of acoustic waves. For the inviscid
fluid, a reduced model involving the pressure only is derived. We restrict the
study to the first order linearization of the fluid advection, but assume that
the permanent flow can be quite important. In this context, to capture the
acoustic streaming effect, where the flow is generated by the acoustic waves,
higher order approximations are needed [5, 26, 14].

Two approaches are examined to analyze the wave dispersion: the periodic
homogenization (PH) and the Floquet-Bloch wave decomposition (FB), cf.  [11],
enabling to analyze waves of lengths comparable with the periodicity size.
Using both these approaches we derive models of viscous and inviscid fluids in
the rigid scaffolds, respecting the advection effect of the permanent fluid
flow upon which perturbations induced by the wave propagation are
superimposed.

Pursuing the first approach for inviscid fluid, a macroscopic model governing
the pressure fluctuations is derived, which involves an advection term related
to the time rate of the pressure. For the viscous fluid, the homogenization of
the N-S equations in the rigid skeleton provides the dynamic permeability of
the effective porous medium. The dispersion phenomenon is remarkable
especially for low frequencies, whereas constant phase velocity characterizes
the asymptotic behaviour for larger wave numbers. In both the cases, the PH
approach captures the wave propagation for wave lengths significantly larger
than the characteristic porosity size corresponding to one period of the
lattice. Using the second approach, the FB decomposition enables to capture
the wave response for wave numbers within the whole first Brillouin zone. We
derive the equations of the cell problems describing the local fluctuations of
the wave polarization. The dispersion analysis leads to quadratic eigenvalue
problems (QEP). In the case of viscous fluids, the QEP with rather complicated
structure appears, that requires a suitable linearization demanding new
variables to be introduced. The computational analysis for the inviscid fluid
model involving 3D scaffolds with variable porosity is performed, showing the
porosity influence on a frequency band gap opening between the two lowest
modes. However, the main emphasis is put on the illustration of the permanent
flow advection influence on the sound speed, in both the inviscid and viscous
fluids.

The plan of the paper is as follows. In Section 2 the model of fluid flow in
porous structures and the acoustic waves are introduced using a decomposition
of the model responses into the steady part and the acoustic fluctuations. For
the inviscid and viscous fluid saturating periodic rigid scaffolds, the
respective homogenized models are derived in Section 3, where also the
pressure plane wave propagation in the homogenized medium is described for the
two fluids. In Section 4, the Floquet-Bloch wave decomposition is applied to
analyze the wave propagation in the considered porous structures; the
generalized quadratic eigenvalue problems resulting from the finite-element
discretization of the weak formulations are introduced and their linearization
is explained. Both the modelling approaches are illustrated on examples of 2D
and 3D periodic scaffolds in Section 5. There the dispersion analysis
performed by solving the eigenvalue problems resulting from the Floquet-Bloch
wave decomposition is compared with the corresponding approximation obtained
by the homogenized model for viscous and inviscid fluids, whereby the size
effects and advection phenomena are explored. Some technical results employed
in the paper are postponed in the Appendix.

##### General notations

Spatial position of a point is specified by its Cartesian coordinates,
$x=(x_{1},x_{2},x_{3})\in{\mathbb{R}}^{3}$, where ${\mathbb{R}}$ is the set of
real numbers. The boldface notation for vectors ${\textbf{{a}}}=(a_{i})$ and
second-order tensors ${\textbf{{b}}}=(b_{ij})$ is used. The second-order
identity tensor is denoted by ${\textbf{{I}}}=(\delta_{ij})$. The fourth-order
elasticity tensor is denoted by ${{\rm I}\kern-1.99997pt{\rm D}}=(D_{ijkl})$.
The superposed dot denotes a derivative with respect to time. The gradient,
divergence and Laplace operators are denoted by $\nabla,\nabla\cdot$ and
$\nabla^{2}$, respectively. When these operators have a subscript referring to
the space variable, it is for indicating that the operator acts relatively at
this space variable, for instance
$\nabla_{y}=(\partial_{i}^{y})=(\partial/\partial y_{i})$. The symbol dot
‘$\cdot$’ denotes the scalar product between two vectors and the symbol colon
‘$:$’ stands for scalar (inner) product of two second-order tensors, e.g.
${\textbf{{A}}}:{\textbf{{B}}}=A_{ij}B_{ij}={\tt{tr}}[{\textbf{{A}}}^{T}{\textbf{{B}}}]=A_{ki}B_{kj}\delta_{ij}$,
where ${\tt{tr}}[\star]$ is the trace of a tensor and superscript $T$ in
$\star^{T}$ is the transposition operator. Operator $\otimes$ designates the
tensor product between two vectors, e.g.
${\textbf{{a}}}\otimes{\textbf{{v}}}=(a_{i}v_{j})$. Standard notations for
functional spaces are adhered. Throughout the paper, $x$ denotes the global
(“macroscopic”) coordinates, while the “local” coordinates $y$ describe
positions within the representative unit cell $Y\subset{\mathbb{R}}^{3}$. The
normal vectors on a boundary of domains $\Omega_{\alpha}$ (or $Y_{\alpha}$)
are denoted by ${\textbf{{n}}}^{\alpha}$, $\alpha=s,f$, to distinguish their
orientation outward to $\Omega_{\alpha}$ (or $Y_{\alpha}$) when dealing with
the solid-fluid interfaces. By
${\textbf{{e}}}({{\textbf{{w}}}})=1/2(\nabla{\textbf{{w}}}+(\nabla{\textbf{{w}}})^{T})$
we denote the strain of a vector field w (displacements, or velocities). The
“tilde”-notation can have various meanings which are explained through the
text and are clear within the particular context. The following standard
functional spaces are used: by $L^{2}(\Omega)$ we refer to square integrable
functions defined in an open bounded domain $\Omega$; by $H^{1}(\Omega)$ we
mean the Sobolev space $W^{1,2}(\Omega)\subset L^{2}(\Omega)$ formed by square
integrable functions including their first generalized derivatives; space
$C_{0}^{\infty}(\Omega)$ is constituted by infinitely differentiable functions
with the compact support, thus, with zero trace on $\partial\Omega$; Bold
notation is used to denote spaces of vector-valued functions, e.g.
${\bf{H}}^{1}(\Omega)$; by subscript # we refer to the _$Y$ -periodic
functions_.

##  2 Model of fluid flow in porous structures and acoustic waves

We introduce formulations for acoustic waves in a fluid saturating channels
$\Omega_{f}$ of a porous two-phase medium situated in an open bounded domain
$\Omega\subset{\mathbb{R}}^{3}$. The solid skeleton occupies domain
$\Omega_{s}\subset\Omega$, whereas viscous, or inviscid slightly compressible
fluids saturate the pores $\Omega_{f}=\Omega\setminus\overline{\Omega_{s}}$
constituting a connected network of channels whose walls
$\Gamma_{\it{fs}}=\overline{\Omega_{s}}\cap\overline{\Omega_{f}}$ are
impermeable for the fluid. Each of the two domains $\Omega_{k}$, with $k=s,f$
are connected, hence also their interface $\Gamma_{\it{fs}}$ is connected. The
fluid can flow through, or be loaded on the external part of
$\partial\Omega_{f}$ denoted by
$\partial_{\rm{ext}}\Omega_{f}=\partial\Omega_{f}\setminus\Gamma_{fs}$.
Obviously, $\partial_{\rm{ext}}\Omega_{f}\subset\partial\Omega$ can be
decomposed according to the specific boundary conditions.

We recall the Navier-Stokes equations; the fluid velocity w and pressure $p$
satisfy

$\begin{split}\rho\left(\frac{{\rm{d}}\,}{{\rm{d}}\,t}{\textbf{{w}}}+{\textbf{{w}}}\cdot\nabla{\textbf{{w}}}\right)&=-\nabla
p+\nabla\cdot{{\rm I}\kern-1.99997pt{\rm
D}}{\textbf{{e}}}({{\textbf{{w}}}})\;,\\\
\frac{\partial{\rho}}{\partial{t}}+\nabla\cdot(\rho{\textbf{{w}}})=0\;,\end{split}$
(2.1)

where $\rho$ is the fluid density, and ${{\rm I}\kern-1.99997pt{\rm
D}}{\textbf{{e}}}({{\textbf{{w}}}})$ represents the viscous stress given by
the velocity strain
${\textbf{{e}}}({{\textbf{{w}}}})=1/2(\nabla{\textbf{{w}}}+(\nabla{\textbf{{w}}})^{T})$
and by the viscosity tensor, ${{\rm I}\kern-1.99997pt{\rm D}}=(D_{ijkl})$ with
$D_{ijkl}=\eta\delta_{ij}\delta_{kl}+\mu(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})$
depending on the 1st and the 2nd viscosity, $\mu$ and $\eta$, respectively.

###  2.1 Assumptions for acoustic wave decomposition and linearization

We shall study acoustic waves propagating in a slightly compressible viscous,
or inviscid fluid, while thermal effects are disregarded (the barotropic
fluid). The fluid flows through the rigid porous structure, whereby the
following assumptions are made:

  * (A1) The total fields w, $p$ and the mass density of the fluid, $\rho$, are split into the “stationary flow” parts $\bar{\textbf{{w}}}$, $\bar{\rho}$ and $\bar{p}$ and the “acoustic fluctuation” parts $\tilde{\textbf{{w}}}$, $\tilde{\rho}$ and $\tilde{p}$, so that

$\begin{split}{\textbf{{w}}}=\bar{\textbf{{w}}}+\tilde{\textbf{{w}}}\;,\quad
p=\bar{p}+\tilde{p}\;,\quad\rho=\bar{\rho}+\tilde{\rho}\;.\end{split}$ (2.2)

  * (A2) The fluid is assumed to be homogeneous and under the stationary flow described by $(\bar{\textbf{{w}}},\bar{p},\bar{\rho})$ is considered as incompressible, thus

$\begin{split}\nabla\bar{\rho}&=0\;,\quad\frac{\partial{}}{\partial{t}}\bar{\textbf{{w}}}=0\;,\quad\frac{\partial{}}{\partial{t}}\bar{\rho}=0\;.\end{split}$
(2.3)

  * (A3) The acoustic perturbations introduced due to the split (2.2) are small, so that mutual multiplications of any two perturbations is neglected. This assumption is important for the linearization.

  * (A4) The acoustic response is barotropic: denoting by $p_{0}$ and $\rho_{0}$ reference state variables, it holds that $p-p_{0}=c^{2}(\rho-\rho_{0})$, where $c^{2}$ is the squared acoustic velocity $c=\sqrt{k_{f}/\rho_{0}}$ with the bulk stiffness $k_{f}=1/\gamma_{f}$, thus, $\gamma_{f}$ is the fluid compressibility for the reference state. Since we consider the steady flow as incompressible, $\bar{\rho}=\bar{\rho}_{0}$, the reference pressure $p_{0}$ can be associated with the one of the steady flow, i.e. $\bar{p}\approx\bar{p}_{0}$. As a consequence, we have $\tilde{p}=c^{2}\tilde{\rho}$. In this paper we adhere to this linear approximation which is adequate in the context of the incompressibility of the stationary flow part, assumption A2). A nonlinear relationship between the pressure and the density, as treated e.g. by [13], is required to handle some “higher order effects”, such as the acoustic streaming, see e.g. [14].

  * (A5) The stationary flow is characterized by a periodic velocity field $\bar{\textbf{{w}}}$. This assumption is needed to analyze the acoustic waves using the Floquet-Bloch theory.

The above assumptions lead the following equations governing the acoustic
fluctuations $(\tilde{\textbf{{u}}},\tilde{p})$ of the velocity and pressure
fields (see A),

$\begin{split}\rho_{0}\left(\frac{{\rm{d}}\,}{{\rm{d}}\,t}\tilde{\textbf{{u}}}+\bar{\textbf{{w}}}\cdot\nabla\tilde{\textbf{{u}}}+\tilde{\textbf{{u}}}\cdot\nabla\bar{\textbf{{w}}}\right)&=-\nabla\tilde{p}+\nabla\cdot{{\rm
I}\kern-1.99997pt{\rm D}}{\textbf{{e}}}({\tilde{\textbf{{u}}}})\;,\\\
\frac{{\rm{d}}\,}{{\rm{d}}\,t}\tilde{p}+{\textbf{{w}}}\cdot\nabla\tilde{p}&=-k_{f}\nabla\cdot\tilde{\textbf{{u}}}\;,\end{split}$
(2.4)

while $(\bar{\textbf{{w}}},\bar{p})$ describing the stationary flow in the
periodic scaffolds satisfy

$\begin{split}\rho_{0}\bar{\textbf{{w}}}\cdot\nabla\bar{\textbf{{w}}}+\nabla\bar{p}-\mu\nabla^{2}{\bar{\textbf{{w}}}}&=\bar{\textbf{{f}}}\;,\\\
\nabla\cdot\bar{\textbf{{w}}}&=0\;,\end{split}$ (2.5)

where f is the volume force.

Concerning Assumption (A5), in general, periodic heterogeneities in the flow
can be caused by different aspects:

  * 1.

presence of periodically distributed obstacles (or flow in a periodic porous
structure) – static fluid;

  * 2.

periodic heterogeneities in the fluid parameters; such a situation is relevant
for a static fluid, i.e. $\bar{\textbf{{w}}}\equiv\bf{0}$, for instance two
fluids, air bubbles…);

  * 3.

fluid flow through a periodic porous structure with non-negligible convection
velocity which is periodically perturbed due to the porous structure.

###  2.2 Inviscid fluid

In this section we consider an inviscid compressible fluid characterized by
constants $\rho_{0}$ and $k_{f}$ only, since the vanishing viscosity cancels
the second right hand side term in (2.4)1 involving ${{\rm
I}\kern-1.99997pt{\rm D}}$.

While for the viscous fluid model the advection velocity $\bar{\textbf{{w}}}$
is computed using (2.5), for inviscid fluids the velocity field
$\bar{\textbf{{w}}}$ is assumed to have vanishing vorticity. Hence, it is
determined by (2.5)2 governing the steady incompressible flow:
$\bar{\textbf{{w}}}=-\nabla\Psi$, where the potential $\Psi$ satisfies

$.\begin{split}-\nabla^{2}\Psi&=0\quad\mbox{ in }\Omega_{f}\;,\\\
\Psi&=\bar{\Psi}\quad\mbox{ on }\partial_{\rm{ext}}\Omega_{f}\;,\\\
{\mbox{\boldmath$\nu$\unboldmath}}\cdot\nabla\Psi&=0\quad\mbox{ on
}\Gamma_{\it{fs}}\;.\end{split}$ (2.6)

The incompressibility condition $\nabla\cdot\bar{\textbf{{w}}}=0$ in
$\Omega_{f}$ and the impermeability of $\Gamma_{\it{fs}}$ are verified.

For stationary fluids, i.e. when $\bar{\textbf{{w}}}\equiv 0$, it is an easy
exercise to eliminate velocity $\tilde{\textbf{{u}}}$ from (2.4), recalling
${{\rm I}\kern-1.99997pt{\rm D}}$ is zero, so that the Helmholtz equation
governing the pressure fluctuations $\tilde{p}$ is obtained. Alternatively,
upon eliminating the pressure, the velocity must satisfy the wave equation
supplemented by the zero vorticity constraint. The acoustic pressure waves in
the free fluid propagate with the phase velocity
$c_{f}=\sqrt{k_{f}/\rho_{0}}$, whereas the acoustic waves in the scaffolds are
governed by either of the following two hyperbolic equations,

$\begin{split}\mbox{ pressure formulation:
}\quad\frac{{\rm{d^{2}}}\,}{{\rm{d}}\,{t}^{2}}\tilde{p}&=c_{f}^{2}\nabla^{2}\tilde{p}\;,\\\
\mbox{ velocity formulation:
}\quad\frac{{\rm{d^{2}}}\,}{{\rm{d}}\,{t}^{2}}\tilde{\textbf{{u}}}&=c_{f}^{2}\nabla(\nabla\cdot\tilde{\textbf{{u}}})\quad\mbox{
and }\nabla\times\tilde{\textbf{{u}}}=\bf{0}\;.\end{split}$ (2.7)

When steady advection is respected, $\bar{\textbf{{w}}}\not\equiv\bf{0}{}$,
similar reduced formulations can be derived. However, it is not
straightforward to eliminate the velocity and to arrive at a pressure
formulation. For this, divergence operator is applied to (2.4)2, which yields,

$\begin{split}\rho_{0}\left(\nabla\cdot\dot{\tilde{\textbf{{u}}}}+2\partial_{k}\bar{w}_{i}\partial_{i}\tilde{u}_{k}+\bar{\textbf{{w}}}\cdot\nabla(\nabla\cdot\tilde{\textbf{{u}}})+\tilde{\textbf{{u}}}\nabla(\nabla\cdot\bar{\textbf{{w}}})\right)=-\nabla^{2}\tilde{p}\;,\end{split}$
(2.8)

where the last term on the left-hand side vanishes due to the
incompressibility of $\bar{\textbf{{w}}}$.

The following lemma is employed to treat the second left-hand side term.

######  Lemma 2.1.

For any open bounded domain $Q\subset\Omega_{f}$, it holds that

$\begin{split}\int_{Q}\partial_{k}\bar{w}_{i}\partial_{i}\tilde{u}_{k}=\int_{\partial
Q}\left(\bar{w}_{i}\partial_{i}\tilde{u}_{k}-\bar{w}_{k}\partial_{i}\tilde{u}_{i}\right)n_{k}+\int_{Q}\nabla\cdot\bar{\textbf{{w}}}\nabla\cdot\tilde{\textbf{{u}}}\;.\end{split}$
(2.9)

The proof is a straightforward consequence of
$\partial_{i}(\partial_{k}\tilde{u}_{i})=\partial_{k}(\nabla\cdot\tilde{\textbf{{u}}})$.

As the consequence of Lemma 2.1 and the incompressibility
$\nabla\cdot\bar{\textbf{{w}}}=0$, the following approximation may be
considered,

$\begin{split}\left\langle{\partial_{i}\bar{w}_{k}\partial_{k}\tilde{u}_{i}}\right\rangle_{Q}&\approx|Q|^{-1}\int_{\partial
Q}\left(\bar{w}_{k}\partial_{k}\tilde{u}_{i}-\bar{w}_{i}\left\langle{\nabla\cdot\tilde{\textbf{{u}}}}\right\rangle_{Q}\right)n_{i}\\\
&\approx\bar{w}_{k}\partial_{k}\left\langle{\nabla\cdot\tilde{\textbf{{u}}}}\right\rangle_{Q}-|Q|^{-1}\int_{\partial
Q}\bar{w}_{i}n_{i}\left\langle{\nabla\cdot\tilde{\textbf{{u}}}}\right\rangle_{Q}\\\
&=\partial_{w}\left\langle{\nabla\cdot\tilde{\textbf{{u}}}}\right\rangle_{Q}\;,\end{split}$
(2.10)

where $\left\langle{\leavevmode\nobreak\ }\right\rangle_{Q}=|Q|^{-1}\int_{Q}$
is the average. Therefore, we can employ the approximate relationship

$\begin{split}\partial_{i}\bar{w}_{k}\partial_{k}\tilde{u}_{i}\approx\bar{\textbf{{w}}}\cdot\nabla(\nabla\cdot\tilde{\textbf{{u}}})\;.\end{split}$
(2.11)

Using (2.11) and upon substituting $\nabla\cdot\tilde{\textbf{{u}}}$ by
(2.4)2, the left hand side in (2.8) can be approximated, so that the acoustic
waves in the pore fluid are governed by the following equation:

$\begin{split}\left(\frac{{\rm{d}}\,}{{\rm{d}}\,t}+\zeta\bar{\textbf{{w}}}\cdot\nabla\right)\left(\dot{\tilde{p}}+\bar{\textbf{{w}}}\cdot\nabla{\tilde{p}}\right)&=c_{f}^{2}\nabla^{2}\tilde{p}\;,\end{split}$
(2.12)

where $\zeta=3$ and $c_{f}=\sqrt{k_{f}/\rho_{0}}$. Further we define
$\theta=(1+\zeta)/2$, thus, $\theta=2$. However, in what follows, we keep the
abstract notation $\zeta$ and $\theta$.

##  3 Homogenization of fluid in periodic porous structures

We are confined to rigid scaffolds saturated by viscous, or inviscid fluids,
so that, at the pore level, the acoustic waves are described by the models
introduced in Section 2. The size of the porous microstructure is expressed by
the scale parameter $\varepsilon=\ell/L$ defined by the ratio of the micro-
and macroscopic characteristic lengths, denoted by $\ell$ and $L$,
respectively. The homogenization procedure applied to derive an effective
model of the rigid-porous medium consists in the asymptotic analysis
$\varepsilon\rightarrow 0$ of the micro-model presented above.

The notation related to the geometry formerly introduced is adhered to,
however, now the superscript ε is appended to respect the dependence on the
scale, i.e. the fluid occupies pores
$\Omega_{f}^{\varepsilon}=\Omega\setminus\overline{\Omega_{s}^{\varepsilon}}$,
bounded by the pore walls
$\Gamma_{\it{fs}}^{\varepsilon}=\overline{\Omega_{s}^{\varepsilon}}\cap\overline{\Omega_{f}^{\varepsilon}}$,
and by the external boundary of the pores
$\partial_{\rm{ext}}\Omega_{f}^{\varepsilon}$.

Figure 1: A 2D scheme of the heterogeneous porous structure. (In 3D, both the
solid and fluid subdomains $\Omega_{s}^{\varepsilon}$ and
$\Omega_{f}^{\varepsilon}$ are connected.) Left: decomposition of the domain
$\Omega$ and its boundary $\partial\Omega$ for a given scale $\varepsilon$
into the solid and fluid parts. Right: the representative unit periodic cell
$Y=\prod_{i=1}^{3}]0,\bar{y}^{i}[$.

###  3.1 Periodic microstructure

We consider a periodic structure of channels $\Omega_{f}^{\varepsilon}$
saturated by one homogeneous fluid. Due to the assumed periodicity of
$\Omega_{f}^{\varepsilon}$, a representative periodic cell
$\mathcal{Z}_{f}^{\varepsilon}$ can be defined, which generates the fluid
domain, see Fig. 1. A periodic cell
$\mathcal{Z}^{\varepsilon}=\prod_{i=1}^{3}]0,\varepsilon\bar{y}^{i}[$ can be
introduced, such that
$\mathcal{Z}^{\varepsilon}=\mathcal{Z}_{f}^{\varepsilon}\cup\mathcal{Z}_{s}^{\varepsilon}\cup\partial_{s}\mathcal{Z}_{f}^{\varepsilon}$,
where $\mathcal{Z}_{s}^{\varepsilon}\subset\mathcal{Z}^{\varepsilon}$ is the
solid skeleton and $\partial_{s}\mathcal{Z}_{f}^{\varepsilon}$ is the fluid-
solid interface. By
$\partial_{\\#}\mathcal{Z}_{f}^{\varepsilon}=\partial\mathcal{Z}_{f}^{\varepsilon}\setminus\partial_{s}\mathcal{Z}_{f}^{\varepsilon}$
we denote the “periodic part” of the boundary. For the purpose of the
homogenization we consider the “unit periodic cell”
$Y=\varepsilon^{-1}\mathcal{Z}^{\varepsilon}$ which consists of the fluid and
solid parts, $Y_{f}=\varepsilon^{-1}\mathcal{Z}_{f}^{\varepsilon}$ and
$Y_{s}=\varepsilon^{-1}\mathcal{Z}_{s}^{\varepsilon}$ , respectively,
accordingly the decomposition of $\mathcal{Z}^{\varepsilon}$, thus,
$Y=Y_{f}\cup Y_{s}\cup\Gamma$, where
$\Gamma=\overline{Y_{f}}\cap\overline{Y_{s}}$ is the interface.

The homogenization procedure is presented below formally, without giving
convergence proofs, however, the derivations of the limit equations can be
followed. For this purpose, basics of the unfolding method are summarized in
B. We use the standard notation; for any $D\subset Y$ we abbreviate
$\sim\kern-10.70009pt\int_{D}=\frac{1}{|Y|}\int_{D}$; note that usually one
may chose $\bar{y}^{i}$, $i=1,2,3$, such that $|Y|=1$. Moreover, the following
spaces are employed: $H_{\\#}^{1}(Y_{f})\subset H^{1}(Y_{f})$ containing only
$Y$-periodic functions. In analogy, a space of $Y$-periodic vector-valued
functions is denoted by ${\bf{H}}_{\\#}^{1}(Y_{f})$.

###  3.2 Homogenization of an inviscid fluid

We recall the notation with the superscript ε indicating dependence on the
scale of pores $\Omega_{f}^{\varepsilon}$. The parameters $\theta=2$,
$\zeta=3$ and $c_{f}$ are fixed. The homogenization procedure based on the
periodic unfolding method [8] will be applied formally using the truncated
asymptotic expansions of the pressure $p^{\varepsilon}$ and the associated
test functions $q^{\varepsilon}$ involved in the weak formulation introduced
below. It arises from (2.12) supplemented by the Dirichlet boundary conditions
given on $\partial_{\rm{ext}}\Omega_{f}^{\varepsilon}$ by
$\tilde{p}^{\partial}$, and non-penetration conditions on
$\Gamma_{fs}^{\varepsilon}$, thus
$\tilde{\textbf{{u}}}\cdot{\mbox{\boldmath$\nu$\unboldmath}}=0$. The specific
choice of the boundary conditions on
$\partial_{\rm{ext}}\Omega_{f}^{\varepsilon}$ has no limiting consequences,
since our aim is to analyze the wave dispersion in an unbounded media,

####  3.2.1 Weak formulation for the acoustic waves with advection in an
inviscid fluid

Recalling Assumption (A5), an advection velocity field
$\bar{\textbf{{w}}}^{\varepsilon}$ involved in (2.12) by virtue of the
projected gradients
$\partial_{w}^{\varepsilon}:=\bar{\textbf{{w}}}^{\varepsilon}\cdot\nabla$ must
be introduced. For this we consider homogenization of the potential flow
(2.6); it is a classical academic problem of the homogenization in the
perforated domains, see e.g. [2, 9], nevertheless, a brief information on
computing $\bar{\textbf{{w}}}^{\varepsilon}$ in response to a given
macroscopic velocity $\bar{\textbf{{w}}}^{0}$ is given in C.

We shall consider the following problem: Find $\tilde{p}\in H^{1}(\Omega_{f})$
such that $\tilde{p}=\tilde{p}^{\partial}$ a.e. on
$\partial_{\rm{ext}}\Omega_{f}$, and

$\begin{split}\int_{\Omega_{f}^{\varepsilon}}q^{\varepsilon}\left(\ddot{\tilde{p}}^{\varepsilon}+\theta\partial_{w}^{\varepsilon}\dot{\tilde{p}}^{\varepsilon}\right)-\int_{\Omega_{f}^{\varepsilon}}\partial_{w}^{\varepsilon}q^{\varepsilon}\left(\theta\dot{\tilde{p}}^{\varepsilon}+\zeta\partial_{w}^{\varepsilon}{\tilde{p}^{\varepsilon}}\right)+c_{f}^{2}\int_{\Omega_{f}^{\varepsilon}}\nabla{\tilde{p}^{\varepsilon}}\cdot\nabla
q^{\varepsilon}=0\;,\end{split}$ (3.1)

holds for a.a. $q^{\varepsilon}\in H^{1}(\Omega_{f}^{\varepsilon})$ such that
$q^{\varepsilon}=0$ a.e. on $\partial_{\rm{ext}}\Omega_{f}^{\varepsilon}$.

Note that (3.1) results from (2.12) which is multiplied by $q$ and integrated
in $\Omega_{f}$, so that the obvious integration by parts leads to the
following boundary integral,

$\begin{split}\mathcal{I}_{\partial\Omega_{f}^{\varepsilon}}(\tilde{p}^{\varepsilon},q^{\varepsilon})=\int_{\partial\Omega_{f}}q^{\varepsilon}\left({\mbox{\boldmath$\nu$\unboldmath}}^{f}\cdot\nabla{\tilde{p}^{\varepsilon}}-\zeta
w_{\nu}^{\varepsilon}\partial_{w}^{\varepsilon}{\tilde{p}}-\theta
w_{\nu}^{\varepsilon}\dot{\tilde{p}}^{\varepsilon}\right)\;,\end{split}$ (3.2)

where
$w_{\nu}^{\varepsilon}=\bar{\textbf{{w}}}^{\varepsilon}\cdot{\mbox{\boldmath$\nu$\unboldmath}}^{f}$.
This integral vanishes on $\partial_{\rm{ext}}\Omega_{f}$ due to the
considered Dirichlet boundary conditions. On the walls $\Gamma_{\it{fs}}$, all
the integrands in (3.2) vanish as well due to the non-penetration condition
which yields
${\mbox{\boldmath$\nu$\unboldmath}}^{f}\cdot\nabla\tilde{p}^{\varepsilon}=0$
by virtue of (2.4)2, since, in general, we assume
${\mbox{\boldmath$\nu$\unboldmath}}^{f}\cdot(\partial_{w}\tilde{\textbf{{u}}}+\tilde{\textbf{{u}}}\cdot\nabla\bar{\textbf{{w}}})=0$.
For polygonal (polyhedral) interfaces this assumption can be verified, since
$\partial_{w}{\mbox{\boldmath$\nu$\unboldmath}}^{f}:=\bar{\textbf{{w}}}\cdot\nabla{\mbox{\boldmath$\nu$\unboldmath}}^{f}=0$.

####  3.2.2 Limit problem

We assume the advection velocity $\bar{\textbf{{w}}}^{\varepsilon}$ to be of
the order $o(\varepsilon)=1$, thus, using the unfolding operator, see B,
${\mathcal{T}}_{\varepsilon}{\left({\bar{\textbf{{w}}}^{\varepsilon}}\right)}=\bar{\textbf{{w}}}(x,y)$.
Also the fluid phase velocity $c_{f}$ is assumed to be independent of
$\varepsilon$, so that the formal asymptotic expansion of
$\tilde{p}^{\varepsilon}$ is

$\begin{split}{\mathcal{T}}_{\varepsilon}{\left({\tilde{p}^{\varepsilon}}\right)}&=p^{0}(x,t)+\varepsilon
p^{1}(x,y,t)\;.\end{split}$ (3.3)

The following sets are employed:

$\begin{split}W_{0}(\Omega)&=H_{0}^{1}(\Omega)=\\{q\in
H^{1}(\Omega)\,|\;q=0\mbox{ on }\partial\Omega\\}\;,\\\ W_{*}(\Omega)&=\\{q\in
H^{1}(\Omega)\,|\;q=p^{\partial}\mbox{ on }\partial\Omega\\}\;,\end{split}$
(3.4)

We may assume the Dirichlet boundary condition represented by
$p^{\partial,\varepsilon}$ converging with $\varepsilon\rightarrow 0$ to
$p^{\partial}$ defined on $\partial\Omega$. Hence, by virtue of the
convergence
${\mathcal{T}}_{\varepsilon}{\left({\tilde{p}^{\varepsilon}}\right)}\rightarrow
p^{0}$ we get $p^{0}\in W_{*}(\Omega)$.

######  Remark  1.

By virtue of the weak convergence of the unfolded pressure gradient, i.e.
${\mathcal{T}}_{\varepsilon}{\left({p^{\varepsilon}}\right)}(\cdot,t)\rightharpoonup\nabla_{x}p^{0}(\cdot,t)+\nabla_{y}p^{1}(\cdot,t)$
weakly in $L^{2}(\Omega\times Y_{y})$, the trace of $p^{0}(\cdot,t)\in
H^{1}(\Omega)$ on $\partial_{p}\Omega$ provides the limit trace
$p^{0}=p^{\partial}$ on $\Omega$.

We define $\partial_{w}^{x}=\bar{\textbf{{w}}}\cdot\nabla_{x}$ and
$\partial_{w}^{y}=\bar{\textbf{{w}}}\cdot\nabla_{y}$. In the limit
$\varepsilon\rightarrow 0$, (3.1) yields the two-scale problem for $p^{0}$ and
$p^{1}$ satisfying (recall B for the integration formula and
$\sim\kern-10.70009pt\int_{D}=\frac{1}{|Y|}\int_{D}$ for any $D\subset Y$)

$\begin{split}\int_{\Omega}q^{0}\ddot{p}^{0}+\theta\int_{\Omega}q^{0}\sim\kern-11.99998pt\int_{Y_{f}}\left(\partial_{w}^{x}\dot{p}^{0}+\partial_{w}^{y}\dot{p}^{1}\right)-\theta\int_{\Omega}\dot{p}^{0}\sim\kern-11.99998pt\int_{Y_{f}}\left(\partial_{w}^{x}q^{0}+\partial_{w}^{y}q^{1}\right)&\\\
-\zeta\int_{\Omega}\sim\kern-11.99998pt\int_{Y_{f}}\left(\partial_{w}^{x}p^{0}+\partial_{w}^{y}p^{1}\right)\left(\partial_{w}^{x}q^{0}+\partial_{w}^{y}q^{1}\right)&\\\
+c_{f}^{2}\int_{\Omega}\sim\kern-11.99998pt\int_{Y_{f}}(\nabla_{x}p^{0}+\nabla_{y}p^{1})\cdot(\nabla_{x}q^{0}+\nabla_{y}q^{1})&=0\;,\end{split}$
(3.5)

for all $q^{0}\in W_{0}(\Omega)$ and $q^{1}\in
L^{2}(\Omega;H_{\\#}^{1}(Y_{f}))$. Due to the incompressibility
$\nabla\cdot\bar{\textbf{{w}}}=0$ and since $\bar{\textbf{{w}}}(x,\cdot)$ is
$Y$-periodic, for any $Y$-periodic function $\psi\in H_{\\#}^{1}(Y_{f})$

$\begin{split}\sim\kern-11.99998pt\int_{Y_{f}}\partial_{w}\psi=\sim\kern-11.99998pt\int_{Y_{f}}\bar{\textbf{{w}}}\cdot\nabla_{y}\psi=0\;,\end{split}$
(3.6)

so that, in the first line of (3.5), both the integrals of
$\partial_{w}^{y}\dot{p}^{1}$ and $\partial_{w}^{y}q^{1}$ vanish. The local
problem is distinguished in (3.5) when putting $q^{0}\equiv 0$. By virtue of
the linearity, the characteristic responses $\pi^{k}\in H_{\\#}^{1}(Y_{f})$,
$k=1,\dots,3$ are introduced, such that

$\begin{split}p^{1}=\pi^{k}\partial_{k}^{x}p^{0}\;.\end{split}$ (3.7)

Using the bilinear form (note that it involves $\zeta=3$),

$\begin{split}a_{w}\left({p},\,{q}\right)&=\sim\kern-11.99998pt\int_{Y_{f}}\left(c_{f}^{2}\nabla_{y}p\cdot\nabla_{y}q-\zeta\partial_{w}^{y}p\partial_{w}^{y}q\right)\\\
&=\sim\kern-11.99998pt\int_{Y_{f}}\left[\left(c_{f}^{2}{\textbf{{I}}}-\zeta\bar{\textbf{{w}}}\otimes\bar{\textbf{{w}}}\right)\nabla_{y}p\right]\cdot\nabla_{y}q\;,\end{split}$
(3.8)

the autonomous problem for computing $\pi^{k}$ is established: Find
$\pi^{k}\in H_{\\#}^{1}(Y_{f})$, $k=1,\dots,3$, such that

$\begin{split}a_{w}\left({\pi^{k}},\,{\psi}\right)&=-a_{w}\left({y_{k}},\,{\psi}\right)\;,\quad\forall\psi\in
H_{\\#}^{1}(Y_{f})\;.\end{split}$ (3.9)

The macroscopic model of the acoustic waves which is presented below involves
the homogenized coefficients
${\mbox{\boldmath$\mathcal{A}$\unboldmath}}=(\mathcal{A}_{ij})$ and the mean
advection velocity ${\textbf{{w}}}^{0}$, describing effective medium
properties,

$\begin{split}\mathcal{A}_{ij}&=a_{w}\left({\pi^{i}+y_{i}},\,{y_{j}}\right)=a_{w}\left({\pi^{i}+y_{i}},\,{\pi^{j}+y_{j}}\right)\;,\\\
{\textbf{{w}}}^{0}&=\sim\kern-11.99998pt\int_{Y_{f}}\bar{\textbf{{w}}}\;,\end{split}$
(3.10)

where the alternative symmetric expression for the _anisotropic acoustic phase
speed_ $\mathcal{A}$ can be obtained due to problem (3.9). Expression (3.10)2
for the mean velocity ${\textbf{{w}}}^{0}$ is consistent with the definition
of field $\bar{\textbf{{w}}}$ by virtue of the homogenization, see (C). In
this context, when ${\textbf{{w}}}^{0}(x)$ varies with the macroscopic
position, all expressions in (3.7)-(3.10) should be understood pointwise for
a.a. $x\in\Omega$.

####  3.2.3 Macroscopic model of acoustic waves

From the two-scale equation (3.5), the macroscopic problem is obtained for
$q^{1}\equiv 0$ and the two-scale function $p^{1}$ is substituted by the split
(3.7), which yields

$\begin{split}&\int_{\Omega}q^{0}\left(\phi_{f}\ddot{p}^{0}+\theta\sim\kern-11.99998pt\int_{Y_{f}}\partial_{w}^{x}\dot{p}^{0}\right)-\int_{\Omega}\dot{p}^{0}\theta\sim\kern-11.99998pt\int_{Y_{f}}\partial_{w}^{x}q^{0}\\\
&+\int_{\Omega}\zeta\sim\kern-11.99998pt\int_{Y_{f}}\left(\partial_{w}^{x}p^{0}+\partial_{w}^{y}\pi^{k}\partial_{k}^{x}p^{0}\right)\partial_{w}^{x}q^{0}+c_{f}^{2}\int_{\Omega}\sim\kern-11.99998pt\int_{Y_{f}}\left(\nabla_{x}p^{0}+\nabla_{y}\pi^{k}\partial_{k}^{x}p^{0}\right)\cdot\nabla_{x}q^{0}=0\;,\end{split}$
(3.11)

hence

$\begin{split}\int_{\Omega}\phi_{f}\ddot{p}^{0}q^{0}+\theta\int_{\Omega}\left({\textbf{{w}}}^{0}\cdot\nabla\dot{p}^{0}q^{0}-\dot{p}^{0}{\textbf{{w}}}^{0}\cdot\nabla
q^{0}\right)+\int_{\Omega}a_{w}\left({\pi^{k}+y_{k}},\,{y_{l}}\right)\partial_{k}^{x}p^{0}\partial_{l}^{x}q^{0}=0\;.\end{split}$
(3.12)

In the last integral, the expression involving the characteristic response
$\pi^{k}$ is substituted using the homogenized $\mathcal{A}$, whereby the
symmetry relationships follow due to the local micro-problem (3.9).

##### Macroscopic problem

Find $p^{0}(\cdot,t)\in W_{*}(\Omega)$, such that

$\begin{split}\int_{\Omega}\phi_{f}\ddot{p}^{0}q^{0}+\int_{\Omega}({\mbox{\boldmath$\mathcal{A}$\unboldmath}}\nabla_{x}p^{0})\cdot\nabla_{x}q^{0}+\theta\int_{\Omega}{\textbf{{w}}}^{0}\cdot(\nabla_{x}\dot{p}^{0}q^{0}-\nabla_{x}q^{0}\dot{p}^{0})=0\;,\end{split}$
(3.13)

for all $q^{0}\in W_{0}(\Omega)$.It is worth to recall that the limit field
$p^{0}$ satisfies the limit Dirichlet boundary condition given by
$p^{\partial}$ on $\partial\Omega$, as pointed above.

From (3.13), the differential equation can be extracted upon integration by
parts in $\Omega$,

$\begin{split}\phi_{f}\ddot{p}^{0}-\nabla\cdot({\mbox{\boldmath$\mathcal{A}$\unboldmath}}\nabla_{x}p^{0})+2\theta{\textbf{{w}}}^{0}\cdot\nabla_{x}\dot{p}^{0}=0\;,\end{split}$
(3.14)

whereby the boundary integrals on $\partial\Omega$ vanish due to $q^{0}\in
W_{0}(\Omega)$.

######  Remark  2.

Besides the gradient projection
${\textbf{{w}}}^{0}\cdot\nabla_{x}\dot{p}^{0}$, the macroscopic stationary
flow represented by ${\textbf{{w}}}^{0}$ influences the anisotropy of the wave
propagation through the coefficients $\mathcal{A}$ describing an anisotropic
diffusion by virtue of (3.8) and (3.9). Let us note that also for a static
fluid, $\bar{\textbf{{w}}}\equiv 0$, although coefficients $\mathcal{A}$
describe homogenized properties of an isotropic diffusion governed by the
Laplace operator, the propagation is anisotropic in general due to the pore
geometry given by $Y_{f}$.

####  3.2.4 Plane wave propagation in the homogenized medium

Acoustic plane waves propagating in the homogenized infinite medium are
expressed using the usual ansatz, defined in terms of the (constant) amplitude
$\bar{p}$, the wave direction n and the wave number $\varkappa$ (the wave
vector can be established,
${\mbox{\boldmath$\kappa$\unboldmath}}=\varkappa{\textbf{{n}}}$),

$\begin{split}p^{0}(x,t)&:=\bar{p}e^{{\rm{i}}\omega
t}e^{-{\rm{i}}\varkappa{\textbf{{n}}}\cdot{\textbf{{x}}}}\;.\\\ \end{split}$
(3.15)

Substituted in (3.14), the following eigenvalue problem is obtained: For a
given $\omega\in{\mathbb{R}}$ find $\varkappa\in{\mathbb{C}}$, such that

$\begin{split}{\textbf{{n}}}\otimes{\textbf{{n}}}:{\mbox{\boldmath$\mathcal{A}$\unboldmath}}\varkappa^{2}+2\varkappa\omega\theta{\textbf{{w}}}^{0}\cdot{\textbf{{n}}}-\phi_{f}\omega^{2}=0\;.\end{split}$
(3.16)

Using abbreviations,

$a_{n}={\textbf{{n}}}\otimes{\textbf{{n}}}:{\mbox{\boldmath$\mathcal{A}$\unboldmath}}\;,\quad
b_{n}=\theta{\textbf{{w}}}^{0}\cdot{\textbf{{n}}}\;,$ (3.17)

alternative relationships $\omega\mapsto\varkappa_{1,2}$ and
$\varkappa\mapsto\omega_{1,2}$ are established,

$\begin{split}\varkappa_{1,2}&=-\frac{\omega}{a_{n}}(b_{n}\mp\sqrt{b_{n}^{2}+a_{n}\phi_{f}})\;,\\\
\omega_{1,2}&=\frac{\varkappa}{\phi}\left(b_{n}\mp\sqrt{b_{n}^{2}+a_{n}\phi_{f}}\right)\;.\end{split}$
(3.18)

Since $a_{n}>0$ due to the positive definiteness of $\mathcal{A}$, (3.18)
presents two waves of different sound speeds $c_{1,2}=\omega/\varkappa_{1,2}$
(for a given $\omega$), propagating in mutually opposite directions. When
$\bar{\textbf{{w}}}\equiv 0$, $b_{n}=0$ and a unique wave number (up to a
sign) is obtained,
$\varkappa=\omega\sqrt{\phi_{f}/a_{n}}=\kappa_{f}c_{f}/\sqrt{a_{n}/\phi_{f}}$,
where $\kappa_{f}=\omega/c_{f}$ denotes the wave number of waves in a free
stationary fluid. Then the effective sound speed of the static fluid in the
rigid scaffolds is given by
$c_{*}=\sqrt{({\textbf{{n}}}\otimes{\textbf{{n}}}:{\mbox{\boldmath$\mathcal{A}$\unboldmath}})/\phi_{f}}$.

###  3.3 Homogenization of a viscous fluid

We consider homogenization of the linearized problem arising from (2.4) which
is now rewritten to respect its dependence on the scale parameter
$\varepsilon$. For any time $t>0$, assuming zero initial conditions, thus
$\tilde{\textbf{{u}}}^{\varepsilon}(x,t=0)=\textbf{{0}}$, and
$\tilde{p}^{\varepsilon}=0$ in $\Omega_{f}^{\varepsilon}$, the acoustic
fluctuations $\tilde{\textbf{{u}}}^{\varepsilon}$ and
$\tilde{p}^{\varepsilon}$ satisfy the following equations:

$\begin{split}\rho_{0}\left(\frac{\partial{\tilde{\textbf{{u}}}^{\varepsilon}}}{\partial{t}}+\bar{\textbf{{w}}}^{\varepsilon}\cdot\nabla\tilde{\textbf{{u}}}^{\varepsilon}+\tilde{\textbf{{u}}}^{\varepsilon}\cdot\nabla\bar{\textbf{{w}}}^{\varepsilon}\right)+\nabla\tilde{p}^{\varepsilon}-\nabla\cdot{{\rm
I}\kern-1.99997pt{\rm
D}}^{\varepsilon}{\textbf{{e}}}({\tilde{\textbf{{u}}}^{\varepsilon}})&={\textbf{{f}}}^{\varepsilon}\;,\quad\mbox{
in }\Omega_{f}^{\varepsilon}\;,\\\
\frac{\partial{\tilde{p}^{\varepsilon}}}{\partial{t}}+\bar{\textbf{{w}}}^{\varepsilon}\cdot\nabla\tilde{p}^{\varepsilon}+k_{f}\cdot\nabla\cdot{\textbf{{u}}}^{\varepsilon}&=0\;,\quad\mbox{
in }\Omega_{f}^{\varepsilon}\;,\\\
\tilde{\textbf{{u}}}^{\varepsilon}&={\textbf{{u}}}^{\partial,\varepsilon}\;,\quad\mbox{
on }\partial_{\rm{ext}}\Omega_{f}^{\varepsilon}\;,\\\
\tilde{\textbf{{u}}}^{\varepsilon}&=\textbf{{0}}\;,\quad\mbox{ on
}\Gamma_{\it{fs}}^{\varepsilon}\;,\end{split}$ (3.19)

where ${\textbf{{u}}}^{\partial,\varepsilon}(x,t)$ is a given function
describing incident waves; we assume existence of an extension
$\tilde{\textbf{{u}}}^{\partial,\varepsilon}$ of
${\textbf{{u}}}^{\partial,\varepsilon}$ to $\Omega_{f}^{\varepsilon}$, such
that $\tilde{\textbf{{u}}}^{\partial,\varepsilon}=0$ on
$\Gamma_{fs}^{\varepsilon}$. The following spaces and the set
$\hat{V}^{\varepsilon}(\Omega_{f}^{\varepsilon},t)$ will be used:

$\begin{split}V_{0}(\Omega_{f}^{\varepsilon})&=\\{{\textbf{{v}}}\in{\bf{H}}^{1}(\Omega_{f}^{\varepsilon})|\;{\textbf{{v}}}=\textbf{{0}}\;,\mbox{
on }\partial\Omega_{f}^{\varepsilon}\\}\;,\\\
\hat{V}^{\varepsilon}(\Omega_{f}^{\varepsilon},t)&=\\{{\textbf{{v}}}\in{\bf{H}}^{1}(\Omega_{f}^{\varepsilon})|\;{\textbf{{v}}}={\textbf{{v}}}^{\partial,\varepsilon}(\cdot,t)\;,\mbox{
on
}\partial\Omega_{f}^{\varepsilon}\;,\quad{\textbf{{v}}}=\textbf{{0}}\;,\mbox{
on }\Gamma_{fs}^{\varepsilon}\\}\;,\\\
Q(\Omega_{f}^{\varepsilon})&=L^{2}(\Omega_{f}^{\varepsilon})\;,\\\
Q^{1}(\Omega_{f}^{\varepsilon})&=H^{1}(\Omega_{f}^{\varepsilon})\;.\end{split}$
(3.20)

Alternatively
$\hat{V}^{\varepsilon}(\Omega_{f}^{\varepsilon},t)=V_{0}(\Omega_{f}^{\varepsilon})+\tilde{\textbf{{u}}}^{\partial,\varepsilon}(\cdot,t)$.
Below we suppress the time dependence in the notation, so that we use
$\hat{V}^{\varepsilon}(\Omega_{f}^{\varepsilon})$.

The weak formulation of the acoustic problem arises from problem (3.19)
whereby the steady state flow $\bar{\textbf{{w}}}^{\varepsilon}$ is assumed to
be known. The fluctuating velocity and pressure fields
$(\widetilde{{\textbf{{u}}}}^{\varepsilon},\widetilde{p}^{\varepsilon})$, such
that for a.a. $t>0$,
$\widetilde{{\textbf{{u}}}}^{\varepsilon}(t,\cdot)\in\hat{V}^{\varepsilon}(\Omega_{f}^{\varepsilon})$,
and $\widetilde{p}^{\varepsilon}(t,\cdot)\in Q^{1}(\Omega_{f}^{\varepsilon})$,
satisfy

$\begin{split}\rho_{0}\int_{\Omega_{f}^{\varepsilon}}{\textbf{{v}}}^{\varepsilon}\cdot\left(\dot{\widetilde{{\textbf{{u}}}}}^{\varepsilon}+\bar{\textbf{{w}}}^{\varepsilon}\cdot\nabla\widetilde{{\textbf{{u}}}}^{\varepsilon}+\widetilde{{\textbf{{u}}}}^{\varepsilon}\cdot\nabla\bar{\textbf{{w}}}^{\varepsilon}\right)&\\\
-\int_{\Omega_{f}^{\varepsilon}}\widetilde{p}^{\varepsilon}\nabla\cdot{\textbf{{v}}}^{\varepsilon}+\int_{\Omega_{f}^{\varepsilon}}{{\rm
I}\kern-1.99997pt{\rm
D}}^{f,\varepsilon}{\textbf{{e}}}({\widetilde{{\textbf{{u}}}}^{\varepsilon}}):{\textbf{{e}}}({{\textbf{{v}}}^{\varepsilon}})&=\int_{\Omega_{f}^{\varepsilon}}\tilde{\textbf{{f}}}^{f}\cdot{\textbf{{v}}}^{\varepsilon}\;,\quad\forall{\textbf{{v}}}^{\varepsilon}\in
V_{0}^{\varepsilon}(\Omega_{f}^{\varepsilon})\;,\\\
\int_{\Omega_{f}^{\varepsilon}}q^{\varepsilon}\left(\dot{\widetilde{p}}^{\varepsilon}+\bar{\textbf{{w}}}^{\varepsilon}\cdot\nabla\widetilde{p}^{\varepsilon}+k_{f}\nabla\cdot\widetilde{{\textbf{{u}}}}^{\varepsilon}\right)&=0\;,\quad\forall
q^{\varepsilon}\in Q^{0}(\Omega_{f}^{\varepsilon})\;.\end{split}$ (3.21)

We pursue the formal procedure of the asymptotic analysis
$\varepsilon\rightarrow 0$ based on the formal unfolding method. The truncated
asymptotic expansions and the viscosity scaling are adopted, however, in
contrast with the steady state flow problem, see e.g. [6, 24], a different
scaling of the viscosity and of the advection velocity must be used,

$\begin{split}\mu^{\varepsilon}&=\varepsilon^{2}\bar{\mu}\;,\quad{{\rm
I}\kern-1.99997pt{\rm D}}^{f,\varepsilon}=\varepsilon^{2}2\bar{\mu}({{\rm
I}\kern-1.99997pt{\rm
I}}-\frac{1}{3}{\textbf{{I}}}\otimes{\textbf{{I}}})\;,\\\
\quad\quad{\mathcal{T}}_{\varepsilon}{\left({\bar{\textbf{{w}}}^{\varepsilon}}\right)}&=\varepsilon\bar{\textbf{{w}}}(x,y)\;,\end{split}$
(3.22)

where ${{\rm I}\kern-1.99997pt{\rm I}}$ is the fourth-order identity tensor
and $\bar{\textbf{{w}}}(x,y)$ is the two-scale solution of steady state flow
problem. The unfolded solutions
$(\widetilde{{\textbf{{u}}}}^{\varepsilon},\widetilde{{\textbf{{w}}}}^{\varepsilon},\widetilde{p}^{\varepsilon})$
are represented by the truncated expansions

$\begin{split}{\mathcal{T}}_{\varepsilon}{\left({\widetilde{{\textbf{{u}}}}^{\varepsilon}(x,t)}\right)}&=\hat{\textbf{{u}}}(x,y,t)\;,\\\
{\mathcal{T}}_{\varepsilon}{\left({\widetilde{p}^{\varepsilon}(x,t)}\right)}&=p^{0}(x,t)+\varepsilon
p^{1}(x,y,t)\;.\end{split}$ (3.23)

The test functions ${\textbf{{v}}}^{\varepsilon}$ and $q^{\varepsilon}$ are
considered in the analogous form, thus

$\begin{split}{\mathcal{T}}_{\varepsilon}{\left({\widetilde{{\textbf{{v}}}}^{\varepsilon}(x)}\right)}&=\hat{\textbf{{v}}}(x,y)\;,\\\
{\mathcal{T}}_{\varepsilon}{\left({q^{\varepsilon}(x)}\right)}&=q^{0}(x)+\varepsilon
q^{1}(x,y)\;.\end{split}$ (3.24)

All the two-scale functions are $Y$-periodic in the second variable $y$ and
for almost all $t>0$

$\begin{split}\hat{\textbf{{u}}}(\cdot,t),\hat{\textbf{{v}}}&\in
L^{2}(\Omega;{\bf{H}}_{\\#0}^{1}(Y_{f}))\;,\\\ p^{1}(\cdot,t),q^{1}&\in
L^{2}(\Omega;H_{\\#}^{1}(Y_{f}))\;,\quad p^{0}(\cdot,t)\in
H^{1}(\Omega)\;.\end{split}$ (3.25)

Note that the non-slip condition of $\hat{\textbf{{u}}}$ on $\Gamma_{\it{fs}}$
is imposed due to the space
${\bf{H}}_{\\#0}^{1}(Y_{f})\subset{\bf{H}}_{\\#}^{1}(Y_{f})$ introduced as
follows, see Section 3.1,

${\bf{H}}_{\\#0}^{1}(Y_{f})=\\{{\textbf{{v}}}\in{\bf{H}}_{\\#}^{1}(Y_{f})|\;{\textbf{{v}}}=\textbf{{0}}\mbox{
on }\partial Y_{f}\setminus\partial Y\\}\;.$ (3.26)

######  Remark  3.

It is important to identify the limit boundary conditions for the homogenized
flow. To respect a prescribed velocity, in the limit, the Dirichlet condition
of the micromodel transforms to a Neumann-type condition associated with the
pressure. We assume the existence of ${\textbf{{u}}}^{\partial}\in
L^{2}(\partial\Omega)$ such that
$\bar{\phi}_{f}{\textbf{{u}}}^{\partial}={\rm{trace}}_{\Omega}\left({\phi_{f}{\textbf{{u}}}^{0}}\right)$,
and

$\begin{split}\int_{\partial\Omega_{f}^{\varepsilon}}\varphi{\textbf{{u}}}^{\partial,\varepsilon}\cdot{\mbox{\boldmath$\nu$\unboldmath}}\rightarrow\int_{\partial\Omega}\varphi\bar{\phi}_{f}{\textbf{{u}}}^{\partial}\cdot{\mbox{\boldmath$\nu$\unboldmath}}\;,\quad\forall\varphi\in
L^{2}(\partial\Omega)\;,\end{split}$ (3.27)

where $\bar{\phi}_{f}$ is the surface porosity while $\phi_{f}\in
H^{1}(\Omega)$ is the volume porosity. In D, we show
${\textbf{{u}}}^{0}\cdot{\mbox{\boldmath$\nu$\unboldmath}}={\textbf{{u}}}^{\partial}\cdot{\mbox{\boldmath$\nu$\unboldmath}}$.
This identifies the desired boundary condition for the limit macroscopic
velocity ${\textbf{{u}}}^{0}$ which, however, is expressed in terms of the
Darcy law involving $\nabla p^{0}$.

The volume forces fluctuations $\tilde{\textbf{{f}}}^{f}$ are assumed to be
independent of the scale $\varepsilon$, such that the associated virtual
powers converge, as follows,

$\begin{split}\int_{\Omega_{f}^{\varepsilon}}\tilde{\textbf{{f}}}^{f}\cdot\widetilde{{\textbf{{v}}}}^{\varepsilon}&\rightarrow\int_{\Omega}{\widehat{\textbf{{f}}}}^{f}\cdot\sim\kern-11.99998pt\int_{Y_{f}}\hat{\textbf{{v}}}\;.\end{split}$
(3.28)

###  3.4 Limit problem

The limit analysis of (3.19) for $\varepsilon\rightarrow 0$ leads to the
following 2-scale problem: For a.a. $t\in]0,T]$, find
$(\hat{\textbf{{u}}},p^{0},p^{1})$ such that $\hat{\textbf{{u}}}(t,\cdot)\in
L^{2}(\Omega;{\bf{H}}_{\\#0}^{1}(Y_{f}))$, $p^{1}\in L^{2}(\Omega\times
Y_{f})$, $p^{0}\in H^{1}(\Omega)$, and

$\begin{split}\rho_{0}\sim\kern-11.99998pt\int_{Y_{f}}\left(\frac{\partial{\hat{\textbf{{u}}}}}{\partial{t}}+\bar{\textbf{{w}}}\cdot\nabla_{y}\hat{\textbf{{u}}}+\hat{\textbf{{u}}}\cdot\nabla_{y}\bar{\textbf{{w}}}\right)\cdot{\textbf{{v}}}&\\\
+\sim\kern-11.99998pt\int_{Y_{f}}\bar{{\rm I}\kern-1.99997pt{\rm
D}}{\textbf{{e}}}_{y}({\hat{\textbf{{u}}}}):{\textbf{{e}}}_{y}({{\textbf{{v}}}})+\sim\kern-11.99998pt\int_{Y_{f}}(\nabla_{x}p^{0})\cdot{\textbf{{v}}}-\sim\kern-11.99998pt\int_{Y_{f}}p^{1}\nabla_{y}\cdot{\textbf{{v}}}&={\textbf{{f}}}\cdot\sim\kern-11.99998pt\int_{Y_{f}}{\textbf{{v}}}\;,\quad\forall{\textbf{{v}}}\in{\bf{H}}_{\\#0}^{1}(Y_{f})\;,\\\
\sim\kern-11.99998pt\int_{Y_{f}}q\nabla_{y}\cdot\hat{\textbf{{u}}}&=0\;,\quad\forall
q\in L^{2}(Y_{f})\;,\\\
\phi_{f}\frac{\partial{p^{0}}}{\partial{t}}+k_{f}\nabla_{x}\cdot{\textbf{{U}}}&=0\quad\mbox{
in }\Omega\;,\\\ \quad\mbox{ where
}{\textbf{{U}}}&=\sim\kern-11.99998pt\int_{Y_{f}}\hat{\textbf{{u}}}\;,\end{split}$
(3.29)

with initial and boundary conditions to be specified.

##### Characteristic responses

For decoupling the scales, the problem is transformed by the Laplace
transformation $\mathcal{L}\\{{a}\\}\mapsto{\underaccent{{*}}{{a}}}$, so that
the following split can be defined:

$\begin{split}{\underaccent{{*}}{{\hat{\textbf{{u}}}}}}&=\lambda{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{k}\rho_{0}^{-1}({\underaccent{{*}}{{f}}}_{k}-\partial_{k}^{x}{\underaccent{{*}}{{p^{0}}}})\;,\\\
{\underaccent{{*}}{{p^{1}}}}&=\lambda{\underaccent{{*}}{{\pi}}}^{k}({\underaccent{{*}}{{f}}}_{k}-\partial_{k}^{x}{\underaccent{{*}}{{p^{0}}}})\;.\end{split}$
(3.30)

We shall employ the following notation:

$\begin{split}\left\langle{{\textbf{{A}}}},\,{{\textbf{{B}}}}\right\rangle_{Y_{f}}&=\sim\kern-11.99998pt\int_{Y_{f}}{\textbf{{A}}}:{\textbf{{B}}}\;,\\\
a_{f}\left({{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)&=\left\langle{\nabla_{y}{\textbf{{u}}}},\,{\nabla_{y}{\textbf{{v}}}}\right\rangle_{Y_{f}}\;,\\\
\tilde{a}_{f}\left({{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)&=2\left\langle{{\textbf{{e}}}_{y}({{\textbf{{u}}}})-\frac{1}{3}\nabla_{y}\cdot{\textbf{{u}}}},\,{{\textbf{{e}}}_{y}({{\textbf{{v}}}})-\frac{1}{3}\nabla_{y}\cdot{\textbf{{v}}}}\right\rangle_{Y_{f}}\\\
&=2\left\langle{\widetilde{{\textbf{{e}}}_{y}({{\textbf{{u}}}})}},\,{\widetilde{{\textbf{{e}}}_{y}({{\textbf{{v}}}})}}\right\rangle_{Y_{f}}\;,\\\
b_{f}\left({\bar{\textbf{{w}}}};\,{{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)&=2\left\langle{\bar{\textbf{{w}}}\otimes{\textbf{{u}}}},\,{{\textbf{{e}}}_{y}({{\textbf{{v}}}})}\right\rangle_{Y_{f}}\;,\\\
c_{f}\left({\bar{\textbf{{w}}}};\,{{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)&=\left\langle{\bar{\textbf{{w}}}\cdot\nabla_{y}{\textbf{{u}}}},\,{{\textbf{{v}}}}\right\rangle_{Y_{f}}\;,\\\
g_{f}\left({\bar{\textbf{{w}}}};\,{{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)&=\left\langle{{\textbf{{u}}}\cdot\nabla_{y}\bar{\textbf{{w}}}},\,{{\textbf{{v}}}}\right\rangle_{Y_{f}}\;,\end{split}$
(3.31)

where
$\widetilde{{\textbf{{e}}}}={\textbf{{e}}}-(1/3){\textbf{{e}}}:{\textbf{{I}}}$
is the deviatoric part of the strain tensor.

Note that for $Y$-periodic fields ${\textbf{{u}}},{\textbf{{v}}}$ vanishing on
$\Gamma$ and satisfying $\nabla_{y}\cdot{\textbf{{u}}}=0$ and
$\nabla_{y}\cdot{\textbf{{v}}}=0$, it holds that
$a_{f}\left({{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)=\tilde{a}_{f}\left({{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)$.
Using the notation introduced above, (3.29)1 can be rewritten alternatively,
as follows:

$\begin{split}\left\langle{\frac{\partial{\hat{\textbf{{u}}}}}{\partial{t}}},\,{{\textbf{{v}}}}\right\rangle_{Y_{f}}-b_{f}\left({\bar{\textbf{{w}}}};\,{\hat{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)+\frac{\bar{\mu}}{\rho_{0}}\tilde{a}_{f}\left({\hat{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)-\frac{1}{\rho_{0}}\left\langle{p^{1}},\,{\nabla_{y}\cdot{\textbf{{v}}}}\right\rangle_{Y_{f}}&=\frac{1}{\rho_{0}}(f_{k}-\partial_{k}^{x}p^{0})\left\langle{1},\,{v_{k}}\right\rangle_{Y_{f}}\;,\\\
\left\langle{\frac{\partial{\hat{\textbf{{u}}}}}{\partial{t}}},\,{{\textbf{{v}}}}\right\rangle_{Y_{f}}+\frac{1}{2}\left(c_{f}\left({\bar{\textbf{{w}}}};\,{\hat{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)-c_{f}\left({\bar{\textbf{{w}}}};\,{{\textbf{{v}}}},\,{\hat{\textbf{{u}}}}\right)\right)+g_{f}\left({\bar{\textbf{{w}}}};\,{\hat{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)&\\\
+\frac{\bar{\mu}}{\rho_{0}}\tilde{a}_{f}\left({\hat{\textbf{{u}}}},\,{{\textbf{{v}}}}\right)-\frac{1}{\rho_{0}}\left\langle{p^{1}},\,{\nabla_{y}\cdot{\textbf{{v}}}}\right\rangle_{Y_{f}}&=\frac{1}{\rho_{0}}(f_{k}-\partial_{k}^{x}p^{0})\left\langle{1},\,{v_{k}}\right\rangle_{Y_{f}}\;.\end{split}$
(3.32)

Now, (3.32) can be transformed and (3.30) substituted, which yields a problem
for
${\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{k}\in{\bf{H}}_{\\#0}^{1}(Y_{f})$
and ${\underaccent{{*}}{{\pi}}}^{k}\in L^{2}(Y_{f})$, such that

$\begin{split}\lambda\left\langle{{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{k}},\,{{\textbf{{v}}}}\right\rangle_{Y_{f}}-b_{f}\left({\bar{\textbf{{w}}}};\,{{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{k}},\,{{\textbf{{v}}}}\right)+\frac{1}{{\rm{Re}\,}}\tilde{a}_{f}\left({{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{k}},\,{{\textbf{{v}}}}\right)-\left\langle{{\underaccent{{*}}{{\pi}}}^{k}},\,{\nabla_{y}\cdot{\textbf{{v}}}}\right\rangle_{Y_{f}}&=\frac{1}{\lambda}\left\langle{1},\,{v_{k}}\right\rangle_{Y_{f}}\quad\forall{\textbf{{v}}}\in{\bf{H}}_{\\#0}^{1}(Y_{f})\;,\\\
\left\langle{q},\,{\nabla_{y}\cdot{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{k}}\right\rangle_{Y_{f}}&=0\quad\forall
q\in L^{2}(Y_{f})\;,\end{split}$ (3.33)

where the $\frac{1}{{\rm{Re}\,}}=\bar{\mu}/\rho_{0}$ is the Reynolds number;
note that, for a given scale parameter
$\varepsilon_{0}>0,$$\frac{1}{{\rm{Re}\,}}=\mu^{\rm{phy}}/(\varepsilon_{0}^{2}\rho_{0})$
can be expressed by the physical viscosity.

By virtue of (3.29)3,4 which are the macroscopic equations, the effective
velocity is computed using the permeability which is defined by (note
$\bar{\textbf{{w}}}(x,y)$ is a two-scale function, in general)

$\begin{split}{\underaccent{{*}}{{\mathcal{K}}}}_{ij}(\lambda,x)&=\frac{\lambda}{\rho_{0}}\sim\kern-11.99998pt\int_{Y_{f}}{\underaccent{{*}}{{w}}}_{i}^{j}=\frac{\lambda}{\rho_{0}}\left\langle{1},\,{{\underaccent{{*}}{{w}}}_{i}^{j}}\right\rangle_{Y_{f}}\\\
&=\frac{1}{\rho_{0}}\left(\lambda^{3}\left\langle{{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{i}},\,{{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{j}}\right\rangle_{Y_{f}}-\lambda^{2}b_{f}\left({\bar{\textbf{{w}}}};\,{{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{i}},\,{{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{j}}\right)+\lambda\frac{1}{{\rm{Re}\,}}\tilde{a}_{f}\left({{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{i}},\,{{{\mbox{\boldmath$\chi$\unboldmath}}}^{j}}\right)\right)\;.\end{split}$
(3.34)

It is possible to define

$\begin{split}{\underaccent{{*}}{{\widehat{\mathcal{K}}}}}_{ij}(\lambda,x)&=\lambda^{2}\sim\kern-11.99998pt\int_{Y_{f}}{\underaccent{{*}}{{w}}}_{i}^{j}=\lambda\rho_{0}{\underaccent{{*}}{{\mathcal{K}}}}_{ij}(\lambda,x)\\\
&=\lambda^{4}\left\langle{{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{i}},\,{{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{j}}\right\rangle_{Y_{f}}-\lambda^{3}b_{f}\left({\bar{\textbf{{w}}}};\,{{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{i}},\,{{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{j}}\right)+\lambda^{2}\frac{1}{{\rm{Re}\,}}\tilde{a}_{f}\left({{\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{i}},\,{{{\mbox{\boldmath$\chi$\unboldmath}}}^{j}}\right)\;,\end{split}$
(3.35)

which will be used in the the plane wave propagation analysis.

It is possible compute ${\mbox{\boldmath$\chi$\unboldmath}}^{k}$ by solving
problem (3.33) transformed in the time domain: for a.a. $t\in]0,\infty[$ let
${\mbox{\boldmath$\chi$\unboldmath}}^{k}(t,\cdot)\in{\bf{H}}_{\\#0}^{1}(Y_{f})$
and $\pi^{k}(t,\cdot)\in L^{2}(Y_{f})$ satisfy

$\begin{split}\frac{\partial{}}{\partial{t}}\left\langle{{\mbox{\boldmath$\chi$\unboldmath}}^{k}},\,{{\textbf{{v}}}}\right\rangle_{Y_{f}}-b_{f}\left({\bar{\textbf{{w}}}};\,{{{\mbox{\boldmath$\chi$\unboldmath}}}^{k}},\,{{\textbf{{v}}}}\right)+\frac{1}{{\rm{Re}\,}}\tilde{a}_{f}\left({{{\mbox{\boldmath$\chi$\unboldmath}}}^{k}},\,{{\textbf{{v}}}}\right)-\left\langle{{\pi}^{k}},\,{\nabla_{y}\cdot{\textbf{{v}}}}\right\rangle_{Y_{f}}&=H(t)\left\langle{1},\,{v_{k}}\right\rangle_{Y_{f}}\quad\forall{\textbf{{v}}}\in{\bf{H}}_{\\#0}^{1}(Y_{f})\;,\\\
\left\langle{q},\,{\nabla_{y}\cdot{{\mbox{\boldmath$\chi$\unboldmath}}}^{k}}\right\rangle_{Y_{f}}&=0\quad\forall
q\in L^{2}(Y_{f})\;,\\\
{\mbox{\boldmath$\chi$\unboldmath}}^{k}(0,y)&=\textbf{{0}}\quad\mbox{ for a.a.
}y\in Y_{f}\;,\end{split}$ (3.36)

where $H(t)$ is the Heaviside function (since $t>0$, $H(t)=1$ in the right
hand side term). We assume zero initial conditions
$\hat{\textbf{{u}}}(t=0,x,y)=\textbf{{0}}$ in $\Omega\times Y_{f}$. Therefore,
the macroscopic velocity is given by

$\begin{split}{\textbf{{U}}}(t,x)&=\int_{0}^{t}\frac{{\rm{d}}\,}{{\rm{d}}\,t}\left(\sim\kern-11.99998pt\int_{Y_{f}}{\mbox{\boldmath$\chi$\unboldmath}}^{k}(t-s,y)\mathrm{d}y\right)(f_{k}(s,x)-\partial_{k}^{x}p^{0}(s,x))\mathrm{d}s\;.\end{split}$
(3.37)

Then, the permeability can be evaluated, as follows:

$\begin{split}{\mathcal{K}}_{ij}(t)&=\frac{1}{\rho_{0}}\frac{{\rm{d}}\,}{{\rm{d}}\,t}\sim\kern-11.99998pt\int_{Y_{f}}\chi_{i}^{j}(t,y)\mathrm{d}y\;.\end{split}$
(3.38)

Note the physical dimensions: $[w_{i}^{j}]=\rm{s}$, hence
$[{\mathcal{K}}_{ij}]=[\rho_{0}^{-1}][w_{i}^{j}]=\frac{\rm{m}^{3}\cdot\rm{s}}{\rm{kg}}=\frac{\rm{m}^{2}}{\rm{Pa}\cdot\rm{s}^{2}}$.

###  3.5 Macroscopic model — time domain

As announced earlier, the macroscopic model is given by (3.29)3,4. Using the
Darcy law

$\begin{split}{\underaccent{{*}}{{{\textbf{{U}}}}}}(\lambda,x)&={\underaccent{{*}}{{{\mbox{\boldmath$\mathcal{K}$\unboldmath}}}}}({\underaccent{{*}}{{{\textbf{{f}}}}}}-\nabla_{x}{\underaccent{{*}}{{p^{0}}}})\;,\\\
{\textbf{{U}}}(t,x)&=\int_{0}^{t}{\mbox{\boldmath$\mathcal{K}$\unboldmath}}(t-s)({\textbf{{f}}}(s)-\nabla_{x}p^{0}(s))\mathrm{d}s\;,\end{split}$
(3.39)

the macroscopic acoustic problems reads: Given $p^{0}(t=0,\cdot)$ by the
initial condition and ${\textbf{{f}}}(t,\cdot)$ in $\Omega$, for a.a. $t>0$
find $p^{0}(t,x)$, such that

$\begin{split}\phi_{f}\frac{\partial{p^{0}}}{\partial{t}}+k_{f}\nabla_{x}\cdot\int_{0}^{t}{\mbox{\boldmath$\mathcal{K}$\unboldmath}}(t-s)({\textbf{{f}}}(s)-\nabla_{x}p^{0}(s))\mathrm{d}s&=0\quad\mbox{
in }\Omega\;,\\\
{\mbox{\boldmath$\nu$\unboldmath}}\cdot{\textbf{{U}}}&={\mbox{\boldmath$\nu$\unboldmath}}\cdot{\textbf{{u}}}^{\partial}\quad\mbox{
on }\partial\Omega\;,\end{split}$ (3.40)

see Remark 3, where U is expressed by (3.39)2.

##### Macroscopic model — frequency domain

In the frequency domain ($\lambda={\rm{i}}\omega$), (3.29)3,4 transform to

$\begin{split}{\rm{i}}\omega\phi_{f}{\underaccent{{\circ}}{{p^{0}}}}+\frac{k_{f}}{{\rm{i}}\omega\rho_{0}}\nabla_{x}\cdot{\underaccent{{\circ}}{{\widehat{{\mbox{\boldmath$\mathcal{K}$\unboldmath}}}}}}({\rm{i}}\omega)({\underaccent{{\circ}}{{{\textbf{{f}}}}}}-\nabla_{x}{\underaccent{{\circ}}{{p^{0}}}})&=0\quad\mbox{
in }\Omega\;.\end{split}$ (3.41)

Recalling the notation formerly introduced, $\kappa_{f}=\omega/c_{f}$
designating the wave number associated with the phase velocity of wave
propagation $c_{f}=\sqrt{k_{f}/\rho_{0}}$ in the unconfined (free) fluid, the
ansatz (3.15) for the pressure plane wave yields the wave number of the
pressure wave in the porous medium,

$\begin{split}\varkappa&=\kappa_{f}\sqrt{\frac{\phi_{f}}{{\underaccent{{\circ}}{{\widehat{{\mbox{\boldmath$\mathcal{K}$\unboldmath}}}}}}({\rm{i}}\omega):{\textbf{{n}}}\otimes{\textbf{{n}}}}}=\sqrt{\frac{-{\rm{i}}\omega\phi_{f}}{k_{f}{\underaccent{{\circ}}{{{\mbox{\boldmath$\mathcal{K}$\unboldmath}}}}}({\rm{i}}\omega):{\textbf{{n}}}\otimes{\textbf{{n}}}}}\;,\end{split}$
(3.42)

where the alternative expressions correspond to the permeabilities introduced
in (3.34) and (3.35), respectively.

##  4 Floquet-Bloch wave decomposition for analysis of fluid acoustics

We shall consider plane waves propagating in an infinite porous medium, such
that $\Omega_{f}^{\varepsilon}$ is unbounded, being generated as a periodic
structure by the representative scaffold cell
$\mathcal{Z}_{f}^{\varepsilon}\subset\mathcal{Z}^{\varepsilon}$; we recall
that $\mathcal{Z}^{\varepsilon}=\varepsilon Y$ is the real size of the
periodic cell $Y$. Since in the following analysis the scale is fixed,
$\varepsilon>0$, for the sake of brevity, we shall drop the superscript
$\varepsilon$ which is related to a given scale of the structure, thus,
$\mathcal{Z}_{f}$ is used to refer to $\mathcal{Z}_{f}^{\varepsilon}$, etc.

By virtue of the Floquet-Bloch theory, the wave response of the flow models
introduced in Section 2 is represented by functions $\tilde{\textbf{{u}}}$ and
$\tilde{p}$ expressed in the decomposed form,

$\begin{split}\tilde{\textbf{{u}}}(x,t)&={\textbf{{u}}}(x)e^{-{\rm{i}}{\mbox{\boldmath$\kappa$\unboldmath}}\cdot{\textbf{{x}}}}e^{{\rm{i}}\omega
t}\;,\\\
\tilde{p}(x,t)&=p(x)e^{-{\rm{i}}{\mbox{\boldmath$\kappa$\unboldmath}}\cdot{\textbf{{x}}}}e^{{\rm{i}}\omega
t}\;,\end{split}$ (4.1)

where ${\mbox{\boldmath$\kappa$\unboldmath}}=\varkappa{\textbf{{n}}}$ is the
wave vector given by the direction of wave propagation n, and the wave number
$\varkappa$. Functions u and $p$ are $\mathcal{Z}^{\varepsilon}$-periodic and,
on the channel walls $\partial_{s}\mathcal{Z}_{f}^{\varepsilon}$, satisfy the
boundary conditions according to the flow model.

To employ the weak formulations for derivation of the dispersion
relationships, it is convenient to consider the test functions, such that they
describe waves propagating in the opposite directions with respect to those of
(4.1), thus

$\begin{split}\tilde{\textbf{{v}}}(x,t)&={\textbf{{v}}}(x)e^{{\rm{i}}{\mbox{\boldmath$\kappa$\unboldmath}}\cdot{\textbf{{x}}}}e^{{\rm{i}}\omega
t}\;,\\\
\tilde{q}(x,t)&=q(x)e^{{\rm{i}}{\mbox{\boldmath$\kappa$\unboldmath}}\cdot{\textbf{{x}}}}e^{{\rm{i}}\omega
t}\;,\end{split}$ (4.2)

where ${\textbf{{v}}},q$ are $\mathcal{Z}$-periodic. As a consequence, the
strain tensor applied to $\tilde{\textbf{{u}}}$, or to another function of the
same form, is expressed in terms of the wave-strain
${\textbf{{g}}}^{{\textbf{{n}}}}({\tilde{\textbf{{u}}}})$, thus,

$\begin{split}{\textbf{{e}}}({\tilde{\textbf{{u}}}})=\left({\textbf{{e}}}({{\textbf{{u}}}})-{\rm{i}}\varkappa{\textbf{{g}}}^{{\textbf{{n}}}}({{\textbf{{u}}}})\right)e^{-{\rm{i}}{\mbox{\boldmath$\kappa$\unboldmath}}\cdot{\textbf{{x}}}}e^{{\rm{i}}\omega
t}\;,\quad\mbox{ where
}{\textbf{{g}}}^{{\textbf{{n}}}}({{\textbf{{u}}}})=\frac{1}{2}({\textbf{{u}}}\otimes{\textbf{{n}}}+{\textbf{{n}}}\otimes{\textbf{{u}}})\;.\end{split}$
(4.3)

In analogy,
${\textbf{{e}}}({\tilde{\textbf{{v}}}})=\left({\textbf{{e}}}({{\textbf{{v}}}})+{\rm{i}}\varkappa{\textbf{{g}}}^{{\textbf{{n}}}}({{\textbf{{v}}}})\right)e^{{\rm{i}}{\mbox{\boldmath$\kappa$\unboldmath}}\cdot{\textbf{{x}}}}e^{{\rm{i}}\omega
t}$.

###  4.1 Floquet-Bloch wave decomposition for inviscid flows

The plane wave propagation in the porous medium is governed by (2.12),
involving the two parameters $\theta=2$ and $\zeta=3$. In analogy with (3.1),
the weak formulation is obtained which describes the wave propagation in the
representative cell $\mathcal{Z}_{f}$,

$\begin{split}\int_{\mathcal{Z}_{f}}\tilde{q}(\ddot{\tilde{p}}+\theta\partial_{w}\dot{\tilde{p}})-\int_{\mathcal{Z}_{f}}\partial_{w}\tilde{q}(\theta\dot{\tilde{p}}+\zeta\partial_{w}{\tilde{p}})+c_{f}^{2}\int_{\mathcal{Z}_{f}}\nabla\tilde{p}\cdot\nabla\tilde{q}&=\mathcal{I}_{\partial\mathcal{Z}_{f}}(\tilde{p},q)\;,\end{split}$
(4.4)

where $\mathcal{I}_{\partial\mathcal{Z}_{f}}(\tilde{p},q)$ is defined
according to (3.2). For the wave ansatz (4.1)2 and (4.2)2 substituted in
(4.4), there the right-hand side integral
$\mathcal{I}_{\partial\mathcal{Z}_{f}}(\tilde{p},q)$ vanishes. The velocity
field is confined by the non-penetration and the free-slip conditions, so that
${\textbf{{u}}}\cdot{\mbox{\boldmath$\nu$\unboldmath}}^{f}=0$ is prescribed on
$\partial_{s}\mathcal{Z}_{f}$. As a consequence, for almost any time $t$, the
admissible acoustic pressure $p(\cdot,t)$ must satisfy

$\begin{split}p(\cdot,t)\in P_{\\#}(\mathcal{Z}_{f}):=\\{q\in
H_{\\#}^{1}(\mathcal{Z}_{f})|q\mbox{ is $\mathcal{Z}$-periodic,
}\nabla(qe^{-{\rm{i}}{\mbox{\boldmath$\kappa$\unboldmath}}\cdot{\textbf{{x}}}})\cdot{\mbox{\boldmath$\nu$\unboldmath}}^{f}=0\mbox{
on }\partial_{s}\mathcal{Z}_{f}\\}\;.\end{split}$ (4.5)

Note that, in the definition of $P_{\\#}(\mathcal{Z}_{f}^{\varepsilon})$, the
zero normal-projected gradient concerns $\tilde{p}$ involving the exponential
part and not only $p$, see (4.1).

To find nontrivial solutions to (4.4), the following dispersion problem must
be solved: Find $p\in P_{\\#}(\mathcal{Z}_{f})$ and
$\varkappa\in{\mathbb{C}}$, such that

$\begin{split}-\omega^{2}\int_{\mathcal{Z}_{f}}pq+{\rm{i}}\omega\theta\left(\int_{\mathcal{Z}_{f}}\bar{\textbf{{w}}}\cdot[\nabla
p-{\rm{i}}\varkappa{\textbf{{n}}}p]q-\int_{\mathcal{Z}_{f}}\bar{\textbf{{w}}}\cdot[\nabla
q+{\rm{i}}\varkappa{\textbf{{n}}}q]p\right)\\\
-\zeta\int_{\mathcal{Z}_{f}}\bar{\textbf{{w}}}\cdot[\nabla
p-{\rm{i}}\varkappa{\textbf{{n}}}p][\nabla
q+{\rm{i}}\varkappa{\textbf{{n}}}q]\cdot\bar{\textbf{{w}}}+c_{f}^{2}\int_{\mathcal{Z}_{f}}[\nabla
p-{\rm{i}}\varkappa{\textbf{{n}}}p]\cdot[\nabla
q+{\rm{i}}\varkappa{\textbf{{n}}}q]=0\;,\end{split}$ (4.6)

for all $q\in P_{\\#}(\mathcal{Z}_{f})$. Alternatively, for a given
$\kappa\in{\mathbb{R}}$, frequency $\omega\in{\mathbb{C}}$ can be computed

##### FEM – discretization of the Bloch wave analysis

Problem (4.6) is solved numerically using the finite element method (FEM). The
matrices resulting from the conforming FEM discretization of integral forms
involved in (4.6) are introduced below, employing a self-explaining notation:

$\begin{split}\int_{\mathcal{Z}_{f}}\nabla p\cdot\nabla
q&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}={\bf{q}}^{T}{\bf{C}}{\bf{p}}\;,\\\
\int_{\mathcal{Z}_{f}}\nabla p\cdot{\textbf{{n}}}q&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}={\bf{q}}^{T}{\bf{Y}}{\bf{p}}\;,\\\
\int_{\mathcal{Z}_{f}}pq&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{q}}^{T}{\bf{M}}{\bf{p}}\;,\\\
\zeta\int_{\mathcal{Z}_{f}}\bar{\textbf{{w}}}\cdot\nabla
p\bar{\textbf{{w}}}\cdot\nabla q&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{q}}^{T}{\bf{G}}_{w}{\bf{p}}\;,\\\
\zeta\int_{\mathcal{Z}_{f}}\bar{\textbf{{w}}}\cdot\nabla
p(\bar{\textbf{{w}}}\cdot{\textbf{{n}}})q&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{q}}^{T}{\bf{R}}_{w}{\bf{p}}\;,\\\
\zeta\int_{\mathcal{Z}_{f}}(\bar{\textbf{{w}}}\cdot{\textbf{{n}}})p(\bar{\textbf{{w}}}\cdot{\textbf{{n}}})q&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{q}}^{T}{\bf{W}}_{w}{\bf{p}}\;,\\\
\theta\int_{\mathcal{Z}_{f}}(\bar{\textbf{{w}}}\cdot{\textbf{{n}}})pq&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{q}}^{T}{\bf{V}}_{w}{\bf{p}}\;,\\\
\theta\int_{\mathcal{Z}_{f}}q(\bar{\textbf{{w}}}\cdot\nabla
p)&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{q}}^{T}{\bf{N}}_{w}{\bf{p}}\;.\end{split}$ (4.7)

Above the column matrices ${\bf{p}}$ and ${\bf{q}}$ represent all degrees of
freedom (DOFs) arising from the discretization of $p$ and $q$, respectively,
considering $\mathcal{Z}_{f}$, whereas the DOFs reduction due to the periodic
conditions on the external channel surfaces
$\partial_{\\#}\mathcal{Z}_{f}=\partial\mathcal{Z}\cap\partial\mathcal{Z}_{f}$
is applied. Using this notation, identity (4.6) is approximated by the
following equation,

$\begin{split}\left[c_{f}^{2}{\bf{C}}+(c_{f}^{2}\varkappa^{2}-\omega^{2}){\bf{M}}+\varkappa^{2}{\bf{W}}_{w}+2\varkappa\omega{\bf{V}}_{w}-{\bf{G}}_{w}\right]{\bf{p}}\\\
+{\rm{i}}\left[\omega({\bf{N}}_{w}-{\bf{N}}_{w}^{T})-\varkappa({\bf{R}}_{w}-{\bf{R}}_{w}^{T})-\varkappa
c_{f}^{2}({\bf{Y}}-{\bf{Y}}^{T})\right]{\bf{p}}=0\;.\end{split}$ (4.8)

Two alternative eigenvalue problems can be considered: a) for a given real
frequency $\omega$, find a complex $\varkappa\in{\mathbb{C}}$ which satisfies
(4.8), or b) for a given real $\varkappa\in{\mathbb{R}}$, find a complex
frequency $\omega\in{\mathbb{C}}$ satisfying (4.8). Both these formulations
lead to a quadratic eigenvalue problem (QEP). To transform it into a standard
eigenvalue problem, we employ the following notations,

$\begin{split}{\bf{Q}}&:={\bf{R}}_{w}+c_{f}^{2}{\bf{Y}}\;,\\\
{\bf{H}}&:={\bf{G}}_{w}-c_{f}^{2}{\bf{C}}\;,\\\
{\bf{U}}&:={\bf{W}}_{w}+c_{f}^{2}{\bf{M}}\;.\end{split}$ (4.9)

Then, using the Choleski decomposition, auxiliary variables ${\bf{s}}$ and
${\bf{r}}$ can be introduced, such that

$\begin{split}{\bf{U}}&={\bf{T}}_{U}^{T}{\bf{T}}_{U}\;,\quad\varkappa^{2}{\bf{U}}{\bf{p}}=\varkappa{\bf{T}}_{U}^{T}{\bf{s}}\;,\quad\varkappa{\bf{T}}_{U}{\bf{p}}={\bf{s}}\;,\\\
{\bf{M}}&={\bf{T}}_{M}^{T}{\bf{T}}_{M}\;,\quad\omega^{2}{\bf{M}}{\bf{p}}=\omega{\bf{T}}_{M}^{T}{\bf{r}}\;,\quad{\bf{T}}_{M}{\bf{p}}={\bf{r}}\;.\end{split}$
(4.10)

Let us first consider the problem for computing $\varkappa$ for a given
$\omega$. Using the decomposition in (4.10)1, the identity (4.8) can be
transformed into the following standard eigenvalue problem for computing
$\lambda_{\varkappa}:=1/\varkappa$ which satisfy,

$\begin{split}\left[\begin{array}[]{ll}{\bf{A}}_{\varkappa}(\omega),&{\bf{T}}_{U}^{T}\\\
{\bf{T}}_{U},&\bf{0}\end{array}\right]\left[\begin{array}[]{c}{\bf{p}}\\\
{\bf{s}}\end{array}\right]=\lambda_{\kappa}\left[\begin{array}[]{ll}{\bf{B}}_{\varkappa}(\omega),&\bf{0}\\\
\bf{0},&{\bf{I}}\end{array}\right]\left[\begin{array}[]{c}{\bf{p}}\\\
{\bf{s}}\end{array}\right]\;,\\\ \mbox{ where
}{\bf{A}}_{\varkappa}(\omega)=2\omega{\bf{V}}_{w}-{\rm{i}}({\bf{Q}}-{\bf{Q}}^{T})\;,\quad{\bf{B}}_{\varkappa}(\omega)=\omega^{2}{\bf{M}}+{\bf{H}}+{\rm{i}}\omega({\bf{N}}_{w}^{T}-{\bf{N}}_{w})\;.\end{split}$
(4.11)

The reason for computing $\lambda_{\varkappa}$ (which is proportional to the
wavelength) rather than $\varkappa$ arises form the regularity of the right
hand side block matrix of this generalized eigenvalue problem, since the left
hand side matrix may become nearly singular. This may happen because of
vanishing ${\bf{V}}_{w}$, when $\bar{\textbf{{w}}}\cdot{\textbf{{n}}}\approx
0$, thus, when the advection flow field $\bar{\textbf{{w}}}$ is almost
orthogonal to the wave direction n.

Alternatively, $\omega$ can be computed in response to $\varkappa$. For this,
in analogy, we employ the decomposition in (4.10)2. Using the notation (4.9),
the following generalized eigenvalue problem can be established which yields a
real $\lambda_{\omega}:=1/\omega$,

$\begin{split}\left[\begin{array}[]{ll}{\bf{A}}_{\omega}(\varkappa),&-{\bf{T}}_{M}^{T}\\\
-{\bf{T}}_{M},&\bf{0}\end{array}\right]\left[\begin{array}[]{c}{\bf{p}}\\\
{\bf{r}}\end{array}\right]=\lambda_{\omega}\left[\begin{array}[]{ll}{\bf{B}}_{\omega}(\varkappa),&\bf{0}\\\
\bf{0},&-{\bf{I}}\end{array}\right]\left[\begin{array}[]{c}{\bf{p}}\\\
{\bf{r}}\end{array}\right]\;,\\\ \mbox{ where
}{\bf{A}}_{\omega}(\varkappa)=2\varkappa{\bf{V}}_{w}+{\rm{i}}({\bf{N}}_{w}-{\bf{N}}_{w}^{T})\;,\quad{\bf{B}}_{\omega}(\varkappa)={\bf{H}}-\varkappa^{2}{\bf{U}}+{\rm{i}}{\varkappa}({\bf{Q}}-{\bf{Q}}^{T})\;.\end{split}$
(4.12)

As in the case of (4.11), the right hand side matrix is regular even for
vanishing $\bar{\textbf{{w}}}$, whereas the left hand side matrix may be not.

######  Remark  4.

All matrices defined in (4.7) are real, whereby
${\bf{C}},{\bf{M}},{\bf{G}},{\bf{V}}$ and ${\bf{W}}$are symmetric, thereby
both the block matrices ${\bf{A}}$ and ${\bf{B}}$ in the QEPs problems (4.11)
and (4.12) are Hermitian. Therefore, their eigenvalues $\lambda$ are real, or
come in complex-conjugate pairs, see e.g. [29]. The definiteness of the
diagonal blocks and, therefore, of the whole block matrices in both the QEP
depend on the respective parameters, i.e. on $\omega$ in (4.11), and on
$\varkappa$ in (4.12). Moreover, all matrices labelled by subscript w depend
on the advection and matrices ${\bf{R}}_{w}$, ${\bf{V}}_{w}$ and
${\bf{W}}_{w}$ may become nearly singular, when
$\bar{\textbf{{w}}}\cdot{\textbf{{n}}}\approx 0$. However, as we demonstrate
in the numerical examples reported in Section 5, up to rather high advection
velocities $w^{0}=|{\textbf{{w}}}^{0}|$, the pressure modes corresponding to
predictions obtained by the homogenization approach are featured by real
eigenvalues $\lambda\in{\mathbb{R}}$, as expected by virtue of the inviscid
fluid properties.

###  4.2 Floquet-Bloch wave decomposition for viscous fluids

In analogy with the Bloch wave analysis for the model of inviscid fluids in an
infinite periodic porous structure, due to the wave decomposition ansatz
(4.1), we consider the reduced problem imposed in the periodic cell
$\mathcal{Z}_{f}$. We consider a steady flow represented by velocity field
$\bar{\textbf{{w}}}$ in the scaffold periodic cell $\mathcal{Z}_{f}$. The weak
formulation for $(\tilde{\textbf{{u}}}(x),\tilde{p}(x))$ given by the ansatz
(4.1) is derived from (2.4) upon multiplying (2.4)1 by $\tilde{\textbf{{v}}}$
and integrating over $\mathcal{Z}_{f}$; this yields

$\begin{split}\int_{\mathcal{Z}_{f}}\rho_{0}\left(\frac{{\rm{d}}\,}{{\rm{d}}\,t}\tilde{\textbf{{u}}}+\bar{\textbf{{w}}}\cdot\nabla\tilde{\textbf{{u}}}+\tilde{\textbf{{u}}}\cdot\nabla\bar{\textbf{{w}}}\right)\cdot\tilde{\textbf{{v}}}&\\\
+\int_{\mathcal{Z}_{f}}{{\rm I}\kern-1.99997pt{\rm
D}}{\textbf{{e}}}({\tilde{\textbf{{u}}}}):{\textbf{{e}}}({\tilde{\textbf{{v}}}})-\int_{\mathcal{Z}_{f}}\tilde{p}\nabla\cdot\tilde{\textbf{{v}}}&=\int_{\partial\mathcal{Z}_{f}}({{\rm
I}\kern-1.99997pt{\rm
D}}({\textbf{{e}}}({\tilde{\textbf{{u}}}})-\tilde{p}{\textbf{{I}}}):({\mbox{\boldmath$\nu$\unboldmath}}^{f}\otimes\tilde{\textbf{{v}}})\;,\\\
\int_{\mathcal{Z}_{f}}\tilde{q}\left(\frac{{\rm{d}}\,}{{\rm{d}}\,t}\tilde{p}+\bar{\textbf{{w}}}\cdot\nabla\tilde{p}+k_{f}\nabla\cdot\tilde{\textbf{{u}}}\right)&=0\;,\end{split}$
(4.13)

to hold for all test functions $(\tilde{\textbf{{v}}},\tilde{q})$ in the form
(4.2) which describes waves running in the opposite direction with respect to
the waves $(\tilde{\textbf{{u}}},\tilde{p})$. Above, the right hand side
boundary integral vanishes due to the $\mathcal{Z}$-periodicity and vanishing
${\textbf{{u}}},{\textbf{{v}}}$ on the solid-fluid interface.

######  Remark  5.

Due to $\nabla\cdot\bar{\textbf{{w}}}=0$, the following integration by parts
can be employed in (4.13),

$\begin{split}\int_{\mathcal{Z}_{f}}q\bar{\textbf{{w}}}\cdot\nabla\tilde{p}&=\frac{1}{2}\int_{\mathcal{Z}_{f}}\left(q\bar{\textbf{{w}}}\cdot\nabla\tilde{p}-\tilde{p}\bar{\textbf{{w}}}\cdot\nabla
q\right)\;,\\\
\int_{\mathcal{Z}_{f}}\bar{\textbf{{w}}}\cdot\nabla\tilde{\textbf{{u}}}\cdot{\textbf{{v}}}&=\frac{1}{2}\int_{\mathcal{Z}_{f}}\left(\bar{\textbf{{w}}}\cdot\nabla\tilde{\textbf{{u}}}\cdot{\textbf{{v}}}-\bar{\textbf{{w}}}\cdot\nabla{\textbf{{v}}}\cdot\tilde{\textbf{{u}}}\right)\;.\end{split}$
(4.14)

Now the wave ansatz (4.1) is substituted in (4.13), so that (4.13) leads to
the following dispersion eigenvalue problem: Given $\varkappa\in{\mathbb{R}}$,
find a complex $\omega\in{\mathbb{C}}$, such that

$\begin{split}\rho_{0}\int_{\mathcal{Z}_{f}}\left({\textbf{{u}}}\cdot\nabla\bar{\textbf{{w}}}+\bar{\textbf{{w}}}\cdot\nabla{\textbf{{u}}}-{\rm{i}}\varkappa({\textbf{{n}}}\cdot\bar{\textbf{{w}}}){\textbf{{u}}}\right)\cdot{\textbf{{v}}}-\int_{\mathcal{Z}_{f}}(\nabla\cdot{\textbf{{v}}}+{\rm{i}}\varkappa{\textbf{{n}}}\cdot{\textbf{{v}}})p&\\\
+\int_{\mathcal{Z}_{f}}\left({{\rm I}\kern-1.99997pt{\rm
D}}{\textbf{{e}}}({{\textbf{{u}}}}):{\textbf{{e}}}({{\textbf{{v}}}})+\varkappa^{2}{{\rm
I}\kern-1.99997pt{\rm
D}}{\textbf{{g}}}^{{\textbf{{n}}}}({{\textbf{{u}}}}):{\textbf{{g}}}^{{\textbf{{n}}}}({{\textbf{{v}}}})\right)&\\\
+{\rm{i}}\varkappa\int_{\mathcal{Z}_{f}}\left({{\rm I}\kern-1.99997pt{\rm
D}}{\textbf{{e}}}({{\textbf{{u}}}}):{\textbf{{g}}}^{{\textbf{{n}}}}({{\textbf{{v}}}})-{{\rm
I}\kern-1.99997pt{\rm
D}}{\textbf{{g}}}^{{\textbf{{n}}}}({{\textbf{{u}}}}):{\textbf{{e}}}({{\textbf{{v}}}})\right)&=-{\rm{i}}\omega\rho_{0}\int_{\mathcal{Z}_{f}}{\textbf{{u}}}\cdot{\textbf{{v}}}\;,\quad\\\
\gamma\int_{\mathcal{Z}_{f}}(\bar{\textbf{{w}}}\cdot\nabla
p-{\rm{i}}\varkappa{\textbf{{n}}}\cdot\bar{\textbf{{w}}}p)q+\int_{\mathcal{Z}_{f}}q(\nabla\cdot{\textbf{{u}}}-{\rm{i}}\varkappa{\textbf{{n}}}\cdot{\textbf{{u}}})&=-{\rm{i}}\omega\gamma\int_{\mathcal{Z}_{f}}pq\;,\end{split}$
(4.15)

for all ${\textbf{{v}}}\in{\bf{H}}_{\\#0}^{1}(\mathcal{Z}_{f})$ and for all
$q\in H_{\\#}^{1}(\mathcal{Z}_{f})$.

##### FEM approximation and matrix notation.

In analogy with the inviscid case, the FEM model is obtained upon discretizing
problem (4.15). The non-slip condition for the velocity is applied on
$\partial_{s}\mathcal{Z}_{f}$ along with the periodicity conditions on
$\partial_{\\#}\mathcal{Z}_{f}$. In addition to column matrices ${\bf{p}}$ and
${\bf{q}}$ representing the pressure field, column matrices ${\bf{u}}$ and
${\bf{v}}$ involve all degrees of freedom (DOFs) arising form the
discretization of u and v, respectively, respecting the boundary conditions.

The matrices obtained by the FE discretization of (4.15) are defined below in
(4.16) and (4.17) — some of the matrix symbols are reused with different
meanings with respect those introduced above in Section 4.1.

The first group of matrices arises from the discretized integrals which do not
depend on the advection velocity,

$\begin{split}\int_{\mathcal{Z}_{f}}{{\rm I}\kern-1.99997pt{\rm
D}}{\textbf{{e}}}({{\textbf{{u}}}}):{\textbf{{e}}}({{\textbf{{v}}}})&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{v}}^{T}{\bf{A}}{\bf{u}}\;,\\\ \int_{\mathcal{Z}_{f}}{{\rm
I}\kern-1.99997pt{\rm
D}}{\textbf{{g}}}^{{\textbf{{n}}}}({{\textbf{{u}}}}):{\textbf{{g}}}^{{\textbf{{n}}}}({{\textbf{{v}}}})&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{v}}^{T}{\bf{S}}{\bf{u}}\;,\\\ \int_{\mathcal{Z}_{f}}{{\rm
I}\kern-1.99997pt{\rm
D}}{\textbf{{e}}}({{\textbf{{u}}}}):{\textbf{{g}}}^{{\textbf{{n}}}}({{\textbf{{v}}}})&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{v}}^{T}{\bf{X}}{\bf{u}}\;,\\\
\int_{\mathcal{Z}_{f}}pq&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{q}}^{T}{\bf{Q}}{\bf{p}}\;,\\\
\int_{\mathcal{Z}_{f}}{\textbf{{u}}}\cdot{\textbf{{v}}}&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{v}}^{T}{\bf{M}}{\bf{u}}\;,\\\
\int_{\mathcal{Z}_{f}}q\nabla\cdot{\textbf{{u}}}&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{q}}^{T}{\bf{B}}{\bf{u}}\;,\\\
\int_{\mathcal{Z}_{f}}q{\textbf{{n}}}\cdot{\textbf{{u}}}&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{q}}^{T}{\bf{N}}_{n}{\bf{u}}\;,\\\ \end{split}$ (4.16)

whereas the second group of matrices depend on $\bar{\textbf{{w}}}$; the
respective matrices vanish for the static fluid, when
$\bar{\textbf{{w}}}\equiv 0$,

$\begin{split}\int_{\mathcal{Z}_{f}}({\textbf{{n}}}\cdot\bar{\textbf{{w}}})pq&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{q}}^{T}{\bf{Q}}_{w}{\bf{p}}\;,\\\
\int_{\mathcal{Z}_{f}}({\textbf{{n}}}\cdot\bar{\textbf{{w}}}){\textbf{{u}}}\cdot{\textbf{{v}}}&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{v}}^{T}{\bf{M}}_{w}{\bf{u}}\;,\\\
\int_{\mathcal{Z}_{f}}q\bar{\textbf{{w}}}\cdot\nabla
p&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}\frac{1}{2}{\bf{q}}^{T}({\bf{Y}}_{w}-{\bf{Y}}_{w}^{T}){\bf{p}}\;,\\\
\int_{\mathcal{Z}_{f}}(\bar{\textbf{{w}}}\cdot\nabla{\textbf{{u}}})\cdot{\textbf{{v}}}&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}\frac{1}{2}{\bf{v}}^{T}({\bf{C}}_{w}-{\bf{C}}_{w}^{T}){\bf{u}}\;,\\\
\int_{\mathcal{Z}_{f}}({\textbf{{u}}}\cdot\nabla\bar{\textbf{{w}}})\cdot{\textbf{{v}}}&\stackrel{{\scriptstyle\rm
FEM}}{{\approx}}{\bf{v}}^{T}{\bf{G}}_{w}{\bf{u}}\;,\end{split}$ (4.17)

where matrices ${\bf{Y}}_{w}$ and ${\bf{C}}_{w}$ were introduced by virtue of
the expressions (4.14). Note that ${\bf{G}}_{w}$ is non-symmetric. Using this
notation, the discretized form of (4.15) can be written, as

$\begin{split}\rho_{0}[{\rm{i}}\omega{\bf{M}}-{\rm{i}}\varkappa{\bf{M}}_{w}+\frac{1}{2}({\bf{C}}_{w}-{\bf{C}}_{w}^{T})+{\bf{G}}_{w}]{\bf{u}}-({\bf{B}}^{T}+{\rm{i}}\varkappa{\bf{N}}_{n}^{T}){\bf{p}}&\\\
+\left[{\bf{A}}+{\rm{i}}\varkappa({\bf{X}}-{\bf{X}}^{T})+\varkappa^{2}{\bf{S}}\right]{\bf{u}}=\textbf{0}\;,\\\
[{\rm{i}}\omega{\bf{Q}}+\frac{1}{2}({\bf{Y}}_{w}-{\bf{Y}}_{w}^{T})-{\rm{i}}\varkappa{\bf{Q}}_{w}]{\bf{p}}+k_{f}({\bf{B}}-{\rm{i}}\varkappa{\bf{N}}_{n})\cdot{\bf{u}}&=\textbf{0}\;.\end{split}$
(4.18)

######  Remark  6.

Formally one can consider the degenerate case of an inviscid fluid, $\mu=0$,
so that matrices ${\bf{A}},{\bf{S}}$ and ${\bf{X}}$ vanish, however, the
irrotational velocity constraint $\nabla\times\tilde{\textbf{{u}}}=0$ should
be imposed to suppress spurious modes, and also the free-slip, non-penetration
conditions must be applied, i.e.
${\mbox{\boldmath$\nu$\unboldmath}}^{f}\cdot{\textbf{{u}}}=0$ on
$\partial_{s}\mathcal{Z}_{f}$.

Two alternative formulations of the dispersion problem arising from (4.18) can
be considered: 1) For a given real frequency $\omega\in{\mathbb{R}}$, compute
a complex wave number $\varkappa\in{\mathbb{C}}$, such that (4.18) holds with
nonvanishing eigenvectors, ${\bf{u}},{\bf{p}}$; 2) For a given real wave
number $\varkappa\in{\mathbb{R}}$, compute a complex frequency
$\omega\in{\mathbb{C}}$ and nonvanishing ${\bf{u}},{\bf{p}}$;

In this paper we use the first formulation, although it leads to the quadratic
eigenvalue problem (QEP) which requires to introduce an additional variable
and to use a substitution enabling to derive a standard generalized eigenvalue
problem. The second alternative seems to be advantageous, since no such
substitutions are needed, however, it appears to be quite cumbersome to
distinguish the correct pressure modes among all the dispersion curves.
Nevertheless, even this second formulation is outlined briefly in Section
4.2.2.

####  4.2.1 Formulation $\omega\mapsto\varkappa$

We consider the problem: For a given $\omega\in{\mathbb{R}}$, find
$\varkappa\in{\mathbb{C}}$, such that (4.18) admits a non-trivial solution
$({\bf{u}},{\bf{p}})$. It is advisable to transformed this QEP into a problem
involving $\varkappa$ in the linear form. To do so, the Choleski decomposition
of ${\bf{S}}={\bf{V}}^{T}{\bf{V}}$ is applied, which yields the substitution
${\bf{z}}:={\rm{i}}\varkappa{\bf{V}}{\bf{u}}$ and
$\varkappa^{2}{\bf{S}}{\bf{u}}=-{\rm{i}}\varkappa{\bf{V}}^{T}{\bf{z}}$.

The following block matrices can be established, where $\gamma=1/k_{f}$,

$\begin{split}{{\bf I}\kern-1.99997pt{\bf
P}}_{w}:=\left(\begin{array}[]{lll}{\bf{A}}+\rho_{0}[\frac{1}{2}({\bf{C}}_{w}-{\bf{C}}_{w}^{T})+{\rm{i}}\omega{\bf{M}}+{\bf{G}}_{w}]\;,&-{\bf{B}}^{T}\;,&\textbf{0}\\\
{\bf{B}}\;,&{\rm{i}}\omega\gamma{\bf{Q}}+\frac{\gamma}{2}({\bf{Y}}_{w}-{\bf{Y}}_{w}^{T})\;,&\textbf{0}\\\
\textbf{0}\;,&\textbf{0}\;,&{\bf{I}}\end{array}\right)\;,\\\ {{\bf
I}\kern-1.99997pt{\bf
N}}_{w}:=\left(\begin{array}[]{lll}\rho_{0}{\bf{M}}_{w}+({\bf{X}}^{T}-{\bf{X}})\;,&{\bf{N}}_{n}^{T}\;,&{\bf{V}}^{T}\\\
{\bf{N}}_{n}\;,&\gamma{\bf{Q}}_{w}\;,&\textbf{0}\\\
{\bf{V}}\;,&\textbf{0}\;,&\textbf{0}\end{array}\right)\;,\quad\quad{{\bf
v}\kern-5.0pt{\bf v}}=\left(\begin{array}[]{l}{\bf{u}}\\\ {\bf{p}}\\\
{\bf{z}}\end{array}\right)\;.\end{split}$ (4.19)

Now (4.18) can be rewritten,

$\begin{split}\frac{1}{{\rm{i}}\varkappa}{{\bf I}\kern-1.99997pt{\bf
P}}_{w}{{\bf v}\kern-5.0pt{\bf v}}={{\bf I}\kern-1.99997pt{\bf N}}_{w}{{\bf
v}\kern-5.0pt{\bf v}}\;,\end{split}$ (4.20)

where ${{\bf I}\kern-1.99997pt{\bf P}}_{w}$ is regular. Unfortunately, because
of the combination of the viscosity and advection effects associated with
velocity $\bar{\textbf{{w}}}$, no special properties can be identified
concerning matrices ${{\bf I}\kern-1.99997pt{\bf P}}_{w}$ and ${{\bf
I}\kern-1.99997pt{\bf N}}_{w}$. Nevertheless, this formulation of the
generalized eigenvalue problem has been implemented to analyze the dispersion
$\omega\mapsto\varkappa$, as reported in Section 5.2. The eigenmodes of
interest are those with the smallest imaginary part $\varkappa_{I}$ (for a
given real $\omega$).

######  Remark  7.

(Balancing the equations) Since $\rho_{0}/\gamma\approx 10^{10}$, the system
should be balanced to improve conditioning of matrix ${{\bf
I}\kern-1.99997pt{\bf N}}$. For this we set,

$\begin{split}\tilde{\bf{u}}:=\rho_{0}^{1/2}{\bf{u}}\;,\quad\tilde{\bf{p}}:=\gamma^{1/2}{\bf{p}}\;,\\\
\tilde{\bf{B}}:=(\rho_{0}\gamma)^{-1/2}{\bf{B}}\;,\quad\tilde{\bf{N}}_{n}:=(\rho_{0}\gamma)^{-1/2}{\bf{N}}_{n}\;,\\\
\tilde{\bf{A}}:={\bf{A}}/\rho_{0}\;,\quad\tilde{\bf{S}}:={\bf{S}}/\rho_{0}\;,\quad\tilde{\bf{X}}:={\bf{X}}/\rho_{0}\;.\end{split}$
(4.21)

Using these substitutions, matrices (4.19) are replaced by

$\begin{split}\tilde{{\bf I}\kern-1.99997pt{\bf
P}}:=\left(\begin{array}[]{lll}\tilde{\bf{A}}+{\rm{i}}\omega{\bf{M}}+\frac{1}{2}({\bf{C}}_{w}-{\bf{C}}_{w}^{T})+{\bf{G}}_{w}\;,&-\tilde{\bf{B}}^{T}\;,&\textbf{0}\\\
\tilde{\bf{B}}\;,&{\rm{i}}\omega{\bf{Q}}+\frac{1}{2}({\bf{Y}}_{w}-{\bf{Y}}_{w}^{T})\;,&\textbf{0}\\\
\textbf{0}\;,&\textbf{0}\;,&{\bf{I}}\end{array}\right)\;,\\\ \tilde{{\bf
I}\kern-1.99997pt{\bf
N}}:=\left(\begin{array}[]{lll}{\bf{M}}_{w}+\tilde{\bf{X}}^{T}-\tilde{\bf{X}}\;,&\tilde{\bf{N}}_{n}^{T}\;,&\tilde{\bf{V}}^{T}\\\
\tilde{\bf{N}}_{n}\;,&{\bf{Q}}_{w}\;,&\textbf{0}\\\
\tilde{\bf{V}}\;,&\textbf{0}\;,&\textbf{0}\end{array}\right)\;,\quad\quad\tilde{{\bf
v}\kern-5.0pt{\bf v}}=\left(\begin{array}[]{l}\tilde{\bf{u}}\\\
\tilde{\bf{p}}\\\ \tilde{\bf{z}}\end{array}\right)\;,\end{split}$ (4.22)

where $\tilde{\bf{S}}=\tilde{\bf{V}}^{T}\tilde{\bf{V}}$, hence
$\tilde{\bf{V}}:={\bf{V}}\rho_{0}^{-1/2}$.

####  4.2.2 Formulation $\varkappa\mapsto\omega$

Alternatively, one may consider a given real wave number
$\varkappa\in{\mathbb{R}}$ and find an $\omega\in{\mathbb{C}}$, such that
(4.18) admits a nontrivial solution $(\omega,({\bf{u}},{\bf{p}}))$. In this
case, (4.18) presents a linear eigenvalue problem for the complex wave
frequency $\omega$, whose imaginary part signifies the wave attenuation. The
following matrices can be established, recalling $\gamma=1/k_{f}$,

$\begin{split}{{\bf I}\kern-1.99997pt{\bf
R}}_{w}:=\left(\begin{array}[]{ll}{\bf{A}}+\varkappa^{2}{\bf{S}}+{\rm{i}}\varkappa({\bf{X}}-{\bf{X}}^{T})\;,&\textbf{0}\\\
\textbf{0}\;,&\textbf{0}\end{array}\right)\;,\\\ {{\bf I}\kern-1.99997pt{\bf
K}}_{w}:=\left(\begin{array}[]{ll}\rho_{0}[\frac{1}{2}({\bf{C}}_{w}-{\bf{C}}_{w}^{T})-{\rm{i}}\varkappa{\bf{M}}_{w}+{\bf{G}}_{w}]\;,&-({\bf{B}}^{T}+{\rm{i}}\varkappa{\bf{N}}_{n}^{T})\\\
{\bf{B}}-{\rm{i}}\varkappa{\bf{N}}_{n}&-{\rm{i}}\varkappa\gamma{\bf{Q}}_{w}+\frac{\gamma}{2}({\bf{Y}}_{w}-{\bf{Y}}_{w}^{T})\end{array}\right)\;,\\\
{{\bf I}\kern-1.99997pt{\bf
M}}:=\left(\begin{array}[]{ll}\rho_{0}{\bf{M}}\;,&\textbf{0}\\\
\textbf{0}&\gamma{\bf{Q}}\end{array}\right)\;,\quad\quad{{\bf
v}\kern-5.0pt{\bf v}}=\left(\begin{array}[]{l}{\bf{u}}\\\
{\bf{p}}\end{array}\right)\;.\end{split}$ (4.23)

Now (4.18) can be rewritten,

$\begin{split}[{{\bf I}\kern-1.99997pt{\bf K}}_{w}+{{\bf I}\kern-1.99997pt{\bf
R}}_{w}]{{\bf v}\kern-5.0pt{\bf v}}=-{\rm{i}}\omega{{\bf I}\kern-1.99997pt{\bf
M}}{{\bf v}\kern-5.0pt{\bf v}}\;,\end{split}$ (4.24)

where ${{\bf I}\kern-1.99997pt{\bf M}}$ is real symmetric positive definite
and ${{\bf I}\kern-1.99997pt{\bf R}}_{w}$ being associated with the fluid
viscosity, is Hermitean, but singular. Matrix ${\rm{i}}{{\bf
I}\kern-1.99997pt{\bf K}}_{w}$ is general, but becomes Hermitean for the case
of a static fluid, i.e. when $\bar{\textbf{{w}}}=0$ so that the nonsymmetric
matrix ${\bf{G}}_{w}$ becomes zero.

######  Remark  8.

As the consequence of the above statements, for static and inviscid fluids,
(4.24) yields real frequencies, thus,
${\mathbb{R}}\ni\varkappa\mapsto\omega\in{\mathbb{R}}$. In contrast, when
$\bar{\textbf{{w}}}\not\equiv 0$, even for inviscid fluids when ${{\bf
I}\kern-1.99997pt{\bf R}}_{w}$ vanishes, (4.24) yields complex frequencies in
general, as discussed in Remark 4. However, in such a case, by virtue of
Remark 6, free slip boundary conditions should be respected (thus, the
displacement-associated matrices modified) and the irrotationality constraint
related to the velocity field should be involved to reduce spurious
oscillation modes.

##  5 Numerical examples

Properties and mutual correspondence of the models introduced in the preceding
sections we be illustrated using examples of 2D and 3D geometries representing
the periodic scaffolds saturated by viscous, or inviscid fluids. For this we
compute dispersion curves and phase velocities characterizing the wave
propagation in such media. Besides comparison of the responses provided by the
two modelling approaches, i.e. the periodic homogenization (PH) and the
Floquet-Bloch wave decomposition (FB), we aim to explore the following
phenomena related to the wave propagation:

  * 1.

the flow advection – the influence of the permanent fluid flow given by the
macroscopic velocity $\bar{\bm{w}}^{0}$ (its size
$\bar{w}^{0}=|\bar{\bm{w}}^{0}|$ and orientation with respect to the wave
direction $\bm{n}$;

  * 2.

the scale effects – the influence of the characteristic size of the
microstructure, as described by the scale parameter $\varepsilon^{0}$.

Numerical solutions of the characteristic problems associated with the
homogenized models, and of the eigenvalue problems, arising from the FB
analysis, were computed using the finite element method implemented in the
software SfePy [7]. In all the examples reported below, meshes with a
sufficient resolution have been used, such that any further uniform refinement
did not change the wave dispersion results for long wave lengths near the
origin of the Brillouin zone. In all figures, results of the numerical
experiments are displayed in the basic SI units, unless specified otherwise.

###  5.1 Inviscid fluid

The inviscid fluid parameters were those of water $\gamma=5\cdot 10^{-10}$
Pa-1, $\rho_{0}=1000$ kg/m3. We considered 3D robocast structures, [16],
produced by sinterization of ceramic fibres. The representative cell
$\mathcal{Z}$ is illustrated in Fig. 2 (left). This domain generates a lattice
of periodic 3D scaffolds. The largest edge of the cell $\mathcal{Z}$
($7.5\cdot 10^{-4}$ m in the $x_{1}$ direction) limits the shortest
wavelengths for the FB analysis, hence, in the reciprocal lattice, the
corresponding 1st irreducible Brillouin zone spans $8378$ m-1. The P1 finite
elements on tetrahedrons were used for the discretization of both the FB and
PH analyses.

|  $d=1.5\cdot 10^{-4}$ m |  $d=2.3\cdot 10^{-4}$ m
---|---|---
|  |
solid fibres $\mathcal{Z}_{s}$ | fluid in $\mathcal{Z}_{f}$ | fluid in $\mathcal{Z}_{f}$
Figure 2: The periodic cell of solid fibres (left, from [16]) and
corresponding finite element meshes ($l=6.5\cdot 10^{-4}$ m, $h=2.75\cdot
10^{-4}$ m) with two radii of the fibres for inviscid flow experiments.

####  5.1.1 Static fluid, $w^{0}=0$

Below we study the influence of the solid fibres diameters on the dispersion
properties using (4.8) for the Bloch wave analysis (the frequencies $\omega$
were computed for given wave numbers $\varkappa$) and (3.18) for the
homogenization-based analysis, simplified for the case of zero convective
velocity, $\bar{\textbf{{w}}}\equiv\bf{0}$. Specifically, (4.8) then presents
a generalized eigenvalue problem for $\lambda:=\omega^{2}$, which satisfies

${\bf{A}}_{0}{\bf{p}}=\omega^{2}{\bf{M}}{\bf{p}}\;\quad\mbox{ with
}{\bf{A}}_{0}:=c_{f}^{2}\left[{\bf{C}}+\varkappa^{2}{\bf{M}}-{\rm{i}}\varkappa({\bf{Y}}-{\bf{Y}}^{T})\right]\;.$
(5.1)

Due to the positive definiteness of both ${\bf{A}}_{0}$ and ${\bf{M}}$, real
positive $\omega^{2}$ can be computed for given $\varkappa\in{\mathbb{R}}$.
Examples of the finite element (FE) meshes of the 3D domains are shown for two
different fibre diameters in Fig. 2 (middle, right). The eigenvalues were
computed using the ARPACK solver (implicitly restarted Lanczos method [18])
through the SciPy package (function eigsh() [30]). The shift-invert mode was
used to accelerate the calculation of the smallest eigenvalues.

We considered incident waves in the direction $\bm{n}=(1,0,0)^{T}$ and
examined influence of the porosity. While increasing fibres radii, thus
reducing the fluid volume fraction, a partial band gap appears and increases,
as illustrated in Fig. 3 in terms of the dispersion curves
($\omega(\varkappa)$ dependence) and phase velocity plots. The dashed black
lines correspond to the homogenization-based asymptotes, see (3.18), while the
thick color lines depict results of the FB analysis, obtained by solving
(5.1).

Further, using the mesh in Fig. 5 (left), fixing the fiber diameter
$d=2.1\cdot 10^{-4}$ m, the wave vector
${\mbox{\boldmath$\kappa$\unboldmath}}=\varkappa\bm{n}$ was modified on a path
in the reciprocal lattice; we traced a part of the boundary of the Brillouin
zone, while fixing $\varkappa_{2}=0$ and varying $\varkappa_{1}$,
$\varkappa_{3}$, see Fig. 4 (left). The dispersion curves (the color lines)
resulting from the FB analysis are displayed in Fig. 4 (right), where the
dashed line represents the homogenization-based asymptotes.

$d=0.1\cdot 10^{-4}$ m
---
|
$d=1.5\cdot 10^{-4}$ m
|
$d=2.3\cdot 10^{-4}$ m
|
Figure 3: The influence of the solid fibres diameters $d$ on the dispersion properties. As the radius of the fibres increases, a band gap appears and increases. The Bloch-based (dots) and homogenization-based (lines) dispersion results are shown using the dispersion curves (left), and the phase velocity curves (right). Units: $\varkappa$: 1/mm, $\omega$: MHz. |
---|---
Figure 4: Left: a path of the wave vectors $\varkappa\bm{n}$ on the Bragg
planes, bounding the Brillouin zone: $\varkappa_{2}=0$ and $\varkappa_{1}$,
$\varkappa_{3}$ vary. Right: the dispersion curves computed by the FB analysis
(dotted color lines), and the corresponding homogenization-based predictions
(the black dashed lines). Units: [$\varkappa$] = 1/mm, [$\omega$] = MHz.

####  5.1.2 Nonstationary inviscid fluid – effects of the advection

We report the influence of the advection magnitude $\bar{w}^{0}$ (the
macroscopic flow velocity) and the incident wave direction $\bm{n}$ on the
dispersion properties. The dispersion analysis was performed for
$\bar{w}^{0}\in(0,300)$ using (4.12) for the FB analysis (the frequencies
$\omega$ were computed for given wave numbers $\varkappa$) and (3.18) for the
homogenization-based analysis. The eigenvalues were computed using the MATLAB®
function eigs() [28, 18] via the MATLAB Engine API for Python.

The macroscopic velocity $\bar{\bm{w}}^{0}=(\bar{w}^{0},0,0)$ had the overall
direction along the $x_{1}$-axis and was specified using its magnitude
$\bar{w}^{0}$. The local velocity field $\bar{\textbf{{w}}}$ was established
by virtue of the homogenization result, see C. Field $\bar{\textbf{{w}}}$ is
illustrated in Fig. 5, in terms of the streamlines distributed in the 3D
computational domain.

$d=2.1\cdot 10^{-4}$ m |
---|---
|
Figure 5: The finite element mesh (left) and inviscid flow streamlines
(right).

Two cases studies were performed for two different directions $\bm{n}$ of the
incident waves:

  * 1.

Study I-A: $\bm{n}=(0,1,0)^{T}$, i.e. perpendicular to the macroscopic flow,
thus $\bm{n}\perp\bar{\bm{w}}^{0}$.

  * 2.

Study I-B: $\bm{n}=(1,0,0)^{T}$, i.e. in the direction of the macroscopic
flow, thus $\bm{n}\parallel\bar{\bm{w}}^{0}$.

In Fig. 6, the dispersion analysis by the FB method is compared with the
homogenization-based approximation, whereby $\bar{w}^{0}\in\\{1,100,300\\}$
m/s. To compare the two possible methods of the FB analysis implementation,
namely the dispersion mappings $\omega\mapsto\varkappa$, problem (4.11), and
$\varkappa\mapsto\omega$, problem (4.12), in Fig. 7, the corresponding curves
are displayed for one fixed macroscopic advection velocity. For the case study
I-B with the advection velocity $\bar{w}^{0}=100$, responses computed using
both the FB analysis implementations, i.e. $\varkappa\mapsto\omega$ and
$\omega\mapsto\varkappa$ are displayed in Fig. 8.

In Fig. 9 we present the dependence of the homogenized coefficients $a_{n}$,
$b_{n}$ involved in (3.17) on the macroscopic advection velocity magnitude
$\bar{w}^{0}$. In Fig. 10 the dependence of the homogenization-based phase
velocities on $\bar{w}^{0}$ is shown: $c_{+}$, $-c_{-}$ in the left subplot
and $|c_{+}|-|c_{-}|$ in the right subplot. As can be seen in the both
figures, in study I-A, the influence of $\bar{w}^{0}$ magnitude is relatively
small (nevertheless increasing with increasing $\bar{w}^{0}$): $a_{n}$,
$b_{n}$ are almost constant and $c_{+}$, $-c_{-}$ relative difference stays
below $10^{-5}$. On the other hand, in study I-B, $a_{n}$, $b_{n}$ vary
significantly and large progressively increasing differences between $c_{+}$
and $-c_{-}$ can be observed.

$\bar{w}^{0}=1$
---
|
$\bar{w}^{0}=100$
|
$\bar{w}^{0}=300$
|
Figure 6: The comparison of the dispersion analysis by the FB method (color
point lines) and the homogenization-based analysis (black lines), various
macroscopic flow velocities $\bar{\bm{w}}^{0}$ we considered: Left: study I-A;
Right: study I-B. $\bar{w}^{0}=100$
---
|
Figure 7: The comparison of the dispersion analysis by the FB method (color point lines) and the homogenization-based analysis (black lines) for various inlet convective velocities $\bar{w}^{0}$, computed $\varkappa$ for given $\omega$. Left: study I-A, right: study I-B. Figure 8: The comparison of $\varkappa\rightarrow\omega$ (Fig. 6) and $\omega\rightarrow\varkappa$ (Fig. 7) dispersion results for the study I-B, $\bar{w}^{0}=100$. Only the first quadrant is shown. |
---|---
|
Figure 9: The dependence of homogenized coefficients $a_{n}$, $b_{n}$ on inlet convective velocities $\bar{w}^{0}$. Left: study I-A, right: study I-B. |
---|---
|
Figure 10: The dependence of the homogenization-based phase velocities
$c_{+}$, $-c_{-}$ on $\bar{w}^{0}$ (top) and the dependence of the difference
$|c_{+}|-|c_{-}|$ on $\bar{w}^{0}$ (bottom). Left: study I-A, right: study
I-B.

###  5.2 Viscous fluid

Because of computational hurdles associated with solving the quadratic
eigenvalue problems (QEP), we consider 2D structures only leading to smaller
discretized models . The viscous fluid parameters are those of water:
$\mu=1.02\cdot 10^{-3}$ Pa$\cdot$s, $\gamma=5\cdot 10^{-10}$ Pa-1,
$\rho_{0}=1000$ kg/m3. The P2 (for velocity variables) and P1 (for pressure
variables) finite elements on triangles were used for the discretization of
both the FB and PH analyses.

####  5.2.1 Homogenization-based analysis

By using (3.42) we studied, for a fixed finite scale parameter
$\varepsilon=10^{-5}$, the influence of the convective velocity magnitude
$\bar{w}^{0}$ on the dispersion properties. The incident wave vector direction
is $\bm{n}=(0,1)^{T}$, whereas the convective flow represented by
$\bar{\bm{w}}^{0}\parallel(1,0)^{T}$ is aligned with the $x_{1}$-axis. The
associated advection field $\bar{\textbf{{w}}}$ is obtained by the standard
reconstruction of the homogenized Stokes flow, see e.g. [25, 1, 10, 32]; for
illustration of $\bar{\textbf{{w}}}$, the streamlines are depicted in Fig 11.
Four velocities were used: $\bar{w}^{0}\in\\{0,1,5,10\\}$ m/s.

Figure 11: The finite element mesh and viscous flow velocity field.

The characteristic responses
${\underaccent{{*}}{{{\mbox{\boldmath$\chi$\unboldmath}}}}}^{k}$ defined by
(3.33) are shown in Fig. 12 for two velocities: $\bar{w}^{0}=0$ and
$\bar{w}^{0}=10$ m/s. The steady advection flow modifies the so-called
correctors (the characteristic responses) associated with the macroscopic
pressure gradients, see (3.30), thus, inducing a non-symmetry of the dynamic
permeability calculated according to (3.34). This effect is demonstrated in
Fig. 13; while for zero advection, the responses of the modes $k=1$ and $k=2$
are symmetric with respect to the coordinate axes (rotation by $\pi/2$), this
symmetry is lost for $\bar{w}^{0}\not=0$. Finally, the dispersion properties,
as depending on $\bar{w}^{0}$, in terms of the $\varkappa(\omega)$ plots and
phase velocity and attenuation plots are shown in Fig. 14.

| $k=1\quad\quad k=2$ | $k=1\quad\quad k=2$
---|---|---
$\Im{{\textbf{{w}}}^{k}}$
$\Re{{\textbf{{w}}}^{k}}$ |  |
|  $|\bar{\textbf{{w}}}|=0$ [m/s] |  $|\bar{\textbf{{w}}}|=10$ [m/s]
Figure 12: Characteristic responses ${{\textbf{{w}}}^{k}}$ for “horizontal”
flow $\bar{\bm{w}}$.

|
---|---
(a) 1st eig.val. of ${\underaccent{{\circ}}{{{\mbox{\boldmath$\mathcal{K}$\unboldmath}}}}}({\rm{i}}\omega)$, $\Re{\lambda_{1}(\omega)}$ | (b) 2nd eig.val. of ${\underaccent{{\circ}}{{{\mbox{\boldmath$\mathcal{K}$\unboldmath}}}}}({\rm{i}}\omega)$, $\Re{\lambda_{2}(\omega)}$

|
---|---
(c) 1st eig.val. of ${\underaccent{{\circ}}{{{\mbox{\boldmath$\mathcal{K}$\unboldmath}}}}}({\rm{i}}\omega)$, $\Im{\lambda_{1}(\omega)}$ | (d) 2nd eig.val. of ${\underaccent{{\circ}}{{{\mbox{\boldmath$\mathcal{K}$\unboldmath}}}}}({\rm{i}}\omega)$, $\Im{\lambda_{2}(\omega)}$

Figure 13: The influence of the convective velocity magnitude $\bar{w}^{0}$ on
the eigenvalues of the homogenized permeability tensor.

|
---|---
(a) $\Re{\varkappa(\omega)}$ | (b) $\Im{\varkappa(\omega)}$

|
---|---
(c) phase velo.: $c_{w}(\omega)=\omega/\Re{\varkappa}$ | (d) attenuation $\omega/\Im{\varkappa}$

Figure 14: The influence of the convective velocity magnitude $\bar{w}^{0}$ on
the wave number, phase velocity and attenuation.

####  5.2.2 Bloch wave analysis compared to homogenization

The homogenization-based dispersion analysis (3.42) results were compared with
the Bloch wave analysis (4.20) as follows. The incoming wave vector direction
was $\bm{n}=(0,1)^{T}$. The convective flow $\bar{\bm{w}}$ had the overall
direction along the $x_{1}$-axis $(1,0)^{T}$, and was specified using its
magnitude $\bar{w}^{0}$ on the inflow (left) boundary. It is illustrated,
together with the used finite element mesh, in Fig 15 in terms of streamlines.

Figure 15: The finite element mesh and viscous flow streamlines.

Two parametric studies of the dispersion were performed:

  * 1.

Study V-A: The finite scale parameter $\varepsilon_{0}$ was modified in the
range $[10^{-4},5.23\cdot 10^{-3}]$ in lockstep with the finite scale
advection velocity $\bar{\textbf{{w}}}^{\varepsilon_{0}}$, so that
$|\bar{\bm{w}}|=1$ in (3.33), and (3.34); recall the scaling (3.22).

  * 2.

Study V-B: The dispersion analyses were performed for the fixed finite scale
parameter $\varepsilon_{0}=10^{-3}$, while changing $\bar{w}^{0}$ in range
$[1,21.5]$.

While in the V-A study, the advection effect is being changed proportionally
to the microstructure size, the V-B study is aimed to explore the advection
effect for a fixed microstructure size.

The eigenvalues in the FB analysis were computed using the MATLAB® function
eigs() [28, 18] via the MATLAB Engine API for Python. Ten eigenvalues
($\frac{1}{i\varkappa}$) with the largest real part, i.e., $\varkappa$ with
the smallest imaginary part, were requested in the study V-A, while only three
eigenvalues in the study V-B.

The dispersion curves and phase velocities obtained by the homogenization-
based and FB-based analyses are compared in Fig. 16 for the study V-A, for two
selected finite scale parameters $\varepsilon_{0}$. Similarly, in Fig. 17, for
the study V-B, the results for two of the $\bar{w}^{0}$ are shown. The
dispersion curves $\omega\mapsto\varkappa$ computed by the homogenization-
based method were always closest to those obtained from the FB analysis with
the smallest imaginary part of eigenvalues $\varkappa$, thus, to the least
attenuated modes.

Since, in this case, the dispersion analyses by the two approaches correspond
well, only those results obtained by the homogenization method are presented.
For the study V-A, the dependence of $\varkappa(\omega)$, phase velocity
$c_{w}(\omega)=\omega/\Re{\varkappa}$ and attenuation $\omega/\Im{\varkappa}$
curves on $\varepsilon_{0}\in[10^{-4},5.23\cdot 10^{-3}]$ is shown in Fig. 18.
In analogy, for the study V-B, the dependence of the dispersion curves on
$\bar{w}^{0}\in[1,21.5]$ for the fixed $\varepsilon_{0}=10^{-3}$ is presented
in Fig. 19.

$\varepsilon_{0}=10^{-4}$
---
|
$\varepsilon_{0}=5.23\cdot 10^{-3}$
|
Figure 16: Study V-A, the comparison of dispersion curves (left) and phase
velocities (right) for homogenization (solid lines) and Bloch-based analyses
(triangles).  $\bar{w}^{0}=1$
---
|
$\bar{w}^{0}=21.5$
|
Figure 17: Study V-B, the comparison of dispersion curves (left) and phase velocities (right) for homogenization (solid lines) and Bloch-based analyses (triangles). |
---|---
|
Figure 18: Study V-A, the dependence of dispersion curves (left) and phase velocities (right) on $\varepsilon_{0}$ obtained by the homogenization analysis. |
---|---
|
Figure 19: Study V-B, the dependence of dispersion curves (left) and phase
velocities (right) on $\bar{w}^{0}$ obtained by the homogenization analysis.

##  6 Conclusion

We derived a macroscopic model governing the acoustic waves in rigid
periodically porous media perfused by permanent flows of Newtonian, or
inviscid fluids. This model extends the results reported in the literature by
considering the advection phenomenon associated with the Navier-Stokes
equations. The dynamic effects described by the model are limited by the
linearization based on the assumptions stated in Section 2.1. Besides the
thermal effects which are neglected by virtue of the restriction to barotropic
fluids, the acoustic linearization assumes the oscillations of the velocity
field magnitudes and of its gradient being small in comparison with the
permanent flow velocity field. For the inviscid fluid with considered
advection permanent flow, a reduced formulation for pressure field only was
proposed.

As the second contribution of the paper, the Floquet-Bloch (FB) theory was
applied to study the dispersion phenomena for models featured by the advection
effects for the viscous and inviscid fluids. In both considered types of
fluid, the analysis leads to a general form of the quadratic eigenvalue
problem (QEP). Linearizations of QEPs for all treated models were suggested.
Numerical solutions were obtained using the finite element method implemented
in the software SfePy [7]. The numerical studies demonstrated expected
relationships between the wave dispersion results obtained by the
homogenization-based and the FB-based analyses. While the first approach
provides a sufficiently accurate approximation for long waves only, the second
approach enables to analyze planar waves propagating in infinite media.
However, further research of numerical methods for the FB should be pursued to
improve robustness and efficiency of the computational algorithms, especially
in the case of the QEP using techniques presented, e.g., in [29].

We believe that the modelling approach presented in this paper can be
developed further towards a more detailed study of the nonlinear phenomena of
acoustics in porous media, such as the acoustic streaming, with a high
applicability in biomedical technologies and smart material design.

##### Acknowledgments

This research was supported by project GACR 17-01618S of the Scientific
Foundation of the Czech Republic and due to the European Regional Development
Fund-Project “Application of Modern Technologies in Medicine and Industry”
(No. CZ.02.1.01/0.0/ 0.0/17 048/0007280), and in part by project LO 1506 of
the Czech Ministry of Education, Youth and Sports.

##  Appendix A Model of acoustic fluctuations

Based on the Assumptions established in Section 2.1, we may consider the
fluctuations in (2.2) being proportional to a small parameter $\alpha>0$, so
that, without changing the notation, (2.2) becomes
${\textbf{{w}}}=\bar{\textbf{{w}}}+\alpha\tilde{\textbf{{w}}}$,
$p=\bar{p}+\alpha\tilde{p}$, and $\rho=\rho_{0}+\alpha\tilde{\rho}$ where
$\bar{\rho}=\rho_{0}$ is assumed to be a constant, see Assumption (A4). Upon
substituting these expansions in (2.1), at order $\alpha^{0}$, this yields the
stationary flow model (2.5). At order $\alpha^{1}$, the momentum equation
(2.1)1 yields (2.4)1 straightforwardly, while from the mass conservation
(2.1)2 we obtain

$\begin{split}{\frac{\partial{\tilde{\rho}}}{\partial{t}}}+\nabla\cdot(\rho_{0}\tilde{\textbf{{u}}}+{\tilde{\rho}\bar{\textbf{{w}}}})&=0\;.\end{split}$

Hence, using the incompressibility (2.5)2 and the assumptions (A5) related to
the barotropic approximation, i.e. $\tilde{p}=c_{f}^{2}\tilde{\rho}$, with
$c_{f}^{2}=(\gamma\rho_{0})^{-1}$, the density $\tilde{\rho}$ can be
eliminated, so that the above equation becomes

$\begin{split}\rho_{0}\gamma\left(\frac{\partial{\tilde{p}}}{\partial{t}}+\bar{\textbf{{w}}}\cdot\nabla\tilde{p}\right)+\rho_{0}\nabla\cdot\tilde{\textbf{{u}}}&=0\;,\end{split}$
(A.1)

thus, (2.4)2 holds. Obviously, terms of higher orders in the perturbation
parameter $\alpha$ could be considered, bringing further equations involving
the first-order approximation handled in this paper.

##  Appendix B Basics of the homogenization by the unfolding operator method

In Section 3, the homogenization results were obtained by the _periodic
unfolding method_ , see [9]. To introduce the periodic unfolding operator, a
domain is needed, containing the “entire” periods $\varepsilon Y$ only:

$\begin{split}\hat{\Omega}^{\varepsilon}&=\mbox{interior}\bigcup_{\zeta\in\Xi^{\varepsilon}}Y_{\zeta}^{\varepsilon}\;,\quad
Y_{\zeta}^{\varepsilon}=\varepsilon(\overline{Y}+\zeta)\;,\\\ \mbox{ where
}\Xi^{\varepsilon}&=\\{\zeta\in{\mathbb{Z}}^{3}\,|\;\varepsilon(\overline{Y}+\zeta)\subset\Omega\\}\;.\end{split}$

For simplicity we may consider such domains $\Omega$ and such subsequences
$\\{\varepsilon_{k}\\}$ for which $\Omega=\hat{\Omega}^{\varepsilon_{k}}$. For
all $z\in{\mathbb{R}}^{3}$, let $[z]$ be the unique integer such that
$z-[z]\in Y$. Since $z=[z]+\\{z\\}$ for all $z\in{\mathbb{R}}^{3}$, for all
$\varepsilon>0$, the unique decomposition holds,

$x=\varepsilon\left(\left[\frac{x}{\varepsilon}\right]+\left\\{\frac{x}{\varepsilon}\right\\}\right)=\xi+\varepsilon
y\quad\forall
x\in{\mathbb{R}}^{3}\;,\quad\xi=\varepsilon\left[\frac{x}{\varepsilon}\right]\;.$
(B.1)

Based on this decomposition, the periodic unfolding operator
$\mathcal{T}_{\varepsilon}\,{}:L^{2}(\Omega;{\mathbb{R}})\rightarrow
L^{2}(\Omega\times Y;{\mathbb{R}})$ is defined as follows: for any function
$v\in L^{1}(\Omega;{\mathbb{R}})$, extended to
$L^{1}({\mathbb{R}}^{3};{\mathbb{R}})$ by zero outside $\Omega$, i.e. $v=0$ in
${\mathbb{R}}^{3}\setminus\Omega$,

${\mathcal{T}}_{\varepsilon}{\left({v}\right)}(x,y)=\left\\{\begin{array}[]{ll}v\left(\varepsilon\displaystyle\left[\frac{x}{\varepsilon}\right]+\varepsilon
y\right)\;,&x\in\hat{\Omega}^{\varepsilon},y\in Y\;,\\\ 0&\mbox{ otherwise
}.\\\ \end{array}\right.$

For product of any $u$ and $v$ the unfolding yields
${\mathcal{T}}_{\varepsilon}{\left({uv}\right)}={\mathcal{T}}_{\varepsilon}{\left({u}\right)}{\mathcal{T}}_{\varepsilon}{\left({v}\right)}$.
The following integration formula holds:

$\int_{\hat{\Omega}^{\varepsilon}}v\,dx=\frac{1}{|Y|}\int_{\Omega\times
Y}{\mathcal{T}}_{\varepsilon}{\left({v}\right)}\,dy\,dx=\int_{\Omega}\sim\kern-11.99998pt\int_{Y}{\mathcal{T}}_{\varepsilon}{\left({v}\right)}\,dy\,dx\quad\forall
v\in L^{1}(\Omega)\;.$

##  Appendix C Local periodic advection velocity field for the inviscid fluid

Here we briefly present the homogenization result of the problem (2.6). For a
given “macroscopic” advection velocity ${\textbf{{w}}}^{0}$, which may depend
on the macroscopic position $x\in\Omega$, the fluctuation of the velocity
potential $\psi(x,\cdot)\in H_{\\#}^{1}(Y_{f})$ for almost all $x\in\Omega$
satisfies

$\begin{split}\int_{Y_{f}}\nabla_{y}\psi\cdot\nabla_{y}q={\textbf{{w}}}^{0}\cdot\int_{Y_{f}}\nabla_{y}q\quad\forall
q\in H_{\\#}^{1}(Y_{f})\;.\end{split}$ (C.1)

Note that the non-penetration condition holds, i.e.
${\mbox{\boldmath$\nu$\unboldmath}}\cdot(\nabla_{y}\psi-{\textbf{{w}}}^{0})=0$
on $\Gamma$, where ${\textbf{{w}}}^{1}:=-\nabla_{y}\psi$ is the fluctuating
part of the velocity. Then, the local advection field $\bar{\textbf{{w}}}$ is
computed, as follows:

$\begin{split}\bar{\textbf{{w}}}(x,y)={\textbf{{w}}}^{0}(x)-\nabla_{y}\psi(x,y)\;.\end{split}$
(C.2)

In the present paper, we assume periodic advection
$\bar{\textbf{{w}}}^{\varepsilon}$, so that ${\textbf{{w}}}^{0}$ is a constant
vector.

##  Appendix D Limit seepage velocity

We prove the convergence of ${\textbf{{u}}}^{\partial,\varepsilon}$ to trace
of the limit seepage velocity ${\textbf{{u}}}^{0}$ on $\partial\Omega$, see
(3.27). Let $\varphi^{\delta}\in C^{\infty}(\Omega)$. We compute limits of the
identity,

$\begin{split}\int_{\partial\Omega_{f}^{\varepsilon}}\varphi^{\delta}{\textbf{{u}}}^{\partial,\varepsilon}\cdot\bm{n}-\int_{\Omega_{f}^{\varepsilon}}\nabla\cdot{\textbf{{u}}}^{\varepsilon}\varphi^{\delta}=\int_{\Omega_{f}^{\varepsilon}}{\textbf{{u}}}^{\varepsilon}\cdot\nabla\varphi^{\delta}\;.\end{split}$
(D.1)

Clearly, the right hand side integral converges, as follows

$\begin{split}\int_{\Omega_{f}^{\varepsilon}}{\textbf{{u}}}^{\varepsilon}\cdot\nabla\varphi^{\delta}&\rightarrow\int_{\Omega}\sim\kern-11.99998pt\int_{Y_{f}}\hat{\textbf{{u}}}\cdot\nabla_{x}\varphi^{\delta}=\int_{\Omega}\nabla_{x}\varphi^{\delta}\cdot(\phi_{f}{\textbf{{u}}}^{0})\\\
&=-\int_{\Omega}\varphi^{\delta}\nabla_{x}\cdot(\phi_{f}{\textbf{{u}}}^{0})+\int_{\partial\Omega}(\phi_{f}{\textbf{{u}}}^{0})\cdot\bm{n}\;.\end{split}$
(D.2)

The first left hand side integral in (D.1) converges by virtue of the given
data ${\textbf{{u}}}^{\partial,\varepsilon}$, thus

$\begin{split}\int_{\partial\Omega_{f}^{\varepsilon}}\varphi^{\delta}{\textbf{{u}}}^{\partial,\varepsilon}\cdot\bm{n}\rightarrow\int_{\partial\Omega}\varphi^{\delta}\bar{\phi}_{f}{\textbf{{u}}}^{\partial}\cdot\bm{n}\;.\end{split}$
(D.3)

The volume integral on the left hand side in (D.1) converges, as follows

$\begin{split}-\int_{\Omega_{f}^{\varepsilon}}\nabla\cdot{\textbf{{u}}}^{\varepsilon}\varphi^{\delta}\rightarrow-\int_{\Omega}\nabla_{x}\varphi^{\delta}\cdot\sim\kern-11.99998pt\int_{Y_{f}}\hat{y}\nabla_{y}\cdot\hat{\textbf{{u}}}=0\;,\end{split}$
(D.4)

since $\nabla_{y}\cdot\hat{\textbf{{u}}}=0$ in $Y_{f}$. Above $\hat{y}$ is the
relative position of $y$ with respect to the barycenter of $Y$. Hence (D.2)
and (D.3) yield

$\begin{split}\int_{\partial\Omega}\varphi^{\delta}\bar{\phi}_{f}{\textbf{{u}}}^{\partial}\cdot\bm{n}=-\int_{\Omega}\varphi^{\delta}\nabla_{x}\cdot(\phi_{f}{\textbf{{u}}}^{0})+\int_{\partial\Omega}\varphi^{\delta}(\phi_{f}{\textbf{{u}}}^{0})\cdot\bm{n}\;.\end{split}$
(D.5)

The assertion (3.27) now follows by $\delta\rightarrow 0$.

## References

  * [1] G. Allaire.  Homogenization of the unsteady Stokes equations in porous media.  In C. Bandle, J. Bemelmans, M. Chipot, M. Grüter, and J. Saint Jean Paulin, editors, Progress in Partial Differential Equations: Calculus of Variations, Applications, volume 296 of Pitman Research Notes in Mathematics Series, pages 109–123. Longman Scientific & Technical, 1992.
  * [2] A. Bensoussan, J.L. Lions, and G. Papanicolaou.  Asymptotic methods in periodic media.  North-Holland, 1978.
  * [3] Thierry Bourbié, O. Coussy, and B. Zinszner.  Acoustics of Porous Media.  Institut français du pétrole publications. Editions Technip, 1987.
  * [4] J.M. Carcione.  Wave fields in real media. Wave propagation in anisotropic, anelastic, porous and electromagnetic media, volume 38 of Handbook of geophysical exploration. Section I. Seismic exploration.  Elsevier, third edition, extended and revised edition, 2014.
  * [5] Clifford E. Chafin.  Wave-flow interactions and acoustic streaming.  2016\.  arXiv:1602.04893.
  * [6] Z. Chen, L. Lyons, and G. Qin.  Derivation of the Forchheimer law via homogenization.  Transport in Porous Media, 44:325–335, 2001.
  * [7] Robert Cimrman, Vladimír Lukeš, and Eduard Rohan.  Multiscale finite element calculations in python using sfepy.  Advances in Computational Mathematics, 45(4):1897–1921, Aug 2019\.
  * [8] D. Cioranescu, A. Damlamian, P. Donato, G. Griso, and R. Zaki.  The periodic unfolding method in domains with holes.  SIAM Journal on Mathematical Analysis, 44(2):718–760, 2012.
  * [9] D. Cioranescu, A. Damlamian, and G. Griso.  The periodic unfolding method in homogenization.  SIAM Journal on Mathematical Analysis, 40(4):1585–1620, 2008.
  * [10] D. Cioranescu, A. Damlamian, G. Griso, and et al.  The Stokes problem in perforated domains by the periodic unfolding method. new trends in continuum mechanics.  Theta Series in Advanced Mathematics, 3:67–80, 2005.
  * [11] M. Collet, M. Ouisse, M. Ruzzene, and M.N. Ichchou.  Floquet-Bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems.  International Journal of Solids and Structures, 48(20):2837–2848, 2011.
  * [12] Jose Diaz-Alban and Nader Masmoudi.  Asymptotic analysis of acoustic waves in a porous medium: Microincompressible flow.  Communications in Partial Differential Equations, 39(11):2125–2167, 2014.
  * [13] R.P. Gilbert and A. Panchenko.  Effective acoustic equations for a two-phase medium with microstructure.  Mathematical and Computer Modelling, 39:1431–1448, 2004.
  * [14] Wu J.  Acoustic streaming and its applications.  Fluids, 3(4):108–125.
  * [15] Hom N. Kandel and Dong Liang.  The long wave fluid flows on inclined porous media with nonlinear Forchheimer’s law.  AIP Advances, 9(9):095302, 2019.
  * [16] Alena Kruisová, Martin Ševčík, Hanuš Seiner, Petr Sedlák, Benito Román-Manso, Pilar Miranzo, Manuel Belmonte, and Michal Landa.  Ultrasonic bandgaps in 3d-printed periodic ceramic microlattices.  Ultrasonics, 82:91–100, Jan 2018.
  * [17] G. Laschet.  Forchheimer law derived by homogenization of gas flow in turbomachines.  Jour. of Comput. and Appl Math., 215:467–476, 2008.
  * [18] R. B. Lehoucq, D. C. Sorensen, and C. Yang.  ARPACK Users’ Guide.  Software, Environments and Tools. Society for Industrial and Applied Mathematics, Jan 1998.
  * [19] N. Masmoudi.  Homogenization of the compressible Navier–Stokes equations in a porous medium.  ESAIM: Control, Optimisation and Calculus of Variations, 8:885–906, 2002.
  * [20] A. Mikelic.  Homogenization of nonstationary Navier-Stokes equations in a domain with a grained boundary.  Annali di Matematica Pura ed Applicata, 158:167–179, 1991.
  * [21] A. Mikelic and L. Paoli.  Homogenization of the inviscid incompressible fluid flow through a 2D porous medium.  Proceedings of the American Mathematical Society Volume 127, Num, 127(7):2019–2028, 1999.
  * [22] E. Miroshnikova.  Some new results in homogenization of flow in porous media with mixed boundary condition.  PhD thesis, Lulea University of Technology, 2016.  Graphic Production.
  * [23] A.N. Norris.  On the viscodynamic operator in Biot’s equations of poroelasticity.  J. Wave-Material Interaction, 1:365–380, 1986.
  * [24] M. Peszyńska and A. Trykozko.  Forchheimer law in computational and experimental studies of flow through porous media at porescale and mesoscale.  In Current Advances in Nonlinear Analysis and Related Topics, volume 32 of GAKUTO International Series, pages 463–482, 2010.
  * [25] D. Poliševski.  Homogenization of Navier-Stokes model: the dependence upon parameters.  Z. angew. Math. Phys., 40:387–394, 1989.
  * [26] Raghu Raghavan.  Theory for acoustic streaming in soft porous matter and its applications to ultrasound-enhanced convective delivery.  Journal of Therapeutic Ultrasound, 6:6, 2018.
  * [27] E. Sanchez-Palencia.  Non-homogeneous media and vibration theory.  Number 127 in Lecture Notes in Physics. Springer, Berlin, 1980.
  * [28] G. W. Stewart.  A krylov–schur algorithm for large eigenproblems.  SIAM Journal on Matrix Analysis and Applications, 23(3):601–614, Jan 2002.
  * [29] F. Tisseur and K. Meerbergen.  The quadratic eigenvalue problem.  SIAM Review, 43:235–286, 2001.
  * [30] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, CJ Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake Vand erPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1. 0 Contributors.  SciPy 1.0: fundamental algorithms for scientific computing in Python.  Nature Methods, page 1–12, Feb 2020.
  * [31] Yan-Feng Wang, Jun-Wei Liang, A-Li Chen, Yue-Sheng Wang, and Vincent Laude.  Wave propagation in one-dimensional fluid-saturated porous metamaterials.  Phys. Rev. B, 99:134304, Apr 2019.
  * [32] R. Zaki.  Homogenization of a Stokes problem in a porous medium by the periodic unfolding method.  Asymptotic Analysis, 79(3–4):229–250, 2012.
  * [33] Y. Zhengan and Z. Hongxing.  Homogenization of a stationary Navier–Stokes flow in porous medium with thin film.  Acta Mathematica Scientia, 28B(4):963–974, 2008.