Flow behavior and heat transfer characteristics in Rayleigh-Bénard laminar convection with fluid-particle interaction

International Journal of Heat and Mass Transfer 146 (2020) 118840

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Flow behavior and heat transfer characteristics in Rayleigh-Bénard laminar convection with fluid-particle interaction Mufeng Chen a, Xiaodong Niu b,⇑, Peng Yu c,d,⇑, Haruhiko Yamasaki e, Hiroshi Yamaguchi f a

College of Physics and Electromechanic Engineering, Longyan University, Longyan 364012, China College of Engineering, Shantou University, Shantou 515063, China c Shenzhen Key Laboratory of Complex Aerospace Flows, Department of Mechanics and Aerospace Engineering, Southern University of Science and Technology, Shenzhen 518055, China d Center for Complex Flows and Soft Matter Research, Southern University of Science and Technology, Shenzhen 518055, China e Osaka Prefecture University, 1-1 Gakuen-cho, Naka-ku, Sakai City, Osaka 599-8531, Japan f Energy Conversion Research Center, Doshisha University, Kyoto 630-0321, Japan b

a r t i c l e

i n f o

Article history: Received 13 April 2019 Received in revised form 27 September 2019 Accepted 4 October 2019

Keywords: Rayleigh-Bénard convection Transition process Heat transfer Fluid-particle interaction Bifurcation

a b s t r a c t The present study experimentally and numerically investigates the Rayleigh-Bénard (R-B) laminar convection of a Newtonian fluid in a square cavity heated from the bottom wall at the initial conditions of constant temperature field and stationary flow field over the range of Rayleigh number Ra
4  104. The effect of a freely-moving particle in the cavity on the flow pattern and heat transfer characteristics is discussed. When the stable multi-stability patterns coexist, i.e. Ra is beyond a certain critical value, the final flow pattern can be manipulated by adjusting the initial position of the particle. The bicellular mode is established if the particle is located at the bottom center of the cavity while the unicellular mode is established if the particle is away from the center. The present results demonstrate that at the same Ra, the thermal performance is mainly determined by the flow pattern but not the presence of a particle. However, one can control the final flow pattern for R-B convection by adjusting the initial position of the particle, which consequently determines the thermal performance. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Rayleigh-Bénard (R-B) convection in a rectangular cavity heated from below is a fundamental heat transfer problem related to many industrial applications and natural phenomena such as heat dissipation of chips, the secure storage of natural liquefied gas, weather forecasting system, permafrost roadbed, and mantle convection. It should be noted that the stability of R-B convection is affected by many factors [1], including the side wall boundary conditions [2,3], the aspect ratio (A) between the horizontal dimension of the cavity to the depth [4–6], and the tilted angle when the cavity is not horizontal [7]. For example, the critical Rayleigh number (Rac) in R-B convection with constant temperature and perfectly insulating boundary conditions of side walls at A = 1 is Rac = 2585 and 5011, respectively [2]. For the perfectly conducting side wall boundary condition, Rac for different A = 2/3, 1, 1.3, 1.5, and 2 are 13380, 5030, 3510, 3090, and 2350, respectively [6]. Besides, the ⇑ Corresponding authors at: College of Engineering, Shantou University, Shantou 515063, China (X. Niu); Shenzhen Key Laboratory of Complex Aerospace Flows, Department of Mechanics and Aerospace Engineering, Southern University of Science and Technology, Shenzhen 518055, China (P. Yu). E-mail addresses: [email protected] (X. Niu), [email protected] (P. Yu). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118840 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

stability of bifurcation was greatly affected when the cavity is tilted. Mizushima and Hara [7] obtained the bifurcation map from multicellular convection to unicellular convection in a rectangular cavity (A = 1 and 2) heated from below when the cavity was tilted gradually from the horizontal direction to 90°. R-B convection has been investigated extensively. Busse [8] reviewed the non-linear behaviors of thermal convection in a layer heated from below. He summarized the dependence of heat transfer characteristics on Rayleigh number (Ra) and Prandtl number (Pr). Although he underlined the stability of convection rolls, he did not mention the multi-stability of the fluid behavior. Cliffe et al. [9] investigated R-B convection as a symmetry-breaking bifurcation problem and obtained the paths of bifurcation points. Goldhirsch et al. [10] simulated R-B convection in a twodimensional rectangular region. They observed the complicated flow patterns and the textural transition in both non-chaotic and chaotic flow regimes, as well as the multi-stability. Mizushima et al. [2,3] studied the sequential transition convection in both constant temperature and perfectly insulating boundary conditions. They found that branches of bifurcation occur, around which the flow evolves from unicellular to bicellular through a multistability process. This multi-stability indicates that the convection evolved from different initial conditions is supposed to reach one

2

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840

of the stable steady states eventually, which makes an important ramification for modeling physical thermal systems. Yamaguchi et al. [11] investigated heat transfer in magnetic fluids (regarded as Newtonian fluids), and found that the Nusselt number (Nu) increases with an increase in Ra, but experiences a sharply drop when the flow transits from unicellular mode to bicellular mode. Bodenschatz et al. [12] reviewed the convection in compressed gases and gas mixtures with Pr
~1. They also covered the convection with rotation about a vertical axis, inclination, and modulation of the vertical acceleration. Cisse et al. [13] studied the instability of R-B convection under high-frequency vibration. They concluded that mechanical vibrations may generate or suppress the convective flow induced by thermal gradient, depending on the angle between the temperature gradient and vibration. Puigjaner et al. [14] investigated bifurcations and stabilities of R-B convection, and found that different stable convective flow patterns can coexist for different ranges of Ra, where the final stable state depended on the initial conditions. Venturi et al. [15] studied the stochastic bifurcations and stabilities of natural convection within twodimensional square enclosures. They obtained the steady-state solutions and primary bifurcations as well as the coexistence of multiple stable steady states. Although significant advances in the underlying physics on R-B convection have been achieved, it is still far from complete, especially with respect to the influence of particles. More specially, heat transfer in the cavity with freely-moving particles may be enhanced, which opens a new way for the design, control and optimization of heat dissipating system [16–22]. The chief reason is that solid particles in thermal fluid highly disturb the temperature field, which consequently affect the efficiency of heat transfer. Takeuchi et al. [17] investigated the heat transfer and particle behaviors in fluidparticle thermal flow. They found that under a low heat conductivity ratio, the particles showed a simple circulating flow, while under a high heat conductivity ratio, particles created oscillation flow around the center region. Takaaki et al. [18] studied convection involving freely-moving particles, and discovered that heat conductivity ratio and particle size had significant effects on heat transfer. In this study, R-B laminar convection of a Newtonian fluid in a square cavity heated from the bottom wall is systematically investigated by experiment and simulation. This kind of convection usually occurs in a large range of phenomena in geophysics, meteorology, oceanography and engineering, e.g., core melt progression in a hypothetical severe accident in a nuclear reactor. Different from the previous studies [16–18], which investigated the heat transfer enhancement by adding many particles with different thermophyiscal properties, the present study is focusing on the effect of a single freely-moving particle on the transition modes of convection from unicellular to bicellular and the corresponding Nu. It shows that the freely-moving particle can change the final flow pattern by adjusting the particle location, which may improve heat transfer performance at moderate Ra. A coupled method, which hybrids a recently developed distribution function correction-based immersed boundary-lattice Boltzmann method (IB-LBM) for flow field with fluid-particle interaction [22] and finite difference method (FDM) for temperature field, is adopted for simulation. Corresponding experiment is conducted to evaluate the heat transfer efficiency. The thermal and flow behaviors and the transition process associated with R-B convection are explored in detail.

