Fourth-order nite dierence scheme for the incompressible
Navier–Stokes equations in a disk
Ming-Chih Lai
∗;†Department of Applied Mathematics; National Chiao Tung University; 1001; Ta Hsueh Road; Hsinchu 30050; Taiwan
SUMMARY
We develop an ecient fourth-order nite dierence method for solving the incompressible Navier– Stokes equations in the vorticity-stream function formulation on a disk. We use the fourth-order Runge– Kutta method for the time integration and treat both the convection and diusion terms explicitly. Using a uniform grid with shifting a half mesh away from the origin, we avoid placing the grid point directly at the origin; thus, no pole approximation is needed. Besides, on such grid, a fourth-order fast direct method is used to solve the Poisson equation of the stream function. By Fourier ltering the vorticity in the azimuthal direction at each time stage, we are able to increase the time step to a reasonable size. The numerical results of the accuracy test and the simulation of a vortex dipole colliding with circular wall are presented. Copyright ? 2003 John Wiley & Sons, Ltd.
KEY WORDS: Navier–Stokes equations; vorticity-stream function formulation; polar co-ordinates; fast Poisson solver; Runge–Kutta method
1. INTRODUCTION
For the past three decades, nite dierence approximation of the incompressible Navier– Stokes equations on Cartesian geometry has been developed and implemented extensively. However, for the uid ow in polar=cylindrical or spherical domain, the relevant literature is still quite limited. The reason may be attributed to introduce some new additional diculties when the governing equations are expressed by those orthogonal curvilinear co-ordinates. For instance, those diculties might include the coupling of velocity components, the co-ordinate singularities, and the restrictive CFL constraint. Let us describe those diculties briey in the following.
∗Correspondence to: M.-C. Lai, Department of Applied Mathematics, National Chiao Tung University, 1001,
TaHsueh Road, Hsinchu 30050, Taiwan.
†E-mail: mclai@math.nctu.edu.tw
Contract=grant sponsor: National Science Council of Taiwan; contract=grant number: NSC-91-2115-M-009-016
The coupling of velocity components is not unusual even in the case of Cartesian geometry. This can be seen through the non-linear convection term in the governing equations. For the case of other orthogonal curvilinear co-ordinates systems, the velocity components are coupled even in the linear diusion terms. This can be explained by when the vector elliptic equation is written in those co-ordinates, the Laplacian acting on the vector eld results a coupled system of elliptic equations. Fortunately, this kind of coupling can be eliminated by applying a similarity transformation to the vector eld as described in Reference [1] for the polar=cylindrical co-ordinates. A detailed summary for how to separate the vector elliptic equations in other dierent co-ordinates can be found in the book of Quartapelle [2].
The second diculty arises from the co-ordinate singularities occurring at the origin=centre-line (r = 0) of the polar=cylindrical domain, and the north and south poles ( = 0; ) of spherical domain. It is important to note that the occurrence of those singularities is due to the representation of the governing equation in those co-ordinates and it is irrelevant to the solution at those singular points. For simplicity, we will just name those singularities as the ‘poles’ in dierent domains. In order to capture the solution behaviour near the poles, some suitable treatments must be employed. For instance, in Reference [3], the authors deduced a new governing equation which is valid at the centerline in polar=cylindrical domain. This traditional technique is to use a regular grid with imposing conditions at the poles. Some natural pole conditions can be found in the book of Canuto et al. [4]. Note that, in nite dierence setting, these pole conditions can be used to approximate numerical boundary values at the poles.
Another technique is to use a uniform grid by shifting a half mesh away from the poles to avoid placing the grid points directly at the poles; thus, no pole conditions are needed References [5, 6]. This simple and elegant approach has been successfully applied to simulate the convective ow over a sphere [5] and compressible Navier–Stokes equations [6]. Using this kind of grid, Lai and Wang have been able to develop an ecient class of second- and fourth-order accurate fast direct solvers for Poisson=Helmholtz equations on 2D polar and spherical domains without pole conditions [7, 8].
The third diculty arising is the severe CFL stability restriction. This is because when using a uniform grid to discretize the governing equations in polar=cylindrical co-ordinates, the azimuthal grid resolution is proportional to the distance from the poles. In other words, the physical mesh spacing of the azimuthal grid points becomes ner and ner as the poles were approached. Thus, near the poles, this leads to an over-resolution in the azimuthal direction which requires choosing a very small time step to full the CFL constraint for the time-dependent problem if an explicit time-integration is used. Several techniques were proposed to overcome this diculty. For instance, Akselvoll and Moin [9] treated all terms with derivatives in the azimuthal direction implicitly in the core region (near the pole) but all other terms are treated explicitly. While in the outer region (away from the pole), the diusive terms and the derivatives in the radial direction are treated implicitly but all other terms are treated explicitly. The idea behind their approach is signicant, however, the treatment of the interface between the two regions and the maintenance of the required accuracy could be quite complex. Another simple treatment introduced by References [5, 10] is to lter out the high wavenumber modes of the dependent variable in the azimuthal direction. This treatment has the same eect to coarsen the grid near the poles thus allowing us to increase the time step to a reasonable size.
In this paper, we shall present an ecient nite dierence scheme for solving the incom-pressible Navier–Stokes equations in a two-dimensional disk. The reason that we choose the
ow problem in a disk domain is because it suces to have all aforementioned numerical dif-culties occurring in higher-dimensional cylindrical and spherical geometries. Note that here we are not trying to compare or compete with the nite element or spectral methods for solving such kind of problems. We simply want to introduce some simple nite dierence treatments for the diculties such as the co-ordinate singularity and the restrictive CFL con-straint. Another main feature of our method is the fourth-order accuracy which is sucient for most of ow computations. The rest of the paper is organized as follows. In Section 2, we formulate the incompressible Navier–Stokes equations using the vorticity-stream function form in polar co-ordinates. The complete fourth-order numerical scheme for solving the equations is described in Section 3. A simple fourth-order fast direct solver for Poisson equation in a disk is reviewed in Section 4. The numerical results are given in Section 5 and followed by some conclusion.
2. INCOMPRESSIBLE NAVIER–STOKES EQUATIONS
The incompressible Navier–Stokes equations has the standard form
@u
@t +u· ∇u +∇p = u (1)
∇ ·u = 0 (2)
where u(x; t) is the uid velocity, p(x; t) the pressure, and the uid viscosity. The rst equation describes the conservation of momentum and the second one is the conservation of mass. In 2D geometry, we can express the Navier–Stokes equations (1) – (2) by so called the vorticity-stream function formulation. By taking the curl of Equation (1) to eliminate the pressure gradient term, we have
@!
@t + J (!; ) = ! (3)
where ! is the non-zero z component of the vorticity, the stream function, dened by u =∇ ×ez, and J (!; ) is the Jacobian determinant. Using this formulation, the velocity
u automatically satises the incompressibility constraint (2). By the denition of !, we can easily derive the relation of ! and as
=−! (4)
Therefore, the original 2D Navier–Stokes equations (1) – (2) with three primitive variables now has an alternate formulation described by (3) – (4) with only two unknown variables.
We are interested in the numerical approximations of Equations (3) – (4) in a unit disk ={0¡r61; 0662}. Thus, it is natural to rewrite the equation by the polar co-ordinates. The non-linear Jacobian describing the vorticity transport now has the form
J (!; ) =1 r @ @ @! @r − @ @r @! @ (5)
and the Laplacian operator is ≡ @ 2 @r2 + 1 r @ @r + 1 r2 @2 @2 (6)
The radial and azimuthal velocity components can be recovered from the stream function by the formulas
ur=1r @ @; u=−@ @r (7)
In this paper, we restrict our attention to the ow inside a unit disk with the no-slip velocity conditions specied on the boundary; that is, ur= u= 0 at r = 1. From the relation of (7),
the boundary conditions for Equations (3) – (4) become
(1; ) = 0; @
@r (1; ) = 0 (8)
So the governing equations consist of Equations (3) – (4) and the boundary conditions (8). It is interesting to note that there are two boundary conditions for the stream function but no boundary condition for the vorticity !. This is troublesome since Equations (3) – (4) are also coupled by the boundary conditions. There are several dierent approaches summarized in Reference [2] to overcome this diculty. In this paper, we will just use an explicit time integration to solve Equations (3) – (4), and transfer one boundary condition of the stream function to the vorticity.
3. NUMERICAL METHOD
3.1. Time integration
We employ a fourth-order Runge–Kutta method for the time integration of Equations (3) – (4) as !1−!n 1 2t + J (!n; n) = !n; 1=−!1; 1(1; ) = 0 (9) !2−!n 1 2t + J (!1; 1) = !1; 2=−!2; 2(1; ) = 0 (10) !3−!n t + J (!2; 2) = !2; 3=−!3; 3(1; ) = 0 (11) K4=−J (!3; 3) + !3 !n+1 =1 3(−! n + !1+ 2!2+ !3) + t 6 K4 (12) n+1=−!n+1; n+1(1; ) = 0
At each Runge–Kutta stage, we discretize both the convection and diusion terms in the vorticity equation explicitly such that the vorticity can be computed directly (without solving any linear system of equations). The spatial derivatives with respect to r and in both convection and diusion terms are approximated by the standard fourth-order centred dierence discretization, see next subsection in detail. Meanwhile, to avoid the severe CFL stability constraint near the origin, increasingly many high wavenumber modes for the vorticity in direction are removed as the grid point toward to the origin. The resultant vorticity is then used as a right-hand side function of Poisson equation for solving the stream function. The scheme is very ecient in terms of computational complexity since it only involves one Poisson equation to be solved at each time stage. In Section 4, we will introduce a fourth-order accurate fast Poisson solver for solving the stream function. The employment of Fourier ltering will be discussed more detail in Section 5.
3.2. Spatial discretization
First, let us choose a M×N computational grid inside the disk domain with the grid point locations
G ={(ri; j) = ((i−1=2)r; (j−1))|16i6M; 16j6N} (13)
where r = 2=(2M + 1) and = 2=N . By the choice of r, we have rM +1= 1. By shifting
a half mesh away from the origin, we avoid placing the grid point directly at the origin; thus, hopefully no pole approximation is needed. Because a ve-point stencil will be used for approximating the derivatives, we dene a larger grid set by adding the computational grid G with a few ghost points along the radial direction as
G ={(ri; j) = ((i−1=2)r; (j−1))| −16i6M + 2; 16j6N} (14)
The unknowns are computed at the grid points of G whereas G consists of all grid points involved. This polar grid G is illustrated in Figure 1 with M = 8 and N = 16.
The computation of the non-linear convection and linear diusion terms in the vorticity equation involves the spatial discretization for the rst and second derivatives in radial and azimuthal directions. Those derivatives can be approximated by the standard fourth-order centred dierence formulas as
@! @r
ij≈
−!i+2; j+ 8!i+1; j−8!i−1; j+ !i−2; j
12r (15) @2! @r2 ij≈
−!i+2; j+ 16!i+1; j−30!ij+ 16!i−1; j−!i−2; j
12r2 (16)
Similarly, we can write down the same formulas for the approximations of @!=@ and @2!=@2.
Since the function dened in a disk is periodic in , the approximation of -derivative does not have any numerical boundary value problem. By this, we mean that !ij= !ij,
where j≡j (mod N ). However, the numerical boundary values for the approximation of r-derivatives are more complicated. First, the centred dierence formulas (15)–(16) are ap-plied to all interior points except the last radial interior point (rM; j) adjacent to the boundary.
Figure 1. A polar grid in a disk based on G with M = 8 and N = 16. The circle denotes the boundary. The vorticity ! are dened on those ‘.’ grid points. The stream function
are dened on ‘.’ points as well as the outer ghost points ‘+’.
At the point (rM; j), we use the following one-sided dierence approximations described in
Reference [11] @! @r M; j≈ 3!M +1; j+ 10!M; j−18!M −1; j+ 6!M −2; j−!M −3; j 12r (17) @2! @r2 M; j≈11!M +1; j−20!M; j+ 6!M −1; j+ 4!M −2; j−!M −3; j 12r2 (18)
Second, due to the ve-point stencil of (15) – (16), the numerical boundary approximations in the vicinity of the origin must also be provided. For instance, when i = 1 in the formula of (15), the numerical approximations of !−1; j and !0j at the corresponding ghost points
(r−1; j); (r0; j) need to be provided. Those approximations can be found in the following.
When we replace r by −r and by + in the Cartesian-polar transformation, the Cartesian co-ordinates of a point remain xed. Therefore, any scalar function satises !(−r; ) = !(r; +) if the domain of the function is extended to negative values of r. So we have
!0j= !(−r=2; j) = !(r=2; j+ ) = !1; j+N=2 (19)
3.3. Vorticity boundary condition
As we mentioned before, the no-slip conditions become two boundary conditions of the stream function as in (8). This causes that the Poisson equation of the stream function is over-determined. One easy way to overcome this diculty is to translate one boundary condition of the stream function to the vorticity, which is rst derived by Thom [12]. A very detailed discussion of the vorticity boundary condition in Cartesian co-ordinates and dierent variants of Thom’s formula can be found in the work of E and Liu [11]. Next, we will derive the fourth-order vorticity boundary condition in polar coordinates.
As in the scheme (9) – (12), the Poisson equation of the stream function is solved with the Dirichlet boundary condition (1; ) = 0. The other Neumann boundary condition of can be used to obtain the boundary condition of the vorticity as follows. Using the one-sided dierence formula (17) to approximate the Neumann condition on the boundary points (rM +1; j), we have @ @r M +1; j≈ 3 M +2; j+ 10 M +1; j−18 M; j+ 6 M −1; j− M −2; j 12r = 0 (20)
Since M +1; j= 0 for all j, we obtain
M +2; j= 6 M; j−2 M −1; j+ 13 M −2; j (21)
Using the fact that M +1; j= 0 and @ =@r|M +1; j= 0, the boundary vorticity becomes
!M +1; j=−( )M +1; j=− @ 2 @r2 M +1; j (22)
Applying the one-sided dierence (18) and using the approximation of (21), we obtain
!M +1; j=−108 M; j+ 27 M −1; j−4 M −2; j
18r2 (23)
Therefore, the vorticity at the boundary can be approximated by the neighbouring values of the stream function.
It is interesting to note that the vorticity boundary condition derived here (under the assumption of (8)) is exactly the same one rst derived by Briley [13] in the case of Cartesian co-ordinates. Although our derivation is for the no-slip boundary conditions, the extension to slip boundary conditions can be derived straightforwardly. In fact, in our numerical test in Section 5, we have implemented the slip boundary conditions for Example 1.
4. REVIEW OF FOURTH-ORDER FAST POISSON SOLVER
In the numerical scheme (9) – (12), we know that at each time stage requires solving one Dirichlet Poisson equation for the stream function. In this section, we review a simple FFT-based fast Poisson solver in a disk [8] which is applied to our numerical computations. This solver is based on truncating the Fourier series expansion, and then solving the dierential equations of Fourier coecients by nite dierence discretization. Note that, this kind of
approach is not a new one, since we can nd it in dierent literature, e.g. Reference [14]. However, our method diers from the one in Reference [14] by two folds: without pole conditions and fourth-order accuracy [8] vs with pole conditions and second-order accuracy [14]. One can see that our method has a very simple treatment for the co-ordinate singularity and it can be applied to dierent boundary conditions as well.
The Poisson equation with Dirichlet boundary in a unit disk can be written as @2u @r2 + 1 r @u @r + 1 r2 @2u @2 = f(r; ); u(1; ) = g() (24)
Due to the periodicity of u in direction, we can approximate the solution by the truncated Fourier series as
u(r; ) = N=2−1
n=−N=2un(r) e
in (25)
where un(r) is the complex Fourier coecient given by
un(r) =N1 N
j=1u(r; j) e
−inj (26)
Here, j is a uniform grid in [0; 2] which is dened exactly the same way as in (14).
The above transformation between the physical space and Fourier space can be eciently performed using the fast Fourier transform (FFT) with O(N log2N ) arithmetic operations. Substituting the expansion in (25) into Equation (24), and equating the Fourier coecients, the un(r) satises the boundary value problem
d2u n dr2 + 1 r dun dr − n2 r2 un= fn; un(1) = gn (27)
Here, the complex Fourier coecients fn(r) and gn are dened in the same manner as
Equations (25) –(26).
Now let us denote the discrete values Ui≈un(ri) and Fi≈fn(ri), where ri is chosen exactly
the same way as in (14). We now discretize Equation (27) using the fourth-order dierence formulas given in (15) – (16) at the interior points i = 1; 2; : : : ; M as
−Ui+2+ 16Ui+1−30Ui+ 16Ui−1−Ui−2
12r2
+−Ui+2+ 8Ui+1−8Ui−1+ Ui−2 12rir −
n2
r2 i
Ui= Fi (28)
This is a penta-diagonal linear system which can be solved by O(M ) arithmetic operations. In order to close the above linear system, we need to supply numerical boundary values such as U−1; U0, UM +1 and UM +2. The value UM +1 is simply given by gn. The inner numerical
which are
U0= un(r0) = un(−r=2) = (−1)nun(r=2) = (−1)nU1 (29)
U−1= un(r−1) = un(−3r=2) = (−1)nun(3r=2) = (−1)nU2 (30)
The outer numerical boundary value UM +2 can be obtained as follows. Since Equation (27)
also holds at the boundary rM +1= 1, we have,
u
n(rM +1) + un(rM +1) = n2un(rM +1) + fn(rM +1) (31)
Substituting the one-sided dierence formulas (17) – (18) into the above equation and after some careful calculations, we obtain a formula for UM +2 in terms of UM +1; UM; UM −1; UM −2
and FM +1. Thus, all the necessary numerical boundary values are derived. We conclude that
this fast solver is spectral accurate in the direction while is fourth-order accurate in the r direction.
5. NUMERICAL RESULTS
In this section, we rst demonstrate the rate of convergence (or the order of accuracy) for our fourth-order scheme (9) – (12) together with the spatial discretization and the vorticity boundary conditions outlined in Section 3. The Poisson equation of the stream function at each Runge–Kutta step is solved by the fourth-order fast direct solver described in Section 4. The second numerical test for our new scheme is to simulate a vortex dipole colliding with a circular wall [15].
Example 1 (Accuracy check)
We start our numerical tests by checking the accuracy of our scheme for a computation in a unit disk ={0¡r61; 0662}. We have taken the exact solution for Navier–Stokes equations as
!(x; y; t) = 2 e−2tcos x cos y; (x; y; t) = e−2tcos x cos y (32)
Using the polar co-ordinate transformation, the above exact solution can be represented by the functions of r and . Notice that, the above example has the non-zero vorticity at the origin (x = y = 0), but the velocity components are zero at this particular point. In fact, we did another numerical test which has non-zero velocity components and got the similar numerical results as in Table I so we omitted here.
In our test, we choose N=2×N grid points in the disk so that there are N=2 points in the radial direction and N points in the azimuthal direction. The time steps for the cases of N = 32; 64; 128; 256 are t = 0:05; 0:025; 0:0125; 0:00625, respectively. The viscosity is chosen as = 0:001. The approximate solution was computed at T = 3. Table I shows the L∞; L1,
and L2 errors for dierent number of grid points. In addition to comparing the errors of the
vorticity (!) and stream function ( ), we also list the errors of the velocity components (ur and u). One can clearly see that the fourth-order accuracy of our scheme can be indeed
Table I. Grid renement analysis at T = 3.
N L∞ error Rate L1 error Rate L2 error Rate
! 32 8.8196E-05 3.5668E-05 4.1421E-05
64 1.4197E-05 2.64 2.7410E-06 3.70 4.1679E-06 3.31 128 6.8231E-07 4.38 1.2070E-07 4.51 1.8914E-07 4.46 256 3.5084E-08 4.28 5.0087E-09 4.59 7.5868E-09 4.64
32 4.8863E-07 6.1366E-07 4.0599E-07
64 1.9435E-08 4.65 1.9195E-08 5.50 1.3896E-08 4.87 128 7.8116E-10 4.64 9.3920E-10 4.35 6.5151E-10 4.41 256 6.1184E-11 3.67 6.8173E-11 3.78 5.0133E-11 3.70 ur 32 7.4158E-04 5.9268E-04 5.0019E-04
64 5.3446E-05 3.79 4.3944E-05 3.75 3.5112E-05 3.83 128 3.5364E-06 3.92 2.9198E-06 3.91 2.2929E-06 3.94 256 2.2664E-07 3.96 1.8705E-07 3.96 1.4597E-07 3.97 u 32 7.6869E-06 2.7102E-06 3.1914E-06
64 5.0719E-07 3.92 1.5524E-07 4.13 1.6942E-07 4.24 128 2.4950E-08 4.35 7.1842E-09 4.43 7.0078E-09 4.60 256 1.1435E-09 4.45 3.7646E-10 4.25 3.2495E-10 4.43
Example 2 (Vortex dipole colliding with a circular wall)
In order to demonstrate the performance of our new fourth-order scheme, we consider the problem of a vortex dipole colliding with a circular wall. This example was taken directly from Reference [15] where the authors used it to test their pseudospectral code. Notice that, we are not trying to compare with the pseudospectral method in Reference [15] but just simply introducing an alternative nite dierence approach to the simulation of the uid dynamical problem in a disk. In particular, our approach should be valid for moderate Reynolds number ows since a fourth-order explicit time stepping procedure is used.
The vortex dipole consists of two point vortices with equal magnitude of strength but dierent signs. Those two vortex centres are initially located at the two points (0:15; 0) and (−0:15; 0) inside a unit disk. The smoothing initial vorticity is chosen as
!(x; y; 0) = 1:5e−20((x−0:15)2+y2)
−1:5e−20((x+0:15)2+y2) (33) whose contour lines can be found in the rst plot of Figure 2. The no-slip boundary conditions are imposed on the wall. The viscosity is chosen as 2×10−5.
The calculation was performed on a M×N = 512×512 polar grid. Since the uniform grid is used, the physical mesh spacings along the grid circles in the core region (the region near the origin) are very small. This has the consequence that the CFL constraint is very restrictive for the time-dependent problem. One simple treatment is to employ the Fourier ltering [5, 10] for the vorticity in the direction inside the core region at each Runge–Kutta time stage. This treatment has the same eect to coarsen the physical mesh near the origin thus allowing us to increase the time step to a reasonable size. The idea is to remove the high wavenumber modes of the vorticity gradually as the grid points toward to the origin. For simplicity, we choose the core region as the half disk r60:5. At the inner grid circle r = r=2, we keep the modes
Figure 2. Evolutionary vorticity contours for a vortex dipole colliding with a circular wall. The positive vorticity denotes by ‘-’ and the negative vorticity
by ‘−’. The magnitude ranges from −2 to 2.
with wavenumber n = 0; 1;−1 and remove the Fourier modes whose wavenumbers are greater than one. At the second grid circle r = 3r=2, we keep two more modes with wavenumber n = 0; 1; 2;−1;−2 and remove the rest of high modes. Continuing in such way, the Fourier modes are no longer to be removed outside the core region. By applying this Fourier ltering, we are able to choose a much larger time step to ensure the numerical stability. In our present computation, we choose the time step t = 2:5×10−3 which is about ten times larger than the one used in Reference [15].
Figure 2 shows the evolutionary contour plots for the vorticity. The solid line denotes the positive vorticity and the dash line denotes the negative values. The magnitude of those plots ranges from−2 to 2. One can see the vortex dipole collides with the circular wall then splits. Each vortex generates another dierent sign of vortex and forms another dipole and moves to the centre. Then they exchange the partners to form two symmetric vortex dipoles which the smaller one moves upward and the larger one moves downward to repeat the colliding process. Figure 3 shows an enlarging part of the vorticity plot at T = 120. One can see that there are four vortex dipoles totally after two complete collisions.
The energy (E) and enstrophy () are the two interesting quantities measured in this test. They are dened by
E =1 2 1 0 2 0 (u2 r+ u2) r dr d; = 1 0 2 0 !2r dr d (34)
We compute those quantities via a discrete approximation of the above integral based on the nodes described in (13). Thus, the computation is almost the same as the usage of the midpoint rule for the double integral integration. Figure 4(a) and (b) are the time evolutionary
Figure 3. Four vortex dipoles were seen after two complete dipole collisions at T = 120.
plots for the energy and enstrophy, respectively. One can see that the energy is decreasing. The enstrophy is oscillating with the maximum occurring around at T = 30 when the dipole completes the rst hitting to the wall and generates another vortices.
As in Reference [15], we also calculate the energy decay rate by the formula dE
dt =− (35)
The relative error of the decay rate is computed by
rel(t) =(E(t + t)−E(t−t))=2t + (t)
−(t) (36)
Figure 4(c) shows the relative error of the energy decay rate which is within 1.6E-03 after some earlier time. The computed error at T = 120 is 4.3E-05 which is about the magnitude of O(t2). So the numerical evidence shows that our scheme does provide the right energy
decay rate.
6. CONCLUSION
Let us summarize our scheme as follows. The Navier–Stokes equations are written in the vorticity-stream function form. We use the fourth-order Runge–Kutta method for the time integration and treat both the convection and diusion terms explicitly. Using a uniform grid with shifting a half mesh away from the origin, we avoid placing the grid point directly at the origin; thus, no pole approximation is needed. Besides, on such grid, a fourth-order fast direct
0 20 40 60 80 100 120 1.8 2 2.2 2.4 2.6 2.8 3x 10 -3 Time Energy Energy 0 20 40 60 80 100 120 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 Time Enstrophy Enstrophy 0 20 40 60 80 100 120 -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02 Time Relative error
Relative error in energy
(a) (b)
(c)
Figure 4. The evolutionary plot of energy, enstrophy, and the relative error of energy decay rate.
method is used to solve the Poisson equation of the stream function. By Fourier ltering the vorticity in the azimuthal direction at each time stage, we are able to increase the time step to a reasonable size. The numerical results of the accuracy test and the simulation of a vortex dipole colliding with circular wall are presented.
REFERENCES
1. Manna M, Vacca A. An ecient method for the solution of the incompressible Navier–Stokes equations in cylindrical geometries. Journal of Computational Physics 1999;151:563–584.
2. Quartapelle L. Numerical Solution of the Incompressible Navier–Stokes Equations. Birkhauser: New York, 1993.
3. Grin MD, Jones E, Anderson JD. A computational uid dynamic technique valid at the centerline for non-axisymmetric problems in cylindrical co-ordinates. Journal of Computational Physics 1979;30:352–360. 4. Canuto C, Hussaini MY, Quarteroni A, Zang TA. Spectral Methods in Fluid Dynamics. Springer: Berlin, 1988;
90–91.
5. Fornberg B, Merrill D. Comparison of nite dierence and pseudospectral methods for convective ow over a sphere. Geophysical Research Letters 1997;24:3245–3248.
6. Mohseni K, Colonius T. Numerical treatment of polar co-ordinate singularities. Journal of Computational Physics 2000;157:787–795.
7. Lai MC. A note on nite dierence discretizations for Poisson equation on a disk. Numerical Methods for Partial Dierential Equations 2001; 17:199–203.
8. Lai MC, Wang WC. Fast direct solvers for Poisson equation on polar and spherical geometries. Numerical Methods for Partial Dierential Equations 2002; 18:56–68.
9. Akselvoll K, Moin P. An ecient method for temporal integration of the Navier–Stokes equations in conned axisymmetric geometries. Journal of Computational Physics 1996;125:454–463.
10. Williamson DL. Linear stability of nite-dierence approximations on a uniform latitude-longitude grid with Fourier ltering. Monthly Weather Review 1976;104:31–41.
11. E W, Liu JG. Vorticity boundary condition and related issues for nite dierence schemes. Journal of Computational Physics 1996;124:368–382.
12. Thom A. The ow past circular cylinders at low speeds. Proceedings of the Royal Society, London, Series A 1933;141:651–666.
13. Briley WR. A numerical study of laminar separation bubbles using the Navier–Stokes equations. Journal of Fluid Mechanics 1971;47:713–736.
14. Iserles A. Numerical Analysis of Dierential Equations. Cambridge University Press, Cambridge, 1996; 256–260.
15. Torres DJ, Coutsias EA. Pseudospectral solution of the two-dimensional Navier–Stokes equations in a disk. SIAM Journal of Scientic Computing 1999; 21:378–403.
