# Virtual Element methods for non-Newtonian shear-thickening fluid flow problems

Paola F. Antonietti, Lourenço Beirão da Veiga, Michele Botti, André Harnist, Giuseppe Vacca

## ► To cite this version:

Paola F. Antonietti, Lourenço Beirão da Veiga, Michele Botti, André Harnist, Giuseppe Vacca. Virtual Element methods for non-Newtonian shear-thickening fluid flow problems. 2026. hal-05448476

HAL Id: hal-05448476

https://hal.science/hal-05448476v1

Preprint submitted on 8 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
Virtual Element methods for non-Newtonian shear-thickening fluid flow problems

P. F. Antonietti*¹, L. Beirão da Veiga†²,⁵, M. Botti‡¹, A. Harnist§³, G. Vacca¶⁴, and M. Verani ||¹

¹MOX-Laboratory for Modeling and Scientific Computing, Dipartimento di Matematica, Politecnico di Milano, Piazza Leonardo da Vinci 32 - 20133 Milano, Italy

²Dipartimento di Matematica e Applicazioni, Università degli Studi di Milano-Bicocca, Via Roberto Cozzi 55 - 20125 Milano, Italy
³Université de technologie de Compiègne, Alliance Sorbonne Université, LMAC, Compiègne, France

⁴Dipartimento di Matematica, Università degli Studi di Bari, Via Edoardo Orabona 4 - 70125 Bari, Italy

⁵IMATI-CNR, Via Ferrata 1 - 27100 Pavia, Italy

January 8, 2026

Abstract

In this work, we present a comprehensive theoretical analysis for Virtual Element discretizations of incompressible non-Newtonian flows governed by the Carreau-Yasuda constitutive law, in the shear-thickening regime ($r > 2$) including both degenerate ($\delta = 0$) and non-degenerate ($\delta > 0$) cases. The proposed Virtual Element method features two distinguishing advantages: the construction of an exactly divergence-free discrete velocity field and compatibility with general polygonal meshes. The analysis presented in this work extends the results of [5], where only shear-thinning behavior ($1 < r < 2$) was considered. Indeed, the theoretical analysis of the shear-thickening setting requires several novel analytical tools, including: an inf-sup stability analysis of the discrete velocity-pressure coupling in non-Hilbertian norms, a stabilization term specifically designed to address the nonlinear structure as the exponent $r > 2$; and the introduction of a suitable discrete norm tailored to the underlying nonlinear constitutive relation. Numerical results demonstrate the practical performance of the proposed formulation.

1 Introduction

Numerous applications, such as polymer processing, additive manufacturing, material deposition, concentrated suspensions, and high-shear biological fluids, as well as various materials science problems, involve fluids that exhibit non-Newtonian behavior. A nonlinear relation between the strain rate and the shear stress characterizes this behavior. A paradigmatic example is the Carreau-Yasuda model, introduced in [58], where the relation between the shear stress $\sigma(\cdot, \epsilon)$ and the strain rate $\epsilon$ is given by

$$ \sigma(\cdot, \epsilon) = \mu(\cdot)(\delta^{\alpha} + |\epsilon|^{\alpha})^{\frac{r-2}{\alpha}} \epsilon, $$

for $\mu$, $\alpha$, $\delta$ and the power-law index $r$ to be specified later on. The Carreau-Yasuda model is a generalization of the Carreau model, which corresponds to the choice $\alpha = 2$. The Carreau-Yasuda model models both shear-thinning (pseudo-plastic) behavior for $1 < r < 2$ and shear-thickening (dilatant) behavior for $r > 2$. When $r = 2$, the model reduces to the standard Newtonian fluid case. The case $\delta = 0$ corresponds to the classical power-law model. From a mathematical and numerical analysis perspective, the shear-thickening regime presents distinct challenges, because the value of the power-law index ($r > 2$) and the possible

<footnote>*paola.antonietti@polimi.it</footnote>

<footnote>† lourenco.beirao@unimib.it</footnote>

<footnote>‡michele.botti@polimi.it</footnote>

<footnote>§andre.harnist@utc.fr</footnote>

<footnote>¶giuseppe.vacca@uniba.it</footnote>

<footnote>||marco.verani@polimi.it</footnote>
presence of the degenerate case $\delta = 0$ lead to a significantly different mathematical structure. Indeed, in such cases, designing robust discretization schemes and deriving stability and a priori estimates are particularly challenging.

Numerical methods for non-Newtonian flows have a long history, beginning with the seminal work of [8], which proposed a Finite Element approximation of a non-Newtonian flow model governed by either the Carreau or the power-law model. Still in the Finite Element framework, sharp error estimates were subsequently established in [53, 10, 9]. In particular, the pioneering studies [10] and [9] derived (in some cases, optimal) velocity and pressure error bounds in appropriate quasi-norms for models incorporating Carreau or power-law models. We also refer to, e.g., [20, 43, 47, 45] for more recent contributions. In practical applications, computational domains often exhibit complex geometric features, necessitating discretization schemes that efficiently accommodate flexible, possibly adapted, grids. Consequently, discretization methods that can support general polyhedral meshes have been recently investigated, including Virtual Element Methods (VEM) for non-Newtonian incompressible fluids in the shear-thinning regime ($1 < r < 2$), see [5], discontinuous Galerkin [51] and hybridizable discontinuous Galerkin [41] methods, and Hybrid High-Order schemes for non-Newtonian fluids governed by the Stokes equations [25] and Navier-Stokes equations with nonlinear convection in [30].

In this work, we focus on the Virtual Element method, originally introduced in [19] for second-order elliptic problems and subsequently generalized to a broad class of differential problems. In the context of fluid flow problems, a key advantage of VEM is its ability to construct divergence-free discrete velocity spaces on general polygonal or polyhedral meshes, thereby eliminating the need for ad hoc stabilization and ensuring mass conservation. For this reason, over the past decade, VEM has been extensively developed for approximating Newtonian fluid flow problems. In [7], a novel stream-function-based VEM formulation for the Stokes problem, relying on a suitably designed stream function space that characterizes the divergence-free subspace of the discrete velocity field, has been proposed and analyzed. Divergence-free virtual elements have been introduced in [15]. Further developments for the Stokes problem can be found in [29, 27, 13, 31, 48, 23, 14, 22]; see also [37, 12], where arbitrary-order pressure-robust VEMs have been investigated. The Virtual Element discretization of the Navier-Stokes equations was first studied in [16] and further investigated in [39, 18, 1, 40]; we also refer to [28, 2, 50], for quasi-Newtonian Stokes flows. Recently, Virtual Element approximations to the coupled Navier-Stokes and heat equations have been analyzed in [4, 17]. Moreover, least-squares Virtual Element discretizations of the Stokes and Navier-Stokes systems have recently been analyzed in [49] and [57], respectively. For a comprehensive review of recent advances in VEM, we refer to the monograph [6].

In this paper, we analyze a Virtual Element discretization for steady incompressible non-Newtonian flow governed by the Carreau-Yasuda constitutive law, explicitly addressing the shear-thickening range ($r > 2$) and both the degenerate ($\delta = 0$) and non-degenerate ($\delta > 0$) cases. To the best of the authors' knowledge, previous $W^{1,r}$ a priori estimates on polytopal methods for non-Newtonian flow have focused on the shear-thinning (pseudoplastic) regime, with the exception of the results in [25], where a Hybrid High Order discretization of non-Newtonian fluids with $r \in (1, \infty)$ is studied, leading to velocity and pressure error bounds of order $\frac{k+1}{r-1}$, being $k$ the polynomial order of the discretization. Notice that here we get a similar dependence of $r$ in the convergence order. However, we provide a better convergence order ($2k/r$) of the $W^{1,r}$ error norm under more regularity of the flux contribution. The Virtual Element formulation is based on employing the divergence-free Virtual Element spaces introduced in [15, 16] and we focus on the extension of the stability and convergence analysis in the case $r > 2$ and, possibly, $\delta = 0$, addressing thereby the theoretical analysis not covered in [5], where only the shear thinning regimes was considered. The proposed method offers two principal advantages: it accommodates general polygonal meshes and leads to a discrete velocity field that is exactly divergence-free. We point out that, the mathematical analysis of the shear-thickening range $r > 2$ and the degenerate case $\delta = 0$ presented in this work requires several novel technical contributions, including the establishment of inf-sup stability of the discrete velocity-pressure coupling in non-Hilbertian norms and the introduction of a stabilization term tailored to the distinct analytical structure that arises when the exponent crosses the threshold $r = 2$. Furthermore, we note that, in the presence of non convex-elements $E$, the definition of the aforementioned Virtual Elements does not guarantee that the local spaces are in the natural Sobolev space $W^{1,r}(E)$ for large values of $r$. As a consequence, and in order to include in our analysis all type of polygonal meshes, we evaluate the stability and the error of the scheme in a discrete norm. Whenever the elements are convex, the standard norm is immediately recovered. We prove that if $r > 2$, and the underlying solution is sufficiently smooth, the order of convergence of the velocity and the pressure is
$\frac{k}{r-1}$ for $\delta \ge 0$. In the non-degenerate case, i.e. $\delta > 0$, we show that the order of convergence is $\frac{2k}{r}$. To the authors' best knowledge, this is the first result for a polytopal discretization of non-Newtonian flows in the degenerate power-law/shear-thickening regime.

The remainder of the manuscript is organized as follows. In Section 2, after introducing the notation used throughout the paper, we present the weak formulation of the continuous problem and discuss its well-posedness. Section 3 describes the proposed divergence-free Virtual Element discretization and presents the well-posedness analysis. The *a priori* error analysis is detailed in Section 4. Section 5 reports numerical experiments illustrating the method's performance. Finally, Section 6 contains a summary of the obtained results and some concluding remarks.

## 2 Model problem

In this section, we introduce the notation used throughout the paper, present the model problem, and discuss its well-posedness.

The vector spaces considered hereafter are over $\mathbb{R}$. We denote by $\mathbb{R}_+$ the set of non-negative real numbers. Given a vector space $V$ with norm $\| \cdot \|_V$, the notation $V' \ denotes its dual space and $V' \langle \cdot, \cdot \rangle_V$ the duality between $V$ and $V'$. The notation $v \cdot w$ and $v \times w$ designates the scalar and vector products of two vectors $v, w \in \mathbb{R}^d$, and $|v|$ denotes the Euclidean norm of $v$ in $\mathbb{R}^d$. The inner product in $\mathbb{R}^{d \times d}$ is defined for $\tau, \eta \in \mathbb{R}^{d \times d}$ by $\tau : \eta := \sum_{i,j=1}^d \tau_{i,j} \eta_{i,j}$ and the induced norm is given by $|\tau| = \sqrt{\tau : \tau}$.

Let $\Omega \subset \mathbb{R}^d$ denote a bounded, connected, polyhedral open set with Lipschitz boundary $\partial \Omega$ and let $n$ be the outward unit normal to $\partial \Omega$. To simplify the exposition, we restrict the presentation to the two-dimensional case, i.e., $d=2$, but the analysis remains valid in the three-dimensional case $d=3$ as well, with minor technical differences. We denote with $x = (x_1, x_2)$ the independent variable. We assume that the boundary is partitioned into two disjoint subsets $\partial \Omega = \Gamma_D \cup \Gamma_N$, with $|\Gamma_D| > 0$, such that a Dirichlet condition is given on $\Gamma_D$ and a Neumann condition on $\Gamma_N$.

Throughout the article, spaces of functions, vector fields, and tensor fields, defined over any $X \subset \bar{\Omega}$ are denoted by italic capitals, boldface Roman capitals, and special Roman capitals, respectively. The subscript $s$ denotes a space of symmetric tensor fields. For example, $L^2(X)$, $L^2(X)$, and $L_s^2(X)$ denote the spaces of square-integrable functions, vector fields, and symmetric tensor fields, respectively. The notation $W^{m,r}(X)$, for $m \ge 0$ and $r \in [1, +\infty]$, with the convention that $W^{0,r}(X) = L^r(X)$, and $W^{m,2}(X) = H^m(X)$, designates the classical Sobolev spaces. The trace map is denoted by $\gamma : W^{1,r}(\Omega) \to W^{1-\frac{1}{r},r}(\partial \Omega)$. Finally, given $\Gamma \subset \partial \Omega$, we denote by $W_{0,\Gamma}^{1,r}(\Omega)$ the subspace of $W^{1,r}(\Omega)$ spanned by functions having zero trace on $\Gamma$. The symbol $\nabla$ denotes the gradient for scalar functions, while $\nabla \cdot \epsilon(\cdot) = \frac{\nabla(\cdot)+\nabla^T(\cdot)}{2}$, and $\nabla \cdot$ denote the gradient, the symmetric gradient operator, and the divergence operator, respectively, whereas $\nabla \cdot$ denotes the vector-valued divergence operator for tensor fields.

We consider the creeping flow of a non-Newtonian fluid occupying $\Omega$ and subjected to a volumetric force $f : \Omega \to \mathbb{R}^d$ and a traction $g : \Gamma_N \to \mathbb{R}^d$ described by the non-linear Stokes equations

$$ 
\begin{aligned}
-\nabla \cdot \sigma(\cdot, \epsilon(\mathbf{u})) + \nabla p &= f && \text{in } \Omega, \\
\nabla \cdot \mathbf{u} &= 0 && \text{in } \Omega, \\
\sigma(\cdot, \epsilon(\mathbf{u}))\mathbf{n} - p\mathbf{n} &= \mathbf{g} && \text{on } \Gamma_N, \\
\mathbf{u} &= \mathbf{0} && \text{on } \Gamma_D,
\end{aligned}
\quad (1) $$

where $\mathbf{u} : \Omega \to \mathbb{R}^d$ and $p : \Omega \to \mathbb{R}$ denote the velocity field and the pressure field, respectively. Non-homogeneous Dirichlet conditions can be considered in place of (1) up to minor modifications.

In this work, we consider the Carreau-Yasuda model, introduced in [58], as a reference model for the non-linear shear stress-strain rate relation, i.e.,

$$ \sigma(\mathbf{x}, \epsilon(\mathbf{v})) = \mu(\mathbf{x})(\delta^{\alpha} + |\epsilon(\mathbf{v})|^{\alpha})^{\frac{r-2}{\alpha}} \epsilon(\mathbf{v}), \quad (2) $$

where $\mu : \Omega \to [\mu_-, \mu_+]$, with $0 < \mu_- < \mu_+ < \infty$, $\alpha \in (0, \infty)$, and $\delta \ge 0$ and $r \in [2, \infty)$. The Carreau-Yasuda law is a generalization of the Carreau model corresponding to the case $\alpha = 2$. The case $\delta = 0$ corresponds to the classical power-law model. The stress-strain law (2) satisfies the following assumption:
**Assumption 1.** The shear stress-strain rate law $\sigma : \Omega \times \mathbb{R}_s^{d \times d} \to \mathbb{R}_s^{d \times d}$ appearing in (1) is a Caratheodory function satisfying $\sigma(\cdot, \mathbf{0}) = \mathbf{0}$ and for a fixed $r \in [2, +\infty)$ there exist real numbers $\delta \in [0, +\infty)$ and $\sigma_c, \sigma_m \in (0, +\infty)$ such that the following conditions hold:

$$ |\sigma(\mathbf{x}, \tau) - \sigma(\mathbf{x}, \eta)| \le \sigma_c (\delta^r + |\tau|^r + |\eta|^r)^{\frac{r-2}{r}} |\tau - \eta|, \quad (\text{Hölder continuity}) \quad (3a) $$

$$ (\sigma(\mathbf{x}, \tau) - \sigma(\mathbf{x}, \eta)) : (\tau - \eta) \ge \sigma_m (\delta^r + |\tau|^r + |\eta|^r)^{\frac{r-2}{r}} |\tau - \eta|^2, \quad (\text{strong monotonicity}) \quad (3b) $$

for almost every $\mathbf{x} \in \Omega$ and all $\tau, \eta \in \mathbb{R}_s^{d \times d}$.

The constants $\sigma_c, \sigma_m > 0$ in (3) for the Carreau-Yasuda model (2) satisfy

$$ \sigma_c = \mu_+(r-1)2^{\left(\frac{1}{a}-\frac{1}{r}\right)^{\Theta}(r-2)} \quad \text{and} \quad \sigma_m = \frac{\mu_-}{r-1}2^{\left[-\left(\frac{1}{a}-\frac{1}{r}\right)^{\Theta}-1\right](r-2)-1}, $$

where $\xi^\oplus := \max(0, \xi)$ and $\xi^\ominus := -\min(0, \xi)$ for all $\xi \in \mathbb{R}$ (cf. [25]). Finally, for further use, we adopt the short-hand notation $a \le b$ to denote the inequality $a \le Cb$, for a constant $C > 0$ that may depend on $\sigma_c, \sigma_m, \delta$ (or related parameters) in Assumption 1 and $r$, but is independent of the discretization parameters. The obvious extensions $a \ge b$ and $a \approx b$ hold.

## 2.1 Weak formulation

In this section, we provide the variational formulation of (1). For $r \in (1, \infty)$, we introduce the conjugate index defined as $r' = \frac{r}{r-1}$ and recall Korn's first inequality (see, e.g., [32, Theorem 1.2] and [42, Theorem 1]): there is $C_K > 0$ depending only on $\Omega$ and $r$ such that for all $\mathbf{v} \in W_{0, \Gamma_D}^{1, r}(\Omega)$,

$$ \|\mathbf{v}\|_{W^{1,r}(\Omega)} \le C_K \|\epsilon(\mathbf{v})\|_{L^r(\Omega)}. \quad (4) $$

Let $r \in [2, \infty)$ be the Sobolev exponent dictated by the non-linear stress-strain law characterizing problem (1) and satisfying Assumption 1. We define the velocity and pressure spaces incorporating the boundary condition on $\Gamma_D$ and the zero-average constraint in the case $\Gamma_D = \partial\Omega$, respectively:

$$ U = W_{0, \Gamma_D}^{1, r}(\Omega) \quad P = \begin{cases} L^{r'}(\Omega) & \text{if } |\Gamma_D| < |\partial\Omega| \\ L_0^{r'}(\Omega) := \{q \in L^{r'}(\Omega) : \int_{\Omega} q = 0\} & \text{if } \Gamma_D = \partial\Omega. \end{cases} $$

Assuming $f \in L^{r'}(\Omega)$ and $g \in L^{r'}(\Gamma_N)$, the weak form of (1) reads: Find $(u, p) \in U \times P$ such that

$$ \begin{aligned} a(\mathbf{u}, \mathbf{v}) + b(\mathbf{v}, p) &= \int_{\Omega} f \cdot \mathbf{v} + \int_{\Gamma_N} g \cdot \mathbf{v} & \quad \forall \mathbf{v} \in U, \\ -b(\mathbf{u}, q) &= 0 & \quad \forall q \in P, \end{aligned} \quad (5) $$

where $a: U \times U \to \mathbb{R}$ and $b: U \times P \to \mathbb{R}$ are defined for all $\mathbf{v}, \mathbf{w} \in U$ and all $q \in L^{r'}(\Omega)$ by

$$ a(\mathbf{w}, \mathbf{v}) := \int_{\Omega} \sigma(\cdot, \epsilon(\mathbf{w})) : \epsilon(\mathbf{v}), \quad b(\mathbf{v}, q) := - \int_{\Omega} (\nabla \cdot \mathbf{v}) q. \quad (6) $$

Introducing the subset $Z = \{\mathbf{v} \in U \text{ s.t. } \nabla \cdot \mathbf{v} = 0\} \subset U$ corresponding to the kernel of $b(\cdot, \cdot)$, problem (5) can be formulated in the equivalent kernel form: Find $\mathbf{u} \in Z$ such that

$$ a(\mathbf{u}, \mathbf{v}) = \int_{\Omega} f \cdot \mathbf{v} + \int_{\Gamma_N} g \cdot \mathbf{v} \quad \forall \mathbf{v} \in Z. \quad (7) $$

## 2.2 Well-posedness

In this section, we report the properties of the functions $a(\cdot, \cdot)$ and $b(\cdot, \cdot)$ defined in (6) and we prove the well-posedness of (5).

**Lemma 1** (Continuity and monotonicity of $a$). For all $\mathbf{u}, \mathbf{v}, \mathbf{w} \in U$, setting $\mathbf{e} = \mathbf{u} - \mathbf{w}$, it holds

$$ |a(\mathbf{u}, \mathbf{v}) - a(\mathbf{w}, \mathbf{v})| \le \sigma_c \left( \delta^r + \|\mathbf{u}\|_{W^{1,r}(\Omega)}^r + \|\mathbf{w}\|_{W^{1,r}(\Omega)}^r \right)^{\frac{r-2}{r}} \|\mathbf{e}\|_{W^{1,r}(\Omega)} \|\mathbf{v}\|_{W^{1,r}(\Omega)}, \quad (8a) $$

$$ a(\mathbf{u}, \mathbf{e}) - a(\mathbf{w}, \mathbf{e}) \ge \sigma_m \|\mathbf{e}\|_{W^{1,r}(\Omega)}^r. \quad (8b) $$
**Proof.** We follow the lines of [25, Lemma 7.3]. Let $u, v, w \in U$ and set $e = u - w$.

(i) **Hölder continuity.** Recalling the Hölder continuity property (3a) and using the generalized Hölder inequality with exponents $(\frac{r}{r-2}, r, r)$, we have

$$ egin{aligned} |a(u, v) - a(w, v)| &\le \int_{\Omega} |[\sigma(\cdot, \epsilon(u)) - \sigma(\cdot, \epsilon(w))] : \epsilon(v)| \\ &\le \sigma_c \left( \int_{\Omega} \delta^r + |\epsilon(u)|^r + |\epsilon(w)|^r \right)^{\frac{r-2}{r}} \|e\|_{W^{1,r}(\Omega)} \|v\|_{W^{1,r}(\Omega)}. \end{aligned} \quad (9) $$

(ii) **Strong monotonicity.** First, we observe that for all $\tau, \eta \in \mathbb{R}_s^{d \times d}$ the triangle inequality implies $2^{1-r}|\tau - \eta|^r \le |\tau|^r + |\eta|^r$. Thus, since $\delta \ge 0$ and $r \ge 2$,

$$ (\delta^r + |\tau|^r + |\eta|^r)^{\frac{r-2}{r}} \ge (2^{1-r} |\tau - \eta|^r)^{\frac{r-2}{r}} \ge |\tau - \eta|^{r-2}. \quad (10) $$

Using Korn's inequality (4) together with (3b) and the inequality in (10), we obtain

$$ egin{aligned} \sigma_m \|e\|_{W^{1,r}(\Omega)}^r &\le C_K^r \left( \int_{\Omega} \sigma_m |\epsilon(u-w)|^r \right) \le \int_{\Omega} \sigma_m |\epsilon(u) - \epsilon(w)|^2 |\epsilon(u) - \epsilon(w)|^{r-2} \\ &\le \int_{\Omega} (\delta^r + |\epsilon(u)|^r + |\epsilon(w)|^r)^{\frac{2-r}{r} + \frac{r-2}{r}} (\sigma(\cdot, \epsilon(u)) - \sigma(\cdot, \epsilon(w))) : \epsilon(u-w) \\ &\le a(u, e) - a(w, e). \end{aligned} $$

□

The following result is needed to infer the existence of a unique pressure $p \in P$ solving problem (5) from the well-posedness of problem (7). For its proof, we refer to [24, Theorem 1].

**Lemma 2 (Inf-sup condition).** For any $r \in [2, \infty)$ there exists a positive constant $\beta(r)$ such that the bilinear form $b(\cdot, \cdot)$ defined in (6) satisfies

$$ \inf_{q \in P} \sup_{w \in U \setminus \{0\}} \frac{b(w, q)}{\|q\|_{L^{r'}(\Omega)} \|w\|_{W^{1,r}(\Omega)}} \ge \beta(r) > 0. \quad (11) $$

We are now ready to establish the well-posedness of problem (5).

**Proposition 3 (Well-posedness).** For any $r \in [2, \infty)$, there exists a unique solution $(u, p) \in U \times P$ to problem (5) satisfying the a priori estimates

$$ \|u\|_{W^{1,r}(\Omega)} \le \sigma_m^{-\frac{1}{r-1}} \left( \|f\|_{L^{r'}(\Omega)} + \|g\|_{L^{r'}(\Gamma_N)} \right)^{\frac{1}{r-1}}, \quad (12a) $$

$$ \|p\|_{L^{r'}(\Omega)} \le \left(1 + \frac{\sigma_c}{\sigma_m}\right) \left( \|f\|_{L^{r'}(\Omega)} + \|g\|_{L^{r'}(\Gamma_N)} \right) + \sigma_c \delta^{r-2} \left( \frac{\|f\|_{L^{r'}(\Omega)} + \|g\|_{L^{r'}(\Gamma_N)}}{\sigma_m} \right)^{\frac{1}{r-1}}. \quad (12b) $$

**Proof.** We focus here on the a priori bounds (12a) and (12b), while for the uniqueness and existence we refer to [11, 21]. Owing to (8b) with $w = 0$, taking $v = u$ in (7), and using the Hölder inequality together with the continuity of the trace map, one has

$$ \sigma_m \|u\|_{W^{1,r}(\Omega)}^r \le a(u, u) = \int_{\Omega} f \cdot u + \int_{\Gamma_N} g \cdot \gamma(u) \le (\|f\|_{L^{r'}(\Omega)} + \|g\|_{L^{r'}(\Gamma_N)}) \|u\|_{W^{1,r}(\Omega)}. $$

From the previous bound, the velocity estimate in (12a) follows. Concerning the estimate of the pressure field: owing to the inf-sup condition (11) and equation (5), it is inferred that

$$ \beta(r) \|p\|_{L^{r'}(\Omega)} \le \sup_{v \in U \setminus \{0\}} \frac{b(v, p)}{\|v\|_{W^{1,r}(\Omega)}} = \sup_{v \in U \setminus \{0\}} \frac{\int_{\Omega} f \cdot v + \int_{\Gamma_N} g \cdot v - a(u, v)}{\|v\|_{W^{1,r}(\Omega)}}. $$

Applying the Hölder inequality, the continuity of $a$ in (8a) with $w = 0$, and the a priori estimate of the velocity (12a), we obtain

$$ egin{aligned} \|p\|_{L^{r'}(\Omega)} &\le \|f\|_{L^{r'}(\Omega)} + \|g\|_{L^{r'}(\Gamma_N)} + \sigma_c \left( \delta^r + \|u\|_{W^{1,r}(\Omega)}^r \right)^{\frac{r-2}{r}} \|u\|_{W^{1,r}(\Omega)} \\ &\le \left(1 + \frac{\sigma_c}{\sigma_m}\right) \left( \|f\|_{L^{r'}(\Omega)} + \|g\|_{L^{r'}(\Gamma_N)} \right) + \sigma_c \delta^{r-2} \sigma_m^{-\frac{1}{r-1}} \left( \|f\|_{L^{r'}(\Omega)} + \|g\|_{L^{r'}(\Gamma_N)} \right)^{\frac{1}{r-1}}. \end{aligned} $$

□
# 3 Virtual Element method

In the present section, we initially review divergence-free Virtual Elements of general order [15, 16]. Afterterwards, we design the computable forms and formulate the discrete problem that approximates the nonlinear equation (5), finally establishing its well-posedness.

## 3.1 Mesh and discrete spaces

Let {$\Omega_h$}$_h$ be a sequence of decompositions of the domain $\Omega \subset \mathbb{R}^2$ into general polytopal elements. For each $E \in \Omega_h$, we denote with $h_E$ the diameter, with $|E|$ the area, with $x_E = (x_{E,1}, x_{E,2})$ the centroid and we set $h = \sup_{E \in \Omega_h} h_E$. We suppose that {$\Omega_h$}$_h$ fulfills the following assumption.

**Assumption 2. (Mesh assumptions).** There exists a positive constant $\rho$ such that for any $E \in \{\Omega_h\}_h$

* $E$ is star-shaped with respect to a ball $B_E$ of radius $\ge \rho h_E$;

* any edge $e$ of $E$ has length $\ge \rho h_E$.

Given $\omega \subset \bar{\Omega}$ and $n \in \mathbb{N}$, we denote by $\mathbb{P}_n(\omega)$ the set of polynomials on $\omega$ of degree less or equal to $n$, with the convention that $\mathbb{P}_{-1}(\omega) = \{0\}$. A natural basis for the space $\mathbb{P}_n(E)$ is the set of normalized monomials $M_n(E) = \{m_\alpha$, with $\alpha = (\alpha_1, \alpha_2) \in \mathbb{N}^2$ such that $|\alpha| \le n\}$, where

$$m_{\alpha} = \prod_{i=1}^{2} \left( \frac{x_i - x_{E,i}}{h_E} \right)^{\alpha_i} \quad \text{and} \quad |\alpha| = \sum_{i=1}^{2} \alpha_i.$$

For any $e$ edge of $\Omega_h$, let $t_e$ and $n_e$ denote the tangent and the normal vectors to the edge $e$ respectively. Moreover, the normalized monomial set $M_n(e)$ is defined analogously as the span of all one-dimensional normalized monomials of degree up to $n$. For any $m \le n$, we denote with

$$\hat{\mathbb{P}}_{n\setminus m}(E) = \operatorname{span}\{m_{\alpha}, \text{ with } m+1 \le |\alpha| \le n\}.$$

For any $E \in \Omega_h$, the $L^2$-projection $\Pi_n^{0,E}: L^2(E) \to \mathbb{P}_n(E)$ is defined such that

$$\int_E q_n(v - \Pi_n^{0,E} v) dE = 0 \quad \text{for all } v \in L^2(E) \text{ and } q_n \in \mathbb{P}_n(E), \quad (13)$$

with obvious extension $\Pi_n^{0,E}: L^2(E) \to [\mathbb{P}_n(E)]^d$ and $\Pi_n^{0,E}: L^2(E) \to [\mathbb{P}_n(E)]^{d \times d}$ for vector and tensor functions, respectively. Moreover, the elliptic projection $\Pi_n^{\nabla,E}: W^{1,2}(E) \to \mathbb{P}_n(E)$ is given by

$$\left\{ \begin{array}{ll} \displaystyle \int_E \nabla q_n \cdot \nabla (v - \Pi_n^{\nabla,E} v) \, dE = 0 & \text{for all } v \in W^{1,2}(E) \text{ and } q_n \in \mathbb{P}_n(E), \\ \\ \displaystyle \int_{\partial E} (v - \Pi_n^{\nabla,E} v) \, ds = 0, & \end{array} \right.$$

with extension for vector fields $\Pi_n^{\nabla,E}: W^{1,2}(E) \to [\mathbb{P}_n(E)]^d$.

We also recall the following useful results:

* Trace inequality with scaling [26]: For any $E \in \Omega_h$ and for any function $v \in W^{1,r}(E)$ it holds

$$\|v\|_{L^r(\partial E)}^r \lesssim h_E^{-1} \|v\|_{L^r(E)}^r + h_E^{r-1} \|\nabla v\|_{L^r(E)}^r. \quad (14)$$

* Polynomial inverse estimate [26, Theorem 4.5.11]: Let $1 \le q, \ell \le \infty$ and let $s \ge 0$, then for any $E \in \Omega_h$

$$\|p_n\|_{W^{s,q}(E)} \lesssim h_E^{2/q - 2/\ell - s} \|p_n\|_{L^\ell(E)} \quad \text{for any } p_n \in \mathbb{P}_n(E). \quad (15)$$

At the global level, given $n \in \mathbb{N}$, $m \in \mathbb{R}_+$, and $l \in [1, +\infty)$, we introduce the piecewise regular spaces

* $\mathbb{P}_n(\Omega_h) = \{q \in L^2(\Omega) \text{ s.t. } q|_E \in \mathbb{P}_n(E) \text{ for all } E \in \Omega_h\}$,

* $W^{m,l}(\Omega_h) = \{v \in L^l(\Omega) \text{ s.t. } v|_E \in W^{m,l}(E) \text{ for all } E \in \Omega_h\}$,

equipped with the broken norm and seminorm

$$\|v\|_{W^{m,l}(\Omega_h)}^l = \sum_{E \in \Omega_h} \|v\|_{W^{m,l}(E)}^l, \quad |v|_{W^{m,l}(\Omega_h)}^l = \sum_{E \in \Omega_h} |v|_{W^{m,l}(E)}^l, \quad \text{if } 1 \le l < \infty. \quad (16)$$
Further we define the operator $\Pi_n^0: L^2(\Omega) \to \mathbb{P}_n(\Omega_h)$ such that $\Pi_n^0|_E = \Pi_n^{0,E}$ for any $E \in \Omega_h$.

Let $k \ge 1$ be the polynomial order of the method. We consider on each polygonal element $E \in \Omega_h$ the “enhanced” virtual space [3, 15, 16, 56]:

$$ 
\begin{align} 
U_h(E) = & \left\{ \mathbf{v}_h \in [C^0(\bar{E})]^2 \text{ s.t. } 
\begin{array}{@{}l@{}}
(i) \ \Delta \mathbf{v}_h + \nabla s \in \mathbf{x}^{\perp} \mathbb{P}_{k-1}(E), \text{ for some } s \in L_0^2(E), \\
(ii) \ \nabla \cdot \mathbf{v}_h \in \mathbb{P}_{k-1}(E), \\
(iii) \ \mathbf{v}_{h|e} \cdot \mathbf{n}_e \in \mathbb{P}_{\max\{2,k\}}(e), \mathbf{v}_{h|e} \cdot \mathbf{t}_e \in \mathbb{P}_k(e) \quad \forall e \in \partial E, \\
(iv) \ (\mathbf{v}_h - \Pi_k^{\nabla,E} \mathbf{v}_h, \mathbf{x}^{\perp} \hat{p}_{k-1})_E = 0 \quad \forall \hat{p}_{k-1} \in \hat{\mathbb{P}}_{(k-1)\setminus(k-3)}(E)
\end{array}
\right\}, 
\end{align} 
$$

where $\mathbf{x}^{\perp} = (x_2, -x_1)$. Next, we summarize the main properties of the space $U_h(E)$.

(P1) **Polynomial inclusion:** $[\mathbb{P}_k(E)]^2 \subseteq U_h(E)$;

(P2) **Degrees of freedom:** the following linear operators $D_U$ constitute a set of DoFs for $U_h(E)$:

D_U 1 the values of $\mathbf{v}_h$ at the vertexes of the polygon $E$,

D_U 2 the edge moments of $\mathbf{v}_h$ for every edge $e \in \partial E$,

$$ \frac{1}{|e|} \int_{e} \mathbf{v}_h \cdot \mathbf{t}_e m_{\alpha} \, ds, \quad \text{for any } m_{\alpha} \in \mathbb{M}_{k-2}(e), $$

$$ \frac{1}{|e|} \int_{e} \mathbf{v}_h \cdot \mathbf{n}_e m_{\alpha} \, ds \quad \text{for any } m_{\alpha} \in \mathbb{M}_{\max\{2,k\}-2}(e), $$

D_U 3 the moments of $\mathbf{v}_h$

$$ \frac{1}{|E|} \int_{E} \mathbf{v}_h \cdot \frac{m_{\alpha}}{h_{E}} (x_2 - x_{2,E}, -x_1 + x_{1,E}) \, dE \quad \text{for any } m_{\alpha} \in \mathbb{M}_{k-3}(E), $$

D_U 4 the moments of $\nabla \cdot \mathbf{v}_h$

$$ \frac{h_E}{|E|} \int_E (\nabla \cdot \mathbf{v}_h) m_{\alpha} \, dE \quad \text{for any } m_{\alpha} \in \mathbb{M}_{k-1}(E) \text{ with } |\alpha| > 0; $$

(P3) **Polynomial projections:** the DoFs $D_U$ allow us to compute the following linear operators:

$$ \Pi_k^{0,E}: U_h(E) \to [\mathbb{P}_k(E)]^2, \quad \Pi_{k-1}^{0,E}: \nabla U_h(E) \to [\mathbb{P}_{k-1}(E)]^{2 \times 2}. $$

The global velocity space $U_h = \{\mathbf{v}_h \in C^0(\Omega) \text{ s.t. } \mathbf{v}_h|_E \in U_h(E) \text{ for all } E \in \Omega_h\}$ is defined by gluing the local spaces with the obvious associated sets of global DoFs.

The discrete pressure space $P_h$ is given by the piecewise polynomial functions of degree $k-1$:

$$ P_h = \{q_h \in P \text{ s.t. } q_{h|E} \in \mathbb{P}_{k-1}(E) \text{ for all } E \in \Omega_h\}. \quad (18) $$

The couple of spaces $(U_h, P_h)$ is well known to be inf-sup stable in the classical Hilbertian setting [15, 16]. The inf-sup stability for $r > 2$ is proven below (Lemma 8). Let us introduce the discrete kernel

$$ Z_h = \{\mathbf{v}_h \in U_h \text{ s.t. } b(\mathbf{v}_h, q_h) = 0 \text{ for all } q_h \in P_h\} \quad (19) $$

and observe that, owing to (ii) in (17) and (18), $\nabla \cdot \mathbf{v}_h = 0$, for all $\mathbf{v}_h \in Z_h$.

Given any $\mathbf{v} \in W^{s,p}(E)$, with $p \in (1, \infty)$ and $s \in \mathbb{R}_+$, $s > 2/p$, we define its approximant $\mathbf{v}_I \in U_h$ as the unique function in $U_h$ that interpolates $\mathbf{v}$ with respect to the DoF set $D_U$. It is easy to check that, whenever $\nabla \cdot \mathbf{v} = 0$, then $\mathbf{v}_I \in Z_h$. Furthermore, the following approximation property is a trivial generalization of the results in [52].

**Lemma 4.** Let $E \in \Omega_h$, $n \in \mathbb{N}$, $\ell \in [1, \infty]$, $s \in \mathbb{R}_+$ and $\mathbf{v} \in W^{s,\ell}(E)$. Then

$$ |\mathbf{v} - \Pi_n^{0,E} \mathbf{v}|_{W^{m,\ell}(E)} \le h_E^{s-m} |\mathbf{v}|_{W^{s,\ell}(E)} \quad \text{for } 0 \le m \le s \le n+1. $$

Furthermore, given $\mathbf{v} \in W^{s,2}(E)$, $s > 1$, let $\mathbf{v}_I \in U_h$ be the interpolant of $\mathbf{v}$ defined above. It holds

$$ |\mathbf{v} - \mathbf{v}_I|_{W^{m,2}(E)} \le h_E^{s-m} |\mathbf{v}|_{W^{s,2}(E)} \quad \text{for } 1 < s \le k+1, m \in \{0, 1\}. $$

The first bound above extends identically to the scalar and tensor-valued cases.
Remark 5. We must note that, by known regularity results on Lipschitz domains, definition (17) guarantees $U_h(E) \subset W^{1,r}(E)$ for all $r \in (1, \infty)$ only if the polygonal element $E$ is convex [33, 46] or has a small Lipschitz constant [38]. Otherwise, if $d=2$ and $r>4$ or $d=3$ and $r>3$, such inclusion does not hold (see, e.g., [44, 54]). This is the reason why, in the following analysis, we make use of discrete norms (see definition (27) below). Additionally, assuming that all the mesh elements are convex, the a priori estimates established in Section 4 below imply error bounds with respect to standard $W^{1,r}$-norms (see Corollary 20 below).

## 3.2 Discrete problem

To formulate the VEM approximation of (5), we aim to construct a discrete version of the nonlinear form $a(\cdot, \cdot)$ in (6) and the approximation of the loading term $f$. For the latter, we define the discrete volumetric force as

$$f_h = \Pi_k^0 f, \qquad (20)$$

owing to the fact that $(f_h, v_h)$ is computable by property (P3). The discrete nonlinear form $a_h(\cdot, \cdot)$ is given as the sum of a consistency term and a stabilizing form that is suited for the non-linearity under consideration. Following [5], we consider a non-linear dofi-dofi stabilization $S(\cdot, \cdot): U_h \times U_h \to \mathbb{R}$ defined as

$$S(v_h, w_h) = \sum_{E \in \Omega_h} S^E(v_h, w_h) \quad \text{for all } v_h, w_h \in U_h, \qquad (21)$$

with $S^E(\cdot, \cdot): U_h(E) \times U_h(E) \to \mathbb{R}$ resembling the nonlinear law in (2), i.e.

$$S^E(v_h, w_h) = \bar{\mu}_E (\delta^\alpha + h_E^{-\alpha} |\chi(v_h)|^\alpha)^{\frac{r-2}{\alpha}} \chi(v_h) \cdot \chi(w_h), \qquad (22)$$

where $\bar{\mu}_E = \Pi_0^{0,E} \mu$ and $\chi: U_h(E) \to \mathbb{R}^{N_E}$, with $N_E$ denoting the dimension of $U_h(E)$, is the function that associates to each $v_h \in U_h(E)$ the vector of the local degrees of freedom in (P2). We remark that, according to the strong monotonicity and Hölder continuity of the Carreau-Yasuda law, any choice of $\alpha \in (0, \infty)$ in (22) give rise to an equivalent stabilization, in the sense that

$$\bar{\mu}_E (\delta^r + h_E^{-r} |\chi(v_h)|^r)^{\frac{r-2}{r}} |\chi(v_h)|^2 \lesssim S^E(v_h, v_h) \lesssim \bar{\mu}_E (\delta^r + h_E^{-r} |\chi(v_h)|^r)^{\frac{r-2}{r}} |\chi(v_h)|^2. \qquad (23)$$

Thus, we define the global form $a: U_h \times U_h \to \mathbb{R}$ such that

$$a_h(v_h, w_h) = \int_{\Omega} \sigma(\cdot, \Pi_{k-1}^0 \epsilon(v_h)) : \Pi_{k-1}^0 \epsilon(w_h) + S((I - \Pi_k^0) v_h, (I - \Pi_k^0) w_h) \quad \forall v_h, w_h \in U_h. \qquad (24)$$

We observe that, owing to property (P3), all projection operators appearing above are computable explicitly in terms of the velocity DOFs.

The virtual element discretization of Problem (5) is given by: Find $(u_h, p_h) \in U_h \times P_h$ such that

$$\begin{align} a_h(u_h, v_h) + b(v_h, p_h) &= \int_{\Omega} f_h \cdot v_h + \int_{\Gamma_N} g \cdot v_h & \forall v_h \in U_h, \tag{25} \\ b(u_h, q_h) &= 0 & \forall q_h \in P_h. \nonumber \end{align}$$

Recalling the definition of the discrete kernel $Z_h$ in (19), the previous problem can also be written in the kernel formulation: Find $u_h \in Z_h$ such that

$$a_h(u_h, v_h) = \int_{\Omega} f_h \cdot v_h + \int_{\Gamma_N} g \cdot v_h \quad \forall v_h \in Z_h. \qquad (26)$$

## 3.3 Well-posedness

This section aims to establish the well-posedness of the discrete problem (25). To do so, first we prove the inf-sup stability of the bilinear form $b(\cdot, \cdot)$ and then we investigate the continuity and monotonicity properties of $a_h(\cdot, \cdot)$.

We define, for all $v \in W^{1,r}(\Omega) \cup U_h$ with $r \in [2, \infty)$, the discrete quantity:

$$\|v\|_r := \|\Pi_{k-1}^0 \epsilon(v)\|_{L^r(\Omega)} + \sum_{E \in \Omega_h} h_E^{2-r} |\chi((I - \Pi_k^0)v)|^r. \qquad (27)$$
Note that, owing to Lemma 7 below, $\|\cdot\|_r$ defines a norm on $U_h$. In Section 4, we will measure the discretization error with respect to the quantity

$$ \|\mathbf{v}\|_{\delta,r}^r := \|\mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{v})\|_{\mathbb{L}^r(\Omega)}^r + S((I - \Pi_k^0)\mathbf{v}, (I - \Pi_k^0)\mathbf{v}), \quad (28) $$

which is not a norm since absolute homogeneity does not hold due to the dependence on $\delta$ of the stabilization term. However, for all $\delta \ge 0$ and $r \ge 2$, the error measure in (28) controls the discrete norm. This is established in the next Lemma.

**Lemma 6.** Given $\delta \ge 0$, the discrete norm defined in (27), and the error measure as in (28), the following inequalities hold:

$$ \|\mathbf{v}\|_r \le \|\mathbf{v}\|_{\delta,r} \lesssim (\delta + \|\mathbf{v}\|_r)^{\frac{r-2}{r}} \|\mathbf{v}\|_r^{\frac{2}{r}}. \quad (29) $$

*Proof.* The first inequality in (29) is a direct consequence of the definition of the stabilization function (22) and the fact that $r \ge 2$. To prove the second inequality, we set $\mathbf{v}^{\perp} = (I - \Pi_k^0)\mathbf{v}$, use the equivalence property (23), recall that $h_E^2 \approx |E|$, and apply the discrete $(\frac{r}{r-2}, \frac{r}{2})$-Hölder inequality, to infer

$$ \begin{aligned} S(\mathbf{v}^{\perp}, \mathbf{v}^{\perp}) &\le \sum_{E \in \Omega_h} (\delta^r + h_E^{-r} |\chi(\mathbf{v}^{\perp})|^r)^{\frac{r-2}{r}} |\chi(\mathbf{v}^{\perp})|^2 \\ &\le \sum_{E \in \Omega_h} (|E|\delta^r + h_E^{2-r} |\chi(\mathbf{v}^{\perp})|^r)^{\frac{r-2}{r}} (h_E^{2-r} |\chi(\mathbf{v}^{\perp})|^r)^{\frac{2}{r}} \\ &\le (|\Omega|\delta^r + \sum_{E \in \Omega_h} h_E^{2-r} |\chi(\mathbf{v}^{\perp})|^r)^{\frac{r-2}{r}} (\sum_{E \in \Omega_h} h_E^{2-r} |\chi(\mathbf{v}^{\perp})|^r)^{\frac{2}{r}} \le (\delta + \|\mathbf{v}\|_r)^{r-2} \|\mathbf{v}\|_r^2. \end{aligned} $$

Moreover, owing to $\delta \ge 0$ and $r \ge 2$, one also has $\|\Pi_{k-1}^0 \epsilon(\mathbf{v})\|_{\mathbb{L}^r(\Omega)} \le (\delta + \|\mathbf{v}\|_r)^{r-2} \|\mathbf{v}\|_r^2$. $\square$

### 3.3.1 Discrete inf-sup condition

We here prove a discrete inf-sup condition analogous to the continuous one. The difference with respect to the analogous condition proved in [5] is that here we make use of the discrete norm in (27).

We recall the following Lemma established in [52] (see also [5, Lemmas 8, 9]).

**Lemma 7.** Let the mesh regularity in Assumption 2 hold. For any $E \in \Omega_h$ we have

$$ |\mathbf{v}_h|_{W^{1,2}(E)} \lesssim |\chi(\mathbf{v}_h)| \lesssim |\mathbf{v}_h|_{W^{1,2}(E)} \quad \text{for all } \mathbf{v}_h \in U_h(E) \text{ s.t. } \Pi_k^{0,E} \mathbf{v}_h = 0. $$

**Lemma 8 (Discrete inf-sup).** Let the mesh regularity assumptions stated in Assumption 2 hold. Then, for any $r \in [2, \infty)$ it exists a constant $\bar{\beta}(r)$, such that

$$ \inf_{q_h \in P_h} \sup_{w_h \in U_h} \frac{b(w_h, q_h)}{\|q_h\|_{L^{r'}(\Omega)} \|w_h\|_r} \ge \bar{\beta}(r) > 0. $$

*Proof.* The proof is a modification of that of Lemma 16 in [5]; to shorten the exposition, we refer here directly to the notation and derivations in that article. We consider the same Fortin operator $\Pi^F$ introduced in the lemma above. The only difference in the present proof is that when showing the continuity of $\Pi^F$ from $W^{1,r}(\Omega)$ into $U_h$, the latter space is equipped with the $\|\cdot\|_r$ norm instead of the $W^{1,r}$ norm.

Let now $w \in W^{1,r}(\Omega)$. First, by definition of the stabilization form and some trivial algebra, then combining a standard inverse estimate for polynomials with the continuity of the $\Pi_{k-1}^0$ operator and recalling Lemma 7, we obtain

$$ \begin{aligned} \|\Pi^F w\|_r^r &= \sum_{E \in \Omega_h} (\|\Pi_{k-1}^0 \epsilon(\Pi^F w)\|_{\mathbb{L}^r(E)}^r + h_E^{2-r} |\chi(I - \Pi_k^0)\Pi^F w|^r) \\ &\lesssim \sum_{E \in \Omega_h} h_E^{2-r} (\|\epsilon(\Pi^F w)\|_{\mathbb{L}^2(E)}^r + |(I - \Pi_k^0)\Pi^F w|_{W^{1,2}(E)}^r). \end{aligned} $$

Using the continuity of $\Pi_k^0$, the “local” continuity of $\Pi^F$ in $W^{1,2}$ (see for instance Lemma 16 in [5]) and a Hölder inequality we get

$$ \|\Pi^F w\|_r^r \lesssim \sum_{E \in \Omega_h} h_E^{2-r} |\Pi^F w|_{W^{1,2}(E)}^r \lesssim \sum_{E \in \Omega_h} h_E^{2-r} |w|_{W^{1,2}(\omega_E)}^r \lesssim \sum_{E \in \Omega_h} |w|_{W^{1,r}(\omega_E)}^r, $$

where $\omega_E$ is the union of all elements sharing at least a vertex with E. The above bound shows the required continuity for the $\Pi^F$ operator. $\square$
### 3.3.2 Properties of the stabilization

We establish the continuity and monotonicity of the local stabilization form.

**Lemma 9** (Hölder continuity and strong monotonicity of $S^E(\cdot, \cdot)$). Let the mesh regularity assumptions stated in Assumption 2 hold. Let $\mathbf{u}_h, \mathbf{w}_h \in \mathbf{U}_h(E)$ and set $\mathbf{e}_h = \mathbf{u}_h - \mathbf{w}_h$. Then, for all $\mathbf{v}_h \in \mathbf{U}_h(E)$, there holds

$$|S^E(\mathbf{u}_h, \mathbf{v}_h) - S^E(\mathbf{w}_h, \mathbf{v}_h)| \lesssim (\delta^r + h_E^{-r} |\chi(\mathbf{u}_h)|^r + h_E^{-r} |\chi(\mathbf{w}_h)|^r)^{\frac{r-2}{r}} |\chi(\mathbf{e}_h)| |\chi(\mathbf{v}_h)|. \quad (30)$$

Moreover, there holds

$$S^E(\mathbf{u}_h, \mathbf{e}_h) - S^E(\mathbf{w}_h, \mathbf{e}_h) \geq S^E(\mathbf{e}_h, \mathbf{e}_h) \geq h_E^{2-r} |\chi(\mathbf{e}_h)|^r. \quad (31)$$

*Proof.* (i) Hölder continuity. Recalling the definition of $S^E(\cdot, \cdot)$ in (22) and employing (3a), bound (30) is derived as follows

$$\begin{aligned} |S^E(\mathbf{u}_h, \mathbf{v}_h) - S^E(\mathbf{w}_h, \mathbf{v}_h)| & \lesssim |(\delta^\alpha + h_E^{-\alpha} |\chi(\mathbf{u}_h)|^\alpha)^{\frac{r-2}{\alpha}} \chi(\mathbf{u}_h) - (\delta^\alpha + h_E^{-\alpha} |\chi(\mathbf{w}_h)|^\alpha)^{\frac{r-2}{\alpha}} \chi(\mathbf{w}_h)| |\chi(\mathbf{v}_h)| \\ & \lesssim (\delta^r + h_E^{-r} |\chi(\mathbf{u}_h)|^r + h_E^{-r} |\chi(\mathbf{w}_h)|^r)^{\frac{r-2}{r}} |\chi(\mathbf{e}_h)| |\chi(\mathbf{v}_h)|. \end{aligned}$$

(ii) Strong monotonicity. We recall the strong monotonicity bound stating that, for all $\mathbf{x}, \mathbf{y} \in \mathbb{R}^n$, it holds

$$|\mathbf{x} - \mathbf{y}|^2 (\delta^r + |\mathbf{x}|^r + |\mathbf{y}|^r)^{\frac{r-2}{r}} \lesssim \left\{ (\delta^\alpha + |\mathbf{x}|^\alpha)^{\frac{r-2}{\alpha}} \mathbf{x} - (\delta^\alpha + |\mathbf{y}|^\alpha)^{\frac{r-2}{\alpha}} \mathbf{y} \right\} \cdot (\mathbf{x} - \mathbf{y}). \quad (32)$$

Inserting $\delta^\alpha \ge 0$, and using a triangle inequality together with the fact that $(x+y)^{\frac{\alpha}{r}} \lesssim x^{\frac{\alpha}{r}} + y^{\frac{\alpha}{r}}$ for all $x, y \in [0, \infty)$, and employing (32) with $\mathbf{x} = \chi(\mathbf{u}_h)$, $\mathbf{y} = \chi(\mathbf{w}_h)$, and recalling the definition of $S^E(\cdot, \cdot)$ in (22) we infer

$$\begin{aligned} h_E^{2-r} |\chi(\mathbf{e}_h)|^r & \lesssim |\chi(\mathbf{e}_h)|^2 (\delta^\alpha + h_E^{-\alpha} |\chi(\mathbf{e}_h)|^\alpha)^{\frac{r-2}{\alpha}} \approx S^E(\mathbf{e}_h, \mathbf{e}_h) \\ & \lesssim |\chi(\mathbf{e}_h)|^2 (\delta^r + h_E^{-r} |\chi(\mathbf{u}_h)|^r + h_E^{-r} |\chi(\mathbf{w}_h)|^r)^{\frac{r-2}{r}} \\ & \lesssim ((\delta^\alpha + h_E^{-\alpha} |\chi(\mathbf{u}_h)|^\alpha)^{\frac{r-2}{\alpha}} \chi(\mathbf{u}_h) - (\delta^\alpha + h_E^{-\alpha} |\chi(\mathbf{w}_h)|^\alpha)^{\frac{r-2}{\alpha}} \chi(\mathbf{w}_h)) \cdot \chi(\mathbf{e}_h) \\ & \lesssim S^E(\mathbf{u}_h, \mathbf{e}_h) - S^E(\mathbf{w}_h, \mathbf{e}_h). \end{aligned} \quad (33)$$

□

Next, we prove the following result that be useful in the sequel.

**Proposition 10.** For all $\mathbf{u}_h, \mathbf{w}_h \in \mathbf{U}_h$,

$$|S(\mathbf{u}_h, \mathbf{w}_h)| \lesssim (\delta^r + S(\mathbf{u}_h, \mathbf{u}_h))^{\frac{r-2}{2r}} S(\mathbf{u}_h, \mathbf{u}_h)^{\frac{1}{2}} S(\mathbf{w}_h, \mathbf{w}_h)^{\frac{1}{r}}. \quad (34)$$

*Proof.* Employing the discrete 3-terms $(\frac{2r}{r-2}, 2, r)$-Hölder inequality together with $h_E^2 \approx |E|$, also recalling (31), we obtain

$$\begin{aligned} |S(\mathbf{u}_h, \mathbf{w}_h)| & \lesssim \sum_{E \in \Omega_h} (\delta^r + h_E^{-r} |\chi(\mathbf{u}_h)|^r)^{\frac{r-2}{r}} |\chi(\mathbf{u}_h)| |\chi(\mathbf{w}_h)| \\ & = \sum_{E \in \Omega_h} (h_E^2 \delta^r + h_E^{-2-r} |\chi(\mathbf{u}_h)|^r)^{\frac{r-2}{2r}} ((\delta^r + h_E^{-r} |\chi(\mathbf{u}_h)|^r)^{\frac{r-2}{r}} |\chi(\mathbf{u}_h)|^2)^{\frac{1}{2}} (h_E^{-2-r} |\chi(\mathbf{w}_h)|^r)^{\frac{1}{r}} \\ & \lesssim \left( \sum_{E \in \Omega_h} (|E| \delta^r + h_E^{-2-r} |\chi(\mathbf{u}_h)|^r) \right)^{\frac{r-2}{2r}} \left( \sum_{E \in \Omega_h} (\delta^r + h_E^{-r} |\chi(\mathbf{u}_h)|^r)^{\frac{r-2}{r}} |\chi(\mathbf{u}_h)|^2 \right)^{\frac{1}{2}} \left( \sum_{E \in \Omega_h} h_E^{-2-r} |\chi(\mathbf{w}_h)|^r \right)^{\frac{1}{r}} \\ & \lesssim (\delta^r + S(\mathbf{u}_h, \mathbf{u}_h))^{\frac{r-2}{2r}} S(\mathbf{u}_h, \mathbf{u}_h)^{\frac{1}{2}} S(\mathbf{w}_h, \mathbf{w}_h)^{\frac{1}{r}}. \end{aligned} \quad (35)$$

□
### 3.3.3 Properties of the discrete viscous term

Hinging on the results of the previous section, we establish the continuity and monotonicity properties of the discrete viscous function $a_h(\cdot, \cdot)$.

**Proposition 11** (Hölder continuity and strong monotonicity of $a_h(\cdot, \cdot)$). Let $u_h$, $w_h \in U_h$ and set $e_h = u_h - w_h \in U_h$. Then for any $v_h \in U_h$ there holds

$$|a_h(\mathbf{u}_h, \mathbf{v}_h) - a_h(\mathbf{w}_h, \mathbf{v}_h)| \lesssim (\delta^r + \|\mathbf{u}_h\|_r^r + \|\mathbf{w}_h\|_r^r)^{\frac{r-2}{r}} \|\mathbf{e}_h\|_r \|\mathbf{v}_h\|_r. \quad (36)$$

Moreover, it holds

$$a_h(\mathbf{u}_h, \mathbf{e}_h) - a_h(\mathbf{w}_h, \mathbf{e}_h) \geq \|\mathbf{e}_h\|_{\delta,r}^r \geq \|\mathbf{e}_h\|_r^r. \quad (37)$$

*Proof.* For the sake of conciseness, for all $v_h \in U_h$, we let $v_h^{\perp} := (I - \Pi_k^0) v_h$.

(i) Hölder continuity. Recalling the definition of $a_h(\cdot, \cdot)$ in (24) we can write

$$|a_h(\mathbf{u}_h, \mathbf{v}_h) - a_h(\mathbf{w}_h, \mathbf{v}_h)| \leq |T_1| + |T_2| \quad (38)$$

where

$$T_1 := \int_{\Omega} (\sigma(\cdot, \Pi_{k-1}^0 \epsilon(\mathbf{u}_h)) - \sigma(\cdot, \Pi_{k-1}^0 \epsilon(\mathbf{w}_h))) : \Pi_{k-1}^0 \epsilon(\mathbf{v}_h), \\ T_2 = S(\mathbf{u}_h^{\perp}, \mathbf{v}_h^{\perp}) - S(\mathbf{w}_h^{\perp}, \mathbf{v}_h^{\perp}).$$

Following the lines of (9) and recalling the definition of the discrete norm in (27) it is inferred that

$$\begin{aligned} |T_1| &\le \int_{\Omega} |(\sigma(\cdot, \Pi_{k-1}^0 \epsilon(\mathbf{u}_h)) - \sigma(\cdot, \Pi_{k-1}^0 \epsilon(\mathbf{w}_h))) : \Pi_{k-1}^0 \epsilon(\mathbf{v}_h)| \\ &\le \left( \int_{\Omega} \delta^r + |\Pi_{k-1}^0 \epsilon(\mathbf{u}_h)|^r + |\Pi_{k-1}^0 \epsilon(\mathbf{w}_h)|^r \right)^{\frac{r-2}{r}} \|\Pi_{k-1}^0 \epsilon(\mathbf{e}_h)\|_{L^r(\Omega)} \|\Pi_{k-1}^0 \epsilon(\mathbf{v}_h)\|_{L^r(\Omega)} \quad (39) \\ &\le (|\Omega|\delta^r + \|\mathbf{u}_h\|_r^r + \|\mathbf{w}_h\|_r^r)^{\frac{r-2}{r}} \|\mathbf{e}_h\|_r \|\mathbf{v}_h\|_r. \end{aligned}$$

Then, applying (30) together with a discrete Hölder inequality with exponents $(\frac{r}{r-2}, r, r)$, we infer

$$\begin{aligned} |T_2| &\lesssim \sum_{E \in \Omega_h} |S^E(\mathbf{u}_h^{\perp}, \mathbf{v}_h^{\perp}) - S^E(\mathbf{w}_h^{\perp}, \mathbf{v}_h^{\perp})| \\ &\lesssim \left( \sum_{E \in \Omega_h} \left( h_E^2 \delta^r + h_E^{2-r} |\chi(\mathbf{u}_h^{\perp})|^r + h_E^{2-r} |\chi(\mathbf{w}_h^{\perp})|^r \right) \right)^{\frac{r-2}{r}} \left( \sum_{E \in \Omega_h} h_E^{2-r} |\chi(\mathbf{e}_h^{\perp})|^r \right)^{\frac{1}{r}} \left( \sum_{E \in \Omega_h} h_E^{2-r} |\chi(\mathbf{v}_h^{\perp})|^r \right)^{\frac{1}{r}} \\ &\lesssim (|\Omega|\delta^r + \|\mathbf{u}_h\|_r^r + \|\mathbf{w}_h\|_r^r)^{\frac{r-2}{r}} \|\mathbf{e}_h\|_r \|\mathbf{v}_h\|_r, \end{aligned} \quad (40)$$

where we used the fact that $h_E^2 \approx |E|$. The proof follows inserting (39) and (40) in (38).

(ii) Strong monotonicity. Recalling (3b) and using the inequality in (10), we get

$$\begin{aligned} \|\Pi_{k-1}^0 \epsilon(\mathbf{e}_h)\|_{L^r(\Omega)}^r &\le \int_{\Omega} (\delta + |\Pi_{k-1}^0 \epsilon(\mathbf{u}_h - \mathbf{w}_h)|)^{r-2} |\Pi_{k-1}^0 \epsilon(\mathbf{u}_h - \mathbf{w}_h)|^2 \\ &\le \int_{\Omega} (\sigma(\cdot, \Pi_{k-1}^0 \epsilon(\mathbf{u}_h)) - \sigma(\cdot, \Pi_{k-1}^0 \epsilon(\mathbf{w}_h))) : \Pi_{k-1}^0 \epsilon(\mathbf{e}_h). \end{aligned}$$

Moreover, from the first inequality in (31) it is readily inferred that

$$S(\mathbf{e}_h^{\perp}, \mathbf{e}_h^{\perp}) \lesssim S(\mathbf{u}_h^{\perp}, \mathbf{e}_h^{\perp}) - S(\mathbf{w}_h^{\perp}, \mathbf{e}_h^{\perp})$$

The assertion follows by summing the previous bounds, recalling the definition of the $\|\cdot\|_{\delta,r}$-norm and using (29). $\square$

### 3.3.4 Main results

We are now ready to prove the well-posedness of the discrete Virtual Element problem (25).

**Theorem 12** (Existence and uniqueness). For any $r \in [2, \infty)$, there exists a unique solution $u_h \in Z_h$ to the discrete problem (26).
Proof. (i) Existence. Let the mesh $\Omega_h$ be fixed. We equip the space $Z_h$ with the $(\cdot, \cdot)_{W^{1,2}}$-inner product and induced norm $\|\cdot\|_{W^{1,2}(\Omega)}$. Owing to the equivalence of norms in finite-dimensional spaces, we have

$$ \|v_h\|_{W^{1,2}(\Omega)} \lesssim C_h \|v_h\|_{r}, \qquad (41) $$

with the positive constant $C_h$ depending on the mesh size $h$. We also define the nonlinear function $\Phi_h: Z_h \to Z_h$ such that

$$ (\Phi_h(v_h), w_h)_{W^{1,2}(\Omega)} := a_h(v_h, w_h), \quad \forall v_h, w_h \in Z_h. $$

The strong monotonicity of $a_h(\cdot, \cdot)$ established in Lemma 11 together with (41) leads to, for any $v_h \in Z_h$,

$$ \lim_{\|v_h\|_{W^{1,2}(\Omega)} \to \infty} \frac{(\Phi_h(v_h), v_h)_{W^{1,2}(\Omega)}}{\|v_h\|_{W^{1,2}(\Omega)}} \ge \lim_{\|v_h\|_{W^{1,2}(\Omega)} \to \infty} \frac{\|v_h\|_{r}^{r}}{\|v_h\|_{W^{1,2}(\Omega)}} \ge C_h^{-1} \lim_{\|v_h\|_{W^{1,2}(\Omega)} \to \infty} \frac{\|v_h\|_{r}^{r-1}}{\|v_h\|_{r}^{r}} \to \infty. $$

By applying [34, Theorem 3.3], the previous result shows that the operator $\Phi_h$ is onto. As a result, there exists $u_h \in Z_h$ such that $\Phi_h(u_h) = z_h$, with $z_h \in Z_h$ defined such that

$$ (z_h, w_h)_{W^{1,2}(\Omega)} = \int_{\Omega} f_h \cdot w_h + \int_{\Gamma_N} g \cdot w_h \quad \forall w_h \in Z_h. $$

Thanks to the definition of $\Phi_h$, this implies that $u_h$ is a solution to the discrete problem (26).

(ii) Uniqueness. Let $u_{h,1}, u_{h,2} \in Z_h$ solve (26). Subtracting (26) for $u_{h,2}$ from (26) for $u_{h,1}$ and then taking $v_h = u_{h,1} - u_{h,2}$ as test function, we obtain

$$ a_h(u_{h,1}, u_{h,1} - u_{h,2}) - a_h(u_{h,2}, u_{h,1} - u_{h,2}) = 0. $$

Hence, using again the strong monotonicity of $a_h(\cdot, \cdot)$ with $e_h = u_{h,1} - u_{h,2}$, we get $\|u_{h,1} - u_{h,2}\|_r^r = 0$, that implies $u_{h,1} = u_{h,2}$. $\square$

The next result is derived by using the discrete inf-sup condition established in Lemma 8 and the equivalence of the discrete problems (25) and (26).

**Corollary 13** (Well-posedness of (25)). For any $r \in [2, \infty)$, there exists a unique solution $(u_h, p_h) \in U_h \times P_h$ to the discrete problem (25).

**Remark 14** (Stability estimates). A priori estimates for the unique discrete velocity and pressure fields $(u_h, p_h)$ solving problem (25) can be obtained by reasoning as in the proof of Proposition 3. The estimate

$$ \|u_h\|_r \lesssim \left( \|f_h\|_{L^{r'}(\Omega)} + \|g\|_{L^{r'}(\Gamma_N)} \right)^{\frac{1}{r-1}} $$

hinges on the monotonicity property in Proposition 11; whereas the estimate for the discrete pressure field follows from the inf-sup condition in Lemma 8 and again Proposition 11 and reads

$$ \|p_h\|_{L^{r'}(\Omega)} \lesssim \left( \|f_h\|_{L^{r'}(\Omega)} + \|g\|_{L^{r'}(\Gamma_N)} \right) + \delta^{r-2} \left( \|f_h\|_{L^{r'}(\Omega)} + \|g\|_{L^{r'}(\Gamma_N)} \right)^{\frac{1}{r-1}}. $$

# 4 A priori error analysis

This section is devoted to the a priori error analysis.

## 4.1 Additional properties of the stress-strain law

We recall some important results regarding the stress-strain relation that are instrumental for the a priori analysis of the scheme. We mainly follow [21, Section 3] and [43, Section 2]. For $r \ge 2$, we introduce, with $a \ge 0$, the shifted functions $\varphi_a(t) = \int_0^t (a+s)^{r-2}s$ ds. The following Lemma provides important properties of the shifted functions $\varphi_a$. We refer the reader to [35, Lemmata 28–32] and [36, Corollary 26] for the detailed proof.

**Lemma 15** (Young type inequalities). Let $r \ge 2$. For all $\varepsilon > 0$ there exists $C_\varepsilon > 0$ only depending on $r$ such that for all $s, t, a \ge 0$ and all $\tau, \eta \in \mathbb{R}^{d \times d}$ there holds

$$ s\varphi'_{a}(t) + t\varphi'_{a}(s) \le \varepsilon\varphi_{a}(s) + C_{\varepsilon}\varphi_{a}(t), \qquad (42a) $$

$$ \varphi_{a+|\tau|}(t) \le \varepsilon\varphi_{a+|\eta|}(|\tau - \eta|) + C_{\varepsilon}\varphi_{a+|\eta|}(t). \qquad (42b) $$
Proof. In the following we only sketch the proof that $C_\epsilon$ does not depend on $a$, thanks to [35, Lemmata 28–32] and [36, Corollary 26]. We define

$$ \Delta(\varphi_a) := \min\{\alpha > 0 : \varphi_a(2t) \le \alpha \varphi_a(t) \quad \forall t \in \mathbb{R}_0^+\}, \quad (43a) $$

$$ \Delta(\{\varphi_a, \varphi_a^*\}) := \max(\Delta(\varphi_a), \Delta(\varphi_a^*)), \quad (43b) $$

where $\varphi_a^*$ is the Fenchel conjugate of $\varphi_a$, i.e., $(\varphi_a')'(\varphi_a'(t)) = t$ for all $t \in \mathbb{R}_0^+$. Now, we only need to show that $\Delta(\{\varphi_a, \varphi_a^*\})$ is bounded independently of $a$. We have for all $t \in \mathbb{R}_0^+$,

$$ \varphi_a(2t) = \int_0^{2t} (a+s)^{r-2}s \, ds = 4 \int_0^t (a+2s)^{r-2}s \, ds \le 2^r \int_0^t (a+s)^{r-2}s \, ds = 2^r \varphi_a(t) \quad (44) $$

so $\Delta(\varphi_a) \le 2^r$. Moreover, using the fact that $\varphi_a'(2t) = (a+2t)^{r-2}2t \ge 2(a+t)^{r-2}t = 2\varphi_a'(t)$, we obtain, $\varphi_a'(2\varphi_a'(2t)) \ge 2\varphi_a'(\varphi_a'(t)) = 2t$, so $2\varphi_a'(t) \ge \varphi_a'(2t)$, thus $2\varphi_a^*(t) \ge \frac{1}{2}\varphi_a^*(2t)$ by integrating, hence, $\varphi_a^*(2t) \le 4\varphi_a^*(t)$ and we obtain $\Delta(\varphi_a^*) \le 4$. Therefore,

$$ \Delta(\{\varphi_a, \varphi_a^*\}) \le \max(2^r, 4) = 2^r. \quad (45) $$

□

The next result showing the equivalence of several quantities is strictly related to the continuity and monotonicity assumptions given in Assumption 1. The proofs of the next lemma can be found in [43, Section 2.3]. The lemma here below applies to any (scalar or tensor valued) function $\sigma$ which satisfies Assumption 1. In the following, with a slight abuse of notation, we will apply such lemma both to the constitutive law $\sigma$ but also to the auxiliary scalar function $\sigma(\tau) = (\delta + |\tau|)^{r-2}\tau$.

**Lemma 16.** Let $\sigma$ satisfy (3) for $r \in [2, \infty)$ and $\delta \ge 0$. Then, uniformly for all $\tau, \eta \in \mathbb{R}_s^{d \times d}$ and all $v, w \in U$ there hold

$$ |\sigma(\cdot, \tau) - \sigma(\cdot, \eta)| \approx (\delta + |\tau| + |\eta|)^{r-2} |\tau - \eta| \approx \varphi'_{\delta+|\tau|}(|\tau - \eta|), \quad (46a) $$

$$ (\sigma(\cdot, \tau) - \sigma(\cdot, \eta)) : (\tau - \eta) \approx (\delta + |\tau| + |\eta|)^{r-2} |\tau - \eta|^2 \approx \varphi_{\delta+|\tau|}(|\tau - \eta|), \quad (46b) $$

where the hidden constants only depend on $\sigma_c$, $\sigma_m$ in Assumption 1 and on $r$.

## 4.2 A priori error estimate: velocity

We start by a simple lemma; the proof is perhaps different from expected since $U_h(E) \subset W^{1,2}(E)$ but $U_h(E) \not\subset W^{1,r}(E)$ for non-convex elements $E$ and large $r$, see Remark 5.

**Lemma 17.** Let the mesh regularity in Assumption 2 hold. Let $E \in \Omega_h$ and $r \in [2, \infty]$. Let $v \in W^{s,r}(E)$, $1 < s \le k+1$, and $v_I \in U_h$ denote the interpolant of $v$ previously introduced. Then it holds

$$ \|\epsilon(v) - \Pi_{k-1}^{0,E} \epsilon(v_I)\|_{L^r(E)} \le h_E^{s-1} |v|_{W^{s,r}(E)}, \quad (47) $$

$$ |\chi((I - \Pi_k^0)v_I)| \le h_E^{s-2/r} |v|_{W^{s,r}(E)}. \quad (48) $$

*Proof.* We start by some trivial manipulation and afterwards apply the first bound in Lemma 4 together with the polynomial inverse estimate (15), obtaining

$$ \begin{aligned} \|\epsilon(v) - \Pi_{k-1}^{0,E} \epsilon(v_I)\|_{L^r(E)} &\le \|\epsilon(v) - \Pi_{k-1}^{0,E} \epsilon(v)\|_{L^r(E)} + \|\Pi_{k-1}^{0,E} \epsilon(v - v_I)\|_{L^r(E)} \\ &\lesssim h_E^{s-1} |v|_{W^{s,r}(E)} + |E|^{1/r-1/2} \|\Pi_{k-1}^{0,E} \epsilon(v - v_I)\|_{L^2(E)}. \end{aligned} $$

We conclude the proof of the first bound by the $L^2(E)$ continuity of $\Pi_{k-1}^{0,E}$, the interpolation estimates in Lemma 4 and finally a Hölder inequality on the element:

$$ \begin{aligned} \|\epsilon(v) - \Pi_{k-1}^{0,E} \epsilon(v_I)\|_{L^r(E)} &\lesssim h_E^{s-1} |v|_{W^{s,r}(E)} + |E|^{1/r-1/2} h_E^{s-1} |v|_{W^{s,2}(E)} \\ &\lesssim h_E^{s-1} |v|_{W^{s,r}(E)}. \end{aligned} $$

In order to deal with the second bound, we first apply Lemma 7, then some obvious manipulations, finally Lemma 4 and an inverse estimate for polynomials. We obtain

$$ \begin{aligned} |\chi(I - \Pi_k^0)v_I| &\lesssim |(I - \Pi_k^0)v_I|_{W^{1,2}(E)} \le |v - v_I|_{W^{1,2}(E)} + |(I - \Pi_k^0)v|_{W^{1,2}(E)} + |\Pi_k^0(v - v_I)|_{W^{1,2}(E)} \\ &\lesssim h_E^{s-1} |v|_{W^{s,2}(E)} + h_E^{-1} \|\Pi_k^{0,E}(v - v_I)\|_{L^2(E)}. \end{aligned} $$
We now recall the $L^2(E)$ continuity of $\Pi_k^{0,E}$, again Lemma 4 and a Hölder inequality:

$$|\chi(I - \Pi_k^0)\nu_I)| \lesssim h_E^{s-1} |\nu|_{W^{s,2}(E)} \lesssim h_E^{s-2/r} |\nu|_{W^{s,r}(E)}.$$

□

Moreover, we state the following interpolation lemma.

**Lemma 18.** Let the mesh regularity in Assumption 2 hold. Let $\nu \in W^{s,r}(\Omega_h)$, with $1 < s \le k+1$ and $r \ge 2$. Let $\nu_I \in U_h$ denote the interpolant of $\nu$ previously introduced. Then it holds

$$\|\mathbf{u} - \mathbf{u}_I\|_{\delta,r}^r \lesssim (\delta + h^{s-1}|\mathbf{u}|_{W^{s,r}(\Omega_h)})^{r-2} h^{2(s-1)} |\mathbf{u}|_{W^{s,r}(\Omega_h)}^2. \quad (49)$$

*Proof.* The proof will be presented briefly, since it essentially makes use of techniques already developed in previous results of this contribution. Let now $e_I = u - u_I$ and, as usual, $e_I^\perp := (I - \Pi_k^0)e_I$. By definition we have

$$\|e_I\|_{\delta,r}^r = \|\Pi_{k-1}^0 e_I\|_{L^r(\Omega_h)}^r + S(e_I^\perp, e_I^\perp) := T_1 + T_2.$$

The first term on the right-hand side is easily bounded by the triangle inequality, polynomial approximation estimates and Lemma 17, obtaining

$$T_1 \lesssim h^{r(s-1)} |\mathbf{u}|_{W^{s,r}(\Omega_h)}^r.$$

By definition of the stabilization form, first by a trivial manipulation, then by a discrete $(\frac{r}{2}, \frac{r}{r-2})$-Hölder inequality, we write (recall $|E| \approx h_E^2$)

$$\begin{aligned}
T_2 &\approx \sum_{E \in \Omega_h} (\delta + h_E^{-1} |\chi(e_I^\perp)|)^{r-2} |\chi(e_I^\perp)|^2 \\
&= \sum_{E \in \Omega_h} (h_E^{2/r} \delta + h_E^{2/r-1} |\chi(e_I^\perp)|)^{r-2} (h_E^{(2-r)/r} |\chi(e_I^\perp)|)^2 \\
&\lesssim \left( \sum_{E \in \Omega_h} h_E^2 \delta^r + h_E^{2-r} |\chi(e_I^\perp)|^r \right)^{\frac{r-2}{r}} \left( \sum_{E \in \Omega_h} h_E^{2-r} |\chi(e_I^\perp)|^r \right)^{\frac{2}{r}} \\
&\lesssim \left( \delta^r + \sum_{E \in \Omega_h} h_E^{2-r} |\chi(e_I^\perp)|^r \right)^{\frac{r-2}{r}} \left( \sum_{E \in \Omega_h} h_E^{2-r} |\chi(e_I^\perp)|^r \right)^{\frac{2}{r}}.
\end{aligned} \quad (50)$$

We now apply Lemma 7 and approximation properties for polynomials

$$\begin{aligned}
T_2 &\lesssim (\delta^r + \sum_{E \in \Omega_h} h_E^{2-r} |e_I^\perp|^r |W^{1,2}(E)|^{\frac{r-2}{r}}) (\sum_{E \in \Omega_h} h_E^{2-r} |e_I^\perp|^r |W^{1,2}(E)|^{\frac{2}{r}}) \\
&\lesssim (\delta^r + \sum_{E \in \Omega_h} h_E^{2-r+r(s-1)} |\mathbf{u}|_{W^{s,2}(E)}^r)^{\frac{r-2}{r}} (\sum_{E \in \Omega_h} h_E^{2-r+r(s-1)} |\mathbf{u}|_{W^{s,2}(E)}^r)^{\frac{2}{r}} \\
&\le (\delta^r + \sum_{E \in \Omega_h} h_E^{r(s-1)} |\mathbf{u}|_{W^{s,r}(E)}^r)^{\frac{r-2}{r}} (\sum_{E \in \Omega_h} h_E^{r(s-1)} |\mathbf{u}|_{W^{s,r}(E)}^r)^{\frac{2}{r}} \\
&\le (\delta + h^{s-1} |\mathbf{u}|_{W^{s,r}(\Omega_h)})^{r-2} (h^{s-1} |\mathbf{u}|_{W^{s,r}(\Omega_h)})^{2}.
\end{aligned}$$

The proof is concluded trivially by combining the bounds above.

□

We now present the main result of this section (see also the important Remark 25).

**Theorem 19.** Let $\mathbf{u}$ be the solution of problem (7) and let $\mathbf{u}_h$ be the solution of problem (26). Assume that $\mathbf{u} \in W^{k_1+1,r}(\Omega_h)$, $\sigma(\cdot, \epsilon(\mathbf{u})) \in W^{k_2,r'}(\Omega_h)$, $f \in W^{k_3+1,r'}(\Omega_h)$ for some positive integers $k_1, k_2, k_3 \le k$. Let the mesh regularity assumptions stated in Assumption 2 hold. Then, we have

$$\|\mathbf{u} - \mathbf{u}_h\|_{\delta,r} \lesssim (\delta + h^{k_1} R_1 + R_4)^{\frac{r-2}{r}} h^{\frac{2k_1}{r}} R_1^{\frac{2}{r}} + h^{\frac{k_2}{r-1}} R_2^{\frac{1}{r-1}} + h^{\frac{k_3+2}{r-1}} R_3^{\frac{1}{r-1}}, \quad (51)$$

where the regularity terms are

$$\begin{aligned}
R_1 &= |\mathbf{u}|_{W^{k_1+1,r}(\Omega_h)}, & R_2 &= |\sigma(\cdot, \epsilon(\mathbf{u}))|_{W^{k_2,r'}(\Omega_h)}, \\
R_3 &= |f|_{W^{k_3+1,r'}(\Omega_h)}, & R_4 &= \|\epsilon(\mathbf{u})\|_{L^r(\Omega_h)}.
\end{aligned} \quad (52)$$
Proof. We set $\xi_h = u_h - u_I$, and as usual, $v^\perp = (I - \Pi_k^0)v$ for any $v \in L^2(\Omega)$. First, by a triangle inequality and Lemma 18 (with $s = k_1 + 1$), we have

$$ \|u - u_h\|_{\delta,r} \le \|u - u_I\|_{\delta,r} + \|\xi_h\|_{\delta,r} \lesssim (\delta + h^{k_1} R_1)^{\frac{r-2}{r}} h^{\frac{2k_1}{r}} R_1^{\frac{2}{r}} + \|\xi_h\|_{\delta,r}. \quad (53) $$

Since $\xi_h \in Z_h$, manipulating (26) and (1) and recalling (24) we have

$$ \begin{aligned} a_h(u_h, \xi_h) - a_h(u_I, \xi_h) &= \int_{\Omega} (-\nabla \cdot \sigma(\cdot, \epsilon(u)) - f) \cdot \xi_h - a_h(u_I, \xi_h) + (f_h, \xi_h) \\ &= \int_{\Omega} \left( \sigma(\cdot, \epsilon(u)) : \epsilon(\xi_h) - \sigma(\cdot, \Pi_{k-1}^0 \epsilon(u_I)) : \Pi_{k-1}^0 \epsilon(\xi_h) \right) - S(u_I^\perp, \xi_h^\perp) + (f_h - f, \xi_h) \\ &=: T_1 + T_2 + T_3. \end{aligned} \quad (54) $$

We next estimate each term on the right-hand side above. In the following $C$ will denote a generic positive constant independent of $h$ that may change at each occurrence, whereas the positive parameter $\theta$ adopted in (58) and (65) will be specified later.

• Estimate of $T_1$. Employing the definition of $L^2$-projection (13) we have

$$ \begin{aligned} T_1 &= \int_{\Omega} (\sigma(\cdot, \epsilon(u)) - \Pi_{k-1}^0 \sigma(\cdot, \Pi_{k-1}^0 \epsilon(u_I))) : \epsilon(\xi_h) \\ &= \int_{\Omega} ((I - \Pi_{k-1}^0)\sigma(\cdot, \epsilon(u))) : ((I - \Pi_{k-1}^0)\epsilon(\xi_h)) \\ &\quad + \int_{\Omega} (\sigma(\cdot, \epsilon(u)) - \sigma(\cdot, \Pi_{k-1}^0 \epsilon(u_I))) : \Pi_{k-1}^0 \epsilon(\xi_h) \\ &=: T_1^A + T_1^B. \end{aligned} \quad (55) $$

We now recall the standard polynomial interpolation result

$$ \|(I - \Pi_{k-1}^{0,E})\sigma(\cdot, \epsilon(u))\|_{L^2(E)} \lesssim h_E^{k_2-1} |\sigma(\cdot, \epsilon(u))|_{\mathbb{W}^{k_2,1}(E)}. \quad (56) $$

Furthermore, combining Lemma 7 with the first line of equation (33), it is easy to check

$$ \|\epsilon(\xi_h^\perp)\|_{L^2(E)} \lesssim |E|^{\frac{1}{2}-\frac{1}{r}} S^E(\xi_h^\perp, \xi_h^\perp)^{\frac{1}{r}}. \quad (57) $$

The term $T_1^A$ can be bounded as follows

$$ \begin{aligned} T_1^A &= \sum_{E \in \Omega_h} \int_E ((I - \Pi_{k-1}^{0,E})\sigma(\cdot, \epsilon(u))) : ((I - \Pi_{k-1}^{0,E})\epsilon(\xi_h)) \\ &\le \sum_{E \in \Omega_h} \|(I - \Pi_{k-1}^{0,E})\sigma(\cdot, \epsilon(u))\|_{L^2(E)} \|(I - \Pi_{k-1}^{0,E})\epsilon(\xi_h)\|_{L^2(E)} \quad (\text{Cauchy-Schwarz ineq.}) \\ &\le C \sum_{E \in \Omega_h} h_E^{k_2-1} |\sigma(\cdot, \epsilon(u))|_{\mathbb{W}^{k_2,1}(E)} \|\epsilon(\xi_h^\perp)\|_{L^2(E)} \quad (\text{Cont of } \Pi_{k-1}^{0,E} \quad (56)) \\ &\le C \sum_{E \in \Omega_h} h_E^{k_2-1} |E|^{\frac{1}{r}} |\sigma(\cdot, \epsilon(u))|_{\mathbb{W}^{k_2,r'}(E)} |E|^{\frac{1}{2}-\frac{1}{r}} S^E(\xi_h^\perp, \xi_h^\perp)^{\frac{1}{r}} \quad ((r', r)\text{-Hölder ineq. } \& (57)) \\ &\le \frac{C}{r'\theta^{r'}} h^{k_2 r'} |\sigma(\cdot, \epsilon(u))|_{\mathbb{W}^{k_2,r'}(\Omega_h)}^{r'} + \frac{\theta^r}{2r} S(\xi_h^\perp, \xi_h^\perp), \quad ((r', r)\text{-Young ineq.}) \end{aligned} \quad (58) $$

where we used the fact that $|E|^{\frac{1}{2}} \approx h_E$. Employing (46a), we obtain

$$ \begin{aligned} T_1^B &\le \sum_{E \in \Omega_h} \int_E |\sigma(\cdot, \epsilon(u)) - \sigma(\cdot, \Pi_{k-1}^{0,E} \epsilon(u_I))| |\Pi_{k-1}^{0,E} \epsilon(\xi_h)| \\ &\le \sum_{E \in \Omega_h} \int_E \varphi'_{\delta+|\Pi_{k-1}^{0,E} \epsilon(u_I)|} (|\Pi_{k-1}^{0,E} \epsilon(u_I) - \epsilon(u)|) |\Pi_{k-1}^{0,E} \epsilon(\xi_h)|. \end{aligned} $$

Employing (42a) we get that for every $\varepsilon > 0$, there exists a positive constant $C_\varepsilon$ such that

$$ T_1^B \le \varepsilon \sum_{E \in \Omega_h} \int_E \varphi_{\delta+|\Pi_{k-1}^{0,E} \epsilon(u_I)|} (|\Pi_{k-1}^{0,E} \epsilon(\xi_h)|) + C_\varepsilon \sum_{E \in \Omega_h} \int_E \varphi_{\delta+|\Pi_{k-1}^{0,E} \epsilon(u_I)|} (|\Pi_{k-1}^{0,E} \epsilon(u_I) - \epsilon(u)|). $$
Using (46b) we obtain (with $\gamma$ denoting the associated uniform hidden constant)

$$ T_1^B \le \gamma\epsilon(\sigma(\cdot, \mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{u}_h)) - \sigma(\cdot, \mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{u}_I)), \mathbf{\Pi}_{k-1}^0 \epsilon(\xi_h)) \\ + C_\epsilon \sum_{E \in \Omega_h} \int_E (\delta + |\mathbf{\Pi}_{k-1}^{0,E} \epsilon(\mathbf{u}_I)| + |\epsilon(\mathbf{u})|)^{r-2} |\epsilon(\mathbf{u}) - \mathbf{\Pi}_{k-1}^{0,E} \epsilon(\mathbf{u}_I)|^2. $$

Notice that the constant $C_\epsilon$ depends only on $\sigma_c, \sigma_m, r$, and $\epsilon$. With respect to $\epsilon$, it may depend on the degree $k$ and the domain $\Omega$. However, given our mesh assumptions, it is independent of the particular mesh or mesh element within the family $\{\Omega_h\}_h$. Using an $(\frac{r}{r-2}, \frac{r}{2})$-Hölder inequality since $r > 2$, from the last equation we get

$$ T_1^B \le \gamma\epsilon(\sigma(\cdot, \mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{u}_h)) - \sigma(\cdot, \mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{u}_I)), \mathbf{\Pi}_{k-1}^0 \epsilon(\xi_h)) \\ + C_\epsilon \sum_{E \in \Omega_h} (|E|\delta^r + \|\mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{u}_I)\|_{\mathbb{L}^r(E)}^r + \|\epsilon(\mathbf{u})\|_{\mathbb{L}^r(E)}^r)^{\frac{r-2}{r}} \|\epsilon(\mathbf{u}) - \mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{u}_I)\|_{\mathbb{L}^r(E)}^2, $$

which, making use of Lemma 17 and taking $\epsilon = \frac{1}{2\gamma}$ becomes (using a discrete $(\frac{r}{r-2}, \frac{r}{2})$-Hölder inequality)

$$ T_1^B \le \frac{1}{2}(\sigma(\cdot, \mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{u}_h)) - \sigma(\cdot, \mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{u}_I)), \mathbf{\Pi}_{k-1}^0 \epsilon(\xi_h)) \\ + Ch^{2k_1} (|\Omega|\delta^r + \|\epsilon(\mathbf{u})\|_{\mathbb{L}^r(\Omega_h)}^r)^{\frac{r-2}{r}} |\mathbf{u}|_{W^{k_1+1,r}(\Omega_h)}^2. \quad (59) $$

Combining (58) and (59) in (55) we infer

$$ T_1 \le \frac{1}{2}(\sigma(\cdot, \mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{u}_h)) - \sigma(\cdot, \mathbf{\Pi}_{k-1}^0 \epsilon(\mathbf{u}_I)), \mathbf{\Pi}_{k-1}^0 \epsilon(\xi_h)) \\ + Ch^{2k_1}(\delta^r + R_4^r)^{\frac{r-2}{r}} R_1^2 + \frac{C}{r'\theta^{r'}} h^{k_2 r'} R_2^{r'} + \frac{\theta^r}{2r} \|\xi_h\|_{\delta,r}^r. \quad (60) $$

• Estimate of $T_2$. Recalling definitions (21) and (22), with some trivial algebra we obtain

$$ T_2 = - \sum_{E \in \Omega_h} S^E(\mathbf{u}_I^{\perp}, \xi_h^{\perp}) \le \sum_{E \in \Omega_h} Ch_E^{2-r}(\delta h_E + |\chi(\mathbf{u}_I^{\perp})|)^{r-2} |\chi(\mathbf{u}_I^{\perp})| |\chi(\xi_h^{\perp})|. $$

Employing (46a) in Lemma 16 to the scalar function $\sigma(\tau) = (\delta h_E + |\tau|)^{r-2}\tau$ (hence $\eta = 0$) we have

$$ T_2 \le C \sum_{E \in \Omega_h} h_E^{2-r} \varphi'_{\delta h_E + |\chi(\mathbf{u}_I^{\perp})|} (|\chi(\mathbf{u}_I^{\perp})|) |\chi(\xi_h^{\perp})|. $$

Employing (42a) we get that for every $\epsilon > 0$ there exists a positive constant $C_\epsilon$ such that

$$ T_2 \le \epsilon \sum_{E \in \Omega_h} h_E^{2-r} \varphi_{\delta h_E + |\chi(\mathbf{u}_I^{\perp})|} (|\chi(\xi_h^{\perp})|) + C_\epsilon \sum_{E \in \Omega_h} h_E^{2-r} \varphi_{\delta h_E + |\chi(\mathbf{u}_I^{\perp})|} (|\chi(\mathbf{u}_I^{\perp})|). $$

We now use (46b) and, denoting with $\gamma$ the hidden constant, we infer

$$ T_2 \le \epsilon\gamma \sum_{E \in \Omega_h} h_E^{2-r} (\delta h_E + |\chi(\mathbf{u}_I^{\perp})| + |\chi(\mathbf{u}_h^{\perp})|)^{r-2} |\chi(\xi_h^{\perp})|^2 \\ + C_\epsilon \sum_{E \in \Omega_h} h_E^{2-r} (\delta h_E + |\chi(\mathbf{u}_I^{\perp})|)^{r-2} |\chi(\mathbf{u}_I^{\perp})|^2 =: T_2^A + T_2^B. \quad (61) $$

Employing (33) (cf. proof of Lemma 9) with $\mathbf{u}_h = \mathbf{u}_h^{\perp}$ and $\mathbf{w}_h = \mathbf{u}_I^{\perp}$ and recalling definition (21) we obtain

$$ T_2^A \le \epsilon C \gamma \sum_{E \in \Omega_h} (S^E(\mathbf{u}_h^{\perp}, \xi_h^{\perp}) - S^E(\mathbf{u}_I^{\perp}, \xi_h^{\perp})) = \epsilon C \gamma (S(\mathbf{u}_h^{\perp}, \xi_h^{\perp}) - S(\mathbf{u}_I^{\perp}, \xi_h^{\perp})). \quad (62) $$

Reasoning as in (50) and employing Lemma 17 we can easily derive

$$ T_2^B \le h^{2k_1}(\delta + h^{k_1} R_1)^{r-2} R_1^2. \quad (63) $$

Combining (62) and (63) in (61) and taking $\epsilon = \frac{1}{2C\gamma}$ we infer

$$ T_2 \le \frac{1}{2}(S(\mathbf{u}_h^{\perp}, \xi_h^{\perp}) - S(\mathbf{u}_I^{\perp}, \xi_h^{\perp})) + Ch^{2k_1}(\delta + h^{k_1} R_1)^{r-2} R_1^2. \quad (64) $$
• Estimate of $T_3$. With analogous arguments to those in (58) we infer

$$
\begin{align*}
T_3 &= \sum_{E \in \Omega_h} \int_E (\Pi_k^{0,E} f - f) \cdot \xi_h^{\perp} & \text{(def. (20) & def. (13))} \\
&\le C \sum_{E \in \Omega_h} h_E^{k_3} |f|_{\overline{W}^{k_3+1,1}(E)} h_E |\chi(\xi_h^{\perp})| & \text{(same reasoning as in (58))} \\
&\le C \sum_{E \in \Omega_h} h_E^{k_3} |E|^{\frac{1}{r}} |f|_{\overline{W}^{k_3+1,r'}(E)} h_E^{\frac{1+r-2}{r}} S^E(\xi_h^{\perp}, \xi_h^{\perp})^{\frac{1}{r}} & \text{(r', r)-Hölder ineq. & Lemma 9} \tag{65} \\
&\le \frac{C}{r'\theta r'} h^{(k_3+2)r'} |f|_{\overline{W}^{k_3+1,r'}(\Omega_h)}^{r'} + \frac{\theta^r}{2r} S(\xi_h^{\perp}, \xi_h^{\perp}) & \text{($|E|^{\frac{1}{2}} \approx h_E$ & $\&$ ($r'$, $r$)-Young ineq.)} \\
&\le \frac{C}{r'\theta r'} h^{(k_3+2)r'} R_3^{r'} + \frac{\theta^r}{2r} \|\xi_h\|_{\delta,r}^r &
\end{align*}
$$

Plugging the estimates in (60), (64) and (65) in (54) and recalling definition (24) we obtain

$$
\begin{align}
\frac{1}{2}(a_h(\mathbf{u}_h, \boldsymbol{\xi}_h) - a_h(\mathbf{u}_I, \boldsymbol{\xi}_h)) \le C h^{2k_1} (\delta + h^{k_1} R_1 + R_4)^{r-2} R_1^2 &+ \frac{C}{r'\theta r'} h^{k_2 r'} R_2^{r'} \nonumber \\
&+ \frac{C}{r'\theta r'} h^{(k_3+2)r'} R_3^{r'} + \frac{\theta^r}{r} \|\boldsymbol{\xi}_h\|_{\delta,r}^r. \tag{66}
\end{align}
$$

We now write with (37) that

$$
\frac{1}{2} (a_h(\mathbf{u}_h, \boldsymbol{\xi}_h) - a_h(\mathbf{u}_I, \boldsymbol{\xi}_h)) \geq \tilde{C} \|\boldsymbol{\xi}_h\|_{\delta,r}^r. \quad (67)
$$

Combining (66) and (67) and taking $\theta = (\frac{\bar{c}r}{2})^{1/r}$ we obtain:

$$
\|\xi_h\|_{\delta,r}^r \le h^{2k_1} (\delta + h^{k_1} R_1 + R_4)^{r-2} R_1^2 + h^{k_2 r'} R_2^{r'} + h^{(k_3+2) r'} R_3^{r'}. \quad (68)
$$

The proof follows by (53) together with the previous bound.

□

According to Remark 5, if in addition to Assumption 2 we assume that all the mesh elements are convex, we can state the error estimate with respect to the standard $W^{1,r}$-norm. Indeed, in this case, due to the elliptic regularity results established in [46] together with Sobolev embeddings, we have $U_h(E) \subset W^{2,2}(E) \subset W^{1,r}(E)$ for all $E \in \Omega_h$ and $r \in [2, \infty)$, which combined with the global continuity of $U_h$ implies $U_h \subset W^{1,r}(\Omega)$. Furthermore, in such case our error estimate in Theorem 19, which is expressed in the $\| \cdot \|_{\delta,r}$ norm, directly translates into equivalent estimates in the classical $W^{1,r}(\Omega)$ norm.

**Corollary 20.** Let **u** be the solution of problem (7) and let **u**<sub>h</sub> be the solution of problem (26). Let all the mesh elements be convex. Under the regularity assumptions of Theorem 19 and with the regularity terms defined as in (52), it holds

$$
\|\mathbf{u} - \mathbf{u}_h\|_{W^{1,r}(\Omega)} \lesssim (\delta + h^{k_1} R_1 + R_4)^{\frac{r-2}{r}} h^{\frac{2k_1}{r}} R_1^{\frac{2}{r}} + h^{\frac{k_2}{r-1}} R_2^{\frac{1}{r-1}} + h^{\frac{k_3+2}{r-1}} R_3^{\frac{1}{r-1}}.
$$

*Proof.* We show the proof only briefly. The key step is deriving that

$$
\|\mathbf{v}_h\|_{W^{1,r}(\Omega)} \lesssim \|\mathbf{v}_h\|_{\mathbf{r}} \quad \forall \mathbf{v}_h \in \mathcal{U}_h, \tag{69}
$$

which can be shown easily, first by the Korn and triangle inequalities,

$$
\begin{align*}
\|\mathbf{v}_h\|_{W^{1,r}(\Omega)} &\lesssim \|\Pi_{k-1}^0 \epsilon(\mathbf{v}_h)\|_{L^r(\Omega_h)} + \|(I - \Pi_{k-1}^0)\epsilon(\mathbf{v}_h)\|_{L^r(\Omega_h)} \\
&\lesssim \|\Pi_{k-1}^0 \epsilon(\mathbf{v}_h)\|_{L^r(\Omega_h)} + |(I - \Pi_k^0)\mathbf{v}_h|_{W^{1,r}(\Omega_h)}
\end{align*}
$$

and then by definition of $\|\cdot\|_r$ norm and the stabilizing form, using Lemma 7 and an inverse estimate on
virtual functions. Afterwards, by the triangle inequality and interpolation estimates combined with (69), we
can write

$$
\|\mathbf{u} - \mathbf{u}_h\|_{W^{1,r}(\Omega)} \lesssim \|\mathbf{u} - \mathbf{u}_I\|_{W^{1,r}(\Omega)} + \|\mathbf{u}_h - \mathbf{u}_I\|_{W^{1,r}(\Omega)} \lesssim C h^k |\mathbf{u}|_{W^{k+1,r}(\Omega)} + \|\mathbf{u}_h - \mathbf{u}_I\|_r.
$$

The proof then follows from the above bound and (29), combined with Theorem 19.

□
We have also the following Corollary, stating a better bound for the $\sigma$ approximation term in Theorem 19, valid in the $\delta > 0$ case.

**Corollary 21.** Under the same assumptions of Theorem 19, if in addition $\delta > 0$, $\sigma(\cdot, \epsilon(\mathbf{u})) \in \mathbb{W}^{k_2,2}(\Omega_h)$ and $f \in \mathbf{W}^{k_3+1,2}(\Omega_h)$, it holds

$$ \|\mathbf{u} - \mathbf{u}_h\|_{\delta,r} \le h^{\frac{2k_1}{r}} (\delta + h^{k_1} R_1 + R_4)^{\frac{r-2}{r}} R_1^{\frac{2}{r}} + h^{\frac{2k_2}{r}} \delta^{\frac{2-r}{r}} \hat{R}_2^{\frac{2}{r}} + h^{\frac{2(k_3+2)}{r}} \delta^{\frac{2-r}{r}} \hat{R}_3^{\frac{2}{r}}, \quad (70) $$

with new regularity terms

$$ \hat{R}_2 = |\sigma(\cdot, \epsilon(\mathbf{u}))|_{\mathbb{W}^{k_2,2}(\Omega_h)}, \quad \hat{R}_3 := |f|_{\mathbf{W}^{k_3+1,2}(\Omega_h)}. \quad (71) $$

*Proof.* The proof follows the same steps as that of Theorem 19, the only modification being the bounds for $T_1^A$ and for $T_3$. We start by noting that from definition (22) and Lemma 7 it follows

$$ S^E(\xi_h^{\perp}, \xi_h^{\perp}) \ge \delta^{r-2} |\chi(\xi_h^{\perp})|^2 \ge \delta^{r-2} |\xi_h^{\perp}|^2_{W^{1,2}(E)}. \quad (72) $$

Furthermore, standard polynomial interpolation results yield

$$ \|(I - \Pi_{k-1}^{0,E})\sigma(\cdot, \epsilon(\mathbf{u}))\|_{L^2(E)} \le C h_E^{k_2} |\sigma(\cdot, \epsilon(\mathbf{u}))|_{\mathbb{W}^{k_2,2}(E)}. \quad (73) $$

We now take the steps from (58) but modify the bound for $\|\epsilon(\xi_h^{\perp})\|_{L^2(E)}$ using (72) instead of (57), and apply (73) instead of (56). We obtain, also using a classical Young inequality,

$$ \begin{aligned} T_1^A &\le C \delta^{(2-r)/2} \sum_{E \in \Omega_h} h_E^{k_2} |\sigma(\cdot, \epsilon(\mathbf{u}))|_{\mathbb{W}^{k_2,2}(E)} S^E(\xi_h^{\perp}, \xi_h^{\perp})^{1/2} \\ &\le \delta^{2-r} \frac{C}{2\theta^r} h^{2k_2} |\sigma(\cdot, \epsilon(\mathbf{u}))|^2_{\mathbb{W}^{k_2,2}(\Omega_h)} + \frac{\theta^r}{2r} S(\xi_h^{\perp}, \xi_h^{\perp}). \end{aligned} \quad (74) $$

The loading term $T_3$ can be handled by analogous modifications. The rest of the proof then follows as in Theorem 19. $\square$

## 4.3 A priori error estimate: pressure

**Theorem 22.** Let $(\mathbf{u}, p)$ be the solution of problem (5) and $(\mathbf{u}_h, p_h)$ the solution of problem (25). Assume that $\mathbf{u} \in \mathbf{W}^{k_1+1,r}(\Omega_h)$, $\sigma(\cdot, \epsilon(\mathbf{u})) \in \mathbb{W}^{k_2,r'}(\Omega_h)$, $f \in \mathbf{W}^{k_3+1,r'}(\Omega_h)$, and $p \in \mathbf{W}^{k_4,r'}(\Omega_h)$, for some $k_1, k_2, k_3, k_4 \le k$. Let the mesh regularity in Assumption 2 hold. Then, we have

$$ \begin{aligned} \|p - p_h\|_{L^{r'}(\Omega_h)} &\le h^{k_2} R_2 + h^{k_3+2} R_3 + (\delta + \|\mathbf{u} - \mathbf{u}_h\|_{\delta,r} + R_4)^{r-2} (\|\mathbf{u} - \mathbf{u}_h\|_{\delta,r} + h^{k_1} R_1) \\ &\quad + (\delta^r + h^{2k_1} (\delta + h^{k_1} R_1 + R_4)^{r-2} R_1^2 + h^{k_2 r'} R_2^{r'} + h^{(k_3+2) r'} R_3^{r'})^{\frac{r-2}{2r}} \\ &\quad \times (h^{2k_1} (\delta + h^{k_1} R_1 + R_4)^{r-2} R_1^2 + h^{k_2 r'} R_2^{r'} + h^{(k_3+2) r'} R_3^{r'})^{\frac{1}{2}} + h^{k_4} R_5, \end{aligned} \quad (75) $$

where the regularity term $R_1, R_2, R_3$, and $R_4$ are defined in (52) and $R_5 = |p|_{\mathbf{W}^{k_4,r'}(\Omega_h)}$.

*Proof.* Let $p_I = \Pi_{k-1}^0 p \in \mathbb{P}_{k-1}(\Omega_h)$ and let $\rho_h = p_h - p_I$. For the sake of brevity also in this proof we employ the notation $\mathbf{v}^{\perp} = (I - \Pi_k^0)\mathbf{v}$ for any $\mathbf{v} \in L^2(\Omega)$. Employing (5) and (25), recalling the definition of the form $b(\cdot, \cdot)$ in (6) and combining item (ii) in (17) with the definition of $L^2$-projection, for all $\mathbf{w}_h \in \mathcal{U}_h$ we get

$$ \begin{aligned} b(\mathbf{w}_h, \rho_h) &= -a_h(\mathbf{u}_h, \mathbf{w}_h) + (f_h, \mathbf{w}_h) + a(\mathbf{u}, \mathbf{w}_h) - (f, \mathbf{w}_h) + b(\mathbf{w}_h, p - p_I) \\ &= -a_h(\mathbf{u}_h, \mathbf{w}_h) + a(\mathbf{u}, \mathbf{w}_h) + (f_h - f, \mathbf{w}_h) \\ &= \int_{\Omega} \left( (I - \Pi_{k-1}^0)\sigma(\cdot, \epsilon(\mathbf{u})) \right) : (I - \Pi_{k-1}^0)\epsilon(\mathbf{w}_h) + (f_h - f, \mathbf{w}_h) \\ &\quad + \int_{\Omega} \left( \sigma(\cdot, \Pi_{k-1}^0 \epsilon(\mathbf{u})) - \sigma(\cdot, \epsilon(\mathbf{u}_h)) \right) : \Pi_{k-1}^0 \epsilon(\mathbf{w}_h) - S(\mathbf{u}_h^{\perp}, \mathbf{w}_h^{\perp}) \\ &=: T_1 + T_2 + T_3 + T_4. \end{aligned} \quad (76) $$

In the first identity above, we have taken $\mathbf{v} = \mathbf{w}_h$ as a test function in the continuous weak problem (5) even if, as observed in Remark 5, $\mathbf{w}_h \in \mathbf{W}^{1,2}(\Omega)$ but $\mathbf{w}_h \notin \mathbf{W}^{1,r}(\Omega)$. This is possible due to the additional regularity assumptions on the forcing term $f \in \mathbf{W}^{1,r'}(\Omega_h) \subset L^2(\Omega)$ and exact stress field
$\sigma(\cdot, \epsilon(\mathbf{u})) \in \mathbb{W}^{1,r'}(\Omega_h) \subset \mathbb{L}^2(\Omega)$. We estimate separately each term in (76). Employing the Hölder inequality with exponents $(r', r)$ and polynomial approximation properties we infer

$$
\begin{align}
T_1 &\le \sum_{E \in \Omega_h} \|(I - \Pi_{k-1}^{0,E}) \sigma(\cdot, \epsilon(\mathbf{u}))\|_{\mathbb{L}^2(E)} \|(I - \Pi_{k-1}^{0,E}) \epsilon(\mathbf{w}_h)\|_{\mathbb{L}^2(E)} \quad \text{(Cauchy-Schwarz ineq.)} \\
&\le \sum_{E \in \Omega_h} h_E^{k_2-1} |\sigma(\cdot, \epsilon(\mathbf{u}))|_{\mathbb{W}^{k_2,1}(E)} \|\epsilon(\mathbf{w}_h^\perp)\|_{\mathbb{L}^2(E)} \quad \text{(Cont of $\Pi_{k-1}^{0,E}$ & (56))} \\
&\le \sum_{E \in \Omega_h} h_E^{k_2-1} |E|^{\frac{1}{r}} |\sigma(\cdot, \epsilon(\mathbf{u}))|_{\mathbb{W}^{k_2,r'}(E)} |E|^{\frac{1}{2}-\frac{1}{r}} S^E(\mathbf{w}_h^\perp, \mathbf{w}_h^\perp)^{\frac{1}{r}} \quad \text{((r', r)-Hölder ineq. & (57))} \\
&\le h^{k_2} |\sigma(\cdot, \epsilon(\mathbf{u}))|_{\mathbb{W}^{k_2,r'}(\Omega_h)} S(\mathbf{w}_h^\perp, \mathbf{w}_h^\perp)^{\frac{1}{r}} \quad \text{($|E|^{\frac{1}{2}} \approx h_E$ & $(r', r)$-Young ineq.)} \\
&\le h^{k_2} R_2 (\delta + \|\mathbf{w}_h\|_r)^{\frac{2}{r}} \|\mathbf{w}_h\|_r^{\frac{2}{r}}. \quad \text{(apply (29))} \tag{77}
\end{align}
$$

On the other hand,

$$
\begin{align}
T_2 &= \sum_{E \in \Omega_h} \int_E (\Pi_k^{0,E} f - f) \cdot w_h^\perp \tag{def. (20) & def. (13)} \\
&\le \sum_{E \in \Omega_h} \| (\Pi_k^{0,E} f - f) \|_{\mathbb{L}^2(E)} \| w_h^\perp \|_{\mathbb{L}^2(E)} \le C \sum_{E \in \Omega_h} h_E^{k_3} |f|_{\mathbb{W}^{k_3+1,r'}(E)} h_E |\chi(w_h^\perp)| \tag{same as in (65)} \\
&\le \sum_{E \in \Omega_h} h_E^{k_3} |E|^{\frac{1}{r}} |f|_{\mathbb{W}^{k_3+1,r'}(E)} h_E^{1+\frac{r-2}{r}} S^E(w_h^\perp, w_h^\perp)^{\frac{1}{r}} \tag{(r', r)-Hölder ineq. & Lemma 9} \\
&\le h^{k_3+2} |f|_{\mathbb{W}^{k_3+1,r'}(\Omega_h)} S(w_h^\perp, w_h^\perp)^{\frac{1}{r}} \tag{$|E|^{\frac{1}{2}} \approx h_E$ & $(r', r)$-Young ineq.) \\
&\le h^{k_3+2} R_3 (\delta + \|w_h\|_r)^{\frac{2}{r}} \|w_h\|_r^{\frac{2}{r}}. \tag{apply (29)} \tag{78}
\end{align}
$$

Employing (3a), the 3-terms ($\frac{r}{r-2}$, $r, r$)-Hölder inequality, the triangle inequality together with the $L^r$-stability of $\Pi_{k-1}^0$ we have

$$
\begin{align}
T_3 &\lesssim \int_{\Omega} (\delta + |\Pi_{k-1}^0 \epsilon(\mathbf{u}_h)| + |\epsilon(\mathbf{u})|)^{r-2} |\Pi_{k-1}^0 \epsilon(\mathbf{u}_h) - \epsilon(\mathbf{u})| |\Pi_{k-1}^0 \epsilon(\mathbf{w}_h)| \nonumber \\
&\lesssim (\delta^r + \|\Pi_{k-1}^0 \epsilon(\mathbf{u}_h)\|_{\mathbb{L}^r(\Omega)}^r + \|\epsilon(\mathbf{u})\|_{\mathbb{L}^r(\Omega)}^r)^{\frac{r-2}{r}} \|\Pi_{k-1}^0 \epsilon(\mathbf{u}_h) - \epsilon(\mathbf{u})\|_{\mathbb{L}^r(\Omega)} \|\Pi_{k-1}^0 \epsilon(\mathbf{w}_h)\|_{\mathbb{L}^r(\Omega)} \nonumber \\
&\lesssim (\delta^r + \|\mathbf{u} - \mathbf{u}_h\|_{\delta,r}^r + \|\epsilon(\mathbf{u})\|_{\mathbb{L}^r(\Omega)}^r)^{\frac{r-2}{r}} (\|\mathbf{u} - \mathbf{u}_h\|_{\delta,r} + \|(I - \Pi_{k-1}^0)\epsilon(\mathbf{u})\|_{\mathbb{L}^r(\Omega)} \|\mathbf{w}_h\|_r \nonumber \\
&\lesssim (\delta^r + \|\mathbf{u} - \mathbf{u}_h\|_{\delta,r}^r + R_4^r)^{\frac{r-2}{r}} (\|\mathbf{u} - \mathbf{u}_h\|_{\delta,r} + h^{k_1} R_1) (\delta + \|\mathbf{w}_h\|_r)^{\frac{r-2}{r}} \|\mathbf{w}_h\|_r^{\frac{2}{r}} \nonumber \\
&=: R_{1,h} (\delta + \|\mathbf{w}_h\|_r)^{\frac{r-2}{r}} \|\mathbf{w}_h\|_r^{\frac{2}{r}}. \tag{79}
\end{align}
$$

Using (34) and Lemma 24, we obtain

$$
\begin{align}
T_4 &\lesssim (\delta^r + S(\mathbf{u}_h^\perp, \mathbf{u}_h^\perp))^{\frac{r-2}{2r}} S(\mathbf{u}_h^\perp, \mathbf{u}_h^\perp)^{\frac{1}{2}} S(\mathbf{w}_h^\perp, \mathbf{w}_h^\perp)^{\frac{1}{r}} \nonumber \\
&\lesssim (\delta^r + S(\mathbf{u}_h^\perp, \mathbf{u}_h^\perp))^{\frac{r-2}{2r}} S(\mathbf{u}_h^\perp, \mathbf{u}_h^\perp)^{\frac{1}{2}} \|\mathbf{w}_h\|_{\delta,r} \nonumber \\
&\lesssim (\delta^r + S(\mathbf{u}_h^\perp, \boldsymbol{\xi}_h^\perp) - S(\mathbf{u}_I^\perp, \boldsymbol{\xi}_h^\perp) + R_h)^{\frac{r-2}{2r}} (S(\mathbf{u}_h^\perp, \boldsymbol{\xi}_h^\perp) - S(\mathbf{u}_I^\perp, \boldsymbol{\xi}_h^\perp) + R_h)^{\frac{1}{2}} \|\mathbf{w}_h\|_{\delta,r} \nonumber \\
&\lesssim (\delta^r + a_h(\mathbf{u}_h, \boldsymbol{\xi}_h) - a_h(\mathbf{u}_I, \boldsymbol{\xi}_h) + R_h)^{\frac{r-2}{2r}} (a_h(\mathbf{u}_h, \boldsymbol{\xi}_h) - a_h(\mathbf{u}_I, \boldsymbol{\xi}_h) + R_h)^{\frac{1}{2}} \|\mathbf{w}_h\|_{\delta,r}, \tag{80}
\end{align}
$$

which, using the bound for $(a_h(\mathbf{u}_h, \boldsymbol{\xi}_h) - a_h(\mathbf{u}_I, \boldsymbol{\xi}_h))$ encased in the proof of Theorem 19, check equation (54), becomes

$$
\begin{align}
T_4 &\lesssim (\delta^r + R'_h + h^{k_2 r'} R_2^{r'} + h^{(k_3+2) r'} R_3^{r'})^{\frac{r-2}{2r}} (R_h + h^{k_2 r'} R_2^{r'} + h^{(k_3+2) r'} R_3^{r'})^{\frac{1}{2}} (\delta + \|\mathbf{w}_h\|_r)^{\frac{r-2}{r}} \|\mathbf{w}_h\|_r^{\frac{2}{r}} \tag{81} \\
&:= R_{2,h} (\delta + \|\mathbf{w}_h\|_r)^{\frac{r-2}{r}} \|\mathbf{w}_h\|_r^{\frac{2}{r}}, \nonumber
\end{align}
$$

where $R'_h = h^{2k_1} (\delta + h^{k_1} R_1 + R_4)^{r-2} R_1^2$. Employing the discrete inf-sup of Lemma 8, and collecting (77),
(78), (79), and (81) in (76), we obtain

$$
\begin{align}
\|\rho_h\|_{L^{r'}(\Omega_h)} &\le \frac{1}{\bar{\beta}(r)} \sup_{\mathbf{w}_h \in \mathcal{U}_h, \|\mathbf{w}_h\|_r=1} b(\mathbf{w}_h, \rho_h) \nonumber \\
&\lesssim \sup_{\mathbf{w}_h \in \mathcal{U}_h, \|\mathbf{w}_h\|_r=1} \left[ (h^{k_2} R_2 + h^{k_3+2} R_3 + R_{1,h} + R_{2,h}) (\delta + \|\mathbf{w}_h\|_r)^{\frac{r-2}{r}} \|\mathbf{w}_h\|_r^{\frac{2}{r}} \right] \tag{82} \\
&= (h^{k_2} R_2 + h^{k_3+2} R_3 + R_{1,h} + R_{2,h}) (\delta + 1)^{\frac{r-2}{r}}. \nonumber
\end{align}
$$

Now the thesis follows from triangular inequality and standard polynomial approximation properties

$$
\|p - p_h\|_{L^{r'}(\Omega_h)} \le \|p - p_I\|_{L^{r'}(\Omega_h)} + \|\rho_h\|_{L^{r'}(\Omega_h)} \le h^{k_4} R_5 + \|\rho_h\|_{L^{r'}(\Omega_h)}. \quad (83)
$$

<figure><img src="image_2.png" /></figure>

**Corollary 23.** Under the same assumptions of Theorem 19, if in addition $\delta > 0$, $\sigma(\cdot, \epsilon(\mathbf{u})) \in \mathbb{W}^{k_2,2}(\Omega_h)$ and $f \in \mathbb{W}^{k_3+1,2}(\Omega_h)$, it holds

$$
\begin{align}
\|p - p_h\|_{L^{r'}(\Omega_h)} &\lesssim h^{k_2} \delta^{\frac{2-r}{2}} \hat{R}_2 + h^{k_3+2} \delta^{\frac{2-r}{2}} \hat{R}_3 + (\delta + \|\mathbf{u} - \mathbf{u}_h\|_{\delta,r} + R_4)^{r-2} (\|\mathbf{u} - \mathbf{u}_h\|_{\delta,r} + h^{k_1} R_1) \nonumber \\
&\quad + (\delta^r + h^{2k_1} (\delta + h^{k_1} R_1 + R_4)^{r-2} R_1^2 + h^{2k_2} \delta^{\frac{2-r}{2}} \hat{R}_2^2 + h^{2(k_3+2)} \delta^{\frac{2-r}{2}} \hat{R}_3^2)^{\frac{r-2}{2r}} \nonumber \\
&\quad \times (h^{2k_1} (\delta + h^{k_1} R_1 + R_4)^{r-2} R_1^2 + h^{2k_2} \delta^{\frac{2-r}{2}} \hat{R}_2^2 + h^{2(k_3+2)} \delta^{\frac{2-r}{2}} \hat{R}_3^2)^{\frac{1}{2}} + h^{k_4} R_5, \tag{84}
\end{align}
$$

where the regularity terms $R_1$ and $R_4$ are from (52), $\hat{R}_2$ and $\hat{R}_3$ from (71), and $R_5 = |p|_{W^{k_4,r'}(\Omega_h)}$.

Proof. The proof follows the same steps as that of Theorem 22 but applying Corollary 21 instead of Theorem 19.
□

**Lemma 24.** Under the same assumptions of Theorem 19, let $\mathbf{u}_I \in \mathcal{U}_h$ be the interpolant of $\mathbf{u}$ (cf. Lemma 4) and $\xi_h = \mathbf{u}_h - \mathbf{u}_I$, then the following holds

$$
S(\mathbf{u}_h^{\perp}, \mathbf{u}_h^{\perp}) \lesssim S(\mathbf{u}_h^{\perp}, \xi_h^{\perp}) - S(\mathbf{u}_I^{\perp}, \xi_h^{\perp}) + h^{2k_1} (\delta + h^{k_1} R_1)^{r-2} R_1^2, \quad (85)
$$

where the regularity term $R_1$ is from (52).

*Proof.* Let $\hat{\sigma}(\mathbf{x}) = (h_E \delta + |\mathbf{x}|)^{r-2} \mathbf{x}$ for any $\mathbf{x} \in \mathbb{R}^{N_E}$. Then simple computations yield

$$
\begin{align}
S(\mathbf{u}_h^{\perp}, \mathbf{u}_h^{\perp}) &\approx \sum_{E \in \Omega_h} h_E^{2-r} (h_E \delta + |\chi(\mathbf{u}_h^{\perp})|)^{r-2} |\chi(\mathbf{u}_h^{\perp})|^2 = \sum_{E \in \Omega_h} h_E^{2-r} \hat{\sigma}(\chi(\mathbf{u}_h^{\perp})) \cdot \chi(\mathbf{u}_h^{\perp}) \nonumber \\
&= \sum_{E \in \Omega_h} h_E^{2-r} (\hat{\sigma}(\chi(\mathbf{u}_h^{\perp})) - \hat{\sigma}(\chi(\mathbf{u}_I^{\perp}))) \cdot \chi(\mathbf{u}_h^{\perp}) + \sum_{E \in \Omega_h} h_E^{2-r} \hat{\sigma}(\chi(\mathbf{u}_h^{\perp})) \cdot \chi(\mathbf{u}_h^{\perp}) \tag{86} \\
&=: T_1 + T_2. \nonumber
\end{align}
$$

In the following $C$ will denote a generic positive constant independent of $h$ that may change at each occurrence,
whereas the parameter $\varepsilon$ adopted in (88) and (89) will be specified later. Using Lemma 15 and Lemma 16
we infer

$$
T_1 \le C \sum_{E \in \Omega_h} h_E^{2-r} \varphi'_{h_E \delta + |\chi(\mathbf{u}_h^\perp)|} (|\chi(\xi_h^\perp)|) |\chi(\mathbf{u}_h^\perp)| \quad \text{(by (46a))}
$$

$$
\le \varepsilon \sum_{E \in \Omega_h} h_E^{2-r} \varphi_{h_E \delta + |\chi(\mathbf{u}_h^\perp)|} (|\chi(\mathbf{u}_h^\perp)|) + C_\varepsilon \sum_{E \in \Omega_h} h_E^{2-r} \varphi_{h_E \delta + |\chi(\mathbf{u}_h^\perp)|} (|\chi(\xi_h^\perp)|) \quad \text{(by (42a))}
$$

$$
\le \gamma \varepsilon \sum_{E \in \Omega_h} h_E^{2-r} (h_E \delta + |\chi(\mathbf{u}_h^\perp)|)^{r-2} |\chi(\mathbf{u}_h^\perp)|^2 + C_\varepsilon \sum_{E \in \Omega_h} h_E^{2-r} (h_E \delta + |\chi(\mathbf{u}_h^\perp)| + |\chi(\mathbf{u}_I^\perp)|)^{r-2} |\chi(\xi_h^\perp)|^2, \quad \text{(by (46b))} \tag{87}
$$

where in the last line $\gamma$ denotes the uniform hidden positive constant from Lemma 16. Applying a consequent bound of (33) (from the second row to the last row) yields

$$
\begin{align}
T_1 &\le \gamma \epsilon S(\mathbf{u}_h^\perp, \mathbf{u}_h^\perp) + C_\epsilon \sum_{E \in \Omega_h} h_E^{2-r} (h_E \delta + |\chi(\mathbf{u}_h^\perp)| + |\chi(\mathbf{u}_I^\perp)|)^{r-2} |\chi(\xi_h^\perp)|^2 \tag{88} \\
&\le \gamma \epsilon S(\mathbf{u}_h^\perp, \mathbf{u}_h^\perp) + C_\epsilon (S(\mathbf{u}_h^\perp, \xi_h^\perp) - S(\mathbf{u}_I^\perp, \xi_h^\perp)). \nonumber
\end{align}
$$
Using analogous arguments we have

$$ T_2 \le C \sum_{E \in \Omega_h} h_E^{2-r} \varphi'_{h_E \delta + |\chi(\mathbf{u}_I^\perp)|} (|\chi(\mathbf{u}_I^\perp)|) |\chi(\mathbf{u}_h^\perp)| \quad (\text{by (46a)}) $$

$$ \le \varepsilon \sum_{E \in \Omega_h} h_E^{2-r} \varphi_{h_E \delta} (|\chi(\mathbf{u}_h^\perp)|) + C_\varepsilon \sum_{E \in \Omega_h} h_E^{2-r} \varphi_{h_E \delta} (|\chi(\mathbf{u}_I^\perp)|) \quad (\text{by (42a)}) $$

$$ \le \gamma \varepsilon \sum_{E \in \Omega_h} h_E^{2-r} (h_E \delta + |\chi(\mathbf{u}_h^\perp)|)^{r-2} |\chi(\mathbf{u}_h^\perp)|^2 + C_\varepsilon \sum_{E \in \Omega_h} h_E^{2-r} (h_E \delta + |\chi(\mathbf{u}_I^\perp)|)^{r-2} |\chi(\mathbf{u}_I^\perp)|^2 \quad (\text{by (46b)}) $$

$$ \le \gamma \varepsilon S(\mathbf{u}_h^\perp, \mathbf{u}_h^\perp) + C_\varepsilon \sum_{E \in \Omega_h} h_E^{2-r} (h_E \delta + |\chi(\mathbf{u}_I^\perp)|)^{r-2} |\chi(\mathbf{u}_I^\perp)|^2, \quad (89) $$

Taking in (88) and (89) $\epsilon = \frac{1}{4\gamma}$, from (86) we obtain

$$ \begin{aligned} S(\mathbf{u}_h^\perp, \mathbf{u}_h^\perp) &\le S(\mathbf{u}_h^\perp, \xi_h^\perp) - S(\mathbf{u}_I^\perp, \xi_h^\perp) + \sum_{E \in \Omega_h} h_E^{2-r} (h_E \delta + |\chi(\mathbf{u}_I^\perp)|)^{r-2} |\chi(\mathbf{u}_I^\perp)|^2 \\ &\le S(\mathbf{u}_h^\perp, \xi_h^\perp) - S(\mathbf{u}_I^\perp, \xi_h^\perp) + h^{2k_1} (\delta + h^{k_1} R_1)^{r-2} R_1^2. \end{aligned} $$

where we used the bound (63) of $T_2^B$.

□

Remark 25 (Orders of convergence). We comment on the orders of convergence of velocity and pressure under enough regularity, i.e., if $k_1 = k_2 = k_3 = k$. Employing Theorem 19 and Corollary 21, for $h$ asymptotically small and for any $r \ge 2$ and $\delta \ge 0$, the order of convergence of the velocity is:

$$ O(u) = \min \left( \frac{k}{r-1}, \frac{2k}{r} \right) = \frac{k}{r-1} \quad \text{valid for } \delta \ge 0, \quad (90) $$

$$ O(u) = \frac{2k}{r} \quad \text{valid for } \delta > 0. $$

In the first bound above we leave the minimum expressed explicitly to shed light on the origin of the leading order $k/(r-1)$. Indeed, this order stems from the "$\sigma$ approximation" term $T_A^1$ (cf. proof of Theorem 19) and therefore, in many situations, it is expected to dominate the estimate only asymptotically, but possibly not for practical mesh sizes (see also Section 5). Furthermore note that such "$\sigma$-approximation term" could be ameliorated by raising the order of the projection $\Pi_{k-1}^0$ appearing in the first addendum (consistency part) of (24), that is using $\Pi_\ell^0$ with $\ell \ge k$. As a consequence, if $\sigma$ is sufficiently regular, a simple modification of bound (58) would lead to the more favorable final bound

$$ O(u) = \min \left( \frac{\ell + 1}{r-1}, \frac{2k}{r} \right) \quad \text{valid for } \delta \ge 0. $$

The above improvement can be achieved by suitably enhancing the virtual space, resulting in a more cumbersome computation of the local discrete forms (but not increasing the size of the global system).

For what concerns the pressure, from Theorem 22 and Corollary 23, for $h$ asymptotically small and for any $r \ge 2$ and $\delta \ge 0$, the order of convergence of the pressure is:

$$ O(p) = O(u). \quad (91) $$

Remark 26 (The role of $\delta$). When $\delta$ is positive but small and comparable to $h$, the pre-asymptotic error reduction rate will be the outcome of the competition among the two bounds appearing in (90). In such cases, although the asymptotic convergence rate will clearly behave as $h^{2k/r}$, we may experience a slower pre-asymptotic error reduction rate, more similar to $h^{k/r-1}$. More precisely, a careful analysis reveals that when $\delta$ is small and comparable to $h$, the term $\delta^{2-r}/r$ in (70) (and (84)) plays a role in the error reduction rate. In particular, we obtain for $h < 1$ that the velocity error is bounded by $Ch^\gamma$, with $C$ independent of $h$, $\delta$ and

$$ \gamma = \max \left( \frac{2k}{r} + \frac{2-r}{r} \frac{\ln(\min(\delta, 1))}{\ln(h)}, \frac{k}{r-1} \right) \quad \text{valid for } \delta > 0, \quad (92) $$

where we used $\delta^{2-r}/r = e^{2-r/r} \ln(\delta) = e^{(2-r/r) \ln(\delta)} \ln(h) = h^{2-r/r} \frac{\ln(\delta)}{\ln(h)}$ and we noted that if $\delta > 1$ then the other terms of (70) dominate the estimate, hence the minimum in (92).
# 5 Numerical Results

In this section, we present three numerical tests to validate the theoretical results of Theorems 19 and 22 (and the associated corollaries) for different values of the parameters $\delta$ and $r$, as well as of the Sobolev regularity indices $k_1, k_2, k_3$ and $k_4$. To compute the VEM error between the exact solution ($u_{ex}, p_{ex}$) and the VEM solution ($u_h, p_h$), we consider the computable error quantities

$$ 
\begin{aligned}
\mathrm{err}(u_h, \|\cdot\|_r) &= \frac{\|u_{ex} - u_h\|_r}{\|u_{ex}\|_r}, \\
\mathrm{err}(u_h, W^{1,r}) &= \frac{\|\nabla u_{ex} - \Pi_{k-1}^0 \nabla u_h\|_{L^r(\Omega)}}{\|\nabla u_{ex}\|_{L^r(\Omega)}}, \\
\mathrm{err}(p_h, L^{r'}) &= \frac{\|p_{ex} - p_h\|_{L^{r'}(\Omega)}}{\|p_{ex}\|_{L^{r'}(\Omega)}}, \\
\mathrm{err}(\sigma, L^{r'}) &= \frac{\|\sigma(\cdot, \epsilon(u_{ex})) - \sigma(\cdot, \Pi_{k-1}^0 \epsilon(u_h))\|_{L^{r'}(\Omega)}}{\|\sigma(\cdot, \epsilon(u_{ex}))\|_{L^{r'}(\Omega)}}.
\end{aligned}
\quad (93) $$

We make use of the $\|\cdot\|_r$ norm, which is bounded by the $\|\cdot\|_{\delta,r}$ norm, in order to have the same error measure in all tests. Furthermore, note that we also include an error measure on the stress $\sigma$, although deriving a theoretical estimate for such a quantity is beyond the scope of the present contribution.

Given a sequence of $N + 1$ meshes with mesh diameters $h_0 > \dots > h_N$, and denoting by $E_h$ any of the error quantities listed in (93), we define the average experimental order of convergence AEOC as

$$ \mathrm{AEOC} = \frac{1}{N} \sum_{n=1}^{N} \frac{\log(E_{h_{n-1}}/E_{h_n})}{\log(h_{n-1}/h_n)}. $$

As a model equation, we consider the Carreau-Yasuda model (2), with $\alpha = 2$ (i.e. corresponding to the Carreau model), $\mu = 1$. In order to verify the apriori error estimates of Section 4, numerical tests are performed with the following values of $r$ and $\delta$:

$$ r = 2.00, 2.25, 2.50, 3.00, \quad \delta = 1, 0. \quad (94) $$

In the forthcoming tests, we consider the scheme (25) with $k = 2$. The nonlinear problem is solved by means of a two-step Picard-type iteration. First, we solve the problem corresponding to $\bar{r}$, defined as the midpoint between 2 and $r$, using as initial guess the solution of the associated linear Stokes problem. Then, the solution obtained for $\bar{r}$ is employed as the initial iterate for a Picard iteration with exponent $r$. An analogous strategy was adopted in [5, Section 5.1] for the case $r \in (1, 2]$. The domain $\Omega$ (specified in each test) is partitioned with the following sequences of polygonal meshes: QUADRILATERAL distorted meshes, RANDOM Voronoi meshes, and CARTESIAN meshes (see Fig. 1). For the generation of the Voronoi meshes we used the code Polymesher [55]. We emphasize that, for the families of meshes under consideration, all mesh elements are convex, therefore, according to Remark 5, the discrete solution satisfies $u_h \in W^{1,r}$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Example of the adopted polygonal meshes.</figcaption></figure>

## 5.1 Test 1. Regular solution

In the first test case, we consider Problem (1) with full Dirichlet boundary conditions (i.e. $\Gamma_D = \partial\Omega$) on the unit square $\Omega = (0, 1)^2$. The load terms $f$ (depending on $r$ and $\delta$ in (2)) and the Dirichlet boundary
conditions are chosen according to the analytical solution

$$ \mathbf{u}_{\text{ex}}(x_1, x_2) = \left[ \begin{array}{l} \sin\left(\frac{\pi}{2}x_1\right) \cos\left(\frac{\pi}{2}x_2\right) \\ -\cos\left(\frac{\pi}{2}x_1\right) \sin\left(\frac{\pi}{2}x_2\right) \end{array} \right], \quad p_{\text{ex}}(x_1, x_2) = -\sin\left(\frac{\pi}{2}x_1\right) \sin\left(\frac{\pi}{2}x_2\right) + \frac{4}{\pi^2}. $$

The domain $\Omega$ is partitioned using the family of QUADRILATERAL distorted meshes and the family of RANDOM meshes. For each mesh family, we consider a mesh sequence with diameters $h = 1/4, 1/8, 1/16, 1/32, 1/64$. In Fig. 2 and Fig. 3, we plot the computed error quantities in (93) for the sequences of aforementioned meshes and parameters $r$ and $\delta$ as in (94). We observe that for $\delta = 1$ (left panels of Fig. 2 and Fig. 3), a convergence rate of order 2 is observed, whereas for $\delta = 0$ (left panels of Fig. 2 and Fig. 3), the plot shows the average experimental orders of convergence AEOC. In order to interpret the results illustrated in Fig. 2 and Fig. 3 with respect to the theoretical estimates established in Section 4, Table 1 reports the expected convergence orders corresponding to the different sources of error derived in Theorems 19 and 22, specifically the two terms appearing on the right in equations (53) and (83), respectively (see also Remark 25). In particular, we report the interpolation errors $\mathbf{u}_{\text{ex}} - \mathbf{u}_{I}$ and $p_{\text{ex}} - p_{I}$, as well as the terms $h^{k/(r-1)}$ and $h^{2k/r}$ appearing in the bounds (68) and (82). To analyze the stress errors $\text{err}(\sigma, L^{r'})$, we also show the quantities $\|\sigma(\cdot, \epsilon(\mathbf{u}_{\text{ex}})) - \sigma(\cdot, \epsilon(\mathbf{u}_{I}))\|_{L^{r'}}$ (denoted by $\|\sigma_{\text{ex}} - \sigma_{I}\|_{L^{r'}}$). It can be observed that, for $\delta = 1$, the inter-

<table><caption>Table 1: Test 1. Expected orders of convergence for the terms appearing in the *a priori* error estimates in Section 4.</caption><thead><tr><th>r</th><th>||<b>u</b><sub>ex</sub> - <b>u</b><sub>I</sub>||<sub>r</sub></th><th>||<b>u</b><sub>ex</sub> - <b>u</b><sub>I</sub>||<sub>W<sup>1,r</sup></sub></th><th>||<i>p</i><sub>ex</sub> - <i>p</i><sub>I</sub>||<sub>L<sup>r'</sup></sub></th><th>||<i>σ</i><sub>ex</sub> - <i>σ</i><sub>I</sub>||<sub>L<sup>r'</sup></sub></th><th>2/(<i>r</i> - 1)</th><th>4/<i>r</i></th></tr></thead><tbody><tr><td>2.00</td><td>2.00</td><td>2.00</td><td>2.00</td><td>2.00</td><td>2.00</td><td>2.00</td></tr><tr><td>2.25</td><td>2.00</td><td>2.00</td><td>2.00</td><td>2.00</td><td>1.60</td><td>1.77</td></tr><tr><td>2.50</td><td>2.00</td><td>2.00</td><td>2.00</td><td>2.00</td><td>1.33</td><td>1.60</td></tr><tr><td>3.00</td><td>2.00</td><td>2.00</td><td>2.00</td><td>2.00</td><td>1.00</td><td>1.33</td></tr></tbody></table>

polation errors dominate all the error quantities defined in (93). For $\delta = 0$ the results are less pronounced compared to the case $\delta = 1$. Let us analyze the velocity errors in the discrete norm. For $r = 2.25$ the averaged rates 1.80 and 1.87 are close to the rate $4/r$. For $r = 2.50$ velocity errors have rates 1.50 and 1.52, which fall between the rates $2/(r-1)$ and $4/r$. For $r = 3.00$ we observe rates 1.13 and 0.85, with the expected rate $2/(r-1)$ nearly attained. Similar rates are observed for the continuous norm. The pressure errors exhibit in general better rates, lying between $4/r$ and 2.
QUADRILATERAL MESHES

<figure><img src="image_4.png" /><figcaption>Figure 2: Test 1. Computed errors defined as in (93) as a function of the mesh size (loglog scale), for the mesh family QUADRILATERAL. Left panel: $\delta = 1$, right panel: $\delta = 0$.</figcaption></figure>
RANDOM MESHES

<figure><img src="image_5.png" /><figcaption>Figure 3: Test 1. Computed errors defined as in (93) as a function of the mesh size (log-log scale), for the mesh family RANDOM. Left panel: $\delta = 1$, right panel: $\delta = 0$.</figcaption></figure>
## 5.2 Test 2. Polynomial solution

To further investigate how the different sources of error combine in the error estimates, we consider Problem (1) on $\Omega = (0, 1)^2$ where the Dirichlet datum and the loading term are chosen in accordance with the exact solution

$$ u_{\text{ex}}(x_1, x_2) = \begin{bmatrix} x_1^2 + x_2^2 + 3x_1 + 5 \\ -2x_1x_2 - x_1^2 - 3x_2 + 7 \end{bmatrix}, \quad p_{\text{ex}}(x_1, x_2) = 0. $$

We notice that $u_{\text{ex}} \in [\mathbb{P}_2(\Omega)]^2 \subseteq U_h$, hence, by Theorem 19, we have $R_1 = 0$. As a consequence the asymptotically dominant contribution to the error arising from the approximation of $\sigma(\cdot, \epsilon(u_{\text{ex}}))$ should be better appreciated (with less or no influence by the other terms). In Fig. 4 we show the error quantities in (93) (for the pressures we plot the absolute errors) for the sequence of QUADRILATERAL meshes and parameters $r$ and $\delta$ in (94).

<figure><img src="image_6.png" /><figcaption>Figure 4: Test 2. Computed errors defined as in (93) as a function of the mesh size (loglog scale), for the mesh family QUADRILATERAL. Left panel: $\delta = 1$, right panel: $\delta = 0$.</figcaption></figure>
As expected, for $r = 2$, we recover the so-called “patch test”, i.e. the discrete solution and the exact solution coincide up to machine precision. For $r > 2$ the average experimental order of convergence AEOC of the $err(u_h, ||\cdot||_r)$ is in good agreement with rates predicted by Corollary 21 for $\delta = 1$ and Theorem 19 for $\delta = 0$, namely $4/r$ and $2/(r - 1)$ respectively (cf. Table 1). The pressure errors exhibit in both cases better rates.

## 5.3 Test 3. Singular solution

The purpose of this test is to assess the performance of the method in the presence of solutions with low Sobolev regularity. To this end, we examine the behavior of the proposed method for the benchmark test

<table><thead><tr><th colspan="5">err(u<sub>h</sub>, ||·||<sub>r</sub>)</th></tr><tr><th>1/h</th><th>2.00</th><th>2.25</th><th>r</th><th>2.50</th><th>3.00</th></tr></thead><tbody><tr><td>4</td><td>7.577256e-04</td><td>1.461758e-03</td><td></td><td>3.500512e-03</td><td>1.119863e-02</td></tr><tr><td>8</td><td>3.772397e-04</td><td>8.214921e-04</td><td></td><td>2.262381e-03</td><td>8.787189e-03</td></tr><tr><td>16</td><td>1.874426e-04</td><td>4.634299e-04</td><td></td><td>1.468054e-03</td><td>6.907361e-03</td></tr><tr><td>32</td><td>9.308978e-05</td><td>2.622244e-04</td><td></td><td>9.550462e-04</td><td>5.434720e-03</td></tr><tr><td>AEOC</td><td>1.0008326e+00</td><td>8.262771e-01</td><td></td><td>6.246411e-01</td><td>3.476817e-01</td></tr><tr><td><sup>4</sup>&frasl;<sub>r<sup>2</sup></sub></td><td>1.00</td><td>0.79</td><td></td><td>0.64</td><td>0.44</td></tr></tbody></table>

<table><thead><tr><th colspan="5">err(u<sub>h</sub>, W<sup>1,r</sup>)</th></tr><tr><th>1/h</th><th>2.00</th><th>2.25</th><th>r</th><th>2.50</th><th>3.00</th></tr></thead><tbody><tr><td>4</td><td>7.576044e-04</td><td>1.457011e-03</td><td></td><td>3.487738e-03</td><td>1.117410e-02</td></tr><tr><td>8</td><td>3.772245e-04</td><td>8.207718e-04</td><td></td><td>2.260225e-03</td><td>8.782263e-03</td></tr><tr><td>16</td><td>1.874407e-04</td><td>4.633219e-04</td><td></td><td>1.467695e-03</td><td>6.906383e-03</td></tr><tr><td>32</td><td>9.308954e-05</td><td>2.622083e-04</td><td></td><td>9.549867e-04</td><td>5.434526e-03</td></tr><tr><td>AEOC</td><td>1.0008251e+00</td><td>8.247422e-01</td><td></td><td>6.229130e-01</td><td>3.466444e-01</td></tr><tr><td><sup>4</sup>&frasl;<sub>r<sup>2</sup></sub></td><td>1.00</td><td>0.79</td><td></td><td>0.64</td><td>0.44</td></tr></tbody></table>

<table><caption>Table 2: Test 3. Computed errors $err(u_h, ||\cdot||_r)$ (top), $err(u_h, W^{1,r})$ (middle), and $err(p_h, L^{r'})$ (bottom) as in (93) for the mesh family CARTESIAN: $\delta = 1$.</caption><thead><tr><th colspan="5">err(p<sub>h</sub>, L<sup>r'</sup>)</th></tr><tr><th>1/h</th><th>2.00</th><th>2.25</th><th>r</th><th>2.50</th><th>3.00</th></tr></thead><tbody><tr><td>4</td><td>1.173547e-01</td><td>1.170544e-01</td><td></td><td>1.238788e-01</td><td>1.354350e-01</td></tr><tr><td>8</td><td>5.832495e-02</td><td>5.822885e-02</td><td></td><td>6.167576e-02</td><td>6.752532e-02</td></tr><tr><td>16</td><td>2.896772e-02</td><td>2.892814e-02</td><td></td><td>3.065185e-02</td><td>3.358190e-02</td></tr><tr><td>32</td><td>1.438465e-02</td><td>1.436589e-02</td><td></td><td>1.522427e-02</td><td>1.668436e-02</td></tr><tr><td>AEOC</td><td>1.0009424e+00</td><td>1.0008820e+00</td><td></td><td>1.0008161e+00</td><td>1.0007010e+00</td></tr><tr><td>||p<sub>ex</sub> - p<sub>I</sub>||<sub>L<sup>r'</sup></sub></td><td>1.00</td><td>1.00</td><td></td><td>1.00</td><td>1.00</td></tr></tbody></table>

introduced in [20, Section 7]. We consider Problem (1) on the square domain $\Omega = (-1, 1)^2$, where the forcing term $f$ (depending on $r$ and $\delta$ in (2)) and the Dirichlet boundary conditions prescribed on $\partial\Omega$ are chosen in accordance with the exact solution

$$u_{\text{ex}}(x_1, x_2) = |x|^{0.01} \left[ \begin{array}{l} x_2 \\ -x_1 \end{array} \right], \quad p_{\text{ex}}(x_1, x_2) = -|x|^\gamma + c_\gamma,$$

where $\gamma = \frac{2}{r} - 1 + 0.01$ and $c_\gamma$ is s.t. $p_{\text{ex}}$ is zero averaged. We note that for all $r \in [2, \infty)$

$$u_{\text{ex}} \in W^{2/r+1,r}(\Omega), \quad \sigma(\cdot, \epsilon(u_{\text{ex}})) \in W^{2/r',r'}(\Omega), \quad f \in W^{2/r'-1,r'}(\Omega), \quad p_{\text{ex}} \in W^{1,r'}(\Omega),$$
therefore, with the notation of Theorem 19 and Theorem 22:

$$k_1 = \frac{2}{r}, \quad k_2 = \frac{2}{r'}, \quad k_3 = \frac{2}{r'} - 2, \quad k_4 = 1.$$

The domain $\Omega$ is partitioned with a sequence of CARTESIAN meshes with diameter $h=1/4, 1/8, 1/16, 1/32$ (see Fig. 1). Table 2 presents the computed errors $err(u_h, \|\cdot\|_r)$ (top panel), $err(u_h, W^{1,r})$ (middle panel), and $err(p_h, L^{r'})$ (bottom panel), cf. (93), and the associated average experimental orders of convergence AEOC, for the values of $r = 2.00, 2.25, 2.50, 3.00$ and $\delta = 1$. The corresponding results for $\delta = 0$ are shown in Table 3. Notice that, according to Theorem 19 and Corollary 20, the expected rate of convergence for the velocity in both the discrete and continuous norm is $2k/r^2$. Finally, for the error $err(\sigma, L^{r'})$, linear convergence is observed.

<table><thead><tr><th colspan="5">err(u<sub>h</sub>, ||·||<sub>r</sub>)</th></tr><tr><th>1/h</th><th colspan="4">r</th></tr><tr><th></th><th>2.00</th><th>2.25</th><th>2.50</th><th>3.00</th></tr></thead><tbody><tr><th>4</th><td>7.577256e-04</td><td>3.757320e-03</td><td>1.784934e-02</td><td>1.572178e-01</td></tr><tr><th>8</th><td>3.772397e-04</td><td>2.139041e-03</td><td>1.114599e-02</td><td>1.064344e-01</td></tr><tr><th>16</th><td>1.874426e-04</td><td>1.222004e-03</td><td>6.981838e-03</td><td>7.427671e-02</td></tr><tr><th>32</th><td>9.308978e-05</td><td>6.988137e-04</td><td>4.376918e-03</td><td>5.255378e-02</td></tr><tr><th>AEOC</th><td>1.008326e+00</td><td>8.0890813-01</td><td>6.759612e-01</td><td>5.269660e-01</td></tr><tr><th><sup>4</sup>/<sub>r<sup>2</sup></sub></th><td>1.00</td><td>0.79</td><td>0.64</td><td>0.44</td></tr></tbody></table>

<table><thead><tr><th colspan="5">err(u<sub>h</sub>, W<sup>1,r</sup>)</th></tr><tr><th>1/h</th><th colspan="4">r</th></tr><tr><th></th><th>2.00</th><th>2.25</th><th>2.50</th><th>3.00</th></tr></thead><tbody><tr><th>4</th><td>7.576044e-04</td><td>3.732379e-03</td><td>1.774680e-02</td><td>1.564706e-01</td></tr><tr><th>8</th><td>3.772245e-04</td><td>2.135423e-03</td><td>1.112967e-02</td><td>1.063065e-01</td></tr><tr><th>16</th><td>1.874407e-04</td><td>1.221485e-03</td><td>6.979270e-03</td><td>7.425567e-02</td></tr><tr><th>32</th><td>9.308954e-05</td><td>6.987394e-04</td><td>4.376515e-03</td><td>5.255063e-02</td></tr><tr><th>AEOC</th><td>1.008251e+00</td><td>8.057564e-01</td><td>6.732348e-01</td><td>5.247039e-01</td></tr><tr><th><sup>4</sup>/<sub>r<sup>2</sup></sub></th><td>1.00</td><td>0.79</td><td>0.64</td><td>0.44</td></tr></tbody></table>

<table><caption>Table 3: Test 3. Computed errors $err(u_h, ||\cdot||_r)$ (top), $err(u_h, W^{1,r})$ (middle), and $err(p_h, L^{r'})$ (bottom) as in (93) for the mesh family CARTESIAN: $\delta = 0$.</caption><thead><tr><th colspan="5">err(p<sub>h</sub>, L<sup>r'</sup>)</th></tr><tr><th>1/h</th><th colspan="4">r</th></tr><tr><th></th><th>2.00</th><th>2.25</th><th>2.50</th><th>3.00</th></tr></thead><tbody><tr><th>4</th><td>1.173547e-01</td><td>1.170279e-01</td><td>1.237627e-01</td><td>1.349763e-01</td></tr><tr><th>8</th><td>5.832495e-02</td><td>5.821788e-02</td><td>6.162369e-02</td><td>6.731919e-02</td></tr><tr><th>16</th><td>2.896772e-02</td><td>2.892321e-02</td><td>3.062604e-02</td><td>3.348035e-02</td></tr><tr><th>32</th><td>1.438465e-02</td><td>1.436359e-02</td><td>1.521138e-02</td><td>1.663401e-02</td></tr><tr><th>AEOC</th><td>1.009424e+00</td><td>1.008788e+00</td><td>1.008117e+00</td><td>1.006832e+00</td></tr><tr><th>||p<sub>ex</sub> - p<sub>I</sub>||<sub>L<sup>r'</sup></sub></th><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td></tr></tbody></table>

# 6 Conclusions

We presented a theoretical analysis of Virtual Element discretizations of incompressible non-Newtonian flows governed by the Carreau–Yasuda constitutive law in the shear-thickening regime ($r > 2$). Our analysis also covers the degenerate limit ($\delta = 0$), which corresponds to the power-law model. The proposed Virtual
Element method is fully compatible with general polygonal meshes and yields an exactly divergence-free discrete velocity field. To carry out the analysis, we introduced novel theoretical tools, including an inf-sup stability bound in non-Hilbertian norms, a suitably tuned stabilization for the case $r > 2$, and discrete norms consistent with the constitutive law. We presented numerical results to demonstrate the theoretical findings and assess the practical performance of the proposed method. The present results extend and complete those in [5], which covered the case $1 < r < 2$ (shear-thinning regime), demonstrating that the VEM provides a robust discretization framework for Carreau-Yasuda non-Newtonian flows in both shear-thickening and shear-thinning regimes.

## Acknowledgments

This research has been partially funded by the European Union (ERC, NEMESIS, project number 101115663). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. The present research is part of the activities of “Dipartimento di Eccellenza 2023-2027”. PA and MV also acknowledge MUR-PRIN/PNRR 2022 grant n. P2022BH5CB, funded by MUR. GV has been partially funded by PRIN2022 n. 2022MBY5JM “FREYA - Fault REactivation: a hYbrid numerical Approach” research grant, and PRIN2022PNRR n. P2022M7JZW “SAFER MESH - Sustainable mAnagement oF watEr Resources: ModEls and numerical MetHods” research grant, funded by the Italian Ministry of Universities and Research (MUR) and by the European Union through Next Generation EU, M4C2. GV acknowledges financial support of INdAM-GNCS through project CUPE53C24001950001 “VEM per la trattazione di problemi definiti su domini parametrici o randomici”. The authors are members of INdAM-GNCS.

## References

[1] D. Adak, D. Mora, S. Natarajan, and A. Silgado. A virtual element discretization for the time dependent Navier-Stokes equations in stream-function formulation. *ESAIM Math. Model. Numer. Anal.*, 55(5):2535-2566, 2021.

[2] D. Adak, D. Mora, and A. Silgado. The Morley-type virtual element method for the Navier-Stokes equations in stream-function form. *Comput. Methods Appl. Mech. Engrg.*, 419:Paper No. 116573, 28, 2024.

[3] P. F. Antonietti, L. Beirão da Veiga, D. Mora, and M. Verani. A stream virtual element formulation of the Stokes problem on polygonal meshes. *SIAM J. Numer. Anal.*, 52(1):386-404, 2014.

[4] P. F. Antonietti, G. Vacca, and M. Verani. Virtual element method for the Navier-Stokes equation coupled with the heat equation. *IMA J. Numer. Anal.*, 43(6):3396-3429, 2023.

[5] P.F. Antonietti, L. Beirão da Veiga, M. Botti, G. Vacca, and M. Verani. A virtual element method for non-newtonian pseudoplastic stokes flows. *Computer Methods in Applied Mechanics and Engineering*, 428:117079, 2024.

[6] P.F. Antonietti, L. Beirão da Veiga, and G. Manzini. *The Virtual Element Method and its Applications*. SEMA SIMAI Springer series 31, Springer International Publishing, 2022.

[7] P.F. Antonietti, L. Beirão da Veiga, D. Mora, and M. Verani. A stream virtual element formulation of the Stokes problem on polygonal meshes. *SIAM J. Numer. Anal.*, 52(1):386-404, 2014.

[8] J. Baranger and K. Najib. Analyse numérique des écoulements quasi-newtoniens dont la viscosité obéit à la loi puissance ou la loi de carreau. *Numer. Math.*, 58(1):35-49, 1990.

[9] J. W. Barrett and W. B. Liu. Quasi-norm error bounds for the finite element approximation of a non-Newtonian flow. *Numer. Math.*, 68(4):437-456, 1994.

[10] J.W. Barrett and W.B. Liu. Finite element error analysis of a quasi-Newtonian flow obeying the Carreau or power law. *Numer. Math.*, 64(4):433-453, 1993.

[11] H. Beirão da Veiga. On the global regularity of shear thinning flows in smooth domains. *J. Math. Anal. Appl.*, 349(2):335-360, 2009.

[12] L. Beirão da Veiga, F. Dassi, D. A. Di Pietro, and J. Droniou. Arbitrary-order pressure-robust DDR and VEM methods for the Stokes problem on polyhedral meshes. *Comput. Methods Appl. Mech. Engrg.*, 397:Paper No. 115061, 31, 2022.
[13] L. Beirão da Veiga, F. Dassi, and G. Vacca. The Stokes complex for virtual elements in three dimensions. *Math. Models Methods Appl. Sci.*, 30(3):477–512, 2020.

[14] L. Beirão da Veiga, F. Dassi, and G. Vacca. Vorticity-stabilized virtual elements for the Oseen equation. *Math. Models Methods Appl. Sci.*, 31(14):3009–3052, 2021.

[15] L. Beirão da Veiga, C. Lovadina, and G. Vacca. Divergence free virtual elements for the Stokes problem on polygonal meshes. *ESAIM Math. Model. Numer. Anal.*, 51(2):509–535, 2017.

[16] L. Beirão da Veiga, C. Lovadina, and G. Vacca. Virtual elements for the Navier-Stokes problem on polygonal meshes. *SIAM J. Numer. Anal.*, 56(3):1210–1242, 2018.

[17] L. Beirão da Veiga, D. Mora, and A. Silgado. A fully-discrete virtual element method for the nonstationary Boussinesq equations in stream-function form. *Comput. Methods Appl. Mech. Engrg.*, 408:115947, 2023.

[18] L. Beirão da Veiga, D. Mora, and G. Vacca. The Stokes complex for virtual elements with application to Navier-Stokes flows. *J. Sci. Comput.*, 81(2):990–1018, 2019.

[19] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo. Basic principles of virtual element methods. *Math. Models Methods Appl. Sci.*, 23(1):199–214, 2013.

[20] L. Belenki, L. C. Berselli, L. Diening, and M. Růžička. On the finite element approximation of $p$-Stokes systems. *SIAM J. Numer. Anal.*, 50(2):373–397, 2012.

[21] L.C. Berselli, L. Diening, and M. Růžička. Existence of strong solutions for incompressible fluids with shear dependent viscosities. *J. Math. Fluid Mech.*, 12(1):101–132, 2010.

[22] T. Bevilacqua, F. Dassi, S. Zampini, and S. Scacchi. BDDC Preconditioners for Virtual Element Approximations of the Three-Dimensional Stokes Equations. *SIAM J. Sci. Comput.*, 46(1):A156–A178, 2024.

[23] T. Bevilacqua and S. Scacchi. BDDC preconditioners for divergence free virtual element discretizations of the Stokes equations. *J. Sci. Comput.*, 92(2):Paper No. 63, 27, 2022.

[24] M. Bogovskiï. Solution of the first boundary value problem for an equation of continuity of an incompressible medium. *Dokl. Akad. Nauk SSSR*, 248(5):1037–1040, 1979.

[25] M. Botti, D. Castanon Quiroz, D. A. Di Pietro, and A. Harnist. A Hybrid High-Order method for creeping flows of non-Newtonian fluids. *ESAIM: Math. Model Numer. Anal.*, 55(5):2045–2073, 2021.

[26] S. C. Brenner and L. R. Scott. *The Mathematical Theory of Finite Element Methods*, volume 15 of *Texts in Applied Mathematics*. Springer, New York, third edition, 2008.

[27] E. Cáceres and G. N. Gatica. A mixed virtual element method for the pseudostress-velocity formulation of the Stokes problem. *IMA J. Numer. Anal.*, 37(1):296–331, 2017.

[28] E. Cáceres, G. N. Gatica, and F. A. Sequeira. A mixed virtual element method for quasi-Newtonian Stokes flows. *SIAM J. Numer. Anal.*, 56(1):317–343, 2018.

[29] A. Cangiani, V. Gyrya, and G. Manzini. The nonconforming virtual element method for the Stokes equations. *SIAM J. Numer. Anal.*, 54(6):3411–3435, 2016.

[30] D. Castanon Quiroz, D.A. Di Pietro, and A. Harnist. A Hybrid High-Order method for incompressible flows of non-newtonian fluids with power-like convective behaviour. *IMA J. Numer. Anal.*, 43(1):144–186, 12 2021.

[31] A. Chernov, C. Marcati, and L. Mascotto. $p$- and $hp$-virtual elements for the Stokes problem. *Adv. Comput. Math.*, 47(2):Paper No. 24, 31, 2021.

[32] P. G. Ciarlet and P. Ciarlet. Another approach to linearized elasticity and a new proof of Korn's inequality. *Math. Models Methods Appl. Sci.*, 15(02):259–271, 2005.

[33] M. Dauge. Stationary Stokes and Navier-Stokes systems on two- or three-dimensional domains with corners. I. Linearized equations. *SIAM J. Math. Anal.*, 20(1):74–97, 1989.

[34] K. Deimling. *Nonlinear functional analysis*. Springer-Verlag, Berlin, 1985.

[35] L. Diening and F. Ettwein. Fractional estimates for non-differentiable elliptic systems with general growth. *Forum Math.*, 20(3):523–556, 2008.

[36] L. Diening and C. Kreuzer. Linear convergence of an adaptive finite element method for the $p$-laplacian equation. *SIAM Journal on Numerical Analysis*, 46(2):614–638, 2008.
[37] D. Frerichs and C. Merdon. Divergence-preserving reconstructions on polygons and a really pressure-robust virtual element method for the Stokes problem. *IMA J. Numer. Anal.*, 42(1):597–619, 2022.

[38] G. P. Galdi, C. G. Simader, and H. Sohr. On the Stokes problem in Lipschitz domains. *Annali di Matematica Pura ed Applicata*, 167(1):147–163, 1994.

[39] G. N. Gatica, M. Munar, and F. A. Sequeira. A mixed virtual element method for the Navier-Stokes equations. *Math. Models Methods Appl. Sci.*, 28(14):2719–2762, 2018.

[40] G. N. Gatica and F. A. Sequeira. An $L^p$ spaces-based mixed virtual element method for the two-dimensional Navier-Stokes equations. *Math. Models Methods Appl. Sci.*, 31(14):2937–2977, 2021.

[41] G.N. Gatica and F. A. Sequeira. Analysis of an augmented HDG method for a class of quasi-newtonian Stokes flows. *J. Sci. Comput.*, 65(3):1270–1308, 2015.

[42] G. Geymonat and P. M. Suquet. Functional spaces for Norton-Hoff materials. *Math. Methods Appl. Sci.*, 8(2):206–222, 1986.

[43] A. Hirn. Approximation of the $p$-Stokes equations with equal-order finite elements. *J. Math. Fluid Mech.*, 15(1):65–88, 2013.

[44] Geng J. and Kilty J. The $L^p$ regularity problem for the Stokes system on Lipschitz domains. *Journal of Differential Equations*, 259(4):1275–1296, 2015.

[45] A. Kaltenbach and M. Růžička. A local discontinuous Galerkin approximation for the $p$-Navier-Stokes system, part i: Convergence analysis. *SIAM Journal on Numerical Analysis*, 61(4):1613–1640, 2023.

[46] R. B. Kellogg and J. E. Osborn. A regularity result for the Stokes problem in a convex polygon. *J. Functional Analysis*, 21(4):397–431, 1976.

[47] C. Kreuzer and E. Süli. Adaptive finite element approximation of steady flows of incompressible fluids with implicit power-law-like rheology. *ESAIM Math. Model. Numer. Anal.*, 50(5):1333–1369, 2016.

[48] F. Lepe and G. Rivera. A virtual element approximation for the pseudostress formulation of the Stokes eigenvalue problem. *Comput. Methods Appl. Mech. Engrg.*, 379:Paper No. 113753, 21, 2021.

[49] Y. Li, C. Hu, and M. Feng. A Least-Squares stabilization virtual element method for the Stokes problem on polygonal meshes. *Int. J. Numer. Anal. Model.*, 19(5):685–708, 2022.

[50] Y. Li, C. Hu, and M. Feng. On stabilized equal-order virtual element methods for the Navier-Stokes equations on polygonal meshes. *Comput. Math. Appl.*, 154:267–286, 2024.

[51] T. Malkmus, M. Růžička, S. Eckstein, and I. Toulopoulos. Generalizations of SIP methods to systems with $p$-structure. *IMA J. Numer. Anal.*, 38(3):1420–1451, 09 2017.

[52] J. Meng, L. Beirão da Veiga, and L. Mascotto. Stability and interpolation properties for Stokes-like virtual element spaces. *J. Sci. Comput.*, 94(56), 2023.

[53] D. Sandri. Sur l'approximation numérique des écoulements quasi-newtoniens dont la viscosité suit la loi puissance ou la loi de carreau. *M2AN Math. Model. Numer. Anal.*, 27(2):131–155, 1993.

[54] G. Savaré. Regularity results for elliptic equations in Lipschitz domains. *Journal of Functional Analysis*, 152(1):176–201, 1998.

[55] C. Talischi, G. H. Paulino, A. Pereira, and I-F. M. Menezes. *PolyMesher*: a general-purpose mesh generator for polygonal elements written in Matlab. *Struct. Multidiscip. Optim.*, 45(3):309–328, 2012.

[56] G. Vacca. An $H^1$-conforming virtual element for Darcy and Brinkman equations. *Math. Models Methods Appl. Sci.*, 28(1):159–194, 2018.

[57] G. Wang and Y. Wang. Least-Squares Virtual Element Method for Stokes Problems on Polygonal Meshes. *J. Sci. Comput.*, 98(2):Paper no. 46, 2024.

[58] K. Yasuda, R.C. Armstrong, and R.E. Cohen. Shear flow properties of concentrated solutions of linear and star branched polystyrenes. *Rheologica Acta*, 20(2):163–178, 3 1981.