Free Surface Flows in Electrohydrodynamics with a Constant Vorticity Distribution

In 1895, Korteweg and de Vries (Philos Mag 20:20, 1895) studied an equation describing the motion of waves using the assumptions of long wavelength and small amplitude. Two implicit assumptions which they also made were irrotational and inviscid fluids. Comparing experiment and observation seems to suggest that these two assumptions are well justified. This paper removes the assumption of irrotationality in the case of electrohydrodynamics with an assumption of globally constant vorticity in the fluid. A study of the effect of vorticity on wave profiles and amplitudes is made revealing some unusual features. The velocity potential is an important variable in irrotational flow; the vertical component of velocity takes place of this variable in our analysis. This allows the bypassing of the Burns condition and also demonstrates that waves exist even for negative values of the vorticity. The linear and weakly nonlinear models are derived.


Introduction
Water waves constitute a very classical problem in hydrodynamics [7]. This problem is traditionally formulated in terms of the velocity potential to achieve some simplifications. In other words, there has always been an implicit assumption of zero vorticity in the flow region. In numerous recent studies, this assumption started to be dropped and the assumption of constant vorticity in the flow region used. One of the pioneering studies was made by Burns [6], and later, Da Silva and Peregrine [8] studied steep and steady waves on finite depth with constant vorticity. More recently, this problem was analyzed mathematically in some two-component systems [12]. The effect of the vorticity on travelling wave solutions (solitary and cnoidal) was investigated in the purely hydrodynamic context in [11] using the qualitative phase space analysis methods. and a Hamiltonian formulation has been reported in [20]. This problem in electrohydrodynamics seems to be still open to the best of our knowledge and the present study should be considered as a further attempt to fill in this gap in the literature.
The current approach to examining flows with constant vorticity in two dimensions is via the use of a stream function, ψ, and its harmonic conjugate, the velocity potential, and ϕ, so u = ∇ϕ +∇ ⊥ ψ. This approach introduces two unnecessary functions which complicates the problem and has the limiting effect in being restricted to fully nonlinear and linear computations. There has been no attempt to undertake a weakly nonlinear analysis which is the purpose of this manuscript. Previous work on constant vorticity models typically has been fully nonlinear; for example, see [19].
The present manuscript is organized as follows. The problem is formulated in Sect. 2. The linear analysis of this problem is performed in Sect. 3, while the weakly nonlinear analysis is presented in Sect. 4. Some numerical predictions of the weakly nonlinear theory are presented in Sect. 5. Finally, the main conclusions and perspectives of this study are outlined in Sect. 6.

Formulation
A two-dimensional fluid in region 1 is considered which is incompressible and inviscid. The vorticity, ω, is constant as is the surface tension σ . Cartesian co-ordinates are introduced, as shown in Fig. 1. Region 1 is defined as −h < y < η(t, x) ∀x ∈ R. The moving pressure distribution P(t, x) is chosen to act along the interface y = η(t, x) and P → 0 as |x| → 0.
In region 2 defined by {(x, y)|x ∈ R, y > η(t, x)}, there is an electric field, E which has no charges and is, therefore, obtained by a potential E = −∇V . The potential is chosen as V (x, −h) = 0 and as the fluid is perfectly conducting this also means that V (x, η(t, x)) = 0. A vertical electric field is set up by imposing: The equation for the electric potential is, therefore, given by: with the condition: which is one of the general boundary conditions derivable from Maxwell's equations. In region 1, the Navier-Stokes equations are used with the stress tensor: where: Here, p is called the electric permittivity. The tensor i j has various names, in the fluids literature, it is known as the Maxwell-stress tensor. It can be shown that: Therefore, the Navier-Stokes equations reduce to the Euler equations: where u = (u, v) is the velocity and P is the pressure in the fluid. The boundary y = −h is taken as impenetrable, so v(x, −h) = 0. The fluid is incompressible and has constant vorticity, ω, and so: These equations can be cross differentiated to obtain a single equation for v: The function v will take the place of the velocity potential when the KdV equation is derived for irrotational fluids. The free surface equation is given by: This gives a boundary condition for v on y = η(t, x). The lower boundary for v is given by: Equations (2.9)-(2.11) are the core of the technique in deriving the free surface profiles. The other boundary condition used is the Young-Laplace equation given by: (2.12) The Euler equations may be simplified using electrostatic equilibrium. The equilibrium condition is: Integrating this equation shows that p = −ρgy + C To compute the value of C, use the Young-Laplace equation to see C = P a − d E 2 0 /2, where P a is the atmospheric pressure. Therefore, now, write: (2.14) The Euler equations now become: The Young-Laplace equation becomes: (2.16)

