Oscillatory flows of electro-thermo-convection in eccentric annulus

International Journal of Heat and Mass Transfer 134 (2019) 920–932

Contents lists available at ScienceDirect

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

Oscillatory flows of electro-thermo-convection in eccentric annulus Tian-Fu Li, Xue-Rao He, Kang Luo, Hong-Liang Yi ⇑ School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001, PR China Key Laboratory of Aerospace Thermophysics, Ministry of Industry and Information Technology, Harbin 150001, PR China

a r t i c l e

i n f o

Article history: Received 11 September 2018 Received in revised form 2 December 2018 Accepted 18 January 2019

Keywords: Electrohydrodynamics Electro-thermo-convection Oscillatory flow Lattice Boltzmann method Eccentric annulus

a b s t r a c t A detailed numerical investigation for oscillatory electro-thermo-convection of dielectric liquids in horizontal eccentric annulus is carried out based on a unified lattice Boltzmann method (LBM). The liquid is under the simultaneous action of a strong unipolar charge injection and a thermal gradient. The selfoscillatory behaviors and heat transfer characteristics for different Rayleigh number Ra, electric Rayleigh number T, Prandtl number Pr, mobility parameter M and radii ratio are studied systematically. The results show that the values of Ra within the range from 103 to 104 have little effect on the periodic oscillatory convection, while the flow develops from steady to periodic oscillatory and finally evolves into chaotic as T increases, and the fundamental frequency and the fluctuation amplitude as well as the heat transfer rate increase with T. Moreover, it is found that the fluid properties have significant impact on the oscillatory convection, and the formation of oscillating flows is strongly related to the value of radius ratio. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction As one of the most important branches of Electrohydrodynamics (EHD), Electro-thermo-hydrodynamic (ETHD) is concerned with the complex interaction between an electric field, a thermal field and a flow field. Such configurations have been extensively studied and applied in industrial systems, such as convection systems [1], boiling and condensing systems [2], drying and evaporating systems [3], solar energy systems [4] and so on. One classic problem in ETHD, named electro-thermo-convection, deals with the flow induced by the joint effect of a unipolar charge injection and a thermal gradient in dielectric liquids. This problem has received considerable attention due to some attractive advantages, like simple design, no moving parts, rapid and smart control, low power consumption [5], etc. Excellent reviews covering various aspects of electro-thermo- convection can be found in [6]. Owing to the complex mathematical model and the strong nonlinear couplings within electro-thermo-convection, most of previous researches in this problem are based on experimental study [7,8] and simple theoretical analysis [9,10]. Therefore, it is necessary to develop some advanced research methods to better understand the physical mechanism of electro-thermo-convection. During the past few decades, the direct numerical simulation has become a popular way to gain some fundamental insights into ⇑ Corresponding author. E-mail address: [email protected] (H.-L. Yi). https://doi.org/10.1016/j.ijheatmasstransfer.2019.01.090 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

many poorly understood electro-convective or electro-thermoconvective phenomena, such as the finite difference method (FDM) [11], the finite element method (FEM) [12], and the finite volume method (FVM) [13–15]. It is notable that most previous numerical studies with this topic are simulations of the simplest parallel-plane electrode configuration. However, compared with the parallel geometry, the concentric configuration is not subjected to the lateral-wall effect, and experimental set-up is more realizable. Moreover, the flow and heat transfer in configurations of complex geometries (e.g. the annular passage) is found in a variety of engineering applications, such as concentrating solar collector receiver design [16]. Accordingly, in this work we adopted the configuration of an eccentric annulus, such configuration (or similar) has only been used in a few recent numerical works [17–20]. Remarkably, all these results above are limited to twodimensional cases. As the two-dimensional simulation is a first approach to the very complex three-dimensional case, it could also provide valuable insight into the physics of problems. Actually, in [21], the complete three-dimensional linear stability problem in the case of unipolar injection between two infinitely long coaxial cylinders of arbitrary radii was considered. The two-dimensional problem with charge diffusion was analyzed by the linear stability analysis [22] and the direct numerical simulation [19]. Besides, the two-dimensional problem with the radius ratio C = 0.5 was carefully investigated by the FVM, and a complete bifurcation diagram in the finite amplitude regime was obtained [20]. Excellent agreements were obtained for the theoretical and numerical predictions

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

of the linear stability criterion. In this work, we also compare the present two-dimensional results with the analytical solutions given in [21], the results fit well, and the details would be shown in Section 3.3. In our previous report, we developed a unified LBM model for EHD and ETHD problems [23,24], and this model has been successfully extended to simulate the electro-thermo-convection in complex geometries [25,26]. In [26], we found that the numerical calculations of electro-thermo-convection gave oscillatory solutions at specific parameter combinations, and such an oscillatory phenomenon has also been observed in the plate-plate configuration for both EHD [27] and ETHD [28] flows. In fact, the similar oscillatory behaviors exist widely in the natural convection [29,30] and the multi-field coupling mixed convection flows [31– 33], it is known that these flows can exhibit interesting dynamical phenomena, and understanding the processes of heat and mass transport in oscillatory flow is a topic of interest for many applications. However, to our knowledge, the research work on oscillatory flows of electro-thermo- convection in complex geometries has never been reported. Since in [26] we only showed the oscillation at a few discrete cases with related parameters, and considering that the oscillatory behavior in electro-thermo-convection is of great significance due to its rich dynamical features and potential applications in engineering, therefore, in the present study, a detailed numerical simulation for oscillatory flows of electro-thermo-convection in eccentric annulus is investigated. The objective of this study is to investigate the effect of the related parameters on the oscillatory behavior in electro-thermo-convection. And the results show that the flow undergoes steady state, periodic oscillation and chaotic state as the electric Rayleigh number and the Rayleigh number increase, and the formation of oscillating flows is strongly related to the value of radius ratio and the fluid properties. Our finding can provide some guidance and reference for the ETHD applications in engineering practice. The remainder of this paper is organized as follows. In Section 2, we state the physical problem to be studied, and the governing equations and boundary conditions are also described. Section 3 is devoted to a brief presentation of the numerical method. Numerical results are provided and discussed in Section 4. Finally, we present our conclusion in Section 5.

