From bf23c7e9bb6a7b3418cecc649cea73592b6570c6 Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 9 Jan 2013 14:38:24 +0000 Subject: [PATCH] compressibleTwoPhaseEulerFoam: Updated thermodynamics --- .../Make/options | 3 + .../compressibleTwoPhaseEulerFoam/TEqns.H | 36 --- .../compressibleTwoPhaseEulerFoam/UEqns.H | 6 +- .../compressibleTwoPhaseEulerFoam/alphaEqn.H | 26 ++- .../compressibleTwoPhaseEulerFoam.C | 3 +- .../createFields.H | 208 +++++------------- .../createRASTurbulence.H | 10 +- .../interfacialCoeffs.H | 4 +- .../interfacialModels/Make/options | 6 +- .../RanzMarshall/RanzMarshall.C | 3 +- .../kineticTheoryModels/Make/options | 1 + .../kineticTheoryModel/kineticTheoryModel.C | 6 +- .../compressibleTwoPhaseEulerFoam/pEqn.H | 32 ++- .../phaseModel/Make/options | 5 +- .../phaseModel/phaseModel/phaseModel.C | 43 ++-- .../phaseModel/phaseModel/phaseModel.H | 44 ++-- .../turbulenceModel/kEpsilon.H | 2 +- .../turbulenceModel/wallFunctions.H | 6 +- .../turbulenceModel/wallViscosity.H | 10 +- 19 files changed, 172 insertions(+), 282 deletions(-) delete mode 100644 applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options index 16ce2e3674..d90ba7bbbd 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options @@ -1,4 +1,5 @@ EXE_INC = \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -IturbulenceModel \ @@ -8,6 +9,8 @@ EXE_INC = \ -Iaveraging EXE_LIBS = \ + -lfluidThermophysicalModels \ + -lspecie \ -lincompressibleTransportModels \ -lcompressiblePhaseModel \ -lcompressibleEulerianInterfacialModels \ diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H deleted file mode 100644 index 8f38ca9d91..0000000000 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H +++ /dev/null @@ -1,36 +0,0 @@ -{ - volScalarField kByCp1("kByCp1", alpha1*(k1/Cp1/rho1 + sqr(Ct)*nut2/Prt)); - volScalarField kByCp2("kByCp2", alpha2*(k2/Cp2/rho2 + nut2/Prt)); - - fvScalarMatrix T1Eqn - ( - fvm::ddt(alpha1, T1) - + fvm::div(alphaPhi1, T1) - - fvm::laplacian(kByCp1, T1) - == - heatTransferCoeff*T2/Cp1/rho1 - - fvm::Sp(heatTransferCoeff/Cp1/rho1, T1) - + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))/Cp1 - ); - - fvScalarMatrix T2Eqn - ( - fvm::ddt(alpha2, T2) - + fvm::div(alphaPhi2, T2) - - fvm::laplacian(kByCp2, T2) - == - heatTransferCoeff*T1/Cp2/rho2 - - fvm::Sp(heatTransferCoeff/Cp2/rho2, T2) - + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))/Cp2 - ); - - T1Eqn.relax(); - T1Eqn.solve(); - - T2Eqn.relax(); - T2Eqn.solve(); - - // Update compressibilities - psi1 = 1.0/(R1*T1); - psi2 = 1.0/(R2*T2); -} diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H index 606b4573e4..a243db6300 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H @@ -12,12 +12,12 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); } else // If not using kinetic theory is using Ct model { - nuEff1 = sqr(Ct)*nut2 + nu1; + nuEff1 = sqr(Ct)*nut2 + thermo1.mu()/rho1; } volTensorField Rc1 ( - "Rc1", + "Rc", (((2.0/3.0)*I)*nuEff1)*tr(gradU1T) - nuEff1*gradU1T ); @@ -52,7 +52,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime); volTensorField gradU2T(fvc::grad(U2)().T()); volTensorField Rc2 ( - "Rc2", + "Rc", (((2.0/3.0)*I)*nuEff2)*tr(gradU2T) - nuEff2*gradU2T ); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H index a818ee2e9e..86d9203dce 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H @@ -1,9 +1,9 @@ -surfaceScalarField alphaPhi1("alphaPhi1", phi1); -surfaceScalarField alphaPhi2("alphaPhi2", phi2); +surfaceScalarField alphaPhi1("alphaPhi", phi1); +surfaceScalarField alphaPhi2("alphaPhi", phi2); { - word scheme("div(phi,alpha1)"); - word schemer("div(phir,alpha1)"); + word scheme("div(phi,alpha)"); + word schemer("div(phir,alpha)"); surfaceScalarField phic("phic", phi); surfaceScalarField phir("phir", phi1 - phi2); @@ -13,7 +13,7 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2); surfaceScalarField alpha1f(fvc::interpolate(alpha1)); surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf()); phir += phipp; - phic += fvc::interpolate(alpha1)*phipp; + phic += alpha1f*phipp; } for (int acorr=0; acorr 0.0) { + surfaceScalarField alpha1f(fvc::interpolate(alpha1)); + ppMagf = fvc::interpolate((1.0/rho1)*rAU1) - *fvc::interpolate - ( - g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax) - ); + *g0*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax); + + // ppMagf = + // fvc::interpolate((1.0/rho1)*rAU1) + // *fvc::interpolate + // ( + // g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax) + // ); alpha1Eqn -= fvm::laplacian ( - (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf, + alpha1f*ppMagf, alpha1, "laplacian(alphaPpMag,alpha1)" ); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C index 855a862b6b..17d569c759 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C @@ -31,6 +31,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" +#include "rhoThermo.H" #include "nearWallDist.H" #include "wallFvPatch.H" #include "fixedValueFvsPatchFields.H" @@ -86,7 +87,7 @@ int main(int argc, char *argv[]) #include "alphaEqn.H" #include "kEpsilon.H" #include "interfacialCoeffs.H" - #include "TEqns.H" + #include "EEqns.H" #include "UEqns.H" // --- Pressure corrector loop diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H index e0cd5328c9..6095ad6845 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H @@ -12,55 +12,43 @@ ) ); + word phase1Name + ( + transportProperties.found("phases") + ? wordList(transportProperties.lookup("phases"))[0] + : "phase1" + ); + + word phase2Name + ( + transportProperties.found("phases") + ? wordList(transportProperties.lookup("phases"))[1] + : "phase2" + ); + autoPtr phase1 = phaseModel::New ( mesh, transportProperties, - "1" + phase1Name ); autoPtr phase2 = phaseModel::New ( mesh, transportProperties, - "2" + phase2Name ); + volScalarField& alpha1 = phase1(); + volScalarField& alpha2 = phase2(); + alpha2 = scalar(1) - alpha1; + volVectorField& U1 = phase1->U(); surfaceScalarField& phi1 = phase1->phi(); - const dimensionedScalar& nu1 = phase1->nu(); - const dimensionedScalar& k1 = phase1->kappa(); - const dimensionedScalar& Cp1 = phase1->Cp(); - dimensionedScalar rho10 - ( - "rho0", - dimDensity, - transportProperties.subDict("phase1").lookup("rho0") - ); - dimensionedScalar R1 - ( - "R", - dimensionSet(0, 2, -2, -1, 0), - transportProperties.subDict("phase1").lookup("R") - ); volVectorField& U2 = phase2->U(); surfaceScalarField& phi2 = phase2->phi(); - const dimensionedScalar& nu2 = phase2->nu(); - const dimensionedScalar& k2 = phase2->kappa(); - const dimensionedScalar& Cp2 = phase2->Cp(); - dimensionedScalar rho20 - ( - "rho0", - dimDensity, - transportProperties.subDict("phase2").lookup("rho0") - ); - dimensionedScalar R2 - ( - "R", - dimensionSet(0, 2, -2, -1, 0), - transportProperties.subDict("phase2").lookup("R") - ); dimensionedScalar pMin ( @@ -69,102 +57,16 @@ transportProperties.lookup("pMin") ); + rhoThermo& thermo1 = phase1->thermo(); + rhoThermo& thermo2 = phase2->thermo(); - Info<< "Reading field T1" << endl; - volScalarField T1 - ( - IOobject - ( - "T1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); + volScalarField& p = thermo1.p(); - Info<< "Reading field T2" << endl; - volScalarField T2 - ( - IOobject - ( - "T2", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); + volScalarField rho1("rho" + phase1Name, thermo1.rho()); + const volScalarField& psi1 = thermo1.psi(); - volScalarField psi1 - ( - IOobject - ( - "psi1", - runTime.timeName(), - mesh - ), - 1.0/(R1*T1) - ); - - volScalarField psi2 - ( - IOobject - ( - "psi2", - runTime.timeName(), - mesh - ), - 1.0/(R2*T2) - ); - psi2.oldTime(); - - Info<< "Reading field p\n" << endl; - volScalarField p - ( - IOobject - ( - "p", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - volScalarField rho1("rho1", rho10 + psi1*p); - volScalarField rho2("rho2", rho20 + psi2*p); - - - Info<< "Reading field alpha1\n" << endl; - volScalarField alpha1 - ( - IOobject - ( - "alpha1", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - volScalarField alpha2 - ( - IOobject - ( - "alpha2", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - scalar(1) - alpha1 - ); + volScalarField rho2("rho" + phase2Name, thermo2.rho()); + const volScalarField& psi2 = thermo2.psi(); volVectorField U ( @@ -179,27 +81,6 @@ alpha1*U1 + alpha2*U2 ); - dimensionedScalar Cvm - ( - "Cvm", - dimless, - transportProperties.lookup("Cvm") - ); - - dimensionedScalar Cl - ( - "Cl", - dimless, - transportProperties.lookup("Cl") - ); - - dimensionedScalar Ct - ( - "Ct", - dimless, - transportProperties.lookup("Ct") - ); - surfaceScalarField phi ( IOobject @@ -226,8 +107,6 @@ alpha1*rho1 + alpha2*rho2 ); - #include "createRASTurbulence.H" - Info<< "Calculating field DDtU1 and DDtU2\n" << endl; volVectorField DDtU1 @@ -248,6 +127,29 @@ Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); + dimensionedScalar Cvm + ( + "Cvm", + dimless, + transportProperties.lookup("Cvm") + ); + + dimensionedScalar Cl + ( + "Cl", + dimless, + transportProperties.lookup("Cl") + ); + + dimensionedScalar Ct + ( + "Ct", + dimless, + transportProperties.lookup("Ct") + ); + + #include "createRASTurbulence.H" + IOdictionary interfacialProperties ( IOobject @@ -297,8 +199,8 @@ if ( !( - dispersedPhase == "1" - || dispersedPhase == "2" + dispersedPhase == phase1Name + || dispersedPhase == phase2Name || dispersedPhase == "both" ) ) @@ -337,7 +239,7 @@ ( IOobject ( - "rAU1", + "rAU" + phase1Name, runTime.timeName(), mesh, IOobject::NO_READ, @@ -387,5 +289,5 @@ ); Info<< "Creating field kinetic energy K\n" << endl; - volScalarField K1("K1", 0.5*magSqr(U1)); - volScalarField K2("K2", 0.5*magSqr(U2)); + volScalarField K1("K" + phase1Name, 0.5*magSqr(U1)); + volScalarField K2("K" + phase2Name, 0.5*magSqr(U2)); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H index 8f05a1d1b3..0a782ef51e 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H @@ -151,7 +151,7 @@ ( IOobject ( - "nut2", + "nut" + phase2Name, runTime.timeName(), mesh, IOobject::NO_READ, @@ -165,13 +165,13 @@ ( IOobject ( - "nuEff1", + "nuEff" + phase1Name, runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), - sqr(Ct)*nut2 + nu1 + sqr(Ct)*nut2 + thermo1.mu()/rho1 ); Info<< "Calculating field nuEff2\n" << endl; @@ -179,11 +179,11 @@ ( IOobject ( - "nuEff2", + "nuEff" + phase2Name, runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), - nut2 + nu2 + nut2 + thermo2.mu()/rho2 ); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H index d6ccf90289..d53bec5ea4 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H @@ -38,12 +38,12 @@ volScalarField heatTransferCoeff volVectorField Ur(U1 - U2); volScalarField magUr(mag(Ur) + residualSlip); - if (dispersedPhase == "1") + if (dispersedPhase == phase1Name) { dragCoeff = drag1->K(magUr); heatTransferCoeff = heatTransfer1->K(magUr); } - else if (dispersedPhase == "2") + else if (dispersedPhase == phase2Name) { dragCoeff = drag2->K(magUr); heatTransferCoeff = heatTransfer2->K(magUr); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options index 92fe0b0d74..f031e05898 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options @@ -1,7 +1,9 @@ EXE_INC = \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I../phaseModel/lnInclude LIB_LIBS = \ - -lcompressiblePhaseModel - + -lcompressiblePhaseModel \ + -lfluidThermophysicalModels \ + -lspecie diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C index 9bc6cadaf6..dc41c062df 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C @@ -72,8 +72,7 @@ Foam::tmp Foam::heatTransferModels::RanzMarshall::K ) const { volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3))); - dimensionedScalar Prb = - phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa(); + volScalarField Prb(phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa()); volScalarField Nu(scalar(2) + 0.6*sqrt(Re)*cbrt(Prb)); return 6.0*phase2_.kappa()*Nu/sqr(phase1_.d()); diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options index 2fcce9913d..7fdabf5730 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options @@ -1,5 +1,6 @@ EXE_INC = \ -I$(LIB_SRC)/foam/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I../phaseModel/lnInclude \ -I../interfacialModels/lnInclude diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C index 7d0767914e..f8185f5272 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C @@ -45,8 +45,8 @@ Foam::kineticTheoryModel::kineticTheoryModel phi1_(phase1.phi()), draga_(draga), - rho1_(phase1.rho()), - nu1_(phase1.nu()), + rho1_(phase1.rho()[0]), //***HGW + nu1_(phase1.nu()()[0]), //***HGW kineticTheoryProperties_ ( @@ -120,7 +120,7 @@ Foam::kineticTheoryModel::kineticTheoryModel ( IOobject ( - "mu1", + "mu" + phase1.name(), U1_.time().timeName(), U1_.mesh(), IOobject::NO_READ, diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H index f7010b1c45..f9a7807d00 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H @@ -1,6 +1,6 @@ { - rho1 = rho10 + psi1*p; - rho2 = rho20 + psi2*p; + rho1 = thermo1.rho(); + rho2 = thermo2.rho(); surfaceScalarField alpha1f(fvc::interpolate(alpha1)); surfaceScalarField alpha2f(scalar(1) - alpha1f); @@ -11,10 +11,10 @@ surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rAU1)); surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rAU2)); - volVectorField HbyA1("HbyA1", U1); + volVectorField HbyA1("HbyA" + phase1Name, U1); HbyA1 = rAU1*U1Eqn.H(); - volVectorField HbyA2("HbyA2", U2); + volVectorField HbyA2("HbyA" + phase2Name, U2); HbyA2 = rAU2*U2Eqn.H(); mrfZones.absoluteFlux(phi1.oldTime()); @@ -38,14 +38,14 @@ surfaceScalarField phiHbyA1 ( - "phiHbyA1", + "phiHbyA" + phase1Name, (fvc::interpolate(HbyA1) & mesh.Sf()) + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1) ); surfaceScalarField phiHbyA2 ( - "phiHbyA2", + "phiHbyA" + phase2Name, (fvc::interpolate(HbyA2) & mesh.Sf()) + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2) ); @@ -96,8 +96,16 @@ //} //else { - surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi1); - surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi2); + surfaceScalarField phid1 + ( + "phid" + phase1Name, + fvc::interpolate(psi1)*phi1 + ); + surfaceScalarField phid2 + ( + "phid" + phase2Name, + fvc::interpolate(psi2)*phi2 + ); pEqnComp1 = fvc::ddt(rho1) + psi1*correction(fvm::ddt(p)) @@ -174,13 +182,15 @@ p = max(p, pMin); - rho1 = rho10 + psi1*p; - rho2 = rho20 + psi2*p; + thermo1.correct(); + thermo2.correct(); + rho1 = thermo1.rho(); + rho2 = thermo2.rho(); K1 = 0.5*magSqr(U1); K2 = 0.5*magSqr(U2); - //***HGW if (thermo.dpdt()) + if (thermo1.dpdt()) { dpdt = fvc::ddt(p); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options index 0ec1139209..e441b0417b 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options @@ -1,6 +1,9 @@ EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude LIB_LIBS = \ - -lincompressibleTransportModels + -lincompressibleTransportModels \ + -lfluidThermophysicalModels \ + -lspecie diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C index d80758ac8c..bafa3060eb 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C @@ -37,35 +37,22 @@ Foam::phaseModel::phaseModel const word& phaseName ) : - dict_ + volScalarField ( - transportProperties.subDict("phase" + phaseName) + IOobject + ( + "alpha" + phaseName, + mesh.time().timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("alpha", dimless, 0) ), name_(phaseName), - nu_ - ( - "nu", - dimensionSet(0, 2, -1, 0, 0), - dict_.lookup("nu") - ), - kappa_ - ( - "kappa", - dimensionSet(1, 1, -3, -1, 0), - dict_.lookup("kappa") - ), - Cp_ - ( - "Cp", - dimensionSet(0, 2, -2, -1, 0), - dict_.lookup("Cp") - ), - rho_ - ( - "rho", - dimDensity, - dict_.lookup("rho") - ), + phaseDict_(transportProperties.subDict(phaseName)), + thermo_(rhoThermo::New(mesh, phaseName)), U_ ( IOobject @@ -79,6 +66,8 @@ Foam::phaseModel::phaseModel mesh ) { + thermo_->validate("phaseModel " + phaseName, "h", "e"); + const word phiName = "phi" + phaseName; IOobject phiHeader @@ -147,7 +136,7 @@ Foam::phaseModel::phaseModel dPtr_ = diameterModel::New ( - dict_, + phaseDict_, *this ); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H index 90b2d2184e..a96779b8a2 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H @@ -36,6 +36,7 @@ SourceFiles #include "dimensionedScalar.H" #include "volFields.H" #include "surfaceFields.H" +#include "rhoThermo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -50,25 +51,18 @@ class diameterModel; \*---------------------------------------------------------------------------*/ class phaseModel +: + public volScalarField { // Private data - dictionary dict_; - //- Name of phase word name_; - //- Kinematic viscosity - dimensionedScalar nu_; + dictionary phaseDict_; - //- Thermal conductivity - dimensionedScalar kappa_; - - //- Heat capacity - dimensionedScalar Cp_; - - //- Density - dimensionedScalar rho_; + //- Thermophysical properties + autoPtr thermo_; //- Velocity volVectorField U_; @@ -116,24 +110,34 @@ public: tmp d() const; - const dimensionedScalar& nu() const + tmp nu() const { - return nu_; + return thermo_->mu()/thermo_->rho(); } - const dimensionedScalar& kappa() const + tmp kappa() const { - return kappa_; + return thermo_->kappa(); } - const dimensionedScalar& Cp() const + tmp Cp() const { - return Cp_; + return thermo_->Cp(); } - const dimensionedScalar& rho() const + const volScalarField& rho() const { - return rho_; + return thermo_->rho(); + } + + const rhoThermo& thermo() const + { + return thermo_(); + } + + rhoThermo& thermo() + { + return thermo_(); } const volVectorField& U() const diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H index 7f53411529..bc9a07b0a8 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H @@ -61,4 +61,4 @@ if (turbulence) #include "wallViscosity.H" } -nuEff2 = nut2 + nu2; +nuEff2 = nut2 + thermo2.mu()/rho2; diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H index c91cce9a00..d85181cba2 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H @@ -4,7 +4,6 @@ scalar Cmu25 = ::pow(Cmu.value(), 0.25); scalar Cmu75 = ::pow(Cmu.value(), 0.75); scalar kappa_ = kappa.value(); - scalar nu2_ = nu2.value(); const fvPatchList& patches = mesh.boundary(); @@ -30,6 +29,8 @@ forAll(patches, patchi) { const fvPatch& currPatch = patches[patchi]; + const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi]; + const scalarField& rho2_ = rho2.boundaryField()[patchi]; if (isA(currPatch)) { @@ -52,7 +53,8 @@ /(kappa_*y[patchi][facei]); G[faceCelli] += - (nut2w[facei] + nu2_)*magFaceGradU[facei] + (nut2w[facei] + mu2_[facei]/rho2_[facei]) + *magFaceGradU[facei] *Cmu25*::sqrt(k[faceCelli]) /(kappa_*y[patchi][facei]); } diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H index b153a36014..9aa032149c 100644 --- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H +++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H @@ -2,13 +2,14 @@ scalar Cmu25 = ::pow(Cmu.value(), 0.25); scalar kappa_ = kappa.value(); scalar E_ = E.value(); - scalar nu2_ = nu2.value(); const fvPatchList& patches = mesh.boundary(); forAll(patches, patchi) { const fvPatch& currPatch = patches[patchi]; + const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi]; + const scalarField& rho2_ = rho2.boundaryField()[patchi]; if (isA(currPatch)) { @@ -20,11 +21,14 @@ // calculate yPlus scalar yPlus = - Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nu2_; + Cmu25*y[patchi][facei]*::sqrt(k[faceCelli]) + /(mu2_[facei]/rho2_[facei]); if (yPlus > 11.6) { - nutw[facei] = nu2_*(yPlus*kappa_/::log(E_*yPlus) - 1); + nutw[facei] = + (mu2_[facei]/rho2_[facei]) + *(yPlus*kappa_/::log(E_*yPlus) - 1); } else {