High-Order Cartesian Cut-Stencil Finite Difference Solutions for Streamfunction-Vorticity Formulation of 2-D Navier-Stokes Equations

Mohammadali Esmaeilzadeh¹, Ronald M. Barron^1,2*

¹Department of Mechanical, Automotive & Materials Engineering, Windsor, Canada.
²Department of Mathematics & Statistics, University of Windsor, Windsor, Canada.

DOI: 10.4236/am.2026.173011                

Abstract

The formulation and implementation of a high-order Cartesian cut-stencil finite difference method (CCST-FDM) to 2-D steady incompressible fluid flow in regular and irregular domains is considered in this paper. The CCST-FDM is capable of simulating flows in complex geometries by employing 1-D quadratic transformation functions to map any (uniform or non-uniform) physical stencil to a uniform computational stencil. In this work, the CCST-FDM is combined with compact high-order (HO) Padé-Hermitian approximations to produce HO CCST-FD schemes. Two different high-order (globally 4^th-order) accurate schemes are formulated for the CCST-FDM. Using the streamfunction-vorticity formulation of the Navier-Stokes equations, low and high-order solutions are computed for lid-driven flows in four different geometries, namely a square cavity and three irregular regions, namely L-shape, skewed quadrilateral and triangular cavities. Results from these geometries are compared to earlier studies for various Reynolds numbers. It is shown that the high-order CCST-FDMs can achieved higher accuracy on coarser grid than other high-order methods.

Keywords

Computational Fluid Dynamics, Cut-Stencil Finite Difference, High-Order Discretization, Streamfunction-Vorticity Formulation

1. Introduction

The equations of fluid motion can be derived from conservation principles for mass, linear momentum and energy, associated with continuity, Newton’s second law and the first law of thermodynamics, respectively. The resulting system of partial differential equations (PDEs) is often collectively referred to as the Navier-Stokes (N-S) equations [1]. Each equation in this system can be modelled as a convection-diffusion equation. For an incompressible flow, the energy equation can be decoupled from the system of equations. Then, the N-S equations for 2-D steady incompressible laminar flow in Cartesian coordinates $\left(x,y\right)$ , in non-dimensional and non-conservative primitive variable form, are given by Equations (1.1), where $u$ , $v$ denote the x and y components of velocity, $p$ is pressure, $\nu$ is kinematic viscosity and $Re=UL/\nu$ is the Reynolds number.

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ (1.1a)

$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+\frac{\partial p}{\partial x}=\frac{1}{Re}\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)$ (1.1b)

$u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+\frac{\partial p}{\partial y}=\frac{1}{Re}\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)$ (1.1c)

The conservation of mass Equation (1.1a) implies the existence of a function $\psi$ , defined by $\frac{\partial \psi }{\partial y}=u$ and $\frac{\partial \psi }{\partial x}=-v$ . Curves $\psi \left(x,y\right)$ = constant represent streamlines of the flow, and the function $\psi$ is referred to as the streamfunction. Taking curl of the velocity yields the vorticity $\left(0,0,\omega \right)$ where $\omega =\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$ .

Replacing $u$ and $v$ in terms of $\psi$ yields Equation (1.2a) for the streamfunction. Finally, taking the curl of the momentum equations (in vector form), one can eliminate the pressure and obtain Equation (1.2b) for vorticity.

$\frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial y^2}=-\omega$ (1.2a)

$u\frac{\partial \omega }{\partial x}+v\frac{\partial \omega }{\partial y}=\frac{1}{Re}\left(\frac{{\partial }^2\omega }{\partial x^2}+\frac{{\partial }^2\omega }{\partial y^2}\right)$ (1.2b)

Since there are only a limited number of practical exact solutions of the N-S equations, researchers have turned to numerical techniques to simulate relevant engineering flows. The ψ-ω formulation of the N-S equations has received considerable attention in the literature, particularly for the development and verification of new numerical schemes [2]-[5]. In this regard, three different numerical approaches, namely finite difference (FD) [1] [6], finite volume (FV) [7] [8] and finite element (FE) [9], each with their own strengths and weaknesses, dominate the discipline.

The primary limitation of the finite difference approach is the difficulties associated with the generation of a structured mesh for domains with complex geometry. However, recent research employing the 2-D Poisson equation has demonstrated that the Cartesian cut-stencil finite difference method (CCST-FDM) can successfully eliminate this limitation [10]. In this paper, the CCST-FDM is utilized to obtain high-order numerical solutions of the coupled nonlinear incompressible N-S equations. CCST-FDM utilizes a 5-point stencil, rather than a cell as in FV and FE, and relies on the construction of a unique 1-D quadratic mapping of each physical stencil, in any arbitrarily shaped domain, to a generic uniform computational stencil. This permits the implementation of standard FD approximations to solve PDEs in any complex shaped region. The formulation is also well-suited for combining CCST-FDM and high-order (HO) Padé-Hermitian methods to create higher-order cut-stencil schemes. This paper implements these schemes on uniform and clustered grid systems to compare solutions for low-order (1^st and 2^nd) and HO CCST-FD schemes applied to lid-driven cavity flows in regular and irregular geometries.

2. Essential Features of the Cartesian Cut-Stencil FDM

2.1. Physical and Computational Stencils: The Mapping

The Cartesian cut-stencil FD method facilitates discretization of PDEs in any arbitrarily shaped complex physical domain. This is accomplished by employing a unique localized mapping of each stencil in the physical domain to a generic computational stencil, as depicted in Figure 1.

Figure 1. Mapping of an arbitrary physical stencil to a uniform computational stencil.

In particular, the mapping illustrated in Figure 1 is defined by:

$x\left(\xi \right)=a_2{\xi }^2+a_1\xi +a_0$ (2.1a)

$y\left(\eta \right)=b_2{\eta }^2+b_1\eta +b_0$ (2.1b)

where the coefficients $a_i$ and $b_i$ $\left(i=0,1,2\right)$ can be expressed in terms of the coordinates of nodes P, W, S, E and N. The central point P and endpoints W, E, S and N are mapped to $p:\left(0,0\right)$ , $w:\left(-1,0\right)$ , $e:\left(1,0\right)$ , $s:\left(0,-1\right)$ and $n:\left(0,1\right)$ on the computational stencil, respectively [10] [11].

The N-S equations, in streamfunction-vorticity form (1.2), can be represented by general 2-D steady convection-diffusion equations for function $\overline{\varphi }$ , where $\overline{\varphi }$ is either $\psi$ or $\omega$ , i.e.,

$-\Gamma {\overline{\nabla }}^2\overline{\varphi }+\overline{U}\frac{\partial \overline{\varphi }}{\partial x}+\overline{V}\frac{\partial \overline{\varphi }}{\partial y}=\overline{S}\left(x,y\right)$ (2.2)

where ${\overline{\nabla }}^2$ , $\left(\overline{U},\overline{V}\right)$ , $\text{Γ}$ , $\overline{S}\left(x,y\right)$ are the Laplacian operator, convective velocity, constant diffusion coefficient and arbitrary source term, respectively. In the CCST-FD approach, the governing equation is written in terms of the computational stencil variables as:

$-\Gamma \left\{\frac{1}{{\left(x^{\prime }\right)}^2}\frac{{\partial }^2\varphi }{\partial {\xi }^2}-\frac{x^{″}}{{\left(x^{\prime }\right)}^3}\frac{\partial \varphi }{\partial \xi }+\frac{1}{{\left(y^{\prime }\right)}^2}\frac{{\partial }^2\varphi }{\partial {\eta }^2}-\frac{y^{″}}{{\left(y^{\prime }\right)}^3}\frac{\partial \varphi }{\partial \eta }\right\}+\frac{U}{x^{\prime }}\frac{\partial \varphi }{\partial \xi }+\frac{V}{y^{\prime }}\frac{\partial \varphi }{\partial \eta }=s\left(\xi ,\eta \right)$ (2.3)

where $\varphi =\varphi \left(\xi ,\eta \right)=\overline{\varphi }\left(x\left(\xi ,\eta \right),y\left(\xi ,\eta \right)\right)$ , etc. As mentioned above, the uniform arms of the computational stencil in both $\xi$ - and $\eta$ -directions allow implementation of standard FD approximations of all derivatives in the mapped form of the model equation. It can be shown that the Jacobian of transformation (2.1) is given by $J_P={x^{\prime }}_p{y^{\prime }}_p$ , which is never zero. The details and important concepts of this method have been discussed in [10] [11], which demonstrated the feasibility and advantages of this formulation for the numerical solution of the model PDE (2.2) for steady and unsteady convection-diffusion phenomena in engineering and applied sciences.

