From c957f2450d5d5b7f638781f8b8d9099163c55f16 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 26 Mar 2014 12:48:37 +0000
Subject: [PATCH] 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<BasicTurbulenceModel>::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();
 }