[rolando]

Hi,
I had some troubles experimenting with moving meshes using an adjustable time step with the icoDyMFoam solver.
If I use a timestep, which is far from the "optimal" one at the beginning of the calculation the timestep doesn´t "converge" at all, but sometimes changes many orders of magnitude from step to step.

Watching the meshPhi fields showed that this fields don´t contain the correct values of mesh movement.

I tried to track the problem and I think it is located in fvc::meshPhi.
As I use the backward timescheme, I can only talk about backwardDdtScheme<type>::meshPhi().

I think the problem is caused when deltaT and deltaT0 have very different size, which may happen at the beginning of a simulation.

The whole calculation in this method seems somehow inconsistent.
The methods return is:
(1 + dT/(dT + dT0))*phi - (dT^2/(dT0*(dT + dT0)))*phi0

If I look at the extreme case with constant meshPhi (phi=phi0) this method produces incorrect results.

The return now is:
(dT0^2+2*dT*dT0-dT^2)/(dT0*(dT+dT0))*phi

If I now say dT = 1.2*dT0 (this is the upper limit, set by setDeltaT.H) the error is about 11%, which is quite a lot.
If dT -> 0 which can also occur, as there is no limiter for lower bounds the error can even "explode".

Would it be better to use the uncorrected value fvMesh::phi() instead of fvc::meshPhi()?
If so, why then even do the above meshPhi correction?

Rolando

[henry]

Does your case run correctly if you use Euler implicit rather than the backward temporal scheme?

Henry

[rolando]

Hi Henry,
after quickly looking at the case, it seems to work with the Euler scheme.

Rolando

[henry]

OK, I will look more closely at the implementation of fvc::meshPhi for the backward scheme; it's not well tested because for moving-mesh problems we usually use Euler because the Courant number is limited by stability issues and Euler is usually accurate enough.

Henry

[rolando]

Thanks Henry.
I think the first order Euler scheme is not the correct choise for what I´m planning to do. I hope to do some LES calculations with a moving grid.
Would it meanwhil be incorrect to use the uncorrected fvMesh::phi() instead of fvc::meshPhi()?

Rolando

[henry]

Rolando,

You can do LES with the Euler implicit scheme if your Courant number is sufficiently small although it is preferable to use either backward or Crank Nicholson with a small amount of off-centering. Using fvMesh::phi() directly is not consistent with using the backward scheme for a case in which the time step changes which is why fvc::meshPhi() was introduced. I will test the backward differencing implementation of this function and post the correction if necessary.

Henry

[rolando]

Hello Hrvoje,
is it possible to get a patched version of OpenFOAM-1.3?
This would help me to continue my work, while waiting for 1.4.

Rolando

[rolando]

Henry,
what do you mean with a sufficiently small Courant number?
I use max Co of about 1.5. Therefore I use an additional nonlinear loop around all equations.

By the way, what small amount of off-centering do you mean?

Rolando

[henry]

Certainly less than Co < 0.5 would be necessary to get good results using the Euler scheme. With a Courant number as high as 1.5 backward may be your only choice because you will probably find Crank Nicholson unstable. If you want to try it use a coefficient of 0.9 which corresponds to using 10% Euler, we find that is pretty much as stable as backward and is useful for equation systems in which sub-cycling precludes the use of backward.

Henry

[rolando]

Henry,
that´s also what I noticed. The Cranck Nicholson scheme is very unstable at that high Co.
But I will try it with the coefficient of 0.9.

Do you mean the nonlinear loop around the equations with sub-cycling?


Hrvoje,
thanks for your eMail.

Rolando

[henry]

No I mean sub-cycling with in the time-step, e.g. in the interFoam codes the gamma equation is sub-cycled and this is consistent with the rest of the equations for Euler and Crank-Nicholson but not with backward.

There is a problem with fvc::meshPhi() for backward and I will fix it and post the fix after testing.

Henry

[rolando]

Thanks a lot Henry,

Rolando

[rolando]

Hello Hrvoje,
I´m still compiling your development version.

1) I had problems with the files skewLinear.C and leastSquaresGrad.C. They both use the method mesh(), which is not yet implemented in the base class MeshObject. I changed the calls from mesh() to mesh_.

2) I get some errors compiling tetDecompositionMotionSolver. I´ll delete all lnInclude´s (I didn´t do so far) and recompile everything. I tell you tomorrow, if I had success.

Rolando

[henry]

Rolando,

I have now checked the backwardDdtScheme<type>::meshPhi function in version 1.3 carefully and tested it and I don't think there is any error, the coefficients are the way they are to be second-order accurate with time-step variation. If you are having stability problems when your Courant number is high and your time-step varies rapidly this may be fundamental to the backward scheme and that you will either need to reduce the time-step, reduce the rate of change of the time-step or perhape use Euler implict during the start-up phase of the calculation and change to backward after the flow has evolved and the time-step stopped changing so rapidly. Play around with these options to see if you can stabilise your calculation.

Henry

[rolando]

Hello Henry,
thanks for that.

Rolando

[rolando]

Hello Hrvoje,
I´m now using your new version of backwardDdtScheme.C and I´m experimenting quite a lot with it.
It works better than the original version, but the problem is not cured.
fvc::meshPhi now produces results that are about the same order of magnitude of what they should be. But they still can differ quite a lot.
They seem to oscillate around the target value.
So the problem is still the same:
- Adjusting the timestep causes a wrong prediction of meshPhi
- This will cause a change in timestep
...

When I start calculation with a fixed and "correct" time step and I switch to adjustable timestep afterwards, the oscillation also occurs, after some delay.

The absolute fluxes in contrast have correct values.

Rolando

[henry]

Rolando,

You seem to be missing the point of fvc::meshPhi: it is a mechanism to generate "effective" mesh fluxes corresponding to the mesh motion such that the order of the temporal discretisation scheme is maintained during mesh motion; they should not necessarily be the same as current time-step mesh fluxes obtained from mesh.phi() except for the Euler scheme. Perhaps the problem you are seeing is not directly associated with fvc::meshPhi but that this affects the relative fluxes which are then used in the time-step evaluation. Try using the absolute fluxes in the time-step evaluation and see if that removes the instability the current approach seems to be exciting.

Henry

[rolando]

Thanks Henry,
if I got you right, it seems to be a stability problem inherent in the backward scheme.
The Euler scheme as posted above works well, even the Crank-Nicholson scheme works well (for this very simple test case).

Evaluating the absolute fluxes helps and the case runs absolutly stable and correct.
But for evaluating flows with moving grids the relative fluxes should be evaluated.
Maybe the way of time-step evaluation can be changed for stabilization reasons. E.g. using a proper damping mechanism or using fvMesh::phi() instead of fvc::meshPhi() for the evaluation.
What do you think about it?

Rolando

[henry]

It isn't clear that there is stability problem inherent in the backward scheme except perhaps in the case of moving meshes with variations in time step derived from the relative fluxes; this is a very special case. If this is so then you need to find a way of avoiding this feed-back mechanism from the high-order mesh fluxes to the time-step and back to the high-order mesh fluxes, one way being using the absolute fluxes or as you suggest approximate relative fluxes for the evalauation of the time-step adjustment.

Henry

[rolando]

Thanks Henry,
I´ll try to find a proper time-step adjustment for that kind of problem.

Rolando