Thermophysics and Aeromechanics, 2016, Vol. 23, No. 5

DOI: 10.1134/S0869864316050061

Turbulent circulation above the surface heat source in a stably stratified environment* 1

A.F. Kurbatskii and L.I. Kurbatskaya

2

1

Khristianovich Institute of Theoretical and Applied Mechanics SB RAS, Novosibirsk, Russia

2

Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Novosibirsk, Russia

E-mail: [email protected] (Received December 27, 2015; in revised form January 12, 2016) The results of the numerical modeling of turbulent structure of the penetrating convection above the urban heat island with a small aspect ratio in a stably stratified medium at rest are presented. The gradient diffusion representations for turbulent momentum and heat fluxes are used, which depend on three parameters ⎯ the turbulence kinetic energy, the velocity of its spectral expenditure, and the dispersion of temperature fluctuations. These parameters are found from the closed differential equations of balance in the RANS approach of turbulence description. The distributions of averaged velocity and temperature fields as well as turbulent characteristics agree well with measurement data. Keywords: turbulence, planetary boundary layer, urban heat island, large-scale circulation, numerical modeling.

Introduction The air circulation above the urban heat island is generated by means of the energy of anthropogenic sources within the urban area (Fig. 1). It becomes the most intense at the nighttime at a clean sky and weak wind. The distributions of velocity and temperature fields as well as the turbulence intensities are the fundamental characteristics reflecting the structure of the night urban heat island. The vertical turbulent thermal plume from the localized source (the heat island) and the circulation related to it develop because of a temperature difference between the heat island and its environment, which has a lower temperature. A typical thermal plume of the urban heat island has a small aspect ratio ( zi /D << 1) that is a small ratio of the mixing height ( zi ) to the linear size of the heat island (D) (Fig. 2). In the measurements [1], the mixing height zi is defined as a height at which the maximum negative difference is reached between the temperature at the plume center and the ambient temperature (density) as shown in Fig. 2. The full-scale measurements in the night boundary layer show the following typical values: zi ≈ 200 m, D ≈ 20 km. The smallness of the ratio zi /D *

The work was financially supported by the Russian Foundation for Basic Research (Grant No. 13-05-00006 and partially No. 14-01-00125) and Program for Basic Research of SB RAS No. II. 2П.

© A.F. Kurbatskii and L.I. Kurbatskaya, 2016

677

A.F. Kurbatskii and L.I. Kurbatskaya Fig. 1. Diagram of the realization in laboratory experiment [1] of the circulation above the urban heat island in a stably stratified urban boundary layer at a weak wind. On the middle fragment on the right ⎯ the heated plate in the form of a disc, which simulates the prototype (on the left); the shadow photograph at the bottom shows a thermal plume above the localized heat source (the heated plate) in the quasi-steady state; the work medium in the laboratory setup is the thermally stably stratified water [1, 2].

for the prototype presents the main limitation of the laboratory experiment. Due to a considerable reduction of the linear scale in laboratory modeling it is difficult to resolve the flow structure in the entire range of heights ⎯ from the surface to the plume top boundary. A typical Reynolds number of the flow to be modeled is by several orders of magnitude less than for the prototype. Therefore, the laboratory modeling is restricted by a reproduction of the large-scale circulation above the heat island inside the urban boundary layer; the processes in the viscous sublayer (in the roughness layer near the ground surface) are not considered. The turbulent motion at the urban heat island center, however, is predominant at a weak wind because of work of a fluctuating buoyancy force. The mechanical action of the urban roughness and the wind shear are less important. To obtain the thermal plume with a small aspect ratio in the laboratory setup [1] (120×124×45 cm) it was required to create large heat fluxes from the heated surfaces (Fig. 1) and very strong temperature gradients in the water, which was used for modeling the turbulent convection, as was first proposed to model the convection in the planetary boundary layer in the work [2]. But the authors of the work [1] failed to obtain a typical region of small values of the ratio zi /D and a small Froude number Fr (the definition of the latter will be given in the following), which are expected for a real heat island in the night urban boundary layer. In the experiment, the values zi /D ≅ 0.09 and Fr = 0.03 were reached, which hardly reached the upper values in the atmosphere. An important consequence of maintaining the above conditions in experiment was, in particular, the fact that the mixing height zi , the velocities and temperatures were the functions of time after the heating was started. After several minutes from the moment of the disc heating start, a quasi-steady state was observed because the heat island intensity ΔTm = Tm − T0 and the heat flux from the surface remained constant in time.

Fig. 2. Diagram of the circulation above the heat island. On the left ⎯ the horizontal velocity distribution, on the right ⎯ the vertical profile of temperature in the ambient medium; zi ⎯ the mixed layer height, ze ⎯ the equilibrium height, Tm ⎯ the temperature on the thermal plume axis, Ta ⎯ the temperature of the ambient medium near the surface, Tp ⎯ the temperature in the ambient medium at the mixed layer height.

678

Thermophysics and Aeromechanics, 2016, Vol. 23, No. 5

