Numerical and experimental investigation of unsteady natural convection in a non-uniformly heated vertical open-ended channel

International Journal of Thermal Sciences 99 (2016) 9e25

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Numerical and experimental investigation of unsteady natural convection in a non-uniformly heated vertical open-ended channel
ne
zo b, c, G.H. Yeoh a, O.A. Tkachenko a, V. Timchenko a, *, S. Giroux-Julien b, C. Me J.A. Reizes a, E. Sanvicente b, M. Fossa d a

School of Mechanical and Manufacturing Engineering, University of New South Wales, Sydney 2052, Australia Centre d'Energ
etique et de Thermique de Lyon (CETHIL UMR 5008), Universit
e de Lyon, Villeurbanne, France Chaire INSA/EDF “habitats et Innovations
energ
etiques”, Villeurbanne, France d Dime, Dept. Mech. Eng., University of Genova, Italy b c

a r t i c l e i n f o

a b s t r a c t

Article history: Received 3 December 2014 Received in revised form 14 July 2015 Accepted 15 July 2015 Available online xxx

Experimental and three-dimensional numerical investigations have been undertaken to improve the understanding of the flow and heat transfer phenomena in narrow, open-ended channels similar to those formed by the double skin façades with PV panels. Non-uniform heating configurations in which heat sources alternated with unheated zones on both walls were studied to simulate opaque PV arrays and glazed panes/windows of the building. Heat transfer and flow measurements were obtained for periodicity 1/15 of the overall height of the heated and unheated zones with a heat input of 220 W/m2. Largeeddy simulations were utilized to capture the flow behaviour, heat transfer rates and possible transition to turbulence. A separation of the thermal plumes from the heated walls of the channel has been observed both experimentally and numerically. In the case of the non-uniformly, both sides heated channel these plumes move across the channel and turn into travelling waves leading to significant mixing thereby allowing increases in the heat transfer rates comparing with one side uniformly heated channel. In addition, the staggered position of the heated and unheated zones on the opposite walls of the channel results in a recirculation zone occurred in the unheated area of one of the walls at the channel entrance. This separation contributes to the periodicity of flow within the channel and the lower average temperatures compared with uniformly heated channel. The improved mixing as the result of the staggered heating on both sides of the channel was confirmed by the parameter termed degree of mixing which was significantly higher in the channel heated from both sides. The results of this study reveal that staggered arrangement of PV panels will produce higher mass flow rate of air entrained into the channel which together with enhanced mixing would further reduce the PV surface temperatures and prolong their life span. © 2015 Elsevier Masson SAS. All rights reserved.

Keywords: Unsteady natural convection Non-uniform heating Large eddy simulation (LES) Building integrated photovoltaic (BIPV) systems

1. Introduction As the operating temperature of photovoltaic (PV) cells increases, their electrical energy output is reduced, whilst their ageing process is accelerated [1e3]. It follows that one of the many significant challenges in the utilization of building integrated photovoltaic (BIPV) systems is controlling and limiting the operating temperatures experienced in PV modules [1,2]. One possible method of cooling the modules is to use natural convection

* Corresponding author. Tel.: þ61 29385 4148. E-mail address: [email protected] (V. Timchenko). http://dx.doi.org/10.1016/j.ijthermalsci.2015.07.029 1290-0729/© 2015 Elsevier Masson SAS. All rights reserved.

generated in the double façade of a building formed by the PV arrays and the outward skin of the building. One embodiment of such a façade is shown in Fig. 1. The narrow channel between them the PV array and the building allows air to circulate by natural convection since the PV arrays are heated by solar radiation. Passive cooling of the PV modules is therefore obtained as noted by Sandberg and Moshfegh [4]. This has been successfully used as a cost-effective means of reducing the operating temperatures [5,6] of PV arrays. Brinkworth et al. [5] have shown that natural convection, particularly in vertical channels, is capable of generating flow rates comparable with those obtained with forced ventilation.

10

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

Nomenclature

Roman symbols A disturbances amplitude (ms1) b width of heated walls (m) Cp heat capacity (J/K) D distance between channel walls (m) g gravitational acceleration (m/s2) h heat transfer coefficient (W/m2 K) H height of the channel (m) F view factor () k thermal conductivity (W/m K) l thickness of insulation (m) L width of channel wall (m) Nu Nusselt number () p pressure (Pa) Pr Prandtl number () Q power (W) q heat flux (Wm2) Ra Rayleigh number ()

t T u Vy

time (s) temperature (K) velocity (m/s) vertical velocity component (m/s)

Greek symbols ε emissivity () m dynamic viscosity (Pa s) r density (kg/m3) s StefaneBoltzmann constant (W/m2 K4) ~ s stress tensor (Pa) Subscripts a air conv convective f foil i, j, k indices in the x, y and z directions respectively ref reference s insulation out outside the channel

Fig. 1. Double façade, with PV modules separated by glass partitions.

An example of an actual double façade is shown in Fig. 1(a). In this configuration, the outer skin of the building is composed of alternating PV modules and glazed panels over its entire height. Similarly, the inner skin has alternating irradiated opaque and shaded portions of wall or windows. The façade in Fig. 1(a) is depicted schematically in Fig. 1(b). Here, the PV modules are represented as localized heat sources and the glazed panes shown as unheated zones on the outer skin, similarly the building envelope has heated and unheated zones which out of step with those on the PV array. Considerable research has been focused on uniformly-heated open-ended channels to gain an understanding of processes within BIPV systems [4,7e10]. Configurations with adjustable channel inclination angles, aspect ratios, and uniform heat flux on boundaries, have been extensively studied. Arrangements such as the staggered heating configuration demonstrated schematically in Fig. 1(c) has been shown experimentally by Fossa et al. [11] to achieve better cooling of the heated zones inside the channel. For a given channel heat flux and channel aspect ratio, the discrete heating arrangement was shown to promote better heat transfer, consequently reducing the maximum wall temperature compared with a uniformly heated wall [11e15]. For non-uniform heating, the

wall temperature distribution has been found to vary periodically in the vertical direction with the alternation of heated and unheated zones [11,12,14,16]. The convective heat transfer coefficient thereby rises sharply in the heated zones, from a maximum value at the leading edge to a minimum at the trailing edge, and decreases on succeeding heated zones. Moreover, optimal arrangements of heat sources along the walls appears to exist, so that criteria has been given by Silva et al. [15,17] on the location and distribution of heat sources that maximizes the global thermal performance of the system. The exact reasons for the improved heat transfer rates by the introduction of alternating heated and unheated zones on both sides of the channel is not fully understood and needs to be investigated. Therefore, to improve understanding of physical processes occurring in such configurations numerical and experimental investigations have been undertaken to determine the effects of the non-uniformity of the wall heat flux distribution on the natural convective flow and heat transfer within a vertical rectangular open ended channel. Unfortunately, since, as may be seen in Fig. 1(b), the heated and unheated zones on the building envelope will change with the geographic location of the building, the seasons and of course the time of day, the study would involve a

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

