Electromagnetic and Thermal Phenomena Modeling of Electrical Discharges in Liquids

Abstract

Electrical discharges in liquids have received lots of attention with respect to their potential applications in various techniques and technical processes. Exemplary, they are useful for water treatment, chemical and thermal processes acceleration, or nanoparticles production. In this paper the special utility of discharges for cold pasteurization of fruit juices is presented. Development of devices for its implementation is a significant engineering problem and should be performed using modeling and simulation techniques to determine the real parameters of discharges. Unfortunately, there is a lack of clear and uniform description of breakdown phenomena in liquids. To overcome this limitation, new methods and algorithms for streamers propagation and breakdown phase analysis are presented in the paper. All solutions were tested in “active area” in the form of liquid material model, placed between two flat electrodes. Electromagnetic and thermal-coupled field analysis were performed to determine all the factors that affect the discharge propagation. Additionally, some circuit models were used to include the power source cooperation with discharge region. In general, presented solutions can be defined as universal and one can use them for numerical simulation of other types of discharges.

1. Introduction

Electrical discharges generated in liquids have many potential advantages and can be used in some low-temperature plasma techniques. Most popular utilities of such a devices are connected with water treatment [1,2], sterilization, and pasteurization of liquid foods [1,3,4] and surface treatment [1,5]. Basic benefits of HVED (high voltage electrical discharges) technique are reduced process time and temperature, minimal degradation of thermosensitive compounds, and minimization of energy consumption in comparison to classic treatment methods [4,6]. A very promising utility of HVED in agri-food industry is microbial inactivation as a result of discharges in liquid media [7,8]. The exact mechanism of the influence of the electric discharges on the inactivation of microorganisms has not yet been fully described [4,9]. Nevertheless, it can be concluded, that there are some important factors leading to a reduction in the number of microbials. These factors include production of ozone, hydrogen peroxide, hydroxyl and superoxide free radicals [4,9]. Additionally, UV radiation generated in the discharge, mechanical stresses, and waves can help reducing some microbials [9,10]. Electric breakdown process results from liquid sample ionization under high voltage applied between electrodes. Generation of the discharge process can be divided into three phases: the pulse initiation, pre-breakdown phase (streamers), and breakdown stage (arc formation) [11,12].

The characteristics of the shape of the electric discharge and the discharge propagation of processed liquid are not recognized in a satisfactory manner. Because of the varied physical parameters of liquids (juices in cases analyzed by authors), it is not possible to use universal devices and programs to implement this technology. Additionally, discharge duration, shape, and other parameters depend on electrodes system geometry and electrical parameters of power source [10]. The development of HVED technique is possible through the usage of modeling methods, facilitating the selection of discharge parameters to a specific demand. There are currently no models that allow precise description of the development of the discharge in liquid environments [13]. Main aim of the work is developing numerical models and procedures to compute pulse discharges in liquids. Mentioned task requires to account models of streamers propagation, electrical breakdown (arc) stage, and models for thermal energy propagation in the time between next discharges. Authors’ calculation procedures that use field and circuit models are described in the paper. Despite of some simplifications, presented models and procedures enable to effective modeling of electrical discharges distributions in liquids.

2. Discharge Models in Liquids

When a high voltage of appropriate value is applied to an electrode system filled with dielectric material, filamentary discharges called “streamers” appear at the first stage of breakdown in gases, liquids, and solids. Environmental conditions of the discharges determine the streamers structures including their shape, speed of propagation, current values, electric field intensity at the streamer’s head, plasma channel radius, and degree of ionization [14]. According to the large number of factors that have influence on electric breakdown process, uniform theory for the development of discharges in liquids cannot be presented. The discharge process depends on the type and purity of the liquid material, electrostatic pressure, viscosity, conductivity, thermal conductivity, and dielectric permittivity [14,15,16]. Therefore, a large number of theories have been proposed [12,14,16]. Discharges in liquids are similar, under certain conditions, to discharge process in gas [16]. In gaseous environments most important factor for discharge process is impact ionization that depend on kinetic energy of electrons in the discharge area [17,18]. In most cases it has been assumed that the discharges can be initiated by the absorption of radiation quantum [19]. These assumptions are generally acceptable for gaseous environments [14,19]. In liquids, however, such correlations cannot be found so far, because of the difficult performance of studies of the mentioned effects [20]. Therefore, the proposed descriptions are burdened with greater uncertainty. The most popular descriptions of discharges development, important from their applicability point of view, are presented below.

