[wschosta]

Hi all.

I've been working on a solver that couples the functionality of interFoam and mhdFoam, the idea is that when finished it will be able to model the effects of instabilities in the boundary region between two fluids (hopefully eventually liquid metals and plasma).

I've run into a bit of a problem with running the code. I've been running tests using simple fluids, eg. Air and Water, and the solver seems to work pretty well and, as far as I can tell, accurate, if a bit slow to process. I just started running tests using more advanced fluids (eg metal and plasma) and now I've run into some fairly difficult problems with the courant number which I use (as is done in interFoam) to control the time step. Basically what happens is the solver will start and then after one or two iterations the courant number will spike to something ridiculous like 2 million which will cause the time step to reduce itself to something along the lines of 1e-99 and crash the program. I've been able to moderate this by increasing the precision of the solver (in the control dict) but I still have a problem with a high courant number forcing my time step to be small (on the order of 1e-17).

Any thoughts on how to improve these time step values? I've done everything I can think of to optimize my code so that it will be able to run in a reasonable amount of time but I'm not sure what to do about this. To give you an idea of how long this is taking, I started a simulation with the adjusted precision, 24 hours later it had advanced through several hundred iterations but only 6.59183724832e-7 of the simulation.

Thanks in advance!

(P.S. I'm running OpenFOAM 1.7.1 and with the first part of my question I was also wondering if anyone could point me in the right direction on printing out the courant number to a file with the time step after each iteration so that I can keep track of how it changes?)

[marupio]

The Courant number limitation is fundamental to the solver design. I'm not sure how those solvers work, but for reacting flows, it ensures an inherent error is negligible. The species transport does not account for reaction sources, but if the Courant number is less than 1, this is negligible.

Are you sure that it is the Courant number limitation you are running into (i.e. maybe it's a bug somewhere...). If so, can you change the problem to improve your Courant number? Coarsen the mesh, reduce the flow speed?

I wrote a solver for biochemical reactions, and it worked beautifully, except the Courant number restricted my simulation to essentially stagnant reactors, otherwise it would have taken years to produce a solution. I had to implement an entirely different solver design that did not have the Courant number limitation. Ultimately, I had to couple transport with reactions.

I was speaking with some combustion guys at the 6th OpenFOAM Workshop at PennState, and they tell me there's a way to get away from the Courant number limitation by iterating several times in a single timestep. Not sure how that works.

[wschosta]

I think you're exactly right, from what I've found by rooting around the code the courant number doesn't affect anything unless it's greater than 1 (or a values set by the user but 1 is sufficient). When it's greater than 1 there's a short bit of code that, if the user has allowed it, takes into account the current deltaT and adjusts it accordingly. I think it works in the other direction too, if the courant number gets extremely close to zero it will increase deltaT.

It's very possible that there's some bug in the program, this project is my first with OpenFOAM and in the two months I've been working on it I've learned a lot and made a lot of revisions. But I've checked several times and I can't seem to find anything that would indicate an issue. And as far as I can tell it's working correctly (I was able to validate with lower speed tests).

I'm not really sure how I can change my test case, it's a pretty well defined problem that I'm trying to model. There might be a way to deal in proportionality, fluid one moving at ten times fluid two but I'm not sure if that's something that I want to take liberties with. I tried coarsening the mesh and I have shown a bit more success in coming up with some more timely results. But it's not a huge change and I need a fine mesh for what I'm trying to do.

In your experience with you biochemical solver, did you find that the solver without the courant number produced results just as accurate as the solver with the courant number? I'm not really looking to rebuild my solver but it's always a possibility and, if I can't get it working like this, may eventually become a necessity.

Iterating several times within a time step? Ha that sounds like crazy talk! I have no idea how that would work, I know the solvers iterate several times but I'm not sure how you would be able to do something like that in such a way that it would affect the courant number.

Thanks for the help!

[marupio]

Yes, my coupled solver (no Courant number limitation) produces results just as good as the ODE solver(Courant number limitation). As for iterating several times within a timestep, the idea is: if you can somehow include the effect of the reaction source term in the transport equation, you can avoid the Courant limitation. I just don't see how you can do that without a "coupled solver".

You may need to speak with someone who is more familiar with CFD solver design to figure out how to get away from the Courant number limitation.

[wschosta]

That's interesting. And something to definitely look into.

it turns out that the scope of the problem set is on the nano-scale in terms of time. So as opposed to the miliseconds I was thinking (and thus billion iterations to do a single worthwhile time step) it's on the order of a thousand time steps. I'm going to keep optimizing and digging through the code because I wouldn't be surprised if I made an error somewhere that's causing this problem. But for now I think it'll be okay.

[tom_flint2012]

Hi,


Trying to revive a long dead thread here.


I'm doing something similar. I'm interested in the interaction between a plasma and liquid metal.



I've combined interfoam and mhdfoam, but I'm not certain of the changes I need to make to the equations. I went through the derivation of the magnetic induction equation, but assuming position dependent electrical and thermal conductivity, and I got a new equation that I implemented in my BEqn.


This new equation has some extra grad(permeability) ^ J and grad(sigma) ^ J terms. I cant seem to find any literature on this, does anyone know how we can derive the magnetic induction equation with position dependent conductivity and permeability?


Thanks in advance,


Tom