How to add temperature to icoFoam
Contents
1 How to add temperature transport to icoFoam
This HOWTO will cover rudimentary methods for altering an existing solver (icoFoam) to solve thermal transport. The main items to be accomplished are first, to copy and test that your installation of OpenFOAM can compile the existing solver correctly. Once that is accomplished, various small changes are necessary to the solver files themselves and these will be detailed below. Finally, new fields have to be added to the initial and boundary conditions and alterations to the fvSchemes file.
2 Copy and recompile icoFoam
This will go through the steps of creating a personal version of icoFoam in the user's subdirectory. The first step is to insure that your installation of OpenFOAM works properly and compiles the unedited solver. First, bring up a console and move to your OpenFOAM installation folder.
cd OpenFOAM
Your particular OpenFOAM installation folder will have a version number following it such as:
cd OpenFOAM-1.7.1
OpenFOAM has organized the solvers separate from the source code of OpenFOAM calling them 'applications'. Inside the 'applications' folder, there are subdirectories including one for solvers.
cd applications/solvers/
A shortcut you can use is stored in the environment variable FOAM_SOLVERS:
cd $FOAM_SOLVERS
Have a look around, but this tutorial is based on icoFoam, which we will copy to our own location:
cd incompressible
If you don't have a solver directory in $WM_PROJECT_USER_DIR do
mkdir -p $WM_PROJECT_USER_DIR/applications/solvers cp -r icoFoam $WM_PROJECT_USER_DIR/applications/solvers/my_icoFoam cd $WM_PROJECT_USER_DIR/applications/solvers/my_icoFoam
Now, a couple of alterations need to be made to the make files in order for everything to compile and not overwrite the original solver. First, rename the primary file to your new solver name and delete the old dependency file:
mv icoFoam.C my_icoFoam.C rm icoFoam.dep
Now go into the Make subdirectory and open the 'files' file with your favorite editor. Change it to read:
my_icoFoam.C EXE = $(FOAM_USER_APPBIN)/my_icoFoam
No changes are necessary for the 'options' file. Delete the old binaries subdirectory:
rm -rf linuxGccDP0pt cd ..
Now, test that the renamed solver (and your installation of OpenFOAM) works:
wmake
If everything worked correctly, your new solver binary should appear in the FOAM_USER_APPBIN directory. Check this with:
ls $FOAM_USER_APPBIN
3 Adding the temperature field
This will show the steps to add another field variable to a solver. Open the my_icoFoam.C (or whatever you have named it) with your favorite text editor.
First, edit the 'Application' to reflect the new name.
Following the flow of the program, one notices that the header file 'createField.H' is called prior to the solution loop. This file was copied with the solver and has the specific information pertaining to what variables will be solved.
Open createFields.H in your favorite editor.
The first items loaded is the kinematic viscosity from the transportProperties dictionary file. We will add a new transport property related to the thermal diffusion which will be denoted as DT. Make the following edits:
Info<< "Reading transportProperties\n" << endl; IOdictionary transportProperties ( IOobject ( "transportProperties", runTime.constant(), mesh, IOobject::MUST_READ, IOobject::NO_WRITE ) ); dimensionedScalar nu ( transportProperties.lookup("nu") ); //Add here... dimensionedScalar DT ( transportProperties.lookup("DT") ); //Done for now...
Aside: What does this code do and how do we learn about it?
Later on, we will need to edit the dictionary file to reflect this change.
Following this in the file there are lines which pertain to the creation of the pressure (p) and velocity (U) fields. We will add a new field for temperature (T). The fastest way to do this is to copy and paste the pressure lines and then edit them appropriately like so:
Info<< "Reading field T\n" <<endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh );
Save these changes. You've completed adding a new field variable, which will be called 'T' to the solver.
4 Adding a new equation to solve
The next step is to add a new equation describing the transport of the temperature. Return to editing the my_icoFoam.C file (or whatever you named it).
Because the temperature transport depends on the velocity field, we will add the equation after the momentum equation is solved (after the PISO loop), but before the time step is written. Edit your file so it looks like this:
# include "continuityErrs.H"
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}
//add these lines...
fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
);
TEqn.solve();
//done adding lines...
runTime.write();
These lines add a new equation for the temperature and make use of the face flux variable, phi, which is already used in the momentum equation solution.
Save your changes and run wmake:
wmake
Your computer should then produce a newly compiled binary.
5 Add a new file for initial and boundary conditions
The next step in this process is to test the new solver. This will first be accomplished by modifying the existing cavity tutorial. Later, we will benchmark the solver against an analytical solution.
First, copy the cavity tutorial files to a new folder:
cd $FOAM_RUN/tutorials/incompressible/icoFoam cp -r cavity $FOAM_RUN/my_icoFoam_cavity cd $FOAM_RUN/my_icoFoam_cavity
We will start by editing the transportProperties dictionary:
cd constant
Open the transportProperties dictionary in your favorite editor and add the following line under the definition of nu:
DT DT [0 2 -1 0 0 0 0] 0.002;
The value doesn't really matter at this point.
Next, we will create a new initial and boundary conditions file for the temperature field:
cd ../0 cp p T
Open the 'T' file in your favorite editor.
For this initial test, we will simply make the moving wall and the fixed walls different temperatures. Edit the file so it looks like the following:
class volScalarField;
object T;
}
//**********************************//
dimensions [0 0 0 1 0 0 0];
internalField uniform 300;
boundaryField
{
movingWall
{
type fixedValue;
value uniform 350;
}
fixedWalls
{
type fixedValue;
value uniform 300;
}
frontAndBack
{
type empty;
}
}
Save those changes. There is one more change to be made before we can run the solver.
6 What to add in fvSchemes and fvSolution
By adding a new equation to solve, we need to tell OpenFOAM what discretization schemes to apply to the equations. This is done in the fvSchemes dictionary. Enter the 'system' subdirectory for this case:
cd ../system
Open the 'fvSchemes' dictionary with your favorite editor.
Now, there are two main items added in the thermal transport equation above: a divergence term and a laplacian. Under each of these sections in the dictionary, we need to add terms which match what was added in the solver source code. Edit your file so it includes the following:
divSchemes
{
default none;
div(phi,U) Gauss linear; //NOTICE: there is no space between the comma and the variables
div(phi,T) Gauss upwind;
}
laplacianSchemes
{
default none;
laplacian(nu,U) Gauss linear corrected;
laplacian((1|A(U)),p) Gauss linear corrected;
laplacian(DT,T) Gauss linear corrected;
}
Alternatively, a scheme can be added to the 'default' field instead of adding specific ones below.
Save these changes.
Next, open the fvSolution dictionary in your favorite editor. It should look like this:
solvers
{
p PCG
{
//information about the pressure solver
};
//add this...
T
{
solver BICCG;
preconditioner DILU;
tolerance 1e-7;
relTol 0;
};
//done editing...
Leave the velocity information intact and save these changes.
Review the controlDict file to remind yourself what it looks like. You are now ready to run the case and test your solver.
cd .. my_icoFoam
Your case should run and display the usual information in the window.
The solution for time-step 0.5 should look like this:
Media:CaseFile.tar.gz
--Media:My_icoFoam_cavity_OF-1.7.0.tar.gz case according OpenFoam-1.7.0
Media:SolverSource.tar.gz --Mike Jaworski 06:49, 25 February 2008 (CET)
7 Benchmarking your new solver
A case study of the Blasius flat-plate flow problem is shown using this solver on another wiki-page. Blasius Flat-Plate Flow Benchmark