2.1. Direct Impact Ionization

Description of the pre-breakdown (streamer) phenomenon in liquids, resulting from direct impact ionization, is similar to the electrical breakdown process in gases [13,14]. Nevertheless, because of the higher liquid densities, electrons free paths are significantly shorter, which also affects the increased electron beam scattering [16]. Acceleration of free electrons to reach required ionization energy in a short time and free paths is difficult in such conditions and few requirements must be fulfilled. The most important is to provide appropriate voltage rise rate on the electrodes. It allows rapid acceleration of free electrons. As shown, among others, in [14,16], this process requires extremely high intensities of the electric field at the level of 220 MV/cm. Pulse rise time should be maintained at picoseconds level [20,21]. The second possibility of impact ionization is to create areas characterized by reduced mass density. In such areas (gas bubbles), electrons free paths between next collisions are longer. In addition, in the areas of lower density, the electrical breakdown voltages are reduced [6,14,16].

2.1.1. Theory of Gas Bubbles

Numerous studies [22,23,24,25] have shown that the liquid phase may be unstable during the impulse discharges. Plasma channels are formed along the electric field gradients [6,16]. Such channels, because of their lower density and different dielectric permittivity, allow easier development of discharges (especially local discharges) leading to the electrical breakdown inside the liquid [16,26]. Models that use gas bubbles theory usually assume the presence of micro-bubbles in non-gaseous media or the self-ignition of gas bubbles. The generation of gas bubbles in liquids may result from local overheating of the medium because of Joule effect in areas characterized by high intensity of electric field (near electrodes), electrostatic expansion of micro-bubbles occurring in liquid or chemical processes [20,27,28]. The process of gas bubble development has been schematically shown in Figure 1. However, in fast discharges (aprox.10 ns), characterized by fast pulses rise times (approx. 150 ps), gas bubbles are not observed. As it was mentioned in [26], gas bubbles can exist in such conditions, but their dimensions are lower than the visible size. It has been generally accepted that the theory of gas bubbles can be used to describe discharges with significant durations, larger than microseconds [16].

2.2.2. Electric Field Dependent Molecular Ionization

Electric field-dependent molecular ionization (field ionization) does not require, in general, the assumption of free electrons and ions occurrence in the space before discharge propagation. In the environment characterized by high values of the electric field intensity, it is possible to separate electrons from inert particles [16,29]. As a result of this process, free electrons and positive ions appear in the inter-electrode space (right panel of Figure 2). Because of the different mobility of electrons and ions, local variations in density of mentioned particles may occur, leading to the generation and propagation of streamers. Such models are commonly used to describe ionization in transformer oils [30,31] and water [32]. The basic disadvantage in the practical utility of this theory is the lack of analytical description of the ionization process in liquids. In practice, Zener models are in use, but they are originally used to describe the tunneling process in semiconductor materials [33].

The ionization process under the influence of the field depends on many factors, such as ionization potential, liquid density, or electric field intensity gradients [16,29]. This mechanism is characterized by many advantages, such as the lack of necessity to take into account the free energy carriers or pollutants (with lower ionization energy) during the initiation of the discharge. However, it should be mentioned that this type of ionization dominates in areas with very high electric field strengths, such as in the vicinity of blades or leader’s head [16,34]. Because of the lack of theoretical description of this process for liquid environments, utility of the proposed model is significantly limited.

2.2.3. Electric Field Dependent Ionic Dissociation

Formally, the ionic dissociation process is similar, in terms of the description, to the above mentioned ionization under the influence of the electric field [16]. To describe this process, one has to assume “active area” filled with a mixture of neutral particles, positive and negative ions. Under the influence of the external electric field, the neutral particles are dissociated, which leads to an increase in the number of ions (Figure 3) [35,36]. This fact reduces the breakdown voltage and provides the development of the discharge. It should be mentioned, however, that the mobility of positive and negative ions is negligible compared to the rate of electrical discharge development. For this reason, utility of this process to describe the development of the discharge is limited [16,35].

2.2.4. Electrostriction Mechanism

The electrostriction mechanism is used to describe discharges characterized by short duration. In the vicinity of electrodes, an inhomogeneous distribution of the high intensity electric field is created. Under conditions of significant gradients of the electric field intensity, in poorly conducting liquids, there are stresses causing deformations and discontinuities of the liquid [37,38]. In these areas, micropores filled with gas are formed, featured by a characteristic extension toward the electric fields force lines. In these structures, the free electrons can obtain energy values that enable the initiation of the collision ionization process [16,38,39]. The forces affecting the liquid as a result of the electrostriction phenomenon are described by the formula (1) for non-polar and (2) for polar liquids [16,38].

