27 #include "surfaceTensionModel.H" 29 #include "surfaceInterpolate.H" 35 #include "alphaContactAngleFvPatchScalarField.H" 55 const phaseModelList& phaseModels
58 tmp<surfaceScalarField> tmpPhi
79 const dictTable& modelDicts
84 const phasePairKey& key = iter.key();
87 if (phasePairs_.found(key))
91 else if (key.ordered())
100 phaseModels_[key.first()],
101 phaseModels_[key.second()]
117 phaseModels_[key.first()],
118 phaseModels_[key.second()]
129 tmp<volScalarField> sumAlphaMoving
134 movingPhaseModels_[0],
135 calculatedFvPatchScalarField::typeName
141 label movingPhasei=1;
142 movingPhasei<movingPhaseModels_.size();
146 sumAlphaMoving.ref() += movingPhaseModels_[movingPhasei];
149 return sumAlphaMoving;
159 forAll(movingPhaseModels_, movingPhasei)
162 movingPhaseModels_[movingPhasei]
163 *movingPhaseModels_[movingPhasei].U();
166 forAll(movingPhaseModels_, movingPhasei)
168 movingPhaseModels_[movingPhasei].URef() += dUm;
175 const PtrList<surfaceScalarField>& alphafs,
183 forAll(movingPhaseModels_, movingPhasei)
186 alphafs[movingPhaseModels_[movingPhasei].index()]
187 *movingPhaseModels_[movingPhasei].phi();
190 forAll(movingPhaseModels_, movingPhasei)
192 movingPhaseModels_[movingPhasei].phiRef() += dphim;
219 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
230 return nHatfv(alpha1, alpha2) & mesh_.Sf();
236 const phaseModel& phase1,
237 const phaseModel& phase2,
238 surfaceVectorField::Boundary& nHatb
241 const volScalarField::Boundary& gbf
242 = phase1.boundaryField();
244 const fvBoundaryMesh&
boundary = mesh_.boundary();
248 if (isA<alphaContactAngleFvPatchScalarField>(gbf[
patchi]))
250 const alphaContactAngleFvPatchScalarField& acap =
251 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
257 mesh_.Sf().boundaryField()[
patchi]
258 /mesh_.magSf().boundaryField()[
patchi]
261 alphaContactAngleFvPatchScalarField::thetaPropsTable::
264 .find(phasePairKey(phase1.name(), phase2.name()));
266 if (tp == acap.thetaProps().end())
269 <<
"Cannot find interface " 270 << phasePairKey(phase1.name(), phase2.name())
271 <<
"\n in table of theta properties for patch " 272 << acap.patch().name()
276 bool matched = (tp.key().first() == phase1.name());
278 scalar theta0 =
degToRad(tp().theta0(matched));
279 scalarField theta(boundary[patchi].size(), theta0);
281 scalar uTheta = tp().uTheta();
286 scalar thetaA =
degToRad(tp().thetaA(matched));
287 scalar thetaR =
degToRad(tp().thetaR(matched));
292 phase1.U()().boundaryField()[
patchi].patchInternalField()
293 - phase1.U()().boundaryField()[
patchi]
295 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
300 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
304 nWall /= (
mag(nWall) + small);
310 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
324 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
332 nHatPatch = a*AfHatPatch +
b*nHatPatch;
334 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
342 const phaseModel& phase1,
343 const phaseModel& phase2
346 tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
348 correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
351 return -
fvc::div(tnHatfv & mesh_.Sf());
367 mesh.time().constant(),
369 IOobject::MUST_READ_IF_MODIFIED,
380 lookup<word>(
"type").find(
"TwoPhase") != string::npos
382 : lookupOrDefault(
"referencePhase", word::null)
388 phaseModel::iNew(*this, referencePhaseName_)
391 phi_(calcPhi(phaseModels_)),
407 cAlphas_(lookupOrDefault(
"interfaceCompression", cAlphaTable())),
416 label movingPhasei = 0;
417 label stationaryPhasei = 0;
418 label anisothermalPhasei = 0;
419 label multiComponentPhasei = 0;
422 phaseModel& phase = phaseModels_[
phasei];
423 movingPhasei += !phase.stationary();
424 stationaryPhasei += phase.stationary();
425 anisothermalPhasei += !phase.isothermal();
426 multiComponentPhasei += !phase.pure();
428 movingPhaseModels_.resize(movingPhasei);
429 stationaryPhaseModels_.resize(stationaryPhasei);
430 anisothermalPhaseModels_.resize(anisothermalPhasei);
431 multiComponentPhaseModels_.resize(multiComponentPhasei);
434 stationaryPhasei = 0;
435 anisothermalPhasei = 0;
436 multiComponentPhasei = 0;
439 phaseModel& phase = phaseModels_[
phasei];
440 if (!phase.stationary())
442 movingPhaseModels_.set(movingPhasei++, &phase);
444 if (phase.stationary())
446 stationaryPhaseModels_.set(stationaryPhasei++, &phase);
448 if (!phase.isothermal())
450 anisothermalPhaseModels_.set(anisothermalPhasei++, &phase);
454 multiComponentPhaseModels_.set(multiComponentPhasei++, &phase);
464 blendingMethods_.insert
477 generatePairsAndSubModels(
"surfaceTension", surfaceTensionModels_);
478 generatePairsAndSubModels(
"aspectRatio", aspectRatioModels_);
486 phaseModel* referencePhasePtr = &
phases()[referencePhaseName_];
493 if (&phaseModels_[
phasei] != referencePhasePtr)
495 referenceAlpha -= phaseModels_[
phasei];
503 mesh_.setFluxRequired(alphai.name());
518 tmp<volScalarField>
rho(movingPhaseModels_[0]*movingPhaseModels_[0].
rho());
522 label movingPhasei=1;
523 movingPhasei<movingPhaseModels_.size();
528 movingPhaseModels_[movingPhasei]
529 *movingPhaseModels_[movingPhasei].rho();
532 if (stationaryPhaseModels_.empty())
538 return rho/sumAlphaMoving();
545 tmp<volVectorField>
U(movingPhaseModels_[0]*movingPhaseModels_[0].
U());
549 label movingPhasei=1;
550 movingPhasei<movingPhaseModels_.size();
555 movingPhaseModels_[movingPhasei]
556 *movingPhaseModels_[movingPhasei].U();
559 if (stationaryPhaseModels_.empty())
565 return U/sumAlphaMoving();
573 if (aspectRatioModels_.found(key))
575 return aspectRatioModels_[key]->E();
581 aspectRatioModel::typeName +
":E",
592 if (surfaceTensionModels_.found(key))
594 return surfaceTensionModels_[key]->sigma();
600 surfaceTensionModel::typeName +
":sigma",
611 if (surfaceTensionModels_.found(key))
613 return surfaceTensionModels_[key]->sigma(patchi);
617 return tmp<scalarField>
628 tmp<volScalarField> tnearInt
658 phaseModels_[key.first()],
659 phaseModels_[key.second()]
673 return PtrList<volScalarField>(phaseModels_.size());
681 if (!phaseModels_[
phasei].incompressible())
705 const phaseModel& phase1
708 tmp<surfaceScalarField> tSurfaceTension
720 const phaseModel& phase2 =
phases()[phasej];
722 if (&phase2 != &phase1)
726 cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
728 if (cAlpha != cAlphas_.end())
730 tSurfaceTension.ref() +=
740 return tSurfaceTension;
748 phaseModels_[
phasei].correct();
755 const PtrList<volScalarField> dmdts = this->dmdts();
757 forAll(movingPhaseModels_, movingPhasei)
759 phaseModel& phase = movingPhaseModels_[movingPhasei];
773 if (
fvOptions().appliesToField(rho.name()))
778 if (dmdts.set(phase.index()))
780 source += dmdts[phase.index()];
783 phase.correctContinuityError(source);
790 bool updateDpdt =
false;
794 phaseModels_[
phasei].correctKinematics();
796 updateDpdt = updateDpdt || phaseModels_[
phasei].thermo().dpdt();
811 phaseModels_[
phasei].correctThermo();
820 phaseModels_[
phasei].correctReactions();
829 phaseModels_[
phasei].correctSpecies();
838 phaseModels_[
phasei].correctTurbulence();
847 phaseModels_[
phasei].correctEnergyTransport();
860 readOK &= phaseModels_[
phasei].read();
882 return vf/vf.
mesh().time().deltaT();
895 return sf/sf.
mesh().time().deltaT();
static const word propertiesName
Default name of the phase properties dictionary.
dimensionedScalar tanh(const dimensionedScalar &ds)
void generatePairs(const dictTable &modelDicts)
Generate pairs.
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.
tmp< volScalarField > E(const phasePairKey &key) const
Return the aspect-ratio for a pair.
virtual void correctTurbulence()
Correct the turbulence.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
static const surfaceScalarField & localRDeltaTf(const fvMesh &mesh)
Return the reciprocal of the local face time-step.
errorManipArg< error, int > exit(error &err, const int errNo=1)
tmp< surfaceScalarField > nHatf(const volScalarField &alpha1, const volScalarField &alpha2) const
Normal to interface between two phases dotted with face areas.
virtual bool read()
Read object.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Unit conversion functions.
const dimensionedScalar & sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
void setMixturePhi(const PtrList< surfaceScalarField > &alphafs, const surfaceScalarField &phim)
Re-normalise the flux of the phases.
static tmp< GeometricField< scalar, fvsPatchField, surfaceMesh > > New(const word &name, const Internal &, const PtrList< fvsPatchField< scalar >> &)
Return a temporary field constructed from name,.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
tmp< volScalarField > sigma(const phasePairKey &key) const
Return the surface tension coefficient for a pair.
tmp< volScalarField > K(const phaseModel &alpha1, const phaseModel &alpha2) const
Curvature of interface between two phases.
Calculate the snGrad of the given volField.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
rhoReactionThermo & thermo
GeometricField< vector, fvPatchField, volMesh > volVectorField
CGAL::Exact_predicates_exact_constructions_kernel K
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
tmp< surfaceScalarField > surfaceTension(const phaseModel &phase) const
Return the surface tension force.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Calculate the first temporal derivative.
virtual void correctContinuityError()
Correct the continuity errors.
virtual void correctKinematics()
Correct the kinematics.
stressControl lookup("compactNormalStress") >> compactNormalStress
static bool enabled(const fvMesh &mesh)
Return true if LTS is enabled.
dimensionedScalar cos(const dimensionedScalar &ds)
phaseSystem(const fvMesh &mesh)
Construct from fvMesh.
const dimensionedScalar & b
Wien displacement law constant: default SI units: [m K].
Calculate the gradient of the given field.
A class for handling words, derived from string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static word groupName(Name name, const word &group)
tmp< volVectorField > U() const
Return the mixture velocity.
static autoPtr< blendingMethod > New(const word &modelName, const dictionary &dict, const wordList &phaseNames)
virtual void correctThermo()
Correct the thermodynamics.
const Type & value() const
Return const reference to value.
static const word null
An empty word.
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)
const dimensionSet dimPressure
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
tmp< surfaceVectorField > nHatfv(const volScalarField &alpha1, const volScalarField &alpha2) const
Normal to interface between two phases.
Calculate the divergence of the given field.
dimensionedScalar pos0(const dimensionedScalar &ds)
virtual void correct()
Correct the fluid properties other than those listed below.
tmp< volScalarField > rho() const
Return the mixture density.
virtual PtrList< volScalarField > dmdts() const
Return the mass transfer rates for each phase.
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
tmp< volScalarField > sumAlphaMoving() const
Return the sum of the phase fractions of the moving phases.
const Mesh & mesh() const
Return mesh.
defineTypeNameAndDebug(combustionModel, 0)
volScalarField sf(fieldObject, mesh)
bool incompressible() const
Return incompressibility.
tmp< volScalarField > byDt(const volScalarField &vf)
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 dimensionSet dimDensity
virtual void correctReactions()
Correct the reactions.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
List< word > wordList
A List of words.
virtual void correctEnergyTransport()
Correct the energy transport e.g. alphat.
tmp< surfaceScalarField > calcPhi(const phaseModelList &phaseModels) const
Calculate and return the mixture flux.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
static const dimensionSet dimSigma
Surface tension coefficient dimensions.
virtual ~phaseSystem()
Destructor.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
virtual bool read()
Read base phaseProperties dictionary.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
phasePairKey()
Construct null.
phaseSystem::phaseModelList & phases
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
void setMixtureU(const volVectorField &Um)
Re-normalise the velocity of the phases.
const doubleScalar e
Elementary charge.
void correctContactAngle(const phaseModel &alpha1, const phaseModel &alpha2, surfaceVectorField::Boundary &nHatb) const
Correction for the boundary condition on the unit normal nHat on.
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
A class for managing temporary objects.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
virtual bool implicitPhasePressure() const
Returns true if the phase pressure is treated implicitly.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
virtual void correctSpecies()
Correct the species mass fractions.
virtual tmp< volScalarField > dmdtf(const phasePairKey &key) const
Return the mass transfer rate for an interface.