New approach to an underwater discharge simulation using elliptic coordinates

In this study, a system of two equations is developed that allows the time evolution of the underwater spark channel to be calculated numerically from a given power input. The proposed mathematical model utilizes the elliptical coordinates. This approach has the advantage of considering the underwater spark as an expanding ellipsoid, which closely corresponds to experimental observations. Similar to spherical or cylindrical models, the proposed method considers only one spatial coordinate as a function of time, which simplifies the analysis. The predictions of this model are compared with the experimental results.


I. INTRODUCTION
An electrical breakdown in liquids is usually accompanied by a formation of underwater spark. The accelerated expanded spark channel generates a powerful acoustic emission in the form of shock waves. The shock wave generation has various industrial applications: in medicine, 1 oil and gas technology for well stimulation, 2,3 and electrical discharge machining. 4 Therefore, a detailed study of underwater discharge has a particular interest.
The efficiency of shock wave generation depends on impedance matching between the characteristic impedance of a high voltage pulse generator and the resistance of a spark channel. As a rule, the characteristic impedance significantly exceeds the spark resistance. The electrical characteristics of capacitor, the inductance of cables and the electrode system put a strong limit on minimization of the characteristic impedance. Therefore, it is desirable to increase the resistance of the spark channel by having a long distance between the electrodes. Here, long distance means that the distance should exceed the streamer thickness before the breakdown. Therefore, a mathematical simulation of underwater discharge generated in long gap attracts a particular interest. In this case, the cylindrical coordinate system is used for the analytical model and the dynamics of the underwater spark is described by the cylinder radius defined as a function of time. [5][6][7][8] The long spark channel develops rapidly into a spherical bubble. The cylindrical approach has two disadvantages: the real underwater spark is not a cylinder and the conversion from the cylindrical to the spherical coordinate system is required to simulate the full process of underwater spark growth. 6 An elliptical coordinate as a two-dimensional orthogonal coordinate system has been identified as a useful tool in many different disciplines in applied mathematics and physics such as electrodynamic applications, 9,10 insect aerodynamics, 11 study of electrostrictive ponderomotive forces in fluid, 12,13 and Kirchhoff vortex. 14 In this study, in an effort to mathematically describe the expansion process of an electrical spark in an uncompressible liquid, we study the Rayleigh method using the elliptic coordinates. The model developed in this paper provides a close approximation of the spark observed in the experiment.

II. SPARK GROWTH MODEL
The proposed model is based on the theoretical analysis of the spark expansion performed in the elliptic coordinates. The elliptical coordinates (μ, ϑ) are related to the Cartesian ones (x, y), with 0 ≤ μ < +∞ and 0 ≤ ϑ ≤ 2π. Surfaces of constant μ are confocal ellipses centered around (x, y) = (0, 0). Constant ϑ gives hyperbolas.

