From d7fcc07675e04359d5678a186678747c2692bbcd Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Thu, 29 Dec 2011 14:01:02 +0000
Subject: [PATCH] compressibleTwoPhaseEulerFoam, multiphaseEulerFoam: correct
 drag coeff boundary field

---
 .../GidaspowErgunWenYu/GidaspowErgunWenYu.C   | 41 +++++++------------
 .../GidaspowSchillerNaumann.C                 | 14 +++----
 .../SchillerNaumann/SchillerNaumann.C         | 14 +++----
 .../dragModels/SyamlalOBrien/SyamlalOBrien.C  | 14 +++----
 .../dragModels/WenYu/WenYu.C                  | 14 +++----
 .../Schaeffer/SchaefferFrictionalStress.C     |  2 +
 .../GidaspowErgunWenYu/GidaspowErgunWenYu.C   | 41 +++++++------------
 .../GidaspowSchillerNaumann.C                 | 14 +++----
 .../SchillerNaumann/SchillerNaumann.C         | 14 +++----
 .../dragModels/SyamlalOBrien/SyamlalOBrien.C  | 14 +++----
 .../dragModels/WenYu/WenYu.C                  | 14 +++----
 .../Schaeffer/SchaefferFrictionalStress.C     |  2 +
 12 files changed, 74 insertions(+), 124 deletions(-)

diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
index bac588c7389..3f659d44ec4 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
@@ -76,34 +76,23 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
     volScalarField bp(pow(beta, -2.65));
     volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3)));
 
-    volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     // Wen and Yu (1966)
-    tmp<volScalarField> tKWenYu(0.75*Cds*phase2_.rho()*Ur*bp/d);
-    volScalarField& KWenYu = tKWenYu();
-
-    // Ergun
-    forAll (beta, cellj)
-    {
-        if (beta[cellj] <= 0.8)
-        {
-            KWenYu[cellj] =
-                150.0*alpha_[cellj]*phase2_.nu().value()*phase2_.rho().value()
-               /sqr(beta[cellj]*d[cellj])
-              + 1.75*phase2_.rho().value()*Ur[cellj]
-               /(beta[cellj]*d[cellj]);
-        }
-    }
-
-    return tKWenYu;
+    return
+    (
+        pos(beta - 0.8)
+       *(0.75*Cds*phase2_.rho()*Ur*bp/d)
+      + neg(beta - 0.8)
+       *(
+           150.0*alpha_*phase2_.nu()*phase2_.rho()/(sqr(beta*d))
+         + 1.75*phase2_.rho()*Ur/(beta*d)
+        )
+    );
 }
 
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
index 99b5cacdf2f..0a624c6b0a2 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
@@ -75,15 +75,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
index 79259cbe058..f68e5379d24 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
@@ -72,15 +72,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::K
 ) const
 {
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
index 4f971b4bb76..541868cf870 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
@@ -73,15 +73,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
 {
     volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
     volScalarField A(pow(beta, 4.14));
-    volScalarField B(0.8*pow(beta, 1.28));
-
-    forAll (beta, celli)
-    {
-        if (beta[celli] > 0.85)
-        {
-            B[celli] = pow(beta[celli], 2.65);
-        }
-    }
+    volScalarField B
+    (
+        neg(beta - 0.85)*(0.8*pow(beta, 1.28))
+      + pos(beta - 0.85)*(pow(beta, 2.65))
+    );
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
 
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
index f1126125bb1..bf279a3990a 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
@@ -75,15 +75,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
index 4320bbad354..d458f07b754 100644
--- a/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+++ b/applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
@@ -149,6 +149,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
         }
     }
 
+    muff.correctBoundaryConditions();
+
     return tmuf;
 }
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
index e34f219ebc0..b52b766de92 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowErgunWenYu/GidaspowErgunWenYu.C
@@ -75,34 +75,23 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowErgunWenYu::K
     volScalarField bp(pow(beta, -2.65));
     volScalarField Re(max(Ur*d/phase2_.nu(), scalar(1.0e-3)));
 
-    volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     // Wen and Yu (1966)
-    tmp<volScalarField> tKWenYu = 0.75*Cds*phase2_.rho()*Ur*bp/d;
-    volScalarField& KWenYu = tKWenYu();
-
-    // Ergun
-    forAll (beta, cellj)
-    {
-        if (beta[cellj] <= 0.8)
-        {
-            KWenYu[cellj] =
-                150.0*phase1_[cellj]*phase2_.nu().value()*phase2_.rho().value()
-               /sqr(beta[cellj]*d[cellj])
-              + 1.75*phase2_.rho().value()*Ur[cellj]
-               /(beta[cellj]*d[cellj]);
-        }
-    }
-
-    return tKWenYu;
+    return
+    (
+        pos(beta - 0.8)
+       *(0.75*Cds*phase2_.rho()*Ur*bp/d)
+      + neg(beta - 0.8)
+       *(
+           150.0*phase1_*phase2_.nu()*phase2_.rho()/(sqr(beta*d))
+         + 1.75*phase2_.rho()*Ur/(beta*d)
+        )
+    );
 }
 
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
index ecd25eb386f..25de6b1dc70 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/GidaspowSchillerNaumann/GidaspowSchillerNaumann.C
@@ -74,15 +74,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::GidaspowSchillerNaumann::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(beta*Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
index 3259b5d0e08..d4d77a2bb35 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SchillerNaumann/SchillerNaumann.C
@@ -71,15 +71,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::K
 ) const
 {
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
index 736be885e0f..07f85b38a1f 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/SyamlalOBrien/SyamlalOBrien.C
@@ -72,15 +72,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::SyamlalOBrien::K
 {
     volScalarField beta(max(phase2_, scalar(1.0e-6)));
     volScalarField A(pow(beta, 4.14));
-    volScalarField B(0.8*pow(beta, 1.28));
-
-    forAll (beta, celli)
-    {
-        if (beta[celli] > 0.85)
-        {
-            B[celli] = pow(beta[celli], 2.65);
-        }
-    }
+    volScalarField B
+    (
+        neg(beta - 0.85)*(0.8*pow(beta, 1.28))
+      + pos(beta - 0.85)*(pow(beta, 2.65))
+    );
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
index 6cfc119c1c4..6319276408b 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/WenYu/WenYu.C
@@ -74,15 +74,11 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::WenYu::K
     volScalarField bp(pow(beta, -2.65));
 
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
-
-    forAll(Re, celli)
-    {
-        if (Re[celli] > 1000.0)
-        {
-            Cds[celli] = 0.44;
-        }
-    }
+    volScalarField Cds
+    (
+        neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+      + pos(Re - 1000)*0.44
+    );
 
     return 0.75*Cds*phase2_.rho()*Ur*bp/phase1_.d();
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C b/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
index 4320bbad354..d458f07b754 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
@@ -149,6 +149,8 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
         }
     }
 
+    muff.correctBoundaryConditions();
+
     return tmuf;
 }
 
-- 
GitLab