On Unsteady Internal Flows of Incompressible Fluids Characterized by Implicit Constitutive Equations in the Bulk and on the Boundary

Long-time and large-data existence of weak solutions for initial- and boundary-value problems concerning three-dimensional flows of incompressible fluids is nowadays available not only for Navier–Stokes fluids but also for various fluid models where the relation between the Cauchy stress tensor and the symmetric part of the velocity gradient is nonlinear. The majority of such studies however concerns models where such a dependence is explicit (the stress is a function of the velocity gradient), which makes the class of studied models unduly restrictive. The same concerns boundary conditions, or more precisely the slipping mechanisms on the boundary, where the no-slip is still the most preferred condition considered in the literature. Our main objective is to develop a robust mathematical theory for unsteady internal flows of implicitly constituted incompressible fluids with implicit relations between the tangential projections of the velocity and the normal traction on the boundary. The theory covers numerous rheological models used in chemistry, biorheology, polymer and food industry as well as in geomechanics. It also includes, as special cases, nonlinear slip as well as stick–slip boundary conditions. Unlike earlier studies, the conditions characterizing admissible classes of constitutive equations are expressed by means of tools of elementary calculus. In addition, a fully constructive proof (approximation scheme) is incorporated. Finally, we focus on the question of uniqueness of such weak solutions.


Dedication
This article is written as a contribution to the celebration of the 100th anniversary of the birth of Olga Aleksandrovna Ladyzhenskaya (March 7, 1922 -January 12,2004) and to honor her scientific achievements.
Olga Ladyzhenskya seems to have been the first to initiate interest in the mathematical community to study incompressible fluid models that go beyond the Navier-Stokes equations.At the International Congress of Mathematicians in Moscow 1966, she presented arguments (see [73] and [75]), based on the kinetic formulation 1 , indicating that the macroscopic relation between the stress and the symmetric part of the velocity gradient should be polynomial.Ladyzhenskaya's model belongs to the class of power-law fluid models (sometimes also called modified or generalized Navier-Stokes fluids) characterized by a power-law index r, where the value r = 2 corresponds to the Navier-Stokes fluid.Ladyzhenskaya was interested in the rigorous analysis of models with r > 2: she has found that for those with r ≥ 5  2 one could prove that the weak solution corresponding to the relevant initial-and boundary-value problem (in the sense of Leray and Hopf [78,64]) not only exists for long-time and large-data but it is unique.This uniqueness result should be contrasted with her counterexample to uniqueness of a weak solution of the Navier-Stokes equations in special time-dependent domains, see [74].She also addressed, particularly in her subsequent studies, other aspects of weak solutions of these equations such as higher temporal and spatial differentiability and long-time behavior (the existence of a global attractor and estimates of its dimension).K. R. Rajagopal together with the second author of this study reviewed Olga Ladyzhenskaya's foundational results concerning the analysis of modified Navier-Stokes equations 2 , achieved during the period 1967-2003, in the second part of their handbook article [86].One of the objectives of this study is to provide a brief review of the results obtained in the mathematical analysis of fluids with nonlinear algebraic relation between the Cauchy stress and the velocity gradient achieved after 2003.The main objective is however to present a novel existence theory.

Formulation of the problem and of the main result
Materials are incompressible if the volume of any measurable subpart of the body remains unchanged during a deformation process.For fluids flowing in a d-dimensional domain 3 Ω, the condition of incompressibility expressed in the terms of the velocity v = (v 1 , . . ., v d ) takes the form (2.1) div v(t, x) = 0 for all t ≥ 0 and x ∈ Ω.
Incompressibility, which should be considered as a useful idealization, implies that the Cauchy stress tensor T is of the form (2.2) where only the part S can be determined experimentally.Homogeneous incompressible fluids are characterized as incompressible fluids in which the density remains unchanged and is equal to a positive constant ρ * .Such fluids automatically fulfil the balance of mass equation.
To conclude, setting Q = (0, T ) × Ω and Γ := (0, T ) × ∂Ω with T > 0, the governing equations for unsteady flows of any homogeneous incompressible fluid flowing in a fixed domain Ω with no outflows and inflows and with initial velocity v 0 take the form Here b stands for the density of external body forces.The second equation in (2.3) comes from the balance of linear momentum once (2.2) is incorporated.The third equation says that the tensor S is symmetric; this implies that the balance of angular momentum is fulfilled.The fourth equation states that all considered flows are internal : the fluid cannot enter or leave Ω.
The system (2.3) is incomplete as, in Q, we have d + 1 + d(d + 1)/2 unknowns v, p and S, but merely d + 1 equations.Also on Γ, we only have one scalar equation, but in fact d boundary conditions are expected.
Taking the scalar product of the second equation in (2.3) and v, integrating the result over Ω and using the remaining equations in (2.3), one obtains, after the integration over (0, t) for any t ∈ (0, T ], the energy identity in the form (see [83,Sect. 4.6] for details) 2 See also chapters in the book by J.-L.Lions [79]. 3 Throughout the whole study, the term domain stands for an open bounded connected set in R d .
where s stands for the projection of the normal traction T n to the tangent plane, i.e., s := −(T n) τ , where z τ := z − (z • n)n, n : ∂Ω → R d being the outer normal to ∂Ω.Note that (T n) τ = (Sn) τ .
The second and third terms on the left-hand side represent two independent dissipation mechanisms: the former is associated with the internal friction inside the fluid, the latter corresponds to the interaction of the flowing fluid with (the inner part of) the boundary.Both terms should be, in accordance with the second law of thermodynamics, non-negative.For the Euler fluid, when S = 0 and also s = 0, both terms vanish.The first term will also vanish if all admissible flows are rigid (i.e.Dv = 0), while the second term is equal to zero if all considered flows are subject to the no slip boundary condition, i.e. v τ = 0 on Γ.For the Navier-Stokes fluid characterized by the constitutive equation (2.5) S = 2ν * Dv, where ν * > 0, we conclude that 4   (2.8) This in conjuction with the energy identity (2.4) guarantees control of ∇v in L 2 (Q), which happens to be a key piece of information to establish long-time existence of a weak solution for any domain Ω, . This is what Leray and Hopf proved, see [78] and [64].
The above constitutive equations (2.5) and (2.6) are linear.There are however many fluids (as is also illustrated in more detail in the next section) exhibiting nonlinear relationship between S and Dv.The same concerns slipping boundary conditions.Following Rajagopal [101,102,103], it is tempting to include all these equations under the umbrella of implicit equations relating S and Dv on the one hand and s and v τ on the other hand.Hence, as nonlinear generalizations of (2.5) and (2.6), we add to the problem (2.3) the following equations where G : R d×d × R d×d → R d×d and g : R d × R d → R d are given continuous functions.Adding (2.9) to the first three equations in (2.3) we obtain a closed system of partial differential equations consisting of d + 1 + d(d + 1)/2 equations for d + 1 + d(d + 1)/2 unknowns v, p and S. Adding (2.10) to the fourth equation in (2.3) we get d equations on the boundary: one in the normal direction and (d − 1) of them are formulated at the tangent plane to ∂Ω.Stated differently, the problem (2.3) together with (2.9) and (2.10) is well-formulated.Motivated by Leray-Hopf's theory for the Navier-Stokes equations, it is natural to ask: Can we formulate conditions on G and g that would allow us to establish the long-time and large-data existence of a weak solution for (2.3), (2.9) and (2.10)?If so, can these conditions be formulated in terms of the tools of elementary calculus so that they are accessible to a broad scientific community?This study provides positive answers to these questions.Before formulating the main result, we give the admissibility conditions on G and g.Here, we closely follow our preceding study [20] focused however on a simpler problem (a "mixed" formulation for problems of parabolic type).Regarding the tensorial function G which determines the material response inside the domain (G3) one of the following holds: . In (G2), we used the following notation for (G Further, A T denotes the transpose tensor to A, i.e. (A T ) ij kl = A kl ij and AB T is the standard tensor multiplication, i.e. (AB T ) ij kl = d m,n=1 A ij mn B kl mn .Also, for any tensor A ∈ R d×d × R d×d , the expression A ≥ 0 means that for any X ∈ R d×d there holds AX : X ≥ 0 (which, written in terms of components, is In addition, if we write A > 0 then we mean that the above inequality is strict for all X = 0. Before formulating similar conditions on g, some comments are in order.First, note that the constitutive equation G(S, Dv) = 0 can be replaced by −G(S, Dv) = 0. Then all inequalities in (G2) and (G3) have the opposite signs except the last inequality in (G2).This ambiguity could be fixed for example by requiring that G is such that the first condition in (G2) holds (compare it with the special case of the Navier-Stokes fluids, see (2.5), when one would consider G(S, Dv) = S − 2ν * Dv and not G(S, Dv) = 2ν * Dv − S).Second, as the null points of G are of our interest, we can require the validity of (G2) only in the neighbourhood G(S, Dv) = 0. Also, the Lipschitz continuity in (G1) is required only to guarantee the existence of the partial derivatives in (G2) almost everywhere.Alternatively, one can assume merely the continuity of G in (G1) and substitute (G2) by See also the statement of Lemma 4.2 below.Finally, the conditions (G1)-(G4) are formulated using elementary tools of calculus (limes superior, partial derivatives).As proved in [20], these conditions are equivalent to the statement that the set of null points of G is a maximal monotone r-coercive graph that passes through the origin (see [20, Definition 3.1 and Lemma 3.2] for details.
Regarding the vectorial function g which determines the relation between the shear stress s and the tangential velocity v τ on the boundary, we assume (in a similar way as above) that (g1 (g3) one of the following conditions holds: (g4) there exist c 1 , c 2 > 0 such that, for all (s, v) ∈ R d × R d fulfilling g(s, v) = 0, the following condition holds: The above comments related to G are applicable to g as well.Now, we are ready to formulate our main result (in a vague way): For arbitrary Ω, T , v 0 and b, and for any G and g fulfilling (G1)-(G4) with r > 2d d+2 and (g1)-(g4) with q > 1, there exists a weak solution to the problem (2.3), (2.9) and (2.10).
Overview of the existence theory for G(S, Dv) = 0 after 2003.We restrict our discussion to the most interesting case d = 3.In the first period, prior 2011, the focus of research were models of the type S = (1 + |Dv| 2 ) r−2 2 Dv following the goal to establish global-in-time existence theory for large data for models with low values of r. (Note that the local-in-time existence of smooth solutions to models of power-law type is addressed in [9].)Let us recall that the problems studied by Olga Ladyzhenskaya concerned the subcritical regime when the velocity itself is an admissible test function in the weak formulation of the balance of linear momentum.This corresponds to the case when r ≥ 11  5 .The method of Lipschitz truncation developed in [53] (see also [44]) for time-independent (stationary) problems covers the case r > 6  5 but was left open for the evolutionary case.For the evolutionary case, the "best" result known around the year 2005, see [52], covered the case r ≥ 8  5 using the L ∞ -truncation technique; the approach is restricted to the spatially periodic problem.The extension for flows in general bounded domains subject to Navier's slip boundary was established in [21], while the no-slip boundary conditions were successfully treated in [120], still for r ≥ 8  5 .Finally, Diening, Růžička and Wolf, see [45], inspired by the works of Kinnunen and Lewis [70], extend the method of Lipschitz truncation to the evolutionary case and proved the existence of a weak solution to the evolutionary problem with no-slip boundary condition for r > 6  5 .The remaining case r ∈ [1, 2d d+2 ] is covered by two somehow contradictory recent results, see [2] and [27], which can be interpreted in the way that the range of possible r's in Theorem 5.1 is optimal.In fact, Abbatiello and Feireisl [2] introduce a novel generalized concept of solution (dissipative solution) and establish its existence theory for r ∈ (1, 2d d+2 ].The theory is developed for a smaller class of possible constitutive relations than considered here.More importantly, their concept of solution does not imply either the validity of weak formulation of balance of linear momentum (see (5.5) in Sect.5) or the validity of G(S, Dv) = 0 almost everywhere in Q.On the other hand, in [27], Burczak, Modena and Székelyhidi show the non-uniqueness of even Leray-Hopf solutions for r < 2d d+2 .From this perspective, it seems also to be reasonable to consider only the case r ≥ 2d d+2 .Note that even for this range of r's Burczak et al. [27] prove the result concerning non-uniqueness of a very weak5 solution.As we are dealing with weak solutions, their result is not applicable to our setting.
Inspired by the foundational works on implicit constitutive relations, see [101,102], the question to develop a robust theory covering the whole class of implicitly constituted incompressible fluids arose.Following initial attempts (see [82,13]), a successful theory covering both polynomial and activated fluids was established in [14,12], even in a broader context than considered here: the r-coercivity condition is generalized in terms of Young's functions in the setting of Orlicz spaces and the constitutive equation (2.9) was allowed to vary with time and space, i.e.G(t, x, S, Dv) = 0 in Q.As discussed in length in [20, pp 2048-2049], there are two shortcomings of the results proved in [14,12] (a non-constructive proof and an a priori assumption concerning the existence of a Borel measurable selection).These shortcomings motivated the development of an alternative approach, see [20].The extension of the approach developed in [20] for problems of parabolic type to problems involving flows of incompressible fluids is one of the main objections of this study.
Note that the assumption (g4) eliminates no-slip and perfect-slip boundary conditions from the analysis presented here.However, as indicated in the above discussion of available results one can incorporate both conditions into the analysis.The case of s = 0 on Γ is in fact easy as the boundary term just vanishes.For no-slip boundary conditions, one needs to change the function space for the velocity and pay attention to differences associated with the reconstruction of the pressure (see [120] and [8] for details).
Structure of the paper.In Sect.3, we illustrate how rich the classes of fluids under consideration are by providing a list of models used in various areas of science (completed by a list of references).In Sect.4, we introduce, in a constructive way, ε-approximations of the constitutive equations and provide a summary of their properties (proved in [20]).In particular, at this approximate level the term − div S leads to Lipschitz continuous uniformly monotone elliptic operator (the nicest one can wish to deal with).After introducing basic function spaces in Section 5 we give a precise formulation of the main theorem including also the precise definition of weak solution to (2.3), (2.9) and (2.10).Here we also recall properties of the Lipschitz approximations of Bochner functions needed in the proof of the main theorem.This forms the content of Sect.6.As this article aims at surveying the results in the field, we give, in Sect.7, a summary of the results that concern similar problems including an additional component that makes the whole problem more complicated.Finally, we comment on the available uniqueness results regarding the studied problem in Sect.8.

Examples of implicit constitutive equations
The purpose of this section is to provide an illustrative list of models and boundary conditions covered by the implicit equations (2.9) and (2.10).The aim is to show that these classes of fluids and boundary conditions are rich and particular models appear in various areas of science and engineering.We first focus on the constitutive equations in the bulk, then we discuss the boundary conditions.
These models fit to the setting characterized by the form As the fluid is incompressible, and consequently the trace of Dv vanishes, one observes that within the class (3.1) one has p = − 1 3 tr T .
In Table 1, we distinguish two special subclasses of (3.1), namely S = 2ν(|Dv| 2 )Dv and S = 2ν(|S| 2 )Dv.The simplest deviation from the Navier-Stokes fluid model represents the power law model that can be described in two equivalent ways as follows: where r ∈ (1, ∞) and ν 0 > 0. Referring to Table 1, we thus observe that the same model is called Ostwald-de Waele's model in chemistry, while it is named Glen's model in geomechanics.Denoting r ′ := r/(r − 1) we also have It also serves as the main motivation for the (r, r ′ )-coercivity assumption (G4).All the models listed in Table 1 describe, for suitable range of parameters, a non-Newtonian phenomenon called shear thinning/shear thickening (the generalized viscosity is decreasing/increasing function of the shear rate).
The constitutive equations of the form (2.9) are also suitable to describe fluids with the activation criteria.Bingham and Herschel-Bulkey fluids [5,62] can be written in the form where also satisfy (G1)-(G4).Activated Euler fluids, see [8], are described by the formula If ν is constant, then the fluid behaves as the Navier-Stokes fluid once |Dv| exceeds the activation parameter δ * .It is straightforward to check that these models fulfil (G1)-(G4).(A large-data analysis of activated Euler fluids, for steady and unsteady flows and for various boundary conditions including complete slip as well as no-slip is developed in [8]. ) The examples discussed above are summed up in the following Table 2, where for brevity we set all physical constants to be 1, except the exponent r related to the (r, r ′ )-coercivity condition (G4).If r does not appear in the equation, then the model leads to (G4) with r = 2.
Examples of two classes of explicit constitutive relations covered by (2.9).The first two lines describe the Navier-Stokes and the power-law fluids with the power-law index r ∈ (1, ∞), r ′ := r/(r−1).In these rows, the formulas in the two columns are equivalent, see (3.2).The formulas in the third and fourth row hold for r ∈ (−∞, ∞); in the range r > 1 the models in the left and right columns behave in the same way for large values of |Dv| and |S|; for r ≥ 1 and r ′ ≥ 1 these formulas satisfy (G2) or (G2*), i.e. the response is monotone.For r < 1 and for r ′ < 1 the formulas in the left and right column behave differently; their response is non-monotone, see [84,77,97,66] for details.The last row describes an activated Euler fluid (left) and a Bingham fluid (right); after the activation takes place, both fluids behave as a Navier-Stokes fluid.
In conclusion, constitutive equation (2.9) covers models designed to describe two non-Newtonian phenomena: shear thinning/shear thickening and the presence of activation criteria in a simple shear flow.(A detailed description of non-Newtonian phenomena is given for example in [86].)Interestingly, (2.9) covers one additional phenomenon, called normal stress differences, which is usually attributed to viscoelastic nature of the fluid; see Perlácová and Průša [97] for more details.
Constitutive equations (boundary conditions) covered by (2.10).Boundary equations are of the same importance as the constitutive equations in the bulk.This assertion can be supported by a recent study [31], where flows are shown to change quantitatively in an essential manner in dependence of the boundary conditions (only linear Navier's slip boundary conditions were tested with their limiting cases (no slip vs complete slip).
Navier proposed a linear constitutive relation (2.5) as the proper boundary condition in [95].Stokes [112] discusses the boundary conditions at length and one variant that he considers concerns a nearly quadratic relation between wall shear stress and the velocity.He states: ". . .when the velocity is not small the tangential force called into action by the sliding of water over the inner surface of the pipe varies nearly as the square of the velocity".Mooney [93] proposed a more general form of slip and introduced a technique that evaluates this relationship.Comprehensive overviews concerning general boundary conditions and slipping mechanisms can be found for example in [58] or [105].More detailed discussions concerning the boundary conditions, their importance, including references to earlier studies are available in [86, A.4] and [85,Sect. 4.6].An overview of basic models of the type (2.10) is given in Table 3.
Examples of two classes of boundary conditions belonging to (2.10).
The first line describes Navier's slip.The second, third and fourth rows describe nonlinear slip of polynomial type with the power-law index q ∈ (1, ∞), q ′ := q/(q − 1).The last row describes the stick-slip boundary condition (right column) and the boundary condition that describes the complete slip before activation and Navier's slip once the activation takes place (left column).
Measurements for molten polymers clearly document that there are nonlinear responses between s and v τ including various activations.For example, Hatzikiriakos in [58] considers models of power-law type, i.e., (3.6) v τ = γ|s| q ′ −2 s, referring to [104] for the value q ′ = 3, to [59] for q ′ = 4 and to [63], [67] for q ′ = 7. Nonlinear responses (3.6) include the material parameter γ that can be a function of other relevant quantities.Besides [58], models of the type (3.6) were studied also in [60], [32], [76], [36].Stick-slip boundary conditions were added to large-data and long-time existence analysis in [17,18].A treatment of complex nonlinear (non-monotone) boundary conditions of stickslip type is presented within the context of analysis of Kolmogorov's two-equation model of turbulence in [19].How the choice of boundary conditions influences the definition of proper function spaces and the subsequent analysis is studied for different boundary conditions within the context of activated Euler's fluids in [8].

ε-approximations of implicit constitutive equations
In our preceding study [20], while studying problems of parabolic type (think of a nonlinear heat equation), we found out the structural assumptions on the implicit function that characterizes the relation between the (heat) flux and the (temperature) gradient which allows us to built a theory parallel ("equivalent" -in the sense specified in [20]) to that of the maximal monotone r-coercive (potentially multi-valued) graphs.Here, we intentionally prefer to avoid using the concept of maximal monotone graphs and we wish to present a theory based only on the assumptions on G and g that require only the knowledge of basic tools of calculus.
The intention of this section is to introduce ε-approximations of the functions G and g and summarize their nice properties: the approximations always lead to L 2 -coercivity of S and Dv (and similarly for s and v τ ), on the ε-approximation level, S is always a function of Dv and upon inserting it into − div S one obtains a Lipschitz continuous uniformly monotone operator.Also here, we follow closely the approach developed in Sect. 4 in [20].Lemma 4.1.Let G satisfy (G1)-(G4) for any r > 1, let g satisfy (g1)-(g4) for any q > 1 and let ε ∈ (0, 1).Then the approximating functions defined by Moreover, there exist two unique functions (single-valued mappings) If, in addition, for any bounded measurable U ⊂ Q and for S ε , D ε : U → R d×d then there exist S ∈ L r ′ (U ) and D ∈ L r (U ) so that (for subsequences) then there exist s ∈ L q ′ (V ) and v ∈ L q (V ) so that (for subsequences) Proof.See [20], Lemma 4.1, Lemma 4.2 and Lemma 4.4.Let Ω ⊂ R d be a domain.We say that Ω is a Lipschitz domain/C 1,1 -domain and we write Ω ∈ C 0,1 /Ω ∈ C 1,1 if, roughly speaking, the boundary ∂Ω can be covered by finite number of overlapping C 0,1 /C 1,1 mappings.For t ∈ (0, T ], we denote Q t := [0, t) × Ω and Γ t := [0, t) × Γ, and we recall that Q := Q T and Γ := Γ T .The abbreviation a.a.stands for almost all, while a.e.stands for almost everywhere.Generic constants, that depend only on the data but are independent of any approximation parameter, are denoted by C and may vary from line to line.
For a Banach space (X, • X ), its dual is denoted by X * .For x ∈ X and x * ∈ X * , the duality is denoted by x * , x X .For r ∈ [1, ∞], we denote (L r (Ω), • r ) and (W 1,r (Ω), • W 1,r (Ω) ) the corresponding Lebesgue and Sobolev spaces with the norms defined in the standard way.Bochner spaces are denoted by L r (0, T ; X).We use the notation L r (Ω; R d ) and L r (Ω; R d×d ) for Lebesgue spaces of vector-or matrix-valued functions, respectively.C ∞ 0 (U ) stands for smooth functions with compact support in an open set U .
Next, we define the function spaces of divergenceless functions with the normal component vanishing on the boundary that are relevant to our setting.For r > 1, we set (5.1) Here, E(Ω) := {u ∈ L 2 (Ω; R d ); div u ∈ L 2 (Ω)}, for which it is known that the trace operator has well defined normal component (in the sense of distributions) and there holds (u • n) ∈ (W 1/2,2 (∂Ω; R d )) * , see e.g.[39].Referring back to (5.1), for any r ∈ [ 2d d+2 , ∞) and z > r, one has , where all the embeddings are continuous and dense.
Finally, we define We are now in a position to precisely formulate our main result.
Theorem 5.1.Let Ω ∈ C 0,1 , T > 0, b ∈ L r ′ (0, T ; V * r ) and v 0 ∈ H be arbitrary.Let G and g be arbitrary functions satisfying (G1) -(G4) with r ∈ ( 2d d+2 , ∞) and (g1)-(g4) with q ∈ (1, ∞).Then there exists a weak solution to (2.3), (2.9) and (2.10) in the following sense: there exist (v, S, s) such that for z := max{r, q, (d+2)r (d+2 the balance of linear momentum is satisfied in a weak sense, i.e. for a.a.t ∈ (0, T ) and for all ϕ ∈ V z (5.5) the constitutive equations (2.9) and (2.Also, for all t ∈ (0, T ) the energy inequality holds, i.e., (5.9) In addition, if Ω ∈ C 1,1 and b ∈ L r ′ (0, T ; V * r ), then there exists a pressure p ∈ L z ′ (0, T ; L z ′ (Ω)) such that As stated in Sect.2, the assumption (g4) eliminates no-slip and perfect-slip boundary conditions from the analysis presented here.However, one can incorporate both conditions into the analysis.The case of s = 0 on Γ is in fact easy as the boundary term just vanishes.For no-slip boundary conditions, one needs to change the function space for the velocity and pay attention to differences associated with the reconstruction of the pressure (see [120] and [8] for details).
In the proof of Theorem 5.1, we use the following powerful convergence result that is a consequence of the properties of the suitably constructed Lipschitz approximations of Bochner functions, see [10].Here, we provide a simplified version (omitting the discussion concerning Lipschitz approximations) suited to the analysis in Sect.6. Lemma 5.2.For any interval I 0 ⊂ (0, T ) and any ball B 0 ⊂ Ω, set Q 0 := I 0 × B 0 .Assume that for δ → 0+ the following convergences hold: which is a weak formulation of , and for every k ∈ N there exists a family {Q δ,k } δ∈(0,1) fulfilling Proof.See [10, Theorem 2.2 and Corollary 2.4], which is adapted to our setting.Compare also with [45,14].
The parameter ε is used to approximate the constitutive equations (2.9) and (2.10) in the same way as presented in Sect.4, see formulas (4.1a) and (4.1b).The motivation for such a choice is that the approximations G ε and g ε possess much better properties in comparison with the properties of G and g, as summarized in Lemma 4.1.In particular, for ε-approximation, we can use (4.5) and thus stay on the level of nonlinear yet uniformly Lipschitz continuous and uniformly monotone operators.Consequently, we come from the (r, r ′ )-coercivity for Dv and S to the L 2 -coercivity for the ε-approximations Dv ε and S ε .The mathematical theory for problems of parabolic type with such nonlinear operators is well known, see [20,Appendix C] for example.Since this type of ε-approximation changes the r-structure to 2-structure, it is suitable to define the following auxiliary numbers µ := min{r, 2}, µ ′ := max{r ′ , 2}, ν := min{r ′ , 2}, ν ′ := max{r, 2}.
We summarize the results proved above (in Subsection 6.2) in the following lemma.
Next, we need show that G(S, Dv) = 0 a.e. in Q and g(s, v τ ) = 0 a.e. on Γ.Here, we have several possibilities.First, if z = r, i.e. if q ≤ r and r ≥ 3d+2 d+2 , one can simply use ϕ := v in (6.37) and therefore, one might mimic the theory developed in [20].Next, if r ≥ 3d+2 d+2 but q > r, one may observe that the choice ϕ := v is admissible in the second and third integral on the left-hand side and also in the term on the right-hand side.Consequently, one might try to generalize the concept of the Gelfand triple and define properly a duality pairing ∂ t v, ϕ and again to mimic the theory from [20].Finally, in the case r < 3d+2 d+2 we cannot set ϕ := v in the second term and therefore the theory from [20] cannot be adapted directly to our case.As the novel method developed here, which is based also on the use of Lemma 5.2, covers also the simple cases discussed above, we just present a unified procedure for all values of r > 2d d+2 .
6.3.2.Identification of nonlinearities on the boundary.By virtue of (6.36j) and Egoroff's theorem, for every η > 0 there exists Γ η so that |Γ \ Γ η | < η and (6.38) It follows from (g4) and Young's inequality that This together with (6.38) implies that {s δ } is a bounded sequence on Γ η .By (6.36h) and Lebesgue's Dominated Convergence Theorem, we conclude that, as Then from Lemma 4.1, g(s, v) = 0 a.e. on Γ η , and letting η → 0 + , we obtain that g(s, v) = 0 a.e. on Γ, and also that for all η, (6.39) 6.3.3.Identification of nonlinearities inside the domain.Identification in Q is not so straightforward, especially due to the lack of proper duality pairing in the convective term and consequently potential failure of the energy equality for the limiting equation.We start with subtracting the weak formulation for v δ (6.31) from the one for v (5.5) and integrate the difference over (0, T ).We deduce that, for all ϕ ∈ L z (0, T ; V z ), (6.40) Next, we localize the above formulation and also omit writing the duality pairing in V r .Indeed, by using the classical theory for r-Stokes problems, we can find 7 F δ and F such that (6.41) F δ → F strongly in L r ′ (Q) 7 We can set F as F := |∇w| r−2 ∇w, where w solves the homogeneous Dirichlet problem − div(|∇w| r−2 ∇w) = −∇π + b, div w = 0 in Ω.
Similarly from b δ we come to F δ .
6.3.4.Energy inequality.Next, we show that (5.9) holds true.For 0 < ǫ ≪ 1 and t ∈ (0, T − ǫ), let η ∈ C 0,1 ([0, T ]) be defined as a piece-wise linear function of three parameters, such that (6.47) We multiply (6.34) by η, and integrate the result over (0, T ) to deduce, after integrating by parts, that The next step is to take the limit as δ → 0 + .We know that (S δ : Dv δ ) ≥ 0 and (s δ • v δ ) ≥ 0 because G(S δ , Dv δ ) = 0 in Q, g(s δ , v δ ) = 0 on Γ and we have (G2 * ) from Lemma 4.2 and an analogous result holds for the function g as well.Therefore, from the above identity we deduce (6.48) For the first term, we can use the weak lower semicontinuity of the norm.For the products S δ : Dv δ and φ δ (|v δ | 2 )s δ • v δ , we use (6.46), (6.39) and the weak convergence of v δ on Γ η .For the duality term, we use (6.1) and (6.36b), and get that Next, we proceed with ǫ, η → 0 + , then , and finally, thanks to v ∈ C w ([0, T ]; H) and the fact that the other terms are well-defined, we obtain the energy inequality (5.9) for any t ∈ (0, T ).
6.4.Attainment of the initial condition.Considering η introduced in (6.47), we multiply (6.31) by ϕη, where ϕ ∈ V z is arbitrary.Integrating the result over (0, T ), we get Next, we apply integration by parts in the first term (using the properties of η and the fact that ϕ is independent of t).Then we take the limit δ → 0 + .Using arguments from Sect.6.3, in paricular the convergence result (6.36c) to take the limit in the convective term, we conclude that Also, taking the limes superior for t → 0 + in the energy inequality (5.9), we obtain that lim sup t→0+ v(t) 2 2 ≤ v 0 2 2 .The last two pieces of information imply the strong convergence in H as claimed in (5.8).6.4.1.Existence of an integrable pressure.In order to reconstruct the pressure, we need to assume that Ω ∈ C 1,1 .The procedure to obtain an integrable pressure for problems with slipping boundary conditions, is explained in [14] in detail.This is why we merely show here a formal estimate concerning the "best" integrability of p.To do so, we assume that there exists an integrable pressure p such that (5.10) holds.Moreover, we may assume that Ω p(t) dx = 0 for a.a.t.Next, we find φ solving the equation ∆φ Using classical theory we know that such a φ exists and satisfies the estimate (6.49) Setting ϕ := ∇φ in (5.10), we have (note that ϕ ∈ V z and also that the term with the time derivative vanishes) (6.50) By use of (6.49), the Hölder inequality and the Trace Theorem we can deduce that After integrating it over (0, T ), thanks to the definition of z ′ and a stronger assumption on b ∈ L z ′ (0, T ; V * z ), all terms on the right-hand side are bounded and we conclude that p ∈ L z ′ (0, T ; L z ′ (Ω)).

Extensions: existence results for related problems (a summary)
The above established result concerns isothermal homogeneous and incompressible fluids; these properties can be seen as limitations and one may wishes to develop a theory in the similar spirit as above for heat-conducting or inhomogeneous or compressible fluids or for fluids that share more of these properties.Also the constitutive equation (2.9) does not cover viscoelastic rate-type or integral models.Similarly, the boundary condition (2.10) does not include dynamic boundary conditions.Below, we provide references that can be relevant to anyone who would like to extend the study in the directions indicated.
When one wishes to include the dependence of the viscosity on the temperature, the system of governing equations has to be completed by the formulations of the balance of energy and the second law of thermodynamics.One also needs to specify the constitutive equation relating the heat flux to the temperature gradient, boundary conditions for the temperature etc.These extensions give rise, within the context of weak solutions, to several concepts of solution.A sound long-time and large-data existence theory for heat-conducting fluids described by incompressible Navier-Stokes-Fourier equations goes back to [50,11].The first existence result for nonlinear models of the power-law type is due to Consiglieri [38] for r ≥ 3d+2 d+2 .Then, a very similar theory was obtained for smaller r's in [23] but merely for explicit models of the type S = 2ν(θ, |Dv| 2 )Dv, where θ is the temperature.The first fully implicit approach with material parameters depending on the temperature was developed in [88], where the authors dealt with a specific activated model with the activation depending on the temperature that also allowed the model to range from activated Euler through the Navier-Stokes regime to a Bingham-type response.
The mathematical theory for unsteady flows of inhomogeneous Navier-Stokes fluids is developed in the book by P.-L.Lions [80].Inhomogeneous isothermal fluids of power-law type are analyzed in [54], while heat-conducting processes for such fluids are treated in [51].Compressible non-Newtonian fluids of power-law type serve represent a completely open field from the point of view of the large-data analysis of relevant initial-and boundary-value problems (see [3] for one of the first attempts and further references).
Regarding viscoelastic rate-type fluids, where for a part of the Cauchy stress one has an additional evolutionary (nonlinear) equation, we distinguish two basic classes: without stress diffusion and with stress diffusion.The first class includes the standard Maxwell, Oldroyd-B or Giesekus models, and the additional equation is of transport type.The long-time and large-data mathematical theory goes back to [81], where Oldroyd-B type models with a corotational time derivative are studied.The type of derivative is however non-physical and simplifies the analysis tremendously.For Giesekus type of models with more general objective derivatives the idea of existence proof is presented in [89] and rigorously proven in the planar case in [16].Regarding the analysis of viescoelastic rate-type fluids with stress diffusion, we refer to [4] where a very robust theory is developed and other relevant studies are cited.
Finally, dynamic boundary conditions (see [58] for their relevance to observations connected with the experiments regarding molten polymers) are the subject of a recent investigation from the point of view of analysis of partial differential equations.The theory for the Stokes system with such dynamic boundary conditions is developed in [1].
Extensions in other directions are possible.Examples include the analysis of rapidly shearthickening fluids (see [57]) or fluids with a priori bounded velocity gradient (see [92,26]).Also, one can consider instead of (2.9) a more general class of incompressible fluids given by the relation G(T , Dv) = 0, which allows one to naturally include naturally fluids with pressure and shear-rate dependent viscosity; here we refer to [21,23,22,25] for further details.

Uniqueness, smoothness, open problems and concluding remarks
We have studied long-time and large-data mathematical properties of unsteady internal flows of incompressible fluids with frictional properties characterized by implicit constitutive equations in the bulk and on the boundary.The developed theory that has origin in the seminal works of O. A. Ladyzhenskaya addresses positively the question concerning existence of weak solutions to large classes of fluids as well as boundary conditions.The structural assumptions characterizing the admissible class of fluids and boundary conditions are expressed in terms of basic tools of calculus and can be checked directly for a given constitutive equation without any deeper knowledge of concepts of the operator theory.Despite broad applicability of the developed theory to many models in different scientific areas, the analysis for fluids satisfying (G1) and (G4) but having non-monotone response (as for example the model computationally tested in [66]) is an open problem.(Note however that non-monotone responses in the boundary conditions can be included, see [19] for details.) The main achievement of this study lies in a novel existence theory.Ladyzhenskaya's interest in these fluids was however motivated by the uniqueness and smoothness of these solutions.Concerning the uniqueness of weak solution, one can easily observe that for models fulfilling (G1)-(G4) we automatically/easily obtain the uniqueness of the velocity field for all r > 1 provided that we neglect the convective term div(v ⊗ v); compare with [20].For the complete model, i.e. for model with the convective term, already Ladyzhenskaya was able to show the uniqueness of a weak solution for the models of the form provided that r ≥ d+2 2 .To date the most general uniqueness result is due to [15], where the authors prove uniqueness of a weak solution in three dimensions (for sufficiently regular data) and a class of models having the r-growth with r ≥ for all (S i , D i ) such that G(S i , D i ) = 0, i = 1, 2, which is fulfilled by the models given in (8.1) above, but it is much more restrictive in comparison with the condition (G2*) (or its equivalent form (G2)) needed in the existence theory.We can slightly strengthen these uniqueness results by considering fluids that behave as the Navier-Stokes fluid prior the activation, i.e. for |Dv| ≤ δ * , where δ * > 0 can be arbitrary, and behave as a power-law fluid with r ≥ 3d+2 d+2 once the activation takes place, i.For this model, one can then establish both the existence and the uniqueness of a weak solution for r sufficiently large (r ≥ 3d+2 d+2 ).We wish to emphasize that the constant δ * can be chosen arbitrarily large.On the other hand, for r < 3d+2 d+2 we have a counterexample to uniqueness thanks to [27] in the class of very weak solutions.Hence, a natural open problem is the (non)uniqueness of a weak solution for smooth data in natural function spaces also for r < 11  5 in dimension three.Finally, concerning smoothness of weak solutions to nonlinear models studied in this paper, there is one striking open problem.Independently of whether one excludes or includes the convective term and independently of the value of the parameter r, it is not clear (even in the situation when we know that there is a unique weak solution) whether, for smooth but large data, there exists a global-in-time C 1,α -solution for any special case of the problem (2.3), (2.9) and (2.10) in three dimensional setting.The regularity theory in two dimensions is available, see e.g.[15,42,68,69,24], the theory in dimension three is however basically untouched.
) d×d and G(0, 0) = 0; (G2) for almost all (S, D) ∈ R d×d × R d×d : ν is a positive constant in the case of Bingham fluids and is a polynomial function of Dv in the case of Herschel-Bulkley fluids.It is proved in [20, Appendix, Example A.3] that Bingham fluids satisfy (G1)-(G4) with r = 2. Following the same line of argument, one can check that Herschel-Bulkley fluids with v(|Dv| 2 ) = (1 + |Dv| 2 )

S
ε : D ε dx dt ≤ U S : D dx dt, then G(S, D) = 0 a.e. in U and S ε : D ε ⇀ S : D weakly in L 1 (U ).Analogously, if for a bounded and measurable V ⊂ Γ and for s