Numerical simulation of conjugate forced turbulent

International Journal of Thermal Sciences 87 (2015) 121e135

Contents lists available at ScienceDirect

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

Numerical simulation of conjugate forced turbulent heat convection with induced natural laminar convection in a 2D inner cavity Diego A. Vasco a, *, Carlos Zambra b, c, Nelson O. Moraga d a

Department of Mechanical Engineering, University of Santiago de Chile, Av. Lib. Bdo. O'Higgins 3363, Santiago 9170022, Chile University Arturo Prat, Av. Arturo Prat 2120, Iquique 1110939, Chile c Centro de Estudios en Alimentos Procesados (CEAP), CONICYT-Regional, GORE Maule, R09I2001 Talca, Chile d Department of Mechanical Engineering, University of La Serena, Av. Benavente 980, La Serena 1720170, Chile b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 10 May 2014 Received in revised form 31 July 2014 Accepted 9 August 2014 Available online

Unsteady natural convection process of a fluid inside an inner walled vessel caused by an external forced turbulent convective flow of air in a rectangular duct has been analyzed numerically for Re ¼ 5000; Re ¼ 22,000. Fluid mechanics and conjugate convective heat transfer, for both inner fluid (water) and outer fluid (air), described in terms of continuity, linear momentum, keε turbulence model and energy equations were predicted by using the finite-volume method (FVM). Results for the streamlines, isotherms, local Nusselt numbers along the cavity walls and the average friction factor are presented. A comparison between obtained results for both Reynolds number is made. In general, it was found that the cooling process is improved up to 79% when Re number is increased from 5000 to 22,000; and the fluid mechanics characteristics inside the vessel is not considerable affected by the velocity inlet of the cooling air. © 2014 Elsevier Masson SAS. All rights reserved.

Keywords: Conjugate heat transfer Turbulence Natural convection Finite-volume method

1. Introduction Conjugate heat transfer processes are found in several relevant engineering applications. Heat and mass transfer processes such as drying, cooling, freezing, quenching and many others thermal treatments are usual in metallurgy, electronics, chemical and food industry. Specially, several conjugate processes in food industry are the sum of physical complexities such as the turbulent flow of external fluids and the heating/cooling/freezing of a liquid food with concurrently natural convection. These complexities must be captured in the mathematical model through the implementation of conjugation conditions and additional transport equations (i.e. keε model) and by its representation with suitable approximations for the variations of the thermo-physical properties. The numerical simulation of conjugate processes implies the solution of several systems of differential equations because of the variation of the thermo-physical properties or through splitting of the domain in subdomains. Therefore the mathematical model should be defined along with conjugation conditions [1]. Through the extension of the computational physical situation a better

* Corresponding author. Tel.: þ56 (51) 27183120. E-mail address: [email protected] (D.A. Vasco). http://dx.doi.org/10.1016/j.ijthermalsci.2014.08.008 1290-0729/© 2014 Elsevier Masson SAS. All rights reserved.

understanding of the interaction between the domains of interest as well as of the surrounding fluid are possible. Moreover, the calculation of properties of interest such as the convective heat coefficient and the friction factor around the domain of interest can be calculated a posteriori [2]. Kim and Viskanta [3] performed one of the first numerical and experimental studies of conjugate processes. These authors analyzed the effect of three differentially heated walls configurations and the ratio of the thermal conductivities on the natural convection process. The experiments show the importance of the effect of the walls on the convective heat transfer and on the Nusselt number. Liaqat and Baytas [4] performed a transient numerical study of the heat transfer process in a cavity with heat conductive walls using FVM and the SIMPLER algorithm of the heat transfer process in a cavity with heat conductive walls. The authors analyzed the effect of the Rayleigh number (107e1012) and the thermal conductivities ratio on the fluid mechanics and heat transfer process. It was found that the obtained results are qualitatively different from those obtained in nonconjugate domains. House et al. [5] analyzed the effect of a heat conductor body inside a cavity in the natural convection process through FVM and the SIMPLER algorithm. The authors found that from a critical size of the conductive body, the effects of the body size and the thermal conductivity ratio become important in the calculated Nusselt number. For smaller sizes than the critical size of

122

D.A. Vasco et al. / International Journal of Thermal Sciences 87 (2015) 121e135

Nomenclature a Cf Cp C1,2,V d dþ D E g G k H L Nu p Re T u,v U, V x, y X, Y

under relaxation factor friction coefficient heat capacity [J/kg K] parameters of dissipation of energy inner characteristic length of the cavity (m) dimensionless characteristic length of the laminar layer outer characteristic length of the cavity (m) dimensionless wall function parameter (E ¼ 9.793) gravitational acceleration [m/s2] generation of turbulent kinetic energy kinetic turbulent energy width of the channel (m) length of the channel (m) Nusselt number [Nu ¼ hD/l] absolute pressure [Pa] Reynolds number [Re ¼ rLu∞/m] temperature [K] velocity components (m/s) dimensionless velocities [U ¼ ui,j/u∞; V ¼ vi,j/u∞] cartesian coordinate position cartesian dimensionless coordinate position [X ¼ x/D; Y ¼ y/D]

Greek symbol b blockage ratio (b ¼ D/H) l thermal conductivity [W/m K] ε dimensionless parameter of the dissipation equation

