diff --git a/etc/controlDict b/etc/controlDict index 647656140..9598ac8b2 100644 --- a/etc/controlDict +++ b/etc/controlDict @@ -788,6 +788,7 @@ DebugSwitches solidBodyMotionFunction 0; solidBodyMotionFvMesh 0; solution 0; + specie 0; spectEddyVisc 0; sphereToCell 0; spherical 0; diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C index 3d7047f86..3896761fb 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.C +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.C @@ -191,7 +191,8 @@ void Foam::ODESolver::solve FatalErrorInFunction << "Integration steps greater than maximum " << maxSteps_ << nl << " xStart = " << xStart << ", xEnd = " << xEnd - << ", x = " << x << ", dxDid = " << step.dxDid + << ", x = " << x << ", dxDid = " << step.dxDid << nl + << " y = " << y << exit(FatalError); } diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C index e1993957e..ad2c75577 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C @@ -196,17 +196,17 @@ Foam::scalar Foam::StandardChemistryModel::omega if (c[si] < c[lRef]) { const scalar exp = R.lhs()[slRef].exponent; - pf *= pow(max(0.0, c[lRef]), exp); + pf *= pow(max(c[lRef], 0.0), exp); lRef = si; slRef = s; } else { const scalar exp = R.lhs()[s].exponent; - pf *= pow(max(0.0, c[si]), exp); + pf *= pow(max(c[si], 0.0), exp); } } - cf = max(0.0, c[lRef]); + cf = max(c[lRef], 0.0); { const scalar exp = R.lhs()[slRef].exponent; @@ -238,23 +238,23 @@ Foam::scalar Foam::StandardChemistryModel::omega if (c[si] < c[rRef]) { const scalar exp = R.rhs()[srRef].exponent; - pr *= pow(max(0.0, c[rRef]), exp); + pr *= pow(max(c[rRef], 0.0), exp); rRef = si; srRef = s; } else { const scalar exp = R.rhs()[s].exponent; - pr *= pow(max(0.0, c[si]), exp); + pr *= pow(max(c[si], 0.0), exp); } } - cr = max(0.0, c[rRef]); + cr = max(c[rRef], 0.0); { const scalar exp = R.rhs()[srRef].exponent; if (exp < 1.0) { - if (cr>small) + if (cr > small) { pr *= pow(cr, exp - 1.0); } @@ -284,9 +284,9 @@ void Foam::StandardChemistryModel::derivatives const scalar T = c[nSpecie_]; const scalar p = c[nSpecie_ + 1]; - for (label i = 0; i < nSpecie_; i++) + forAll(c_, i) { - c_[i] = max(0.0, c[i]); + c_[i] = max(c[i], 0.0); } omega(c_, T, p, dcdt); @@ -366,7 +366,7 @@ void Foam::StandardChemistryModel::jacobian { if (c_[si] > small) { - kf *= el*pow(c_[si] + vSmall, el - 1.0); + kf *= el*pow(c_[si], el - 1.0); } else { @@ -412,7 +412,7 @@ void Foam::StandardChemistryModel::jacobian { if (c_[si] > small) { - kr *= er*pow(c_[si] + vSmall, er - 1.0); + kr *= er*pow(c_[si], er - 1.0); } else { @@ -741,6 +741,9 @@ Foam::scalar Foam::StandardChemistryModel::solve deltaTMin = min(this->deltaTChem_[celli], deltaTMin); + this->deltaTChem_[celli] = + min(this->deltaTChem_[celli], this->deltaTChemMax_); + for (label i=0; i::omega if (c[si] < c[lRef]) { const scalar exp = R.lhs()[slRef].exponent; - pf *= pow(max(0.0, c[lRef]), exp); + pf *= pow(max(c[lRef], 0.0), exp); lRef = si; slRef = s; } else { const scalar exp = R.lhs()[s].exponent; - pf *= pow(max(0.0, c[si]), exp); + pf *= pow(max(c[si], 0.0), exp); } } - cf = max(0.0, c[lRef]); + cf = max(c[lRef], 0.0); { const scalar exp = R.lhs()[slRef].exponent; @@ -274,23 +274,23 @@ Foam::scalar Foam::TDACChemistryModel::omega if (c[si] < c[rRef]) { const scalar exp = R.rhs()[srRef].exponent; - pr *= pow(max(0.0, c[rRef]), exp); + pr *= pow(max(c[rRef], 0.0), exp); rRef = si; srRef = s; } else { const scalar exp = R.rhs()[s].exponent; - pr *= pow(max(0.0, c[si]), exp); + pr *= pow(max(c[si], 0.0), exp); } } - cr = max(0.0, c[rRef]); + cr = max(c[rRef], 0.0); { const scalar exp = R.rhs()[srRef].exponent; if (exp < 1) { - if (cr>small) + if (cr > small) { pr *= pow(cr, exp - 1); } @@ -334,14 +334,14 @@ void Foam::TDACChemistryModel::derivatives // efficiencies for (label i=0; ic_[simplifiedToCompleteIndex_[i]] = max(0.0, c[i]); + this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0); } } else { for (label i=0; inSpecie(); i++) { - this->c_[i] = max(0.0, c[i]); + this->c_[i] = max(c[i], 0.0); } } @@ -416,7 +416,7 @@ void Foam::TDACChemistryModel::jacobian this->c_ = completeC_; for (label i=0; ic_[simplifiedToCompleteIndex_[i]] = max(0.0, c[i]); + this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0); } } else @@ -456,7 +456,7 @@ void Foam::TDACChemistryModel::jacobian { if (this->c_[si] > small) { - kf *= el*pow(this->c_[si] + vSmall, el - 1); + kf *= el*pow(this->c_[si], el - 1); } else { @@ -514,7 +514,7 @@ void Foam::TDACChemistryModel::jacobian { if (this->c_[si] > small) { - kr *= er*pow(this->c_[si] + vSmall, er - 1); + kr *= er*pow(this->c_[si], er - 1); } else { @@ -797,6 +797,9 @@ Foam::scalar Foam::TDACChemistryModel::solve this->nSpecie_ = mechRed_->nSpecie(); } deltaTMin = min(this->deltaTChem_[celli], deltaTMin); + + this->deltaTChem_[celli] = + min(this->deltaTChem_[celli], this->deltaTChemMax_); } // Set the RR vector (used in the solver) diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C index b47b71ecd..2a896711b 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -58,6 +58,7 @@ Foam::basicChemistryModel::basicChemistryModel(basicThermo& thermo) mesh_(thermo.p().mesh()), chemistry_(lookup("chemistry")), deltaTChemIni_(readScalar(lookup("initialChemicalTimeStep"))), + deltaTChemMax_(lookupOrDefault("maxChemicalTimeStep", great)), deltaTChem_ ( IOobject diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H index 1d321912a..8443e3f41 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -80,6 +80,9 @@ protected: //- Initial chemical time step const scalar deltaTChemIni_; + //- Maximum chemical time step + const scalar deltaTChemMax_; + //- Latest estimation of integration step volScalarField::Internal deltaTChem_; diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C index f903f7914..e56812a97 100644 --- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C +++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C @@ -387,9 +387,9 @@ jacobian { if (exp < 1.0) { - if (c2[si]>small) + if (c2[si] > small) { - kf *= exp*pow(c2[si] + vSmall, exp - 1.0); + kf *= exp*pow(c2[si], exp - 1.0); } else { diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C index c32221bf4..3eaa07994 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C +++ b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C @@ -125,16 +125,7 @@ Foam::scalar Foam::ReversibleReaction const scalarField& c ) const { - scalar Kc = this->Kc(p, T); - - if (mag(Kc) > vSmall) - { - return kfwd/Kc; - } - else - { - return 0; - } + return kfwd/max(this->Kc(p, T), 1e-6); }