26 #include "multiphaseSystem.H" 27 #include "alphaContactAngleFvPatchScalarField.H" 51 const Foam::scalar Foam::multiphaseSystem::convertToRad =
57 void Foam::multiphaseSystem::calcAlphas()
64 alphas_ += level*
phases()[i];
70 void Foam::multiphaseSystem::solveAlphas()
74 PtrList<surfaceScalarField> alphaPhiCorrs(
phases().size());
88 "phi" + alpha1.name() +
"Corr",
93 "div(phi," + alpha1.name() +
')' 102 phaseModel& phase2 =
phases()[phasej];
105 if (&phase2 == &phase)
continue;
109 cAlphaTable::const_iterator cAlpha
111 cAlphas_.find(phasePairKey(phase.name(), phase2.name()))
114 if (cAlpha != cAlphas_.end())
118 (
mag(phi_) +
mag(phir))/mesh_.magSf()
126 "div(phir," + alpha2.name() +
',' + alpha1.name() +
')' 137 surfaceScalarField::Boundary& alphaPhiCorrBf =
138 alphaPhiCorr.boundaryFieldRef();
145 if (!alphaPhiCorrp.coupled())
150 forAll(alphaPhiCorrp, facei)
152 if (phi1p[facei] < 0)
154 alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
178 const scalar rDeltaT = 1.0/mesh_.time().deltaTValue();
203 mesh_.time().timeName(),
219 alphaPhic += upwind<scalar>(mesh_, phi_).
flux(phase);
226 mesh_.time().timeName(),
238 mesh_.time().timeName(),
246 if (phase.divU().valid())
252 if (dgdt[celli] > 0.0)
254 Sp[celli] -= dgdt[celli];
255 Su[celli] += dgdt[celli];
257 else if (dgdt[celli] < 0.0)
261 *(1.0 - alpha[celli])/
max(alpha[celli], 1
e-4);
268 const phaseModel& phase2 =
phases()[phasej];
271 if (&phase2 == &phase)
continue;
273 if (phase2.divU().valid())
279 if (dgdt2[celli] < 0.0)
283 *(1.0 - alpha2[celli])/
max(alpha2[celli], 1
e-4);
287 *alpha[celli]/
max(alpha2[celli], 1
e-4);
289 else if (dgdt2[celli] > 0.0)
291 Sp[celli] -= dgdt2[celli];
306 phase.alphaPhi() += alphaPhic;
308 Info<< phase.name() <<
" volume fraction, min, max = " 309 << phase.weightedAverage(mesh_.V()).value()
310 <<
' ' <<
min(phase).value()
311 <<
' ' <<
max(phase).value()
317 Info<<
"Phase-sum volume fraction, min, max = " 318 << sumAlpha.weightedAverage(mesh_.V()).value()
319 <<
' ' <<
min(sumAlpha).value()
320 <<
' ' <<
max(sumAlpha).value()
347 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
358 return nHatfv(alpha1, alpha2) & mesh_.Sf();
368 void Foam::multiphaseSystem::correctContactAngle
371 const phaseModel& phase2,
372 surfaceVectorField::Boundary& nHatb
375 const volScalarField::Boundary& gbf
376 = phase1.boundaryField();
378 const fvBoundaryMesh&
boundary = mesh_.boundary();
382 if (isA<alphaContactAngleFvPatchScalarField>(gbf[
patchi]))
384 const alphaContactAngleFvPatchScalarField& acap =
385 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
391 mesh_.Sf().boundaryField()[
patchi]
392 /mesh_.magSf().boundaryField()[
patchi]
398 .find(phasePairKey(phase1.name(), phase2.name()));
400 if (tp == acap.thetaProps().end())
403 <<
"Cannot find interface " 404 << phasePairKey(phase1.name(), phase2.name())
405 <<
"\n in table of theta properties for patch " 406 << acap.patch().name()
410 bool matched = (tp.key().first() == phase1.name());
412 scalar theta0 = convertToRad*tp().theta0(matched);
413 scalarField theta(boundary[patchi].size(), theta0);
415 scalar uTheta = tp().uTheta();
420 scalar thetaA = convertToRad*tp().thetaA(matched);
421 scalar thetaR = convertToRad*tp().thetaR(matched);
426 phase1.U()().boundaryField()[
patchi].patchInternalField()
427 - phase1.U()().boundaryField()[
patchi]
429 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
434 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
438 nWall /= (
mag(nWall) + SMALL);
444 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
458 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
466 nHatPatch = a*AfHatPatch +
b*nHatPatch;
468 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
476 const phaseModel& phase1,
477 const phaseModel& phase2
480 tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
482 correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
485 return -
fvc::div(tnHatfv & mesh_.Sf());
512 cAlphas_(
lookup(
"interfaceCompression")),
523 mesh_.setFluxRequired(alphai.name());
538 const phaseModel& phase1
541 tmp<surfaceScalarField> tSurfaceTension
548 mesh_.time().timeName(),
555 dimensionSet(1, -2, -2, 0, 0),
563 const phaseModel& phase2 =
phases()[phasej];
565 if (&phase2 != &phase1)
567 phasePairKey key12(phase1.name(), phase2.name());
569 cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
571 if (cAlpha != cAlphas_.end())
573 tSurfaceTension.ref() +=
583 return tSurfaceTension;
590 tmp<volScalarField> tnearInt
597 mesh_.time().timeName(),
620 const Time& runTime = mesh_.time();
629 tmp<volScalarField> trSubDeltaT;
637 PtrList<volScalarField> alpha0s(
phases().size());
638 PtrList<surfaceScalarField> alphaPhiSums(
phases().size());
658 "phiSum" + alpha.name(),
670 subCycleTime alphaSubCycle
672 const_cast<Time&>(runTime),
675 !(++alphaSubCycle).end();
695 alpha.timeIndex() = runTime.timeIndex();
698 alpha.oldTime() = alpha0s[
phasei];
699 alpha.oldTime().timeIndex() = runTime.timeIndex();
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
const double e
Elementary charge.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
multiphaseSystem::phaseModelList & phases
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
static tmp< volScalarField > localRSubDeltaT(const fvMesh &mesh, const label nAlphaSubCycles)
Calculate and return the reciprocal of the local sub-cycling.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
Calculate the matrix for the laplacian of the field.
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
const Type & value() const
Return const reference to value.
Calculate the snGrad of the given volField.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
DimensionedField< scalar, volMesh > Internal
Type of the internal field from which this GeometricField is derived.
CGAL::Exact_predicates_exact_constructions_kernel K
surfaceScalarField phir(phic *interface.nHatf())
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Calculate the first temporal derivative.
dimensionedScalar pos(const dimensionedScalar &ds)
stressControl lookup("compactNormalStress") >> compactNormalStress
static bool enabled(const fvMesh &mesh)
Return true if LTS is enabled.
dimensionedScalar cos(const dimensionedScalar &ds)
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
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.
Calulate the matrix for the first temporal derivative.
static const volScalarField & localRDeltaT(const fvMesh &mesh)
Return the reciprocal of the local time-step.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
faceListList boundary(nPatches)
Calculate the field for explicit evaluation of implicit and explicit sources.
void solve()
Solve for the mixture phase-fractions.
label readLabel(Istream &is)
const dimensionSet & dimensions() const
Return dimensions.
Calculate the divergence of the given field.
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
defineTypeNameAndDebug(combustionModel, 0)
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.
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const phaseModel & phase() const
Return the phase.
dimensioned< scalar > mag(const dimensioned< Type > &)
volScalarField divU(fvc::div(fvc::absolute(phi, U)))
Field< vector > vectorField
Specialisation of Field<T> for vector.
MULES: Multidimensional universal limiter for explicit solution.
const surfaceScalarField & alphaPhi() const
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.
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
virtual ~multiphaseSystem()
Destructor.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Calculate the matrix for implicit and explicit sources.