present study concerns that the effect of a freely-moving particle on the R-B convection transition process in a square cavity heated from the bottom. Fig. 1 shows the two-dimensional fluid-particle interaction system in the present study, where the particle boundary is described by a closed curve U which divides the fluid region X into the exterior domain X1 and the interior domain X2. According to the immersed boundary method (IBM) of Peskin [23], a fixed Eulerian mesh for the fluid and a set of Lagrangian points for the solid boundaries are adopted. The effect of the particle boundary U to the surrounding fluid can be considered by introducing a forcing term (F) and a heat source term (q) into the momentum equation and the energy equation, respectively. The working fluid is the suspension of Nano-Mn-Zn ferrite in Kerosene base fluid (TS-50 K). The previous study has shown that this fluid can be regarded as a Newtonian fluid without an applied magnetic field [24]. In this work, the Boussinesq approximation is adopted. Consequently, the governing equations for an incompressible viscous thermal Newtonian fluid with fluid-particle interaction can be expressed as follows [20–24]:

r  u ¼ 0;

ð1Þ

@ðqf 0 uÞ þ r  ðqf 0 uuÞ ¼ rp þ gf r2 u  qf 0 bðT  T 0 Þg þ F; @t

ð2Þ

@T þ u  rT ¼ kf r2 T þ q; @t

ð3Þ

where qf0, u, and T are reference density, velocity, and temperature of the fluid, respectively; p, gf, b, g, kf and T0 are pressure, dynamical viscosity, expansion coefficient under Boussinesq approximation, gravitational acceleration, thermal diffusivity and the reference temperature of the fluid, respectively; F and q are the body force density and the heat source term of the fluid, respectively, which represent the effects of the immersed particle boundary on the surrounding fluid under the frame of IBM [21–26]. The corresponding discrete values of F(rij, t) and q(rij, t) on the Eulerian mesh points can be allocated from those at the Lagrangian mesh points by introducing the Dirac delta function. The detailed procedure will be given in Sections 2.3.2 and 2.3.3 respectively. The particle dynamics will be introduced in Section 2.3.4. The dimensionless parameters governing the flow behavior and the thermal characteristic are Ra and Pr, which are given as

2. Method 2.1. Problem description of fluid-particle interaction system A particle presented in a Newtonian fluid has an influence on the flow field and temperature field in a thermal system. The

Fig. 1. Illustration of IBM for the fluid-particle interaction problems.

3

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840

Ra ¼

gf bg DTL3 ; Pr ¼ ; kf g f qf kf

ð4Þ

where L is the side length of the cavity, g is the magnitude of g, DT is the temperature difference between the hot and cold walls. In the present study, Ra is no more than 4  104, which ensure the steady flow for a pure fluid. However, it is not necessarily true when considering unpredictable perturbations, such as particles [17,18,22]. It’s reported that the particle size, volume fraction and heat conductivity affect the final flow state. 2.2. Experimental arrangement An image of the experimental setup in the present study is displayed in Fig. 2. The schematic of the experimental setup is shown in Fig. 3. The working section is a cubic container with size of 8 mm  8 mm  8 mm, which is fully filled with the working fluid. The working section is sandwiched between two acrylic conductors (8 mm  8 mm  9 mm in length, width and height). The working section and the acrylic conductors are enclosed by the acrylic plates of 1.0 mm in thickness. The whole system is then tightly wrapped by the glass-wood plywood insulation layer and the styrofoam insulation layer, which ensures perfectly thermal insulation. The bottom wall of the container is heated through an acrylic conductor by a silicon code heater and the constant high temperature is maintained by a heater controller connected to a Peltier device. The top wall is cooled by a Peltier device through another acrylic conductor and the constant low temperature is maintained by a Peltier controller. Both the bottom and top walls of the container are maintained at the constant temperatures and the top wall temperature is lower than that of the bottom wall. The experiment setup is placed in a temperature controlled room, which insures the initial temperature for the working fluid is uniform and constant. The temperatures (T1, T2, T3 and T4) along the center line of the acrylic conductors and fluid container of four different locations are measured by the thermistor temperature sensors, and the obtained data are recorded and processed by a temperature data logger. The corresponding locations of measurement points are indicated in the sketch, with L1 = L3 = L4 = L6 = 2.5 mm and L2 = L5 = 4 mm. Nu is defined as follows,

Nu ¼

hf L ; kf

ð5Þ

where kf is the thermal conductivity of the fluid, hf is the heat convection coefficient. From the energy conservation law, at the steady

state, heat convection carrying by the moving fluid between the bottom and top walls is equal to heat conduction from the bottom or top wall. Thus, hf can be calculated by the following equation,

ka

@T ¼ hf ðT t  T b Þ; @n

ð6Þ

where ka is the thermal conductivity of acrylic conductor, Tt and Tb are the top and bottom wall temperatures, respectively, and @T/on is the temperature gradient. Because the whole system is perfectly thermally insulated at all vertical surfaces, one dimensional heat transfer condition can be assumed along the vertical direction. Thus, the temperature distribution in the acrylic conductor is linear along the vertical direction, and the temperatures at the top and bottom walls can be calculated as,

T t ¼ T 2 þ L3  ðT 2  T 1 Þ=L2 ;

ð7Þ

T b ¼ T 3  L4  ðT 4  T 3 Þ=L5 :

ð8Þ

The temperature gradient can be calculated as

@T @T ¼ ðT 1  T 2 Þ=L2 or ¼ ðT 3  T 4 Þ=L5 : @n @n

ð9Þ

Table 1 lists the thermophysical properties of the working fluid (TS-50 K, Nano-Mn-Zn ferrite in Kerosene base fluid) in which qf0 denotes the fluid density at the reference temperature T0. It is worth mentioning that the working fluid shows good Newtonian behavior without imposing magnetic field [24,32,33]. The density and average diameter of the nanoparticles (Mn-Zn ferrite) are 5.00  103 kg/m3 and 9.7 nm, respectively. The mass concentration of the nanoparticles is 49.86 wt% and the number of nanoparticles per unit volume is about 5.45  1023. The volume fraction of the nanoparticles can thus be calculated, which is 14%. The fluid properties such as density, specific heat, viscosity, and thermal conductivities in Table 1 are the functions of the volume fraction of nanoparticles. However, since the particle volume fraction is fixed and the fluid properties are known, no correlations for the variations in thermal conductivity and dynamic viscosity with nanoparticle volume fraction are provided in the present study. Table 2 shows the properties of the added particle. The material of the particle is alumina. The alumina particle is hollow and its apparent density is the same as that of the working fluid. 2.3. Numerical methods 2.3.1. The numerical model The 2D simulation is performed in the present study. Under the present flow and heat transfer configuration, the 2D simulation is

T1 T2

T3 T4

Fig. 2. Experimental apparatus (left) and testing section (right).

