**1. Introduction**

Diffusers are ubiquitous in engineering applications. Usually simple in design, they serve to increase the static pressure of a flow by reducing its velocity, albeit often with significant losses. It has a simple design, but it develops complex three-dimensional flow features. A stable separation bubble forms early in the expanding section of the diffuser and spreads across one of the two expanding walls of the diffuser. The flow eventually reattaches in a straight exhaust duct where further pressure recovery occurs. Because of its simple three-dimensional geometry and the existence of a high quality velocity dataset, this diffuser has become a popular test case for measurement of mild separation and validating numerical simulations. However, accurate measurement and prediction of the pressure-driven separation in diffusers has always been challenging in fluid mechanics. In an experimental setup of the diffuser, the first challenge is achieving spanwise homogeneity of turbulence due to the presence of sidewalls producing strong backflow. The point of separation and the extent of the flow reversal zone in diffusers are particularly sensitive to the inlet condition [1]. In a numerical simulation of the diffuser, the presence of an adverse pressure gradient and the formation of an unsteady separation bubble make this flow very sensitive and difficult to predict with numerical means.

Two canonical laboratory incompressible diffuser flows have emerged as standard test-flows in the past: the diffuser studied by Azad [2] and the diffuser studied by Obi et al. [3]. The diffuser studied by Obi et al. [3] is selected in this paper since it has several desirable features. Firstly, the fully developed channel flow is utilized as the inlet condition. Flow in a duct with smooth parallel walls was fully investigated by Kim et al. [4]. For validating the simulation of a separating flow in diffusers. it is vitally significant to know accurately the statistics of the upstream flow because the separation bubble is so sensitive to upstream conditions. Secondly, the flow in this diffuser has rich flow physics, such as pressure-driven separation from a smooth wall, subsequent reattachment, and redevelopment of the downstream boundary layer, which are challenges for large eddy simulation (LES) with subgrid scale mode (SGS) models. Thirdly, based on the mean center velocity *Uc* and the inlet duct of height *H*, the Reynolds number of the inlet channel flow is 7488. The corresponding Reynolds number based on the friction velocity is 407. A direct numerical simulation (DNS) of a diffuser flow is feasible at this Reynolds number, thus the DNS data can be used as supplementary benchmarks for assessing different SGS models [5]. Moreover, this Reynolds number is high enough that the flow is not sensitive to this parameter [6].

Actually, the diffuser studied by Obi et al. has been applied as a test flow in Reynolds-averaged Navier–Stokes (RANS) simulations, LES studies and DNS. Obi et al. [3] demonstrated that the standard *k* − *-* model fails to accurately predict the separation bubble in the diffuser. Durbin [7] proposed the *k* − *-* − *ν*2 model for separated flow and validated his new model in asymmetric plane diffuser with satisfactory accuracy in comparison with measurement. Iaccarino [8] investigated the turbulent flow in an asymmetrical two-dimensional diffuser using three commercial CFD codes: CFX, Fluent, and star-CD. El Behery and Hamed [9] presented a comparative study of turbulence models performance for the separating flow in a planar asymmetric diffuser, in which the steady RANS equations for turbulent incompressible fluid flow and six turbulence closures are used. Schneider et al. [10] performed the LES and RANS calculations for two asymmetrical three-dimensional diffusers and validated LES with wall function delivering the results within the accuracy of experimental data. Abe and Ohtsuka [11] investigated high Reynolds-number complex turbulent flows in a 3D diffuser using LES and hybrid LES/RANS models. Kaltenbach et al. [6] elaborately studied the flow in a planar asymmetric diffuser using LES with the dynamic Smagorinsky model. The overall good agreemen<sup>t</sup> of simulation results and measurement is obtained and it is found that the SGS model plays a significant role for both mean momentum and turbulent kinetic energy balances. Schluter et al. [12] compared LESs with no subgrid model, the standard Smagorinsky model, the dynamic Smagorinsky model and the dynamic localization model in a separated plane diffuser and demonstrated the dynamic localization model performing the best agreemen<sup>t</sup> with measurement. Kobayashi et al. [13] tested their own coherent structures model through some complex geometries in which the asymmetric plane diffuser is included. Taghinia et al. [14] developed a one-equation subgrid scale model with a variable eddy-viscosity coefficient for LES and validated the model in complex separated and reattaching turbulent flows, i.e., the turbulent flow through an asymmetric planar diffuser. Shuai and Agarwal [15] proposed a new one-equation eddy-viscosity model from the two-equation k-kL model. The new turbulence model is used to simulate an asymmetric plane diffuser and improved the accuracy of the flow simulations compared to the one-equation Spalart–Allmaras model. DNS of turbulent flow through an asymmetric plane diffuser was performed by Ohta and Kajishima [5] using a high-order finite difference method. The DNS results from this type of field can be used as one of the benchmarks for numerical simulation schemes and SGS models.