the body, the heat transfer inside the cavity, in comparison with the classic natural convection process, can be improved (worsen) by the presence of a conductive body with a thermal conductivity ratio lower (greater) than unity. Yeong-Ha et al. [6] performed a transient MVF -based natural convection in a square cavity with a heat generating body inside. The movement of the fluid is leaded by two mean temperature gradients; a temperature gradient along the cavity because of the isothermal walls and another one caused by the heat source. These authors studied the effect of both temperature gradients and the dimensionless numbers Rayleigh and Prandtl numbers (melted sodium, air and water) on the heat transfer and fluid mechanics. Later, Yeong-Ha and Jung [7] extended the analysis to a three-dimensional domain. In this case, the effect of the Rayleigh number, for a constant Pr number (melted sodium), and the heat generation on the fluid mechanics and heat transfer was studied. The main input of this work was the analysis of the impact of the third dimension in the conjugate heat transfer process and on the two-dimensional distribution of the Nusselt number in the isothermal walls. Another contribution to the original investigation of House et al. [5] was made by KumarDas and Kumar-Reddy [8], in which the authors analyzed the effect of the Rayleigh number (103e106), the thermal conductivity ratio (l*) and the inclination angle of the cavity (15 e90 ). It was found that for Ra>103 there is a critical point where the average Nusselt number changes its relative magnitude. Under the critical point, a higher l* enhances the heat transfer and above this point a lower l* causes further improvement in the heat transfer than the value obtained with a higher l*. Moreover, in this study the sudden change of the thermo-physical properties at the interface was considered through a harmonic mean formulation [9]. This is remarked due to the importance of the required treatment by the

h q k m r t

apparent viscosity [Pa s] dimensionless temperature [q ¼ (TTC)/(THTC)] constant of the universal profile of velocity dynamic viscosity [Pa s] density [kg/m3] dimensionless time [t ¼ u∞t/D]

Subscripts 0 reference value 1, 2 differential subscripts of turbulence standard keε turbulence model a air c characteristic turbulent velocity or distance of the standard keε model C cold f friction H hot i, j cartesian component l laminar m inner fluid (water) s solid tkε turbulent property T top U, D distance of body edges from the inlet (Lu) and outlet (LD) planes w wall ∞ average Superscripts k iteration number

thermo-physical properties of the fluid and the solid when coupling SIMPLE-based algorithms are implemented [10]. Considering the several applications of the conjugate processes, other works have been focused on the effect of the position of heat sources [11] and mass source bodies [12,13] in the heat transfer inside the cavity. The turbulent flow around objects is of great interest in aeronautical, environmental and food industry. The pioneering experimental study performed by Lee [14] is one of the more referenced works in the area. In this study, the author measured the fluctuation of the pressure for a turbulent flow around a squared object with Re ¼ 1.8 105 and analyzed the effect of the position of the object respect to the main direction of the flow on the vortex formation and on the drag coefficient. Lyn and Rodi [15] and Lyn et al. [16] provided a set of experimental measurements of the velocity and turbulent quantities for the turbulent flow (Re ¼ 22,000) around a square object. This experimental work, along with many other investigations, has provided a strong basis for computational investigations of turbulent flows around objects focused as well on the turbulence modeling. Kamawura and Kawashima [17] validated their results with the keε standard model for low Reynolds numbers with the experimental results obtained by Lyn and coworkers [16]. The authors found a good agreement between the calculated Strouhal number and drag coefficient with those values obtained experimentally. In a set of previous numerical works by FVM, conjugated laminar mixed heat transfer with natural convection/solidification of a fluid in an inner vessel contained in a cavity has been studied. In the first work [18], water cooling considering opposite outlet position of the cooling fluid (Re ¼ 200, 500) respect inlet was studied. Moraga et al. [19] studied the effect of the Reynolds


number and the rheological behavior of the solidifying fluid in a conjugate heat transfer process. Lately, the effect of the inlet/outlet position in the cavity of the cooling fluid on the solidification process of a non-Newtonian fluid was analyzed [20]. LemusMondaca et al. [21] studied the unsteady natural convection driven sterilization of non-Newtonian Ostwald-de Waele carboxymethyl-cellulose solution inside a square container caused by external natural convection of a non-Newtonian. In all these studies, food industry processes were the main subjects of interest. The main purpose of the present work is to continue with the aforementioned investigations by predicting the fluid mechanics and heat transfer in a conjugate two-dimensional process of a forced turbulent flow of air (Re ¼ 5000 and Re ¼ 22,000) with an induced natural convection of a Newtonian fluid (water) inside an inner container. A study with two Re numbers (Re ¼ 5000 and 22,000) was performed to emphasize the effect of the outer flow mechanics in the fluid dynamics and natural convective heat transfer of the inner fluid. The transient friction factor and Nusselt number at the outer walls of the inner vessel are calculated during the different time scales of the cooling process. The work starts with the description of the physical situation and the mathematical model, including the constitutive equations that describe the thermo-physical properties of the inner fluid. In the third section the computational implementation of the FVM code is described. In the fourth section, the validation and verification of the numerical code by comparison of the obtained results for turbulence and natural convection with those experimentally reported is presented, along with the grid size and time step analyses. Finally, the obtained results for the conjugate forced turbulent convection with natural convection of water are shown and discussed.

2. Physical and mathematical model In food industry applications such as cooling and freezing of liquid food in coolers the external flow of air is turbulent in order to reduce cooling times and space requirements, while the inner flow in the vessel is laminar due to the vessel dimensions. Fig. 1 represents the two-dimensional conjugate domain of the forced turbulent convection with natural convection of a Newtonian fluid. This physical domain has been studied experimentally [15,16] and numerically [22]. The dimensions of the channel in Fig. 1a are defined according to the characteristic length of the cavity D (0.044 m) and a blockage ratio (b ¼ D/H) of 0.07. To avoid the three-

123

dimensional effects, the channel side along the lateral direction is assumed to be very long. The distances of the obstacle edges from the inlet and outlet planes of the channel are set to LU ¼ 4D and LD ¼ 14D, respectively. The cooling fluid (air) gets into the channel with a constant velocity U∞ and a constant temperature (T∞). The surfaces of the channel are considered adiabatic. The Reynolds number is calculated on the inlet velocity, the air properties and the characteristic length D. For the present study, two Reynolds number values are chosen. The value of Re ¼ 22,000 is chosen because the results for the asymptotic isothermal case can be compared with reliable experimental data [16], meanwhile a lower Reynolds number value of 5000 is used because it is considered a limiting value for which a flow can exhibit a turbulent flow motion. The immersed body (Fig. 1b) is a vessel with thin solid walls (d ¼ 0.038 m) and filled with a fluid, initially at TH > T∞. Inside the vessel density variation originated and maintained by temperature gradients is important that is why the inner fluid is subjected to the buoyancy forces and therefore laminar natural convection during the cooling process is considered. The governing equations are formulated for both incompressible fluids (water and air) and the solid walls. The variation of the thermo-physical properties with temperature is only considered important for the water contained in the inner vessel. The twodimensional mathematical model for each fluid and the vessel walls are specified below. 2.1. Water Density is considered an important function of temperature but not a function of pressure, that is the reason why water is considered an incompressible fluid, nevertheless the continuity equation is expressed by:

vr vðruÞ vðrvÞ þ þ ¼0 vt vx vy

(1a)

the laminar NaviereStokes equations are given by:

vðrm uÞ vðrm uÞ vðrm uÞ vP þu þv þ VðmVuÞ ¼ 0 vt vx vy vx

(2a)

vðrm vÞ vðrm vÞ vðrm vÞ vP þu þv þ Vðmm VvÞ rm g ¼ 0 vt vx vy vy

(2b)

and the energy equation is:

v ðrCpÞm T v ðrCpÞm T v ðrCpÞm T þu þv Vðlm VTÞ ¼ 0 vt vx vy

(3)

2.2. Air The thermophysical properties of air are considered constant with temperature, moreover air is considered incompressible in this application due to the Mach number is much lower than 0.3, therefore the continuity equation corresponds to Eq. (1b):

vu vv þ ¼0 vx vy

Fig. 1. Two-dimensional conjugate physical domain: forced turbulent convection of air in a channel (a) with natural convection of a fluid inside the inner container of side D, detailed in (b).

(1b)

The flow is considered turbulent and the keε model with the coefficients proposed by Lauder et al. [24] is selected instead of for instance keu or LES models because is the most used and tested of all turbulent models. The keε model has been successfully applied to a large number of different flow situations such as plane jets, mixing layers, boundary layer flows and internal duct flows [25]. In

124


the momentum equations, the velocity components and pressure are expressed as average values and the viscosity of air is augmented by the turbulent viscosity (mtkε):

ra

vu vu vu vP þu þv þ V½ðma þ mtkε ÞVu ¼ 0 vt vx vy vx

(4a)

ra

vv vv vv vP þu þv þ V½ðma þ mtkε ÞVv ¼ 0 vt vx vy vy

(4b)

The air turbulent equations for the kinetic turbulent energy (k) and its dissipation (ε) are expressed considering the keε standard model [23]:

vε vε vε ma þ mtkε þu þv V Vε vt vx vy sε 2 ε ε C2 ra ¼0 þ C1 G k k

ra

(5)

vk vk vk ma þ mtkε þu þv V ra Vk þ G ra ε ¼ 0 vt vx vy sk

(6)

Wall functions to the keε standard model should be added. The hypothesis establishes that in two points located out of the laminar layer the velocity component of the velocity parallel to the surface follows a logarithmic law and the turbulence is in local equilibrium, that is, the production of turbulent kinetic energy is equal to its dissipation. Thus, the magnitude of the velocity component parallel to the wall Uc and to a distance dc is given by:

Uc 1 þ ¼ ln Edc k Uf

(13)

with the von Karman constant of the universal profile of velocity k ¼ 0.418 and the friction velocity Uf and dþ c given by the Eqs. (14) and (15), respectively:

Uf ¼

qffiffiffiffiffiffiffiffiffi tw=r

(14)

dþ c ¼

Uf dc ml

(15)

where ml is the laminar viscosity. The terms for the production and dissipation of turbulent energy are expressed as follows:

where the turbulent viscosity is defined as:

k2 mtkε ¼ Cu ra ε

(7)

( "

) 2 # vu 2 vv vv vu 2 þ þ þ vx vy vx vy

(8)

The parameters of the dissipation Eq. (5) are those proposed by Launder & Spalding [24]:

C1 ¼ 1:44; C2 ¼ 1:92; Cu ¼ 0:09; sk ¼ 1:0; sε ¼ 1:0; PrT ¼ 0:85 (9) Finally, the energy equation in terms of the average properties of velocity and temperature is given by:

ðrCpÞa

vT vT vT þu þv vt vx vy

!

m V la þ Cpa tkε VT ¼ 0 PrT

(10)

2.3. Vessel walls The energy equation of the walls of the vessel (w), made of stainless steel, are giving by Eq. (11):

ðrCpÞs

vT ls V2 T ¼ 0 vt

(11)

2.4. Boundary and conjugation conditions The mathematical model is stated together with a set of additional restraints, the boundary and conjugation conditions. No-slip boundary conditions for the velocities are considered at the walls of the channel. For the energy equation at the walls of the channel, the boundary conditions consider adiabatic walls in the channel with null heat flux:

vT

¼0 vy y¼0;y¼H

(16)

Finally, the conjugation conditions for the inner (17) and outer walls (18) of the vessel are denoted by:

and the generation of turbulent kinetic energy (G) is given by:

G ¼ ðma þ mtkε Þ 2

Uf2 Uf3 kc ¼ pffiffiffiffiffiffi; εc ¼ kdc Cu

(12)

ls

vTs

vTm

¼ l m vx w vx w

(17a)

ls

vTs

vTm

¼ lm

vy w vy w

(17b)

Ts jw ¼ Tm jw

(17c)

vTs

vT a

ls ¼ la

vx w vx

w

vTs

vT a

ls ¼ la

vy w vy

w

(18a)

(18b)

Ts jw ¼ T a w

(18c)

2.5. Thermo-physical properties Water is the liquid assumed to be inside the inner vessel. The fourth order polynomial that describes the density of water as a function of the temperature ( C) in the interval between 273 and 293 K is given by Ref. [26]:

rðTÞ ¼ 6:621 107 T 4 þ 8:78 105 T 3 0:00894T 2 h . i þ 0:06732T þ 999:84 kg m3

(19)

Heat capacity, Eq. (20), the thermal conductivity, Eq. (21) and the Newtonian viscosity, in Eq. (22) are expressed in terms of the absolute temperature in the interval between 273 and 283 K according to [27]:


