An efficient compact difference scheme for solving the streamfunction formulation of the incompressible Navier–Stokes equations

Abstract

Recently, a new paradigm for solving the steady Navier–Stokes equations using a streamfunction–velocity formulation was proposed by Gupta and Kalita [M.M. Gupta, J.C. Kalita, A new paradigm for solving Navier–Stokes equations: streamfunction–velocity formulation, J. Comput. Phys. 207 (2005) 52–68], which avoids difficulties inherent in the conventional streamfunction–vorticity and primitive variable formulations. It is discovered that this formulation can reached second-order accurate and obtained accuracy solutions with little additional cost for a couple of fluid flow problems.

In this paper, an efficient compact finite difference approximation, named as five point constant coefficient second-order compact (5PCC-SOC) scheme, is proposed for the streamfunction formulation of the steady incompressible Navier–Stokes equations, in which the grid values of the streamfunction and the values of its first derivatives (velocities) are carried as the unknown variables. The derivation approach is simple and can be easily used to derive compact difference schemes for other similar high order elliptic differential equations. Numerical examples, including the lid driven cavity flow problem and a problem of flow in a rectangular cavity with the hight–width ratio of 2, are solved numerically to demonstrate the accuracy and efficiency of the newly proposed scheme. The results obtained are compared with ones by different available numerical methods in the literature. The present scheme not only shows second-order accurate, but also proves more effective than the existing second-order compact scheme of the streamfunction formulation in the aspect of computational cost.

Introduction

Numerical solution of incompressible Navier–Stokes equations is the most significant area in computational fluid dynamics (CFD) related fields in science and engineering. The two-dimensional (2D) incompressible flows are usually solved with the conventional primitive variable or the streamfunction–vorticity formulation of the Navier–Stokes equations. Over the past couple of decades, these formulations have been utilized by a large number of authors to test methods proposed for solving numerically a variety of fluid flow or heat transfer problems [2], [3], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [22]. The primitive variable form of Navier–Stokes equations can accurately describe the fluid flow phenomena, but a major difficult comes from the pressure term in the momentum equations and must be implicitly updated for the incompressibility to be satisfied. The streamfunction–vorticity form of Navier–Stokes equations has been used to avoid handling the pressure variable for several decades. The vorticity formulation is obtained by taking the rotational of the momentum equations. The difficulty with the vorticity formulation is the lack of the simple physical boundary conditions for the vorticity field at the no-slip boundaries. When the vorticity transport equation is numerically solved, as mentioned by some researchers [27], [28], [45], a variety of numerical approximations for the boundary values of vorticity has to be carried out. In [45], a method of vorticity–streamfunction dynamics has been developed by Ben-Artzi, Fishelov and Trachtenberg. An important corollary of this method is that the difficulty of determining the vorticity on the boundary was already avoided. Later, in order to circumvent the pitfalls associated with the vorticity values at the boundary, the streamfunction–velocity or the pure streamfunction formulation-based methodology for the solution of the 2D incompressible fluid flows has been utilized by some authors [27], [28], [29], [30], [31], [46], [47]. The main advantage of this type of formulation is that the boundary conditions of streamfunction and velocity are generally known and are easy to implement computationally, thus the computational schemes proposed are very efficient and accurate. The scheme developed in this paper belongs to this category.

We consider the following streamfunction–velocity formulation of the 2D steady Navier–Stokes equations representing incompressible fluid flows:

$$\frac{{\partial }^4\psi }{\partial x^4}+2\frac{{\partial }^4\psi }{\partial x^2\partial y^2}+\frac{{\partial }^4\psi }{\partial y^4}=\mathit{Re}\left[u\left(\frac{{\partial }^3\psi }{\partial x^3}+\frac{{\partial }^3\psi }{\partial x\partial y^2}\right)+v\left(\frac{{\partial }^3\psi }{\partial x^2\partial y}+\frac{{\partial }^3\psi }{\partial y^3}\right)\right]$$