With increasing computing power, LES is more widely used in three-dimensional, unsteady complex flows. In LES technique, turbulent motions are separated into large-scale and small-scale contributions by using a filtering operation in which the large-scale fluid motions are directly calculated, whereas the unresolved SGS motions are modeled. The outcome of an actual LES therefore depends on the quality of the SGS model, as well as the accuracy of the numerical method used to solve the equations for the resolved scales. The important interactions between the resolved large eddies and unresolved SGS eddies, in the case of fully developed isotropic homogeneous turbulent flows, can be reduced to be seen as that of the energy transfers by omitting a portion of information included in the

SGS eddies, for instance, the structural information associated to the anisotropy. It has been shown [16] that the energy transfer between SGS eddies and large eddies mainly exhibit two mechanisms: a forward energy transfer from large eddies to the SGS eddies and a backward transfer to the resolved scales, which, it seems, is much weaker in intensity. The SGS models are responsible for describing the desired dissipation or energy production effects.

In the present study, we focus on the influence of SGS models for LES of the separated turbulent flow in an asymmetric diffuser on the relative coarse mesh through using a high-order finite difference method to solve the equations for the resolved scales. Four existing SGS models, namely the Smagorinsky model [17], the dynamic Smagorinsky model [18,19], the one-equation model [20] and the one-equation dynamic model [21], and a new dynamic one-equation SGS model are investigated in diffuser flow. The Smagorinsky model, which can give a moderately accurate magnitude of energy transfer, is a simple and robust eddy viscosity model; however, it uses a priori model parameter. In the dynamic Smagorinsky model, the model coefficient is dynamically decided from the resolved scales to the unresolved SGS ranges based on an assumption of the scale invariance. However, the clipping procedure or averaging operation is still required in a homogeneous direction or in the global volume of domain. The one-equation model directly uses the information concerned with the SGS motions, i.e., the SGS kinetic energy, to decide the eddy viscosity, in which a transport equation of the SGS kinetic energy is solved and local equilibrium hypothesis of previous models can be removed. The one-equation dynamic model is a modification of the standard one-equation model, in which the production of SGS kinetic energy and the energy loss in large-scale portion are treated with different dynamic mechanisms, i.e., the production term in SGS kinetic energy transport equation is solved using a SGS model based on the resolved scales such as the dynamic Smagorinsky model, while the eddy viscosity in the filtered equation of motion is determined from information of unresolved scales—the transport equation of SGS kinetic energy. Since the test filtering operation is required in the one-equation dynamic model, all undesirable features and inconsistencies related to the test filtering operation are retained. The computation cost of solving an additional transport equation, as well as a dynamic procedure over test filter in the one-equation dynamic model is relatively huge. Thus, these properties still limit its application for the simulation of turbulent flows in complex geometries. It is necessary to continue the search for a new one-equation dynamic model that performs as well as the one-equation dynamic model, while the new model does not require any test filter and is not more expensive in terms of computational cost than the standard one-equation model.

Therefore, the goal of this study is to develop a new one-equation dynamic SGS model and examine the performance of the new SGS model and four existing SGS models on the flow prediction of pressure-driven separation in a diffuser by using a high-order finite difference method to solve the equations for the resolved scales. In addition, we examine the choice of mesh resolution of an asymmetric plane diffuser.

## **2. SGS Models**

