Benjamin A Toms^{13*}  
^{1}Summer Undergraduate Research Fellow, Fire Research Division, Engineering Laboratory, National Institute of Standards and Technology, Gaithersburg, Maryland, 20899, USA  
^{2}Department of Civil Engineering and Environmental Sciences, University of Oklahoma, Norman, Oklahoma, 73019, USA  
^{3}Department of Atmospheric and Geographic Sciences, University of Oklahoma, Norman, Oklahoma, 73019,USA  
Corresponding Author :  Benjamin A Toms Department of Civil Engineering and Environmental Sciences University of Oklahoma, Norman, Oklahoma, 73019, USA Tel: 303 483 8268 Email: [email protected] 
Received May 22, 2015; Accepted June 19, 2015; Published June 26, 2015  
Citation: Toms BA (2015) Largeeddy Simulation of Flow Over a Backward Facing Step: Assessment of Inflow Boundary Conditions, Eddy Viscosity Models and Wall Functions. J Appl Mech Eng 4: 169. doi:10.4172/21689873.1000169  
Copyright: © 2015 Toms BA. This is an openaccess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.  
Related article at Pubmed, Scholar Google 
Visit for more related articles at Journal of Applied Mechanical Engineering
Keywords 
Largeeddy simulation; Backward facing step; Turbulent boundary conditions; Wall functions; Computational fluid dynamics 
Introduction 
The Fire Dynamics Simulator (FDS) is a largeeddy simulation (LES) code used to model lowmach flows driven by combustion heat release and buoyancy [1,2]. FDS has recently been applied to urban canopy modeling [3] and wind engineering. As such, it is important to study geometrically influenced flows, including the appropriate specification of boundary conditions and resultant downstream flow evolution. 
The backward facing step is a commonly used geometry in the study of turbulent flows.Extensive research on this flow has been conducted using direct numerical simulation (DNS), ReynoldsAveraged Navier Stokes (RANS), and LES. DNS has offered computational results most consistently resemblant of experimental data (e.g. the comparison of LES results of Panjwani et al. [4] to the DNS results of Le et al. [5] and experimental results of Jovic and Driver [6]. However, studies performed by authors such as Saric et al. [7] have suggested optimally conditioned LES code may simulate the flow as accurately as DNS. 
The backward facing step has been utilized to study such cases as wall heat transfer, pressure drops following kinematic interaction with complex geometry, and influence of geometry on turbulence related mixing of fluid properties [810]. Although the foundation of research relating to the backward facing step is meritorious, parametric study focusing on influence of commonly varied computational domain characteristics (e.g. wall boundary conditions, inlet turbulence, etc.) is limited. The number of studies using FDS is further limited. As such, the goal of the present work is to confirm the capability of FDS to model geometrically complex flows and to assess the importance of boundary conditions and subgrid models on the results. 
Le et al. [5] first modeled flow over a backward facing step using DNS with a Reynolds number of 5,100. To verify the results, Jovic and Driver (J and D) [6] supplemented the DNS simulation with an identically proportioned wind tunnel experiment. Together, the data sets from these two studies have provided the baseline for analysis of recent simulations of flow over a backward facing step. 
Panjwani et al. [4] used LES to study the influence of eddy viscosity subgridscale (SGS) modeling, discretization schemes, and grid refinement on LES accuracy in comparison to the DNS of Le et al. [5]. Effects of SGS modeling were determined to be minimal, whereas the discretization schemes were found to have a substantial impact. 
Aider et al. [11] studied the influence of inlet turbulent boundary conditions on downstream flow characteristics. Inlet turbulence boundary conditions were found to have significant impact on flow development. These results supplement those of Isomoto and Honami who concluded the maximum turbulence intensity to have an inverse relationship with reattachment length, and are further confirmed by Dol et al. who determined integration of a pulsation frequency on inlet kinematics influence downstream eddy development. Namely, Aider et al. [11] found higher values of rootmeansquare (RMS) velocity led to a more rapid destabilization of the shear layer that develops downstream of the step. This resulted in a reduction in reattachment length and a faster redevelopment of the mean stream wise velocity profile, . (Throughout this document, an over line is used to denote a mean quantity.) 
Sarwar et al. [12] used FDS to compare SGS eddy viscosity models. The constant Smagorinsky model [13] performed the best, although Germano [14,15] and Deardorff [16] models, are nearly equivalent in accuracy, were found to perform better than the Vreman model [17]. To avoid explicit specification of inlet turbulence conditions, Sarwar et al. [12] created an extremely long inlet section to allow turbulence to develop. One of the goals of the present study is to show that synthetic turbulent inflow boundary conditions can be successfully used to greatly reduce the size of the inlet domain and hence the calculation cost. 
For the present study, grid resolution, inlet turbulence definition, SGS eddy viscosity modeling, and wall models are analyzed using FDS 6.1.0 to determine the most influential variables on flow development. Simulations with varying grid resolutions are compared to the J and D experimental data. Using adequate grid resolution, the effects of inlet RMS velocity magnitude are studied. The effectiveness of the Deardorff and dynamic Smagorinsky SGS eddy viscosity models are compared. Influences of the wall condition are determined through a comparison of an LES loglawbased wall function with the noslip and freeslip wall conditions. 
Methodology 
Computational domain 
A schematic view of the computational domain is provided in Figure 1. The dimensions of the channel, based on step height h=.0098 m, are defined as follows: length of channel, L_{x}=24 h; length of inlet, L_{i}=4 h; width of channel, L_{y}=4 h; height of channel inlet, L_{iz}=5 h; height of channel poststep, L_{z}=6 h. The expansion ratio, defined as Liz/Lz, is thus 1.2. The inlet is split into three subinlets to permit localized variation of inlet turbulence. The simulation uses a stepheight Reynolds number of 5,100 based on a free stream velocity of U_{0}=7.2 m/s. 
To give the reader a qualitative sense of the flow field, Figure 2 presents instantaneous contours of velocity magnitude. For quantitative results, virtual measurement devices are placed throughout the channel to collect data relating to flow characteristics such as velocity, turbulence RMS velocity, and friction velocity. These virtual measurement devices are placed into lines—four vertical and one horizontalwith a device in the volumetric center of each grid cell. A vertical line device is placed within the inlet region at a location of x=−3h, and three line devices are placed in the poststep region at locations of 4h, 6h, and 10h. The poststep vertical line devices are placed accordingly to sample the recirculation, reattachment, and recovery regions. A horizontal array of virtual anemometers is used to sample velocity data directly adjacent to the bottom wall of the channel in the poststep region (0h to 20h). 
Boundary conditions 
The inlet profile is established using experimental data provided by J and D [6], while , , and inlet freestream turbulence are set to zero. Based on the profile provided by J and D, the prestep boundary layer depth is equivalent to h [6]. The summed vertical span of the two bottom inlets is equivalent to the inlet boundary layer height. Turbulent eddies are injected using the Synthetic Eddy Method of Jarrin [18]. Eddies are injected at random locations in the bottom two inlets, advected with the flow over a distance equal to the maximum eddy length scale l_{0}, and recycled at the inlet. The inlet maximum eddy length scale, number of eddies, and RMS velocity is predefined. The Reynolds stress and eddy length scale are assumed to be isotropic. 
The boundary conditions for velocity and pressure on the top of the domain are “mirror”, that is, zero normal gradients. The span wise boundaries are periodic. The outlet boundary is “open”, meaning the outlet pressure is specified based on the direction of the flow [1]. 
In an attempt to replicate the J and D experiment, the inlet boundary conditions are modified to resemble the experimental inlet data as closely as possible. This is complicated by inlet turbulence profile limitations of FDS. Currently, FDS allows specification of turbulence parameters unique to each inlet (a “vent” in FDS nomenclature), but does not allow direct designation of an inlet turbulence profile. The turbulence parameterization that results in kinematic profiles and turbulence statistics most representatives of J and D data is used as the baseline case throughout the study. The RMS velocity is 1.0 m/s for the bottom inlet (Inlet A), 0.5 m/s for the middle inlet (Inlet B), and 0.0 m/s for the top inlet (Inlet C). For both inlets A and B, the maximum eddy length scale is 0.03 m and number of injected eddies is 100. The number of eddies is set high enough to ensure correct inlet statistics. Figure 3 compares J and D data to the resultant inlet mean velocity and turbulence statistics profiles using a grid resolution of h/δz=10. 
Of note, J and D provided error associated with the experimental measurements – these values have been incorporated as error bars into subsequent plots. J and D designated the error associated with velocity measurements, Reynolds stress and turbulence intensity measurements, and friction coefficients as, respectively, 2%, 15%, and ± 0.0005. 
Results and Discussion 
Grid resolution 
Grid resolutions of h/δz= 5, 10, and 20 are compared to determine the minimum grid resolution requirement. Case 3(ii) is used as the baseline case for grid resolution analysis (Table 1). For validation, in the following sections the longitudinal friction coefficient and pressure coefficient profiles, mean velocity profiles, and turbulence intensity and Reynolds stress profiles are compared to J and D experimental data [6]. 
Friction coefficient and reattachment length: The reattachment length is defined as the location of zero vertical gradient of longitudinal velocity at the wall: . If two locations of occur due to a secondary recirculation region [19], the point farther from the step is designated the reattachment point. The reattachment lengths for the 5, 10, and 20 resolution cases are, respectively, x/h=9.9, 5.95, 5.275. Averaged flow quantities have been shown to be dependent on reattachment length [20]. Therefore, due to the significant difference between the h/δz=5 and J and D reattachment lengths, all statistics are poorly represented by the lowest resolution case. 
The friction coefficient is defined as 
(1) 
where 
(2) 
and is therefore proportional to the vertical gradient of at the wall. Figure 4a shows poststep longitudinal Cf profiles for each resolution. The h/δz=10 case correctly displays the qualitative characteristics of the longitudinal Cf profile, but does not realize the magnitude of the sharp dip within the recirculation region. In the recirculation region, the magnitude of C_{f} is most accurately captured by the h/δz=20 case, although the recirculation length is too short and thus the qualitative characteristics are skewed toward the step. The h/δz=5 case incorrectly models both the qualitative and quantitative characteristics of the C_{f} longitudinal profile [21]. 
There is a direct relationship between grid resolution and C_{f} magnitude within the recirculation region. This is likely the result of enhanced sampling of the nearwall region with increased resolution. For the given computational domain, DNS simulations have shown deviations from the LES loglaw viscous sub layer depth (z+ <12) [5], necessitating enhanced grid resolution to correctly resolve the nearwall profile. The h/δz= 20 case is the only resolution capable of sampling data within both the viscous sub layer and buffer layer in the mean. 
Pressure coefficient: The coefficient of pressure is defined as 
(3) 
Where p is the pressure and p_{0} is the background pressure. The pressure coefficient is thus proportional to localized pressure perturbations. Figure 4b presents the computational C_{p} data compared to J and D data. 
The h/δz=5 case performs poorly. The h/δz=10 case is qualitatively correct, but results in the pressure coefficient profile being shifted downstream relative to the J and D data. The h/δz=20 case is within 5% of the J and D data at the reattachment location. However, the recirculation length for the h/δz=20 case is shorter than that of J and D, so the quantitative accuracy in the streamwise extremity of the recirculation region is questionable. All grid resolutions result in similar maximum magnitude of pressure perturbation within the recirculation region. 
Wallnormal velocity and Reynolds stress profiles: Figure 5 shows wallnormal velocity and Reynolds stress profiles for the grid resolution study. The h/δz=5 case performs poorly for all kinematic variables and will not be discussed further. The profile of at h/ δz=20 matches J and D data more closely than the h/δz=10 case. The profiles for both cases identically match J and D data. 
With increasing resolution, the spike in near the wall is better captured. The values of the h/δz= 10 and 20 cases are within the margins of error of the J and D data, except for within the recovery region. The and (Reynolds stress) values are too large in both cases. 
Minimum grid resolution requirement: The h/δz=5 grid resolution performs poorly. The h/δz=10 grid resolution adequately models the qualitative characteristics of the flow. Although the h/δz=10 and h/δz=20 grid resolution simulation results differ quantitatively, both cases are within experimental error of data provided by J and D. Therefore, the h/δz=10 grid resolution is adequate both qualitatively and quantitatively, and is utilized throughout the remainder of the study. 
Inlet turbulence modification 
The significance of inlet turbulence is studied through a comparison of a set of inlet turbulence cases. Refer to Table 1 for parameterizations of each case. In the first case, no turbulence is injected. For the second case, turbulence of the same RMS velocity magnitude is injected into both lower inlets. For the third case, the RMS velocity magnitude of the middle inlet (inlet B) is half that of the bottom inlet (inlet A). These cases are henceforth referred to as the no turbulence, two inlet, and three inlet cases respectively. Neither the maximum eddy length scale nor the number of injected eddies is varied. 
The reattachment lengths for the no turbulence, two inlet, and three inlet turbulence modification cases are, respectively, x/h=11.25, 5.55, 5.95. This aligns with previous research suggesting inlet RMS velocity magnitude and reattachment length are inversely related [11]. The no turbulence case drastically overestimates the reattachment length, and therefore significantly differs from J and D data for all studied flow variables. 
The longitudinal C_{f} and C_{p} profiles for the two and three inlet turbulence cases are similar, although the recirculation region associated pressure gradient is larger for the three inlet case. Figure 6 shows and values for all three cases. As discussed in Sec. 3.1.3, the values for the three inlet case falls within experimental error bounds of the J and D data, except within the recovery region. The two inlet case exhibits excessive values within the recirculation and recovery regions. Both cases exhibit excessive and values throughout the channel, although the three inlet case is more accurate. 
Turbulence model 
In this section, the Deardorff and dynamic Smagorinsky eddy viscosity models [1] are compared. Of note, in the first layer of cells adjacent to the wall, FDS uses the constant Smagorinsky model with Van Driest damping [1]. This is to avoid inconsistencies related to test filtering operations near the wall. The results are shown in Figure 7. 
There is no distinguishable difference between the poststep C_{f} profiles, reattachment lengths, C_{p} profiles, or profiles. The reattachment length for both models is x/h=5.95. All subsequently discussed differences between the models are small. The dynamic Smagorinsky model shows more negative w within the recirculation region and within the region near the reattachment point. Additionally, the turbulence intensity within the recirculation region is larger for the dynamic Smagorinsky model than for the Deardorff model. In contrast, the Deardorff model exhibits larger turbulence intensity within the recirculation region. Neither model captures the spike in near the wall at the inlet. 
Aside from these small differences, the Deardorff and dynamic Smagorinsky eddy viscosity models result in similar flow characteristics. 
Wall model 
The LES loglaw wall model, noslip wall condition, and freeslip wall condition are also compared (Figure 8). For the LES loglaw wall model, the stream wise velocity component in the grid cell closest to the wall, u(δz/2), is sampled from the following piecewise function: 
(4) 
Where u+=u/u∗ and z+=z/δν. Here, u∗ is friction velocity and δνis the viscous lengthscale; κ= 0.41 is the von Kármán constant and B=5.2 for a smooth wall. 
For the noslip wall condition, wallnormal gradients of streamwise velocity are calculated using a finite difference. For the freeslip wall condition, the gradient of streamwise velocity in the wallnormal direction is zero at the wall and hence there is no stress at the wall. 
The reattachment points for the loglaw model, noslip condition, and freeslip condition are x/h=5.95, 6.05, 5.85 respectively. As shown in Figure 9a, both the LES loglaw model and noslip condition Cf profiles qualitatively match that of J and D. However, the LES loglaw model results are more accurate. Neither model captures the dip in Cf within the recirculation region. The freeslip Cf is zero throughout the channel due to . There is no discernible difference in Cp profiles between the LES loglaw model and nonslip condition, although the freeslip condition results in an excessive negative pressure perturbation within the recirculation region (plot not shown). 
As expected, the freeslip condition is unable to resolve the nearwall velocity gradient within the inlet and recovery regions. Based on the u(z) profile, the noslip condition most accurately models the mean gradient of stream wise velocity near the wall. There are no significant differences in the profiles. 
There are no significant differences in the or profiles. The profiles, shown in Figure 9b, experience the largest variations. The nearwall spike in in the inlet and recovery regions is not captured by any of the wall model cases. However, as previously discussed in 3.1.3, this is likely more dependent on grid resolution than the wall model. With all three models, is slightly over predicted near the surface in the recirculation and reattachment regions. 
For the given computational domain, the noslip condition and LES loglaw model perform similarly. In the mean, the virtual measurement device closest to the wall is always within the linear sub layer of the LES wall model (z+<12). Therefore, for both the loglaw model and noslip condition, is calculated by a simple finite difference. Of note for the h/δz=10 resolution used in this study, within the inlet region the z+ location of the virtual measurement device closest to the wall is always within the log layer. Additionally, within the recovery region, the device records occasional excursions into the log layer. The loglaw model should, therefore, still be used for this case. 
In addition to the placement of the virtual measurement device closest to the wall residing within the linear sub layer of the LES wall model, Le et al. concluded the vertical profile of velocity within the recirculation deviates from the loglaw profile defined by the LES wall model. This is supported by present computational data, as shown in Figure 8. u+ is observed to be independent of z+ within the recirculation region, leading to significant deviations from the LES wall model specified u+ profile. This leads to the noslip condition being of similar, or greater, accuracy than the LES wall model for the given computational domain. Refer to Le et al. [5] for further details. 
Conclusion 
Largeeddy simulation of flow over a backward facing step is performed to study the effects of grid resolution, inlet turbulence conditions, eddy viscosity models, and wall boundary conditions. Computational results are compared to experimental data of Jovic and Driver [6]. 
The key finding of the present work is that, aside from grid resolution, specification of the mean profile and turbulence intensity at the inlet boundary is the most influential factor governing quality of the downstream flow profiles. A grid resolution of h/δz=10 is adequate, although nearwall kinematic features and locations of turbulence RMS velocity extreme are more accurately represented with twice the resolution. Based on present FDS limitations, we conclude that with grid resolutions of h/δz ≥ 10 the reattachment length may be predicted to within 10% (Table 1). The Deardorff and dynamic Smagorinsky reddy viscosity models produce very similar results. The noslip wall condition and LES loglaw model result in similarly accurate results due to the virtual measurement device closest to the wall residing within the wall model viscous sub layer (z+ <12) in the mean in the poststep region. However, the LES loglaw model is suggested to be the optimal wall boundary condition due to occasional instantaneous sampling in the loglaw region. 
Of note, the h/δz= 20 case performed worse for reattachment length in comparison to experimental data. The authors hypothesize this may be a manifestation of turbulence inlet profile specification limitations of FDS 6.1. However, the current study aims to present the simulations as is to compare the effects of various parameterization modifications, without modifying the LES code. Thus, while the authors recognize there is developmental work to be done on the Fire Dynamics Simulator that is not the aim of the current study. 
For future work, inlet Reynolds stress specification should be studied to determine if this would be beneficial for matching experimental data near the inlet. Additionally, for the present study, specification of vertical profiles of inlet turbulence statistics is limited by the current FDS requirement that turbulence parameters are unique to each inlet. Modification of FDS to permit specification of inlet turbulence profiles should be considered. 
Acknowledgement 
This work was supported by the NIST Summer Undergraduate Research Fellowship (SURF) program. 
References 

Table 1 
Figure 1  Figure 2  Figure 3  Figure 4  Figure 5 
Figure 6  Figure 7  Figure 8  Figure 9 