Large data existence theory for three-dimensional unsteady flows of rate-type viscoelastic fluids with stress diffusion

We prove that there exists a weak solution to a system governing an unsteady flow of a viscoelastic fluid in three dimensions, for arbitrarily large time interval and data. The fluid is described by the incompressible Navier-Stokes equations for the velocity $v$, coupled with a diffusive variant of a combination of the Oldroyd-B and the Giesekus models for a tensor $\mathbb{B}$. By a proper choice of the constitutive relations for the Helmholtz free energy (which, however, is non-standard in the current literature, despite the fact that this choice is well motivated from the point of view of physics) and for the energy dissipation, we are able to prove that $\mathbb{B}$ enjoys the same regularity as $v$ in the classical three-dimensional Navier-Stokes equations. This enables us to handle any kind of objective derivative of $\mathbb{B}$, thus obtaining existence results for the class of diffusive Johnson-Segalman models as well. Moreover, using a suitable approximation scheme, we are able to show that $\mathbb{B}$ remains positive definite if the initial datum was a positive definite matrix (in a pointwise sense). We also show how the model we are considering can be derived from basic balance equations and thermodynamical principles in a natural way.


Introduction
We aim to establish a global-in-time and large-data existence theory, within the context of weak solutions, to a class of homogeneous incompressible rate-type viscoelastic fluids flowing in a closed three-dimensional container. The studied class of models can be seen as the Navier-Stokes system coupled with a viscoelastic ratetype fluid model that shares the properties of both Oldroyd-B and Giesekus models and is completed with a diffusion term. Such models are frequently encountered in the theory of non-Newtonian fluid mechanics, see [13,12] and further references cited in [12].
In order to precisely formulate the problems investigated in this study, we start by introducing the necessary notation. For a bounded domain Ω ⊂ R 3 with the Lipschitz boundary ∂Ω and a time interval of the length T > 0, we set Q := (0, T ) × Ω for a time-space cylinder and Σ := (0, T ) × ∂Ω for a part of its boundary. The symbol n denotes the outward unit normal vector on ∂Ω and, for any vector z, the vector z τ denotes the projection of the vector to a tangent plane on ∂Ω, i.e., z τ := z − (z · n)n. Then, for a given density of the external body forces f : Q → R 3 , a given initial velocity v 0 : Ω → R 3 and a given initial extra stress tensor B 0 : Ω → R 3×3 >0 (here R 3×3 >0 denotes the set of symmetric positive definite (3 × 3)-matrices), we look for a vector field v : Q → R 3 , a scalar field p : Q → R and a positive definite matrix field B : Q → R 3×3 >0 solving the following system in Q: div v = 0, (1.1) The parameters γ ∈ (0, 1), ν, λ, σ > 0, δ 1 , δ 2 ≥ 0 and a ∈ R are given numbers.
The main result of this study can be stated as: Let v 0 and B 0 be such that the initial total energy is bounded. Then, for sufficiently regular f , there exists a global-in-time weak solution to (1.1)-(1.6).
Although the above result is stated vaguely, we would like to emphasize that we are going to establish the long-time existence of a weak solution for large data and for three-dimensional flows. A more precise and rigorous version of the above result including the correct function spaces and the properly defined weak formulation is stated in the Theorem below.
We complete the introductory part by providing the physical background relevant to the studied problem and by recalling earlier results relevant to the problem (1.1)-(1.6) analyzed here. In (1.7), • v stands for the material time derivative of v, i.e., • v = ∂ t v + (v · ∇)v. Defining similarly the material time derivative of a tensor B as which is supposed to hold true in Q and which is completed by the initial conditions (1.5), (1.6) fulfilled in Ω and by the boundary conditions (1.4) on Σ that take the form: v · n = 0, (1.12) We provide several comments regarding (1.8)-(1.11) as well as the boundary conditions (1.12)-(1.14). The Navier slip boundary condition (1.13) is considered here just for simplicity; note that for smooth domains, namely if Ω ∈ C 1,1 , we can introduce the pressure p as an integrable function, e.g., by using an additional layer of approximation as in [4], see also [9,8] or [2] which discuss the treatment of the pressure in evolutionary models subject to the Navier boundary condition. Nevertheless, since we shall always deal with formulation without the pressure (see the Definition), we can also treat the Dirichlet boundary condition, as well as very general implicitly specified boundary conditions see, e.g., [24,5,6] or [2]. The Neumann boundary condition for B is considered here only for simplicity and without any specific physical meaning.
A further aspect, which makes the above system more complicated than the Navier-Stokes equation is the form of the Cauchy stress tensor T as in (1.8). The term −pI + 2νDv corresponds to the standard Newtonian fluid flow model with a constant kinematic viscosity ν. The next part of the Cauchy stress, which depends linearly on B, appears in all the viscoelastic rate-type fluid models -see, e.g., [20, (7.20b), (8.20e)], [15, (6.43e)] or [12, (43a)]. On the other hand, the addition of the term 2aµγ(B 2 − B) is, to our best knowledge, considered here for the first time. The fact that we require that γ is positive (and strictly less than 1) plays a key role in the analysis of the problem, as will be shown below. Note that the linearization of T with respect to B when B is close to the identity I yields and we recover the standard form of T (after possible redefinition of the pressure).
The quantity B takes into account the elastic responses of the fluid and the equation (1.11) describes its evolution in the current configuration (Eulerian coordinates), just as the velocity v. It is frequent to call the tensor µ(B − I) the extra stress or conformation tensor and to denote it by τ . More importantly, since the material derivative of B is not objective, it must be "corrected" and this is the reason, why in (1.11) the derivative ⋄ B appears. The parameter a in the definition of ⋄ B determines the type of the objective derivative. The case a = 1 leads to the upper convected Oldroyd derivative, that has favourable physical properties and that leads to a clear interpretation of B within the thermodynamical framework developed in [26], see also [27,22,23,21]. Next, the case a = 0 leads to the corrotational Jaumann-Zaremba derivative and this is the only case for which the analysis is much simpler than in other cases. Furthermore, if a ∈ [−1, 1], one obtains the entire class of Gordon-Schowalter derivatives. However, it turns out that the physical properties of these derivatives are irrelevant for the analysis presented below (except the case a = 0), therefore we may take any a ∈ R. For a = 1 and λ = 0 we distinguish two cases: if δ 1 > 0 and δ 2 = 0 we obtain the classical Oldroyd-B model while if δ 1 = 0 and δ 2 > 0 we get the Giesekus model. Next, by considering a ∈ [−1, 1], we obtain the class of Johnson-Segalman models. If we further let λ > 0, we are introducing diffusive variants of the previous models. It has been observed that including the diffusion term in (1.11) is physically reasonable, see, e.g., [13] or [12] and references therein. However, up to now, it has been unknown what precise form should the diffusion term take and also whether it actually helps in the analysis of the model. Our main result provides a partial answer to this question, namely: for γ ∈ (0, 1) and with the diffusion term being of the form ∆B (or more generally, a linear second order operator), the global existence of a weak solution is available.
The reader familiar with the equations describing flows of the standard Oldroyd-B viscoleastic rate-type fluid can identify two deviations in the set of equations (1.9)-(1.11) studied hereafter. We provide a few comments on these differences.
The first deviation concerns the incorporation of the stress diffusion term, i.e. the term −∆B, into the equations. Following the pioneering work of [13] it is clear that a quantity related to |∇B| 2 has to be added into the list of underlying dissipation mechanisms. On the other hand, the precise form in which stress diffusion should appear depends on the choice of a thermodynamical approach and specific assumptions. In fact, using the thermodynamical concepts as in [20] or [12], one can derive models, where the stress diffusion term takes the form −B∆B − ∆BB, −B 1 2 ∆BB 1 2 etc., however, we would prefer −∆B simply because it coincides with the form proposed by [13], and, from the perspective of PDE analysis and numerical approximation, one prefers to deal with stress diffusion that leads to a linear operator.
The second deviation from usual viscoelastic models consists in the presence of the term (B 2 − B) in the Cauchy stress tensor, see (1.8). This term arises if we slightly modify energy storage mechanism and apply the thermodynamic approach as developed in [20]. In what follows, we shall give a clear interpretation and a thermodynamic derivation of our model.

