Axisymmetric buoyant far-field plume in a quiescent unstratified environment

Underlying Flow Regime 1-06

Test Case

Brief Description of the Study Test Case

The experiments used in this UFR are those of George et al. [3] which were conducted in 1974 at the Factory Mutual Research Corporation and were subsequently repeated by Shabbir & George [34] at the University of Buffalo.

Heated air is discharged through a circular orifice into ambient air that is at rest.
The plume source temperature is 300°C and the ambient air is 29°C.
The source has diameter, D = 6.35 cm.
The hot air is discharged at a velocity of U0 = 67 cm/s with a approximately a top-hat profile.
Temperature and velocity fluctuations at the inlet are less than 0.1%.
George et al. [3] present experimentally measured profiles

of both mean and fluctuating components of the temperature and axial velocity in the self-similar region at x/D = 8, 12 and 16 above the source.

Test Case Experiments

The experiments used in this UFR are those of George et al. [3] which were conducted in 1974 at the Factory Mutual Research Corporation and were subsequently repeated by Shabbir & George [34] at the University of Buffalo.

The general arrangement is shown in Figure 4. Compressed air is passed through a set of heaters and porous mesh screens before exiting through a nozzle into the enclosure. The nozzle is stated as a 15:1 contraction in [3], a 12:1 contraction in [sg92] and appears to be different again in a drawing of the arrangement in [3] (see Figure 5). It resulted in a velocity profile through the exit which was uniform to within 2% outside the wall boundary layer. The velocity and temperature fluctuations at the exit were measured to be very low, less than 0.1% in [3] and 0.5% in [34]. The temperature of the source was 300°C and the ambient environment 29°C. Both were controlled to an accuracy of within 1°C. The discharge velocity was 67 cm/s, as calculated from the measured heat flux. These source conditions corresponded to Reynolds number, Re₀ = 870, and densimetric Froude number, Fr₀ = 1.23 [#sdfootnote1sym ¹]. There was no evidence of laminar flow behaviour at a position two inlet diameters downstream from the source. The effective origin of the plume, x₀, was found to be at the same location as the exit (see [3] for details of how this was determined).

The screen enclosure around the plume exit was 2.44 × 2.44 metres in cross-section and 2.44 metres high (there is, presumably, an error in [3] which suggests that the enclosure is 2.44 × 2.44 × 2.44 mm). In the later Shabbir & George experiments, a 2 × 2 × 5 metre enclosure was used. The purpose of the screens was to minimize the effect of cross-draughts and other disturbances affecting the flow. Two-wire probes were used by George et al. [3] to record velocities and temperature.

[#sdfootnote1anc 1]The densimetric Froude number is calculated here from the source and ambient temperatures, the exit velocity and source diameter given by George et al. [gat77], using Equation (1). However, George et al. [gat77] stated that the densimetric Froude number was 1.4. It is unclear how they determined this value. Using the approach taken by Chen & Rodi [cr80] in which the source density instead of the ambient density is used to make the density difference dimensionless, and Froude number is defined using the square of the expression given in Equation (1), this gives a Froude number of 0.80.

Figure 5 Schematic of the George et al. [3] experiments, from Shabir & George [11]

Figure 6 Schematic of the plume generator used in the experiments, from George et al. [3]

George et al. [3] reported that measurement errors, stemming from directional ambiguity of the hot wire and its thermal inertia, were around 3% for the velocity and lower for other mean and RMS values. The frequency response of the hot wires was estimated to be around 300 Hz compared to the frequency of the energy-containing eddies at around 50 Hz and the Kolmogorov microscale at 1 Khz. It was noted that measurement errors were likely to be higher on the outer edge of the plume where the velocity fluctuations were higher.

In their review of plume experiments, Chen & Rodi [1] noted that the data from George et al. differed significantly from earlier measurements by Rouse et al. [64]. However, they considered it to be more reliable due to its use of more sophisticated instrumentation. George [40], describes an experimental program at the University of Buffalo that was set up following publication of the original George et al. [3] paper to investigate possible causes of differences in experimental plume results. Possible sources of errors discussed included:

ambient thermal stratification
the size of the enclosure
the use of porous screens used to minimise disturbances from the far-field affected the plume source.
hot wire measurement errors

The most significant concern was ambient thermal stratification. One of the features of buoyant plumes in neutral environments is that the integral of the buoyancy across the whole cross-section of the plume, F, should remain constant and equal to the buoyancy added at the source, F₀. George [40] discussed how thermal stratification involving small temperature differences of the order of 1°C across a 3 metre vertical span would be sufficient to cause F to decrease to only 50% of the source value. This would be likely to cause differences in measured temperature and velocity plume profiles.

In the initial experiments of George et al. [3], the thermal stratification was not strictly controlled. However, results from later experiments published in the PhD thesis of Shabbir [32] (reproduced in [34] and [40]), which conserved buoyancy to within 10%, are in good agreement with the earlier results from George et al. [3]. This suggests that, perhaps fortunately, ambient thermal stratification did not contaminate the George et al. [3] results significantly.

A summary of the original results from George et al. [3] and those reproduced later by Shabbir & George [11] is presented in Table 3. Also shown are the recommended values from Chen & Rodi's review [1] and other studies. The parameters given in Table 3 relate to the following empirical formulae for the mean vertical velocity:

W=F_{0}^{-{1/3}}z^{-{1/3}}f_{w}

and effective buoyancy acceleration:

g{\frac {\Delta \rho }{\rho }}=F_{0}^{2/3}z^{-{5/3}}f_{T}

where $f_{w}$ and $f_{T}$ are Gaussian functions:

f_{w}=A_{w}\exp \left[-B_{w}\left({\frac {r}{x}}\right)^{2}\right]\qquad f_{T}=A_{T}\exp \left[-B_{T}\left({\frac {r}{x}}\right)^{2}\right]

The parameters, $l_{\Delta T/2}$ and $l_{w/2}$ are the dimensionless half-widths of the plume, as defined by the location where the normalized buoyancy or mean velocity falls to half its centreline value. The RMS temperature and axial velocity fluctuations normalized by their centreline mean values are denoted, $\left({\overline {t^{2}}}\right)^{1/2}/\Delta T_{c}$ and $\left({\overline {w^{2}}}\right)^{1/2}/W_{c}$ , respectively.

As noted earlier, Dai&nbspet al. [10][37] [38][39][41] disputed the accuracy of the George et al. [3] experiments and suggested that they had made measurements too near the source, before the plume had reached a fully-developed state. Their arguments are disregarded by Shabbir & George [11] [34].

Table 3

CFD Methods

Van Maele & Merci: Description of CFD Work

Numerical Methods

Van Maele & Merci [2] used the finite-volume-based commercial CFD code, Fluent, to simulate the plume experiments of George et al. [3]. For the discretization of the convective terms in the momentum, turbulence and energy equations a second-order upwind scheme was used. Diffusion terms were discretized using second-order central differences and the SIMPLE algorithm was used for pressure-velocity coupling. The flow was treated as axisymmetric and elliptic calculations were performed used a Cartesian grid arrangement.

The low-Mach-number form of the Favre-averaged Navier-Stokes equations were used. In this weakly-compressible approach, the density is treated as only a function of temperature and not pressure. Pressure only affects the flow field through the pressure-gradient term in the momentum equations. The ideal gas law is used to link the mean density, ${\overline {\rho }}$ , to mean temperature, T as follows:

p_{*}={\overline {\rho }}RT

where p_* is taken as constant and equal to the atmospheric pressure. The low-Mach-number approximation implies that the effect of the mean kinetic energy and the work done by viscous stresses and pressure are negligible in the energy equation.

Turbulence Modelling

Two turbulence models were used by Van Maele & Merci [2]: the standard k – ε model of Jones & Launder [65] and the realizable k – ε model of Shih et al. [66]. In the former model, the eddy viscosity is given by:

\mu ={\overline {\rho }}c_{\mu }{\frac {k^{2}}{\varepsilon }}

where c_μ is a constant equal to 0.09 and the standard k and ε equations are written:

{\frac {\partial }{\partial x_{j}}}\left({\overline {\rho }}U_{j}k\right)={\frac {\partial }{\partial x_{j}}}\left[\left(\mu +{\frac {\mu }{\sigma _{k}}}\right){\frac {\partial k}{\partial x_{j}}}\right]+P_{k}+G-{\overline {\rho }}\varepsilon \ \ \ \ (14)

{\frac {\partial }{\partial x_{j}}}\left({\overline {\rho }}U_{j}\varepsilon \right)={\frac {\partial }{\partial x_{j}}}\left[\left(\mu +{\frac {\mu }{\sigma _{\varepsilon }}}\right){\frac {\partial \varepsilon }{\partial x_{j}}}\right]+c_{\varepsilon 1}P_{k}{\frac {\varepsilon }{k}}-c_{\varepsilon 2}{\overline {\rho }}{\frac {\varepsilon ^{2}}{k}}+S_{\varepsilon B}\ \ \ \ (15)

where c_ε1 = 1.44, c_ε2 = 1.92, σ_k = 1.0, σ_ε = 1.3 and P_k is the production term due to mean shear. The terms G and S_εB are source terms related to the influence of buoyancy on the k and ε equations. The treatment of these gterms is discussed below.

The Shih et al. [66] model involves two changes to the standard k – ε model. Firstly, c_μ is made a function of strain and vorticity invariants to ensure that the model always returns positive normal Reynolds stresses and satisfies the Schwarz inequality for the turbulent shear stresses. The function form of c_μ is given by:

c_{\mu }=\left(A_{0}+A_{s}U^{(\ast {})}{\frac {k}{\varepsilon }}\right)^{-1}\ \ \ \ \ \ \ \ (16)

where:

U^{(\ast {})}={\sqrt {S_{\mathit {ij}}S_{\mathit {ij}}+\Omega _{\mathit {ij}}\Omega _{\mathit {ij}}}}\ \ \ S_{\mathit {ij}}={\frac {1}{2}}\left({\frac {\partial U_{i}}{\partial x_{j}}}+{\frac {\partial U_{j}}{\partial x_{i}}}\right)\ \ \Omega _{\mathit {ij}}={\frac {1}{2}}\left({\frac {\partial U_{i}}{\partial x_{j}}}-{\frac {\partial U_{j}}{\partial x_{i}}}\right)\ \ (17)

A_{s}={\sqrt {6}}\cos \phi \ \ \phi ={\frac {1}{3}}\arccos \left({\sqrt {6}}W\right)\ \ \ W=2^{3/2}{\frac {S_{\mathit {ij}}S_{\mathit {jk}}S_{\mathit {ki}}}{S^{3}}}\ \ S={\sqrt {2S_{\mathit {ij}}S_{\mathit {ij}}}}\ \ \ \ (18)

and A₀ is a constant equal to 4.04.

Secondly, a different ε-equation is used to resolve the problem of the round-jet/plane-jet anomaly (see Pope [67]):

{\frac {\partial }{\partial x_{j}}}\left({\overline {\rho }}U_{j}\varepsilon \right)={\frac {\partial }{\partial x_{j}}}\left[\left(\mu +{\frac {\mu }{\sigma _{\varepsilon }}}\right){\frac {\partial \varepsilon }{\partial x_{j}}}\right]+c_{\varepsilon 1}{\overline {\rho }}S{\frac {\varepsilon }{k}}-c_{\varepsilon 2}{\overline {\rho }}{\frac {\varepsilon ^{2}}{k+{\sqrt {\nu \varepsilon }}}}+S_{\varepsilon B}\ \ \ \ (19)

c_{\varepsilon 1}={\mathit {max}}\left[0.43,\eta /(\eta +5)\right]\ \ \ \ \eta =Sk/\varepsilon \ \ \ \ \ \ \ \ (20)

where S is the strain-rate invariant as before, c_ε2 = 1.9, σ_k = 1.0 and σ_ε = 1.3.

The Shih et al. [66] model was developed for high Reynolds number turbulent flows and therefore a zonal or wall-function approach must be used to bridge the viscous sub-layer near walls. Compared to the standard k–ε model, it has been shown to produce improved behaviour in a number of free shear flows, boundary-layer flows and a backward-facing step flow [66]. One of the major weaknesses of the standard k–ε model is that it produces too much turbulent kinetic energy at stagnation points [68]. The Shih et al. model should in principle suffer less from this weakness since the functional form of c_μ should reduce the over-production of k. However, its overall performance in stagnating flows will depend on the type of wall model used.

Production due to Buoyancy, G

The term G in the k-equation relates to the influence of buoyancy on the turbulent kinetic energy, and is given by:

G={\overline {\rho u_{j}}}g_{j}\ \ \ \ \ \ \ \ \ \ \ \ (21)

where g_j is the gravitational acceleration vector. In stably stratified flows, where the temperature increases with height, G is negative. Conversely, in unstably stratified flows, where temperature decreases with height, G is positive and acts to increase k. The unknown density-velocity correlation, ${\overline {\rho u_{j}}}$ , must be modelled. The most common approximation of this term is the so-called Boussinesq Simple Gradient Diffusion Hypothesis (SGDH):

{\overline {\rho u_{j}}}=-{{\frac {\mu _{t}}{\sigma _{t}}}{\frac {1}{\overline {\rho }}}{\frac {\partial {\overline {\rho }}}{\partial x_{j}}}}\ \ \ \ \ \ \ \ \ \ (22)

The production due to buoyancy using SGDH is then as follows:

G=-{{\frac {\mu _{t}}{\sigma _{t}}}{\frac {1}{{\overline {\rho }}^{2}}}{\frac {\partial {\overline {\rho }}}{\partial x_{j}}}\rho _{\infty }g_{j}}\ \ \ \ \ \ \ \ (23)

In their paper, Van Maele & Merci [2] erroneously included an additional pressure-gradient term $\left(\partial P/\partial x_{j}\right)$ in Equation (23) related to the pressure-work rather than buoyancy (see Wilcox [69]). Since the term is negligible in incompressible flows, such as the buoyant plumes considered here, it has therefore been ignored. The ratio of the reference density to the mean density, $\rho _{\infty }/{\overline {\rho }}$ , appears in Equation (23) due to the use of a non-Boussinesq approach and Favre-averaging, which are discussed later. Van Maele & Merci [2] assumed that σ_t was constant and equal to 0.85.

Instead of writing the buoyancy production in terms of the density-velocity correlation, ${\overline {\rho u_{j}}}$ , the equation can be written in terms of the heat flux, ${\overline {u_{j}t^{\prime }}}$ :

{\overline {u_{j}t^{\prime }}}=-{{\frac {\mu _{t}}{\sigma _{t}}}{\frac {\partial T}{\partial x_{j}}}}\ \ \ \ \ \ \ \ \ \ (24)

and the G term is then written:

G=-{\beta {\frac {\mu _{t}}{\sigma _{t}}}{\frac {\partial T}{\partial x_{j}}}g_{j}}\ \ \ \ \ \ \ \ \ \ (25)

where t′ is the temperature fluctuation, T is the mean temperature and β is the volumetric expansion coefficient, $\beta =-\left(1/{\overline {\rho }}\right)\partial {\overline {\rho }}/\partial T$ . Other equivalent expressions can also be formulated using the ideal gas law and the assumption that density is only a function of temperature, not pressure (the low-Mach-number approximation). The conversion from mean density to temperature gradients is then as follows:

{\frac {\partial {\overline {\rho }}}{\partial x_{j}}}={\frac {\partial {\overline {\rho }}}{\partial T}}{\frac {\partial T}{\partial x_{j}}}={\frac {-P_{\ast {}}}{RT^{2}}}{\frac {\partial T}{\partial x_{j}}}={\frac {\overline {\rho }}{T}}{\frac {\partial T}{\partial x_{j}}}\ \ \ \ \ \ (26)

The SGDH model predicts zero density-velocity correlation or heat flux components ( ${\overline {\rho u_{j}}}=0$ or ${\overline {u_{j}t^{\prime }}}=0$ ) in situations were the density or temperature gradients are zero in that direction. However, as Ince & Launder [70] noted, in a simple shear flow in which there are only cross-stream temperature gradients, the heat flux in the streamwise direction actually exceeds that in the cross-stream direction. This shortcoming of the SGDH model was confirmed by the analysis of Shabbir & Taulbee [33], who showed that the model significantly underestimates the magnitude of the heat flux in vertical buoyant plumes. The underprediction of ${\overline {\rho u_{j}}}$ or ${\overline {u_{j}t^{\prime }}}$ by the SGDH model leads to an overly-small production term, G, and hence a turbulent kinetic energy, k, which is too small, producing too little mixing in the modelled plume. The study by Yan & Holmstedt [53] provides a clear example of how the k – ε model with SGDH produces buoyant plumes which are too narrow and with overly high temperatures and velocities in the core of the flow.

Van Maele & Merci [2] examined a different model for G based on the the Generalized-Gradient Diffusion Hypothesis (GGDH) of Daly & Harlow [52]. This was first used in the context of practical CFD calculations with the k – ε model by Ince & Launder [70], and is written as follows:

{\overline {\rho u_{j}}}=-{{\frac {3}{2}}{\frac {c_{\mu }}{\sigma _{t}}}{\frac {k}{\varepsilon }}\left({\overline {u_{j}u_{k}}}{\frac {\partial {\overline {\rho }}}{\partial x_{k}}}\right)}\ \ \ \ \ \ \ \ \ \ (27)

which Van Maele & Merci [2] expressed as follows:

G=-{{\frac {3}{2}}{\frac {\mu _{t}}{\sigma _{t}{\overline {\rho }}^{2}k}}\left({\overline {u_{j}u_{k}}}{\frac {\partial {\overline {\rho }}}{\partial x_{k}}}\right)}\rho _{\infty }g_{j}\ \ \ \ \ \ \ \ (28)

again with σ_t = 0.85. As previously, Van Maele & Merci [2] included a pressure-gradient term in the above equation related to pressure-work but this can effectively be ignored in the present application. The advantage of the GGDH approach is that transverse density gradients affect the production term.

In their plume simulations, Van Maele & Merci [2] used a slightly modified form of the above relation. They replaced the normal stress in the streamwise (vertical) direction, ${\overline {\mathit {ww}}}$ , with the turbulent kinetic energy, k. They justified this on the basis that the k – ε model gives poor predictions of normal stresses in plumes. Experimental measurements indicate that the streamwise normal stress is approximately twice the magnitude of the transverse components (i.e. ${\overline {\mathit {ww}}}\approx 2{\overline {\mathit {uu}}}$ ) whereas the k – ε model predicts them to be roughly equal. Since k can be approximated from $k\approx {\frac {1}{2}}\left({\overline {\mathit {ww}}}+2{\overline {\mathit {uu}}}\right)$ and the k – ε model predicts ${\overline {\mathit {ww}}}\approx {\overline {\mathit {uu}}}$ , they suggest that it is more appropriate to use k rather than ${\overline {\mathit {ww}}}$ , to artificially increase the stress to a more realistic value. This ad-hoc correction may not be appropriate in more complex flows where the gravitational vector is not aligned to the Cartesian axes.

A simplification frequently made to the buoyancy treatments described above is to assume that the mean density is equal to the reference density, ${\overline {\rho }}\approx \rho _{\infty }$ , an approach known as the Boussinesq approximation. For the SGDH written in terms of temperature gradients, this gives:

G=-{\beta {\frac {\mu _{t}}{\sigma _{t}}}{\frac {\partial T}{\partial x_{j}}}g_{j}}\ \ \ \ \ \ \ \ \ \ (29)

Van Maele & Merci [2] examined the effect of this simplification on the prediction of the George et al. [3] buoyant plume experiments.

Buoyancy Source Term in the ε–Equation, S_εB

The buoyancy source term, S_εB, in the ε–equation is given by:

S_{\varepsilon B}=c_{\varepsilon 1}\left(1-c_{\varepsilon 3}\right){\frac {\varepsilon }{k}}G\ \ \ \ \ \ \ \ (30)

Unlike other model constants in the k – ε model, there is still some controversy over the best value or formula for c_ε3. Different approaches have been proposed by different researchers, partly depending on whether the flows are horizontal or vertical and whether there is stable or unstable stratification. For a review of the performance of various models, see Rodi [71], Markatos et al. [72] or Worthy et al. [73]. In their paper, Van Maele & Merci [2] provided a summary of the values proposed previously in 20 published papers and, based on analysis of these studies, used a constant value for c_ε3 of 0.8.

Favre-Averaging

Throughout their paper, Van Maele & Merci [2] refer to Favre-averaged mean velocity, enthalpy and temperature ~T

( ${\tilde {U}}$ , ${\tilde {h}}$ and ${\tilde {T}}$ ) instead of the perhaps more familiar Reynolds-averaged values, (U, H and T). The Favre average of a variable ${\phi }$ , denoted ${\tilde {\phi }}$ is calculated from: