From 2d0f459ea5a33eafa6fe423d5ff6395c20a38635 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Wed, 31 Jul 2019 15:40:04 +0100 Subject: [PATCH] chemistryModel, reactions, ODESolver: Added "li" argument to all functions which evaluate the reaction rate In chemistryModel "li" is set to the current cell index but for other reacting systems it should be set to the current index of the element for which the reaction system is being evaluated. In the ODESolver "li" is the current index of the element for which the ODE system is being solved if there is a list of related systems being solved, otherwise it can be set to 0. --- src/ODE/ODESolvers/Euler/Euler.C | 6 +- src/ODE/ODESolvers/Euler/Euler.H | 2 + src/ODE/ODESolvers/EulerSI/EulerSI.C | 8 ++- src/ODE/ODESolvers/EulerSI/EulerSI.H | 2 + src/ODE/ODESolvers/ODESolver/ODESolver.C | 11 ++-- src/ODE/ODESolvers/ODESolver/ODESolver.H | 25 ++++--- src/ODE/ODESolvers/RKCK45/RKCK45.C | 16 +++-- src/ODE/ODESolvers/RKCK45/RKCK45.H | 2 + src/ODE/ODESolvers/RKDP45/RKDP45.C | 18 ++--- src/ODE/ODESolvers/RKDP45/RKDP45.H | 2 + src/ODE/ODESolvers/RKF45/RKF45.C | 16 +++-- src/ODE/ODESolvers/RKF45/RKF45.H | 4 +- .../ODESolvers/Rosenbrock12/Rosenbrock12.C | 10 +-- .../ODESolvers/Rosenbrock12/Rosenbrock12.H | 2 + .../ODESolvers/Rosenbrock23/Rosenbrock23.C | 10 +-- .../ODESolvers/Rosenbrock23/Rosenbrock23.H | 2 + .../ODESolvers/Rosenbrock34/Rosenbrock34.C | 12 ++-- .../ODESolvers/Rosenbrock34/Rosenbrock34.H | 2 + src/ODE/ODESolvers/SIBS/SIBS.C | 9 +-- src/ODE/ODESolvers/SIBS/SIBS.H | 2 + src/ODE/ODESolvers/SIBS/SIMPR.C | 7 +- src/ODE/ODESolvers/Trapezoid/Trapezoid.C | 8 ++- src/ODE/ODESolvers/Trapezoid/Trapezoid.H | 2 + .../adaptiveSolver/adaptiveSolver.C | 7 +- .../adaptiveSolver/adaptiveSolver.H | 2 + src/ODE/ODESolvers/rodas23/rodas23.C | 12 ++-- src/ODE/ODESolvers/rodas23/rodas23.H | 2 + src/ODE/ODESolvers/rodas34/rodas34.C | 18 ++--- src/ODE/ODESolvers/rodas34/rodas34.H | 2 + src/ODE/ODESolvers/seulex/seulex.C | 16 +++-- src/ODE/ODESolvers/seulex/seulex.H | 2 + src/ODE/ODESystem/ODESystem.H | 10 ++- .../StandardChemistryModel.C | 33 ++++++---- .../StandardChemistryModel.H | 17 +++-- .../TDACChemistryModel/TDACChemistryModel.C | 35 ++++++---- .../TDACChemistryModel/TDACChemistryModel.H | 17 +++-- .../TDACChemistryModel/reduction/DAC/DAC.C | 10 +-- .../TDACChemistryModel/reduction/DAC/DAC.H | 5 +- .../TDACChemistryModel/reduction/DRG/DRG.C | 12 ++-- .../TDACChemistryModel/reduction/DRG/DRG.H | 5 +- .../reduction/DRGEP/DRGEP.C | 11 ++-- .../reduction/DRGEP/DRGEP.H | 5 +- .../TDACChemistryModel/reduction/EFA/EFA.C | 10 +-- .../TDACChemistryModel/reduction/EFA/EFA.H | 5 +- .../TDACChemistryModel/reduction/PFA/PFA.C | 9 +-- .../TDACChemistryModel/reduction/PFA/PFA.H | 5 +- .../chemistryReductionMethod.H | 7 +- .../noChemistryReduction.C | 7 +- .../noChemistryReduction.H | 7 +- .../TDACChemistryModel/tabulation/ISAT/ISAT.C | 8 ++- .../TDACChemistryModel/tabulation/ISAT/ISAT.H | 2 + .../chemistryTabulationMethod.H | 3 +- .../noChemistryTabulation.H | 1 + .../EulerImplicit/EulerImplicit.C | 13 ++-- .../EulerImplicit/EulerImplicit.H | 5 +- .../chemistrySolver/chemistrySolver.H | 7 +- .../noChemistrySolver/noChemistrySolver.C | 7 +- .../noChemistrySolver/noChemistrySolver.H | 7 +- .../chemistryModel/chemistrySolver/ode/ode.C | 9 +-- .../chemistryModel/chemistrySolver/ode/ode.H | 5 +- .../pyrolysisChemistryModel.C | 66 +++++++++---------- .../pyrolysisChemistryModel.H | 25 +++---- .../solidChemistryModel/solidChemistryModel.H | 22 ++++--- .../solidArrheniusReactionRate.H | 10 ++- .../solidArrheniusReactionRateI.H | 10 ++- .../IrreversibleReaction.C | 27 +++++--- .../IrreversibleReaction.H | 17 +++-- .../NonEquilibriumReversibleReaction.C | 33 ++++++---- .../NonEquilibriumReversibleReaction.H | 17 +++-- .../reaction/Reactions/Reaction/Reaction.C | 26 +++++--- .../reaction/Reactions/Reaction/Reaction.H | 23 +++++-- .../Reactions/ReactionProxy/ReactionProxy.C | 17 +++-- .../Reactions/ReactionProxy/ReactionProxy.H | 17 +++-- .../ReversibleReaction/ReversibleReaction.C | 29 ++++---- .../ReversibleReaction/ReversibleReaction.H | 17 +++-- .../ArrheniusReactionRate.H | 10 ++- .../ArrheniusReactionRateI.H | 10 ++- .../ChemicallyActivatedReactionRate.H | 10 ++- .../ChemicallyActivatedReactionRateI.H | 32 +++++---- .../FallOffReactionRate/FallOffReactionRate.H | 10 ++- .../FallOffReactionRateI.H | 32 +++++---- .../JanevReactionRate/JanevReactionRate.H | 10 ++- .../JanevReactionRate/JanevReactionRateI.H | 12 ++-- .../LandauTellerReactionRate.H | 10 ++- .../LandauTellerReactionRateI.H | 10 ++- .../LangmuirHinshelwoodReactionRate.H | 10 ++- .../LangmuirHinshelwoodReactionRateI.H | 10 ++- .../MichaelisMentenReactionRate.H | 10 ++- .../MichaelisMentenReactionRateI.H | 10 ++- .../infiniteReactionRate.H | 12 ++-- .../infiniteReactionRateI.H | 12 ++-- .../powerSeries/powerSeriesReactionRate.H | 10 ++- .../powerSeries/powerSeriesReactionRateI.H | 10 ++- .../thirdBodyArrheniusReactionRate.H | 10 ++- .../thirdBodyArrheniusReactionRateI.H | 16 +++-- 95 files changed, 713 insertions(+), 415 deletions(-) diff --git a/src/ODE/ODESolvers/Euler/Euler.C b/src/ODE/ODESolvers/Euler/Euler.C index 883928b73a..bf1f2a1b3f 100644 --- a/src/ODE/ODESolvers/Euler/Euler.C +++ b/src/ODE/ODESolvers/Euler/Euler.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -68,6 +68,7 @@ Foam::scalar Foam::Euler::solve ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y @@ -93,10 +94,11 @@ void Foam::Euler::solve ( scalar& x, scalarField& y, + const label li, scalar& dxTry ) const { - adaptiveSolver::solve(odes_, x, y, dxTry); + adaptiveSolver::solve(odes_, x, y, li, dxTry); } diff --git a/src/ODE/ODESolvers/Euler/Euler.H b/src/ODE/ODESolvers/Euler/Euler.H index 179cd83d18..4821d9255e 100644 --- a/src/ODE/ODESolvers/Euler/Euler.H +++ b/src/ODE/ODESolvers/Euler/Euler.H @@ -97,6 +97,7 @@ public: ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y @@ -107,6 +108,7 @@ public: ( scalar& x, scalarField& y, + const label li, scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/EulerSI/EulerSI.C b/src/ODE/ODESolvers/EulerSI/EulerSI.C index c1773083a4..f4dcc466cf 100644 --- a/src/ODE/ODESolvers/EulerSI/EulerSI.C +++ b/src/ODE/ODESolvers/EulerSI/EulerSI.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -78,12 +78,13 @@ Foam::scalar Foam::EulerSI::solve ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y ) const { - odes_.jacobian(x0, y0, dfdx_, dfdy_); + odes_.jacobian(x0, y0, li, dfdx_, dfdy_); for (label i=0; i= 0) diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H index cdf54061d0..3450c3b19d 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.H +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.H @@ -183,35 +183,42 @@ public: inline void resizeMatrix(scalarSquareMatrix& m) const; - //- Solve the ODE system as far as possible up to dxTry - // adjusting the step as necessary to provide a solution within - // the specified tolerance. + //- Solve the ODE system from the current state xStart, y + // and the optional index into the list of systems to solve li + // as far as possible up to dxTry adjusting the step as necessary + // to provide a solution within the specified tolerance. // Update the state and return an estimate for the next step in dxTry virtual void solve ( scalar& x, scalarField& y, + const label li, scalar& dxTry ) const; - //- Solve the ODE system as far as possible up to dxTry - // adjusting the step as necessary to provide a solution within - // the specified tolerance. - // Update the state and return an estimate for the next step in dxTry + //- Solve the ODE system from the current state xStart, y + // and the optional index into the list of systems to solve li + // as far as possible up to step.dxTry adjusting the step as necessary + // to provide a solution within the specified tolerance. + // Update the state and return the step actually taken in step.dxDid + // and an estimate for the next step in step.dxTry virtual void solve ( scalar& x, scalarField& y, + const label li, stepState& step ) const; - //- Solve the ODE system from xStart to xEnd, update the state - // and return an estimate for the next step in dxTry + //- Solve the ODE system from the current state xStart, y + // and the optional index into the list of systems to solve li + // to xEnd and return an estimate for the next step in dxTry virtual void solve ( const scalar xStart, const scalar xEnd, scalarField& y, + const label li, scalar& dxEst ) const; diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.C b/src/ODE/ODESolvers/RKCK45/RKCK45.C index 7c9881b0fa..5c591e3b1d 100644 --- a/src/ODE/ODESolvers/RKCK45/RKCK45.C +++ b/src/ODE/ODESolvers/RKCK45/RKCK45.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -114,6 +114,7 @@ Foam::scalar Foam::RKCK45::solve ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y @@ -124,21 +125,21 @@ Foam::scalar Foam::RKCK45::solve yTemp_[i] = y0[i] + a21*dx*dydx0[i]; } - odes_.derivatives(x0 + c2*dx, yTemp_, k2_); + odes_.derivatives(x0 + c2*dx, yTemp_, li, k2_); forAll(yTemp_, i) { yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]); } - odes_.derivatives(x0 + c3*dx, yTemp_, k3_); + odes_.derivatives(x0 + c3*dx, yTemp_, li, k3_); forAll(yTemp_, i) { yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]); } - odes_.derivatives(x0 + c4*dx, yTemp_, k4_); + odes_.derivatives(x0 + c4*dx, yTemp_, li, k4_); forAll(yTemp_, i) { @@ -146,7 +147,7 @@ Foam::scalar Foam::RKCK45::solve + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]); } - odes_.derivatives(x0 + c5*dx, yTemp_, k5_); + odes_.derivatives(x0 + c5*dx, yTemp_, li, k5_); forAll(yTemp_, i) { @@ -155,7 +156,7 @@ Foam::scalar Foam::RKCK45::solve *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]); } - odes_.derivatives(x0 + c6*dx, yTemp_, k6_); + odes_.derivatives(x0 + c6*dx, yTemp_, li, k6_); forAll(y, i) { @@ -178,10 +179,11 @@ void Foam::RKCK45::solve ( scalar& x, scalarField& y, + const label li, scalar& dxTry ) const { - adaptiveSolver::solve(odes_, x, y, dxTry); + adaptiveSolver::solve(odes_, x, y, li, dxTry); } diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.H b/src/ODE/ODESolvers/RKCK45/RKCK45.H index 17666c8e61..a6a25db712 100644 --- a/src/ODE/ODESolvers/RKCK45/RKCK45.H +++ b/src/ODE/ODESolvers/RKCK45/RKCK45.H @@ -117,6 +117,7 @@ public: ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y @@ -127,6 +128,7 @@ public: ( scalar& x, scalarField& y, + const label li, scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.C b/src/ODE/ODESolvers/RKDP45/RKDP45.C index d43550aed3..604f24c679 100644 --- a/src/ODE/ODESolvers/RKDP45/RKDP45.C +++ b/src/ODE/ODESolvers/RKDP45/RKDP45.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -119,6 +119,7 @@ Foam::scalar Foam::RKDP45::solve ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y @@ -129,21 +130,21 @@ Foam::scalar Foam::RKDP45::solve yTemp_[i] = y0[i] + a21*dx*dydx0[i]; } - odes_.derivatives(x0 + c2*dx, yTemp_, k2_); + odes_.derivatives(x0 + c2*dx, yTemp_, li, k2_); forAll(yTemp_, i) { yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]); } - odes_.derivatives(x0 + c3*dx, yTemp_, k3_); + odes_.derivatives(x0 + c3*dx, yTemp_, li, k3_); forAll(yTemp_, i) { yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]); } - odes_.derivatives(x0 + c4*dx, yTemp_, k4_); + odes_.derivatives(x0 + c4*dx, yTemp_, li, k4_); forAll(yTemp_, i) { @@ -151,7 +152,7 @@ Foam::scalar Foam::RKDP45::solve + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]); } - odes_.derivatives(x0 + c5*dx, yTemp_, k5_); + odes_.derivatives(x0 + c5*dx, yTemp_, li, k5_); forAll(yTemp_, i) { @@ -160,7 +161,7 @@ Foam::scalar Foam::RKDP45::solve *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]); } - odes_.derivatives(x0 + dx, yTemp_, k6_); + odes_.derivatives(x0 + dx, yTemp_, li, k6_); forAll(y, i) { @@ -169,7 +170,7 @@ Foam::scalar Foam::RKDP45::solve } // Reuse k2_ for the derivative of the new state - odes_.derivatives(x0 + dx, y, k2_); + odes_.derivatives(x0 + dx, y, li, k2_); forAll(err_, i) { @@ -187,10 +188,11 @@ void Foam::RKDP45::solve ( scalar& x, scalarField& y, + const label li, scalar& dxTry ) const { - adaptiveSolver::solve(odes_, x, y, dxTry); + adaptiveSolver::solve(odes_, x, y, li, dxTry); } diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.H b/src/ODE/ODESolvers/RKDP45/RKDP45.H index 6068de2b66..a8e7aa59df 100644 --- a/src/ODE/ODESolvers/RKDP45/RKDP45.H +++ b/src/ODE/ODESolvers/RKDP45/RKDP45.H @@ -120,6 +120,7 @@ public: ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y @@ -130,6 +131,7 @@ public: ( scalar& x, scalarField& y, + const label li, scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/RKF45/RKF45.C b/src/ODE/ODESolvers/RKF45/RKF45.C index fae956314a..a6f72e3bd5 100644 --- a/src/ODE/ODESolvers/RKF45/RKF45.C +++ b/src/ODE/ODESolvers/RKF45/RKF45.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -115,6 +115,7 @@ Foam::scalar Foam::RKF45::solve ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y @@ -125,21 +126,21 @@ Foam::scalar Foam::RKF45::solve yTemp_[i] = y0[i] + a21*dx*dydx0[i]; } - odes_.derivatives(x0 + c2*dx, yTemp_, k2_); + odes_.derivatives(x0 + c2*dx, yTemp_, li, k2_); forAll(yTemp_, i) { yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]); } - odes_.derivatives(x0 + c3*dx, yTemp_, k3_); + odes_.derivatives(x0 + c3*dx, yTemp_, li, k3_); forAll(yTemp_, i) { yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]); } - odes_.derivatives(x0 + c4*dx, yTemp_, k4_); + odes_.derivatives(x0 + c4*dx, yTemp_, li, k4_); forAll(yTemp_, i) { @@ -147,7 +148,7 @@ Foam::scalar Foam::RKF45::solve + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]); } - odes_.derivatives(x0 + c5*dx, yTemp_, k5_); + odes_.derivatives(x0 + c5*dx, yTemp_, li, k5_); forAll(yTemp_, i) { @@ -156,7 +157,7 @@ Foam::scalar Foam::RKF45::solve *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]); } - odes_.derivatives(x0 + c6*dx, yTemp_, k6_); + odes_.derivatives(x0 + c6*dx, yTemp_, li, k6_); // Calculate the 5th-order solution forAll(y, i) @@ -183,10 +184,11 @@ void Foam::RKF45::solve ( scalar& x, scalarField& y, + const label li, scalar& dxTry ) const { - adaptiveSolver::solve(odes_, x, y, dxTry); + adaptiveSolver::solve(odes_, x, y, li, dxTry); } diff --git a/src/ODE/ODESolvers/RKF45/RKF45.H b/src/ODE/ODESolvers/RKF45/RKF45.H index aa4f579a70..eae458378d 100644 --- a/src/ODE/ODESolvers/RKF45/RKF45.H +++ b/src/ODE/ODESolvers/RKF45/RKF45.H @@ -1,4 +1,4 @@ -/*---------------------------------------------------------------------------*\ + /*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org @@ -122,6 +122,7 @@ public: ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y @@ -132,6 +133,7 @@ public: ( scalar& x, scalarField& y, + const label li, scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C index ab6b48a113..8e309b4ce6 100644 --- a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C +++ b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -94,12 +94,13 @@ Foam::scalar Foam::Rosenbrock12::solve ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y ) const { - odes_.jacobian(x0, y0, dfdx_, dfdy_); + odes_.jacobian(x0, y0, li, dfdx_, dfdy_); for (label i=0; i 1) diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H index 27ac8bdedd..77433c727e 100644 --- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H +++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H @@ -82,6 +82,7 @@ public: ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y @@ -93,6 +94,7 @@ public: const ODESystem& ode, scalar& x, scalarField& y, + const label li, scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/rodas23/rodas23.C b/src/ODE/ODESolvers/rodas23/rodas23.C index deff872201..f2c532620f 100644 --- a/src/ODE/ODESolvers/rodas23/rodas23.C +++ b/src/ODE/ODESolvers/rodas23/rodas23.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -100,12 +100,13 @@ Foam::scalar Foam::rodas23::solve ( const scalar x0, const scalarField& y0, + const label li, const scalarField& dydx0, const scalar dx, scalarField& y ) const { - odes_.jacobian(x0, y0, dfdx_, dfdy_); + odes_.jacobian(x0, y0, li, dfdx_, dfdy_); for (label i=0; i jacRedo_) { - odes_.jacobian(x, y, dfdx_, dfdy_); + odes_.jacobian(x, y, li, dfdx_, dfdy_); jacUpdated = true; } @@ -299,7 +301,7 @@ void Foam::seulex::solve for (k=0; k<=kTarg_+1; k++) { - bool success = seul(x, y0_, dx, k, ySequence_, scale_); + bool success = seul(x, y0_, li, dx, k, ySequence_, scale_); if (!success) { @@ -427,7 +429,7 @@ void Foam::seulex::solve if (theta_ > jacRedo_ && !jacUpdated) { - odes_.jacobian(x, y, dfdx_, dfdy_); + odes_.jacobian(x, y, li, dfdx_, dfdy_); jacUpdated = true; } } diff --git a/src/ODE/ODESolvers/seulex/seulex.H b/src/ODE/ODESolvers/seulex/seulex.H index 497b3b5326..f6e138503a 100644 --- a/src/ODE/ODESolvers/seulex/seulex.H +++ b/src/ODE/ODESolvers/seulex/seulex.H @@ -108,6 +108,7 @@ class seulex ( const scalar x0, const scalarField& y0, + const label li, const scalar dxTot, const label k, scalarField& y, @@ -150,6 +151,7 @@ public: ( scalar& x, scalarField& y, + const label li, stepState& step ) const; }; diff --git a/src/ODE/ODESystem/ODESystem.H b/src/ODE/ODESystem/ODESystem.H index c7390d2894..2d7dd68057 100644 --- a/src/ODE/ODESystem/ODESystem.H +++ b/src/ODE/ODESystem/ODESystem.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -67,19 +67,25 @@ public: virtual label nEqns() const = 0; //- Calculate the derivatives in dydx + // for the current state x and y + // and optional index into the list of systems to solve li virtual void derivatives ( const scalar x, const scalarField& y, + const label li, scalarField& dydx ) const = 0; //- Calculate the Jacobian of the system - // Need by the stiff-system solvers + // for the current state x and y + // and optional index into the list of systems to solve li. + // Need by stiff-system solvers virtual void jacobian ( const scalar x, const scalarField& y, + const label li, scalarField& dfdx, scalarSquareMatrix& dfdy ) const = 0; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C index 7ad38cd52d..bf1d52ed2d 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C @@ -103,9 +103,10 @@ Foam::StandardChemistryModel:: template void Foam::StandardChemistryModel::omega ( - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalarField& dcdt ) const { @@ -116,7 +117,7 @@ void Foam::StandardChemistryModel::omega { const Reaction& R = reactions_[i]; - R.omega(p, T, c, dcdt); + R.omega(p, T, c, li, dcdt); } } @@ -125,9 +126,10 @@ template Foam::scalar Foam::StandardChemistryModel::omegaI ( const label index, - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -137,7 +139,7 @@ Foam::scalar Foam::StandardChemistryModel::omegaI ) const { const Reaction& R = reactions_[index]; - scalar w = R.omega(p, T, c, pf, cf, lRef, pr, cr, rRef); + scalar w = R.omega(p, T, c, li, pf, cf, lRef, pr, cr, rRef); return(w); } @@ -147,6 +149,7 @@ void Foam::StandardChemistryModel::derivatives ( const scalar time, const scalarField& c, + const label li, scalarField& dcdt ) const { @@ -158,7 +161,7 @@ void Foam::StandardChemistryModel::derivatives c_[i] = max(c[i], 0); } - omega(c_, T, p, dcdt); + omega(p, T, c_, li, dcdt); // Constant pressure // dT/dt = ... @@ -197,6 +200,7 @@ void Foam::StandardChemistryModel::jacobian ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt, scalarSquareMatrix& J ) const @@ -227,8 +231,8 @@ void Foam::StandardChemistryModel::jacobian { const Reaction& R = reactions_[ri]; scalar kfwd, kbwd; - R.dwdc(p, T, c_, J, dcdt, omegaI, kfwd, kbwd, false, dummy); - R.dwdT(p, T, c_, omegaI, kfwd, kbwd, J, false, dummy, nSpecie_); + R.dwdc(p, T, c_, li, J, dcdt, omegaI, kfwd, kbwd, false, dummy); + R.dwdT(p, T, c_, li, omegaI, kfwd, kbwd, J, false, dummy, nSpecie_); } // The species derivatives of the temperature term are partially computed @@ -318,7 +322,7 @@ Foam::StandardChemistryModel::tc() const { const Reaction& R = reactions_[i]; - R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef); + R.omega(pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef); forAll(R.rhs(), s) { @@ -410,7 +414,10 @@ Foam::StandardChemistryModel::calculateRR } const Reaction& R = reactions_[ri]; - const scalar omegai = R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef); + const scalar omegai = R.omega + ( + pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef + ); forAll(R.lhs(), s) { @@ -461,7 +468,7 @@ void Foam::StandardChemistryModel::calculate() c_[i] = rhoi*Yi/specieThermo_[i].W(); } - omega(c_, Ti, pi, dcdt_); + omega(pi, Ti, c_, celli, dcdt_); for (label i=0; i::solve while (timeLeft > small) { scalar dt = timeLeft; - this->solve(c_, Ti, pi, dt, this->deltaTChem_[celli]); + this->solve(pi, Ti, c_, celli, dt, this->deltaTChem_[celli]); timeLeft -= dt; } diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H index 28091d7d43..98b0935910 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H @@ -154,9 +154,10 @@ public: //- dc/dt = omega, rate of change in concentration, for each species virtual void omega ( - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalarField& dcdt ) const; @@ -166,9 +167,10 @@ public: virtual scalar omegaI ( label iReaction, - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -227,6 +229,7 @@ public: ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt ) const; @@ -234,15 +237,17 @@ public: ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt, scalarSquareMatrix& J ) const; virtual void solve ( - scalarField &c, - scalar& T, scalar& p, + scalar& T, + scalarField& c, + const label li, scalar& deltaT, scalar& subDeltaT ) const = 0; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C index 986b8af0b9..a874fb7003 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C @@ -149,9 +149,10 @@ Foam::TDACChemistryModel::~TDACChemistryModel() template void Foam::TDACChemistryModel::omega ( - const scalarField& c, // Contains all species even when mechRed is active - const scalar T, const scalar p, + const scalar T, + const scalarField& c, // Contains all species even when mechRed is active + const label li, scalarField& dcdt ) const { @@ -170,7 +171,7 @@ void Foam::TDACChemistryModel::omega scalar omegai = R.omega ( - p, T, c, pf, cf, lRef, pr, cr, rRef + p, T, c, li, pf, cf, lRef, pr, cr, rRef ); forAll(R.lhs(), s) @@ -204,9 +205,10 @@ template Foam::scalar Foam::TDACChemistryModel::omega ( const Reaction& R, - const scalarField& c, // Contains all species even when mechRed is active - const scalar T, const scalar p, + const scalar T, + const scalarField& c, // Contains all species even when mechRed is active + const label li, scalar& pf, scalar& cf, label& lRef, @@ -215,8 +217,8 @@ Foam::scalar Foam::TDACChemistryModel::omega label& rRef ) const { - const scalar kf = R.kf(p, T, c); - const scalar kr = R.kr(kf, p, T, c); + const scalar kf = R.kf(p, T, c, li); + const scalar kr = R.kr(kf, p, T, c, li); const label Nl = R.lhs().size(); const label Nr = R.rhs().size(); @@ -314,6 +316,7 @@ void Foam::TDACChemistryModel::derivatives ( const scalar time, const scalarField& c, + const label li, scalarField& dcdt ) const { @@ -345,7 +348,7 @@ void Foam::TDACChemistryModel::derivatives } } - omega(this->c_, T, p, dcdt); + omega(p, T, this->c_, li, dcdt); // Constant pressure // dT/dt = ... @@ -398,6 +401,7 @@ void Foam::TDACChemistryModel::jacobian ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt, scalarSquareMatrix& J ) const @@ -451,6 +455,7 @@ void Foam::TDACChemistryModel::jacobian p, T, this->c_, + li, J, dcdt, omegaI, @@ -464,6 +469,7 @@ void Foam::TDACChemistryModel::jacobian p, T, this->c_, + li, omegaI, kfwd, kbwd, @@ -651,7 +657,7 @@ Foam::scalar Foam::TDACChemistryModel::solve if (reduced) { // Reduce mechanism change the number of species (only active) - mechRed_->reduceMechanism(c, Ti, pi); + mechRed_->reduceMechanism(pi, Ti, c, celli); nActiveSpecies += mechRed_->NsSimp(); nAvg++; scalar timeIncr = clockTime_.timeIncrement(); @@ -672,7 +678,12 @@ Foam::scalar Foam::TDACChemistryModel::solve // Solve the reduced set of ODE this->solve ( - simplifiedC_, Ti, pi, dt, this->deltaTChem_[celli] + pi, + Ti, + simplifiedC_, + celli, + dt, + this->deltaTChem_[celli] ); for (label i=0; i::solve } else { - this->solve(c, Ti, pi, dt, this->deltaTChem_[celli]); + this->solve(pi, Ti, c, celli, dt, this->deltaTChem_[celli]); } timeLeft -= dt; } @@ -713,7 +724,7 @@ Foam::scalar Foam::TDACChemistryModel::solve Rphiq[Rphiq.size()-1] = pi; } label growOrAdd = - tabulation_->add(phiq, Rphiq, rhoi, deltaT[celli]); + tabulation_->add(phiq, Rphiq, celli, rhoi, deltaT[celli]); if (growOrAdd) { this->setTabulationResultsAdd(celli); diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H index a1935d4d82..30400754d0 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H @@ -173,9 +173,10 @@ public: //- dc/dt = omega, rate of change in concentration, for each species virtual void omega ( - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalarField& dcdt ) const; @@ -184,9 +185,10 @@ public: virtual scalar omega ( const Reaction& r, - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -215,6 +217,7 @@ public: ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt ) const; @@ -222,15 +225,17 @@ public: ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt, scalarSquareMatrix& J ) const; virtual void solve ( - scalarField& c, - scalar& T, scalar& p, + scalar& T, + scalarField& c, + const label li, scalar& deltaT, scalar& subDeltaT ) const = 0; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C index 50435dbcfe..b0d9e55144 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -270,9 +270,10 @@ Foam::chemistryReductionMethods::DAC::~DAC() template void Foam::chemistryReductionMethods::DAC::reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ) { scalarField& completeC(this->chemistry_.completeC()); @@ -303,10 +304,11 @@ void Foam::chemistryReductionMethods::DAC::reduceMechanism forAll(this->chemistry_.reactions(), i) { const Reaction& R = this->chemistry_.reactions()[i]; + // for each reaction compute omegai scalar omegai = this->chemistry_.omega ( - R, c1, T, p, pf, cf, lRef, pr, cr, rRef + R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef ); // Then for each pair of species composing this reaction, diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.H index e173c5d5c4..ff9603698b 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.H @@ -147,9 +147,10 @@ public: //- Reduce the mechanism virtual void reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ); }; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C index 894de70016..3cb832af6b 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -69,9 +69,10 @@ Foam::chemistryReductionMethods::DRG::~DRG() template void Foam::chemistryReductionMethods::DRG::reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ) { scalarField c1(this->nSpecie_+2, 0.0); @@ -102,11 +103,12 @@ void Foam::chemistryReductionMethods::DRG::reduceMechanism forAll(this->chemistry_.reactions(), i) { const Reaction& R = this->chemistry_.reactions()[i]; + // For each reaction compute omegai scalar omegai = this->chemistry_.omega ( - R, c1, T, p, pf, cf, lRef, pr, cr, rRef - ); + R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef + ); // Then for each pair of species composing this reaction, diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.H index 7b42184061..c92c98da23 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.H @@ -83,9 +83,10 @@ public: //- Reduce the mechanism virtual void reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ); }; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C index 6e520b6ae5..a500d7c733 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -115,9 +115,10 @@ template void Foam::chemistryReductionMethods::DRGEP:: reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ) { scalarField& completeC(this->chemistry_.completeC()); @@ -153,8 +154,8 @@ reduceMechanism // for each reaction compute omegai scalar omegai = this->chemistry_.omega ( - R, c1, T, p, pf, cf, lRef, pr, cr, rRef - ); + R, p,T, c1, li, pf, cf, lRef, pr, cr, rRef + ); omegaV[i] = omegai; // then for each pair of species composing this reaction, diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.H index a806bd9f0c..267cd624d1 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.H @@ -154,9 +154,10 @@ public: //- Reduce the mechanism virtual void reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ); }; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C index 18ad092deb..f7a13446c7 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -94,9 +94,10 @@ Foam::chemistryReductionMethods::EFA::~EFA() template void Foam::chemistryReductionMethods::EFA::reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ) { scalarField& completeC(this->chemistry_.completeC()); @@ -132,10 +133,11 @@ void Foam::chemistryReductionMethods::EFA::reduceMechanism forAll(this->chemistry_.reactions(), i) { const Reaction& R = this->chemistry_.reactions()[i]; + // for each reaction compute omegai this->chemistry_.omega ( - R, c1, T, p, pf, cf, lRef, pr, cr, rRef + R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef ); scalar fr = mag(pf*cf)+mag(pr*cr); scalar NCi(0.0),NHi(0.0),NOi(0.0),NNi(0.0); diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.H index a21010dc60..0702a11f14 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.H @@ -80,9 +80,10 @@ public: //- Reduce the mechanism virtual void reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ); }; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C index 492eeeb8cf..f27abc902b 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -72,9 +72,10 @@ Foam::chemistryReductionMethods::PFA::~PFA() template void Foam::chemistryReductionMethods::PFA::reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ) { scalarField& completeC(this->chemistry_.completeC()); @@ -110,7 +111,7 @@ void Foam::chemistryReductionMethods::PFA::reduceMechanism // for each reaction compute omegai scalar omegai = this->chemistry_.omega ( - R, c1, T, p, pf, cf, lRef, pr, cr, rRef + R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef ); // then for each pair of species composing this reaction, diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.H index 50da03ddd3..1daec070e3 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.H @@ -82,9 +82,10 @@ public: //- Reduce the mechanism virtual void reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ); }; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H index 1efde6e03b..b4cf137eb3 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -149,9 +149,10 @@ public: //- Reduce the mechanism virtual void reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ) = 0; }; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.C index 719b6fb801..2ae43317b6 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -53,9 +53,10 @@ template void Foam::chemistryReductionMethods::none:: reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ) { NotImplemented; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.H index 234fee7cee..dc7a372462 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -76,9 +76,10 @@ public: //- Reduce the mechanism virtual void reduceMechanism ( - const scalarField &c, + const scalar p, const scalar T, - const scalar p + const scalarField& c, + const label li ); }; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C index f04c3dc47e..fe99195ca8 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -341,6 +341,7 @@ void Foam::chemistryTabulationMethods::ISAT::computeA ( scalarSquareMatrix& A, const scalarField& Rphiq, + const label li, const scalar rhoi, const scalar dt ) @@ -374,7 +375,7 @@ void Foam::chemistryTabulationMethods::ISAT::computeA // A = C(psi0,t0)/(I-dt*J(psi(t0+dt))) // where C(psi0,t0) = I scalarField dcdt(speciesNumber + 2, Zero); - this->chemistry_.jacobian(runTime_.value(), Rcq, dcdt, A); + this->chemistry_.jacobian(runTime_.value(), Rcq, li, dcdt, A); // The jacobian is computed according to the molar concentration // the following conversion allows the code to use A with mass fraction @@ -536,6 +537,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT::add ( const scalarField& phiq, const scalarField& Rphiq, + const label li, const scalar rho, const scalar deltaT ) @@ -614,7 +616,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT::add // Compute the A matrix needed to store the chemPoint. label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2; scalarSquareMatrix A(ASize, Zero); - computeA(A, Rphiq, rho, deltaT); + computeA(A, Rphiq, li, rho, deltaT); chemisTree().insertNewLeaf ( diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.H index c3d6dee564..06c899ae8b 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.H @@ -168,6 +168,7 @@ class ISAT ( scalarSquareMatrix& A, const scalarField& Rphiq, + const label li, const scalar rho, const scalar dt ); @@ -232,6 +233,7 @@ public: ( const scalarField& phiq, const scalarField& Rphiq, + const label li, const scalar rho, const scalar deltaT ); diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H index a9de9182ba..f04b31430b 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2016-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -161,6 +161,7 @@ public: ( const scalarField& phiQ, const scalarField& RphiQ, + const label li, const scalar rho, const scalar deltaT ) = 0; diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H index d51c8aefba..3e7d746101 100644 --- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H +++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H @@ -105,6 +105,7 @@ public: ( const scalarField& phiq, const scalarField& Rphiq, + const label li, const scalar rho, const scalar deltaT ) diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C index 3a0898787b..18124e9947 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -91,9 +91,10 @@ void Foam::EulerImplicit::updateRRInReactionI template void Foam::EulerImplicit::solve ( - scalarField& c, - scalar& T, scalar& p, + scalar& T, + scalarField& c, + const label li, scalar& deltaT, scalar& subDeltaT ) const @@ -124,8 +125,10 @@ void Foam::EulerImplicit::solve scalar pf, cf, pr, cr; label lRef, rRef; - const scalar omegai = - this->omegaI(i, c, T, p, pf, cf, lRef, pr, cr, rRef); + const scalar omegai + ( + this->omegaI(i, p, T, c, li, pf, cf, lRef, pr, cr, rRef) + ); scalar corr = 1; if (eqRateLimiter_) diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H index 2e3900c5e1..41306c53f8 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H @@ -107,9 +107,10 @@ public: //- Update the concentrations and return the chemical time virtual void solve ( - scalarField& c, - scalar& T, scalar& p, + scalar& T, + scalarField& c, + const label li, scalar& deltaT, scalar& subDeltaT ) const; diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H index ea8b8987c6..722cbb10e7 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -70,9 +70,10 @@ public: //- Update the concentrations and return the chemical time virtual void solve ( - scalarField &c, - scalar& T, scalar& p, + scalar& T, + scalarField& c, + const label li, scalar& deltaT, scalar& subDeltaT ) const = 0; diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C index 24556cdc32..12d582b855 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -50,9 +50,10 @@ Foam::noChemistrySolver::~noChemistrySolver() template void Foam::noChemistrySolver::solve ( + scalar&, + scalar&, scalarField&, - scalar&, - scalar&, + const label li, scalar&, scalar& ) const diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H index 303c7f333f..704b689a17 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,9 +74,10 @@ public: //- Update the concentrations and return the chemical time virtual void solve ( - scalarField& c, - scalar& T, scalar& p, + scalar& T, + scalarField& c, + const label li, scalar& deltaT, scalar& subDeltaT ) const; diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C index 49085bfa14..3ea09388f4 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -49,9 +49,10 @@ Foam::ode::~ode() template void Foam::ode::solve ( - scalarField& c, - scalar& T, scalar& p, + scalar& T, + scalarField& c, + const label li, scalar& deltaT, scalar& subDeltaT ) const @@ -73,7 +74,7 @@ void Foam::ode::solve cTp_[nSpecie] = T; cTp_[nSpecie+1] = p; - odeSolver_->solve(0, deltaT, cTp_, subDeltaT); + odeSolver_->solve(0, deltaT, cTp_, li, subDeltaT); for (int i=0; iYs_.size() + nGases_), RRg_(nGases_), - Ys0_(this->nSolids_), - cellCounter_(0) + Ys0_(this->nSolids_) { // create the fields for the chemistry sources forAll(this->RRs_, fieldi) @@ -184,17 +183,16 @@ template Foam::scalarField Foam:: pyrolysisChemistryModel::omega ( - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, const bool updateC0 ) const { scalar pf, cf, pr, cr; label lRef, rRef; - const label celli = cellCounter_; - scalarField om(nEqns(), 0.0); forAll(this->reactions_, i) @@ -207,7 +205,7 @@ pyrolysisChemistryModel::omega scalar omegai = omega ( - R, c, T, p, pf, cf, lRef, pr, cr, rRef + R, p, T, c, li, pf, cf, lRef, pr, cr, rRef ); scalar rhoL = 0.0; forAll(R.lhs(), s) @@ -226,7 +224,7 @@ pyrolysisChemistryModel::omega if (updateC0) { - Ys0_[si][celli] += sr*omegai; + Ys0_[si][li] += sr*omegai; } } forAll(R.grhs(), g) @@ -245,9 +243,10 @@ Foam::scalar Foam::pyrolysisChemistryModel::omega ( const Reaction& R, - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -258,14 +257,12 @@ Foam::pyrolysisChemistryModel::omega { scalarField c1(nSpecie_, 0.0); - label celli = cellCounter_; - for (label i=0; i::omega const scalar exp = R.lhs()[si].exponent; kf *= - pow(c1[si]/Ys0_[si][celli], exp) - *(Ys0_[si][celli]); + pow(c1[si]/Ys0_[si][li], exp) + *(Ys0_[si][li]); } return kf; @@ -288,9 +285,10 @@ Foam::scalar Foam::pyrolysisChemistryModel:: omegaI ( const label index, - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -301,7 +299,7 @@ omegaI { const Reaction& R = this->reactions_[index]; - scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef); + scalar w = omega(R, p, T, c, li, pf, cf, lRef, pr, cr, rRef); return(w); } @@ -311,7 +309,8 @@ void Foam::pyrolysisChemistryModel:: derivatives ( const scalar time, - const scalarField &c, + const scalarField& c, + const label li, scalarField& dcdt ) const { @@ -320,7 +319,7 @@ derivatives dcdt = 0.0; - dcdt = omega(c, T, p); + dcdt = omega(p, T, c, li); // Total mass concentration scalar cTot = 0.0; @@ -354,6 +353,7 @@ jacobian ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt, scalarSquareMatrix& dfdc ) const @@ -377,13 +377,13 @@ jacobian } // length of the first argument must be nSolids - dcdt = omega(c2, T, p); + dcdt = omega(p, T, c2, li); for (label ri=0; rireactions_.size(); ri++) { const Reaction& R = this->reactions_[ri]; - scalar kf0 = R.kf(p, T, c2); + scalar kf0 = R.kf(p, T, c2, li); forAll(R.lhs(), j) { @@ -434,14 +434,13 @@ jacobian // calculate the dcdT elements numerically scalar delta = 1.0e-8; - scalarField dcdT0 = omega(c2, T - delta, p); - scalarField dcdT1 = omega(c2, T + delta, p); + scalarField dcdT0 = omega(p , T - delta, c2, li); + scalarField dcdT1 = omega(p, T + delta, c2, li); for (label i=0; imesh().V()[celli]; if (this->reactingCells_[celli]) @@ -505,7 +502,7 @@ calculate() c[i] = rhoi*this->Ys_[i][celli]*delta; } - const scalarField dcdt = omega(c, Ti, pi, true); + const scalarField dcdt = omega(pi, Ti, c, celli, true); forAll(this->RRs_, i) { @@ -570,8 +567,6 @@ Foam::pyrolysisChemistryModel::solve { if (this->reactingCells_[celli]) { - cellCounter_ = celli; - scalar rhoi = rho[celli]; scalar pi = p[celli]; scalar Ti = T[celli]; @@ -590,7 +585,7 @@ Foam::pyrolysisChemistryModel::solve while (timeLeft > small) { scalar dt = timeLeft; - this->solve(c, Ti, pi, dt, this->deltaTChem_[celli]); + this->solve(pi, Ti, c, celli, dt, this->deltaTChem_[celli]); timeLeft -= dt; } @@ -608,7 +603,7 @@ Foam::pyrolysisChemistryModel::solve } // Update Ys0_ - dc = omega(c0, Ti, pi, true); + dc = omega(pi, Ti, c0, celli, true); } } @@ -654,9 +649,10 @@ Foam::pyrolysisChemistryModel::gasHs template void Foam::pyrolysisChemistryModel::solve ( - scalarField &c, - scalar& T, scalar& p, + scalar& T, + scalarField& c, + const label li, scalar& deltaT, scalar& subDeltaT ) const diff --git a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H index a63f7f2cf3..e0c3f91c1c 100644 --- a/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H +++ b/src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H @@ -87,9 +87,6 @@ private: //- List of accumulative solid concentrations mutable PtrList Ys0_; - //- Cell counter - label cellCounter_; - public: @@ -125,9 +122,10 @@ public: //- dc/dt = omega, rate of change in concentration, for each species virtual scalarField omega ( - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, const bool updateC0 = false ) const; @@ -137,9 +135,10 @@ public: virtual scalar omega ( const Reaction& r, - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -155,9 +154,10 @@ public: virtual scalar omegaI ( label iReaction, - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -204,6 +204,7 @@ public: ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt ) const; @@ -211,15 +212,17 @@ public: ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt, scalarSquareMatrix& dfdc ) const; virtual void solve ( - scalarField &c, - scalar& T, scalar& p, + scalar& T, + scalarField& c, + const label li, scalar& deltaT, scalar& subDeltaT ) const; diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H index de42392156..44b95483bf 100644 --- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H +++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H @@ -126,9 +126,10 @@ public: //- dc/dt = omega, rate of change in concentration, for each species virtual scalarField omega ( - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, const bool updateC0 = false ) const = 0; @@ -137,9 +138,10 @@ public: virtual scalar omega ( const Reaction& r, - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -154,9 +156,10 @@ public: virtual scalar omegaI ( label iReaction, - const scalarField& c, - const scalar T, const scalar p, + const scalar T, + const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -204,6 +207,7 @@ public: ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt ) const = 0; @@ -211,15 +215,17 @@ public: ( const scalar t, const scalarField& c, + const label li, scalarField& dcdt, scalarSquareMatrix& dfdc ) const = 0; virtual void solve ( - scalarField &c, - scalar& T, scalar& p, + scalar& T, + scalarField& c, + const label li, scalar& deltaT, scalar& subDeltaT ) const = 0; diff --git a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H index dcab54cb81..f6c64f3f41 100644 --- a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H +++ b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H @@ -101,14 +101,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -123,6 +125,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -132,7 +135,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; diff --git a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H index a9dd9d0795..0c9035e3e2 100644 --- a/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H +++ b/src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H @@ -57,7 +57,8 @@ inline Foam::scalar Foam::solidArrheniusReactionRate::operator() ( const scalar, const scalar T, - const scalarField& + const scalarField&, + const label li ) const { if (T < Tcrit_) @@ -75,7 +76,8 @@ inline Foam::scalar Foam::solidArrheniusReactionRate::ddT ( const scalar p, const scalar T, - const scalarField& + const scalarField&, + const label li ) const { if (T < Tcrit_) @@ -101,6 +103,7 @@ inline void Foam::solidArrheniusReactionRate::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const {} @@ -110,7 +113,8 @@ inline Foam::scalar Foam::solidArrheniusReactionRate::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; diff --git a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C index eb77b4fa5b..8eb9fa851e 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C +++ b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -99,10 +99,11 @@ Foam::scalar Foam::IrreversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return k_(p, T, c); + return k_(p, T, c, li); } @@ -122,7 +123,8 @@ Foam::scalar Foam::IrreversibleReaction const scalar kfwd, const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; @@ -144,7 +146,8 @@ Foam::scalar Foam::IrreversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; @@ -166,10 +169,11 @@ Foam::scalar Foam::IrreversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return k_.ddT(p, T, c); + return k_.ddT(p, T, c, li); } @@ -189,6 +193,7 @@ Foam::scalar Foam::IrreversibleReaction const scalar p, const scalar T, const scalarField& c, + const label li, const scalar dkfdT, const scalar kr ) const @@ -231,10 +236,11 @@ void Foam::IrreversibleReaction const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const { - k_.dcidc(p, T, c, dcidc); + k_.dcidc(p, T, c, li, dcidc); } @@ -253,10 +259,11 @@ Foam::scalar Foam::IrreversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return k_.dcidT(p, T, c); + return k_.dcidT(p, T, c, li); } diff --git a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H index 0d8473ad7c..34c8586bbe 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H +++ b/src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H @@ -146,7 +146,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Reverse rate constant from the given forward rate constant @@ -156,7 +157,8 @@ public: const scalar kfwd, const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Reverse rate constant @@ -165,7 +167,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; @@ -176,7 +179,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Temperature derivative of reverse rate @@ -186,6 +190,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, const scalar dkfdT, const scalar kr ) const; @@ -202,6 +207,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -211,7 +217,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; diff --git a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C index 87002858d0..228b940f91 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C +++ b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -124,10 +124,11 @@ Foam::NonEquilibriumReversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return fk_(p, T, c); + return fk_(p, T, c, li); } @@ -148,10 +149,11 @@ Foam::NonEquilibriumReversibleReaction const scalar, const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return rk_(p, T, c); + return rk_(p, T, c, li); } @@ -171,10 +173,11 @@ Foam::NonEquilibriumReversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return rk_(p, T, c); + return rk_(p, T, c, li); } @@ -194,10 +197,11 @@ Foam::NonEquilibriumReversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return fk_.ddT(p, T, c); + return fk_.ddT(p, T, c, li); } @@ -218,11 +222,12 @@ Foam::NonEquilibriumReversibleReaction const scalar p, const scalar T, const scalarField& c, + const label li, const scalar dkfdT, const scalar kr ) const { - return rk_.ddT(p, T, c); + return rk_.ddT(p, T, c, li); } @@ -260,10 +265,11 @@ void Foam::NonEquilibriumReversibleReaction const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const { - fk_.dcidc(p, T, c, dcidc); + fk_.dcidc(p, T, c, li, dcidc); } @@ -282,10 +288,11 @@ Foam::scalar Foam::NonEquilibriumReversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return fk_.dcidT(p, T, c); + return fk_.dcidT(p, T, c, li); } diff --git a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H index ae107ce08a..7c3b36ed4b 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H +++ b/src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H @@ -134,7 +134,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Reverse rate constant from the given formard rate constant @@ -143,7 +144,8 @@ public: const scalar kfwd, const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Reverse rate constant. @@ -153,7 +155,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; @@ -164,7 +167,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Temperature derivative of backward rate @@ -173,6 +177,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, const scalar dkfdT, const scalar kr ) const; @@ -189,6 +194,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -198,7 +204,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C index 02ba9e9561..1c75b240e7 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C +++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C @@ -205,6 +205,7 @@ void Foam::Reaction::ddot const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& d ) const { @@ -217,6 +218,7 @@ void Foam::Reaction::fdot const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& f ) const { @@ -229,6 +231,7 @@ void Foam::Reaction::omega const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcdt ) const { @@ -237,7 +240,7 @@ void Foam::Reaction::omega scalar omegaI = omega ( - p, T, c, pf, cf, lRef, pr, cr, rRef + p, T, c, li, pf, cf, lRef, pr, cr, rRef ); forAll(lhs_, i) @@ -261,6 +264,7 @@ Foam::scalar Foam::Reaction::omega const scalar p, const scalar T, const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -272,8 +276,8 @@ Foam::scalar Foam::Reaction::omega scalar clippedT = min(max(T, this->Tlow()), this->Thigh()); - const scalar kf = this->kf(p, clippedT, c); - const scalar kr = this->kr(kf, p, clippedT, c); + const scalar kf = this->kf(p, clippedT, c, li); + const scalar kr = this->kr(kf, p, clippedT, c, li); pf = 1; pr = 1; @@ -375,6 +379,7 @@ void Foam::Reaction::dwdc const scalar p, const scalar T, const scalarField& c, + const label li, scalarSquareMatrix& J, scalarField& dcdt, scalar& omegaI, @@ -387,7 +392,7 @@ void Foam::Reaction::dwdc scalar pf, cf, pr, cr; label lRef, rRef; - omegaI = omega(p, T, c, pf, cf, lRef, pr, cr, rRef); + omegaI = omega(p, T, c, li, pf, cf, lRef, pr, cr, rRef); forAll(lhs_, i) { @@ -402,8 +407,8 @@ void Foam::Reaction::dwdc dcdt[si] += sr*omegaI; } - kfwd = this->kf(p, T, c); - kbwd = this->kr(kfwd, p, T, c); + kfwd = this->kf(p, T, c, li); + kbwd = this->kr(kfwd, p, T, c, li); forAll(lhs_, j) { @@ -504,7 +509,7 @@ void Foam::Reaction::dwdc { // This temporary array needs to be cached for efficiency scalarField dcidc(beta.size()); - this->dcidc(p, T, c, dcidc); + this->dcidc(p, T, c, li, dcidc); forAll(beta, j) { @@ -538,6 +543,7 @@ void Foam::Reaction::dwdT const scalar p, const scalar T, const scalarField& c, + const label li, const scalar omegaI, const scalar kfwd, const scalar kbwd, @@ -550,8 +556,8 @@ void Foam::Reaction::dwdT scalar kf = kfwd; scalar kr = kbwd; - scalar dkfdT = this->dkfdT(p, T, c); - scalar dkrdT = this->dkrdT(p, T, c, dkfdT, kr); + scalar dkfdT = this->dkfdT(p, T, c, li); + scalar dkrdT = this->dkrdT(p, T, c, li, dkfdT, kr); scalar sumExp = 0.0; forAll(lhs_, i) @@ -582,7 +588,7 @@ void Foam::Reaction::dwdT // For reactions including third-body efficiencies or pressure dependent // reaction, an additional term is needed - scalar dcidT = this->dcidT(p, T, c); + scalar dcidT = this->dcidT(p, T, c, li); dcidT *= omegaI; // J(i, indexT) = sum_reactions nu_i dqdT diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H index 92beddb174..33b1d7fa7e 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H +++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H @@ -210,6 +210,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& d ) const; @@ -219,6 +220,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& f ) const; @@ -228,6 +230,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcdt ) const; @@ -237,6 +240,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalar& pf, scalar& cf, label& lRef, @@ -252,7 +256,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const = 0; //- Reverse rate constant from the given forward rate constant @@ -261,7 +266,8 @@ public: const scalar kfwd, const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const = 0; //- Reverse rate constant @@ -269,7 +275,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const = 0; @@ -282,6 +289,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarSquareMatrix& J, scalarField& dcdt, scalar& omegaI, @@ -298,6 +306,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, const scalar omegaI, const scalar kfwd, const scalar kbwd, @@ -312,7 +321,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const = 0; //- Temperature derivative of reverse rate @@ -321,6 +331,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, const scalar dkfdT, const scalar kr ) const = 0; @@ -336,6 +347,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const = 0; @@ -344,7 +356,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const = 0; diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C index 4002e4af76..4a26822ba7 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C +++ b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C @@ -106,7 +106,8 @@ Foam::scalar Foam::ReactionProxy::kf ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { NotImplemented; @@ -120,7 +121,8 @@ Foam::scalar Foam::ReactionProxy::kr const scalar kfwd, const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { NotImplemented; @@ -133,7 +135,8 @@ Foam::scalar Foam::ReactionProxy::kr ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { NotImplemented; @@ -146,7 +149,8 @@ Foam::scalar Foam::ReactionProxy::dkfdT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { NotImplemented; @@ -160,6 +164,7 @@ Foam::scalar Foam::ReactionProxy::dkrdT const scalar p, const scalar T, const scalarField& c, + const label li, const scalar dkfdT, const scalar kr ) const @@ -184,6 +189,7 @@ void Foam::ReactionProxy::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const { @@ -196,7 +202,8 @@ Foam::scalar Foam::ReactionProxy::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { NotImplemented; diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H index 6ee6dd249d..61fb2d0b63 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H +++ b/src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H @@ -107,7 +107,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Reverse rate constant from the given forward rate constant @@ -116,7 +117,8 @@ public: const scalar kfwd, const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Reverse rate constant @@ -124,7 +126,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; @@ -135,7 +138,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Temperature derivative of reverse rate @@ -144,6 +148,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, const scalar dkfdT, const scalar kr ) const; @@ -159,6 +164,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -167,7 +173,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; }; diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C index cf0db201ec..d850c6940d 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C +++ b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -99,10 +99,11 @@ Foam::scalar Foam::ReversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return k_(p, T, c); + return k_(p, T, c, li); } @@ -122,7 +123,8 @@ Foam::scalar Foam::ReversibleReaction const scalar kfwd, const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return kfwd/max(this->Kc(p, T), rootSmall); @@ -144,10 +146,11 @@ Foam::scalar Foam::ReversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return kr(kf(p, T, c), p, T, c); + return kr(kf(p, T, c, li), p, T, c, li); } @@ -166,10 +169,11 @@ Foam::scalar Foam::ReversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return k_.ddT(p, T, c); + return k_.ddT(p, T, c, li); } @@ -189,6 +193,7 @@ Foam::scalar Foam::ReversibleReaction const scalar p, const scalar T, const scalarField& c, + const label li, const scalar dkfdT, const scalar kr ) const @@ -233,10 +238,11 @@ void Foam::ReversibleReaction const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const { - k_.dcidc(p, T, c, dcidc); + k_.dcidc(p, T, c, li, dcidc); } @@ -255,10 +261,11 @@ Foam::scalar Foam::ReversibleReaction ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - return k_.dcidT(p, T, c); + return k_.dcidT(p, T, c, li); } diff --git a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H index f923eacdd2..b80ca64449 100644 --- a/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H +++ b/src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H @@ -143,7 +143,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Reverse rate constant from the given formard rate constant @@ -152,7 +153,8 @@ public: const scalar kfwd, const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Reverse rate constant. @@ -162,7 +164,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; @@ -173,7 +176,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Temperature derivative of backward rate @@ -182,6 +186,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, const scalar dkfdT, const scalar kr ) const; @@ -198,6 +203,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -207,7 +213,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H index 5d22494103..021d68cf1c 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H @@ -97,14 +97,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -119,6 +121,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -128,7 +131,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Write to stream diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H index 1af507b9f0..698f5b1394 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H @@ -56,7 +56,8 @@ inline Foam::scalar Foam::ArrheniusReactionRate::operator() ( const scalar p, const scalar T, - const scalarField& + const scalarField&, + const label ) const { scalar ak = A_; @@ -79,7 +80,8 @@ inline Foam::scalar Foam::ArrheniusReactionRate::ddT ( const scalar p, const scalar T, - const scalarField& + const scalarField&, + const label ) const { scalar ak = A_; @@ -110,6 +112,7 @@ inline void Foam::ArrheniusReactionRate::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const {} @@ -119,7 +122,8 @@ inline Foam::scalar Foam::ArrheniusReactionRate::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H index ef93894d25..ea9bffb0e7 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H @@ -107,14 +107,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -132,6 +134,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -141,7 +144,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Write to stream diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H index 175ac1cc67..782fbfe239 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -91,11 +91,12 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - const scalar k0 = k0_(p, T, c); - const scalar kInf = kInf_(p, T, c); + const scalar k0 = k0_(p, T, c, li); + const scalar kInf = kInf_(p, T, c, li); const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf; return k0*(1/(1 + Pr))*F_(T, Pr); @@ -111,14 +112,15 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - const scalar k0 = k0_(p, T, c); - const scalar kInf = kInf_(p, T, c); + const scalar k0 = k0_(p, T, c, li); + const scalar kInf = kInf_(p, T, c, li); const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf; - return (1/(1 + Pr))*F_(T, Pr)*k0_.ddT(p, T, c); + return (1/(1 + Pr))*F_(T, Pr)*k0_.ddT(p, T, c, li); } @@ -132,6 +134,7 @@ inline void Foam::ChemicallyActivatedReactionRate const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const { @@ -139,8 +142,8 @@ inline void Foam::ChemicallyActivatedReactionRate if (M > small) { - const scalar k0 = k0_(p, T, c); - const scalar kInf = kInf_(p, T, c); + const scalar k0 = k0_(p, T, c, li); + const scalar kInf = kInf_(p, T, c, li); const scalar Pr = k0*M/kInf; const scalar F = F_(T, Pr); @@ -171,20 +174,21 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { const scalar M = thirdBodyEfficiencies_.M(c); if (M > small) { - const scalar k0 = k0_(p, T, c); - const scalar kInf = kInf_(p, T, c); + const scalar k0 = k0_(p, T, c, li); + const scalar kInf = kInf_(p, T, c, li); const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf; const scalar F = F_(T, Pr); const scalar dPrdT = - Pr*(k0_.ddT(p, T, c)/k0 - kInf_.ddT(p, T, c)/kInf - 1/T); + Pr*(k0_.ddT(p, T, c, li)/k0 - kInf_.ddT(p, T, c, li)/kInf - 1/T); const scalar dFdT = F_.ddT(Pr, F, dPrdT, T); return (-dPrdT/(1 + Pr) + dFdT/F); diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H index 2b10cd5078..a92afd38a0 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H @@ -104,14 +104,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -128,6 +130,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -136,7 +139,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Write to stream diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H index ea6d49dd1b..d07410c0f3 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -82,11 +82,12 @@ Foam::FallOffReactionRate::operator() ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - const scalar k0 = k0_(p, T, c); - const scalar kInf = kInf_(p, T, c); + const scalar k0 = k0_(p, T, c, li); + const scalar kInf = kInf_(p, T, c, li); const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf; return kInf*(Pr/(1 + Pr))*F_(T, Pr); @@ -99,14 +100,15 @@ Foam::FallOffReactionRate::ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { - const scalar k0 = k0_(p, T, c); - const scalar kInf = kInf_(p, T, c); + const scalar k0 = k0_(p, T, c, li); + const scalar kInf = kInf_(p, T, c, li); const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf; - return (Pr/(1 + Pr))*F_(T, Pr)*kInf_.ddT(p, T, c); + return (Pr/(1 + Pr))*F_(T, Pr)*kInf_.ddT(p, T, c, li); } @@ -116,6 +118,7 @@ inline void Foam::FallOffReactionRate::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const { @@ -123,8 +126,8 @@ inline void Foam::FallOffReactionRate::dcidc if (M > small) { - const scalar k0 = k0_(p, T, c); - const scalar kInf = kInf_(p, T, c); + const scalar k0 = k0_(p, T, c, li); + const scalar kInf = kInf_(p, T, c, li); const scalar Pr = k0*M/kInf; const scalar F = F_(T, Pr); @@ -151,20 +154,21 @@ Foam::FallOffReactionRate::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { const scalar M = thirdBodyEfficiencies_.M(c); if (M > small) { - const scalar k0 = k0_(p, T, c); - const scalar kInf = kInf_(p, T, c); + const scalar k0 = k0_(p, T, c, li); + const scalar kInf = kInf_(p, T, c, li); const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf; const scalar F = F_(T, Pr); const scalar dPrdT = - Pr*(k0_.ddT(p, T, c)/k0 - kInf_.ddT(p, T, c)/kInf - 1/T); + Pr*(k0_.ddT(p, T, c, li)/k0 - kInf_.ddT(p, T, c, li)/kInf - 1/T); const scalar dFdT = F_.ddT(Pr, F, dPrdT, T); return (dPrdT/(Pr*(1 + Pr)) + dFdT/F); diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H index d8d8b067c6..5678cead58 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H @@ -100,14 +100,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -121,6 +123,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -129,7 +132,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Write to stream diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H index 22db3fa181..82f4309ad3 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -59,7 +59,8 @@ inline Foam::scalar Foam::JanevReactionRate::operator() ( const scalar p, const scalar T, - const scalarField& + const scalarField&, + const label ) const { scalar lta = A_; @@ -93,7 +94,8 @@ inline Foam::scalar Foam::JanevReactionRate::ddT ( const scalar p, const scalar T, - const scalarField& + const scalarField&, + const label ) const { scalar lta = A_; @@ -142,6 +144,7 @@ inline void Foam::JanevReactionRate::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const {} @@ -151,7 +154,8 @@ inline Foam::scalar Foam::JanevReactionRate::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H index fb70d20c76..299ab83d17 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H @@ -99,14 +99,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -120,6 +122,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -128,7 +131,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Write to stream diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H index d7f7e3d7dd..967ab2ed36 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H @@ -62,7 +62,8 @@ inline Foam::scalar Foam::LandauTellerReactionRate::operator() ( const scalar p, const scalar T, - const scalarField& + const scalarField&, + const label ) const { scalar lta = A_; @@ -102,7 +103,8 @@ inline Foam::scalar Foam::LandauTellerReactionRate::ddT ( const scalar p, const scalar T, - const scalarField& + const scalarField&, + const label ) const { scalar lta = A_; @@ -157,6 +159,7 @@ inline void Foam::LandauTellerReactionRate::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const {} @@ -166,7 +169,8 @@ inline Foam::scalar Foam::LandauTellerReactionRate::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H index d98419b3fc..f59536e425 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H @@ -102,14 +102,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -123,6 +125,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -131,7 +134,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Write to stream diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H index a6b2b20c26..7849f750c2 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H @@ -68,7 +68,8 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::operator() ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { const scalar c0m = pow(c[r_[0]], m_[0]); @@ -90,7 +91,8 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { const scalar c0m = pow(c[r_[0]], m_[0]); @@ -125,6 +127,7 @@ inline void Foam::LangmuirHinshelwoodReactionRate::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const {} @@ -134,7 +137,8 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H index 21f9efde34..add5efc2a3 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H @@ -104,14 +104,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -125,6 +127,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -133,7 +136,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Write to stream diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H index 7e0bdd655c..8327972a80 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H @@ -46,7 +46,8 @@ inline Foam::scalar Foam::MichaelisMentenReactionRate::operator() ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return Vmax_/(Km_ + c[s_]); @@ -57,7 +58,8 @@ inline Foam::scalar Foam::MichaelisMentenReactionRate::ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; @@ -76,6 +78,7 @@ inline void Foam::MichaelisMentenReactionRate::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const { @@ -87,7 +90,8 @@ inline Foam::scalar Foam::MichaelisMentenReactionRate::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H index 19affef813..aaaab6dcd2 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -84,14 +84,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -105,6 +107,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -113,7 +116,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Write to stream diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H index 01d89958ba..e79e2e3b73 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -47,7 +47,8 @@ inline Foam::scalar Foam::infiniteReactionRate::operator() ( const scalar p, const scalar, - const scalarField& + const scalarField&, + const label ) const { return (1); @@ -57,7 +58,8 @@ inline Foam::scalar Foam::infiniteReactionRate::ddT ( const scalar p, const scalar, - const scalarField& + const scalarField&, + const label ) const { return (0); @@ -76,6 +78,7 @@ inline void Foam::infiniteReactionRate::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const {} @@ -85,7 +88,8 @@ inline Foam::scalar Foam::infiniteReactionRate::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H index d666d9a4b1..4a3ecaa959 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H @@ -100,14 +100,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -121,6 +123,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -129,7 +132,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Write to stream diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H index c0fe1fae47..2d845a0d52 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H @@ -59,7 +59,8 @@ inline Foam::scalar Foam::powerSeriesReactionRate::operator() ( const scalar p, const scalar T, - const scalarField& + const scalarField&, + const label ) const { scalar lta = A_; @@ -86,7 +87,8 @@ inline Foam::scalar Foam::powerSeriesReactionRate::ddT ( const scalar p, const scalar T, - const scalarField& + const scalarField&, + const label ) const { scalar lta = A_; @@ -124,6 +126,7 @@ inline void Foam::powerSeriesReactionRate::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const {} @@ -133,7 +136,8 @@ inline Foam::scalar Foam::powerSeriesReactionRate::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return 0; diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H index e8026eb0c3..7afc90bb23 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H @@ -97,14 +97,16 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; inline scalar ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Third-body efficiencies (beta = 1-alpha) @@ -121,6 +123,7 @@ public: const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const; @@ -129,7 +132,8 @@ public: ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const; //- Write to stream diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H index 261ff65b94..eca08c3475 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -76,12 +76,13 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::operator() ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return thirdBodyEfficiencies_.M(c) - *ArrheniusReactionRate::operator()(p, T, c); + *ArrheniusReactionRate::operator()(p, T, c, li); } @@ -89,12 +90,13 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::ddT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return thirdBodyEfficiencies_.M(c) - *ArrheniusReactionRate::ddT(p, T, c); + *ArrheniusReactionRate::ddT(p, T, c, li); } @@ -103,6 +105,7 @@ inline void Foam::thirdBodyArrheniusReactionRate::dcidc const scalar p, const scalar T, const scalarField& c, + const label li, scalarField& dcidc ) const { @@ -118,7 +121,8 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::dcidT ( const scalar p, const scalar T, - const scalarField& c + const scalarField& c, + const label li ) const { return -1.0/T;