From 3087fcbb1e6bf82956f34f741bd7a8e76ca449f2 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 9 Feb 2015 21:25:36 +0000
Subject: [PATCH] constTransport: Avoid /0 in + and - operators averaging 1/Pr
 when the number of moles = 0 Resolves bug-report
 http://www.openfoam.org/mantisbt/view.php?id=1348

---
 .../specie/transport/const/constTransportI.H  | 56 +++++++++++++------
 1 file changed, 40 insertions(+), 16 deletions(-)

diff --git a/src/thermophysicalModels/specie/transport/const/constTransportI.H b/src/thermophysicalModels/specie/transport/const/constTransportI.H
index c0c15372..050d6d79 100644
--- a/src/thermophysicalModels/specie/transport/const/constTransportI.H
+++ b/src/thermophysicalModels/specie/transport/const/constTransportI.H
@@ -203,15 +203,27 @@ inline Foam::constTransport<Thermo> Foam::operator+
         static_cast<const Thermo&>(ct1) + static_cast<const Thermo&>(ct2)
     );
 
-    scalar molr1 = ct1.nMoles()/t.nMoles();
-    scalar molr2 = ct2.nMoles()/t.nMoles();
+    if (mag(ct1.nMoles()) + mag(ct2.nMoles()) < SMALL)
+    {
+        return constTransport<Thermo>
+        (
+            t,
+            0,
+            ct1.rPr_
+        );
+    }
+    else
+    {
+        scalar molr1 = ct1.nMoles()/t.nMoles();
+        scalar molr2 = ct2.nMoles()/t.nMoles();
 
-    return constTransport<Thermo>
-    (
-        t,
-        molr1*ct1.mu_ + molr2*ct2.mu_,
-        1.0/(molr1/ct1.rPr_ + molr2/ct2.rPr_)
-    );
+        return constTransport<Thermo>
+        (
+            t,
+            molr1*ct1.mu_ + molr2*ct2.mu_,
+            1.0/(molr1/ct1.rPr_ + molr2/ct2.rPr_)
+        );
+    }
 }
 
 
@@ -227,15 +239,27 @@ inline Foam::constTransport<Thermo> Foam::operator-
         static_cast<const Thermo&>(ct1) - static_cast<const Thermo&>(ct2)
     );
 
-    scalar molr1 = ct1.nMoles()/t.nMoles();
-    scalar molr2 = ct2.nMoles()/t.nMoles();
+    if (mag(ct1.nMoles()) + mag(ct2.nMoles()) < SMALL)
+    {
+        return constTransport<Thermo>
+        (
+            t,
+            0,
+            ct1.rPr_
+        );
+    }
+    else
+    {
+        scalar molr1 = ct1.nMoles()/t.nMoles();
+        scalar molr2 = ct2.nMoles()/t.nMoles();
 
-    return constTransport<Thermo>
-    (
-        t,
-        molr1*ct1.mu_ - molr2*ct2.mu_,
-        1.0/(molr1/ct1.rPr_ - molr2/ct2.rPr_)
-    );
+        return constTransport<Thermo>
+        (
+            t,
+            molr1*ct1.mu_ - molr2*ct2.mu_,
+            1.0/(molr1/ct1.rPr_ - molr2/ct2.rPr_)
+        );
+    }
 }