Convergence Characteristics of Upwind Method for Modified Artificial Compressibility Method
 Author: Lee Hyungro, Lee Seungsoo
 Organization: Lee Hyungro; Lee Seungsoo
 Publish: International Journal Aeronautical and Space Sciences Volume 12, Issue4, p318~330, 30 Dec 2011

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 NavierStokes equations. The cellcentered 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’skω 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 NavierStokes equations , Upwind method , Artificial compressibility method

1. Introduction
Steady solutions of fluid flows can be obtained by solving the Euler equations, or the NavierStokes equations, with time marching methods. The time marching methods, however, are not applicable to the incompressible NavierStokes 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 NavierStokes 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 builtin 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 NavierStokes 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 CourantFriedrichsLewy (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 NavierStokes (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 nondimensionalized governing equations are as follows:
where
and
are the inviscid vector and the viscous flux vector respectively. The source term,
is defined by:
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 withI = 0 while axisymmetric problems can be solved withI = 1. Here,represents the outward normal vector of the computational cell. All of the vectors are defined as
where
p, u, v, andT 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.V_{n} is the normal velocity component. The stress tensor are defined bywhere
v is the kinematic viscosity. Γ is the preconditioning matrix and is defined as follows: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 cellcentered fi nite volume method is applied to Eq. (1). The semidiscretized equation is found to be:
where the residual term is defined as
The total flux vector,
is defined as follows
In Eq. (8),
and
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:
where Δ
Q =Q_{R} ?Q_{L} , andand
represent the flux computed using the right state,
Q_{R} and left state,Q_{L} on each side of the computational cell faces. In addition, the monotone upstreamcentered schemes for conservation laws (Van Leer, 1979) and van Albada’s limiter are applied in order to obtain higherorder spatial accuracy. The inviscid flux vector normal to the cell surface,is defined by:
The Jacobian matrix,
A _{Γ} is defined byThe eigenvalues of the preconditioned system (1) are found to be:
where
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:where
X _{Γ} is modal matrix derived from the preconditioned system. The diagonal matrix,A _{Γ} is given by:The modal matrix is defined by:
Th us, the determinant of
X _{Γ} is: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 ofX _{Γ} approaches zero. Whenα_{u} = 2.0,β should be chosen to satisfy the following equation for the entire computational domain: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 timeintegration method, the approximated factorizationalternate 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):where the blending function,
F _{2} is defined by: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:where
q_{app} andq_{exact} are the numerical and exact solutions respectively. Here,V_{i} andV 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:
where
A _{Γ} andB _{Γ} are the Jacobian matrices of the Cartesian flux vectors,E andF . That is: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 Ctype 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. Fromthe 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 withthe 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
threedimensional when Re>220. In this paper, numerical simulations are performed with Reynolds numbers less than 220. An Otype 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: 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’sk ？ω 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 thatthe 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×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 isk ？ω SSTV, which is a modified version of thek ？ω SST model, while QT2D uses thek ？ω 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.

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]

[Fig. 1.] Grid convergence with variation of β (αu = 0.0).

[Table 1.] Grid size used for circular cylinder

[Fig. 2.] Pressure coefficients for a circular cylinder (αu = 0.0)

[Fig. 3.] Pressure coefficients for NACA 0012 problem.

[Fig. 4.] Convergence histories for NACA 0012 problem (AOA = 0.0 degree).

[Fig. 5.] Determinant contours of modal matix.

[Fig. 6.] Convergence histories with variation of angle of attack (AOA) (αu = 5.0).

[Fig. 7.] Streamlines around a sphere with variation of Reynolds numbers (αu = 0.0, β = 5.0).

[Fig. 8.] Drag coefficients of a sphere versus the Reynods numbers (β = 5.0).

[Fig. 9.] Convergence histories for a sphere problem (Re = 100).

[Fig. 10.] 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)

[Fig. 12.] Skin friction coefficient versus axial location of plate (αu = 0.0, β = 5.0).

[Fig. 13.] Convergence histories for flatplate problem.

[Fig. 14.] Comparison of pressure/skin friction coefficients (β = 5.0).

[Fig. 15.] Convergence histories for a bump problem.