Coupled electromechanical-hydrodynamic modeling of a low-frequency vibratory cavitator

1. Introduction

Hydrodynamic cavitation is increasingly applied to water and wastewater treatment because cavity collapse intensifies mixing and physicochemical degradation. Reviews [1]-[3] summarize advanced-oxidation mechanisms, reactor configurations, pollutant removal, and energy-efficiency constraints. For perforated restrictions, Maynes et al. [4] and Malavasi et al. [5] showed that head loss and cavitation inception depend on the free-area ratio, hole geometry, plate thickness, and Reynolds number. Ezoji et al. [6] demonstrated that oscillating perforated plates exhibit amplitude-dependent added mass and nonlinear damping.

Electromagnetic oscillators must be modeled jointly with their mechanical load. Hung and Lim [7] experimentally confirmed the influence of voltage and air gap on current, force, and plunger displacement, whereas Naseradinmousavi and Nataraj [8] formulated nonlinear actuator dynamics with position-dependent inductance. The Rayleigh-Plesset equation [9] remains the classical basis for reduced bubble-growth and collapse descriptions. Aftanaziv et al. [10] developed an electromagnetic vibratory cavitator for low-frequency liquid treatment.

The authors’ previous study [11] proposed a continuous-flow device with two oppositely moving perforated grates, but prescribed the electromagnetic forces harmonically and represented liquid loading by generalized viscous damping. A recent hydrodynamic-cavitation generator driven by a linear electromagnetic piston [12] further demonstrates the relevance of coupled actuator-fluid analysis. However, a compact model combining a variable-gap electromagnetic drive, oscillating perforated-grate hydrodynamics, and finite-rate vapor-content evolution has not yet been developed for the considered continuous-flow cavitator. The present study therefore evaluates the effects of excitation frequency on grate motion, local pressure, vapor content, current, power demand, and hydraulic dissipation.

Experimental studies involving Aftanaziv and co-workers have investigated vibrocavitation of brewing wastewater [13], frequency-dependent treatment of model contaminated water [14], and gas-assisted ultrasonic treatment of fat-and-oil-plant wastewater [15]. These studies establish relevant treatment experience, but do not validate the dynamic predictions of the present device. Accordingly, the discussion below distinguishes existing experimental evidence from planned model validation and examines the changes required for viscous and non-Newtonian liquids.

2. Research methodology

The computational diagram of the cavitator is shown in Fig. 1. The body of mass $m_1$ is connected to the stationary base by the stiffness $k_1$ and damping $c_1$, whereas two identical armature-rod-grate assemblies $m_2=m_3=m$ are connected to the body by elastic and dissipative elements $k$ and $c_s$. The electromagnetic drives act in opposite global directions. The absolute displacements $x_1$, $x_2$, and $x_3$ are shown in Fig. 1. The local displacements of the moving assemblies toward the corresponding stators are defined as $y_2=x_1-x_2$ and $y_3=x_3-x_1$. Under the synchronous symmetric operating conditions considered in Section 3 of this research, the local motions satisfy $y_2=y_3=y$, while the global grate motions are antiphase and their dynamic reactions on the cavitator body are mutually compensated. Therefore, both moving assemblies can be represented by one local coupled subsystem.

Fig. 1Computational diagram of the coupled cavitator model: m1 – cavitator body; m2, m3 – armature-rod-grate assemblies; x1, x2, x3 – absolute displacements; Fem,j and Fh,j – electromagnetic and hydrodynamic forces; Qh,j – liquid flow rates through the perforated grates

The reduced mechanical equation is written as:

where the relative flow rate through the perforated grate, nonlinear pressure-loss increment, added liquid mass, and mixture properties are:

Eq. (2) accounts for hydraulic inertance, linear and quadratic losses, and the dependence of the liquid-vapor mixture properties on the instantaneous vapor fraction. The electromagnetic subsystem is voltage-driven and includes the variable air gap, magnetic saturation, and motion-induced electromotive force:

where $L_d=\partial \lambda /\partial i$, ${\lambda }_g=\partial \lambda /\partial g$, and the static magnetic force $F_{em,0}=F_{em}\left(\frac{U_b}{R_e},g_{eq}\right)$ is balanced by the initial spring preload.

For long-term frequency-response simulations, the cavitation dynamics are represented by a regularized pressure-based vapor-fraction model:

where, ${\alpha }_{eq}$ is the pressure-dependent equilibrium vapor fraction, $\mathrm{\Delta }{p}_s$ controls the smoothness of the cavitation transition, and ${\tau }_g$ and ${\tau }_c$ are effective vapor-cloud growth and collapse times, respectively. Eq. (4) is a phenomenological reduced-order closure rather than a direct integration of the Rayleigh-Plesset equation [9].

All pressures in Eq. (4) are absolute. Because the regularized closure introduces a smooth transition, the condition $p_{\mathrm{m}\mathrm{i}\mathrm{n}}\le p_v\left(T\right)$ is used only as a pressure-based indicator that the cavitation threshold has been crossed, rather than as a discontinuous phase-change criterion.

The adopted magnetic, hydraulic, and vapor-fraction coefficients require experimental identification for the present device. The calculations below therefore constitute model predictions; numerical solver agreement does not establish experimental validity.

3. Results and discussion

3.1. Numerical response for water purification process

Numerical simulation was performed for the coupled electromechanical-hydrodynamic model derived in Section 2. The calculations were implemented in Python 3.13.5 using NumPy 2.3.5 and SciPy 1.17.0. The reported results were obtained using the implicit Radau solver, while selected operating points were recalculated with the BDF method to verify numerical consistency. The maximum integration step was restricted to approximately one hundredth of the excitation period. Each frequency case was simulated for at least 35 oscillation cycles or 0.70 s, whichever was longer, and only the final eight cycles were used to calculate steady-state amplitudes, root-mean-square values, phase shifts, and mean powers.

The processing temperature and pressure selected in Table 1 correspond to the operating range proposed for the cavitator, namely approximately 75-80°C and 0.35-0.40 MPa. In the present simulations, the chamber pressure $p_c=$ 0.380 MPa was interpreted as an absolute pressure. The quantities $A_p$, $A_h$, $i_b$, and $F_{em,0}$ were calculated from the corresponding geometric and electrical parameters. Frequencies below 41 Hz were excluded because the predicted minimum air gap approached the imposed limit $g_{\mathrm{m}\mathrm{i}\mathrm{n}}=$ 0.8 mm, indicating the onset of electromagnetic pull-in or mechanical contact beyond the validity of the no-contact model.

Table 1Parameters adopted for the numerical simulation

Symbol	Value	Unit	Symbol	Value	Unit
$m_1$	29.5	kg	${\mu }_l$	3.55×10-4	Pa·s
$m_2=m_3=m$	0.85	kg	$A_p$	7.854×10-3	M2
$k_1$	1.9×104	N/m	$A_h$	4.241×10-4	m2
$c_1$	50	N·s/m	$\epsilon$	0.054	–
$k_2=k_3=k$	8.8×104	N/m	$t_p/d_h$	1.5	–
$c_s$	55	N·s/m	$Q_0$	6.0×10-3	M3/s (21.6 m3/h)
$g_{eq}$, $g_{min}$	5.0; 0.8	mm	$\chi$	0.15	–
$g_{Fe}$	0.8	mm	$K_0$	1.8	–
$N$	700	–	$R_h$	1.0×106	Pa·s/m3
$A_g$	2.0×10-3	m²	$l_h$	0.012	m
$I_s$	4.5	A	$l_p$	0.010	m
$L_l$	0.012	H	$C_p$	2.8	–
$R_e$	35	Ω	$A_e$	1.5×10-3	M2
$U_b$, $U_a$	84; 140	V	$c_{\tau }$	20	N·s/m
$i_b$	2.40	A	${\alpha }_{max}$	0.06	–
$F_{em,0}$	98.84	N	$k_{\alpha }$	2.5	–
$T$	80	°C	$\Delta p_s$	4.0	kPa
$p_c$	0.380	MPa (abs.)	${\tau }_g; {\tau }_c$	0.5; 1.5	ms
${\rho }_l$, ${\rho }_v$	971.8; 0.293	kg/m³	$f$	41-100	Hz
$p_v$	47.4	kPa	$y\left(0\right)$	0	m
$D_p$, $t_p$	100; 4.5	mm	$\dot{y}\left(0\right)$	0	m/s
$n_h$	60	–	$i\left(0\right)$	2.4	A
$d_h$	3.0	mm	$\mathrm{\alpha }\left(0\right)$	0	–

