Commit 24c6772a authored by mattijs's avatar mattijs
Browse files

Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

parents 1aea998b 3d039770
chemFoam.C
EXE = $(FOAM_APPBIN)/chemFoam
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude\
-I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
EXE_LIBS = \
-lfiniteVolume \
-lcompressibleRASModels \
-lreactionThermophysicalModels \
-lbasicThermophysicalModels \
-lchemistryModel \
-lODE \
-lthermophysicalFunctions \
-lspecie
{
forAll(Y, specieI)
{
volScalarField& Yi = Y[specieI];
solve
(
fvm::ddt(rho, Yi) - chemistry.RR(specieI),
mesh.solver("Yi")
);
}
}
\ No newline at end of file
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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
chemFoam
Description
Solver for chemistry problems
- designed for use on single cell cases to provide comparison against
other chemistry solvers
- single cell mesh created on-the-fly
- fields created on the fly from the initial conditions
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "hCombustionThermo.H"
#include "turbulenceModel.H"
#include "psiChemistryModel.H"
#include "chemistrySolver.H"
#include "OFstream.H"
#include "thermoPhysicsTypes.H"
#include "basicMultiComponentMixture.H"
#include "cellModeller.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createSingleCellMesh.H"
#include "createFields.H"
#include "readInitialConditions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readControls.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "solveChemistry.H"
{
#include "YEqn.H"
#include "hEqn.H"
#include "pEqn.H"
}
#include "output.H"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info << "Number of steps = " << runTime.timeIndex() << endl;
Info << "End" << nl << endl;
return(0);
}
// ************************************************************************* //
// write base thermo fields - not registered since will be re-read by
// thermo package
Info<< "Creating base fields for time " << runTime.timeName() << endl;
{
volScalarField Ydefault
(
IOobject
(
"Ydefault",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("Ydefault", dimless, 1)
);
Ydefault.write();
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("p", dimPressure, p0)
);
p.write();
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("T", dimTemperature, T0)
);
T.write();
}
if (mesh.nCells() != 1)
{
FatalErrorIn(args.executable())
<< "Solver only applicable to single cell cases"
<< exit(FatalError);
}
Info<< "Reading initial conditions.\n" << endl;
IOdictionary initialConditions
(
IOobject
(
"initialConditions",
runTime.constant(),
runTime,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
scalar p0 = readScalar(initialConditions.lookup("p"));
scalar T0 = readScalar(initialConditions.lookup("T"));
#include "createBaseFields.H"
Info<< nl << "Reading thermophysicalProperties" << endl;
autoPtr<psiChemistryModel> pChemistry(psiChemistryModel::New(mesh));
psiChemistryModel& chemistry = pChemistry();
scalar dtChem = refCast<const psiChemistryModel>(chemistry).deltaTChem()[0];
hsCombustionThermo& thermo = chemistry.thermo();
basicMultiComponentMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
runTime,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
thermo.rho()
);
volScalarField& p = thermo.p();
volScalarField& hs = thermo.hs();
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
runTime,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedVector("zero", dimVelocity, vector::zero),
p.boundaryField().types()
);
#include "createPhi.H"
Info << "Creating turbulence model.\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
OFstream post(args.path()/"chemFoam.out");
post<< "# Time" << token::TAB << "Temperature [K]" << token::TAB
<< "Pressure [Pa]" << endl;
Info<< "Constructing single cell mesh" << nl << endl;
labelList owner(6, 0);
labelList neighbour(0);
pointField points(8);
points[0] = vector(0, 0, 0);
points[1] = vector(1, 0, 0);
points[2] = vector(1, 1, 0);
points[3] = vector(0, 1, 0);
points[4] = vector(0, 0, 1);
points[5] = vector(1, 0, 1);
points[6] = vector(1, 1, 1);
points[7] = vector(0, 1, 1);
const cellModel& hexa = *(cellModeller::lookup("hex"));
faceList faces = hexa.modelFaces();
fvMesh mesh
(
IOobject
(
fvMesh::defaultRegion,
runTime.timeName(),
runTime,
IOobject::NO_READ
),
xferMove<Field<vector> >(points),
faces.xfer(),
owner.xfer(),
neighbour.xfer()
);
List<polyPatch*> patches(1);
patches[0] = new emptyPolyPatch("boundary", 6, 0, 0, mesh.boundaryMesh());
mesh.addFvPatches(patches);
{
if (constProp == "volume")
{
hs[0] = u0 + p[0]/rho[0] + integratedHeat;
}
else
{
hs[0] = hs0 + integratedHeat;
}
}
\ No newline at end of file
runTime.write();
Info<< "Sh = " << Sh
<< ", T = " << thermo.T()[0]
<< ", p = " << thermo.p()[0]
<< ", " << Y[0].name() << " = " << Y[0][0]
<< endl;
post<< runTime.value() << token::TAB << thermo.T()[0] << token::TAB
<< thermo.p()[0] << endl;
{
thermo.correct();
rho = thermo.rho();
if (constProp == "volume")
{
p[0] = rho0*R0*thermo.T()[0];
rho[0] = rho0;
}
}
\ No newline at end of file
if (runTime.controlDict().lookupOrDefault("suppressSolverInfo", false))
{
lduMatrix::debug = 0;
}
Switch adjustTimeStep(runTime.controlDict().lookup("adjustTimeStep"));
scalar maxDeltaT(readScalar(runTime.controlDict().lookup("maxDeltaT")));
word constProp(initialConditions.lookup("constantProperty"));
if (constProp == "pressure" || constProp == "volume")
{
Info << constProp << " will be held constant." << nl
<< " p = " << p[0] << " [Pa]" << nl
<< " T = " << thermo.T()[0] << " [K] " << nl
<< " rho = " << rho[0] << " [kg/m3]" << nl
<< endl;
}
else
{
FatalError << "in initialConditions, unknown constantProperty type "
<< constProp << nl << " Valid types are: pressure volume."
<< abort(FatalError);
}
word fractionBasis(initialConditions.lookup("fractionBasis"));
if ((fractionBasis != "mass") && (fractionBasis != "mole"))
{
FatalError << "in initialConditions, unknown fractionBasis type " << nl
<< "Valid types are: mass or mole."
<< fractionBasis << abort(FatalError);
}
label nSpecie = Y.size();
PtrList<gasThermoPhysics> specieData(Y.size());
forAll(specieData, i)
{
specieData.set
(
i,
new gasThermoPhysics
(
dynamic_cast<const reactingMixture<gasThermoPhysics>&>
(thermo).speciesData()[i]
)
);
}
scalarList Y0(nSpecie, 0.0);
scalarList X0(nSpecie, 0.0);
dictionary fractions(initialConditions.subDict("fractions"));
if (fractionBasis == "mole")
{
forAll(Y, i)
{
const word& name = Y[i].name();
if (fractions.found(name))
{
X0[i] = readScalar(fractions.lookup(name));
}
}
scalar mw = 0.0;
const scalar mTot = sum(X0);
forAll(Y, i)
{
X0[i] /= mTot;
mw += specieData[i].W()*X0[i];
}
forAll(Y, i)
{
Y0[i] = X0[i]*specieData[i].W()/mw;
}
}
else // mass fraction
{
forAll(Y, i)
{
const word& name = Y[i].name();
if (fractions.found(name))
{
Y0[i] = readScalar(fractions.lookup(name));
}
}
scalar invW = 0.0;
const scalar mTot = sum(Y0);
forAll(Y, i)
{
Y0[i] /= mTot;
invW += Y0[i]/specieData[i].W();
}
const scalar mw = 1.0/invW;
forAll(Y, i)
{
X0[i] = Y0[i]*mw/specieData[i].W();
}
}
scalar hs0 = 0.0;
forAll(Y, i)
{
Y[i] = Y0[i];
hs0 += Y0[i]*specieData[i].Hs(T0);
}
hs = dimensionedScalar("hs", dimEnergy/dimMass, hs0);
thermo.correct();
rho = thermo.rho();
scalar rho0 = rho[0];
scalar u0 = hs0 - p0/rho0;
scalar R0 = p0/(rho0*T0);
scalar integratedHeat = 0.0;
if (adjustTimeStep)
{
runTime.setDeltaT(min(dtChem, maxDeltaT));
Info<< "deltaT = " << runTime.deltaT().value() << endl;
}
dtChem = chemistry.solve
(
runTime.value() - runTime.deltaT().value(),
runTime.deltaT().value()
);
scalar Sh = chemistry.Sh()()[0]/rho[0];
integratedHeat += Sh*runTime.deltaT().value();
......@@ -118,8 +118,8 @@ void writeVTKFields
forAll(values, fieldI)
{
Info<< " writing field " << fieldNames[fieldI] << endl;
os << nl << fieldNames[fieldI] << ' ' << pTraits<Type>::nComponents << ' '
<< values[fieldI].size() << " float" << nl;
os << nl << fieldNames[fieldI] << ' ' << pTraits<Type>::nComponents
<< ' ' << values[fieldI].size() << " float" << nl;
label offset = 0;
forAll(agePerTrack, trackI)
{
......
......@@ -67,13 +67,15 @@ cleanCase()
rm -rf forces* > /dev/null 2>&1
rm -rf sets > /dev/null 2>&1
rm -rf system/machines > /dev/null 2>&1
(cd constant/polyMesh && \
rm -rf \
allOwner* cell* face* meshModifiers* \
owner* neighbour* point* edge* \
cellLevel* pointLevel* refinementHistory* surfaceIndex* sets \
> /dev/null 2>&1 \
)
if [ -d "constant/polyMesh" ]; then
(cd constant/polyMesh && \
rm -rf \
allOwner* cell* face* meshModifiers* \
owner* neighbour* point* edge* \
cellLevel* pointLevel* refinementHistory* surfaceIndex* sets \
> /dev/null 2>&1 \
)
fi
(cd constant && \
rm -rf \
cellToRegion cellLevel* pointLevel* \
......
......@@ -8,10 +8,10 @@
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 2 of the License, or (at your
option) any later version.
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
......@@ -19,8 +19,7 @@ License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
......
......@@ -8,10 +8,10 @@
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 2 of the License, or (at your
option) any later version.
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
......@@ -19,8 +19,7 @@ License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::cachedRandom
......
......@@ -8,10 +8,10 @@
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 2 of the License, or (at your
option) any later version.
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
......@@ -19,8 +19,7 @@ License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA