From 994761a422b4724f356b584b6a82957a97eed735 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Wed, 12 Jan 2011 13:17:39 +0000
Subject: [PATCH] ENH: Updated cell value corrections for reacting particles

---
 .../ReactingMultiphaseParcel.C                | 44 +------------------
 .../Templates/ReactingParcel/ReactingParcel.C | 25 +++++++++--
 2 files changed, 23 insertions(+), 46 deletions(-)

diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
index 98e357950d..e610f956b6 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C
@@ -151,48 +151,8 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::cellValueSourceCorrection
     const label cellI
 )
 {
-    scalar massCell = this->massCell(cellI);
-
-    scalar addedMass = 0.0;
-    forAll(td.cloud().rhoTrans(), i)
-    {
-        addedMass += td.cloud().rhoTrans(i)[cellI];
-    }
-
-    this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
-
-    scalar massCellNew = massCell + addedMass;
-    this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
-
-    scalar CpEff = 0;
-    if (addedMass > ROOTVSMALL)
-    {
-        forAll(td.cloud().rhoTrans(), i)
-        {
-            scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
-            CpEff += Y*td.cloud().thermo().carrier().Cp(i, this->Tc_);
-        }
-    }
-    const scalar Cpc = td.CpInterp().psi()[cellI];
-    this->Cpc_ = (massCell*Cpc + addedMass*CpEff)/massCellNew;
-
-    this->Tc_ += td.cloud().hsTrans()[cellI]/(this->Cpc_*massCellNew);
-
-    if (this->Tc_ < td.cloud().constProps().TMin())
-    {
-        WarningIn
-        (
-            "void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection"
-            "("
-                "TrackData&, "
-                "const scalar, "
-                "const label"
-            ")"
-        )   << "Limiting observed temperature in cell " << cellI << " to "
-            << td.cloud().constProps().TMin() <<  nl << endl;
-
-        this->Tc_ = td.cloud().constProps().TMin();
-    }
+    // Re-use correction from reacting parcel
+    ReactingParcel<ParcelType>::cellValueSourceCorrection(td, dt, cellI);
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
index 2c05144d9d..b5199b3465 100644
--- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C
@@ -75,28 +75,43 @@ void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
     const label cellI
 )
 {
-    scalar massCell = this->massCell(cellI);
-
     scalar addedMass = 0.0;
     forAll(td.cloud().rhoTrans(), i)
     {
         addedMass += td.cloud().rhoTrans(i)[cellI];
     }
 
-    this->rhoc_ += addedMass/td.cloud().pMesh().cellVolumes()[cellI];
+    if (addedMass < ROOTVSMALL)
+    {
+        return;
+    }
+
+    scalar massCell = this->massCell(cellI);
+
+    const scalar V = td.cloud().pMesh().cellVolumes()[cellI];
+    this->rhoc_ += addedMass/V;
 
     scalar massCellNew = massCell + addedMass;
     this->Uc_ += td.cloud().UTrans()[cellI]/massCellNew;
 
-    scalar CpEff = 0;
+    scalar CpEff = 0.0;
+    scalar CsEff = 0.0;
+    scalar Csc = 0.0;
     if (addedMass > ROOTVSMALL)
     {
         forAll(td.cloud().rhoTrans(), i)
         {
             scalar Y = td.cloud().rhoTrans(i)[cellI]/addedMass;
             CpEff += Y*td.cloud().composition().carrier().Cp(i, this->Tc_);
+
+            scalar W = td.cloud().composition().carrier().W(i);
+            CsEff += td.cloud().rhoTrans(i)[cellI]/V*Y/W;
+            scalar Yc = td.cloud().composition().carrier().Y(i)[cellI];
+            Csc += massCell/V*Yc/W;
         }
     }
+    CsEff = (massCell*Csc + addedMass*CsEff)/massCellNew;
+
     const scalar Cpc = td.CpInterp().psi()[cellI];
     this->Cpc_ = (massCell*Cpc + addedMass*CpEff)/massCellNew;
 
@@ -117,6 +132,8 @@ void Foam::ReactingParcel<ParcelType>::cellValueSourceCorrection
 
         this->Tc_ = td.cloud().constProps().TMin();
     }
+
+    this->pc_ = CsEff*specie::RR*this->Tc_;
 }