Electroconvection of Thin Liquid Crystals: Model Reduction and Numerical Simulations

TL;DR

This paper develops a finite element method to simulate electroconvection in thin liquid crystals between electrodes, comparing infinite and slim electrode configurations, and analyzing parameter effects on convection sustainability.

Contribution

It introduces a novel numerical approach for electroconvection modeling in thin liquid crystals with different electrode geometries and studies parameter influences on the phenomena.

Findings

Slim electrodes promote sustained electroconvection.

Finite element method effectively simulates electroconvection.

Parameter variations significantly affect convection behavior.

Abstract

We propose a finite element method for the numerical simulation of electroconvection of thin liquid crystals. The liquid is located in between two concentric circular electrodes which are either assumed to be of infinite height or slim. Each configuration results in a different nonlocal electro-magnetic model defined on a two dimensional bounded domain. The numerical method consists in approximating the surface charge density, the liquid velocity and pressure, and the electric potential in the two dimensional liquid region. Finite elements for the space discretization coupled with standard time stepping methods are put forward. Unlike for the infinite electrodes configuration, our numerical simulations indicate that slim electrodes are favorable for electroconvection to occur and are able to sustain the phenomena over long period of time. Furthermore, we provide a numerical study on the influence of the three main parameters of the system: the Rayleigh number, the Prandtl number and the electrodes aspect ratio.

Figures

Tables

Table 1. Table 1 : Infinite electrodes configuration: Number of quadrilaterals # ​ 𝒯 # 𝒯 \#\mathcal{T} and subsequent uniform refinements used in the coarse subdivision for different values of aspect ratio α 𝛼 \alpha .

$\alpha$	0.1	0.2	0.33	0.452	0.56	0.6446	0.9
$\#𝒯$	4	5	6	9	12	15	29
$\#$ uniform ref.	5	5	5	5	5	5	6

Table 2. Table 2 : Parameter settings for infinite electrode simulations. Here “DoFs” stands for degrees of freedom.

$\alpha$	$\tau$	DoFs for ${\Psi }_{0,\Omega }$	DoFs for $q$	DoFs for $𝒖$ and $p$
0.1	0.001	4,224	4,224	33,280/4,224
0.2	0.001	5,280	5,280	41,600/5,280
0.33	0.001	7,392	7,392	58,240/7,392
0.452	0.001	9,504	9,504	74,880/9,504
0.56	0.001	12,672	12,672	99,840/12,672
0.6446	0.001	15,840	15,840	124,800/15,840
0.8	0.001	120,640	120,640	957,696/120,640

Table 3. Table 3 : Parameter settings for slim electrode simulations. Here “DoFs” stands for degrees of freedom.

$\alpha$	$\tau$	$M$	DoFs for ${\Psi }_{0,\Omega }$	DoFs for $q$	DoFs for $𝒖$ and $p$
$\le 0.452$	0.001	3	9,009	6,272	49,664/6,272
0.56, 0.6446	0.001	3	24,960	24,960	198,144/24,960
0.8	0.001	3	99,072	99,072	789,504/99,072

Equations

$$\nabla \times \mathbf{E}=0,$$

$$\mathbf{E}=-\nabla \Psi .$$

$$\nabla \cdot \mathbf{E}=0\text{or}\Delta \Psi =0\text{in}{D}_s.$$

$$\Psi =V\text{on}{\mathcal{K}}_i\text{and}\Psi =0\text{on}{\mathcal{K}}_o.$$

$$\underset{z\to \infty }{\lim }|\Psi (x,y,z)|=0,(x,y)\in {\mathbb{R}}^2,$$

$$\underset{z\to \infty }{\lim }|\Psi (x,y,z)|<\infty ,(x,y)\in \Omega ,$$

$$G(\bar{x},\bar{y},r):=(\bar{x}-r,\bar{x}+r)\times (\bar{y}-r,\bar{y}+r)\times (-s-r,s+r)$$

$${\epsilon }_0(\mathbf{E}(\bar{x},\bar{y},s)-\mathbf{E}(\bar{x},\bar{y},-s))\cdot {\mathbf{e}}_3=2q+O(s).$$

$$\rho ({\partial }_t+\mathbf{u}\cdot \nabla )\mathbf{u}+\eta \Delta \mathbf{u}+\nabla p=\mathbf{L},\nabla \cdot \mathbf{u}=0,\text{both in }{\Omega }_s\times {\mathbb{R}}_{+},$$

$$\underset{z\downarrow 0}{\lim }\frac{\partial }{\partial z}\Psi -\underset{z\uparrow 0}{\lim }\frac{\partial }{\partial z}\Psi =-\frac{2q}{{\epsilon }_0}\text{in }\Omega ,$$

$$\Delta \Psi =0,in{D}^{+},$$

$$\begin{array}{rl}& \Psi =V\text{leavevmode}\text{on}\text{leavevmode}{\mathcal{K}}_i^{+},\text{leavevmode}\Psi =0\text{leavevmode}\text{on}\text{leavevmode}{\mathcal{K}}_o^{+},\text{leavevmode}\text{and }\underset{z\downarrow 0}{\lim }{\partial }_z\Psi =-\frac{q}{{\epsilon }_0}\text{leavevmode}\text{on}\text{leavevmode}\Omega ,\end{array}$$

$$\begin{array}{rl}{\rho }_{\Omega }({\partial }_t+{\mathbf{u}}_{\Omega }\cdot {\nabla }_{\Omega }){\mathbf{u}}_{\Omega }+{\eta }_{\Omega }{\Delta }_{\Omega }{\mathbf{u}}_{\Omega }+{\nabla }_{\Omega }{p}_{\Omega }& ={\mathbf{L}}_{\Omega },\text{in }\Omega \times {\mathbb{R}}_{+},\\ {\nabla }_{\Omega }\cdot {\mathbf{u}}_{\Omega }& =0,\text{in }\Omega \times {\mathbb{R}}_{+},\end{array}$$

$${\mathbf{u}}_{\Omega }=0\text{on}\partial \Omega \times {\mathbb{R}}_{+}.$$

$${\nabla }_{\Omega }q\cdot \mathbf{\nu }=-{\sigma }_{\Omega }{\mathbf{E}}_{\Omega }\cdot \mathbf{\nu },$$

$${\partial }_{t}q+{\nabla }_{\Omega }\cdot ({\mathbf{u}}_{\Omega }q-{\sigma }_{\Omega }{\mathbf{E}}_{\Omega })=0,\text{in }\Omega \times {\mathbb{R}}_{+}.$$

$${\partial }_{t}q+{\nabla }_{\Omega }\cdot ({\mathbf{u}}_{\Omega }q)-{\sigma }_{\Omega }{\Delta }_{\Omega }{\Psi }_{\Omega }=0,\text{in }\Omega \times {\mathbb{R}}_{+},$$

