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
PDF
HTML XML
119
Downloads
488
Views
Citations

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

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

Share and Cite:

So, R. (2022) Effects of Numerical Methods on the Calculation of Developing Plane Channel Flow. *Journal of Applied Mathematics and Physics*, **10**, 2086-2104. doi: 10.4236/jamp.2022.106142.

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,
$\stackrel{\xaf}{uv}$,
$\stackrel{\xaf}{uw}$, and
$\stackrel{\xaf}{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

$-\stackrel{\xaf}{{u}_{i}{u}_{j}}=2{\nu}_{t}{S}_{ij}-\frac{2}{3}k{\delta}_{ij}$. (1)

In Equation (1),
$-\stackrel{\xaf}{{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}*,

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

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

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
$-\stackrel{\xaf}{{u}_{i}{u}_{j}}$, but also the coupling between
$-\stackrel{\xaf}{{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 \stackrel{\xaf}{{u}_{i}{u}_{j}}}{\partial {x}_{j}}$. (3)

Here, the Einstein summation convention is followed, *U _{i}* is the velocity in the

${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,
$-\stackrel{\xaf}{{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}^{*}\stackrel{\u02dc}{P}\right){b}_{ij}+{C}_{2}\epsilon \left({b}_{ik}{b}_{kj}-\Pi {\delta}_{ij}/3\right)-\left({P}_{ij}-2\stackrel{\u02dc}{P}{\delta}_{ij}/3\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}+{\beta}_{1}\left({D}_{ij}-2\stackrel{\u02dc}{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(\stackrel{\xaf}{{u}_{i}{u}_{j}}-2k{\delta}_{ij}/3\right)/2k$ is the anisotropic tensor,
$\stackrel{\u02dc}{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}^{*}\stackrel{\u02dc}{P}\right){b}_{ij}-{C}_{2}\epsilon \left({b}_{ik}{b}_{kj}-\Pi {\delta}_{ij}/3\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}+{\epsilon}^{*}\left({P}_{ij}-2\stackrel{\u02dc}{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 \stackrel{\xaf}{{u}_{i}{u}_{k}}}{\partial {x}_{l}}\right){n}_{k}{n}_{j}+\frac{\partial}{\partial {x}_{l}}\left(\nu \frac{\partial \stackrel{\xaf}{{u}_{j}{u}_{k}}}{\partial {x}_{l}}\right){n}_{k}{n}_{i}\right]\\ \text{\hspace{0.17em}}\text{\hspace{0.05em}}+\frac{1}{3}\frac{\partial}{\partial {x}_{m}}\left(\nu \frac{\partial \stackrel{\xaf}{{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{\stackrel{\xaf}{{u}_{i}{u}_{j}}+\stackrel{\xaf}{{u}_{i}{u}_{k}}{n}_{k}{n}_{j}+\stackrel{\xaf}{{u}_{j}{u}_{k}}{n}_{k}{n}_{i}+{n}_{i}{n}_{j}\stackrel{\xaf}{{u}_{i}{u}_{k}}{n}_{l}{n}_{k}}{1+3\stackrel{\xaf}{{u}_{k}{u}_{l}}{n}_{l}{n}_{k}/2k}\right]$. (8)

In these equations, *n _{i}*

${D}_{ij}^{t}=\frac{\partial}{\partial {x}_{k}}\left[{C}_{s}\frac{k}{\epsilon}\left(\stackrel{\xaf}{{u}_{i}{u}_{l}}\frac{\partial \stackrel{\xaf}{{u}_{j}{u}_{k}}}{\partial {x}_{l}}+\stackrel{\xaf}{{u}_{j}{u}_{l}}\frac{\partial \stackrel{\xaf}{{u}_{k}{u}_{i}}}{\partial {x}_{l}}+\stackrel{\xaf}{{u}_{k}{u}_{l}}\frac{\partial \stackrel{\xaf}{{u}_{i}{u}_{j}}}{\partial {x}_{l}}\right)\right]$. (9)

Since the exact and modelled terms are of order *y*^{3} 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}\stackrel{\xaf}{{u}_{i}{u}_{j}}\frac{\partial \epsilon}{\partial {x}_{i}}\right)+{C}_{\epsilon 1}\frac{\epsilon}{k}\stackrel{\u02dc}{P}-{C}_{\epsilon 2}\frac{\epsilon \stackrel{\u02dc}{\epsilon}}{k}+\xi $, (10)

where

$\xi ={f}_{w2}\left(-N\frac{\epsilon \stackrel{\u02dc}{\epsilon}}{k}+M\frac{{\stackrel{\xaf}{\epsilon}}^{2}}{k}-L\frac{\epsilon}{k}\stackrel{\u02dc}{P}\right)$, (11)

is the near-wall correcting function and decays to zero away from the wall,
$\stackrel{\u02dc}{\epsilon}=\epsilon -2\nu {\left(\partial \sqrt{k}/\partial {x}_{i}\right)}^{2}$ and
$\stackrel{\xaf}{\epsilon}=\epsilon -2\nu k/{x}_{2}^{2}$ are reduced *ε*, and *x*_{2} is the normal coordinate. The damping function in (11) is given by
${f}_{w2}=\mathrm{exp}\left[-{\left({\mathrm{Re}}_{t}/40\right)}^{2}\right]$ and the model constants are *C*_{1} = 3.4, *C _{2}* = 4.2,
${C}_{1}^{\ast}$ = 1.8,
${C}_{3}^{\ast}$ = 1.3,

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.

${\stackrel{\xaf}{S}}_{u}=-\rho \left(\frac{\stackrel{\xaf}{{u}_{E}^{2}}-\stackrel{\xaf}{{u}_{W}^{2}}}{2}\right)\Delta y\Delta z-\rho \left(\frac{{\stackrel{\xaf}{uv}}_{N}-{\stackrel{\xaf}{uv}}_{S}}{2}\right)\Delta x\Delta z-\rho \left(\frac{{\stackrel{\xaf}{uw}}_{T}-{\stackrel{\xaf}{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
$\stackrel{\xaf}{{u}_{e}^{2}}=\left(\stackrel{\xaf}{{u}_{E}^{2}}+\stackrel{\xaf}{{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 {\stackrel{\xaf}{{u}_{i}{u}_{j}}}^{p}}{\partial {x}_{j}}$, (13)

where
${\stackrel{\xaf}{{u}_{i}{u}_{j}}}^{p}=\stackrel{\xaf}{{u}_{i}{u}_{j}}+2{\nu}_{t}{S}_{ij}$ is the pseudo stress tensor and *ν _{t}* is evaluated from

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,

$-\stackrel{\xaf}{{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]\\ \text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}+\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

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,

${\stackrel{\xaf}{S}}_{u}=-\rho \left(\stackrel{\xaf}{{u}_{e}^{2}}-\stackrel{\xaf}{{u}_{w}^{2}}\right)\Delta y\Delta z-\rho \left(\stackrel{\xaf}{u{v}_{n}}-\stackrel{\xaf}{u{v}_{s}}\right)\Delta x\Delta z-\rho \left(\stackrel{\xaf}{u{w}_{t}}-\stackrel{\xaf}{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,

$\stackrel{\xaf}{{u}_{p}^{2}}={\left(\frac{\Sigma {A}_{nb}\stackrel{\xaf}{{u}_{nb}^{2}}+\stackrel{\xaf}{S}}{{A}_{p}}\right)}_{P}-{\left({\Gamma}_{U}^{11}\frac{{U}_{e}-{U}_{w}}{\Delta x}\right)}_{P}$, (18a)

$\stackrel{\xaf}{u{v}_{p}}={\left(\frac{\Sigma {A}_{nb}\stackrel{\xaf}{u{v}_{nb}}+\stackrel{\xaf}{S}}{{A}_{p}}\right)}_{P}-{\left({\Gamma}_{U}^{12}\frac{{U}_{n}-{U}_{s}}{\Delta y}\right)}_{P}$, (18b)

$\stackrel{\xaf}{u{w}_{p}}={\left(\frac{\Sigma {A}_{nb}\stackrel{\xaf}{u{w}_{nb}}+\stackrel{\xaf}{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
$\stackrel{\xaf}{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,

$\stackrel{\xaf}{{u}_{e}^{2}}={\langle \frac{\Sigma {A}_{nb}\stackrel{\xaf}{{u}_{nb}^{2}}+\stackrel{\xaf}{S}}{{A}_{p}}\rangle}_{e}-{\langle {\Gamma}_{U}^{11}\rangle}_{e}\frac{{U}_{E}-{U}_{P}}{\delta {x}_{e}}$, (19a)

$\stackrel{\xaf}{u{v}_{n}}={\langle \frac{\Sigma {A}_{nb}\stackrel{\xaf}{u{v}_{nb}}+\stackrel{\xaf}{S}}{{A}_{p}}\rangle}_{n}-{\langle {\Gamma}_{U}^{12}\rangle}_{n}\frac{{U}_{N}-{U}_{P}}{\delta {y}_{n}}$, (19b)

$\stackrel{\xaf}{u{w}_{t}}={\langle \frac{\Sigma {A}_{nb}\stackrel{\xaf}{u{w}_{nb}}+\stackrel{\xaf}{S}}{{A}_{p}}\rangle}_{t}-{\langle {\Gamma}_{U}^{13}\rangle}_{t}\frac{{U}_{T}-{U}_{P}}{\delta {z}_{t}}$, (19c)

where $\langle \text{\hspace{0.05em}}\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 \text{\hspace{0.05em}}\rangle $ terms are reduced to

${\langle \frac{\Sigma {A}_{nb}\stackrel{\xaf}{{u}_{nb}^{2}}+\stackrel{\xaf}{S}}{{A}_{p}}\rangle}_{e}=\frac{1}{2}\left[{\left(\frac{\Sigma {A}_{nb}\stackrel{\xaf}{{u}_{nb}^{2}}+\stackrel{\xaf}{S}}{{A}_{p}}\right)}_{E}+{\left(\frac{\Sigma {A}_{nb}\stackrel{\xaf}{{u}_{nb}^{2}}+\stackrel{\xaf}{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}\stackrel{\xaf}{{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}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\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}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\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}\stackrel{\xaf}{{u}_{nb}^{2}}+\stackrel{\xaf}{S}}{{A}_{p}}\rangle}_{w}-{\langle \frac{\Sigma {A}_{nb}\stackrel{\xaf}{{u}_{nb}^{2}}+\stackrel{\xaf}{S}}{{A}_{p}}\rangle}_{e}\right]\Delta y\Delta z$, (22a)

${S}_{U}^{12}=\left[{\langle \frac{\Sigma {A}_{nb}\stackrel{\xaf}{u{v}_{nb}}+\stackrel{\xaf}{S}}{{A}_{p}}\rangle}_{s}-{\langle \frac{\Sigma {A}_{nb}\stackrel{\xaf}{u{v}_{nb}}+\stackrel{\xaf}{S}}{{A}_{p}}\rangle}_{n}\right]\Delta x\Delta z$, (22b)

${S}_{U}^{13}=\left[{\langle \frac{\Sigma {A}_{nb}\stackrel{\xaf}{u{w}_{nb}}+\stackrel{\xaf}{S}}{{A}_{p}}\rangle}_{b}-{\langle \frac{\Sigma {A}_{nb}\stackrel{\xaf}{u{w}_{nb}}+\stackrel{\xaf}{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 \stackrel{\xaf}{{u}_{P}^{2}}+2\left({\gamma}_{1}-{f}_{w1}{\gamma}^{*}\right)\rho {k}_{P}\begin{array}{c}\text{\hspace{0.05em}}\\ \text{\hspace{0.05em}}\\ \text{\hspace{0.05em}}\end{array}\\ \text{\hspace{0.17em}}\text{\hspace{0.05em}}+\frac{1}{3}\rho \left(1-{f}_{w1}\right){C}_{1}^{*}\stackrel{\xaf}{{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(\stackrel{\xaf}{{u}_{P}^{2}}\right)}^{2}}{{k}_{P}}\}\frac{\Delta x\Delta y\Delta z}{{A}_{P}^{\stackrel{\xaf}{{u}^{2}}}}\end{array}$ (23a)

$\begin{array}{c}{\Gamma}_{U}^{12}=\{-{\beta}_{1}\rho \stackrel{\xaf}{{u}_{P}^{2}}+\left(1-{\alpha}_{1}+{\alpha}^{*}{f}_{w1}\right)\rho \stackrel{\xaf}{{v}_{P}^{2}}+\left({\gamma}_{1}-{f}_{w1}{\gamma}^{*}\right)\rho {k}_{P}\begin{array}{c}\text{\hspace{0.05em}}\\ \text{\hspace{0.05em}}\\ \text{\hspace{0.05em}}\end{array}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{2}\rho {C}_{3}^{*}{\Pi}_{P}^{1/2}{k}_{P}-\rho \left(1-{f}_{w1}\right){C}_{1}^{*}\frac{{\left(\stackrel{\xaf}{uv}\right)}^{2}}{{k}_{P}}\}\frac{\Delta x\Delta y\Delta z}{{A}_{P}^{\stackrel{\xaf}{uv}}}\end{array}$ (23b)

$\begin{array}{c}{\Gamma}_{U}^{13}=\{-{\beta}_{1}\rho \stackrel{\xaf}{{u}_{P}^{2}}+\left(1-{\alpha}_{1}+{\alpha}^{*}{f}_{w1}\right)\rho \stackrel{\xaf}{{w}_{P}^{2}}+\left({\gamma}_{1}-{f}_{w1}{\gamma}^{*}\right)\rho {k}_{P}\begin{array}{c}\text{\hspace{0.05em}}\\ \text{\hspace{0.05em}}\\ \text{\hspace{0.05em}}\end{array}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{2}\rho {C}_{3}^{*}{\Pi}_{P}^{1/2}{k}_{P}-\rho \left(1-{f}_{w1}\right){C}_{1}^{*}\frac{{\left(\stackrel{\xaf}{uw}\right)}^{2}}{2{k}_{P}}\}\frac{\Delta x\Delta y\Delta z}{{A}_{P}^{\stackrel{\xaf}{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*^{−5}.

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*^{−5}. 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,

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

Figure 3. Comparison of the calculated *C _{f}*

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

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

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.

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

Journals Menu

Contact us

+1 323-425-8868 | |

customer@scirp.org | |

+86 18163351462(WhatsApp) | |

1655362766 | |

Paper Publishing WeChat |

Copyright © 2024 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.