1. Introduction
1.1. General description of electrohydrodynamic flow
Electrohydrodynamics (EHD) is concerned with the interaction between an electric field and a flow field. Such configurations have broad applications in a range of industrial and biological devices. EHD effects can be used to enhance the heat transfer efficiency (Jones Reference Jones1978; Allen & Karayiannis Reference Allen and Karayiannis1995), to design microscale EHD pumps (Bart et al. Reference Bart, Tavrow, Mehregany and Lang1990; Darabi et al. Reference Darabi, Rada, Ohadi and Lawler2002), to fabricate diagnostic devices and drug delivery systems (Chakraborty et al. Reference Chakraborty, Liao, Adler and Leong2009) and DNA microarrays (Lee et al. Reference Lee, Cho, Huh, Ko, Lee, Jang, Lee, Kang and Choi2006), and to design new strategies for active flow control (Bushnell & McGinley Reference Bushnell and McGinley1989). Physically, EHD flow is characterized by a strong nonlinear interaction between the velocity field, the electric field and space charges: the electric force results in flow motion, which in turn affects the charge transport. The intricate nature of this nonlinearity defies a fundamental understanding of EHD flow. Moreover, as we will see, there still remains a mismatch or discrepancy between experimental observations and theoretical analysis.
One classic problem in EHD, named electroconvection, deals with the convective motions induced by unipolar charge injection into a dielectric liquid (of very low conductivity) that fills the gap between two parallel rigid plane electrodes. The Coulomb force acting on the free charge carriers tends to destabilize the system. Electroconvection is often compared to Rayleigh–Bénard convection (RBC) because of their similar geometry and convection patterns. Moreover, RBC is known to be analogous to the Taylor–Couette (TC) flow in the gap between two concentric rotating cylinders, where thermal energy transport in RBC corresponds to the transport of angular momentum in TC flow (Bradshaw Reference Bradshaw1969; Grossmann & Lohse Reference Grossmann and Lohse2000). In the linear regime of RBC, the flow is destabilized by the buoyancy force caused by the continued heating of the lower wall (an analogous role is played by centrifugal force in TC flow). As the thermal gradient exceeds a critical value, chaotic motion sets in. In EHD flow, the destabilizing factor is the electric force, acting in the wall-normal direction. However, the analogy between the two flows ends as soon as nonlinearities arise, especially where diffusive effects are concerned: in RBC, molecular diffusion constitutes the principal dissipative mechanism, whereas in EHD flow, it is the ion drift velocity
$KE$
(with
$K$
being the ionic mobility) that diffuses perturbations in the fluid. It is well known that RBC is of a supercritical nature, i.e. transition from the hydrostatic state to the finite-amplitude state occurs continuously as the controlling parameter, i.e. the Rayleigh number, is increased. For EHD flow, on the other hand, the bifurcation is subcritical, characterized by (i) an abrupt jump in motion amplitude from zero to a finite value, as a critical parameter is crossed, and (ii) the existence of a hysteresis loop. It is interesting to mention an analogy between EHD flow and polymeric flow: polymeric flow shows a hysteresis loop as well, as the first bifurcation is considered. In fact, the counterpart of EHD flow, i.e. magnetohydrodynamic (MHD) flow, has been compared to polymeric flow in Ogilvie & Proctor (Reference Ogilvie and Proctor2003).
Most studies in the EHD literature address electroconvection in the hydrostatic condition, i.e. without cross-flow. In this work we also investigate the EHD stability properties in the presence of cross-flow. Our interest is twofold. First, the potential of this flow configuration resides in the possibility of using the electric field to create large-scale rollers for flow manipulation; turbulent drag reduction designed in the spirit of Schoppa & Hussain (Reference Schoppa and Hussain1998) and investigated by Soldati & Banerjee (Reference Soldati and Banerjee1998) in the nonlinear regime is an example of this type. Secondly, EHD with cross-flow has been applied to wire–plate electrostatic precipitators, but, owing to the complex nature of the chaotic interaction between wall turbulence and the electric field, our current understanding of such flows is rather limited. Nonlinear EHD simulations with a cross-flow component have been reported in Soldati & Banerjee (Reference Soldati and Banerjee1998). More relevant to our linear problem is the unipolar-injection-induced instabilities in plane-parallel flows studied by Atten & Honda (Reference Atten and Honda1982) and Castellanos & Agrait (Reference Castellanos and Agrait1992). The former work focused on so-called electroviscous effects, defined by an increase of viscosity due to the applied electric field compared to the canonical channel flow. The latter work found that, at high Reynolds numbers, the destabilizing mechanism is linked to inertia, while, at sufficiently low Reynolds numbers, EHD instabilities are dominant. In this article, we will not only address the modal stability problem of EHD channel flow, as those two previous studies did, but also take into account transient effects, discussed briefly below, of the high-
$\mathit{Re}$
channel flow in the presence of an electric field. Our results should be interesting to those researching flow instability and transition to turbulence, especially for high-
$\mathit{Re}$
flow. The results will also shed light on the study of flow control strategies using EHD effects.
1.2. Stability of electrohydrodynamic flow
The endeavour to understand stability and transition to turbulence in EHD flow dates back to the 1970s, when Schneider & Watson (Reference Schneider and Watson1970) and Atten & Moreau (Reference Atten and Moreau1972), among the first, performed a linear stability analysis on the flow of dielectric liquids confined between two parallel electrodes with unipolar injection of charges. The mechanism for linear instability could be explained via the formation of an electric torque engendered by the convective motion when the driving electric force is sufficiently strong to overcome viscous diffusion. It was established in Atten & Moreau (Reference Atten and Moreau1972) that, in the weak-injection limit,
$C\ll 1$
, where
$C$
is the charge injection level, the flow is characterized by the criterion
$T_{c}C^{2}\approx 220.7$
, where
$T_{c}$
is the linear stability criterion for the stability parameter
$T$
, defined in the mathematical modelling section (§ 2.2). In the case of space-charge-limited (SCL) injection,
$C\rightarrow \infty$
, they found
$T_{c}\approx 160.75$
. However, according to Lacroix, Atten & Hopfinger (Reference Lacroix, Atten and Hopfinger1975) and Atten & Lacroix (Reference Atten and Lacroix1979), the experimentally determined stability criterion was notably different from the theoretical calculations. In the experiments performed by Atten & Lacroix (Reference Atten and Lacroix1979), the linear criterion was found to be
$T_{c}\approx 100$
in the case of SCL, which is far lower than the theoretically predicted value. It was argued then that this disagreement might be due to neglect of charge diffusion (Atten Reference Atten1976). We will address this discrepancy in the SCL case in this paper, and confirm that charge diffusion is indeed an important factor influencing the linear stability criterion in this case.
The first nonlinear stability analysis was performed by Félici (Reference Félici1971), who assumed a two-dimensional, a priori hydraulic model for the velocity field in the case of weak injection between two parallel plates. It was found that within the interval
$[T_{nl},T_{c}]$
, where
$T_{nl}$
is the nonlinear stability criterion for
$T$
, two solutions exist, namely, a stable state and an unstable finite-amplitude state. This finding corroborated the fact that the bifurcation in the unipolar injection problem is of a subcritical nature and that the flow has a hysteresis loop, as experimentally verified by Atten & Lacroix (Reference Atten and Lacroix1979). Physically, this subcritical bifurcation is related to the formation of a region of zero charge (Pérez & Castellanos Reference Pérez and Castellanos1989). Later, this simple hydraulic model was extended to three-dimensional hexagonal convective cells for the SCL case by Atten & Lacroix (Reference Atten and Lacroix1979), and it was shown that the most unstable hydrodynamic mode consists of hexagonal cells with the interior liquids flowing towards the injector. The nonlinear stability criterion for three-dimensional hexagonal cells, according to Atten & Lacroix (Reference Atten and Lacroix1979), was
$T_{nl}\approx 90$
in the experiments, but theoretical studies produced
$T_{nl}\approx 110$
.
Most of the previous linear stability analyses of EHD flow focus on the most unstable mode of the linear system, which is insufficient for a comprehensive flow analysis. In fact, theoretically, the linear stability analysis is linked to the characteristics of the linearized Navier–Stokes (NS) operator
$\boldsymbol{{\mathcal{L}}}$
, which, in the case of shear flows (in this paper, the cross-flow case), may be highly non-normal, i.e.
$\boldsymbol{{\mathcal{L}}}^{+}\boldsymbol{{\mathcal{L}}}\neq \boldsymbol{{\mathcal{L}}}\boldsymbol{{\mathcal{L}}}^{+}$
(with
$\boldsymbol{{\mathcal{L}}}^{+}$
denoting the adjoint of
$\boldsymbol{{\mathcal{L}}}$
) or, expressed differently, the eigenvectors of the linear operator are mutually non-orthogonal (see Trefethen et al.
Reference Trefethen, Trefethen, Reddy and Driscoll1993; Schmid & Henningson Reference Schmid and Henningson2001). For a normal operator (
$\boldsymbol{{\mathcal{L}}}^{+}\boldsymbol{{\mathcal{L}}}=\boldsymbol{{\mathcal{L}}}\boldsymbol{{\mathcal{L}}}^{+}$
), the dynamics of the perturbations is governed by the most unstable mode over the entire time horizon. In contrast, a non-normal operator has the potential for large transient amplification of the disturbance energy in the early linear phase, even though the most unstable mode is stable. The theory of non-modal stability analysis (Farrell & Ioannou Reference Farrell and Ioannou1996; Schmid Reference Schmid2007), the main tool to be used in this work, has been applied successfully to explain processes active during transition to turbulence in several shear flows. The fact that the bifurcation of EHD flow is subcritical, a trait often observed in shear flows governed by non-normal linearized operators, tempts one to think that the discrepancy between the experimental value
$T_{c}\approx 100$
and the theoretical value
$T_{c}\approx 161$
in the SCL regime of EHD flow might be examined in the light of non-modal stability theory. In fact, it seems surprising that this type of stability analysis has so far only rarely been applied to EHD flows, except for the work of Atten (Reference Atten1974) in the case of hydrostatic flow. The method we employ here is different from Atten’s quasi-stationary approach: non-modal stability theory, based on solving the initial-value problem, seeks the maximum disturbance energy growth over the entire time horizon when considering all admissible initial conditions and identifies the optimal initial condition for achieving this maximum energy growth. In Atten (Reference Atten1974), a quasi-stationary approach was taken that proposed that disturbances grow rapidly, compared to the time variation of the thickness of the unipolar layer; however, transient energy growth due to the non-normality of the linearized operator in hydrostatic EHD has been found to be rather limited in this work. This is in contrast with EHD Poiseuille flow, where non-normality is prevalent and should be considered from the outset.
The present paper extends the work by Martinelli, Quadrio & Schmid (Reference Martinelli, Quadrio and Schmid2011) and is organized as follows. In § 2, we present the mathematical model, the governing equations and the framework of the linear stability analysis. In § 3, numerical details are given; and a code validation is provided in the appendix. We then present in § 4 the results of the modal and non-modal stability analysis, and in § 5 the energy analysis. Finally, in § 6, we summarize our findings and conclude with a discussion.
2. Problem formulation
2.1. Mathematical modelling
We consider the planar geometry sketched in figure 1, where the Cartesian coordinate system used in this work is (
$x,y,z$
) or (
$\mathbf{1}_{x},\mathbf{1}_{y},\mathbf{1}_{z}$
) as the streamwise, wall-normal and spanwise directions, respectively. The two flat electrodes extend infinitely in the
$x$
and
$z$
directions, and the applied voltage only varies in the
$y$
direction. The distance between the two electrodes is
$2L^{\ast }$
. The dimensional variables and parameters are denoted with a superscript
$\text{}^{\ast }$
. The electric field satisfies the reduced Maxwell equations. The charges are generated through electrochemical reactions on the charge-injecting electrode (Alj et al.
Reference Alj, Denat, Gosse, Gosse and Nakamura1985). Since the electric conductivity is very low, conduction currents are negligible even in the presence of large electric fields. Therefore, magnetic effects in the Maxwell equations can be neglected (Melcher Reference Melcher1981; Castellanos Reference Castellanos1998), leading to the quasi-electrostatic limit of the Maxwell equations
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}^{\ast }\times \boldsymbol{E}^{\ast }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}^{\ast }\boldsymbol{\cdot }\boldsymbol{D}^{\ast }=Q^{\ast }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial Q^{\ast }}{\partial t^{\ast }}+\boldsymbol{{\rm\nabla}}^{\ast }\boldsymbol{\cdot }\boldsymbol{J}^{\ast }=0, & \displaystyle\end{eqnarray}$$
$\boldsymbol{E}^{\ast }$
is the electric field,
$\boldsymbol{D}^{\ast }={\it\epsilon}^{\ast }\boldsymbol{E}^{\ast }$
denotes the electric displacement,
${\it\epsilon}^{\ast }$
stands for the fluid permittivity, which we assume constant here,
$Q^{\ast }$
represents the charge density and
$\boldsymbol{J}^{\ast }$
is the current density. Considering (2.1a
), it is a well-known practice to define a potential field
${\it\phi}^{\ast }$
according to
$\boldsymbol{E}^{\ast }=-\boldsymbol{{\rm\nabla}}^{\ast }{\it\phi}^{\ast }$
. Combining the first two equations, (2.1a
) and (2.1b
), we can write the governing equation for
${\it\phi}^{\ast }$
as 
$$\begin{eqnarray}{\rm\nabla}^{\ast 2}{\it\phi}^{\ast }=-\frac{Q^{\ast }}{{\it\epsilon}^{\ast }}.\end{eqnarray}$$
Sketch of the electroconvection problem with coordinate system (
$x,y,z$
). In the non-hydrostatic case, a flow rate is induced along the streamwise (
$x$
) direction.
The current density
$\boldsymbol{J}^{\ast }$
arises from several sources. By modelling the EHD flow with only one ionic species in a perfectly insulating fluid (conductivity
${\it\sigma}^{\ast }=0$
), one can express
$\boldsymbol{J}^{\ast }$
as (Castellanos Reference Castellanos1998)
$$\begin{eqnarray}\boldsymbol{J}^{\ast }=K^{\ast }\boldsymbol{E}^{\ast }Q^{\ast }+\boldsymbol{U}^{\ast }Q^{\ast }-D_{{\it\nu}}^{\ast }\boldsymbol{{\rm\nabla}}^{\ast }Q^{\ast },\end{eqnarray}$$
where the first term accounts for the drift of ions (with respect to the fluid) under the effect of the electric field, moving at the relative velocity
$K^{\ast }\boldsymbol{E}^{\ast }$
, with
$K^{\ast }$
as the ionic mobility, the second term represents the convection of ions due to the fluid velocity
$\boldsymbol{U}^{\ast }$
, and the last term takes into account the charge diffusion, with
$D_{{\it\nu}}^{\ast }$
as the diffusion coefficient. Since the work of Pérez & Castellanos (Reference Pérez and Castellanos1989), the vast body of literature, with the exception of Kourmatzis & Shrimpton (Reference Kourmatzis and Shrimpton2012) for turbulent EHD flow, neglects the charge diffusion term because of its small value when compared to the drift terms. However, we will show that, even though the numerical value of
$D_{{\it\nu}}^{\ast }$
is very small, its impact on the flow dynamics is undeniable.
The flow field is incompressible, viscous and Newtonian and governed by the NS equations, which, in vector notation, read
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}^{\ast }\boldsymbol{\cdot }\boldsymbol{U}^{\ast }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}^{\ast }\frac{\partial \boldsymbol{U}^{\ast }}{\partial t^{\ast }}+{\it\rho}^{\ast }(\boldsymbol{U}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}^{\ast })\boldsymbol{U}^{\ast }=-\boldsymbol{{\rm\nabla}}^{\ast }P^{\ast }+{\it\mu}^{\ast }{\rm\nabla}^{\ast 2}\boldsymbol{U}^{\ast }+\boldsymbol{F}_{q}^{\ast }, & \displaystyle\end{eqnarray}$$
$\boldsymbol{U}^{\ast }$
is the velocity field,
$P^{\ast }$
the pressure,
${\it\rho}^{\ast }$
the density,
${\it\mu}^{\ast }={\it\rho}^{\ast }{\it\nu}^{\ast }$
the dynamic viscosity (
${\it\nu}^{\ast }$
the kinematic viscosity) and
$\boldsymbol{F}_{q}^{\ast }$
the volumetric density of electric force, which expresses the coupling between the fluid and the electric field. In general,
$\boldsymbol{F}_{q}^{\ast }$
can be written as 
$$\begin{eqnarray}\boldsymbol{F}_{q}^{\ast }=Q^{\ast }\boldsymbol{E}^{\ast }-\frac{1}{2}|\boldsymbol{E}^{\ast }|^{2}\boldsymbol{{\rm\nabla}}^{\ast }{\it\epsilon}^{\ast }+\boldsymbol{{\rm\nabla}}^{\ast }\left[\frac{|\boldsymbol{E}^{\ast }|^{2}}{2}{\it\rho}^{\ast }\frac{\partial {\it\epsilon}^{\ast }}{\partial {\it\rho}^{\ast }}\right],\end{eqnarray}$$
where the three terms on the right-hand side represent, respectively, the Coulomb force, the dielectric force and the electrostrictive force. The Coulomb force is commonly the strongest force when a DC voltage is applied. As we assume an isothermal and homogeneous fluid, the permittivity
${\it\epsilon}$
is constant in space. As a result, the dielectric force is zero (however, it would be dominant in the case of an AC voltage). The electrostrictive force can be incorporated into the pressure term of the NS equation, as it can be expressed as the gradient of a scalar field. Therefore, the only remaining term of interest in our formulation is the Coulomb force.
The system is supplemented by suitable boundary conditions. In our problem, we assume periodic boundary conditions in the wall-parallel directions. The no-slip and no-penetration conditions for the velocities are assumed at the channel walls. For the potential field, we require Dirichlet conditions on both walls, on the injector
${\it\phi}^{\ast }(L^{\ast })={\it\phi}_{0}^{\ast }$
and the collector
${\it\phi}^{\ast }(-L^{\ast })=0$
in order to fix the potential drop
${\rm\Delta}{\it\phi}_{0}^{\ast }$
between the electrodes. The injection mechanism is autonomous and homogeneous, meaning that the charge density is constant on the injector, not influenced by the nearby electric field and has a zero wall-normal flux of charge on the collector, i.e.
$Q^{\ast }(L^{\ast })=-Q_{0}^{\ast }$
and
$(\partial Q^{\ast }/\partial y^{\ast })(-L^{\ast })=0$
. Owing to the homogeneity in the wall-parallel directions, there is no requirement for boundary conditions in the
$x$
and
$z$
directions.
2.2. Non-dimensionalized governing equations
In the no-cross-flow case, as we are interested in the effect of the electric field on the flow dynamics, we non-dimensionalize the full system with the characteristics of the electric field, i.e.
$L^{\ast }$
(half-distance between the electrodes),
${\rm\Delta}{\it\phi}_{0}^{\ast }$
(voltage difference applied to the electrodes) and
$Q_{0}^{\ast }$
(injected charge density). Accordingly, the time
$t^{\ast }$
is non-dimensionalized by
$L^{\ast 2}/(K^{\ast }{\rm\Delta}{\it\phi}_{0}^{\ast })$
, the velocity
$\boldsymbol{U}^{\ast }$
by
$K^{\ast }{\rm\Delta}{\it\phi}_{0}^{\ast }/L^{\ast }$
, the pressure
$P^{\ast }$
by
${\it\rho}_{0}^{\ast }K^{\ast 2}{\rm\Delta}{\it\phi}_{0}^{\ast 2}/L^{\ast 2}$
, the electric field
$\boldsymbol{E}^{\ast }$
by
${\rm\Delta}{\it\phi}_{0}^{\ast }/L^{\ast }$
and the electric density
$Q^{\ast }$
by
$Q_{0}^{\ast }$
. Therefore, the non-dimensional equations read
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{U}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{U}}{\partial t}+(\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{U}=-\boldsymbol{{\rm\nabla}}P+\frac{M^{2}}{T}{\rm\nabla}^{2}\boldsymbol{U}+CM^{2}Q\boldsymbol{E}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial Q}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[(\boldsymbol{E}+\boldsymbol{U})Q]=\frac{1}{\mathit{Fe}}{\rm\nabla}^{2}Q, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\rm\nabla}^{2}{\it\phi}=-CQ, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{E}=-\boldsymbol{{\rm\nabla}}{\it\phi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}M=\frac{({\it\epsilon}^{\ast }/{\it\rho}_{0}^{\ast })^{1/2}}{K^{\ast }},\quad T=\frac{{\it\epsilon}^{\ast }{\rm\Delta}{\it\phi}_{0}^{\ast }}{K^{\ast }{\it\mu}^{\ast }},\quad C=\frac{Q_{0}^{\ast }L^{\ast 2}}{{\rm\Delta}{\it\phi}_{0}^{\ast }{\it\epsilon}^{\ast }},\quad \mathit{Fe}=\frac{K^{\ast }{\rm\Delta}{\it\phi}_{0}^{\ast }}{D_{{\it\nu}}^{\ast }}.\end{eqnarray}$$
Additionally, the non-dimensional boundary conditions are
$\boldsymbol{U}(\pm 1)=0$
,
${\it\phi}(1)=1$
,
${\it\phi}(-1)=0$
,
$Q(1)=-1$
and
$(\partial Q/\partial y)(-1)=0$
.
Various dimensionless groups appear in the equations as written above. First,
$M$
is the ratio between the hydrodynamic mobility
$({\it\epsilon}^{\ast }/{\it\rho}_{0}^{\ast })^{1/2}$
and the true ion mobility
$K^{\ast }$
. Gases usually take on a value of
$M$
less than 0.1 and liquids have values of
$M$
greater than 1 (Castellanos & Agrait Reference Castellanos and Agrait1992). Second, Taylor’s parameter
$T$
represents the ratio of the Coulomb force to the viscous force. It is the principal stability parameter, assuming a similar role as the Rayleigh number in RBC. Next,
$C$
measures the injection level. When
$C\gg 1$
, the system is in a strong-injection regime, and when
$C\ll 1$
, it is in a weak-injection regime. Finally,
$\mathit{Fe}$
is the reciprocal of the charge diffusivity coefficient. The factor
$M^{2}/T$
appearing in (2.6b
) can be interpreted as the ratio between the charge relaxation time
$L^{\ast 2}/(K^{\ast }{\rm\Delta}{\it\phi}_{0}^{\ast })$
by drift and the momentum relaxation time
$L^{\ast 2}/{\it\nu}^{\ast }$
. This mathematical model for EHD flow has been assumed and studied in many previous investigations of linear stability and turbulence analyses for a dielectric liquid subject to unipolar injection of ions (Lacroix et al.
Reference Lacroix, Atten and Hopfinger1975; Traoré & Pérez Reference Traoré and Pérez2012; Wu et al.
Reference Wu, Traoré, Vázquez and Pérez2013), except that the diffusion term in (2.6c
) is usually neglected (excluding the study of Kourmatzis & Shrimpton Reference Kourmatzis and Shrimpton2012).
2.3. Linear stability problem
The linear problem is obtained by decomposing the flow variable as a sum of base state and perturbation, i.e.
$\boldsymbol{U}=\bar{\boldsymbol{U}}+\boldsymbol{u}$
,
$P=\bar{P}+p$
,
$\boldsymbol{E}=\bar{\boldsymbol{E}}+\boldsymbol{e}$
,
$\boldsymbol{D}=\bar{\boldsymbol{D}}+\boldsymbol{d}$
,
$Q=\bar{Q}+q$
and
${\it\phi}=\bar{{\it\phi}}+{\it\varphi}$
. For the vector fields, we have
$\boldsymbol{u}=(u,v,w)$
and
$\boldsymbol{e}=(e_{1},e_{2},e_{3})$
along the three Cartesian coordinate directions. After substituting the decompositions into the governing equations (2.6a–e
), subtracting from them the governing equations for the base states and retaining the terms of first order, the linear system reads
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\bar{\boldsymbol{U}}+(\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p+\frac{M^{2}}{T}{\rm\nabla}^{2}\boldsymbol{u}+CM^{2}(q\bar{\boldsymbol{E}}+\bar{Q}\boldsymbol{e}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial q}{\partial t}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[(\bar{\boldsymbol{E}}+\bar{\boldsymbol{U}})q+(\boldsymbol{e}+\boldsymbol{u})\bar{Q}]=\frac{1}{\mathit{Fe}}{\rm\nabla}^{2}q, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\rm\nabla}^{2}{\it\varphi}=-Cq, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{e}=-\boldsymbol{{\rm\nabla}}{\it\varphi}, & \displaystyle\end{eqnarray}$$
$\boldsymbol{u}(\pm 1)=0$
,
${\it\varphi}(\pm 1)=0$
,
$q(1)=0$
and
$(\partial q/\partial y)(-1)=0$
.
2.3.1. Base states
The base states are the solutions to (2.6a–e
) in the case of no time dependence. Owing to the periodicity in the wall-parallel directions, we can reduce the shape of the base states as functions of
$y$
only, that is,
$\bar{\boldsymbol{U}}=\bar{U}(y)\mathbf{1}_{y}$
and
$\bar{\boldsymbol{E}}=\bar{E}(y)\mathbf{1}_{y}$
. For the base flow
$\bar{U}(y)$
, we are interested in the hydrostatic and pressure-driven Poiseuille flows, which, after non-dimensionalization, are given by
$$\begin{eqnarray}\bar{U}(y)=0,\quad \bar{U}(y)=\mathit{Re}\frac{M^{2}}{T}(1-y^{2})=(1-y^{2}),\end{eqnarray}$$
respectively, in which the (electric) Reynolds number is defined as
$\mathit{Re}=T/M^{2}=K^{\ast }{\rm\Delta}{\it\phi}_{0}^{\ast }/{\it\nu}^{\ast }$
(in order to enforce the same constant flow rate). It is a passive parameter in the hydrostatic case, but becomes a free parameter in the presence of high-
$\mathit{Re}$
cross-flow, in which, consequently,
$M$
would be the passive parameter. Therefore, in the Poiseuille flow case, we modify the governing equation (2.8b
) by substituting the relation
$\mathit{Re}=T/M^{2}$
to obtain
$$\begin{eqnarray}\frac{\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\bar{\boldsymbol{U}}+(\bar{\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p+\frac{1}{\mathit{Re}}{\rm\nabla}^{2}\boldsymbol{u}+\frac{CT}{\mathit{Re}}(q\bar{\boldsymbol{E}}+\bar{Q}\boldsymbol{e}).\end{eqnarray}$$
By doing so, it is more obvious to identify the effects of
$T$
and
$C$
on the electric force term. The parameter
$\mathit{Re}=K^{\ast }{\rm\Delta}{\it\phi}_{0}^{\ast }/{\it\nu}^{\ast }$
here coincides with the canonical hydrodynamic equivalent
$\mathit{Re}_{h}=U^{\ast }L^{\ast }/{\it\nu}^{\ast }$
because of the electric scaling we chose. However, in a general sense, the two may not necessarily be identical. The non-dimensional quantity
$$\begin{eqnarray}\frac{\mathit{Re}}{\mathit{Re}_{h}}=\frac{K^{\ast }{\rm\Delta}{\it\phi}_{0}^{\ast }}{U^{\ast }L^{\ast }}=\frac{L^{\ast }/U^{\ast }}{L^{\ast 2}/(K^{\ast }{\rm\Delta}{\it\phi}_{0}^{\ast })}\end{eqnarray}$$
relates the eddy turnover time and the charge relaxation time by the drift. According to the equality
$T=\mathit{Re}\,M^{2}$
, when
$\mathit{Re}$
is near linear criticality at 5772 and
$T$
is around
$10^{2}$
,
$M\approx 0.1$
. This implies that the working liquid is gas. Moreover, in contrast to the nonlinear constitutive modelling for polymers in viscoelastic flow, the base flow is not modified under the influence of the base electric field, even though the coupling between
$\boldsymbol{U}$
and
$Q$
is nonlinear in (2.6c
). This is because the directions of the base flow and the base electric field are perpendicular. Nevertheless, the base pressure gradient in the wall-normal direction is no longer zero.
The base electric field
$\bar{\boldsymbol{E}}(y)$
can be solved from (2.6c–e
), recast into an equation for
$\bar{{\it\phi}}$
only, which reads
$$\begin{eqnarray}\bar{{\it\phi}}^{\prime }\bar{{\it\phi}}^{\prime \prime \prime }+(\bar{{\it\phi}}^{\prime \prime })^{2}+\frac{1}{\mathit{Fe}}\bar{{\it\phi}}^{\prime \prime \prime \prime }=0,\end{eqnarray}$$
where prime
$\text{}^{\prime }$
denotes the spatial derivative with respect to the
$y$
direction. The boundary conditions are
$\bar{{\it\phi}}(1)=1$
,
$\bar{{\it\phi}}(-1)=0$
,
$\bar{{\it\phi}}^{\prime \prime }(1)=-C$
and
$\bar{{\it\phi}}^{\prime \prime \prime }(-1)=0$
. Analytical solutions to this fourth-order ordinary differential equation can be obtained by observing that the equation can be transformed into a Riccati equation; alternatively, as we do here, a simple numerical integration combined with a nonlinear gradient method provides us with the required
$\bar{{\it\phi}}(y)$
profile. The Poiseuille base flow and the base states of the electric and charge fields are shown in figure 2.
The base states: (a)
$\bar{U}$
, (b)
$\bar{{\it\phi}}$
and (c)
$\bar{Q}$
.
2.3.2. Matrix representation
In linear stability analysis, it is a common practice to rewrite the fluid system (2.8a,b
) in terms of the wall-normal velocity
$v$
and the wall-normal vorticity
${\it\eta}=\partial _{z}u-\partial _{x}w$
by eliminating the pressure term. For the electric field, the three equations (2.8c–e
) can be reduced to one for
${\it\varphi}$
. Therefore, the governing equations (2.8a–e
) become, in terms of a
$v{-}{\it\eta}{-}{\it\varphi}$
formulation,
$$\begin{eqnarray}\displaystyle \frac{\partial {\rm\nabla}^{2}v}{\partial t} & = & \displaystyle \left[-\bar{U}\frac{\partial }{\partial x}{\rm\nabla}^{2}+\bar{U}^{\prime \prime }\frac{\partial }{\partial x}+\frac{M^{2}}{T}{\rm\nabla}^{4}\right]\!v\nonumber\\ \displaystyle & & \displaystyle +\,M^{2}\left[-\bar{{\it\phi}}^{\prime \prime \prime }\left({\rm\nabla}^{2}-\frac{\partial ^{2}}{\partial y^{2}}\right){\it\varphi}+\bar{{\it\phi}}^{\prime }\left({\rm\nabla}^{2}-\frac{\partial ^{2}}{\partial y^{2}}\right){\rm\nabla}^{2}{\it\varphi}\right]\!,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\partial {\it\eta}}{\partial t} & = & \displaystyle -\bar{U}\frac{\partial }{\partial x}{\it\eta}-\bar{U}^{\prime }\frac{\partial v}{\partial z}+\frac{M^{2}}{T}{\rm\nabla}^{2}{\it\eta},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\partial {\rm\nabla}^{2}{\it\varphi}}{\partial t} & = & \displaystyle \bar{{\it\phi}}^{\prime }\frac{\partial {\rm\nabla}^{2}{\it\varphi}}{\partial y}+\bar{{\it\phi}}^{\prime \prime \prime }\frac{\partial {\it\varphi}}{\partial y}+2\bar{{\it\phi}}^{\prime \prime }{\rm\nabla}^{2}{\it\varphi}-\bar{U}\frac{\partial {\rm\nabla}^{2}{\it\varphi}}{\partial x}-\bar{{\it\phi}}^{\prime \prime \prime }v+\frac{1}{\mathit{Fe}}{\rm\nabla}^{4}{\it\varphi},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}v(\pm 1)=0,\quad v^{\prime }(\pm 1)=0,\\ {\it\eta}(\pm 1)=0,\\ {\it\varphi}(\pm 1)=0,\quad {\it\varphi}^{\prime \prime }(1)=0,\quad {\it\varphi}^{\prime \prime \prime }(-1)=0.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
For compactness, we write
${\bf\gamma}=(v,{\it\eta},{\it\varphi})^{\text{T}}$
, and the linearized system, recast in matrix notation, becomes
$$\begin{eqnarray}\left(\begin{array}{@{}ccc@{}}{\rm\nabla}^{2} & \pmb{{\it 0}} & \pmb{{\it 0}}\\ \pmb{{\it 0}} & \unicode[STIX]{x1D644} & \pmb{{\it 0}}\\ \pmb{{\it 0}} & \pmb{{\it 0}} & {\rm\nabla}^{2}\end{array}\right)\frac{\partial }{\partial t}\left(\begin{array}{@{}c@{}}v\\ {\it\eta}\\ {\it\varphi}\end{array}\right)=\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D647}_{os} & \pmb{{\it 0}} & \unicode[STIX]{x1D647}_{v{\it\varphi}}\\ \unicode[STIX]{x1D647}_{c} & \unicode[STIX]{x1D647}_{sq} & \pmb{{\it 0}}\\ \unicode[STIX]{x1D647}_{{\it\varphi}v} & \pmb{{\it 0}} & \unicode[STIX]{x1D647}_{{\it\varphi}{\it\varphi}}\end{array}\right)\left(\begin{array}{@{}c@{}}v\\ {\it\eta}\\ {\it\varphi}\end{array}\right),\end{eqnarray}$$
where
$\unicode[STIX]{x1D644}$
denotes the identity matrix and the submatrices
$\unicode[STIX]{x1D647}_{os}$
,
$\unicode[STIX]{x1D647}_{v{\it\varphi}}$
,
$\unicode[STIX]{x1D647}_{c}$
,
$\unicode[STIX]{x1D647}_{sq}$
,
$\unicode[STIX]{x1D647}_{{\it\varphi}v}$
and
$\unicode[STIX]{x1D647}_{{\it\varphi}{\it\varphi}}$
can be easily deduced from (2.13a–c
). To represent the system even more compactly, we can rewrite the linearized problem (2.15) as
$$\begin{eqnarray}\boldsymbol{{\mathcal{A}}}\frac{\partial {\bf\gamma}}{\partial t}=\boldsymbol{{\mathcal{B}}}{\bf\gamma}\quad \Longrightarrow \quad \frac{\partial {\bf\gamma}}{\partial t}=\boldsymbol{{\mathcal{L}}}{\bf\gamma},\end{eqnarray}$$
where
$\boldsymbol{{\mathcal{L}}}=\boldsymbol{{\mathcal{A}}}^{-1}\boldsymbol{{\mathcal{B}}}$
represents the linearized NS operator for EHD flow.
Since the flow is homogeneous in the wall-parallel directions, the perturbations are assumed to take on a wave-like shape. Moreover, as we consider a linear problem with a steady base flow, it is legitimate to examine the frequency response of the linear system for each frequency individually. These two simplifications lead to
$$\begin{eqnarray}\boldsymbol{f}(x,y,z,t)=\hat{\boldsymbol{f}}(y,t)\exp (\text{i}{\it\alpha}x+\text{i}{\it\beta}z)=\tilde{\boldsymbol{f}}(y)\exp (-\text{i}{\it\omega}t)\exp (\text{i}{\it\alpha}x+\text{i}{\it\beta}z),\end{eqnarray}$$
where
$\boldsymbol{f}$
could represent any flow variable in
$(\boldsymbol{u},p,\boldsymbol{e},q,{\it\varphi})^{\text{T}}$
,
$\hat{\boldsymbol{f}}(y,t)$
and
$\tilde{\boldsymbol{f}}(y)$
are the shape functions,
${\it\alpha}$
and
${\it\beta}$
are the real-valued streamwise and spanwise wavenumbers, and the complex-valued
${\it\omega}$
is the circular frequency of the perturbation, with its real part
${\it\omega}_{r}$
representing the phase speed and its imaginary part
${\it\omega}_{i}$
representing the growth rate of the linear perturbation. Upon substitution of the above expression into the linear problem (2.16), we arrive at an eigenvalue problem for the
$v{-}{\it\eta}{-}{\it\varphi}$
formulation that reads
$$\begin{eqnarray}-\!\text{i}{\it\omega}\tilde{{\bf\gamma}}=\boldsymbol{{\mathcal{L}}}\tilde{{\bf\gamma}},\end{eqnarray}$$
where
$-\text{i}{\it\omega}$
is the eigenvalue and
$\tilde{{\bf\gamma}}$
is the corresponding eigenvector. Both formulations, (2.16) and (2.18), would be relevant as discussed in a recent review by Schmid & Brandt (Reference Schmid and Brandt2014). The least unstable eigenvalues obtained from the eigenproblem formulation (2.18) would determine the asymptotic behaviour of the linear system, while the initial-value problem formulation (2.16) could be used to examine the dynamics of the fluid system evolving over a finite time scale.
2.3.3. Energy norm
In our calculation of the non-modal transient growth, we define the total energy density of the perturbation contained in a control volume
${\it\Omega}$
as
$$\begin{eqnarray}\displaystyle \int _{{\it\Omega}}{\mathcal{E}}^{\ast }\text{d}V^{\ast } & = & \displaystyle \int _{{\it\Omega}}({\mathcal{E}}_{k}^{\ast }+{\mathcal{E}}_{{\it\varphi}}^{\ast })\text{d}V^{\ast }=\int _{{\it\Omega}}\frac{1}{2}({\it\rho}_{0}^{\ast }\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\boldsymbol{u}^{\ast }+\boldsymbol{e}^{\ast }\boldsymbol{\cdot }\boldsymbol{d}^{\ast })\text{d}V^{\ast }\nonumber\\ \displaystyle & = & \displaystyle \int _{{\it\Omega}}\frac{1}{2}({\it\rho}_{0}^{\ast }(u^{\ast 2}+v^{\ast 2}+w^{\ast 2})+{\it\epsilon}^{\ast }|\boldsymbol{{\rm\nabla}}^{\ast }{\it\varphi}^{\ast }|^{2})\text{d}V^{\ast }.\end{eqnarray}$$
The perturbed electric energy
${\mathcal{E}}_{{\it\varphi}}^{\ast }$
follows the definition in Castellanos (Reference Castellanos1998). In terms of the
$v{-}{\it\eta}{-}{\it\varphi}$
formulation, the non-dimensionalized energy norm in spectral space becomes
$$\begin{eqnarray}\displaystyle \int _{{\it\Omega}}{\mathcal{E}}\text{d}V & = & \displaystyle \frac{1}{2}\cdot \frac{1}{2}\int \hat{{\bf\gamma}}^{\dagger }\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D644}+\displaystyle \frac{1}{k^{2}}\unicode[STIX]{x1D63F}_{1}^{\dagger }\unicode[STIX]{x1D63F}_{1} & \pmb{{\it 0}} & \pmb{{\it 0}}\\ \pmb{{\it 0}} & \displaystyle \frac{1}{k^{2}}\unicode[STIX]{x1D644} & \pmb{{\it 0}}\\ \pmb{{\it 0}} & \pmb{{\it 0}} & M^{2}(k^{2}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D63F}_{1}^{\dagger }\unicode[STIX]{x1D63F}_{1})\end{array}\right)\hat{{\bf\gamma}}\text{d}y\nonumber\\ \displaystyle & = & \displaystyle \int _{{\it\Omega}}\hat{{\bf\gamma}}^{\dagger }\unicode[STIX]{x1D648}\hat{{\bf\gamma}}\text{d}y,\end{eqnarray}$$
where the superscript dagger
$\text{}^{\dagger }$
denotes the complex conjugate,
$k^{2}={\it\alpha}^{2}+{\it\beta}^{2}$
, and
$\unicode[STIX]{x1D63F}_{1}$
represents the first-derivative matrix with respect to the wall-normal direction (likewise for
$\unicode[STIX]{x1D63F}_{2}$
and
$\unicode[STIX]{x1D63F}_{3}$
below). The positive definite matrix
$\unicode[STIX]{x1D648}$
allows us to work in the
$L_{2}$
-norm. To do so, we apply a Cholesky decomposition to the weight matrix according to
$\unicode[STIX]{x1D648}=\unicode[STIX]{x1D641}^{\dagger }\unicode[STIX]{x1D641}$
and define
$\hat{{\bf\xi}}=\unicode[STIX]{x1D641}\hat{{\bf\gamma}}$
to arrive at
$$\begin{eqnarray}\int _{{\it\Omega}}{\mathcal{E}}\text{d}V=\int _{{\it\Omega}}\hat{{\bf\gamma}}^{\dagger }\unicode[STIX]{x1D648}\hat{{\bf\gamma}}\text{d}y=\int _{{\it\Omega}}\hat{{\bf\gamma}}^{\dagger }\unicode[STIX]{x1D641}^{\dagger }\unicode[STIX]{x1D641}\hat{{\bf\gamma}}\text{d}y=\int _{{\it\Omega}}\hat{{\bf\xi}}^{\dagger }\hat{{\bf\xi}}\text{d}y=\Vert \hat{{\bf\xi}}\Vert _{2},\end{eqnarray}$$
where
$\Vert \cdot \Vert _{2}$
represents the
$L_{2}$
-norm and, accordingly, the eigenvalue problem (2.18) becomes
$$\begin{eqnarray}-\!\text{i}{\it\omega}(\unicode[STIX]{x1D641}\hat{{\bf\gamma}})=\unicode[STIX]{x1D641}\boldsymbol{{\mathcal{L}}}\unicode[STIX]{x1D641}^{-1}(\unicode[STIX]{x1D641}\hat{{\bf\gamma}}).\end{eqnarray}$$
Therefore, once the linear operator is redefined as
$\boldsymbol{{\mathcal{L}}}_{L_{2}}=\unicode[STIX]{x1D641}\boldsymbol{{\mathcal{L}}}\unicode[STIX]{x1D641}^{-1}$
, we can conveniently use the
$L_{2}$
-norm and its associated inner product for all computations. The transient growth
$G$
, defined as the maximum energy growth over all possible initial conditions
$\hat{{\bf\xi}}_{0}$
, is given below in the
$L_{2}$
-norm,
$$\begin{eqnarray}G(t)=\max _{\hat{{\bf\xi}}_{0}}\frac{\Vert \hat{{\bf\xi}}(t)\Vert _{2}}{\Vert \hat{{\bf\xi}}(0)\Vert _{2}}=\max _{\hat{{\bf\xi}}_{0}}\frac{\Vert \boldsymbol{{\mathcal{T}}}\hat{{\bf\xi}}(0)\Vert _{2}}{\Vert \hat{{\bf\xi}}(0)\Vert _{2}}=\Vert \boldsymbol{{\mathcal{T}}}\Vert _{2}=\Vert \text{e}^{t\unicode[STIX]{x1D641}\boldsymbol{{\mathcal{L}}}\unicode[STIX]{x1D641}^{-1}}\Vert _{2},\end{eqnarray}$$
where
$\boldsymbol{{\mathcal{T}}}$
is the linear evolution operator, i.e. the solution to (2.16).
The parameters that are to be investigated include the injection level
$C$
, the mobility parameter
$M$
, the charge diffusion coefficient
$\mathit{Fe}$
, the Taylor parameter
$T$
, the Reynolds number
$\mathit{Re}$
, and the streamwise and spanwise wavenumbers
${\it\alpha}$
and
${\it\beta}$
.
3. Numerical method and validation
3.1. Numerical method
To discretize the eigenvalue problem (2.18), we use a spectral method based on collocation points chosen as the roots of Chebyshev polynomials. The MATLAB suite for partial differential equations by Weideman & Reddy (Reference Weideman and Reddy2000) is used for differentiation and integration.
To impose the boundary condition, we employ the boundary boarding technique (Boyd Reference Boyd2001), in which selected rows of the linear matrices are replaced directly by the boundary conditions. When solving the eigenvalue problem via the MATLAB routine eig with the above boundary condition enforced, we find that the eigenvalues converge for a sufficient number
$N$
of collocation points (see figure 12 and table 4 in the validation section in appendix A) and approach the pure hydrodynamic results as electric effects become negligible (see figure 13 and table 6). The corresponding eigenvectors, however, are incorrect, since they do not satisfy the proper boundary conditions (not shown). To overcome this difficulty, we employ an iterative technique to obtain the eigenvector associated with a specified eigenvalue. In the generalized eigenvalue problem (2.16), a desired eigenvalue
${\it\omega}$
(and its corresponding eigenvector) is targeted by applying the spectral transformation
$$\begin{eqnarray}\boldsymbol{{\mathcal{S}}}=(\boldsymbol{{\mathcal{B}}}-{\it\omega}\boldsymbol{{\mathcal{A}}}_{1})\setminus \boldsymbol{{\mathcal{A}}}_{1},\end{eqnarray}$$
where
$\boldsymbol{{\mathcal{A}}}_{1}=-\text{i}\boldsymbol{{\mathcal{A}}}$
and
$\boldsymbol{{\mathcal{S}}}$
will be processed by an iterative routine (Saad Reference Saad2011).
3.2. Validation
The stability problem for EHD flow is exceedingly challenging from a numerical point of view, which warrants a careful and thorough validation step, before results about stability characteristics, modal and non-modal solutions and physical mechanisms are produced. To conserve the clarity of the paper’s structure, we postpone the validation steps to the appendix A.
4. Results of stability analysis
4.1. Electrohydrodynamics without cross-flow
As mentioned earlier, the parameter
$T$
plays the main role in determining the flow instability. The critical
$T_{c}$
denotes the minimum value of
$T$
within the linear regime, above which infinitesimal disturbances can grow exponentially in time;
$T_{c}$
will vary with the flow parameters. In the case of no cross-flow, the effects of
$\mathit{Fe}$
,
$T$
,
$M$
and
$C$
on the flow stability are investigated. As has already been assumed, the flow will be confined to the SCL regime, implying a large value for
$C$
.
We display the neutral stability curve in figure 3 for different
$\mathit{Fe}$
at
$C=50$
,
$M=100$
,
$T=155$
,
${\it\alpha}=2.5$
and
${\it\beta}=0$
. In the case without cross-flow, one does not need to distinguish between the
$x$
and
$z$
axes, since neither is preferred by the base flow
$\bar{\boldsymbol{U}}=0$
; thus, we simply set
${\it\beta}=0$
. As mentioned in the validation section, results for
$\mathit{Fe}=10^{7}$
are very close to previous investigations. Even though the diffusion coefficient is small, it plays an important role in determining the critical
$T_{c}$
, as shown in figure 3(a,b). For example, for
$\mathit{Fe}=10^{3}$
the critical
$T_{c}$
declines to 140. In fact, the value of
$\mathit{Fe}$
could fall within the range
$10^{3}{-}10^{4}$
for real liquids (Pérez & Castellanos Reference Pérez and Castellanos1989), when
$\mathit{Fe}$
is non-dimensionalized in the same way as presented here. Physically, the effect of diffusion will smooth out sharp gradients in the flow. Unlike the unidirectional electric field pointing in the wall-normal direction, the diffusion effects act equally in all directions. With charge diffusion considered in the model, the discontinuous separatrix is blurred in the nonlinear phase (Pérez & Castellanos Reference Pérez and Castellanos1989). The physical mechanism of how charge diffusion influences the critical stability parameter
$T_{c}$
will be discussed by using an energy analysis (see § 5.1). In addition, the transient growth of disturbance energy has been discussed in Atten (Reference Atten1974) using the quasi-stationary method; transient energy growth has been confirmed as a minor factor in that work. This is also confirmed in our computations, as presented in figure 3(c): specifically, the figure shows that disturbance energy growth
$G$
reaches a value of approximately 3 at
$T=155$
for stable flows (
$\mathit{Fe}>10^{3}$
).
Effect of
$\mathit{Fe}$
: (a) neutral stability curves for various
$\mathit{Fe}$
: (b)
$T_{c}$
as a function of
$\mathit{Fe}$
; and (c) transient energy growth versus time. The parameters are
$C=50$
,
$M=100$
,
$T=155$
,
${\it\alpha}=2.5$
,
${\it\beta}=0$
and
$N=250$
. The direction of the arrow indicates increasing
$\mathit{Fe}$
.
The role of
$M$
in EHD is analogous to that of the Prandtl number in RBC. In figure 4(a), it is shown that the variation of
$M$
exerts no influence on the linear stability criterion,
$T_{c}=159.58$
at
$C=100$
,
$\mathit{Fe}=10^{5}$
,
${\it\alpha}=2.57$
and
${\it\beta}=0$
; the same finding has been reported in Atten & Moreau (Reference Atten and Moreau1972). For the transient dynamics, however, the same conclusion does not hold, as evidenced in figure 4(b). The plot describes a trend of increasing
$G_{max}$
with smaller
$M$
. The slopes at the final time are slightly different for each
$M$
, indicating that the asymptotic growth rates differ slightly (while the linear stability criterion remains the same).
Effect of
$M$
: (a) neutral stability curve; and (b) transient energy growth versus time for different
$M$
and
$T=155$
. The parameters are
$C=100$
,
$\mathit{Fe}=10^{5}$
,
${\it\alpha}=2.57$
,
${\it\beta}=0$
and
$N=250$
without cross-flow. The direction of the arrow indicates increasing
$M$
.
Figure 5 depicts the influence on
$C$
, which measures the intensity of charge injection. Atten & Moreau (Reference Atten and Moreau1972) and Atten & Lacroix (Reference Atten and Lacroix1979) reported a dependence of the critical value
$T_{c}$
on the parameter
$C$
. In figure 5(a), we see that, in the SCL regime, increasing
$C$
will yield lower
$T_{c}$
. This result can be understood from a physical argument. Increasing the intensity of charge injection will lead to a higher concentration of charges between the electrodes. The linear instability mechanism, as discussed above, relies on the formation of an electric torque due to convective motions. With higher charge concentration, the electric torque is stronger. Therefore, a lower voltage difference is required, which amounts to stating that a lower
$T$
will be sufficient to form an electric torque of comparable strength. But as we are in the SCL regime (with a value of
$C=50$
considered very large), a rise of
$C$
to 200 only yields a minor decrease in
$T_{c}$
. In contrast, the transient dynamics of the perturbation energy
$G$
appears not to be influenced by a change in
$C$
for early times; for example, see the time interval
$t\in [0,3]$
in figure 5(b).
Effect of
$C$
: (a) neutral stability curve for four different injection levels
$C$
within the SCL regime; and (b) transient energy growth for different
$C$
at
$T=155$
and
${\it\alpha}=2.57$
. The parameters are
$M=100$
,
$\mathit{Fe}=10^{4}$
,
${\it\beta}=0$
and
$N=250$
, without cross-flow. The direction of the arrow indicates increasing
$C$
.
4.2. Electrohydrodynamics with cross-flow
When cross-flow is considered, the property of the linearized system changes due to the presence of a base shear in the flow. Especially, this shear will render the linearized operator ‘more non-normal’. We first note that, in the modal stability analysis, Squire’s theorem still holds for EHD Poiseuille flow, that is, a two-dimensional instability will be encountered first. This can be easily verified by a perfect analogy with standard viscous theory (Schmid & Henningson Reference Schmid and Henningson2001). Moreover, there are two sets of scales in the EHD problem with cross-flow. To study the influence of the cross-flow on the electric and the charge fields, the values of
$M$
and
$T$
are kept in the vicinity of the values in the previous section: the scale of the electric field will be considered primarily, whereas, when we examine the effects of the electric field exerted on the canonical Poiseuille flow, we take the value of the free parameter
$\mathit{Re}$
around 5772, i.e. the linear stability criterion for pressure-driven flow; the latter choice introduces a scale based on the hydrodynamics. In both cases, we will enforce the relation
$\mathit{Re}=T/M^{2}$
, which results in the Reynolds number
$\mathit{Re}$
being rather low in the former case (denoted as the low-
$\mathit{Re}$
case) and relatively high in the latter case (referred to as the high-
$\mathit{Re}$
case).
4.2.1. Electrohydrodynamics: low
$\mathit{Re}$
We have demonstrated that non-modal effects in hydrostatic EHD are not significant. In the presence of cross-flow, given that the Reynolds number in this section is considered small, we expect the non-normality to be rather moderate as well. For this reason, we will mainly focus on the modal stability characteristics for the low-
$\mathit{Re}$
case.
In figure 6(a), it is observed that the symmetry of the hydrostatic EHD spectrum is now broken due to the presence of cross-flow. The most unstable perturbation travels at a positive phase speed
$u_{p}={\it\omega}_{r}/{\it\alpha}=2.256/2.57=0.8778$
, induced by cross-flow convection (the centreline velocity of the cross-flow is 1, as we set
$\mathit{Re}=T/M^{2}$
in (2.9)). In figure 6(b), we show the neutral stability curve for
$C=50$
,
$M=100$
,
${\it\beta}=0$
and
$N=250$
, which can be directly compared to the results in figure 3. Note that, since
$\mathit{Re}=T/M^{2}$
is enforced, the Reynolds number
$\mathit{Re}$
is not identical for each point, but generally small. We see that, with cross-flow, the critical
$T_{c}$
decreases compared to the no-cross-flow case; this indicates that the flow is more unstable in the presence of a low-
$\mathit{Re}$
cross-flow compared to the results in figure 3(a). To investigate the reason behind this destabilization, we again resort to an energy analysis in § 5.1.2. Previously, an energy analysis for EHD with cross-flow has been studied in Castellanos & Agrait (Reference Castellanos and Agrait1992).
(a) Spectra for the case with and without cross-flow for
$C=50$
,
$M=100$
,
$\mathit{Fe}=10^{5}$
,
$T=160$
,
$\mathit{Re}=T/M^{2}=0.016$
(for cross-flow),
${\it\alpha}=2.57$
,
${\it\beta}=0$
and
$N=250$
. (b) Effect of
$\mathit{Fe}$
on the neutral stability curve. The parameters are
$C=50$
,
$M=100$
,
${\it\beta}=0$
and
$N=250$
, with cross-flow. The direction of the arrow indicates increasing
$\mathit{Fe}$
.
Even though varying
$M$
has no effect on the linear stability when
$\bar{U}=0$
, as has been discussed briefly in the previous section, in the presence of cross-flow, changing
$M$
does influence the linear stability. This is displayed in figure 7(a), where we see that the effects of
$M$
are only discernible when
$M$
is small. We will discuss this issue further in the energy analysis section (§ 5.1.2). Considering non-normal linear stability, transient energy growth
$G$
is still small, even though slightly higher than in the no-cross-flow case.
Effect of
$M$
: (a) neutral stability curve with cross-flow for
$C=100$
,
$\mathit{Fe}=10^{5}$
,
${\it\alpha}=2.57$
,
${\it\beta}=0$
and
$N=250$
; and (b) transient energy growth for different
$M$
and
$T=145$
. The direction of the arrow indicates increasing
$M$
.
4.2.2. Electrohydrodynamics: high
$\mathit{Re}$
In this section, we consider the flow governed by the inertial scale, i.e. in the high-
$\mathit{Re}$
regime. To discuss the results more properly, the Reynolds number
$\mathit{Re}$
will be the free parameter, and the governing momentum equation is given by (2.10); see § 2.3. The modal stability is examined in figure 8. In panel (a), changes in the spectrum due to the additional electric field are visible. It appears that the core modes, wall modes and centre modes do not change appreciably, except that the growth rate of the most unstable mode increases. In panel (b), we plot the neutral stability curve for varying
$T$
. The pure hydrodynamic linear stability limit
$\mathit{Re}=5772.2$
is recovered by considering a minute value for
$T$
such as
$T=10^{-8}$
(we could have taken
$T=0$
, but to be compatible with (2.8) and the discussion based on that equation in other literature, we assign to
$T$
a negligibly small value). With increasing
$T$
, the system becomes more unstable, as the critical linear stability criterion becomes smaller. The reason for this is obviously due to the effect of the electric field transferring energy into the velocity fluctuations, while at the same time modifying the canonical channel flow; see table 3 in the energy analysis section (§ 5.1.3) at
$C=100$
,
$\mathit{Fe}=10^{5}$
,
$\mathit{Re}=5500$
,
${\it\alpha}=1$
and
${\it\beta}=0$
.
(a) Comparison of two spectra for
$T=10^{-8}$
and
$T=200$
at
$C=100$
,
$\mathit{Re}=5600$
,
$\mathit{Fe}=10^{5}$
,
${\it\alpha}=1$
and
${\it\beta}=0$
. (b) Effect of
$T$
on the neutral stability curve with cross-flow at
$C=100$
,
$\mathit{Fe}=10^{5}$
,
$M=\sqrt{T/\mathit{Re}}$
and
${\it\beta}=0$
. The direction of the arrow indicates increasing
$T$
. (c) Effect of
$\mathit{Fe}$
on the neutral stability curve with cross-flow at
$C=100$
,
$T=100$
,
$M=\sqrt{T/\mathit{Re}}$
and
${\it\beta}=0$
. The direction of the arrow indicates increasing
$\mathit{Fe}$
.
We also investigate the effect of charge diffusion
$\mathit{Fe}$
on the flow stability. The results concerning the neutral stability curve are shown in figure 8(c) at
$C=100$
,
$T=100$
,
$M=\sqrt{T/\mathit{Re}}$
and
${\it\beta}=0$
. For small charge diffusion (large
$\mathit{Fe}$
), the critical Reynolds number
$\mathit{Re}$
is only slightly affected by changes in
$\mathit{Fe}$
. Only when
$\mathit{Fe}=10^{3}$
does the critical Reynolds number
$\mathit{Re}$
drop noticeably, though the destabilization effect is still small. It can thus be concluded that charge diffusion has only a small influence on the dynamics of EHD cross-flow at high
$\mathit{Re}$
. This is due to the inertial scale we are considering. As we have seen in the hydrostatic EHD flow, the effect of charge diffusion is significant, considering the electric scale, i.e. at relatively small (or zero) Reynolds numbers.
It is well known that in high-
$\mathit{Re}$
Poiseuille flow the two-dimensional Orr mechanism is not the principal mechanism for perturbation energy growth over a finite time horizon. The flow is expected to become turbulent within a short time interval, even though the asymptotic growth rate of the linear system is negative. The non-normal nature of the linearized NS operator for channel flow – in physical terms, due to the base flow modulation by spanwise vorticity tilting into the wall-normal direction – suggests that transient disturbance growth during the early phase should be considered primarily.
Effect of
$T$
: (a) transient energy growth for
$C=100$
,
$\mathit{Re}=5200$
,
$\mathit{Fe}=10^{5}$
,
${\it\alpha}=0$
and
${\it\beta}=2$
; (b) contours of transient growth
$G$
in the
${\it\alpha}{-}{\it\beta}$
plane at
$C=100$
,
$\mathit{Re=5000}$
,
$\mathit{Fe}=10^{5}$
and
$T=100$
; and (c) transient growth
$G$
as a function of
${\it\beta}$
at
$C=100$
,
$\mathit{Fe}=10^{5}$
and
${\it\alpha}=0$
. The direction of the arrow indicates increasing
$T$
.
In figure 9(a), we present the transient growth
$G$
for different
$T$
. Mainly, the effect of increasing
$T$
is to enhance transient growth. The optimal initial condition that achieves maximum transient growth is shown in figure 9(b). The optimal wavenumbers for the pure hydrodynamic case, independent of
$\mathit{Re}$
, are found to be
${\it\alpha}=0$
and
${\it\beta}=2.05$
, suggesting streamwise-independent vortices as the most amplified structures (Schmid & Henningson Reference Schmid and Henningson2001). For high-
$\mathit{Re}$
EHD with cross-flow, the maximum transient growth is still found to favour streaks (
${\it\alpha}=0$
), but with a different optimal spanwise wavenumber of
${\it\beta}=2.36$
at
$C=100$
,
$\mathit{Re}=5000$
,
$\mathit{Fe}=10^{5}$
and
$T=100$
(see figure 9
b). Interestingly, for a different value of
$T$
, i.e. a different amount of potential drop across the electrodes, the optimal wavenumber would be different. For example, at
$T=50$
the optimal
${\it\beta}=2.18$
, while at
$T=10$
the optimal
${\it\beta}=2.08$
, as shown in figure 9(c). It seems that, for smaller
$T$
, approaching the regime of pure hydrodynamics, the optimal
${\it\beta}$
converges towards
${\it\beta}=2.05$
. The independence of maximum transient growth on
$\mathit{Re}$
for various
${\it\beta}$
still holds in the limit of high-
$\mathit{Re}$
EHD flow, as indicated by the dashed lines connecting the peaks of the two curves for
$\mathit{Re}=5000$
and
$\mathit{Re}=3000$
in figure 9(c). In the nonlinear regime of EHD Poiseuille flow, the influence of the electric field on the streaks has been reported in Soldati & Banerjee (Reference Soldati and Banerjee1998). These authors reported that the spanwise spacing of the low-speed streaks is approximately
$105\pm 15$
in wall units, which is different from the average spanwise spacing of the streaks in Poiseuille flow, which is 100 in wall units (see e.g. Butler & Farrell Reference Butler and Farrell1993). Thus, to some extent, our results that the spanwise spacing of the streaks changes in the linear EHD cross-flow agree qualitatively with these findings. These authors also found that the cross-flow is weakened by the electric field. This does not stand in conflict with the current results, because enhanced transient growth due to the electric field, as found here, only indicates that, in the linear phase, transition to turbulence is more rapid when compared to canonical channel flow; no conclusions can be drawn for the flow behaviour in the nonlinear regime. To more fully understand how the electric field influences streaks and streamwise vortices in the nonlinear phase of transition, a more comprehensive study of the role played by EHD in the formation and dynamics of a self-sustaining cycle (Jiménez & Pinelli Reference Jiménez and Pinelli1999) is called for.
Effect of
$T$
. (a) The optimal initial condition for
$v$
as a function of
$T$
and
${\it\beta}$
at
$C=100$
,
$\mathit{Fe}=10^{5}$
,
$\mathit{Re}=5000$
and
${\it\alpha}=0$
. (b) The optimal initial conditions for velocity and potential at
$C=100$
,
$\mathit{Re}=5000$
,
$\mathit{Fe}=10^{5}$
,
$T=100$
,
${\it\alpha}=0$
and
${\it\beta}=2.36$
. (c) Velocity vectors of the optimal initial condition in the cross-stream plane; same parameters as in panel (b). (d) Contours of the optimal response for
${\it\varphi}$
; same parameters as in panel (b).
To further investigate the effect of
$T$
, we plot in figure 10 the optimal initial conditions that achieve
$G_{max}$
in a given finite time for parameters
$C=100$
,
$\mathit{Re}=5000$
,
$\mathit{Fe}=10^{5}$
,
${\it\alpha}=0$
and different values of
$T$
with its corresponding optimal
${\it\beta}$
. In figure 10(a), the optimal initial conditions for
$v$
for various
$T$
are presented. The symmetry of the optimal
$v$
with respect to the flow centreline
$y=0$
, when the flow is close to the pure-hydrodynamics limit, is broken due to the action of the electric field in the wall-normal direction as
$T$
increases. Since the electrode with higher potential is at
$y=1$
in our case, the optimal initial conditions for
$w$
and
$v$
are tilted towards
$y=1$
(see also figure 10(b), which additionally shows the optimal initial condition for
${\it\varphi}$
). In figure 10(c), the formation of streamwise vortices in the
$y{-}z$
plane is shown; their centres are shifted upwards by the electric field. In figure 10(d), the optimal response of
${\it\varphi}$
, taking the form of waves in the spanwise direction, is displayed. Recalling that the non-modal transient growth is due to base-flow modulations arising from the tilting of spanwise into wall-normal vorticity, we can state that the variation of the optimal spanwise wavenumber for different
$T$
is the direct result of the three-dimensional nature of the non-normal linearized operator under the influence of a constant electric field pointing in the wall-normal coordinate direction.
5. Results of energy analysis
5.1. Asymptotic energy analysis
The dynamics of the disturbance energy (of the velocity fluctuations) in the limit of an infinite time horizon is examined in this section. The governing equation for the energy evolution is obtained by multiplying the linearized equation (2.8b
) by the complex conjugate velocity
$v_{i}^{\dagger }$
, i.e.
$$\begin{eqnarray}v_{i}^{\dagger }\frac{\partial v_{i}}{\partial t}+v_{i}^{\dagger }v_{j}\frac{\partial \bar{U}_{i}}{\partial x_{j}}+v_{i}^{\dagger }\bar{U}_{j}\frac{\partial v_{i}}{\partial x_{j}}=-v_{i}^{\dagger }\frac{\partial p}{\partial v_{i}}+v_{i}^{\dagger }\frac{M^{2}}{T}\frac{\partial ^{2}v_{i}}{\partial x_{j}^{2}}+v_{i}^{\dagger }M^{2}\left(\frac{\partial \bar{{\it\phi}}}{\partial x_{i}}\frac{\partial ^{2}{\it\varphi}}{\partial x_{j}x_{j}}+\frac{\partial ^{2}\bar{{\it\phi}}}{\partial x_{j}x_{j}}\frac{\partial {\it\varphi}}{\partial x_{i}}\right),\end{eqnarray}$$
taking the complex conjugate of the obtained equation, and averaging the two equations, which leaves us with
$$\begin{eqnarray}\displaystyle \frac{\partial {\mathcal{E}}}{\partial t} & = & \displaystyle -\frac{1}{2}(v_{i}^{\dagger }v_{j}+v_{j}^{\dagger }v_{i})\frac{\partial \bar{U}_{i}}{\partial x_{j}}-\frac{M^{2}}{T}\frac{\partial v_{i}^{\dagger }}{\partial x_{j}}\frac{\partial v_{i}}{\partial x_{j}}-\frac{M^{2}}{2}\frac{\partial \bar{{\it\phi}}}{\partial x_{i}}\left(\frac{\partial {\it\varphi}^{\dagger }}{\partial x_{j}}\frac{\partial v_{i}}{\partial x_{j}}+\frac{\partial {\it\varphi}}{\partial x_{j}}\frac{\partial v_{i}^{\dagger }}{\partial x_{j}}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{M^{2}}{2}\frac{\partial ^{2}\bar{{\it\phi}}}{\partial x_{i}\partial x_{j}}\left(v_{i}^{\dagger }\frac{\partial {\it\varphi}}{\partial x_{j}}+v_{i}\frac{\partial {\it\varphi}^{\dagger }}{\partial x_{j}}\right)-\frac{M^{2}}{2}\frac{\partial ^{3}\bar{{\it\phi}}}{\partial x_{i}\partial x_{i}\partial x_{j}}({\it\varphi}v_{j}^{\dagger }+{\it\varphi}^{\dagger }v_{j})\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\partial }{\partial x_{j}}\left[-\frac{1}{2}v_{i}v_{i}^{\dagger }\bar{U}_{j}-\frac{1}{2}(v_{i}^{\dagger }p+v_{i}p^{\dagger }){\it\delta}_{ij}+\frac{M^{2}}{2T}\left(v_{i}^{\dagger }\frac{\partial v_{i}}{\partial x_{j}}+v_{i}\frac{\partial v_{i}^{\dagger }}{\partial x_{j}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.\frac{M^{2}}{2}\frac{\partial \bar{{\it\phi}}}{\partial x_{i}}\left(\frac{\partial {\it\varphi}}{\partial x_{j}}v_{i}^{\dagger }+\frac{\partial {\it\varphi}^{\dagger }}{\partial x_{j}}v_{i}\right)+\frac{M^{2}}{2}\frac{\partial ^{2}\bar{{\it\phi}}}{\partial x_{i}x_{i}}\left({\it\varphi}v_{j}^{\dagger }+{\it\varphi}^{\dagger }v_{j}\right)\right],\end{eqnarray}$$
where
${\mathcal{E}}=v_{i}^{\dagger }v_{i}/2$
is the perturbation energy density of the hydrodynamic part of the spectral space. The terms in the square brackets are the transport terms, which, in the case of periodic as well as no-slip and no-penetration boundary conditions, exert no influence on the energy balance. Therefore, after integrating the above equation over the control volume
${\it\Omega}$
, we obtain
$$\begin{eqnarray}\displaystyle \int _{{\it\Omega}}\frac{\partial {\mathcal{E}}}{\partial t}\text{d}V & = & \displaystyle -\underbrace{\int _{{\it\Omega}}\frac{1}{2}(v_{i}^{\dagger }v_{j}+v_{j}^{\dagger }v_{i})\frac{\partial \bar{U}_{i}}{\partial x_{j}}\text{d}V}_{\mathsf{Pr}}-\underbrace{\int _{{\it\Omega}}\frac{M^{2}}{T}\frac{\partial v_{i}^{\dagger }}{\partial x_{j}}\frac{\partial v_{i}}{\partial x_{j}}\text{d}V}_{\mathsf{VD}}\nonumber\\ \displaystyle & & \displaystyle -\,\underbrace{\int _{{\it\Omega}}\frac{M^{2}}{2}\frac{\partial \bar{{\it\phi}}}{\partial x_{i}}\left(\frac{\partial {\it\varphi}^{\dagger }}{\partial x_{j}}\frac{\partial v_{i}}{\partial x_{j}}+\frac{\partial {\it\varphi}}{\partial x_{j}}\frac{\partial v_{i}^{\dagger }}{\partial x_{j}}\right)\text{d}V}_{\mathsf{VE1}}\nonumber\\ \displaystyle & & \displaystyle -\,\underbrace{\int _{{\it\Omega}}\frac{M^{2}}{2}\frac{\partial ^{2}\bar{{\it\phi}}}{\partial x_{i}\partial x_{j}}\left(v_{i}^{\dagger }\frac{\partial {\it\varphi}}{\partial x_{j}}+v_{i}\frac{\partial {\it\varphi}^{\dagger }}{\partial x_{j}}\right)\text{d}V}_{\mathsf{VE2}}\nonumber\\ \displaystyle & & \displaystyle -\,\underbrace{\int _{{\it\Omega}}\frac{M^{2}}{2}\frac{\partial ^{3}\bar{{\it\phi}}}{\partial x_{i}\partial x_{i}\partial x_{j}}({\it\varphi}v_{j}^{\dagger }+{\it\varphi}^{\dagger }v_{j})\text{d}V}_{\mathsf{VE3}}.\end{eqnarray}$$
Since the boundary conditions are periodic in the wall-parallel coordinate directions, it is legitimate to consider the control ‘volume’
${\it\Omega}$
in only the
$y$
direction, that is,
${\it\Omega}=[-1,1]$
and
$\text{d}V=\text{d}y$
. The first term on the right-hand side of (5.3) represents energy production from the mean shear (
$\mathsf{Pr}$
), which is zero in the hydrostatic case; the second term describes viscous dissipation (
$\mathsf{VD}$
); the third to fifth terms are the energy transfer terms between the velocity fluctuation field and the electric field (
$\mathsf{VE1}$
,
$\mathsf{VE2}$
and
$\mathsf{VE3}$
, respectively). The Einstein summation convention does not apply for the subscripts of
$\mathsf{VE1}_{ij}$
; the term
$\mathsf{VE1}_{21}$
, for example, represents
$\int _{{\it\Omega}}(M^{2}/2)(\partial \bar{{\it\phi}}/\partial y)[(\partial {\it\varphi}^{\dagger }/\partial x)(\partial v/\partial x)+(\partial {\it\varphi}/\partial x)(\partial v^{\dagger }/\partial x)]\text{d}V$
.
As has been discussed and verified for polymeric flows in Zhang et al. (Reference Zhang, Lashgari, Zaki and Brandt2013), the time variation of the normalized perturbation energy density
$R_{e}$
should be equal to twice the asymptotic growth rate of linear disturbances, i.e.
$$\begin{eqnarray}R_{e}=\frac{\displaystyle \int _{{\it\Omega}}\displaystyle \frac{\partial {\mathcal{E}}}{\partial t}\text{d}V}{\displaystyle \int _{{\it\Omega}}{\mathcal{E}}\text{d}V}=2{\it\omega}_{i},\end{eqnarray}$$
where
${\it\omega}_{i}$
denotes the growth rate of the least stable mode. We will validate this relation in the following sections and use it as an a posteriori check for our results.
5.1.1. Electrohydrodynamics without cross-flow
We apply the energy analysis for hydrostatic EHD flow with different values of charge diffusion coefficients
$\mathit{Fe}$
to probe how the electric field interacts with the velocity fluctuations. Quantitative results are listed in table 1, with the notation
$\mathsf{VD}=\mathsf{VD}_{11}+\mathsf{VD}_{12}+\mathsf{VD}_{21}+\mathsf{VD}_{22}$
and, likewise,
$\mathsf{VE}=\mathsf{VE1}_{21}+\mathsf{VE1}_{22}+\mathsf{VE2}+\mathsf{VE3}$
. There is no spanwise dependence, as
${\it\beta}=0$
. Immediately, one can make several direct observations. First, viscous dissipation is always negative for the hydrodynamics. Secondly, since the EHD flow is hydrostatic, there is no production from the mean shear,
$\mathsf{Pr}=0$
. The only terms that can lead to growth in the hydrodynamic disturbance energy density
${\mathcal{E}}$
are linked to the energy transfer from the electric field,
$\mathsf{VE}$
. The most efficient mechanism seems to be related to the term
$\mathsf{VE1}_{21}$
, which represents the interaction between the streamwise perturbed electric field and the wall-normal velocity shear under the constant effect of the wall-normal base electric field. The term
$\mathsf{VE3}$
is even negative, indicating that the electric field can absorb energy from the perturbed hydrodynamic field by an out-of-phase configuration between
${\it\varphi}$
and
$v$
(in the energy-budget equation for the perturbed electric field, one would find exactly the same term with opposite sign). Regarding the effect of charge diffusion, with increasing
$1/\mathit{Fe}$
(increasing charge diffusion) from the right column to the left in table 1, the total energy transfer
$\mathsf{VE}$
diminishes, but, at the same time, the hydrodynamic diffusion is also dissipating less energy into heat. Furthermore, even though
$\mathsf{VE}$
(and thus
$\mathsf{VE1}_{22}$
,
$\mathsf{VE2}$
and
$\mathsf{VE3}$
) decreases with rising charge diffusion, the primary mechanism of energy transfer
$\mathsf{VE1}_{21}$
transfers more energy from the electric field to the hydrodynamic fluctuations, together with less dissipation, leading to an unstable flow for the chosen parameters. Therefore, it seems that the effect of charge diffusion is to catalytically enhance the efficiency of the most productive energy transfer mechanism between the perturbed electric field and the hydrodynamics, expressed by the term
$\mathsf{VE1}_{21}$
, resulting in a lower dissipation. As a consequence, increasing charge diffusion leads to a more unstable flow.
Energy budget for modal instability of hydrostatic EHD flow for different values of the charge diffusion. The results are normalized as in (5.4). The parameters are
$C=50$
,
$M=100$
,
$T=160$
,
${\it\alpha}=2.57$
and
${\it\beta}=0$
for hydrostatic flow.
It is instructive to assess the effect of
$M$
on the linear stability (see § 4.1) with the help of the energy-budget equation (5.3). For a vanishing time derivative of the energy density, the factor
$M^{2}$
on the right-hand side can be eliminated for the hydrostatic case
$\bar{U}=0$
. In the case of cross-flow, however, by the same reasoning
$M$
will have an influence on the linear stability criterion.
5.1.2. Electrohydrodynamics with low-
$\mathit{Re}$
cross-flow
Table 2 shows the results for
$C=50$
,
$M=100$
,
$T=160$
,
$\mathit{Re}=T/M^{2}=0.016$
,
${\it\alpha}=2.57$
and
${\it\beta}=0$
with cross-flow for different
$\mathit{Fe}$
. Compared to the case without cross-flow, the results are quite similar. However, it is interesting to note that the fluctuation energy production from the mean shear
$\mathsf{Pr}$
, even though rather small, is negative, indicating that the perturbed flow field transfers energy to the base flow. Recalling the results in figure 7 of § 4.2.1, a change of
$M$
does not have a strong effect on the rate of change of the disturbance energy density
${\mathcal{E}}$
since
$\mathsf{Pr}$
is very small.
Energy budget for modal instability of EHD flow with low-
$\mathit{Re}$
cross-flow for different values of the charge diffusion
$\mathit{Fe}$
. The results are normalized as in (5.4). The parameters are
$C=50$
,
$M=100$
,
$T=160$
,
$\mathit{Re}=T/M^{2}=0.016$
,
${\it\alpha}=2.57$
and
${\it\beta}=0$
with cross-flow.
Energy budget for modal instability of EHD flow with a high-
$\mathit{Re}$
cross-flow for different values of the stability parameter
$T$
. The results are normalized as in (5.4). The parameters are
$C=100$
,
$\mathit{Fe}=10^{5}$
,
$\mathit{Re}=5500$
,
${\it\alpha}=1$
and
${\it\beta}=0$
with cross-flow.
In the case of EHD flow with a weak cross-flow (
$\mathit{Re}=0.016$
), the main mechanism for transferring energy into the hydrodynamic subsystem is still based on the potential difference across the two electrodes – the same as for the no cross-flow case. This can be confirmed by inspecting table 2:
$\mathsf{VE1}_{21}$
is the dominant energy transfer term.
5.1.3. Electrohydrodynamics with high-
$\mathit{Re}$
cross-flow
The energy analysis for the EHD Poiseuille flow at
$C=100$
,
$\mathit{Fe}=10^{5}$
,
$\mathit{Re}=5500$
,
${\it\alpha}=1$
and
${\it\beta}=0$
is summarized in table 3. Note that the production
$\mathsf{Pr}$
is diminishing with increasing
$T$
. On the other hand,
$\mathsf{VE}$
increases with larger values of
$T$
, compensating and exceeding the decrease of
$\mathsf{Pr}$
at higher
$T$
. This is consistent with the results in figure 8: higher values of
$T$
yield a more unstable flow. However, the principal mechanism underlying the flow instability is still linked to the production
$\mathsf{Pr}$
. The electric field only assumes a secondary role in destabilizing the flow, at least for the parameters considered in this case. Unlike the hydrostatic case where
$\mathsf{VE1}_{21}$
is responsible for the dominant energy transfer, in the presence of cross-flow
$\mathsf{VE1}_{22}$
becomes the most efficient agent transferring fluctuation energy
${\mathcal{E}}$
between the electric field and the perturbed velocity field.
5.2. Transient energy analysis
To investigate the cause for the increase of non-modal growth with
$T$
(recall figure 9), we formulate and perform an energy analysis for the initial-value problem of (2.16). We consider the energy density evolution over a finite time horizon following Butler & Farrell (Reference Butler and Farrell1992),
$$\begin{eqnarray}\displaystyle \frac{1}{|{\it\Omega}|}\int _{{\it\Omega}}\frac{\partial {\mathcal{E}}}{\partial t}\text{d}V & = & \displaystyle \frac{1}{|{\it\Omega}|}\int _{-1}^{1}\int _{0}^{a}\int _{0}^{b}\frac{\partial }{\partial t}\left(\frac{u^{2}+v^{2}+w^{2}}{2}\right)\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{|{\it\Omega}|}\int _{-1}^{1}\int _{0}^{a}\int _{0}^{b}\left[-\!\underbrace{uv\frac{\text{d}U}{\text{d}y}}_{\mathsf{Pr}}-\underbrace{\frac{1}{\mathit{Re}}\frac{\partial u_{i}}{\partial x_{j}}\frac{\partial u_{i}}{\partial x_{j}}}_{\mathsf{VD}}\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\,\underbrace{\frac{T}{\mathit{Re}}\frac{\partial \bar{{\it\phi}}}{\partial y}\frac{\partial {\it\varphi}}{\partial x_{j}}\frac{\partial v}{\partial x_{j}}}_{\mathsf{VE1}}-\underbrace{\frac{T}{\mathit{Re}}\frac{\partial ^{2}\bar{{\it\phi}}}{\partial y^{2}}v\frac{\partial {\it\varphi}}{\partial y}}_{\mathsf{VE2}}-\underbrace{\frac{T}{\mathit{Re}}\frac{\partial ^{3}\bar{{\it\phi}}}{\partial y^{3}}{\it\varphi}v}_{\mathsf{VE3}}\right]\text{d}z\text{d}x\text{d}y,\end{eqnarray}$$
where
$|{\it\Omega}|=2ab$
,
$a=2{\rm\pi}/{\it\alpha}$
and
$b=2{\rm\pi}/{\it\beta}$
. In the above equation, we label, as before, the first term on the right-hand side as
$\mathsf{Pr}$
(production from the mean shear), the second term as
$\mathsf{VD}$
(viscous dissipation), the third to fifth terms collectively as
$\mathsf{VE}$
(energy density transfer between the perturbed velocity field and the perturbed electric field) and the sum of all five terms as
$\mathsf{Total}$
. In this temporal evolution problem, the initial condition is the optimal one, following the procedure in § 4.2.2.
The results of our energy analysis are presented in figure 11. In panel (a) and its inset, the pure hydrodynamic result is shown, where the production
$\mathsf{Pr}$
counteracts the viscous dissipation
$\mathsf{VD}$
. In panel (b) for EHD cross-flow, we observe that the term
$\mathsf{VE}$
is insignificant, even though at
$T=100$
; this is in contrast to both the linear modal stability criterion (see figure 8) and the overall non-modal transient growth (see figures 9 and 10), where its effect is not negligible. Furthermore, production
$\mathsf{Pr}$
increases by a factor of 2–3 compared to the pure hydrodynamic flow. These results seem to indicate that, concerning the non-modal analysis, the effect of the additional electric field on the canonical channel flow is incidental, i.e. the perturbation velocity energy is only indirectly influenced by the electric field; in fact, the electric field does not induce a substantial energy transfer directly into velocity fluctuations at all and its effect is to enhance the lift-up mechanism, and therefore the production.
Energy analysis over a finite time horizon at
$C=100$
,
$\mathit{Re}=5000$
,
$\mathit{Fe}=10^{5}$
and
${\it\alpha}=0$
for: (a)
$T=10^{-8}$
and
${\it\beta}=2.05$
; and (b)
$T=100$
and
${\it\beta}=2.36$
.
Examining more closely the inset in figure 11(b), we see that the term
$\mathsf{VE}$
surpasses
$\mathsf{Pr}$
only at the very beginning of the time horizon. This is due to the high-
$\mathit{Re}$
regime we are investigating. As discussed earlier,
$\mathit{Re}=T/M^{2}$
represents the ratio of the momentum relaxation time
$L^{\ast 2}/{\it\nu}^{\ast }$
to the charge relaxation time
$L^{\ast 2}/(K^{\ast }{\rm\Delta}{\it\phi}_{0}^{\ast })$
. With the maximum
$\mathsf{Total}$
energy achieved at
$t_{max}\approx 144$
in figure 11(b), we can estimate the time horizon in the inset by observing that
$144/5000=0.029$
, a value close to the time scale depicted in the inset of panel (b). Moreover, the minimal energy growth due to the electric force validates our previous observation that the purely EHD-induced non-normality is rather small (see appendix B for a direct proof of this statement via an input–output formulation). In the case of other complex flows at high Reynolds numbers, a similar conclusion can be drawn; for instance, in viscoelastic flows (Zhang et al.
Reference Zhang, Lashgari, Zaki and Brandt2013; Brandt Reference Brandt2014), polymer stretching cannot induce disturbance growth when the fluid inertia is prevalent.
6. Discussion and conclusions
In this article, we have presented a comprehensive linear stability analysis of charge-injection-induced EHD flows between two plate electrodes, covering the hydrostatic as well as the cross-flow case, employing modal as well as non-modal tools. We intend to examine whether the linear framework is sufficient for describing the transition to turbulence of EHD flow in the early phase of perturbation evolution. It is hoped that the results presented above and summarized below would help to understand better the EHD flow instability and its transition mechanism and shed light on its flow physics as well as flow control design.
6.1. Electrohydrodynamics without cross-flow
In the hydrostatic case, the often-omitted charge diffusion is taken into account and found to have a non-negligible effect, particularly on the critical linear stability parameter
$T_{c}$
with SCL injection – a finding running contrary to a common assumption in previous studies. In those studies, a linear stability analysis predicts a critical value of
$T_{c}\approx 161$
in the strong-injection limit. This result is reproduced in our computations for a negligible value of
$1/\mathit{Fe}$
, but even for a moderate amount of charge diffusion the flow quickly becomes more unstable. Hence, we suggest that charge diffusion be accounted for in linear stability analyses and numerical simulations whenever the real physics indicates charge diffusion that cannot be neglected, as it improves not only the model of the flow physics but also the robustness of the numerics as well. In fact, the common use of total variation diminishing (TVD) schemes (Harten Reference Harten1983) in direct numerical simulations of EHD flow, which introduces artificial numerical diffusion, seems unnecessary when true physical charge diffusion could be included.
The longstanding discrepancy of the critical stability parameter
$T_{c}$
between the experimental and theoretical values, however, could not be resolved by our analysis: even though
$T_{c}$
in the SCL limit drops to 140 at
$\mathit{Fe}=10^{3}$
(a physical value according to Pérez & Castellanos (Reference Pérez and Castellanos1989)), a substantial gap remains to the experimentally measured parameter of
$T_{c}\approx 110$
. Motivated by research in subcritical channel flow, we examine other mechanisms for early transition to turbulence, specifically, transient growth due to the non-normality of the linearized EHD operator. Non-modal stability theory has been successfully applied to the variety of wall-bounded shear flows in an attempt to explain aspects of the transition process. In the case of hydrostatic EHD flow, our calculations seem to indicate that transient energy growth, as defined in (2.19), is not significant, reaching gains of
${\sim}10$
at most: rather, the flow instability is dictated by the asymptotic growth rate of the least stable mode. These results seem to indicate that below the critical
$T_{c}$
the significant energy growth observed in real EHD flow is not of a linear nature, otherwise the linear framework would succeed in detecting it. It might be hypothesized that the major energy growth mechanism in subcritical hydrostatic EHD follows a nonlinear route; nevertheless, it is only after performing a full nonlinear simulation of subcritical EHD flow that can one conclude whether its energy growth mechanism is truly nonlinear or not. Besides, these results also seem to shed some light on the flow control of hydrostatic EHD. It is now well established that, in canonical channel flow, where the perturbation energy is found to grow linearly in the early phase, a linear flow control strategy is sufficient to abate the perturbation development (Kim Reference Kim2003; Kim & Bewley Reference Kim and Bewley2007). Owing to the limited early perturbation energy growth in hydrostatic EHD, we thus suggest that different flow control methods be examined and applied in addressing the flow control of such flow. There might exist another possibility for the limited transient growth. As discussed by Atten (Reference Atten1974), the correct prediction of the linear stability criterion might require a closer comparison between the experimental conditions and the mathematical model. In light of this, one may suggest a re-examination of the charge creation and transport processes, as the current charge creation model does not seem to accommodate any efficient energy transfer from the electric to the flow field, during the linear phase.
6.2. Electrohydrodynamics with cross-flow
The flow instability and the transition to turbulence in canonical or complex channel (for example, EHD, MHD or polymeric) flows are currently not well understood. The study of complex channel flow, serving as a supplement to the investigation of canonical flows, focuses on the flow modification under the influence of external fields, for example, electric field, magnetic field or polymer stress field. The study of such flows will not only improve our understanding of these particular flow configurations, but also, more importantly, help us to better understand, during the linear, transition and turbulent phases, the dynamics of the important flow structures, for instance, streak formation and attenuation, by probing the interaction between the fluids (or flow structures) and the external fields. For example, from the point of view of flow control, research on polymer turbulence drag reduction reveals the mechanism of how the auto-generation cycles of turbulence are modified in the presence of polymer molecules (Dubief et al.
Reference Dubief, White, Terrapon, Shaqfeh, Moin and Lele2004). This has led to an even broader picture of the dynamics of turbulence. Similarly, in the case of EHD, our goal is to understand how the flow changes in response to the electric effects and to provide a physical interpretation. Below we present the results of EHD cross-flow. We differentiated low-
$\mathit{Re}$
and high-
$\mathit{Re}$
cases. For low-
$\mathit{Re}$
flow, the effects of
$M$
and
$\mathit{Fe}$
are similar to those of the hydrostatic flow, with the linear stability criterion being smaller for low-
$\mathit{Re}$
cross-flow when compared to hydrostatic flow.
The high-
$\mathit{Re}$
case is more interesting. In both modal and non-modal stability analyses, the canonical channel flow becomes more unstable, once an electric field is applied between the two electrodes. From an input–output and an energy analysis, we found, however, that energy growth directly related to the electric field is not significant and that the effect of the electric field on the flow instability is indirect. In general, in high-
$\mathit{Re}$
channel flow, the maximum transient growth is achieved by vortices aligning along the streamwise coordinate direction and generating streamwise streaks via an efficient energy growth mechanism known as lift-up. These optimal streamwise vortices are symmetric with respect to the channel centreline for standard Poiseuille flow. In contrast to other complex flows, in EHD flows the electric field, which always points in the wall-normal direction, actively participates in the formation of the streamwise rolls by accelerating the downward-moving fluid (note that, in our setting, the injector is at
$y=1$
). Consequently, this yields stronger transient growth via the lift-up mechanism, when compared to the common channel flow. In other words, the electric field provides wall-normal momentum. As has been discussed in Landahl (Reference Landahl1980) and recently reviewed by Brandt (Reference Brandt2014), the presence of wall-normal momentum will cause any three-dimensional, asymptotically stable or unstable shear flow to exhibit energy growth during a transient phase. In the present study, the role of the electric field is to provide the shear flow with such a source of wall-normal momentum and to strengthen the lift-up mechanism for EHD flow with high-
$\mathit{Re}$
cross-flow. Besides, we also find that the optimal wavenumbers for maximum transient growth increase under a stronger electric effect. Since the electric field will help to establish streamwise vortices, it may constitute a good actuator for drag reduction techniques, using the two-dimensional rolls together with a flow control strategy as described in Schoppa & Hussain (Reference Schoppa and Hussain1998) and Soldati & Banerjee (Reference Soldati and Banerjee1998).
Acknowledgements
F.M. was supported by the Italian Ministry for University and Research under grant PRIN 2010. The authors would like to thank Emanuele Bezzecchi for his initial input. M.Z. would like to thank Professor Luca Brandt of the Royal Institute of Technology (KTH), Sweden, and Dr P. Jordan of Université de Poitiers, France.
Appendix A. Code validation
We first perform a resolution check to examine the convergence of the results. The parameters in this case are
$C=50$
,
$\mathit{Fe}=2000$
,
$\mathit{Re}=6000$
,
$T=100$
,
$M=\sqrt{T/\mathit{Re}}=0.129$
,
${\it\alpha}=1$
and
${\it\beta}=0$
. The eigenspectra for four different grid resolutions
$N$
are shown in figure 12(a). The most unstable modes in these cases are listed in table 4. Satisfactory convergence, with increasing
$N$
, is observed.
Code validation. (a) Resolution check for EHD flow with cross-flow at
$C=50$
,
$\mathit{Fe}=2000$
,
$\mathit{Re}=6000$
,
$T=100$
,
$M=\sqrt{T/\mathit{Re}}=0.1291$
,
${\it\alpha}=1$
and
${\it\beta}=0$
. (b) Eigenvector component
$v$
for the most unstable mode (normalized to have the same peak value for the two codes). The parameters are
$C=50$
,
$M=10^{-11}$
,
$\mathit{Fe}=2000$
,
$\mathit{Re}=6000$
,
$T=M^{2}\,\mathit{Re}$
,
${\it\alpha}=1$
,
${\it\beta}=0$
and
$N=250$
for the EHD code; and
$\mathit{Re}=6000$
,
${\it\alpha}=1$
,
${\it\beta}=0$
and
$N=250$
for the hydrodynamic stability code.
Code validation. Resolution check for the most unstable eigenvalue of EHD flow with cross-flow; with the same parameters as in figure 12.
Secondly, the EHD eigenvector, using (3.1), is examined against a verified, pure hydrodynamic stability code employing the same spectral collocation method and solving the Orr–Sommerfeld–Squire system (see Schmid & Henningson Reference Schmid and Henningson2001), as shown in figure 12(b). The parameters for the EHD code are
$C=50$
,
$M=10^{-11}$
,
$\mathit{Fe}=2000$
,
$\mathit{Re}=6000$
,
$T=M^{2}\,\mathit{Re}$
,
${\it\alpha}=1$
,
${\it\beta}=0$
and
$N=250$
. The parameters for the hydrodynamic stability code are
$\mathit{Re}=6000$
,
${\it\alpha}=1$
,
${\it\beta}=0$
and
$N=250$
. We see that the iteratively solved EHD eigenvector is the same as the pure hydrodynamic one, which is solved directly by the MATLAB routine eig. For the computation of the transient amplification
$G$
in (2.23), it is legitimate to include only the first several, most unstable modes (Schmid & Henningson Reference Schmid and Henningson2001), i.e. eigenmodes corresponding to eigenvalues with imaginary part smaller than a certain
${\it\omega}_{i}^{c}$
are discarded (see table 5 for a validation of this approach). The reason for a minor increase of
$G$
, as more modes are included, is due to the newly incorporated eigenvectors, not because of an insufficiently refined grid.
Values of
$G$
versus different cutoff growth rates
${\it\omega}_{i}^{c}$
. The first two columns represent hydrostatic EHD flow at
$C=50$
,
$\mathit{Fe}=10^{5}$
,
$M=100$
,
$T=155$
,
${\it\alpha}=2.5$
,
${\it\beta}=0$
and
$N=250$
. The last two columns represent EHD flow with cross-flow at
$C=100$
,
$\mathit{Fe}=10^{5}$
,
$\mathit{Re}=5000$
,
$T=100$
,
${\it\alpha}=0$
,
${\it\beta}=2.36$
and
$N=250$
. Refer to the text for the definition of
${\it\omega}_{i}^{c}$
.
With the eigenvalue problem reliably solved as shown above, we present validation for the specific flows considered here. In the case of hydrostatic flow, our results for
$\mathit{Fe}=10^{7}$
, approximating the case of zero charge diffusion,
$T_{c}=160.67$
and
${\it\alpha}=2.57$
at
$C=50$
(see figure 3
a in § 4.1), are very close to the linear stability criterion reported in Atten & Moreau (Reference Atten and Moreau1972),
$T_{c}=160.75$
and
${\it\alpha}=2.569$
in the case of
$C\rightarrow \infty$
, where a coupled flow and electric system with neglected charge diffusion has been considered. This additionally implies that a value of
$C$
higher than 50 can well approximate the space-charge limit.
In the presence of cross-flow, since there exist no quantitative results for eigenvalues and eigenvectors of the EHD problem in the literature, we partially verify our results by examining the pure hydrodynamic limit of the EHD linearized problem, i.e. with electric effects being very small. This comparison is made with the stability code. The parameters are identical to the ones chosen above for the comparison of the eigenvectors. The Poiseuille base flow is
$\bar{U}=1-y^{2}$
in both codes. It is obvious that, with these selected parameters, the governing equation (2.13c
) for
${\it\varphi}$
is void of the coupling with
$v$
, since
$\unicode[STIX]{x1D647}_{v{\it\varphi}}$
in (2.15) is negligible. Therefore, the hydrodynamic equations for
$v$
and
${\it\eta}$
in (2.13a
) and (2.13b
) must reproduce the results of the stability code. This match is shown in figure 13. In panel (a), the spectra of the two codes are seen to collapse, even in the intersection region of the three eigenbranches, which is known to be sensitive due to the high non-normality of the linearized system (Schmid & Henningson Reference Schmid and Henningson2001). Additionally, the eigenmodes (blue squares) in panel (a) not matched by the hydrodynamic modes (red asterisks) are the supplementary eigenvalues linked to the presence of an electric field. The most unstable eigenvalue is shown in table 6. In panel (b), transient growth using an eigenvector expansion with
$n=71$
eigenmodes is shown. A quantitative comparison of the maximum transient growth
$G_{max}$
and its corresponding time
$t_{max}$
is presented in table 6. Agreement up to the fourth digit is achieved. The computations of
$t_{max}$
and
$G_{max}$
involve
$n=71$
eigenfunctions, each one solved with the iterative method. Even though each individual mode may be prone to small inaccuracies, figure 13(b) illustrates that transient growth (a multi-modal phenomenon) can be reliably and robustly computed using the eigenvector expansion outlined above.
Comparison between EHD Poiseuille flow and canonical Poiseuille flow: (a) eigenvalue spectrum; and (b) transient amplification of initial energy. The parameters are
$C=50$
,
$M=10^{-11}$
,
$\mathit{Fe}=2000$
,
$\mathit{Re}=6000$
,
$T=M^{2}\,\mathit{Re}$
,
${\it\alpha}=1$
,
${\it\beta}=0$
and
$N=250$
for the EHD code; and
$\mathit{Re}=6000$
,
${\it\alpha}=1$
,
${\it\beta}=0$
and
$N=250$
for the purely hydrodynamic stability code.
Code validation. The parameters are the same as in figure 13.
Results from an input–output analysis at
$C=100$
,
$\mathit{Fe}=10^{5}$
,
$\mathit{Re}=5000$
,
$T=100$
,
${\it\alpha}=0$
and
${\it\beta}=2.36$
.
Appendix B. Input–output formulation
An input–output formulation can reveal additional information on prevalent instability mechanisms by considering different types of forcings (input) and responses (output) (Jovanović & Bamieh Reference Jovanović and Bamieh2005). To demonstrate that the transient growth due to perturbative
${\it\varphi}$
is small, we compare the full responses to perturbations consisting of (i) all variables
$v,{\it\eta}$
and
${\it\varphi}$
, (ii) both
$v$
and
${\it\eta}$
, and (iii) only
${\it\varphi}$
. We thus define for these three cases different input filters
$\unicode[STIX]{x1D63D}$
, where
$\unicode[STIX]{x1D63D}_{in1}=\unicode[STIX]{x1D644}_{3N\times 3N}$
for the first case, while for the second and third cases we have
$$\begin{eqnarray}\unicode[STIX]{x1D63D}_{in2}=\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D644}_{N\times N} & \pmb{{\it 0}}\\ \pmb{{\it 0}} & \unicode[STIX]{x1D644}_{N\times N}\\ \pmb{{\it 0}} & \pmb{{\it 0}}\end{array}\right),\quad \unicode[STIX]{x1D63D}_{in3}=\left(\begin{array}{@{}c@{}}\pmb{{\it 0}}\\ \pmb{{\it 0}}\\ \unicode[STIX]{x1D644}_{N\times N}\end{array}\right).\end{eqnarray}$$
The output filter
$\unicode[STIX]{x1D63E}_{out}$
is
$\unicode[STIX]{x1D644}_{3N\times 3N}$
for all three cases: we examine the flow response in all velocities and the electric field. Consequently, the energy weight matrices should be redefined with
$\unicode[STIX]{x1D648}_{out}=\unicode[STIX]{x1D63E}_{out}\unicode[STIX]{x1D648}\unicode[STIX]{x1D63E}_{out}^{\text{T}}$
and
$\unicode[STIX]{x1D648}_{in}=\unicode[STIX]{x1D63D}_{in}^{\text{T}}\unicode[STIX]{x1D648}\unicode[STIX]{x1D63D}_{in}$
. After applying a Cholesky decomposition to these energy weight matrices, we obtain
$\unicode[STIX]{x1D641}_{out}$
and
$\unicode[STIX]{x1D641}_{in}$
for a formulation based on the
$L_{2}$
-norm. Finally, the maximum transient growth
$G$
over a finite time interval is given by
$$\begin{eqnarray}\displaystyle G(t) & = & \displaystyle \max _{{\bf\gamma}_{0}}\frac{\Vert {\bf\gamma}_{out}(t)\Vert _{E_{out}}}{\Vert {\bf\gamma}_{in}(0)\Vert _{E_{in}}}=\max _{{\bf\gamma}_{0}}\frac{\Vert \boldsymbol{{\mathcal{T}}}{\bf\gamma}_{in}(0)\Vert _{E_{out}}}{\Vert {\bf\gamma}_{in}(0)\Vert _{E_{in}}}=\max _{{\bf\gamma}_{0}}\frac{\Vert \unicode[STIX]{x1D641}_{out}\boldsymbol{{\mathcal{T}}}{\bf\gamma}_{in}(0)\Vert _{2}}{\Vert \unicode[STIX]{x1D641}_{in}{\bf\gamma}_{in}(0)\Vert _{2}}\nonumber\\ \displaystyle & = & \displaystyle \max _{{\bf\gamma}_{0}}\frac{\Vert \unicode[STIX]{x1D641}_{out}\boldsymbol{{\mathcal{T}}}\unicode[STIX]{x1D641}_{in}^{-1}\unicode[STIX]{x1D641}_{in}{\bf\gamma}_{in}(0)\Vert _{2}}{\Vert \unicode[STIX]{x1D641}_{in}{\bf\gamma}_{in}(0)\Vert _{2}}=\Vert \unicode[STIX]{x1D641}_{out}\boldsymbol{{\mathcal{T}}}\unicode[STIX]{x1D641}_{in}^{-1}\Vert _{2}\nonumber\\ \displaystyle & = & \displaystyle \Vert \unicode[STIX]{x1D641}_{out}\unicode[STIX]{x1D63E}_{out}\text{e}^{t\boldsymbol{{\mathcal{L}}}}\unicode[STIX]{x1D63D}_{in}\unicode[STIX]{x1D641}_{in}^{-1}\Vert _{2}.\end{eqnarray}$$
We report the transient growth results for the above three cases in table 7 at
$C=100$
,
$\mathit{Fe}=10^{5}$
,
$\mathit{Re}=5000$
,
$T=100$
,
${\it\alpha}=0$
and
${\it\beta}=2.36$
. We observe that perturbations solely in
${\it\varphi}$
(case (iii)) exhibit transient growth two orders smaller than in the other two cases. For cases (i) and (ii) the transient growth characteristics are nearly identical, which suggests that the non-normality of the linear operator is mainly related to the hydrodynamics.
