Optimal error estimate for a space-time discretization for incompressible generalized Newtonian fluids: The Dirichlet problem

In this paper we prove optimal error estimates for {solutions with natural regularity} of the equations describing the unsteady motion of incompressible shear-thinning fluids. We consider a full space-time semi-implicit scheme for the discretization. The main novelty, with respect to previous results, is that we obtain the estimates directly without introducing intermediate semi-discrete problems, which enables the treatment of homogeneous Dirichlet boundary conditions.


Introduction
In this paper we study a space-time discretization of the unsteady system describing the motion of homogeneous (for simplicity the density ρ is set equal to 1), incompressible shear-thinning fluids under homogeneous Dirichlet boundary conditions. We prove optimal error estimates (cf. Section 2.4) for solutions possessing a natural regularity, extending the results in [5] to the case of homogeneous Dirichlet boundary conditions. Our method differs from most previous investigations in as much as we use no intermediate semi-discrete problems to prove our result. We restrict ourselves to the three-dimensional setting, however, all results can be easily adapted to the general setting in d-dimensions.
More precisely, we consider for a bounded polyhedral domain Ω ⊂ R 3 and a finite time interval I := (0, T ), for some given T > 0, the system ∂ t u − div S(Du) + [∇u]u + ∇q = f in I × Ω, div u = 0 in I × Ω, u(0) = u 0 in Ω, (1.1) where the vector field u = (u 1 , u 2 , u 3 ) ⊤ is the velocity, the scalar q is the kinematic pressure, the vector f = (f 1 , f 2 , f 3 ) ⊤ is the external body force and u 0 is the initial velocity. We assume that the extra stress tensor S has (p, δ)-structure for some p ∈ (1,2], and δ ∈ [0, ∞) (cf. Section 2.2). For the convective term we use the notation ([∇u]u) i = 3 j=1 u j ∂ j u i , i = 1, 2, 3, while Du := 1 2 (∇u + ∇u ⊤ ) denotes the symmetric part of the gradient ∇u. For smooth enough solutions (u, q) the variational formulation of (1.1) reads (∂ t u(t), v) + (S(Du(t)), Dv) + b(u(t), u(t), v) − (q(t), divv) = (f (t), v) ,  for all v ∈ V := W 1,p 0 (Ω) 3 , η ∈ Q := L p ′ 0 (Ω) and almost every t ∈ I. We used the notation b(u, v, w) : for the convective term to have a stable space-discretization. Note that b(·, ·, ·) is skew-symmetric with respect to the last two arguments, i.e., b(u, v, w) = −b (u, w, v) and that for solenoidal functions it holds b(u, v, w) = ([∇v]u, w). We perform an error analysis for the semi-implicit space-time discretization, which for given h > 0 and M ∈ N reads: for u 0 ) is the backward difference quotient with κ := T M . Here V h ⊂ V , Q h ⊂ Q are appropriate stable finite element spaces with mesh size h > 0. The precise setup can be found in Section 2.3.

Preliminaries and main results
In this section we introduce the notation, the setup, and recall some technical results which will be needed in the proof of the main result.
2.1. Function spaces. We use c, C to denote generic constants, which may change from line to line, but are not depending on the crucial quantities. We write f ∼ g if and only if there exist constants c, C > 0 such that c f ≤ g ≤ C f .
We will use the customary Lebesgue spaces (L p (Ω), . p ) and Sobolev spaces (W k,p (Ω), . k,p ), k ∈ N. We do not distinguish between scalar, vector-valued or tensor-valued function spaces in the notation if there is no danger of confusion. However, we denote vector-valued functions by small boldfaced letters and tensor-valued functions by capital boldfaced letters. If a norm is considered on a set M different than Ω we indicate this in the respective norms as . p,M , . k,p,M . We equip W 1,p 0 (Ω) with the gradient norm ∇ . p . We denote by |M | the 3-dimensional Lebesgue measure of a measurable set M . The mean value of a locally integrable function f over a measurable set M ⊂ Ω is denoted by f M := − M f dx = 1 |M|Ḿ f dx. By L p 0 (Ω) we denote the space of functions f ∈ L p (Ω) with f Ω = 0. Moreover, we use the notation (f, g) :=Ώ f g dx, whenever the right-hand side is well-defined. We use the following notation , for the most often used function spaces.

