## Contents

## 1 New reference

If you want to check out the latest work from Tobias Holzmann about the usage of the PIMPLE algorithm, you can checkout his book »Mathematics, Numerics, Derivations and OpenFOAM(R)«. Here, you will find an own chapter about the whole stuff. In addition reference [21] on www.holzmann-cfd.de explain the pressure-velocity coupling very well.

## 2 Storyboard

First I wanted to write a nice story like »Ignaz« but that would only copy the good attention of the author and of course I am not a good story-teller. Therefore I will go directly to the topic and will tell you how the PIMPLE algorithm is working in OpenFOAM, how to set up the keywords and which abilities you have to control the algorithm. The test case it the pitzDaily case in the pimpleFoam folder of incompressible fluids. Only the Courant number was decreased to 3.

## 3 Case A

In the first case we are using the tutorial case and change the fvSolution file. In the pimple control directory we will remove all entries to get only a dummy body.

PIMPLE { }

After starting the simulation you can see that the solver is still running but after a while it will crash. The following questions occur at that case:

- why is the solver working without complaining that keywords are missing
- in which mode is the solver running when there are no keywords

**Question 1**:

- If you don't insert any keywords into the PIMPLE control dictionary, the solver is still running using the default settings. After calling the pimpleControl the solver initialize the necessary keywords with default values. These can be checked out in the
*pimpleControl.C*file:

121 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // 122 123 Foam::pimpleControl::pimpleControl(fvMesh& mesh) 124 : 125 solutionControl(mesh, "PIMPLE"), 126 nCorrPIMPLE_(0), 127 nCorrPISO_(0), 128 corrPISO_(0), 129 turbOnFinalIterOnly_(true), 130 converged_(false) 131 { 132 read();

- As you can see there are several keywords that can be used if you are using the PIMPLE algorithm. These keywords are set »by default« to the values you can see at the code above. After initialize the values, the constructor will call the »read()« function which will check out the user defined values and overwrite the above values with user values or other values which are necessary that the solver is working. If there are no values in the PIMPLE control available the solver is using the default one. The functions is given below:

