Weak Versus Strong Wall Boundary Conditions for the Incompressible Navier-Stokes Equations

The pressure-velocity formulation of the incompressible Navier-Stokes equations is solved using high-order finite difference operators satisfying a summation-by-parts property. Two methods for imposing Dirichlet boundary conditions (one strong and one weak) are presented and proven stable using the energy method. Additionally, novel diagonal-norm second-derivative finite difference operators are derived with highly improved boundary accuracy. Accuracy and convergence measurements are presented and verified against theoretical expectations. Numerical experiments also show that subtle effects close to solid walls are more efficiently captured with strong boundary condition imposition methods rather than weak (less degrees of freedom required).


Introduction
Multiple numerical methods for the incompressible Navier-Stokes equations have been suggested in the past.Among the popular ones are fractional step projection methods [1][2][3], where an approximation of the velocity field is computed and then corrected in a following step to impose the divergence free condition.Used together with staggered grids [4] this type of method has achieved great success and is probably the most commonly used today.However, although flexible (in principle any method can be used for the spatial discretization) and efficient (the velocity and pressure computations are decoupled), it exhibits some issues in terms of temporal accuracy and boundary conditions.Another common approach is to rewrite the problem in terms of the streamfunction and vorticity of the velocity field, which results in the streamfunction-vorticity formulation [5].This approach has achieved

Continuous Analysis
The dimensionless incompressible Navier-Stokes equations are given by where V is the velocity field, p the pressure, Re the Reynolds number, G the boundary data and H the initial data.In two spatial dimensions we have V = (u, v) T , where u and v are the velocities in the x-and y-direction, respectively.The boundary operator B and boundary data G determine the boundary conditions.The spatial domain is ∈ R 2 and its boundary ∂ .We refer to the first equation in (1) as the momentum equation and the second equation as the divergence free condition.
Taking the divergence of the momentum equation leads to where φ = ∇ • V.Here subscripts denote partial differentiation.Substituting the divergence free condition φ = 0 leads to the pressure Poisson equation (PPE), Note that (3) by itself is not enough to replace the divergence free condition.This is seen by substituting (3) in (2) and not imposing φ = 0. We get i.e. the divergence satisfies an advection-diffusion equation.For φ to be non-growing boundary conditions for the divergence are required, here we use homogeneous Dirichlet boundary conditions.See [19] for a derivation of maximally dissipative boundary conditions to this problem.Additionally, for the divergence to be zero for all time the initial divergence has to be zero.Finally, we arrive at the velocity-pressure formulation of (1): Note that the divergence free condition in (1) has been replaced by the PPE and a homogeneous Dirichlet boundary condition for the divergence.

Energy of the Momentum Equation
Let (u, v) = uv denote the L 2 -inner product of two real-valued functions u, v and ||u|| 2 = (u, u) the corresponding norm.Taking the inner product of V and the momentum Eq. in (5) (the energy method) leads to where E is an energy defined by E = ||u|| 2 + ||v|| 2 and the boundary terms are given by where n = (n (x) , n (y) ) is the normal vector to the boundary.To arrive at (6) the following identities are used: From ( 6) it is clear that the energy is decreasing if the boundary terms are taken care by suitable (i.e. the correct number of) boundary conditions such that BT ≤ 0 (for homogeneous data) and the divergence free condition ∇ • V = 0 is fulfilled.

Characteristic Boundary Conditions
In this paper we use characteristic boundary conditions (CBC) as non-reflecting outflow boundary conditions and as interface conditions for the backwards facing step problem in Sect.4.3.They are derived for the incompressible Navier-Stokes equations in [19], look there for more details.Boundary characteristic variables to the incompressible Navier-Stokes equations are given by where , λ 3 = 0, and and V n = un x + vn y denotes the normal velocity and V t = −un y + vn x the tangential velocity.Since λ 2,5 < 0 the CBC are where g 2 and g 5 are boundary data.With these boundary conditions the energy Eq. ( 6) becomes where With g 2,5 = 0 it is clear that these boundary conditions leads to additional energy dissipation.