CpðTÞ ¼ 1:785 104 T 3 0:1915T 2 þ 67:95T 3755:9½J=kgK (20) lðTÞ ¼ 4:2365 109 T 3 1:144 105 T 2 þ 7:196 103 T

Cf ¼

2m g_ w ru2∞

3:8208 102 ½Pa s T 252:33

vu

vy y¼H

(21)

g_ w jy¼H ¼

(22)

with the velocity gradient being calculated with a second order accuracy by:

( ! " !

ujy¼Hþ1 ujy¼H Dyjy¼Hþ1 ujy¼Hþ2 ujy¼Hþ1 vu

¼ vy y¼H Dyjy¼Hþ1 Dyjy¼Hþ1 þ Dyjy¼Hþ2 Dyjy¼Hþ2

The thermo-physical properties of the cooling fluid (air, m ¼ 2) and the walls of the vessel (stainless steel, w) are considered constants, as given by expressions (23) and (24), respectively.

. rw ¼ 7750 kg m3 ; Cpw ¼ 0:420 J=kg K; lw ¼ 11:2 W=mK (23) rair

. ¼ 1:293 kg m3 ; Cpair ¼ 1:005 kJ=kgK;

lair ¼ 0:0243 W=mK; mair ¼ 1:33 105 Pa s

(24)

2.6. Nusselt number and friction factor calculation The heat transfer in the outer walls of the vessel is calculated through the local and average Nusselt number as a function of the non-dimensional temperature gradients between the vessel walls and the cooling fluid using the Eqs. (25) and (26):

hY vq

Nuw ðYÞ ¼ ¼

la vX w

(25)

Z Nuw ¼

Nuw ðYÞdY

T T∞ TH T∞

(30)

ujy¼Hþ1 ujy¼H Dyjy¼Hþ1

!#

þ O Dy2

) (31)

3. Solution procedure The governing coupled non-linear second order partial differential Eqs. (1)e(6), (10) and (11) with the no-slippery boundary conditions, thermal boundary conditions (12) and wall functions (13e16) were solved with the finite volume method, FVM, by using the SIMPLE algorithm [28]. A fifth power law was used to calculate the convective terms while the diffusion terms were determined by using linear interpolation functions for the dependent variables between the nodes. The resulting systems of algebraic equations were solved through the block Gauss-Seidel (backward/forward) iteration with regard to planes in Y-direction and an alternating direction line-by-line iteration in each plane (TDMA algorithm). Successive under-relaxation was used during the iterative process in the calculation of the improved values for the dependent primitive variables. The calculations were accomplished by using the supercomputing infrastructure of the NLHPC (Intel Xeon X5550 quad-core, 2.67 GHz; 24 GB memory). The INTEL Fortran compiler IFORT under Linux platform was used for the compilation of the Fortran90 in-house made source code. All variables were defined as double-precision floating-point data types. The under-relaxations factors used for the velocity components, turbulence variables, temperature and pressure were:

(26)

au;v ¼ 0:3; ak;ε ¼ 0:1; aT ¼ 0:5; ap ¼ 0:1

(27)

The convergence criteria applied to stop the velocity, turbulence and temperature calculations at each time step was based in a maximum value for the difference in the calculated value at two successive iterations, defined as follows:

The non-dimensional temperature is defined by:

q¼

(29)

in this expression ðg_ w Þ is the local deformation rate given by:

0:6326½W=mK

mðTÞ ¼

125

Numerically, the local Nusselt number is calculated from an implicit second order error expression:

qw1 qw DXw1 qwþ2 qwþ1 Nuw ðYÞ ¼ DXw DX þ DXw2 DXwþ2 w1 qwþ1 qw þ O DX 2 DXwþ1

(32)

k

5 k k1

5

ui;j uk1 i;j 10 ; vi;j vi;j 10

(33a)

k 4 k k1

4 10 ; ε

ε

ki;j kk1 i;j i;j i;j 10

(33b)

k k1

106

Ti;j Ti;j

(33c)

(28) where qwþ1, qwþ2 are non-dimensional values of temperature of the two nearest nodes perpendicular to the outer walls of the vessel at qw, which value is obtained through an energy balance at the interface. The local friction factor is calculated through the Eq. (29) in the outer walls of the vessel and the duct walls.

4. Validations In order to prove the correctness and accuracy of the obtained results by our MVF code, two different experimental physical situations were replicated, the laminar natural convection of water in


a squared cavity [29] and the turbulent forced convective flow around a square cylinder submerged in a duct [16]. Along with the validation of the numerical code, the grid size and time step analysis is performed.

4.1. Grid size study The first validation and grid size study of the code was performed by comparing the obtained results with those experimentally and numerically reported by Banaszek et al. [29] for laminar convection of water. The numerical results were obtained by using non-uniform adaptive grids of 42 42 and 62 62 nodes and a constant time step (Dt ¼ 0.5 s). The authors implemented the finite element method to study the natural convection in a cavity with L ¼ 0.038 m. In this numerical experiment, the right and left walls were kept at constant temperature TH (10 C) and TC (0 C), respectively; meanwhile the remaining walls were considered to be adiabatic. In Figs. 2 and 3 a comparison of the isotherms, streamlines and the vertical component of the velocity along the horizontal axis between the calculated and experimental values at time t ¼ 2000 s is shown. In general, the obtained results are in good agreement with the numerical results and experimental observations reported by Banaszek et al. [29]. It can be noticed the importance of the effect of the natural convection in the heat transfer process. The second process of interest is the turbulent forced convective flow around a square cylinder. For the verification of our code for turbulent flows, the numerical study through FVM and the SIMPLE coupling algorithm performed by Raisee et al. [22] was also replicated. In this case, the computational domain consists of a square cylinder with a blockage ratio (b ¼ D/H) of 0.07, which is located in the middle of the channel. The Reynolds number based

v (mm/s)

126

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8

Grid 2 Grid1 Experimental

0.0

0.2

0.4

0.6

0.8

1.0

X Fig. 3. Comparison of the vertical component of the velocity along the horizontal axis after t ¼ 2000 s (Grid 1: 42 42 nodes, grid 2: 62 62 nodes).

