From 0696e9aa441f29defef1a9ec5a8354c4d22b0b9c Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Wed, 20 May 2009 14:50:08 +0100
Subject: [PATCH] corrections to chemical enthalpy transfer

---
 .../ReactingMultiphaseParcel.C                     |  7 +++++--
 .../ReactingMultiphaseParcelI.H                    | 14 +++++++++++++-
 2 files changed, 18 insertions(+), 3 deletions(-)

diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 7c2fd37583..e5d8e25d5f 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -174,9 +174,12 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::cellValueSourceCorrection
     const scalar cpc = td.cpInterp().psi()[cellI];
     this->cpc_ = (massCell*cpc + addedMass*cpEff)/massCellNew;
 
-    const scalar fCarrier = -1.0/td.constProps().hRetentionCoeff();
+    const scalar fCarrier =
+        (1.0 - td.cloud().constProps().hRetentionCoeff())
+       /td.cloud().constProps().hRetentionCoeff();
+
     const scalar dh =
-        td.cloud().hsTrans()[cellI] + fCarrier*td.cloud().hcTrans()[cellI];
+        td.cloud().hsTrans()[cellI] - fCarrier*td.cloud().hcTrans()[cellI];
 
     this->Tc_ += dh/(this->cpc_*massCellNew);
 }
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
index ef40cdc994..e70b922ae2 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H
@@ -39,7 +39,19 @@ constantProperties
     (
         dimensionedScalar(this->dict().lookup("hRetentionCoeff")).value()
     )
-{}
+{
+    if ((hRetentionCoeff_ < 0) || (hRetentionCoeff_ > 1))
+    {
+        FatalErrorIn
+        (
+            "ReactingMultiphaseParcel<ParcelType>::constantProperties::"
+            "constantProperties"
+        )   << "hRetentionCoeff must be in the range 0 to 1" << nl
+            << exit(FatalError) << endl;
+    }
+
+    hRetentionCoeff_ = max(1e-06, hRetentionCoeff_);
+}
 
 
 template<class ParcelType>