Applying a filter with scale Δ, and assuming the filtering operations are commuting with the operations of differentiation, the filtered Navier–Stokes (N-S) equation for LES of incompressible flows can be given as:

$$\frac{\partial \overline{u}\_i}{\partial t} + \frac{\partial \overline{u}\_i \overline{u}\_j}{\partial \mathbf{x}\_j} = -\frac{1}{\rho} \frac{\partial \overline{p}}{\partial \mathbf{x}\_i} + \nu \frac{\partial^2 \overline{u}\_i}{\partial \mathbf{x}\_j \partial \mathbf{x}\_j} - \frac{\partial \mathbf{r}\_{ij}}{\partial \mathbf{x}\_j}, \quad \frac{\partial \overline{u}\_i}{\partial \mathbf{x}\_i} = 0. \tag{1}$$

where the SGS stress tensor is defined as *τij* = *uiuj* − *uiuj*. A index notation is utilized to represent the components, in which the coordinates *x*1, *x*2 and *x*3 are denoted as streamwise, wall-normal and spanwise directions, respectively, and the corresponding velocity components are presented as *u*, *v* and *w*. Throughout this paper, summation convention is implied for repeated index.

As Boussinesq proposed, the SGS stress can be modeled as:

$$-\frac{\partial \mathbf{r}\_{ij}^{a}}{\partial \mathbf{x}\_{i}} \equiv -\frac{\partial}{\partial \mathbf{x}\_{i}} \left(\boldsymbol{\tau}\_{ij} - \frac{1}{3}\boldsymbol{\tau}\_{kk}\delta\_{ij}\right) = \frac{\partial}{\partial \mathbf{x}\_{i}} \left(2\boldsymbol{\nu}\_{\text{sys}}\boldsymbol{\mathsf{S}}\_{ij}\right),\tag{2}$$

in which *τaij* is the anisotropic component of *<sup>τ</sup>ij*; *Sij* is the characteristic filtered rate-of-strain tensor and is defined as (*∂ui*/*∂xj* + *<sup>∂</sup>uj*/*∂xi*)/2; and the coefficient of proportionality *<sup>ν</sup>sgs* is the SGS viscosity that remains to be evaluated by a SGS model.

#### *2.1. Smagorinsky Model (SM)*

Using the local equilibrium assumption, i.e., the dissipation rate of SGS energy is in balance with the production rate, the Smagorinsky model can be obtained for LES:

$$\nu\_{\text{sgs}} = (\mathbb{C}\_s \Lambda)^2 |\overline{\mathbb{S}}| \,\tag{3}$$

in which *Cs* is a constant and taken to be 0.1 in this study. |*S*| denotes the norm of the rate-of-strain tensor and is defined as 2*SijSij*. For wall-resolved LES, to make sure that modeled SGS stresses exhibit the near-wall behavior, a damping function *fs* has to be incorporated with the Smagorinsky model:

$$\nu\_{s\mathcal{S}^s} = (\mathbb{C}\_s f\_s \Lambda)^2 |\overline{\mathcal{S}}|. \tag{4}$$

The van Driest function is employed as the damping function:

$$f\_s = 1 - \exp\left(-\frac{y^+}{A^+}\right),\tag{5}$$

