Effects of Numerical Methods on the Calculation of Developing Plane Channel Flow

Ronald M. C. So

Mechanical Engineering Department, The Hong Kong Polytechnic University, Hong Kong, China.

DOI: 10.4236/jamp.2022.106142                

Abstract

In wall-bounded turbulent flow calculations, the past focus has been directed to the modelling of the Reynolds-stress gradients. Not much attention has been paid to the effects of the numerical methods used to calculate these terms and the modelled equations. Discrepancies between model calculations and measurements are quite often attributed to incorrect modelling, while the suitability and accuracy of the numerical methods used are seldom scrutinized. Instead, alternate near-wall and Reynolds-stress models are proposed to remedy the incorrect turbulent flow calculations. On the other hand, if care is not taken in the numerical treatment of the Reynolds-stress gradient terms, physically unrealistic results and solution instability could occur. Previous studies by the author and his collaborators on the effects of numerical methods have shown that some of the more commonly used numerical methods could enhance numerical stability in the solution procedure but would introduce considerable inaccuracy to the results. The flow cases chosen to demonstrate these inaccuracies are a backstep flow and flow in a square duct, where flow complexities are present. The current investigation attempts to show that the above-mentioned effects of numerical methods could also occur in the calculation of a developing plane channel flow, where flow complexities are absent. In addition, this study shows that the results thus obtained lead to a predicted skin friction coefficient that is influenced more by the numerical method used than by the turbulence model invoked. Together, these results show that numerical treatment of the Reynolds-stress gradients in the equations play an important role, even for a developing plane channel flow.

Keywords

Numerical Representation of Pressure-Gradient Terms, Numerical Representation of Stress-Gradient Terms, Numerical Grid Effects, Numerical Methods Effects, Calculation of Plane Channel Flow

1. Introduction

Regardless of whether laminar or turbulent flows are considered, pressure and stress gradient terms are always present in the Navier-Stokes or Reynolds equations. Therefore, using a finite difference scheme to solve these equations, difficulties associated with numerically resolving the pressure-velocity coupling are well known (Harlow and Welch (1965) [1]; Patankar (1980) [2]; and van Doormal and Raithby (1984) [3]). When central differences are used on an ordinary grid system to approximate the pressure gradient terms, the converged results reveal the presence of oscillatory pressure and velocity fields. This problem was recognized by Harlow and Welch (1965) [1]. In order to prevent these oscillations from occurring, they suggested calculating the velocity variables at grid locations staggered with respect to those of the pressure and all other scalar variables. This staggered grid arrangement was later adopted by Patankar and Spalding (1972) [4] to solve three-dimensional parabolic flows using a finite volume scheme and incorporated into the SIMPLE algorithm of Patankar (1980) [2]. In the staggered grid arrangement, the velocity is driven by the pressure variables adjacent to it. Thus, a strong linkage between the velocity and pressure field is created. On the other hand, when all the variables are calculated at the same grid location, the velocities at grid point P are affected by pressure at every other grid location except that at P. Consequently, a physically incorrect pressure field could result during each iteration of the governing equations, and the final converged solution could give rise to an oscillatory pressure field. This phenomenon is commonly known as the checkerboard problem (Patankar (1980) [2]).

When turbulent flows are solved using Reynolds-stress models, a similar problem exists due to the presence of the Reynolds-stress gradient terms in the Reynolds equations. Like the pressure variable, the Reynolds stresses appear as first-order differentials in the coordinate directions. If central differencing is used to approximate them, a problem similar to that due to velocity-pressure decoupling occurs as explained above. Again, this leads to an oscillatory solution field, and usually, the iteration procedure diverges. This problem was recognized by Pope and Whitelaw (1976) [5]. Drawing on the pressure-velocity decoupling analogy, they suggested calculating the turbulent shear stresses, $\overline{uv}$, $\overline{uw}$, and $\overline{vw}$, at control volumes located halfway between the control volumes for the velocity variables. Here, u, v and w are the components of the fluctuating velocity along x, y, and z coordinates, respectively. In this way, they were able to create a linkage between the Reynolds stresses and their associated strain rates. Without this remedy, they were unable to obtain converged solutions.

The reasoning behind staggering the control volumes for the Reynolds shear stresses can be explained by considering the Boussinesq approximation assumed for the Reynolds stress and the mean strain-rate tensor, which is usually written as

$-\overline{u_i{u}_j}=2{\nu }_t{S}_{ij}-\frac{2}{3}k{\delta }_{ij}$. (1)

In Equation (1), $-\overline{u_i{u}_j}$ is the Reynolds-stress tensor, $S_{ij}=\left(\partial U_i/\partial x_j+\partial U_j/\partial x_i\right)/2$ is the mean strain-rate tensor, ${\nu }_t$ is an isotropic eddy viscosity, U_i, u_i and x _i are the ith components of the mean and fluctuating velocity, and coordinates, respectively, and k is the turbulent kinetic energy. According to this expression, $\overline{u_i{u}_j}$ is driven by S_ij. Therefore, in order to provide the physical coupling, $\overline{uv}$ needs to be calculated at control volumes halfway between the control volumes for U and V, $\overline{uw}$ between U and W, and $\overline{vw}$ between V and W. There is no need to stagger the normal stresses because, if the velocities are computed at staggered grid locations, $\overline{u^2}$, $\overline{v^2}$, and $\overline{w^2}$ are driven by the differences between the neighboring velocities U, V, and W, respectively, as desired.

When turbulent flows are solved using Reynolds-stress models, the problem associated with the velocity-pressure decoupling is compounded by the velocity-Reynolds-stress decoupling; both of which need to be addressed. However, if the velocity and the Reynolds shear stress variables are calculated at staggered control volumes as suggested by Harlow and Welch (1965) [1] and Pope and Whitelaw (1976) [5], seven different control volumes need to be specified within the code; three for the velocity variables, three for the Reynolds shear stresses, and one common control volume for all the other scalar variables, such as the pressure and the Reynolds normal stresses. This, in turn, could make the implementation of the Reynolds-stress models into a computer code a very difficult task, especially if a generalized, curvilinear coordinate system is used.

Furthermore, the absence of turbulent diffusion terms in the mean momentum equations renders the solution procedure less stable. This problem is usually avoided by the introduction of a turbulent kinematic viscosity, ν_t, when two-equation models are used. However, since Reynolds-stress models are anisotropic models, this approach is not suitable. Therefore, one has to resort to other means to extract diffusion like terms and introduce them into the momentum equations. The two issues of velocity-Reynolds-stress decoupling and turbulent diffusion are independent of each other and may introduce considerable errors to the calculation of complex turbulent flows if not handled properly. These two issues were examined by Aksoy et al. (1995) [6] in their study of a square duct flow, and by Aksoy and So (1995) [7] in their study of a backstep flow. In these studies, they used four different numerical methods (listed in Table 1) to evaluate the stress-gradient terms in the governing equations. Among the four numerical methods tested, not one could give consistently stable calculations with fairly accurate results for different flow types. Their performance or lack thereof could be attributed to turbulence effects generated by flow and geometric complexity. In order to gain further understanding of this insufficiency, it is necessary to scrutinize the performance of these numerical methods used to evaluate the pressure- and stress-gradient terms in a simple flow, where turbulence created by complex geometry and high shear are essentially absent. The present objective is to investigate whether the four different numerical methods used to treat the pressure- and stress-gradient terms in the governing

