Abstract: 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.
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 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.
GOVERNING EQUATIONS
The non-dimensional form of Navier-Stokes equations in finite-volume form with the ACM and the Boussinesq approximation is given by:
| (1) |
| (2) |
| (3) |
| (4) |
At which for each cell interface, ζ, uN, uP, ΔL,
FN, a,
| 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
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:
| (5) |
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:
| (6) |
The eigenvalues and corresponding eigenvectors of Jacobian matrix A are:
| (7) |
| (8) |
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 (
| (9) |
| (10) |
| (11) |
To evaluate virtual acoustic wave strength matrix αK, following system of equations is solved:
| (12) |
| (13) |
Finally, the normal convective flux at any cell side is obtained as:
| (14) |
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:
| (15) |
The viscous and thermal conduction terms were discretized using the secondary cells as shown in Fig. 2:
| (16) |
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:
| (17) |
| Fig. 2: | Schematic representation of primary and secondary cells to evaluate viscous terms |
| (18) |
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):
| (19) |
| (20) |
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:
| (21) |
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 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 |
CONCLUSION
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.
NOMENCLATURE
| 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 |
| 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 | θ |
| Normal, parallel | N, P |
| Environment | ∞ |
| Right, Left | R, L |
| Wall | W |
| Mean | - |