## 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.* Reference Bart, Tavrow, Mehregany and Lang1990; Darabi *et al.* Reference Darabi, Rada, Ohadi and Lawler2002), to design devices and drug delivery systems (Chakraborty *et al.* Reference Chakraborty, Liao, Adler and Leong2009) and DNA microarrays (Lee *et al.* Reference Lee, Cho, Huh, Ko, Lee, Jang, Lee, Kang and Choi2006) and to design new strategies for active flow control (Bushnell & Mcginley Reference Bushnell and Mcginley1989). 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 Reference Laohalertdecha, Naphon and Wongwises2007), such as in boiling and condensing systems (Cotton *et al.* Reference Cotton, Robinson, Shoukri and Chang2005), drying and evaporating systems (Rashidi *et al.* Reference Rashidi, Bafekr, Masoodi and Languri2017) and solar energy systems (Ghalamchi *et al.* Reference Ghalamchi, Kasaeian, Ghalamchi, Fadaei and Daneshazarian2017). An excellent review covering various aspects of electro-thermo-convection (ETC) can be found in Léal *et al.* (Reference Léal, Miscevic, Lavieille, Amokrane, Pigache, Topin, Nogarède and Tadrist2013). 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.* Reference Togun, Ab Dulrazzaq, Kazi, adarudin, Kadhum and Sadeghinezhad2014). In the literature, there are few numerical works (Fernandes *et al.* Reference Fernandes, Lee, Alapati and Yong2012; Wu *et al.* Reference Wu, Traoré, Zhang, Pérez and Vázquez2016; Hassen *et al.* Reference Hassen, Oztop, Kolsi, Borjini and Abu-Hamdeh2017; Li *et al.* Reference Li, He, Luo and Yi2019; He *et al.* Reference He, Chai, Wang, Ma and Shi2021; Ma *et al.* Reference Ma, Zhang, Wang, He and Li2022) 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 (Reference Worraker and Richardson1979) 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 (Reference Castellanos, Atten and Velarde1984) 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 (Reference Rodriguez-Luis, Castellanos and Richardson1986) investigated the onset of steady convection under charge injection strengths of any finite magnitude. Pontiga & Castellanos (Reference Pontiga and Castellanos1994) 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 (Reference Roberts1969) initially performed a linear instability study on the ETC by assuming the variation in dielectric constant as a function of temperature. Pontiga, Castellanos & Richardson (Reference Pontiga, Castellanos and Richardson1992) 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 (Reference Taraut and Smorodin2012) found that the flow system evolves from the monotonic mode to the oscillatory mode via a transition in frequency. Mordvinov & Smorodin (Reference Mordvinov and Smorodin2012) 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.* (Reference Guan, He, Wang, Song, Zhang and Wu2021) 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 Reference Félici1971; Lacroix, Atten & Hopfinger Reference Lacroix, Atten and Hopfinger1975; Zhang Reference Zhang2016) 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 (Reference Richardson1980), Agrait & Castellanos (Reference Agrait and Castellanos1990), Fernandes *et al.* (Reference Fernandes, Lee, Park and Suh2013) and recently Wu *et al.* (Reference Wu, Vázquez, Traoré and Pérez2014) 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.*(Reference Huang, Wang, Guan, Du, Deepak Selvakumar and Wu2021) 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 (Reference Yoo1999) 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 (Reference Mizushima, Hayashi and Adachi2001) 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 (Reference Serrano-Aguilera, Blanco-Rodríguez and Parras2021) conducted abundant numerical investigations on dual solutions for values of the Prandtl number between 0.01 and 1 and Rayleigh numbers between 10

^{2}and 5 × 10

^{6}. 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.*(Reference Traoré, Pérez, Koulova and Romat2010) 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 (Reference Atten, Lacroix and Malraison1980) 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 (Reference Malraison and Atten1982) 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.* Reference Malraison, Atten, Berge and Dubois1983), and the fractal dimension of the chaotic attractor was calculated based on the Grassberger–Procaccia method (Grassberger & Procaccia Reference Grassberger and Procaccia1983). 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.* (Reference Atten, Caputo, Malraison and Gagne1984) 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 (Reference Chicón, Pérez and Castellanos2001) applied the algorithm proposed by Kantz & Schreiber (Reference Kantz and Schreiber2004) 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.* (Reference Wolf, Swift, Swinney and Vastano1985) has been widely used recently (Feng *et al.* Reference Feng, Zhang, Vazquez and Shu2021; Huang *et al.* Reference Huang, Wang, Guan, Du, Deepak Selvakumar and Wu2021). Moreover, Li *et al.* (Reference Li, Su, Luo and Yi2020) 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 (Reference Pérez, Aguilar and Theofilis2017). 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*is located in the centre, and the aspect ratio is fixed:

_{i}*A*= 1.25. The inner and outer cylinders are kept at different constant temperatures

*θ*and

_{hot}*θ*(Δ

_{cold}*θ*

_{0}=

*θ*−

_{hot}*θ*), respectively, and a radial direct current electric field is established by applying a voltage difference Δ

_{cold}*ϕ*

_{0}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.*Reference Laohalertdecha, Naphon and Wongwises2007), which means that the injected charge density

*q*

_{0}at the inner cylinder remains constant and uniform.

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 Reference Castellanos1998). We introduce dimensionless quantities, denoted with a star, as follows:

where *L* = *R _{o}* −

*R*is the length scale. If we drop the superscript stars for clarity, the resulting governing equations of this ETC system are (Traoré

_{i}*et al.*Reference Traoré, Pérez, Koulova and Romat2010)

In the above equations, ** u** = (

*u*,

_{x}*u*) is the velocity,

_{y}**= (**

*E**E*,

_{x}*E*) is the electric field and

_{y}

*e**is the vertical direction. The scalars*

_{y}*ρ*,

*θ*,

*ϕ*and

*q*denote the fluid density, temperature, electric potential and charge density, respectively. Six dimensionless parameters appear in governing equations, defined as

*a*–

*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^{−4} and 10^{−3} (Pérez & Castellanos Reference Pérez and Castellanos1989).

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

*a*,

*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*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

_{o}*a*,

*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 (Reference He and Zhang2021) 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

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

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 (Reference He and Zhang2021). Therefore, (2.11)–(2.16) can be written in a compact manner as

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 Reference Schmid, Henningson and Jankowski2002) $\boldsymbol{n^{\prime}} = \hat{n}(x,y)\,{\textrm{e}^{ - \textrm{i}\omega t}}$ is assumed, leading to an eigenvalue problem that reads as

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

*ω*representing the growth rate of the linear perturbation. If

_{i}*ω*= 0 when

_{i}*ω*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 Reference Crawford1991). If

_{r}*ω*≠ 0 when

_{r}*ω*vanishes, the steady solution is unstable to an oscillatory disturbance (this behaviour was named ‘overstability’ by Eddington Reference Eddington1920), and it makes a transition to a periodic solution with the angular velocity

_{i}*ω*. Such a transition is called a Hopf bifurcation.

_{r} To solve the eigenvalue problem, the matrix-free method (Theofilis Reference Theofilis2011) 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.* Reference Luo, Wu, Yi and Tan2016*a*,Reference Luo, Wu, Yi and Tan*b*; Lu, Liu & Wang Reference Lu, Liu and Wang2019), 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 Reference Guo, Shi and Wang2000) 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.* Reference Pérez, Aguilar and Theofilis2017) and thermal convective flows (Jiang *et al.* Reference Jiang, Luo, Zhang, Wu and Yi2022), under D2Q9 (two-dimensional nine-velocity) velocity discretization, and are formulated as follows:

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

*a*,

*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.*Reference Pérez, Aguilar and Theofilis2017; Jiang

*et al.*Reference Jiang, Luo, Zhang, Wu and Yi2022) as

*a*)\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}

*b*)\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}

*a*)\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}

*b*)\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}

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

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 (Reference Chai and Shi2008) is used to solve $\phi ^{\prime}$, and the corresponding lattice Boltzmann equation (LBE) can be given as

