Convergence analysis of an augmented dual-mixed fnite element method with skew-symmetric Nitsche stabilization for the generalized Oseen equations with mixed boundary conditions

Introduction

The generalized Oseen equations arise from the linearization of both the steady and unsteady Navier–Stokes equations, which describe the motion of viscous incompressible fluids in a wide range of physical settings. The Navier–Stokes equations play a fundamental role in numerical analysis and scientific computing, with broad applications in engineering, geophysics, and biomedical sciences, among other fields (see, e.g., [30, 35, 46]). The numerical solution of these equations remains challenging due to the nonlinear convective term, particularly at high Reynolds numbers, where convection may dominate diffusion. In this context, generalized Oseen-type problems naturally arise at each step of fixed-point or Newton iterations applied to the Navier–Stokes equations (see, e.g., [45, 52]).

Classical mixed finite element methods have become widely used tools for the numerical approximation of fluid flows, as they allow for the simultaneous approximation of primal variables, such as velocity and pressure, together with additional variables of physical interest, often referred to as dual variables, such as the pseudostress tensor or the velocity gradient, depending on the formulation. Furthermore, other relevant quantities, such as vorticity, can be accurately recovered through suitable postprocessing techniques [20, 25, 38]. The well-posedness of these formulations typically relies on an inf-sup condition, also known as the Ladyzhenskaya–Babuška–Brezzi (LBB) condition, which must be satisfied by the discrete approximation spaces and may impose compatibility restrictions between the spaces associated with the different unknowns. An alternative approach consists of incorporating residual-based stabilization terms derived from the governing equations of the model problem, such as the conservation of momentum and mass equations, with the aim of allowing more general choices of finite element approximation spaces. This type the formulation is knows in the literature as augmented formulations (for an incomplete list, see for example [7, 26, 34, 39]).

The imposition of boundary conditions in finite element formulations can be handled in different ways. For instance, a natural boundary condition in a primal formulation becomes an essential condition in the corresponding dual-mixed formulation [20, 38]. Consequently, such a condition, as in the case of Neumann boundary conditions, must be incorporated into the corresponding continuous function space, which may complicate the construction and practical implementation of suitable finite element subspaces satisfying the prescribed boundary conditions. One way to avoid imposing this constraint directly on the approximation space is to enforce the essential boundary condition weakly by means of a Lagrange multiplier [10, 36, 37, 40]. However, this approach introduces an additional unknown and requires an additional finite element approximation space, together with a discrete inf-sup condition that must be satisfied. Another approach consists of homogenizing the essential boundary conditions by means of a suitable lifting. Although this technique allows the problem to be reformulated in terms of homogeneous boundary conditions, its discrete implementation may require constructing a lifting compatible with the finite element approximation spaces [14, 15, 16]. For problems involving interface conditions between different subdomains, mortar finite element methods [18] provide a suitable alternative, particularly when the meshes do not match across the interface. These methods are also of interest in elasticity, especially for contact problems involving two deformable bodies [42].

Another approach for the weak imposition of essential boundary conditions is the classical penalty method [9, 11], which incorporates a mesh-dependent penalty term into the variational formulation without introducing additional unknowns. However, this method is not consistent, and its convergence properties depend on the choice of the penalty parameter (see, e.g., [29, 43]). An alternative is Nitsche’s method, originally introduced in [49], which weakly enforces essential boundary conditions while preserving consistency (see, e.g., [29, 51]). Under suitable assumptions, Nitsche’s method allows optimal convergence rates to be achieved without increasing the number of discrete unknowns [17, 49]. Both symmetric and nonsymmetric formulations are available, the latter including skew-symmetric variants. While the classical symmetric formulation requires a sufficiently large penalty parameter to ensure coercivity, the stability of a penalty-free nonsymmetric Nitsche formulation has also been established for Poisson and convection–diffusion problems [23]. Furthermore, like the classical penalty method, Nitsche’s method avoids introducing additional boundary unknowns and the associated saddle-point structure arising in Lagrange multiplier formulations.

In recent years, there has been growing interest in the use of Nitsche’s method for the weak imposition of nonstandard boundary conditions in different variational formulations. Among the numerous contributions involving primal formulations, we first mention [28], where Nitsche’s method is employed to impose nonhomogeneous Robin boundary conditions for the Poisson problem. Also for the Poisson problem, but considering Dirichlet boundary conditions, a symmetric Nitsche formulation combined with least-squares and finite cell stabilization is proposed in [47], within the framework of unfitted finite element methods. In the context of linear elasticity, a Nitsche-based finite element approximation of frictionless unilateral contact conditions of Signorini type is studied in [22], with particular emphasis on the convergence analysis arising in shape optimization.