complex, transient convection problem equally difficult to model numerically as it would be experimentally. A configuration which does not change in time had to be employed so as reduce the computing requirements. Further, in order to understand the reasons for an improved performance found by Fossa and his coworkers [11,12,14,16], the same configuration, as that employed by these researchers, namely that shown in Fig. 1(c), was adopted in this paper. Experiments were performed in an open-ended channel built at CETHIL research centre in Lyon, France, with the surface temperatures and velocity measurements at different locations of the channel being obtained. In parallel, a three-dimensional computational model, using large eddy simulation (LES) has been developed to simulate the fluid flow and heat transfer in the open-ended channel to better understand the flow behaviour and possible transition from the laminar to turbulent regimes in the complex non-uniformly heating configuration. In a preliminary two-dimensional numerical study it has been shown that when a natural convection apparatus of the type shown in Fig. 1(c) is located in a laboratory, complex interactions between the flow in the room and that in the channel occur. Indeed, the large-scale flow structures ingested into the channel contribute greatly to the unsteady flow within the channel. The effect of disturbances introduced by the external environment is modelled by small random fluctuations introduced at the channel inlet, thereby leading to better agreement with experimental results. Within the channel, thermal radiation contributes significantly to the heat transfer rates, even at low emissivities and wall temperatures [8]. An approximate but, numerically efficient model of wall-to-wall radiation is developed and solved alongside the equations for fluid flow and convective heat transfer. In order to determine the reasons for the improved performance when discrete heat sources as shown in Fig. 1(c) are used, the main aspects investigated in this study concern the unsteady character of the flow and its structure and the characterization of the transition stage. Additionally, the enhancement of the mixing in the flow, which might be expected in such a configuration, is also studied. These questions are addressed through the examination of the mean flow quantities of the fluid flow obtained from experimentally measured particle image velocimetry (PIV) data and numerical results. Typical temporal and spatial temperature distribution and flow dynamics in terms of the frequency spectra of velocity and temperature fluctuations are also obtained experimentally and numerically. Finally, the changes in the mean cross-sectional temperature distribution as a function of height are used to illustrate the improved mixing as the result of the staggered heating on both sides of the channel. 2. Experimental details An open-ended channel built at CETHIL research centre in Lyon, France, shown in Fig. 2 was used in this study. The channel consisted of two 1.5 m  0.7 m parallel walls 0.1 m apart with, as may be seen in Fig. 2(b), 30 chamfers cut at the channel inlet. The lateral sides were closed by plexiglass sheets. The vertical walls were insulated with 12 cm thick polyurethane blocks with thermal conductivity of 0.027 W/(m,K). In order to mimic the heat transmitted from the solar cells, the channel walls were covered on the inside with 15 independently controlled stainless steel foil heaters, each 10 cm wide and 50 mm thick. The emissivity of each foil was 0.092 and the thermal conductivity 13 W/(m,K). As may be seen in Fig. 2(b) the channel walls were vertical with the channel inlet placed 0.75 m above the floor and the channel outlet 0.75 m below an artificial ceiling.

11

Standard k-type thermocouples were used to measure the temperatures and a Particle Image Velocimetry (PIV) system with appropriate software was employed to obtain the values of velocity at different locations of the channel. The 205 thermocouples employed are divided in 5 groups. Each group of 41 thermocouples has its own cold junction whose temperature was measured by two platinum resistance thermometers (type Pt100). 75 thermocouples were placed every 2 cm on the vertical centreline of the wide sides (z/L ¼ 0.5) in contact with the rear of the heating foils. The first thermocouple was placed 1 cm from the leading edge of each heater as shown on the expanded portion in Fig. 3(a). This type of arrangement gives five measurements of temperature per heating band and limits the risk of disturbing the flow. Wall temperatures were also measured with 16 thermocouples fitted, as shown in Fig. 3(a), on one of the walls between the central plane and the lateral edge of the channel and at five different vertical levels in order to monitor the horizontal variation of the heated surface temperature and evaluate the lateral conduction losses in the z-direction. A total of 30 thermocouples, one per heating strip, were also embedded 4 cm deep in the insulating material of each wall, as may be seen in Fig. 3(a), on the line z/ L ¼ 0.5, so as to allow an estimate of the heat losses through the insulation. Particles of Di-Ethyl-Hexyl-Sebacat (DEHS) with diameters less than 1 mm and mean density of r ¼ 912 kg/m3 were injected into the channel for the measurement of velocities. The acquisition of the PIV images was done utilizing a charge-coupled device (CCD) camera with a resolution of 10-bits and a frame grabber whose acquisition frequency was 11 Hz. Each image consisted of an array of 1024  1024 pixels, with a square pixel of 7.4 mm. Each field of view represents a section of the flow of 200  200 mm2, with the camera located at a distance of 0.9 m from the object plane. The experimental measurement techniques as well as experimental methodology are detailed by Sanvicente et al. [10,18]. 3. Computational model 3.1. Conservation equations The in-house computer code developed and validated in Lau et al. [9,19,20] has been used to simulate three-dimensional natural convective flows in a non-uniformly heated open ended channel such as that shown in Fig. 1(c) was adopted. Since Lau et al. [9] have shown that Reynolds-averaged NaviereStokes (RANS) methods do not yield results which are in reasonable agreement with experimental data, whereas Large Eddie Simulations (LES) do [20], the LES approach with the Vreman subgrid model [21] suggested by Lau et al. [9] was adopted here. In natural convective flows, the changes in density may be large and are taken into account by means of Favre-averaging. For any arbitrary variable 4(x,t), the Favre~ tÞ is obtained by averaged variable fðx;

~ tÞ ¼ rðx; tÞ$fðx; tÞ ; fðx; rðx; tÞ

(1)

in which the curly overbar (~) represents the Favre-averaging operator and the overbar (e) represents a mean quantity. Since the Favre averaged variable defined in Eq. (1) contains two other averages, the equations of continuity the conservation equation need to be averaged in the same way. The time over which these averages are taken has to be sufficiently short so as to allow the averaged parameters to remain functions of time. This operation is often called “filtering” as the high frequencies are cut from the solution. It follows that whilst the large, low frequency structures are simulated, the high frequency small structures, have to be

12

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

Fig. 2. Schematic of experimental apparatus at Lyon, CETHIL UMR 5008, France (a) Photograph of the experimental cell; (b) Schematic of geometry and materials.

Fig. 3. Measurement arrangements (a) Locations of the thermocouples; (b) Laser and PIV camera arrangement.

modelled. After averaging, the conservation of mass, momentum and energy equations in a three-dimensional Cartesian coordinate frame are

vr vruej þ ¼ 0; vt vxj

(2)

 vr$uej vruei uej sij vMij  vp vf þ ¼ þ þ þ r  rref gj ; vxj vxi vt vxj vxi

(3)

vT~ vT~ v r$Cp $ þ rCp uei ¼ vt vxi vxi

m$Cp vT~ Pr vxi

! þ

vtu;T ; vxi

(4)

in which r and T are density and temperature respectively, while rref is the reference density at the reference temperature Tref, ui is the velocity component in the ith direction, p is the dynamic pressure, Pr is the molecular Prandtl number, m is the dynamic viscosity and sij is the stress tensor which, using the Stokes' hypothesis can be written as

0

