Weak-strong uniqueness and energy-variational solutions for a class of viscoelastoplastic fluid models

We study a model for a fluid showing viscoelastic and viscoplastic behavior, which describes the flow in terms of the fluid velocity and an internal stress. This stress tensor is transported via the Zaremba--Jaumann rate, and it is subject to two dissipation processes: one induced by a nonsmooth convex potential and one by stress diffusion. We show short-time existence of strong solutions as well as their uniqueness in a class of Leray--Hopf type weak solutions satisfying the tensorial component in the sense of an evolutionary variational inequality. The global-in-time existence of such generalized solutions has been established in a previous work. We further study the limit when stress diffusion vanishes. In this case, the above notion of generalized solutions is no longer suitable, and we introduce the concept of energy-variational solutions, which is based on an inequality for the relative energy. We derive general properties of energy-variational solutions and show their existence by passing to the non-diffusive limit in the relative energy inequality satisfied by generalized solutions for non-zero stress diffusion.


Introduction
In a three-dimensional bounded domain Ω ⊂ R 3 , we consider the flow of an incompressible viscoelastoplastic fluid governed by the equations where (0, T ) is a time interval with T ∈ (0, ∞].Equation (1a) describes the evolution of the fluid flow (with constant density ρ = 1) in terms of the Eulerian velocity v v v : Ω × (0, T ) → R 3 and pressure p : Ω × (0, T ) → R, subject to an external forcing f f f : Ω × (0, T ) → R 3 .The fluid stress T = ηS + 2µ(∇v v v) sym − pI decomposes into the classical term 2µ(∇v v v) sym − pI for Newtonian fluids and an additional stress S : Ω × (0, T ) → R 3×3 .Here (∇v v v) sym = 1 2 (∇v v v + ∇v v v ⊤ ) is the rateof-strain tensor, and µ > 0, η ≥ 0 denote fixed constants.Since p shall describe the physical pressure, that is, the (negative of the) spherical part of the fluid stress, we assume S to be a symmetric deviatoric tensor field, that is, S = S ⊤ and Tr S = 0.By equation (1b), the extra stress S is transported via the Zaremba-Jaumann rate where Moreover, S is subject to a diffusion process induced by the term γ∆S with γ ≥ 0 and an additional nonlinear dissipation due to the subdifferential ∂ P(S) of the nonsmooth convex potential P. The system is completed by no-slip conditions for v v v and homogeneous Neumann conditions for S (in the case γ > 0) on the boundary as well as initial conditions.
The Zaremba-Jaumann derivative (2) describes only one possible choice for the objective stress rate, and there exist other objective tensor derivatives used for the description of stress evolution.However, as we explain below, the model ( 1) is motivated from geodynamics, where the Zaremba-Jaumann derivative is a common choice (see [21,12,13,22]).Moreover, it is crucial in the mathematical analysis of (1) since it (formally) guarantees the identity This property can be used to reveal information on the evolution of the total quadratic energy which consists of the kinetic energy associated with v v v and the stored elastic energy associated with S.More precisely, smooth solutions (v v v, S) to (1) formally satisfy the energy-dissipation balance for all t ∈ [0, T ).This shows that the total energy is dissipated by three processes: the direct fluid viscosity with parameter µ > 0, the nonsmooth stress-dissipation potential, and the stress diffusion with parameter γ ≥ 0.
The main non-standard feature of the system (1) is the occurrence of the set-valued subdifferential ∂ P(S) in (1b), the meaning of which will be specified in the following.Observe that the energy-dissipation balance (3) suggests to examine (1) in an L 2 framework and we thus seek stress tensors S = S(t) taking values in the space L 2 δ (Ω) = {S ∈ L 2 (Ω) 3×3 : S = S ⊤ , TrS = 0}.The dissipation potential P is now defined to be a convex and lower semicontinuous function P : L 2 δ (Ω) → [0, ∞] that satisfies P(O) = 0, where O denotes the zero tensor.By definition, the convex subdifferential ∂ P of P is then given by S → ∂ P(S) := G ∈ L 2 δ (Ω) : P( S) ≥ P(S) + Ω G : ( S − S) dx for all S ∈ L 2 δ (Ω) .
Such nonsmooth dissipation potentials allow to include plastic effects in the model, and related viscoelastoplastic fluid models are used in geodynamics to describe rock deformation in the lithosphere; see [21,12], where Here a > 0 is a constant, and the yield stress σ yield > 0 determines the transition to plastic behavior.One readily verifies that P defined in this way has the above properties.Further examples for possible choices of P can be found in [10].
The beginning of the mathematical analysis of viscoelastic fluid models, also using objective derivatives different from the Zaremba-Jaumann rate (2), can be dated back to the middle 1980s; see [14,25,24,8,23] for example.Since all objective derivatives come along with strong nonlinearities, the first result on global existence of weak solutions was only established several years later by Lions and Masmoudi [20], who studied the system (1) for γ = 0 and a quadratic dissipation potential such that ∂ P(S) = aS for some a ≥ 0. In this case, (1b) becomes a transport equation, and existence can be deduced from the propagation of compactness in L 2 .This tool is no longer available when P is nonlinear and nonsmooth.For such potentials, large-data global existence can be achieved for diffusive regularizations of the tensorial transport equation, as recently demonstrated in [6,2,10].The article [10] considers problem (1) with γ > 0 and proves global existence of generalized solutions compsed of a weak formulation for (1a) and a variational inequality for (1b) (cf.Definition 3.1 below).In the present article, we continue this analysis in two directions.
In a first part, we investigate the case γ > 0, complementing and refining the existence analysis of generalized solutions in [10].Here, our main results are the short-time existence of strong solutions (see Theorem 3.4) as well their uniqueness among generalized solutions (see Theorem 3.5).The presence of the nonsmooth dissipation potential P renders the construction of (local) strong solutions a non-trivial and interesting question.Indeed, some care has to be taken to derive a priori estimates compatible with the nonsmoothness of P, and our construction strongly relies on the fact that γ > 0. Concerning our uniqueness result, we note that the generalized solutions as considered here comprise the family of Leray-Hopf weak solutions of the Navier-Stokes equations (for η = 0), and hence the uniqueness of generalized solutions seems to be out of reach.The main step in the proof of the weak-strong uniqueness principle mentioned above is the derivation of an evolutionary inequality for the relative energy of the form for t ∈ (0, T ), see (38) below, which allows to compare generalized solutions (v v v, S) with (more regular) test functions (ṽ v v, S).Here F contains more or less the (nonlinear) differential operator associated with problem (1) applied to the test function (ṽ v v, S), and W (K ) consists of terms describing the relative dissipation as well as terms arising from the nonlinearities in (1).It further depends on the non-negative function K , which we call the regularity weight since it determines the class of admissible test functions (ṽ v v, S) such that (5) is meaningful.
In the second part we study the problem without stress diffusion, that is, when γ = 0. Our interest in this case stems, among other things, from the original models used in geodynamics [21,12], where stress diffusion does not appear.If γ = 0, problem (1) reduces to the system As mentioned above, stress diffusion serves as a regularization making existence results accessible by means of parabolic theory.When γ = 0, the energy estimate (3) no longer controls the gradient of S in L 2 loc (Ω × [0, T )), and it becomes unclear how to pass to the limit along approximate solutions in the term S(∇v v v) skw −(∇v v v) skw S in (6b) associated with the Zaremba-Jaumann derivative.In this article, we provide a framework allowing us to treat the case γ = 0.It relies on the rather weak notion of energy-variational solutions to (6) (cf.Def.4.1 below), which is based on an inequality for the relative energy R(v v v, S |ṽ v v, S), adapting the aforementioned relative energy inequality (5) for generalized solutions when γ > 0. Under a convexity assumption, we may pass to the limit γ → 0 in this inequality using weak lower semicontinuity arguments.This allows us to construct an energy-variational solution to (6) as the limit of a sequence of generalized solutions (v v v γ , S γ ) to (1) for γ > 0, see Theorem 4.2.
The idea to base a solution concept on a relative energy estimate goes back to Lions [19,Def. 4.1], who introduced the notion of dissipative solutions for the incompressible Euler equations.A dissipative solution does not fulfill the differential equation in a distributional sense, but rather satisfies a relative energy inequality with respect to any sufficiently smooth test function.After the seminal work by Lions, this concept has been adopted in other contexts as well, e.g., in the context of viscous electro-magneto-hydrodynamics [1], liquid crystals [16], and nematic electrolytes [3].
Another generalized solution concept, which is often used in the context of fluid dynamics, is the notion of so-called measure-valued solutions [9].Measure-valued solutions carry more information than dissipative solutions, but this is achieved by increasing the degrees of freedom: In every point in time and space the solution carries an infinite dimensional measure.The expectation of every measure-valued solution fulfills the dissipative formulation [5], which is identified as a desirable quantity in the case of liquid crystals [15].Moreover, the concept of dissipative solution is amenable from a numerical point of view.In the case of anisotropic fluids, a structure-preserving finite element scheme was proven to converge to a dissipative solution, but the convergence to a measure-valued solution seems to be out of reach [3,Rem. 3.7].Especially the high-order regularizations used to prove convergence to a measure-valued solution [16] are not amenable from a numerical standpoint.Moreover, in contrast to the set of weak solutions, the set of dissipative solutions is convex and weakly- * closed in case of convex energy and dissipation potential.This can be used to define selection criteria choosing a certain dissipative solution in order to achieve uniqueness [17], see also Remark 4.4 below.
The article at hand extends previous concepts of dissipative solutions in two directions.
Firstly, we do not fix the regularity weight K in the relative energy estimate (5) defining energy-variational solutions.This allows to consider different classes of test functions and lets us derive different results for different choices of the regularity weight K .For one class, we preserve the weak formulation of the Navier-Stokes-like part in the generalized formulation (see Prop. 4.3), and for another choice, we can deduce the convexity of the solution set (see Rem. 4.4).As mentioned above, this convexity property may be used to select a unique physically relevant solution [18,17].Secondly, we refrain from dropping the term W (K ) in ( 5), which for suitable K is non-negative and usually estimated by 0. By choosing K such that the relevant terms in W (K ) are convex and lower semicontinuous, they can be kept when passing to the limit γ → 0, thus making the solution concept more selective (see Remark 4.2).The structure of this article is as follows.We first introduce the general notation and prepare some auxiliary results in Section 2. In Section 3 we recall the existence theorem for generalized solutions to (1) in the case γ > 0 and show the local-in-time existence of strong solutions.Subsequently, we prove that strong solutions are unique in the class of generalized solutions.For this purpose, we derive a suitable relative energy inequality.In Section 4 we introduce the notion of energy-variational solutions to (1) in the case γ = 0, which is based on a similar relative energy inequality, and we show their existence by an approximation with generalized solutions for γ > 0. Subsequently, we derive general properties of energy-variational solutions.
In the appendix we explain how to infer the present notion of generalized solutions from the slightly different notion introduced in [10].