Table 1. Designation of different numerical methods.

equations would also give rise to similar numerical errors in the calculation of a simple developing plane channel flow.

To enhance comprehension, the present paper is organized into the following sections: 2) Objective, 3) Numerical Methods and Their Characteristics, 4) Near-Wall Modelling, 6) Numerical Treatment of the Stress-Gradient Terms, 7) Numerical Solution, 8) Results and Discussion, and 9) Conclusion.

2. Objective

The objective of the present study is focused on studying the effect of velocity-pressure and velocity-Reynolds-stress decoupling on the calculation results. Therefore, the present investigation limits the examination to the effects of numerical methods on the velocity-Reynolds-stress-gradient problem, and its associated effects on the near-wall Reynolds-stress modelling of simple turbulent flows. In particular, the ability to correctly replicate the development of the skin friction coefficient C_f is compared. Therefore, the nature of the present investigation is different from those where the emphasis is concentrated on the more general effects of physical and numerical influences on turbulent flow calculations. Some relevant studies on the more general physical and numerical influences are represented by Launder and Spalding (1974) [8], Gunzburger and Labovsky (2012) [9], Junqueire-Juniot et al. (2013) [10], Ahsan (2014) [11], Nguyen and Kempf (2017) [12], and Shahjada Tarafder and Al Mursaline (2019) [13].

The developing plane channel flow is free of turbulence complexities created by rapid change of flow geometry and the high shear thus created. Therefore, differences in the calculated results, if any, given by the numerical methods tabulated in Table 1 could be attributed to the effects of the method itself, and not to the turbulence models invoked for the stress-gradient terms in the governing turbulent flow equations. Validation of current calculations are achieved by comparing them with the channel flow expirements of Laufer (1951) [14] and the direct numerical simulation (DNS) data of Kim et al. (1987) [15].

3. Numerical Methods and Their Characteristics

Various numerical methods, which differ in the way they treat the stress-gradient terms in the momentum equations, are available. Merits and/or defects of each method had been thoroughly discussed in So (1996) [16] with the calculations of swirling and non-swirling flow modelling in mind. In the discussion, four methods (refer to Table 1) were identified as more commonly used and they had been adoped by Aksoy et al. (1995) [7] and Aksoy and So (1995) [8] in their investigations of square duct flow and backstep flow, respectively. These methods are labeled as: Method 1: explicit source method, Method 2: isotropic viscosity method, Method 3: anisotropic viscosity method, and Method 4: semi-coupled method. These methods and their characterists are briefly descriped below, while their mathematical details are provided in a later section.

Method 1 It evaluates the Reynolds-stress gradients in a straightforward manner and treats them as explicit source terms. It pays no attention to the absence of turbulent diffusion and the decoupling of the Reynolds stresses and the mean velocity variables. This scheme is expected to be stable for a wide range of simple flows, where streamline curvature and/or other external body force effects are absent or relatively small.

Method 2 It adopts Equation (1) directly and is similar to the calculations carried out by Amano and Goel (1985) [17] and Demuren (1992) [18]. It is commonly known as the isotropic viscosity method, or Method 2 as designated in Aksoy et al. (1995) [6]. The Reynolds stresses in the momentum equations are replaced by Equation (1) with ${\nu }_t$ defined by k and its dissipation rate, ε. As a result, the explicit coupling between the Reynolds stresses and the momentum equations is removed. The Reynolds stresses are still obtained from their respective transport equations and the normal stresses are used to evaluate k and ε. Therefore, this method not only models $-\overline{u_i{u}_j}$, but also the coupling between $-\overline{u_i{u}_j}$ and the mean momentum equations.

Method 3 It assumes anisotropic viscosity and replaces ${\nu }_t$ in Equation (1) by ${\nu }_t^n$, whose components: ${\nu }_t^1$, ${\nu }_t^2$ ${\nu }_t^3$, ${\nu }_t^4$, ${\nu }_t^5$, and ${\nu }_t^6$ are obtained from the 11, 22, 33, 12, 13, and 23 components of the Reynolds stress tensor, respectively. The method brings in stability to the solution procedure by introducing implicit diffusion to the momentum equations and yet without completely sacrificing the anisotropic behavior of the Reynolds stresses with the associated mean strain rates; thus, strongly couples the Reynolds stresses to the velocity variables. However, the calculated anisotropic viscosities could become negative or very large during calculations; thus, the method requires the use of limiters, which are usually applied as multiples of ${\nu }_t$ (So et al. (1988) [19]; Lai et al. (1991) [20]). This could lead to numerical errors in the results.

Method 4 It makes use of the SIMPLE algorithm of Patankar (1980) [2], Rhie and Chow (1983) [21], and Peric et al. (1988) [22] directly. Also, it suggests a linear interpolation technique on a collocated grid arrangement where the velocity variable is coupled to the pressure variable in such a way that the occurrence of an oscillatory pressure or velocity field is avoided. Therefore, one common control volume can be used for all the variables. This special interpolation technique has been successfully applied to the solution of vastly different laminar and turbulent flows using finite-volume method by Rhie and Chow (1983) [21], Peric et al. (1988) [22], Miller and Schmidt (1988) [23], Majumdar (1988) [24], and Lien and Leschziner (1993) [25]. Along similar line, Obi et al. (1991) [26] proposed a numerical treatment for the stress-gradient terms by extracting apparent viscosities from the discretized momentum equations and suggested to couple the velocity and the Reynolds-stress variables through the use of a special interpolation technique analogous to that introduced by Peric et al. (1988) [22] for the velocity and pressure variables. In this way, they were able to calculate all flow variables at collocated grid locations; thus making the method especially favorable when three-dimensional flows are studied in a generalized curvilinear coordinate system. This treatment is termed the semi-coupled method since only a portion of the Reynolds-stress equations are coupled to the momentum equations.

4. Near-Wall Modelling

The effects of near-wall turbulence models on flow in a square duct have been thoroughly investigated by Aksoy and So (1994) [27]. Even then, improvements on the overall replication of the corner flow, including the formation of corner vortices, are not as satisfactory. The effects of numerical methods on the calculations of complex turbulent flows in a square duct has also been carried out by Aksoy et al. (1995) [6], and the results showed that calculations of the corner vortices are influenced more by the numerical methods used to evaluate the stress-gradient terms. In their study, they found that Method 3 and Method 4, the semi-coupled method, erroneously predict an absence of corner vortices, thus giving rise to a zero vorticity field everywhere. Numerical methods were also found to have an effect on the calculation of a backstep flow by Aksoy and So (1995) [7], who used Methods 2, 3 and 4 to evaluate the stress-gradient terms in the Reynolds-stress equations. According to their study, iterations with Method 3 divereged rapidly, therefore, it was not possible to obtain a converged solution with Method 3. Furthermore, Methods 2 and 3 could give rise to calculation instability. It might be possible to obtain converged solutions if different numerical algorithms or turbulence models were used. The turbulence field in the corner of a square duct and behind a backstep are quite complex; therefore, it is not clear whether the errors observed in these calculations are the result of inapppropriate turbulence models adopted or are given rise by the numerical methods used in modelling the stress-gradient terms in the Reynolds-stress equations.

5. The Near-Wall Reynolds-Stress Model

Incompressible turbulent flows are considered. The flows are described by the Reynolds-averaged equations, which can be written in Cartesian tensor form as

$\frac{\partial U_i}{\partial x_i}=0$ (2)

$\frac{D{U}_i}{Dt}=-\frac{1}{\rho }\frac{\partial P}{\partial x_i}+\nu \frac{{\partial }^2{U}_j}{\partial x_i\partial x_j}-\frac{\partial \overline{u_i{u}_j}}{\partial x_j}$. (3)

Here, the Einstein summation convention is followed, U_i is the velocity in the ith direction, x_i is the coordinate system with x₁ denoting the stream direction, x₂ denotes the normal to x₁ in the flow plane and x₃ is normal to the flow plane, P is the mean pressure, ρ is the fluid density and ν is the kinematic viscosity. If a Reynolds-stress model is invoked, the equation governing the transport of the Reynolds stresses can be symbolically written as

$C_{ij}=D_{ij}^{\nu }+D_{ij}^t+P_{ij}+{\Pi }_{ij}-{\epsilon }_{ij}$, (4)

where the terms from left to right represent the convection, molecular diffusion, turbulent diffusion, production, velocity-pressure-gradient correlation and viscous dissipation of the Reynolds stress tensor, $-\overline{u_i{u}_j}$. The pressure-strain model of Speziale et al. (1991) [28] is adopted as the high-Reynolds-number model for ${\Pi }_{ij}$ and a near-wall correction, ${\Pi }_{ij}^w$, is formulated so that the pressure-diffusion near a wall is also accounted for. Thus modified, the model for ${\Pi }_{ij}$ can be written as

$\begin{array}{c}{\Pi }_{ij}=-\left(C_1\epsilon +C_1^{*}\tilde{P}\right)b_{ij}+C_2\epsilon \left(b_{ik}{b}_{kj}-\Pi {\delta }_{ij}/3\right)-\left(P_{ij}-2\tilde{P}{\delta }_{ij}/3\right)\\ +{\beta }_1\left(D_{ij}-2\tilde{P}{\delta }_{ij}/3\right)-2\left({\gamma }_1+C_3^{*}{\Pi }^{1/2}/2\right)k{S}_{ij}+{\Pi }_{ij}^w\end{array}$ (5)

where $b_{ij}=\left(\overline{u_i{u}_j}-2k{\delta }_{ij}/3\right)/2k$ is the anisotropic tensor, $\tilde{P}=P_{ii}/2$ is the production of k, $\Pi =b_{ij}{b}_{ij}$, and the near-wall correction is derived to be

$\begin{array}{c}{\Pi }_{ij}^w=f_{w1}[\left(C_1\epsilon +C_1^{*}\tilde{P}\right)b_{ij}-C_2\epsilon \left(b_{ik}{b}_{kj}-\Pi {\delta }_{ij}/3\right)\\ +{\epsilon }^{*}\left(P_{ij}-2\tilde{P}{\delta }_{ij}/3\right)+2{\gamma }^{*}k{S}_{ij}]+{\Pi }_{ij}^P\end{array}$ (6)

with

$\begin{array}{c}{\Pi }_{ij}^P=\frac{1}{3}\left[\frac{\partial }{\partial x_l}\left(\nu \frac{\partial \overline{u_i{u}_k}}{\partial x_l}\right)n_k{n}_j+\frac{\partial }{\partial x_l}\left(\nu \frac{\partial \overline{u_j{u}_k}}{\partial x_l}\right)n_k{n}_i\right]\\ +\frac{1}{3}\frac{\partial }{\partial x_m}\left(\nu \frac{\partial \overline{u_l{u}_k}}{\partial x_m}\right)n_k{n}_l{n}_i{n}_j\end{array}$ (7)

The model for the dissipation rate tensor is given by

${\epsilon }_{ij}=\frac{2}{3}\epsilon {\delta }_{ij}+f_{w1}\frac{\epsilon }{k}\left[-\frac{2}{3}k{\delta }_{ij}+\frac{\overline{u_i{u}_j}+\overline{u_i{u}_k}{n}_k{n}_j+\overline{u_j{u}_k}{n}_k{n}_i+n_i{n}_j\overline{u_i{u}_k}{n}_l{n}_k}{1+3\overline{u_k{u}_l}{n}_l{n}_k/2k}\right]$. (8)

In these equations, n_i = (0, 1, 0) is the wall unit normal and f_w1 is a damping function defined as $f_{w1}=\exp \left[-{\left({\mathrm{Re}}_t/200\right)}^2\right]$, where ${\mathrm{Re}}_t=k^2/\nu \epsilon$ is the turbulent Reynolds number. The high-Reynolds-number model of Hanjalic and Launder (1972) [29] is adopted for $D_{ij}^t$, or

$D_{ij}^t=\frac{\partial }{\partial x_k}\left[C_s\frac{k}{\epsilon }\left(\overline{u_i{u}_l}\frac{\partial \overline{u_j{u}_k}}{\partial x_l}+\overline{u_j{u}_l}\frac{\partial \overline{u_k{u}_i}}{\partial x_l}+\overline{u_k{u}_l}\frac{\partial \overline{u_i{u}_j}}{\partial x_l}\right)\right]$. (9)

Since the exact and modelled terms are of order y³ and higher near a wall, no near-wall correction for $D_{ij}^t$ is necessary. The ε transport equation is also formulated to allow integration to the wall by appropriately modifying the high-Reynolds-number modeled equation. Its final form is given by

$\frac{D\epsilon }{Dt}=\frac{\partial }{\partial x_j}\left(\nu \frac{\partial \epsilon }{\partial x_j}\right)+\frac{\partial }{\partial x_j}\left(C_{\epsilon }\frac{k}{\epsilon }\overline{u_i{u}_j}\frac{\partial \epsilon }{\partial x_i}\right)+C_{\epsilon 1}\frac{\epsilon }{k}\tilde{P}-C_{\epsilon 2}\frac{\epsilon \tilde{\epsilon }}{k}+\xi$, (10)