on the inlet velocity and the obstacle side (D) is 22,000. The Strouhal number is defined as (St ¼ fD/uin), where f is the shedding frequency of the vortexes. The values of the Strouhal number reported by Raisee and co-workers (0.115 and 0.126) obtained with the unsteady linear and non-linear keε turbulence model respectively are compared with the experimental value (0.132) obtained by Lyn et al. [16]. Implementing a non-uniform adaptive grid of 122 117 nodes and a time step of 104, a Strouhal number of 0.120 was obtained with the keε standard turbulence model, i.e. a relative error of 9.1% respect to the experimental value reported by Lyn and co-workers. In Fig. 4, the time-averaged U and V velocity profiles obtained by time-averaging process over one period are presented and compared with the counterpart experimental data of Lyn et al. [16]. The locations are chosen to be on the top of the immersed body (x/D ¼ 0.0) and within the recirculation region (x/D ¼ 1.5). At the first location there is a good agreement between the calculated and the experimental values for both velocity

Fig. 2. Comparison of the distribution of isotherms (up) and streamlines (down) obtained by Banaszek et al. [29] (left) and those obtained in the present work (right) after t ¼ 2000 s.

D.A. Vasco et al. / International Journal of Thermal Sciences 87 (2015) 121e135 2

2

Present study Experimental

1.75 1.5

1.5 1 0.75 0.5

0.5

a

0.25

U/Uin

0 -0.6 -0.35 -0.1 2

0.15

0.4

0.65

0.25

0.9

1.15

1.5

0.0

0.1

0.2

-0.1

0

0.3

2


1.75 1.5 1.25

y/D

1 0.75

0.5

0.5

c

0.25 0 -0.25

y/D

1.25 0.75

-0.1

1.4

b

V/Uin

0


1.75

1

y/D

1.25

y/D

0.75


1.75

1.25 1

127

U/Uin 0

0.25

0.5

0.75

1

1.25

d

0.25 V/Uin

0 1.5

-0.3

-0.2

0.1

Fig. 4. Comparison of predicted stream-wise and cross-stream velocity profiles at (x/D ¼ 0.0) (a, b) and at (x/D ¼ 1.5) (c, d) with experimental data of Lyn et al. [16] (b ¼ 7% and Re ¼ 22,000).

components. Although the correspondence between the results in the second location is not as good as the observed for the first locations, the shape of the predicted velocity profile denotes the correctness of the predicted results. According to Fig. 3, a grid with 42 42 nodes was fairly good to describe the process of natural convection of water that takes place in the inner vessel. The calculated Strouhal number for different grid sizes used in the calculations is shown in Table 1. According to these values a grid size of 122 117 nodes is good enough to capture the fluid dynamics of the turbulent process with Re ¼ 22,000. The problem of the conjugate forced turbulent convection of air with natural convection of water was built based on the contribution of the two previously validated works available in the literature. Therefore, to assure the correctness of the predicted heat transfer results a higher nodes concentration was kept in the vessel region, explained why a grid size of 153 138 nodes for the whole domain was chosen, with 50 50 nodes located in the inner vessel to capture the natural convection.

4.2. Time step study Based on the numerical study performed by Raisee et al. [22], it is known that low time steps are required to obtain correct results for the turbulent forced convection of air. That is why the effect of the time step (103 and 104 s), with a grid size of 122 117 nodes, on the Strouhal number was analyzed. The results obtained of St ¼ 0.1199 and 0.1198 were obtained with a time step of 104 s and 103 s, respectively. Moreover as shown in the previous section, a time step of 0.5 s is good enough to describe the natural convection of water inside the vessel. According to this, a time step of 103 s is appropriated to capture the fluid dynamic and heat transfer in the coupled case.

Table 1 Calculated Strouhal number (St) obtained for the forced turbulent convective flow around a body (Re ¼ 22,000). Grid 47 47 87 82 105 102 122 117 172 162 Experimental size Lyn et al. [16] St 0.008 0.125 0.116 0.120 0.119 0.132

5. Results and discussion The analysis and discussion of the results is divided in two sections. The first one gathers the results dealing with the fluid mechanics of the surrounding air and those corresponding to water inside the vessel. The transient behavior of streamlines, velocity profiles and the average friction factors calculated around the vessel and along the walls of the channel are reported. The second section comprises those results related to heat transfer during the cooling process of the fluid contained the inner vessel. These results include the transient behavior of isotherms, temperature profiles and cooling rates in the inner vessel and finally the Nusselt number calculated at the outer walls of the vessel. 5.1. Fluid mechanics The transient behavior of the streamlines of air along with selected isotherms of the water in the vessel during the convective process is shown in Fig. 5 for Re ¼ 5000; in the time interval between 1150 s and 1200 s, and in Fig. 6 during the time interval between 600 s and 650 s for Re ¼ 22,000. Those time intervals correspond to the last 50 s of the cooling process of water contained in the inner vessel. The cooling time process is diminished up to 79% when the Reynolds number is increased from 5000 to 22,000. In Fig. 5 can be observed the vortex formation on both top and bottom surfaces. Then, how each vortex moves (Fig. 5a,d) toward the backside of the vessel, grows (Fig. 5d,e) and sheds (Fig. 5c,f). When the Reynolds number becomes higher the behavior of the fluid mechanics around the vessel is quite different. The vortexes formed on the top and bottom surfaces of the vessel do not grow alternatively on the right surface as in the preceding case, but they grow at the same time (Fig. 6b,c) and the shedding occurs when one of the growing vortexes is big enough (Fig. 6c,e). It can be noticed that when a vortex sheds is replaced by upcoming vortex from the vessel surface (Fig. 6a,f). It is important to notice how the outer fluid motion induces the heat transfer process in fluid contained in the inner vessel. The U and V velocity profiles along the horizontal axis at the middle of the duct are shown in Fig. 7 for both Reynolds numbers. After the fluid cross the vessel region, the vertical component of velocity acquires non-zero values, pointing downwards when Re ¼ 5000 and upwards when Re ¼ 22,000. In both cases a wavy

128


