1. Introduction
The construction and mathematical analysis of finite element approximations of models of nonNewtonian fluids has been a subject of active research in recent years. Some of the most general results in this direction concern the convergence of mixed finite element approximations of models of incompressible fluids with implicit constitutive laws relating the Cauchy stress tensor to the symmetric velocity gradient (cf.
[15], [33] and [18]). Motivated by the groundbreaking contributions of Cohen, Dahmen and DeVore [12, 13] and Binev, Dahmen and DeVore [8] concerning the convergence of adaptive algorithms for linear elliptic problems, progress, albeit much more limited in both scope and extent, has also been made on the analysis of adaptive finite element approximations of implicitly constituted nonNewtonian fluid flow models (cf. [26]).Upon decomposing the Cauchy stress tensor into its traceless part, called the deviatoric stress tensor or shearstress tensor, and its diagonal part, called the mean normal stress, models of incompressible fluids typically involve the velocity of the fluid, , its pressure, , and the shearstress tensor, . For Newtonian fluids the shearstress tensor is a scalar multiple of the symmetric velocity gradient. The finite element approximation of Newtonian fluids is therefore usually performed in the velocitypressure formulation. For nonNewtonian fluids on the other hand the situation is more involved, because the shearstress tensor exhibits nonlinear dependence as a function of the symmetric velocity gradient, and the functional relationship between the shearstress tensor and the symmetric velocity gradient may even be completely implicit and multivalued. For powerlaw fluids, such as the ones considered in this work, the shearstress tensor exhibits powerlaw type growth as a function of the symmetric velocity gradient, the simplest instance of which results in an Laplace type operator in the balance of linear momentum equation, with a powerlaw exponent ; for , corresponding to a Newtonian fluid, the operator is linear, the Laplace operator. From a mathematical point of view, in the presence of a convection term in the balance of linear momentum equation in the model, the lower the value of the more difficult the problem is to analyse. The existence of solutions for small values of was first proved in [19], where an Acerbi–Fusco type Lipschitz truncation was used in conjunction with Minty’s method from monotone operator theory; thus, weak solutions were shown to exist for in space dimensions.
Finite element approximations of problems with powerlaw rheology have been extensively studied, including stabilised (or variationalmultiscale) methods (cf. [10, 1], for example) and local discontinuous Galerkin methods (see, [27], for example). The relevant literature is vast and it is beyond the scope of this work to provide an exhaustive survey of the various contributions; the interested reader may wish to consult [29], for example. Concerning implicitlyconstituted models, in the recent papers [15, 33] the convergence of generic infsup stable velocity/pressurebased mixed finite element methods was proved for , while convergence for the full range, , was shown to be achievable only in the case the finite element methods where the velocity space consists of pointwise divergencefree functions. The reason for this dichotomy is that in the case of velocity approximations that are discretely divergencefree only, as is the case in generic infsup stable mixed finite element methods, the finite element approximation of the convection term does not vanish when tested with
, and it needs to be skewsymmetrized (cf.
[34]) for this to happen. While in the case of the Navier–Stokes equations (corresponding to ) membership of the velocity field to the natural function space for weak solutions, , ensures that the convection term and its skewsymmetric modification can be bounded by the same expression using Hölder’s inequality, this is not the case for the powerlaw model under consideration here for entire range for which weak solutions to the problem are known to exist. In fact, in the case of nonNewtonian powerlaw models the natural function space for the velocity field is , and while the original convection term can be bounded in terms of the norm for all , for the skewsymmetric modification of the convection term, whose use is essential so as to be able to derive an energy inequality for discretely divergencefree velocity fields, this can only be achieved for the limited range . This was precisely the bottleneck encountered in [15] for discretely divergencefree velocity approximations, resulting in the reduction of the range of from the maximal range for which weak solutions are known to exist, toThe advantage of pointwise divergencefree finite element methods over discretely divergencefree finite element methods is therefore that, besides the physical consistency they provide, there is no need to rewrite the convection term in a skewsymmetric form. The topic of divergencefree finite element spaces has been treated extensively in the literature, most commonly presenting pairs of spaces for which the divergence of the velocity space is a subspace of, or equal to, the pressure space. For example, the early Scott–Vogelius element [32] (analysed recently in [24]) uses conforming piecewise polynomials of degree for the velocity, while discontinuous piecewise polynomials of degree are used for the pressure. The stability of this pair requires either special meshes, or a highenough degree (for example, is needed in [24], and is only allowed in very special cases such as those described in Remark 3). Another possibility is to relax the continuity requirements and consider a discontinuous Galerkin method, as was done, for example, in [11], or to relax only the tangential continuity of the approximate velocity on faces of elements while still preserving its continuity in the direction of the normal to faces of elements, thus using conforming methods, as was the case in [31], for example. In this latter case the viscous term (defined as the divergence of the shearstress needs to be modified, for stability reasons, by adding terms controlling the jumps and averages of the velocity into the formulation, with, obviously undesirable, extra complications if the viscous term in the balance of linear momentum equation has a more complex structure, as is the case for the powerlaw model considered herein.
The recent works [2, 3] offer a way of preserving the advantages of a pointwise divergencefree approximation to the velocity field while working with the, computationally simplest, lowestorder conforming velocity/pressure pair, namely, . The key idea in those works can be summarised as follows: the discrete continuity equation contains a stabilising term based on the jumps of the discrete pressure. As the jumps of the pressure are constant along element faces, there exists a unique Raviart–Thomas field such that its normal component is equal to the jumps. This field can be built at no extra computational cost, and then the continuity equation can be rewritten as a standard continuity equation, but for a modified velocity field, which is now solenoidal. The finite element method then involves replacing the original discrete velocity field with the new, now solenoidal, modified velocity field in the convection term. This facilitates the proofs of stability and convergence of the resulting finite element method without the need to rewrite the convection term in a skewsymmetric form. Our aim here is to apply this idea to a problem in nonNewtonian fluid mechanics. As a first step in this direction, we have chosen an explicit constitutive law with powerlaw rheology. Even though this is the simplest constitutive law, it has been shown experimentally to faithfully reproduce many situations of physical interest (see the discussion in [22], and the experimental results in, e.g., [25]); we therefore believe that it is a representative model for exemplifying the applicability of the proposed method in a mathematically nontrivial and physically relevant setting. Since the convection term does not need to be rewritten in a skewsymmetric form, the resulting method can now be proved to be stable and convergent to a weak solution for the whole range of the powerlaw index for which weak solutions to the model are known to exist. In addition, the sequence of numerical approximations is shown to converge strongly, and this strong convergence result is, to the best of our knowledge, a new contribution even in the, very special, Newtonian case ().
The rest of the manuscript is organised as follows. A section on preliminaries, containing the necessary notational conventions, basic definitions and results, the finite element spaces, the lifting operator, the definition of the stabilising form, and properties of the discrete Lipschitz truncation method that we use, are presented following this Introduction. An important ingredient enabling the use of the discrete Lipschitz truncation technique is a discrete infsup condition that is given in the Appendix. The finite element method is presented in Section 3, where we also show a uniform boundedness result for the sequence of approximations. Based on this and results pertaining to the discrete Lipschitz truncation, in Section 4 the convergence of the discrete solution to a weak solution of the model problem is proved using a compactness argument. Finally, some conclusions are drawn and potential future extensions are indicated.
2. Preliminaries
2.1. Notation and the problem of interest
We use standard notation for Sobolev spaces. In particular, for , and , we denote by the closure of with respect to the norm, and by the space of functions in with zero mean value. The norm in is denoted by ; when we shall use the simpler notation , and the inner product in will be denoted by . For , the norm (seminorm) in is denoted by (). Moreover, the space is the dual of with duality pairing denoted by . We also denote by the space of functions in whose distributional divergence belongs to , and by the set of elements in whose normal trace on
is zero. In the above inner products and norms we do not make a distinction between scalar and vector or tensorvalued functions.
Let be an open, bounded, polyhedral domain with a Lipschitz boundary. In this work we treat the problem with powerlaw rheology: given and a righthand side , find the velocity , the pressure , and the shearstress tensor satisfying
(2.1) 
There are many possible choices for the constitutive law, linking and the velocity . In this work we have chosen the powerlaw description where , where is a reference viscosity. In order to simplify matters we will suppose that , but we should keep in mind that, to maintain physical consistency this reference value should be kept. Similarly, in physically realistic models the gradient of the velocity is usually replaced by the symmetric velocity gradient . The results obtained in this paper can be extended, with minor modifications based on Korn’s inequality, to that case as well, so for the sake of simplicity of the exposition we shall proceed with the constitutive relation (with ) instead of .
In order to state the weak formulation of (2.1) we need to present a few additional ingredients associated with the exponent in the constitutive law relating and . For , let be its conjugate given by the relation , and let us define the critical exponent as follows:
(2.2) 
Remark 1.
With the definition (2.2) of , the space is continuously embedded in if and in , for every , if (see, e.g., [9, Corollary 9.14]). Then, in particular, is continuously embedded in and there exists a such that
(2.3) 
Moreover, the value of exhibits two different regimes, as can be seen in Figure 1, where its range of values is depicted. We will distinguish between and . The latter case occurs for and the maximum value of is attained when , at which point we have the following values:
(2.4) 
With this choice of stress tensor , the weak formulation of (2.1) is as follows: find and such that
(2.5)  
(2.6) 
2.2. Finite element spaces and preliminary results
Let be a shaperegular family of triangulations of consisting of closed simplices of diameter . To avoid technical difficulties we will suppose that the family of triangulations is quasiuniform. For reasons that will become apparent later, in the proof of convergence of the finite element method we will distinguish between the cases and . To cover the latter case (and for that purpose only) we need to make the following assumption on the mesh:
Assumption (A1). The triangulation is the result of performing one (for ), or two (for ), red refinement(s) of a, coarser, shaperegular triangulation .
We will denote the (closed) elements contained in (referred to, in some instances, as macroelements) by .
Remark 3.
raises the question whether the space itself is stable on carefully constructed meshes. Some result are known in this direction. For example, in two space dimensions, this pair is infsup stable on Powell–Sabin meshes [37]. In the recent work [23] local infsup stability is proved for this element in barycentricallyrefined meshes (also known as the Alfeld split [28], and a Hsieh–Clough–Tocher triangulation, see the references quoted in [37, P. 461]). This is then used to build enriched elements that are divergencefree (although the velocity space contains quadratic face bubbles with constant divergence). For threedimensional meshes, for the Alfeld split the lowest order infsup stable pair is (cf. [36]), while for the Powell–Sabin split the lowest order infsup stable pair is [38]. However, for the case considered in this paper, that is, taking the pair on general shaperegular meshes, stabilisation is a necessity. In addition, it is important to note that the papers cited above concern the Newtonian case only and are mostly focused on the Stokes equations. The analysis of some of those alternatives in the case of nonNewtonian flow models treated in the present work has not been carried out so far, and it will constitute a topic of future research.
Remark 4.
(i) By letting , clearly, , where does not depend on . In fact, for and for .
(ii) Under , for every , a facet of , there exists at least one node of that belongs to the interior of . In fact, this last remark is the main reason why has been made on the meshes. In particular, could also result from first making a barycentric refinement of each facet of and then building a conforming triangulation of . For ease of exposition we shall simply adopt in what follows.
In the triangulation we shall use the following notation:

: the set of all facets (edges in and faces in ) of the triangulation , with diameter . The set of internal facets is denoted by and those on the boundary of are denoted by , so ;

for every we denote by the set of facets of whose interior lies in the interior of ;

for and we define the neighbourhoods
(2.8) 
for each facet and every piecewise regular function , we denote by the jump of across ;

for we denote by the space of polynomials defined on of total degree smaller than, or equal to, , and introduce the following finite element spaces:
(2.9) (2.10) (2.11)
Using the finite element spaces defined in (2.9)–(2.11), we denote by
the Scott–Zhang interpolation operator and by
, the projections defined by (see, e.g., [17]):(2.12)  
(2.13) 
These operators satisfy ([17]):
(2.14)  
(2.15) 
The following result, whose proof can be carried out using the techniques presented in [17, Lemma 2.23], will be fundamental in the derivation (and analysis) of the proposed finite element method: for every there exists a constant , independent of , such that
(2.16) 
for all , all , and all .
We now recall three inequalities that will be useful in what follows. Let , and . The following local trace inequality is a corollary of the multiplicative trace inequality proved in [17, Lemma 12.15]:
(2.17) 
In addition, we recall the following local inverse inequality (see, e.g., [17, Lemma 12.1]): for all and all , there exists a constant , independent of , such that
(2.18) 
for every polynomial function defined on . A global version of this inequality can also be derived using the quasiuniformity of the mesh family. Finally, for , a set of indices , and any vector , the following inequality holds (see [14, Proposition 3.4(a)] for its proof):
(2.19) 
Finally, we note that under the spaces and satisfy the following discrete infsup condition: for any there exists a constant , independent of , such that for all the following inequality holds:
(2.20) 
The proof of this result, to the best of our knowledge, has not been given previously and thus we report it in the Appendix. It is based on the construction of a Fortin operator satisfying
(2.21)  
(2.22) 
In addition, (2.20) guarantees the existence of a nontrivial subspace of discretely divergencefree functions
(2.23) 
2.3. Results linked to the discrete Lipschitz truncation
In the convergence proof given below we will need the following two results. These are known as discrete Lipschitz truncation and divergencefree discrete Lipschitz truncation, respectively. Their proofs are omitted since they are essentially a rewriting of Corollary 17 and the proof on pages 1006–1007 in [15] (see also [35, Lemmas 2.29 and 2.30]).
Lemma 5.
Let . Let us suppose that for all and weakly in as . Then, there exist

a double sequence such that for all ;

a double sequence of open sets , of the form
(2.24) where denotes the collection of some elements of the mesh ;

a double sequence with for all and all ;
satisfying

in for all and all ;

there exists a such that
(2.25) 
there exists a such that
(2.26) 
for any fixed ,
(2.27) as .
Lemma 6.
Let and assume that Assumption (A1) is satisfied. Let be a sequence such that for all and such that weakly in as . Furthermore, let be the sequence of Lipschitz truncations given by Lemma 5. Then, there exists a double sequence such that

for all and all ;

there exists a such that
(2.28) 
for any fixed the following convergences hold (up to a subsequence, if necessary):
(2.29) as , for all .
2.4. The stabilising bilinear form and the lifting operator
The finite element method studied in this work is based on the pair . Since this pair is not infsup stable some form of stabilisation is needed. In this work our proposal is to use the following stabilising bilinear form
(2.30) 
where the stabilisation parameter is defined as follows:
(2.31) 
The behaviour of is depicted in Figure 1. It can be observed there that the stabilisation gets stronger as . The reason for this behaviour will become clear when we perform the convergence analysis in Section 4.
Remark 7.
Thanks to and the infsup condition (2.20) it can be expected to have stability of (where is the finite element approximation of the pressure). The stabilisation is then built with the aim of controlling . More precisely, using (2.16), (2.19) (or the inverse inequality (2.18)), and the definition of the bilinear form we see that there exists a constant such that
(2.32) 
for all , where
(2.33) 
It will be useful in what follows to observe that for we have .
Another important ingredient in the definition of the method is a lifting of the pressure jumps defined with the help of the lowest order Raviart–Thomas basis functions. To define this, for each we choose a unique normal vector . Its orientation is of no importance, but it needs to point outwards of if . Moreover, for each such that , we denote the node in opposite by . Using this unique normal vector, we introduce the lowest order Raviart–Thomas basis function defined as
(2.34) 
and extended by zero outside . In this definition, the sign of the function depends on whether the normal vector points in or out of . Thanks to its definition, satisfies the following: for every the normal component of is given by (with the obvious abuse of notation considering that is not defined at the boundary of ):
(2.35) 
With the help of these Raviart–Thomas basis functions, we define the following operator, which will be fundamental in the definition of the finite element method:
(2.36) 
Since the velocity is bounded in , then it is bounded in as well. In the finite element method proposed in Section 3, we will consider a modified velocity built with the help of the mapping just defined. The following result states that the stability just mentioned is preserved by the operator .
Lemma 8.
There exists a constant , independent of , such that
(2.37) 
for all .
Proof.
Thanks to the embedding (2.3) and denoting
(2.38) 
the following bound follows
(2.39) 
To bound the second term on the righthand side of this inequality we start by noticing that the definition of (cf. (2.34)) gives for each such that . So, let and let be the unique macroelement such that . Then, using the mesh regularity and the CauchySchwarz inequality we get
Hence, squaring, summing over all the elements, and using the mesh regularity gives
(2.40) 
Thus, using the inverse inequality (2.18) we arrive at
(2.41) 
To complete the proof we only need to make sure that the exponent of in (2.41) is not negative. Let . If then and . So, . If then and so
Since in the whole range of values for we have , the proof is complete. ∎
3. The finite element method
The finite element method studied in this work reads as follows: find such that
(3.1)  
(3.2) 
for all , where is defined by (2.36) and the stabilising bilinear form is defined in (2.30).
Remark 9.
(i) The main differences between (3.1), (3.2) and a standard Galerkin method are twofold: first, the stabilising term involving the jumps of the discrete pressure are added to the formulation to compensate for the fact that the pair does not satisfy the discrete infsup condition. Additionally, and perhaps more significantly, the convection velocity has been replaced by the modified version . In Lemma 10 this modified velocity will be proved to be solenoidal, which allows us to analyse the finite element method without the need to rewrite the convection term in a skewsymmetric form. This will lead to a convergence result valid in the whole range .
(ii) As can be expected, the power of in the stabilisation parameter depends strongly on the value of . Two important remarks are in order:

for all ;

for all we have .
Thus, there is always a positive power of multiplying the jump terms of the pressure involved in the definition of and , but the stabilisation becomes stronger as .
3.1. Existence of a solution and a priori bounds
Before exploring the stability of the scheme, we present the following a priori result concerning qualitative properties of and , whenever solves (3.1), (3.2).
Lemma 10.

is discretely divergencefree with respect to the coarse space , that is,
(3.3) 
on , and
(3.4)
Proof.
The proof of (i) is a consequence of the fact that the stabilisation vanishes on the coarse space , that is, for all and all . For (ii), we can follow similar arguments as those presented in [6, Lemma 3.8] and [7, Lemma 3] (see also [2, Theorem 3] for a different proof).
The following result states the existence of a solution to the discrete problem (3.1), (3.2). In addition, it provides uniform a priori bounds for the sequence of solutions as .