Adding a Scalar Transport Equation to a Custom Solver in OpenFOAM 12

Introduction

OpenFOAM is an incredibly versatile open-source CFD toolbox that allows you to modify existing solvers or create entirely new ones. In this tutorial, I’ll walk you through the steps of taking an existing solver (foamRun introduced in OpenFOAM 12), duplicating it, renaming it, and adding a scalar transport equation. This process will help you understand the solver structure, the compilation process, and how to integrate new fields into your own custom solver. Prerequisites are:

• OpenFOAM 12 installed and sourced properly.
• A basic understanding of the OpenFOAM directory structure and C++ code organization.
• Familiarity with the standard OpenFOAM solvers and how they operate.

1. Setting Up Your Working Environment

Starting from the OpenFOAM 12 root directory, navigate to the solver directory:

cd $WM_PROJECT_DIR/applications/solvers

Here, you’ll find various solver families. We’ll work with a solver called foamRun. First, copy the foamRun solver into a directory of your choice where you can make modifications:

cp -r foamRun $FOAM_USER_APPBIN/myFoamRunSolver

This creates a copy that you can customize without affecting the system-level installations.

2. Understanding the Make Folder

Inside your newly copied solver folder, you’ll see a Make directory containing two files: files and options.

“files” file:
This file lists which source files should be compiled into the final executable and what the resulting executable name should be. It typically includes .C files and may contain directives that specify linking rules.

Currently, it looks something like this (from the original solver):

EXE = $(FOAM_APPBIN)/foamRun

We need to rename the solver to myFoamRun to avoid conflicts with any existing bash commands or executables. Inside files, change:

EXE = $(FOAM_APPBIN)/foamRun

to

EXE = $(FOAM_USER_APPBIN)/myFoamRun

The $(FOAM_USER_APPBIN) variable points to the user’s local application binary directory, ensuring that your custom solver’s executable myFoamRun is placed in a user-level directory rather than the system-level directory. This keeps your custom developments separate, organized, and avoids overwriting system executables.

Also, ensure that the main source file in files is renamed from foamRun.C to myFoamRun.C. This step ensures consistency in naming and uniqueness.

options file:

The options file specifies additional compile and link directives, such as including libraries and setting up include paths. It often remains unchanged for simple modifications, as it already references the necessary OpenFOAM libraries required by the solver. For now, no changes are needed in the options file.

3. Reviewing Supporting Files: setDeltaT.C and setDeltaT.H

Before diving into the main solver code, let’s look at two files that handle time-step adjustments: setDeltaT.C and setDeltaT.H.

setDeltaT.C and setDeltaT.H:
These files provide functions (setDeltaT and adjustDeltaT) that dynamically control the time-step size. They look at conditions like whether the solver is transient and what maximum time-step (maxDeltaT) is allowed. Initially, they ensure that the solver doesn’t use an unstable or excessively large time-step. Think of these as small utility functions that your solver calls before and during the simulation loop to maintain stability and optimal runtime conditions. We do not need to modify these files for our scalar transport addition. They serve as good examples of how OpenFOAM neatly separates solver logic from utility functions.

4. Renaming and Inspecting the Main Solver File

Now, rename the main solver file from foamRun.C to myFoamRun.C:

mv foamRun.C myFoamRun.C

In this step, we’ll examine the original solver file foamRun.C closely. We will explain its contents line by line or in logical groups. Then, we’ll show the modified version and highlight where and how the scalar transport equation is integrated.

Original foamRun.C Breakdown:

Header and License Block:

*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright ...
     \\/     M anipulation  |
-------------------------------------------------------------------------------
...
\*---------------------------------------------------------------------------*/

This is the standard OpenFOAM license and header information. It states the distribution and license terms. Every official OpenFOAM file starts with a similar block.

Application and Description Sections:

Application
    foamRun

Description
    Loads and executes an OpenFOAM solver module...

Explanation:

• Application: Declares the executable’s name (foamRun).
• Description: Explains what this solver or application does. In this case, it dynamically loads and runs a chosen solver. The solver name can be set in the controlDict or via command-line arguments.

Usage Examples and Additional Notes:

Usage
    foamRun [OPTION]

  -solver <name>
    Solver name

  -libs '("lib1.so" ... "libN.so")'
    Specify additional libraries

Example usage:
    foamRun -solver fluid

