Analysis and approximation of a vorticity-velocity-pressure formulation for the Oseen equations

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. 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.


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,8] as well as stabilised and least-squares schemes [2,9] for Navier-Stokes; also several mixed and augmented methods for Brinkman [4,5,7], and a number of other discretisations specifically designed for Stokes flows [6,22,23,25,30].
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 velocitypressure formulation) include for instance [10], that presents a least-squares method for Navier-Stokes equations with vorticity-based first-order reformulation, and whose analysis exploits the elliptic theory of Agmon-Douglas-Nirenberg. 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 [27], where a suitable preconditioner is also proposed. In [32] 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 div least-squares minimisation problem based on the stress-velocity-pressure formulation was introduced in [13]. 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 [15], also including vorticity and total pressure in the formulation.
Discontinuous Galerkin (DG) methods have also been used to solve the Oseen problem, as for example, in [18,19] 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 [16] for the Oseen equations. The family of DG methods we propose here has resemblance with those schemes, but exploits 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 [24] 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 ≥ 0 for the velocity field, Nédélec elements or 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. The method produces exactly divergencefree 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 [19,20], we present another scheme involving the discontinuous Galerkin discretisation of the curl-curl and grad-div operators. We prove the well-posedness of the DG scheme and derive error estimates under some solution regularity assumptions.
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 Section 2. The finite element discretisation is constructed in Section 3, where we also derive the stability and convergence bounds. In Section 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 Section 5 with a set of numerical tests that illustrate the properties of the proposed numerical schemes in a variety of scenarios.

Statement and solvability of the continuous problem
Oseen problem in terms of velocity-vorticity-pressure. A standard backward Euler timediscretisation 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 (see [26,29]): where ν > 0 is the kinematic fluid viscosity, σ > 0 is inversely proportional to the time-step, β is an adequate approximation of velocity to be made precise below, and the vector of external forces 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 + 1 2 |u| 2 , where P is the actual fluid pressure.
The structure of (2.1) suggests to introduce the rescaled vorticity vector ω := √ ν curl 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 Γ ⊂ ∂Ω, whereas a non-homogeneous tangential velocity u Σ and a fixed Bernoulli pressure p Σ are set on the remainder of the boundary Σ = ∂Ω \ Γ. Therefore, system (2.1) can be recast in the form where n stands for the outward unit normal on ∂Ω. Should the boundary Σ have zero measure, the additional condition (p, 1) Ω,0 = 0 is required to enforce uniqueness of the Bernoulli pressure. where the operator γ t is the tangential trace operator on Γ, defined by: γ t (θ) = θ ×n. Let us endow H and Q with their natural norms. For the space Z however, we consider the following viscosity-weighted norm: From now on, we will assume that the data are regular enough: β ∈ L ∞ (Ω) 3 and f ∈ L 2 (Ω) 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: where the bilinear forms a : and the linear functionals F : H → R, and G : Z → R are specified as follows for all u, v ∈ H, ω, θ ∈ Z, and q ∈ Q.
Solvability analysis. In order to analyse the variational formulation (2.3), let us introduce the Kernel of the bilinear form b 2 (·, ·) and its classical characterisation and let us recall that b 2 satisfies the inf-sup condition: with an inf-sup constant β 2 > 0 only depending on Ω (see e.g. [24]).
We will now address the well-posedness of (2.3). To that end, it is enough to study its reduced counterpart, defined on X × Z: Find (u, ω) ∈ X × Z such that The equivalence between (2.3) and (2.5) is established in the following result, whose proof follows [26] and it is basically a direct consequence of the inf-sup condition (2.4).
The abstract setting that will permit the analysis of (2.3) is stated in the following general result [24, Theorem 1.2].
Theorem 2.1 Let A : X × X → R be a bounded bilinear form and G : X → R a bounded functional, both defined on the Hilbert space (X , ·, · X ). If there exists α > 0 such that

6)
and sup then there exists a unique solution x ∈ X to the problem Furthermore, there exists C > 0 (independent of x) such that Lemma 2.2 Let us assume that 2 β 2 ∞,Ω νσ < 1, (2.8) and let us define the bilinear form A(·, ·) specified as Then, there exist α 1 , α 2 > 0 such that where X := X × Z, endowed with the corresponding product norm, is a Hilbert space.
Proof. For all (v, θ) ∈ X , we have that: As a consequence of the previous lemmas, we have the following result.
Theorem 2.2 Let us assume (2.8). Then, the variational problem (2.5) admits a unique solution (u, ω) ∈ X × Z. Moreover, there exists C > 0 such that Proof. It suffices to verify the hypotheses of Theorem 2.1. First, we define the linear functional which is bounded on X × Z. Thus, the proof follows from Lemmas 2.2 and 2.3.
The following result establishes the corresponding stability estimate for the Bernoulli pressure.
Corollary 1 Let (u, ω) ∈ X × Z, be the unique solution of (2.5), with u and ω satisfying (2.11). In addition, let p ∈ Q be the unique pressure provided by Lemma 2.1, so that (u, ω, p) ∈ H × Z × Q is the unique solution of (2.3). Then, there exists C > 0 such that

Proof. Combining the inf-sup condition (2.4) with the first equation in (2.3) gives the bound
which together with (2.11), and the boundedness of F , a, b 1 and c, complete the proof.
Remark 2.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 [14]. Establishing inf-sup conditions for the off-diagonal bilinear forms in the original nonsymmetric formulation is, however, much more involved (see e.g. [28]).
Remark 2.2 Assumption (2.8) holds provided one chooses σ appropriately. As this parameter represents the inverse of the timestep, the aforementioned relation constitutes then a CFL-type condition. We also note that this bound for β coincides with the hypotheses that yield solvability of least-squares formulations for the Oseen problem analysed in [13,15].

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.
Defining the discrete problem. Let {T h (Ω)} h>0 be a shape-regular family of partitions of the polyhedral regionΩ, by tetrahedrons T of diameter h T , with mesh size h := max{h T : T ∈ T h (Ω)}. In what follows, given an integer k ≥ 0 and a subset S of R 3 , P k (S) will denote the space of polynomial functions defined locally in S and being of total degree ≤ k.
Moreover, for any T ∈ T h (Ω), we introduce the local Nédélec space where R k+1 (T ) is a subspace of P k+1 (T ) 3 composed by homogeneous polynomials of degree k + 1, and being orthogonal to x. With these tools, let us define the following finite element subspaces: The proposed Galerkin scheme approximating (2.3) reads as follows: 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 div H h ⊆ Q h . Moreover, it is well-known that the following discrete inf-sup condition holds (see [24]): We again resort to a reduced version of the problem, now defined on the product space and its equivalence with (3.4) is once more a direct consequence of the inf-sup condition (3.6).
In order to establish the well-posedness of (3.7), we will employ the following discrete version of Theorem 2.1.
Theorem 3.1 Assume (2.8). Let k ≥ 0 be an integer and let X h and Z h be given by (3.5) and (3.1), respectively. Then, there exists a unique (u h , ω h ) ∈ X h × Z h solution of the discrete scheme (3.7). Moreover, there exist positive constantsĈ 1 ,Ĉ 2 > 0 independent of h such that and where (u, ω) ∈ X × Z is the unique solution to (2.5).
Proof. Let us define X h := X h × Z h and reuse the forms A(·, ·) and G(·) as in the proof of Lemma 2.2. The next step consists in proving that A(·, ·) satisfies the discrete version of the inf-sup conditions (2.6)-(2.7), as in Lemmas 2.2 and 2.3. In order to assert (2.6), we consider (u h , ω h ) ∈ X h , and definẽ Then, repeating exactly the same steps used in the proof of Lemma 2.2 the discrete version of (2.6) follows. Regarding the discrete version of (2.7), we once again repeat the same arguments given in the proof of Lemma 2.3. Finally, the Céa estimate follows from classical arguments.
We now state the stability and an adequate approximation property of the discrete pressure.
In addition, let p h ∈ Q h be the unique discrete Bernoulli pressure provided by Lemma 3.1, so that Then, there exist positive constants C 1 ,C 2 > 0, independent of h and ν, such that 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.
A priori error estimates. Let us introduce for a given s > 1/2, the Nédeléc global interpolation operator R h : H s (curl; Ω) ∩ Z → Z h . From [1] we known that for all θ ∈ H s (curl; Ω) with s > 1/2, there exists C > 0 independent of h, such that On the other hand, for the Raviart-Thomas interpolation Π h : H s (Ω) 3 ∩ H → H h , with s > 0, we recall (see e.g. [24]) that there exists C > 0, independent of h, such that for all s > 0: Finally we recall that the orthogonal projection from L 2 (Ω) onto the finite element subspace Q h , here denoted P h , satisfies the following error estimate for all s > 0: These operators fulfil the following commuting diagram The following result summarises the error analysis for our mixed finite element scheme (3.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.
Preliminaries. Apart from the definitions laid out at the beginning of Section 3, let us denote by E h the set of internal faces, by F Σ h the set of external faces on Σ and by F Γ h the set of external faces on Γ. We set We denote by h e the diameter of each face e. Let T + and T − be two adjacent elements of T h and let n + (respectively n − ) be the outward unit normal vector on ∂T + (respectively ∂T − ). For a vector field u, we denote by u ± the trace of u from the interior of T ± . We define jumps and averages and adopt the convention that for boundary faces 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 Section 3, that isQ h := Q h .
Discrete formulation and solvability analysis. Multiplying each equation in (2.2) by suitable functions, the resulting DG scheme consists in finding where u w h , u p h , ω h and p h are numerical fluxes, which approximate the traces of u h , ω h and p h on the boundary. The fluxes ω h and u w h are related to the curl-curl operator and are defined by whereas the fluxes u p h and p h are associated with the grad-div operator and defined by (4.5) The parameters C 11 , A 11 and D 11 are positive stabilisation parameters, and following [20] 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 C 11 and C 11 (respectively D 11 , D 11 and A 11 , A 11 ).
We then proceed to integrate by parts equations (4.1) and (4.3), and then summing up over all T ∈ T h , we obtain the following DG scheme: where the forms a, c and d are the same in (2.3), whileb 1 ,b 2 , j and e are defined, respectively, by: In addition, the linear functionals F , G and L associated with the source terms are defined as: By integration by parts and as a consequence of the identity: it follows that the formb 1 can be written as: Similarly, using the following identity: the formb 2 can be recast, after integration by parts, as follows Note that, differently from the conforming method introduced in Section 3, the discrete velocity generated by scheme (4.9) is not necessarily divergence-free.
To simplify the exposition of the analysis of the method, we will write the mixed scheme (4.9) in the following equivalent form: where Let us now show the existence and uniqueness of solution to formulation (4.9).

Proposition 1
The DG method (4.9) with the numerical fluxes given by (4.4)-(4.5) defines a unique approximate solution Proof. Since the problem is linear and finite dimensional, it suffices to show that if f = 0, p Σ = 0 and u Σ = 0, then (u h , ω h , p h ) = (0, 0, 0). To this end, take v h = u h , θ h = ω h and q h = p h in (4.9), summing up the three equations, we obtain: It follows that Using Young's inequality, we can assert that Therefore, in particular we have that: 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 L 2 -norm. For this, we introduce the following semi-norm | · | A : |(u, ω, p)| 2 A := σ u 2 0,Ω + ω 2 0,Ω + |u| 2 j + |p| 2 e . (4.14) We shall suppose that the exact solution (u, p) satisfies the following regularity u ∈ H 1+s (Ω) 3 , and p ∈ H s (Ω), s ≥ 1. And we will also employ the following norm We define e u = u − u h , e ω = ω − ω h and e p = p − p h . Let us denote by Π H (respectively Π Z and Π Q ) the L 2 -projection onto H (respectively Z and Q), and let us split the errors in the following manner e u = ξ u + η u , e ω = ξ ω + η ω and e p = ξ p + η p , where the numerical and approximation errors are defined by: We recall the following standard approximation properties (see for instance [17]). And as a consequence, we have the following result.
Lemma 4.2 Suppose that the analytical solution (u, p) of (2.1) satisfies (4.15) and we set ω = √ ν curl u. 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.16). Next we can state that Recalling that C 11 , A 11 and D 11 are as in (4.6)-(4.8) and using (4.17), we end up with the bound and similarly we obtain Proof. The bounds associated with the formb 2 can be proved exactly with the same arguments as [20,Section 3.3]. Let us now deal with the termb 1 . We observe that due to the properties of the L 2 -projection, we can write T ξ u · curl θ h dx = 0.
Using Cauchy-Schwarz's inequality, we then readily obtain , and the desired estimate follows from the inverse inequality and Lemma 4.2.
Similarly, using (4.10) and again the properties of the L 2 -projections we have and then, we simply have to use (4.6) and Lemma 4.1.
The estimates for the forms j and e are obtained in a similar way. Indeed, using once again Cauchy-Schwarz's inequality and Lemma 4.2, it follows that and proceeding analogously as before, we get Finally, concerning the term c(ξ ω , v h ) we can assert that and exploiting the previous bounds, we obtain the corresponding estimate with where C A and C are positive constants independent of the meshsize.
Proof. We begin with the estimate (4.18). A direct application of the definition of the A−seminorm in combination with Lemma 4.2 gives Concentrating on the projection of the errors, we can exploit the Galerkin orthogonality to obtain Due to the orthogonality of the L 2 -projections, we have that a(ξ u , η u ) = 0 and d(ξ ω , η ω ) = 0. Then, from the the definition of the form A it follows that Note that all these terms can be controlled using Lemma 4.3. Indeed we have where C 1 = C b1 + (C b1 + C b2 + C j ) + (C b2 + C e ) + Cω σ . Next we only need to estimate c(η ω , η u ) in (4.21). To do so we use Young's inequality Substituting (4.22) and (4.23) back into (4.21), we obtain the bound and thanks to assumption (4.13), we can arrive at The error estimate in (4.18) is then obtained by combining estimates (4.20), (4.24) and using triangle inequality.
Let us now estimate each of the terms T i , i = 1, . . . , 7 in (4.26). Using the properties of the L 2 −projection, Cauchy-Schwarz's inequality and the inverse inequality, we obtain the bounds Then, using (4.24) and (4.25), we can deduce that Furthermore, since T ξ z · ∇η p dx = 0, we get from (4.11) the following estimates Then, using (4.24) and (4.25), we can infer that T 2 ≤ Ch min{s, k+1} (u, p) s e p 0,Ω .
Next, using Cauchy-Schwarz's inequality together with Lemma 4.2 and (4.24), we get 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.18) and the fact that z ∈ H 1 0 (Ω) 3 . Finally, the pressure estimate follows after putting all individual bounds back into (4.26).

