Turbulence Modeling Using Z-F Model and RSM for Flow Analysis in Z-SHAPE Ducts

Turbulent ﬂow in Z-shape duct conﬁguration is investigated using Reynolds stress model (RSM) and ζ - f model and compared to experimental results. Both RSM and ζ - f models are based on steady-state RANS solutions. The focus was on regions where the RSM has over- or underpredicted the ﬂow when compared to the experimental results and on regions where there are ﬂow separations and high turbulence. The performance of predicting the ﬂow reattachment length in each model is studied as well. RSM has shown the mean ﬂow velocity proﬁle results match reasonably well with the experiment. Advanced ζ - f turbulence model is introduced as user-deﬁned function (UDF) code and applied to the Z-shape duct. It is found that the turbulent kinetic energy production in ζ equation is much easier to reproduce accurately. Both mean velocity gradient and local turbulent stress terms are also much easier to be resolved properly. The current research has found that not only ζ - f model takes less time to complete the simulation but also the mean ﬂow velocity proﬁle results are in better agreement with the experimental data than the RSM although both are coupled steady-state RANS. ζ - f model numerically resolved both the ﬂow separation and reattachment regions better than the RSM. The current numerical results from ζ - f model are attractive and encouraging for wall-bounded ﬂow applications where ﬂow separation and ﬂow reattachment are important for the ﬂow mechanism.


Introduction
Duct and pipe systems with various curvature shapes have widespread industrial applications, e.g., HVAC systems, heat exchangers, and gas and fluid transport lines (such as water, gas, and oil). Laminar flow through the curved ducts is well studied, and the physics are well understood while the turbulent flow is more complex, and the flow pattern is not recognized universally. is is especially true in duct shapes that are nonconventional such as Z-shape that is considered in this study. e Z-shape ducts can be found in almost all ventilation and air-conditioning duct systems as well as in many other industrial duct and pipe systems.
In open literature, research on turbulent flow in Z-shape ducts is very limited. Salehi et al. [1,2] studied experimentally and numerically the turbulent flow in Z-shape and U-shape ducts. e separation distance, which is the duct length between the two elbows, was varied. ey concluded that the two-equation turbulence models were unable to predict the flow accurately; however, RSM with enhanced wall treatment was able to predict elbow loss coefficients (<15% of error). Sleiti et al. [3] measured detailed velocity profiles in Z-shape with turning radii r/D � 1.5 to understand the flow physics for geometries with bends and provided experimental data to verify the results of computational fluid dynamics. ey compared the measured velocity profiles to RSM and large-eddy simulation (LES) predictions for the effect of separation distance. It was concluded that RSM was able to predict the velocity profiles accurately. However, LES model was unable to predict the velocity profiles. Challenges and opportunities in predicting complex flow in duct fitting losses using different CFD turbulence models were presented in [4]. e authors found that two-equation RANS turbulence models could only predict the flow trends and advanced turbulence models are needed to predict such complex flow.
Several researchers studied other duct geometries with curvatures and secondary flows in stationary and even rotational frame. Rectangular helical ducts were studied in [5] using realizable k-ε eddy-viscosity model (EVM). Due to the curvature of the ducts, as fluid flows through curved tubes, a centrifugal force is generated and a secondary flow induced by the centrifugal force has enhanced the heat transfer rate. Taylor et al. [6] studied laminar and turbulent flow experimentally in 90°-bend duct with a curvature ratio of δ � 3.7. Afterward, Sudo et al. [7,8] worked on 90°-bend ducts (square-sectioned) and pipes (circular-sectioned) with a curvature ratio of δ � 4. According to their results, boundary layer separation did not happen in the bend due to small values of curvature ratios. To investigate the effect of curvature ratio, Ono et al. [9] studied the water flow through two elbows with δ � 2 and 3. Later, Tan et al. [10] performed an independent experiment and showed that boundary layer separation would occur by increasing curvature ratio. Tunstall and Harvey [11] observed a secondary flow pattern in a bent pipe for Re � 4 × M10N ∧ 4, which differed from classical twin vortex pattern. e new flow pattern was characterized by single swirl secondary flow that was unsteady bistable flow and its direction changed from clockwise to counterclockwise randomly at low-frequency timescale. Hellstrom et al. [12] investigated the flow field in an elbow with a curvature ratio of δ � 2 and Reynolds numbers between Re � 2×M10N ∧ 4 and Re � 1.15×M10N ∧ 5. Applying proper orthogonal decomposition (POD) method to experimental snapshot data, they found that the Dean motion is not the most energetic structure of the flow. Instead, a separated secondary flow with random rotation direction in clockwise and counterclockwise direction, called the "swirl switching" mode, is the dominant structure. e flow of hydrogen in complex geometry of tube bundles coaxially and guided radially is investigated in [13] experimentally and numerically using two-equation CFD models. Due to the drastic change in the flow direction and its inherently three-dimensional nature and sharp curvatures, optimizing the geometry resulted in 15% reduction in the average pressure drop.
Rutten et al. [14] studied numerically turbulent flows in two 90°bends with δ � 2 and 3 using LES. ey observed the separation of flow and swirl switching for Re � 27000 and obtained the Strouhal number of 0.2-0.3, independent of curvature ratio. In [15], a comparison between EVM and RSM turbulence models in predicting flow and heat transfer in rotating rib-roughened channels was conducted. e authors concluded that RSM is superior to EVMs in predicting such flow. e effect of high-density ratios and rotation was studied in [16][17][18] using RSM, where the Coriolis and centrifugal forces in U-shape geometry and parallel rotation induced secondary flows in several locations not seen at low-density ratios and rotations.
Several other researchers summarized in Table 1 used RANS, LES, and even DNS modeling to study the flow characteristics in curved ducts. ζ-f model, in comparison with models mentioned in Table 1, is a version of the eddyviscosity model, which solves a transport equation for the velocity scales ratio instead of the equation for stresses, which improves the numerical stability of the model and improvements for nonequilibrium wall flows as will be detailed in the below sections. e model is more robust and less sensitive to nonuniformities and clustering of the computational grid.
It has been noted that none of the studies found in open literature were performed on Z-shape ducts with different separation distances, which is the subject of the current research study.
Turbulent flow in Z-shape ducts is rather complicated and needs to be studied and analyzed thoroughly. Experimental work in [3] showed that such flow is a strong function of the separation distance, and therefore, more studies are needed to understand the turbulence behavior. e accurate numerical simulation of the flow through Z-shape duct depends on the correct prediction of boundary layer transition characteristics. Reynolds stress turbulence model considers more flow physics and models the turbulence redistribution close to the wall which plays an important role within the transition process. Compared to other turbulence models, ζ-f model has greater numerical stability and robustness advantages for computing the flow field. is proposed new version of eddy-viscosity model is based on modification of the elliptical relaxation method implemented on solving steady-state transport equation. Instead of turbulent kinetic energy, k, that was implemented in standard k-ε model, the ζ-f model uses ζ to evaluate eddy viscosity. ζ is referred to velocity scale ratio defined by υ 2 /k, where υ 2 is "wall-normal" velocity scale and k is the turbulence kinetic energy. f is the elliptic function used to model the anisotropic wall effects. Typically, f is solved by elliptic equation based on Helmholtz formulation, which will be described in the next section. υ 2 also represents the velocity fluctuation normal to the streamlines.
is new concept was mainly to make use of the right scaling for better representation of the damping characteristics of turbulent transport near to the wall.
In this paper, the authors will address specific turbulence modeling issues applicable to the flow in Z-shape ducts and will study and analyze the flow turbulence using Reynolds stress model (RSM) and ζ-f model.