Since the mapping is confined to a 3-point stencil (in each direction), discretization of Equation (2.3) allows implementation of 3-point finite difference approximations (about p on the computational stencil) of the derivatives appearing in the equation. These approximations lead to a FD equation of the form

$a_p{\varphi }_p+a_w{\varphi }_w+a_e{\varphi }_e+a_s{\varphi }_s+a_n{\varphi }_n=$ 𝓈_p (2.4)

for which different discretization schemes produce different coefficients and source terms.

2.2. Discretization

Various FD approximations for the derivatives appearing in Equation (2.3) can be derived by expressing their value at p as a linear combination of values at the stencil endpoints. First and second derivatives of $\varphi$ (for m = 1 or 2) can be written as,

${\frac{{\partial }^m\varphi }{\partial {\xi }^m}|}_p={\displaystyle \sum _{nb}c{0}_{nb}^{\left(m\right)}{\varphi }_{nb}}+{\displaystyle \sum _{nb}c{1}_{nb}^{\left(m\right)}{\frac{\partial \varphi }{\partial \xi }|}_{nb}}+{\displaystyle \sum _{nb}c{2}_{nb}^{\left(m\right)}{\frac{{\partial }^2\varphi }{\partial {\xi }^2}|}_{nb}}$ (2.5)

where nb represents nodes p, w and e. A similar expression holds for $\eta$ -derivatives with nb representing nodes p, s and n. The coefficients in (2.5) are determined by the standard procedure of using Taylor series expansions about node P:$\left(x_p,y_p\right)$ , by noting that, for example, ${\overline{\varphi }}_E=\overline{\varphi }\left(x_P+\Delta x_2,y_P\right)$ . Applying this on the computational stencil, where $\Delta \xi =1$ , gives

${\varphi }_e={\varphi }_p+{\frac{\partial \varphi }{\partial \xi }|}_p+\frac{1}{2}{\frac{{\partial }^2\varphi }{\partial {\xi }^2}|}_p+\frac{1}{3!}{\frac{{\partial }^3\varphi }{\partial {\xi }^3}|}_p+\cdots$ . Similar series expansions can be derived

for ${\varphi }_w,{\varphi }_s,{\varphi }_n$ .

2.2.1. Low-Order Schemes

For 1^st- and 2^nd-order schemes, the diffusion terms in Equation (2.3) are discretized using 2^nd-order accurate 3-point central approximations, derived from Equation (2.5) by taking m = 2 and setting $c{1}_{nb}^{\left(2\right)}=0$ and $c{2}_{nb}^{\left(2\right)}=0$ . However, an upwind scheme for the convective terms is often necessary in numerical analysis of fluid flows to accurately capture the problem physics, especially in convection-dominated flows, and to stabilize the numerical solution since central differencing may not yield a correct or stable solution. Thus, various upwind differencing schemes have been proposed for the convection-diffusion equation [1] [12]. The 1^st-order upwind scheme for the first derivatives in the convection terms of Equation (2.3) is obtained from Equation (2.5) by setting m = 1, $c{1}_{nb}^{\left(1\right)}=0$ and $c{2}_{nb}^{\left(1\right)}=0$ . A 2^nd-order upwind scheme can be developed by taking m = 1 and $c{2}_{nb}^{\left(1\right)}=0$ . This 2^nd-order scheme is actually a Pade-Hermitian scheme involving the first derivatives at w, e, s and n. Using the resulting finite difference formulae to approximate the derivatives in Equation (2.3), and collecting like terms, yields a FD equation of the form (2.4). The coefficients $a_p,a_w,a_e,a_s,a_n$ and source term 𝓈_p are recorded in the Appendix, Equations (A.1a-A.1b). The derivatives at w, e, s and n appearing in the source terms are approximated by central FD formulae.

2.2.2. High-Order Schemes

The order of discretization is one of the criteria used to evaluate the solution accuracy in a FD method. Thus, higher-order FD methods have generated intense interest for many applications [13] [14]. Researchers have developed two main techniques for HO discretization [15]-[18]—a wider stencil which suffers from low-order accuracy at nodes adjacent to boundaries along with more complicated matrix calculations, and the compact Padé-Hermitian FDM, which employs narrower stencils. For example, 4^th-order Padé-Hermitian approximations for 1^st- and 2^nd-derivatives in the $\xi$ -direction are derived from Equation (2.5) by setting $c{2}_{nb}^{\left(1\right)}=0$ for the 1^st-derivative, and $c{1}_{nb}^{\left(2\right)}=0$ for the 2^nd-derivative, yielding the following 4^th-order approximations:

${\frac{\partial \varphi }{\partial \xi }|}_p\approx \frac{3}{4}\left({\varphi }_e-{\varphi }_w\right)-\frac{1}{4}\left({\frac{\partial \varphi }{\partial \xi }|}_e+{\frac{\partial \varphi }{\partial \xi }|}_w\right)$ (2.6)

${\frac{{\partial }^2\varphi }{\partial {\xi }^2}|}_p\approx \frac{6}{5}\left({\varphi }_w-2{\varphi }_p+{\varphi }_e\right)-\frac{1}{10}\left({\frac{{\partial }^2\varphi }{\partial {\xi }^2}|}_e+{\frac{{\partial }^2\varphi }{\partial {\xi }^2}|}_w\right)$ (2.7)

with similar expressions in the $\eta$ -direction. This combination of CCST-FDM and Padé-Hermitian scheme is referred to as CCST-HOM1 throughout this study.

An alternative compact 4^th-order approximation for the 2^nd-derivative of $\varphi$ (in the $\xi$ -direction), can be expressed as:

${\frac{{\partial }^2\varphi }{\partial {\xi }^2}|}_p\approx 2\left({\varphi }_w-2{\varphi }_p+{\varphi }_e\right)+\frac{1}{2}\left({\frac{\partial \varphi }{\partial \xi }|}_w-{\frac{\partial \varphi }{\partial \xi }|}_e\right)$ (2.8)

which is obtained by setting m = 2 and $c{2}_{nb}^{\left(2\right)}=0$ in Equation (2.5). The scheme which approximates the 2^nd-derivative using (2.8) is referred to as CCST-HOM2 throughout this study. The 4^th-order accurate approximation of the 1^st-derivative (for the convection term) in the CCST-HOM2 scheme is the same as in CCST-HOM1, i.e., Equation (2.6). These schemes have been applied to a Poisson equation with known solution, on domains with regular and irregular boundary nodes, confirming their global 4^th-order accuracy [11] [19]. The coefficients $a_p,a_w,a_e,a_s,a_n$ and source terms 𝓈_p for these 4^th-order schemes are recorded in the Appendix, Equations (A.2a-A.2c). The derivatives at w, e, s and n appearing in the source terms are approximated by central FD formulae. Both HO CCST formulations are constructed on a 5-point stencil, and no special treatment is required at nodes adjacent to the boundaries. Further details of these HO schemes, such as derivative approximations at endpoints w, e, s and n of the computational stencil, are provided in [11] [19].

2.3. Boundary Conditions for CCST-FDM

The treatment of boundary nodes in CCST-FDM is dependent on the type of condition imposed at the boundary node, which may be either a Dirichlet or Neumann condition, expressed in mixed form as $\delta \widehat{n}\cdot \nabla \overline{\varphi }+\gamma \overline{\varphi }=\overline{G}\left(x,y\right)$, where $\widehat{n}=\left(n_x,n_y\right)$ is the unit normal vector at the boundary, $\delta =0$ yields a Dirichlet condition and $\gamma =0$ corresponds to a Neumann condition. In $\xi ,\eta$ variables, this can be written as

$\delta \widehat{n}\cdot \left(\frac{1}{x^{\prime }}\frac{\partial \varphi }{\partial \xi },\frac{1}{y^{\prime }}\frac{\partial \varphi }{\partial \eta }\right)+\gamma \varphi =G\left(\xi ,\eta \right)$. (2.9)

In the case of a Dirichlet condition, the value at the boundary node is known so this value can be directly used in the calculation on the corresponding stencil adjacent to the boundary. For the Neumann boundary condition, as indicated in Figure 2, a one-sided differencing scheme with accuracy related to either 2^nd-order or HO discretization, can be used to approximate the 1^st-derivative at the boundary node.