On the other hand, concerning mixed formulations, a Nitsche-based mixed finite element method is proposed in [5] for the Poisson problem with mixed boundary conditions, in which the Neumann boundary condition is weakly enforced through a symmetric Nitsche formulation. The discrete scheme employs Brezzi–Douglas–Marini (BDM) elements, and the approach is further extended to the Poisson problem with Signorini boundary conditions. In a different context, two strategies inspired by contact mechanics are proposed in [41] for the weak imposition of nonlinear seepage boundary conditions in Richards’ equation. The first approach relies on a strongly consistent penalty term inspired by Nitsche’s method, whereas the second employs a hybridized formulation in which the boundary pressure is introduced as an additional unknown.

For fluid flow problems, three variants of Nitsche’s method are studied in [8] for the Stokes equations with slip boundary conditions. These schemes combine the weak imposition of essential boundary conditions through Nitsche terms with a stabilized velocity-pressure formulation based on equal-order $P_k^n \times P_k^n$ Lagrange finite elements. For the Brinkman problem, a symmetric Nitsche scheme based on a velocity-pressure formulation with grad-div stabilization is proposed in [6] to weakly impose Dirichlet and tangential boundary conditions. In the context of porous media flow, two mixed finite element formulations for the Darcy problem with Dirichlet and Neumann boundary conditions are analyzed in [24], using Raviart–Thomas elements conforming in $H(\operatorname{div}; \Omega)$. The first formulation is consistent and belongs to the class of Nitsche-type methods, whereas the second follows a penalty approach. Further applications of Nitsche’s method to fluid flow problems can be found, among others, in [12, 13, 19].

More recently, a skew-symmetric Nitsche scheme based on an augmented dual-mixed finite element formulation has been proposed in [27] for the Brinkman problem with mixed Dirichlet–Neumann boundary conditions. The augmented formulation incorporates additional terms derived from the governing equations, allowing greater flexibility in the choice of finite element approximation spaces without requiring the discrete inf-sup compatibility conditions of the underlying mixed formulation.

Motivated by the preceding developments, in this work we build upon the frameworks introduced in [16, 15] to propose and analyze a new skew-symmetric Nitsche scheme for the generalized Oseen equations with non-solenoidal convective fields, mixed Dirichlet–Neumann boundary conditions, and non-null source terms. The method is based on an augmented dual-mixed pseudostress–velocity formulation, whose well-posedness at the continuous level is established by means of a suitable lifting of the Neumann datum. At the discrete level, the essential boundary condition is weakly enforced through consistent skew-symmetric Nitsche terms, using Raviart–Thomas elements for the pseudostress tensor and Lagrange elements for the velocity. We establish the existence and uniqueness of the discrete solution and prove the stability of the proposed scheme for homogeneous Dirichlet boundary conditions and general Neumann data. Furthermore, in contrast to [27], we derive optimal convergence rates in an energy norm under weaker regularity assumptions on the exact solution. The proposed scheme does not introduce additional unknowns into the augmented formulation and can be readily implemented using the finite element library FEniCS [2]. Finally, numerical experiments in two and three dimensions confirm the theoretical convergence rates and illustrate the performance of the method for non-solenoidal convective fields and convection-dominated flow regimes.

The remainder of this paper is organized as follows. In § 2, we introduce the model problem and the corresponding functional setting. We then recall the pseudostress–velocity formulation and establish the well-posedness of the resulting augmented variational formulation under suitable assumptions on the stabilization parameters. In § 3, we introduce the finite element approximation spaces and formulate the discrete scheme, in which the essential boundary condition is weakly enforced through suitable Nitsche terms. We also establish the existence and uniqueness of the discrete solution. § 4 is devoted to the error analysis. We first establish the error orthogonality property and derive suitable approximation properties of the interpolation operators, which are subsequently used to derive an a priori error estimate in the energy norm. Finally, in § 5, we present numerical experiments that confirm the theoretical convergence rates and illustrate the performance of the proposed method.

Model problem and preliminary results

Function spaces

