Stability and bifurcation of annular electro-thermo-convection

Abstract

We numerically investigated the global linear instability and bifurcations in electro-thermo-convection (ETC) of a dielectric liquid confined in a two-dimensional (2-D) concentric annulus subjected to a strong unipolar injection. Seven kinds of solutions exist in this ETC system due to the complex bifurcations, i.e. saddle-node, subcritical and supercritical Hopf bifurcations. These bifurcation routes constitute at most four solution branches. Global linear instability analysis and energy analysis were conducted to explain the instability mechanism and transition of different solutions and to predict the local instability regions. The linearized lattice Boltzmann method (LLBM) for global linear instability analysis, first proposed by Pérez et al. (Theor. Comput. Fluid Dyn., vol. 31, 2017, pp. 643–664) to analyse incompressible flows, was extended here to solve the whole set of coupled linear equations, including the linear Navier–Stokes equations, the linear energy equation, Poisson's equation and the linear charge conservation equation. A multiscale analysis was also performed to recover the macroscopic linearized Navier–Stokes equations from the four different discrete lattice Boltzmann equations (LBEs). The LLBM was validated by calculating the linear critical value of 2-D natural convection; it has an error of 1.39% compared with the spectral method. Instability with global travelling wave behaviour is a unique behaviour in the annulus configuration electrothermohydrodynamic system, which may be caused by the baroclinity. Finally, the chaotic behaviour was quantitatively analysed through calculation of the fractal dimension and Lyapunov exponent.

1. Introduction

Electrohydrodynamics (EHD) studies the complex interaction between an electric field and a dielectric fluid of very low conductivity. This configuration has broad applications in a range of industrial and biological devices. The EHD effects can be used to design microscale EHD pumps (Bart et al. 1990; Darabi et al. 2002), to design devices and drug delivery systems (Chakraborty et al. 2009) and DNA microarrays (Lee et al. 2006) and to design new strategies for active flow control (Bushnell & Mcginley 1989). As one of the most important applications, EHD combining the thermal effect, called electrothermohydrodynamics, has been applied to enhance heat transfer efficiency in the past 30 years (Laohalertdecha, Naphon & Wongwises 2007), such as in boiling and condensing systems (Cotton et al. 2005), drying and evaporating systems (Rashidi et al. 2017) and solar energy systems (Ghalamchi et al. 2017). An excellent review covering various aspects of electro-thermo-convection (ETC) can be found in Léal et al. (2013). It is notable that most previous studies on this topic concerned the simplest parallel-plane-electrode configuration. However, compared with the parallel geometry, the concentric configuration is not subjected to the lateral-wall effect, and its experimental set-up is more realizable, i.e. the concentrating solar collector receiver design (Togun et al. 2014). In the literature, there are few numerical works (Fernandes et al. 2012; Wu et al. 2016; Hassen et al. 2017; Li et al. 2019; He et al. 2021; Ma et al. 2022) on ETC adopting the configuration of a concentric annulus, which solely focused on the heat transfer enhancement. The fundamental fluid mechanics in this flow have not been clearly revealed, potentially inhibiting a more extensive application of ETC. Therefore, in this work, we systematically study the ETC in a concentric annulus, including its global linear instability analysis, bifurcations, routes of transition to chaos and heat transfer enhancement, to improve our general understanding of this complex flow. In the following, we summarize the works on linear stability analysis and bifurcations analysis of ETC and concisely discuss the position of the current work in the literature.

1.1. Linear instability analysis in ETC

As there is no work in the literature on the linear stability analysis of ETC in an annulus, we summarize the well-studied configuration of two parallel-plane electrodes. More specifically, we separately discuss ETC induced by the electric potential difference plus a destabilizing thermal gradient or a stabilizing thermal gradient. In the parallel-plane ETC with a destabilizing thermal gradient, Worraker & Richardson (1979) first considered that a thermal gradient is superimposed upon a unipolar charge injection of arbitrary strength by assuming that the ionic mobility varied linearly with temperature but that the electric permittivity remained constant. Castellanos, Atten & Velarde (1984) subsequently analysed the combined effect of Coulombic and polarization forces in the case of a weak unipolar injection by considering temperature-induced linear variations in both mobility and permittivity. Rodriguez-Luis, Castellanos & Richardson (1986) investigated the onset of steady convection under charge injection strengths of any finite magnitude. Pontiga & Castellanos (1994) comprehensively studied the flow stability in a parallel-plane ETC based on the force balance and noted the effect of residual conductivity on the stability of the liquid layer.

In the parallel-plane ETC with a stabilizing thermal gradient, Roberts (1969) initially performed a linear instability study on the ETC by assuming the variation in dielectric constant as a function of temperature. Pontiga, Castellanos & Richardson (1992) showed that the overstability in electrothermohydrodynamics (ETHD) is caused by the restoring forces provided by a thermal gradient. The indirect coupling between the charge density and the thermal field through mobility and permittivity induces the oscillation. Taraut & Smorodin (2012) found that the flow system evolves from the monotonic mode to the oscillatory mode via a transition in frequency. Mordvinov & Smorodin (2012) showed that for relatively weak heating from above (low Rayleigh number Ra, the ratio of the buoyancy to the viscous force), the initial perturbation either decays in an oscillatory manner (low electric Rayleigh number T, the ratio of Coulomb force to the viscous force) or grows monotonically (large T). For a large Ra, oscillatory growth of the initial perturbation occurs at large T. More recently, Guan et al. (2021) observed a two-stage bifurcation for overstability near the threshold Rayleigh number with a significant change in phase and amplitude.

1.2. Bifurcations and chaos in ETC

Subcritical bifurcation (Félici 1971; Lacroix, Atten & Hopfinger 1975; Zhang 2016) is an intrinsic property for EHD flow, which is characterized by an abrupt jump in the strength of the flow motion from zero to a finite value when finite-amplitude perturbations are present. Richardson (1980), Agrait & Castellanos (1990), Fernandes et al. (2013) and recently Wu et al. (2014) studied the influences of the injection strength, the ratio of the inner cylinder diameter D* to the gap width L* (A = D*/L*), and the mobility number M on T_c in the annulus geometry. Huang et al. (2021) numerically investigated the second bifurcation for higher parameter flows. However, the types of bifurcation for natural convection in an annulus are not fixed and rely on the Prandtl number and A. Yoo (1999) reported the occurrence of dual solutions (two kinds of steady solutions under the same parameter) at Rayleigh numbers larger than a critical value at A = 1.25. The hysteresis phenomenon (indicative of a subcritical bifurcation) occurs for fluids of 0.3 ≤ Pr ≤ 0.4, but it is not observed for 0.5 ≤ Pr ≤ 1 (saddle-node bifurcation). Mizushima, Hayashi & Adachi (2001) researched the influence of A on the bifurcation type and clarified the origin of the dual solutions from the bifurcation analysis. Serrano-Aguilera, Blanco-Rodríguez & Parras (2021) conducted abundant numerical investigations on dual solutions for values of the Prandtl number between 0.01 and 1 and Rayleigh numbers between 10² and 5 × 10⁶. Steady flow transition to oscillation (Hopf bifurcation) was also studied systematically. In the case where both electrical and thermal effects are considered, Traoré et al. (2010) gave an analytical model with a plane-electrode configuration to understand the appearance of subcritical or supercritical bifurcations for ETHD depending on the value of the Prandtl number and mobility parameter M. As we can see, for the annulus configuration, the bifurcation of the ETC predicted is rather complex. One of the focuses of this work is to further study the bifurcation in ETC and clarify the flow mechanisms underneath.

After the emergence of the ETC, as the strength of the electric field and thermal field further increase, various successive flow bifurcations can take place, such as periodic, quasiperiodic and even chaotic motions. The chaotic characteristics of the electro-convection (EC) have been confirmed in some experiments. Atten, Lacroix & Malraison (1980) showed that the EC became time dependent and chaotic above the instability threshold, and one characteristic frequency and its subharmonic were found in the power spectra of the total current fluctuations. Malraison & Atten (1982) found two types of chaotic behaviours in the power spectra of intensity fluctuations, namely, exponential decay in viscous-dominated flow and power-law decay in inertially dominated flow. Furthermore, the trajectories in an n-dimensional phase space have been reconstructed based on the experimentally obtained time series of the total current fluctuations with a time-delay technique (Malraison et al. 1983), and the fractal dimension of the chaotic attractor was calculated based on the Grassberger–Procaccia method (Grassberger & Procaccia 1983). However, this method becomes inapplicable when the logarithm of the integral correlation function does not show a defined slope. In such a case, Atten et al. (1984) concluded that the fractal dimension of the EC flow seems to increase without limit. By extracting time series from numerical simulation, Chicón, Pérez & Castellanos (2001) applied the algorithm proposed by Kantz & Schreiber (2004) to obtain a linear function S(Δn) whose coefficient of proportionality was the desired largest Lyapunov exponent. However, in these algorithms, some free parameters must be introduced, such as the embedding dimension and delay time. To characterize the chaotic motions in the EC in a more systematic way and avoid assigning free parameters, the method of Wolf et al. (1985) has been widely used recently (Feng et al. 2021; Huang et al. 2021). Moreover, Li et al. (2020) first introduced this method to investigate the transition process from laminar to chaotic flow in ETC.

