Project 1- Creating your own openFOAM solver
Problem statement
The problem under consideration is to add an additional scalar transport equation (with out diffusion) which is covected through the velocity field.
`\frac{\partial s }{\partial t}+\grad.(mathbf u s)`
This requires us to modify the icofam solver to solve for the velocity and then solve for the scalar variable (s) using the above equation.
Introduction:
icoFoam is transient solver for incompressible and laminar flows. It solves the Navier-Stokes equation using the PISO(Pressure Iimplicit with Spliting of Operators) algorithm. The need of Pressure-Velocity Coupling alogrithms like SIMPLE and PISO arise from the fact that there is no Transport Equation for Pressure.
Procedure:
copy the icofoam solver from the location: openfoam/applications/solver/incompressible/icofoam to the run folder and rename it as scalarfoam. Similarlt rename the icofoam.C to scalarfoam.C
open the folder scalarfoam in the visual studio and the screen would like this
The following are the changes you need to add to the scalarfoam.C foam to add the transport equation:
Passive Scalar Transport Equation:
fvScalarMatrix sEqn
(
fvm::ddt(s) //Time derivative
+ fvm::div(phi, s) //convective term
- fvm::laplacian(su, s) //diffusion term
);
sEqn.relax(); //Considering under relaxing factor
sEqn.solve(); //Calling in Public member function SolveThe above part of the code creates an eqution with the type scalar matrix and the object name as sEqn. ddt is the time derivative term, div is for the divergence term and laplacian is for the diffusion term in the general equation. Once we create this equation we need to solve it so we use the object member function solve and relax to apply the relaxation factor while solving. it is important to note that all the terms here are treated implicitly.
The sequation sEqn needs to be declared in the while loop where the PISO solver is being solved, but after the momontum correction is performed. This is done and is shown in the code below:
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
icoFoam
Group
grpIncompressibleSolvers
Description
Transient solver for incompressible, laminar flow of Newtonian fluids.
\heading Solver details
The solver uses the PISO algorithm to solve the continuity equation:
\f[
\div \vec{U} = 0
\f]
and momentum equation:
\f[
\ddt{\vec{U}}
+ \div \left( \vec{U} \vec{U} \right)
- \div \left(\nu \grad \vec{U} \right)
= - \grad p
\f]
Where:
\vartable
\vec{U} | Velocity
p | Pressure
\endvartable
\heading Required fields
\plaintable
U | Velocity [m/s]
p | Kinematic pressure, p/rho [m2/s2]
\endplaintable
\*---------------------------------------------------------------------------*/
#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();
}
//Passive Scalar Transport Eqn
fvScalarMatrix sEqn
(
fvm::ddt(s)
+ fvm::div(phi, s)
- fvm::laplacian(su, s)
);
sEqn.relax(); //Considering under relaxing factor
sEqn.solve(); //Calling in Public member function Solve
runTime.write();
runTime.printExecutionTime(Info);
}
Info<< "End\n" << endl;
return 0;
}
Declaring the scalar variable 's' and the diffucion coefficient (su)
In order to use the variable 'su' and s i.e the scalar variable we need to create the variables at the runtime in the createFields.H file
The file looks like this:
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
);
#include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());
Info<< "Reading field s\n" << endl;
volScalarField s
(
IOobject
(
"s",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);we have declared a IO dictionary of the transportproperties and have read the values of su and nu from this dictionary object which is of type dimviscosity.
Similary inorder to set the value of s i.e. transport variable we have declared another object of type volume scalar field.
Next thing is to make the Make folder with the following corrections:
The file contains the following piece of code"
scalarfoam.C
EXE = $(FOAM_USER_APPBIN)/scalarfoam
Here: EXE is used to create an executable and it is stored with the name of scalarfoam at the user_app_bin location this is like an local directory where we can store all the user created executables which can then be accessed by the user.
The options file contains the following piece of information:
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools
In order to create the executable user need to type: wmake in the command terminal of the visual studio, as shown in the figure below:
The executable is created and now we are good to go to use the solver created.
Copy the test case named cavity from the OpenFoam Tutorials folder
$FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity/
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.
We will create an additional s file in the 0 folder which contains the information about 's' at various boundaries:
0/s File contents:
This file is a copy of pressure with specify changes such as:
dimension is zero since s is a scalar quantity
Internal field value is also uniformly zero.
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| / F ield | OpenFOAM: The Open Source CFD Toolbox |
| / O peration | Version: v2106 |
| / A nd | Website: www.openfoam.com |
| / M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
arch "LSB;label=32;scalar=64";
class volScalarField;
location "0";
object s;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
movingWall
{
type zeroGradient;
}
fixedWalls
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //0-P-File:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2112 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
movingWall
{
type zeroGradient;
}
fixedWalls
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //
0-U File
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2112 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
movingWall
{
type fixedValue;
value uniform (1 0 0);
}
fixedWalls
{
type noSlip;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //
System folder files:
Control dict:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
// ************************************************************************* //
fvSchemes:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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 upwind; // additionally added line to evaluate the div(Us) gausssian method and LinearUpwind Method // ******
}
laplacianSchemes
{
default Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //
fvsolution:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
}
// writing the fv solution i.e. for the scalar variable s which has to be transported by the velocity field U
s
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-08;
relTol 0;
minIter 1;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //
setfields:
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| / F ield | OpenFOAM: The Open Source CFD Toolbox |
| / O peration | Version: v2106 |
| / A nd | Website: www.openfoam.com |
| / M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object setFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defaultFieldValues
(
volScalarFieldValue s 0
);
regions
(
sphereToCell
{
centre (0.05 0.05 0);
radius 0.03;
fieldValues
(
volScalarFieldValue s 1
);
}
);
// ************************************************************************* //The set fields is used to create an patch of circle with radius 0.03 and centre at (X,Y)=(0.05 , 0.05) and make the value of the scalar field intially at t=0 to 1 as shown in the figure:
The above image is at t 0.01 second which shows the variation of the scalar field 's'. There is a patch where the value of s is defined equal to 1 and then its getting advected due to the zero diffusion in the transport properties.
Zero folder:
transportproperties file:
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| / F ield | OpenFOAM: The Open Source CFD Toolbox |
| / O peration | Version: v2106 |
| / A nd | Website: www.openfoam.com |
| / M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.001;
su [0 2 -1 0 0 0 0] 0.001;
// ************************************************************************* //
here i have declared a small value of diffusion coefficient to see the diffusion of the scalar 's' into the domain with the advection due to the the velocity field.
The video (Video 1)is attached along with the report to visulalize the effect of the diffusion.
Another test case was run by changeing the value of diffusion coefficient 'su' in the transport properties to zero.
constant -transportroperties
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| / F ield | OpenFOAM: The Open Source CFD Toolbox |
| / O peration | Version: v2106 |
| / A nd | Website: www.openfoam.com |
| / M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.001;
su [0 2 -1 0 0 0 0] 0;
// ************************************************************************* //
The output after the simulation is as follows:
We can see that there is no diffusion and the scalar variable is getting convected at the very first time step. The animated video of the pure advection is attached with the file name video 2.