Electromechanical coupling mechanisms at a plasma–liquid interface

Cold atmospheric plasmas (CAPs) have been the focus in recent years due to their unique ability to generate a mixture of highly reactive chemical species under ambient conditions. The number of applications where the use of CAP has been explored is growing rapidly, and promising results have been obtained spanning the domains of agriculture,^1,2 plasma medicine,^3–5 and materials processing.^6,7 In many of these applications, the discharge inevitably interacts with some form of liquid, be it moisture on freshly harvested vegetables, right through to the deliberate treatment of a liquid volume to create plasma activated liquid (PAL). Using PAL offers the intriguing possibility of retaining some chemical activity initiated by the plasma for several days,^8–10 which is a major advantage for most applications. Unfortunately, this advantage comes at a price, as the introduction of a plasma–water interface vastly increases the complexity of the physiochemical processes at play.

A variety of CAP sources have been used to activate liquids, these can be categorized based on whether the plasma is in direct contact with the liquid surface or not. In non-contact sources, the plasma is generated some distance away from the interface, and the reactive species generated by the plasma are carried to the interface by diffusion, convection, or both.^11,12 Such sources include surface barrier discharges (SBDs),^11,13 gliding arc discharges,^14,15 and plasma jets, where the plasma terminates before arriving at the interface.^12,16 In direct-contact sources, on the other hand, the plasma physically contacts the liquid surface, creating a complex gas–water–plasma interface. This leads to a higher dose of short-lived species such as ions and radicals reaching water and thus a higher activation in a shorter period of time.^17,18 Examples of such discharge configurations include pin-water discharges,^19,20 water film dielectric barrier discharge (DBD),²¹ and some configurations of plasma jets.²² 

Whether PAL is generated using a contact or a non-contact source, there is potential for the plasma treatment to induce a flow in the liquid phase.²³ For non-contact sources, this flow can be explained by the mechanical coupling between the gas phase flow and liquid surface, where the gas flow can be that from a plasma jet or from the electrohydrodynamic (EHD) forces in an SBD. This gas flow sets the liquid's free surface in motion, leading to a flow in the bulk liquid due to shear stresses, which was reported in numerical and experimental studies.^24,25 For contact sources, however, there are multiple processes at play, one of which is the mechanical coupling between the gas flow, generated by electrohydrodynamic (EHD) forces and the liquid phase.^26,27 The mechanical coupling typically leads to the formation of a dip at the interface as it pushes the liquid downwards.^12,24 Another process at play is the electric coupling, represented by electric surface stresses that can drive a flow in the bulk liquid phase, which typically leads to the formation of a cone or a hump at the interface by pulling the surface upwards.^20,28 Flows driven by electric surface stresses are the main concept behind electrosprays and EHD atomizers,^29,30 where a strong electric field is applied to a liquid, causing its surface to form a Taylor cone, which breaks into a jet of droplets. Considering its importance for many applications, this process has been studied extensively experimentally and numerically.^31–34 In the context of plasma–liquid interaction, previous works have focused on the pre-breakdown stage of the discharge and the initial deformation of the liquid before the discharge ignites,^20,28,35 where it has been reported that the surface of the water forms a hump for values of the applied voltage below the critical stability limit, for higher voltages, a Taylor cone was observed at the interface.²⁵ 

Few quantitative studies of the complex interplay between the mechanical and the electrical processes have been performed, despite its importance for the transport of reactive species across the liquid interface. For example, Mitsugi et al. have studied the plasma jet-induced flow in a water sample using particle image velocimetry (PIV),³⁶ where they reported that the velocity of the bulk liquid is of the order of 10 mm s⁻¹ and that surface pressure waves were driving the surface flow in the radial direction. Lai et al. used PIV and shadowgraphy imaging to study flow instabilities at a plasma–liquid interface,³⁷ where their analysis indicated the presence of a Kelvin–Helmholtz perturbation that causes sharp velocity shear, inducing vortices in the shearing layers.

The objective of this work was to quantify the influence of these processes in driving the liquid flow in plasma activated water (PAW) reactors, where the plasma contacts the liquid's interface. This was achieved by implementing a numerical model to describe a pin-water discharge configuration operating in air. The model was experimentally validated and was able to capture the EHD induced flow in the gas phase and the electric surface stresses exerted by the plasma's sheath on the interface. The remainder of this manuscript is organized as follows: Sec. II describes the experimental setup used to validate the numerical model, which is described in Sec. III. Section IV presents the experimental and numerical results of the study, and a conclusion is provided in Sec. V.

The experimental setup employed in this study was used to generate plasma in direct contact with a liquid surface and facilitate the measurement of the velocity fields in both the gas and liquid phases using particle image velocimetry (PIV).

A pin-water discharge configuration operating in air was setup to expose 125 ml volume of water. A 2 mm diameter tungsten pin mas mounted vertically 5 mm above the water surface. The pin was connected to a custom-made high voltage power supply, capable of providing a sinusoidal signal in the kV range at 10 s of a kHz frequency range. The ground electrode was a copper wire inserted into the liquid volume at the container wall, approximately 35 mm from the discharge contact point. The discharge voltage and current were measured using a Tektronix DPO 5054 oscilloscope and a Pearson 2877 current probe. The dissipated power was calculated by averaging the product of the voltage and current over 450 cycles in real-time on the oscilloscope.

Two samples of water were considered in this study, a de-ionized sample with a conductivity of 2 μS cm⁻¹, measured using a Hanna Instruments HI 1285-5 multiparameter probe and a tap sample with a conductivity of 300 μS cm⁻¹. For each test, 125 ml volume of liquid was pipetted into a plastic container measuring 50 × 50 × 75 mm in width × length × depth. All samples were plasma treated at a constant power of approximately 15 W. To achieve a dissipated power of 15 W, peak to peak voltages of 11 and 8 kV were required for the de-ionized and tap water samples, respectively. During operation, it was apparent that plasma formation caused distortion of the applied voltage waveform due to the loading of the high voltage transformer. To ensure comparable behavior between the model and the experiment, the recorded voltage waveforms were used as inputs to the model. A schematic diagram of the discharge configuration is shown in Fig. 1(a).

While the discharge was in operation, PIV measurements were performed to obtain the velocity vector fields in both the gas and the liquid phase. The full operating principle of the PIV technique is described elsewhere,³⁸ here it is explained in brief. The entire experimental setup was placed within a large sealed chamber, measuring 1 × 1 × 1.5 m in height × width × length. The chamber had two windows, one enabled a laser light sheet to be projected across the experimental setup and the other enabled an observation to be made perpendicular to the laser sheet, as shown in Fig. 1(b). The chamber was filled with oil droplets, with a nominal diameter of 1 μm, to scatter the laser light. To obtain velocity measurements in the liquid volume, melamine particles with a nominal diameter of 10 μm were used. When the discharge was in operation, a 30 mJ pulsed Nd:YLF laser operating with a pulse duration of 100 ns at 527 nm was used to produce two consecutive laser sheets with a controllable time delay. Simultaneously, a high-speed camera (Phantom Micro 340) was used to capture a single laser pulse in each frame. Using a cross-correlation procedure, the 2D velocity vector field was calculated for each pair of camera frames.³⁹ Considering that the Stokes number of the seeding particles for both the gas and liquid phase measurements was less than 0.1, the velocity of the particle seeds was within 1% of the background flow.^40,41

Considering the significant difference between the liquid velocity and gas velocity, the time delay between the laser sheets was set to 50 μs for the gas phase and set to 5200 μs for the liquid phase at laser pulse repletion frequency of 400 and 185 Hz, respectively. Consequently, it was not possible to obtain gas and liquid velocity measurements simultaneously. For each experimental run, a total of 500 frames were captured and averaged to calculate the time-averaged velocity field after a given time, corresponding to 1.25 s for the gas phase and approximately 2.7 s for the liquid phase. Due to the scattering of laser light from the liquid surface, the PIV processing interrogation area was set 0.5 mm below the interface when calculating the liquid velocity field. For the same reason, in the gas phase, the interrogation was set 0.5 mm above the interface, as labeled in Fig. 1(a). It was not possible to perform PIV measurements for the de-ionized liquid sample due to severe surface perturbation resulting from plasma interaction. A video showing the unstable surface for the de-ionized sample case can be found in the supplementary material.

A two-dimensional time dependent axisymmetric computational model was used in this study, consisting of two components: a mechanical component for capturing the fluid flow in both the gas and liquid phases and an electric component for capturing the plasma dynamics. The model was designed to represent the experimental conditions as close as possible; therefore, all the parameters of the model such as the geometry, and the liquid sample's conductivity were chosen to match with those used in the experiments. The inputs to the model were the voltage waveforms from the experiments. The outputs of the model were the velocity fields in both phases, in addition to the electric field and the species densities in the plasma. The computed velocity fields were compared to those measured using PIV in order to validate the computational model. Despite the difference in the geometry between the experimental setup and the axisymmetric nature of the model, such a difference will not have a significant influence on the comparison given that the plasma–liquid contact area is much smaller than the surface area of the gas-liquid interface and is situated in the center of the container well away from the walls. The domain of the computational model and the relation between its components is shown in Fig. 2.

Initially, the mechanical component of the model was run to compute the induced flow by the electric field just before breakdown, in addition to the initial surface structure and induced surface charge on the surface. Those parameters served as the initial conditions for the electric component, which was solved for the time-averaged EHD force fields and electric surface stresses at the interface, which were subsequently used as inputs to the mechanical component, in order to solve for the plasma-induced flow.

In formulating the model, it was assumed that the plasma treatment time of the water samples is not long enough to alter their conductivity. This assumption was justified as the model followed the first few seconds of the plasma treatment. Experimentally, it was shown that a significant change in the conductivity of the treated water occurs on the order of minutes of treatment.⁴² 

The mechanical component of the model solved for the flow field in both phases, the plasma-free electric field, the interfacial surface charge from the liquid phase, and the free surface (the interface) deformation. For the fluid flow, the mass continuity equation and Navier–Stokes equation were solved in both phases, as given by Eqs. (1) and (2), respectively,

In Eqs. (1) and (2), ρ is the fluid’s mass density (kg m⁻³), which was calculated by the model in air assuming the ideal gas law and was set to 998 in water meaning that water was treated as an incompressible fluid, $u⇀$ is the velocity field in both phases (m s⁻¹), p is the pressure (Pa) that was obtained from the ideal gas law in the gas phase, and was computed as a variable in the liquid phase, I is the identity matrix, $τ¯¯$ is the stress tensor given in Eq. (3), $g⇀$ is gravity's acceleration (m s⁻²). The gravity acceleration term was only applied to the liquid phase. Finally, F_EHD is the time-averaged force density due to the EHD forces exerted by the charged species in plasma in the gas phase (N m⁻³). This term was set to 0 when the mechanical component was initially solved. To close the system, an explicit form of the stress tensor is needed. The stress tensor is given by Eq. (3),

In Eq. (3), ζ is the dynamic viscosity of the fluid (Pa s), which was set to 1.81 × 10⁻⁵ and 1 × 10⁻³ for air and water, respectively, ε is the electric permittivity of the medium that was set to ε₀ (F m⁻¹) for air and 81ε₀ for water, $E⇀$ is the electric field (V m⁻¹) that was determined from Eq. (4).

Equation (3) has two types of stresses, the first is the viscous stress represented by the first two terms, the second is the electric field stress represented by Maxwell's stress tensor appearing as the last two terms in Eq. (3).⁴³ Considering that F_EHD accounted for the electric forces in the gas phase, the Maxwell stress tensor was accounted for in the liquid phase only. To calculate the electric field, the Laplace equation was solved for the electric potential everywhere in the computational domain, as given in Eq. (4),

In Eq. (4), V is the electric potential (V), and the electric field was computed from the potential by $E⇀=−∇V$⁠. It should be noted that the electric field was also solved in the mechanical component to account for any field-induced deformation of the interface prior to breakdown. When an electric field is applied in the liquid phase, where positive and negative ions exist, counterions are attracted from the bulk to the interface to shield the electric field, forming surface charge. To compute that surface charge, Eq. (5) was solved. The surface charge density on the interface was treated as a boundary condition when solving Eq. (4),

In Eq. (5), σ is the surface charge density (C m⁻²), $n^$ is the normal unit vector on the surface, μ is the ionic mobility in water (m² V⁻¹ s⁻¹), n_i is the initial ion density in the liquid (m⁻³). Considering that the ions mobility and density appear as a product in the equation, the value of their product μn_i was set to match the experimental conductivity for the investigated water samples.

The imposed boundary conditions on the interface enforced the continuity of the flow field, as given in Eqs. (6a)–(6c),

In Eqs. (6a)–(6c), the “[ ]” brackets represent the “jump” across the surface from the liquid phase into the gas phase. Equation (6a) is a no-permeability condition, Eq. (6b) is a no-slip condition (meaning the velocity components are continuous along the interface), and Eq. (6c) states that the jump in the normal stress is equal to surface forces. The first two terms in the right-hand side of Eq. (6c) represent surface tension forces, where γ is the surface tension coefficient (N m⁻¹), which was set equal to 0.073 for the air–water interface. Revisiting the electric field stress given in Eq. (3); considering that the width of a typical electric double layer at the air–water interface is in the nm range,²³ i.e., much smaller that the dimensions of the model, it was possible to treat the electric stresses as a surface force instead of a volumetric force as described by Eq. (3). Consequently, the last two terms in the right-hand side of Eq. (6c) represent the electric “surface” stress, and the “volumetric” electric stress is no longer accounted for in Eq. (3).

