lemmadefinition theoremdefinition corollarydefinition propositiondefinition
assumptiondefinition remarkdefinition exampledefinition conjecturedefinition
lemma theorem corollary proposition remark conjecture assumption

# Variational methods for fluid-structure interaction and porous media

B. Benešová, M. Kampschulte, S. Schwarzacher  Department of Mathematics and
Physics, Charles University Prague<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

In this work we consider a poroelastic flexible material that may deform
largely which is situated in an incompressible fluid driven by the Navier-
Stokes equations in two or three space dimensions. By a variational approach
we show existence of weak solutions for a class of such coupled systems. We
consider the unsteady case, this means that the PDE for the poroelastic solid
involves the Frechet derivative of a non-convex functional as well as (second
order in time) inertia terms.

###### Key words and phrases:

Mathematics for continuum mechanics, Fluid poroelastic structure interactions,
Minimizing movements, Navier-Stokes equations, Elastic solids, Hyperbolic
evolutions, Coupled systems of PDEs

##  1\. Introduction

_Poromechanics_ or _mechanics of porous media_ has been a lively area of
research in engineering and continuum mechanics (see e.g. the monographs [18,
16]) as the possible applications of the proposed theories are numerous,
ranging from the mechanics of tissue and biological materials in general (see
e.g. [27]) to the recently important breathing through masks. In many of these
cases, the poroelastic matrix is not isolated but instead immersed into a
fluid which may flow through the elastic structure. Naturally, then, the
porous medium and the free fluid are interacting with each other. Indeed on
the mutual interface they both exert pressure and stresses to each other and
on the top of that the fluid may enter the pores of the structure and flow
through it. Such a setting is referred to as _fluid-pororelastic structure
interaction (FPSI)_ and will be the subject of the present work. The main
result of this paper is the well-posedness for a class of FPSI problems
allowing large deformations and involving incompressible fluids driven by the
Navier-Stokes equations in two or three space dimensions.

In order to predict the response of the FPSI-system, we need to devise a model
for the free fluid, the poroelastic structure and prescribe suitable interface
conditions. Each of the tasks can be approached on several scales, ranging
from the microscopic to the macroscopic one. For the former the skeleton is
modeled in detail and a (yet another) fluid-structure interaction problem
between the skeleton and the circumflowing fluid is set up. More common in the
engineering literature, however, is a more _macroscopic approach_ where the
porous medium is assumed to be saturated by the fluid, i.e. they form a kind
of “mixture” in each material point and an _ad hoc_ form of the stress tensor
is prescribed for the porous medium [12]. Such an approach covers the early
and famous models for the porous medium due to Darcy [17], Brinkman [10] or
Biot [7] and will be also the one we follow in this contribution. Let us also
note that a justification of the proposed models by homogenization is highly
relevant indeed and some results concerning the aforementioned standard models
for porous media have already been obtained, see e.g. [1, 15, 24].

Following the macroscopic approach, we will assume that the porous medium is
fully saturated by the fluid and can be completely described by the fluid
velocity $v$ and the deformation of the solid $\eta$ (see Section 2 for a more
detailed modeling description). We will take a _variational point of view_ for
modeling as well as for the analysis by defining the constitutive behavior of
the mechanical system by prescribing its _stored energy_ as well as
_dissipation functional_ (dissipation pseudo-potential). Such an approach has
been advocated by many authors (e.g. [22] for solids or [36, 37, 35] for
fluids) as it may simplify the modeling. On the top of that, we argue in
Section 3 that the variational approach is essential also from the point of
view of _mathematical analysis_ if a _large-strain, non-linear_ model is to be
used for the involved solid. Indeed, in such a case, the Helmholtz free energy
cannot be chosen as a convex functional for physical restrictions and thus the
resulting differential operator is not monotone.

In the case of FPSI, the variational approach may also facilitate the
derivation of boundary conditions at the interface of the poroelastic medium
itself and the fluid into which it is immersed, see e.g. [20]. Prescribing the
above mentioned energy and dissipation for the whole system then automatically
yields a balance of stresses at this boundary and, additionally, imposes a
“weak continuity” condition for the fluid velocity over the porous media
interface, as we ask for its derivatives to be globally integrable.
Nonetheless, generalizations are, of course, possible. For this see in
particular the notes at the end of Section 2.

While the engineering literature on models in poroelasticity and FPSI seems to
be rich, the mathematical literature studying well-posedness of the FPSI is,
to the authors’ knowledge, rather scarce. Interaction of the linear or non-
linear poroelastic medium with Stokes flow has been studied in [31, 2, 8].
Moreover, the coupling of the Stokes and Darcy flow has been studied (see e.g.
[34] for an overview) but here, intrinsically, the region occupied by the
porous medium is fixed a-priori. Similarly, the interaction of the Biot-
poroelastic medium with a Navier-Stokes flow was studied in [11], but again
with the interface wall between the porous solid and fluid fixed. Thus, well-
posedness of a FPSI problem featuring a fully non-linear poroelastic solid
that may undergo large deformations coupled to a Navier-Stokes fluid has not
been addressed in literature so far. In this contribution, we work exactly in
this setting and we establish _existence of weak solutions_. In order to do
so, it is crucial to approach the modeling as well as the subsequent analysis
in an _energetic (or variational)_ way to access methods from the calculus of
variations. We provide an overview of the analytical details of the approach
in Section 3 but would like to emphasize here that working the proposed
framework not only allows to prove the result presented here but also paves
the road for studying further FPSI models. Indeed, models including dependence
of the stresses on e.g. the porosity, the porous pressure of the fluid or more
general coupling conditions are accessible to the analysis, though a careful
investigation will be needed.

This work is built up as follows: In Section 2 we present the model considered
in this work and formulate the main results. In Section 3 we overview the
ingredients of the variational approach and in Section 4 we present the proof
of the main result, Theorem 2.

##  2\. Modeling and main results

To set up the porous media problem, we consider a regular enough111The
regularity of the boundary of $\Omega$ as well as $Q$ only fully comes into
play when discussing contact of the solid with $\partial\Omega$. This is
discussed in a bit more detail in [6], and will be thoroughly discussed in a
forthcoming paper involving the second author. For the purpose of the paper at
hand, we ignore this aspect completely and just require enough regularity to
meaningfully talk about boundary values. spatial domain
$\Omega\subset\mathbb{R}^{n}$ in which both the porous medium an the
circumlying fluid are contained.

$Q$$\Omega$$\Omega$$\eta$$v$$\partial_{t}\eta\circ\eta^{-1}$$\partial_{t}\eta$
Figure 1. A short illustration of the geometry involved in the model.

As already announced, we assume that the poroelastic medium is a perfectly
homogenized mixture of the solid and the fluid in each material point while
outside the porous medium the fluid is found in its “pure” state. Thus, the
fluid is allowed to flow through the whole of $\Omega$, including the part
already occupied by the solid. The state of the solid will be described by its
deformation, a map $\eta:Q\to\Omega$; i.e. from some regular enough reference
domain $Q\subset\mathbb{R}^{n}$ to the container $\Omega$. We choose to fix
“boundary values” $\eta|_{P}=\gamma$ for some subset $P\subset\overline{Q}$,
but technically this is optional. The deformation $\eta$ additionally fixes
the geometry of the porous medium $\eta(Q)$.

The fluid, on the other hand, is described by its velocity
$v:\Omega\to\mathbb{R}^{n}$ and pressure $p:\Omega\to\mathbb{R}^{n}$ with
$\operatorname{div}v=0$. In order to simplify the discussion, we will assume
no-slip conditions at the boundary of $\Omega$, i.e. $v|_{\partial\Omega}=0$.
As the fluid is assumed incompressible, the pressure will later “disappear”
from the weak formulation thanks to the choice of appropriate test functions.

Having defined all the quantities involved, we can now turn back to modeling
by prescribing the stored energy. We assume an additive decomposition of the
energy for the fluid and the solid. However, as the fluid is assumed to be
purely Newtonian and incompressible, no energy can be stored in the fluid.
Thus, the stored energy consists just of the solid one denoted as $E(\eta)$.
The overall energy is then given by the sum

$\displaystyle\int_{\Omega}\frac{\rho_{f}}{2}\left|{v}\right|^{2}dy+\int_{Q}\frac{\rho_{s}}{2}\left|{\partial_{t}\eta}\right|^{2}dx+E(\eta).$

Here $\rho_{f}$ and $\rho_{s}$ are given mass densities of the fluid and solid
respectively.

To complete the modeling we need to prescribe a form of $E(\eta)$. For
simplicity, we stick to the the Saint Venant-Kirchhoff energy enriched with a
term introducing resistance of the solids to infinite compression and a
higher-order gradient regularisation:

(2.1) $\displaystyle E(\eta)$
$\displaystyle:=\begin{cases}\int_{Q}\frac{1}{8}|\nabla\eta^{T}\nabla\eta-I|^{2}_{\mathcal{C}}+\frac{1}{(\det\nabla\eta)^{a}}+\frac{1}{q}\left|{\nabla^{2}\eta}\right|^{q}dx&\text{if
$\det\nabla\eta>0$ a.e. in $Q$ }\\\ +\infty&\text{otherwise}\end{cases}$

where we used the notation
$|\nabla\eta^{T}\nabla\eta-I|^{2}_{\mathcal{C}}:=\big{(}\mathcal{C}(\nabla\eta^{T}\nabla\eta-I)\big{)}\cdot\big{(}\nabla\eta^{T}\nabla\eta-I\big{)}$,
with $\mathcal{C}$ being a fourth order positive definite tensor of elastic
constants, $q>n$ and $a>\frac{qn}{q-n}$.

######  Remark  (Non-simple materials).

Let us stress that the proposed energy functional for the solid depends not
only on the _deformation gradient_ , which is the standard setting, but it
also includes a term depending on the _second gradient of the deformation_.
Thus, the proposed model belongs to the class of so-called _non-simple_
materials that have been proposed by Toupin [33] and later further developed
in other works, e.g. [28, 32, 21]. In the mathematical literature concerned
with the large-deformation regime this concept is often employed to gain
additional, necessary regularity. Also the present work relies on it. Indeed,
at the heart of our variational approach (see also Section 3) is that we will
be able to solve a time-discrete variant of the set of PDEs in (2.5) by
solving a minimization problem and then relying on the fact that this
minimizer solves the associate Euler-Langrage equation. However, as our energy
functional takes infinite values, the relation between the functional and the
Euler-Lagrange equation is not straighforward. Indeed, the minimizer may lie
almost at the boundary of the domain of $E(\eta)$ so that taking variations at
the minimizer would not be possible. This is indeed a well-known problem in
mathematical elasticity [4] that has not found a satisfactory solution to
date. Let us just mention that even in the one-dimensional setting explicit
examples have been found in [5] that illustrate this exact phenomenon: a
minimizer of the functional exists but it does not satisfy the Euler-Lagrange
equation.

To the authors’ knowledge, as of today, the only setting where it is known
that this phenomenon does not occur is the one of _non-simple materials_.
Indeed, Healey and Krömer have found in [23] that if the indices $q$ and $a$
in $E(\eta)$ are chosen appropriately, then every deformation of a finite
energy has a uniform lower bound for its Jacobian. Thus, the minimizer is well
separated from the boundary of the domain in $E(\eta)$ and variations at the
minimizer can be performed. ∎