Turbulence Models
e most accurate method to simulate the turbulent flow is the direct numerical simulation (DNS) in which the numerical grid is small enough to resolve the smallest eddies. However, the computational cost of DNS is such high that it is not practical for most of the engineering problems. RANS method is popular in many engineering applications; however, the models are not universal, i.e., they cannot predict many complicated flows such as stagnation point flow, boundary layer transition from laminar to turbulence, separation of flow, and curvature effects. Smagorinsky [27] introduced an alternative approach called large-eddy simulation (LES), which resolves the large-scale vortices and models the small-scale ones. eoretically, this approach is more accurate than the RANS modeling. Some for the most recent applications of LES are demonstrated by Liu et al. [28]. ey have investigated the interaction of cavitation-vortex-turbulence flow over a hydrofoil based on modified Schnerr-Sauer model. Compared to experiment, mean flow average velocity profile predicted by LES is in good agreement. Liu and Tan [29] conducted a numerical study of the tip leakage vortex characteristics in a mixed flow pump based on shear stress transport (SST) k-ε method. eir predictions show good agreement with the experiment.
In computational expenses, LES categorized between RANS models and DNS and tries to cover the disadvantages of these approaches. LES does not require very fine meshes as DNS; however, the grid should be finer than that those typically used for RANS calculations. One of the great RANS models is ζ-f model which is proposed by Hanjalic et al. [30] and Popovac and Hanjalic [31] and will be used in the current study. In the following sections, description of the turbulence models used in this study is provided.