4

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840

8

1

6

L1 L2 L3

9

4 5

8mm

2

7

L4 L5 L6

10

6

T1 T2

12

8

3 4

11

T3 T4

8mm 1

Peltier controller

2

Temperature data logger

3

Heater controller

4

Peltier device

5

Thermistor sensor

6

Heat sink

7

Styrofoam insulation

8

Acrylic conductor

9

Acrylic plate

10 Glass-wool

insulation

11 Testing

points

12 Working

fluid

Fig. 3. Sketch of experimental installation (left) and testing section (right).

U = (0, 0), θc = 0 (Cold)

Table 1 Cavity length and properties of working fluid (TS-50 K). Length of the cavity L [mm] Reference density qf0 [kg/m3] Viscosity gf [Pa s]

8

Thermal conductivity kf [W/(m K)]

0.175

1.40  103 1.68  102

Specific heat Cpf [J/(Kg K)] Reference temperature T0 [K] Gravity acceleration g [m/s2] Expansion coefficient b [1/K]

1.39  103 278

θ0 = 0.5

g

9.8 6.90  104

U

U

(0, 0) particle

y Table 2 Properties of the particle.

(0, 0)

x

Diameter D [mm]

0.5

Specific heat Cpp [J/(Kg・K)]

0.78  103

Thermal conductivity kp [W/(m K)] Density qp [kg/m3]

U = (0, 0), θh = 1 (Hot)

17 1.40  103

able to capture the main physics. Also, the 3D simulation is timeconsuming, especial for Ra close to the bifurcation point. Thus, to compromise the efficiency and accuracy, the 2D simulation is adopted. The computational domain is presented schematically in Fig. 4. In numerical simulation, the non-dimensional parameters are adopted, which are defined as,

u x y T  T0 U ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; X ¼ ; Y ¼ ; h ¼ ;t ¼ t L L DT gbDTL

rffiffiffiffiffiffiffiffiffiffiffiffi gbDT ; L

Fig. 4. Schematic of the computational domain.

Z Nu ¼ 0:5 0

1

@h dXjY¼0 þ @Y

Z 0

1

@h dXjY¼1 @Y

 ð11Þ

2.3.2. The IB-LBM for flow field In the present study, the lattice BGK model of D2Q9 is chosen as the flow solver. The evolution equation of LBM with the forcing term can be written as [22,25–30]: eq

ð10Þ

Numerical simulations are performed in a two-dimension square cavity of 1  1. The gravitational acceleration g acts in the negative y direction. The whole cavity is filled with the working fluid. The non-slip and impermeable wall boundary conditions U = (0, 0) are specified at the four side walls. The bottom and top wall temperatures are kept at hh = 1 and hc = 0, respectively. The two side walls are kept perfectly insulated, which gives oh/on = 0 there. For the cases with a particle, the fluid in the whole cavity is stationary U = (0, 0) and at a uniform constant temperature h0 = 0.5 initially (at time t = 0). However, to obtain the specified final flow state, the initial conditions for the cases without a particle may be different, which will be provided later. Once the temperature field is obtained by simulation, Nu of the whole cavity can be calculated as an average of those along the hot and cold walls,

f a ðr þ ea dt; t þ dtÞ ¼ f a ðr; tÞ 

f a ðr; tÞ  f a ðr; tÞ

sf

1 noneq þ dtð1  Þf ðr; tÞ þ dtGa ; 2sf a;b

ð12Þ

where fa (r, t) and feq a (r, t) denote the density distribution function and the equilibrium density distribution function for the discrete noneq velocity, f a;b ðr; tÞ is the non-equilibrium part of distribution function, r is the spatial vector of the Euler mesh point, dt denotes the time step, Ga is the discrete form of the buoyancy force, and sf is the relaxation time of density distribution function. eq f a ðr; tÞ can be expressed as,

  3 9 3 eq f a ðr; tÞ ¼ xa qf ðr; tÞ 1 þ 2 ðea  uÞ þ 4 ðea  uÞ2  2 u2 ; c 2c 2c

ð13Þ

where c ¼ dx=dt ¼ 1 is the unit lattice velocity, with dx the lattice pffiffiffi spacing and c ¼ cs 3, where cs is the lattice sound speed, ea and

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840

xa are the discrete lattice velocity and weight function, respectively, which depend on the lattice velocity model. The discrete external force Ga in Eq. (12) is defined as, #  " 1 3ðea  uÞ 9ðea  uÞ2 Ga ¼ xa 1  þ ea  G; 2sf c2 2c4

ð14Þ

where G is the buoyancy force with G = q0b(T  T0)g. The macro density and velocity are the collective behaviors of the distribution functions and described as follows:

8 P > < qf ¼ f a ; a P P noneq > q u ¼ f a ea þ 0:5dt½ f a;b ðr; tÞea þ G; : f a

ð15Þ

a

dependence on the relaxation time. Also, the mutual restraint between the relaxation times for density and temperature imposes a restriction on simulation of the high Ra problem. Therefore, FDM is adopted to discretize the energy equation. The continuous heat flux boundary condition is imposed at the solid-fluid interface [19–21], which reads

kp

@T @T j ¼ kf j ; @n Cþ @n C

sf

Q ðXl ; tÞ ¼ kp ð16Þ

where Ma is the Mach number. noneq The non-equilibrium part of distribution function f a;b ðr; tÞ is obtained by taking account for the effect of the non-slip boundary condition [22] due to the particle immersed in the flow field, which represents the fluid-particle interaction term. Firstly, to model the effect of the immersed boundary, the density distribution functions at the neighboring Eulerian mesh points are interpolated to the Lagrangian points through the Dirac delta function, which read

f ðXl ; tÞ ¼

X

2

f a ðrij ; tÞDij ðrij ; tÞh ;

ð17Þ

ij

where rij denotes the spatial vector of Euler mesh point with i, j the mesh index in the x and y directions, respectively, fa (Xl, t) is the density distribution functions at the Lagrangian boundary point Xl, and h is equal to the space of the mesh dx. The Dirac delta function [23–26] is written as follows:

Dij ðrij  Xl Þ ¼  dh ðaÞ ¼

1 2

h

 dh

   yij  Y l xij  X l dh h h

0:25½1 þ cosð0:5paÞ; jaj < 2; 0;

others:

ð18Þ

ð19Þ

where xij and yij are the coordinates of rij in the x and y directions, respectively; Xl and Yl are the coordinates of Xl in the x and y directions, respectively; and a is the arbitrary variable of Dirac delta function dh. Secondly, the equilibrium density distribution functions eq f a ðXl ; tÞ at the Lagrangian points are calculated from the non-slip boundary condition. The new set of non-equilibrium density distrinoneq bution functions f a;b ðXl ; tÞ on the boundary points can be written as: noneq

f a;b

eq

ðXl ; tÞ ¼ f a ðXl ; tÞ  f a ðXl ; tÞ;

ð20Þ



where a denotes the opposite direction of a. To reflect the boundary effects on the flow field, the obtained non-equilibrium distribution functions at the Lagrangian boundary points are allocated to the neighboring Eulerian mesh points, which give

X noneq  noneq F r ij ; t ¼ f a;b ðrij ; tÞ ¼ f a;b ðXl ; tÞDðrij  Xl ÞDsl ;