Shows how to run this application, what arguments it accepts, and how you might specify a solver or libraries from the command line. Also gives sample configurations for controlDict.

Include Statements:

#include "argList.H"
#include "solver.H"
#include "pimpleSingleRegionControl.H"
#include "setDeltaT.H"

• argList.H: Handles command-line arguments and parsing.
• solver.H: Provides the base class and functionality for loading and using various solver modules.
• pimpleSingleRegionControl.H: Implements the PIMPLE loop control structure (a hybrid of PISO and SIMPLE algorithms).
• setDeltaT.H: Functions to set and adjust the time-step according to solver constraints.

using namespace Foam;

Brings the main OpenFOAM namespace into scope for convenience.

Main Function and Argument Handling:

int main(int argc, char *argv[])
{
    argList::addOption
    (
        "solver",
        "name",
        "Solver name"
    );

    #include "setRootCase.H"
    #include "createTime.H"

• argList::addOption(“solver”, …): Adds a -solver command-line option.
• “setRootCase.H”: Sets up paths to the case directory, reading from environment variables.
• “createTime.H”: Initializes the runTime object, which manages simulation time steps, controlDict, and data I/O.

Reading and Setting the Solver Name:

word solverName
    (
        runTime.controlDict().lookupOrDefault("solver", word::null)
    );

    args.optionReadIfPresent("solver", solverName);

    if (solverName == word::null)
    {
        args.printUsage();
        FatalErrorIn(args.executable())
            << "solver not specified ..."
            << exit(FatalError);
    }
    else
    {
        solver::load(solverName);
    }

Retrieves the solver keyword from the controlDict.

• If -solver was specified on the command line, it overrides the controlDict setting.
• If no solver is found, it prints usage and exits with a fatal error. Otherwise, it dynamically loads the specified solver library.

Mesh Creation and Solver Instantiation:

#include "createMesh.H"

    autoPtr<solver> solverPtr(solver::New(solverName, mesh));
    solver& solver = solverPtr();

• “createMesh.H”: Builds the computational mesh from the provided case setup.
• solver::New(…): Uses the chosen solver’s factory method to create a solver instance, passing in the mesh. solverPtr is an autoPtr, a smart pointer managing memory. solver is a reference to the solver object.

PIMPLE Control Setup and Time-Step Adjustment:

pimpleSingleRegionControl pimple(solver.pimple);
setDeltaT(runTime, solver);

• Constructs the PIMPLE control object that manages the outer iteration loop for steady/transient simulations.
• Calls setDeltaT to initialize the time-step size.

Main Time Loop:

Info<< nl << "Starting time loop\n" << endl;

    while (pimple.run(runTime))
    {
        pimple.read();

        solver.preSolve();

        adjustDeltaT(runTime, solver);

        runTime++;

        Info<< "Time = " << runTime.userTimeName() << nl << endl;

        // PIMPLE corrector loop
        while (pimple.loop())
        {
            solver.moveMesh();
            solver.motionCorrector();
            solver.fvModels().correct();
            solver.prePredictor();
            solver.momentumPredictor();
            solver.thermophysicalPredictor();
            solver.pressureCorrector();
            solver.postCorrector();
        }

        solver.postSolve();

        runTime.write();

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

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

    return 0;
}

• pimple.run(runTime): Determines if the simulation should continue.
• pimple.read(): Reads updated PIMPLE settings if they’ve changed.
• solver.preSolve()/postSolve(): Hooks for solver-specific operations before and after each iteration.
• adjustDeltaT(runTime, solver): Adjusts time-step if needed.
• runTime++: Increments time and moves to the next time-step.
• PIMPLE corrector loop: Executes multiple correction steps per time-step if required, improving solution accuracy and stability.
• After completing the loops, it writes out results, prints execution times, and repeats until pimple.run(runTime) returns false (end of simulation).

This original file does not include a scalar transport equation. It focuses on loading a solver and running a simulation loop with the PIMPLE algorithm, leaving the actual physics (like momentum, pressure, energy) to the dynamically loaded solver module.

---

Modified Version with Scalar Transport Equation:

In the modified version, we add a separate header file (e.g., tracerEqn.H) and include it at an appropriate point in the time loop. The modifications would look like this:

// After setting and adjusting the time-step and before writing results, we include the scalar transport equation:
#include "tracerEqn.H"

What Does tracerEqn.H Contain?
tracerEqn.H defines and solves the scalar transport equation. It might look like this:

5. Introducing the Scalar Transport Equation

To integrate a scalar transport equation (such as a tracer concentration) into your solver, you need to modify the solver loop to include an equation that represents the scalar’s evolution over time.

Where:

1. Temporal term:
   • This represents the rate of change of the scalar T over time.
   • Represented in the code as fvm::ddt(T).
   • In steady-state, this term is typically zero.
2. Convection term:
   • Describes the transport of scalar T due to the flow field velocity.
   • Represented in the code as fvm::div(phi, T), where phi is the volumetric flux.
3. Diffusion term:
   • Accounts for the diffusion of scalar T within the medium.
   • DT the thermal diffusivity (or a similar diffusion coefficient).
   • Represented in the code as fvm::laplacian(DT, T).

Create a new file called tracerEqn.H and place it in the same directory as myFoamRun.C. This file will define the scalar equation. For example:

fvScalarMatrix tracerEqn
(
    fvm::ddt(tracer)            // Time derivative of tracer field
  + fvm::div(phi, tracer)       // Convective transport of tracer by the flow field
  - fvm::laplacian(DT, tracer)  // Diffusive spreading of tracer
);

tracerEqn.solve();

Info << "Min/Max Tracer: " << min(tracer).value() << "/ " << max(tracer).value() << endl;

The scalar transport equation typically has three main parts:

1. Temporal derivative (fvm::ddt(tracer)): This represents how the scalar changes over time.
2. Convection term (fvm::div(phi, tracer)): This represents how the scalar is carried by the flow field.
3. Diffusion term (fvm::laplacian(DT, tracer)): This accounts for the scalar spreading out due to diffusion.

DT is the diffusivity of the scalar, and phi is the volumetric flux field. Together, these terms form a partial differential equation that OpenFOAM discretizes and solves at each time-step.]

