A Hybrid High-Order method for creeping flows of non-Newtonian fluids

In this paper, we design and analyze a Hybrid High-Order discretization method for the steady motion of non-Newtonian, incompressible fluids in the Stokes approximation of small velocities. The proposed method has several appealing features including the support of general meshes and high-order, unconditional inf-sup stability, and orders of convergence that match those obtained for scalar Leray-Lions problems. A complete well-posedness and convergence analysis of the method is carried out under new, general assumptions on the strain rate-shear stress law, which encompass several common examples such as the power-law and Carreau-Yasuda models. Numerical examples complete the exposition.


Introduction
In this paper, we design and analyze a Hybrid High-Order (HHO) discretization method for the steady motion of a non-Newtonian, incompressible fluid in the Stokes approximation of small velocities. Notable applications include ice sheet dynamics [32 32], mantle convection [44 44], chemical engineering [33 33], and biological fluids rheology [28 28, 37 37]. We focus on fluids with shear-rate-dependent viscosity, whose behavior is characterized by a nonlinear strain rate-shear stress function. Physical interpretations and discussions of non-Newtonian fluid models can be found, e.g., in [8 8, 40 40]. Typical examples that are frequently used in the applications include the power-law and Carreau-Yasuda model, covered by the present analysis.
The earliest investigations of fluids with shear-dependent viscosity date back to the pioneering work of Ladyzhenskaya [36 36]. For a detailed mathematical study of the well-posedness and regularity of the continuous problem, see also [3 3, 7 7, 24 24, 39 39, 41 41] and references therein. Early results on the numerical analysis of non-Newtonian fluid flow problems were given in [2 2, 30 30, 43 43]. Later, these results were improved in [6 6] and [31 31] by proving error estimates that are optimal for fluids with shear thinning behavior (described by a power-law exponent ≤ 2). In [6 6], the authors considered a conforming inf-sup stable finite element discretization, while in [31 31] a low-order scheme with local projection stabilization was proposed. In both works, the use of Orlicz functions is instrumental to unify the treatment of the shear thinning and shear thickening cases (also called pseudoplastic and dilatant, respectively; cf. Example 4 4).
: Ω → R . Its flow is governed by the generalized Stokes problem, which consists in finding the velocity field : Ω → R and the pressure field : Ω → R such that −∇· (·, ∇ s ) + ∇ = in Ω, (1a) where ∇· denotes the divergence operator applied to vector or tensor fields, ∇ s is the symmetric part of the gradient operator ∇ applied to vector fields, and, denoting by R × s the set of square, symmetric, real-valued × matrices, : Ω × R × s → R × s is the strain rate-shear stress law. In what follows, we formulate assumptions on that encompass common models for non-Newtonian fluids and state a weak formulation for problem (1 1) that will be used as a starting point for its discretization.
Assumption 1 (Strain rate-shear stress law). Let a real number ∈ (1, ∞) be fixed, denote by −1 ∈ (1, ∞) the conjugate exponent of , and define the singular exponent of bỹ ≔ min( , 2) ∈ (1, 2]. ( The strain rate-shear stress law satisfies : Moreover, there exist real numbers de ∈ [0, ∞) and hc , sm ∈ (0, ∞) such that, for all , ∈ R × s and almost every ∈ Ω, we have the Hölder continuity property These relations are reminiscent of the ones used in [17 17] in the context of scalar Leray-Lions problems. The advantage of assumptions (3c 3c)-(3d 3d), expressed in terms of the singular index˜, is that they enable a unified treatment of the cases < 2 and ≥ 2 in the proofs of Lemma 18 18, Theorem 11 11, Lemma 20 20, and Theorem 12 12 below.
Indeed, let ∈ R × s be such that | | × > 0. Using the strong monotonicity (3d 3d) (with = 0), the Cauchy-Schwarz inequality, and the Hölder continuity (3c 3c) (again with = 0), we infer that where : Ω → [ − , + ] is a measurable function with − , + ∈ (0, ∞) corresponding to the local flow consistency index, ∈ [0, ∞) is the degeneracy parameter, : Ω → [ − , + ] is a measurable function with − , + ∈ (0, ∞) expressing the local transition flow behavior index, and ∈ (1, ∞) is the flow behavior index. The Carreau-Yasuda law is a generalization of the Carreau law (corresponding to − = + = 2) that takes into account the different local levels of flow behavior in the fluid. The degenerate case = 0 corresponds to the power-law model. Non-Newtonian fluids described by constitutive laws with a ( , , , )-structure exhibit a different behavior according to the value of . If > 2, then the fluid shows shear thickening behavior and is called dilatant. Examples of dilatant fluids are wet sand and oobleck. The case < 2, on the other hand, corresponds to pseudoplastic fluids having shear thinning behavior, such as blood. Finally, if = 2, then the fluid is Newtonian and (1 1) becomes the classical (linear) Stokes problem. We show in Appendix A A that the strain rate-shear stress law (5 5) is an -power-framed function with de = , if ≥ 2, where ⊕ ≔ max(0, ) and ≔ − min(0, ) denote, respectively, the positive and negative parts of a real number . As a consequence, it matches Assumption 1 1.

Weak formulation
From this point on, we omit both the integration variable and the measure from integrals, as they can be in all cases inferred from the context. We define the following velocity and pressure spaces embedding, respectively, the homogeneous boundary condition and the zero-average constraint: Assuming ∈ (Ω, R ), the weak formulation of problem (1 1) reads: Find ( , ) ∈ × such that where the function : × → R and the bilinear form : × (Ω, R) → R are defined such that, for all , ∈ and all ∈ (Ω, R), Remark 5 (Mass equation). The test space in (6b 6b) can be extended to (Ω, R) since, for all ∈ , the divergence theorem and the fact that with Ω denoting the unit vector normal to Ω and pointing out of Ω. Remark 6 (Well-posedness and a priori estimates). It can be checked that, under Assumption 1 1, the continuous problem (6 6) admits a unique solution ( , ) ∈ × ; see, e.g., [31 31,Section 2.4], where slightly stronger assumptions are considered. For future use, we also note the following a priori bound on the velocity: where K > 0 comes from the Korn inequality given at (33 33) below. To prove (8 8), use the strongmonotonicity (3d 3d) of , sum (6a 6a) written for = to (6b 6b) written for = , and use the Hölder inequality together with the Korn inequality (33 33) to write where |Ω| is the measure of Ω, that is, Observing that ∇ s , we obtain, enumerating the cases for the maximum and summing the corresponding bounds, ∇ s (Ω,R × ) ≤ (2 2−˜N ) Combining this inequality with (9 9) gives (8 8).

Mesh and notation for inequalities up to a multiplicative constant
We define a mesh as a couple M ℎ ≔ (T ℎ , F ℎ ), where T ℎ is a finite collection of polyhedral elements such that ℎ = max ∈ T ℎ ℎ with ℎ denoting the diameter of , while F ℎ is a finite collection of planar faces with diameter ℎ . Notice that, here and in what follows, we use the three-dimensional nomenclature also when = 2, i.e., we speak of polyhedra and faces rather than polygons and edges. It is assumed henceforth that the mesh M ℎ matches the geometrical requirements detailed in [19 19,Definition 1.7]. In order to have the boundedness property (14 14) for the interpolator, we additionally assume that the mesh elements are star-shaped with respect to every point of a ball of radius uniformly comparable to the element diameter; see [19 19,Lemma 7.12] for the Hilbertian case. Boundary faces lying on Ω and internal faces contained in Ω are collected in the sets F b ℎ and F i ℎ , respectively. For every mesh element ∈ T ℎ , we denote by F the subset of F ℎ containing the faces that lie on the boundary of . For every face ∈ F ℎ , we denote by T the subset of T ℎ containing the one (if ∈ F b ℎ ) or two (if ∈ F i ℎ ) elements on whose boundary lies. Finally, for each mesh element ∈ T ℎ and face ∈ F , denotes the (constant) unit vector normal to pointing out of .
Our focus is on the ℎ-convergence analysis, so we consider a sequence of refined meshes that is regular in the sense of [19 19, Definition 1.9] with regularity parameter uniformly bounded away from zero. The mesh regularity assumption implies, in particular, that the diameter of a mesh element and those of its faces are comparable uniformly in ℎ and that the number of faces of one element is bounded above by an integer independent of ℎ.
To avoid the proliferation of generic constants, we write henceforth (resp., ) for the inequality ≤ (resp., ≥ ) with real number > 0 independent of ℎ, of the constants de , hc , sm in Assumption 1 1, and, for local inequalities, of the mesh element or face on which the inequality holds. We also write to mean and . The dependencies of the hidden constants are further specified when needed.

Discrete spaces and norms
Let an integer ≥ 1 be fixed. The HHO space of discrete velocity unknowns is The interpolation operator ℎ : 1,1 (Ω, R ) → ℎ maps a function ∈ 1,1 (Ω, R ) on the vector of discrete unknowns ℎ defined as follows: For all ∈ T ℎ , we denote by and the restrictions of ℎ and ℎ to , respectively and, for all ℎ ∈ ℎ , we let ( , ( ) ∈ F ) ∈ denote the vector collecting the discrete unknowns attached to and its faces. Furthermore, for all ℎ ∈ ℎ , we define the broken polynomial field ℎ ∈ P (T ℎ , R ) obtained patching element unknowns, that is, We define on ℎ the 1, (Ω, R )-like strain seminorm · ,ℎ such that, for all ℎ ∈ ℎ , The following boundedness property for can be proved adapting the arguments of [19 19,Proposition 6.24] and requires the star-shaped assumption on the mesh elements: For all ∈ T ℎ and all ∈ 1, ( , R ), where the hidden constant depends only on , the mesh regularity parameter, , and .
The discrete velocity and pressure are sought in the following spaces, which embed, respectively, the homogeneous boundary condition for the velocity and the zero-average constraint for the pressure: By the discrete Korn inequality proved in Lemma 15 15 below, · ,ℎ is a norm on ℎ,0 (the proof is obtained reasoning as in [19 19, Corollary 2.16]).

HHO scheme
In this section, after introducing the discrete counterparts of the viscous and pressure-velocity coupling terms, we state the discrete problem along with the main results.

Local symmetric gradient reconstruction
For all ∈ T ℎ , we define the local symmetric gradient reconstruction G s, : The global symmetric gradient reconstruction G s,ℎ : is obtained patching the local contributions, that is, for all ℎ ∈ ℎ , we set

Discrete viscous function
The discrete counterpart of the function defined by (7 7) is a ℎ : In the above definition, recalling (4 4), is a stabilization parameter such that while the stabilization function s ℎ : where the local contributions are assumed to satisfy the following assumption.

An example of viscous stabilization function
Taking inspiration from the scalar case (cf., e.g., [18 18, Eq. (4.11c)]), a local stabilization function that matches Assumption 2 2 can be obtained setting, for all , ∈ , where, denoting by P (F , R ) the space of vector-valued broken polynomials of total degree ≤ on F , the boundary residual operator : with velocity reconstruction r +1 : Above, ∇ ss denotes the skew-symmetric part of the gradient operator ∇ applied to vector fields and ⊗ is the tensor product such that, for all = ( (25 25)). The local stabilization function defined by (25 25) satisfies Assumption 2 2.

Pressure-velocity coupling
For all ∈ T ℎ , we define the local divergence reconstruction D : We have the following characterization of D : For all ∈ , as can be checked writing (15 15) for = I . Taking the trace of (16 16), it is inferred that, for all ∈ T ℎ and all ∈ 1,1 ( , R ), D ( ) = (∇· ). The pressure-velocity coupling is realized by the bilinear form b ℎ :

Discrete problem and main results
The discrete problem reads: Remark 9 (Discrete mass equation). The space of test functions in (29b 29b) can be extended to P (T ℎ , R) since, for all ℎ ∈ ℎ,0 , the divergence theorem together with the fact that = 0 for all ∈ F b ℎ and Remark 10 (Efficient implementation). When solving the system of nonlinear algebraic equations corresponding to (29 29) by, e.g., the Newton algorithm, all element-based velocity unknowns and all but one pressure unknown per element can be locally eliminated at each iteration by static condensation. As all the computations are local, this procedure is an embarrassingly parallel task which can fully benefit from multi-thread and multi-processor architectures. This implementation strategy has been described for the linear Stokes problem in [22 22, Section 6.2]. After further eliminating the boundary unknowns by strongly enforcing the boundary condition (1c 1c), we end up solving, at each iteration of the nonlinear solver, a linear system of size card(F i ℎ ) + −1 −1 + card(T ℎ ). Concerning the interplay between the static condensation strategy and the performance of -multilevel linear solvers, we refer to [11 11].
In what follows, we state the main results for the HHO scheme (29 29). The proofs are postponed to Section 7 7.
Theorem 11 (Well-posedness). There exists a unique solution ( ℎ , ℎ ) ∈ ℎ,0 × ℎ to the discrete problem (29 29). Additionally, the following a priori bounds hold: Proof. See Section 7.2 7.2. (6 6) and (29 29), respectively. Assume the additional regularity Then, under Assumptions 1 1 and 2 2, where we have set, for the sake of brevity, Notice that, owing to the presence of higher-order terms in the right-hand sides of (31 31), higher convergence rates may be observed before attaining the asymptotic ones; see Section 5 5. The asymptotic order of convergence for the velocity coincides with the one proved in [17 17, Theorem 3.2] for HHO discretizations of scalar Leray-Lions problems. We refer to [20 20] for recent improvements on these estimates depending on the degeneracy of the problem.

Numerical examples
In this section, we evaluate the numerical performance of the HHO method on a complete panel of numerical test cases. We focus on the ( , 0, 1, )-Carreau-Yasuda law (5 5)

Trigonometric solution
We begin by considering a manufactured solution to problem (1 1) in order to assess the convergence of the method. We take Ω = (0, 1) 2 and exact velocity and pressure given by, respectively, The volumetric load and the Dirichlet boundary condition are inferred from the exact solution. Considering = 1 and ∈ {1.5, 1.75, . . . , 2.75}, this solution matches the assumptions required in Theorem 12 12 for = 1, except the case = 1.5 for which (·, ∇ s ) ∉ 1, (Ω, R × s ). We consider the HHO scheme for = 1 on three mesh families, namely Cartesian orthogonal, distorted triangular, and distorted Cartesian; see Figure 1 1. Overall, the results are in agreement with the theoretical predictions, and in some cases the expected asymptotic orders of convergence are exceeded. Specifically, for ≠ 2, the convergence rates computed on the last refinement surpass in some cases the theoretical ones. As noticed in Remark 13 13, this suggests that the asymptotic order is still not attained. A similar phenomenon has been observed on certain meshes for the -Laplace problem; see [17 17, Section 3.5.2] and [21 21, Section 3.7]. In some cases, we observe a better convergence for the velocity on distorted triangular meshes than on Cartesian meshes. This phenomenon possibly results from the combination of two factors: on one hand, the improved robustness of HHO methods with respect to elongated elements when compared to classical discretization methods; on the other hand, the fact that unstructured triangular meshes have more elements than Cartesian meshes for a given meshsize and lack privileged directions, which reduces mesh bias. Further investigation is postponed to a future work.

Lid-driven cavity flow
We next consider the lid-driven cavity flow, a well-known problem in fluid mechanics. The domain is the unit square Ω = (0, 1) 2 , and we enforce a unit tangential velocity = (1, 0) on the top edge (of equation 2 = 1) and wall boundary conditions on the other edges. This boundary condition is incompatible with the formulation (6 6), even generalized to non-homogeneous boundary conditions, since ∉ 1, (Ω, R ). However, this is a very classical test that demonstrates the quality of the method. We consider a low Reynolds number Re ≔ 2 = 1. For ∈ {1.25, 2, 2.75}, we solve the discrete problem on Cartesian and distorted triangular meshes (cf. Figure 1 1) of approximate size 128 × 128 for = 1, and 16 × 16 for = 5. This choice is meant to compare the low-order version of the method on a fine mesh with the high-order version on a very coarse one. The corresponding total number of degrees of freedom is: 130048 for the fine Cartesian mesh with = 1; 5760 for the coarse Cartesian mesh with = 5; 298676 for the fine triangular mesh with = 1; and 14196 for the coarse triangular mesh with = 5. In the left column of Figure 4 4 we display the velocity magnitude, while in the right column we plot the horizontal component 1 of the velocity along the vertical centreline 1 = 1 2 (resp., vertical component 2 along the horizontal centreline 2 = 1 2 ). The lines corresponding to = 1 on the fine mesh and to = 5 on the coarse mesh are perfectly superimposed, regardless of the mesh family and of the value of . This shows that, despite the lack of regularity of the exact solution, high-order versions of the scheme on very coarse meshes deliver similar results as low-order versions on very fine grids. Furthermore, we observe significant differences in the behavior of the flow according to , coherent with the expected physical behavior. In particular, the viscous effects increase with , as reflected by the size of the central vortex.

Discrete Korn inequality
We prove in this section a discrete counterpart of the following Korn inequality (see [29 29, Theorem 1]) that will be needed in the analysis: There is K > 0 depending only on Ω and such that for all ∈ , We start by recalling the following preliminary result concerning the node-averaging interpolator (sometimes called Oswald interpolator). Let ℎ be a matching simplicial submesh of M ℎ in the sense of [19 19, Definition 1.8]. The node-averaging operator av,ℎ : P (T ℎ , R ) → P ( ℎ , R ) ∩ 1, (Ω, R ) is such that, for all ℎ ∈ P (T ℎ , R ) and all Lagrange node of ℎ , denoting by the set of simplices sharing , For all ∈ F i ℎ , denote by 1 , 2 ∈ T ℎ the elements sharing , taken in an arbitrary but fixed order. We define the jump operator such that, for any function This definition is extended to boundary faces ∈ F b ℎ by setting [ ] ≔ | . Proposition 14 (Boundedness of the node-averaging operator). For all ℎ ∈ P (T ℎ , R ), it holds where F V, collects the faces whose closure has non-empty intersection with . Using the local inverse where we have used the fact that ℎ − ≤ ℎ − along with inequality (35 35) to pass to the second line, and we have exchanged the sums after setting T V, ∈ T ℎ : ∩ ≠ ∅ for all ∈ F ℎ to pass to the third line. Observing that max ∈ F ℎ card(T V, ) 1 (since, for any ∈ F ℎ , card(T V, ) is bounded by the left-hand side of [19 19, Eq. (4.23)] written for any ∈ T ℎ to which belongs), (34 34) follows.
The following inequalities between sums of powers will be often used in what follows without necessarily recalling this fact explicitly each time. Let an integer ≥ 1 and a real number ∈ (0, ∞) be given. Then, for all 1 , . . . , ∈ (0, ∞), we have . Gathering the above cases yields (36 36).

Lemma 15. (Discrete Korn inequality)
We have, for all ℎ ∈ ℎ,0 , recalling the notation (12 12), Proof. Let ℎ ∈ ℎ,0 . Using a triangle inequality followed by (36 36), we can write where we have used the continuous Korn inequality (33 33) to pass to the second line, we have inserted ±∇ s,ℎ ℎ into the first norm and used a triangle inequality followed by (36 36) to pass to the third line, and we have invoked the bound (34 34) to conclude. Observing that, for any ∈ F ℎ , | [ ℎ ] | ≤ ∈ T | − | by a triangle inequality, and using (36 36), we can continue writing where we have exchanged the sums over faces and elements and recalled definition (13a 13a) to conclude. This proves the bound for the second term in the left-hand side of (37 37

Well-posedness and convergence analysis
In this section, after studying the stabilization function s ℎ , we prove the main results stated in Section 4.3 4.3.

Properties of the stabilization function
Lemma 16 (Consistency of s ). For any ∈ T ℎ and any s satisfying Assumption 2 2, it holds, for all ∈ +2, ( , R ) and all ∈ , where the hidden constant is independent of ℎ, , and .
Proof. The proof adapts the arguments of [19 19, Propositon 2.14]. Using the polynomial consistency property (22 22), we can write where we have used the Hölder continuity (23 23) and observed that, by the consistency property (22 22), s ( ( +1 ), ( +1 )) = 0 to pass to the second line, we have used the boundedness property (21 21) to pass to the third line, the boundedness (14 14) of to pass to the fourth line, and the ( + 2, , 1)approximation property (11a 11a) of +1 to conclude.
In what follows, we will need generalized versions of the continuous and discrete Hölder inequalities, recalled hereafter for the sake of convenience. Let ⊂ R be measurable, ∈ N * , and let , 1 , . . . , ∈ (0, ∞] be such that Proposition 17 (Properties of s ℎ ). Let s ℎ be given by (20 20) with, for all ∈ T ℎ , s satisfying Assumption 2 2. Then it holds, for all ℎ ∈ ℎ , Finally, for any ∈ ∩ +2, (T ℎ , R ), it holds Above, the hidden constants are independent of ℎ and of the arguments of s ℎ .

Well-posedness
In this section, after proving Hölder continuity and strong monotonicity properties for the discrete viscous function a ℎ and the inf-sup stability of the pressure-velocity coupling bilinear form b ℎ , we prove Theorem 11 11.

Stability of the pressure-velocity coupling
Lemma 19 (Inf-sup stability of b ℎ ). It holds, for all ℎ ∈ ℎ , with hidden constant depending only on , , , Ω, and the mesh regularity parameter.
(i) Fortin operator. We need to prove that the following properties hold for any ∈ 1, (Ω, R ): Property (49a 49a) is obtained by raising both sides of (14 14) to the power , summing over ∈ T ℎ , then taking the th root of the resulting inequality. The proof of (49b 49b) is given, e.g., in [19 19,Lemma 8.12].

Proof of Theorem 11 11
Proof of Theorem 11 11. (i) Existence. Denote by , * ℎ the dual space of ℎ and let ℎ : ℎ,0 → , * ℎ be such that, for all ℎ ∈ ℎ,0 , Here and in what follows, ·, · denotes the appropriate duality pairing as inferred from its arguments. Define the following subspace of ℎ,0 spanned by vectors of discrete unknowns with zero discrete divergence: and consider the following problem: Existence of a solution to this problem for a fixed ℎ can be proved adapting the arguments of [18 18,Theorem 4.5]. Specifically, equip ℎ with an inner product (·, ·) ,ℎ (which need not be further specified), denote by · ,ℎ the induced norm, and let ℎ : ℎ → ℎ be such that, for all ℎ ∈ ℎ , The strong monotonicity (43b 43b) of a ℎ yields, for any ℎ ∈ ℎ such that ℎ ,ℎ ≥ de , where denotes the constant (possibly depending on ℎ) in the equivalence of the norms · ,ℎ and · ,ℎ (which holds since ℎ is finite-dimensional). This shows that ℎ is coercive hence, by [16 16, Theorem 3.3], surjective. Let now ℎ ∈ ℎ be such that ( ℎ , ℎ ) ,ℎ = ∫ Ω · ℎ for all ℎ ∈ ℎ . By the surjectivity of ℎ , there exists ℎ ∈ ℎ such that ℎ ( ℎ ) = ℎ which, by definition of ℎ and ℎ , is a solution to the discrete problem (53 53).

Error estimate
In this section, after studying the consistency of the viscous and pressure-velocity coupling terms, we prove Theorem 12 12.

Consistency of the viscous function
Lemma 20 (Consistency of a ℎ ). Let ∈ ∩ +2, (T ℎ , R ) be such that (·, ∇ s ) ∈ 1, (Ω, R × s ) ∩ ( +1) (˜−1), (T ℎ , R × s ). Define the viscous consistency error linear form E a,ℎ ( ; ·) : ℎ → R such that, for all ℎ ∈ ℎ , Then, under Assumptions 1 1 and 2 2, we have Proof. Letˆℎ ℎ and ℎ ∈ ℎ,0 . Expanding a ℎ according to its definition (18 18) in the expression (55 55) of E a,ℎ , inserting ± ∫ Ω (·, ∇ s ) : G s,ℎ ℎ + ∫ Ω ℎ (·, ∇ s ) : G s,ℎ ℎ , and rearranging, we obtain where have used the definition (10 10) of ℎ together with the fact that G s,ℎ ℎ ∈ P (T ℎ , R × s ) in the cancellation. We proceed to estimate the terms in the right-hand side. For the first term, we start by noticing that ∑︁ as a consequence of the continuity of the normal trace of (·, ∇ s ) together with the single-valuedness of across each interface ∈ F i ℎ and of the fact that = 0 for every boundary face ∈ F b ℎ . Using an element by element integration by parts on the first term of T 1 along with the definitions (17 17 where we have used the definition (10 10) of ℎ together with the fact that ∇ s,ℎ ℎ ∈ P −1 ( ) to cancel the term in the first line, and we have inserted (58 58) and rearranged to conclude. Therefore, applying the Hölder inequality together with the bound ℎ ≤ ℎ , we infer where the conclusion follows using the (( + 1) (˜− 1), )-trace approximation properties (11b 11b) of along with ℎ ≤ ℎ for the first factor and the definition (13 13) of the · ,ℎ -norm for the second.
For the second term, using the Hölder inequality and again (41a 41a), we get We estimate the first factor as follows: where we have used the Hölder continuity (3c 3c) of in the first bound, the ( ; −˜,˜−1 )-Hölder inequality (39 39) in the second, the boundedness of Ω along with (41a 41a) and the commutation property (16 16) of G s,ℎ in the third, and we have concluded invoking the ( + 1, , 0)-approximation property (11a 11a) of . Plugging this estimate into (60 60), we get Finally, using the fact that ≤ hc together with the consistency (42 42) of s ℎ and the norm equivalence (41a 41a), we obtain for the third term Plug the bounds (59 59), (61 61), and (62 62) into (57 57) and pass to the supremum to conclude.

Consistency of the pressure-velocity coupling bilinear form
Lemma 21 (Consistency of b ℎ ). Let ∈ 1, (Ω, R) ∩ ( +1) (˜−1), (T ℎ , R). Let E b,ℎ ( ; ·) : ℎ → R be the pressure consistency error linear form such that, for all ℎ ∈ ℎ , Then, we have that Proof. Let ℎ ∈ ℎ,0 . Integrating by parts element by element, we can reformulate the first term in the right-hand side of (63 63) as follows: where the introduction of in the boundary term is justified by the fact that the jumps of vanish across interfaces by the assumed regularity and that = 0 on every boundary face ∈ F b ℎ . On the other hand, expanding, for each ∈ T ℎ , D according to its definition (27 27), we get Summing (65 65) and (66 66) and observing that the first terms in parentheses cancel out by the definition (10 10) of since ∇· ∈ P −1 ( , R) ⊂ P ( , R) for all ∈ T ℎ , we can write where we have used the Hölder inequality along with ℎ ≥ ℎ whenever ∈ F in the second line and the (( + 1)(˜− 1), )-trace approximation property (11b 11b) of together with the bound ℎ ≤ ℎ and the definition (13 13) of the · ,ℎ -norm to conclude. Passing to the supremum yields (64 64).
Step 3. Error estimate for the pressure. Using the Hölder continuity (43a 43a) of a ℎ , we have, for all where the first factor is estimated as in (69 69). Thus, using the inf-sup condition (48 48), we can write where we have used the definition (67 67) of the consistency error together with equation (29a 29a) to pass to the second line, (71 71) to pass to the third line (recall that $ denotes here the supremum in the left-hand side of (68 68)), and the bounds (68 68) and (31a 31a) (proved in Step 2) to conclude.

A Power-framed functions
In the following theorem, we introduce the notion of power-framed function and discuss sufficient conditions for this property to hold.
Remark 23 (Notation). The boldface notation for the elements of is reminiscent of the fact that Theorem 22 22 is used with = R × s in Corollary 24 24 to characterize the Carreau-Yasuda law as an -power-framed function and in Lemma 8 8 with = R to study the local stabilization function s . Theorem 22 22. Let ∈ be such that (73 73) holds, and , ∈ . By symmetry of inequalities (74 74) and the fact that is continuous, we can assume, without loss of generality, that > > 0. where, to pass to the second line, we have removed negative contributions if < 2 and used the fact that ( − ) −1 ≤ de + + if ≥ 2, to pass to the third line we have used the fact that ↦ → −2 is non-increasing if < 2, and the fact that ≤ otherwise, while the conclusion follows from the definition of sm . This shows that is non-decreasing. Hence, for all ∈ [ , ∞), ( ) ≥ ( ) = 0, i.e.