ð21Þ

i;j

where Dsl is the arc length of the boundary element. 2.3.3. FDM for energy equation The energy equation can be solved by LBM with double distribution functions of density and temperature. However, the convergence of LBM has a limitation for large Pr and has a great

ð22Þ

where kp and kf are the thermal conductivities of particle and fluid, respectively, C is the interface of fluid and particle, with the subscripts + and  denote the particle and the fluid sides, respectively, n is the outward normal unit vector on the solid boundary. The heat flux across the interface can be expressed as [20,21]

The relaxation parameter sf can be defined as,

pffiffiffiffiffiffiffiffi MaL 3Pr p ffiffiffiffiffiffi þ 0:5; ¼ cdt Ra

5

# #

"

" @T

1 @T

1 ¼ k f @n Cþ qCpp @n C qCpf

ð23Þ

where the square brackets [ ] denotes the jump function of an arbitrary variable u across the interface, which is

½/ðX l Þ ¼ lim ½/ðX l þ enÞ  /ðX l  enÞ; e!0þ

ð24Þ

The heat source term on interface is then allocated to that on Euler mesh point, which is expressed as

qðrij ; tÞ ¼

X

QðXl ; tÞDðrij  Xl ÞDsl ;

ð25Þ

i;j

2.3.4. Forces subjected to the particle In the fluid-particle interaction problem, besides the hydrodynamic force (F) of fluid interacting with the particle, three other forces act on the suspended particle, which are the buoyancy force, particle-wall interaction force (denoted as Fp-w), and inertia force (FIner). The total force (FTotal) [22,31] acting on the particle can be expressed as

!

X qf m g þ Fpw þ FTotal ¼ 1  FðXl ; tÞDsl þ FIner qp p l

ð26Þ

where mp is the mass of the particle. F(Xl, t) is calculated by a momentum exchange method [31]. The particle-wall interaction force Fp-w is obtained by the Lennard-Jones potential [31], which can be expressed as pw

F

  8   < 2:4ePjmax 2 rp 14  rp 8 Rp rwj ; R 6 r þ f p p;j j¼1 Rp;j Rp;j rp 2 ¼ ; : 0; Rp;j > r p þ f

ð27Þ

where rp is the radius of the particle and e = r2p, Rp is the coordinate of the particle, the subscript j is the index of the discrete point on the cavity wall, rwj is the position of the discrete point j, Rp,j = |Rp - rwj|, and f is the lattice unit. The inertia force (FIner) can be expressed as

FIner ¼ mf

Vt  Vtdt ; dt

ð28Þ

where mf is the mass of fluid based on the volume occupied by the particle, V is the velocity of the particle. 3. Results and discussions Before the discussion on R-B laminar convection for the special boundary conditions as shown in Fig. 4, a brief introduction of R-B convection for random initial conditions is carried out for understanding the transition process and the stability of such convection. It is noted that both the initial condition as well as unpredictable perturbations can lead to a notable difference in

6

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840

the final flow patterns. The bifurcation diagrams for the current configuration have been discussed in detail [14,15]. At small Ra, conduction is dominated. With increasing Ra, the flow initially experiences a transition from conduction to convection (onset of R-B convection) at Rac = 2585.02 through a pitchfork bifurcation. The flow then develops into unicellular pattern after this first bifurcation point, denoted as P1 in Fig. 5. The flow experiences the second bifurcation at Ra2 = 6742.31. A bicellular mode may be observed after this second bifurcation point, denoted as P2 in Fig. 5. Initially, this bicellular mode is unstable and the flow eventually reaches the unicellular mode under any perturbations. The initially unstable bicellular mode becomes stable at Ra3 = 11279, which is denoted as P3 in Fig. 5. When Ra > Ra3, the unicellular and bicellular modes coexist, and the final mode depends on the initial condition in the cavity. The flow experiences another bifurcation at Ra4 = 21000, which yields another unstable tricellular mode. In the present study, we focus on the unicellular and bicellular modes. The convection circulation cells can be either clockwise or anticlockwise, which result in different flow patterns. However, the heat transfer performance does not affect by the flow direction, which is not distinguished in the present study. Although the multi-stability patterns may coexist under the same flow configuration, the final flow pattern for a certain boundary condition with the same initial flow condition is deterministic. A systematic analysis and understanding of final flow pattern corresponding to special initial conditions as well as unpredictable perturbations is still lacking. Therefore, in this section, we focus our attention on the R-B laminar convection at the initial condition of constant temperature field and stationary flow field as shown in Fig. 4, which usually happens at the start-up stage of cooling devices. The transition process of finial flow pattern and the corresponding heat transfer characteristic are studied in detail. The final flow pattern has a significant effect on heat transfer performance. A freely-moving particle, which can be regarded as unpredictable perturbation, is introduced into the cavity. The effect of the fluidsolid interaction on the final flow state, and consequently the heat transfer characteristic, is discussed. 3.1. Validation and grid independence study For numerical simulation, the validation and the grid independence study are performed. The case considered is the R-B laminar

convection in a square cavity without any particles. The boundary conditions for the numerical simulation are shown in Fig. 4. The grid independence study is conducted on three different meshes of 100  100, 160  160, and 200  200. The three critical Rayleigh numbers Rac, Ra2, and Ra3 are obtained on the three meshes and listed in Table 3. For the cases of searching Rac, the initial conditions are U = (0, 0) and h0 = 0.5, i.e., stationary fluid with a uniform constant temperature field. However, the initial condition for the cases of searching Ra2 is different. First, the flow and temperature fields of Ra = 20000 are calculated with the initial conditions U = (0, 0) and h0 = 0.5, which yields a stable bicellular mode. These flow and temperature fields are then set as the initial conditions for the cases of searching Ra2. At a given Ra within the range of Ra2 < Ra < Ra3, the flow initially evolves from the initial conditions of Ra = 20000 and reaches a quasi-steady state, which is still a bicellular mode. This quasisteady bicellular mode maintains for quite a long period and eventually transits to a stable unicellular mode. Our results indicate that the Ra–Nu relationship for this quasi-steady bicellular mode agrees well with that for the unstable bicellular mode. Thus, this intermediate quasi-steady bicellular mode is the unstable bicellular mode indeed. The final transition to a stable unicellular mode is actually due to the accumulated disturbance introduced by the long time LBM operation. Both Rac and Ra2 indicate the transitions from conduction to convection. At the small departures from the onset of transition and Ra > Racri (Racri = Rac or Ra2), Nu is the linear function of square root of (Ra  Racri) [32], i.e.,

Nu 

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ra  Racri

ð29Þ