Preliminaries
In this section we prepare the notation used throughout this article as well as some helpful inequalities.We further introduce the basic regularity hypotheses on the data assumed throughout this manuscript.

Notations
The third-order tensor A ⊗ a a a is defined by (A ⊗ a a a) jkℓ = A jk a a a ℓ .By A ⊤ and Tr A we denote the transpose and the trace of A, and R 3×3 δ := {A ∈ R 3×3 : A ⊤ = A, TrA = 0} denotes the class of symmetric deviatoric matrices.Moreover, 0 0 0 ∈ R 3 and O ∈ R 3×3 denote the zero vector and the zero tensor, respectively.
The symbol Ω always denotes a bounded Lipschitz domain in R 3 , and points in (x,t) ∈ Ω × (0, T ), T > 0, consist of a spatial variable x ∈ Ω and a time variable t ∈ (0, T ).By ∂ t u and ∂ j u := ∂ x j u, j = 1, 2, 3, we denote time and spatial partial derivatives of a (sufficiently regular) function u.We further write ∇ and ∆ for the gradient and the Laplace operator.For a vector , the symmetric and skew-symmetric parts of ∇v v v are denoted by We let • X denote the norm of a Banach space X .When the dimension is clear from the context, we do not distinguish between X and its n-fold Cartesian product X n .We further write X * for the dual space of X , and ϕ, x denotes the duality pairing of ϕ ∈ X * and x ∈ X .When X is a Hilbert space, we sometimes write (x, y) for the inner product of two elements x, y ∈ X .
For an interval I ⊂ R, the class of continuous X -valued functions is denoted by C 0 (I; X ).For the associated Bochner-Lebesgue spaces we write L q (I; X ), and we define W 1,q (I; X ) := {u ∈ L q (I; X ) : ∂ t u ∈ L q (I; X )}.As above, we set H 1 (I; X ) := W 1,2 (I; X ), and the classes L q loc (I; X ) and H 1 loc (I; X ) contain all functions that, when restricted to any compact subinterval J ⊂ I, belong to L q (J; X ) or H 1 (J; X ), respectively.If I = (0, T ), we simply write C 0 (0, T ; X ) = C 0 (I; X ) and L q (0, T ; X ) = L q (I; X ).Moreover, for functions u on Ω × I, we sometimes abbreviate u(t) := u( • ,t) for t ∈ I.
We let C ∞ 0,σ (Ω) := ϕ ∈ C ∞ 0 (Ω) 3 : ∇• ϕ = 0 denote the class of smooth solenoidal vector fields, and the associated Lebesgue and Sobolev spaces are given by where the conditions ∇•v v v = 0 and v v v| ∂ Ω • n n n = 0 have to be understood in a weak sense; see [11, Theorem III.2.3] for example.If I ⊂ R is an interval, we further set Moreover, we introduce the spaces of symmetric deviatoric fields In view of the (formal) energy dissipation law (3), for γ > 0 we seek solutions (v v v, S) of (1) in the natural energy space LH T × X T where Usually, the time T ∈ (0, ∞] will be fixed and we simply write for the solution space.We further need function spaces obeying additionally a Serrin-type regularity criterion, which are defined by and we introduce the space Functionals.Let X be a Banach space and P : X → [0, ∞] be convex and proper (that is, P ≡ ∞) .Then ∂ P denotes the subdifferential of P defined by for u ∈ X .When X is a Hilbert space, then we can identify X * with X such that ∂ P(u) ⊂ X .Moreover, the convex conjugate of P is denoted by On the state space and we let R : Furthermore, given a functional K : Q → [0, ∞], we let