Linear Theory
The scaling for the linear theory is: (3. 2) The equations become: Here: Expand according to:û The set of linear equations now becomes: (3.23)

General Dispersion Relation
To obtain a dispersion relation set P 1 = 0 and write all perturbations in the form: Solving the equation for v 1 and V 1 shows that: Only p 1 at the surface is required, so it is possible to set y = 0 in Eq. (3.19) to find: Inserting everything into the linearised Young-Laplace equations shows: The phase velocity, c, is given by c = ξ/k and solving for the phase velocity shows that: It can be seen that setting = 0 reduces to the dispersion relation in [14] and setting E b = 0 results in the dispersion relation in [15]. It was noted in [14] that to have a linear wave profile the parameters had to satisfy the inequality 4B ≥ E 2 b ; however, with the inclusion of positive vorticity, this is no longer the case.
on the branch for which c(0) > 0. One can also show that for large k, c(k) ∼ √ k as k → ∞ which shows the existence of a minimum in the dispersion relation. The beginning point at k = 0 can be seen to satisfy the equation: (3.31)

Free Surface Profiles
Consider a moving pressure distribution moving with non-dimensional speed U . Then, a frame of reference moving with speed U is selected, so all time derivative terms may be dropped. The horizontal velocity component is expanded as: The equations which are changed are then: The method of derivation is very similar to that of the derivation of the dispersion relation will be omitted. The perturbations will be expressed at: The free surface is given by: Choosing a U below the minimum of dispersion relation, one uses Eq. (3.37). To compute the waves for which U is above the minimum of the dispersion, one must include Rayleigh viscosity μ in the following way: For 0 < μ < 1, it gives waves of the form (Fig. 4a). Figure 4a, b show the free surface profiles under a moving pressure distribution for B = E b = 2 and = 1 .

