UFR 1-07 Test Case

Front Page

Description

Test Case Studies

Evaluation

Best Practice Advice

References

Unsteady Near-Field Plumes

Underlying Flow Regime 1-07

Test Case Study

Brief Description of the Study Test Case

• A summary of the boundary conditions is shown in Figure 8.
• A gas mixture mainly composed of helium is discharged through a circular orifice into ambient air.
• The gas is composed of 96.4% helium, 1.7% acetone and 1.9% oxygen by volume.
• The molecular weight of the gas released is 5.45 g/mol ±2.7%.
• The mixture is discharged at a temperature of T_He = 11°C ±3°C and the air is at T_air = 13°C ±3°C.
• The circular plume source has diameter, D = 1 metre.
• The helium is discharged at a Reynolds-averaged velocity V₀ = 0.325 m/s ±1.3% and a Favre-averaged velocity of approximately 0.339 m/s.
• The flow through the orifice is laminar.
• The ambient pressure is 80.9 kPa ±0.4 kPa.
• The measurements include:
  1. Time-history of vertical velocity at a point 0.5 m from the centreline and 0.5 m above the inlet, used to estimate the puffing frequency
  2. Measurements on a vertical plane through the plume from the plume source to a distance of one orifice diameter of:
     • Reynolds-averaged and Favre-averaged mean axial and radial velocities
     • Reynolds-averaged and Favre-averaged shear stresses, normal stresses and turbulent kinetic energy
     • Favre-averaged helium concentrations
  3. Movies of helium concentration and velocities
  4. Profiles of the mean and RMS velocities, and mean and RMS helium concentrations at six measurement positions (0.1, 0.2, 0.3, 0.4, 0.5 and 0.6 m downstream of the plume source)

Item 1 is available in the O‘Hern et al. [4] paper, Items 2 and 3 can be obtained by contacting the authors of the study. and Item 4 is presented by Chung & Devaud [39].

Test Case Experiments

The experiments selected for this UFR are those undertaken by O‘Hern et al. [4] at the Fire Laboratory for Accreditation of Models by Experimentation (FLAME) facility at Sandia National Laboratories, Albuquerque, New Mexico, in the late 1990s/early 2000s. The aim of these experiments was to examine the characteristics of turbulent buoyant plumes and provide data that could be used to help validate LES models suitable for modelling fires.

Figure 8  Boundary conditions for the O‘Hern et al. [4] experiments.

The experimental arrangement is shown in Figures 8 and 9. The main chamber has dimensions 6.1 × 6.1 × 7.3 metres and converges to a square chimney outlet at the top with nominal dimensions of 2.4 m on each side. The plume source is located in the centre of the chamber 2.45 m off the floor. Air is directed through a series of diverters, screens and honeycombs to form an annular low-velocity inlet flow surrounding the helium plume. A relatively large plume source (diameter, D = 1 m) was chosen to ensure that the plume would be fully turbulent. This is surrounded by a 0.51 m wide sheet of steel which simulates the ground plane. Air is drawn into the helium plume passing over this sheet flowing radially inwards. The experiments were designed specifically to mimic an unconfined plume on an infinite ground plane with negligible wind effects. Extensive CFD simulations were performed to help design the facility and to ensure that any separation bubble formed by the vertical flow of air around the 0.51 m ground plane did not disturb the plume.

The helium flowed through a diffuser, a series of perforated plates and three layers of honeycomb before being released through the orifice. The honeycomb immediately upstream of the orifice suppressed turbulence and flow visualization suggested that the inflow conditions were laminar. A detailed study of the inlet flow characteristics also found that the inlet velocity profile was uniform to within 6% [64]. Within just a few centimetres downstream of the inlet, observations suggested that the plume had become fully-turbulent. To ensure that the flow had reached a quasi-steady state, the helium was released for a couple of minutes before recordings were taken. Particle Image Velocimetry (PIV) was conducted using around 11,500 images spanning 70 puff cycles while Planar Laser-Induced Fluorescence (PLIF) analyses were performed on approximately 2,300 images, covering 33 puffs. The experiments were repeated 10 times and the inlet velocity was on average 0.325 m/s ±1.3% [4]. The acetone and oxygen needed to be added into the helium released in order for laser fluorescence. As a consequence, the molecular weight of the mixture was 5.45 g/mol ±2.7% compared to the pure helium value of 4.00 g/mol.