Thus, Racri can be obtained by a linear extrapolating method. First, Nus at three different Ra close to each bifurcation point are calculated, i.e. Ra = 2600, 2700, 2800 for Rac and Ra = 6800, 6900, 7000 for Ra2. A linear fit can then be established between Nu2 and Ra for each bifurcation point. Eventually, Racri of each bifurcation point can be determined by applying the fact of Nu = 1 when conduction is dominated. However, Ra3 cannot be determined by the above method. In the present study, the dichotomy method is applied to search the third bifurcation point. Initially, the stable bicellular mode is obtained by setting Ra = 2  104 with the initial conditions of U = (0, 0) and h0 = 0.5. The obtained flow and temperature fields are then used as the initial conditions for the cases of searching Ra3. The initial searching range is 1  104
Ra
1.3  104, with the lower bound for the unicellular mode and the upper bound for the stable bicellular mode. The flow field of each searching is being monitored until the final state reaches. The time period for each simulation should be sufficiently long because the flow may maintain an unstable bicellular mode for quite a long time before reaching the final unicellular mode. Results indicate that the maximum differences in Rac, Ra2, and Ra3 between meshes 160  160 and 200  200 are less than 0.5%. Thus, to compromise accuracy and computational time, mesh 160  160 is chosen to perform the final simulation. Table 3 also compares the present results with those reported in the literature [2,15]. The good agreement indicates the validation of the present numerical method. Fig. 5 shows variation of Nu with Ra for the validation case. Both experimental and simulation results are reported. The agreement between the present numerical and experimental results further validates the present numerical and experimental methods. It is worth mentioning that the experiment at each Ra is 

Fig. 5. Variation of the Nusselt number with the Rayleigh number.

repeated five times and the reported Nu are the average value (x) of those obtained from the five experiments. Nu for each experimental run (xk, with k = 1, 2, 3, 4 and 5) and the average value

7

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840 Table 3 The comparison of Rac, Ra2 and Ra3 at different mesh size. Ra

Mesh

Rac Ra2 Ra3

(x ¼ 1n

n P

Reference

100  100

160  160

200  200

2.546  103 6.653  103 1.120  104

2.581  103 6.714  103 1.127  104

2.583  103 6.726  103 1.128  104

xk , with n = 5) at different Ra are recorded in Table 4. The

k¼1

uncertainty analysis for present experimental results is performed, which is evaluated by the arithmetic mean of standard deviation of each result,

rx

sP ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 2 k¼1 ðxk  xÞ ¼ nðn  1Þ

ð30Þ

The measurement uncertainty factors rx at different Ra are summarized in Table 4. The maximum rx is 6.68% at Ra = 8.99  103, while the minimum one is 3.31% at Ra = 2.50  104. The small uncertainty factor ensures good repeatability of experimental runs, and thus the experimental data are reliable. 3.2. Stability of Rayleigh-Bénard laminar convection without any particles The stability of R-B laminar convection in a square cavity without any particles is further investigated in this section. Under the present initial condition as shown in Fig. 4, the dependence of Nu on Ra in the range of Ra
4  104 for both experiments and simulations is shown in Fig. 5. The blue solid line with the empty circle symbols and the red solid line with the solid square symbols correspond to the results for steady unicellular and bicellular modes obtained from the simulation, respectively. The ‘’ symbols denote the experimental results. We do not visualize the flow field in the experiment since the working fluid is not transparent. Thus, we can only distinguish the flow patterns by the numerical simulations. However, by the detailed comparison of Nu between the numerical and experimental results in Fig. 5, the flow patterns in the experiment can be identified. The typical flow and temperature fields in different flow regimes are shown in Fig. 6, and the corresponding data points are indicated by the arrows in Fig. 5. The results discussed here are restricted to the present initial condition. Over a

2.585  103 [2,15] 6.742  103 [2,15] 1.128  104 [15]

certain range of Ra, flow from different initial condition may evolve to different final state with different heat transfer performance. When Ra is very small, the buoyance-induced flow is weak and conduction is dominated. As expected, Nu is equal to 1 at this range of Ra (Fig. 5). When Ra is beyond a critical value, which is usually defined as the critical Rayleigh number (Rac), convection appears. Rac predicted by the present study is marked by the star symbol (the point P1) in Fig. 5. Rac obtained in the present simulation is 2581.4, which accords well with those of Mizushima et al. (Rac = 2585.0) [2] and Yamaguchi et al. (Rac = 2580.7) [32]. Natural convection starts when Ra > Rac. Initially, a full cavity scale circulation (unicellular mode) is established (refer to Fig. 6a). When Ra > Ra3, two circulations (bicellular mode) are observed (refer to Fig. 6c). Fig. 5 shows that Nu increases with Ra for either unicellular mode or bicellular mode. However, a sharp drop in Nu is observed when the flow transits from the unicellular mode to bicellular mode (marked as P3 and Ra3). Both experimental observation and numerical simulation capture this transition behavior, which verifies the reliability of the present study. Ra3 obtained by the present simulation is 1.127  104, which is consistent with that obtained by Venturi et al. [15]. It is worth mentioning that when Ra2 < Ra < Ra3, the flow may experience a bicellular-structure stage under certain specified initial condition [15]. For example, if the initial flow field is assigned as that for steady state at Ra = 4  104, the bicellular mode can be observed during the transient development (refer to Fig. 6b). This bicellular mode exists for a rather long period and the corresponding Nu keeps almost constant before the flow reaches the final unicellular state. This bicellular mode is an unstable state [15]. Any perturbations due to the numerical approximation may cause the flow to deviate to the stable unicellular mode. Thus, this flow pattern is referred to the unstable bicellular mode, which is denoted by the dash red line with the empty square symbols in Fig. 5. Interestingly, the dash red line follows well with the trend of the solid red line. Nu for the unstable bicellular mode decreases with decreasing Ra. If Ra continuously decreases until Ra < Ra2, the flow

Table 4 Nu for each experimental run and the average value at different Ra. Ra  103 Nu

x1

x2

x3

x4

x5



x

rx

6.84 7.65 8.99 10.31 11.50 12.80 13.50 16.42 17.35 20.92 22.46 25.00 27.23 30.78 34.54 36.78 39.00

1.708 1.934 1.943 2.143 2.392 1.678 1.965 2.033 2.336 2.495 2.493 2.624 2.851 2.877 3.083 3.142 3.122

1.848 1.965 2.081 2.282 2.612 1.732 1.972 2.245 2.398 2.548 2.586 2.704 2.873 2.939 3.165 3.183 3.169

1.857 2.088 2.165 2.297 2.274 1.825 2.098 2.265 2.444 2.617 2.611 2.712 2.964 3.021 3.177 3.218 3.209

1.942 2.149 2.261 2.357 2.362 1.849 2.107 2.372 2.503 2.623 2.688 2.756 2.998 3.123 3.199 3.241 3.215

2.051 2.267 2.321 2.433 2.475 2.028 2.142 2.408 2.521 2.827 2.691 2.826 3.164 3.151 3.322 3.379 3.337

1.8812 2.0806 2.1542 2.3024 2.4230 1.8224 2.0568 2.2646 2.4404 2.6220 2.6138 2.7244 2.9700 3.0222 3.1892 3.2326 3.2104

5.67% 6.10% 6.68% 4.79% 5.71% 6.00% 3.68% 6.56% 3.40% 5.64% 3.66% 3.31% 5.57% 5.23% 3.86% 4.02% 3.58%

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840

Isotherms

Streamlines

8

(a) Ra = 8×103, Stable

(b) Ra = 8×103, Unstable

(c) Ra = 2×104, Stable

Fig. 6. Streamlines and isotherms for stable unicellular, unstable bicellular, and stable bicellular flow patterns obtained by the present numerical method.