For those cases where only one grid line, either horizontal or vertical, passes from the boundary node, the node is referred to as an “irregular” boundary node. Thus, the dashed line in Figure 2(a) is not a real grid line, and the subscripts $e{s}^{*}$ and $es{s}^{*}$ in Equation (2.10) refer to computational nodes corresponding to the fictitious physical nodes inserted at ${\text{ES}}^{\text{*}}$ and ${\text{ESS}}^{\text{*}}$ . Distances $L_1$ and $L_2$ , from ${\text{ES}}^{\text{*}}$ to $\text{S}$ and $\text{ES}$ respectively, are introduced to employ weighted averaging

for calculation of ${\varphi }_{e{s}^{*}}=\frac{L_2{\varphi }_s+L_1{\varphi }_{es}}{L_1+L_2}$ , and similarly for ${\varphi }_{es{s}^{*}}=\frac{L_2{\varphi }_{ss}+L_1{\varphi }_{ess}}{L_1+L_2}$ . The

case of a regular boundary node (see Figure 2(b)) is recovered by setting $L_2=0$ .

Research has shown that a k^th-order scheme for the interior nodes along with a (k − 1)-order scheme for the boundary nodes approximation can retain the global accuracy of the scheme used at the interior nodes [20] [21]. In this paper, for the HO CCST-FDM, a third-order accurate formulation is implemented at boundary nodes with a Neumann boundary condition. Thus, on the computational stencil, the mixed boundary condition is discretized as:

Figure 2. Schematic of the boundary node E with Neumann condition: (a) irregular boundary node, (b) regular boundary node.

$\begin{array}{l}\delta \left[{\left(\frac{n_x}{x^{\prime }}\right)}_e\left(\frac{5{\varphi }_e-4{\varphi }_p-{\varphi }_w}{2}-2{\frac{\partial \varphi }{\partial \xi }|}_p\right)+{\left(\frac{n_y}{y^{\prime }}\right)}_e\left(\frac{5{\varphi }_e-4{\varphi }_{e{s}^{*}}-{\varphi }_{es{s}^{*}}}{2}-2{\frac{\partial \varphi }{\partial \eta }|}_{e{s}^{*}}\right)\right]\\ +\gamma {\varphi }_e={G\left(\xi ,\eta \right)|}_e\end{array}$ (2.10)

Equation (2.10) can be rewritten to calculate the boundary value ${\varphi }_e$ explicitly. The first derivatives at points $p$ and $e{s}^{*}$ are calculated depending on the location of center point $p$ , and updated through the iterative solution process [11] [19].

3. CCST-FDM Formulation of Streamfunction-Vorticity Equations

The streamfunction Equation (1.2a) is a Poisson equation which, in $\xi ,\eta$ variables, is obtained from the general convection-diffusion Equation (2.3) by setting $\varphi =\psi$ , $U=0$ , $V=0$ , $\text{Γ}=1$ and $s=\omega$ . The vorticity equation is obtained from (2.3) by taking $\varphi =\omega$ , $U=u$ , $V=v$ , $\text{Γ}=R{e}^{-1}$ and $s=0$ .

3.1. Low-Order Discretization of ψ-ω Equations

Following the discussion in section 2.2.1 and using the above settings for the streamfunction equation, the coefficients and source term in the FD Equation (2.4) for the streamfunction, for the low-order hybrid schemes, are given by Equations (A.1a-A.1b). Similarly, the low-order coefficients in the FD equation for vorticity are obtained by inserting the appropriate quantities indicated above into Equations (A.1a-A.1b).

3.2. High-Order Discretizations of ψ-ω Equations

High-order discretization schemes for the streamfunction-vorticity formulation of the Navier-Stokes equations have been extensively discussed in the literature. Erturk [22] compared wide and compact 4^th-order schemes. High-order wide finite difference formulations use a nine-point stencil [23]-[25], but struggle with complex domains, requiring various alternative stencils based on the location of the stencil central point. Moreover, these stencils typically use transformation functions $x=x\left(\xi ,\eta \right)$ and $y=y\left(\xi ,\eta \right)$ for complex domains, complicating the procedure compared to the 1-D mappings in CCST-FDM.

Applying the settings mentioned above in the convection-diffusion Equation (2.3), the CCST-HOM1 discretized forms (2.4) of the streamfunction and vorticity Equations (1.2) have coefficients recorded in Equation (A.2a) with the value λ = 3/5 and the source term given by (A.2b). Alternatively, Equation (A.2a) with λ = 1 provides the coefficients for the CCST-HOM2 discretization, with source term given by (A.2c).

4. CCST-FDM Results for Lid-Driven Cavity Flows

Lid-driven cavity flows, a key fluid mechanics benchmark problem in CFD, will be used to assess the capabilities and results of the high-order CCST-FD methods, HOM1 and HOM2. Results for square, L-shape, skewed and triangular cavities are compared to published numerical results to verify the methods used in the CCST-FDM coding. A non-uniform grid is used to focus on areas with high gradients that occur near the boundaries of the cavity, supporting the CCST-FDM’s effectiveness in reducing computational time and resources. A logarithmic clustering function

$\delta \left(x\right)=0.5+0.5\left\{\frac{\ln \left|\frac{B-1+2x}{B+1-2x}\right|}{\ln \left|\frac{B+1}{B-1}\right|}\right\}$ (4.1)

with parameter B = 1.15 is used to achieve the desired clustering [11]. The fluid motion in the cavity is driven by the top wall, which slides in the x-direction at a constant speed $u=u_{lid}$ , while the other walls are held fixed.

4.1. Boundary Conditions for Streamfunction and Vorticity

The boundary conditions for the streamfunction are well-known, i.e., $\psi$ = constant on all walls and $\partial \psi /\partial y=u_{lid}$ on the sliding top wall. We also note that the no-slip condition implies that $\partial \psi /\partial y=u=0$ on the stationary walls and $\partial \psi /\partial x=-v=0$ on all walls.

Boundary conditions for vorticity, on the other hand, have been a topic of considerable discussion since the 1930s, when Thom [26] [27] applied the no-slip condition to propose the 1^st-order accurate approximation formula for the wall

vorticity ${\omega }_S=\frac{2\left({\psi }_{S+1}-{\psi }_S\right)}{{\left(\text{Δ}h\right)}^2}$ , where ${\psi }_S$ and ${\psi }_{S+1}$ refer to the streamfunction

value on the wall and at the first node at distance $\text{Δ}h$ normal to the wall, respectively. However, some studies argue that no physical boundary condition is defined for vorticity and consequently, no numerical boundary condition should be imposed on vorticity [3] [23] [28]. Despite this, most studies still employ various differencing expressions with different orders of accuracy for boundary vorticity calculations [29]-[32]. Orszag and Israeli [33] have provided a summary of high-order approximation methods for boundary vorticity.

In the present study, the vorticity on the cavity walls is computed using the compact finite difference method. Noting that Equation (1.2a) is also valid on the boundaries yields Dirichlet boundary conditions for vorticity. Referring to the discussion in section 2.3, we take $\varphi =\omega$ , $\delta =0$ , $\gamma =1$ and $G=-{\nabla }^2\psi$ , where ${\nabla }^2$ is the transformed Laplacian shown in Equation (2.3). $G$ is evaluated at each iteration, as new values of $\psi$ become available. As discussed above, the cut-stencil transformation facilitates use of non-uniform stencils and regular or irregular boundary nodes.

For example, if the west node W is a boundary node, then ${{\partial }^2\psi /\partial y^2|}_w=0$ since $\psi$ is constant on the cavity walls, and Equation (1.2a) reduces to

${\omega }_w=-{\frac{{\partial }^2\psi }{\partial x^2}|}_w=-{\left\{\frac{1}{{\left(x^{\prime }\right)}^2}\frac{{\partial }^2\psi }{\partial {\xi }^2}-\frac{x^{″}}{{\left(x^{\prime }\right)}^3}\frac{\partial \psi }{\partial \xi }\right\}|}_w={\left\{\frac{1}{{\left(x^{\prime }\right)}^2}\frac{{\partial }^2\psi }{\partial {\xi }^2}\right\}|}_w$ , since

$v_w={\left\{-\frac{1}{x^{\prime }}\frac{\partial \psi }{\partial \xi }\right\}|}_w=0$ . A similar analysis can be applied to obtain a boundary condition for vorticity on the lid (a north node N), noting that $u_n={\left\{\frac{1}{y^{\prime }}\frac{\partial \psi }{\partial \eta }\right\}|}_n=u_{lid}$ .

Using Equation (2.5) to express the second derivatives and expanding in Taylor series gives 2^nd-order approximations for ${\omega }_w$ and ${\omega }_{lid}$ :

${\omega }_w=-\left[\frac{-6{\psi }_w+6{\psi }_p}{{\left({x^{\prime }}_w\right)}^2}-\frac{2}{{\left({x^{\prime }}_w\right)}^2}\left(-{x^{\prime }}_p{v}_p\right)\right]$ (4.2a)