2. Problem statement and governing equations Fig. 1 shows the schematic diagram for electro-thermoconvection in a two-dimensional eccentric annulus due to charge injection and thermal buoyant force. The system consists of an outer cylinder with radius Ro, within which an inner cylinder with radius Ri is located on the vertical centerline. The eccentricity C is defined as C = e/(Ro
Ri), where e is the eccentricity along vertical axis (positive upwards). The inner and outer cylinders are kept at different constant temperature h1 and h0 (h1 > h0), respectively; and a radial DC electric field is established by applying a voltage D/0 to this two electrodes. In this work, we assume that the fluid is incompressible, Newtonian and perfectly insulating, and the free charges only generate at the inner electrode, then injected into the bulk of liquid, i.e. unipolar injection. It is also assumed that the injection is autonomous and homogenous, which means the injected charge density q0 remains constant and uniform. Flow motion will be induced under the combined effect of the electric field and the temperature gradient, and finally formulate the electro-thermo-convection. The governing equations describing this problem include the Navier-Stokes equations, the energy equation, the equations for the electric field and the charge conservation equation. To further

921

Fig. 1. Schematic diagram for electro-thermo-convection in a two-dimensional eccentric annulus.

simplify the discussion, the viscous dissipation is assumed to be negligible. Moreover, the magnetic effects and Joule heating can be safely neglected since the electric current through the dielectric liquid is rather small. Following the previous assumptions and under the Boussinesq approximation, the above equations can be written as [34,35]:

ru¼0

ð1aÞ

@ ðquÞ þ r  ðquuÞ ¼
rp þ r  ðlruÞ þ f b @t

ð1bÞ

r2 / ¼
q=e

ð1cÞ

E ¼
r/

ð1dÞ

@q þ r  ðqKE
Drq þ quÞ ¼ 0 @t

ð1eÞ

@h þ r  ðuhÞ ¼ r  ðvrhÞ @t

ð1fÞ



 where the vector u ¼ ux ; uy and E ¼ Ex ; Ey denote the fluid velocity and the electric field, respectively. The scalars t, q, p, q, / and h, in turn, represent the time, fluid density, pressure, charge density, electric potential and temperature. The symbols l, e, K, D and v stand for the dynamic viscosity, electrical permittivity, ionic mobility, charge diffusivity and thermal-diffusion coefficient. The last term fb in (1b) is the total force including the electrical force f e and the buoyant force f t , and the electrical force can be further divided into the Coulomb force, dielectric force and electrostrictive force, that is:

fb ¼ ft þ fe

      1 1 @e ¼ qgb h
href þ qE
E2 re þ r qE2 2 2 @q h

ð2Þ

where b denotes the thermal expansion coefficient, g is the gravity acceleration and href is the reference temperature. In this study, the electrical permittivity e and ionic mobility K are assumed to be constant and independent of temperature and electric field. As a result, the dielectric force
12 E2 re vanishes, and the electrostrictive term h i 1 r qE2 ð@ e=@ qÞh can be viewed as a gradient of scalar, so it is 2

922

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

conventionally incorporated into the pressure term. Then the final   form for fb reads as f b ¼ qgb h
href þ qE. The above equations can be nondimensionalized using the characteristic scales H = Ro
Ri, D/0 , eD/0 =H2 , l/qH and Dh0 = h1
h0 for length, electric potential, charge density, fluid velocity and temperature, respectively. The detailed dimensionless governing equations can be expressed as:

Besides, to characterize the progress of flow and the heat transfer, the maximum fluid velocity and the Nusselt numbers will be output:

V max ¼ max Nui ¼ Ri ln

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2x þ u2y

ð5Þ

      Ro @h Ro @h ; Nu ¼ R ln o o Ri @r r¼Ri Ri @r r¼Ro

ð6Þ

ru¼0

ð3aÞ

_ @u T2 Ra þ r  ðuuÞ ¼
r p þr2 u þ 2 qC  E þ h  ey @t Pr M

ð3bÞ

in which Nui and Nuo denote local Nusselt numbers at the surfaces of inner and outer cylinders. Then the mean Nusselt number on each cylinder is computed by integration:

r2 / ¼
qC

ð3cÞ

Nui ¼

E ¼
r/ @q þr @t

ð3dÞ



 T E þ u q ¼ ar2 q M2

ð3eÞ

@h 1 þ u  rh ¼ r2 h @t Pr

ð3fÞ

_

In Eq. (3b), p is the modified pressure including an extra component from the electrostrictive force [36] and ey ¼ ½0; 1 stand for a unit vector in the y direction. Then the system is governed by six nondimensional parameters defined as:

Ra ¼

gbDhH3

mv

 1=2

;T ¼

eD/0 m q H2 1 e ; Pr ¼ ; C ¼ 0 ; M ¼ K q lK v eD/0

;a ¼

1 2p

Z 2p 0

Nui ðaÞda;

Nuo ¼

1 2p

Z 2p 0

Nuo ðaÞda

ð7Þ

where a is the angle, see Fig. 1. Because of the energy conservation law, the values of Nui and Nu0 should be equal to each other. Slight difference may be induced by the numerical method, so the final mean Nusselt number is calculated as the average value of Nui   and Nu0 , i.e. Nuav e ¼ Nui þ Nuo =2. The non-slip boundary conditions for fluid velocity are applied on the walls. Besides, the associated non-dimensional boundary conditions are summarized as:

D K D/0 ð4Þ

in which the Rayleigh number Ra represent the ratio of the buoyant force to the viscous force and the electric Rayleigh number T is defined as the ratio of the Coulomb force to the viscous force; these two are the main parameters governing the system. Both of the Prandtl number Pr and the non-dimensional mobility parameter M are physical parameters determined by the fluid properties, for typical dielectric liquids, their values generally obey Pr  5 [37] and M  3 [36]. The parameter C represents the measurement of the injection strength, for strong injection the C = 10 is frequently used, which can be viewed as a good approximation of the space charge limited (SCL) injection [38]. a is the non-dimensional charge-diffusion number with its typical value in the range from 10
4 to 10
3 [39], in this work, a constant a = 10
3 is selected.

6

2

3

5

1

0

7

4 (a)

Fig. 3. Comparison of radial dimensionless temperature distribution between our numerical results and experimental results of Ref. [47] for Pr = 0.706, C =
0.623, R* = 0.385 and Ra = 4.93  104.

8 (b)

Fig. 2. (a) Spatial discretization of D2Q9 lattice; (b) illustration of the NEES for the treatment of a curved boundary.

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

at the inner cylinder, ux ¼ uy ¼ 0, / ¼ 1, q ¼ 1 and h ¼ 1; at the outer cylinder, ux ¼ uy ¼ 0, / ¼ 0, @q=@r ¼ 0 and h ¼ 0, And the zero-field for all independent variables serves as the initial conditions. 3. Numerical solution In this section, a unified lattice Boltzmann method (LBM) developed in our previous work [23,24] is constructed to compute the

923

flow field, electric field, charge density distribution and the thermal field of the electro-thermo-convection. As a novel and powerful numerical approach for convection-diffusion problems, the LBM possesses the advantages of simple calculation procedure, efficient implementation for parallel computation, easy treatment of complex geometries [40,41], etc. Instead of solving the macroscopic governing equations directly just like most of the conventional fluid dynamic solvers, the basic idea of LBM is to solve the lattice Boltzmann equations (LBEs) that can recover the macroscopic equations in some limit. Corresponding to the macroscopic govern-

Fig. 4. Comparison between numerical and analytical solutions of the hydrostatic regime with three radius ratios R* = 0.1, 0.3 and 0.5, C = 10: (a) Charge density; (b) electric field strength.

Fig. 5. Oscillatory behaviors for Ra = 104, T = 400, Pr = 10, M = 10: (a) Time-evolution of Vmax; (b) time-evolution of Nuave; (c) phase spaces trajectory of velocity at the sampling point A; (d) Fourier frequency spectrum of Vmax.

924

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

ing equations (Eqs. (3a)–(3f)), four consistent LBEs are formulized and will be presented later, a Chapman-Enskog expansion reproducing the macroscopic equations from these LBEs had been performed in our previous work, further detailed description of the present LBM model can be found in [23] and [24]. 3.1. The lattice Boltzmann model Compared to the conventional partial differential equations (PDEs) based methods, the LBM is another powerful numerical tool for both simple and complex flows, which can naturally track the transient evolution of fluid field without the assumption of continuity medium [42]. The LBM is most frequently used to compute flow field, however, its mesoscopic idea of modelling can be extended and applied to simulate other macroscopic physical fields. Taking the body force term into consideration, the generalized evolution equation for the velocity field can be expressed as:

f ðx þ cDt; t þ DtÞ
f ðx; tÞ ¼ X þ Dt  F

ð8Þ

in which f(x,t) is the distribution function for the fluid particles with microscopic velocity c at position x and time t, Dt is the time step, X is the collision term and F is the body force. We adopt the commonly used D2Q9 velocity discretization scheme [43] together with the single relaxation time Bhatnagar-Gross-Krook (LBGK) model [44] to facilitate the implementation of our model, as shown in Fig. 2(a), the nine velocity vectors are given by:

8 j¼0 > < ð0; 0Þ j ¼ 1—4 cj ¼ cðcos½ð j
1Þp=2; sin½ð j
1Þp=2Þ > : pffiffiffi 2cðcos½ð2j
1Þp=4; sin½ð2j
1Þp=4Þ j ¼ 5—8

ð9Þ

where the streaming speed c is defined as c = Dx/Dt, Dx is the lattice size. The weight function xj with the jth velocity direction is given by:

8 j¼0 > < 4=9 xj ¼ 1=9 j ¼ 1—4 > : 1=36 j ¼ 5—8

ð10Þ

Fig. 6. Evolution of charge density distribution, temperature field and stream function contours (from top to bottom) at the wave crest and the wave trough of Nuave in Fig. 5 (b): (a) Fields at the wave trough; (b) comparison between the wave crest (solid line) and the wave trough (dash-dotted). Ra = 104, T = 400, Pr = 10, M = 10.

925

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

Then the collision operator can be expressed as:

h i eq X ¼
f j ðx; t Þ
f j ðx; t Þ =sm

where sm is the non dimensional relaxation time and librium distribution function which can be given as:

sm ¼

3m 1 þ c 2 Dt 2

f j ¼ qxj eq

ð11Þ eq fj

is the equi-

i   1h eq f j x þ cj Dt; t þ Dt
f j ðx; tÞ ¼
f j ðx; t Þ
f j ðx; tÞ þ Dt  F j

sm

ð14Þ ð12Þ

!  2 cj  u cj  u u2 1þ 2 þ
2 cs 2c4s 2cs

pffiffiffi For D2Q9 model, sound velocity cs ¼ c= 3, then the LBE for the velocity field can be written as:

ð13Þ

In this study, the split-forcing scheme [45] is used to calculate the body forces, for each discretization direction it can be:

   
1 cj qE þ qgb h
href F j ¼ xj 1
c2s 2sm

ð15Þ

Fig. 7. Partial enlarged detail of the upper right part of temperature fields, Ra = 104, T = 400, Pr = 10, M = 10.

Fig. 8. Oscillatory behaviors for different values of Ra, T = 400, Pr = 10, M = 10: Time-evolutions of Nuave for (a) 103  Ra  104, (b) Ra = 5  104 and (c) Ra = 105; (b) the phase spaces trajectory of velocity at the sampling point A, Ra from 103 to 104.

926

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

The mass density of fluid and macroscopic velocity are given by:

X q ¼ f j; j

X Dt qu ¼ cj f j þ ðqE þ qgÞ 2 j

ð16Þ

in which gj, hj and lj are the distribution functions for the electric potential, charge density distribution and temperature, respectively. The corresponding equilibrium distribution functions g eq j , eq

eq

hj , and lj are defined by:

Similarly, the LBEs for electric potential (17a), charge density (17b) and temperature field (17c) can be written as following:

i   1h g j x þ cj Dt; t þ Dt
g j ðx; tÞ ¼
g j ðx; tÞ
g eq j ðx; t Þ þ Dt  Sj