Throughout this paper, let $\Omega\subset\mathbb{R}^{n}$, $n \in\{2,3\}$, be a bounded Lipschitz domain with boundary $\Gamma:= \Gamma_D \cup\Gamma_N$, where $\Gamma_D$ and $\Gamma_N$ are nonempty and disjoint. We denote by $\boldsymbol{n}$ the outward unit normal vector to $\Gamma$. We adopt standard notation for Sobolev spaces, inner products, and norms (see, e.g., [33]). In particular, for any open set $\mathcal{O} \subset\mathbb{R}^{n}$, we denote by $L^2(\mathcal{O})$ the space of square-integrable functions on $\mathcal{O}$. Moreover, $\lVert\cdot\rVert_{m,\mathcal{O}}$ and $|\cdot|_{m,\mathcal{O}}$ denote the standard norm and seminorm, respectively, in $H^m(\mathcal{O})$, $H^m(\mathcal{O})^n$, and $H^m(\mathcal{O})^{n \times n}$, for $m \ge0$, with the convention $H^0(\mathcal{O})=L^2(\mathcal{O})$. On the other hand, given a scalar field $q$ and a vector field $\mathbf{v}=(v_i)_{i=1,\ldots,n}$, we set the laplacian, gradient and divergence operators, respectively by

Furthermore, the Laplacian of a vector field is defined componentwise, whereas the divergence of a tensor field is taken row-wise. We also introduce the Hilbert spaces $H(\operatorname{div};\Omega):=\{\mathbf{v}\in L^2(\Omega)^n:\nabla\cdot\mathbf{v}\in L^2(\Omega)\}$ and $H:=\{\boldsymbol{\tau}\in L^2(\Omega)^{n\times n}:\nabla\cdot\boldsymbol{\tau}\in L^2(\Omega)^n\}$, as well as the closed subspace $V:=\{\boldsymbol{\tau}\in H:\boldsymbol{\tau}\mathbf{n}=0\text{ on }\Gamma_N\}$, where the boundary condition is understood in the sense of normal traces. In all these cases, we consider the standard $H(\operatorname{div})$-norm $\|\cdot\|_{\operatorname{div};\Omega}$. For vector fields, we introduce the following Hilbert spaces:

We denote by $H^{-1/2}_{00}(E):=(H^{1/2}_{00}(E))'$ with the topological dual of the Lions–Magenes space $H^{1/2}_{00}(E)$, endowed with the duality pairing $\langle\cdot,\cdot\rangle_E$, where $E\in\{\Gamma_D,\Gamma_N\}$ (for details, see [50]).

Problem model

In this work, we consider the generalized Oseen problem: Find a velocity vector field $\mathbf{u}$ and a pressure scalar field $p$ such as

where $\nu>0$ is the kinematic viscosity, $\theta>0$ is a reaction coefficient, $\boldsymbol{\alpha}\in W^{1,\infty}(\Omega)^n$ denotes a convective velocity field, not necessarily solenoidal, $I$ is the identity matrix of $\mathbb{R}^{n\times n}$, $\mathbf{g}_D\in H^{1/2}_{00}(\Gamma_D)$, $\mathbf{g}_N\in H^{-1/2}_{00}(\Gamma_N)$, and $\mathbf{f}\in L^2(\Omega)^n$ is a body force vector. Moreover, we assume the following conditions on the convective term $\boldsymbol{\alpha}$:

(H1) $\|\nabla\cdot\boldsymbol{\alpha}\|_{\infty,\Omega}<\varepsilon$ in $\Omega$, where $\varepsilon>0$ is small;

(H2) $\boldsymbol{\alpha}\cdot\mathbf{n}\ge0$ on $\Gamma_N$.

Remark 1. Typically, the reaction parameter $\theta$ arises from a time discretization and is given by $\frac{1}{\Delta t}$, where $\Delta t$ denotes the time-step size. In turn, $\boldsymbol{\alpha}$ is commonly obtained from a linearization of the Navier–Stokes equations and is therefore regarded as a known velocity field. Since the underlying flow is incompressible, it is natural to assume that $\boldsymbol{\alpha}$ is nearly divergence-free, which motivates assumption (H1). Furthermore, assumption (H2) is naturally satisfied by the boundary conditions arising in many applications.

Next, at the continuous level, we follow the approach in [15, 16] to derive a mixed formulation of problem $(P)$. To this end, we introduce the pseudostress tensor $\boldsymbol{\sigma}:=\nu\nabla\mathbf{u}-pI$ as an auxiliary unknown. In addition, given a tensor $\boldsymbol{\tau}$, we denote its deviator by $\boldsymbol{\tau}^{d}:=\boldsymbol{\tau}-\frac{\operatorname{tr}(\boldsymbol{\tau})}{n}I$. With these definitions, problem $(P)$ can be equivalently reformulated as follows: Find a tensor field $\boldsymbol{\sigma}$ and a vector field $\mathbf{u}$ such that

The following lemma, whose proof is straightforward, establishes the equivalence between formulations $(P)$ and $(2.2)$.