Dirichlet Boundary Conditions
If we use Dirichlet boundary conditions for no-slip walls and inflow boundaries, then an energy estimate can easily be obtained if the divergence free condition is fulfilled.Assuming zero data or a wall (V = G = 0) the boundary terms (7) are zero and the energy Eq. ( 6) becomes 3 Discrete Analysis

Definitions
The spatial discretization is done using high order finite difference SBP operators.Here we include a short description of the operators required for this problem.For more details on SBP finite differences, see [8,9,20,21].The two-dimensional domain is discretized using m points in each dimension (the extension to rectangular domains is straightforward): where x and y are vectors containing the grid points.For an equidistant grid we have where h x,y are the grid sizes and x l,r , y l,r are the domain limits.Let v denote a semidiscrete solution vector so that v T = [v (1) , v (2) , . . ., v (m) ], where m ] contains the solution points at x i along the y-axis, see Fig. 1.Here we have three such solution vectors: u, v and p corresponding to the discrete values of x-velocity, y-velocity and pressure respectively.Note that u, v and p previously denoted the space-continuous The incompressible Navier-Stokes equations require first and second derivative SBP operators, the following definitions are central: x, using a pth-order accurate narrow-stencil, is said to be a pth-order accurate narrow-stencil diagonal-norm first derivative SBP operator if H is diagonal and positive definite and Q + Q T = B = diag(−1, 0, . . ., 0, 1).
using a pth-order accurate narrow-stencil, is said to be a pth-order accurate narrow-stencil diagonalnorm second derivative SBP operator if H is diagonal and positive definite, M is symmetric and positive semi-definite, S approximates the first derivative operator at the boundaries and B = diag(−1, 0, . . ., 0, 1).
The matrix B can be expressed in terms of the frequently used vectors as B = −e 1 e T 1 + e m e T m .Additionally, the matrix S is often written in terms of one sided first-derivative operators d 1 and d m , approximating the first derivative at the left and right boundary respectively, as B S = −e 1 d 1 + e m d m .The extension of the one-dimensional SBP operators to two dimensions is done using the following relationships: where ⊗ denotes the Kronecker product and I m is the m × m identity matrix.The discrete Laplace operator is given by We also define the following discrete inner products and corresponding norms: where H = H x H y = (H ⊗ H ) and H = H 0 0 H .In this paper we consider two types of narrow-stencil diagonal-norm SBP operators: traditional and optimal.The main difference is in the distribution of grid points.The traditional operators are designed to act on standard equidistant Cartesian grids.And the optimal operators on non-equidistant grids where the distribution of grid points near the boundaries are chosen to maximize the accuracy of the operators.See [12,18] for more details.Note that we do not map the PDE in the usual sense, we simply discretize it on different grids and use operators designed for those grids.In the present study we use 6th and 8th order accurate operators of both types.

Compatible Second Derivative Operators
The derivation of the discrete PPE (see Sect. 3.5) necessitates that the discrete Laplace operators in the momentum equation are given by the discrete divergence of the discrete gradient operator, i.e.

D (w)
The superscript in (21) indicates that this is a wide operator (wide internal stencil) constructed from the wide second derivative operator However, since central first derivative operators do not resolve the highest frequency mode (the π-mode) that can exist on the grid, the induced second derivative operator D (w) 2 will not either.This can lead to odd-even decoupled oscillations in the solution [22].To rectify this we combine D (w) 2 with artificial dissipation to obtain narrow-stencil second derivative operators that are compatible.The following definition, first introduced in [22], is central: 2 be pth-order accurate narrow-stencil diagonal-norm first-and second-derivative SBP operators.If and the remainder R ( p) is positive semi-definite, D 1 and 2 are called compatible.

