
[Sponsors] 
March 28, 2007, 07:08 
Hi,
I had some troubles exper

#1 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
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*dT0dT^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 

March 28, 2007, 07:15 
Does your case run correctly i

#2 
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 14 
Does your case run correctly if you use Euler implicit rather than the backward temporal scheme?
Henry 

March 28, 2007, 07:26 
Hi Henry,
after quickly looki

#3 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
Hi Henry,
after quickly looking at the case, it seems to work with the Euler scheme. Rolando 

March 28, 2007, 07:30 
OK, I will look more closely a

#4 
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 14 
OK, I will look more closely at the implementation of fvc::meshPhi for the backward scheme; it's not well tested because for movingmesh problems we usually use Euler because the Courant number is limited by stability issues and Euler is usually accurate enough.
Henry 

March 28, 2007, 07:42 
Thanks Henry.
I think the fir

#5 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
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 

March 28, 2007, 07:51 
Rolando,
You can do LES wit

#6 
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 14 
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 offcentering. 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 

March 28, 2007, 07:52 
Hello Hrvoje,
is it possible

#7 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
Hello Hrvoje,
is it possible to get a patched version of OpenFOAM1.3? This would help me to continue my work, while waiting for 1.4. Rolando 

March 28, 2007, 07:59 
Henry,
what do you mean with

#8 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
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 offcentering do you mean? Rolando 

March 28, 2007, 08:06 
Certainly less than Co < 0.5 w

#9 
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 14 
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 subcycling precludes the use of backward.
Henry 

March 28, 2007, 08:21 
Henry,
thatīs also what I not

#10 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
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 subcycling? Hrvoje, thanks for your eMail. Rolando 

March 28, 2007, 10:09 
No I mean subcycling with in

#11 
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 14 
No I mean subcycling with in the timestep, e.g. in the interFoam codes the gamma equation is subcycled and this is consistent with the rest of the equations for Euler and CrankNicholson 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 

March 28, 2007, 10:14 
Thanks a lot Henry,
Rolando

#12 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
Thanks a lot Henry,
Rolando 

March 28, 2007, 12:18 
Hello Hrvoje,
Iīm still compi

#13 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
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 

March 28, 2007, 14:22 
Rolando,
I have now checked

#14 
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 14 
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 secondorder accurate with timestep variation. If you are having stability problems when your Courant number is high and your timestep varies rapidly this may be fundamental to the backward scheme and that you will either need to reduce the timestep, reduce the rate of change of the timestep or perhape use Euler implict during the startup phase of the calculation and change to backward after the flow has evolved and the timestep stopped changing so rapidly. Play around with these options to see if you can stabilise your calculation. Henry 

March 29, 2007, 04:48 
Hello Henry,
thanks for that.

#15 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
Hello Henry,
thanks for that. Rolando 

March 30, 2007, 04:31 
Hello Hrvoje,
Iīm now using y

#16 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
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 

March 30, 2007, 04:50 
Rolando,
You seem to be mis

#17 
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 14 
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 timestep 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 timestep evaluation. Try using the absolute fluxes in the timestep evaluation and see if that removes the instability the current approach seems to be exciting. Henry 

March 30, 2007, 05:34 
Thanks Henry,
if I got you ri

#18 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
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 CrankNicholson 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 timestep 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 

March 30, 2007, 05:40 
It isn't clear that there is s

#19 
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 14 
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 feedback mechanism from the highorder mesh fluxes to the timestep and back to the highorder mesh fluxes, one way being using the absolute fluxes or as you suggest approximate relative fluxes for the evalauation of the timestep adjustment.
Henry 

March 30, 2007, 05:49 
Thanks Henry,
Iīll try to fin

#20 
Member
Rolando Maier
Join Date: Mar 2009
Posts: 89
Rep Power: 9 
Thanks Henry,
Iīll try to find a proper timestep adjustment for that kind of problem. Rolando 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Bug caused by CrankNicholson scheme  cbeck  OpenFOAM Bugs  8  March 28, 2009 02:48 
waves caused by stones  gamego  Main CFD Forum  5  October 1, 2008 10:33 
Forces caused by fluids  anja  OpenFOAM Running, Solving & CFD  20  May 8, 2006 06:37 
Adiabatic compression caused by gas ram  Chris  FLUENT  0  December 13, 2005 05:42 
problem caused by the different materials  limingtiger  CDadapco  1  November 21, 2005 09:02 