Let us just remark that the second term in $E(\eta)$ assures that the stored
energy blows up whenever $\det({\nabla\eta})\to 0$ and as such penalizes
compression (in fact we get a uniform energy dependent lower bound on the
determinant, see [23]). This is a crucial requirement non-linear elasticity
(see e.g. [13]). We also stress that such a requirement alone rules out
convexity of the stored energy, albeit when $E(\eta)$ is of the form of (2.1)
also the first term is not convex.

Much more challenging for the analysis than the non-convex terms in the energy
is the fact, that the admissible set of deformations has to be chosen
carefully; in particular it is bound to be a non-convex set. This is related
to the physical expectation that any admissible deformation will be injective.
Indeed, while it might be possible to make sense of the equations of the
problem even in the case that two different parts of the solid occupy the same
Eulerian space, physically such a situation is of course nonsensical. Now
observe that locally the blow-up of the determinant already guarantees
injectivity of any deformation of finite energy. But in addition to the local
injectivity we further have to take into account the important non-convex
restriction of global injectivity. What thus turns out to be a good choice for
the admissible set of $\eta$ is the following set that is in coherence with
the celebrated _Ciarlet-Nečas condition_ proposed in [14]:

(2.2) $\displaystyle\mathcal{E}:=\left\\{\eta\in
W^{2,q}(Q;\Omega)\,:\,E(\eta)<\infty,\,\left|{\eta(Q)}\right|=\int_{Q}\det\nabla\eta\,dx\right\\}.$

Here, the finite energy guarantees local injectivity and together with the
last condition we have that any $C^{1}$-local homeomorphism is globally
injective except for possible touching at the boundary. Hence in the following
we will construct deformations in the set (2.2). In fact the set $\mathcal{E}$
only includes deformations for which the dissipation function $R$ of the solid
deformation, that we introduce next, is also well defined.

As for the dissipation function, we also assume an additive splitting of the
fluid and solid dissipation, respectively. In addition to that we introduce a
phenomenological coupling term, a _drag term_ , that introduces additional
dissipation if the velocities of the fluid and solid are not equal. For the
purposes of this paper we stick to the simplest possible setting when the drag
force is proportional to the difference of the solid and fluid velocity (see
[29]); this is also the setting within which, upon further simplifications,
the law of Darcy or the equations of Brinkman are derived [29]. However since
one of these velocities is given in Lagrangian and the other in Eulerian
representation, we need to transform one of them into the proper frame, the
easier of the two options being the Lagrangian. Thus, we obtain for the the
dissipation functional:

$\displaystyle\mathcal{D}=\underbrace{\frac{\nu}{2}\left\|{\varepsilon}v\right\|_{{\Omega}}^{2}}_{\text{fluid
dissipation}}+\underbrace{R(\eta,\partial_{t}\eta)}_{\text{solid
disspation}}+\underbrace{A(\eta,\partial_{t}\eta-v\circ\eta)}_{\text{drag
potential}}$

where

$A(\eta,\phi):=\int_{Q}\frac{1}{2}\phi^{T}\cdot a(\nabla\eta)\cdot\phi\,dx$

and $a:\\{M\in\mathbb{R}^{n\times n}:\det M>0\\}\to\mathbb{R}^{n\times n}$ is
a smooth function mapping into the set of positive semi-definite symmetric
matrices. Note that the drag potential is quadratic, so that the drag force
depends linearly, through a factor depending on $\nabla\eta$, on the
difference of the two velocities. This quite general dependence on
$\nabla\eta$ is necessary to model anisotropy in the correct frame independent
way.222If we assume the simpler case of an isotropic solid with regards to
drag, we can instead use
$A(\eta,b)=\frac{a_{0}}{2}\left\|b\right\|_{{Q}}^{2}$, but doing so does not
fundamentally change any of the later computations. Note also that for
physical reasons the drag does only depend on the relative and not on the
absolute velocities. Finally, we need to prescribe the solid dissipation, for
which we choose the following form which leads to an viscoelastic solid of
Kelvin-Voigt type:

(2.3) $\displaystyle R(\eta,\partial_{t}\eta)$
$\displaystyle:=\int_{Q}|(\nabla\partial_{t}\eta)^{T}\nabla\eta+(\nabla\eta)^{T}(\nabla\partial_{t}\eta)|^{2}dx=\int_{Q}|\partial_{t}(\nabla\eta^{T}\nabla\eta)|^{2}dx.$

Here, let us notice that the dissipation of the solid depends on the state.
This is indeed necessary to assure frame indifference [3].

With all quantities in hand, we can now write down the (formal) energy
equality

$\displaystyle\underbrace{E(\eta(T))+\frac{\rho_{s}}{2}\left\|\partial_{t}\eta(T)\right\|_{{Q}}^{2}+\frac{\rho_{f}}{2}\left\|v(T)\right\|_{{\Omega}}^{2}}_{\text{energy
at final time
$T$}}+\underbrace{\int_{0}^{T}\left[\nu\left\|{\varepsilon}v\right\|_{{\Omega}}^{2}+2R(\eta,\partial_{t}\eta)+2A(\eta,\partial_{t}\eta-v\circ\eta)\right]dt}_{\text{dissipated
energy}}$
$\displaystyle=\underbrace{E(\eta_{0})+\frac{\rho_{s}}{2}\left\|b\right\|_{{Q}}^{2}+\frac{\rho_{f}}{2}\left\|v_{0}\right\|_{{\Omega}}^{2}}_{\text{energy
at the initial time $0$}}+\underbrace{\int_{0}^{T}\left\langle
f,v\right\rangle_{\Omega}dt.}_{\text{power of the external forces}}$

Here $f$ is an external force acting on the fluid. Forces acting on the solid
can be considered in a similar fashion. We use the subscript “0” to denote the
initial data for the solid deformation and the fluid velocity given at $t=0$,
while $b$ is used for the initial solid velocity.

Once the stored energy and the dissipation function have been specified, we
embark on giving a weak and a strong formulation of the considered problem.
This formulation needs to be closed using boundary conditions on the interface
between the pure fluid and the porous media. We do so in the least invasive
way, by requiring weak continuity of the fluid (in the sense that $v\in
W^{1,2}(\Omega;\mathbb{R}^{n})$) and stress free boundary condition for the
solid deformation.

Now, we start with the weak formulation and fix the space of test functions,
which consists of the pairs $(\phi,\xi)\in W^{2,q}(Q;\mathbb{R}^{n})\times
W^{1,2}_{0}(\Omega;\mathbb{R}^{n})$ with $\phi|_{P}=0$. Taking the derivatives
of stored energy and dissipation with respect to position and velocity
respectively, in directions $\phi$ and $(\phi,\xi)$ gives us the corresponding
forces and we end up with:

######  Definition 2.1.

We call the pair $(\eta,v)\in L^{\infty}(0,T;\mathcal{E})\cap
W^{1,2}(0,T;W^{1,2}(Q))\times L^{2}(0,T;W^{1,2}_{0}(\Omega;\mathbb{R}^{n}))$
with $\eta|_{P}=0$ and $\eta(0)=\eta_{0}$ satisfying

(2.4) $\displaystyle\left\langle b,\phi(0)\right\rangle_{Q}+\left\langle
v_{0},\xi(0)\right\rangle_{\Omega}$
$\displaystyle=\int_{0}^{T}\rho_{s}\left\langle\partial_{t}\eta,\partial_{t}\phi\right\rangle_{Q}+\left\langle
DE(\eta),\phi\right\rangle+\left\langle
D_{2}R(\eta,\partial_{t}\eta),\phi\right\rangle+\nu\left\langle{\varepsilon}v,{\varepsilon}\xi\right\rangle_{\Omega}\
$ $\displaystyle\quad+\left\langle
D_{2}A(\eta,\partial_{t}\eta-v\circ\eta),\phi-\xi\circ\eta\right\rangle+\rho_{f}\left\langle
v,\partial_{t}\xi-v\cdot\nabla\xi\right\rangle_{\Omega}-\left\langle
f,\xi\right\rangle_{\Omega}dt,$

for all $(\phi,\xi)\in C^{\infty}([0,T]\times Q)\times
C^{\infty}(0,T;C^{\infty}_{0}(\Omega))$, with $\phi(T)=0$, $\phi|_{P}=0$,
$\operatorname{div}\xi=0$ and $\xi(T)=0$ a weak solution of the FPSI-problem
given by (2.1) and (2.3) with respective initial and boudnary values.

The corresponding strong formulation reads

