Electroconvection in a dielectric liquid between two concentric half-cylinders with rigid walls: Linear and nonlinear analysis

PHYSICAL REVIEW E 97, 023106 (2018) Electroconvection in a dielectric liquid between two concentric half-cylinders with rigid walls: Linear and nonlinear analysis Pedro A. Vázquez* Departamento de Física Aplicada III, ETSII, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain Alberto T. Pérez Departamento de Electrónica y Electromagnetismo, Facultad de Física, Universidad de Sevilla, Avda. de Reina Mercedes s/n, 41013 Sevilla, Spain Philippe Traoré Institut PPRIME, Département Fluide-Thermique-Combustion, Boulevard Pierre et Marie Curie, BP 30179, 86962 Futuroscope-Chasseneuil, France Jian Wu† School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001, People’s Republic of China (Received 2 October 2017; published 20 February 2018) We study the linear stability and nonlinear behavior of the electroconvection between two concentric halfcylinders with no-slip conditions on all boundaries. The no-slip condition makes impossible to apply the standard modal approach. Hence, we apply a finite element technique similar to the one we have used in a previous paper about the electroconvection in a rectangular enclosed domain. When compared to the classical problem of electroconvection between two full concentric cylinders, the linear criterion is higher, due to the viscous shear introduced by the lateral walls. As a consequence, the structure of the eigenmodes is very different. There is a repulsion between modes with the same symmetry, forcing pairs of modes to cross each other repeatedly. For inner injection and small value of the ratio between the inner and outer radii the no-slip condition changes the nature of the bifurcation from subcritical to supercritical, while it is always subcritical for outer injection. To understand this behavior, we perform a modal analysis using the eigenfunctions obtained from the linear stability analysis as modal basis. We show that the supercritical branch is originated by the nonorthogonality of the modes when no-slip boundary conditions are imposed. These mechanism explains the previously unexplained appearance of the supercritical branch in the enclosed rectangular configuration. DOI: 10.1103/PhysRevE.97.023106 I. INTRODUCTION We study the linear stability and the nonlinear bifurcation of the electrohydrodynamic (EHD) convection in a dielectric liquidbetweentwohalf-cylinderswithno-slipboundarycondi- tionsin thestrong injectionregime.EHD isan interdisciplinary science dealing with the interaction of fluids and electric fields. It has a great number of applications in industrial processes: EHD pumping [1–3], heat transfer enhancement [4–6], turbulent mixing [7], EHD spraying and atomization [8,9], and biomicroelectromechanical systems (bio-MEMS) and nanotechnology applications [10–12]. One of the classical problems in EHD is the convection generated between two infinite cylindrical and concentric electrodeswith radiiRi<R o.The spacebetweenthecylinders is filled with a perfectly dielectric liquid. An electric potential difference is applied between the cylinders. When the applied voltage is above a critical value, injection on one or both of the *pv[email protected] †[email protected] electrodes occurs. We will consider only unipolar injection, inner if the injecting electrode is the one with radius Riand outer if it is the one with radius Ro. The injected charge having the same polarity as the injecting electrode, the ions are pushed away from the electrode, setting the liquid into motion. The linear stability analysis of the full 3D problem was performed by Agräit and Castellanos [13]. In a recent paper, Fernandes et al. [14] studied the linear stability of this problem, both numerically and analytically, although they did not considered the nonlinear regime. In a previous paper [15], the authors studied numerically the linear stability and the nonlinear bifurcation of the full cylinder problem simplified to a 2D configuration. The bifurcation is subcritical, that is, once the critical voltage is attained, the stationarybranchisunstable.Thevelocityofthefluidincreases until, in the final steady state, the maximum velocity of the liquid is greater than the drift velocity of the ions. The charge distribution is then controlled by the velocity field and regions void of charge appear. This nonlinear bifurcation is similar to what happens in one of the most classical problems in EHD, the convection generated by two infinite parallel electrodes immersed in a perfectly dielectric liquid and subjected to a 2470-0045/2018/97(2)/023106(15) 023106-1 ©2018 American Physical Society VÁZQUEZ, PÉREZ, TRAORÉ, AND WU PHYSICAL REVIEW E 97, 023106 (2018) electric voltage applied between the electrodes. However, in two recent papers [16,17] the authors considered a variation of this two-plate problem, with the infinite domain substituted by a finite closed rectangular domain with no-slip boundary conditionsappliedonalltheboundaries.Inthefirstpaper,using numerical simulations, the authors showed that the bifurcation may become supercritical for some values of the aspect ratio of thedomain.In this case,steady convectiverollswith maximum velocity smaller than the drift velocity of the ions can be formed. In the second paper the authors studied the linear stability of the enclosed problem, and showed that the structure of the modes is quite different from that of the infinite case. Here, we study a similar modified version of the twocylinder configuration. We consider the electroconvection between two concentric half-cylinders, of radii Ri<R owith no-slip boundary conditions on all the boundaries. An electric voltage is applied between the half-cylinders, with the inner cylinder connected to a higher voltage and the outer one grounded for inner injection, and the reversed way for outer injection. We will assume that unipolar injection takes place in both cases. We study the linear stability of this problem, along with the nature of the nonlinear bifurcation. We restrict our analysis to the strong injection regime, although the main conclusions are applicable to arbitrary injection strength. The change of geometry from the planar one introduces two new elements influencing the system: the injection direction and the radius ratio. We will see that the behavior of the system is different when the injection of charge occurs at the inner or the outer electrodes. This is especially true for the change of the nature of the bifurcation from subcritical to supercritical. This different behavior offers some clues of the explanation of this change of nature. Theresults ofthe linearstabilityanalysisaresimilartothose obtained in Pérez et al. [17]. The no-slip boundary conditions increases the threshold value for injection and changes the structure of the modes. The modes with the same symmetry repel each other as the ratio between the inner and outer radii changes. For the nonlinear aspects, the presence of no-slip walls can change the nature of the bifurcation, as it was discussed by the authors in Wu et al. [16] in a rectangular domain with top and bottom flat electrodes. With free-slip lateral walls the bifurcation is always subcritical. When no-slip boundary conditions are imposed the bifurcation becomes supercriticalforsomerangesofvaluesof.Inthehalf-cylinder configurations,the ratiobetween theradii ofthe innerand outer cylinders plays the role of the aspect ratio in the rectangular configuration. When free-slip conditions are imposed on the flat walls the bifurcation is always subcritical, as it is the case withfullconcentriccylinders[15].Whenno-slipconditionsare imposed,andtheinjectionofionsoccursontheinnerelectrode, thebifurcationmaybecomesupercriticalforsmallvaluesofthe radiiratio.Asthisratioisincreaseditbecomessubcriticalagain and remains that way. However, when the injection occurs at the outer electrode, the bifurcation is subcritical for all the values of the radii ratio that we have explored. We have performed a modal analysis of the evolution of the flow pattern using the eigenmodes for each configuration computed in the linear stability analysis as a basis. In the free-slip case the modes are essentially orthogonal, while in the no-slip configuration they are not. We show that, for inner injection, the supercritical nature of the bifurcation can be explained by this nonorthogonality, as several modes can extract energy from the electric field and their final amplitudes areofsimilar values.The bifurcationbecomessubcriticalwhen the second mode is excited. This mode modifies the charge distribution and throws the system away from the supercritical branchinto theupperone.Inthe outerinjection case,themodes are still nonorthogonal, but the critical value of the stability parameter is much higher than in the inner injection regime. As a consequence, the modes grow so much that the velocity field takes control of the distribution of electric charge and a region void of charge is created, which is a characteristic of the subcritical bifurcation. This mechanism explains also the previously unexplained appearance of the supercritical branch in the closed rectangular configuration described in reference Pérez et al. [17]. The paper is organized as follows. In Sec. II we state the physical domain, the governing equations, the boundary conditions of both free- and rigid-walls cases and inner and outer injection configurations. The relevant nondimensional parameters are also set. Section III is devoted to the description of the method dealing with the linear stability analysis, as well as analyzing the results. As it was the case in Pérez et al. [17], the no-slip boundary conditions on the flat walls prevents the use of the standard modal analysis used by Agräit and Castellanos [13]. A semianalytical technique based on finite elements is applied. In Sec. IV we examine the nonlinear behavior. First, as a test of the numerical technique, we present the linear stability criteria obtained from the simulations and compare them to the values obtained from the eigenvalue analysis. After that, we discuss the nature of the bifurcation for each configuration. Finally, we perform a modal analysis to understand the nonlinear behavior for each configuration. The main conclusions are summarized in Sec. V. II. PROBLEM FORMULATION A. Basic equations and boundary conditions We consider a dielectric liquid confined in the 2D region between two-half cylinders, of radius Riand Ro, with Ri<R o (see Fig. 1). The liquid is assumed to be incompressible, Newtonian, isothermal, and perfectly insulating, with mass density ρ, kinematic viscosity ν, and absolute permittivity ε. When inner injection is considered, a stationary electric potential a is applied to the inner cylinder and the outer is grounded. The applied potentials are reversed when the outer injection case is studied. The electric field acts upon these charges and, if the applied potential is high enough, put the liquid in motion. We will make the assumption of unipolar injection, where only one of the electrodes injects charges and only one type of charge carriers is present. We will further assume that the density of injected charge on the injecting electrode, qi, is constant (autonomous injection). The collecting electrode is assumed to behave as an open boundary for the ions. These are common assumptions in EHD problems [18]. The constitutive law for the current density is J=qKE−D∇q+qu,(1) where qis the volumetric charge density, uis the fluid velocity, and Eis the electric field. The first term in Eq. (1)isthe 023106-2 ELECTROCONVECTION IN A DIELECTRIC LIQUID … PHYSICAL REVIEW E 97, 023106 (2018) FIG. 1. Computational domain with dimensional boundary conditions for inner injection. For outer injection the boundary conditions on the electric potential are interchanged and the boundary condition for the charge density is applied at the outer electrode. drift current, Kbeing the ionic mobility of the charge carriers. The second term is the contribution from diffusion where D is the charge diffusion coefficient. The last term is due to the convection of the electric charge by the velocity of the liquid u. In EHD problems, usually diffusion is negligible compared to drift [18]. Hence, the charge conservation equation becomes: ∂q ∂t +∇·[q(KE+u)]=0.(2) As boundary condition for the charge we impose q=qion the injecting electrode. Figure 1shows the dimensional boundary conditions for inner injection. The equations describing the problem are the Poisson equation, the Navier-Stokes equations, and the charge conservation equation. To build the nondimensional equations we choose the following scales: x,y ∼d=(Ro−Ri),t∼d2/Ka,u∼Ka/d, ∼a,E∼a/d, q ∼εa/d2, p∼ε2 a/d2. (3) Wehaveusedthedistancebetweenthecylinders,d=Ro−Ri, to define the spatial scale. With these choices, the nondimensional equations of the problem become ∂u ∂t +(u·∇)u=−M2∇p+M2 Tu+M2qE, ∇·u=0,  =−q, E=−∇, ∂q ∂t +∇·[q(u+E)]=0.(4) As the fluid is considered to be isothermal, the dielectric force is not relevant. The electrostriction force, being the gradient of a scalar magnitude, can be included in the gradient of the pressure. Therefore, the force term in the first equation in Eq. (4) is restricted to the Coulomb force [18]. The nondimensional boundary conditions are r=¯ Ri= 1−:=1(inner),u x,uy=0,q=C(inner), =0(outer), r=¯ Ro=1 1−:=0(inner),u x,uy=0, =1(outer),q=C(outer), θ=0,π :∂/∂y =0ux,uy=0 (rigid walls) ∂ux/∂y,uy=0 (free walls). (5) Here, ¯ Ri=Ri/d,¯ Ro=Ro/d,and=Ri/Ro.Figure2shows the geometry and boundary conditions of the problem in nondimensional form. In the rest of the paper we will use the terms lateral walls or flat sides to describe the θ=0,π boundaries. The nondimensional parameters appearing in the equations and boundary conditions are T=εa ρνK ,C=qid2 εa ,M=1 Kε ρ1/2 ,=Ri Ro .(6) The parameter Tis the ratio between the electric and viscous terms in Navier-Stokes equation. It acts as a kind of electric Rayleigh number and will stand for the stability parameter. The mobility number Mis the ratio between the so called hydrodynamic mobility and the ionic mobility [19]. It only depends on the fluid properties. The injection number C measures the injection strength. In this paper we will only consider the strong injection case (C=10). III. LINEAR STABILITY ANALYSIS We study here the linear stability of the hydrostatic state. We apply the same technique detailed in Pérez et al. [17], adapted to the half-cylinders geometry. For small values of the applied voltage viscosity prevents any motion of the liquid. As the voltage is increased, the liquid is eventually put in motion when a certain threshold value is attained. We study the linear stability of the system using the standard method of small perturbations. The steady state is perturbed with small perturbations and we study how they grow. If all perturbations decay, the system is stable. The smallest value of the stability parameter for which the perturbations grow in time is the critical value Tc. 023106-3 VÁZQUEZ, PÉREZ, TRAORÉ, AND WU PHYSICAL REVIEW E 97, 023106 (2018) FIG. 2. Computational domain with nondimensional boundary conditions for inner injection. For outer injection the boundary conditions on the electric potential are interchanged and the boundary condition for the charge density is applied at the outer electrode. A. Hydrostatic solution The hydrostatic solution corresponds to a steady state characterized by the liquid at rest. The hydrostatic profiles of the charge and the electric field can be computed from equations Eq. (4) imposing that the temporal derivatives are zero andthe velocity U0=0. Theequations for thehydrostatic state are −∇P0+Q0E0=0, 0=−Q0, E0=−∇0, ∇·(Q0E0)=0.(7) The solutions only depend on the radial coordinate, hence the electric field can be expressed as E0=E0(r)erand the charge profile as Q0=Q0(r). The equations are 1 r d(rE0) dr =Q0, d(rQ0E0) dr =0.(8) The boundary conditions are inner injection: Q0(r=¯ Ri)=C, ¯ Ro ¯ Ri E0dy =1, outer injection: Q0(r=¯ Ro)=C, ¯ Ri ¯ Ro E0dy =1.(9) The solutions have the form E0(r)=RcEs r1+C EsRcr2−R2 c, Q0(r)=C 1+C EsRcr2−R2 c .(10) Here, Rc=¯ Ri(¯ Ro) stands for inner (outer) injection. Esis the value of E0at the injector. It is defined by Eq. (9), which is an implicit equation in Es. This equation can be solved by numerical integration for each value of C.Esis positive for inner injection and negative for outer injection. B. Perturbation equations We look for perturbations of the base magnitudes. We express the physical magnitudes as u(x,y,t)=u(x,y,t), p(x,y,t)=P0(x,y)+p(x,y,t), (x,y,t)=0(x,y)+(x,y,t), E(x,y,t)=E0(x,y)+E(x,y,t), q(x,y,t)=Q0(x,y)+q(x,y,t).(11) Theprimedquantitiesaresupposedtobemuchsmallerthanthe corresponding base quantities. Introducing these expressions into Eq. (4), taking into account Eqs. (7), and keeping only the first order terms in the primed quantities leads to the following set of equations: ∂u ∂t =−M2∇p+M2 Tu+M2(Q0E+qE0), ∇·u=0, =−q, E=−∇, ∂q ∂t +∇·[Q0(u+E)+qE0]=0.(12) The standard way to proceed is to write the primed equations as a product of an exponential function depending on time and a function depending on the spatial variables alone, ⎛ ⎜ ⎝ u p  q ⎞ ⎟ ⎠=⎛ ⎜ ⎝ U(x,y) (x,y) V(x,y) Q(x,y) ⎞ ⎟ ⎠eσt.(13) The symbol σdenotes the growth rate, which can be complex. When Im(σ)=0, the principle of exchange of stabilities applies. As it is discussed in Atten and Moreau [20]forthe plane-plane geometry, it is not possible to prove rigorously that this principle applies in the free walls case. But the numerical results show that this is actually the case. Here, we will assume that Im(σ)=0 for both free and rigid walls situations. Then, for the marginal state that corresponds to the onsetofthe motionitis σ=0. Thecomparisonwith theresults issued from the direct numerical computation will validate this assumption. In the free walls case, when symmetric boundary conditions are applied at the flat sides, the spatial functions can be developed in terms of normal modes, adapting the 3D analysis in Agräit and Castellanos [13] to the 2D geometry. As an example, the charge perturbation would be expressed, in polar coordinates,asQ(r,θ)=f(r)eimθ,thenumbermbeingrelated to the number of rolls. The no-slip boundary conditions on the flat sides prevent us from applying the same procedure in the rigid walls case, and a different strategy based on a numerical computation must be envisaged. 023106-4 ELECTROCONVECTION IN A DIELECTRIC LIQUID … PHYSICAL REVIEW E 97, 023106 (2018) We introduce Eqs. (13) into Eqs. (12). We are interested in the margin stability, σ=0, so the time derivative does not appear. In this way, the number Mdisappears too and will have no influence on the linear stability. The pressure can be rescaled to include the parameter Tin its definition. We also drop for convenience the primes from the magnitudes. Then the stability equations can be rewritten as ∇p−u=T(−Q0∇V+qE0), ∇·u=0, V =−q, ∇·[Q0(u−∇V)+qE0]=0.(14) The boundary conditions are now r=¯ Ri= 1−:V=0,u x,uy=0,inner injection : q=0, r=¯ Ro=1 1−:V=0,u x,uy=0,outer injection : q=0, θ=0,π :∂V/∂y =0,u x,uy=0 (rigid walls), ∂ux/∂y,uy=0 (free walls). (15) Equations (14) along with boundary conditions Eqs. (15) define an eigenvalue problem. The secular determinant leads to a relation between the relevant parameters of the form F(T,C,)=0.(16) We take Tas the stability parameter, that is, for Cand fixed, the different values of the eigenvalue Tgive the excitation thresholds of the different modes. The smallest of these eigenvalues, Tc, will be the threshold for the linear stability. The generalized eigenvalue problem is solved using finite elements. To build the matrices, we use DOLFIN [21], an interface to the FEniCS package [22]. This is a framework for automated solution of partial differential equations based on the Finite Element method. The package allows the construction of the relevant matrices using a high level formulation of the weak problem. It also handles the computation of eigenvalues and eigenfunctions through the SLEPc package [23], a software library for the solution of large scale sparse eigenvalue problems. We refer to Pérez et al. [17]formore details. C. Results We have considered the case C=10, both with free and rigid walls, for both inner and outer injection. We have used triangular meshes built over a 200 ×100 structured grid, that is,with40000triangularelements.Themeshesarethinnernear the injecting half-cylinder, and coarser next to the collecting one. 1. Free walls case We use the free walls case as a validation of our method. In this configuration the standard modal expansion can be used to study the linear stability. Let us point out that, when free-slip boundary conditions are imposed on the flat sides in Fig. 2, the linear stability problem becomes equivalent to the complete cylinders configuration. Agräit and Castellanos [13] studied the linear stability of the 3D configuration with whole cylinders by applying the classical modal expansion. We have adapted their study to the 2D case. We have solved the resulting eigenvalue problem using the bvp4c function in MATLAB. This function implements a collocation method with C1 piecewise cubic polynomials [24]. Figure 3shows the critical value Tcas a function of for inner and outer injection obtained with the standard modal analysis and the one developed here for purpose based on finite elements. The results obtained with the two different methods are in excellent agreement. In the strong injection regime under consideration, inner injection is always less stable than outer injection for all . In addition, as approaches 1, the critical value for both injection configurations tends towards the critical value for infinitely large plate-plate case, Tp c=164.1 [20]. The results of the FE method for very close to 1 are not provided. The reason is that, as →1, the computational FIG. 3. Linear threshold computed with Finite Elements (FE) and modal analysis for inner and outer injection for C=10. The results computed with FE for rigid walls are also plotted. As approaches 1, the value of Tcgoes to the critical value corresponding to the plane case, Tp c, for all the cases. The number of rolls corresponding at each hump is indicated for the inner injection branch in the free walls case. 023106-5 VÁZQUEZ, PÉREZ, TRAORÉ, AND WU PHYSICAL REVIEW E 97, 023106 (2018) ( a()b ) FIG. 4. Eigenvalues of the first four modes vs. for the free walls case for inner injection (a) and outer injection (b). The numbering in the legend corresponds to the number of rolls observed in the half-cylinders configuration. Only a portion of each mode is showed. domain becomes very large, considering ¯ Ri=/(1 −) and ¯ Ro=1/(1 −). Thus, a larger number of elements is required and the computation time increases dramatically. Anyway, for values of close to 1 the results are similar to those obtained for a planar geometry. InFig.3wecanobservesomehumps inbothbranches. Each humpcorresponds toa differentnumber ofvelocityrolls forthe velocityfield(relatedtothemvalueinthemodalanalysis[13]). The number of rolls observed in the half-cylinder geometry are indicated up to 4 for the inner injection case. Figure 4shows the eigenvalues of the first four modes for each value of .The cusps in the lower branch of Fig. 3corresponds to the crossing of two different modes as displayed in Fig. 4. Note that we have plotted only a portion of each mode. If we prolong the branches all the modes cross each other, although at values of Thigher than the critical value for the corresponding value of . We have plotted in Fig. 5the eigenfunctions of the most unstable modes in the inner injection case for =0.05 and =0.15. The former corresponds to a point inside the first hump in Fig. 3, while the latter lies inside the second hump. The stream function for the first mode is antisymmetric with respect to the x=0 line, while the second mode is symmetric. This pattern repeats itself with higher modes. Also, all the velocity rolls are equivalent, in the sense that they correspond to the modal factor eimθ in the modal expansion. They only adapt their sizes to the total number of rolls present and the space available between the electrodes. The crossing of the modes for the free walls case is similar to the behavior in the planar injection problem discussed in Pérez et al. [17]. However, there is a fundamental difference in the interpretation of the results. In the planar case, the free walls configuration is a way of modeling two infinite parallel planes. In this problem, all the modes have the same minimum value of the linear threshold Tp c=164.1. In a real experiment with a sufficiently large domain, the system would choose this absolute minimum value as the linear stability threshold. In the cylinder configuration the situation is different. Indeed, if we couldperforman experimentkeeping the2D configuration, the value of the linear criteria Tcwould be fixed by the geometric parameter , as plotted in Fig. 3. (a ( )b) FIG. 5. Stream function contour plots of the eigenfunctions of the most unstable modes for =0.05 (a) and =0.15 (b) for the free-walls case. The number of rolls is related with the number mof the corresponding mode in the modal expansion. 023106-6 ELECTROCONVECTION IN A DIELECTRIC LIQUID … PHYSICAL REVIEW E 97, 023106 (2018) (a ()b) FIG. 6. First four eigenvalues as a function of for rigid walls. Plot (a) corresponds to inner injection, and plot (b) to outer injection. The number of rolls of the dominant mode is showed up to four. 2. Rigid walls case Figure 3also shows the critical value Tcas a function of  for rigid walls for both inner and outer injection, along with the values for the free wall configurations. In both injection configurations the linear criterion of the free walls case is the lowest one, especially for small values of . The difference becomes negligible as increases above 0.4. This behavior is due to the viscous shear introduced by the lateral walls in the rigid case. This damping affects only the convective cells in the vicinity of the flat sides. Hence, its influence on the global flow decreases as the number of cells increases, that is, as increases. In both curves, humps similar to those present in the free walls configuration are visible, although the cusps are displaced. Each hump corresponds also to a different number of rolls. As approaches 1, the value of Tctends to Tp c=164.1, the critical value for two infinite parallel plates. The same behavior has already been found for the free walls case. We emphasize that the numbering of the modes does not correspond to the number of rolls appearing in the flow. As increases, the first mode develops two lateral smaller rolls. They grow in size and eventually, when the first mode is again the most unstable one, a pattern of three fully developed convective rolls is observed. Then two smaller lateral rolls start to develop and the process repeats itself. A similar behavior is exhibited by the second mode (see Fig. 7). Figure 6plots the value of Tas a function of for the first four most unstable modes for the rigid walls case, for inner and outer injection. Note that the structure of the Tversus  diagram is completely different from that of Fig. 4. The first and second modes cross each other several times, but they don’t cross the other modes. It is the same for the third and fourth modes, and other higher ones. This behavior is similar to what was observed in the 2D rectangular cavity configuration described in Pérez et al. [17], and it can be explained by a similar argument. Figure 7shows the stream functions of the most unstable mode for =0.05 and =0.20 with inner injection. The former corresponds to the first mode and the latter is related to the second one. Both are the most unstable modes for the corresponding values of . The modes with the same parity have the same symmetry with respect to the x=0 axis: antisymmetric for the odd modes and symmetric for the even ones. The paths of the modes with the same symmetry (a ( )b) FIG. 7. Stream function contour plots of the eigenfunctions of the most unstable modes for =0.05 (a) and =0.20 (b) for the rigid walls case and inner injection. The first one is antisymmetric and the second one is symmetric. The counter rolls in the corners appear because of the no-slip boundary conditions for velocity. 023106-7 VÁZQUEZ, PÉREZ, TRAORÉ, AND WU PHYSICAL REVIEW E 97, 023106 (2018) ( a()b ) FIG. 8. Linear stability criteria as a function of for C=10 for inner (a) and outer (b) injection in the case of rigid walls. The discrete points correspond to values computed from the numerical simulations, while the solid lines have been computed with the FEM-based stability analysis. do not cross each other, so mode 3 repels mode 1, and mode 4 repels mode 2, forcing modes 1 and 2 to cross each other. The mathematical reason is related to the conditions needed to get degenerate eigenvalues in the eigenvalue problem. In general, in a two-parameter system, two paths of symmetrybreakingbifurcationpointsintersectonlyiftheybreakdifferent symmetries (or if the solutions belong to different symmetry classes), whereas the path of bifurcation points that break the same symmetry cannot cross one another. For a more complete explanation we refer the reader to our previous paper [17], and also to other works on the Rayleigh-Bénard-Marangoni problem, where this behavior has also been reported (see for example refs. [25–27]). The free walls case is a degenerate problem, and the modes can cross each other without any restriction. It is worth to not that, while for the free walls configuration all the convective rolls are equivalent, this is not the case for the rigid walls case. The new convective roll in the central region in the plot (b) of Fig. 7is not a replica of the rolls in plot (a) of the figure. Indeed, this is due to the presence of the counter-rotating rolls close to the corners, induced by the no-slip boundary conditions for velocity. IV. NONLINEAR ANALYSIS In this section, we consider the nonlinear aspects of the problem, in particular the nature of the bifurcation and how it changes with the values of the parameter. The complexity of the mathematical problem makes very difficult to perform an analytic or semianalytic study, similar to what we have done for the linear stability. Therefore, we have solved numerically the full problem defined by the set of Eqs. (4) along with the associated boundary conditions Eqs. (5). A. Numerical methods We solve the set of Eqs. (4) using our in-house code Oracle. This code implements a second order in time and space finite volume method [28]. The computational domain is discretized with a structured grid consisting of nonorthogonal quadrilaterals. All the meshes have 150 cells along the radial direction, but the number of cells along the azimuthal direction increases with the value of . For example, the mesh is 330 ×150 for =0.05 and 628 ×150 for =0.5. In all the cases, the mesh is thinner near the injecting electrode. We use a colocated arrangement, i.e., all variables are stored at the center of each control volume. The central differencing (CD) scheme and the improved deferred correction (IDC) scheme [29] are used to compute the convective and diffusive fluxes in the Navier- Stokes equations, respectively. For time integration, we use a second-order semi-implicit three time levels (I3L) scheme [28]. The velocity-pressure coupling algorithm is undertaken by the SIMPLE algorithm [30], as the fluid is assumed to be incompressible. The charge transport equation is hyperbolic, therefore it requires the use of special methods [31,32]. To prevent the development of spurious numerical oscillations in the charge density distribution, it is recommended to use nonoscillating, nondiffusive, and bounded schemes [33,34]. In this study, we have chosen a TVD scheme, the third-order smooth monotonic algorithm for real transport (SMART) scheme of Gaskell and Lau [35]. We refer to Refs. [31,32] for additional details. The improved least-squares approach [36] is used to compute all gradients including the electric field. The simulations related to the projection on the modal basis have been made with finite elements using the COMSOL multiphisics software. The numerical computations of the charge transport equation have been made with the transport of diluted species module. This module includes a crosswind stabilization scheme by Do Carmo and Galeao [37]. Although more diffusive than the SMART scheme described above, it controls well the undershoots and overshoots, so that the values of electric charge density remain positive at all times. 023106-8 ELECTROCONVECTION IN A DIELECTRIC LIQUID … PHYSICAL REVIEW E 97, 023106 (2018) (a ()b) FIG. 9. Bifurcation diagrams for inner injection, rigid walls, =0.5(a) and =0.05(b). The bifurcation is subcritical in the first case and supercritical in the second one. B. Results All the simulations in this paper have been performed with C=10 and M=60. 1. Linear stability As a validation, we reobtain from the direct numerical simulations the linear stability criteria discussed previously, for inner and outer injection in the case of no-slip walls. We start from the hydrostatic solution of Eqs. (4) and boundary conditions Eqs. (5) for a given value of . Then we start a simulation computing all the physical magnitudes for different values of Tuntil the onset of fluid motion. Then we obtain the critical value Tcfor each from the growth factor of the perturbations. The procedure is more detailed in Wu et al. [15]. Figure 8plots the computed thresholds for several values of for inner [Fig. 8(a)] and outer [Fig. 8(b)] injection. The solid lines correspond to the values obtained with the Finite Element (FEM)analysisdetailed inSec. III. Theagreement is verygood for both cases. The cusps observed in the plot correspond to the crossing of modes described in Sec. III. These cusps are clearer in the inner injection configuration. 2. Inner injection In Wu et al. [15] we discussed the non linear phenomena associated to two coaxial complete cylinders. The bifurcation in that case was subcritical. This means that, once the value of Tis greater than the linear threshold Tcand the motion starts, the system jumps to a steady state where the maximum velocity of the fluid is always greater than the maximum drift velocity induced by the electric field. This behavior can also be obtained in the half cylinder configuration. Plot (a) in Fig. 9 shows the bifurcation diagram for =0.5 and inner injection. For T>T c=124.34, the system jumps to the upper branch. If Tis further increased, the maximum velocity increases as indicated by the vertical upwards arrow. The figure depicts the typical hysteresis associated to subcritical bifurcations. Once we reach the upper branch, if we start decreasing Tthe system stays on the upper branch even for values T<T c. This goes on down to a nonlinear critical value Tf(about 89.0 for this value of ). If T<T fthe motion stops and the system goes back to the hydrostatic state. However, the nature of the bifurcation can change for some values of . Plot (b) in Fig. 9depicts the bifurcation diagram for inner injection and =0.05. In that case, we can see that above the linear critical value Tc1=92.63 a supercritical stable branch appears. Along this branch we find steady states where the maximum velocity of the fluid is lower than the maximum drift velocity. As the value of Tincreases we follow this supercritical branch, up to a new critical value Tc2= 104.5, where the system jumps to a finite amplitude branch similar to that one observed for =0.5. In this branch the maximum velocity of the liquid is greater than the maximum drift velocity. If we decrease now the value of Twe follow this branch down to a non linear threshold Tf≃69.0. Figure 10 shows the electric charge distribution for the stationary states corresponding to points B(a), C(b), and D(c) in Fig. 9.The first point is on the supercritical branch, the velocity of the fluid is not high enough to overcome the drift velocity, and there is no region void of charge. Points Cand Dare on the upper branch. Here, the electric charge distribution in the bulk is mainly controlled by the velocity of the fluid. As a consequence, regions void of charge are formed in the bulk. Our simulations show that the supercritical branch only exists for small values of . We can define three critical values of parameter T:Tc1is the linear criterion to start the motion of the liquid; Tc2the critical value at which the system jumps from the supercritical to the upper branch; Tfis the non linear criterionatwhichthemotionoftheliquidstopswhenwefollow down the upper branch. Figure 11 plots these critical values for severalvaluesof.WehavenotfoundanyvaluesofTc2beyond =0.125. The bifurcation is subcritical for higher values of the radius ratio and supercritical for smaller ones. In all cases, the values of Tfare smaller than the corresponding Tc1,asit is expected from the bifurcation structure. 023106-9