Large-eddy simulation of the flow in Z-Shape duct

Abstract A numerical study is conducted of the transient flow in a z-shape three-dimensional duct using large-eddy simulation (LES) with dynamic eddy viscosity subgrid-scale (SGS) model with a fully structured grid system. The numerical results of the velocity profiles are quantitatively validated against experimental data for z-shape ducts with various lateral separation distance configurations. The framework of the current LES model has been studied and discussed and the performance of LES in predicting the flow in z-shape ducts as a function of separation distances was evaluated. LES predictions of the mean flow velocity profiles are in good agreement (within the experimental uncertainty) with experimental data for the investigated wide range of l/d configurations. This is attributed to the well-resolved large-scale flow structures. Some slight over-predictions and under-predictions were found at certain separation distances. These numerical errors are due to the limited modeling approach to predict small eddies structures with the current SGS model. The main key features of the flow after the first elbow are also identified as separation and re-attachment regions. Some discrepancies are identified for lateral separation distances at sections x/d = 3 and x/d = 5 inside the flow transition regions. These discrepancies are believed to be inherited from the upstream flow numerical errors that arise in the non-uniform flow mixing regions. The potential remedy includes applying finer mesh resolution and/or higher order spatial discretization to accurately resolve the local velocity gradients and the complex flow structures.

ABOUT THE AUTHOR Mohammed Karbon is currently a PhD candidate in the Department of Mechanical and Industrial Engineering at Qatar University. Karbon expertise is in Computational Fluid Dynamics and Turbulence Modeling working under the supervision of Professor Ahmad K. Sleiti. He has many years of industrial work experience before joining the group of Prof. Ahmad. Karbon is also familiar with several computational skills and commercial software packages including ANSYS FLUENT and others.

PUBLIC INTEREST STATEMENT
Although the computational modeling of incompressible fluid flow has been making advancement in recent years, accurate numerical predictions of high Reynolds number turbulent flow remain challenging for wall-bounded flows, especially in Z-shape geometries. In an attempt to address this challenge, this research presents the modeling framework and detailed numerical results using transient large-eddy Simulation (LES) approach for flow in Z-shape duct configuration. This paper reports a comprehensive validation of the current state-of-the-art LES simulated results against experimental data. The effect of the separation distance between the two elbows of the Z-shape duct on mean-flow velocity at varying duct lateral separation distances is studied extensively. This study calls for more attention to the use of transient LES method to capture high Reynolds number flow physics in complicated geometries.

Introduction
Ducts and duct systems with various curvature shapes have widespread industrial applications, e. g., HVAC systems, heat exchangers, 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 detailed flow physics remain unsolved fully. It is especially true in duct shapes that are non-conventional such as Z-shape that is considered in this study. The Z-shape ducts can be found in all ventilation and air-conditioning duct systems as well as in many other industrial applications.
In the open literature, research on turbulent flow in Z-shape ducts is very limited.  and , studied experimentally and numerically the turbulent flow in Z-shape and U-shape ducts. They defined a separation distance as the duct length separating the center-points of the elbows and this distance was systematically varied. Zerolength pressure loss coefficients were predicted using five different two-equations eddy viscosity models including the standard k-ε, the Realizable k-ε, RNG k-ε, standard k-ω, SST k-ω models, and the Reynolds Stress Model (RSM). Their results have shown that the two-equation turbulence models predicted incorrect trends, while RSM with enhanced wall treatment was generally able to correctly predict elbow loss coefficients with less than 15% of error. Sleiti et al., (2017) measured detailed velocity profiles in Z-shape and U-shape ducts with turning radii r/D = 1.5 to study the physics of the flow in complex geometries and to provide data to verify the accuracy of computational fluid dynamics (CFD) modeling predictions. They compared the measured velocity profiles to RSM and LES predictions for the effect of separation distance. Their results showed that RSM and LES predicted the velocity trends and magnitudes correctly, with more accurate predictions credited to RSM as LES error was very high (16%) at sections that are close to elbows. This high LES error predicted in  is attributed to using the steady-state LES modeling. In the present study, the transient LES modeling approach is used in an attempt to conduct more accurate investigations and to report more LES data for different separation distances to further evaluate LES modeling approach.
Several researchers studied other duct geometries with curvatures (Sleiti & Kapat, 2005) and with other strong flow conditions (Sleiti & Kapat, 2006a;2006b;Sleiti, 2009;Sleiti, 2011;Sleiti et al., 2013). Taylor et al. (1982) studied laminar and turbulent experimentally in 90°-bend duct with a curvature ratio of δ = 3.7. Afterward, Sudo et al. (1998;2001) worked on 90°-bend ducts (squaresectioned) and ducts (circular-sectioned) with 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. (2011) studied the water flow through two elbows with δ = 2 and 3. Later, Tan et al. (2014) performed an independent experiment and showed that boundary layer separation would occur by increasing curvature ratio. Hellström et al. (2013) investigated the flow field in an elbow with curvature ratio of δ = 2 and Reynolds numbers between Re = 20,000 and Re = 115,000. 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. Rütten et al. (2005) studied numerically turbulent flows in two 90° bends with δ = 2 and 3 using LES. They observed the separation of flow and swirl switching for Re = 27,000 and obtained the Strouhal number of 0.2-0.3, independent of curvature ratio. In (Sleiti & Kapat, 2006c), the authors investigated 3d U-bend ducts using Standard k-ε, Realizable k-ε (RKE), RNG k-ε, SST k-Omega (k-ω) and RSM. They concluded that all models over predicted the reattachment length except that the SKE model, RSM and RNG k-ε models have satisfactory results at recirculation zones but require more computational time. Several other researchers used RANS, LES and even DNS modeling to study the flow characteristics in curved ducts. However, none of the studies were performed on Zshape ducts with different separation distances.
Turbulent flow in Z-shape ducts is rather complicated and needs to be studied and analyzed thoroughly. Experimental work in  showed that such flow is a strong function of the separation distance and therefore more studies are needed to understand the turbulence behavior. In this paper, specific turbulence modeling issues applicable to the flow in Z-shape ducts, not reported before, are addressed and analyzed using transient LES simulations. Specific details related to numerical mesh, time step, discretization and solution are provided in context with the problem being considered, i.e. Z-shape ducts.

