27 #include "surfaceTensionModel.H" 29 #include "surfaceInterpolate.H" 37 #include "alphaContactAngleFvPatchScalarField.H" 60 const phaseModelList& phaseModels
63 tmp<surfaceScalarField> tmpPhi
84 const dictTable& modelDicts
89 const phasePairKey& key = iter.key();
92 if (phasePairs_.found(key))
96 else if (key.ordered())
105 phaseModels_[key.first()],
106 phaseModels_[key.second()]
122 phaseModels_[key.first()],
123 phaseModels_[key.second()]
134 tmp<volScalarField> sumAlphaMoving
139 movingPhaseModels_[0],
140 calculatedFvPatchScalarField::typeName
146 label movingPhasei=1;
147 movingPhasei<movingPhaseModels_.size();
151 sumAlphaMoving.ref() += movingPhaseModels_[movingPhasei];
154 return sumAlphaMoving;
164 forAll(movingPhaseModels_, movingPhasei)
167 movingPhaseModels_[movingPhasei]
168 *movingPhaseModels_[movingPhasei].U();
171 forAll(movingPhaseModels_, movingPhasei)
173 movingPhaseModels_[movingPhasei].URef() += dUm;
180 const PtrList<surfaceScalarField>& alphafs,
188 forAll(movingPhaseModels_, movingPhasei)
191 alphafs[movingPhaseModels_[movingPhasei].index()]
192 *movingPhaseModels_[movingPhasei].phi();
195 forAll(movingPhaseModels_, movingPhasei)
197 movingPhaseModels_[movingPhasei].phiRef() += dphim;
224 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
235 return nHatfv(alpha1, alpha2) & mesh_.Sf();
241 const phaseModel& phase1,
242 const phaseModel& phase2,
243 surfaceVectorField::Boundary& nHatb
246 const volScalarField::Boundary& gbf
247 = phase1.boundaryField();
249 const fvBoundaryMesh&
boundary = mesh_.boundary();
253 if (isA<alphaContactAngleFvPatchScalarField>(gbf[
patchi]))
255 const alphaContactAngleFvPatchScalarField& acap =
256 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
262 mesh_.Sf().boundaryField()[
patchi]
263 /mesh_.magSf().boundaryField()[
patchi]
266 alphaContactAngleFvPatchScalarField::thetaPropsTable::
269 .find(phasePairKey(phase1.name(), phase2.name()));
271 if (tp == acap.thetaProps().end())
274 <<
"Cannot find interface " 275 << phasePairKey(phase1.name(), phase2.name())
276 <<
"\n in table of theta properties for patch " 277 << acap.patch().name()
281 bool matched = (tp.key().first() == phase1.name());
283 scalar theta0 =
degToRad(tp().theta0(matched));
284 scalarField theta(boundary[patchi].size(), theta0);
286 scalar uTheta = tp().uTheta();
291 scalar thetaA =
degToRad(tp().thetaA(matched));
292 scalar thetaR =
degToRad(tp().thetaR(matched));
297 phase1.U()().boundaryField()[
patchi].patchInternalField()
298 - phase1.U()().boundaryField()[
patchi]
300 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
305 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
309 nWall /= (
mag(nWall) + small);
315 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
329 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
337 nHatPatch = a*AfHatPatch +
b*nHatPatch;
339 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
347 const phaseModel& phase1,
348 const phaseModel& phase2
351 tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
353 correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
356 return -
fvc::div(tnHatfv & mesh_.Sf());
372 mesh.time().constant(),
374 IOobject::MUST_READ_IF_MODIFIED,
381 referencePhaseName_(lookupOrDefault(
"referencePhase", word::null)),
386 phaseModel::iNew(*this, referencePhaseName_)
389 phi_(calcPhi(phaseModels_)),
405 cAlphas_(lookupOrDefault(
"interfaceCompression", cAlphaTable())),
414 label movingPhasei = 0;
415 label stationaryPhasei = 0;
416 label anisothermalPhasei = 0;
417 label multiComponentPhasei = 0;
420 phaseModel& phase = phaseModels_[
phasei];
421 movingPhasei += !phase.stationary();
422 stationaryPhasei += phase.stationary();
423 anisothermalPhasei += !phase.isothermal();
424 multiComponentPhasei += !phase.pure();
426 movingPhaseModels_.resize(movingPhasei);
427 stationaryPhaseModels_.resize(stationaryPhasei);
428 anisothermalPhaseModels_.resize(anisothermalPhasei);
429 multiComponentPhaseModels_.resize(multiComponentPhasei);
432 stationaryPhasei = 0;
433 anisothermalPhasei = 0;
434 multiComponentPhasei = 0;
437 phaseModel& phase = phaseModels_[
phasei];
438 if (!phase.stationary())
440 movingPhaseModels_.set(movingPhasei++, &phase);
442 if (phase.stationary())
444 stationaryPhaseModels_.set(stationaryPhasei++, &phase);
446 if (!phase.isothermal())
448 anisothermalPhaseModels_.set(anisothermalPhasei++, &phase);
452 multiComponentPhaseModels_.set(multiComponentPhasei++, &phase);
462 blendingMethods_.insert
475 generatePairsAndSubModels(
"surfaceTension", surfaceTensionModels_);
476 generatePairsAndSubModels(
"aspectRatio", aspectRatioModels_);
484 phaseModel* referencePhasePtr = &
phases()[referencePhaseName_];
491 if (&phaseModels_[
phasei] != referencePhasePtr)
493 referenceAlpha -= phaseModels_[
phasei];
501 mesh_.setFluxRequired(alphai.name());
516 tmp<volScalarField>
rho(movingPhaseModels_[0]*movingPhaseModels_[0].
rho());
520 label movingPhasei=1;
521 movingPhasei<movingPhaseModels_.size();
526 movingPhaseModels_[movingPhasei]
527 *movingPhaseModels_[movingPhasei].rho();
530 if (stationaryPhaseModels_.empty())
536 return rho/sumAlphaMoving();
543 tmp<volVectorField>
U(movingPhaseModels_[0]*movingPhaseModels_[0].
U());
547 label movingPhasei=1;
548 movingPhasei<movingPhaseModels_.size();
553 movingPhaseModels_[movingPhasei]
554 *movingPhaseModels_[movingPhasei].U();
557 if (stationaryPhaseModels_.empty())
563 return U/sumAlphaMoving();
571 if (aspectRatioModels_.found(key))
573 return aspectRatioModels_[key]->E();
579 aspectRatioModel::typeName +
":E",
590 if (surfaceTensionModels_.found(key))
592 return surfaceTensionModels_[key]->sigma();
598 surfaceTensionModel::typeName +
":sigma",
609 if (surfaceTensionModels_.found(key))
611 return surfaceTensionModels_[key]->sigma(patchi);
615 return tmp<scalarField>
626 tmp<volScalarField> tnearInt
656 phaseModels_[key.first()],
657 phaseModels_[key.second()]
671 return PtrList<volScalarField>(phaseModels_.size());
677 return PtrList<volScalarField>(phaseModels_.size());
685 if (!phaseModels_[
phasei].incompressible())
709 const phaseModel& phase1
712 tmp<surfaceScalarField> tSurfaceTension
724 const phaseModel& phase2 =
phases()[phasej];
726 if (&phase2 != &phase1)
730 cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
732 if (cAlpha != cAlphas_.end())
734 tSurfaceTension.ref() +=
744 return tSurfaceTension;
752 phaseModels_[
phasei].correct();
759 const PtrList<volScalarField> dmdts = this->dmdts();
761 forAll(movingPhaseModels_, movingPhasei)
763 phaseModel& phase = movingPhaseModels_[movingPhasei];
777 if (
fvModels().addsSupToField(rho.name()))
782 if (dmdts.set(phase.index()))
784 source += dmdts[phase.index()];
787 phase.correctContinuityError(
source);
794 bool updateDpdt =
false;
798 phaseModels_[
phasei].correctKinematics();
800 updateDpdt = updateDpdt || phaseModels_[
phasei].thermo().dpdt();
815 phaseModels_[
phasei].correctThermo();
824 phaseModels_[
phasei].correctReactions();
833 phaseModels_[
phasei].correctSpecies();
842 phaseModels_[
phasei].correctTurbulence();
851 phaseModels_[
phasei].correctEnergyTransport();
858 if (mesh_.changing())
872 forAll(movingPhases(), movingPhasei)
874 phaseModel& phase = movingPhases()[movingPhasei];
876 const volVectorField::Boundary& UBf = phase.U()().boundaryField();
878 FieldField<fvsPatchField, scalar> phiRelBf
880 MRF_.relative(mesh_.Sf().boundaryField() & UBf)
883 surfaceScalarField::Boundary& phiBf = phase.phiRef().boundaryFieldRef();
889 isA<fixedValueFvsPatchScalarField>(phiBf[patchi])
890 && !isA<movingWallVelocityFvPatchVectorField>(UBf[
patchi])
903 const tmp<volScalarField>&
divU,
905 nonOrthogonalSolutionControl&
pimple 908 forAll(movingPhases(), movingPhasei)
910 phaseModel& phase = movingPhases()[movingPhasei];
912 volVectorField::Boundary& Ubf = phase.URef().boundaryFieldRef();
913 surfaceVectorField::Boundary& UfBf = phase.UfRef().boundaryFieldRef();
917 if (Ubf[patchi].fixesValue())
919 Ubf[
patchi].initEvaluate();
925 if (Ubf[patchi].fixesValue())
934 correctBoundaryFlux();
938 PtrList<surfaceScalarField>
alphafs(phaseModels_.size());
939 forAll(movingPhases(), movingPhasei)
941 phaseModel& phase = movingPhases()[movingPhasei];
949 phi_ += alphafs[
phasei]*(mesh_.Sf() & phase.Uf());
955 movingPhases()[0].
U(),
967 setMixturePhi(alphafs, phi_);
978 forAll(phaseModels_, phasei)
980 readOK &= phaseModels_[
phasei].read();
1002 return vf/vf.
mesh().time().deltaT();
1015 return sf/sf.
mesh().time().deltaT();
static const word propertiesName
Default name of the phase properties dictionary.
pressureReference & pressureReference
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.
fluidReactionThermo & thermo
pimpleNoLoopControl & pimple
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 PtrList< volScalarField > d2mdtdps() const
Return the mass transfer pressure implicit coefficients.
virtual bool read()
Read object.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
PtrList< surfaceScalarField > alphafs(phases.size())
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Unit conversion functions.
const dimensionSet dimPressure
void setMixturePhi(const PtrList< surfaceScalarField > &alphafs, const surfaceScalarField &phim)
Re-normalise the flux of the phases.
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
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
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
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.
void correctPhi(const volScalarField &p_rgh, const tmp< volScalarField > &divU, const pressureReference &pressureReference, nonOrthogonalSolutionControl &pimple)
const dimensionSet dimless
dimensionedScalar det(const dimensionedSphericalTensor &dt)
GeometricField< vector, fvPatchField, volMesh > volVectorField
CGAL::Exact_predicates_exact_constructions_kernel K
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
void correctBoundaryFlux()
Correct fixed-flux BCs to be consistent with the velocity BCs.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
void CorrectPhi(surfaceScalarField &phi, const volVectorField &U, const volScalarField &p, const RAUfType &rAUf, const DivUType &divU, const pressureReference &pressureReference, nonOrthogonalSolutionControl &pcorrControl)
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
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
Calculate the first temporal derivative.
virtual void correctContinuityError()
Correct the continuity errors.
const dimensionSet dimTime
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.
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)
const dimensionSet dimDensity
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)
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)
tmp< volScalarField > divU
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.
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
virtual void correctReactions()
Correct the reactions.
Foam::fvModels & fvModels
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
virtual void correctEnergyTransport()
Correct the energy transport e.g. alphat.
tmp< surfaceScalarField > calcPhi(const phaseModelList &phaseModels) const
Calculate and return the mixture flux.
fvModels source(alpha1, mixture.thermo1().rho())
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.
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.
A class for managing temporary objects.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
void makeRelative(surfaceScalarField &phi, const volVectorField &U)
Make the given flux relative.
virtual void meshUpdate()
Update the fluid properties for mesh changes.
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.