Abstract
We introduce a family of mixed methods and discontinuous Galerkin discretisations designed to numerically solve the Oseen equations written in terms of velocity, vorticity, and Bernoulli pressure. The unique solvability of the continuous problem is addressed by invoking a global inf-sup property in an adequate abstract setting for non-symmetric systems. The proposed finite element schemes, which produce exactly divergence-free discrete velocities, are shown to be well-defined and optimal convergence rates are derived in suitable norms. This mixed finite element method is also pressure-robust. In addition, we establish optimal rates of convergence for a class of discontinuous Galerkin schemes, which employ stabilisation. A set of numerical examples serves to illustrate salient features of these methods.
Similar content being viewed by others
1 Introduction
The Oseen equations stem from linearisation of the steady (or alternatively from the backward Euler time-discretisation of the transient) Navier–Stokes equations. Of particular appeal to us is their formulation in terms of fluid velocity, vorticity vector, and pressure. A diversity of discretisation methods is available to solve incompressible flow problems using these three fields as principal unknowns. Some recent examples include spectral elements [3, 9, 30], cell-based pressure schemes [12], as well as stabilised and least-squares schemes [2, 10] for Navier–Stokes; also several mixed and augmented methods for Brinkman [4,5,6, 8], and a number of other discretisations specifically designed for Stokes flows [7, 24, 25, 27, 34].
Both the implementation and the analysis of numerical schemes for Navier–Stokes are typically based on the Oseen linearisation. A few related contributions (not only restricted to the velocity–pressure formulation) include for instance [11], that presents a least-squares method for Navier–Stokes equations with vorticity-based first-order reformulation, and whose analysis appeals to the elliptic theory of Agmon–Douglas–Nirenberg, and the proposed conforming finite element methods exhibit optimal order of accuracy for diverse boundary conditions. We also mention the non-conforming exponentially accurate least-squares spectral method for Oseen equations proposed in [32], where a suitable preconditioner is also proposed. In [36] the authors introduce a velocity–vorticity-pressure least-squares finite element method for Oseen and Navier–Stokes equations with velocity boundary conditions. They derive error estimates and reported a degeneracy of the convergence for large Reynolds numbers. A least-squares minimisation problem based on the stress–velocity–pressure formulation was introduced in [15]. The study shows that the corresponding homogeneous least-squares functional is elliptic and continuous in suitable norms. Several first-order Oseen-type systems are analysed in [17], also including vorticity and pressure in the formulation.
Discontinuous Galerkin (DG) methods have also been used to solve the Oseen problem, as for example, in [20, 21] for the case of Dirichlet boundary conditions. Compared with conforming finite elements, discretisations based on DG methods have a number of attractive, and well-documented features. These include high order accuracy, being amenable for hp-adaptivity, relatively simple implementation on highly unstructured meshes, and superior robustness when handling rough coefficients. We also mention the a priori error analysis of hybridisable DG schemes introduced in [18] for the Oseen equations. The family of DG methods we propose here has resemblance with those schemes, but concentrates on a three-field formulation described below.
This paper is concerned with mixed non-symmetric variational problems which will be analysed using a global inf-sup argument. To do this, we conveniently restrict the set of equations to the space of divergence-free velocities, and apply results from [26, Sect. 2] in order to prove that the equivalent resulting non-symmetric saddle-point problem is well-posed. For the numerical approximation, we first consider Raviart–Thomas elements of order \(k\ge 0\) for the velocity field, Nédélec elements of order k for the vorticity, and piecewise polynomials of degree k without continuity restrictions, for the Bernoulli pressure. We prove unique solvability of the discrete problem by adapting the same tools utilised in the analysis of the continuous problem. In addition, the proposed family of Galerkin finite element methods turns out to be optimally convergent, under the common assumptions of enough regularity of the exact solutions to the continuous problem. An appealing feature of this mixed finite element method is that it produces exactly divergence-free approximations of the velocity by construction; thus preserving, at the discrete level, an essential constraint of the governing equations. Next, inspired by the methods presented in [21, 22], we present another scheme involving the discontinuous Galerkin discretisation of the \(\mathop {\mathbf {curl}}\nolimits \)-\(\mathop {\mathbf {curl}}\nolimits \) and \(\mathrm {grad}\)-\(\mathrm {div}\) operators. An advantage is related to the robustness with respect to rough coefficients and the relaxation of inter-element continuity conditions. We prove the well-posedness of the DG scheme and derive error estimates under appropriate regularity assumptions. The novelties of our contribution are mainly related to the formulation of Oseen equations in terms of vorticity, in considering H(div)-conforming velocities, in proposing mixed finite element and discontinuous Galerkin schemes based on such formulation, and in establishing their convergence properties.
We have structured the contents of the paper in the following manner. Notation-related preliminaries are stated in the remainder of this section. We then present the model problem as well as the three-field weak formulation and its solvability analysis in Sect. 2. The finite element discretisation is constructed in Sect. 3, where we also derive the stability and convergence bounds. In Sect. 4, we present the mixed DG formulation for the model problem. The well-posedness of the method and the error analysis are established in the same section. We close in Sect. 5 with a set of computational tests that illustrate the properties of the proposed numerical schemes in a variety of scenarios.
Let \(\varOmega \) be a bounded domain of \(\mathbb {R}^3\) with Lipschitz boundary \(\partial \varOmega \). Moreover, we assume that \(\partial \varOmega \) admits a disjoint partition \(\partial \varOmega =\varGamma \cup \varSigma \). For any \(s\ge 0\), the symbol \(\left\| \cdot \right\| _{s,\varOmega }\) denotes the norm of the Hilbertian Sobolev spaces \(\mathrm {H}^s(\varOmega )\) or \(\mathrm {H}^s(\varOmega )^3\), with the usual convention \(\mathrm {H}^0(\varOmega ):=\mathrm {L}^2(\varOmega )\). For \(s\ge 0\), we recall the definition of the space
that we will endow with the norm \(\left\| \varvec{\theta }\right\| ^2_{\mathrm {H}^s(\mathop {\mathbf {curl}}\nolimits ;\varOmega )}=\left\| \varvec{\theta }\right\| _{s,\varOmega }^2+\left\| \mathop {\mathbf {curl}}\nolimits \varvec{\theta }\right\| ^2_{s,\varOmega }\), and will denote \(\mathrm {H}(\mathop {\mathbf {curl}}\nolimits ;\varOmega )=\mathrm {H}^0(\mathop {\mathbf {curl}}\nolimits ;\varOmega )\).
2 Statement and Solvability of the Continuous Problem
2.1 Oseen Problem in Terms of Velocity–Vorticity–Pressure
A standard backward Euler time-discretisation of the classical Navier–Stokes equations, or a linearisation of the steady version of the problem combined with standard curl-div identities, leads to the following set of equations, known as the Oseen equations:
where \(\nu >0\) is the kinematic fluid viscosity, \(\sigma >0\) is inversely proportional to the time-step, \(\varvec{\beta }\) is an adequate approximation of velocity to be made precise below, and the vector of external forces \(\varvec{f}\) also absorbs the contributions related to previous time steps, or to fixed states in the linearisation procedure of the steady Navier–Stokes equations. As usual in this context, in the momentum equation we have conveniently introduced the Bernoulli (also known as dynamic) pressure \(p := P + \frac{1}{2}\vert \varvec{u}\vert ^2\), where P is the actual fluid pressure.
The structure of (2.1) suggests to introduce the rescaled vorticity vector \(\varvec{\omega }:=\sqrt{\nu }\mathop {\mathbf {curl}}\nolimits \varvec{u}\) as a new unknown. Furthermore, in this study we focus on the case of zero normal velocities and zero tangential vorticity trace imposed on a part of the boundary \(\varGamma \subset \partial \varOmega \), whereas a non-homogeneous tangential velocity \(\varvec{u}_{\varSigma }\) and a fixed Bernoulli pressure \(p_\varSigma \) are set on the remainder of the boundary \(\varSigma = \partial \varOmega \setminus \varGamma \). Therefore, system (2.1) can be recast in the form
where \(\varvec{n}\) stands for the outward unit normal on \(\partial \varOmega \). Should the boundary \(\varSigma \) have zero measure, the additional condition \( (p,1)_{\varOmega ,0} = 0\) is required to enforce uniqueness of the Bernoulli pressure. Note that even if the formulation and its analysis are restricted to the case of constant coefficients, the theory might be extended to the case of piecewise constant model coefficients.
2.2 Defining a Weak Formulation
Let us introduce the following functional spaces
where \(\gamma _{t}\) is the tangential trace operator on \(\varGamma \), defined by: \(\gamma _{t}(\varvec{\theta })=\varvec{\theta }\times \varvec{n}\). Let us endow \(\mathrm {H}\) and \(\mathrm {Q}\) with their natural norms. For the space \(\mathrm {Z}\) however, we consider the following viscosity-weighted norm:
From now on, we will assume that the data are regular enough: \(\varvec{\beta }\in \mathrm {L}^{\infty }(\varOmega )^3\) and \(\varvec{f}\in \mathrm {L}^2(\varOmega )^3\). We proceed to test (2.2) against adequate functions and to impose the boundary conditions in such a manner that we end up with the following formulation: Find \((\varvec{u},\varvec{\omega },p)\in \mathrm {H}\times \mathrm {Z}\times \mathrm {Q}\) such that
where the bilinear forms \(a:\mathrm {H}\times \mathrm {H}\rightarrow \mathbb {R}\), \(b_1:\mathrm {H}\times \mathrm {Z}\rightarrow \mathbb {R}\), \(d:\mathrm {Z}\times \mathrm {Z}\rightarrow \mathbb {R}\), \(b_2:\mathrm {H}\times \mathrm {Q}\rightarrow \mathbb {R}\), \(c:\mathrm {Z}\times \mathrm {H}\rightarrow \mathbb {R}\), and the linear functionals \(F:\mathrm {H}\rightarrow \mathbb {R}\), and \(G:\mathrm {Z}\rightarrow \mathbb {R}\) are specified as follows
for all \(\varvec{u},\varvec{v}\in \mathrm {H}\), \(\varvec{\omega },\varvec{\theta }\in \mathrm {Z}\), and \(q\in \mathrm {Q}\).
The continuity of these bilinear and linear functionals is stated in the following lemma, whose proof is obtained by rather standard arguments. In particular, the estimate for \(c(\cdot ,\cdot )\) is proven using the assumption on \(\varvec{\beta }\) and the fact that \(\Vert \varvec{\theta }\times \varvec{\beta }\Vert _{0,\varOmega }\le 2\Vert \varvec{\beta }\Vert _{\infty ,\varOmega }\Vert \varvec{\theta }\Vert _{0,\varOmega }\).
Lemma 1
The following estimates hold true,
where C is a positive constant.
2.3 Solvability Analysis
In order to analyse the variational formulation (2.3), let us introduce the kernel of the bilinear form \(b_2(\cdot ,\cdot )\) and its classical characterisation
and let us recall that \(b_2\) satisfies the inf-sup condition:
with an inf-sup constant \(\beta _2>0\) only depending on \(\varOmega \) (see e.g. [26, Sect. 2.4.1]).
We will now address the well-posedness of (2.3). To that end, it is enough to study its reduced counterpart, defined on \(\mathrm {X}\times \mathrm {Z}\): Find \((\varvec{u},\varvec{\omega })\in \mathrm {X}\times \mathrm {Z}\) such that
The equivalence between (2.3) and (2.12) is established in the following result, whose proof follows [28, Sect. I.4.1] and it is basically a direct consequence of the inf-sup condition (2.11).
Lemma 2
If \((\varvec{u},\varvec{\omega },p)\in \mathrm {H}\times \mathrm {Z}\times \mathrm {Q}\) is a solution of (2.3), then \(\varvec{u}\in \mathrm {X}\) and \((\varvec{u},\varvec{\omega })\in \mathrm {X}\times \mathrm {Z}\) also solve (2.12). Conversely, if \((\varvec{u},\varvec{\omega })\in \mathrm {X}\times \mathrm {Z}\) is a solution of (2.12), then there exists a unique \(p\in \mathrm {Q}\) such that \((\varvec{u},\varvec{\omega },p)\in \mathrm {H}\times \mathrm {Z}\times \mathrm {Q}\) solves (2.3).
The abstract setting that will permit the analysis of (2.3) is stated in the following general result [26, Theorem 1.2].
Theorem 1
Let \({\mathcal {A}}: {\mathcal {X}}\times {\mathcal {X}} \rightarrow \mathbb {R}\) be a bounded bilinear form and \({\mathcal {G}}: {\mathcal {X}}\rightarrow \mathbb {R}\) a bounded functional, both defined on the Hilbert space \(({\mathcal {X}},\langle \cdot ,\cdot \rangle _{{\mathcal {X}}})\). If there exists \(\alpha > 0\) such that
and
then there exists a unique solution \(x \in {\mathcal {X}}\) to the problem
Furthermore, there exists \(C>0\) (independent of x) such that
Lemma 3
Let us assume that
and let us define the bilinear form \({\mathcal {A}}(\cdot ,\cdot )\) specified as
Then, there exist \(\alpha _1,\alpha _2>0\) such that
and
where \({\mathcal {X}}:=\mathrm {X}\times \mathrm {Z}\), endowed with the corresponding product norm, is a Hilbert space.
Proof
As a consequence of Lemma 1, we have that the bilinear form \({\mathcal {A}}(\cdot ,\cdot )\) is bounded and therefore the condition (2.16) readily follows.
Concerning the satisfaction of the inf-sup condition (2.17), for a given \((\varvec{u},\varvec{\omega })\in {\mathcal {X}}\), we can define
where \(\hat{c}>0\) is a constant to be chosen later. We can then immediately assert that
where we have used the bound \(\Vert \varvec{\omega }\times \varvec{\beta }\Vert _{0,\varOmega }\le 2\Vert \varvec{\beta }\Vert _{\infty ,\varOmega }\Vert \varvec{\omega }\Vert _{0,\varOmega }\). Choosing \(\hat{c}=1/(4\sigma )\) and employing (2.15), we arrive at
with C independent of \(\nu \). On the other hand, by construction we realise that \(\Vert {{\tilde{\varvec{\theta }}}}\Vert _{\mathrm {Z}}=\Vert \varvec{\omega }\Vert _{\mathrm {Z}}\) and \(\Vert {\tilde{\varvec{v}}}\Vert _{0,\varOmega }\le C\hat{c}(\Vert \varvec{u}\Vert _{0,\varOmega }+\Vert \varvec{\omega }\Vert _{\mathrm {Z}})\), and consequently
which finishes the proof. \(\square \)
Lemma 4
Suppose that the bound (2.15) is satisfied. Then,
Proof
For all \((\varvec{v},\varvec{\theta })\in {\mathcal {X}}\), we have that:
\(\square \)
As a consequence of the previous lemmas, we have the following result.
Theorem 2
Let us assume (2.15). Then, the variational problem (2.12) admits a unique solution \((\varvec{u},\varvec{\omega })\in \mathrm {X}\times \mathrm {Z}\). Moreover, there exists \(C>0\) such that
Proof
It suffices to verify the hypotheses of Theorem 1. First, we define the linear functional
which is bounded on \(\mathrm {X}\times \mathrm {Z}\) as a consequence of the last two estimates in Lemma 1. Thus, the proof follows from Lemmas 3 and 4. \(\square \)
The following result establishes the corresponding stability estimate for the Bernoulli pressure.
Corollary 1
Let \((\varvec{u},\varvec{\omega })\in \mathrm {X}\times \mathrm {Z}\), be the unique solution of (2.12), with \(\varvec{u}\) and \(\varvec{\omega }\) satisfying (2.18). In addition, let \(p\in Q\) be the unique pressure provided by Lemma 2, so that \((\varvec{u},\varvec{\omega },p)\in \mathrm {H}\times \mathrm {Z}\times \mathrm {Q}\) is the unique solution of (2.3). Then, there exists \(C>0\) such that
Proof
Combining the inf-sup condition (2.11) with the first equation in (2.3) gives the bound
which together with (2.18) and Lemma 1, complete the proof.\(\square \)
Remark 1
An alternative analysis for the nonsymmetric variational problem (2.3) can be carried out using a fixed-point argument that allows a symmetrisation of the mixed structure. The resulting weak form could then be analysed using classical tools for saddle-point problems, for instance, following the similar treatment in [16, Sect. 3.3]. Establishing inf-sup conditions for the off-diagonal bilinear forms in the original nonsymmetric formulation is, however, much more involved (see e.g. [33, Sect. 3]).
Remark 2
Assumption (2.15) holds provided one chooses \(\sigma \) appropriately. As this parameter represents the inverse of the timestep, the aforementioned relation constitutes then a CFL-type condition. We have observed experimentally that the bound for \(\varvec{\beta }\) is not sharp, but we also stress that the same condition coincides with the hypotheses that yield solvability of least-squares formulations for the Oseen problem analysed in [15, 17]. Further investigations would be then required.
3 Finite Element Discretisation
In this section we introduce a Galerkin scheme for (2.3) and analyse its well-posedness by establishing suitable assumptions on the finite element subspaces involved. Error estimates are also derived.
3.1 Defining the Discrete Problem
Let \(\{{\mathcal {T}}_{h}(\varOmega )\}_{h>0}\) be a shape-regular family of partitions of the polyhedral region \({\bar{\varOmega }}\), by tetrahedrons T of diameter \(h_T\), with mesh size \(h:=\max \{h_T:\; T\in {\mathcal {T}}_{h}(\varOmega )\}\). In what follows, given an integer \(k\ge 0\) and a subset S of \(\mathbb {R}^3\), \({\mathcal {P}}_k(S)\) will denote the space of polynomial functions defined locally in S and being of total degree \(\le k\).
Moreover, for any \(T\in {\mathcal {T}}_{h}(\varOmega )\), we introduce the local Nédélec space
where \(R_{k+1}(T)\) is a subspace of \({\mathcal {P}}_{k+1}(T)^3\) composed by homogeneous polynomials of degree \(k+1\), and being orthogonal to \(\varvec{x}\). With these tools, let us define the following finite element subspaces:
where \(RT_k(T)={\mathcal {P}}_{k}(T)^3\oplus {\mathcal {P}}_{k}(T)\varvec{x}\) is the Raviart–Thomas space defined locally in \(T\in {\mathcal {T}}_h(\varOmega )\).
The proposed Galerkin scheme approximating (2.3) reads as follows: Find \((\varvec{u}_h,\varvec{\omega }_h,p_h)\in \mathrm {H}_h\times \mathrm {Z}_h\times \mathrm {Q}_h\) such that
3.2 Solvability and Stability of the Discrete Problem
The analysis of the Galerkin formulation will follow the same arguments exploited in the continuous setting. Let us then consider the discrete kernel of \(b_2\):
where the characterisation is indeed possible thanks to the inclusion \(\mathop {\mathrm {div}}\nolimits \mathrm {H}_h\subseteq \mathrm {Q}_h\). Moreover, it is well-known that the following discrete inf-sup condition holds (see [26, Sect. 4.2]):
We again resort to a reduced version of the problem, now defined on the product space \(\mathrm {X}_h\times \mathrm {Z}_h\). Find \((\varvec{u}_h,\varvec{\omega }_h)\in \mathrm {X}_h\times \mathrm {Z}_h\) such that
and its equivalence with (3.4) is once more a direct consequence of the inf-sup condition (3.6).
Lemma 5
If \((\varvec{u}_h,\varvec{\omega }_h,p_h)\in \mathrm {H}_h\times \mathrm {Z}_h\times \mathrm {Q}_h\) is a solution of (3.4), then \(\varvec{u}_h\in \mathrm {X}_h\), and \((\varvec{u}_h,\varvec{\omega }_h)\in \mathrm {X}_h\times \mathrm {Z}_h\) is also a solution of (3.7). Conversely, if \((\varvec{u}_h,\varvec{\omega }_h)\in \mathrm {X}_h\times \mathrm {Z}_h\) is a solution of (3.7), then there exists a unique \(p_h\in \mathrm {Q}_h\) such that \((\varvec{u}_h,\varvec{\omega }_h,p_h)\in \mathrm {H}_h\times \mathrm {Z}_h\times \mathrm {Q}_h\) is a solution of (3.4).
In order to establish the well-posedness of (3.7), we will employ the following discrete version of Theorem 1.
Theorem 3
Assume (2.15). Let \(k\ge 0\) be an integer and let \(\mathrm {X}_h\) and \(\mathrm {Z}_h\) be given by (3.5) and (3.1), respectively. Then, there exists a unique \((\varvec{u}_h,\varvec{\omega }_h)\in \mathrm {X}_h\times \mathrm {Z}_h\) solution of the discrete scheme (3.7). Moreover, there exist positive constants \(\hat{C}_1,\,\hat{C}_2>0\) independent of h such that
and
where \((\varvec{u},\varvec{\omega })\in \mathrm {X}\times \mathrm {Z}\) is the unique solution to (2.12).
Proof
Let us define \({\mathcal {X}}_h:=\mathrm {X}_h\times \mathrm {Z}_h\) and reuse the forms \({\mathcal {A}}(\cdot ,\cdot )\) and \({\mathcal {G}}(\cdot )\) as in the proof of Lemma 2. The next step consists in proving that \({\mathcal {A}}(\cdot ,\cdot )\) satisfies the discrete version of the inf-sup conditions (2.13)–(2.14), as in Lemmas 3 and 4. In order to assert (2.13), we consider \((\varvec{u}_h,\varvec{\omega }_h)\in {\mathcal {X}}_h\), and define
Then, repeating exactly the same steps used in the proof of Lemma 3 the discrete version of (2.13) follows. Regarding the discrete version of (2.14), we once again repeat the same arguments given in the proof of Lemma 4. Finally, the Céa estimate follows from classical arguments. \(\square \)
Remark 3
Note that since the discrete kernel satisfies \(\mathrm {X}_h\subset \mathrm {X}\), the infimum in the Céa estimate (3.9) can be taken over \((\varvec{v}_h,\varvec{\theta }_h)\in \mathrm {H}_h\times \mathrm {Z}_h\), which implies in particular that the mixed finite element scheme does not depend on the pressure (see e.g. [26, pp. 58–59]). Implications of pressure-robust schemes as this one are reviewed in e.g. [29].
We now state the stability and an adequate approximation property of the discrete pressure.
Corollary 2
Let \((\varvec{u}_h,\varvec{\omega }_h)\in \mathrm {X}_h\times \mathrm {Z}_h\) be the unique solution of (3.7), with \(\varvec{u}_h\) and \(\varvec{\omega }_h\) satisfying (3.8). In addition, let \(p_h\in Q_h\) be the unique discrete Bernoulli pressure provided by Lemma 5, so that \((\varvec{u}_h,\varvec{\omega }_h,p_h)\in \mathrm {H}_h\times \mathrm {Z}_h\times \mathrm {Q}_h\) is the unique solution of (3.4). Then, there exist positive constants \(\bar{C}_1,\,{\bar{C}}_2>0\), independent of h and \(\nu \), such that
and
Proof
The result follows using the same arguments considered in the proof of Corollary 1, but using the discrete inf-sup condition (3.6). We omit further details. \(\square \)
3.3 A Priori Error Estimates
Let us introduce for a given \(s>1/2\), the Nédeléc global interpolation operator \({{\mathcal {R}}}_{h}:\mathrm {H}^s(\mathop {\mathbf {curl}}\nolimits ;\varOmega )\cap \mathrm {Z}\rightarrow \mathrm {Z}_h\). From [1, Proposition 5.6] we know that for all \(\varvec{\theta }\in \mathrm {H}^{s}(\mathop {\mathbf {curl}}\nolimits ;\varOmega )\) with \(s>1/2\), there exists \(C>0\) independent of h, such that
On the other hand, for the Raviart–Thomas interpolation \(\varPi _h:\mathrm {H}^s(\varOmega )^{3}\cap \mathrm {H}\rightarrow \mathrm {H}_h\), with \(s>0\), we recall (see e.g. [26, Theorem 3.6 and Lemma 3.19]) that there exists \(C>0\), independent of h, such that for all \(s>0\):
Finally we recall that the orthogonal projection from \(\mathrm {L}^2(\varOmega )\) onto the finite element subspace \(\mathrm {Q}_h\), here denoted by \(\varPi _{\mathrm {Q}}\), satisfies the following error estimate for all \(s>0\):
The following result summarises the error analysis for our mixed finite element scheme (3.4).
Theorem 4
Assume (2.15). Let \(k\ge 0\) be an integer and let \(\mathrm {H}_h,\mathrm {Z}_h\) and \(\mathrm {Q}_h\) be given by (3.1), (3.2), and (3.3). Let \((\varvec{u},\varvec{\omega },p)\in \mathrm {H}\times \mathrm {Z}\times \mathrm {Q}\) and \((\varvec{u}_h,\varvec{\omega }_h,p_h)\in \mathrm {H}_h\times \mathrm {Z}_h\times \mathrm {Q}_h\) be the unique solutions to the continuous and discrete problems (2.3) and (3.4), respectively. Assume that \(\varvec{u}\in \mathrm {H}^s(\varOmega )^3\), \(\mathop {\mathrm {div}}\nolimits \varvec{u}\in \mathrm {H}^s(\varOmega )\), \(\varvec{\omega }\in \mathrm {H}^{1+s}(\varOmega )^3\) and \(p\in \mathrm {H}^s(\varOmega )\), for some \(s>1/2\). Then, there exist positive constants \(\hat{C}\) and \({\tilde{C}}\), independent of h, such that
and
Proof
The proof follows from (3.9), (3.10), and standard interpolation estimates satisfied by the operators \({{\mathcal {R}}}_h\), \(\varPi _h\) and \(\varPi _{\mathrm {Q}}\) (see (3.11), (3.12) and (3.13), respectively). \(\square \)
4 Discontinuous Galerkin Method
In this section, we propose and analyse a DG method for (2.2). We provide solvability and stability of the discrete scheme by introducing suitable numerical fluxes. A priori error estimates are also derived.
4.1 Preliminaries
Apart from the definitions laid out at the beginning of Sect. 3, let us denote by \({{\mathcal {E}}}_h\) the set of internal faces, by \({{\mathcal {F}}}_h^{\varSigma }\) the set of external faces on \(\varSigma \) and by \({{\mathcal {F}}}_h^{\varGamma }\) the set of external faces on \(\varGamma \). We set \({{\mathcal {F}}}_h={{\mathcal {E}}}_h \cup {{\mathcal {F}}}_h^{\varSigma }\cup {{\mathcal {F}}}_h^{\varGamma }\). We denote by \(h_e\) the diameter of each face e. Let \(T^+\) and \(T^-\) be two adjacent elements of \({\mathcal {T}}_h\) and let \(\varvec{n}^+\) (respectively \(\varvec{n}^-\)) be the outward unit normal vector on \(\partial T^+\) (respectively \(\partial T^-\)). For a vector field \(\varvec{u}\), we denote by \(\varvec{u}^\pm \) the trace of \(\varvec{u}\) from the interior of \(T^\pm \). We define jumps
and averages
and adopt the convention that for boundary faces \(e\in {{\mathcal {F}}}_h^\varSigma \cup {{\mathcal {F}}}_h^\varGamma \), we set \({\left[ \!\left[ \varvec{v}\right] \!\right] }_T=\varvec{v}\times \varvec{n}\), \({\left[ \!\left[ \varvec{v}\right] \!\right] }_N=\varvec{v}\cdot \varvec{n}\), \({[\![q]\!]}=q\varvec{n}\), \(\{\!\{\varvec{v}\}\!\}=\varvec{v}\) and \(\{\!\{q\}\!\}=q\).
Now, for \(k\ge 0\), suitable finite dimensional spaces for vorticity and velocity that remove the restriction of continuity are defined by:
and we remark that the space for pressure approximation will coincide with the one used in Sect. 3, that is \({\tilde{\mathrm {Q}}}_h := \mathrm {Q}_h\).
4.2 Discrete Formulation and Solvability Analysis
Multiplying each equation in (2.2) by suitable functions, the resulting DG scheme consists in finding \((\varvec{u}_h,\, \varvec{\omega }_h,\,p_h)\in {\tilde{\mathrm {H}}}_h \times {{\tilde{Z}}}_h\times {\tilde{\mathrm {Q}}}_h\), such that for any test functions \((\varvec{v}_h,\, \varvec{\theta }_h,\,q_h)\in {\tilde{\mathrm {H}}}_h \times {{\tilde{Z}}}_h\times {\tilde{\mathrm {Q}}}_h\) and for all elements T in the partition \({\mathcal {T}}_h\)
where \(\widehat{\varvec{u}}_h^w,\, \widehat{\varvec{u}}_h^p, \,\widehat{\varvec{\omega }}_h\) and \({{\widehat{p}}}_h\) are numerical fluxes, which approximate the traces of \(\varvec{u}_h\), \(\varvec{\omega }_h\) and \(p_h\) on the boundary. The fluxes \(\widehat{\varvec{\omega }}_h\) and \(\widehat{\varvec{u}}_h^w\) are related to the \(\mathop {\mathbf {curl}}\nolimits \)-\(\mathop {\mathbf {curl}}\nolimits \) operator and are defined by
whereas the fluxes \(\widehat{\varvec{u}}_h^p\) and \({{\widehat{p}}}_h\) are associated with the \(\mathop {\mathrm {grad}}\nolimits \)-\(\mathop {\mathrm {div}}\nolimits \) operator and defined by
The parameters \(C_{11}\), \(A_{11}\) and \(D_{11}\) are positive stabilisation parameters, and following [22] we choose
where \(c_{11},d_{11},a_{11}>0\). Moreover, we suppose that \(C_{11}\) (respectively \(D_{11}\) and \(A_{11}\)) have a uniform positive bound above and below denoted by \(\overline{C_{11}}\) and \(\underline{C_{11}}\) (respectively \(\overline{D_{11}}\), \(\underline{D_{11}}\) and \(\overline{A_{11}}\), \(\underline{A_{11}}\)).
We then proceed to integrate by parts equations (4.1) and (4.3), and then summing up over all \(T\in {\mathcal {T}}_h\), we obtain the following DG scheme: Find \((\varvec{u}_h,\varvec{\omega }_h,p_h)\in {{\widetilde{\mathrm {H}}}}_h\times {\widetilde{\mathrm {Z}}}_h\times {{\widetilde{\mathrm {Q}}}}_h\) such that
where the forms a, c and d are the same in (2.3), while \({{\tilde{b}}}_1\), \({{\tilde{b}}}_2\), j and e are defined, respectively, by:
In addition, the linear functionals \({\widetilde{F}}\), \({\widetilde{G}}\) and \({\widetilde{L}}\) associated with the source terms are defined as:
By integration by parts and as a consequence of the identity:
it follows that the form \({{\tilde{b}}}_{1}\) can be written as:
Similarly, using the following identity:
the form \({{\tilde{b}}}_{2}\) can be recast, after integration by parts, as follows
To simplify the exposition of the analysis of the method, we will write the mixed scheme (4.9) in the following equivalent form: Find \((\varvec{u}_h,\varvec{\omega }_h,p_h)\in {{\widetilde{\mathrm {H}}}}_h\times {\widetilde{\mathrm {Z}}}_h\times {{\widetilde{\mathrm {Q}}}}_h\) such that
where
and
Let us now show the existence and uniqueness of solution to formulation (4.9).
Proposition 1
Let \(k\ge 0\) be an integer. The DG method (4.9) with the numerical fluxes given by (4.4)–(4.5) defines a unique approximate solution \((\varvec{u}_h,\varvec{\omega }_h,p_h)\in {{\widetilde{\mathrm {H}}}}_h\times {\widetilde{\mathrm {Z}}}_h\times {{\widetilde{\mathrm {Q}}}}_h\) provided that
Proof
Since the problem is linear and finite dimensional, it suffices to show that if \(\varvec{f}=\mathbf{0 }\), \(p_{\varSigma }=0\) and \(\varvec{u}_{\varSigma }=\mathbf{0 }\), then \((\varvec{u}_h,\varvec{\omega }_h,p_h)=(\mathbf{0 },\mathbf{0 },0)\). To this end, take \(\varvec{v}_h=\varvec{u}_h\), \(\varvec{\theta }_h=\varvec{\omega }_h\) and \(q_h=p_h\) in (4.9), summing up the three equations, we obtain:
It follows that
where we define
Using Young’s inequality, we can assert that
Therefore, in particular we have that:
which, owing to the assumption (4.13), implies that \(\varvec{u}_h=\mathbf{0 }\), \(\varvec{\omega }_h=\mathbf{0 }\) and \({[\![p_h]\!]}=0\) on \({{\mathcal {E}}}_h\), \(p_h=0\) on \({{\mathcal {F}}}_h^{\varSigma }\). The first equation in (4.9) then becomes
and then \(\nabla p_h=\mathbf{0 }\). Employing this result, together with \({[\![p_h]\!]}=0\) on \({{\mathcal {E}}}_h\), \(p_h=0\) on \({{\mathcal {F}}}_h^{\varSigma }\) or the fact that \(p_h\) has zero mean value if \(\varSigma \) have zero measure, we conclude that \(p_h=0\). \(\square \)
4.3 A Priori Error Bounds
Let us now present and discuss a priori error bounds for the proposed DG method. The proof involves two steps. The first one consists in establishing an error estimate in the natural semi-norm. In the second step, we prove the error estimate for the pressure in the \(\mathrm {L}^2\)-norm. For this, we introduce the following semi-norm \(\vert \cdot \vert _{{\mathcal {A}}}\):
For the analysis, we will also employ the following norm
We define \(\varvec{e}_{\varvec{u}}=\varvec{u}-\varvec{u}_h\), \(\varvec{e}_{\varvec{\omega }}=\varvec{\omega }-\varvec{\omega }_h\) and \(e_p=p-p_h\). Let us denote by \(\varvec{\varPi }_{{{\widetilde{\mathrm {H}}}}}\) (respectively \(\varvec{\varPi }_{{\widetilde{\mathrm {Z}}}}\) and \(\varPi _{{{\widetilde{\mathrm {Q}}}}}\)) the \(\mathrm {L}^{2}\)-projection onto \({{\widetilde{\mathrm {H}}}}_h\)\(({\mathcal {P}}_{k+1}^3)\) (respectively \({\widetilde{\mathrm {Z}}}_h\)\(({\mathcal {P}}_{k}^3)\) and \({{\widetilde{\mathrm {Q}}}}_h\)\(({\mathcal {P}}_{k})\)), and let us split the errors in the following manner
where the numerical and approximation errors are defined by:
We recall the following standard approximation properties (see for instance [19]).
Lemma 6
Let \(\varvec{v}\in \mathrm {H}^{1+r}(\varOmega )\), \(r\ge 0\). Let \(\varPi \) the projection operator such that \(\varPi \varvec{v}=\varvec{v}\) for all \(\varvec{v}\in {\mathcal {P}}_{k}(T)\), \(k\ge 0\). Then we have
As a consequence, we have the following result.
Lemma 7
Let \((\varvec{u},\varvec{\omega },p)\in \mathrm {H}\times \mathrm {Z}\times \mathrm {Q}\) be the unique solution to the continuous problem (2.3). Assume that \(\varvec{u}\in \mathrm {H}^{1+s}(\varOmega )^3\), \(\varvec{\omega }\in \mathrm {H}^s(\varOmega )^3\) and \(p\in \mathrm {H}^s(\varOmega )\), for some \(s\ge 1\). Then, we have:
where \(C_{a}\), \(C_{d}\), \(C_{j}\) and \(C_{e}\) are positive constants independent of the meshsize.
Proof
The two first estimates are a simple consequence of (4.15). Next we can state that
and
Recalling that \(C_{11}\), \(A_{11}\) and \(D_{11}\) are as in (4.6)–(4.8) and using (4.16), we end up with the bound
and similarly we obtain
where \(C_j\) (respectively \(C_e\)) depends on \(c_{11} \), \(a_{11}\) and \(\nu \) (respectively \(d_{11}\)). \(\square \)
We next concentrate on obtaining bounds for the forms \({{\tilde{b}}}_{1}\), \({{\tilde{b}}}_{2}\), c and j.
Lemma 8
Let \((\varvec{u},\varvec{\omega },p)\in \mathrm {H}\times \mathrm {Z}\times \mathrm {Q}\) be the unique solution to the continuous problem (2.3). Assume that \(\varvec{u}\in \mathrm {H}^{1+s}(\varOmega )^3\), \(\varvec{\omega }\in \mathrm {H}^s(\varOmega )^3\) and \(p\in \mathrm {H}^s(\varOmega )\), for some \(s\ge 1\). Then one has
where \(C_{b_{1}}\), \(C_{b_{2}}\), \(C_j\), \(C_e\) and \(C_{\varvec{\omega }}\) are positive constants independent of the meshsize.
Proof
The bounds associated with the form \({{\tilde{b}}}_{2}\) can be proved exactly with the same arguments as [22, Sect. 3.3]. Let us now deal with the term \(\tilde{b}_{1}\). We observe that due to the properties of the \(\mathrm {L}^2\)-projection onto \({{{\widetilde{\mathrm {H}}}}_h}\)\(({\mathcal {P}}_{k+1}^3)\), we can write
Using Cauchy–Schwarz’s inequality, we then readily obtain
and the desired estimate follows from the inverse inequality and Lemma 7.
Similarly, using (4.10) and again the properties of the \(\mathrm {L}^{2}\)-projection onto \({{\widetilde{\mathrm {Z}}}_h}\)\(({\mathcal {P}}_{k}^3)\), we have
thus
and then, we simply have to use (4.6) and Lemma 6.
The estimates for the forms j and e are obtained in a similar way. Indeed, using once again Cauchy–Schwarz’s inequality and Lemma 7, it follows that
and proceeding analogously as before, we get
Finally, concerning the term \(c(\varvec{\xi }_{\varvec{\omega }},\varvec{v}_h)\) we can assert that
and exploiting the previous bounds, we obtain the corresponding estimate with \(C_{\varvec{\omega }}=\frac{2C_{d}}{\sqrt{\nu }}\Vert \varvec{\beta }\Vert _{\infty }\). \(\square \)
Theorem 5
Assume (4.13). Let \(k\ge 0\) be an integer and let \((\varvec{u},\varvec{\omega },p)\in \mathrm {H}\times \mathrm {Z}\times \mathrm {Q}\) be the unique solution to the continuous problem (2.3). Assume that \(\varvec{u}\in \mathrm {H}^{1+s}(\varOmega )^3\), \(\varvec{\omega }\in \mathrm {H}^s(\varOmega )^3\) and \(p\in \mathrm {H}^s(\varOmega )\), for some \(s\ge 1\). Then, the mixed DG approximation \((\varvec{u}_h,\varvec{\omega }_h,p_h)\in {{\widetilde{\mathrm {H}}}}_h\times {\widetilde{\mathrm {Z}}}_h\times {{\widetilde{\mathrm {Q}}}}_h\) defined by (4.9), satisfies the following a priori error bounds
where \(C_{{\mathcal {A}}}\) and C are positive constants independent of the meshsize.
Proof
We begin with the estimate (4.17). A direct application of the definition of the \({\mathcal {A}}-\)seminorm in combination with Lemma 7 gives
Concentrating on the projection of the errors, we can exploit the Galerkin orthogonality to obtain
Due to the orthogonality of the \(\mathrm {L}^{2}\)-projections, we have that \(a(\varvec{\xi }_{\varvec{u}},\varvec{\eta }_{\varvec{u}})=0\) and \(d(\varvec{\xi }_{\varvec{\omega }},\varvec{\eta }_{\varvec{\omega }})=0\). Then, from the definition of the form \({\mathcal {A}}\) it follows that
Note that all these terms can be controlled using Lemma 8. Indeed we have
where \(C_{1}=C_{b_{1}}+(C_{b_{1}}+C_{b_{2}}+C_{j})+(C_{b_{2}}+C_{e})+\frac{C_{\varvec{\omega }}}{\sigma }\). Next we only need to estimate \(c(\varvec{\eta }_{\varvec{\omega }},\varvec{\eta }_{\varvec{u}})\) in (4.20). To do so we use Young’s inequality
Substituting (4.21) and (4.22) back into (4.20), we obtain the bound
and thanks to assumption (4.13), we can arrive at
where \(C'=C_{1}(\min \{\frac{1}{2},1-\frac{2\Vert \varvec{\beta }\Vert _{\infty }^2}{\nu \sigma }\})^{-1}\). The error estimate in (4.17) is then obtained by combining estimates (4.19), (4.23) and using triangle inequality.
We now turn to the estimate of the \(\mathrm {L}^{2}\)-norm of error in the pressure (4.18). Since \(e_{p}\in \mathrm {L}^{2}_{0}(\varOmega )\), we can find \(\varvec{z}\in \mathrm {H}^{1}_{0}(\varOmega )^3\) such that (see for instance [28, Chapter I, Corollary 2.4])
where \(\kappa >0\) is the inf-sup constant. Therefore, we infer from (4.12) that
where we have used that \(j(e_{\varvec{u}},\varvec{z})=0\) for \(\varvec{z}\in \mathrm {H}^{1}_{0}(\varOmega )^3\). Using the Galerkin orthogonality, we obtain
Therefore
Next, by the definition of the form \({\mathcal {A}}\), we have the relation
and then we can write
The first term in the right-hand side of the above inequality can be easily estimated by using Lemma 8. Indeed, it follows by choosing \(\varvec{v}_h=\varvec{\xi }_{\varvec{z}}\) that
Let us now estimate each of the terms \(T_i\), \(i=1,\ldots ,7\) in (4.25). Using the properties of the \(\mathrm {L}^2-\)projection, Cauchy–Schwarz’s inequality and the inverse inequality, we obtain the bounds
Then, using (4.23) and (4.24), we can deduce that
Furthermore, since \(\int _{T} \varvec{\xi }_{\varvec{z}}\cdot \nabla \eta _p\, d\varvec{x}=0\), we get from (4.11) the following estimates
Then, using (4.23) and (4.24), we can infer that
Next, using Cauchy–Schwarz’s inequality and again (4.23) together with (4.24). we get
where \(\mu _{h}=\frac{h\Vert \varvec{\beta }\Vert _{\infty }}{\sqrt{\nu }}\). Similarly, we have:
The terms \(T_5\), \(T_6\) and \(T_7\) can be readily estimated, much in the same way as before, using the error bound in (4.17) and the fact that \(\varvec{z}\in \mathrm {H}^{1}_{0}(\varOmega )^3\).
Finally, the pressure estimate follows after putting all individual bounds back into (4.25). \(\square \)
Note that, differently from the conforming method introduced in Sect. 3, the discrete velocity generated by scheme (4.9) is not necessarily divergence-free, and as a consequence the method is not pressure-robust. A recent remedy in the context of DG methods can be found in [31].
5 Numerical Tests
We present a set of examples to confirm numerically the convergence rates anticipated in Theorems 4 and 5. We stress that whenever \(\varGamma =\partial \varOmega \), the zero-mean condition enforcing the uniqueness of the Bernoulli pressure is implemented using a real Lagrange multiplier (which amounts to adding one row and one column to the corresponding matrix system). Linear solves are performed with the direct method SuperLU.
5.1 Test 1: Experimental Convergence in 2D
For our first example we produce the error history associated with the proposed mixed finite element and mixed DG approximations. Let us consider the following closed-form solutions to the Oseen equations defined on the unit square domain \(\varOmega =(0,1)^2\):
The exact velocity has zero normal component on the whole boundary, and the exact vorticity is employed to impose a non-homogeneous vorticity trace. In this example we are assuming that \(\varGamma =\partial \varOmega \), and the exact Bernoulli pressure fulfils the null-average condition. We consider the model parameters \(\nu = 0.1\) and \(\sigma = 10\), and the convecting velocity \(\varvec{\beta }\) is taken as the exact velocity solution, which in particular satisfies the bound (2.15).
Experimental convergence in 2D. Lowest-order mixed finite element approximation of velocity magnitude (a), vorticity (b), and Bernoulli pressure (c) on the unit square
Test 1A On a sequence of uniformly refined meshes we compute errors between the exact and approximate solutions, measured in the norms involved in the convergence analysis of Sect. 3. The obtained error history is reported in Table 1, where the rightmost column displays the \(\ell ^\infty -\)norm of the nodal values of the velocity divergence projected to the space \(\mathrm {Q}_h\), all approaching machine precision. The asymptotic \(O(h^{k+1})\) decay of the error observed for each field variable confirms the overall optimal convergence predicted by Theorem 4. Sample approximate solutions generated with the lowest order method on a coarse mesh are portrayed in Fig. 1.
Test 1B With the purpose of verifying that the mixed finite element scheme is pressure-robust and also performs well for different viscosity values, we modify the exact pressure to be \(p(x,y) = 1000(x^4-y^4)\) and rerun Test 1A for decreasing values of \(\nu \) (down to \(\nu =\)1e-6). Sample values of the error history are collected in Table 2, where we only show the results for the lowest-order method. Even if the pressure error shows optimal convergence rates, the error values are quite large but these do not affect the error decay of the velocity, suggesting pressure-robustness. We also see that smaller viscosity values do not pollute the approximations. In view of (3.9) and Remark 3, perhaps a clearer numerical confirmation of pressure-robustness can be observed by choosing exact velocity and exact vorticity belonging to the finite element subspaces (for instance, we can use constant values \(\varvec{u}= (0,0)\), \(\varvec{\omega }= 0\)), while keeping the exact pressure \(p(x,y) = x^4-y^4\). As discussed in [29, Sect. 3], pressure robustness of the method will imply that the velocity and vorticity errors will be practically zero whereas the pressure error will decay with \(O(h^{k+1})\). Table 3 shows that this is precisely the case, where we have taken \(\nu = 0.01\).
Test 1C An analogous example is now carried out to confirm numerically the convergence rates of the DG methods defined by (4.9). The same model parameters and closed-form solutions as in Test 1A are used, and the stabilisation constants in (4.6)–(4.8) take the values \(a_{11}=c_{11}=\sigma \) and \(d_{11}=\nu \). The results collected in Table 4 indicate that the DG scheme converges optimally when we measure errors in the energy \({\mathcal {A}}\)-seminorm (4.14) and in the \(\mathrm {L}^2-\)norm of the pressure. Here, however, we do not expect divergence-free approximate velocities. Moreover, for the DG case we do not have pressure-robustness. Nevertheless, as evidenced in Table 5, the errors in the \({\mathcal {A}}\)-norm and the pressure errors exhibit optimal decay regardless of the viscosity values.
5.2 Test 2: Transient Flow in an Open Cavity
In this example we illustrate a more complex problem involving the non-stationary behaviour of the flow in an open 2D cavity. The main compartment of the domain consists on a rectangle \((0,1.2)\times (0,1)\) whereas smaller rectangles \((0.25,0.45)\times (-0.1,0)\) and \((1.2,1.3)\times (0.7,0.9)\) play the role of inlet and outlet channels. The domain is discretised into an unstructured mesh of 35433 triangular elements. Normal velocities and a compatible trace vorticity are imposed on the whole boundary \(\varGamma =\partial \varOmega \) according to
that is, parabolic inlet and outlet profiles together with slip velocities elsewhere on \(\partial \varOmega \). We set a fluid viscosity of \(\nu = 0.001\) and use as initial velocity the solution adapted from the previous test \(\varvec{u}_0(x,y) =[\sin (\pi /1.3 x)^2\sin (\pi /1.1 (y+0.1))^2\cos (\pi /1.1 (y+0.1)), -\frac{1}{3}\sin (2/1.3\pi x)\sin (\pi /1.1 (y+0.1))^3]^T\). The parameter \(\sigma = 10\) indicates a timestep of \(\varDelta t = 0.1\), and a backward Euler discretisation implies that we take \(\varvec{f}= \sigma {\hat{\varvec{u}}}\), where \({\hat{\varvec{u}}}\) denotes the velocity approximation at the previous iteration. The convective velocity \(\varvec{\beta }= {\hat{\varvec{u}}}\) therefore needs to be updated at each iteration. At least for the initial solution we have that the convecting velocity satisfies the assumption (2.15). The simulation is run until \(T_{\text {final}} = 4\varDelta t\) and we present in Fig. 2 two snapshots of the numerical solutions at \(t = \varDelta t\) and \(t = T_{\text {final}}\), computed with a second-order DG scheme, and using the stabilisation parameters \(a_{11}=c_{11}=\sigma \) and \(d_{11}=\nu \). From the velocity plots (including a line integral convolution visualisation), we can evidence the formation of a main vortex on the centre of the domain plus smaller recirculation areas that emerge on the top left and bottom left corners, together with a preferential path joining the inlet and outlet boundaries.
Test 2: Second-order DG approximation of the transient flow patterns in an open cavity after one timestep (a–c) and after four time steps (d–f), using \(\varDelta t = 0.1\) and \(\nu = 0.001\)
5.3 Test 3: Lid Driven Cavity Flow
For this classical benchmark problem we consider zero external forces and concentrate on the case where flow recirculation occurs by Dirichlet conditions only. First we consider the two-dimensional case, where on the top lid of the unit square (at \(y=1\)) we set a unidirectional velocity of unit magnitude, whereas no-slip velocity and zero tangential vorticity are imposed on the remaining sides of the boundary. We set the parameters \(\nu =0.001,\sigma =50\) and employ a structured mesh of 4096 elements. The initial velocity is computed from a Stokes solution (setting both \(\sigma \) and \(\varvec{\beta }\) to zero). We compare the results obtained with our lowest-order FE scheme against the benchmark data from [13] (produced with a spectral method applied to a vorticity-based formulation). Figure 3a shows the generated velocity profile at \(t=20\), having all the flow features expected for this regime. The solid lines in Fig. 3b–e portray cuts of the solution on the mid-lines of the domain, whereas the circle markers indicate benchmark data. The approximate vorticity has been rescaled with \(\nu ^{-1/2}\) to reflect the overall agreement with the results reported in [13].
Test 3. Line integral convolution visualisation of velocity for the 2D cavity flow benchmark with \(\nu =0.001\) (a), cuts of the vertical velocity and vorticity on the line \(y=0.5\) (b, c); and horizontal velocity and vorticity on the line \(x=0.5\) (d, e). The circle markers indicate benchmark values from [13]. Lowest-order approximation of velocity (f) for the 3D case with \(\nu =0.0025\), shown at \(t=20\); and profile of the horizontal velocity on the line \(y=0.5,z=0.5\) (g), and of the vertical velocity on the line \(x=0.5,y=0.5\) (h). The asterisks indicate benchmark values from [23] and all numerical approximations for this test were obtained with the lowest-order mixed finite element method
We also test the 3D implementation and formulation by conducting the same benchmark on the unit cube \(\varOmega = (0,1)^3\). Again, boundary \(\varSigma \) is the top plate (defined by \(z=1\)), where we set tangential velocity of magnitude one, and on \(\varGamma =\partial \varOmega \setminus \varSigma \) we consider no-slip velocities and zero tangential vorticity. The fluid viscosity is now \(\nu = 0.0025\) and a structured mesh of 58752 tetrahedral elements is employed. Once again we focus on a non-stationary regime with a backward Euler scheme, now using \(\sigma = 10\), and proceed to update the convective velocity and the right-hand side using the velocity approximation at the previous time iteration. The velocity field for a converged solution after 200 time steps is shown in Fig. 3f. We can observe the expected asymmetric vortex forming parallel to the xz plane (also the generation of high pressure near the corners where the Dirichlet velocity datum has a discontinuity). We have also compared our results with the benchmark values obtained in [23] (using multiquadric differential quadratures) for a Reynolds number of 400: the solid lines in Fig. 3g–h show velocity profiles captured on the plane \(y=0.5\), concentrating on the vertical and horizontal centrelines, where we also include the data from [23] (in asterisks) showing a reasonable match (we have rotated the data, as in their tests the unit velocity is imposed on the face \(y=1\)).
Test 4. Transients at early (left), moderate (middle column) and advanced (right) times, for the Kelvin–Helmholtz instability problem computed with the first-order DG scheme (4.9). Velocity magnitude (a–c), scalar vorticity (d–f), and Bernoulli pressure (g–i)
5.4 Test 4: Kelvin–Helmholtz Mixing Layer
We close this section with a benchmark test related to the well-known vortex formation mechanisms known as the Kelvin–Helmholtz instability problem. The setup of the test follows the specifications in [35] (see also [14]), the transient Oseen equations are solved on the unit square \(\varOmega =(0,1)^2\) and the bottom and top walls constitute \(\varGamma \), where we impose a free-slip velocity condition and zero vorticity. The left and right walls are regarded as a periodic boundary. The initial velocity is
with perturbation scaling \(c_n = 0.001\), reference velocity \(u_\infty = 1\), \(w_{\infty }^a = 8\pi \), \(w_{\infty }^b = 20\pi \), \(\delta _0 = 1/28\). The characteristic time is \(\bar{t}=\delta _0/u_{\infty }\), the Reynolds number is Re\(=10000\), and the kinematic viscosity is \(\nu = \delta _0u_{\infty }/\text {Re}\). We use a structured mesh of 128 segments per side, representing 131072 triangular elements, and we solve the problem using our first-order DG scheme, setting again the stabilisation constants to \(a_{11}=c_{11}=\sigma = 1/\varDelta t\) and \(d_{11}=\nu \), where the timestep is taken as \(\varDelta t = \bar{t}/20\). The specification of this problem implies that the solutions will be quite sensitive to the initial perturbations present in the velocity, which will amplify and consequently vortices will appear. We proceed to compute numerical solutions until the dimensionless time \(t = 7\), and present in Fig. 4 sample solutions at three different simulation times. For visualisation purposes we zoom into the region \(0.25 \le y \le 0.75\), where all flow patterns are concentrated.
References
Alonso, A., Valli, A.: An optimal domain decomposition preconditioner for low-frequency time harmonic Maxwell equations. Math. Comput. 68, 607–631 (1999)
Amara, M., Capatina-Papaghiuc, D., Trujillo, D.: Stabilized finite element method for Navier–Stokes equations with physical boundary conditions. Math. Comput. 76(259), 1195–1217 (2007)
Amoura, K., Azaïez, M., Bernardi, C., Chorfi, N., Saadi, S.: Spectral element discretization of the vorticity, velocity and pressure formulation of the Navier–Stokes problem. Calcolo 44(3), 165–188 (2007)
Anaya, V., Bendahmane, M., Mora, D., Ruiz-Baier, R.: On a vorticity-based formulation for reaction-diffusion-Brinkman systems. Netw. Heterog. Media 13(1), 69–94 (2018)
Anaya, V., Gatica, G.N., Mora, D., Ruiz-Baier, R.: An augmented velocity–vorticity–pressure formulation for the Brinkman problem. Int. J. Numer. Methods Fluids 79(3), 109–137 (2015)
Anaya, V., Mora, D., Oyarzúa, R., Ruiz-Baier, R.: A priori and a posteriori error analysis for a mixed scheme for the Brinkman problem. Numer. Math. 133(4), 781–817 (2016)
Anaya, V., Mora, D., Ruiz-Baier, R.: An augmented mixed finite element method for the vorticity–velocity–pressure formulation of the Stokes equations. Comput. Methods Appl. Mech. Eng. 267, 261–274 (2013)
Anaya, V., Mora, D., Ruiz-Baier, R.: Pure vorticity formulation and Galerkin discretization for the Brinkman equations. IMA J. Numer. Anal. 37(4), 2020–2041 (2017)
Azaïez, M., Bernardi, C., Chorfi, N.: Spectral discretization of the vorticity, velocity and pressure formulation of the Navier–Stokes equations. Numer. Math. 104(1), 1–26 (2006)
Bochev, P.B.: Analysis of least-squares finite element methods for the Navier–Stokes equations. SIAM J. Numer. Anal. 34(5), 1817–1844 (1997)
Bochev, P.B.: Negative norm least-squares methods for the velocity–vorticity–pressure Navier–Stokes equations. Numer. Methods PDEs 15(2), 237–256 (1999)
Bonelle, J., Ern, A.: Analysis of compatible discrete operator schemes for the Stokes equations on polyhedral meshes. IMA J. Numer. Anal. 35(4), 1672–1697 (2015)
Botella, O., Peyret, R.: Benchmark spectral results on the lid-driven cavity flow. Comput. Fluids 27(4), 421–433 (1998)
Burman, E., Ern, A., Fernández, M.A.: Fractional-step methods and finite elements with symmetric stabilization for the transient Oseen problem. ESAIM Math. Model. Numer. Anal. 51(2), 487–507 (2017)
Cai, Z., Chen, B.: Least-squares method for the Oseen equation. Numer. Methods PDEs 32, 1289–1303 (2016)
Caucao, S., Mora, D., Oyarzúa, R.: A priori and a posteriori error analysis of a pseudostress-based mixed formulation of the Stokes problem with varying density. IMA J. Numer. Anal. 36(2), 947–983 (2016)
Chang, C.L., Yang, S.-Y.: Analysis of the \([L^2, L^2, L^2]\) least-squares finite element method for incompressible Oseen-type problems. Int. J. Numer. Anal. Model. 4(3–4), 402–424 (2007)
Çeşmelioğlu, A., Cockburn, B., Nguyen, N.C., Peraire, J.: Analysis of HDG methods for Oseen equations. J. Sci. Comput. 55(2), 392–431 (2013)
Ciarlet, P.: The Finite Element Method for Elliptic Problems. North-Holland, Amsterdam (1978)
Cockburn, B., Kanschat, G., Schötzau, D.: The local discontinuous Galerkin method for linear incompressible fluid flow: a review. Comput. Fluids 34(4–5), 491–506 (2005)
Cockburn, B., Kanschat, G., Schötzau, D.: The local discontinuous Galerkin method for the Oseen equations. Math. Comput. 73(246), 569–593 (2003)
Cockburn, B., Kanschat, G., Schötzau, D., Schwab, C.: Local discontinuous Galerkin methods for the Stokes system. SIAM J. Numer. Anal. 40(1), 319–343 (2002)
Ding, H., Shu, C., Yeo, K.S., Xu, D.: Numerical computation of three-dimensional incompressible viscous flows in the primitive variable form by local multiquadric differential quadrature method. Comput. Methods Appl. Mech. Eng. 195, 516–533 (2006)
Duan, H.-Y., Liang, G.-P.: On the velocity–pressure–vorticity least-squares mixed finite element method for the 3D Stokes equations. SIAM J. Numer. Anal. 41(6), 2114–2130 (2003)
Dubois, F., Salaün, M., Salmon, S.: First vorticity–velocity–pressure numerical scheme for the Stokes problem. Comput. Methods Appl. Mech. Eng. 192(44–46), 4877–4907 (2003)
Gatica, G.N.: A Simple Introduction to the Mixed Finite Element Method. Theory and Applications. Springer Briefs in Mathematics. Springer, Cham (2014)
Gatica, G.N., Gatica, L.F., Márquez, A.: Augmented mixed finite element methods for a vorticity-based velocity–pressure–stress formulation of the Stokes problem in 2D. Int. J. Numer. Methods Fluids 67(4), 450–477 (2011)
Girault, V., Raviart, P.A.: Finite Element Methods for Navier–Stokes Equations. Theory and Algorithms. Springer, Berlin (1986)
John, V., Linke, A., Merdon, C., Neilan, M., Rebholz, L.G.: On the divergence constraint in mixed finite element methods for incompressible flows. SIAM Rev. 59(3), 492–544 (2017)
Kreeft, J., Gerritsma, M.: Mixed mimetic spectral element method for Stokes flow: a pointwise divergence-free solution. J. Comput. Phys. 240, 284–309 (2013)
Lederer, P.L., Lehrenfeld, C., Schöberl, J.: Hybrid discontinuous Galerkin methods with relaxed H(div)-conformity for incompressible flows. Part I. SIAM J. Numer. Anal. 56(4), 2070–2094 (2018)
Mohapatra, S., Ganesan, S.: A non-conforming least squares spectral element formulation for Oseen equations with applications to Navier–Stokes equations. Numer. Funct. Anal. Optim. 37(10), 1295–1311 (2016)
Nicolaides, R.A.: Existence, uniqueness and approximation for generalized saddle point problems. SIAM J. Numer. Anal. 19, 349–357 (1982)
Salaün, M., Salmon, S.: Low-order finite element method for the well-posed bidimensional Stokes problem. IMA J. Numer. Anal. 35(1), 427–453 (2015)
Schroeder, P.W., John, V., Lederer, P.L., Lehrenfeld, C., Lube, G., Schöberl, J.: On reference solutions and the sensitivity of the 2D Kelvin–Helmholtz instability problem. Comput. Math. Appl. 77(4), 1010–1028 (2019)
Tsai, C.-C., Yang, S.-Y.: On the velocity–vorticity–pressure least-squares finite element method for the stationary incompressible Oseen problem. J. Comput. Appl. Math. 182(1), 211–232 (2005)
Acknowledgements
This work has been supported by CNRS through the PEPS programme; by CONICYT - Chile through FONDECYT project 11160706; by DIUBB, Universidad del Bío-Bío through projects 165608–3/R and 171508 GI/VC; by CONICYT-Chile through the project AFB170001 of the PIA Program: Concurso Apoyo a Centros Científicos y Tecnológicos de Excelencia con Financiamiento Basal. by DIDULS, Universidad de La Serena through project PR17151; and by Unversidad de Córdoba through the programme “Estrategias para la Sostenibilidad de Grupos de Investigación”, project FCB-03-17.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Anaya, V., Bouharguane, A., Mora, D. et al. Analysis and Approximation of a Vorticity–Velocity–Pressure Formulation for the Oseen Equations. J Sci Comput 80, 1577–1606 (2019). https://doi.org/10.1007/s10915-019-00990-7
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10915-019-00990-7
Keywords
- Oseen equations
- Vorticity-based formulation
- Mixed finite elements
- Exactly divergence-free velocity
- Discontinuous Galerkin schemes
- Numerical fluxes
- A priori error bounds