Hi,
This is mainly addressed to Henry:
(apologies for lengthiness...no need to urgently reply to this though!)

I had to give up...I've been trying to understand the very innocent seeming BISO loop in the mhdFoam.C solver. My own MHD codes have not used BISO type corrections, but I nevertheless thought it wouldn't be something I'd have to spend much time on in FOAM, however I simply cannot figure it out to my absolute satisfaction.

Naively I would think that because BPISO adjusts "phiB" it also indirectly destroys the accuracy of the momentum equation, yet the momentum equation is not resolved in the BISO loop, only "phiB" gets corrected. It seems analagous to an incompressible PISO loop because it is enforcing div(B)=0, as incompressibly one enforces div(U)=0 in deriving the pressure correction equation. Only the big difference is that we know div(B)=0 is automatically guaranteed by a physically accurate solution, so it should not need extreme correction as the pressure requires in PISO.

Now, for compressible MHD flow BISO is unchanged because one still wants div(B)=0, so I currently have a BISO loop in my compressible MHD code that has the same form as the incompressible BPISO loop. Will that be OK?

Those are the general remarks. My specific query has to do with velocity-magnetic field coupling.

Taking incompresssible PISO as a lesson on how to do all this in FOAM, I would expect there to be a velocity correction after the updating of the flux, as indeed there is in the outer loop of ico-PISO. But in mhd-BISO there is only a magnetic field flux correction

"phiB=phiB-pBEqn.flux();"

no velocity correction follows. But the B-transport equation is re-solved, which seems logical to me because after all div(B)=0 should be already satisfied by a well-behaved solution, so we only require weak coupling here, and hence we can forgo velocity and pressure correction. But this it not the rigorous rationale that I seek!

The question I have is: has this been analytically proven somewhere to lead to an overall consistent segregated solver? Or is the BISO correction loop a seat-of-the-pants rough correction attempt and not analytically rigorous. In relation to this, the appearance of the scalar field "pB" seems to be entirely a dummy variable to temporarily store the error in div(B)=0. I'm assuming that's why in the hartmann case (from the tutorial package) it is initialized to zero everywhere.

But then why is the coefficient "1.0/Ueqn.A()" used in the BISO pBEqn for the correction to phiB? Has this got soemthing to do with my worry about recoupling to the velocity somehow? e.g., one can derive a PISO type equation for the magnetic field correction if the momentum equation was rewritten with the grad(magSqr(B)) terms dropped off, to get a psuedo-velocity. But that sort of approach would seem to me naively to require explicit velocity correction because it would be changing B or phiB at least, based upon the coefficients in the discretized UEqn, which I suspect is what you are trying to avoid? So again that leads me to ask why does "1.0/UEqn.A()" make an appearance in the BISO loop pBEqn?

A related question is does the "solve" on the "ddt(B)+..." B-transport equation change the velcotiy field, which afterall appears in that equation? The question is, does the FOAM solver always solve for all fields that appear in PDE expresions that are "solve"d? If so, then this confounds the accuracy of the solution to the momentum and mass continuity equations which may have repurcussions on the enxt time step. Has the segregated solver emthod been rigorously justified in these cases of multiply coupled sets of equations, or am I better off waiting until an OpenFOAm contributer delivers a robust block matrix solver set of OOP features?

Currently I'm thinking that for MHD I'll have to treat the equations as explicitly solved (as a set) and therefore requiring a strong CFL limit so that none of the coupled fields gets too far out of synch on a given time step.

Clearly I have not fully analyzed the mhd-BISO method fully or carefully, otherwise I wouldn't be baffled by all this. But I cannot contribute much to OpenFOAM without understanding some of these things that are not covered in Hrvoje's thesis.

So far whenever I've worried about FOAM it has been my understanding that's been astray, not the standard solvers. So once again I beg your patient indulgence.

Yours FOAMcerely,
Blair.

The first thing I must stress is that I am not an MHD expert.
I wrote mhdFoam one Sunday afternoon after buying a small
book on MHD from the Imperial College library which looked
interesting; it cost me 50p :-). What struck me most about
the presentation of the B equation in the book was the
similarity to the momentum equation so I decided to create an
algorithm based on PISO to solve it. I have not read any
papers on MHD or numerical solution of MHD problems so I have
no idea if the approach I have taken is common or used at all,
it just seemed like a good idea to me.

