Viscoelastic rate-type fluids with stress-diffusion

# PDE analysis of a class of thermodynamically compatible viscoelastic rate-type fluids with stress-diffusion

Miroslav Bulíček Mathematical Institute, Charles University, Faculty of Mathematics and Physics, Sokolovská 83, 186 75 Prague 8, Czech Republic Josef Málek Mathematical Institute, Charles University, Faculty of Mathematics and Physics, Sokolovská 83, 186 75 Prague 8, Czech Republic Vít Pruša Mathematical Institute, Charles University, Faculty of Mathematics and Physics, Sokolovská 83, 186 75 Prague 8, Czech Republic  and  Endre Süli Mathematical Institute, University of Oxford, Woodstock Road, Oxford OX2 6GG, United Kingdom
###### Abstract.

We establish the long-time existence of large-data weak solutions to a system of nonlinear partial differential equations. The system of interest governs the motion of non-Newtonian fluids described by a simplified viscoelastic rate-type model with a stress-diffusion term. The simplified model shares many qualitative features with more complex viscoelastic rate-type models that are frequently used in the modeling of fluids with complicated microstructure. As such, the simplified model provides important preliminary insight into the mathematical properties of these more complex and practically relevant models of non-Newtonian fluids. The simplified model that is analyzed from the mathematical perspective is shown to be thermodynamically consistent, and we extensively comment on the interplay between the thermodynamical background of the model and the mathematical analysis of the corresponding initial-boundary-value problem.

###### 2000 Mathematics Subject Classification:
Primary 35Q35, 76A05, 76A10
The authors acknowledge the support of the ERC-CZ project LL1202, financed by MŠMT

## 1. Introduction

The main goal of this study is to establish the long-time existence of large-data weak solutions to the following simplified set of governing equations encountered in the mechanics of incompressible non-Newtonian fluids. This simple model shares many standard properties with more complex viscoelastic rate-type fluid models with stress-diffusion, used in applications.

For any given Lipschitz domain , , , and for any , we set and we seek the functions satisfying, in , the system of PDEs:

 (1.1) div→v =0, (1.2) ϱ(∂t→v+div(→v⊗→v))−divT =0, (1.3) T =−pI+2νD−σ(∇b⊗∇b), (1.4) ν1∂tb+ν1div(b→v)+μ(b2−b)−2σb2Δb =0,

together with the following boundary conditions, where denotes the unit outward normal vector on :

 (1.5) →v=→0 and ∇b⋅→n=0 on% ∂Ω×(0,T),

and the initial conditions:

 (1.6) →v(0,⋅)=→v0(⋅) and b(0,⋅)=b0(⋅) in Ω.

Here, is the velocity, stands for the dyadic product (), is the Cauchy stress, is the spherical stress (the modified pressure), is a scalar quantity characterizing the volumetric elastic changes exhibited by the fluid; denotes the symmetric part of the velocity gradient, i.e., ; is the viscosity, is the density, and , and are other material constants, all of which are positive.

Provided further that (the definitions of all function spaces appearing in the paper are given in Section 4)

 (1.7) →v0∈L20,div(Ω)d,b0∈W1,2(Ω),b0>0  % a.e. in Ω,b0∈L∞(Ω) and 1b0∈L∞(Ω),

we will prove, see Section 5, the following statement:

For any Lipschitz domain , , , , , , and for any satisfying (1.7), there exists a couple solving the problem (1.1)–(1.6) in a weak sense.

The precise definition of weak solution is given in Section 4.

Note that if in the above equations the problem splits into two separate problems: a transport equation with damping for and the standard Navier–Stokes equations for . Since the work of Leray [30] on the incompressible Navier–Stokes equations, see also [22, 8], the question of long-time existence of large-data weak solutions has also been explored for more general classes of viscous fluids; further references will be given in the next section. It is then natural to attempt to advance this program by exploring how large the class of fluids might be for which long-time and large-data existence of weak solutions can be established. This task is particularly interesting if one considers fluid models that include terms that are associated with elastic properties of the material. This paper aims to contribute to the study of this question.

