Compressible fluids interacting with 3D visco-elastic bulk solids

We consider the physical setup of a three-dimensional fluid-structure interaction problem. A viscous compressible gas or liquid interacts with a nonlinear, visco-elastic, three-dimensional bulk solid. The latter is described by a hyperbolic evolution with a non-convex elastic energy functional. The fluid is modelled by the compressible Navier--Stokes equations with a barotropic pressure law. Due to the motion of the solid, the fluid domain is time-changing. Our main result is the long-time existence of a weak solution to the coupled system until the time of a collision. The nonlinear coupling between the motions of the two different matters is established via the method of minimising movements. The motion of both the solid and the fluid is chosen via an incrimental minimization with respect to dissipative and static potentials. These variational choices together with a careful construction of an underlying flow map for our approximation then directly result in the pressure gradient and the material time derivatives.


Introduction
The interaction of fluids with elastic strucures has attracted the interest of scientists from various fields due to the immense potential for applications ranging from hydro-and aero-elasticity [Dow15] over bio-mechanics [BGN14] to hydrodynamics [Cha15].In last two decades there has also been a huge development in the understanding of the mathematical models -typically highly coupled systems of nonlinear partial differential equations (PDEs).The interaction of incompressible viscous fluids with elastic structures has been studied intensively.Results concerning the existence of weak solutions to the coupled system (which exist as long as the moving part of the structure does not touch the fixed part of the fluid boundary) can be found, for instance, in [CDEG05, Gra08, LR14, MC13, Mv14, MS19].The first result in [CDEG05] is concerned with a flexible elastic plate located on one part of the fluid boundary.The shell equation is linearised and the shell is assumed to be one-dimensional.The existence of a weak solution to the incompressible Navier-Stokes equations coupled with a plate in flexion is proved in [Gra08].In [MC13] the incompressible Navier-Stokes equations are studied in a cylindrical wall.The movement of the latter is modelled by the one-dimensional cylindrical linearised Koiter shell model.The elastodynamics of the cylinder wall in [Mv14] is governed by the one-dimensional linear wave equation modelling a thin structural layer, and by the two-dimensional equations of linear elasticity modelling a thick structural layer.The interaction with a linear-elastic shell of Koiter-type in a general geometric set-up (where the middle surface of the shell serves as the mathematical boundary of the threedimensional fluid domain) is studied in [LR14].The result has recently been extended to original (and fully nonlinear) Koiter model, cf.[MS19].What all these results have in common is that the incompressible Navier-Stokes equations are coupled to a lower dimensional hyperbolic equation for the structure.We are interested in models where fluid domain and structure have the same dimension which is the case for the interaction with an elastic bulk.A major difficulty is that the associated elastic energy involves a non-convex functional.Consequently, it is not possible to apply fixed point arguments to obtain a solution and one is forced to use variational methods instead.A corresponding approach has been developed recently by the second and third author together with B. Benešová in [BKS20].They provide an approach to the fully coupled system which is based on De Giorgi's celebrated minimising movement method.The key idea in [BKS20] is that the second order time derivative of the structure displacement is discretised in two steps (yielding a velocity scale and an acceleration scale) and to use a Lagrangian approximation of the material time derivative on the time-discrete level.The result is the existence of a global-in-time weak solution to coupled system describing the interaction of an incompressible fluid with a three-dimensional visco-elastic bulk solid.The situation in the compressible case is completely different and the analysis of problems from fluid-structure interaction is still at the beginning.A first result has been achieved by the first and third author in [BS18].They prove the existence of a weak solution to a coupled system describing the motion of a three-dimensional compressible fluid interacting with a two-dimensional linear elastic shell of Koiter-type.Recently, they extended the result to heat-conducting fluids and fully nonlinear shell models in [BS21].Results on the existence of local strong solutions appeared recently in [MRT21,MT20,Mit20,ST20].In this paper we aim at the natural next step and consider the interaction of a compressible fluid with a three-dimensional visco-elastic bulk solid in order to arrive at a compressible counterpart of the result in [BKS20].
Let us present the model in detail.The fluid together with the elastic structure are both confined to a fixed container -a bounded Lipschitz domain Ω ⊂ R n with n ≥ 2 (where n = 3 is the most interesting case, but also n = 2 has physical relevance).The deformation of the solid is described by the deformation function η : Q → Ω, where Q ⊂ R n is the reference configuration of the solid, which is assumed to be a bounded Lipschitz domain as well.We denote by M the part of Q which is mapped to the contact interface between the fluid and the solid and set P := ∂Q \ M .For a given deformation, we then deal with a fluid domain Ω η := Ω \ η(Q), which will be time-dependent.In Ω η we observe the flow of a viscous compressible fluid subject to the volume force f f : I × Ω η → R n .We seek the velocity field v : I × Ω η → R n and the density ̺ : I × Ω η → R solving the system in I × M, (1.3) v(t, x) = 0 in I × ∂Ω \ η(P ), (1.4) in Ω η(0) .(1.5) Here, p(̺) is the pressure which is assumed to follow the γ-law, that is p(̺) ∼ ̺ γ for large ̺, where γ > 1 (see Section 2.1 for the precise assumptions).Further, we suppose Newton's rheological law with strictly positive viscosity coefficients µ, λ.The balance of linear momentum for the solid is given by where ̺ s > 0 is the density of the shell and f s : I × Q → R 3 a given external force.We suppose that the first Piola-Kirchhoff stress tensor σ can be derived from underlying energy and dissipation potentials; that is for some energy-and dissipation functionals E and R. Prototypical example examples are given by R(η, ∂ t η) = ˆQ |∂ t ∇η ⊤ ∇η + ∇η ⊤ ∂ t ∇η| 2 dy, and for some q > n and a > nq The latter is defined provided det ∇η > 0 a.e. in Q and we set E(η) = ∞ otherwise.Here C is the positive definite tensor of elastic constants.The general assumptions for E and R are collected in Section 2.1.Equation (1.6) is supplemented with the boundary conditions σ(t, x)ν(x) = S(∇v) − p(̺)I ν(t, η(t, x)) in I × M, η(t, x) = γ in I × P. (1.7)Here ν(x) is the unit normal to M while ν(t, η(t, x)) = cof(∇η(t, x))ν(x) is the normal transformed to the actual configuration and γ : P → Ω is a given function.Finally, we assume the initial conditions η(0, •) = η 0 , ∂ t η(0, •) = η 1 in Q, (1.8) where η 0 , η 1 : Q → R are given functions.We aim to prove the existence of a weak solution to (1.1)-(1.8),where the precise formulation can be found in Section 2.3.A simplified versions of our main result reads as follows and we refer to Theorem 2.11 for the complete statement.
Theorem 1.1.Under natural assumptions on the data there exists a weak solution (η, v, ̺) to (1.1)-(1.8)which satisfies the energy inequality dy + E(η(t)), (1.9) for all s ∈ I, where H is the pressure potential related to p by p(̺) = H ′ (̺)̺ − H(̺).The interval of existence is of the form I = (0, T ), where T is the first time of collision or ∞ if there is none.
To prove Theorem 1.1 we aim to apply a variational approach in the spirit of [BKS20], where the same problem was solved in the easier case of an incompressible fluid.This approach is based on a time-delayed approximation: One seeks a continuous solution for which the inertial terms (i.e.∂ tt η and the material derivative Dv Dt ) occur in form of a difference quotient involving a time step h > 0. The resulting equation is of gradient-flow type, which allows us to construct solutions using the minimizing movements scheme.On a formal level, this method is well suited to the compressible case: The flow map, that is constructed to obtain the material derivative, can directly be used to transport the density.It turns out that the variational nature of the scheme allows us to generate the pressure term directly from the pressure potential which we include in the functional we are minimising in each step.The material derivative of the density also occurs as a difference quotient which leads to the desired equation of continuity in the limit.
Let us now explain how to make these ideas rigorous.When solving the compressible Navier-Stokes equations it is common to work with an artificial viscosity ε in (1.1) and add ε∆̺ on the right-hand side to make it a proper parabolic equation.As it turns out the limit ε → 0 can only be performed if the integrability of the density (which results from the parameter γ in (1.2)) is large and thus outside the realm of physical interest.Consequently, a second regularisation (δ-level with artificial pressure) is used and one adds an additional term δ̺ β (with β sufficiently large) to the pressure in (1.2).This approach has been introduced in [FNP01].In order to solve the regularised problem with ε, δ > 0 fixed, we aim at a variational approach as described above.Due to the additional term ε∆̺, we need to modify the mass transport from a simple update of the densities via the flow map.This change also necessitates a rather specific choice of the discretised inertial term in order to obtain the correct energy inequality (see Section 4.1 for details).Additionally, when combined with the pressure potential, this perturbation creates in the limit of the time-step h → 0 the term in the momentum equation (as well as a similar term in the energy inequality).This corresponds to a regularisation used similarly in [FNP01] and many subsequent papers.In order to recover this term, it turns out that we have to exclude the vacuum for the limit passage h → 0, see equation (4.29).Excluding the hypothetical vacuum is a big open problem for the compressible Navier-Stokes system even in presence of positive ε.To overcome this difficulty, we split the regularisation of v, that was already present in [BKS20], on the h-level into its own κ-level: For fixed κ > 0 we add a higher order dissipation to the momentum equation.This gives sufficient regularity of the velocity which ultimately yields a minumum principle for the equation of continuity.
For fixed κ, ε, δ > 0 we are able to apply the ideas just explained and to obtain a solution to the approximate problem by the minimising movement method (with the acceleration scale limit τ → 0 in Section 3.3 and the velocity scale limit h → 0 in Section 4).Eventually, we pass to the limit with respect to the regularisation parameters κ, ε and δ in Section 5.For technical reasons this has to be done in three independent steps.The limit κ → 0 is rather straightforward as the density remains compact for ε > 0. The compactness of the density becomes critical in the subsequence limit procedures where we pass to the limit in ε and δ respectively.This is done via the method of effective viscous flux which is due to Lions [Lio98] (with important extensions by Feireisl et al. [FNP01]).This method has been extended to the setting of variable domains in [BS18].
It is important to note that in each of these approximations, our approximate solutions are already confined to the set of admissible states.In particular, on each level the solid deformation η is injective, the solid and the fluid move with the same velocities at the interface and the total mass of the fluid is conserved at all times.Additionally, an energy inequality holds on each level.It mirrors the physical energy inequality and will be our main tool to obtain a priori estimates.

