27 #include "surfaceTensionModel.H" 28 #include "surfaceInterpolate.H" 36 #include "alphaContactAngleFvPatchScalarField.H" 58 const phaseModelList& phaseModels
61 tmp<surfaceScalarField> tmpPhi
82 tmp<volScalarField> sumAlphaMoving
87 movingPhaseModels_[0],
88 calculatedFvPatchScalarField::typeName
95 movingPhasei<movingPhaseModels_.size();
99 sumAlphaMoving.ref() += movingPhaseModels_[movingPhasei];
102 return sumAlphaMoving;
112 forAll(movingPhaseModels_, movingPhasei)
115 movingPhaseModels_[movingPhasei]
116 *movingPhaseModels_[movingPhasei].U();
119 forAll(movingPhaseModels_, movingPhasei)
121 movingPhaseModels_[movingPhasei].URef() += dUm;
128 const PtrList<surfaceScalarField>& alphafs,
136 forAll(movingPhaseModels_, movingPhasei)
139 alphafs[movingPhaseModels_[movingPhasei].index()]
140 *movingPhaseModels_[movingPhasei].phi();
143 forAll(movingPhaseModels_, movingPhasei)
145 movingPhaseModels_[movingPhasei].phiRef() += dphim;
172 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
183 return nHatfv(alpha1, alpha2) & mesh_.Sf();
189 const phaseModel& phase1,
190 const phaseModel& phase2
193 tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
199 phase1.U()().boundaryField(),
201 tnHatfv.ref().boundaryFieldRef()
205 return -
fvc::div(tnHatfv & mesh_.Sf());
221 mesh.time().constant(),
223 IOobject::MUST_READ_IF_MODIFIED,
230 referencePhaseName_(lookupOrDefault(
"referencePhase", word::null)),
235 phaseModel::iNew(*this, referencePhaseName_)
238 phi_(
"phi", calcPhi(phaseModels_)),
261 label movingPhasei = 0;
262 label stationaryPhasei = 0;
263 label anisothermalPhasei = 0;
264 label multiComponentPhasei = 0;
267 phaseModel& phase = phaseModels_[
phasei];
268 movingPhasei += !phase.stationary();
269 stationaryPhasei += phase.stationary();
270 anisothermalPhasei += !phase.isothermal();
271 multiComponentPhasei += !phase.pure();
273 movingPhaseModels_.resize(movingPhasei);
274 stationaryPhaseModels_.resize(stationaryPhasei);
275 anisothermalPhaseModels_.resize(anisothermalPhasei);
276 multiComponentPhaseModels_.resize(multiComponentPhasei);
279 stationaryPhasei = 0;
280 anisothermalPhasei = 0;
281 multiComponentPhasei = 0;
284 phaseModel& phase = phaseModels_[
phasei];
285 if (!phase.stationary())
287 movingPhaseModels_.set(movingPhasei++, &phase);
289 if (phase.stationary())
291 stationaryPhaseModels_.set(stationaryPhasei++, &phase);
293 if (!phase.isothermal())
295 anisothermalPhaseModels_.set(anisothermalPhasei++, &phase);
299 multiComponentPhaseModels_.set(multiComponentPhasei++, &phase);
307 if (this->
found(
"interfaceCompression"))
309 generateInterfacialValues(
"interfaceCompression", cAlphas_);
313 generateInterfacialModels(surfaceTensionModels_);
321 phaseModel* referencePhasePtr = &
phases()[referencePhaseName_];
328 if (&phaseModels_[
phasei] != referencePhasePtr)
330 referenceAlpha -= phaseModels_[
phasei];
338 mesh_.schemes().setFluxRequired(alphai.name());
353 tmp<volScalarField>
rho(movingPhaseModels_[0]*movingPhaseModels_[0].
rho());
357 label movingPhasei=1;
358 movingPhasei<movingPhaseModels_.size();
363 movingPhaseModels_[movingPhasei]
364 *movingPhaseModels_[movingPhasei].rho();
367 if (stationaryPhaseModels_.empty())
373 return rho/sumAlphaMoving();
380 tmp<volVectorField>
U(movingPhaseModels_[0]*movingPhaseModels_[0].
U());
384 label movingPhasei=1;
385 movingPhasei<movingPhaseModels_.size();
390 movingPhaseModels_[movingPhasei]
391 *movingPhaseModels_[movingPhasei].U();
394 if (stationaryPhaseModels_.empty())
400 return U/sumAlphaMoving();
408 if (surfaceTensionModels_.found(key))
410 return surfaceTensionModels_[key]->sigma();
416 surfaceTensionModel::typeName +
":sigma",
427 if (surfaceTensionModels_.found(key))
429 return surfaceTensionModels_[key]->sigma(patchi);
433 return tmp<scalarField>
444 tmp<volScalarField> tnearInt
469 const phaseInterfaceKey& key
483 return PtrList<volScalarField>(phaseModels_.size());
489 return PtrList<volScalarField>(phaseModels_.size());
521 const phaseModel& phase1
524 tmp<surfaceScalarField> tSurfaceTension
536 const phaseModel& phase2 =
phases()[phasej];
538 if (&phase2 != &phase1)
540 const phaseInterface interface(phase1, phase2);
542 if (cAlphas_.found(interface))
544 tSurfaceTension.ref() +=
554 return tSurfaceTension;
562 phaseModels_[
phasei].correct();
569 const PtrList<volScalarField> dmdts = this->dmdts();
571 forAll(movingPhaseModels_, movingPhasei)
573 phaseModel& phase = movingPhaseModels_[movingPhasei];
587 if (
fvModels().addsSupToField(rho.name()))
592 if (dmdts.set(phase.index()))
594 source += dmdts[phase.index()];
597 phase.correctContinuityError(source);
604 bool updateDpdt =
false;
608 phaseModels_[
phasei].correctKinematics();
610 updateDpdt = updateDpdt || phaseModels_[
phasei].thermo().dpdt();
625 phaseModels_[
phasei].correctThermo();
634 phaseModels_[
phasei].correctReactions();
643 phaseModels_[
phasei].correctSpecies();
652 phaseModels_[
phasei].correctTurbulence();
661 phaseModels_[
phasei].correctEnergyTransport();
668 if (mesh_.changing())
682 forAll(movingPhases(), movingPhasei)
684 phaseModel& phase = movingPhases()[movingPhasei];
688 FieldField<fvsPatchField, scalar> phiRelBf
690 MRF_.relative(mesh_.Sf().boundaryField() & UBf)
699 isA<fixedValueFvsPatchScalarField>(phiBf[patchi])
700 && !isA<movingWallVelocityFvPatchVectorField>(UBf[
patchi])
713 const tmp<volScalarField>&
divU,
715 nonOrthogonalSolutionControl&
pimple 718 forAll(movingPhases(), movingPhasei)
720 phaseModel& phase = movingPhases()[movingPhasei];
727 if (Ubf[patchi].fixesValue())
729 Ubf[
patchi].initEvaluate();
735 if (Ubf[patchi].fixesValue())
744 correctBoundaryFlux();
748 PtrList<surfaceScalarField>
alphafs(phaseModels_.size());
749 forAll(movingPhases(), movingPhasei)
751 phaseModel& phase = movingPhases()[movingPhasei];
759 phi_ += alphafs[
phasei]*(mesh_.Sf() & phase.Uf());
765 movingPhases()[0].
U(),
777 setMixturePhi(alphafs, phi_);
788 forAll(phaseModels_, phasei)
790 readOK &= phaseModels_[
phasei].read();
812 return vf/vf.
mesh().time().deltaT();
825 return sf/sf.
mesh().time().deltaT();
static const word propertiesName
Default name of the phase properties dictionary.
pressureReference & pressureReference
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
#define forAll(list, i)
Loop across all elements in list.
virtual void correctTurbulence()
Correct the turbulence.
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
pimpleNoLoopControl & pimple
static const surfaceScalarField & localRDeltaTf(const fvMesh &mesh)
Return the reciprocal of the local face time-step.
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.
PtrList< surfaceScalarField > alphafs(phases.size())
void correctContactAngle(const volScalarField &alpha1, const volScalarField &alpha2, const volVectorField::Boundary &Ubf, const dimensionedScalar &deltaN, surfaceVectorField::Boundary &nHatbf)
Correct the contact angle for the two volume fraction fields.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
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
tmp< volScalarField > K(const phaseModel &alpha1, const phaseModel &alpha2) const
Curvature of interface between two phases.
virtual bool incompressible() const
Return whether the phase is incompressible.
Calculate the snGrad of the given volField.
void correctPhi(const volScalarField &p_rgh, const tmp< volScalarField > &divU, const pressureReference &pressureReference, nonOrthogonalSolutionControl &pimple)
GeometricBoundaryField< vector, fvPatchField, volMesh > Boundary
Type of the boundary field.
const dimensionSet dimless
GeometricField< vector, fvPatchField, volMesh > volVectorField
CGAL::Exact_predicates_exact_constructions_kernel K
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
virtual volScalarField & p()=0
Pressure [Pa].
virtual const rhoThermo & thermo() const
Return the thermophysical model.
void correctBoundaryFlux()
Correct fixed-flux BCs to be consistent with the velocity BCs.
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.
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)
virtual tmp< volScalarField > dmdtf(const phaseInterfaceKey &key) const
Return the mass transfer rate for an interface.
tmp< volVectorField > U() const
Return the mixture velocity.
const dimensionSet dimDensity
virtual void correctThermo()
Correct the thermodynamics.
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)
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.
word name(const complex &)
Return a string representation of a complex.
tmp< volScalarField > byDt(const volScalarField &vf)
tmp< volScalarField > sigma(const phaseInterfaceKey &key) const
Return the surface tension coefficient for an interface.
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 tmp< volScalarField > rho() const
Return the density field.
virtual void correctEnergyTransport()
Correct the energy transport e.g. alphat.
tmp< surfaceScalarField > calcPhi(const phaseModelList &phaseModels) const
Calculate and return the mixture flux.
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.
phaseSystem::phaseModelList & phases
dimensioned< scalar > mag(const dimensioned< Type > &)
void setMixtureU(const volVectorField &Um)
Re-normalise the velocity of the phases.
const doubleScalar e
Elementary charge.
void CorrectPhi(surfaceScalarField &phi, const volVectorField &U, const volScalarField &p, const RAUfType &rAUf, const DivUType &divU, const pressureReference &pressureReference, nonOrthogonalSolutionControl &pcorrControl)
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.