Note also that if then not only is there a diffusion term present in the equation for , but also the Korteweg stress appears in the expression for the Cauchy stress featuring in the equation for the balance of linear momentum. This structure of the governing equations results from a careful derivation of the model based on the thermodynamical approach established for viscoelastic rate-type fluids in [42] and further refined in [39] (see also [36]), and extended to rate-type fluids with stress-diffusion in [37]. This thermodynamical approach automatically guarantees that the resulting model is consistent with the laws of thermodynamics. From the point of view of PDE analysis, the approach readily provides the relevant a priori estimates upon which the analysis is based, and which would be otherwise completely nontrivial to discover solely from the PDE system (1.1)–(1.4).

Thus, the second goal of this paper is to highlight how the analysis of PDEs for complex systems such as (1.1)–(1.4) is related to the thermodynamical approach used in the derivation of the model and to show that the problem (1.1)–(1.6) therefore merits investigation.

The paper is structured as follows. In the next section, we place the model under consideration into a hierarchy of phenomenological models. In Section 3, we sketch its derivation. As was noted above, this thermodynamical approach automatically guarantees the consistency with the second law of thermodynamics and provides directly the a priori estimates available for the system. In Section 4, we re-derive these a priori estimates, but now from a purely PDE analytical point of view, and then, after introducing the appropriate function spaces, we precisely state the main result of the paper. The proof of the theorem is contained in Section 5. The existence of solutions to the finite-dimensional approximating system, upon which the proof of the main result is based, is given in the Appendix.

## 2. Incompressible non-Newtonian fluids: a brief overview

For incompressible fluids, the Cauchy stress tensor decomposes as

 (2.1) T=−ϕI+S,

where is a scalar unknown quantity and is related to the symmetric part of the velocity gradient via the so-called constitutive relation. Three classes of constitutive relations are frequently used; they are of the following forms:

 (2.2) G(S,D) =O, (2.3) =O, (2.4) =O.

Here stands for an arbitrary continuous tensorial function and signifies an objective time derivative.

We will give examples and discuss briefly the usefulness of these classes. Then we shall look at the associated initial-boundary-value problems for internal flows from the perspective of the long-time existence of large-data weak solutions.

Besides the (linear, i.e. Newtonian) Navier–Stokes fluid for which the constitutive relation is given by with , the class (2.2) includes for example the Bingham fluid, described by , as well as various generalizations of activated or non-activated (stress) power-law fluids that can be characterized by

 2ν(|S|2,|D|2)(|D|−d∗)+D=2α(|S|2,|D|2)(|S|−τ∗)+S,

which involves two non-negative activation parameters and , but only one is assumed to be positive in a given model. The possibility of (2.2) to incorporate nonlinear relations between and and to include sudden changes from one type of response to another one (within the considered class) makes the class (2.2) suitable for modeling even mixing type phenomena, and the class (2.2) is therefore very popular in many areas (see [38, 5] and [36] for details and further references).

On the other hand, fluids described by the class (2.2) are not capable of capturing fundamental phenomena such as stress relaxation or nonlinear creep observed in most fluids with complex microstructure. If these phenomena are of interest, then the class (2.3) is a suitable choice.

As was already mentioned above in connection with (2.3), denotes an objective time derivative. The set of objective derivatives includes the upper-convected Oldroyd derivative defined by

the Jaumann–Zaremba (corotational) derivative specified by

or the Gordon–Schowalter derivative defined by

Assuming that , we see that the appropriate choice of material parameter and the objective derivative lead to the popular Maxwell, Oldroyd-B and Johnson–Segalman models respectively:

 τ▽A+A=2ν1D with ν=0 and τ=ν1E,τ▽A+A=2ν1D with ν>0 and τ=ν1E,τ□A+S=2aD with ν>0 and a∈[−1,1].

Although these models are widely used, there are several subtle issues regarding their physical underpinnings. These include the ambiguity with respect to the choice of the objective derivative, the possibility of the derivation of the model at a purely macroscopic level, the consistency of the models with the second law of thermodynamics, the extension of the models to the compressible setting, and the inclusion of thermal effects into the models.

