21_simpleSol

Important

Visit https://aerosand.cc for the latest updates.

0. Preface

In the previous discussion, we examined the SIMPLE algorithm in detail and briefly looked at some code snippets of the SIMPLE algorithm in OpenFOAM.

Here, we will implement the code based on our discussion of the algorithm and understand the usage of some code.

This section primarily discusses:

Implementation of the SIMPLE algorithm
Computing the cavity test case

1. Governing Equations

The governing equations are as follows:

Continuity equation (mass equation):

\nabla\cdot U = 0

Momentum equation:

\nabla \cdot (UU) = \nabla\cdot(\nu\nabla U)-\nabla p

The following assumptions are still applied:

Steady state
Viscous term is simplified
Gravity is neglected
Density has been accounted for

2. Project Preparation

Run the following commands in the terminal to create the project:

terminal

ofsp
foamNewApp ofsp_21_simpleSol
cd ofsp_21_simpleSol
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity debug_case
code .

2.1. Documentation File

Provide a documentation file for the project:

README.md

## About

A solver that simply reproduces the SIMPLE algorithm.

## Bio

aerosand @ aerosand

## Caution

Pay attention to the OpenFOAM version.

## Deploy

Ensure the solver files are complete and include the debug_case folder

Enter the following commands in the terminal:

./Allclean
./Allrun


## Event

@20260312

- Added scripts #done

2.2. Script Files

We will streamline the previous scripts.

Run the following command in the terminal to create the scripts:

terminal

code {Allclean,Allrun}

The cleaning script Allclean is as follows:

Allclean

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#!/bin/sh
cd "${0%/*}" || exit                                # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
#------------------------------------------------------------------------------

cd debug_case

# Function: initialize the 0 folder
init_0() {
    if [ -d "0.orig" ] && [ -d "0" ]; then
        echo "Both 0 and 0.orig exist. Removing 0 and restoring from 0.orig..."
        rm -rf 0
        cp -a 0.orig 0
    elif [ -d "0" ] && [ ! -d "0.orig" ]; then
        echo "Backing up 0 to 0.orig..."
        cp -a 0 0.orig
    elif [ -d "0.orig" ] && [ ! -d "0" ]; then
        echo "Restoring 0 from 0.orig...".
        cp -a 0.orig 0
    else
        echo "Warning: Neither 0 nor 0.orig exists. Skipping."
    fi
    echo ">>>>>>>>>>>>> Initial condition done."
}

# Ensure initial conditions
init_0

# Clean the case
. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
cleanCase0

The run script Allrun is as follows:

Allrun

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
#!/bin/sh
cd "${0%/*}" || exit                                # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
#------------------------------------------------------------------------------

cd debug_case

restore0Dir

runApplication blockMesh

runApplication $(getApplication)

Do not worry about the scripts; as the discussion deepens, we will continuously expand the script content and try different writing styles.

Note

Use chmod +x {Allclean,Allrun} to grant execution permissions to the scripts.

These scripts can continue to be used in the future unless otherwise noted.

3. Project Implementation

We will implement the SIMPLE algorithm with the simplest code possible.

3.1. Main Source Code

Considering the discussion in 20_simple, implement the main framework of the SIMPLE algorithm in the main source code:

ofsp_21_simpleSol.C

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
#include "fvCFD.H" // Refer to 07_firstApp, will not be repeated

#include "simpleControl.H" // Header file containing the SIMPLE algorithm control
// This header file belongs to the finiteVolume library; no additional linking is needed in Make


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

int main(int argc, char *argv[])
{
    #include "setRootCase.H" // Refer to 13_commandLine, will not be repeated
    #include "createTime.H" // Refer to 14_time, will not be repeated

    #include "createMesh.H" // Refer to 15_mesh, will not be repeated
    #include "createControl.H" // Create algorithm control; can be automatically created based on usage
    #include "createFields.H" // User-provided field inclusion


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

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

    while (simple.loop()) // SIMPLE algorithm loop
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // --- Pressure-velocity SIMPLE corrector
        {
            #include "UEqn.H" // Momentum predictor
            #include "pEqn.H" // Pressure and momentum correction
        }

        runTime.write();

        runTime.printExecutionTime(Info);
    }

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

    return 0;
}

It can be seen that the main algorithm framework is the same as discussed in the previous section.

3.2. Field Inclusion

The field inclusion createFields.H is as follows:

createFields.H

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
/*
    # Basic fields
    # transportProperties
    # Pressure reference
*/

// # Basic fields

Info<< "Reading field p\n" << endl; // Pressure field
volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading field U\n" << endl; // Velocity field
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

#include "createPhi.H" // Mass flux field

// # transportProperties

Info<< "Reading transportProperties\n" << endl;

IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

dimensionedScalar nu
(
    "nu",
    dimViscosity,
    transportProperties
);


