zhyang-dev/OpenFOAM-6
1
0
Fork 0
mirror of https://github.com/OpenFOAM/OpenFOAM-6.git synced 2025-12-08 06:57:46 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

lagrangian: Rewrite of integration schemes and transfer terms

Browse Source 
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
This commit is contained in:
Will Bainbridge
2017-09-20 16:06:48 +01:00
parent 8743c2fa26
commit a5806207c3
10 changed files with 163 additions and 256 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -28,59 +28,27 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type> template<class Type>
Foam::Analytical<Type>::Analytical Foam::integrationSchemes::Analytical<Type>::Analytical()
(
const word& phiName,
const dictionary& dict
)
:
IntegrationScheme<Type>(phiName, dict)
{}
template<class Type>
Foam::Analytical<Type>::Analytical(const Analytical& is)
:
IntegrationScheme<Type>(is)
{} {}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type> template<class Type>
Foam::Analytical<Type>::~Analytical() Foam::integrationSchemes::Analytical<Type>::~Analytical()
{} {}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type> template<class Type>
typename Foam::IntegrationScheme<Type>::integrationResult Foam::scalar Foam::integrationSchemes::Analytical<Type>::factor
Foam::Analytical<Type>::integrate
( (
const Type& phi,
const scalar dt, const scalar dt,
const Type& alphaBeta, const scalar Beta
const scalar beta
) const ) const
{ {
typename IntegrationScheme<Type>::integrationResult retValue; return mag(Beta*dt) > SMALL ? (1 - exp(- Beta*dt))/Beta : dt;
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;
} }

28
src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -25,7 +25,11 @@ Class
Foam::Analytical Foam::Analytical
Description 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 Foam
{ {
namespace integrationSchemes
{
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class Analytical Declaration Class Analytical Declaration
@ -56,11 +62,8 @@ public:
// Constructors // Constructors
//- Construct from components //- Construct
Analytical(const word& phiName, const dictionary& dict); Analytical();
//- Copy constructor
Analytical(const Analytical& is);
//- Construct and return clone //- Construct and return clone
virtual autoPtr<IntegrationScheme<Type>> clone() const virtual autoPtr<IntegrationScheme<Type>> clone() const
@ -78,19 +81,14 @@ public:
// Member Functions // Member Functions
//- Perform the integration //- Return the integration factor
virtual typename IntegrationScheme<Type>::integrationResult integrate virtual scalar factor(const scalar dt, const scalar Beta) const;
(
const Type& phi,
const scalar dt,
const Type& alphaBeta,
const scalar beta
) const;
}; };
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace integrationSchemes
} // End namespace Foam } // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

32
src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -28,47 +28,27 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type> template<class Type>
Foam::Euler<Type>::Euler Foam::integrationSchemes::Euler<Type>::Euler()
(
const word& phiName,
const dictionary& dict
)
:
IntegrationScheme<Type>(phiName, dict)
{}
template<class Type>
Foam::Euler<Type>::Euler(const Euler& is)
:
IntegrationScheme<Type>(is)
{} {}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type> template<class Type>
Foam::Euler<Type>::~Euler() Foam::integrationSchemes::Euler<Type>::~Euler()
{} {}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type> template<class Type>
typename Foam::IntegrationScheme<Type>::integrationResult Foam::scalar Foam::integrationSchemes::Euler<Type>::factor
Foam::Euler<Type>::integrate
( (
const Type& phi,
const scalar dt, const scalar dt,
const Type& alphaBeta, const scalar Beta
const scalar beta
) const ) const
{ {
typename IntegrationScheme<Type>::integrationResult retValue; return dt/(1 + Beta*dt);
retValue.value() = (phi + alphaBeta*dt)/(1.0 + beta*dt);
retValue.average() = 0.5*(phi + retValue.value());
return retValue;
} }

28
src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -25,7 +25,11 @@ Class
Foam::Euler Foam::Euler
Description 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 Foam
{ {
namespace integrationSchemes
{
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class Euler Declaration Class Euler Declaration
@ -56,11 +62,8 @@ public:
// Constructors // Constructors
//- Construct from components //- Construct
Euler(const word& phiName, const dictionary& dict); Euler();
//- Copy constructor
Euler(const Euler& is);
//- Construct and return clone //- Construct and return clone
virtual autoPtr<IntegrationScheme<Type>> clone() const virtual autoPtr<IntegrationScheme<Type>> clone() const
@ -75,19 +78,14 @@ public:
// Member Functions // Member Functions
//- Perform the integration //- Return the integration factor
virtual typename IntegrationScheme<Type>::integrationResult integrate virtual scalar factor(const scalar dt, const scalar Beta) const;
(
const Type& phi,
const scalar dt,
const Type& alphaBeta,
const scalar beta
) const;
}; };
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace integrationSchemes
} // End namespace Foam } // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

