26 #include "multiphaseSystem.H" 27 #include "alphaContactAngleFvPatchScalarField.H" 32 #include "surfaceInterpolate.H" 41 const Foam::scalar Foam::multiphaseSystem::convertToRad =
47 void Foam::multiphaseSystem::calcAlphas()
52 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
54 alphas_ += level*iter();
60 void Foam::multiphaseSystem::solveAlphas()
62 PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
65 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
67 phaseModel& phase = iter();
75 "phi" + alpha1.name() +
"Corr",
80 "div(phi," + alpha1.name() +
')' 87 forAllIter(PtrDictionary<phaseModel>, phases_, iter2)
89 phaseModel& phase2 = iter2();
92 if (&phase2 == &phase)
continue;
98 cAlphas_.
find(interfacePair(phase, phase2))
101 if (cAlpha != cAlphas_.
end())
105 (
mag(phi_) +
mag(phir))/mesh_.magSf()
113 "div(phir," + alpha2.name() +
',' + alpha1.name() +
')' 124 phase.correctInflowOutflow(alphaPhiCorr);
128 1.0/mesh_.time().deltaT().value(),
150 mesh_.time().timeName(),
159 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
161 phaseModel& phase = iter();
164 alphaPhi += upwind<scalar>(mesh_, phi_).
flux(phase);
165 phase.correctInflowOutflow(alphaPhi);
178 Info<< phase.name() <<
" volume fraction, min, max = " 179 << phase.weightedAverage(mesh_.V()).value()
180 <<
' ' <<
min(phase).value()
181 <<
' ' <<
max(phase).value()
189 Info<<
"Phase-sum volume fraction, min, max = " 190 << sumAlpha.weightedAverage(mesh_.V()).value()
191 <<
' ' <<
min(sumAlpha).value()
192 <<
' ' <<
max(sumAlpha).value()
197 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
199 phaseModel& phase = iter();
201 alpha += alpha*sumCorr;
230 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
241 return nHatfv(alpha1, alpha2) & mesh_.Sf();
251 void Foam::multiphaseSystem::correctContactAngle
254 const phaseModel& phase2,
255 surfaceVectorField::Boundary& nHatb
258 const volScalarField::Boundary& gbf
259 = phase1.boundaryField();
261 const fvBoundaryMesh&
boundary = mesh_.boundary();
265 if (isA<alphaContactAngleFvPatchScalarField>(gbf[
patchi]))
267 const alphaContactAngleFvPatchScalarField& acap =
268 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
274 mesh_.Sf().boundaryField()[
patchi]
275 /mesh_.magSf().boundaryField()[
patchi]
278 alphaContactAngleFvPatchScalarField::thetaPropsTable::
280 acap.thetaProps().find(interfacePair(phase1, phase2));
282 if (tp == acap.thetaProps().end())
285 <<
"Cannot find interface " << interfacePair(phase1, phase2)
286 <<
"\n in table of theta properties for patch " 287 << acap.patch().name()
291 bool matched = (tp.key().first() == phase1.name());
293 scalar theta0 = convertToRad*tp().theta0(matched);
294 scalarField theta(boundary[patchi].size(), theta0);
296 scalar uTheta = tp().uTheta();
301 scalar thetaA = convertToRad*tp().thetaA(matched);
302 scalar thetaR = convertToRad*tp().thetaR(matched);
307 phase1.U().boundaryField()[
patchi].patchInternalField()
308 - phase1.U().boundaryField()[
patchi]
310 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
315 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
319 nWall /= (
mag(nWall) + SMALL);
325 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
339 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
347 nHatPatch = a*AfHatPatch +
b*nHatPatch;
349 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
357 const phaseModel& phase1,
358 const phaseModel& phase2
361 tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
363 correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
366 return -
fvc::div(tnHatfv & mesh_.Sf());
382 "transportProperties",
390 phases_(
lookup(
"phases"), phaseModel::iNew(U.
mesh())),
409 sigmas_(
lookup(
"sigmas")),
410 dimSigma_(1, 0, -2, 0, 0),
411 cAlphas_(
lookup(
"interfaceCompression")),
412 Cvms_(
lookup(
"virtualMass")),
422 interfaceDictTable dragModelsDict(
lookup(
"drag"));
432 *phases_.lookup(iter.key().first()),
433 *phases_.lookup(iter.key().second())
440 const phaseModel& phase1 = iter1();
444 const phaseModel& phase2 = iter2();
446 if (&phase2 != &phase1)
448 scalarCoeffSymmTable::const_iterator
sigma 450 sigmas_.find(interfacePair(phase1, phase2))
453 if (
sigma != sigmas_.end())
455 scalarCoeffSymmTable::const_iterator cAlpha
457 cAlphas_.find(interfacePair(phase1, phase2))
460 if (cAlpha == cAlphas_.end())
463 <<
"Compression coefficient not specified for " 465 << phase1.name() <<
' ' << phase2.name()
466 <<
") for which a surface tension " 467 "coefficient is specified" 483 tmp<volScalarField>
trho = iter()*iter().rho();
486 for (++iter; iter != phases_.end(); ++iter)
488 rho += iter()*iter().rho();
500 tmp<scalarField> trho = iter().boundaryField()[
patchi]*iter().rho().value();
503 for (++iter; iter != phases_.end(); ++iter)
505 rho += iter().boundaryField()[
patchi]*iter().rho().value();
516 tmp<volScalarField> tmu = iter()*(iter().rho()*iter().nu());
519 for (++iter; iter != phases_.end(); ++iter)
521 mu += iter()*(iter().rho()*iter().nu());
533 tmp<scalarField> tmu =
534 iter().boundaryField()[
patchi]
535 *(iter().rho().value()*iter().nu().value());
538 for (++iter; iter != phases_.end(); ++iter)
541 iter().boundaryField()[
patchi]
542 *(iter().rho().value()*iter().nu().value());
545 return tmu/
rho(patchi);
551 const phaseModel& phase
554 tmp<volScalarField> tCvm
561 mesh_.time().timeName(),
568 dimensionSet(1, -3, 0, 0, 0),
576 const phaseModel& phase2 = iter();
578 if (&phase2 != &phase)
580 scalarCoeffTable::const_iterator Cvm
582 Cvms_.find(interfacePair(phase, phase2))
585 if (Cvm != Cvms_.end())
587 tCvm.ref() += Cvm()*phase2.rho()*
phase2;
591 Cvm = Cvms_.find(interfacePair(phase2, phase));
593 if (Cvm != Cvms_.end())
595 tCvm.ref() += Cvm()*phase.rho()*
phase2;
607 const phaseModel& phase
610 tmp<volVectorField> tSvm
617 mesh_.time().timeName(),
624 dimensionSet(1, -2, -2, 0, 0),
632 const phaseModel& phase2 = iter();
634 if (&phase2 != &phase)
636 scalarCoeffTable::const_iterator Cvm
638 Cvms_.find(interfacePair(phase, phase2))
641 if (Cvm != Cvms_.end())
643 tSvm.ref() += Cvm()*phase2.rho()*phase2*phase2.DDtU();
647 Cvm = Cvms_.find(interfacePair(phase2, phase));
649 if (Cvm != Cvms_.end())
651 tSvm.ref() += Cvm()*phase.rho()*phase2*phase2.DDtU();
657 volVectorField::Boundary& SvmBf =
658 tSvm.ref().boundaryFieldRef();
665 isA<fixedValueFvsPatchScalarField>
667 phase.phi().boundaryField()[
patchi]
682 autoPtr<dragCoeffFields> dragCoeffsPtr(
new dragCoeffFields);
686 const dragModel& dm = *iter();
694 dm.phase1()*dm.phase2(),
695 dm.residualPhaseFraction()
701 mag(dm.phase1().U() - dm.phase2().U()),
714 isA<fixedValueFvsPatchScalarField>
716 dm.phase1().phi().boundaryField()[
patchi]
724 dragCoeffsPtr().insert(iter.key(), Kptr);
727 return dragCoeffsPtr;
733 const phaseModel& phase,
737 tmp<volScalarField> tdragCoeff
744 mesh_.time().timeName(),
751 dimensionSet(1, -3, -1, 0, 0),
757 dragModelTable::const_iterator dmIter = dragModels_.begin();
758 dragCoeffFields::const_iterator dcIter = dragCoeffs.begin();
762 dmIter != dragModels_.end() && dcIter != dragCoeffs.end();
768 &phase == &dmIter()->phase1()
769 || &phase == &dmIter()->phase2()
772 tdragCoeff.ref() += *dcIter();
782 const phaseModel& phase1
785 tmp<surfaceScalarField> tSurfaceTension
792 mesh_.time().timeName(),
799 dimensionSet(1, -2, -2, 0, 0),
807 const phaseModel& phase2 = iter();
809 if (&phase2 != &phase1)
811 scalarCoeffSymmTable::const_iterator
sigma 813 sigmas_.find(interfacePair(phase1, phase2))
816 if (
sigma != sigmas_.end())
818 tSurfaceTension.ref() +=
829 return tSurfaceTension;
836 tmp<volScalarField> tnearInt
843 mesh_.time().timeName(),
854 max(tnearInt(),
pos0(iter() - 0.01)*
pos0(0.99 - iter()));
863 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
868 const Time& runTime = mesh_.time();
877 PtrList<volScalarField> alpha0s(phases_.size());
878 PtrList<surfaceScalarField> alphaPhiSums(phases_.size());
881 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
883 phaseModel& phase = iter();
899 "phiSum" + alpha.name(),
913 subCycleTime alphaSubCycle
915 const_cast<Time&>(runTime),
918 !(++alphaSubCycle).end();
924 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
932 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
934 phaseModel& phase = iter();
944 alpha.oldTime() = alpha0s[
phasei];
945 alpha.oldTime().timeIndex() = runTime.timeIndex();
963 PtrList<entry> phaseData(
lookup(
"phases"));
966 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
968 readOK &= iter().read(phaseData[phasei++].
dict());
971 lookup(
"sigmas") >> sigmas_;
972 lookup(
"interfaceCompression") >> cAlphas_;
973 lookup(
"virtualMass") >> Cvms_;
dimensionedScalar tanh(const dimensionedScalar &ds)
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
dimensionedScalar acos(const dimensionedScalar &ds)
#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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
const double e
Elementary charge.
virtual bool read()
Read object.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
IOobject(const word &name, const fileName &instance, const objectRegistry ®istry, readOption r=NO_READ, writeOption w=NO_WRITE, bool registerObject=true)
Construct from name, instance, registry, io options.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< volScalarField > trho
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
autoPtr< multiphaseSystem::dragCoeffFields > dragCoeffs(fluid.dragCoeffs())
Calculate the snGrad of the given volField.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
GeometricField< vector, fvPatchField, volMesh > volVectorField
CGAL::Exact_predicates_exact_constructions_kernel K
tmp< volScalarField > dragCoeff(const phaseModel &phase, const dragCoeffFields &dragCoeffs) const
Return the sum of the drag coefficients for the given phase.
Area-weighted average a surfaceField creating a volField.
const surfaceScalarField & alphaPhi() const
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
tmp< surfaceScalarField > surfaceTension(const phaseModel &phase) const
stressControl lookup("compactNormalStress") >> compactNormalStress
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
static autoPtr< dragModel > New(const dictionary &interfaceDict, const phaseModel &phase1, const phaseModel &phase2)
const phaseModel & phase() const
Return the phase.
dimensionedScalar cos(const dimensionedScalar &ds)
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
void limit(const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phi, surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su, const scalar psiMax, const scalar psiMin, const bool returnCorr)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Calculate the gradient of the given field.
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Calculate the face-flux of the given field.
friend class const_iterator
Declare friendship with the const_iterator.
tmp< volVectorField > U() const
Return the mixture velocity.
const Type & value() const
Return const reference to value.
tmp< volVectorField > Svm(const phaseModel &phase) const
Return the virtual-mass source for the given phase.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
faceListList boundary(nPatches)
void solve()
Solve for the mixture phase-fractions.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
label readLabel(Istream &is)
tmp< volScalarField > rho() const
Return the mixture density.
Calculate the divergence of the given field.
dimensionedScalar pos0(const dimensionedScalar &ds)
tmp< volScalarField > nu() const
Return the mixture laminar viscosity.
bool read()
Read base transportProperties dictionary.
autoPtr< dragCoeffFields > dragCoeffs() const
Return the drag coefficients for all of the interfaces.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
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.
label timeIndex() const
Return the time index of the field.
const dimensionedScalar mu
Atomic mass unit.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
IOdictionary(const IOobject &)
Construct given an IOobject.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
surfaceScalarField phic(mixture.cAlpha() *mag(phi/mesh.magSf()))
#define WarningInFunction
Report a warning using Foam::Warning.
const surfaceScalarField & phi() const
Constant access the mixture flux.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const Time & time() const
Return time.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
MULES: Multidimensional universal limiter for explicit solution.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
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
const objectRegistry & db() const
Return the local objectRegistry.
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux(phi, alpha1, alphaScheme))
multiphaseSystem(const volVectorField &U, const surfaceScalarField &phi)
Construct from components.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
const fvMesh & mesh() const
Constant access the mesh.
tmp< volScalarField > Cvm(const phaseModel &phase) const
Return the virtual-mass coefficient for the given phase.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
friend class const_iterator
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.