Fast direct solver for the biharmonic
equation on a disk and its application
to incompressible flows
Ming-Chih Lai
a,*, Hsi-Chi Liu
baDepartment of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road,
Hsinchu 300, Taiwan
bDepartment of Mathematics, National Chung Cheng University, Minghsiung, Chiayi 621, Taiwan
Abstract
We develop a simple and efficient FFT-based fast direct solver for the biharmonic equation on a disk. The biharmonic equation is split into a coupled system of harmonic problems. We first use the truncated Fourier series expansion to derive a set of coupled singular ODEs, then we solve those singular equations by second-order finite difference discretizations. Using a radial grid with shifting a half mesh away from the origin, we can handle the coordinate singularity easily without pole conditions. The Sherman– Morrison formula is then applied to solve the resultant linear system in a cost-efficient way. The computational complexity of the method consists of O(MN log2N) arithmetic
operations for M· N grid points. The numerical accuracy check and some applications to the incompressible Navier–Stokes flows inside a disk are conducted.
2004 Elsevier Inc. All rights reserved.
Keywords: Biharmonic equation; Polar coordinates; Sherman–Morrison formula; FFT; Vorticity stream function formulation
0096-3003/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.04.064
*
Corresponding author.
E-mail address:mclai@math.nctu.edu.tw(M.-C. Lai).
1. Introduction
The biharmonic problem arises from several applications in solid mechanics and fluid mechanics. In solid mechanics, finding the displacement of the bend-ing of elastic plates involves solvbend-ing the biharmonic equation. In fluid mechan-ics, the stream function of incompressible Stokes flow in two-dimensional space is the solution of a biharmonic equation as well.
The biharmonic problem (of the first kind) on a two-dimensional domain X has the form
D2uðx; yÞ ¼ f ðx; yÞ in X; ð1:1Þ u¼ gðx; yÞ; ou
on¼ hðx; yÞ on oX; ð1:2Þ
where the differential operator D is the well-known Laplacian defined by D¼ o2
ox2þo 2
oy2and D
2
u = D(D u). As mentioned in[2], direct discretization of this fourth-order biharmonic equation (1.1) leads to a very ill-conditioned linear equations. Thus, most of the iterative methods require a large number of iter-ations in order to obtain some satisfactory solutions. One popular approach to avoid solving such ill-conditioned matrix equations is to introduce an auxiliary variable v(x,y) = Du(x,y) and to split the biharmonic equation(1.1)into a cou-pled system of Poisson equations as
Duðx; yÞ ¼ vðx; yÞ; Dvðx; yÞ ¼ f ðx; yÞ: ð1:3Þ The boundary conditions(1.2)are still the same. One can easily see that under this formulation, there are two boundary conditions for the solution u but no boundary conditions for v. The functions u and v are coupled through the boundary conditions implicitly which turns out to be the main difficulty of solving such problem.
Throughout this paper, we are interested in the domain of a unit disk X= {(x,y) : x2+ y2< 1}. Therefore, it is natural to apply the polar coordinate transformation x = r cos h, y = r sin h to the equations. For simplicity, we use the same notations to represent the functions both in Cartesian and polar coor-dinates. The coupled system of the biharmonic equation of u(r, h) can be writ-ten as o2u or2þ 1 r ou orþ 1 r2 o2u oh2¼ vðr; hÞ; 0 < r < 1; 0 6 h < 2p; ð1:4Þ o2v or2þ 1 r ov orþ 1 r2 o2v oh2¼ f ðr; hÞ; 0 < r < 1; 0 6 h < 2p; ð1:5Þ uð1; hÞ ¼ gðhÞ; ou orð1; hÞ ¼ hðhÞ: ð1:6Þ
Various approaches for the numerical solution of the boundary value prob-lem(1.4)–(1.6), along with the applications to the steady incompressible flow inside circular geometries have been developed in the literature. Those include the spectral method[6], the integral equation method[4], and the spectral/dif-ference method[1,5]. In this paper, we shall develop a fast direct method sim-ilar to [1]but differs with the treatments of coordinate singularity (with pole conditions[1]vs. without pole conditions) and the boundary conditions (global integral conditions[1]vs. local difference approximations). Besides, our result-ant linear equations can be solved in an efficient algorithm. Unlike those papers aforementioned, we apply the present biharmonic solver to study the unsteady incompressible Navier–Stokes flows.
Our method is a FFT-based fast direct solver for the biharmonic equation
(1.4)–(1.6). We first use the truncated Fourier series expansion to derive a set of coupled singular ODEs, then we solve those singular equations by second-order finite difference discretizations. Using a radial grid with shifting a half mesh away from the origin, we can handle the coordinate singularity easily without pole conditions. The Sherman–Morrison formula is then applied to solve the resultant linear system in a cost-efficient way. The computational complexity of the method consists of O(MN log2 N) arithmetic operations
for M· N grid points.
The rest of the paper is organized as follows. In Section 2, we present our fast direct solver for the biharmonic equation(1.4)–(1.6). We then apply this solver to develop a numerical scheme for the unsteady incompressible Na-vier–Stokes flows inside a disk in Section 3. The numerical accuracy check and some test applications have been performed in Section 4.
2. FFT-based fast biharmonic solver 2.1. Fourier mode equations
Since the solution u in Eqs.(1.4) and (1.5)is periodic in h, we can approx-imate it by the truncated Fourier series as
uðr; hÞ ¼ X
N =21
k¼N =2
ukðrÞeikh; ð2:1Þ
where uk(r) is the complex Fourier coefficient given by
ukðrÞ ¼
1 N
XN j¼1
hj= 2jp/N, and N is the number of grid points along a circle. The functions
v(r, h), f(r, h), g(h), and h(h) are defined in the same manner as Eqs.(2.1) and (2.2). The above transformation between the physical space and Fourier space can be efficiently performed using the fast Fourier transform(FFT) with O(N log2N) arithmetic operations.
Substituting those expansions into Eqs. (1.4)–(1.6), we reduce the original PDE to a set of singular ODEs. This common approach is known as the sep-aration of variables in the solution of the linear PDEs. The kth Fourier coef-ficients uk(r) and vk(r) now satisfy the boundary value problems
u00k þu 0 k r k2 r2uk¼ vkðrÞ; ð2:3Þ v00kþv 0 k r k2 r2vk ¼ fkðrÞ; ð2:4Þ ukð1Þ ¼ gk; u0kð1Þ ¼ hk; ð2:5Þ
where the prime denotes the derivative with respect to r. The remaining task is to solve those coupled Fourier mode equations for ukand vkby second-order
finite difference discretizations.
2.2. Spatial discretization and boundary conditions
Throughout this paper, we adapt a radial grid by shifting half mesh width away the origin as
ri¼ ði 1=2ÞDr; i¼ 1; 2; . . . M; M þ 1; ð2:6Þ
where the mesh width Dr = 2/(2M + 1). Under such grid, we have rM+1= 1. The
advantage of this grid is that we do not place grid points directly at the origin; thus, as we shall see, the numerical boundary value near the origin is not needed. This radial grid has been intensively used to develop the efficient Pois-son solvers in 2D polar[9,7]and 3D cylindrical[10]geometries, and to simulate the compressible Navier–Stokes and Euler equations[11].
For convenience of presentation, we simply denote the discrete values Ui uk(ri), Vi vk(ri), and Fi fk(ri). Using second-order centered difference
approximations to discretize Eqs. (2.3) and (2.4), we obtain the difference equations Uiþ1 2Uiþ Ui1 Dr2 þ 1 ri Uiþ1 Ui1 2Dr k2 r2 i Ui¼ Vi; ð2:7Þ
Viþ1 2Viþ Vi1 Dr2 þ 1 ri Viþ1 Vi1 2Dr k2 r2 i Vi¼ Fi ð2:8Þ
for the index 1 6 i 6 M.
In order to close the linear system, the numerical boundary values U0, V0,
UM+1and VM+1should be supplied. Choosing rias described in(2.6), we have
r1= Dr/2. When the index of i = 1, the coefficients of U0and V0in the difference
equations(2.7) and (2.8) equal to zero, respectively; thus, no approximations for U0and V0are actually needed. This is the advantage of using a shifted grid
(2.6). The outer numerical boundary value UM+1 is known by the boundary
value gk. However, there is no direct given value for VM+1. One simple way
to obtain the numerical boundary value VM+1is by the local finite difference
approximation as follows.
First, let us impose Eq.(2.7)at the index i = M + 1; that is, we have VMþ1¼ UMþ2 2UMþ1þ UM Dr2 þ 1 rMþ1 UMþ2 UM 2Dr k2 r2 Mþ1 UMþ1; ð2:9Þ
where UM+2 is the ghost value outside the computational domain.
Approxi-mating the boundary condition u0
kð1Þ ¼ hk by the second-order difference
for-mula, we have UMþ2 UM
2Dr ¼ hk: ð2:10Þ
Therefore, the ghost value UM+2can be obtained by the formula UM+2= UM+
2Drhk. Substituting the value of UM+2, UM+1(=gk), and rM+1= 1 into Eq.(2.9),
we derive the numerical boundary value for VM+1as
VMþ1¼ 2 Dr2UMþ 2 Dr2 k 2 gkþ 2 Drþ 1 hk: ð2:11Þ
This is also explained that the solutions of Uiand Viare coupled through the
boundary conditions.
2.3. Efficient solver for the resultant linear system
After multiplying Dr2 in both sides of Eqs. (2.7) and (2.8), the resultant 2M· 2M linear system has the form
T D E T U V ¼ 0 F : ð2:12Þ
Here, T is an M· M tridiagonal matrix T ¼ d1 s1 l2 d2 s2 lM1 dM1 sM1 lM dM 2 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 5 ; ð2:13Þ
whose entries are di¼ 2 k2 ði 1=2Þ2; si¼ 1 þ 1 2ði 1=2Þ; li¼ 1 1 2ði 1=2Þ ð2:14Þ for 1 6 i 6 M. The matrix D is the M· M diagonal matrix with D = Dr2I. The matrix E has the only one nonzero entry (EM,M= 2sM/Dr2) which
repre-sents the coupling of the boundary condition of V with U. One should notice that this coupling comes from the local approximation(2.11). Solving this sys-tem by the direct Gaussian elimination needs at least O(M2) arithmetic opera-tions. However, one can immediately recognize that the matrix E is as simple as a rank-one matrix. Therefore, we can uncouple the linear equations by ignoring the nonzero entry of the matrix E (that is, to replace E by a zero matrix) and solve the new linear system by inverting the tridiagonal matrix T. In such way, the computational cost has been cut to O(M) operations. This idea is com-pletely described by the Sherman–Morrison formula as follows.
Lemma 1 (Sherman–Morrison formula [15]). If By = b and Bz = a, then (B + abt)x = b has the solution
x¼ y b
t
y
1þ btzz: ð2:15Þ
This formula can be easily verified by direct substitution.
In our implementation of the above formula, the matrices B and abt are chosen as B¼ T D 0 T ; abt¼ 0 0 E 0 : ð2:16Þ
The 2M column vectors a and b are chosen as a¼ ð0; 0; . . . ; 0; 2sM=Dr2Þ
t
where the nonzero entries of a and b are at 2Mth and Mth components, respec-tively. By those choices, we can easily check E = abt. The inversion of the matrix B involves solving the linear system of the tridiagonal matrix T which needs only O(M) operations. One should also notice that the matrix T is the resultant matrix of the Fourier mode equations of Dirichlet Poisson problem which is completely solvable[9].
2.4. Summary of the algorithm
Let us close the section by summarizing the algorithm and the operation counts in the following three steps:
1. Compute the Fourier coefficients for the right-hand side function and the boundary conditions using FFT described in (2.2). This requires O(MN log2N) arithmetic operations.
2. Solve the coupled tridiagonal linear system(2.12) resulting from (2.7) and (2.8)by the Sherman–Morrison formula. This requires O(MN) operations. 3. Convert the Fourier coefficients by inverse FFT(2.1)to obtain the solution,
which requires O(MN log2N) operations.
The overall operation count is thus O(MN log2N) for M· N grid points.
3. Incompressible Navier–Stokes solver on a disk 3.1. Vorticity stream function formulation
The incompressible Navier–Stokes equations has the standard form ou
otþ u ru þ rp ¼ 1
ReDu; ð3:1Þ
r u ¼ 0; ð3:2Þ
where u(x, t) is the fluid velocity, p(x, t) the pressure, and Re is the Reynolds number. The first 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(3.1) and (3.2)by so called the vorticity stream func-tion formulafunc-tion. By taking the curl of Eq.(3.1)to eliminate the pressure gra-dient term, we have
ox
ot þ J ðx; wÞ ¼ 1
ReDx; ð3:3Þ
where x is the nonzero vorticity of the z component, w the stream function, defined by u = ez· $w, and J(x,w) is the Jacobian determinant. Note that,
the velocity u automatically satisfies the incompressibility constraint(3.2). Now using the definition of x, it yields the relation of x and w by
Dw¼ x: ð3:4Þ
Therefore, the original 2D Navier–Stokes equations(3.1) and (3.2)with three primitive variables now has an alternate formulation described by (3.3) and (3.4)with only two unknown variables.
We are interested in the numerical approximations of Eqs. (3.3) and (3.4)
in a unit disk geometry X = {0 <r 6 1, 0 6 h 6 2p}; thus, the polar coordi-nates is used. The nonlinear Jacobian describing the vorticity transport is writ-ten as Jðx; wÞ ¼1 r ow or ox oh ow oh ox or : ð3:5Þ
The radial and azimuthal velocity components can be recovered from the stream function by the formulas
ur¼ 1 r ow oh; uh¼ ow or: ð3:6Þ
The vorticity can be written as x¼ouh or þ uh r 1 r our oh: ð3:7Þ
We restrict our attention to the flow inside a unit disk with some particular velocity specified on the boundary as ur= 0 and uh= h(h) at r = 1. (This
partic-ular boundary condition corresponds to the non-normal flow condition. In addition, if h(h) = 0, the velocity is no-slip at the boundary.) From the relation of(3.6), the above boundary conditions become
wð1; hÞ ¼ 0; ow
orð1; hÞ ¼ hðhÞ: ð3:8Þ
So the complete governing equations include Eqs. (3.3) and (3.4) and the boundary conditions(3.8). Again, one can easily see that there are two bound-ary conditions for the stream function w but no boundbound-ary condition for the vorticity x. This is exactly the same situation as the case of solving the coupled system of biharmonic equations (1.4)–(1.6). It should not be surprising since the 2D Navier–Stokes equations can be actually formulated to a time-depend-ent biharmonic problem of the stream function w. This pure stream function formulation can be obtained by simply substituting the equation (3.4) into
3.2. Time integration
We employ a second-order IMEX (implicit–explicit) backward integration scheme for Eqs.(3.3) and (3.4)as
3xnþ1 4xnþ xn1 2Dt þ ½2J ðx n;wn Þ J ðxn1;wn1Þ ¼ 1 ReDx nþ1; ð3:9Þ Dwnþ1¼ xnþ1; ð3:10Þ wnþ1ð1; hÞ ¼ 0; ow nþ1 or ð1; hÞ ¼ hðhÞ: ð3:11Þ
The superscript on a variable represents the time step index where Dt is the time step. One can easily see that the above time integration scheme has local trun-cation error O(Dt2). Here, we treat the nonlinear convection term explicitly and the linear viscous term implicitly so that at each time step we need to solve a coupled system of Poisson problems just like Eqs. (1.4)–(1.6). Therefore, the fast biharmonic solver described in the previous section can be applied without much modifications.
Recently, the first author has introduced a finite difference scheme for Eqs.
(3.3) and (3.4) which uses the Runge-Kutta method as a time integrator and treats the convection and viscous terms explicitly[8]. The method involves solv-ing a ssolv-ingle Poisson problem for the stream function at each time stage. How-ever, the price to be paid for such simplicity is that the time step has to be chosen very small in order to guarantee the numerical stability. To have a rea-sonable time step size, a Fourier filtering must be implemented to the vorticity near the center at each time stage.
There are other numerical schemes for the unsteady Navier–Stokes equa-tions(3.3) and (3.4)on a disk in the literature. For instance, Torres and Cout-sias[12]have implemented a pseudospectral method with a third-order IMEX backward differencing for time integration to Eqs.(3.3) and (3.4). The pesudo-spectral method involves expanding the vorticity and the stream function in a truncated Chebyshev–Fourier series in r h directions. Therefore, it leads to solve the similar singular Fourier mode equations(2.3) and (2.4)by Chebyshev method. In order to keep the spectral accuracy, some complicated precondi-tioning techniques must be employed for different modes which makes solving the resultant linear equations quite complex. On the other hand, as we dis-cussed before, our resultant linear equations can be solved in a simple and effi-cient way without any preconditioning techniques.
3.3. Spatial discretization
In the computation of the nonlinear convection term, we need to compute the first derivatives of w and x. This can be easily approximated by the
second-order centered difference method. Here, we use the same M· N grid as in the biharmonic problem, that is,
ðri;hjÞ ¼ ðði 1=2ÞDr; jDhÞ; ð3:12Þ
where Dr = 2/(2M + 1) and Dh = 2p/N. Again, by shifting a half mesh in radial direction, we avoid placing grid points directly at the origin.
Let the discrete values of the scalar function w be denoted by wi,j w(ri,hj).
Then, the first derivatives of r and h can be approximated by ow or ij ¼wiþ1;j wi1;j 2Dr ; ð3:13Þ ow oh i;j ¼wi;jþ1 wi;j1 2Dh : ð3:14Þ
Since the function is periodic in h, the approximation of h-derivative does not run into any trouble. However, at i = 1, the numerical boundary value w0,jmust be provided. One appropriate choice of the value is w0;j¼ w1;jþN
2. This
is because if we replacer by r and h by h + p in the Cartesian-polar transfor-mation, the Cartesian coordinates of a point remain fixed. Therefore, any sca-lar function satisfies w(r, h) = w(r, h + p) if the domain of the function is extended to negative values of r. The same spatial discretization is applied to the vorticity x.
3.4. Vorticity boundary conditions
After finishing one time step in our scheme (3.9) and (3.10), we obtain the vorticity xn+1and the stream function wn+1at those interior grid points. How-ever, the boundary vorticity is needed for the approximation of the convection term near the boundary. This can be derived using ThomÕs formula[16]as fol-lows. We first approximate the second boundary condition of (3.11) at the boundary rM+1= 1 by
wnþ1Mþ2;j wM ;jnþ1
2Dr ¼ hj; ð3:15Þ
where wnþ1Mþ2;jis a ghost value outside the computational domain. Therefore, we have wnþ1Mþ2;j¼ wM ;jnþ1þ 2Drhj. Substituting the value of wnþ1Mþ2;jand using the fact
of wnþ1Mþ1;j¼ 0 for all j, we can compute the boundary vorticity by the discrete Laplacian as
xnþ1Mþ1;j¼ Dwnþ1Mþ1;j¼2w
nþ1
M ;j þ 2Drhj
The readers who are interested in the vorticity boundary conditions and the re-lated issues can refer to[3].
4. Numerical results
In this section, we first perform the accuracy check for our numerical schemes to the biharmonic and Navier–Stokes equations on a disk. Then we demonstrate the numerical validness of our Navier–Stokes solver by simulating the moving wall and tripole formation problems. In particular, those problems are picked to test our correct treatments near the boundary (moving wall prob-lem) and near the center (tripole formation).
4.1. Accuracy check for the biharmonic solver
We start our numerical tests by checking the accuracy of our FFT-based fast biharmonic solver on a disk. We simply test two exact solutions of Eq.(1.1)as
u1ðr; hÞ ¼ 1 4ð1 r 2Þð1 þ r cos hÞ; f 1ðr; hÞ ¼ 0; ð4:1Þ u2ðr; hÞ ¼ erðcos hþsin hÞ; f2ðr; hÞ ¼ 4u2ðr; hÞ: ð4:2Þ
The first solution is chosen as the same one used in[6].
Table 1shows the L1errors for our test problems. Here, we fix the number
of grid points in the azimuthal direction as N = 64 and vary the number of grid points M in the radial direction. It is clear that both convergent rates approach two. Therefore, our method is indeed second-order accurate for the biharmonic equation.
4.2. Accuracy check for the Navier–Stokes solver
In this example, we check the accuracy of our scheme for the Navier–Stokes equations. We have taken the exact solution for the Navier–Stokes equations as
Table 1
L1errors for the biharmonic equation
M L1error (u1) Rate L1error (u2) Rate
16 2.4736E04 – 1.2272E03 –
32 6.3760E05 1.96 3.1760E04 1.92
64 1.6186E05 1.98 8.0719E05 1.98
128 4.0786E06 1.99 2.0340E05 1.99
xðx; y; tÞ ¼ 2e2t=Recos x cos y; wðx; y; tÞ ¼ e2t=Recos x cos y: ð4:3Þ
The functions are described in Cartesian coordinates for simplicity of presen-tation. The actual computations are all in polar coordinates.
In our test, we use M· N grid points in the disk so that there are M points in the radial direction and N points in the azimuthal direction. The Reynolds number is Re = 20. The time step is chosen as Dt = 0.01 and the approximate solutions were computed up to T = 2.Table 2shows the L1errors for different
number of grid points. One can easily see that the second-order accuracy has been achieved for both the stream function and the vorticity.
Table 2
L1errors for the Navier–Stokes equations
M· N L1error (w) Rate L1error (x) Rate
16· 32 3.2481E04 – 1.6555E03 – 32· 64 8.3609E05 1.96 6.0361E04 1.46 64· 128 2.1269E05 1.97 1.7393E04 1.80 128· 256 5.3372E06 1.99 4.1738E05 1.88 T=4 T=8
4.3. Moving-wall problem
The moving-wall problem is a flow problem generated by the tangential mo-tion of the boundary of the disk. Thus, the radial boundary velocity is always kept to be ur(1,h) = 0. Here, we choose the same azimuthal (tangential)
bound-ary velocity as in[4]
uh¼ cos h sin h: ð4:4Þ
This problem is chosen to test the numerical treatments near the boundary for our Navier–Stokes solver since the vorticity is generated by the movement of the boundary. One can expect that the maximal vorticity occurs near the boundary. In our run, we use a 128· 128 grid and the time step Dt = 0.01. The compu-tations were computed up to time T = 8.Figs. 1and2show the vorticity con-tours and the streamlines at time T = 4 and 8 for two different Reynolds numbers Re = 100 and 300, respectively. One can see the flow is breaking into four quadrants which is because the azimuthal velocity(4.4)changes sign four times on the boundary. Besides, the gradients of the vorticity contours of the case Re = 300 are steeper than the one observed at Re = 100 which reflects the more diffusion of the vortex in the latter case.
T=4
T=8
4.4. Tripole formation
In order to verify the validness of our numerical treatment near the center, we perform the similar numerical simulation of the tripole formation as in[13]. It has been shown in laboratory and computer experiments[14]that a shielded monopolar vortex, when perturbed, will produce a tripolar vortex. The physics and characteristics of such tripole formation are well explained in [14]. This example serves a perfect test for the coordinate singularity treatment if we place the monopolar vortex at the center.
As in[13], the initial vorticity profile is given by xðrÞ ¼ 1 a 2 r q a eðr=qÞa; ð4:5Þ T=0 T=15 T=30 T=45
Fig. 3. Vorticity contour plots of the tripole formation from T = 0 to 45: Ô-Õ positive values, Ô.Õ negative values.
where the parameter a controls the steepness of the vorticity gradients and q controls the size of the monopole. One can easily see from(4.5), the maximal vorticity occurs at r = 0. In our test, the parameters are chosen as a = 3, q = 1, and the Reynolds number Re = 2000. We extend the radial computational do-main to r = 4 where the no-slip conditions are imposed on the boundary. We use a 128· 128 grid and the time step Dt = 0.005.
To speed up the tripole formation, a random perturbation has been added to the vorticity distribution of(4.5)in the neighboring region where the vortic-ity changes sign. Fig. 3andFig. 4show the formation process of the tripolar vortex from T = 0 to 105. As the time evolves, the initial vortex structure soon breaks the symmetry. The central (positive) vortex becomes elliptical and along the longer sides the negative vorticity organizes into two satellite vortices. This steady configuration then rotates about the center of the positive vortex core in the sense of the positive vorticity. The time evolution of the peak vorticity of
T=60 T=75
T=90 T=105
Fig. 4. Vorticity contour plots of the tripole formation from T = 60 to 105: Ô-Õ positive values, Ô.Õ negative values.
the positive and negative vortices are shown in Fig. 5. One can see the peak vorticity decreases due to the diffusion effect.
Acknowledgement
The authors were supported in part by the National Science Council of Tai-wan under research grant NSC-92-2115-M-009-012.
References
[1] S.C.R. Dennis, M. Ng, P. Nguyen, Numerical solution for the steady motion of a viscous fluid inside a circular boundary using integral conditions, J. Comput. Phys. 108 (1993) 142–152. [2] L.W. Ehrlich, M.M. Gupta, Some difference schemes for the biharmonic equation, SIAM J.
Numer. Anal. 12 (1975) 773–789.
[3] E. Weinan, J.-G. Liu, Vorticity boundary condition and related issues for finite difference schemes, J. Comput. Phys. 124 (1996) 368.
[4] L. Greengard, M.C. Kropinski, An integral equation approach to the incompressible Navier– Stokes equations in two dimensions, SIAM J. Sci. Comput. 20 (1) (1998) 318.
[5] W. Huang, T. Tang, Pseudospectral solutions for steady motion of a viscous fluid inside a circular boundary, Appl. Numer. Math. 33 (2000) 167–173.
0 10 20 30 40 50 60 70 80 90 100 –0.4 –0.2 0 0.2 0.4 0.6 0.8 1
Fig. 5. Time evolution of the peak vorticity for the tripole formation: Ô-Õ peak positive vorticity, Ô.Õ peak negative vorticity.
[6] A. Karageorghis, T. Tang, A spectral domain decomposition approach for steady Navier– Stokes problems in circular geometries, Comput. Fluids 25 (1996) 541–549.
[7] M.-C. Lai, A simple compact fourth-order Poisson solver on polar geometry, J. Comput. Phys. 182 (2002) 337–345.
[8] M.-C. Lai, Fourth-order finite difference scheme for the incompressible Navier–Stokes equations in a disk, Int. J. Numer. Meth. Fluids 42 (2003) 909–922.
[9] M.-C. Lai, W.-C. Wang, Fast direct solvers for Poisson equation on polar and spherical geometries, Numer. Meth. Partial Differ. Eq. 18 (2002) 56–68.
[10] M.-C. Lai, W.-W. Lin, W. Wang, A fast spectral/difference method without pole conditions for Poisson-type equations in cylindrical and spherical geometries, IMA J. Numer. Anal. 22 (2002) 537–548.
[11] K. Mohseni, T. Colonius, Numerical treatment of polar coordinate singularities, J. Comput. Phys. 157 (2000) 787–795.
[12] D.J. Torres, E.A. Coutsias, Pseudospectral solution of the two-dimensional Navier–Stokes equations in a disk, SIAM J. Sci. Comput. 21 (1) (1999) 378–403.
[13] R. Verzicco, P. Orlandi, A finite-difference scheme for three-dimensional incompressible flows in cylindrical coordinates, J. Comput. Phys. 123 (1996) 402–414.
[14] P. Orlandi, G.F. van Heijst, Numerical simulation of tripolar vortices in 2D flow, Fluid Dyn. Res. 9 (1992) 179–206.
[15] J.C. Strikwerda, Finite Difference Scheme and Partial Differential Equations, 1989. [16] A. Thom, The flow past circular cylinders at low speeds, Proc. Roy. Soc. London A 141 (1933)
