Buoyancy-opposed wall jet

Application Challenge 3-01 © copyright ERCOFTAC 2004

Overview of CFD Simulations

Introduction

A brief overview of the computational approach is given in paragraph 1.1 of this report. The preliminary computational studies described in issue 1 had considered only the isothermal case. This enabled an early evaluation of the different turbulence models considered, as well as providing verification that the implementation of the computational methods was correct. The more extensive studies, reported here, covered both isothermal and non-isothermal cases at two different velocity ratios (V_chan / V_jet 0.077 and 0.15). Moreover, two alternative wall functions have been used with all the high-Reynolds-number models. For the ratio V_chan / V_jet of 0.077, LES computations have also been carried out by Prof. Laurence’s group.

RANS Turbulence Modelling

The four turbulence models used in more detail are:

Model 1: Eddy-Viscosity scheme with wall functions (KEWF)

This model applies the standard eddy-viscosity method, wherein transport equations are solved for the turbulence energy, κ and its dissipation rate ε and from the values of κ and ε values of the turbulent viscosity μ_t are obtained. When the turbulent viscosity is multiplied by the mean strain rate, one obtains the Reynolds stresses. This model, however, cannot directly handle the low Reynolds number region near the wall and necessitates the use of a wall function to define the behaviour in the near-wall control volumes. In addition to the conventional wall-function strategies, that assume a logarithmic variation for the near-wall velocity, a novel approach is also tested here, in which the near-wall velocity variation is obtained from the analytical solution of the near-wall form of the momentum equation. Both these strategies are further discussed later on in this section.

Model 2: Low-Reynolds-number Eddy-Viscosity scheme (LRKE)

One need not, necessarily, be reliant on wall functions to compute the flow behaviour near to a rigid wall. If one adapts the turbulence model such that it can properly compute the flow in regions of low turbulent Reynolds numbers, one may solve the transport equations through the viscous sublayer, right up to the wall, without the need for a truncated grid and wall functions. The low Reynolds-number eddy-viscosity scheme used here is essentially that of Launder and Sharma (1974). There is a small modification to the transport equation for the dissipation ε, from that originally proposed. Further details are provided in Media:Appendix_A_issue2e.pdf

Model 3: Basic Second Moment Closure (Basic2mc)

The second pair of turbulence closures solve transport equations for the four non-zero Reynolds stresses <μ_iμ_j>. In producing exact transport equations for these terms, one arrives at some terms which can be implemented exactly and others terms which require some form of turbulence model in order to close the equations. The first closure considered here is generally referred to as the ‘Basic Model’ adopting the Launder, Reece and Rodi (1975) proposal to handle the pressure strain j_ij, assumes isentropic dissipation ε_ij, and used the generalised gradient diffusion hypothesis of Daly and Harlow (1970) for the diffusion d_ij. Two additional features are needed in order to handle the effects of the rigid wall correctly. Firstly, since this is essentially a high Reynolds number scheme, wall functions are required as described above for the KEWF model. Secondly, one needs to adapt the pressure-strain model j_ij, introducing a wall-reflection component ω^w_ij. The form adopted here is that of Craft and Launder (1992), which is a development of the earlier and better known Gibson and Launder (1978) model, but which has been designed to handle flows impinging on, as well as parallel to the wall.

Model 4: Two-Component Limit Second Moment Closure (TCL2mc)

The fourth turbulence model considered here is a more advanced form of second-moment closure. Workers at UMIST over the last decade have evolved a model for the pressure-strain j_ij and dissipation ε_ijprocesses (Kidger (2000)). This model exhibits two features: it is implicitly realisable, in that it cannot produce non-physical values for the Reynolds stresses, particularly when the turbulence approaches a two-component state as occurs near a wall or free-surface. Its second feature is that it eliminates the need for a wall-reflection process, or more importantly, the need to define a wall-normal vector y. In its first decade of testing this TCL (i.e. two-component limit) model has been found to yield superior results to the Basic model for both wall jets and flows which exhibit buoyant effects, both of which are very relevant to the present study (Craft et al, 1996; Craft and Launder, 2001).

