From f17071cf7040ad7250561e7d435c2ad27e252d5e Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Fri, 13 Feb 2015 12:41:34 +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  | 58 +++++++++++++------
 1 file changed, 41 insertions(+), 17 deletions(-)

diff --git a/src/thermophysicalModels/specie/transport/const/constTransportI.H b/src/thermophysicalModels/specie/transport/const/constTransportI.H
index c0c15372eb..669623abab 100644
--- a/src/thermophysicalModels/specie/transport/const/constTransportI.H
+++ b/src/thermophysicalModels/specie/transport/const/constTransportI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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_)
+        );
+    }
 }