Remark 2. The previous result shows that, from the solution of the strong problem $(2.2)$ formulated in terms of the pseudostress and velocity, the pressure $p$ of the original problem $(P)$ can be recovered. Moreover, by means of a straightforward postprocessing procedure, other quantities of interest can also be obtained. In particular, since $\boldsymbol{\sigma}^{d}=\nu\nabla\boldsymbol{u}$, the velocity gradient and the vorticity tensor are given, respectively, by $\nabla\boldsymbol{u}=\nu^{-1}\boldsymbol{\sigma}^{d}$ and $\boldsymbol{\omega}=\frac{1}{2\nu}\left(\boldsymbol{\sigma}^{d}-(\boldsymbol{\sigma}^{d})^{t}\right)$. For further details, we refer to [25].

The following lemma will be instrumental in the analysis of this work.

By testing with suitable functions, integrating by parts, and adding appropriate residual terms from the momentum equations in (2.2), we obtain the following augmented mixed formulation, which is derived using arguments similar to those in [15, 16]: Find $(\boldsymbol{\sigma},\boldsymbol{u})\in\mathbf{H}\times\mathbf{Q}$, with $\boldsymbol{\sigma}\boldsymbol{n}=-\boldsymbol{g}_{N}$ on $\Gamma_{N}$, such that

for all $(\boldsymbol{\tau},\boldsymbol{v})\in\mathbf{V}\times\mathbf{Q}$, where $A:(\mathbf{H}\times\mathbf{Q})\times(\mathbf{H}\times\mathbf{Q})\longrightarrow\mathbb{R}$ is the bilinear form defined by

for all $(\boldsymbol{\sigma},\boldsymbol{u}),(\boldsymbol{\tau},\boldsymbol{v})\in\mathbf{H}\times\mathbf{Q}$, and $F:\mathbf{H}\times\mathbf{Q}\longrightarrow\mathbb{R}$ is the linear form defined by

for all $(\boldsymbol{\tau},\boldsymbol{v})\in\mathbf{H}\times\mathbf{Q}$, with $\kappa_{1},\kappa_{2}$ and $\kappa_{3}$, positive parameters which will be specified later.

The following technical results will be used to establish the coercivity of the bilinear form $A$ on the subspace $\mathbf{V}\times\mathbf{Q}$ of $\mathbf{H}\times\mathbf{Q}$, as well as the existence and uniqueness of the solution to the proposed discrete scheme.

Then, for all $(\boldsymbol{\tau},v) \in\mathbf{H} \times Q$ there is a positive constant $C$, such that

Moreover, there exists another positive constant $C$, depending on the physical parameters and the domain $\Omega$, such that

for all $(\boldsymbol{\tau},v) \in\mathbf{V} \times Q$.

Proof. Let $(\boldsymbol{\tau},v) \in\mathbf{H} \times Q$. Then, using the definition of the bilinear form $A(\cdot,\cdot)$ (cf. (6)), Hölder and Young’s inequalities, (H1)–(H2) and Lemma 2, we obtain

and the estimate (9) follows from the conditions on the parameters given in (8). Moreover, the second lower bound (10) follows by restricting the previous estimate to $\mathbf{V} \times Q$ and applying [37 Lemma 4]. □

The following result guarantees the well-posedness of the variational formulation (5) by means of standard lifting and homogenization arguments. For completeness, we provide the details.

A skew-symmetric Nitsche scheme

Let $\{\mathcal{T}_h\}_{h>0}$ be a shape-regular family of partitions of $\overline{\Omega}$ into triangles for $n=2$ and tetrahedra for $n=3$. Let $h_T$ denote the diameter of an element $T\in\mathcal{T}_h$, and define $h:=\max_{T\in\mathcal{T}_h}h_T$. Moreover, let $E$ be an edge for $n=2$ or a face for $n=3$, and $\mathcal{E}_N$ denote the set of all edges/faces of $\mathcal{T}_h$ lying on $\Gamma_N$. For $T\in\mathcal{T}_h$, we define the local tensor-valued Raviart–Thomas space of order $k\geq0$ by

where $\tau_i$ denotes the $i$-th row of the tensor $\tau$, and $\operatorname{RT}_k(T):=\mathbb{P}_k(T)^n\oplus\mathbf{x}\mathbb{P}_k(T)$ is the classical Raviart–Thomas space of order $k$ on $T$ (see, e.g., [31 Chapter 3] and [38 Section 3.3]). Here, $\mathbb{P}_k(\mathcal{O})$ denotes the space of polynomials of total degree at most $k$ on $\mathcal{O}\subset\mathbb{R}^n$.