Turbulence models
The 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 so high that it is not practical for most of the engineering problems. Reynolds Averaged Navier-Stokes (RANS) method is widely used in many engineering applications; however, RANS 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, curvature effects and so on. Smagorinsky (Smagorinsky, 1963) introduced an alternative approach; LES, which resolves the large-scale vortices and models the small-scale ones. Compared to LES, RANS equations are formulated to account for averaged flow quantities with a whole range of turbulence scales being modeled. Therefore, the RANS-based modeling approach requires less computational resources with increased turn-around time. Since LES is required to resolve the energycontaining turbulent eddies in both spatial and temporal domains, the computational resources involved are significantly intensive, particularly for high Reynolds number industrial flow applications. LES computational expense is usually found in between RANS models and DNS approaches. For wall-bounded flow applications, it is possible to employ coarse mesh in the near wall region and combine with near wall function to achieve lowered computational cost. For this reason, it is critical to justify carefully the near wall modeling approach including the type of discretization schemes.

Large eddies simulation model
For incompressible flow, applying the filter operator to the governing equations of conservation of mass and momentum (Navier Stokes equations) (ANSYS INC., 2013), the filtered formulations of the resolved field can be obtained as follows: where the bar means the filtered field and σ ij represents the viscous stress tensor defined as: and τ ij is the Sub-Grid-Scale (SGS) stress tensor obtained from filtering and defined as follows: The SGS stresses that are obtained from the filtering procedure are not known, and they should be modeled using appropriate assumptions. Smagorinsky (Smagorinsky, 1963) proposed to use Boussinesq's hypothesis or eddy viscosity assumption to obtain the SGS stress tensor as follows: where σ t is called SGS eddy viscosity. The term τ kk is isotropic part of the SGS stresses that is not modeled, and it is assumed to be added to the filtered static pressure term. � S ij is the strain rate tensor of resolved scales that is defined as below: Employing the idea of the "mixing-length model" in the RANS modeling approach, Smagorinsky (1963) proposed to use the following model for SGS eddy-viscosity: where L s is the mixing length for SGSs, Δ is the local grid scale and C s is Smagorinsky constant.
The drawback of the Smagorinsky model is that C s is flow-dependent and it is not constant. Researchers reported different values for C S from C S ¼ 0:065 (Kim & Moin, 1982) to C S ¼ 0:250 (Wille & Jones, 1996). In addition, Smagorinsky model is too dissipative and it is not appropriate for some physics such as boundary layer transition. Germano et al. (1991) and Lilly (1992) developed a model that calculated the C s according to the resolved eddies of the flow, which is called dynamic Smagorinsky-Lilly Model. Nicoud & Nicoud (1999) present another subgrid-scale model called, Wall-Adapting Local Eddy-Viscosity (WALE) model in which the eddy viscosity is calculated from the following equation; Where where V is the volume of the computational cell.
WALE subgrid-scale model takes into account the effect of both strain rate tensor and rotation rate tensor of resolved eddies. In addition, the asymptotic near wall behavior of eddies (y 3 ) is modeled correctly. The accurate predictions of wall-bounded turbulent flow are determined by an accurate representation of the flow in the near wall region. Chapman (1979) has studied the impact of grid size required in the LES computational domain with respect to Reynolds number. He estimated the total number of grids is proportional to N total ~ Re c 1.8 . This implies LES modeling is far from realistic at high Reynolds numbers as the extremely high computational cost makes it close to DNS simulation. In particular, the turbulent boundary layer developed near the wall is dominated by the presence of multiscale flow structures with high energetic and dynamic motions. Such scales are becoming smaller as the Reynolds number is increased. This raises the demand for increasing total grid size with increased Reynolds number and therefore near wall problem imposes significant challenges for LES implementation. In the present work, dynamic Smagorinsky wall model is used to obtain the correct SGS viscosity as part of the near wall modeling strategy. The details of the near wall numerical mesh are discussed in the subsequent sections.

