On well-posedness of a velocity-vorticity formulation of the Navier-Stokes equations with no-slip boundary conditions

We study well-posedness of a velocity-vorticity formulation of the Navier--Stokes equations, supplemented with no-slip velocity boundary conditions, a no-penetration vorticity boundary condition, along with a natural vorticity boundary condition depending on a pressure functional. In the stationary case we prove existence and uniqueness of a suitable weak solution to the system under a small data condition. The topic of the paper is driven by recent developments of vorticity based numerical methods for the Navier--Stokes equations.

1. Introduction. The evolution of an incompressible, viscous Newtonian fluid is governed by the Navier-Stokes equations (NSE), which for a given bounded, connected domain Ω ⊂ R 3 with a piecewise, smooth, Lipschitz boundary ∂Ω, an end time T , an initial condition u 0 : Ω → R 3 and force field f : (0, T ] × Ω → R 3 read: find a velocity field u : [0, T ] × Ω → R 3 and a pressure field p : (0, T ] × Ω → R such that ∂u ∂t + (u · ∇)u − ν∆u + ∇p = f in (0, T ] × Ω, where ν > 0 is the kinematic viscosity. The system must also be equipped with appropriate boundary conditions. For suitably smooth solutions, equations (1) can be re-formulated in other equivalent forms, in particular through the introduction of new variables such as vorticity, streamfunction, rate of strain tensor, Bernoulli pressure, etc., see, e.g., [14,18,23]. While different formulations often provide useful insights into both physical and mathematical properties of NSE solutions, it is believed that in general they do not include more information than that already contained in (1). The situation changes, however, if one is interested in numerical solutions. It is well known that by applying a discretization method to different equivalent formulations of the equations, one commonly obtains non-equivalent discrete problems with strikingly different numerical properties and, for particular problems, appropriate discretizations of certain formulations outperform others. For this reason, the literature on computational methods for (1) often considers formulations different from the primitive variables "convective" formulation in (1). Among these, we can mention the conservative form, skew-symmetric form, vorticity-streamfunction formulation, streamfunction formulation, rotational form, and EMAC form. We refer the reader to [6,15,20,28] for a description of these formulations and their discretizations.
For numerical analysis of discretization methods, such as a finite element method, one typically needs discrete counterparts of fundamental a priori bounds and wellposedness results. However due to the reasons outlined above, for discrete solutions such results often cannot be obtained by resorting to the primitive variables that describe the equations in (1). This is the main motivation of this work where, through direct arguments, we show well-posedness of a particular formulation of NSE based on vorticity dynamics. Our analysis overcomes two major (technical) difficulties: (i) The implicit enforcement of the relation w = rot u between two independent variables, vorticity and velocity of the fluid. In fact, our arguments extend to the more general case when this identity is no longer true. (ii) Imposing the practically important no-slip conditions on solid walls on the velocity u leads to natural boundary conditions for the vorticity w, which involve boundary values of the pressure variable.
The two observations above have two immediate consequences: first, the energy estimate for u does not necessarily imply any bound for w; and, second, one has to look for more regular weak solutions to give sense to the vorticity boundary conditions on solid walls. A key new technical result in this paper is the wellposedness and regularity of a Stokes-type problem with non-standard boundary conditions. The problem was introduced in [3,13] but, to the best of our knowledge, it has not been further explored in the literature. Velocity-vorticity formulations have been widely used in the numerical literature and have been found to provide greater accuracy in simulations (compared to more common formulations), especially for flows where boundary effects are critical [29,22,24,11,17,30,19]. Such formulations directly build on the equation of the vorticity dynamics, There are several ways to close the system of equations for w and u. One common way is to supplement (2) with the vector Poisson problem for the velocity, ∆u = rot w. In this paper, we consider the formulation that couples (2) with the rotational form of the momentum equation, We are motivated to follow this path by the recent development of numerical methods for this coupling [26,21,4,27,19]. As these references show, these schemes provide superior long-time numerical stability and more accurate fulfillment of discrete conservation laws (in addition to other well-known benefits of vorticity based numerical formulations).
The system of velocity-vorticity equations should be supplemented with suitable boundary conditions for both u and w. The most practically important are no-slip, no-penetration conditions on solid boundaries of the fluid domain, i.e., u = 0. The corresponding vorticity boundary conditions are not straightforward to obtain and were actually the subject of a long and controversial discussion in the literature; see, e.g., [25] for a brief overview. The normal condition w · n = 0 on ∂Ω is the immediate consequence of u = 0 and w = rot u for u smooth up to the boundary. Here and further n is the unit normal vector on ∂Ω. To define tangential vorticity boundary conditions, we follow [25], where the following physically consistent boundary condition on solid walls was derived, This boundary condition can be obtained by taking the tangential component of (3) on ∂Ω, and recalling that ∆u = rot w and u = 0 on ∂Ω. It can be shown, cf. [25], that for f = 0, condition (4) characterizes the production of the streamwise and spanwise vorticity on no-slip boundaries due to the tangential pressure gradient and vorticity intensification, depending on the shape of ∂Ω. The vorticity boundary condition (4) was used in [25,19] to produce numerical simulations with excellent results. We shall see that the condition in (4) is a natural boundary condition, which enters the weak formulation as a functional. Based on the numerical success of the velocity-vorticity formulation (2)-(3) with boundary conditions u = 0, w · n = 0 on ∂Ω, and (4), a theoretical study of the underlying PDE system is warranted. Hence, the purpose of this work is to initiate the study of the velocity-vorticity system equipped with no-slip velocity boundary conditions, and the corresponding vorticity boundary conditions. We shall focus on a stationary version of the problem: with α ≥ 0. For this system we will show existence and uniqueness of solutions, provided that the problem data is sufficiently small, or α is sufficiently large.
We introduced the non-negative parameter α in the system, with the aim of covering the case when an additional zero order term results from the approximation of time derivative by a finite difference. That is, if a temporal discretization is applied to (1) in rotation form coupled with (2) and boundary conditions, then the system (5) could be considered as a single time step of the resulting semidiscrete system (with right hand sides appropriately modified). One, of course, would be also interested in passing to the limit of time step turning to zero. However, this seems beyond our reach, and so we cannot provide an existence result for the full time dependent problem. This, in particular, is due to the presence of the pressure in the vorticity boundary conditions (4). Nevertheless, the partial analysis that we provide here comes with its own set of unique features. As mentioned above, to overcome them, we need to analyze a class of Stokes-type problems with nonstandard boundary conditions, and Helmholtz-Weyl decompositions of the space L r (Ω), r = 2.
Remark 1 (solenoidal vorticity). We note that (5) explicitly enforces the divergence free condition for vorticity with the help of the Lagrange multiplier η in the vorticity equation. There are several reasons for doing this. First, we find that this helps us in the analysis. Second, this is motivated by numerical experience, since a common finite element method for (2) does not yield divergence free discrete solutions and enforcing it as a constraint benefits the accuracy of the recovered vorticity field. Finally, all the arguments in this paper extend now in a straightforward way to a more generic system of equations, where rot f on the right-hand side of the vorticity equation in (5) is replaced by a sufficiently regular source term f w , not necessary divergence free, which can be interpreted as modeling external sources of vorticity production.
We also note the alternative way of enforcing div w = 0 by re-writing the vorticity equation (2) in the form In this equation, η has the physical meaning of the helical density (rather than an auxiliary variable) and D(w)u = u · ∇w − 1 2 (rot w) × u. When (6) is used as the vorticity equation in (5), the system is called the velocity-vorticity-helicity (VVH) formulation of the NSE [26]. It is unique in the fact that it naturally enforces the mathematical constraint that the vorticity be divergence-free, and also that it solves for helicity directly (via helical density). The analysis of the paper can be extended to the system (5), where the vorticity equation is written in the vorticity-helical density form (6), but we shall not elaborate further on this.
Our presentation is organized as follows. In section 2 we provide some preliminary results along with some new results concerning the Stokes problem with nonstandard boundary conditions. In section 3 we show that system (5) is well posed. The analysis follows carefully the dependence of all estimates on problem parameters, since such dependence is of interest for understanding the properties of numerical solutions. In particular, we are interested in the dependence on α, since it comes from the time-discretization and so is "user-controlled". Conclusions and potential directions for future work are provided in section 4.

