From d89a8d3daaa04aad02c2fa49052e66e24b2d79a3 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Thu, 1 Oct 2020 15:24:14 +0100 Subject: [PATCH] compressibleMultiphaseInterFoam: updated mixing of thermophysical properties Thermodynamic properties are now mass-fraction mixed Transport properties remain volume-fraction mixed --- .../compressibleMultiphaseInterFoam/TEqn.H | 1 + .../createFields.H | 12 +- .../multiphaseMixtureThermo.C | 171 +++++++----------- .../multiphaseMixtureThermo.H | 12 +- .../phaseModel/phaseModel.C | 13 +- .../phaseModel/phaseModel.H | 13 ++ .../compressibleMultiphaseInterFoam/pEqn.H | 13 +- .../laminar/damBreak4phase/0/p_rgh | 1 + .../laminar/damBreak4phase/system/fvSchemes | 2 +- 9 files changed, 106 insertions(+), 132 deletions(-) diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H index 66274441c3..546b232d44 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H @@ -13,5 +13,6 @@ TEqn.relax(); TEqn.solve(); + mixture.correctThermo(); mixture.correct(); } diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H index f8728472a0..aae7d1782a 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H @@ -31,17 +31,7 @@ volVectorField U Info<< "Constructing multiphaseMixtureThermo\n" << endl; multiphaseMixtureThermo mixture(U, phi); -volScalarField rho -( - IOobject - ( - "rho", - runTime.timeName(), - mesh, - IOobject::READ_IF_PRESENT - ), - mixture.rho() -); +const volScalarField& rho = mixture.rho(); dimensionedScalar pMin("pMin", dimPressure, mixture); diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C index 2a0b2b382e..d5b8aefdb0 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C @@ -67,7 +67,7 @@ Foam::multiphaseMixtureThermo::multiphaseMixtureThermo const surfaceScalarField& phi ) : - psiThermo::composite(U.mesh(), word::null), + rhoThermo::composite(U.mesh(), word::null), phases_(lookup("phases"), phaseModel::iNew(p_, T_)), mesh_(U.mesh()), @@ -118,25 +118,36 @@ Foam::multiphaseMixtureThermo::multiphaseMixtureThermo // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -void Foam::multiphaseMixtureThermo::correct() +void Foam::multiphaseMixtureThermo::correctThermo() { forAllIter(PtrDictionary, phases_, phasei) { phasei().correct(); } +} + +void Foam::multiphaseMixtureThermo::correct() +{ PtrDictionary::iterator phasei = phases_.begin(); + rho_ = phasei()*phasei().thermo().rho(); psi_ = phasei()*phasei().thermo().psi(); mu_ = phasei()*phasei().thermo().mu(); alpha_ = phasei()*phasei().thermo().alpha(); for (++phasei; phasei != phases_.end(); ++phasei) { + rho_ += phasei()*phasei().thermo().rho(); psi_ += phasei()*phasei().thermo().psi(); mu_ += phasei()*phasei().thermo().mu(); alpha_ += phasei()*phasei().thermo().alpha(); } + + forAllIter(PtrDictionary, phases_, phasei) + { + phasei().Alpha() = phasei()*phasei().thermo().rho()/rho_; + } } @@ -144,7 +155,7 @@ void Foam::multiphaseMixtureThermo::correctRho(const volScalarField& dp) { forAllIter(PtrDictionary, phases_, phasei) { - phasei().thermo().rho() += phasei().thermo().psi()*dp; + phasei().thermo().rho() += phasei().thermo().psi()*dp; } } @@ -198,11 +209,11 @@ Foam::tmp Foam::multiphaseMixtureThermo::he { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp the(phasei()*phasei().thermo().he(p, T)); + tmp the(phasei().Alpha()*phasei().thermo().he(p, T)); for (++phasei; phasei != phases_.end(); ++phasei) { - the.ref() += phasei()*phasei().thermo().he(p, T); + the.ref() += phasei().Alpha()*phasei().thermo().he(p, T); } return the; @@ -242,13 +253,14 @@ Foam::tmp Foam::multiphaseMixtureThermo::he tmp the ( - phasei().boundaryField()[patchi]*phasei().thermo().he(T, patchi) + phasei().Alpha().boundaryField()[patchi]*phasei().thermo().he(T, patchi) ); for (++phasei; phasei != phases_.end(); ++phasei) { the.ref() += - phasei().boundaryField()[patchi]*phasei().thermo().he(T, patchi); + phasei().Alpha().boundaryField()[patchi] + *phasei().thermo().he(T, patchi); } return the; @@ -259,11 +271,11 @@ Foam::tmp Foam::multiphaseMixtureThermo::hs() const { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp ths(phasei()*phasei().thermo().hs()); + tmp ths(phasei().Alpha()*phasei().thermo().hs()); for (++phasei; phasei != phases_.end(); ++phasei) { - ths.ref() += phasei()*phasei().thermo().hs(); + ths.ref() += phasei().Alpha()*phasei().thermo().hs(); } return ths; @@ -278,11 +290,11 @@ Foam::tmp Foam::multiphaseMixtureThermo::hs { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp ths(phasei()*phasei().thermo().hs(p, T)); + tmp ths(phasei().Alpha()*phasei().thermo().hs(p, T)); for (++phasei; phasei != phases_.end(); ++phasei) { - ths.ref() += phasei()*phasei().thermo().hs(p, T); + ths.ref() += phasei().Alpha()*phasei().thermo().hs(p, T); } return ths; @@ -322,13 +334,14 @@ Foam::tmp Foam::multiphaseMixtureThermo::hs tmp ths ( - phasei().boundaryField()[patchi]*phasei().thermo().hs(T, patchi) + phasei().Alpha().boundaryField()[patchi]*phasei().thermo().hs(T, patchi) ); for (++phasei; phasei != phases_.end(); ++phasei) { ths.ref() += - phasei().boundaryField()[patchi]*phasei().thermo().hs(T, patchi); + phasei().Alpha().boundaryField()[patchi] + *phasei().thermo().hs(T, patchi); } return ths; @@ -339,11 +352,11 @@ Foam::tmp Foam::multiphaseMixtureThermo::ha() const { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp tha(phasei()*phasei().thermo().ha()); + tmp tha(phasei().Alpha()*phasei().thermo().ha()); for (++phasei; phasei != phases_.end(); ++phasei) { - tha.ref() += phasei()*phasei().thermo().ha(); + tha.ref() += phasei().Alpha()*phasei().thermo().ha(); } return tha; @@ -358,11 +371,11 @@ Foam::tmp Foam::multiphaseMixtureThermo::ha { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp tha(phasei()*phasei().thermo().ha(p, T)); + tmp tha(phasei().Alpha()*phasei().thermo().ha(p, T)); for (++phasei; phasei != phases_.end(); ++phasei) { - tha.ref() += phasei()*phasei().thermo().ha(p, T); + tha.ref() += phasei().Alpha()*phasei().thermo().ha(p, T); } return tha; @@ -402,13 +415,14 @@ Foam::tmp Foam::multiphaseMixtureThermo::ha tmp tha ( - phasei().boundaryField()[patchi]*phasei().thermo().ha(T, patchi) + phasei().Alpha().boundaryField()[patchi]*phasei().thermo().ha(T, patchi) ); for (++phasei; phasei != phases_.end(); ++phasei) { tha.ref() += - phasei().boundaryField()[patchi]*phasei().thermo().ha(T, patchi); + phasei().Alpha().boundaryField()[patchi] + *phasei().thermo().ha(T, patchi); } return tha; @@ -419,11 +433,11 @@ Foam::tmp Foam::multiphaseMixtureThermo::hc() const { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp thc(phasei()*phasei().thermo().hc()); + tmp thc(phasei().Alpha()*phasei().thermo().hc()); for (++phasei; phasei != phases_.end(); ++phasei) { - thc.ref() += phasei()*phasei().thermo().hc(); + thc.ref() += phasei().Alpha()*phasei().thermo().hc(); } return thc; @@ -466,52 +480,15 @@ Foam::tmp Foam::multiphaseMixtureThermo::THE } -Foam::tmp Foam::multiphaseMixtureThermo::rho() const -{ - PtrDictionary::const_iterator phasei = phases_.begin(); - - tmp trho(phasei()*phasei().thermo().rho()); - - for (++phasei; phasei != phases_.end(); ++phasei) - { - trho.ref() += phasei()*phasei().thermo().rho(); - } - - return trho; -} - - -Foam::tmp Foam::multiphaseMixtureThermo::rho -( - const label patchi -) const -{ - PtrDictionary::const_iterator phasei = phases_.begin(); - - tmp trho - ( - phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi) - ); - - for (++phasei; phasei != phases_.end(); ++phasei) - { - trho.ref() += - phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi); - } - - return trho; -} - - Foam::tmp Foam::multiphaseMixtureThermo::Cp() const { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp tCp(phasei()*phasei().thermo().Cp()); + tmp tCp(phasei().Alpha()*phasei().thermo().Cp()); for (++phasei; phasei != phases_.end(); ++phasei) { - tCp.ref() += phasei()*phasei().thermo().Cp(); + tCp.ref() += phasei().Alpha()*phasei().thermo().Cp(); } return tCp; @@ -528,13 +505,14 @@ Foam::tmp Foam::multiphaseMixtureThermo::Cp tmp tCp ( - phasei().boundaryField()[patchi]*phasei().thermo().Cp(T, patchi) + phasei().Alpha().boundaryField()[patchi]*phasei().thermo().Cp(T, patchi) ); for (++phasei; phasei != phases_.end(); ++phasei) { tCp.ref() += - phasei().boundaryField()[patchi]*phasei().thermo().Cp(T, patchi); + phasei().Alpha().boundaryField()[patchi] + *phasei().thermo().Cp(T, patchi); } return tCp; @@ -545,11 +523,11 @@ Foam::tmp Foam::multiphaseMixtureThermo::Cv() const { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp tCv(phasei()*phasei().thermo().Cv()); + tmp tCv(phasei().Alpha()*phasei().thermo().Cv()); for (++phasei; phasei != phases_.end(); ++phasei) { - tCv.ref() += phasei()*phasei().thermo().Cv(); + tCv.ref() += phasei().Alpha()*phasei().thermo().Cv(); } return tCv; @@ -566,13 +544,14 @@ Foam::tmp Foam::multiphaseMixtureThermo::Cv tmp tCv ( - phasei().boundaryField()[patchi]*phasei().thermo().Cv(T, patchi) + phasei().Alpha().boundaryField()[patchi]*phasei().thermo().Cv(T, patchi) ); for (++phasei; phasei != phases_.end(); ++phasei) { tCv.ref() += - phasei().boundaryField()[patchi]*phasei().thermo().Cv(T, patchi); + phasei().Alpha().boundaryField()[patchi] + *phasei().thermo().Cv(T, patchi); } return tCv; @@ -581,16 +560,7 @@ Foam::tmp Foam::multiphaseMixtureThermo::Cv Foam::tmp Foam::multiphaseMixtureThermo::gamma() const { - PtrDictionary::const_iterator phasei = phases_.begin(); - - tmp tgamma(phasei()*phasei().thermo().gamma()); - - for (++phasei; phasei != phases_.end(); ++phasei) - { - tgamma.ref() += phasei()*phasei().thermo().gamma(); - } - - return tgamma; + return Cp()/Cv(); } @@ -600,21 +570,7 @@ Foam::tmp Foam::multiphaseMixtureThermo::gamma const label patchi ) const { - PtrDictionary::const_iterator phasei = phases_.begin(); - - tmp tgamma - ( - phasei().boundaryField()[patchi]*phasei().thermo().gamma(T, patchi) - ); - - for (++phasei; phasei != phases_.end(); ++phasei) - { - tgamma.ref() += - phasei().boundaryField()[patchi] - *phasei().thermo().gamma(T, patchi); - } - - return tgamma; + return Cp(T, patchi)/Cv(T, patchi); } @@ -622,11 +578,11 @@ Foam::tmp Foam::multiphaseMixtureThermo::Cpv() const { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp tCpv(phasei()*phasei().thermo().Cpv()); + tmp tCpv(phasei().Alpha()*phasei().thermo().Cpv()); for (++phasei; phasei != phases_.end(); ++phasei) { - tCpv.ref() += phasei()*phasei().thermo().Cpv(); + tCpv.ref() += phasei().Alpha()*phasei().thermo().Cpv(); } return tCpv; @@ -643,7 +599,8 @@ Foam::tmp Foam::multiphaseMixtureThermo::Cpv tmp tCpv ( - phasei().boundaryField()[patchi]*phasei().thermo().Cpv(T, patchi) + phasei().Alpha().boundaryField()[patchi] + *phasei().thermo().Cpv(T, patchi) ); for (++phasei; phasei != phases_.end(); ++phasei) @@ -661,11 +618,11 @@ Foam::tmp Foam::multiphaseMixtureThermo::CpByCpv() const { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp tCpByCpv(phasei()*phasei().thermo().CpByCpv()); + tmp tCpByCpv(phasei().Alpha()*phasei().thermo().CpByCpv()); for (++phasei; phasei != phases_.end(); ++phasei) { - tCpByCpv.ref() += phasei()*phasei().thermo().CpByCpv(); + tCpByCpv.ref() += phasei().Alpha()*phasei().thermo().CpByCpv(); } return tCpByCpv; @@ -682,7 +639,8 @@ Foam::tmp Foam::multiphaseMixtureThermo::CpByCpv tmp tCpByCpv ( - phasei().boundaryField()[patchi]*phasei().thermo().CpByCpv(T, patchi) + phasei().Alpha().boundaryField()[patchi] + *phasei().thermo().CpByCpv(T, patchi) ); for (++phasei; phasei != phases_.end(); ++phasei) @@ -700,14 +658,14 @@ Foam::tmp Foam::multiphaseMixtureThermo::W() const { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp tW(phasei()*phasei().thermo().W()); + volScalarField rW(phasei().Alpha()/phasei().thermo().W()); for (++phasei; phasei != phases_.end(); ++phasei) { - tW.ref() += phasei()*phasei().thermo().W(); + rW += phasei().Alpha()/phasei().thermo().W(); } - return tW; + return 1/rW; } @@ -718,18 +676,19 @@ Foam::tmp Foam::multiphaseMixtureThermo::W { PtrDictionary::const_iterator phasei = phases_.begin(); - tmp tW + scalarField rW ( - phasei().boundaryField()[patchi]*phasei().thermo().W(patchi) + phasei().Alpha().boundaryField()[patchi]/phasei().thermo().W(patchi) ); for (++phasei; phasei != phases_.end(); ++phasei) { - tW.ref() += - phasei().boundaryField()[patchi]*phasei().thermo().W(patchi); + rW += + phasei().Alpha().boundaryField()[patchi] + /phasei().thermo().W(patchi); } - return tW; + return 1/rW; } @@ -1006,6 +965,8 @@ void Foam::multiphaseMixtureThermo::solve() { solveAlphas(cAlpha); } + + correct(); } diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H index 0c61f0a291..76bf51257d 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H @@ -39,7 +39,6 @@ SourceFiles #include "volFields.H" #include "surfaceFields.H" #include "rhoThermo.H" -#include "psiThermo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -52,7 +51,7 @@ namespace Foam class multiphaseMixtureThermo : - public psiThermo::composite + public rhoThermo::composite { public: @@ -229,6 +228,9 @@ public: return rhoPhi_; } + //- Correct the thermodynamics of each phase + virtual void correctThermo(); + //- Update properties virtual void correct(); @@ -366,12 +368,6 @@ public: // Fields derived from thermodynamic state variables - //- Density [kg/m^3] - virtual tmp rho() const; - - //- Density for patch [kg/m^3] - virtual tmp rho(const label patchi) const; - //- Heat capacity at constant pressure [J/kg/K] virtual tmp Cp() const; diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C index 7f75ad869a..a60a069105 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2013-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2020 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -50,6 +50,17 @@ Foam::phaseModel::phaseModel p_(p), T_(T), thermo_(nullptr), + Alpha_ + ( + IOobject + ( + IOobject::groupName("Alpha", phaseName), + p.mesh().time().timeName(), + p.mesh() + ), + p.mesh(), + dimensionedScalar(dimless, 0) + ), dgdt_ ( IOobject diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H index 0ab051b7bc..fb5bb5ea66 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H @@ -61,6 +61,7 @@ class phaseModel const volScalarField& p_; const volScalarField& T_; autoPtr thermo_; + volScalarField Alpha_; volScalarField dgdt_; @@ -129,6 +130,18 @@ public: return thermo_(); } + //- Return const-access to phase mass-fraction + const volScalarField& Alpha() const + { + return Alpha_; + } + + //- Return access to phase mass-fraction + volScalarField& Alpha() + { + return Alpha_; + } + //- Return const-access to phase divergence const volScalarField& dgdt() const { diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H index 7265f327d9..fc3c1d017b 100644 --- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H @@ -68,7 +68,7 @@ phase ) { - tmp hmm + tmp p_rghEqnCompi ( (max(phase(), scalar(0))/phase().thermo().rho()) *p_rghEqnComps[phasei] @@ -76,11 +76,11 @@ if (phasei == 0) { - p_rghEqnComp = hmm; + p_rghEqnComp = p_rghEqnCompi; } else { - p_rghEqnComp.ref() += hmm; + p_rghEqnComp.ref() += p_rghEqnCompi; } phasei++; @@ -90,6 +90,9 @@ if (pimple.finalNonOrthogonalIter()) { + p = max(p_rgh + mixture.rho()*gh, pMin); + p_rgh = p - rho*gh; + phasei = 0; forAllIter ( @@ -113,11 +116,9 @@ } } - p = max(p_rgh + mixture.rho()*gh, pMin); - // Update densities from change in p_rgh mixture.correctRho(p_rgh - p_rgh_0); - rho = mixture.rho(); + mixture.correct(); // Correct p_rgh for consistency with p and the updated densities p_rgh = p - rho*gh; diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/p_rgh b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/p_rgh index 983378041e..edb028a263 100644 --- a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/p_rgh +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/0/p_rgh @@ -42,6 +42,7 @@ boundaryField { type totalPressure; p0 uniform 1e5; + rho thermo:rho; } defaultFaces diff --git a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes index e2cb769546..957f0ac871 100644 --- a/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes +++ b/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/fvSchemes @@ -34,7 +34,7 @@ divSchemes div(rhoPhi,T) Gauss upwind; div(rhoPhi,K) Gauss upwind; div(phi,p) Gauss upwind; - div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; + div(((thermo:rho*nuEff)*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes