BUG: Corrected Pilch-Erdman break-up model

src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C

@@ -36,7 +36,7 @@ Foam::PilchErdman<CloudType>::PilchErdman
 :
     BreakupModel<CloudType>(dict, owner, typeName),
     B1_(0.375),
-    B2_(0.236)
+    B2_(0.2274)
 {
     if (!this->defaultCoeffs(true))
     {
@@ -90,59 +90,80 @@ bool Foam::PilchErdman<CloudType>::update
     scalar& massChild
 )
 {
-    scalar semiMass = nParticle*pow3(d);
+    // Weber number - eq (1)
-    scalar We = 0.5*rhoc*sqr(Urmag)*d/sigma;
+    scalar We = rhoc*sqr(Urmag)*d/sigma;
+
+    // Ohnesorge number - eq (2)
     scalar Oh = mu/sqrt(rho*d*sigma);
 
-    scalar Wec = 6.0*(1.0 + 1.077*pow(Oh, 1.6));
+    // Critical Weber number - eq (5)
+    scalar Wec = 12.0*(1.0 + 1.077*pow(Oh, 1.6));
 
     if (We > Wec)
     {
-        // We > 1335, wave crest stripping
+        // We > 2670, wave crest stripping - eq (12)
         scalar taubBar = 5.5;
 
-        if (We < 1335)
+        if (We < 2670)
         {
-            if (We > 175.0)
+            if (We > 351)
             {
-                // sheet stripping
+                // sheet stripping - eq (11)
-                taubBar = 0.766*pow(2.0*We - 12.0, 0.25);
+                taubBar = 0.766*pow(We - 12.0, 0.25);
             }
-            else if (We > 22.0)
+            else if (We > 45)
             {
-                // Bag-and-stamen breakup
+                // bag-and-stamen breakup  - eq (10)
-                taubBar = 14.1*pow(2.0*We - 12.0, -0.25);
+                taubBar = 14.1*pow(We - 12.0, 0.25);
             }
-            else if (We > 9.0)
+            else if (We > 18)
             {
-                // Bag breakup
+                // bag breakup - eq (9)
-                taubBar = 2.45*pow(2.0*We - 12.0, 0.25);
+                taubBar = 2.45*pow(We - 12.0, 0.25);
             }
-            else if (We > 6.0)
+            else if (We > 12)
             {
-                // Vibrational breakup
+                // vibrational breakup - eq (8)
-                taubBar = 6.0*pow(2.0*We - 12.0, -0.25);
+                taubBar = 6.0*pow(We - 12.0, -0.25);
+            }
+            else
+            {
+                // no break-up
+                taubBar = GREAT;
             }
         }
 
         scalar rho12 = sqrt(rhoc/rho);
 
-        scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*taubBar*taubBar);
+        // velocity of fragmenting drop - eq (20)
+        scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*sqr(taubBar));
+
+        // maximum stable diameter - eq (33)
         scalar Vd1 = sqr(1.0 - Vd/Urmag);
         Vd1 = max(Vd1, SMALL);
-        scalar Ds = 2.0*Wec*sigma/(Vd1*rhoc*sqr(Urmag));
+        scalar dStable = Wec*sigma/(Vd1*rhoc*sqr(Urmag));
-        scalar A = Urmag*rho12/d;
 
-        scalar taub = taubBar/A;
+        if (d < dStable)
+        {
+            // droplet diameter already stable = no break-up
+            // - do not update d and nParticle
+            return false;
+        }
+        else
+        {
+            scalar semiMass = nParticle*pow3(d);
 
+            // invert eq (3) to create a dimensional break-up time
+            scalar taub = taubBar*d/(Urmag*rho12);
+
+            // update droplet diameter according to the rate eq (implicitly)
             scalar frac = dt/taub;
-
+            d = (d + frac*dStable)/(1.0 + frac);
-        // update the droplet diameter according to the rate eq. (implicitly)
-        d = (d + frac*Ds)/(1.0 + frac);
 
             // correct the number of particles to conserve mass
             nParticle = semiMass/pow3(d);
         }
+    }
 
     return false;
 }

src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Class
     Foam::PilchErdman
 
 Description
-    secondary breakup model
+    Particle secondary breakup model, based on the reference:
 
     @verbatim
     Pilch, M. and Erdman, C.A.
@@ -35,6 +35,24 @@ Description
     Int. J. Multiphase Flows 13 (1987), 741-757
     @endverbatim
 
+    The droplet fragment velocity is described by the equation:
+
+    \f[
+        V_d = V sqrt(epsilon)(B1 T + B2 T^2)
+    \f]
+
+    Where:
+        V_d     : fragment velocity
+        V       : magnitude of the relative velocity
+        epsilon : density ratio (rho_carrier/rho_droplet)
+        T       : characteristic break-up time
+        B1, B2  : model input coefficients
+
+    The authors suggest that:
+        compressible flow   : B1 = 0.75*1.0; B2 = 3*0.116
+        incompressible flow : B1 = 0.75*0.5; B2 = 3*0.0758
+
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef PilchErdman_H