This experiment was in fact limited by the first several minutes of a quasi-steady state until the induced circulation reached the setup side walls and the latter started affecting the circulation. The numerical model of the urban heat island is formulated with regard for the abovenoted limitations in laboratory modeling. In the present study, the RANS approach is used for modeling the turbulent circulation above the heat island and reproducing the velocity and temperature fields as well as the turbulence structure. This approach includes the averaged equations of thermohydro-dynamics. For turbulent momentum and heat fluxes, the gradient diffusion equations are formulated, which depend on three parameters: the turbulence kinetic energy, the rate of its spectral expenditure (dissipation), and the dispersion of turbulent temperature fluctuations. These parameters are determined from the solution of closed differential equations of balance [3], which are given in the Appendix. The main difference of the present heat island model from other theoretical and numerical models lies in the description of turbulent convection. So, for example, the two-dimensional model [4] accounts for the differences in the urban roughness and the ambient medium, however, the turbulence is handled in a parameterized form, which does not enable a complete consideration of the structure of the evolving thermal plume and the obtaining of the distributions in space and time for the main characteristics of the heat island. Such distributions are important for modeling the atmospheric diffusion of admixtures. 1. Governing equations The penetrative turbulent convection initiated by a constant heat flux H 0 from the surface of a plate of diameter D (Figs. 1 and 2) models in the laboratory experiment a prototype of an urban heat island with a small aspect ratio ( zi /D << 1 ) at a weak wind in a stably stratified ambient medium: ( ∂Ta /∂ z ) = γ z > 0. The medium motion is considered as possessing the axial symmetry. Based on the processing of measurement data [1] for

(

the horizontal (radial) velocity the scale wD = βgDH0 /ρ0c p

)

has been chosen, where β is

the coefficient of thermal expansion of the medium, g is the gravity acceleration, cp is the specific heat. The quantity wD slightly differs from the convective velocity scale introduced by Deardorff [2] because the diameter D of the heated plate is used as the basic controlled length parameter. The quantity wD N / g β is taken as the reference temperature value, and the quantity D/ wD is used as the time scale. The Froude number Fr = wD /( N ⋅ D), the quantity N = βg (∂T /∂z )0 

1/ 2

is the Brünt−Väisälä frequency.

The governing equations of fluid dynamics for describing the circulation above the heat island of a small relative aspect ratio can be taken in the hydrostatic approach [5]. Without regard for the Coriolis force and radiation, these equations have the following dimensionless form for the averaged values of the velocity and temperature in the Boussinesq approximation: h

∂U 1 ∂ ∂ ∂T ∂ (u 2 ) ∂ (uw) u 2 − v 2 + ⋅ ( rU ) U + Fr UW = Fr −1  dz− − + + ∂ t r ∂r ∂z ∂r ∂r ∂z r z  1 ∂ ∂U ∂ 2U + Re −1  ⋅ r +  r ∂r ∂r ∂ z2 

∂ ∂ rU + Fr rW = 0, ∂r ∂z

  , 

(1)

(2)

679

A.F. Kurbatskii and L.I. Kurbatskaya

 1 ∂ ∂T ∂ 2T ∂T 1 ∂ ∂ + ( rU T ) + Fr (W T ) = (Re⋅ Pr)−1  ⋅ r + 2 ∂r r ∂r ∂z  r ∂r ∂r ∂ z

  1 ∂  ∂ wθ  , (3)  +  − ⋅ r uθ − ∂z    r ∂r

here U is the averaged horizontal velocity, W is the averaged vertical velocity, u is the horizontal turbulent velocity, w is the vertical turbulent velocity, v is the azimuthal turbulent velocity, T is the averaged temperature, θ is the turbulent fluctuation of temperature, Re = = (wD D)/ν is the Reynolds number, λ is the thermal diffusivity coefficient, ν is the kinematic viscosity, Pr = ν /λ is the Prandtl number, h is the given height of the stratified layer of the medium. In equations (1)−(3) and in the following, the lower-case letters supplied with a bar denote the averaged values of quantities, and the capital letters denote the turbulent fluctuations of quantities. 2. Turbulent fluxes of heat and momentum

Equations (1)−(3) include as the unknown quantities the components of turbulent stresses ui u j (the normal and shear stresses) and the vector of the turbulent heat flux uiθ . One can obtain a physically correct description of the effect of the stratification on the circulation turbulent structure above the heat island in the approach of the locally equilibrium turbulence [3], and the turbulent fluxes are parameterized by the gradient diffusion relations. Such an approach minimizes the complexity of the description of turbulence of stratified flows of the environment and the costs of the numerical realization of the model. The exact equation for the vector of the turbulent heat flux uiθ has the form [6]:

Duiθ ∂ = ∂xj Dt −u jθ

 ∂u ∂θ ν ⋅θ i + λ ⋅ ui  ∂xj ∂xj 

 ∂ ∂T − ui u jθ − ui u j −  ∂xj ∂xj 

∂U i θ ∂ p ∂ u ∂θ − ⋅ − ( λ +ν ) i ⋅ − β gi θ 2 . ∂ x j ρ ∂ xi ∂ xk ∂ xk

(4)

For the case of large turbulent Reynolds numbers, equation (4) in closed form has the form Duiθ ∂ − Dt ∂ xk

 ∂ uk θ ∂ uiθ   ∂U i E ∂T + uk uα − u jθ −   = −ui u j csθ  ui uα ε x x x ∂ ∂ ∂ ∂xj j α α    

−c1θ

uiθ 2ττ θ

+ c2θ u jθ

∂U i + c3θ gi βθ 2 − gi βθ 2 , ∂xj

(5)

where E = (1/ 2 ) ui u i is the kinetic energy of turbulence, and ε is the rate of its dissipation,

