The Characteristics of the Electro-Thermo-Convective Flow of a Dielectric Liquid Analyzed through the Electric Nusselt Number ^†

Abstract

This paper presents a fundamental study of electro-thermo-convective flows within a layer of dielectric liquid subjected to both an electric field and a thermal gradient. A low-conductivity liquid enclosed between two horizontal electrodes and subjected to unipolar charge injection is considered. The interplay between electric and thermal fields ignites complex physical interactions within the flows, all governed by a set of coupled electro-thermo-hydrodynamic equations. These equations include Maxwell, Navier–Stokes, and energy equations and are solved numerically using an in-house code based on the finite volume method. Electro-thermo-convective flows are driven by two dimensionless instability criteria: Rayleigh number Ra and the stability parameter T, and also by the dimensionless mobility parameter M and Prandtl number Pr. The electric Nusselt number (Ne) analogue to the Nusselt number (Nu) in pure thermal problems serves as an indicator to monitor the shift from a thermo- to an electro-convective flow and its eventual evolution into unsteady, and, later, chaotic flow. This change in regime is observed by tracking the electric Nusselt number’s behavior as a function of the stability parameter (T), for different values of the non-dimensional parameters (M, Ra, and Pr). The important role of mobility parameter M for the development of the flow is shown. The flow structure during different development stages in terms of the number of convective cells is also discussed.

1. Introduction

The combined effects of an electric field and a thermal gradient applied simultaneously to a dielectric liquid containing a space charge lead to very complex physical interactions in the flow. The electroconvection (EC) and electro-thermo-convection (ETC) in dielectric liquids are problems that have been intensively investigated in recent years [1,2,3,4]. Their study is essential for applications with electric fields in heat transfer and flow stability control processes, as well as for a thorough understanding of systems subjected to various natural convective forces (atmospheric convection, solar magneto-convection).

The space charge in a dielectric liquid can be created by unipolar injection that induces EC under the action of the Coulomb force. In the case of a planar electrode geometry, EC presents an analogy with the classical Rayleigh–Benard problem of a horizontal fluid layer heated from below. In both cases, the flow is driven by the body force (thermal buoyancy force for Rayleigh–Benard convection (RBC) and Coulomb force for EC). In the stability analysis, for weak driving forces, both systems remain at rest due to the viscous damping. Above some critical values, Ra_c and T_c, the flow motion takes place and evolves into a regular pattern of convective cells.

It is well known that the linear bifurcation of RBC is supercritical, while the linear bifurcation of EC is subcritical [4,5]. The subcritical bifurcation was first explained by Felici [6] for weak injection, later for SCL (space-charge-limited) injection in [7] and confirmed experimentally in [8,9]. The subcritical bifurcation is a representative characteristic of EC, featured by a linear stability criterion T_c, a non-linear stability criterion T_f, and a hysteresis loop between these two criteria [8,10]. Recently, Luo et al. showed that a lattice Boltzmann model (LBM) could predict the linear and finite amplitude stability criteria of the subcritical bifurcation in the EC flow [11,12,13] for both 2D and 3D flow scenarios.

In order to explain the different regimes of instability in a pure EC problem, it is convenient to introduce the electric Nusselt number (Ne), which is defined as the ratio of the effective current and the current existing without liquid motion. It can be considered as the analog of the Nusselt number (Nu) for a pure thermal problem but has not been studied yet in an ETC problem.

In this article, we consider the convective phenomena in an insulating liquid layer enclosed between two parallel electrodes and subjected to the combined effects of an electric field and a thermal gradient. The effects of both fields simultaneously applied to the fluid layer leads to very complex physical interactions in the flow; these have received much attention by many authors [14,15,16]. When the space charge only results from ion injection, the coupling between the conservation equations of momentum, electric charge, and energy is ensured via the Coulomb and buoyancy forces. This coupling results from the direct interaction between the velocity, temperature, and charge perturbations as well as from the interaction between the velocity and the electric field. Due to the complexity of the mathematical model, the problem of ETC flows can be investigated mainly by numerical simulations; different methods have been developed and are shown in Refs. [17,18,19,20]. For example, an ETHD system with solid–liquid phase change in both rectangular and annular cavities have been analyzed by LBM simulations in Ref. [20]. For SCL injection, EC produces a significant increase in charge transfer [21]. Additional studies on ETHD deal with the evaluation of heat transfer enhancement from numerical [21,22,23,24] and experimental [25,26] point of view.