Figure 9  Schematic of the Sandia FLAME facility showing the laser-light sheet that bisects the plume and two video cameras that record the PIV and PLIF images.

The Reynolds number based on the inlet diameter and velocity, and the helium mixture properties was $Re=DV_{0}/\nu =3200\pm 0.6$% and the Richardson number was $Ri=\left(\rho _{\infty }-\rho _{0}\right)gD/\left(\rho _{\infty }V_{0}^{2}\right)=76\pm 6.5$%, where $\rho _{\infty }$ is the air density and $\left.\rho _{0}\right.$ the plume fluid density.

The PIV and PLIF measurements produced simultaneous time-resolved velocity and mass fraction data. The data was used to calculate density-weighted Favre-averaged statistics in addition to the more usual Reynolds or time-averaged statistics. Interestingly, the difference between the Favre- and Reynolds-averaged quantities was found to be less than the uncertainty in the data throughout the flow field [4].

The puffing frequency of the plume was analysed from the time-history of the vertical velocity at a point in space 0.5 m above the inlet and 0.5 m radially from the centreline. The recorded mean measured frequency was 1.37 Hz which compares well with the empirical correlation of $f=\left.0.8V_{0}{Ri}^{0.38}/D\right.$ from Cetegen & Kaspar [18] for helium-air plumes with Ri < 100, which gives a frequency of 1.35 Hz, and the empirical correlation of $f=\left.1.5D^{-1/2}\right.$ from Cetegen & Ahmed [25] for fire plumes which gives a frequency of 1.5 Hz.

O‘Hern et al. [4] discussed in some detail the dynamics of the unsteady plume and the role of the Rayleigh-Taylor instability in producing bubble and spike flow structures. Figure 10, taken from their paper, shows four snapshots of the plume where the spike and bubble structures are identified with arrows and the location of the large coherent puffing vortex is indicated with a circle.

Figure 10 Four snapshots of the helium plume of O‘Hern et al. [4] taken 115 ms apart. The left-hand-side of each image shows the mass fraction field from the PLIF, the right-hand-side shows the corresponding PIV vector field overlaid with the general plume outline. The development and movement of a large toroidal vortex is indicated by circular arrows. The spike and bubble structures characteristic of Rayleigh-Taylor instability are indicated by straight arrows.

Details of the uncertainties in the experiments are discussed at length in their paper. These include measurement errors due to the effects of out-of-plane motion and improper choice of peak correlation in the cross-correlation analysis of the PIV measurements, and the influence of film response, image registration and laser-sheet intensity normalization in the PLIF measurements. Overall, the uncertainties are estimated to be ±18% for the difference between the plume and air density, ±5% for the air density, ±20% for the velocities and ±30% for the turbulence statistics [2].

CFD Methods

DesJardin et al. [1]: Description of CFD Work

Governing Equations

Desjardin et al. [1] used the fully-compressible form of the Favre-averaged Navier Stokes equations. Transport equations were solved for the Favre-averaged momentum, species mass fraction and energy:

${\frac {\partial \left(\rho U_{i}\right)}{\partial t}}+{\frac {\partial \left(\rho U_{i}U_{j}+p\delta _{\mathit {ij}}\right)}{\partial x_{j}}}={\frac {\partial }{\partial x_{j}}}\left[\mu \left({\frac {\partial U_{i}}{\partial x_{j}}}+{\frac {\partial U_{j}}{\partial x_{i}}}\right)-{\frac {2}{3}}\mu {\frac {\partial U_{k}}{\partial x_{k}}}\delta _{\mathit {ij}}\right]+{\frac {\partial \tau _{u_{i}u_{j}}}{\partial x_{j}}}+\rho g_{i}$

$\left(12\right)$

${\frac {\partial \left(\rho Y\right)}{\partial t}}+{\frac {\partial \left(\rho U_{j}Y\right)}{\partial x_{j}}}={\frac {\partial }{\partial x_{j}}}\left({\frac {\mu }{\mathit {Sc}}}{\frac {\partial Y}{\partial x_{j}}}\right)+{\frac {\partial \tau _{u_{j}Y}}{\partial x_{j}}}$

