From 2d017d6cc46f4cd8a94c3939dede930ff4d51d05 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Tue, 8 Sep 2009 19:29:37 +0100
Subject: [PATCH] Hacking in a sympletic, leapfrog integrator for interparticle
 forces. Conserves kinetic energy correctly.

---
 .../InteractingKinematicCloud.C                  | 16 +++++++++++++++-
 .../InteractingKinematicParcel.C                 |  3 +--
 2 files changed, 16 insertions(+), 3 deletions(-)

diff --git a/src/lagrangian/intermediate/clouds/Templates/InteractingKinematicCloud/InteractingKinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/InteractingKinematicCloud/InteractingKinematicCloud.C
index d063b2c1ff..c05b0fe4d3 100644
--- a/src/lagrangian/intermediate/clouds/Templates/InteractingKinematicCloud/InteractingKinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/InteractingKinematicCloud/InteractingKinematicCloud.C
@@ -254,10 +254,24 @@ void Foam::InteractingKinematicCloud<ParcelType>::evolve()
         resetSourceTerms();
     }
 
-    this->collision().collide();
+    const scalar deltaT = mesh_.time().deltaT().value();
+
+    forAllIter(typename InteractingKinematicCloud<ParcelType>, *this, iter)
+    {
+        ParcelType& p = iter();
+        p.U() += 0.5*deltaT*p.f()/p.mass();
+    }
 
     Cloud<ParcelType>::move(td);
 
+    this->collision().collide();
+
+    forAllIter(typename InteractingKinematicCloud<ParcelType>, *this, iter)
+    {
+        ParcelType& p = iter();
+        p.U() += 0.5*deltaT*p.f()/p.mass();
+    }
+
     postEvolve();
 }
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.C
index 702abded1e..48214ecd01 100644
--- a/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/InteractingKinematicParcel/InteractingKinematicParcel.C
@@ -167,14 +167,13 @@ const Foam::vector Foam::InteractingKinematicParcel<ParcelType>::calcVelocity
         mass*td.cloud().forces().calcCoupled(cellI, dt, rhoc_, rho, Uc_, U);
     const vector FNonCoupled =
         mass*td.cloud().forces().calcNonCoupled(cellI, dt, rhoc_, rho, Uc_, U);
-    const vector& FCollision = f();
 
     // New particle velocity
     //~~~~~~~~~~~~~~~~~~~~~~
 
     // Update velocity - treat as 3-D
     const scalar As = this->areaS(d);
-    const vector ap = Uc_ + (FCoupled + FNonCoupled + FCollision + Su)/(utc*As);
+    const vector ap = Uc_ + (FCoupled + FNonCoupled + Su)/(utc*As);
     const scalar bp = 6.0*utc/(rho*d);
 
     IntegrationScheme<vector>::integrationResult Ures =