[djbungee]

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 OpenFOAM-group, the CFL-number 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 time-step 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 time-step loop do not decrease, the time-step 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 bof-session 'Hybrid RANS/LES and DES'.

Best regards, Ulf.

 /*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Application
    transientSimpleOodles

Description
    Incompressible LES solver using transient SIMPLE.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "incompressible/transportModel/transportModel.H"
#include "LESModel.H"
#include "IFstream.H"
#include "OFstream.H"
#include "Random.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"//#include "createMeshNoClear.H"
    #include "createFields.H"
    #include "initContinuityErrs.H"
    #include "readTransportProperties.H"

    Info<< "\nStarting time loop\n" << endl;

    for (runTime++; !runTime.end(); runTime++)
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "readPISOControls.H"
        #include "CourantNo.H"
    scalar eqnResidual = 1;
    scalar convergenceCriterion = 1.E-20;
    scalar maxResidual = convergenceCriterion + 1;
    scalar MaxResidualChange = 0;
    scalar MaxResidualPrev = 0;
    int IIcount=0;
        // Pressure-velocity SIMPLE corrector
        //for (int corr=0; corr<nCorr; corr++)
    while(maxResidual > convergenceCriterion)
        {
         // Momentum predictor
         eqnResidual = 1; 
         MaxResidualPrev=maxResidual;
         maxResidual = 0;
         IIcount++;

         tmp<fvVectorMatrix> UEqn
           (
                  fvm::ddt(U)
        + fvm::div(phi, U)
        + sgsModel->divDevBeff(U)
        );

         UEqn().relax();

        eqnResidual = solve
          (
           UEqn() == -fvc::grad(p)
           ).initialResidual();
        
        //  Retain maximum residual for global convergence check
        maxResidual = max(eqnResidual, maxResidual);

            p.boundaryField().updateCoeffs();
        
        volScalarField rUA = 1.0/UEqn().A();
        U = rUA*UEqn().H();
        UEqn.clear();
        
        phi = (fvc::interpolate(U) & mesh.Sf())
          + fvc::ddtPhiCorr(rUA, U, phi);
        
        adjustPhi(phi, U, p);
        
        // Non-orthogonal pressure corrector loop
        if (nCorr != 1) //Ulf aus pimple
          {
        p.storePrevIter();
          }
        //////--------------------------------------------------------------------------------------
        scalar StoppPressCorr =1;
        int CountPressCorr =0;
        scalar PrevPressEqnResidual=1.E20;
        //for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
        while ( (StoppPressCorr > convergenceCriterion) && (CountPressCorr < 10) )
          {
        CountPressCorr++;
        fvScalarMatrix pEqn
          (
           fvm::laplacian(rUA, p) == fvc::div(phi)
           );
        
        pEqn.setReference(pRefCell, pRefValue);

        eqnResidual = pEqn.solve().initialResidual();
        StoppPressCorr = eqnResidual;

        // Stop loop if residual increases, i.e. set loop counter to 10
        if (eqnResidual > PrevPressEqnResidual) CountPressCorr=10;
        PrevPressEqnResidual = eqnResidual;

        if (StoppPressCorr <= convergenceCriterion)
          {
            phi -= pEqn.flux();
            maxResidual = max(eqnResidual, maxResidual);
          }
        if (CountPressCorr == 10)
          {
            // Do not save maxResidual so that the outer, U-based maxResidual is not influenced if convergence is not reached 
            phi -= pEqn.flux();
          }
          }
        ////----------------------------------------------------------------------------------------
        
#           include "continuityErrs.H"

        // Explicitly relax pressure for momentum corrector
        p.relax();

        // Ulf, nach unten geschoben: DAS SOLLTE HINTER MOMENTUM CORRECTOR:  sgsModel->correct();
        
        // Momentum corrector
        U -= rUA*fvc::grad(p);
        U.correctBoundaryConditions();
        
        sgsModel->correct();
        // Info << "   ->  maxResidual(U,p) " << maxResidual << endl;


        // CONVERGENCE CHECK 

        MaxResidualChange=fabs(MaxResidualPrev - maxResidual);
        if (IIcount > 100 ) MaxResidualChange = MaxResidualChange * 100./IIcount;
        // If residual does not change anymore, stop loop 
        if ( (MaxResidualChange/maxResidual < 0.001) && (IIcount!=1) )
          {
        Info << "   ->  maxResidual(U,p) " << maxResidual << endl;
        Info << "   ->  MaxResidualChange (U,p)   " <<MaxResidualChange <<endl;
        Info << "   ->  MaxResidualChange/maxResidual   " <<MaxResidualChange/maxResidual <<endl;
        Info << "   ->  NOT reached convergence criterion " << convergenceCriterion << " on time step " << runTime.timeName() << " after " << IIcount <<" cycles " <<endl;
                Info << "   ->  maxResidual (U,p) set to 0 to end loop due to low MaxResidualChange" <<endl;
        maxResidual=0;
          }
        
        // FINALLY REACHED CONVERGENCE
        if (maxResidual < convergenceCriterion)
          {
        Info << "   ->  reached convergence criterion " << convergenceCriterion << " on time step " << runTime.timeName() << " after " << IIcount <<" cycles " <<endl;
          }

        }
    
        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return(0);
}


// ************************************************************************* //

[mvoss]

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

[djbungee]

Dear all,

as we were not able to present the energy-spectra (Engergy vs. wave-number) for the computations of decay of isotropic turbulence / box turbulence we made on the OpenFOAM-Stammtisch on Friday for different settings they are given here. The Fourier-transformation 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 under-relaxation (no under-relaxation: RELAX=1, and 50% under-relaxation: 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=1e-20). And we varied the time-step 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 time-step size. The magenta-colored curves are even
congruent
- With coupling of equations on a time-step 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 Rhie-and-Chow 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.

[charlie]

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:

Originally Posted by djbungee

- 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.

This is what I would expect.

Originally Posted by djbungee

- For 'small' time steps (blue curves), this is also true.

I wouldn't necessarily expect this. My feeling is that if the time step is fine enough, explicit and implicit time integration should become equivalent... but I'm really not sure about this.

Originally Posted by djbungee

- Without relaxation (green and magenta) the influence of not coupling the equations
decreases with decreasing time-step size. The magenta-colored curves are even
congruent

This is what I would expect: Relaxation "slows down" the progress that the solution makes between iterations, so the effect of "insufficient iterations" per time step is magnified with relaxation. That the magenta curves are identical is also what I would expect: with such a fine time step the system becomes equivalent to an explicit time integration and only one iteration per time step is sufficient (related to above point).

Originally Posted by djbungee

- With coupling of equations on a time-step the results get 'improved' (are closer
to each other), see the dashed lines.

I would also expect this and it is nice to see.

Originally Posted by djbungee

-> Coupling on a time step (inner loops, controlled by residuals) is required.

I agree entirely.

Originally Posted by djbungee

Another effect, but a totally different issue, is the dependence of the results on relaxation. Does this relate to the Rhie-and-Chow interpolation? Has anyone tried to evaluate that?

This throws me totally! I would expect (or maybe hope?) that the solution should never depend on the under-relaxation assuming the solution is converged on each time step. (Do I remember rightly that Fabian reported this behaviour too?). So, for my understanding of things, the dashed blue and magenta curves (for example) should be identical...

So, being sufficiently rattled, I decided to repeat your experiment with ELAN (our in-house code, also a co-located, FV, "transient SIMPLE" type). My setup was identical except for two things:

• Different initialisation field and molecular viscosity
• A peculiar bug meant that only one iteration per time step doesn't work. Could be considered a feature in light of this study . I therefore was lazy and ran the relevant cases with two iterations per time step rather than fix the bug.

The results are therefore not directly comparable, however similar trends should be observed if the solvers behave the same way. The spectra are attached, and I used your colour coding and nomenclature for ease of comparison.

spectraTimeStepRelax.png

So, what are my conclusions?

• I also see that insufficient inner iterations affects the results significantly. I also see this effect magnified by large time steps and under-relaxation, whereby the effect of the latter appears much stronger than that of the former.
• For a given time step and when sufficient iterations per time step are computed, the results are identical irrespective of the under-relaxation.

Hence, I can sleep well tonight rest assured that my intuition is not too far off the mark . However, it shows that the URF dependency shown by your OpenFOAM simulations isn't a necessary feature, and I think something must be wrong somewhere. I unfortunately have no idea where to start looking for an explanation .

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.

[djbungee]

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 . 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 do-while 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 Arbitrary-Euler-Lagrange () 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, LL-k-epsilon-based turbulence model for hybrid RANS/LES

[charlie]

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 under-relaxation 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.

[alberto]

Hello,

I have some comments, questions and a suggestion

Originally Posted by djbungee

Dear all,
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.

This is true only if the equations do not converge in that time step, meaning you are using an under-relaxation which is actually too aggressive and prevents the evolution of the solution in the number of prescribed iterations. Under-relaxation should not be applied at all in PISO-like algorithm for this reason, and in SIMPLE-like algorithm you must allow a sufficient number of iterations so that the solution actually is not affected by the under-relaxation.

- For 'small' time steps (blue curves), this is also true.

The same comment I wrote above applies.

- Without relaxation (green and magenta) the influence of not coupling the equations
decreases with decreasing time-step size. The magenta-colored curves are even
congruent

This is also known, and it is one of the basis of the PISO algorithm. You don't use under-relaxation but you choose an appropriate size of the time-step, typically relying on the CFL condition.

- With coupling of equations on a time-step the results get 'improved' (are closer
to each other), see the dashed lines.

That's because you allow the solution to reach its actual value at the end of the time step, through your sub-iterations.

-> Coupling on a time step (inner loops, controlled by residuals) is required.

This is surely true if you want to use large time steps, and I agree a convergence check is requires at each time step.

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 time-step adaption procedure, which has been already successfully used in MFIX (www.mfix.org), with a SIMPLE-like procedure.
The approach they use is very effective at speeding up the simulation, and the idea is the following:

• The SIMPLE method is used, solving the equations in a segregated fashion at each iteration
• At each SIMPLE corrector loop, the behaviour of the residuals is checked, based on what happened in the previous few iterations.
  • If the residuals stall or increase,
    • the SIMPLE loop is interrupted
    • the solution is reset to the value at the previous time-step
    • the time-step is reduced of a factor (0.8, for example)
    • a new SIMPLE loop with the new time-step is performed
    • The same procedure is applied to the new loop
  • If the solution converges fast, in a low number of iterations
    • The time step is increased of a factor (0.2) for the next time step. This does not represent a risk of divergence, since the checks on the residuals would capture the divergence and recover the previous solution before the divergence causes a crash in the solver

This procedure has various advantages:

• It ensures the convergence criteria are always satisfied
• It does not require a fixed number of sub-iterations per time-step, since the number is based on the convergence criteria
• It significantly reduces the computational time, by trying to use always (almost ) the maximum possible time step compatible with the convergence criteria
• It automatically detects divergence, and recovers the solution without interrupting the computation

There is a disadvantage too: computing statistics becomes trickier

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

Best,

[djbungee]

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 .

Originally Posted by alberto

This is true only if the equations do not converge in that time step,
meaning you are using an under-relaxation which is actually too aggressive and prevents the evolution of the solution in the number of prescribed iterations.
Under-relaxation should not be applied at all in PISO-like algorithm for this reason, and in SIMPLE-like algorithm you must allow a sufficient number of iterations so that the solution actually is not affected by the under-relaxation.

That is what I have written, but the concept does not go far or general enough. You have to ensure a coupled convergence. And, and as mentioned by us, the number of iterations was NOT fixed but only after convergence (initial residual 'small enough') the solution is proceeded by one time step. This is independent of SIMPLE, SIMPLEX, SIMPLER or PISO or whatever.

Originally Posted by alberto

My question is: where the solutions in the different cases actually converged? If so, what criterion was adopted to define the convergence?

Please take a look at the source code where you can see the rather simple convergence check. It was the first development step and not very elegant, however, effective and that is one of the features, that I marked as one that can easily be improved.

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 pressure-based procedure. We rather think that this comes from the way pressure-velocity decoupling on non-staggered grids is avoided by employing a Rhie-and-Chow interpolation.

Best regards, Ulf.

[alberto]

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

Originally Posted by djbungee

You have to ensure a coupled convergence.

My question about convergence was on this point. Thanks for the clarification.

And, and as mentioned by us, the number of iterations was NOT fixed but only after convergence (initial residual 'small enough') the solution is proceeded by one time step. This is independent of SIMPLE, SIMPLEX, SIMPLER or PISO or whatever.

Well, in the pure PISO algorithm you do not have the convergence check at all inside the time step, and you're simply driven by your Courant numbers, so it makes some difference. Of course this does not ensure anything about coupled convergence.

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!

Yes, you're bringing residuals at extremely low levels (10e-20), which is much lower than what is usually done.

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.

OK.

This should not depend on the way you implement your pressure-based procedure. We rather think that this comes from the way pressure-velocity decoupling on non-staggered grids is avoided by employing a Rhie-and-Chow interpolation.

Yes, it surely can be. The standard Rhie-Chow interpolation is sensitive to under-relaxation, and might give checkerboarding when small time steps or small under-relaxation factors are used.

To stabilize the interpolation in those cases, a modified Rhie-Chow interpolation is used, including the effect of the body-force 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

Best,

[alberto]

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: 221-239, 2003

In summary

• They consider the solution proposed in S. Majumdar, Role of underrelaxation in momentum interpolation for calculation of flow with nonstaggered grids, Numerical Heat Transfer, 13(1): 125-132, 1998. This solutions consists in replacing the fluxes interpolated from the velocities at the previous iteration with the fluxes already obtained from the previous pressure-correction step.
• Shen at al show this is not sufficient for the SIMPLEC method, which becomes not convergent when their correction is applied, since the p-v coupling becomes inconsistent.
• They propose a solution, function of a parameter beta, used to weight the contribution of the previous step solution in the interpolation formula, claiming that, if beta is fixed, the solution is not dependent on the value of URF_u (URF_p does not matter, since the optimal value for SIMPLEC is 1).
• In the paper they show that, however, the improved algorithm becomes, at least slightly, grid-dependent. In addition beta influences the number of iterations to obtain convergence (1500 with beta = 0.03, 700 with beta = 0.04, 900 with beta 0.07...with a non-monotonic behaviour).

Best,

[djbungee]

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.

[alberto]

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.

You're welcome. If you can keep me posted on this when you have some update or question (take all the time you need, I don't want to sound pushy ), it would be great.

Best,

[charlie]

Also from me thanks to Alberto for this very comprehensive advice - this will surely help us to overcome this problem!

Best regards,

Charlie.

[a.carnarius]

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 non-linear 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

[djbungee]

Dear Angelo, dear all,

just one comment on your statement:

Originally Posted by a.carnarius

this is a very interesting topic as I was really surprised to see, that OpenFOAM seems to have problems with such fundamental numerical schemes.

Actually, this is not a problem of OpenFOAM (and I am sorry if I should be responsible for causing this impression ). There might still be solvers that already have remedies for the dependence on relaxation available (cf. Alberto's post) and pimpleFoam already has inner loops implemented (but with a fixed number of steps). So it is our task to get this knowledge / information / source code together and make it available to the community.

Thanks to all of you for this very fruitful discussion.

Best regards, Ulf.

[alberto]

Originally Posted by a.carnarius

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.

The dependence of the solution on the relaxation factor when Rhie-Chow interpolation is used is not something alarming anyway, since it is in most of the cases of order of magnitude smaller than the discretization error, according to Ferziger and Peric (see their book, pp 198 and following):

• In a lid-driven cavity case, they compare the difference between the solutions, after residuals were reduced of five orders of magnitude, with URF_u = 0.9 and URF_0.5, obtaining
  • 0.5% difference on a 8 x 8 CV grid - Dicretization error 53%
  • 0.002% difference on a 128 x 128 CV grid - Discretization error 0.4%
• In a buoyancy-driven case, they repeated the test, obtaining
  • 0.23% on a 8x8 CV grid
  • 0.0008% on a 128x128 CV grid

In all the cases, the error is significantly lower than the discretization error. Probably this comparison should be considered in the calculations discussed in this thread.

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 grid-dependent.

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.

Actually the fact the error goes to zero with the time step is essential. If it doesn't, you do not have time-convergence, and your scheme should not be used.

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.

Well, for time-resolved simulations like LES should be to have it done right, the current PISO implementation in OpenFOAM is OK. The CFL condition can be often less restrictive than the characteristic lenght scale of your flow.

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.

I will post some results of test calculations on this topics today or on monday for further discussion.

That would be great!
What test case are you using?

Best,
Alberto

[djbungee]

Dear all,

what you mentioned, Alberto, is right:

Originally Posted by alberto

In all the cases, the error is significantly lower than the discretization error. Probably this comparison should be considered in the calculations discussed in this thread.

Moreover, please be aware that the dependence of the results on relaxation due to employing a Rhie-and-Chow interpolation was only an assumption.

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:

• only the relaxation for velocity causes the difference for with and without relaxation as initially described.
• turbulence and (explicit) pressure relaxation do not influence the results in the presented implementation.
• doing a final loop after convergence on a time step without relaxation of momentum improves the situation (only) somewhat.

Therefore, the recommendation for this procedure and also, e.g., for pimpleFoam can only be not to apply relaxation on the momentum equation.

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.

[alberto]

Hello and thank you for the update

Originally Posted by djbungee

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.

Good. I expected that, since in your case you do not have any reason why this should make a difference. In this sense the code behaved correctly.

• only the relaxation for velocity causes the difference for with and without relaxation as initially described.
• turbulence and (explicit) pressure relaxation do not influence the results in the presented implementation.
• doing a final loop after convergence on a time step without relaxation of momentum improves the situation (only) somewhat.

This seems consistent with what observed in some of the references we talked about before.

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?

I did not try anything yet. What we could do is to start with a very simple case without any submodel. Probably the lid-driven cavity test case could be a good one, with a laminar flow condition. It would also allow to study the problem both with steady and unsteady solvers, and the computational requirement is low.

What do you think?

Best,

[niklas]

Moved the thread to this forum instead.
It seems more appropriate.

[djbungee]

Dear Niklas,

Originally Posted by niklas

Moved the thread to this forum instead.
It seems more appropriate.

thanks. However, there might be a bug and we might further develop the code... so it fits in either forum.

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.