To incorporate tracerEqn.H into your solver, include it in myFoamRun.C at an appropriate location in the time loop, typically after your other primary equations are solved but before writing the results.

6. Compiling Your New Solver

Once you’ve completed all modifications (renamed files, updated EXE path in Make/files, introduced tracerEqn.H and included it in myFoamRun.C), it’s time to compile.

First, ensure your OpenFOAM environment is loaded:

source /opt/openfoam12/etc/bashrc

This command sets up all the required environment variables, paths, and aliases to compile and run OpenFOAM utilities and solvers. If you’ve defined a custom alias for sourcing OpenFOAM, use that instead.

Before compiling, clean any stale build files:

wclean

The wclean command removes previously compiled objects and ensures that you start the compilation process from a clean state. This is helpful when switching code branches, changing compilers, or modifying solver source files to avoid leftover binaries causing unexpected behavior.

Now, compile with:

wmake

The wmake tool is a wrapper around the underlying compiler commands, making it straightforward to rebuild OpenFOAM applications and libraries. After running wmake, you should see your new myFoamRun solver in the directory specified by FOAM_USER_APPBIN.

7. Setting Up Your Case

For testing the newly developed solver, let’s consider a prepared pitzDaily case. Before running the simulation, ensure the following steps are completed:

1. Adjusting the Diffusivity:
   In the constant/transportProperties file, add or modify the DT (diffusivity) entry for the tracer field. This specifies how quickly the tracer diffuses through the flow field.
2. Initializing the Tracer Field:
   In the 0 directory, create an initial conditions file for tracer—similar to how U (velocity) or p (pressure) fields are defined. Set its internal field and boundary conditions according to your problem setup. For example, you may start with a uniform tracer concentration or a specified gradient.

Once these preparations are complete, run the solver from the case directory:

myFoamRun

If everything is set correctly, you will observe your tracer field evolving alongside the primary flow variables (e.g., velocity and pressure) as the simulation progresses.

Conclusion

Congratulations! You’ve successfully taken an existing solver in OpenFOAM, customized it, and introduced a scalar transport equation. This journey gave you insight into the solver’s internal structure, file organization, and the flexibility of OpenFOAM’s framework. With your new solver running the pitzDaily case, you can now visualize how the tracer diffuses and interacts with the flow field.

For even easier access and collaboration, both the solver and the test case are available on GitHub, ready for you to download, experiment with, and improve.

Happy coding and happy simulating! 🎉🚀