Fig. 5. Transient behavior of the streamlines around the vessel during the final instants of the cooling process of water with a forced convective flow of air (Re ¼ 5000).

behavior of the vertical component of the velocity is noticed, being practically independent of time when Re ¼ 5000. The horizontal component of velocity presents several inflexion points after the vessel zone. This can be noticed clearly after a distance equals to approximately four and six times the vessel length for Re ¼ 5000 and Re ¼ 22,000, respectively. The average friction factor calculated on the surface of the duct is practically independent of time and it is almost 50% lower when the Reynolds number is increased from 5000 (Cf ¼ 5.73 105) to 22,000 (Cf ¼ 2.84 105). Fig. 8 shows the transient behavior of the average friction coefficient calculated at the vessel. The apparently chaotic behavior observed when Re ¼ 5000 (Fig. 8a) is explained by the more complex vortex forming and shedding process, showed in Fig. 5. When the inlet velocity is increased, Re ¼ 22,000, an oscillating behavior is observed since the vortex forming and shedding process occurs at a higher frequency (Fig. 8b). The chosen scale of time does not allow to distinguish the differences observed when Re ¼ 5000.

The water contained in the vessel is subjected to natural convection promoted by the temperature gradients during the cooling process. Fig. 9 shows the transient vertical component of velocity along the horizontal and vertical axis for both Reynolds numbers. The V-velocity component profiles along the horizontal axis share common features (Fig. 9a,b). There are two velocity peaks pointing downwards near to the vessel walls and the vertical component of velocity acquires an almost null value between those peaks. The classical natural convection problem of a fluid contained in a differentially heated cavity shows two velocity peaks near to the walls but pointing to opposite directions when the Rayleigh number is high (106, Pr ¼ 0.71) [30]. This is because of the asymmetric boundary conditions applied to the vertical walls. In our case, there is an outer fluid taking heat from the water contained in the vessel so such as almost-symmetric behavior along X ¼ 0 is expected. The symmetry is not perfect because the fluid mechanics out of the vessel is not symmetric. The U-velocity component profiles along the vertical axis are shown in Fig. 9c and d. When Re ¼ 5,000, at


129

Fig. 6. Transient behavior of the streamlines around the vessel during the final instants of the cooling process of water with a forced convective flow of air (Re ¼ 22,000).

Fig. 7. Transient U and V velocity profiles along the horizontal axis at the middle of the duct for two Reynolds numbers.

Fig. 8. Transient behavior of the average friction coefficient calculated at the vessel walls during the cooling process of water with a forced convective flow of air (a. Re ¼ 5000; b. Re ¼ 22,000).

130


Fig. 9. U and V velocity profiles in the vessel during the cooling process of water with a forced convective flow of air.

t ¼ 20 s the U-velocity component along the vertical axis is low since temperature variations inside the vessel are not important enough to promote buoyancy. At higher times, t ¼ 400 s and 600 s, it can be noticed two velocity peaks pointing right and two to left, showing certain symmetry respect to the origin. When Re becomes higher, it can be observed an almost-symmetric behavior respect to the origin except for the higher U-velocity peak pointing left near to the top of the vessel. The U-velocity component of velocity along Yaxis is lower for Re ¼ 22,000 than for Re ¼ 5000.

During the cooling process the streamlines of the water contained tend to form two main vortexes along both vertical sides of the vessel (Fig. 10). When the Reynolds number is lower (Re ¼ 5000), two secondary vortexes in the middle of the vessel near to the top wall can be noticed during the early times of the cooling process. After that, both main vortexes prevail changing their form. When Re ¼ 22,000, the streamlines behavior is more complex. Both secondary vortexes are present during most of the cooling process being important even at high cooling times

Fig. 10. Transient behavior of the streamlines in the vessel during the cooling process of water.


131

Fig. 11. Transient behavior of dimensionless isotherms out of the vessel region at Y ¼ H/2 during the cooling process of water with a forced convective flow of air.

(t ¼ 400 s). At the end of the cooling process (t > 600 s), both secondary vortexes vanish and a new secondary vortex appears near to the left low corner. During the cooling process of water in a differentially heated cavity the streamlines adopt a different behavior (Fig. 2) of that observed at the end of the cooling process shows in Fig. 10. 5.2. Heat transfer The transient behavior of the isotherms out of the vessel region at Y]H/2 is shown in Fig. 11 for Re ¼ 5000 and Re ¼ 22,000. When Re ¼ 5,000, the isotherms show an oscillatory behavior at every selected time step. This behavior is explained by the vortex shedding from the east wall of the vessel. When the Reynolds number becomes higher (Re ¼ 22,000) this oscillatory behavior is not noticeable because vortex-shedding frequency becomes

higher, moreover the temperature gradient at the east wall is higher and the temperature profile becomes almost a step-shaped function. The isotherms of the water contained in the inner vessel are shown in Figs. 12 and 13 for Re ¼ 5000 and Re ¼ 22,000, respectively. In both cases, it can be observed that the thin walls of the vessel made of stainless steel offer a lower resistance to heat transfer than the inner fluid. Therefore the heat is extracted first from the walls and then from the left lower corner of the vessel. For Re ¼ 5,000, during most of the cooling process isotherms exhibit a similar behavior characterized by hotter regions which rise to the top of the vessel maintaining a characteristic shape (Re ¼ 5000: 200e400 s; Re ¼ 22,000: 500e600 s) until they are split (Re ¼ 5000: 500 s; Re ¼ 22,000: 400 s). After t ¼ 1100 s and t ¼ 650 s for Re ¼ 5000 and Re ¼ 22,000, respectively, there is a change in tendency when isotherms q ¼ 0.40 and q ¼ 0.45 appears in the

Fig. 12. Transient behavior of dimensionless isotherms in the vessel during the cooling process of water with a forced convective flow of air (Re ¼ 5000).

132


Fig. 13. Transient behavior of dimensionless isotherms in the vessel during the cooling process of water with a forced convective flow of air (Re ¼ 22,000).