$${\int }_{\Omega }q(t)={\int }_{\Omega }q(0),t>0.$$

$$\hat{\mathbf{x}}:=\frac{\mathbf{x}}{d},\hat{t}:=\frac{{\sigma }_{\Omega }}{{\epsilon }_{0}d}t,\hat{\Omega }:=\{\hat{\mathbf{x}}:\mathbf{x}\in \Omega \},\text{and}{\hat{D}}^{+}:=\{\hat{\mathbf{x}}:\mathbf{x}\in D^{+}\},$$

$$\hat{\Psi }(\hat{\mathbf{x}},\hat{t})$$

$$\hat{\mathbf{u}}(\hat{\mathbf{x}},\hat{t})$$

$${\partial }_{t}q+\mathbf{u}\cdot {\nabla }_{\Omega }q-{\Delta }_{\Omega }{\Psi }_{\Omega }=0,in\Omega \times {\mathbb{R}}_{+},$$

$$-\Delta \Psi =0,in{D}^{+}\times {\mathbb{R}}_{+},$$

$$\begin{array}{rl}& \Psi =1\text{leavevmode}\text{on}\text{leavevmode}{\mathcal{K}}_i^{+},\text{leavevmode}\Psi =0\text{leavevmode}\text{on}\text{leavevmode}{\mathcal{K}}_o^{+},\text{leavevmode}\text{and}\text{leavevmode}\underset{z\downarrow 0}{\lim }{\partial }_z\Psi =-q\text{leavevmode}\text{on}\text{leavevmode}\Omega \end{array}$$

$${\partial }_t\mathbf{u}+(\mathbf{u}\cdot {\nabla }_{\Omega })\mathbf{u}-\mathcal{P}{\Delta }_{\Omega }\mathbf{u}+{\nabla }_{\Omega }p=-\mathcal{R}\mathcal{P}q{\nabla }_{\Omega }{\Psi }_{\Omega },{\nabla }_{\Omega }\cdot \mathbf{u}=0,$$

$$\mathcal{P}:=\frac{{ϵ}_0{\eta }_{\Omega }}{{\rho }_{\Omega }{\sigma }_{\Omega }d},and\mathcal{R}:=\frac{{ϵ}_0^2{V}^2}{{\eta }_{\Omega }{\sigma }_{\Omega }}.$$

$$R_i=\frac{\alpha }{1-\alpha },\text{and}{R}_o=\frac{1}{1-\alpha }.$$

$$\Delta {\Psi }_K=0\; in{D}^{+}$$

$${\Psi }_K(x,y,z)=\eta (x,y):=\frac{\ln (\sqrt{x^2+y^2}/(1-\alpha ))}{\ln (\alpha )}.$$

$$-\Delta {\Psi }_0=0\; in{D}^{+}$$

$${\Psi }_0=0\text{leavevmode}\text{on}\text{leavevmode}{\mathcal{K}}_i^{+}\cup {\mathcal{K}}_o^{+},\text{and}\underset{z\downarrow 0}{\lim }{\partial }_z{\Psi }_0=-q\text{on }\Omega .$$

Full text

1 Introduction

The electrical convection, or electroconvection in short, is the appearance of pairs of vortices caused by an electric field across a thin-layer of electrically charged fluid. We refer to  Morris et al. [27], Daya et al. [12], Deyirmenjian et al. [16, 17], Daya et al. [13, 14, 15], Tsai et al. [29, 31] for experimental studies in different geometries and to Kim et al. [25], Kwak et al. [26] for applications of electroconvection in bio-technologies.

In this work, we follow the experimental setting in Daya et al. [13] and Tsai et al. [30]. The thin film of liquid crystal is at room temperature and the crystals are arranged in layers with their long axis aligned with the normal of the film plane (semectic-A phase). Under such arrangement, the fluid molecules cannot move across the film thickness, thereby restricted to a planar motion. The fluid has a poor conductivity in the direction perpendicular to the long axis of the liquid crystal molecules. As a consequence, the magnetic fields generated by the low current inside the fluid are neglected. The thin-layered film is placed in between two concentric electrodes as illustrated in Figure 1. Two different settings are considered: the electrodes extend to infinity in the fluid normal direction (infinite case) and the electrodes have negligible thickness (slim case).

The mathematical models are derived in the limiting case of vanishing thickness. The resulting surface charge density conservation relation and the mass and momentum conservation relations for the incompressible liquid are derived on a two dimensional bounded domain $\Omega$. However, the electric field (or potential) remains described in all $\mathbb{R}^{3}$. To fully take advantage of a reduced modeling setting, we resort to equivalent nonlocal representations of the restriction to $\Omega$ of the electric potential. Depending on the electrodes configuration, these nonlocal representations involve either the spectral Laplacian defined via spectral expansions (infinite electrodes) or the integral Laplacian defined using a Fourier transform (slim electrodes). We refer to Bonito et al. [8] for a survey of numerical methods associated with these two operators.

Numerical simulations of electroconvection for slim electrodes are obtained in Tsai et al. [30], Tsai [32], Jabbour [24]. Second order backward differentiation schemes for the time evolution coupled with spectral methods for the space discretization of all the quantities are advocated. The space discretizations are based on standard continuous, piecewise polynomial finite elements: Taylor-Hood elements for the fluid dynamics and continuous piecewise (bi-)linear finite elements for the electric potential and the surface charge density. At this point, it is worth mentioning that the discretization of the nonlocal operators only requires finite element solvers for standard reaction-diffusion elliptic partial differential equations, thereby making the overall algorithm easily implementable. The evolution of the velocity-pressure in the liquid is approximated using a second order backward difference method. In our experimental settings and unlike for the fluid dynamics, the surface charge density conservation law may be convection dominated and prone to numerical instabilities. To circumvent this issue, we propose an explicit second order strongly stability preserving scheme (SSP) (see Gottlieb and Gottlieb [19]) coupled with a second order entropy viscosity stabilization from Guermond and Popov [22] for its approximation.

We assess numerically the advantages of the two different electrodes configurations. We also provide a numerical study on the effect of three critical nondimensional parameters in the electroconvection process. They are the Rayleigh number $\mathcal{R}$ representing the ratio between the electric forcing and viscous dissipation, the Prandtl number $\mathcal{P}$ measuring the ratio between the charge relaxation time and the viscous relaxation time, and the circular electrodes aspect ratio $\alpha$.

The outline of this paper is as follows. We begin with the derivation of the full mathematical model in Section 2 followed with its two dimensional model reduction in Section 3. The numerical algorithms for the approximations of the fluid dynamics, the surface charge density and electric potential are described in Section 4. Their general performances together with numerical assessments of favorable conditions to obtain electroconvection are discussed in Section 5. When using the same experimental configurations as in Tsai et al. [30], Tsai [32], our findings are in good agreement. This is also detailed in Section 5. Conclusions are drawn in Section 6.

2 Mathematical models

We describe below the derivation of the electroconvection model from Maxwell system. We consider two settings: infinite and slim electrodes. The former corresponds to the simulation reported in Daya et al. [13, 14], Tsai et al. [30] while the later is related to the analysis in Constantin et al. [11]. Worth mentioning, vanishing charge densities on $\partial\Omega$ are considered in Constantin et al. [11] allowing for smoother charge densities and simplifications of the mathematical model. We do not make such assumption incompatible with the requirement of charge conservation and refer to Remark 3.1 for further details.

We recall that our experimental setting consists of a thin layer of liquid crystal in an annular region. We assume that the liquid crystal molecules can only move along the $xy-$plane and polarized in the $z$ direction without being able to change polarization. This corresponds to the so-called Smectic-A phase, which manifests itself when the liquid crystal are at room temperature $24\pm 2^{\circ}$C. At this point, we remark that the assumptions made on the state of the liquid crystals are critical for the subsequently derived model reduction. For instance, if instead in a Smectic-C phase, the long axis of liquid crystal molecules pertain a fixed angle with $z$ direction. The possibility of model reduction in this context is open for further investigations.

2.1 Geometry

We denote by $\Omega:=\{{\bm{x}}\in\mathbb{R}^{2}\ :\ R_{i}<\lvert{\bm{x}}\rvert<R_{o}\}$, $0<R_{i}<R_{o}<\infty$, the annular region. The liquid crystal is confined in the domain $\Omega_{s}:=\Omega\times(-s,s)$, where $s>0$ stands the film thickness. We also denote by $K_{i}$ and $K_{o}$ the regions in $\Omega$ occupied by the inner and outer electrodes respectively: $K_{i}:=\{{\bm{x}}\in\mathbb{R}^{2}\ :\ |{\bm{x}}|<R_{i}\}$ and $K_{o}:=\{{\bm{x}}\in\mathbb{R}^{2}\ :\ |{\bm{x}}|>R_{o}\}$.

As already mentioned, we consider two experimental settings: $K_{i}\times\{0\}$ and $K_{o}\times\{0\}$ (slim electrodes) or $K_{i}^{\infty}:=K_{i}\times\mathbb{R}$ and $K_{o}^{\infty}:=K_{o}\times\mathbb{R}$ (infinite electrodes). Generically, we use the notation $\mathcal{K}_{o}$ and $\mathcal{K}_{i}$ to denote either the slim or infinite electrodes. Furthermore, we set $D_{s}:=\mathbb{R}^{3}\setminus(\overline{\Omega}_{s}\cup\overline{\mathcal{K}}_{i}\cup\overline{\mathcal{K}}_{o})$, with $D:=D_{0}$, to denote the free space.

2.2 Electro-magnetism

The liquid crystals considered exhibit poor conductivity in the xy-plane, thereby resulting in low currents inside the liquid film and negligible magnetic field. Without magnetic effects, the electric field satisifes

[TABLE]

which guarantees the existence of a potential function $\Psi:\mathbb{R}^{3}\rightarrow\mathbb{R}$ satisfying

[TABLE]

The relation between the potential $\Psi$ and the surface charge density $q:\partial\Omega_{s}\rightarrow\mathbb{R}$ (units: $\mathrm{C}\text{\,}{\mathrm{m}}^{-2}$) is derived from the Gauss law. Because there is no charge in the free space $D_{s}$, it directly implies that

[TABLE]

In addition, the electric potential $\Psi$ is set to the appropriate voltages on the electrodes

[TABLE]

Behavior of the electric potential at $z\to\infty$ is required to close the system. In the case of slim electrodes, one assumes that

[TABLE]

while

[TABLE]

is imposed in the infinite electrodes case.

For the liquid region $\Omega_{s}$, we recall that even with low conductivity, the charges are located at the boundary of $\Omega_{s}$. To derive a relation between the electric field and the surface charge density, we consider a “Gaussian pillboxes” (see Figure 2 and also Griffiths [21])

[TABLE]

centered at $(\bar{x},\bar{y})\in\Omega$ with $r=s$. Applying the Gauss law on $G(\bar{x},\bar{y},r)$ yields

[TABLE]

Here ${\bm{e}}_{3}=(0,0,1)$ and $\varepsilon_{0}$ is the permittivity of the free space. Notice that the constant 2 is due to contribution from the top and bottom sides of the film. Also, the permittivity of the liquid crystal is assumed to be isotropic in the $xy$ directions and therefore the contributions from the other sides are of $O(s)$ in (2.5).

2.3 Fluid Dynamics

In the Smectic-A phase, the liquid crystal behaves like a two-dimensional incompressible fluid. The Navier-Stokes system is used to model the relation between the pressure $p(t):\Omega_{s}\rightarrow\mathbb{R}$ and the fluid velocity ${\bm{u}}(t):\Omega_{s}\rightarrow\mathbb{R}^{3}$

[TABLE]

where $\rho$ is the fluid density, $\eta$ is the shear viscosity and ${\bm{L}}:=q\nabla\Psi$ is the Lorentz force induced by the motion of the electric current. The velocity is two dimensional ${\bm{u}}\cdot{\bm{e}}_{3}=0$ in $\Omega_{s}$ (and in particular on the top and bottom sides). At the electrodes $\partial\Omega_{s}\cap\partial\mathcal{K}_{i/o}$, we restrict further the motion and consider a no-slip conditions ${\bm{u}}=0$.

2.4 Small Thickness Limiting Model

We now consider the limiting model when the thickness $s$ tends to [math]. Relation (2.5) between the electric field ${\bm{E}}=-\nabla\Psi$ and the surface charge density $q$ becomes

[TABLE]

Notice that by symmetry, relations (2.1), (2.7) and (2.2) reduces to a system of partial differential equations on the half plane $D^{+}:=D\cap\{z\geq 0\}$. It reads

[TABLE]

together with the boundary conditions

[TABLE]

where $\mathcal{K}_{i/o}^{+}:=\mathcal{K}_{i/o}\cap\{z\geq 0\}$.

Thanks to the two dimensional assumption on the fluid motion (Smectic-A phase), the Navier-Stokes system (2.6) reduces to

[TABLE]

where $p_{\Omega}(t):=p(t)|_{\Omega}:\Omega\rightarrow\mathbb{R}$, ${\bm{u}}_{\Omega}(t):\Omega\rightarrow\mathbb{R}^{2}$ denotes the first two components of ${\bm{u}}$ restricted to $\Omega$ and $\rho_{\Omega}$, $\eta_{\Omega}$ are the two-dimensional fluid mass density (units: $\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$) and two-dimensional fluid shear viscosity (units: $\mathrm{Pa}\text{\,}\mathrm{s}\text{\,}\mathrm{m}$) respectively. Also ${\bm{L}}_{\Omega}$ is the projection of the Lorenz force ${\bm{L}}$ on the plane supporting $\Omega$, i.e. ${\bm{L}}_{\Omega}=q\nabla_{\Omega}\Psi_{\Omega}$. In addition, the boundary conditions on the velocity become