General hypotheses
Throughout this manuscript, we investigate (1) under the following assumptions.
Hypothesis 2.3.We let Ω ⊂ R 3 be a bounded Lipschitz domain and T ∈ (0, ∞].Moreover, P : L 2 δ (Ω) → [0, ∞] denotes a convex, lower semicontinuous functional that satisfies P(O) = 0.In particular, P is weakly lower semicontinuous in L 2 δ (Ω).For the remaining data we assume the regularity 3 Generalized and strong solutions in the case of stress diffusion In this section, we investigate system (1) for γ > 0, that is, when stress diffusion is present.We begin by introducing the notions of strong and generalized solutions, and show that for sufficiently regular data with P(S 0 ) < ∞ a strong solution exists at least locally in time.The global-in-time existence of generalized solutions has been established in [10].The second main result of this section concerns a weak-strong stability estimate for the relative energy between a generalized and a strong solution.An important implication of this estimate is that any generalized solution coincides with the strong solution starting from the same initial data as long as the latter exists (see Theorem 3.5).The basis of the stability estimate is an inequality involving the relative energy between a generalized solution and an arbitrary sufficiently regular competitor taking the role of a test function, see Proposition 3.6.We will take up this energyvariational inequality in Section 4, where it provides a framework for passing to the limit γ → 0.