s/

g eq j ðx; t Þ ¼ xj / ( eq hj ðx; t Þ

¼ qxj

ð17aÞ 



hj x þ cj Dt; t þ Dt
hj ðx; tÞ ¼
  lj x þ cj Dt; t þ Dt
lj ðx; t Þ ¼


1h

sq

1h

sh

hj ðx; tÞ


i

eq hj ðx; t Þ

eq

lj ðx; t Þ
lj ðx; tÞ

i

ð17bÞ

ð18aÞ )
2 cj ðKE þ uÞ
c2s ðKE þ uÞ2 cj ðKE þ uÞ 1þ þ c2s 2c4s ð18bÞ

" lj ¼ hxj 1 þ eq

#  2 cj  u cj  u u2 þ
2c4s 2c2s c2s

ð18cÞ

ð17cÞ

Fig. 9. Oscillatory behaviors for different values of T, Ra = 104, Pr = 10, M = 10: Time-evolutions of Nuave for (a) 100  T  400 and (b) T = 500; the phase spaces trajectory of velocity at the sampling point A for (c) T = 100, (d) T = 200, (e) T = 300 and (f) T = 500.

927

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

The source term in (17a) are formulated as Sj ¼ xj cq=e, and s/ , sq and sh are the non-dimensional relaxation time for each LBEs, which can be calculated by:

s/ ¼

3c 1 3D 1 3v 1 þ ; sq ¼ 2 þ ; sh ¼ 2 þ c Dt 2 c 2 Dt 2 c Dt 2

ð19Þ

Finally, the macroscopic electric potential, charge density distribution and temperature are given by:

X /¼ gj ; E ¼ j

X X 1 X cj g j ; q ¼ hj ; h ¼ lj : 2 s / c s Dt j j j

ð20Þ

3.2. Boundary condition treatment In the application of LBM, some distribution functions at boundary nodes are unknown, which need to be derived using the LB techniques. In this study, the 2nd order non-equilibrium extrapolation scheme (NEES) proposed by Guo et al. [46] is adopted to treat the mesoscopic boundaries for all LBEs. The advantage of the NEES for the coupled multi-field problems lies in that the same functional module can be applied to all independent variables. Fig. 2 (b) shows the illustration of the NEES for the treatment of a curved boundary, in which the curved wall separates the whole region into three parts: the flow region with open circle, the solid wall with light grey circle and the curved boundary with the black circle. In the NEES, the unknown distribution function at xb can be divided into an equilibrium component and a non-equilibrium part as following:

Rj ðxb ; t Þ ¼ Req ðxb ; t Þ þ Rneq ðxb ; t Þ j j

ð21Þ

where Req ðxb ; t Þ is the equilibrium component, which can be comj ðx b ; t Þ puted using the macroscopic boundary conditions, and Rneq j is the non-equilibrium part, which can be obtained by a 2nd interpolation approximation: Rneq ðxb ; t Þ ¼ j

8     eq if D P 0:75 < Rj xf ; t
Rj xf ; t h  h    i   i : D R xf ; t
Req xf ; t þ ð1
DÞ R xff ; t
Req xff ; t if D < 0:75 j j j j

ð22Þ where xf and xff are neighboring fluid nodes, see Fig. 2(b). For more details description of the NEES, we refer the interested readers to Ref. [46].

means the fluid keeps rest while charge transfer solely by drift due to electric field, the flow motion arises only when the Coulomb force is sufficiently strong to overcome the damping action of the viscous force. There are analytical solutions of the hydrostatic state available for the concentric cylinder (C = 0), which may be expressed as [21]:

 
1=2 qðr Þ ¼ Ae r2 þ Be ;

Er ðr Þ ¼

1=2 Ae  2 r þ Be r

ð23Þ

where Ae and Be are two constants depending on the ratio between the radii of cylinders R⁄ and the injection strength C. For C = 10 and various value of R⁄, the values of Ae and Be can be found in [21]. In Fig. 4(a) and (b) we have displayed the charge density and the electric field strength in the radial direction for C = 10, R⁄ = 0.1, 0.3 and 0.5, respectively. A very good agreement between our numerical solutions and the analytical ones is obtained. More comprehensive validation of present model can be found in our upfront paper [25] and [48]. A uniform Cartesian grid is used to discretize the whole computational domain, and the grid independence of the solution scheme has been checked with additional simulations on a series of grids (with cells from 100  100 to 500  500). Considering simulated accuracy and computational cost, the grid of 400  400 is sufficiently fine to ensure the grid independent solution and therefore will be employed for all subsequent simulations. 4. Results and discussion After code validation and the grid independence test, a series of carefully selected computations are then performed to investigate the oscillatory flows of electro-thermo-convection, and the results will be presented and discussed in this section. In all numerical simulations, C = 10 and a = 10
3 are fixed, the eccentricity C is kept to be 0.4 and the radii ratio R⁄ is kept to be 0.3 unless specifically stated. The sampling point A with coordinate (
0.1Ro, 0.4Ro) which located at the top area between the inner and outer cylinders (see Fig. 1) is randomly chosen to track the unsteady flow characteristics of this problem. First, we present the typical characteristics of the oscillatory behavior of electro-thermo-convection in eccentric annulus. Then, by changing the main control parameters Ra and T, we investigate the effects of electric field and temperature gradient and the transition from periodic oscillatory to chaotic flows. Furthermore, we analyze the effects of fluid properties on

3.3. Model validation For the purpose of comparison, we conduct simulations to validate the present LBM model for electro-thermo-convection from two aspects, namely, pure thermal convection and pure electro convection. First of all, we conduct simulations for natural convection in eccentric annulus to validate the thermal module of our numerical method, the benchmark case with Pr = 0.706, C =
0.623, Ra = 4.93  104, and the radius ratio R⁄ = 0.385 (R⁄ = Ri/Ro) is considered. The experimental results of Kuehn and Goldstein [47] with the same configuration and parameters are used for comparison. The radial dimensionless temperature distributions of natural convection for three different u are plotted in Fig. 3, where u is the angle measured from top to line joining centers of inner and outer cylinders with center of outer cylinder as origin. It can be seen that our numerical results are in a fairly good agreement with the experimental data. Moreover, the hydrostatic solution is selected to validate the Poisson’ equation for electric potential and the charge conservation, and their coupling to the Navier-Stokes equations, i.e., the electrostatic module. And the analytical solutions of the hydrostatic state are used for comparison. The hydrostatic solution

