Solutions of 2D and 3D Stokes laws using multiquadrics method
D.L. Young*, S.C. Jane, C.Y. Lin, C.L. Chiu, K.C. Chen
Department of Civil Engineering and Hydrotech Research Institute, National Taiwan University, Taipei 10617, Taiwan Received 27 July 2002; revised 16 February 2003; accepted 10 April 2003
Available online 5 March 2004
Abstract
In this paper, velocity – vorticity formulation and the multiquadrics method (MQ) with iterative scheme are used to solve two (2D) and three-dimensional (3D) steady-state incompressible Stokes cavity flows. The method involves solving of Laplace type vorticity equations and Poisson type velocity equations. The solenoidal velocity and vorticity components are obtained by iterative procedures through coupling of velocity and vorticity fields. Both the Poisson type velocity equations and the Laplace type vorticity equations are solved using the MQ, which renders a meshless (or meshfree) solution. Here, the results of 2D Stokes flow problems in a typical square cavity and a circular cavity are presented and compared with other model results. Besides utilizing the MQ to solve the 3D Stokes cubic cavity flow problem, we are also obtaining promising results for the accuracy of the velocity and vorticity. The MQ model has been found to be very simple and powerful for analyzing the 2D and 3D internal Stokes flow problems.
q2004 Elsevier Ltd. All rights reserved.
Keywords: Velocity – vorticity formulation; Stokes flow; Iterative solution; Meshless; Multiquadrics method; Two-dimensional; Square cavity; Circular cavity; Three-dimensional; Cubic cavity
1. Introduction
Stokes flow is a classical problem in fluid dynamics. After many years of development, three formulations for numerical analysis of the Stokes flow are most known: vorticity – stream vector function or vorticity – vector poten-tial approach, primitive variable (velocity – pressure) approach, and vorticity – velocity approach. Using the vorticity – velocity approach, we can transfer the governing equations of primitive variable Stokes equations into a system of Laplace and Poisson equations for the com-ponents of vorticity and velocity fields. By doing so, the direct computation of pressure that is rather difficult to deal with is avoided.
Among various methods for computational fluid dynamics (CFD), the finite difference method (FDM), the finite element method (FEM) and the boundary element method (BEM) are the three most famous numerical schemes adapted for mesh generation. Although the field of computational science was developed widely and quickly, there still exist some unresolved problems for
these numerical methods. The following drawbacks of the conventional numerical schemes with mesh are observed: the awkward treatment of irregular boundary in FDM; the storage of huge data in FEM; the difficulty of the singularities and fundamental solutions in BEM; etc. In this study, an innovative scheme—the multiquadrics method (MQ) or the so-called Kansa’s method [1,2] is used to deal with the above-mentioned shortcomings and also to provide the meshfree scheme.
MQ scheme is a truly scattered, grid free (or meshless) scheme for representing surfaces and bodies in an arbitrary number of dimensions. The radial basis functions (RBF) depend only upon distances between pairs of points. So, it is very powerful for dealing with irregular domain problems. MQ scheme is a very promising method not only for very accurate interpolation of functions, but also for their partial derivatives, divergences, curls, gradients, and integrals. So, error analysis and solving of differential and integral equations become much easier than with the use of the conventional numerical methods. By assuming the suitable interpolation functions—the radial basis functions, there will be no fundamental solutions or singular integration problems involved, such as in the BEM. Madych and Nelson [3] also have proved the theoretical justification for MQ’s performance and gave us confidence to justify our numerical 0955-7997/$ - see front matter q 2004 Elsevier Ltd. All rights reserved.
doi:10.1016/j.enganabound.2003.04.002
www.elsevier.com/locate/enganabound
* Corresponding author. Tel.: þ886-2-2362-6114; fax: þ886-2-3366-5866.
results. Recently, Cheng et al. [4]also have demonstrated the exponential convergence for the MQ collocation method, if the free shape parameter is chosen properly.
MQ was used to simulate the shallow water problems by Hon et al. [5]. They used the partial derivatives of radial basis functions and scattered collocation points to solve the shallow water equation problems in irregular topographic water bodies. As for the size of the region, the resolving way of collocation, and the rank of matrix, Kansa and Hon[6] discussed the ill conditioning of the problem applied to elliptic partial differential equations. They explore several techniques to improve the condition of the coefficient matrix to get rid of the ill conditioning of matrix inversion and achieve the more accurate results. The shape parameter in the radial basis functions is another issue in hot debate nowadays. Kansa and Carlson[7], for example, investigate the usage of variable shape parameters to improve the accuracy of multiquadric interpolation. An adaptive collo-cation method to solve the Burger’s equation that is a nonlinear partial differential equation is proposed by Hon and Mao[8]. More descriptions of the MQ scheme to solve various linear and nonlinear partial differential equations can be referred to Fedoseyev et al.[9].
Sometimes, in order to increase the accuracy of the result, more points may be required to represent the configuration and characteristics of a region. Therefore, solving large radial basis function interpolation problems encounter some difficulty in terms of computational cost and magnitude of storage and operation. Consequently, Beatson et al.[10]propose a way to resolve this problem. They adopt suitable approximate basis function with precondition and employ the GMRES iterative method to resolve these problems.
For large matrix vector problem, Beatson and Newsam [11 – 13]proposed to use the domain decomposition method to divide the domain into a few subdomains. The techniques are based on the hierarchical and multipole expansion for the calculation of many-body potentials. This method provides the feature for fast, storage-efficient computation of the matrix vector product.
Mei-Duy and Tran-Cong[14]used the MQ RBF to solve the 2D Navier – Stokes equations including planar Poi-seuille, driven cavity and natural convection flows. They claimed that employment of the indirect radial basis function networks method (IRBFN) [15] would render better results with relatively low collocation densities. However, since the streamfunction and vorticity formu-lation are used, only 2D flow problems are solved. Extension to 3D problems would be very difficult.
In this paper, we combine the vorticity – velocity formulation with the MQ scheme to solve steady-state, both 2D and 3D Stokes flow problems. Three numerical examples, a square cavity problem and a circular cavity problem for 2D and a cubic cavity flow problem for 3D are employed to check the validity and efficiency of the MQ scheme. As far as the more complex Navier – Stokes flow
problems are concerned, this paper only provides a preliminary work. MQ scheme really offers an efficient and powerful as well as meshless numerical method to the realm of CFD and deserves further research and application.
2. Governing equations
The governing equations of Stokes flows for the velocity – vorticity formulation in steady-state can be derived from the Navier – Stokes equations by neglecting the inertia and are given as follows
72v~¼ ~0 ð1Þ
72~u ¼ 27 £ ~v ð2Þ
Here ~u is the velocity vector, and ~vis the vorticity vector. The vorticity vectorv~is defined by the following
~
v¼ 7 £~u ð3Þ
In two-dimensions, if ðu; vÞ are the horizontal velocity vectors and v~k is the associated vertical vorticity vector,
then the vorticity transport Eq. (1) can be expressed as
72v¼ 0 ð4Þ
Eq. (2) that governs the u; v-velocity becomes as follows 72u ¼ 2›v
›y ð5Þ
72v ¼ ›v
›x ð6Þ
In steady-state conditions, solution of Eq. (4), in conjunc-tion with the Poisson Eqs. (5) and (6), gives the velocity and vorticity distributions all over the domain.
In three-dimensions, if ~u ¼ ðu; v; wÞ are the velocity vectors andv~¼ ðj;h;zÞ are the associated vorticity vectors, then the individual vorticity transport Eq. (1) can be described as
72j¼ 0 ð7Þ
72h¼ 0 ð8Þ
72z¼ 0 ð9Þ
Similarly, Eq. (2) can be deduced for 3D velocity and illustrated by the following
72u ¼ 2 ›z ›y 2 ›h ›z ð10Þ 72v ¼ 2 ›j ›z 2 ›z ›x ð11Þ 72w ¼ 2 ›h ›x 2 ›j ›y ð12Þ Similarly, being in steady-state conditions, solution of Eqs. (7) – (9), in conjunction with the Poisson Eqs. (10) – (12),
gives the velocity and vorticity distributions all over the domain in 3D Stokes flow problems.
A solution is sought in the domain V that satisfies the boundary conditions, ~u ¼ ~uG ð13Þ and ~ v¼ ð7 £~uÞlG ð14Þ on the boundary. 3. Numerical formulation
3.1. Two-dimensional numerical formulation Consider the 2D Poisson’s equation
72F¼ f ðx; yÞ ðx; yÞ [ domain ð15Þ
F¼ gðx; yÞ ðx; yÞ [ boundary ð16Þ
Let F¼X m i¼1 ai ffiffiffiffiffiffiffiffiffi r2 i þ c2 q ð17Þ Where riis the Euclidean distance between the field point
ðx; yÞ and collocation point ðxi; yiÞ; and c is a free shape
parameter to be determined later. In general, 0:1 # c # 1:0 is recommended in our study. The determination of c is still an attractive research topic.ai’s are the coefficients to be
determined later by the point collocation method. f ðx; yÞ and gðx; yÞ are given functions. In this study,F¼ u; or v; or v;
and f ðx; yÞ ¼ 2ð›v=›yÞ; or ð›v=›xÞ; or 0. We choose m collocation points {ðxi; yiÞ}m1 and use Eq. (5) for the
u-velocity. The other unknowns, v andvcan be obtained in a similar fashion. Let u ¼X m i¼1 aui ffiffiffiffiffiffiffiffiffi r2 i þ c2 q ð18Þ v¼X m i¼1 avi ffiffiffiffiffiffiffiffiffi r2 i þ c2 q ð19Þ ›2 ›x2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 q ¼ ðy 2 yiÞ 2þ c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 p 3 ð20Þ Similarly, ›2 ›y2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 q ¼ ðx 2 xiÞ 2þ c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 p 3 ð21Þ
For the vorticity derivative ofv; › ›y Xm i¼1 avi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 q ¼X m i¼1 avi ðy 2 yiÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 p ð22Þ
Substituting Eqs. (20) – (22) into Eq. (5), we obtain Xm i¼1 aui ðy 2 yiÞ 2þ c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 p 3 0 B @ 1 C A þX m i¼1 aui ðx 2 xiÞ 2þ c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 p 3 ¼X m i¼1 avi 2ðy 2 yiÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 p Also Xm i¼1 aui ðx 2 xiÞ 2þ ðy 2 y iÞ2þ 2c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 p 3 0 B @ 1 C A ¼X m i¼1 avi 2ðy 2 yiÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 p
By collocation method, we obtain n equations from the interior points and m 2 n boundary points.
Xm i¼1 aui ðxj2 xiÞ 2þ ðy j2 yiÞ2þ 2c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxj2 xiÞ2þ ðyj2 yiÞ2þ c2 q 3 0 B @ 1 C A ¼X m i¼1 avi 2ðyj2 yiÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxj2 xÞ2þ ðyj2 yiÞ2þ c2 q j ¼ 1; …; n ð23Þ
We choose ðm 2 nÞ boundary points {ðxi; yiÞ}mnþ1: From Eq.
(16), we obtain Xm i¼1 aui ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðxj2 xiÞ2þ ðyj2 yiÞ2þ c2 q ¼ gðxj; yjÞ j ¼ n þ 1; …; m ð24Þ
Solving m £ m system of Eqs. (23) and (24) by iteration, we obtain the coefficients au
i: After that, u is found from
Eq. (18) for all points located in the interior domain or on the boundary.
3.2. Three-dimensional numerical formulation
In a similar way, we consider the 3D Poisson’s equation 72Fk¼ fkðx; y; zÞ ðx; y; zÞ [ domain ð25Þ
Fk¼ gkðx; y; zÞ ðx; y; zÞ [ boundary ð26Þ Let Fk¼ Xm i¼1 ai ffiffiffiffiffiffiffiffiffi r2 i þ c2 q ð27Þ where riis the Euclidean distance between the space point
ðx; y; zÞ and collocation point ðxi; yi; ziÞ; and c is a free shape
parameter to be determined. fkðx; y; zÞ is equal to 27 £ ~vor
~0; andFkis equal to~u or ~v: We choose m collocation points
{ðxi; yi; ziÞ}m1: Now, we will illustrate Eq. (10) for u-velocity
only. The other unknowns, v; w;j;h; andzcan be obtained in a similar way. Let u ¼X m i¼1 aui ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 q ð28Þ z¼X m i¼1 azi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 q ð29Þ ›2 ›x2 Xm i¼1 aui ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 q ¼X m i¼1 aui ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ ðz 2 ziÞ þ c2 p 3 ð30Þ Similarly, ›2 ›y2 Xm i¼1 aui ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 q ¼X m i¼1 aui ðx 2 xiÞ2þ ðz 2 ziÞ2þ c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 p 3 ð31Þ ›2 ›z2 Xm i¼1 aui ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 q ¼X m i¼1 aui ðx 2 xiÞ2þ ðy 2 yiÞ2þ c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 p 3 ð32Þ
For vorticity derivative ofz; › ›y Xm i¼1 azi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ22 ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 q ¼X m i¼1 azi y 2 yi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 p ð33Þ Similarly › ›z Xm i¼1 ahi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ22 ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 q ¼X m i¼1 ahi z 2 zi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx 2 xiÞ2þ ðy 2 yiÞ2þ ðz 2 ziÞ2þ c2 p ð34Þ
Substituting Eqs. (30) – (34) into Eq. (10), and using the collocation method, we obtain n equations from the internal points and ðm 2 nÞ equations from the boundary points. For the – internal points,
Xm i¼1 aui 2ðxj2 xiÞ2þ 2ðyj2 yiÞ2þ 2ðzj2 ziÞ2þ 3c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxj2 xiÞ2þ ðyj2 yiÞ2þ ðzj2 ziÞ2þ c2 q 3 0 B @ 1 C A ¼ 2 X m i¼1 azi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðyj2 yiÞ ðxj2 xiÞ2þ ðyj2 yiÞ2þ ðzj2 ziÞ2þ c2 q 2 6 4 2X m i¼1 ahi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðzj2 ziÞ ðxj2 xiÞ2þ ðyj2 yiÞ2þ ðzj2 ziÞ2þ c2 q 3 7 5 j ¼ 1;…;n ð35Þ
Similarly for the v and w components, we find Xm i¼1 avi 2ðxj2 xiÞ2þ 2ðyj2 yiÞ2þ 2ðzj2 ziÞ2þ 3c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxj2 xiÞ2þ ðyj2 yiÞ2þ ðzj2 ziÞ2þ c2 q 3 0 B @ 1 C A ¼ 2 X m i¼1 aji ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðzj2 ziÞ ðxj2 xiÞ2þ ðyj2 yiÞ2þ ðzj2 ziÞ2þ c2 q 2 6 4 2X m i¼1 azi ðxj2 xiÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxj2 xiÞ2þ ðyj2 yiÞ2þ ðzj2 ziÞ2þ c2 q 3 7 5 j ¼ 1;…;n ð36Þ Xm i¼1 awi 2ðxj2 xiÞ2þ 2ðyj2 yiÞ2þ 2ðzj2 ziÞ2þ 3c2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxj2 xiÞ2þ ðyj2 yiÞ2þ ðzj2 ziÞ2þ c2 q 3 0 B @ 1 C A ¼ 2 X m i¼1 ahi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðxj2 xiÞ ðxj2 xiÞ2þ ðyj2 yiÞ2þ ðzj2 ziÞ2þ c2 q 2 6 4 2X m i¼1 aji ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðyj2 yiÞ ðxj2 xiÞ2þ ðyj2 yiÞ2þ ðzj2 ziÞ2þ c2 q 3 7 5 j ¼ 1;…;n ð37Þ
We choose ðm 2 nÞ boundary points {ðxi;yi;ziÞ}mnþ1: From
Eq. (26), we obtain Xm i¼1 aki ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxj2 xiÞ2þ ðyj2 yiÞ2þ ðzj2 ziÞ2þ c2 q ¼ gkðxj;yj;zjÞ j ¼ n þ 1;…;m; k ¼ u;v;w ð38Þ
Solving the m £ m system of Eqs. (35) – (38) by iteration, we obtainak
i; the velocity coefficients. Similarly, we use the
equation, and then get the coefficient an
i of the vorticity,
where n is chosen fromj;h; andz:
It is noticed that in order to reduce the size of the coupling of velocity and vorticity fields, we choose the iterative scheme to solve the velocity and vorticity in sequence. The vorticity provides the internal source for the velocity; while the solution of vorticity needs the boundary condition for vorticity from which it is obtained by taking the curl of the velocity field. The details are elaborated in Section 4.
4. Solution procedures
The MQ procedure is applied to solve the Stokes flows via the velocity – vorticity formulation. In most incompres-sible viscous flow problems, the essential boundary conditions are the prescribed velocities. The vorticity boundary conditions are not known a priori in general, and should be determined iteratively from the computation of the velocity field. In the present model, the velocity Poisson equations are firstly solved to obtain the velocity and the associated vorticity boundary conditions. After that, the vorticity is obtained by solving the vorticity Laplace equation. The solution procedures are described in the following steps
1. Solve the velocity Poisson equations using the MQ (Kansa’s) method by taking known vorticity components.
2. Calculate the unknown boundary velocity or velocity flux values.
3. Calculate the velocity distribution and velocity deriva-tives at all nodal points.
4. Determine the new vorticity boundary values using the definition of vorticity.
5. Solve the vorticity Laplace equation using the MQ (Kansa’s) method.
6. Calculate the unknown vorticity values throughout the domain.
7. Calculate the derivatives of the vorticity vectors to be used in the velocity Poisson equations.
8. Check for the convergence of the velocity and vorticity components in the present iteration, for example,
ukþ1i 2 u k i
# 1024
9. If convergence criterion is satisfied, then stop. Otherwise replace the old vorticity by a new one and then go to step 1.
The equations of the full matrix are then solved by using the Gaussian elimination technique. We have encountered no ill-conditioning problem of the full matrix as big as the size of 3375 nodes without the introduction of precondition
or domain decomposition methods as commonly adopted in the literature.
5. Numerical examples
The present MQ scheme has been applied to test three problems to verify the feasibility and accuracy of the method. The test problems are the classical ‘driven flow in a square cavity’ and ‘flow in a circular cavity’ problems for which many analytical and numerical model results are available in the literature for 2D Stokes flows. A cubic cavity flow on the other hand is employed to illustrate the feasibility of the 3D Stokes flows.
Example 1. The first model problem consists of a 2D square cavity filled with incompressible viscous fluid moving at a constant velocity on the top lid. No slip or impervious boundary conditions are imposed on any wall, with the velocity at the upper wall set to unity. For the present computation, 33 boundary points and 88 internal points are used. In this example, we use 11 £ 11 nodes to collocate the region, (1,0) £ (0,1). The compu-tational domain is divided into two regions, boundary and interior regions. The value of the shape function c is taken as 0.39.
The square cavity flow with configuration and boundary conditions is illustrated in Fig. 1. The velocity vector is shown inFig. 2, and the distribution of vorticity is depicted inFig. 3, respectively. We use two sets of collocation points for 11 £ 11 nodes and 13 £ 13 nodes respectively to express the x-component velocity profile ðuÞ on the vertical center-line of the cavity, and compared with the accuracy of the two collocation sets. The diagram ofFig. 4reveals almost no difference between the two collocations for the velocity
profiles. Fig. 5 displays the velocity profile ðvÞ of the y-component on the horizontal centerline for steady-state Stokes flow by 11 £ 11 collocation points, and Fig. 6 exhibits the velocity profile ðuÞ of the x-component on the vertical centerline of the cavity. The results are compared with Burggraf series solution[16], a FEM solution[17], a BEM solution [18], and a DRBEM (direct method) solution [19]. Good agreement is observed even in such small sets of collocation points. As can be seen inFig. 6, the present MQ (Kansa’s) method is as accurate as the other numerical methods, even if very few collocation points are used.
Example 2. The second model problem consists of a recirculating flow in a 2D circular cavity[16,20]. The radius of the circular cavity is assumed to be unity. The configuration and boundary conditions of the problem is depicted in Fig. 7. In the upper half of the boundary, the velocity uu¼ 1 in a contour clockwise sense and ur¼ 0
are prescribed, and in the lower half, uu¼ ur¼ 0 are
imposed. For this computation, a total of 612 boundary points and internal points are used in this study.
Fig. 2. Velocity vectors for Stokes flow in a square cavity.
Fig. 3. Vorticity distribution for Stokes flow in a square cavity.
Fig. 4. Velocity profile ðuÞ on the centerline at x ¼ 0:5 of a square cavity; comparison of mesh convergence.
The velocity vector is observed in Fig. 8, and the distribution of vorticity contour is shown in Fig. 9. The present computations give very good results as compared with the other solutions.Fig. 10depicted the x-component velocity profile ðuÞ on the vertical centerline, and Fig. 11 describes the y-component velocity profile ðvÞ on the horizontal centerline. The results are compared with Burggraf series solution [16] and Hwu et al. analytical solution[20], and both also show good agreements.
Example 3.The third model problem is a 3D cubic cavity, with the top plate moving with unit velocity in the x-direction ðz ¼ 1; 0 # x # 1; 0 # y # 1; u ¼ 1Þ: Others are set with no-slip and impervious Dirichlet conditions.
The configuration and boundary conditions of the problem are shown inFig. 12. We used regular collocation points in the cube for three different uniform meshes of size 11 £ 11 £ 11, 13 £ 13 £ 13, and 15 £ 15 £ 15 in this simulation. A mesh sensitivity study using the three meshes of size 11 £ 11 £ 11, 13 £ 13 £ 13, and 15 £ 15 £ 15 indicates that the mesh 15 £ 15 £ 15 predicts velocity and vorticity with the acceptable error.
The velocity vectors which take x-component velocity u and z-component velocity w are shown in Fig. 13. The distribution of vorticity which takes the x – z plane at y ¼ 0:5 is depicted inFig. 14. From another aspect, velocity vector Fig. 6. Velocity profile ðuÞ on the centerline at x ¼ 0:5 of a square cavity.
Fig. 7. Circular cavity flow problem with boundary conditions.
Fig. 8. Velocity vectors for Stokes flow in a circular cavity.
at x – y plane at z ¼ 0:5 for Stokes flow in the cubic cavity is displayed inFig. 15. The velocity affects the distribution of the intensity of the vorticity, so we look more closely at the vorticity profile for more cross sections. Fig. 16 expresses the jvorticity intensity which is located at z ¼ 0:5 of x–y plane; the results indicate the symmetric characteristics with respective to x- and y-directions. The effect of the wall makes vorticity distribution more concentrated on the edges of the cavity. This consequence is conspicuously revealed inFigs. 17 and 18, respectively. Fig. 19 shows the x-component velocity profile ðuÞ on the vertical centerline of the cavity, and Fig. 20 shows the z-component velocity profile ðwÞ on the horizontal
centerline for steady-state Stokes flow of the cubic cavity. The results are compared with those obtained using meshless BEM method (Tsai et al. [21]) and the results of a traditional BEM – FEM (Young et al.[22]). The result of our present method with 11 £ 11 £ 11 mesh is not as good as the results obtained by the meshless BEM (Tsai et al. [21]). However, it is conjectured that adding more point density, for example by using meshes of size 13 £ 13 £ 13 or 15 £ 15 £ 15, the quality of our results could be improved, and also can capture the effect of wall at driven direction with increased node density. Nevertheless, the MQ (Kansa’s) scheme is very easy to implement as compared to other numerical methods.
Fig. 10. Comparison of u-velocity profile at x ¼ 0 for a circular cavity.
Fig. 11. Comparison of v-velocity profile at y ¼ 0 for a circular cavity.
Fig. 12. Cubic cavity flow problem with boundary conditions.
Fig. 13. Velocity vectors of x – z plane at y ¼ 0:5 for Stokes flow in a cubic cavity.
In order to make more quantitative comparisons with analytic or other numerical solutions, the following criterion of root mean squared (RMS) error analysis is established. The RMS errors of x-component velocity profile ðuÞ on the vertical centerline, as well as y-component velocity profile ðvÞ on the horizontal centerline are both considered. The base used for the square cavity is the Burggraf series solution [16], while the analytic solution of Hwu et al. [20] is employed for the circular cavity. The RMS errors for 2D square cavity and circular
cavity are shown in Table 1. The magnitudes are in the range of 8 £ 1023 to 4 £ 1024. For the 3D cubic cavity flow problem, the base used for calculating RMS is Tsai et al.[21]meshless results since there are no analytic solutions available. The corresponding RMS error is illustrated in Table 2. The values lie between 6 £ 1022 and 8 £ 1028. In general, the present MQ (Kansa’s) method renders reasonable and comparable accuracy as compared with exact and other numerical schemes. Fig. 14.h-component vorticity distribution of x – z plane at y ¼ 0:5 for
Stokes flow in a cubic cavity.
Fig. 15. Velocity vectors of x – y plane at z ¼ 0:5 for Stokes flow in a cubic cavity.
Fig. 16. j-component vorticity distribution of x – y plane at z ¼ 0:5 for Stokes flow in cubic cavity.
Fig. 17.h-component vorticity distribution of y – z plane at x ¼ 0:5 for Stokes flow in cubic cavity.
6. Conclusions
Stokes flows are investigated in 2D square and 2D circular cavity and a 3D cubic cavity. In comparing with analytical solutions and other numerical solutions, the results show that the present scheme is efficient and accurate for using very few collocation points with the merits of meshfree and adaptive mesh advantages. In this investigation, the velocity – vorticity formulation is employed. It mainly deals with the coupled calculation of the vorticity and velocity equations. It means that the accuracy of both results depends on the mesh size of each
other, so we can increase the number of points at some high gradient places, such as near the singular points to raise the resolution of the vorticity accuracy.
The demonstrated good performance of the MQ (Kansa’s) scheme shows that it is a powerful tool for the numerical solution of incompressible viscous fluid flows. Although 2D and 3D Stokes flow problems are primarily solved, we expect the MQ (Kansa’s) scheme will have great potential to solve other 2D and 3D incompressible viscous flow problems that are currently under study.
The main limitation of the MQ scheme is the huge storage of data coming from the operator matrix. If the numerical model is extended to 3D Navier – Stokes problems, the storage of data would become very much larger than that of the 2D model. Therefore, it deserves further research to overcome this problem. Finally, the MQ Fig. 18.z-component vorticity distribution of y – z plane at x ¼ 0:5 for
Stokes flow in cubic cavity.
Fig. 19. Velocity profile ðuÞ on the centerline at x ¼ 0:5 of a cubic cavity.
Fig. 20. Velocity profile ðwÞ on the centerline at y ¼ 0:5 of a cubic cavity. Table 1
The RMS error for square and circular cavities
Shape Variable
Nodes u v
Square 121 0.00299 0.000492
Circular 612 0.00882 0.00116
Table 2
The RMS error for cubic cavity
Mesh size Variable
u w
11 £ 11 £ 11 0.0583 0.00698
13 £ 13 £ 13 0.0177 0.00643
(Kansa’s) is worthy to investigate the more challenging 2D and 3D Navier – Stokes equations especially for the higher Reynolds number flows.
Acknowledgements
This work was supported by the National Science Council, Taiwan. It is greatly appreciated. Also, we would like to express our sincere thanks to Professor A.H.-D. Cheng for his valuable discussions of this paper. Finally, we would like to thank Professor A. S. Muleshkov for his helpful suggestions and revision of this manuscript.
References
[1] Kansa EJ. Multiquadrics—a scattered data approximation scheme with application to computational fluid-dynamics. I. Surface approxi-mations and partial derivative estimates. Comput Math Appl 1990; 19(8/9):127 – 45.
[2] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics.II. Solutions to parabolic, hyperbolic and elliptic partial differential equations. Comput Math Appl 1990;19(8/9):147 – 61.
[3] Madych WR, Nelson S. Bounds on multivariate polynomials and exponential error estimates for multiquadric interpolation. J. Approx Theory 1992; 70:94 – 114.
[4] Cheng AH-D, Golberg MA, Kansa EJ, Zammito G. Exponential convergence and H-c multiquadric collocation method for partial differential equations. Numer Meth Part Differ Eq. 2003;19(5):571– 94. [5] Hon YC, Cheung KF, Mao XZ, Kansa EJ. Multiquadric solution for
shallow water equations. J Hydraul Engng 1999;125(5):524 – 33. [6] Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with
multiquadric radial basis functions: applications to elliptic partial differential equations. Comput Math Appl 2000;39:123– 37. [7] Kansa EJ, Carlson RE. Improved accuracy of multiquadric
interp-olation using variable shape parameters. Comput Math Appl 1992; 24(12):99– 120.
[8] Hon YC, Mao XZ. An efficient numerical scheme for Burger’s equation. Appl Math Comput 1998;95:37 – 50.
[9] Fedoseyev AI, Friedman MJ, Kansa EJ. Continuation for nonlinear elliptic partial differential equations discretized by the multiquadrics method. Int J Bifurcation Chaos 2000;10(2):481 – 92.
[10] Beatson RK, Cherrie JB, Mouat CT. Fast fitting of radial basis functions: method based on preconditioned GMRES iteration. Adv Comput Math 1999;11:253– 70.
[11] Beatson RK, Newsam GN. Fast evaluation of radial basis functions— I. Comput Math Appl 1992;24(12):7 – 19.
[12] Beatson RK, Newsam GN. Fast evaluation of radial basis functions: moment-based methods. SIAM J Sci Comput 1998;19(5):1428 – 49. [13] Beatson RK, Newsam GN. Fast solution of the radial basis function
interpolation equations: domain decomposition methods. SIAM J Sci Comput 2000;22(5):1717 – 40.
[14] Mai-Duy N, Tran-Cong T. Numerical solution of Navier – Stokes equations using multiquadrics radial basis function networks. Int J Numer Meth Fluids 2001;37:65– 86.
[15] Mai-Duy N, Tran-Cong T. Numerical solution of differential equations using multiquadrics radial basis function networks. Neural Netw 2001;14:185 – 99.
[16] Burggraf OR. Analytical and numerical studies of the structure of steady separated flows. J Fluid Mech 1966;124:113– 51.
[17] Young DL, Lin QH. Application of finite element method to 2-D flows. Proceedings of the Third National Conference On Hydraulic Engineering, Taipei, Taiwan; 1986. p. 223 – 42.
[18] Young DL, Yang SK. Simulation of two-dimensional steady Stokes flows by the velocity – vorticity boundary element method. Proceed-ings of the 20th National Conference of Mechanics, Taipei, Taiwan; 1996. p. 422 – 9.
[19] Eldho TI, Young DL. Solution of the velocity – vorticity Navier – Stokes equations using dual reciprocity boundary element method. Proceedings of the Fourteenth Engineering Mechanics Conference, ASCE, Texas, USA; 2000.
[20] Hwu TY, Young DL, Chen YY. Chaotic advections for Stokes flows in circular cavity. J Engng Mech 1997;123:774– 82.
[21] Tsai CC, Young DL, Cheng AH-D. Meshless BEM for three-dimensional Stokes flows. CMES: Comput Model Engng Sci 2002; 3(1):117– 28.
[22] Young DL, Liu YH, Eldho TI. Three-dimensional Stokes flow solution using combined boundary element and finite element methods. Chin J Mech 1999;15:169 – 76.
