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()
79 PtrList<surfaceScalarField> alphaPhiCorrs(
phases().size());
90 "phi" + alpha1.name() +
"Corr",
95 "div(phi," + alpha1.name() +
')' 104 phaseModel& phase2 =
phases()[phasej];
107 if (&phase2 == &phase)
continue;
111 cAlphaTable::const_iterator cAlpha
113 cAlphas_.find(phasePairKey(phase.name(), phase2.name()))
116 if (cAlpha != cAlphas_.end())
120 (
mag(phi_) +
mag(phir))/mesh_.magSf()
128 "div(phir," + alpha2.name() +
',' + alpha1.name() +
')' 139 phase.correctInflowOutflow(alphaPhiCorr);
159 const scalar rDeltaT = 1.0/mesh_.time().deltaTValue();
184 mesh_.time().timeName(),
200 alphaPhi += upwind<scalar>(mesh_, phi_).
flux(phase);
201 phase.correctInflowOutflow(alphaPhi);
208 mesh_.time().timeName(),
220 mesh_.time().timeName(),
228 if (phase.divU().valid())
234 if (dgdt[celli] > 0.0)
236 Sp[celli] -= dgdt[celli];
237 Su[celli] += dgdt[celli];
239 else if (dgdt[celli] < 0.0)
243 *(1.0 - alpha[celli])/
max(alpha[celli], 1
e-4);
250 const phaseModel& phase2 =
phases()[phasej];
253 if (&phase2 == &phase)
continue;
255 if (phase2.divU().valid())
261 if (dgdt2[celli] < 0.0)
265 *(1.0 - alpha2[celli])/
max(alpha2[celli], 1
e-4);
269 *alpha[celli]/
max(alpha2[celli], 1
e-4);
271 else if (dgdt2[celli] > 0.0)
273 Sp[celli] -= dgdt2[celli];
290 Info<< phase.name() <<
" volume fraction, min, max = " 291 << phase.weightedAverage(mesh_.V()).value()
292 <<
' ' <<
min(phase).value()
293 <<
' ' <<
max(phase).value()
299 Info<<
"Phase-sum volume fraction, min, max = " 300 << sumAlpha.weightedAverage(mesh_.V()).value()
301 <<
' ' <<
min(sumAlpha).value()
302 <<
' ' <<
max(sumAlpha).value()
310 alpha += alpha*sumCorr;
337 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
348 return nHatfv(alpha1, alpha2) & mesh_.Sf();
358 void Foam::multiphaseSystem::correctContactAngle
361 const phaseModel& phase2,
362 surfaceVectorField::Boundary& nHatb
365 const volScalarField::Boundary& gbf
366 = phase1.boundaryField();
368 const fvBoundaryMesh&
boundary = mesh_.boundary();
372 if (isA<alphaContactAngleFvPatchScalarField>(gbf[
patchi]))
374 const alphaContactAngleFvPatchScalarField& acap =
375 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
381 mesh_.Sf().boundaryField()[
patchi]
382 /mesh_.magSf().boundaryField()[
patchi]
385 alphaContactAngleFvPatchScalarField::thetaPropsTable::
388 .find(phasePairKey(phase1.name(), phase2.name()));
390 if (tp == acap.thetaProps().end())
393 <<
"Cannot find interface " 394 << phasePairKey(phase1.name(), phase2.name())
395 <<
"\n in table of theta properties for patch " 396 << acap.patch().name()
400 bool matched = (tp.key().first() == phase1.name());
402 scalar theta0 = convertToRad*tp().theta0(matched);
403 scalarField theta(boundary[patchi].size(), theta0);
405 scalar uTheta = tp().uTheta();
410 scalar thetaA = convertToRad*tp().thetaA(matched);
411 scalar thetaR = convertToRad*tp().thetaR(matched);
416 phase1.U()().boundaryField()[
patchi].patchInternalField()
417 - phase1.U()().boundaryField()[
patchi]
419 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
424 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
428 nWall /= (
mag(nWall) + SMALL);
434 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
448 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
456 nHatPatch = a*AfHatPatch +
b*nHatPatch;
458 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
466 const phaseModel& phase1,
467 const phaseModel& phase2
470 tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
472 correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
475 return -
fvc::div(tnHatfv & mesh_.Sf());
502 cAlphas_(
lookup(
"interfaceCompression")),
513 mesh_.setFluxRequired(alphai.name());
528 const phaseModel& phase1
531 tmp<surfaceScalarField> tSurfaceTension
538 mesh_.time().timeName(),
545 dimensionSet(1, -2, -2, 0, 0),
553 const phaseModel& phase2 =
phases()[phasej];
555 if (&phase2 != &phase1)
557 phasePairKey key12(phase1.name(), phase2.name());
559 cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
561 if (cAlpha != cAlphas_.end())
563 tSurfaceTension.ref() +=
573 return tSurfaceTension;
580 tmp<volScalarField> tnearInt
587 mesh_.time().timeName(),
610 const Time& runTime = mesh_.time();
619 tmp<volScalarField> trSubDeltaT;
627 PtrList<volScalarField> alpha0s(
phases().size());
628 PtrList<surfaceScalarField> alphaPhiSums(
phases().size());
648 "phiSum" + alpha.name(),
660 subCycleTime alphaSubCycle
662 const_cast<Time&>(runTime),
665 !(++alphaSubCycle).end();
685 alpha.timeIndex() = runTime.timeIndex();
688 alpha.oldTime() = alpha0s[
phasei];
689 alpha.oldTime().timeIndex() = runTime.timeIndex();
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)
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.
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
const surfaceScalarField & alphaPhi() const
GeometricField< scalar, fvPatchField, volMesh > volScalarField
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)
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)
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.
Calulate the matrix for the first temporal derivative.
const Type & value() const
Return const reference to value.
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)
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.
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)
surfaceScalarField phic(mixture.cAlpha() *mag(phi/mesh.magSf()))
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
MULES: Multidimensional universal limiter for explicit solution.
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.