From quasi-incompressible to semi-compressible fluids

A new concept of semi-compressible fluids is introduced for slightly compressible visco-elastic fluids (typically rather liquids than gasses) where mass density variations are negligible in some sense, while being directly controlled by pressure which is very small in comparison with the elastic bulk modulus. The physically consistent fully Eulerian models with specific dispersion of pressure-wave speed are devised. This contrasts to the so-called quasi-incompressible fluids which are described not physically consistently and, in fact, only approximate ideally incompressible ones in the limit. After surveying and modifying models for the quasi-incompressible fluids, we eventually devise some fully convective models complying with energy conservation and capturing phenomena as pressure-wave propagation with wave-length (and possibly also pressure) dependent velocity.


1.
Introduction. It is well known that propagation of pressure waves (briefly Pwaves, also referred to as longitudinal waves or, if waves are in acoustic frequencies, as sound) cannot be modeled by ideally incompressible fluid models described by incompressible Navier-Stokes equations. On the other hands, fully compressible fluid models which counts with the mass density variations are analytically and computationally very complicated (with e.g. only weak-strong uniqueness known), which is usually not necessary in liquids. For example, water density increases only by 5 % even when so extremely compressed as in the 11 km deep Mariana Trench in the Pacific ocean.
This paper focuses on modelling of liquids where mass density variation can be neglected in a certain way or, in another way, directly controlled by pressure which is very small in comparison with elastic bulk modulus. This elastic modulus is very large (in comparison with occuring pressure) but not +∞. The goal is to allow propagation of P-waves while not to drive the model too far from the incompressible Navier-Stokes one.
In fluids, typically the rheological models are very different in the deviatoric part and in the spherical (volumetric) part of the strain tensor. The usual Newtonian ideally incompressible fluids use the linear Stokes rheology in the deviatoric part and while the spherical part is ideally rigid. The essence of viscoelastic fluids is involvement of some elasticity into purely viscous/rigid rheology. In their simplest variant, viscoelastic fluids use the Kelvin-Voigt rheology in the spherical part and/or the Maxwell in the deviatoric part. Here we will be focused on the former option.
The viscoelastic fluids with a constant viscosity which flow like a liquid but behave like an elastic solid when stretched out have been invented by D.V. Boger [13] in the context of polymer solutions, cf. [36] for a historical survey. This primarily concerns the shear part but we will use this idea for the spherical (volumetric) part where it is even more natural in most of fluids. The spherical elasticity response (i.e. compressibility) is characterized by the bulk elastic modulus K, while the shear modulus is zero in fluids so that shear waves cannot propagate (in contrast to solids). Together with the mass density , this modulus determines the speed of the P-waves as K/ . E.g. water has approximately (cf. Remarks 2 and 8 below) K = 2.15 GPa and is thus much more compressible comparing e.g. with steel which has K ∼ 100 − 200 GPa.
One of a prominent application is in geophysics when modelling propagation of seismic P-waves through fluidic parts of the Earth, i.e. outer core (composed with molten iron and nickel) and oceans (composed from water with solution of salts). In mathematical literature, viscoelastic fluids with elasticity in the spherical part (called quasi-incompressible or sometimes also quasi-compressible) are often considered only as a regularization (or a numerical "pressure stabilization") of incompressible fluids, cf. e.g. [25,26,46] or, even with twisted energetics (due to omitting Temam's force 2 (div v) v in (4) below) also e.g. in [22,23] or [49,Sec. 4.2], referring (although not explicitly) to the ignored mass-density variation or partly ignored convection, cf. the discussion in Sect. 5 below. Sometimes, quasi-incompressibility may refer to mixtures and then it means that density is a function of concentration only, cf. e.g. [29,39].
The goal of this article is to modify the quasi-incompressible models devised as only approximate models for the ideally incompressible one towards physically sensible models. We also focus on identifying dispersive character of the specific models, which seems rather ignored in literature. The dispersion of velocities of P-waves is mainly manifested in solids but it is reported in fluids too, e.g. in the Earth's core [24,35].
In particular, the devised models to which we will focus will exhibit some of (or all) the following attributes: (a) propagation and dispersion of P-waves are controlled in a certain way, (b) the energy balance is preserved at least formally, but in some models even rigorously, (c) the pressure is well defined in a reasonable sense also on the boundary, (d) the equations are consistently written in Eulerian coordinates (i.e. the model is fully convective), (e) in some models, uniqueness of weak solutions holds even in the physically relevant 3-dimensional cases.
We begin with presentation of a rather simple model satisfying (a)-(b) in Sect. 2. In Sect. 3 and 4, we will improve it to satisfy also (c) by including some gradient terms with a conservative or a dissipative character, respectively. Eventually, in Sect. 5, we arrive to fully convective models which satisfy (a)-(d) and, in the case of a multipolar variant, also (e). We will call such new fluidic models satisfying (a)-(d) as semi-compressible fluids.

2.
Simple quasi-incompressible models. The ultimate departure point is the fully compressible Navier-Stokes system which can be rigorously derived from Boltzmann kinetic theory. In term of the velocity field v and the density ρ, this system is where with K v > 0 and G v > 0 a bulk and a shear viscosity coefficients and with the small-strain rate e = e(v) := 1 2 ∇v + 1 2 ∇v and with sph e = (tr e)I/d and dev e = e−(tr e)I/d are the spherical (volumetric) and the deviatoric parts of e, respectively; here " tr " denotes the trace of a matrix. We can also write D ijkl = K v δ ij δ kl + G v (δ ik δ jl +δ il δ jk −2δ ij δ kl /d) with " δ " denoting the Kronecker symbol. The pressure π is assumed to be a function of the density ρ.
The quasi-incompressible model is usually derived by considering a smallperturbation ansatz ( > 0 being a fixed constant -i.e. a homogeneous fluid): where π is the pressure and K the elastic bulk modulus, both with the physical dimension Pa=J/m 3 . Let us remark that 1/K is called compressibility. The mass conservation (1a) written in the form . ρ + div(ρv) = 0 then results to . π/K + div((1+ π K v) = 0 and (1c) results to .
This approximation of the incompressibility constraint, used already e.g. in [22,23], however twists the energy balance, which is desired to be improved by a suitable modification of the approximate system. In term of the velocity field v and the pressure π with > 0 constant, the governing system in its simple variant is considered as .
Here I ∈ R d×d denotes the unit matrix. Beside, g is a given bulk force, e.g. the gravitational force, and K is the mentioned bulk modulus. The force − 2 (div v) v in (4a) is related with the convective term (v·∇)v in the only quasi-(but not fully) incompressible model and was invented (for = 1) by R. Temam [52] (cf. also [53, Ch. III, Sect. 8]) and used also e.g. in [11,21,25,26]. A physical meaning in the energy balance of this term is the variation of the kinetic energy 2 |v| 2 within the volume changes div v and is needed to gain the correct energetics, cf. (7) below. Simultaneously, this forces vanishes in the incompressible limit. Recently, some justification by pullback/pushforward geometrical arguments has been provided in [54], based on [48]. We will consider the system (4) on a bounded domain Ω ⊂ R d evolving in time up to T > 0 a finite fixed time horizon. We abbreviate I := [0, T ].
The system (4) is to be completed by the boundary conditions, say a combination of the Dirichlet and the Navier conditions: where v n = v · n is the normal velocity (as a scalar) with n = n(x) denoting the unit outward normal to Γ, while (·) t denotes the tangential component of a vector, i.e. e.g. v t = v − v n n is the tangential velocity (a vector). Furthermore, σ in (5) is the stress De(v) − πI occurring (4a), and with κ the "friction" coefficient on the boundary acting in the Navier boundary condition (5). Eventually, considering the initial-value problem, we also prescribe the initial conditions The energetics of the model (4)-(5) is revealed by testing (4a) by v. We use (4b) together with the calculus for the convective term which, when used for v = v, gives rise to the bulk force 1 2 (div v)v. Let us note that, in (7), we exploited that is constant. Moreover, we use the calculus of the pressure gradient tested by v: The mentioned test thus gives (at least formally) the energy balance Ω 2 |v(t)| 2 kinetic energy dissipation by a boundary "friction" We will use the standard notation concerning the Lebesgue and the Sobolev spaces, namely L p (Ω; R n ) for Lebesgue measurable functions Ω → R n whose Euclidean norm is integrable with p-power, and W k,p (Ω; R n ) for functions from L p (Ω; R n ) whose all derivative up to the order k have their Euclidean norm integrable with p-power. We also write briefly H k = W k,2 . Moreover, for a Banach space X and for I = [0, T ], we will use the notation L p (I; X) for the Bochner space of Bochner measurable functions I → X whose norm is in L p (I), and H 1 (I; X) for functions I → X whose distributional derivative is in L 2 (I; X). Furthermore, C w (I; X) will denote the Banach space of weakly continuous functions I → X.
The qualification of g uses L 6/5 (Ω; R d ) ⊂ H 1 (Ω; R d ) * and the a-priori estimates (10) are obtained from (9) The former estimate in (11) needs d ≤ 3 and uses the embedding H 1 (Ω) ⊂ L 6 (Ω) and the interpolation L ∞ (I; L 2 (Ω; R d ))∩L 2 (I; L 6 (Ω; R d )) ⊂ L 10/3 (I×Ω; R d ) so that both (v·∇)v and 1 2 (div v)v belong to L 5/4 (I×Ω; R d ). Let us remark that the estimate (10) still does not yield enough integrability of the pressure to allow for uniqueness in the physically relevant 3-dimensional case, while only 2-dimensional case works by estimation like in Proposition 1(iii) below. Therefore, it is not much improvement of the incompressible model where the uniqueness problem has been recognized extremely difficult by involving into seven millennium problems by the Clay institute.
Remark 1 (Energy conservation alternatively). Using the calculus (v·∇)v = div(v ⊗ v) − (div v)v, the force balance (4a) can also be written as as pointed out in [11,Rem.1.2]. Formally, the energetics could be balanced also by other ways, e.g. by the additional pressure instead of the force One could also combine (12) and (13) with suitable weights 2 and −1 to obtain Although the desired energy is conserved in (13) and (14) and the latter case even uses the enhanced pressure 2 |v| 2 + π like in (69a) below, both these variants are physically not justified and seems not relevant because (13) is, in fact, equivalent by replacing the force − 2 (div v) v in (4a) by (∇v) v which does not vanish for incompressible limit.
3. Towards dispersion through conservative gradient terms. The P-waves in real media usually does not propagate with a constant speed, but the speed may depend on their frequency (or equivalently on their wave length). This is referred to as a dispersion and some dispersion is mostly a real effect thus desirable to reflect it in the model in some ways. If the speed increases (or decreases) with the wave length, we speak about normal (or anomalous) dispersion, respectively. Already due to the attenuation by viscous Kelvin-Voigt-like rheology in Sect. 2 controlled by K v causes a normal dispersion, cf. [37,Rem. 6.4.13], but usual modelling ansatz is that K v > 0 is rather very small to allow propagation of waves which have also high frequencies without any essential attenuation. For some fluidic applications where P-waves should propagate for long distances (say thousands of kilometers as in Earth's outer core or oceans), it is thus desirable to have another mechanisms to control dispersion at disposal, not directly causing atenuation. A certain way to realize it is a Sobolev-type regularization of the incompressibility constraint where H is an elastic bulk "hyper-modulus" having the physical dimension Pa/m 2 =J/m 5 . The parabolic equation of the type " . π−∆ . π = f " is sometimes referred to as a Sobolev equation, which is why (15) is called a Sobolev-type regularization here. The term ∆ . π in (15) has been considered in [49,Sect.4.3] but without the K-term for merely numerical purposes.
Sketch of the proof. We can use some approximation, say by the Faedo-Galerkin method, using a nested sequences of finite-dimensional subspaces of H 1 (Ω; R d ) and of H 1 (Ω) whose union is dense in these Sobolev spaces to be used for the approximation of v and of π, respectively. Let us index them by k ∈ N and denote by (v k , π k ) the approximate solutions created by this way. In particular, the discrete variant of (15) leads to the integral identity for all π from the mentioned finite-dimensional subspace of H 1 (Ω) so that π = π k will be a legitimate test to be used for the calculus (17) which then holds for (v k , π k ) too.
The local-in-time existence of (v k , π k ) is due to the standard arguments from the theory of ordinary differential systems together with the subsequent prolongation on the whole time interval I due to the L ∞ (I)-estimates which follows from (18) which holds also for (v k , π k ). These estimates are obtained from (18) by Young's and Gronwall's inequalities and by using the embedding H 1 (Ω; R d ) ⊂ L 6 (Ω; R d ) which holds for d ≤ 3. More specifically, we can see the a-priori estimates By the Banach selection principle, we can select a subsequence converging weakly* to a limit (u, π). This limit is a weak solution to the mentioned initial-boundaryvalue problem for (4a) and (15). The limit passage in the Galerkin approximation for k → ∞ is simply due to the mentioned weak* convergence. For the only two nonlinear terms in (4a) one can use the strong convergence of v k in L 10/3− (Ω; R d ) with > 0 arbitrarily small. This follows (through interpolation) by the Aubin-Lions compact-embedding theorem, using also a uniform bound of a Hahn-Banach extension of . v k on L 1 (I; B) with an arbitrary Banach space B, here on the space from the former inclusion in (11) . π k bounded in L 2 (I; H 1 (Ω)). This is inherited in the limit . π ∈ L 2 (I; H 1 (Ω)). For a smooth Ω, we can use H 2 -regularity of the ∆-operator so that, from ∆ . π = H( . π/K + div v) ∈ L 2 (I; L 2 (Ω)), we obtain also . π ∈ L 2 (I; H 2 (Ω)). As for the uniqueness in the 2-dimensional case, we just expand the estimation which, in the incompressible case, are based on the Gagliardo-Nirenberg inequality Abbreviating v 12 := v 1 −v 2 and π 12 := π 1 −π 2 for two weak solutions (v 1 , π 1 ) and (v 2 , π 2 ), one has the estimate where the constant C is from the Gagliardo-Nirenberg inequality. Then taking > 0 sufficiently small and using Korn's and Gronwall's inequalities, we can see that v 12 = 0 and π 12 = 0. For this, it is also important that v 2 ∈ L 4 (I×Ω; R d ), which follows from the interpolation L ∞ (I; L 2 (Ω)) ∩ L 2 (I; H 1 (Ω)) ⊂ L 4 (I; L 4 (Ω)) holding for d = 2 again by the Gagliardo-Nirenberg inequality.
In comparison with Sect. 2, the nondissipative-gradient model brings an additional mechanism contributing to dispersion not related with any attenuation. To see a more specific character of the dispersion-attenuation behind (15), let us neglect the nonlinear convective terms in (4a) and consider a 1-dimensional linearized situation. We thus consider the simple system Equivalently, we can also write it as ..
In particular, we eliminate v and consider only the displacement u and the pressure to have When subtracting these two equations, we can eliminate also π by using 1 K π x − 1 H π xxx = −u xx and obtain a dispersive equation The term ε .. u xx arisen in (24) is sometimes called micro-inertia; for various models involving such micro-inertia and its role in dispersion we refer in particular to [7,Chap. 6] and [2,27,38,40]. Analogously, εD . u xxxx might be called micro-viscosity. A conventional way to calculate dispersion and attenuation in linear media is to use the ansatz u = e i(wt+x/λ) with the angular frequency w = ω +iγ considered complex with ω, γ ∈ R and the real-valued wavelength λ; here i = √ −1 denotes the imaginary unit. Then ..
. u xxxx = iεDwu/λ 4 , and Ku xx = −Ku/λ 2 , so that the one-dimensional dispersive wave equation (24) yields the algebraic condition When substituting w = ω + iγ so that w 2 = ω 2 − γ 2 + 2iωγ, we obtain two algebraic equations for the real and the imaginary part each, called Kramers-Kronig's relations. More specifically, here Rather exceptionally, these algebraic relations can explicitly be solved: From the latter equation we can read that γ = D/(2 λ 2 ) and then the former equation yields . When substituting η from (25), we obtain the (real part of the) angular frequency Then, realizing that the speed of waves is v = ωλ and substituting η from (25), we obtain which gives a normal dispersion for sufficiently long waves, namely having a length λ > λ crit with the critical wave length solving the equation In terms of the so-called angular wavenumber k = 2π/λ, (28) can be written as cf. Fig. 1. Let us recall that the adjective "normal" for dispersion means that waves with longer lengths propagate faster than those with shorter lengths. Waves with the length λ crit or shorter are so fast attenuated that they cannot propagate at all. Let us note that the coefficient γ determining the attenuation is here inversely proportional to the square power of the wave length. This also reveals the dispersion/attenuation for the simple model from Sect. 2 when putting ε = 0 into (28) and (29); in particular λ crit = D/ √ 4 K. The other extreme case is for the inviscid situation D → 0+ where λ crit → 0+ so that also waves with ultra-high frequencies can propagate. Let us note that the coefficient γ determining the attenuation is here inversely proportional to the square power of the wave length, i.e. naturally the attenuation rises for high frequencies.
From (27), we can also read the so-called group velocity of waves, which is the velocity with which the overall envelope shape of the wave's amplitudes propagates through space and which is defined as dω/dk with k = 2π/λ. Here, (27) with λ = 2π/k yields In this case, we obtain the group velocity: There is an alternative way how to incorporate a conservative gradient term into the model: instead of enhancing (4b) to (15), we can enhance the pressure term in (4a). More specifically, the system (4) is enhanced for .
with some > 0 presumably small, having the physical dimension in meters and determining a certain length scale for possible spacial pressure variations. Modifying still the boundary conditions (16) as with σ = De(v) − (π− 2 ∆π)I, the previous calculus (17) now modifies as This gives again the energy balance (18) only with H replaced by 2 /K and with c replaced by c/ 2 . The Proposition 1(i) and (iii) holds for the boundary-value problem (31)-(32) with the initial condition (6), too. Although the model (31)-(32) is energetically and analytically very similar to (4a) with (15), it differs as far as the dispersive character. Indeed, instead of (23) we now have . v − Dv xx + π x − 2 π xxx = 0, and .
Instead of (25), we now have w 2 − iDw/λ 2 − K/λ 2 − 2 K/λ 4 = 0 so that the Kramers-Kronig's relations look as The In terms of the angular wavenumber k = 2π/λ, we now have, If the dissipative effects do not dominate (i.e. D > 0 is small), (37) gives an anomalous dispersion, i.e. higher-frequency waves (i.e. with longer wave length) propagate faster than waves with lower frequencies, cf. Fig. 2   Remark 2 (Pressure-dependent compressibility). Actually, the compressibility and thus also the speed of P-waves may also depend on the pressure, cf. e.g. [32, Table V] or [42] for the case of seawater. It is not difficult to make K dependent on π in a general continuous manner by replacing the term 1 2K |π| 2 by a general function φ(π), which leads to K(π) = π/φ (π). Indeed, the calculus (17) thus modifies by using π . π/K(π) = φ (π) . π = ∂ ∂t φ(π). Remark 3 (A generalization of the boundary conditions). As the pressure is well defined on the boundary, one can think about a penetrable wall with a flux depending of pressure, which could be modeled by the boundary condition v n = v · n = f (π). Yet, it would expand the energy balance (9) by the boundary flux t 0 Γ 2 |v| 2 f (π)dSdt which seems to corrupt the a-priori estimation even if f (π) would be supposed bounded. 4. Towards dispersion through dissipative gradient terms. An interesting modification has been devised by A.P. Oskolkov [46] who considered (rather for analytical purposes) a parabolic regularization of the incompressibility constraint of the type: where the bulk "hyper-modulus" H has now the physical dimension Pa s/m 2 . Ignoring the convective term, this was later developed under the name "mixed quasicompressibility methods" in [49,Ch. 5]. It actually combines the simple timeregularization (4b) with a mere elliptic regularization div v = 1 H ∆π as used e.g. in [11,18,49].
It can be derived from a fully compressible model when introducing some diffusion into the convective mass transport (1a), i.e.
for some hyper-modulus H. Although it does not seem deducible from the Boltzmann equation by usual procedure toward gas-dynamic limit [4], such a convectiondiffusion transport is occasionally considered, devised by H. Brenner [14,15], cf. also [5,31], although sometimes rather as an artificial regularization only, cf. [30,Eq. (3.172)]. Similarly, [19] devised a diffusion of the deviatoric part of stress in the incompressible case as a regularization. This falls into a general concept of a parabolic perturbation of the (here mass-) conservation law. For a thermodynamical justification see also [55], opposing [47]. Then the small-perturbation ansatz (3) would lead, instead of (4b), to (39). The physical dimension of is m 2 /s. Vaguely speaking, dividing by a "characteristic" velocity of the flow and a "characteristic" length of the system, we obtain a dimensionalless Péclet number expressing dominance of either the convective or the diffusive transport phenomena. Diffusion processes are always dissipative, so it is not surprising that this "regularization" of the incompressibility condition is dissipative, too. The calculus (8) now enhances as  (39)). Let d ≤ 3 and the assumptions (19) hold. Then: (i) the system (4a) and (39) with the initial and boundary conditions (6) and (16) has a weak solution v ∈ C w (I; L 2 (Ω; R d )) ∩ L 2 (I; H 1 (Ω; R d )) and π ∈ C w (I; H 1 (Ω)) ∩ H 1 (I; L 2 (Ω)). (ii) If Ω is smooth, then even π ∈ L 2 (I; H 2 (Ω)) so that, in particular π(t) ∈ C(Ω) for a.a. t ∈ I. (iii)If d = 2, then the weak solution is unique.
Sketch of the proof. Again we use the Faedo-Galerkin method as in the proof of Proposition 1. From (42) written for the approximate solutions, we obtain a-priori estimates In fact, (43a) is the former estimate (21a). Further a-priori estimate follows by testing (39) by . π, using that (43a) yields div v k ∈ L 2 (I×Ω). Thus we obtain π k L ∞ (I;H 1 (Ω)) ≤ C, and . π k L 2 (I×Ω) ≤ C .
By comparison, we obtain also an a-priori information for a Hahn-Banach extension of . v k in L 2 (I; H 1 (Ω; R d ) * )+L 5/4 (I×Ω; R d ), and then, by the Aubin-Lions theorem used for such Hahn-Banach extension as in [50,Sect.8.4], we can pass to the limit as in the proof of Proposition 1.
As to (ii), assuming also Ω to be smooth so that H 2 -regularity of the Laplacian is at disposal, from ∆π = H(div v + . π/K) ∈ L 2 (I×Ω) we obtain the L 2 (I; H 2 (Ω; R d ))estimate of π.
As for (iii), we just modify the estimate (22) appropriately: more specifically, . v − Dv xx + π x = 0, and Eliminating v, we can write it as ..
When subtracting them, we can eliminate π by using 1 . u xx and obtain a dispersive equation Using again the ansatz u = e i(wt+x/λ) , we have εDu xxxx = εDu/λ 4 , and the onedimensional dispersive wave equation (46) yields the algebraic condition Considering again w = ω + iγ, the Kramers-Kronig's relations (26) now looks as From the latter equation we can read that γ = (D+ε )/(2 λ 2 ) and then, realizing that the speed of waves v = ω 2 λ 2 , the former equation in (48) yields Like (28), we can see again the normal dispersion with the effect that high-frequency waves cannot propagate at all. Again, γ determines the attenuation which naturally rises for higher frequencies.
Remark 4 (Multipolar fluids). Another enhancement by dissipative gradient terms can exploit the concept of the 2nd-grade nonsimple fluids, devised by E. Fried and M. Gurtin [33,34] and earlier, even more generally and nonlinearly as multipolar fluids, by J. Nečas at al. [6,44,45]. More specifically, the 2nd-grade multipolar variant of the model (4a) with (39) reads as with H a positive-definite tensor of hyper-viscous moduli. This needs enhancement of the boundary condition (16), say v n = 0, σ n − div S (H∇e(v) n) t + κv t = 0, and H∇e(v): .
( n ⊗ n) = 0 (50) where σ = De(v)−πI and where "div S " is the surface divergence defined as div S (·) = tr ∇ S (·) with tr(·) denoting the trace and ∇ S denoting the surface gradient given by ∇ S v = (I − n⊗ n)∇v = ∇v − ∂v ∂ n n. This model gives a new term H∇e(v): . ∇e(v) into the dissipation rate in the energetics (42), which improves the estimate (43a) by replacing H 1 (Ω; R d ) by H 2 (Ω; R d ). The analysis of the dispersion is more involved: (44) now looks as . v − Dv xx + π x + hv xxxx = 0, and Moreover, one can also think about combination of this dissipative multipolar concept with the conservative-gradient models as (31). This would lead to the system The energetics of (52) when completed by the boundary conditions (32) can be revealed by using again the calculus (33). This gives again the energy balance (18) only with H replaced by 2 /K and with c replaced by c/ 2 , and here in addition with the dissipation rate H∇e(v): . ∇e(v). The estimate like (67) below will yield uniqueness. The strain-rate gradient versus pressure gradient are gradienttheoretical concepts known as Aifantis [1][2][3] versus Eringen [28], respectively. In (52a), we combined both concepts.

5.
Truly convective "semi-compressible" models. All the above models may still be referred under the name "quasi" because they are not fully mechanically consistent when mixing convective time derivative in (4a) or in (31a) or (49a), i.e. the Eulerian description, with the partial time derivative in (4b), (15), (39), or (49b), i.e. the Lagrangian description. Now we want to improve (some of) these models to be formulated consistently fully convectively, i.e. in Eulerian coordinates. It is also likely that the criticism [47] of the diffusive term in (39) or (49) because of lack of Galilean invariance is thus suppressed so that the justification of these diffusive terms in [55] is fully in effect.
In fact, the derivation of the quasi-incompressible model (4b) from the fullycompressible equation (1a) under the small perturbation ansatz (3) was not precise. In fact, from (1a) when using (3), we have Instead of (4b), it should lead rather to when assuming |π| K to replace the coefficient 1/(K+π) by 1/K. Mechanically, one can say that replacing the extensive variable (density) ρ in (1) by the intensive variable (pressure) π in (4), should lead to replacement the transport equation (1a) for the extensive variable ρ by a transport equation (54) for the intensive variable π; let us remind that intensive variables are those whose magnitude is independent of the size of the system (volume) whereas the extensive ones are those whose magnitude is additive for subsystems.
However, replacement of (1a) by (54) brings analytical difficulties because ∇π, similarly as in the nonconvective variant (4a), is not a-priori estimated in the model (4a) with (54). Even the by-part integration Ω (v·∇π) π dx = − Ω π(v·∇ π) + π(div v) π dx which could eliminate ∇π from a weak formulation seems problematic because then one would need a strong convergence in ∇v or in π, which does not seem directly at disposal. Similarly, the convective modification of the model (4a) with (15) from Section 3 seems analytically problematic.
Yet, some of the above discussed models bear the physically relevant convective modification which is simultaneously amenable for analysis. In particular, for (39), the force-equilibrium equation (4a) is to be now slightly modified to obtain the model Beside the "hydrostatic" pressure π, we now also see an "internal" pressure due to the elastic energy π 2 /(2K) in (55a). This contribution due to internal energy (here elastic but it can also be e.g. magnetic or chemical) is a thermodynamically justified effect, cf. e.g. [41], which disappears only in ideally incompressible models. Here, it arises when using the Green formula so that the calculus (41) now enhances as The boundary conditions can again be considered as in (16). The energetics (42) remains the same because the additional terms 1 2K π 2 in (57) just cancel with the additional pressure terms coming by the test of (55a) by v and using the enhanced boundary condition (32). (55)). Let d ≤ 3, , K, H, c > 0, κ ≥ 0, D be symmetric positive definite, and (19) hold. Then the system (55) with the boundary conditions (16) possesses a weak solution v ∈ C w (I; L 2 (Ω; R d )) ∩ L 2 (I; H 1 (Ω; R d )) and π ∈ C w (I; H 1 (Ω)).
The convergence in the additional nonlinear term v k ·∇π k in (55b) is easy because ∇π k converge weakly* in L 2 (I×Ω; R d ) due to (43b) while v k converge strongly in L 2 (I×Ω; R d ) due to (43a) and the Aubin-Lions compact-embedding theorem using also (58a). In addition, we need strong convergence of π k in the additional pressure 1 2K π 2 k , which is again simple when using the Aubin-Lions theorem with (58b) with the only technicality that . π k in (58b) is to be understood as a Hahn-Banach extension or the estimate (58b) is to be weakened by using only seminorms, cf. [50,Sect.8.4].
One should note that the rigorous energy conservation behind the model (55) is not granted by our estimates because the inertial force . v does not enjoy enough regularity to be eligible for being tested by v. Also the uniqueness for d = 3 is not obvious, cf. Remark 9. For these reasons, one may still consider the fully convective variant of the multipolar model (49), which looks as together with the boundary conditions (16) combined with (50), i.e. here v n = 0, .
( n ⊗ n) = 0, and This allows for improving the estimate of . π from (58b). Another important attribute of the semi-compressible convective multipolar model (59) is that it allows for uniqueness even in 3-dimensional situations.
Proof. First, testing (55a) in its Galerkin approximation by the Galerkin approximation v k of v and using the calculus (57) for (59b), we obtain the energy balance (62) for the approximate solutions. Of course, (59b) is understood here in its Galerkin approximation for which the test by π k used in (57) is legitimate. From this, one can read the a-priori estimates Second, by testing (59b) by . π k , we obtain Using the embedding H 2 (Ω) ⊂ L ∞ (Ω) so that due to (63a), by the Gronwall inequality we improve the estimate (63b) as π k L ∞ (I;H 1 (Ω)) ∩ H 1 (I;L 2 (Ω)) ≤ C .
Third, let us realize that, by (65), div(π 2 k I/2) = π k ∇π k is bounded in L ∞ (I; L 3 (Ω; R d )). Therefore, we can test (59a) in its Galerkin approximation by The first right-hand side term is to be estimated as Ω g · (Ω;R d ) ; here the assumption g ∈ L 2 (I×Ω; R d ) is employed. When estimating the last two term as with some C depending on and involving also a constant from the Korn inequality on Ω, and when using again (64), by the Gronwall inequality we still obtain the estimate The limit passage of selected weakly* convergent subsequences towards weak solutions (v, π) to (59) is then easy. Then, from (59b) we can also see that ∆π = H K ( . π + v·∇π) + Hdiv v ∈ L 2 (I×Ω) due to (65). As to the uniqueness, enhancing the calculus (22), we have for a.a. time instances t ∈ I (with t omitted in the following formulas for notational simplicity) that where v 12 := v 1 −v 2 and π 12 := π 1 −π 2 for two weak solutions (v 1 , π 1 ) and (v 2 , π 2 ). The particular terms on the right-hand side can be estimated by Hölder's and Young's inequalities. In particular, by the Gagliadro-Nirenberg inequality together with Korn's inequality for the 2nd-grade nonsimple materials [37, Sect. 5.2], we have v 12 for any 1 ≤ r ≤ 6. Here we will use it for r = 12/5. Thus, we can estimate Taking > 0 small enough, the -terms on the right-hand sides of (68) can be absorbed in the left-hand side of (67) and then, using that t → e(v 1 (t)) 2 L 6 (Ω;R d×d ) , t → π 1 (t)+π 2 (t) 2 L 4 (Ω) and t → v 1 (t) 2 L ∞ (Ω;R d ) + ∇π 1 (t) 2 L 2 (Ω;R d ) are L 1 (I), we treat it by Gronwall inequality.
Eventually, (iii) follows by the H 2 -regularity of the ∆-operator on smooth domains the left-hand side of (59b) is valued in L 2 (Ω). Note that . π ∈ L 2 (I; L 2 (Ω)) due to (65) while v·∇π and div v are in a smaller space due to (63). Similarly, we can use the H 4 -regularity for the bi-harmonic operator div 2 (H∇e(v)) and that . v ∈ L 2 (I; L 2 (Ω; R d )) due to (66), and to obtain the L 2 (I; H 4 (Ω; R d )) regularity of v from (59a).
Remark 6 (Bernoulli principle). Using the calculus (v·∇)v = (∇ × v) × v + 1 2 ∇|v| 2 , the model (55) can be written in the form This recovers the celebrated Daniel Bernoulli's principle [8] which states that an increase in the speed of a fluid occurs simultaneously with a decrease in pressure and vice versa. The "Bernoulli" pressure in (69) is thus the sum of the so-called dynamic pressure 2 |v| 2 , the pressure related with the internal energy of the fluid 1 2K π 2 , and the "hydrostatic" pressure π. The model (59) bears an analogous form and is thus consistent with the Bernoulli principle, too.
cf. e.g. [37,Formula (7.7.27)]. This gives the same energetics as the model (52). However, the limit passage in the Galerkin approximation would now need also the strong convergence of ∇π due to the Korteweg stress in (70a), which seems difficult. In principle, a higher-gradient "multipolar-type" term like −div 2 ∇ 2 π instead of 0 in (70b) with correspondingly enhanced boundary conditions would help. A certain variant of this model even with as a variable subjected to the continuity equation, combined also with the diffusion in (59b) and amenable for existence of solutions, has been developed in [16]. For a different, incompressible variant with fixed combining incompressible model with a convective-diffusive transport of a scalar pressure-like variable, amenable for existence of solutions, see [17].
Remark 8 (Pressure-dependent compressibility III). As in Remark 2, K = K(π) can be considered, too. Now, (53) suggests to take it so that K (0) = 1, cf. also Remark 5. We thus imitate the factor 1/(K+π) occurring in (53) at least for small pressures |π| while allowing still for keeping positivity of K(π) for bigger negative pressures. Actually, dependence of K on π allows for modeling dependence of the sound speed on the pressure, which is a well observed effect. More specifically, (53) suggests the expansion for small pressure as v p (π) = K+π = K 1 + π K ∼ v p,0 1 + π 2K with v p,0 = K .
Yet, actual fluids may exhibit even more pronounced dependence on pressure. E.g. for water with K = 2.15 GPa and v p,0 = 1500 m/s this formula gives v p (π) ∼ 1500 + 0.348 × 10 −6 π while the measurements as in [20,Tab. IV] shows v = 1521.46 m/s for π = 0 and v = 1538.12 m/s for π =100 bar =10 MPa at temperature 20C, which gives v(π) = v 0 +2.146×10 −6 π. Even more pronounced dependence is in [43,Tab.5], reporting the compressibility at 0 • C in 10 6 m 3 /J=1/MPa is 5.0792×10 −4 −1.3389× 10 −6 π for π in MPa. In general, sound speed provides high accuracy experimental data exploited to for state equation from which the pressure dependence of the bulk modulus (or compressibility) is determined, in the particular case of water cf. [56] for the Intl. Assoc. for the Properties of Water and Steam 1995 standard.
Remark 9 (Uniqueness in the model (55)). It is known that the incompressible Navier-Stokes equation enjoys uniqueness and full regularity if the pressure is in the class L 1 (I; L 3/2 (Ω)), cf. [9] or also e.g. [10]. Here, in the model (55), the overall pressure 1 2K π 2 + π is controlled even in L 1 (I; L 3 (Ω)), and the additional force 2 (div v) v has a bilinear form similar to the convective term (v·∇)v. There seems a chance to augment the uniqueness proof for the whole system (55) involving also the other bilinear term 1 K v·∇π in (55b). Remark 10 (Normal dispersion by convective conservative gradient terms). To model dispersion in waves that can effectively propagate on long distances, one should use rather conservative gradient terms than dissipative ones in order not to make a substantial attenuation of the waves. For the normal dispersion, this can be realized through (15), cf. (28). But the convective variant of (15) seems problematic for analysis. Thus designing a convective model which would realize normal dispersion through a conservative gradient term seems to be a challenging open problem.