My main aim was to create a method to solve both for the
cell-centre B as well as face-fluxes for B as PISO does for
U and apply the same kind of Rhie and Chow type approch to avoid
staggering on a collocated mesh. This requires the solution of a
"pressure" equation for B so I introduced a fictitious magnetic
field pressure and used it to derive a flux equation from the
div(B)=0 constraint. The consequence of this is that the fluxes
exactly obey this constraint but the fictitious magnetic field
pressure is non-zero, in fact it holds the non-convergence error
which tends to zero as the system of equations are iterated to
convergence. However, it need not actually reach zero, it will
hold a form of discretisation error representing the difference
between the magnetic field fluxes and the face-interpolate of the
cell-centre magnetic field. Tests have shown this error to be
small and hence the fictitious magnetic field pressure to be small.
In conclusion the fictitious magnetic field pressure is only
introduced to obtain an efficent solution algorithm.

> Now, for compressible MHD flow BISO is unchanged because one
> still wants div(B)=0, so I currently have a BISO loop in my
> compressible MHD code that has the same form as the
> incompressible BPISO loop. Will that be OK?

That sounds correct although I am not sure what form the coupled
term takes for compressible flow; does density not feature in the
B equation at all?

> Taking incompresssible PISO as a lesson on how to do all
> this in FOAM, I would expect there to be a velocity
> correction after the updating of the flux, as indeed there
> is in the outer loop of ico-PISO. But in mhd-BISO there is
> only a magnetic field flux correction

That's because the fictitious magnetic field pressure does not
feature in the B-equation and so cannot be used to correct B in
the same way the pressure is used to correct U although looking
at it now I think I could come up with a way.

> no velocity correction follows. But the B-transport
> equation is re-solved

I don't follow you, where is the B-transport eqn resolved?

> The question I have is: has this been analytically proven
> somewhere to lead to an overall consistent segregated
> solver?

Nope, it was just my bit-of-fun on a Sunday afternoon which
happens to work really well for the cases we have run.

> Or is the BISO correction loop a seat-of-the-pants
> rough correction attempt and not analytically rigorous.

It's not analytically rigorous but more than seat-of-the-pants
stuff; I do have quite a bit of experience with numerics and
solution algorithms.

> But then why is the coefficient "1.0/Ueqn.A()" used in the
> BISO pBEqn for the correction to phiB?

It should be 1.0/BEqn.A() but in fact for most cases they will
be identical. However, for consistency I will change this.

> A related question is does the "solve" on the "ddt(B)+..."
> B-transport equation change the velcotiy field, which
> afterall appears in that equation?

No, it's a segregated solver.

> Has the segregated solver emthod been
> rigorously justified in these cases of multiply coupled sets
> of equations, or am I better off waiting until an OpenFOAm
> contributer delivers a robust block matrix solver set of OOP
> features?

Segregated solvers are fine and correct and efficient if they
converge. Block-solvers are best for cases where segregated
solvers either don't converge or convergence is very slow.
Does the segregated solver not work for your cases?

> Currently I'm thinking that for MHD I'll have to treat the
> equations as explicitly solved (as a set) and therefore
> requiring a strong CFL limit so that none of the coupled
> fields gets too far out of synch on a given time step.

Why do you think so? Why not try the segregated approach to
see how well it performs for your cases? If you are worried
about better coupling between the variables why not try putting
an outer-loop around the equation-set.

Henry

Thanks Henry,
You're a gem.

It is correct that density does not feature in the B-equation, it is derived purely from Maxwell's equations! It is, if you like, a reduction of Maxwell's equations to a single equation, with appropriate "MHD Approximations" such as the weak form of Ohm's law, charge neutrality, and neglect of displacement current, to get rid of all explicit references to E and current density J.

