Posted on 16th August 2019 in CFD

Adding a passive scalar transport equation to pimpleFoam

In this post we will provide a step-to-step guide to implement a passive scalar transport equation within the pimpleFoam solver in OpenFOAM, in order to extend its capabilities to problems like forced turbulent convection heat transfer and generic turbulent convection-diffusion transport.

OpenFOAM is a an extremely popular open-source CFD solver including a broad range of built-in models and solvers. PimpleFOAM is its go-to solver for unsteady turbulent problems. PimpleFOAM relies on the PIMPLE algorithm, which is an extension of the unsteady PISO algorithm with the addition of a SIMPLE-like outer correction loop in order to relax the limitation of a maximum unity Courant number typical of the PISO-based solvers. More details about the PIMPLE algorithm can be found here.

The pimpleFoam solver implemented in OpenFOAM is suitable for the simulation of transient, turbulent incompressible flows of Newtonian fluids in isothermal conditions. Our goal is to add an additional passive scalar transport equation within the solver. The final implementation can be found in this GitHub repository. Both the code within the respository and this step-by-step guide are based on OpenFOAM v6.1 and they should be usable with any reasonably recent OpenFOAM distribution.

Since heat transfer is by far the most popular application for such a solver, we will indicate our passive scalar with a capital T, standing for temperature. The transport equation we seek to implement in pimpleFoam is

(1)   \begin{equation*}\frac{\partial T}{\partial t}+U_j \frac{\partial T}{\partial x_j}= \frac{\partial}{\partial x_j}\left(\alpha \frac{\partial T}{\partial x_j} - \overline{\theta u_j} \right) \end{equation*}

where \alpha=\nu/Pr and Pr are the molecular thermal diffusivity and the molecular Prandtl number of the fluid, respectively. If T is a generic passive scalar, then \alpha represents its molecular diffusivity and Pr corresponds to its molecular Schmidt number. The term \overline{\theta u_j} represents the turbulent flux of the passive scalar T and appears as an additional unclosed term in the standard convection-diffusion equation as a consequence of the averaging operation within the RANS turbulence modelling framework. A similar unclosed term representing the subgrid-scale turbulent flux of T would appear in the equation as a consequence of the spatial filtering operation performed in the LES turbulence modelling framework. This term requires a suitable model in order to close and solve Equation 1. In our implementation we will use a Simple Gradient Diffusion Hypothesis, which expresses the turbulent scalar flux as a function of the scalar gradient as

(2)   \begin{equation*}\overline{\theta u_j}=-\alpha_t\frac{\partial T}{\partial x_j} \end{equation*}


where \alpha_t is the turbulent thermal diffusivity (or the turbulent diffusivity for a generic passive scalar) which will be evaluated following the so-called Reynolds analogy

(3)   \begin{equation*}\alpha_t=\frac{\nu_t}{Pr_t}\end{equation*}

In Equation 3 Pr_t is the turbulent Prandtl number (which corresponds to the turbulent Schmidt number for a generic passive scalar) and is a assumed to be a constant to be supplied by the user. Usually Pr_t is taken in the range 0.8-1.0.

Combining Equations 2 and 3 together and plugging the result into Equation 1 leads to the final version of the transport equation to be implemented in pimpleFoam, i.e.:

(4)   \begin{equation*}\frac{\partial T}{\partial t}+U_j \frac{\partial T}{\partial x_j}=   \frac{\partial}{\partial x_j}\left(\alpha_{eff}  \frac{\partial T}{\partial  x_j} \right) \end{equation*}

where \alpha_{eff}=\alpha+\alpha_t is the effective thermal diffusivity, accounting for the effects of both molecular and turbulent transport. In the case of a generic passive scalar different from the temperature, \alpha_{eff} represents the sum of the molecular and of the turbulent mass diffusivities.

Our goal is now to implement Equation 4 in the native pimpleFoam solver to create a new solver that will be called (not surprisingly) passiveScalarPimpleFoam. Instead of creating a new solver from scratch, it is much more convenient to start from the existing implementation of pimpleFoam and simply plug-in the necessary lines of code to include our additional passive scalar equation in the solver’s algorithm. Firstly, move to the folder where you want to store the source code for your new solver, which usually is $WM_PROJECT_USER_DIR/applications/solvers:
cd $WM_PROJECT_USER_DIR/applications/solvers

From within this folder create a copy of the source code for pimpleFoam with the following command:
cp -r $FOAM_APP/solvers/incompressible/pimpleFoam .
rename the folder as passiveScalarSimpleFoam
mv pimpleFoam passiveScalarPimpleFoam
move inside the folder you have just created
cd passiveScalarPimpleFoam
and remove stuff we are not going to need for the new implementation, i.e:
rm -r SRFPimpleFoam

At this point, we are ready to start to convert the existing implementation of pimpleFoam into our own passiveScalarPimpleFoam. The first step to perform is to rename any file in the source code which has “pimpleFoam” in its name to “passiveScalarPimpleFoam”. In this case, there is only one such file, i.e. pimpleFoam.C, and this operation could be done manually. For more complex codes, it is advisable to do this in an automated way to save time and avoid possible issues due to oversights, so we are going to use the following command useful to bulk-rename files within the Linux shell
rename 's/pimpleFoam/passiveScalarPimpleFoam/' *
which should have renamed your pimpleFoam.C to passiveScalarPimpleFoam.C

Before performing any other operation, we will have a look at the content of the Make folder:
cd Make
This folder contains instructions for the compiler in order to compile our code and build the executable file associated with it. Inside the folder you will see two files, i.e. files and options. The name of the source files to be compiled and the path and name for the executable that will be created are specified within files, whilst options contains the names and paths of the header files included in our code as well as the names and paths to the libraries to be dynamically linked to our code. No changes are needed to the options file in this case, since we are going to use the same header files and the same libraries that are used to build the native pimpleFoam solver. Within files, on the other hand, we will have to specify the path to the folder where we want to generate the new executable file. Currently this is $FOAM_APPBIN, which is located inside the OpenFOAM installation folder and is meant to contain the executables native to the current OpenFOAM distribution. It is good practice to place the executables generated by the user in a different location, and the standard practice is to use $FOAM_USER_APPBIN for this purpose. Therefore, using your favourite text editor, inside files replace:
$(FOAM_APPBIN)
with
$(FOAM_USER_APPBIN)
and save the changes.

In addition, we want to change the name of the .C source file and of the executable from pimpleFoam to passiveScalarPimpleFoam. Again we could do this manually, but what we actually want to do is to make sure that any occurrence of pimpleFoam within every single file in our source code is replaced by passiveScalarPimpleFoam. It is much safer and quicker to leverage the power of the Linux shell to perform the search and replace operation automatically. In order to do so move back from the Make folder to the main folder of our source code:
cd ..
and issue the following command
find . -type f -exec sed -i 's/pimpleFoam/passiveScalarPimpleFoam/g' {} +
Finally, double check within Make/files that the two pimpleFoam occurences have now been replaced with their passiveScalarPimpleFoam counterparts.

We are now ready to start with the actual implementation of the new solver. By going back to the derivation of Equation 4, we can see that we are going to need two additional scalar fields in our solver: the first one is, of course, the passive scalar T that we want to solve for. The second one is the the turbulent thermal diffusivity \alpha_t. Somehow, we have to tell OpenFOAM that these two new fields exist. More technically, both these variables have type volScalarField so we will have to invoke the constructor of the volScalarField class in order to instantiate them. In all OpenFOAM solvers, basic fields are created within the createFields.H file. Open this file and add the following lines at the top of the file to create the field associated with T:
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