1.2.
Thermodynamical derivation of the model. Viscoelastic models with (nonlinear) stress diffusion, but without the term B 2 in the stress tensor are derived, e.g., in [20] and [12] even in the temperature-dependent case. Here, we will briefly explain the approach in a simplified isothermal setting (sufficient for the purpose of this study), referring to the cited works for the derivation in a complete thermal setting and for more details.
First, we postulate the constitutive equation for the Helmholtz free energy in the form , where µ > 0 and γ ∈ [0, 1] is a parameter interpolating between two forms of the energy. The choice γ = 0 would lead to a standard Oldroyd-B diffusive model. To our best knowledge, the case γ > 0 was not considered before in literature. The term 1 2 γ|B−I| 2 , which is newly included in ψ is obviously convex with the minimum at B = I and depends only on tr B and on tr(BB), i.e., on invariants of B, therefore it does not violate any of the basic principles of continuum physics. Moreover, such an addition does not affect the first three terms in the asymptotic expansion of ψ near I, on the logarithmic scale. To see this, let H denote the Hencky logarithmic tensor satisfying e H = B (which exists due to the positive definiteness of B). Using Jacobi's identity, we compute that ). On the other hand, we easily get ) and we see that for B being close to identity, the form of ψ is almost independent of the choice of parameter γ and the second part of ψ in (1.15) can be just understood as a correction for large values of B.
Next, we show how the constitutive equation for T (see (1.8)) appears naturally if we start with the choice of the Helmholtz free energy (1.15) and require that the form of the equation for B is given by (1.11). For the derivation, we followed the approach developed in [20] that stems from the balance equations and requires the knowledge of how the material stores the energy, but we simplify the derivation presented there by assuming that the density is constant (in fact we set for simplicity ̺ = 1 and hence div v = 0) and the flow is isothermal, i.e., the temperature θ is constant as well. Under these assumptions the balance equations of continuum physics (for linear and angular momenta, energy and for formulation of the second law of thermodynamics) take the form where e is the (specific) internal energy, η is the entropy, ξ is the rate of entropy production, T is the Cauchy stress tensor and the quantities j e , j η represent the internal and the entropy fluxes, respectively. Since the quantities ψ, e, θ and η are related through the thermodynamical identity e = ψ + θη, we can easily deduce from above identities that To evaluate the last term, we rewrite (1.11) as Next, it follows from (1.15) that where J is defined as Consequently, taking the inner product of (1.17) with J we observe that (since BJ = JB, the term with Wv vanishes) (1.18) To evaluate the terms on the last line, we use the symmetry and the positive definiteness of the matrix B to obtain Similarly, we obtain Thus, using (1.18)-(1.20) in (1.16), we conclude that (1.21) Hence, assuming that the fluxes fulfil and setting (compare with (1.8)) which gives the nonnegative rate of the entropy production. Moreover, we have seen how the form of the Cauchy stress tensor T in (1.8) is dictated by the second line in (1.21). Furthermore, we can also see in (1.23) (and also in the last line of (1.19)) how the choice of the free energy (1.15) affects the entropy production due to the presence of the diffusive term ∆B in (1.3).
1.3. The concept of weak solution and energy (in)equality. In order to introduce the proper concept of weak solution, we first derive the basic energy estimates based on the observations from the previous section. First, taking the scalar product of (1.10) and v, we deduce the kinetic energy identity and replacing the term T · Dv from the equation (1.16), and using then also (1.22) and (1.23), we finally obtain Integrating the above identity over Ω, using integration by parts and the boundary conditions (1.12)-(1.14), we obtain Definition. Let T > 0 and assume that Ω ⊂ R 3 is a Lipschitz domain. Let The initial conditions are satisfied in the following sense Moreover, we say that the solution satisfies the energy inequality if, for all t ∈ (0, T ): (1.28) In the above definition we used the following notation. By L p (Ω) and W n,p (Ω), 1 ≤ p ≤ ∞, n ∈ N, we denote the usual Lebesgue and Sobolev space, with their usual norms denoted as · p and · n,p , respectively. The trace operator that maps W 1,p (Ω) into L q (∂Ω), for certain q ≥ 1, will be denoted by T . Further, we set W −1,p ′ (Ω) = (W 1,p (Ω)) * , where p ′ = p/(p − 1). We shall use the same notation for the function spaces of scalar-, vector-, or tensor-valued functions, but we will distinguish the functions themselves using different fonts such as a for scalars, a for vectors and A for tensors. Also, we do not specify the meaning of the duality pairing ·, · , assuming that it is clear from the context. Moreover, for certain subspaces of vector valued functions, we shall use the following notation: Occasionally, we shall denote the standard inner products in L 2 (Ω) and L 2 (∂Ω) as (·, ·) and (·, ·) ∂Ω , respectively. The Bochner spaces of mappings from (0, T ) to a Banach space X will be denoted as L p (0, T ; X) with the norm · L p (0,T ;X) = denotes a space of weakly continuous functions, i.e., for every f ∈ C weak (0, T ; X) and every g ∈ X * there holds lim The symbol R 3×3 sym denotes the set of symmetric 3 × 3 real matrices. Furthermore, by R 3×3 >0 we denote the subset of R 3×3 sym which consists of positive definite matrices, i.e., those which satisfy
Let us briefly explain the main difficulties connected with the analysis of the system (1.9)-(1.13) and our ideas how to solve them. In the standard models where γ = 0, to get an a priori estimate for B, the appropriate test function to take in (1.11) is I − B −1 . Then, using (1.9) and (1.10) tested by v, one can eliminate the problematic terms, such as B · Dv coming from the objective derivative. However, the non-negative quantity to be controlled, which comes from the diffusion term, turns out to be just |B − 1 2 ∇BB − 1 2 | 2 and this provides little to no information. In particular, the terms ∇vB appearing in (1.11) are going to be just integrable and it is unclear if one can show strong convergence of B. Instead, one would like to test also by B to achieve control over |∇B| 2 . But this is not possible, since the resulting term ∇vB · B cannot be estimated without some serious simplifications (such as boundedness of ∇v, two or one dimensional setting or small data). Quite remarkably, this problem is solved simply by adding 1 2 γ|B − I| 2 into the constitutive form for ψ. More precisely, considering γ ∈ (0, 1), we observe that the appropriate test function in (1.11) is in fact (1 − γ)(I − B −1 ) + γ(B − I). Indeed, the terms from the objective derivative cancel again due to the presence of γ(B 2 − B) in T. But now, we also get γ|∇B| 2 under control, which is much better information than in the case γ = 0 and it will imply compactness of all the terms appearing in (1.10) and (1.11). We have seen above that such a modification of ψ, and consequently of T, is not ad-hoc and that it rests on solid physical grounds.
The second and also the last major difficulty which we will encounter is how one can justify testing of (1.11) by B −1 on the approximate (discrete level), where B −1 might not even exist. This we overcome by designing a delicate approximation scheme, which takes into account the smallest eigenvalue of B, and also by noting that testing (1.11) only by B yields sufficiently strong a priori estimates for the initial limit passage (in the Galerkin approximation of B).
Up to now, there have been no results on global existence of weak solutions to Oldroyd-B models in three dimensions, including either the standard, or diffusive variants. The closest result so far is probably [25,Theorem 4.1], however there it is assumed that δ 2 > 0 and λ = 0 (Giesekus model), whereas we treat also the case δ 2 = 0, but with λ > 0 (diffusive Oldroyd-B or Giesekus model). Moreover, in [25], only the weak sequential stability of a hypothetical approximation is proved. We, on the other hand, provide the complete existence proof, including the construction of approximate solutions (which, in viscoelasticity, is generally a non-trivial task). In the article [17], Lions and Masmoudi prove the global existence in three dimensions, but only for a = 0 (corrotational case), which is known to be much easier. The local in time existence of regular solutions for the non-diffusive variants of the models above (λ = 0) is proved in the pioneering work [14,Theorem 2.4.]. There, also the global existence for small data is shown. In two dimensions, the problem is solved in [11] in the case λ > 0, δ 1 > 0, δ 2 = 0 (diffusive Oldroyd-B model). There are also global large data existence results in three dimensions for slightly different classes of diffusive rate-type viscoelastic models, but under some simplifying assumptions. For example, in [7] and [3], the authors consider the case where B = bI. This assumption, however, turns (1.11) into a much simpler scalar equation. Moreover, note that if B = bI, then the equations (1.10) and (1.11) decouple (which is not the case in [7] and [3] since there the considered constitutive relation for T is more complicated than here). Furthermore, in [18], the authors consider yet another class of Peterlin viscoelastic models with stress diffusion and prove existence of a global two-or three-dimensional solution. However, the free energy associated with these models depends only on the trace of the extra stress tensor. This is a significant simplification, which can even be seen as unphysical. See also [10] for various modifications of Oldroyd-B viscoelastic models, for which an existence theory is available. Finally, in [1] (see also [16]), the global existence of a weak solution is shown for a certain regularized Oldroyd-B model (including a cut-off or nonlinear p−Laplace operator in the diffusive term in B). Thus, one might argue that since the case γ > 0 could be also seen as a regularization of the original model, we are just proving an existence of a solution to another regularization. However, this argument is not, in our opinion, correct for several reasons. First of all, the "regularization" γ > 0 does not touch the equation (1.11) at all. Second, it is not obvious why the nonlinear term γ(B − I) 2 should have any regularization effect. And, perhaps most importantly, we already showed in Section 1.2 that the model with γ > 0 is physically well founded and worthy of studying in its own right.
Since the topic is quite new and unexplored, we decided, for brevity and clarity of presentation to consider only the isothermal case. However, we believe that the framework and ideas presented here are robust enough to provide an existence analysis also for the full thermodynamical model if the evolution of the internal energy is described correctly. This is the subject of our forthcoming study.
Remark. Finally, we close this section with several concluding remarks on possible extensions, but we do not provide their proofs in this paper.
(i) The Theorem holds also in arbitrary dimensions d > 3 (in d ≤ 2, it is known), however with worse function spaces for the time derivatives and better for the test functions. Indeed, the only dimension-specific argument in the proof below is in the derivation of interpolation inequalities, which are then used to estimate ∂ t v and ∂ t B. Moreover, all of the non-linear terms in (1.25), (1.26) are integrable for arbitrary d if the test functions are smooth. In addition, if d = 2, then we can prove the existence of a weak solution satisfying even the energy equality, i.e., (1.28) holds with the equality sign. (ii) When Ω has C 1,1 boundary, then, in addition, there exists a pressure p ∈ L 5 3 (Q), which appears in (1.2). Then, the test functions in (1.25) need not be divergence-free if we include the term Ω p div ϕ in (1.25). This follows in a standard way, using the Helmholtz decomposition of v (see, e.g., [2] for details). (iii) It is possible to replace (1.12), (1.13) by the no-slip boundary condition v = 0 on ∂Ω. Then, we only need to change the space W 1,2 n to W 1,2 0 , and so on. However, then it seems that the pressure p can be only obtained as a distribution (see [2]).