// # Pressure reference

label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);

3.3. Momentum Predictor

The momentum predictor is in UEqn.H, with the following code:

UEqn.H

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
// Momentum predictor

  fvVectorMatrix UEqn
  (
      fvm::div(phi, U)
    - fvm::laplacian(nu, U)
  );

  UEqn.relax(); // Equation relaxation, to be discussed later

  solve(UEqn == -fvc::grad(p));

Regarding the relaxation of the momentum equation, this will be discussed later. Readers can comment out this relaxation code and compare the differences in calculation results.

3.4. Pressure and Momentum Correction

The pressure and momentum correction is in pEqn.H, with the following code:

pEqn.H

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
{
    volScalarField rAU(1.0/UEqn.A()); // A^-1 
    volVectorField HbyA(rAU*UEqn.H()); // HbyA
    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA)); // phiHbyA

    fvScalarMatrix pEqn
    (
        fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
    );
	// Pressure correction equation
	// Similarly, fvm discretization returns a matrix, with the discretization operation targeting the parameter p, which is the unknown quantity to be solved
	// fvc discretization returns a field, which serves as a known quantity in the equation


    pEqn.setReference(pRefCell, pRefValue); // Set the pressure reference

    pEqn.solve(); // Solve the pressure correction equation

    // Explicitly relax pressure for momentum corrector
    p.relax(); // Field relaxation, to be discussed later

    // Momentum corrector
    U = HbyA - rAU*fvc::grad(p); // Solve for the corrected velocity
    
    // The corrected pressure and velocity participate in the next SIMPLE loop iteration
}

Pressure also undergoes relaxation, which will be discussed later. Readers can try commenting it out to observe the differences in calculation results.

Tip

Note that for ease of understanding, non-orthogonal correction is not performed here. Since the cavity test case used later has a simple mesh, this has no impact.

3.5. Project Make

As discussed above, this solver does not use additional libraries, so no extra linking specifications are needed in the project Make.

The Make/files content is as follows:

Make/files

1
2
3
ofsp_21_simpleSol.C

EXE = $(FOAM_USER_APPBIN)/ofsp_21_simpleSol

The Make/options content is as follows:

Make/options

1
2
3
4
5
6
7
EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

EXE_LIBS = \
    -lfiniteVolume \
    -lmeshTools

3.6. Compilation

Run the following command in the terminal to compile the entire project:

terminal

wclean
wmake

Compilation is successful with no issues.

3.7. Test Case

We adjust the copied test case to test the above solver.

Modify the control dictionary controlDict as follows:

controlDict

1
2
3
4
5
6
...
application     ofsp_21_simpleSol;
...
endTime         1.0; // Extend the time to compare calculation effects with and without relaxation
// If relaxation is not performed, the calculation may diverge and terminate after some time
...

Modify the solution dictionary fvSolution as follows:

fvSolution

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}


SIMPLE // Specify the dictionary parameters for the SIMPLE algorithm
{
    nNonOrthogonalCorrectors 0;
    consistent      yes;
    pRefCell        0;
    pRefValue       0;

    residualControl
    {
        p               1e-2;
        U               1e-3;
    }
}

relaxationFactors // Specify the dictionary parameters for relaxation
{
    fields // Field relaxation
    {
        p               0.3;   // Pressure field relaxation factor 0.3
        ".*"       0.7;   // Use regular expression to match all p_rgh variants
    }
    
    equations // Equation relaxation
    {
        U               0.9; // 0.9 is more stable but 0.95 more convergent
        ".*"            0.9; // 0.9 is more stable but 0.95 more convergent
    }
}

Other files remain unchanged.

3.7. Execution

Since the solver project compiles without issues, it can be run directly using the scripts.

Run the project:

terminal

./Allclean
./Allrun

Post-process visualization:

terminal

paraFoam -case debug_case

The calculation results can be viewed in ParaView.

Readers can modify parts of the code to explore the effects of different statements.

4. Summary

Through the implementation and discussion of the solver project, I believe we now have a relatively comprehensive understanding of the SIMPLE algorithm and solver implementation.

In the next section, we will discuss the PISO algorithm.

This section has completed the following discussions:

Implementation of the SIMPLE algorithm
Computing the cavity test case

Support us

Tip

Hopefully, the sharing here can be helpful to you.

If you find this content helpful, your comments or donations would be greatly appreciated. Your support helps ensure the ongoing updates, corrections, refinements, and improvements to this and future series, ultimately benefiting new readers as well.

The information and message provided during donation will be displayed as an acknowledgment of your support.

Support

AliPay

Copyright @ 2026 Aerosand
Course (text, images, etc.)：CC BY-NC-SA 4.0
Code derived from OpenFOAM：GPL v3
Other code：MIT License

Last updated on May 6, 2026 • Aerosand

20_simple 22_piso