may transit from the bicellular mode to unstable conduction state [15]. This corresponds to another pitchfork bifurcation point denoted as P2 associated with the critical Rayleigh number Ra2. The present simulation indicates that Ra2 is around 6.714  103, which agrees well with the prediction of Mizushima (Ra2 = 6.742  103) [2]. Fig. 6a and b show the streamlines and isotherms at the stable unicellular mode and the unstable bicellular state at Ra = 8  103, respectively. The flow direction is indicated by the arrow. The direction of convective circulation is determined by chance according to the linear stability theory, and there is an equal probability that the circulation is clockwise or anticlockwise. The streamlines in Fig. 6a are diagonal symmetry. However, the streamlines in Fig. 6b show somewhat asymmetrical feature, which indicates that the flow is being deviated from the unstable bicellular mode due to numerical perturbations. Detailed comparison of isotherms in Fig. 6a and b indicates that the temperature gradients near the bottom and top walls for the bicellular mode is weaker than those for the unicellular mode. Thus, Nu for the unicellular mode is larger than that for the bicellular mode at Ra = 8  103. Actually, the sharp drop in Nu at P3 shown in Fig. 5 is due to the transition from the unicellular mode to the bicellular mode. When Ra > Ra3, the stable bicellular mode is obtained. Streamlines and isotherms at Ra = 2  104 are shown in Fig. 6c. The heated fluid moves upwards from the center of the bottom wall. When the fluid imparts the top wall, it splits into two branches, moving leftwards and rightwards, respectively. These two branches then turn the corns, move downwards along the sidewalls, turn the corns again, and eventually meet at the center of the bottom wall, which form two symmetric circulations. Over the range of Ra3 < Ra
4  104, the present initial condition always leads to the bicellular mode. The stability analysis for random initial conditions has demonstrated that the stable unicellular and bicellular modes coexist in this range of Ra [15]. At the same Ra, Nu for the unicellular mode is different from that for the bicellular mode, which indicates different heat transfer performance. Thus, in the followings, the effect of perturbation on the final flow state under the present initial condition is discussed. Particularly, a freely-particle is introduced into the cavity. The effect of its initial position on the final flow state is investigated in detail.

3.3. Stability of Rayleigh-Bénard laminar convection with a moving particle In this section, a small particle with a nondimensional radius of rp = 1/32 is put in the cavity. The effects of the different initial position of the particle are studied. The initial temperature for the experiment and simulation are higher than the reference temperature. Thus, the density of the fluid is slightly lower than that of the particle. Since the fluid is stationary initially, the particle sinks down to the bottom of the cavity. When Ra < Rac, the flow is conduction-dominated. Thus, the particle is stationary and its effect of on the flow and heat transfer behavior is negligible. When Ra > Rac, the particle moves along with the convection circulation. This can be confirmed by examining the particle Stokes number (St). St depends on the response times of the particle and the thermal system, which reads,

St ¼

sp ¼

sp tf

;

ð31Þ

qp D2 gf ; tf ¼ ; 18gf qf u2c

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi uc ¼ bg DTL ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Rakf gf =L;

qf

ð32Þ

ð33Þ

where sp and tf are the response times of the particle and the thermal system, respectively; uc is the characteristic velocity of natural convection. When 2581 < Ra < 4  104, the range of St falls in 4.2  103 < St < 6.5  102, which is less than 1. Over this range of St, the particle is a good follower to the thermal flow system, which starts moving with the onset of convection flow. The flow becomes unsteady with the moving particle and eventually approaches a periodic state. Nu is calculated by averaging the instantaneous Nusselt number over the five periods after the flow reaches the final periodic state. We restrict our discussion on the cases with Ra > Rac in this section. In the simulations, five different initial particle positions at the bottom are considered, which are PA(rp, rp), PB(0.25, rp), PC(0.4, rp), PD(0.45, rp), PE(0.5, rp). The flow and temperature fields,

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840

center breaks this symmetry, which strengthens heat transfer near its adjacent sidewall. The buoyance driven flow initially develops along the sidewall, which eventually results in a unicellular mode. Thus, the particle can be regarded as dynamic interference, which triggers different flow path towards the final flow state. The achievement of unicellular mode indicates that the particle located away from the bottom center can stabilize the natural convection mode by shifting a higher convection mode (bicellular) to a lower one (unicellular). Fig. 8 shows the trajectories of the moving particle at three different Ra. The trajectory follows the flow direction. All the trajectories shown in Fig. 8 are in the clockwise direction, but the trajectories may follow the anticlockwise direction also. At Ra = 7  103, the flow is the unicellular mode. Fig. 8a shows the particle circles around the whole cavity and the trajectory follows almost the same route at every circle. At high Ra = 2  104 and 4  104, either the unicellular mode or the bicellular mode occurs, depending on the initial position of the particle. The particle moves around the whole cavities in the unicellular mode but limits to either left or right half of the cavity in the bicellular mode. For the unicellular mode at Ra = 2  104 and 4  104, the radius of the particle trajectory at each circle shrinks gradually before it reaches the final state. The detailed comparison indicates that with increasing Ra, the trajectory in the unicellular mode is stretched in the anti-diagonal direction and the trajectory radius of the final state decreases. At a higher Ra, convection becomes stronger and the flow tends to stretch the vortex into an elliptical shape. The trajectory of the particle follows the streamline and gradually shifts the route to the streamline with the largest velocity. Similar phenomena can be observed in the bicellular mode from Fig. 8(b and c), where the trajectory in an elliptical shape at Ra = 2.0  104 is further stretched into a diamond shape at Ra = 4.0  104. It is worth mentioning that the trajectory stretches in the diagonal direction if the particle moves in the anticlockwise direction. Fig. 9 shows the corresponding temperature profile along the y direction at X = 0.25, 0.5, and 0.75 for Ra = 7.0  103 (unicellular mode) and 2.0  104 (bicellular mode). The present configuration results in the overall heat transfer in the positive y direction. If

(a) Unicellular

the particle trajectory, and Nu at different Ra are examined. For a low Ra (Ra < Ra3), the initial position of the particle does not affect the final flow pattern. A full cavity scale circulation (unicellular mode) is always established. At a high Ra (Ra > Ra3), however, convection circulation with a particle may reach either unicellular mode or bicellular mode eventually, which is determined by the initial position of the particle. Fig. 7a and b show the different isotherm patterns at Ra = 2  104 due to different initial positions of PB(0.25, rp) and PE(0.5, rp) respectively. The particle initially located at PB(0.25, rp) leads to the unicellular mode. As shown in Fig. 7a, the presence of the particle near the bottom corner enhances heat transfer there, resulting in a much asymmetrical structure near the left side wall at the beginning. Convection initially starts along the left side wall, carrying the particle to move along the walls, which causes the unicellular mode. The particle initially located at PE(0.5, rp) leads to the bicellular mode. The presence of the particle at the bottom center enhances heat transfer there, triggering an upward plume as shown in Fig. 7b. When the upward flow imparts the top wall, it naturally splits into two branches. Thus, the flow eventually forms the bicellular mode. The particle follows one branch (either left or right branch) and always circles within that half region of the cavity. At a high Ra > Ra3, the present simulation indicates that the particle located at PA(rp, rp), PB(0.25, rp), or PC(0.4, rp) results in the unicellular mode, while the particle located at PD(0.45, rp) or PE (0.5, rp) leads to the bicellular mode. Thus, we can speculate that the particle settled down near the bottom center region is supposed to trigger the bicellular mode, while the particle settled down away from the bottom center is expected to trigger the unicellular mode. For the cavity without a particle, the present symmetry initial condition with the stationary flow field and uniform temperature field always brings about the bicellular mode. The presence of the particle at the bottom can be regarded as a disturbance of the initial condition. The particle located at the bottom center promotes the formation of the thermal plume there. The induced flow thus remains symmetric initially, which triggers the bicellular mode. However, the particle away from the bottom