Fig. 10. Map of flow patterns for electro-thermo-convection in a two-dimensional eccentric annulus: square, steady flow; circle, oscillatory flow; triangle, chaotic flow. The range is 5  102  Ra  5  105, 100  T  700, other parameters: Pr = 10, M = 10, R* = 0.3 and C= 0.4.

928

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

the formation of oscillating flows. Finally, the influences of the radius ratio on flow motion and heat transfer are studied. 4.1. Oscillatory flow characteristics Fig. 5(a) and (b) show the temporal evolution of the maximum fluid velocity Vmax and the mean Nusselt number Nuave at Ra = 104, T = 400, Pr = 10, M = 10. It is observed clearly that the fully developed flow is unsteady, and the time series is displaying what appears to be sustained single mode oscillatory behavior. Meanwhile, to better understand oscillatory behaviors of the current electro-thermo- convection, the phase spaces trajectory of velocity

at the sampling point A as well as the Fourier frequency spectrum of maximum velocity analysis are presented in Fig. 5(c) and (d), respectively. It can be seen that the phase portrait shows a slender eight-shaped closed orbit which represent a mono-periodic solution, and the flow is periodic self-sustained oscillation with fundamental frequency f1 = 11.9. In order to visualize the oscillatory behaviors, in Fig. 6 we have presented the evolution of the simulated charge density distributions, temperature fields and stream function contours at the special point of time, i.e. the wave crest and the wave trough of Nuave in Fig. 5(b). As we can see, the charge density distributions is consist of three stronger electro-plumes in the upper region and a

Fig. 11. The influences of Pr on the oscillatory convection, Ra = 104, T = 400, M = 10, Pr from 5 to 120: (a) Temporal evolution of the mean Nusselt number Nuave; (b) the phase spaces trajectory of velocity at the sampling point A.

Fig. 12. The influences of M on the oscillatory convection, Ra = 104, T = 400, Pr = 10: Time-evolution of Vmax and Nuave for (a) M = 5 and (b) M = 15; Fourier frequency spectrum of Vmax for (c) M = 5 and (d) M = 15.

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

weaker electro-plume in the lower region, and the thermal plumes always share the same shape and position with electro-plumes. The flow motion is characterized by four pairs of counterrotating cells, two smaller at the top and two larger ones at the bottom. Besides, in order to highlight the differences between the fields at different time points, a comparison is made in Fig. 6 (b2) and (b3), where the solid line represents the field at the wave crest, and the dash-dotted line represents the field at the wave trough. It is clear that the upper part of the annulus shows more complicated oscillatory behavior compared with the lower part. In fact, the charge density distributions and the flow fields at the bottom are almost not changed throughout the oscillation period. For clearer visualization of the non-stationary nature of the fields, take the temperature fields as example, we have shown the partial enlarged detail of the upper right part of Fig. 6(b2). As shown in Fig. 7, within one cycle, the thermal-plumes at both sides swing up and down under the influence of Coulomb buoyant force, and then the oscillatory behaviors repeat consistently. Accordingly, the smaller counter-rotating cells expand and shrink periodically.

929

time-average Nusselt number for sufficiently long time is found to be 3.18, which is lower than that for T = 400. From the results above, we can see that the stable flow may bifurcate into different patterns as T and Ra increase. In order to clarify the parameter ranges of different flow regimes, a series of numerical simulations in a wider range of Ra and T (5  102  Ra  5  105, 100  T  700) are carried out. Fig. 10 shows the transition diagram for Pr = 10, M = 10, R⁄ = 0.3 and C = 0.4. It is seen that the solution region is divided into three sub-regions: steady (marked by square), oscillatory (marked by circle) and chaotic (marked by triangle). The steady flow occurs only at lower T, and the region expands gradually with the increase of T. While the onset of chaotic motion occurs only at higher T, and the region expands gradually with the decrease of T. The region between them is oscillatory, note that for Ra  105, no oscillatory flow is observed.

4.2. Influences of Rayleigh number Ra and electric Rayleigh number T In this part, the influences of the Rayleigh number Ra and the electric Rayleigh number T on the oscillatory flows of the electrothermo-convection are investigated. Fig. 8(a) and (b) illustrate the time evolution of the average Nusselt number and the phase spaces trajectory of velocity at the sampling point A for T = 400, Pr = 10 and M = 10 with different values of Ra (Ra = 1  103, 5  103 and 1  104). It is found that the Rayleigh number has only a little effect on the oscillatory behaviors within the current range, the phase spaces trajectories for different Ra are very similar to each other and the fluctuation amplitudes of heat transfer show little difference. Actually, the Fourier frequency spectrum of maximum velocity analysis indicate that the fundamental frequencies for Ra = 1  103, 5  103 and 1  104 are 12.6, 12.4 and 11.9, respectively, which means that the increase of Ra will results in a slight decrease in fundamental frequency. However, on increasing the Rayleigh number to Ra = 5  104 the sub-critical flow state will be broken, and the time-periodic oscillatory convection gradually bifurcates to a non-periodic solution after some unpredictable time, see Fig. 8(c). If the Rayleigh number is further increased, the influence of thermal buoyant force becomes more and more strong, and the convection will directly turns into chaotic at high Ra (Ra = 105), see Fig. 8(d). Fig. 9 shows the influences of electric field on the oscillatory behaviors of electro-thermo-convection for Ra = 104, Pr = 10 and M = 10. As we can see from Fig. 9(a), for the low value T = 100, the flow is in a steady state and the heat transfer rate remain at 1.62; while for higher value of T, the flow begins to show obvious periodic oscillatory character, the fundamental frequency and the fluctuation amplitude increase with increasing of T, and the time-average Nusselt number in one cycle of oscillation for T = 200, 300 and 400 are 2.57, 3.03 and 3.28, respectively, which means the heat transfer is strengthen due to the increase of electric field strength. Moreover, the phase spaces trajectory of velocity at the sampling point A for different values of T is presented in Fig. 9 (c)–(f). Combined with Fig. 5(c), we can see that the attractor morphology of steady stationary solution eventually reaches a point at T = 100. However, for higher values of T, the phase portraits showing a variety of closed orbits with different shapes and locations, which means all of these flow oscillations are mono-periodic solutions, while with different oscillatory behaviors. As T increases to 500, the regular temporal pattern disintegrates, ultimately nonperiodic chaotic convection is developed after intermittent steady flow, see Fig. 9(b). And the very irregular phase space trajectory in Fig. 9(f) also clearly depicts its chaotic behavior. By calculating, the

