Advanced
Convergence Characteristics of Upwind Method for Modified Artificial Compressibility Method
Convergence Characteristics of Upwind Method for Modified Artificial Compressibility Method
International Journal of Aeronautical and Space Sciences. 2011. Dec, 12(4): 318-330
Copyright ©2011, The Korean Society for Aeronautical Space Science
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/bync/ 3.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.
  • Received : July 15, 2011
  • Accepted : December 07, 2011
  • Published : December 30, 2011
Download
PDF
e-PUB
PubReader
PPT
Export by style
Share
Article
Author
Metrics
Cited by
TagCloud
About the Authors
Hyungro Lee
Seungsoo Lee
slee@inha.ac.kr

Abstract
This paper investigates the convergence characteristics of the modified artificial compressibility method proposed by Turkel. In particular, a focus is mode on the convergence characteristics due to variation of the preconditioning factor ( αu ) and the artificial compressibility ( β ) in conjunction with an upwind method. For the investigations, a code using the modified artificial compressibility is developed. The code solves the axisymmetric incompressible Reynolds averaged Navier-Stokes equations. The cell-centered finite volume method is used in conjunction with Roe’s approximate Riemann solver for the inviscid flux, and the central difference discretization is used for the viscous flux. Time marching is accomplished by the approximated factorizationalternate direction implicit method. In addition, Menter’s k-ω shear stress transport turbulence model is adopted for analysis of turbulent flows. Inviscid, laminar, and turbulent flows are solved to investigate the accuracy of solutions and convergence behavior in the modified artificial compressibility method. The possible reason for loss of robustness of the modified artificial compressibility method with αu >1.0 is given.
Keywords
1. Introduction
Steady solutions of fluid flows can be obtained by solving the Euler equations, or the Navier-Stokes equations, with time marching methods. The time marching methods, however, are not applicable to the incompressible Navier-Stokes equations because the continuity equation has no time derivative, unlike the momentum equations. To overcome the difficulty associated with the lack of the continuity equation, Chorin (1967) proposed an artificial compressibility method where an artificial time derivative term of pressure is added to the continuity equation with an artificial compressibility parameter ( β ). With the modification, various numerical techniques, originally developed for hyperbolic equations, can be applied to the incompressible Navier-Stokes equations.
After Chorin’s pioneering work, many researchers applied the artificial compressibility method to various incompressible flow analyses. Peyret and Taylor (1983), and Rahman and Siikonen (2008) solved the steady flow problems using the artificial compressibility method. Merkle and Athavale (1987), and Rogers and Kwak (1990) analyzed unsteady flows using the dual time stepping method with the artificial compressibility method. Turkel (1987) suggested a modification of the artificial compressibility method, where artificial time derivatives are added to the momentum equations as well as the continuity equation. This modification introduced an additional parameter ( αu ). The two arbitrary parameters ( αu and β ) play a significant role in the stability and convergence rates of solutions. Choice of values, however, requires considerable numerical experiments as well as users’ intuition. Some investigations of stability and convergence related to these parameters were conducted by Michelassi et al. (1996), Kiris et al. (2006), and Esfahanian and Akbarzadeh (2009). Turkel proposed using αu = 2.0, which makes the condition number equal to 1. As Malan et al. (2002a, b) pointed out, however, when αu = 2.0 the modified artificial compressibility method became less robust.
Most research concerned with the modified artificial compressibility method has so far been based on the central difference method with the 2nd order and the 4th order artificial dissipation. Only a few studies adopted upwind methods for the artificial compressibility method. Pan and Chakravarthy (1989) pointed out that the solution with an upwind method was no longer independent of the artifi cial compressibility parameter, and the built-in dissipation associated with the upwind method could contaminate the solution. However, they argued that accurate solutions could be obtained for a suitable choice of β .
The objective of this paper is to investigate convergence characteristics of the modified artificial compressibility method with an upwind method. To achieve this, an analysis code is developed based on Roe’s approximate Riemann solver. The following section contains the numerical method used in the paper. Also, Eigenstructure of the incompressible Navier-Stokes equations with the modified artificial compressibility method is given in details. Next, the accuracy and the convergence rates of various incompressible flows are compared. The convergence characteristics with combinations of the two parameters are also examined. In the numerical simulations, the Courant-Friedrichs-Lewy (CFL) number is kept constant in order to exclude eff ects of the CFL number on the convergence characteristics. Conclusions will be drawn from the numerical investigation.
2. Numerical Method
- 2.1 Governing equations
The axisymmetric Reynolds averaged Navier-Stokes (RANS) equations and two equation turbulence model equations for the incompressible flows are chosen as the governing equations. The modified artificial compressibility method of Turkel is adopted in order to use a time marching method. The non-dimensionalized governing equations are as follows:
Lager Image
where
Lager Image
and
Lager Image
are the inviscid vector and the viscous flux vector respectively. The source term,
Lager Image
is defined by:
Lager Image
where S is the source term due to the turbulence model and the second term is related to axisymmetry. Two dimensional problems can be solved with I = 0 while axisymmetric problems can be solved with I = 1. Here,
Lager Image
represents the outward normal vector of the computational cell. All of the vectors are defined as
Lager Image
where p, u, v, and T are the pressure, the axial velocity component, the radial velocity component, and the temperature respectively. Also, k , and ω represent the turbulence kinetic energy and the specific dissipation rate. Vn is the normal velocity component. The stress tensor are defined by
Lager Image
where v is the kinematic viscosity. Γ is the preconditioning matrix and is defined as follows:
Lager Image
In Eq. (5), the artifi cial compressibility, β , and the preconditioning factor, αu , are user defined parameters. By setting αu = 0.0, Eq. (1) becomes the original artifi cial compressibility method. Turkel (1992) suggested that a constant β be used for a good convergence rate.
- 2.2 Numerical schemes
The cell-centered fi nite volume method is applied to Eq. (1). The semi-discretized equation is found to be:
Lager Image
where the residual term is defined as
Lager Image
The total flux vector,
Lager Image
is defined as follows
Lager Image
In Eq. (8),
Lager Image
and
Lager Image
represent the inviscid and viscous flux vectors respectively. The viscous flux vectors are evaluated with the gradient theorem, which is equivalent to the central difference in the Cartesian coordinate system. Th e inviscid flux vectors are computed with Roe’s approximation Riemann solver (Roe, 1981). The numerical flux vector of Roe’s method is given by:
Lager Image
where Δ Q = QR ? QL , and
Lager Image
and
Lager Image
represent the flux computed using the right state, QR and left state, QL on each side of the computational cell faces. In addition, the monotone upstream-centered schemes for conservation laws (Van Leer, 1979) and van Albada’s limiter are applied in order to obtain higher-order spatial accuracy. The inviscid flux vector normal to the cell surface,
Lager Image
is defined by:
Lager Image
The Jacobian matrix, A Γ is defined by
Lager Image
The eigenvalues of the preconditioned system (1) are found to be:
Lager Image
where
Lager Image
Lager Image
Lager Image
From Eq. (12)-(15), αu and β are closely connected to wave speed, which has signifi cant impact on convergence rate, particularly αu = 2.0, λ 2, 3 become ∓ β so that they are independent of the flow velocity. Also, the condition number, |λ 2 / λ 3 |, becomes one. For this reason, Turkel proposed using αu = 2.0 for the convergence rate. The matrix | A Γ | represented in Eq. (9) is defined by:
Lager Image
where X Γ is modal matrix derived from the preconditioned system. The diagonal matrix, A Γ is given by:
Lager Image
The modal matrix is defined by:
Lager Image
Th us, the determinant of X Γ is:
Lager Image
Equation (19) indicates that the determinant can be zero when αu >1.0. This causes Roe’s numerical flux vector to become unbounded. The upwind method fails as the determinant of X Γ approaches zero. When αu = 2.0, β should be chosen to satisfy the following equation for the entire computational domain:
Lager Image
Turkel (1987) derived the same expression as Eq. (20) from the symmetrizability requirement of the system. When central difference methods are used, this manifests as loss of robustness, which was observed by Malan et al. (2002a, b). Most boundary condition methods use characteristic information to determine the solution variables along the computational boundaries. When the determinant of X Γ becomes zero, therefore, the boundary condition methods fail. When αu = 2.0, therefore, the value of β should be large enough to satisfy Eq. (20) for the entire computational domain.
For the time-integration method, the approximated factorization-alternate direction implicit method (Beam and Warming, 1982) is employed. Detailed information of time discretization can be found in Beam and Warming’s study.
- 2.3 Turbulence model
Menter’s k −ω shear stress transport (SST) turbulence model (Menter, 1994) is used to compute the turbulent flows. Menter’s model uses the turbulent kinetic energy ( k ) and the specific dissipation rate (ω) as turbulence variables. The eddy viscosity is evaluated by Eq. (21):
Lager Image
where the blending function, F 2 is defined by:
Lager Image
Other details of the parameters and source term of the k −ω SST turbulence model can be found in Menter’s study.
3. Computational Results
- 3.1 Inviscid flow around a circular cylinder
As the first computational example, an inviscid flow around a circular cylinder is calculated. This test case is chosen to verify the accuracy of the upwind method and the grid convergence of the modified artifficial compressibility method in conjunction with the upwind method. The
Lager Image
Grid convergence with variation of β (αu = 0.0).
simulations are performed,increasing grid size in variation of β . The grid sizes used are listed in Table 1 . The errors are defined as:
Lager Image
where qapp and qexact are the numerical and exact solutions respectively. Here, Vi and V are the cell volume and the total volume respectively, of the computational domain. Figure 1 shows the result of the grid convergence test. The values of β in this figure cover the entire range where the converged solution could be obtained. Figure 2 presents the computed surface pressure coefficient distributions as well as the potential solution. All solutions with different β ’s overlap with each other as well as with the exact solution. The numerical solutions converge into the exact solutions as the mesh size decreases for the range of β . In other words,
Grid size used for circular cylinder
Lager Image
Grid size used for circular cylinder
the accuracy of the solution is guaranteed when a converged solution is obtained (Lax equivalence theorem). Moreover, Roe’s method satisfies the consistency when applied to the artificial compressibility method. The numerical dissipation for the first order upwind method for the Cartesian coordinate system can be expressed as:
Lager Image
where A Γ and B Γ are the Jacobian matrices of the Cartesian flux vectors, E and F . That is:
Lager Image
Pressure coefficients for a circular cylinder (αu = 0.0)
Lager Image
Pressure coefficients for NACA 0012 problem.
Lager Image
From Eq. (24), it is clear that the numerical dissipation vanishes as Δ x , Δ y → 0. Similarly, the numerical dissipation of high order methods vanishes as the mesh size decreases.
Lager Image
Convergence histories for NACA 0012 problem (AOA = 0.0 de-gree).
- 3.2 Inviscid flow around a NACA0012 airfoil
To compare convergence rate with the combinations of αu and β , an inviscid flow around an NACA0012 airfoil is solved. Simulations are performed with a C-type grid of 257× 65. The CFL number is set to 5.0 for all simulations. Figure 3 shows the surface pressure coefficient distributions at angle of attack (AOA) of 3.0 degrees. From Fig. 3 , it is clear that the pressure coefficients of the converged solutions are highly accurate, regardless of αu . Also, a result of a panel code is presented for comparison, indicating the accuracy of the artificial compressibility methods.
The convergence histories with various αu and β combinations are plotted in Fig. 4 . All computations are performed with an AOA of 0.0 degree. In the figure, the horizontal lines indicate that the computations fail. From
Lager Image
Determinant contours of modal matix.
the figure, it is clear that there exists an optimal value for convergence. For this case, β of 5.0 gives the best convergence. When Figs. 4 b and c are compared with Fig. 4 a, the modified artificial compressibility method has better convergence characteristics than the original method with the value of β greater than 1. As predicted with Eq. (20), the modified artificial compressibility method with αu = 2.0 suffers loss of robustness for a small value of β . Figure 5 depicts the contour plots of the determinant of the modal matrix with various combinations of αu and β . There are no singular regions in the entire computational domain for β = 0.5, and β = 0.5 with αu = 0.0. However, singular regions are found near the leading edge of the airfoil, with αu = 2.0, and β = 0.5. Figure 6 shows the convergence histories at three angles of attack; 0.0, 1.25, and 3 degrees. All the computations are made with
Lager Image
Convergence histories with variation of angle of attack (AOA) (αu = 5.0).
the optimal values of β = 5.0. As can be seen from the figure, the modified artificial compressibility method with αu = 2.0 performs best when compared with other αu .
- 3.3 Laminar flow around a sphere
The second example of examinining convergence characteristics is a laminar flow around a sphere. It is wellknown that flow around a sphere is steady and axisymmetric when Re<220. However, the flow becomes unsteady and
Lager Image
Streamlines around a sphere with variation of Reynolds num-bers (αu = 0.0, β = 5.0).
three-dimensional when Re>220. In this paper, numerical simulations are performed with Reynolds numbers less than 220. An O-type computational grid of 129× 57 is used for the computations. The CFL number is set to 3.0 for all computations. Figure 7 depicts the streamlines around the sphere at four different Reynolds numbers; 30, 50, 100, and 150. The size of the separation bubble behind the sphere increases commensurately to the Reynolds number. The drag coefficients computed with different values of αu are presented in Fig. 8 . The drag coefficient measured by Roos and Willmarth (1971) and the drag coefficients computed by Mittal (1999) and Sheard et al. (2003) are also presented for comparison. The results of the present study match well with results of previous research.
Figure 9 depicts convergence histories of the sphere problem with various combinations of αu and β . The Reynolds number of all the results in this figure is 100. Unlike the previous example, the best convergence rate is achieved with the different values of β for each values of αu . It confirms the conclusions from the previous example that there exists an optimal value of β for best convergence rate and that the modified artificial compressibility method with αu = 2.0 has the best convergence characteristics at the optimal value of the compressibility factor. In particular, Fig. 9 c indicates that the computations with α u = 2.0 fail with values of β that are less than 1.0 as in the previous problem. Figure 10 depicts the convergence histories at various Reynolds numbers. This figure confirms that the modified method, with αu = 2.0, yields the highest convergence rate while the original method gives the lowest convergence rate.
Lager Image
Drag coefficients of a sphere versus the Reynods numbers (β = 5.0).
Lager Image
Convergence histories for a sphere problem (Re = 100).
Lager Image
Convergence histories with variation of the Reynolds number (β = 5.0).
Lager Image
Velocity profiles at several axial locations of plate (αu = 0.0, β = 5.0)
- 3.4 Turbulent flat plate flow
A turbulent flow over a flat plate flow is selected to study the convergence characteristics of the modified artificial compressibility method for turbulent flows. The size of the computational grid is 111×81. The first grid spacing above the wall satisfies the requirement of y + <1 so that the turbulent boundary layer can be correctly resolved. For all the computations, the CFL number of 5.0, the preconditioning factor, αu of 0.0, and the artificial compressibility factor, β of 0.5 are used. Although the results with the original method ( αu = 0.0) are only represented here, there is little difference in the results with αu = 0.0, and αu = 1.0. However, the converged solution could not be obtained with αu = 2.0 with β = 0.5, as in laminar computations. Figure 11 exhibits the turbulent boundary layer profiles at three different locations: Re x = 2.70×10 6 , 7.62×10 6 , and 1.03×10 7 . The experimental data of Kline et al. (1969), the numerical result of Ryu et al. (2006), and the law of wall are also presented for comparison. The result of Ryu et al. (2006) was obtained with a compressible RANS solver (QT2D) and Menter’s k −ω SST turbulence model. Excellent agreement can be observed from the figure. Figure 12 compares the computed skin friction coefficient with Ryu et al. (2006)’s result, the experimental data, and the 1/7th law and Schlichting’s formula (1979). It is noted that the computation is performed with the assumption of fullyturbulent flow.
In Fig. 13 , convergence histories are compared for the combinations of αu and β . For the original artificial compressibility method, the convergence rates with β = 0.5 and 1.0 are relatively high. A similar trend with αu = 1.0 can be observed from the figure. It is interesting to notice that
Lager Image
Skin friction coefficient versus axial location of plate (αu = 0.0, β = 5.0).
Lager Image
Convergence histories for flat-plate problem.
the optimal value of β , to achieve good convergence, is less than that of the inviscid or laminar problems. This difference is thought to come from high aspect ratio of computational cells near the solid wall, for the resolution of the turbulent boundary layer. The stability characteristics of the artificial compressibility method associated with high aspect ratio cells would be an interesting research topic. Figure 13 c shows that converged solutions cannot be not obtained with β <1.0 because of the zero determinant problem mentioned earlier. It is especially noticeable that the convergence with β = 5.0 suddenly stalls after 13,000 iterations. The result also confirms that the modified artificial compressibility method with αu = 2.0 becomes less robust.
- 3.5 Turbulent flow passing over a bump
As the final example, the turbulent flow over a bump is selected. This example is one validation case of turbulence
Lager Image
Comparison of pressure/skin friction coefficients (β = 5.0).
Lager Image
Convergence histories for a bump problem.
modeling resource at the Langley Research Center website. A grid of 353×161 is downloaded from the website and used for the computation. The first grid spacing above the wall satisfies the grid resolution requirement of y + <1. The Reynolds number based on the length of the bump is 3.0×10 6 . The CFL number of 5.0 is used for the simulations. The results of the numerical simulations are compared with CFL3D (Langley Research Center, 2011) and QT2D, where the freestream Mach number is 0.2. The turbulence model in CFL3D is k −ω SST-V, which is a modified version of the k −ω SST model, while QT2D uses the k −ω SST model. Figure 14 compares the pressure coefficient and the skin friction coefficient with the other results. The results match well with each other.
Figure 15 represents convergence histories with combinations of β and αu . The figure shows that there exists an optimal value of β for good convergence, and the modified artificial compressibility method with αu = 2.0 are the least robust.
4. Conclusions
In this study, an incompressible RANS code, which uses the modified artificial compressibility method of Turkel and the upwind method of Roe, was developed; the convergence characteristics were studied for various incompressible flow problems. It is confirmed numerically and analytically that the accuracy of the solution can be assured with the upwind method and the artificial compressibility method. The modified artificial compressibility method has superior convergence characteristics compared to the original method. However, the modified artificial compressibility method with αu = 2.0 can fail when used with a low value of β for both the upwind method and the central difference method. It is shown that this loss of robustness comes from the fact that the determinant of the modal matrix can be zero when αu >1.0.
References
Beam R. M. , Warming R. F. (1982) Implicit Numerical Methods for the Compressible Navier-Stokes and Euler Equations. Rhode-Saint-Genese Von Karman Institute for Fluid Dynamics. Belgium
Chorin A. J. (1967) A numerical method for solving incompressible viscous flow problems. Journal of Computational Physics 2 12 - 26    DOI : 10.1016/0021-9991(67)90037-X
Esfahanian V. , Akbarzadeh P. (2009) Advanced investigation on design criteria of a robust, artificial compressibility and local preconditioning method for solving the inviscid incompressible flows. Proceeding of the Third International Conference on Modeling Simulation, and Applied Optimization Sharjah, UAE.
Kline S. J. , Coles D. E. , Hirst E. A. US Air Force. (1969) Office of Scientific Research. Thermosciences Division. Computation of Turbulent Boundary Layers--1968 AFOSR-IFP-Stanford Conference Proceedings. Stanford: Stanford University.
Kiris C. , Housman J. , Kwak D. , Groth C. , Zingg D. W. (2006) Comparison of artificial compressibility methods. Computational Fluid Dynamics 2004. Springer Berlin. Heidelberg, Germany 475 - 480
Langley Research Center (2011) Turbfulence Modeling Resource. http://turbmodels.larc.nasa.gov/
Malan A. G. , Lewis R. W. , Nithiarasu P. (2002a) An improved unsteady, unstructured, artificial compressibility, finite volume scheme for viscous incompressible flows. Part I. Theory and implementation. International Journal for Numerical Methods in Engineering 54 695 - 714    DOI : 10.1002/nme.447
Malan A. G. , Lewis R. W. , Nithiarasu P. (2002b) An improved unsteady, unstructured, artificial compressibility, finite volume scheme for viscous incompressible flows. Part II. Application. International Journal for Numerical Methods in Engineering 54 715 - 729    DOI : 10.1002/nme.443
Menter F. R. (1994) Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal 32 1598 - 1605    DOI : 10.2514/3.12149
Merkle C. L. , Athavale M. (1987) Time-accurate unsteady incompressible flow algorithms based on artificial compressibility. AIAA 8th Computational Fluid Dynamics Conference Honolulu, Hawaii. AIAA Paper No. 87-1137.
Michelassi V. , Migliorini F. , Martelli F. (1996) Preconditioned scalar approximate factorization method for incompressible fluid flows. International Journal of Computational Fluid Dynamics 7 311 - 325    DOI : 10.1080/10618569608940768
Mittal R. (1999) A Fourier-Chebyshev spectral collocation method for simulating flow past spheres and spheroids. International Journal for Numerical Methods in Fluids 30 921 - 937
Pan D. , Chakravarthy S. (1989) Unified formulation for incompressible flows. AIAA 27th Aerospace Science Meeting Reno, NV. AIAA Paper No. 89-0122.
Peyret R. , Taylor T. D. (1983) Computational Methods for Fluid Flow. Springer-Verlag. New York 358 -
Rahman M. M. , Siikonen T. (2008) An artificial compressibility method for viscous incompressible and low Mach number flows. International Journal for Numerical Methods in Engineering 75 1320 - 1340    DOI : 10.1002/nme.2302
Roe P. L. (1981) Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics 43 357 - 372    DOI : 10.1016/0021-9991(81)90128-5
Rogers S. E. , Kwak D. (1990) Upwind differencing scheme for the time-accurate incompressible Navier-Stokes equations. AIAA Journal 28 253 - 262    DOI : 10.2514/3.10382
Roos F. W. , Willmarth W. W. (1971) Some experimental results on sphere and disk drag. AIAA Journal 9 285 - 291    DOI : 10.2514/3.6164
Ryu S. , Lee S. , Kim B. (2006) A study of local preconditioning method for compressible low speed flows. Journal of the Korea Institute of Military Science and Technology 9 152 - 160
Schlichting H. (1979) Boundary-Layer Theory. 7th ed. McGraw-Hill. New York
Sheard G. J. , Thompson M. C. , Hourigan K. (2003) From spheres to circular cylinders: the stability and flow structures of bluff ring wakes. Journal of Fluid Mechanics 492 147 - 180    DOI : 10.1017/S0022112003005512
Turkel E. (1987) Preconditioned methods for solving the incompressible and low speed compressible equations. Journal of Computational Physics 72 277 - 298    DOI : 10.1016/0021-9991(87)90084-2
Turkel E. (1992) Review of Preconditioning Methods for Fluid Dynamics. Institute for Computer Applications in Science and Engineering. Hampton ICASE Report No. 92-74.
Van Leer B. (1979) Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method. Journal of Computational Physics 32 101 - 136    DOI : 10.1016/0021-9991(79)90145-1