ENH: Updates to the lagrangian intermediate lib coupling terms

This commit is contained in:
andy
2011-02-28 16:25:11 +00:00
parent 4ae3d1c0e5
commit c4a4055d3a
9 changed files with 32 additions and 35 deletions

View File

@ -127,7 +127,7 @@ void Foam::KinematicParcel<ParcelType>::calc
// ~~~~~~
// Calculate new particle velocity
scalar Cud = 0.0;
scalar Spu = 0.0;
vector U1 =
calcVelocity
(
@ -142,7 +142,7 @@ void Foam::KinematicParcel<ParcelType>::calc
mass0,
Su,
dUTrans,
Cud
Spu
);
@ -154,7 +154,7 @@ void Foam::KinematicParcel<ParcelType>::calc
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Update momentum transfer coefficient
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
td.cloud().UCoeff()[cellI] += np0*Spu;
}
@ -179,7 +179,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
const scalar mass,
const vector& Su,
vector& dUTrans,
scalar& Cud
scalar& Spu
) const
{
typedef typename TrackData::cloudType cloudType;
@ -202,13 +202,14 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
const vector ap = Uc_ + (Feff.Su() + Su)/(Feff.Sp() + ROOTVSMALL);
const scalar bp = Feff.Sp()/mass;
Cud = bp;
Spu = Feff.Sp();
IntegrationScheme<vector>::integrationResult Ures =
td.cloud().UIntegrator().integrate(U, dt, ap, bp);
vector Unew = Ures.value();
// note: Feff.Sp() and Fc.Sp() must be the same
dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su());
// Apply correction to velocity and dUTrans for reduced-D cases

View File

@ -303,7 +303,7 @@ protected:
const scalar mass, // mass
const vector& Su, // explicit particle momentum source
vector& dUTrans, // momentum transfer to carrier
scalar& Cud // linearised drag coefficient
scalar& Spu // linearised drag coefficient
) const;

View File

@ -353,7 +353,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
// ~~~~~~
// Calculate new particle velocity
scalar Cud = 0;
scalar Spu = 0;
vector U1 =
this->calcVelocity
(
@ -368,11 +368,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
mass0,
Su,
dUTrans,
Cud
Spu
);
dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
// Accumulate carrier phase source terms
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -418,7 +416,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Update momentum transfer coefficient
td.cloud().UCoeff()[cellI] += np0*0.5*(mass0 + mass1)*Cud;
td.cloud().UCoeff()[cellI] += np0*Spu;
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*dhsTrans;

View File

@ -362,7 +362,7 @@ void Foam::ReactingParcel<ParcelType>::calc
// ~~~~~~
// Calculate new particle velocity
scalar Cud = 0.0;
scalar Spu = 0.0;
vector U1 =
this->calcVelocity
(
@ -377,7 +377,7 @@ void Foam::ReactingParcel<ParcelType>::calc
mass0,
Su,
dUTrans,
Cud
Spu
);
dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
@ -399,7 +399,7 @@ void Foam::ReactingParcel<ParcelType>::calc
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Update momentum transfer coefficient
td.cloud().UCoeff()[cellI] += np0*0.5*(mass0 + mass1)*Cud;
td.cloud().UCoeff()[cellI] += np0*Spu;
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*dhsTrans;

View File

@ -234,7 +234,7 @@ void Foam::ThermoParcel<ParcelType>::calc
// ~~~~~~
// Calculate new particle velocity
scalar Cud = 0.0;
scalar Spu = 0.0;
vector U1 =
this->calcVelocity
(
@ -249,7 +249,7 @@ void Foam::ThermoParcel<ParcelType>::calc
mass0,
Su,
dUTrans,
Cud
Spu
);
@ -261,7 +261,7 @@ void Foam::ThermoParcel<ParcelType>::calc
td.cloud().UTrans()[cellI] += np0*dUTrans;
// Update momentum transfer coefficient
td.cloud().UCoeff()[cellI] += np0*mass0*Cud;
td.cloud().UCoeff()[cellI] += np0*Spu;
// Update sensible enthalpy transfer
td.cloud().hsTrans()[cellI] += np0*dhsTrans;

View File

@ -28,9 +28,9 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class CloudType>
Foam::scalar Foam::NonSphereDragForce<CloudType>::Cd(const scalar Re) const
Foam::scalar Foam::NonSphereDragForce<CloudType>::CdRe(const scalar Re) const
{
return 24.0/Re*(1.0 + a_*pow(Re, b_)) + Re*c_/(Re + d_);
return 24.0*(1.0 + a_*pow(Re, b_)) + Re*c_/(1 + d_/(Re + ROOTVSMALL));
}
@ -93,7 +93,7 @@ Foam::NonSphereDragForce<CloudType>::~NonSphereDragForce()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CloudType>
Foam::forceSuSp Foam::NonSphereDragForce<CloudType>::calcNonCoupled
Foam::forceSuSp Foam::NonSphereDragForce<CloudType>::calcCoupled
(
const typename CloudType::parcelType& p,
const scalar dt,
@ -104,8 +104,7 @@ Foam::forceSuSp Foam::NonSphereDragForce<CloudType>::calcNonCoupled
{
forceSuSp value(vector::zero, 0.0);
const scalar ReCorr = max(Re, 1e-6);
value.Sp() = mass*0.75*muc*Cd(ReCorr)*ReCorr/(p.rho()*sqr(p.d()));
value.Sp() = mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d()));
return value;
}

View File

@ -97,8 +97,8 @@ protected:
// Private Member Functions
//- Drag coefficient
scalar Cd(const scalar Re) const;
//- Drag coefficient multiplied by Reynolds number
scalar CdRe(const scalar Re) const;
public:
@ -140,7 +140,7 @@ public:
// Evaluation
//- Calculate the non-coupled force
virtual forceSuSp calcNonCoupled
virtual forceSuSp calcCoupled
(
const typename CloudType::parcelType& p,
const scalar dt,

View File

@ -28,15 +28,15 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class CloudType>
Foam::scalar Foam::SphereDragForce<CloudType>::Cd(const scalar Re) const
Foam::scalar Foam::SphereDragForce<CloudType>::CdRe(const scalar Re) const
{
if (Re > 1000.0)
{
return 0.424;
return 0.424*Re;
}
else
{
return 24.0/Re*(1.0 + 1.0/6.0*pow(Re, 2.0/3.0));
return 24.0*(1.0 + 1.0/6.0*pow(Re, 2.0/3.0));
}
}
@ -76,7 +76,7 @@ Foam::SphereDragForce<CloudType>::~SphereDragForce()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class CloudType>
Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcNonCoupled
Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcCoupled
(
const typename CloudType::parcelType& p,
const scalar dt,
@ -87,8 +87,7 @@ Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcNonCoupled
{
forceSuSp value(vector::zero, 0.0);
const scalar ReCorr = max(Re, 1e-6);
value.Sp() = mass*0.75*muc*Cd(ReCorr)*ReCorr/(p.rho()*sqr(p.d()));
value.Sp() = mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d()));
return value;
}

View File

@ -49,8 +49,8 @@ class SphereDragForce
{
// Private Member Functions
//- Drag coefficient
scalar Cd(const scalar Re) const;
//- Drag coefficient multiplied by Reynolds number
scalar CdRe(const scalar Re) const;
public:
@ -91,8 +91,8 @@ public:
// Evaluation
//- Calculate the non-coupled force
virtual forceSuSp calcNonCoupled
//- Calculate the coupled force
virtual forceSuSp calcCoupled
(
const typename CloudType::parcelType& p,
const scalar dt,