Project 1- Creating your own openFOAM solver
Adding Scalar Transport Equation in IcoFoam solver
1. Introduction
The aim of this project is to modify the existing IcoFoam solver to solve a scalar transport equation as given below.
dsdt+∇.(us)=0">dsdt+∇.(us)=0dsdt+∇.(us)=0 ....................................eqn(i)
In this case the scalar transport equation doesn't consist of any diffusive term and source term.
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) IcoFoam.C was renamed as scalarFoam.C. AlsoThe code given below was added to the existing IcoFoam solver file IcoFoam.C. The following code consist of temporal term and convective term only. The code was added after the PISO loop in the file.
#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();
}
fvScalarMatrix sEquations
(
fvm::ddt(s)
+ fvm:: div(phi,s)
- fvm:: laplacian(su,s)
);
sEquations.relax();
sEquations.solve();
runTime.write();
runTime.printExecutionTime(Info);
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
ii) In createFields.H file the scalar quantity S 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 S(scalar quantity).
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 su
(
"su",
dimViscosity,
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 s\n" << endl;
volScalarField s
(
IOobject
(
"s",
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());
iii) the changes were made in the Make>>files directory as given below.
scalarFoam.C
EXE = $(FOAM_USER_APPBIN)/scalarFoam$FOAM_USER_APPBIN saves the executable in the user defined directory.
With the above changes, the solver was compiled using the command wmake and it was successfully compiled.
3. Changes in the Tutorial
$FOAM_TUTORIALS/basic/scalarTransportFoam/pitzDaily case was used to test out our solver and renamed it to scalarpitzDaily.
i) In the folder 0, /as copied p folder was replaced by the scalar S as a dimensionless quantity.
ii) Since, no diffusive term was there to solve, the diffusive coefficient was also removed from the constant folder.
iii) In the fvSchemes folder the interpolation scheme for convective term ∇ · (us) was selected as linearUpwind and the gradient scheme was linear (Gauss cell based) as given below.
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
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,s) Gauss linear;
}
laplacianSchemes
{
default Gauss linear orthogonal;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //iv) In fvSolutions the solver for p was used by make a copy of p with scalar quantity S. As both quantity are scalar, hence same solver was used for solving S.
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
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;
}
s
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-06;
relTol 0.05;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //v) A setFieldsDict file was also used to visually enclose the scalar quantity inside a small circle using the following code below.
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object setFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defaultFieldValues
(
volScalarFieldValue s 0
);
regions
(
sphereToCell
{
centre (0.05 0.05 0);
radius 0.03;
fieldValues
(
volScalarFieldValue s 1
);
}
);
// ************************************************************************* //
vii) The changes were also made in the controlDict file. The solver name was changed and kept as scalarFoam for application. Other parameters were kept same as shown below.
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application scalarFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 1;
deltaT 1e-2;
writeControl timeStep;
writeInterval 1;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
viii) The case was run by the using the following commands:
- blockMesh ( to generate the mesh)
- setFields (to generate the scalar quantity enclosed in a circle/box)
- scalarFoam (to run the solver)
4. Animation of ScalarTransport simulation
After the case was run, Though diffusion term wasn't present in the eqn(i) but due to the numerical diffusion present in the discretization scheme, the scalar quantity can been seen getting diffused later in the simulation.
i uploaded snap shots because getting some error in saving animations
error during deffussion
before that there is no diffussion present in the simulation then i slightly changed reduced writing interval time in controldict then i got deffussion