2.2.
Basic properties of the extra stress tensor. For a tensor P ∈ R 3×3 we denote its symmetric part by P sym := 1 2 (P+ P ⊤ ) ∈ R 3×3 sym := {A ∈ R 3×3 | P = P ⊤ }. The scalar product between two tensors P, Q is denoted by P · Q, and we use the notation |P| 2 = P·P. We assume that the extra stress tensor S has (p, δ)-structure, which will be defined now. A detailed discussion and full proofs of the following results can be found in [12,21].
sym ), satisfies S(P) = S P sym , and S(0) = 0. Moreover, we assume that S has (p, δ)-structure, i.e., there exist p ∈ (1, ∞), δ ∈ [0, ∞), and constants C 0 , C 1 > 0 such that are satisfied for all P, Q ∈ R 3×3 with P sym = 0 and all i, j, k, l = 1, . . . , 3. The constants C 0 , C 1 , and p are called the characteristics of S. Remark 2.3. We would like to emphasize that, if not otherwise stated, the constants in the paper depend only on the characteristics of S, but are independent of δ ≥ 0.
Another important tool are shifted N-functions 1 {ϕ a } a≥0 , cf. [12,13,21]. Defining for t ≥ 0 a special N-function ϕ by with some constants C i > 0, i = 0, 1. Next, the shifted N-functions are defined for t ≥ 0 by Note that ϕ a (t) ∼ (δ + a + t) p−2 t 2 and that the complementary function satisfies (ϕ a ) * (t) ∼ ((δ + a) p−1 + t) p ′ −2 t 2 . Moreover, the N-functions ϕ a and (ϕ a ) * satisfy the ∆ 2 -condition 2 uniformly with respect to a ≥ 0, i.e., ∆ 2 (ϕ a ) ≤ c 2 max{2,p} and ∆ 2 ((ϕ a ) * ) ≤ c 2 max{2,p ′ } , respectively. We will use also Young's inequality: for all ε > 0 there exists c ε > 0, such that for all s, t, a ≥ 0 it holds Closely related to the extra stress tensor S with (p, δ)-structure is the function In the following lemma we recall several useful results, which will be frequently used in the paper. The proofs of these results and more details can be found in [12,21,13,2].
The constants depend only on the characteristics of S.
(ii) For all ε > 0, there exist a constant c ε > 0 (depending only on ε > 0 and on the characteristics of S) such that for all u, v, w ∈ X we have with constants depending only on p.
The proof of the main result uses the following modification of Gronwall's lemma.
The following result will be used frequently in the sequel.
Proof. The assertion is proved in [7, Lemma 3.1] in the special case τ m = t m , m = 1, . . . , M . The general case follows exactly in the same way.
For the spatial discretization we denote by T h a family of shape-regular triangulations, consisting of 3-dimensional closed simplices K. We denote by h K the diameter of K and by ρ K the supremum of the diameters of inscribed balls. We assume that T h is non-degenerate, i.e., there exists a constant γ 0 > 0 such that max K∈T h hK ρK ≤ γ 0 . The global mesh-size h is defined by h := max K∈T h h K . Let S K denote the neighborhood of K, i.e., S K is the union of all simplices of T h intersecting K. By the assumptions we obtain that |S K | ∼ |K| and that the number of patches S K to which a simplex belongs are bounded uniformly in both h > 0 and K ∈ T h .
We denote by P k (T h ), with k ∈ N 0 := N ∪ {0}, the space of scalar or vectorvalued functions, which are polynomials of degree at most k on each K ∈ T h . Given a triangulation T h of Ω with the above properties and given r 0 , r 1 , s 0 ∈ N 0 , with r 0 ≤ r 1 , we define For the weak formulation of discrete problems we use the following function spaces We also need some numerical interpolation operators. Rather than working with a specific interpolation operator we make the following assumptions: 19. We assume that r 0 = 1 and that there exists a linear projection operator Π div h : X → X h which (a) is locally W 1,1 -stable, i.e., for all w ∈ X and K ∈ T h there holds Assumption 2.21. We assume that Y h contains the constant functions, i.e. that R ⊂ Y h , and that there exists a linear projection operator The existence of a projection operator Π div h as in Assumption 2.19 is known (among others) for the Taylor-Hood, the Crouzeix-Raviart, and the MINI element in dimensions two and three; the Clément interpolation operator satisfies Assumption 2.21. For a discussion and consequences of these assumptions we refer to [2, Sec. 3.2], [7,Appendix], and [17,Sec. 4,5]. We collect in the next two propositions the properties of the projection operators, which are relevant for the present paper. (i) Let F(Dv) ∈ W 1,2 (Ω). Then, there exists a constant c = c(p, r 1 , γ 0 ) such that (ii) Let r ∈ (1, ∞). Then, there exists a constant c = c(r, r 1 , γ 0 ) such that for all v ∈ W 1,r (Ω) and all w ∈ Proof. The first two assertions are proved, e.g., in [2,17]. The last two assertions are proved in [7].
Proposition 2.23. Let Π Y h satisfy Assumption 2.21. Let ψ be an N-function satisfying the ∆ 2 -condition. Then, there exists a constant c = c(γ 0 , ∆ 2 (ψ)) such that for all sufficiently smooth functions and all K ∈ T h there holdŝ Proof. The assertions are proved in [17].