Since Eqs. (6a)–(6c) allow for a non-zero velocity to exist on the interface, this means that the interface may deform due to its flow velocity. To capture the deformation of the interface, the arbitrary Lagrangian–Eulerian (ALE) method was implemented.⁴⁴ This method splits the coordinates used in the model into spatial coordinates (r,z) and reference coordinates (R,Z), as the location of the interface evolved, a smoothing equation was solved to map the reference coordinates to the spatial coordinates, which is given by Eq. (7),

In Eq. (7), r and z (mm) are the spatial coordinates (the instantaneous coordinates as the simulation was solved), while R and Z (mm) are the reference coordinates (the initial coordinates at the start time of the model). Equation (7) was connected to Eqs. (6a)–(6c) through Eqs. (8a) and (8b), where $usurf$ and $wsurf$ are the velocity components in the r and z directions, respectively (m s⁻¹),

For the remaining boundaries, a symmetry boundary condition for all variables solved for was imposed on the boundary r = 0. For the boundaries at z = 30 mm and r = 25 mm, the normal viscous stresses and electric fields were set to 0, which represented an open boundary. A no-slip boundary condition was imposed on the pin's surface for the fluid flow, and a time varying potential was imposed for solving Eq. (4). For boundaries making container of the water sample, a no-slip boundary condition was imposed, and the ground electrode was imposed at the bottom boundary.

The electric component of the model accounted for the generation of the plasma and how it affected the electric field. The electric component was based on a plasma fluid model consisting of a system of mass continuity equations following the evolution of 6 species in space and time and 13 reactions associated with them. The species were e⁻, N₂⁺, O₂⁺, O₂⁻, O₂, and N₂ and the reactions included in the model and their rate coefficients are given in Table I. For every species, Eqs. (9a) and (9b) were solved. In addition to one additional conservation equation for the electron energy density, as given by Eq. (10),

In Eqs. (9) and (10), n_k is the density of the kth species (m⁻³), $Γ⇀k$ is the flux of the kth species (m⁻² s⁻¹), k_lj is the two-body reaction rate coefficient (m³ s⁻¹), μ_k is the mobility of the kth species (m² V⁻¹ s⁻¹), D_k is the diffusion coefficient of the kth species (m² s⁻¹), n_en is the electron energy density (eV m⁻³), μ_en and D_en are the mobility and the diffusion coefficients, respectively, for the electron energy density, and $θl$ is the electron's energy cost of the lth reaction (eV).

The transport coefficient of ions and neutrals, namely, the mobility and the diffusion coefficient, were set equal to those in other work.⁴⁵ The mobilities of the neutrals were set to 0. For electrons, the transport coefficients, μ_en and D_en were all calculated by BOLSIG+.⁴⁷ 

Since the plasma alters the electric field, Eqs. (4) and (5) were updated with a plasma term to become Eqs. (11) and (12),

The term on the right-hand side of Eq. (11) is the space charge density of the plasma (C m⁻³), while the second term on the right-hand side of Eq. (12) is the net charge flux to the interface from the plasma.

While the electric component is solved, the instantaneous EHD forces and the instantaneous electric surface stresses were integrated in time.²⁷ After the electric component is finished, the time-integrated quantities were divided by the period of the waveform, thus giving the time-averaged EHD forces and the time-averaged electric surface stresses, as given in Eqs. (13) and (14). Both terms are subsequently used as an input to the mechanical component to account for the plasma-induced flow,

In Eqs. (13) and (14), T is the period of the waveform (s) and $F⇀E$ is the time-average electric surface stress pressure (N  m⁻²). After these variables were evaluated, they were introduced in Eqs. (2) and (6c) in the mechanical component of the model. This approach is applicable as long as the deformation of the surface during a single period is negligible.

All the plasma parameters were solved only in the gas phase. A symmetry boundary condition was imposed at r = 0 boundary for all variables solved for. For the outer boundaries of the gas phase, a continuity condition was imposed. For the boundary defining the pin electrode, a wall loss condition was imposed, which means all the diffusive flux to the wall is lost, in addition to electric-field driven flux when it is directed toward the electrode.

Initially, the mechanical component of the model was run until the surface arrived its steady state structure. Then, the electric component was run for one period of the applied waveform, and the mechanical component was subsequently run (taking the plasma effects into account) until reaching a steady state flow pattern, which took few seconds of solution time. The time-averaged power dissipated at the pin electrode as calculated by the model was 10 W for the de-ionized sample case and 9 W for the tap sample case.