$$u={\psi }_y\text{,}v=-{\psi }_x$$

where ψ is the streamfunction, u and v are the velocities and Re is the nondimensional Reynolds number.

When the streamfunction–velocity formulation (1), which is a fourth-order partial differential equation, is solved by finite differences, a uniform grid with 13 points must be needed to obtain a classical second-order finite difference (FD) approximation. This difference discretization using 13 grid points needs to be amended at grid points near the boundaries and must bring about difficulties for the solution of the resulting linear systems [27]. In order to overcome the above drawback, Gupta and Kalita [27] proposed a new paradigm for solving the steady Navier–Stokes equations: streamfunction–velocity formulation and derived two second-order compact finite difference schemes that carry streamfunction and its first derivatives (velocities) as the unknown variables for this equation using the Mathematica code. They used a biconjugate gradient method to obtain the numerical solutions of the fluid flow problems and solved the cavity flow with a grid size of 161 × 161 for 3200 ⩽ Re ⩽ 10,000 and a problem of flow in a rectangular cavity with the hight–width ratio of 2 with a grid size of 81 × 161 for 100 ⩽ Re ⩽ 1500. Similar a second-order compact finite difference formulation for the streamfunction–velocity formulation (1) was presented by Tian et al. [33]. They used a multigrid method to obtain the numerical solutions for the lid driven cavity flow problem for Reynolds numbers up to Re = 3200 with a maximum of 257 × 257 grid mesh.

In Fig. 1, we denote the placement of nine points, and number the mesh points (x_i, y_j), (x_i + h, y_j), (x_i, y_j + h), (x_i − h, y_j), (x_i, y_j − h), (x_i + h, y_j + h), (x_i − h, y_j + h), (x_i − h, y_j − h) and (x_i + h, y_j − h) as 0, 1, 2, 3, 4, 5, 6, 7 and 8, respectively, where h denotes the mesh size of x- and y-directions. In writing the finite difference approximations a single subscript ‘k’ denotes the corresponding function value at the mesh point numbered ‘k’.

Using the Mathematica code, Gupta and Kalita [27] proposed two second-order compact difference approximations for the solution of the streamfunction formulation (1), which are given by

$$28{\psi }_0-8\sum _{k=1}^4{\psi }_k+\sum _{k=5}^8{\psi }_k=3h(v_1-v_3-u_2+u_4)+\frac{1}{2}{\mathit{Reh}}^2\left(v_0\sum _{k=1}^4{u}_k-u_0\sum _{k=1}^4{v}_k\right)$$

and

$$-112{\psi }_0+(32-2{\mathit{Rehu}}_0){\psi }_1+(32-2{\mathit{Rehv}}_0){\psi }_2+(32+2{\mathit{Rehu}}_0){\psi }_3+(32+2{\mathit{Rehv}}_0){\psi }_4-[4-\mathit{Reh}(u_0+v_0)]{\psi }_5-[4+\mathit{Reh}(u_0-v_0)]{\psi }_6-[4+\mathit{Reh}(u_0+v_0)]{\psi }_7-[4-\mathit{Reh}(u_0-v_0)]{\psi }_8=12h(u_2-u_4-v_1+v_3)+2{\mathit{Reh}}^2[u_0(v_1+v_3)-v_0(u_2+u_4)]$$

respectively.

In [33], another second-order compact difference scheme for the solving the streamfunction formulation (1) was given by

$$-112{\psi }_0+(32+4{\mathit{Rehu}}_0){\psi }_1+(32+4{\mathit{Rehv}}_0){\psi }_2+(32-4{\mathit{Rehu}}_0){\psi }_3+(32-4{\mathit{Rehv}}_0){\psi }_4-[4-\mathit{Reh}(u_0+v_0)]{\psi }_5-[4+\mathit{Reh}(u_0-v_0)]{\psi }_6-[4+\mathit{Reh}(u_0+v_0)]{\psi }_7-[4-\mathit{Reh}(u_0-v_0)]{\psi }_8=12h(u_2-u_4-v_1+v_3)$$