2.4.
Main results. Before we formulate the main result, proving optimal convergence rates for the error between the solution u of the continuous problem (1.1) and the discrete solution {u m h } M m=0 of the space-time discretization (1.4), we discuss the existence of solutions of (1.1) and (1.4).

2.5.
Comparison with previous results and observation on the requested regularity. Here we compare the new result with previous ones on the discretization of generalized non-Newtonian fluids (and general parabolic equations and systems). We also discuss briefly the regularity we are assuming on the continuous solution.
Concerning previously proved error-estimates we can mainly recall the following facts: (i) Problem (1.1), in the case of space periodic boundary conditions, has been treated in [5,14,15,20]. In [5] the same optimal error estimate as in Theorem 2.24 is proved under slightly stronger assumptions on the regularity of the solution u of (1.1), for p ∈ ( 3 2 , 2]. Problem (1.1) without the convective term [∇u]u, in the case of homogeneous Dirichlet boundary conditions, has been treated in [19] for p ∈ ( 6 5 , ∞). It is shown there that (2.27) holds with the right-hand side replaced by c (h min{2, 4 p } + κ 2 ). The proofs of these results are based on intermediate semi-discrete problems, for which a certain regularity has to be proved, to obtain the desired optimal convergence rates. This in fact limits the results in [5,14,15,20] to the case of space periodic boundary conditions. Here, we avoid such (technical) problems by proving the error estimate directly without using intermediate semi-discrete problems. The same perspective is also taken in the recent papers [7,11]. In [11] problem (1.1) is treated for a nonlinear operator S depending on the full gradient ∇u, having (p(·, ·), δ)-structure, with a variable exponent p(·, ·), but without convective term and without the solenoidality condition. In this situation the error estimate (2.27) with the right-hand side replaced by c (h 2αx + κ 2αt ) is proved if the variable exponent belongs to C αx,αt (I × Ω) for appropriate α x , α t ∈ (0, 1] and an additional CFL-condition κ r ≤ c inf K∈T h h k , with some r ≥ 1+2αt 2αx , is satisfied. In [7] problem (1.1) is treated without convective term and without the solenoidality condition. There the error estimate (2.27) is proved under the same conditions as in Theorem 2.24.
(ii) In the recent paper [9] a different approach is used to treat the unsteady p-Laplace problem, i.e., problem (1.1) without convective term and without the solenoidality constraint. By using the L 2 -projection operator instead of the Scott-Zhang operator (cf. [22]), satisfying Assumption 2.19, the error estimate is proved in [9] without a coupling condition between h and κ. The removal of the coupling is due to the usage of the L 2 -projection which in the treatment of the time derivative does not produce terms needing a coupling (cf. estimate (3.4)). However, the treatment of the nonlinear operator S is subtle and requires delicate estimates, which result in the different error estimate (2.28) compared to the error estimate (2.27). The estimate (2.28) is proved under the assumption that F(∇u) ∈ L 2 (I; N αx,2 (Ω)) ∩ N αt,2 (I; L 2 (Ω)) u ∈ L 2 (I; N αx,2 (Ω)) for 1 2 < α t ≤ 1 and 0 < α x ≤ 1, where N α,2 are appropriate Nikol'skiȋ spaces with differentiability α ∈ (0, 1] (cf. [9]). Even for α x = α t = 1 estimate (2.28) differs from (2.27) since it contains time averages of the error rather than a point-wise error. The usage of limited regularity in the time-variable and time averaging of the error is motivated by a similar analysis performed for stochastic parabolic equations in [10].
(iii) In [23,24] the convergence of a fully implicit space-time discretization (without convergence rate but also with no assumptions of smoothness of the limiting problem) of the problem (1.1) in the case of homogeneous Dirichlet boundary conditions is proved for p > 11 5 and even for p > 6 5 for a regularized problem. The convergence of the same numerical scheme (1.4) towards a weak solution has been recently proved in [6] for p > 11 5 . In fact, in [6] the convergence of a general quasi non-conforming Rothe-Galerkin scheme in the context of evolution problems with Bochner pseudo-monotone operators is proved (cf. [1] for the treatment of a conforming Rothe-Galerkin scheme in the context of evolution problems with Bochner pseudo-monotone operators).
Let us now discuss the natural regularity assumption (2.25). The assumption (2.25) 1 is natural in the sense that it is the one obtained from the extra stress tensor S if formally tested with −∆u and ∂ 2 t u. The existence of solutions satisfying (2.25) 1 is proved rigorously for problem (1.1) in the case of periodic boundary conditions locally in time in [4,16], for p > 7 5 and large data. The situation for Dirichlet boundary conditions is more complicated. The existence of a locally in time unique strong solution u ∈ L r (I ′ ; W 2,r (Ω)) ∩ L p (I ′ ; V div ) ∩ W 1,r (I ′ ; L r (Ω)), q ∈ L r (I ′ ; W 1,r (Ω)) ∩ L r (I ′ ; L r 0 (Ω)) for any r ∈ (5, ∞) and some I ′ := (0, T ′ ) with T ′ ∈ (0, T ) is proved in [8] for large data. This regularity implies, using parabolic embedding theory (cf. [16,Appendix]), that F(Du) ∈ L 2 (I ′ ; W 1,2 (Ω)) and that u, ∇u ∈ C(I ′ × Ω). However, in [8] it is not proved that this solution also satisfies ∂ t F(Du) ∈ L 2 (I ′ ; L 2 (Ω)). Nevertheless, one can show that the solution from [8] also satisfies ∂ t F(Du) ∈ L 2 (I ′ ; L 2 (Ω)), using the following auxiliary result.
Proof. This result is proved using ideas from [4,16]. Using the Galerkin method and the theory of monotone operators one shows that there exists a unique weak solutions v ∈ L ∞ (I; L 2 (Ω)) ∩ L p (I; V div ) of (2.30). Moreover, the regularity of the data allows us to show that we can take the time derivative of the Galerkin equations and test with the time derivative of the Galerkin solution. Straightforward manipulations show that this produces an estimate, showing after a limiting process in the Galerkin parameter, the additional regularity stated above.
Next, by using Proposition 2.29 with G = u ⊗ u (where u is the solution from [8]), the monotonicity of S, and the above regularity for u imply that the solution from [8] satisfies also ∂ t F(Du) ∈ L 2 (I ′ ; L 2 (Ω)). Consequently, the unique solution (u, q) from [8] satisfies (2.25) with I replaced by I ′ .
It is useful to formulate the consequences of the assumption F(Du) ∈ W 1,2 (I×Ω) in terms of Bochner-Sobolev spaces. From [16,Thm. 33] and standard embedding results it follows that F(Du) ∈ C(I; L 3 (Ω)).

Proof of the main result
In this section we prove the error estimates from Theorem 2.24. To this end we need to derive an equation for the error and to use the discrete Gronwall lemma 2.10 together with approximation properties due to the regularity of the solution and the properties of the extra stress tensor S.
In the error equation we want to use the test function u m h − Π div h u(t m ), which belongs to the space V h (0). Thus, it is enough to consider test functions v h from V h (0) in (1.4). For such test functions we can replace the discrete pressure q m h by an arbitrary function from Q h . Thus, it follows from (1.4) 1 that for all v h ∈ V h (0), µ h ∈ Q h and m = 1, . . . , M there holds we get that for all v h ∈ V h (0) and m = 1, . . . , M , there holds which we can re-write also as follows: since we are averaging locally constant terms. We subtract from this equation the retarded averages over I m of equation (1.2) and obtain the equation for the error We now discuss and estimate the terms in (3.3) separately, to arrive finally to the estimate (3.25). First, note that the projection operator Π div h has the same properties as the operator P h considered in [7] and that the solution u of (1.2) and the solution treated in [7] possess exactly the same regularity. Thus, the first two terms on the left-hand side can be treated exactly as in [7]. Consequently, [7,Lemmas 3.7,3.9] yield the following estimates: and The term with the external force is treated slightly differently compared to [7,Lemma 3.10]. This is due to the fact that to apply Gronwall's lemma 2.10 we need an estimate involving the L p -norm of the gradient of the error. To shorten the notation in the following computations we denote the error for m = 1, . . . , M by dt.
Proof. Using that together with Hölder's inequality, Young's inequality and the embedding W 1,p (Ω) ֒→ L 2 (Ω), valid for p ≥ 6 5 , we get Thus, we can re-write the integrand in the term to be estimated in (3.8) as follows To the latter we add and subtract, in the order, the terms b(u m−1 (3.9) In view of the skew-symmetry of b(·, ·, ·) we have I m 1 (t) = 0, for all m = 1, . . . , M and t ∈ I m .
Let us now discuss the last term, namely the one including the pressure.    , Π div h u(t m )) an extension of the validity of the error estimate for p ≤ 3 2 with the present technique is impossible, even with further regularity assumptions on u.