Preliminaries
2.1.Assumptions.The assumptions on the solid and its dissipation are identical to those in [BKS20]: Assumption 2.1 (Elastic energy).We assume that q > n and E : W 2,q (Q; Ω) → R satisfies: S1 Lower bound: There exists a number E min > −∞ such that S2 Lower bound in the determinant: For any E 0 > 0 there exists S5 Existence of derivatives: For finite values E has a well defined derivative which we will formally denote by Furthermore, on any sublevel-set of E, DE is bounded and continuous with respect to strong W 2,q (Q; R n )-convergence.S6 Monotonicity and Minty-type property: Definition 2.2 (Domain of definition).The set of admissible deformations in W 2,q (Q; Ω) (injective a.e. and satisfying the Dirichlet boundary condition) can be expressed as where γ : P → ∂Ω is a fixed injective function of sufficient regularity so that E is non-empty.By a slight abuse of notation, we use L ∞ (I; E) to denote the set of functions η : R2 Homogeneity of degree two: The dissipation is homogeneous of degree two in its second argument , i.e., R(η, λb) = λ 2 R(η, b) ∀λ ∈ R. In particular, this implies R(η, b) ≥ 0 and R(η, 0) = 0. R3 Energy-dependent Korn-type inequality: Fix E 0 > 0. Then there exists a constant exists and is weakly continuous in its two arguments.Due to the homogeneity of degree two this in particular implies Additionally, we need to state some assumptions on the potential energy of the density, which will also result in the pressure response.
P3: p has γ-growth with γ > 2n(n−1) 3n−2 , i.e. there exists a > 0 and γ > 2n(n−1) 3n−2 such that We will also consider the regularised pressure p δ given by p δ (̺) := p(̺) + δ̺ β + δ̺ 2 for δ > 0, where β > max{4, γ}.The function p δ clearly inherits P1 and P2 but has β-growth instead of γ-growth.Additionally, one also has that p ′′ δ (̺) ≥ 2δ for all ̺ > 0. Given the pressure through the function p the potential energy of the density can be described by the fluid potential For the pressure potential H we also obtain a regularised version given by H δ (̺) := ̺ ´̺ c p δ (z) During the minimising movement approach in Sections 3 and 4 only H δ appears.Some of its properties, which follow directly from P1-P3 above, are listed in the following lemma.
Recall that the moving part of the boundary of Ω η(t) is given by η(t)| M : M → Ω.The corresponding function spaces for variable domains are defined as follows.
Definition 2.6.(Function spaces) For I = (0, T ), T > 0, and η ∈ C(I × Q; Ω) defining a changing domain Ω(t) := Ω \ η(t, Q) we define for 1 ≤ p, r ≤ ∞ Function spaces of vector-or matrix valued functions are defined accordingly.We now give a definition of convergence in variable domains.Convergence in Lebesgue spaces follows from an extension by zero.
(a) We say that a sequence (g i ) ⊂ L p (I, L q (Ω ηi )) converges to g in L p (I, L q (Ω η )) strongly with respect to (η i ), in symbols (b) Let p, q < ∞.We say that a sequence (g i ) ⊂ L p (I, L q (Ω ηi )) converges to g in L p (I, L q (Ω η )) weakly with respect to (η i ), in symbols g i ⇀ η g in L p (I, L q (Ω η )), if (c) Let p = ∞ and q < ∞.We say that a sequence (g i ) ⊂ L ∞ (I, L q (Ω ηi )) converges to g in L ∞ (I, L q (Ω η )) weakly * with respect to (η i ), in symbols Next we state a compactness lemma from [BS18, Lemma 2.8] which gives a variant of the classical result by Aubin-Lions for PDEs in variables domains.It allows to pass to the limit in the product of two weakly converging sequences provided one is more regular in space and the other one in time.Let us list the required assumptions.(A1) The sequence (η i ) ⊂ C(I × Q; Ω) satisfies η i → η uniformly in I × Q and for some α > 0 (A2) Let (v i ) be a sequence such that for some p, s ∈ [1, ∞) we have (A3) Let (r i ) be a sequence such that for some m, b ∈ [1, ∞) we have Assume further that (∂ t r i ) is bounded in the sense of distributions, i.e., there is c > 0 and m ∈ N such that uniformly in i for all φ ∈ C ∞ c (I × Ω ηi ).In [BS18, Lemma 2.8] the corresponding version of (A3) assumes m = 2.But the very same argument is also valid in the general case as in the classical Aubin-Lions argument.
Lemma 2.8.Let (η i ), (v i ) and (r i ) be sequences satisfying (A1)-(A3) where Then there is a subsequence with Corollary 2.9.In the case r i = v i we find that 2.3.Weak solutions and the main theorem.In this session we make the concept of a weak solution to (1.1)-(1.8)rigorous.We begin with the function spaces to which the triple (η, v, ̺) belongs.
• For the solid deformation η : I × Q → Ω we consider the space • Given η ∈ Y I , for the fluid velocity v : I × Ω η → R n we define the space  for all ψ ∈ C ∞ (I × R 3 ) and we have ̺(0) = ̺ 0 .
• The energy inequality is satisfied in the sense that Here, we abbreviated Remark 2.10.Due to the bulk setting, it is possible to extend the fluid variables onto the fixed domain Ω using their solid counterparts.We will do so quite often for the velocity, where it is convenient to set v(t, η(t, y)) := ∂ t η(t, y) for all y ∈ Q, as this results in the complete Eulerian velocity v ∈ L 2 (I; W 1,2 0 (Ω; R n )).For the density, one can similarly consider pushing the solid density forward onto Ω using η.However in practice this is less useful, as the resulting function will still have a jump at the interface between solid and fluid.Instead, it is often more convenient to keep the inertial effects of the solid in their Lagrangian description and to think of ̺ as extended by 0 onto Ω, as this allows for the removal of most of the time-dependent domains in the weak solution.
We are finally in the position to state our main result in a complete form.
Theorem 2.11.Assume that we have There is a weak solution (η, v, ̺) ∈ Y I × X I η × Z I η to (1.1)-(1.8) in the sense of (2.3)-(2.5).Here, we have I = (0, T * ), where T * < T only if the time T * is the time of the first contact of the free boundary of the solid body either with itself or ∂Ω (i.e.η(T * ) ∈ ∂E).
As will be apparent by the analysis we will show that the renormalised continuity equation in the sense of DiPerna and Lions is satisfied, cf.[DL89,Lio98].
2.4.The damped continuity equation in variable domains.Here we present some results concerning the continuity equation in moving domains.Some of the results are direct consequences of our previous papers [BS18,BS21], but Lemma 2.14 below is new and contains some improvements.We assume that the moving boundary is prescribed by a function η : I ×Q → Ω.For a given function in Q for a.a.t ∈ I and some ε > 0 we consider the equation a) There is at most one weak solution ̺ to (2.7) such that for large values of s and θ(0) = 0.Then the following holds: Proof.The proof follows by the lines of [BS18, Thm.3.1], where much stronger conditions on the regularity of η and also sightly more on w is assumed.One easily checks that these assumptions are only used there to prove the existence of a solution, which we do not claim here.The statements from a)-c) do not require it.Note that the assumption Lemma 2.14.Let the assumptions of Lemma 2.13 be satisfied and suppose additionally that w ∈ L 2 (I; Proof.Let us initially suppose that η, w and ̺ are sufficiently smooth and ̺ is stricly positive such that all the following manipulations are justified.Ad a).Multiplying (2.7) by θ ′ (̺) shows where Similarly, choosing θ(z) = z −m , we have Taking the m-th root shows Ad (b).Multiplying (2.7) by ∆̺ shows The first and last terms are estimated by , where ξ > 0 is arbitrary and where we used the trace theorem which requires together with an interpolation argument.We also have . Absorbing the ξ-terms and applying Gronwall's lemma proves where the constant depends on Let us now remove the regularity assumptions on η, w and ̺.We can regularise η and w by smooth approximation as follows.First, we extend w by ∂ t η • η −1 to Ω and regularise w by a standard smooth approximation in space-time.This yields a smooth sequence (w) ξ which converges to w in the , where (η 0 ) ξ is a regularisation of η 0 in space.The function (η) ξ (•, x) does indeed exist for all given x ∈ Q on the interval [0, T ] by the Picard-Lindelöff theorem as ∇w ξ ∞ is uniformly bounded (in dependence of ξ).By its equation one directly deduces that (η) ξ is smooth.Now [BS21, Theorem 3.3] applies to the regularised problem and we obtain a solution ̺ ξ with ̺ ξ ∈ C 1 (I × Ω η ) and ∇ 2 ̺ ξ ∈ C(I × Ω η ).Also ̺ ξ is stricly positive.One easily checks that the estimates derived above do not depend on ξ.It remains to pass to the limit ξ → 0. Since η 0 ∈ W 2,q (Q; R n ) by assumption and q > n we have (η 0 ) ξ → η 0 in W 1,∞ (Q; R n ).Hence, (using the ODE, the properties of (w) ξ and Gronwall's lemma) one deduces that (η) ξ → η as ξ → 0 uniformly in Q.Actually, for this purpose convergence of w in the L 2 (I; W 1,∞ (Ω; R n )) is sufficient.Since we have convergence in L 2 (I; W 2,∞ (Ω; R n )) we have similarly ∇(η) ξ → ∇η as ξ → 0 uniformly in Q.Also, passing to the limit with the ODE implies that These convergences are sufficient to pass to the limit in the estimate.Since (2.7) is linear the limit procedure ξ → 0 in (2.8) is straightforward and the limit is indeed the unique solution (see Lemma 2.13 a)).