Fig. 13. Time-evolution of Vmax for Ra = 104, T = 400, Pr = 10, M = 10: (a) R* = 0.1; (b) R* = 0.2; (c) R* = 0.4.

930

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

Fig. 14. Distributions of charge density, temperature fields and stream functions (from left to right) under different flow state for Ra = 104, T = 400, Pr = 10, M = 10: (a) steady and (b) chaotic state for R* = 0.1; (c) oscillatory and (d) chaotic state for R* = 0.2; (e) chaotic state for R* = 0.4.

4.3. Influences of fluid property parameters Pr and M Now we study the influences of fluid properties on oscillatory electro-thermo-convection, which can be characterized by the Prandtl number Pr and the non-dimensional mobility parameter

M. Fig. 11 present the temporal evolution of the mean Nusselt number Nuave and the phase spaces trajectory of velocity at the sampling point A at Ra = 104, T = 400 and M = 10 with different Pr. From Fig. 11(a) we can see that with the increase of Pr, the fluctuation amplitudes of heat transfer are decreased while the heat

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

transfer rates are increased obviously. However, the fundamental frequencies of oscillatory flow for Pr = 5, 40, 80 and 120 are 12.5, 12.4, 12.4 and 12.5, respectively, which means that the Pr has little influence on the oscillating frequency. Besides, the phase spaces trajectories for different Pr almost remain the same, see Fig. 11 (b). Therefore, the simulated results indicate that the increasing Pr has obvious influence on the heat transfer enhancement, but almost no influence on the flow motion. In Fig. 12, we have summarized the influences of M on the oscillatory convection for Ra = 104, T = 400 and Pr = 10. The temporal evolution of the maximum fluid velocity Vmax and the mean Nusselt number Nuave at M = 5 is presented in Fig. 12(a). It is note that for M = 5 the convection is chaotic, and this can be also seen in the Fourier frequency spectrum of Vmax (Fig. 12(c)) with broad-banded spectra. However, as mentioned above, for M = 10 the oscillatory convection changes from chaotic into periodic. If the value of M is further increased, the fluctuation signal evolves from periodic oscillation into quasi-periodic oscillation as shown in Fig. 12(b), with the Fourier transform showing a dominant peak at frequency of 6.27, has a number of other discrete peaks at a range of frequencies, see Fig. 12(d). Combined with Figs. 5 and 11, it is easy to find that the increasing M will lead to not only the enhancement of heat transfer, but also the increase of flow intensity, and the change of flowing structure as well. 4.4. Influence of the radii ratio R⁄ We also investigate the influence of the radii ratio of inner and outer cylinders on oscillatory behaviors and the results will be given in this section. Varying the radii ratio R⁄ from 0.1 to 0.4, while other parameters are fixed at the T = 400, Ra = 104, Pr = 10 and M = 10. Fig. 13 shows the temporal evolutions of the Vmax for R⁄ = 0.1, 0.2 and 0.4, respectively. It is clear that the flow intensity decreases with the increases of R⁄. Specially, for R⁄ = 0.1, the maximum fluid velocity Vmax reaches a plain after an initial transient and the flow appeared to briefly in a stable state. The corresponding distributions of charge density, temperature fields and stream functions are provided in Fig. 14(a). Only two electro and thermal plumes are clearly observed, and the flowing structure with two pairs of counter-rotating cells is strictly symmetric about the vertical centerline. Remarkably, in the lower region, the thermal convection and electroconvection behave to be a competition mechanism. Therefore flow is potential unstable and ultimately turn into a chaotic state, the lower electro and thermal plumes disappear and the streamline evolves into two twisted vortexes, as shown in Fig. 14(b). When the value of R⁄ is up to 0.2, the maximum fluid velocity presents periodic oscillation with time varying before the flow ends with a chaotic convection, see Fig. 13(b). To intuitively illustrate this transformation of flow pattern, the snapshots of the periodic oscillatory solution and the chaotic solution are presented in Fig. 14(c) and (d), respectively. Compared with the steady result for R⁄ = 0.1, for the oscillatory solution the upper convection rolls increase their size and strength. However, for the chaotic solution, a new electro plume imping towards the outer cylinder appears in the upper region, accompanied by the generation of an extra pairs of convection cells. Finally, for R⁄ = 0.4 the convection is totally chaotic, and the relatively low frequency non-periodic fluctuation of the Vmax is shown in Fig. 13(c). Meanwhile, the instantaneous results are presented in Fig. 14(e), where the increase in the number of electro and thermal plumes and convection cells are observed, and the flow is more complicated. To illustrate the influences of radii ratio R⁄ on the heat transfer, the temporal evolution of Nuave for the above simulations is provided in Fig. 15. Combined with Fig. 13, it is interesting to find that despite decreased flow intensity, the heat transfer rates increase

931

Fig. 15. Time-evolution of Nuave for Ra = 104, T = 400, Pr = 10, M = 10, R* from 0.1 to 0.4.

with the radii ratio. For R⁄ from 0.1 to 0.4, the corresponding time-average Nusselt numbers for the nondimensional time t’ at the same interval of 1.5–5 are 1.97, 2.72, 3.28 and 3.43, respectively. This may be explained by the fact that a narrower gap between the cylinders will results in a larger impact surface (see Fig. 14), which would contributes to heat transfer enhancement.

