From 9ee3bbc943d4dfcbef766a9ec574d39f8ae1c92e Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 28 Nov 2017 10:23:52 +0000 Subject: [PATCH] integrationSchemes: Corrections to coupled/non-coupled force splitting The integration splitting implemented in commit a5806207 has been shown to be incorrect in some cases. A new procedure has been implemented which can correctly split the implicit-explicit integral into a number of pieces, in order to calculate the contribution of each. This is intended for integrating coupled and non-coupled particle momentum and heat transfers. However, currently there is only ever one implicit coefficient used in these transfers (there is no implicit non-coupled contribution). The evaluation has therefore been short-cutted to only do the integration with respect to the coupled contributions. The splitting functionality has been retained in case additional separate implicit coefficients are required in the future. This change was made with help from Timo Niemi, VTT This resolves bug report https://bugs.openfoam.org/view.php?id=2666 --- .../integrationScheme/Euler/Euler.C | 10 ++++ .../integrationScheme/Euler/Euler.H | 4 ++ .../integrationScheme/analytical/analytical.C | 18 ++++++- .../integrationScheme/analytical/analytical.H | 3 ++ .../integrationScheme/integrationScheme.H | 48 ++++++++++++++----- .../integrationSchemeTemplates.C | 35 ++++++++++---- .../KinematicParcel/KinematicParcel.C | 30 ++++++++++-- .../Templates/ThermoParcel/ThermoParcel.C | 8 ++-- 8 files changed, 125 insertions(+), 31 deletions(-) diff --git a/src/lagrangian/intermediate/integrationScheme/Euler/Euler.C b/src/lagrangian/intermediate/integrationScheme/Euler/Euler.C index a0fa45f7b..3380e1f36 100644 --- a/src/lagrangian/intermediate/integrationScheme/Euler/Euler.C +++ b/src/lagrangian/intermediate/integrationScheme/Euler/Euler.C @@ -61,4 +61,14 @@ Foam::scalar Foam::integrationSchemes::Euler::dtEff } +Foam::scalar Foam::integrationSchemes::Euler::sumDtEff +( + const scalar dt, + const scalar Beta +) const +{ + return sqr(dt)/(1 + Beta*dt); +} + + // ************************************************************************* // diff --git a/src/lagrangian/intermediate/integrationScheme/Euler/Euler.H b/src/lagrangian/intermediate/integrationScheme/Euler/Euler.H index 21c417efa..5ae87cbda 100644 --- a/src/lagrangian/intermediate/integrationScheme/Euler/Euler.H +++ b/src/lagrangian/intermediate/integrationScheme/Euler/Euler.H @@ -79,6 +79,10 @@ public: //- Return the integration effective time step virtual scalar dtEff(const scalar dt, const scalar Beta) const; + + //- Return the integral of the effective time step (using an Euler + // integration method) + virtual scalar sumDtEff(const scalar dt, const scalar Beta) const; }; diff --git a/src/lagrangian/intermediate/integrationScheme/analytical/analytical.C b/src/lagrangian/intermediate/integrationScheme/analytical/analytical.C index ed9362d67..4c9dc9bfd 100644 --- a/src/lagrangian/intermediate/integrationScheme/analytical/analytical.C +++ b/src/lagrangian/intermediate/integrationScheme/analytical/analytical.C @@ -57,7 +57,23 @@ Foam::scalar Foam::integrationSchemes::analytical::dtEff const scalar Beta ) const { - return mag(Beta*dt) > SMALL ? (1 - exp(- Beta*dt))/Beta : dt; + return + mag(Beta*dt) > SMALL + ? (1 - exp(- Beta*dt))/Beta + : dt; +} + + +Foam::scalar Foam::integrationSchemes::analytical::sumDtEff +( + const scalar dt, + const scalar Beta +) const +{ + return + mag(Beta*dt) > SMALL + ? dt/Beta - (1 - exp(- Beta*dt))/sqr(Beta) + : sqr(dt); } diff --git a/src/lagrangian/intermediate/integrationScheme/analytical/analytical.H b/src/lagrangian/intermediate/integrationScheme/analytical/analytical.H index 5bfef1df4..ea1176b4c 100644 --- a/src/lagrangian/intermediate/integrationScheme/analytical/analytical.H +++ b/src/lagrangian/intermediate/integrationScheme/analytical/analytical.H @@ -79,6 +79,9 @@ public: //- Return the integration effective time step virtual scalar dtEff(const scalar dt, const scalar Beta) const; + + //- Return the integral of the effective time step + virtual scalar sumDtEff(const scalar dt, const scalar Beta) const; }; diff --git a/src/lagrangian/intermediate/integrationScheme/integrationScheme/integrationScheme.H b/src/lagrangian/intermediate/integrationScheme/integrationScheme/integrationScheme.H index ead1eb98b..01c41e6a7 100644 --- a/src/lagrangian/intermediate/integrationScheme/integrationScheme/integrationScheme.H +++ b/src/lagrangian/intermediate/integrationScheme/integrationScheme/integrationScheme.H @@ -34,7 +34,8 @@ Description The methods are defined in terms of the effective time-step \f$\Delta t_e\f$ by which the explicit rate is multipled. The effective time-step is - a function of the actual time step and the implicit coefficient. + a function of the actual time step and the implicit coefficient, which must + be implemented in each derived scheme. \f[ \Delta t_e = f(\Delta t, B) @@ -49,14 +50,20 @@ Description 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 + \frac{d \phi_i}{d t} = \alpha_i - \beta_i \phi \f] \f[ - \Delta \phi_i = (\alpha_i - \beta_i \phi^n) \Delta t_e + \Delta \phi_i = \alpha_i \Delta t - + \beta_i \int_0^{\Delta t} \phi d t + \f] + \f[ + \Delta \phi_i = (\alpha_i - \beta_i \phi^n) \Delta t - + (A - B \phi^n) \int_0^{\Delta t} t_e dt \f] - Note that the effective time step \f$\Delta t_e\f$ is not split up. + These partial calculations are defined in terms of the integral of the + effective time-step, \f$\int_0^{\Delta t} t_e dt\f$, which is also + implemented in every derivation. SourceFiles integrationScheme.C @@ -124,9 +131,19 @@ public: // Member Functions - //- Perform the full integration + //- Perform the integration explicitly template - Type Delta + static Type explicitDelta + ( + const Type& phi, + const scalar dtEff, + const Type& Alpha, + const scalar Beta + ); + + //- Perform the integration + template + Type delta ( const Type& phi, const scalar dt, @@ -134,18 +151,23 @@ public: const scalar Beta ) const; - //- Perform a stage of the integration + //- Perform a part of the integration template - static Type delta + Type partialDelta ( const Type& phi, - const scalar dtEff, - const Type& alpha, - const scalar beta - ); + const scalar dt, + const Type& Alpha, + const scalar Beta, + const Type& alphai, + const scalar betai + ) const; //- Return the integration effective time step virtual scalar dtEff(const scalar dt, const scalar Beta) const = 0; + + //- Return the integral of the effective time step + virtual scalar sumDtEff(const scalar dt, const scalar Beta) const = 0; }; diff --git a/src/lagrangian/intermediate/integrationScheme/integrationScheme/integrationSchemeTemplates.C b/src/lagrangian/intermediate/integrationScheme/integrationScheme/integrationSchemeTemplates.C index 7202f2c12..0c7365ba0 100644 --- a/src/lagrangian/intermediate/integrationScheme/integrationScheme/integrationSchemeTemplates.C +++ b/src/lagrangian/intermediate/integrationScheme/integrationScheme/integrationSchemeTemplates.C @@ -28,15 +28,15 @@ License // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -inline Type Foam::integrationScheme::Delta +inline Type Foam::integrationScheme::explicitDelta ( const Type& phi, - const scalar dt, + const scalar dtEff, const Type& Alpha, const scalar Beta -) const +) { - return delta(phi, dtEff(dt, Beta), Alpha, Beta); + return (Alpha - Beta*phi)*dtEff; } @@ -44,12 +44,29 @@ template inline Type Foam::integrationScheme::delta ( const Type& phi, - const scalar dtEff, - const Type& alpha, - const scalar beta -) + const scalar dt, + const Type& Alpha, + const scalar Beta +) const { - return (alpha - beta*phi)*dtEff; + return explicitDelta(phi, dtEff(dt, Beta), Alpha, Beta); +} + + +template +inline Type Foam::integrationScheme::partialDelta +( + const Type& phi, + const scalar dt, + const Type& Alpha, + const scalar Beta, + const Type& alphai, + const scalar betai +) const +{ + return + explicitDelta(phi, dt, alphai, betai) + - explicitDelta(phi, sumDtEff(dt, Beta), Alpha, Beta)*betai; } diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index e654b56d3..26206c454 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -179,6 +179,8 @@ const Foam::vector Foam::KinematicParcel::calcVelocity const forceSuSp Fncp = forces.calcNonCoupled(p, ttd, dt, mass, Re, mu); const scalar massEff = forces.massEff(p, ttd, mass); + /* + // Proper splitting ... // Calculate the integration coefficients const vector acp = (Fcp.Sp()*td.Uc() + Fcp.Su())/massEff; const vector ancp = (Fncp.Sp()*td.Uc() + Fncp.Su() + Su)/massEff; @@ -186,12 +188,32 @@ const Foam::vector Foam::KinematicParcel::calcVelocity const scalar bncp = Fncp.Sp()/massEff; // Integrate to find the new parcel velocity - const scalar dtEff = cloud.UIntegrator().dtEff(dt, bcp + bncp); - const vector deltaUcp = integrationScheme::delta(U_, dtEff, acp, bcp); - const vector deltaUncp = integrationScheme::delta(U_, dtEff, ancp, bncp); + const vector deltaUcp = + cloud.UIntegrator().partialDelta + ( + U_, dt, acp + ancp, bcp + bncp, acp, bcp + ); + const vector deltaUncp = + cloud.UIntegrator().partialDelta + ( + U_, dt, acp + ancp, bcp + bncp, ancp, bncp + ); + const vector deltaT = deltaUcp + deltaUncp; + */ + + // Shortcut splitting assuming no implicit non-coupled force ... + // Calculate the integration coefficients + const vector acp = (Fcp.Sp()*td.Uc() + Fcp.Su())/massEff; + const vector ancp = (Fncp.Su() + Su)/massEff; + const scalar bcp = Fcp.Sp()/massEff; + + // Integrate to find the new parcel velocity + const vector deltaU = cloud.UIntegrator().delta(U_, dt, acp + ancp, bcp); + const vector deltaUncp = ancp*dt; + const vector deltaUcp = deltaU - deltaUncp; // Calculate the new velocity and the momentum transfer terms - vector Unew = U_ + deltaUcp + deltaUncp; + vector Unew = U_ + deltaU; dUTrans -= massEff*deltaUcp; diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C index d44526e75..5d0dfec41 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C @@ -286,12 +286,12 @@ Foam::scalar Foam::ThermoParcel::calcHeatTransfer ancp /= m*Cp_; // Integrate to find the new parcel temperature - const scalar dtEff = cloud.TIntegrator().dtEff(dt, bcp); - const scalar deltaTcp = integrationScheme::delta(T_, dtEff, acp, bcp); - const scalar deltaTncp = integrationScheme::delta(T_, dtEff, ancp, 0); + const scalar deltaT = cloud.TIntegrator().delta(T_, dt, acp + ancp, bcp); + const scalar deltaTncp = ancp*dt; + const scalar deltaTcp = deltaT - deltaTncp; // Calculate the new temperature and the enthalpy transfer terms - scalar Tnew = T_ + deltaTcp + deltaTncp; + scalar Tnew = T_ + deltaT; Tnew = min(max(Tnew, cloud.constProps().TMin()), cloud.constProps().TMax()); dhsTrans -= m*Cp_*deltaTcp;