t1t1footnotetext: A.B. is partially supported by NSF grant DMS-1817691. P.W.
is partially supported by DMS-1254618, DMS-1817691 and the TAMU Hagler
Institute for Advanced Study HEEP Graduate Fellowship.

# Electroconvection of Thin Liquid Crystals: Model Reduction and Numerical
Simulations

Andrea Bonito<EMAIL_ADDRESS>Peng Wei<EMAIL_ADDRESS>Department of
Mathematics, Texas A&M University, College Station, TX 77843, USA;
<EMAIL_ADDRESS>
###### 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.

###### keywords:

Electroconvection, Fractional Operators, Finite Elements, Rayleigh number,
Prandlt number

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

|
---|---
(a) | (b)
Figure 1: Experimental settings. The liquid film is located in between two
concentric electrodes. (a) the two electrodes extend to infinity (infinite
case) and (b) the two electrodes have negligible thickness (slim case). Notice
that in both cases, the outer electrode is extended to infinity in the
$xy$-plane (not pictured).

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

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

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

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

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

(2.1)
$\nabla{\cdot}{\bm{E}}=0\qquad\textrm{or}\qquad\Delta\Psi=0\qquad\textrm{in}\quad
D_{s}.$

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

(2.2)
$\Psi=V\quad\textrm{on}\quad\mathcal{K}_{i}\qquad\textrm{and}\qquad\Psi=0\quad\textrm{on}\quad\mathcal{K}_{o}.$

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

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

while

(2.4) $\lim_{z\to\infty}|\Psi(x,y,z)|<\infty,\qquad(x,y)\in\Omega,$

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])

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

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

(2.5)
$\varepsilon_{0}({\bm{E}}(\bar{x},\bar{y},s)-{\bm{E}}(\bar{x},\bar{y},-s))\cdot{\bm{e}}_{3}=2q+O(s).$

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

Figure 2: Gaussian pillbox enclosing a small region of the fluid.

###  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}$

(2.6)
$\rho(\partial_{t}+{\bm{u}}\cdot\nabla){\bm{u}}+\eta\Delta{\bm{u}}+\nabla
p={\bm{L}},\quad\nabla\cdot{\bm{u}}=0,\quad\textrm{both in
}\Omega_{s}\times\mathbb{R}_{+},$

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 $0$.
Relation (2.5) between the electric field ${\bm{E}}=-\nabla\Psi$ and the
surface charge density $q$ becomes

(2.7) $\lim_{z\downarrow 0}\frac{\partial}{\partial z}\Psi-\lim_{z\uparrow
0}\frac{\partial}{\partial z}\Psi=-\frac{2q}{\varepsilon_{0}}\quad\textrm{in
}\Omega,$

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

(2.8) $\Delta\Psi=0,\quad\hbox{ in }D^{+},$

together with the boundary conditions

(2.9) $\begin{split}&\Psi=V\leavevmode\nobreak\
\textrm{on}\leavevmode\nobreak\ \mathcal{K}_{i}^{+},\leavevmode\nobreak\
\quad\Psi=0\leavevmode\nobreak\ \textrm{on}\leavevmode\nobreak\
\mathcal{K}_{o}^{+},\leavevmode\nobreak\ \textrm{and }\lim_{z\downarrow
0}\partial_{z}\Psi=-\frac{q}{\varepsilon_{0}}\leavevmode\nobreak\
\textrm{on}\leavevmode\nobreak\ \Omega,\end{split}$

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

(2.10)
$\begin{split}\rho_{\Omega}(\partial_{t}+{\bm{u}}_{\Omega}\cdot\nabla_{\Omega}){\bm{u}}_{\Omega}+\eta_{\Omega}\Delta_{\Omega}{\bm{u}}_{\Omega}+\nabla_{\Omega}p_{\Omega}&={\bm{L}}_{\Omega},\qquad\textrm{in
}\Omega\times\mathbb{R}_{+},\\\
\nabla_{\Omega}\cdot{\bm{u}}_{\Omega}&=0,\quad\qquad\textrm{in
}\Omega\times\mathbb{R}_{+},\end{split}$

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

(2.11)
${\bm{u}}_{\Omega}=0\qquad\textrm{on}\quad\partial\Omega\times\mathbb{R}_{+}.$

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

$\nabla_{\Omega}q\cdot{\bm{\nu}}=-\sigma_{\Omega}{\bm{E}}_{\Omega}\cdot{\bm{\nu}},$

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

$\partial_{t}q+\nabla_{\Omega}\cdot({\bm{u}}_{\Omega}q-\sigma_{\Omega}{\bm{E}}_{\Omega})=0,\qquad\textrm{in
}\Omega\times\mathbb{R}_{+}.$