1 e e v u v u j B C 2 g vuei g@ iþ fij ¼ mðTÞ s d ; A  mðTÞ 3 vxj vxi vxj ij

(5)

in which dij is the Kronecker delta. The dynamic viscosity which appears in Eqs. (4) and (5) as a function of temperature is calculated using Sutherland's law [22].

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

The tensor Mij in Eq. (3) is taken as evðr$tui ;uj Þ=vxj , in which the subgrid stress tensor tui,uj for the Favre-averaged momentum field ei uej. is given by ðug i uj Þ  u 3.1.1. Sub-grid scale model (SGS) The terms Mij and tu,T in Eqs. (3) and (4) represent the sub-grid terms obtained when the Favre-averaging operation is applied to the conservation equations. These terms need to be appropriately modelled to account for the effect of small-scale eddies on the large-scale structures. The thermal flux vector due to subgrid-scale (SGS) motion, tu,T, is correlated by the SGS Prandtl number Prsgs in accordance with

tu;T

(6)

in which Prsgs is 0.4 for air [23]. The unresolved turbulent momentum stress tensor Mij is modelled as

Mij ¼ tui ;uj

! 1 1f e  $tuk ;uk dij ¼ 2msgs Sij  Skk dij ; 3 3

(7)

in which msgs needs to be modelled and Seij ¼ 1=2ðvuei =vxj þ vuej =vxi Þ. In the current investigation, the SGS model developed by Vreman [21] is adopted. It has been found by Lau et al. [19] to reasonably accurately predict flows in natural convection and in particular predict transition to turbulence quite accurately. Follow Lau et al. [19], the subgrid viscosity in Eq. (7) is given by

msgs

sffiffiffiffiffiffiffiffiffiffiffi Bb ¼ rCv ; aij aij

(8)

vuej ; b ¼ D2m $ami amj ; vxi ij

Bb ¼ b11 $b22  b212 þ b11 b33  b213 þ b22 b33  b223 ;

3.3. Boundary conditions

  fij þ pdij ¼ 0 V s

(9)

(10)

Dm is filter width and, as suggested by Vreman [21], Cv is set to 0.07225. 3.2. Computational domain Fig. 4 is a schematic illustration of the computational domain and the mesh adopted in this study. This geometry is intended to simulate the test section used by Vareilles [13,24] and GirouxeJulien et al. [13] and in which the experimental results presented in this paper were obtained. As shown in Fig. 4(a), regions outside the channel were included in the computational domain below the channel inlet and above the channel outlet. Including these features in the numerical model was necessary as the flow is significantly affected by the presence of the floor and the false ceiling [7]. These zones were open on all sides, with the upper region was delimited by the horizontal solid wall 0.75 m above the outlet of the channel. However, Lau et al. [7] found that only half the region below the channel needed to be included for satisfactory numerical solutions, so that only half the physical region was included below the channel. For transient numerical calculation, the computational domain has been meshed using a rectangular mesh. The mesh employed in the present investigation is the same as that used by Lau et al. [7,25] who had performed a thorough mesh convergence study for the

(11)

was imposed at the openings, thereby allowing natural entrainment of ambient fluid. A no-slip boundary condition is enforced on all solid walls. As demonstrated by Sanvicente et al. [10], natural convective flows in open-ended channels are highly sensitive to ambient conditions which lead to external disturbances being introduced through the inlet section and possibly amplified in the channel. This was supported by a preliminary two-dimensional numerical study. Unfortunately, due to the very large computational demands it was not possible to extend this study to three dimensions, or for the same reason include it as part of the present investigation. Thus, the disturbance of the inflow needs to be modelled. In this work, the same process as developed by Lau et al. [28] was employed to simulate the effect of external disturbances on the flow. A fluctuating velocity boundary condition,

unþ1 ¼ unijk þ Aε; ijk

in which

aij ¼

open ended study for an open-ended channel configuration with uniform heating studied by Miyamoto et al. [26]. With the internal channel dimensions of 0.1 m by 1.54 m by 0.7 m, the mesh consists of 41  330  79 elements. It is uniform in the x and y directions with step sizes of 2 mm and 5 mm, respectively, and non-uniform in z direction with mesh size increasing from 2.5 mm at the boundaries to 12 mm at the centre of the channel. The chamfers at the channel inlet, as shown in Fig. 4(b), were approximated by rectangular volumes.

For the fluid flow, a stress-free boundary condition [27].

!

m vT~ g ei T~ z sgs $ ; ¼ r ðu i TÞ  u Prsgs vxi

13

1  ε  1

(12)

was imposed on all mesh cells on the lower boundary of the inlet in Fig. 4(a). Here, unijk is the horizontal velocity perpendicular to the face of any mesh cell (i, j, k) on the boundary of the inlet, calculate at the end of the nth time step, A is the amplitude of the disturbance and ε is a random number. The boundary condition for the velocity at the mesh cell (i, j, k) at the (n þ 1)th time step, is then unþ1 . ijk The problem is fully specified once the thermal boundary conditions have been established. The electrical power input per unit area to each of the foils, qelec, is given. The convective flux qin was estimated by subtracting the radiation heat flux, qrad, together with the heat flux lost to the atmosphere through the insulation, qcond from qelec. Radiation heat transfer had to be taken into account as it has been shown to be significant in natural convective flows in open-ended channels [25]. It was assumed that the electric heat flux is uniform all over the heating foils. 3.3.1. Conduction through insulation A one-dimensional model of the heat transfer is employed to estimate the heat loss through the insulation namely,

qcond ¼

Tf  T∞ l ks

þ h1

;

(13)

out

in which l is the thickness of the insulation and ks its thermal conductivity, hout is the coefficient of convective heat transfer from the insulation to the ambient, and the subscripts f and ∞ refer to the foil and to the ambient respectively. The heat transfer coefficient hout has been assumed to have a value of 5 Wm2 K1, which had been used previously by Lau et al. [29].

14

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

Fig. 4. Computational domain: (a) vertical projection, (b) 3D isometric view, (c) slices of 3D mesh.

Fig. 5. Wall-to-wall radiation view-factors diagram. (a) View-factors within the open-ended channel; (b) Notation used in the view factor calculation.

Heat losses to the ambient were calculated for the insulated walls (yez planes) as well as the lateral, plexiglass walls (xey planes). The thickness of the plexiglass sheeting was 5 mm, and its thermal conductivity was 0.2 W/(m K). For polyurethane insulation, thickness is 12 cm and thermal conductivity is 0.027 W/(m K). The ambient temperature, T∞, was set to 18.4  C, the mean temperature in the laboratory.

The radiosity is calculated for each heater; thus, a linear system of equations is solved each time when the radiative heat transfer is calculated:

sTj4 

k  X m¼1

3.3.2. Radiation between walls Similarly to Fossa et al. [11], a simplified calculation for the wallto-wall radiation has been developed to calculate the radiation losses through the radiosity method (Fig. 5). In this study, radiative heat fluxes are calculated from view factors and temperatures which change with time as the flow develops.

 k   q X 1  εm j 4 Fjm sTm ¼  Fjm qm ; εj m¼1 εm

(14)