The Weakly Nonlinear Free Surface
To obtain a weakly nonlinear model of the phenomena, scale according to: where c 0 = √ gh. The scaled equations are then: where = hω/c 0 . The Young-Laplace equation becomes: where: The term, F E , is the ratio of a velocity to c 0 , which shows that there is a natural velocity occurring which is given by, U = d E 2 0 /ρ; for this reason, F E should be referred to as the electric Froude number. The next step is to make the transformation: To obtain the required PDE for the free surface, the KdV scaling is required. This is α = β = ε 1. The speed of propagation, c, of the (linear) waves is unknown at this point; the following co-ordinate transformation is used: (4.6) Dropping bars and hats, the equations then become: (4.14) It can be noted that the combination of Eqs. (4.12) and (4.13) can be combined into: Expand according to: (4.21) The O(1) equations are: with boundary conditions: The equation ∂ 2 y v 0 = 0 has solution: Setting y = 0 shows that A 0 = −c∂ X η 0 , and so, u 0 = cη 0 . The equation ∂ y p 0 = 0 shows that p 0 = η 0 . The equation −(c + y)∂ X u 0 − v 0 = −∂ X p 0 can be evaluated at y = 0. Using the previous solutions yields: which gives the following expression for c: (4.26) The usual way to obtain this expression is to evaluate the Burns condition, which, in this case, is evaluating the integral: (4.27) The method presented here bypasses the evaluation of increasingly complicated integrals with simple substitution. The next order equations are: To keep the electric term in the Young-Laplace equation, a scaling has to be made on F E of the form, F E =F E ε 1/4 . To progress, one finds the expression for p 1 by integrating (4.31) and using (4.34) to obtain: To obtain the required equation for the free surface, the function v 1 (T , X , 0) is needed. Solving Eq. (4.32) for v 1 shows that: Setting y = 0 in Eq. (4.30) allows A 1 (T , X ) to be computed: showing that: Therefore, this gives v 1 (T , X , 0) to be: To compute V 1 , use the Fourier representation to obtain: and to obtain: Rη 0 e −|k|y e ik X dk. (4.41) Then: where H (·) denotes the Hilbert transform. From the properties of Hilbert transforms: Inserting v 1 into Eq. (4.33) yields the equation of interest: This is the celebrated Benjamin equation, first discovered in [4]. Travelling wave profiles are examined by η = η(x − V t) to get the equation: (4.45) where F = V /c 0 . Setting = 0 which it follows that c = 1 reduces to the equation in [14]. In the paper [15], a weakly nonlinear equation for constant vorticity is derived comparable to the one with E b = 0.
The derived equation (4.44) belongs to the family of the so-called Benjamintype (integro-differential) equations [4]. Its localized travelling wave solutions are described by the ordinary differential equation (4.45) with one non-local term. The travelling wave solution to the classical Benjamin equation was first analyzed in [5]. The stability of its solitary waves was recently studied in [9].

Numerical Methods and Results
To numerically generate the solitary wave solutions to Eq. (4.45), the classical Fouriertype pseudo-spectral discretization is employed on a sufficiently large domain. The solution decays below the machine accuracy to annihilate the effect of implied periodic boundary conditions. In other words, formally, a periodic BVP is solved, but the repeated value is actually zero in agreement with the decaying properties of solitary waves. The discrete problem for spectral coefficients is solved using the classical Petviashvili iteration as it was described in [10] for the classical Benjamin equation.
To implement the Petviashvili scheme, Eq. (4.45) is written in the following form: where the linear L and nonlinear N (·) operators are defined as: Then, the Petviashvili iteration for Eq. (5.1) reads: taking into account that the mapping N (·) is a homogeneous function of degree two of its argument. Finally, the stabilizing factor γ n is defined as: For more details on the Petviashvili iteration, refer [3,10].

Numerical Results
Equation (4.45) was discretized in space using the standard Fourier-type pseudospectral method; N = 1024 modes were used in the computations. Localized travelling waves are of interest to the study, and the computational domain was taken to be T = [−10, 10], assuming periodic boundary conditions ensured by the choice of basis functions. The Petviashvili iterations were stopped when the L ∞ norm of the difference of two successive iterations became smaller than = 5 × 10 −15 . The initial guess was a localized bump of negative polarity. The convergence of the method was marginally dependent on the choice of initial guess, which only influenced the total number of iterations. In any case, from the end user perspective, the code ran virtually instantaneously. It was noticed that the number of iterations increased with the electric Froude number F E . The computation of an oscillating travelling wave for F E = 1.27 took 117 Petviashvili iterations, while, for F E = 0.5, the method needed only 59 iterations to converge. The dimensionless physical and numerical parameters used in this computation are reported in Table 1. At this point, it is pertinent to mention that for some values of parameters, it was possible to compute the periodic travelling wave solutions, which was not the initial goal. One such solution is reported in Fig. 5. Even if this solution appears to be a single Fourier mode (such as a sine wave), the Fourier power spectrum shown in the bottom panel of Fig. 5 shows that it is actually a superposition of infinite number of modes with exponentially decaying amplitude. There are only a finite number of modes have  Table 1).
The dependence of the localized travelling wave solutions (such as depicted in Figs. 6, 7) on the vorticity parameter is also studied. For this purpose, a series of computations were performed with varying parameter between its lowest positive value for which the travelling solution exists and the maximal value fixed here to be = 5.0 . The result is presented in Fig. 8. The dependence turns out to be nonmonotonic. Though, for high values of the vorticity parameter, the positive amplitude seems to decrease and the negative amplitude seems to increase monotonically. To illustrate the shape of solutions at two extremes on Fig. 8c, a typical solution near the lowest positive vorticity value is depicted in Fig. 9a, while the solution corresponding to the highest value is depicted in Fig. 9b.