To address some of these issues Rajagopal and Srinivasa [42] provided a simple, yet general method for the derivation of thermodynamically consistent rate-type fluid models. The method is based on the concept of natural configuration and on the knowledge of constitutive relations for two scalar quantities: the Helmholtz free energy (characterizing how the material stores energy) and the rate of entropy production (characterizing how the material dissipates energy). Using the method, it is possible to derive new thermodynamically compatible classes of non-linear viscoelastic rate-type fluid models, and to investigate the conditions under which these models reduce to the standard Oldroyd-B, Maxwell and Burgers models. More recently, the approach has been refined and extended to compressible fluids and thermal processes (see [39, 36, 23]).

Rate-type fluid models with stress-diffusion, c.f. (2.4), are popular in the modelling of shear and vorticity banding phenomena; see, for example, the reviews [12, 17, 14]. The ad hoc addition of the regularizing term to the rate-type equation (2.3) raises doubts about the consistency of the resulting model with the second law of thermodynamics and calls for the specification of appropriate boundary conditions for . These issues are addressed in the recent study [37]. Here, in Section 3, we sketch the derivation under several simplifying assumptions.

As to the question of existence of solutions, for any time interval and any reasonable set of data, a satisfactory theory, within the context of weak solutions, has been established during the last fifty years for a general class of fluids given by (2.2) with a continuous monotone function (see [30, 8, 25, 27, 28, 31, 4, 2, 34, 35, 19, 7, 43, 13, 3, 5, 6, 40]).

Regarding the long-time and large-data existence theory for rate-type fluids in three dimensions much less is known (see [32, 24, 41] or [29]).

One expects that the addition of the term to a rate-type equation will automatically improve its mathematical properties and a well-founded theory for a class of fluids of the type (2.4) will therefore emerge. However, most of the results that are in place concern two-dimensional (or steady) flows, see [15, 1, 10, 9, 33, 16] or require stronger regularization, see [26].

Unlike the above two-dimensional results, the general analytical approach developed and presented in Sections 4 and 5 concerns -dimensional flows. Although the tensor is ultimately “replaced” here by a scalar quantity , the analysis is performed in such a way that it avoids properties which are unlikely to be true for a general tensor .

## 3. Derivation of the model

We denote by the superscript the material derivative, i.e., for a scalar function ; if stands for a longer expression, we write instead of in order to avoid confusion. For vector- or tensor-valued quantities, the same notation is used and the above formula is applied to each component.

The fundamental balance equations of mass, linear momentum and energy as well as the formulation of the second law of thermodynamics are the following:

 (3.1) ˙ϱ =−ϱdiv→v, (3.2) ϱ˙→v =divT,T=TT, (3.3) ϱ˙e =T:D−div→je, (3.4) ϱ˙η =ϱζ−div→jη with % ζ≥0,

where is the density, is the velocity, is the specific internal energy, is the specific entropy, is the symmetric part of the velocity gradient, is the Cauchy stress tensor, and are the energy and entropy fluxes, and stands for the specific rate of entropy production.

Introducing the Helmholtz free energy by

 ψ:=e−θη,

where stands for the (positive) temperature, the equations (3.3) and (3.4) lead to

 T:D−ϱ˙ψ−div(→je−θ→jη)−ϱη˙θ−∇θ⋅→jη=ϱθζ with ζ≥0.

In what follows, we restrict ourselves to isothermal processes (referring the reader to [37] for the development of models where both mechanical and thermal processes are included). Consequently, is constant and the last identity reduces to

 (3.5) T:D−ϱ˙ψ−div(→je−θ→jη)=ξ with ξ≥0,

where denotes the rate of dissipation. Note that (3.5) simplifies further to provided that , which we do not require here, however.

The approach that we will exploit is based on the concept of natural configuration that splits the total response described by the deformation tensor between the current and initial configuration into the purely elastic (reversible, non-dissipative) part that operates between the natural and current configuration and the dissipative part that maps from the reference to the natural configuration, i.e. . In analogy with the relations and , we set and define . We also set .

In order to derive the constitutive relations for the Cauchy stress tensor and the energy and entropy fluxes we start with postulating the constitutive relation for the Helmholtz free energy in the form:

 (3.6) ψ=ψ0(ϱ)+μ2ϱ(trBkp(t)−3−lndetBkp(t))+σ2ϱ|∇trBkp(t)|2=:ψ0(ϱ)+πϱ,

with , constant, and

 π:=μ2(trBkp(t)−3−lndetBkp(t))+σ2|∇trBkp(t)|2.

The following identities (see [39] for their derivation) hold:

 (3.7) ˙¯¯¯¯¯¯¯¯¯¯¯¯Bkp(t) =LBkp(t)+Bkp(t)LT−2Fkp(t)Dkp(t)FTkp(t), (3.8) ˙¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯trBkp(t) =2Bkp(t):D−2Ckp(t):Dkp(t), (3.9) ˙¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯lndetBkp(t) =2I:D−2I:Dkp(t), (3.10) 12˙¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯|∇trBkp(t)|2 =∇trBkp(t)⋅(∇˙¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯trBkp(t)−(∇→v)∇trBkp(t)).

Setting

 b:=trBkp(t)

and using the formula

 ∇b⋅∇˙b=div(˙b∇b)−Δb˙b,

the above identities (3.6), (3.8) and (3.10) take the form

 (3.11) ψ =ψ0(ϱ)+πϱ with π=μ2(b−3−lndetBkp(t))+σ2|∇b|2, (3.12) ˙b =2Bkp(t):D−2Ckp(t):Dkp(t), (3.13) ˙¯¯¯¯¯¯¯¯¯¯¯¯|∇b|22 =−((∇b⊗∇b)+2ΔbBkp(t)):D+div(˙b∇b)+2ΔbCkp(t):Dkp(t).

Consequently, using also (3.1) and (3.9), we compute

 ϱ˙ψ =ϱψ′0(ϱ)˙ϱ+˙π−π˙ϱϱ =(μ(Bkp(t)−I)−σ((∇b⊗∇b)+2ΔbBkp(t))):D −(μ(Ckp(t)−I)−2σΔbCkp(t)):Dkp(t)+div(σ˙b∇b)−(ϱ2ψ′0(ϱ)−π)div→v.

Inserting the result into (3.5) we obtain

 ξ =(T−μ(Bkp(t)−I)+σ((∇b⊗∇b)+2ΔbBkp(t))+(ϱ2ψ′0(ϱ)−π)I):D +(μ(Ckp(t)−I)−2σΔbCkp(t)):Dkp(t)−div(→je−θ→jη+σ˙b∇b) with ξ≥0.

Finally, setting111Note that the first line states how the entropy flux is related to the energy flux , while in the second and third lines we have introduced new notation to simplify the formulae.

 →jη =→je+σ˙b∇bθ, pNSth :=ϱ2ψ′0(ϱ), Tel :=−(pNSth−π)I+μ(Bkp(t)−I)−σ((∇b⊗∇b)+2ΔbBkp(t)),

we arrive at

 (3.14) ξ =(T−Tel):D+(μ(Ckp(t)−I)−2σΔbCkp(t)):Dkp(t) with ξ≥0,

which provides key information in the derivation of an entire hierarchy of models, ranging from (compressible/incompressible) Euler and Navier–Stokes fluids, through Maxwell and Oldroyd-B fluids, to diffusive Maxwell and diffusive Oldroyd-B fluids. In addition, the procedure also clarifies how these fluids store and dissipate energy and consequently what quantities are a priori controlled by the initial data (which then indicates what a priori information/estimates one can use in the analysis of the model).

To illustrate the procedure and to derive a simplified (toy) model that is then analyzed in the second part of this paper, we rewrite (3.14) in the following way. Using the subscript to signify the deviatoric (traceless) part of a tensor, i.e. , noticing that

 A:B=Aδ:Bδ+trAtrB3

