From a5806207c3a836bf73abfdf455ffc0d43e08a4dc Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Wed, 20 Sep 2017 16:06:48 +0100 Subject: [PATCH] lagrangian: Rewrite of integration schemes and transfer terms The integration of force and heat transfer onto the particle is facilitated by a run-time-selectable integration scheme. These schemes were written to generate the value at the end of an intregration step and also an average value over the step from which the total transfer was computed. The average value in the Euler scheme was implemented incorrectly, which resulted in the momentum and heat transfer processes being non-conservative. Implementing the average correctly, however, would have inteoduced a number of trancendental functions which would have negated the purpose of the Euler scheme as the cheap and stable option. The schemes have been rewritten to generate changes over the step, rather than the final value. This change is then used to calculate the transfers. Regardless of the scheme, this formulation is guaranteed to be conservative, and the Euler scheme remains computationally inexpensive. This change was made with help from Timo Niemi, VTT This resolves bug report https://bugs.openfoam.org/view.php?id=2666 --- .../IntegrationScheme/Analytical/Analytical.C | 44 +----- .../IntegrationScheme/Analytical/Analytical.H | 28 ++-- .../IntegrationScheme/Euler/Euler.C | 32 +--- .../IntegrationScheme/Euler/Euler.H | 28 ++-- .../IntegrationScheme/IntegrationScheme.C | 48 ++++-- .../IntegrationScheme/IntegrationScheme.H | 148 +++++++----------- .../IntegrationScheme/IntegrationSchemeNew.C | 12 +- .../makeIntegrationSchemes.C | 4 +- .../KinematicParcel/KinematicParcel.C | 21 +-- .../Templates/ThermoParcel/ThermoParcel.C | 54 +++---- 10 files changed, 163 insertions(+), 256 deletions(-) diff --git a/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C b/src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C index 5d95ad26e..43b90b5a1 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 84f8c98b6..5326e9699 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 cd483d8e1..be2e65e96 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 efd91eea2..ea8adb214 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 50d05a2ee..4a86a9924 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 88326f8c2..9ddca46bf 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 b85a24c61..5e5cf065c 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 efbd10d2f..e16fca0bf 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 832d0dca9..0d2ee7b2e 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 580a3ac1a..40e5666a4 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; }