Proof of the Theorem
Throughout the proof, we shall simplify notation by assuming λ = µ = ν = σ = 1 and refer to Section 1.2 for a detailed computation for general parameters. To shorten all formulae, we also denote The general scheme of the proof is the following: In order to invert the matrix B and to avoid problems with low integrability in the objective derivative, we introduce the special cut-off function where Λ(A) denotes a minimal eigenvalue of A (whose spectrum is real due to its symmetry) 1 . Since eigenvalues of a matrix depend continuously on its entries, the function ρ ε is continuous. Moreover, for any positive definite matrix A there holds ρ ε (A) → 1 as ε → 0 + . We construct a solution by an approximation scheme with parameters k, l and ε, where k, l ∈ N correspond to the Galerkin approximation for v and B, respectively, and ε corresponds to the presence of the cut-off function ρ ε in certain terms. The first limit we take is l → ∞, which corresponds to the limit in the equation for B. This way, the limiting object B is infinite-dimensional and, using the properties of ρ ε , we prove that B −1 exists. With the help of this information, we derive the energy estimates that are uniform with respect to all the parameters. Next, we let ε → 0 + in order to remove the truncation function and finally we take k → ∞, which corresponds to the limiting procedure in the equation for the velocity v.
Then for fixed k, l ∈ N and ε ∈ (0, 1), we look for the functions v k,l ε , B k,l ε of the form where c k,l,ε i , d k,l,ε j , i = 1, . . . , k, j = 1, . . . , l, are unknown functions of time, and we require that v k,l ε , B k,l ε (and consequently the functions c k,l,ε i (t) and d k,l,ε j (t)) satisfy the following system of (k + l) ordinary differential equations in time interval (0, T ): for j = 1, . . . , l.
Here, B ε 0 is defined as Since B 0 (x) ∈ R 3×3 >0 for almost every x ∈ Ω, we have that Λ(B 0 (x)) > 0 for almost all x ∈ Ω. Consequently, using the fact B 0 ∈ L 2 (Ω), we obtain, as ε → 0 + , that Note also that the initial conditions blow up as t → T * − . But using the estimate presented below (see (3.8) valid for all t ∈ (0, T * )), this will be seen never to happen. Thus, we can set T * = T .

