26 #include "multiphaseSystem.H" 27 #include "alphaContactAngleFvPatchScalarField.H" 52 const Foam::scalar Foam::multiphaseSystem::convertToRad =
58 void Foam::multiphaseSystem::calcAlphas()
65 alphas_ += level*
phases()[i];
71 void Foam::multiphaseSystem::solveAlphas()
84 mesh_.time().timeName(),
90 forAll(stationaryPhases(), stationaryPhasei)
92 alphaVoid -= stationaryPhases()[stationaryPhasei];
106 upwind<scalar>(mesh_, phi_).interpolate(phase)
112 PtrList<surfaceScalarField> alphaPhiCorrs(
phases().size());
113 forAll(stationaryPhases(), stationaryPhasei)
115 phaseModel& phase = stationaryPhases()[stationaryPhasei];
123 - upwind<scalar>(mesh_, phi_).flux(phase)
127 forAll(movingPhases(), movingPhasei)
129 phaseModel& phase = movingPhases()[movingPhasei];
138 fvc::flux(phi_, phase,
"div(phi," + alpha.name() +
')')
149 if (&phase2 == &phase)
continue;
153 cAlphaTable::const_iterator cAlpha
155 cAlphas_.find(phasePairKey(phase.name(), phase2.name()))
158 if (cAlpha != cAlphas_.end())
162 (
mag(phi_) +
mag(phir))/mesh_.magSf()
170 "div(phir," + alpha2.name() +
',' + alpha.name() +
')' 181 phase.correctInflowOutflow(alphaPhiCorr);
191 min(alphaVoid.primitiveField(), phase.alphaMax())(),
199 forAll(stationaryPhases(), stationaryPhasei)
201 fixedAlphaPhiCorrs.insert(stationaryPhases()[stationaryPhasei].index());
206 forAll(movingPhases(), movingPhasei)
208 phaseModel& phase = movingPhases()[movingPhasei];
212 alphaPhi += upwind<scalar>(mesh_, phi_).
flux(phase);
213 phase.correctInflowOutflow(alphaPhi);
220 mesh_.time().timeName(),
233 if (phase.divU().valid())
239 if (dgdt[celli] > 0.0)
241 Sp[celli] -= dgdt[celli];
242 Su[celli] += dgdt[celli];
244 else if (dgdt[celli] < 0.0)
248 *(1 - alpha[celli])/
max(alpha[celli], 1
e-4);
255 const phaseModel& phase2 =
phases()[phasej];
258 if (&phase2 == &phase)
continue;
260 if (phase2.divU().valid())
266 if (dgdt2[celli] < 0.0)
270 *(1 - alpha2[celli])/
max(alpha2[celli], 1
e-4);
274 *alpha[celli]/
max(alpha2[celli], 1
e-4);
276 else if (dgdt2[celli] > 0.0)
278 Sp[celli] -= dgdt2[celli];
301 Info<< phase.name() <<
" fraction, min, max = " 302 << phase.weightedAverage(mesh_.V()).value()
303 <<
' ' <<
min(phase).value()
304 <<
' ' <<
max(phase).value()
313 mesh_.time().timeName(),
319 forAll(movingPhases(), movingPhasei)
321 sumAlphaMoving += movingPhases()[movingPhasei];
324 Info<<
"Phase-sum volume fraction, min, max = " 325 << (sumAlphaMoving + 1 - alphaVoid)().weightedAverage(mesh_.V()).value()
326 <<
' ' <<
min(sumAlphaMoving + 1 - alphaVoid).value()
327 <<
' ' <<
max(sumAlphaMoving + 1 - alphaVoid).value()
331 forAll(movingPhases(), movingPhasei)
333 movingPhases()[movingPhasei] *= alphaVoid/sumAlphaMoving;
360 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
371 return nHatfv(alpha1, alpha2) & mesh_.Sf();
381 void Foam::multiphaseSystem::correctContactAngle
384 const phaseModel& phase2,
385 surfaceVectorField::Boundary& nHatb
388 const volScalarField::Boundary& gbf
389 = phase1.boundaryField();
391 const fvBoundaryMesh&
boundary = mesh_.boundary();
395 if (isA<alphaContactAngleFvPatchScalarField>(gbf[
patchi]))
397 const alphaContactAngleFvPatchScalarField& acap =
398 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
404 mesh_.Sf().boundaryField()[
patchi]
405 /mesh_.magSf().boundaryField()[
patchi]
408 alphaContactAngleFvPatchScalarField::thetaPropsTable::
411 .find(phasePairKey(phase1.name(), phase2.name()));
413 if (tp == acap.thetaProps().end())
416 <<
"Cannot find interface " 417 << phasePairKey(phase1.name(), phase2.name())
418 <<
"\n in table of theta properties for patch " 419 << acap.patch().name()
423 bool matched = (tp.key().first() == phase1.name());
425 scalar theta0 = convertToRad*tp().theta0(matched);
426 scalarField theta(boundary[patchi].size(), theta0);
428 scalar uTheta = tp().uTheta();
433 scalar thetaA = convertToRad*tp().thetaA(matched);
434 scalar thetaR = convertToRad*tp().thetaR(matched);
439 phase1.U()().boundaryField()[
patchi].patchInternalField()
440 - phase1.U()().boundaryField()[
patchi]
442 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
447 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
451 nWall /= (
mag(nWall) + small);
457 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
471 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
479 nHatPatch = a*AfHatPatch +
b*nHatPatch;
481 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
489 const phaseModel& phase1,
490 const phaseModel& phase2
493 tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
495 correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
498 return -
fvc::div(tnHatfv & mesh_.Sf());
525 cAlphas_(
lookup(
"interfaceCompression")),
536 mesh_.setFluxRequired(alphai.name());
551 const phaseModel& phase1
554 tmp<surfaceScalarField> tSurfaceTension
561 mesh_.time().timeName(),
568 dimensionSet(1, -2, -2, 0, 0),
576 const phaseModel& phase2 =
phases()[phasej];
578 if (&phase2 != &phase1)
580 phasePairKey key12(phase1.name(), phase2.name());
582 cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
584 if (cAlpha != cAlphas_.end())
586 tSurfaceTension.ref() +=
596 return tSurfaceTension;
603 tmp<volScalarField> tnearInt
610 mesh_.time().timeName(),
633 const Time&
runTime = mesh_.time();
642 tmp<volScalarField> trSubDeltaT;
650 PtrList<volScalarField> alpha0s(
phases().size());
651 PtrList<surfaceScalarField> alphaPhiSums(
phases().size());
671 "phiSum" + alpha.name(),
683 subCycleTime alphaSubCycle
685 const_cast<Time&>(runTime),
688 !(++alphaSubCycle).end();
702 if (phase.stationary())
continue;
710 alpha.timeIndex() = runTime.timeIndex();
713 alpha.oldTime() = alpha0s[
phasei];
714 alpha.oldTime().timeIndex() = runTime.timeIndex();
725 if (phase.stationary())
continue;
727 phase.alphaRhoPhiRef() =
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)
#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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
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.
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
GeometricField< scalar, fvPatchField, volMesh > volScalarField
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Calculate the first temporal derivative.
tmp< surfaceScalarField > surfaceTension(const phaseModel &phase) const
stressControl lookup("compactNormalStress") >> compactNormalStress
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
const phaseModel & phase() const
Return the phase.
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].
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Calculate the face-flux of the given field.
static word groupName(Name name, const word &group)
Calculate the matrix for the first temporal derivative.
const Type & value() const
Return const reference to value.
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)
Calculate the divergence of the given field.
dimensionedScalar pos0(const dimensionedScalar &ds)
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.
PtrList< surfaceScalarField > alphafs(phases.size())
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
virtual surfaceScalarField & alphaPhiRef()=0
Access the volumetric flux of the phase.
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)
surfaceScalarField phic(mixture.cAlpha() *mag(phi/mesh.magSf()))
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
MULES: Multidimensional universal limiter for explicit solution.
const doubleScalar e
Elementary charge.
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
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux(phi, alpha1, alphaScheme))
multiphaseSystem(const volVectorField &U, const surfaceScalarField &phi)
Construct from components.
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
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.