or

(2.12)
$\partial_{t}q+\nabla_{\Omega}\cdot({\bm{u}}_{\Omega}q)-\sigma_{\Omega}\Delta_{\Omega}\Psi_{\Omega}=0,\qquad\textrm{in
}\Omega\times\mathbb{R}_{+},$

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

(2.13) $\int_{\Omega}q(t)=\int_{\Omega}q(0),\qquad t>0.$

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

$\hat{\mathbf{x}}:=\frac{\mathbf{x}}{d},\;\hat{t}:=\frac{\sigma_{\Omega}}{\varepsilon_{0}d}t,\quad\widehat{\Omega}:=\\{\hat{\mathbf{x}}\
:\mathbf{x}\in\Omega\\},\quad\textrm{and}\quad\widehat{D}^{+}:=\\{\hat{\mathbf{x}}\
:\mathbf{x}\in D^{+}\\},$

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

$\displaystyle\widehat{\Psi}(\hat{\mathbf{x}},\hat{t})$
$\displaystyle:=\frac{1}{V}\Psi(\mathbf{x},t),\quad\widehat{\Psi}_{\widehat{\Omega}}=\widehat{\Psi}|_{\widehat{\Omega}},\qquad\hat{q}(\hat{\mathbf{x}},\hat{t}):=\frac{d}{\epsilon_{0}V}q(\mathbf{x},t),$
$\displaystyle\widehat{{\bm{u}}}(\hat{\mathbf{x}},\hat{t})$
$\displaystyle:=\frac{\epsilon_{0}}{\sigma_{\Omega}}{\bm{u}}(\mathbf{x},t),\qquad\hat{p}(\hat{\mathbf{x}},\hat{t}):=\frac{\epsilon_{0}^{2}}{\sigma_{\Omega}^{2}\rho_{\Omega}}p(\mathbf{x},t).$

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

(2.14)
$\partial_{t}{q}+{\bm{u}}\cdot\nabla_{\Omega}q-\Delta_{\Omega}\Psi_{\Omega}=0,\qquad\hbox{
in }{\Omega}\times\mathbb{R}_{+},$

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

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

with boundary conditions

(2.16) $\begin{split}&\Psi=1\leavevmode\nobreak\
\textrm{on}\leavevmode\nobreak\ \mathcal{K}_{i}^{+},\leavevmode\nobreak\
\Psi=0\leavevmode\nobreak\ \textrm{on}\leavevmode\nobreak\
\mathcal{K}_{o}^{+},\leavevmode\nobreak\ \textrm{and}\leavevmode\nobreak\
\lim_{z\downarrow 0}\partial_{z}\Psi=-q\leavevmode\nobreak\
\textrm{on}\leavevmode\nobreak\ \Omega\end{split}$

(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

(2.17)
$\partial_{t}{\bm{u}}+({\bm{u}}\cdot\nabla_{\Omega}){\bm{u}}-{\mathcal{P}}\Delta_{\Omega}{\bm{u}}+\nabla_{\Omega}p=-{\mathcal{R}}{\mathcal{P}}q\nabla_{\Omega}\Psi_{\Omega},\qquad\nabla_{\Omega}\cdot{\bm{u}}=0,$

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

(2.18)
$\displaystyle{\mathcal{P}}:=\frac{\epsilon_{0}\eta_{\Omega}}{\rho_{\Omega}\sigma_{\Omega}d},\quad\hbox{
and
}\quad{\mathcal{R}}:=\frac{\epsilon_{0}^{2}V^{2}}{\eta_{\Omega}\sigma_{\Omega}}.$

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

(2.19) $R_{i}=\frac{\alpha}{1-\alpha},\qquad\textrm{and}\qquad
R_{o}=\frac{1}{1-\alpha}.$

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

Figure 3: Domain and boundaries, cross section view in the infinite electrodes
case; compare with Figure 4.

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

(3.1) $\Delta\Psi_{K}=0\hbox{ in }D^{+}\quad$

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

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

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

$-\Delta\Psi_{0}=0\hbox{ in }D^{+}$

together the boundary conditions

$\Psi_{0}=0\leavevmode\nobreak\ \textrm{on}\leavevmode\nobreak\
\mathcal{K}_{i}^{+}\cup\mathcal{K}_{o}^{+},\quad\textrm{and}\quad\lim_{z\downarrow
0}\partial_{z}\Psi_{0}=-q\quad\textrm{on }\Omega.$

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$

(3.3) $(-\Delta_{\Omega})^{\frac{1}{2}}\Psi_{0,\Omega}=q\quad\hbox{ in
}\quad\Omega,$

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

(3.4)
$(-\Delta_{\Omega})^{\frac{1}{2}}v:=\sum_{j=0}^{\infty}\lambda_{j}^{\frac{1}{2}}\left(\int_{\Omega}v\phi_{j}\right)\phi_{j}.$

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

(3.5) $\Psi_{\Omega}=\Psi_{0,\Omega}+\eta$

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

$-(\nabla_{\Omega}\cdot\nabla_{\Omega})(-\Delta_{\Omega})^{-\frac{1}{2}}=(-\Delta_{\Omega})(-\Delta_{\Omega})^{-\frac{1}{2}}=(-\Delta_{\Omega})^{\frac{1}{2}},$

to simplify the charge relation (2.14) as

$\partial_{t}{q}+{{\bm{u}}}\cdot\nabla_{\Omega}q+(-\Delta_{\Omega})^{\frac{1}{2}}q=0\quad\hbox{
in }\quad{\Omega}\times\mathbb{R}_{+}.$

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

$\lim_{t\to\infty}\int_{\Omega}q(t)=0,$

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.

Figure 4: Domain and boundaries, cross section view in the slim electrodes
case; compare with Figure 3.

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

(3.6) $\Psi=\Psi_{0}+\Psi_{K}.$

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

(3.7) $\Delta\Psi_{K}=0\quad\hbox{ in }\mathbb{R}^{2}\times\mathbb{R}_{+}$

together with the boundary condition

(3.8) $\Psi_{K}=\eta,\qquad\textrm{on }\mathbb{R}^{2},$

and decay condition

(3.9) $\lim_{z\to+\infty}\Psi_{K}=0.$

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

$-\Delta_{\Omega}\eta=0\quad\textrm{in }\Omega,\quad\eta=0\quad\textrm{on
}K_{o},\quad\textrm{and}\quad\eta=1\quad\textrm{on }K_{i}.$

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

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

(3.10) $\begin{split}\left\\{\begin{array}[]{rll}-\Delta\Psi_{0}&=0,&\hbox{ in
}D^{+},\\\ \partial_{z}\Psi_{0}&=-q-\partial_{z}\Psi_{K},&\hbox{ on
}\Omega,\\\ \Psi_{0}&=0,&\hbox{ on }K_{i}\cup
K_{o},\end{array}\right.\end{split}$

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

(3.11) $\Psi_{K}(r,\theta,z)=\int_{0}^{\infty}e^{-kz}J_{0}(kr)A_{0}(k)dk.$

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

$A_{0}(k)=k\int_{0}^{\infty}\eta(r)J_{0}(kr)rdr=\frac{1}{k^{2}\ln\alpha}\left(J_{0}(R_{0}k)-J_{0}(R_{i}k)\right),$

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

$\begin{split}\lim_{z\downarrow
0}\partial_{z}\Psi_{K}(r,\theta,z)&=-\int_{0}^{\infty}kJ_{0}(kr)A_{0}(k)dk\\\
&=-\frac{1}{\ln\alpha}\int_{0}^{\infty}(J_{0}(R_{o}k)-J_{0}(R_{i}k))J_{0}(rk)dk\\\
&=-\frac{1}{\ln\alpha}\left(\frac{1}{R_{o}}{}_{2}F_{1}(\frac{1}{2},\frac{1}{2};1;\frac{r^{2}}{R_{o}^{2}})-\frac{1}{r}{}_{2}F_{1}(\frac{1}{2},\frac{1}{2};1;\frac{R_{i}^{2}}{r^{2}})\right)\end{split}$

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

$\lim_{z\downarrow 0}\partial_{z}\Psi_{K}=g,$

where

(3.12)
$g(r):=-\frac{1}{\ln\alpha}\left(\frac{1}{R_{o}}{}_{2}F_{1}(\frac{1}{2},\frac{1}{2};1;\frac{r^{2}}{R_{o}^{2}})-\frac{1}{r}{}_{2}F_{1}(\frac{1}{2},\frac{1}{2};1;\frac{R_{i}^{2}}{r^{2}})\right).$

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

$\lim_{\beta\to
1^{-}}\frac{{}_{2}F_{1}(a,b;a+b;\beta)}{-\ln(1-\beta)}=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}.$

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.

$\widetilde{\Psi}_{0,\Omega}=\Psi_{0}\qquad\textrm{on}\quad\Omega,\qquad\textrm{and}\qquad\widetilde{\Psi}_{\Omega}=0\quad\textrm{on}\quad\mathbb{R}^{2}\setminus\Omega.$

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

(3.13)
$(-\Delta_{\mathscr{F}})^{\frac{1}{2}}\widetilde{\Psi}_{0,\Omega}=q+g\qquad\textrm{in}\quad\Omega,$

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.

(3.14)
$\mathscr{F}((-\Delta_{\mathscr{F}})^{\frac{1}{2}}v)(\zeta):=\lvert\zeta\rvert\mathscr{F}(v)(\zeta).$

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

$\frac{1}{2}\int_{\mathbb{R}^{2}}\int_{0}^{\infty}|\nabla\Phi(x,y,z)|^{2}\leavevmode\nobreak\
dz\leavevmode\nobreak\
dxdy-\int_{\Omega}(q(x,y)+g(x,y))\Psi(x,y,0)\leavevmode\nobreak\ dxdy$

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

$J(w):=\frac{1}{2}\int_{\mathbb{R}^{2}}|\xi||\mathscr{F}(w)|^{2}d\xi\
-\frac{1}{\pi}\int_{\Omega}(q(x,y)+g(x,y))w(x,y)\leavevmode\nobreak\ dxdy$

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

(3.15) $\Psi_{\Omega}=\Psi_{0,\Omega}+\eta$

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

(4.1)
$\Psi_{0,\Omega}=\frac{2}{\pi}\int_{-\infty}^{\infty}e^{s}\Phi(s;q)\,ds.$

Here $\Phi(s;q)$ solves

(4.2) $\Phi(s;q)-e^{2s}\Delta\Phi(s;q)=q\qquad\textrm{in
}\Omega,\quad\textrm{with}\quad\Phi(s;q)=0\qquad\textrm{on }\partial\Omega.$

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

(4.3)
$\Psi_{0,\Omega}\approx\Psi^{k}:=\frac{2}{\pi}k\sum_{j=-N_{-}}^{N_{+}}e^{s_{j}}\Phi(s_{j};q),$

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}$.

