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

1.1 Mathematical and physical background

The system (1.1)–(1.4) can be rewritten into a more concise form once one recognizes some physical quantities. First of all, let

denote the symmetric and antisymmetric parts of the velocity gradient ∇v, respectively. Then, looking at the equation (1.2), we see that (1.2) is obtained from a general form of the balance of linear momentum, namely

once we set the density ϱ = 1 and require that the Cauchy stress tensor 𝕋 has the form

In (1.7), v̇ stands for the material time derivative of v, i.e., v̇ = ∂_tv + (v ⋅ ∇)v. Defining similarly the material time derivative of a tensor 𝔹 as

we can recognize the presence of a general objective derivative in (1.3). Namely, defining

we can rewrite the system (1.1)–(1.3) into a more familiar form 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:

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) (and in general all boundary conditions allowing the fluid to slip ever so slightly) has recently attracted lot of attention. It was well documented that in certain situations the Navier slip boundary conditions are more appropriate than no slip boundary conditions, we refer e.g. to [12, 20, 23, 25] or [36] and references therein. In addition, it was shown that the Navier slip boundary condition can be understood as an asymptotic limit of no slip boundary conditions in case we consider rough and highly oscillating boundary, see e.g. [1, 6, 9]. Furthermore, for the classical Navier-Stokes equation or the Stokes equation, we can say that the available mathematical theory for no-slip boundary condition has been already “re-proven” for Navier boundary conditions, see e.g. [3] for the existence analysis, [2, 4, 30] for regularity theory for the Stokes system and [7] for a conditional regularity result for Navier-Stokes system. The key difference and also the main mathematical advantage of the Navier slip boundary conditions is, that for smooth domains, namely if Ω ∈ 𝓒^1,1, we can introduce the pressure p as an integrable function, e.g., by using an additional layer of approximation as in [11], see also [15, 16] or [8] 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. [12, 13, 36] or [8]. The Neumann boundary condition for 𝔹 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 𝕋 as in (1.8). The term −p𝕀+2ν𝔻v corresponds to the standard Newtonian fluid flow model with a constant kinematic viscosity ν. The next part of the Cauchy stress, which depends linearly on 𝔹, appears in all the viscoelastic rate-type fluid models - see, e.g., [32, (7.20b), (8.20e)], [24, (6.43e)] or [19, (43a)]. On the other hand, the addition of the term 2aμ β(𝔹² − 𝔹) 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 𝕋 with respect to 𝔹 when 𝔹 is close to the identity 𝕀 yields

and we recover the standard form of 𝕋 (after possible redefinition of the pressure).