${\omega }_{lid}=-\left[\frac{-6{\psi }_{lid}+6{\psi }_p}{{\left({y^{\prime }}_{lid}\right)}^2}+\left(\frac{4}{{y^{\prime }}_{lid}}-\frac{{y^{″}}_{lid}}{{\left({y^{\prime }}_{lid}\right)}^2}\right)u_{lid}+\frac{2}{{\left({y^{\prime }}_{lid}\right)}^2}\left({y^{\prime }}_p{u}_p\right)\right]$ (4.2b)

4.2. Lid-Driven Flow in a Square Cavity

High-order CCST-FDM solutions for the lid-driven cavity flow in a square domain are considered in this section. Verification of the proposed high-order formulae, introduced in section 3.2, is associated with demonstrating the capability of these HO formulations to accurately capture the flow features and quantitative parameters with a coarser mesh compared to that used for 2^nd-order schemes.

Initial computations using high-order central differencing approximation of the convective terms failed to converge for grids of 65 × 65 and 81 × 81 nodes, even for $Re=400$ (based on lid speed and cavity length), for both CCST-HOM1 and CCST-HOM2. Many researchers have discussed the cause of this failure in convective-dominated flows, and have provided different remedies to resolve this issue [1] [6] [8] [12] [34] [35]. However, computations for these cases using HO upwind formulations show relatively good agreement on key qualitative and quantitative features with those in the literature. The computations start with evaluation of the transformation metrics for the grid, which remain fixed throughout the iterative solution procedure. An initialization of the solution variables is followed by a sweep of all nodes in the cavity, first for vorticity then for streamfunction. Velocity components in the vorticity equation are computed from derivatives of the streamfunction and updated at each iteration. Upwinding (backward or forward) or central differencing is determined by the sign of the velocity components at each node. The successive under-relaxation (SUR) method is used to perform the computations for the CCST-HOMs and the optimum value of the relaxation factor $\sigma$ can be determined for each grid size and each HO formulation. The iterations were considered to have converged when the absolute difference between two successive iterations fell below 1 × 10⁻¹². For the 65 × 65 grid, the optimum value of $\sigma$ was found to be 0.30 for both CCST-HOMs. For the 81 × 81 grid, the optimum value of $\sigma$ was 0.35 for HOM1 and 0.465 for HOM2. Further details can be found in [11]. Employing a HO upwind formulation at $Re=1000$ , on an 81 × 81 grid and CCST-HOM1, the absolute value of streamfunction ($\psi$ ) at the center of the primary vortex, the bottom left and bottom right vortices are 1.21E−1, 2.27E−4 and 1.87E−3, respectively. The corresponding values using CCST-HOM2 on the same grid are 1.22E−1, and 2.41E−4 and 1.87E−3.

Comparison of cut-stencil solutions in Table 1 with studies by Gupta and Kalita [24], Pandit [36] and Nishida and Satofuka [37] indicate close agreement on nearly all key features of the solution, except for the value of streamfunction at the bottom right vortex predicted by CCST-HOM1 with a non-uniform 65 × 65 grid. This issue has been resolved by the CCST-HOM2 solution with the same grid, or by the CCST-HOM1 formulation with a finer 81 × 81grid.

Table 1. CCST-HOMs vs. other HO solutions for square cavity flow ($Re=1000$ ).

[Study] (Grid)	Primary vortex			B.L.* vortex			B.R.** vortex
	Abs. $\left(\psi \right)$	Location $\left(x,y\right)$	Abs. $\left(\omega \right)$	Abs. $\left(\psi \right)$	Location $\left(x,y\right)$	Abs. $\left(\omega \right)$	Abs. $\left(\psi \right)$	Location $\left(x,y\right)$	Abs. $\left(\omega \right)$
[24] (81 × 81)	0.117	(0.5250, 0.5625)	N. A.	2.02E−4	(0.0875, 0.0750)	N.A.	1.70E−3	(0.8625, 0.1125)	N.A.
[36] (61 × 61)	0.118	(0.5266, (0.5532)	N. A.	2.31E−4	(0.0840, 0.0840)	N.A.	1.72E−3	(0.8577, 0.1092)	N.A.
[37] (65 × 65)	0.118	(0.5313, 0.5625)	2.05692	2.24E−4	(0.0781, 0.0781)	3.134E−1	1.67E−3	(0.8594, 0.1094)	1.125
[37] (129 × 129)	0.119	(0.5313, 0.5625)	2.06616	2.32E−4	(0.0859, 0.0781)	3.671E−1	1.72E−3	(0.8594, 0.1094)	1.144
CCST-HOM1 (65 × 65)	0.120	(0.5239, 0.5714)	2.08386	2.30E−4	(0.0839, 0.0728)	3.209E−1	1.81E−3	(0.8648, 0.1213)	1.175
CCST-HOM1 (81 × 81)	0.119	(0.5382, 0.5572)	2.07696	2.35E−4	(0.0794, 0.0794)	3.432E−1	1.78E−3	(0.8590, 0.1081)	1.121
CCST-HOM2 (65 × 65)	0.118	(0.5239, 0.5714)	2.06372	2.30E−4	(0.0839, 0.0728)	3.161E−1	1.76E−3	(0.8648, 0.1081)	1.031
CCST-HOM2 (81 × 81)	0.119	(0.5382, 0.5572)	2.06270	2.33E−4	(0.0794, 0.0794)	3.415E−1	1.74E−3	(0.8590, 0.1081)	1.127

*: Bottom left, **: Bottom right.

Analysis of the numerical results of Gupta and Kalita [24] reveals that their high-order solution was not able to capture the streamfunction value at the bottom left corner as accurately as other studies. Furthermore, the streamfunction values at both corner vortices reported in [37] with a 65 × 65 grid are somewhat different from other studies, such as [2] [36] [38]. The streamfunction and vorticity contours for the CCST-HOM2 solution with a non-uniform 65 × 65 grid are shown in Figure 3, providing good qualitative comparison between the results from other numerical studies and the 4^th-order solution of the present study.

Figure 3. Contours of CCST-HOM2 solution for lid-driven square cavity flow ($Re=1000$ , non-uniform 65 × 65 grid, 4225 active nodes): (a) streamfunction, (b) vorticity.

4.3. Lid-Driven Flow in an L-Shape Cavity

The lid-driven cavity flow in an L-shape domain, as depicted in Figure 4(a), has been addressed in other numerical studies to verify domain decomposition methods [39] [40] and discretization techniques [39]. The same transformation functions and mapping procedure that have been used above for square domains can be employed for this semi-irregular domain, demonstrating the generality of the CCST-FD procedure to solve PDEs in any type of domain without requiring any changes in the mapping procedure. In the CCST-FD approach, any non-rectangular physical flow region is enclosed in a rectangle. Nodes within the flow region are identified as active nodes. The portion of the bounding rectangle outside of the flow region is labeled as inactive and shown as a blank region in the following contour plots.

Comparisons of the 2^nd-order CCST-FD solutions and other published numerical results show good agreement for the key features of the solution, where previous numerical studies have used much finer meshes than those in the present work. The sum of relative differences for extremes of $\psi$ in the present study on a 129 × 129 grid with 12545 active nodes, compared to the values in [39] on a 256 × 256 grid, are 12.02%, 2.71% and 2.08% for 2^nd-order, CCST-HOM1 and CCST-HOM2, respectively. These differences are reduced to 6.61%, 1.75% and 0.79%, respectively, when a 161 × 161 grid (19,521 active nodes) is employed. The contours of streamfunction for the CCST-HOM2 solution for the lid-driven cavity flow in an L-shape domain, using the 129 × 129 grid (12,545 active nodes), are plotted in Figure 4(b).

Figure 4. L-shape domain: a) schematic, b) streamfunction contours of CCST-HOM2 solution ($Re=1000$ , 129 × 129 grid, 12,545 active nodes).

4.4. Lid-Driven Flow in Skewed Quadrilateral Cavities

Typically, irregular domains contain non-uniform cut stencils which are created by intersection of only one grid line (either horizontal or vertical) and the boundary of the domain of interest. The domains in this study are skewed quadrilateral and triangular cavities found in other studies, which can be used to verify solutions obtained by the CCST-FDMs.

Two skew angles, $\alpha =45˚$ and $\alpha =135˚$ , as depicted in Figure 5(a), Figure 5(b), are considered. Oosterlee et al. [39] proposed a FVM solution using curvilinear coordinates, while Erturk and Dursun [41] addressed skewed cavity flow with the FDM and a non-orthogonal skewed grid system. Their grid systems required the use of 2-D transformation functions, leading to more complexity in the equations compared to the CCST-FDM setup.

Figure 5. Schematic of domain for skewed lid-driven cavities: (a) $\alpha =45˚$ , (b) $\alpha =135˚$ , and isosceles right triangular lid-driven cavities: (c) left-side aligned, (d) right-side aligned.

Although not shown here, the data for $Re=100$ indicates that all cut-stencil FD results align well with other studies, but CCST-FDM can predict solutions of comparable accuracy on a coarser mesh [11]. No significant differences are seen between the 2^nd-order and HO solutions or between the grid sizes used in the CCST-FD solution.

At higher Reynolds numbers, point-Jacobi iterations converge if the convective derivatives in the vorticity equation are upwinded, though this reduces the solution accuracy. More accurate results at this Reynolds number can be obtained using central differencing for the convective terms, with convergence achieved through successive under-relaxation (SUR). The optimum SUR factor is found by adjusting $\sigma$ within a specific range. The results of using central schemes are reported in Table 2, indicating that better agreement between the present CCST-FD solutions and other numerical results is achieved when central differencing is used for the convective derivatives for each formulation. The capability of the high-order formulations to produce values closer to those from other numerical studies is more apparent for this higher value of Reynolds number ($Re=1000$ , $\alpha =45˚$ ) and, between the two high-order CCST-FD formulations, CCST-HOM2 appears to yield slightly better agreement with other published data. Additionally, the point-Jacobi method converges for CCST-HOM2 on a finer 273 × 113 mesh (18,193 active nodes) with central differencing of the convective terms. The contours of streamfunction obtained from CCST-HOM2, as depicted in Figure 6, shows good qualitative agreement with those from other numerical studies, such as [41].

Table 2. Comparison of CCST-FDM solutions for skewed lid-driven cavity flow ($Re=1000$ , $\alpha =45˚$ ).

Study (Grid)	Properties of vortex centre-1			Properties of vortex centre-2
	Abs. $\left(\psi \right)$	Location $\left(x,y\right)$	Abs. $\left(\omega \right)$	Abs. $\left(\psi \right)$	Location $\left(x,y\right)$	Abs. $\left(\omega \right)$
[39] (256 × 256)	5.3523E−2	(1.3128, 0.5745)	N.A.	1.0039E−2	(0.7775, 0.4005)	N.A.
[41] (513 × 513)	5.3423E−2	(1.3148, 0.5745)	6.955	1.0024E−2	(0.7780, 0.3991)	6.269E−1
[43] (320 × 320)	5.2553E−2	(1.3120, 0.5745)	N.A.	1.0039E−2	(0.7766, 0.3985)	N.A.
[44] (120 × 120)	5.4690E−2	(1.3100, 0.5700)	N.A.	1.0170E−2	(0.7760, 0.3980)	N.A.
2nd-order cut-stencil (137 × 57 − 4617*)	4.6779E−2 $\left(0.400\le \sigma \le 0.660\right)$ i $\left({\sigma }_{\text{Opt}.}=0.565\right)$ ii	(1.3500, 0.5875)	7.961	9.45E−3	(0.7875, 0.3875)	5.756E−1
Cut-stencil CCST-HOM1 (137 × 57 − 4617*)	5.6464E−2 $\left(0.450\le \sigma \le 0.550\right)$ $\left({\sigma }_{\text{Opt}.}=0.512\right)$	(1.2875, 0.5625)	6.510	1.01E−2	(0.7750, 0.4000)	6.454E−1
Cut-stencil CCST-HOM2 (137 × 57 − 4617*)	5.3191E−2 $\left(0.470\le \sigma \le 0.656\right)$ $\left({\sigma }_{\text{Opt}.}=0.656\right)$	(1.3125, 0.5625)	6.992	1.01E−2	(0.7750, 0.4000)	6.404E−1
2nd-order cut-stencil (273 × 113 – 18,193*)	5.1531E−2 $\left(0.500\le \sigma \le 0.890\right)$ $\left({\sigma }_{\text{Opt}.}=0.890\right)$	(1.3186, 0.5750)	7.214	9.87E−3	(0.7750, 0.3938)	6.216E−1
Cut-stencil CCST-HOM1 (273 × 113 – 18,193*)	5.4278E−2 $\left(0.650\le \sigma \le 0.840\right)$ $\left({\sigma }_{\text{Opt}.}=0.830\right)$	(1.3000, 0.5625)	6.849	1.00E−2	(0.7688, 0.3938)	6.389E−1
Cut-stencil CCST-HOM2 (273 × 113 – 18,193*)	5.3438E−2	(1.3063, 0.5688)	6.968	1.00E−2	(0.7750, 0.3938)	6.270E−1

*: Number of active nodes. (i): Study range of under-relaxation factor, (ii): Optimum value of under-relaxation factor.

Figure 6. Streamfunction contours of CCST-HOM2 formulation for skewed lid-driven cavity ($Re=1000$ , $\alpha =45˚$ , 273 × 113 grid, 18,193 active nodes).

Contours of the streamfunction predicted by CCST-HOM2 for $Re=1000$ and $\alpha =135˚$ , using the finer mesh of 18193 active nodes, is plotted in Figure 7. For this irregular domain, the absolute values of streamfunction at the center of each main vortex are 9.4407E−2 and 4.30E−3 for the 2^nd-order accurate formulation when a grid of 4617 active nodes is used. If CCST-HOM2 is employed with the same grid size, these values are 9.4036E−2 and 3.69E−3. Quantitative and qualitative comparison with the contours in Erturk and Dursun [41] show more evidence of the accuracy of the methods and algorithms developed in this research, particularly noting that the data in [41] has been obtained with a much finer mesh (~260,000 nodes) than those used in the present study.

4.5. Lid-Driven Flow in Triangular Cavities

The CCST-FD solutions for lid-driven cavity flow in the domains shown in Figure 5(c), Figure 5(d) are considered in this section, to further demonstrate the capability of the CCST-FDM to produce accurate solutions of fluid flow problems in any irregular shaped domains.

Paramane and Sharma [35] used FVM to simulation lid-driven cavity flows for different cross-sections, such as skewed, trapezoidal or triangular cross-sections. Erturk and Gokcol [42] obtained solutions for general triangular cavities using FDM. They defined a transformation which maps an arbitrary triangle in the physical domain to an isosceles right triangle as the computational domain. The

Figure 7. Streamfunction contours of CCST-HOM2 formulation for skewed lid-driven cavity ($Re=1000$ , $\alpha =135˚$ , 273 × 113 grid, 18,193 active nodes).

mapping function used in [42] is solely applicable to the transformation of triangular domains. However, different mapping functions can be defined for different cross-sectional geometries, e.g. [41] [42] [45]. Additionally, the 2-D transformation functions in [42] produce cross derivatives in the mapped form of the streamfunction-vorticity equations, making them more complex and hence more difficult to discretize. It is important to note that in the CCST-FDM approach, the mapping process is performed with the same transformation functions for every geometrical cross-section, also generating the mapped form of the governing equations in a relatively simpler form. Flow inside the triangular cavity shown in Figure 5(c) has also been studied using other numerical schemes such as FEM [46] and the Lattice Boltzmann method [47].

For the domain depicted in Figure 5(c), at Reynolds number $Re=500$ , comparison of the results from CCST-HOM2 and Erturk and Gokcol [42] shows that the relative differences of $\psi$ and $\omega$ at the primary vortex are equal to 0.91% and 1.52%, respectively, for the 65 × 65 grid with 2145 active nodes in the current simulation, while a mesh of (512*512)/2 nodes is used in [42]. These differences are 3.41% and 4.34% for the 2^nd-order CCST-FD formulation on a 101 × 101 grid with 5151 active nodes. This provides further strong evidence that the CCST-FDMs, especially the HO methods, can accurately simulate complex flow phenomena on a coarse mesh in irregular domains. Additionally, comparisons between the CCST-FDMs and other numerical studies in this domain and for the same $Re$ , such as [46] [47], show the same level of agreement, despite these studies employing much finer grids.