|
---|---
(a) |  (b)
Figure 5: Polygonal approximation $\Omega_{h}$ of $\Omega$ for an aspect ratio
$\alpha=0.33$. (a) The coarse initial subdivision; (b) Approximation resulting
in two successive uniform refinements and placing the boundary vertices on the
boundary of $\Omega$.

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.

$\alpha$ | 0.1 | 0.2 | 0.33 | 0.452 | 0.56 | 0.6446 | 0.9
---|---|---|---|---|---|---|---
$\\#\mathcal{T}$ | 4 | 5 | 6 | 9 | 12 | 15 | 29
$\\#$ uniform ref. | 5 | 5 | 5 | 5 | 5 | 5 | 6
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$.

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

$V_{h}:=\\{v\in C^{0}(\overline{\Omega}_{h})\ :\ v|_{T}\circ F_{T}\textrm{ is
bilinear for all }T\in\mathcal{T}_{h}\quad\textrm{and}\quad
v|_{\partial\Omega_{h}}=0\\}.$

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

(4.4)
$\displaystyle\Psi_{0,h}:=\Psi^{k}_{0,h}:=\frac{2k}{\pi}\sum_{j=-N_{-}}^{N_{+}}e^{s_{j}}\Phi_{h}(s_{j};q)\in
V_{h},$

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

(4.5)
$\int_{\Omega}\Phi_{h}(s_{j};q)\varphi_{h}+e^{2s_{j}}\int_{\Omega}\nabla\Phi_{h}(s_{j};q)\cdot\nabla\varphi_{h}=\int_{\Omega}q\varphi_{h},\hbox{\
\ for all }\varphi_{h}\in V_{h}.$

Returning to the decomposition (3.5), we arrive at

(4.6) $\Psi_{h}:=\Psi_{0,h}+\pi_{h}\eta\in V_{h},$

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

(4.7)
$\frac{1}{\pi}\int_{-\infty}^{\infty}e^{\frac{s}{2}}\left(\int_{\Omega}(\Phi(s;\Psi_{0,\Omega})+\Psi_{0,\Omega})\varphi\right)ds=\int_{\Omega}(q+g)\varphi$

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})$

(4.8)
$e^{s}\int_{\mathbb{R}^{2}}\Phi(s;\Psi_{0,\Omega})\varphi+\int_{\mathbb{R}^{2}}\nabla\Phi(s;\Psi_{0,\Omega})\cdot\nabla\varphi=-e^{s}\int_{\Omega}\Psi_{0,\Omega}\varphi.$

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