where

$\xi =f_{w2}\left(-N\frac{\epsilon \tilde{\epsilon }}{k}+M\frac{{\bar{\epsilon }}^2}{k}-L\frac{\epsilon }{k}\tilde{P}\right)$, (11)

is the near-wall correcting function and decays to zero away from the wall, $\tilde{\epsilon }=\epsilon -2\nu {\left(\partial \sqrt{k}/\partial x_i\right)}^2$ and $\bar{\epsilon }=\epsilon -2\nu k/x_2^2$ are reduced ε, and x₂ is the normal coordinate. The damping function in (11) is given by $f_{w2}=\exp \left[-{\left({\mathrm{Re}}_t/40\right)}^2\right]$ and the model constants are C₁ = 3.4, C₂ = 4.2, $C_1^{\ast }$ = 1.8, $C_3^{\ast }$ = 1.3, α₁ = 0.4125, β₁ = 0.2125, γ₁ = 0.01667, α* = −0.29, γ* = 0.065, C_s = 0.11, C_ε = 0.12, C_ε1 = 1.5, C_ε2 = 1.83, L = 2.25, M = 0.5 and N = 0.57. Further details and the rationale underlying the assumed constants for this near-wall Reynolds-stress model are provided in So et al. (1996) [30].

6. Numerical Treatment of the Stress-Gradient Terms

In the present investigation, four different methods of treating the stress gradient terms in the mean flow equations are examined. These are: Method 1: explicit source method, Method 2: isotropic viscosity method, Method 3: anisotropic viscosity method and Method 4: semi-coupled method. The numerical mathematics for each are briefly summarized below. For the sake of conciseness, only the source term of the x-momentum equation is presented.

Method 1: Explicit Source Method

With this approach, the stress gradient term in Equation (3) is considered to be the source term. Thus, integrating the x-momentum equation around the control volume shown in Figure 1 gives

Figure 1. Typical control volume and grid system.

${\bar{S}}_u=-\rho \left(\frac{\overline{u_E^2}-\overline{u_W^2}}{2}\right)\Delta y\Delta z-\rho \left(\frac{{\overline{uv}}_N-{\overline{uv}}_S}{2}\right)\Delta x\Delta z-\rho \left(\frac{{\overline{uw}}_T-{\overline{uw}}_B}{2}\right)\Delta x\Delta y$, (12)

where the Reynolds stresses at the control-volume faces are obtained by linear interpolation from the neighboring grids, such as $\overline{u_e^2}=\left(\overline{u_E^2}+\overline{u_P^2}\right)/2$ if the cell face is located midway between the grids. It can be seen from this equation that the x-momentum equation at grid P is affected by the Reynolds-stresses at the neighboring grids but not at grid P. Therefore, similar to the problem which arises as a result of pressure-velocity decoupling (Patankar (1980) [2]), an oscillatory variation of the Reynolds stresses could result during each iteration of the mean momentum equation. This could lead to a physically unrealistic solution or the iterative process could diverge.

Method 2: Isotropic Viscosity Method

With this method the momentum equation is cast into a form given by,

$\frac{D{U}_i}{Dt}=-\frac{1}{\rho }\frac{\partial P}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\left(\nu +{\nu }_t\right)\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)\right]-\frac{\partial {\overline{u_i{u}_j}}^p}{\partial x_j}$, (13)

where ${\overline{u_i{u}_j}}^p=\overline{u_i{u}_j}+2{\nu }_t{S}_{ij}$ is the pseudo stress tensor and ν_t is evaluated from k and ε. In essence, a turbulent diffusion term is added to the mean momentum equation. This term is treated implicitly together with the viscous diffusion term. However, the same term is subtracted from the source term, which is treated explicitly. When the iterations converge, the two terms cancel out, and a correct solution is obtained. In this respect, the method is a deferred-corrector method.

Method 3: Anisotropic Viscosity Method

With this approach, Equation (1) is still invoked. However, unlike Equation (1), each component of the stress is assumed to have a different eddy viscosity. In this way, six different eddy viscosities are defined using an expression of the form,

$-\overline{u_i{u}_j}=2{\nu }_t^n{S}_{ij}-\frac{2}{3}k{\delta }_{ij}$, $n=1,2,\cdots ,6$. (14)

Substituting this expression into the x-component of Equation (3) and making use of the continuity equation, the x-momentum equation can be cast into the following form,

$\begin{array}{c}\frac{DU}{Dt}=-\frac{1}{\rho }\frac{\partial P}{\partial x}+\frac{\partial }{\partial x}\left[\left(\nu +{\nu }_t^1\right)\frac{\partial U}{\partial x}\right]+\frac{\partial }{\partial y}\left[\left(\nu +{\nu }_t^4\right)\frac{\partial U}{\partial y}\right]\\ +\frac{\partial }{\partial z}\left[\left(\nu +{\nu }_t^6\right)\frac{\partial U}{\partial z}\right]+S_u\end{array}$ (15)

where the source term S_u is given by

$S_u=\frac{\partial }{\partial x}\left[\left(\nu +{\nu }_t^1\right)\frac{\partial U}{\partial x}\right]+\frac{\partial }{\partial y}\left[\left(\nu +{\nu }_t^4\right)\frac{\partial U}{\partial x}\right]+\frac{\partial }{\partial z}\left[\left(\nu +{\nu }_t^6\right)\frac{\partial U}{\partial x}\right]-\frac{2}{3}\frac{\partial k}{\partial x}$, (16)

and the anisotropic viscosities ${\nu }_t^n$ are calculated from Equation (14). For most flows, the strain-rate terms could vanish or assume an arbitrary sign in different regions of the flow field. Consequently, ${\nu }_t^n$ obtained from Equation (14) could become negative or very large during the calculations. A negative diffusion or an extremely large ${\nu }_t^n$ could become detrimental in the iterative procedure. Therefore, in order to prevent the iterative scheme from diverging, this method requires the use of limits as upper and lower bounds on the calculated viscosities. This in turn introduces numerical errors to the results. The upper and lower limits are usually applied as multiples of ν_t, such as $L_l{\nu }_t\le {\nu }_t^n\le L_u{\nu }_t$, where L_l and L_u are usually chosen to be of the order of 0.1 and 10, respectively.

Method 4: Semi-Coupled Method

The approach of Obi et al. (1991) [26] is based on the reasoning that the turbulent Reynolds stresses are coupled to the associated strain rates by an expression similar to Equation (1), and that a strong link between the Reynolds stresses and these strain rates needs to be formed in order to avoid a diverging iterative solution. Considering the x-component of Equation (3) and treating the stress gradient term as the source, S_u, the integration of this equation around the control volume shown in Figure 2 gives rise to the following source term,