48
src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -28,22 +28,7 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type> template<class Type>
Foam::IntegrationScheme<Type>::IntegrationScheme Foam::IntegrationScheme<Type>::IntegrationScheme()
(
const word& phiName,
const dictionary& dict
)
:
phiName_(phiName),
dict_(dict)
{}
template<class Type>
Foam::IntegrationScheme<Type>::IntegrationScheme(const IntegrationScheme& is)
:
phiName_(is.phiName_),
dict_(is.dict_)
{} {}
@ -54,6 +39,35 @@ Foam::IntegrationScheme<Type>::~IntegrationScheme()
{} {}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Type Foam::IntegrationScheme<Type>::delta
(
const Type& phi,
const scalar dt,
const Type& Alpha,
const scalar Beta
) const
{
return delta(phi, dt, Beta, Alpha, Beta);
}
template<class Type>
Type Foam::IntegrationScheme<Type>::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" #include "IntegrationSchemeNew.C"

148
src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -25,7 +25,35 @@ Class
Foam::IntegrationScheme Foam::IntegrationScheme
Description 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 SourceFiles
IntegrationScheme.C IntegrationScheme.C
@ -51,79 +79,6 @@ namespace Foam
template<class Type> template<class Type>
class IntegrationScheme 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: public:
//- Runtime type information //- Runtime type information
@ -136,22 +91,16 @@ public:
( (
autoPtr, autoPtr,
IntegrationScheme, IntegrationScheme,
dictionary, word,
( (),
const word& phiName, ()
const dictionary& dict
),
(phiName, dict)
); );
// Constructors // Constructors
//- Construct from components //- Construct
IntegrationScheme(const word& phiName, const dictionary& dict); IntegrationScheme();
//- Copy constructor
IntegrationScheme(const IntegrationScheme& is);
//- Construct and return clone //- Construct and return clone
virtual autoPtr<IntegrationScheme<Type>> clone() const = 0; virtual autoPtr<IntegrationScheme<Type>> clone() const = 0;
@ -159,7 +108,7 @@ public:
// Selectors // Selectors
//- Return a reference to the selected radiation model //- Select an integration scheme
static autoPtr<IntegrationScheme> New static autoPtr<IntegrationScheme> New
( (
const word& phiName, const word& phiName,
@ -173,14 +122,27 @@ public:
// Member Functions // Member Functions
//- Perform the Integration //- Perform the full integration
virtual integrationResult integrate Type delta
( (
const Type& phi, const Type& phi,
const scalar dt, 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 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 \ defineTemplateRunTimeSelectionTable \
( \ ( \
IntegrationScheme<Type>, \ IntegrationScheme<Type>, \
dictionary \ word \
); );
@ -205,7 +167,7 @@ public:
\ \
defineNamedTemplateTypeNameAndDebug(SS<Type>, 0); \ defineNamedTemplateTypeNameAndDebug(SS<Type>, 0); \
\ \
IntegrationScheme<Type>::adddictionaryConstructorToTable<SS<Type>> \ IntegrationScheme<Type>::addwordConstructorToTable<SS<Type>> \
add##SS##Type##ConstructorToTable_; add##SS##Type##ConstructorToTable_;

12
src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationSchemeNew.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -41,20 +41,20 @@ Foam::IntegrationScheme<Type>::New
Info<< "Selecting " << phiName << " integration scheme " Info<< "Selecting " << phiName << " integration scheme "
<< schemeName << endl; << schemeName << endl;
typename dictionaryConstructorTable::iterator cstrIter = typename wordConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(schemeName); wordConstructorTablePtr_->find(schemeName);
if (cstrIter == dictionaryConstructorTablePtr_->end()) if (cstrIter == wordConstructorTablePtr_->end())
{ {
FatalErrorInFunction FatalErrorInFunction
<< "Unknown integration scheme type " << "Unknown integration scheme type "
<< schemeName << nl << nl << schemeName << nl << nl
<< "Valid integration scheme types are:" << nl << "Valid integration scheme types are:" << nl
<< dictionaryConstructorTablePtr_->sortedToc() << nl << wordConstructorTablePtr_->sortedToc() << nl
<< exit(FatalError); << exit(FatalError);
} }
return autoPtr<IntegrationScheme<Type>>(cstrIter()(phiName, dict)); return autoPtr<IntegrationScheme<Type>>(cstrIter()());
} }
// ************************************************************************* // // ************************************************************************* //

4
src/lagrangian/intermediate/IntegrationScheme/makeIntegrationSchemes.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -32,6 +32,8 @@ License
namespace Foam namespace Foam
{ {
using namespace integrationSchemes;
makeIntegrationScheme(scalar); makeIntegrationScheme(scalar);
makeIntegrationSchemeType(Euler, scalar); makeIntegrationSchemeType(Euler, scalar);
makeIntegrationSchemeType(Analytical, scalar); makeIntegrationSchemeType(Analytical, scalar);

21
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
View File

@ -177,26 +177,27 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
// Momentum source due to particle forces // Momentum source due to particle forces
const forceSuSp Fcp = forces.calcCoupled(p, ttd, dt, mass, Re, mu); const forceSuSp Fcp = forces.calcCoupled(p, ttd, dt, mass, Re, mu);
const forceSuSp Fncp = forces.calcNonCoupled(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); const scalar massEff = forces.massEff(p, ttd, mass);
// New particle velocity // New particle velocity
//~~~~~~~~~~~~~~~~~~~~~~ //~~~~~~~~~~~~~~~~~~~~~~
// Update velocity - treat as 3-D const vector acp = (Fcp.Sp()*td.Uc() + Fcp.Su())/massEff;
const vector abp = (Feff.Sp()*td.Uc() + (Feff.Su() + Su))/massEff; const vector ancp = (Fncp.Sp()*td.Uc() + Fncp.Su() + Su)/massEff;
const scalar bp = Feff.Sp()/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<vector>::integrationResult Ures = vector Unew = U_ + deltaUcp + deltaUncp;
cloud.UIntegrator().integrate(U_, dt, abp, bp);
vector Unew = Ures.value(); dUTrans -= massEff*deltaUcp;
// note: Feff.Sp() and Fc.Sp() must be the same Spu = dt*Fcp.Sp();
dUTrans += dt*(Feff.Sp()*(Ures.average() - td.Uc()) - Fcp.Su());
// Apply correction to velocity and dUTrans for reduced-D cases // Apply correction to velocity and dUTrans for reduced-D cases
const polyMesh& mesh = cloud.pMesh(); const polyMesh& mesh = cloud.pMesh();

54
src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
View File

@ -263,57 +263,41 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
const scalar d = this->d(); const scalar d = this->d();
const scalar rho = this->rho(); 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 // Calc heat transfer coefficient
scalar htc = cloud.heatTransfer().htc(d, Re, Pr, kappa, NCpW); scalar htc = cloud.heatTransfer().htc(d, Re, Pr, kappa, NCpW);
if (mag(htc) < ROOTVSMALL && !cloud.radiation()) // Calculate the integration coefficients
{ const scalar bcp = htc*As/(m*Cp_);
return const scalar acp = bcp*td.Tc();
max scalar ancp = Sh;
(
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);
if (cloud.radiation()) if (cloud.radiation())
{ {
tetIndices tetIs = this->currentTetIndices(); const tetIndices tetIs = this->currentTetIndices();
const scalar Gc = td.GInterp().interpolate(this->coordinates(), tetIs); const scalar Gc = td.GInterp().interpolate(this->coordinates(), tetIs);
const scalar sigma = physicoChemical::sigma.value(); const scalar sigma = physicoChemical::sigma.value();
const scalar epsilon = cloud.constProps().epsilon0(); const scalar epsilon = cloud.constProps().epsilon0();
// Assume constant source ancp += As*epsilon*(Gc/4.0 - sigma*pow4(T_));
scalar s = epsilon*(Gc/4.0 - sigma*pow4(T_));
ap += s/htc;
} }
ancp /= m*Cp_;
// Integrate to find the new parcel temperature // Integrate to find the new parcel temperature
IntegrationScheme<scalar>::integrationResult Tres = const scalar deltaTcp =
cloud.TIntegrator().integrate(T_, dt, ap*bp, bp); cloud.TIntegrator().delta(T_, dt, bcp, acp, bcp);
const scalar deltaTncp =
cloud.TIntegrator().delta(T_, dt, bcp, ancp, 0);
scalar Tnew = // Calculate the new temperature and the transfer terms
min scalar Tnew = T_ + deltaTcp + deltaTncp;
( Tnew = min(max(Tnew, cloud.constProps().TMin()), cloud.constProps().TMax());
max
(
Tres.value(),
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; return Tnew;
} }