$\Omega=\\{(x,y)\in\mathbb{R}^{2}\ :\
\frac{\alpha}{1-\alpha}<\sqrt{x^{2}+y^{2}}<\frac{1}{1-\alpha}\\}$

and truncate the domain of integration to the disks

$\Omega^{M}(s):=\left\\{\begin{array}[]{ll}\\{(1+e^{-\frac{s}{2}}(M+1))(x,y)\;:\;\sqrt{x^{2}+y^{2}}<1/(1-\alpha)\\},&e^{-\frac{s}{2}}\geq
1,\\\
\\{(M+2)(x,y)\;:\;\sqrt{x^{2}+y^{2}}<1/(1-\alpha)\\},&e^{-\frac{s}{2}}<1,\end{array}\right.$

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

$\frac{k}{\pi}\sum_{j=-N_{-}}^{N_{+}}e^{\frac{s_{j}}{2}}\int_{\Omega_{h}}(\Phi_{h}(s_{j};\Psi_{0,h})+\Psi_{0,h})\varphi_{h}=\int_{\Omega_{h}}(q+g)\varphi_{h},\qquad\forall\varphi_{h}\in
V_{h},$

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

$\begin{split}e^{s_{j}}\int_{\Omega_{h}^{M}(s_{j})}\Phi_{h}(s_{j};\Psi_{0,h})\varphi_{h}+\int_{\Omega^{M}_{h}(s_{j})}\nabla\Phi_{h}(s_{j};\Psi_{0,h})&\cdot\nabla\varphi_{h}\\\
&=-e^{s_{j}}\int_{\Omega_{h}}\Psi_{0,h}\varphi_{h},\end{split}$

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

(4.9) $\Psi_{h}:=\Psi_{0,h}+\pi_{h}\eta\in V_{h},$

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

  * 1.

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

$r_{j}:=e^{jh_{0}}/(1-\alpha),j=1,2,\cdots,\lceil
M/h\rceil\quad\textrm{with}\quad h_{0}=\ln(R^{M}(s)(1-\alpha))h/M;$

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

|
---|---
|
(a) | (b)
Figure 6: Polygonal approximation $\Omega_{h}^{M}(s)$ (first row) and zoom of
the liquid region $\Omega_{h}$ (second row) with $M=3$ and $s=1$ associated to
$\Omega$ with aspect ratio $\alpha=0.33$. The approximation of the liquid
domain $\Omega_{h}$ is in gray. (a) Initial subdivision and (b) Three
successive iterations of the refinement procedure.

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,

$-\Delta_{\Omega}\Psi_{\Omega}=-\Delta_{\Omega}\Psi_{0,\Omega}.$

Hence, the surface charge density conservation relation reduces to

$\partial_{t}q+{\bm{u}}\cdot\nabla_{\Omega}q-\Delta_{\Omega}\Psi_{0,\Omega}=0\qquad\textrm{in
}\Omega\times\mathbb{R}_{+}.$

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

$Q_{h}:=\\{v\in C^{0}(\overline{\Omega}_{h})\ :\ v|_{T}\circ F_{T}\textrm{ is
bilinear for all }T\in\mathcal{T}_{h}\textrm{ and }\int_{\Omega_{h}}v=0\\}.$

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

(4.10)
$\begin{split}&\int_{\Omega_{h}}\partial_{t}q_{h}(t)\varphi_{h}+\int_{\Omega_{h}}{\bm{u}}_{h}(t)\cdot\nabla_{\Omega}q_{h}(t)\varphi_{h}\\\
&\qquad+\int_{\Omega_{h}}\nabla_{\Omega}\Psi_{0,h}\cdot\nabla_{\Omega}\varphi_{h}-\int_{\partial\Omega_{h}}\nabla_{\Omega}\Psi_{0,h}\cdot{\bm{\nu}}_{h}\varphi_{h}=0,\qquad\forall\varphi_{h}\in
Q_{h},\end{split}$

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

(4.11)
$\begin{split}\int_{\Omega_{h}}\mu_{h}^{(1)}\varphi_{h}=&\int_{\Omega_{h}}q_{h}^{n}\varphi_{h}-\tau\int_{\Omega_{h}}{\bm{u}}_{h}^{n}\cdot\nabla_{\Omega}q_{h}^{n}\varphi_{h}\\\
&-\tau\int_{\Omega_{h}}\nabla_{\Omega}\Psi_{0,h}\cdot\nabla_{\Omega}\varphi_{h}+\tau\int_{\partial\Omega_{h}}\nabla_{\Omega}\Psi_{0,h}\cdot{\bm{\nu}}_{h}\varphi_{h},\end{split}$

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

(4.12) $q_{h}^{n+1}:=\frac{1}{2}(q_{h}^{n}+\mu_{h}^{(2)}).$

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

$\displaystyle{\bm{W}}_{h}$ $\displaystyle:=\\{{\bm{v}}\in
C^{0}(\overline{\Omega}_{h})^{2}\ :\ {\bm{v}}|_{T}\circ
F_{T}\in{\mathcal{Q}}_{2}^{2},\leavevmode\nobreak\ \forall
T\in{\mathcal{T}}_{h},\leavevmode\nobreak\
{\bm{v}}|_{\partial\Omega_{h}}=0\\},$ $\displaystyle X_{h}$
$\displaystyle:=\\{\theta\in C(\overline{\Omega}_{h}):\
\int_{\Omega_{h}}\theta=0,\ \theta|_{T}\circ F_{T}\in{\mathcal{Q}}_{1},\forall
T\in{\mathcal{T}}_{h}\\},$

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

(4.13)
$\begin{split}\int_{\Omega_{h}}{\bm{u}}_{h}^{n+1}\cdot{\bm{v}}_{h}&+\frac{2}{3}\tau{\mathcal{P}}\int_{\Omega_{h}}\nabla{\bm{u}}_{h}^{n+1}\cdot\nabla{\bm{v}}_{h}-\frac{2}{3}\tau\int_{\Omega_{h}}\nabla{\cdot}{\bm{v}}_{h}p_{h}^{n+1}\\\
&=\int_{\Omega_{h}}\left(\frac{4}{3}{\bm{u}}_{h}^{n}-\frac{1}{3}{\bm{u}}_{h}^{n-1}\right)\cdot{\bm{v}}_{h}-\frac{2}{3}\tau\int_{\Omega_{h}}{\bm{v}}_{h}\cdot({\bm{u}}_{h}^{n}\cdot\nabla){\bm{u}}_{h}^{n}\\\
&\qquad\qquad\qquad-\frac{2}{3}\tau{\mathcal{R}}{\mathcal{P}}\int_{\Omega_{h}}q_{n}^{n+1}\nabla_{\Omega}\Psi_{h}^{n}\cdot{\bm{v}}_{h},\end{split}$

and

$\int_{\Omega_{h}}\nabla{\cdot}{\bm{u}}_{h}^{n+1}\theta_{h}=0$

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}$.