Fig. 2(a) shows the dependence of the minimum local absolute pressure and maximum vapor volume fraction on the excitation frequency. At 41 Hz, the minimum pressure decreases to 37.84 kPa, which is 9.56 kPa below the saturated vapor pressure of water at 80 °C. The corresponding maximum vapor fraction is 5.43 %. At 50 Hz, the minimum pressure equals 41.20 kPa and the maximum vapor fraction reaches 4.79 %.

The pressure-threshold crossing occurs at approximately 60-61 Hz. At 59 Hz, the minimum pressure is 46.64 kPa, whereas at 61 Hz it rises to 47.93 kPa. Above this range, the instantaneous pressure no longer falls below $p_v$, although a small modeled vapor fraction remains because the regularized closure in Eq. (4) gives a nonzero equilibrium fraction above the saturation pressure; finite relaxation additionally produces phase lag. At frequencies above approximately 78-80 Hz, the maximum vapor fraction becomes lower than 0.1 %. At 100 Hz, the predicted value is only 8.9×10^-4%, and cavitation is practically absent.

These results indicate that the cavitation intensity is governed not only by the nominal liquid-flow rate but also by the grate velocity and acceleration. The lower-frequency operating regime generates larger volumetric displacement of liquid through the holes and a stronger inertial pressure reduction. Thus, the enhanced mathematical model additionally identifies an approximate pressure-based cavitation interval of 41-60 Hz for the adopted parameter set.

Fi2(b) presents the electromechanical response at 50 Hz. In Fig. 2(b), $A_y$, $I_a$, and $F_a$ denote the half peak-to-pg. eak amplitudes of the displacement, current, and electromagnetic force, respectively, while overbars indicate cycle-averaged values. Although the strongest pressure reduction is predicted near 41 Hz, the corresponding grate amplitude reaches approximately 3.04 mm and the minimum air gap decreases to about 1.13 mm. At 50 Hz, the minimum instantaneous gap increases to approximately 2.59 mm, providing a margin of approximately 1.79 mm above the imposed lower limit $g_{\mathrm{m}\mathrm{i}\mathrm{n}}=0.8$ mm. The steady-state displacement amplitude is 2.205 mm, the root-mean-square grate velocity is 0.491 m/s, and the maximum acceleration is 221.1 m/s², or approximately 22.54$g$.

The instantaneous coil current varies from 0.70 to 5.00 A, and its root-mean-square value is 2.81 A. The electromagnetic force changes from approximately 9.3 to 270.7 N. However, the fundamental harmonic of this force has an amplitude of only approximately 133.5 N. Therefore, the maximum instantaneous force must not be used as an equivalent amplitude of an ideal sinusoidal excitation.

Frequency-domain analysis of the steady-state signals gives the following phase shifts relative to the applied voltage: ${\phi }_i\approx$ –29.6°, ${\phi }_{F_{em}}\approx$ –59.7°, ${\phi }_y\approx$ –149.8°. The substantial current and force phase lags confirm that the electromagnetic drive cannot be represented adequately by a mechanical force that is assumed to be perfectly synchronized with the supply voltage.