Fig. 14. Transient behavior of the dimensionless temperature at different positions (X1: -4/11; 0, X2: 4/11; 0, Y1: 0; 4/11, Y2: 0; 4/11) during the cooling process of water with a forced convective flow of air.

inner vessel. It is in this temperature interval where density inversion of water takes place. The cooling rate curves shown in Fig. 14a and b represent the transient response of the dimensionless temperature at five referential positions in the vessel: the center of the vessel, which coincides with the origin of the coordinate axis system, two points

along the x-axis X1 (4/11; 0), X2 (4/11; 0) and two more along the y-axis Y1 (0; 4/11), Y2 (0; 4/11). Taking as reference the temperature at the center of the vessel and at X1, the cooling rate is increased 73% and 79% when the Reynolds number of the external flow of air is increased to 22,000, respectively. The lack of symmetry of the cooling process can be observed noticing the gap

Fig. 15. Transient behavior of the dimensionless temperature in the inner vessel.


133

Fig. 16. Transient behavior of the Nusselt number calculated at the outer west (a, b) and east (c, d) vessel walls during the cooling process of water with a forced convective flow of air.

between the corresponding curves X1, X2 and Y1, Y2. In a process with symmetry, these curves would be overlapping. When the Reynolds number is increased such a gap gets broader. Fig. 15 depicts the transient behavior of dimensionless temperature along the vertical axis at X ¼ 0, for Re ¼ 5000 (Fig. 15a) and Re ¼ 22,000 (Fig. 15b). Even the shape of the temperature profiles are similar for both Reynolds numbers, it can be observed how at the same time the temperature profile is lower when Re ¼ 22,000 (Fig. 15b). In both cases, it can be noticed how the temperature gradient becomes higher at the north wall of the vessel as a consequence of buoyancy. To account for the relative importance of convective heat transfer, the transient behavior of the local Nusselt number was calculated at the outer surfaces of the vessel. At the west wall, where the vessel is cooled directly by the airflow, the local Nusselt number is higher than in the other surfaces (Fig. 16a,b). The local Nusselt number (Nuw) along the west wall is described by a parabolic curve with a minimum just in the center of the wall where a stagnation point is present. When Re ¼ 5,000, Nuw calculated at the east wall (Fig. 16c) increases from a lower value (Nuw z 5) at the vessel bottom to a higher value (Nuw z 30) at the vessel top. Convection reigns over conduction as long as Y is increased. The behavior of Nuw is practically constant along the east side of the cavity when Re ¼ 22,000, but there are important variations close to the vessel corners (Fig. 16d). When Re ¼ 5,000, four main regions of the Nuw behavior at the north wall of the vessel can be noticed (Fig. 17a). First, Nuw decreases steeply near to the left corner, where flow velocity decreases abruptly. Then, Nuw rises in the region where two vortexes are attached to the cavity walls, as it can be seen in Fig. 18, later Nuw decreases where the vortex is shed to the west wall, and finally Nuw increases near to the right corner of the vessel. At the south wall, three regions can be distinguished (Fig. 17c). First Nuw decreases steeply near to the left corner, then, it

rises slowly where the fluid is recirculating, and finally Nuw decreases where streamlines are detached from the wall surface. At a higher Reynolds number, the variation of Nuw along the south (Fig. 17d) and north (Fig. 17b) walls are quite similar, Nuw decreases steeply near to the left corner of the wall then, there is a sudden change observed like a peak, Nuw increases along the wall and finally decreases at the right corner of the wall. 6. Conclusions The study of a heat transfer process through a conjugate model allows a more realistic characterization of fluid mechanics and heat transfer in the domain of interest. Instead of imposing typical boundary conditions, difficult to find for turbulent and steady heat convective flows, local heat flux and friction factors can be obtained at each wall, after solutions for velocity and temperature distributions have been found. Nusselt number and friction factor calculations require the use of implicit second order error expressions and a higher concentration of nodes near to the vessel walls for the calculation of gradients. During earlier times of the cooling process, outer airflow dynamics governs the convergence of the FVM simulations, and hence a lower time step is required (Dt ¼ 103 s) for the conjugate process in comparison with the non-conjugate natural convection process of water (Dt ¼ 0.5 s). The cooling rate of the water contained in the vessel is improved by up to 79% when the Reynolds number is increased from 5000 to 22,000. The lower thermal resistance offered by the vessel material controls the concentric shape of the isotherms, since first the walls of the vessel become an isothermal surface and then the heat is extracted almost uniformly from the water contained in it. The isotherms of water inside the vessel is similar for R ¼ 5000 and Re ¼ 22,000. During most of the cooling process, the isotherms

134


Fig. 17. Transient behavior of the Nusselt number calculated at the outer north (a, b) and south (c, d) vessel walls during the cooling process of water with a forced convective flow of air.

exhibit a similar behavior for both Reynolds numbers that is characterized by hotter regions, which rise to the top of the vessel maintaining a characteristic shape until they are split at the top of the vessel. The behavior of fluid mechanics in the water contained in the vessel changes when density inversion of water takes place. The local Nusselt number calculated at the outer surfaces of the vessel shows a higher prevalence of convection over conduction, mainly at the west wall, where fluid impacts the vessel directly. The local Nusselt number variation at the north and south walls, especially at the north wall where an inflexion point can be observed, is more important for Re ¼ 5000 than for Re ¼ 22,000. This is explained by the higher vortex shedding frequency appearing when Re ¼ 22,000.

Even it is simpler to implement a third order boundary condition, with a constant convective heat transfer coefficient, at the vessel walls, it is found, that the Nusselt number, and therefore the heat transfer coefficient, is highly dependent of time and position, therefore to assume a constant third order boundary condition could lead to wrong predictions of fluid mechanics and heat transfer in the inner fluid. Acknowledgments D.A. Vasco acknowledges to CONICYT-Chile for the support received by FONDECYT project 11130168. Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02)

References

Fig. 18. Streamlines behavior around the vessel for t ¼ 600 s and Re ¼ 22,000.