At the instant prior to breakdown, the applied electric field was strong enough to cause a small “hump” to form on the interface. For the applied voltages and the pin-interface distance used in experiments, the system was under the critical stability limit for behaving as an electric spray, thus no jetting was observed experimentally, which was consistent with the model.

To form the hump on the surface, the flow in the liquid phase was in a direction that pumped water against gravity. This can be seen from Fig. 3, which shows the liquid’s flow direction in Fig. 3(a), and the axial velocity component in the liquid phase in Fig. 3(b).

Figure 3(a) shows that the direction of the water flow was from the bulk to the interface at the outer boundary of the container, then parallel to the surface to form the hump at the symmetry axis and then back into the bulk. The formation of the hump caused the hydrostatic pressure to increase locally at the symmetry axis, which created a pressure gradient that drove the liquid flow downwards at the symmetry axis, as indicated in Fig. 3(a). The flow of the liquid at the interface dragged the gas, causing it to be sucked toward the symmetry axis, where the build-up of the pressure created a stream in the positive z direction as can be seen close to the pin electrode.

To quantify the flow field, Figs. 3(a) and 3(b) show the axial velocity component in the liquid phase. The velocity at the surface is zero as the stagnation pressure is developed in the gas side of the interface. Moving down in the liquid, the velocity was directed toward the bottom of the container and peaked at 8 and 2.1 mm s⁻¹ for the de-ionized and tap water samples, respectively. The difference in the peak velocities stemmed from using different amplitudes to drive the discharge at the same operating power.

The behavior of the discharge in a pin-water configuration has been widely studied and the findings in this study are in line with the past findings. The discharge was observed to ignite in the positive cycle of the applied voltage, originating at the pin electrode, where the electric field had its maximum value due to the sharp geometry at the pin's tip. A plasma cloud formed and subsequently propagated in the gap between the pin electrode and the interface as a streamer, forming a conductive plasma channel between the pin electrode and the liquid. The two-dimensional profile of the plasma channel as it formed and propagated is shown in Figs. 4(a) and 4(b). The model showed that the estimated propagation velocity of the streamer was approximately 15 km s⁻¹, and the radius of the streamer was slightly above 1 mm. Both parameters are consistent with experimental parameters reported for similar configurations.^48–50 

Figure 4(b) also shows that as the plasma channel approaches the interface, the discharge mode transitions from a streamer-like mode into a glow-like mode, where a plasma is ignited between two conductors, leading the electron density to increase throughout the plasma channel reaching electron density of 10¹⁹–10²⁰ m⁻³, which is consistent with experimental reports on similar discharges.⁴⁸ 

As soon as the plasma channel was established, its interaction with the interface varied between the two investigated cases. For the de-ionized sample, the plasma expanded along the interface as a surface streamer, creating a branched plasma channel parallel to the surface and thus increasing the contact area. For the more conductive tap water sample, however, this effect was less profound. For both investigated cases, the area of the plasma–water interface was less than 1% of the total surface area of the air–water interface, thus supporting the assumption that the influence of the difference in the geometry between the model and the experimental setup on the results was minimal.

The difference in the behavior at the interface can be explained in terms of the conductivity of the liquid samples. As stated in Sec. II, the conductivity of the de-ionized sample was 2 μS cm⁻¹ while that of the tap sample was 300 μS cm⁻¹. Critically, the conductivity of the liquid controls affects its charge relaxation time, which defines the characteristic time required for the liquid to shield an external electric field. For an ideal dielectric, the relaxation time is infinity as dielectrics do not shield the external electric field, instead they polarize. The relaxation time of a perfect conductor is zero. The charge relaxation times were calculated to be approximately 3.6 μs for the de-ionized sample and 24 ns for the tap sample.

The plasma's propagation along the water's surface occurred on a tens of nanoseconds timescale. Consequently, the surface discharge behaved liked a dielectric barrier discharge (DBD) in the de-ionized sample case, where the plasma deposits surface charge on the interface, leading to a weakened electric field and causing a further radial spread along the surface. In contrast, when a tap water sample was used, the plasma's propagation time was comparable to that of the charge relaxation time of the sample, thus the deposited surface charge started to decay as soon as it was deposited on the interface, making the discharge behave as a hybrid between a DBD and a glow discharge.

Figure 4(c) shows the surface charge density for both cases 0.5 μs after the streamer arrived at the interface. As can be seen, the deposited surface charge for the de-ionized sample is an order of magnitude higher than that of the tap sample. Another observation is that the surface charge distribution in the tap sample case is narrower than that in the de-ionized case. This is because of surface charge neutralization, which weakened the radial propagation of the discharge.