ARTICLE scitation.org/journal/adv
Two focuses of ellipses and hyperbolas are at (−a, 0) and (a, 0) to the plane. A set of shapes are drawn in Fig. 1, which represent the elliptic coordinate system. In the frame of the presenting model, we further assume that the trajectory of every point on the spark surface is a hyperbola corresponding to a coordinate ϑ. Subsequently, the velocity vector ⃗ v(μ, ϑ) is always normal to the elliptic channel surface and have to be considered as a function of time with respect to the coordinate μ(t) (Fig. 1).
This theoretical analysis will be simplified by neglecting the stochasticity of the streamer growth, the kink instability of the spark channel, and the presence of electrodes, by assuming that at time t = 0, an elliptical plasma channel 2a long and 2b thick exists. It is also assumed that the effect of the magnetic field on the underwater spark dynamics is negligible. For t > 0, the electrical energy is deposited to the spark and the spark channel expands due to the high temperature and pressure caused by Joule heating. The physical parameter distributions in the plasma column are assumed to be homogeneous. These assumptions allow us to consider the spark channel as a prolate expanding spheroid (ellipsoid) with the semimajor radius a and semiminor radius b. The radius a is considered to be equal to the sum of a half of interelectrode distance l and the small parameter r el , where r el reflects the processes of vapor bubble formation on the anode and a minor streamer expansion in the direction opposite to its growth from the cathode before the breakdown. The prolate spheroid is invariant under rotation, and a two dimensional model may be considered.
We use the approach and conventions by Rayleigh, who analyzed the dynamics of a spherical bubble in an infinite incompressible liquid. 5,15,16 Then, at a given instant, the rate of mass of liquid flowing through any spheroid surface Sw(μ) must be constant. In time ∆t, the mass of liquid flowing with velocity ⃗ v(μ(t), ϑ) across a surface corresponding to μ(t) is where ρ 0 is the water density. Equating this to the liquid flow through the spark wall surface Ssp(μ(t)) gives (4) where υsp(μ(t), ϑ) is the spark wall velocity. Performing differentiation of both parts of Eq. (4) by angle ϑ, we can express the fluid velocity as Since φ = −∇v, it is easy to obtain the velocity potential distribution of the liquid from formula (5), where KL is the Lame coefficient, The surface area of the prolonged spheroid 11 in formulas (5) and (6) is Substituting Eqs. (1) and (2) into the expression for the velocity v = √ Vx 2 + Vy 2 and after some transformations, we obtain the following formula for the velocity of spark channel expansion: Using Eqs. (7)-(9), the velocity potential (6) can be written as Substituting expressions (9) and (10) into the Bernoulli equation for unsteady potential flow, and determining the constant from the boundary condition at infinity after some transformation, we find the equation in elliptic coordinates for spark expansion, we receive a system of two differential equations, where the volume of the spark channel V(t) is expressed in the prolonged spheroid approximation, γ = 1.26 is the ratio of specific heats, and N(t) is the power deposited in the discharge. The volume bounded by the prolonged spheroid is Substituting Eqs. (1) and (2) into (14) and taking into account that the coordinate Rx lies on the X axis (ϑ = 0) and the coordinate Ry lies on the Y axis (ϑ = π/2), we obtain

III. EXPERIMENTAL SETUP
The experimental setup was similar to that used in experiments described in Refs. 18 and 19. Figure 2 provides its illustration. Electrical underwater discharge was generated by discharging a capacitor bank via an electrode system mounted on a laboratory tank filled with tap water. Spark was generated between two pin electrodes made of tungsten wires with a diameter of 0.8 mm, and the electrode gap was 11 mm. The inductance and resistance of the electrode system were measured with Keysight U1733C: 0.7 μH and 30 mΩ, respectively. These values are negligible and are not taken into account in calculations. The electrodes were connected to a capacitor bank Cs with a capacitance of 0.8 μF by a high voltage coaxial cable. The capacitor bank was charged up to 17.5 kV. A triggerable spark-gap switch was used to apply the accumulated electrical energy to the underwater electrical discharge. 17 Video sequences of the expanding spark channel were obtained using a high-speed Phantom v710 camera equipped with a Nikon 200 mm f/4D IF-ED AF micro-objective lens with a focal distance of 200 mm. A low camera resolution of 128 × 128 pixels was chosen with the aim of giving the maximum frame rate (210 526 fps).
The camera is positioned in front of the tank window, as shown in Fig. 2. The underwater spark was visualized using a shadowgraph method. A strong light source COOLH Dedocool tungsten light head was used for this purpose. It projects a powerful beam of light through the electrode region in the camera direction. Temporal voltage and current waveforms were measured using a PVM-1 2000:1 voltage divider (North Star Research Co., maximum DC/pulsed voltage 40/60 kV and maximum frequency 80 MHz) and a Pearson model 4997 probe (maximum peak current 20 kA and maximum frequency 15 MHz), respectively. The voltage divider is placed between the high voltage coaxial cable and electrode system. The measured signals were recorded using a Tektronix MD04054C oscilloscope. Wolfram Mathematica 11.3 software was used to model and simulate the temporal behavior of an expanding spark channel.

IV. EXPERIMENT AND SIMULATION RESULTS
The preceding discussion outlines the mathematical model that will be integrated to determine the time evolution of underwater spark. The result of calculation will be compared with the expansion data of underwater spark obtained from the experiment and shown in the sequence of eight images of spark expansion taken from supplementary material Video 1 (Fig. 3). Initial conditions had to be chosen for the variables to be integrated in the model equations. A logical choice for the initial pressure is P(0) = P∞, where P∞ is the atmospheric pressure. The initial conditions for μ(0) andμ(0) (the overdots denote derivatives with respect to time) have to be chosen carefully, taking elliptic coordinates. For this purpose, let us consider the velocity of spark expansion along the X-axis (ϑ = 0) obtained from Eq. (1), and for initial time, we obtain a * sinh(μ(0)) = b, where b is a minor radius of the ellipse determined as an initial radius of streamer just before the breakdown. Subsequently, the expression for Vμ (0) is The approximate value of radius b before the breakdown is determined from the image in supplementary material Video 1 and is approximately equal to 3 × 10 −4 m. As was mentioned above, the vapor bubble is growing on the anode. The velocity of its expansion is measured as a difference of spark channel radii observed on two adjacent frames divided by a period T f between them, The streamer growing from the cathode to the anode reveals small expansion in the direction opposite to its growth. This expansion velocity and vapor bubble growth are approximately equal, Vx(0) ≈ 8 m/s. After substitution of Vx(0) in Eq. (18), we receive the initial value Vμ = 26 667 a.u./s. It should be noted that the system of differential equations (12) and (13) is not sensitive for the measured value Vx(0). However, a zero value, as it is often assumed for spherical models, leads to the numerical results different from the experimental data. The electrical energy deposited into the spark determines its dynamic behavior and is described by the power function N(t) in Eq. (13). The power function is expressed as a product of the experimentally measured voltage and current, as shown in Fig. 4, (20) The measured current and voltage waveforms were smoothed in the Origin program with the purpose of accelerating numerical simulation.
Since the breakdown stage is not considered in the present work, we start calculations from the breakdown moment characterized by an abrupt decrease in the voltage, as shown in Fig. 4. Substituting the initial conditions and numerically solving the system of two differential equations (12) and (13), we receive the function μ(t) defined on the time interval from 0 μs to 165 μs. The upper time limit was chosen arbitrary and meets two requirements: it has to exceed the period of energy deposition and be less than the period of spark observation on supplementary material Video 1. Equations (1) and (2) parametrically represent the spark channel shape in x and y spatial coordinates through the function μ(t). Considering angle ϑ as a parameter and calculating (1) and (2) for different times t, the set of spark channel shapes was received and is shown in Fig. 5. At the initial moment 0 μs, the simulated spark is a long thin ellipse according to model assumptions. It expands more vertically than horizontally with time and approaches a circle. Substituting ϑ = 0 in Eq. (1) and ϑ = π/2 in Eq.
(2), we obtain the spark major radius Rx(t) and minor radius Ry(t) depending on time. The curves Rx(t) and Ry(t) are presented in Fig. 6.
Taking the time derivative of Eqs. (1) and (2) and after the substitution of ϑ = 0 in Eq. (1) and ϑ = π/2 in Eq. (2), we obtain the velocity of spark expansion along the coordinate axes X and Y, The temporal behavior of the calculated velocity profiles during spark channel expansion is presented in Fig. 7. The maximum expansion velocity reaches 580 m/s approximately in 2 μs after the discharge initiation in the vertical direction. The horizontal expansion is approximately five times less. After reaching maximum, the expansion velocity is gradually slowing down. Figure 8 shows the pressure in the spark channel calculated using the system of differential equations (12) and (13). The maximum pressure coincides in time with the maximum expansion velocity.
Our calculation method is verified with the available experimental data. Figure 3 shows examples of the frame series taken from supplementary material Video 1. The channel of underwater spark is black in color. The experimentally observed spark channel has a tendency to twist into the structure with a helical topology [ Fig. 3(b)]. This signifies that a current-driven kink instability has taken place. 13,14 After the current termination, the radially expanding force suppresses the helicity and the spark expands as a prolate spheroid. During this expansion process, the spark broadens out and asymptotically approaches a spherical shape. The experimental radii Rx(t) and Ry(t) were obtained from supplementary material Video 1.
The major radius is defined as with Lsp being the maximum spark length on the frame corresponding to time t. The minor radius of spark is defined as with dsp being the spark thickness measured in the middle of the interelectrode gap on the frame corresponding to time t. The measured radii vs time are presented as scatter curves in Fig. 6. The surface of the spark channel is not homogeneous and has a fine structure with ripples and folds [Figs. 3(b)-3(h)]. Therefore, some of the data points representing the spark expansion do not lie on smoothly ascending curves. During the spark channel expansion, the resultant elliptical volume develops gradually into a spherical bubble. At moment t, the semiaxes of an elliptical channel become approximately equal and transfer from the elliptical coordinate system to a spherical coordinate system. 6 Consider a new variable Rtr expressing the initial radius of the expanding bubble in the spherical coordinates. Then, In accordance with Eqs. (1) and (2), the ratio of radius Rtr to the half of interelectrode distance a is With a sufficiently high accuracy, the ratio y(t) x(t) can be chosen as 0.9 (in general, R y (t) R x (t) < 1). Then, after calculation, the initial radius Rtr for the expanding spherical bubble is 0.11 m in the case considered in this paper. The full simulation model of a cavitation bubble generated by the spark lies outside the scope of the present paper and can be found in Ref. 6. The ratio in (25) is chosen arbitrary   FIG. 9. Dependence of the initial radius of the spherical model divided by a half of interelectrode distance on a ratio between the semiaxes Ry(t) and Rx(t).
according to the required accuracy. The radius Rtr is not a function of time and depends only on the interelectrode distance. The general dependence of the ratio of bubble radius Rtr to the half of interelectrode distance a is plotted in Fig. 9.
There are some discrepancies between the modeling results and experimental data. These discrepancies could mainly be explained by the assumptions employed for our mathematical model. The difference between the experimental data and simulation results is distinct in two regions. First, the calculated semimajor radius Rx(t) is slightly less, approximately 70 μs after breakdown, when the expansion velocity is decreasing. This disparity can be explained by the violation of the initial conditions determined by Eqs. (9) and (10) for velocity υ(μ, ϑ) and potential φ(μ, ϑ) functions at the coordinates ϑ = 0 and ϑ = π. This violation is caused by the presence of the electrodes which is not taken into account in the considered model. Second, the calculated minor radius expansion significantly exceeds during the first 80 μs. A reason explains this: the magnetic field generated by the current passing through the spark squeezes the spark channel. When the current is terminated, the simulated curve Rx(t) equalizes the experimental curve and coincides with it. In the literature, the effect of magnetic constriction is supposed to be negligible. 5 Nevertheless, the magnetic field can substantially affect the spark channel dynamics 18,19 as it can be seen in supplementary material Video 1 and Figs. 3(b) and 3(c). Therefore, the factor of the magnetic field on the spark expansion has to be studied, which is not in the frame of the presented paper. Nevertheless, the considered model generally corresponds to the dynamics of the observed underwater spark.

V. CONCLUSIONS
We have developed a new mathematical model, which will enable one to numerically investigate the behavior of fast transient cavities generated in water by high voltage underwater spark discharges in uncompressible liquid. The elliptical coordinates are used for analysis. The main difference with the previous reported models lies in considering the spark channel as an expanding prolate ellipsoid. The temporal development of the spark channel ARTICLE scitation.org/journal/adv was obtained using this model and compared with the experiment. A minor discrepancy from the experimental results was revealed in simulated spark dynamics. The reasonable explanation of this lies on the effect of magnetic field on the spark channel and the effect of electrodes not taking into account in the model. The proposed mathematical model provides a reasonable agreement with experimental measurements.

SUPPLEMENTARY MATERIAL
See the supplementary material for a video of an underwater expanding electric spark.