in which the non-dimensional constant *A*<sup>+</sup> ≈ 25 and *y*<sup>+</sup> = *yuτ*/*<sup>ν</sup>*. For the separated flow in an asymmetric diffuser, it would be challenging to determine the friction velocity *uτ* in the vicinity of a separation point due to *uτ* = 0. This is why few researchers try to use the Smagorinsky model for LES to simulate the diffuser flow. Even though the Smagorinsky model is applied in diffuser, its results are far from satisfactory. For example, Schluter et al. [12] found the flow separation was predicted on the upper wall instead on the inclined wall using the Smagorinsky model. In this study, an interesting method that *uτ* is determined by a linear interpolation is used in Smagorinsky model.

The two-dimensional schematic of Obi et al. diffuser is shown in Figure 1, in which the expansion ratio and expansion part of the diffuser is 4.7 and 21*H*, respectively. A Cartesian coordinate system with the origin on the upper wall where the inlet channel wall and the deflected wall form a corner is used to define x in the streamwise direction and y in the downward wall-normal direction. Then, *y*<sup>+</sup> in diffuser can be approximately defined as (see Appendix A):

$$y^+(x,y) = \begin{cases} \text{Re}\_{\text{7}} 0 \cdot \min(y, 1-y) / H, & x \le 0 \\ \text{Re}\_{\text{7}} 0 \cdot (1 - 3.7 / 4.7 \cdot x / H) \cdot \min[y, Y(x) - y] / H, & 0 < x \le 21H \\ \text{Re}\_{\text{7}} 0 / 4.7 \cdot \min(y, 4.7H - y) / H, & x > 21H. \end{cases} \tag{6}$$

Here, *Reτ*0 denotes the Reynolds number based on the friction velocity found in the inlet of duct. *<sup>Y</sup>*(*x*) means the local width of expansion in diffuser.

**Figure 1.** Schematic of an asymmetric plane diffuser.

#### *2.2. Standard Dynamic Smagorinsky Model (DSM)*

Applying the test filter on the grid-filtered N-S equations, the Germano identity can be defined as

$$L\_{\rm ij} = T\_{\rm ij} - \widetilde{\mathbf{r}}\_{\rm ij} = \overline{\overline{\mathbf{u}}\_{\rm i}} \overline{\overline{\mathbf{u}}\_{\rm j}} - \overline{\overline{\overline{\mathbf{u}}}} \,\, \overline{\overline{\mathbf{u}}}\_{\rm j} \tag{7}$$

where *Lij* can be calculated based on the resolved scales, *Tij* = *u iuj* − *u iu j* represents the residual turbulent stress at a test-filter scale Δ , and can be given as

$$T\_{ij} - \frac{1}{3} \delta\_{ij} T\_{kk} = -2C\tilde{\Lambda}^2 |\widetilde{\mathfrak{F}}| \widetilde{\mathfrak{F}}\_{ij}.\tag{8}$$

On substituting Equations (2) and (8) into Equation (7) and assuming Δ and C are constant inside the test filter, an equation for determining C is obtained:

$$L\_{ij} - \frac{1}{3} \delta\_{ij} L\_{kk} = -2\mathbb{C}\overline{\Delta}^2 M\_{ij\star} \tag{9}$$

where

$$M\_{i\bar{j}} = \frac{\widetilde{\Delta}^2}{\overline{\Delta}^2} |\widetilde{\overline{\mathcal{S}}}| \widetilde{\overline{\mathcal{S}}}\_{i\bar{j}} - |\overline{\overline{\mathcal{S}}}| \overline{\overline{\mathcal{S}}}\_{i\bar{j}}.\tag{10}$$

Minimization of the error of Equation (9) over all independent tensor components [19], as well as over some averaging region of statistical homogeneity, leads to

$$C = -\frac{1}{2\overline{\Delta}^2} \frac{\langle L\_{ij}M\_{ij} \rangle}{\langle M\_{ij}M\_{ij} \rangle}.\tag{11}$$

#### *2.3. Standard One-Equation Model (OM)*

The SGS eddy viscosity *<sup>ν</sup>sgs* is represented as *<sup>ν</sup>sgs* = *<sup>C</sup>ν*Δ*νksgs* by the dimensional analysis, where *ksgs* = (*uiui* − *uiui*)/2 = *<sup>τ</sup>ii*/2, i.e., the SGS kinetic energy. *ksgs* is evaluated by solving an evolution equation

$$\begin{split} \frac{\partial k\_{\rm sgs}}{\partial t} + \overline{u}\_{\dot{f}} \frac{\partial k\_{\rm sgs}}{\partial x\_{\dot{f}}} &= \left. - \right| \tau\_{\rm lj} \overline{S}\_{\rm ij} - \mathbb{C}\_{\rm r} \frac{k\_{\rm sgs}^{3/2}}{\Delta} - \varepsilon\_{\rm w} \\ &+ \quad \frac{\partial}{\partial x\_{\dot{f}}} \left[ \left( \mathbb{C}\_{\rm d} \Delta\_{\rm V} \sqrt{k\_{\rm sgs}} + \nu \right) \frac{\partial k\_{\rm sgs}}{\partial x\_{\dot{f}}} \right], \end{split} \tag{12}$$

