From 7ec7cd88e83697ea45be4d989470493f3f7684f6 Mon Sep 17 00:00:00 2001 From: Yeongdo Park Date: Tue, 30 Oct 2018 13:14:46 -0400 Subject: [PATCH] copy of chemFoam, base for test application --- .gitignore | 6 ++ testApp/Make/files | 3 + testApp/Make/options | 15 +++++ testApp/YEqn.H | 12 ++++ testApp/chemFoam.C | 93 +++++++++++++++++++++++++++ testApp/createBaseFields.H | 57 +++++++++++++++++ testApp/createControls.H | 2 + testApp/createFieldRefs.H | 4 ++ testApp/createFields.H | 84 +++++++++++++++++++++++++ testApp/createSingleCellMesh.H | 46 ++++++++++++++ testApp/hEqn.H | 14 +++++ testApp/output.H | 11 ++++ testApp/pEqn.H | 16 +++++ testApp/readControls.H | 3 + testApp/readInitialConditions.H | 108 ++++++++++++++++++++++++++++++++ testApp/setDeltaT.H | 6 ++ testApp/solveChemistry.H | 3 + 17 files changed, 483 insertions(+) create mode 100644 testApp/Make/files create mode 100644 testApp/Make/options create mode 100644 testApp/YEqn.H create mode 100644 testApp/chemFoam.C create mode 100644 testApp/createBaseFields.H create mode 100644 testApp/createControls.H create mode 100644 testApp/createFieldRefs.H create mode 100644 testApp/createFields.H create mode 100644 testApp/createSingleCellMesh.H create mode 100644 testApp/hEqn.H create mode 100644 testApp/output.H create mode 100644 testApp/pEqn.H create mode 100644 testApp/readControls.H create mode 100644 testApp/readInitialConditions.H create mode 100644 testApp/setDeltaT.H create mode 100644 testApp/solveChemistry.H diff --git a/.gitignore b/.gitignore index 1193b2d..23da90b 100644 --- a/.gitignore +++ b/.gitignore @@ -18,3 +18,9 @@ lnInclude /laminarReactingFoam/Make/* !/laminarReactingFoam/Make/file !/laminarReactingFoam/Make/option + +/testApp/Make/* +!/testApp/Make/file +!/testApp/Make/option + +canteraTest diff --git a/testApp/Make/files b/testApp/Make/files new file mode 100644 index 0000000..33e8507 --- /dev/null +++ b/testApp/Make/files @@ -0,0 +1,3 @@ +chemFoam.C + +EXE = canteraTest diff --git a/testApp/Make/options b/testApp/Make/options new file mode 100644 index 0000000..d08ab34 --- /dev/null +++ b/testApp/Make/options @@ -0,0 +1,15 @@ +EXE_INC = \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ + -I$(LIB_SRC)/ODE/lnInclude\ + -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lchemistryModel \ + -lfiniteVolume \ + -lmeshTools \ + -lcantera diff --git a/testApp/YEqn.H b/testApp/YEqn.H new file mode 100644 index 0000000..edb4072 --- /dev/null +++ b/testApp/YEqn.H @@ -0,0 +1,12 @@ +{ + forAll(Y, specieI) + { + volScalarField& Yi = Y[specieI]; + + solve + ( + fvm::ddt(rho, Yi) - chemistry.RR(specieI), + mesh.solver("Yi") + ); + } +} \ No newline at end of file diff --git a/testApp/chemFoam.C b/testApp/chemFoam.C new file mode 100644 index 0000000..a68d7d0 --- /dev/null +++ b/testApp/chemFoam.C @@ -0,0 +1,93 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\/ 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 . + +Application + chemFoam + +Description + Solver for chemistry problems, designed for use on single cell cases to + provide comparison against other chemistry solvers, that uses a single cell + mesh, and fields created from the initial conditions. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "psiReactionThermo.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[]) +{ + argList::noParallel(); + + #define CREATE_MESH createSingleCellMesh.H + #define NO_CONTROL + #include "postProcess.H" + + #include "setRootCase.H" + #include "createTime.H" + #include "createSingleCellMesh.H" + #include "createFields.H" + #include "createFieldRefs.H" + #include "readInitialConditions.H" + #include "createControls.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; +} + + +// ************************************************************************* // diff --git a/testApp/createBaseFields.H b/testApp/createBaseFields.H new file mode 100644 index 0000000..0762f27 --- /dev/null +++ b/testApp/createBaseFields.H @@ -0,0 +1,57 @@ +// 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(); +} + diff --git a/testApp/createControls.H b/testApp/createControls.H new file mode 100644 index 0000000..6aeb2ff --- /dev/null +++ b/testApp/createControls.H @@ -0,0 +1,2 @@ +Switch adjustTimeStep(runTime.controlDict().lookup("adjustTimeStep")); +scalar maxDeltaT(readScalar(runTime.controlDict().lookup("maxDeltaT"))); diff --git a/testApp/createFieldRefs.H b/testApp/createFieldRefs.H new file mode 100644 index 0000000..ef34180 --- /dev/null +++ b/testApp/createFieldRefs.H @@ -0,0 +1,4 @@ +scalar dtChem = refCast(chemistry).deltaTChem()[0]; +basicMultiComponentMixture& composition = thermo.composition(); +PtrList& Y = composition.Y(); +volScalarField& p = thermo.p(); diff --git a/testApp/createFields.H b/testApp/createFields.H new file mode 100644 index 0000000..62be3ce --- /dev/null +++ b/testApp/createFields.H @@ -0,0 +1,84 @@ + if (mesh.nCells() != 1) + { + FatalErrorInFunction + << "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 pChemistry(psiChemistryModel::New(mesh)); + + psiChemistryModel& chemistry = pChemistry(); + + psiReactionThermo& thermo = chemistry.thermo(); + thermo.validate(args.executable(), "h"); + + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + runTime, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + thermo.rho() + ); + + volScalarField Rspecific + ( + IOobject + ( + "Rspecific", + runTime.timeName(), + runTime, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar + ( + "zero", + dimensionSet(dimEnergy/dimMass/dimTemperature), + 0.0 + ) + ); + + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + runTime, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedVector("zero", dimVelocity, Zero) + ); + + #include "createPhi.H" + + OFstream post(args.path()/"chemFoam.out"); + post<< "# Time" << token::TAB << "Temperature [K]" << token::TAB + << "Pressure [Pa]" << endl; diff --git a/testApp/createSingleCellMesh.H b/testApp/createSingleCellMesh.H new file mode 100644 index 0000000..9b9d0f3 --- /dev/null +++ b/testApp/createSingleCellMesh.H @@ -0,0 +1,46 @@ +Info<< "Constructing single cell mesh" << nl << endl; + +labelList owner(6, label(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::READ_IF_PRESENT + ), + xferMove>(points), + faces.xfer(), + owner.xfer(), + neighbour.xfer() +); + +List patches(1); + +patches[0] = new emptyPolyPatch +( + "boundary", + 6, + 0, + 0, + mesh.boundaryMesh(), + emptyPolyPatch::typeName +); + +mesh.addFvPatches(patches); diff --git a/testApp/hEqn.H b/testApp/hEqn.H new file mode 100644 index 0000000..4dd593e --- /dev/null +++ b/testApp/hEqn.H @@ -0,0 +1,14 @@ +{ + volScalarField& h = thermo.he(); + + if (constProp == "volume") + { + h[0] = u0 + p[0]/rho[0] + integratedHeat; + } + else + { + h[0] = h0 + integratedHeat; + } + + thermo.correct(); +} diff --git a/testApp/output.H b/testApp/output.H new file mode 100644 index 0000000..5620fe6 --- /dev/null +++ b/testApp/output.H @@ -0,0 +1,11 @@ + 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; + diff --git a/testApp/pEqn.H b/testApp/pEqn.H new file mode 100644 index 0000000..a919ff1 --- /dev/null +++ b/testApp/pEqn.H @@ -0,0 +1,16 @@ +{ + rho = thermo.rho(); + if (constProp == "volume") + { + scalar invW = 0.0; + forAll(Y, i) + { + invW += Y[i][0]/specieData[i].W(); + } + + Rspecific[0] = 1000.0*constant::physicoChemical::R.value()*invW; + + p[0] = rho0*Rspecific[0]*thermo.T()[0]; + rho[0] = rho0; + } +} diff --git a/testApp/readControls.H b/testApp/readControls.H new file mode 100644 index 0000000..c56206d --- /dev/null +++ b/testApp/readControls.H @@ -0,0 +1,3 @@ +runTime.controlDict().lookup("adjustTimeStep") >> adjustTimeStep; + +maxDeltaT = readScalar(runTime.controlDict().lookup("maxDeltaT")); diff --git a/testApp/readInitialConditions.H b/testApp/readInitialConditions.H new file mode 100644 index 0000000..8e340a7 --- /dev/null +++ b/testApp/readInitialConditions.H @@ -0,0 +1,108 @@ + word constProp(initialConditions.lookup("constantProperty")); + if ((constProp != "pressure") && (constProp != "volume")) + { + 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 specieData(Y.size()); + forAll(specieData, i) + { + specieData.set + ( + i, + new gasHThermoPhysics + ( + dynamic_cast&> + (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 h0 = 0.0; + forAll(Y, i) + { + Y[i] = Y0[i]; + h0 += Y0[i]*specieData[i].Hs(p[0], T0); + } + + thermo.he() = dimensionedScalar("h", dimEnergy/dimMass, h0); + thermo.correct(); + + rho = thermo.rho(); + scalar rho0 = rho[0]; + scalar u0 = h0 - p0/rho0; + scalar R0 = p0/(rho0*T0); + Rspecific[0] = R0; + + scalar integratedHeat = 0.0; + + 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; diff --git a/testApp/setDeltaT.H b/testApp/setDeltaT.H new file mode 100644 index 0000000..46d9f7b --- /dev/null +++ b/testApp/setDeltaT.H @@ -0,0 +1,6 @@ +if (adjustTimeStep) +{ + runTime.setDeltaT(min(dtChem, maxDeltaT)); + Info<< "deltaT = " << runTime.deltaT().value() << endl; +} + diff --git a/testApp/solveChemistry.H b/testApp/solveChemistry.H new file mode 100644 index 0000000..cb4258c --- /dev/null +++ b/testApp/solveChemistry.H @@ -0,0 +1,3 @@ + dtChem = chemistry.solve(runTime.deltaT().value()); + scalar Sh = chemistry.Sh()()[0]/rho[0]; + integratedHeat += Sh*runTime.deltaT().value();