1.3. The position and structure of the current work

The present work was motivated to investigate the bifurcation phenomena and transition routes to chaos in the ETHD system with an annulus configuration. Instability mechanisms can be explained by global linear instability analysis and energy analysis. The implementation of global linear instability analysis is based on a linearized lattice Boltzmann method (LLBM), which was first proposed by Pérez, Aguilar & Theofilis (2017). This method is extended to solve the coupled linear Navier–Stokes equations, linear thermal equation, linear charge potential equation and linear charge density equation in this work. With the instability and bifurcation mechanisms being illustrated clearly, we have also discussed the heat transfer enhancement in the ETHD flow in the annulus, which will help the expansion the applications of ETHD in engineering.

The remainder of the present paper is organized as follows. In § 2, we state the physics problem to be studied, the non-dimensional governing equations, derivations of our mathematical method and the basics of the global linear instability analysis. Section 3 is devoted to a detailed presentation of the model treatment. The numerical results are presented and discussed in § 4. The conclusions are presented in the last section.

2. Problem statement and governing equations

2.1. Governing equations

Figure 1 shows a schematic diagram of ETC in a two-dimensional concentric annulus due to charge injection and thermal buoyant force. The system consists of an outer cylinder with radius R_o, within which an inner cylinder with radius R_i is located in the centre, and the aspect ratio is fixed: A = 1.25. The inner and outer cylinders are kept at different constant temperatures θ_hot and θ_cold (Δθ ₀ = θ_hot − θ_cold), respectively, and a radial direct current electric field is established by applying a voltage difference Δϕ ₀ to the two electrodes. In this work, we assume that the fluid is incompressible, Newtonian, perfectly insulating and under the Boussinesq approximation. The free charges are only generated at the inner electrode and are then injected into the bulk of the liquid. A unipolar injection of the free charges is considered. It is also assumed that the injection is autonomous and homogenous (Laohalertdecha et al. 2007), which means that the injected charge density q ₀ at the inner cylinder remains constant and uniform.

Figure 1. Schematic diagram of ETC in a two-dimensional concentric annulus with A = 1.25. Red, yellow and blue are different regions divided.

The governing equations describing this problem include the Navier–Stokes equations, the temperature equation, the Poisson equation describing the electric potential and the charge conservation equation describing the charge transport. To further simplify the modelling, the viscous thermal dissipation, the magnetic effects and Joule heating are assumed to be negligible (Castellanos 1998). We introduce dimensionless quantities, denoted with a star, as follows:

(2.1)\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {x_i^\ast= \dfrac{x}{L},\quad u_i^\ast= \dfrac{u}{{\nu /L}},\quad {t^\ast } = \dfrac{t}{{{L^2}/\nu }},\quad {q^\ast } = \dfrac{q}{{{q_0}}},\quad {\phi^\ast } = \dfrac{\phi }{{\varDelta {\phi_0}}},}\\ {E_i^\ast= \dfrac{{{E_i}}}{{\varDelta {\phi_0}/L}},\quad {\theta^\ast } = \dfrac{{\theta - {\theta_{ref}}}}{{\varDelta {\theta_\textrm{0}}}},\quad {\rho^\ast } = \dfrac{\rho }{{{\rho_0}}},\quad {p^\ast } = \dfrac{p}{{{\rho_0}{{(\nu /L)}^2}}}} \end{array}}, \right\}\end{equation}

where L = R_o − R_i is the length scale. If we drop the superscript stars for clarity, the resulting governing equations of this ETC system are (Traoré et al. 2010)

(2.2)\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = \textrm{0,}\end{gather}

(2.3)\begin{gather}\frac{{\partial \boldsymbol{u}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{uu}) =- \boldsymbol{\nabla }p + {\mathrm{\nabla }^\textrm{2}}\boldsymbol{u} + \frac{{Ra}}{{Pr}}\theta {\boldsymbol{e}_y} + {\left( {\frac{T}{M}} \right)^\textrm{2}}qC\boldsymbol{E},\end{gather}

(2.4)\begin{gather}\frac{{\partial \theta }}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\theta \boldsymbol{u}) = \frac{1}{{Pr}}{\mathrm{\nabla }^\textrm{2}}\theta ,\end{gather}

(2.5)\begin{gather}{\mathrm{\nabla }^\textrm{2}}\phi =- Cq,\end{gather}

(2.6)\begin{gather}\boldsymbol{E} =- \boldsymbol{\nabla }\phi ,\end{gather}

(2.7)\begin{gather}\frac{{\partial q}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left[ {\left( {\frac{T}{{{M^2}}}\boldsymbol{E} + \boldsymbol{u}} \right)q} \right] = \frac{T}{{{M^2}}}\alpha {\mathrm{\nabla }^\textrm{2}}q.\end{gather}

In the above equations, u = (u_x, u_y) is the velocity, E = (E_x, E_y) is the electric field and e_y is the vertical direction. The scalars ρ, θ, ϕ and q denote the fluid density, temperature, electric potential and charge density, respectively. Six dimensionless parameters appear in governing equations, defined as

(2.8a–f) \begin{gather} Ra = \frac{{g\beta \varDelta {\theta _0}{L^3}}}{{\nu \chi }}\textrm{,}\quad Pr = \frac{\nu }{\chi },\quad T = \frac{{\varepsilon \varDelta {\phi _\textrm{0}}}}{{\mu K}},\quad C = \frac{{{q_0}{L^2}}}{{\varepsilon \varDelta {\phi _\textrm{0}}}},\quad M = \frac{1}{K}{\left( {\frac{\varepsilon }{{{\rho_0}}}} \right)^{1/2}},\notag\\ \alpha = \frac{D}{{K\varDelta {\phi _\textrm{0}}}}, \end{gather}

where the symbols μ, χ, $\varepsilon$, K, β and D represent the dynamic viscosity, thermal diffusivity, electrical permittivity, ionic mobility, thermal expansion coefficient and charge diffusion coefficient, respectively. The Rayleigh number Ra and electric Rayleigh number T are the two driving parameters of the system. They represent the ratio of the buoyancy and Coulomb force to the viscous force, respectively. The Prandtl number Pr is defined as the ratio of the momentum diffusivity to the thermal diffusivity. The injection parameter C measures the strength of the charge injection. The dimensionless mobility parameter M is defined as the ratio of the so-called hydrodynamic mobility and the true mobility of the ions; α is the dimensionless ion diffusion number and the typical values of α are between 10⁻⁴ and 10⁻³ (Pérez & Castellanos 1989).

In addition, to characterize the heat transfer in the flow, Nusselt numbers are defined

(2.9a,b)\begin{equation}N{u_i} = {R_i}\ \textrm{ln}\left( {\frac{{{R_o}}}{{{R_i}}}} \right){\left. {\frac{{\partial \theta }}{{\partial r}}} \right|_{r = {R_i}}},\quad N{u_o} = {R_o}\ \textrm{ln}\left( {\frac{{{R_o}}}{{{R_i}}}} \right){\left. {\frac{{\partial \theta }}{{\partial r}}} \right|_{r = {R_o}}},\end{equation}

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

(2.10a,b)\begin{equation}\overline {N{u_i}} = \frac{\textrm{1}}{{\textrm{2}{\rm \pi} }}\int_\textrm{0}^{\textrm{2}{\rm \pi} } {N{u_i}(\delta )\,\textrm{d}\delta } ,\quad \overline {N{u_o}} = \frac{\textrm{1}}{{\textrm{2}{\rm \pi} }}\int_\textrm{0}^{\textrm{2}{\rm \pi} } {N{u_o}(\delta )\,\textrm{d}\delta } ,\end{equation}

where δ is the angle (see figure 1). Theoretically, the values of $\overline {N{u_i}}$ and $\overline {N{u_o}}$ should be equal to each other due to energy conservation. However, a slight difference may be induced by the numerical error, so the final mean Nusselt number is defined as the average value of $\overline {N{u_i}}$ and $\overline {N{u_o}}$, i.e. $Nu = (\overline {N{u_i}} + \overline {N{u_o}} )/2$.

2.2. Linearization and global stability analysis

The linear problem is formulated following the work of He & Zhang (2021) based on Reynolds’ decomposition $\boldsymbol{N} = \bar{\boldsymbol{N}} + \boldsymbol{n^{\prime}}$, where N represents any flow variable discussed above, $\bar{\boldsymbol{N}}$ is the base state and $\boldsymbol{n^{\prime}}$ is the perturbation. After substituting the decompositions into the governing equations, subtracting from them the governing equations for the base states and retaining the terms of the first order, the linear system reads

(2.11)\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u^{\prime}} = \textrm{0,}\end{gather}