The time delayed problem
Following [BKS20] we begin constructing a solution to a time-delayed problem, where the material derivative in the momentum equation and its solid counterpart are discretised at level h > 0, but all other variables are already continuous.We initially solve only on the interval [0, h] and then iterate this in the next section to a time-delayed solution on I by decomposing it into intervals [0, h], [h, 2h], . . . .As it is common in the literature on the compressible Navier-Stokes system we use an artificial diffusion in the continuity equation and an approximate pressure.As it turns out this alone does not yield sufficient regularity to pass to the limit in the time-step h, which we are going to do in the next section.In order to overcome this problem we add additional terms and set for κ > 0 and k 0 ∈ N large enough In the model for the bulk we replace E by E κ and R by R κ respectively.Additionally, we add κ ´Ω(t) |∇ k0 v| 2 dx to the dissipation of the fluid.Similar terms are also used in [BKS20] with h instead of κ.Hence there they disappear already the limit h → 0, simultaneously with turning the differential quotient into the material derivative.In our case this is split into two limit procedures.
• The time delayed momentum equation holds, that is we have • The approximate equation of continuity holds, that is we have 2 Due to the dissipation of density in the equation of continuity, we cannot simply work with the previous fluid velocity.Additionally, we need to deal with the transport effects inherent in the Eulerian description of the fluid.As a result, the natural quantity w(t) turns out to be ̺(t − h)v(t − h) transported forward along the flow to the domain Ω 0 at time t = 0.This will become more apparent in the subsequent section, when we extend the problem to [0, T ] and perform the limit h → 0. For the purpose of the current section, all terms involving w can be thought of as a given force term.
Theorem 3.1.Suppose that there are Then there is a triple which solves the time delayed problem in the sense of (3.1)-(3.3).The rest of this section is devoted to the proof of Theorem 3.1.Some aspects are reminiscent to Theorem 4.2 and its predecessors Theorems 2.2 (for the FSI) and 3.5 (for the time delayed terms) in [BKS20] to which we refer when possible.The main effort is to understand to contribution of ̺ and its behaviour.In order to construct a solution we now splitting [0, h] again into small time-steps of length τ ≪ h and discretise in time.For each of these steps we solve a stationary minimisation problem (the iterative problem), prove a discrete version of the energy inequality and finally pass to the limit τ → 0 with the resulting piecewise constant and affine approximations.
3.1.The iterative approximation.Assume that τ, h, κ, ε, δ > 0 are fixed and that the forces f f and f s as well as where we require For the fluid, we require a regularised condition for mass transport, that is ̺ and v are related through Here and in the future Ψ v := id + τ v is a helpful shorthand and (id − τ ε∆) −1 is to be understood as the solution operator to the respective Neumann problem on Ω k , i.e. (id Clearly, there is a unique solution ̺k ∈ W 2,β (Ω k ) since ̺ k can be assumed to be in L β (Ω k ) and Ω k is a C 1,α -domain. 3Furthermore, due to the fact that v uniquely determines ̺ through (3.6) we can rewrite (3.4) as a minimisation problem in (η, v) only.Specifically, we can write Finally, we update Φ k to Φ k+1 by setting Lemma 3.2.Suppose that η k and ̺ k are given, where Then the minimisation problem (3.4)-(3.6)has a solution (η k+1 , v k+1 , ̺ k+1 ), where Proof.We argue by the direct method in the calculus of variation.The functional is clearly welldefined by the choice of the function spaces.Using Young's inequality to estimate the force terms we can also show that it is bounded from below in each term and inserting (η k , 0, ̺k ) as a candidate shows that the minimiser must have a finite value.Using (3.8) we can rewrite the problem as a minimisation in (η, v) only.Coercivity in the relevant functional spaces is now obvious, recalling assumption 2.1 S4.However, we still need to verify that any limit obeys the lower bounds.
A standard application of the minimum principle 4 to (3.7) implies that inf is part of the functional, we know that ∇v C α (Ω k ) is uniformly bounded along any minimising sequence.Thus det(I + τ ∇v) is trivially bounded from above, which gives us inf Ω k ̺ > 0 for any ̺ of finite energy in the functional.
Additionally, by the growth condition H3 of Lemma 2.5 and (3.8) we find the following bound on the determinant in dependence of the lower bound of ̺k .Combining this with the C α bound on ∇v and choosing β sufficiently large, we get a nonzero lower bound on det(I + τ ∇v) by [BKS20, Prop.2.24 (S2)], which ultimately goes back to [HK09].Note that together with the uniform C 1,α -bounds on η, this also implies that id +τ v is a diffeomorphism up to the boundary of the fluid domain, which in turn implies that there cannot be a collision.The final point which needs clarification is lower semi-continuity, which for most terms does not differ from the analysis in the incompressible case and for these we refer to [BKS20, Propositions 2.13 and 4.3].For the new term Ũ δ,ε ̺ k however, weak convergence in W k0,2 (Ω k ) implies strong convergence in C 1 (Ω k ).Since the functional Ũ δ,ε ̺ k is continuous on C 1 (Ω) (as det(I + τ ∇v) is strictly positive for all velocities for which the functional is of finite value) the proof is complete.
Next we calculate the corresponding Euler-Lagrange equation, which will give us the discrete, time-delayed momentum-balance.Assuming that (η k+1 , v k+1 ) is in the interior of the admissible set 3 In fact, we can infer from Lemma 3.2 that ̺k ∈ W 3,β (Ω k ), but note that this cannot be iterated for increasing k, as the Ψv k+1 appearing in the definition of ̺ k+1 presents an upper limit to regularity. 4 Here we use the assumption that ρ k is strictly bounded from above and below.Due to the Neumann boundary conditions, we find that ρk satisfies the same upper and lower bounds as ρ k .
(i.e. the solid has no collisions and det The first integral is the discrete version of the pressure term and the second is a dissipation related error that will be shown to vanish in the limit τ → 0. Furthermore we note that a short calculation reveals

The discrete energy inequality.
A key in the analysis is the energy inequality.We start with a discrete energy inequality for the solution of the minimisation problem (3.4)-(3.6). 5 Lemma 3.4.Suppose that η 0 , v 0 and ̺ 0 , as well as w, ζ, f f and f s are given, where Note that this inequality is non-optimal; one would expect an additional factor of two in front of all dissipation terms.Using a more involved argument, one can indeed derive an improved inequality at this point.But since there is no qualitative improvement on the resulting estimates and this inequality will be replaced by a more correct, continuous energy inequality afterwards, there is no need to do so here.
Proof.We compare the value of the actual minimizer in (3.4) with the value at (η, v, ̺) = (η k , 0, ̺k ), where ̺k := (id − τ ε∆) −1 ̺ k .Thus we get (3.10) We are now going to estimate the error between Next we construct piecewise constant and piecewise continuous interpolations of all our quantities by setting as well as Ω (τ ) (t) = Ω k for τ k ≤ t < τ (k + 1).Note that from this point on, we extend v to all of Ω, using the corresponding solid velocity, cf.Remark 2.10.Lemma 3.4 implies the following corollary for the interpolated quantities by setting t = τ N : Corollary 3.5.Under the assumptions of Lemma 3.4 we have Absorbing the forcing terms into the left-hand side by means of Young's inequality we obtain the following estimates They are uniform in τ but may depend on the other parameters h, κ, ε and δ.Note also that (3.12) implies for some α > 0 using Sobolev's embedding.
Let us finally deduce uniform bounds for det ∇Φ k .Note the eventhough the general idea stays the same as in [BKS20, Prop.4.6], the estimate is slightly different since we cannot use that div v = 0 in the compressible regime.
Lemma 3.6 (Bounds on Φ).There exists constants c, C > 0 such that for all small enough h and τ as well as where M l are the polynomials of order l stemming from the expansion of the determinant.Estimating the arithmetric mean by the geometric mean yields In order to get a similar bound form below we use (1 + a) −1 ≤ 1 + 2|a| for |a| ≪ 1 as well as uniform boundedness of ) , which follows from restricting (3.14) to a single term.We then can infer similarly to the above (3.17) In total we arrive at (3.15).Arguing in the same way, we can show uniformly in k, τ and h using also (3.12).
Similarly to the above we can also control second order derivative of Φ k .It holds such that, using (3.16), uniformly in k by (3.12) provided we choose k 0 large enough.Again c is independent of τ and h.
We only have to prove that the terms involving the density converge to their correct counterparts, that is as τ → 0. The limit in the remaining terms can be performed as in [BKS20, Section 4.1].The convergence in (3.32) follows directly from (3.22) and (3.24), whereas (3.30) and (3.31) require strong convergence of the density.Using (3.6) and (3.7) we can write for all ψ ∈ W 1,2 (Ω k ).Now we choose a parabolic cylinder J × B such that 2B ⋐ Ω η (τ ) for all τ small enough.This is possible due to (3.24).We obtain for ψ Here the quantity o(τ ) is such that o(τ )/τ vanishes in the L 2 (J; L ∞ (B))-norm as a consequence of (3.21).Using (3.6) we may deduce from the regularity of v (τ ) in (3.11)-(3.13)as well as the uniform lower bound for det ∇Φ k from (3.15)) that uniformly in τ .Combining this with (3.23) we conclude that and This together with (3.29) yields Using (3.13) and arbitrariness of J × B, the convergences in (3.34) and (3.35) even hold in L q ((0, h)× Ω (τ ) ) with q = 4+3β 3 > β.This in combination with (3.24) and (3.29) is enough to prove the convergenc (3.30) Taking into account also (3.21), (3.26) and (3.27) proves (3.31) (note also the uniform lower bound for det ∇Φ k from (3.15)).We thus conclude that (3.1) holds.Now we take a look at the continuity equation.We use a test-function ψ ∈ C ∞ c ((τ, h − τ ) × Ω) and obtain similarly to (3.33) which we denote by I = II + III.The term I on the left-hand side can be rewritten as with o(τ ) as above.Due to (3.26), (3.34), the smoothness of ψ and (3.25) we find On the other hand, for all ψ ∈ C ∞ c ((0, h) × Ω).We have shown (3.2).Finally, for the energy inequality (3.3), we note that we cannot simply pass to the limit in Lemma 3.4, as this is missing a factor 2 in front of the dissipation terms.Instead, we obtain it as in [BKS20, Lemma 4.8] by testing (3.1) with (∂ t η, v).The only substantial addition here is the pressure term which reads as 2) and Reynold's transport theorem (due to Lemma 2.14 (b) the density is smooth enough to rigorously perform these computations).For the inertial term we now have by Young's inequality and a change of variables.This also shows the need for the factor √ det ∇Φ −1 in the equations.

