A compact streamfunction-velocity scheme for the 2-D unsteady incompressible Navier-Stokes equations in arbitrary curvilinear coordinates ＊

Jian-xin Qiu, Bo Peng, Zhen-fu Tian

Department of Mechanics and Engineering Science, Fudan University, Shanghai 200433, China

Abstract: A streamfunction-velocity formulation-based compact difference method is suggested for solving the unsteady incompressible Navier-Stokes equations in the arbitrary curvilinear coordinates, in which the streamfunction and its first derivatives as the unknown variables are utilized. Numerical examples, involving the boundary layer problem, a constricted channel flow, driven polar cavity flow and trapezoidal cavity flow problem, are solved by the present method. Numerical results demonstrate the accuracy of the proposed scheme and exhibit the numerical capability to simulate the flow problems on geometries beyond rectangular. For driven polar cavity flow problem, the results show that the flow for Re=5000 is not steady but time-periodic, and the critical Reynold number (Rec ) for the occurrence of a Hopf bifurcation is given.

Key words: Streamfunction-velocity formulation, arbitrary curvilinear coordinates, compact scheme, unsteady, incompressible flow

Introduction

Over the past few years, as a significant alternative computational technique to 2-D Navier-Stokes(N-S) equations solvers, the streamfunction-velocity or pure streamfunction approach has captured more attentions. The vorticity formulation (vorticityvelocity and streamfunction-vorticity)[1-5], due to the absence of the physical boundary conditions for the vorticity, the corresponding numerical approximation for the vorticity boundary has to be supplemented. By comparison with vorticity formulation, the boundary conditions of velocity and streamfunction in the streamfunction-velocity or pure streamfunction formulation are generally known and are not difficult to implement computationally. However, the governing equation based on the streamfunction-velocity formulation is a fourth order nonlinear partial differential equation, such that a uniform grid with 13 grid points is carried out to acquire a classical second-order finite difference (FD) scheme when the streamfunction-velocity formulation is approached via the FD method. Whence the FD scheme using 13 points has to be amended at grid points in the vicinity of the boundaries and there are further difficulties with the solution of the resulting linear systems. To eliminate these pitfalls, a variety of compact schemes with respect to the streamfunction-velocity formulation have been proposed and are very efficient and accurate[6-12]. These compact schemes are confined invariably to only rectangular regions, so that they could not fully exploit the advantages associated with non-rectangular geometric problem, that of smaller scales in the regions of large gradients in particular.As a matter of fact, many flow problems, such as large gradient boundary layer problems, have uneven distribution in various places. Hence, it is extremely indispensable to establish effective numerical algorithms for simulating the flow in geometries beyond rectangular. Fortunately, some relevant researches have begun to be performed for the streamfunctionvelocity or pure streamfunction formulation of the incompressible N-S equations based on the irregular planar domains. Yu and Tian[13]addressed a compact streamfunction-velocity scheme for a steady 2-D incompressible N-S equations via the polar coordinates transformation based on five-point stencil.Pandit[14]developed a new methodology based on streamfunction velocity formulation for a steady biharmonic equation on irregular geometries, and the gradient of streamfunction (or velocity) is replaced through the second-order central difference scheme.Shortly afterwards, Sen et al.[15]presented an implicit compact scheme for 2-D unsteady flows in which the fourth-order pure streamfunction equation is transformed by a conformal mapping. These methods mentioned above[13-15]have the following several characters. On the one hand, the methods in Refs.[14-15] have been developed based on a conformal transformation rather than the arbitrary coordinate transformation, so that it may be invalid for the non-conformal transformation cases. The compact schemes in Refs. [13-14] merely happen to be applied to steady problems. On the other hand, these discrete difference schemes are based on the nine-point compact stencil except for the literature[13]. As mentioned in Ref. [7], they are called as nine-point variable coefficient second-order compact scheme, i.e.,the influences coefficients of the finite difference formulation for the streamfunction are connected to Reynolds number, grid size and streamfunction gradient, thus making the discrete algebra systems have stronger nonlinearity. For five-point compact scheme, the coefficient matrix is constant and gives rise to a strictly dominant block tri-diagonal system of equations. In the present paper, we extend this five-point compact scheme to the arbitrary curvilinear coordinates. In addition, the element involving the discrete coefficient matrix contains the unknown streamfunction gradient so that it is variable coefficient matrix in each iteration step. Meanwhile,the coefficient matrix is not completely diagonally dominant matrix and may become ill-conditioned matrix during the iterations. As a result, one can hardly obtain convergence results using traditional iterative method.