$$F\approx \frac{{\epsilon }_0}{2}\frac{\left(\epsilon -1\right)\left(\epsilon +2\right)}{3}\nabla E^2,$$

(1)

$$F\approx \frac{\alpha \epsilon {\epsilon }_0}{2}\nabla E^2,$$

(2)

where: ε₀—vacuum dielectric permittivity; ε—dielectric permittivity of the medium in which the discharge occurs; α—a coefficient determined empirically for polar liquids

The analysis of the given equations shows that the electrostriction mechanism is intensified in the conditions of fast-changing electromagnetic fields with significant gradients. Therefore, this process can be used to describe discharges of nanosecond duration, in the vicinity of pierce electrodes [39,40]. For longer discharges, the hydrostatic pressure reduces the formation of micropores [40].

2.2. Material and Method

In the previous part of the paper, most popular models of electrical discharges formation and propagation in liquids are detailed. Presented models cannot be, in general, used interchangeably, and utility of specific description depends on electrical parameters of the discharge and material parameters of the region, where breakdown process occurs.

Main topic of the paper is developing universal models and calculating procedure for the complex analysis of high voltage electric discharges. Models consider the pre-breakdown and breakdown stage of the discharge. Additionally, calculations of thermal phenomena in the time between the discharges are realized.

Before development of numerical models, some tests were performed using physical model of discharge system. Based on measurements of electrical parameters of discharges, it was possible to select appropriate model of ionization. These tasks are detailed in Section 3 of the paper.

In Section 4, general concept of HVED modeling is presented. It is assumed, that the mentioned issue can be solved using coupled electromagnetic and thermal fields. Both field and circuit models are in use.

Electrostatic calculations, described in Section 5, are used to analyze the streamers propagation in models of liquids by using FEM (ANSYS) system. Electric field distribution and gas bubble models are used to determine streamers geometry in the pre-breakdown stage. Utility of additional analytical equations enables to determine speed of streamers propagations.

In the Section 6, procedure for analysis of breakdown (arc) stage of discharge is detailed. FEM (ANSYS) numerical model with plasma channel is used. Plasma conductivity and cross section of the channel are determined by analytical equations. Momentary values of voltage drop and current in plasma channel are computed using coupled circuit model (composed of resistor and capacitor) of power supply. Such solution enables to determine transient characteristics of discharge, including power supply capacitor discharge.

Section 7 presents the analysis of temperature distribution in the time between the discharges. The computational fluid dynamics procedures are used to guarantee the precise analysis of convection heat transfer in liquid environment. Based on the calculations results obtained in this stage, it is possible to determine the minimal break time between discharges.

3. Preliminary Research of Discharge Parameters in Juices

In the previous section of the paper, most popular mechanisms describing the nature of the high voltage electrical discharges in liquids are detailed. These mechanisms usually describe specific electrical discharges, classified primarily, because of the speed of their propagation (fast and slow discharges) [27,39]. In order to select the model for further analysis of electrical breakdown process in juices, basic parameters of physical model are measured. All presented measurements are performed to determine the energy values of pulsed discharges and specific discharge periods. The test stand, consisting of a high voltage pulse generator, working chamber, and a control system, is used. Sample of juice is placed between the two flat electrodes. Figure 4 shows the design of working chamber (Figure 4a) and block diagram of the test stand (Figure 4b).

The AS3-mini analyzer equipped with an oscilloscope to record voltage and current waveforms is used. During the measurements of high voltage pulses parameters, the initial voltage was kept at 20 kV. Recorded waveforms are presented in the Figure 5. Results show the oscillatory character of the voltage values of the generated pulses. Its initial amplitude depends on the applied voltage triggering the pulse.

Tests results were used to determine the shape of the generated pulses, together with all the specific quantities (amplitudes, rise and fall times). They were used to model the electrical discharge process. According to the measurements results, it was assumed, that the theory of gas bubbles can be used to describe the discharge phenomena in the analyzed juice model [14,16,39]. In pursuance of the description of this process, the initiation of the discharge occurs as a result of the gas bubble generation because of the Joule heat effect or large field intensity gradient in vicinity of electrodes. Depending on the initial voltage and energy value, it is possible to develop “slow” (propagation speed within several km/s) or fast (propagation velocity from several to several tens km/s) discharges [27,39]. Based on the measurements of current and voltage values in the tested device, the approximate propagation velocity (Figure 5) of the discharge was determined based on the current value of the rise time.

Estimated propagation speed was at 3470 m/s, which allowed us to classify the discharge as “slow.” Under these conditions, it has been assumed that the electrical discharge is a result of ionization process in the gas phase filling the plasma channel [19,41]. This fact is important in the discharge numerical modeling point of view. As it was mentioned above, the discharge parameters in liquids are determined by the value of the electric strength, the distance between the electrodes, their shape, temperature, and material parameters of the liquid (juice in analyzed case). When appropriate number of free electrons and positive ions occur in the discharge area, one or more plasma channels (streamers) can be created. Almost all electrons will flow through the channel characterized by minimal resistance value [40,41]. Streamers starts at the cathode toward the anode. When the anode area is reached by plasma channel, the impedance of the discharge area decreases dramatically, which leads to a significant increment of the current. This stage, called the breakdown, is characterized by the occurrence of a plasma channel with almost cylindrical cross-section and duration depending on the capacity in the system [28,42]. The discharge energy can be simply defined as:

$$W=\frac{C}{2}\left(U_i^2-U_f^2\right),$$

(3)

where: U_i—voltage before breakdown, U_f—voltage after breakdown, C—capacitance

The formula was used during determination of energy balance of the discharge presented in the next section of the paper.

4. The Concept of Discharges Modeling in Liquids

Discharge model, designed to analyze the phenomena occurring in devices using HVED for the preservation of liquid food products, should enable analysis of coupled electromagnetic and thermal fields. Primarily, the model was used during the design of the laboratory discharges generator and planning future tests to implement this technique in industry. In particular, the model should allow:

-

Selection of the electrodes system in terms of their dimensions and distance;

-

Determination of the voltage values to initiate the discharge;

-

Determination of discharges energy;

-

Determination of the duration of discharges;

-

Determination of the temperature distribution after discharge in order to analyze the time required to initiate next discharge. This task was aimed to eliminate the possibility of subsequent electrical discharges in the same channel;

-

Determination of the power supply circuit parameters based on the required voltage and energy of discharges.

Characterization of all above mentioned parameters requires advanced models and calculating procedures using both field and circuit models. Field computations, carried out in the commercial FEM system, were used to compute the electrodes system variables, streamers propagation scheme, and temperature distributions. Power source parameters have been determined using circuit models. The calculation algorithm is shown in Figure 6. Both field and circuit models have been extended with Author’s procedures and algorithms, characterized in subsequent section of this paper.

5. Modeling of the Discharge Development

Numerical calculations of discharges in liquids are divided into three basic stages, allowing for the implementation of the modeling objectives described in the previous section:

-

Determination of the streamers propagation (pre-breakdown);

-

Determination of discharge breakdown (arc) stage parameters (voltage, current and energy values);

-

Modeling of the discharge system including the power source to select appropriate capacitor used during discharge process.

The analysis presented in this paper uses a model consisting of two flat electrodes of 50 mm diameter. Active area (model of liquid juice), wherein the discharge occurs, was placed between the electrodes. Distance between the electrodes was 60 mm. The geometry of the model is shown in Figure 7a.

The gas bubble theory was used to determine streamers propagation in the active area. This assumption was introduced, because of the results of preliminary analysis presented in previous section of the paper. According to the obtained results, one can classify the discharge as “slow” [39] and theoretical description of gas bubble theory is applicable in such conditions [14]. The static models implemented in ANSYS program were used to determine electric field distributions. According to the proposed procedure, streamers development was computed using a special “virtual” surfaces arranged in parallel to electrodes (Figure 7b). Adopted “virtual” surfaces form a kind of differential division of the active area. All streamers heads can be located only on these surfaces and streamers spread are computed using a “frame by frame” method. So that, resolution of calculation results depends directly on the number of surfaces. The calculation procedure has been detailed further in the paper. Computations of streamers development could not be performed using only commercial FEM system (ANSYS). Some additional procedures were implemented in Authors’ subprograms. Mentioned procedures were used to compute electric field gradients, voltage values at streamers heads, and to generate geometry of streamers paths and gas bubbles.

According to the description of the HVED technique, a series of several tens to several hundred discharges is used [4,9]. Therefore, the procedure for determining the streamers propagation have to be performed many times, regarding the number of analyzed impulses. The initial conditions for each impulse, in the form of variable physical parameters of the liquid environment, have been introduced based on the temperature distribution that occurs after the previous pulse was generated (Figure 6). During the modeling of the first discharge, it was assumed that the active area is characterized by homogeneous temperature, conductivity, and concentration of charges. The solution of the electric field intensity equation and the electric potential (4) for the electrostatic problem in such a system always leads to results, where the electric field intensity is homogeneous throughout the analyzed area and the electric potential distribution is characterized by iso-surfaces parallel to the electrodes (Figure 8a).

$$div\left(\epsilon E\right)=4\pi \rho$$

(4)

where: E—electric field intensity; ε—dielectric permittivity; ρ—charge density

Because of homogeneous electric field distribution, discharge initiation point cannot be clearly determined. In order to initiate first discharge, the random number generator was used, implemented in a MathCAD [43] program. Node numbers (of the finite element mesh) lying on the cathode surface were exported to this software. Random node number was generated and re-entered into the ANSYS program to initiate the discharge. The gas bubble model with spherical geometry of 10 μm diameter was created around the selected node (Figure 8b). According to the Equation (5), the intensity of the electric field in the bubble is higher than in its direct vicinity [39]. The phenomenon was considered during calculations by introducing the value resulting from the formula (5) as an additional source of electric field in the numerical model (Figure 8c) [17].

$$E_B=\frac{3{\epsilon }_r}{2{\epsilon }_r+1}{E}_{ia},$$

(5)

where: E_ia—the intensity of the electric field in the vicinity of the gas bubble; ε_r—electrical dielectric permittivity

The plasma channel can be developed only if the local electric field intensity on the streamer’s head exceeds the critical value specified for breakdown of the dielectric material (juice model in the analyzed case) (E *) [17]. This condition was used to determine the streamer’s propagation path, by determination of electric field intensity gradients (Figure 9a). However, in many cases, the main discharge streamer can be branched. The condition for creation of new streamers may be presented basing on the field fluctuation criterion, leading to formula (6) [17].

$$E_i>E_{\ast }-{\delta }_i,$$

(6)

where: E_i—the gradient of the electric field intensity in the direction in which a new streamer can propagate; δ_i—random variable, taking into account the fluctuations associated with the local heterogeneity of the critical values of the breakdown voltage (randomness of the discharge).

Computations of streamers propagation scheme were performed by the analysis of electric field intensity gradients distributions between individual nodes located on the successive “virtual” surfaces of the discharge area (Figure 7). In cases, where gradients were characterized by similar values (differences lower than 10%) (Figure 9b), the condition (6) was analyzed and a new streamer was introduced (Figure 9c,d).

The proposed procedure uses static calculations. Considerable disadvantage of this solution is the inability to use it to transient analysis of the pre-breakdown stage. In order to enable such calculations, the Equation (7) [17,44] was used to determine the average speed of streamers propagation, based on the analysis of local values of the electric field intensity.

$$v(E)=\frac{h}{\tau }{e}^{-\left(\frac{E_{\ast }}{g}\right)}{e}^{\left(\frac{E}{g}\right)}$$

(7)

where: τ—time step; h—distance (between subsequent surfaces of the differential division of the discharge space); g—width of the streamer; E—electric strength, E_*—breakdown field strength

Equation (7) was used to determine the time, after which, the streamer head reaches subsequent “virtual” surfaces dividing the area of the discharge propagation (Figure 7b). Such computations can be performed for all calculation steps to determine the temporal characteristics of the discharge propagation path.

Procedure shown in Figure 9 was used to determine the locations of all nodes, where streamers heads occur. Around these nodes, models of the gas bubbles were created, similarly to the above described modeling of the first node located on the surface of the electrode (Figure 8b). The occurrence of the plasma channel between the gas bubble models is taken into account. The channel was simulated in the form of a cylinder with a diameter of 10 μm (Figure 9c,d). Inside the channel, the model of the plasma of conductivity determined by Equation (8) was assumed [45,46,47]. This formula is useful for poorly ionized plasma. During calculations movements of plasma particles was neglected [46].

$${\sigma }_M=0.85\frac{n_e\cdot e^2{\left(8kT\right)}^{1/2}}{m_e{\left(\pi m_e\right)}^{1/2}\cdot n_o\cdot s_a},$$

(8)

where: T—temperature, k-Boltzmann constant, m_e—electron mass, n_e—concentration of electrons, n_o—concentration of electrically neutral particles, e—atomic unit of charge, s_a—active cross-section in which collisions occur