|
---|---
(a) | (b)
Figure 7: The initial condition of charge density $q_{0}=-g$ with $g$ defined
by (3.12) for $\alpha=0.33$. (a) Initial surface charge density distribution;
(b) Cross-section view of the initial surface charge density with
$r=\sqrt{x^{2}+y^{2}}$.

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

$E_{k}:=\frac{1}{2}\int_{\Omega}\lvert{\bm{u}}\rvert^{2}\ d{\bm{x}},$

and the circulation energy

$E_{curl}:=\int_{\Omega}|\nabla\times{\bm{u}}|^{2}\ d{\bm{x}}.$

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.

|
---|---
(a) | (b)
Figure 8: Electroconvection for $\mathcal{P}=10$, $\mathcal{R}=100$, and $\alpha=0.33$ at time $40$. (a) Numerical approximation of the velocity field ${\bm{u}}$; (b) Numerical approximation of the electric surface charge density distribution $q$. |
---|---
(a) | (b)
Figure 9: Electroconvection for $\mathcal{P}=10$, $\mathcal{R}=100$, and
$\alpha=0.56$ at time $40$. (a) Numerical approximation of the velocity field
${\bm{u}}$; (b) Numerical approximation of the electric surface charge density
distribution $q$.

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.

$\alpha$ | $\tau$ | DoFs for $\Psi_{0,\Omega}$ | DoFs for $q$ | DoFs for ${\bm{u}}$ 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 2: Parameter settings for infinite electrode simulations. Here “DoFs” stands for degrees of freedom. $\alpha$ | $\tau$ | $M$ | DoFs for $\Psi_{0,\Omega}$ | DoFs for $q$ | DoFs for ${\bm{u}}$ and $p$
---|---|---|---|---|---
$\leq 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
Table 3: Parameter settings for slim electrode simulations. Here “DoFs” stands
for degrees of freedom.

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

|
---|---
(a) | (b)
Figure 10: Comparison of the Kinetic energy (a) and circulation enegry (b) for
the infinite electrode model with $\mathcal{P}=10$, $\alpha=0.33$ and for
several values of $\mathcal{R}$.

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.