Details of the experimental setup used for validation and comparison
The details of the test set up and the measurement data used for validation and comparison purposes are provided in Sleiti et al., 2017) and summarized here for the convenience of the reader. The experimental setup shown in Figure 1 is implemented specifically to study the velocity profiles and pressure losses associated with closecoupled round five-gore elbows. The measurements of pressure loss and volumetric flow rate through the ductwork and fittings were performed in accordance with ASHRAE Standard 120-2008. The elbow pressure loss experiments were preceded by a series of tests designed to evaluate the friction factor of straight ducts. A bellmouth was mounted at the entrance of the ductwork to ensure uniform inlet flow. Pressure taps soldered to the ducts were employed to measure the pressure loss at specific distances prescribed in Standard 120-2008.
The turning radii of the elbows were fixed at r/D = 1.5. The ducts upstream and downstream of the test section were each 1.20 m (4 ft) in length and were connected by slip couplings. For close coupled elbows the intermediate length, measured from the exit plane of the upstream elbow to the entrance plane of the downstream elbow, was varied from 0 m (0 ft) to 3.05 m (10 ft) in increments of 0.60 m (2 ft).
A 30 hp centrifugal fan provided air flow through the test apparatus in the forced draft mode. A multiple-nozzle chamber in compliance with the requirements of Standard 120-2008 was used to measure the volume flow rate through the test setup. Pressure loss measurements over the test section were performed using a liquid-filled micromanometer having a measurement accuracy of ±0.025 mm (0.001 in). The static gage pressure upstream and downstream of the test section was measured with electronic manometers having a readability of ±0.25 mm (0.01 in). The nozzle pressure loss was measured using a micromanometer having a scale readability of ±0.025 mm (0.001 in.). Static pressures were measured in the chamber (to assess air density) using an electronic manometer having the scale readability of ±0.25 mm (0.01 in.). A VFD was used to control the flow rate through the test section. The air temperature in the nozzle chamber was measured using a mercury thermometer having a scale readability of ±0.6°C (1°F). The dry-bulb temperature and wet-bulb temperature of the ambient air was measured using an aspirated psychrometer, with an accuracy of ±0.6°C (1°F). Ambient pressures were measured with a Fortintype barometer, with an accuracy of ±0.25 mm (0.01 in.) of mercury. At least eight test points, evenly spaced over the range of test velocities, were obtained for each test apparatus.
For the Z-shape close-coupled intermediate section lengths and apparatus configurations, detailed velocity profile measurements were performed using a five-hole directional velocity probe on the 305 mm (12 in.) diameter test apparatus. One traverse plane was located one duct diameter upstream of the first elbow. Another measurement plane was situated one duct diameter downstream of the second elbow. Likewise, detailed velocity profile measurements were performed at various axial locations in the straight section between the close-coupled elbows. The specific locations of the velocity profile measurements are presented in Table 1 as a function of intermediate length x=D ð Þ, measured relative to the exit plane of the upstream elbow. The five-hole probe was first located at the centerline of the duct. Thereafter, the probe was re-positioned in radial increments of 25.4 mm (1 in.) on either side of the centerline; refer to Figure 2. The sensing tip of the five-hole probe was located with an accuracy of ±1.25 mm (0.05 in.), relative to the duct centerline. For every traverse performed, air velocities were measured in two mutually perpendicular planes. The two planes divided the cross section of the circular duct into four equal quadrants.

Computational domain and boundary conditions
The experimental data from  are used in this study for comparison purposes. The Z-shape duct with diameter, D = 12 inch (304.8 mm) is shown in Figure 3. The geometry consists of three parts, a 15D first horizontal section as the inlet section, a vertical separation distance section where the separation distance between the two elbows was varied systemically, namely L/D was investigated for five cases (2, 4, 6, 8 and 10) and a 10D second horizontal section as the outlet section. D is the diameter of the duct and L is the separation distance length. The inlet boundary condition is set to "velocity inlet" with a value such that the flow Reynolds number, Re ¼ 3:5 � 10 5 (same as experimental Re). The length of the entrance of the duct in the experimental study of  was elongated to 15D to ensure that the flow becomes fully developed before reaching the first elbow of the duct. For the present numerical study, the same 15D entrance section length was used to ensure that the flow becomes fully developed. The outlet condition is set as "pressure outlet" with zero-gauge pressure. All walls are assumed "non-slip" condition. A computational grid with a structured hexahedral cell type is constructed using commercial software ANSYS Workbench 16. Total mesh size is about 2 million grid nodes. As shown in Figure 4, the mesh is refined in the bending geometry by a further 30% compared to the bulk mesh region. The mesh at near wall region is also refined to keep y+ within the required range.

Numerical schemes
ANSYS FLUENT 16 solver is used for the CFD simulation. The spatial discretization is done using node-based gradient finite volume approach. The pressure equation is discretized with a standard scheme. Second-order bounded central differencing scheme is employed for the momentum equation. Temporal discretization used second-order accurate scheme. The convergence criteria are set such that the momentum equation residuals reduce to less than 10 −4 and scalars to less than 10 −6 . The LES solution is advancing with a time-step of 10 −5 sec. The eddy-viscosity, dynamic Smogorinsky-Lilly sub-grid scale model is employed. The simulations were performed on a supercomputer using 72 parallel processors. Higher number of iterations per time step (∆t) at the very beginning of the simulation was used, as many as 100 iterations for the first time step. The number of iterations per time step was gradually reduced from 100 to 15 as simulation progressed as it is required to converge the solution at every time-step due to coupling of the Navier-Stokes equations. NITA with PISO pressure-velocity coupling methods in order to speed up the simulations.

Grid refinement study
Mesh convergence analysis for LES simulations is not exactly the same as in RANS since the resolved variables are filtered. So, the convergence analysis has to be done carefully since the filter is changing for each mesh, for example, see (Weinman et al., 2006;You et al., 2010). LES requires mesh and time-step sizes sufficiently fine to resolve the energy-containing eddies. The mesh resolution determines the fraction of turbulent kinetic energy directly resolved. Turbulent kinetic energy peaks at l 0 integral length scale. This scale must be sufficiently resolved by having few cells in each direction to resolve an eddy with a length scale l (half length, ∆ = l/2). In the current study, the aim is to resolve 80% of the turbulent kinetic energy, so we need to resolve the eddies whose sizes are larger than roughly half the size of the integral length scale l 0 . To achieve this, approximately 5 cells are placed across the integral length scale l 0 . Then, contours showing the ratio l0/∆ were plotted at different cross sections. The upper values of l0/∆ were removed so that the well-resolved areas do not appear and the not-so-well resolved regions could be identified easily. Critical under-resolved regions near the elbows of the flow domain were remeshed.
In LES, unlike in RANS, both the wall-normal and wall-parallel (stream-wise & spanwise) spacing must also decrease to resolve smaller eddies, so the density of grid points should, ideally, increase in all three directions as a wall is approached. The LES wall functions are used in the current study and the near-wall mesh resolution was taken such that yp+ is about 10. The combination of the subgrid-scale (SGS) WALE model was used that give correct levels of SGS viscosity in the near-wall regions.