and denoting the mean normal stress by , i.e.,

 m:=13trT,

we rewrite (3.14) as

 (3.15) ξ=(Tδ−[Tel]δ):Dδ+(m−trTel3)div→v+(μ[Ckp(t)]δ−2σΔb[Ckp(t)]δ):[Dkp(t)]δ+(μ(b−3)−2σbΔb)trDkp(t)3 with ξ≥0,

where

 [Tel]δ =μ[Bkp(t)]δ−σ((∇b⊗∇b)δ+2Δb[Bkp(t)]δ), trTel3 =−pNSth+13[μ(b−3)−σ|∇b|2−2σbΔb]+π with π=μ2(b−3−lndetBkp(t))+σ2|∇b|2.

Referring to [36] for details, we wish to emphasize that the structure of (3.15) is rich enough to incorporate the compressible and incompressible Euler and Navier–Stokes, Maxwell, Oldroyd-B, Giesekus, diffusive Maxwell, diffusive Oldroyd-B, and diffusive Giesekus models.

Here, we make three simplifications: the fluid is assumed to be incompressible, the density is taken to be uniform and the elastic part of the deformation is supposed to be purely spherical222This is reminiscent of the response leading to the Euler fluid.. This means that

 (3.16) div→v=0,ϱ is a positive constant and [Ckp(t)]δ=O.

Then,

 Ckp(t)=[Ckp(t)]δ+trCkp(t)3I=O+trBkp(t)3I=b3I,

where is the zero tensor; note that also vanishes. As a consequence of these simplifications, (3.12) and (3.15) reduce to

 (3.17) ˙b=−23btrDkp(t),

and

 (3.18)

where

 [Tel]δ=−σ(∇b⊗∇b)δ.

Requiring that

 (3.19) Tδ−[Tel]δ=2νDwith ν>0,μ(b−3)−2σbΔb=2ν1trDkp(t)3with ν1>0,

we obtain

 (3.20) ξ=2ν|D|2+2ν1|trDkp(t)|29.

Referring then to (3.17) we deduce that (3.19) leads to

 (3.21) T=mI+2νD−σ(∇b⊗∇b)δ=ϕI+2νD−σ(∇b⊗∇b),ν1˙bb+μ(b−3)−2σbΔb=0,

and (3.20) implies that

 (3.22) ξ=2ν|D|2+ν12∣∣∣˙bb∣∣∣2.

Note that and consequently cannot be determined constitutively because of the incompressibility constraint. Consequently, it is an additional unknown quantity. In what follows we relabel by , as this symbol is used in most studies concerning incompressible fluids. To conclude, recalling also the definition of the material time derivative, the system of governing equations takes the form

 (3.23) div→v =0, (3.24) ϱ(∂t→v+div(→v⊗→v))−divT =0, (3.25) T =−pI+2νD−σ(∇b⊗∇b), (3.26) ν1∂tb+ν1div(b→v)+μ(b2−3b)−2σb2Δb =0.

We wish to emphasize that despite the simplification (see (3.16)), the reduced model (3.23)–(3.26) retains all of the material parameters of the original problem in the framework. It is also worth emphasizing that not only does the diffusion term appear in the equation for the mean normal elastic stress , but also the Korteweg stress is present in the constitutive relation for the Cauchy stress.

Setting

 ~b:=b3 and ~σ:=9σ

we can rewrite (3.23)–(3.26) as follows:

 (3.27) div→v =0, (3.28) ϱ(∂t→v+div(→v⊗→v))−divT =0, (3.29) T =−pI+2νD−~σ(∇~b⊗∇~b), (3.30) ν1∂t~b+ν1div(~b→v)+3μ(~b2−~b)−2~σ~b2Δ~b =0,

which, upon relabeling , and coincides with the system (1.1)–(1.4) that is analyzed in the next sections.

We conclude this section by showing that the choice of the constitutive relations for and directly yields the relevant a priori estimates in a straightforward way. Note first that with the simplifications (3.16) and relabeling ( by and by ) the constitutive relation (3.11) reduces to

 (3.31) ϱψ=ϱψ0(ϱ)+μ2(b−1−lnb)+σ2|∇b|2.