The transition from laminar to chaotic flow in the electro-thermal convection of a dielectric fluid for different Rayleigh numbers Ra was investigated in Ref. [27], using a lattice Boltzmann method. By increasing Ra, the different means of transition to chaos were observed and the chaotic behavior was quantitatively analyzed.

Some preliminary results of electro-thermo-convective flow and its eventual evolution into unsteady and later chaotic flow from the point of view of the electric Nusselt number were obtained in Ref. [28]. Here, we present a detailed and extended study of the behavior of the electric Nusselt number (Ne) on a large time scale as a function of the stability parameter T at various values of the following non-dimensional parameters: M = 5, 10, 20, 40; Ra = 0, 2.10³, 5.10³, 10⁴ and for Pr = 10. The main purpose of the present study is to explore the subcritical feature of ETC from the point of view of the electrical Nusselt number Ne and to investigate the effect of mobility M on Ne. The entire set of ETC equations is solved by a direct numerical simulation based on the finite volume method [29,30].

2. Statement of the Problem

2.1. Basic Equations and Parameters

We consider an incompressible and perfectly insulating liquid of density ρ, permittivity ε, the ionic mobility K, constant kinematic viscosity υ and constant thermal diffusivity κ, enclosed between two parallel planar electrodes of length L. The liquid layer of depth H is subjected to a potential difference ΔV = V₀ − V₁ and a thermal gradient $\Delta \theta ={\theta }_1-{\theta }_0$ (see Figure 1).

The z-axis is taken perpendicular to the electrodes and the gravity is considered as acting in the negative z-direction. In this study, as the liquid is perfectly insulating, we do not take into account the ion creation due to the dissociation process that can occur in conductive liquids. Here, only the homogeneous and autonomous injection of the unipolar charge at the emitter electrode is considered. This means that q = q₀ at z = 0 at all times, i.e., that the injector and, therefore, the injection rate are not influenced by the electric field nor by the liquid motion. These injected ions are supposed to have a given ionic mobility K. The lower electrode is chosen as the source of the ions, which are then injected into the fluid and collected at the upper electrode. The liquid is heated from below. The general set of equations for an incompressible fluid expressing the conservation of mass and momentum (Navier–Stokes equations) including electrical and buoyancy effects, the energy balance under Boussinesq assumption, the charge density conservation, Gauss theorem, and the definition of the electric field towards electric potential V, takes the following form:

$$\nabla .\overrightarrow{u}=0$$

(1)

