...

Ding2004.pdf

by user

on
Category: Documents
2

views

Report

Comments

Transcript

Ding2004.pdf
Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744
www.elsevier.com/locate/cma
Simulation of incompressible viscous flows past
a circular cylinder by hybrid FD scheme and meshless
least square-based finite difference 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 efficient method for simulating the two-dimensional steady and unsteady incompressible
flows. The method is a hybrid approach, which combines the conventional finite difference (FD) scheme and the meshfree least square-based finite difference (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 flow domain to take advantage of its high computational efficiency. Numerical simulation is
carried out for the steady flow with Reynolds numbers of 10, 20 and 40, and unsteady flow 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 efficiency.
Ó 2003 Elsevier B.V. All rights reserved.
Keywords: Incompressible viscous flows; 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 fluid 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-specified 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 scientific and engineering problems. Despite of their preliminary
successes, they do have some drawbacks, i.e., in general, they are less efficient than the conventional
numerical methods such as finite difference (FD), finite element (FE) and finite 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
finite difference 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 efficiency is unavoidably impaired. In this regard, the mesh-free methods achieve the
geometry flexibility at the cost of computational efficiency. Despite of recent advances in the fields of meshfree methodology, it is certainly the bottleneck in the industrial applications.
To handle flows with complex geometry effectively and efficiently is always a tough task. Traditional FD/
FV methods enjoy the computational efficiency, but they confront difficulties 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 flow
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-fitted or unstructured mesh,
this method is thought to be more efficient. 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 difficulty 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 flow 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 fit 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 finite difference methods, to deal with flow problems of complex geometry. The main
drawback of using this kind of mesh is that the information exchange between different grids is based on the
interpolation. If the solution exhibits intensive variation near the interface between grids, the interpolation
process may introduce significant numerical errors or cause the break-down of the computation.
To improve the above methods, a hybrid method, which combines the conventional finite difference
scheme and the mesh-free least square-based finite difference (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 flow domain to take advantage
of its high computational efficiency. The present method is validated by simulating the flow field around a
circular cylinder at various Reynolds numbers. The simulated flow patterns include both the steady and
unsteady state. Numerical experiments are also carried out to examine the efficiency 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 firstly give a brief review on the ‘‘truly’’ mesh-free least square-based finite difference
(MLSFD) method. Then, we show how the MLSFD and conventional FD schemes are combined to solve
the partial differential 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 finite difference 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 first-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 difficulty is not easy to be resolved because little is known about the effects 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 coefficient matrix has good properties such as positive, symmetric, and definite. Therefore, the
problem of singular coefficient matrix is circumvented by means of using more supporting points (more
than the unknowns). More detailed descriptions about MLSFD method can be found in [10]. 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
finite-difference approximation [10].
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 efficient 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-fitted 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, efficient 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
efficiency 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 identified for the convenience of the
numerical discretization in the next step. To figure out this concept, a detailed local pattern of node distribution is demonstrated in Fig. 3. In this figure, 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 different 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 filled symbols,
which separate the two regions. In fact, the nodes in this region are on the Cartesian mesh. However, two
different 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 filled circles serves
as the boundary for the mesh-free regions, and the layer of filled 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 different regions.
3. Flow past a circular cylinder
The capabilities of the present method are demonstrated by simulating the laminar flow over a circular
cylinder. This flow 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 flow exhibits vastly different 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 ffi 49, the
flow 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 flow 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 flow problem was often used as a
prototype of unsteady separated flows. In the present study, we performed numerical simulation at a series
of Reynolds number from 10 to 200 with various flow patterns in the steady and unsteady state. The
efficiency improvement was also examined using present method.
732
H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744
Fig. 4. Configuration of flow past a circular cylinder.
3.1. Stream function-vorticity formulation
The flow 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-flow velocity is specified 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 flow
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-field according to literature [16]. For the steady flow
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-flow boundary. The out-flow
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 influence 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 fluid is equal
to the cylinder velocity, zero. Therefore, the stream-function must satisfy two boundary conditions: a notangential flow condition and a no-normal flow 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 [17] and Tezduyar et al. [18] 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 finite difference 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-fitted 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 flow
field, which serves as an artificial initiator for the numerical simulation, is provided.
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
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 flow field eventually settle into a periodic
oscillatory pattern.
3.3. Node generation
In principle, the concept of mesh-free least square-based finite difference discretization works well for
any ‘‘mesh’’ system, in which the nodes can be either regularly or irregularly distributed. This makes
MLSFD method possess geometry flexibility. 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 finite difference scheme is employed to do the spatial discretization, is appropriately generated. Since the flow field 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 different 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 coefficients
For viscous flow, 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-fitted local orthogonal coordinate
* *
system ð s ; n Þ, the tangential momentum equation at the surface of a stationary body can be simplified as
1 op
*
*
¼ t s ðr xÞ;
q os
ð12Þ
which can be further simplified 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
flow 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 coefficients 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 flow
around an impulsively started circular cylinder for various Reynolds numbers. Unless otherwise mentioned,
the basic method consists of a second-order MLSFD and central difference scheme. In the MLSFD scheme,
16 supporting points are employed for every centre node, and central difference scheme is applied on the
uniform Cartesian mesh. The number of the time steps used is related to both the flow complexity and the
‘‘mesh’’ size. A four-stage Runge–Kutta method is used for temporal discretization. To demonstrate
the improved efficiency 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 flow past a cylinder is the vortex
shedding behind the cylinder and the periodic variation of the flow field with moderate Reynolds number.
The ability of the method to simulate the unsteady flow was demonstrated by simulating the flow field with
Reynolds number of 100 and 200.
4.1. Simulation of steady flow 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 efficiency 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 final 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 significant reduction in computational cost, and achieves significant improvement on the efficiency as compared with pure MLSFD method. As pointed out before, the
efficiency improvement is mainly due to the less summation and multiplication required by the central
difference scheme within one time step. Therefore, the extent of efficiency 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 efficiency improved. This effect is also illustrated in Fig. 6 by
efficiency 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 flow reaches its final 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 coefficient Cd ,
as well as the results from other researchers [19–22], are listed in Table 2. All these flow 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
Efficiency comparison between MLSFD and present method for Re ¼ 20
Total node number
Method
Percent of mesh-free
region (%)
CPU-time (s)
Efficiency 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. Efficiency improvement versus meshfree region for Re ¼ 20.
4.2. Simulation of unsteady flow 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 flow field 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 coefficient (Cd ) for Re ¼ 10, 20 and 40
Re
Source
Lsep
hsep
Cd
10
Dennis et al. [19]
Takami et al. [20]
Tuann et al. [21]
Fornberg [22]
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. [19]
Takami et al. [20]
Tuann et al. [21]
Fornberg [22]
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. [19]
Takami et al. [20]
Tuann et al. [21]
Fornberg [22]
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 flow past blunt
bodies, the drag and lift coefficients at the surface of body are two very important parameters. The timeevolution of these two characteristic parameters illustrate the variation of the flow field. 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 coefficients show obvious periodic oscillations for both Re ¼ 100 and
200 cases. This implies the periodic variation of flow field. From Figs. 12 and 13, it can also be found that
the lift coefficient oscillates with larger amplitude than the drag coefficient, and the drag coefficient varies
twice as fast as the lift coefficient. These phenomena are consistent with those observed by other
researchers. The reason lies on that the drag coefficient is affected 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 coefficients 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 coefficients for Re ¼ 200.
sides of the cylinder. In Tables 3 and 4, the present results are quantitatively compared with the numerical
results of Belov [23], Braza [24] and Liu et al. [25] 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 coefficients for Re ¼ 100 and 200
Drag (Cd )
Re ¼ 100
Re ¼ 200
Belov et al. [23]
Braza et al. [24]
Liu et al. [25]
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. [23]
Braza et al. [24]
Liu et al. [25]
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. [23]
Braza et al. [24]
Liu et al. [25]
Present
–
0.160
0.164
0.164
0.193
0.200
0.192
0.196
Table 4
Lift coefficients 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 fluid flow phenomenon in the wake region. It is connected to the flow
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 coefficient
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. [25] for Re ¼ 100 and 200. Our value of Strouhal number for Re ¼ 200 differs 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 flow 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 effective for solving large scale flow problems.
From the viewpoint of methodology, the present approach greatly improves the computational efficiency
as compared with MLSFD method, and it is easy to deal with complex geometry. In the future, simple
effective 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 flows are to be considered.
References
[1] T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Int. J. Numer. Meth. Engng. 37 (1994) 229–256.
[2] W. Liu, S. Jun, Y. Zhang, Reproducing kernel particle methods, Int. J. Numer. Methods Fluids 20 (1995) 1081–
1106.
[3] I. Babuska, J. Melenk, The partition of unity method, Int. J. Numer. Meth. Engng. 40 (1997) 727–758.
[4] C.A. Duarte, J.T. Oden, Hp clouds-a meshless method to solve boundary-value problems, TICAM Report 95-05, 1995.
[5] E. O
nate, S. Idelsohn, O.C. Zienkiewicz, R.L. Taylor, A finite point method in computational mechanics. Application to
convective transport and fluid flow, Int. J. Numer. Meth. Engng. 39 (1996) 3839–3866.
[6] S.N. Atluri, T. Zhu, New meshless local Petrov–Galerkin (MLPG) approach in computational mechanics, Comput. Mech. 22 (2)
(1998) 117–127.
[7] T. Liszka, An interpolation method for an irregular net of nodes, Int. J. Numer. Meth. Engng. 20 (1984) 1599–1612.
[8] E.J. Kansa, Multiquadrics––a scattered data approximation scheme with applications to computational fluid-dynamics. I. Surface
approximations and partial derivative estimates, Comput. Math. Appl. 19 (6–8) (1990) 127–145.
[9] E.J. Kansa, 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. 19(6–8) 147–161 (1990) 29;
H. Takami, H.B. Keller, Steady two-dimensional viscous flow of an Incompressible fluid past a circular cylinder, Phys. Fluids 12
(Suppl. II) (1969) II-51.
[10] H. Ding, C. Shu, K.S. Yeo, Development of least square-based two-dimensional finite difference schemes and their application to
simulate natural convection in a cavity, Comput. Fluids 33 (2004) 137–154.
[11] D. Calhoun, R.J. Leveque, A Cartesian grid finite-volume method for the advection–diffusion equation in irregular geometries, J.
Comput. Phys. 157 (2000) 143–180.
[12] R.B. Pember, J.B. Bell, An adaptive cartesian grid method for unsteady compressible flow in irregular regions, J. Comput. Phys.
120 (1995) 278–304.
[13] J. Falcovitz, G. Alfandary, G. Hanoch, A two-dimensional conservation laws scheme for compressible flows with moving
boundaries, J. Comput. Phys. 138 (1997) 83–102.
[14] M. Hinatsu, J.H. Ferziger, Numerical computation of unsteady incompressible flow in complex geometry using a composite
multigrid technique, Int. J. Numer. Methods Fluids 13 (1991) 971–997.
[15] C.Y. Perng, R.L. Street, A coupled multigrid-domain-splitting technique for simulating incompressible flows in geometrically
complex domains, Int. J. Numer. Methods Fluids 13 (1991) 269–286.
[16] M. Behr, D. Hastreiter, S. Mittal, T.E. Tezduyar, Incompressible flow past a circular cylinder: dependence of the computed flow
field on the location of the lateral boundaries, Comput. Methods Appl. Mech. Engng. 123 (1995) 309–316.
[17] 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.
[18] 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.
[19] S.C.R. Dennis, G.Z. Chang, Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100, J. Fluid
Mech. 42 (1970) 471.
[20] H. Takami, H.B. Keller, Steady two-dimensional viscous flow of an incompressible fluid past a circular cylinder, Phys. Fluids 12
(Suppl. II) (1969) II-51.
[21] S.Y. Tuann, M.D. Olson, Numerical studies of the flow around a circular cylinder by a finite element method, Comput. Fluids 6
(1978) 219.
[22] B. Fornberg, A numerical study of steady viscous flow past a circular cylinder, J. Fluid Mech. 98 (1980) 819.
[23] A. Belov, L. Martinelli, A. Jameson, A new implicit algorithm with multigrid for unsteady incompressible flow calculation, AIAA
95-0049 January, 9 (1995).
744
H. Ding et al. / Comput. Methods Appl. Mech. Engrg. 193 (2004) 727–744
[24] M. Braza, P. Chassaing, H. Ha Minh, Numerical study and physical analysis of the pressure and velocity fields in the near wake of
a circular cylinder, J. Fluid Mech. 165 (1986) 79.
[25] C. Liu, X. Zheng, C.H. Sung, Preconditioned multigrid methods for unsteady incompressible flows, J. Comput. Phys. 139 (1998)
35.
Fly UP