Project 3- Adding a temperature dictionary to the incompressible icoFoam case to solve the energy equation

Project 3 : Adding a temperature dictionary to the incompressible icoFoam case to solve the energy equation

1. Introduction

The aim of this project is to add a temperature dictionary to the incompressible icoFoam case to solve the energy equation as given below.

∂T∂t+(∇.U)T=k∇2T">∂T∂t+(∇.U)T=k∇2T∂T∂t+(∇.U)T=k∇2T .....................(i)

where kappa "k" is the heat diffusion coefficient.

The default icoFoam solver only solves for continuity and momentum equation. In this project eqn(i) for temperature will also be solved.

2. Changes in the solver:

IcoFoam solver was copied to a user defined folder to modify it.

The changes that were made in the icoFoam solver are given below.

i) The code given below was added to the existing icoFoam solver file icoFoam.C. The code was added after the PISO loop in the file.

  
//Energy equation for temperature

fvScalarMatrix TEqn
        (
            fvm::ddt(T)                   //Temporal term
           +fvm::div(phi,T)               //Convective term
           -fvm::laplacian(kappa,T)       //Diffusive term
        );

        TEqn.solve();

ii) In createFields.H file temperature T was added as an IO object with the following code. The code given below is similar to Pressure. Hence same code was copied and P(pressure) was replaced with T.

volScalarField T
(
    IOobject
    (
        "T",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

iii) Again in createFields.H file a dictionary reader for kappa was created as given below.

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

iv) icoFoam.C was renamed as icoTempFoam.C. Also the changes were made in the Make>>files directory.

With the above changes, the solver was compiled using the command wmake and it was successfully compiled.

3. Changes in the case:

$FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity case was used to test out our solver and renamed it to cavityTemp.

i) In the folder 0, T was added with the following boundary conditions.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 300;

boundaryField
{
    movingWall
    {
        type            fixedValue;
        value           uniform 400;
    }

    fixedWalls
    {
        type            fixedValue;
        value           uniform 300;
    }

    frontAndBack
    {
        type            empty;
    }
}


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

ii) Heat diffusion coefficient was set as : kappa [ 0 2 -1 0 0 0 0 ] 0.001

iii) Added a div(phi,T) in divSchemes and laplacian(kappa,T) in laplacianSchemes in the fvSchemes folder.

iv) Added a smoothSolver linear solver for temperature in fvSolution and a tolerance of 1e-07 and a relTol of 0 as shown below.

    T
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-07;
        relTol          0;
    }

v) The name of the solver was changed from icoFoam to icoTempFoam in controlDict file.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LaTeX Example: Project Report
\documentclass{article}

% Language setting
% Replace `english' with e.g. `spanish' to change the document language
\usepackage[english]{babel}

% Set page size and margins
% Replace `letterpaper' with`a4paper' for UK/EU standard size
\usepackage[letterpaper,top=2cm,bottom=2cm,left=3cm,right=3cm,marginparwidth=1.75cm]{geometry}

% Useful packages
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{subfigure}

\title{\textbf{Project-3:To add a temperature dictionary to the incompressible {\fontfamily{qcr}\selectfont icoFoam} solver to solve the energy equation}}


\begin{document}
\maketitle

\begin{abstract}
In this project we are going to add the temperature equation into the existing {\fontfamily{qcr}\selectfont icoFoam} solver and modify it into the {\fontfamily{qcr}\selectfont icoTempFoam}. Steps for the same are described here.
Then this {\fontfamily{qcr}\selectfont icoTempFoam} is tested on the modified tutorial case {\fontfamily{qcr}\selectfont cavityTemp}.The details of the modification of tutorial case are explained in this report. 
The results of the successful implementation of the modified solver are presented. 
\end{abstract}

\section{Introduction to the new solver {\fontfamily{qcr}\selectfont icoTempFoam}}
The {\fontfamily{qcr}\selectfont icoFoam} solver solves the incompressible flow equations for mass and momentum:
\begin{equation}
    \nabla . U = 0 
\end{equation}
\begin{equation}
    \frac{\partial U}{\partial t} + \nabla . U = - \nabla p + \nabla.(\nu \nabla U) + g 
\end{equation}
Our goal is to add an equation for temperature as:
\begin{equation}
    \frac{\partial T}{\partial t} + (U.\nabla) T  =  \nabla.(\kappa \nabla T)  
\end{equation}
where $\kappa $ is the heat diffusion coefficient,
\newline into the existing {\fontfamily{qcr}\selectfont icoFoam} solver which then will become the {\fontfamily{qcr}\selectfont icoTempFoam} solver, as the name suggests the temperature equations are added to the {\fontfamily{qcr}\selectfont icoFoam} solver.
\section{Steps in the creation \& compilation of the new solver \\{\fontfamily{qcr}\selectfont icoTempFoam} .}

