ENH: Cloud tracking - handle case where tracking stalls due to stepFraction

src/lagrangian/basic/particle/particle.C

@@ -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
@@ -34,6 +34,8 @@ const Foam::scalar Foam::particle::trackingCorrectionTol = 1e-5;
 
 const Foam::scalar Foam::particle::lambdaDistanceToleranceCoeff = 1e3*SMALL;
 
+const Foam::scalar Foam::particle::minStepFractionTol = 1e5*SMALL;
+
 namespace Foam
 {
     defineTypeNameAndDebug(particle, 0);

src/lagrangian/basic/particle/particle.H

@@ -307,6 +307,9 @@ public:
         //  for the denominator and numerator of lambda
         static const scalar lambdaDistanceToleranceCoeff;
 
+        //- Minimum stepFraction tolerance
+        static const scalar minStepFractionTol;
+
 
     // Constructors
 

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

@@ -271,6 +271,8 @@ bool Foam::KinematicParcel<ParcelType>::move
     scalar tEnd = (1.0 - p.stepFraction())*trackTime;
     const scalar dtMax = tEnd;
 
+    bool moving = true;
+
     while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
     {
         // Apply correction to position for reduced-D cases
@@ -281,22 +283,36 @@ bool Foam::KinematicParcel<ParcelType>::move
         // Set the Lagrangian time-step
         scalar dt = min(dtMax, tEnd);
 
-        // Remember which cell the parcel is in since this will change if
-        // a face is hit
+        // Cache the parcel current cell as this will change if a face is hit
         const label cellI = p.cell();
 
         const scalar magU = mag(U_);
-        if (p.active() && magU > ROOTVSMALL)
+        if (p.active())
         {
             const scalar d = dt*magU;
             const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
-            dt *=
-                dCorr/d
-               *p.trackToFace(p.position() + dCorr*U_/magU, td);
+            dt *= dCorr/d;
+
+            if (moving && (magU > ROOTVSMALL))
+            {
+                dt *= p.trackToFace(p.position() + dCorr*U_/magU, td);
+            }
         }
 
         tEnd -= dt;
-        p.stepFraction() = 1.0 - tEnd/trackTime;
+
+        scalar newStepFraction = 1.0 - tEnd/trackTime;
+
+        if
+        (
+            mag(p.stepFraction() - newStepFraction)
+          < particle::minStepFractionTol
+        )
+        {
+            moving = false;
+        }
+
+        p.stepFraction() = newStepFraction;
 
         // Avoid problems with extremely small timesteps
         if (dt > ROOTVSMALL)