which was theoretically derived by the authors of [22–24]. For meeting the correct asymptotic behavior to the wall, an additional modification was proposed by Okamoto & Shima [24] for the characteristic length Δ*ν* and additional dissipation term *εw*, i.e.,

$$
\Delta\_V = \frac{\overline{\Delta}}{1 + \mathcal{C}\_k \overline{\Delta}^2 \overline{\mathcal{S}}^2 / k\_{\mathcal{S}\mathcal{S}^s}} \, \tag{13}
$$

$$\varepsilon\_{\rm w} = 2\nu \frac{\partial \sqrt{k\_{\rm sgs}}}{\partial x\_{\rm j}} \frac{\partial \sqrt{k\_{\rm sgs}}}{\partial x\_{\rm j}}.\tag{14}$$

Non-dimensional constants associated with one-equation model are as follow: *Cν* = 0.05, *Cε* = 0.835, *Cd* = 0.10, and *Ck* = 0.08.

#### *2.4. One-Equation Dynamic Model (ODM)*

Kajishima and Nomachi [21] considered that the dynamic procedure (that is, Section 2.2 ) is suitable for the determination of energy transfer from resolved scales to SGS components since this transfer is regarded to be local and instantaneous, whereas the energy loss in large-scale portion is considered from the historic effect of SGS turbulence due to the transport. Thus, the dynamic procedure is applied to the production term in the transport equation of *ksgs*, while the SGS eddy viscosity in the filtered N-S equation is determined by using *ksgs* of Equation (12). The first term in the right-hand side of Equation (12) is given by

$$-\tau\_{i\bar{j}}\overline{\mathfrak{S}}\_{i\bar{j}} = \mathbb{C}\_{s}\overline{\mathfrak{A}}^{2}|\overline{\mathfrak{S}}|^{3} \,. \tag{15}$$

where *Cs* is calculated by the dynamic procedure of Equation (11). It should be noted that *Cs* in the ODM model does not require any averaging over homogeneous direction or ad hoc clipping technique to avoid negative Smagorinsky constants. The negative value for *Cs* is acceptable and it results in the decrease in *ksgs*. However, the computation of the test filter to evaluate the Germano Identity *Lij* and the parameter *Mij* is still required.

#### *2.5. One-Equation Vreman Model (OVM)*

Since the test filtering operation is required in the ODM model, all undesirable features and inconsistencies related to the test filtering operation are retained. The computation cost of solving an additional transport equation, as well as a dynamic procedure over test filter in the ODM model is relatively huge. Thus, these properties still limit its application for the simulation of turbulent flows in complex geometries. It is necessary to continue the search for a new one-equation dynamic model that performs as well as the ODM model, while not requiring any test filter and not being more expensive in terms of computational cost than the standard one-equation model.

Vreman [25] proposed a SGS model for the LES of turbulent shear flows (henceforth, referred to as Vreman model). Vreman model is essentially not more complicated than the Smagorinsky model, but is able to adequately handle turbulent as well as transitional flow. The model is expressed in first-order derivative, and does not involve test filtering, averaging, or clipping procedures. From two test cases of a transitional and turbulent mixing layer at high Reynolds number and a turbulent channel flow, Vreman model is found to be more accurate than Smagorinsky model and as good as the standard dynamic Smagorinsky model. Because of these desirable properties, it seems to be particularly pertinent for incorporating Vreman model into the ODM model to develop a new one-equation dynamic model that is suitable for LES of turbulent flows in complex geometries. In the new model, the production term in the transport equation of *ksgs* (Equation (12)) is represented as:

$$-\left.\tau\_{ij}\overline{\mathbf{S}}\_{ij} = \mathcal{C}\_{vm}^{+} \sqrt{\frac{B\_{\beta}}{\alpha\_{ij}\alpha\_{ij}}} |\overline{\mathbf{S}}|^{2} \,. \tag{16}$$

where