\subsection{To copy and rename the files from existing {\fontfamily{qcr}\selectfont icoFoam} solver}

Copy the {\fontfamily{qcr}\selectfont icoFoam} solver source code from the  {\fontfamily{qcr} \selectfont \$FOAM\_SOLVERS/incompressible/icoFoam} \newline using copy command to the {\fontfamily{qcr} \selectfont \$WM\_PROJECT\_USER\_DIR/my\_solvers}. \newline If  {\fontfamily{qcr} \selectfont my\_solvers}  directory is not already created then make it and proceed further i.e.
\newline
\newline {\fontfamily{qcr} \selectfont cd \$WM\_PROJECT\_USER\_DIR/my\_solvers} 
\newline {\fontfamily{qcr} \selectfont cp -r \$FOAM\_SOLVERS/incompressible/icoFoam icoTempFoam}
\newline
\newline Rename the {\fontfamily{qcr}\selectfont icoFoam.C} file to the {\fontfamily{qcr}\selectfont icoTempFoam.C} by using rename command.
\newline {\fontfamily{qcr} \selectfont cd   icoTempFoam}
\newline {\fontfamily{qcr} \selectfont mv icoFoam.C icoTempFoam.C}
\newline Keep the file {\fontfamily{qcr} \selectfont createFields.H} as it is.

\subsection{To modify the file {\fontfamily{qcr}\selectfont Make/files}}
Go to the {\fontfamily{qcr}\selectfont Make} directory and into the file {\fontfamily{qcr}\selectfont files} rename {\fontfamily{qcr}\selectfont icoFoam.C} to {\fontfamily{qcr}\selectfont icoTempFoam.C}. 
\newline Also change the {\fontfamily{qcr}\selectfont EXE=\$(FOAM\_APPBIN)/icoFoam}  to {\fontfamily{qcr}\selectfont EXE=\$(FOAM\_USER\_APPBIN)/icoTempFoam}  as follows:
\begin{verbatim}
    icoTempFoam.C
    EXE = $(FOAM_USER_APPBIN)/icoTempFoam
\end{verbatim}
Do not make any changes to the {\fontfamily{qcr}\selectfont Make/options}.

\subsection{To add {\fontfamily{qcr}\selectfont fvScalarMatrix TEqn} into the {\fontfamily{qcr}\selectfont icoTempFoam.C}}

The {\fontfamily{qcr}\selectfont icoTempFoam.C} file is modified by adding the {\fontfamily{qcr}\selectfont fvScalarMatrix TEqn} into it. The modified file {\fontfamily{qcr}\selectfont icoTempFoam.C} is as follows:
\begin{verbatim}
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "pisoControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Transient solver for incompressible, laminar flow"
        " of Newtonian fluids."
    );
    #include "postProcess.H"
    #include "addCheckCaseOptions.H"
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"
    pisoControl piso(mesh);
    #include "createFields.H"
    #include "initContinuityErrs.H"
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    Info<< "\nStarting time loop\n" << endl;
    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;
        #include "CourantNo.H"
    
    // Momentum predictor
        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );
    
    if (piso.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));
        }
    
    // --- PISO loop
        while (piso.correct())
        {
            volScalarField rAU(1.0/UEqn.A());
            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
      
      surfaceScalarField phiHbyA
            (
                "phiHbyA",
                fvc::flux(HbyA)
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );
      
      adjustPhi(phiHbyA, U, p);
      
      // Update the pressure BCs to ensure flux consistency
            constrainPressure(p, U, phiHbyA, rAU);
      
      // Non-orthogonal pressure corrector loop
            while (piso.correctNonOrthogonal())
            {
                // Pressure corrector
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );
                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
                if (piso.finalNonOrthogonalIter())
                {
                    phi = phiHbyA - pEqn.flux();
                }
            }
            #include "continuityErrs.H"
            U = HbyA - rAU*fvc::grad(p);
            U.correctBoundaryConditions();
        }       
        //Temperature equation 
        fvScalarMatrix TEqn
        (
              fvm::ddt(T)
            + fvm::div(phi, T)
            - fvm::laplacian(kappa, T)
        );
        {
            TEqn.solve();
            T.correctBoundaryConditions();
        }
        runTime.write();
        runTime.printExecutionTime(Info);
    }
    Info<< "End\n" << endl;
    return 0;
}
// ************************************************************************* //
\end{verbatim}

\subsection{To add dictionary reader for kappa ($\kappa$)  and IOobject for {\fontfamily{qcr}\selectfont volScalarField T} into the {\fontfamily{qcr}\selectfont createFields.H}}
To read the value of kappa ($\kappa$) from the {\fontfamily{qcr}\selectfont constant/transportProperties}, the entry of kappa is added below the entry of nu. Also, to read and write the vaue of {\fontfamily{qcr}\selectfont volScalarField T} the IOobject is written into the {\fontfamily{qcr}\selectfont createFields.H} as shown below:
\begin{verbatim}
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
);

dimensionedScalar kappa
(  
    "kappa",
    dimensionSet (0,2,-1,0,0,0,0),
    transportProperties
);
Info<< "Reading field p\n" << endl;
volScalarField p
(    IOobject
    (   
        "p",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(   IOobject
    (  
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);
Info<< "Reading field T\n" << endl;
volScalarField T
(   IOobject
    (   "T",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);
#include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());
\end{verbatim}
\subsection{Compilation of the newly modified {\fontfamily{qcr}\selectfont icoTempFoam} solver.}
The {\fontfamily{qcr}\selectfont icoTempFoam} solver is compiled by running the {\fontfamily{qcr}\selectfont wmake} command into the \\{\fontfamily{qcr}\selectfont \$WM\_PROJECT\_USER\_DIR/my\_solver/icoTempFoam} directory.
\\ The output message into the terminal after successful compilation of the solver is similar to the shown below. 

\subsection{To check the {\fontfamily{qcr}\selectfont icoTempFoam} solver in {\fontfamily{qcr}\selectfont \$FOAM\_USER\_APPBIN} }
The {\fontfamily{qcr}\selectfont \$FOAM\_USER\_APPBIN} directory is checked for the presence of newly compiled  solver is compiled {\fontfamily{qcr}\selectfont icoTempFoam} by the command {\fontfamily{qcr}\selectfont cd \$FOAM\_USER\_APPBIN} as shown below. 

\newpage
\section{Steps to modify the {\fontfamily{qcr}\selectfont icoFoam/cavity/cavity} tutorial case \& to test it for the {\fontfamily{qcr}\selectfont icoTempFoam} solver .}

\subsection{To copy the existing tutorial case to {\fontfamily{qcr}\selectfont \$FOAM\_RUN} directory and rename it to {\fontfamily{qcr}\selectfont cavityTemp}}

To copy the tutorial case from existing {\fontfamily{qcr}\selectfont icoFoam/cavity/cavity} case to {\fontfamily{qcr}\selectfont \$FOAM\_RUN} directory and rename it to {\fontfamily{qcr}\selectfont cavityTemp}. The following sets of command are used.
\newline
\newline {\fontfamily{qcr} \selectfont cd \$FOAM\_RUN} 
\newline {\fontfamily{qcr} \selectfont cp -r \$FOAM\_TUTORIALS/incompressible/icoFoam/cavity/cavity cavityTemp}
\newline
\newline Rename the {\fontfamily{qcr}\selectfont icoFoam.C} file to the {\fontfamily{qcr}\selectfont icoTempFoam.C} by using rename command.
\newline {\fontfamily{qcr} \selectfont cd   icoTempFoam}
\newline {\fontfamily{qcr} \selectfont mv icoFoam.C icoTempFoam.C}
\newline Keep the file {\fontfamily{qcr} \selectfont createFields.H} as it is.

\subsection{To add the {\fontfamily{qcr}\selectfont 0/T} file to specify the initial and boundary conditions for Temperature.}
The {\fontfamily{qcr}\selectfont 0/p} file is copied and renamed as {\fontfamily{qcr}\selectfont 0/T}. The object entry is changed to {\fontfamily{qcr}\selectfont T} while dimensions of temperature are set as  {\fontfamily{qcr}\selectfont [0 0 0 1 0 0 0]}. The {\fontfamily{qcr}\selectfont internalField} value is set to {\fontfamily{qcr}\selectfont uniform 300} and {\fontfamily{qcr}\selectfont boundaryField} type are changed to {\fontfamily{qcr}\selectfont fixedValue}. The details of the {\fontfamily{qcr}\selectfont 0/T} file are as given below.
\begin{verbatim}
\*---------------------------------------------------------------------------*/
FoamFile
{   version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions      [0 0 0 1 0 0 0];
internalField   uniform 300;
boundaryField
{
    movingWall
    {
        type            fixedValue;
        value           uniform 400;    
    }
    fixedWalls
    {
        type            fixedValue;
        value           uniform 300;
    }
    frontAndBack
    {
        type            empty;
    }
}
// ************************************************************************* //
\end{verbatim}
\newpage
\subsection {To add the {\fontfamily{qcr}\selectfont kappa} entry into the {\fontfamily{qcr}\selectfont constant/transportProperties}.}
The {\fontfamily{qcr}\selectfont constant/transportProperties} file is modified by adding the value of kappa as 0.001 . The details are as follows:
\begin{verbatim}
\*---------------------------------------------------------------------------*/
FoamFile
{   version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu              0.01;
kappa           [ 0 2 -1 0 0 0 0 ] 0.001;
// ************************************************************************* //
\end{verbatim}

\subsection{To change the solver name into the  {\fontfamily{qcr}\selectfont system/controlDict} file.}
In the {\fontfamily{qcr}\selectfont systrm/controlDict} file the solver name is changed to {\fontfamily{qcr}\selectfont icoTempFoam}. The rest of the details are kept as it is except {\fontfamily{qcr}\selectfont endTime} which is increased to 1.
\begin{verbatim}
\*---------------------------------------------------------------------------*/
FoamFile
{   version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application     icoTempFoam;
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         1;
deltaT          0.005;
writeControl    timeStep;
writeInterval   20;
purgeWrite      0;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;
// ************************************************************************* //
\end{verbatim}

\subsection {To add the divergence and laplacian scheme for temperature into the  {\fontfamily{qcr}\selectfont system/fvSchemes} file.}
To divergence  scheme for temperature i.e. {\fontfamily{qcr}\selectfont div(phi, T)   Gauss linear} is added into the list of divergence schemes whereas the laplacian scheme for temperature i.e. {\fontfamily{qcr}\selectfont laplacian(kappa, T)   Gauss linear orthogonal} is added into the list of laplacian schemes into the {\fontfamily{qcr}\selectfont system/fvSchemes} file.

\begin{verbatim}


\*---------------------------------------------------------------------------*/
FoamFile
{    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
    default         Euler;
}
gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
}
divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
    div(phi,T)      Gauss linear;
}
laplacianSchemes
{
    default         Gauss linear orthogonal;
    laplacian(kappa,T)         Gauss linear orthogonal;
}
interpolationSchemes
{
    default         linear;
}
snGradSchemes
{
    default         orthogonal;
}
// ************************************************************************* //
\end{verbatim}

\subsection{To add the solver settings for temperature in the {\fontfamily{qcr}\selectfont system/fvSolution} file }

The {\fontfamily{qcr}\selectfont smoothSolver} is set as solver for temperature into the {\fontfamily{qcr}\selectfont system/fvSolution} file as shown below.
\begin{verbatim}



\*---------------------------------------------------------------------------*/
FoamFile
{   version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
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;
    }
    T
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-07;
        relTol          0;
    }
}
PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}
// ************************************************************************* //
\end{verbatim}

\subsection{To test the solver {\fontfamily{qcr}\selectfont icoTempFoam} by running it on the above modified {\fontfamily{qcr}\selectfont cavityTemp} tutorial case. }
Steps to solve the {\fontfamily{qcr}\selectfont cavityTemp} tutorial case with {\fontfamily{qcr}\selectfont icoTempFoam} solver are as follows
\begin{itemize}
\item Run {\fontfamily{qcr}\selectfont blockMesh} command to generate mesh.
\item Run {\fontfamily{qcr}\selectfont checkMesh} command to check generated mesh
\item Run {\fontfamily{qcr}\selectfont icoTempFoam} to solve the tutorial case.
\item Post-process the results in the {\fontfamily{qcr}\selectfont ParaView} using appropriate command.
\end{itemize}
The terminal output for solution with {\fontfamily{qcr}\selectfont icoTempFoam} is as given below.

\newpage
\section{Results of the tutorial case {\fontfamily{qcr}\selectfont cavityTemp} for the test run of {\fontfamily{qcr}\selectfont icoTempFoam} solver}

\subsection{The initial \& final temperature contours} 

\begin{figure}[h]
    \centering
    \includegraphics[scale=3]{t ini.png.png}{ \label{fig:}}
    \caption{Initial temperature contours}
\end{figure}
\begin{figure}[h]
    \centering
    \includegraphics[scale=3]{t fin.png.png}{ \label{fig:}}
    \caption{Final temperature contours}
\end{figure}
\newpage
 \subsection{Observations \& Conclusion}
 \begin{enumerate}
\item  The Pressure and Velocity contours are similar to the results of  {\fontfamily{qcr}\selectfont icoFoam/cavity/cavity} tutorial case, as expected.
\item  The temperature contours show the sign of convection and diffusion as the temperature is more convected on the downstream side of the {\fontfamily{qcr}\selectfont movingWall} with respect to the upstream side. Also the temperature is slightly diffused to the sides and bottom.
\item Thus the we have successfully added the temperature equation to the {\fontfamily{qcr}\selectfont icoFoam} solver  and  successfully tested it on the test case {\fontfamily{qcr}\selectfont cavityTemp}. 
 \end{enumerate}
 
\end{document}

Results

temperature