30#include "alphaContactAngleFvPatchScalarField.H"
39#include "surfaceInterpolate.H"
52void Foam::multiphaseMixtureThermo::calcAlphas()
57 for (
const phaseModel& phase : phases_)
59 alphas_ += level * phase;
73 psiThermo(
U.
mesh(), word::null),
74 phases_(lookup(
"phases"), phaseModel::iNew(p_, T_)),
108 sigmas_(lookup(
"sigmas")),
109 dimSigma_(1, 0, -2, 0, 0),
116 rhoPhi_.setOriented();
127 for (phaseModel& phase : phases_)
132 auto phasei = phases_.cbegin();
149 for (phaseModel& phase : phases_)
151 phase.thermo().rho() += phase.thermo().psi()*dp;
158 auto phasei = phases_.cbegin();
173 for (
const phaseModel& phase : phases_)
175 if (!phase.thermo().incompressible())
187 for (
const phaseModel& phase : phases_)
189 if (!phase.thermo().isochoric())
205 auto phasei = phases_.cbegin();
225 auto phasei = phases_.cbegin();
249 auto phasei = phases_.cbegin();
259 phasei().boundaryField()[patchi]*
phasei().thermo().he(
p,
T, patchi);
268 auto phasei = phases_.cbegin();
309 auto phasei = phases_.cbegin();
327 auto phasei = phases_.cbegin();
329 tmp<scalarField>
trho
337 phasei().boundaryField()[patchi]*
phasei().thermo().rho(patchi);
346 auto phasei = phases_.cbegin();
366 auto phasei = phases_.cbegin();
376 phasei().boundaryField()[patchi]*
phasei().thermo().Cp(
p,
T, patchi);
385 auto phasei = phases_.cbegin();
405 auto phasei = phases_.cbegin();
415 phasei().boundaryField()[patchi]*
phasei().thermo().Cv(
p,
T, patchi);
424 auto phasei = phases_.cbegin();
444 auto phasei = phases_.cbegin();
446 tmp<scalarField> tgamma
454 phasei().boundaryField()[patchi]
455 *
phasei().thermo().gamma(
p,
T, patchi);
464 auto phasei = phases_.cbegin();
484 auto phasei = phases_.cbegin();
486 tmp<scalarField> tCpv
494 phasei().boundaryField()[patchi]
495 *
phasei().thermo().Cpv(
p,
T, patchi);
504 auto phasei = phases_.cbegin();
524 auto phasei = phases_.cbegin();
526 tmp<scalarField> tCpByCpv
534 phasei().boundaryField()[patchi]
535 *
phasei().thermo().CpByCpv(
p,
T, patchi);
544 auto phasei = phases_.cbegin();
568 return mu(patchi)/
rho(patchi);
574 auto phasei = phases_.cbegin();
592 auto phasei = phases_.cbegin();
594 tmp<scalarField> tkappa
602 phasei().boundaryField()[patchi]*
phasei().thermo().kappa(patchi);
611 PtrDictionary<phaseModel>::const_iterator
phasei = phases_.begin();
617 talphaEff.ref() +=
phasei()*
phasei().thermo().alphahe();
629 PtrDictionary<phaseModel>::const_iterator
phasei = phases_.begin();
631 tmp<scalarField> talphaEff
633 phasei().boundaryField()[patchi]
640 phasei().boundaryField()[patchi]
641 *
phasei().thermo().alphahe(patchi);
653 auto phasei = phases_.cbegin();
659 tkappaEff.ref() +=
phasei()*
phasei().thermo().kappaEff(alphat);
672 auto phasei = phases_.cbegin();
674 tmp<scalarField> tkappaEff
676 phasei().boundaryField()[patchi]
683 phasei().boundaryField()[patchi]
684 *
phasei().thermo().kappaEff(alphat, patchi);
696 auto phasei = phases_.cbegin();
702 talphaEff.ref() +=
phasei()*
phasei().thermo().alphaEff(alphat);
715 auto phasei = phases_.cbegin();
717 tmp<scalarField> talphaEff
719 phasei().boundaryField()[patchi]
726 phasei().boundaryField()[patchi]
727 *
phasei().thermo().alphaEff(alphat, patchi);
736 auto phasei = phases_.cbegin();
752 tmp<surfaceScalarField> tstf
758 "surfaceTensionForce",
759 mesh_.time().timeName(),
785 <<
"Cannot find interface " << interfacePair(
alpha1,
alpha2)
786 <<
" in list of sigma values"
821 !(++alphaSubCycle).
end();
859 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
880void Foam::multiphaseMixtureThermo::correctContactAngle
890 const fvBoundaryMesh&
boundary = mesh_.boundary();
894 if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
896 const alphaContactAngleFvPatchScalarField& acap =
897 refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
903 mesh_.Sf().boundaryField()[patchi]
904 /mesh_.magSf().boundaryField()[patchi]
913 <<
"Cannot find interface " << interfacePair(
alpha1,
alpha2)
914 <<
"\n in table of theta properties for patch "
915 << acap.patch().name()
919 const bool matched = (tp.key().first() ==
alpha1.
name());
921 const scalar theta0 =
degToRad(tp().theta0(matched));
924 const scalar uTheta = tp().uTheta();
929 const scalar thetaA =
degToRad(tp().thetaA(matched));
930 const scalar thetaR =
degToRad(tp().thetaR(matched));
935 U_.boundaryField()[patchi].patchInternalField()
936 - U_.boundaryField()[patchi]
938 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
943 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
947 nWall /= (
mag(nWall) + SMALL);
953 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
967 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
975 nHatPatch = a*AfHatPatch +
b*nHatPatch;
977 nHatPatch /= (
mag(nHatPatch) + deltaN_.value());
989 tmp<surfaceVectorField> tnHatfv = nHatfv(
alpha1,
alpha2);
994 return -
fvc::div(tnHatfv & mesh_.Sf());
1001 tmp<volScalarField> tnearInt
1008 mesh_.time().timeName(),
1016 for (
const phaseModel& phase : phases_)
1019 max(tnearInt(),
pos0(phase - 0.01)*
pos0(0.99 - phase));
1026void Foam::multiphaseMixtureThermo::solveAlphas
1031 static label nSolves(-1);
1034 const word alphaScheme(
"div(phi,alpha)");
1040 PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
1043 for (phaseModel&
alpha : phases_)
1062 for (phaseModel&
alpha2 : phases_)
1078 1.0/mesh_.time().deltaT().value(),
1079 geometricOneField(),
1102 mesh_.time().timeName(),
1115 for (phaseModel&
alpha : phases_)
1125 mesh_.time().timeName(),
1137 mesh_.time().timeName(),
1150 if (dgdt[celli] < 0.0 &&
alpha[celli] > 0.0)
1152 Sp[celli] += dgdt[celli]*
alpha[celli];
1153 Su[celli] -= dgdt[celli]*
alpha[celli];
1155 else if (dgdt[celli] > 0.0 &&
alpha[celli] < 1.0)
1157 Sp[celli] -= dgdt[celli]*(1.0 -
alpha[celli]);
1162 for (
const phaseModel&
alpha2 : phases_)
1170 if (dgdt2[celli] > 0.0 &&
alpha2[celli] < 1.0)
1172 Sp[celli] -= dgdt2[celli]*(1.0 -
alpha2[celli]);
1173 Su[celli] += dgdt2[celli]*
alpha[celli];
1175 else if (dgdt2[celli] < 0.0 &&
alpha2[celli] > 0.0)
1177 Sp[celli] += dgdt2[celli]*
alpha2[celli];
1184 geometricOneField(),
1204 Info<<
"Phase-sum volume fraction, min, max = "
1205 << sumAlpha.weightedAverage(mesh_.V()).value()
1206 <<
' ' <<
min(sumAlpha).value()
1207 <<
' ' <<
max(sumAlpha).value()
CGAL::Exact_predicates_exact_constructions_kernel K
surfaceScalarField phic(mixture.cAlpha() *mag(alphaPhic/mesh.magSf()))
MULES: Multidimensional universal limiter for explicit solution.
const volScalarField & alpha1
const volScalarField & alpha2
const dimensionSet & dimensions() const
Return dimensions.
dimensioned< Type > weightedAverage(const DimensionedField< scalar, GeoMesh > &weightField) const
Calculate and return weighted average.
void setOriented(const bool oriented=true) noexcept
Set the oriented flag.
DimensionedField< scalar, volMesh > Internal
The internal field type from which this GeometricField is derived.
Internal & ref(const bool updateAccessTime=true)
Return a reference to the dimensioned internal field.
GeometricBoundaryField< vector, fvsPatchField, surfaceMesh > Boundary
Type of boundary fields.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
const word & name() const noexcept
Return the object name.
dimensionedScalar deltaT() const
Return time step.
T & first()
Return the first element of the list.
const_iterator cend() const noexcept
Return const_iterator to end traversing the constant UList.
volScalarField alpha_
Laminar thermal diffusivity [kg/m/s].
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const surfaceScalarField & nHatf() const
tmp< volScalarField > alphaEff() const
Effective thermal turbulent diffusivity of mixture [kg/m/s].
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
virtual word thermoName() const
Return the name of the thermo physics.
tmp< surfaceScalarField > surfaceTensionForce() const
virtual tmp< volScalarField > Cv() const
Heat capacity at constant volume [J/kg/K].
virtual tmp< volScalarField > rho() const
Density [kg/m^3].
virtual tmp< volScalarField > W() const
Molecular weight [kg/kmol].
virtual tmp< volScalarField > gamma() const
Gamma = Cp/Cv [].
virtual bool incompressible() const
Return true if the equation of state is incompressible.
virtual tmp< volScalarField > Cpv() const
Heat capacity at constant pressure/volume [J/kg/K].
virtual tmp< scalarField > THE(const scalarField &h, const scalarField &p, const scalarField &T0, const labelList &cells) const
Temperature from enthalpy/internal energy for cell-set.
virtual tmp< volScalarField > kappa() const
Thermal diffusivity for temperature of mixture [J/m/s/K].
tmp< volScalarField > rCv() const
Return the phase-averaged reciprocal Cv.
virtual tmp< volScalarField > Cp() const
Heat capacity at constant pressure [J/kg/K].
virtual void correct()
Update properties.
virtual volScalarField & he()
Enthalpy/Internal energy [J/kg].
virtual tmp< volScalarField > hc() const
Chemical enthalpy [J/kg].
virtual tmp< volScalarField > nu() const
Kinematic viscosity of mixture [m^2/s].
void correctRho(const volScalarField &dp)
Update densities for given pressure change.
virtual tmp< volScalarField > alphahe() const
Thermal diffusivity for energy of mixture [kg/m/s].
virtual tmp< volScalarField > CpByCpv() const
Heat capacity ratio [].
void solve()
Solve for the mixture phase-fractions.
virtual bool isochoric() const
Return true if the equation of state is isochoric.
const Time & time() const noexcept
Return time registry.
volScalarField mu_
Dynamic viscosity [kg/m/s].
volScalarField psi_
Compressibility [s^2/m^2].
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
A class for managing temporary objects.
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const volScalarField & mu
const tmp< volScalarField > & tCv
const volScalarField & Cv
const tmp< volScalarField > & tCp
const volScalarField & Cp
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
word alpharScheme("div(phirb,alpha)")
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Calculate the divergence of the given field.
Calculate the face-flux of the given field.
Calculate the gradient of the given field.
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
Calculate the snGrad of the given volField.
volScalarField alphaEff("alphaEff", turbulence->nu()/Pr+alphat)
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
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)
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
const dimensionedScalar kappa
Coulomb constant: default SI units: [N.m2/C2].
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
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.
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
GeometricField< vector, fvPatchField, volMesh > volVectorField
const dimensionSet dimless
Dimensionless.
dimensionedScalar pos0(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
GeometricField< scalar, fvPatchField, volMesh > volScalarField
dimensionedScalar tanh(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Field< vector > vectorField
Specialisation of Field<T> for vector.
static constexpr const zero Zero
Global zero (0)
dimensionedScalar cbrt(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
tmp< volScalarField > trho
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux(phi, alpha1, alphaScheme))
Unit conversion functions.