27 #include "alphaContactAngleFvPatchScalarField.H" 36 #include "surfaceInterpolate.H" 46 const Foam::scalar Foam::multiphaseMixtureThermo::convertToRad =
52 void Foam::multiphaseMixtureThermo::calcAlphas()
57 forAllIter(PtrDictionary<phaseModel>, phases_, phase)
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),
126 forAllIter(PtrDictionary<phaseModel>, phases_, phasei)
137 for (++phasei; phasei != phases_.end(); ++
phasei)
148 forAllIter(PtrDictionary<phaseModel>, phases_, phasei)
159 word name =
phasei().thermo().thermoName();
161 for (++ phasei; phasei != phases_.end(); ++
phasei)
163 name +=
',' +
phasei().thermo().thermoName();
176 ico &= phase().thermo().incompressible();
189 iso &= phase().thermo().incompressible();
206 for (++phasei; phasei != phases_.end(); ++
phasei)
229 for (++phasei; phasei != phases_.end(); ++
phasei)
253 for (++phasei; phasei != phases_.end(); ++
phasei)
269 for (++phasei; phasei != phases_.end(); ++
phasei)
310 for (++phasei; phasei != phases_.end(); ++
phasei)
326 tmp<scalarField>
trho 331 for (++phasei; phasei != phases_.end(); ++
phasei)
347 for (++phasei; phasei != phases_.end(); ++
phasei)
370 for (++phasei; phasei != phases_.end(); ++
phasei)
386 for (++phasei; phasei != phases_.end(); ++
phasei)
409 for (++phasei; phasei != phases_.end(); ++
phasei)
425 for (++phasei; phasei != phases_.end(); ++
phasei)
443 tmp<scalarField> tgamma
448 for (++phasei; phasei != phases_.end(); ++
phasei)
452 *
phasei().thermo().gamma(p, T, patchi);
465 for (++phasei; phasei != phases_.end(); ++
phasei)
483 tmp<scalarField> tCpv
488 for (++phasei; phasei != phases_.end(); ++
phasei)
492 *
phasei().thermo().Cpv(p, T, patchi);
505 for (++phasei; phasei != phases_.end(); ++
phasei)
523 tmp<scalarField> tCpByCpv
528 for (++phasei; phasei != phases_.end(); ++
phasei)
532 *
phasei().thermo().CpByCpv(p, T, patchi);
545 for (++phasei; phasei != phases_.end(); ++
phasei)
565 return mu(patchi)/
rho(patchi);
575 for (++phasei; phasei != phases_.end(); ++
phasei)
591 tmp<scalarField> tkappa
596 for (++phasei; phasei != phases_.end(); ++
phasei)
612 for (++phasei; phasei != phases_.end(); ++
phasei)
630 phasei().boundaryField()[patchi]
634 for (++phasei; phasei != phases_.end(); ++
phasei)
638 *
phasei().thermo().alphahe(patchi);
654 for (++phasei; phasei != phases_.end(); ++
phasei)
656 tkappaEff.ref() +=
phasei()*
phasei().thermo().kappaEff(alphat);
671 tmp<scalarField> tkappaEff
673 phasei().boundaryField()[patchi]
677 for (++phasei; phasei != phases_.end(); ++
phasei)
681 *
phasei().thermo().kappaEff(alphat, patchi);
697 for (++phasei; phasei != phases_.end(); ++
phasei)
716 phasei().boundaryField()[patchi]
720 for (++phasei; phasei != phases_.end(); ++
phasei)
724 *
phasei().thermo().alphaEff(alphat, patchi);
737 for (++phasei; phasei != phases_.end(); ++
phasei)
749 tmp<surfaceScalarField> tstf
755 "surfaceTensionForce",
756 mesh_.time().timeName(),
762 "surfaceTensionForce",
763 dimensionSet(1, -2, -2, 0, 0),
773 const phaseModel& alpha1 =
phase1();
778 for (; phase2 != phases_.end(); ++
phase2)
780 const phaseModel& alpha2 =
phase2();
782 sigmaTable::const_iterator sigma =
783 sigmas_.find(interfacePair(alpha1, alpha2));
785 if (sigma == sigmas_.end())
788 <<
"Cannot find interface " << interfacePair(alpha1, alpha2)
789 <<
" in list of sigma values" 808 const Time& runTime = mesh_.time();
810 const dictionary& alphaControls = mesh_.solverDict(
"alpha");
812 scalar cAlpha(
readScalar(alphaControls.lookup(
"cAlpha")));
824 !(++alphaSubCycle).end();
828 rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
862 return gradAlphaf/(
mag(gradAlphaf) + deltaN_);
873 return nHatfv(alpha1, alpha2) & mesh_.Sf();
883 void Foam::multiphaseMixtureThermo::correctContactAngle
885 const phaseModel& alpha1,
886 const phaseModel& alpha2,
887 surfaceVectorField::Boundary& nHatb
890 const volScalarField::Boundary& gbf
891 = alpha1.boundaryField();
893 const fvBoundaryMesh&
boundary = mesh_.boundary();
897 if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
899 const alphaContactAngleFvPatchScalarField& acap =
900 refCast<const alphaContactAngleFvPatchScalarField>(gbf[
patchi]);
906 mesh_.Sf().boundaryField()[
patchi]
907 /mesh_.magSf().boundaryField()[
patchi]
910 alphaContactAngleFvPatchScalarField::thetaPropsTable::
912 acap.thetaProps().find(interfacePair(alpha1, alpha2));
914 if (tp == acap.thetaProps().end())
917 <<
"Cannot find interface " << interfacePair(alpha1, alpha2)
918 <<
"\n in table of theta properties for patch " 919 << acap.patch().name()
923 bool matched = (tp.key().first() == alpha1.name());
925 scalar theta0 = convertToRad*tp().theta0(matched);
926 scalarField theta(boundary[patchi].size(), theta0);
928 scalar uTheta = tp().uTheta();
933 scalar thetaA = convertToRad*tp().thetaA(matched);
934 scalar thetaR = convertToRad*tp().thetaR(matched);
939 U_.boundaryField()[
patchi].patchInternalField()
940 - U_.boundaryField()[
patchi]
942 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
947 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
951 nWall /= (
mag(nWall) + small);
957 theta += (thetaA - thetaR)*
tanh(uwall/uTheta);
971 b2[facei] =
cos(
acos(a12[facei]) - theta[facei]);
979 nHatPatch = a*AfHatPatch +
b*nHatPatch;
981 nHatPatch /= (
mag(nHatPatch) + deltaN_.
value());
989 const phaseModel& alpha1,
990 const phaseModel& alpha2
993 tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
995 correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef());
998 return -
fvc::div(tnHatfv & mesh_.Sf());
1005 tmp<volScalarField> tnearInt
1012 mesh_.time().timeName(),
1023 max(tnearInt(),
pos0(phase() - 0.01)*
pos0(0.99 - phase()));
1030 void Foam::multiphaseMixtureThermo::solveAlphas
1035 static label nSolves=-1;
1038 word alphaScheme(
"div(phi,alpha)");
1044 PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
1047 forAllIter(PtrDictionary<phaseModel>, phases_, phase)
1049 phaseModel& alpha = phase();
1056 phi_.name() + alpha.name(),
1068 forAllIter(PtrDictionary<phaseModel>, phases_, phase2)
1070 phaseModel& alpha2 =
phase2();
1072 if (&alpha2 == &alpha)
continue;
1086 1.0/mesh_.time().deltaT().value(),
1087 geometricOneField(),
1110 mesh_.time().timeName(),
1123 forAllIter(PtrDictionary<phaseModel>, phases_, phase)
1125 phaseModel& alpha = phase();
1128 alphaPhi += upwind<scalar>(mesh_, phi_).
flux(alpha);
1135 mesh_.time().timeName(),
1147 mesh_.time().timeName(),
1160 if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
1162 Sp[celli] += dgdt[celli]*alpha[celli];
1163 Su[celli] -= dgdt[celli]*alpha[celli];
1165 else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
1167 Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
1174 const phaseModel& alpha2 =
phase2();
1176 if (&alpha2 == &alpha)
continue;
1182 if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
1184 Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
1185 Su[celli] += dgdt2[celli]*alpha[celli];
1187 else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
1189 Sp[celli] += dgdt2[celli]*alpha2[celli];
1196 geometricOneField(),
1205 Info<< alpha.name() <<
" volume fraction, min, max = " 1206 << alpha.weightedAverage(mesh_.V()).value()
1207 <<
' ' <<
min(alpha).value()
1208 <<
' ' <<
max(alpha).value()
1216 Info<<
"Phase-sum volume fraction, min, max = " 1217 << sumAlpha.weightedAverage(mesh_.V()).value()
1218 <<
' ' <<
min(sumAlpha).value()
1219 <<
' ' <<
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.
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.
Info<< "Creating turbulence model\"<< endl;tmp< volScalarField > talphaEff
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
virtual tmp< volScalarField > hc() const
Chemical enthalpy [J/kg].
word alpharScheme("div(phirb,alpha)")
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
tmp< volScalarField > trho
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
Calculate the snGrad of the given volField.
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
const dimensionedScalar kappa
Coulomb constant: default SI units: [N.m2/C2].
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.
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
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.
const Type & value() const
Return const reference to value.
List< label > labelList
A List of labels.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
faceListList boundary(nPatches)
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
label readLabel(Istream &is)
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(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil) - psi) *pSat, rhoMin);# 1 "/home/ubuntu/OpenFOAM-6/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()
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 dimensionedScalar mu
Atomic mass unit.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
virtual tmp< volScalarField > kappaEff(const volScalarField &alphat) const
Effective thermal diffusivity of mixture [J/m/s/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)
volScalarField & h
Planck constant.
surfaceScalarField phic(mixture.cAlpha() *mag(phi/mesh.magSf()))
virtual tmp< volScalarField > W() const
Molecular weight [kg/kmol].
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const scalarList W(::W(thermo))
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
MULES: Multidimensional universal limiter for explicit solution.
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.
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 alphaEff("alphaEff", turbulence->nu()/Pr+alphat)
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 [J/m/s/K].
virtual tmp< volScalarField > Cpv() const
Heat capacity at constant pressure/volume [J/kg/K].
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux(phi, alpha1, alphaScheme))
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
#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 [J/m/s/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.