(2.12)\begin{gather}\frac{{\partial \boldsymbol{u^{\prime}}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u^{\prime}\bar{U}} + \boldsymbol{\bar{U}u^{\prime}}) =- \boldsymbol{\nabla }p^{\prime} + {\mathrm{\nabla }^2}\boldsymbol{u^{\prime}} + \frac{{Ra}}{{Pr}}\theta ^{\prime}{\boldsymbol{e}_y} + C{\left( {\frac{T}{M}} \right)^2}(q^{\prime}\bar{\boldsymbol{E}} + \bar{Q}\boldsymbol{e^{\prime}}),\end{gather}

(2.13)\begin{gather}\frac{{\partial \theta ^{\prime}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u^{\prime}}\bar{\varTheta } + \bar{\boldsymbol{U}}\theta ^{\prime}) = \frac{1}{{Pr}}{\mathrm{\nabla }^2}\theta ^{\prime},\end{gather}

(2.14)\begin{gather}{\mathrm{\nabla }^2}\phi ^{\prime} =- Cq^{\prime},\end{gather}

(2.15)\begin{gather}\boldsymbol{e^{\prime}} =- \boldsymbol{\nabla }\phi ^{\prime},\end{gather}

(2.16)\begin{gather}\frac{{\partial q^{\prime}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left[ {\left( {\frac{T}{{{M^2}}}\boldsymbol{e^{\prime}} + \boldsymbol{u^{\prime}}} \right)\bar{Q} + \left( {\frac{T}{{{M^2}}}\bar{\boldsymbol{E}} + \bar{\boldsymbol{U}}} \right)q^{\prime}} \right] = \frac{T}{{{M^2}}}\alpha {\mathrm{\nabla }^2}q^{\prime},\end{gather}

where $\bar{\boldsymbol{U}}$, $\bar{\varTheta }$, $\bar{\varPhi }$ and $\bar{Q}$ are the base states of velocity, temperature, electrical potential and charge density, respectively, and $\boldsymbol{u^{\prime}}$, $\theta ^{\prime}$, $\phi ^{\prime}$, $q^{\prime}$ are the corresponding perturbation fields. The boundary condition for the fluctuations on the two impermeable isothermal walls given by

(2.17)\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {\boldsymbol{u^{\prime}}{|_{r = {R_i}}} = 0,\quad \theta^{\prime}{|_{r = {R_i}}} = 0,\quad \phi^{\prime}{|_{r = {R_i}}} = 0,\quad q^{\prime}{|_{r = {R_i}}} = 0,}\\ {\boldsymbol{u^{\prime}}{|_{r = {R_o}}} = 0,\quad \theta^{\prime}{|_{r = {R_o}}} = 0,\quad \phi^{\prime}{|_{r = {R_o}}} = 0,\quad {{\left. {\dfrac{{\partial q^{\prime}}}{{\partial r}}} \right|}_{r = {R_o}}} = 0.} \end{array}} \right\}\end{equation}

In the linear stability analysis, it is a common practice to rewrite the fluid system in the form of a matrix; see the work of He & Zhang (2021). Therefore, (2.11)–(2.16) can be written in a compact manner as

(2.18)\begin{equation}\frac{{\partial \boldsymbol{n^{\prime}}}}{{\partial t}} = \boldsymbol{Mn^{\prime}},\end{equation}

where M represents the linearized Navier–Stokes operator for the ETC in the annulus and $\boldsymbol{n^{\prime}} = {({v^\prime },\eta ^{\prime},\phi ^{\prime},\theta ^{\prime})^T}$, $v^{\prime}$ and $\eta ^{\prime}$ are the wall-normal velocity and wall-normal vorticity, respectively. In the temporal stability analysis, a wave-like solution (Schmid, Henningson & Jankowski 2002) $\boldsymbol{n^{\prime}} = \hat{n}(x,y)\,{\textrm{e}^{ - \textrm{i}\omega t}}$ is assumed, leading to an eigenvalue problem that reads as

(2.19)\begin{equation}- \textrm{i}\omega \hat{n} = \boldsymbol{M}\hat{n}, \end{equation}

where −iω is the eigenvalue and $\hat{n}$ is the corresponding eigenvector. The complex-valued ω is the circular frequency of the perturbation, with its real part ω_r representing the phase speed and its imaginary part ω_i representing the growth rate of the linear perturbation. If ω_i = 0 when ω_r vanishes, it is said that an exchange of stabilities occurs, and the flow makes a transition to another steady state. Such a transition can be classified into pitchfork bifurcations, saddle-node bifurcations and transcritical bifurcations according to bifurcation theory (Crawford 1991). If ω_r ≠ 0 when ω_i vanishes, the steady solution is unstable to an oscillatory disturbance (this behaviour was named ‘overstability’ by Eddington 1920), and it makes a transition to a periodic solution with the angular velocity ω_r. Such a transition is called a Hopf bifurcation.

To solve the eigenvalue problem, the matrix-free method (Theofilis 2011) is engaged by using time-stepping methods, which employ temporal integrations of the perturbations using the linear Navier–Stokes equations, ${\partial _t}\boldsymbol{n^{\prime}} = \boldsymbol{Mn^{\prime}}$, from 0 to t. This is equivalent to applying the exponential operator $\textrm{exp}(\int_0^t \boldsymbol{M} \,\textrm{d}t)$ to the perturbation $\boldsymbol{n^{\prime}}$. Thus, the eigenvalue problem to be solved using the time-stepping method is defined by the exponential operator of the linearized version of the equations. The leading eigenvalue of the exponential operator to which the power method converges is usually the most unstable eigenvalue of the linearized operator.

2.3 Linearized lattice Boltzmann methods

To discover the full bifurcation diagram of the annular ETC, we need to simultaneously solve the fully coupled nonlinear equations (2.2)–(2.7) and linearized equations (2.11)–(2.16). Instead of directly solving the governing equations at the macroscopic level, here, we perform both the direct numerical simulation (DNS) and global linear stability analysis in the framework of the mesoscopic lattice Boltzmann method (LBM). As the LBM for the DNS of this problem has been well studied (Luo et al. 2016a,b; Lu, Liu & Wang 2019), we focus on the derivation of a LLBM solver for the linearized governing equations. In our LLBM model, a standard Bhatnagar–Gross–Krook scheme (Guo, Shi & Wang 2000) was employed for all fields, including the perturbation flow, perturbation electric potential, perturbation charge density and perturbation temperature.

For the perturbation flow and perturbation temperature, the LLBM equations have been built for the global linear stability analysis of incompressible flows (Pérez et al. 2017) and thermal convective flows (Jiang et al. 2022), under D2Q9 (two-dimensional nine-velocity) velocity discretization, and are formulated as follows:

(2.20)\begin{gather}{f_i}(\boldsymbol{x} + {\boldsymbol{e}_i}\varDelta t,t + \varDelta t) - {f_i}(\boldsymbol{x},t) =- \frac{1}{{{\tau _v}}}[\,{f_i}(\boldsymbol{x},t) - f_i^{eq}(\boldsymbol{x},t)] + \varDelta t{F_i}(\boldsymbol{x},t),\end{gather}

(2.21)\begin{gather}{l_i}(\boldsymbol{x} + {\boldsymbol{e}_i}\varDelta t,t + \varDelta t) - {l_i}(\boldsymbol{x},t) =- \frac{1}{{{\tau _\theta }}}[{l_i}(\boldsymbol{x},t) - l_i^{eq}(\boldsymbol{x},t)],\end{gather}

where f and l are the distribution functions of the perturbation fluid density and perturbation temperature, respectively. The relaxation times ${\tau _v}$ and ${\tau _\theta }$ are defined as

(2.22a,b)\begin{equation}{\tau _v} = \frac{{3\nu }}{{{c^2}\varDelta t}} + \frac{1}{2},\quad {\tau _\theta } = \frac{{3\chi }}{{{c^2}\mathrm{\Delta }t}} + \frac{1}{2}.\end{equation}

The corresponding equilibrium distributions $f_i^{eq}$ and $l_i^{eq}$, and the source terms F_i in (2.20) and (2.21) are formulated (Pérez et al. 2017; Jiang et al. 2022) as

(2.23a)\begin{gather}f_i^{eq} = \frac{{\rho ^{\prime}}}{{\bar{\rho }}}\bar{f}_i^{eq}\textrm{ + }\bar{\rho }{w_i}\left[ {\frac{{{\boldsymbol{e}_i} \boldsymbol{\cdot} \boldsymbol{u^{\prime}}}}{{c_s^2}} + \frac{{({\boldsymbol{e}_i} \boldsymbol{\cdot} \boldsymbol{u^{\prime}})({\boldsymbol{e}_i}\boldsymbol{\cdot }\bar{\boldsymbol{U}})}}{{c_s^4}} - \frac{{\boldsymbol{u^{\prime}} \boldsymbol{\cdot} \bar{\boldsymbol{U}}}}{{c_s^2}}} \right],\end{gather}

(2.23b)\begin{gather}\bar{f}_i^{eq} = \bar{\rho }{w_i}\left[ {1 + \frac{{{\boldsymbol{e}_i}\boldsymbol{\cdot }\bar{\boldsymbol{U}}}}{{c_s^2}} + \frac{{{{({\boldsymbol{e}_i}\boldsymbol{\cdot }\bar{\boldsymbol{U}})}^2}}}{{c_s^4}} - \frac{{{{\bar{\boldsymbol{U}}}^2}}}{{c_s^2}}} \right],\end{gather}