The finite element subspaces of $\mathbf{H}$ and $\mathbf{Q}$ employed in this work are defined, respectively, by:

with $k\geq0$, and we equip the product space $\mathbf{H}_{h,k}\times\mathbf{Q}_{h,k}$ with the norm:

To formulate the discrete scheme, we assume, for simplicity, that the Neumann datum satisfies $\mathbf{g}_N\in L^2(\Gamma_N)^n$. Under this assumption, the duality pairing involving $\mathbf{g}_N\in\mathbf{H}_{00}^{-1/2}(\Gamma_N)$ and functions in $\mathbf{H}_{00}^{1/2}(\Gamma_N)$ coincides with the $L^2(\Gamma_N)^n$ inner product. We now introduce the following discrete scheme, incorporating skew-symmetric Nitsche’s terms to impose the essential boundary conditions of problem (2): Find $(\boldsymbol{\sigma}_h,\mathbf{u}_h)\in\mathbf{H}_{h,k}\times\mathbf{Q}_{h,k}$ such that

for all $(\boldsymbol{\tau}_h,\mathbf{v}_h)\in\mathbf{H}_{h,k}\times\mathbf{Q}_{h,k}$, where $A_N:(\mathbf{H}_{h,k}\times\mathbf{Q}_{h,k})\times(\mathbf{H}_{h,k}\times\mathbf{Q}_{h,k})\longrightarrow\mathbb{R}$ is the bilinear form defined by

for all $(\boldsymbol{\sigma}_h,\mathbf{u}_h),(\boldsymbol{\tau}_h,\mathbf{v}_h)\in\mathbf{H}_{h,k}\times\mathbf{Q}_{h,k}$, and $F:\mathbf{H}_{h,k}\times\mathbf{Q}_{h,k}\longrightarrow\mathbb{R}$ is the linear form defined by

for all $(\boldsymbol{\tau}_h,\mathbf{v}_h)\in\mathbf{H}_{h,k}\times\mathbf{Q}_{h,k}$, with $\kappa_1$, $\kappa_2$ and $\kappa_3$ positive parameters that satisfy the conditions (8) and $\gamma_N>0$ being the Nitsche parameter.

The following result guarantees the existence and uniqueness of a solution to the discrete problem (15).

Proof. Using the definition of the bilinear form $A_N(\cdot,\cdot)$ and the norm $\lVert\!(\cdot,\cdot)\!\rVert$ (cf. (16) and (14), respectively), together with Lemma 3, there exists a positive constant $C_A := \min\left\{\frac{\kappa_2}{2\theta},\theta\right\}$ such that

for all $(\tau_h,v_h) \in H_{h,k}\times Q_{h,k}$, which proves the result.

Remark 3. In the case of homogeneous Dirichlet data, $g_D=0$, we obtain the stability estimate

where the constant $C$ depends on $\theta$ and $\lVert\alpha\rVert_{\infty,\Omega}$, but is independent of the viscosity $\nu$. For nonhomogeneous Dirichlet data, the present theoretical framework does not provide a stability estimate for the Nitsche scheme. Establishing such an estimate requires additional arguments and may therefore constitute an interesting direction for future work.

Error analysis

We begin by establishing the consistency of the discrete formulation (15).

Proof. The claim follows directly from (15), integration by parts, and the extra regularity assumptions of $(\sigma,u)\in H^s(\Omega)^{n\times n}\cap H\mathbin{\times}H^{1+s}(\Omega)^n\cap Q$, which satisfies (2).

Now, we introduce some interpolation operators, together with their corresponding approximation properties, which will be instrumental in establishing the optimal convergence of the proposed method. For each $T\in\mathcal{T}_h$, we define $P_T^k:L^2(T)^n\longrightarrow\mathbb{P}_k(T)^n$ as the $L^2$-orthogonal projector. Then, the following approximation estimates hold (see [33 Theorem 18.16]):

for all $v\in H^1(T)^n$. To estimate the pseudostress error, we consider, for $s\in(1/2,k+1]$, the global and local Raviart–Thomas interpolation operators $\Pi_h^k:H^s(\Omega)^{n\times n}\cap H\longrightarrow H_{h,k}$ and $\Pi_T^k(\tau):=\Pi_h^k(\tau)|_T$, respectively, characterized by the following identities (for details, see [20, 38]):

and

From (20), it follows that

for all edge/face $E$ of the element $T \in\mathcal{T}_h$. Moreover, for $s \in(0,k+1]$ the interpolation operator satisfies the following estimate (see [37] for details):

for all $\boldsymbol{\tau}\in H^s(\Omega)^{n\times n}\cap\mathbb{H}$ with $\nabla\cdot\boldsymbol{\tau}\in H^s(\Omega)^n$.