where *g*, ** e**, ${\tau _\phi }$,

*R*and

*D*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

_{a}*ω*and

*ϖ*being expressed as

*a*–

*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

*a*,

*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*> 0) is chosen to be

_{a}*D*= 1/2 to balance the evolution speed and numerical stability. The coefficient

_{a}*β*has been derived to be 1/2 by the Chapman–Enskog expansion (Chai & Shi Reference Chai and Shi2008). Finally, the electric potential and charge density can be evaluated as

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 Reference Shi and Guo2009; Chai & Shi Reference Chai and Shi2020; Zhao

*et al.*Reference Zhao, Wu, Chai and Shi2020; Chen

*et al.*Reference Chen, Chai, Shang and Shi2021; Wang

*et al.*Reference Wang, Wei, Li, Chai and Shi2021

*a*; Chai, Shi & Zhan Reference Chai, Shi and Zhan2022), the nonlinear charge mobility term can be treated as a cross-diffusion term (Wang

*et al.*Reference Wang, Wei, Li, Chai and Shi2021

*a*) or a convective term (Chai

*et al.*Reference Chai, Shi and Zhan2022); 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.*(Reference Wang, Wei, Li, Chai and Shi2021

*a*). Under D2Q9 (two-dimensional nine-velocity) velocity discretization, we have

where $h_i^{eq}$ is the equilibrium distribution function, *G _{i}* is the term to recover the cross-diffusion term and

*S*is the correction term to eliminate the additional terms appearing in the recovered macroscopic equations (Wang

_{i}*et al.*Reference Wang, Wei, Li, Chai and Shi2021

*a*,Reference Wang, Yang, Wang, Chai and Wei

*b*). We further have

*a*–

*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}

*a*)\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}

*b*)\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

Theoretically, we need to recover the linearized governing equations by Chapman–Enskog (CE) expansion. Pérez *et al.* (Reference Pérez, Aguilar and Theofilis2017) 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.* Reference Jiang, Luo, Zhang, Wu and Yi2022) with the CE expansion. Besides, the CE analysis of Poisson-type perturbation charge potential equation can be found in the work of Chai & Shi (Reference Chai and Shi2008). 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

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*

^{2}(

*N*being the number of samples), and is defined as

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 Reference Grassberger and Procaccia1983).

The algorithm of Wolf *et al.* (Reference Wolf, Swift, Swinney and Vastano1985) 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.* (Reference Wolf, Swift, Swinney and Vastano1985). Generally, for a given time series {*x* _{1}, *x* _{2}, …, *x _{k}*, …}, its reconstructed phase space (dimension

*N*) can be described as

in which *n _{dim}* is the number of embedding dimensions,

*τ*represents the delay time. We follow the evolution in time of the initial point

_{L}*Y*(

*t*

_{0}) of the trajectory and its nearest point

*Y*

_{0}(

*t*

_{0}) belonging to a different portion of the trajectory. Then, the distance between these two points changes continuously over time. At time

*t*

_{1}, the distance between these two points changes from

*L*(

*t*

_{0}) = |

*Y*(

*t*

_{0}) −

*Y*

_{0}(

*t*

_{0})| to

*L*(

*t*

_{1}) = |

*Y*(

*t*

_{1}) −

*Y*

_{0}(

*t*

_{1})|, then another point

*Y*

_{1}(

*t*

_{1}) of the trajectory is selected to keep the distance |

*Y*

_{1}(

*t*

_{1}) −

*Y*

_{0}(

*t*

_{1})| =

*L*(

*t*

_{0}). Repeat the above process until

*Y*(

_{E}*t*) reaches the end of the time series; at this time, the total number of iterations is

_{E}*E*. If the value of

*E*is large enough, then the maximum Lyapunov exponent can be obtained by

It should be noted that the corresponding dynamical system is chaotic when the maximum Lyapunov exponent *λ _{max}* is positive (Rempel & Chian Reference Rempel and Chian2007). 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.* Reference Latt2021), 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 Reference Guo, Zheng and Shi2002*a*,Reference Guo, Zheng and Shi*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(*a*1,*a*2) 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*= 2472. Our numerical results are in fairly good agreement with the data of Serrano-Aguilera

_{c}*et al.*(Reference Serrano-Aguilera, Blanco-Rodríguez and Parras2021) (

*Ra*= 2438 by a spectral method), with a relative error of 1.39%.

_{c} 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

^{−3}. At the same parameter condition, Wu

*et al.*(Reference Wu, Vázquez, Traoré and Pérez2014) 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(

*b*1) shows the electrohydrostatic state at

*T*= 120 and figure 2(

*b*2) gives the leading mode for

*k*= 8, which is similar to the results of Pérez

*et al.*(Reference Pérez, Vázquez, Wu and Traoré2014). After the emergence of EC, due to the competition between the ionic velocity and fluid velocity (Castellanos, Atten & Perez Reference Castellanos, Atten and Perez1987), the flow is characterized by eight pairs of charge-void regions (where the value of

*q*is almost zero in figure 2

*c*1). 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.*(Reference Huang, Wang, Guan, Du, Deepak Selvakumar and Wu2021), 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

^{−3}(the eigenfunction

*q*of the most unstable mode is shown in figure 2

*c*2).

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.*(Reference Zhang, Martinelli, Wu, Schmid and Quadrio2015, 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*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 (Reference Castellanos, Pérez and Atten1989) in their DNS of such flows.

_{c}## 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^{−3}. 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^{4}] 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.

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

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 (Reference Yoo1999) 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*= 15 423; however, the P1-C0-S solution remains stable.

_{c} 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:

(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*decreases as_{r}*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*).(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*is less than the critical finite amplitude at_{r}*Ra*>*Ra*_{c}_{2}= 7066. Figure 7(*c*) shows the distribution of the P3-C2-S solution after this transition.(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 7*d*) through a phase of chaos.(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*rises with the decline in_{r}*Ra,*and new charge-void regions occur at a special*Ra*at the bottom of the annulus. Therefore, the P5-C4-S and P3-C2-S solutions bifurcate to the P5-C5-S state at_{f}*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}.

#### 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:

(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 7*a*). 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.(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.(iii) The P6-C6-S solution transitions to the P3-C2-S state (in figure 7

*c*) 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.

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

**. Due to the symmetry of the boundary condition, the two-dimensional base flow admits the symmetry property about the**

*M**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:

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.* (Reference Serrano-Aguilera, Blanco-Rodríguez and Parras2021), 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*= 15 423 and the corresponding distribution of eigenfunctions

_{c}*θ*and

*q*are shown in figures 10(

*a*1) and 10(

*a*2). 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.

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.* (Reference Zhang, Martinelli, Wu, Schmid and Quadrio2015) 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.

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:

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*), the third to fifth terms are the energy transfer terms between the velocity fluctuation field and the electric field (

_{d}*K*) 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 (

_{e}*K*). Summing the right-hand side terms, a positive sign means that the terms play a destabilizing role on the flow, and

_{b}*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*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,

_{d}*K*,

_{e}*K*and

_{b}*K*;

_{v}*K*is the leading term for destabilizing,

_{b}*K*is the second important term and production

_{e}*K*can be neglected due to its weak influence. The most efficient electric mechanism seems to be related to the term

_{v}*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 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*is found around the jet region of the plume at the top of the annulus (in figure 12

_{b}*a*2). The electric field both stabilizes and destabilizes this flow in different areas, as shown in figure 12(

*a*3), where the main destabilizing region of

*K*is opposite to

_{e}*K*.

_{b} 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*and the buoyancy term

_{e}*K*. Figure 11(

_{b}*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 d

*K/*d

*t*concentrates in the boundary layers around the outside circle, which mainly comes from the contribution of

*K*. 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(

_{e}*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*for

_{b}*T*= 95, 96, 97 and 98 at

*Ra*= 5000 and

*t*> 1.7, i.e. the model instability phase. The leading role of

*K*is gradually replaced by

_{b}*K*in this range, and

_{e}*K*is the largest destabilizing term when

_{e}*T*>

*T*(critical electric Rayleigh number with

_{energy}*K*=

_{b}*K*);

_{e}*T*∈ [95, 98] is less than the linear critical value

_{energy}*T*= 122. Figure 13(

_{c}*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*cannot offset the destabilizing effect of

_{d}*K*and

_{e}*K*, which causes the sign of d

_{b}*K/*d

*t*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.

### 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*

_{1}= 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*

_{1},

*t*

_{2}and

*t*

_{3}marked in figure 14(

*a*2). 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.

There is another instability type, the global travelling wave from the P5-C5-U or P6-C6-U solutions. Figures 14(*b*1) and 14(*b*2) 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* _{1} = 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* _{1}, *t* _{2} and *t* _{3} (labelled in figure 14*b*2) are shown in figure 15(*b*). It is clear that charge-void regions and thermal plumes rotate clockwise. Huang *et al.* (Reference Huang, Wang, Guan, Du, Deepak Selvakumar and Wu2021) 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

*λ*were adopted to quantitatively test the chaos based on our previous work (Li

_{max}*et al.*Reference Li, Su, Luo and Yi2020). The transition route from steady to chaos through periodic and quasi-periodic phases for ETC has been investigated by Li

*et al.*(Reference Li, Su, Luo and Yi2020). 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* _{1} = 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*> 1 and

_{c}*λ*

_{max}> 0 indicate a steady transition to chaos directly without suffering a periodic and quasiperiodic phase. Around this transition criterion

*Ra*of chaos, the variation in

_{c}*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*

_{1},

*t*

_{2}and

*t*

_{3}are shown in figure 16, where regular flow patterns are broken and the Nusselt number is slightly larger than 3.

## 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.* (Reference Pérez, Aguilar and Theofilis2017) 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.

## Acknowledgments

The authors are grateful to Dr W. Lei for his help with the derivation of equations in this paper.

## Funding

This work is supported by the National Natural Science Foundation of China (grant nos 52076055, 52276057). The financial support from the Ministry of Education, Singapore, is also acknowledged (the WBS no. R-265-000-689-114).

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. The CE expansion for LLBM

Here, we provided the CE analysis of perturbation charge density equation. In (2.31), *h _{i}* and $h_i^{eq}$ satisfy

*a*–

*c*)\begin{equation}\sum\limits_i {{h_i}} = \sum\limits_i {h_i^{eq}} = q^{\prime},\quad \sum\limits_i {{\boldsymbol{e}_i}h_i^{eq}} = \bar{Q}\boldsymbol{u^{\prime}}\textrm{ + }q^{\prime}\bar{\boldsymbol{U}}\textrm{,}\quad \sum\limits_i {{\boldsymbol{e}_i}{\boldsymbol{e}_i}h_i^{eq}} = q^{\prime}c_s^2\boldsymbol{I}.\end{equation}

The corresponding cross-diffusion term *G _{i}* and the correction term

*S*satisfy

_{i}*a*–

*d*)\begin{equation}\sum\limits_i {{G_i}} = 0,\quad \sum\limits_i {{\boldsymbol{e}_i}{G_i}} = \boldsymbol{G},\quad \sum\limits_i {{S_i}} = 0,\quad \sum\limits_i {{\boldsymbol{e}_i}{S_i}} = \boldsymbol{S}.\end{equation}

With an expansion of the distribution function *h _{i}* and terms

*G*and

*S*in time and space, we obtain

*a*)\begin{gather}{h_i} = h_i^{eq} + \varepsilon h_i^{(1)} + {\varepsilon ^2}h_i^{(2)},\end{gather}

*b*)\begin{gather}{G_i} = \varepsilon G_i^{(1)},\quad {S_i} = \varepsilon S_i^{(1)}.\end{gather}

Combining (A1) and (A3), we have

*a*)\begin{gather}\sum\limits_i {h_i^{(n)}} = \textrm{0}\quad \textrm{(}n \ge \textrm{1),}\end{gather}

*b*)\begin{gather}\sum\limits_i {G_i^{(1)}} = 0\textrm{,}\quad \sum\limits_i {{\boldsymbol{e}_i}G_i^{(1)}} = {\boldsymbol{G}^{(1)}},\quad \sum\limits_i {S_i^{(1)}} = 0\textrm{,}\quad \sum\limits_i {{\boldsymbol{e}_i}S_i^{(1)}} = {\boldsymbol{S}^{(1)}}.\end{gather}

By applying the Taylor expansion to (2.31), we obtain

where ${D_i} = {\partial _t} + {\boldsymbol{e}_i}\boldsymbol{\cdot }\boldsymbol{\nabla } = \varepsilon {\partial _{{t_1}}} + {\varepsilon ^2}{\partial _{{t_2}}} + \varepsilon {\boldsymbol{e}_i}\boldsymbol{\cdot }{\boldsymbol{\nabla }_1}$. Setting ${D_{i1}} = {\partial _{t1}} + {\boldsymbol{e}_i}\boldsymbol{\cdot }{\boldsymbol{\nabla }_1}$, then we have ${D_i} = \varepsilon {D_{i1}} + {\varepsilon ^2}{\partial _{{t_2}}}$. Inserting (A3) into (A5), and by the scale analysis, we can obtain the following equations for *ε* ^{1} and *ε* ^{2}:

*a*)\begin{gather}{\varepsilon ^1}\quad{D_{i1}}h_i^{eq} =- \frac{1}{{{\tau _q}\varDelta t}}h_i^{(1)} + G_i^{(1)} + S_i^{(1)},\end{gather}

*b*)\begin{gather}{\varepsilon ^2}\quad{\partial _{{t_2}}}h_i^{eq} + {D_{i1}}h_i^{(1)} + \frac{{\varDelta t}}{2}D_{i1}^2h_i^{eq} =- \frac{1}{{{\tau _q}\varDelta t}}h_i^{(2)}.\end{gather}

Inserting (A6*a*) into (A6*b*), we can obtain

Summing (A6*a*) and (A7) over *i*, we can obtain

*a*)\begin{gather}\frac{{\partial q^{\prime}}}{{\partial {t_1}}} + {\boldsymbol{\nabla }_1}\boldsymbol{\cdot }{\bf (}\bar{Q}\boldsymbol{u^{\prime}}\textrm{ + }q^{\prime}\bar{\boldsymbol{U}}) = 0,\end{gather}

*b*)\begin{gather}\frac{{\partial q^{\prime}}}{{\partial {t_\textrm{2}}}} + \sum\limits_i {{D_{i1}}\left[ {\left( {1 - \frac{1}{{2{\tau_q}}}} \right)h_i^{(1)}} \right]} + \frac{{\varDelta t}}{2}\sum\limits_i {{\boldsymbol{e}_i}\boldsymbol{\cdot }{\boldsymbol{\nabla }_1}[G_i^{(1)} + S_i^{(1)}]} = 0.\end{gather}

The third term in (A8*b*) can be expressed as

Consider the second term in (A8*b*),

Using (A6*a*), we have

Substituting (A11) from (A8*b*), we have

Summing (A12) multiplied by *ε* ^{2} and (A8*a*) multiplied by *ε*, we can obtain

and setting $(T/{M^2})\alpha = c_s^2({\tau _q} - 1/2)\varDelta t$, we can recover (A14*e*).

Finally, the LLBM equations can recover the linearized macroscopic equations (A14*a*–*e*) by CE expansion

*a*)\begin{gather}\frac{{\partial \rho ^{\prime}}}{{\partial t}}+\boldsymbol{\nabla }\boldsymbol{\cdot }(\bar{\rho }\boldsymbol{u^{\prime}}+\rho ^{\prime}\bar{\boldsymbol{U}}) = \textrm{0,}\end{gather}

*b*)\begin{gather}\frac{{\partial \bar{\rho }\boldsymbol{u^{\prime}}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\bar{\rho }\boldsymbol{u^{\prime}\bar{U}} + \bar{\rho }\boldsymbol{\bar{U}u^{\prime}}) =- \boldsymbol{\nabla }p^{\prime} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\mu \boldsymbol{\nabla }\boldsymbol{u^{\prime}}) + {\bar{\rho }_0}\boldsymbol{g}\beta (\theta - {\theta _\textrm{0}}) + q\boldsymbol{E},\end{gather}

*c*)\begin{gather}\frac{{\partial \theta ^{\prime}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u^{\prime}}\bar{\varTheta } + \bar{\boldsymbol{U}}\theta ^{\prime}) = \boldsymbol{\nabla }\boldsymbol{\cdot }(\chi \boldsymbol{\nabla }\theta ^{\prime}),\end{gather}

*d*)\begin{gather}{\mathrm{\nabla }^2}\phi ^{\prime} =- Cq^{\prime},\quad \boldsymbol{e^{\prime}} =- \boldsymbol{\nabla }\phi ^{\prime},\end{gather}

*e*)\begin{gather}\frac{{\partial q^{\prime}}}{{\partial t}} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\bar{Q}\boldsymbol{u^{\prime}} + q^{\prime}\bar{\boldsymbol{U}}) = \frac{T}{{{M^2}}}\boldsymbol{\nabla }\boldsymbol{\cdot }(\alpha {\mathrm{\nabla }^2}q^{\prime} - \boldsymbol{e^{\prime}}\bar{Q} - \bar{\boldsymbol{E}}q^{\prime}).\end{gather}

It is noted that (A14*b*) is obtained by subtracting the moment equation of the base state and eliminating the high-order perturbation term (Pérez *et al.* Reference Pérez, Aguilar and Theofilis2017). As the LLBM is an artificial compressible solver for incompressible ETC, under incompressible conditions, $\bar{\rho } = \rho ^{\prime} = \textrm{constant}$ and constant parameters, after non-dimensionalization, we can obtained non-dimensional incompressible governing equations (2.11)–(2.16).

## Appendix B. Validations for DNSs

In this section, we first present the result of the grid independency test. The mean Nusselt number is computed at *T* = 100 and *Ra* = 1 × 10^{5} under different grid resolutions, as shown in table 2. The relative errors between resolutions of 400 × 400 and 500 × 500 are less than 0.2%. Therefore, a resolution of 400 × 400 is applied in the following numerical simulations considering both the accuracy and the computational efficiency.

Then, we verify the thermal module and electric module. The benchmark case for natural convection at *Pr* = 0.706 and *Ra* = 4.7 × 10^{4} is considered. The radial temperature distributions for three different angles *δ* = 0°, 90° and 180° and *Nu* at the inner and outer surfaces are plotted in figures 17(*a*) and 17(*b*). Our results have a rather small error compared with the result obtained by the spectral code (Serrano-Aguilera *et al.* Reference Serrano-Aguilera, Blanco-Rodríguez and Parras2021). Finally, the hydrostatic solution is selected to validate the solving of Poisson’s equation for the electric potential and charge conservation by LBM. There are analytical solutions of hydrostatic state available for the concentric cylinder, which may be expressed as

*a*–

*c*)\begin{equation}{\boldsymbol{u}_s} = [0,0],\quad {q_s}(R) = {A_e}{[\delta ({R^2} + {B_e})]^{ - 1/2}},\quad {E_s}(R) = \frac{{\delta {A_e}}}{R}{[\delta ({R^2} + {B_e})]^{1/2}},\end{equation}

where *δ* = 1 for the inner injections; *A _{e}* and

*B*are two constants depending on $\varGamma = {R_0}/{R_1}$. In the case of strong injection

_{e}*C*= 10, the values for

*A*and

_{e}*B*can be found in Agrait & Castellanos (Reference Agrait and Castellanos1990). In figure 18(

_{e}*a*,

*b*), we validate our numerical steady solutions by DNS against the analytical solutions equation (A1). The difference is 0.8% for

*q*and 0.183% for

*E*in terms of the root-mean-square error for

_{r}*R*= 0.1, 0.3 and 0.5.