Modelling of the Turbulent Heat Fluxes

When the linear eddy-viscosity model is used for momentum transport, the scalar fluxes are consistently represented by

{\overline {u_{i}\theta }}=-{\nu _{t} \over {\sigma _{\theta }}}{{\partial \Theta } \over {\partial x_{i}}}

With a second-moment closure, however, rather than the differential transport equations for ${\overline {u_{i}\theta }}$ , the so-called generalized gradient diffusion hypothesis (GGDH) is used instead:

{\overline {u_{i}\theta }}=-c_{\theta }\ {k \over \varepsilon }\ {\overline {u_{i}u_{j}}}\ {{\partial \Theta } \over {\partial x_{j}}}

The experience at UMIST is that this form is usually adequate as long as the influential turbulent scalar and momentum fluxes are horizontal. It is emphasized, however, that in horizontally directed shear flows affected by buoyancy where the turbulent fluxes are vertical, it would be essential to use a complete second-moment closure, as described in Craft et al (1996). Indeed, in situations of extreme stable stratification third-moment closure is essential, Craft and Launder (2002), Kidger (2000).

Wall Function Strategies

In near-wall flows, the turbulence field undergoes rapid changes across a very thin near-wall sub-layer, known as the viscous sub-layer. Integration of the mean flow equations from the fully turbulent outer region to the wall is possible, but requires the use of specially extended models of turbulence, known as low-Reynolds number models, and also very fine near-wall meshes. The latter requirement makes this approach very expensive, especially in three-dimensional flows. A more economical alternative, the wall-function approximation, has consequently been developed and has become the most widely used practice for the mathematical modelling of the effects of near-wall turbulence. Its popularity lies in the fact that it does not need to resolve the viscous sub-layer. The near-wall grid node is located far enough from the wall to be in the fully turbulent region. The semi-empirical co-relations are then used to account for the wall effects on turbulence, resulting in estimates for the wall shear stress and also for the generation and dissipation rates of turbulence over the wall-adjacent control volumes. The main drawback of such strategies, at least until recently, has been the assumption that the near-wall region is one of uniform shear stress and with the production and dissipation rate of k in balance. This has limited the range of flows for which the use of wall functions can result in reliable predictions. A more refined alternative has been recently developed by the UMIST group, Craft et al (2002), which relies on assumptions made at a deeper, more generally valid level. In the present work this alternative has been used in combination with the three high-Reynolds-number models, in addition to the standard wall function. Both the standard and the Analytical (UMIST-A) wall functions are briefly outlined below.

Standard Wall-Function

Figure 11. Prescribed semi-logarithmic velocity and temperature profiles

In this widely used approach, the variation in velocity between the near-wall node and the wall is assumed to follow the log-law. The total shear stress is assumed to remain constant and equal to the wall shear stress. The turbulent shear stress is assumed to be constant (also equal to the wall shear stress) up to the edge of the viscous sub-layer where it is assumed to abruptly fall to zero. Outside the viscous sub-layer, the dissipation rate of turbulence is assumed to be proportional to the inverse of the wall distance, and within the viscous sub-layer it is assumed to remain constant. Moreover the dimensionless thickness of the viscous sub-layer is also assumed to be independent of the mean flow field. Inevitably this rather large set of assumptions makes this approach unsuitable for many flows, including those affected by buoyancy.

Analytical Wall Function, Craft et al (2002).

Figure 12. Prescribed variation of turbulent viscosity

The overall strategy is the same as in the standard wall function, namely the use of a large wall-adjacent control-volume and the wall shear stress and the generation and dissipation rates of turbulence estimated from certain assumptions. The approximations involved, however, are not as limiting as those of the standard wall-function. The wall shear stress and the appropriate thermal boundary condition are obtained through the analytical solution of simplified near-wall versions of the transport equations for the wall-parallel momentum and the thermal energy respectively. These analytical solutions are obtained by prescribing the variation of the turbulent viscosity and also by assuming that convective transport and the wall-parallel pressure gradient do not change across the near-wall control volumes. Over an inner viscosity-dominated sub-layer, determined by the value of the dimensionless wall distance y* (º ykP1/2/n), the turbulent viscosity is assumed to be zero. Outside the viscosity-dominated sub-layer, the turbulent viscosity is taken to be proportional to the distance from the edge of this viscosity-dominated sub-layer. The dimensionless thickness of the viscosity-dominated sub-layer is determined with reference to fully developed pipe flows. The analytical velocity variation and the prescribed turbulent viscosity lead to an expression for the local generation rate of turbulence, which can then be integrated over the wall-adjacent control volume. Further refinements have also been added to this approach. In heated flows the molecular viscosity is allowed to vary with temperature across the viscosity-dominated sub-layer. The analytical solution of the momentum equation also takes into account the buoyancy force, which is introduced through the analytical temperature variation. The dissipation rate of turbulence over the near-wall control volume is made sensitive to the mean flow development through the introduction of a parameter based on the ratio between the wall shear stress and the shear stress at the edge of the viscosity-dominated sub-layer. The wall-parallel momentum flux across faces of wall-adjacent control volumes is evaluated from the analytical velocity variation. In high Prandtl number fluids an effective molecular Prandtl number is introduced, to correct for the fact that the conduction sub-layer is thinner than the viscosity-dominated sublayer. These refinements make the analytical wall function applicable over a wider range of flows than the standard approach.

The equations involved together with the physical basis of the assumptions are included in Appendix B.

LES

LES computations have been carried out using two different codes. One is a research code, Saturne, developed at EDF and the other is the commercial code Star-CD. Both codes use second-order centered schemes in time and space and they employ pressure correction methods for the velocity-pressure coupling. Both the classical, Smagorinsky, and the dynamic, Germano et al sub-grid models have been tested.

Computational grid

A computational domain was chosen which covered the total flow measurement region (Figure 1). This necessitated having a cut-out in the flow to handle the region of the thick splitter plate. The inlet profiles of the two fluid streams were not reported in the experimental work, and the strategy adopted here has been to assume a ‘top-hat’ velocity profile both at entry to the splitter plate and at the bottom of the test section where the counter-flowing main stream enters. It is assumed that the solver will yield an acceptable boundary layer profile by the time the jet exits into the main cavity at the bottom of the splitter plate and where the channel flow reaches the mixing region.

The above strategy, however, has shortcomings, because of the presence of the splitter and the fact that the TEAM code uses a structured staggered grid. At the edges of the computational domain, the boundaries are arranged in such a way as to make the application of the appropriate boundary conditions easy since all dependent variables have a node on the boundary. However, this arrangement becomes somewhat complicated as a result of having the splitter plate within the computational domain since there will not be any nodes on the boundary on the splitter surface. For that reason, the strategy adopted was to use two distinct calculations, the first to compute the developing profile of the inlet jet between the splitter and the outer wall and the second to resolve the flow in the primary domain taking the values from the former to define the inlet jet.

Table 2. Isothermal Flow Computations
V_ch/V_jet	T_ch °C	T_jet °C	Turbulence Model	Wall Function	Flow Domain		Grids
					Y_L	Y_U
O.077	42	42	Low-Re κ-ε	NA	1.3m	0.4m	130x280 103x165
0.077	42	42	High-Re κ-ε	Standard	1.3m	0.4m	102x280 75x165
0.077	42	42	High-Re κ-ε	Analytical	1.3m	0.4m	102x280 75x165
0.077	42	42	Low-Re κ-ε	NA	1m	1m	130x300
0.077	42	42	High-Re κ-ε	Standard	1m	1m	102x300
0.077	42	42	High-Re κ-ε	Analytical	1m	1m	102x300
0.077	42	42	Basic 2mc	Standard	1.3m	0.4m	102x280 75x165
0.077	42	42	Basic 2mc	Analytical	1.3m	0.4m	102x280 75x165
0.077	42	42	TLC 2mc	Standard	1.3m	0.4m	102x280 75x165
0.077	42	42	TLC 2mc	Analytical	1.3m	0.4m	102x280 75x165
0.077	42	42	LES	NA	1m	0.4m	58x142x27 59x141x40 601.6x10³†
0.15	42	42	Low-Re κ-ε	NA	1.3m	0.4m	130x280 103x165
0.15	42	42	High-Re κ-ε	Standard	1.3m	0.4m	102x280 75x165
0.15	42	42	High-Re κ-ε	Analytical	1.3m	0.4m	102x280 75x165
0.15	42	42	Basic 2mc	Standard	1.3m	0.4m	75x165
0.15	42	42	Basic 2mc	Analytical	1.3m	0.4m	75x165
0.15	42	42	TLC 2mc	Standard	1.3m	0.4m	75x165
0.15	42	42	TLC 2mc	Analytical	1.3m	0.4m	75x165

YL and YU are the lengths of the computational domain below and above the splitter plates respectively

† Unstructured Mesh

Table 3. Non-Isothermal Flow Computations
V_ch/V_jet	T_ch °C	T_ch °C	Turbulence Model	Wall Function	Flow Domain		Grids
					Y_L	Y_U
0.077	34	42	Low-Re κ-ε	NA	0.6m	1.8m	130x280 130x150
0.077	34	42	High-Re κ-ε	Standard	0.6m	1.8m	102x280 75x150
0.077	34	42	High-Re κ-ε	Analytical	0.6m	1.8m	102x280 75x150
0.15	34	42	Low-Re κ-ε	NA	0.6m	1.8m	130x280 103x150
0.15	34	42	High-Re κ-ε	Standard	0.6m	1.8m	102x280 75x150
0.15	34	42	High-Re κ-ε	Analytical	0.6m	1.8m	102x280 75x150
0.15	34	42	Low-Re κ-ε	NA	0.6m	2.4m	130x300 103x150
0.15	34	42	High-Re κ-ε	Standard	0.6m	2.4m	102x300 75x150
0.15	34	42	High-Re κ-ε	Analytical	0.6m	2.4m	102x300 75x150

Computations have been carried out with different lengths of the computational domain, both below (Y_L) and above (Y_U) the splitter plate, in order to ensure that the boundary locations did not affect the computations. In the isothermal case the lower entry plane was located 1.3m below the splitter plate, so as not to influence the predicted penetration length of the jet. In non-isothermal flows the upper exit boundary was moved to 1.8 m above the splitter plate to prevent back flow across the exit boundary.

Different values for the turbulent kinetic energy and its dissipation rate have also been tried at the lower inlet boundary. The same case (isothermal, V_chan /V_jet=0.077) was computed with either k_in = 0.01V_in² and ε_in = rc_μκ_in²/(50μ), or κ_in = 0.025V_in² and ε_in = κ_in^3/2 / l_in, where l_in is prescribed from DNS data for channel flow. The resulting predictions were identical.

With regard to the mesh spacing adjacent to the wall, with the low-Re model a maximum dimensionless near-wall mesh size (y*) of 1 at the jet exit, based on the near-wall turbulent kinetic energy, has been used. In the high-Re computations, mesh sizes have been used which result in y* values of between 65 and 160 at the jet exit, with no discernible differences between the predictions. On the wall opposite the jet, mesh sizes have been used which result in y* values of the order 5 to 10 at the bottom of the domain and 40 to 80 at the top, again with no discernible differences between the predictions.

In the LES computations three grids have been tested, two of which are fully structured and one that involves local grid refinement.

The cases computed are tabulated in Tables 2 and 3. In all cases the jet Reynolds number is 4,000.

Contributors: Jeremy Noyce - Magnox Electric

Site Design and Implementation: Atkins and UniS