Definition of generalized and strong solutions
For the definition of generalized solutions, recall the definition of the energy space X from (7).
Definition 3.1 (Generalized solution).We call a couple (v v v, S) a generalized solution of system (1) if (v v v, S) ∈ X = LH T × X T and the following holds true: 1.The velocity field satisfies the weak formulation , and the partial energy inequality is satisfied for almost all t ∈ (0, T ).
2. The extra stress tensor S satisfies the evolutionary variational inequality for all S ∈ Z q T , q ∈ (2, ∞), and a.a.t ∈ (0, T ).In the article [10], existence for (1) was shown based on a slightly weaker form of the notion of generalized solutions than introduced above.The main difference lies in the fact that the first term in (18), that is, the partial relative energy at time t, does not appear in the notion used in [10].The above version including this term is essentially equivalent (cf.Lemma A.1) but better suited for the relative energy methods used in the present paper.Theorem 3.2 (Existence of generalized solutions [10]).Assume that Hypothesis 2.3 is satisfied.Then there exists a generalized solution (v v v, S) ∈ X in the sense of Definition 3.1.
For the most part, Theorem 3.2 was established in [10].Details on how to infer the version above, which is formulated with the current, upgraded version of generalized solutions, are provided in Appendix A. Remark 3.1.Setting S = 0 in (18) yields the partial energy inequality Since we have (∇v v v) sym : S = ∇v v v : S by the symmetry of S, summation of ( 17) and ( 19) further yields the total energy-dissipation inequality where the energy functional E was introduced in (10).
• S ∈ Z q T for some q ∈ (1, ∞), S(0) = S 0 , and for all T ∈ H 1 δ (Ω) and a.a.t ∈ (0, T ) it holds Lemma 2.1 ensures that the first and the second term in the first line and the term in the second line of inequality (21) are in L 1 loc ([0, T )) and thus, in particular, finite a.e. in (0, T ).
Moreover, the velocity field v v v of a strong solution (v v v, S) in the sense of Definition 3.3 is a weak solution to (1a) that satisfies a Serrin condition.This is in accordance with the well-known notion of strong solutions to the classical Navier-Stokes equations as used in [27] for example.In particular, due to this regularity condition, it is not necessary to additionally assume that v v v satisfies the partial energy inequality (17) since the corresponding energy equality is satisfied automatically, that is, As a consequence of these observations, we infer the following consistency property.
Remark 3.2 (Strong solutions are generalized solutions).Any strong solution (v v v, S) in the sense of Definition 3.3 is a generalized solution in the sense of Definition 3.1.For the velocity component, this follows from the fact that strong solutions satisfy the energy equality (22).Inequality (18) for given S ∈ Z q T is obtained upon integrating inequality (21) at time s and with the choice T := S(s) in time from s = 0 to s = t, where one uses the identity Let us finally point out the relation between the variational inequality ( 21) and the differential inclusion (1b).For this purpose let P : . By definition, its convex subdifferential ∂ P maps elements of H 1 δ (Ω) to subsets of H 1 δ (Ω) * , and for almost all t ∈ (0, T ), inequality (21) can be written as a differential inclusion in H 1 δ (Ω) * , namely with the understanding that −γ∆ S, T := γ Ω ∇ S • • • ∇T dx.In Section 3.2 we will construct local-in-time strong solutions enjoying the extra regularity ∆ S(t), ∂ t S(t) ∈ L 2 δ (Ω), which satisfy the inclusion (23) in the L 2 δ (Ω) sense as well (cf.Remark 3.4).Remark 3.3 (Energy equality).We note that for a strong solution (v v v, S), we may recover the energy equality (3).Indeed, the inclusion (23) implies that there exists G with G(t) ∈ ∂ P(S(t)) for a.a.t ∈ (0, T ) such that the equality Thanks to the extra regularity, we may test this identity with S. Integrating the resulting equality in time leads to the partial energy equality for a.e.t ∈ (0, T ), where we used the Fenchel identity G, S = P(S)+ P * (G) for G ∈ ∂ P(S).Furthermore, we may use v v v as a test function in (16), which leads to the partial energy equality (22) for a.a.t ∈ (0, T ).Summing up the two partial energy equalities, we arrive at the energy equality (3).

Local existence of strong solutions
In the following, we show that under an additional regularity hypothesis on the data, problem (1) has a strong solution on a time interval (0, T ) provided T > 0 is small enough.For simplicity, we only consider the case f f f ≡ 0 here, but the result equally holds for forcings Theorem 3.4 (Local existence of strong solutions).Let γ > 0, µ > 0 and η ≥ 0. Suppose that Ω ⊂ R 3 is a bounded domain with C 2 boundary, and let f f f ≡ 0. In addition to the basic hypotheses (15) assume that Then there exists T ∈ (0, ∞) such that problem (1) has at least one strong solution (v v v, S) with initial data (v v v 0 , S 0 ).This solution enjoys the following additional regularity properties: Proof.The main point in the proof is to derive sufficiently strong a priori bounds (cf.( 27) below) for suitably regular solutions of (1) on some short time interval [0, T ] in the case that P is Fréchet differentiable with globally Lipschitz continuous derivative.For the actual construction of a strong solution on (0, T ) × Ω, one may then follow the roadmap in [10]: One first regularizes the potential P : L 2 δ (Ω) → [0, ∞] by its Moreau envelope and then performs a Galerkin approximation for the regularized potentials similar to [10, Section 4].Along the approximate sequence, the a priori estimates can be derived rigorously.Moreover, they allow us to pass to the limit along the approximate solutions (up to a subsequence), not only in the weak formulation for (v v v, S) but also in the strong formulation (21) for the internal stress tensor, and imply the asserted regularity (25).Let us note that for the velocity component, in the Galerkin approximation one should here use a Galerkin basis composed of eigenfunctions of the Stokes operator −P∆ on L 2 σ (Ω) with P denoting the Helmholtz projector P : L 2 (Ω) 3 → L 2 σ (Ω) (rather than the smooth basis functions introduced in [10] for constructing weak solutions).The eigenfunctions of the Stokes operator only belong to H 1 0,σ (Ω) ∩ H 2 (Ω) 3 in general, but this regularity suffices for our purpose.
In the rest of the proof, we will show that, in addition to the fundamental energy estimate with , regular solutions of (1) enjoy the following bound: for Here, the quantity C 27 can be chosen to be non-decreasing in each of its arguments, while the time T (•) > 0 is non-increasing in each of its arguments.
For the derivation of ( 27), we first test the Navier-Stokes equations (1a) with Using Agmon's inequality as well as the Poincaré and Young inequalities and the elliptic estimate ), the integral terms can be estimated as follows: To proceed, let us recall that we may assume without loss of generality that ∂ P(S) is globally Lipschitz continuous.As mentioned above, the rigorous version of this step is to be carried out with the Moreau envelope of the nonsmooth dissipation potential, cf.[10].We may then test the evolution law (1b) for the tensorial component with ∂ t S to find The integral I 3 involving the convective term is estimated, using Agmon's inequality and the bound For I 4 we estimate, using the Gagliardo-Nirenberg-Sobolev inequality, The estimate for I 5 is immediate In combination, we deduce where we used the regularity of P and S and the chain rule to rewrite the term involving the dissipation potential.
Adding up inequalities (29) and (30) gives , we find that the function where Since, by hypothesis, ψ(0) is finite, we may use a classical ODE comparison argument to deduce the existence of a time When P ≡ 0, the evolution law for the tensorial component implies that any strong solution (v v v, S) with the regularity (25) further satisfies S ∈ L 2 (0, T ; H 2 (Ω) 3×3 ).The following remark shows that this conclusion continues to hold true for nonsmooth potentials P given by an integral functional.
Let us first provide the formal argument demonstrating this assertion.The estimate combined with (26) and (27) shows that the term Assuming for the moment that P ∈ C 2 and testing (32) with DP(S) (= G) yields where the first term on the left-hand side is non-negative thanks to the positive semidefiniteness of the Hesse form of the convex function P. We thus deduce the a priori bound which further yields γ∆ S L 2 (0,T ;L 2 (Ω)) ≤ C F L 2 (0,T ;L 2 (Ω)) .
To make the above reasoning rigorous, we may argue as follows.Following the proof of Theorem 3.4, we first construct local-in-time strong solutions using the regularized functions P ε,λ := ρ λ * P ε with λ , ε ∈ (0, 1].Here, P ε denotes the Moreau envelope of the nonsmooth function P, that is, which by construction satisfies the pointwise bound P ε (T) ≤ 1 2ε |T| 2 , while ρ λ denotes a standard (non-negative) mollifier.The functions P ε,λ are C 2 and convex and satisfy the λ -uniform bound 0 whenever λ ∈ (0, 1].(The fact that possibly P ε,λ (0) = 0 does not affect the construction of solutions.)At the level of the strong solutions associated with P ε,λ , estimate (33) and the L 2 bound for γ∆ S can be derived rigorously.The control (34) and the a priori bounds allow to pass to the limit λ → 0 in the strong formulation ( 21) for all T ∈ L 2 δ (Ω), where the term γ Ω ∇S(t) • • • ∇(S(t) − T) dx is to be replaced by −γ Ω ∆S(t) : S(t) − T dx.Finally, in the limit ε → 0, we obtain a strong solution on (0, T ) with the additional regularity (31) satisfying the differential inclusion (1b) in the L 2 δ sense.Since our weak-strong uniqueness principle (Theorem 3.5 below) implies that strong solutions are unique, this shows that already the strong solutions obtained in Theorem 3.4 must have these regularity properties if the dissipation potential is in integral form.

Relative energy inequality and weak-strong stability
The aim of this subsection is to establish the following weak-strong stability result, which compares generalized solutions (v v v, S) with strong solutions (ṽ v v, S) to (1) in terms of the relative energy R(v v v, S|ṽ v v, S) defined in (11).It implies the uniqueness of strong solutions in the class of generalized solutions.Theorem 3.5 (Weak-strong stability).Let (ṽ v v, S) be a strong solution according to Def. 3.3 with initial data (ṽ v v 0 , S0 ).Choose p, q, r, s ∈ (2, ∞) satisfying ( 8) such that ṽ for a.e.t ∈ (0, T ), where for a constant C = C(Ω, p, q, r, s, µ, γ) > 0. In particular, if the initial data (v v v 0 , S 0 ) and (ṽ v v 0 , S0 ) coincide, then ṽ v v ≡ v v v and S = S.
In the above formulation, we call the function K the regularity weight since it measures the minimal regularity of the function (ṽ v v, S) such that both sides of the relative energy inequality (35) remain finite.
We will deduce Theorem 3.5 from a suitable relative energy inequality that compares a generalized solution with any sufficiently smooth function.For its formulation, let us (formally) introduce the operator A = A (1) , A (2) by Note that the interpolation inequality of Lemma 2.1 shows that where X and Z were introduced in ( 7) and ( 9).Now we can state the above-mentioned relative energy inequality.Actually, we derive a family of relative energy inequalities where the regularity weight K is not fixed in advance.Observe that the choice of K influences the class of admissible test functions, which must belong to the set D K defined in (12).Proposition 3.6.Let (v v v, S) ∈ X be a generalized solution according to Definition 3.1, and let K : Q → [0, ∞] be a given functional.Then (v v v, S) fulfills the relative energy inequality

R(v v v(t), S(t)|ṽ v v(t), S(t))
for a.e.t ∈ (0, T ) and all (ṽ v v, S) ∈ D K ∩ Z, where W := W (K ) denotes the relative dissipationlike quantity which depends on the regularity weight K .
Observe that by Lemma 2.1 the quantity W (K ) satisfies W (K ) (v v v, S |ṽ v v, S) ∈ L 1 loc ([0, T )) for all (v v v, S), (ṽ v v, S) as above.
If t 0 P( S) dτ = +∞, the asserted inequality is trivially satisfied.Thus, we may henceforth assume that P( S) ∈ L 1 (0,t).Since (v v v, S) is a generalized solution, it fulfills (17) and (18).Given φ ∈ C([0,t]), we transform these two inequalities according to Lemma 2.2 into their weak formulation in time and add the weak form ( 16) of the Navier-Stokes equations with the test function Φ = −φ ṽ v v (which is admissible as a test function due to a classical approximation argument using Lemma 2.1) to obtain Observing that we can rewrite inequality (40), upon integration by parts in time, as where we used the canceling of several terms (notably of the coupling term Ω η(∇v v v) sym : S dx = Ω η∇v v v : S dx), the fact that v v v and ṽ v v are solenoidal vector fields, and the antisymmetry of (∇ṽ v v) skw .Choosing φ (τ) = ϕ(τ)e − τ 0 K (ṽ v v, S) d τ and invoking Lemma 2.2 yields the asserted inequality (38).
For the derivation of ( 35) from ( 38) it suffices to show that for a.e.t ∈ (0, T ) we have for K given by (36).To verify inequality (42), we first show, for a.e.t ∈ (0, T ), the two relations Note that, as a consequence of the weak solution property of (ṽ v v, S) and the extra regularity of ṽ v v, for a.e.t ∈ (0, T ) we have the equality which immediately gives (43).Inequality (44) follows upon choosing T = S(t) ∈ H 1 δ (Ω) as a test function in the strong variational inequality (21) satisfied by (ṽ v v, S).We are thus left to prove that if K is given by (36) for some sufficiently large constant C < ∞.Using Lemma 2.1, we may estimate where (r, s) and (p, q) fulfill ( 8).Similarly, we estimate the terms stemming from the Zaremba-Jaumann rate by Hölder's, Gagliardo-Nirenberg's, and Young's inequality as where α = 3/p.Hence, (42) holds for K as in (36) with C = max{C 1 , 2C 2 }.