5. Conclusions Systematic numerical investigations of oscillatory electrothermo-convection in an eccentric annulus filled with insulating liquid have been carried out for the first time. The liquid is under the simultaneous action of the buoyant force and the Coulomb forces due to a strong unipolar charge injection. A unified LBM model developed recently is adopted to solve the whole set of macroscopic governing equations for this coupled multi-field problem. Meanwhile, phase space trajectories and Fourier frequency spectrum method of nonlinear analysis are also implemented to investigate the oscillatory convection. After code validation, we present the evolution of charge density distributions, temperature fields and streamlines for Ra = 104, T = 400, Pr = 10 and M = 10 to show the typical characteristics of the oscillatory behavior. The phase portrait and Fourier frequency analysis of velocity indicate that the convection is mono-periodic self-sustained oscillation with fundamental frequency f1 = 11.9. Then, by changing the governing parameters Ra, T, Pr and M, we investigate the effects of electric field, temperature gradient and fluid properties on the oscillatory behaviors of electro-thermoconvection. The results show that the values of Ra within the range from 103 to 104 have little effect on either heat transfer or flow motion of the periodic oscillatory convection. For even higher values of Ra, the balance between the Coulomb force and the buoyant force will be broken, and the flow ultimately turns into chaotic directly at Ra = 105. However, it is found that the value of T has significant impact on the oscillatory convection. As T increases, the flow develops from steady to periodic oscillatory, and finally evolves into chaotic. For the periodic oscillatory flow, the fundamental frequency and the fluctuation amplitude as well as the heat transfer rate increase with increasing of T. The results also show that the increasing Pr will lead to a decrease of fluctuation amplitude and improve the heat transfer, while has little influence on the oscillating frequency. Moreover, the flow undergoes chaotic state, periodic oscillation and quasi-periodic oscillation when M is varied from 5 to 15, accompanied by the increase of heat transfer rate and flow intensity.

932

T.-F. Li et al. / International Journal of Heat and Mass Transfer 134 (2019) 920–932

Furthermore, we also study the effects of the radius ratio on the oscillatory flow characteristics. It is shown that the formation of oscillating flows is strongly related to the value of radius ratio. For current parameter conditions, almost all of these systems under consideration eventually evolve into chaotic flow, except for the case with R⁄ = 0.3, in which the stable periodic oscillatory solution is obtained. Besides, with the increase of R⁄, the number of electro plumes and convection rolls increase, and an enhancement of heat transfer can be observed. Conflict of interest The authors declared that there is no conflict of interest. Acknowledgements The authors are grateful to Prof. Jian Wu for helpful suggestions. This work is supported by the National Natural Science Foundation of China (Grant No. 51776054), and partially by financial support from National Postdoctoral Program for Innovative Talents of China (Grant No. AUGA4150000117). References [1] S. Laohalertdecha, P. Naphon, S. Wongwises, A review of electrohydrodynamic enhancement of heat transfer, Renew. Sustain. Energy Rev. 11 (2007) 858–876. [2] J. Cotton, A.J. Robinson, M. Shoukri, J.S. Chang, A two-phase flow pattern map for annular channels under a DC applied voltage and the application to electrohydrodynamic convective boiling analysis, Int. J. Heat Mass Transfer 48 (2005) 5563–5579. [3] J.E. Bryan, J. Seyed-Yagoobi, Electrohydrodynamically enhanced convective boiling: relationship between electrohydrodynamic pressure and momentum flux rate, J. Heat Transfer 122 (1999) 266–277. [4] M. Ghalamchi, A. Kasaeian, M. Ghalamchi, et al., Optimizing of solar chimney performance using electrohydrodynamic system based on array geometry, Energy Convers. Manage. 135 (2017) 261–269. [5] J. Seyed-Yagoobi, J.E. Bryan, Enhancement of heat transfer and mass transport in single-phase and two-phase flows with electrohydrodynamics, Adv. Heat Transfer 33 (1999) 95–186. [6] L. Léal, M. Miscevic, P. Lavieille, et al., An overview of heat transfer enhancement methods and new perspectives: Focus on active methods using electroactive materials, Int. J. Heat Mass Transfer 61 (2013) 505–524. [7] J.L. Fernández, R. Poulter, Radial mass flow in electrohydrodynamicallyenhanced forced heat transfer in tubes, Int. J. Heat Mass Transfer 30 (1987) 2125–2136. [8] J.S. Paschkewitz, D.M. Pratt, The influence of fluid properties on electrohydrodynamic heat transfer enhancement in liquids under viscous and electrically dominated flow conditions, Exp. Therm. Fluid Sci. 21 (2000) 187–197. [9] G.F. Smirnov, V.G. Lunev, Heat Transfer during condensation of vapour of dielectric in electric field, Appl. Electr. Phenomena USSR 2 (1978) 37–42. [10] J. Trommelmans, J.P. Janssens, F. Maelfait, et al., Electric field heat transfer augmentation during condensation of nonconducting fluids on a horizontal surface, IEEE Trans. Industry Appl. 2 (1985) 530–534. [11] A. Castellanos, P. Atten, Numerical modeling of finite amplitude convection of liquids subjected to unipolar injection, IEEE Trans. Industry Appl. 5 (1987) 825–830. [12] P.A. Vázquez, A. Castellanos, Numerical simulation of EHD flows using Discontinuous Galerkin finite element methods, Computers Fluids 84 (2013) 270–278. [13] P. Traoré, A.T. Pérez, D. Koulova, et al., Numerical modelling of finiteamplitude electro-thermo-convection in a dielectric liquid layer subjected to both unipolar injection and temperature gradient, J. Fluid Mech. 658 (2010) 279–293. [14] K. Dantchi, T. Philippe, R. Hubert, et al., Numerical simulations of electrothermo-convection and heat transfer in 2d cavity, J. Electrost. 71 (2013) 341– 344. [15] J. Wu, P. Traoré, A finite-volume method for electro-thermoconvective phenomena in a plane layer of dielectric liquid, Numer. Heat Transfer, Part A: Appl. 68 (2015) 471–500. [16] H. Togun, T. Abdulrazzaq, S.N. Kazi, et al., A review of studies on forced, natural and mixed heat transfer to fluid and nanofluid flow in an annular passage, Renew. Sustain. Energy Rev. 39 (2014) 835–856. [17] J. Wu, P. Traoré, M. Zhang, et al., Charge injection enhanced natural convection heat transfer in horizontal concentric annuli filled with a dielectric liquid, Int. J. Heat Mass Transfer 92 (2016) 139–148.

