How to add a passive scalar to your OpenFOAM simulations

How to add a passive scalar to your OpenFOAM® simulations

The purpose of this post is to show how to add a passive scalar to your simulations without modifying a given solver. This feature is, for example, very useful when trying to track the concentration of inert species.

The tutorial Cavity

The Cavity tutorial is part of the icoFoam solver tutorials. The tutorial is located here

[pastacode lang=”bash” manual=”%24FOAM_TUTORIALS%2Fincompressible%2FicoFoam%2Fcavity%2Fcavity%2F” message=”” highlight=”” provider=”manual”/]

The Cavity tutorial presents the modeling of the flow driven by a moving lid in a square cavity (lid-driven cavity flow). This fluid mechanics problem is a typical validation test case for CFD codes. The problem is considered 2D and the boundary conditions are shown in the figure below. The temporal evolution of the velocity field is also represented.

Reminder about the icoFoam solver

The icoFoam solver of OpenFOAM is an unsteady, incompressible, one phase solver adapted to laminar flows. The solver is based on the PISO loop (and therefore necessarily unsteady). The solver algorithm is described in detail here.

Changing the cavity tutorial - Adding a passive scalar

In order to illustrate the possibility of adding a passive scalar to an OpenFOAM simulation, it is assumed that there is an inert chemical species (called s, for example a rare gas) and modeled by a scalar field whose evolution is governed by the convection / diffusion equation

$\frac{\partial s}{\partial t} + \overrightarrow{U}\overrightarrow{\nabla s} = \nabla^2 (D s)$

The concentration of the species is assumed to be zero in the cavity at t = 0 sec and set at 1 on the upper boundary condition representing the moving lid.

The definition of the scalar field is done in the controlDict file located in the /system directory:

[pastacode lang=”cpp” manual=”functions%0A%7B%0A%09scalar1%0A%09%7B%0A%09type%20%20%20%20%20%20%20%20%20%20%20%20scalarTransport%3B%0A%09functionObjectLibs%20(%22libsolverFunctionObjects.so%22)%3B%0A%09enabled%20true%3B%0A%09writeControl%20outputTime%3B%20%2F%2F%20write%20scalar%20field%20results%0A%09%0A%09field%20s%3B%20%2F%2Fname%20of%20scalar%20field%0A%09nCorr%201%3B%20%2F%2Fnumber%20of%20corrector%20loop%0A%09D%200.0001%3B%20%2F%2Fdifussion%20coefficient%0A%2F%2F%20%20%20%20schemesField%20U%3B%20%20%2F%2F%20Name%20of%20field%20to%20use%20when%20looking%20up%20schemes%20from%20fvSchemes%0A%0A%09log%20yes%3B%0A%09%7D%0A%7D%3B” message=”Scalar Definition” highlight=”” provider=”manual”/]

In this case, schemesField U being commented the discretization scheme used for the convection term has to to defined in fvScheme the (if not OpenFOAM® uses that defined for U). We choose to use the second order scheme linearUpwind (grad(s) being calculated with the scheme defined in the section gradSchemes, ie a centered scheme). Note also that the scheme used here is not as stable as possible and more complex cases would require an upwind scheme (1st order) for stability reason.

[pastacode lang=”cpp” manual=”div(phi%2Cs)%20%20%20%20%20Gauss%20linearUpwind%20grad(s)%3B” message=”Divergence scheme for passive scalar field” highlight=”” provider=”manual”/]

Once the fvSchemes file is correctly adapted, the resolution algorithms of the matrix system need to be set in the fvSolution file:

[pastacode lang=”cpp” manual=”%20%20%20%20s%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20solver%20%20%20%20%20%20%20%20%20%20PBiCGStab%3B%0A%20%20%20%20%20%20%20%20preconditioner%20%20DILU%3B%0A%20%20%20%20%20%20%20%20tolerance%20%20%20%20%20%20%201e-08%3B%0A%20%20%20%20%20%20%20%20relTol%20%20%20%20%20%20%20%20%20%200%3B%0A%20%20%20%20%20%20%20%20minIter%09%20%20%20%20%201%3B%0A%20%20%20%20%7D” message=”Linear-solver for passive scalar s” highlight=”” provider=”manual”/]

Finally, it only remains to define the boundary conditions used for the field s by creating a file “s” in the directory /0:

[pastacode lang=”cpp” manual=”%2F*——————————–*-%20C%2B%2B%20-*———————————-*%5C%0A%7C%20%3D%3D%3D%3D%3D%3D%3D%3D%3D%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%7C%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%7C%0A%7C%20%5C%5C%20%20%20%20%20%20%2F%20%20F%20ield%20%20%20%20%20%20%20%20%20%7C%20OpenFOAM%3A%20The%20Open%20Source%20CFD%20Toolbox%20%20%20%20%20%20%20%20%20%20%20%7C%0A%7C%20%20%5C%5C%20%20%20%20%2F%20%20%20O%20peration%20%20%20%20%20%7C%20Version%3A%20%20plus%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%7C%0A%7C%20%20%20%5C%5C%20%20%2F%20%20%20%20A%20nd%20%20%20%20%20%20%20%20%20%20%20%7C%20Web%3A%20%20%20%20%20%20www.OpenFOAM.com%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%7C%0A%7C%20%20%20%20%5C%5C%2F%20%20%20%20%20M%20anipulation%20%20%7C%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%7C%0A%5C*—————————————————————————*%2F%0AFoamFile%0A%7B%0A%20%20%20%20version%20%20%20%20%202.0%3B%0A%20%20%20%20format%20%20%20%20%20%20ascii%3B%0A%20%20%20%20class%20%20%20%20%20%20%20volScalarField%3B%0A%20%20%20%20object%20%20%20%20%20%20s%3B%0A%7D%0A%2F%2F%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20*%20%2F%2F%0A%0Adimensions%20%20%20%20%20%20%5B0%200%200%200%200%200%200%5D%3B%0A%0AinternalField%20%20%20uniform%200%3B%0A%0AboundaryField%0A%7B%0A%20%20%20%20movingWall%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20type%20%20%20%20%20%20%20%20%20%20%20%20fixedValue%3B%0A%20%20%20%20%20%20%20%20value%20%20%20%20%20%20%20%20%20%20%20uniform%201.0%3B%0A%20%20%20%20%7D%0A%0A%20%20%20%20fixedWalls%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20type%20%20%20%20%20%20%20%20%20%20%20%20zeroGradient%3B%0A%20%20%20%20%7D%0A%0A%20%20%20%20frontAndBack%0A%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20type%20%20%20%20%20%20%20%20%20%20%20%20empty%3B%0A%20%20%20%20%7D%0A%7D%0A%0A%2F%2F%20*************************************************************************%20%2F%2F” message=”Boundary condition for passive scalar field” highlight=”” provider=”manual”/]

Results

The two images below show the passive scalar evolution for the modified cavity tutorial of OpenFOAM. Two different values of diffusion coefficients are tested:

$D = 0.0001 m^2/s$ and $D = 0.0005 m^2/s$

Conclusions

It is easily possible to add a passive scalar to a simulation under OpenFOAM® by following the steps previously described:

Define the scalar field in the controlDict file
Define the discretization schemes in the fvSchemes file
Set the linear solver in the fvSolution file
Create a file with the name of the scalar in directory 0 and set the boundary conditions

Note that the procedure can be repeated to add as many passive scalars as required !