Energy-variational solutions for vanishing stress diffusion
In this section, we investigate problem (6), that is, (1) for γ = 0, where stress diffusion is not present.In this case, the basic energy estimates no longer provide weak sequential compactness in L 1 loc (Ω × [0, T )) for the nonlinear term S(∇v v v) skw − (∇v v v) skw S coming from the Zaremba-Jaumann derivative.Thus, we have to further extend the concept of generalized solutions introduced in Def.3.1.Here, we use a concept in the spirit of [15,17], which is reminiscent of the dissipative solutions introduced by P.-L.Lions in [19].In the literature, the term dissipative is, however, employed for various kinds of paradigms, and here we prefer the term energy-variational since the basis of our solution concept is a relative energy inequality (similar to (38)), which can be interpreted as a variation of the energy-dissipation mechanism with respect to functions in the domain D K of the regularity weight K .We show existence of energy-variational solutions and derive some important properties inherent in this solution concept.

The concept of energy-variational solutions
We introduce the notation where As explained above, the notion of energy-variational solutions to ( 6) is based on an adaptation of the relative energy inequality (38) to the case γ = 0. To emphasize the γ-dependence in the quantities appearing in (38), we set for γ ≥ 0, and W γ = W (K ) γ = W (K ) for γ > 0, where A (2) and W (K ) are given in (37b) and (39), respectively.Observe that in the limiting case γ = 0 the quantity W (K ) is not well defined by (39) for general (v v v, S) ∈ X 0 and (ṽ v v, S) ∈ Z 0 since ∇S appears in the advective term.Therefore, we formally integrate by parts and define We now define energy-variational solutions by imposing a relative energy inequality analogous to (38).Recall that the regularity weight K was not fixed in (38), which is reflected in the following solution concept.We emphasize though that the definition to follow is only meaningful for sufficiently "nice" weights K .
is fulfilled for a.e.t ∈ (0, T ) and all (ṽ v v, S) ∈ D K ∩ Z 0 .As before, we call the function K the regularity weight.
Remark 4.1 (Choice of the regularity weight).While the assumption K (0 0 0, O) = 0 is not needed to define such a solution concept, it is a natural hypothesis which ensures the classical energy inequality, i.e., inequality (20) with γ = 0. Note that this inequality directly follows by letting (ṽ v v, S) ≡ (0 0 0, O) in (47), which is admissible when K (0 0 0, O) = 0.This in turn guarantees that P(S) ∈ L 1 loc ([0, T )) and hence, that the terms in (47) are well defined for all (ṽ v v, S) ∈ D K ∩ Z 0 .Observe that standard interpolation estimates show that for (v v v, S) ∈ X 0 and (ṽ v v, S) ∈ Z 0 ∩ D K , see also the estimates in the proof of Theorem 4.2 below.Notice that the classical energy inequality encodes some basic information concerning the long-time behavior of solutions, cf.[15, Thm.6.2].Remark 4.2 (Comparison to previous formulations).The above formulation differs from previous notions of dissipative solutions (cf.[19]), which were based on similar relative energy inequalities.In those definitions, the regularity weight K was fixed and chosen such that the last three lines of (46) are non-negative.Additionally, these terms are usually just estimated from below by zero.By keeping the mentioned terms in the inequality (47), the present solution concept is more selective than previous versions.

