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
\\ / 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<class Type>
Foam::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)
Foam::integrationSchemes::Analytical<Type>::Analytical()
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
Foam::Analytical<Type>::~Analytical()
Foam::integrationSchemes::Analytical<Type>::~Analytical()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
typename Foam::IntegrationScheme<Type>::integrationResult
Foam::Analytical<Type>::integrate
Foam::scalar Foam::integrationSchemes::Analytical<Type>::factor
(
const Type& phi,
const scalar dt,
const Type& alphaBeta,
const scalar beta
const scalar Beta
) const
{
typename IntegrationScheme<Type>::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;
}

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

@ -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<IntegrationScheme<Type>> clone() const
@ -78,19 +81,14 @@ public:
// Member Functions
//- Perform the integration
virtual typename IntegrationScheme<Type>::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
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@ -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<class Type>
Foam::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)
Foam::integrationSchemes::Euler<Type>::Euler()
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
Foam::Euler<Type>::~Euler()
Foam::integrationSchemes::Euler<Type>::~Euler()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
typename Foam::IntegrationScheme<Type>::integrationResult
Foam::Euler<Type>::integrate
Foam::scalar Foam::integrationSchemes::Euler<Type>::factor
(
const Type& phi,
const scalar dt,
const Type& alphaBeta,
const scalar beta
const scalar Beta
) const
{
typename IntegrationScheme<Type>::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);
}

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

@ -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<IntegrationScheme<Type>> clone() const
@ -75,19 +78,14 @@ public:
// Member Functions
//- Perform the integration
virtual typename IntegrationScheme<Type>::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
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@ -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<class Type>
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_)
Foam::IntegrationScheme<Type>::IntegrationScheme()
{}
@ -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"

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

@ -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 Type>
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<IntegrationScheme<Type>> clone() const = 0;
@ -159,7 +108,7 @@ public:
// Selectors
//- Return a reference to the selected radiation model
//- Select an integration scheme
static autoPtr<IntegrationScheme> 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<Type>, \
dictionary \
word \
);
@ -205,7 +167,7 @@ public:
\
defineNamedTemplateTypeNameAndDebug(SS<Type>, 0); \
\
IntegrationScheme<Type>::adddictionaryConstructorToTable<SS<Type>> \
IntegrationScheme<Type>::addwordConstructorToTable<SS<Type>> \
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
\\ / 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<Type>::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<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
\\ / 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);

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
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<vector>::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();

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 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<scalar>::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;
}