in which j and m are indices of heaters on the walls (see), qj is the radiative heat flux emitted by jth heater, Fjm is the view factor of pair of heaters, s is the StefaneBoltzmann constant, and εj is emissivity of jth heater. Radiation from and to the lateral sides of the channel, as well as from and to the outside was neglected in the 4 was used. model, and simplified formula qm ¼ εm =1  εm sTm

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

View factors for each pair of heaters are evaluated from a reference parallel plates equation [30] as function of their size and location:

P2 P2 F12 ¼

l¼1

k¼1

P2 P2 j¼1

i¼1 ð1Þ

ðiþjþkþlÞ

ðx2  x1 Þðy2  y1 Þ

  G xi ; yj ; hk ; xl

;

(15)

in which G(xi,yj,hk,xl) is the integration function of a geometry of parallel plates given by.

9 8 0 > > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi = < 1 B y  h 2 G¼ @ðy  hÞ ðx  xÞ þ z2 tan1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ > 2p ; : ðx  xÞ2 þ z2 > 9 8 > > qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi = < xx þ ðx  xÞ ðy  hÞ2 þ z2 tan1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  > ; : ðy  hÞ2 þ z2 > 1  z2  C  ln ðx  xÞ2 þ ðy  hÞ2 þ z2 A 2

(16)

3.4. Numerical methodology The conservation equations for variable-property Newtonian flows were solved numerically using an explicit two-step predictorcorrector method. This algorithm has been demonstrated to effectively handle the changes of density in natural convective flow. Details of this numerical scheme including its full derivation and implementation can be found in Lau et al. [9] and for brevity will not be included here. The solution was advanced in time using a second-order AdamsBashforth time integration scheme for the predictor stage and a second-order quasi Crank-Nicolson integration for the corrector stage. The pressure correction, which takes into account the compressibility of air and which was incorporated in both the predictor and corrector stages, involves the solution of pressure correction Poisson equations. The equations are discretised through the finite volume formulation on a collocated grid. A fourth-order central differencing scheme was employed to approximate the convective fluxes while a second-order central differencing scheme was adopted for all other spatial derivatives. To accurately simulate the flow, it is important to recognize the range of length and time scales that are required to be resolved for valid large eddy simulation calculations. An adequate grid needs to be developed. Lau et al. [29] performed a grid convergence study on a very similar problem and, as stated above, the same grid was used here. Numerical stability was ensured by employing that the time step used resulted in a maximum CouranteFriedricheLewis criterion of 0.35. The convergence in each case for the pressure correction steps at each time-step was less than 105. Simulations were continued for approximately 160,000 time steps in order to have sufficient data to satisfactorily capture the turbulent fluctuations and allow reliable statistics to be obtained. 4. Results and discussion The effects of the non-uniformity of the wall heat flux distribution on the temperature and the velocity fields in the open ended channel, as well as on the heat transfer efficiency, are explored both numerically and experimentally for the Case V1 shown in Fig. 6. The aspect ratio of the channel was D/H ¼ 1/15 and each of the heated

15

and unheated regions had a height of a/H ¼ 1/15 with an average electrical heat input to each of the heated zones of qelec ¼ 220 W/ m2. This configuration was compared with the case of one uniformly heated wall with qelec ¼ 220 W/m2 wall, designated as VREF [18] in Fig. 6. It should be noted that the total electric power input to the channel was the same in both cases as is shown in the last column of Table 1. However, since the power supplied to sides S1 and S2 in Fig. 6, were not the same in the two cases, the electric power supplied to each to the two sides are also shown in Table 1. Since the contribution of heat input on S1 and S2 for the non-uniformly heating configuration is summarized in Table 1. As mentioned above, a preliminary two-dimensional numerical study of a channel which was confined in a closed laboratory such as that was used for the experiments, showed that external disturbance significantly affected the flow. In particular, the large flow structures generated by the hot air emanating from the apparatus and which circulated around the room were from time ‘sucked back’ into the inlet of the channel and completely altered the flow patterns. As a consequence, following Lau et al. [28], the process depicted by Eq. (12) has been adopted. As a first step therefore, the value of A in Eq. (12) needed to be established so as to allow a reasonable comparison with the experimental results. Lau et al. [28] had used A ¼ 0.02 ms1 and obtained results which compared well with experimental data for the VREF case. Although the same experimental apparatus was used, it was not clear that the same value of A would be appropriate when the channel was heated non-uniformly and asymmetrically on both sides. Computational results were obtained with imposed disturbances A ¼ 0, 0.02 and 0.05 ms1. Since the temperatures on the heated walls depend on the details of the flow as well as the radiation and heat loss calculations, the values of the time averaged temperature calculated numerically on the centreline of each of the heated walls were compared with those obtained experimentally to determine the value of A which should be used. As may be seen in Fig. 7(a), the results with no disturbance from the surroundings, A ¼ 0 ms1, led to a significant overestimate of the time-averaged temperature of the walls in comparison with the experimental data. Similarly, as is shown in Fig. 7(b), the in numerical results with no disturbance on S2 wall, A ¼ 0 ms1, also appreciably overestimated the magnitude of the time averaged temperature on the centreline relative to the experimental data. When A ¼ 0.05 ms1 was used, as may be seen in Fig. 7(a), the temperature on the lower part of the S1 wall lower than the experimental results and those calculated with A ¼ 0.02 ms1. However, above the first heater on S1, both numerically calculated temperature distributions with As0 ms1 almost coincided and were in good agreement with the experimental data. Surprisingly the numerical results obtained As0 on S2, coincided over the whole height of the heater and, except for the region close to the exit, were in good agreement with experiment. As the result, in concord with Lau et al. [28], unless otherwise stated, the results presented in this paper are the data obtained with A ¼ 0.02 ms1 as the disturbance. 4.1. Heat transfer characteristics at the walls For the non-uniformly heated configuration, Case V1, as would have been expected and may be seen in Fig. 7, the time-averaged wall temperature distribution has a spatial periodicity corresponding to the oscillation between heated and unheated zones. The surface temperature rises sharply in the heated zones, from the minimum value at the leading edge to the maximum one at the trailing edge, and decreases along the unheated zones. Similar

16

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

Fig. 6. Configurations studied.

Table 1 Heat flux and total power values for the cases studied.

Experimental Numerical

2

Case

qelec (W/m )

Qelec

V1 VREF V1

220 220 220

123 231 123

S1

(W)

Qelec 108 0 108

S2

(W)

Qelec

S1þS2

(W)

231 231 231

temperature structures were observed for non-uniformly heated channels by other researchers [11e14]. Case V1 yielded approximately 12  C lower values of the experimentally and numerically measured time averaged temperatures on the centreline of wall S1, than those obtained experimentally for uniformly heated configuration, Case VREF, indicated by dotted line in Fig. 7(a). Since there is one fewer heater and one more unheated zone on side S2 than side S1, the time averaged temperatures at the centreline of S2 shown on Fig. 7(b) for case V1 are 1.5  C lower than those on wall S1. It is interesting to note in Fig. 7(a) that the experimentally obtained temperature rise on the lowest heater on side S1, follows almost exactly the experimentally determined temperature distribution for the VREF configuration over the same height. However, as may also be seen in Fig. 7(a) none of the numerically evaluated temperature distributions on the surface follow the experimentally found temperature distributions. Experimentally there is gradual temperature rise with height on the first heater on S1. The numerically predicted temperature on the same heater, with A ¼ 0.02 ms1, has a very sharp rise follow by a significant fall in temperature which is followed by another, albeit smaller temperature rise. As mentioned above, the chamfers were introduced on the experimental apparatus to prevent separation at the inlet. Since the in-house computer code