Inertial problem (h → 0)
The aim of the present section is to pass to the limit h → 0 in (3.1)-(3.3)and to obtain a solution to the regularised system (where κ, ε and δ are fixed).We beginn with a definition of the latter.
We introduce the function spaces ), which replace the spaces Y I , X I η and Z I η (defined in Section 2.3) on the κ and ε-level respectively.A weak solution to the regularised system is a triple (η, v, ̺) ∈ Y I k0 × X I η,k0 × ẐI η,ε that satisfies the following.
• The momentum equation holds in the sense that • The continuity equation holds in the sense that holds in I × Ω η and we have ̺(0) = ̺ 0 as well as ∂ νη ̺ = 0.
• The energy inequality is satisfied in the sense that Here, we abbreviated Theorem 4.1.Assume that we have for some α ∈ (0, 1) Furthermore suppose that ̺ 0 is strictly positive.Then there is a solution Here, we have I = (0, T ), where T ∈ (0, ∞) is arbitrary.

4.1.
A priori analysis.We proceed as follows.First we need to produce an h-approximation on the whole interval I.For h ≪ 1 we decompose the interval I into subintervals (0, h), (h, 2h), . . . .For any given h we obtain from Theorem 3.1 the existence of a solution to (3.1)-(3.3) in (0, h).We will then use the resulting η(h), ̺(h), ∂ t η and a suitably modified version of v • Φ(t) • Φ(h) −1 as η 0 , ̺ 0 ζ and w on (h, 2h).We can prove that these are valid initial data, because of the energy inequality (3.3).We repeat this procedure on the following time-intervals to get a global solution (η h , v h , ̺ h ) which solves a variant of (3.1) on I. Specifically, given a solution (η 3) on the interval [0, h] (note that later this will correspond to [(l − 1)h, lh] in I, but for now, we will consider each of these intervals as [0, h] and distinguish them by the index l), we apply Theorem 3.1 with for all t ∈ [0, h] and as w(t) for all t ∈ [0, h] (defined on Ω \ η . We also write as usual Ω (h) These assignements are obvious, except for w.Here one should keep in mind that w(t) needs to be defined on the new initial fluid domain Ω 0 , but needs to correspond to the quantity √ ̺v for the matching fluid particle at an earlier time.We thus employ the flow map to move it there.As this flow map was only defined with respect to the previous reference configuration, we take a slight detour there, resulting in Φ(t) • Φ(h) −1 and finally we need to correct the distortion due to the extra term ε∆̺ in the continuity equation such that, in particular, From the energy inequality (3.3) applied to the previous solution, we can conclude that our new initial data fulfills the conditions for Theorem 3.1.The solutions constructed by this theorem will then be denoted by (η With these solutions in hand, we can now construct an h-approximation on the whole interval I. Specifically, we set as well as a redefined flow map for t ∈ [(l − 1)h, lh] and so on.Here the new flow map needs to be read in the style of semi-groups, i.e.Φ (h) s (t) maps the fluid domain at time t to the fluid domain at time t + s, in such a way that we follow the fluid flow.Note that in concert with the uniform (only κ-dependent) bounds on det ∇Φ(t) derived in the last section in (3.15) and the C 1 -regularity of η, this also implies that Φ(t) is a diffeomorphism up to the boundary and thus there cannot be a collision.
If we now translate (3.1)-(3.3)into this new notation, then we have that the h-approximation fulfills the following: • The momentum equation (3.1) translates to • The continuity equation (3.2) holds unchanged, i.e.
• the energy balance holds in the sense that for a.a.t 1 ∈ I.Most of this is a straigthforward replacement of the new definitions together with a telescope argument for the energy balances.Note that in the intertial term of the momentum-equation, the different flow maps from the equation and the definition of w combine to Φ h , as do their Jacobian determinants.
Additionally, we recover the expected properties of the flow map: Corollary 4.2.Let the assumptions of Theorem 4.1 be valid.For any t, t Proof.The first set of assertions follows directly from the definition and the properties of the shorttime solutions.Additionally, we can calculate s (t) using the contunuity equation.Combining the two relations above results in the final assertion.
From (4.6) we obtain the following uniform bounds: which are uniform in h.As a consequence we have the following convergences for some α ∈ (0, 1) after taking a (non-relabelled) subsequence.Taking further Lemma 2.14 into account we have Using this together with (4.14) in the equation of continuity we also have Moreover, ̺ (h) stays bounded and strictly positive, that is for some ̺, ̺ > 0 which do not depend on h.Finally, passing to the limit in (3.15)-(3.16)and using (4.14) yields uniform bounds for Φ (h) .In particular, we have for some c, c > 0 independent of h and s (however diverging with T → ∞ or κ → 0) as well as for a.a.s as h → 0.

4.2.
Strong convergence.To pass to the limit in (4.4)-(4.6),we need to improve some of the convergence from the last subsection to strong convergences.First we will argue similarly to [BKS20, Lem.4.14, Prop.4.15] and use (4.4) to estimate Lemma 4.3 (W −m,2 -estimates).There exists a constant C > 0 independent of h such that for an m ∈ N large enough ≤ C, (4.23) Proof.We need the following key estimate to correct the flow map.For any Lipschitz continuous function a : I × Ω → R and any s ∈ [0, h] we have where we used that ´Ω(t+r) ̺ (h) (t) dx is constant (this follows immediately from (3.2)) and the fact that ffl t t−h ´Ω(t) ̺ (h) (t + r) v (h) 2 dx dr is uniformly bounded by the energy estimate.By Lip t , Lip x and Lip t,x we denote the Lipschitz constants with respect to time, space and space-time respectively.This estimate varies from the corresponding estimate in [BKS20] by the necessary inclusion of the density, as the flow is no longer volume preserving.Now we consider (with ̺ (h) per convention extended by 0 to all of Ω) where the first term can be estimated using equation (4.4).For the second one we perform a change of variables with Φ h (t − h) in its first term to obtain Now the absolute value of the first integral can be estimated by Finally we need to prove strong convergence of the density.For this purpose we choose a parabolic cylinder J × B such that 2B ⋐ Ω (h) for all t ∈ J and all h small enough.This is possible due to (4.13).From the continuity equation (4.5) we obtain uniformly in h.This, in combination with the gradient estimate from (4.8), yields Using (4.20) and the arbitrariness of J × B we obtain for some p > β.Using Theorem 2.13 (b) with θ(s) = s 2 (which is admissible by approximation) we obtain Now We apply Theorem 2.13 (b) to the limit version of the continuity equation (that is, equation (4.2)) which results in a counterpart of (4.27).On account of (4.26) all terms converge except for ´t 0 ´Ωη (h) 2ε|∇̺ (h) | 2 dx ds, which yields convergence of the latter.Hence we obtain ). (4.28) 4.3.Derivation of the material derivative.We now take the inertial term from (4.4) and shift its second half in time by a change of variables.This gives us (note the symmetry of the √ ̺-terms) Using the fundamental theorem of calculus and Corollary 4.2 the integrand can be rewritten as6 using in particular Corollary 4.2.
We can now deal with each of these terms one after the other.Using (Φ As h → 0 we have (IV) h → 0 on account of (4.18), (4.20) and For the first term we have a W −m,2 -estimate given in (4.24).The second term can be estimated in L 1 (I; L 2 (Ω (h) )) by (4.18), (4.19), (4.20) and (4.21).We conclude that uniformly in h.Now, Lemma 2.8 yields  For this purpose we decompose Here [•] ξ denotes a regularisation in space defined by an extension to the whole space 7 and a mollification with a smooth kernel.In the above we use [. . .] as a shorthand for For fixed ξ we can use smoothness of [•] ξ to conclude that the first term vanishes as h → 0. This is a consequence of (4.25) and the a priori estimates.Also, the second term converges to zero as ξ → 0 (uniformly with respect to h, recall (4.17), (4.18) and (4.20)) by standard properties of the mollification.The last term converges to the expected limit as we have seen in (4.30).In conclusion, we have as h → 0, which finishes the proof of (4.1).

Removal of the remaining approximation parameters
In this section we pass to the limit in the approximate equations.For technical reasons the limits κ → 0, ε → 0 and δ → 0 have to be performed independently from each other.The limit κ → 0 is rather straightforward as the density remains compact for ε > 0. For the greater part of this section we study the limit ε → 0 and only highlight the difference in the δ-limit.5.1.The limit system for κ → 0. Recalling the definition of the function spaces from Section 2.3 we seek a triple (η, v, ̺) ∈ Y I × X I η × ẐI η that satisfies the following.
• The continuity equation holds in the sense that Theorem 5.1.Assume that we have for some α ∈ (0, 1) Furthermore suppose that ̺ 0 is strictly positive.There is a solution (η, v, ̺) ∈ Y I × X I η × ẐI η to (5.1)-(5.3).Here, we have I = (0, T * ), where T * < T only if the time T * is the time of the first contact of the free boundary of the solid body either with itself or ∂Ω (i.e.η(T * ) ∈ ∂E ).
Passing to a subsequence we obtain for some α ∈ (0, 1) η (κ) ⇀ * η in L ∞ (I; W 2,q (Q; Ω)), (5.9) (5.17) (5.18) Clearly, the κ-terms in (4.1) vanishes as κ → 0 as a consequence of (5.13), (5.14) and (5.16).Arguing as in [BKS20, Prop.2.23] one can use assumption S6 to deduce η (κ) → η in L q (I; W 2,q (Q; Ω)) (5.19) such that for a.e.t ∈ I DE(η (κ) (t)) → DE(η(t)) in W −2,q (Q; Ω). (5.20) In order to pass to the limit in various terms in the equations we are concerned with the compactness of ̺ (κ) .Due to (5.18) and (4.2) we can apply Corollary 2.9 to conclude (5.21) In combination with (5.7) this can be improved to for some p > β.It is easy to see that (5.22) allows to pass to the limit in all nonlinear terms of (4.1) and (4.2) except of Due to (5.15) we obtain the expected limit as κ → 0 provided we have This can be proved exactly as in (4.28) using Theorem 2.13 (b).Finally we complete proof of Theorem (5.1), by extending the solution to long times: Assume that I = [0, T ) is a maximal interval of existence with T < ∞.Then using the energy-inequality, we conclude existence of limits η(T ), ∂ t η(T ), ̺(T ) and v(T ).Now η(T ) has to have a collision, which proves the theorem.otherwise we could construct an extended solution by applying the theorem with these as initial data, which would be a contradiction.5.2.The limit system for ε → 0. We wish to establish the existence of a weak solution (η, v, ̺) to the system with artificial pressure in the following sense: We define )) as the function space for the density, whereas the other function spaces are defined in Section 2.3.A weak solution is a triple (η, v, ̺) ∈ Y I × X I η × Z I η that satisfies the following.• The momentum equation holds in the sense that •η(t) in Q and φ(t) = 0 on P , where Ω(t) := Ω\η(t, Q).Moreover, we have (̺v)(0) = q 0 , η(0) = η 0 and ∂ t η(0) = η 1 we well as • The energy inequality is satisfied in the sense that Here, we abbreviated Theorem 5.2.Assume that we have for some α ∈ (0, 1) and s > 0 Furthermore suppose that ̺ 0 is strictly positive.There is a solution (η, v, ̺) ∈ Y I × X I η × Z I η to (5.24)-(5.26).Here, we have I = (0, T * ), where T * < T only if the time T * is the time of the first contact of the free boundary of the solid body either with itself or ∂Ω (i.e.η(T * ) ∈ ∂E ).
Lemma 5.3.Under the assumptions of Theorem 5.2 the continuity equation holds in the renormalized sense as specified in Definition 2.12.
For a given ε we obtain a solution (η ε , v ε , ̺ ε ) to (5.1)-(5.3)by Theorem 5.1.In particular, we have for any 0 ≤ t 1 ≤ T , where Ω (ε) := Ω \ η (ε) (t, Q).We deduce the bounds Finally, we deduce from the equation of continuity (5.2) (using the renormalized formulation from Theorem 2.13 (b) with θ(z) = z 2 and testing with ψ ≡ 1)) that (5.31) Note that all estimates are independent of ε.Hence, we may take a subsequence such that for some α ∈ (0, 1) we have (5.38) Arguing as in [BKS20, Prop.2.23] one can again benefit from assumption S6 to deduce We observe that the a priori estimates (5.29) imply uniform bounds of Therefore, we may apply Lemma 2.8 with the choice v i ≡ v (ε) , r i = ̺ (ε) , p = s = 2, b = β and m sufficiently large to obtain where a ∈ (1, 2β β+1 ) and q ∈ (1, We apply Lemma 2.8 once more with the choice β+1 and m sufficiently large to find that We also obtain (5.44) for all q < 6β β+6 .At this stage of the proof the pressure is only bounded in L 1 , so we have to exclude its concentrations.The standard approach only works locally, where the moving shell is not seen and we obtain the following Lemma (see [BS18, Lemma 6.3] for details).
We still have to exclude concentrations of the pressure at the boundary, which is the object of the following lemma.
Lemma 5.5.Let ξ > 0 be arbitrary.There is a measurable set A ξ ⋐ I × Ω(•) such that we have for all ε ≤ ε 0 (ξ) Proof.The proof is exactly as in [BS18, Lemma 6.4] (which is inspired by [Kuk09]) provided we know that uniformly in ε.This follows from (5.34) due to the trace theorem.
We connect Lemma 5.4 and Lemma 5.5 to obtain the following corollary.
In order to exclude concentrations of the pressure at the moving boundary we need the assumption γ > 12 7 .
Proof.The proof is exactly as in [BS18, Lemma 7.4] provided we know that for all q < 4 ∂ t η (ε) ∈ L 2 (I; L q (M )) uniformly in ε.This follows from (5.62) due to the trace theorem.
(5.73) Similarly to Corollary 5.6 we have the following.
Using (5.72) and the convergences (5.60)-(5.66)we can pass to the limit in (5.24) and (5.25) and obtain It remains to show strong convergence of ̺ (δ) .We define the L ∞ -truncation Here T is a smooth concave function on R such that T (z) = z for z ≤ 1 and T (z) = 2 for z ≥ 3. we clearly have T k (̺ (δ) ) ⇀ T 1,k in C w (I; L p (R 3 )) ∀p ∈ [1, ∞), (5.78) ) div v (δ) ⇀ T 2,k in L 2 (I × R 3 ), (5.79) for some limit functions T 1,k and T 2,k .Now we have to show that ˆI×Ω η (δ)  p δ (̺ (δ) ) − (λ + 2µ) div v (δ) T k (̺ (δ) ) dx dt −→ ˆI×Ωη p − (λ + 2µ) div v T 1,k dx dt. (5.80) For this step we are able to use the theory established in [Lio98] on a local level.As in [BS18, Subsection 7.1] one can first prove a localised version of (5.80) and then use Lemma 5.8 and Corollary 5.9 to deduce the global version.The next aim is to prove that ̺ is a renormalized solution (in the sense of Definition 2.12).In order to do so it suffices to use the continuity equation and (5.80) again on the whole space.Following line by line the arguments from [BS18, Subsection 7.2] we have ∂ t T 1,k + div T 1,k v + T 2,k = 0 (5.81) in the sense of distributions on I × R n .Note that we extended ̺ by zero to R n .The next step is to show lim sup δ→0 ˆI×Ω |T k (̺ (δ) ) − T k (̺)| q dx dt ≤ C, (5.82) where C does not depend on k and q > 2 will be specified later.The proof of (5.82) follows exactly the arguments from the classical setting with fixed boundary (see [FNP01]) using (5.80) and the uniform bounds on v (δ) (with the only exception that we do not localise).Using (5.82) and arguing as in [BS18,Sec. 7.2] we obtain the renormalised continuity equation.As in [BS18,Sec. 7.3] we can use the latter one to show strong convergence of the density and as in the end of Theorem 5.1 we then extend the existence interval until the first collision, which finishes the proof.

