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.

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:

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 with I = 0 while axisymmetric problems can be solved with I = 1. Here,

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

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. V_n is the normal velocity component. The stress tensor are defined by

where 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.

The cell-centered fi nite volume method is applied to Eq. (1). The semi-discretized 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, and

and

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 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,

is defined by:

The Jacobian matrix, A_Γ is defined by

The 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, |λ₂ / λ₃|, 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 of X_Γ 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 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.

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₂ is defined by:

Other details of the parameters and source term of the k？ω SST turbulence model can be found in Menter’s study.

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 and q_exact are the numerical and exact solutions respectively. Here, V_i 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:

where A_Γ and B_Γ are the Jacobian matrices of the Cartesian flux vectors, E and F. 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.

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.

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.

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⁶, 7.62×10⁶, and 1.03×10⁷. 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.

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⁶. 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.

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.