$$B\_{\beta} = \beta\_{11}\beta\_{22} - \beta\_{12}^2 + \beta\_{11}\beta\_{33} - \beta\_{13}^2 + \beta\_{22}\beta\_{33} - \beta\_{23}^2. \tag{17}$$

$$
\beta\_{i\bar{j}} = \Delta\_m^2 \alpha\_{mi} \alpha\_{mj}, \quad \alpha\_{i\bar{j}} = \frac{\partial \overline{u}\_{\bar{j}}}{\partial x\_i} \tag{18}
$$

$$\mathcal{C}\_{vm}^{+} = \begin{cases} \mathcal{C}\_{vm} |\overline{\Omega}| / |\overline{\mathbb{S}}| , & |\overline{\Omega}| / |\overline{\mathbb{S}}| < 1, \\\mathcal{C}\_{vm} , & |\overline{\Omega}| / |\overline{\mathbb{S}}| \ge 1. \end{cases} \tag{19}$$

The eddy viscosity is obtained from solving Equations (12) and (16) simultaneously. We call the new SGS eddy viscosity model as one-equation Vreman model. Note that a fixed coefficient *Cvm*(= 0.025) is used in Vreman model, while in our new one-equation dynamic model a modified coefficient *<sup>C</sup>*+*vm* is introduced for taking excess SGS production rate where both |Ω| and |*S*| are large. Actually, to further improve the performance of Vreman model, Park et al. [26] and You and Moin [27,28] proposed different dynamic procedures to dynamically determine the parameter *Cvm*, in which single-level or double-level test filters are employed.

## **3. Computational Method**

In this study, to make much more persuasive comparison of SGS models performance for diffuser flow and expel the other factors, which can interface the observation of comparison, the aforementioned SGS models in LES and a DNS are numerically implemented through the same set of numerical methods in an asymmetric planar diffuser. The results of DNS are used as benchmarks for SGS models as well. The configuration of the diffuser, as shown in Figure 2, and Reynolds number *Rec* = 7488 based on the mean center velocity *Uc* found in the inlet duct of height *H* match previous experiments performed by Obi et al. [3,29].

**Figure 2.** Three-dimensional schematic of Obi et al. diffuser.

#### *3.1. Domain Size and Boundary Condition*

There are two computational regions sandwiched by two walls, called the driver and the spatially developing region, as given in Figure 2. The expansion ratio of the diffuser is 4.7, the inclination of the wall (the lower wall) is 10 degrees and the spanwide extend is 7.68 H. The driver region length is 15.84 H and the upstream channel, expansion part and downstream extension of diffuser are 3 H, 21 H and 39.36 H, respectively. Both the joins between the expansion part and upstream/downstream parts are rounded with radius of 11.43 H, which are the same as in experiment. A Cartesian coordinate system with the origin on the upper wall where the inlet channel wall and the deflected wall form a corner is used to define x in the streamwise direction, y in the downward wall-normal direction, and z in the spanwise direction.

The flow in two regions are numerically computed simultaneously with periodic and no-slip boundary conditions for the velocity in spanwise direction and wall-normal direction, respectively. There, it should be noted that the boundary condition in the spanwise direction between simulations and the experiment performed by Obi et al. [3] is different, i.e., periodic (homogeneous) boundary condition in simulations but closed (rectangular duct) boundary condition in the experiment. To some extent, this kind of difference will affect the agreemen<sup>t</sup> between the simulations and the experiment, which is specifically discussed in the following. The fully-developed turbulent channel flow that is treated as inflow for the spatially-developing region are generated by the driver region with the periodic boundary condition in streamwise direction. The pressure gradients for driver region in streamwise direction are controlled to keep the flow rate a constant. Velocity profiles at outflow boundary condition for spatially-developing region are determined by a convective outflow condition

$$\frac{d\mu}{dt} = -u\_{\rm env} \frac{d\mu}{d\mathbf{x}'}\tag{20}$$

in which convective velocity *ucvn* is initialized with a constant bulk velocity to keep the flow rate in spatially-developing region. Finally, the pressure gradients at spatially-developing region are obtained automatically.