(2.5) $\displaystyle\left\\{\begin{aligned}
\rho_{s}\partial_{t}^{2}\eta&=\operatorname{div}\sigma&\text{ in }&[0,T]\times
Q\\\
\operatorname{div}\sigma&=DE(\eta)+D_{2}R(\eta,\partial_{t}\eta)+a(\nabla\eta)\cdot(\partial_{t}\eta-v\circ\eta)&\text{
in }&[0,T]\times Q\\\ \operatorname{div}v&=0&\text{ in }&[0,T]\times\Omega\\\
\rho_{f}(\partial_{t}v+v\cdot\nabla v)&=\nu\Delta v+F-f-\nabla p&\text{ in
}&[0,T]\times\Omega\\\
F&=(\det\nabla(\eta^{-1}))a(\nabla\eta\circ\eta^{-1})\cdot(v-\partial_{t}\eta\circ\eta^{-1})&\text{
in }&[0,T]\times\eta(Q)\\\ F&=0&\text{ in
}&[0,T]\times\Omega\setminus\eta(Q)\\\ \sigma\cdot n&=0&\text{ in
}&[0,T]\times\partial Q\setminus P\\\ \eta&=\gamma&\text{ in
}&[0,T]\times\partial P\\\ v&=0&\text{ in }&[0,T]\times\partial\Omega\\\
v(0)&=v_{0}&\text{ in }&\Omega\\\ \eta(0)=\eta_{0}&\text{ and
}\partial_{t}\eta(0)=b&\text{ in }&Q,\end{aligned}\right.$

where $\sigma$ is the Piola-Kirchhoff stress tensor derived from $E$ and $F$
is the Eulerian drag from the solid acting on the fluid. In deriving the
strong formulation from the energy contribution, the drag term is used twice
(it appears in the PDE for the velocity $v$ as well as the deformation $\eta$)
but deriving both constituents from the same term guarantees that these two
forces are indeed opposite and equal. Notice also that the boundary value of
$\sigma$ is of Neumann type and due to a balance of stresses on the porous
medium - fluid interface; in the weak formulation this is already encoded
implicitly.

The goal of this paper is to prove the following theorem:

######  Theorem  (Existence for the porous media problem).

Given $T>0$, $f\in L^{2}([0,T]\times\Omega;\mathbb{R}^{n})$,
$\eta_{0}\in\mathcal{E}$, $b\in L^{2}(Q;\mathbb{R}^{n})$, $v_{0}\in
L^{2}(\Omega;\mathbb{R}^{n})$ there exists time $T^{*}\in(0,T]$ and a weak
solution (in the sense of (2.4)) to (2.5). More precisely we show the
existence of a pair

$(\eta,v)\in C^{0}([0,T^{*}];W^{2,q}(Q;\mathbb{R}^{n}))\cap
W^{1,2}([0,T^{*}];W^{1,2}(Q;\mathbb{R}^{n}))\times
W^{1,2}([0,T^{*}]\times\Omega;\mathbb{R}^{n})$

with initial data $\eta(0)=\eta_{0}$, $\partial_{t}\eta_{0}=b$, $v(0)=v_{0}$
(in the weak sense), satisfying (2.4). Additionally this solution satisfies
the physical energy inequality. Here $T^{*}$ either equals $T$ or is the the
first time of (self-) contact, where the free part of $\partial\eta(T^{*},Q)$
either touches itself or $\partial\Omega$.

Theorem 2 is proved in Section 4. However, as the method of proof has the
potential to be extended to more general models (see below) we first present,
for a better understanding of the techniques, a respective overview in Section
3.

### Limitations and generalizations

The model considered in this paper for which we prove existence of weak
solutions in Theorem 2 includes many important features that have not been
considered in mathematical FPSI yet: namely, the model is _fully dynamic
including inertia_ and the response of the solid is _fully nonlinear allowing
for large structural deformations_. On the other hand, we have made several
simplifying assumptions including the following ones:

  1. (1)

Non-saturated/non-zero volume porous media: In the model at hand, the volume
fraction of the fluid inside the porous medium is always fixed to one. This
greatly simplifies the discussion as it in particular implies
$\operatorname{div}v=0$ everywhere, even over the interface. It is however
possible to generalize this condition and assume divergence free flow only on
a part of the domain. In particular in [6] we already dealt with a global
velocity field that is only divergence free outside the solid.

  2. (2)

Dependence of the solid stored energy on fluid variables: In the model we
assume that the stored energy of the solid does not depend on the properties
of the saturating fluid. Generalizations in this directions are known in
modeling and include e.g. Biot’s model where the stored energy depends on the
pressure of the fluid in pores or models that depend on the volume fraction of
the fluid and solid in each point.

  3. (3)

Compressibility of the fluid: The fluid is assumed incompressible and
Newtonian in the present model. Nonetheless, in view of possible applications
like the breathing through masks, incorporating compressibility seems an
important task. This will call for changes in the analysis, but possibly also
modelling, in particular if the pressure additionally enters the solid stored
energy like in Biot’s model. (See [9] for an example of the variational
strategy in the context of fluid-structure interaction with a compressible
fluid).

  4. (4)

Drag force: We assume the most basic form of the drag force that depends
linearly on the difference between the fluid and solid velocity. It is also
just the drag force where the coupling between the fluid and solid velocities
happens in the porous medium model. One could think of a drag force that
depends also on other variables included in the model (like volume fractions)
or other coupling possibilities.

  5. (5)

Interface conditions and boundary layers: The interface condition at the
transition between the porous medium and the circumlying fluid currently
consists of continuity of the fluid velocity and a stress free boundary
conditions for the solid. Generalizations are thinkable for instance an
additional boundary layer described by an additional dissipation functional.

  6. (6)

Geometry of the porous medium: In many cases the porous medium can be well
described as a lower-dimensional structure, thus an appropriate membrane,
shell or plate model is applicable. Applications include cell walls or the
aforementioned breathing through masks. The variational approach presented
here is well suited to handle also these situations; for the classical FSI
this situation is currently investigated [25] so that generalizations to the
FPSI setting seem feasible.

It seems plausible that the discussed aspects (or other) could be incorporated
into the presented model in a future work. Let us however remark that some
limitations of the given model seem not to be removable with the current state
of art. Most significantly for the solid model. In particular it is essential
to include regularizations featuring the second gradient (see Remark 2).

##  3\. Variational strategy - Overview of the proof of Theorem 2

We are interested in proving existence of weak solutions to the system (2.5)
in the sense of Definition 2.1. As is standard, we will first prove existence
of approximations of (2.5) and then pass to the limit. To enable the limit
passage a crucial prerequisite is the availability of suitable a-priori
estimates, which are usually based on an appropriate energy (in)equality.
Indeed, multiplying formally the balances of the respective momenta in (2.5)
by the velocity and integrating over $[0,T]$ leads to

$\displaystyle\underbrace{E(\eta(T))+\frac{\rho_{s}}{2}\left\|\partial_{t}\eta(T)\right\|_{{Q}}^{2}+\frac{\rho_{f}}{2}\left\|v(T)\right\|_{{\Omega}}^{2}}_{\text{energy
at final time
$T$}}+\underbrace{\int_{0}^{T}\left[\nu\left\|{\varepsilon}v\right\|_{{\Omega}}^{2}+2R(\eta,\partial_{t}\eta)+2A(\eta,\partial_{t}\eta-v\circ\eta)\right]dt}_{\text{dissipated
energy}}$
$\displaystyle=\underbrace{E(\eta_{0})+\frac{\rho_{s}}{2}\left\|\partial_{t}\eta(0)\right\|_{{Q}}^{2}+\frac{\rho_{f}}{2}\left\|v_{0}\right\|_{{\Omega}}^{2}}_{\text{energy
at the initial time $0$}}+\underbrace{\int_{0}^{T}\left\langle
f,v\right\rangle_{\Omega}dt}_{\text{power of the external forces}},$

which, on a more abstract level, takes the compact form

$\displaystyle
E_{\text{st}}(T)+E_{\text{kin}}(T)+\int_{0}^{T}W_{\text{diss}}(t)dt=E_{\text{st}}(0)+E_{\text{kin}}(0)+\int_{0}^{T}W_{\text{ext}}(t)dt$

where we distinguished the the following four contributions: The stored energy
$E_{\text{st}}$, the kinetic energy $E_{\text{kin}}$, energy lost through
dissipation $W_{\text{diss}}$ and work done by external forces
$W_{\text{ext}}$. The kinetic energy and the external forces will generally
each always have a similar form, independent of the considered material. On
the other hand, the other two carry the constitutive information on material
modelling. In this work, we stick to the form (2.1) and (2.3) of these
contributions, but the method of proof is to a large extent independent of
those, so we shall not use the specific form within this section.

We will prove Theorem 2 by _time-discretization_ ; this is particularly suited
as the geometry between the porous medium and the circumlying fluid changes
during the evolution. Thus, even if we do not deal with the situation here, a
time-discretization method has the potential to handle complex transition
conditions between the porous medium and the fluid. More importantly, however,
a time-discretization allows us to set up a _variational problem_ to prove
_existence of approximate solutions_. This is of great importance as
standardly used fixed point methods do not function well if the included
operators are not monotone (in some sense) and/or the set of admissible
functions is not convex.

Nonetheless, the selected time-discretization method needs to be chosen
carefully as the balance of momentum for the solid porous medium (i.e. first
equation in (2.5)) is _hyperbolic featuring a non-monotone differential
operator stemming from the non-convex stored energy_. Let us explain the
difficulty on a simple toy model involving a single unit mass particle with
position $x(t)\in\mathbb{R}^{n}$ and a _non-convex_ potential energy
$E(x(t))$. Hence, we seek the solution to the following hyperbolic ODE:

$\partial_{t}^{2}x=-\nabla E(x)$

with initial data $x(0)=x_{0}$ and $\partial_{t}x(0)=x^{*}$. The naive ansatz
is to consider a time-discretization with step-size $\tau$ and simply set up a
backward Euler discretization as follows:

(3.1) $\displaystyle 0=\nabla
E(x_{k+1})+\frac{\tfrac{x_{k+1}-x_{k}}{\tau}-\tfrac{x_{k}-x_{k-1}}{\tau}}{\tau}.$

Assuming that existence of solutions to (3.1) can be proved, one would like to
derive an energy inequality from (3.1) to obtain a-priori estimates. Mimicking
the continuous case, the most straighforward way to do so seems to test (3.1)
with the discretized time derivative. This yields

$0=\left\|\tfrac{x_{k+1}-x_{k}}{\tau}\right\|^{2}-\left\langle\tfrac{x_{k}-x_{k-1}}{\tau},\tfrac{x_{k+1}-x_{k}}{\tau}\right\rangle+\left\langle\nabla
E(x_{k+1}),\tfrac{x_{k+1}-x_{k}}{\tau}\right\rangle,$

where the second term can be estimated using Young’s inequality. Thus, we need
to concentrate on the last term; in the time-continuous case we would apply a
chain rule here. However in the discrete setting a chain rule is generally not
available, only if $E$ would be convex the following inequality

$\left\langle\nabla E(x_{k+1}),\tfrac{x_{k+1}-x_{k}}{\tau}\right\rangle\geq
E(x_{k+1})-E(x_{k}),$

that is sometimes called the _discrete chain rule_ could be applied. In the
_non-convex setting_ this inequality is false in general and this path to
a-priori estimates seems closed.333In some cases, e.g. in the case of
$\lambda$-convexity, it is still possible to derive a variant of the discrete
chain rule with some additional error terms on the right hand side. This can
lead the way to an a-priori estimate and a resulting existence proof but is
much more problem specific than the general method we present here.

The solution to this quandry has been found in [6] by noting that rather a
_two-scale_ approximation is needed in (3.1). Thus we consider two parameters
$\tau,h$ and write

(3.2) $0=\nabla
E(x_{k+1})+\frac{\tfrac{x_{k+1}-x_{k}}{\tau}-\tfrac{x_{k}-x_{k-1}}{\tau}}{h}.$

If $h$ is fixed and $t<h$, (3.2) can be viewed as a $\tau$-time discretization
of the following gradient flow problem under forcing

(3.3) $\displaystyle\nabla E(x(t))=-\frac{\partial_{t}x(t)-x^{*}}{h},\quad
x(0)=x_{0},$

where $x^{*}=\partial_{t}x(0)$. Indeed, if we could pass to the limit $\tau\to
0$ in (3.2) we construct a solution $x^{h}$ to (3.3) on $[0,h]$. Now using the
known values of $\partial_{t}x(t-h)$ in place of $x^{*}$ for the next interval
of length $h$ we can iteratively prolong this process to get what we will call
a _time-delayed solution_ , satisfying

(3.4) $\displaystyle\nabla
E(x(t))=-\frac{\partial_{t}x(t)-\partial_{t}x(t-h)}{h}$

for $t\in[0,T]$. We can test this time-delayed solution with
$\partial_{t}x(t)$ and find the hyperbolic a-priori estimate

$\displaystyle E(x(b))-E(x(a))$
$\displaystyle=-\int_{a}^{b}\left\langle\frac{\partial_{t}x(t)-\partial_{t}x(t-h)}{h},\partial_{t}x(t)\right\rangle\,dt$
$\displaystyle\quad\leq-\frac{1}{2}\fint_{b-h}^{b}\left|{\partial_{t}x(t)}\right|^{2}dt+\frac{1}{2}\fint_{a-h}^{a}\left|{\partial_{t}x(t)}\right|^{2}dt,$

whenever the solution was constructed over $[a-h,b]$.444Actually in the
construction procedure the hyperbolic a-priori estimate already needs to be
used for each interval of length $h$ to guarantee the admissibility of the
initial/right hand data of the next.

The above gives us a good estimate on $E(x(t))$ and an averaged time-
derivative, independent of $h$ which allows for sending $h\to 0$ in (3.4).
Thus this idea paves the path to Theorem 2.

On the level of the $\tau$-discretization we deal with the gradient flow
problem (3.4). For these variational methods are well established. We will
explain this in the following. Let us consider the rescaled problem with $h=1$
so that we obtain

$\nabla E(x(t))=-\partial_{t}x(t)-f,\quad x(0)=x_{0},$

with some given time dependent force $f$, which includes the previously
constructed $\partial_{t}x(t-h)=\partial_{t}x(t-1)$. We may discretize this by
means of the backward Euler method to get

(3.5) $\displaystyle 0=\nabla E(x_{k+1})+\frac{x_{k+1}-x_{k}}{\tau}+f_{k}.$

where $f_{k}$ is a corresponding time discretization of $f$. Now one realizes
that (3.5) is actually the Euler-Lagrange equation to the minimization problem

$\mathcal{F}_{k}(x)\longrightarrow\text{min},$

where $\mathcal{F}_{k}:x\mapsto
E(x)+\frac{1}{2\tau}\left|{x-x_{k}}\right|^{2}+(x-x_{k})\cdot f_{k}$. Thus,
one constructs $x_{k+1}$ using the minimization problem. Now, since one
specifically has a minimizer, instead of using the equation one can compare
the values of $F_{k}(x)$ at $x_{k+1}$ and $x_{k}$ in order to derive a-priori
estimates. This leads to

$\displaystyle
E(x_{k+1})+\frac{\tau}{2}\left|{\frac{x_{k+1}-x_{k}}{\tau}}\right|^{2}+\tau\frac{x_{k+1}-x_{k}}{\tau}f=\mathcal{F}_{k}(x_{k+1})\leq\mathcal{F}_{k}(x_{k})=E(x_{k}),$

and summing up yields the energy inequality (up to a factor on the
dissipation) and, in turn, a-priori estimates. This method, when one resorts
to minimization, has been known as the _minimizing movements approximation_ of
gradient flows and is due to De Giorgi [19]. Since then it has been widely
used in mathematical solid mechanics since it allows to cope with the non-
convexities of $E$ (see e.g. [26]) and is well suitable also in our case.

######  Remark  (Minimization in the hyperbolic case).

A more direct attempt at a proof of the main theorem would be to generalize
the minimizing movements by minimizing the functional

(3.6)
$\displaystyle\mathcal{F}_{k}(x):=E(x)+\frac{1}{2}\left|{\tfrac{x-x_{k}}{\tau}-\tfrac{x_{k}-x_{k-1}}{\tau}}\right|^{2},$

which is easily checked to have (3.1) as Euler-Lagrange equation. But then the
same arguments apply to establishing estimates via the equation and comparing
the values of $F_{k}(x)$ at $x_{k+1}$ and $x_{k}$ instead yields

$\displaystyle
E(x_{k+1})+\frac{\tau^{2}}{2}\left|{\frac{\tfrac{x_{k+1}-x_{k}}{\tau}-\tfrac{x_{k}-x_{k-1}}{\tau}}{\tau}}\right|^{2}=\mathcal{F}_{k}(x_{k+1})\leq\mathcal{F}_{k}(x_{k})=E(x_{k})+\frac{1}{2}\left|{\tfrac{x_{k}-x_{k-1}}{\tau}}\right|^{2}$

with a term on the right hand side that turns out to have entirely the wrong
scaling to estimate.555This is not surprising, as we are comparing a proper
approximately inertial solution with one that suddenly stops. A better
competitor might be the “straight continuation” $x_{k}+\tau(x_{k}-x_{k+1})$,
but then the estimate again requires convexity to deal with the energy-term.
Thus, the energy estimate cannot be obtained. ∎

The approach explained here turns out to be admissible for infinite-
dimensional spaces instead of $\mathbb{R}^{n}$ and even coupling between
Eulerian and Lagrangian coordinates. For that Eulerian-Lagrangian coupling
however additional difficulties appear that require some novel ideas on its
own. This will be discussed in the forthcoming sections.

######  Remark  (Numerical use of the method).

Since numerous numerical schemes for minimization (over discrete spaces) are
available the above methodology might also be attractive for computational
mathematics. The idea would here be to do a two-scale approximation: This
means that once $x_{k+1}^{\ell-1},x_{k}^{\ell-1}$ and $x^{\ell}_{k}$ are
constructed. Then we can define $x_{k+1}^{\ell}$ as the minimizer of

$\displaystyle\mathcal{F}_{k}^{\ell}(x):=E(x)+\frac{\tau}{2h}\left|{\frac{x-x_{k}^{\ell}}{\tau}-\frac{x_{k+1}^{\ell-1}-x_{k}^{\ell-1}}{\tau}}\right|^{2}.$

In order to pass to the limit it is in general unavoidable to use the
hyperbolic structure on a time-continuous level. This means that first
$\tau\to 0$ and only afterwards $h\to 0$. The question is how much smaller
does $\tau$ needs to be? One observes quickly, that in case $E$ is convex
$\tau$ and $h$ can be chosen arbitrarily. Hence, the smallness of $\tau$ in
relation should depend on the non-convexity of the assumed energies. In a
forthcoming paper we hope to investigate this issue further. ∎

##  4\. Proof of Theorem 2

In this section we prove Theorem 2, i.e. we prove existence of weak solutions
of the FPSI-problem (2.5). We follow the general strategy as outlined in
Section 3 and as devised in [6]. First we show existence for solutions to a
time-delayed, parabolic problem and then we use these solutions to construct
the weak solution to the actual problem. These two parts will be the topic of
the next two subsections.

###  4.1. Existence for the time-delayed problem

The key to finding a solution to the time-delayed problem is to understand
that on short time-scales, the time-delayed problem is ultimatively parabolic.
Instead of treating the delayed velocities $\partial_{t}\eta(t-h)$ and
$v(t-h)$ as non-local (in time) parts of the equation, we can treat them as a
fixed given data by considering times $t<h$ only. This mindset allows us to
state and prove the following Proposition 4.1 during this subsection.

Additionally, in this subsection and in Proposition 4.1, we need to introduce
regularizations of the involved stored energy and dissipation function
depending on the following parameters: We take $k_{0}\in\mathbb{N}$ in such a
way that $W^{k_{0},2}(\Omega;\mathbb{R}^{n})$ embeds into
$C^{1}(\Omega;\mathbb{R}^{n})\cap W^{2,q}(Q;\mathbb{R}^{n})$. Further we
choose $a_{0}\in(0,1)$ appropriately to be fixed later.

######  Proposition  (Existence for short time time-delayed solutions).

Given $h>0$ sufficiently small, $f,w\in
L^{2}([0,h]\times\Omega;\mathbb{R}^{n})$, $\zeta\in L^{2}([0,h]\times
Q;\mathbb{R}^{n})$, $\eta_{0}\in\mathcal{E}$ there exist

$\eta\in C^{0}([0,h];\overline{\mathcal{E}})\cap
W^{1,2}([0,h];W^{k_{0},2}(Q;\mathbb{R}^{n}))\text{ with }\eta(t)|_{P}=\gamma$

and $v\in W_{0}^{k_{0},2}([0,h]\times\Omega;\mathbb{R}^{n})$ with
$\operatorname{div}v=0$ such that

$\displaystyle\int_{0}^{h}\langle
DE(\eta),\phi\rangle+h^{a_{0}}\left\langle\nabla^{k_{0}}\eta,\nabla^{k_{0}}\phi\right\rangle_{Q}+\langle
D_{2}R(\eta,\partial_{t}\eta),\phi\rangle+h\left\langle\nabla^{k_{0}}\partial_{t}\eta,\nabla^{k_{0}}\phi\right\rangle_{Q}$
$\displaystyle+\langle
D_{2}A(\eta,\partial_{t}\eta-v\circ\eta),\phi-\xi\circ\eta\rangle+\nu\left\langle{\varepsilon}v,{\varepsilon}\xi\right\rangle_{\Omega}+h\left\langle\nabla^{k_{0}}v,\nabla^{k_{0}}\xi\right\rangle_{\Omega}$
$\displaystyle+\rho_{s}\left\langle\tfrac{\partial_{t}\eta-\zeta}{h},\phi\right\rangle_{Q}+\rho_{f}\left\langle\tfrac{v\circ\Phi-w}{h},\xi\circ\Phi\right\rangle_{\Omega}-\left\langle
f,\xi\right\rangle_{\Omega}dt=0$

for all $\phi\in L^{2}([0,h];W^{2,q}(Q;\mathbb{R}^{n}))$ with $\phi_{P}=0$ and
$\xi\in L^{2}([0,h];W^{1,2}_{0}(\Omega;\mathbb{R}^{n}))$ with
$\operatorname{div}\xi=0$. Here $\Phi:[0,h]\times\Omega\to\Omega$ is the flow
map of $v$, i.e. a family of volume preserving diffeomorphisms
$\Phi(t,.):\Omega\to\Omega$ such that $\Phi(0,.)=\textrm{id}$ and
$\partial_{t}\Phi(t,y)=v(t,\Phi(t,y))$.

In accordance with the regularizing parameters, we define

$\displaystyle
E_{h}(\eta):=E(\eta)+\frac{h^{a_{0}}}{2}\left\|\nabla^{k_{0}}\eta\right\|_{{Q}}^{2}\text{
and
}R_{h}(\eta,b):=R(\eta,b)+\frac{h}{2}\left\|\nabla^{k_{0}}b\right\|_{{Q}}^{2}$

as a shorthand and note that

$\displaystyle\left\langle DE_{h}(\eta),\phi\right\rangle=\left\langle
DE(\eta),\phi\right\rangle+h^{a_{0}}\left\langle\nabla^{k_{0}}\eta,\nabla^{k_{0}}\phi\right\rangle_{Q}$

and similarly for $\left\langle
D_{2}R_{h}(\eta,\partial_{t}\eta),\phi\right\rangle$. These regularizations
will help us at several points, both with establishing the energy inequality
as well as the flow map of the fluid.

As explained previously, the key to deal with such problems in a manner
consistent with our energy-considerations, is De Giorgi’s method of minimizing
movements [19]. We thus start the proof by minimizing the functional

$\displaystyle\mathcal{F}_{k}^{(\tau)}(\eta,v)$
$\displaystyle:=E_{h}(\eta)+\tau\left[R_{h}(\eta_{k}^{(\tau)},\partial_{k}^{\tau}\eta)+A(\eta_{k}^{(\tau)},\partial_{k}^{\tau}\eta-v\circ\eta_{k}^{(\tau)})+\frac{\nu}{2}\left\|{\varepsilon}v\right\|_{{\Omega}}^{2}+\frac{h}{2}\left\|\nabla^{k_{0}}v\right\|_{{\Omega}}^{2}\right]$
$\displaystyle+\frac{\tau}{2h}\left[\rho_{s}\left\|\partial_{k}^{\tau}\eta-\zeta_{k}^{(\tau)}\right\|_{{Q}}^{2}+\rho_{f}\left\|v\circ\Phi_{k}^{(\tau)}-w_{k}^{(\tau)}\right\|_{{\Omega}}^{2}\right]-\tau\left\langle
f_{k}^{(\tau)},v\right\rangle_{\Omega}$

in the class of all admissible $\eta,v$, i.e. those such that
$\eta\in\mathcal{E}_{h}=\mathcal{E}\cap W^{k_{0},2}(Q;\mathbb{R}^{n})$ and
$v\in W_{0,\mathrm{div}}^{1,2}(\Omega;\mathbb{R}^{n})$, i.e. $v\in
W_{0}^{1,2}(\Omega;\mathbb{R}^{n})$ with $\operatorname{div}v=0$. Here
$\partial_{k}^{\tau}\eta:=\smash{\tfrac{\eta-\eta_{k}^{(\tau)}}{\tau}}$ is a
discrete derivative and
$\smash{f_{k}^{(\tau)}}:=\smash{\fint_{k\tau}^{(k+1)\tau}fdt}$ is a time-
discretization of $f$. In the same way we define $\smash{\zeta_{k}^{(\tau)}}$
and $\smash{w_{k}^{(\tau)}}$. For any given
$\eta_{k}^{(\tau)},\Phi_{k}^{(\tau)}$, this minimizer will yield the next
step, i.e. the pair $(\eta_{k+1}^{(\tau)},v_{k+1}^{(\tau)})$. From this we
then also construct $\Phi_{k+1}^{(\tau)}:=(\textrm{id}+\tau
v_{k+1}^{(\tau)})\circ\Phi_{k}^{(\tau)}$.666Note that, as expected for a
parabolic problem, the rate variables, i.e. the velocities
$\partial_{k}^{\tau}\eta^{(\tau)}_{k}$ and $v_{k}^{(\tau)}$ are almost
entirely discarded. Only the latter occurs indirectly in $\Phi_{k}^{(\tau)}$.
Similarly we are not using the initial data for these velocities. Only later
will they all reappear through proper choice of $\zeta$ and $w$. This again
emphasizes that these rate variables are associated with the $h$-scale which
is kept constant for the time-delayed problem.

The existence of a minimizer to $\mathcal{F}_{k}^{(\tau)}$ follows by the
direct method from the calculus of variations. Indeed, by the coercivity of
$\smash{\mathcal{F}_{k}^{(\tau)}}$, we find a bounded sequence
$(\eta_{k}^{l},v_{k}^{l})\in\mathcal{E}_{h}\times\smash{W^{1,2}_{0,\mathrm{div}}(\Omega;\mathbb{R}^{n})}$,
such that

$\mathcal{F}_{k}^{(\tau)}(\eta_{k}^{l},v_{k}^{l})\to\mathrm{inf}_{(\nu,\eta)\in\mathcal{E}_{h}\times
W^{1,2}_{0,\mathrm{div}}(\Omega;\mathbb{R}^{n})}\mathcal{F}_{k}^{(\tau)}(\eta,v)$

Here, the boundedness of sequence is particularly due to the added
regularizing terms, though existence of minimizers could be shown even if they
were omitted, as well as the highest order term in the elastic stored energy.

By the Banach-Alaoglu theorem, we conclude the existence of a pair
$(\eta_{k}^{\tau},v_{k}^{\tau})\in W^{k_{0},2}(Q;\mathbb{R}^{n})\cap
W^{2,q}(Q;\mathbb{R}^{n})\times W^{1,2}(\Omega,\mathbb{R}^{n})$ such that

$(\eta_{k}^{l},v_{k}^{l})\rightharpoonup(\eta_{k}^{\tau},v_{k}^{\tau})\in
W^{k_{0},2}(Q;\mathbb{R}^{n})W^{2,q}(Q;\mathbb{R}^{n})\times
W^{1,2}(\Omega,\mathbb{R}^{n}).$

By compact embeddings, owing to the regularizing terms, we can show that
$\eta_{k}^{\tau}$ fulfills the Ciarlet-Nečas condition and thus lies in
$\mathcal{E}_{h}$ 777This can be show even under a weaker type of convergence;
see the original paper [14]., due to linearity $v_{k}^{\tau}$ is divergence
free and satisfies the zero boundary condition.

It remains to show that $\mathcal{F}_{k}^{\tau}$ is weakly lower
semicontinuous with respect to the weak convergence above. This is
nevertheless standard, due to the convexity of the highest order terms, the
strong convergence of $\eta_{k}^{l}\to\eta_{k}^{\tau}$ in
$C^{1}(Q;\mathbb{R}^{n})$. Thus, $(\eta_{k}^{\tau},v_{k}^{\tau})$ is the
sought minimizer.

The reason we prefer a minimization to other approaches is that it allows us
to immediately access a discrete energy inequality, without having to rely on
additional properties of the energy, such as convexity.

For this we compare the minimal value, i.e. that of
$(\eta_{k+1}^{(\tau)},v_{k+1}^{(\tau)})$ to that of “standing still”, i.e.
$(\eta_{k}^{(\tau)},0)$ in the same functional. This removes all dissipative
terms on one side of the inequality and gives us

$\displaystyle\phantom{{}={}}E_{h}(\eta_{k+1}^{(\tau)})+\tau\left[R_{h}(\eta_{k}^{(\tau)},\partial_{k}^{\tau}\eta_{k+1}^{(\tau)})+A(\eta_{k}^{(\tau)},\partial_{k}^{\tau}\eta_{k+1}^{(\tau)}-v_{k+1}^{(\tau)}\circ\eta_{k}^{(\tau)})+\frac{\nu}{2}\left\|{\varepsilon}v_{k+1}^{(\tau)}\right\|_{{\Omega}}^{2}+\frac{h}{2}\left\|\nabla^{k_{0}}v^{(\tau)}_{k+1}\right\|_{{\Omega}}^{2}\right]$
$\displaystyle+\frac{\tau}{2h}\left[\rho_{s}\left\|\partial^{\tau}_{k}\eta_{k+1}^{(\tau)}-\zeta_{k}^{(\tau)}\right\|_{{Q}}^{2}+\rho_{f}\left\|v_{k+1}^{(\tau)}\circ\Phi_{k}^{(\tau)}-w_{k}^{(\tau)}\right\|_{{\Omega}}^{2}\right]-\tau\left\langle
f_{k}^{(\tau)},v_{k+1}^{(\tau)}\right\rangle_{\Omega}$ $\displaystyle\leq
E_{h}(\eta_{k}^{(\tau)})+\frac{\tau}{2h}\left[\rho_{s}\left\|\zeta_{k}^{(\tau)}\right\|_{{Q}}^{2}+\rho_{f}\left\|w_{k}^{(\tau)}\right\|_{{\Omega}}^{2}\right]$

A telescope argument and the usual weighted Young’s inequality then gives us
an energy estimate as well.

$\displaystyle
E_{h}(\eta_{N}^{(\tau)})+\\!\sum_{k=0}^{N-1}\tau\left[R_{h}(\eta_{k}^{(\tau)},\partial^{\tau}_{k}\eta_{k+1}^{(\tau)})+A(\eta_{k}^{(\tau)},\partial_{k}^{\tau}\eta_{k+1}^{(\tau)}-v_{k+1}^{(\tau)}\circ\eta_{k}^{(\tau)})+\frac{\nu}{2}\left\|{\varepsilon}v_{k+1}^{(\tau)}\right\|_{{\Omega}}^{2}\\!+\frac{h}{2}\left\|\nabla^{k_{0}}v^{(\tau)}_{k+1}\right\|_{{\Omega}}^{2}\right]$
$\displaystyle+\sum_{k=0}^{N-1}C\tau\left[\rho_{s}\left\|\partial^{\tau}_{k}\eta_{k+1}^{(\tau)}\right\|_{{Q}}^{2}+\rho_{f}\left\|v_{k+1}^{(\tau)}\circ\Phi_{k}^{(\tau)}\right\|_{{\Omega}}^{2}\right]$
$\displaystyle\leq
E_{h}(\eta_{k}^{(\tau)})+C\sum_{k=0}^{n-1}\frac{\tau}{2h}\left[\rho_{s}\left\|\smash{\zeta_{k}^{(\tau)}}\right\|_{{Q}}^{2}+\rho_{f}\left\|\smash{w_{k}^{(\tau)}}\right\|_{{\Omega}}^{2}+\left\|\smash{f_{k}^{(\tau)}}\right\|_{{\Omega}}^{2}\right]$

This is the first of many similar energy estimates related to the physical
energy inequality, which we will make use of.

Note also that in each step this also allows us to derive a uniform estimate
on the distance from the initial data in the form

$\displaystyle\left\|\eta_{0}-\smash{\eta_{N}^{(\tau)}}\right\|_{{Q}}\leq\sum_{k=0}^{N-1}\tau\left\|\partial^{\tau}_{k}\eta^{(\tau)}_{k+1}\right\|_{{Q}}\leq\sqrt{\sum_{k=0}^{N-1}\tau}\sqrt{\sum_{k=0}^{N-1}\tau\left\|\partial^{\tau}_{k}\eta^{(\tau)}_{k+1}\right\|_{{Q}}^{2}}\leq
C\sqrt{h}.$

Since a straightforward calculation (see [6, Prop. 2.7]) shows that any given
state of finite energy has a minimum $L^{2}$-distance to any (self-)coliding
state of finite energy, this also tells us that we can indeed choose $h$ small
enough to avoid such a collision.

Next, we define the follwoing piecewise constant and piecewise affine (in
time) approximations:

$\displaystyle\bar{\eta}^{(\tau)}(t,x)$
$\displaystyle:=\eta_{k+1}^{(\tau)}(x)$ for $\displaystyle t\in[\tau
k,\tau(k+1)),x\in Q$ $\displaystyle\underline{\eta}^{(\tau)}(t,x)$
$\displaystyle:=\eta_{k}^{(\tau)}(x)$ for $\displaystyle t\in[\tau
k,\tau(k+1)),x\in Q$ $\displaystyle\hat{\eta}^{(\tau)}(t,x)$
$\displaystyle:=((k+1)-t/\tau)\eta_{k}^{(\tau)}(x)-(t/\tau-k)\eta_{k+1}^{(\tau)}(x)$
for $\displaystyle t\in[\tau k,\tau(k+1)),x\in Q$ $\displaystyle
v^{(\tau)}(t,y)$ $\displaystyle:=v_{k}^{(\tau)}(y)$ for $\displaystyle
t\in[\tau k,\tau(k+1)),y\in\Omega$ $\displaystyle\Phi^{(\tau)}(t,y)$
$\displaystyle:=\Phi_{k}^{(\tau)}(y)$ for $\displaystyle t\in[\tau
k,\tau(k+1)),y\in\Omega$ in particular we have
$\displaystyle\partial_{t}\hat{\eta}^{(\tau)}(t,x)$
$\displaystyle=\partial^{\tau}_{k}\eta_{k+1}^{(\tau)}(x)$ for $\displaystyle
t\in(\tau k,\tau(k+1)),x\in\Omega$

Plugging these quantities into the energy estimate above, we get that there is
a constant $C$ independent of $\tau$ such that

$\displaystyle\sup_{t\in[0,h]}E_{h}(\bar{\eta}^{(\tau)}(t))\leq C,$
$\displaystyle\quad\sup_{t\in[0,h]}\left\|\bar{\eta}^{(\tau)}(t)\right\|_{{W^{k_{0},2}(Q)}}\leq
C,$
$\displaystyle\quad\sup_{t\in[0,h]}\left\|\hat{\eta}^{(\tau)}(t)\right\|_{{W^{k_{0},2}(Q)}}\leq
C$
$\displaystyle\int_{0}^{h}\left\|\partial_{t}\hat{\eta}^{(\tau)}(t)\right\|_{{W^{k_{0},2}(Q)}}^{2}dt\leq
C$
$\displaystyle\quad\int_{0}^{h}\left\|v^{(\tau)}\right\|_{{W^{k_{0},2}(\Omega)}}^{2}dt\leq
C,$
$\displaystyle\int_{0}^{h}A(\eta,v^{(\tau)}\circ\bar{\eta}^{(\tau)}-\partial_{t}\hat{\eta}^{(\tau)})dt\leq
C.$

Thus we can find a (non-relabelled) subsequence of $\tau$’s and $\eta\in
W^{1,2}([0,T];W^{k_{0},2}(Q;\mathbb{R}^{n}))$ and $v\in
L^{2}([0,T];W^{k_{0},2}(\Omega;\mathbb{R}^{n}))$ such that

$\displaystyle\bar{\eta}^{(\tau)},\underline{\eta}^{(\tau)}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\eta\text{
in }L^{\infty}([0,T];W^{k_{0},2}(Q;\mathbb{R}^{n}))\quad\text{ and
}\quad\hat{\eta}^{(\tau)}\rightharpoonup\eta\text{ in
}W^{1,2}([0,T];W^{k_{0},2}(Q;\mathbb{R}^{n}))$ $\displaystyle
v^{(\tau)}\rightharpoonup v\text{ in
}L^{2}([0,T];W^{k_{0},2}(\Omega;\mathbb{R}^{n})).$

The next question we are going to have to deal with is how to construct the
flow map $\Phi$. In each step, the flow map is updated by concatenation. Not
only is this in itself an inherently non-linear procedure, but also for any
non-zero time, the number of those steps goes to infinity when we send $\tau$
to zero. On the other hand, each of these steps represents a change of scale
$\tau$. This hints at an exponential structure and this is indeed what we
find. As often (e.g. when working with conservation laws) the Lipschitz
continuity of the flow-velocity is essential to have a well defined Lagrangian
flow map. Here we particularly rely on the regularization of the flow-velocity
is used at the $\tau$ level.

######  Lemma  (Existence of a flow map).

Assume that $v^{(\tau)}\in L^{2}([0,h];C^{0,1}(\Omega;\mathbb{R}^{n}))$
uniformly in $\tau$, such that $\operatorname{div}v^{(\tau)}=0$. Then the maps
$\Phi^{(\tau)}$ are uniformly Lipschitz continuous in space, independent of
$\tau$ (but not of $h$). Additionally we have in the limit that

$\displaystyle\det\nabla\Phi^{(\tau)}\to 1.$

###### Proof.

Let us illustrate the proof for the Jacobian determinant. Here we are using an
estimate involving an expansion of the determinant, the inequality between
arithmetric and geometric mean, as well as the fact that
$(1+a/N)^{N}\nearrow\exp(a)$:

$\displaystyle\phantom{{}={}}\det\nabla\Phi_{N}^{(\tau)}=\prod_{k=1}^{N}\left[\det\left(I+\tau\nabla
v_{k}^{(\tau)}\right)\right]\circ\Phi_{k-1}^{(\tau)}$
$\displaystyle=\prod_{k=1}^{N}\bigg{[}1+\tau\underbrace{\operatorname{tr}\left(\nabla
v_{k}^{(\tau)}\right)}_{=\operatorname{div}v_{k}^{(\tau)}=0}+\sum_{l=2}^{n}\tau^{l}M_{l}\left(\nabla
v_{k}^{(\tau)}\right)\bigg{]}\circ\Phi_{k-1}^{(\tau)}$
$\displaystyle\leq\prod_{k=1}^{N}\left[1+\sum_{l=2}^{n}c\tau^{l}\operatorname{Lip}\left(\nabla
v_{k}^{(\tau)}\right)^{l}\right]\leq\left(1+\frac{1}{N}\sum_{k=1}^{N}\sum_{l=2}^{n}\tau^{l}c\operatorname{Lip}\left(v_{k}^{(\tau)}\right)^{l}\right)^{N}$
$\displaystyle\leq\exp\left(C\tau\sum_{k=1}^{N}\tau\operatorname{Lip}\left(v_{k}^{(\tau)}\right)^{2}\right),$

where $M_{l}$ includes all terms of order $l$ in the polynomial expansion of
the determinant (e.g. $M_{1}(A)=\operatorname{tr}(A)$, …,
$M_{n-1}(A)=\operatorname{tr}(\operatorname{cof}(A))$, $M_{n}(A)=\det A$, cf.
e.g. [6, Lemma A.1.]). Now the sum in the last term is nothing but $\tau$
times the $L^{2}([0,h];C^{0,1}(\Omega;\mathbb{R}^{n}))$-norm of $v^{(\tau)}$,
something on which we have $\tau$-independent bounds by assumption. Similar
arguments also derive a bound from below, which gives us
$\det\nabla\Phi_{N}^{(\tau)}\to 1$ uniformly. Furthermore applying a very
similar argument to $\left|{\smash{\nabla\Phi^{(\tau)}_{N}}}\right|$ gives us
a $\tau$-independent bound on the Lipschitz constant of $\Phi$. ∎

In particular, since $\Phi^{(\tau)}$ is piecewise constant, this allows us to
use a variant of the Arzela-Ascoli theorem to find a limit
$\Phi:[0,h]\times\Omega\to\Omega$ (after a subsequence), which has to be a
volume preserving diffeomorphism with

$\displaystyle\partial_{t}\Phi(t,y)=v(t,\Phi(t,y))\text{ for all
}t\in[0,h],y\in\Omega$

i.e. a flow map.

With this, we have one half of the proof in hand, namely constructing the
objects that are going to form our solution. The other half is showing that
these objects actually are a solution, i.e. that they satisfy the right
equation.

For this we can follow along the same steps as the construction. Apart from
giving us an energy inequality, the other benefit of starting with a
minimization is that any such minimizer has to satisfy the Euler-Lagrange
equation.

$\displaystyle\left\langle DE_{h}(\eta_{k+1}),\phi\right\rangle+\left\langle
D_{2}R_{h}(\eta_{k},\partial^{\tau}\eta_{k+1}),\phi\right\rangle+\left\langle
D_{2}A(\eta_{k},\partial^{\tau}\eta_{k+1}-v_{k+1}\circ\eta_{k}),\phi-\xi\circ\eta_{k}\right\rangle$
$\displaystyle\quad+\nu\left\langle\varepsilon
v_{k+1},\varepsilon\xi\right\rangle_{\Omega}+h\left\langle\nabla^{k_{0}}v_{k+1},\nabla^{k_{0}}\xi\right\rangle_{\Omega}+\rho_{s}\left\langle\tfrac{\partial^{\tau}\eta-\zeta_{k}}{h},\phi\right\rangle_{Q}+\rho_{f}\left\langle\tfrac{v_{k+1}\circ\Phi_{k}-w_{k}}{h},\xi\right\rangle_{\Omega}+\left\langle
f,\xi\right\rangle_{\Omega}=0$

for any pair $(\phi,\xi)$ such that
$(\eta_{k+1}+\varepsilon\phi,v_{k+1}+\varepsilon\xi/\tau)$ is also admissible
for our minimization888Note that we scale $\phi$ and $\xi$ differently in
$\tau$, so that they both behave like a change of position. This is not
required at this point, but convenient to do for the limit passage. for some
$\varepsilon_{0}$ and any $\varepsilon\in[-\varepsilon_{0},\varepsilon_{0}]$,
$\phi\in W^{1,2}(Q;\mathbb{R}^{n})$, $\xi\in
W_{0}^{1,2}(\Omega;\mathbb{R}^{n})$ with $\phi|_{P}=0$ and
$\operatorname{div}\xi=0$. Here, we point out that, as we start with a
configuration without contact we know due to the bounds in the energy estimate
that, for $h$ small enough, every minimizer will be contact-free. Thus, indeed
variations in all directions are possible.

Next we replace all occurences of the discrete approximations with the
corresponding time dependent functions, giving us

$\displaystyle\left\langle
DE_{h}(\bar{\eta}^{(\tau)}),\phi\right\rangle+\left\langle
D_{2}R_{h}(\underline{\eta}^{(\tau)},\partial_{t}\hat{\eta}^{(\tau)}),\phi\right\rangle+\left\langle
D_{2}A(\underline{\eta}^{(\tau)},\partial_{t}\hat{\eta}^{(\tau)}-v^{(\tau)}\circ\underline{\eta}^{(\tau)}),\phi-\xi\circ\underline{\eta}^{(\tau)}\right\rangle$
$\displaystyle\quad+\nu\left\langle\varepsilon
v^{(\tau)},\varepsilon\xi\right\rangle_{\Omega}+h\left\langle\nabla^{k_{0}}v^{(\tau)},\nabla^{k_{0}}\xi\right\rangle_{\Omega}$
$\displaystyle\quad+\rho_{s}\left\langle\tfrac{\partial_{t}\hat{\eta}^{(\tau)}-\zeta_{k}}{h},\phi\right\rangle_{Q}+\rho_{f}\left\langle\tfrac{v^{(\tau)}\circ\Phi^{(\tau)}-w_{k}}{h},\xi\circ\Phi^{(\tau)}\right\rangle_{\Omega}+\left\langle
f,\xi\right\rangle_{\Omega}=0$

for almost all times $t$.

Note in particular that there is now no longer any direct occurrence of $\tau$
in this equation. Thus one can form the integral over $[0,h]$ and pass to the
limit $\tau\to 0$. Due to the added regularizing terms this is
straightforward. Indeed, the highest terms (added through linearization) are
linear. For all other terms (except the
$\frac{1}{\mathrm{det}(\nabla\eta)^{a}}$ for which we rely on uniform
continuity) we exploit compact embeddings and Nemitskii continuity (see e.g.
[30]) and thus obtain:

$\displaystyle\int_{0}^{h}\left\langle
DE_{h}(\eta),\phi\right\rangle+\left\langle
D_{2}R_{h}(\eta,\partial_{t}\eta),\phi\right\rangle+\left\langle
D_{2}A(\eta,\partial_{t}\eta-v\circ\eta),\phi-\xi\circ\eta\right\rangle$
$\displaystyle\quad+\nu\left\langle\varepsilon
v,\varepsilon\xi\right\rangle_{\Omega}+h\left\langle\nabla^{k_{0}}v,\nabla^{k_{0}}\xi\right\rangle_{\Omega}+\rho_{s}\left\langle\tfrac{\partial_{t}\eta-\zeta_{k}}{h},\phi\right\rangle_{Q}+\rho_{f}\left\langle\tfrac{v\circ\Phi-w}{h},\xi\circ\Phi\right\rangle_{\Omega}+\left\langle
f,\xi\right\rangle_{\Omega}dt=0$

Finally we need to show that the new solution also satisfies the energy
inequality. There is a way to do so directly from the approximation using De
Giorgi’s methods.999Note that our previous energy estimate was off by a factor
2 in all dissipation terms, so we cannot simply apply lower-semicontinuity to
this. However in the situation that we are studying, due to the regularizer,
we can directly show that the energy inequality holds true for all solutions.

For this, we test the equation directly with the pair $(\partial_{t}\eta,v)$,
which we are allowed due to our regularizing terms and we obtain

$\displaystyle E_{h}(\eta(h))+\int_{0}^{h}\left\langle
D_{2}R_{h}(\eta,\partial_{t}\eta),\partial_{t}\eta\right\rangle+2A(\eta,\partial_{t}\eta-v\circ\eta)+\nu\left\|\varepsilon
v\right\|_{{\Omega}}^{2}+h\left\|\nabla^{k_{0}}v\right\|_{{\Omega}}^{2}dt$
$\displaystyle\quad+\fint_{0}^{h}\rho_{s}\left\langle\partial_{t}\eta-\zeta,\partial_{t}\eta\right\rangle_{Q}+\rho_{f}\left\langle
v\circ\Phi-w,v\circ\Phi\right\rangle_{\Omega}dt+\int_{0}^{h}\left\langle
f,v\right\rangle_{\Omega}dt=E_{h}(\eta(0))$

Now we apply a consequence of Young’s inequality, $\left\langle
a-b,a\right\rangle=\left|{a}\right|^{2}-\left\langle
b,a\right\rangle\geq\frac{1}{2}\left|{a}\right|^{2}-\frac{1}{2}\left|{b}\right|^{2}$,
to the two inertial terms. Together with the conservation of volume implying
$\left\|v\circ\Phi\right\|_{{\Omega}}^{2}=\left\|v\right\|_{{\Omega}}^{2}$,
this results in the energy inequality:

$\displaystyle
E_{h}(\eta(h))+\int_{0}^{h}2R_{h}(\eta,\partial_{t}\eta)+2A(\eta,\partial_{t}\eta-v\circ\eta)+\nu\left\|\varepsilon
v\right\|_{{\Omega}}^{2}+h\left\|\nabla^{k_{0}}v\right\|_{{\Omega}}^{2}dt$
$\displaystyle+\fint_{0}^{h}\frac{\rho_{s}}{2}\left\|\partial_{t}\eta\right\|_{{Q}}^{2}+\frac{\rho_{f}}{2}\left\|v\right\|_{{\Omega}}dt=E_{h}(\eta(0))+\fint_{0}^{h}\frac{\rho_{s}}{2}\left\|\zeta\right\|_{{Q}}^{2}+\frac{\rho_{f}}{2}\left\|w\right\|_{{\Omega}}dt+\int_{0}^{h}\left\langle
f,v\right\rangle_{\Omega}dt.$

###  4.2. Convergence to the full problem

Taking a look at the solutions constructed in the previous section, we can use
the energy inequality to see that a solution on the interval $[0,h]$ can be
turned into admissible initial and right hand side data with $\eta(h)$ as the
new $\eta_{0}$, $\partial_{t}\eta$ as $\zeta$,
$v(t)\circ\Phi(t)\circ\Phi(h)^{-1}$ as $w$ (to correct for the flow) and so
on. So we can apply the result again to construct a solution on the interval
$[h,2h]$ and so on. We do so and combine these solutions into a single pair of
functions $\eta^{(h)}:[0,T]\times Q\to\Omega$ and
$v^{(h)}:[0,T]\times\Omega\to\mathbb{R}^{n}$. These then fulfill the full
time-delayed equation:

(4.1) $\displaystyle\begin{aligned}
&\int_{0}^{T}\rho_{s}\left\langle\tfrac{\partial_{t}\eta^{(h)}(t)-\eta^{(h)}(t-h)}{h},\phi\right\rangle_{Q}+\rho_{f}\left\langle\tfrac{v^{(h)}(t)-v^{(h)}(t-h)\circ\Phi_{-h}^{(h)}(t)}{h},\xi\right\rangle_{\Omega}\,dt\\\
&\quad+\int_{0}^{T}\left\langle
DE_{h}(\eta^{(h)}),\phi\right\rangle_{Q}+\left\langle
D_{2}A(\eta,\partial_{t}\eta^{(h)}-v^{(h)}\circ\eta^{(h)}),\phi-\xi\circ\eta^{(h)}\right\rangle_{Q}\\\
&\quad+\int_{0}^{T}\left\langle
D_{2}R_{h}(\eta^{(h)},\partial_{t}\eta^{(h)}),\phi\right\rangle_{Q}+\nu\left\langle\varepsilon
v^{(h)},\varepsilon\xi\right\rangle_{\Omega}+h\left\langle\nabla^{k_{0}}v^{(h)},\nabla^{k_{0}}\xi\right\rangle\,dt=\int_{0}^{T}\left\langle
f,\xi\right\rangle_{\Omega}dt\end{aligned}$

Special care needs to be taken in constructing all the fluid related terms
along the way. For the flow map we take $\Phi_{s}^{(h)}(t,.)$ to be the map
that maps a fluid particle’s position at time $t$ to its position at time
$t+s$. Such a map is constructed from composing the flow maps on each of the
subintervals, by first using their inverse to get to a multiple of $h$ and
then moving forward again.

Or in other words, if $\Phi^{0},\Phi^{1},\dots$ denote the flow maps on the
respective intervals $[0,h],[h,2h],\dots$ then we first define for
$t\in[lh,(l+1)h]$

$\displaystyle\Phi_{t}^{(h)}(0):=\Phi^{l}(t-lh)\circ\Phi^{l-1}(h)\circ\dots\circ\Phi^{0}(h)$

and then for arbitrary $t,t+s\in[0,T]$ using the invertibility of the flow
maps

$\displaystyle\Phi^{(h)}_{s}(t):=\Phi^{(h)}_{t+s}(0)\circ\Phi^{(h)}(t)^{-1}.$

In practice of course we are mostly interested in the case when $t$ and $t+s$
are at most a distance $h$ apart and most of the back and forth in this
definition cancels.

In particular these maps are flow maps for $v^{h}$, so we have

$\displaystyle\partial_{s}\Phi_{s}^{(h)}(t,y)=v^{(h)}(t+s,\Phi_{s}^{(h)}(t,y))\text{
for all }t,t+s\in[0,T],y\in\Omega$

as well as $\Phi_{0}^{(h)}(t,.)=\textrm{id}$ and
$\det\nabla\Phi_{s}^{(h)}(t,.)=1$.

Additionally, a simple iteration of the energy-inequality yields the following
energy inequality for the problem on $[0,T]$:

$\displaystyle
E_{h}(\eta^{(h)}(t_{0}))+\int_{0}^{t_{0}}2R_{h}(\eta^{(h)},\partial_{t}\eta^{(h)})+2A(\eta,\partial_{t}\eta^{(h)}-v^{(h)}\circ\eta^{(h)})+\nu\left\|\varepsilon
v^{(h)}\right\|_{{\Omega}}^{2}+h\left\|\nabla^{k_{0}}v^{(h)}\right\|^{2}dt$
$\displaystyle\quad+\fint^{t_{0}}_{t_{0}-h}\frac{\rho_{s}}{2}\left\|\partial_{t}\eta^{(h)}\right\|_{{Q}}^{2}+\frac{\rho_{f}}{2}\left\|v^{(h)}\right\|_{{\Omega}}^{2}dt+\int_{0}^{t_{0}}\left\langle
f,v^{(h)}\right\rangle_{\Omega}dt\leq
E_{h}(\eta_{0})+\frac{\rho_{s}}{2}\left\|b\right\|_{{Q}}^{2}+\frac{\rho_{f}}{2}\left\|v_{0}\right\|_{{\Omega}}^{2}.$

Similarly to the estimate for the discrete approximation, this again gives us
a uniform bound on $\left\|\eta_{0}-\eta^{(h)}(t_{0})\right\|_{{Q}}$ in terms
of $t_{0}$, which we can use to initially choose $T$ small enough to avoid
collisions.

As a consequence, we conclude that there is a constant $C$ independent of $h$
such that

$\displaystyle\sup_{t\in[0,T]}\Big{(}E_{h}(\eta^{(h)}(t))+\left\|\eta^{(h)}(t)\right\|_{{W^{2,q}(Q)}}+h^{a_{0}}\left\|\eta^{(h)}(t)\right\|_{{W^{k_{0},2}(Q)}}\Big{)}\leq
C,$
$\displaystyle\int_{0}^{T}\Big{(}\left\|\partial_{t}\eta^{(h)}(t)\right\|_{{W^{1,2}(Q)}}^{2}+h\left\|\partial_{t}\eta^{(h)}(t)\right\|_{{W^{k_{0},2}(Q)}}^{2}+\left\|v^{(h)}\right\|_{{W^{1,2}(\Omega)}}^{2}+h\left\|v^{(h)}\right\|_{{W^{k_{0},2}(\Omega)}}^{2}\Big{)}dt\leq
C,$
$\displaystyle\int_{0}^{T}A(\eta,v^{(h)}\circ\eta^{(h)}-\partial_{t}\eta^{(h)})dt\leq
C.$

These can again be used to pick weak*-converging sub-sequences and a limit
pair $(\eta,v)$. Our goal is again to prove convergence of the equation and
these weak convergences are indeed enough to do so for all the terms which are
unchanged from before. However this time, we have to additionaly deal with the
two $h$-dependent inertial terms, of which in particular the fluid-term
requires some attention as weak convergence is not enough to go to the limit.

Interestingly in this case, the correct quantities to consider are time-
averages. Consider first the solid. Then we can define the rolling average of
the momentum as

$\displaystyle
m^{(h)}(t,x)=\rho_{s}\fint_{t-h}^{t}\partial_{t}\eta^{(h)}(s,x)ds$

and note that its time derivative

$\displaystyle\partial_{t}m^{(h)}(t,x)=\rho_{s}\frac{\partial_{t}\eta^{(h)}(t)-\partial_{t}\eta^{(h)}(t-h)}{h}$

is something already occuring in the equation and over which we thus already
have some weak control.

The same is true for the fluid, where we however have to be a bit more careful
and instead of the physically correct quantity $\fint\rho_{f}v\circ\Phi$, we
instead have to consider a similarly “straightened” version

$\displaystyle m^{(h)}(t,x)=\rho_{f}\fint_{t-h}^{t}v(s,x)ds$

For these averaged momenta, since we have $L^{2}(W^{1,2})$-bounds from the
energy estimate and bounds on their time-derivatives in negative spaces from
the equation, we can apply a version of the Aubin-Lions-lemma (see e.g. [30,
Sec. 7.3]) to obtain their strong $L^{2}$-convergence (in space-time). In
particular, for the fluid this implies

$\displaystyle\int_{0}^{T}\left\langle\frac{v^{(h)}(t)-v^{(h)}(t-h)\circ\Phi_{-h}^{(h)}(t)}{h},\xi(t)\right\rangle
dt=-\int_{0}^{T}\left\langle
v^{(h)}(t),\frac{\xi(t+h)\circ\Phi_{h}^{(h)}(t)-\xi(t)}{h}\right\rangle dt$
$\displaystyle=-\int_{0}^{T}\left\langle
v^{(h)}(t),\fint_{0}^{h}\partial_{s}\xi(t+s)\circ\Phi_{s}^{(h)}(t)ds\right\rangle
dt$ $\displaystyle=-\int_{0}^{T}\left\langle
v^{(h)}(t),\fint_{0}^{h}\left[\partial_{t}\xi(t+s)+v^{(h)}(t+s)\cdot\nabla\xi(t+s)\right]\circ\Phi_{s}^{(h)}(t)ds\right\rangle
dt$ $\displaystyle\to-\int_{0}^{T}\left\langle
v(t),\partial_{t}\xi+v\cdot\nabla\xi\right\rangle dt.$

Next note that a Minty-type argument (compare [6, Lem. 3.9]) allows us to deal
with the second order term in $DE(\eta^{(h)})$ while the convergence of the
first order terms results from compact embeddings and the boundedness of the
determinant. In total implies the weak convergence of $DE(\eta^{(h)})$ to
$DE(\eta)$. Conclusively we then have

$\displaystyle\int_{0}^{T}\left\langle DE(\eta),\phi\right\rangle+\left\langle
D_{2}R(\eta,\partial_{t}\eta),\phi\right\rangle+\left\langle
D_{2}A(\eta,\partial_{t}\eta-v\circ\eta),\phi-\xi\circ\eta\right\rangle$
$\displaystyle\quad+\nu\left\langle\varepsilon
v,\varepsilon\xi\right\rangle-\rho_{s}\left\langle\partial_{t}\eta,\partial_{t}\phi\right\rangle-\rho_{f}\left\langle
v,\partial_{t}\xi-v\cdot\nabla\xi\right\rangle+\left\langle f,\xi\right\rangle
dt=0$

for all $\phi\in C^{\infty}([0,T]\times Q;\mathbb{R}^{n})$ and $\xi\in
C_{0}^{\infty}([0,T]\times\Omega;\mathbb{R}^{n})$ with $\phi|_{[0,T]\times
P}=0$ and $\operatorname{div}\xi=0$, as desired.

This initially holds only on a small interval on which collisions are avoided
a-priori, but by standard arguments, we can extend this interval until we
reach the original $T$ or an actual collision.

Finally, we can spare some thoughts on the pressure, which can be found
through the equation. For this we need the so called Bogovskiĭ-operator which
is a bounded linear operator $\mathcal{B}:W^{k,p}(\Omega)\to
W_{0}^{k+1,p}(\Omega;\mathbb{R}^{n})$ such that
$\operatorname{div}\mathcal{B}\psi=\psi$. With this we can define the operator

$\displaystyle P(\psi):=$ $\displaystyle\int_{0}^{T}\Big{(}\left\langle
D_{2}A(\eta,\partial_{t}\eta-v\circ\eta),\mathcal{B}(\psi)\circ\eta\right\rangle-\nu\left\langle\varepsilon
v,\varepsilon\mathcal{B}(\psi)\right\rangle$
$\displaystyle\quad+\rho_{f}\left\langle
v,\partial_{t}\mathcal{B}(\psi)-v\cdot\nabla\mathcal{B}(\psi)\right\rangle-\left\langle
f,\mathcal{B}\psi\right\rangle\Big{)}dt$

for which a short calculation reveals that it is a bounded distribution in
space-time which gives the pressure of the fluid. Indeed, adding $\nabla
p(\xi):=-P(\operatorname{div}(\xi))$ to the equation allows us to test with
testfunctions of non-zero divergence and shows that the solution
(distributionally) satisfies the Navier-Stokes equations.

### Declarations

This work has been supported by the Primus research programme PRIMUS/19/
SCI/01, the University Centre UNCE/SCI/023 of Charles University as well as
the grant GJ19-11707Y of the Czech national grant agency (GAČR). Further S. S.
wishes to thank the University of Vienna for their kind hospitality in winter
2020/21. There is no conflict of interest.

## References

  * Allaire [1990] Grégoire Allaire.  Homogenization of the Navier-Stokes equations in open sets perforated with tiny holes. II. Noncritical sizes of the holes for a volume distribution and a surface distribution of holes.  _Arch. Rational Mech. Anal._ , 113(3):261–298, 1990.  doi: 10.1007/BF00375066.
  * Ambartsumyan et al. [2019] Ilona Ambartsumyan, Vincent J Ervin, Truong Nguyen, and Ivan Yotov.  A nonlinear stokes–biot model for the interaction of a non-newtonian fluid with poroelastic media.  _ESAIM: Mathematical Modelling and Numerical Analysis_ , 53(6):1915–1955, 2019.
  * Antman [1998] Stuart S. Antman.  Physically unacceptable viscous stresses.  _Zeitschrift für angewandte Mathematik und Physik_ , 49(6):980–988, 1998.
  * Ball [2002] John M Ball.  Some open problems in elasticity.  In _Geometry, mechanics, and dynamics_ , pages 3–59. Springer, 2002\.  doi: doi.org/10.1007.
  * Ball and Mizel [1987] John M Ball and Victor J Mizel.  One-dimensional variational problems whose minimizers do not satisfy the euler-lagrange equation.  In _Analysis and thermomechanics_ , pages 285–348. Springer, 1987\.
  * Benešová et al. [2020] Barbora Benešová, Malte Kampschulte, and Sebastian Schwarzacher.  A variational approach to hyperbolic evolutions and fluid-structure interactions.  _arXiv:2008.04796 [math]_ , August 2020.
  * Biot [1941] Maurice A Biot.  General theory of three-dimensional consolidation.  _Journal of applied physics_ , 12(2):155–164, 1941.
  * Bociu et al. [2020] Lorena Bociu, Sunčica Čanić, Boris Muha, and Justin T Webster.  Multilayered poroelasticity interacting with stokes flow.  _arXiv preprint arXiv:2011.12602_ , 2020.
  * Breit et al. [2021] Dominic Breit, Malte Kampschulte, and Sebastian Schwarzacher.  Compressible fluids interacting with 3D visco-elastic bulk solids.  _arXiv:2108.03042 [math]_ , August 2021.
  * Brinkman [1949] Hendrik C Brinkman.  A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles.  _Flow, Turbulence and Combustion_ , 1(1):27–34, 1949.
  * Cesmelioglu [2017] Aycil Cesmelioglu.  Analysis of the coupled navier–stokes/biot problem.  _Journal of Mathematical Analysis and Applications_ , 456(2):970–991, 2017.
  * Chapelle and Moireau [2014] Dominique Chapelle and Philippe Moireau.  General coupling of porous flows and hyperelastic formulations—from thermodynamics principles to energy balance and compatible time schemes.  _European Journal of Mechanics-B/Fluids_ , 46:82–96, 2014\.
  * Ciarlet [2000] Philippe G. Ciarlet.  _Mathematical elasticity. Vol. III_ , volume 29 of _Studies in Mathematics and its Applications_.  North-Holland Publishing Co., Amsterdam, 2000.  ISBN 0-444-82891-5.  Theory of shells.
  * Ciarlet and Nečas [1987] Philippe G. Ciarlet and Jindřich Nečas.  Injectivity and self-contact in nonlinear elasticity.  _Arch. Rat. Mech. Anal._ , 97(3):171–188, Sep 1987.  doi: 10.1007/BF00250807.
  * Conca [1988] Carlos Conca.  The stokes sieve problem.  _Communications in applied numerical methods_ , 4(1):113–121, 1988.
  * Coussy [2004] Olivier Coussy.  _Poromechanics_.  John Wiley & Sons, 2004.
  * Darcy [1856] Henry Darcy.  _Les fontaines publiques de la ville de Dijon: exposition et application…_ Victor Dalmont, 1856.
  * De Boer [2006] Reint De Boer.  _Trends in continuum mechanics of porous media_ , volume 18.  Springer Science & Business Media, 2006.
  * De Giorgi [1993] Ennio De Giorgi.  New problems on minimizing movements.  _Ennio de Giorgi: Selected Papers_ , pages 699–713, 1993.
  * Dell’Isola et al. [2009] Francesco Dell’Isola, Angela Madeo, and Pierre Seppecher.  Boundary conditions at fluid-permeable interfaces in porous media: A variational approach.  _International Journal of Solids and Structures_ , 46(17):3150–3164, 2009.
  * Fried and Gurtin [2006] Eliot Fried and Morton E Gurtin.  Tractions, balances, and boundary conditions for nonsimple materials with application to liquid flow at small-length scales.  _Archive for Rational Mechanics and Analysis_ , 182(3):513–554, 2006.
  * Halphen and Son Nguyen [1975] Bernard Halphen and Quoc Son Nguyen.  Sur les matériaux standard généralisés.  _Journal de Mécanique_ , 14:39–63, 1975.
  * Healey and Krömer [2009] Timothy J. Healey and Stefan Krömer.  Injective weak solutions in second-gradient nonlinear elasticity.  _ESAIM: Control, Optimisation and Calculus of Variations_ , 15(4):863–871, 2009.
  * Hornung [1996] Ulrich Hornung.  _Homogenization and porous media_ , volume 6.  Springer Science & Business Media, 1996.
  * Kampschulte et al. [2020] Malte Kampschulte, Sebastian Schwarzacher, and Gianmarco Sperone.  Unrestricted deformations of thin elastic structures interacting with fluids.  _in preperation_ , 2020.
  * Kružík and Roubíček [2019] Martin Kružík and Tomáš Roubíček.  _Mathematical methods in continuum mechanics of solids_.  Springer, 2019.
  * Ogden and Holzapfel [2006] Ray W Ogden and Gerhard A Holzapfel.  _Mechanics of biological tissue_.  Springer, 2006.
  * Podio-Guidugli and Vianello [2010] Paolo Podio-Guidugli and Maurizio Vianello.  Hypertractions and hyperstresses convey the same mechanical information.  _Continuum Mechanics and Thermodynamics_ , 22(3):163–176, 2010.
  * Rajagopal [2007] Kumbakonam R Rajagopal.  On a hierarchy of approximate models for flows of incompressible fluids through porous solids.  _Mathematical Models and Methods in Applied Sciences_ , 17(02):215–252, 2007.
  * Roubíček [2005] Tomáš Roubíček.  _Nonlinear partial differential equations with applications_ , volume 153 of _International Series of Numerical Mathematics_.  Birkhäuser Verlag, Basel, 2005.  ISBN 978-3-7643-7293-4; 3-7643-7293-1.
  * Showalter [2005] Ralph E Showalter.  Poroelastic filtration coupled to stokes flow.  In _Control theory of partial differential equations_ , pages 243–256. Chapman and Hall/CRC, 2005.
  * Šilhavỳ [1985] Miroslav Šilhavỳ.  Phase transitions in non-simple bodies.  _Archive for rational mechanics and analysis_ , 88(2):135–161, 1985.
  * Toupin [1962] R Toupin.  Elastic materials with couple-stresses.  _Archive for rational mechanics and analysis_ , 11(1):385–414, 1962.
  * Wilbrandt [2019] Ulrich Wilbrandt.  _Stokes–Darcy Equations: Analytic and Numerical Analysis_.  Springer, 2019.
  * Ziegler and Wehrli [1987a] H. Ziegler and C. Wehrli.  On a principle of maximal rate of entropy production.  _J. Non-Equilib. Thermodyn._ , 12(3):229–243, 1987a.  doi: 10.1515/jnet.1987.12.3.229.
  * Ziegler [1963] Hans Ziegler.  Some extremum principles in irreversible thermodynamics with application to continuum mechanics.  In _Progress in Solid Mechanics, Vol. IV_ , pages 91–193. North-Holland, Amsterdam, 1963.
  * Ziegler and Wehrli [1987b] Hans Ziegler and Christoph Wehrli.  The derivation of constitutive relations from the free energy and the dissipation function.  _Adv. Appl. Mech._ , 25:183–238, 1987b.  doi: 10.1016/S0065-2156(08)70278-3.