employed orthogonal meshes, in computational model the chamfers were approximated by a series of steps shown in Fig. 4(b) which possibly led to larger separations on the wall S1 and resulted in discrepancy between numerical and experimental results near the inlet. At the top of the heater the experimental of the temperature value is almost the same as that generated by the simulation. The temperature curves generated with A ¼ 0 and A ¼ 0.02 ms1, in Fig. 7(a) follow similar patterns on either side of the A ¼ 0.02 ms1 temperature curve. At the top of the S1 wall the numerically obtained data over predicts the experimentally determined temperature measurements, but this is due to losses through the ends of the apparatus not modelled here. However, the agreement between the numerical simulations and the experimental measurements over the remainder of the S1 wall is excellent. Whilst the agreement between the temperature distribution from the numerical calculations and the experiment over the first unheated zone on S2 is very good, there is a discrepancy over the heater immediately above that zone. This discrepancy diminishes with height with almost exact agreement for 0.27  (y/H)  0.53. For (y/H) > 0.53, the numerically calculated temperature drop in the unheated zones on S2 is lower than the experimental result and, as a consequence, the temperature at the trailing edge of the heaters in the simulations is higher that determined by experiment. As shown in Table 1, cases V1 and VREF have the same power input. But, since the power supplied to VREF is dissipated to the fluid through the S1 wall only, the temperatures on that wall, as may be seen in Fig. 7(a), are substantially higher than either S1 or S2 in the case V1 which has power supplied to both walls.

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

A ¼ 0.02 ms1;

A ¼ 0.05 ms1; Experimental:

17

V1; j j : j j VREF.

Fig. 7. Time average temperature on the centreline of the heated walls. Numerical V1:

A ¼ 0;

Since in both the VREF [28] and the V1 cases there is good agreement between the numerical and experimental timeaveraged temperature results, a similarly close agreement between experimental and numerical data may be expected between the heat transfer coefficients, hconv, defined as

4.2. Time-averaged flow behaviour in the non-uniformly heated channel

hconv ¼ qconv =ðTwall  Tinlet Þ;

(17)

in which

  .    dx1  qconv ¼ ka ðTÞ Tf  Tadj

(18)

with ka being the thermal conductivity of air, dx1 the mesh size in the x direction at the wall and the subscript adj the mesh in the fluid immediately adjacent to a heated wall. The good agreement between the numerical and experimental data is shown in Fig. 8. As the air in the channel is gradually heated, the capacity of the fluid to carry heat away from the solid surface is reduced. In the case VREF with uniform heating the heat transfer coefficient, as expected, gradually decreases with height of the channel. On contrary, in the case V1, due to the disruption of the thermal boundary layer, the heat transfer coefficient decreases along each heated zone but, on the average it remains higher than in the uniformly heated on one side case VREF which also contributes the a lower temperatures along the wall.

Streamwise velocity profiles at different heights of the channel for the V1 configuration are shown in Fig. 9(a). Similarly to the temperature distributions shown Fig. 7, numerically produced time-average velocity profiles are in reasonable agreement with experimental data, particularly in the upper part of the channel. The time-averaged velocity distribution near the inlet section at y/H ¼ 0.18 has an almost flat shape across the channel width. The presence of localized small peaks close to wall S1 is linked to the distribution of the heated zones. As the flow moves upwards, the shape of the velocity profiles evolves with large maximums in the vicinity of each wall. In all cases, the higher the cross section is taken, the larger the magnitude of the velocity peaks relative to the velocity near the centre of the passage. Since the S1 wall receives 14% more electric power than the S2 wall, as may be seen in Fig. 9(a) at all levels, the maximum velocity near the S1 wall is higher than that near the S2 wall. It should be noted that the axial velocity distributions presented in Fig. 9(a) are for the centre plane, (z/L) ¼ 0.5 as these were the only experimental velocity profiles. However, these distributions cannot be taken as fully representative of a whole xez cross section at any level as is shown in Fig. 9(b). Wave like structures approximately parallel to the adiabatic seem to form. The difference between the maximum and minimum velocities grows with height, however, their positions remain approximately constant in the z direction. This is particularly obvious in the corners of the channel

18

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

Fig. 8. Comparison of time-averaged coefficient of heat transfer on the z/L ¼ 0.5 line on heated wall S1; Numerical: V1, VREF; Experimental:  V1, VREF.

as is shown in Fig. 9(b) for (y/H)  0.36, where particularly high velocities are generated. Corner vortices are causing local “hot spots” in the four corners at each cross-section thereby heating the air in the vicinity to a higher temperature than in its surroundings, as is shown by the temperature distributions in Fig. 9(b). This leads to locally higher body forces with the consequence that larger velocities result as the flow develops. Since the average absolute temperature of the fluid leaving the channel is only about 1.7% higher than its value at the inlet, the average velocities at the two cross-sections would not vary significantly. Indeed, the average velocity on the line 0  x  1 at z/ L ¼ 0.5, y/H ¼ 0.18, 0.36, 0.54, 0.8 and 0.95 of the numerically obtained data in Fig. 9(a) are all approximately 0.36 ms1. The velocities averaged over the whole cross-section at each height using the data presented in Fig. 9(b) are 0.364, 0.364, 0.365, 0.366, 0.367 ms1, respectively. These indicate that the velocity averaged on the centre plane is representative of the whole velocity profile. However, as may be seen in Fig. 9(a), the experimental obtained velocity profiles were lower than the distributions generated numerically numerical data. The estimated error in the values of velocities obtained experimentally was evaluated at ±15% [18], so that the error in the average velocity should also be ±15%. Unlike the numerically obtained data, the values of the experimentally determined vertical

