Subscribe Now Subscribe Today
Research Article

Characteristics-Based Finite-Volume Solution for Natural Convection Around a Horizontal Cylinder

Seyed Esmail Razavi, Farzad Barar and Vahid Farhangmehr
Facebook Twitter Digg Reddit Linkedin StumbleUpon E-mail

In the current investigation, a Characteristics-Based (CB) finite-volume solution was developed to improve the convergence rate and the solution stability of the numerical simulation of incompressible viscous flows. The Artificial Compressibility Method (ACM) as a pre-conditioner was taken into account in our study. A fifth-order Runge-Kutta scheme has been used for marching in time. Local time stepping and implicit residual smoothing were applied as convergence acceleration techniques. The convective fluxes have been calculated by a Roe-like high-order flux splitting scheme. The viscous and thermal conduction terms were found by a second-order technique. This method admitted high Courant-Friedrichs-Levy (CFL) numbers and revealed a good convergence rate and accuracy within a wide range of Rayleigh numbers from 1 to 10000 for natural convection flow around a horizontal circular cylinder which is chosen as a case study. The obtained results were compared with available benchmark data in literature and showed a high consensus with them.

Related Articles in ASCI
Similar Articles in this Journal
Search in Google Scholar
View Citation
Report Citation

  How to cite this article:

Seyed Esmail Razavi, Farzad Barar and Vahid Farhangmehr, 2008. Characteristics-Based Finite-Volume Solution for Natural Convection Around a Horizontal Cylinder. Journal of Applied Sciences, 8: 1905-1911.

DOI: 10.3923/jas.2008.1905.1911



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 low-speed 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 third-order 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 higher-order finite-volume 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 hyperbolic-dominated 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 second-order technique. Having wide range of stability and high convergence rate, fifth-order Runge-Kutta scheme was used for time discretization. Because of convergence and accuracy problems of natural convection in comparison with other case studies of low-speed flows, it was chosen for convergence and accuracy analysis of presented scheme.


The non-dimensional form of Navier-Stokes equations in finite-volume form with the ACM and the Boussinesq approximation is given by:




At which for each cell interface, ζ, uN, uP, ΔL, FN, a, are respectively angle between normal and horizontal lines, normal and parallel components of velocity, length, the normal convective flux, the pseudo-sound 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 no-slip 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 far-field 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 Roe-like scheme was used to evaluate convective flux at any cell side, while the first-order derivatives of left and right cells were averaged for viscous and thermal conduction terms. In our work, a similar approach for incompressible Navier-Stokes 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 quasi-linear 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 pseudo-shear wave which travels with fluid velocity in the direction of propagation (). Second one introduces a pseudo-entropy wave which travels with the same velocity as first eigenvalue. Third and forth ones represent virtual acoustic disturbances which travel with velocities uN ± a. It should be mentioned that eigenvalues and eigenvectors are real and linearly independent. A Roe-like 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 WR, WL. 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 fifth-order Runge-Kutta 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 steady-state was reduced by using a variable time step which depends on the local flow changes and the grid spacing as:



For model validation, a series of numerical tests were conducted with different Rayleigh numbers around a horizontal circular cylinder. The incompressible laminar two-dimensional 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 steady-state 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 low-speed 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, far-field 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


In present research, a CB second-order 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 fifth-order Runge-Kutta 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 second-order technique. Results obtained by this scheme were compared with benchmark results, from which substantial consensus was observed.


Dimensionless flux Jacobian matrix A
Courant-Friedrichs-Levy number CFL
Cylinder diameter (m) D
Dimensionless right eigenvector E
Dimensionless normal convective flux vector FN
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 (m2 sec-1) α'
Angle between normal and horizontal lines at cell face ζ
Dimensionless primary and secondary cell area Ω, Ω'
Kinematic viscosity (m2 sec-1) v
Time (sec) τ
Dimensional eigenvalue λ
Density (kg m-3) ρ
Roe-like weighted averaging parameter ω
Dimensionless temperature θ
Subscripts and superscripts
Normal, parallel N, P
Right, Left R, L
Wall W
Mean -
Drikakis, D. and P.K. Smolarkiewicz, 2001. On spurious vertical structures. J. Comput. Phys., 172: 309-325.
Direct Link  |  

Drikakis, D., P.A. Govatsos and D.E. Papatonis, 1994. A characteristic based method for incompressible flows. Int. J. Numerical Meth. Fluids, 19: 667-685.
Direct Link  |  

Jameson, A., 1995. Analysis and design of numerical schemes for gas dynamics 1: Artificial diffusion, upwind biasing, limiters and their effect on accuracy and multi-grid convergence. Int. J. Comp. Flu. Dyn., 4: 171-218.
Direct Link  |  

Kuehn, T.H. and R.J. Goldstein, 1980. Numerical solution to the Naveir-Stokes equations for laminar natural convection about a horizontal isothermal circular cylinder. Int. J. Heat Mass Transfer, 23: 971-980.
Direct Link  |  

Kwak, D., C. Kiris and C.S. Kim, 2005. Computational challenges of viscous incompressible flows. Comput. Fluids, 34: 283-299.
CrossRef  |  

Mallinger, F. and D. Drikakis, 2002. Instability in, three-dimensional, unsteady stenotic flows. Int. J. Heat Fluid Flow, 23: 657-663.
Direct Link  |  

Morgan, V.T., 1975. The Overall Convective Heat Transfer from Smooth Circular Cylinders. Academic Press, Inc., New York.

Neofytou, P., 2007. Revision of the characteristic-based scheme for incompressible flows. J. Comput. Phys., 222: 475-484.
CrossRef  |  Direct Link  |  

Roe, P.L., 1981. Approximate Riemann solvers, parameter vectors and difference schemes. J. Comput. Phys., 43: 357-372.
CrossRef  |  Direct Link  |  

Saitoh, T., T. Sajik and K. Maruhara, 1993. Bench mark solutions to natural convection heat transfer problems around a horizontal circular cylinder. Int. J. Heat Mass Transfer, 36: 1251-1259.
Direct Link  |  

Shapiro, E. and D. Drikakis, 2006. Non-conservative and conservative formulations of characteristics numerical reconstructions for incompressible flows. Int. J. Numer. Meth. Eng., 66: 1466-1482.
Direct Link  |  

Shapiro, E. and D. Drikakis, 2005. Artificial compressibility, characteristics-based schemes for variable density, incompressible, multi-species flows. Part I. Derivation of different formulations and constant density limit. J. Comput. Phys., 210: 548-607.
Direct Link  |  

Volpe, G., 1993. Performance of compressible flow codes at low mach numbers. AIAA. J., 31: 49-56.
Direct Link  |  

Wang, P., M. Kahawita and T. Naguyen, 1990. Numerical computational of the natural convection flow about a horizontal cylinder using splines. Numer. Heat Transfer, 17: 191-215.
Direct Link  |  

Zhao, Y. and B. Zhang, 2000. A high-order characteristics upwind FV method for incompressible flow and heat transfer simulation on unstructured grids. Comput. Meth. Applied Mech. Eng., 190: 733-756.
Direct Link  |  

©  2020 Science Alert. All Rights Reserved