(2.24a)\begin{gather}l_i^{eq} = \frac{{\theta ^{\prime}}}{{\bar{\varTheta }}}\bar{l}_i^{eq} + \bar{\varTheta }{w_i}\left[ {\frac{{{\boldsymbol{e}_i}\boldsymbol{\cdot }\boldsymbol{u^{\prime}}}}{{c_s^2}} + \frac{{({\boldsymbol{e}_i}\boldsymbol{\cdot }\boldsymbol{u^{\prime}})({\boldsymbol{e}_i}\boldsymbol{\cdot }\bar{\boldsymbol{U}})}}{{c_s^4}} - \frac{{\boldsymbol{u^{\prime}}\boldsymbol{\cdot }\bar{\boldsymbol{U}}}}{{c_s^2}}} \right],\end{gather}

(2.24b)\begin{gather}\bar{l}_i^{eq} = \bar{\varTheta }{w_i}\left[ {1 + \frac{{{\boldsymbol{e}_i}\boldsymbol{\cdot }\bar{\boldsymbol{U}}}}{{c_s^2}} + \frac{{{{({\boldsymbol{e}_i}\boldsymbol{\cdot }\bar{\boldsymbol{U}})}^2}}}{{c_s^4}} - \frac{{{{\bar{\boldsymbol{U}}}^2}}}{{c_s^2}}} \right],\end{gather}

(2.25)\begin{gather}{F_i} = {w_i}\left( {1 - \frac{1}{{2{\tau_v}}}} \right)\frac{{{\boldsymbol{e}_i}\{ \boldsymbol{g}\beta [\bar{\rho }\theta ^{\prime} + \rho ^{\prime}(\bar{\varTheta } - {{\bar{\varTheta }}_{ref}})] + q^{\prime}\bar{\boldsymbol{E}} + \bar{Q}\boldsymbol{e^{\prime}}\} }}{{c_s^2}},\end{gather}

and the macroscopic quantities are calculated through (2.26), as follows:

(2.26) \begin{gather}\rho ^{\prime} = \sum\limits_i {{f_i}} ,\;\;\bar{\rho }\boldsymbol{u^{\prime}} + \rho ^{\prime}\bar{U} = \sum\limits_i {{\boldsymbol{e}_i}{\,f_i} + \frac{{\Delta t}}{2}\{ \boldsymbol{g}\beta [\bar{\rho }\theta ^{\prime} + \rho ^{\prime}(\bar{\varTheta } - {{\bar{\varTheta }}_{ref}})] + q^{\prime}\bar{\boldsymbol{E}} + \bar{Q}\boldsymbol{e^{\prime}}\} } ,\notag\\ \theta ^{\prime} = \sum\limits_i {{l_i}.}\end{gather}

The electrical part equations, including the Poisson equation for perturbation electric potential $\phi ^{\prime}$ and Nernst–Plank equation for the perturbation charge density $q^{\prime}$, are also modelled by LLBMs. Here, the improved lattice Boltzmann model for the Poisson equation proposed by Chai & Shi (2008) is used to solve $\phi ^{\prime}$, and the corresponding lattice Boltzmann equation (LBE) can be given as

(2.27)\begin{equation}{g_i}(\boldsymbol{x} + {\boldsymbol{e}_i}\varDelta t,t + \varDelta t) - {g_i}(\boldsymbol{x},t) =- \frac{1}{{{\tau _\phi }}}[{g_i}(\boldsymbol{x},t) - g_i^{eq}(\boldsymbol{x},t)] + \varDelta t{\varpi _i}R{D_a},\end{equation}

where g, e, ${\tau _\phi }$, R and D_a are distribution functions, microscopic velocity, relaxation time, source term and artificial diffusion coefficient, respectively. The D2Q5 (two-dimensional five-velocity) velocity discretization is adopted for the electric potential, with the equilibrium distribution function $g_i^{eq}$ and the corresponding weight coefficients ω and ϖ being expressed as

(2.28a–c) \begin{gather}g_i^{eq}(\boldsymbol{x},t) = \left\{ {\begin{array}{@{}ll@{}} {({\omega_0} - 1.0)\phi^{\prime}}&{i = 0}\\ {{\omega_i}\phi^{\prime}}&{i = 1 - 4} \end{array}} \right.,\quad {\omega _i} = \left\{ \begin{array}{@{}ll@{}} 0&{i = 0}\\ {1/4}&{i = 1 - 4} \end{array}\right.,\notag\\ {\varpi_i} = \left\{ {\begin{array}{@{}ll@{}} 0&{i = 0}\\ {1/4}&{i = 1 - 4} \end{array}}. \right.\end{gather}

The relaxation time ${\tau _\phi }$ and source term R can be computed from

(2.29a,b)\begin{equation}{\tau _\phi } = \frac{{{D_a}}}{{\beta {c^2}\varDelta t}} + \frac{1}{2},\quad R = Cq.\end{equation}

The artificial diffusion coefficient D_a (D_a > 0) is chosen to be D_a = 1/2 to balance the evolution speed and numerical stability. The coefficient β has been derived to be 1/2 by the Chapman–Enskog expansion (Chai & Shi 2008). Finally, the electric potential and charge density can be evaluated as

(2.30)\begin{equation}\phi ^{\prime} = \frac{1}{{1 - {\omega _0}}}\sum\limits_{i = 1}^4 {{g_i}} .\end{equation}

The LLBM for perturbation charge density $q^{\prime}$ is more complex because (2.16) is inherently nonlinear as the electric field E is a function of q. It is noticed that (2.16) has the form of a convection–diffusion equation (CDE). In the adoption of LBM for solving the CDE (Shi & Guo 2009; Chai & Shi 2020; Zhao et al. 2020; Chen et al. 2021; Wang et al. 2021a; Chai, Shi & Zhan 2022), the nonlinear charge mobility term can be treated as a cross-diffusion term (Wang et al. 2021a) or a convective term (Chai et al. 2022); these two models show almost no difference in this problem. Here, we can build the LLBM for $q^{\prime}$ by treating the nonlinear charge mobility term inspired by the work of Wang et al. (2021a). Under D2Q9 (two-dimensional nine-velocity) velocity discretization, we have

(2.31)\begin{equation}{h_i}(\boldsymbol{x} + {\boldsymbol{e}_i}\varDelta t,t + \varDelta t) - {h_i}(\boldsymbol{x},t) =- \frac{1}{{{\tau _q}}}[{h_i}(\boldsymbol{x},t) - h_i^{eq}(\boldsymbol{x},t)] + \varDelta t{G_i} + \varDelta t{S_i},\end{equation}

where $h_i^{eq}$ is the equilibrium distribution function, G_i is the term to recover the cross-diffusion term and S_i is the correction term to eliminate the additional terms appearing in the recovered macroscopic equations (Wang et al. 2021a,b). We further have

(2.32a–c)\begin{gather}h_i^{eq} = \frac{{q^{\prime}}}{{\bar{Q}}}\bar{h}_i^{eq} + \bar{Q}{w_i}\frac{{{\boldsymbol{e}_i} \boldsymbol{\cdot} \boldsymbol{u^{\prime}}}}{{c_s^2}},\quad \bar{h}_i^{eq} = \bar{Q}{w_i}\left[ {1 + \frac{{{\boldsymbol{e}_i} \boldsymbol{\cdot} \bar{\boldsymbol{U}}}}{{c_s^2}}} \right],\quad {\tau _q} = \frac{{\alpha T\nu }}{{{M^2}c_s^2\mathrm{\Delta }t}} + \frac{1}{2},\end{gather}

(2.33a)\begin{gather}\boldsymbol{G} =- \left( {\frac{T}{{{M^2}}}} \right)\frac{{\bar{Q}\boldsymbol{\nabla }\phi ^{\prime} + q^{\prime}\boldsymbol{\nabla }\bar{\varPhi }}}{{{\tau _q}\varDelta t}},\quad {G_i} = \frac{{{w_i}{\boldsymbol{e}_i}}}{{c_s^2}} \boldsymbol{\cdot} \boldsymbol{G}\end{gather}

(2.33b)\begin{gather}\boldsymbol{S} = \left( {1 - \frac{1}{{2{\tau_q}}}} \right){\partial _t}(\boldsymbol{u^{\prime}}\bar{Q}\boldsymbol{+ \bar{U}}q^{\prime}),\quad {S_i} = \frac{{{w_i}{\boldsymbol{e}_i}}}{{c_s^2}}\boldsymbol{\cdot }\boldsymbol{S}, \end{gather}

and the charge density is evaluated from

(2.34)\begin{equation}q^{\prime} = \sum\limits_j {{h_j}} .\end{equation}

Theoretically, we need to recover the linearized governing equations by Chapman–Enskog (CE) expansion. Pérez et al. (2017) first recovered the linearized Navier–Stokes equation without an external force from the LBE by a CE analysis. Recently, we have extended Pérez's model for the global instability analysis of thermal convection (Jiang et al. 2022) with the CE expansion. Besides, the CE analysis of Poisson-type perturbation charge potential equation can be found in the work of Chai & Shi (2008). Finally, the CE analysis of perturbation charge density equation is provided in Appendix A.

2.4 Fractal dimension and Lyapunov exponent

In this paper, the fractal dimension and Lyapunov exponent are adopted to further test the chaos quantitatively. The fractal dimension d_c measures the degrees of freedom associated with hydrodynamics, given by

(2.35)\begin{equation}{d_c} = \mathop {\lim }\limits_{s \to 0} \frac{{\log ({C_p}(s))}}{{\log (s)}},\end{equation}

where correlation function C_p(s) represents the minimum number of hyper cubes (of size s) covering a given set of points in a p-dimensional space, divided by N ² (N being the number of samples), and is defined as

(2.36)\begin{equation}{C_p}(s) = \mathop {\lim }\limits_{N \to \infty } \frac{1}{{{N^2}}}\sum\limits_{\scriptstyle i,j = 1 \atop \scriptstyle i \ne j}^N {h(s - |{\boldsymbol{X}_i} - {\boldsymbol{X}_j}|)} ,\end{equation}

in which h is the Heaviside function, ${\boldsymbol{X}_i}$ and ${\boldsymbol{X}_j}$ are a pair of given data points. The straight part of the curves for C_p(s) as a function of s on log–log scale will tend to have a constant slope as p increases. In particular, the slope of the periodic flow is close to 1, the slope of the quasi-periodic flow is close to 2 and the slope of the chaos is a real number greater than 2 (Grassberger & Procaccia 1983).

The algorithm of Wolf et al. (1985) is used for estimating the maximum Lyapunov exponent from time traces of the temperature at sampling point P. Since the implementation of this algorithm is relatively complicated, only a simple description is given here. For more details of this algorithm, interested readers can refer to Wolf et al. (1985). Generally, for a given time series {x ₁, x ₂, …, x_k, …}, its reconstructed phase space (dimension N) can be described as

(2.37)\begin{equation}Y({t_i}) = [x({t_i}),x({t_i} + {\tau _L}), \ldots ,x({t_i} + ({n_{dim}} - 1){\tau _L})],\quad i = 1,2, \ldots ,N,\end{equation}

in which n_dim is the number of embedding dimensions, τ_L represents the delay time. We follow the evolution in time of the initial point Y(t ₀) of the trajectory and its nearest point Y ₀(t ₀) belonging to a different portion of the trajectory. Then, the distance between these two points changes continuously over time. At time t ₁, the distance between these two points changes from L(t ₀) = |Y(t ₀) − Y ₀(t ₀)| to L(t ₁) = |Y(t ₁) − Y ₀(t ₁)|, then another point Y ₁(t ₁) of the trajectory is selected to keep the distance |Y ₁(t ₁) − Y ₀(t ₁)| = L(t ₀). Repeat the above process until Y_E(t_E) reaches the end of the time series; at this time, the total number of iterations is E. If the value of E is large enough, then the maximum Lyapunov exponent can be obtained by

(2.38)\begin{equation}{\lambda _{max}} = \frac{1}{{{t_E} - {t_0}}}\sum\limits_{k = 0}^E {\frac{{L({t_k})}}{{L({t_0})}}} .\end{equation}

It should be noted that the corresponding dynamical system is chaotic when the maximum Lyapunov exponent λ_max is positive (Rempel & Chian 2007). The larger the positive exponent, the more chaotic the system is.

2.6. Numerical method

Both the DNS and the solving of linear Navier–Stokes equations of the ETC are performed using the computational fluid dynamics solver Palabos (Latt et al. 2021), which is well known for its efficient parallelization and second-order accuracy. In the current work, the time step Δt is determined by the Courant–Friedrichs–Lewy condition with the target Courant number being 0.005. The mesoscopic boundary conditions of non-equilibrium extrapolation (Guo, Zheng & Shi 2002a,b) are used.

3. Validation

The grid independence test and basic validation results of DNS are shown in Appendix B. Here, the accuracy of the solution of the linear equations (2.17)–(2.22) using the LLBM solver was tested by calculating the linear stability criteria. Figure 2(a1,a2) shows the spatial distribution of the base flow and the leading global mode for natural convection (without the electric effect) at Pr = 0.0733, Ra = 1608.8 and A = 1.25. It was detected that the change of sign of the exponent ω_i occurs approximately at the critical value Ra_c = 2472. Our numerical results are in fairly good agreement with the data of Serrano-Aguilera et al. (2021) (Ra_c = 2438 by a spectral method), with a relative error of 1.39%.

Figure 2. Base flows (first row) and their corresponding most unstable mode (second row). (a) Natural convection at Pr = 0.0733, Ra = 1608.8 and A = 1.25; (b) electrohydrostatic state at T = 120, k = 8 and A = 2; (c) electroconvection state at T = 120 and A = 2.

Then, we verified the linear critical T_c of pure electric convection for both the first bifurcation and the second bifurcation. Figure 3 shows the neutral curve of the first bifurcation for modes k = 6, 7, 8 and 9 as a function of A at C = 10, M = 10 and α = 10⁻³. At the same parameter condition, Wu et al. (2014) obtained the linear criterion as T_{c 1} = 122 for mode k = 8 with A = 2 and α = 0. We find that the current T_{c 1} = 119.8 at A = 2 is lower than the linear stability criterion. Figure 2(b1) shows the electrohydrostatic state at T = 120 and figure 2(b2) gives the leading mode for k = 8, which is similar to the results of Pérez et al. (2014). After the emergence of EC, due to the competition between the ionic velocity and fluid velocity (Castellanos, Atten & Perez 1987), the flow is characterized by eight pairs of charge-void regions (where the value of q is almost zero in figure 2c1). As T continues to increase, the second bifurcation occurs. There has been no linear instability analysis of the second instability, and only DNSs were used to calculate the criterion by Huang et al. (2021), who obtained T_{c 2} = 213 for eight cells under parameters A = 2, M = 5, C = 10 and α = 0. We obtain the linear results as T_{c 2} = 242 with α = 10⁻³ (the eigenfunction q of the most unstable mode is shown in figure 2c2).

Figure 3. Linear critical values T_{c 1} of modes k = 6, 7, 8 and 9 as a function of A/(A + 2) at α = 10⁻³.

We find that the current T_{c 1} and T_{c 2} are lower and higher than the linear stability criterion obtained from the finite volume method, which neglected the charge diffusion effect. Here, we attribute this discrepancy to the effects of charge diffusion. Zhang et al. (2015, 2016) analysed the linear and weakly nonlinear stability of EC of dielectric liquids subjected to strong unipolar injection with the charge diffusion effect considered. It was found that the charge diffusion can destabilize the linearized EC but stabilizes the flow in an early phase of the nonlinear development of the disturbance, which means that weaker charge diffusion leads to greater T_{c 1} and smaller T_{c 2}, consistent with the comparison here. The quantitative effects of α on T_c for the first and second bifurcations are shown in figure 4. This dual effect of charge diffusion has also been discussed by Castellanos, Pérez & Atten (1989) in their DNS of such flows.

Figure 4. Effect of α on the criteria T_c of (a) the first bifurcation at A = 2 and M = 10; (b) the second bifurcation at A = 2 and M = 5.

4. Results and discussion

4.1. DNS

These parameters are fixed for all following cases: A = 1.25, C = 10, Pr = 10, M = 5 and α = 10⁻³. To study the effect of the radial electric field on the pure thermal convection in the concentric annulus, the first problem to be addressed is the determination of steady-state solutions in a Ra–T map, as shown in figure 5. To obtain this map, we visited the different positions of Ra ∈ [0, 3 × 10⁴] and T ∈ [0, 220] by dividing the map into 30 and 22 points, respectively. Therefore, this study aimed to provide a steady solution for 660 cases in the range of Ra and T investigated. Here, Ra is fixed and we increase T from zero; we solve for the steady-state solution by means of the Newton–Raphson iteration using as the initial guess the final-time result obtained from the previous simulations at lower T. As shown in figure 5, the solution region is divided into three flow states (seven subregions). The number of thermal plumes ‘P×’, the number of charge-void regions ‘C×’ and the stable/unstable flow state ‘S/U’ constitute the name of every subregion, i.e. the P1-C0-S solution represents steady convection with a main thermal plume at the top of the annulus without a charge-void region. Chaos is represented by the yellow dot region.

Figure 5. Map of flow patterns for ETC in a two-dimensional concentric annulus.

From the dimensionless momentum equation (2.2), the total destabilizing effect from electric and thermal forces along the radial direction is

(4.1)\begin{equation}{F_r} = |{C{T^2}\,Pr\,q\boldsymbol{E}+{M^2}Ra\theta \,\textrm{cos}(\delta )} |.\end{equation}

Compared with the plate–plate configuration, the electrothermo-convective phenomenon in the cylindrical configuration is much more complex. This is partially due to the differences in the directions of the driving forces and the basic flow patterns. From (4.1), the Rayleigh–Bénard thermal instability can be enhanced by destabilizing the electric field occurring in the top part of the thermally unstable region (red shadow region with δ < 45°); hydrodynamic instability can occur in the vertical section (yellow shadow regions with δ near 90°), which is called ‘natural convection heated from sidewalls’; and the instability caused by the baroclinic torque at the bottom of the annulus (red shadow region), as shown in figure 1. Therefore, instability of ETC in the annulus may be a combination of these mechanisms.

4.1.1. Bifurcations and multiple solutions with Ra

We keep T constant but different from zero (select three cases: T = 50, 100 and 150) and vary the Rayleigh number in such a way as to cross the critical line in a direction parallel to the Ra axis. Yoo (1999) first researched the saddle-node bifurcation for natural convection in an annulus when Pr > 0.5. The P1-C0-S solution tends to be excited more easily at zero-field initial conditions, and the P2-C1-S solution branch can be obtained with the help of artificial numerical disturbances, which were introduced during a short initial period (t < 0.01) to enhance the separation of the boundary layer from a point other than the top of the inner cylinder for one plume case. Saddle-node bifurcation occurs at Ra > Ra_s = 2024 for the case of T = 50, where the P1-C0-S solution branch and P2-C1-S solution branch are independent of each other, as shown in the bifurcation diagram in figure 6(a). As Ra continues to rise, the P2-C1-S solution transitions to an oscillatory state at Ra_c = 15 423; however, the P1-C0-S solution remains stable.

Figure 6. The mean Nu changes over Ra for three different T; (a) T = 50 (dotted lines represent the oscillation), (b) T = 100 with red coloured lines (increasing Ra), blue coloured lines (decreasing Ra) and T = 150 with green coloured lines (green points represent the chaos state).

Figure 6(b) shows the relation between Nu and Ra at T = 100 and T = 150. The saddle-node bifurcation disappears and is replaced by subcritical bifurcation at these two electric Rayleigh numbers. At T = 150, the P5-C5-S solution remains steady until Ra > Ra_c = 12 435, just like the case shown in figure 7(a). Then, infinite local maximum and minimum amplitudes of Nu indicate that steady convection transitions to the disorder state. Complex transition processes and hysteresis phenomena occur in the case of T = 100, and we summarize them as follows:

1. (i) The P5-C5-S solution possesses a symmetrical distribution along the vertical centreline at Ra < Ra_{c 1} = 5468, as shown in figure 7(a). The value of F_r decreases as Ra increases when δ > 90° from (4.1). At Ra = Ra_{c 1}, the strength of the electric plume does not maintain the existence of a charge-void region in the bottom of the annulus, and the P5-C5-S solution bifurcates to the P5-C4-S solution, as shown in figure 7(b).

2. (ii) The intensity of the main plume P_m (one plume with maximum flow intensity) continually increases with Ra and squeezes the room taken by secondary thermal plumes P_s (other plumes except P_m), so charge-void regions C_s (corresponding to P_s) move along the −y-direction for the P5-C4-S solution. Identically, these two charge-void regions will disappear once the magnitude of F_r is less than the critical finite amplitude at Ra > Ra_{c 2} = 7066. Figure 7(c) shows the distribution of the P3-C2-S solution after this transition.

3. (iii) Different from the former two processes, no plume suffers from the inhibition of the stabilizing buoyancy force for the P3-C2-S solution. Two large charge-void regions C_m (corresponding to P_s) are stable and fixed at their location with increasing Ra (<Ra_{c 3}). Figure 6(a) indicates that multiple plume structures are more unstable than a single plume. When Ra > Ra_{c 3} = 17 323, the P3-C2-S solution is broken and then convection transitions to the P1-C0-S state (in figure 7d) through a phase of chaos.

4. (iv) The hydrodynamic behaviours of this ETHD system with continually decreasing Ra can be understood as a reverse process of increasing Ra. The magnitude of F_r rises with the decline in Ra, and new charge-void regions occur at a special Ra_f at the bottom of the annulus. Therefore, the P5-C4-S and P3-C2-S solutions bifurcate to the P5-C5-S state at Ra_{f 1} = 2506 and Ra_{f 2} = 712, respectively. The P1-C0-S state first transitions to the P3-C2-S solution at Ra_{f 3} = 2820 and then returns to the P5-C5-S solution at Ra_{f 2}.

Figure 7. Four solutions at Ra = 5000 and T = 100. Panels (a–d) are the finite-amplitude solutions corresponding to the P5-C5-S, P5-C4-S, P3-C2-S and P1-C0-S solution branches, respectively.

4.1.2. Bifurcations and heat transfer enhancement with T

The next research question is investigation of the hydrodynamics of this ETC by varying T instead of Ra. We consider the case of Ra = 5000; to more clearly depict the subcritical feature, the bifurcation diagram is drawn in figure 8. Three hysteresis loops can be observed. The explanations are as follows:

1. (i) The P1-C0-S solution can be obtained from the zero-field initial condition at T = 0. Increasing T until T > T_{c 1} = 122, the instability exchange and the flow transition to the P5-C5-S solution branch (in figure 7a). An interesting transition occurs at T = 187, where the P5-C5-S solution loses its stability and all charge-void regions travel clockwise. We discuss this instability mechanism in the next section. However, the P5-C5-S state bifurcation to a P6-C6-S state is shown in the enlarged inset of figure 8 at T = 192. Then, the P6-C6-S solution transitions to the travelling wave state when T > T_{c 5} = 201.

2. (ii) Saddle-node bifurcation is proven in the ETHD system with a low electric Rayleigh number, as shown in figure 6(a), and we also obtain a P2-C1-S solution branch at T = 0. This P2-C1-S state is more stable than P1-C0-S, which bifurcates to the P5-C5-S solution at T > T_{c 2} = 129.

3. (iii) The P6-C6-S solution transitions to the P3-C2-S state (in figure 7c) with the disappearance of four charge-void regions at T < T_{f 1} = 99. These two charge-void regions cannot be maintained when T < T_{f 2} = 62 and the flow returns to the P1-C0-S state. Moreover, the P3-C2-S solution can be a saddle-node bifurcation solution branch, which bifurcates to the P5-C5-S solution at T_{c 3} = 111.

Figure 8. The mean Nu over T at Ra = 5000. The shaded area shows the enhanced heat transfer caused by the P6-C6-S solution compared with the P1-C0-S solution. The blue lines represent decreasing T continually.

In summary, this ETHD system consists of saddle-node bifurcations, subcritical bifurcations and supercritical bifurcations. Two conclusions can be drawn about the heat transfer from these determined bifurcation routes. First, the Nusselt number decreases with increasing Ra when thermal plumes exist in the lower part of the annulus. Second, the heat transfer efficiency is connected to the number of plumes under the same drive parameters, as shown in figures 6 and 8. Six plumes had the largest Nusselt number Nu = 2.90; however, the Nusselt number was only equal to 1.70 when only one plume existed, as shown in figure 7. Moreover, we can enhance the heat transfer in a determined manner: with the zero-field initial condition, we first increase the electric Nusselt number until T > T_{c 1} and then decrease T to the value that we need.

4.2. Global linear instability analysis

The transition of multiple plume solutions (more than two plumes) is explained by analysing the magnitude of F_r at the bottom of the annulus successfully. However, it cannot explain the transition behaviour of the single and two-plume solutions (P1-C0-S and P2-C1-S) on the saddle-node bifurcation branches shown in figures 6(a) and 8. In this section, we carry out a global linear instability analysis by analysing the spectrum of the eigenvalues of the operator M. Due to the symmetry of the boundary condition, the two-dimensional base flow admits the symmetry property about the x = 0 line, i.e. Φ_b(x, y, t) = Φ_b(−x, y, t) with Φ = [u, v, θ, q, ϕ]. Solutions φ of the linearized Navier–Stokes equations admit the following two types of symmetry with respect to the x = 0 plane:

(4.2)\begin{gather}\textrm{S-mode}\quad\varphi (x,y,t) = \varphi ( - x,y,t),\end{gather}

(4.3)\begin{gather}\textrm{A-mode}\quad \varphi (x,y,t) =- \varphi ( - x,y,t).\end{gather}

Figure 9(a,b) illustrates the behaviour of eigenvalues in figure 6(a) with T = 50 and Ra ∈ [14 500, 15 300]. Eigenpairs are extracted by adding a random symmetry/anti-symmetry perturbation about the x = 0 line as the initial condition. We only give the S-mode (leading mode) solution here because the A-mode has a too rapid damping rate with time to be extracted. The P1-C0-S solution reflects a rapid decline with monotonic instability in this parameter range. However, this instability of the P1-C0-S solution will be replaced with overstability because the flow is thermally dominant, where the magnitude of buoyancy is approximately 4.5 times larger than the electric force from (4.1). Therefore, instability of the P1-C0-S solution is similar to the case of pure thermal convection discussed by Serrano-Aguilera et al. (2021), who proved that the S-mode is the most unstable mode and the bifurcation is Hopf type, which is consistent with our conclusion. For the P2-C1-S solution, it shows an overstability with a critical angular frequency ω_r = 12.07 at Ra_c = 15 423 and the corresponding distribution of eigenfunctions θ and q are shown in figures 10(a1) and 10(a2). The eigenfunction θ has its maxima at an approximate angle of 45° in the lower region of the ascending jets, so these parts will be first destabilized and oscillate around the steady state. From these results, the A-mode is more stable than the S-mode and steady convection losses its stability through a Hopf bifurcation.

Figure 9. Variation in the growth rate ω_i and the frequency ω_r with Ra and T. (a,b) Varying Ra from 14 500 to 15 300 with T = 50; (c,d) varying T between 114 and 122 with Ra = 5000. The red line represents the P1-C0-S solution, the blue line is the P2-C1-S solution, the solid line is the S-mode and the dotted line is the A-mode.

Figure 10. (a) Critical mode of the P2-C1-S solution at Ra = 15 423 and T = 50; (b,c) are the (s)-mode and (a)-mode of the P1-C0-S solution at Ra = 5000 and T = 121.8, respectively.

For the case of saddle-node bifurcation solution branches in figure 8 with T ∈ [114, 122] and Ra = 5000, figures 9(c) and 9(d) give the eigenvalues of the A-mode and S-mode for the P1-C0-S and P2-C1-S solutions with T, respectively. We first focus on the eigenvalue of the P1-C0-S solution. Different from the conclusions above, the A-mode has a larger increasing rate with T than the S-mode and it becomes the leading unstable mode instead of the S-mode when T > 115. Moreover, the leading A-mode plays a role of monotonic instability, while the S-mode reflects the overstability. Instability exchange occurs at T_c = 122, and the distribution of eigenfunctions θ and q of these two modes are shown in figures 10(b) and 10(c). The same conclusion emerges for the P2-C1-S solution at this range of T. Therefore, the strength ratio of buoyancy and electric force tend to influence the symmetry of the leading mode and the bifurcation types.

4.3. Energy analysis

To gain a deeper understanding of instability mechanisms from the global linear instability results above, an energy analysis method was used. The energy analysis was first applied to the EHD flow by Zhang et al. (2015) in local analysis. Here, we use this method to explain the instability mechanism for the ETC. The hydrodynamic perturbation kinetic energy is derived by multiplying the linearized (2.18) by the complex conjugate of the perturbation velocity $u_i^\textrm{ + }$, i.e.

(4.4) \begin{align} & u_i^ + \frac{{\partial {u_i}}}{{\partial t}} + u_i^ + {u_j}\frac{{\partial {{\bar{U}}_i}}}{{\partial {x_j}}} + u_i^ + {\bar{U}_j}\frac{{\partial {u_i}}}{{\partial {x_j}}}\notag\\ &\quad =- u_i^ + \frac{{\partial p}}{{\partial {x_i}}} + u_i^\textrm{ + }\frac{{{\partial ^\textrm{2}}{u_i}}}{{\partial x_j^\textrm{2}}} + u_i^ + {\left( {\frac{T}{M}} \right)^\textrm{2}}\left( {\frac{{\partial \bar{\phi }}}{{\partial {x_i}}}\frac{{{\partial^\textrm{2}}\phi }}{{\partial x_j^\textrm{2}}} + \frac{{\partial \phi }}{{\partial {x_i}}}\frac{{{\partial^\textrm{2}}\bar{\phi }}}{{\partial x_j^\textrm{2}}}} \right) + u_i^\textrm{ + }\frac{{Ra}}{{Pr}}\theta \boldsymbol{\cdot }{\boldsymbol{e}_y},\end{align}

where u_i, ϕ, θ and p are perturbation values. By taking the complex conjugate of the obtained equation and averaging the two equations, the rate of the perturbation energy density is obtained. After discarding transport terms (which only redistribute the energy inside the domain and do not contribute to the net energy rate in the case of no-penetration boundary conditions) and integrating the equation in the domain, we obtain the following:

(4.5) \begin{align} \dfrac{{\partial K}}{{\partial t}} &= \underbrace{{ - \dfrac{1}{2}\int_\varOmega {(u_i^ + {u_j} + u_j^ + {u_i})\dfrac{{\partial {{\bar{U}}_i}}}{{\partial {x_j}}}} }}_{{{K_v}}}\underbrace{{ - \int_\varOmega {\dfrac{{\partial u_i^ + }}{{\partial {x_j}}}\dfrac{{\partial {u_i}}}{{\partial {x_j}}}} }}_{{{K_d}}}\underbrace{{ + \dfrac{1}{2}\int_\varOmega {\dfrac{{Ra}}{{Pr}}(\theta u_i^++ {\theta ^ + }{u_i})\boldsymbol{\cdot }{\boldsymbol{e}_y}} }}_{{{K_b}}}\notag\\ &\quad \underbrace{{\begin{array}{@{}c@{}} + \dfrac{\textrm{1}}{2}\displaystyle\int_\varOmega {{\left( {\dfrac{T}{M}} \right)}^2}\Bigg[ - \dfrac{{\partial \bar{\phi }}}{{\partial {x_i}}}\Bigg( \dfrac{{\partial {\phi^\textrm{ + }}}}{{\partial {x_j}}}\dfrac{{\partial {u_i}}}{{\partial {x_j}}} + \dfrac{{\partial \phi }}{{\partial {x_j}}} \dfrac{{\partial u_i^ + }}{{\partial {x_j}}} \Bigg) \\ - \dfrac{{{\partial^2}\bar{\phi }}}{{\partial {x_i}\partial {x_j}}}\left( {u_i^ + \dfrac{{\partial \phi }}{{\partial {x_j}}} + {u_i}\dfrac{{\partial {\phi^\textrm{ + }}}}{{\partial {x_j}}}} \right) - \dfrac{{{\partial^3}\bar{\phi }}}{{\partial {x_i}\partial {x_i}\partial {x_j}}}(\phi u_j^ ++ {\phi^\textrm{ + }}{u_j}) \vphantom{\dfrac{{\partial u_i^ + }}{{\partial {x_j}}}}\Bigg] \end{array}}}_{{{K_e}}}\notag\\ &\quad \underbrace{{\begin{array}{@{}c@{}}+ \dfrac{1}{2}\displaystyle\int_\varOmega \dfrac{\partial }{{\partial {x_j}}}\Bigg[ - ({{\bar{U}}_j}{u_i}u_i^ + ) - (\,pu_i^ ++ {p^ + }{u_i}){\delta_{ij}} + \Bigg( {u_i^ + \dfrac{{\partial {u_i}}}{{\partial {x_j}}} + {u_i}\dfrac{{\partial u_i^ + }}{{\partial {x_j}}}} \Bigg)\\ + {{\left( {\dfrac{T}{M}} \right)}^2}\dfrac{{\partial \bar{\phi }}}{{\partial {x_i}}}\left( {\dfrac{{\partial \phi }}{{\partial {x_j}}}u_i^ ++ \dfrac{{\partial {\phi^\textrm{ + }}}}{{\partial {x_j}}}{u_i}} \right) \Bigg] \end{array}}}_{{Transport\;Terms}}, \end{align}

where $K = u_i^ + {u_i}/2$ is the perturbation energy density. The first term on the right-hand side of (4.5) represents the interaction between the disturbance velocity and the base flow (K_v), the second term indicates the viscous dissipation of the disturbance energy (K_d), the third to fifth terms are the energy transfer terms between the velocity fluctuation field and the electric field (K_e) and the last term describes the energy exchanged between the fluctuating flow field and the thermal field, which is the rate of work done by buoyancy (K_b). Summing the right-hand side terms, a positive sign means that the terms play a destabilizing role on the flow, and vice versa.

Figure 11(a) shows the time evolution of various perturbed kinetic energies for the P2-C1-S solution at Ra_c = 15 423 and T = 50. Immediately, one can make several direct observations. First, viscous dissipation K_d is always negative for the hydrodynamics. The terms that can lead to growth in the hydrodynamic disturbance energy density are linked to the energy transfer from the electric field, K_e, K_b and K_v; K_b is the leading term for destabilizing, K_e is the second important term and production K_v can be neglected due to its weak influence. The most efficient electric mechanism seems to be related to the term K_{e 12}, which represents the interaction between the x-direction perturbed electric field and the vertical direction velocity under the constant effect of the vertical direction base flow field. Figure 11(c) gives the local contribution of the perturbed energy budget over the angle δ from 0° to 180°, which indicates that the destabilized electric field enhanced the Rayleigh–Bénard instability at the top of the annulus for δ < 60°. There are two main destabilized areas: δ ∈ [0°, 11.4°] with a peak value at δ = 0° and δ ∈ [32.4°, 46.8°] with a peak value at δ = 40.8°, where plume cycles originate and return to the inner cylinder, respectively. When δ > 70°, all energy terms tend to zero, so these regions are stable.

Figure 11. Temporal evolution and spatial contributions to the local growth rate of the perturbed kinetic energy through mechanisms of various budget terms. Panels (a,c) are the P2-C1-S solution at Ra = 15 423 and T = 50, respectively; panels (b,d) are the P1-C0-S solution at T = 120 and Ra = 5000, respectively.

Figure 12(a) shows the local spatial structures of various budget terms for this P2-C1-S solution. The viscous dissipation term K_d is proportional to the magnitude of the velocity gradient, so the boundary layer and the excessive region of the jets are the main perturbation energy dissipation (not given here). The energy budget representing buoyancy K_b is found around the jet region of the plume at the top of the annulus (in figure 12a2). The electric field both stabilizes and destabilizes this flow in different areas, as shown in figure 12(a3), where the main destabilizing region of K_e is opposite to K_b.

Figure 12. Contributions to the local growth rate of the perturbed kinetic energy through mechanisms of various budget terms: (a) P2-C1-S solution at Ra = 15 423, T = 50; (b) P1-C0-S solution at Ra = 5000 and T = 120.

Figure 11(b) shows the time series of the energy terms for the P1-C0-S solution at Ra = 5000 and T = 120 (<T_c), in which a transient energy rise is found at t ∈ [0, 0.05] by adding a random perturbation of q as the initial condition. Different from the overstability shown in figure 11(a), the pitchfork bifurcation predicts a steady exponential decay for both the electric term K_e and the buoyancy term K_b. Figure 11(d) gives the local contributions of the perturbed energy budget, and there are destabilizing areas δ ∈ [6.2°, 29.9°] with a peak value of δ = 22.2° and a stabilized region at δ < 6.2°. The local instability region in figure 12(b) illustrates that the total destabilizing energy dK/dt concentrates in the boundary layers around the outside circle, which mainly comes from the contribution of K_e. These destabilizing regions tend to excite a pair of charge-void regions that are symmetric along the vertical midplane. Therefore, the P1-C0-S solution bifurcates to the P3-C2-S solution and then continues to bifurcate the P5-C5-S state, as shown in the inset of figure 11(d).

Here, we sought a quantitative energy guideline. Figure 13(a) shows the time evolution of the two main destabilizing terms K_e and K_b for T = 95, 96, 97 and 98 at Ra = 5000 and t > 1.7, i.e. the model instability phase. The leading role of K_b is gradually replaced by K_e in this range, and K_e is the largest destabilizing term when T > T_energy (critical electric Rayleigh number with K_b = K_e); T_energy ∈ [95, 98] is less than the linear critical value T_c = 122. Figure 13(b) indicates the energy budget terms for different electric Rayleigh numbers. We obtained these terms by using the most unstable mode with a kinetic density k = 10 as the initial condition and then evolving it over time. The magnitude of all terms increases with T, and the most efficient terms seem to relate to K_{e 11}. The only stabilizing term K_d cannot offset the destabilizing effect of K_e and K_b, which causes the sign of dK/dt to change at T = 122. Therefore, the electric destabilizing mechanism may be associated with the A-mode, while thermal mechanisms correspond to the S-mode.

Figure 13. (a) Time evolution of K_b and K_e for four different T = 98, 97, 96 and 95 for the P1-C0-S solution with Ra = 5000; (b) is the perturbation kinetic energy contributions for the P1-C0-S solution at T ∈ [114, 122] and Ra = 5000.

4.4. Transition to oscillation and disorder

Energy analysis predicted the local instability region for the P2-C1-U solution at Ra_c and T = 50. Therefore, we first investigated the instability features of this solution after Hopf bifurcation at Ra = 18 000. Figure 14(a) shows the velocity time series and the corresponding Fourier frequency spectrum at sampling point P = (0.1, 0.2). The amplitude of the velocity fluctuation increases gradually and eventually maintains a stable periodic oscillation. The dominant frequency of the oscillatory flow (f ₁ = 2.03) intuitively illustrates this oscillation. In figure 15(a), we show snapshots of charge-void regions (green line with q = 0.05) and the temperature distribution at the special time points t ₁, t ₂ and t ₃ marked in figure 14(a2). The two-plume structure is still maintained but oscillates with time along the vertical direction, and there is one pair of local maximum (Nu = 2.46) and minimum (Nu = 2.37) amplitudes in the time series. Moreover, the P1-C0-U solution in figure 5 has the same instability behaviours.

Figure 14. Time evolution of u for three different unstable solutions (t ₁, t ₂ and t ₃) and the Fourier frequency spectrum of temperature at the sampling point. (a) One-frequency oscillation at Ra = 5000, T = 180; (b) one-frequency oscillation for the P2-C1-S solution at Ra = 18 000, T = 50.

Figure 15. Distributions of temperature fields and the charge-void region (green line with q = 0.05) at three different times t ₁, t ₂ and t ₃.

There is another instability type, the global travelling wave from the P5-C5-U or P6-C6-U solutions. Figures 14(b1) and 14(b2) display the time evolution of the velocity u and its Fourier spectrum for the P5-C5-U solution at T = 180 and Ra = 5000. The flow shows a regular mono-periodic oscillation characterized by a main frequency f ₁ = 6.02 and its harmonics. To better understand this unsteady flow, three instantaneous snapshots of the charge-void regions and temperature distribution corresponding to times t ₁, t ₂ and t ₃ (labelled in figure 14b2) are shown in figure 15(b). It is clear that charge-void regions and thermal plumes rotate clockwise. Huang et al. (2021) found that all charge-void regions were located on the clockwise side and then swung to the anticlockwise side in one period for the second instability of the pure electroconvection. The baroclinic effect caused by the stabilizing thermal effect at the bottom of the annulus may be responsible for the difference between the ETC and EC. Moreover, part of the Coulomb force is used to drive the fluid's circumferential motion, which reduces the radial convective intensity. Therefore, the mean Nusselt is lower than the steady state (i.e. Nu = 4.2 at T = 175) and slightly oscillatory about Nu at approximately 3.61.

Finally, the routes of steady flow transition to chaos were investigated. In addition to the basic analysis methods above, i.e. the bifurcation diagram analysis and the frequency spectral analysis, the fractal dimension d_c and Lyapunov exponent λ_max were adopted to quantitatively test the chaos based on our previous work (Li et al. 2020). The transition route from steady to chaos through periodic and quasi-periodic phases for ETC has been investigated by Li et al. (2020). Therefore, we focused on another route that occurs widely in our present ETHD system, such as that in figures 6(b) and 5 with Ra > 24 000.

Figure 16 gives the time evolution of u at Ra = 12 500 and T = 150; to reduce the computational time, this series used the P5-C5-S solution at Ra = 12 400 as the initial condition. Through a steady increasing phase occurring at t < 7, irregular oscillation was observed for a large time t >10. A continuous broadband frequency spectrum (not shown here) with a main frequency f ₁ = 2.73 indicates a typical chaos. Before computing the Lyapunov exponents, the nonlinear equations (2.1) are first integrated for a long time (t > 50) to decay all initial transitions and let the flow enter the chaotic state without such a disturbance. The values of λ _max and d_c are listed in table 1. The maximum Lyapunov exponents and fractal dimension are zero at Ra = 12 400. At Ra = 12 430, d_c > 1 and λ _max > 0 indicate a steady transition to chaos directly without suffering a periodic and quasiperiodic phase. Around this transition criterion Ra_c of chaos, the variation in Ra is not significant, and these Lyapunov dimensions are found to be smaller than 4. The distributions of the temperature fields and the charge-void regions with Ra = 12 500 at t ₁, t ₂ and t ₃ are shown in figure 16, where regular flow patterns are broken and the Nusselt number is slightly larger than 3.

Figure 16. Time evolution of u at Ra = 12 500 and T = 150 for three different unstable solutions (t ₁, t ₂ and t ₃). Chaos from the P5-C5-S solution.

Table 1. Summary of the flow state, fractal dimension and maximum Lyapunov exponent for different electric Rayleigh numbers T and Rayleigh number Ra.

5. Conclusion

Global linear instability and bifurcations in ETC of a dielectric liquid confined between a two-dimensional concentric annulus subjected to a strong unipolar injection were investigated numerically. A LLBM was first proposed by Pérez et al. (2017) and extended to solve the whole set of coupled linear equations, including the linear Navier–Stokes equations, the linear energy equation, Poisson's equation, and the linear charge conservation equation. A multiscale analysis was also performed to recover the macroscopic linearized Navier–Stokes equations from four different discrete LBEs. It was validated by calculating the linear critical value of two-dimensional natural convection and the LLBM had an error of 1.39% compared with the spectral method.

Seven kinds of solutions exist in this ETC system due to the complex bifurcations, including saddle-node, subcritical and Hopf bifurcations. These bifurcation routes constitute at most four solution branches. Force magnitude analysis, global linear instability analysis and asymptotic energy analysis were used to explain the instability mechanism and transition of different solutions and predict the local instability regions. Instability with global travelling wave behaviour is a unique behaviour in the annulus configuration ETHD system, which may be caused by the baroclinity. Finally, the chaotic behaviour was quantitatively analysed through the calculation of the fractal dimension and Lyapunov exponent.