velocity averaged over the line 0  x  1 at z/L ¼ 0.5 varied widely at the various heights mentioned above and were 0.32 ± 0.05, 0.31 ± 05, 0.26 ± 0.04, 0.30 ± 0.05 and 0.34 ± 0.05 ms1, respectively. The difference between the highest and the lowest average velocity found experimentally is 0.08 ms1, which is just inside the estimated error of 0.09 ms1. Thus all experimentally determined average values are plausible. The difference between the average value of all the experimentally determined velocities, 0.31 ms1 and the averaged numerically obtained velocity, 0.365 ms1, is 16% which, once again, is acceptable. It is also observed in Fig. 9(b) that there is a significant reduction of the velocity, particularly in the ‘core region’ (0.2  x/D  0.7) as the fluid rises in the channel. This is linked to a significant difference in the temperature which directly leads to large density differences between the heated and unheated portions of the fluid at any cross-section. The fluid near one of the heated wall becomes increasingly hotter as it rises, so that its velocity grows as the fluid rises. However, as there is little mixing, the fluid in the core regions in Fig. 9(b) at any level is hardly hotter than the fluid at inlet. As a result, the buoyancy forces in the central regions are low and the fluid in these areas has a low velocity. In fact this part of the flow is “pulled” up by shear stresses between the faster moving fluid near the heated walls and the slower travelling fluid in the core region. This local reduction of velocity in the core region has already been observed in natural convection in a uniformly heated right vertical circular channel [31,32]. Indeed, in extreme cases a stationary persistent separated region may occur in the central region of the tube [32]. The authors relate this phenomenon to a ‘competition’ between the shear forces generated by the relatively large velocities differences between the fluid near the heated walls and the slower moving fluid in the core and the force needed to support the ‘weight’ of fluid in the low temperature, therefore high density central region of the channel. In this non-uniformly heated configuration in which heated and unheated zones alternate and are not symmetrically located on each heated wall, the reduction of the velocity in the core region, as shown in Fig. 9 is particularly pronounced in the experimental results between y/H ¼ 0.36 and y/ H ¼ 0.80, whereas in the numerical results this phenomenon continues to the outlet. Beyond y/H ¼ 0.54, both experimental and numerical results indicate better mixing, with the streamwise velocity in the core region increased. 4.3. Instantaneous flow and temperature fields in the nonuniformly heated channel A separation of the thermal plumes formed as the result of heat transfer to the fluid from the heated foils has been previously observed both experimentally and numerically for uniformly heated channel for the VREF configuration [7]. However, as may be seen in the three numerically obtained snapshots at three times shown in Fig. 10 different behaviour can be observed for the nonuniform heating V1 configuration. For both the V1 and VREF configurations the plumes rise on both walls case V1, but only along the one heated wall in the VREF design. In the V1 case, these plumes move across the channel, impinge on the opposite wall and interact with the plumes rising on the opposite wall leading to significant mixing between the hot fluid at the both walls and the cool fluid in the central regions of the channel. However in the VREF case, despite having doubled the power supply, such large scale motions only occur towards the exit of the channel. Further, in the V1 arrangement, the plumes turn into travelling waves that continuously keep changing shape within the thickening boundary layer and appear to be ‘flap’ from one side of the channel to another. This seems to enhance mixing and produces a

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

Fig. 9. Streamwise time-averaged velocity Vy for case V1 at different heights of the channel:

Numerical;

19

Experimental.

20

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

Fig. 10. Instantaneous contours of the temperature and vertical velocity at mid-plane z/L ¼ 0.5.

strong coupling between instantaneous temperature and velocity fields can be clearly also be seen in Fig. 14, where the contour plots of temperature are also shown. There appears to be much less coupling between the velocity and temperature fields in the VREF case. This indicates that much less mixing occurs in the VREF case as is clearly seen in the three snapshots in Fig. 10. It can be observed in the V1 arrangement that areas near unheated zones on both S1 and S2 are also regions of lower velocity. As a consequence, jet like structures are generated. These grow and travel downstream and towards the opposite wall. The interaction between the two boundary layers brings cool fluid to the walls which results in higher coefficient of heat transfer as shown in Fig. 8. The degree of mixing, Md, defined as 0

Md ¼

0

smax  sj ; 0

smax

(19)

is one way of quantifying mixing. In Eq. (19) the non-dimensional standard deviation, s0j is defined as

0

sj ¼

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u  2 u u∬ 1  DT dA t A DT A

(20)

in which

DT ¼

∬ A DTdA ; DT ¼ T  Tref A

(21)

and A is the cross-section area of the duct. It should be noted that here the temperature T is the time average temperature. The variation of Md as a function of height in both the V1 and VREF arrangements is shown in Fig. 11. Except for the first y/ H ¼ 0.05, the V1 configuration yielded better mixing that the VREF setup over the rest of the height of the channel. This means that the air in the vicinity of the heated surfaces was cooler in the V1 than in the VREF configuration, leading to higher heat transfer coefficients in the V1 arrangement, as is clearly shown in Fig. 8. Surprisingly, at

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

Fig. 11. Degree of mixing vs height within the channel for cases V1 and VREF.

the very top, the value of Md in both configurations was approximately the same. As mentioned by Lau et al. [7] and seen in Fig. 10(b), relatively “hot” fluid is drawn in from the outside of the channel at the outlet from the VREF arrangement. The intermittent nature of this inflow, as may be seen in Fig. 10(b), leads to the appearance of large structures of hot fluid being ejected from the heated wall, thereby producing a more uniform temperature at exit when averaged over time. This is the reason for the large gradient in Md towards the ends of the channel. This inflow is not present in the case of the V1 configuration, so that, after an initial very sharp increase in Md there is little further increase in its value further downstream. The need for introducing external disturbance is clearly illustrated in Fig. 12(a and b) in which there is a large difference

21

between a typical numerically generated snapshot of the flow and a similar snapshot obtained experimentally. From comparison of instantaneous isotherms and streamlines obtained numerically with visualization of the velocity field obtained using PIV it is clear that the case with imposed disturbance represents more adequately the evolution of the dynamic behaviour of natural convective flows within the channel such as seen in Fig. 12. In the V1 arrangement, the unsteady structures within the fluid flow are generated at recirculation zones at the inlet and lower portion of the channel, particularly on the unheated zone at the inlet of the S2 wall as it has seven heaters only. This is similar to the VREF case in which a large separation occurred on the both heated and unheated walls [7]. The process of intermittent growth and decay of the separated structures in the V1 case, shown in Fig. 13(a), results in what could be called ‘flapping’. Four such large structures, which are more or less fixed in the z direction exist and lead to the non-uniform appearance of the streamlines on the S2 wall in Fig. 13(a). Similar, but, smaller and less well-defined structures also appear on the face S1 in Fig. 13(a). This also leads to less visible nonuniformities in the streamline distribution on the S1 wall as well. The non-uniformities generated at the inlet in the z direction persist all the way to the outlet as is clearly seen in Fig. 13(b and c) and result in the non-uniform temperature distributions on the heated walls shown in Fig. 13(c). The fact that the structures are reasonably fixed in the z direction has already been discussed above in connection with the non-uniformities shown in Fig. 9(b). The numerically calculated fluctuating separated zone at the leading edge of the S1 wall is not closed, so that fluid enters and leaves this region. The reduced fluid velocity together with constant power supply to the heated wall leads to the fluid near the leading edge becoming quite hot so, as may be seen in Figs. 13(c) and 7(a), the temperature of the boundary near the leading edge of S1 also increases. However, at y/H z 0.015 the flow reattaches. This brings cold fluid to the surface of S1 and leads to a sharp fall in the wall temperature as may be seen in Fig. 7(a). At y/H z 0.04 as shown in Fig. 7(a), the fluid near the boundary has been heated to sufficiently high temperature to allow the wall temperature to increase again, albeit at a slightly lower rate than the experimentally determined rate. The separation at the leading edge of the S1 face is

Fig. 12. Instantaneous contours of vertical velocity near the outlet at mid-plane z/L ¼ 0.5. (a) 0 disturbances; (b) Experimental (c) 0.02 disturbances.

22

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

Fig. 13. Effect of unheated zone: instantaneous, at time 100.56 s. (a) Wall temperature T; (b) Vertical velocity Vy; (c) Stream traces coloured by temperature.