Global existence of energy-variational solutions
In Proposition 3.6 we have seen that the generalized solutions for γ > 0 satisfy the relative energy inequality (38).For suitable regularity weights K we may perform the limit γ → 0 in this inequality to obtain the following existence result.Theorem 4.2 (Existence of energy-variational solutions).Under the assumptions stated in Hypothesis 2.3, there exists an energy-variational solution to the initial-boundary value problem (6) of type K given by for a suitable constant C = C(Ω, µ) > 0.
Proof.In order to prove the existence of energy-variational solutions, for each γ > 0 we consider a generalized solution (v v v γ , S γ ) ∈ X to (1), which exists due to Theorem 3.2.By Remark 3.1 the generalized solution (v v v γ , S γ ) fulfills the energy inequality for a.e.t ∈ (0, T ).Additionally, we may find a γ-independent bound on the time derivative of the velocity field.Indeed, since (v v v γ , S γ ) satisfies the weak formulation (16), for every Φ ∈ C ∞ 0,σ (Ω × (0, T )) we can estimate By the aforementioned bounds and a classical diagonalization argument (see [27, Proof of Theorem 3.1.1]for example), we infer the existence of a limit function (v v v, S) such that for each t ∈ (0, T ) we have for an appropriately chosen subsequence γ → 0. Next we conclude that via these convergences, we may pass to the limit in the relative energy inequality (38).Proposition 3.6 guarantees that (38) holds for all non-negative functionals K : Q → [0, ∞], all (ṽ v v, S) ∈ D K ∩ Z and a.e.t ∈ (0, T ).We may reformulate this inequality using Lemma 2.2 and observing that W γ ≥ W 0 .This leads to for all φ ∈ C([0, T ]).We now show that, when choosing K as in (48), we may pass to the limit γ → 0 in each of the terms on the left-hand side of (51) via the convergence properties from (50).Since φ ′ ≤ 0, the first term in ( 51) is convex and continuous on L 2 (Ω × (0,t)) as a function of (v v v γ , S γ ).Hence, it is weakly lower semicontinuous, and its convergence follows from (50c) and (50d).
In order to pass to the limit in the second integral term in (51), we split where The convergence of the part of the integral involving the term C (v v v γ , ṽ v v) follows as usual since the strong convergence (50c) combined with (50a) implies the weak convergence of v v v γ • ∇v v v γ .Concerning the remaining part involving Q(v v v γ , S γ |ṽ v v, S), we assert that (for fixed ṽ is chosen large enough.Since Q is the sum of a quadratic form Q 1 and an affine part, it suffices to show that Q ≥ 0, as this implies the non-negativity of Q 1 and hence its convexity (cf.[4,Prop. 3.71]).To show the non-negativity of Q (and its continuity) we use the Sobolev-Poincaré inequality ṽ v v L 6 (Ω) ≤ C ∇ṽ v v L 2 (Ω) and estimate, choosing C sufficiently large, Combined with the bound this shows that Q is non-negative and continuous.Hence, we may use weak lower semicontinuity to take the limit γ → 0 in the integral term t 0 φ Q(v v v γ , S γ |ṽ v v, S)e t s K (ṽ v v, S) dτ ds.For the third term in (51), note that the induced functional T → t 0 φ P(T)e t s K (ṽ v v, S) dτ ds is convex and, by Fatou's lemma, lower semicontinuous on L 2 (0,t; L 2 (Ω)).This implies its weak lower semicontinuity, whence lim inf For the term involving the operator A γ , we first estimate the γ-dependent term as for all t ∈ (0, T ).Due to the bound (49) on t 0 γ ∇ S γ 2 L 2 (Ω) ds and the regularity of S, the right-hand side of the last inequality tends to zero as γ → 0. Since the remaining terms do not explicitly depend on γ and (v v v γ , S γ ) occurs at most linearly, we conclude convergence of the last term on the left-hand side of (51).
In total, we can pass to the limit inferior with γ → 0 in (51), which implies inequality (47) by invoking Lemma 2.2.Therefore, (v v v, S) is an energy-variational solution to (6).
Since K (ṽ v v, O) = 0 for all ṽ v v ∈ C ∞ 0,σ (Ω), we further conclude from Proposition 4.3 below that (6a) is satisfied in the weak sense.Remark 4.3.An argument similar to the proof of Theorem 4.2 shows that the set of energyvariational solutions of type K given by ( 48) is weakly- * sequentially compact in X 0 in the following sense: For every sequence of energy-variational solutions (v v v j , S j ) there exists an energy-variational solution (v v v, S) ∈ X 0 and a subsequence (also denoted by for all t ∈ (0, T ).Indeed, a sequence (v v v j , S j ) of energy-variational solutions of type K satisfies the relative energy inequality (47) with (ṽ v v, S) = (0 0 0, O), that is, the energy inequality for a.a.t ∈ (0, T ).Moreover, Proposition 4.3 below ensures that the Navier-Stokes component holds in the weak sense.In the same way as above, we further deduce that (∂ t v v v j ) is bounded in L 4/3 (0,t; (H 1 0,σ ) * ) for all t ∈ (0, T ), and we conclude the asserted convergence properties as well as v v v j → v v v in L 2 (0,t; L 2 (Ω)) for all t ∈ (0, T ).In nearly the same way as for the limit γ → 0 above, we can now perform the limit j → ∞ to infer that (v v v, S) is an energy-variational solution of type K .

Further properties
Here we collect general properties of energy-variational solutions.Most importantly, we show that energy-variational solutions are subject to the weak formulation of the Navier-Stokes equations (6a) if the regularity weight K vanishes on C ∞ 0,σ (Ω) × {O}, where O ∈ L 2 δ (Ω) denotes the zero tensor.Thus, this is in particular the case for the solutions obtained in Theorem 4.2.Notice that this property is to be expected for problem (6) since the lack of weak sequential compactness in L 1 loc (Ω × [0, T )) only occurs in the tensor component (6b) but not in (6a).Proposition 4.3.Suppose that the regularity weight Then every energy-variational solution of type K is a weak solution of the Navier-Stokes component in the sense of (16).
Proof.The hypothesis on K implies the identity e for all ṽ v v ∈ C ∞ 0,σ (Ω × [0, T )) and for a.e.t ∈ (0, T ).Applying Lemma 2.2 and observing that 0 (ṽ v v, O), S = − Ω η(∇ṽ v v) sym : S dx, we find We now integrate by parts in time the second term in the last line of (52), regroup terms, let ṽ ) and α > 0, and multiply the inequality by 1/α to infer In the limit α → ∞ the two terms with the prefactor 1 α disappear.Thus, appealing to Lemma 2.2 we deduce where we used the fact that ũ u u(•, T ) = 0. Since the left-hand side of the last inequality is a linear functional of ũ u u with the latter being allowed to vary in the linear space C ∞ 0,σ (Ω × [0, T )), we infer the asserted equation (16).
Another quite natural property is that energy-variational solutions are monotonic in the type in the following manner.Proof.This can readily be seen by expressing (47) in a weak form with Lemma 2.2, using the test function φ (s) = ϕ(s)e − s 0 (L (ṽ v v, S)−K (ṽ v v, S)) dτ for ϕ ∈ W ((0, T )), and using Lemma 2.2 again to return to a pointwise version of (47) with K replaced with L .
for s ∈ (2, ∞) and r ∈ (3, ∞) such that 2/s + 3/r = 1.Similarly to (45), we observe that W 0 (v v v, S |ṽ v v, S) ≥ 0 for all (v v v, S) ∈ X 0 and all (ṽ v v, S) ) and non-negative, it is convex in (v v v, S).This implies that the solution set of energy-variational solutions of type K s is convex.Additionally, the solution set is weakly- * closed (cf.Remark 4.3) and bounded and therewith compact in the weak- * topology of X 0 .This convexity and compactness properties may be used to select an energy-variational solution with maximal dissipation [18].Such a selected maximally dissipative solution can not only be argued to be physically more reasonable since they minimize the energy along all energy-variational solutions, but they also turn out to be analytically favorable since the solution concept is well posed (see [17,18]).