ζ-f Model.
Originally, the ζ-f model is based on v 2 − f model which is proposed by Durbin [32] by replacing ζ � v 2 /k. Both ζ-f model and v 2 − f model are considered similar to the standard k-ε model, except both of them incorporated some near-wall turbulence anisotropy and nonlocal pressure strain effects.
e new ζ-f formulation is introduced and described as follows: where ere is no need to make use of wall blending functions as this method resolved the transport equations all the way to solid wall . Compared to v 2 − f model, ζ can be computed more efficiently with the new ζ-f arrangement [30]. It is found that the turbulent kinetic energy production P term in ζ equation is much easier to reproduce accurately. Both mean velocity gradient and local turbulent stress terms are also much easier to be resolved properly. In the near-wall region, ζ is found proportional to the distance squared, ζ ∝ y 2 . e wall boundary condition for the ζ-f model is simply described as follows: (2) e equation above is much simple compared to v 2 − f model. Together, these combined benefits provide robustness improvement and higher efficiency for computation, which is the main goal of this dissertation study.
e following elliptical f formulation is derived using quasilinear model for equation (1) by applying P � 0 and X � 0. Both L and Τ are defined as length and timescale, respectively. For simplicity reason, (3) represents the elliptic relaxation equation for f. It is noted that the last term can be neglected since C 4 /3-C 5 � 0.008: Due to special treatment of near-wall turbulence anisotropy and nonlocal pressure strain effects, ζ-f model yields promising results for many flow phenomena including separated flow. e model is considered a four-equation model which solves transport equations for turbulent kinetic energy k, turbulent dissipation rate ε, velocity scale ratio ζ, and elliptic function f parameters as follows [30,31]: where e turbulent eddy viscosity, timescale, and length scale are computed from e closure coefficients are as follows: C μ is the same as v 2 − f model of Durbin. e parameter values above are recommended for numerical stability reason based on Durbin's realizability constraints.

Reynolds Stress Model.
e model solves additional transport equations for six independent Reynolds stresses, requiring one equation for turbulent dissipation; isotopic eddy-viscosity assumption is avoided. e model has shown excellent results with complex flows. Models that solve transport equations are called second-moment closures. e RSM has the advantage to adopt high Reynolds number turbulent flows, generally superior over one-and twoequation models. It is important to implement a secondmoment closure model when simulating stress induced secondary flows in ducts [33]. e boundary condition of RSM requires user-defined Reynolds stresses and dissipation rate. In which case, there are equations that could be used to have an estimation value from the characteristic length or the hydraulic diameter [33]. Wall functions available with RSM are standard, nonequilibrium, and enhanced wall functions. e standard wall function is investigated by Launder and Spalding [34]. It is the default option in FLUENT and is good to adopt with a wide range of bounded flows. However, it has the disadvantage of poor predictions in cases with separation and reattachment flows such as U-bends. e nonequilibrium model produces better results in cases with flow separations due to it taking an account of pressure gradients. e nonequilibrium wall function can be used with RSM and k-ε models [33].

Numerical Simulation Approach
e Z-shape duct with diameter D � 12 is shown in Figure 1. e geometry consists of three parts, a 15D first horizontal section, a 10D vertical section, and a 10D second horizontal section where D is the diameter of the pipe. Inlet type is set to "velocity inlet" with a value such that the flow Reynolds number becomes R e � 3.5 × 10 5 . e outlet was set to pressure outlet with zero-gauge pressure. A computational grid with hexahedral cell type is produced with 2.1 × 10 6 cells and it is depicted in Figure 2.
A UDF script was developed for the new turbulence model interface with FLUENT solver. e custom source terms functions for transport equations from equations (3)-(7) are defined in the UDF script. is includes turbulent kinetic energy, k, turbulent dissipation rate, ε, "wallnormal" velocity scale, v 2 , and elliptic function, f.

Mesh Independence Study
Several mesh resolution variations were applied to the computational domain. e purpose of this mesh independence study is to ensure the flow inside the pipe is resolved accurately and the fluid flow physics is captured sufficiently. Since multiple simulations are required to be completed on the final optimal mesh resolution, a mesh independence study is desired for adopting the coarser converged mesh from multiple meshes to save computational time than suggesting a fine mesh that would produce the same results but will be computational demanding. is process requires several mesh resolutions to be studied and attempted systematically. e final optimal mesh resolution is determined based on the comparison of the converged numerical results from the previous coarse mesh. Such comparison process would determine the final optimal mesh is adequate to resolve the turbulence flow physics.
As shown in Figure 3, the mesh is refined in the bending geometry by a further 30% compared to the bulk mesh region. e mesh at the near-wall region is also refined to keep y + as low as possible.
e range of y+ in all cases studied in this paper was between 0.6 and 1. Results of turbulence kinetic energy for a single case were analyzed at section x/D � 5 at duct centerline r/D � 0 for different grid sizes as can be seen in Table 2. Medium mesh was adopted since TKE is within acceptable tolerance (±2%) range and it has less 2.1 × 10 6 number of cells than the fine mesh which will save computational time.