I liked your ico-mhdFoam code, which is why I haven't yet tried any more drastic departures for the compressible MHD case. I did not mean at all to take my ignorance of FOAM and CFD to impute that you have less than expert experience! Your numerics do seem soundly reasoned to me, and like you say, there is proof in the pudding of some sort, at least for 'mhdFoam'.

I hope (for a thrill) that I can add to OpenFOAM's MHD features before some other real MHD experts get in on the act, but if someone does take over and release a great MHD solver then I'll gladly take a back seat. My strengths are not in numerics and algorithms, but rather in weakly ionized plasma properties, and so I hope to add some code to allow a user to model nearly real partial plasmas with temperature dependent gas conductivity and so forth.

I'll take the rest of your notes on board and conduct further experimentation. I'll try to analyse the robustness of the numerical algorithms if I can, time and patience willing, but that is not my forte.

One of the main potential problenms I see for the codes I'm developing are the often vastly different time scales involved in MHD systems, e.g., when strong B-field gradients exist the magnetic field diffusion can cause numerical instabilities. That's a problem for someone trying to simulate electromagnetically driven shocks for example, because you need really strong initial fields of 5 to 100 Tesla maybe, to generate well-formed MHD shocks, while the gas velocity is still subsonic. So at some stage I will have to query the notice board on how to go about calculating all the waves speeds (acoustic as well as Alfven, fast MHD shock, and various signal speeds associated with all forms of diffusion etc.) from the primitive variables over a mesh, and how to include them in a modified generalized "CourantNo.H" code to get a maximum time step for guaranteed stability, then compare with less restrained CFL number and see how the code performs.

I do not have a vast amount of experience with these issues, so all FOAMer's please don't hold back making some suggections...

If anyone else is reading this, please do not count on me to deliver a new FOAM MHD code, because I've been forced at present to do this all on my weekend time, as my research has suddenly shifted directions (by necessity not desire) away from MHD and towards fiddly experimental work in materials science.

-Blair.

Hi,

I am trying to use the mhdFoam and I have two questions:

1. What is the 'pB' ?
2. How I can impose 'dB/dx'?

Thank you,
Rita

Hi,

Here are some precisions about question 2:

It seems that the boundary conditions we can impose on the magnetic field B are only of Neuman type (B constant). Is it possible to impose a Dirichlet type condition (given rotational of B, or imposed current)?

thank u,
Rita

[ed_teller]

Hi everyone,

First i would like to give my compliments to Weller and Hrv Jasak for their outstanding work in developing such a powerful (and extendable) tool as OpenFOAM.

My main issue in this message is:

1) Would it be difficult to develop/implement a Solver for the compressible MHD CFD problem in a hot dense plasma, with radiation transport assuming a variable opacity (as a function of state vars in plasma) and multi spectral blackbody radiation?

2) What could be used as a coupling variable, between CFD and radiation transport process ( I assume the possibility of ocurrence of a Marshak wave in the problem) for the sake of numerical stability?

3) Does the moving mesh ability of OpenFOAM allows me to keep track of a boundary between two different media in case of the ocurrence of a Richtmeyer-Meshkov instability (in case of a cylindrical cumulative implosion problem, driven by MHD Lorentz force) or a Rayleigh Taylor Instability ( in case of a radiation shock dominated spherical implosion process, driven by an outer spherical shell ablation) ?

Thanks in advance :-)

[giovanni10]

Hi everyone,
I am just a beginner in openFoam.
At the moment, I am trying to add the term JxB in sonicFoam. I have seen that it is descriping by two terms -div(phiB,2.0*DBU*B) and grad(DBU*magSqr(B)) in MHDFoam which I put them in UEqn.H in my new Folder my_sonicFoam. I added the magnetic Field in creatField.H in OpenFoam/solvers/my_sonicFoam, DBU and #include "createPhiB.H" and I typed wmake and I have errors.
I would appreciated if you could help me.
Best wishes

[giovanni10]

Does anyone know if it is possible to use LES ( LARGE EDDY SIMULATION ) in mhdFoam? Is there any tutorial?
THANKS !

[levka]

Hello guys,
do you know where can i find mhd solver with electric potential formulation that utilize Four Step Projection Method (for computing incompressible conducting flow at high Ha ) ?