[shariq12]

Greetings of the time,
I am trying to code a standard SIMPLE algorithm within OpenFOAM. Here's the logic,

Start by computing U* (Initialise); I do this outside the main while loop as:

fvVectorMatrix UEqn
     (
        // div(U * phi) - div (nu * grad(phi)) = -div(p) 

        fvm::div(phi,U) - fvm::laplacian(nu,U) == -fvc::grad(p)
    );
    UEqn.solve();

With this, I then enter the while loop, compute the scalar matrix of coefficient A, and interpolate its inverse to the face of the volume:

volScalarField A = UEqn.A(); 
surfaceScalarField A_inv = fvc::interpolate(1/A);

Next, ignoring the correction from neighbouring cell centres, I define corrections as:
U' = -grad(p)' / A


Now we know, U = U* + U', i.e.
div(U) = 0 // Continuity equation
div(U* + U') = 0
div(grad(p)' / A) = div(U*) // This equation drives the pressure to correct
//velocities

I do this using,

fvScalarMatrix pEqn
(
    fvm::laplacian(A_inv,p) == fvc::div(phi)
);
pEqn.setReference(pRefCell,pRefValue); // Provided by the user
pEqn.solve();

Finally, I correct the velocities and fluxes using

U = U - 1/A * (fvc::grad(p));    // U = U* + U'
phi = fvc::interpolate(U) & mesh.Sf();

Is this implementation correct? When testing it on a simple 2D case, I can get convergence with appropriate velocity values. However, my pressure seems to be off and has oscillations.

Can someone please point out if there's a flaw in my logic?

Thank You

[dlahaye]

Correct. The Rhie-Chow interpolation is there to avoid pressure oscillations expected is cell-centered discretizations.

[shariq12]

Hey,

Thanks for the input on the odd-odd and even-even checkerboard pattern that is well-known when both the velocities and pressures are stored in a collocated fashion.
However, I faced issues (mismatch between OF SIMPLE algorithm) with the previous implementation primarily because I was not updating the pressure field correctly.
I managed to fix this. Here's how I implemented the standard SIMPLE algorithm (This is not the same as OF implementation, which is similar to SIMPLEC with a few differences):

1. Start by declaring the appropriate declaration in the includeFields.H; this includes the volVectorField U and volScalarField p.
   
   
2. Initialise the guess values (^*) with the initial conditions (Note: we want to compute the corrections denoted by ^')
   
   

   // Initilise U* , p* , phi* and p' NOTE: U = U* + U' and so on.
   volVectorField U_star(U);
   volScalarField p_star(p);
   surfaceScalarField phi_star(phi);
   volScalarField p_cor(p);

3. Enter the time loop and solve the momentum equations with the guess values, (Note: we are solving for U_star; this will not satisfy the continuity equation).
   
   

   fvVectorMatrix UEqn
   (   
   // div(U * phi) - div (nu * grad(phi)) = -div(p) , Note in this case phi is u,v,w scalar components of velocities.
      fvm::div(phi_star,U_star) - fvm::laplacian(nu,U_star) == 
      fvc::grad(p_star) 
   );
   // Solve for U*
   UEqn.solve();

4. Obtain the scalar matrix of the coefficient and interpolate it to the faces.
   
   

   // a_p * U_p' = sum(i->n){a_i * U_i'} + sources - div(p')
   // U_p' = sum(i->n){a_i * U_i'} / a_p + sources / a_p - div(p') / a_p
   // For Std SIMPLE ignore contri from neighbours
   // i.e. U_p' = div(p') / a_p.
   volScalarField A = UEqn.A();
   surfaceScalarField A_inv = fvc::interpolate(1/A);
   phi_star = fvc::interpolate(U_star) & mesh.Sf();

5. Assemble and solve the pressure correction equation (i.e. p').
   
   

   // Now div(U) = 0
   // div(U* + U') = 0
   // div(U*) + div(U') = 0
   // div(U') = -div(U*)
   // div(-1/a_p * div(p')) = -div(U*)
   fvScalarMatrix pEqn
   (
      fvm::laplacian(A_inv,p_cor) == fvc::div(phi_star)
   );
   
   // Solve for p'
   pEqn.setReference(pRefCell,pRefValue);
   pEqn.solve();

6. Correct the fields and make sure to apply relaxations.
   
   

   p = p_star + (alpha * p_cor);
   
   // Correct U as U = U* + U' , with U' = div(-1/a_p * div(p')) also apply relaxation.
   // U' is computed using the method highlighted in step 4.
   U = U_star - 1/A * (fvc::grad(p_cor));
   U = alpha * U + (1 - alpha) * U_old;
   
   // Correct phi now
   phi = fvc::interpolate(U) & mesh.Sf();

7. Update p* and U* as p and U respectively, ensuring to correct the BCs and saving the present U to U_old (this is to apply relaxation in the next it).
   
   

   // Set p* = p ; U* = U ; phi* = phi and repeat
   p_star = p;
   U_star = U;
   phi_star = phi;
   
   
   // Correct BCs
   U.correctBoundaryConditions();
   p.correctBoundaryConditions();
   
   // Save for relaxation
   U_old = U;
   
   runTime.write();

   
   I hope this helps someone. Cheers!!

[Tobermory]

Is that any different to what simpleFoam does? (I didn't take the time to read carefully through each of your steps)

[shariq12]

Originally Posted by Tobermory

Is that any different to what simpleFoam does? (I didn't take the time to read carefully through each of your steps)

Thank you for your interest,

Indeed, the standard SIMPLE algorithm implemented in simpleFoam, is quite different from the standard SIMPLE algorithm developed by Patankar found here (https://en.wikipedia.org/wiki/SIMPLE_algorithm).

Let be briefly highlight the difference

Initial step is same for both in which we start by solving for U* (guess velocity); this however will not satisfy the continuity. This is where things start to differ, in the classic SIMPLE algorithm we explicitly define a correction U' to correct the velocity field i.e. (U = U* + U'); and this corrected velocity(U) is what we substitute in continuity to derive the correction U'; mathematically:

Notice how because of the presence of U* we get a sink component that is known from the initial guess, intiuitively this could be pictured as continuity error in each control volume i.e.

The idea now is to come up with a formulation of U' that can be substituted in the equation above, the correction (U') can very simply be shown to be

Now based on the requirement we can ignore the contribution from H' (or neighbouring cells) this leaves with just the pressure correction p' which can be substituted back in the velocity correction equation and by doing this we can solve for p' (Notice that p' is driven by the presence of ); once this is done we can correct the pressure as p = p* + p' and similary velocity can be corrected from the equation above.

In the simpleFoam implementation, we do not define an explicit correction for p, i.e., we solve for the pressure equation directly. This is similar to the SIMPLEC algorithm (can be found here https://en.wikipedia.org/wiki/SIMPLEC_algorithm); mathematically,

Solving this gives the corrected pressure p; the next step is to correct the velocity field, and this is where simpleFoam's implementation differs from the SIMPLEC algorithm. The pressure is used to update the velocity field as

This process is then repeated in a loop. Clearly, the pressure here is driven by H. Additionally, notice that we have not explicitly included a velocity predictor step of the form (U = U* + U') which is typically seen in SIMPLE. In fact, we are solving for the corrected pressure field first (using the velocities field which would not satisfy the continuity) and then correcting the velocity fields using the same pressure to obtain corrected velocity. This is computationally cheaper when compared to SIMPLEC, where the corrected pressure field is used to solve the momentum equation within the same loop before applying the standard pressure correction to update the velocity. And it is more robust when compared to the SIMPLE algorithm because we are directly solving for the corrected pressure field. Additionally, it's is quite efficient as we do not need to store p*.

Let me know if this helps; I can share my implementation for all three algorithms if required. Cheerrss

[Tobermory]

Thanks for your explanation - I was wondering why you were interested in reimplementing the SIMPLE scheme when there is already a lot of info in the public domain on OpenFOAM's implementation (eg https://openfoamwiki.net/index.php/SimpleFoam).

However - I understand now - the first time I read through simpleFoam, trying to compare against Ferziger & Peric (which is well worth reading as well as Patankar), I was also a bit puzzled at the differences. But I realised that they make sense when you realise that you want the same code to be able to run SIMPLE and SIMPLEC. Your previous post should help others in the future - good job.