The above three second-order compact finite difference schemes discretize the 2D streamfunction Eq. (1) at grid point (x_i, y_j) using the eight neighboring values of ψ, the nearest neighbors values of −ψ_x (=v) and values of ψ_y (=u). Corresponding the fourth order compact approximations for −ψ_x (=v) and ψ_y (=u) are given, respectively, by

$$\frac{1}{6}{v}_1+\frac{4}{6}{v}_0+\frac{1}{6}{v}_3=\frac{{\psi }_3-{\psi }_1}{2h}$$

$$\frac{1}{6}{u}_2+\frac{4}{6}{u}_0+\frac{1}{6}{u}_4=\frac{{\psi }_2-{\psi }_4}{2h}$$

Scheme (3) may be called as the nine point constant coefficient second-order compact (9PCC-SOC) scheme; i.e., the influence coefficients of the FD formula for the streamfunction are constant. Schemes (4), (5) may be classified as the nine point variable coefficient second-order compact (9PVC-SOC) schemes; i.e., the influence coefficients of the FD formula for the streamfunction are connected to Re, u, v and h.

It is easily found that the coefficients of the schemes (4), (5) depend on the first derivatives (velocities) of the unknown variables ψ and Reh. This indicates that the resulting algebra systems of the equations arising from schemes (4), (5) are strongly nonlinear. So, the direct solvers for the solution of the resulting systems of algebraic equations can only be used for small values of Reh and the iterative methods either converge very slowly or diverge. In contrary, the algebraic systems associated with the nine point constant coefficient second-order compact scheme (3) is linear for the streamfunction ψ. However, the coefficient matrix arising from the discretization (3) is not diagonally dominant and conventional iterative methods such as Gauss–Seidel method may not obtain the convergence results. Hence to solve the system of equation associated with scheme (3), Gupta and Kalita [2] used the hybrid biconjugate gradient stabilized method (BiCGStab) [24], [23], [34].

In this paper, we develop a new second-order compact FD formulation for solving the streamfunction–velocity formulation of the incompressible Navier–Stokes equations (1), (2). The formulation also uses the streamfunction and its first derivatives (velocities) as unknown variables. The new scheme for the streamfunction ψ is established on a uniform grid using the four nearest neighbors values of ψ and is the constant coefficient second-order compact scheme, named as the five point constant coefficient second-order compact (5PCC-SOC) scheme. The present scheme gives rise to a diagonally dominant block tri-diagonal system of equations, which results in more sparse matrix than one resulted from Gupta and Kalita’s scheme (3). This will significantly decrease the computational complexity. As mentioned in Ref. [29], the idea of using the streamfunction formulation is not novel, what was missing for many years was an efficient implementation of it. We use a version of multigrid called FAS (full approximation scheme) to obtain the numerical solutions of the fluid flow problems. We employ the newly proposed 5PCC-SOC scheme and the above three second-order compact schemes (3), (4), (5) mentioned with multigrid techniques to simulate the two-dimensional square driven cavity flow for Reynolds numbers up to Re = 7500 with a maximum of 129 × 129 grid mesh. Our computational results are compared with those obtained by using the 9PCC-SOC scheme [27] and the 9PVC-SOC schemes [27], [33], and well established results from various other schemes in the literature. It is discovered that our new formulation successfully provides second-order accuracy solutions and high computational efficiency.

The remainder of this paper is organized as follows. In Section 2, a second-order compact FD scheme for the streamfunction formulation of Navier–Stokes equations is proposed using a new derivation approach. In Section 3, the solution of algebraic systems associated with the proposed FD approximations is described. Numerical experiments for two test problems are performed to validate the accuracy and efficiency of the newly derived compact difference scheme in Section 4. Conclusions are included in Section 5.