9

= 120s

= 180s

= 80s

= 120s

= 180s

(b) Bicellular

= 80s

Fig. 7. The isotherms for unicellular and bicellular modes at Ra = 2  104.

10

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840

(a) Ra = 7×103, unicellular

(b) Ra = 2×104, two modes

(c) Ra = 4×104, two modes

Fig. 8. Trajectories of the particle at different Rayleigh number.

the flow direction is parallel to the heat transfer direction, the temperature changes slightly along the y direction (i.e. in the region 0.3 < Y < 0.7). However, if the flow direction is normal to the heat transfer direction or the flow velocity is small, the temperature varies significantly along the y direction (i.e. in the region near the bottom or top wall). Fig. 10 plots the Nu dependence upon Ra with and without a freely-moving particle. The blue and red solid lines denote the numerical results for the stable unicellular and bicellular modes in the cavity without a particle, respectively. The red dash line denotes the numerical results for the unstable bicellular mode in the cavity without a particle. The green lines with the empty and solid symbols represent the numerical results for the sable unicellular and bicellular modes in the cavity with a freely-moving particle, respectively. The ‘+’ symbols denote the experimental results in the cavity with a freely-moving particle. The particle in the experiment settles down by chance. Since the working fluid is nontransparent, the initial position of the particle in the experiment is unknown. Although we cannot directly observe the detailed flow pattern, the dependence of Nu on Ra obtained from the experiment and the comparison between the numerical and experimental results can be applied to speculate the flow pattern in the experiment. In this case, we repeat the experiment at the same Ra for three times and Nu for each experimental run is calculated. The reported Nu for the experiment is the average value of those for the same flow pattern at the same Ra. The uncertain analysis is also performed by using the same method for the experiment without a

Fig. 9. Temperature profile along the y direction at X = 0.25, 0.5, and 0.75.

Fig. 10. Nu versus Ra in the cavity with and without a freely-moving particle.

particle. The measurement uncertainty factor rx is somewhat larger than that for the experiment without a particle, but still within the acceptable range. The less experimental run for each Ra and the unsteady flow feature may contribute to a larger rx . Fig. 10 indicates that the experimental results compare well with the numerical results. For either the unicellular or the bicellular mode, Nu increases with Ra for the cavities with and without a particle. For the same mode, the particle has a negligible effect on Nu at the same Ra and the largest error in Nu between the cavities with and without a particle is within 1.5%. Over the low range of Ra, Nu for the unicellular mode is higher than that of bicellular mode at the same Ra. The unicellular mode and the bicellular mode show the equal Nu at the point P4 with Ra4 = 2.944  104. When Ra > Ra4, Nu for the unicellular mode is smaller than that of bicellular mode at the same Ra. When Ra > Ra3, the stable unicellular and bicellular modes coexist. Ra3 for the cavity with a particle is higher than that without a particle. The presence of the particle may break up the symmetry of the flow configuration and thus the bicellular mode can only occur at higher Ra. However, the exact value of Ra3 with a particle is not provided in the present study because the initial position of the particle affects the final flow mode, which thus results in different value of Ra3. The above results indicate that the thermal performance of the present flow configuration depends on the final flow pattern. Fig. 10 shows that the high thermal performance can be achieved

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840

Fig. 11. Variation of Nu with time in the cavity with and without a particle.

by the unicellular mode when Ra < Ra4 but the bicellular mode when Ra > Ra4. Thus, the high thermal performance can be obtained by adjusting the initial position of the particle. For Ra < Ra3, the presence of the particle always leads to the unicellular mode, which is favorable for high Nu. To achieve high Nu over the range of Ra3 < Ra < Ra4, the particle initially located away from the center of the cavity bottom is preferred. When Ra > Ra4, high Nu can be obtained by putting the particle at the center of the cavity bottom. Fig. 11 compares the time evolution of Nu at Ra = 7  103, 2  104 and 4  104 in the cavity with and without a particle. For the cavity without a particle, Nu eventually reaches a steady state. For the cavity with a particle, Nu demonstrates a periodic state finally. At the same Ra, the time to reach the final state for the cavities with and without a particle is almost the same. The flow reaches the final state earlier with increasing Ra. The period for the fluctuation of Nu is shorter at a higher Ra, which corresponds to the time for the particle to complete a circle in the cavity. The moving particle influences the instantaneous Nu in the cavity. However, at the same Ra and flow pattern, the average Nu in the cavity with a particle is almost the same as that in the cavity without a particle. The present results indicate that the presence of the particle only affects the final flow pattern, depending on its initial location. The thermal performance is mainly determined by the flow pattern. The present study focuses on adding a small freely-moving particle (D/L = 1/16) in the R-B laminar convection flow and investigating its influence on flow pattern and heat transfer characteristic. We use the spherical particle in experiment and the circular particle in simulation because they are the most common shapes. We also test the particle with either smaller or larger diameter (D/L = 1/20 and 1/10), and the similar flow and heat transfer behaviors are observed. However, if the particle size is too small, the interference to the flow is not strong enough to trigger the bifurcation. On the contrary, a too large particle may not well follow the fluid motion. The range of appropriate particle size to trigger the bifurcation is difficult to determine because it is Ra-dependent. 4. Conclusions In this study, the R-B laminar convection of a Newtonian fluid in a square cavity heated from the bottom wall is systematically investigated by experiment and simulation over the range of

11