Comparison of results for the CCST-FDM and other numerical studies for cavity flow inside the triangular cross-section shown in Figure 5(c), for $Re=1000$ , is presented in Table 3. This data demonstrates that there is good agreement between CCST-HOM1 and other numerical solutions for all grid sizes studied, with even better agreement when CCST-HOM2 is employed. The relative difference of $\psi$ and $\omega$ at the primary vortex between the CCST-HOM2 results and reported data in [42] can be seen as 0.28% and 0.41%, respectively, for the coarser mesh studied in Table 3. These values increase to 5.50% and 5.23% for $\psi$ and $\omega$ at the primary vortex, respectively, when the 2^nd-order CCST-FD formulation is employed on the finer grid (126 × 126 grid, 8001 active nodes). The data in Table 3 shows that more accurate solutions in irregular shaped domains can be obtained by applying the high-order formulations. Additionally, the benefit of central differencing of the convective terms can be achieved by employing SUR as the iterative method when necessary. This supports the proposition that the CCST-FDM is a powerful numerical formulation with known features of traditional well-established numerical methods and, even though it is a pure finite difference method, it can be used to solve complex PDEs in irregular shaped domains.

Table 3. Comparison of solutions for left-side aligned triangular cavity flow ($Re=1000$ ).

Study (Grid)	Properties of primary vortex center
	Abs. $\left(\psi \right)$	Location $\left(x,y\right)$	Abs. $\left(\omega \right)$
[42] (512 × 512)/2	0.05306	(0.6094, 0.8691)	7.022
[46] (100 subspaces)	0.05325	(0.6081, 0.8678)	6.997
[47] (300 × 300)	N.A.	(0.6050, 0.8650)	N.A.
2nd-order CCST-FDM (81 × 81 − 3321*)	0.04647 $\left(0.350\le \sigma \le 0.575\right)$ i $\left({\sigma }_{\text{Opt}.}=0.575\right)$ ii	(0.6500, 0.8875)	8.036
CCST-HOM1 (81 × 81 − 3321*)	0.05614 $\left(0.325\le \sigma \le 0.515\right)$ $\left({\sigma }_{\text{Opt}.}=0.515\right)$	(0.5875, 0.8625)	6.565
CCST-HOM2 (81 × 81 − 3321*)	0.05291 $\left(0.470\le \sigma \le 0.665\right)$ $\left({\sigma }_{\text{Opt}.}=0.665\right)$	(0.6125, 0.8625)	7.051
2nd-order CCST-FDM (101 × 101 − 5151*)	0.04865 $\left(0.425\le \sigma \le 0.665\right)$ $\left({\sigma }_{\text{Opt}.}=0.665\right)$	(0.6400, 0.8800)	7.669
CCST-HOM1 (101 × 101 − 5151*)	0.05510 $\left(0.400\le \sigma \le 0.600\right)$ $\left({\sigma }_{\text{Opt}.}=0.600\right)$	(0.6000, 0.8600)	6.703
CCST-HOM2 (101 × 101 − 5151*)	0.05305 $\left(0.400\le \sigma \le 0.700\right)$ $\left({\sigma }_{\text{Opt}.}=0.700\right)$	(0.6100, 0.8700)	7.032
2nd-order CCST-FDM (126 × 126 − 8001*)	0.05014 $\left(0.550\le \sigma \le 0.790\right)$ $\left({\sigma }_{\text{Opt}.}=0.790\right)$	(0.6320, 0.8800)	7.389
CCST-HOM1 (126 × 126 − 8001*)	0.05449 $\left(0.475\le \sigma \le 0.700\right)$ $\left({\sigma }_{\text{Opt}.}=0.700\right)$	(0.6000, 8640)	6.813
CCST-HOM2 (126 × 126 − 8001*)	0.05310 $\left(0.600\le \sigma \le 0.725\right)$ $\left({\sigma }_{\text{Opt}.}=0.725\right)$	(0.6080, 8720)	7.015

*: Number of active nodes. (i): Study range of under-relaxation factor, (ii): Optimum value of under-relaxation factor.

Figure 8 shows streamfunction contours for $Re=500$ and $Re=1000$ , for lid-driven cavity flow in a left-side aligned right triangle, employing CCST-HOM2 and using the finest grid size studied for each Reynolds number. The contours of streamfunction for $Re=500$ and $Re=1000$ , for lid-driven cavity flow in a right-side aligned right triangle employing CCST-HOM2 and using the finest grid size studied for each Reynolds number, are illustrated in Figure 9.

(a) (b)

Figure 8. Streamfunction contours for CCST-HOM2 formulation for left-side aligned right triangle lid-driven cavity flow: (a) $Re=500$ , 101 × 101 grid, 5151 active nodes, (b) $Re=1000$ , 126 × 126 grid, 8001 active nodes.

(a) (b)

Figure 9. Streamfunction contours for CCST-HOM2 formulation for right-side aligned right triangle lid-driven cavity flow: (a) $Re=500$ , 101 × 101 grid, 5151 active nodes, (b) $Re=1000$ , 126 × 126 grid, 8001 active nodes.

A detailed comparison of the above results can be found in [11]. For example, for lid-driven cavity flow inside the right-side aligned triangle (Figure 5(d)), for $Re=1000$ , the relative differences for streamfunction and vorticity for the 2^nd-order CCST-FD solution compared to results in [42] are 2.68% and 2.30%, respectively, for the 126 × 126 grid with 8001 active nodes. For the same grid, CCST-HOM1 reduces these differences to 0.13% and 0.26% for streamfunction and vorticity, respectively. This comparison shows relative differences reduced to 0.03% and 0.12% for streamfunction and vorticity, respectively, for the CCST-HOM2 solution.

5. Conclusions

The simulation of lid-driven cavity flow, as a fluid flow benchmark problem, has been explored using the Cartesian Cut-Stencil Finite Difference Method (CCST-FDM). The CCST-FDM is simple to formulate and is capable of solving PDEs in any irregular and complex shaped domain. This capability of CCST-FDM derives from the unique mapping of each finite difference stencil in the physical domain to a generic uniform computational stencil. In this research, the CCST-FDM, which exhibits 2^nd-order accuracy in its basic formulation, has been combined with Padé-Hermitian formulae to construct compact high-order accurate schemes.

Solutions based on the 2^nd-order accurate formulation of the CCST-FDM were initially compared to other available numerical results, confirming that the cut-stencil formulations can be applied to coupled sets of non-linear PDEs to solve complex fluid flows such as those described by the streamfunction-vorticity equations. Good agreement with relatively fine mesh results from the literature was obtained using the 2^nd-order accurate cut-stencil formulation with a coarser mesh. Two high-order formulations of the CCST-FDM were tested, initially in a regular square domain, and results were compared to high-order results in other numerical studies. Solutions using both CCST-HOM1 and CCST-HOM2 show good agreement with data from other researchers.

The CCST-FDM solutions of the streamfunction-vorticity equations for cavity flow in several irregular shaped domains were also compared to results in the literature. Unlike other finite difference studies for irregular shaped domains, the CCST-FDM does not require specialized grid generation for each different geometry of the domain. For all the domains, the quadratic transformation functions introduced in Section 2 provide the uniform computational stencil on which the numerical solution algorithm is based. The high-order CCST-FDM formulations developed in this study show excellent agreement with comparative studies, while the 2^nd-order formulation shows good agreement, even with relatively coarser meshes.

The point-Jacobi method was chosen as the default iterative scheme and, for cases when the flow is dominated by convection, central differencing of the convective terms led to divergence of the iterations. As expected, upwinding schemes were able to resolve this issue and provide a converged solution for convective-dominant flow. Additionally, successive under-relaxation gave a converged solution for central differencing of the convective terms, yielding more accurate results compared to the lower order upwind approximation.

Acknowledgements

The authors gratefully acknowledge funding to Professor Barron from the Discovery Grant program of the Natural Sciences and Engineering Research Council of Canada (NSERC).

Appendix: Scheme Coefficients and Source Terms

Discretization of the transformed convection-diffusion Equation (2.3) yields a finite difference equation of the form

$a_p{\varphi }_p+a_w{\varphi }_w+a_e{\varphi }_e+a_s{\varphi }_s+a_n{\varphi }_n=$ 𝓈_p (2.4)

for which different discretization schemes produce different coefficients and source terms.

A. LOW-ORDER HYBRID SCHEMES:

a) Convection: k^th-order upwind (k = 1, 2)

Backward: $\alpha ,\beta =+1$ , Forward: $\alpha ,\beta =-1$

2^nd-order central: $\alpha =0,\beta =0$

b) Diffusion: 2^nd-order central

c) Evaluation at stencil endpoints (w, e, s, n in source term): 2^nd-order central

d) Coefficients in (2.4):

$a_p=2\Gamma \left[\frac{1}{{\left({x^{\prime }}_p\right)}^2}+\frac{1}{{\left({y^{\prime }}_p\right)}^2}\right]+\frac{k\alpha U}{{x^{\prime }}_p}+\frac{k\beta V}{{y^{\prime }}_p}$

$a_w=-\Gamma \left[\frac{1}{{\left({x^{\prime }}_p\right)}^2}+\frac{{x^{″}}_p}{2{\left({x^{\prime }}_p\right)}^3}\right]-{\left(1+\alpha \right)}^k\frac{U}{2{x^{\prime }}_p}$

$a_e=-\Gamma \left[\frac{1}{{\left({x^{\prime }}_p\right)}^2}-\frac{{x^{″}}_p}{2{\left({x^{\prime }}_p\right)}^3}\right]+{\left(1-\alpha \right)}^k\frac{U}{2{x^{\prime }}_p}$ (A.1a)

$a_s=-\Gamma \left[\frac{1}{{\left({y^{\prime }}_p\right)}^2}+\frac{{y^{″}}_p}{2{\left({y^{\prime }}_p\right)}^3}\right]-{\left(1+\beta \right)}^k\frac{V}{2{y^{\prime }}_p}$

$a_n=-\Gamma \left[\frac{1}{{\left({y^{\prime }}_p\right)}^2}-\frac{{y^{″}}_p}{2{\left({y^{\prime }}_p\right)}^3}\right]+{\left(1-\beta \right)}^k\frac{V}{2{y^{\prime }}_p}$

e) Source term in (2.4):

For k = 1:

𝓈_p$=s_p$

For k = 2:

𝓈_p$=s_p+\frac{\alpha U}{2{x^{\prime }}_p}\left[\left(1+\alpha \right){\frac{\partial \varphi }{\partial \xi }|}_w-\left(1-\alpha \right){\frac{\partial \varphi }{\partial \xi }|}_e\right]+\frac{\beta V}{2{y^{\prime }}_p}\left[\left(1+\beta \right){\frac{\partial \varphi }{\partial \eta }|}_s-\left(1-\beta \right){\frac{\partial \varphi }{\partial \eta }|}_n\right]$ (A.1b)

B. HIGH-ORDER SCHEMES

1. HOM1 (4^th-order)

a) Convection: Equation (2.6)

b) Diffusion: Equation (2.6) for first derivatives, Equation (2.7) for second derivatives.

c) Evaluation at stencil endpoints (w, e, s, n in source term): 2^nd-order central

d) Coefficients in (2.4), with $\lambda =3/5$ :

$a_p=4\lambda \Gamma \left[\frac{1}{{\left({x^{\prime }}_p\right)}^2}+\frac{1}{{\left({y^{\prime }}_p\right)}^2}\right]$

$a_w=-\Gamma \left[\frac{2\lambda }{{\left({x^{\prime }}_p\right)}^2}+\frac{3{x^{″}}_p}{4{\left({x^{\prime }}_p\right)}^3}\right]-\frac{3U}{4{x^{\prime }}_p}$

$a_e=-\Gamma \left[\frac{2\lambda }{{\left({x^{\prime }}_p\right)}^2}-\frac{3{x^{″}}_p}{4{\left({x^{\prime }}_p\right)}^3}\right]+\frac{3U}{4{x^{\prime }}_p}$ (A.2a)

$a_s=-\Gamma \left[\frac{2\lambda }{{\left({y^{\prime }}_p\right)}^2}+\frac{3{y^{″}}_p}{4{\left({y^{\prime }}_p\right)}^3}\right]-\frac{3V}{4{y^{\prime }}_p}$

$a_n=-\Gamma \left[\frac{2\lambda }{{\left({y^{\prime }}_p\right)}^2}-\frac{3{y^{″}}_p}{4{\left({y^{\prime }}_p\right)}^3}\right]+\frac{3V}{4{y^{\prime }}_p}$

e) Source term in (2.4):

𝓈_p$=s_p+\left[\frac{U}{4{x^{\prime }}_p}+\frac{\Gamma {x^{″}}_p}{4{\left({x^{\prime }}_p\right)}^3}\right]\left[{\frac{\partial \varphi }{\partial \xi }|}_w+{\frac{\partial \varphi }{\partial \xi }|}_e\right]+\left[\frac{V}{4{y^{\prime }}_p}+\frac{\Gamma {y^{″}}_p}{4{\left({y^{\prime }}_p\right)}^3}\right]\left[{\frac{\partial \varphi }{\partial \eta }|}_s+{\frac{\partial \varphi }{\partial \eta }|}_n\right]$

$-\frac{\Gamma }{10{\left({x^{\prime }}_p\right)}^2}\left[{\frac{{\partial }^2\varphi }{\partial {\xi }^2}|}_w+{\frac{{\partial }^2\varphi }{\partial {\xi }^2}|}_e\right]-\frac{\Gamma }{10{\left({y^{\prime }}_p\right)}^2}\left[{\frac{{\partial }^2\varphi }{\partial {\eta }^2}|}_s+{\frac{{\partial }^2\varphi }{\partial {\eta }^2}|}_n\right]$ (A.2b)

2. HOM2 (4^th-order)

a) Convection: Equation (2.6)

b) Diffusion: Equation (2.6) for first derivatives, Equation (2.8) for second derivatives.

c) Evaluation at stencil endpoints (w, e, s, n in source term): 2^nd-order central

d) Coefficients in (2.4): Eqns. (A.2a), with $\lambda =1$ .

e) Source term in (2.4):

𝓈_p$=s_p+\frac{\Gamma }{2{\left({x^{\prime }}_p\right)}^2}\left[{\frac{\partial \varphi }{\partial \xi }|}_w-{\frac{\partial \varphi }{\partial \xi }|}_e\right]+\left[\frac{U}{4{x^{\prime }}_p}+\frac{\Gamma {x^{″}}_p}{4{\left({x^{\prime }}_p\right)}^3}\right]\left[{\frac{\partial \varphi }{\partial \xi }|}_w+{\frac{\partial \varphi }{\partial \xi }|}_e\right]$