3.2.
Limit l → ∞. In this part, we simplify the notation and denote the approximating solution, constructed in the previous section, by (v l , B l ) := (v k,l ε , B k,l ε ). We start by proving estimates independent of l. Since B l (t) and v l (t) belong for almost all t to the linear hull of {W j } l j=1 and {w i } k i=1 , respectively, we can use v l instead of w i in (3.1) and B l instead of W j in (3.2) to deduce, where we used the integration by parts formula and the facts that div v l = 0 and T v · n = 0. Next, it follows from the definition of ρ ε , R and S that Here, the notation C(ε) emphasizes that the constant C depends on ε; we keep this notation in what follows. Summing (3.4) and (3.5) and using the estimate (3.6) to bound the term on the right-hand side, we obtain with the help of Hölder's, Young's and Korn's inequalities that After integrating over (0, T ) with respect to time, we obtain the following bound: where the last inequality follows from the continuity of the projections P k and Q l and from the assumptions on data, namely that v 0 Next, we focus on the estimate for time derivatives. First, it follows from L 2orthonormality of the bases and the estimate (3.7) that Then, since v l is a linear combination of and we can deduce from (3.1) that Finally, it follows from (3.2) and (3.7) that (3.10) Using (3.7), (3.9)-(3.11) and Banach-Alaoglu's theorem, we can find subsequences (which we do not relabel) and corresponding weak limits (denoted with the subscript ε), such that, for l → ∞, we get v l ⇀ v ε weakly in L 2 (0, T ; W 1,2 n,div ), weakly * in L ∞ (0, T ; W 1,∞ (Ω)), (3.14) T v l ⇀ T v ε weakly in L 2 (0, T ; L 2 (∂Ω)), (3.15) B l ⇀ B ε weakly in L 2 (0, T ; W 1,2 (Ω)), (3.16) weakly in L 2 (0, T ; W −1,2 (Ω)). (3.17) Moreover, it follows from (3.12), (3.14), (3.16), (3.17) and from the Aubin-Lions lemma that for some further subsequences, we have (3.13) strongly in L 2 (Q) and a.e. in Q, Using the convergence results (3.12)- (3.20), it is rather standard to let l → ∞ in (3.1)-(3.2). This way, for almost all t ∈ (0, T ), we obtain for i = 1, . . . , k, and for all A ∈ W 1,2 (Ω). Moreover, from (3.16) and (3.17), we get B ε ∈ C(0, T ; L 2 (Ω)) and it is standard to show that B ε (0, ·) = B ε 0 and v ε (0, ·) = P k v 0 .
3.3. Limit ε → 0. In this part we consider the solutions (v ε , B ε ) constructed in the preceding section for ε ∈ (0, 1) and we study their behaviour as ε → 0 + . To do so, we first have to derive estimates that are uniform with respect to ε. Following the ideas used before in the derivation of the model, we wish to test (3.22) by the function . This test function, however, contains B −1 ε and we need to justify that it exists (for any ε ∈ (0, 1)).

