From c957f2450d5d5b7f638781f8b8d9099163c55f16 Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 26 Mar 2014 12:48:37 +0000 Subject: [PATCH 1/4] continuousGasKEpsilon: Omega now consistent with Lahey 2005 paper --- .../continuousGasKEpsilon/continuousGasKEpsilon.C | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.C b/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.C index b6eef980..b3b453d6 100644 --- a/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.C +++ b/src/TurbulenceModels/phaseIncompressible/RAS/continuousGasKEpsilon/continuousGasKEpsilon.C @@ -124,9 +124,17 @@ void continuousGasKEpsilon::correctNut() const transportModel& liquid = fluid.otherPhase(gas); volScalarField thetal(liquidTurbulence.k()/liquidTurbulence.epsilon()); - volScalarField thetag((1.0/(18*liquid.nu()))*sqr(gas.d())); - volScalarField expThetar(exp(min(thetal/thetag, scalar(50)))); - volScalarField omega(sqr(expThetar - 1)/(sqr(expThetar) - 1)); + volScalarField rhodv(gas.rho() + fluid.virtualMass(gas).Cvm()*liquid.rho()); + volScalarField thetag((rhodv/(18*liquid.rho()*liquid.nu()))*sqr(gas.d())); + volScalarField expThetar + ( + min + ( + exp(min(thetal/thetag, scalar(50))), + scalar(1) + ) + ); + volScalarField omega((1 - expThetar)/(1 + expThetar)); nutEff_ = omega*liquidTurbulence.nut(); } From 962bb79676c3a56cd1620ac79c705fade3bc423d Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 26 Mar 2014 12:49:18 +0000 Subject: [PATCH 2/4] driftFluxFoam: Changed the laminar viscosity used in the k-epsilon model that of the mixture --- .../multiphase/driftFluxFoam/driftFluxFoam.C | 1 - .../multiphase/driftFluxFoam/kEpsilon.H | 20 ++++++++++--------- .../multiphase/driftFluxFoam/wallFunctions.H | 4 ++-- .../multiphase/driftFluxFoam/wallViscosity.H | 6 +++--- 4 files changed, 16 insertions(+), 15 deletions(-) diff --git a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C index da03517e..8ab48d7b 100644 --- a/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C +++ b/applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C @@ -87,7 +87,6 @@ int main(int argc, char *argv[]) #include "alphaEqnSubCycle.H" twoPhaseProperties.correct(); - Info<< average(twoPhaseProperties.mu()) << endl; #include "UEqn.H" diff --git a/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H b/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H index 66216b87..ddcc8cf2 100644 --- a/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H +++ b/applications/solvers/multiphase/driftFluxFoam/kEpsilon.H @@ -10,8 +10,6 @@ if (turbulence) dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), 0); dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL); - volScalarField divU(fvc::div(phi)); - tmp tgradU = fvc::grad(U); volScalarField G(mut*(tgradU() && dev(twoSymm(tgradU())))); tgradU.clear(); @@ -21,7 +19,7 @@ if (turbulence) Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin) ); - volScalarField muc(twoPhaseProperties.nucModel().nu()*rho2); + volScalarField mul(twoPhaseProperties.mu()); #include "wallFunctions.H" @@ -32,12 +30,12 @@ if (turbulence) + fvm::div(rhoPhi, epsilon) - fvm::laplacian ( - mut/sigmaEps + muc, epsilon, + mut/sigmaEps + mul, epsilon, "laplacian(DepsilonEff,epsilon)" ) == C1*G*epsilon/(k + kMin) - - fvm::SuSp(C1*(1.0 - C3)*Gcoef + (2.0/3.0*C1)*rho*divU, epsilon) + - fvm::SuSp(C1*(1.0 - C3)*Gcoef, epsilon) - fvm::Sp(C2*rho*epsilon/(k + kMin), epsilon) ); @@ -56,12 +54,12 @@ if (turbulence) + fvm::div(rhoPhi, k) - fvm::laplacian ( - mut/sigmak + muc, k, + mut/sigmak + mul, k, "laplacian(DkEff,k)" ) == G - - fvm::SuSp(Gcoef + 2.0/3.0*rho*divU, k) + - fvm::SuSp(Gcoef, k) - fvm::Sp(rho*epsilon/(k + kMin), k) ); @@ -75,6 +73,10 @@ if (turbulence) mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin); #include "wallViscosity.H" -} -muEff = mut + twoPhaseProperties.mu(); + muEff = mut + mul; +} +else +{ + muEff = mut + twoPhaseProperties.mu(); +} diff --git a/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H b/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H index b9ff8481..2a064b47 100644 --- a/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H +++ b/applications/solvers/multiphase/driftFluxFoam/wallFunctions.H @@ -33,7 +33,7 @@ if (isA(curPatch)) { const scalarField& mutw = mut.boundaryField()[patchi]; - const scalarField& mucw = muc.boundaryField()[patchi]; + const scalarField& mulw = mul.boundaryField()[patchi]; scalarField magFaceGradU ( @@ -55,7 +55,7 @@ /(kappa_*y[patchi][facei]); G[faceCelli] += - (mutw[facei] + mucw[facei]) + (mutw[facei] + mulw[facei]) *magFaceGradU[facei] *Cmu25*::sqrt(k[faceCelli]) /(kappa_*y[patchi][facei]); diff --git a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H index 2b54f6c2..d73e0a5d 100644 --- a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H +++ b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H @@ -12,7 +12,7 @@ if (isA(curPatch)) { scalarField& mutw = mut.boundaryField()[patchi]; - const scalarField& mucw = muc.boundaryField()[patchi]; + const scalarField& mulw = mul.boundaryField()[patchi]; forAll(curPatch, facei) { @@ -20,12 +20,12 @@ scalar yPlus = Cmu25*y[patchi][facei]*::sqrt(k[faceCelli]) - /(mucw[facei]/rho2.value()); + /(mulw[facei]/rho2.value()); if (yPlus > 11.6) { mutw[facei] = - mucw[facei]*(yPlus*kappa_/::log(E_*yPlus) - 1); + mulw[facei]*(yPlus*kappa_/::log(E_*yPlus) - 1); } else { From 840078c02e073350059dad00fe2f8638174e2245 Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 26 Mar 2014 14:31:09 +0000 Subject: [PATCH 3/4] driftFluxFoam: corrected density to get kinematic viscosity for yPlus --- applications/solvers/multiphase/driftFluxFoam/wallViscosity.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H index d73e0a5d..b75db346 100644 --- a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H +++ b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H @@ -20,7 +20,7 @@ scalar yPlus = Cmu25*y[patchi][facei]*::sqrt(k[faceCelli]) - /(mulw[facei]/rho2.value()); + /(mulw[facei]/rho[facei]); if (yPlus > 11.6) { From 43841790346c8f52a0188899eb6189b671a08b8e Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 26 Mar 2014 15:19:39 +0000 Subject: [PATCH 4/4] Corrected rho used for yPlus --- applications/solvers/multiphase/driftFluxFoam/wallViscosity.H | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H index b75db346..632432ff 100644 --- a/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H +++ b/applications/solvers/multiphase/driftFluxFoam/wallViscosity.H @@ -13,6 +13,7 @@ { scalarField& mutw = mut.boundaryField()[patchi]; const scalarField& mulw = mul.boundaryField()[patchi]; + const scalarField& rhow = rho.boundaryField()[patchi]; forAll(curPatch, facei) { @@ -20,7 +21,7 @@ scalar yPlus = Cmu25*y[patchi][facei]*::sqrt(k[faceCelli]) - /(mulw[facei]/rho[facei]); + /(mulw[facei]/rhow[facei]); if (yPlus > 11.6) {