|
---|---
(a) | (b)
Figure 11: Comparison of the evolution the energies $E_{k}$ and $E_{curl}$
versus time for the slim and infinite electrodes configurations with
$\alpha=0.33$, ${\mathcal{R}}=800$ and ${\mathcal{P}}=10$. The slim
configuration energies are significantly larger and remain large during the
entire evolution indicating a sustained electroconvection phenomena.

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

|
---|---
(a) | (b)
Figure 12: Kinetic energy (a) and $E_{curl}$ (b) evolutions over time for
different Prandtl number for $\alpha=0.33$ and $\mathcal{R}=100$. Increasing
the Prandtl number decreases the time for the electroconvection to develop. We
observe that electroconvection cannot occur in this setting before $t=5$
irrespectively of the Prandtl number. Compare with Figure 15.

###  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$.

|
---|---
(a) | (b)
Figure 13: Kinetic energy (a) and circulation energy (b) for slim electrodes
with ${\mathcal{P}}=10$, $\alpha=0.33$ and different values of
${\mathcal{R}}$.

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].

|
---|---
(a) | (b)
Figure 14: Critical Rayleigh number (a) and number of vortex pairs (b) in the
slim electrode case for $\mathcal{P}=10$ and several aspect ratio $\alpha$.
Compare with the results provided in Tsai et al. [30]. The uncertainty
intervals in part (a) are due to the increment used in the critical Rayleigh
number exploration.

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.

|
---|---
(a) | (b)
Figure 15: Comparison of the Kinetic energy (a) and circulation energy (b) for
the slim electrode model with $\mathcal{P}=10$, $\alpha=0.33$ and for several
values of $\mathcal{R}$.

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