The time-averaged EHD force field for both investigated cases showed a similar structure with different amplitudes. The radial and the axial force fields are shown in Fig. 5, where it can be seen that the maximum values of the force fields coincided with the positions of the plasma sheath, namely, the tip of the pin electrode, the contact area between the plasma and the liquid, and in some cases, the radial sheath of the plasma channel. This affiliation can be explained by Eq. (13), as the plasma sheaths are the only regions where a non-zero charge density exists and as a result a strong electric field exists. Therefore, they are the only regions where the EHD forces are exerted.

Focusing on the radial force, Figs. 5(a) and 5(b) show that the maximum values of the force, directed radially inward and outward, are at the tip of the pin electrode. This is a consequence of the difference in the widths of the sheaths during the positive and the negative cycles of the waveform. During the positive cycle, the plasma channel was at a lower potential than the electrode; therefore, the electric field was directed radially outward, and the direction of the radial force was radially outward as a result. During the negative cycle, the plasma channel had a higher potential than the pin electrode and the direction of the radial force was radially inward. Considering that streamer propagation occurred in the positive cycle, the negative cycle started at a higher plasma density than the positive cycle, which led to a shorter Debye length and a narrower sheath, thus explaining why the time-averaged radial force is negative on the surface of the pin electrode.

The other region of interest in the computational domain is the contact area between the plasma and the liquid, where the time-averaged radial force has a local maximum in the radially outward direction for both samples. This force was exerted during the positive cycle where the plasma potential was higher than its surrounding, which meant that the electric field was in the radially outward direction. During the negative cycle, the electric field changed its direction, however, because the plasma density was low at the contact area, the difference between the plasma potential and its surroundings was low, which led to a weak electric field. Consequently, the positive cycle behavior dominated the time-averaged EHD force in the contact area at the plasma–liquid interface.

The main difference between Figs. 5(a) and 5(b) is the presence of a stronger radial force at the side sheath of the plasma channel in Fig. 5(a) compared to Fig. 5(b). As explained in Sec. IV B, the liquid in the de-ionized sample [Fig. 5(a)] behaved as a lossy dielectric. This meant that the surface charge accumulation on the interface weakened the electric field in the sheath at the contact area, thus the plasma potential was primarily influenced by the pin electrode. Indeed, the variation of the plasma potential across the plasma channel did not deviate significantly from the pin electrode's potential at any given time. For the tap sample [Fig. 5(b)], the liquid behaved as a conductor with a finite conductivity, which meant that the ground electrode submerged in the water influenced the plasma potential, such that the plasma potential in the tap sample case was lower than that in the de-ionized sample case. This led to a weaker radial electric field and a lower time-averaged radial force at the side sheath of the plasma channel.

The time-averaged axial force for the de-ionized and tap water samples are shown in Figs. 5(c) and 5(d), respectively. Similar to the radial force, the maximum values of the axial force, whether it is upwards or downwards, exist at the tip of the pin electrode. The positive cycle contributed to the downward force, while the negative cycle contributed to the upward force. For the same reasons explained earlier, the contribution of the negative cycle extended over a smaller region than that of the positive cycle, leading to the force structure observed in Figs. 5(c) and 5(d).

The sheath at the plasma–liquid contact area shows a significant difference between Figs. 5(c) and 5(d), where there is a strong force component in the tap sample case but not in the de-ionized sample case. The absence of the strong force in the de-ionized sample case can be attributed to the accumulation of surface charge on the interface, which attenuated the axial electric field in the sheath leading to a weaker axial force. The high conductivity of the tap sample meant that significant charge accumulation could not occur on the surface, which allowed for a strong electric field to exist over longer periods of time, which contributed to the time-averaged EHD force.

Despite it not being obvious in Fig. 5, it should be noted that non-zero time-averaged force components exist in the plasma channel. The origin of this time-averaged force was the instantaneous force exerted by the streamer head as it propagated in the gap between the pin electrode and the liquid surface. As soon as a plasma channel was established in the gap, the instantaneous force dropped to zero as the plasma channel was charge neutral. The time-averaged force components in the plasma channel were fairly uniform, where the axial component had approximate values of 400 and 1800 N m⁻³ for the tap sample and the de-ionized sample, respectively, in the downward direction. For the radial component, the approximate values were 20 and 80 N m⁻³ at a radius of 0.2 mm in the outward direction for the tap and de-ionized samples, respectively. The higher values obtained for the de-ionized sample can be attributed to the higher voltage needed to drive the discharge for the de-ionized case.