τθ = θ 2 /2εθ is the time scale of the turbulent temperature field, τ = E /ε is the time scale of the turbulent velocity field. In equation (5), the dissipation of the turbulent heat flux uiθ is assumed to be equal to zero in view of the considerations of the absence of an isotropic tensor of the first rank. It is assumed within the framework of this assumption that c3θ = c2θ ; the numerical values of coefficients c1θ , c2θ , c3θ will be given in the following. A fully explicit gradient model for the vector of the turbulent heat flux can be derived from equation (5) in the approach of the weakly equilibrium turbulence [3]. In the present study, a simplified version of the model for the vector of turbulent heat flux is used, which is obtained in the approach of the equilibrium turbulence from equation (5) originally in the implicit form because uiθ enters the right-hand side of expression (6): 680

Thermophysics and Aeromechanics, 2016, Vol. 23, No. 5

 ∂U i  ∂T 2 + (1 − c2θ )u jθ (6)  + (1 − c2θ ) gi βθ  . ∂xj ∂ x j  c1θ    The obtaining of fully explicit approachs for turbulent fluxes of the momentum and heat presents in the general three-dimensional case a difficult task [7]. In the two-dimensional case, the matrix of linear algebraic equations (6) may be inverted for turbulent fluxes of the heat and momentum by using the symbolic algebra. Such a strategy will enable the obtaining of fully explicit expressions of the gradient diffusion for turbulent fluxes [3, 8]. To write equation (6) in an explicit form one employs simple models of the Boussinesq gradient diffusion for the fluxes of the momentum ( −ui u j ) and heat ( uiθ ) entering the right-hand side of (6): −uiθ =

2

 

ττ θ  ui u j

−ui u j = 2 K m Sij −

2 ∂T Eδ ij , − u jθ = K h , 3 ∂xj

(7)

where Sij = (∂U i /∂ x j + ∂U j /∂ xi ) /2 is the mean strain tensor, K m = cμ E 2 /ε is the eddy diffusivity for momentum, K h = cλ 2R ⋅ E 2 /ε is the eddy diffusivity for heat, R = τθ /τ

is

the parameter of the ratio of time scales of the temperature ( τ θ ) and dynamic (τ) turbulent fields of temperature and velocity, respectively. The numerical values of coefficients in (6) were obtained at the solution of “limiting” problems of stratified flows [9, 10]: cμ = 0.09,

cλ = 0.095, c1θ = 3.28, c2θ = с3θ = 0.40, and R = 0.6 [11]. The substitution of (7) into (6) leads to a fully explicit expression for the vector of the turbulent heat flux: −uiθ = cλ

E2

ε

2R

∂T 2R E − ⋅ [{2 K m + (1 − c2θ ) K h }Sij + ∂ xi c1θ ε

+(1 − c2θ ) K h Ωij ]

∂T 1 − c2θ + ∂xj c1θ

2R

Ε g βθ 2 , ε i

(8)

where Ωij = (∂U i /∂ x j − ∂U j /∂ xi )/2 is the mean rotation tensor. The gradient model (8) expresses explicitly the flux vector −uiθ in terms of mean gradients, the eddy diffusivities for momentum and heat as well as the dispersion of temperature fluctuations θ 2 . Equation (8) implies the expressions for the normalized heat fluxes ⎯ the vertical one

( )

( wθ )

and the hori-

zontal one uθ :

− wθ = cλ

E2

ε

2R

 ∂W ∂U  ∂T ∂W  ∂T 2R E  − ⋅  K m  Fr + −   + (1 − c2θ ) K h Fr ∂z ∂ r  ∂ r c1θ ε   ∂r ∂ z  −

−uθ = cλ

E2

ε

2R

1 − c2θ c1θ

2R

E

ε

θ 2⋅ Fr −1 ,

 ∂U ∂T ∂W 2R E  − ⋅  Km  + Fr ∂r ∂ ∂r c1θ ε  z 

(9)

 ∂U  ∂T .   + (1 − c2θ ) K h ∂ z  ∂ z 

(10)

The substitution of (9) and (10) in (3) leads to a closed form of the equation for the mean temperature. One may note that the Boussinesq model (7) for the heat flux yields a zero horizontal heat flux for the flows with the zero longitudinal gradient of the mean temperature. In expression (10), when the gradient ∂T /∂ r vanishes, the heat flux may be different from zero 681

A.F. Kurbatskii and L.I. Kurbatskaya

(because of the second item in (10)) due to the interaction of turbulent vortices with the mean temperature gradient, which is normal to the main flow direction. For the normal turbulent stresses entering the right-hand side of equation (1), a Boussinesq simple model is employed, which preserves some anisotropy of normal stresses: 2 ∂U 2 ∂W 2 U u 2 = E − 2Km , w2 = E − 2 K m Fr , v2 = E − 2Km . (11) 3 3 ∂r ∂z 3 r The substitution of (11) in (1) enables one to write equation (1) in closed form. Expressions (11) yield with regard for the continuity equation (2) as a sum the turbulence kinetic energy. For the shear turbulent stress, the gradient diffusion model (7) has the form:  ∂U ∂W  (12) −uw = K m  + Fr . ∂r   ∂z

The vertical averaged velocity W is found in the numerical model by a quadrature from the continuity equation (2). 3. Initial and boundary conditions