2.
Preliminaries. In this work we assume that the fluid occupies a bounded, simply connected domain Ω ⊂ R 3 , with ∂Ω ∈ C 3 . We will follow standard notation concerning function spaces. In particular, we will denote the L 2 (Ω) inner product and norm by (·, ·) and · , respectively. All other norms will be clearly labeled.
The natural function spaces for velocity, vorticity and pressure in the presence of solid walls are (Ω) denotes the subspace of L 2 (Ω) consisting of functions with zero mean. The Nečas inequality implies that there exists β 1 > 0 such that 2.1. The operators rot and div. Let us here recall some well-known results about the operators rot and div, and explore some of their consequences that shall be useful in the sequel. The assumptions on Ω yield that W = H 0 (div) ∩ H(rot) both algebraically and topologically and, as shown in [14, section 3.5], we additionally have that If Ω is convex, then C 1 = 1 in (9).
Since Ω is assumed to be bounded, simply connected, and have a smooth boundary, we have the following result about the solvability of the equation rot v = u, see [14, Theorems 3.5 and 3.6].
Notice that if v ∈ H(rot), then rot v ∈ H(div), and so the condition n · rot v = 0 on ∂Ω is well-defined. We therefore define the space Owing to Lemma 2.1, we can characterize this space.
Proof. The result of the Lemma is given in [13, Remark 2.1], but we prove it here for completeness. For a given v ∈ H(rot), let w = rot v in Lemma 2.1 and note that rot(u−v) = 0 implies that u−v = ∇q for some q ∈ H 1 (Ω). The reverse implication, i.e. u + ∇q ∈ H(rot) for all q ∈ H 1 (Ω), u ∈ H 1 (Ω) satisfying the conditions in (11), is also straightforward since u × n = 0 on ∂Ω implies n · rot u = 0 on ∂Ω.
We also define the subspace of H(rot) consisting of divergence free functions with zero normal component on ∂Ω as Notice that if v ∈ H 0 (rot), we have in particular that (v, ∇q) = 0 ∀ q ∈Ċ(Ω), wherė C(Ω) denotes the set of C ∞ (Ω) functions with compact support. Hence div v = 0 inĊ(Ω) and so, owing to a density argument, div v = 0 in L 2 (Ω). Therefore, v ∈ H 0 (rot) implies v ∈ H(div) and the condition n · v = 0 holds in H − 1 2 (∂Ω), and so we conclude n · v = 0 a.e. on ∂Ω. Therefore, for v ∈ H 0 (rot), we get due to the regularity of ∂Ω that v ∈ H 1 (Ω), and because of the boundary condition we have v ∈ W. Summarizing, we have shown that Owing to the fact that we are assuming ∂Ω ∈ C 3 we have the following Helmholtz-Weyl decomposition of u ∈ L r (Ω), r ∈ (1, ∞), see [10, Theorems III.1.2 and III.2.3]: where q ∈ W 1,r (Ω)/R is unique, The existence of the Helmholtz-Weyl decomposition (13) implies the well-posedness of the following problem [10]: The space H 0 r (div) from the Helmholtz-Weyl decomposition admits the following characterization for simply connected domains Ω for which ∂Ω ∈ C 2 , see [5, Corollary 2.3].
This characterization yields the following result.
Corollary 1 (gradient estimate). In the setting of Lemma 2.3, if u ∈ H 0 r (div) is such that rot u ∈ H 0 r (div), then u ∈ W 1,r (Ω) and ∇u L r ≤ c rot u L r . Proof. Note that rot u ∈ H 0 r (div) implies rot u = rot v with v as in (15), satisfying ∇v L r ≤ c rot u L r .
Define e = u − v, and notice that we have div e = div v, e · n = 0 on ∂Ω, and rot e = 0. This, in particular, implies that e = ∇ψ. Therefore, the function ψ solves the Neumann problem ∆ψ = − div v in Ω with n · ∇ψ = 0 on ∂Ω. Using a regularity result for this problem, see [ where we used that ∇v L r ≤ c rot u L r .
2.2. The Stokes problem with nonstandard boundary conditions. Consider first the classical Stokes problem, supplemented with no-slip boundary conditions: Recall the following regularity result, [10, Lemma IV.6.1], if ∂Ω ∈ C 2 and F ∈ L r (Ω) with r ∈ (1, ∞), then the solution to (16) satisfies u ∈ W 2,r (Ω), p ∈ W 1,r (Ω)/R, and Our analysis in section 3 will also require regularity results for the following non-standard Stokes-type problem, with a ∈ X and div a = 0. The problem with a = 0, α = 0 was discussed in [13].
We will now establish the necessary well-posedness and regularity for (18). A weak formulation reads: 3 2 (Ω). In light of the decomposition given in Lemma 2.2 the variables u and p can be decoupled and we can consider, instead, the following weak formulation: find for all ψ ∈ H 0 (rot), q ∈ W 1,3 (Ω). To see this, it suffices to set, in (19), ψ = v ∈ H 0 (rot) and ψ = ∇q with q ∈ W 1,3 (Ω) and recall that (12) implies that H 0 (rot) ⊂ L 3 (Ω).
We will now prove the well-posedness of (19) and (20). This is done as an auxiliary step towards showing that (5) is well-posed for α large enough (if the other data is fixed). For this reason, we need to make sure that certain constants in our estimates are independent of α ≥ 0. Since this includes two extreme cases: α = 0 and α → +∞, it is helpful to introduce the following parameter, where C P is the optimal constant in the Poincaré inequality v ≤ C P ∇v , v ∈ W. We introduce it in the definition of α + so that it has the proper (physical) scaling. The definition of α + allows us to obtain the following simple bound: To check the last inequality in (22) it is helpful to note the trivial bound for non-negative reals a, b, c, d. Notice, in addition, that by considering two cases α ≤ C −2 P ν and α > C −2 P ν and also using Poincare's or Young's inequality, we obtain that v ∇v ≤ α These two inequalities will help us to obtain conditions to guarantee that problems (19) and (20) then problems (19) and (20) are well-posed and have the same unique solution.
Proof. We begin the proof by noting that all bilinear forms in (19) and (20) are well-defined and continuous on the corresponding spaces. In particular, we have The Gagliardo-Nirenberg interpolation inequality [1,Theorem 4.17] provides where in the last step we applied (9). Now (26)-(27) yields the bound for u ∈ H 0 (rot), where the constant C depends only on Ω. In view of estimate (26) for 0 ≤ α < α + and (28) for α = α + , respectively, we see that there is a constant C * = C * (Ω) such that, whenever (25) holds, we have, for every u ∈ H 0 (rot), that where the constant c > 0 depends only on Ω. With these estimates at hand, we can now show the well posedness of each problem. We begin with (20) since it is simpler. Notice that, in light of (29) and the Banach-Nečas-Babuška theorem [8, Theorem 2.6], the u-problem is well posed. Now since u ∈ H 0 (rot) is uniquely defined, estimate (26) together with the well posedness of (14) show that the p-problem is well-defined as well. In addition, estimates (22) and (29) yield From the second equation in (20), (26) and (30) we conclude that, if (25) holds, we have where, for the last inequality, we used that ν − 1 2 α + . We now proceed with (19). For that we begin by noticing that both spaces H(rot) ∩ L 3 (Ω) and W 1, 3 2 (Ω)/R are Banach and reflexive. Hence the well-posedness of (19) follows from the theory of saddle point problems as detailed, for instance, in [8,Theorem 2.34]. Indeed, inequality (29) gives the inf-sup property and nondegeneracy of the u-form over the kernel of the p-form, which happens to coincide with H 0 (rot) ∩ L 3 (Ω). On the other hand, the decomposition (11) and the Helmholtz-Weyl decomposition of L 3 (Ω) yield the inf-sup property of the p-form: It remains to show that the solutions coincide, but this is immediate upon choosing appropriate test functions.
Remark 2 (large α). Notice that, due to the definition of α + , for any given viscosity ν and vector a ∈ X with div a = 0, there exists α large enough so that condition (25) is satisfied.
We now establish a regularity result for the weak solution of (19) and (20). Although relatively straightforward, such a result is not found in [13] nor seemingly anywhere else in the literature.
Proof. By definition, the solution to (20) satisfies (u, p) ∈ H 0 (rot) × W 1, 3 2 (Ω). Thanks to the embedding (12), we have that div u = 0 in Ω and u · n = 0 on ∂Ω. We want to apply the regularity result in (10) and so it remains to show rot u ∈ H 1 (Ω) together with a suitable bound on rot u H 1 .
The estimates above and the second equation in (20) imply the claimed regularity for the pressure: ∇p ∈ L 2 (Ω). Moreover, (30) and (34) yield the estimate with a constant that depends only on Ω.
It remains to show the regularity of u. Since ∇p, g, a × rot u, αu ∈ L 2 (Ω) we have rot rot u ∈ L 2 (Ω). Recalling that n · rot u = 0 on ∂Ω we have rot u ∈ H(rot) ∩ H 0 (div) and so, owing to the regularity of ∂Ω, rot u ∈ H 1 (Ω), as we intended to show. Moreover, we have the estimate where C only depends on Ω.
This completes the proof.
For given α ≥ 0 and a ∈ X, with a solenoidal and satisfying (25) we denote by ∆ −1 a : L 2 (Ω) → H 2 (Ω)∩W the (velocity-)solution operator to the Stokes problem (18). Owing to Lemma 2.5, ∆ −1 a is well-defined. We will denote by ∆ −1 0 the solution operator with a = 0, α = 0. This operator has the following properties: a w = 0 on ∂Ω, and for allw ∈ W, s.t. div w = 0, The first three are obvious. For the fourth one we note that w ∈ W, s.t. div w = 0 yields w ∈ H 0 (rot). Now the identity follows by taking g = w and ψ = w in (20). The last one follows by taking g = w and ψ = ∆ −1 a w in (20).

2.3.
Additional forms and their associated bounds. As a last preliminary step, we define some forms that will be needed for the analysis of the velocityvorticity formulation (5). We first define the bilinear form f bc : The form f bc is well-defined on Q × W and continuous, since for any smooth P and χ it holds that We also note the identities which implies that, for any a ∈ X with div a = 0 such that ∆ −1 a is well defined, we have Finally, define the trilinear form b : and note that, whenever div u = 0, this form is skew-symmetric, i.e. b(u, v, v) = 0. We will utilize the following bounds for the nonlinear terms that arise in our analysis.
with the last bound requiring div u = div w = 0.
Proof. The constant M will depend only on Ω. Each of the bounds in this lemma will hold with a potentially different constant that depends only on Ω, and we take M to be the maximum of these constants. Estimate (38) follows from Hölder's inequality, Sobolev inequalities and (27). For (40), to bound the first term, we use the fact that div u = 0, Hölder's inequality and the embedding H 1 → L 6 to obtain where in the last step we used (35) and (27). Similarly, for the second term, we obtain Adding these bounds gives the result claimed in (40). For (39), Hölder and Sobolev inequalities along with (27) and (35) provide which proves the stated result.
3. Weak formulation of the velocity-vorticity problem in the presence of solid walls and its well-posedness. In this section, we will prove that the velocity-vorticity formulation with no-slip velocity boundary conditions and the corresponding vorticity boundary conditions (5) is well posed. Our method of proof will be as follows: we will obtain a priori bounds which, through a fixed point argument, will allow us to show existence. In addition, provided the data satisfies further restrictions, we will be able to show uniqueness. The highlight of our approach lies in that our arguments hinge on tools developed in section 2.1 and the problems studied in section 2.2.
The weak formulation of (5) is obtained in a standard way, that is multiplying by suitable test functions and integrating by parts. It reads: given f ∈ L 2 (Ω), find (div w, λ) = (div u, π) = 0, (42) We note that the term ∂Ω (f ×n)·χ ds coming from the vorticity boundary condition (4) cancels, when we integrate by parts in (rot f , χ).

3.1.
A priori estimates and uniqueness. We begin the analysis of (41)-(44) with a priori bounds on solutions.
Lemma 3.1 (a priori bounds). Assume that system (41)-(44) has solutions. Then they satisfy where the norm · * was defined in (23) and α + in (21). If, in addition, the forcing term f satisfies where the constant C * is the same as in (25) of Lemma 2.4, then they verify with K 2 = K 2 (ν, α + , f , Ω) and K 3 = K 3 (ν, α + , f , Ω). The constants K 2 and K 3 , however, remain uniformly bounded as α → ∞, while all the other data of the problem remain fixed.
Proof. We begin with the velocity bound. Setting v = u in (41) and using (22) yields estimate (45).
Notice now that, if (46) holds, then (45) implies that a = u satisfies condition (25) of Lemma 2.4 and, consequently, the operator ∆ −1 u is well defined. We set χ = ∆ −1 u w in (44) and observe that, due to (37) and (35), the η-term and boundary functional vanish, thus (44) reduces to which when compared with (35) yields To obtain the velocity-pressure part of estimate (47) we employ a bootstrapping argument. Hölder's inequality, the embedding H 1 → L 6 , the bound (45) on the velocity, and the L 2 -bound on the vorticity yield the existence of a constant C(Ω) that depends only on the domain Ω, for which and, consequently w × u ∈ L 3 2 (Ω). We now apply the regularity result for the Stokes problem (17) using where we used (45) and the fact that α ≤ α + . The estimate given above shows that the velocity part of the solution to (41)-(44) satisfies We now invoke that, for every r < ∞, we have the embedding W 2, 3 2 → L r to obtain where, in the last step, we used (48) and the L 2 -estimate on the vorticity. This shows that F ∈ L 7 4 (Ω) with the bound F Using, once again, (17) yields that u ∈ W 2, 7 4 (Ω) with the estimate u W 2, 7 4 ≤ ν −1 C 2 (ν, α + , f , Ω). Finally, we use the embedding W 2, 7 4 → L ∞ to assert that This gives us that F ∈ L 2 (Ω) with the estimate so that invoking, one last time, (17) we obtain Let us now bound the H 1 -norm of the vorticity and the L 2 -norm of η. Setting χ = w in (44) yields α w 2 + ν rot w 2 + ν div w 2 = b(w, u, w) + (f , rot w) − f bc (P, w).
Applying (9) and then using the continuity of the functionals on the right hand side we obtain To control the trilinear term, we use (47) and the embedding H 2 → W 1,6 to obtain and, as a consequence, follows. An application of (7) yields the desired bound for η.
It remains then to show that K 2 and K 3 are bounded as α grows large. To see this, we first observe that for α sufficiently large we have α + = α, and so we need to study the dependence on α + . In the course of the proof of estimates (47), we obtained that Since K 1 does not depend on α + , it then is sufficient to show that K 2 remains bounded as α → ∞. From (49) we have where we successively applied the definitions of C j , j = 0, 1, 2. The only power of α + that appears is negative and so we are able to conclude.
Remark 3 (large α). If all the problem data besides α is kept fixed, the a priori bounds of Lemma 3.1 will all hold if α is taken sufficiently large. Indeed, this follows from the definition of α + , as for a given f it is always possible to find α large enough so that condition (46) is satisfied.
Remark 4 (large α). As in Lemma 3.1, if all the data besides α is fixed, Theorem 3.2 implies that uniqueness can be obtained by taking α sufficiently large.

3.2.
Existence. To prove the existence of solutions, we will utilize the following fixed point theorem, referred to as Shaefer's fixed point theorem in [9] and as Leray-Schauder's fixed point theorem in [12]. find y λ ∈ Y satisfying y λ = λF (y λ ), 0 ≤ λ ≤ 1, are uniformly bounded. Then the problem y * = F (y * ) has a solution y * ∈ Y .
We will proceed by constructing a compact map whose fixed points are solutions of (41)-(44), then consider the family of fixed point problems, and finally apply Lemma 3.3.
Proof. The bounds follow standard arguments. The first bound repeats the proof of (45), while the second is (17) for r = 2. Once this is established, we invoke (9) to conclude the last one.
We are now ready to show the existence of solutions.
Remark 5 (well-posedness). We comment that the conditions of both Theorem 3.2 and Theorem 3.5 are satisfied if α is large enough. Thus, in this case, we conclude that the problem is well-posed. 4. Conclusions and outlook. We have proven well-posedness of a steady velocityvorticity system with no-slip velocity boundary conditions, a corresponding zeronormal condition for vorticity on the boundary, and a natural boundary condition for vorticity involving a pressure functional. The well-posedness result provides a mathematical foundation for numerical methods, which exploit vorticity equations with these boundary conditions, and suggests a possible framework for their numerical analysis.
Significant technical difficulties arose due to the use of vorticity, vorticity boundary conditions, and the boundary conditions containing a pressure functional, and were overcome by a long, technical analysis that combined and extended analyses from [13,14,5]. These results are important for at least two reasons in addition to that stated above: First, there are very few analytical results known for fluid problems involving vorticity boundary conditions. Second, several results in section 2 are new, and can aid in future works for related PDEs which use non-standard boundary conditions.
Moving forward, these results will allow for improved analysis of numerical methods for NSE in velocity-vorticity formulations, and thus likely also to improved algorithms. To date, analysis of numerical schemes for these systems has been limited to 2D or the 3D steady case, and without physically derived boundary conditions, e.g. [2,21,19].