diff --git a/applications/solvers/combustion/fireFoam/YEEqn.H b/applications/solvers/combustion/fireFoam/YEEqn.H index b84d655a0d..1c9c2176c5 100644 --- a/applications/solvers/combustion/fireFoam/YEEqn.H +++ b/applications/solvers/combustion/fireFoam/YEEqn.H @@ -16,7 +16,7 @@ tmp> mvConvection forAll(Y, i) { - if (i != inertIndex) + if (i != inertIndex && composition.active(i)) { volScalarField& Yi = Y[i]; diff --git a/applications/solvers/combustion/reactingFoam/YEqn.H b/applications/solvers/combustion/reactingFoam/YEqn.H index 246f7119ab..7c3f0f117b 100644 --- a/applications/solvers/combustion/reactingFoam/YEqn.H +++ b/applications/solvers/combustion/reactingFoam/YEqn.H @@ -16,7 +16,7 @@ tmp> mvConvection forAll(Y, i) { - if (i != inertIndex) + if (i != inertIndex && composition.active(i)) { volScalarField& Yi = Y[i]; diff --git a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H index 52e44e2c47..4edbce9b42 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H @@ -17,7 +17,7 @@ tmp> mvConvection forAll(Y, i) { - if (i != inertIndex) + if (i != inertIndex && composition.active(i)) { volScalarField& Yi = Y[i]; diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H index 1fbb98445b..555b5bc0bc 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H @@ -17,7 +17,7 @@ tmp> mvConvection forAll(Y, i) { - if (i != inertIndex) + if (i != inertIndex && composition.active(i)) { volScalarField& Yi = Y[i]; diff --git a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H index 7259a1ca2a..db6628936c 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H @@ -16,7 +16,7 @@ tmp> mvConvection forAll(Y, i) { - if (i != inertIndex) + if (i != inertIndex && composition.active(i)) { volScalarField& Yi = Y[i]; diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H index 6498fd73f2..7c7cf9a4c0 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H @@ -16,7 +16,7 @@ tmp> mvConvection forAll(Y, i) { - if (i != inertIndex) + if (i != inertIndex && composition.active(i)) { volScalarField& Yi = Y[i]; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C index efa486a0e4..081e5c893a 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -142,6 +142,12 @@ Foam::MultiComponentPhaseModel::YiEqn this->name() ) ) + || ( + !this->thermo_->composition().active + ( + this->thermo_->composition().species()[Yi.member()] + ) + ) ) { return tmp(); diff --git a/bin/tools/CleanFunctions b/bin/tools/CleanFunctions index dde53bafc3..3df910567d 100644 --- a/bin/tools/CleanFunctions +++ b/bin/tools/CleanFunctions @@ -57,6 +57,7 @@ cleanCase() rm -rf processor* > /dev/null 2>&1 rm -rf postProcessing > /dev/null 2>&1 + rm -rf TDAC > /dev/null 2>&1 rm -rf probes* > /dev/null 2>&1 rm -rf forces* > /dev/null 2>&1 rm -rf graphs* > /dev/null 2>&1 diff --git a/src/thermophysicalModels/chemistryModel/Make/files b/src/thermophysicalModels/chemistryModel/Make/files index ef8c6d84c8..d0782c33b9 100644 --- a/src/thermophysicalModels/chemistryModel/Make/files +++ b/src/thermophysicalModels/chemistryModel/Make/files @@ -6,6 +6,9 @@ chemistryModel/psiChemistryModel/psiChemistryModels.C chemistryModel/rhoChemistryModel/rhoChemistryModel.C chemistryModel/rhoChemistryModel/rhoChemistryModels.C +chemistryModel/TDACChemistryModel/reduction/makeChemistryReductionMethods.C +chemistryModel/TDACChemistryModel/tabulation/makeChemistryTabulationMethods.C + chemistrySolver/chemistrySolver/makeChemistrySolvers.C LIB = $(FOAM_LIBBIN)/libchemistryModel diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C new file mode 100644 index 0000000000..de13a49d01 --- /dev/null +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C @@ -0,0 +1,843 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "TDACChemistryModel.H" +#include "UniformField.H" +#include "clockTime.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::TDACChemistryModel::TDACChemistryModel +( + const fvMesh& mesh, + const word& phaseName +) +: + chemistryModel(mesh, phaseName), + NsDAC_(this->nSpecie_), + completeC_(this->nSpecie_,0.0), + reactionsDisabled_(this->reactions_.size(), false), + completeToSimplifiedIndex_(this->nSpecie_,-1), + simplifiedToCompleteIndex_(this->nSpecie_), + specieComp_(this->nSpecie_) +{ + basicMultiComponentMixture& composition = this->thermo().composition(); + + // Store the species composition according to the species index + speciesTable speciesTab = composition.species(); + + const HashTable>& specComp = + dynamicCast&>(this->thermo()) + .specieComposition(); + + forAll(specieComp_,i) + { + specieComp_[i] = specComp[this->Y()[i].name()]; + } + + mechRed_ = chemistryReductionMethod::New + ( + *this, + *this + ); + + // When the mechanism reduction method is used, the 'active' flag for every + // species should be initialized (by default 'active' is true) + if (mechRed_->active()) + { + forAll(this->Y(), i) + { + IOobject header + ( + this->Y()[i].name(), + mesh.time().timeName(), + mesh, + IOobject::NO_READ + ); + + // Check if the species file is provided, if not set inactive + // and NO_WRITE + if (!header.headerOk()) + { + composition.setInactive(i); + this->Y()[i].writeOpt() = IOobject::NO_WRITE; + } + } + } + + tabulation_ = chemistryTabulationMethod::New + ( + *this, + *this + ); + + if (mechRed_->log()) + { + cpuReduceFile_ = logFile("cpu_reduce.out"); + nActiveSpeciesFile_ = logFile("nActiveSpecies.out"); + } + + if (tabulation_->log()) + { + cpuAddFile_ = logFile("cpu_add.out"); + cpuRetrieveFile_ = logFile("cpu_retrieve.out"); + } + + if (mechRed_->log() || tabulation_->log()) + { + cpuSolveFile_ = logFile("cpu_solve.out"); + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::TDACChemistryModel::~TDACChemistryModel() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::TDACChemistryModel::omega +( + const scalarField& c, // Contains all species even when mechRed is active + const scalar T, + const scalar p, + scalarField& dcdt +) const +{ + const bool reduced = mechRed_->active(); + + scalar pf, cf, pr, cr; + label lRef, rRef; + + dcdt = Zero; + + forAll(this->reactions_, i) + { + if (!reactionsDisabled_[i]) + { + const Reaction& R = this->reactions_[i]; + + scalar omegai = omega + ( + R, c, T, p, pf, cf, lRef, pr, cr, rRef + ); + + forAll(R.lhs(), s) + { + label si = R.lhs()[s].index; + if (reduced) + { + si = completeToSimplifiedIndex_[si]; + } + + const scalar sl = R.lhs()[s].stoichCoeff; + dcdt[si] -= sl*omegai; + } + forAll(R.rhs(), s) + { + label si = R.rhs()[s].index; + if (reduced) + { + si = completeToSimplifiedIndex_[si]; + } + + const scalar sr = R.rhs()[s].stoichCoeff; + dcdt[si] += sr*omegai; + } + } + } +} + + +template +Foam::scalar Foam::TDACChemistryModel::omega +( + const Reaction& R, + const scalarField& c, // Contains all species even when mechRed is active + const scalar T, + const scalar p, + scalar& pf, + scalar& cf, + label& lRef, + scalar& pr, + scalar& cr, + label& rRef +) const +{ + const scalar kf = R.kf(p, T, c); + const scalar kr = R.kr(kf, p, T, c); + + const label Nl = R.lhs().size(); + const label Nr = R.rhs().size(); + + label slRef = 0; + lRef = R.lhs()[slRef].index; + + pf = kf; + for (label s=1; s SMALL) + { + pf *= pow(cf, exp - 1); + } + else + { + pf = 0; + } + } + else + { + pf *= pow(cf, exp - 1); + } + } + + label srRef = 0; + rRef = R.rhs()[srRef].index; + + // Find the matrix element and element position for the rhs + pr = kr; + for (label s=1; sSMALL) + { + pr *= pow(cr, exp - 1); + } + else + { + pr = 0; + } + } + else + { + pr *= pow(cr, exp - 1); + } + } + + return pf*cf - pr*cr; +} + + +template +void Foam::TDACChemistryModel::derivatives +( + const scalar time, + const scalarField& c, + scalarField& dcdt +) const +{ + const bool reduced = mechRed_->active(); + + const scalar T = c[this->nSpecie_]; + const scalar p = c[this->nSpecie_ + 1]; + + if (reduced) + { + // When using DAC, the ODE solver submit a reduced set of species the + // complete set is used and only the species in the simplified mechanism + // are updated + this->c_ = completeC_; + + // Update the concentration of the species in the simplified mechanism + // the other species remain the same and are used only for third-body + // efficiencies + for (label i=0; ic_[simplifiedToCompleteIndex_[i]] = max(0.0, c[i]); + } + } + else + { + for (label i=0; inSpecie(); i++) + { + this->c_[i] = max(0.0, c[i]); + } + } + + omega(this->c_, T, p, dcdt); + + // Constant pressure + // dT/dt = ... + scalar rho = 0; + for (label i=0; ic_.size(); i++) + { + const scalar W = this->specieThermo_[i].W(); + rho += W*this->c_[i]; + } + + scalar cp = 0; + for (label i=0; ic_.size(); i++) + { + // cp function returns [J/(kmol K)] + cp += this->c_[i]*this->specieThermo_[i].cp(p, T); + } + cp /= rho; + + // When mechanism reduction is active + // dT is computed on the reduced set since dcdt is null + // for species not involved in the simplified mechanism + scalar dT = 0; + for (label i=0; inSpecie_; i++) + { + label si; + if (reduced) + { + si = simplifiedToCompleteIndex_[i]; + } + else + { + si = i; + } + + // ha function returns [J/kmol] + const scalar hi = this->specieThermo_[si].ha(p, T); + dT += hi*dcdt[i]; + } + dT /= rho*cp; + + dcdt[this->nSpecie_] = -dT; + + // dp/dt = ... + dcdt[this->nSpecie_ + 1] = 0; +} + + +template +void Foam::TDACChemistryModel::jacobian +( + const scalar t, + const scalarField& c, + scalarSquareMatrix& dfdc +) const +{ + const bool reduced = mechRed_->active(); + + // If the mechanism reduction is active, the computed Jacobian + // is compact (size of the reduced set of species) + // but according to the informations of the complete set + // (i.e. for the third-body efficiencies) + + const scalar T = c[this->nSpecie_]; + const scalar p = c[this->nSpecie_ + 1]; + + if (reduced) + { + this->c_ = completeC_; + for (label i=0; ic_[simplifiedToCompleteIndex_[i]] = max(0.0, c[i]); + } + } + else + { + forAll(this->c_, i) + { + this->c_[i] = max(c[i], 0.0); + } + } + + dfdc = Zero; + + forAll(this->reactions_, ri) + { + if (!reactionsDisabled_[ri]) + { + const Reaction& R = this->reactions_[ri]; + + const scalar kf0 = R.kf(p, T, this->c_); + const scalar kr0 = R.kr(kf0, p, T, this->c_); + + forAll(R.lhs(), j) + { + label sj = R.lhs()[j].index; + if (reduced) + { + sj = completeToSimplifiedIndex_[sj]; + } + scalar kf = kf0; + forAll(R.lhs(), i) + { + const label si = R.lhs()[i].index; + const scalar el = R.lhs()[i].exponent; + if (i == j) + { + if (el < 1) + { + if (this->c_[si] > SMALL) + { + kf *= el*pow(this->c_[si] + VSMALL, el - 1); + } + else + { + kf = 0; + } + } + else + { + kf *= el*pow(this->c_[si], el - 1); + } + } + else + { + kf *= pow(this->c_[si], el); + } + } + + forAll(R.lhs(), i) + { + label si = R.lhs()[i].index; + if (reduced) + { + si = completeToSimplifiedIndex_[si]; + } + const scalar sl = R.lhs()[i].stoichCoeff; + dfdc(si, sj) -= sl*kf; + } + forAll(R.rhs(), i) + { + label si = R.rhs()[i].index; + if (reduced) + { + si = completeToSimplifiedIndex_[si]; + } + const scalar sr = R.rhs()[i].stoichCoeff; + dfdc(si, sj) += sr*kf; + } + } + + forAll(R.rhs(), j) + { + label sj = R.rhs()[j].index; + if (reduced) + { + sj = completeToSimplifiedIndex_[sj]; + } + scalar kr = kr0; + forAll(R.rhs(), i) + { + const label si = R.rhs()[i].index; + const scalar er = R.rhs()[i].exponent; + if (i == j) + { + if (er < 1) + { + if (this->c_[si] > SMALL) + { + kr *= er*pow(this->c_[si] + VSMALL, er - 1); + } + else + { + kr = 0; + } + } + else + { + kr *= er*pow(this->c_[si], er - 1); + } + } + else + { + kr *= pow(this->c_[si], er); + } + } + + forAll(R.lhs(), i) + { + label si = R.lhs()[i].index; + if (reduced) + { + si = completeToSimplifiedIndex_[si]; + } + const scalar sl = R.lhs()[i].stoichCoeff; + dfdc(si, sj) += sl*kr; + } + forAll(R.rhs(), i) + { + label si = R.rhs()[i].index; + if (reduced) + { + si = completeToSimplifiedIndex_[si]; + } + const scalar sr = R.rhs()[i].stoichCoeff; + dfdc(si, sj) -= sr*kr; + } + } + } + } + + // Calculate the dcdT elements numerically + const scalar delta = 1e-3; + + omega(this->c_, T + delta, p, this->dcdt_); + for (label i=0; inSpecie_; i++) + { + dfdc(i, this->nSpecie_) = this->dcdt_[i]; + } + + omega(this->c_, T - delta, p, this->dcdt_); + for (label i=0; inSpecie_; i++) + { + dfdc(i, this->nSpecie_) = + 0.5*(dfdc(i, this->nSpecie_) - this->dcdt_[i])/delta; + } + + dfdc(this->nSpecie_, this->nSpecie_) = 0; + dfdc(this->nSpecie_ + 1, this->nSpecie_) = 0; +} + + +template +void Foam::TDACChemistryModel::jacobian +( + const scalar t, + const scalarField& c, + scalarField& dcdt, + scalarSquareMatrix& dfdc +) const +{ + jacobian(t, c, dfdc); + + const scalar T = c[this->nSpecie_]; + const scalar p = c[this->nSpecie_ + 1]; + + // Note: Uses the c_ field initialized by the call to jacobian above + omega(this->c_, T, p, dcdt); +} + + +template +template +Foam::scalar Foam::TDACChemistryModel::solve +( + const DeltaTType& deltaT +) +{ + const bool reduced = mechRed_->active(); + + basicMultiComponentMixture& composition = this->thermo().composition(); + + // CPU time analysis + const clockTime clockTime_= clockTime(); + clockTime_.timeIncrement(); + scalar reduceMechCpuTime_ = 0; + scalar addNewLeafCpuTime_ = 0; + scalar solveChemistryCpuTime_ = 0; + scalar searchISATCpuTime_ = 0; + + // Average number of active species + scalar nActiveSpecies = 0; + scalar nAvg = 0; + + CompType::correct(); + + scalar deltaTMin = GREAT; + + if (!this->chemistry_) + { + return deltaTMin; + } + + const volScalarField rho + ( + IOobject + ( + "rho", + this->time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + this->thermo().rho() + ); + + const scalarField& T = this->thermo().T(); + const scalarField& p = this->thermo().p(); + + scalarField c(this->nSpecie_); + scalarField c0(this->nSpecie_); + + // Composition vector (Yi, T, p) + scalarField phiq(this->nEqns()); + + scalarField Rphiq(this->nEqns()); + + forAll(rho, celli) + { + const scalar rhoi = rho[celli]; + scalar pi = p[celli]; + scalar Ti = T[celli]; + + for (label i=0; inSpecie_; i++) + { + c[i] = rhoi*this->Y_[i][celli]/this->specieThermo_[i].W(); + c0[i] = c[i]; + phiq[i] = this->Y()[i][celli]; + } + phiq[this->nSpecie()]=Ti; + phiq[this->nSpecie()+1]=pi; + + // Initialise time progress + scalar timeLeft = deltaT[celli]; + + // Not sure if this is necessary + Rphiq = Zero; + + clockTime_.timeIncrement(); + + // When tabulation is active (short-circuit evaluation for retrieve) + // It first tries to retrieve the solution of the system with the + // information stored through the tabulation method + if (tabulation_->active() && tabulation_->retrieve(phiq, Rphiq)) + { + // Retrieved solution stored in Rphiq + for (label i=0; inSpecie(); i++) + { + c[i] = rhoi*Rphiq[i]/this->specieThermo_[i].W(); + } + + searchISATCpuTime_ += clockTime_.timeIncrement(); + } + // This position is reached when tabulation is not used OR + // if the solution is not retrieved. + // In the latter case, it adds the information to the tabulation + // (it will either expand the current data or add a new stored poin). + else + { + clockTime_.timeIncrement(); + if (reduced) + { + // Reduce mechanism change the number of species (only active) + mechRed_->reduceMechanism(c, Ti, pi); + nActiveSpecies += mechRed_->NsSimp(); + nAvg++; + } + + reduceMechCpuTime_ += clockTime_.timeIncrement(); + + // Calculate the chemical source terms + while (timeLeft > SMALL) + { + scalar dt = timeLeft; + if (reduced) + { + // completeC_ used in the overridden ODE methods + // to update only the active species + completeC_ = c; + + // Solve the reduced set of ODE + this->solve + ( + simplifiedC_, Ti, pi, dt, this->deltaTChem_[celli] + ); + + for (label i=0; isolve(c, Ti, pi, dt, this->deltaTChem_[celli]); + } + timeLeft -= dt; + } + + solveChemistryCpuTime_ += clockTime_.timeIncrement(); + + // If tabulation is used, we add the information computed here to + // the stored points (either expand or add) + if (tabulation_->active()) + { + forAll(c, i) + { + Rphiq[i] = c[i]/rhoi*this->specieThermo_[i].W(); + } + + Rphiq[Rphiq.size()-2] = Ti; + Rphiq[Rphiq.size()-1] = pi; + tabulation_->add(phiq, Rphiq, rhoi); + } + + addNewLeafCpuTime_ += clockTime_.timeIncrement(); + + // When operations are done and if mechanism reduction is active, + // the number of species (which also affects nEqns) is set back + // to the total number of species (stored in the mechRed object) + if (reduced) + { + this->nSpecie_ = mechRed_->nSpecie(); + } + deltaTMin = min(this->deltaTChem_[celli], deltaTMin); + } + + // Set the RR vector (used in the solver) + for (label i=0; inSpecie_; i++) + { + this->RR_[i][celli] = + (c[i] - c0[i])*this->specieThermo_[i].W()/deltaT[celli]; + } + } + + if (mechRed_->log() || tabulation_->log()) + { + cpuSolveFile_() + << this->time().timeOutputValue() + << " " << solveChemistryCpuTime_ << endl; + } + + if (mechRed_->log()) + { + cpuReduceFile_() + << this->time().timeOutputValue() + << " " << reduceMechCpuTime_ << endl; + } + + if (tabulation_->active()) + { + // Every time-step, look if the tabulation should be updated + tabulation_->update(); + + // Write the performance of the tabulation + tabulation_->writePerformance(); + + if (tabulation_->log()) + { + cpuRetrieveFile_() + << this->time().timeOutputValue() + << " " << searchISATCpuTime_ << endl; + + cpuAddFile_() + << this->time().timeOutputValue() + << " " << addNewLeafCpuTime_ << endl; + } + } + + if (reduced && nAvg && mechRed_->log()) + { + // Write average number of species + nActiveSpeciesFile_() + << this->time().timeOutputValue() + << " " << nActiveSpecies/nAvg << endl; + } + + if (Pstream::parRun()) + { + List active(composition.active()); + Pstream::listCombineGather(active, orEqOp()); + Pstream::listCombineScatter(active); + + forAll(active, i) + { + if (active[i]) + { + composition.setActive(i); + } + } + } + + forAll(this->Y(), i) + { + if (composition.active(i)) + { + this->Y()[i].writeOpt() = IOobject::AUTO_WRITE; + } + } + + return deltaTMin; +} + + +template +Foam::scalar Foam::TDACChemistryModel::solve +( + const scalar deltaT +) +{ + // Don't allow the time-step to change more than a factor of 2 + return min + ( + this->solve>(UniformField(deltaT)), + 2*deltaT + ); +} + + +template +Foam::scalar Foam::TDACChemistryModel::solve +( + const scalarField& deltaT +) +{ + return this->solve(deltaT); +} + + +// ************************************************************************* // diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H new file mode 100644 index 0000000000..52a8f0860f --- /dev/null +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H @@ -0,0 +1,280 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Class + Foam::TDACChemistryModel + +Description + Extends chemistryModel by adding the TDAC method. + + References: + \verbatim + Contino, F., Jeanmart, H., Lucchini, T., & D’Errico, G. (2011). + Coupling of in situ adaptive tabulation and dynamic adaptive chemistry: + An effective method for solving combustion in engine simulations. + Proceedings of the Combustion Institute, 33(2), 3057-3064. + + Contino, F., Lucchini, T., D'Errico, G., Duynslaegher, C., + Dias, V., & Jeanmart, H. (2012). + Simulations of advanced combustion modes using detailed chemistry + combined with tabulation and mechanism reduction techniques. + SAE International Journal of Engines, + 5(2012-01-0145), 185-196. + + Contino, F., Foucher, F., Dagaut, P., Lucchini, T., D’Errico, G., & + Mounaïm-Rousselle, C. (2013). + Experimental and numerical analysis of nitric oxide effect on the + ignition of iso-octane in a single cylinder HCCI engine. + Combustion and Flame, 160(8), 1476-1483. + + Contino, F., Masurier, J. B., Foucher, F., Lucchini, T., D’Errico, G., & + Dagaut, P. (2014). + CFD simulations using the TDAC method to model iso-octane combustion + for a large range of ozone seeding and temperature conditions + in a single cylinder HCCI engine. + Fuel, 137, 179-184. + \endverbatim + +SourceFiles + TDACChemistryModelI.H + TDACChemistryModel.C + +\*---------------------------------------------------------------------------*/ + +#ifndef TDACChemistryModel_H +#define TDACChemistryModel_H + +#include "chemistryModel.H" +#include "chemistryReductionMethod.H" +#include "chemistryTabulationMethod.H" +#include "OFstream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class TDACChemistryModel Declaration +\*---------------------------------------------------------------------------*/ + +template +class TDACChemistryModel +: + public chemistryModel +{ + // Private member data + + // Mechanism reduction + label NsDAC_; + scalarField completeC_; + scalarField simplifiedC_; + Field reactionsDisabled_; + Field