Regarding a Track of Proper Liquid-gas Interfacial Charge Density Calculation for EHD Solution

[Nazmul Hassan]

Dear Concerns,

I am conducting an axisymmetric VOF simulation of electrically influenced bubble bursting using Basilisk. I use:

#include "ehd/implicit.h"
#include "ehd/stress.h"

The liquid is electrically conducting, while the gas conductivity is zero. I currently associate rhoe with the liquid VOF phase using:

event defaults (i = 0)
{
 rhoe.inverse = false;
 f.tracers = list_append (f.tracers, rhoe);
}

I do not include tracer.h, and rhoe is not registered in the ordinary BCG tracer list.

The electrical properties are updated using:

#define perm(T) (1./(((T)/BETA_12) + (1. - (T))))
#define cond(T) (COND_1*(T))

event properties (i++) {
 face vector cf[];
 face_fraction (f, cf);
 foreach_face() {
 double ff = clamp ((f[] + f[-1])/2., 0., 1.);
 double fc = clamp (cf.x[], 0., 1.);
 epsilon.x[] = perm(ff)*fm.x[];
 K.x[] = cond(fc)*fm.x[];
 }
}

For a fresh simulation, I initialise:

foreach()
{
 rhoe[] = 0.;
 phi[] = VOLTnd*(X0 + L0 - x)/L0;
 foreach_dimension()
 u.x[] = 0.;
}
 boundary ({rhoe, phi});
boundary ((scalar *){u});

The simulation remains numerically stable, but I observe a persistent rhoe contour corresponding to the initial bubble meniscus, even after the liquid gas interface has moved considerably. The old contour is visible in cells where the current volume fraction is (f=0).

For example, at (t_sim=0.54), close to the first jet droplet pinch-off, my integrated diagnostics are:

qabs = 66.7774
qgas = 47.2055
qint = 19.3510
qliquid = 0.2194

The classifications used were:

if (f[] > 1e-4 && f[] < 1. - 1e-4)
 qint += fabs(rhoe[])*vol;
 if (f[] <= 1e-12)
 qgas += fabs(rhoe[])*vol;
else if (f[] >= 1. - 1e-12)
 qliquid += fabs(rhoe[])*vol;

Here, qgas is the integral of fabs(rhoe) over cells satisfying (f\le10^{-12}). Therefore, approximately 70.7% of the absolute charge is located in cells classified as pure gas. Multiplying the output by f removes the silhouette visually, but I understand that this would only mask the issue and would not correct the charge evolution.

For visualisation, I also output:

rhoe_liquid[] = rhoe[]*f[];

This hides the old contour, but I understand that it is only a display variable and does not correct the charge used by the electrical solver.

I would be grateful for guidance on the following questions:

1. Is it valid to register rhoe in f.tracers while using ehd/implicit.h, or is rhoe subsequently reconstructed by the implicit electric-potential step in a way that overrides the VOF-constrained transport?
2. Is rhoe intended to represent a volumetric free-charge density distributed across interface cells, rather than the physical surface-charge density?
3. If I require the interfacial free surface-charge density, should I calculate it from the electric-displacement jump, σs​ = n⋅(ϵg * ​Eg​ − ϵl​ * El​) where the normal points from liquid to gas? If so, what is the recommended way to evaluate the two one-sided electric fields on an adaptive VOF interface in Basilisk?
4. How should interfacial charge be transported conservatively with the moving VOF interface while also accounting for Ohmic conduction?
5. Is charge in gas cells immediately adjacent to the reconstructed interface expected numerically? If so, what diagnostic should be used to distinguish the numerical representation of surface charge from an unphysical residual charge left behind by the moving interface?

Any advice about the correct Basilisk formulation or a suitable validation case would be greatly appreciated.

Kind regards,
MD Nazmul Hasan
MSc in Advanced Mechanical Engineering
Queen Mary University of London