INTRODUCTION
The numerical simulation of incompressible flows is one of the problems in computational fluid dynamics which has received continuous attention in recent years. Although researchers use the compressible flow schemes for incompressible ones there was some difficulties in simulation of incompressible flows. Stiffness of the convective terms in the governing equations of lowspeed flows causes convergence and accuracy problems and increases sensitivity to boundary conditions (Volpe, 1993). To overcome these problems CB method for incompressible flows was developed by Drikakis et al. (1994) which is based on the definition of primitive variables as functions of their values on the characteristics and formation of compatibility equations in characteristic directions. Afterwards, this method has been developed and applied successfully in computation of different incompressible flow cases (Shapiro and Drikakis, 2005, 2006; Drikakis and Smolarkiewicz, 2001; Mallinger and Drikakis, 2002). Neofytou (2007) proposed revisions in the derivation of the compatibility equations which are used to derive the CB formulas for the calculations of the primitive flow valuables within the iterative process of the scheme for both the 2D and 3D versions which is directly related to the convergence rate. They used a thirdorder of accuracy formula to interpolate characteristic variables on the cell faces. It is found that over simplifying assumptions have been made in the original CB scheme for the derivation of the compatibility equations. At the same time, Zhao and Zhang (2000) have developed a general higherorder finitevolume unstructured grid method for incompressible flows based on the CB method and ACM. The above references show the CB scheme has been successfully used in a number of applications from accuracy and convergence rate viewpoints.
In this study, CB method of Roe (1997) which originally developed for the simulation
of compressible Euler equations was extended to estimate convective fluxes of
incompressible viscous equations. Presented scheme is based on the propagation
of virtual acoustic waves which makes it able to damp numerical oscillations
generated by discontinuities of flux calculations at cell interfaces. This means
that the scheme requires no special boundary treatment. By the aid of ACM, governing
equations took a hyperbolicdominated nature which is necessary for implementation
of CB method. It has been shown that the ACM requires a fast convergence scheme
at each time step in order to satisfy the incompressibility condition (Kwak
et al., 2004). To achieve faster convergence rate, convergence acceleration
techniques such as local time stepping and implicit residual smoothing were
utilized. Viscous and thermal conduction terms were discretized by a common
secondorder technique. Having wide range of stability and high convergence
rate, fifthorder RungeKutta scheme was used for time discretization. Because
of convergence and accuracy problems of natural convection in comparison with
other case studies of lowspeed flows, it was chosen for convergence and accuracy
analysis of presented scheme.
GOVERNING EQUATIONS
The nondimensional form of NavierStokes equations in finitevolume form with
the ACM and the Boussinesq approximation is given by:
Where,
At which for each cell interface, ζ, u_{N}, u_{P}, ΔL,
F_{N}, a, are
respectively angle between normal and horizontal lines, normal and parallel
components of velocity, length, the normal convective flux, the pseudosound
speed and the normal unit vector or the direction of virtual wave propagation
respectively as shown in Fig. 1.

Fig. 1: 
Schematic representation of the primary cell to evaluate normal
convective flux 
An algebraic grid with clustering near the surface was used. The noslip and
constant surface temperature, i.e., u = 0, v = 0, θ = 1 were imposed at
solid boundary. The pressure on the solid wall was found by the extrapolation
using a first order polynomial. The inflow and outflow portions of farfield
boundary were distinguished by At
the inflow region, the pressure was extrapolated from interior domain. Velocity
and temperature were set to environment values and at the outflow region, the
pressure was fixed and other flow variables were interpolated from interior
domain.
Convective flux treatment: A Roelike scheme was used to evaluate convective
flux at any cell side, while the firstorder derivatives of left and right cells
were averaged for viscous and thermal conduction terms. In our work, a similar
approach for incompressible NavierStokes equations based on ACM was used. Considering
Fig. 1, one has:
where, K is attributed to each cell side, namely (1), (2), (3) and (4). By
converting the differential form of the governing equations into quasilinear
form with normal convective flux at any cell side, one has:
The eigenvalues and corresponding eigenvectors of Jacobian matrix A are:
where, first eigenvalue and eigenvector represent a change in velocity direction
and can be thought as a pseudoshear wave which travels with fluid velocity
in the direction of propagation ().
Second one introduces a pseudoentropy wave which travels with the same velocity
as first eigenvalue. Third and forth ones represent virtual acoustic disturbances
which travel with velocities u_{N} ± a. It should be mentioned
that eigenvalues and eigenvectors are real and linearly independent. A Roelike
weighted averaging parameter at cell interface averages flow variables as:
To evaluate virtual acoustic wave strength matrix α_{K}, following
system of equations is solved:
Finally, the normal convective flux at any cell side is obtained as:
The order of accuracy of the scheme depends on choosing W_{R}, W_{L}.
As shown in Fig. 1 for cell side AB, they are obtained as:
The viscous and thermal conduction terms were discretized using the secondary
cells as shown in Fig. 2:
where, φ stands for u, v, θ and K is attributed to each cell side,
namely DA, AB, BC and CD. For example for face AB one has:

Fig. 2: 
Schematic representation of primary and secondary cells to
evaluate viscous terms 
Time discretization: For time discretization, the fifthorder RungeKutta
scheme due to its high accuracy and wide range of stability was used (Jameson,
1995):
where, Q is the right hand side of Eq. 19. The number of
iterations to achieve steadystate was reduced by using a variable time step
which depends on the local flow changes and the grid spacing as:
RESULTS AND DISCUSSION
For model validation, a series of numerical tests were conducted with different Rayleigh numbers around a horizontal circular cylinder. The incompressible laminar twodimensional equations of natural convection were solved using the presented approach. Air was taken as a working fluid with constant Prandtl number of 0.7. The local and average angular Nusselt numbers based on cylinder diameter was determined to compare obtained results with other available data in literature.
To ensure that the solution had no dependency to grid, six levels of grids
namely, 60x60 and 70x70 for Ra = 60, 100x100 and 110x110 for Ra = 400, 120x120
and 130x130 for Ra = 2000 were tested and the results were demonstrated in Fig.
3.

Fig. 3: 
Achieving the grid independency 

Fig. 4: 
Independency of results to the artificial compressibility
parameter 
The grid independency was achieved in a 60x60 grid for 10≤Ra<100, 100x100
for 100≤Ra≤1000 and 120x120 for 1000<Pa≤10000. In fact, the difference
in local Nusselt number distribution appeared to be negligible for two finest
successive levels of grids. This has leaded us to use mentioned grids in all
subsequent calculations considering relative cost of computation with achievable
accuracy.
In order to study the effect of ACM parameter on the solution accuracy, the
angular distribution of local Nusselt number for different values of it are
plotted at Ra = 60 in Fig. 4. This Fig. 4
reveals that the steadystate results are not sensitive to ACM parameter, although
our findings illustrated that its value influences stability and convergence
rate of the solution. For this particular case study ACM parameter of 10 shows
better behavior.

Fig. 5: 
A typical comparison of convergence rates between current
work and flux averaging method for Ra = 1000 

Fig. 6: 
Comparison of average Nusselt number versus different Rayleigh
numbers 
In Fig. 5, the current scheme resulted in a faster convergence
rate with very negligible numerical oscillations that proves the capability
of provided scheme to handle lowspeed flows without any special boundary treatment.
Figure 6 shows substantial agreement between average Nusselt
number calculated by current scheme, flux averaging method of Jameson and benchmark
results. At the same grid size, farfield boundary and ACM parameter, the Jameson
flux averaging failed to converge at Ra>1000, in which the virtual acoustic
wave interactions lead to high and low frequency waves.

Fig. 7: 
Comparison of local Nusselt number for Ra = 1000 and 10000 
The lower frequencies postpone the convergence, but the higher frequencies
violate the numerical stability. Because in the proposed upwind method physical
behavior of flow was related to the mathematical characteristics of the governing
equations, convergence could be archived without any use of artificial numerical
dissipation.
The angular distribution of local Nusselt number obtained from current work
was compared with other researcher’ data at Ra = 1000 and 10000 in
Fig. 7. The peak is observed at the front stagnation point, which was anticipated
since the downstream growth of a thermal boundary layer could certainly increase
the thermal resistance.
Figure 8 displays the streamlines and isotherms for different
Rayleigh numbers. From all these frames, it can be observed that the growth
of thermal boundary layer over the surface of cylinder is markedly high. From
the lower stagnation point, the hot fluid rises up due to the gravity hence
the thickness of the thermal boundary layer is expected to grow. But this phenomenon
is not very straightforward. There appears deepness on the isotherms near to
the surface of the cylinder. At low Rayleigh numbers, the induced flow moves
upward making a turn around the cylinder and its increase resulted in the intensification
of induced flow. Initially the isotherms are approximately concentric circles
indicating pure conduction which then become asymmetrical circles due to convection.
Obtained results showed good agreement with other available results (Saitoh
et al., 1993).

Fig. 8: 
Streamlines and isotherms for (a): Ra = 20, (b): Ra = 600,
(c): Ra = 2000 
CONCLUSION
In present research, a CB secondorder scheme was developed to solve incompressible flows. For convective flux treatment, the propagation of virtual acoustic waves was taken into account which results in high convergence rate with less numerical oscillations without any requirement to damping or special boundary treatment. Governing equations were discretised in time by fifthorder RungeKutta scheme. Local time stepping and residual smoothing was utilized to accelerate convergence. Also the ACM was employed to couple continuity and momentum equations. Viscous and thermal conduction terms were discretized by a secondorder technique. Results obtained by this scheme were compared with benchmark results, from which substantial consensus was observed.
NOMENCLATURE
Dimensionless flux Jacobian matrix 
A 
CourantFriedrichsLevy number 
CFL 
Cylinder diameter (m) 
D 
Dimensionless right eigenvector 
E 
Dimensionless normal convective flux vector 
F_{N} 
Gravity (m sec^{2}) 
g 
Grashof number 
Gr 
Nusselt number 
Nu 
Pressure (Pa) 
P 
Dimensionless pressure 
p 
Prandtl number 
Pr 
Dimensionless viscous and thermal conduction components 
R, S 
Rayleigh number 
Ra 
Reynolds number 
Re 
Temperature (K) 
T 
Dimensionless time 
t 
Velocity components (m sec^{1}) 
U, V 
Dimensionless velocity components 
u,v 
Dimensionless state vector 
W 
Physical coordinates (m) 
X, Y 
Dimensionless coordinates 
y 
Source term 
x,Z 
Greek symbols
Artificial compressibility parameter 
β 
Thermal expansion coefficient (1/K) 
β' 
Dimensionless virtual wave strength 
α 
Thermal diffusivity (m^{2} sec^{1}) 
α' 
Angle between normal and horizontal lines at cell face 
ζ 
Dimensionless primary and secondary cell area 
Ω, Ω' 
Kinematic viscosity (m^{2} sec^{1}) 
v 
Time (sec) 
τ 
Dimensional eigenvalue 
λ 
Density (kg m^{3}) 
ρ 
Roelike weighted averaging parameter 
ω 
Dimensionless temperature 
θ 
Subscripts and superscripts
Normal, parallel 
N, P 
Environment 
∞ 
Right, Left 
R, L 
Wall 
W 
Mean 
 