In the present study, we establish a new fivepoint constant coefficient compact stream-functionvelocity scheme for the 2-D unsteady N-S equations based on the arbitrary curvilinear coordinates. The coefficient matrix is constant and strictly pentadiagonal dominant which can ensure stability for the algebraic system utilizing conventional iterative method and also decrease significantly the computational complexity. To check the proposed scheme,we are concerned with an analytical solution for the boundary layer, a constricted channel problem, driven polar cavity flow and trapezoidal cavity flow problem.

1. The mathematical formulation and numerical method

1.1 The mathematical formulation

Considering the streamfunction-vorticity (ψ-ζ)formulation of the unsteady incompressible N-S equations in the Cartesian coordinates (x,y) is

whereu,vare the velocities,ψ,ζare the streamfunction, vorticity respectively,Sis the source term, andReis the Reynolds number. For the sake of convenience, Eqs. (1), (2) can be expressed as the unified formulation

whereΦ1,Φ2stand for streamfunctionψand vorticityζ. The coefficients and the source terms are given by

Now, we use the coordinate transformation to extend this formulation to a general coordinate system(ξ,)η. The relationship of the coordinate transformation is following:

Under this transformation, the above Eq. (3) is transformed as

The coefficients are given by

And then, eliminating vorticityζ, the streamfunction-velocity formulation of N-S equation in the general curvilinear coordinate system is obtained by

where the coefficientsA,B,C,D,E,F,G,H,K,L,M,N,P,Q,Rgiven explicitly as follows:

In addition, the relations between the velocities(u,v) and the first derivations of the streamfunction(ξψ,ηψ) are given by

From Eq. (11), we can easily obtain the boundary conditions for the first derivations of the streamfunction. Using Eq. (10), we can determine the values of the velocities byξψ,ηψto update those coefficients in Eq. (9).

1.2 The discretization scheme

For convenience, standard central difference operators with first and second derivatives at each internal nodal point (ξi,ηj) are given by

wherehξ,hηdenote the spatial increment inx-,y-directions, 0, 1, 2, 3, 4, 5, 6, 7, 8 stand for grid points, respectively.

In Eq. (9), all derivative terms of the streamfunctionψare replaced by the first and second derivative of the streamfunction gradientsξψ,ηψexcept for the fourth derivative terms ∂4ψ/∂ξ4and∂4ψ/∂η4. And the newly governing Eq. (9) is discretized by using the standard central difference operators as follows

where the following Stephenson finite difference operators are employed to approximate the difference operators

And the following Padé schemes are chosen to approach the stream function gradientsξψ,ηψ

Meanwhile, using second-order Crank-Nicolson scheme with a temporal increment Δtfor approaching the time derivative, and the discrete scheme of Ref. [9]can be written by

The coefficients of Eq. (16) are listed as follows:

and

1.3 The algebraic system

Next, we discuss the solution of algebraic system associated with the newly proposed algorithm, and the matrix form of the system can be written as follows

(1) Begin withn,kψ.

(3) Predictψn+1,k+1using Eq. (18).

(4) Correctψn+1,k+1,ψn+1,k+1using Eq. (15).

ξ

η

0convergence criterion, then. Otherwise,setk=k+1 and go to step 3.

(6) March in time withn←n+1 and go to step 2.

2. Stability analysis

Now we prove the stability for constant coefficients of the finite difference scheme (16) by von Neumann linear stability analysis. Taking into account the complexity of the general curvilinear coordinate system, we investigate only the stability analysis of the current scheme with a conformal transformation in here. By means of linearizing, Eq. (9) is simplified as following

Herein,a＞0,b,c,dare assumed to be independent ofψand its derivativesξψ,ηψ. Ulteriorly, weighted thought is applied to Eq. (19)

Theorem 1:The finite difference scheme (16) is unconditionally stable for 0.5≤ω≤1 in the von Neumann sense whend≤0.

