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().deltaTValue(),
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
547 for (
const phaseModel&
phase2 : phases_)
554 auto iterCvm = Cvms_.cfind(interfacePair(phase,
phase2));
562 iterCvm = Cvms_.cfind(interfacePair(
phase2, phase));
566 tCvm.
ref() += iterCvm()*phase.rho()*
phase2;
577 const phaseModel& phase
587 dimensionSet(1, -2, -2, 0, 0),
592 for (
const phaseModel&
phase2 : phases_)
599 auto Cvm = Cvms_.cfind(interfacePair(phase,
phase2));
607 Cvm = Cvms_.cfind(interfacePair(
phase2, phase));
617 tSvm.ref().boundaryFieldRef();
620 forAll(phase.phi().boundaryField(), patchi)
624 isA<fixedValueFvsPatchScalarField>
626 phase.phi().boundaryField()[patchi]
630 SvmBf[patchi] =
Zero;
645 const multiphaseEuler::dragModel& dm = *iter();
653 dm.phase1()*dm.phase2(),
654 dm.residualPhaseFraction()
660 mag(dm.phase1().U() - dm.phase2().U()),
669 forAll(dm.phase1().phi().boundaryField(), patchi)
673 isA<fixedValueFvsPatchScalarField>
675 dm.phase1().phi().boundaryField()[patchi]
683 dragCoeffsPtr().set(iter.key(), Kptr);
686 return dragCoeffsPtr;
692 const phaseModel& phase,
704 dragModelTable::const_iterator dmIter = dragModels_.begin();
705 dragCoeffFields::const_iterator dcIter =
dragCoeffs.begin();
709 dmIter.good() && dcIter.good();
715 &phase == &dmIter()->phase1()
716 || &phase == &dmIter()->phase2()
719 tdragCoeff.ref() += *dcIter();
739 tSurfaceTension.ref().setOriented();
741 for (
const phaseModel&
phase2 : phases_)
752 tSurfaceTension.ref() +=
762 return tSurfaceTension;
777 for (
const phaseModel& phase : phases_)
789 for (phaseModel& phase : phases_)
803 PtrList<volScalarField> alpha0s(phases_.size());
804 PtrList<surfaceScalarField> alphaPhiSums(phases_.size());
807 for (phaseModel& phase : phases_)
838 subCycleTime alphaSubCycle
843 !(++alphaSubCycle).
end();
849 for (
const phaseModel& phase : phases_)
858 for (phaseModel& phase : phases_)
862 phase.alphaPhi() = alphaPhiSums[
phasei];
888 PtrList<entry> phaseData(lookup(
"phases"));
891 for (phaseModel& phase : phases_)
893 readOK &= phase.read(phaseData[
phasei++].
dict());
896 lookup(
"sigmas") >> sigmas_;
897 lookup(
"interfaceCompression") >> cAlphas_;
898 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) to Type reference.
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.
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 the 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.
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, IOobjectOption::registerOption regOpt, const Mesh &mesh, const dimensionSet &dims, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return tmp field (NO_READ, NO_WRITE) from name, mesh, dimensions and patch type. [Takes current timeN...
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 a time name for the given scalar time value 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
static autoPtr< T > New(Args &&... args)
Construct autoPtr with forwarding arguments.
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.
Do not request registration (bool: false)
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