Mean flow velocity profiles
This section presents and discusses the validation of numerical predictions compared to experimental data from Sleiti et al. (2017). Figure 5-10 show the distributions of mean flow velocity profile at each selected cross-section in the fluid domain at different duct lateral separation distances, L/D. The mean flow velocity profiles in the Figures, V* is the ration of the local mean flow velocity normalized by velocity at the center of each section. The normalized radial distance, r* is defined as the ratio of radial distance from the center of the duct to the duct radius. Such normalized radial locations have a value of 0 at the center and 1 at the duct inner wall. Each profile section is compared using iso-line along both N-S (North-South) and E-W (East-West) orientations as shown in Figure 2.
Different duct lateral separation distance configurations, L/D are validated against experimental data, namely, L/D = 2, 4, 6, 8, and 10. The total number of fluid sectional planes (x/D) used for each configuration is varied depending on the different L/D configurations. For example, there are total of five sections used for the longest configuration for the case of L/D = 10 (i.e. X/D = 1, 3, 5, 7, 9), whereas velocity profile data are only available at a single section for the shortest configuration of L/D = 2. The locations of these sections (x/D) are given in Figure 3. The "inlet" section is referred to the fluid plane before the first 90 deg elbow turn. The "outlet" section is referred to the fluid plane after the second 90 deg elbow turn. The remaining x/D = 1 to 9 sections are located along the straight duct in between elbows. Detailed LES velocity profiles compared to experimental profiles at each section are presented and discussed below for different L/D configurations. This is done in an attempt to fully assess the LES capabilities in predicting such flow and to explain the discrepancies. After that, the LES velocity contours are plotted and analyzed followed by kinetic energy contours and discussion.

Velocity profiles at the Inlet section (before the first elbow)
The "Inlet" section represents the section right before the first elbow, see Figure 3. As shown in Figure 5, the flow has become fully developed, similar to the experiment flow profile for all L/D configurations. Compared to experiment data, the predicted LES mean flow velocity results at inlet section near the duct center are in excellent agreement (within 2% difference). The sharp ends in the computational velocity profile may indicate that the grid resolution close to the wall is not sufficient; however, as the experimental uncertainty of the velocities was about 4%, the difference is considered acceptable and within the experimental uncertainty. The near center bulk flow structures are well resolved with the current LES method to match a fully developed flow profile. Towards the near-wall region inside the sub-viscous layer, the flow velocities are also in excellent agreement (within 3%). However, the current LES flow velocities are slightly overpredicted in the narrow region between the flow boundary layer and bulk flow. Such a slight overshoot issue could possibly be due to a lack of accurate scheme implementation to capture small and medium-sized eddies structures with the current SGS model.

Velocity profiles at x/D = 1 section
At x/D = 1, after the flow makes the first turn, the flow velocity distribution becomes highly nonuniform due to the curvature, which causes local flow separation and re-attachment regions, Figure 6. The flow pattern has demonstrated high asymmetry 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. Similar patterns are found consistently for all L/D configurations. The local reduction of mean flow velocity occurs inside the flow separation region and local increase of flow velocity occurs inside the flow re-attachment region. The re-attachment region occurs when the local mean flow velocity ratio is found greater than 1.0.
It is noted that the profiles for all L/D configurations show a very similar pattern despite the differences in separation distances, L/D. The flow distribution after the first elbow is largely dominated by the elbow turning radius which in this case is the same for all L/D configurations; therefore, the velocity distributions remain unchanged with different lateral separation distance. In general, the predicted flow distribution trends are similar to the experimental data. The current LES results are slightly under-predicted the local separation regions consistently for all L/D configurations as shown by EW profile. In addition, the EW profiles showed that LES predicted flow reattachment region located near r* = 0.2-1, while experimental data are located near r* = 0.18-1. On the other hand, flow separation region as indicated by V*<1.0 is found in the region within r* <0.
NS profiles showed a similar distribution in the near-wall region as at the "inlet" section. However, flow velocities are drastically reduced near the center of the duct at this section due to loss of flow momentum, dominated by the flow separation region. Both experiment and LES results demonstrated similar trends with LES prediction slightly underperforming near flow separation region. It is also believed the reduction of velocities are attributed to flow swirling in the circumferential direction.

