33 #include "surfaceInterpolate.H" 39 #include "alphaContactAngleFvPatchScalarField.H" 43 void Foam::multiphaseMixture::calcAlphas()
48 for (
const phase& ph : phases_)
50 alphas_ += level * ph;
68 "transportProperties",
71 IOobject::MUST_READ_IF_MODIFIED,
76 phases_(lookup(
"phases"), phase::iNew(
U,
phi)),
121 sigmas_(lookup(
"sigmas")),
122 dimSigma_(1, 0, -2, 0, 0),
129 rhoPhi_.setOriented();
141 auto iter = phases_.cbegin();
143 tmp<volScalarField>
trho = iter()*iter().rho();
146 for (++iter; iter != phases_.cend(); ++iter)
148 rho += iter()*iter().rho();
158 auto iter = phases_.cbegin();
160 tmp<scalarField>
trho = iter().boundaryField()[patchi]*iter().rho().value();
163 for (++iter; iter != phases_.cend(); ++iter)
165 rho += iter().boundaryField()[patchi]*iter().rho().value();
175 auto iter = phases_.cbegin();
177 tmp<volScalarField> tmu = iter()*iter().rho()*iter().nu();
180 for (++iter; iter != phases_.cend(); ++iter)
182 mu += iter()*iter().rho()*iter().nu();
192 auto iter = phases_.cbegin();
194 tmp<scalarField> tmu =
196 iter().boundaryField()[patchi]
197 *iter().rho().value()
203 for (++iter; iter != phases_.cend(); ++iter)
207 iter().boundaryField()[patchi]
208 *iter().rho().value()
220 auto iter = phases_.cbegin();
222 tmp<surfaceScalarField> tmuf =
226 for (++iter; iter != phases_.cend(); ++iter)
246 return nu_.boundaryField()[patchi];
260 tmp<surfaceScalarField> tstf
266 "surfaceTensionForce",
267 mesh_.time().timeName(),
280 const phase&
alpha1 = iter1();
284 for (++iter2; iter2 != phases_.cend(); ++iter2)
286 const phase&
alpha2 = iter2();
293 <<
"Cannot find interface " << interfacePair(
alpha1,
alpha2)
294 <<
" in list of sigma values" 342 !(++alphaSubCycle).
end();
363 for (phase& ph : phases_)
392 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
413 void Foam::multiphaseMixture::correctContactAngle
423 const fvBoundaryMesh&
boundary = mesh_.boundary();
427 if (isA<alphaContactAngleFvPatchScalarField>(gb1f[patchi]))
429 const alphaContactAngleFvPatchScalarField& acap =
430 refCast<const alphaContactAngleFvPatchScalarField>(gb1f[patchi]);
432 correctBoundaryContactAngle(acap, patchi,
alpha1,
alpha2, nHatb);
434 else if (isA<alphaContactAngleFvPatchScalarField>(gb2f[patchi]))
436 const alphaContactAngleFvPatchScalarField& acap =
437 refCast<const alphaContactAngleFvPatchScalarField>(gb2f[patchi]);
439 correctBoundaryContactAngle(acap, patchi,
alpha2,
alpha1, nHatb);
445 void Foam::multiphaseMixture::correctBoundaryContactAngle
447 const alphaContactAngleFvPatchScalarField& acap,
454 const fvBoundaryMesh&
boundary = mesh_.boundary();
460 mesh_.Sf().boundaryField()[patchi]
461 /mesh_.magSf().boundaryField()[patchi]
464 const auto tp = acap.thetaProps().cfind(interfacePair(
alpha1,
alpha2));
469 <<
"Cannot find interface " << interfacePair(
alpha1,
alpha2)
470 <<
"\n in table of theta properties for patch " 471 << acap.patch().name()
475 const bool matched = (tp.key().first() ==
alpha1.
name());
477 const scalar theta0 =
degToRad(tp().theta0(matched));
480 const scalar uTheta = tp().uTheta();
485 const scalar thetaA =
degToRad(tp().thetaA(matched));
486 const scalar thetaR =
degToRad(tp().thetaR(matched));
491 U_.boundaryField()[patchi].patchInternalField()
492 - U_.boundaryField()[patchi]
494 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
499 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
503 nWall /= (
mag(nWall) + SMALL);
509 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
523 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
531 nHatPatch = a*AfHatPatch +
b*nHatPatch;
533 nHatPatch /= (
mag(nHatPatch) + deltaN_.value());
543 tmp<surfaceVectorField> tnHatfv = nHatfv(
alpha1,
alpha2);
548 return -
fvc::div(tnHatfv & mesh_.Sf());
555 tmp<volScalarField> tnearInt
562 mesh_.time().timeName(),
570 for (
const phase& ph : phases_)
572 tnearInt.ref() =
max(tnearInt(),
pos0(ph - 0.01)*
pos0(0.99 - ph));
579 void Foam::multiphaseMixture::solveAlphas
584 static label nSolves(-1);
587 const word alphaScheme(
"div(phi,alpha)");
593 PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
596 for (phase&
alpha : phases_)
615 for (phase&
alpha2 : phases_)
631 1.0/mesh_.time().deltaT().value(),
655 mesh_.time().timeName(),
664 for (phase&
alpha : phases_)
679 <<
alpha.weightedAverage(mesh_.V()).value()
689 Info<<
"Phase-sum volume fraction, min, max = " 690 << sumAlpha.weightedAverage(mesh_.V()).value()
691 <<
' ' <<
min(sumAlpha).value()
692 <<
' ' <<
max(sumAlpha).value()
697 for (phase&
alpha : phases_)
712 PtrList<entry> phaseData(lookup(
"phases"));
715 for (phase& ph : phases_)
717 readOK &= ph.read(phaseData[
phasei++].
dict());
720 readEntry(
"sigmas", sigmas_);
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)
tmp< surfaceScalarField > muf() const
Return the face-interpolated dynamic laminar viscosity.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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.
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.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Unit conversion functions.
word alpharScheme("div(phirb,alpha)")
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< volScalarField > trho
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
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.
GeometricField< vector, fvPatchField, volMesh > volVectorField
CGAL::Exact_predicates_exact_constructions_kernel K
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
multiphaseMixture(const volVectorField &U, const surfaceScalarField &phi)
Construct from components.
#define forAll(list, i)
Loop across all elements in list.
tmp< volScalarField > rho() const
Return the mixture density.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
const dimensionedScalar e
Elementary charge.
dimensionedScalar cos(const dimensionedScalar &ds)
bool read()
Read base transportProperties dictionary.
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.
void correct()
Correct the mixture properties.
const Time & time() const noexcept
Return time registry.
Calculate the face-flux of the given field.
dimensionedScalar cbrt(const dimensionedScalar &ds)
void setOriented(bool on=true) noexcept
Set the oriented flag: on/off.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
tmp< volScalarField > nu() const
Return the kinematic laminar viscosity.
Calculate the divergence of the given field.
dimensionedScalar pos0(const dimensionedScalar &ds)
static word timeName(const scalar t, const int precision=precision_)
Return time name of given scalar time formatted with the given precision.
tmp< volScalarField > mu() const
Return the dynamic laminar viscosity.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil) - psi) *pSat, rhoMin);# 1 "/home/chef2/andy/OpenFOAM/release/v2212/OpenFOAM-v2212/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
void solve()
Solve for the mixture phase-fractions.
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 word & name() const noexcept
Return const reference to name.
surfaceScalarField phic(mixture.cAlpha() *mag(alphaPhic/mesh.magSf()))
const dimensionedScalar mu
Atomic mass unit.
tmp< surfaceScalarField > nuf() const
Return the face-interpolated dynamic laminar viscosity.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
messageStream Info
Information stream (stdout output on master, null elsewhere)
virtual bool read()=0
Read transportProperties dictionary.
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.
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
A class for managing temporary objects.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux(phi, alpha1, alphaScheme))
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)
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)
forAllConstIters(mixture.phases(), phase)
tmp< surfaceScalarField > surfaceTensionForce() const
static constexpr const zero Zero
Global zero (0)
const volScalarField & alpha1