Relating to Equation (8), it is necessary to know the plasma temperature to determine its conductivity. In general, the presented algorithm for streamers development stage, does not require to perform thermal calculations. This fact results directly from the analysis of bibliography [39,44,47], where it has been proved that in pre-breakdown stage (before electrical breakdown), the value of dissipated energy is relatively small in comparison to the energy released during the breakdown (arc) stages of discharge. Using the equations for isothermally heated body (9), it was possible to develop an approximate dynamic temperature characteristic (Figure 10a) for the discharge propagation. The obtained temperatures were used to compute the non-linear value of plasma conductivity (Figure 10b) [48]. The proposed solution is simplified because of the assumed constant power value during the streamer conductivity computations. To consider mentioned effects it is necessary to perform iterative calculations.

$$t\left(\tau \right)\approx P_p\cdot \frac{\tau }{c_{pl}\cdot m_{pl}}+t_p,$$

(9)

where: Pp—power dissipated in plasma; c_pl—specific heat of the plasma; m_pl—plasma mass; t_p—initial temperature

Plasma conductivity values were used to determine the charge distribution (10). All steps described above were used to create the computational algorithm for the propagation of electrical discharges in liquids, as shown in Figure 11.

$$\frac{\partial \rho }{\partial \tau }+div\left({\sigma }_M\cdot E\right)=0,$$

(10)

The electric potential distributions in the model for the first discharge have been shown in Figure 12. In terms of quality, the results confirm the proper operation of the algorithm.

The algorithm allows to calculate the density of currents in streamers and the power loses during the initiation stage of the discharge. Current densities in different time steps are presented in Figure 13.

One can use the algorithm, to perform full calculations of the pre-breakdown, both in geometrical and energy point of view.

6. Modeling of the Breakdown (Arc) Stage

Further analyses were carried out for the breakdown stage, where the plasma channel occurs between the electrodes (channel “connects” the electrodes). It was assumed, during this stage, that all “secondary” streamers (streamers that not reach the electrode surface) were neglected (Figure 14). To compute the current and voltage values in the time domain, it was necessary to determine the electrical conductivity of the plasma, where charge carriers are in the form of electrons and ionized particles. It was required to establish the influence of thermodynamic parameters (mainly temperature and pressure) on the conductivity value. For highly ionized plasma, the conductivity value can be computed using formulas (11), and (8) for weakly ionized plasma [39,46].

$${\sigma }_D=0.591\frac{{\left(4\pi \epsilon \right)}^2{\left(kT\right)}^{3/2}}{e^2{m}^{1/2}}\frac{1}{ln\theta },$$

(11)

where: T—temperature, k—Boltzmann constant, m—mass of the particle, θ—electron concentration, e—elementary energy

It was assumed that in the presented system there is a plasma characterized by medium ionization degree. Electrical conductivity of the plasma was estimated as the average value from formulas (8) and (11) [39,47]. Figure 15 presents the conductivity values as a function of the plasma temperature. Because of the mentioned relation, arc stage analysis requires coupled electromagnetic and thermal calculations.

Dimensions of the plasma channel were determined automatically. In the case of a circular arc, the equation for radius of the plasma channel can be specified as [47]:

$$r_{ł}={\left(\frac{2}{{\pi }^2}\right)}^{1/4}\cdot {\left(\frac{1}{\sigma \rho h}\right)}^{1/4}\cdot {\left(2g\frac{{\rho }_a}{\rho }\right)}^{1/8}\cdot z^{1/8}\cdot I^{1/2},$$

(12)

where: I—current, z—arc length, g—gravity, h—specific enthalpy of the arc plasma, ρ—arc plasma density

According to formula (12), to determine the plasma channel radius, it is necessary to input discharge current value and iterative calculations must be performed to use it. The values of discharge currents were introduced on the basis of the results obtained in previous iteration steps of the conducted analyzes. In the first step, the discharge diameter was assumed as 10 μm (the value was the same as in the pre-breakdown modeling). The arc plasma of the given conductivity was introduced to the numerical model in the form of a solid body (Figure 14b). The physical phenomena occurring in the plasma have been omitted. Nevertheless, this simplification is acceptable in relation to determination of active power of plasma channel. Static calculations of the electromagnetic problem were carried out to compute the discharge conduction current (Figure 16a,b) and the heat sources distribution (Figure 16c). These computations were performed iteratively, with the iteration frequency depending on the thermal analysis results discussed further in the paper.

