[linyanx]

Dear Foamers,

I am currently studying about the metal pad instability. There is a good solver named epotFoam created by Alessandro.

It works on a standard MHD solving procedure. Everything is good but I found one problem while doing the test to see the divergence of the electric current:

Because the MHD's elliptic equation is solved demanding to the charge conservation(divergence operation•J=0). I wrote a test phase into the eqotFoam's original code to plot the divergence of the jn.

//Interpolating cross product u x B over mesh faces
surfaceScalarField psiub = fvc::interpolate(U ^ B0) & mesh.Sf();

//Poisson equation for electric potential
fvScalarMatrix PotEEqn
(
fvm::laplacian(PotE)
==
fvc::div(psiub)
);

//Reference potential
PotEEqn.setReference(PotERefCell, PotERefValue);

//Solving Poisson equation
PotEEqn.solve();

//Computation of current density at cell faces
surfaceScalarField jn = -(fvc::snGrad(PotE) * mesh.magSf()) + psiub;

//Current density at face center
surfaceVectorField jnv = jn * mesh.Cf();

//Interpolation of current density at cell center
volVectorField jfinal = fvc::surfaceIntegrate(jnv) - (fvc::surfaceIntegrate(jn) * mesh.C());

//Update current density distribution and boundary condition
jfinal.correctBoundaryConditions();

//Lorentz force computation
lorentz = sigma* (jfinal ^ B0);

volScalarField test3 = fvc::div(jn);

test3.write();

I tried to run this solver in huntFlow case and let it run for 2 time steps to see the divergence of ‘jn’. In principle, the result is supposed to be uniform zero around the field. But in fact, it was not.The divergence of the jn is far bigger than it is supposed to be. Therefore, there seems to somewhere wrong to me.

Could anybody please correct if my thought for the divergence test is not right? Or there is something I am missing?

Appreciate your time and effort in advance!

Best regards,
Linyan