Compliance with Ethical Standards
Funding.This research was partly funded by: (i) Primus Research Programme of Charles University, Prague (grant number PRIMUS/19/SCI/01).(ii) The program GJ19-11707Y of the Czech national grant agency (GA ČR).
Conflict of Interest.The authors declare that they have no conflict of interest.
27) in addition to (3.26).Combining (3.26) and (3.27) shows dt consists of T /h terms uniformly bounded by a multiple of h, due to the energy-inequality and by applying (4.25) with s = h and a = b For the second integral we have by Corollary 4.2 and the various estimates ∈ E and v(t) = 0 on ∂Ω for a.a.t ∈ I. • The continuity equation holds in the sense that ˆI d dt ˆΩ(t) ̺ψ dx dt − ˆI ˆΩ(t) ̺∂ t ψ + ̺v • ∇ψ dx dt = 0 (5.25) for all ψ ∈ C ∞ (I × Ω) and we have ̺(0) = ̺ 0 .
The remaining term involving (III) h is more critical since v(h), which is based on our a-priori estimates only weakly converging, appears in a product with itself.However, we may use the discrete time derivative for the momentum.It holds • The momentum equation holds in the sense that for a.a.t ∈ I and we have ̺(0) = ̺ 0 .•Theenergy inequality is satisfied in the sense that− ˆI ∂ t ψ E δ dt + ˆI ψ 7Since ∂Ω (h) is a Lipschitz boundary uniformly in time by (4.7) standard results on the extension of Sobolev functions apply.