[TABLE]

In $\Omega$, the charge density flux and electric field are assumed to satisfy the Ohm’s law, i.e., across any closed curves with outside pointing normal ${\bm{\nu}}$ we have

[TABLE]

where $\sigma_{\Omega}$ stands for the two-dimensional fluid electrical conductivity (units: $\mathrm{S}$) and ${\bm{E}}_{\Omega}$ denotes the first two component of ${\bm{E}}$ (projection onto the plan supporting $\Omega$). With this assumption, the conservation of surface charge density reads

[TABLE]

or

[TABLE]

using the electric potential $\Psi_{\Omega}:=\Psi|_{\Omega}$. The surface charge density $q$ in (2.12) is defined up to a constant fixed by the total charge conservation relation

[TABLE]

In summary, the electric potential $\Psi$, the surface charge density $q$, the velocity and pressure $({\bm{u}},p)$ are related by the system of differential equations (2.8), (2.10), (2.12), and (2.13), which is supplemented with the boundary conditions (2.9) and (2.11). From now on, we drop the subscript $\Omega$ on ${\bm{u}}_{\Omega}$, $p_{\Omega}$.

2.5 Nondimensional Model

To sort out the effects of the different parameters present in the model, we follow Tsai [32] and rewrite the governing equations (2.8), (2.7), (2.12), and (2.10) using the rescaled variables

[TABLE]

where $d:=R_{o}-R_{i}$ is the distance between the two electrodes. Similarly for the electrodes, we set $\widehat{\mathcal{K}}_{i/o}^{+}:=\{\hat{\mathbf{x}}\ :\mathbf{x}\in\mathcal{K}_{i/o}^{+}\}$. In addition, we define the rescaled functions

[TABLE]

In order to simplify the notations, we omit the notation $\hat{.}$ and for now on only consider the rescaled variables. In particular, we write

[TABLE]

with $\int_{\Omega}q(t)=\int_{\Omega}q(0)$ for the conservation of surface charge density relation,

[TABLE]

with boundary conditions

[TABLE]

(together with (2.3) in the case of slim electrodes and (2.4) in the case of infinite electrodes) for the relations between the electric potential and the surface charge density, and

[TABLE]

in $\Omega\times\mathbb{R}_{+}$ with boundary conditions ${\bm{u}}=0$ on $\partial\Omega$ for the Navier-Stokes system. The two dimensionless parameters ${\mathcal{P}}$ and $\mathcal{R}$ appearing in (2.17) are the Prandtl and Rayleigh numbers. They are given by

[TABLE]

The Prandtl number indicates the fluid viscous relaxation ability relative to its charge relaxation ability while the Rayleigh relates the electric with the dissipation forces. In Section 5, we propose a numerical study determining what range of parameters allows for electroconvection. We also introduce a geometric characteristic parameter $\alpha:=R_{i}/R_{o}\in(0,1)$ so that

[TABLE]

As we shall see in Section 5, this parameter affects the number of vortices during the electroconvection process.

3 Model Reduction

The electroconvection system (2.14)-(2.17) is mainly two dimensional (i.e. defined on $\Omega$) except for (2.15), where the electric potential must be computed in the entire free space $D^{+}$. However, notice that only its trace in $\Omega$ is required in (2.14) and (2.17). This, together with the two different electrodes configurations (infinite and slim), is exploited in Sections 3.1 and 3.2 to replace (2.15) with nonlocal problems in $\Omega$.

3.1 Infinite Electrodes

We assume that the two electrodes $K_{i}$ and $K_{o}$ extend to infinity along the $z$ directions, refer to Figure 3 for an illustration. This is the setting considered for instance in Constantin et al. [11]. However, because we do not impose vanishing charge densities on $\partial\Omega$, our model do not reduce to the one analyzed in Constantin et al. [11]. Worse, it appears that less regular charge densities are to be expected (see Remark 3.1).

We decompose $\Psi$ in two parts $\Psi=\Psi_{0}+\Psi_{K}$. The component $\Psi_{K}$ accounts for the voltage set by the electrodes and is defined as the solution to

[TABLE]

together with $\Psi_{K}=1$ on $\mathcal{K}_{i}^{+}$, $\Psi_{K}=0$ on $\mathcal{K}_{o}^{+}$ and $\lim_{z\downarrow 0}\partial_{z}\Psi_{K}=0$ on $\Omega$. Notice that $\Psi_{K}$ is independent of the $z$ variable. In fact, its exact expression is given by

[TABLE]

The second part $\Psi_{0}=\Psi-\Psi_{K}$ depends on the charge density $q$ and solves

[TABLE]

together the boundary conditions

[TABLE]

Following Stinga and Torrea [28], we realize that the trace $\Psi_{0,\Omega}:=\Psi_{0}|_{\Omega}$ satisfies the following non-local partial differential equation on $\Omega$

[TABLE]

where $(-\Delta_{\Omega})^{\frac{1}{2}}$ is the spectral Laplacian defined using a spectral expansion as described now. We denote by $\{\phi_{j}\}_{j=1}^{\infty}\subset H^{1}_{0}(\Omega)$ an $L_{2}(\Omega)$-orthonormal basis of eigenfunctions of $-\Delta$ and by $\{\lambda_{j}\}_{j=1}^{\infty}$ the associated eigenvalues. The spectral Laplacian is then defined for $v\in H^{1}_{0}(\Omega)$ via the relation

[TABLE]

Returning to the expression for $\Psi$, we find that its restriction to $\Omega$, $\Psi_{\Omega}:=\Psi|_{\Omega}$ satisfies

[TABLE]

where the expression of $\eta$ is given in  (3.2). The reduced model system, defined on $\Omega$, consists of (2.14), (2.17), (3.3) and (3.5).

Remark 3.1

We have already mentioned that compared with Constantin et al. [11], we do not impose $q|_{\partial\Omega}=0$ but rather $\int_{\Omega}q=0$. In particular, the charge density does not belong to $H^{1}_{0}(\Omega)$, the domain of the operator $(-\Delta_{\Omega})^{\frac{1}{2}}$, and we cannot use the relation

[TABLE]

to simplify the charge relation (2.14) as

[TABLE]

The former is the starting point of the analysis proposed in Constantin et al. [11]. We did not make such assumption because it implies

[TABLE]

which is incompatible with the surface charge density conservation required in our context.

Remark 3.2

The decomposition (3.5) of $\Psi_{\Omega}$ happens to correspond to the definition of the fractional laplacian with non vanishing boundary condition proposed in Antil et al. [2].

3.2 Slim electrodes

Instead of assuming infinite electrodes in the $z$ direction, we consider electrodes with negligible height as illustrated in Figure 4.

In this setting as well, we separate the voltage and charge contributions in the electric potential

[TABLE]

In view of (2.15), (2.16) and the decay relation $\lim_{z\to+\infty}\Psi=0$, the first part $\Psi_{K}$ is defined as the solution of

[TABLE]

together with the boundary condition

[TABLE]

and decay condition

[TABLE]

Here $\eta:\mathbb{R}^{2}\rightarrow\mathbb{R}$ matches the voltage imposed by the electrodes and is extended harmonically in $\Omega$:

[TABLE]

Its exact expression on $\Omega$ matches (3.2).

The second component $\Psi_{0}=\Psi-\Psi_{K}$ satisfies

[TABLE]

together with the decay condition $\lim_{z\to+\infty}\Psi_{0}=0$. It depends on the value of $\partial_{z}\Psi_{K}|_{\Omega}$ we determine now. Standard techniques based on a separation of variables $\Psi_{K}(x,y,z)=X(x,y)Z(z)$ and the representation of $X(x,y)$ in terms of Bessel’s functions of the first kind (cf. Evans [18, Chapter 4.1] and Tsai [32]) reveal that together with the axy-symmetry property of $\Psi_{K}|_{\Omega}=\eta$, the solution to (3.7), (3.8), (3.9) expressed in cylindrical coordinates $(r,\theta,z)$ takes the form

[TABLE]

The function $A_{0}$ is determined using the boundary condition $\Psi_{K}(r,\theta,0)=\eta(r)$ and the orthogonality properties of $J_{0}$:

[TABLE]

where we used the expression (3.2) of $\eta$ to derive the second equality. We now take advantage of Formula 6.574.1 in Gradshteyn and Ryzhik [20] to compute

[TABLE]

where the function ${}_{2}F_{1}$ denotes the hyper-geometric function, see Abramowitz and Stegun [1, Chapter 15]. In short, we write

[TABLE]

where

[TABLE]

Notice that $g\in L^{2}(\Omega)$ thanks to the property

[TABLE]

As a consequence, the system (3.10) for $\Psi_{0}$ is uniquely determined by $q$ and $g$.

We denote by $\Psi_{0,\Omega}$, the trace of $\Psi_{0}$ on $\Omega$, and define $\widetilde{\Psi}_{0,\Omega}$ to be its zero extension to $\mathbb{R}^{2}$, i.e.

[TABLE]

It turns out that $\widetilde{\Psi}_{0,\Omega}$ satisfies the following nonlocal problem

[TABLE]

where the fractional operator $(-\Delta_{\mathscr{F}})^{\frac{1}{2}}$ is the integral fractional laplacian defined via the Fourier transform $\mathscr{F}$ in $\mathbb{R}^{2}$ as detailed now. For $v\in H^{1}(\mathbb{R}^{2})$, the integral fractional Laplacian is a pseudo-differential operator with symbol $\lvert\zeta\rvert$, i.e.

[TABLE]

The claim that the restriction to $\Omega$ of the solution to (3.10) is the unique solution to (3.13) directly follows from an argument provided in Caffarelli and Silvestre [10], where $\partial_{z}\Psi_{0}$ is imposed on $\mathbb{R}^{2}$ rather than only in $\Omega$ as in our context. Following the argument provided in Section 3.2 in Caffarelli and Silvestre [10], it is easy to see that the unique solution $\Psi$ to (3.10) minimize the energy

[TABLE]

over all functions $\Phi:\mathbb{R}^{3}\rightarrow\mathbb{R}$ vanishing in $\mathbb{R}^{3}\setminus\Omega$ and when $z\to\infty$. Furthermore, its restriction $\Phi(z,y,0)$ also minimizes the energy

[TABLE]

over all functions $w:\mathbb{R}^{2}\rightarrow\mathbb{R}$ vanishing on $\mathbb{R}^{2}\setminus\Omega$. Since the solution to (3.13) is the unique minimizer of $J(.)$, the claim follows. Details are available in the Appendix A of Wei [33].

Combining the relations obtained for $\Psi_{K}$ and $\Psi_{0}$, we find that the restriction of $\Psi_{\Omega}=\Psi|_{\Omega}$ is given by

[TABLE]

where $\Psi_{0,\Omega}$ satisfies (3.13) with $g$ as in (3.12) and $\eta$ is given by (3.2).

The reduced model for the slim electrodes case defined on $\Omega$ consists of (2.14), (2.17) and (3.15).

4 Numerical Algorithms

In this section we detail the numerical algorithms advocated to approximate the fluid dynamic (2.17), the surface charge density convection (2.14), and the two non-local problems for the electric field (3.5) or (3.15) depending on the assumption made on the electrodes. In fact, the time marching algorithm proposed consists of three sub-steps. First, the electric potential is approximated in $\Omega$ using the (previous) surface charge density. Second, the surface charge density approximation is updated using the electric potential and (previous) fluid velocity. Third, the fluid velocity (and pressure) is updated with a Lorentz force computed using the surface charge density and electric potential. We detail each step separately below. Notice that the initial surface charge density and velocity are given allowing the algorithm to start.

4.1 Approximation of the Electric Potential (3.3) or (3.13)

We discuss the two different configurations separately.

Infinite Electrodes

We start with the simpler case of infinite electrodes, where the electric potential satisfies (3.3) involving the spectral Laplacian. Recall that $\Psi_{\Omega}=\Psi_{0,\Omega}+\eta$ with $\eta$ given by (3.2).

We adopt the numerical procedure put forward in Bonito and Pasciak [4, 5], which relies on the integral representation (Balakrishnan formula) of the electric potential $\Psi_{0,\Omega}$

[TABLE]

Here $\Phi(s;q)$ solves

[TABLE]

An exponentially convergent sinc quadrature, see Bonito and Pasciak [4, 5], Bonito et al. [6], is used to approximate the undefined integral. It reads

[TABLE]

where $k>0$ is the quadrature step size, $s_{j}:=jk$, $N_{+}:=\left\lceil\frac{\pi^{2}}{k^{2}}\right\rceil,\hbox{ and }N_{-}:=\left\lceil\frac{\pi^{2}}{k^{2}}\right\rceil$. In the numerical simulations proposed in Section 5, we fix $k=0.04$ so that $N_{-}=N_{+}=62$.

The solutions $\Phi(s_{j},q)$, $j=-N_{-},...,N_{+}$ of the subproblems (4.2) are in turn approximated by continuous piecewise linear functions. We start with a coarse polygonal approximation of $\Omega$ subdivided into quadrilaterals as in Figure 5(a). This polygonal approximation is then uniformly refined using quad-refinement but placing the boundary vertices on the exact boundary of $\Omega$. This gives rise to a polygonal domain $\Omega_{h}$ and a partition $\mathcal{T}_{h}$ made of quadrilaterals without hanging nodes, see for instance Figure 5(b). Here $h$ denotes the maximum diameter of elements on $\mathcal{T}_{h}$.

The number of quadrilaterals in the coarse subdivision depends on the aspect ratio in order to maximize the quality of the subdivision, see GridGenerator::hyper_shell documentation in Arndt et al. [3]. Furthermore, the final resolution $h$ required for the simulations presented in Section 5 depends on the aspect ratio $\alpha$ as well. Indeed, we shall see in Section 5.3 that larger aspect ratios yield more pairs of vortices during electroconvection, thereby increasing the resolution required. The number of uniform refinements performed on the coarse subdivisions for each aspect ratio $\alpha$ are given in Table 1.

For each $T\in\mathcal{T}_{h}$, we denote by $F_{T}:[0,1]^{2}\rightarrow T$ the bilinear (invertible) map from the unit square to $T$. This map is instrumental to define the finite element space

[TABLE]

With these notations, the finite element approximation of $\Psi_{0,\Omega}^{k}$ in (4.3) is given by

[TABLE]

where $\Phi_{h}(s_{j};q)\in V_{h}$ solves

[TABLE]

Returning to the decomposition (3.5), we arrive at

[TABLE]

where $\pi_{h}$ stands for the $L^{2}$ projection onto $V_{h}$.

It is worth mentioning, that the $N_{-}+N_{+}+1$ finite element problems (4.5) are mutually independent making the parallel implementation straightforward. Also, the algorithm consists of an outer loop (4.4) gathering the contributions of the finite element solutions at each quadrature point $s_{j}$. In particular, it only requires the implementation of a classical finite element solver for diffusion-reaction problems.

Slim electrodes

We now consider the case of slim electrodes. In this case, the electric potential on $\Omega$ is given by $\Psi_{\Omega}=\Psi_{0,\Omega}+\eta$ where $\Psi_{0,\Omega}$ satisfies (3.13). Unfortunately, an integral representation of $\Psi_{0,\Omega}$ like (4.1) is not available for the integral Laplacian. Instead, we follow Bonito et al. [7] and take advantage of an integral representation of the action of $(-\Delta_{{\mathscr{F}}})^{\frac{1}{2}}$ by multiplying (3.13) by a smooth function compactly supported in $\Omega$ and integrating over $\Omega$. Theorem 4.1 in Bonito et al. [7] guarantees that $\Psi_{0,\Omega}$ satisfies

[TABLE]

for all smooth functions $\varphi$ compactly supported in $\Omega$. Here $\Phi(s;\Psi_{0,\Omega})\in H^{1}(\mathbb{R}^{2})$ solves for all $\varphi\in H^{1}(\mathbb{R}^{2})$

[TABLE]

Similarly to the infinite electrode case, we can now use a sinc quadrature and standard finite element methods to provide an approximation of (4.7). The additional caveat stem from the fact that the domain of integration in (4.8) is $\mathbb{R}^{2}$. To circumvent this issue, we recall that

[TABLE]

and truncate the domain of integration to the disks

[TABLE]

where $M>1$ is a given truncation parameter.

We postpone for the moment the discussion regarding the automatic subdivision of $\Omega^{M}(s)$ but denote by $\Omega^{M}_{h}(s)$ its polygonal approximation and by $V^{M}_{h}(s)$ the associated finite element space based on continuous piecewise bilinear finite elements vanishing on $\partial\Omega^{M}_{h}(s)$. The approximation $\Psi_{0,h}\in V_{h}$ of $\Psi_{0,\Omega}$ satisfying (4.7) is then given by the relations

[TABLE]

where $\Phi_{h}(s_{j};\Psi_{0,h}):=\Phi_{h}^{k,M}(s_{j};\Psi_{0,h})\in V_{h}^{M}(s_{j})$ solves

[TABLE]

for all $\varphi_{h}\in V_{h}^{M}(s_{j})$.

At this point, a few comments are in order. The numerical procedure described above yields an approximation of the bilinear form (4.7) associated with the relation (3.13) for the electric potential but not the approximation $\Psi_{0,h}$ of $\Psi_{0,\Omega}$ directly. However, this is sufficient in the context of residual based iterative methods (such as conjugate gradient implemented in the present context), where only the action of the operator is required to compute residuals. The accumulation point of the iterative algorithm is $\Psi_{0,h}\in V_{h}$.

For the numerical experiments provided in Section 5, we take $k=0.04$, $N_{-}=62$, $N_{+}=124$ and $M=3$. It is shown in Bonito et al. [7] the proposed numerical method to approximate $\Psi_{0,\Omega}$ is exponentially convergent in the sinc quadrature parameter $k$ and in truncation parameter $M$. It is of optimal order (depending on the regularity of the exact electric potential $\Psi_{0,\Omega}$) in the space discretization. We refer to Bonito et al. [7] for additional information.

Returning to the decomposition (3.15), we define

[TABLE]

where $\eta$ is given by (3.2).

We end the discussion by describing the automatic generation of polygonal approximation $\Omega_{h}^{M}(s)$ of $\Omega^{M}(s)$ and its associated subdivision constrained to be an extension of a subdivisons of $\Omega$. We denote by $R^{M}(s)$ the radius of the truncated domain $\Omega^{M}(s)$. We start with a coarse subdivision made of 25 quadrilaterals illustrated in Figure 6(a). To increase the accuracy while keeping the complexity under control, given $0<h\leq 1$, the refinement procedure consists of

refine the fluid domain $\Omega_{h}$ uniformly until all quadtilaterals have diameters smaller than $h$, with possible additional refinements on $\Omega_{h}^{M}(s)$ to keep the number of hanging nodes to a maximum of one by edge. During this refinement process, newly created vertices at the boundary of $\Omega_{h}$ are placed on the boundary of $\Omega$.

• 2.

An exponential grading is performed outside $\Omega_{h}$, i.e. vertices on the azimuthal directions are placed at radii $r_{j}$ with

[TABLE]

We refer to Figure 6(b) for an illustration.

As in the case of infinite elecrodes configuration, the final resolution depends on the aspect ratio $\alpha$ to accommodate for the number of vortices during electroconvection. In the simulations presented in Section 5, we performed 5 uniform refinements for $\alpha\in\{0.1,0.2,0.33,0.452\}$, 6 for $\alpha\in\{0.56,0.6446\}$ and $7$ for $\alpha=0.8$.

4.2 Approximation of the charge density (2.14)

We now discuss the approximation of the surface charge density satisfying (2.14). Because $\Psi_{K}$ is harmonic in $\Omega$ in both electrodes configurations,

[TABLE]

Hence, the surface charge density conservation relation reduces to

[TABLE]