Ra
4  104. Specially, a freely-moving particle is introduced in the cavity and its influences on the flow pattern and heat transfer characteristics are examined in detail. A coupled approach, with immersed boundary-lattice Boltzmann method (IB-LBM) for flow field and finite difference method (FDM) for temperature field, is adopted for numerical simulation. The present numerical and experimental results agree with each other, which are also consistent with those reported in the literature. Although the multi-stability patterns may coexist for R-B convection, the final flow pattern for a certain boundary condition with the same initial flow condition is uniquely identified. Under the present boundary and initial conditions, the unicellular mode is established when Rac < Ra < Ra3 while the bicellular mode is established when Ra > Ra3. A sharp drop in Nu is observed when the flow transits from the unicellular mode to bicellular mode at Ra = Ra3. A freely-moving particle inside the cavity affects the final flow pattern, depending on its initial position. When Ra < Ra3, the unicellular mode is always achieved with the presence of the particle. When Ra > Ra3, the initial position of the particle determines the final flow pattern. The unicellular mode is established if the particle is initially located away from the center of the cavity bottom. However, the particle located at the center of the cavity bottom leads to the bicellular mode. It is worth noting that Ra3 for the presence of the particle is somewhat larger than that for the cavity without a particle. The Nu – Ra curves indicate that for the same flow pattern, Nu increases with Ra. For the same Ra and flow pattern, the presence of particle has a negligible effect on the thermal performance. When Ra < Ra4, Nu for the unicellular mode is larger than that for the bicellular mode. However, the situation is opposite when Ra > Ra4. The present results demonstrate that the particle can be applied to control the final flow pattern for R-B convection, which consequently manipulates the thermal performance. The present working fluid in the experiment is a magnetic fluid, which is treated as a simple Newtonian fluid in the numerical simulation. It is a natural extension to further consider the flow and thermal characteristics of the present flow configuration with an external magnetic field. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgements This project is supported by the National Natural Science Foundation for the Youth of China (Grant No. 11902135), the National Natural Science Foundation of China (Grant No. 11772179), the Shenzhen Fundamental Research Funds (Grant No. ZDSYS201802081843517), Scientific Research Foundation of Longyan University (Grant No. LB2018036), and JSPS KAKENHI Grant in Aid for Scientific Research (C) (Grant No. 26420126). References [1] G.Z. Gershuni, E.M. Zhukhovitskii, Convective Stability of Incompressible Fluids, first ed., Keterpress Enterprises, Moscow, 1976, pp. 117–143. [2] J. Mizushima, Onset of the thermal convection in a finite two-dimensional box, J. Phys. Soc. Japan 64 (7) (1995) 2420–2432. [3] J. Mizushima, A. Takahiro, Sequential transitions of the thermal convection in a square cavity, J. Phys. Soc. Japan 66 (1) (1997) 79–90. [4] M.S. Chana, P.G. Daniels, Onset of Rayleigh-Bénard convection in a rigid channel, J. Fluid Mech. 199 (1989) 257–279. [5] E.L. Koschmieder, Bénard Cells and Taylor Vortices, Cambridge University Press, 1993. [6] W. Velte, Stabilitätsverhalten und verzweigung stationärer Lösungen der navier-stokesschen gleichungen, Arch. Rational Mech. Anal. 16 (2) (1964) 97– 125.

12

M. Chen et al. / International Journal of Heat and Mass Transfer 146 (2020) 118840

[7] J. Mizushima, Y. Hara, Routes to unicellular convection in a tilted rectangular cavity, J. Phys. Soc. Japan 69 (8) (2000) 2371–2374. [8] F.H. Busse, Non-linear properties of thermal convection, Rep. Progr. Phys. 41 (12) (1978) 1929. [9] K.A. Cliffe, K.H. Winters, The use of symmetry in bifurcation calculations and its application to the Bénard problem, J. Comput. Phys. 67 (2) (1986) 310–326. [10] I. Goldhirsch, R.B. Pelz, S.A. Orszag, Numerical simulation of thermal convection in a two-dimensional finite box, J. Fluid Mech. 199 (1989) 1–28. [11] H. Yamaguchi, I. Kobori, Y. Uehata, Heat transfer in natural convection of magnetic fluids, J. Thermophys. Heat Transfer 13 (4) (1999) 501–507. [12] E. Bodenschatz, W. Pesch, G. Ahlers, Recent developments in Rayleigh-Bénard convection, Annu. Rev. Fluid Mech. 32 (1) (2000) 709–778. [13] I. Cissé, G. Bardan, A. Mojtabi, Rayleigh Bénard convective instability of a fluid under high-frequency vibration, Int. J. Heat Mass Transfer 47 (19–20) (2004) 4101–4112. [14] D. Puigjaner, J. Herrero, F. Giralt, C. Simó, Stability analysis of the flow in a cubical cavity heated from below, Phys. Fluids 16 (10) (2004) 3639–3655. [15] D. Venturi, X. Wan, G.E. Karniadakis, Stochastic bifurcation analysis of Rayleigh-Bénard convection, J. Fluid Mech. 650 (2010) 391–413. [16] A. Ueyama, S. Moriya, M. Nakamura, T. Kajishima, Immersed boundary method for liquid-solid two-phase flow with heat transfer, Trans. JSME Series B 77 (775) (2011) 803–814. [17] S. Takeuchi, T. Takaaki, K. Takeo, Effect of temperature gradient within a solid particle on the rotation and oscillation modes in solid-dispersed two-phase flows, Int. J. Heat Fluid Flow 43 (2013) 15–25. [18] T. Tsutsumi, S. Takeuchi, T. Kajishima, Heat transfer and particles behaviors in dispersed two-phase flow with different heat conductivities for liquid and solid, Flow, Turbul. Combust. 92 (2014) 103–119. [19] L. Li, C. Chen, R.W. Mei, F.K. James, Conjugate heat and mass transfer in the lattice Boltzmann equation method, Phys. Rev. E 89 (4) (2014) 043308. [20] Y. Hu, D.C. Li, S. Shu, X.D. Niu, Simulation of steady fluid-solid conjugate heat transfer problems via immersed boundary-lattice Boltzmann method, Comput. Math. Appl. 70 (9) (2015) 2227–2237.

[21] Y. Hu, C.D. Li, Full Eulerian lattice Boltzmann model for conjugate heat transfer, Phys. Rev. E 92 (6) (2015) 063305. [22] M.F. Chen, X.D. Niu, H. Yamaguchi, C. Shu, A lattice Boltzmann modeling fluidstructure interaction problems and its applications in natural convections in a square cavity with particles suspended inside, Adv. Appl. Math. Mech. 10 (2) (2018) 275–300. [23] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (3) (1977) 220. [24] H. Yamaguchi, Engineering Fluid Mechanics, Springer, Netherlands, 2008, pp. 497–537. [25] C. Shu, W.W. Ren, L.M. Yang, Novel immersed boundary methods for thermal flow problems, Int. J. Numer. Method Heat Fluid Flow 23 (2013) 124–142. [26] Y. Peng, C. Shu, Y.T. Chew, Simplified thermal lattice Boltzmann model for incompressible thermal flows, Phys. Rev. E 68 (2003) 02671–2678. [27] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1) (1998) 329–364. [28] E. Moeendarbary, T.Y. Ng, M. Zangeneh, Dissipative particle dynamics: introduction, methodology and complex fluid applications – a review, Int. J. Appl. Mech. 1 (04) (2009) 737–763. [29] A. Xu, W. Shyy, T. Zhao, Lattice Boltzmann modeling of transport phenomena in fuel cells and flow batteries, Acta Mech. Sin. 33 (3) (2017) 555–574. [30] Z.H. Chai, T.S. Zhao, A pseudopotential-based multiple-relaxation-time lattice Boltzmann model for multicomponent/multiphase flows, Acta Mech. Sin. 28 (4) (2012) 983–992. [31] X.D. Niu, C. Shu, Y.T. Chew, Y. Peng, A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows, Phys. Lett. A 354 (3) (2006) 173–182. [32] H. Yamaguchi, I. Kobori, Y. Uehata, Natural convection of magnetic fluid in a rectangular box, J. Magnet. Magnet. Mater. 201 (1999) 264–267. [33] Y. Xuan, Q. Li, Heat transfer enhancement of nanofluids, Int. J. Heat Fluid Flow 21 (1) (2000) 58–64.