ENH: Improvements to particle integration methods

src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C

@@ -61,7 +61,7 @@ Foam::Analytical<Type>::integrate
 (
     const Type& phi,
     const scalar dt,
-    const Type& alpha,
+    const Type& alphaBeta,
     const scalar beta
 ) const
 {
@@ -69,9 +69,18 @@ Foam::Analytical<Type>::integrate
 
     const scalar expTerm = exp(min(50, -beta*dt));
 
-    retValue.average() =
-        alpha + (phi - alpha)*(1 - expTerm)/(beta*dt + ROOTVSMALL);
-    retValue.value() =  alpha + (phi - alpha)*expTerm;
+    if (beta > ROOTVSMALL)
+    {
+        const Type alpha = alphaBeta/beta;
+        retValue.average() = alpha + (phi - alpha)*(1 - expTerm)/(beta*dt);
+        retValue.value() =  alpha + (phi - alpha)*expTerm;
+    }
+    else
+    {
+        retValue.average() = phi;
+        retValue.value() =  phi;
+    }
+
 
     return retValue;
 }

src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H

@@ -83,7 +83,7 @@ public:
         (
             const Type& phi,
             const scalar dt,
-            const Type& alpha,
+            const Type& alphaBeta,
             const scalar beta
         ) const;
 };

src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C

@@ -61,12 +61,12 @@ Foam::Euler<Type>::integrate
 (
     const Type& phi,
     const scalar dt,
-    const Type& alpha,
+    const Type& alphaBeta,
     const scalar beta
 ) const
 {
     typename IntegrationScheme<Type>::integrationResult retValue;
-    retValue.value() = (phi + beta*dt*alpha)/(1.0 + beta*dt);
+    retValue.value() = (phi + alphaBeta*dt)/(1.0 + beta*dt);
     retValue.average() = 0.5*(phi + retValue.value());
 
     return retValue;

src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H

@@ -80,7 +80,7 @@ public:
         (
             const Type& phi,
             const scalar dt,
-            const Type& alpha,
+            const Type& alphaBeta,
             const scalar beta
         ) const;
 };

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C

@@ -62,7 +62,7 @@ Foam::IntegrationScheme<Type>::integrate
 (
     const Type& phi,
     const scalar dt,
-    const Type& alpha,
+    const Type& alphaBeta,
     const scalar beta
 ) const
 {

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H

@@ -185,7 +185,7 @@ public:
         (
             const Type& phi,
             const scalar dt,
-            const Type& alpha,
+            const Type& alphaBeta,
             const scalar beta
         ) const;
 };

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

@@ -192,21 +192,20 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     const parcelType& p = static_cast<const parcelType&>(*this);
     const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
     const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
-    forceSuSp Feff = Fcp + Fncp;
-    Feff.Sp() += ROOTVSMALL;
+    const forceSuSp Feff = Fcp + Fncp;
 
 
     // New particle velocity
     //~~~~~~~~~~~~~~~~~~~~~~
 
     // Update velocity - treat as 3-D
-    const vector ap = Uc_ + (Feff.Su() + Su)/Feff.Sp();
+    const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/mass;
     const scalar bp = Feff.Sp()/mass;
 
     Spu = Feff.Sp();
 
     IntegrationScheme<vector>::integrationResult Ures =
-        td.cloud().UIntegrator().integrate(U, dt, ap, bp);
+        td.cloud().UIntegrator().integrate(U, dt, abp, bp);
 
     vector Unew = Ures.value();
 

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

@@ -333,7 +333,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     // Integrate to find the new parcel temperature
     IntegrationScheme<scalar>::integrationResult Tres =
-        td.cloud().TIntegrator().integrate(T, dt, ap, bp);
+        td.cloud().TIntegrator().integrate(T, dt, ap*bp, bp);
 
     scalar Tnew = max(Tres.value(), td.cloud().constProps().TMin());