The following result establishes a new interpolation error estimate for $\Pi_h^k$ on the boundary $\Gamma_N$, which will be used to control the Nitsche term in the convergence analysis of the scheme.

Proof. We begin with the case $s\in(1/2,k+1]$. For $E\in\mathcal{E}_N$, let $P_E^k:L^2(E)^n\longrightarrow\mathbb{P}_k(E)^n$ be the $L^2$-orthogonal projector. Then, by the definition of $P_E^k$ and the characterization (20) of $\Pi_h^k$ on each $E\in\mathcal{E}_N$, we have

Now, let $T\in\mathcal{T}_h$ be the unique element such that $E\in\mathcal{E}(T)\cap\mathcal{E}_N$. We define $Q_0:=\frac{1}{|T|}\int_T\boldsymbol{\tau}$. It is clear that $Q_0\boldsymbol{n}\in\mathbb{P}_0(E)^n\subseteq\mathbb{P}_k(E)^n$. Hence, by the best-approximation property of the orthogonal projector $P_E^k$ and (25), we have

To proceed with the proof, we recall the following fractional trace inequality [3 Lemma 2.1]: For all $\boldsymbol{v}\in H^s(\Omega)^n$, we have

and the fractional Poincaré inequality (for details, see [33 Lemma 3.26] or [32 Lemma 7.1]): For all $\boldsymbol{\tau}\in H^s(\Omega)^{n\times n}$, we get

Combining these estimates with (26), we obtain

The case $s=1$ is a direct consequence of estimate [37 (25)]. For the remaining case, when $s\in(1/2,k+1]$, let $r\in\mathbb{N}$ be such that $1\leq r\leq k$ and $r<s\leq r+1$. By [44 Lemma 3.7], there exists a positive constant $C$ and a matrix $Q_r\in\mathbb{P}_r(T)^{n\times n}$, such that

If define $\overline{Q}_r:=Q_r+\frac{1}{|T|}\int_T(\boldsymbol{\tau}-Q_r)$, we have $\int_T\overline{Q}_r=\int_T\boldsymbol{\tau}$, and therefore $\overline{Q}_r-\boldsymbol{\tau}\in L_0^2(T)^{n\times n}$. Next, as $h_T<C$, by Poincaré inequality on $L_0^2(T)^{n\times n}$ (cf. [33 Lemma 3.27]) and (27), we have

In turn, since $r \le k$ and $\overline{Q}_r \in\mathbb{P}_r(T)^{n\times n} \subset\mathbb{RT}_k(T)$, we have $\Pi_h^k\overline{Q}_r=\overline{Q}_r$. Finally, applying the result for $s=1$ together with (4.10) and (4.11), we obtain

which completes the proof.

Finally, for the velocity approximations, we consider the Lagrange interpolant $I_h^k:H^s(\Omega)^n\longrightarrow Q_{h,k}$, which satisfies: For all $T\in\mathcal{T}_h$

for all $v\in H^s(T)^n$ with $l=0,1$ and $\frac{n}{2}<s\leq k+2$ (for details, see [20, 21, 33]).

Theorem 8 provides an optimal a priori error estimate in the energy norm (3.3), with a constant $C$ that is independent of $h$. With the aim of making the dependence on the physical parameters at least partially explicit, the following result is established under assumptions on the parameters of problem $(P)$ that are typically satisfied in practice, namely:

Proof. We proceed as in the proof of [27 Theorem 5]. To this end, we define $\eta^\sigma:=\sigma-\Pi_h^k(\sigma)$ and $\eta^u:=u-I_h^k(u)$. Recalling the definition of the bilinear form $A_N(\cdot,\cdot)$ (cf. (3.5)) and applying (4.5) together with the Cauchy–Schwarz inequality, we obtain

Now, using (2.4), Young’s inequality, and a trace inequality, we obtain

for all $(\tau_h,v_h)\in H_{h,k}\times Q_{h,k}$.

In turn, combining the local trace inequality (cf. [33 Lemma 12.15])

and Young’s inequality with the mesh regularity and (4.1)–(4.2), we deduce

Thus, using these last estimates, we obtain

On the other hand, from the definition of the deviator $(\eta^\sigma)^d$, we get

Then, applying Lemma 7, (4.6), (4.12), (4.15), (4.18), and employing again the mesh regularity, we deduce that

Connecting (4.14), (4.17) and (4.19), we get

for all $(\tau_h,v_h)\in H_{h,k}\times Q_{h,k}$.

On the other hand, from (23), (29), the mesh regularity, and the local and global trace inequalities, we have

