- Email: [email protected]
Direct numerical simulation of Rayleigh-Bénard-Marangoni convection of cold water near its density maximum in a cylindrical pool
International Journal of Thermal Sciences 124 (2018) 327–337
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Direct numerical simulation of Rayleigh-Bénard-Marangoni convection of cold water near its density maximum in a cylindrical pool
MARK
You-Rong Lia,∗, Lu Jiaa, Li Zhanga, Yu-Peng Hub, Jia-Jia Yua a Key Laboratory of Low-grade Energy Utilization Technologies and Systems of Ministry of Education, College of Power Engineering, Chongqing University, Chongqing 400044, China b Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang 621999, China
A R T I C L E I N F O
A B S T R A C T
Keywords: Rayleigh-Bénard-Marangoni convection Density inversion Flow pattern Heat transfer Cylindrical pool
In order to understand the influence of Marangoni effect on Rayleigh-Bénard convection, a series of threedimensional numerical simulations of Rayleigh-Bénard- Marangoni convection of cold water near its maximum density in a cylindrical pool heated from below and with a free surface at the top were carried out by using the finite volume approach. Results show that the flow pattern is strongly influenced by the Rayleigh number, as well as by the Biot number. When the Biot number is small, the critical Rayleigh number for the onset of convection is low due to the strong Marangoni effect on the free surface. As the Biot number increases, the Marangoni effect becomes weaker. Therefore, the critical Rayleigh number becomes much larger and more flow patterns can be observed at different density inversion parameters. Furthermore, the average Nusselt number is influenced by the Rayleigh number, the Biot number and the density inversion parameter. Specifically, the average Nusselt number increases with the increase of the Rayleigh number and the decrease of the density inversion parameter. Meanwhile, as the Biot number goes up, the average Nusselt number increases at first and then decreases. In addition, there exists a sudden increase of the Nusselt number when the flow pattern transits from a relative stable state to a violent state, which means the flow pattern can also influence the heat transfer ability.
1. Introduction For a cavity heated from the bottom and with a free surface at the top, when the temperature difference between the bottom and top is large enough, the flow driven by the buoyancy and the surface tension gradient will happen, which is named as Rayleigh- Bénard-Marangoni (R-B-M) convection [1]. Because of its wide existence in industrial applications, such as cold storage air conditioning, crystal growth, welding process etc. [2–4], this convection has been widely concerned, especially in recent decades. The pioneering work on Bénard convection was carried out by Bénard [5]. According to his experiment, it was found that the fluid keeps stationary at a small temperature difference; however, the hexagon flow pattern appears when the temperature difference exceeds a certain value. After that, Block [6] pointed out that the hexagonal cell pattern may be caused by the surface tension gradient. Ondarcuhu et al. [7] investigated the effect of the temperature difference on the flow pattern of Bénard-Marangoni convection in a rectangular cavity by experiments. Results indicated that when the temperature difference is small, the flow pattern appears as the pentagon and square patterns. ∗
With the gradual increase of the temperature difference, the pentagon pattern disappears at first, and then reproduces. Dijkstra [8] numerically analyzed the stability of two-dimensional steady flow patterns in rectangular containers, and obtained the bifurcation diagram of flow patterns. Similar studies of a volatile fluid in a sealed and confined rectangular container were carried out by Qin et al. [9,10]. Then, Xu et al. [11] studied R-B-M convection in a cylinder heated laterally and cooled at top by using higher order finite difference scheme. They found that the critical Grashof number increases with the increase of the Marangoni number. Nezar and Rahal [12] numerically investigated the effect of gravity and surface tension gradient on the flow patterns. It was concluded that with the increase of the thickness of liquid layer, the flow pattern would bifurcate from hexagonal cells to stationary longitudinal rolls. Rahal et al. [13] conducted a series of experiments to study the effect of the Marangoni number, Biot number and Prandtl number on the free surface deformation, and found that the flow will bifurcate to a time-dependent flow when the Biot number is large enough. Sakhy et al. [14] carried out a series of three-dimensional numerical simulations on R-B-M convection in a cylindrical container heated from below by a non-uniform heat flux, the results showed that
Corresponding author. E-mail address: [email protected] (Y.-R. Li).
http://dx.doi.org/10.1016/j.ijthermalsci.2017.10.034 Received 13 May 2017; Received in revised form 5 September 2017; Accepted 26 October 2017 1290-0729/ © 2017 Elsevier Masson SAS. All rights reserved.
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
γ γT λ Г μ ν θ Θ ρ σ τ
Nomenclature Bi d ez g h Ma Nu P Pr q r R Ra T V z Z
Biot number depth of cylindrical pool, m z-directional unit vector gravitational acceleration, m/s2 heat exchange coefficient on the free surface, W/(m2·K) Marangoni number Nusselt number dimensionless pressure Prandtl number exponent in Eq. (1) radial coordinate or radius, m dimensionless radius Rayleigh number temperature, °C dimensionless velocity vector axial coordinate dimensionless axial coordinate
Subscripts 0 ave cri e h m t
Greek symbols α
thermal expansion coefficient, (°C)-q surface tension coefficient, N/(m·K) thermal conductivity, W/(m·K) aspect ratio dynamic viscosity, kg/(m·s) kinematic viscosity, m2/s azimuthal coordinate, rad dimensionless temperature density, kg/m3 surface tension, N/m dimensionless time
initial average critical environment hot wall density inversion point top free surface
thermal diffusivity, m2/s
the flow patterns depend on the Rayleigh number, Marangoni number, Biot number and the ratio of thermal conductivities. Boeck [15] investigated the effect of homogeneous magnetic field on Bénard convection in a rectangular domain, and the transition from stationary hexagonal pattern to time-dependent convection was explored for the fluid with zero Prandtl number. Guo and Narayanan [16] performed a linear stability analysis of R-B-M convection of a silicone oil in the annular cavity. It was found that the depth of the cavity rather than the width of the gap has an important effect on the Marangoni convection. Furthermore, the critical wavelength generally decreases with the decrease of the scaled gap. For all reports mentioned above, the density of working fluid was supposed to be a linear function of temperature. However, cold water is a typical fluid which has the maximum density at the temperature around 4 °C, this special property makes the flow become more complicated [17]. Sasaguchi et al. [18] carried out numerical analyses on the influence of eccentricity on natural convection of cold water in a rectangular cavity and found that when the initial water temperature is 4 °C, the relative positon of the cylinder and rectangular cavity has a significant effect on the average Nusselt number. When the initial water temperature exceeds 4 °C,a downward flow pattern can be observed at first, and then changes to upward flow pattern gradually. Osorio et al. [19] numerically simulated natural convection of cold water in a square cavity and discussed the effect of the density inversion parameter on the flow pattern. Furthermore, Nishimura et al. [20] investigated ReB convection of cold water in a rectangular enclosure whose aspect ratio was 1.25. Results showed that the flow pattern and the temperature field are symmetric; however, a sinking cell is formed at the center of the enclosure because of the influence of the density inversion. Lin and Nansteel [21] simulated natural convection of cold water in a square enclosure heated and cooled on opposite vertical walls, and found that there is a multicellular flow pattern. Furthermore, the flow pattern is independent on the Rayleigh number at 103≤Ra≤106. Yuan and Li et al. [22,23] performed a series of numerical simulations on natural convection of cold water in the horizontal annulus and the horizontal eccentric annulus, and discussed the effect of Rayleigh number and density inversion parameter on the flow patterns. Li et al. [24] also found some new flow patterns in ReB convection of cold water at a fixed aspect ratio. Hu et al. [25] performed numerical simulations of cold water in complex configurations and obtained the heat transfer
correlations in different cavities. However, when both the buoyancy and the surface tension gradient are taken into account, the formation of the flow pattern is very complicated and different. Therefore, the aim of present work is to obtain various flow patterns of R-B-M convection of cold water in a cylindrical container. Furthermore, we also focus on the effects of Rayleigh number, Biot number and density inversion parameter on the heat transfer and evolution of flow patterns. 2. Model formulation 2.1. Basic assumption and governing equations The physical model is schematically shown in Fig. 1. We consider a cylindrical pool with the aspect ratio Г = r0/d = 1.5, which is filled with cold water near its maximum density. Where, r0 is the radius and d is the depth of the cylindrical pool. The sidewall is regarded as perfectly adiabatic. The bottom of the pool is kept at constant temperature Th, while the top is a free surface with a variation temperature. The existence of free surface means that there is the heat exchange between cold water and the environment. The heat transfer coefficient and the
Fig. 1. Physical model and the coordinate system.
328
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
environment temperature are marked as h and Te, respectively. In order to simplify, some assumptions are introduced in the present model: (1) The no-slip boundary condition on all solid-liquid boundaries is applied. (2) The cold water is incompressible Newtonian fluid and the flow is laminar. (3) Viscous dissipation is neglected. (4) The physical properties are taken as constant except for the density ρ and surface tension σ. Furthermore, the density is assumed to be a nonlinear function of temperature [26] as
ρ (T ) = ρm (1 − γ T − Tm q )
Table 1 Mesh dependence test at Bi = 5, Θm = 0.5 and Ra = 104. Mesh 20R 40R 60R 80R
4.1958 4.2309 4.2472 4.2527
1.34 0.51 0.13 –
Table 2 Comparison of average Nusselt number at different Rayleigh numbers.
where ρm = 999.972 kg/m is the maximum density at temperature Tm = 4.029325 °C, the thermal expansion coefficient γ = 9.297173 × 10−6(°C)−q, q = 1.894816. The surface tension can be expressed as
σ = σ0 − γT (T − T0)
Ra
(2)
0.5 1.0 1.5 2.0
where σ0 = σ(T0), γT = −∂σ/∂T, and T0 is the initial temperature. Using d, d2/ν, ν/d and μν/d2 as scale quantities for length, time, velocity and pressure, the dimensionless governing equations in the cylindrical coordinate system can be expressed as follows:
Nuave
× × × ×
4
10 104 104 104
Deviation (%)
Present
Ref. [27]
2.11 2.33 2.45 2.56
2.07 2.29 2.41 2.53
1.89 1.71 1.63 1.17
(3)
∇⋅V =0
∂τ V+V ⋅∇V= − ∇P +
∇2 V+ (Ra/Pr )
∂τ Θ + V ⋅∇Θ = (1/ Pr ) ∇2 Θ
Θ−Θm q eZ
2.2. Boundary conditions (4) Under above assumptions, the boundary conditions are: At the bottom wall (0 ≤ R < Г, 0≤θ < 2π, Z = 0):
(5)
where V(VR, Vθ, VZ) is the dimensionless velocity vector, Θ=(T-Tt)/(ThTt) the dimensionless temperature. Θm=(Tm-Tt)/(Th-Tt) is the density inversion parameter which defines the location of the maximum density temperature within fluid layer. ez is the unit vector in the Z direction. Meanwhile, two important dimensionless parameters, Ra and Pr, are defined respectively as: Ra = gγd3(Th−Tt)q/(αν), Pr=ν/α
Vθ = VR = VZ = 0, Θ = 1
Bi z T (z ) = (Th − Te ) + Th 1 + Bi d
Vθ = VR = VZ = 0, ∂R Θ = 0. At the free surface (0 ≤ R < Г, 0≤θ < 2π, Z = 1):
(6a-b)
(12a-d) where Ma is the Marangoni number, which is defined as:
Ma =
γT d (Th − Tt ) . μα
(13)
The average Nusselt number on the bottom wall is defined as (7)
Nuave =
1 πΓ 2
2π
∫0 ∫0
Γ
∂Θ RdRdθ. ∂Ζ
(14)
2.3. Numerical procedure and validation (8a-b) The non-uniform computational mesh is shown in Fig. 2, the finer grids cluster close to the boundaries for the reason that the temperature gradient and velocity gradient are large. Finite volume approach is applied to discretize the fundamental governing equations. The convective terms and the diffusion terms are discretized through using the QUICK scheme and the second-order central difference scheme, respectively. The SIMPLE algorithm is adopted to handle the pressure-
It is worth pointing out that the temperature difference ΔT is defined as ΔT = Th−Tt rather than ΔT = TheTe in this paper for the reason that the attention is the flow of cold water in the pool. Therefore, according to Eq. (8), the temperature difference can be expressed as
Bi (Th − Te ) 1 + Bi
(11a-b)
VZ = 0, ∂Z Θ + BiΘ + 1 = 0, ∂Z VR = −Ma∂R Θ , ∂Z Vθ = −(Ma/ R) ∂θ Θ,
where Bi is the Biot number, which is defined as Bi = hd/λ. Therefore, the temperature at the free surface is
Bi (Th − Te ) + Th 1 + Bi
(10a-b)
At the cylinder sidewall (R = Г, 0≤θ < 2π, 0 ≤ Z ≤ 1):
At a small temperature difference, as considered by Touihri et al. [27], the imposed thermal conditions lead to a conductive state without any flow. In this case, the temperature profile is linear along the vertical direction,
ΔT =
Deviation (%)
(1) 3
Tt = T (d ) = −
80θ × 10Z 120θ × 20Z 160θ × 30Z 200θ × 40Z
× × × ×
Nu
(9)
Fig. 2. The cylindrical mesh model, R-θ plane (left), R-Z plane (right).
329
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 3. Contours of the axial velocity in the Z = 1/2 plane at Bi = 5 and Θm = 0.5. Solid lines correspond to positive value; dashed lines correspond to negative value.
Fig. 4. Flow patterns obtained from different routes.
Table 3 Stability range of all stable convective branches at Bi = 5 and Θm = 0.5. Flow pattern
Corresponding figure
Range of stability
upward one-torus downward two-roll single-roll upward two-roll four-spoke dipole
Fig. Fig. Fig. Fig. Fig. Fig.
140 ≤ Ra≤650 160 ≤ Ra≤104 0.1 × 104≤Ra≤4.0 × 104 0.8 × 104≤Ra≤4.0 × 104 1.2 × 104≤Ra≤4.0 × 104 160 ≤ Ra≤500
3(a) 3(b) 3(c) 3(d) 4(a) 3(b)
Fig. 6. The variation of the average Nusselt number with Rayleigh number at Bi = 5 and Θm = 0.5.
velocity coupling, and the fully implicit scheme is used for the temporal discretization. The dimensionless time step is 10−4. At each time step, the solution is considered to be convergent when the residual errors of both the momentum and energy equations are less than 10−5. In order to verify the independence of the grid, a set of numerical simulations with four different meshes were performed at Bi = 5, Θm = 0.5 and Ra = 104. The average Nusselt numbers on the bottom wall are listed in Table 1. It can be found that the deviation of the average Nusselt numbers is only 0.13% between two denser grids. Comparing with results for grid independence in Refs. [29,34], this deviation can be completely negligible. Therefore, the grid of 60R × 160θ × 30Z is an appropriate choice considering both the
Fig. 5. Schematic diagram of stability range and the flow pattern transition at Bi = 5.
330
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 7. Flow pattern obtained from the conductive state at Θm = 0.1 and Bi = 100. Upper plot: contours of the axial velocity in Z = 1/2 plane; solid lines correspond to positive value and dashed lines correspond to negative value. Lower plot: isothermal surfaces of Θ = 0.3, 0.5 and 0.7.
moment, the hot-rising flow is located in the center of the cylindrical pool while the cold-descending flow is near the sidewall. Therefore, the flow pattern is named as upward one-torus. When the Rayleigh number is between 550 and 2340, the flow transits to an axisymmetric structure, which is named as downward two-roll. This flow pattern has a cold-descending flow in the center and two hot-rising flows that are distributed symmetrically on both sides, as shown in Fig. 3(b). When 2.35 × 103≤Ra≤2.17 × 104, a single-roll pattern with two approximate symmetrical sectors is observed, as shown in Fig. 3(c). When the Rayleigh number exceeds 2.18 × 104, the flow evolves into an upward two-roll, as shown in Fig. 3(d), which has the opposite flow direction with the flow pattern observed at 550 ≤ Ra≤2340. Li et al. [29] had determined that the critical Rayleigh number of RB convection onset for cold water at Θm = 0.5 is Racri = 1.73 × 104, which is much larger than that of R-B-M convection onset in a cylindrical pool with the free surface at the same condition. It hints that the top solid wall can remarkably enhance the stability of the fluid layer. On the other hand, the nonuniform temperature distribution on the free surface leads to the surface tension gradient. This surface tension gradient depends on the Biot number and is very large at small Biot number. Meanwhile, it can help to overcome the stabilization effect of viscosity, as considered by Block [6]. As a result, the critical Rayleigh number of R-B-M convection is smaller than that of R-B convection.
numerical accuracy and computational expenses. To validate the numerical method, some simulations were performed at the same condition with Ref. [27]. When Bi = 100 and Ma = 0, the comparison of the average Nusselt numbers is given in Table 2. Obviously, the obtained results are very close to those in the reference and the maximum deviation is less than 1.9%. Furthermore, the same flow pattern is also observed. As a consequence, the numerical method is accurate enough to simulate Rayleigh-Bénard-Marangoni convection of cold water near its maximum density in a cylindrical pool.
3. Results and discussion 3.1. Effect of Rayleigh number Firstly, we analyze the effect of the Rayleigh number on the flow pattern and the heat transfer characteristics in the case of Bi = 5. A series of numerical simulations at the Rayleigh number range from 100 to 4.0 × 104 were carried out when the conductive solution was used as the initial condition. As shown in Fig. 3, four typical flow patterns orderly appear with the increase of the Rayleigh number. When Ra < 160 at Θm = 0.5, the fluid keeps stationary and no flow appears in the cylindrical pool. When Rayleigh number increases to the critical value of Racri = 160, the fluid begins to destabilize. The conductive solution bifurcates to a steady one-torus flow pattern, which is the same with the centrosymmetric flow pattern observed by Ma et al. [28]. At this
3.1.1. Evolution from upward one-torus pattern The simulation result shows this flow pattern remains stable at 331
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 8. Evolution processes of various flow pattern.
Fig. 9. Flow pattern obtained from the conductive state at Θm = 0.5 and Bi = 100. Upper plot: contours of the axial velocity in Z = 1/2 plane; solid lines correspond to positive value and dashed lines correspond to negative value. Lower plot: isothermal surfaces of Θ = 0.3, 0.5 and 0.7.
140 ≤ Ra≤650. However, the critical Rayleigh number of the flow transition from the one-torus pattern to the conductive state is smaller than that of the flow transition from the conductive state to the onetorus pattern. Therefore, there exists the transition hysteresis of the flow pattern, which is also found by Buell and Catton [30]. When the Rayleigh number is between 700 and 2850, the upward one-torus pattern turns into a downward two-roll pattern. Furthermore, it bifurcates to a single-roll pattern at 0.3 × 104≤Ra≤2.3 × 104. Finally, when the Rayleigh number is above 2.4 × 104, the upward two-roll pattern appears.
3.1.2. Evolution from downward two-roll pattern When the downward two-roll pattern converged at Ra = 1.5 × 103 is used as an initial condition, it's found that this flow pattern is stable at 160 ≤ Ra≤ × 104. When the Rayleigh number is below 160, it decays to the conductive state. When the Rayleigh number is beyond 1.2 × 104, this pattern bifurcates to a new flow pattern that is called as four-spoke pattern, as shown in Fig. 4(a).
3.1.3. Evolution from single-roll pattern The single-roll pattern remains stable at a wide range of 332
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 10. Evolution from flabellum pattern to two-tori pattern at Ra = 8.0 × 104. Fig. 11. Flow pattern obtained from the conductive state at Θm = 0.7 and Bi = 100. Upper plot: contours of the axial velocity at Z = 1/3 plane; solid lines correspond to positive value and dashed lines correspond to negative value. Lower plot: isothermal surfaces of Θ = 0.3, 0.5 and 0.7.
rather than one-torus pattern. In summary, because of the strong surface tension gradient at Bi = 5, the flow patterns are relatively simple and mainly affected by the Rayleigh number. Specifically, four distinct flow patterns can be obtained when the conductive solution is used as the initial condition, including upward one-torus, downward two-roll, single-roll and upward two-roll. Furthermore, two new flow patterns, four-spoke and dipole patterns, are also observed when the different flow states are used as initial conditions. The stability range of the different flow patterns and schematic diagram on flow pattern transitions are given in Table 3 and Fig. 5, respectively. It can be found that there exist a variety of the flow patterns at the same Rayleigh number and the transition of the flow patterns is abundant. The variation of the average Nusselt number on the bottom with the Rayleigh number is shown in Fig. 6. It can be found that the average Nusselt number strongly depends on the Rayleigh number. Specifically, for all the flow patterns, the average Nusselt number will increase simultaneously with the Rayleigh number. Besides, at a constant Rayleigh number, the average Nusselt number of different flow patterns is also different, which means that the flow pattern can also influence heat transfer ability. With the increase of the Biot number, the free surface temperature Tt approaches gradually the environment temperature Te. Therefore, the effect of the thermocapillary force on the free surface becomes weaker. When the density inversion parameter is small, the density gradient direction in the fluid is almost same as that of the common fluid. The simulation results at Θm = 0.1 and Bi = 100 are displayed in Fig. 7. The critical Rayleigh number for the convection onset is 2750 if the conductive solution is used as the initial condition. When the Rayleigh number is over the threshold, the upward one-torus state appears first. At this time, the intensity of flow is very weak, as shown
Table 4 Stability range of all stable convective branches. Θm
Flow pattern
Corresponding figure
Range of stability
0.1
upward one-torus sadman four-spoke X two-roll downward one-torus one-torus dice
Fig. 7(a) Fig. 7(b) Fig. 7(c) Fig. 7(d) Figs. 7(e) and 8(b) Fig. 8(a)
2.75 × 103≤Ra≤0.55 × 104 0.3 × 104≤Ra≤0.85 × 104 0.6 × 104≤Ra≤6.8 × 104 4.7 × 104≤Ra≤8.0 × 104 0.55 × 104≤Ra≤8.0 × 104 2.0 × 104≤Ra≤8.0 × 104
Fig. 8(c)
0.3 × 104≤Ra≤0.58 × 104
0.5
two-tori three-spot flabellum
Fig. 9(a) and (d) Fig. 9(b) Fig. 9(c)
1.84 × 104≤Ra≤1.25 × 105 1.85 × 104≤Ra≤6.5 × 104 1.85 × 104≤Ra≤7.8 × 104
0.7
two-tori seven-cell eight-cell
Fig. 11(a) and (d) Fig. 11(b) Fig. 11(c)
1.75 × 105≤Ra≤3.5 × 105 1.83 × 105≤Ra≤3.5 × 105 1.9 × 105≤Ra≤3.5 × 105
103≤Ra≤4 × 104 when the flow state at Ra = 0.5 × 104 is used as an initial condition. When the Rayleigh number varies from 150 to 900, it transits into the downward two-roll pattern. If Ra≤150, the system decays to the conductive state. 3.1.4. Evolution from upward two-roll pattern In this case, the simulation results show that the upward two-roll pattern exists at a wide range of the Rayleigh number from 0.8 × 104 to 4 × 104. When Ra < 0.7 × 104, it bifurcates to the single-roll pattern. If the Rayleigh number is less than 500, the pattern turns into a new flow state that is named as the dipole pattern, as depicted in Fig. 4(b), 333
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 12. Schematic diagram of stability range and the flow pattern transition at Bi = 100.
flow pattern is also detected by Borońska and Tuckerman [31,32], which is a secondary bifurcation flow pattern from the star pattern. When 0.66 × 104≤Ra≤7 × 104, the flow bifurcates into the spoke
in Fig. 7(a). When the Rayleigh number is between 3200 and 6500, the flow bifurcates to a new flow pattern, which looks like a sadman. Therefore, it is called as the sadman pattern, as shown in Fig. 7(b). This 334
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
pattern. When the Rayleigh number is continually increased to 2.0 × 104≤Ra≤8.0 × 104, the upward one-torus pattern evolves into the downward one-torus pattern, and the evolution process is displayed in Fig. 8(a). 3.1.6. Evolution from sadman pattern When the sadman pattern converged at Ra = 0.5 × 104 is used as the initial condition, it is found that this flow pattern is stable at 0.3 × 104≤Ra≤0.85 × 104. When the Rayleigh number increases from 0.85 × 104 to 4.0 × 104, it turns into the two-roll pattern, as shown in Fig. 8(b). When Ra > 4.0 × 104, the aberrant two-roll pattern is observed. 3.1.7. Evolution from spoke pattern The spoke patterns include the four-spoke pattern and the X pattern. The simulation result shows that the four-spoke pattern is stable in the range of 0.6 × 104≤Ra≤6.8 × 104. When Ra > 7.0 × 104, it evolves into the X pattern. On the contrary, the spoke pattern decays to a new dice pattern at 0.3 × 104≤Ra≤0.58 × 104, which is also observed by Hu et al. [33], and this evolution process is displayed in Fig. 8(c). The X pattern is stable at 4.7 × 104≤Ra≤8.0 × 104, and transits into the four-spoke pattern and the dice pattern with the decrease of the Rayleigh number. 3.1.8. Evolution from two-roll pattern The two-roll flow pattern exists at the range of the Rayleigh number between 0.56 × 104 and 8.0 × 104 when the result at Ra = 7.5 × 104 is used as the initial condition. With the decrease of the Rayleigh number, this flow pattern becomes more regular. At Ra < 0.55 × 104, the sadman pattern appears and it evolves into the upward one-torus pattern at Ra≤0.3 × 104. Therefore, there is no new flow pattern during this evolution. When the density inversion parameter is increased to Θm = 0.5 at Bi = 100, there are only four typical flow patterns when the conductive state is used as the initial condition, as shown in Fig. 9. In this case, the critical Rayleigh number for the convection onset increases to 1.84 × 104. When the Rayleigh number is over this value, the two-tori flow pattern appears firstly, as shown in Fig. 9(a). When the Rayleigh number varies from 2.2 × 104 to 3.4 × 104, the three-spot pattern is observed, which is consist of three hot-rising flows near the sidewall and a cold-descending flow in the center, as shown in Fig. 9(b). When the Rayleigh number is between 3.5 × 104 and 4.8 × 104, a new flow pattern stably exists, which looks like a flabellum. This flow pattern contains four flow cells near the surrounding wall and a hot-rising flow in the center of the cylindrical pool, as shown in Fig. 9(c). Furthermore, in a wide range of 5.0 × 104≤Ra≤1.25 × 105, the two-tori flow pattern appears again, but the flow becomes more intense, as shown in Fig. 9(d). Finally, the flow becomes oscillatory at Ra > 1.25 × 105. When the two-tori pattern is used as the initial condition, there is no flow pattern evolution regardless of the variation of the Rayleigh number, which means this flow pattern is relatively stable.
Fig. 13. The variation of the average Nusselt number with Rayleigh number.
pattern. When the Rayleigh number is small, it appears as the regular four-spoke pattern, like that in Fig. 7(c). However, with the further increase of the Rayleigh number, the region of the hot-rising flow becomes smaller, which leads to a ‘X’ flow pattern, as depicted in Fig. 7(d). When the Rayleigh number varies from 7.2 × 104 to 8.0 × 104, an aberrant two-roll pattern is observed, as shown in Fig. 7(e). At last, the flow becomes oscillatory state at Ra > 8.0 × 104.
3.1.9. Evolution from spot pattern In this case, it's found that the three-spot pattern is stable at 1.85 × 104≤Ra≤6.5 × 104. When the Rayleigh number exceeds 6.5 × 104, it turns into the two-tori pattern. On the other hand, when the result at Ra = 4.5 × 104 is used as the initial condition, the flabellum pattern keeps stable at 1.85 × 104≤Ra≤7.8 × 104. When Ra≥8.0 × 104, it evolves into the two-tori pattern, the evolution process as shown in Fig. 10. When Ra > 1.25 × 105, the flow becomes oscillatory. When the density inversion parameter is approaching 1, the position of the maximum density in the fluid layer moves to the bottom of the pool, which means that the flow occurs only in a thin liquid layer near the bottom wall. At Θm = 0.7, the distribution of axial velocity at Z = 1/3 plane and isothermal surfaces of Θ = 0.3, 0.5 and 0.7 are
3.1.5. Evolution from upward one-torus pattern By increasing or decreasing the Rayleigh number, the evolution from the upward one-torus flow pattern is exhibited when the flow state at Ra = 3 × 103 is used as the initial condition. It is found that this flow pattern remains stable at 2.75 × 103≤Ra≤5.5 × 103. When the Rayleigh number exceeds 0.6 × 104, it bifurcates into the sadman 335
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 14. Schematic diagram of flow patterns with the Biot number.
shown in Fig. 11(b). When 2.5 × 105≤Ra≤3.1 × 105, the seven-cell pattern evolves into the eight-cell pattern, like that in Fig. 11(c). Finally, when the Rayleigh number exceeds 3.2 × 105, the flow bifurcates into the two-tori pattern again, as shown in Fig. 11(d). When Θm = 0.7, there is no new flow pattern to be found when the primary bifurcation patterns are used as initial conditions. Since the similar secondary bifurcation results have been discussed at Θm = 0.5, it will not be repeated here. Comparing the results at Bi = 100 with those at Bi = 5, it is found that more flow patterns are observed at the large Biot number. The stability range of each flow pattern at different density inversion parameters is listed in Table 4. The schematic diagram of the flow pattern transition as a function of Rayleigh number is shown in Fig. 12, it can be found that the flow pattern becomes more abundant with the decrease of the density inversion parameter. In addition, the average Nusselt number as a function of the Rayleigh number is also depicted in Fig. 13. It is obvious that the average Nusselt number gets smaller as the density inversion parameter increases, which means the heat transfer ability becomes weaker.
3.2. Effect of Biot number When Bi = 0, the free surface can be regarded as perfectly adiabatic. In this case, there is no temperature difference between the bottom and the free surface because the sidewall is also adiabatic. Therefore, no matter how much the Rayleigh number or the density inversion parameter is, the fluid layer is always static. When Bi > 0, the heat exchange from the free surface leads to a vertical temperature gradient. If the Rayleigh number is large enough, the R-B-M convection will appear. Fig. 14 qualitatively shows the evolution of the flow pattern with the Biot number at different Rayleigh numbers. As shown in Fig. 14(a) at Ra = 103, with the increase of the Biot number, it's found that two flow patterns can be observed regardless of the density inversion parameter. In detail, when Θm = 0.1 and Bi ≤ 18, Θm = 0.5 and Bi ≤ 15 or Θm = 0.9 and Bi ≤ 11, the tworoll flow pattern first appears. However, when the Biot number is above each critical value, the two-roll pattern evolves to the one-torus pattern. Continuing to increase the Biot number, the flow pattern remains the same. When Θm = 0.1 and Bi ≥ 75, Θm = 0.5 and Bi ≥ 65 or Θm = 0.9 and Bi ≥ 60, there is no flow happening in the fluid layer. When Ra = 5 × 103, the simulation result is similar to that at Ra = 103. Specifically, at Θm = 0.1, when Bi ≤ 28, the flow state that bifurcates from the conductive state leads to a single-roll pattern. When the Biot number is between 29 and 75, the flow turns into the two-roll pattern. Finally, when Bi > 75, it appears as a sadman pattern and remains constant. At Θm = 0.5, the single-roll pattern and the two-roll pattern can be observed at Bi ≤ 42 and 42 < Bi ≤ 70, respectively. When Bi > 70, the fluid layer becomes static. When Θm = 0.9, the single-roll pattern exists at Bi ≤ 35 while the two-roll pattern is
Fig. 15. The variation of average Nusselt numbers with the Biot number.
labeled in Fig. 11. In this case, only three stable states can be observed when the Rayleigh number is large enough. The critical Rayleigh number for convection onset at Θm = 0.7 further increases to 181000. When the Rayleigh number exceeds this value, the two-tori flow pattern appears first, which is same with the result at Θm = 0.5, as depicted in Fig. 11(a). However, the flow region becomes obviously thinner. When the Rayleigh number varies from 1.9 × 105 to 2.45 × 105, a new flow pattern with six cells near the sidewall and one cell in the center is observed. Meanwhile, each cell appears as the flow pattern with the hot-rising flow surrounded by the cold-descending flow. Therefore, we name it as seven-cell pattern, as 336
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
obtained at 35 < Bi ≤ 65. When Bi > 65, no flow happens. Fig. 15 shows the variation of the average Nusselt number with the Biot number at Ra = 103. We can see that the Biot number has a slight impact on the Nusselt number. Besides, with the increase of the Biot number, the average Nusselt number appears to increase at first and then decrease. Furthermore, it is also obvious that the average Nusselt number decreases with the increase of the density inversion parameter. When Ra = 5 × 103, the variation of the average Nusselt number with the Rayleigh number is the same with that at Ra = 103, as shown in Fig. 15(b). However, there is a sudden increase of the average Nusselt number at the point of the flow pattern evolution. As explained in Ref. [34], the main reason for this phenomenon is that the flow pattern bifurcates from a relatively stable state to a violent state, which means the flow pattern can also influence the heat transfer characteristics.
3892–3895. [8] H.A. Dijkstra, On the structure of cellular solution in Rayleigh-Bénard-Marangoni flows in small-aspect-ratio containers, J Fluid Mech 243 (1992) 73–102. [9] T.R. Qin, Z. Tuković, R.O. Gigoriev, Buoyancy-thermocapillary convection of volatile fluids under their vapors, Int J Heat Mass Transf 80 (2015) 38–49. [10] T.R. Qin, Z. Tuković, R.O. Gigoriev, Buoyancy-thermocapillary convection of volatile fluids under atmospheric conditions, Int J Heat Mass Transf 75 (2014) 284–301. [11] B. Xu, X. Ai, B.Q. Li, Bénard-Marangoni instability in an open vertical cylinder with lateral heating, ASME Int Mech Eng Congr Expo (2005) 963–968, http://dx.doi. org/10.1115/IMECE2005-81857. [12] D. Nezar, S. Rahal, Numerical simulation of convective instabilities in a liquid layer submitted to an inclined gradient of temperature, Energy Procedia 36 (2013) 380–385. [13] S. Rahal, P. Cerisier, H. Azuma, Bénard-Marangoni convection in a small circular container: influence of the Biot and Prandtl numbers on pattern dynamics and free surface deformation, Exp Fluids 43 (2007) 547–554. [14] R.E. Sakhy, K. El Omari, Y. Le Guer, S. Blancher, Rayleigh-Bénard-Marangoni convection in an open cylindrical container heated by a non-uniform flux, Int J Therm Sci 86 (2014) 198–209. [15] T. Boeck, Low-Prandtl-number Bénard-Marangoni convection in a vertical magnetic field, Theor Comput Fluid Dyn 23 (2009) 509–524. [16] W. Guo, R. Narayanan, Onset of Rayleigh-Marangoni convection in a cylindrical annulus heated from below, J Colloid Interface Sci 314 (2007) 727–732. [17] H. Mlaouah, T. Tsuji, Y. Nagano, A study of Non-Boussinesq effect on transition of thermally induced flow in a square cavity, Int J Heat Fluid Flow 18 (1997) 100–106. [18] K. Sasaguchi, K. Kuwabara, K. Kusano, H. Kitagawa, Transient cooling of water around a cylinder in a rectangular cavity - a numerical analysis of the effect of the cylinder, Int J Heat Mass Transf 41 (1998) 3149–3156. [19] A. Osorio, R. Avila, J. Cervantes, On the natural convection of water near its density inversion in an inclined square cavity, Int J Heat Mass Transf 47 (2004) 4491–4495. [20] T. Nishimura, A. Wake, E. Fukumori, Natural convection of water near the density extremum for a wide range of Rayleigh numbers, Numer Heat Transf A-Appl 27 (1995) 433–449. [21] D.S. Lin, M.W. Nansteel, Natural convection heat transfer in a square enclosure containing water near its density maximum, Int J Heat Mass Transf 30 (1987) 2319–2329. [22] X.F. Yuan, Y.R. Li, L. Peng, Y. Liu, Numerical study on natural convection of water near its density maximum in horizontal annulus, J Supercond Nov Magn 23 (2010) 1105–1109. [23] Y.R. Li, X.F. Yuan, C.M. Wu, Y.P. Hu, Natural convection of water near its density maximum between horizontal cylinder, Int J Heat Mass Transf 54 (2011) 2550–2559. [24] Y.R. Li, Y.Q. Ouyang, Y.P. Hu, Pattern formation of Rayleigh-Bénard convection of cold water near its density maximum in a vertical cylindrical container, Phys Rev E 86 (2012) 046323. [25] Y.P. Hu, Y.R. Li, C.M. Wu, Comparison investigation on natural convection of cold water near its density maximum in annular enclosures with complex configurations, Int J Heat Mass Transf 72 (2014) 572–584. [26] B. Gebhart, J.C. Mollendorf, New density relation for pure and saline water, Deep Sea Res 24 (1977) 831–848. [27] R. Touihri, A. El Gallaf, D. Henry, H. Ben Hadid, Instabilities in a cylindrical cavity heated from below with a free surface. I. Effect of Biot and Marangoni numbers, Phys Rev E 84 (2011) 056302. [28] D.J. Ma, D.J. Sun, X.Y. Yin, Multiplicity of steady states in cylindrical RayleighBénard convection, Phys Rev E 74 (2006) 037302. [29] Y.R. Li, Y.P. Hu, Y.Q. Ouyang, C.M. Wu, Flow state multiplicity in Rayleigh-Bénard convection of cold water with density maximum in a cylinder of aspect ratio 2, Int J Heat Mass Transf 86 (2015) 244–257. [30] J.C. Buell, I. Catton, The effect of wall conduction on the stability of a fluid in a right circular cylinder heated from below, J Heat Transf-Trans ASME 105 (1983) 255–260. [31] K. Borońska, L.S. Tuckerman, Extreme multiplicity in cylindrical Rayleigh-Bénard convection. I. Time dependence and oscillations, Phys Rev E 81 (2010) 036320. [32] K. Borońska, L.S. Tuckerman, Extreme multiplicity in cylindrical Rayleigh-Bénard convection. II. Bifurcation diagram and symmetry classification, Phys Rev E 81 (2010) 036321. [33] Y.P. Hu, Y.R. Li, C.M. Wu, Aspect ratio dependence of Rayleigh-Bénard convection of cold water near its density maximum in vertical cylindrical containers, Int J Heat Mass Transf 97 (2016) 932–942. [34] Y.R. Li, H. Zhang, L. Zhang, C.M. Wu, Three-dimensional numerical simulation of double-diffusive Rayleigh-Bénard convection in a cylindrical enclosure of aspect ratio 2, Int J Heat Mass Transf 98 (2016) 472–483.
4. Conclusions A series of three-dimensional numerical simulations on RayleighBénard-Marangoni convection of cold water near the maximum density in a cylindrical pool heated from below and with the heat exchange at the top free surface were carried out. From the simulation results, the following conclusions can be drawn: (1) At a small Biot number, the critical Rayleigh number for the convection onset is very small because of the influence of Marangoni effect. With the increases of the Biot number and the density inversion parameter, the critical Rayleigh number increases greatly. At a large Biot number, flow patterns become more abundant. (2) During the flow pattern bifurcation, the hysteresis phenomenon is certified. Besides, at the same Rayleigh number, there coexists a variety of the flow patterns. The Biot number also influences the flow patterns, but the effect degree is not as obvious as that of the Rayleigh number. (3) The average Nusselt number strongly depends on the Rayleigh number, but also is influenced by the Biot number and the density inversion parameter. Moreover, the effect of the flow patterns is relatively slight. Acknowledgments This work is supported by National Natural Science Foundation of China (Grant No. 51376199). References [1] B. Trouette, E. Chenier, F. Doumenc, C. Delcarte, B. Guerrier, Transient RayleighBenard-Marangoni solutal convection, Phys Fluids 24 (2012) 196–208. [2] E.M. Alawadhi, A solidification process with free convection of water in an elliptical enclosure, Energy Conv Manag 50 (2009) 360–364. [3] S.G. Egorov, I.F. Chervonyi, R.N. Volyar, The effect of preheating inductor electromagnetic field frequency on temperature gradient in the process of silicon singlecrystal growth, Inorg Mater 45 (2009) 596–598. [4] A.P. Mackwood, R.C. Crafer, Thermal modelling of laser welding and related processes: a literature review, Opt Laser Technol 37 (2005) 99–115. [5] H. Benard, Les tourbillons cellulaires dan une nappe liquid transportant dela chaleur par convection en regime permanent, Ann Chim Phys 23 (1901) 62–144. [6] M.J. Block, Surface tension as the cause of Bénard cells and surface deformation in a liquid film, Nature 178 (1956) 650–651. [7] T. Ondarcuhu, G.B. Mindlin, H.L. Mancini, C.P. Garcia, Dynamical patterns in Bénard Marangoni convection in a square container, Phys Rev Lett 70 (1993)
337
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Direct numerical simulation of Rayleigh-Bénard-Marangoni convection of cold water near its density maximum in a cylindrical pool
MARK
You-Rong Lia,∗, Lu Jiaa, Li Zhanga, Yu-Peng Hub, Jia-Jia Yua a Key Laboratory of Low-grade Energy Utilization Technologies and Systems of Ministry of Education, College of Power Engineering, Chongqing University, Chongqing 400044, China b Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang 621999, China
A R T I C L E I N F O
A B S T R A C T
Keywords: Rayleigh-Bénard-Marangoni convection Density inversion Flow pattern Heat transfer Cylindrical pool
In order to understand the influence of Marangoni effect on Rayleigh-Bénard convection, a series of threedimensional numerical simulations of Rayleigh-Bénard- Marangoni convection of cold water near its maximum density in a cylindrical pool heated from below and with a free surface at the top were carried out by using the finite volume approach. Results show that the flow pattern is strongly influenced by the Rayleigh number, as well as by the Biot number. When the Biot number is small, the critical Rayleigh number for the onset of convection is low due to the strong Marangoni effect on the free surface. As the Biot number increases, the Marangoni effect becomes weaker. Therefore, the critical Rayleigh number becomes much larger and more flow patterns can be observed at different density inversion parameters. Furthermore, the average Nusselt number is influenced by the Rayleigh number, the Biot number and the density inversion parameter. Specifically, the average Nusselt number increases with the increase of the Rayleigh number and the decrease of the density inversion parameter. Meanwhile, as the Biot number goes up, the average Nusselt number increases at first and then decreases. In addition, there exists a sudden increase of the Nusselt number when the flow pattern transits from a relative stable state to a violent state, which means the flow pattern can also influence the heat transfer ability.
1. Introduction For a cavity heated from the bottom and with a free surface at the top, when the temperature difference between the bottom and top is large enough, the flow driven by the buoyancy and the surface tension gradient will happen, which is named as Rayleigh- Bénard-Marangoni (R-B-M) convection [1]. Because of its wide existence in industrial applications, such as cold storage air conditioning, crystal growth, welding process etc. [2–4], this convection has been widely concerned, especially in recent decades. The pioneering work on Bénard convection was carried out by Bénard [5]. According to his experiment, it was found that the fluid keeps stationary at a small temperature difference; however, the hexagon flow pattern appears when the temperature difference exceeds a certain value. After that, Block [6] pointed out that the hexagonal cell pattern may be caused by the surface tension gradient. Ondarcuhu et al. [7] investigated the effect of the temperature difference on the flow pattern of Bénard-Marangoni convection in a rectangular cavity by experiments. Results indicated that when the temperature difference is small, the flow pattern appears as the pentagon and square patterns. ∗
With the gradual increase of the temperature difference, the pentagon pattern disappears at first, and then reproduces. Dijkstra [8] numerically analyzed the stability of two-dimensional steady flow patterns in rectangular containers, and obtained the bifurcation diagram of flow patterns. Similar studies of a volatile fluid in a sealed and confined rectangular container were carried out by Qin et al. [9,10]. Then, Xu et al. [11] studied R-B-M convection in a cylinder heated laterally and cooled at top by using higher order finite difference scheme. They found that the critical Grashof number increases with the increase of the Marangoni number. Nezar and Rahal [12] numerically investigated the effect of gravity and surface tension gradient on the flow patterns. It was concluded that with the increase of the thickness of liquid layer, the flow pattern would bifurcate from hexagonal cells to stationary longitudinal rolls. Rahal et al. [13] conducted a series of experiments to study the effect of the Marangoni number, Biot number and Prandtl number on the free surface deformation, and found that the flow will bifurcate to a time-dependent flow when the Biot number is large enough. Sakhy et al. [14] carried out a series of three-dimensional numerical simulations on R-B-M convection in a cylindrical container heated from below by a non-uniform heat flux, the results showed that
Corresponding author. E-mail address: [email protected] (Y.-R. Li).
http://dx.doi.org/10.1016/j.ijthermalsci.2017.10.034 Received 13 May 2017; Received in revised form 5 September 2017; Accepted 26 October 2017 1290-0729/ © 2017 Elsevier Masson SAS. All rights reserved.
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
γ γT λ Г μ ν θ Θ ρ σ τ
Nomenclature Bi d ez g h Ma Nu P Pr q r R Ra T V z Z
Biot number depth of cylindrical pool, m z-directional unit vector gravitational acceleration, m/s2 heat exchange coefficient on the free surface, W/(m2·K) Marangoni number Nusselt number dimensionless pressure Prandtl number exponent in Eq. (1) radial coordinate or radius, m dimensionless radius Rayleigh number temperature, °C dimensionless velocity vector axial coordinate dimensionless axial coordinate
Subscripts 0 ave cri e h m t
Greek symbols α
thermal expansion coefficient, (°C)-q surface tension coefficient, N/(m·K) thermal conductivity, W/(m·K) aspect ratio dynamic viscosity, kg/(m·s) kinematic viscosity, m2/s azimuthal coordinate, rad dimensionless temperature density, kg/m3 surface tension, N/m dimensionless time
initial average critical environment hot wall density inversion point top free surface
thermal diffusivity, m2/s
the flow patterns depend on the Rayleigh number, Marangoni number, Biot number and the ratio of thermal conductivities. Boeck [15] investigated the effect of homogeneous magnetic field on Bénard convection in a rectangular domain, and the transition from stationary hexagonal pattern to time-dependent convection was explored for the fluid with zero Prandtl number. Guo and Narayanan [16] performed a linear stability analysis of R-B-M convection of a silicone oil in the annular cavity. It was found that the depth of the cavity rather than the width of the gap has an important effect on the Marangoni convection. Furthermore, the critical wavelength generally decreases with the decrease of the scaled gap. For all reports mentioned above, the density of working fluid was supposed to be a linear function of temperature. However, cold water is a typical fluid which has the maximum density at the temperature around 4 °C, this special property makes the flow become more complicated [17]. Sasaguchi et al. [18] carried out numerical analyses on the influence of eccentricity on natural convection of cold water in a rectangular cavity and found that when the initial water temperature is 4 °C, the relative positon of the cylinder and rectangular cavity has a significant effect on the average Nusselt number. When the initial water temperature exceeds 4 °C,a downward flow pattern can be observed at first, and then changes to upward flow pattern gradually. Osorio et al. [19] numerically simulated natural convection of cold water in a square cavity and discussed the effect of the density inversion parameter on the flow pattern. Furthermore, Nishimura et al. [20] investigated ReB convection of cold water in a rectangular enclosure whose aspect ratio was 1.25. Results showed that the flow pattern and the temperature field are symmetric; however, a sinking cell is formed at the center of the enclosure because of the influence of the density inversion. Lin and Nansteel [21] simulated natural convection of cold water in a square enclosure heated and cooled on opposite vertical walls, and found that there is a multicellular flow pattern. Furthermore, the flow pattern is independent on the Rayleigh number at 103≤Ra≤106. Yuan and Li et al. [22,23] performed a series of numerical simulations on natural convection of cold water in the horizontal annulus and the horizontal eccentric annulus, and discussed the effect of Rayleigh number and density inversion parameter on the flow patterns. Li et al. [24] also found some new flow patterns in ReB convection of cold water at a fixed aspect ratio. Hu et al. [25] performed numerical simulations of cold water in complex configurations and obtained the heat transfer
correlations in different cavities. However, when both the buoyancy and the surface tension gradient are taken into account, the formation of the flow pattern is very complicated and different. Therefore, the aim of present work is to obtain various flow patterns of R-B-M convection of cold water in a cylindrical container. Furthermore, we also focus on the effects of Rayleigh number, Biot number and density inversion parameter on the heat transfer and evolution of flow patterns. 2. Model formulation 2.1. Basic assumption and governing equations The physical model is schematically shown in Fig. 1. We consider a cylindrical pool with the aspect ratio Г = r0/d = 1.5, which is filled with cold water near its maximum density. Where, r0 is the radius and d is the depth of the cylindrical pool. The sidewall is regarded as perfectly adiabatic. The bottom of the pool is kept at constant temperature Th, while the top is a free surface with a variation temperature. The existence of free surface means that there is the heat exchange between cold water and the environment. The heat transfer coefficient and the
Fig. 1. Physical model and the coordinate system.
328
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
environment temperature are marked as h and Te, respectively. In order to simplify, some assumptions are introduced in the present model: (1) The no-slip boundary condition on all solid-liquid boundaries is applied. (2) The cold water is incompressible Newtonian fluid and the flow is laminar. (3) Viscous dissipation is neglected. (4) The physical properties are taken as constant except for the density ρ and surface tension σ. Furthermore, the density is assumed to be a nonlinear function of temperature [26] as
ρ (T ) = ρm (1 − γ T − Tm q )
Table 1 Mesh dependence test at Bi = 5, Θm = 0.5 and Ra = 104. Mesh 20R 40R 60R 80R
4.1958 4.2309 4.2472 4.2527
1.34 0.51 0.13 –
Table 2 Comparison of average Nusselt number at different Rayleigh numbers.
where ρm = 999.972 kg/m is the maximum density at temperature Tm = 4.029325 °C, the thermal expansion coefficient γ = 9.297173 × 10−6(°C)−q, q = 1.894816. The surface tension can be expressed as
σ = σ0 − γT (T − T0)
Ra
(2)
0.5 1.0 1.5 2.0
where σ0 = σ(T0), γT = −∂σ/∂T, and T0 is the initial temperature. Using d, d2/ν, ν/d and μν/d2 as scale quantities for length, time, velocity and pressure, the dimensionless governing equations in the cylindrical coordinate system can be expressed as follows:
Nuave
× × × ×
4
10 104 104 104
Deviation (%)
Present
Ref. [27]
2.11 2.33 2.45 2.56
2.07 2.29 2.41 2.53
1.89 1.71 1.63 1.17
(3)
∇⋅V =0
∂τ V+V ⋅∇V= − ∇P +
∇2 V+ (Ra/Pr )
∂τ Θ + V ⋅∇Θ = (1/ Pr ) ∇2 Θ
Θ−Θm q eZ
2.2. Boundary conditions (4) Under above assumptions, the boundary conditions are: At the bottom wall (0 ≤ R < Г, 0≤θ < 2π, Z = 0):
(5)
where V(VR, Vθ, VZ) is the dimensionless velocity vector, Θ=(T-Tt)/(ThTt) the dimensionless temperature. Θm=(Tm-Tt)/(Th-Tt) is the density inversion parameter which defines the location of the maximum density temperature within fluid layer. ez is the unit vector in the Z direction. Meanwhile, two important dimensionless parameters, Ra and Pr, are defined respectively as: Ra = gγd3(Th−Tt)q/(αν), Pr=ν/α
Vθ = VR = VZ = 0, Θ = 1
Bi z T (z ) = (Th − Te ) + Th 1 + Bi d
Vθ = VR = VZ = 0, ∂R Θ = 0. At the free surface (0 ≤ R < Г, 0≤θ < 2π, Z = 1):
(6a-b)
(12a-d) where Ma is the Marangoni number, which is defined as:
Ma =
γT d (Th − Tt ) . μα
(13)
The average Nusselt number on the bottom wall is defined as (7)
Nuave =
1 πΓ 2
2π
∫0 ∫0
Γ
∂Θ RdRdθ. ∂Ζ
(14)
2.3. Numerical procedure and validation (8a-b) The non-uniform computational mesh is shown in Fig. 2, the finer grids cluster close to the boundaries for the reason that the temperature gradient and velocity gradient are large. Finite volume approach is applied to discretize the fundamental governing equations. The convective terms and the diffusion terms are discretized through using the QUICK scheme and the second-order central difference scheme, respectively. The SIMPLE algorithm is adopted to handle the pressure-
It is worth pointing out that the temperature difference ΔT is defined as ΔT = Th−Tt rather than ΔT = TheTe in this paper for the reason that the attention is the flow of cold water in the pool. Therefore, according to Eq. (8), the temperature difference can be expressed as
Bi (Th − Te ) 1 + Bi
(11a-b)
VZ = 0, ∂Z Θ + BiΘ + 1 = 0, ∂Z VR = −Ma∂R Θ , ∂Z Vθ = −(Ma/ R) ∂θ Θ,
where Bi is the Biot number, which is defined as Bi = hd/λ. Therefore, the temperature at the free surface is
Bi (Th − Te ) + Th 1 + Bi
(10a-b)
At the cylinder sidewall (R = Г, 0≤θ < 2π, 0 ≤ Z ≤ 1):
At a small temperature difference, as considered by Touihri et al. [27], the imposed thermal conditions lead to a conductive state without any flow. In this case, the temperature profile is linear along the vertical direction,
ΔT =
Deviation (%)
(1) 3
Tt = T (d ) = −
80θ × 10Z 120θ × 20Z 160θ × 30Z 200θ × 40Z
× × × ×
Nu
(9)
Fig. 2. The cylindrical mesh model, R-θ plane (left), R-Z plane (right).
329
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 3. Contours of the axial velocity in the Z = 1/2 plane at Bi = 5 and Θm = 0.5. Solid lines correspond to positive value; dashed lines correspond to negative value.
Fig. 4. Flow patterns obtained from different routes.
Table 3 Stability range of all stable convective branches at Bi = 5 and Θm = 0.5. Flow pattern
Corresponding figure
Range of stability
upward one-torus downward two-roll single-roll upward two-roll four-spoke dipole
Fig. Fig. Fig. Fig. Fig. Fig.
140 ≤ Ra≤650 160 ≤ Ra≤104 0.1 × 104≤Ra≤4.0 × 104 0.8 × 104≤Ra≤4.0 × 104 1.2 × 104≤Ra≤4.0 × 104 160 ≤ Ra≤500
3(a) 3(b) 3(c) 3(d) 4(a) 3(b)
Fig. 6. The variation of the average Nusselt number with Rayleigh number at Bi = 5 and Θm = 0.5.
velocity coupling, and the fully implicit scheme is used for the temporal discretization. The dimensionless time step is 10−4. At each time step, the solution is considered to be convergent when the residual errors of both the momentum and energy equations are less than 10−5. In order to verify the independence of the grid, a set of numerical simulations with four different meshes were performed at Bi = 5, Θm = 0.5 and Ra = 104. The average Nusselt numbers on the bottom wall are listed in Table 1. It can be found that the deviation of the average Nusselt numbers is only 0.13% between two denser grids. Comparing with results for grid independence in Refs. [29,34], this deviation can be completely negligible. Therefore, the grid of 60R × 160θ × 30Z is an appropriate choice considering both the
Fig. 5. Schematic diagram of stability range and the flow pattern transition at Bi = 5.
330
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 7. Flow pattern obtained from the conductive state at Θm = 0.1 and Bi = 100. Upper plot: contours of the axial velocity in Z = 1/2 plane; solid lines correspond to positive value and dashed lines correspond to negative value. Lower plot: isothermal surfaces of Θ = 0.3, 0.5 and 0.7.
moment, the hot-rising flow is located in the center of the cylindrical pool while the cold-descending flow is near the sidewall. Therefore, the flow pattern is named as upward one-torus. When the Rayleigh number is between 550 and 2340, the flow transits to an axisymmetric structure, which is named as downward two-roll. This flow pattern has a cold-descending flow in the center and two hot-rising flows that are distributed symmetrically on both sides, as shown in Fig. 3(b). When 2.35 × 103≤Ra≤2.17 × 104, a single-roll pattern with two approximate symmetrical sectors is observed, as shown in Fig. 3(c). When the Rayleigh number exceeds 2.18 × 104, the flow evolves into an upward two-roll, as shown in Fig. 3(d), which has the opposite flow direction with the flow pattern observed at 550 ≤ Ra≤2340. Li et al. [29] had determined that the critical Rayleigh number of RB convection onset for cold water at Θm = 0.5 is Racri = 1.73 × 104, which is much larger than that of R-B-M convection onset in a cylindrical pool with the free surface at the same condition. It hints that the top solid wall can remarkably enhance the stability of the fluid layer. On the other hand, the nonuniform temperature distribution on the free surface leads to the surface tension gradient. This surface tension gradient depends on the Biot number and is very large at small Biot number. Meanwhile, it can help to overcome the stabilization effect of viscosity, as considered by Block [6]. As a result, the critical Rayleigh number of R-B-M convection is smaller than that of R-B convection.
numerical accuracy and computational expenses. To validate the numerical method, some simulations were performed at the same condition with Ref. [27]. When Bi = 100 and Ma = 0, the comparison of the average Nusselt numbers is given in Table 2. Obviously, the obtained results are very close to those in the reference and the maximum deviation is less than 1.9%. Furthermore, the same flow pattern is also observed. As a consequence, the numerical method is accurate enough to simulate Rayleigh-Bénard-Marangoni convection of cold water near its maximum density in a cylindrical pool.
3. Results and discussion 3.1. Effect of Rayleigh number Firstly, we analyze the effect of the Rayleigh number on the flow pattern and the heat transfer characteristics in the case of Bi = 5. A series of numerical simulations at the Rayleigh number range from 100 to 4.0 × 104 were carried out when the conductive solution was used as the initial condition. As shown in Fig. 3, four typical flow patterns orderly appear with the increase of the Rayleigh number. When Ra < 160 at Θm = 0.5, the fluid keeps stationary and no flow appears in the cylindrical pool. When Rayleigh number increases to the critical value of Racri = 160, the fluid begins to destabilize. The conductive solution bifurcates to a steady one-torus flow pattern, which is the same with the centrosymmetric flow pattern observed by Ma et al. [28]. At this
3.1.1. Evolution from upward one-torus pattern The simulation result shows this flow pattern remains stable at 331
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 8. Evolution processes of various flow pattern.
Fig. 9. Flow pattern obtained from the conductive state at Θm = 0.5 and Bi = 100. Upper plot: contours of the axial velocity in Z = 1/2 plane; solid lines correspond to positive value and dashed lines correspond to negative value. Lower plot: isothermal surfaces of Θ = 0.3, 0.5 and 0.7.
140 ≤ Ra≤650. However, the critical Rayleigh number of the flow transition from the one-torus pattern to the conductive state is smaller than that of the flow transition from the conductive state to the onetorus pattern. Therefore, there exists the transition hysteresis of the flow pattern, which is also found by Buell and Catton [30]. When the Rayleigh number is between 700 and 2850, the upward one-torus pattern turns into a downward two-roll pattern. Furthermore, it bifurcates to a single-roll pattern at 0.3 × 104≤Ra≤2.3 × 104. Finally, when the Rayleigh number is above 2.4 × 104, the upward two-roll pattern appears.
3.1.2. Evolution from downward two-roll pattern When the downward two-roll pattern converged at Ra = 1.5 × 103 is used as an initial condition, it's found that this flow pattern is stable at 160 ≤ Ra≤ × 104. When the Rayleigh number is below 160, it decays to the conductive state. When the Rayleigh number is beyond 1.2 × 104, this pattern bifurcates to a new flow pattern that is called as four-spoke pattern, as shown in Fig. 4(a).
3.1.3. Evolution from single-roll pattern The single-roll pattern remains stable at a wide range of 332
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 10. Evolution from flabellum pattern to two-tori pattern at Ra = 8.0 × 104. Fig. 11. Flow pattern obtained from the conductive state at Θm = 0.7 and Bi = 100. Upper plot: contours of the axial velocity at Z = 1/3 plane; solid lines correspond to positive value and dashed lines correspond to negative value. Lower plot: isothermal surfaces of Θ = 0.3, 0.5 and 0.7.
rather than one-torus pattern. In summary, because of the strong surface tension gradient at Bi = 5, the flow patterns are relatively simple and mainly affected by the Rayleigh number. Specifically, four distinct flow patterns can be obtained when the conductive solution is used as the initial condition, including upward one-torus, downward two-roll, single-roll and upward two-roll. Furthermore, two new flow patterns, four-spoke and dipole patterns, are also observed when the different flow states are used as initial conditions. The stability range of the different flow patterns and schematic diagram on flow pattern transitions are given in Table 3 and Fig. 5, respectively. It can be found that there exist a variety of the flow patterns at the same Rayleigh number and the transition of the flow patterns is abundant. The variation of the average Nusselt number on the bottom with the Rayleigh number is shown in Fig. 6. It can be found that the average Nusselt number strongly depends on the Rayleigh number. Specifically, for all the flow patterns, the average Nusselt number will increase simultaneously with the Rayleigh number. Besides, at a constant Rayleigh number, the average Nusselt number of different flow patterns is also different, which means that the flow pattern can also influence heat transfer ability. With the increase of the Biot number, the free surface temperature Tt approaches gradually the environment temperature Te. Therefore, the effect of the thermocapillary force on the free surface becomes weaker. When the density inversion parameter is small, the density gradient direction in the fluid is almost same as that of the common fluid. The simulation results at Θm = 0.1 and Bi = 100 are displayed in Fig. 7. The critical Rayleigh number for the convection onset is 2750 if the conductive solution is used as the initial condition. When the Rayleigh number is over the threshold, the upward one-torus state appears first. At this time, the intensity of flow is very weak, as shown
Table 4 Stability range of all stable convective branches. Θm
Flow pattern
Corresponding figure
Range of stability
0.1
upward one-torus sadman four-spoke X two-roll downward one-torus one-torus dice
Fig. 7(a) Fig. 7(b) Fig. 7(c) Fig. 7(d) Figs. 7(e) and 8(b) Fig. 8(a)
2.75 × 103≤Ra≤0.55 × 104 0.3 × 104≤Ra≤0.85 × 104 0.6 × 104≤Ra≤6.8 × 104 4.7 × 104≤Ra≤8.0 × 104 0.55 × 104≤Ra≤8.0 × 104 2.0 × 104≤Ra≤8.0 × 104
Fig. 8(c)
0.3 × 104≤Ra≤0.58 × 104
0.5
two-tori three-spot flabellum
Fig. 9(a) and (d) Fig. 9(b) Fig. 9(c)
1.84 × 104≤Ra≤1.25 × 105 1.85 × 104≤Ra≤6.5 × 104 1.85 × 104≤Ra≤7.8 × 104
0.7
two-tori seven-cell eight-cell
Fig. 11(a) and (d) Fig. 11(b) Fig. 11(c)
1.75 × 105≤Ra≤3.5 × 105 1.83 × 105≤Ra≤3.5 × 105 1.9 × 105≤Ra≤3.5 × 105
103≤Ra≤4 × 104 when the flow state at Ra = 0.5 × 104 is used as an initial condition. When the Rayleigh number varies from 150 to 900, it transits into the downward two-roll pattern. If Ra≤150, the system decays to the conductive state. 3.1.4. Evolution from upward two-roll pattern In this case, the simulation results show that the upward two-roll pattern exists at a wide range of the Rayleigh number from 0.8 × 104 to 4 × 104. When Ra < 0.7 × 104, it bifurcates to the single-roll pattern. If the Rayleigh number is less than 500, the pattern turns into a new flow state that is named as the dipole pattern, as depicted in Fig. 4(b), 333
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 12. Schematic diagram of stability range and the flow pattern transition at Bi = 100.
flow pattern is also detected by Borońska and Tuckerman [31,32], which is a secondary bifurcation flow pattern from the star pattern. When 0.66 × 104≤Ra≤7 × 104, the flow bifurcates into the spoke
in Fig. 7(a). When the Rayleigh number is between 3200 and 6500, the flow bifurcates to a new flow pattern, which looks like a sadman. Therefore, it is called as the sadman pattern, as shown in Fig. 7(b). This 334
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
pattern. When the Rayleigh number is continually increased to 2.0 × 104≤Ra≤8.0 × 104, the upward one-torus pattern evolves into the downward one-torus pattern, and the evolution process is displayed in Fig. 8(a). 3.1.6. Evolution from sadman pattern When the sadman pattern converged at Ra = 0.5 × 104 is used as the initial condition, it is found that this flow pattern is stable at 0.3 × 104≤Ra≤0.85 × 104. When the Rayleigh number increases from 0.85 × 104 to 4.0 × 104, it turns into the two-roll pattern, as shown in Fig. 8(b). When Ra > 4.0 × 104, the aberrant two-roll pattern is observed. 3.1.7. Evolution from spoke pattern The spoke patterns include the four-spoke pattern and the X pattern. The simulation result shows that the four-spoke pattern is stable in the range of 0.6 × 104≤Ra≤6.8 × 104. When Ra > 7.0 × 104, it evolves into the X pattern. On the contrary, the spoke pattern decays to a new dice pattern at 0.3 × 104≤Ra≤0.58 × 104, which is also observed by Hu et al. [33], and this evolution process is displayed in Fig. 8(c). The X pattern is stable at 4.7 × 104≤Ra≤8.0 × 104, and transits into the four-spoke pattern and the dice pattern with the decrease of the Rayleigh number. 3.1.8. Evolution from two-roll pattern The two-roll flow pattern exists at the range of the Rayleigh number between 0.56 × 104 and 8.0 × 104 when the result at Ra = 7.5 × 104 is used as the initial condition. With the decrease of the Rayleigh number, this flow pattern becomes more regular. At Ra < 0.55 × 104, the sadman pattern appears and it evolves into the upward one-torus pattern at Ra≤0.3 × 104. Therefore, there is no new flow pattern during this evolution. When the density inversion parameter is increased to Θm = 0.5 at Bi = 100, there are only four typical flow patterns when the conductive state is used as the initial condition, as shown in Fig. 9. In this case, the critical Rayleigh number for the convection onset increases to 1.84 × 104. When the Rayleigh number is over this value, the two-tori flow pattern appears firstly, as shown in Fig. 9(a). When the Rayleigh number varies from 2.2 × 104 to 3.4 × 104, the three-spot pattern is observed, which is consist of three hot-rising flows near the sidewall and a cold-descending flow in the center, as shown in Fig. 9(b). When the Rayleigh number is between 3.5 × 104 and 4.8 × 104, a new flow pattern stably exists, which looks like a flabellum. This flow pattern contains four flow cells near the surrounding wall and a hot-rising flow in the center of the cylindrical pool, as shown in Fig. 9(c). Furthermore, in a wide range of 5.0 × 104≤Ra≤1.25 × 105, the two-tori flow pattern appears again, but the flow becomes more intense, as shown in Fig. 9(d). Finally, the flow becomes oscillatory at Ra > 1.25 × 105. When the two-tori pattern is used as the initial condition, there is no flow pattern evolution regardless of the variation of the Rayleigh number, which means this flow pattern is relatively stable.
Fig. 13. The variation of the average Nusselt number with Rayleigh number.
pattern. When the Rayleigh number is small, it appears as the regular four-spoke pattern, like that in Fig. 7(c). However, with the further increase of the Rayleigh number, the region of the hot-rising flow becomes smaller, which leads to a ‘X’ flow pattern, as depicted in Fig. 7(d). When the Rayleigh number varies from 7.2 × 104 to 8.0 × 104, an aberrant two-roll pattern is observed, as shown in Fig. 7(e). At last, the flow becomes oscillatory state at Ra > 8.0 × 104.
3.1.9. Evolution from spot pattern In this case, it's found that the three-spot pattern is stable at 1.85 × 104≤Ra≤6.5 × 104. When the Rayleigh number exceeds 6.5 × 104, it turns into the two-tori pattern. On the other hand, when the result at Ra = 4.5 × 104 is used as the initial condition, the flabellum pattern keeps stable at 1.85 × 104≤Ra≤7.8 × 104. When Ra≥8.0 × 104, it evolves into the two-tori pattern, the evolution process as shown in Fig. 10. When Ra > 1.25 × 105, the flow becomes oscillatory. When the density inversion parameter is approaching 1, the position of the maximum density in the fluid layer moves to the bottom of the pool, which means that the flow occurs only in a thin liquid layer near the bottom wall. At Θm = 0.7, the distribution of axial velocity at Z = 1/3 plane and isothermal surfaces of Θ = 0.3, 0.5 and 0.7 are
3.1.5. Evolution from upward one-torus pattern By increasing or decreasing the Rayleigh number, the evolution from the upward one-torus flow pattern is exhibited when the flow state at Ra = 3 × 103 is used as the initial condition. It is found that this flow pattern remains stable at 2.75 × 103≤Ra≤5.5 × 103. When the Rayleigh number exceeds 0.6 × 104, it bifurcates into the sadman 335
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
Fig. 14. Schematic diagram of flow patterns with the Biot number.
shown in Fig. 11(b). When 2.5 × 105≤Ra≤3.1 × 105, the seven-cell pattern evolves into the eight-cell pattern, like that in Fig. 11(c). Finally, when the Rayleigh number exceeds 3.2 × 105, the flow bifurcates into the two-tori pattern again, as shown in Fig. 11(d). When Θm = 0.7, there is no new flow pattern to be found when the primary bifurcation patterns are used as initial conditions. Since the similar secondary bifurcation results have been discussed at Θm = 0.5, it will not be repeated here. Comparing the results at Bi = 100 with those at Bi = 5, it is found that more flow patterns are observed at the large Biot number. The stability range of each flow pattern at different density inversion parameters is listed in Table 4. The schematic diagram of the flow pattern transition as a function of Rayleigh number is shown in Fig. 12, it can be found that the flow pattern becomes more abundant with the decrease of the density inversion parameter. In addition, the average Nusselt number as a function of the Rayleigh number is also depicted in Fig. 13. It is obvious that the average Nusselt number gets smaller as the density inversion parameter increases, which means the heat transfer ability becomes weaker.
3.2. Effect of Biot number When Bi = 0, the free surface can be regarded as perfectly adiabatic. In this case, there is no temperature difference between the bottom and the free surface because the sidewall is also adiabatic. Therefore, no matter how much the Rayleigh number or the density inversion parameter is, the fluid layer is always static. When Bi > 0, the heat exchange from the free surface leads to a vertical temperature gradient. If the Rayleigh number is large enough, the R-B-M convection will appear. Fig. 14 qualitatively shows the evolution of the flow pattern with the Biot number at different Rayleigh numbers. As shown in Fig. 14(a) at Ra = 103, with the increase of the Biot number, it's found that two flow patterns can be observed regardless of the density inversion parameter. In detail, when Θm = 0.1 and Bi ≤ 18, Θm = 0.5 and Bi ≤ 15 or Θm = 0.9 and Bi ≤ 11, the tworoll flow pattern first appears. However, when the Biot number is above each critical value, the two-roll pattern evolves to the one-torus pattern. Continuing to increase the Biot number, the flow pattern remains the same. When Θm = 0.1 and Bi ≥ 75, Θm = 0.5 and Bi ≥ 65 or Θm = 0.9 and Bi ≥ 60, there is no flow happening in the fluid layer. When Ra = 5 × 103, the simulation result is similar to that at Ra = 103. Specifically, at Θm = 0.1, when Bi ≤ 28, the flow state that bifurcates from the conductive state leads to a single-roll pattern. When the Biot number is between 29 and 75, the flow turns into the two-roll pattern. Finally, when Bi > 75, it appears as a sadman pattern and remains constant. At Θm = 0.5, the single-roll pattern and the two-roll pattern can be observed at Bi ≤ 42 and 42 < Bi ≤ 70, respectively. When Bi > 70, the fluid layer becomes static. When Θm = 0.9, the single-roll pattern exists at Bi ≤ 35 while the two-roll pattern is
Fig. 15. The variation of average Nusselt numbers with the Biot number.
labeled in Fig. 11. In this case, only three stable states can be observed when the Rayleigh number is large enough. The critical Rayleigh number for convection onset at Θm = 0.7 further increases to 181000. When the Rayleigh number exceeds this value, the two-tori flow pattern appears first, which is same with the result at Θm = 0.5, as depicted in Fig. 11(a). However, the flow region becomes obviously thinner. When the Rayleigh number varies from 1.9 × 105 to 2.45 × 105, a new flow pattern with six cells near the sidewall and one cell in the center is observed. Meanwhile, each cell appears as the flow pattern with the hot-rising flow surrounded by the cold-descending flow. Therefore, we name it as seven-cell pattern, as 336
International Journal of Thermal Sciences 124 (2018) 327–337
Y.-R. Li et al.
obtained at 35 < Bi ≤ 65. When Bi > 65, no flow happens. Fig. 15 shows the variation of the average Nusselt number with the Biot number at Ra = 103. We can see that the Biot number has a slight impact on the Nusselt number. Besides, with the increase of the Biot number, the average Nusselt number appears to increase at first and then decrease. Furthermore, it is also obvious that the average Nusselt number decreases with the increase of the density inversion parameter. When Ra = 5 × 103, the variation of the average Nusselt number with the Rayleigh number is the same with that at Ra = 103, as shown in Fig. 15(b). However, there is a sudden increase of the average Nusselt number at the point of the flow pattern evolution. As explained in Ref. [34], the main reason for this phenomenon is that the flow pattern bifurcates from a relatively stable state to a violent state, which means the flow pattern can also influence the heat transfer characteristics.
3892–3895. [8] H.A. Dijkstra, On the structure of cellular solution in Rayleigh-Bénard-Marangoni flows in small-aspect-ratio containers, J Fluid Mech 243 (1992) 73–102. [9] T.R. Qin, Z. Tuković, R.O. Gigoriev, Buoyancy-thermocapillary convection of volatile fluids under their vapors, Int J Heat Mass Transf 80 (2015) 38–49. [10] T.R. Qin, Z. Tuković, R.O. Gigoriev, Buoyancy-thermocapillary convection of volatile fluids under atmospheric conditions, Int J Heat Mass Transf 75 (2014) 284–301. [11] B. Xu, X. Ai, B.Q. Li, Bénard-Marangoni instability in an open vertical cylinder with lateral heating, ASME Int Mech Eng Congr Expo (2005) 963–968, http://dx.doi. org/10.1115/IMECE2005-81857. [12] D. Nezar, S. Rahal, Numerical simulation of convective instabilities in a liquid layer submitted to an inclined gradient of temperature, Energy Procedia 36 (2013) 380–385. [13] S. Rahal, P. Cerisier, H. Azuma, Bénard-Marangoni convection in a small circular container: influence of the Biot and Prandtl numbers on pattern dynamics and free surface deformation, Exp Fluids 43 (2007) 547–554. [14] R.E. Sakhy, K. El Omari, Y. Le Guer, S. Blancher, Rayleigh-Bénard-Marangoni convection in an open cylindrical container heated by a non-uniform flux, Int J Therm Sci 86 (2014) 198–209. [15] T. Boeck, Low-Prandtl-number Bénard-Marangoni convection in a vertical magnetic field, Theor Comput Fluid Dyn 23 (2009) 509–524. [16] W. Guo, R. Narayanan, Onset of Rayleigh-Marangoni convection in a cylindrical annulus heated from below, J Colloid Interface Sci 314 (2007) 727–732. [17] H. Mlaouah, T. Tsuji, Y. Nagano, A study of Non-Boussinesq effect on transition of thermally induced flow in a square cavity, Int J Heat Fluid Flow 18 (1997) 100–106. [18] K. Sasaguchi, K. Kuwabara, K. Kusano, H. Kitagawa, Transient cooling of water around a cylinder in a rectangular cavity - a numerical analysis of the effect of the cylinder, Int J Heat Mass Transf 41 (1998) 3149–3156. [19] A. Osorio, R. Avila, J. Cervantes, On the natural convection of water near its density inversion in an inclined square cavity, Int J Heat Mass Transf 47 (2004) 4491–4495. [20] T. Nishimura, A. Wake, E. Fukumori, Natural convection of water near the density extremum for a wide range of Rayleigh numbers, Numer Heat Transf A-Appl 27 (1995) 433–449. [21] D.S. Lin, M.W. Nansteel, Natural convection heat transfer in a square enclosure containing water near its density maximum, Int J Heat Mass Transf 30 (1987) 2319–2329. [22] X.F. Yuan, Y.R. Li, L. Peng, Y. Liu, Numerical study on natural convection of water near its density maximum in horizontal annulus, J Supercond Nov Magn 23 (2010) 1105–1109. [23] Y.R. Li, X.F. Yuan, C.M. Wu, Y.P. Hu, Natural convection of water near its density maximum between horizontal cylinder, Int J Heat Mass Transf 54 (2011) 2550–2559. [24] Y.R. Li, Y.Q. Ouyang, Y.P. Hu, Pattern formation of Rayleigh-Bénard convection of cold water near its density maximum in a vertical cylindrical container, Phys Rev E 86 (2012) 046323. [25] Y.P. Hu, Y.R. Li, C.M. Wu, Comparison investigation on natural convection of cold water near its density maximum in annular enclosures with complex configurations, Int J Heat Mass Transf 72 (2014) 572–584. [26] B. Gebhart, J.C. Mollendorf, New density relation for pure and saline water, Deep Sea Res 24 (1977) 831–848. [27] R. Touihri, A. El Gallaf, D. Henry, H. Ben Hadid, Instabilities in a cylindrical cavity heated from below with a free surface. I. Effect of Biot and Marangoni numbers, Phys Rev E 84 (2011) 056302. [28] D.J. Ma, D.J. Sun, X.Y. Yin, Multiplicity of steady states in cylindrical RayleighBénard convection, Phys Rev E 74 (2006) 037302. [29] Y.R. Li, Y.P. Hu, Y.Q. Ouyang, C.M. Wu, Flow state multiplicity in Rayleigh-Bénard convection of cold water with density maximum in a cylinder of aspect ratio 2, Int J Heat Mass Transf 86 (2015) 244–257. [30] J.C. Buell, I. Catton, The effect of wall conduction on the stability of a fluid in a right circular cylinder heated from below, J Heat Transf-Trans ASME 105 (1983) 255–260. [31] K. Borońska, L.S. Tuckerman, Extreme multiplicity in cylindrical Rayleigh-Bénard convection. I. Time dependence and oscillations, Phys Rev E 81 (2010) 036320. [32] K. Borońska, L.S. Tuckerman, Extreme multiplicity in cylindrical Rayleigh-Bénard convection. II. Bifurcation diagram and symmetry classification, Phys Rev E 81 (2010) 036321. [33] Y.P. Hu, Y.R. Li, C.M. Wu, Aspect ratio dependence of Rayleigh-Bénard convection of cold water near its density maximum in vertical cylindrical containers, Int J Heat Mass Transf 97 (2016) 932–942. [34] Y.R. Li, H. Zhang, L. Zhang, C.M. Wu, Three-dimensional numerical simulation of double-diffusive Rayleigh-Bénard convection in a cylindrical enclosure of aspect ratio 2, Int J Heat Mass Transf 98 (2016) 472–483.
4. Conclusions A series of three-dimensional numerical simulations on RayleighBénard-Marangoni convection of cold water near the maximum density in a cylindrical pool heated from below and with the heat exchange at the top free surface were carried out. From the simulation results, the following conclusions can be drawn: (1) At a small Biot number, the critical Rayleigh number for the convection onset is very small because of the influence of Marangoni effect. With the increases of the Biot number and the density inversion parameter, the critical Rayleigh number increases greatly. At a large Biot number, flow patterns become more abundant. (2) During the flow pattern bifurcation, the hysteresis phenomenon is certified. Besides, at the same Rayleigh number, there coexists a variety of the flow patterns. The Biot number also influences the flow patterns, but the effect degree is not as obvious as that of the Rayleigh number. (3) The average Nusselt number strongly depends on the Rayleigh number, but also is influenced by the Biot number and the density inversion parameter. Moreover, the effect of the flow patterns is relatively slight. Acknowledgments This work is supported by National Natural Science Foundation of China (Grant No. 51376199). References [1] B. Trouette, E. Chenier, F. Doumenc, C. Delcarte, B. Guerrier, Transient RayleighBenard-Marangoni solutal convection, Phys Fluids 24 (2012) 196–208. [2] E.M. Alawadhi, A solidification process with free convection of water in an elliptical enclosure, Energy Conv Manag 50 (2009) 360–364. [3] S.G. Egorov, I.F. Chervonyi, R.N. Volyar, The effect of preheating inductor electromagnetic field frequency on temperature gradient in the process of silicon singlecrystal growth, Inorg Mater 45 (2009) 596–598. [4] A.P. Mackwood, R.C. Crafer, Thermal modelling of laser welding and related processes: a literature review, Opt Laser Technol 37 (2005) 99–115. [5] H. Benard, Les tourbillons cellulaires dan une nappe liquid transportant dela chaleur par convection en regime permanent, Ann Chim Phys 23 (1901) 62–144. [6] M.J. Block, Surface tension as the cause of Bénard cells and surface deformation in a liquid film, Nature 178 (1956) 650–651. [7] T. Ondarcuhu, G.B. Mindlin, H.L. Mancini, C.P. Garcia, Dynamical patterns in Bénard Marangoni convection in a square container, Phys Rev Lett 70 (1993)
337