Proof:LetwhereΦnis the amplitude at time levelnandis imaginary unit,θξ=2πhξ/K1,θη=2πhη/K2are the phase angles inx-,y- directions with wavelengthsK1,K2respectively. Then according to Eqs. (14),(15), we have

Table 1 The L∞-errors and convergence rates of streamfunction for various λ, and Re

Using relations (21)-(24) in the scheme (16), we obtain the amplification factor (G) associated with the finite difference scheme (16)

Note from Eq. (26) thatA1≤0,A2≥0, thusA1A2≤0, and the above inequality is apparently satisfied which completes the required proof.

3. Numerical results

In this section, the accuracy, effectiveness and robustness of our scheme are tested through using the analytical issues and three classical flow problems with complicated geometry. It includes the following four test problems: the boundary layer problem, a constricted channel problem, driven polar cavity flow and trapezoidal cavity flow problem.

3.1 Boundary layer problem

Firstly, to check the accuracy of the proposed scheme, we consider a problem governed by N-S equations with a constructed analytical solution[2], the analytical solution of this problem with respect to streamfunction is expressed by

and the source termSin Eq. (2) is given by

The boundary conditions are given by analytical solution. This problem has a steep boundary layer atx=0,y=0, and the thickness of boundary layer becomes thinner as the value ofpincreases. Similar to Ref. [14], we takepas 10 in our calculation. In order to capture preferably boundary layer information, the following transformation is given by

whereλξ,λη(0≤λξ,λη≤1) are the stretching parameter inx-,y-directions which determine the degree of clustering near the corresponding boundaries, respectively. TheL∞-errors and the convergence rates for streamfunction with various stretching parameterλand Reynolds number are summarized in Table 1. Thereinto, the rate of convergence is estimated by

whereE2h,Ehrepresent the residual errors for the coarse mesh (2h) and refined mesh (h) in the computational domain, respectively. And the results obtained by our proposed scheme and Pandit[14]have been compared. It is clear from Table 1 that the convergence rates of the present scheme and Pandit's scheme[14]achieve second order accuracy as expected.Furthermore, comparing the maximum absolute errors(L∞-error), one can see that the accuracy of our proposed scheme for the uniform grid (λ=0) with any Reynolds numbers is better than the scheme of Ref. [14], meanwhile, it also illustrates that our scheme is more accurate than the scheme of Ref. [14]not only on the uniform grid but also on the nonuniform grid.

Fig. 1 Schematic diagram of laboratory flume

3.2 Constricted channel problem

The incompressible fluid flow problem in a 2-D rigid step-down constricted channel, which has been investigated numerically in Refs. [16-18], is a basic challenge to obtain the recirculation region immediately downstream of the constriction for appropriate Reynolds numbers in numerical simulation. As shown in Fig. 1(b), with the controlling parameterτincreases, the channel boundary can be gradually varied from a smooth constriction to one with a very sharp corner, the coordinate may become singularity when the controlling parameter reaches a certain critical value. Thusτis chosen in the range of 0＜τ≤1 for avoiding the coordinate singularity which appears at (ξ0,η0)=(0.56,1.05) in this problem[16]. Furthermore, to resolve the flow information accurately at the sharp corner, a large number of mesh points are arranged in the vicinity of corner. Therefore, we consider the following mapping function of independent variables

herein,Γ=cosh(2ξ)+cos(2η), and the constantλ,μare defined as

wherea,bare the outlet, inlet heights of the constricting channel and are taken asa=1,b=2 respectively. In order to exclude the influence of the location of the inlet and outlet boundaries on the flow development, the inlet and outlet boundary are set atx=-10,x=30 far away from the constricted corner respectively. And then, as shown in Fig. 1(a),the corresponding boundary conditions for the above problem are listed as following:

Fig. 2 Streamlines and vorticity contours for different Reynolds numbers with τ=0.9

It is worthy pointing out that, combining (30),(31),ywill close to (bη/ 2τ) (aη/ 2τ) whenξappro- aches negative (positive) infinity. The qualitative results in terms of the streamlines and post processed vorticity contours for the constricted channel problem withτ=0.9,Re=250, 500 are displayed in Fig. 2. The shape and value of the streamline and vorticity contour line are in good agreement with the results of Refs. [16, 18]. And the centreline velocity profiles for various Reynolds number and the controlling parameters are exhibited in Fig. 3. It is noticed that velocity profile is almost overlap for different Reynolds numbers and controlling parameters at the upstream of the throat. When the flow encounters the constriction in the vicinity of the throat, the centerline velocity profile occurs significant jump and the changed degrees increase with the controlling parameter (0＜τ≤1) increases when fixedRe=500. Forτ=0.8, the changed degree of velocity profile is negative correlation withRe(Re≤1000), the main reason is that the width of the channel becomes more and more narrow with decreasing inReor increasing inτ. In addition,the quantitative comparisons of the horizontal position of the separation and reattachment point on the upper wall vortex with various Reynolds numbers forτ=0.9 are tabulated in Table 2, it is clear that the current computational results are in good agreement with those results of Refs. [17-18].

Fig. 3 Centerline velocity profile for the constricted channel problem

Table 2 Comparison of the separation and reattachment points on the upper wall vortex for τ=0.9

3.3 Driven polar cavity flow

A fluid flow problem in a polar cavity, which was first researched by Fuchs and Tillmark[19]both experimental and numerical method, is nowadays being often employed to validate the capability of numerical algorithms for physical domains having circular arcs as boundaries. It is usually more difficult to capture accurately the flow character for some numerical algorithms than the square cavity flow problem because of its irregular geometrical shape.The schematic diagram of the polar cavity flow and boundary conditions are shown in Fig. 4, and the abbreviations V1, VL2 and VR2, VL3 refer to the primary vortex, left and right secondary vortex,tertiary vortex of the cavity flow, respectively.Meanwhile, we consider the following polar coordinate transformation

Fig. 4 Schematic diagram of the polar cavity flow and boundary conditions

In our calculation, the termination condition which the flow reaches steady states is that the absolute error tolerancewas achieved. Firstly, we present the comparisons of the steady-state streamfunction contour between our numerical and experimental results of Ref. [19] forRe=55, 350 in Fig. 5. Good agreement has been achieved between the current numerical and the experimental results. To illustrate the richness of the flow patterns, the computed streamlines and post processed vorticity contours for 0≤Re≤3000 are displayed in Fig. 6, it can be seen that two symmetrical secondary vortices about the radial line atη=π/2 are generated on the top left and right corner forRe=0. With increasing inRe, the two secondary vortices grow gradually in size, and the growth degree for the top left secondary vortexis larger than the top right one, at the same time, the primary vortex moves also towards the right side of the polar cavity center gradually. The small size tertiary vortex on the top left corner is captured forRe≥2 000, As Reynolds number further increases,the small size tertiary vortex also gradually becomes larger, and the position of primary vortex almost becomes stabilization. The corresponding vorticity contours withRe≤3000 are also shown in the right Fig. 6. All of the phenomena are accordant with those of Refs. [13, 20]. ForRe=5000, the stable state results are given in Ref. [20] based on the steady N-S equations. However, We ran our code on the grids of size 81×81, 129×129 and 257×257 for computing the flow inRe=5000. Comparisons of the total energy,amplitude and frequency (f) for different grids withRe=5000 at dimensionless timet=1 000-1100 are shown in Table 3. The total energy (E) and amplitude (A) are calculated by

Fig. 5 Comparison of the steady-state streamlines

