30 #include "alphaContactAngleFvPatchScalarField.H" 35 #include "surfaceInterpolate.H" 46 void Foam::multiphaseSystem::calcAlphas()
51 for (
const phaseModel& phase : phases_)
53 alphas_ += level * phase;
59 void Foam::multiphaseSystem::solveAlphas()
61 PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
64 for (phaseModel& phase : phases_)
85 for (phaseModel&
phase2 : phases_)
89 if (&
phase2 == &phase)
continue;
93 const auto cAlpha = cAlphas_.cfind(interfacePair(phase,
phase2));
105 const word phirScheme
122 1.0/mesh_.time().deltaT().value(),
144 mesh_.time().timeName(),
153 for (phaseModel& phase : phases_)
168 Info<< phase.
name() <<
" volume fraction, min, max = " 170 <<
' ' <<
min(phase).value()
171 <<
' ' <<
max(phase).value()
179 Info<<
"Phase-sum volume fraction, min, max = " 180 << sumAlpha.weightedAverage(mesh_.V()).value()
181 <<
' ' <<
min(sumAlpha).value()
182 <<
' ' <<
max(sumAlpha).value()
187 for (phaseModel& phase : phases_)
219 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
240 void Foam::multiphaseSystem::correctContactAngle
250 const fvBoundaryMesh&
boundary = mesh_.boundary();
256 isA<multiphaseEuler::alphaContactAngleFvPatchScalarField>
265 const multiphaseEuler::alphaContactAngleFvPatchScalarField
275 mesh_.Sf().boundaryField()[patchi]
276 /mesh_.magSf().boundaryField()[patchi]
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 const scalar theta0 =
degToRad(tp().theta0(matched));
296 scalar uTheta = tp().uTheta();
301 const scalar thetaA =
degToRad(tp().thetaA(matched));
302 const scalar thetaR =
degToRad(tp().thetaR(matched));
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());
361 tmp<surfaceVectorField> tnHatfv = nHatfv(
phase1,
phase2);
382 "transportProperties",
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())
438 for (
const phaseModel&
phase1 : phases_)
440 for (
const phaseModel&
phase2 : phases_)
452 <<
"Compression coefficient not specified for phase pair (" 454 <<
") for which a surface tension coefficient is specified" 466 auto iter = phases_.cbegin();
468 tmp<volScalarField>
trho = iter()*iter().rho();
471 for (++iter; iter != phases_.cend(); ++iter)
473 rho += iter()*iter().rho();
483 auto iter = phases_.cbegin();
485 tmp<scalarField>
trho = iter().boundaryField()[patchi]*iter().rho().value();
488 for (++iter; iter != phases_.cend(); ++iter)
490 rho += iter().boundaryField()[patchi]*iter().rho().value();
499 auto iter = phases_.cbegin();
501 tmp<volScalarField> tmu = iter()*(iter().rho()*iter().nu());
504 for (++iter; iter != phases_.cend(); ++iter)
506 mu += iter()*(iter().rho()*iter().nu());
516 auto iter = phases_.cbegin();
518 tmp<scalarField> tmu =
519 iter().boundaryField()[patchi]
520 *(iter().rho().value()*iter().nu().value());
523 for (++iter; iter != phases_.cend(); ++iter)
526 iter().boundaryField()[patchi]
527 *(iter().rho().value()*iter().nu().value());
530 return tmu/
rho(patchi);
536 const phaseModel& phase
539 tmp<volScalarField> tCvm
546 mesh_.time().timeName(),
554 for (
const phaseModel&
phase2 : phases_)
561 auto iterCvm = Cvms_.cfind(interfacePair(phase,
phase2));
569 iterCvm = Cvms_.cfind(interfacePair(
phase2, phase));
573 tCvm.
ref() += iterCvm()*phase.rho()*
phase2;
584 const phaseModel& phase
587 tmp<volVectorField> tSvm
594 mesh_.time().timeName(),
601 dimensionSet(1, -2, -2, 0, 0),
607 for (
const phaseModel&
phase2 : phases_)
614 auto Cvm = Cvms_.cfind(interfacePair(phase,
phase2));
622 Cvm = Cvms_.cfind(interfacePair(
phase2, phase));
632 tSvm.ref().boundaryFieldRef();
635 forAll(phase.phi().boundaryField(), patchi)
639 isA<fixedValueFvsPatchScalarField>
641 phase.phi().boundaryField()[patchi]
645 SvmBf[patchi] =
Zero;
656 autoPtr<dragCoeffFields> dragCoeffsPtr(
new dragCoeffFields);
660 const multiphaseEuler::dragModel& dm = *iter();
668 dm.phase1()*dm.phase2(),
669 dm.residualPhaseFraction()
675 mag(dm.phase1().U() - dm.phase2().U()),
684 forAll(dm.phase1().phi().boundaryField(), patchi)
688 isA<fixedValueFvsPatchScalarField>
690 dm.phase1().phi().boundaryField()[patchi]
698 dragCoeffsPtr().set(iter.key(), Kptr);
701 return dragCoeffsPtr;
707 const phaseModel& phase,
711 tmp<volScalarField> tdragCoeff
718 mesh_.time().timeName(),
725 dimensionSet(1, -3, -1, 0, 0),
731 dragModelTable::const_iterator dmIter = dragModels_.begin();
732 dragCoeffFields::const_iterator dcIter =
dragCoeffs.begin();
736 dmIter.good() && dcIter.good();
742 &phase == &dmIter()->phase1()
743 || &phase == &dmIter()->phase2()
746 tdragCoeff.ref() += *dcIter();
759 tmp<surfaceScalarField> tSurfaceTension
766 mesh_.time().timeName(),
773 dimensionSet(1, -2, -2, 0, 0),
778 tSurfaceTension.ref().setOriented();
780 for (
const phaseModel&
phase2 : phases_)
791 tSurfaceTension.ref() +=
801 return tSurfaceTension;
808 tmp<volScalarField> tnearInt
815 mesh_.time().timeName(),
823 for (
const phaseModel& phase : phases_)
835 for (phaseModel& phase : phases_)
849 PtrList<volScalarField> alpha0s(phases_.size());
850 PtrList<surfaceScalarField> alphaPhiSums(phases_.size());
853 for (phaseModel& phase : phases_)
884 subCycleTime alphaSubCycle
889 !(++alphaSubCycle).
end();
895 for (
const phaseModel& phase : phases_)
904 for (phaseModel& phase : phases_)
908 phase.alphaPhi() = alphaPhiSums[
phasei];
934 PtrList<entry> phaseData(lookup(
"phases"));
937 for (phaseModel& phase : phases_)
939 readOK &= phase.read(phaseData[
phasei++].
dict());
942 lookup(
"sigmas") >> sigmas_;
943 lookup(
"interfaceCompression") >> cAlphas_;
944 lookup(
"virtualMass") >> Cvms_;
tmp< volScalarField > nu() const
Return the mixture laminar viscosity.
ITstream & lookup(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return an entry data stream. FatalIOError if not found, or not a stream. ...
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)
void correctInflowOutflow(surfaceScalarField &alphaPhi) const
Ensure that the flux at inflow/outflow BCs is preserved.
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
dimensionedScalar acos(const dimensionedScalar &ds)
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &f1)
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
virtual bool read()
Read object.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
bool found(const Key &key) const
Same as contains()
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
const word & name() const noexcept
Return the object name.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Unit conversion functions.
Type & refCast(U &obj)
A dynamic_cast (for references). Generates a FatalError on failed casts and uses the virtual type() m...
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
autoPtr< multiphaseSystem::dragCoeffFields > dragCoeffs(fluid.dragCoeffs())
Calculate the snGrad of the given volField.
const volScalarField & alpha2
GeometricBoundaryField< vector, fvsPatchField, surfaceMesh > Boundary
Type of boundary fields.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
const dimensionSet dimless
Dimensionless.
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T. FatalIOError if not found, or if the number of tokens is incorrect.
Lookup type of boundary radiation properties.
GeometricField< vector, fvPatchField, volMesh > volVectorField
CGAL::Exact_predicates_exact_constructions_kernel K
const surfaceScalarField & phi() const
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
const volVectorField & DDtU() const
Area-weighted average a surfaceField creating a volField.
tmp< volScalarField > rho() const
Return the mixture density.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
#define forAll(list, i)
Loop across all elements in list.
const surfaceScalarField & alphaPhi() const
GeometricField< scalar, fvPatchField, volMesh > volScalarField
const volVectorField & U() const
tmp< volScalarField > Cvm(const phaseModel &phase) const
Return the virtual-mass coefficient for the given phase.
const dimensionedScalar e
Elementary charge.
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.
const Time & time() const noexcept
Return time registry.
Calculate the face-flux of the given field.
bool set(const Key &key, T *ptr)
Assign a new entry, overwrites existing.
dimensionedScalar cbrt(const dimensionedScalar &ds)
virtual bool write(const bool writeOnProc=true) const
Write using setting from DB.
label timeIndex() const noexcept
Return current time index.
tmp< volScalarField > dragCoeff(const phaseModel &phase, const dragCoeffFields &dragCoeffs) const
Return the sum of the drag coefficients for the given phase.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
void solve()
Solve for the mixture phase-fractions.
Calculate the divergence of the given field.
dimensionedScalar pos0(const dimensionedScalar &ds)
tmp< volVectorField > Svm(const phaseModel &phase) const
Return the virtual-mass source for the given phase.
static word timeName(const scalar t, const int precision=precision_)
Return time name of given scalar time formatted with the given precision.
bool read()
Read base transportProperties dictionary.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
static autoPtr< dragModel > New(const dictionary &dict, const phaseModel &phase1, const phaseModel &phase2)
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
const word & name() const noexcept
Return const reference to name.
surfaceScalarField phic(mixture.cAlpha() *mag(alphaPhic/mesh.magSf()))
const dimensionedScalar mu
Atomic mass unit.
const dimensionedScalar & rho() const
#define WarningInFunction
Report a warning using Foam::Warning.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
Single incompressible phase derived from the phase-fraction. Used as part of the multiPhaseMixture fo...
const word & name() const
messageStream Info
Information stream (stdout output on master, null elsewhere)
Internal & ref(const bool updateAccessTime=true)
Same as internalFieldRef()
Field< vector > vectorField
Specialisation of Field<T> for vector.
MULES: Multidimensional universal limiter for explicit solution.
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
dimensioned< Type > weightedAverage(const DimensionedField< scalar, GeoMesh > &weightField) const
Calculate and return weighted average.
A class for managing temporary objects.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
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: [].
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
Defines the attributes of an object for which implicit objectRegistry management is supported...
Single incompressible phase derived from the phase-fraction. Used as part of the multiPhaseMixture fo...
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
dimensionedScalar deltaT() const
Return time step.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
autoPtr< dragCoeffFields > dragCoeffs() const
Return the drag coefficients for all of the interfaces.
forAllConstIters(mixture.phases(), phase)
tmp< surfaceScalarField > surfaceTension(const phaseModel &phase) const
static constexpr const zero Zero
Global zero (0)
const volScalarField & alpha1