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& phase1 = iter();
78 "phi" + alpha1.name() +
"Corr",
83 "div(phi," + alpha1.name() +
')' 90 forAllIter(PtrDictionary<phaseModel>, phases_, iter2)
92 phaseModel& phase2 = iter2();
95 if (&phase2 == &phase1)
continue;
101 cAlphas_.
find(interfacePair(phase1, phase2))
104 if (cAlpha != cAlphas_.
end())
108 (
mag(phi_) +
mag(phir))/mesh_.magSf()
116 "div(phir," + alpha2.name() +
',' + alpha1.name() +
')' 127 surfaceScalarField::Boundary& alphaPhiCorrBf =
128 alphaPhiCorr.boundaryFieldRef();
135 if (!alphaPhiCorrp.coupled())
140 forAll(alphaPhiCorrp, facei)
142 if (phi1p[facei] < 0)
144 alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
152 1.0/mesh_.time().deltaT().value(),
174 mesh_.time().timeName(),
183 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
185 phaseModel& phase1 = iter();
188 alphaPhi += upwind<scalar>(mesh_, phi_).
flux(phase1);
201 Info<< phase1.name() <<
" volume fraction, min, max = " 202 << phase1.weightedAverage(mesh_.V()).value()
203 <<
' ' <<
min(phase1).value()
204 <<
' ' <<
max(phase1).value()
212 Info<<
"Phase-sum volume fraction, min, max = " 213 << sumAlpha.weightedAverage(mesh_.V()).value()
214 <<
' ' <<
min(sumAlpha).value()
215 <<
' ' <<
max(sumAlpha).value()
244 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
255 return nHatfv(alpha1, alpha2) & mesh_.Sf();
265 void Foam::multiphaseSystem::correctContactAngle
267 const phaseModel& phase1,
268 const phaseModel& phase2,
269 surfaceVectorField::Boundary& nHatb
272 const volScalarField::Boundary& gbf
273 = phase1.boundaryField();
275 const fvBoundaryMesh&
boundary = mesh_.boundary();
279 if (isA<alphaContactAngleFvPatchScalarField>(gbf[
patchi]))
281 const alphaContactAngleFvPatchScalarField& acap =
282 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
288 mesh_.Sf().boundaryField()[
patchi]
289 /mesh_.magSf().boundaryField()[
patchi]
294 acap.thetaProps().find(interfacePair(phase1, phase2));
296 if (tp == acap.thetaProps().end())
299 <<
"Cannot find interface " << interfacePair(phase1, phase2)
300 <<
"\n in table of theta properties for patch " 301 << acap.patch().name()
305 bool matched = (tp.key().first() == phase1.name());
307 scalar theta0 = convertToRad*tp().theta0(matched);
310 scalar uTheta = tp().uTheta();
315 scalar thetaA = convertToRad*tp().thetaA(matched);
316 scalar thetaR = convertToRad*tp().thetaR(matched);
321 phase1.U().boundaryField()[
patchi].patchInternalField()
322 - phase1.U().boundaryField()[
patchi]
324 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
329 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
333 nWall /= (
mag(nWall) + SMALL);
339 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
353 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
361 nHatPatch = a*AfHatPatch +
b*nHatPatch;
363 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
371 const phaseModel& phase1,
372 const phaseModel& phase2
375 tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
377 correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
380 return -
fvc::div(tnHatfv & mesh_.Sf());
396 "transportProperties",
404 phases_(
lookup(
"phases"), phaseModel::iNew(U.
mesh())),
423 sigmas_(
lookup(
"sigmas")),
424 dimSigma_(1, 0, -2, 0, 0),
425 cAlphas_(
lookup(
"interfaceCompression")),
426 Cvms_(
lookup(
"virtualMass")),
436 interfaceDictTable dragModelsDict(
lookup(
"drag"));
446 *phases_.lookup(iter.key().first()),
447 *phases_.lookup(iter.key().second())
454 const phaseModel& phase1 = iter1();
458 const phaseModel& phase2 = iter2();
460 if (&phase2 != &phase1)
462 scalarCoeffSymmTable::const_iterator
sigma 464 sigmas_.find(interfacePair(phase1, phase2))
467 if (
sigma != sigmas_.end())
469 scalarCoeffSymmTable::const_iterator cAlpha
471 cAlphas_.find(interfacePair(phase1, phase2))
474 if (cAlpha == cAlphas_.end())
477 <<
"Compression coefficient not specified for " 479 << phase1.name() <<
' ' << phase2.name()
480 <<
") for which a surface tension " 481 "coefficient is specified" 497 tmp<volScalarField>
trho = iter()*iter().rho();
500 for (++iter; iter != phases_.end(); ++iter)
502 rho += iter()*iter().rho();
514 tmp<scalarField> trho = iter().boundaryField()[
patchi]*iter().rho().value();
517 for (++iter; iter != phases_.end(); ++iter)
519 rho += iter().boundaryField()[
patchi]*iter().rho().value();
530 tmp<volScalarField> tmu = iter()*(iter().rho()*iter().nu());
533 for (++iter; iter != phases_.end(); ++iter)
535 mu += iter()*(iter().rho()*iter().nu());
547 tmp<scalarField> tmu =
548 iter().boundaryField()[
patchi]
549 *(iter().rho().value()*iter().nu().value());
552 for (++iter; iter != phases_.end(); ++iter)
555 iter().boundaryField()[
patchi]
556 *(iter().rho().value()*iter().nu().value());
559 return tmu/
rho(patchi);
565 const phaseModel&
phase 568 tmp<volScalarField> tCvm
575 mesh_.time().timeName(),
582 dimensionSet(1, -3, 0, 0, 0),
590 const phaseModel& phase2 = iter();
592 if (&phase2 != &phase)
594 scalarCoeffTable::const_iterator Cvm
596 Cvms_.find(interfacePair(phase, phase2))
599 if (Cvm != Cvms_.end())
601 tCvm.ref() += Cvm()*phase2.rho()*
phase2;
605 Cvm = Cvms_.find(interfacePair(phase2, phase));
607 if (Cvm != Cvms_.end())
609 tCvm.ref() += Cvm()*phase.rho()*
phase2;
621 const phaseModel& phase
624 tmp<volVectorField> tSvm
631 mesh_.time().timeName(),
638 dimensionSet(1, -2, -2, 0, 0),
646 const phaseModel& phase2 = iter();
648 if (&phase2 != &phase)
650 scalarCoeffTable::const_iterator Cvm
652 Cvms_.find(interfacePair(phase, phase2))
655 if (Cvm != Cvms_.end())
657 tSvm.ref() += Cvm()*phase2.rho()*phase2*phase2.DDtU();
661 Cvm = Cvms_.find(interfacePair(phase2, phase));
663 if (Cvm != Cvms_.end())
665 tSvm.ref() += Cvm()*phase.rho()*phase2*phase2.DDtU();
671 volVectorField::Boundary& SvmBf =
672 tSvm.ref().boundaryFieldRef();
679 isA<fixedValueFvsPatchScalarField>
681 phase.phi().boundaryField()[
patchi]
696 autoPtr<dragCoeffFields> dragCoeffsPtr(
new dragCoeffFields);
700 const dragModel& dm = *iter();
708 dm.phase1()*dm.phase2(),
709 dm.residualPhaseFraction()
715 mag(dm.phase1().U() - dm.phase2().U()),
728 isA<fixedValueFvsPatchScalarField>
730 dm.phase1().phi().boundaryField()[
patchi]
738 dragCoeffsPtr().insert(iter.key(), Kptr);
741 return dragCoeffsPtr;
747 const phaseModel& phase,
751 tmp<volScalarField> tdragCoeff
758 mesh_.time().timeName(),
765 dimensionSet(1, -3, -1, 0, 0),
771 dragModelTable::const_iterator dmIter = dragModels_.begin();
772 dragCoeffFields::const_iterator dcIter = dragCoeffs.begin();
776 dmIter != dragModels_.end() && dcIter != dragCoeffs.end();
782 &phase == &dmIter()->phase1()
783 || &phase == &dmIter()->phase2()
786 tdragCoeff.ref() += *dcIter();
796 const phaseModel& phase1
799 tmp<surfaceScalarField> tSurfaceTension
806 mesh_.time().timeName(),
813 dimensionSet(1, -2, -2, 0, 0),
821 const phaseModel& phase2 = iter();
823 if (&phase2 != &phase1)
825 scalarCoeffSymmTable::const_iterator
sigma 827 sigmas_.find(interfacePair(phase1, phase2))
830 if (
sigma != sigmas_.end())
832 tSurfaceTension.ref() +=
843 return tSurfaceTension;
850 tmp<volScalarField> tnearInt
857 mesh_.time().timeName(),
867 tnearInt.ref() =
max(tnearInt(),
pos(iter() - 0.01)*
pos(0.99 - iter()));
876 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
881 const Time& runTime = mesh_.time();
890 PtrList<volScalarField> alpha0s(phases_.size());
891 PtrList<surfaceScalarField> alphaPhiSums(phases_.size());
894 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
896 phaseModel& phase = iter();
912 "phiSum" + alpha.name(),
926 subCycleTime alphaSubCycle
928 const_cast<Time&>(runTime),
931 !(++alphaSubCycle).end();
937 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
945 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
947 phaseModel& phase = iter();
957 alpha.oldTime() = alpha0s[
phasei];
958 alpha.oldTime().timeIndex() = runTime.timeIndex();
976 PtrList<entry> phaseData(
lookup(
"phases"));
979 forAllIter(PtrDictionary<phaseModel>, phases_, iter)
981 readOK &= iter().read(phaseData[phasei++].
dict());
984 lookup(
"sigmas") >> sigmas_;
985 lookup(
"interfaceCompression") >> cAlphas_;
986 lookup(
"virtualMass") >> Cvms_;
tmp< volScalarField > nu() const
Return the mixture laminar viscosity.
fvsPatchField< scalar > fvsPatchScalarField
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.
tmp< surfaceScalarField > surfaceTension(const phaseModel &phase) const
const objectRegistry & db() const
Return the local objectRegistry.
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< volVectorField > U() const
Return the mixture velocity.
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< dragCoeffFields > dragCoeffs() const
Return the drag coefficients for all of the interfaces.
autoPtr< multiphaseSystem::dragCoeffFields > dragCoeffs(fluid.dragCoeffs())
const Type & value() const
Return const reference to value.
tmp< volScalarField > Cvm(const phaseModel &phase) const
Return the virtual-mass coefficient for the given phase.
Calculate the snGrad of the given volField.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
label size() const
Return number of elements in list.
GeometricField< vector, fvPatchField, volMesh > volVectorField
CGAL::Exact_predicates_exact_constructions_kernel K
Area-weighted average a surfaceField creating a volField.
surfaceScalarField phir(phic *interface.nHatf())
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
const surfaceScalarField & phi() const
Constant access the mixture flux.
tmp< volScalarField > dragCoeff(const phaseModel &phase, const dragCoeffFields &dragCoeffs) const
Return the sum of the drag coefficients for the given phase.
dimensionedScalar pos(const dimensionedScalar &ds)
surfaceScalarField alphaPhi(IOobject("alphaPhi", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), phi *fvc::interpolate(alpha1))
stressControl lookup("compactNormalStress") >> compactNormalStress
static autoPtr< dragModel > New(const dictionary &interfaceDict, const phaseModel &phase1, const phaseModel &phase2)
dimensionedScalar cos(const dimensionedScalar &ds)
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
const fvMesh & mesh() const
Constant access the mesh.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Calculate the gradient of the given field.
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
const Time & time() const
Return time.
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.
label timeIndex() const
Return the time index of the field.
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.
tmp< volVectorField > Svm(const phaseModel &phase) const
Return the virtual-mass source for the given phase.
bool read()
Read base transportProperties dictionary.
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.
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)
#define WarningInFunction
Report a warning using Foam::Warning.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const phaseModel & phase() const
Return the phase.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
MULES: Multidimensional universal limiter for explicit solution.
const surfaceScalarField & alphaPhi() const
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
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)
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
multiphaseSystem(const volVectorField &U, const surfaceScalarField &phi)
Construct from components.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
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.