[1] S. Dorfman, Conjugate Problems in Convective Heat Transfer, CRC Press, Nueva York, 2010. [2] C. Lamnatou, E. Papanicolaou, V. Belessiotis, N. Kyriakis, Finite-volume modeling of heat and mass transfer during convective drying of porous bodies e non-conjugate and conjugate formulations involving the aerodynamic effects, Renew. Energy 35 (7) (2010) 1391e1402. [3] D. Kim, R. Viskanta, A laser-doppler study of effects of wall conductance on natural convection in differently oriented square cavities, J. Fluid Mech. 144 (1) (1984) 153e176. [4] A. Liaqat, A. Baytas, Conjugate natural convection in a square enclosure containing volumetric sources, Int. J. Heat Mass Transf. 44 (17) (2001) 3273e3280. [5] J. House, C. Beckermann, T. Smith, Effect of a centered conducting body on a natural convection heat transfer in an enclosure, Numer. Heat Transf. Part A 18 (2) (1990) 213e225. [6] M. Yeong-Ha, M. Jung-Jung, Y. Soo-Kim, Numerical study on transient heat transfer and fluid flow of natural convection in an enclosure with a heatgenerating conducting body, Numer. Heat Transf. Part A 35 (4) (1999) 415e433.

D.A. Vasco et al. / International Journal of Thermal Sciences 87 (2015) 121e135 [7] M. Yeong-Ha, M. Jung-Jung, A numerical study on three dimensional conjugate heat transfer of natural convection and conduction in a differentially heated cubic enclosure with a heat-generating cubic conducting body, Int. J. Heat Mass Transf. 43 (23) (2000) 4229e4248. [8] M. Kumar-Das, S. Kumar-Reddy, Conjugate natural convection heat transfer in an inclined square cavity containing a conducting block, Int. J. Heat Mass Transf. 49 (25e26) (2006) 4987e5000. [9] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington DC, 1980. [10] X. Chen, P. Han, A note on the solution of conjugate heat transfer problems using simple-like algorithms, Int. J. Heat Fluid Flow. 21 (4) (2000) 463e467. [11] G. Kuznetsov, M. Sheremet, Conjugate natural convection in an enclosure with local heat sources, Comput. Therm. Sci. 1 (3) (2009) 341e360. [12] G. Kuznetsov, M. Sheremet, Conjugate heat transfer in an enclosure under the condition of internal mass transfer and in the presence of the local heat source, Int. J. Heat Mass Transf. 52 (1e2) (2009) 1e8. [13] G. Kuznetsov, M. Sheremet, Conjugate natural convection in an enclosure with a heat source of constant heat transfer rate, Int. J. Heat Mass Transf. 54 (1e3) (2011) 260e268. [14] B. Lee, The effect of turbulence on the surface pressure field of a square prism, J. Fluid Mech. 69 (2) (1975) 263e282. [15] D. Lyn, W. Rodi, The flapping shear layer formed by the flow separation from the forward corner of a square cylinder, J. Fluid Mech. 267 (1) (1994) 353e376. [16] D. Lyn, S. Einva, W. Rodi, J. Park, A laser-doppler velocimetry study of ensemble averaged characteristics of the turbulence near wake of a square Cylinder, J. Fluid Mech. 304 (1) (1995) 285e319. [17] H. Kawamura, N. Kawashima, An application of a near-wall keε model to the turbulence flow with transpiration and to the oscullatory flow around a square cylinder, in: Second International Symposium on Turbulence, Heat and Mass Transfer, Delft, Holland, 1997, pp. 379e388. [18] N. Moraga, J. Riquelme, L. Jauriat, Unsteady conjugate water/air mixed convection in a square cavity, Int. J. Heat Mass Transf. 52 (23e24) (2009) 5512e5524.

135

[19] N.O. Moraga, M. Andrade, D.A. Vasco, Unsteady conjugate mixed convection phase change of a power-law non newtonian fluid in a square cavity, Int. J. Heat Mass Transf. 53 (15e16) (2010) 3308e3318. [20] N.O. Moraga, R.A. Lemus-Mondaca, Numerical conjugate air mixed convection/ non-newtonian liquid solidification for various cavity configurations and rheological models, Int. J. Heat Mass Transf. 54 (23e24) (2011) 5116e5125. [21] R.A. Lemus-Mondaca, N.O. Moraga, J. Riquelme, Unsteady 2D conjugate natural non-newtonian convection with non-newtonian liquid sterilization in square cavity, Int. J. Heat Mass Transf. 61 (2013) 73e81. [22] M. Raisee, A. Jafari, H. Babaei, H. Iacovides, Two-dimensional prediction of time dependent, turbulent flow around a square cylinder confined in a channel, Int. J. Numer. Methods Fluids 62 (11) (2010) 1232e1263. [23] B.E. Launder, B.I. Sharma, Application of the energy dissipation model of turbulence to the calculation of flow near a spinning disc, Lett. Heat Mass Transf. 1 (2) (1974) 131e138. [24] B.E. Launder, D.B. Spalding, Mathematical Models of Turbulence, Academic Press, Londres, 1972. [25] C.J. Chen, S.Y. Jaw, Fundamentals of Turbulence Modeling, Taylor and Francis, Washington, 1998. [26] M. Giangi, T.A. Kowalewski, F. Stella, E. Leonardi, Natural convection during ice formation: numerical simulation vs. experimental results, CAMES 7 (2000) 321e342. [27] A. Zografos, W. Martin, J. Sunderland, Equation of properties as a function of temperature for seven fluids, J. Comput. Methods Appl. Mech. Eng. 61 (2) (1987) 477e493. [28] S.V. Patankar, D.B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Int. J. Heat Mass Transf. 15 (10) (1972) 1787e1806. [29] J. Banaszek, Y. Jaluria, T.A. Kowalewski, M. Rebow, Semi-implicit FEM analysis of natural convection in freezing water, Numer. Heat Transf. Part A 36 (5) (1999) 449e472. [30] T. Fusegi, J.M. Hyun, K. Kuwahara, B. Farouk, A numerical study of threedimensional natural convection in a differentially heat cubical enclosure, Int. J. Heat Mass Transf. 34 (6) (1991) 1543e1557.