diff --git a/applications/solvers/dfLowMachFoam/createFields.H b/applications/solvers/dfLowMachFoam/createFields.H index 5c4ccc7e..e4c5ba9b 100644 --- a/applications/solvers/dfLowMachFoam/createFields.H +++ b/applications/solvers/dfLowMachFoam/createFields.H @@ -22,6 +22,7 @@ volScalarField rho thermo.rho() ); + Info<< "Reading field U\n" << endl; volVectorField U ( @@ -54,12 +55,34 @@ autoPtr turbulence ) ); +Info<< "Creating field dpdt\n" << endl; +volScalarField dpdt +( + IOobject + ( + "dpdt", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("dpdt",p.dimensions()/dimTime, 0) +); + + Info<< "Creating reaction model\n" << endl; autoPtr> combustion ( CombustionModel::New(thermo, turbulence()) ); Info<< "end Creating reaction model\n" << endl; + + +//const word combModelName(combustion->type()); +const word combModelName(mesh.objectRegistry::lookupObject("combustionProperties").lookup("combustionModel")); +Info << "Combustion Model Name is confirmed as "<< combModelName << endl; + dfChemistryModel* chemistry = combustion->chemistry(); PtrList& Y = chemistry->Y(); const word inertSpecie(chemistry->lookup("inertSpecie")); @@ -69,28 +92,22 @@ const label inertIndex(chemistry->species()[inertSpecie]); chemistry->correctThermo(); Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl; -Info<< "Creating field dpdt\n" << endl; -volScalarField dpdt -( - IOobject - ( - "dpdt", - runTime.timeName(), - mesh - ), - mesh, - dimensionedScalar(p.dimensions()/dimTime, 0) -); +//for dpdt Info<< "Creating field kinetic energy K\n" << endl; volScalarField K("K", 0.5*magSqr(U)); multivariateSurfaceInterpolationScheme::fieldTable fields; + +if(combModelName!="flareFGM") +{ forAll(Y, i) { fields.add(Y[i]); } fields.add(thermo.he()); +} + const scalar Sct = chemistry->lookupOrDefault("Sct", 1.); volScalarField diffAlphaD @@ -152,4 +169,4 @@ const Switch splitting = CanteraTorchProperties.lookupOrDefault("splittingStrate #ifdef USE_LIBTORCH const Switch log_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("log", false); const Switch torch_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("torch", false); -#endif \ No newline at end of file +#endif diff --git a/applications/solvers/dfLowMachFoam/dfLowMachFoam.C b/applications/solvers/dfLowMachFoam/dfLowMachFoam.C index ec90b78e..437866bb 100644 --- a/applications/solvers/dfLowMachFoam/dfLowMachFoam.C +++ b/applications/solvers/dfLowMachFoam/dfLowMachFoam.C @@ -148,17 +148,24 @@ int main(int argc, char *argv[]) end = std::clock(); time_monitor_flow += double(end - start) / double(CLOCKS_PER_SEC); - #include "YEqn.H" + if(combModelName!="ESF" && combModelName!="flareFGM" ) + { + #include "YEqn.H" - start = std::clock(); - #include "EEqn.H" - end = std::clock(); - time_monitor_E += double(end - start) / double(CLOCKS_PER_SEC); + start = std::clock(); + #include "EEqn.H" + end = std::clock(); + time_monitor_E += double(end - start) / double(CLOCKS_PER_SEC); - start = std::clock(); - chemistry->correctThermo(); - end = std::clock(); - time_monitor_corrThermo += double(end - start) / double(CLOCKS_PER_SEC); + start = std::clock(); + chemistry->correctThermo(); + end = std::clock(); + time_monitor_corrThermo += double(end - start) / double(CLOCKS_PER_SEC); + } + else + { + combustion->correct(); + } Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl; @@ -189,6 +196,8 @@ int main(int argc, char *argv[]) runTime.write(); + Info << "output time index " << runTime.timeIndex() << endl; + Info<< "========Time Spent in diffenet parts========"<< endl; Info<< "Chemical sources = " << time_monitor_chem << " s" << endl; Info<< "Species Equations = " << time_monitor_Y << " s" << endl; diff --git a/src/dfCombustionModels/FGM/baseFGM/baseFGM.C b/src/dfCombustionModels/FGM/baseFGM/baseFGM.C new file mode 100644 index 00000000..670ac0b0 --- /dev/null +++ b/src/dfCombustionModels/FGM/baseFGM/baseFGM.C @@ -0,0 +1,610 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2018 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 "baseFGM.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::combustionModels::baseFGM::baseFGM +( + const word& modelType, + ReactionThermo& thermo, + const compressibleTurbulenceModel& turb, + const word& combustionProperties +) +: + laminar(modelType, thermo, turb, combustionProperties), + buffer_(this->coeffs().lookupOrDefault("buffer", false)), + scaledPV_(this->coeffs().lookupOrDefault("scaledPV", false)), + incompPref_(this->coeffs().lookupOrDefault("incompPref", -10.0)), + ignition_(this->coeffs().lookupOrDefault("ignition", false)), + combustion_(this->coeffs().lookupOrDefault("combustion", false)), + solveEnthalpy_(this->coeffs().lookupOrDefault("solveEnthalpy", false)), + flameletT_(this->coeffs().lookupOrDefault("flameletT", false)), + psi_(const_cast(dynamic_cast(thermo).psi())), + Wt_ + ( + IOobject + ( + "Wt", + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("Wt",dimensionSet(1,0,0,0,-1,0,0),28.96) + ), + Cp_ + ( + IOobject + ( + "Cp", + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("Cp",dimensionSet(0,2,-2,-1,0,0,0),1010.1) + ), + Z_ + ( + IOobject + ( + "Z", + this->mesh().time().timeName(), + this->mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh() + ), + Zvar_ + ( + IOobject + ( + "Zvar", + this->mesh().time().timeName(), + this->mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh() + ), + He_ + ( + IOobject + ( + "He", + this->mesh().time().timeName(), + this->mesh(), + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("He",dimensionSet(0,2,-2,0,0,0,0),0.0) + ), + Hf_ + ( + IOobject + ( + "Hf", + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("Hf",dimensionSet(0,2,-2,0,0,0,0),1907.0) + ), + c_ + ( + IOobject + ( + "c", + this->mesh().time().timeName(), + this->mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh() + ), + cvar_ + ( + IOobject + ( + "cvar", + this->mesh().time().timeName(), + this->mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh() + ), + Zcvar_ + ( + IOobject + ( + "Zcvar", + this->mesh().time().timeName(), + this->mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh() + ), + omega_c_ + ( + IOobject + ( + "omega_c", + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("omega_c",dimensionSet(1,-3,-1,0,0,0,0),0.0) + ), + chi_c_ + ( + IOobject + ( + "chi_c", + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("chi_c",dimensionSet(0,0,-1,0,0,0,0),0.0) + ), + chi_Z_ + ( + IOobject + ( + "chi_Z", + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("chi_Z",dimensionSet(0,0,-1,0,0,0,0),0.0) + ), + chi_Zc_ + ( + IOobject + ( + "chi_Zc", + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("chi_Zc",dimensionSet(0,0,-1,0,0,0,0),0.0) + ), + YCO2_ + ( + IOobject + ( + "YCO2", + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("YCO2",dimless,0.0) + ), + YCO_ + ( + IOobject + ( + "YCO", + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("YCO",dimless,0.0) + ), + YH2O_ + ( + IOobject + ( + "YH2O", + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar("YH2O",dimless,0.0) + ), + speciesNames_(this->coeffs().lookup("speciesName")), + Y_(speciesNames_.size()), + cOmega_c_(omega_c_), + ZOmega_c_(omega_c_), + WtCells_ (Wt_.primitiveFieldRef()), + CpCells_ (Cp_.primitiveFieldRef()), + ZCells_(Z_.primitiveFieldRef()), + ZvarCells_(Zvar_.primitiveFieldRef()), + HCells_(He_.primitiveFieldRef()), + HfCells_(Hf_.primitiveFieldRef()), + cCells_(c_.primitiveFieldRef()), + cvarCells_(cvar_.primitiveFieldRef()), + ZcvarCells_(Zcvar_.primitiveFieldRef()), + omega_cCells_(omega_c_.primitiveFieldRef()), + cOmega_cCells_(cOmega_c_.primitiveFieldRef()), + ZOmega_cCells_(ZOmega_c_.primitiveFieldRef()), + chi_cCells_(chi_c_.primitiveFieldRef()), + chi_ZCells_(chi_Z_.primitiveFieldRef()), + chi_ZcCells_(chi_Zc_.primitiveFieldRef()), + YCO2Cells_(YCO2_.primitiveFieldRef()), + YCOCells_(YCO_.primitiveFieldRef()), + YH2OCells_(YH2O_.primitiveFieldRef()), + ZMax_(1.0), + ZMin_(0.0), + ZvarMax_(0.25), + ZvarMin_(0.0), + cMax_(1.0), + cMin_(0.0), + Ycmaxall_(1.0), + cvarMax_(0.25), + cvarMin_(0.25), + ZcvarMax_(0.25), + ZcvarMin_(-0.25), + rho_(const_cast(this->mesh().objectRegistry::lookupObject("rho"))), + p_(this->thermo().p()), + T_(this->thermo().T()), + U_(this->mesh().objectRegistry::lookupObject("U")), + dpdt_(this->mesh().objectRegistry::lookupObject("dpdt")), + phi_(this->mesh().objectRegistry::lookupObject("phi")), + TCells_(T_.primitiveFieldRef()), + ignBeginTime_(this->coeffs().lookupOrDefault("ignBeginTime", GREAT)), + ignDurationTime_(this->coeffs().lookupOrDefault("ignBeginTime", 0.0)), + reactFlowTime_(0.0), + x0_(this->coeffs().lookupOrDefault("x0", 0.0)), + y0_(this->coeffs().lookupOrDefault("y0", 0.0)), + z0_(this->coeffs().lookupOrDefault("z0", 0.05)), + R0_(this->coeffs().lookupOrDefault("R0", 0.04)), + Sct_(this->coeffs().lookupOrDefault("Sct", 0.7)), + Sc_(this->coeffs().lookupOrDefault("Sc", 1.0)), + bufferTime_(this->coeffs().lookupOrDefault("bufferTime", 0.0)), + relaxation_(this->coeffs().lookupOrDefault("relaxation", false)), + DpDt_(this->coeffs().lookupOrDefault("DpDt", false)) +{ + if(incompPref_ < 0.0) //local pressure used to calculate EOS + { + Info<< "Equation of State: local pressure used" << nl << endl; + } + else //constant pressure used to calculate EOS + { + Info<< "Equation of State: constant pressure used: " + << incompPref_ << " Pa" << nl << endl; + } + + //initialize species fields + forAll(Y_, speciesI) + { + + Y_.set + ( + speciesI, + new volScalarField + ( + IOobject + ( + "Y_" + speciesNames_[speciesI], + this->mesh().time().timeName(), + this->mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + this->mesh(), + dimensionedScalar(dimless, Zero) + ) + ); + + } + + // add fields + fields_.add(Z_); + fields_.add(Zvar_); + fields_.add(c_); + fields_.add(cvar_); + fields_.add(Zcvar_); + fields_.add(He_); + +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::combustionModels::baseFGM::~baseFGM() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +template +void Foam::combustionModels::baseFGM::transport() +{ + buffer(); + + tmp tmut(this->turbulence().mut()); + const volScalarField& mut = tmut(); + + tmp tmu(this->turbulence().mu()); + const volScalarField& mu = tmu(); + + //scalarUW used for cEqn, cvarEqn, ZEqn, ZvarEqn,ZcvarEqn to ensure the convergence when buffer_ is true + tmp > scalarUWConvection + ( + fv::convectionScheme::New + ( + this->mesh(), + fields_, + phi_, + this->mesh().divScheme("div(phi,scalarUW)") + ) + ); + + + // Solve the mixture fraction transport equation + if(buffer_) + { + Info<< "UW schemes used for scalars" << endl; + + fvScalarMatrix ZEqn + ( + fvm::ddt(rho_,Z_) + +scalarUWConvection->fvmDiv(phi_, Z_) + -fvm::laplacian( mut/Sct_ + mu/Sc_ , Z_) + ); + if(relaxation_) + { + ZEqn.relax(); + } + ZEqn.solve(); + } + + else + { + Info<< "TVD schemes used for scalars" << endl; + + fvScalarMatrix ZEqn + ( + fvm::ddt(rho_,Z_) + +fvm::div(phi_, Z_) + -fvm::laplacian(mut/Sct_ + mu/Sc_ , Z_) + ); + + if(relaxation_) + { + ZEqn.relax(); + } + ZEqn.solve(); + } + Z_.min(ZMax_); + Z_.max(ZMin_); + + // solve the total enthalpy transport equation + if(solveEnthalpy_) + { + fvScalarMatrix HEqn + ( + fvm::ddt(rho_,He_) + +( + buffer_ + ? scalarUWConvection->fvmDiv(phi_, He_) + : fvm::div(phi_, He_) + ) + +( + DpDt_ + ? - dpdt_ - ( U_ & fvc::grad(p_) ) - fvm::laplacian( mut/Sct_ + mu/Sc_, He_) + : - fvm::laplacian( mut/Sct_ + mu/Sc_, He_) + ) + ); + + if(relaxation_) + { + HEqn.relax(); + } + HEqn.solve(); + } + + + // Solve the mixture fraction variance transport equation + + fvScalarMatrix ZvarEqn + ( + fvm::ddt(rho_,Zvar_) + +( + buffer_ + ? scalarUWConvection->fvmDiv(phi_, Zvar_) + : fvm::div(phi_, Zvar_) + ) + -fvm::laplacian( mut/Sct_+mu/Sc_, Zvar_) + -(2.0*mut/Sct_*(fvc::grad(Z_) & fvc::grad(Z_))) + +(2.0*rho_*chi_Z_) + ); + + if(relaxation_) + { + ZvarEqn.relax(); + } + + ZvarEqn.solve(); + + Zvar_.min(ZvarMax_); + Zvar_.max(ZvarMin_); + + + if(combustion_ && reactFlowTime_ > 0.0) + { + // At initial time, cMax is set as Ycmaxall when unscaled PV employed + if(this->mesh().time().timeIndex() == 1 && !scaledPV_) cMax_ = Ycmaxall_; + + // Solve the progress variable transport equation + fvScalarMatrix cEqn + ( + fvm::ddt(rho_, c_) + +( + buffer_ + ? scalarUWConvection->fvmDiv(phi_, c_) + : fvm::div(phi_, c_) + ) + -fvm::laplacian( mut/Sct_ + mu/Sc_, c_) + -omega_c_ + ); + if(relaxation_) + { + cEqn.relax(); + } + cEqn.solve(); + c_.min(cMax_); + c_.max(cMin_); + + // Solve the progress variable variance transport equation + fvScalarMatrix cvarEqn + ( + fvm::ddt(rho_,cvar_) + +( + buffer_ + ? scalarUWConvection->fvmDiv(phi_, cvar_) + : fvm::div(phi_, cvar_) + ) + -fvm::laplacian( mut/Sct_ + mu/Sc_, cvar_) + -(2.0*mut/Sct_*(fvc::grad(c_) & fvc::grad(c_))) + +2.0*(rho_*chi_c_) + -2.0*(cOmega_c_-omega_c_*c_) + ); + + if(relaxation_) + { + cvarEqn.relax(); + } + cvarEqn.solve(); + cvar_.min(cvarMax_); + cvar_.max(cvarMin_); + + + // Solve the covariance transport equation + fvScalarMatrix ZcvarEqn + ( + fvm::ddt(rho_,Zcvar_) + +( + buffer_ + ? scalarUWConvection->fvmDiv(phi_, Zcvar_) + : fvm::div(phi_, Zcvar_) + ) + -fvm::laplacian( mut/Sct_+mu/Sc_, Zcvar_) + -(1.0*mut/Sct_*(fvc::grad(c_) & fvc::grad(Z_))) + -(1.0*mut/Sct_*(fvc::grad(Z_) & fvc::grad(c_))) + +(2.0*rho_*chi_Zc_) + -1.0*(ZOmega_c_-omega_c_*Z_) + ); + + if(relaxation_) + { + ZcvarEqn.relax(); + } + + ZcvarEqn.solve(); + + Zcvar_.min(ZcvarMax_); + Zcvar_.max(ZcvarMin_); + + } + + +} + +template +void Foam::combustionModels::baseFGM::initialiseFalmeKernel() +{ + + reactFlowTime_= this->mesh().time().value()-ignBeginTime_; + + Info<< "Time = " << this->mesh().time().timeName() + << " reactFlowTime = " << reactFlowTime_ << nl << endl; + + if(ignition_ && reactFlowTime_ > 0.0 && reactFlowTime_ < ignDurationTime_) + { + const vectorField& centre = this->mesh().cellCentres(); + const scalarField x = centre.component(vector::X); + const scalarField y = centre.component(vector::Y); + const scalarField z = centre.component(vector::Z); + + forAll(cCells_,celli) + { + scalar R = Foam::sqrt(magSqr(x[celli]-x0_) + magSqr(y[celli]-y0_) + + magSqr(z[celli]-z0_)); + if(R <= R0_) {cCells_[celli] = 1.0;} + } + + Info<< "Flame initialisation done"<< endl; + } + + +} + +template +Foam::tmp +Foam::combustionModels::baseFGM::R(volScalarField& Y) const +{ + return laminar::R(Y); +} + +template +Foam::tmp +Foam::combustionModels::baseFGM::Qdot() const +{ + return volScalarField::New + ( + this->thermo().phasePropertyName("Qdot"), + laminar::Qdot() + ); +} + + +template +bool Foam::combustionModels::baseFGM::buffer() +{ + if(this->mesh().time().value() - this->mesh().time().startTime().value() < bufferTime_) buffer_ = true; + else buffer_ = false; + + return buffer_; +} + + +// ************************************************************************* // diff --git a/src/dfCombustionModels/FGM/baseFGM/baseFGM.H b/src/dfCombustionModels/FGM/baseFGM/baseFGM.H new file mode 100644 index 00000000..19eb378c --- /dev/null +++ b/src/dfCombustionModels/FGM/baseFGM/baseFGM.H @@ -0,0 +1,295 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2019 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::combustionModels::baseFGM + +Description + Partially stirred reactor turbulent combustion model. + + This model calculates a finite rate, based on both turbulence and chemistry + time scales. Depending on mesh resolution, the Cmix parameter can be used + to scale the turbulence mixing time scale. + +SourceFiles + baseFGM.C + +\*---------------------------------------------------------------------------*/ + +#ifndef baseFGM_H +#define baseFGM_H + +#include "../laminar/laminar.H" +#include "uniformDimensionedFields.H" +#include "CanteraMixture.H" +//#include "basicThermo.H" +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace combustionModels +{ + +/*---------------------------------------------------------------------------*\ + Class baseFGM Declaration +\*---------------------------------------------------------------------------*/ + +template +class baseFGM +: + public laminar +{ + // Protected Data + protected: + + //- buffer switch + Switch buffer_; + + //- scaled PV switch + Switch scaledPV_; + + //- incompressible pressure reference + scalar incompPref_; + + //- ignition switch + Switch ignition_; + + //- combustion switch + Switch combustion_; + + //- enthalpy solve witch + Switch solveEnthalpy_; + + //- T-lookup switch + Switch flameletT_; + + //- psi=\rho/p + volScalarField& psi_; + + //- air mass/mole(28.96kg/mol) + volScalarField Wt_; + + //- heat capacity(m^2*s^{-2}*K^{-1}) + volScalarField Cp_; + + //- mixture fraction Z + volScalarField Z_; + + //- mixture fraction variance Zvar + volScalarField Zvar_; + + //- enthalpy + volScalarField He_; + + //- formation enthalpy(m^2*s^{-2}) + volScalarField Hf_; + + //- progress variable c + volScalarField c_; + + //- progress variable c variance + volScalarField cvar_; + + //- co-variance Z, c + volScalarField Zcvar_; + + //- chemistry source reaction rate(kg*m^{-3}*s^{-1}) + volScalarField omega_c_; + + //- progress variable dissipation rate (s^{-1}) + volScalarField chi_c_; + + //- mixture fraction dissipation rate (s^{-1}) + volScalarField chi_Z_; + + //-covariance Z-c dissipation rate (s^{-1}) + volScalarField chi_Zc_; + + //- default species needed to lookup + volScalarField YCO2_; + volScalarField YCO_; + volScalarField YH2O_; + + //- self-defined species name + hashedWordList speciesNames_; + + //- self-defined speices mass fraction + PtrList Y_; + + //- copy of omega_c_ + volScalarField cOmega_c_; + volScalarField ZOmega_c_; + + + //- internal fields + scalarField& WtCells_; + scalarField& CpCells_; + scalarField& ZCells_; + scalarField& ZvarCells_; + scalarField& HCells_; + scalarField& HfCells_; + scalarField& cCells_; + scalarField& cvarCells_; + scalarField& ZcvarCells_; + scalarField& omega_cCells_; + scalarField& cOmega_cCells_; + scalarField& ZOmega_cCells_; + scalarField& chi_cCells_; + scalarField& chi_ZCells_; + scalarField& chi_ZcCells_; + + scalarField& YCO2Cells_; + scalarField& YCOCells_; + scalarField& YH2OCells_; + + //- max and min Z + scalar ZMax_; + scalar ZMin_; + + //- max and min Z variance + scalar ZvarMax_; + scalar ZvarMin_; + + //- max and min c + scalar cMax_; + scalar cMin_; + scalar Ycmaxall_; + + //- max and min c variance + scalar cvarMax_; + scalar cvarMin_; + + //- max and min Zc variance + scalar ZcvarMax_; + scalar ZcvarMin_; + + //- fluid density + volScalarField& rho_; + + //- fluid pressure + const volScalarField& p_; + + //- fluid temperature + volScalarField& T_; + + //- fluid velocity + const volVectorField& U_; + + // pressure derivate + const volScalarField& dpdt_; + + //- fluid flux + const surfaceScalarField& phi_; + + scalarField& TCells_; + + //- ignition begin time + scalar ignBeginTime_; + + //- ignition duration time + scalar ignDurationTime_; + + //- reactFlowTime_=currentTime - ignBeginTime_ + scalar reactFlowTime_; + + //- igntion origin point + scalar x0_; + scalar y0_; + scalar z0_; + + //- ignition radius + scalar R0_; + + //- turbulent schmidt number + scalar Sct_; + + //- laminar schmidt number + scalar Sc_; + + scalar bufferTime_; + + //- relaxation switch + Switch relaxation_; + + //- Dpdt switch + Switch DpDt_; + + multivariateSurfaceInterpolationScheme::fieldTable fields_; + + public: + //- Runtime type information + TypeName("baseFGM"); + + // Constructors + + //- Construct from components + baseFGM + ( + const word& modelType, + ReactionThermo& thermo, + const compressibleTurbulenceModel& turb, + const word& combustionProperties + ); + + + //- Destructor + virtual ~baseFGM(); + + + // Member Functions + + //- transport equation + virtual void transport(); + + //- initialize flame kernel + virtual void initialiseFalmeKernel(); + + //- Fuel consumption rate matrix. + virtual tmp R(volScalarField& Y) const; + + //- Heat release rate [kg/m/s^3] + virtual tmp Qdot() const; + + //- reaturn buffer_ + virtual bool buffer(); + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace combustionModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "baseFGM.C" +#endif + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dfCombustionModels/FGM/baseFGM/baseFGMs.C b/src/dfCombustionModels/FGM/baseFGM/baseFGMs.C new file mode 100644 index 00000000..c5d60d16 --- /dev/null +++ b/src/dfCombustionModels/FGM/baseFGM/baseFGMs.C @@ -0,0 +1,41 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2018 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 "baseFGM.H" +#include "makeCombustionTypes.H" + +#include "basicThermo.H" + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +makeCombustionTypes(baseFGM, basicThermo); + +} + +// ************************************************************************* // diff --git a/src/dfCombustionModels/FGM/flameletTableSolver/readThermChemTables.H b/src/dfCombustionModels/FGM/flameletTableSolver/readThermChemTables.H new file mode 100644 index 00000000..dbea81fd --- /dev/null +++ b/src/dfCombustionModels/FGM/flameletTableSolver/readThermChemTables.H @@ -0,0 +1,199 @@ + Info<< "Reading " << H_fuel << "\n" << endl; + Info<< "Reading " << H_ox << "\n" << endl; + + z_Tb5={new double[NZL]{}}; + sl_Tb5={new double[NZL]{}}; + th_Tb5={new double[NZL]{}}; + tau_Tb5={new double[NZL]{}}; + kctau_Tb5={new double[NZL]{}}; + + Info<< "Reading laminar flame properties\n" << endl; + + for(int jj = 0; jj < NZL; jj++) + { + f1=fscanf( table, "%lf %lf %lf %lf %lf\n", + &z_Tb5[jj],&sl_Tb5[jj],&th_Tb5[jj], + &tau_Tb5[jj],&kctau_Tb5[jj] ); + } + + fNZ=fscanf(table, "%d",&NZ); + z_Tb3={ new double[NZ]{} }; + for(int ii = 0; ii < NZ; ii++) + { + f2=fscanf( table, "%lf\n", &z_Tb3[ii] ); //load z to z_Tb3 + } + + fNC=fscanf(table, "%d",&NC); + c_Tb3={ new double[NC]{} }; + for(int ii = 0; ii < NC; ii++) + { + f3=fscanf( table, "%lf\n", &c_Tb3[ii] ); ///load c to c_Tb3 + } + + fNGZ=fscanf(table, "%d",&NGZ); + gz_Tb3={ new double[NGZ]{} }; + for(int ii = 0; ii < NGZ; ii++) + { + f4=fscanf( table, "%lf\n", &gz_Tb3[ii] ); //load z-varance to gz_Tb3 + } + + fNGC=fscanf(table, "%d",&NGC); + gc_Tb3={ new double[NGC]{} }; + for(int ii = 0; ii < NGC; ii++) + { + f5=fscanf( table, "%lf\n", &gc_Tb3[ii] ); //load c-variance to gc_Tb3 + } + + fNZC=fscanf(table, "%d",&NZC); + gzc_Tb3={ new double[NZC]{} }; + for(int ii = 0; ii < NZC; ii++) + { + f6=fscanf( table, "%lf\n", &gzc_Tb3[ii] ); //load zc-vovariance to gzc_Tb3 + } + + fNS_NY=fscanf(table, "%d %d",&NS,&NY); + + omgc_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }; + cOc_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }; + ZOc_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }; + cp_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }; + hiyi_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }; + mwt_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }; + Tf_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }; + nu_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }; + Ycmax_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }; + + Yi01_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }, + Yi02_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }, + Yi03_Tb3={ new double[NZ*NC*NGZ*NGC*NZC]{} }; + + Info << "NS is read as: " << NS << endl; + + if(NS == 8) + { + scaledPV_ = true; + Info<< "=============== Using scaled PV ===============" + << "\n" << endl; + } + else if(NS == 9) + { + scaledPV_ = false; + Info<< "=============== Using unscaled PV ===============" + << "\n" << endl; + } + else + { + WarningInFunction << "Number of columns wrong in flare.tbl !!!" + << "\n" << endl; + } + + const char* fmt8{"%lf %lf %lf %lf %lf %lf %lf %lf"}; + const char* fmt9{"%lf %lf %lf %lf %lf %lf %lf %lf %lf"}; + + //Note：if nYis is changed，fmt_nYis is also needed to change！ + const char* fmt_nYis{"%lf %lf %lf"}; + + + Info<< "Reading turbulent flame properties\n" << endl; + + int count = 0; + for(int ii = 0; ii < NZ; ii++) + { + for(int jj = 0; jj < NC; jj++) + { + for(int kk = 0; kk < NGZ; kk++) + { + for(int ll = 0; ll < NGC; ll++) + { + for(int mm = 0; mm < NZC; mm++) + { + if(scaledPV_) + { + f7=fscanf + ( + table,fmt8,&omgc_Tb3[count],&cOc_Tb3[count],&ZOc_Tb3[count],&cp_Tb3[count], + &mwt_Tb3[count],&hiyi_Tb3[count],&Tf_Tb3[count],&nu_Tb3[count] + ); + } + else + { + f8=fscanf + ( + table,fmt9,&omgc_Tb3[count],&cOc_Tb3[count],&ZOc_Tb3[count],&cp_Tb3[count], + &mwt_Tb3[count],&hiyi_Tb3[count],&Tf_Tb3[count],&nu_Tb3[count],&Ycmax_Tb3[count] + ); + } + + if(NY > 0) + { + + f9=fscanf( + table,fmt_nYis,&Yi01_Tb3[count],&Yi02_Tb3[count],&Yi03_Tb3[count] + ); + + //- loop speciesNames + /* forAll(Yi_Tb3,speciesI) + { + f9=fscanf( + table,fmt_nYis,&Yi_Tb3[speciesI][count] + ); //再读入'H2O','CO','CO2'的数据 + } + */ + } + + count++; + } + } + } + } + if(ii%10 == 0) + { + prog = 100.0*ii/(NZ-1); + Info<< " Progress -- " << prog << "%" << endl; + } + } + + //- findthe maxmum PV value from Ycmax_Tb3 + const scalar Ycmaxall{*std::max_element(Ycmax_Tb3,Ycmax_Tb3+count+1)}; + cMaxAll_ = Ycmaxall; + + if(!scaledPV_) + { + Info<< "\nunscaled PV -- Ycmaxall = "< (1.0-small)) gvar = 0.0; + else gvar = var/(mean*(1.-mean)) ; + } + else + { + if(mean < small || mean > (1.0-small)) gvar = 0.0; + else gvar = var/(mean*(Ycmax-mean)) ; + } + + gvar = min(1.0,gvar); + gvar = max(smaller,gvar); + + return gvar; + +} + +double Foam::tableSolver::cal_gcor +( + double Z, + double c, + double Zvar, + double cvar, + double Zcvar +) +{ + double gcor; + + if(abs(Zcvar) < 1e-12) + { + gcor=0.0; + } + else if(abs(cvar) < 1e-12 || abs(Zvar) < 1e-12) + { + gcor=0.0; + } + else + { + gcor = (Zcvar - Z * c) /(Foam::sqrt(Zvar * cvar)); + } + + gcor = fmin(1.0,gcor); + gcor = fmax(-1.0,gcor); + + return gcor; +} + +int Foam::tableSolver::locate_lower +( + int n, + double array[], + double x +) +{ + int i,loc; + + loc = 0; + if(x <= array[loc]) return loc; + + loc = n-2; + if(x >= array[loc+1]) return loc; + + for(i=0; i= array[i] && x < array[i+1]) + { + loc = i; + return loc; + } + } + + return loc; +} + +double Foam::tableSolver::intfac +( + double xx, + double low, + double high +) +{ + double fac; + + if(xx <= low) fac = 0.0; + else if(xx >= high) fac = 1.0; + else fac = (xx-low)/(high-low); + + return fac; +} + +double Foam::tableSolver::interp1d +( + int n1, + int loc_z, + double zfac, + double table_1d[] +) +{ + int i1,j1; + double factor,result; + + result =0.0; + for(i1=0; i1<2; i1++) + { + factor = (1.0-zfac+i1*(2.0*zfac-1.0)); + + if(i1 == 1) j1 = loc_z+1; + else j1 = loc_z; + + result = result + factor*table_1d[j1]; + } + + return result; +} + +double Foam::tableSolver::interp5d +( + int nz, int nc, int ngz, int ngc, int ngcz, + int loc_z, int loc_c, int loc_gz, int loc_gc, + int loc_gcz,double zfac, double cfac, double gzfac, + double gcfac,double gczfac, double table_5d[] +) +{ + int i1,i2,i3,i4,i5,j1,j2,j3,j4,j5,loc; + double factor, result; + + result =0.0; + for(i1=0; i1<2; i1++) + { + if(i1 == 1) j1 = loc_z+1; + else j1 = loc_z; + + for(i2=0; i2<2; i2++) + { + if(i2 == 1) j2 = loc_c+1; + else j2 = loc_c; + + for(i3=0; i3<2; i3++) + { + if(i3 == 1) j3 = loc_gz+1; + else j3 = loc_gz; + + for(i4=0; i4<2; i4++) + { + if(i4 == 1) j4 = loc_gc+1; + else j4 = loc_gc; + + for(i5=0; i5<2; i5++) + { + factor = (1.0-zfac+i1*(2.0*zfac-1.0)) + *(1.0-cfac+i2*(2.0*cfac-1.0)) + *(1.0-gzfac+i3*(2.0*gzfac-1.0)) + *(1.0-gcfac+i4*(2.0*gcfac-1.0)) + *(1.0-gczfac+i5*(2.0*gczfac-1.0)); + + if(i5 == 1) j5 = loc_gcz+1; + else j5 = loc_gcz; + + loc = j1*nc*ngz*ngc*ngcz + +j2*ngz*ngc*ngcz + +j3*ngc*ngcz + +j4*ngcz + +j5; + result = result + factor*table_5d[loc] ; + } + } + } + } + } + + return result; +} + +double Foam::tableSolver::interp2d +( + int n1, int n2, int loc_z, int loc_gz, + double zfac,double gzfac, double table2d[] +) +{ + int i1,i2,j1,j2,loc; + double factor,result; + + result = 0.0; + for(i1=0; i1<2; i1++) + { + if(i1 == 1) j1 = loc_z+1; + else j1 = loc_z; + + for(i2=0; i2<2; i2++) + { + factor = (1.0-zfac+i1*(2.0*zfac-1.0)) + *(1.0-gzfac+i2*(2.0*gzfac-1.0)); + + if(i2 == 1) j2 = loc_gz+1; + else j2 = loc_gz; + + loc = j1*n2+j2; + result = result + factor*table2d[loc]; + } + } + + return result; + +} + + +double Foam::tableSolver::sdrFLRmodel +( + double cvar, double uSgs_pr, double filterSize, + double sl, double dl, double tau, double kc_s,double beta +) +{ + double KaSgs,c3,c4,kc,DaSgs,theta_5; + //- The model parameter K_c^* is also obtained from the laminar flame calculation [Kolla et al., 2009] + //- and this parameter varies with Z for partially premixed flames + kc=kc_s*tau; + theta_5=0.75; + if(uSgs_pr < 1.0E-6) uSgs_pr = smaller; + + KaSgs = Foam::pow((uSgs_pr/sl),1.5) * Foam::pow((dl/filterSize),0.5); + DaSgs = sl*filterSize/uSgs_pr/dl; + + c3=1.5/(1.0/Foam::sqrt(KaSgs)+1.0); + c4=1.1/Foam::pow((1.0+KaSgs),0.4); + + return (1-Foam::exp(-theta_5*filterSize/dl))*cvar/beta + *( 2.0*kc*sl/dl + (c3 - tau*c4*DaSgs) + *(2.0/3.0*uSgs_pr/filterSize) ); +} + + +double Foam::tableSolver::sdrLRXmodel +( + double Csdr, double nut, double delta, double var +) +{ + return Csdr*nut/sqr(delta)*var; +} + +double Foam::tableSolver::lookup1d +( + int n1, double list_1[], double x1, double table_1d[] +) +{ + int loc_1 = locate_lower(n1,list_1,x1); + double fac_1 = intfac(x1,list_1[loc_1],list_1[loc_1+1]); + return interp1d(n1,loc_1,fac_1,table_1d); +} + +double Foam::tableSolver::lookup2d +( + int n1, double list_1[], double x1, + int n2, double list_2[], double x2, + double table_2d[] +) +{ + int loc_1 = locate_lower(n1,list_1,x1); + double fac_1 = intfac(x1,list_1[loc_1],list_1[loc_1+1]); + int loc_2 = locate_lower(n2,list_2,x2); + double fac_2 = intfac(x2,list_2[loc_2],list_2[loc_2+1]); + + return interp2d(n1,n2,loc_1,loc_2,fac_1,fac_2,table_2d); +} + +double Foam::tableSolver::lookup5d +( + int n1, double list_1[], double x1, + int n2, double list_2[], double x2, + int n3, double list_3[], double x3, + int n4, double list_4[], double x4, + int n5, double list_5[], double x5, + double table_5d[] +) +{ + int loc_1 = locate_lower(n1,list_1,x1); + double fac_1 = intfac(x1,list_1[loc_1],list_1[loc_1+1]); + int loc_2 = locate_lower(n2,list_2,x2); + double fac_2 = intfac(x2,list_2[loc_2],list_2[loc_2+1]); + int loc_3 = locate_lower(n3,list_3,x3); + double fac_3 = intfac(x3,list_3[loc_3],list_3[loc_3+1]); + int loc_4 = locate_lower(n4,list_4,x4); + double fac_4 = intfac(x4,list_4[loc_4],list_4[loc_4+1]); + int loc_5 = locate_lower(n5,list_5,x5); + double fac_5 = intfac(x5,list_5[loc_5],list_5[loc_5+1]); + + return interp5d(n1,n2,n3,n4,n5,loc_1,loc_2,loc_3,loc_4,loc_5, + fac_1,fac_2,fac_3,fac_4,fac_5,table_5d); +} + + +double Foam::tableSolver::RANSsdrFLRmodel +( + double cvar, double epsilon, double k, double nu, + double sl, double dl, double tau, double kc_s,double rho +) +{ + double beta=6.7, C3, C4, Kc, Ka; + + Ka=(dl/sl)/Foam::sqrt(nu/epsilon); + Kc=kc_s*tau; + C3=1.5*Foam::sqrt(Ka)/(1+Foam::sqrt(Ka)); + C4=1.1/Foam::pow((1+Ka),0.4); + + return rho/beta*cvar*( (2.0*Kc-tau*C4)*sl/dl + C3*epsilon/k ); +} + + + +} // End Foam namespace \ No newline at end of file diff --git a/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.H b/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.H new file mode 100644 index 00000000..4133f9f8 --- /dev/null +++ b/src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.H @@ -0,0 +1,192 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2021 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + + 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::tableSolver + +Description + class for the interface between table look-up and combustion model. + +SourceFiles + tableSolver.C +\*---------------------------------------------------------------------------*/ + +#ifndef tableSolver_H +#define tableSolver_H + + +#include "hashedWordList.H" +#include "labelList.H" +#include "scalarList.H" +#include "scalarField.H" +#include "scalar.H" +#include "dimensionedScalar.H" +#include "Switch.H" +#include "PtrList.H" + +namespace Foam +{ + +class tableSolver +{ + //- protected data + public: + double small; + double smaller; + double smallest; + double T0; + + + //- create table + FILE *table; + + //- self-defined species name + hashedWordList speciesNames_; + + //- enthalpy solve witch + Switch& scaledPV_; + + //- T-lookup switch + Switch flameletT_; + + //- unscaled max C + scalar& cMaxAll_; + //PtrList Yi_Tb3; + + int NZL,NZ,NC,NGZ,NGC,NZC,NY,NS; + + double Hox,Hfu,prog; + + int fHox_fu, fNZL, fNZ, fNC, fNGZ, fNGC, fNZC, fNS_NY; + + int f1, f2, f3, f4, f5, f6, f7, f8, f9, f10; + + //- fuel and oxidizer enthalpy + const dimensionedScalar H_fuel; + const dimensionedScalar H_ox; + + + double *z_Tb5, //- interpolation z_space + *sl_Tb5, //- laminar flame speed + *th_Tb5, //- laminar flame thickness + *tau_Tb5, //- heat release parameter + *kctau_Tb5; //- for calculating chi_c + + double *z_Tb3, //- z space + *c_Tb3, //- c space + *gz_Tb3, //- z-variance space + *gc_Tb3, //- c-variance space + *gzc_Tb3; //- z-c co-variance space + + + double *omgc_Tb3, //- chemistry reaction source rate + *cOc_Tb3, //- c*omega_c + *ZOc_Tb3, //- z*omega_c + *cp_Tb3, //- heat capacity + *hiyi_Tb3, //- formation enthalpy + *mwt_Tb3, //- mixture mole mass fraction + *Tf_Tb3, //- flame temperature + *nu_Tb3, //- viscosity + *Ycmax_Tb3; //- maximum progress variable + + double *Yi01_Tb3, //- mass fraction of H2O + *Yi02_Tb3, //- mass fraction of CO + *Yi03_Tb3; //- mass fraction of CO2 + + + double *d2Yeq_Tb2; + + + public: + + //- Constructor + tableSolver(wordList speciesNames, Switch& scaledPV_, Switch flameletT_, scalar& cMaxAll_); + + //- Member function + + //- solve variance + double cal_gvar(double mean, double var, double Ycmax = -1.0); + + //- Compute correlation coefficient between Z and c + double cal_gcor(double Z, double c, double Zvar, double cvar, double Zcvar); + + //- Locate lower index for 1D linear interpolation + int locate_lower(int n, double array[], double x); + + //- Compute 1D linear interpolation factor + double intfac(double xx, double low, double high); + + //- 1D linear interpolation + double interp1d(int n1, int loc_z, double zfac, double table_1d[]); + + //- 5D linear interpolation + double interp5d(int nz, int nc, int ngz, int ngc, int ngcz + ,int loc_z, int loc_c, int loc_gz, int loc_gc, int loc_gcz + ,double zfac, double cfac, double gzfac, double gcfac + ,double gczfac, double table_5d[]); + + //- 2D linear interpolation + double interp2d(int n1, int n2, int loc_z, int loc_gz, double zfac + ,double gzfac, double table2d[]); + + + //- SDR model for progress variable + double sdrFLRmodel(double cvar, double uSgs_pr, double filterSize + ,double sl, double dl, double tau, double kc_s + ,double beta); + + + //- SDR model for passive scalar + double sdrLRXmodel(double Csdr, double nut, double delta, double var); + + + //- 1D table look-up + double lookup1d(int n1, double list_1[], double x1, double table_1d[]); + + //- 2D table look-up + double lookup2d(int n1, double list_1[], double x1,int n2, + double list_2[], double x2, double table_2d[]); + + + //- 5D table look-up + double lookup5d(int n1, double list_1[], double x1, + int n2, double list_2[], double x2, + int n3, double list_3[], double x3, + int n4, double list_4[], double x4, + int n5, double list_5[], double x5, + double table_5d[]); + + + //- RANS SDR model for progress variable + double RANSsdrFLRmodel(double cvar, double epsilon, double k, double nu, + double sl, double dl, double tau, double kc_s,double rho); + + //- Destructor + virtual ~tableSolver(); + +}; + + + +} // End Foam namespace + +#endif /* TABLESOLVER_H_ */ diff --git a/src/dfCombustionModels/FGM/flareFGM/flareFGM.C b/src/dfCombustionModels/FGM/flareFGM/flareFGM.C new file mode 100644 index 00000000..62c1b7f0 --- /dev/null +++ b/src/dfCombustionModels/FGM/flareFGM/flareFGM.C @@ -0,0 +1,457 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2018 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 "flareFGM.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::combustionModels::flareFGM::flareFGM +( + const word& modelType, + ReactionThermo& thermo, + const compressibleTurbulenceModel& turb, + const word& combustionProperties +) +: + baseFGM(modelType, thermo, turb, combustionProperties), + tableSolver( + baseFGM::speciesNames_, + baseFGM::scaledPV_, + baseFGM::flameletT_, + baseFGM::Ycmaxall_ + ) +{ + //- retrieval data from table + retrieval(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::combustionModels::flareFGM::~flareFGM() +{} + +template +void Foam::combustionModels::flareFGM::correct() +{ + //- initialize flame kernel + baseFGM::initialiseFalmeKernel(); + + //- solve transport equation + baseFGM::transport(); + + //update enthalpy using lookup data + if(!(this->solveEnthalpy_)) + { + this->He_ = this->Z_*(H_fuel-H_ox) + H_ox; + } + + //- retrieval data from table + retrieval(); + +} + +template +void Foam::combustionModels::flareFGM::retrieval() +{ + + tmp tk(this->turbulence().k()); + volScalarField& k = const_cast(tk()); + scalarField& kCells =k.primitiveFieldRef(); + + tmp tepsilon(this->turbulence().epsilon()); + volScalarField& epsilon = const_cast(tepsilon()); + const scalarField& epsilonCells =epsilon.primitiveFieldRef(); + + tmp tmu = this->turbulence().mu(); + volScalarField& mu = const_cast(tmu()); + scalarField& muCells = mu.primitiveFieldRef(); + + + //- calculate reacting flow solution + const scalar Zl{z_Tb5[0]}; + const scalar Zr{z_Tb5[NZL-1]}; + + dimensionedScalar TMin("TMin",dimensionSet(0,0,0,1,0,0,0),200.0); + dimensionedScalar TMax("TMax",dimensionSet(0,0,0,1,0,0,0),3000.0); + + forAll(this->rho_, celli) + { + + this->chi_ZCells_[celli] = 1.0*epsilonCells[celli]/kCells[celli] *this->ZvarCells_[celli]; + + this->chi_ZcCells_[celli] = 1.0*epsilonCells[celli]/kCells[celli] *this->ZcvarCells_[celli]; + + if(this->ZCells_[celli] >= Zl && this->ZCells_[celli] <= Zr + && this->combustion_ && this->cCells_[celli] > this->small) + { + + double kc_s = this->lookup1d(NZL,z_Tb5,this->ZCells_[celli],kctau_Tb5); + double tau = this->lookup1d(NZL,z_Tb5,this->ZCells_[celli],tau_Tb5); + double sl = this->lookup1d(NZL,z_Tb5,this->ZCells_[celli],sl_Tb5); + double dl = this->lookup1d(NZL,z_Tb5,this->ZCells_[celli],th_Tb5); + + this->chi_cCells_[celli] = + this->RANSsdrFLRmodel(this->cvarCells_[celli],epsilonCells[celli], + kCells[celli],muCells[celli]/this->rho_[celli], + sl,dl,tau,kc_s,this->rho_[celli]); + + } + else + { + this->chi_cCells_[celli] = 1.0*epsilonCells[celli]/kCells[celli]*this->cvarCells_[celli]; + } + + // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = + + double gz{cal_gvar(this->ZCells_[celli],this->ZvarCells_[celli])}; + double gcz{cal_gcor(this->ZCells_[celli],this->cCells_[celli],this->ZvarCells_[celli],this->cvarCells_[celli],this->ZcvarCells_[celli])}, + Ycmax{-1.0},cNorm{},gc{}; + + if(scaledPV_) + { + cNorm = this->cCells_[celli]; + } + else + { + Ycmax = this->lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,0.0, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,0.0, + NZC,gzc_Tb3,0.0, + Ycmax_Tb3); + Ycmax = max(this->smaller,Ycmax); + cNorm = this->cCells_[celli]/Ycmax; + } + + gc = cal_gvar(this->cCells_[celli],this->cvarCells_[celli],Ycmax); + + this->WtCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + mwt_Tb3); + + muCells[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + nu_Tb3)*this->rho_[celli]; + + // -------------------- Yis begin ------------------------------ + if(NY > 0) + { + this->YH2OCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + Yi01_Tb3); + this->YCOCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + Yi02_Tb3); + this->YCO2Cells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + Yi03_Tb3); + } + + // -------------------- Yis end ------------------------------ + + if(this->ZCells_[celli] >= Zl && this->ZCells_[celli] <= Zr + && this->combustion_ && this->cCells_[celli] > this->small) + { + this->omega_cCells_[celli] = + lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + omgc_Tb3) + + ( + scaledPV_ + ? this->chi_ZCells_[celli]*this->cCells_[celli] + *lookup2d(NZ,z_Tb3,this->ZCells_[celli], + NGZ,gz_Tb3,gz,d2Yeq_Tb2) + : 0.0 + ); + + this->cOmega_cCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + cOc_Tb3); + + this->ZOmega_cCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + ZOc_Tb3); + + } + else + { + this->omega_cCells_[celli] = 0.0; + this->cOmega_cCells_[celli] = 0.0; + this->ZOmega_cCells_[celli] = 0.0; + } + + if(flameletT_) + { + this->TCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + Tf_Tb3); + } + else + { + this->CpCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + cp_Tb3); + + this->HfCells_[celli] = lookup5d(NZ,z_Tb3,this->ZCells_[celli], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + hiyi_Tb3); + + this->TCells_[celli] = (this->HCells_[celli]-this->HfCells_[celli])/this->CpCells_[celli] + + this->T0; + } + + this->omega_cCells_[celli] = this->omega_cCells_[celli]*this->rho_[celli]; + this->cOmega_cCells_[celli] = this->cOmega_cCells_[celli]*this->rho_[celli]; + this->ZOmega_cCells_[celli] = this->ZOmega_cCells_[celli]*this->rho_[celli]; + + } + + + //-----------update boundary--------------------------------- + + forAll(this->rho_.boundaryFieldRef(), patchi) + { + fvPatchScalarField& pZ = this->Z_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pZvar = this->Zvar_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pH = this->He_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pWt = this->Wt_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pCp = this->Cp_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pHf = this->Hf_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pT = this->T_.boundaryFieldRef()[patchi]; + fvPatchScalarField& prho_ = this->rho_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pomega_c = this->omega_c_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pcOmega_c = this->cOmega_c_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pZOmega_c = this->ZOmega_c_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pc = this->c_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pcvar = this->cvar_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pZcvar = this->Zcvar_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pchi_Z = this->chi_Z_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pchi_c = this->chi_c_.boundaryFieldRef()[patchi]; + fvPatchScalarField& pchi_Zc = this->chi_Zc_.boundaryFieldRef()[patchi]; + + tmp tmuw = this->turbulence().mu(patchi); + scalarField& pmu = const_cast(tmuw()); + + fvPatchScalarField& pk = k.boundaryFieldRef()[patchi]; + fvPatchScalarField& pepsilon = epsilon.boundaryFieldRef()[patchi]; + + + forAll(prho_, facei) + { + + pchi_Z[facei] = 1.0*pepsilon[facei]/pk[facei] *pZvar[facei]; + + pchi_Zc[facei]= 1.0*pepsilon[facei]/pk[facei] *pZcvar[facei]; + + if(pZ[facei] >= Zl && pZ[facei] <= Zr + && this->combustion_ && pc[facei] > this->small) + { + + double kc_s = lookup1d(NZL,z_Tb5,pZ[facei],kctau_Tb5); + double tau = lookup1d(NZL,z_Tb5,pZ[facei],tau_Tb5); + double sl = lookup1d(NZL,z_Tb5,pZ[facei],sl_Tb5); + double dl = lookup1d(NZL,z_Tb5,pZ[facei],th_Tb5); + + pchi_c[facei] = + RANSsdrFLRmodel(pcvar[facei],pepsilon[facei], + pk[facei],pmu[facei]/prho_[facei], + sl,dl,tau,kc_s,prho_[facei]); + } + else + { + pchi_c[facei] = 1.0*pepsilon[facei]/pk[facei] *pcvar[facei]; + + } + + // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = + + double gz{cal_gvar(pZ[facei],pZvar[facei])}; + double gcz{cal_gcor(pZ[facei],pc[facei],pZvar[facei],pcvar[facei],pZcvar[facei])}, + Ycmax{-1.0},cNorm{},gc{}; + + if(scaledPV_) + { + cNorm = pc[facei]; + } + else + { + Ycmax = lookup5d(NZ,z_Tb3,pZ[facei], + NC,c_Tb3,0.0, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,0.0, + NZC,gzc_Tb3,0.0, + Ycmax_Tb3); + Ycmax = max(this->smaller,Ycmax); + cNorm = pc[facei]/Ycmax; + } + + gc = cal_gvar(pc[facei],pcvar[facei],Ycmax); + + pWt[facei] = lookup5d(NZ,z_Tb3,pZ[facei], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + mwt_Tb3); + + pmu[facei] = lookup5d(NZ,z_Tb3,pZ[facei], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + nu_Tb3)*prho_[facei]; + + if(pZ[facei] >= Zl && pZ[facei] <= Zr + && this->combustion_ && pc[facei] > this->small) + { + pomega_c[facei] = + lookup5d(NZ,z_Tb3,pZ[facei], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + omgc_Tb3) + + ( + scaledPV_ + ? pchi_Z[facei]*pc[facei] + *lookup2d(NZ,z_Tb3,pZ[facei], + NGZ,gz_Tb3,gz,d2Yeq_Tb2) + : 0.0 + ); + + pcOmega_c[facei] = lookup5d(NZ,z_Tb3,pZ[facei], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz,cOc_Tb3); + + pZOmega_c[facei] = lookup5d(NZ,z_Tb3,pZ[facei], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz,ZOc_Tb3); + } + else + { + pomega_c[facei] = 0.0; + pcOmega_c[facei] = 0.0; + pZOmega_c[facei] = 0.0; + } + + if(flameletT_) + { + pT[facei] = lookup5d(NZ,z_Tb3,pZ[facei], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + Tf_Tb3); + } + else + { + pCp[facei] = lookup5d(NZ,z_Tb3,pZ[facei], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + cp_Tb3); + + pHf[facei] = lookup5d(NZ,z_Tb3,pZ[facei], + NC,c_Tb3,cNorm, + NGZ,gz_Tb3,gz, + NGC,gc_Tb3,gc, + NZC,gzc_Tb3,gcz, + hiyi_Tb3); + + pT[facei] = (pH[facei]-pHf[facei])/pCp[facei] + + this->T0; + } + + pomega_c[facei] = pomega_c[facei]*prho_[facei]; + pcOmega_c[facei] = pcOmega_c[facei]*prho_[facei]; + pZOmega_c[facei] = pZOmega_c[facei]*prho_[facei]; + + } + + } + + this->T_.max(TMin); + this->T_.min(TMax); + + if(this->mesh().time().timeIndex() > 0) + { + dimensionedScalar p_operateDim("p_operateDim", dimensionSet(1,-1,-2,0,0,0,0),this->incompPref_); + + if(this->incompPref_ > 0.0) + { + this->rho_ = p_operateDim*this->psi_; + } + else + { + this->rho_ = this->p_*this->psi_; + } + } + + dimensionedScalar R_uniGas("R_uniGas",dimensionSet(1,2,-2,-1,-1,0,0),8.314e3); + this->psi_ = this->Wt_/(R_uniGas*this->T_); +} + +// ************************************************************************* // diff --git a/src/dfCombustionModels/FGM/flareFGM/flareFGM.H b/src/dfCombustionModels/FGM/flareFGM/flareFGM.H new file mode 100644 index 00000000..897a1e48 --- /dev/null +++ b/src/dfCombustionModels/FGM/flareFGM/flareFGM.H @@ -0,0 +1,118 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2019 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::combustionModels::flareFGM + +Description + Partially stirred reactor turbulent combustion model. + + This model calculates a finite rate, based on both turbulence and chemistry + time scales. Depending on mesh resolution, the Cmix parameter can be used + to scale the turbulence mixing time scale. + +SourceFiles + flareFGM.C + +\*---------------------------------------------------------------------------*/ + +#ifndef flareFGM_H +#define flareFGM_H + +#include "baseFGM.H" +#include "tableSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace combustionModels +{ + +/*---------------------------------------------------------------------------*\ + Class flareFGM Declaration +\*---------------------------------------------------------------------------*/ + +template +class flareFGM +: + public baseFGM, + public tableSolver +{ + +public: + + //- Runtime type information + TypeName("flareFGM"); + + + // Constructors + + //- Construct from components + flareFGM + ( + const word& modelType, + ReactionThermo& thermo, + const compressibleTurbulenceModel& turb, + const word& combustionProperties + ); + + //- Disallow default bitwise copy construction + flareFGM(const flareFGM&); + + + //- Destructor + virtual ~flareFGM(); + + + // Member Operators + + //- Correct combustion rate + virtual void correct(); + + //- retrieval data from table + virtual void retrieval(); + + + + //- Disallow default bitwise assignment + void operator=(const flareFGM&) = delete; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace combustionModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "flareFGM.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dfCombustionModels/FGM/flareFGM/flareFGMs.C b/src/dfCombustionModels/FGM/flareFGM/flareFGMs.C new file mode 100644 index 00000000..67b3035f --- /dev/null +++ b/src/dfCombustionModels/FGM/flareFGM/flareFGMs.C @@ -0,0 +1,41 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2018 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 "flareFGM.H" +#include "makeCombustionTypes.H" + +#include "basicThermo.H" + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +makeCombustionTypes(flareFGM, basicThermo); + +} + +// ************************************************************************* // diff --git a/src/dfCombustionModels/Make/files b/src/dfCombustionModels/Make/files index 6cf11d1a..7dc602c6 100644 --- a/src/dfCombustionModels/Make/files +++ b/src/dfCombustionModels/Make/files @@ -7,6 +7,10 @@ CombustionModel/CombustionModel/CombustionModels.C PaSR/PaSRs.C EDC/EDCs.C +FGM/flameletTableSolver/tableSolver.C +FGM/baseFGM/baseFGMs.C +FGM/flareFGM/flareFGMs.C + laminar/laminars.C noCombustion/noCombustions.C diff --git a/src/dfCombustionModels/Make/options b/src/dfCombustionModels/Make/options index e063e7ed..fa21fbee 100644 --- a/src/dfCombustionModels/Make/options +++ b/src/dfCombustionModels/Make/options @@ -10,8 +10,21 @@ EXE_INC = -std=c++14 \ -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(DF_SRC)/lagrangian/intermediate/lnInclude \ + -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ + -I$(DF_SRC)/lagrangian/spray/lnInclude \ + -I$(LIB_SRC)/lagrangian/spray/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(DF_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \ + -I$(DF_SRC)/thermophysicalModels/SLGThermo/lnInclude \ + -I$(LIB_SRC)/regionModels/regionModel/lnInclude \ + -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ -I$(DF_SRC)/dfCanteraMixture/lnInclude \ -I$(DF_SRC)/dfChemistryModel/lnInclude \ + -I$(LIB_SRC)/ODE/lnInclude \ + -I$(LIB_SRC)/functionObjects/field/lnInclude \ -I$(CANTERA_ROOT)/include \ $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include,) \ $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include/torch/csrc/api/include,) \ @@ -20,12 +33,22 @@ EXE_INC = -std=c++14 \ LIB_LIBS = \ -lcompressibleTransportModels \ -lturbulenceModels \ + -llagrangian \ -lfiniteVolume \ -lmeshTools \ -L$(DF_LIBBIN) \ -ldfCompressibleTurbulenceModels \ + -ldfLagrangianIntermediate \ + -ldfLagrangianTurbulence \ + -ldfLagrangianSpray \ + -ldfFluidThermophysicalModels \ + -ldfThermophysicalProperties \ + -ldfSLGThermo \ + -lregionModels \ + -lsurfaceFilmModels \ -ldfCanteraMixture \ -ldfChemistryModel \ + -lODE \ $(CANTERA_ROOT)/lib/libcantera.so \ $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libtorch.so,) \ $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libc10.so,) \