3.3.1.
Estimates for the inverse matrix -still ε-dependent. First, we prove that Λ(B ε ) ≥ ε. For this purpose, let z ∈ R 3 be arbitrary and consider 3 (3.22). Due to the properties of B ε (see (3.16)), we know that A belongs to L 2 (0, T ; W 1,2 (Ω)) and we can use it as a test function in (3.22). Upon inserting A into (3.22), we integrate the result over (0, τ ) with some fixed τ ∈ (0, T ). We evaluate all terms in (3.22) separately. For the time derivative, we have where, for the last equality, the definition of B ε 0 was used. Furthermore, we obtain integrating by parts and using the fact that div v ε = 0 and T v ε = 0. Since we also observe, that Hence, we get Consequently, inserting A of the form (3.24) into (3.22), we see that the right-hand side is identically zero. Therefore, relations (3.25), (3.26), (3.27) and (3.28) yield which implies (3.29) B ε z · z ≥ ε|z| 2 for every z ∈ R 3 and a.e. in Q.
Thus, we have the following estimate for the minimal eigenvalue of B ε : Therefore, the inverse matrix B −1 ε is well defined and satisfies Furthermore, since we conclude from (3.7) and (3.30), that Hence, the inverse of B ε exists and B −1 ε ∈ L 2 (0, T ; W 1,2 (Ω)).

