zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

updates for parcel calc routined

Browse Source
This commit is contained in:
andy
2009-03-13 19:17:25 +00:00
parent 5e76511069
commit 55da1f976a
10 changed files with 144 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
View File

@ -102,6 +102,10 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
mu_(mu),
g_(g),
interpolationSchemes_(particleProperties_.subDict("interpolationSchemes")),
forcesDict_(particleProperties_.subDict("forces")),
forceGravity_(forcesDict_.lookup("gravity")),
forceVirtualMass_(forcesDict_.lookup("virtualMass")),
forcePressureGradient_(forcesDict_.lookup("pressureGradient")),
dispersionModel_
(
DispersionModel<KinematicCloud<ParcelType> >::New

27
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
View File

@ -133,6 +133,21 @@ class KinematicCloud
dictionary interpolationSchemes_;
// Forces to include in particle motion evaluation
//- Dictionary of forces
dictionary forcesDict_;
//- Gravity
Switch forceGravity_;
//- Virtual mass
Switch forceVirtualMass_;
//- Pressure gradient
Switch forcePressureGradient_;
// References to the cloud sub-models
//- Dispersion model
@ -252,6 +267,18 @@ public:
inline const dictionary& interpolationSchemes() const;
// Forces to include in particle motion evaluation
//- Return reference to the gravity force flag
inline Switch forceGravity() const;
//- Return reference to the virtual mass force flag
inline Switch forceVirtualMass() const;
//- Return reference to the pressure gradient force flag
inline Switch forcePressureGradient() const;
// Sub-models
//- Return reference to dispersion model

22
src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
View File

@ -108,6 +108,28 @@ Foam::KinematicCloud<ParcelType>::interpolationSchemes() const
}
template<class ParcelType>
inline Foam::Switch Foam::KinematicCloud<ParcelType>::forceGravity() const
{
return forceGravity_;
}
template<class ParcelType>
inline Foam::Switch Foam::KinematicCloud<ParcelType>::forceVirtualMass() const
{
return forceVirtualMass_;
}
template<class ParcelType>
inline Foam::Switch Foam::KinematicCloud<ParcelType>::
forcePressureGradient() const
{
return forcePressureGradient_;
}
template<class ParcelType>
inline const Foam::DispersionModel<Foam::KinematicCloud<ParcelType> >&
Foam::KinematicCloud<ParcelType>::dispersion() const

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

@ -27,6 +27,8 @@ License
#include "KinematicParcel.H"
#include "dimensionedConstants.H"
#include "fvcGrad.H"
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
template<class ParcelType>
@ -69,7 +71,7 @@ void Foam::KinematicParcel<ParcelType>::calc
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
vector dUTrans = vector::zero;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
const vector U1 = calcVelocity(td, dt, vector::zero, mass(), Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -88,7 +90,7 @@ void Foam::KinematicParcel<ParcelType>::calc
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 3. Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
this->U() = U1;
U_ = U1;
}
@ -98,27 +100,51 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
TrackData& td,
const scalar dt,
const vector& Fx,
const scalar mass,
scalar& Cud,
vector& dUTrans
vector& dUTrans,
) const
{
// Return linearised term from drag model
Cud = td.cloud().drag().Cu(U_ - Uc_, d_, rhoc_, rho_, muc_);
// Initialise total force (per unit mass)
vector Ftot = vector::zero;
// Gravity force
if (td.cloud().forceGravity())
{
Ftot += td.g()*(1 - rhoc_/rho_);
}
// Virtual mass force
if (td.cloud().forceVirtualMass())
{
// Ftot += td.constProps().Cvm()*rhoc_/rho_*d(Uc - U_)/dt;
}
// Pressure gradient force
if (td.cloud().forcePressureGradient())
{
Ftot += rhoc_/rho_*(U_ & fvc::grad(Uc_))
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Set new particle velocity
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Update velocity - treat as 3-D
const vector ap = Uc_ + (1 - rhoc_/rho_)/(Cud + VSMALL)*td.g();
const vector ap = Uc_ + (Ftot + Fx)/(Cud + VSMALL);
const scalar bp = Cud;
vector Unew = td.cloud().UIntegrator().integrate(U_, dt, ap, bp).value();
// Calculate the momentum transfer to the continuous phase
// - do not include gravity impulse
dUTrans = -mass()*(Unew - U_ - dt*td.g());
// TODO: USE AVERAGE PARTICLE MASS
dUTrans = -mass*(Unew - U_ - dt*td.g());
return Unew;
}

8
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
View File

@ -93,6 +93,9 @@ public:
//- Minimum particle mass [kg]
const scalar minParticleMass_;
//- Virtual mass coefficient []
const scalar Cvm_;
public:
@ -106,6 +109,9 @@ public:
//- Return const access to the minimum particle mass
inline scalar minParticleMass() const;
//- Return const access to the virtual mass coefficient
inline scalar Cvm() const;
};
@ -228,6 +234,8 @@ protected:
(
TrackData& td,
const scalar dt,
const vector& Fx,
const scalar mass,
scalar& Cud,
vector& dUTrans
) const;