The calculated mean electrical input power is 322.7 W per coil. The corresponding copper loss is 276.2 W, while approximately 46.5 W is converted into mechanical power. The corresponding model-based electromechanical conversion efficiency is approximately 14.4 % per drive, excluding losses in the power converter. For two identical electromagnetic drives, the total electrical input is therefore approximately 645 W. Of the mechanical power delivered to each moving assembly, approximately 33.3 W is dissipated through the modeled hydraulic resistance and 13.3 W through structural damping.

The hydrodynamic-force amplitude, evaluated as one-half of its peak-to-peak range, is approximately 99.4 N at 50 Hz. Its fundamental harmonic is shifted by approximately 174.5° relative to the grate velocity and therefore acts predominantly as a damping force. The corresponding equivalent linear hydrodynamic damping evaluated from the fundamental components is approximately 138 N·s/m. Together with the assumed structural damping of 55 N·s/m, this gives approximately 193 N·s/m.

Fig. 2Predicted responses of the coupled model

a) Frequency dependence of the minimum local pressure and maximum vapor volume fraction

b) Normalized steady-state displacement, coil current, and electromagnetic force at $f=$ 50 Hz

c) Local absolute pressure and vapor volume fraction during three steady-state oscillation cycles at 50 Hz

As shown in Fig. 2(c), at 50 Hz, the local absolute pressure varies from approximately 41.2 to 175.2 kPa. The pressure remains below the saturated vapor pressure for approximately 15.8 % of each oscillation period. Due to the finite vapor-cloud collapse time, the vapor fraction remains above 1 % for approximately 30.4 % of the period.

The vapor fraction reaches its maximum after the instant of minimum pressure rather than simultaneously with it. This phase delay reflects the finite rate of vapor-cloud development.

The mean velocity through the grate holes is approximately 14.15 m/s and varies from 12.28 to 16.16 m/s during the cycle. The corresponding hole Reynolds number, $R{e}_h={\rho }_l\left|u_h\right|d_h/{\mu }_l$, varies approximately within 1.01×10⁵$\le R{e}_h\le$ 1.33×10⁵, confirming strongly turbulent through-hole flow and supporting the use of a quadratic pressure-loss term.

4. Experimental evidence and planned validation

Shevchuk et al. [13] experimentally investigated brewing yeast-containing wastewater using an electromagnetic vibratory cavitator and found that combined vibrocavitation and nitrogen treatment improved biological and organic decontamination. Aftanaziv et al. [14] measured chemical oxygen demand during low-frequency vibrocavitation of model contaminated water. For the reported uric-acid-to-water ratio of 1:10 under nitrogen, organic-contaminant destruction after one hour reached 70 %, 84 %, and 75 % at 25, 37, and 75 Hz, respectively. These published results demonstrate that treatment performance depends on excitation frequency, but the 37 Hz optimum for that apparatus cannot be transferred to the present geometry and operating conditions.

Vashkurak et al. [15] measured chemical oxygen demand and microbial counts during ultrasonic treatment of fat-and-oil-plant wastewater in different gas environments, demonstrating the influence of gas composition on treatment efficiency. That study concerns aqueous wastewater and ultrasonic cavitation; it is not an experiment on pure oil or on the present low-frequency drive. Collectively, these studies support the technological feasibility of cavitation treatment, but their process-level measurements cannot establish the accuracy of the displacement, current, local pressure, or vapor-fraction predictions reported here.

The theoretical predictions developed in this article will undergo dedicated experimental validation in subsequent research using the laboratory equipment developed at Lviv Polytechnic National University (see Fig. 3). Synchronized measurements of grate displacement, coil voltage and current, chamber pressure, local pressure near the holes, flow rate, and temperature will be compared with predicted amplitudes, phases, pressure minima, and mean electrical power. Optical observations will assess cavity onset and cloud timing; quantitative vapor-fraction validation will require a calibrated measurement method. Parameters will be identified using a subset of operating points and tested at separate frequencies, with repeated measurements and uncertainty estimates. Water tests will precede investigations of viscous and non-Newtonian liquids.

Fig. 3Experimental prototype of the vibratory cavitator