3.3.2.
Estimates independent of (ε, k). At this point, we can test (3.22) with J ε defined in (3.23). This way, we obtain Next, we evaluate all terms. Here, we follow very closely the procedure developed in Section 1.2, see the derivation of (1.18) and consequent identities. Since where ψ is defined in (1.15), it is clear that Next, recalling (1.19), we get and due to the fact that B ε J ε = J ε B ε we also have where we used the fact that the trace of Dv ε is identically zero. Hence, using A := J ε (defined in (3.23)) in (3.22) and taking into account the above identities, we deduce that Similarly as in previous section, replacing w i in (3.21) by v ε , we get Thus, summing (3.31) and (3.32) and integrating the result with respect to time t ∈ (0, τ ), we deduce the identity where, for the last inequality we used the continuity of P k , the definition of B ε 0 and the fact that ψ(I) = 0.
Hence Λ(B k ) ≥ 0 and det B k ≥ 0 a.e. in Q. Therefore, using (3.39) again and the continuity of ψ, there exists (still possibly infinite) limit However, since ψ ≥ 0, Fatou's lemma implies that, for almost every t ∈ (0, T ), we have Thus, we deduce that If there existed a set E ⊂ Q of a positive measure, where Λ(B k ) = 0, then also − ln det B k = ∞ on that set, which contradicts (3.40). Thus, we have Then it directly follows from the continuity of Λ, that ρ ε (B ε ) → 1 a.e. in Q. Then, since ρ ε (B ε ) ≤ 1, we further get, by Vitali's theorem, that Using the established convergence results, it is easy to let ε → 0 + in (3.21) and (3.22) and obtain, for almost all t ∈ (0, T ), that Furthermore, we can take the limit in the estimates (3.33), (3.35), (3.36) and (3.37) using either the weak lower semi-continuity of norms or, in the terms which depend on B ε , e.g. Q ρ ε (B ε )|B ε | 2 , we apply (3.41) to conclude the pointwise limit and then use Fatou's lemma. Thus, inequalities (3.33), (3.35), (3.36) and (3.37) continue to hold in the same form, but for (v k , B k ) instead of (v ε , B ε ) and with 1 instead of ρ ε (B ε ). In particular, for almost all t ∈ (0, T ), we have The attainment of initial conditions is standard (see the last section for details in a more complicated case).

