diff --git a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C index 5d95ad26ea..43b90b5a17 100644 --- a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C +++ b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -28,59 +28,27 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -Foam::Analytical::Analytical -( - const word& phiName, - const dictionary& dict -) -: - IntegrationScheme(phiName, dict) -{} - - -template -Foam::Analytical::Analytical(const Analytical& is) -: - IntegrationScheme(is) +Foam::integrationSchemes::Analytical::Analytical() {} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template -Foam::Analytical::~Analytical() +Foam::integrationSchemes::Analytical::~Analytical() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -typename Foam::IntegrationScheme::integrationResult -Foam::Analytical::integrate +Foam::scalar Foam::integrationSchemes::Analytical::factor ( - const Type& phi, const scalar dt, - const Type& alphaBeta, - const scalar beta + const scalar Beta ) const { - typename IntegrationScheme::integrationResult retValue; - - const scalar expTerm = exp(min(50, -beta*dt)); - - if (beta > ROOTVSMALL) - { - const Type alpha = alphaBeta/beta; - retValue.average() = alpha + (phi - alpha)*(1 - expTerm)/(beta*dt); - retValue.value() = alpha + (phi - alpha)*expTerm; - } - else - { - retValue.value() = phi + alphaBeta*dt; - retValue.average() = 0.5*(phi + retValue.value()); - } - - return retValue; + return mag(Beta*dt) > SMALL ? (1 - exp(- Beta*dt))/Beta : dt; } diff --git a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H index 84f8c98b6e..5326e9699a 100644 --- a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H +++ b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,11 @@ Class Foam::Analytical Description - Analytical integration + Analytical integration scheme + + \f[ + \Delta \phi = (A - B \phi^n) \frac{1}{B} (1 - e^{- B \Delta t}) + \f] \*---------------------------------------------------------------------------*/ @@ -38,6 +42,8 @@ Description namespace Foam { +namespace integrationSchemes +{ /*---------------------------------------------------------------------------*\ Class Analytical Declaration @@ -56,11 +62,8 @@ public: // Constructors - //- Construct from components - Analytical(const word& phiName, const dictionary& dict); - - //- Copy constructor - Analytical(const Analytical& is); + //- Construct + Analytical(); //- Construct and return clone virtual autoPtr> clone() const @@ -78,19 +81,14 @@ public: // Member Functions - //- Perform the integration - virtual typename IntegrationScheme::integrationResult integrate - ( - const Type& phi, - const scalar dt, - const Type& alphaBeta, - const scalar beta - ) const; + //- Return the integration factor + virtual scalar factor(const scalar dt, const scalar Beta) const; }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +} // End namespace integrationSchemes } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C index cd483d8e11..be2e65e966 100644 --- a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C +++ b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -28,47 +28,27 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -Foam::Euler::Euler -( - const word& phiName, - const dictionary& dict -) -: - IntegrationScheme(phiName, dict) -{} - - -template -Foam::Euler::Euler(const Euler& is) -: - IntegrationScheme(is) +Foam::integrationSchemes::Euler::Euler() {} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template -Foam::Euler::~Euler() +Foam::integrationSchemes::Euler::~Euler() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -typename Foam::IntegrationScheme::integrationResult -Foam::Euler::integrate +Foam::scalar Foam::integrationSchemes::Euler::factor ( - const Type& phi, const scalar dt, - const Type& alphaBeta, - const scalar beta + const scalar Beta ) const { - typename IntegrationScheme::integrationResult retValue; - retValue.value() = (phi + alphaBeta*dt)/(1.0 + beta*dt); - retValue.average() = 0.5*(phi + retValue.value()); - - return retValue; + return dt/(1 + Beta*dt); } diff --git a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H index efd91eea25..ea8adb2144 100644 --- a/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H +++ b/src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,11 @@ Class Foam::Euler Description - Euler-implicit integration + Euler-implicit integration scheme. + + \f[ + \Delta \phi = (A - B \phi^n) \frac{\Delta t}{1 + B \Delta t} + \f] \*---------------------------------------------------------------------------*/ @@ -38,6 +42,8 @@ Description namespace Foam { +namespace integrationSchemes +{ /*---------------------------------------------------------------------------*\ Class Euler Declaration @@ -56,11 +62,8 @@ public: // Constructors - //- Construct from components - Euler(const word& phiName, const dictionary& dict); - - //- Copy constructor - Euler(const Euler& is); + //- Construct + Euler(); //- Construct and return clone virtual autoPtr> clone() const @@ -75,19 +78,14 @@ public: // Member Functions - //- Perform the integration - virtual typename IntegrationScheme::integrationResult integrate - ( - const Type& phi, - const scalar dt, - const Type& alphaBeta, - const scalar beta - ) const; + //- Return the integration factor + virtual scalar factor(const scalar dt, const scalar Beta) const; }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +} // End namespace integrationSchemes } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C index 50d05a2ee3..4a86a99245 100644 --- a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C +++ b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -28,22 +28,7 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -Foam::IntegrationScheme::IntegrationScheme -( - const word& phiName, - const dictionary& dict -) -: - phiName_(phiName), - dict_(dict) -{} - - -template -Foam::IntegrationScheme::IntegrationScheme(const IntegrationScheme& is) -: - phiName_(is.phiName_), - dict_(is.dict_) +Foam::IntegrationScheme::IntegrationScheme() {} @@ -54,6 +39,35 @@ Foam::IntegrationScheme::~IntegrationScheme() {} +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +Type Foam::IntegrationScheme::delta +( + const Type& phi, + const scalar dt, + const Type& Alpha, + const scalar Beta +) const +{ + return delta(phi, dt, Beta, Alpha, Beta); +} + + +template +Type Foam::IntegrationScheme::delta +( + const Type& phi, + const scalar dt, + const scalar Beta, + const Type& alpha, + const scalar beta +) const +{ + return (alpha - beta*phi)*factor(dt, Beta); +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #include "IntegrationSchemeNew.C" diff --git a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H index 88326f8c2b..9ddca46bf4 100644 --- a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H +++ b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,35 @@ Class Foam::IntegrationScheme Description - Top level model for Integration schemes + Base for a set of schemes which integrate simple ODEs which arise from + semi-implcit rate expressions. + + \f[ + \frac{d \phi}{d t} = A - B \phi + \f] + + The methods are defined in terms of the factor \f$f\f$ by which the result + of a trivial explicit integration is multipled. + + \f[ + \Delta \phi = (A - B \phi^n) f(\Delta t, B) + \f] + + This class also facilitates integration in stages. If the explicit and + implicit coefficients, \f$A\f$ and \f$B\f$, are a summation of differing + contributions, \f$\sum \alpha_i\f$ and \f$\sum \beta_i\f$, then the + integration can be split up to detemine the effect of each contribution. + + \f[ + \frac{d \phi}{d t} = \left( \sum \alpha_i \right) - + \left( \sum \beta_i \right) \phi + \f] + \f[ + \Delta \phi_i = (\alpha_i - \beta_i \phi^n) f(\Delta t, B) + \f] + + Note that the sum of implicit coefficients \f$B\f$ is needed for each stage + of the integral, whilst the sum of explicit coefficients \f$A\f$ is not. SourceFiles IntegrationScheme.C @@ -51,79 +79,6 @@ namespace Foam template class IntegrationScheme { - -public: - - //- Helper class to supply results of integration - class integrationResult - { - //- Integration value - Type value_; - - //- Average value across integration step - Type average_; - - - public: - - //- Constructor - integrationResult() - : - value_(Zero), - average_(Zero) - {} - - - // Member functions - - // Access - - //- Return const access to the value - Type value() const - { - return value_; - } - - //- Return const access to the average - Type average() const - { - return average_; - } - - - // Edit - - //- Return access to the value for changing - Type& value() - { - return value_; - } - - //- Return access to the average for changing - Type& average() - { - return average_; - } - }; - - -private: - - // Private data - - //- Name of the Integration variable - const word& phiName_; - - //- Reference to the dictionary - const dictionary& dict_; - - - // Private Member Functions - - //- Disallow default bitwise assignment - void operator=(const IntegrationScheme&); - - public: //- Runtime type information @@ -136,22 +91,16 @@ public: ( autoPtr, IntegrationScheme, - dictionary, - ( - const word& phiName, - const dictionary& dict - ), - (phiName, dict) + word, + (), + () ); // Constructors - //- Construct from components - IntegrationScheme(const word& phiName, const dictionary& dict); - - //- Copy constructor - IntegrationScheme(const IntegrationScheme& is); + //- Construct + IntegrationScheme(); //- Construct and return clone virtual autoPtr> clone() const = 0; @@ -159,7 +108,7 @@ public: // Selectors - //- Return a reference to the selected radiation model + //- Select an integration scheme static autoPtr New ( const word& phiName, @@ -173,14 +122,27 @@ public: // Member Functions - //- Perform the Integration - virtual integrationResult integrate + //- Perform the full integration + Type delta ( const Type& phi, const scalar dt, - const Type& alphaBeta, + const Type& Alpha, + const scalar Beta + ) const; + + //- Perform a stage of the integration + Type delta + ( + const Type& phi, + const scalar dt, + const scalar Beta, + const Type& alpha, const scalar beta - ) const = 0; + ) const; + + //- Return the factor + virtual scalar factor(const scalar dt, const scalar Beta) const = 0; }; @@ -197,7 +159,7 @@ public: defineTemplateRunTimeSelectionTable \ ( \ IntegrationScheme, \ - dictionary \ + word \ ); @@ -205,7 +167,7 @@ public: \ defineNamedTemplateTypeNameAndDebug(SS, 0); \ \ - IntegrationScheme::adddictionaryConstructorToTable> \ + IntegrationScheme::addwordConstructorToTable> \ add##SS##Type##ConstructorToTable_; diff --git a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationSchemeNew.C b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationSchemeNew.C index b85a24c61b..5e5cf065c6 100644 --- a/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationSchemeNew.C +++ b/src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationSchemeNew.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -41,20 +41,20 @@ Foam::IntegrationScheme::New Info<< "Selecting " << phiName << " integration scheme " << schemeName << endl; - typename dictionaryConstructorTable::iterator cstrIter = - dictionaryConstructorTablePtr_->find(schemeName); + typename wordConstructorTable::iterator cstrIter = + wordConstructorTablePtr_->find(schemeName); - if (cstrIter == dictionaryConstructorTablePtr_->end()) + if (cstrIter == wordConstructorTablePtr_->end()) { FatalErrorInFunction << "Unknown integration scheme type " << schemeName << nl << nl << "Valid integration scheme types are:" << nl - << dictionaryConstructorTablePtr_->sortedToc() << nl + << wordConstructorTablePtr_->sortedToc() << nl << exit(FatalError); } - return autoPtr>(cstrIter()(phiName, dict)); + return autoPtr>(cstrIter()()); } // ************************************************************************* // diff --git a/src/lagrangian/intermediate/IntegrationScheme/makeIntegrationSchemes.C b/src/lagrangian/intermediate/IntegrationScheme/makeIntegrationSchemes.C index efbd10d2fe..e16fca0bfd 100644 --- a/src/lagrangian/intermediate/IntegrationScheme/makeIntegrationSchemes.C +++ b/src/lagrangian/intermediate/IntegrationScheme/makeIntegrationSchemes.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -32,6 +32,8 @@ License namespace Foam { + using namespace integrationSchemes; + makeIntegrationScheme(scalar); makeIntegrationSchemeType(Euler, scalar); makeIntegrationSchemeType(Analytical, scalar); diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index 832d0dca91..0d2ee7b2ea 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -177,26 +177,27 @@ const Foam::vector Foam::KinematicParcel::calcVelocity // Momentum source due to particle forces const forceSuSp Fcp = forces.calcCoupled(p, ttd, dt, mass, Re, mu); const forceSuSp Fncp = forces.calcNonCoupled(p, ttd, dt, mass, Re, mu); - const forceSuSp Feff = Fcp + Fncp; const scalar massEff = forces.massEff(p, ttd, mass); // New particle velocity //~~~~~~~~~~~~~~~~~~~~~~ - // Update velocity - treat as 3-D - const vector abp = (Feff.Sp()*td.Uc() + (Feff.Su() + Su))/massEff; - const scalar bp = Feff.Sp()/massEff; + const vector acp = (Fcp.Sp()*td.Uc() + Fcp.Su())/massEff; + const vector ancp = (Fncp.Sp()*td.Uc() + Fncp.Su() + Su)/massEff; + const scalar bcp = Fcp.Sp()/massEff; + const scalar bncp = Fncp.Sp()/massEff; - Spu = dt*Feff.Sp(); + const vector deltaUcp = + cloud.UIntegrator().delta(U_, dt, bcp + bncp, acp, bcp); + const vector deltaUncp = + cloud.UIntegrator().delta(U_, dt, bcp + bncp, ancp, bncp); - IntegrationScheme::integrationResult Ures = - cloud.UIntegrator().integrate(U_, dt, abp, bp); + vector Unew = U_ + deltaUcp + deltaUncp; - vector Unew = Ures.value(); + dUTrans -= massEff*deltaUcp; - // note: Feff.Sp() and Fc.Sp() must be the same - dUTrans += dt*(Feff.Sp()*(Ures.average() - td.Uc()) - Fcp.Su()); + Spu = dt*Fcp.Sp(); // Apply correction to velocity and dUTrans for reduced-D cases const polyMesh& mesh = cloud.pMesh(); diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C index 580a3ac1a9..40e5666a46 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C @@ -263,57 +263,41 @@ Foam::scalar Foam::ThermoParcel::calcHeatTransfer const scalar d = this->d(); const scalar rho = this->rho(); + const scalar As = this->areaS(d); + const scalar V = this->volume(d); + const scalar m = rho*V; // Calc heat transfer coefficient scalar htc = cloud.heatTransfer().htc(d, Re, Pr, kappa, NCpW); - if (mag(htc) < ROOTVSMALL && !cloud.radiation()) - { - return - max - ( - T_ + dt*Sh/(this->volume(d)*rho*Cp_), - cloud.constProps().TMin() - ); - } - - htc = max(htc, ROOTVSMALL); - const scalar As = this->areaS(d); - - scalar ap = td.Tc() + Sh/(As*htc); - const scalar bp = 6.0*htc/max(rho*d*Cp_, ROOTVSMALL); - + // Calculate the integration coefficients + const scalar bcp = htc*As/(m*Cp_); + const scalar acp = bcp*td.Tc(); + scalar ancp = Sh; if (cloud.radiation()) { - tetIndices tetIs = this->currentTetIndices(); + const tetIndices tetIs = this->currentTetIndices(); const scalar Gc = td.GInterp().interpolate(this->coordinates(), tetIs); const scalar sigma = physicoChemical::sigma.value(); const scalar epsilon = cloud.constProps().epsilon0(); - // Assume constant source - scalar s = epsilon*(Gc/4.0 - sigma*pow4(T_)); - - ap += s/htc; + ancp += As*epsilon*(Gc/4.0 - sigma*pow4(T_)); } + ancp /= m*Cp_; // Integrate to find the new parcel temperature - IntegrationScheme::integrationResult Tres = - cloud.TIntegrator().integrate(T_, dt, ap*bp, bp); + const scalar deltaTcp = + cloud.TIntegrator().delta(T_, dt, bcp, acp, bcp); + const scalar deltaTncp = + cloud.TIntegrator().delta(T_, dt, bcp, ancp, 0); - scalar Tnew = - min - ( - max - ( - Tres.value(), - cloud.constProps().TMin() - ), - cloud.constProps().TMax() - ); + // Calculate the new temperature and the transfer terms + scalar Tnew = T_ + deltaTcp + deltaTncp; + Tnew = min(max(Tnew, cloud.constProps().TMin()), cloud.constProps().TMax()); - Sph = dt*htc*As; + Sph = dt*m*Cp_*bcp; - dhsTrans += Sph*(Tres.average() - td.Tc()); + dhsTrans -= m*Cp_*deltaTcp; return Tnew; }