Now, defining $(\xi_h^\sigma,\xi_h^u):=(\Pi_h^k(\sigma)-\sigma_h,I_h^k(u)-u_h)$, and using the coercivity of $A_N(\cdot,\cdot)$ (cf. (17)) with $C_A^{-1}=\frac{2\theta}{\kappa_2^2}$, and Lemma 6, we obtain

and applying (37), we conclude that

The proof is completed noting that

and then combining the above estimate with (38) and (39).

Remark 4. An important difference between Theorem 8 and the corresponding result in [27] lies in the regularity assumptions imposed on the exact solution $(\sigma,u)$. More precisely, thanks to Lemma 7, the stronger assumption $(\sigma,u)\in H^{k+1}(\Omega)^{n\times n}\cap H\times H^{k+2}(\Omega)^n\cap Q$ required in [27], can be substantially relaxed, leading to weaker regularity requirements that are essentially minimal for the present analysis.

Numerical experiments

In this section, we present a series of numerical experiments in two and three dimensions to assess the performance of the proposed method and confirm the convergence rates established in Theorem 8. The numerical scheme was implemented using the open-source finite element library FEniCS [2], and the resulting linear systems were solved using the direct solver MUMPS [4]. The numerical solutions in two dimensions were visualized using ParaView [1], whereas the convergence curves were generated using Gnuplot [53]. Moreover, in accordance with Lemma 3, we use the following parameters in all numerical examples:

Robustness to non-solenoidal convection

The aim of this example is to assess the performance of the proposed numerical method when the convective field is not solenoidal. To this end, we consider the following manufactured solution of problem (P):

on the domain $\Omega:= (0,1) \times(0,1)$. Dirichlet and Neumann boundary conditions are prescribed on $\Gamma_D := \{(x,y) \in\partial\Omega: x = 0 \vee y = 0\}$ and $\Gamma_N := \partial\Omega\setminus\Gamma_D$, respectively. The parameters of the generalized Oseen equations are set to $\theta= 1$, $\nu= 10^{-1}$, and the convective field is defined by $\boldsymbol{\alpha}(x,y) := \boldsymbol{u}(x,y) + \varepsilon(x/4,y/4)$, where $\varepsilon> 0$ is the parameter appearing in assumption (H1). Since $\nabla\cdot\boldsymbol{u} = 0$, we have $\nabla\cdot\boldsymbol{\alpha} = \varepsilon/2 < \varepsilon$, and hence (H1) is satisfied. Finally, the boundary data $\boldsymbol{g}_D$ and $\boldsymbol{g}_N$, together with the source term $\boldsymbol{f}$, are determined from the exact solution given in (40).

Figure 1 illustrates the convergence behavior of the proposed method on a sequence of uniformly refined meshes for $\varepsilon\in\{1, 10^{-3}, 10^{-6}, 0\}$ with $k = 0$ and $k = 1$. We observe that the numerical performance remains essentially unaffected as $\varepsilon$ decreases towards zero, corresponding to the transition from non-solenoidal to solenoidal convective fields. Moreover, the observed convergence rates are in agreement with the theoretical predictions of Theorem 8. Similarly, Figure 2 shows the convergence behavior of the individual errors in the pseudostress, velocity, and pressure, for which the expected convergence rates are observed across all values of $\varepsilon$, providing further numerical evidence of the robustness of the proposed method with respect to non-solenoidal convection.

An example with boundary layers

Inspired by the numerical example presented in [48], we consider the domain $\Omega:= (0,1) \times(0,1)$, with parameters $\boldsymbol{\alpha} = (1,1)$, $\nu= 10^{-2}$ and $\theta= 1$. The source term $\boldsymbol{f}$ is chosen so that the exact solution of problem (P) reads

where

with $\lambda= \frac{1}{2\nu}$. We impose homogeneous Dirichlet boundary conditions, $\boldsymbol{g}_D = 0$, on $\Gamma_D := \{(x,y) \in\partial\Omega: x = 0 \vee y = 0\}$, while the Neumann datum $\boldsymbol{g}_N$ on $\Gamma_N := \partial\Omega\setminus\Gamma_D$ is determined from the exact solution given in (41).

The first component of the velocity field $\boldsymbol{u}$ exhibits a boundary layer along $y = 1$, particularly near the corner $(1,1)$, whereas the second component exhibits a similar behavior along $x = 1$. Figure 3 illustrates the convergence behavior of the error in the energy norm $\lVert\!\lVert\!\lVert(\boldsymbol{\sigma} - \boldsymbol{\sigma}_h,\boldsymbol{u} - \boldsymbol{u}_h)\rVert\!\rVert\!\rVert$, for $k = 0$ and $k = 1$. As expected, the presence of boundary layers affects the convergence behavior on coarse meshes. However, as the mesh is sufficiently refined, the expected convergence rates are recovered, in agreement with the

