[chriss85]

Dear Foamers,

I'm currently trying to model an ablation process caused by radiation and heat conduction. I'm calculating the ablation density rate at a patch and add it to the rho equation as a field with values in the cells next to the patch. My solver is based on the pressure-based sonicFoam solver using OF 2.3.x.

I now wanted to do some studies on the required mesh and time step size for generating a correct solution and I'm noticing weird behavior.

I created a one dimensional case with an ablating wall on one side and a fixed pressure on the other side. Temperature is also fixed, so that incident radiation falls on the wall and causes ablation (and keeps the temperature steady).

After some time the solution converges, so that the velocity, pressure and temperature fields don't change anymore, however the pressure is constant over the whole domain and the velocity is zero.

This is not what I have been expecting though: If density is added in each time step, the pressure at the wall should be higher and there should be a net flow away from the wall.

I have tried to look at the terms of the density equation and it looks as if my source term is not increasing the density anymore, even though it is not zero. Multiplying this ablation rate with the time step also produces a relevant density that should be added.

Now I don't know how the boundary conditions for rho are generated in the code. Can anyone explain to me how the boundary values for rho are calculated and if there could be a problem here?

I have changed the mesh size to get enough spatial precision compared to smaller meshes.

I noticed that when I increase the PIMPLE corrector iterations, the solution becomes unsteady and temporal oscilations occur in pressure and velocity.

This works:
nOuterCorrectors 5;
nCorrectors 6;

While this doesn't:
nOuterCorrectors 10;
nCorrectors 10;

The same happens when the Courant number is too high, 0.1 works while 1 produces an oscillating behavior. I'm currently using upwind schemes.
In the case of a converging solution, I also see oscillating behaviour before, but it converges to a constant solution.

Any idea why? From what I've read in the forum more iterations (and ofc lower Courant numbers) should help with the stability.

I have also thought about creating boundary conditions for p or U to model the ablation, but I don't know of a reasonable approximation to find the correct impulse generated by the ablated mass. Putting it in the first cell should produce a pressure increase and thus a flow away from the wall (although the pressure decrease due to cooling acts against it).


I'm attaching some images with results so you can see what I mean. The ablating wall is on the right.

Is there anything obvious I'm missing? I'd be grateful if someone can provide me with a few hints