also responsible for a similar numerically generated temperature distribution on the S1 surface in the VREF configuration. The above distribution in temperatures of the wall leads to a large change in the distribution of temperatures in the fluid. This is the reason for the ‘dip’ in the degree of mixing for both the V1 and VREF arrangements seen in Fig. 11 between 0.02 < y/H < 0.05. Typical patterns of evolution of fluctuations in time of velocity and temperature are shown in Fig. 14 at approximately at midheight of the channel (y/H ¼ 0.54, z/H ¼ 0.5) near the wall S1 at x/D ¼ 0.15 and the wall S2 at x/D ¼ 0.85. The fluctuation in both temperature and velocity on both positions are in phase with each other. However, the fluctuations at the two points are clearly not in phase with each other. Fig. 15 depicts the power spectral density (PSD), evaluated at the above-mentioned points for both numerical and experimental velocity values. Since there were only 340 s seconds of data obtained numerically after 3 months running on the fastest cluster available to the authors, the lowest frequency which could be obtained with any reliability was of the order of 101 Hz. On the other hand the experimental results obtained with 7000 double images taken [18]

allowed frequencies as low as approximately 102 Hz to be evaluated. On the basis of the FFT results shown in this Fig. 15, the values of frequency obtained from both experimental and numerical data near wall S1 in the range between 0.1 Hz and 1 Hz can be usefully compared. Meanwhile, near wall S2 the dominant frequency along all height of the channel remains 0.1 Hz, with another significant frequency between 0.2 Hz and 0.3 Hz. As can be seen the two PSD curves in Fig. 15, obtained numerically and represented by continuous lines are similar to the experimental measurements for Hz > 0.1; in particular, for the wall S2. The fall off rate on at x/ D ¼ 0.15 is lower for the numerical results than that obtained experimentally. 4.4. Induced mass flow rate estimation The mean mass flow rates per unit width of the channel (kg/ (m,s)) are presented in Fig. 16(a). Since no temperatures were measured experimentally within the channel, the experimental mass flow rate per unit width of the channel, m_ exp was evaluated with

Fig. 14. Instantaneous temperature T and vertical velocity Vy. (a) x/D ¼ 0.15: near wall S1; (b) x/D ¼ 0.85: near wall S2.

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

Fig. 15. Power spectral density in proximity to the walls (a) S1; (b) S2. dxd Experimental;

  1 m_ exp ¼ $A$Vy $r Tfl L

(22)

in which the density of the air, rðTfl Þ was based on the mean temperature Tfl , between the S1 and S2 wall temperatures [18]. This led to an estimated ±21% uncertainty [18] which is significantly greater than the ±15% uncertainty of the average velocity at a particular cross-section, discussed above. Further, also as discussed above, there is a significant variation in the values of the experimentally measured velocity at various cross sections, so that the average of the average velocity at each cross section, namely 0.31 ms1, was used calculating of the mass flow rate in the V1 case. The same procedure was used in the VREF arrangement. Unlike the experimental results, the mass flow rate at a height y, m_ nuy , calculated from the numerical results, used the temperature

23

0.02 disturbances V1.

of the air at the relevant points in evaluating the density of air so that

m_ nuy ¼

1 ∬ rðTðx; y; zÞÞVy dxdz L A

(23)

As continuity has to be satisfied, it was found that for both configurations m_ nuy had the same values at all positions in the channel. Since there are thermal plumes on either side of the channel in the V1 configuration rather just one, as shown in Fig. 16(a) both experimentally and numerically, for the same total power input, the mass flow rate for the V1 Case is higher than that in VREF Case.

Fig. 16. Global characteristics: comparison of different configurations. (a) Mean mass flow rates per unit width of the channel (kg/(m,s)), cases V1-VREF; (b) Mean NuD and Ra*ðD=aÞ numbers at wall S1; cases V1, V2, VREF for the range of the heat loads 18, 47, 100, 225, 475 W/m2) [18].

24

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25

4.5. Overall heat transfer characteristics In order to generalise the results, the wall heat transfer data have been recast in terms of dimensionless Nusselt, NuD , defined as

NuD ¼

D H

ZH 0

Fconv ðyÞ dy kðyÞ$ðTwall ðyÞ  Tinlet Þ

(24)

and the modified Rayleigh number, Ra*ðD=aÞ; which on the assumption that the air is a perfect gas so that the coefficient of thermal expansion is equal to the inverse of the absolute temperature, defined as

Ra*ðD=aÞ ¼

gD5 aHTinlet

ZH 0

Fconv ðyÞ dy; nðyÞkðyÞaðyÞ

(25) References

with both integrals being taken along the z/L ¼ 0.5 line. In the above expressions the temperature at the inlet of the channel was taken as a reference temperature whilst the thermo-physical properties were evaluated at the wall temperature Twall(y). The integrations in Eqs. (24) and (25) were performed over heated areas of the wall, only. This means that in the VREF Case values on the S1 wall only, needed to be used, whereas in the V1 Case values on the S1 and S2 walls were used alternatively as the integration was performed. NuD as a function of Ra*ðD=aÞ is plotted in Fig. 16(b) which also includes a series of other experimental and numerical data found in the literature. As may have been expected the present results are in good agreement with the previously obtained experimental and numerical data. It has been shown above that the temperatures of the walls in the V1 Case are lower than that of the S1 wall at the same height in the VREF case, so that the channel average Nusselt number is larger for the periodically heated configuration. In particular, the numerically obtained NuD at the same heat input is 34% larger in Case V1 than in Case VREF. It is of practical interest to have a global correlation regardless the size of alternated heated zones and non-heated zones. Sanvicente [18] proposed that

NuD ¼ 2:11$ðRa  ðD=aÞÞ0:14

The power spectral density generated from the numerically axial velocity fluctuations was in reasonable agreement with the results obtained experimentally in 0.1e1 Hz range. Despite running the program for more than three months on the fastest available computer cluster, frequencies below 0.1 Hz could not be identified. The channel heated from both side, albeit non-uniformly, has two plumes driven by natural convection, whereas the uniformly heated channel has only one. It was shown that the mixing, as measured by the parameter termed degree of mixing was more vigorous in the channel heated form both sides. This was confirmed by the higher wall temperatures in the case of the channel heated from one side. The channel heated from both sides also generated a larger mass flow rate.

(26)

is an appropriate correlation for this task. As may be seen in Fig. 16(b), the numerically generated data for both Case V1 and Case VREF almost exactly fall on the correlation line. 5. Conclusions The effects of the non-uniformity of the wall heat flux distribution on the natural convection flow and heat transfer within an open-ended vertical channel has been studied numerically and experimentally. Investigations were conducted for a configuration consisting of asymmetric alternately, heated and unheated zones on opposite sides of the channel: each zone having a height to channel height ratio of 1/15. The channel dimensionless width to channel height ratio of 1/15 and the heat flux was 220 W=m2 . Numerical results have been obtained using an in-house finitevolume solver and the predictions have been compared with the corresponding experimental results obtained in the CETHIL research centre in Lyon, France. As in previous numerical studies, it has been found necessary to introduce fluctuations at the inlet of the computational domain in order to reasonably match experimental data. It has been found here that the intermittently fluctuating separations at the inlet lead to the growth of flow structures which promote mixing within the channel.