The quantity 𝔹 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 μ(𝔹-𝕀) the extra stress or conformation tensor and to denote it by τ. More importantly, since the material derivative of 𝔹 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 𝔹 within the thermodynamical framework developed in [38], see also [33, 34, 35, 39]. 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 ∈ ℝ. For a = 1 and λ = 0 we distinguish two cases: if δ₁ > 0 and δ₂ = 0 we obtain the classical Oldroyd-B model while if δ₁ = 0 and δ₂ > 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., [21] or [19] 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 Δ 𝔹 (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 −Δ𝔹, into the equations. Following the pioneering work of [21] it is clear that a quantity related to ∣∇𝔹∣² 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 [32] or [19], one can derive models, where the stress diffusion term takes the form −BΔB−ΔBB,−B12ΔBB12 etc., however, we would prefer −Δ𝔹 simply because it coincides with the form proposed by [21], 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 (𝔹²-𝔹) 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 [32]. 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 𝔹² in the stress tensor are derived, e.g., in [32] and [19] 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 12 β∣𝔹-𝕀∣², which is newly included in ψ is obviously convex with the minimum at 𝔹 = 𝕀 and depends only on tr 𝔹 and on tr (𝔹𝔹), i.e., on invariants of 𝔹, 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 𝕀, on the logarithmic scale. To see this, let ℍ denote the Hencky logarithmic tensor satisfying e^ℍ = 𝔹 (which exists due to the positive definiteness of 𝔹). Using Jacobi’s identity, we compute that

On the other hand, we easily get

hence we also have

and we see that for 𝔹 being close to identity, the form of ψ is almost independent of the choice of parameter β and the second part of ψ in (1.5) can be just understood as a correction for large values of 𝔹.

Next, we show how the constitutive equation for 𝕋 (see (1.8)) appears naturally if we start with the choice of the Helmholtz free energy (1.5) and require that the form of the equation for 𝔹 is given by (1.11). For the derivation, we followed the approach developed in [32] 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 divv = 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, 𝕋 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

we can easily deduce from above identities that

To evaluate the last term, we rewrite (1.11) as

Next, it follows from (1.5) that

where 𝕁 is defined as

Consequently, taking the inner product of (1.17) with 𝕁 we observe that (since 𝔹𝕁 = 𝕁 𝔹, the term with 𝕎v vanishes)

To evaluate the terms on the last line, we use the symmetry and the positive definiteness of the matrix 𝔹 to obtain

Similarly, we obtain

Thus, using (1.19)–(1.21) in (1.16), we conclude that

Hence, assuming that the fluxes fulfil

and setting (compare with (1.8))

the identity (1.22) reduces to (noticing that −p𝕀 ⋅ 𝔻v = −p div v = 0)

which gives the nonnegative rate of the entropy production. Moreover, we have seen how the form of the Cauchy stress tensor 𝕋 in (1.8) is dictated by the second line in (1.22). Furthermore, we can also see in (1.24) (and also in the last line of (1.20)) how the choice of the free energy (1.5) affects the entropy production due to the presence of the diffusive term Δ 𝔹 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 𝕋 ⋅ 𝔻v from the equation (1.16), and using then also (1.23) and (1.24), we finally obtain

Integrating the above identity over Ω, using integration by parts and the boundary conditions (1.12)–(1.14), we obtain

The identity (1.27) indicates the proper choice of the function spaces for the solution (v, 𝔹) and the form of the (weak) formulation of the solution to (1.1)–(1.6).

1.4 Notation

In order to formulate the definition of a weak solution conveniently, let us fix some notation. By L^p(Ω) and W^n,p(Ω), 1 ≤ p ≤ ∞, 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 𝓣. 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 𝔸 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²(Ω) and L²(∂Ω) 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 ∥⋅∥Lp(0,T;X)=(∫0T∥⋅∥Xp)1p. If X = L^q(Ω), or X = W^k,q(Ω), we will write just ∥⋅∥_L^pL^q, or ∥⋅∥_L^pW^k,q, respectively. The space 𝓒_weak(0, T; X) ⊂ L^∞(0, T; X) denotes a space of weakly continuous functions, i.e., for every f ∈ 𝓒_weak(0, T; X), every t₀ ∈ [0, T] and every g ∈ X^* there holds

The symbol Rsym3×3 denotes the set of symmetric 3 × 3 real matrices. Furthermore, by R>03×3 we denote the subset of Rsym3×3 which consists of positive definite matrices, i.e., those which satisfy

3.1 Galerkin approximation

Following e.g., [31, Appendix A.4], we know that there exists a basis {wi}i=1∞ of Wn,div3,2, which is orthonormal in L²(Ω) and orthogonal in Wn,div3,2 . Moreover, the projection Pk:L2(Ω)→span{wi}i=1k, defined as^[2]

is continuous in L²(Ω) and also in Wn,div3,2 independently of k, i.e.,

for all φ ∈ Wn,div3,2 , where the constant C is independent of k. Furthermore, by the standard embedding, we also have that Wn,div3,2 ↪ W^2,6(Ω) ↪ W^1,∞(Ω). Similarly, we construct the basis {Wj}j=1∞ of W^1,2(Ω), which is L²-orthonormal, W^1,2-orthogonal and the projection

is continuous in L²(Ω) and in W^1,2(Ω) independently of l.

Then for fixed k, l ∈ ℕ and ε ∈ (0, 1), we look for the functions vεk,l,Bεk,l of the form

where cik,l,ε,djk,l,ε,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 cik,l,ε(t) and djk,l,ε(t)) satisfy the following system of (k+l) ordinary differential equations in time interval (0, T):

Due to the L²-orthonormality of the bases {wi}i=1∞ and {Wj}j=1∞, the system (3.1)–(3.2) can be rewritten as a nonlinear system of ordinary differential equations for cik,l,ε and djk,l,ε, where i = 1, …, k and j = 1, …, l, and we equip this system with the initial conditions

Here, B0ε is defined as

Since 𝔹₀(x) ∈ R>03×3 for almost every x ∈ Ω, we have that Λmbda(𝔹₀(x)) > 0 for almost all x ∈ Ω. Consequently, using the fact 𝔹₀ ∈ L²(Ω), we obtain, as ε → 0₊, that

Note also that the initial conditions (3.3) can be rewritten as vεk,l(0)=Pkv0 and Bεk,l(0)=QlB0ε.

For the system (3.1)–(3.3), Carathéodory’s theorem can be applied and therefore there exists T^* > 0 and absolutely continuous functions cik,l,ε,djk,l,ε satisfying (3.3) and (3.1)–(3.2) almost everywhere in (0,T^*). If T^* is the maximal time, for which the solution exists, and T^* < T, then at least one of the functions cik,l,ε,djk,l,ε must 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 (vl,Bl):=(vεk,l,Bεk,l). We start by proving estimates independent of l. Since 𝔹_l(t) and v_l(t) belong for almost all t to the linear hull of {Wj}j=1l and {wi}i=1k, respectively, we can use v_l instead of w_i in (3.1) and 𝔹_l instead of 𝕎_j in (3.2) to deduce,

where we used the integration by parts formula and the facts that div v_l = 0 and 𝓣v ⋅ n = 0. Next, it follows from the definition of ρ_ε, ℝ and 𝕊 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

Next, we focus on the estimate for time derivatives. First, it follows from L²-orthonormality of the bases and the estimate (3.7) that

Then, since v_l is a linear combination of {wi}i=1k⊂W1,∞(Ω), we can estimate

and we can deduce from (3.1) and f∈L2(0,T;Wn,div−1,2) that

Finally, it follows from (3.2) and (3.7) that

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