$\left(13\right)$

${\frac {\partial \left(\rho e\right)}{\partial t}}+{\frac {\partial \left(\rho U_{j}h\right)}{\partial x_{j}}}={\frac {\partial }{\partial x_{j}}}\left[{\frac {\mu c_{p}}{\mathit {Pr}}}\left(1-{\frac {\mathit {Pr}}{\mathit {Sc}}}\right){\frac {\partial T}{\partial x_{j}}}+{\frac {\mu }{\mathit {Sc}}}{\frac {\partial h}{\partial x_{j}}}\right]+{\frac {\partial \tau _{u_{j}h}}{\partial x_{j}}}+{\frac {\partial \tau _{u_{j}u_{k}u_{k}}}{\partial x_{j}}}+{\frac {\partial \left(U_{k}\tau _{u_{j}u_{k}}\right)}{\partial x_{j}}}+\rho U_{k}g_{k}$

$\left(14\right)$

where ρ is the density, U_i the velocity components, p the pressure, Y the species mass fraction, e the total energy and h the enthalpy (N.B. all of these parameters are Favre-averaged quantities). For the diffusion of helium into air, the molecular Schmidt and Prandtl numbers were set to values of Sc = 0.2 and Pr = 0.7. Thermodynamic properties, such as c_p, were evaluated based on the mixture composition using the Chemkin libraries. The molecular viscosity, μ, was determined from Sutherland's law for pure air.

Turbulence Modelling

The Smagorinsky model was used for the SGS stresses in the momentum equation, $\tau _{u_{i}u_{j}}$ , and a simple Boussinesq gradient-diffusion model was used for the SGS stress terms in the species mass fraction and energy equations, $\tau _{u_{j}Y}$ and $\tau _{u_{j}h}$:

$\tau _{u_{i}u_{j}}=2\rho C_{R}\Delta ^{2}\left|S\right|\left(S_{\mathit {ij}}-{\frac {1}{3}}S_{\mathit {kk}}\delta _{\mathit {ij}}\right)$

$\left(15\right)$

$\tau _{u_{j}Y}=2\rho C_{Y}\Delta ^{2}\left|S\right|\left({\frac {\partial Y}{\partial x_{j}}}\right)$

$\left(16\right)$

$\tau _{u_{j}h}=2\rho C_{h}\Delta ^{2}\left|S\right|\left({\frac {\partial h}{\partial x_{j}}}\right)$

$\left(17\right)$

where $S_{\mathit {ij}}=1/2\left(\partial U_{i}/\partial x_{j}+\partial U_{j}/\partial x_{i}\right)$ is the strain rate and $\left|S\right|=\left(2S_{\mathit {ij}}S_{\mathit {ij}}\right)^{1/2}$ is the strain-invariant. The filter width was taken as twice the cube-root of the local computational cell volume, $\Delta =2(\delta x\delta y\delta z)^{1/3}$ , where $\delta x$, $\delta y$, and $\delta z$, are the cell widths. This is twice as large as the commonly used value and was chosen in an effort to minimize numerical errors. The three modelling “constants”, C_R, C_Y and C_h were calculated dynamically [65][66] using a second explicit filter that was twice the size of the first implicit filter. To ensure numerical stability, the constants were locally smoothed using explicit filtering with a filter function described by Fureby [67]. The magnitude of the effective viscosity or diffusivity was also clipped to be always greater than or equal to zero ( $\mu +\mu _{t}\geqslant 0$ and $D+D_{t}\geqslant 0$ ). DesJardin et al. [1] noted that the advantage of this approach over other approaches is that it allows for some backscatter that may be important for laminar to turbulent transition. Triple correlations appearing in the energy equation, $\tau _{u_{j}u_{k}u_{k}}$ , were modelled using an approach proposed by Ragab et al. [68][69], for details see [1]. DesJardin et al. [1] also presented results obtained without using the SGS model (i.e. a “no-model” approach).

Numerical Methods

DesJardin et al. [1] used a finite-volume treatment where the mass, momentum and energy equations equations were differenced in time using a fourth-order Runge-Kutta scheme. Convective terms were discretized using a blend of a fifth-order ENO scheme for the first two stages of the Runge-Kutta integration and ninth-order upwind-biased scheme for the final two stages. This combination of high-order schemes was chosen to prevent dispersive errors (undershoots and overshoots), minimize numerical dissipation and avoid odd-even decoupling errors (chequerboarding) in regions of the flow where the Mach number was small. Where the flow was aligned to the grid, their approach should have provided up to ninth-order accuracy for the momentum equations and fifth-order accuracy for mass, energy and species mass fraction. Diffusion terms were discretized using fourth-order central differences. Near the boundaries, the differencing schemes used for the convection and diffusion terms were of lower order accuracy.

To avoid having to use a time-step limited by the acoustic wave speed, DesJardin et al. [1] used the Pressure Gradient Scaling (PGS) method of O‘Rouke et al. [70][71]. This approach decomposes the pressure into two parts comprising thermodynamic and hydrostatic components. The thermodynamic component, which contains acoustic information, is pre-multiplied by a scaling factor which artificially reduces the acoustic wave speeds. Details of this technique are given in the Appendix of DesJardin et al.‘s paper [1].

Boundary Conditions

The vertical boundaries and top outlet planes were assumed to be open, allowing for flow to be entrained into or exit the domain. On the inlet plane, the inlet velocity for the helium was U_p = 0.351 m/s, whereas the experimental Favre-averaged velocity was 0.339 m/s [4]. A small axial coflow velocity of 0.01 m/s was specified outside the plume whereas in the experiments there was a fixed ground plane. The cross-stream velocities were set to zero on the inlet plane. A considerable amount of detail on the treatment used to avoid acoustic waves reflecting back from open boundaries into the domain and contaminating the solution was provided in an Appendix to their paper [1]. The non-reflective pressure relation used at open boundaries was based on an approach developed by Rudy & Strikwerda [72][73]. In the simulations, the gas released was pure helium with a molecular weight of 4.0 g/mol, whilst in the experiments, the gas released had a molecular weight of 5.4 g/mol.

The inlet velocity, the co-flow and the gas density differed slightly compared to the experiments since the simulations and the experiments were undertaken concurrently, and the final measured conditions differed from those originally planned.

It was not found necessary to superimpose turbulent fluctuations on the inlet velocity to obtain transition to turbulence. Tests found that using different prescribed inlet turbulence intensities did not affect the resulting flow behaviour.

Grid Used

Rather than model the same geometry as used in the experiments, it was assumed that the experimental conditions represented an unconfined plume. The computational domain was then constructed as a cube with sides of length 4 metres. DesJardin et al. [1] commented that this domain size was found necessary for the plume to be unaffected by the presence of the domain boundaries, due to the large quantity of air drawn into the plume in each puffing cycle. Note that in comparison, the central chamber of the FLAME facility, where the experiments were conducted, comprised a cube of sides 6.1 m and the ground plane around the plume source extended 0.51 m radially outwards from the perimeter of the source orifice (see Figure 9). Two grids were used comprising 80 × 80 × 80 (= 512k nodes) and 136 × 136 × 136 (${\approx }$ 2.5M nodes). The grids were refined near the centreline and the base of the plume resulting in minimum and maximum grid spacings of 2.8 cm and 13.1 cm for the coarse grid and 1.6 cm and 7.8 cm for the fine grid. A plot showing a cross-section through the mesh is given in their paper. It is not clear how the circular inlet orifice was modelled using the structured mesh although from their plot of the grid it would appear that a stair-stepped or sawtoothed approach was probably taken.

Time-Averaging

Calculations were run for 10 seconds of physical time to allow for initial transients to move downstream and for the flow to develop. Another 10 seconds of physical time were then used to collect at least 2000 realizations of the flow field for constructing mean and RMS values.

Calculations were performed using 128 processors with typical run times of 5.5 hours/processor for every second of physical time (a total CPU time of ~14,000 hours).

Discussion

The CFD methodology employed by DesJardin et al. [1] appears to have been performed to a high standard. Details of the modelling and numerical techniques used in their work were recorded clearly in their paper.

A number of studies have shown that upwind-biased numerical schemes should be avoided when using LES, since they can lead to excessive numerical dissipation [74][75]. DesJardin et al. [1] recognized that the use of upwind-biased convection schemes could introduce some undesirable numerical errors. They noted that whilst pure central-differencing schemes are commonly used in non-reacting flows, some degree of upwinding is necessary to stabilise the solution when there are strong scalar gradients. They chose a high-order upwind scheme to minimize unwanted dissipation and noted that they did not expect the flow predictions to be sensitive to the details of the discretization scheme.