${\bar{S}}_u=-\rho \left(\overline{u_e^2}-\overline{u_w^2}\right)\Delta y\Delta z-\rho \left(\overline{u{v}_n}-\overline{u{v}_s}\right)\Delta x\Delta z-\rho \left(\overline{u{w}_t}-\overline{u{w}_b}\right)\Delta x\Delta y$. (17)

If the cell face Reynolds stresses in Equation (17) are obtained using a linear interpolation of the values at the neighboring grid points, then Equation (17) reduces to Equation (12) and the explicit-source method is recovered. Instead, the present method calculates the associated strain rates from the discretized Reynolds-stress equations which can be written in the following form,

$\overline{u_p^2}={\left(\frac{\Sigma A_{nb}\overline{u_{nb}^2}+\bar{S}}{A_p}\right)}_P-{\left({\Gamma }_U^{11}\frac{U_e-U_w}{\Delta x}\right)}_P$, (18a)

$\overline{u{v}_p}={\left(\frac{\Sigma A_{nb}\overline{u{v}_{nb}}+\bar{S}}{A_p}\right)}_P-{\left({\Gamma }_U^{12}\frac{U_n-U_s}{\Delta y}\right)}_P$, (18b)

$\overline{u{w}_p}={\left(\frac{\Sigma A_{nb}\overline{u{w}_{nb}}+\bar{S}}{A_p}\right)}_P-{\left({\Gamma }_U^{13}\frac{U_t-U_b}{\Delta z}\right)}_P$, (18c)

Figure 2. Staggered control volume for U in the x-y plane.

where the A’s are coefficients, Γ’s are diffusivities and $\bar{S}$ is the linearized source term. The symbols and nomenclature adopted in this section are identical to those used by Patankar (1980) [2]. Similar expressions can now be written for the discretized Reynolds stresses at the cell faces,

$\overline{u_e^2}={\langle \frac{\Sigma A_{nb}\overline{u_{nb}^2}+\bar{S}}{A_p}\rangle }_e-{\langle {\Gamma }_U^{11}\rangle }_e\frac{U_E-U_P}{\delta x_e}$, (19a)

$\overline{u{v}_n}={\langle \frac{\Sigma A_{nb}\overline{u{v}_{nb}}+\bar{S}}{A_p}\rangle }_n-{\langle {\Gamma }_U^{12}\rangle }_n\frac{U_N-U_P}{\delta y_n}$, (19b)

$\overline{u{w}_t}={\langle \frac{\Sigma A_{nb}\overline{u{w}_{nb}}+\bar{S}}{A_p}\rangle }_t-{\langle {\Gamma }_U^{13}\rangle }_t\frac{U_T-U_P}{\delta z_t}$, (19c)

where $\langle \rangle$ denotes linear interpolation. When the Reynolds stresses are evaluated from the expressions given by Equation (19), they are directly coupled with the mean velocity at the neighboring grids in the same way if a staggered grid system were used. If the cell faces lie midway between the grid points, the $\langle \rangle$ terms are reduced to

${\langle \frac{\Sigma A_{nb}\overline{u_{nb}^2}+\bar{S}}{A_p}\rangle }_e=\frac{1}{2}\left[{\left(\frac{\Sigma A_{nb}\overline{u_{nb}^2}+\bar{S}}{A_p}\right)}_E+{\left(\frac{\Sigma A_{nb}\overline{u_{nb}^2}+\bar{S}}{A_p}\right)}_P\right]$, (20a)

${\langle {\Gamma }_U^{11}\rangle }_e=\frac{1}{2}\left[{\langle {\Gamma }_U^{11}\rangle }_E+{\langle {\Gamma }_U^{11}\rangle }_P\right]$, (20b)

etc. Substituting the expressions given by Equation (19) into the right hand side of Equation (17),

$\begin{array}{c}\overline{S_U}=\frac{{\langle {\Gamma }_U^{11}\rangle }_e\Delta y\Delta z}{\delta x_e}\left(U_E-U_P\right)-\frac{{\langle {\Gamma }_U^{11}\rangle }_w\Delta y\Delta z}{\delta x_w}\left(U_E-U_P\right)+S_U^{11}\\ +\frac{{\langle {\Gamma }_U^{12}\rangle }_n\Delta x\Delta z}{\delta y_n}\left(U_N-U_P\right)-\frac{{\langle {\Gamma }_U^{12}\rangle }_w\Delta x\Delta z}{\delta y_s}\left(U_P-U_S\right)+S_U^{12}\\ +\frac{{\langle {\Gamma }_U^{13}\rangle }_t\Delta x\Delta y}{\delta z_t}\left(U_T-U_P\right)-\frac{{\langle {\Gamma }_U^{13}\rangle }_s\Delta y\Delta z}{\delta z_b}\left(U_P-U_B\right)+S_U^{13}\end{array}$ (21)

where

$S_U^{11}=\left[{\langle \frac{\Sigma A_{nb}\overline{u_{nb}^2}+\bar{S}}{A_p}\rangle }_w-{\langle \frac{\Sigma A_{nb}\overline{u_{nb}^2}+\bar{S}}{A_p}\rangle }_e\right]\Delta y\Delta z$, (22a)

$S_U^{12}=\left[{\langle \frac{\Sigma A_{nb}\overline{u{v}_{nb}}+\bar{S}}{A_p}\rangle }_s-{\langle \frac{\Sigma A_{nb}\overline{u{v}_{nb}}+\bar{S}}{A_p}\rangle }_n\right]\Delta x\Delta z$, (22b)

$S_U^{13}=\left[{\langle \frac{\Sigma A_{nb}\overline{u{w}_{nb}}+\bar{S}}{A_p}\rangle }_b-{\langle \frac{\Sigma A_{nb}\overline{u{w}_{nb}}+\bar{S}}{A_p}\rangle }_t\right]\Delta x\Delta y$, (22c)

During the computation of the mean momentum equations, the diffusion-like terms in Equation (21) are treated implicitly and they directly contribute to the coefficients $A_E$, $A_W$, $A_N$, $A_S$, $A_T$, $A_B$, and $A_P$. On the other hand, $S_U^{11}$, $S_U^{12}$, and $S_U^{13}$ are calculated explicitly as source terms.

The apparent-viscosities extracted from the discretized Reynolds-stress equations depend on the second-order closure adopted. If the present near-wall Reynolds-stress model is invoked at grid P, the coefficients for the x-momentum equation can be determined and they are given by the following expressions:

$\begin{array}{c}{\Gamma }_U^{11}=\{2\left[1-\frac{2}{3}\left({\alpha }_1-{\alpha }^{*}{f}_{w1}+{\beta }_1\right)\right]\rho \overline{u_P^2}+2\left({\gamma }_1-f_{w1}{\gamma }^{*}\right)\rho k_P\begin{array}{c}\\ \\ \end{array}\\ +\frac{1}{3}\rho \left(1-f_{w1}\right)C_1^{*}\overline{u_P^2}+\rho C_3^{*}{\Pi }_P^{1/2}{k}_P-\frac{1}{2}\rho \left(1-f_{w1}\right)C_1^{*}\frac{{\left(\overline{u_P^2}\right)}^2}{k_P}\}\frac{\Delta x\Delta y\Delta z}{A_P^{\overline{u^2}}}\end{array}$ (23a)

$\begin{array}{c}{\Gamma }_U^{12}=\{-{\beta }_1\rho \overline{u_P^2}+\left(1-{\alpha }_1+{\alpha }^{*}{f}_{w1}\right)\rho \overline{v_P^2}+\left({\gamma }_1-f_{w1}{\gamma }^{*}\right)\rho k_P\begin{array}{c}\\ \\ \end{array}\\ +\frac{1}{2}\rho C_3^{*}{\Pi }_P^{1/2}{k}_P-\rho \left(1-f_{w1}\right)C_1^{*}\frac{{\left(\overline{uv}\right)}^2}{k_P}\}\frac{\Delta x\Delta y\Delta z}{A_P^{\overline{uv}}}\end{array}$ (23b)

$\begin{array}{c}{\Gamma }_U^{13}=\{-{\beta }_1\rho \overline{u_P^2}+\left(1-{\alpha }_1+{\alpha }^{*}{f}_{w1}\right)\rho \overline{w_P^2}+\left({\gamma }_1-f_{w1}{\gamma }^{*}\right)\rho k_P\begin{array}{c}\\ \\ \end{array}\\ +\frac{1}{2}\rho C_3^{*}{\Pi }_P^{1/2}{k}_P-\rho \left(1-f_{w1}\right)C_1^{*}\frac{{\left(\overline{uw}\right)}^2}{2k_P}\}\frac{\Delta x\Delta y\Delta z}{A_P^{\overline{uw}}}\end{array}$ (23c)

For numerical stability, it is important that extracted viscosities given by Equation (23) are positive. It can be seen from these equations that both near the wall ( $f_{w1}=1$ ) and away from the wall ( $f_{w1}=0$ ), they assume positive values. If a convergent solution field cannot be achieved with the set of coefficients given above, the terms which act to decrease the value of the apparent diffusion coefficients may be removed from these expressions. However, in the present study, it was possible to obtain physically realistic, convergent solutions with the set of coefficients given above for all cases examined. The converged results are found to be independent of the expressions chosen for the apparent viscosities.

When Obi et al. (1991) [26] introduced this new method, they applied it to the numerical solution of the turbulent flow over a backward-facing step using the high-Reynolds-number second-order model of Gibson and Launder (1978) [31]. However, their computational results indicated the presence of unrealistic velocity fields in the neighborhood of the reattachment point close to the wall. This point will be discussed further below.

For ease of reference later, the different methods discussed above are designated as in Table 1. In this table, the SIMPLEC method of van Doormal and Raithby (1984) [3] is an improved version of SIMPLE. It is used to improve the convergence rate of the SIMPLE algorithm when three-dimensional computations are carried out.

7. Numerical Soultion

All computations of square duct flow by Aksoy et al. (1995) [6] and backstep flow by Aksoy and So (1995) [7] were carried out using the SIMPLE/SIMPLEC algorithm, and the convective/diffusive terms are discretized with the power-law scheme described in Patankar (1980) [2]. A three-dimensional elliptic code has been used to calculate the fully-developed turbulent flow inside the square duct while a two-dimanional code was used in the case of backstep flow. Periodic flow conditions are applied in the axial direction. Calculations are carried out using 4 × 91 × 91 grids in the axial (x), the normal (y) and the transverse (z) directions, respectively. With this grid resolution, the results were found to be grid independent for both Reynolds numbers attempted. The grids were placed such that at least 5 grids are located within y⁺ < 5 with the first grid point at y⁺ ≃ 1, and 25 grids are located in the region 5 < y⁺ < 65. The developing flow at Re = 250,000 is solved using an iterative, forward marching solution procedure based on the parabolized equations as described by Patankar and Spalding (1972) [4]. The same grid spacing (91 × 91) is adopted in the cross sectional plane. The forward step size is varied progressively from 0.4 percent of the hydraulic diameter near the inlet to 4 percent until the flow becomes fully developed. The iterations are carried out until the normalized residual sources reduce to 10⁻⁵.

A two-dimensional (2-D) elliptic code is used to calculate the backstep flow. The governing equations are discretized without the unsteady terms, put into a tridiagonal form and solved by using the SIMPLE/SIMPLEC algorithm until the normalized residual sources for the momentum equations and mass reduce to 10⁻⁵. Special attention is given to the placement of the grids near the wall. Again, a minimum of 5 grids are placed within y⁺ < 5 with the first grid point at y⁺ ≃ 1, and about 30 grids are located in the region 5 < y⁺ < 75. Different grid sizes of 51 × 41, 91 × 81, 131 × 151 in the x and y directions are tested for grid independence and a 91 × 81 grid is found to be sufficient. This 2-D elliptic algorithm was found to be quite effective in its calculation of backstep flow. Therefore, it is chosen as the algorithm of choice for the calculation of a developing plane channel flow in the present study.

8. Results and Discussion

Methods 1 - 4 had also been used by Aksoy et al. (1995) [6] to study a square duct flow and by Aksoy and So (1995) [7] to simulate a backstep flow. Therefore, their findings will be drawn on in the following discussion to bring out the indiviual effects of these four numerical methods on the calculation of complex flows. These calculations and the present developing plane channel flow calculation are all carried out using the Reynolds-stress model discussed above. It is anticipated that through this analysis, the merits or lack thereof of each method can be demonstrated.

The results of the developing plane channel flow calculations are discussed below together with the performance of the four numerical methods. Calculations are carried out at Re = 7900 and 30,800, based on the maximum velocity in the fully developed region and the channel half width. These Re’s are the same as the experimental measurements of Laufer (1951) [14] and the DNS data of Kim et al. (1987) [15]. Figure 3 shows the development of the skin-friction coefficient C_f along the channel wall. Here, C_f is defined with respect to the bulk velocity. The DNS data and the measurements are also shown for comparison. Predictions from Methods 1, 2 and 3 give almost identical results and agree well with both sets of data. However, Method 4 underpredicts C_f at both Re by about 20% in the fully-developed state. This is a significant under prediction; thus indicating that in comparison with Methods 1, 2 and 3, Method 4 is not suitable for the calculation of a simple developing channel flow.

