From 55da1f976a3ec5168081de206a01353163fdfb5f Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Fri, 13 Mar 2009 19:17:25 +0000
Subject: [PATCH] updates for parcel calc routined
---
.../Templates/KinematicCloud/KinematicCloud.C | 4 ++
.../Templates/KinematicCloud/KinematicCloud.H | 27 ++++++++++++
.../KinematicCloud/KinematicCloudI.H | 22 ++++++++++
.../KinematicParcel/KinematicParcel.C | 36 +++++++++++++---
.../KinematicParcel/KinematicParcel.H | 8 ++++
.../KinematicParcel/KinematicParcelI.H | 11 ++++-
.../Templates/ReactingParcel/ReactingParcel.C | 43 ++++++++-----------
.../Templates/ThermoParcel/ThermoParcel.C | 39 +++++++++--------
.../Templates/ThermoParcel/ThermoParcel.H | 1 +
.../Kinematic/DragModel/DragModel/DragModel.C | 2 +-
10 files changed, 144 insertions(+), 49 deletions(-)
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index 6584495b637..6df424c523a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -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
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index 8a03a73cd7e..77a61b02d29 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -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
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index f129df9bdb0..adae027b878 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -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
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 7fc3a62b792..c9df5a57768 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -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;
}
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
index c60a2164984..b6df66d414c 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H
@@ -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;
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
index 3e8d49b7ad4..4a6be04a971 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H
@@ -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>
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 8e8b9f9d68a..65d8bc1d8fd 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -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;
-
- // New specific heat capacity of mixture - using old temperature
- scalar cp1 = td.cloud().composition().cp(0, Y_, pc_, T0);
+ scalar dhTrans = ShHT + ShPC;
- // 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
)
{
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
index 61cbdf9535d..da7a31200ee 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C
@@ -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();
}
diff --git a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
index bc66a2a7d2c..8506bb7c0ad 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H
@@ -210,6 +210,7 @@ protected:
TrackData& td,
const scalar dt,
const label cellI,
+ const scalar Sh,
scalar& htc,
scalar& dhTrans
);
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/DragModel/DragModel/DragModel.C b/src/lagrangian/intermediate/submodels/Kinematic/DragModel/DragModel/DragModel.C
index 273e93b0e17..f1b98e609c2 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/DragModel/DragModel/DragModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/DragModel/DragModel/DragModel.C
@@ -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);
--
GitLab