11
src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
View File

@ -33,7 +33,8 @@ inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
)
:
rho0_(dimensionedScalar(dict.lookup("rho0")).value()),
minParticleMass_(dimensionedScalar(dict.lookup("minParticleMass")).value())
minParticleMass_(dimensionedScalar(dict.lookup("minParticleMass")).value()),
Cvm_(dimensionedScalar(dict.lookup("Cvm")).value())
{}
@ -103,6 +104,14 @@ Foam::KinematicParcel<ParcelType>::constantProperties::minParticleMass() const
}
template <class ParcelType>
inline Foam::scalar
Foam::KinematicParcel<ParcelType>::constantProperties::Cvm() const
{
return Cvm_;
}
// * * * * * * * * * * * trackData Member Functions * * * * * * * * * * * * //
template <class ParcelType>

43
src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
View File

@ -82,35 +82,33 @@ void Foam::ReactingParcel<ParcelType>::calc
scalar Cud = 0.0;
vector dUTrans = vector::zero;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
// 2. Calculate heat transfer
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
scalar dhHT = 0.0;
// TODO: T1 no longer used - return dhHT instead??????????????????????????????
// scalar T1 = calcHeatTransfer(td, dt, cellI, htc, dhHT);
calcHeatTransfer(td, dt, cellI, htc, dhHT);
const vector U1 = calcVelocity(td, dt, vector::zero, mass0, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// 3. Calculate phase change
// 2. Calculate phase change
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Mass transfer from particle to carrier phase
scalarList dMassPC(td.cloud().gases().size(), 0.0);
const scalar dMassPCTot = calcPhaseChange(td, dt, cellI, T0, 1.0, dMassPC);
// Update particle component mass fractions
updateMassFraction(mass0, dMassPC, Y_);
// Enthalpy change due to change in particle composition (sink)
scalar dhPC = -dMassPCTot*td.cloud().composition().L(0, Y_, pc, T0);
scalar ShPC = -dMassPCTot*td.cloud().composition().L(0, Y_, pc, T0);
// Enthalpy change due to species released into the carrier (source)
scalar HEff = td.cloud().composition().H(0, Y_, pc, T0);
dhPC += dMassPCTot*HEff;
ShPC += dMassPCTot*HEff;
// Update particle component mass fractions
updateMassFraction(mass0, dMassPC, Y_);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
// 3. Calculate heat transfer
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
scalar ShHT = 0.0;
scalar T1 = calcHeatTransfer(td, dt, cellI, ShPC, htc, ShHT);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -121,13 +119,10 @@ void Foam::ReactingParcel<ParcelType>::calc
scalar mass1 = mass0 - dMassPCTot;
// Total enthalpy transfer from the particle to the carrier phase
scalar dhTrans = dhHT + dhPC;
scalar dhTrans = ShHT + ShPC;
// New specific heat capacity of mixture - using old temperature
scalar cp1 = td.cloud().composition().cp(0, Y_, pc_, T0);
// New particle temperature - using average mass over the time interval
scalar T1 = T0 + dhTrans/(0.5*(mass0 + mass1)*cp1);
// New specific heat capacity of mixture - using new temperature
scalar cp1 = td.cloud().composition().cp(0, Y_, pc_, T1);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -203,7 +198,7 @@ Foam::scalar Foam::ReactingParcel<ParcelType>::calcPhaseChange
const scalar dt,
const label cellI,
const scalar T,
const scalar YPhase,
const scalar YPhase, // TODO: NEEDED?????????????????????????????????????????
scalarList& dMassMT
)
{

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

@ -72,18 +72,18 @@ void Foam::ThermoParcel<ParcelType>::calc
scalar dhTrans = 0.0;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 2. Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 3. Calculate heat transfer - update T
// 2. Calculate heat transfer - update T
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar htc = 0.0;
const scalar T1 = calcHeatTransfer(td, dt, cellI, htc, dhTrans);
const scalar T1 = calcHeatTransfer(td, dt, cellI, 0.0, htc, dhTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 3. Calculate velocity - update U
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalar Cud = 0.0;
const vector U1 = calcVelocity(td, dt, vector::zero, Cud, mass0, dUTrans);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -108,7 +108,7 @@ void Foam::ThermoParcel<ParcelType>::calc
// 5. Set new particle properties
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
this->U() = U1;
this->T() = T1;
T_ = T1;
}
@ -119,6 +119,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
TrackData& td,
const scalar dt,
const label cellI,
const scalar Sh,
scalar& htc,
scalar& dhTrans
)
@ -142,9 +143,12 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
this->muc_
);
//- Assuming diameter = diameter at start of time step
scalar Ap = this->areasS();
// Determine ap and bp coefficients
scalar ap = Tc_;
scalar bp = htc;
scalar ap = Tc_ + Sh/(htc*Ap + ROOTVSMALL);
scalar bp = 6.0*htc/(this->rho_*this->d_*cp_);
if (td.cloud().radiation())
{
// Carrier phase incident radiation field
@ -156,19 +160,18 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
// Helper variables
const scalar sigma = radiation::sigmaSB.value();
const scalar epsilon = td.constProps().epsilon0();
const scalar epsilonSigmaT3 = epsilon*sigma*pow3(T_);
ap = (htc*Tc_ + 0.25*epsilon*G[cellI])/(htc + epsilonSigmaT3);
bp += epsilonSigmaT3;
const scalar D = epsilon*sigma*pow3(T_)/(htc + ROOTVSMALL) + 1.0;
ap += 0.25*epsilon*G[cellI]/(htc + ROOTVSMALL);
ap /= D;
bp *= D;
}
bp *= 6.0/(this->rho_*this->d_*cp_);
// Integrate to find the new parcel temperature
IntegrationScheme<scalar>::integrationResult Tres =
td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
// Using average parcel temperature for enthalpy transfer calculation
dhTrans = dt*this->areaS()*htc*(Tres.average() - Tc_);
dhTrans = dt*Ap*htc*(Tres.average() - Tc_);
return Tres.value();
}

1
src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
View File

@ -210,6 +210,7 @@ protected:
TrackData& td,
const scalar dt,
const label cellI,
const scalar Sh,
scalar& htc,
scalar& dhTrans
);

2
src/lagrangian/intermediate/submodels/Kinematic/DragModel/DragModel/DragModel.C
View File

@ -75,7 +75,7 @@ Foam::scalar Foam::DragModel<CloudType>::Cu
{
const scalar magUr = mag(Ur);
const scalar Re = rhoc*magUr*d/(mu + SMALL);
const scalar Re = rhoc*magUr*d/(mu + ROOTVSMALL);
const scalar Cd = this->Cd(Re);