theoretical predictions of Theorem 8. A similar behavior is observed in Figure 4 for the individual errors in the pseudostress tensor, velocity, and pressure, with the expected convergence rates recovered for both $k = 0$ and $k = 1$. Figure 5 displays the isovalues of the computed velocity magnitude $\lvert\mathbf{u}_h \rvert$ and pressure $p_h$, whereas Figure 6 shows the four components of the pseudostress tensor $\boldsymbol{\sigma}_h$. All results were obtained with $k = 1$ on a mesh consisting of 1,048,576 elements. Only the numerical solutions are presented, since their isovalues are visually indistinguishable from those of the corresponding exact solutions.

A manufactured solution in three dimensions

The final experiment assesses the performance of the proposed scheme for a three-dimensional problem with a non-solenoidal convective field. To this end, we consider the manufactured solution of problem (P) given by $\mathbf{u}(x,y,z) = e^{-y}(xz,2z,z^2/2)$ and $p(x,y,z) = (x - 0.5)(y - 0.5)z^2$, defined on the unit cube $\Omega:= (0,1)^3$. We set $\theta= 1.2$, $\nu= 10^{-2}$, and define the convective field by $\boldsymbol{\alpha} = \mathbf{u} + 10^{-3}(0,0,z)$. Although $\nabla\cdot\mathbf{u} = 0$, the convective field satisfies $\nabla\cdot\boldsymbol{\alpha} = 10^{-3}$, and hence $\lVert\nabla\cdot\boldsymbol{\alpha} \rVert_{\infty,\Omega} = 10^{-3}$. Accordingly, we choose $\varepsilon= 1.001 \times10^{-3}$ to satisfy assumption (H1). Furthermore, Dirichlet and Neumann boundary conditions are prescribed on $\Gamma_D := \partial\Omega\setminus\Gamma_N$ and $\Gamma_N := \{(x,y,z) \in\partial\Omega: y = 1\}$, respectively. The boundary data $\mathbf{g}_D$ and $\mathbf{g}_N$, as well as the source term $\mathbf{f}$, are determined from the exact solution $(\mathbf{u},p)$.

Figure 7 illustrates the optimal convergence rates in the energy norm established in Theorem 8, for $k = 0$ and $k = 1$, on a sequence of tetrahedral meshes with decreasing mesh size $h$. Similarly, Figure 8 shows the convergence behavior of the individual errors in the pseudostress tensor, velocity, and pressure, for which the expected convergence rates are observed.

Conclusions

In this work, we have proposed and analyzed a new skew-symmetric Nitsche scheme based on an augmented dual-mixed finite element formulation for the generalized Oseen equations with mixed boundary conditions and non-solenoidal convective fields. By introducing the pseudostress tensor as a new unknown, the classical velocity-pressure formulation is recast as a dual-mixed pseudostress-velocity formulation, in which the boundary conditions are treated differently from those in the primal formulation. In particular, the Dirichlet boundary condition for the velocity is naturally incorporated through integration by parts, whereas the prescribed normal trace of the pseudostress tensor is weakly enforced through consistent Nitsche terms. We have established the well-posedness of the continuous formulation by means of a suitable lifting of the Neumann datum and proved the existence and uniqueness of the discrete solution under appropriate assumptions on the augmentation parameters, without requiring a sufficiently large Nitsche stabilization parameter. Furthermore, the stability of the discrete scheme has been established for homogeneous Dirichlet boundary conditions and general Neumann data. We have also derived optimal a priori error estimates in

an energy norm under relatively mild regularity assumptions on the exact solution, thereby extending the analysis to non-solenoidal convective fields. The theoretical results have been supported by numerical experiments in two and three dimensions, implemented using the open-source finite element library FEniCS. The numerical results confirm the theoretical convergence rates and illustrate the performance of the proposed scheme for non-solenoidal convective fields and solutions exhibiting boundary layers. As future research directions, we envisage extending the present framework to the Navier–Stokes equations and developing adaptive finite element procedures based on a posteriori error estimators, with the aim of improving the local accuracy of the numerical approximation while controlling the computational cost.

Acknowledgements

The second author was partially supported by Dirección de Investigación of the Universidad Católica de la Santísima Concepción through project DIREG 01/2025. The third author was partially supported by ANID-Chile through grant FONDECYT de Iniciación No 11241197.