5. Numerical response for water purification process

The present predictions apply to water at 80 °C. For illustration, the adopted mean flow rate and hole geometry give a mean hole velocity of 14.15 m/s. A hypothetical Newtonian oil with density 850 kg/m³ and dynamic viscosity 0.10 Pa·s at its operating temperature would have a mean hole Reynolds number of approximately 361 at the same imposed flow rate, compared with approximately 1.16×10⁵ for water. These are scaling estimates, not additional simulations. Low Reynolds number does not itself exclude cavity formation: experiments comparing water and oil in the same nozzle reported cavities in laminar oil flow and showed the importance of released noncondensable gas [16].

In the present formulation, viscosity enters the hydraulic load through coefficients that are held fixed for the water calculations; changing only the tabulated viscosity would therefore not predict its effect on motion. For viscous Newtonian liquids, linear hydraulic resistance and shear-related damping must depend on viscosity, and the loss and local-pressure coefficients must be identified for the appropriate flow regime. Higher viscosity increases viscous dissipation at prescribed motion and generally attenuates motion near resonance, but displacement, phase shifts, and pressure minima require recalculation of the coupled system. The outcome also depends on whether flow rate or available pumping pressure is maintained. The short holes require allowance for entrance losses and unsteady flow. Density, viscosity, and vapor pressure must correspond to the operating temperature. For oils, gas release must be distinguished from vapor formation; a separate gas-content description is needed where degassing is significant [16].

For many non-Newtonian liquids, apparent viscosity depends on shear rate, while yield stress or viscoelasticity may introduce additional effects. Shear thinning can reduce local viscous resistance relative to that estimated using zero-shear viscosity, whereas shear thickening can increase it. Experiments on shear-thickening suspensions demonstrate that rheology also changes bubble deformation and collapse [17]. A measured constitutive law, such as Carreau-Yasuda for shear thinning or Herschel-Bulkley for yield-stress behavior, must therefore be incorporated into the hydraulic-load description; significant viscoelasticity requires stress evolution as well. The vapor-fraction closure and its growth and collapse times must also be reassessed. Consequently, neither the predicted 41-60 Hz pressure-threshold interval nor the operating compromise near 50 Hz can be assumed valid for other fluids without further modeling and experiments.

6. Conclusions

A compact coupled electromechanical-hydrodynamic model of the low-frequency cavitator has been developed. The formulation accounts for the variable air gap, magnetic saturation, motion-induced electromotive force, liquid added mass, nonlinear pressure losses, and finite-rate vapor-content evolution. For the adopted parameter set, the model predicts a pressure-based cavitation interval of approximately 41-60 Hz, whereas cavitation becomes practically negligible at 100 Hz. At 50 Hz, the predicted grate amplitude is 2.21 mm, the minimum local pressure is 41.2 kPa, the peak vapor fraction is 4.8 %, the maximum acceleration reaches approximately 22.5g, and the RMS coil current is 2.81 A. The current and electromagnetic-force phase lags relative to the voltage are approximately 29.6° and 59.7°, confirming that representation of the electromagnetic excitation by a prescribed in-phase harmonic force is inadequate. The predicted total electrical input to the two drives is approximately 0.65 kW, while the dissipative component of the hydrodynamic load corresponds to an equivalent damping of about 138 N·s/m. Compared with the previous force-driven model [11], the predicted 50 Hz displacement amplitude is approximately 21 % lower because the actual fundamental electromagnetic-force amplitude is limited by current dynamics, variable air gap, and magnetic saturation. Thus, operation near 50 Hz appears to provide a practical compromise between pressure reduction, grate displacement, and air-gap safety for the adopted parameter set.

Published experiments on related cavitation devices support the feasibility of the treatment principle, while direct quantitative validation of the present coupled model remains to be performed. The theoretical predictions will be experimentally tested in subsequent research. Application to viscous oils and non-Newtonian liquids requires fluid-specific hydraulic coefficients, appropriate rheological relations, and reassessment of the vapor and gas dynamics.