Convergence Characteristics of Upwind Method for Modified Artificial Compressibility Method

  • cc icon
  • 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.


  • KEYWORD

    Computational fluid dynamics , Incompressible Navier-Stokes equations , Upwind method , Artificial compressibility method

  • 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:

    image

    where

    image

    and

    image

    are the inviscid vector and the viscous flux vector respectively. The source term,

    image

    is defined by:

    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,

    image

    represents the outward normal vector of the computational cell. All of the vectors are defined as

    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

    image

    where v is the kinematic viscosity. Γ is the preconditioning matrix and is defined as follows:

    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:

    image

    where the residual term is defined as

    image

    The total flux vector,

    image

    is defined as follows

    image

    In Eq. (8),

    image

    and

    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:

    image

    where ΔQ = QR?QL, and

    image

    and

    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,

    image

    is defined by:

    image

    The Jacobian matrix, AΓ is defined by

    image

    The eigenvalues of the preconditioned system (1) are found to be:

    image

    where

    image
    image
    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:

    image

    where XΓ is modal matrix derived from the preconditioned system. The diagonal matrix, AΓ is given by:

    image

    The modal matrix is defined by:

    image

    Th us, the determinant of XΓ is:

    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:

    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):

    image

    where the blending function, F2 is defined by:

    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

    simulations are performed,increasing grid size in variation of β. The grid sizes used are listed in Table 1. The errors are defined as:

    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,

    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:

    image

    where AΓ and BΓ are the Jacobian matrices of the Cartesian flux vectors, E and F. That is:

    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.

       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

    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. 4b and c are compared with Fig. 4a, 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

    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

    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. 9c 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.

       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: Rex = 2.70×106, 7.62×106, and 1.03×107. 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

    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 13c 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

    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×106. 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.

  • 1. Beam R. M., Warming R. F. (1982) Implicit Numerical Methods for the Compressible Navier-Stokes and Euler Equations. Rhode-Saint-Genese google
  • 2. Chorin A. J. (1967) A numerical method for solving incompressible viscous flow problems. [Journal of Computational Physics] Vol.2 P.12-26 google doi
  • 3. 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] google
  • 4. Kline S. J., Coles D. E., Hirst E. A. (1969) Office of Scientific Research. [Computation of Turbulent Boundary Layers--1968 AFOSR-IFP-Stanford Conference Proceedings.] google
  • 5. Kiris C., Housman J., Kwak D., Groth C., Zingg D. W. (2006) Comparison of artificial compressibility methods. Computational Fluid Dynamics 2004. P.475-480 google
  • 6. (2011) Turbfulence Modeling Resource. google
  • 7. 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] Vol.54 P.695-714 google doi
  • 8. 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] Vol.54 P.715-729 google doi
  • 9. Menter F. R. (1994) Two-equation eddy-viscosity turbulence models for engineering applications. [AIAA Journal] Vol.32 P.1598-1605 google doi
  • 10. Merkle C. L., Athavale M. (1987) Time-accurate unsteady incompressible flow algorithms based on artificial compressibility. [AIAA 8th Computational Fluid Dynamics Conference] google
  • 11. Michelassi V., Migliorini F., Martelli F. (1996) Preconditioned scalar approximate factorization method for incompressible fluid flows. [International Journal of Computational Fluid Dynamics] Vol.7 P.311-325 google doi
  • 12. Mittal R. (1999) A Fourier-Chebyshev spectral collocation method for simulating flow past spheres and spheroids. [International Journal for Numerical Methods in Fluids] Vol.30 P.921-937 google
  • 13. Pan D., Chakravarthy S. (1989) Unified formulation for incompressible flows. [AIAA 27th Aerospace Science Meeting] google
  • 14. Peyret R., Taylor T. D. (1983) Computational Methods for Fluid Flow. P.358 google
  • 15. 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] Vol.75 P.1320-1340 google doi
  • 16. Roe P. L. (1981) Approximate Riemann solvers, parameter vectors, and difference schemes. [Journal of Computational Physics] Vol.43 P.357-372 google doi
  • 17. Rogers S. E., Kwak D. (1990) Upwind differencing scheme for the time-accurate incompressible Navier-Stokes equations. [AIAA Journal] Vol.28 P.253-262 google doi
  • 18. Roos F. W., Willmarth W. W. (1971) Some experimental results on sphere and disk drag. [AIAA Journal] Vol.9 P.285-291 google doi
  • 19. 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] Vol.9 P.152-160 google
  • 20. Schlichting H. (1979) Boundary-Layer Theory. google
  • 21. 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] Vol.492 P.147-180 google doi
  • 22. Turkel E. (1987) Preconditioned methods for solving the incompressible and low speed compressible equations. [Journal of Computational Physics] Vol.72 P.277-298 google doi
  • 23. Turkel E. (1992) Review of Preconditioning Methods for Fluid Dynamics. google
  • 24. Van Leer B. (1979) Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method. [Journal of Computational Physics] Vol.32 P.101-136 google doi
  • [Fig. 1.] Grid convergence with variation of β (αu = 0.0).
    Grid convergence with variation of β (αu = 0.0).
  • [Table 1.] Grid size used for circular cylinder
    Grid size used for circular cylinder
  • [Fig. 2.] Pressure coefficients for a circular cylinder (αu = 0.0)
    Pressure coefficients for a circular cylinder (αu = 0.0)
  • [Fig. 3.] Pressure coefficients for NACA 0012 problem.
    Pressure coefficients for NACA 0012 problem.
  • [Fig. 4.] Convergence histories for NACA 0012 problem (AOA = 0.0 de-gree).
    Convergence histories for NACA 0012 problem (AOA = 0.0 de-gree).
  • [Fig. 5.] Determinant contours of modal matix.
    Determinant contours of modal matix.
  • [Fig. 6.] Convergence histories with variation of angle of attack (AOA) (αu = 5.0).
    Convergence histories with variation of angle of attack (AOA) (αu = 5.0).
  • [Fig. 7.] Streamlines around a sphere with variation of Reynolds num-bers (αu = 0.0, β = 5.0).
    Streamlines around a sphere with variation of Reynolds num-bers (αu = 0.0, β = 5.0).
  • [Fig. 8.] Drag coefficients of a sphere versus the Reynods numbers (β = 5.0).
    Drag coefficients of a sphere versus the Reynods numbers (β = 5.0).
  • [Fig. 9.] Convergence histories for a sphere problem (Re = 100).
    Convergence histories for a sphere problem (Re = 100).
  • [Fig. 10.] Convergence histories with variation of the Reynolds number (β = 5.0).
    Convergence histories with variation of the Reynolds number (β = 5.0).
  • [Fig. 11.] Velocity profiles at several axial locations of plate (αu = 0.0, β = 5.0)
    Velocity profiles at several axial locations of plate (αu = 0.0, β = 5.0)
  • [Fig. 12.] Skin friction coefficient versus axial location of plate (αu = 0.0, β = 5.0).
    Skin friction coefficient versus axial location of plate (αu = 0.0, β = 5.0).
  • [Fig. 13.] Convergence histories for flat-plate problem.
    Convergence histories for flat-plate problem.
  • [Fig. 14.] Comparison of pressure/skin friction coefficients (β = 5.0).
    Comparison of pressure/skin friction coefficients (β = 5.0).
  • [Fig. 15.] Convergence histories for a bump problem.
    Convergence histories for a bump problem.