39 void Foam::pimpleControl::read() 40 { 41 solutionControl::read(false); 42 43 // Read solution controls 44 const dictionary& pimpleDict = dict(); 45 nCorrPIMPLE_ = pimpleDict.lookupOrDefault<label> ("nOuterCorrectors", 1); 46 nCorrPISO_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1); 47 turbOnFinalIterOnly_ = 48 pimpleDict.lookupOrDefault<Switch>("turbOnFinalIterOnly", true); 49 }

- As you can see here, this function is looking in the
**fvSolutions**file for some entries in the PIMPLE control structure:

nOuterCorrctors (nCorrPIMPLE) is set by default to 1 nCorrectors (nCorrPISO) is set by default to 1 nNonOrthogonalCorrectors (corrPISO) is set by default to 0 (not shown here) turbOnFinalIterOnly is set by default to true Residual control is set to false

**Question 2**:

- If you are running the PIMPLE solver without any keywords it is operating in PISO mode and therefore you need very small time steps. Thats the reason why the simulation will blow up with these settings (Co = 3). That the solver is operating in PISO mode can be shown in the source code too, and will report it on the screen:

134 if (nCorrPIMPLE_ > 1) 135 { . . . 156 } 157 else 158 { 159 Info<< nl << algorithmName_ << ": Operating solver in PISO mode" << nl 160 << endl; 161 }

**Residuals**if you are using a empty PIMPLE directory

- As you can see the residuals are not converging and the solution blow up. The reason is due to the high Co number (Co = 3) and the fact that without keywords the PIMPLE algorithm is working in PISO mode. To solve the problem you have to reduce the Co number but that is not is not the aim to introduce merged PISO - SIMPLE algorithm.

**Legend**:

- Iteration at the bottom means one PISO iteration; that means in each iteration you go on with the time. In other words every iteration solves a new time step.

## 4 Case B

In the second case we are using the standard entries in the **fvSolutions** file used in the tutorial.

PIMPLE { nCorrectors 2; }

Now we correct the pressure field twice and that means that we calculate the following equation twice:

In that case the pressure - momentum coupling is stable enough to work and the solver will not blow up.

**Advantage**- While you are using nCorrectors with a number higher than 1 the PIMPLE algorithm is still working in PISO mode with the advantage that you calculate your pressure twice, this is also available in the normal PISO solver but is still not the aim of the PIMPLE solver.

- The residual of the pressure is much better and the solution is stable in that moment - not tested what happen if you solve more seconds.

**Legend:**- For the momentum one iteration means »one time step«, in other words that means in that case we were solving round about 65 time steps
- For the pressure we solve the equation given above two times and therefore it is not anymore possible to say that one iteration is one time step. In that case two iterations are necessary to go on in time.

## 5 Case C

Now we will use the MERGED PISO - SIMPLE algorithm and therefore we are using the advantage of the SIMPLE mode idea (relaxation). What does that mean? In fact it is really easy. If you are using the PISO mode you are strict limited to the time step (Co < 1). That mean that it is very expensive to solve a real time problem in transient manner for long time especially in very complex geometries. SIMPLE is developed for steady-state condition to reach that state very fast but (normally) does not contain time derivation and therefore no time information. To provide a stable simulation you have to introduce relaxation factors that limit the values of the new values.

In other words, the new value is a combination of the old value and the new one which is calculated. This make the simulation stable but not time conservative. It allows you to reach the steady-state condition very fast; time information is not needed.

**Miss understanding**

- A lot of people miss understand the function of the PIMPLE algorithm. The problem in the case we are using here is, that this case is very simple. Imagine a very complex simulation in which it is very straight forward to simulate small time steps. If you want to have bigger time steps you have to increase your Co number which leads you to unstable convergence. At that moment you are in a discrepancy because you want go on faster in time step to reduce computational effort but it is not possible because you are limited in time step due to the Co number.

### 5.1 Case C.1

Now we will use the following **fvSolutions** file:

PIMPLE { nCorrectors 1; nOuterCorrectors 2; }

Compared to Case B there is not a big difference. The only difference is that we calculate the pressure - momentum coupling in one time step twice:

- Calculate momentum
- Calculate pressure
- Recalculate momentum from new pressure
- Recalculate pressure from new momentum

- 1.

- 2.

- correct U with

- now you use the and for the second loop

- 3.

- 4.

- correct U with

During this calculation it has not a big advantage and is like Case A calculated twice and then move on in time step (with the only difference, that the turbulence is calculated at the end. That means in PIMPLE Iteration 2 due to the fact that *turbOnFinalIterOnly = true*). You can assume that the calculation will blow up which is shown in the residual plot.

**Residuals**if you are using the PIMPLE algorithm with the settings above- As you can see it is much worse compared to Case 1. At that moment there is none advantage of the pimple mode that we are using but we are using it in a wrong way.

### 5.2 Case C.2

Now we want to generate a complete PIMPLE case. The **fvSolutions** file looks like that:

PIMPLE { nNonOrthogonalCorrectors 0; nCorrectors 1; nOuterCorrectors 2; // Set to 50 relaxationFactors { fields { p 0.3; pFinal 1; } equations { "U|k|epsilon" 0.3; "(U|k|epsilon)Final" 1; } }

We introduced the relaxation factors to obtain a smooth convergence. The problem now is that we use the relaxation factors in the first pimple loop but after we reach the last one (here iteration 2 because *nOuterCorrectors = 2*) the relaxation factor will be 1 again to ensure time consistence but this makes no sense in the following case. Therefore we will change the nOuterCorrectors to 50. Now we calculate 50 times the pressure - momentum coupling. Therefore 49 times with relaxation (like SIMPLE) and the last time without relaxation. For cases that are very unstable this is perfect to get a smooth convergence. In other words we are using the idea of SIMPLE mode (under relaxation) between the time steps to converge the solution at these time step to a certain limit. The limit is set in a way to be sure the solution is accurate in time and space.

**Residuals**if you are using the right PIMPLE algorithm with relaxation- You can see that the solution is stable and the convergence is very good. Further more you can see the pimple loops; here from zero to fifty iteration, 100 to 150 and so on. After you finished the pimple loop we go on in time and do the same thing again. If you introduce the
*nCorrectors > 1*you will calculate the pressure field n times (more hints to that later).

**Do not forget**

- At the moment we have 50 Iterations per time step and therefore the computational effort could be bigger as before. To avoid that you should always add residual controls to stop the pimple loop after reaching some special criterion » Case C.3

### 5.3 Case C.3

At last, we introduce the residual control to exit the pimple loop after we reached a good solution during one time step. The **fvSolutions** file looks like that:

PIMPLE { nNonOrthogonalCorrectors 0; nCorrectors 1; nOuterCorrectors 50; residualControl { U { tolerance 1e-5; relTol 0; } p { tolerance 5e-4; relTol 0; } } } relaxationFactors { fields { p 0.3; pFinal 1; } equations { "U|k|epsilon" 0.3; "(U|k|epsilon)Final" 1; } }

**Residuals**if you are using the right PIMPLE algorithm with relaxation and residualControl- As you can see now, the time iteration is faster because after the pimple loop will reach the residual control it will go to the next time step.

smoothSolver: Solving for Ux, Initial residual = 2.07804e-06, Final residual = 2.07804e-06, No Iterations 0 smoothSolver: Solving for Uy, Initial residual = 3.42321e-06, Final residual = 3.42321e-06, No Iterations 0 GAMG: Solving for p, Initial residual = 0.000407312, Final residual = 3.14679e-08, No Iterations 8 time step continuity errors : sum local = 5.66964e-12, global = 3.09312e-13, cumulative = -8.61807e-06 smoothSolver: Solving for epsilon, Initial residual = 0.00512886, Final residual = 2.08382e-06, No Iterations 2 smoothSolver: Solving for k, Initial residual = 0.0660137, Final residual = 1.08191e-06, No Iterations 3 '''PIMPLE: converged in 28 iterations''' ExecutionTime = 13.84 s ClockTime = 14 s

- Although we set 50 iterations for the pimple loop it will exit the loop after the residual control settings are true. That saves computational effort.

## 6 Some hints

- the pitzDaily tutorial is to easy to show the advantage of PIMPLE in a good manner

- PIMPLE is very effective in big cases to increase the Co number

- PIMPLE speed up your simulation and could reach a stable solution due to under relaxation

- For high turbulent flow you should add turbOnFinalIterOnly = false to your PIMPLE control

- Set nOuterCorrectors up to a very high number (~ 50 to 200) and control your pimple loop with the residual control

- Set nCorrectors to a small number (1 to 3; you have to play around)

- If you set nCorrectors higher it could be possible to leave the pimple loop earlier and prevent extra calculation but could also do the opposite

- nNonOrthogonalCorr means that you correct the pressure field more often with the new calculated one:

- till nNonOrtogonalCorr is reached

- The optimum relaxation factors depend on the simulation case you are running
- The converged solution is in depended of the relaxation factors which are applied
- Due to the fact that you have bigger time steps it could happen that some physics are not resolved (extreme fast pressure drop etc.)
- At the beginning of an transient simulation it could happen that after 100 pimple loops the solution is still not converged but this will change during the simulation time

- simple heat transfer case has at the beginning 200 pimple loops and is still not converging
- after the tenth time step it converges and the pimple loops decrease significantly till it get a steady state pimple iteration of round about 3 - 5

**Please feel free to change this site if I got something wrong**