30 #include "alphaContactAngleFvPatchScalarField.H" 60 void Foam::multiphaseSystem::calcAlphas()
67 alphas_ += level*
phases()[i];
73 void Foam::multiphaseSystem::solveAlphas()
86 mesh_.time().timeName(),
92 forAll(stationaryPhases(), stationaryPhasei)
94 alphaVoid -= stationaryPhases()[stationaryPhasei];
108 upwind<scalar>(mesh_, phi_).interpolate(phase)
114 PtrList<surfaceScalarField> alphaPhiCorrs(
phases().size());
115 forAll(stationaryPhases(), stationaryPhasei)
117 phaseModel& phase = stationaryPhases()[stationaryPhasei];
125 - upwind<scalar>(mesh_, phi_).flux(phase)
129 forAll(movingPhases(), movingPhasei)
131 phaseModel& phase = movingPhases()[movingPhasei];
151 if (&
phase2 == &phase)
continue;
155 cAlphaTable::const_iterator cAlpha
157 cAlphas_.find(phasePairKey(phase.name(),
phase2.
name()))
160 if (cAlpha != cAlphas_.end())
183 phase.correctInflowOutflow(alphaPhiCorr);
193 min(alphaVoid.primitiveField(), phase.alphaMax())(),
201 forAll(stationaryPhases(), stationaryPhasei)
203 fixedAlphaPhiCorrs.insert(stationaryPhases()[stationaryPhasei].index());
208 forAll(movingPhases(), movingPhasei)
210 phaseModel& phase = movingPhases()[movingPhasei];
215 phase.correctInflowOutflow(
alphaPhi);
222 mesh_.time().timeName(),
235 if (phase.divU().valid())
241 if (dgdt[celli] > 0.0)
243 Sp[celli] -= dgdt[celli];
244 Su[celli] += dgdt[celli];
246 else if (dgdt[celli] < 0.0)
260 if (&
phase2 == &phase)
continue;
268 if (dgdt2[celli] < 0.0)
278 else if (dgdt2[celli] > 0.0)
280 Sp[celli] -= dgdt2[celli];
303 Info<< phase.name() <<
" fraction, min, max = " 304 << phase.weightedAverage(mesh_.V()).value()
305 <<
' ' <<
min(phase).value()
306 <<
' ' <<
max(phase).value()
315 mesh_.time().timeName(),
321 forAll(movingPhases(), movingPhasei)
323 sumAlphaMoving += movingPhases()[movingPhasei];
326 Info<<
"Phase-sum volume fraction, min, max = " 327 << (sumAlphaMoving + 1 - alphaVoid)().weightedAverage(mesh_.V()).value()
328 <<
' ' <<
min(sumAlphaMoving + 1 - alphaVoid).value()
329 <<
' ' <<
max(sumAlphaMoving + 1 - alphaVoid).value()
333 forAll(movingPhases(), movingPhasei)
335 movingPhases()[movingPhasei] *= alphaVoid/sumAlphaMoving;
362 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
383 void Foam::multiphaseSystem::correctContactAngle
393 const fvBoundaryMesh&
boundary = mesh_.boundary();
399 isA<reactingMultiphaseEuler::alphaContactAngleFvPatchScalarField>
408 const reactingMultiphaseEuler::
409 alphaContactAngleFvPatchScalarField
416 mesh_.Sf().boundaryField()[patchi]
417 /mesh_.magSf().boundaryField()[patchi]
420 reactingMultiphaseEuler::alphaContactAngleFvPatchScalarField::
421 thetaPropsTable::const_iterator tp =
425 if (tp == acap.thetaProps().end())
428 <<
"Cannot find interface " 430 <<
"\n in table of theta properties for patch " 431 << acap.patch().name()
435 bool matched = (tp.key().first() ==
phase1.
name());
437 scalar theta0 =
degToRad(tp().theta0(matched));
440 scalar uTheta = tp().uTheta();
445 const scalar thetaA =
degToRad(tp().thetaA(matched));
446 const scalar thetaR =
degToRad(tp().thetaR(matched));
451 phase1.
U()().boundaryField()[patchi].patchInternalField()
452 -
phase1.
U()().boundaryField()[patchi]
454 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
459 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
463 nWall /= (
mag(nWall) + SMALL);
469 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
483 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
491 nHatPatch = a*AfHatPatch +
b*nHatPatch;
493 nHatPatch /= (
mag(nHatPatch) + deltaN_.value());
505 tmp<surfaceVectorField> tnHatfv = nHatfv(
phase1,
phase2);
537 cAlphas_(
lookup(
"interfaceCompression")),
548 mesh_.setFluxRequired(alphai.
name());
568 tSurfaceTension.ref().setOriented();
578 cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
580 if (cAlpha != cAlphas_.end())
582 tSurfaceTension.ref() +=
592 tSurfaceTension->setOriented();
594 return tSurfaceTension;
633 tmp<volScalarField> trSubDeltaT;
641 List<volScalarField*> alphaPtrs(
phases().size());
642 PtrList<surfaceScalarField> alphaPhiSums(
phases().size());
670 subCycleTime alphaSubCycle
675 !(++alphaSubCycle).
end();
689 if (phase.stationary())
continue;
702 if (phase.stationary())
continue;
704 phase.alphaRhoPhiRef() =
707 phase.clamp_range(SMALL, 1 - SMALL);
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)
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &f1)
multiphaseSystem::phaseModelList & phases
virtual tmp< volScalarField > divU() const =0
Return the phase dilatation rate (d(alpha)/dt + div(alpha*phi))
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.
const word & name() const noexcept
Return the object name.
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
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.
static tmp< volScalarField > localRSubDeltaT(const fvMesh &mesh, const label nAlphaSubCycles)
Calculate and return the reciprocal of the local sub-cycling.
Type & refCast(U &obj)
A dynamic_cast (for references) to Type reference.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
Calculate the matrix for the laplacian of the field.
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.
DimensionedField< scalar, volMesh > Internal
The internal field type from which this GeometricField is derived.
CGAL::Exact_predicates_exact_constructions_kernel K
const surfaceScalarField & phi() const
dimensionedScalar sigma("sigma", dimMass/sqr(dimTime), transportProperties)
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Calculate the first temporal derivative.
const volVectorField & U() const
static word groupName(StringType base, const word &group)
Create dot-delimited name.group string.
const dimensionedScalar e
Elementary charge.
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()))
Class to represent a system of phases and model interfacial transfers between them.
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.
Calculate the matrix for the first temporal derivative.
Calculate the field for explicit evaluation of implicit and explicit sources.
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.
label find(const T &val) const
Find index of the first occurrence of the value.
static tmp< GeometricField< scalar, fvsPatchField, surfaceMesh > > New(const word &name, IOobjectOption::registerOption regOpt, const Mesh &mesh, const dimensionSet &dims, const word &patchFieldType=fvsPatchField< 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)
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.
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
defineTypeNameAndDebug(combustionModel, 0)
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.
const word & name() const noexcept
Return const reference to name.
surfaceScalarField phic(mixture.cAlpha() *mag(alphaPhic/mesh.magSf()))
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.
PtrList< surfaceScalarField > alphafs(phases.size())
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
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.
A class for managing temporary objects.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
const PtrDictionary< phaseModel > & phases() const
Return the phases.
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...
const Boundary & boundaryField() const noexcept
Return const-reference to the boundary field.
Do not request registration (bool: false)
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Calculate the finiteVolume matrix for implicit and explicit sources.
tmp< surfaceScalarField > surfaceTension(const phaseModel &phase) const
static constexpr const zero Zero
Global zero (0)
const volScalarField & alpha1