It is clear that the differences between 81×81,129×129 and 257×257 grids about the average total energy, amplitude and frequency are no more than 1%.The partial time trace of total energy, Fourier frequency spectrum and phase-space trajectory at the monitoring point (5/4,(2π+1)/4 are displayed in Fig. 7 forRe=5000 with 129×129 grid. It indicates that the flow is not steady but time-periodic, which is characterized by the fact that the frequency spectrum shows a spike (f=0.2224) and its harmonics. We can also see from Fig. 7(c) that the phase trajectory ofuversusvconverges to a steady limit cycle,which manifests further that the flow is strictly periodic. In fact, the occurrence of a Hopf bifurcation is believed to take place at the Reynolds number lower than 3625 through our elaborate calculation. Finally the temporal evolution of streamlines forRe=5000 on 129×129 grid are plotted during one periodT0in Fig. 8.

Fig. 6 The streamfunction (a) and vorticity (b) contours with 129×129 grid for different Reynolds numbers

Table 3 Comparison of the total energy, amplitude and frequency for different grid with Re=5000

Additionally, the quantitative data associated with the character of the primary, secondary and tertiary vortices for the driven polar cavity flow problem with 0≤Re≤3000 are tabulated in Table 4.It can be noted that our results with other existing literature results[13,20]are in excellent agreement when Reynolds number is less than 2 000. As forRe=3000, there is a slight difference between the secondary and tertiary vortices.

Fig. 7 Unsteady states at Re=5000 with 129×129 grid

3.4 Driven trapezoidal cavity flow

In this section, to simulate fluid flow within isosceles trapezoidal cavity in Fig. 9, the following nonorthogonal transformation, which will map the trapezoidal domain in (x,y) plane to the rectangular domain in (ξ,)ηplane, is introduced by

Fig. 8 The temporal evolution of streamlines on grid 129×129 for Re=5000

Fig. 9 Schematic diagram of the trapezoidal cavity flow and boundary conditions

In Table 5, the minimum values of streamfunction (ψmin) the central locations (x0,y0) of the primary vortex and post processed vorticity (ω0) at these locations are compared with the results of Mcquain et al.[21]and Zhang et al.[22]As seen from Table 5, the present computational results are in good agreement with those results of Refs. [21-22].Meanwhile, streamline patterns for the trapezoid withRe=100, 500 are shown in Fig. 10, it is clear that the shapes of the streamline agree excellently with the available results of Ref. [21]. Figure 11 shows the computed streamlines patterns for 100≤Re≤1000.ForRe=100, the central eddy is seen to occupy the whole region of the trapezoidal cavity. As the Reynolds number is increased to 400, the two secondary circulations occur at the bottom corner regions and the left secondary eddy is greater than the right one in size. With increase in Reynolds number,the two secondary eddies coalesce into one. AsReis increased further, the secondary eddy gains a more significant size, it can be seen from Fig. 11(d) that the primary vortex is compressed and splits into two individual vortices.

Table 4 The results of primary, secondary and tertiary vortices with different Reynolds numbers

4. Conclusions

In the present paper, a temporal and spatial second-order accurate compact difference scheme for solving the time-dependent incompressible N-S equations based on streamfunction-velocity (pure streamfunction) formulation is proposed in the arbitrary curvilinear coordinates. The main advantages of the current algorithm are that it is established based on the arbitrary curvilinear coordinates and the coefficient matrix based on five-point compact stencil is constant and satisfies strictly penta-diagonal dominant which can ensure stability for the discrete algebraic system utilizing conventional iterative method. Compared with the established numerical and experimental results, our present algorithm has been verified that it is quite effective and the numerical results agree excellently with the available results in literatures.

Table 5 Validation of present simulations for trapezoidal cavity with Re=100,500

Fig. 10 Streamline patterns for the trapezoid

We apply current approach to four problems which encompass an analytical solution for the boundary layer, constricted channel problem, driven polar cavity flow and trapezoidal cavity flow problem.The robustness and effectiveness of the scheme are demonstrated in the above cases. The relevant tabular data imply that the present compact scheme achieves the theoretical convergence accuracy as expected and is slightly better than numerical results for other available second order schemes. For the driven polar cavity flow problem, the relevant data and graphs show that our current computed solutions are in excellent agreement with the existing numerical results for 0≤Re≤3000. ForRe=5000, as displayed in Fig. 7, the frequency spectrum shows a spike(f=0.2224) and its harmonics, the phase trajectories ofuversusvconverge to a steady limit cycle. The facts show that the corresponding flow withRe=5000 is not steady buttime-periodic flow with a period close to 4.5. Meanwhile, the occurrence of a Hopf bifurcation is believed to take place at the Reynolds number lower than 3 625 through our elaborate calculation. For the driven trapezoidal cavity flow case, which is the nonorthogonal transformation problem, our computational results are in excellent agreement with the established numerical results.

Fig. 11 Streamline and vorticity contours for the trapezoid with various Reynolds numbers

Acknowledgement

This work was supported by the China Postdoctoral Science Foundation (Grant No. 2014M550211).