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);
176 Info<< phase.name() <<
" volume fraction, min, max = " 177 << phase.weightedAverage(mesh_.V()).value()
178 <<
' ' <<
min(phase).value()
179 <<
' ' <<
max(phase).value()
187 Info<<
"Phase-sum volume fraction, min, max = " 188 << sumAlpha.weightedAverage(mesh_.V()).value()
189 <<
' ' <<
min(sumAlpha).value()
190 <<
' ' <<
max(sumAlpha).value()
195 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
197 phaseModel& phase = iter();
199 alpha += alpha*sumCorr;
228 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
239 return nHatfv(alpha1, alpha2) & mesh_.Sf();
249 void Foam::multiphaseSystem::correctContactAngle
252 const phaseModel& phase2,
253 surfaceVectorField::Boundary& nHatb
256 const volScalarField::Boundary& gbf
257 = phase1.boundaryField();
259 const fvBoundaryMesh&
boundary = mesh_.boundary();
263 if (isA<alphaContactAngleFvPatchScalarField>(gbf[
patchi]))
265 const alphaContactAngleFvPatchScalarField& acap =
266 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
272 mesh_.Sf().boundaryField()[
patchi]
273 /mesh_.magSf().boundaryField()[
patchi]
276 alphaContactAngleFvPatchScalarField::thetaPropsTable::
278 acap.thetaProps().find(interfacePair(phase1, phase2));
280 if (tp == acap.thetaProps().end())
283 <<
"Cannot find interface " << interfacePair(phase1, phase2)
284 <<
"\n in table of theta properties for patch " 285 << acap.patch().name()
289 bool matched = (tp.key().first() == phase1.name());
291 scalar theta0 = convertToRad*tp().theta0(matched);
292 scalarField theta(boundary[patchi].size(), theta0);
294 scalar uTheta = tp().uTheta();
299 scalar thetaA = convertToRad*tp().thetaA(matched);
300 scalar thetaR = convertToRad*tp().thetaR(matched);
305 phase1.U().boundaryField()[
patchi].patchInternalField()
306 - phase1.U().boundaryField()[
patchi]
308 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
313 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
317 nWall /= (
mag(nWall) + small);
323 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
337 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
345 nHatPatch = a*AfHatPatch +
b*nHatPatch;
347 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
355 const phaseModel& phase1,
356 const phaseModel& phase2
359 tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
361 correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
364 return -
fvc::div(tnHatfv & mesh_.Sf());
380 "transportProperties",
388 phases_(
lookup(
"phases"), phaseModel::iNew(U.
mesh())),
407 sigmas_(
lookup(
"sigmas")),
408 dimSigma_(1, 0, -2, 0, 0),
409 cAlphas_(
lookup(
"interfaceCompression")),
410 Cvms_(
lookup(
"virtualMass")),
420 interfaceDictTable dragModelsDict(
lookup(
"drag"));
430 *phases_.lookup(iter.key().first()),
431 *phases_.lookup(iter.key().second())
438 const phaseModel& phase1 = iter1();
442 const phaseModel& phase2 = iter2();
444 if (&phase2 != &phase1)
446 scalarCoeffSymmTable::const_iterator
sigma 448 sigmas_.find(interfacePair(phase1, phase2))
451 if (
sigma != sigmas_.end())
453 scalarCoeffSymmTable::const_iterator cAlpha
455 cAlphas_.find(interfacePair(phase1, phase2))
458 if (cAlpha == cAlphas_.end())
461 <<
"Compression coefficient not specified for " 463 << phase1.name() <<
' ' << phase2.name()
464 <<
") for which a surface tension " 465 "coefficient is specified" 481 tmp<volScalarField>
trho = iter()*iter().rho();
484 for (++iter; iter != phases_.end(); ++iter)
486 rho += iter()*iter().rho();
498 tmp<scalarField> trho = iter().boundaryField()[
patchi]*iter().rho().value();
501 for (++iter; iter != phases_.end(); ++iter)
503 rho += iter().boundaryField()[
patchi]*iter().rho().value();
514 tmp<volScalarField> tmu = iter()*(iter().rho()*iter().nu());
517 for (++iter; iter != phases_.end(); ++iter)
519 mu += iter()*(iter().rho()*iter().nu());
531 tmp<scalarField> tmu =
532 iter().boundaryField()[
patchi]
533 *(iter().rho().value()*iter().nu().value());
536 for (++iter; iter != phases_.end(); ++iter)
539 iter().boundaryField()[
patchi]
540 *(iter().rho().value()*iter().nu().value());
543 return tmu/
rho(patchi);
549 const phaseModel& phase
552 tmp<volScalarField> tCvm
559 mesh_.time().timeName(),
566 dimensionSet(1, -3, 0, 0, 0),
574 const phaseModel& phase2 = iter();
576 if (&phase2 != &phase)
578 scalarCoeffTable::const_iterator Cvm
580 Cvms_.find(interfacePair(phase, phase2))
583 if (Cvm != Cvms_.end())
585 tCvm.ref() += Cvm()*phase2.rho()*
phase2;
589 Cvm = Cvms_.find(interfacePair(phase2, phase));
591 if (Cvm != Cvms_.end())
593 tCvm.ref() += Cvm()*phase.rho()*
phase2;
605 const phaseModel& phase
608 tmp<volVectorField> tSvm
615 mesh_.time().timeName(),
622 dimensionSet(1, -2, -2, 0, 0),
630 const phaseModel& phase2 = iter();
632 if (&phase2 != &phase)
634 scalarCoeffTable::const_iterator Cvm
636 Cvms_.find(interfacePair(phase, phase2))
639 if (Cvm != Cvms_.end())
641 tSvm.ref() += Cvm()*phase2.rho()*phase2*phase2.DDtU();
645 Cvm = Cvms_.find(interfacePair(phase2, phase));
647 if (Cvm != Cvms_.end())
649 tSvm.ref() += Cvm()*phase.rho()*phase2*phase2.DDtU();
655 volVectorField::Boundary& SvmBf =
656 tSvm.ref().boundaryFieldRef();
663 isA<fixedValueFvsPatchScalarField>
665 phase.phi().boundaryField()[
patchi]
680 autoPtr<dragCoeffFields> dragCoeffsPtr(
new dragCoeffFields);
684 const dragModel& dm = *iter();
692 dm.phase1()*dm.phase2(),
693 dm.residualPhaseFraction()
699 mag(dm.phase1().U() - dm.phase2().U()),
712 isA<fixedValueFvsPatchScalarField>
714 dm.phase1().phi().boundaryField()[
patchi]
722 dragCoeffsPtr().insert(iter.key(), Kptr);
725 return dragCoeffsPtr;
731 const phaseModel& phase,
735 tmp<volScalarField> tdragCoeff
742 mesh_.time().timeName(),
749 dimensionSet(1, -3, -1, 0, 0),
755 dragModelTable::const_iterator dmIter = dragModels_.begin();
756 dragCoeffFields::const_iterator dcIter = dragCoeffs.begin();
760 dmIter != dragModels_.end() && dcIter != dragCoeffs.end();
766 &phase == &dmIter()->phase1()
767 || &phase == &dmIter()->phase2()
770 tdragCoeff.ref() += *dcIter();
780 const phaseModel& phase1
783 tmp<surfaceScalarField> tSurfaceTension
790 mesh_.time().timeName(),
797 dimensionSet(1, -2, -2, 0, 0),
805 const phaseModel& phase2 = iter();
807 if (&phase2 != &phase1)
809 scalarCoeffSymmTable::const_iterator
sigma 811 sigmas_.find(interfacePair(phase1, phase2))
814 if (
sigma != sigmas_.end())
816 tSurfaceTension.ref() +=
827 return tSurfaceTension;
834 tmp<volScalarField> tnearInt
841 mesh_.time().timeName(),
852 max(tnearInt(),
pos0(iter() - 0.01)*
pos0(0.99 - iter()));
861 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
866 const Time&
runTime = mesh_.time();
875 PtrList<volScalarField> alpha0s(phases_.size());
876 PtrList<surfaceScalarField> alphaPhiSums(phases_.size());
879 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
881 phaseModel& phase = iter();
897 "phiSum" + alpha.name(),
911 subCycleTime alphaSubCycle
913 const_cast<Time&>(runTime),
916 !(++alphaSubCycle).end();
922 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
930 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
932 phaseModel& phase = iter();
942 alpha.oldTime() = alpha0s[
phasei];
943 alpha.oldTime().timeIndex() = runTime.timeIndex();
961 PtrList<entry> phaseData(
lookup(
"phases"));
964 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
966 readOK &= iter().read(phaseData[phasei++].
dict());
969 lookup(
"sigmas") >> sigmas_;
970 lookup(
"interfaceCompression") >> cAlphas_;
971 lookup(
"virtualMass") >> Cvms_;
void limit(const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phi, surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su, const PsiMaxType &psiMax, const PsiMinType &psiMin, const bool returnCorr)
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
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)
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
Return 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...
const doubleScalar e
Elementary charge.
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
Return 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.