The construction and mathematical analysis of finite element approximations of models of non-Newtonian 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.,  and ). Motivated by the groundbreaking contributions of Cohen, Dahmen and DeVore [12, 13] and Binev, Dahmen and DeVore  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 non-Newtonian fluid flow models (cf. ).
Upon decomposing the Cauchy stress tensor into its traceless part, called the deviatoric stress tensor or shear-stress 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 shear-stress tensor, . For Newtonian fluids the shear-stress tensor is a scalar multiple of the symmetric velocity gradient. The finite element approximation of Newtonian fluids is therefore usually performed in the velocity-pressure formulation. For non-Newtonian fluids on the other hand the situation is more involved, because the shear-stress tensor exhibits nonlinear dependence as a function of the symmetric velocity gradient, and the functional relationship between the shear-stress tensor and the symmetric velocity gradient may even be completely implicit and multi-valued. For power-law fluids, such as the ones considered in this work, the shear-stress tensor exhibits power-law 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 power-law 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 , 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 power-law rheology have been extensively studied, including stabilised (or variational-multiscale) methods (cf. [10, 1], for example) and local discontinuous Galerkin methods (see, , 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 , for example. Concerning implicitly-constituted models, in the recent papers [15, 33] the convergence of generic inf-sup stable velocity/pressure-based 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 divergence-free functions. The reason for this dichotomy is that in the case of velocity approximations that are discretely divergence-free only, as is the case in generic inf-sup stable mixed finite element methods, the finite element approximation of the convection term does not vanish when tested with
, and it needs to be skew-symmetrized (cf.) 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 skew-symmetric modification can be bounded by the same expression using Hölder’s inequality, this is not the case for the power-law model under consideration here for entire range for which weak solutions to the problem are known to exist. In fact, in the case of non-Newtonian power-law 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 skew-symmetric modification of the convection term, whose use is essential so as to be able to derive an energy inequality for discretely divergence-free velocity fields, this can only be achieved for the limited range . This was precisely the bottleneck encountered in  for discretely divergence-free velocity approximations, resulting in the reduction of the range of from the maximal range for which weak solutions are known to exist, to
The advantage of pointwise divergence-free finite element methods over discretely divergence-free finite element methods is therefore that, besides the physical consistency they provide, there is no need to rewrite the convection term in a skew-symmetric form. The topic of divergence-free 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  (analysed recently in ) 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 high-enough degree (for example, is needed in , 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 , 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 , for example. In this latter case the viscous term (defined as the divergence of the shear-stress 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 power-law model considered herein.
The recent works [2, 3] offer a way of preserving the advantages of a pointwise divergence-free approximation to the velocity field while working with the, computationally simplest, lowest-order -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 skew-symmetric form. Our aim here is to apply this idea to a problem in non-Newtonian fluid mechanics. As a first step in this direction, we have chosen an explicit constitutive law with power-law 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 , and the experimental results in, e.g., ); 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 skew-symmetric form, the resulting method can now be proved to be stable and convergent to a weak solution for the whole range of the power-law 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 inf-sup 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.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 tensor-valued functions.
Let be an open, bounded, polyhedral domain with a Lipschitz boundary. In this work we treat the problem with power-law rheology: given and a right-hand side , find the velocity , the pressure , and the shear-stress tensor satisfying
There are many possible choices for the constitutive law, linking and the velocity . In this work we have chosen the power-law 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:
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
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:
With this choice of stress tensor , the weak formulation of (2.1) is as follows: find and such that
In order for the variational formulation (2.5), (2.6) to be meaningful it is necessary that with , which necessitates that , and under this condition the existence of a solution to (2.5), (2.6) has been proved (see ). Thus, for the rest of this work we will assume that .
2.2. Finite element spaces and preliminary results
Let be a shape-regular family of triangulations of consisting of closed simplices of diameter . To avoid technical difficulties we will suppose that the family of triangulations is quasi-uniform. 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, shape-regular triangulation .
We will denote the (closed) elements contained in (referred to, in some instances, as macro-elements) by .
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 inf-sup stable on Powell–Sabin meshes . In the recent work  local inf-sup stability is proved for this element in barycentrically-refined meshes (also known as the Alfeld split , and a Hsieh–Clough–Tocher triangulation, see the references quoted in [37, P. 461]). This is then used to build enriched elements that are divergence-free (although the velocity space contains quadratic face bubbles with constant divergence). For three-dimensional meshes, for the Alfeld split the lowest order inf-sup stable pair is (cf. ), while for the Powell–Sabin split the lowest order inf-sup stable pair is . However, for the case considered in this paper, that is, taking the pair on general shape-regular 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 non-Newtonian flow models treated in the present work has not been carried out so far, and it will constitute a topic of future research.
(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
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)
the Scott–Zhang interpolation operator and by, the projections defined by (see, e.g., ):
These operators satisfy ():
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
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]:
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
for every polynomial function defined on . A global version of this inequality can also be derived using the quasi-uniformity 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):
Finally, we note that under the spaces and satisfy the following discrete inf-sup condition: for any there exists a constant , independent of , such that for all the following inequality holds:
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
In addition, (2.20) guarantees the existence of a non-trivial subspace of discretely divergence-free functions
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 divergence-free 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  (see also [35, Lemmas 2.29 and 2.30]).
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
where denotes the collection of some elements of the mesh ;
a double sequence with for all and all ;
in for all and all ;
there exists a such that
there exists a such that
for any fixed ,
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
for any fixed the following convergences hold (up to a subsequence, if necessary):
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 inf-sup stable some form of stabilisation is needed. In this work our proposal is to use the following stabilising bilinear form
where the stabilisation parameter is defined as follows:
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.
Thanks to and the inf-sup 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
for all , where
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
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 ):
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:
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 .
There exists a constant , independent of , such that
for all .
Thanks to the embedding (2.3) and denoting
the following bound follows
To bound the second term on the right-hand 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 macro-element such that . Then, using the mesh regularity and the Cauchy-Schwarz inequality we get
Hence, squaring, summing over all the elements, and using the mesh regularity gives
Thus, using the inverse inequality (2.18) we arrive at
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
(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 inf-sup 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 skew-symmetric 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
is discretely divergence-free with respect to the coarse space , that is,
on , and
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).