Results and Discussion
is section presents and discusses the validation of numerical predictions compared to the experimental data from 4 Journal of Engineering Sleiti et al. [3]. Section x/D � 1 represents the section right after first elbow transition. Based on basic principles of fluid mechanics, after flow making the first turn at the elbow, the flow velocity distribution becomes highly nonuniform due to local flow separation. erefore, at x/D � 1 section, such flow pattern has shown highly asymmetric in the E-W orientation. Higher velocity was found near the outboard (negative) radial locations of the section and lower velocities within the inner (positive) radial location. e mean flow velocity reduces due to flow separation and increases due to flow reattachment. Figures 4 to 8 show the results of the normalized mean flow velocity profiles, and Figure 9 shows the normalized velocity contours at different separation distances, x/D. Figure 4 illustrates the results of normalized mean flow velocity profiles at x/D � 1. e velocity profiles at the E-W direction show that ζ-f turbulence model results agree well with experiment. Within the flow separation region, the ζ-f model prediction is found much closer to experiment trend compared to RSM. RSM has poorly underpredicted the trend within this flow separation region. e ζ-f model also predicts well near the flow reattachment region. Similarly, RSM is also found to predict well at the flow reattachment region. RSM region of lower velocities is larger than ζ-f (Figure 9), which explains the RSM underprediction. e velocity profiles at NS (y) direction show that ζ-f and RSM results are in good agreement with experiment.
When the flow enters the straight pipe section, flow distribution is found transitioning from nonuniform type to uniform as flow moves downstream towards outlet. e following sections are analyzed at different locations along the straight pipe locations at x/D � 3, 5, and 7. Figure 5 depicts the results of normalized velocity profiles at x/D � 3. e velocity profiles at EW (z) direction indicate that ζ-f turbulence model produces very close results to experiment. Moreover, the results of RSM are in good agreement with experiment at most locations. RSM overpredicts at the flow reattachment region. Similarly, the velocity profiles at NS (y) direction demonstrate that the predictions of ζ-f models are better. Figure 6 displays the results of normalized velocity profiles at x/D � 5. e velocity profiles at both EW (z) and NS (y) directions show that ζ-fmodel and RSM capture the flow physics very well.
At further downstream in the straight pipe, it is found that the flow velocity becomes much uniform due to fully turbulent flow development. Figure 7 shows results at the        shedding in the flow separation region near bending elbow geometry is found to be the key flow mechanism dictating downstream flow development. Turbulence eddy dissipation downstream accompanied with the flow mixing mechanism is found to produce complex flow structures with highly nonhomogeneous and nonisotropic behaviors.
e ζ-f model is newly developed in this thesis and has not been tested before for this Z-shape pipe. erefore, it is believed the benchmarking effort of this research would contribute to further exploration and use of this ζ-f model in similar or other applications. e predictions of RSM showed good accuracy. But there are some discrepancies in the results. e RSM is efficient inside the boundary layer but has excessive diffusion in the separated region.

Conclusion
Turbulent flow predictions studied for Z-shape duct configuration using Reynolds stress model (RSM) and ζ-f model. Both RSM and ζ-f models used in this research study are based on steady-state RANS solutions. e focus was on regions where the RSM has over-or underpredicted the flow when compared to the experimental results and on regions where there are flow separations and high turbulence. e performance of predicting the flow reattachment length in each model is studied as well.
e current RSM numerical analysis and validation results using Z-shape ducts have shown that the mean flow velocity profile results match reasonably well with experiment. However, RSM of the current study is coupled to steady-state RANS equations and therefore assumed the temporal effect of turbulence flow is minimal and turbulence flow is steady state/frozen state.
Advanced ζ-f turbulence model is developed and written with UDF code and is applied to the Z-shape duct. is newly developed ζ-f model is never investigated before for this duct configuration application. e ζ-f model is used instead of v 2 − f model because ζ can be computed more efficiently with the new ζ-f arrangement. It is found that the turbulent kinetic energy production in ζ equation is much easier to reproduce accurately. Both mean velocity gradient and local turbulent stress terms are also much easier to be resolved properly. e current research has found that ζ-f model not only takes less time to complete the simulation but also the mean flow velocity profile results found better agreement with experimental model than RSM although both are coupled steady-state RANS.
e ζ-f model numerically resolved both the flow separation and reattachment regions better than RSM. e current numerical results from ζ-f model are attractive and encouraging for wallbounded flow applications where flow separation and flow reattachment are important for the flow mechanism.

Data Availability
Experimental data used in this research article were obtained from Ahmad Sleiti, Steve Idem, and Mohamed Salehi (2016)

Conflicts of Interest
e authors declare that they have no conflicts of interest.