27 #include "alphaContactAngleFvPatchScalarField.H" 37 #include "surfaceInterpolate.H" 49 void Foam::multiphaseMixtureThermo::calcAlphas()
54 forAllIter(PtrDictionary<phaseModel>, phases_, phase)
56 alphas_ += level*phase();
70 psiThermo(U.
mesh(), word::null),
71 phases_(
lookup(
"phases"), phaseModel::iNew(p_, T_)),
105 sigmas_(
lookup(
"sigmas")),
106 dimSigma_(1, 0, -2, 0, 0),
123 forAllIter(PtrDictionary<phaseModel>, phases_, phasei)
134 for (++phasei; phasei != phases_.end(); ++
phasei)
145 forAllIter(PtrDictionary<phaseModel>, phases_, phasei)
156 word name =
phasei().thermo().thermoName();
158 for (++ phasei; phasei != phases_.end(); ++
phasei)
160 name +=
',' +
phasei().thermo().thermoName();
173 ico &= phase().thermo().incompressible();
186 iso &= phase().thermo().incompressible();
203 for (++phasei; phasei != phases_.end(); ++
phasei)
225 for (++phasei; phasei != phases_.end(); ++
phasei)
248 for (++phasei; phasei != phases_.end(); ++
phasei)
264 for (++phasei; phasei != phases_.end(); ++
phasei)
283 for (++phasei; phasei != phases_.end(); ++
phasei)
305 for (++phasei; phasei != phases_.end(); ++
phasei)
328 for (++phasei; phasei != phases_.end(); ++
phasei)
344 for (++phasei; phasei != phases_.end(); ++
phasei)
363 for (++phasei; phasei != phases_.end(); ++
phasei)
385 for (++phasei; phasei != phases_.end(); ++
phasei)
408 for (++phasei; phasei != phases_.end(); ++
phasei)
424 for (++phasei; phasei != phases_.end(); ++
phasei)
463 for (++phasei; phasei != phases_.end(); ++
phasei)
479 tmp<scalarField>
trho 484 for (++phasei; phasei != phases_.end(); ++
phasei)
500 for (++phasei; phasei != phases_.end(); ++
phasei)
522 for (++phasei; phasei != phases_.end(); ++
phasei)
538 for (++phasei; phasei != phases_.end(); ++
phasei)
560 for (++phasei; phasei != phases_.end(); ++
phasei)
576 for (++phasei; phasei != phases_.end(); ++
phasei)
593 tmp<scalarField> tgamma
598 for (++phasei; phasei != phases_.end(); ++
phasei)
602 *
phasei().thermo().gamma(T, patchi);
615 for (++phasei; phasei != phases_.end(); ++
phasei)
632 tmp<scalarField> tCpv
637 for (++phasei; phasei != phases_.end(); ++
phasei)
641 *
phasei().thermo().Cpv(T, patchi);
654 for (++phasei; phasei != phases_.end(); ++
phasei)
671 tmp<scalarField> tCpByCpv
676 for (++phasei; phasei != phases_.end(); ++
phasei)
680 *
phasei().thermo().CpByCpv(T, patchi);
693 for (++phasei; phasei != phases_.end(); ++
phasei)
714 for (++phasei; phasei != phases_.end(); ++
phasei)
735 return mu(patchi)/
rho(patchi);
745 for (++phasei; phasei != phases_.end(); ++
phasei)
761 tmp<scalarField> tkappa
766 for (++phasei; phasei != phases_.end(); ++
phasei)
782 for (++phasei; phasei != phases_.end(); ++
phasei)
784 talphaEff.ref() +=
phasei()*
phasei().thermo().alphahe();
798 tmp<scalarField> talphaEff
800 phasei().boundaryField()[patchi]
804 for (++phasei; phasei != phases_.end(); ++
phasei)
808 *
phasei().thermo().alphahe(patchi);
824 for (++phasei; phasei != phases_.end(); ++
phasei)
826 tkappaEff.ref() +=
phasei()*
phasei().thermo().kappaEff(alphat);
841 tmp<scalarField> tkappaEff
843 phasei().boundaryField()[patchi]
847 for (++phasei; phasei != phases_.end(); ++
phasei)
851 *
phasei().thermo().kappaEff(alphat, patchi);
867 for (++phasei; phasei != phases_.end(); ++
phasei)
869 talphaEff.ref() +=
phasei()*
phasei().thermo().alphaEff(alphat);
884 tmp<scalarField> talphaEff
886 phasei().boundaryField()[patchi]
890 for (++phasei; phasei != phases_.end(); ++
phasei)
894 *
phasei().thermo().alphaEff(alphat, patchi);
907 for (++phasei; phasei != phases_.end(); ++
phasei)
919 tmp<surfaceScalarField> tstf
923 "surfaceTensionForce",
933 const phaseModel& alpha1 = phase1();
938 for (; phase2 != phases_.end(); ++phase2)
940 const phaseModel& alpha2 = phase2();
942 sigmaTable::const_iterator sigma =
943 sigmas_.find(interfacePair(alpha1, alpha2));
945 if (sigma == sigmas_.end())
948 <<
"Cannot find interface " << interfacePair(alpha1, alpha2)
949 <<
" in list of sigma values" 968 const Time& runTime = mesh_.time();
970 const dictionary& alphaControls = mesh_.solverDict(
"alpha");
972 scalar cAlpha(alphaControls.lookup<scalar>(
"cAlpha"));
984 !(++alphaSubCycle).end();
988 rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
1022 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
1033 return nHatfv(alpha1, alpha2) & mesh_.Sf();
1043 void Foam::multiphaseMixtureThermo::correctContactAngle
1045 const phaseModel& alpha1,
1046 const phaseModel& alpha2,
1047 surfaceVectorField::Boundary& nHatb
1050 const volScalarField::Boundary& gbf
1051 = alpha1.boundaryField();
1053 const fvBoundaryMesh&
boundary = mesh_.boundary();
1057 if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
1059 const alphaContactAngleFvPatchScalarField& acap =
1060 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
1066 mesh_.Sf().boundaryField()[
patchi]
1067 /mesh_.magSf().boundaryField()[
patchi]
1070 alphaContactAngleFvPatchScalarField::thetaPropsTable::
1072 acap.thetaProps().find(interfacePair(alpha1, alpha2));
1074 if (tp == acap.thetaProps().end())
1077 <<
"Cannot find interface " << interfacePair(alpha1, alpha2)
1078 <<
"\n in table of theta properties for patch " 1079 << acap.patch().name()
1083 bool matched = (tp.key().first() == alpha1.name());
1085 scalar theta0 =
degToRad(tp().theta0(matched));
1086 scalarField theta(boundary[patchi].size(), theta0);
1088 scalar uTheta = tp().uTheta();
1093 scalar thetaA =
degToRad(tp().thetaA(matched));
1094 scalar thetaR =
degToRad(tp().thetaR(matched));
1099 U_.boundaryField()[
patchi].patchInternalField()
1100 - U_.boundaryField()[
patchi]
1102 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
1107 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
1111 nWall /= (
mag(nWall) + small);
1117 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
1131 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
1139 nHatPatch = a*AfHatPatch +
b*nHatPatch;
1141 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
1149 const phaseModel& alpha1,
1150 const phaseModel& alpha2
1153 tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
1155 correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef());
1158 return -
fvc::div(tnHatfv & mesh_.Sf());
1165 tmp<volScalarField> tnearInt
1178 max(tnearInt(),
pos0(phase() - 0.01)*
pos0(0.99 - phase()));
1185 void Foam::multiphaseMixtureThermo::solveAlphas
1190 static label nSolves=-1;
1197 phic =
min(cAlpha*phic,
max(phic));
1199 PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
1202 forAllIter(PtrDictionary<phaseModel>, phases_, phase)
1204 phaseModel& alpha = phase();
1211 phi_.name() + alpha.name(),
1223 forAllIter(PtrDictionary<phaseModel>, phases_, phase2)
1225 phaseModel& alpha2 = phase2();
1227 if (&alpha2 == &alpha)
continue;
1241 1.0/mesh_.time().deltaT().value(),
1242 geometricOneField(),
1265 mesh_.time().timeName(),
1278 forAllIter(PtrDictionary<phaseModel>, phases_, phase)
1280 phaseModel& alpha = phase();
1283 alphaPhi += upwind<scalar>(mesh_, phi_).
flux(alpha);
1290 mesh_.time().timeName(),
1302 mesh_.time().timeName(),
1315 if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
1317 Sp[celli] += dgdt[celli]*alpha[celli];
1318 Su[celli] -= dgdt[celli]*alpha[celli];
1320 else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
1322 Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
1329 const phaseModel& alpha2 = phase2();
1331 if (&alpha2 == &alpha)
continue;
1337 if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
1339 Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
1340 Su[celli] += dgdt2[celli]*alpha[celli];
1342 else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
1344 Sp[celli] += dgdt2[celli]*alpha2[celli];
1351 geometricOneField(),
1360 Info<< alpha.name() <<
" volume fraction, min, max = " 1361 << alpha.weightedAverage(mesh_.V()).value()
1362 <<
' ' <<
min(alpha).value()
1363 <<
' ' <<
max(alpha).value()
1371 Info<<
"Phase-sum volume fraction, min, max = " 1372 << sumAlpha.weightedAverage(mesh_.V()).value()
1373 <<
' ' <<
min(sumAlpha).value()
1374 <<
' ' <<
max(sumAlpha).value()
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)
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.
const dimensionedScalar & kappa
Coulomb constant: default SI units: [N.m2/C2].
errorManipArg< error, int > exit(error &err, const int errNo=1)
virtual tmp< volScalarField > gamma() const
Gamma = Cp/Cv [].
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].
tmp< volScalarField > tCv
const volScalarField & Cv
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
virtual tmp< volScalarField > hc() const
Enthalpy of formation [J/kg].
static tmp< GeometricField< scalar, fvsPatchField, surfaceMesh > > New(const word &name, const Internal &, const PtrList< fvsPatchField< scalar >> &)
Return a temporary field constructed from name,.
virtual tmp< volScalarField > ha() const
Absolute enthalpy [J/kg].
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< volScalarField > trho
Calculate the snGrad of the given volField.
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
dimensionedScalar det(const dimensionedSphericalTensor &dt)
rhoReactionThermo & thermo
GeometricField< vector, fvPatchField, volMesh > volVectorField
T & first()
Return the first element of the list.
DimensionedField< scalar, volMesh > Internal
Type of the internal field from which this GeometricField is derived.
CGAL::Exact_predicates_exact_constructions_kernel K
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
virtual tmp< volScalarField > Cp() const
Heat capacity at constant pressure [J/kg/K].
virtual tmp< volScalarField > CpByCpv() const
Heat capacity ratio [].
GeometricField< scalar, fvPatchField, volMesh > volScalarField
void correctRho(const volScalarField &dp)
Update densities for given pressure change.
void solve()
Solve for the mixture phase-fractions.
stressControl lookup("compactNormalStress") >> compactNormalStress
dimensionedScalar cos(const dimensionedScalar &ds)
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
virtual tmp< volScalarField > alphahe() const
Thermal diffusivity for energy of mixture [kg/m/s].
virtual tmp< volScalarField > rho() const
Density [kg/m^3].
const dimensionedScalar & b
Wien displacement law constant: default SI units: [m K].
Calculate the gradient of the given field.
tmp< volScalarField > rCv() const
Return the phase-averaged reciprocal Cv.
A class for handling words, derived from string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Calculate the face-flux of the given field.
virtual tmp< scalarField > THE(const scalarField &h, const scalarField &T0, const labelList &cells) const
Temperature from enthalpy/internal energy for cell-set.
const Type & value() const
Return const reference to value.
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux(phi, alpha1, alphaScheme))
List< label > labelList
A List of labels.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
faceListList boundary(nPatches)
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Calculate the divergence of the given field.
dimensionedScalar pos0(const dimensionedScalar &ds)
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-8/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
const label nAlphaSubCycles(alphaControls.lookup< label >("nAlphaSubCycles"))
defineTypeNameAndDebug(combustionModel, 0)
virtual tmp< volScalarField > nu() const
Kinematic viscosity of mixture [m^2/s].
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
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.
Internal & ref()
Return a reference to the dimensioned internal field.
const word alphaScheme(mesh.divScheme(divAlphaName)[1].wordToken())
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
virtual tmp< volScalarField > kappaEff(const volScalarField &alphat) const
Effective thermal diffusivity of mixture [W/m/K].
multiphaseMixtureThermo(const volVectorField &U, const surfaceScalarField &phi)
Construct from components.
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)
const dimensionedScalar & mu
Atomic mass unit.
virtual tmp< volScalarField > W() const
Molecular weight [kg/kmol].
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const scalarList W(::W(thermo))
scalar Cp(const scalar p, const scalar T) const
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
word alpharScheme("div(phirb,alpha)")
virtual tmp< volScalarField > hs() const
Sensible enthalpy [J/kg].
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
dimensioned< scalar > mag(const dimensioned< Type > &)
const dimensionedScalar & h
Planck constant.
Field< vector > vectorField
Specialisation of Field<T> for vector.
MULES: Multidimensional universal limiter for explicit solution.
const doubleScalar e
Elementary charge.
virtual tmp< volScalarField > Cv() const
Heat capacity at constant volume [J/kg/K].
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
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.
tmp< surfaceScalarField > surfaceTensionForce() const
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
virtual tmp< volScalarField > alphaEff(const volScalarField &alphat) const
Effective thermal diffusivity of mixture [W/m/K].
virtual tmp< volScalarField > Cpv() const
Heat capacity at constant pressure/volume [J/kg/K].
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
virtual tmp< volScalarField > kappa() const
Thermal diffusivity for temperature of mixture [W/m/K].
virtual bool isochoric() const
Return true if the equation of state is isochoric.
virtual word thermoName() const
Return the name of the thermo physics.
virtual volScalarField & he()
Enthalpy/Internal energy [J/kg].
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
friend class const_iterator
virtual bool incompressible() const
Return true if the equation of state is incompressible.
virtual void correct()
Update properties.