While the electric component of the model was integrated in time, the electric surface stresses across the interface were integrated in time. The time-averaged electric surface stresses across the interface are shown in Fig. 6, which indicates that the electric surface stresses in the de-ionized sample case were an order of magnitude higher than those for the tap sample case. It also shows that the direction of the electric surface stresses for the de-ionized sample was opposite to that of the tap sample. Finally, the electric surface stresses in the de-ionized sample case extended over a larger area than that of the tap sample case.

The tangential electric surface stress of the de-ionized sample is much larger than that of the tap sample case, as Fig. 6(a) shows. This stress is caused by the tangential electric field along the interface. As is known from electromagnetism, the tangential electric field of a perfect conductor is zero.⁵¹ Considering that the tap sample has a significantly higher conductivity than the de-ionized sample, the resultant tangential electric field was much lower, in the order of 1 kV cm⁻¹ at its maximum. In contrast, the tangential electric field in the de-ionized sample was on the order of 10 kV cm⁻¹ at its maximum, which led to larger tangential electric surface stress at the interface.

For the normal electric surface stress, the explanation of the differences between the two cases comes from the deposited surface charge on the interface. As explained in Sec. IV B, the de-ionized sample maintained the surface charge for a longer time due to its lower conductivity. The accumulated surface charge on the interface weakened the electric field in the plasma side of the interface. On the liquid side, however, the electric field was enhanced due to the presence of the surface charge, as highlighted in Fig. 6(c). Considering that the electric stress is a function of the square of electric field, the intensified electric field in the liquid medium on the order of 20 kV cm⁻¹ at its maximum caused the electric stress in the liquid to be comparable to that in the plasma, which was on the order of 60 kV cm⁻¹ at its maximum. Since the dielectric constant of water is 81 times that of air, the normal electric stress shown in Fig. 6(b) had a negative jump as when moving from the liquid phase to the plasma phase. For the tap water, the accumulation of the surface charge was small thus the electric field in the plasma side of the interface was stronger than that in the liquid side, as shown in Fig. 6(d), which was less than 0.05 kV cm⁻¹ at its maximum. This led to a positive jump when moving from the liquid phase into the plasma phase. These observations explain the difference in the polarity between the de-ionized and tap water samples, it also explains the difference in the amplitude of electric surface stresses between the de-ionized sample and the tap sample. Finally, considering that the water behaved as a dielectric in the de-ionized sample case, the plasma extended over a larger area in comparison to the tap sample, which led to the electric field stress for the de-ionized case to extend over a larger area than that of the tap sample case.

After the EHD force fields and electric surface stresses on the interface were computed by the electric component of the model, they were used as an input to the mechanical component to account for the plasma-induced flow, which was solved until the flow arrived a steady state. This allowed for the comparison of the model predicted flow field to that of the time-averaged PIV measurements. Figure 7 shows a comparison between the computed and measured velocity magnitude in the air and water regions, respectively. Reasonable agreement between the computed and measured flow fields can be observed. As can be seen in Fig. 7(a), there is a gap in the measured data extending from the pin tip to the liquid surface. Unfortunately, this is an unavoidable consequence of the plasma emission masking scattered laser light from the seeding particles, making it impossible to calculate the velocity vectors within the most intense parts of the discharge. Figure 7 also shows that the air flow in the computed and the measured cases eventually go up, with the transition point occurring at a closer distance to the plasma in the measurement. This is a result of the higher temperature of the flowing gas in comparison to the ambient air, which sets buoyancy in action when the flow velocity is low, which was not taken into consideration in the model.

Figure 7 also indicates that the computed velocity magnitude in the liquid phase is slightly higher than the measured profile, despite having a similar structure. Considering that the model is deterministic and based on time-averaged force fields, it was unable to capture some of the flow instabilities observed experimentally on plasma–liquid interfaces.³⁷ The absence of such instabilities meant that momentum transfer across the surface was overestimated, which led to a higher velocity in the liquid phase.

In the gas phase, the EHD axial force induced a “jet” flow directed toward the liquid surface, with an average velocity in the region between the pin electrode and the interface of 1.5 ms⁻¹. The jet flow impinged on the surface and turned into a radial flow parallel to it. The radial EHD force in the contact area gives a further push to the radial flow. The radial gas flow dragged the interface through viscous stress coupling across the surface, which induced a shear flow in the bulk liquid on the order of 12 mm s⁻¹. The direction of the radial flow in the liquid phase was parallel to that in the gas phase. The induced axial flow in the liquid phase was directed in the positive z direction in the order of 15 mm s⁻¹, transporting liquid from the bulk to the surface to provide the mass flow to support the radial flow at the surface. Comparing Figs. 7 to 3, it becomes evident that the shear flow in the liquid phase is dominated by the gas flow after breakdown, in comparison to electric surface stresses before breakdown. The contribution of the electric surface stresses after breakdown on the surface to the overall flow was found to be insignificant in the tap sample case. This was inferred by comparing the solution of the model when both mechanisms (EHD flow and electric surface stresses) are accounted for, to the solution when the EHD flow only is taken into account. The liquid flow was found to be less than 10% different between the two cases.