$+\frac{\Gamma }{2{\left({y^{\prime }}_p\right)}^2}\left[{\frac{\partial \varphi }{\partial \eta }|}_s-{\frac{\partial \varphi }{\partial \eta }|}_n\right]+\left[\frac{V}{4{y^{\prime }}_p}+\frac{\Gamma {y^{″}}_p}{4{\left({y^{\prime }}_p\right)}^3}\right]\left[{\frac{\partial \varphi }{\partial \eta }|}_s+{\frac{\partial \varphi }{\partial \eta }|}_n\right]$ (A.2c)

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Hoffmann, K.A. and Chiang, S.T. (2004) Computational Fluid Dynamics. Vol. 1. 4th Edition, Engineering Education System, Wichita, Kansas, USA.
[2]	Ghia, U., Ghia, K.N. and Shin, C.T. (1982) High Resolutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method. Journal of Computational Physics, 48, 387-411.[CrossRef]
[3]	Spotz, W.F. (1998) Accuracy and Performance of Numerical Wall Boundary Conditions for Steady, 2D, Incompressible Streamfunction Vorticity. International Journal for Numerical Methods in Fluids, 28, 737-757.[CrossRef]
[4]	Yu, P.X., Tian, Z.F. and Zhang, H. (2017) A Rational High-Order Compact Difference Method for the Steady-State Stream Function-Vorticity Formulation of the Navier–Stokes Equations. Computers & Mathematics with Applications, 73, 1461-1484.[CrossRef]
[5]	Yu, Q., Xu, H., Liao, S. and Yang, Z. (2019) A Novel Homotopy-Wavelet Approach for Solving Stream Function-Vorticity Formulation of Navier-Stokes Equations. Communications in Nonlinear Science and Numerical Simulation, 67, 124-151.[CrossRef]
[6]	Ferziger, J.H. (1981) Numerical Methods for Engineering Application. Wiley.
[7]	Patankar, S.V. (1980) Numerical Heat Transfer and Fluid Flow. Series in Computational Methods in Mechanics and Thermal Sciences, Hemisphere Publishing.
[8]	Versteeg, H.K. and Malalasekera, W. (2007) An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Pearson Education Ltd.
[9]	Reddy, J.N. (1984) An Introduction to the Finite Element Method. McGraw-Hill.
[10]	Esmaeilzadeh, M., Barron, R.M. and Balachandar, R. (2020) Numerical Solution of Partial Differential Equations in Arbitrary Shaped Domains Using Cartesian Cut-Stencil Finite Difference Method. Part I: Concepts and Fundamentals. Numerical Mathematics: Theory, Methods and Applications, 13, 881-907.[CrossRef]
[11]	Esmaeilzadeh, M. (2016) A Cartesian Cut-Stencil Method for the Finite Difference Solution of PDEs in Complex Domains. Ph.D. Dissertation, University of Windsor.
[12]	Choo, S.M. and Chung, S.K. (2000) High-Order Perturbation-Difference Scheme for a Convection-Diffusion Problem. Computer Methods in Applied Mechanics and Engineering, 190, 721-732.[CrossRef]
[13]	Liu, D., Kuang, W. and Tangborn, A. (2009) High-Order Compact Implicit Difference Methods for Parabolic Equations in Geodynamo Simulation. Advances in Mathematical Physics, 2009, Article 568296.[CrossRef]
[14]	Zhang, W. and Jiang, J. (2017) A New Family of Fourth-Order Locally One-Dimensional Schemes for the Three-Dimensional Wave Equation. Journal of Computational and Applied Mathematics, 311, 130-147.[CrossRef]
[15]	Adam, Y. (1975) A Hermitian Finite Difference Method for the Solution of Parabolic Equations. Computers & Mathematics with Applications, 1, 393-406.[CrossRef]
[16]	Chu, P.C. and Fan, C. (2000) A Three-Point Sixth-Order Staggered Combined Compact Difference Scheme. Mathematical and Computer Modelling, 32, 323-340.[CrossRef]
[17]	Kong, L., Zhu, P., Wang, Y. and Zeng, Z. (2019) Efficient and Accurate Numerical Methods for the Multidimensional Convection-Diffusion Equations. Mathematics and Computers in Simulation, 162, 179-194.[CrossRef]
[18]	Patel, K.S. and Mehra, M. (2020) Fourth Order Compact Scheme for Space Fractional Advection-Diffusion Reaction Equations with Variable Coefficients. Journal of Computational and Applied Mathematics, 380, Article 112963.[CrossRef]
[19]	Esmaeilzadeh, M. and Barron, R.M. (2022) Numerical Solution of Partial Differential Equations in Arbitrary Shaped Domains Using Cartesian Cut-Stencil Finite Difference Method. Part II: Higher-Order Schemes. Numerical Mathematics: Theory, Methods and Applications, 15, 819-850.[CrossRef]
[20]	Zhong, X. (1998) High-Order Finite-Difference Schemes for Numerical Simulation of Hypersonic Boundary-Layer Transition. Journal of Computational Physics, 144, 662-709.[CrossRef]
[21]	Zhong, X. and Tatineni, M. (2003) High-Order Non-Uniform Grid Schemes for Numerical Simulation of Hypersonic Boundary-Layer Stability and Transition. Journal of Computational Physics, 190, 419-458.[CrossRef]
[22]	Erturk, E. (2009) Comparison of Wide and Compact Fourth-Order Formulations of the Navier-Stokes Equations. International Journal for Numerical Methods in Fluids, 60, 992-1010.[CrossRef]
[23]	Li, M., Tang, T. and Fornberg, B. (1995) A Compact Fourth-Order Finite Difference Scheme for the Steady Incompressible Navier-Stokes Equations. International Journal for Numerical Methods in Fluids, 20, 1137-1151.[CrossRef]
[24]	Gupta, M.M. and Kalita, J.C. (2005) A New Paradigm for Solving Navier-Stokes Equations: Streamfunction-Velocity Formulation. Journal of Computational Physics, 207, 52-68.[CrossRef]
[25]	Pandit, S.K., Kalita, J.C. and Dalal, D.C. (2008) A Fourth-Order Accurate Compact Scheme for the Solution of Steady Navier-Stokes Equations on Non-Uniform Grids. Computers & Fluids, 37, 121-134.[CrossRef]
[26]	Thom, A. (1928) An Investigation of Fluid Flow in Two Dimensions. Reports and Memoranda, No. 1194. Aeronautical Research Committee.
[27]	Thom, A. (1933) The Flow Past Circular Cylinders at Low Speeds. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 141, 651-669.[CrossRef]
[28]	Gresho, P.M. (1991) Incompressible Fluid Dynamics: Some Fundamental Formulation Issues. Annual Review of Fluid Mechanics, 23, 413-453.[CrossRef]
[29]	Sahin, M. and Owens, R.G. (2003) A Novel Fully Implicit Finite Volume Method Applied to the Lid-Driven Cavity Problem—Part I: High Reynolds Number Flow Calculations. International Journal for Numerical Methods in Fluids, 42, 57-77.[CrossRef]
[30]	Benjamin, A.S. and Denny, V.E. (1979) On the Convergence of Numerical Solutions for 2-D Flows in a Cavity at Large Re. Journal of Computational Physics, 33, 340-358.[CrossRef]
[31]	Briley, W.R. (1971) A Numerical Study of Laminar Separation Bubbles Using the Navier-Stokes Equations. Journal of Fluid Mechanics, 47, 713-736.[CrossRef]
[32]	Gupta, M.M. (1991) High Accuracy Solutions of Incompressible Navier-Stokes Equations. Journal of Computational Physics, 93, 343-359.[CrossRef]
[33]	Orszag, S.A. and Israeli, M. (1974) Numerical Simulation of Viscous Incompressible Flows. Annual Review of Fluid Mechanics, 6, 281-318.[CrossRef]
[34]	Yu, W. Q. Tao, D. S. Zhang, Q. W. W, B. (2001) Discussion on Numerical Stability and Boundedness of Convective Discretized Scheme. Numerical Heat Transfer, Part B: Fundamentals, 40, 343-365.[CrossRef]
[35]	Paramane, S.B. and Sharma, A. (2008) Consistent Implementation and Comparison of FOU, CD, SOU, and QUICK Convection Schemes on Square, Skew, Trapezoidal, and Triangular Lid-Driven Cavity Flow. Numerical Heat Transfer, Part B: Fundamentals, 54, 84-102.[CrossRef]
[36]	Pandit, S.K. (2008) On the Use of Compact Streamfunction-Velocity Formulation of Steady Navier-Stokes Equations on Geometries Beyond Rectangular. Journal of Scientific Computing, 36, 219-242.[CrossRef]
[37]	Nishida, H. and Satofuka, N. (1992) Higher-Order Solutions of Square Driven Cavity Flow Using a Variable-Order Multi-Grid Method. International Journal for Numerical Methods in Engineering, 34, 637-653.[CrossRef]
[38]	Erturk, E. (2009) Discussions on Driven Cavity Flow. International Journal for Numerical Methods in Fluids, 60, 275-294.[CrossRef]
[39]	Oosterlee, C.W., Wesseling, P., Segal, A. and Brakkee, E. (1993) Benchmark Solutions for the Incompressible Navier-Stokes Equations in General Co-Ordinates on Staggered Grids. International Journal for Numerical Methods in Fluids, 17, 301-321.[CrossRef]
[40]	Perng, C.Y. and Street, R.L. (1991) A Coupled Multigrid‐domain‐splitting Technique for Simulating Incompressible Flows in Geometrically Complex Domains. International Journal for Numerical Methods in Fluids, 13, 269-286.[CrossRef]
[41]	Erturk, E. and Dursun, B. (2007) Numerical Solutions of 2-D Steady Incompressible Flow in a Driven Skewed Cavity. ZAMM—Journal of Applied Mathematics and Mechanics, 87, 377-392.[CrossRef]
[42]	Erturk, E. and Gokcol, O. (2007) Fine Grid Numerical Solutions of Triangular Cavity Flow. The European Physical Journal Applied Physics, 38, 97-105.[CrossRef]
[43]	Shklyar, A. and Arbel, A. (2003) Numerical Method for Calculation of the Incompressible Flow in General Curvilinear Co-Ordinates with Double Staggered Grid. International Journal for Numerical Methods in Fluids, 41, 1273-1294.[CrossRef]
[44]	Louaked, M., Hanich, L. and Nguyen, K.D. (1997) An Efficient Finite Difference Technique for Computing Incompressible Viscous Flows. International Journal for Numerical Methods in Fluids, 25, 1057-1082.[CrossRef]
[45]	McQuain, W.D., Ribbens, C.J., Wang, C.Y. and Watson, L.T. (1994) Steady Viscous Flow in a Trapezoidal Cavity. Computers & Fluids, 23, 613-626.[CrossRef]
[46]	Ahmed, M. and Kuhlmann, H.C. (2012) Flow Instability in Triangular Lid-Driven Cavities with Wall Motion Away from a Rectangular Corner. Fluid Dynamics Research, 44, Article 025501.[CrossRef]
[47]	Munir, F.A., Che Sidik, N.A., Mohd, M.I., et al. (2011) Application of Lattice Boltzmann Method in Predicting Flow of Shear Driven Cavities. Journal of Mechanical Engineering and Technology, 3, 55-70.