Energy-variational solutions versus strong solutions
In this subsection, we prove two results attempting to further justify the concept of energyvariational solutions.First, in Proposition 4.6 we show that any hypothetical strong solution is unique in the class of energy-variational solutions.Second, in Proposition 4.8 we prove that any energy-variational solution enjoying some additional regularity is a (unique) strong solution.Observe that both results only hold for strong solutions that belong to D K .In this sense, the regularity weight K also determines the minimal regularity of a strong solution for comparison with an energy-variational solution.
We show next that if there exists an energy-variational solution of type K and a strong solution in the class D K , both emanating from the same initial data, then these solutions coincide.This statement and its proof are parallel to the one given in Theorem 3.5 for γ > 0. In particular, we also derive a corresponding weak-strong stability estimate.Proposition 4.6 (Weak-strong uniqueness).Assume that K is a regularity weight such that W 0 = W (K ) 0 defined in ( 46) is non-negative, i.e., W 0 (v v v, S |ṽ v v, S) ≥ 0 for all (v v v, S) ∈ X 0 and (ṽ v v, S) ∈ Z 0 .Let (v v v, S) ∈ X 0 be an energy-variational solution of type K with initial data (v v v 0 , S 0 ), and let (ṽ v v, S) ∈ Z 0 ∩ D K be a strong solution in the sense of Def.4.5 with initial data (ṽ v v 0 , S0 ).Then the inequality holds.Especially, if the initial conditions coincide, i.e., (v v v 0 , S 0 ) = (ṽ v v 0 , S0 ), then every energyvariational solution of type K coincides with the (hypothetical) strong solution (ṽ v v, S).
Proof.From the assumption that W 0 is non-negative, the inequality (47) also holds without this term.Similarly to the proof of Theorem 3.5, we observe that A (1) (ṽ v v, S),v v v − ṽ v v = 0, P(S) − P( S) + A (2) 0 (ṽ v v, S), S − S ≥ 0 for a.e.t ∈ (0, T ).Hence, all terms in the second line of (47) may be estimated from below by zero, which implies (55).
Corollary 4.7.The energy-variational solution from Theorem 4.2 also fulfills the weak-strong uniqueneness property, that is, if there exists a strong solution in the sense of Def.4.5 it coincides with the energy-variational solution from Theorem 4.2.
Proof.According to Proposition 4.4, an energy-variational solution of type K given in (48) is also an energy-variational solution solution of type K s given in (53), which fulfills the assumptions of Proposition 4.6.This ensures that the weak-strong uniqueness also holds for the smaller set of energy-variational solutions of Theorem 4.2, as long as the strong solution belongs to D K s for some s ∈ (2, ∞).Clearly, this is the case for all strong solutions in the sense of Def.4.5.
To infer that sufficiently regular energy-variational solutions are already strong solutions, we assume that P is given in integral form P(S) = Ω P(S) dx.This leads to approximation properties that are sufficient to show that energy-variational solutions with more regularity, that Inserting this into (47), multiplying the resulting relation by 1/α, and sorting the different terms according to the appearing exponent of α, we end up with αR(v v v, S,u u u, T, α) + Passing to the limit α → 0 in (59), we infer from the boundedness of all terms in R and the continuity property lim α → 0 K (v v v + αu u u, S +αT) = K (v v v, S), which holds for K = K s given in (53), that where the regularity of (v v v, S) ∈ D K ∩ Z 0 guarantees that G ∈ L 1 (0,t; L 2 δ (Ω)).By choosing T = Te − t • K (v v v,S) dτ for (0 0 0, T) ∈ D K ∩ Z 0 and employing inequality (58), we infer from (61) that as ε → 0. In particular, the strong convergence statement implies v v v ε (t) L 2 (Ω) → v v v(t) L 2 (Ω) for almost every t ∈ (0, T ), at least for a suitable subsequence.This allows to pass to the limit ε → 0 in (17) (with (v v v, S) replaced with (v v v ε , S ε )), where for the last term the identity Ω ηS : ∇v v v dx = − Ω η(∇• S) •v v vdx can be used.In the end, (v v v, S) also satisfies (17).
Finally, observe that in [10] the boundary ∂ Ω was assumed to be of class C 1,1 , which allowed for the construction of a suitable extension of the boundary conditions.Since we consider homogeneous boundary conditions (1c) here, we can use the trivial extension, whence Lipschitz boundary is sufficient.
The following lemma shows that inequality (64) implies the stronger inequality (18), which also takes into account the values of S and S at time t.
Proof.The assertion is obtained similarly as in the proof of [10,Prop. 3.3].We therefore only sketch the argument.Extending S by zero for t < 0, we define for κ > 0 T is arbitrary.We observe that Sκ,δ ∈ Z q T satisfies Sκ,δ (0) = S(0), Sκ,δ (t) = S κ (t) and, as κ, δ ↓ 0, the sequence { Sκ,δ } approximates S (in a suitable sense).To infer ineq.(18), we insert Sκ,δ as a test function in ineq. (64),take the limit δ ↓ 0 and then send κ ↓ 0. Let us only point out how to pass to these limits in the terms involving a time derivative, since the remaining integrals can be handled as in [10,Prop. 3.3] where the last term in the second equation vanishes as δ → 0. Thus, combining these two identities and treating the remaining terms as in the proof of [10,Prop. 3.3], we deduce for General notations.If a a a = (a a a j ), b b b = (b b b j ) ∈ R 3 are two vectors, their inner product and their tensor product are denoted by a a a • b b b = a a a j b b b j and a a a ⊗ b b b with (a a a ⊗ b b b) jk = a a a j b b b k , respectively.Here and in what follows, we tacitly use Einstein's summation convention.For the inner product of two second-order tensors A = (A jk ), B = (B jk ) ∈ R 3×3 , we write A : B = A jk B jk , and the inner product of two third-order tensors C

Remark 3 .
4 (H 2 regularity for S).Suppose that the dissipation potential P takes the form of an integral P(S) = Ω P(S(x)) dx for a lower semicontinuous, convex function P : R 3×3 δ → [0, ∞] with P(O) = 0.Then, under the hypotheses of Theorem 3.4, strong solutions (v v v, S) satisfying (25) enjoy the additional regularity

Proposition 4 . 4 (
Monotonicity).If (v v v, S) is an energy-variational solution of type K , and if K (ṽ v v, S) ≤ L (ṽ v v, S) for all (ṽ v v, S) ∈ D L ⊂ D K and a.a.t ∈ (0, T ), then (v v v, S) is an energyvariational solution of type L .

Remark 4 . 4 (
Convex solution set).By Proposition 4.4, the energy-variational solution established in Theorem 4.2 is also an energy-variational solution of type K s given by