[sheaker]

Hello dear Foamers. I'm using openFoam ext 1.6 to do some internal combustion engine simulations with sonicTurbDyMEngineFoam . To keep Courant number below 1, timestep must be set to 1.5e^-7s. For a simple geometry, it takes 90000sec to simulate 720 CA degree. Geometry example:
https://www.youtube.com/watch?v=THJvHIk0n0Y
https://www.youtube.com/watch?v=jVls41R8j_c

I decided to add PIMPLE loop instead of SIMPLE for sonicTurbDyMEngineFoam to increase maximum Courant Number and to gain some experience in modifying solvers. It shouldn’t be very complicated for the beginners.
First of all I have changed sonicTurbDyMEngineFoam.C from

  [...]
#       include "rhoEqn.H"

        // --- SIMPLE loop
        for (int corr=1; corr<=nCorr; corr++)
        {
#           include "UEqn.H"

#           include "hEqn.H"

#           include "pEqn.H"
        }

        turbulence->correct();
  [...]

To

  [...]
//        PIMPLE LOOP
        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
            {
#           include "rhoEqn.H"
#           include "UEqn.H"
#           include "hEqn.H"
            // --- PISO loop
            for (int corr=1; corr<=nCorr; corr++)
                {
#               include "pEqn.H"
                }

            turbulence->correct();
            }        
    [...]

There is nothing new here – only additional loops.

Then I have changed from

  #include <readPISOcontrols.H>

to

  #include <readPIMPLEcontrols.H>

in two files: sonicTurbDyMEngineFoam.C and readControls.H to read PIMPLE instead of PISO sets.


Is that solver properly prepared? I have compiled it (with different name) but with warnings about no used variables like nNonOrthoCorr, etc. . I think those warnings comes from nowhere because those warnings appears also when I compile original solver.

In case from sonicTurbDyMEngineFoam I have only modify fvSolutions from:

  PISO                                                     
  {
      // SIMPLE correctors
      nOuterCorrectors   2;
   
      // PISO correctors
      nCorrectors        3;
   
      momentumPredictor yes;
      transonic         no;
      nNonOrthogonalCorrectors 0;
  }
   
  relaxationFactors
  {
      U          0.7;
      p          0.25;
      h          0.5;
      k          0.7;
      epsilon    0.7;
  }

To

  PIMPLE
  {
      nNonOrthogonalCorrectors 0;
      nCorrectors          3;
      nOuterCorrectors    10;
   
      residualControl
      {
          U
          {
                  tolerance  1e-8;
                  relTol      0.01;
          }
          p
          {
                  tolerance  1e-7;
                  relTol      0;
          }
                  k
                  {
                  tolerance  1e-8;
                  relTol      0.01;
                  }
                  epsilon
                  {
                  tolerance  1e-8;
                  relTol      0.01;
                  }
                  h
                  {
                  tolerance  1e-8;
                  relTol      0.01;
                  }
       }
  }
   
  relaxationFactors
  {
      fields
      {
          p      0.25;
          pFinal   1;
      }
      equations
      {
          U    0.7;
                  h             0.5;
                  k             0.7;
                  epsilon 0.7;
          "(U|h|k|epsilon)Final"   1;
      }
  }

Is that case properly prepared? I can run it but after 10CA, when maximum velocity occure, the following error appears:

[...]
    From function specieThermo<thermo>::T(scalar f, scalar T0, scalar (specieThermo<thermo>::*F)(const scalar) const, scalar (specieThermo<thermo>::*dFdT)(const scalar) const) const
    in file /home/sheaker/OpenFOAM/OpenFOAM-1.6-ext/src/thermophysicalModels/specie/lnInclude/specieThermoI.H at line 73
    Maximum number of iterations exceeded.  Rescue by HJ

    From function specieThermo<thermo>::T(scalar f, scalar T0, scalar (specieThermo<thermo>::*F)(const scalar) const, scalar (specieThermo<thermo>::*dFdT)(const scalar) const) const
    in file /home/sheaker/OpenFOAM/OpenFOAM-1.6-ext/src/thermophysicalModels/specie/lnInclude/specieThermoI.H at line 73
    Maximum number of iterations exceeded.  Rescue by HJ

    From function specieThermo<thermo>::T(scalar f, scalar T0, scalar (specieThermo<thermo>::*F)(const scalar) const, scalar (specieThermo<thermo>::*dFdT)(const scalar) const) const
    in file /home/sheaker/OpenFOAM/OpenFOAM-1.6-ext/src/thermophysicalModels/specie/lnInclude/specieThermoI.H at line 73
    Maximum number of iterations exceeded.  Rescue by HJ
BiCGStab:  Solving for epsilon, Initial residual = 3.61942e-05, Final residual = 6.94297e-13, No Iterations 1
bounding epsilon, min: 1.68981e-20 max: 2.56723e+07 average: 70394.6
BiCGStab:  Solving for k, Initial residual = 0.0717658, Final residual = 1.22763e-12, No Iterations 1
bounding k, min: -3.38742 max: 544.101 average: 5.35043
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
Floating point exception (core dumped)
sheaker@sheakerLAPTOP:~/sonicTDMEFoamPIMPLEtest$

I had similar errors when I set Courant number to 2 and ran simulation with SIMPLE loop. So it looks like adding PIMPLE loop decrease stability... Or rather I have done something wrong. I believe some of You are experienced in coding and could give me a hint, maybe.

I wish You all good time. Thank You for reading.
sheaker

[piu58]

Without analysing all the material:

The need of a very small time step is often caused by a unbalanced mesh: The Co number is very different in different parts of the mesh. The solution log gives max and min Co number so you may look at this. A relation of 10 may be acceptable, much higher quotients needs to work at the mesh.

[sheaker]

Thank You Uwe for Your replay. The mesh is pretty simple, the simulation is not steady-state and in different timesteps different velocities occur in some regions so it is rather impossible to keep maxCo / minCo = 10 without auto remeshing.
After all I would really like to understand how to successfully add the PIMPLE algorithm because I think there is a problem with my implementation or case setup.

Have a nice day.
sheaker