Next, by taking the scalar product of (3.28) with we obtain

 (3.32) ∂t(ϱ|→v|22)+div(ϱ|→v|22→v)−div(T→v)+T:D=0.

Using the equation (3.5) then leads to

 (3.33) ∂t(ϱ|→v|22)+ϱ˙ψ+ξ+div(ϱ|→v|22→v−T→v+→je−θ→jη)=0.

Inserting (3.20) and (3.31) in (3.33), and recalling that and is constant, we arrive at

 (3.34) 12∂t(ϱ|→v|2+σ|∇b|2+μ(b−lnb))+2ν|D|2+ν12∣∣∣˙bb∣∣∣2+div(12(ϱ|→v|2+μ(b−lnb)+σ|∇b|2)→v−T→v+σ˙b∇b)=0.

The next section starts with the derivation of (3.34) from a PDE-analytic point of view.

## 4. Function spaces and statement of the main result

In this section, we provide the precise statement of the main theorem regarding the long-time existence of large-data weak solutions to the problem (1.1)–(1.7).

To motivate the choice of function spaces, we first briefly derive the basic energy estimate. For these formal calculations, we assume that all quantities are well defined; in particular, we assume that is strictly positive. Taking the scalar product of (1.2) with and using the incompressibility constraint (1.1), we obtain the following identity (for the evolution of the kinetic energy):

 12∂t(ϱ|→v|2)+12div(ϱ→v|→v|2)−div(T→v)+T:D=0.

Employing further (1.3) and again (1.1) we obtain

 (4.1) 12∂t(ϱ|→v|2)+2ν|D|2=σ(∇b⊗∇b):D−12div(ϱ→v|→v|2)+div(T→v).

Next, dividing (1.4) by and multiplying the result by , we obtain

 (4.2) ν1b2|∂tb +div(b→v)|2=(−μ(1−b−1)+2σΔb)(∂tb+div(b→v)) =−μ(∂t(b−lnb)+div(→v(b−lnb)))−σ∂t|∇b|2+2σdiv(∂tb∇b) +2σdiv(∇b(div(b→v)))−2σ∇b⋅∇(div(b→v)).

Using the fact that several times, the last term on the right-hand side can be rewritten as

 ∇b⋅∇(div(b→v)) =∇b⋅∇(→v⋅∇b) =(∇b⊗∇b):D+12div(→v|∇b|2).

Using this expression, (4.2) takes the form

 (4.3) ∂t(σ|∇b|2+μ(b−lnb))+ν1b2|∂tb+div(b→v)|2=−2σ(∇b⊗∇b):D +div(2σ(∂tb+div(b→v))∇b−(μ(b−lnb)+σ|∇b|2)→v).

Finally, dividing (4.3) by 2 and adding the result to (4.1), we deduce the energy identity

 (4.4) 12∂t(ϱ|→v|2+σ|∇b|2+μ(b−lnb))+2ν|D|2+ν12b2|∂tb+div(b→v)|2=div(T→v+σ(∂tb+div(b→v))∇b−12(ϱ|→v|2+μ(b−lnb)+σ|∇b|2)→v).

Hence, integration of the result over , using integration by parts and the boundary conditions (1.5), we deduce that

 supt∈(0,T)∫Ω(ϱ|→v|2+σ|∇b|2+μ(b−lnb))dx+∫Q4ν|D|2+ν1b2|∂tb+div(b→v)|2dxdt≤∫Ω(ϱ|→v0|2+σ|∇b0|2+μ(b0−lnb0))dx.

This natural energy estimate thus motivates the minimal assumptions on the data and and simultaneously indicates the correct choices of function spaces for the problem. In addition, as it will be seen later, the unknown also satisfies a version of minimum/maximum principle and therefore we shall also require that is bounded and uniformly positive on .