The problem of the circulation development above the heat island is solved in the axisymmetric statement. The integration region represents a cylinder of the given height h. The heated plate of diameter D lies at the center of the lower cylinder surface (Figs. 1 and 2). The external boundary lies at the distance of 1.5D from the cylinder axis. At the initial moment of time, the medium is at rest and is stably stratified. According to the results of the laboratory experiment, a large-scale penetrative convection above the heat island is simulated without resolving the flow in a viscous sublayer; the noslip condition is used for the radial velocity U on the lower surface. The given heat flux from the source ( 0 < r/D ≤ 0.5 ) − ρ cp λ ( ∂T /∂ z ) z =0 = H 0

(13)

is employed as a boundary condition when writing a finite difference analog of equation (3) for the temperature in near-boundary nodes of the difference grid. Outside the heater (0.5 < r/D ≤ ≤ 1.5), the temperature is assumed equal to the given surface temperature T0 : T

z =0

= T0 ,

(14)

which agrees with the laboratory experiment conditions (the wall is thermally conductive). The boundary conditions for the second moments of the velocity and temperature fields (the turbulence kinetic energy E, its dissipation rate ε, and the dispersion of temperature fluctuations θ 2 determined by equations (A.1)−(A.4) of the Appendix) are formulated at the first computational layer from the surface (z = z1) with regard for the conditions of thermal stability of the medium separately for the regions of the unstable stratification (above the heat source) and the stable stratification of the medium (outside the heat source). For the unstable stratification, the boundary values for E and ε above the heat source (at 0 < r/D ≤ 0.5) are taken in the form [12]:

E1 = u*2

2/3 2/3   zi  z1    7 + 0.52  + 0.85 1 + 3  , | L| | L|     

(15)

3/ 2 u*2  2/3  1 + 0.5 z / | L | . ( ) i  k z1 

(16)

ε1 = 682

Thermophysics and Aeromechanics, 2016, Vol. 23, No. 5

Condition (16) agrees well with the data of the Kansas experiment [13] and satisfies the limit conditions of neutral stratification ( ε = u*3 /k z1 ) and strongly unstable stratification (ε =

(

= 0.35βgH0). The quantity L = u*3 / k β g H 0 /ρ c p

) in (15) and (16) is the Monin−Obukhov scale,

k ≈ 0.40 is the Karman constant, u* is the turbulent friction velocity, g is the acceleration due to

(

gravity. The Monin−Obukhov scale L = −u*3 / β gk ⋅ wθ

)

determines the height above which

the work of the fluctuating buoyancy force is the predominant mechanism of the generation of the turbulence kinetic energy, and the production of this energy by the velocity shear is small. The boundary condition of the dispersion θ 2 at 0 < r/D ≤ 0.5 at the first mesh point ( z = z1 ) has been obtained from equation (A.4) in the approach of the local balance “production = destruction” with the use of expression (9) for the turbulent heat flux:

θ12

(

2

)

3  Fr 2  H 0 /ρ c p β gD/wD 2 R /( cλ E1 )   = . 3  1 + 2 R (1 − c2θ ) / ( c1θ cλ ε1 )  H 0 /ρ c p β g ( D/wD )  

{

(

)

}

(17)

Outside the heater (0.5 < r/D ≤ 1.5), one specifies a “background” level of the temperature dispersion with respect to its value above the heated plate at the first computational layer (z = z1):

θ 2 ( 0.5 < r/D ≤ 1.5 )1 = 0.01⋅θ 2 ( 0 < r/D ≤ 0.5 )1 .

(18)

One may note that the variation of the numerical coefficient in (18) has affected insignificantly the behavior of the θ 2 profile at z > z1. In the work [14], a non-iterative approach was proposed for finding the surface turbulent scales of the velocity u* and temperature θ* from the Monin−Obukhov similarity hypotheses for the bottom layer of the atmospheric boundary layer [15]. In the problem under study, the scale θ* will not be required, and the formula for its computation is not presented. The turbulent friction velocity u* at the distance z = z1 from the surface (the first computational layer) is found by the formula (19) u*2 = СDs ⋅ U 2 ⋅ f m ( Ri B , z / z0 ) , Ri B = β g z

ΔT U2

(20)

,

where RiB is the volumetric Richardson number, СDs = (u*2 / U 2 ) = k 2 / ln ( z / z0 ) 

2

is

the friction coefficient for neutral conditions, z0 is the surface roughness parameter. The form of the function f m must be compatible with the logarithmic velocity profile. The necessary requirement of the compatibility is f m ( Ri B → 0, z / z0 ) = 1. For the stability function f m , the analytic approximations were obtained in the work [14]: − in the case of an unstable stratification ( Ri B < 0 )

(

f m(1) = 1 − bRi B

(1 + c | Ri

B

))

|1/ 2 ,

(21)

− in the case of a stable stratification ( Ri B > 0 ) 2

f m(2) = 1/ (1 + b1Ri B ) ,

(22)

where b = 2b 1 = 9.4, c = 7.4 CDs b ( z / z0 )1/ 2 . 683

A.F. Kurbatskii and L.I. Kurbatskaya

Thus, the turbulent friction velocity u* at the distance z = z1 from the surface may be computed by formula (19) with f m(1) above the heat source (0 < r/D ≤ 0.5) and with f m(2) outside the heat source (0.5 < r/D ≤ 1.5). Outside the surface heat source (0.5 < r/D ≤ 1.5), the boundary values for the kinetic energy of turbulence and its dissipation rate ε have the form [16]:

E1 = cμ−1/ 2u*2 ,

ε1 =

u*3

(23)

/ k z1 (1 + 4 z1/L ) .

(24)