Velocity profiles at x/D = 3 and x/D = 5
When the flow convects along the straight duct section, the flow distribution becomes more uniform downstream of the duct. Figures 7 and 8   more symmetrical along EW direction. This is expected as the flow transitioning from non-uniform to uniform pattern further downstream at x/D = 5 location including upstream flow separation region. Flow structures are expected to be re-attached to the wall as the flow moves towards downstream of the duct. Similarly, LES also predicts the same local flow velocity deficit region near r* = −0.8 to 0 locations as shown in EW profile; however, the profiles are under-predicted compared to experiment. The current LES method is found to have some discrepancies in capturing the flow transitioning from separation to re-attachment. Numerical errors from upstream flow transitioning at x/D = 3 are also found extending towards downstream of the duct as indicated by the discrepancies appeared at x/D = 5. NS profiles show a distribution that is similar to the previous one at x/D = 1 with more flow recovery near the duct center. The peaks found between duct center and wall are attributed to the flow swirling effect in the circumferential direction. Both experiment and LES demonstrate similar trends. The LES results are closely following the experimental distribution along the NS direction.

Velocity profiles at x/D = 7 and x/D = 9: Downstream In-between Elbows
Further downstream in the straight duct, two sectional planes were studied and analyzed at x/ D = 7 and 9 (Figures 9 and 10). The furthest downstream x/D = 9 section is located right before the second elbow. It is found that the flow velocity in these two sections becomes much uniform than  upstream regions. Such uniform patterns are deemed to fully turbulent flow development profile as expected near downstream of the straight duct. Such patterns are also found consistently for all L/D cases. Both LES results and experiment are in good agreement for both EW and NS profiles. However, EW profiles showed some slight discrepancies and such over-predictions results are deemed to numerical errors accumulated from upstream flow structures. In general, NS profile shows that large-scale flow structures are sufficiently resolved by current LES methodology as the results showed acceptable agreement with experiment.

Velocity profile at the "Outlet" section: after second elbow
At the "outlet" section after the second 90 deg elbow, the mean flow velocity distributions demonstrate higher velocities at negative radial locations and lower velocity at positive radial locations, Figure 11. This is attributed to the presence of local flow separation region near r* = 0.2-0.9 location. Such distributions are consistently found in all L/D configurations. Overall, EW profile of both LES and experiment are in fairly good agreement for all L/D configurations. The longest lateral separation L/D = 10 produces the most satisfying results that match the experiment. In general, both LES predictions and experiment data demonstrate the local flow separation region near r* = 0.2-0.9 location with LES slightly less accurate in most cases. Such local flow separation region is expected as flow is turning at 90 deg around the elbow, it is expected the local flow swirling effect generates such local flow circulation. On the other hand, NS profile is also validated as shown in the figures below and both LES and experiment showed excellent agreement for all L/ D cases.

Mean flow velocity contours plots
LES data are known to include a large number of insightful flow variables, this section discusses, in further details, the two-dimensional contours of mean flow velocities at each corresponding fluid plane previously analyzed at x/D = inlet, 1, 3, 5, and outlet locations for different L/D configurations. Despite there are no experimental data available to fully validate these LES contours, it is essential to visualize the footprint of the entire fluid planes as only part of the planes were validated in the previous section. From Figure 12, the "Inlet" section contours for different L/D consistently showed a much uniform and higher velocity distribution near the center of the duct, except for L/D = 6 configuration, where slightly higher velocities are found near the lower region. However, the differences between red and orange colors scale represented by the legend is only about 1 m/s in velocity magnitude. Therefore, such velocity difference is considered small and only narrowly misses the target velocity.
Velocity contours after the first elbow at x/D = 1 location are found highly non-uniform with an asymmetrical pattern. Such distribution is deemed to the presence of local flow separation and reattachment near the opposite side of the duct walls. High velocities are found in flow re-attachment within the outboard of the duct, whereas low velocities are found in the flow separation region within the inboard of the duct. Both flow regions dominated almost half of the duct circular area on opposite sides. This pattern is found consistent for all L/D configurations.
Flow distribution near upstream of the straight duct is also compared quantitatively for all L/D configurations at x/D = 3 section. Local flow is found recovering from separation as indicated by yellow color scale at about 11.5 m/s near inboard of the straight duct despite such local velocity is still lower than inlet velocity. Flow is also found accelerating near the outboard wall of the straight duct. The iso-color distribution in the fluid plane distributed circumferentially also indicates the presence of flow swirling around the straight duct in the circumferential direction. These secondary flow structures present in the straight duct due to swirling effect.
The downstream x/D = 5 contours represent the mean flow distribution near further downstream of straight duct located before the second elbow. This sectional circular area is the same as the previous x/D = 3. As can be seen, the velocity becomes more uniform as the flow is transitioning from an upstream non-uniform pattern. Flow swirling effect is reduced. Such distributions are also found consistently for all L/D configurations. Last, the "outlet" contours located after the second elbow demonstrated that the velocity distribution is highly non-uniform in the NS direction. The low-velocity region is found in the local flow circulation region near the upper (North) region as shown in Figure 12. Such local flow re-circulations are found due to flow separation after the flow is turning 90 deg elbow.

