Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 www.elsevier.com/locate/cma Simulation of incompressible viscous ﬂows past a circular cylinder by hybrid FD scheme and meshless least square-based ﬁnite diﬀerence method H. Ding, C. Shu *, K.S. Yeo a, D. Xu a b Department of Mechanical Engineering, National University of Singapore 10 Kent Ridge Crescent, Singapore 117576, Singapore b Institute of High Performance Computing, 1 Science Park Road, #01-01 The Capricorn Singapore 117528, Singapore Received 5 May 2003; received in revised form 10 October 2003; accepted 5 November 2003 Abstract In this paper, we describe an eﬃcient method for simulating the two-dimensional steady and unsteady incompressible ﬂows. The method is a hybrid approach, which combines the conventional ﬁnite diﬀerence (FD) scheme and the meshfree least square-based ﬁnite diﬀerence (MLSFD) method. In this method, the MLSFD method is adopted to deal with the complex geometry for its ‘‘truly’’ mesh-free property, and it is responsible for the spatial discretization in the region around the complex geometry which represents the stationary solid obstacle. Correspondingly, conventional FD scheme is applied in the rest of the ﬂow domain to take advantage of its high computational eﬃciency. Numerical simulation is carried out for the steady ﬂow with Reynolds numbers of 10, 20 and 40, and unsteady ﬂow with Reynolds numbers of 100 and 200. The obtained numerical results agree very well with computational and experimental data available in the literature. Compared with the fully MLSFD method, this new approach greatly improves the computational eﬃciency. Ó 2003 Elsevier B.V. All rights reserved. Keywords: Incompressible viscous ﬂows; Hybrid FD scheme; MLSFD 1. Introduction In the past decades, a group of so-called mesh-free or meshless methods have become one of the hottest areas of research in computational ﬂuid dynamics [1–10]. As their name implies, one common characteristic of all these methods is that they do not require the traditional mesh to construct the numerical formulation. Mesh-free methods possess a number of interesting properties. For example, they require node generation instead of mesh generation. In other words, there is no pre-speciﬁed connectivity or relationships among the nodes, thus the computational costs associated with mesh generation are highly reduced. Another attractive property of mesh-free methods is the computational ease of adding and subtracting nodes from the pre-existing nodes. The computational advantages of mesh-free method suggest that they have * Corresponding author. E-mail address: [email protected] (C. Shu). 0045-7825/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.11.002 728 H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 potentials in solving a broad class of scientiﬁc and engineering problems. Despite of their preliminary successes, they do have some drawbacks, i.e., in general, they are less eﬃcient than the conventional numerical methods such as ﬁnite diﬀerence (FD), ﬁnite element (FE) and ﬁnite volume (FV) methods. In the numerical formulation of mesh-free methods, the approximation of the unknown function or its derivatives is constructed on a set of nodes within the local support. Usually, the local support has a circle shape for the two-dimensional case and a sphere shape for the three-dimensional case. Compared with the ﬁnite diﬀerence method on the Cartesian mesh, the construction of unknown function approximation requires much more nodes by mesh-free methods to achieve the same order of accuracy. As a result, the bandwidth of the matrix derived from the resultant algebraic equations is greatly expanded. Therefore, the computational eﬃciency is unavoidably impaired. In this regard, the mesh-free methods achieve the geometry ﬂexibility at the cost of computational eﬃciency. Despite of recent advances in the ﬁelds of meshfree methodology, it is certainly the bottleneck in the industrial applications. To handle ﬂows with complex geometry eﬀectively and eﬃciently is always a tough task. Traditional FD/ FV methods enjoy the computational eﬃciency, but they confront diﬃculties in mesh generation, when dealing with the complex geometry, especially when three-dimensional problems are considered. To improve the drawbacks for the traditional FD/FV methods, Cartesian mesh method was proposed to solve the ﬂow problems in irregular regions [11–13]. This method uses the Cartesian mesh, and the irregular geometries are treated as a specialized boundary embedded in the mesh. Since the computational cost used for the generation of Cartesian mesh is negligible as compared with that required by body-ﬁtted or unstructured mesh, this method is thought to be more eﬃcient. However, in order to ensure the accuracy of the solution by Cartesian mesh method, some special treatments are required for the spatial discretization in the vicinity of the embedded boundary. Another diﬃculty is the time step restriction arising from the small irregular cut cells. This may be the reason why the main applications of Cartesian mesh methods up to date mainly focus on the compressible ﬂow problems [12,13], in which slip boundary condition is applied. In the applications of the Cartesian mesh methods, only one single type of mesh (Cartesian mesh) is used to cover the irregular computational domains. On the other hand, some other researchers suggested using the composite mesh formed by a set of regular grids [14,15]. For example, cylindrical, spherical or nonorthogonal grids are used to ﬁt the geometry of the obstacle bodies while the Cartesian grids are adopted in the rest of the solution domain. This type of composite mesh provides a possible way for the traditional numerical methods like ﬁnite diﬀerence methods, to deal with ﬂow problems of complex geometry. The main drawback of using this kind of mesh is that the information exchange between diﬀerent grids is based on the interpolation. If the solution exhibits intensive variation near the interface between grids, the interpolation process may introduce signiﬁcant numerical errors or cause the break-down of the computation. To improve the above methods, a hybrid method, which combines the conventional ﬁnite diﬀerence scheme and the mesh-free least square-based ﬁnite diﬀerence (FD/MLSFD) approach, is proposed. In the present study, the MLSFD method is adopted for the spatial discretization in the region around the complex geometry. Accordingly, conventional FD scheme is applied in the rest of the ﬂow domain to take advantage of its high computational eﬃciency. The present method is validated by simulating the ﬂow ﬁeld around a circular cylinder at various Reynolds numbers. The simulated ﬂow patterns include both the steady and unsteady state. Numerical experiments are also carried out to examine the eﬃciency improvement using present method as compared with that using MLSFD method solely. 2. MLSFD method and its combination with conventional FD scheme In this section, we ﬁrstly give a brief review on the ‘‘truly’’ mesh-free least square-based ﬁnite diﬀerence (MLSFD) method. Then, we show how the MLSFD and conventional FD schemes are combined to solve the partial diﬀerential equations. H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 729 As its name implies, MLSFD was inspired from the development of traditional ﬁnite diﬀerence method. Both methods were originated from the Taylor series expansion. However, MLSFD uses multi-dimensional Taylor series expansions to construct the approximation of the derivatives of unknown functions while conventional FD scheme employs the one-dimensional Taylor series expansion. As shown in Fig. 1, the functional value at the supporting knot can be approximated by the functional value and its derivatives at the reference knots by using the 2D Taylor series expansion, 2 2 2 o/ o/ 1 1 o/ 2 o / 2 o / / ¼ /0 þ Dx þ Dy þ ðDxÞ þ ðDyÞ þ ðDxÞðDyÞ ox 0 oy 0 2 ox2 0 2 oy 2 0 ox oy 0 3 3 3 3 1 1 1 o/ 1 o/ 3 o / 3 o / 2 2 ðDyÞ ðDxÞ ðDxÞðDyÞ þ ðDxÞ þ þ ðDyÞ þ þ 3 3 2 6 ox 0 6 oy 0 2 ox oy 0 2 ox oy 2 0 ð1Þ Suppose that Eq. (1) is truncated to the third-order derivatives. Then it has nine unknowns. Among them, two are the ﬁrst-order derivatives, three are the second-order derivatives, and four are the third-order derivatives. Application of Eq. (1) at nine neighbouring nodes leads to the following matrix form: Du ¼ Sdu ð2Þ where DuT ¼ ½ /1 /0 ; . . . ; /9 /0 T du ¼ o/ ox 2 2 2 3 3 3 3 o/ o/ o/ o/ o/ o/ o/ o/ ; ; ; ; ; ; ; ; 2 2 3 3 2 oy 0 ox 0 oy 0 ox oy 0 ox 0 oy 0 ox oy 0 ox oy 2 0 0 ST ¼ b s1 ; . . . ; s9 c99 sTj ¼ 2 2 3 3 2 2 Dxj ; Dyj ; 12ðDxj Þ ; 12ðDyj Þ ; ðDxj ÞðDyj Þ; 16ðDxj Þ ;16ðDyj Þ ; 12ðDxj Þ ðDyÞj ; 12ðDxj ÞðDyj Þ The square matrix S contains all the geometric information about the distribution of the supporting nodes. Studying the structure of matrix S, it can be found that the pattern of local distribution of the supporting nodes will determine whether the matrix is singular or ill-conditioned for inversion. This difﬁculty is not easy to be resolved because little is known about the eﬀects of node distribution on the conditioning of the matrix. For numerical implementation, except for a few special cases, such as the case where all the supporting nodes are located on the same straight line, we would hence have to check and ensure that the matrix S is well-conditioned at every node in the computational domain. This can be done Reference knot Supporting knots Non-supporting knots Fig. 1. Supporting knots around a reference knot. 730 H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 by a trial-and-error process. However, it will greatly increase the computational cost. In the present MLSFD method, the least square technique is introduced to overcome this drawback. The least-squares technique allows an optimized approximation derived from an over-determined set of equations, and the resultant coeﬃcient matrix has good properties such as positive, symmetric, and deﬁnite. Therefore, the problem of singular coeﬃcient matrix is circumvented by means of using more supporting points (more than the unknowns). More detailed descriptions about MLSFD method can be found in . We only give the summary in this section. The explicit expression for the optimal derivative approximation by least-square technique is du ¼ ðST SÞ1 ST Du ð3Þ Usually, the row number of matrix S is larger than its column number. If they are of the same value and matrix S is invertible, the derivative vector derived from Eq. (3) has the equivalent entries achieved by Eq. (2). It should be noted that the use of least-square technique does not degrade the order of accuracy of ﬁnite-diﬀerence approximation . From the formulation shown above, we can see that if we approximate a second-order derivative with the second-order of accuracy by MLSFD method, at least 10 nodes are required (including the reference node itself). But for the traditional FD scheme on a uniform mesh, the use of three nodes can achieve the same order of accuracy. More supporting points imply larger computational costs. In this sense, MLSFD method is much less eﬃcient than conventional FD schemes. However, MLSFD method can be easily applied at scattered points, while the conventional FD schemes are only applied along the mesh lines on a structured grid. When the problems with complex geometry are confronted, conventional FD scheme usually requires a body-ﬁtted mesh, which is not necessary for the MLSFD method. From the end-userÕs standpoint, it would be very attractive if a method is simple, eﬃcient and able to deal with complex and large scale problem. If we can deliberately combine the MLSFD and FD methods and let them deal with their appropriate jobs respectively, we should be able to reach a balance point between the computational eﬃciency and the ability of dealing with complex geometry. This is the motivation of present study. To fully exploit the performance of present method, a special composite mesh is designed to suit for our purpose. In the composite mesh, a Cartesian mesh is generated as the background mesh, and the nodes around the solid body are generated with the consideration of the geometrical description to make the meshing easier. A typical example of such composite meshes is shown in Fig. 2. During the mesh generation, the nodes among the interface between the meshes need to be identiﬁed for the convenience of the numerical discretization in the next step. To ﬁgure out this concept, a detailed local pattern of node distribution is demonstrated in Fig. 3. In this ﬁgure, the symbol of circle represents the group of nodes at which MLSFD method is used to do numerical discretization. This group of nodes is called as ‘‘mesh-free region’’. The symbol of diamond represents the other group at which the conventional FD scheme is used, Fig. 2. Two diﬀerent ways of knot generation for present study. H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 731 Fig. 3. Structure of ‘‘exchange layers’’. and they are named as ‘‘FD region’’. The nodes among the interface have been marked by ﬁlled symbols, which separate the two regions. In fact, the nodes in this region are on the Cartesian mesh. However, two diﬀerent discretization methods are employed on them as indicated by the circle and diamond symbols. These interface nodes play a very important role in the current scheme. They act as the ‘‘exchange layers’’ and bridge the mesh-free and FD regions. Unlike the methods applied on the overlapped mesh, the information exchange between MLSFD and FD methods does not depend on the interpolation. Instead, the information exchange will take place naturally in the ‘‘exchange layers’’. The layer of ﬁlled circles serves as the boundary for the mesh-free regions, and the layer of ﬁlled diamonds serves as the boundary for the FD regions. The algorithm is similar to the domain-decomposition though there is no explicit boundary line to divide the solution domain. The concept of ‘‘exchange layers’’ is workable due to the contribution of truly mesh-free scheme. Because there is no requirement of connection between the nodes in the domain, the local support of ‘‘mesh-free’’ nodes near or inside the exchange layers can include the points in the FD region generated by Cartesian mesh. As discussed previously, MLSFD method used in present scheme is a truly mesh-free method. Therefore, this design of inter-woven exchange layer ensures the smooth transmitting of information between the two diﬀerent regions. 3. Flow past a circular cylinder The capabilities of the present method are demonstrated by simulating the laminar ﬂow over a circular cylinder. This ﬂow has been computed extensively and used as a benchmark to examine the accuracy of new numerical methods for a long time. It is well-known that the ﬂow exhibits vastly diﬀerent patterns as the Reynolds number Re ¼ U1t D changes, where U1 is the free-stream velocity, D is the cylinder diameter, and t is the kinematic viscosity. At small Reynolds numbers, i.e., from zero up to approximately Recritical ﬃ 49, the ﬂow maintains a stable pattern with a pair of symmetric counter-rotating vortices behind the cylinder. At moderate Reynolds number, i.e., 50 < Re < 190, though the ﬂow remains laminar and two-dimensional, vortex shedding, also known as the Karman vortex street, can be observed. Because the Karman vortexstreet developing behind the cylinder displays clearly time periodicity, this ﬂow problem was often used as a prototype of unsteady separated ﬂows. In the present study, we performed numerical simulation at a series of Reynolds number from 10 to 200 with various ﬂow patterns in the steady and unsteady state. The eﬃciency improvement was also examined using present method. 732 H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 Fig. 4. Conﬁguration of ﬂow past a circular cylinder. 3.1. Stream function-vorticity formulation The ﬂow domain and boundary conditions are depicted in Fig. 4. The non-dimensional governing equations, expressed in terms of vorticity and stream function, are written as ox ox ox 1 þu þv ¼ ot ox oy Re o2 x o2 x þ 2 ; ox2 ox o 2 w o2 w þ ¼ x; ox2 oy 2 ð4Þ ð5Þ where w denotes stream function, x represents vorticity, and u, v denote the components of velocity in the x and y direction, which can be calculated from the stream function u¼ ow ; oy v¼ ow ox ð6Þ 3.2. Boundary conditions and initial conditions In this study, the in-ﬂow velocity is speciﬁed with a free stream velocity U , which is equivalent to impose the boundary condition for streamfunction with w ¼ U y and vorticity with x ¼ 0. For the unsteady ﬂow with Re > 49, the top and bottom boundaries are located at a transversal distance of 16 times of the cylinder diameter, which are assumed far enough to be the far-ﬁeld according to literature . For the steady ﬂow with Re < 49, transversal distance is chosen as 10 times of the cylinder diameter to save the computational cost. The boundary condition imposed at these places is the same as the in-ﬂow boundary. The out-ﬂow boundary is located at a distance of 30 times of the cylinder diameter downstream of the rear of the cylinder. The corresponding boundary condition is imposed in such a manner that it does not inﬂuence the vortices formed in the near wake of the cylinder, which is of the form ox ¼0 ox ð7Þ H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 733 At the surface of the cylinder, the non-slip boundary condition applies, i.e. the velocity of the ﬂuid is equal to the cylinder velocity, zero. Therefore, the stream-function must satisfy two boundary conditions: a notangential ﬂow condition and a no-normal ﬂow condition, i.e. w ¼ constant; ow ¼0 on ð8Þ In the multiply connected domain, the value of stream function at the surface of stationary body is an unknown constant. However, this unknown constant may vary with time, and must be determined and updated at each time step during the period of computation. To calculate the value of stream function on the wall, Tezduyar  and Tezduyar et al.  suggested the use of the single-value condition of the pressure for the multiply connected domain. The pressure single-value condition is derived from the fact that pressure is a scalar. For the present case, it can be formulated as Z 2p op R dh ¼ 0 ð9Þ oh 0 As shown in Fig. 5, R is the radius of the cylinder. The expression of op on the cylinder wall can be derived oh from the momentum equation under the cylindrical coordinates, which is written as op o3 w o2 w ¼ R 2 ð10Þ oh cylinder wall or3 or Substituting Eq. (10) into Eq. (9) gives Z 2p 3 o w o2 w R 3 þ 2 dh ¼ 0 or or 0 ð11Þ The derivatives in Eq. (11) can be discretized by the one-side ﬁnite diﬀerence schemes and expressed in terms of the stream function at the cylinder wall and interior knots. Then, the constant value of stream function at the cylinder surface can be calculated and updated at each time step. There is no explicit boundary condition available for vorticity x at the surface of cylinder, and it is directly computed from the stream-function. Fig. 5. Local body-ﬁtted coordinate system. 734 H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 In this study, all the simulations are carried out in a time-marching way. An unsymmetrical initial ﬂow ﬁeld, which serves as an artiﬁcial initiator for the numerical simulation, is provided. qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 2 2 w ¼ ðx xc Þ þ ðy yc Þ where xc and yc are the position of the center of cylinder in the x and y coordinates, respectively. For the Reynolds numbers less than the critical value (Recritical ¼ 49), the introduced perturbation is gradually dissipated by viscosity. For the Reynolds numbers greater than the critical value, the introduced perturbation will trigger the vortex shedding process, and make the ﬂow ﬁeld eventually settle into a periodic oscillatory pattern. 3.3. Node generation In principle, the concept of mesh-free least square-based ﬁnite diﬀerence discretization works well for any ‘‘mesh’’ system, in which the nodes can be either regularly or irregularly distributed. This makes MLSFD method possess geometry ﬂexibility. However, in our cases,, only the circular cylinder was considered. For simplicity and convenience, the nodes in the neighborhood of the circular cylinder are generated under the cylindrical coordinate system. Beyond this neighborhood, Cartesian mesh, where the second order ﬁnite diﬀerence scheme is employed to do the spatial discretization, is appropriately generated. Since the ﬂow ﬁeld in the ‘‘FD region’’ does not vary intensively, the mesh size in the Cartesian mesh can be relatively coarser accordingly. The interface between two ‘‘meshes’’ can be of any shape, depending on the convenience and requirement of the end-users. For example, as shown in Fig. 2, two diﬀerent types of interfaces are generated with regard to the present problem concerning circular cylinder: one has the shape of square, and the other has the shape of circle. Of the two types, the former one is comparatively convenient in programming, and the latter one has more evenly distributed nodes, which implies similar discretization error in the ‘‘exchange layers’’. Both are suitable for the solution of current cases. 3.4. Calculation of lift and drag coeﬃcients For viscous ﬂow, the forces exerted on the body come from two sources: surface pressure distribution and surface friction.The surface pressure distribution can be derived from the tangential momentum equation at the surface of the body. As shown in Fig. 5, for the body-ﬁtted local orthogonal coordinate * * system ð s ; n Þ, the tangential momentum equation at the surface of a stationary body can be simpliﬁed as 1 op * * ¼ t s ðr xÞ; q os ð12Þ which can be further simpliﬁed as 1 op ox ¼ t q os on ð13Þ The term on the right-hand side of Eq. (13) is the rate of vorticity creation per unit length at the solid surface, which can be calculated after implementing the boundary conditions. For the particular case of ﬂow around a circular cylinder, the drag and lift can be obtained from Z ox FL ¼ r lr lx cos h dh on ð14Þ Z ox lr lx sin h dh FD ¼ r on H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 735 where r is the radius of the circular cylinder. The non-dimensional drag and lift coeﬃcients of the body are then given by CL ¼ FD ; 2 r qU1 CD ¼ FD 2 r qU1 ð15Þ 4. Results and discussion In this study, the mixed MLSFD and the conventional FD schemes are combined to simulate ﬂow around an impulsively started circular cylinder for various Reynolds numbers. Unless otherwise mentioned, the basic method consists of a second-order MLSFD and central diﬀerence scheme. In the MLSFD scheme, 16 supporting points are employed for every centre node, and central diﬀerence scheme is applied on the uniform Cartesian mesh. The number of the time steps used is related to both the ﬂow complexity and the ‘‘mesh’’ size. A four-stage Runge–Kutta method is used for temporal discretization. To demonstrate the improved eﬃciency of the proposed method, simulations of steady case are also carried out and compared quantitatively with the hybrid simulations. Current computation was performed on a Compaq Alfa machine. The numerical results achieved by the present method are also compared with previous numerical studies. One of the most attractive features in the numerical simulation of ﬂow past a cylinder is the vortex shedding behind the cylinder and the periodic variation of the ﬂow ﬁeld with moderate Reynolds number. The ability of the method to simulate the unsteady ﬂow was demonstrated by simulating the ﬂow ﬁeld with Reynolds number of 100 and 200. 4.1. Simulation of steady ﬂow at low Reynolds numbers An important consideration regarding to the use of present method is the running-time required by the numerical simulation. To estimate the improvement of the computational eﬃciency by present method as compared with the MLSFD method, numerical experiment was carried out for Reynolds number of 20. The reported data are obtained on a ‘‘mesh’’, which was divided into two regions: mesh-free region and FD region, as shown in Fig. 3. The numerical simulation was maintained for a period of about 30 nondimensional time length to reach a ﬁnal steady state, requiring 20,000 time steps with Dt ¼ 0:0015. Table 1 presents the computational cost required to reach the steady-state solution by MLSFD and present hybrid methods. It can be seen in Table 1 that when the mesh-free region holds about 33.54% area of the whole domain, present method only uses half of the running time required for the pure MLSFD method. This fact indicates that the present method leads to signiﬁcant reduction in computational cost, and achieves signiﬁcant improvement on the eﬃciency as compared with pure MLSFD method. As pointed out before, the eﬃciency improvement is mainly due to the less summation and multiplication required by the central diﬀerence scheme within one time step. Therefore, the extent of eﬃciency improvement mainly depends on the size of the area where the conventional FD scheme is applied. It can be observed in Table 1, that the smaller of the mesh-free region, the greater the eﬃciency improved. This eﬀect is also illustrated in Fig. 6 by eﬃciency improvement versus percentage of mesh-free area. Numerical simulations were also carried out for two other small Reynolds numbers: 10 and 40, respectively. Fig. 7 illustrates the streamlines when ﬂow reaches its ﬁnal steady state. In all cases, a pair of vortices develops behind the cylinder and is perfectly aligned. This is consistent with previous observation. Some quantitative parameters for the recirculating region, such as the length of the recirculating region, Lsep from the rearmost point of the cylinder to the end of the wake, separation angle hsep and drag coeﬃcient Cd , as well as the results from other researchers [19–22], are listed in Table 2. All these ﬂow parameters agree well with the results of previous studies for all three Reynolds numbers studied. 736 H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 Table 1 Eﬃciency comparison between MLSFD and present method for Re ¼ 20 Total node number Method Percent of mesh-free region (%) CPU-time (s) Eﬃciency improvement (%) 10 148 Present MLSFD 33.54 100 56.5 113.04 50.02 8364 Present MLSFD 43.88 100 54.72 89.49 38.85 7488 Present MLSFD 51.38 100 51.61 77.00 32.97 7037 Present MLSFD 65.21 100 57.37 74.56 23.06 5716 Present MLSFD 71.02 100 47.73 59.85 20.25 Fig. 6. Eﬃciency improvement versus meshfree region for Re ¼ 20. 4.2. Simulation of unsteady ﬂow at moderate Reynolds numbers An understanding of the vortex shedding process behind a circular cylinder poses a challenge to both basic research and general applications. There are many studies on this topic in the literature [23–25]. In this section, a number of comparisons are made between the present study and other numerical work for Reynolds numbers of Re ¼ 100 and 200, which represent the typical cases in terms of validation of new numerical approaches. Figs. 8 and 9 show the time-dependent behavior of streamlines for both Re ¼ 100 and 200. It is clear that the periodicity of ﬂow ﬁeld has been successfully revealed in both plots. With the periodic vortex shedding from the rear part of the cylinder, the value of stream-function at the cylinder surface also experiences H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 737 Fig. 7. Streamlines for Re ¼ 10, 20 and 40. Table 2 Comparison of length of the recirculating region (Lsep ), separation angle (hsep ) and drag coeﬃcient (Cd ) for Re ¼ 10, 20 and 40 Re Source Lsep hsep Cd 10 Dennis et al.  Takami et al.  Tuann et al.  Fornberg  Present 0.252 0.249 0.25 – 0.252 29.6 29.3 29.7 – 30.0 2.85 2.80 3.18 – 3.07 20 Dennis et al.  Takami et al.  Tuann et al.  Fornberg  Present 0.94 0.935 0.9 0.91 0.93 43.7 43.7 44.1 – 44.1 2.05 2.01 2.25 2.00 2.18 40 Dennis et al.  Takami et al.  Tuann et al.  Fornberg  Present 2.35 2.32 2.1 2.24 2.20 53.8 53.6 54.8 – 53.5 1.522 1.536 1.675 1.498 1.713 periodic variation as shown in Figs. 10 and 11. It can be observed that the magnitude of the variation grows up as the Reynolds number increases, though the value is still comparatively small. For the ﬂow past blunt bodies, the drag and lift coeﬃcients at the surface of body are two very important parameters. The timeevolution of these two characteristic parameters illustrate the variation of the ﬂow ﬁeld. It can be observed 738 H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 Fig. 8. Time-evolution of streamlines for Re ¼ 100. H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 Fig. 9. Time-evolution of streamlines for Re ¼ 200. 739 740 H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 0.02 Stream function 0.01 0 -0.01 -0.02 100 110 120 Time Fig. 10. The time-evolution of stream function value on the cylinder surface for Re ¼ 100. 0.02 Stream function 0.01 0 -0.01 -0.02 -0.03 100 110 120 Time Fig. 11. The time-evolution of stream function value on the cylinder surface for Re ¼ 200. in Figs. 12 and 13 that lift and drag coeﬃcients show obvious periodic oscillations for both Re ¼ 100 and 200 cases. This implies the periodic variation of ﬂow ﬁeld. From Figs. 12 and 13, it can also be found that the lift coeﬃcient oscillates with larger amplitude than the drag coeﬃcient, and the drag coeﬃcient varies twice as fast as the lift coeﬃcient. These phenomena are consistent with those observed by other researchers. The reason lies on that the drag coeﬃcient is aﬀected by vortex shedding process from both H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 0.4 1.38 0.2 1.37 0 741 1.36 -0.2 Cl -0.6 1.34 -0.8 1.33 -1 1.32 -1.2 1.31 -1.4 Cl Cd -1.6 1.3 1.29 -1.8 -2 Cd 1.35 -0.4 100 110 120 1.28 130 Time 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 -1.4 -1.6 -1.8 -2 -2.2 -2.4 1.8 1.7 1.6 1.5 1.4 Cd Cl Fig. 12. The time-evolution of lift and drag coeﬃcients for Re ¼ 100. 1.3 1.2 Cl Cd 90 100 1.1 110 120 1 Time Fig. 13. The time-evolution of lift and drag coeﬃcients for Re ¼ 200. sides of the cylinder. In Tables 3 and 4, the present results are quantitatively compared with the numerical results of Belov , Braza  and Liu et al.  for Re ¼ 100 and 200, respectively. It can be observed that excellent agreement has been achieved. 742 H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 Table 3 Drag coeﬃcients for Re ¼ 100 and 200 Drag (Cd ) Re ¼ 100 Re ¼ 200 Belov et al.  Braza et al.  Liu et al.  Present – 1.364 ± 0.015 1.350 ± 0.012 1.325 ± 0.008 1.19 ± 0.042 1.40 ± 0.05 1.31 ± 0.049 1.327 ± 0.045 Lift (Cl ) Re ¼ 100 Re ¼ 200 Belov et al.  Braza et al.  Liu et al.  Present – ±0.25 ±0.339 ±0.28 ±0.64 ±0.75 ±0.69 ±0.60 Strouhal number (St) Re ¼ 100 Re ¼ 200 Belov et al.  Braza et al.  Liu et al.  Present – 0.160 0.164 0.164 0.193 0.200 0.192 0.196 Table 4 Lift coeﬃcients for Re ¼ 100 and 200 Table 5 Strouhal number for Re ¼ 100 and 200 To determine whether we have achieved the proper shedding frequency, the Strouhal number (St) is used as a measure of the oscillating ﬂuid ﬂow phenomenon in the wake region. It is connected to the ﬂow parameter by the following relationship St ¼ fs D=U1 ð16Þ where fs is the shedding frequency. Since D and U1 are non-dimensionalized as unity, St is equivalent to the shedding frequency. Thus, the Strouhal number is obtained as the inverse for the time period of vortex shedding. In the present study, Strouhal numbers are deduced from the time-evolution of the lift coeﬃcient Cl . It can be seen from Table 5 that, the St number computed using our algorithm was very close to that reported by Liu et al.  for Re ¼ 100 and 200. Our value of Strouhal number for Re ¼ 200 diﬀers from that reported by them only about 2%, and is of the same value for Re ¼ 100. 5. Conclusions In this paper, a hybrid FD/MLSFD method was developed for solving the two-dimensional stream function-vorticity formulation of N–S equations. The method is validated by simulating the ﬂow past a circular cylinder at various Reynolds numbers. The present numerical results were in good agreement with H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744 743 other experimental and numerical studies. It indicates that the combination of conventional FD scheme and MLSFD method is eﬀective for solving large scale ﬂow problems. From the viewpoint of methodology, the present approach greatly improves the computational eﬃciency as compared with MLSFD method, and it is easy to deal with complex geometry. In the future, simple eﬀective tools will be developed on top of this basic framework to handle more complicated problems. It is clear that the above practice and experiments provide only a starting point for the potential end-user, especially when the simulations of multi-dimensional external ﬂows are to be considered. References  T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Int. J. Numer. Meth. Engng. 37 (1994) 229–256.  W. Liu, S. Jun, Y. Zhang, Reproducing kernel particle methods, Int. J. Numer. Methods Fluids 20 (1995) 1081– 1106.  I. Babuska, J. Melenk, The partition of unity method, Int. J. Numer. Meth. Engng. 40 (1997) 727–758.  C.A. Duarte, J.T. Oden, Hp clouds-a meshless method to solve boundary-value problems, TICAM Report 95-05, 1995.  E. O nate, S. Idelsohn, O.C. Zienkiewicz, R.L. Taylor, A ﬁnite point method in computational mechanics. Application to convective transport and ﬂuid ﬂow, Int. J. Numer. Meth. Engng. 39 (1996) 3839–3866.  S.N. Atluri, T. Zhu, New meshless local Petrov–Galerkin (MLPG) approach in computational mechanics, Comput. Mech. 22 (2) (1998) 117–127.  T. Liszka, An interpolation method for an irregular net of nodes, Int. J. Numer. Meth. Engng. 20 (1984) 1599–1612.  E.J. Kansa, Multiquadrics––a scattered data approximation scheme with applications to computational ﬂuid-dynamics. I. Surface approximations and partial derivative estimates, Comput. Math. Appl. 19 (6–8) (1990) 127–145.  E.J. Kansa, Multiquadrics––a scattered data approximation scheme with applications to computational ﬂuid-dynamics––II. Solutions to parabolic, hyperbolic, and elliptic partial diﬀerential equations, Comput. Math. Appl. 19(6–8) 147–161 (1990) 29; H. Takami, H.B. Keller, Steady two-dimensional viscous ﬂow of an Incompressible ﬂuid past a circular cylinder, Phys. Fluids 12 (Suppl. II) (1969) II-51.  H. Ding, C. Shu, K.S. Yeo, Development of least square-based two-dimensional ﬁnite diﬀerence schemes and their application to simulate natural convection in a cavity, Comput. Fluids 33 (2004) 137–154.  D. Calhoun, R.J. Leveque, A Cartesian grid ﬁnite-volume method for the advection–diﬀusion equation in irregular geometries, J. Comput. Phys. 157 (2000) 143–180.  R.B. Pember, J.B. Bell, An adaptive cartesian grid method for unsteady compressible ﬂow in irregular regions, J. Comput. Phys. 120 (1995) 278–304.  J. Falcovitz, G. Alfandary, G. Hanoch, A two-dimensional conservation laws scheme for compressible ﬂows with moving boundaries, J. Comput. Phys. 138 (1997) 83–102.  M. Hinatsu, J.H. Ferziger, Numerical computation of unsteady incompressible ﬂow in complex geometry using a composite multigrid technique, Int. J. Numer. Methods Fluids 13 (1991) 971–997.  C.Y. Perng, R.L. Street, A coupled multigrid-domain-splitting technique for simulating incompressible ﬂows in geometrically complex domains, Int. J. Numer. Methods Fluids 13 (1991) 269–286.  M. Behr, D. Hastreiter, S. Mittal, T.E. Tezduyar, Incompressible ﬂow past a circular cylinder: dependence of the computed ﬂow ﬁeld on the location of the lateral boundaries, Comput. Methods Appl. Mech. Engng. 123 (1995) 309–316.  T.E. Tezduyar, Finite element formulation for the vorticity-stream function form of the incompressible Euler equations on the multiply connected domains, Comput. Methods Appl. Mech. Engng. 73 (1989) 331–339.  T.E. Tezduyar, R. Glowinski, J. Liou, Petrov–Galerkin methods on multiply connected domains for the vorticity-stream function formulation of the incompressible Navier–Stokes equations, Int. J. Numer. Methods Fluids 8 (1988) 1269–1290.  S.C.R. Dennis, G.Z. Chang, Numerical solutions for steady ﬂow past a circular cylinder at Reynolds numbers up to 100, J. Fluid Mech. 42 (1970) 471.  H. Takami, H.B. Keller, Steady two-dimensional viscous ﬂow of an incompressible ﬂuid past a circular cylinder, Phys. Fluids 12 (Suppl. II) (1969) II-51.  S.Y. Tuann, M.D. Olson, Numerical studies of the ﬂow around a circular cylinder by a ﬁnite element method, Comput. Fluids 6 (1978) 219.  B. Fornberg, A numerical study of steady viscous ﬂow past a circular cylinder, J. Fluid Mech. 98 (1980) 819.  A. Belov, L. Martinelli, A. Jameson, A new implicit algorithm with multigrid for unsteady incompressible ﬂow calculation, AIAA 95-0049 January, 9 (1995). 744 H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744  M. Braza, P. Chassaing, H. Ha Minh, Numerical study and physical analysis of the pressure and velocity ﬁelds in the near wake of a circular cylinder, J. Fluid Mech. 165 (1986) 79.  C. Liu, X. Zheng, C.H. Sung, Preconditioned multigrid methods for unsteady incompressible ﬂows, J. Comput. Phys. 139 (1998) 35.