With respect to the de-ionized sample case, despite the absence of PIV data, the computed EHD forces and the electric surface stresses were input to the mechanical component to quantify their contribution to the predicted flow. The model showed that when both mechanisms were considered, the surface experienced large deformations and became unstable, which was consistent with the experimental observations. When the model was solved with EHD flow only, the liquid phase was stable and the deformation of the surface was minimal. To quantify this effect, the model was solved using the EHD force only until steady state flow in the liquid phase was established. At time 0, the electric surface stress was “switched on.” Figure 8 shows the deformation of the surface and the axial velocity on the symmetry axis over the first millisecond of switching on the electric surface stress. In this scenario, data were only captured over a short operating period as the large deformation of the interface required the electric component of the model to be updated to account for the new interface structure.

Figure 8 indicates that the electric surface stresses dominated the liquid flow in the de-ionized water case, and that it is driving the instabilities at the interface in the first minutes of plasma generation. The surface instabilities cause the flow to become highly turbulent at the interface. Considering that the Sherwood number, which is a dimensionless parameter describing the rate of mass transfer across an interface, increases as the Reynolds number increases; a turbulent interface leads to a larger mass transfer rate across the interface, which would result in a faster activation of the treated liquid.⁵² 

It was observed experimentally that longer plasma exposures led to an increasingly stable flow in the liquid phase. This may be explained by considering the conductivity of the liquid, which is known to gradually increase during plasma exposure due to the increase of NO₂⁻ and NO₃⁻ in the liquid.⁴² This acts to transform the liquid from a leaky dielectric into a poor conductor, thus the de-ionized water sample becomes more like the tap water case.

A two-dimensional axisymmetric model was developed to describe the mechanical interaction between the plasma and the treated liquid in a pin-water electrode system. The model solved for the plasma parameters, the flow fields in the gas and the liquid phase, and the deformation of the surface. Particle image velocimetry measurements were used to capture the velocity flow field in both the gas and the liquid, providing data to validate the computational model. Close agreement was observed between the model's predictions and the experimental measurements. The model was solved for two water different samples: a de-ionized sample with a low conductivity of 2 μS cm⁻¹ and a tap sample, which had a conductivity of 300 μS cm⁻¹.

The model was used to investigate two mechanisms that contribute to driving the liquid flow, those were the EHD-induced flow in the background gas and the electric surface stresses across the interface. The time-averaged EHD force field was analyzed, and it was reported that the primary difference between the two investigated cases was in the axial component at the plasma–liquid contact area, where the force for the tap sample case was significantly larger than that of the de-ionized case. This was attributed to the strong electric field being maintained in the sheath for most of the applied period of the waveform as a result of the short lifetime of the surface charge on the interface. The electric surface stresses were shown to be an order of magnitude larger for the de-ionized sample in comparison to that of the tap sample. This was explained by the presence of the longer lifetime of the surface charge on the interface, which weakened the electric field in the air side of the plasma side of the interface and strengthened the electric field in the liquid side, leading to larger electric surface stresses.

The influence of the two mechanisms on the liquid flow was investigated. It was shown that the electric surface stresses played a minimal role in the tap sample case, while in the de-ionized sample case, they were responsible for large deformation of the interface, which significantly affected the liquid flow. These findings indicate that the electric stresses are the dominant mechanism driving the liquid flow when the liquid behaves as a dielectric, while the EHD-induced gas flow is the dominant mechanism driving the liquid flow when the liquid behaves as a conductor.

The large surface deformations reported for the de-ionized sample case led to strong turbulence at the interface. Considering that the mass transport rate across an interface is proportional to the Reynolds number, the findings of this study imply that the plasma activation of a liquid with a low conductivity occurs as a higher rate that those with a lower conductivity, leading to faster activation of the treated liquid.

See the supplementary material for a video recording of the water's interface to show instability at the interface arising when plasma treating the de-ionized water sample as described in the manuscript.

J.L.W. and M.I.H. would like to acknowledge support of Engineering and Physical Sciences Research Council (ESPRC) Grant Nos. EP/S025790/1, EP/N021347/1, and EP/T000104/1.

The data that support the findings of this study are available from the corresponding author upon reasonable request.