Turbulent kinetic energy contours
The instantaneous flow velocity components u, v, and w can be decomposed in the following manner: where � u and � v and � w represent mean flow (timeaveraged) quantities, while u 0 , v 0 and w 0 are fluctuating components of velocity due to turbulence. The resolved turbulent kinetic energy per unit mass, k can be expressed with these fluctuating components.
The process of turbulent kinetic energy production, transport and dissipation can be expressed as the following: where Ñ � T 0 is the turbulence transport of k P is the production of k ε is the dissipation of k The phenomena of each process that governs the motion of the turbulent flow can be investigated using turbulent kinetic energy budget. It is expected that the production term to be greater than the dissipation term in the region where k is increasing near the elbow region. This implies the local boundary layer becomes more turbulent in that region. On the other hand, the boundary layer becomes less turbulent in the region where the production term is found less than the dissipation term.
In Figure 13, the LES prediction of k of cross-sectional planes at the centre of the duct are demonstrated for different L/D configurations. As can be seen, k consistently increases to the peak with maximum value >20 J/kg near the turning radius at both elbows for all cases. Along the vertical section, k reduces gradually until it reaches the second elbow and increases gradually once again near the corner transition. The sheared flow is developed with the presence of turbulent mixing enhancement. It is believed that when k is cascading down due energy transfer, large-scale flow structures are dissipated by viscous forces at the Kolmogorov scale. Compared to other longer lateral separation distance configuration, the shortest L/D = 2 design seems to have higher k in the region between elbows. This phenomenon is largely due to limited vertical length in streamwise direction to facilitate the viscous dissipation of turbulent structures.

λ2-criterion Iso-surfaces
The 3D turbulent flow structures' visualization is detected using λ 2 iso-surfaces post-processing technique as demonstrated in Figure 14. The λ 2 -criterion looks for a local pressure minimum that entails a vortex. The representation was done by removing both the effects from unsteady straining and viscosity terms. The close-up view near the first elbow in Figure 14 highlights the turbulent structures shed around the corner upstream location. The flow structures can be seen dissipating downstream of the elbow as such structures are convecting in the vertical duct. It is believed the viscous dissipation mechanism is dominant near downstream while upstream is dominated by turbulent production mechanism. Therefore, in the downstream region, the smaller eddies are seen due to dissipation. Also, the turbulence structures in the upstream region are found to have a larger scale in size than those in downstream. The larger size turbulent structures also dictate the swirling effect of the flow. Therefore, high vorticity (>100 1/s) flow structures are dominating in the elbow region near the wall.

Conclusions
The complex transient flow in Z-shape ducts was simulated using LES and compared to the experimental data for Re = 350,000. The flow and turbulence behavior as a function of the lateral separation distance (L/D) in the Z-shape duct are investigated in this research study for a wide range of L/D from 2 to 10. L/D = 2 cases represent a Z-shape duct where the two elbows are connected directly without any additional lateral duct. The following has been concluded: • LES predictions of the mean flow velocity profiles are in good agreement (within the experimental uncertainty) with experimental data for the investigated wide range of L/D configurations. Such good trend predictions are attributed to the well-resolved large-scale flow structures. However, there are some slight over-predictions and under-predictions when comparing absolute values in these flow regions. It is believed that such numerical errors are deemed to the limited modeling approach to predict small eddies structures with the current SGS model.
• The main key features of the flow after the first elbow are also identified as separation and reattachment regions.
• Some bigger discrepancies are identified for lateral separation distances at sections x/D = 3 and x/D = 5 inside the flow transition region from non-uniform to uniform distribution. Some of those discrepancies are believed to be inherited from the upstream flow numerical errors that arise in the non-uniform flow mixing regions. Although it has not been studied, the potential remedy includes applying finer mesh resolution and/or higher order spatial discretization to accurately resolve the local velocity gradients and the complex flow structures. It is expected the downstream duct mean flow distribution accuracy can be improved when upstream discrepancy is reduced.

Funding
The publication of this article was funded by the Qatar National Library.