3.4.
Limit k → ∞. Since we start from the same a priori estimates as in the previous section, we follow, step by step, the procedure developed when taking the limit ε → 0 + . The only difference is that the term ρ ε (B ε ) is not present. Thus, using the density of {w i } ∞ i=1 in W 3,2 n,div , we obtain, after letting k → ∞, for almost all t ∈ (0, T ), that Moreover, from the weak lower semi-continuity of norms, we obtain the energy inequality (1.28) for almost all t ∈ (0, T ). Furthermore, the same argument as above implies that B is positive definite a.e. in Q. Now observe that, by Hölder's inequality and (3.35), all the terms in (3.42) except the first one, are integrable for every ϕ ∈ L 4 (0, T ; W 1,2 n,div ) ֒→ L 4 (0, T ; L 6 (Ω)). Indeed, for example for the non-linear terms, we get Hence, the functional ∂ t v can be uniquely extended to ∂ t v ∈ L 4 3 (0, T ; W −1,2 n,div ) and we can use the density argument to conclude (1.25). Analogously, we obtain (1.26). Hence, it remains to show that (1.28) holds for all t ∈ (0, T ) and that the initial data fulfil (1.27).
Next, we notice that the function ψ is convex on the convex set R 3×3 >0 . Indeed, evaluating the second Fréchet derivative of ψ, we get >0 , which is obviously a positive definite operator for any γ ∈ [0, 1] and consequently, ψ must be convex on R 3×3 >0 . Further, we integrate (1.28) over (t 1 , t 1 + δ), where t 1 ∈ (0, T ), and divide the result by δ. Using also an elementary inequality Finally, we let δ → 0 + . The limit on the right hand side is standard and consequently, if we show that 1 2 v(t 1 ) then (1.28) will hold for all t ∈ (0, T ). To show it, we notice that due to (3.43) (3.45) v(t) ⇀ v(t 1 ) weakly in L 2 (Ω) as t → t 1 , B(t) ⇀ B(t 1 ) weakly in L 2 (Ω) as t → t 1 , Consequently, due to the weak lower semicontinuity and the convexity of ψ we also have for all t ∈ (0, T ) Ω |v(t)| 2 + ψ(B(t)) ≤ C.
3.4.2. Attainment of initial conditions. First, it is standard to show from the construction and from the weak continuity (3.45), that for arbitrary ϕ, A ∈ L 2 (Ω) there holds Next, using the convexity of ψ and (3.46) (and consequently weak lower semicontinuity of the corresponding integral) and letting t → 0 + in (1.28), we deduce that  which finishes the proof of (1.27) and consequently also the proof of the Theorem.