Lagrangian radiation models and fvDOM: Apply updates

Applying patches provided by Timo Niemi
Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1636

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C

@@ -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
@@ -210,6 +210,9 @@ void Foam::ThermoParcel<ParcelType>::calc
     // Sum Ni*Cpi*Wi of emission species
     scalar NCpW = 0.0;
 
+    // Store T for consistent radiation source
+    const scalar T0 = this->T_;
+
     // Calculate new particle temperature
     this->T_ =
         this->calcHeatTransfer
@@ -255,7 +258,7 @@ void Foam::ThermoParcel<ParcelType>::calc
         if (td.cloud().radiation())
         {
             const scalar ap = this->areaP();
-            const scalar T4 = pow4(this->T_);
+            const scalar T4 = pow4(T0);
             td.cloud().radAreaP()[cellI] += dt*np0*ap;
             td.cloud().radT4()[cellI] += dt*np0*T4;
             td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
@@ -304,7 +307,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     htc = max(htc, ROOTVSMALL);
     const scalar As = this->areaS(d);
 
-    scalar ap = Tc_ + Sh/As/htc;
+    scalar ap = Tc_ + Sh/(As*htc);
     scalar bp = 6.0*(Sh/As + htc*(Tc_ - T_));
     if (td.cloud().radiation())
     {
@@ -313,8 +316,11 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
         const scalar sigma = physicoChemical::sigma.value();
         const scalar epsilon = td.cloud().constProps().epsilon0();
 
-        ap = (ap + epsilon*Gc/(4.0*htc))/(1.0 + epsilon*sigma*pow3(T_)/htc);
-        bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T_)));
+        // Assume constant source
+        scalar s = epsilon*(Gc/4.0 - sigma*pow4(T_));
+
+        ap += s/htc;
+        bp += 6.0*s;
     }
     bp /= rho*d*Cp_*(ap - T_) + ROOTVSMALL;
 

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C

@@ -482,7 +482,8 @@ Foam::tmp<Foam::volScalarField> Foam::radiation::fvDOM::Rp() const
                 IOobject::NO_WRITE,
                 false
             ),
-            4.0*a_*physicoChemical::sigma //absorptionEmission_->a()
+            // Only include continuous phase emission
+            4*absorptionEmission_->aCont()*physicoChemical::sigma
         )
     );
 }
@@ -494,12 +495,15 @@ Foam::radiation::fvDOM::Ru() const
 
     const DimensionedField<scalar, volMesh>& G =
         G_.dimensionedInternalField();
+
     const DimensionedField<scalar, volMesh> E =
         absorptionEmission_->ECont()().dimensionedInternalField();
-    const DimensionedField<scalar, volMesh> a =
-        a_.dimensionedInternalField();
 
-    return  a*G - E;
+    // Only include continuous phase absorption
+    const DimensionedField<scalar, volMesh> a =
+        absorptionEmission_->aCont()().dimensionedInternalField();
+
+    return a*G - E;
 }
 
 

src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

@@ -30,7 +30,6 @@ License
 
 using namespace Foam::constant;
 
-
 const Foam::word
 Foam::radiation::radiativeIntensityRay::intensityPrefix("ILambda");
 
@@ -225,9 +224,15 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
               + fvm::Sp(k*omega_, ILambda_[lambdaI])
             ==
                 1.0/constant::mathematical::pi*omega_
-              * (
-                    k*blackBody_.bLambda(lambdaI)
-                  + absorptionEmission_.ECont(lambdaI)/4
+               *(
+                    // Remove aDisp from k
+                    (k - absorptionEmission_.aDisp(lambdaI))
+                   *blackBody_.bLambda(lambdaI)
+
+                  + absorptionEmission_.ECont(lambdaI)
+
+                    // Add EDisp term from parcels
+                  + absorptionEmission_.EDisp(lambdaI)
                 )
             );
         }
@@ -240,8 +245,14 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
            ==
                1.0/constant::mathematical::pi*omega_
              * (
-                   k*blackBody_.bLambda(lambdaI)
-                 + absorptionEmission_.ECont(lambdaI)/4
+                   // Remove aDisp from k
+                   (k - absorptionEmission_.aDisp(lambdaI))
+                  *blackBody_.bLambda(lambdaI)
+
+                 + absorptionEmission_.ECont(lambdaI)
+
+                   // Add EDisp term from parcels
+                 + absorptionEmission_.EDisp(lambdaI)
                )
             );
         }