28 #include "BlendedInterfacialModel.H" 29 #include "dragModel.H" 30 #include "virtualMassModel.H" 31 #include "liftModel.H" 32 #include "wallLubricationModel.H" 33 #include "turbulentDispersionModel.H" 50 template<
class BasePhaseSystem>
54 const phasePairKey& key
57 if (dragModels_.found(key))
59 return dragModels_[key]->K();
63 return tmp<volScalarField>
69 dragModel::typeName +
":K",
70 this->mesh_.time().timeName(),
84 template<
class BasePhaseSystem>
88 const phasePairKey& key
91 if (dragModels_.found(key))
93 return dragModels_[key]->Kf();
97 return tmp<surfaceScalarField>
103 dragModel::typeName +
":K",
104 this->mesh_.time().timeName(),
118 template<
class BasePhaseSystem>
122 const phasePairKey& key
125 if (virtualMassModels_.found(key))
127 return virtualMassModels_[key]->K();
131 return tmp<volScalarField>
137 virtualMassModel::typeName +
":K",
138 this->mesh_.time().timeName(),
152 template<
class BasePhaseSystem>
158 phaseSystem::phasePairTable,
163 const phasePair& pair(phasePairIter());
176 if (!pair.phase1().stationary())
181 eqn += dmdt21*pair.phase2().U() -
fvm::Sp(dmdt21, eqn.psi());
184 if (!pair.phase2().stationary())
189 eqn -= dmdt12*pair.phase1().U() -
fvm::Sp(dmdt12, eqn.psi());
197 template<
class BasePhaseSystem>
204 BasePhaseSystem(mesh)
206 this->generatePairsAndSubModels
212 this->generatePairsAndSubModels
218 this->generatePairsAndSubModels
224 this->generatePairsAndSubModels
227 wallLubricationModels_
230 this->generatePairsAndSubModels
232 "turbulentDispersion",
233 turbulentDispersionModels_
243 const phasePair& pair(this->phasePairs_[dragModelIter.key()]);
261 dragModelIter()->Kf()
268 virtualMassModelTable,
273 const phasePair& pair(this->phasePairs_[virtualMassModelIter.key()]);
281 virtualMassModelIter()->K()
291 virtualMassModelIter()->Kf()
300 template<
class BasePhaseSystem>
308 template<
class BasePhaseSystem>
313 autoPtr<phaseSystem::momentumTransferTable> eqnsPtr
320 forAll(this->movingPhases(), movingPhasei)
322 const phaseModel& phase = this->movingPhases()[movingPhasei];
339 *Kds_[dragModelIter.key()] = dragModelIter()->K();
340 *Kdfs_[dragModelIter.key()] = dragModelIter()->Kf();
347 const phasePair& pair(this->phasePairs_[KdIter.key()]);
351 if (!iter().stationary())
363 virtualMassModelTable,
368 *Vms_[virtualMassModelIter.key()] = virtualMassModelIter()->K();
369 *Vmfs_[virtualMassModelIter.key()] = virtualMassModelIter()->Kf();
376 const phasePair& pair(this->phasePairs_[VmIter.key()]);
380 const phaseModel& phase = iter();
381 const phaseModel& otherPhase = iter.otherPhase();
383 if (!phase.stationary())
398 + this->MRF_.DDt(Vm, U - otherPhase.U());
404 addMassTransferMomentumTransfer(eqns);
410 template<
class BasePhaseSystem>
415 autoPtr<phaseSystem::momentumTransferTable> eqnsPtr
422 forAll(this->movingPhases(), movingPhasei)
424 const phaseModel& phase = this->movingPhases()[movingPhasei];
434 PtrList<fvVectorMatrix> UgradUs(this->phaseModels_.size());
437 const phaseModel& phase = this->phaseModels_[
phasei];
439 if (!phase.stationary())
460 const phasePair& pair(this->phasePairs_[VmIter.key()]);
464 const phaseModel& phase = iter();
465 const phaseModel& otherPhase = iter.otherPhase();
467 if (!phase.stationary())
469 *eqns[phase.name()] -=
472 UgradUs[phase.index()]
473 - (UgradUs[otherPhase.index()] & otherPhase.U())
480 addMassTransferMomentumTransfer(eqns);
486 template<
class BasePhaseSystem>
490 PtrList<surfaceScalarField>
AFfs(this->phaseModels_.size());
496 const phasePair& pair(this->phasePairs_[KdfIter.key()]);
500 this->addField(iter(),
"AFf", Kf,
AFfs);
508 const phasePair& pair(this->phasePairs_[VmfIter.key()]);
516 if (this->fillFields_)
525 template<
class BasePhaseSystem>
529 const PtrList<volScalarField>& rAUs
532 PtrList<surfaceScalarField>
phiFs(this->phaseModels_.size());
543 const phasePair& pair(this->phasePairs_[liftModelIter.key()]);
564 wallLubricationModelTable,
565 wallLubricationModels_,
566 wallLubricationModelIter
571 pair(this->phasePairs_[wallLubricationModelIter.key()]);
593 const phaseModel& phase = this->phaseModels_[
phasei];
605 this->addField(phase,
"phiF", pPrimeByAf*snGradAlpha1,
phiFs);
607 const bool implicitPhasePressure =
608 this->mesh_.solverDict(phase.volScalarField::name()).
609 template lookupOrDefault<Switch>
611 "implicitPhasePressure",
615 if (implicitPhasePressure)
617 this->addField(phase,
"DByAf", pPrimeByAf, DByAfs_);
624 turbulentDispersionModelTable,
625 turbulentDispersionModels_,
626 turbulentDispersionModelIter
630 pair(this->phasePairs_[turbulentDispersionModelIter.key()]);
651 if (DByAfs_.found(pair.phase1().name()))
653 this->addField(pair.phase1(),
"DByAf", DByA1f, DByAfs_);
657 if (this->fillFields_)
666 template<
class BasePhaseSystem>
670 const PtrList<surfaceScalarField>& rAUfs
673 PtrList<surfaceScalarField>
phiFfs(this->phaseModels_.size());
679 const phasePair& pair(this->phasePairs_[VmfIter.key()]);
687 - rAUfs[iter().index()]*
Vmf 690 + iter.otherPhase().DUDtf()
706 const phasePair& pair(this->phasePairs_[liftModelIter.key()]);
712 rAUfs[pair.phase1().index()]*
Ff,
719 - rAUfs[pair.phase2().index()]*
Ff,
727 wallLubricationModelTable,
728 wallLubricationModels_,
729 wallLubricationModelIter
734 pair(this->phasePairs_[wallLubricationModelIter.key()]);
740 rAUfs[pair.phase1().index()]*
Ff,
747 - rAUfs[pair.phase2().index()]*
Ff,
756 const phaseModel& phase = this->phaseModels_[
phasei];
768 this->addField(phase,
"phiFf", pPrimeByAf*snGradAlpha1,
phiFfs);
770 const bool implicitPhasePressure =
771 this->mesh_.solverDict(phase.volScalarField::name()).
772 template lookupOrDefault<Switch>
774 "implicitPhasePressure",
778 if (implicitPhasePressure)
780 this->addField(phase,
"DByAf", pPrimeByAf, DByAfs_);
787 turbulentDispersionModelTable,
788 turbulentDispersionModels_,
789 turbulentDispersionModelIter
793 pair(this->phasePairs_[turbulentDispersionModelIter.key()]);
814 if (DByAfs_.found(pair.phase1().name()))
816 this->addField(pair.phase1(),
"DByAf", DByAf1, DByAfs_);
820 if (this->fillFields_)
829 template<
class BasePhaseSystem>
833 const PtrList<volScalarField>& rAUs
836 PtrList<surfaceScalarField> phiKdPhis(this->phaseModels_.size());
842 const phasePair& pair(this->phasePairs_[KdIter.key()]);
857 if (this->fillFields_)
867 return phiKdPhis.xfer();
871 template<
class BasePhaseSystem>
875 const PtrList<surfaceScalarField>& rAUfs
878 PtrList<surfaceScalarField> phiKdPhifs(this->phaseModels_.size());
884 const phasePair& pair(this->phasePairs_[KdfIter.key()]);
892 - rAUfs[iter().index()]*Kf
899 if (this->fillFields_)
909 return phiKdPhifs.xfer();
913 template<
class BasePhaseSystem>
917 const PtrList<volScalarField>& rAUs
920 PtrList<volVectorField> KdUByAs(this->phaseModels_.size());
926 const phasePair& pair(this->phasePairs_[KdIter.key()]);
934 - rAUs[iter().index()]*K*iter.otherPhase().U(),
940 if (this->fillFields_)
945 return KdUByAs.xfer();
949 template<
class BasePhaseSystem>
953 const PtrList<volScalarField>& rAUs,
954 const bool includeVirtualMass
957 PtrList<surfaceScalarField> ddtCorrByAs(this->phaseModels_.size());
960 PtrList<surfaceScalarField> phiCorrs(this->phaseModels_.size());
963 const phaseModel& phase = this->phaseModels_[
phasei];
976 const phaseModel& phase = this->phaseModels_[
phasei];
985 tmp<surfaceScalarField> phiCorrCoeff =
pos0(alphafBar - 0.99);
987 surfaceScalarField::Boundary& phiCorrCoeffBf =
988 phiCorrCoeff.ref().boundaryFieldRef();
995 !this->mesh_.boundary()[
patchi].coupled()
996 || isA<cyclicAMIFvPatch>(this->mesh_.boundary()[
patchi])
999 phiCorrCoeffBf[
patchi] = 0;
1016 if (includeVirtualMass)
1021 const phasePair& pair(this->phasePairs_[VmIter.key()]);
1025 const phaseModel& phase = iter();
1026 const phaseModel& otherPhase = iter.otherPhase();
1034 phiCorrs[phase.index()]
1035 + this->
MRF().absolute(otherPhase.phi())
1037 - phiCorrs[otherPhase.index()]
1045 return ddtCorrByAs.xfer();
1049 template<
class BasePhaseSystem>
1052 const PtrList<volScalarField>& rAUs
1055 Info<<
"Inverting drag systems: ";
1060 PtrList<PtrList<volScalarField>> KdByAs(phases.size());
1061 PtrList<PtrList<surfaceScalarField>> phiKds(phases.size());
1068 new PtrList<volScalarField>(phases.size())
1074 new PtrList<surfaceScalarField>(phases.size())
1081 const phasePair& pair(this->phasePairs_[KdIter.key()]);
1083 const label phase1i = pair.phase1().index();
1084 const label phase2i = pair.phase2().index();
1127 for (
label i = 0; i < phases.size(); ++ i)
1129 for (
label j = i + 1; j < phases.size(); ++ j)
1131 KdByAs[i][j] /= KdByAs[i][i];
1132 phiKds[i][j] /= phiKds[i][i];
1133 for (
label k = i + 1;
k < phases.size(); ++
k)
1135 KdByAs[j][
k] -= KdByAs[j][i]*KdByAs[i][
k];
1136 phiKds[j][
k] -= phiKds[j][i]*phiKds[i][
k];
1143 for (
label i = 1; i < phases.size(); ++ i)
1145 detKdByAs *= KdByAs[i][i];
1146 detPhiKdfs *= phiKds[i][i];
1148 Info<<
"Min cell/face det = " <<
gMin(detKdByAs.primitiveField())
1149 <<
"/" <<
gMin(detPhiKdfs.primitiveField()) <<
endl;
1153 for (
label i = 1; i < phases.size(); ++ i)
1155 if (!phases[i].stationary())
1157 for (
label j = 0; j < i; j ++)
1159 phases[i].URef() -= KdByAs[i][j]*phases[j].U();
1160 phases[i].phiRef() -= phiKds[i][j]*phases[j].phi();
1164 for (
label i = phases.size() - 1; i >= 0; i --)
1166 if (!phases[i].stationary())
1168 for (
label j = phases.size() - 1; j > i; j --)
1170 phases[i].URef() -= KdByAs[i][j]*phases[j].U();
1171 phases[i].phiRef() -= phiKds[i][j]*phases[j].phi();
1173 phases[i].URef() /= KdByAs[i][i];
1174 phases[i].phiRef() /= phiKds[i][i];
1180 template<
class BasePhaseSystem>
1183 const PtrList<surfaceScalarField>& rAUfs
1186 Info<<
"Inverting drag system: ";
1191 PtrList<PtrList<surfaceScalarField>> phiKdfs(phases.size());
1198 new PtrList<surfaceScalarField>(phases.size())
1205 const phasePair& pair(this->phasePairs_[KdfIter.key()]);
1207 const label phase1i = pair.phase1().index();
1208 const label phase2i = pair.phase2().index();
1234 for (
label i = 0; i < phases.size(); ++ i)
1236 for (
label j = i + 1; j < phases.size(); ++ j)
1238 phiKdfs[i][j] /= phiKdfs[i][i];
1239 for (
label k = i + 1;
k < phases.size(); ++
k)
1241 phiKdfs[j][
k] -= phiKdfs[j][i]*phiKdfs[i][
k];
1247 for (
label i = 1; i < phases.size(); ++ i)
1249 detPhiKdfs *= phiKdfs[i][i];
1251 Info<<
"Min face det = " <<
gMin(detPhiKdfs.primitiveField()) <<
endl;
1255 for (
label i = 1; i < phases.size(); ++ i)
1257 if (!phases[i].stationary())
1259 for (
label j = 0; j < i; j ++)
1261 phases[i].phiRef() -= phiKdfs[i][j]*phases[j].phi();
1265 for (
label i = phases.size() - 1; i >= 0; i --)
1267 if (!phases[i].stationary())
1269 for (
label j = phases.size() - 1; j > i; j --)
1271 phases[i].phiRef() -= phiKdfs[i][j]*phases[j].phi();
1273 phases[i].phiRef() /= phiKdfs[i][i];
1279 template<
class BasePhaseSystem>
1287 template<
class BasePhaseSystem>
Class which models interfacial momenum transfer between a number of phases. Drag, virtual mass...
A simple container for copying or transferring objects of type <T>.
surfaceScalarField Vmf("Vmf", fluid.Vmf())
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Type gMin(const FieldField< Field, Type > &f)
HashPtrTable< fvVectorMatrix > momentumTransferTable
virtual autoPtr< phaseSystem::momentumTransferTable > momentumTransferf()
As momentumTransfer, but for the face-based algorithm.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
virtual void partialEliminationf(const PtrList< surfaceScalarField > &rAUfs)
As partialElimination, but for the face-based algorithm. Only solves.
virtual bool read()
Read base phaseProperties dictionary.
Ostream & endl(Ostream &os)
Add newline and flush stream.
MomentumTransferPhaseSystem(const fvMesh &)
Construct from fvMesh.
Calculate the snGrad of the given volField.
label k
Boltzmann constant.
A HashTable specialization for hashing pointers.
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
PtrList< surfaceScalarField > AFfs(fluid.AFfs())
GeometricField< vector, fvPatchField, volMesh > volVectorField
CGAL::Exact_predicates_exact_constructions_kernel K
dimensionedScalar posPart(const dimensionedScalar &ds)
virtual Xfer< PtrList< surfaceScalarField > > phiKdPhifs(const PtrList< surfaceScalarField > &rAUfs) const
As phiKdPhis, but for the face-based algorithm.
virtual Xfer< PtrList< surfaceScalarField > > phiKdPhis(const PtrList< volScalarField > &rAUs) const
Return the explicit drag force fluxes for the cell-based algorithm.
Area-weighted average a surfaceField creating a volField.
virtual Xfer< PtrList< surfaceScalarField > > phiFs(const PtrList< volScalarField > &rAUs)
Return the explicit force fluxes for the cell-based algorithm, that.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Calculate the first temporal derivative.
PtrList< surfaceScalarField > phiFs(fluid.phiFs(rAUs))
virtual const HashPtrTable< surfaceScalarField > & DByAfs() const
Return the phase diffusivities divided by the momentum coefficients.
virtual Xfer< PtrList< surfaceScalarField > > AFfs() const
Return implicit force coefficients on the faces, for the face-based.
Calculate the face-flux of the given field.
static word groupName(Name name, const word &group)
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
virtual Xfer< PtrList< surfaceScalarField > > ddtCorrByAs(const PtrList< volScalarField > &rAUs, const bool includeVirtualMass=false) const
Return the flux corrections for the cell-based algorithm. These.
Calculate the matrix for the first temporal derivative.
PtrListDictionary< phaseModel > phaseModelList
volVectorField F(fluid.F())
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Calculate the divergence of the given field.
dimensionedScalar pos0(const dimensionedScalar &ds)
const dimensionSet dimForce
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
tmp< volScalarField > byDt(const volScalarField &vf)
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
const dimensionSet dimDensity
Calculate the matrix for the divergence of the given field and flux.
virtual ~MomentumTransferPhaseSystem()
Destructor.
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1) *mesh.magSf())
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
fvMatrix< vector > fvVectorMatrix
virtual Xfer< PtrList< volVectorField > > KdUByAs(const PtrList< volScalarField > &rAUs) const
Return the explicit part of the drag force for the cell-based.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
tmp< GeometricField< Type, fvPatchField, volMesh > > average(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Area-weighted average a surfaceField creating a volField.
fluid fillFields("rAU", dimTime/dimDensity, rAUs)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
A class for managing temporary objects.
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
surfaceScalarField Ff(fluid.Ff())
PtrList< surfaceScalarField > phiFfs(fluid.phiFfs(rAUfs))
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Calculate the matrix for implicit and explicit sources.
virtual Xfer< PtrList< surfaceScalarField > > phiFfs(const PtrList< surfaceScalarField > &rAUfs)
As phiFs, but for the face-based algorithm.
dimensionedScalar negPart(const dimensionedScalar &ds)
virtual void partialElimination(const PtrList< volScalarField > &rAUs)
Solve the drag system for the velocities and fluxes.
virtual autoPtr< phaseSystem::momentumTransferTable > momentumTransfer()
Return the momentum transfer matrices for the cell-based algorithm.
const dimensionSet dimVelocity