In the computations, the boundary condition (24) has been used in the form ε1 = u*3 / k z1. According to [17], the magnitude of the Monin−Obukhov scale amounts to about the height of a stable boundary layer and corresponds to z1 << L for the conditions of the problem to be solved. The symmetry conditions

∂U ∂ r = ∂ E ∂ r = ∂ε ∂ r = ∂T ∂ r = ∂θ 2 ∂ r = 0

(25)

are specified both on the thermal plume axis (r/D = 0) and at the outer boundary (r = 1.5D), and at the top boundary (z = h), the condition of the absence of fluxes across the boundary is set:

∂U ∂ z = ∂ E ∂ z = ∂ε ∂ z = ∂θ 2 ∂ z = 0.

(26)

The boundary condition for the temperature at the top boundary is formulated in the form of the equality of the vertical gradient of temperature in two last nodes of the difference grid

( ∂T /∂ z ) z = z

J −1

= ( ∂T /∂ z ) z = z . J

(27)

This condition determines the temperature value at the upper boundary, which corresponds to the thermal stratification of the medium in the laboratory experiment ⎯ the linear law of the growth of the ambient medium temperature with increasing height. The numerical modeling was carried out for two realizations of the experiment [1]. The numerical results for the velocity field (case I) were obtained at a “weak” stratification, and the characteristics of the temperature field (case II) were obtained at a “strong” stratification: H0 = 1.81 W/cm2, wD = 1.51 cm/s, γz = 1.4°C/cm, Fr = 0.088, Re = 4500, h = 16 cm. For case I, the data have been presented in the work [1] on the measurements of mean velocities and standard deviations for turbulent fluctuations of the radial and vertical velocities, and for case II ⎯ the distributions of the mean temperature and dispersion of the temperature turbulent fluctuations. In section 5, the numerical results are compared with experimental data for the above cases I and II. 4. Numerical realization

The system of equations (1)−(3) together with the balance differential equations (A.1)−(A.4) for the second moments E, ε, and θ 2 (see Appendix) is solved under the given initial and boundary conditions by a finite-difference method using a semi-implicit scheme of alternating directions (the second upwind scheme [18]) on a difference grid with staggered 4 nodes. To reach the quasi-steady state of circulation above the heat island about 2.8⋅10 time 4 steps were required for the situation with a “weak” stratification (case I) and about 5⋅10 time steps for the situation with a “strong” stratification (case II). The computations were done on the grid with 25 nodes along the radial coordinate, and the grid had 120 nodes along the vertical 684

Thermophysics and Aeromechanics, 2016, Vol. 23, No. 5

for the layer height h = 20 cm for case I and 135 nodes for the layer height h = 16 cm for case II. The number of grid nodes in the radial direction could be varied by varying the linear size of the integration region. The time step was chosen from the condition of preserving the numerical solution accuracy. 5. Results of numerical modeling. Comparison with experiment

In experiment of [1], a quasi-steady state of thermal circulation formed after several minutes after the heat supply start from the heater to the ambient stably stratified medium. At this state, the heat island intensity  Tm = Tm − T0 and the surface heat flux remained invariable with time ( Tm and T0 are the surface temperatures at the center of the island and outside it, respectively). The experiment was limited by the first several minutes of the quasi-steady state until the medium circulation induced by the heat island reached the lateral sides of the setup and the influence of the latter manifested itself. The numerical results presented below correspond to such a quasi-steady state of the circulation above the heat island. 5.1. Velocity field

Figure 3 shows two large-scale vortices (see Fig. 1) with the counter-clockwise rotation of the left vortex and the clockwise rotation of the right vortex, which form the main ascending flow at the center, which extends up to the entrainment layer ( z / zi ≈ 1 ), and the descending flow on periphery. The thermal plume height is suppressed here by a stable stratification of environment, the lateral motion and the plume turbulence are increasing. Figure 4 shows the radial velocity profiles. Near the surface, the velocity of a flow coming from the heat island periphery grows in the direction towards the center, reaches its maximum at about r/D ≈ 0.25, and then decreases down to zero at the heat island center. The outflow velocity above increases with the distance from the center and reaches its maximum at about r/D ≈ 0.25.

Fig. 3. Streamline contours above the heat island. Experiment [1] (a) and computation (RANS approach) (b) 2 at H0 = 0.65 W/cm , Fr = 0.077, Re = 8280, (∂T/∂z)a = 0.5 °C/cm.

685

A.F. Kurbatskii and L.I. Kurbatskaya

Fig. 4. Profiles of the horizontal (radial) mean velocity U at different altitudes above the heat island. Experiment [1] (a) and computation (the RANS approach) (b) 2 at H0 = 0.65 W/cm , Fr = 0.077, Re = 8280, (∂T/∂z)a = 0.5 °C/cm; a ⎯ z/zi = 0.05 (1), 0.26 (2), 0.47 (3), 0.68 (4), 0.89 (5), 1.11 (6), 1.32 (7), b ⎯ z/zi = 0.06 (1), 0.3 (2), 0.6 (3), 0.75 (4), 0.97 (5), 1.12 (6), 1.32 (7).