## References

  * Abramowitz and Stegun [1964] M. Abramowitz and I. A. Stegun.  _Handbook of mathematical functions: with formulas, graphs, and mathematical tables_ , volume 55.  Courier Corporation, 1964.
  * Antil et al. [2017] H. Antil, J. Pfefferer, and S. Rogovs.  Fractional operators with inhomogeneous boundary conditions: Analysis, control, and discretization.  _arXiv preprint arXiv:1703.05256_ , 2017.
  * Arndt et al. [2017] D. Arndt, W. Bangerth, D. Davydov, T. Heister, L. Heltai, M. Kronbichler, M. Maier, J.-P. Pelteret, B. Turcksin, and D. Wells.  The deal.II library, version 8.5.  _Journal of Numerical Mathematics_ , 2017.  doi: 10.1515/jnma-2016-1045.
  * Bonito and Pasciak [2015] A. Bonito and J. Pasciak.  Numerical approximation of fractional powers of elliptic operators.  _Mathematics of Computation_ , 84(295):2083–2110, 2015.
  * Bonito and Pasciak [2016] A. Bonito and J. E. Pasciak.  Numerical approximation of fractional powers of regularly accretive operators.  _IMA Journal of Numerical Analysis_ , page drw042, 2016.
  * Bonito et al. [2017a] A. Bonito, W. Lei, and J. E. Pasciak.  On sinc quadrature approximations of fractional powers of regularly accretive operators.  _arXiv preprint arXiv:1709.06619_ , 2017a.
  * Bonito et al. [2017b] A. Bonito, W. Lei, and J. E. Pasciak.  Numerical approximation of the integral fractional laplacian.  _arXiv preprint arXiv:1707.04290_ , 2017b.
  * Bonito et al. [2018] A. Bonito, J. P. Borthagaray, R. H. Nochetto, E. Otarola, and A. J. Salgado.  Numerical methods for fractional diffusion.  _Computing and Visualization in Science_ , 19(5-6):19–46, 2018.
  * Brenner and Scott [2007] S. Brenner and R. Scott.  _The mathematical theory of finite element methods_ , volume 15.  Springer Science & Business Media, 2007.
  * Caffarelli and Silvestre [2007] L. Caffarelli and L. Silvestre.  An extension problem related to the fractional laplacian.  _Communications in partial differential equations_ , 32(8):1245–1260, 2007.
  * Constantin et al. [2016] P. Constantin, T. Elgindi, M. Ignatova, and V. Vicol.  On some electroconvection models.  _Journal of Nonlinear Science_ , pages 1–15, 2016.
  * Daya et al. [1997] Z. A. Daya, S. W. Morris, and J. R. De Bruyn.  Electroconvection in a suspended fluid film: a linear stability analysis.  _Physical Review E_ , 55(3):2682, 1997.
  * Daya et al. [1998] Z. A. Daya, V. Deyirmenjian, S. W. Morris, and J. R. De Bruyn.  Annular electroconvection with shear.  _Physical review letters_ , 80(5):964, 1998.
  * Daya et al. [1999] Z. A. Daya, V. Deyirmenjian, and S. W. Morris.  Electrically driven convection in a thin annular film undergoing circular couette flow.  _Physics of Fluids_ , 11(12):3613–3628, 1999\.
  * Daya et al. [2002] Z. A. Daya, V. Deyirmenjian, and S. W. Morris.  Sequential bifurcations in sheared annular electroconvection.  _Physical Review E_ , 66(1):015201, 2002.
  * Deyirmenjian et al. [1997] V. Deyirmenjian, Z. A. Daya, and S. W. Morris.  Weakly nonlinear analysis of electroconvection in a suspended fluid film.  _Physical Review E_ , 56(2):1706, 1997.
  * Deyirmenjian et al. [2005] V. Deyirmenjian, Z. A. Daya, and S. W. Morris.  Codimension-two points in annular electroconvection as a function of aspect ratio.  _Physical Review E_ , 72(3):036211, 2005.
  * Evans [2010] L. C. Evans.  _Partial differential equations_.  American Mathematical Society, 2010.
  * Gottlieb and Gottlieb [2003] S. Gottlieb and L.-A. J. Gottlieb.  Strong stability preserving properties of runge–kutta time discretization methods for linear constant coefficient operators.  _Journal of Scientific Computing_ , 18(1):83–109, 2003.
  * Gradshteyn and Ryzhik [2014] I. S. Gradshteyn and I. M. Ryzhik.  _Table of integrals, series, and products_.  Academic press, 2014.
  * Griffiths [2005] D. J. Griffiths.  Introduction to electrodynamics, 2005.
  * Guermond and Popov [2017] J.-L. Guermond and B. Popov.  Invariant domains and second-order continuous finite element approximation for scalar conservation equations.  _SIAM Journal on Numerical Analysis_ , 55(6):3120–3146, 2017.
  * Guermond et al. [2018] J.-L. Guermond, M. Q. de Luna, B. Popov, C. E. Kees, and M. W. Farthing.  Well-balanced second-order finite element approximation of the shallow water equations with friction.  _SIAM Journal on Scientific Computing_ , 40(6):A3873–A3901, 2018.
  * Jabbour [2015] J. A. Jabbour.  _A matrix-free numerical bifurcation method applied to sheared annular electroconvection_.  Master thesis, Modelling and Computational Science, University of Ontario, 2015.
  * Kim et al. [2010] S. J. Kim, Y.-A. Song, and J. Han.  Nanofluidic concentration devices for biomolecules utilizing ion concentration polarization: theory, fabrication, and applications.  _Chemical Society Reviews_ , 39(3):912–922, 2010\.  doi: 10.1039/b822556g.
  * Kwak et al. [2013] R. Kwak, V. S. Pham, K. M. Lim, and J. Han.  Shear flow of an electrically charged fluid by ion concentration polarization: Scaling laws for electroconvective vortices.  _Physical Review Letters_ , 110(11), 2013.  doi: 10.1103/physrevlett.110.114501.
  * Morris et al. [1991] S. W. Morris, J. R. De Bruyn, and A. May.  Patterns at the onset of electroconvection in freely suspended smectic films.  _Journal of statistical physics_ , 64(5-6):1025–1043, 1991.
  * Stinga and Torrea [2010] P. R. Stinga and J. L. Torrea.  Extension problem and harnack’s inequality for some fractional operators.  _Communications in Partial Differential Equations_ , 35(11):2092–2122, 2010.
  * Tsai et al. [2004] P. Tsai, Z. A. Daya, and S. W. Morris.  Aspect-ratio dependence of charge transport in turbulent electroconvection.  _Physical review letters_ , 92(8):084503, 2004\.
  * Tsai et al. [2007] P. Tsai, Z. A. Daya, V. B. Deyirmenjian, and S. W. Morris.  Direct numerical simulation of supercritical annular electroconvection.  _Physical Review E_ , 76(2):026305, 2007.
  * Tsai et al. [2008] P. Tsai, S. W. Morris, and Z. A. Daya.  Localized states in sheared electroconvection.  _EPL (Europhysics Letters)_ , 84(1):14003, 2008\.
  * Tsai [2007] P.-C. Tsai.  _The Route to Chaos and Turbulence in Annular Electroconvection_.  University of Toronto, 2007.
  * Wei [2019] P. Wei.  _Numerical approximation of time dependent fractional diffusion with drift: Numerical analysis and applications to surface quasi-geostrophic dynamics and electroconvection_.  PhD thesis, Texas A&M University, 2019.