Project 1- Creating your own openFOAM solver
Objective
Create an OpenFoam solver from an existing solver icoFoam by adding passive scalar transport equation.
Check the compiled solver to test cavity case and analysis the results.
Theory -
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 and rename it to scalarFoam :
~/OpenFOAM-v2106/applications/solvers/incompressible/icoFoam scalarFoam
Followed by appropriate changes to the name of the .C file
scalar2Foam.C
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 SolveHere we solve for scalar term s. The LHS of sEqn is stored in an fvScalarMatrix, inside which all the terms (time derivative, convection and diffusion) are treated implicitly. Following this, we apply relaxation factor based on user input and finally solve it(equate to 0).
On LHS side,
First term is a time derivative of scalar s.
Second term is a linearized convection term i.e. volume intergral of divergence of the vector field is converted into to surface integral of that vector, which is flux phi.
Last term is diffusion term with diffusion cofficient su and the diffusion term is subtracted from the LHS.
scalar2Foam.C - This file contains the source code of the solver
/*---------------------------------------------------------------------------*
========= |
\ / 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
scalarFoam
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
Adding passive Scalar Transport equation
f[
ddt{scalar{s}}
+ div left( vec{phi} scalar{s} right)
- div left(scalar su grad scalar{s} right)
= - grad s
f]
*---------------------------------------------------------------------------*/
#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 loopn" << 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<< "Endn" << endl;
return 0;
}
// ************************************************************************* //createFields.H - This file contains all IOObjects.
Here we add an additional volScalarField s and dimensionedScalar su
nfo<< "Reading transportPropertiesn" << 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 pn" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field Un" << 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 sn" << endl;
volScalarField s
(
IOobject
(
"s",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Make/Files -
scalar2Foam.C //Souce File
EXE = $(FOAM_USER_APPBIN)/scalar2Foam //Particular Path to store the excutable and Exceutable NameMake/Options
EXE_INC =
-I$(LIB_SRC)/finiteVolume/lnInclude
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS =
-lfiniteVolume
-lmeshToolsNow run the wmake command to compile the solver.
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.
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
`(∂s)/(∂t)+∇(Us)=∇^2(Ds)`
setFieldsDict
/*--------------------------------*- 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
);
}
);
// ************************************************************************* //
0/s - 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;
}
}
// ************************************************************************* //
consatnt/transportProperties
/*--------------------------------*- 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.01;
su [0 2 -1 0 0 0 0] 0.01;
// ************************************************************************* //system/fvSchemes
/*--------------------------------*- 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 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;
}
laplacianSchemes
{
default Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //system/fvSolutions
Here we must specify particular Linear Solver that is to be used for scalar matrix.Moreover, for scalar s we cannot use PCG as its already been specified for the pressure term.
PCG is meant for symmetric matrices.
PBiCG is for asymmetric matrices.
smoothSolver for any of them.
/*--------------------------------*- 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 fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
s
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-08;
relTol 0;
minIter 1;
}
pFinal
{
$p;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-05;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //Post-Processing
blockMesh - Creates the mesh as specified in blockMeshDict.
checkMesh - A utility which analyzes and evaluates mesh statistics and quality parameters.
setFields - It sets the field as specified in setFieldsDict.
scalarFoam - Executable to run the simulation case.