[1] S. Kurtz, K. Whitfield, D. Miller, J. Joyce, J. Wohlgemuth, M. Kempe, et al., Evaluation of high-temperature exposure of rack-mounted photovoltaic modules, in: 34th IEEE Photovoltaic Specialists Conference, 2009. [2] S. Kurtz, K. Whitfield, G. TamizhMani, M. Koehl, D. Miller, J. Joyce, et al., Evaluation of high-temperature exposure of photovoltaic modules, Prog. Photovolt Res. Appl. 19 (2011) 954e965. [3] H.A. Zondag, Flat-plate PV-Thermal collectors and systems: a review, Renew. Sustain. Energy Rev. 12 (2008) 891e959. [4] M. Sandberg, B. Moshfegh, Buoyancy-induced air flow in photovoltaic facades: effect of geometry of the air gap and location of solar cell modules, Build. Environ. 37 (2002) 211e218. [5] B.J. Brinkworth, B.M. Cross, R.H. Marshall, H. Yang, Thermal regulation of photovoltaic cladding, Sol. Energy 61 (1997) 169e178. [6] M. Posnansky, A. Eckmanns, Practical results with cogeneration of electricity and heat in building integrated PV power systems, in: 13th EPSEC, Nice, 1995.
ne
zo, S. Giroux-Julien, M. Fossa, E. Sanvicente, et [7] G.E. Lau, V. Timchenko, C. Me al., Numerical and experimental investigation of unsteady natural convection in a vertical open-ended channel, Comput. Therm. Sci. 4 (2012) 443e456. [8] G.E. Lau, V. Timchenko, J.A. Reizes, M. Fossa, G.H. Yeoh, Natural convection in an asymmetrically-heated open-ended channel: a three-dimensional computational study, in: ASME 2013 Summer Heat Transfer Conference, HT2013, Minneapolis, MN, USA, 2013. [9] G.E. Lau, G.H. Yeoh, V. Timchenko, J.A. Reizes, Large-eddy simulation of turbulent natural convection in vertical parallel-plate channels, Numer. Heat Transf. Part B Fundam. 59 (2011) 259e287.
ne
zo, H. Bouia, Transitional natural [10] E. Sanvicente, S. Giroux-Julien, C. Me convection flow and heat transfer in an open channel, Int. J. Therm. Sci. 63 (2013) 87e104.
ne
zo, E. Leonardi, Experimental natural convection on vertical [11] M. Fossa, C. Me surfaces for building integrated photovoltaic (BIPV) applications, Exp. Therm. Fluid Sci. 32 (2008) 980e990. [12] L.H. Chen, H.Z. Tian, Y.Z. Li, D. Zhang, Experimental study on natural convective heat transfer from a vertical plate with discrete heat sources mounted on the back, Energy Convers. Manag. 47 (2006) 3447e3455.
ne
zo, J. Vareilles, H. Pabiou, M. Fossa, E. Leonardi, [13] S. Giroux-Julien, C. Me Natural convection in nonuniformly heated channel with application to photovoltaic facades, Comput. Therm. Sci. 1 (2009) 231e258. [14] J. Hern
andez, B. Zamora, Effects of variable properties and non-uniform heating on natural convection flows in vertical channels, Int. J. Heat Mass Transf. 48 (2005) 793e807. [15] A.K. da Silva, G. Lorenzini, A. Bejan, Distribution of heat sources in vertical open channels with natural convection, Int. J. Heat Mass Transf. 48 (2005) 1462e1469. [16] A.A. Dehghan, M. Behnia, Numerical investigation of natural convection in a vertical slot with two heat source elements, Int. J. Heat Fluid Flow 17 (1996) 474e482. [17] A.K. da Silva, S. Lorente, A. Bejan, Optimal distribution of discrete heat sources on a wall with natural convection, Int. J. Heat Mass Transf. 47 (2004) 203e214. [18] E. Sanvicente, Experimental Investigation of Thermal and Fluid Dynamical Behavior of Flows in Open-ended Channels: Application to Building Integrated Photovoltaic (BiPV) Systems, National Institute of Applied Sciences e INSA Lyon, 2013 (PhD thesis). [19] G.E. Lau, G.H. Yeoh, V. Timchenko, J.A. Reizes, Large-eddy simulation of natural convection in an asymmetrically-heated vertical parallel-plate channel: assessment of subgrid-scale models, Comput. Fluids 59 (2012) 101e116. [20] G.E. Lau, G.H. Yeoh, V. Timchenko, R.K.K. Yueen, Natural convection in a PVintegrated double-skin façade using large-eddy simulation, Procedia Eng. 14 (2011) 3277e3284. [21] A.W. Vreman, An eddy-viscosity subgrid-scale model for turbulent shear flow: algebraic theory and applications, Phys. of Fluids (1994-present) 16 (2004) 3670e3681. [22] W. Sutherland, The viscosity of gases and molecular force, Philos. Mag. (S.5) 36 (1893) 507e531.

O.A. Tkachenko et al. / International Journal of Thermal Sciences 99 (2016) 9e25 [23] D.G. Barhaghi, L. Davidson, Natural convection boundary layer in a 5:1 cavity, Phys. Fluids (1994-present) 19 (2007).
[24] J. Vareilles, Etude des transferts de chaleur dans un canal vertical
rentiellement chauffe
: application aux enveloppes photovoltaïquesdiffe
Claude Bernard Lyon 1, 2007 (PhD thesis). thermiques, Universite [25] G.E. Lau, G.H. Yeoh, V. Timchenko, J.A. Reizes, Numerical investigation of passive cooling in open vertical channels, Appl. Therm. Eng. 39 (2012) 121e131. [26] M. Miyamoto, Y. Katoh, J. Kurima, H. Sasaki, Turbulent free convection heat transfer from vertical parallel plates, in: 8th International Heat Transfer Conference, 1986, pp. 1593e1598. [27] T.C. Papanastasiou, N. Malamataris, K. Ellwood, A new outflow boundary condition, Int. J. Numer. Methods Fluids 14 (1992) 587e608.

25

[28] G.E. Lau, Natural Convection in Building-integrated Photovoltaic Systems: a Computational Study, UNSW, 2013 (PhD thesis).
ne
zo, et al., [29] G.E. Lau, E. Sanvicente, G.H. Yeoh, V. Timchenko, M. Fossa, C. Me Modelling of natural convection in vertical and tilted photovoltaic applications, Energy Build. 55 (2012) 810e822. [30] J.R. Ehlert, T.F. Smith, View factors for perpendicular and parallel, rectangular plates, J. Thermophys. Heat Transf. 7 (1993) 173e174. [31] T.F. Ayinde, S.A.M. Said, M.A. Habib, Experimental investigation of turbulent natural convection flow in a channel, Heat Mass Transf. 42 (2006) 169e177. [32] A. Litvak, G.L. Morrison, Study of flow characteristics in thermosyphon solar system, in: Eighth Australasian Fluid Mechanics Conference, University of Newcastle, New South Wales, Australia, 1983, 11C.5e11C.8.