The difference of more than 20% in the density of the gas released in their model compared to that released in the experiments is unfortunate and may have affected their results. Similarly, small errors may have been introduced by using a 0.01 m/s co-flow, and an inlet velocity of 0.351 m/s instead of the experimental Favre-averaged value of 0.339 m/s [4].

Tieszen et al. [2]: Description of CFD Work

Governing Equations

Tieszen et al. [2] used a low-Mach-number code developed at Stanford University by Pierce [76][77]. In the low-Mach-number limit, acoustic interactions, compressibility effects and viscous dissipation effects are neglected. The pressure is decomposed into a background pressure, p₀, and a flow-induced perturbation pressure, δp. The background component is assumed to be constant in space and time and is used in the energy equation and the state equation (the ideal gas law). For details of the equations solved, see [76].

Turbulence Modelling

The basic turbulence model used was the same as that employed by DesJardin et al. [1] (see above). This comprised the Smagorinsky LES model with coefficients determined using the dynamic procedure of Germano et al. [65] and the least-squares approach of Lilly [66]. Details of the averaging procedure used to smooth the dynamic constants are not provided in the paper by Tieszen et al. [2]. The test filter used in the dynamic procedure was twice the width of the grid filter. Subgrid-scale effects in the species mass fraction and energy equations were modelled using a gradient-diffusion approach.

Numerical Methods

The cylindrical form of the governing equations were solved and a structured mesh was used. Velocities were stored at staggered locations with respect to density and other scalars in both space and time. An energy-conservative, second-order central differencing scheme was used for convection in the momentum equations and an upwind-biased QUICK differencing scheme was used for the scalar equations. An iterative semi-implicit approach was used in time similar to the Crank-Nicolson scheme.

Boundary Conditions

For the boundary conditions, Tieszen et al. [2] stated simply that open boundaries were used on all domain surfaces except the floor and inlet. They noted that the inlet treatment differed slightly to that of DesJardin et al. [1] in that the flow in the diffuser is not specified directly using a Dirichlet condition but allowed to develop, ignoring the presence of the honeycomb at the inlet. Although this is not clear, it suggests that rather than impose a flat top-hat profile at the inlet, the flow was allowed to develop for some distance upstream of the orifice before discharging into the main flow domain.

Grid Used

Tieszen et al. [2] did not describe the overall size of the cylindrical domain used in their study, although from one of their figures (copied in Figure 14, below) it appears that it extended at least 1.6 metres radially and 2.7 metres axially. A structured mesh was used and grid-sensitivity studies were undertaken using three different mesh densities: a coarse mesh comprising in total around 250k nodes with 52 × 64 × 80 nodes in the radial × azimuthal × axial directions, a medium mesh of approximately 1M nodes (104 × 64 × 160), and a fine mesh with 4M nodes (208 × 64 × 320). In each case, the same number of nodes was used in the azimuthal direction. No details are given regarding any clustering of nodes and a plot of the mesh was not provided in their paper.

Time-Averaging

Time-averaging was performed once the flow had established a quasi-steady puffing mode. No further details were provided regarding the time-period over which averaging was performed.

Discussion

The description of their work in Tieszen et al.‘s paper [2] was not comprehensive and some important details, such as the size of the flow domain, the smoothing of the dynamic Smagorinsky constant, the time-averaging and the nature of the boundary conditions were not provided. Since the principal author, Tieszen, was a co-author of both the computational work of DesJardin et al. [1] and the experiments of O‘Hern et al. [4], it is reasonable to assume that appropriate choices were made for these aspects of the modelling and that their description was omitted simply in order to keep the paper concise. Tieszen et al.‘s study provides useful information on grid-dependence issues which are critical to understand for industrial LES. The main objective of the paper was to outline plans for a buoyancy-modified SGS model, which is discussed later.

Xin [3]: Discription of CFD Work

Numerical Methods