Notice that for a given electric potential $\Psi_{0,\Omega}$, the above partial differential equation for $q$ is a standard transport equation. It is approximated with an explicit Runge-Kutta 2 method in time and standard finite elements in space.

We recall $\mathcal{T}_{h}$ denotes a partition of a polygonal approximation $\Omega_{h}$ of the liquid domain $\Omega$ (see Section 4.1). We introduce the finite element space for the surface charge density

[TABLE]

Hence, given an approximation $\Phi_{0,h}\in V_{h}$ of $\Phi_{0,\Omega}$ determined as discussed in Section 4.1, the approximation $q_{h}(t)\in Q_{h}$ of $q(t)$ is defined as the solution to

[TABLE]

where ${\bm{\nu}}_{h}$ is the outside pointing normal to $\Omega_{h}$ (defined almost everywhere). This relation is obtained upon multiplying the surface charge density equation and integrating by parts the electric potential term. We suplement (4.10) with an approximated initial condition $\pi_{h}q(x,y,0)$, where $\pi_{h}$ denotes the $L^{2}(\Omega)$ projection onto $Q_{h}$.

A Strongly Stability Preserving (SSP) two stage Runge-Kutta scheme proposed by Gottlieb and Gottlieb [19] is advocated for the time discretization. Set $q_{h}^{0}=\pi_{h}q(0)$ the $L^{2}$ projection of a given surface charge density onto $Q_{h}$ and $t_{n}:=n\tau$ for $n=0,1,2,\cdots$ for a time step parameter $\tau>0$. Given approximation $q_{h}^{n}\in V_{h}$, $\Psi_{0,h}^{n}\in V_{h}$ (given by (4.6) or (4.9) with $q$ replaced by $q_{h}^{n}$) and ${\bm{u}}_{h}^{n}$ (see below) of the charge density $q(t_{n})$, electric potential $\Psi(t_{n})$ and fluid velocity ${\bm{u}}(t_{n})$ respectively, we approximate $q(t_{n+1})$ by $q_{h}^{n+1}\in Q_{h}$ as follows. In the first stage we seek $\mu_{h}^{(1)}\in Q_{h}$ satisfying

[TABLE]

for all $\varphi_{h}\in Q_{h}$. In the second stage, we find $\mu_{h}^{(2)}\in Q_{h}$ solving (4.11) but with $q_{h}^{n}$ replaced by $\mu_{h}^{(1)}$. Then, we set

[TABLE]

It is well documented that the finite element approximations of transport equation might be polluted by spurious oscillations. To circumvent this issue, we include the smoothness-based second-order maximum principle preserving viscosity method proposed in Guermond and Popov [22]. The numerical parameters used in the artificial viscosity are those recommended in (5.4) in Guermond et al. [23].

4.3 Approximation of the Fluid Dynamic (2.17)

We discretize the fluid dynamic using backward differentiation scheme of order 2 (BDF-2) coupled with Taylor-Hood finite element approximations for the space discretization; see Brenner and Scott [9]. Given a subdivision $\mathcal{T}_{h}$ of $\Omega_{h}$ constructed as in Section 4.1, the finite element spaces for the velocity and pressure are defined by

[TABLE]

where $\mathcal{Q}^{i}$, $i=1,2$, stands for the space of polynomial of (total) degree $i$.

We start with ${\bm{u}}_{h}^{0}:=\pi_{h}{\bm{u}}(0)$, the $L^{2}$ projection of a given initial velocity ${\bm{u}}(0)$ onto ${\bm{W}}_{h}$. We assume that at time $t=t_{n}:=n\tau$ (for the same time stepping parameter $\tau$ used for the surface charge density approximation), we have obtained $\Psi_{h}^{n}\in V_{h}$ given by (4.6) or (4.9) with $q$ replaced by $q_{h}^{n}$ and $q_{h}^{n+1}\in Q_{h}$ given by (4.12). The approximation $({\bm{u}}_{h}^{n+1},p_{h}^{n+1})\in{\bm{W}}_{h}\times X_{h}$ of $({\bm{u}}(t_{n+1}),p(t_{n+1}))$, $t_{n+1}:=(n+1)\tau$, is then defined as satisfying

[TABLE]

and

[TABLE]

for all $({\bm{v}}_{h},\theta_{h})\in{\bm{W}}_{h}\times X_{h}$.

5 Numerical Simulations

In this section we discuss the numerical results for two different settings with a particular emphasis on the effects of the Prandtl number $\mathcal{P}$, the Rayleigh number $\mathcal{R}$, and the aspect ratio $\alpha$.

The initial setting is common to all experiments. The fluid is always starting at rest. For the surface charge density, we always start with an initial configuration in equilibrium with the electric potential associated with the slim electrode configuration. More precisely, the initial charge density $q_{0}$ is chosen so that $\Psi=\Psi_{K}$ (i.e. $\Psi_{0,\Omega}=0$) in the decomposition (3.6). This corresponds to setting $q_{0}=-g$, where $g$ is given by  (3.12). We provide in Figure 7 an illustration of the initial surface charge density for $\alpha=0.33$. It is worth mentioning that we numerically observed the need of such compatibility between $q_{0}$ and $\Psi_{0}$ in the slim electrode setting but have not investigated this analytical question further. Notice that this configuration is unstable as the electric charges are aggregated near the outer boundary where the voltage is minimal. In the simulation below, we break the symmetry by adding a Gaussian white noise of magnitude no larger than $10^{-4}$.

For comparison, we consider the same initial surface charge density in the infinite electrode case.

In order to detect when and whether electroconvection occurs, we monitor two quantities: the kinetic energy

[TABLE]

and the circulation energy

[TABLE]

Because of the added white noise to the initial surface charge density, the two quantities $E_{k}$ and $E_{curl}$ evolves initially. We declare the electroconvection phenomena to occur when both $E_{k}$ and $E_{curl}$ undergo a relative change greater than $0.1\%$ compared to their respective initial values before time $t=20$ (we ran our simulation further and did not observe any changes after that). We define the critical Rayleigh number $\mathcal{R}_{c}$ as the threshold value for which fluids with lower Rayleigh numbers electroconvection do not occur (the Lorentz force is not strong enough to overcome the electric and viscous dissipation). For Rayleigh numbers above ${\mathcal{R}}_{c}$, the axis-symmetric distribution of charge density is broken and vortices appear. The critical vortex pair number ${\mathcal{N}}_{c}$ is the number of pairs of vortices during electroconvection. As an illustration, we provide in Figures 8 and  9 the numerical approximation of sustained electroconvection with ${\mathcal{N}}_{c}=4$ and ${\mathcal{N}}_{c}=8$ respectively. Similar structures but with different ${\mathcal{N}}_{c}$ are observed for other different aspect ratios.