Numerical tests
We present a set of examples to confirm numerically the convergence rates anticipated in Theorem 3.2 and Theorem 4.1 . We stress that whenever Γ = ∂Ω, the zero-mean condition enforcing the uniqueness of the Bernoulli pressure is implemented using a real Lagrange multiplier (which amounts to add one row and one column to the corresponding matrix system). Linear solves are performed with the direct method SuperLU. u(x, y) = sin(πx) 2 sin(πy) 2 cos(πy) − 1 3 sin(2πx) sin(πy) 3 , ω(x, y) = √ ν curl u, p(x, y) = x 4 − y 4 .
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 Γ = ∂Ω, and the exact Bernoulli pressure fulfils the null-average condition. We consider the model parameters ν = 0.1 and σ = 10, and the convecting velocity β is taken as the exact velocity solution, which in particular satisfies the bound (2.8). 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 Section 3. The obtained error history is reported in Table 1, where the rightmost column displays the ∞ −norm of the nodal values of the velocity divergence projected to the space 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 3.2. Sample approximate solutions generated with the lowest order method on a coarse mesh are portrayed in Figure 1.
An analogous test 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 are used, and the stabilisation constants in (4.6)-(4.8) take the values a 11 = c 11 = σ and d 11 = ν. The results collected in Table 2 indicate that the DG scheme converges optimally when we measure errors in the energy A−seminorm (4.14) and in the L 2 −norm of the pressure. Here however, we do not expect divergencefree approximate velocities.
otherwise, that is, parabolic inlet and outlet profiles together with slip velocities elsewhere on ∂Ω. We set a fluid viscosity of ν = 0.001 and use as initial velocity the solution adapted from the previous test u 0 (x, y) = [sin(π/1.3x) 2 sin(π/1.1(y+0.1)) 2 cos(π/1.1(y+0.1)), − 1 3 sin(2/1.3πx) sin(π/1.1(y+0.1)) 3 ] T . The parameter σ = 10 indicates a timestep of ∆t = 0.1, and a backward Euler discretisation implies that we take f = σû, whereû denotes the velocity approximation at the previous iteration. The convective velocity β =û therefore needs to be updated at each iteration. At least for the initial solution we have that the convecting velocity satisfies the assumption (2.8). The simulation is run until T final = 4∆t and we present in Figure 2 two snapshots of the numerical solutions at t = ∆t and t = T final , computed with a second-order DG scheme, and using the stabilisation parameters a 11 = c 11 = σ and d 11 = ν. 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 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 conditions and zero tangential vorticity are imposed on the remaining sides of the boundary. We set the parameters ν = 0.001, σ = 50 and employ a structured mesh of 4096 elements. The initial velocity is computed from a Stokes solution (setting both σ and β to zero). We compare the results obtained with our lowest-order FE scheme against the benchmark data from [11] (produced with a spectral method applied to a vorticity-based formulation). The approximate vorticity has been rescaled with ν −1/2 to reflect the overall agreement with the results reported in [11].
We also test the 3D implementation and formulation by conducting the same benchmark on the unit cube Ω = (0, 1) 3 . Again, boundary Σ is the top plate (defined by z = 1), where we set tangential velocity of magnitude one, and on Γ = ∂Ω\Σ we consider no-slip velocity and zero tangential vorticity. The fluid viscosity is now ν = 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 σ = 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 Figure 3(f). 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 [21] (using multiquadric differential quadratures) for a Reynolds number of 400: the solid lines in Figures 3(g)-3(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 [21] (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: 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 [31] (see also [12]) with perturbation scaling c n = 0.001, reference velocity u ∞ = 1, w a ∞ = 8π, w b ∞ = 20π, δ 0 = 1/28. The characteristic time ist = δ 0 /u ∞ , the Reynolds number is Re= 10000, and the kinematic viscosity is ν = δ 0 u ∞ /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 = σ = 1/∆t and d 11 = ν, where the timestep is taken as ∆t =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 Figure 4 sample solutions at three different simulation times. For visualisation purposes we zoom into the region 0.25 ≤ y ≤ 0.75, where all flow patterns are concentrated. In addition, a qualitative comparison against benchmark data from [31] is presented in terms of the temporal evolution of the enstrophy E(t) (here we rescale ω h with √ ν to match again the real vorticity). We also record the evolution of the palinstrophy P (t), a quantity that encodes the dissipation process. These quantities are defined,  and we remark that for the palinstrophy we use the discrete gradient associated with the DG discretisation. We show these quantities in Figure 5, where also include results from [31] that correspond to coarse and fine mesh solutions of the Navier-Stokes equations using a high order scheme based on Brezzi-Douglas-Marini elements.