The narrow discrete Laplace operator D (n)
L , constructed according to (19) using 2 , is given by where The remainders R ( p) for p = 6, 8 were first derived for constant coefficient traditional operators in [22].To find them for the optimal operators we take inspiration from [23], where narrow-stencil variable coefficient second-derivative SBP operators were derived, and make the ansatz 8 . ( The internal stencils in are narrow-stencil approximations of the ith derivative using the same non-equidistant grid point distribution as D 1 .Note that R ( p) is symmetric and positive semi-definite by construction and that the term −H −1 R ( p) in ( 22) constitutes artificial dissipation (since it is an even derivative approximation in the interior).To achieve non-stiff second derivative SBP operators, we zero out the first few rows in D ( p) i .For more details see [22,23].See also [24], where a similar technique is used to construct general artificial dissipation operators for SBP operators.
The optimal SBP operators, including the new second-derivative operators, are available online at https://github.com/guer7/sbp.

Remark
The compatible second derivative SBP operators can be made more or less narrow by including or excluding terms in (25), each term corresponds to canceling out one additional value on each side of the inner stencil.The choices in (25) are those that yield second-derivative operators with minimal stencil width, referred to as narrow-stencil second-derivative operators.Numerical experiments have shown that excluding the terms corresponding to the lowest order derivatives in (25) has some favorable properties.Compared to the minimally narrow scheme, this choice was found to be slightly more accurate and allow for slightly larger time steps.Presently, the reason for this is not fully understood.It may be related to the accuracy of the boundary closures, but a thorough investigation is out of scope of the current work.The numerical results presented in Sect. 4 utilize this slightly wider operator, since it leads to a more efficient scheme.

Semi-discrete Momentum Equation
Let define the semi-discrete velocity field.Ignoring boundary conditions for now, a consistent semi-discrete approximation of the momentum equation is given by where the right-hand-side is given by Here are the skew-symmetric split advection operators (see (8)) and ū, v, ūx and vy denote diagonal matrices with the elements of u, v, D x u and D y v respectively on the diagonal.

Discrete Energy Estimate
To make the stability analysis more readable, we only present boundary terms corresponding to the western and eastern boundaries.The southern and northern boundaries are done analogously.Taking the inner product of w and ( 27) and using Definitions 1 and 2 gives where is the discrete energy, the dissipation (the second row in (32) is due to the artificial dissipation AD ( p) ), the boundary terms where ūw,e denote diagonal matrices with the elements of u w,e on the diagonal and the remainder terms proportional to the discrete divergence.The energy Eq. ( 30) is the discrete equivalent to the continuous Eq. ( 6).To obtain a discrete energy estimate we need to fulfill two things: (i) impose the correct number of boundary conditions such that BT ≤ 0 (with zero boundary data); and (ii) ensure that the discrete divergence converges towards zero as the grid is refined, so that RT = 0 as h → 0.

Semi-discrete Pressure Poisson Equation
Regarding the PPE one has two options: (i) derive the continuous PPE (3) first and then discretize; or (ii) discretize the momentum equation first and then derive a discrete PPE.
To ensure that the semi-discrete momentum equation and PPE are compatible, we take the second approach.
As in the continuous analysis, the first step in deriving a semi-discrete PPE is to multiply (27) by D x D y from the left (take the discrete divergence of the semi-discrete momentum equation).We get where φ = D x u + D y v denotes the discrete divergence and The following lemma is one of the main results of this paper: The norm of the discrete divergence φ is non-increasing if the semi-discrete momentum Eq. ( 27) holds, the discrete divergence satisfies and the pressure satisfies the discrete PPE Proof Taking the inner product of φ and (35) and using Definitions 1 and 2 lead to where If φ w,e,s,n = 0 and D

(w)
L p = F the right-hand-side of (39) is zero and Lemma 3.1 shows that the H -norm of the discrete divergence will decrease over time if the boundary condition (37) holds.For the discrete divergence to be zero at all times, we also require that the initial condition satisfies Remark The derivation of the discrete PPE leads to a wide Laplace operator in (38).One may consider using the narrow operator instead since it includes dissipation that does not allow for π-mode oscillations.However, numerical experiments have shown that this has no clear benefits in practice.Compared to the narrow options, the scheme with the wide operators is slightly more accurate and allows for slightly larger time steps.As for the Laplace operator in the momentum equation, the differences are minor and an extensive study of the reasons for this is beyond the scope of the current work.But we point out that the derivation of the discrete PPE presented in Sect.3.5 directly leads to a wide Laplace operator.Adding the dissipative term to construct the narrow operator makes the discrete PPE in one sense inconsistent with the discrete momentum equation.Which could be a reason for the small decrease in accuracy.

Semi-Discrete Problem
We are now ready to formulate the semi-discrete problem: Here D(u, v, p) is given by ( 28), F is given by (36) and B, Ĝ and Ĥ define the discretized boundary and initial conditions.In this paper we consider two methods to solve (43) with Dirichlet boundary conditions: projection and SAT.Both methods augment the momentum equation and PPE to impose the boundary conditions, resulting in provably stable systems of ODEs.In Sects.3.7 and 3.8 the two methods are presented separately, including procedures for deriving and imposing consistent pressure boundary conditions.

Remark
We use CBC as outflow and interface conditions for the numerical experiments on the backwards facing step problem (presented in Sect. 4. 3).An SBP-SAT discretization of these conditions including energy estimates is presented in [19].

Projection Method
In this section SBP-P discretizations of the momentum and pressure Poisson equations are presented.In Sect.3.7.3 the method is summarized step-by-step for wall boundary conditions.

Momentum Equation
The Dirichlet velocity and divergence boundary conditions can be represented discretely as where Here L w,e w = g w,e corresponds to the Dirichlet boundary conditions and L div w = 0 to the divergence boundary condition on the western and eastern boundaries.A consistent semidiscrete approximation of the momentum equation with (44) imposed using the projection method is given by where û v = ŵ = Pw + (I − P) g.Here is a projection operator, g satisfies L g = g and I is the identity matrix.Note that P is self-adjoint with respect to the inner product (•, •) H , i.e.
and that g never has to be computed explicitly since it only appears as L g in the equations (the same is true for gt ).See [13] for more details on the projection method.
The following lemma is one of the main results in this paper: Lemma 3.
where ūw,e denote diagonal matrices with the elements of ûw,e on the diagonal.Since L ŵ = L Pw = 0, and thus ûw,e = vw,e = 0, we get BT = 0.
Note that with the projection method both the velocity and divergence boundary conditions are imposed simultaneously, thus solving (47) with a pressure satisfying D L p = F is sufficient to ensure that the divergence is zero.

Remark
The projection method can be used to project zero divergence everywhere (not only on the boundaries), this may be useful as pre-or post-processing step.Then the divergence boundary operator becomes (51)

Pressure Equation
If no pressure boundary data is available, consistent pressure values have to be derived numerically.Here we extract the gradient of the pressure on the boundaries from the normal component of the momentum equation.For example, on the western boundary we have The validity of these boundary conditions has been a topic of discussion for a long time, see for example [25][26][27][28][29].In our view no consistency requirements are violated as long as the velocities in (52) satisfy the divergence boundary conditions.The Neumann pressure boundary conditions are given by Lp = g p , where The PPE with Neumann boundary conditions imposed using the projection method is given by where σ is a positive scalar, L gp = g p and The term σ H (I − P)( p − gp ) in ( 54) is required to remove zero eigenvalues due to the projection.A similar term is used in [13] to rectify inconsistent initial and boundary data for PDEs.The parameter σ should be chosen proportional to (h x h y ) −1 so that the term σ H (I − P)( p − gp ) does not go to zero as we grid refine.From numerical experiments (not presented here) we found that is a good choice in terms of accuracy and robustness.By rearranging (54) we get the system where Note that A N is singular (even with σ = 0).Since any constant vector p is a solution to A N p = 0 (since D (w) L p = Lp = 0), any solution to (57) is unique up to an arbitrary constant.To avoid this singularity we impose the additional constraint that the average pressure across the whole domain is zero.A simple numerical formulation of the resulting system is where 1 is a vector with all components equal to one [30].

Summary
We summarize the projection approach by describing the procedure for computing the righthand side of (47) with wall boundary conditions:

SAT Method
In this section SBP-SAT discretizations of the momentum and pressure Poisson equations are presented.In Sect.3.8.3 the method is summarized step-by-step for wall boundary conditions

Momentum Equation
The ODE system (27) augmented with Dirichlet boundary conditions imposed using SAT is given by where and Here σ (u,v) w,e ∈ R m×m are chosen so that an energy estimate is obtained.Note that the pressure terms in (61) are necessary to obtain an energy estimate due to the pressure terms in (33).
The following lemma is one of the main results of this paper: Lemma 3. 3 The system (60) is a stable approximation of the momentum equation with Dirichlet boundary conditions if the discrete divergence free condition is fulfilled (so that RT = 0), Proof Taking the inner product of w and (60) with zero data and using Definitions 1 and 2 leads to boundary terms Choosing yields BT = 0. Thus, if RT = 0, an energy estimate is obtained and ( 60) is a stable ODE system.
To impose the divergence boundary conditions with SAT we employ the same technique as in [19].With this approach the boundary values of the velocity in the normal direction are modified so that the divergence free condition (37) is fulfilled exactly.For example, using Matlab index notation, on the western boundary we strongly enforce (66)

Pressure Equation
An issue with the SAT method for Dirichlet velocity boundary conditions is the pressure data required in (61).As a consequence (52) can not be used directly as pressure boundary data since it only provides the normal derivative.Here we use a numerical trick where the current pressure is used to compute the boundary pressure values so that the normal derivative equals g given by (52).For example, using Matlab notation, on the western boundary we use The PPE with Dirichlet boundary conditions imposed using SAT is given by where We get the system where n . (71)

Summary
We summarize the SAT approach by describing the procedure for computing the right-hand side of (60) with wall boundary conditions: backwards-facing step problem and compare the results to experimental benchmark data.
To integrate the ODE systems (47) and (60) in time we use the classical 4th order explicit Runge-Kutta method with time step t ∝ h 2 Re , where h is the minimum spatial grid interval.The pressure systems (59) and (70) are solved using LU-factorization.

Convergence Study
An analytical solution to the incompressible Navier-Stokes equations in two dimensions is the Taylor-Green vortex flow: where (x 0 , y 0 ) are the starting coordinates of the vortices, θ is the angle at which they are moving and u ∞ is their speed.We use a square domain x, y ∈ [−1, 1] with total degrees of freedom N = m 2 where m is the number of grid points in each dimension.The problem parameters in (72) are chosen to Re = 100, u ∞ = 1, (x 0 , y 0 ) = (0, 0) and θ = π 4 .The initial and boundary data are given by the analytical solution (72) (for this problem we use exact pressure data).The errors of the approximations are measured in the H -norms given by (20) of the difference between approximated and analytical solution, denoted (w, p) , at t = 1.The convergence rate is estimated as where N 1 = N 2 are the degrees of freedom of two separate solutions and d is the spatial dimension, here d = 2.The theoretical rate of convergence for this problem is min ( p, r + 2), where p is the internal stencil order and r the boundary closure order of the least accurate SBP operator in the scheme [31].With a diagonal norm the closure order of the first derivative operator is half that of the internal order, r = p/2.Furthermore, the wide-stencil second derivative operator D (w) 2 yields a scheme with one order of accuracy less at the boundaries, r = p/2 − 1 [9].The decreased boundary order of the second derivative operators is not rectified by the artificial dissipation discussed in Sect.3.2, and thus it limits the convergence rates of the full schemes to min ( p, p/2 + 1) for both the traditional and optimal SBP operators.We get for the 6th and 8th order operators ( p = 6, 8) the theoretical convergence rates 4 and 5 respectively.
In Fig. 2 the velocity and pressure errors are plotted versus grid resolution for Dirichlet boundary conditions imposed using projection and SAT.The results show that projection and SAT are very similar in terms of accuracy and achieve close to or slightly above the expected convergence rates.Furthermore, the results demonstrate a significant difference in efficiency between the traditional and optimal SBP operators.For example, the optimal 8th order SBP operators yield a more accurate solution with 51×51 grid points than the traditional 8th order SBP operators with 151 × 151 grid points.We also emphasize that the gain in accuracy when Fig. 2 Error versus grid resolution for the Taylor-Green vortex flow problem solved with traditional and optimal 6th and 8th order SBP operators.Top row: Dirichlet boundary conditions imposed using the projection method.Bottom row: Dirichlet boundary conditions imposed using SAT.The dashed lines indicate the theoretical convergence rates 4 and 5 for the 6th and 8th order SBP operators respectively.The labels indicate the convergence rate of each method and operator, estimated by (73) using the two finest grid resolutions going from 6th to 8th order operators is much higher with the optimal operators compared to the traditional, this is particularly clear for the pressure.

Lid-Driven Cavity
The lid-driven cavity is one of the most common benchmark problems for viscous incompressible flows.The problem setup consists of a two-dimensional unit square domain with walls at each boundary where the top wall is moving at a constant tangential velocity, see Fig. 3.The velocity field and pressure are initiated from zero and simulated until the velocity field has approximately converged towards a steady-state solution.We define the stopping criteria as where w n is the velocity field at time step n and = 10 −8 .To avoid temporal discontinuities the non-zero boundary condition at the top wall is smoothly brought into the domain using a hyperbolic tangent function of time.All results presented here are obtained using the optimal 8th order SBP operators with m grid points in each dimension.In Fig. 5 the streamlines of steady state solutions to the lid-driven cavity problem with Dirichlet boundary conditions imposed using projection (left column) and SAT (right column) are plotted.With the projection method the secondary corner vortices are visible already at m = 31, whereas with the SAT method the corner vortices require significantly more degrees of freedom to develop, even with m = 151 the vortex in the bottom left corner is not clear.This suggests that strong enforcement of physical boundary conditions (projection) is preferable over weak enforcement methods (SAT) when effects close to the boundaries are of primary interest.We hypothesize that this is a consequence of these effects requiring a clear distinction between boundary and inner domain to develop.Therefore, since the accuracy of the boundary conditions with SAT scales with the spatial resolution, one has to grid refine until this distinction is sufficiently large.Grid refinement with projection this is not required To verify the accuracy of the SAT and projection methods for the lid-driven cavity problem we compare the results with m = 151 to benchmark solutions provided in [32].In Fig.

Remark
We do not stretch the grid by mapping the PDE, which is usually done with finite difference methods to better resolve the boundary layer.The non-equidistant grid point distribution inherent to the optimal operators is in this case enough to capture the boundary effects.

Backward-Facing Step
The third problem we consider is the backward-facing step problem.The domain consists of an inflow pipe connected to an outflow with larger radius, forming a step where the flow is disrupted and interesting effects occur, see Fig. 6.Here we use H i = 0.5 and H o = 0.97 which determines the expansion ratio of the step to H o H i = 1.94.To solve this problem with the SBP finite difference operators used here the computational domain has to be split into at least three rectangular blocks (one for the inflow and at least two for the outflow).Here we couple the blocks together by weakly imposing (using SAT) that the ingoing characteristic variables are equal on each side of the interfaces.This procedure has proven to be accurate and stable in the past, see for example [17].
Since we solve a local PPE in each block we benefit greatly in terms of efficiency by utilizing many smaller blocks as compared to a few larger ones.Therefore we set a fixed block size in the outflow channel of width 4 and height 0.5 or 0.47 (for top and bottom block respectively) with 168 × 21 grid points in each direction.The length of the outflow channel is determined by the number of such blocks.The inflow consists of one block of width 1 and height 0.5 with 42 × 21 grid points in each direction.As boundary conditions for the inflow we specify the normal velocity to be the theoretical solution of fully developed laminar flow in an infinitely long pipe, given by where V max is the velocity in the center of the pipe and R is its radius.Homogeneous characteristic boundary conditions imposed using SAT are used at the outflow.The Dirichlet boundary conditions (inflow and walls) are imposed using the projection method.In Fig. 7 the boundary and interface conditions are presented for a five-block domain.As for the liddriven cavity problem we use the optimal 8th order SBP operators and a stopping criteria given by (74).In Fig. 8 the velocity field of the steady-state solution to the backwards-facing step problem with Re = 600 is shown.The primary and secondary vortices are clearly visible in the streamlines plot.Additionally, one can see the tendencies towards a parabolic flow in the outflow channel.
We verify the accuracy by comparing the reattachment and detachment lengths, denoted x 1,2,3 in Fig. 6, of our results to experimental data provided in [33].For these results the outflow pipe is chosen sufficiently long so that the considered lengths are unaffected.In Fig. 9 the comparison is presented for Reynold's numbers up to Re = 1000 with good agreement up until the secondary vortex is developed at around Re = 400.present study.Circles: experimental data from [33] As pointed out in [33], above approximately this Reynold's number the flow becomes three-dimensional and thus the modeling error of the two-dimensional equations becomes dominant.

Conclusions
A high-order SBP finite difference method is developed for the velocity-pressure formulation of the incompressible Navier-Stokes equations.Two boundary procedures for imposing Dirichlet boundary conditions are presented based on the projection and SAT methods.The resulting schemes are proven to be stable using the energy method.Furthermore, new optimal second-derivative SBP operators are derived and shown to be highly efficient compared to traditional SBP operators.
The methods are evaluated on three benchmark problems: the Taylor-Green vortex flow, the lid-driven cavity and the backwards-facing step.The accuracy and convergence properties of the schemes are verified for the traditional and novel optimal SBP operators.For the liddriven cavity and backwards-facing step problems good agreement with benchmark data is achieved.Additionally, imposing wall boundary conditions strongly (projection) rather than weakly (SAT) is shown to be more efficient in the sense that corner vortices for the lid-driven cavity problem requires less degrees of freedom to appear.
Interesting future extensions include expanding the methods to curvilinear grids in three spatial dimensions and proving multiblock stability using the projection method.

Fig. 1
Fig. 1 Computational domain showing the orientation of the solution vectors

e w = (e 1 ⊗
I m ), e e = (e m ⊗ I m ), e s = (I m ⊗ e 1 ), e n = (I m ⊗ e m ), ) with L w,e = e w,e 0 0 e w,e and L div = d w D 1 e w d e D 1 e e .(

Fig. 3 Fig. 4
Fig. 3 Boundary conditions for the lid-driven cavity problem.The primary vortex in the center of the cavity and the secondary vortices in the bottom corners are indicated

Fig. 5
Fig. 5 Steady-state solutions to the lid-driven cavity problem with Re = 100 and m = 31, 91, 151, where m × m is the number of degrees of freedom.Left column: using projection to impose Dirichlet boundary conditions.Right column: using SAT to impose Dirichlet boundary conditions

Fig. 6 Fig. 7
Fig. 6 Domain for the backward-facing step problem.The in-and outflow velocity profiles and recirculation zones are indicated

Fig. 8
Fig. 8 Steady-state solution to the backward-facing step problem for Re = 600 with expansion ratio 1.94.The number of unknowns are 29106.Top: streamlines of velocity field.Bottom: profiles of the x-component of the velocity, the dashed lines indicate zero velocity