[18] W. Hassen, H.F. Oztop, L. Kolsi, et al., Analysis of the electro-thermoconvection induced by a strong unipolar injection between two concentric or eccentric cylinders, Numer. Heat Transfer, Part A: Appl. 71 (2017) 789–804. [19] D.V. Fernandes, H.D. Lee, S. Alapati, et al., Numerical simulation of the electroconvective onset and complex flows of dielectric liquid in an annulus, J. Mech. Sci. Technol. 26 (2012) 3785–3793. [20] J. Wu, P.A. Vázquez, P. Traoré, et al., Finite amplitude electroconvection induced by strong unipolar injection between two coaxial cylinders, Phys. Fluids 26 (2014) 124105. [21] N. Agrait, A. Castellanos, Linear convective patterns in cylindrical geometry for unipolar injection, Phys. Fluids A: Fluid Dynam. 2 (1990) 37–44. [22] D.V. Fernandes, H.D. Lee, S. Park, et al., Electrohydrodynamic instability of dielectric liquid between concentric circular cylinders subjected to unipolar charge injection, J. Mech. Sci. Technol. 27 (2013) 461–467. [23] K. Luo, J. Wu, H.L. Yi, et al., Lattice Boltzmann model for Coulomb-driven flows in dielectric liquids, Phys. Rev. E 93 (2016) 023309. [24] K. Luo, J. Wu, H.L. Yi, et al., Lattice Boltzmann modelling of electro-thermoconvection in a planar layer of dielectric liquid subjected to unipolar injection and thermal gradient, Int. J. Heat Mass Transfer 103 (2016) 832–846. [25] K. Luo, J. Wu, H.L. Yi, et al., Numerical investigation of heat transfer enhancement in electro-thermo- convection in a square enclosure with an inner circular cylinder, Int. J. Heat Mass Transfer 113 (2017) 1070–1085. [26] K. Luo, T. Li, J. Wu, et al., Electro-thermo-convective flow of a dielectric liquid due to nonautonomous injection of charge by an elliptical electrode, Int. J. Heat Mass Transfer 127 (2018) 373–384. [27] J. Wu, P. Traore, C. Louste, et al., Effect of the mobility parameter on the oscillatory electroconvection of dielectric liquids subject to strong unipolar charge injection, IEEE Trans. Industry Appl. 50 (2014) 2306–2313. [28] A.N. Mordvinov, B.L. Smorodin, Electroconvection under injection from cathode and heating from above, J. Exp. Theor. Phys. 114 (2012) 870–877. [29] P.H. Kao, R.J. Yang, Simulating oscillatory flows in Rayleigh-Benard convection using the lattice Boltzmann method, Int. J. Heat Mass Transfer 50 (2007) 3315– 3328. [30] N. Williamson, S.W. Armfield, M.P. Kirkpatrick, et al., Bifurcation of natural convection flow in an inclined differentially heated closed square cavity, Comput. Therm. Sci.: Int. J. 7 (2015) 417–425. [31] L. Martínez-Suástegui, C. Treviño, J.C. Cajas, Thermal nonlinear oscillator in mixed convection, Phys. Rev. E 84 (2011) 046310. [32] I.E. Sarris, S.C. Kakarantzas, A.P. Grecos, et al., MHD natural convection in a laterally and volumetrically heated square cavity, Int. J. Heat Mass Transfer 48 (2005) 3443–3453. [33] J. Wang, M. Yang, Y. Zhang, Onset of double-diffusive convection in horizontal cavity with Soret and Dufour effects, Int. J. Heat Mass Transfer 78 (2014) 1023– 1031. [34] T. Jones, Electrohydrodynamically enhanced heat transfer in liquids – a review, Adv. Heat Transfer 14 (1978) 107–148. [35] A. Castellanos, Entropy production and the temperature equation in electrohydrodynamics. Dielectrics and electrical insulation, IEEE Trans. Dielectrics Electr. Insulat. 10 (2003) 22–26. [36] P. Atten, Electrohydrodynamic instability and motion induced by injected space charge in insulating liquids, IEEE Trans. Dielectrics Electr. Insulat. 3 (1996) 1–17. [37] W.M. Rohsenow, Y.I. Cho, Handbook of Heat Transfer, McGraw-Hill, New York, 1998. [38] J.C. Lacroix, P. Atten, E.J. Hopfinger, Electro-convection in a dielectric liquid layer subjected to unipolar injection, J. Fluid Mech. 69 (1975) 539–563. [39] A.T. Pérez, A. Castellanos, Role of charge diffusion in finite-amplitude electroconvection, Phys. Rev. A 40 (1989) 5844. [40] Q. Li, K.H. Luo, Q.J. Kang, et al., Lattice Boltzmann methods for multiphase flow and phase-change heat transfer, Progr. Energy Combust. Sci. 52 (2016) 62–105. [41] P. Lallemand, L.S. Luo, Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546. [42] C.K. Aidun, J.R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439–472. [43] Y. Qian, D. d’Humières, P. Lallemand, Lattice BGK models for Navier-Stokes equation, EPL (Europhysics Letters) 17 (1992) 479. [44] Z. Guo, B. Shi, N. Wang, Lattice BGK Model for Incompressible Navier-Stokes equation, J. Comput. Phys. 165 (2000) 288–306. [45] Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E 65 (2002) 046308. [46] Z. Guo, C. Zheng, B. Shi, An extrapolation method for boundary conditions in lattice Boltzmann method, Phys. Fluids 14 (2002) 2007–2010. [47] T.H. Kuehn, R.J. Goldstein, An experimental study of natural convection heat transfer in concentric and eccentric horizontal cylindrical annuli, J. Heat Transfer 100 (1978) 635–640. [48] K. Luo, H.L. Yi, H.P. Tan, et al., Unified lattice Boltzmann method for electric field-space charge coupled problems in complex geometries and its applications to annular electroconvection, IEEE Trans. Industry Appl. 53 (2017) 3995–4007.