Moving Pressure Effect
In the previous section, the moving pressure effect was set to zero to concentrate on the free wave propagation. In this section, the effect of the moving pressure distribution is studied. To fix the ideas, consider the following localized pressure distribution:   Table 1 This additional ingredient poses a slight technical problem, however: the nonlinearity N (η) ceases to be a homogeneous function in η, since the pressure term has the degree zero, while the rest of terms has the degree two. Consequently, this addition forces the algorithm to switched to the extended Petviashvili method [2]. The convergence of the extended Petviashvili method is slower than in its classical counterpart. Indeed, for zero vorticity case (i.e., = 0) the iterative procedure converges up to machine accuracy with about ∼ 200 iterations (cf. ∼ 60 iterations above). When the vorticity effects are included, the extended method takes more than 2500 iterations to converge to the same accuracy. In the same time, the classical method diverges. All parameters   Table 2 used in this study are reported in Table 2. The fully converged results of the numerical computations are shown in Fig. 10, where profiles were obtained for two different values of the vorticity parameter . It was not possible find the localized solutions in this regime. Thus, the periodic profiles to which the extended Petviashvili method converged. These solutions appeared already in the free wave propagation regime, cf. Fig. 5. In particular, it may seen on Fig. 5 that the negative vorticity seems to reduce the effect of the moving pressure. To check this preliminary conclusion, a series of additional experiments was performed for several values of the vorticity parameter and moving pressure amplitude a . The result of these computations is shown in Fig. 11. In particular, it can see that the trough amplitude is essentially a linear function of a (except at small values of this parameter). This conclusion is intuitively clear.

Conclusions and Discussion
This paper has not been the first to study the properties of the Benjamin equation, and previous works such as [1,16,18] where the emphasis was usually been on irrotational flow. It has, to the authors knowledge, been the first to study it in reference to constant  Table 2 vorticity. The goal of the paper is to study the affect of vorticity on waves rather than a study of the properties of the Benjamin equation itself.
In the present manuscript, the propagation of free surface electrohydrodynamic waves in the presence of non-zero, but constant vorticity distribution was investigated. The problem was analyzed from the linear and weakly nonlinear points of view. The linear analysis allowed us to get rid of Burns's condition. The weakly nonlinear approach allowed us to compute solitary wave solutions by solving numerically the non-local ODE which describes them. The non-local effects are described by a linear term involving the Hilbert transform of the free surface excursion derivative η . It turns out that the dynamics of weakly nonlinear electrohydrodynamic waves is described by a generalized Benjamin equation [4]. So far, it appeared as a model equation for internal capillary-gravity waves [4,5]. In our study, it serves to predict the shape of coherent structures in electrohydrodynamic flows with constant vorticity. We notice that the authors studied earlier the effect of viscosity on the electrohydrodynamic waves in [13].
Concerning the perspectives, in future works, we would like to consider internal waves with differing constant vorticities and more general vorticity distributions. The equation: allows the vorticity, ω, to be some a priori prescribed function of t, x, y. Finally, the unsteady simulations have to be performed to understand better the dynamics of solutions discussed hereinabove. We suspect also that the derived Benjamin-type Eq. (4.45) possesses also other types of travelling wave solutions such as multi-pulsed solitary waves which were computed in [10] in the context of internal waves.