The distributions of the vertical velocity across the plume at different altitudes, which are observed in experiments, and computed ones are shown in Fig. 5. The width of the ascending flow region is approximately constant in the interval of altitudes z/ zi = 0.2−0.6. One observes a stationary wave at the equilibrium altitude (z/zi ≈ 0.7). This is a consequence of the excess of the thermal plume above its equilibrium altitude in a stably stratified ambient medium, which can be seen in Fig. 6, which shows the computed velocity vector field and the temperature field isotherms. The plume upper part acquires the “hat” shape, as in the quasi-steady state on the shadow photograph of Fig. 1. The effect of the plume central part swelling will be considered below in section 5.2. The structure of the thermal plume turbulence is shown in Fig. 7 by the distributions 1/ 2   of the root-mean-square fluctuations of the horizontal  σ u = u 2 /wD  and vertical  

( )

( )

 2  σw = w  686

1/2

 /wD  turbulent velocities at the plume center. The data of laboratory measurements 

Thermophysics and Aeromechanics, 2016, Vol. 23, No. 5

Fig. 5. Normalized vertical velocity W/(wD Fr) above the heat source.

The experiment of [1] is on the left, the computation (RANS approach) is on the right; the stationary wave is fixed at an equilibrium height (z/zi ≈ 0.7) due to the excess of the plume central part (see Fig. 6).

Fig. 6. Pattern of the velocity vector field and the field of isotherms above the heat source in the quasi-steady state obtained as a result of present computations. 687

A.F. Kurbatskii and L.I. Kurbatskaya

Fig. 7. Distributions of the horizontal (σu /wD) and vertical (σw /wD) dispersion of turbulent velocity on the plume axis (r/D = 0) above the heat island 2 at H0 = 0.65 W/cm , Fr = 0.077, Re = 8280, (∂T/∂z)a = 0.5 °C/cm. 1, 2 ⎯ experiment [1], 3, 4 ⎯ the computation using the three-parameter RANS approach, 5 ⎯ the computed distribution of turbulence intensity ui ui ; 2 1/2

2 1/2

2 1/2

2 1/2

1/2

1 ⎯ u  /wD, r/D = 0 [1]; 2 ⎯ w  /wD, r/D = 0 [1]; 3 ⎯ u  /wD, r/D = 0; 4 ⎯ w  /wD, r/D = 0; 5 ⎯ (2E) .

and the numerical modeling results show that the high values of σ u /wD and σ w /wD inside the mixing layer rapidly decay with the altitude growth above the entrainment layer ( z > zi ). Both the measured and computed profiles σ w /wD have a maximum near the upper boundary of the mixing layer ( z / zi ≈ 1 ). The domain of the computed maximum value σ u /wD extends along the horizontal outside the maximum value of σ w /wD due to a horizontal divergent flow caused by a stable stratification in the plume upper part. It was noticed in the work [1] that the finite sizes if the experimental setup nevertheless affected the circulation and suppressed the horizontal motion and, thus, reduced the level of the σ u /wD values. The presented results of computing the turbulent structure of thermal plume show that in the problem under consideration, simple parameterizations of turbulent fluxes of the momentum and heat (7), (9), and (10) not only describe correctly the features of the distributions of σ u /wD and σ w /wD but also reflect fairly well their anisotropic character. Figure 8 shows the turbulent friction velocity u* ( r , z1 ) distribution versus the radial coordinate r for the quasi-steady state of turbulent circulation above the heat island. The results of

Fig. 8. Turbulent friction velocity u* near the surface at the height z = z1 as a function of the radial coordinate r in the quasi-steady state of turbulent circulation above the heat island (shown by a thick black line on the abscissa axis). 1 ⎯ RANS computation for a stratified medium, 2 ⎯ the computation for a neutrally stratified medium.

688

Thermophysics and Aeromechanics, 2016, Vol. 23, No. 5 Fig. 9. Normalized temperature profiles in different sections above the heat island. 2

H0 = 1.81 W/cm , Fr = 0.088, Re = 4500, (∂T/∂z)a = 1.4 °C/cm; experiment [1] (a) and computation (b) at r/D = 0 (1), 0.2 (2), 0.4 (3), 1.0 (4).

computation by formulas (19)−(22) are shown for the roughness parameter z0 = 0.01 cm. Outside the heat island (r/D > 0.5), the friction velocity differs little from the friction velocity for neutral conditions (Ri = 0) at a stable stratification of the ambient medium (∂T/∂z = 0.5°C/cm). 5.2. Temperature field

The vertical profiles of the dimensionless temperature (Fig. 9) point to the presence of a region of strong positive buoyancy in the plume lower part and the region of a weak negative buoyancy near the height of mixing ( z / zi ≈ 1 ). The laboratory measurements, as the full-scale observations, register a universal shape of the dimensionless vertical distribution of temperature, which is independent of the Re and Fr numbers.

The vertical temperature distributions above the heat island, which were measured in the laboratory experiment (Fig. 10a) and computed (Fig. 10b), are shown in Fig. 10. The measurement data presented in Fig. 10a correspond to a “strong” stratification (case II) at the moment of time equal to 8 minutes. The computed results shown in Fig. 10b correspond to the same conditions. The computed vertical distribution of temperature at r/D = 0 and r/D = 0.20 are similar to the measured values and point to a good mixing in the lower and central part of the plume. Such a character of the temperature distribution along the height refers to those real night planetary boundary layers,

Fig. 10. Temperature profiles in different sections above the heat island. 2

H0 = 1.81 W/cm , Fr = 0.088, Re = 4500, (∂T/∂z)a = 1.4 °C/cm; experiment [1] (a) and computation (b) at r/D = 0 (1), 0.2 (2), 0.4 (3), 1.0 (4); initial profile (5).

689

A.F. Kurbatskii and L.I. Kurbatskaya Fig. 11. Vertical profile of the dispersion of temperature turbulent fluctuations above the heat island on the plume axis (r/D = 0). 2

H0 = 1.81 W/cm , Fr = 0.088, Re = 4500, (∂T/∂z)a = 1.4 °C/cm; 1 ⎯ experiment [1], 2 ⎯ the results of computation in the RANS approach.

in which the unstable (convective) conditions are predominant, which arise because of the ascending heat flux from the urbanized surface and small motion velocities of the ambient air. In the absence of a strong heat flux directed upwards, or considerable velocities of the ambient air motion, the approximately neutral or slightly stable conditions may prevail inside the urban boundary layer. It is seen in Fig. 10a and 10b that the temperature profiles inside the plume have a typical “swelling”: the temperature inside the plume proves to be lower than the temperature outside it at the same altitude, which points to the region of a negative buoyancy arising because of the plume elevation at the center. This elevation is due to the intersection effect revealed in the vertical temperature profiles the urban town−rural. The elevation height is maximum on the plume axis and drops with the increasing distance from its center. Such a character of the vertical distribution of temperature corresponds to the shadowgraph pictures of the thermal plume in the laboratory experiment in Fig. 1: the plume has a dome-shaped upper part in the form of a “hat” (see also Fig. 6). The turbulent structure of temperature field is shown in Fig. 11 by the distribution of

(

root-mean-square fluctuations of σ T /TD , where TD = H 0 / ρ c p wD

)

is the convective scale of

the temperature field. The profile has a typical form with decreasing from the maximum value near the surface to the minimum value at z / zi ≈ 0.85 both in the experiment and in the RANS approach with the parameters E = (1/ 2)ui ui , ε, and θ 2 computed from the balance equations.

The second maximum near the entrainment layer ( z / zi ≈ 1 ) is due to a large gradient of the generation velocity θ 2 , which is confirmed by the vertical distribution of temperature in Figs. 9 and 10. The difference in the profiles of temperature dispersion in the upper part of the mixing layer ( z / zi ≈ 1 ) is related to the known defect of using the simplified approximation of gradient type of turbulent diffusion processes (the third moments of thermal and hydrodynamic fields) in the RANS approach of the turbulence modeling with the closure at the level of second moments. Conclusions

A mathematical model has been formulated for the urban heat island. The description of the thermal plume turbulence in the RANS approach including the second-order moments of thermal and hydrodynamic fields as the sought quantities enables the reproduction of the structural peculiarities of the penetrative turbulent convection above the heat island, including such fine effects as the intersection of the vertical profiles of the thermal plume temperature with the formation of a negative buoyancy region, which testifies to the development of a dome-shaped form of the plume upper part in the “hat” form (Fig. 1). 690

Thermophysics and Aeromechanics, 2016, Vol. 23, No. 5

The obtained values of the vertical dispersion of the turbulent velocity agree fairly well with experimental data, whereas the computed values of the radial dispersion of the turbulent velocity exceed the experimental analog within the mixing layer ( 0 < z/ zi ≤ 1 ), which may be related to the influence of vertical walls of the laboratory setup [1], whose effect probably suppressed in the quasi-steady state the horizontal motion and reduced the level of the horizontal dispersion of turbulent velocity. In the problem under consideration, a simple model of gradient diffusion for turbulent flows not only describes correctly the features of the distributions of the vertical and horizontal dispersions of the turbulent velocity but also reflects satisfactorily their anisotropic character. Appendix

In the Appendix, the balance equations for the quantities E , ε and θ 2 in the cylindrical coordinate system are presented: K  ∂ E  ∂E 1 ∂ ∂ 1 ∂   + ⋅ ( rUE ) + Fr (WE ) = ⋅ r  Re−1 + m  + ∂t r ∂ r ∂z σ E  ∂ r  r ∂ r   +

∂   −1 K m  ∂ E  r  Re + + P +G −ε,  ∂ z   σ E  ∂ z 

(A.1)

K  ∂ε  ∂ε 1 ∂ ∂ 1 ∂   + ⋅ ( rU ε ) + Fr (W ε ) −  r  Re−1 + m   + ∂t r ∂ r ∂z r ∂ r   σ ε  ∂ r  −

where

∂   −1 K m  ∂ε  r  Re +   = −Ξ, σ E  ∂ z  ∂ z  

  ∂U ∂U 2  ∂U  ∂W 2  ∂W ∂W  P =  2Km − E +  2 K m Fr − E  Fr + Km  + Fr  ∂r 3  ∂r  ∂z 3  ∂ z ∂r    ∂z

(A.2) 2

is

the generation of the turbulence kinetic energy by the velocity shear, G = β g wθ is the generation of the turbulence kinetic energy by the work of the fluctuating buoyancy force, σ E and σ ε are the turbulent Prandtl numbers (σ Е = 1, σ ε = 1.3). The function Ξ=

 1 ε2 1 ε ∂U ∂W  Ξ0 + Ξ1 ( P + G ) + Ξ 2  Fr −1 ⋅ uθ + wθ  2 E 2 E z ∂ ∂z  

(A.3)

includes in comparison with its standard representation [19] in the form of a sum of the terms linear with respect to the mean velocity shear ∂U i /∂ x j of the buoyancy vector β gi the quadratic term with coefficient Ξ 2 to account for the effects of the interaction between the mean velocity shear and the thermal field anisotropy including the night and transitional periods of the planetary boundary layer [20]. The action of the velocity and buoyancy mean shear is taken into account by the first-order terms in (A.3). The coefficients in (A.3) have the numerical values selected at the solution of various problems of the planetary boundary layer [20]: Ξ 0 = 3.8, Ξ1 = 2.4, Ξ3 = 0.3. The balance equation for the dispersion of temperature fluctuations has the following closed form: 691

A.F. Kurbatskii and L.I. Kurbatskaya

∂θ 2 1 ∂ ∂ 1 ∂   E  ∂θ 2  + ⋅ rU θ 2 + Fr W θ 2 = ⋅ r cθ 2 θ 2  + ∂t ∂z r ∂r r ∂r   ε  ∂r   

(

+

)

( )

E 2  ∂θ 2  ∂   ∂T ∂T 1 ε 2 − 2 wθ − ⋅ θ ,  r cθ 2 w   − 2uθ ε ∂z   ∂r ∂z R E  ∂ z  

(A.4)

here R = τ θ /τ = (1/2 ) θ 2 /εθ / ( E /ε ) = 0.6 [11], cθ 2 = 0.11 [10]. References 1. J. Lu, P. Arya, W.H. Snyder, and R.E. Lawson, A laboratory study of the urban heat island in a calm and stably stratified environment. Part I, II, J. Appl. Meteor., 1997, Vol. 36, No. 10, P. 1377−1402. 2. G.E. Willis and J.W. Deardorff, A laboratory model of the unstable planetary boundary layer, J. Atmos. Sci., 1974, Vol. 31, P. 1297−1307. 3. A.F. Kurbatskiy and L.I. Kurbatskaya, Three-parameter model of turbulence for the atmospheric boundary layer over an urbanized surface, Izvestiya. Atmospheric and Oceanic Phys., 2006, Vol. 42, No. 4, P. 439−455. 4. D.W. Byun and S.P.S. Araya, A two-dimensional mesoscale numerical model an urban mixed layer. I. Model formulation, surface energy budget, and mixed layer dynamic, Atmospheric Environment, 1990, Vol. 24A, No. 4, P. 829−844. 5. R.A. Pielke, Mesoscale meteorological modeling. 2nd edition, Academic Press, New York, 2002. 6. A.F. Kurbatsky, Modeling of Non-local Turbulent Transfer of Momentum and Heat, Nauka, Novosibirsk, 1988. 7. C.G. Speziale, Modeling of turbulent transport equations, in: T.B. Gatski, M.Y. Hussaini, J.L. Lumley (Eds.), Simulation and Modeling of Turbulent Flows, Oxford University Press, Oxford, 1996. 8. Y. Cheng, V.M. Canuto, and A.M. Howard, An improved model for the turbulent PBL, J. Atmos. Sci., 2002, Vol. 59, No. 9, P. 1550−1565. 9. T.P. Sommer and R.M.C. So, On the modeling of homogeneous turbulence in a stably stratified flow, Phys. Fluids, 1995, Vol. 7, No. 11, P. 2766−2777. 10. L.H. Jin, R.M.C. So, and T.B. Gatski, Equilibrium states of turbulent homogeneous buoyant flows, J. Fluid Mech., 2003, Vol. 482, P. 207−233. 11. C. Beguier, I. Dekeyser, and B.E. Launder, Ratio of scalar and velocity dissipation time scales in shear flow temperature, Phys. Fluids, 1978, Vol. 21, No. 3, P. 307−310. 12. Y.A. Panofsky, Y. Tennekes, D.Y. Lenshow, and J.C. Wyngaard, The characteristics of turbulent velocity components in the surface layer under convective conditions, Boundary Layer Meteor., 1977, Vol. 11, P. 353−361. 13. J.A. Businger, J.C. Wyngaard, Y. Izumi, and E.F. Bradley, Flux profile relationship in the atmospheric surface layer, J. Atmos. Sci., 1971, Vol. 28, P. 181−189. 14. J.-F. Louis, A parametric model of vertical eddy fluxes in the atmosphere, Boundary-Layer Meteor., 1979, Vol. 17, No. 2, P. 187−202. 15. A.S. Monin and A.M. Yaglom, Statistical Fluid Mechanics. Vol. 1: Mechanics of Turbulence, MIT Press, Cambridge, Mass., 1971. 16. J.C. Andre, G. de Moor, F. Laccarere, G. Therry, and R. du Vachat, Modeling the 24-hour evolution of the mean and turbulent structures of the planetary boundary layer, J. Atmos. Sci. 1978. Vol. 35. P. 1861−1883. 17. P.G. Duynkerke, Application of the e-ε turbulence closure model to the neutral and stable atmospheric boundary layer, J. Atmos. Sci., 1988, Vol. 45, No. 5, P. 865−879. 18. P. Roache, Computational Fluid Dynamics, Hermosa, Albuquerque, New Mexico, 1976. 19. O. Zeman and J.L. Lumley, Buoyancy effects in entraining turbulent boundary layers: a second-order closure study, in: F. Durst, B.E. Launder, F.W. Schmidt, and J.H. Whitelaw (Eds.), Turbulent Shear Flows I, Springer, Berlin, Heidelberg, 1979, P. 295−306. 20. A.A.A. Andren, TKE dissipation model for the atmospheric boundary layer, Boundary Layer Meteor., 1992, Vol. 56, P. 207−221.

692