We summarize the numerical parameters using in all the simulation below in Table 2 for the infinite electrodes configuration and in Table 3 for the slim electrodes configuration.

5.1 Comparison between the infinite and slim electrodes models

For this comparison, we set the aspect ratio to $\alpha=0.33$ and the Prandtl number to ${\mathcal{P}}=10$. When the Rayleigh number is ${\mathcal{R}}=100$, electroconvection is observed in the slim electrodes configuration but not in the infinite electrodes configuration, which seems to requires more energy to trigger electroconvection. In fact, electroconvection is observed in the infinite electrodes configuration for ${\mathcal{R}}\geq 250$, see kinetic and circulation energies in Figure 10.

Moreover, we find that even when electroconvection occurs in the infinite electrodes configuration, it cannot be sustained as in the slim electrode configuration. To substantiate this fact, we set ${\mathcal{P}}=10$ and ${\mathcal{R}}=800$ and compare the energies in Figure 11 for the two configurations. We observe that the energies in the slim case are not only significantly larger, they remain large as time evolves unlike in the infinite electrode configuration.

The difference between the two models is striking. In the presence of infinite electrodes model, electroconvection occurs as well but at much larger Rayleigh number ${\mathcal{R}}$. Still, even when convective flows appears, they do not persist and quickly disappear, see Figure 11. In contrast, the capability of maintaining a stable electroconvection phenomena indicates that the slim electrodes model is more adequate for electroconvection. Therefore, from now on we only consider the slim electrodes case. Worth mentioning, both models predict the same the number of pairs of vortices ${\mathcal{N}}_{c}$, which seems indicating that ${\mathcal{N}}_{c}$ depends mainly on the geometry as discussed later in Section 5.3.

5.2 Effect of the Prandtl number for the slim electrodes case

The Prandtl number ${\mathcal{P}}$ is the dimensionless ratio between the charge and viscous relaxation times. To understand its influence in the electroconvection phenomena, we fix $\alpha=0.33$, $\mathcal{R}=100$ and let $\mathcal{P}$ vary from $0.01$ to $1,000$. In Figure 12 we report the kinetic energies and the circulation energies of the fluid. We observed that the electroconvection occurs for value of $\mathcal{P}\in[0.1,1000]$ which influences the activation time but not the long term behavior. Worth noting, the energy plots in Figure 12 seem also to indicate that electroconvection cannot occur before $t=5$ irrespectively of the value of the Prandtl number. For smaller Prandtl number (${\mathcal{P}}=0.01$ as shown in Figure 12), electroconvection is not observed numerically. The existence of a critical Prandtl number is in accordance with the nonlinear analysis provided in Tsai et al. [30]. For the rest of this section we will be focusing on the value ${\mathcal{P}}=10$ for the purpose of demonstrating convective flow.

5.3 Effect of the Geometry for the slim electrodes case

The geometry of the annulus domain $\Omega$ is characterized by the aspect ratio $\alpha=R_{i}/R_{o}$; see (2.19). It turns out $\alpha$ affects the critical Rayleigh number $\mathcal{R}_{c}$ after which electroconvection occurs as well as the number of vortex pairs ${\mathcal{N}}_{c}$.

We fix ${\mathcal{P}}=10$ and determine for aspect ratio $\alpha\in[0.1,0.8]$ the corresponding critical Rayleigh number ${\mathcal{R}}_{c}$. The latter is determined by running independent simulations starting from a relatively small ${\mathcal{R}}$ at which electroconvection does not occur, and continuously increasing the value ${\mathcal{R}}$ by one unit each time until electroconvection occurs. This allows us to determine the critical value ${\mathcal{R}}_{c}$ up to $1$ unit. The system is considered steady when $E_{k}$ and $E_{curl}$ remains within $0.1\%$ relative difference throughout the simulation time. For example, Figure 13(a) depicts the evolution of $E_{k}$ for $\alpha=0.33$ with $\mathcal{R}$ varying from $80$ to $87$. The critical Rayleigh number satisfies $82<{\mathcal{R}}_{c}<83$.

In Figure 14(a) we report all the critical Rayleigh values for various $\alpha$ and compare them with Tsai et al. [30]. Although not strictly matching the simulation results in Fig. 5 of Tsai et al. [30], our results are in good agreement. In fact, in all cases, our numerical uncertainty intervals encompasses the theoretically predicted ranges reported in Fig. 5(a) of Tsai et al. [30].

The number of vortex pairs ${\mathcal{N}}_{c}$ are strongly influenced by the geometry. To document this we fix $\mathcal{P}=10$ and set, for each aspect ratio considered, the Rayleigh number at the critical value $\mathcal{R}_{c}$. In Figure 14(b) we compare the number of vortex pairs ${\mathcal{N}}_{c}$ obtained by our algorithm with results from Tsai et al. [30]. They match for all aspect ratio considered except for $\alpha=0.56$ where they differ by one. However, our predicted number of vortex pair in this case matches the experimental data Figure 5(b) of Tsai et al. [30].

5.4 Effect of the Rayleigh number for the slim electrodes case

From the definition (2.18), we realize that ${\mathcal{R}}\propto V^{2}$. Increasing the Rayleigh number corresponds to a stronger electric field and thus a stronger Lorentz force. In Section 5.3, we have already discussed the influence of the geometry on the critical Rayleigh value. We now set $\alpha=0.33$ and ${\mathcal{P}}=10$ and complete the investigation by increasing the value of $\mathcal{R}$ up to $10\mathcal{R}_{c}$. The values of $E_{k}$ and $E_{curl}$ are reported in Figure 15. We observe that larger Rayleigh numbers result in faster activation of the convection. This is similar as for the Prandtl number discussed in Section 5.2 but in this case there does not seem to be a limiting time for before which electroconvection cannot occur. Consequently fluids with larger Rayleigh numbers develop a stable electroconvection at earlier times. The strength of convection is also stronger for larger Rayleigh number.

6 Conclusions

In this paper, we have derived a mathematical model for electrically driven convection in an annular two dimensional fluid. The two different electrodes configurations are considered: infinite and slim. Depending on the electrodes configuration, nonlocal representations of the electric potential are derived on the liquid domain. This together with the surface density charge conservation relation and the Navier-Stokes system for the fluid dynamics yield a system of partial differential equations defined only in the two dimensional and bounded liquid region.

The proposed numerical methods take advantage of this dimension reduction and only require discretization tools readily available in most finite element codes. Our numerical simulations reveal that the slim electrodes configuration is more favorable for electroconvection: it requires less energy and is able to sustain the effect. It is therefore the configuration chosen to provide a numerical study of the three nondimensional parameters describing the electroconvection system.

We observe the Prandtl number must be above a critical value for electroconvection occurs. However, values above such critical number only affects the electroconvection activation time but not the long term behavior. For several geometries, we also exhibit critical Rayleigh numbers above which electroconvection occurs.