Xin [3] used the open-source LES-based CFD code: Fire Dynamics Simulator (FDS) version 3 which is developed and maintained by the U.S. National Institute for Standards & Technology (NIST). The code solves the low-Mach-number equations on a staggered Cartesian grid. The spatial discretization is second-order accurate and an explicit second-order predictor-corrector method is used in time.

Turbulence Modelling

FDS uses the Smagorinsky SGS model with constants of C_s = 0.2, Pr = 0.7 and Sc = 1.0. The model in this version of FDS was implemented in a rather unusual way. The diffusion term in the momentum equation is usually written:

${\frac {\partial }{\partial x_{j}}}\left[\left(\mu +\mu _{t}\right)S_{\mathit {ij}}\right]$

$\left(18\right)$

where μ is the molecular viscosity, μ_t is the SGS viscosity and S_ij is the strain-rate. In this version of FDS, the model was implemented as follows:

${\frac {\partial }{\partial x_{j}}}\left[{\mathit {max}}\left(\mu ,\mu _{t}\right)S_{\mathit {ij}}\right]$

$\left(19\right)$

i.e. the effective viscosity was taken as either the molecular viscosity or the SGS viscosity, whichever was largest.

Boundary Conditions

All of the boundaries were treated as openings except for the floor. At these openings the total specific pressure was set to zero if the flow was entering the domain or to u²/2 if the flow was leaving (see [78] for details). Air entering the domain was assigned ambient conditions. For the plume source, Xin [3] used an inlet flow velocity of 0.351 m/s which was superimposed with 1% random noise. This is the same mean velocity as that used by DesJardin et al. [1], despite the the experimental Reynolds and Favre-averaged velocities being 0.325 m/s and 0.339 m/s [4].

Grid Used

A rectangular computational domain was used with dimensions 2 × 2 × 6 metres. This corresponds to half the width and one-and-a-half times the height of the domain used by DesJardin et al. [1]. Since FDS uses a Cartesian grid, the circular inlet was modelled as a stair-stepped or sawtoothed shape. To help reduce unphysical vorticity being produced on surfaces with stepped boundaries, the ‘sawtooth’ model can be activated in FDS [78]. It is not clear from the description given in Xin's paper [3] whether or not this model was used.

The Cartesian grid was composed of grid with cells of sides 2.5 cm. This gives a mesh of 80 × 80 × 240 nodes with in total approximately 1.5M nodes.

Time-Averaging

Xin [3] stated that calculations were run for 20 seconds during which time there were more than 12 puffing cycles. During each simulation 1000 data samples were taken which were used for time-averaging. It is unclear from the description whether a period of time was allowed for the flow to develop before sampling took place.

Discussion

Details of the modelling and numerical techniques are only briefly described by Xin [3]. This reflects the fact that it was only a conference paper and not a peer-reviewed journal article.

The computational domain used was relatively small, only half the width of that used in the earlier study by DesJardin et al. [1]. In their paper, DesJardin et al. noted that a relatively large computational domain was necessary for the plume development to be free from the influences of the domain boundaries. It is possible therefore that some boundary effects could have influenced Xin's results. Given the strong dependence of LES flow predictions on the grid resolution, the choice of a small domain may have been driven by the desire to maximise the grid refinement at the possible expense of introducing some boundary effects.

There have been various different values proposed for C_s: for example, a value of 0.17 for homogeneous isotropic turbulence [56], 0.1 for mixing layers [79], and 0.065 for channel flows [80]. The constant has to be reduced to obtain the correct asymptotic flow behaviour near walls and the optimum value of C_s is also a function of the computational grid resolution [81]. It is reasonable to surmise that the Smagorinsky constant is not really a constant at all but a parameter dependent on the flow and the grid resolution. In Xin's simulations, the Smagorinsky constant was taken as the default FDS value of 0.2 throughout the whole flow domain. This may have lead to excessive damping of the turbulent structures in some regions of the flow. The sensitivity of LES predictions to the Smagorinsky constant was examined by Chung & Devaud [39], who also used the FDS code. They found that values of C_s between 0.15 and 0.20 produced the best overall agreement with the experiments of O'Hern et al. [4].

Footnotes

Front Page

Description

Test Case Studies

Evaluation

Best Practice Advice

References

Contributed by: Simon Gant — UK Health & Safety Laboratory

© copyright ERCOFTAC 2010