Before we formulate the main theorem of the paper, we briefly recall the notations for the relevant function spaces. The standard Lebesgue and Sobolev spaces are denoted by and , respectively, and are equipped with the usual norms denoted by and , respectively. Since we also need to deal with vector- and tensor-valued functions, we will emphasize this character by writing or if necessary. Moreover, since the velocity is a solenoidal function, we define

 L20,div :=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯{→v∈C∞0(Ω;Rd):div→v=0}∥⋅∥2, W1,20,div :=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯{→v∈C∞0(Ω;Rd):div→v=0}∥⋅∥1,2.

Having defined the proper function spaces, we notice that the assumptions on the data, see (1.7), are a natural consequence of (4.4). We can now formulate our main theorem.

###### Theorem 4.1.

Let be a Lipschitz domain and . Assume that , , , and are positive constants and that satisfies (1.7). Then, there exists a couple such that

 (4.5) →v ∈L∞(0,T;L20,div)∩L2(0,T;W1,20,div), (4.6) b ∈L∞(0,T;W1,2(Ω)),b,b−1∈L∞(Q), (4.7) ∂t→v ∈L1(0,T;(W1,20,div∩Wd+1,2(Ω)d)∗), (4.8) ∂tb ∈L1(Q)

and

 (4.9) (∂tb+div(b→v))∈L2(Q),Δb∈L2(Q),

and this couple solves (1.1)–(1.6) in the following sense: for almost all and all we have (recall that denotes the symmetric part of )

 (4.10) ⟨∂t(ϱ→v),→w⟩+∫Ω(2νD−σ∇b⊗∇b−ϱ→v⊗→v):∇→wdx=0;

for almost all and all we have

 (4.11) ∫Ω(ν1(∂tb+div(b→v))b2+μ(1−1b))w+2σ∇b⋅∇wdx=0,

and the initial data are attained in the following sense:

 (4.12) limt→0+(∥→v(t)−→v0∥2+∥b(t)−b0∥1,2)=0.

Moreover, the following energy inequality holds for all :

 (4.13) 12∫Ω(ϱ|→v(t)|2+σ|∇b(t)|2+μ(b(t)−lnb(t)))dx+∫t0∫Ω2ν|D|2+ν12b2|∂tb+div(b→v)|2dxdτ≤12∫Ω(ϱ|→v0|2+σ|∇b0|2+μ(b0−lnb0))dx.

Notice here that (1.1) is automatically satisfied since ; the identity (4.10) represents the weak formulation of (1.2), (1.3) and (1.5), where we have omitted the presence of the pressure by choosing divergence-free test functions. Finally, (4.11) is a weak formulation of (1.4) and (1.5), which is obtained by division of (1.4) by , multiplying the result by a test function , integrating over and in the elliptic term using integration by parts together with (1.5). Observe also that because of the assumed regularity (4.6) and (4.9), we could simply say that (1.4) holds almost everywhere in , but we would lose the information (1.5), which is encoded in (4.11). Nevertheless, (1.4) can be directly obtained from (4.11) by setting with an arbitrary and noting that (here we use the fact that and integration by parts)

 ∫Ω∇b⋅∇(b2φ)dx =∫Ω2b|∇b|2φ+b2∇b⋅∇φdx=∫Ω2b|∇b|2φ−div(b2∇b)φdx =−∫Ωb2Δbφdx.

## 5. Proof of Theorem 4.1

For the sake of simplicity we set the parameters and equal to 2 and , and equal to 1, and we also divide (1.4) by 2 in order to avoid the presence of all constants. Note that such a simplification does not change the assertion of the theorem, as the proof can be straightforwardly extended to any other values of these parameters thanks to their positivity. We will also frequently use the symbol to denote a generic positive constant, whose value may change from line to line, and can depend only on the data. In instances when also depends on the indices of the approximating problem appearing in the proof, this will be clearly indicated.

### 5.1. Galerkin approximation

First, we introduce the following cut-off function

 Tn(s):=min{n,max{n−1,s}} for s∈R.

Next, we find a basis of , which is orthogonal in and a basis of , which is orthonormal in