Based on the conduction currents density (J) and the distribution of heat sources (p_V), it was possible to compute the resistance of the arc (13). The proposed method uses the total active power (P) and current (I) values, obtained by adding its elementary values in individual finite elements in the discharge model (i = 1 … N). The computed resistance values were necessary to design circuit models that allow to determine operational characteristics of discharge model equipped with power source.

$$R=\frac{P}{I^2}=\frac{{\displaystyle \sum _{i=1}^N{p}_{Vi}\cdot V_i}}{{\left[{\displaystyle \sum _{i=1}^N\left(J_i\cdot F_i\right)}\right]}^2},$$

(13)

where: P—active power value, p_Vi—Joule heat of i-th finite element, V_i—volume of i-th finite element, J_i—current density in i-th finite element, F_i—external area if i-th finite element

Thermal calculations were performed as transient, using previously computed heat sources and boundary conditions characteristic for conduction and convection heat transfer. In the example, the temperature distributions in time 0–0.001 s are shown in Figure 17.

After thermal calculations in any time step, it is necessary to change the material parameters (primarily the plasma conductivity) and the supply voltage in the breakdown stage. In current algorithm, voltage source was modeled as the RC circuit. Resistance values determined based on the field modeling (13) and the capacitance (the value assumed as equal to the power source capacitor—Figure 4) were used to compute the time constant of RC circuit and source voltage values in any time step of the analysis (14). Figure 18 shows the voltage over time characteristic, obtained by the conducted analysis.

$$V\left(\tau \right)=V_0{e}^{-\frac{\tau }{RC}}$$

(14)

where: V₀—initial voltage, τ—time

The voltage value on the electrodes, after thermal calculations, was introduced as the initial condition in the next cycle of static electromagnetic calculations. Presented algorithm can be realized iteratively until whole energy stored in the capacitor before discharge initiation is dissipated. From the moment of capacitor discharge, no current will flow in the arc channel. During the break time between next discharges, temperature distribution problem was analyzed.

7. Modeling of Temperature Distribution in the Time between Discharges

The temperature distribution after the discharge extinction, can be used to re-analyze the generation of the next discharge, carried out by using the principles given in Section 5 of this paper. For example, time required for the proposed solutions, despite the significant time consumption, enables a full analysis of discharges in liquid environments.

It is worth to point out that for simplifications of numerical analysis, two different heat transfer models were used during discharge phase and break between the next discharges analysis:

-

The “static” approach without fluid movements modeling (15) was used to compute heat transfer process during breakdown (in the pre-breakdown (streamer) stage, thermal analysis was not performed). The special criterion value of thermal conductivity (16) was used to include thermal convection heat transfer in the fluid [49];

-

The CFD model was used for heat and mass transfer analysis in the time between next discharges. Coupled equations for continuity of flow (18), momentum (19) (here presented in one dimension of Cartesians system), and energy (20) balances were computed in Ansys program. A typical turbulent k-ε model was in use. Temperature field distribution in the last step of discharge analysis was used as initial condition for CFD modeling.

$$\frac{\partial t}{\partial \tau }=\frac{{\mathrm{p}}_{\mathrm{V}}}{\mathrm{c}\rho }+\frac{1}{c\rho }[\nabla (\lambda \nabla \mathrm{t}\left)\right]$$

(15)

$${\lambda }_e=\left(1+\frac{C_{1}R{a}_s^n}{R{a}_s+C_2}\right){\lambda }_s$$

(16)

where: t—temperature, τ—time, c—specific heat, ρ—mass density, λ—thermal conductivity, λ_e—equivalent thermal conductivity (including fluid velocity), Ra_s—the Rayleigh number, C₁, C₂, n—factors depended on fluid flow type (laminar, turbulent).

$$\frac{\partial \rho }{\partial \tau }+\nabla \cdot \left(\rho \overrightarrow{v}\right)=0,$$

(17)

$$\rho \frac{d{w}_x}{d\tau }=\rho g_x-\frac{\partial p}{\partial x}+\mu \left(\frac{{\partial }^2{w}_x}{\partial x^2}+\frac{{\partial }^2{w}_x}{\partial y^2}+\frac{{\partial }^2{w}_x}{\partial z^2}\right)+\frac{1}{3}\mu \frac{\partial }{\partial x}\left(\frac{\partial w_x}{\partial x}+\frac{\partial w_y}{\partial y}+\frac{\partial w_z}{\partial z}\right)$$