In Figure 4, predicted Reynolds normal stresses are compared with DNS data. Again, results from Method 4 differ from those deduced by using the other three methods. However, all four methods yield the same result for the shear stress as shown in Figure 5. The discrepancies noted in Figure 3 and Figure 4 for Method 4 are believed to be the result of numerical errors introduced by the interpolations that need to be implemented as part of this method. In other words, the interpolations have an undue influence on the calculated results. The calculations with Method 1 can be considered more accurate since no limitations, interpolations, or deferred-corrector type terms are introduced. The results with Method 3 are obtained by forcing all components of the calculated anisotropic viscosities to remain within a lower bound of 0.1ν_t and a higher bound of 10ν_t. This means that use of lower and upper bounds in Method 3 is required in order to prevent the iterations from diverging due to the existence of negative or excessively large anisotropic viscosities during the computations.

Figure 3. Comparison of the calculated C_f for plane channel flow.

Figure 4. Comparison of the calculated normal stresses for plane channel flow.

Figure 5. Comparison of the calculated shear stress for plane channel flow.

In general, the four methods give essentially the same results for all variables examined except for C_f. Therefore, in one way or another, all four numerical methods used to treat the stress gradient terms have some defects in the calculation of a 2-D developing channel flow. Overall, it could be concluded that the best performing numerical method is given by Method 4, in spite of its underprediction of C_f.

Finally, the relative effectiveness of these four numerical methods on the calculations of complex flows such as a square duct flow (Aksoy et al. (1995) [6]) and a backstep flow (Aksoy and So (1995) [7]) are examined here for comparison. Since the secondary flows in a square duct are generated primarily by the gradients of the Reynolds stresses (Aksoy et al. (1995) [6]), they are very sensitive to the type of numerical methods and turbulence models adopted. The following discussion is based on the fact that the same turbulence models and numerical methods are invoked in this study and the other two on complex flows. Methods 1 and 4 are found to give identical and physically correct results. The interpolations incorporated as part of Method 4 does not seem to introduce any significant errors to the calculations. Although it is known to be unstable for most flows, Method 1 can still be used to study the three-dimensional turbulent flow inside a square duct. Since Method 4 introduces implicit diffusion terms to the momentum equations by extracting them from the discretized Reynolds-stress equations and resolves the velocity-Reynolds-stress decoupling through a special interpolation technique, it can be used to solve different types of complex flows. The limitations applied to the anisotropic viscosities with Method 3 may actually result in the introduction of an isotropic viscosity into the momentum equations for most of the flow domain, thus leading to incorrect predictions. Furthermore, this method gives rise to zero secondary flow in the case of a square duct. Method 2 also predicts a zero secondary flow field. Therefore, Methods 2 and 3 are not suitable in spite of the fact that the mean velocities and Reynolds stresses are calculated correctly. The results presented for the square duct flow show that, even if the Reynolds stress aniostropy is captured by the turbulence model, unless the stress gradient terms in Equation (3) are treated correctly, there is no guaratee that a secondary flow will form in the calculation. On the other hand, the improved calculations of a backstep flow using Method 3 compared to using Method 4 could be explained as follows. For a backstep flow, the secondary flow is driven by the pressure gradient term in Equation (3). Therefore, even if the stress-gradient term in Equation (3) is not treated properly by Method 3, the net effect on the calculation of the recirculation region is insignificant.

In the case of a backstep flow (Aksoy and So (1995) [7]), the following conclusions can be drawn. Converged solutions could be obtained only with Methods 3 and 4. Even then, the C_f distribution downstream of the step is found to be very sensitive to the numerical method adopted. This last point is particularly true for Method 3 because anisotropic viscosities bounds need to be enforced as part of Method 3. The behavior of the calculated velocities using Methods 3 and 4 are compared near the reattachment point downstream of the backstep. It was found that only the velocity obtained with Method 4 shows a forward flow near the wall in the recirculation region; thus, indicating the presence of kinked streamlines. This is contrary to the findings of Lasher and Taulbee (1992) [32], who attributed the unrealistic behavior to the numerical method adopted and not to the pressure-strain model. In conclusion, Aksoy and So (1995) [7] found that only Method 3 can give rise to physically correct results that agree better with DNS data than those obtained using Method 4.

9. Conclusion

When second-order models are used to calculate turbulent flows, the type of numerical methods used to treat the stress-gradient terms in the mean momentum equations becomes important from both the stability and accuracy point of view. Methods 1 - 4 have been tested using the near-wall Reynolds-stress models of So et al. (1996) [30] in the calculation of the turbulent characteristics in a developing plane channel flow. For this simple flow, all four methods are stable. With the exception of Method 4, they all give approximately identical results except for the development of C_f. with Method 4, the calculated C_f along the channel wall is in error by ~20%. This error could be a consequence of the interpolations implemented in this method. Therefore, this investigation and those carried out by Aksoy et al. (1995) [6] and Aksoy and So (1995) [7] for complex flows show that Method 4 might work well for complex flows. However, it is not appropriate for a simple developing plane channel flow calculation because the error in C_f is too large to ignore (Figure 3). In other words, any method that works for complex flows is not necessarily suitable for simple flows.

Acknowledgements

Administrative and library support is granted to Ronald M.C. So by the Department of Mechanical Engineering, The Hong Kong Polytechnic University, during the course of this study, which is gratefully acknowledged. With their help, he was able to get hold of some key references cited in this paper.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Harlow, F.H. and Welch, J.E. (1965) Numerical Calculation of Time Dependent Viscous Incompressible Flow of Fluid with Free Surface. Physics of Fluids, 8, 2182-2189. https://doi.org/10.1063/1.1761178
[2]	Patankar, S.V. (1980) Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York.
[3]	van Doormal, J.P. and Raithby, G.D. (1984) Enhancements of the SIMPLE Method for Predicting Incompressible Fluid Flows. Numerical Heat Transfer, 7, 147-163. https://doi.org/10.1080/01495728408961817
[4]	Patankar, S.V. and Spalding, D.B. (1972) A Calculation Procedure for Heat, Mass and Momentum Transfer in 3D Parabolic Flows. International Journal of Heat and Mass Transfer, 15, 1787-1806. https://doi.org/10.1016/0017-9310(72)90054-3
[5]	Pope, S.B. and Whitelaw, J.H. (1976) The Calculation of Near Wake Flows. Journal of Fluid Mechanics, 73, 9-32. https://doi.org/10.1017/S0022112076001213
[6]	Aksoy, H., So, R.M.C. and Kotb, N.A. (1995) An Assessment of Different Numerical Techniques on the Calculation of Turbulent Square Duct Flows. Proceedings of the Tenth Symposium on Turbulent Shear Flows, Vol. 1, 7.13-7.18.
[7]	Aksoy, H. and So, R.M.C. (1995) Effects of Numerical Methods and Near-Wall Reynolds Stress Modelling on the Prediction of Backstep Flow. In: Otugen, M.V., et al., Eds., Separated and Complex Flows, ASME International, New York, Vol. 217, 265-272.
[8]	Launder, B.E. and Spalding, D.B. (1974) The Numerical Computation of Turbulent Flows. Computer Methods in Applied Mechanics and Engineering, 3, 269-289. https://doi.org/10.1016/0045-7825(74)90029-2
[9]	Drikakis, D. (2003) Advances in Turbulent Flow Computations Using High-Resolution Methods. Progress in Aerospace Sciences, 39, 405-424. https://doi.org/10.1016/S0376-0421(03)00075-7
[10]	Gunzburger, M. and Labovsky, A. (2012) High Accuracy Method for Turbulent Flow Problems. Mathematical Models and Methods in Applied Sciences, 22, Article ID: 1250005. https://doi.org/10.1142/S0218202512500054
[11]	Junqueire-Junior, C., Azevedo, J.L.F., Scalabrin, L.C. and Basso, E. (2013) A Study of Physical and Numerical Effects of Dissipation on Turbulent Flow Simulations. Journal of Aerospace Technology and Management, 5, 145-168. https://doi.org/10.5028/jatm.v5i2.179
[12]	Ahsan, M. (2014) Numerical Analysis of Friction Factor for a Fully Developed Turbulent Flow Using k-ε Turbulence Model with Enhanced Wall Treatment. Beni-Suef University Journal of Basic and Applied Sciences, 3, 269-277. https://doi.org/10.1016/j.bjbas.2014.12.001
[13]	Shahjada Tarafder, Md. and Al Mursaline, M. (2019) Numerical Analysis of Turbulent Flow around Two-Dimensional Bodies Using Non-Orthogonal Body-Fitted Mesh. International Journal of Applied Mechanics and Engineering, 24, 387-410. https://doi.org/10.2478/ijame-2019-0024
[14]	Laufer, J. (1951) Investigation of Turbulent Flow in a Two-Dimensional Channel. NACA Report No. 1053.
[15]	Kim, J., Moin, P. and Moser, R.D. (1987) Turbulence Statistics in Fully Developed Channel Flow at Low Reynolds Number. Journal of Fluid Mechanics, 177, 133-186. https://doi.org/10.1017/S0022112087000892
[16]	So, R.M.C. (1996) Development of Advanced Turbulence Models for Swirling and Non-Swirling Three-Dimensional Confined Flows. Final Report Submitted to Knolls Atomic Power Laboratory, Schenectady, Report No. ASU/TFML-96-10.
[17]	Amano, R.S. and Goel, P. (1985) Computations of Turbulent Flow beyond Backward-Facing Steps Using Reynolds-Stress Closure. AIAA Journal, 23, 1356-1361. https://doi.org/10.2514/3.9092
[18]	Demuren, A.O. (1992) Multigrid Acceleration and Turbulence Models for Computations of 3D Turbulent Jets in Cross Flow. International Journal of Heat and Mass Transfer, 35, 2783-2794. https://doi.org/10.1016/0017-9310(92)90299-8
[19]	So, R.M.C., Lai, Y.G., Hwang, B.C. and Yoo, G.J. (1988) Low-Reynolds-Number Modeling of Flows over a Backward-Facing Step. Journal of Applied Mathematics and Physics (ZAMP), 39, 13-27. https://doi.org/10.1007/BF00945719
[20]	Lai, Y.G., So, R.M.C. and Zhang, H.S. (1991) Turbulence Driven Secondary Flows in a Curved Pipe. Theoretical and Computational Fluid Dynamics, 3, 163-180. https://doi.org/10.1007/BF00271800
[21]	Rhie, C.M. and Chow, W.L. (1983) Numerical Study of the Turbulent Flow past an Airfoil with Trailing Edge Separation. AIAA Journal, 21, 1525-1532. https://doi.org/10.2514/3.8284
[22]	Peric, M., Kessler, R. and Scheuerer, G. (1988) Comparison of Finite-Volume Numerical Methods with Staggered and Collocated Grids. Computers and Fluids, 16, 389-403. https://doi.org/10.1016/0045-7930(88)90024-2
[23]	Miller, T.F. and Schmidt, F.W. (1988) Use of Pressure Weighted Interpolation Method for the Solution of the Incompressible Navier-Stokes Equations on a Non-Staggered Grid System. Numerical Heat Transfer, 14, 213-233. https://doi.org/10.1080/10407788808913641
[24]	Majumdar, S. (1988) Role of Underrelaxation in Momentum Interpolation for Calculation of Flow with Non-Staggered Grids. Numerical Heat Transfer, 13, 125-132. https://doi.org/10.1080/10407788808913607
[25]	Lien, F.S. and Leschziner, M.A. (1993) A General Non-Orthogonal Collocated Finite Volume Algorithm for Turbulent Flow at All Speeds Incorporating Second-Moment Turbulence-Transport Closure, Part 1: Computational Implementation. Computer Methods in Applied Mechanics and Engineering, 114, 123-148. https://doi.org/10.1016/0045-7825(94)90165-1
[26]	Obi, S., Peric, M. and Scheuerer, G. (1991) Second-Moment Calculation Procedure for Turbulent Flows with Collocated Variable Arrangement. AIAA Journal, 29, 585-590. https://doi.org/10.2514/3.10624
[27]	Aksoy, H. and So, R.M.C. (1994) Near-Wall Turbulence Models and Their Application to Flow in a Square Duct. In: Hussaini, M.Y., et al., Eds., Transition, Turbulence and Combustion, Kluwer Academic Publishers, Dordrecht, Vol. 2, 23-37. https://doi.org/10.1007/978-94-011-1034-1_3
[28]	Speziale, C.G., Sarkar, S. and Gatski, T.B. (1991) Modeling the Pressure-Strain Correlation of Turbulence: An Invariant Dynamical Systems Approach. Journal of Fluid Mechanics, 227, 245-272. https://doi.org/10.1017/S0022112091000101
[29]	Hanjalic, K. and Launder, B.E. (1972) A Reynolds Stress Model of Turbulence and its Application to Thin Shear Flows. Journal of Fluid Mechanics, 52, 609-638. https://doi.org/10.1017/S002211207200268X
[30]	So, R.M.C., Aksoy, H., Sommer, T.P. and Yuan, S.P. (1996) Modeling Reynolds-Number Effects in Wall-Bounded Turbulent Flows. Journal of Fluids Engineering, 118, 260-267. https://doi.org/10.1115/1.2817372
[31]	Gibson, M.M. and Launder, B.E. (1978) Ground Effects on Pressure Fluctuations in the Atmospheric Boundary Layer. Journal of Fluid Mechanics, 86, 491-511. https://doi.org/10.1017/S0022112078001251
[32]	Lasher, W.C. and Taulbee, D.B. (1992) On the Computation of Turbulent Backstep Flow. International Journal of Heat and Fluid Flow, 13, 30-40. https://doi.org/10.1016/0142-727X(92)90057-G