Similarly, add the following lines to createFields.H before the #include “createMRF.H” statement located at the bottom of the file to instantiate the turbulent thermal diffusivity field:
Info<< "Reading field alphat\n" << endl;
volScalarField alphat
(
IOobject
(
"alphat",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

In the original pimpleFoam solver, the working fluid is completely defined by its molecular viscosity \nu. In our extended solver, we also need to define the molecular thermal diffusivity, which will be done by imposing the molecular Prandtl number of the fluid. In addition, we will also have to specify a value for the the turbulent Prandtl number Pr_t, which is needed in Equation 3 to evaluate \alpha_t. In order to do this, replace the line
singlePhaseTransportModel laminarTransport(U, phi);
in createFields.H with
#include "readTransportProperties.H"

Now, create a readTransportProerties.H file in your text editor, add the following content to the file:
singlePhaseTransportModel laminarTransport(U, phi);
// Laminar Prandtl number
dimensionedScalar Pr("Pr", dimless, laminarTransport);
// Turbulent Prandtl number
dimensionedScalar Prt("Prt", dimless, laminarTransport);

and save it within the main folder of the source code.

Finally, we have to code the actual transport equation for the passive scalar, i.e. Equation 4; thanks to OpenFOAM’s high-level syntax, the coded equation will resemble very closely its “textbook” counterpart. Create a TEqn.H file in your text editor, add the following content to the file
{
alphat = turbulence->nut()/Prt;
alphat.correctBoundaryConditions();
volScalarField alphaEff("alphaEff", turbulence->nu()/Pr + alphat);
fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(alphaEff, T)
);
TEqn.relax();
TEqn.solve();
}

and save it within the main folder of the source code.

The last remaining step is to include our new equation within the solver’s algorithm. Since the passive scalar has no impact on the resolved flow-field we will solve for its transport equation outside of the PIMPLE loop. In order to do this open passiveScalarPimpleFoam.C, add the following line:
#include "TEqn.H"
after the closed curly bracket which ends the PIMPLE loop and immediately above the runTime.write() line and save the changes.

Now it is time to compile the source code to build the executable associated with our solver. If you are unsure about any of the steps above double-check them against the source code provided in the GitHub repository. When you are ready to go, just issue the command
wmake
from the main source code folder. When the compilation is complete (hopefully without errors) type
passiveScalarPimpleFoam
in your terminal. OpenFOAM should execute and throw a FOAM_FATAL_ERROR due to a missing controlDict file. If this is the outcome that you get … congratulations, you have just correctly compiled your solver!

Now it is time to put our brand new solver to the test. In order to do this, we will use a modified version of the existing pitzDaily tutorial to include the transport of a passive scalar. To do so download the passiveScalarPitzDaily folder from the GitHub repository folder to your OpenFOAM $FOAM_RUN folder. This test case mimics the classical isothermal pitzDaily tutorial which is shipped with any OpenFOAM distribution. The case has been expanded by including the transport of a passive scalar which has a uniform 0 value at the inlet and a uniform value of 1 at the walls. The initial and boundary conditions for the “new” scalars T and \alpha_t are specified within the T and alphat files within the 0 folder. The molecular Prandtl number of the working fluid is 0.71, which is equal to that of air under standard conditions, and is specified in the transportProperties file within the constant folder. In the same file, a value of 0.85 is imposed for the turbulent Prandtl number.

To run the case go into the case folder:
cd passiveScalarPitzDaily
and execute the solver you have just created:
passiveScalarPimpleFoam > log &

The run should take a couple of minutes on a standard laptop. To check the log file during the run type
tail -f log

Once the run is complete, launch paraFoam
paraFoam
to analyse the results. The velocity field is not affected by the presence of the passive scalar (hence the “passive” in its name…) and should look similar to the the figure below.

passiveScalarPitzDaily test case: velocity magnitude contours

The temperature field shown below illustrates the heating up of the fluid, which enters the domain with a temperature equal to 0, due to the presence of the “hot” walls.

passiveScalarPitzDaily test case: temperature contours for Pr=0.71

Lastly, we are going to investigate the effect of the thermal molecular diffusivity of the fluid on the calculated temperature field by lowering the fluid Prandtl number to 0.01. To run the simulation at this lower Prandtl number, first create a “clone” of the previous case in your run folder
cd $FOAM_RUN
foamCloneCase passiveScalarPitzDaily passiveScalarPitzDaily_lowPr
cd passiveScalarPitzDaily_lowPr

Now modify the value of the Pr parameter in the transportProperties files within the constant folder from 0.71 to 0.01 and run the case
passiveScalarPimpleFoam > log &

This lower value of Pr corresponds to a much higher thermal diffusivity compared to the previous case. Such a low Prandtl value is typical of fluids like liquid metals, which are excellent heat conductors. This should result in a much thicker thermal boundary layer at the walls, which is indeed the effect that is observed in the temperature contours showed below.

passiveScalarPitzDaily test case: temperature contours for Pr=0.01

This concludes our tutorial. Try to test the new solver with your own cases and to modify it, if needed. If you found this post helpful feel free to share it and to ask questions/give us a feedback in the comments section.


Leave a Reply