From e54db499035a668624a71dde27c750ca6ee936c9 Mon Sep 17 00:00:00 2001 From: minzhang0929 <78373564+minzhang0929@users.noreply.github.com> Date: Thu, 19 Jan 2023 16:58:30 +0800 Subject: [PATCH 1/8] Update files --- src/dfCombustionModels/Make/files | 4 ++++ 1 file changed, 4 insertions(+) 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 From 13dc006c25b4a538b86057b5da60e6580ffb8a36 Mon Sep 17 00:00:00 2001 From: minzhang0929 <78373564+minzhang0929@users.noreply.github.com> Date: Thu, 19 Jan 2023 16:59:29 +0800 Subject: [PATCH 2/8] Delete options --- src/dfCombustionModels/Make/options | 34 ----------------------------- 1 file changed, 34 deletions(-) delete mode 100644 src/dfCombustionModels/Make/options diff --git a/src/dfCombustionModels/Make/options b/src/dfCombustionModels/Make/options deleted file mode 100644 index e063e7ed..00000000 --- a/src/dfCombustionModels/Make/options +++ /dev/null @@ -1,34 +0,0 @@ --include $(GENERAL_RULES)/mplibType - -EXE_INC = -std=c++14 \ - -Wno-old-style-cast \ - $(if $(LIBTORCH_ROOT),-DUSE_LIBTORCH,) \ - $(if $(PYTHON_INC_DIR),-DUSE_PYTORCH,) \ - -I$(LIB_SRC)/transportModels/compressible/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ - -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ - -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(DF_SRC)/dfCanteraMixture/lnInclude \ - -I$(DF_SRC)/dfChemistryModel/lnInclude \ - -I$(CANTERA_ROOT)/include \ - $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include,) \ - $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include/torch/csrc/api/include,) \ - $(PYTHON_INC_DIR) - -LIB_LIBS = \ - -lcompressibleTransportModels \ - -lturbulenceModels \ - -lfiniteVolume \ - -lmeshTools \ - -L$(DF_LIBBIN) \ - -ldfCompressibleTurbulenceModels \ - -ldfCanteraMixture \ - -ldfChemistryModel \ - $(CANTERA_ROOT)/lib/libcantera.so \ - $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libtorch.so,) \ - $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libc10.so,) \ - $(if $(LIBTORCH_ROOT),-rdynamic,) \ - $(if $(LIBTORCH_ROOT),-lpthread,) - From f70cccb284f5a14f857630a2aaf6c603eec68b0a Mon Sep 17 00:00:00 2001 From: minzhang0929 <78373564+minzhang0929@users.noreply.github.com> Date: Thu, 19 Jan 2023 21:59:48 +1300 Subject: [PATCH 3/8] Add files via upload --- src/dfCombustionModels/Make/options | 57 +++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 src/dfCombustionModels/Make/options diff --git a/src/dfCombustionModels/Make/options b/src/dfCombustionModels/Make/options new file mode 100644 index 00000000..fa21fbee --- /dev/null +++ b/src/dfCombustionModels/Make/options @@ -0,0 +1,57 @@ +-include $(GENERAL_RULES)/mplibType + +EXE_INC = -std=c++14 \ + -Wno-old-style-cast \ + $(if $(LIBTORCH_ROOT),-DUSE_LIBTORCH,) \ + $(if $(PYTHON_INC_DIR),-DUSE_PYTORCH,) \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -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,) \ + $(PYTHON_INC_DIR) + +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,) \ + $(if $(LIBTORCH_ROOT),-rdynamic,) \ + $(if $(LIBTORCH_ROOT),-lpthread,) + From 5f49ac1aed3168587d91f3803b51738a3313c8fc Mon Sep 17 00:00:00 2001 From: minzhang0929 <78373564+minzhang0929@users.noreply.github.com> Date: Thu, 19 Jan 2023 22:01:08 +1300 Subject: [PATCH 4/8] Add files via upload --- src/dfCombustionModels/FGM/baseFGM/baseFGM.C | 610 ++++++++++++++++++ src/dfCombustionModels/FGM/baseFGM/baseFGM.H | 295 +++++++++ src/dfCombustionModels/FGM/baseFGM/baseFGMs.C | 41 ++ .../flameletTableSolver/readThermChemTables.H | 199 ++++++ .../FGM/flameletTableSolver/tableSolver.C | 392 +++++++++++ .../FGM/flameletTableSolver/tableSolver.H | 192 ++++++ .../FGM/flareFGM/flareFGM.C | 457 +++++++++++++ .../FGM/flareFGM/flareFGM.H | 118 ++++ .../FGM/flareFGM/flareFGMs.C | 41 ++ 9 files changed, 2345 insertions(+) create mode 100644 src/dfCombustionModels/FGM/baseFGM/baseFGM.C create mode 100644 src/dfCombustionModels/FGM/baseFGM/baseFGM.H create mode 100644 src/dfCombustionModels/FGM/baseFGM/baseFGMs.C create mode 100644 src/dfCombustionModels/FGM/flameletTableSolver/readThermChemTables.H create mode 100644 src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.C create mode 100644 src/dfCombustionModels/FGM/flameletTableSolver/tableSolver.H create mode 100644 src/dfCombustionModels/FGM/flareFGM/flareFGM.C create mode 100644 src/dfCombustionModels/FGM/flareFGM/flareFGM.H create mode 100644 src/dfCombustionModels/FGM/flareFGM/flareFGMs.C 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); + +} + +// ************************************************************************* // From cd65cfa64021543afe4ebd5d3ac06685b4ae8d48 Mon Sep 17 00:00:00 2001 From: minzhang0929 <78373564+minzhang0929@users.noreply.github.com> Date: Thu, 19 Jan 2023 17:03:25 +0800 Subject: [PATCH 5/8] Delete createFields.H --- .../solvers/dfLowMachFoam/createFields.H | 155 ------------------ 1 file changed, 155 deletions(-) delete mode 100644 applications/solvers/dfLowMachFoam/createFields.H diff --git a/applications/solvers/dfLowMachFoam/createFields.H b/applications/solvers/dfLowMachFoam/createFields.H deleted file mode 100644 index 5c4ccc7e..00000000 --- a/applications/solvers/dfLowMachFoam/createFields.H +++ /dev/null @@ -1,155 +0,0 @@ -#include "createRDeltaT.H" - -Info<< "Reading thermophysical properties\n" << endl; - -fluidThermo* pThermo = new hePsiThermo(mesh, word::null); -fluidThermo& thermo = *pThermo; -thermo.validate(args.executable(), "ha"); - -const volScalarField& psi = thermo.psi(); -volScalarField& p = thermo.p(); -volScalarField& T = thermo.T(); -volScalarField rho -( - IOobject - ( - "rho", - runTime.timeName(), - mesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - thermo.rho() -); - -Info<< "Reading field U\n" << endl; -volVectorField U -( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -#include "compressibleCreatePhi.H" - -pressureControl pressureControl(p, rho, pimple.dict(), false); - -mesh.setFluxRequired(p.name()); - -Info<< "Creating turbulence model\n" << endl; -autoPtr turbulence -( - compressible::turbulenceModel::New - ( - rho, - U, - phi, - thermo - ) -); - -Info<< "Creating reaction model\n" << endl; -autoPtr> combustion -( - CombustionModel::New(thermo, turbulence()) -); -Info<< "end Creating reaction model\n" << endl; -dfChemistryModel* chemistry = combustion->chemistry(); -PtrList& Y = chemistry->Y(); -const word inertSpecie(chemistry->lookup("inertSpecie")); -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) -); - -Info<< "Creating field kinetic energy K\n" << endl; -volScalarField K("K", 0.5*magSqr(U)); - -multivariateSurfaceInterpolationScheme::fieldTable fields; -forAll(Y, i) -{ - fields.add(Y[i]); -} -fields.add(thermo.he()); - -const scalar Sct = chemistry->lookupOrDefault("Sct", 1.); -volScalarField diffAlphaD -( - IOobject - ( - "diffAlphaD", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar(dimEnergy/dimTime/dimVolume, 0) -); -volVectorField hDiffCorrFlux -( - IOobject - ( - "hDiffCorrFlux", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedVector(dimensionSet(1,0,-3,0,0,0,0), Zero) -); -volVectorField sumYDiffError -( - IOobject - ( - "sumYDiffError", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedVector("sumYDiffError", dimDynamicViscosity/dimLength, Zero) -); - -IOdictionary CanteraTorchProperties -( - IOobject - ( - "CanteraTorchProperties", - runTime.constant(), - mesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) -); -const Switch splitting = CanteraTorchProperties.lookupOrDefault("splittingStrategy", false); -#ifdef USE_PYTORCH - const Switch log_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("log", false); - const Switch torch_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("torch", false); -#endif -#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 From bc8f4a1fd02e4ab56a32a10ddcd144caf57d026c Mon Sep 17 00:00:00 2001 From: minzhang0929 <78373564+minzhang0929@users.noreply.github.com> Date: Thu, 19 Jan 2023 17:04:14 +0800 Subject: [PATCH 6/8] Delete dfLowMachFoam.C --- .../solvers/dfLowMachFoam/dfLowMachFoam.C | 239 ------------------ 1 file changed, 239 deletions(-) delete mode 100644 applications/solvers/dfLowMachFoam/dfLowMachFoam.C diff --git a/applications/solvers/dfLowMachFoam/dfLowMachFoam.C b/applications/solvers/dfLowMachFoam/dfLowMachFoam.C deleted file mode 100644 index ec90b78e..00000000 --- a/applications/solvers/dfLowMachFoam/dfLowMachFoam.C +++ /dev/null @@ -1,239 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 . - -Application - rhoPimpleFoam - -Description - Transient solver for turbulent flow of compressible fluids for HVAC and - similar applications, with optional mesh motion and mesh topology changes. - - Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and - pseudo-transient simulations. - -\*---------------------------------------------------------------------------*/ - -#include "dfChemistryModel.H" -#include "CanteraMixture.H" -#include "hePsiThermo.H" - -#ifdef USE_PYTORCH -#include -#include -#include //used to convert -#endif - -#ifdef USE_LIBTORCH -#include -#include "DNNInferencer.H" -#endif - -#include "fvCFD.H" -#include "fluidThermo.H" -#include "turbulentFluidThermoModel.H" -#include "pimpleControl.H" -#include "pressureControl.H" -#include "localEulerDdtScheme.H" -#include "fvcSmooth.H" -#include "PstreamGlobals.H" -#include "basicThermo.H" -#include "CombustionModel.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ -#ifdef USE_PYTORCH - pybind11::scoped_interpreter guard{};//start python interpreter -#endif - #include "postProcess.H" - - // #include "setRootCaseLists.H" - #include "listOptions.H" - #include "setRootCase2.H" - #include "listOutput.H" - - #include "createTime.H" - #include "createMesh.H" - #include "createDyMControls.H" - #include "initContinuityErrs.H" - #include "createFields.H" - #include "createRhoUfIfPresent.H" - - double time_monitor_flow=0; - double time_monitor_chem=0; - double time_monitor_Y=0; - double time_monitor_E=0; - double time_monitor_corrThermo=0; - double time_monitor_corrDiff=0; - label timeIndex = 0; - clock_t start, end; - - turbulence->validate(); - - if (!LTS) - { - #include "compressibleCourantNo.H" - #include "setInitialDeltaT.H" - } - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - Info<< "\nStarting time loop\n" << endl; - - while (runTime.run()) - { - timeIndex ++; - - #include "readDyMControls.H" - - if (LTS) - { - #include "setRDeltaT.H" - } - else - { - #include "compressibleCourantNo.H" - #include "setDeltaT.H" - } - - runTime++; - - Info<< "Time = " << runTime.timeName() << nl << endl; - - // --- Pressure-velocity PIMPLE corrector loop - while (pimple.loop()) - { - if (splitting) - { - #include "YEqn_RR.H" - } - if (pimple.firstPimpleIter() || moveMeshOuterCorrectors) - { - // Store momentum to set rhoUf for introduced faces. - autoPtr rhoU; - if (rhoUf.valid()) - { - rhoU = new volVectorField("rhoU", rho*U); - } - } - - if (pimple.firstPimpleIter() && !pimple.simpleRho()) - { - #include "rhoEqn.H" - } - - start = std::clock(); - #include "UEqn.H" - end = std::clock(); - time_monitor_flow += double(end - start) / double(CLOCKS_PER_SEC); - - #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(); - chemistry->correctThermo(); - end = std::clock(); - time_monitor_corrThermo += double(end - start) / double(CLOCKS_PER_SEC); - - Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl; - - // --- Pressure corrector loop - - start = std::clock(); - while (pimple.correct()) - { - if (pimple.consistent()) - { - #include "pcEqn.H" - } - else - { - #include "pEqn.H" - } - } - end = std::clock(); - time_monitor_flow += double(end - start) / double(CLOCKS_PER_SEC); - - if (pimple.turbCorr()) - { - turbulence->correct(); - } - } - - rho = thermo.rho(); - - runTime.write(); - - Info<< "========Time Spent in diffenet parts========"<< endl; - Info<< "Chemical sources = " << time_monitor_chem << " s" << endl; - Info<< "Species Equations = " << time_monitor_Y << " s" << endl; - Info<< "U & p Equations = " << time_monitor_flow << " s" << endl; - Info<< "Energy Equations = " << time_monitor_E << " s" << endl; - Info<< "thermo & Trans Properties = " << time_monitor_corrThermo << " s" << endl; - Info<< "Diffusion Correction Time = " << time_monitor_corrDiff << " s" << endl; - Info<< "sum Time = " << (time_monitor_chem + time_monitor_Y + time_monitor_flow + time_monitor_E + time_monitor_corrThermo + time_monitor_corrDiff) << " s" << endl; - Info<< "============================================"< Date: Thu, 19 Jan 2023 22:04:45 +1300 Subject: [PATCH 7/8] Add files via upload --- .../solvers/dfLowMachFoam/createFields.H | 172 ++++++++++++++++++ 1 file changed, 172 insertions(+) create mode 100644 applications/solvers/dfLowMachFoam/createFields.H diff --git a/applications/solvers/dfLowMachFoam/createFields.H b/applications/solvers/dfLowMachFoam/createFields.H new file mode 100644 index 00000000..e4c5ba9b --- /dev/null +++ b/applications/solvers/dfLowMachFoam/createFields.H @@ -0,0 +1,172 @@ +#include "createRDeltaT.H" + +Info<< "Reading thermophysical properties\n" << endl; + +fluidThermo* pThermo = new hePsiThermo(mesh, word::null); +fluidThermo& thermo = *pThermo; +thermo.validate(args.executable(), "ha"); + +const volScalarField& psi = thermo.psi(); +volScalarField& p = thermo.p(); +volScalarField& T = thermo.T(); +volScalarField rho +( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + thermo.rho() +); + + +Info<< "Reading field U\n" << endl; +volVectorField U +( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +#include "compressibleCreatePhi.H" + +pressureControl pressureControl(p, rho, pimple.dict(), false); + +mesh.setFluxRequired(p.name()); + +Info<< "Creating turbulence model\n" << endl; +autoPtr turbulence +( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) +); + +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")); +const label inertIndex(chemistry->species()[inertSpecie]); + + +chemistry->correctThermo(); +Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl; + +//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 +( + IOobject + ( + "diffAlphaD", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar(dimEnergy/dimTime/dimVolume, 0) +); +volVectorField hDiffCorrFlux +( + IOobject + ( + "hDiffCorrFlux", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedVector(dimensionSet(1,0,-3,0,0,0,0), Zero) +); +volVectorField sumYDiffError +( + IOobject + ( + "sumYDiffError", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedVector("sumYDiffError", dimDynamicViscosity/dimLength, Zero) +); + +IOdictionary CanteraTorchProperties +( + IOobject + ( + "CanteraTorchProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) +); +const Switch splitting = CanteraTorchProperties.lookupOrDefault("splittingStrategy", false); +#ifdef USE_PYTORCH + const Switch log_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("log", false); + const Switch torch_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("torch", false); +#endif +#ifdef USE_LIBTORCH + const Switch log_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("log", false); + const Switch torch_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("torch", false); +#endif From 41aa57e2ae927c594671b5802da8de35225ef453 Mon Sep 17 00:00:00 2001 From: minzhang0929 <78373564+minzhang0929@users.noreply.github.com> Date: Thu, 19 Jan 2023 22:05:19 +1300 Subject: [PATCH 8/8] Add files via upload --- .../solvers/dfLowMachFoam/dfLowMachFoam.C | 248 ++++++++++++++++++ 1 file changed, 248 insertions(+) create mode 100644 applications/solvers/dfLowMachFoam/dfLowMachFoam.C diff --git a/applications/solvers/dfLowMachFoam/dfLowMachFoam.C b/applications/solvers/dfLowMachFoam/dfLowMachFoam.C new file mode 100644 index 00000000..437866bb --- /dev/null +++ b/applications/solvers/dfLowMachFoam/dfLowMachFoam.C @@ -0,0 +1,248 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Application + rhoPimpleFoam + +Description + Transient solver for turbulent flow of compressible fluids for HVAC and + similar applications, with optional mesh motion and mesh topology changes. + + Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and + pseudo-transient simulations. + +\*---------------------------------------------------------------------------*/ + +#include "dfChemistryModel.H" +#include "CanteraMixture.H" +#include "hePsiThermo.H" + +#ifdef USE_PYTORCH +#include +#include +#include //used to convert +#endif + +#ifdef USE_LIBTORCH +#include +#include "DNNInferencer.H" +#endif + +#include "fvCFD.H" +#include "fluidThermo.H" +#include "turbulentFluidThermoModel.H" +#include "pimpleControl.H" +#include "pressureControl.H" +#include "localEulerDdtScheme.H" +#include "fvcSmooth.H" +#include "PstreamGlobals.H" +#include "basicThermo.H" +#include "CombustionModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +#ifdef USE_PYTORCH + pybind11::scoped_interpreter guard{};//start python interpreter +#endif + #include "postProcess.H" + + // #include "setRootCaseLists.H" + #include "listOptions.H" + #include "setRootCase2.H" + #include "listOutput.H" + + #include "createTime.H" + #include "createMesh.H" + #include "createDyMControls.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "createRhoUfIfPresent.H" + + double time_monitor_flow=0; + double time_monitor_chem=0; + double time_monitor_Y=0; + double time_monitor_E=0; + double time_monitor_corrThermo=0; + double time_monitor_corrDiff=0; + label timeIndex = 0; + clock_t start, end; + + turbulence->validate(); + + if (!LTS) + { + #include "compressibleCourantNo.H" + #include "setInitialDeltaT.H" + } + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + timeIndex ++; + + #include "readDyMControls.H" + + if (LTS) + { + #include "setRDeltaT.H" + } + else + { + #include "compressibleCourantNo.H" + #include "setDeltaT.H" + } + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + if (splitting) + { + #include "YEqn_RR.H" + } + if (pimple.firstPimpleIter() || moveMeshOuterCorrectors) + { + // Store momentum to set rhoUf for introduced faces. + autoPtr rhoU; + if (rhoUf.valid()) + { + rhoU = new volVectorField("rhoU", rho*U); + } + } + + if (pimple.firstPimpleIter() && !pimple.simpleRho()) + { + #include "rhoEqn.H" + } + + start = std::clock(); + #include "UEqn.H" + end = std::clock(); + time_monitor_flow += double(end - start) / double(CLOCKS_PER_SEC); + + 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(); + 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; + + // --- Pressure corrector loop + + start = std::clock(); + while (pimple.correct()) + { + if (pimple.consistent()) + { + #include "pcEqn.H" + } + else + { + #include "pEqn.H" + } + } + end = std::clock(); + time_monitor_flow += double(end - start) / double(CLOCKS_PER_SEC); + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + rho = thermo.rho(); + + 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; + Info<< "U & p Equations = " << time_monitor_flow << " s" << endl; + Info<< "Energy Equations = " << time_monitor_E << " s" << endl; + Info<< "thermo & Trans Properties = " << time_monitor_corrThermo << " s" << endl; + Info<< "Diffusion Correction Time = " << time_monitor_corrDiff << " s" << endl; + Info<< "sum Time = " << (time_monitor_chem + time_monitor_Y + time_monitor_flow + time_monitor_E + time_monitor_corrThermo + time_monitor_corrDiff) << " s" << endl; + Info<< "============================================"<