(18)

$$\frac{\partial }{\partial \tau }\left(\rho h\right)+\nabla \cdot \left(\rho h\overrightarrow{v}\right)=\nabla \left[\left(\lambda +{\lambda }_t\right)\nabla T\right]+p_V,$$

(19)

where: ρ—fluid mass density; p—hydrostatic pressure; ρg—mass forces; μ—dynamic viscosity; h—enthalpy; λ—thermal conductivity; λt—thermal conductivity for turbulent flow; T—temperature; p_V—volumetric density of heat sources and loses.

Because of the fact that the electrical permittivity and resistivity values of analyzed liquid depend on the temperature value (both parameters decrease in temperature function), calculations of thermal energy propagation in the time between the discharges seems to be a very important factor [50]. If the break time between initiation of the next discharge is too short, a large heterogeneity of temperature distribution is observed in the liquid sample. In such conditions (Equations (4) and (5)), streamers will probably propagate in the channel of the previous discharge or its vicinity. Transient thermal analysis between discharges enables to determine minimal time interval, after which, next discharge will propagate randomly, independently from the pattern of previous breakdown process. Exemplary temperature and fluid velocity distributions computed after first discharge (presented in Figure 16 and Figure 17) are presented in Figure 19.

Performed CFD analysis were used to determine the electric field intensity distributions in the analyzed workload (sample of liquid juice). Temperature distributions were introduced as initial conditions of the electrostatic problem. Such approach enables to consider the temperature dependence of material parameters that are particularly important for electromagnetic simulations (ε, ρ). Based on the presented approach, it was possible to compute break time between the discharges, to guarantee random propagations of streamers and new discharge channel. In Figure 20 electric field intensities (a, b) and the gradients of electric field intensity (c, d) in vicinity of electrodes (1) are shown for the break time value of 0.5 s (a, c) and 1 s (b, d) from previous discharge. In the figure geometry of previous discharge channel is shown (2).

Simulations results, presented above, show that break time of 0.5 s is too short for initiation of next discharge. In this case, heterogeneous distribution of temperature in the analyzed fluid leads to the conclusion that the next discharge will propagate in the channel of previous one (Figure 20c). After 1 s, temperature distribution and electric field intensity distribution are more homogeneous, and random distribution of the discharge channel is probable in this case.

8. Conclusions

The article deals with modelling and simulation issues of pulse discharges occurring in liquid environments, used in many techniques. As it was mentioned in the introduction, there is a lack of universal and accurate models to simulate all stages of discharges (initiation, pre-breakdown, and arc). Most popular models of discharges in liquids (especially dielectric) were shortly described. Because of the lack of information about the parameters of discharges in juices, preliminary tests were performed. Based on the authors’ laboratory tests, most important electrical parameters of discharges were performed. Additionally, specific transient characteristic of discharge propagations were determined. Considering the obtained results, discharge was classified as “slow.” In such type of discharge, gas bubble theory model was used to describe the phenomena of streamers propagation paths. In the paper, the numerical method and algorithm for pre-breakdown were characterized. Subsequently, analysis of electric potential gradients in vicinity of the streamer’s heads were used to determine the streamer position. The proposed method, based on gradients determination and probability distribution, was used to establish total discharge path and time characteristic of discharge growth. The mentioned phase of discharge is described in Section 5 of the paper. It is worth to highlight that the proposed method enables to determine energy values based on momentary values of current and voltage.

In the Section 6, authors’ solutions in the scope of modelling of breakdown (arc) stage of discharge were detailed. For the given geometry of discharge path (determined in previous steps of the analysis), the total energy generated in the arc plasma was established using iterative calculations. Plasma conductivity and the cross section of plasma channel were computed using analytical formulas described in Section 6. Coupled analysis of electromagnetic and thermal fields (in FEM system) were performed to determine the total power of the arc and transient temperature characteristics. Circuit models were used to simulate the power source in the form of RC system. Different capacity values of series capacitor enable to form different arc current characteristics. The CFD models were used in the last stage, where temperature distribution in the time between next discharges was computed. Based on the results, it was possible to consider nonlinear material parameters for the next discharge distribution modelling.

Results presented in the paper concern discharges in juices and other liquids. All algorithms can be classified as original Authors’ solutions. It has to be highlighted, that similar solutions are rarely presented in scientific papers. Obviously, discussed models and algorithms, can be extended to obtain results characterized by higher accuracy. However, obtained results are promising and are characterized by high versatility.