Suggested unsteady, implicit solver stable with arbitrarily large time steps
1 Attachment(s)
Dear all,
as already discussed in the forum "Questions & Answers related to OpenFOAM®", in the article thread "pimpleFoam, transientSimpleFoam, simpleFoam, turbFoam, oodles, transientSimpleOodles", in the XING OpenFOAMgroup, the CFLnumber in most of the transient solvers provided with the OpenFOAM library is not required, and to limit SIMPLE to steady simulations is also wrong. Moreover, in these solvers, you do actually not get convergence on an (arbitrarily) large time step and, therefore, experience convergence problems for larger time steps. What I consider important in order to achieve an implicit procedure with arbitrary timestep sizes is that all discretization is done in an implicit manner (fvm::....) and that the sequentially solved equations are solved for _several_ times in a single time step to reach a converged _SYSTEM_ of equations on a time step. Below and attached you will find an early development version of our code with fixed convergence criteria and some features that still need some improvements. However, the basic idea should become apparent. With this code you can achieve large time steps, and if the initial residuals in the inner loop below the timestep loop do not decrease, the timestep is probably physically too 'coarse' or too large and should be decreased. Suggestions and discussions are welcome. We want to discuss this on the Stammtisch in Berlin tomorrow, for anyone interested in the bofsession 'Hybrid RANS/LES and DES'. Best regards, Ulf. Code:
/**\ 
thanks for posting these additional (but VERY necessary) snippet. i will try to set them up in the baratropicCavitatingFoam and report back.
neewbie P.S.: nice workshop in Berlin ;) 
addon: Suggested unsteady, implicit solver stable with arbitrarily large time steps
1 Attachment(s)
Dear all,
as we were not able to present the energyspectra (Engergy vs. wavenumber) for the computations of decay of isotropic turbulence / box turbulence we made on the OpenFOAMStammtisch on Friday for different settings they are given here. The Fouriertransformation binary did not work on our laptop and we did not get a fortran compiler in a few minutes due to a lack of online access that day. Instead, we presented curves for the temporal evolution of point probes for instantaneous velocity at a given point. Significant differences became apparent. These results are confirmed here and now. Please take a look at the picture attached to this post. What do we see in the attached diagram? Energy vs. wave number at a given point of time / after simulating 2 time units. The basics of the results, how the Fourier transformation was done, and the comparability to available experimental data is of no interest for the following conclusions. The results itself are not very good, admittedly, as a very coarse grid was used. The simulations are done with the source code attached to the first post in this thread. We applied different underrelaxation (no underrelaxation: RELAX=1, and 50% underrelaxation: RELAX=0.5). Then we computed one loop per time step ('one loop') and a variable number of loops on a time step until a certain residual was reached ('more loops, MaxRes=1e20). And we varied the timestep size using DT=0.2 and DT=0.01. Thus, we get 8 simulations. What do we conclude?  For 'large' time steps with relaxation (red curves), the results significantly depend on whether you couple the equations by more loops on a time step or not.  For 'small' time steps (blue curves), this is also true.  Without relaxation (green and magenta) the influence of not coupling the equations decreases with decreasing timestep size. The magentacolored curves are even congruent  With coupling of equations on a timestep the results get 'improved' (are closer to each other), see the dashed lines. > Coupling on a time step (inner loops, controlled by residuals) is required. Another effect, but a totally different issue, is the dependence of the results on relaxation. Does this relate to the RhieandChow interpolation? Has anyone tried to evaluate that? Comments are welcome, but, however, please not on the curves themselves. We know that we need much finder grids for decent results especially when applying full central differencing. Anyway, the effects are the same as described here. We have checked that in other cases. Best regards, Ulf. 
1 Attachment(s)
Dear Ulf,
Thanks for posting your results and conclusions in detail, especially since I couldn't make the workshop and find the topic very important. First of all, I agree with the procedure you have followed for the new transient SIMPLE solver (multiple iterations of the equation system within each time step). To save us having to compare the source code, could you possibly outline the differences between your solver and the transientSimpleFoam and pimpleFoam solvers? Your investigation with isotropic turbulence is very interesting. I find the results partly in line with my intuition, and partly perplexing: Quote:
Quote:
Quote:
Quote:
Quote:
Quote:
So, being sufficiently rattled, I decided to repeat your experiment with ELAN (our inhouse code, also a colocated, FV, "transient SIMPLE" type). My setup was identical except for two things:
Attachment 2282 So, what are my conclusions?
One thing's for sure, this URF dependency is highly undesirable and it needs to be fixed. Has anybody seen this with one of the standard solvers? All the best, Charlie. 
Dear Charlie, dear all,
so we do or did have the same expectations, perfect. The dependence of results on under relaxation is nothing new, I guess. Angelo mentioned that in Berlin and searching for 'under relaxation Rhie and Chow' e.g. in google gives many links to related publications. There the effect is described and workarounds are given. However, I am not that deep in the OpenFOAM source code. As you mentioned, we could do that in ELAN, where my roots are  this is now your task :D. However, there you did not observe this effect  lucky you! Regarding your request to give details about the source code changes I will come back to that later. Basically it is only a dowhile loop on the time step that checks for the initial residuals, and, of course, that is very similar as it is done in ELAN. I remember when I implemented the grid deformation algorithm in ELAN together with the ArbitraryEulerLagrange (:cool:) formulation for all equations that many adjustments were necessary regarding residuals and convergence check  but that is too long ago for me to compare exactly :(. Best regards, Ulf. Remark: For this simulation we used OpenFOAM, version 1.6, July 2009 ('unpatched'), employing our modified, LLkepsilonbased turbulence model for hybrid RANS/LES 
Hi Ulf,
Thanks for the reply. Nice that we agree on how the solver should behave (of course!). :) The problem that the results depend on the underrelaxation factor should definitely be pursued further, as it is a serious limitation. Let's stay on the case and invite contributions / advice from the OpenFOAM community! Best, Charlie. 
Hello,
I have some comments, questions and a suggestion :D Quote:
Quote:
Quote:
Quote:
Quote:
My question is: where the solutions in the different cases actually converged? If so, what criterion was adopted to define the convergence? The suggestion is on the implementation of a different kind of timestep adaption procedure, which has been already successfully used in MFIX (www.mfix.org), with a SIMPLElike procedure. The approach they use is very effective at speeding up the simulation, and the idea is the following:
This algorithm was shown to be very successful in multiphase simulations with very stiff problems, like in granular flows, and can be quite easily implemented changing time class and with minor changes to the solvers :cool: Best, 
Dear Alberto,
thanks for your confirmation! As I see: we succeeded, as we wanted to show something obvious and something that one would expect :). A success like that is a decade ago when teaching CFD  in engineering life this is different and all the engineers in the teams in my department and myself never get results we expect, either in crash, statics, dynamics, solid or fluid mechanics or whatever ;). Quote:
Quote:
The residual in for all equations in the simulations was down to a very low level (except for the simulations with a fixed number of inner iterations, of course), therefore the systems of equations can definitively considered as converged! The number of iterations on each time step naturally increases with increased under relaxation, but the level of convergence was equal for all the simulations we did (we have only shown a small excerpt). Nonetheless, the results differed dependent on the relaxation factor, although we did have equally converged solutions. This should not depend on the way you implement your pressurebased procedure. We rather think that this comes from the way pressurevelocity decoupling on nonstaggered grids is avoided by employing a RhieandChow interpolation. Best regards, Ulf. 
Hi,
sorry, I didn't take a deep look to the code before... my bad! The forum makes it quite hard to read with the scrolling box :p Quote:
Quote:
Quote:
Quote:
Quote:
To stabilize the interpolation in those cases, a modified RhieChow interpolation is used, including the effect of the bodyforce term, pressure gradient and other source terms in the interpolation. You can do something similar in OpenFOAM (and they do that for some compressible solver like buoyantPisoFoam). You might want to give it a try to check if it plays a role also in your case. In your code you do not have any body force, so you simply need to modify the momentum predictor as follows UEqn() == fvc::reconstruct(fvc::snGrad(p)*mesh.magSf()) Some additional change is made to the correction step too, since also the velocity correction is reconstructed from the correction to the flux (see pEqn.H in buoyantPisoFoam) There is an additional term coming from the unsteady term that brings a dependency on the relaxation factors (see for example Chio, Numerical Heat Transfer, Part A, 36:545 ± 550, 1999), which, to my knowledge is ignored in OpenFOAM (and in other codes). In the paper they show that the error that you make when you neglect the term dependent on the relaxation factors in the interpolation is important when dt = 0.001 and URF_u = 0.7. I hope this helps a bit, and I'm very interested in the topic (Very happy you pointed it out and shared your results with us. Thanks!), since I've been fighting with a related interpolation problem for a work on multiphase flows I'm doing :D Best, 
Maybe you know it already, but some interesting considerations are in
W. Z. Shen, J. A. Michelsen, N. N. Sorensen, J. N. Sorensen, An Improved SIMPLEC method on collocated grids for steady and unsteady flow computations, Numerical Heat Transfer, Part B, 43: 221239, 2003 In summary

Dear Alberto,
thank you very much for all the information and suggestions! Wow, I am impressed and see that you are very deep in these issues. We will try what you suggested and have to handle all the information you have given. Please give us some time. We will either come back with some more questions or show the improvements we have achieved. Again, thank you very much and best regards, Ulf. 
Quote:
Best, 
Also from me thanks to Alberto for this very comprehensive advice  this will surely help us to overcome this problem!
Best regards, Charlie. 
Dear all,
this is a very interesting topic as I was really surprised to see, that OpenFOAM seems to have problems with such fundamental numerical schemes. Concerning the discussions above, I have some comments: 1. In agreement to Ulf, I am am very sure that all the pressure correction schemes we are talking about (SIMPLE, SIMPLER, PISO) can be applied to steady as well as unsteady flows. There is no reason to restrict them to only one type of time integration. 2. Independent of how the pressure is calculated, there are two reasons why we have to iterate in an inner loop over the system of equations on each time step: a: We solve the coupled equations in a segregated way, which means that the coefficients of one equation are depending on the other equations. b: There is a nonlinear convective term in the momentum equations. This has to be linearised, which requires an iterative procedure even for time explicit schemes. If we do not apply the inner loop, we get an error, which somehow scales with the time step. I am not sure, whether this error tends to zero with decreasing time step. But this might be an academic problem, which is not of major interest, as long as the error is small for small time steps. However, in the above problem (decay of isotropic turbulence) , a Courant number (CFL) of 0.2, which is quite small, still gives a significant error. Note, that this behaviour is very different from explicit time integration, where CFL<1 is required for numerical stability, even if CFL>>1 would be sufficient for numerical accuracy!!! One should also keep in mind that the main advantage of time implicite schemes is that you don't have to pay attention to the Courant number. You simply choose the time step which is appropriate for the physical representation of your problem. For many industrial applications this leads to a very high Courant number. Thus you definitly need the inner loops. 3. Relaxation for the momentum equation means that we take only a part of the new solution of u_i. This is then used in all other equations. Relaxation for the momentum equations can ALWAYS be applied independent of the pressure correction scheme. 4. Relaxation for the SIMPLE algorithm means that the pressure of the last iteration is only corrected by a part of the pressure correction p'. For the correction of the velocities u_i, the FULL pressure correction p' HAS to be applied to achieve mass conservation. For PISO or SIMPLER, you never use relaxation for the pressure, as these schemes are constructed such that you can use the full pressure correction. It is very important not to mix up the above points in order to separate the influence of all the numerical parameters. From my point of view, inner iteration loops are essential for general applicability, no matter which pressure correction scheme is used. This allows the use of much larger time steps and the solution should be independent of the relaxation factor. Some comments on the Rhie&Chow Interpolation: From my knowledge, there are two flaws of the standard formulation. First, the calculation of the face velocities depends on the relaxation factor of the momentum equation. This holds for steady as well as unsteady problems. There exist different cures for this, simple ones as well as more complicated ones. Second, the damping propertie of the Rhie&Chow Interpolation decreases with small time steps. Modifications to avoid this are quite complex. I will post some results of test calculations on this topics today or on monday for further discussion. Best regards, Angelo 
... not a problem of OpenFOAM
Dear Angelo, dear all,
just one comment on your statement: Quote:
Thanks to all of you for this very fruitful discussion. Best regards, Ulf. 
Quote:
Of course, not having the dependency would be better. However I don't think "moving" the dependency from the URF to an arbitrary parameter, as done in one of the paper I read is any better, especially if this makes the method griddependent. :D Quote:
Quote:
For RANS, I agree with you, but the PIMPLE implementation is OK for that. If you are familiar with commercial codes, there is a famous one that iterates for a given number of times, exiting if residuals fall below a certain value, but it does not do absolutely anything if that criterion is not respected, which is not very different from what OpenFOAM does in PIMPLE loops, excluding the fixed number of iterations. Quote:
What test case are you using? Best, Alberto 
1 Attachment(s)
Dear all,
what you mentioned, Alberto, is right: Quote:
Implementing reconstruct for the pressure gradient in the momentum equation (fvc::reconstruct(fvc::snGrad(p)*mesh.magSf())) as suggested earlier in this thread and alternatively or additionally when updating the momentum after pressure iteration did not change the results at all. Some additional test and implementations yielded the following:
Actually, I am not satisfied with these results and suspect that something is overlooked. Now I am really interested, what you all find out. Has anyone checked this for pimpleFoam? Best regards, Ulf. 
Hello and thank you for the update :)
Quote:
Quote:
Quote:
What do you think? Best, 
Moved the thread to this forum instead.
It seems more appropriate. 
Dear Niklas,
Quote:
Hopefully, other contributers will find the thread again easily  if not, that would be most inappropriate ;). Therefore, taking the history of a thread it does not really matter where it resides, I think. It is more important that contributers find it again. For forums or groups that I moderate I suggest and ask before I move a post. Of course, if a thread is totally or obviously off topic or violating rules it is moved or deleted. Well, never mind. Thanks and best regards, Ulf. 
All times are GMT 4. The time now is 20:58. 