$${\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+(\overrightarrow{u}.\nabla )\overrightarrow{u}=-\nabla \tilde{p}+\Delta \overrightarrow{u}+RTCq\overrightarrow{E}+{\displaystyle \frac{Ra}{Pr}}\theta {\overrightarrow{e}}_z$$

(2)

$${\displaystyle \frac{\partial \theta }{\partial t}}+\overrightarrow{u}.\nabla \theta ={\displaystyle \frac{1}{Pr}}\Delta \theta$$

(3)

$${\displaystyle \frac{\partial q}{\partial t}}+\nabla .\left(q(\overrightarrow{u}+R\overrightarrow{E})\right)=0$$

(4)

$$\Delta V=-Cq, \overrightarrow{E}=-\nabla V$$

(5)

where $\overrightarrow{u}$ is the velocity, p the static pressure, q the electric charge density in the bulk, $\overrightarrow{E}$ the electric field, θ the absolute temperature and V the electric potential, $q\overrightarrow{E}$ is the Coulomb force. The following non-dimensional scales are introduced: H for length, the applied voltage ΔV = V₀ − V₁ for electric potential, ΔV/H for the electric field, q₀ for the charge density, ε₀KΔV²/H³ for the current density, υ/H for the velocity, ρ₀υ²/H² for the pressure, θ for the temperature and H²/υ for time (ρ₀, K₀ and ε₀ are, respectively, the density, the ionic mobility and the permittivity at the reference temperature θ₀).

The system is essentially governed by the following non-dimensional parameters:

$$Ra={\displaystyle \frac{g\beta \Delta \theta H^3}{\nu \kappa }}, T={\displaystyle \frac{{\epsilon }_0\Delta V}{{\rho }_0\nu K_0}}, C={\displaystyle \frac{q_0{H}^2}{{\epsilon }_0\Delta V}}, M={\displaystyle \frac{1}{K_0}}{\left({\displaystyle \frac{{\epsilon }_0}{{\rho }_0}}\right)}^{1/2}, Pr={\displaystyle \frac{\nu }{\kappa }}, R={\displaystyle \frac{T}{M^2}},$$

where Ra is the Rayleigh number (β thermal expansion coefficient), T is the electric instability parameter, defined as the ratio of the driving force and the damping viscous force, C is a measure of the injection level (q₀ the injected charge), g is the gravity, Pr is the Prandtl number, and R is the electrical Reynolds number. M is the non-dimensional mobility parameter and it only depends on the physical properties of the fluid. The name “mobility”, given by Felici in 1969 [6], considering the conversion of the electrical energy (1/2εE²) into mechanical energy (1/2ρw‘), developed the concept of an electrohydrodynamic mobility (ε/ρ)^1/2. Experimental measurements showed that w’~1/3 (ε/ρ)^1/2/E [3]. The ratio between the EHD mobility and the true ion mobility K is the parameter M; for dielectric liquids its typical values are greater than three [4].

2.2. Initial and Boundary Conditions

The boundary conditions on velocity, temperature, electric potential, and charge density are chosen, as shown in Figure 2. On the horizontal and lateral walls, we apply no-slip boundary conditions. We consider here the fluid initially at rest with simultaneous heating and injection at t = 0.

The no-slip condition at the impermeable, thermally, and electrically perfectly conducting electrodes is assumed and the horizontal metallic electrodes are assumed to be rigid. The injection at the emitter is autonomous, which means that the injection rate is not influenced by any perturbation in the bulk and, therefore, that q = 1 at z = 0.

2.3. Numerical Method

The set of coupled partial differential Equations (1)–(5) are integrated using the finite volume method [29]. The numerical procedure based on the augmented Lagrangian method. The calculation is fully transient. Special care is taken for the transport equation (4) for charge density q (because of its hyperbolic nature). The SMART algorithm was utilized in this way.

The boundary condition of the charge density on the collector requires some comments because of the hyperbolic nature of the charge conservation equation. From a mathematical point of view, this implies that a boundary condition on the injector is enough to determine the solution. However, in the numerical procedure, the discretized charge density transport equation involves the value of q at some neighboring nodes so that interior nodes next to the boundaries will require the value of q at those boundaries. This is where boundary conditions come into play, as they provide pre-defined values of q on the boundaries to close the system.

In summary, boundary conditions play a crucial role in solving the charge density transport equation by providing the necessary values of q at the boundaries. They ensure the closure of the system and enable accurate modeling and simulation of the behavior of the system near the boundaries. The selection and implementation of appropriate boundary conditions are important considerations in numerical simulations and modeling studies.

The overall algorithm ensures that nothing originates from the collector, which is both physically and mathematically consistent with the principle of the charge conservation equation.

The code has been thoroughly validated and numerical results have been compared with the analytical solution (hydrostatic solution) as well as with some results from linear and non-linear stability analyses [29,30]. Grid independency tests were conducted in order to select the optimal grid to undertake the computation.

3. Results and Discussion

The electro-thermo-convective flow depends on two instability criteria: the Rayleigh number Ra and the electrical parameter T, as well as on the non-dimensional parameters Pr and M number. By changing the governing parameters Ra and T, we investigate the influence of the mobility parameter M on the electro-thermo-convective flow from the point of view of the electric Nusselt number. The transition from steady to unsteady flow is also analyzed.

The global effect of convective heat transfer is usually presented by the Nusselt number Nu, which is the ratio of the mean heat flux to the flux that would exist without convection for the same temperature difference. In analogy with Nu, we introduce the electric Nusselt number Ne, given as Ne = I/I₀, where I is the electric current and I₀ is the electric current without flow motion (note that the current is transported mainly by electro-convection induced by the ion injection from the electrode at z = 0 and the dependence between the current I and the typical liquid velocity is non-linear).

In this section, we present the electro-thermo-convection (ETC) in the classical geometry of a dielectric liquid layer lying between two parallel planar electrodes and subjected to a strong unipolar injection of charges (C = 10). We consider a rectangular cavity of height H = 1 and length L = 2 in the case of heating and a strong charge injection from the lower electrode. The electric current can be computed by integrating the current density on the entire emitter or the collecting electrodes. Here, Ne is calculated on the injecting electrode.

The time evolution of the electrical Nusselt number for the lower instability parameter value (T = 120) is given in Figure 3a. At the beginning of the injecting and heating, we observe a brusque increase in the Ne number, which reaches an asymptotic value that shows a steady-state behavior. The steady state does not depend on the initial conditions. For Tc < T < 300, Ne reaches stationary values, whereas for T > 300, Ne fluctuates around the average values (Figure 3b), which is the sign of the existence of a transition in the flow from a steady laminar regime to an unsteady periodic one.

The flow moves from a motionless state to a convective flow and finally to a periodic one once the system loses its linear stability. This change in regime foreshadows the later appearance of chaotic flow for certain high values of T. It can be explained by the fact that the system with higher T values (T > 300) is essentially in the inertial state of motion and that the flow will probably never reach a steady state.

The results in the next paragraphs are presented for C = 10, Pr = 10 and for the values of Ra above the critical Rayleigh number (Ra_c~1708).

3.1. Influence of Rayleigh Number Ra on Electric Nusselt Number Ne

Figure 4 shows the evolution of the electric Nusselt number Ne with T in the case of ETC for different values of the Ra number (one is slightly above the critical value Ra = 2000 and the others are for higher Ra values Ra = 5000 and 10,000) and for the pure EC (Ra = 0).

It can be seen that the Ne number does not depend significantly on the Ra number for high enough values of T. We have already demonstrated a similar behavior for the Nusselt number (Nu) in Refs. [31,32]; for large values of T, the electrical forces dominate the buoyancy forces. In the case of pure EC, when T < T_c (T_c = 163.4 [11]), the perturbation does not trigger the flow instability; the perturbed flow reverses to the hydrostatic state. If T decreases after the EC vortices are formed, they are maintained until T = T_f (T_f = 108.5 [11]), when the system returns to the hydrostatic state. As soon as T is greater than T_c (T > T_c) there is a sudden jump in the Ne for EC due to the motion of the liquid initially at rest, with a maximum velocity that exceeds the ion drift velocity. For ETC, when T < T_c we observe greater values for the Ne, which can be explained by the non-linear coupling phenomena between the buoyancy and Coulomb forces; the motion of the liquid induced by thermo-convection drags more electric charge and then causes an increase in the electric current. This means that the increase in the electric Nusselt number Ne is due to thermo-convection. For Ra = 2000 and T > 300, the Ne values are close to the ones obtained in the case of pure EC. At Ra = 10,000 the situation is different; in the initial stage, the flow is thermo-convective with a two-roll structure, (it does not change), but later, at high values of T (T > 280), it changes drastically. There is a jump, and it becomes an electro-thermo-convective flow structure with four rolls.

3.2. Influence of Mobility Parameter M and T

In this section, we study the electro-thermo-convection driven by buoyancy and Coulomb forces for a large set of M parameter values. Numerical simulations are made for larger values of T to those corresponding with the saturation zone and to the fully developed turbulence flows. Figure 5 presents the variation in the electric Nusselt number Ne as a function of T under different mobility parameters.

The shape of the curve for M = 5 is similar to those obtained in Ref. [9] for pure electro-convection; Ne grows with T then reaches its saturation value. For other values of M, we can see that Ne increases with T until T~350 but does not reach the saturation value because of the influence of thermo-convection. For high values of T (T > 350), Ne decreases due to the transition from viscous to non-linear flow. For a high value of M (M > 40), i.e., liquids with very low ion mobility, the influence of T seems less important and the flow remains thermo-convective. For higher values of T (T > 450), Ne increases again. The increases in velocity and turbulent mixing can explain this.

In Figure 6, the variation of the electric Nusselt number Ne with M is presented for different values of the instability criteria T. We can see the fast linear growth of Ne, which is similar to the results obtained in Ref. [9].

Fom theoretical predictions and experimental observations [9], the saturation value of Ne depends on the physical property of the liquids and is given as a function of parameter M by the following expression:

$$Ne={\displaystyle \frac{I}{I_0}}\approx {\left({\displaystyle \frac{M}{3}}\right)}^{1/2}$$

(6)

Figure 7 presents the numerical results for the saturation value of Ne taken from Figure 5 and those obtained from formula (6) for pure electrohydrodynamic problems [9]. As one can see, the shape of both curves is similar, the difference between them being due to the influence of thermo-convection.

3.3. Flow Transition Process Described by the Maximum Fluid Velocity

Figure 8 presents the bifurcation diagram for thermo- and electro-convection expressed with the maximum fluid velocity in the bulk V_max versus T for different values of parameter M. Initially, V_max increases slowly with T. Then, for certain values of T, the slope of the curve changes and continues to increase. The change in slope depends on the values of M and corresponds to the transition from thermo- to electro-convection. For smaller values of M (M = 5 and M = 10, see Figure 8a), the values of V_max are much higher than for the other represented M (see Figure 8b), but the behavior of the curves is similar.

Figure 9 shows the evolution of the maximum velocity magnitude in time for different driving parameters M and T and Ra = 2000 and Pr = 10. Figure 9a plots the velocity magnitude for M = 10 and three different electric Rayleigh numbers, T = 150, 220 and 380. For T = 150 the flow remains inactive at the early stage; the velocity magnitude then grows quickly with an exponentially growing rate. Finally, the system reaches a steady state with two thermo-convective roll patterns. The curve for T= 220 has a slightly different behavior; it first develops as for T = 150 but after a certain time (~10) there is a sudden upward jump and it reaches a steady state with four convective cells. This means that the electro-convection dominates, but the flow rests in a steady state. For T = 380, the Vmax curve grows rapidly; the flow becomes strongly unsteady and looks periodic.

Figure 9b depicts the evolution in time of the maximum velocity magnitude at different settings of the driving parameters M and T corresponding to the values at which the slopes of the curves change (see Figure 8). The evolution of peak velocity Vmax also experiences two different stages after the initial disturbance. The first one is the exponential growth stage, which is the same for the smaller values of T (see Figure 9a). However, the second stage is quite different; Vmax grows rapidly, the flow becomes strongly unsteady and looks chaotic. In simpler terms, the flow transitions from a stationary state to a thermo-convective flow, followed by an electro-convective flow, and eventually reaches a chaotic state once the system loses its linear stability. This means that the system with high T is essentially in the inertial state of motion. Due to the electric force and destabilizing buoyancy force, the system exhibits electro-thermo-convective vortices and transitions to a chaotic flow field, where the flow fluctuates irregularly in time.

In Figure 10, we present the charge density distribution and the stream function at a steady convective state for the values of the electrical instability parameters T, in which the flow passes from a thermo-convective pattern with two rolls to an electro-convective one with four rolls.

In Figure 10, we can observe the formation of electro-convective cells with charge-depleted central regions. This absence of charges in the core of the cell is a consequence of the non-uniform distribution of the electric field and the rotational fluid motion induced by the Coulomb force. This well-established phenomenon is a hallmark of Coulomb-driven electro-convective flows and has been consistently reported in numerous experimental and numerical investigations [4,5,12,30].

In Figure 11, the transition from a thermo- to an electro-convective regime is shown in detail by the stream function. For M = 5 we see that the two-roll patterns for T < T_c change to four-roll patterns from T = 220, as in the pure electro-convective case. The same change in the roll patterns also occurs for M = 10, when four rolls appear at T approximatively equal to 220 and for M = 20, in which the change appears at T = 250. For M = 40 and M = 60, the two-roll pattern turns to a four-cell one at a higher level of T—T~320, and T~380, respectively.

The change in the number of rolls corresponds to a change in the slope of the curve (see Figure 8 and Figure 9a for T = 150 and 220). In all cases, the change in the slope of the curves Ne(T) always occurs when the number of rolls passes from two to four; in all cases Ne seems to saturate for the highest values of T. Such a transition from two to four convective cells has already been observed [31,32] in another study. We must also note that a change from two to four cells for M > 40 appears for higher values of T, which means that the flow remains thermo-convective up to T > 320. In the last column, the patterns are already in an unsteady and chaotic regime.

4. Conclusions

This work presents a numerical study of electro-thermo-convective phenomena in an insulating liquid layer between two parallel plates from the point of view of the electric Nusselt number. Numerical simulations for a large range of M and T values corresponding to the saturation zone and a fully developed chaotic flow were performed. It was confirmed once again that electrical forces dominate buoyancy forces as far as T > T_c. For higher values of T, and for a fixed M, the change in the number of rolls has an influence on the shape of the Ne(T) curves. It can be clearly seen that when the number of rolls changes, the slope of the maximum velocity with T changes as well. The transition from thermal to electro-convective flow was found and studied. We show the transition from thermo- to electro-convective flow, and, later, to unsteady flow as a function of parameter M and T through the maximum fluid velocity. For lower values of M (M = 5 and M = 10), electro-convection dominates for the values of T close to the linear instability criteria and Ne(T) grows fast. For high values of M, the electrical forces exceed the buoyancy forces for high T values. When T reaches the region in which the convective flow moves to the inertial state of motion and finally to chaos, Ne increases with T for a fixed M. Some of the results on Ne obtained here are in good agreement with the results presented in Ref. [9].