Reduce the stripped parcel mass due to evaporation

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1295

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C

@@ -68,42 +68,47 @@ void Foam::SprayParcel<ParcelType>::calc
     const CompositionModel<reactingCloudType>& composition =
         td.cloud().composition();
 
-    // check if parcel belongs to liquid core
+    // Check if parcel belongs to liquid core
     if (liquidCore() > 0.5)
     {
-        // liquid core parcels should not experience coupled forces
+        // Liquid core parcels should not experience coupled forces
         td.cloud().forces().setCalcCoupled(false);
     }
 
-    // get old mixture composition
+    // Get old mixture composition
     const scalarField& Y0(this->Y());
     scalarField X0(composition.liquids().X(Y0));
 
-    // check if we have critical or boiling conditions
+    // Check if we have critical or boiling conditions
     scalar TMax = composition.liquids().Tc(X0);
     const scalar T0 = this->T();
     const scalar pc0 = this->pc_;
     if (composition.liquids().pv(pc0, T0, X0) >= pc0*0.999)
     {
-        // set TMax to boiling temperature
+        // Set TMax to boiling temperature
         TMax = composition.liquids().pvInvert(pc0, X0);
     }
 
-    // set the maximum temperature limit
+    // Set the maximum temperature limit
     td.cloud().constProps().setTMax(TMax);
 
-    // store the parcel properties
+    // Store the parcel properties
     this->Cp() = composition.liquids().Cp(pc0, T0, X0);
     sigma_ = composition.liquids().sigma(pc0, T0, X0);
-    scalar rho0 = composition.liquids().rho(pc0, T0, X0);
+    const scalar rho0 = composition.liquids().rho(pc0, T0, X0);
     this->rho() = rho0;
+    const scalar mass0 = this->mass();
     mu_ = composition.liquids().mu(pc0, T0, X0);
 
     ParcelType::calc(td, dt, cellI);
 
     if (td.keepParticle)
     {
-        // update Cp, sigma, density and diameter due to change in temperature
+        // Reduce the stripped parcel mass due to evaporation
+        // assuming the number of particles remains unchanged
+        this->ms() -= this->ms()*(mass0 - this->mass())/mass0;
+
+        // Update Cp, sigma, density and diameter due to change in temperature
         // and/or composition
         scalar T1 = this->T();
         const scalarField& Y1(this->Y());
@@ -125,7 +130,7 @@ void Foam::SprayParcel<ParcelType>::calc
         {
             calcAtomization(td, dt, cellI);
 
-            // preserve the total mass/volume by increasing the number of
+            // Preserve the total mass/volume by increasing the number of
             // particles in parcels due to breakup
             scalar d2 = this->d();
             this->nParticle() *= pow3(d1/d2);
@@ -136,7 +141,7 @@ void Foam::SprayParcel<ParcelType>::calc
         }
     }
 
-    // restore coupled forces
+    // Restore coupled forces
     td.cloud().forces().setCalcCoupled(true);
 }
 
@@ -163,7 +168,7 @@ void Foam::SprayParcel<ParcelType>::calcAtomization
     scalar R = specie::RR/Wc;
     scalar Tav = atomization.Taverage(this->T(), this->Tc());
 
-    // calculate average gas density based on average temperature
+    // Calculate average gas density based on average temperature
     scalar rhoAv = this->pc()/(R*Tav);
 
     scalar soi = td.cloud().injectors().timeStart();
@@ -171,14 +176,14 @@ void Foam::SprayParcel<ParcelType>::calcAtomization
     const vector& pos = this->position();
     const vector& injectionPos = this->position0();
 
-    // disregard the continous phase when calculating the relative velocity
+    // Disregard the continous phase when calculating the relative velocity
     // (in line with the deactivated coupled assumption)
     scalar Urel = mag(this->U());
 
     scalar t0 = max(0.0, currentTime - this->age() - soi);
     scalar t1 = min(t0 + dt, td.cloud().injectors().timeEnd() - soi);
 
-    // this should be the vol flow rate from when the parcel was injected
+    // This should be the vol flow rate from when the parcel was injected
     scalar volFlowRate = td.cloud().injectors().volumeToInject(t0, t1)/dt;
 
     scalar chi = 0.0;
@@ -234,7 +239,7 @@ void Foam::SprayParcel<ParcelType>::calcBreakup
     scalar R = specie::RR/Wc;
     scalar Tav = td.cloud().atomization().Taverage(this->T(), this->Tc());
 
-    // calculate average gas density based on average temperature
+    // Calculate average gas density based on average temperature
     scalar rhoAv = this->pc()/(R*Tav);
     scalar muAv = this->muc();
     vector Urel = this->U() - this->Uc();
@@ -316,7 +321,7 @@ Foam::scalar Foam::SprayParcel<ParcelType>::chi
     const scalarField& X
 ) const
 {
-    // modifications to take account of the flash boiling on primary break-up
+    // Modifications to take account of the flash boiling on primary break-up
 
     typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
     const CompositionModel<reactingCloudType>& composition =
@@ -333,7 +338,7 @@ Foam::scalar Foam::SprayParcel<ParcelType>::chi
     {
         if (pv >= 0.999*pAmb)
         {
-            // liquid is boiling - calc boiling temperature
+            // Liquid is boiling - calc boiling temperature
 
             const liquidProperties& liq = composition.liquids().properties()[i];
             scalar TBoil = liq.pvInvert(p0);
@@ -368,10 +373,10 @@ void Foam::SprayParcel<ParcelType>::solveTABEq
     scalar r2 = r*r;
     scalar r3 = r*r2;
 
-    // inverse of characteristic viscous damping time
+    // Inverse of characteristic viscous damping time
     scalar rtd = 0.5*TABCmu*mu_/(this->rho()*r2);
 
-    // oscillation frequency (squared)
+    // Oscillation frequency (squared)
     scalar omega2 = TABComega*sigma_/(this->rho()*r3) - rtd*rtd;
 
     if (omega2 > 0)
@@ -383,7 +388,7 @@ void Foam::SprayParcel<ParcelType>::solveTABEq
         scalar y1 = this->y() - We;
         scalar y2 = this->yDot()/omega;
 
-        // update distortion parameters
+        // Update distortion parameters
         scalar c = cos(omega*dt);
         scalar s = sin(omega*dt);
         scalar e = exp(-rtd*dt);
@@ -402,7 +407,7 @@ void Foam::SprayParcel<ParcelType>::solveTABEq
     }
     else
     {
-        // reset distortion parameters
+        // Reset distortion parameters
         this->y() = 0;
         this->yDot() = 0;
     }