sixDoFRigidBodyMotion: Time integration now switches between symplectic and Crank-Nicolson

For explicit motion (and the first iteration of iterative motion
correction) the 2nd-order symplectic motion integrator is used.

For iterative correction a form of lagged Crank-Nicolson is used in
which the current time-step values correspond to the current iteration.
This converges to a 2nd-order implicit solution.

src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -204,18 +204,14 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     const pointPatch& ptPatch = this->patch();
 
     // Store the motion state at the beginning of the time-step
+    bool firstIter = false;
     if (curTimeIndex_ != t.timeIndex())
     {
         motion_.newTime();
         curTimeIndex_ = t.timeIndex();
+        firstIter = true;
     }
 
-    // Patch force data is valid for the current positions, so
-    // calculate the forces on the motion object from this data, then
-    // update the positions
-
-    motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
-
     dictionary forcesDict;
 
     forcesDict.add("type", forces::typeName);
@@ -228,6 +224,11 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
 
     f.calcForcesMoment();
 
+    // Patch force data is valid for the current positions, so
+    // calculate the forces on the motion object from this data, then
+    // update the positions
+    motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
+
     // Get the forces on the patch faces at the current positions
 
     if (lookupGravity_ == 1)
@@ -243,9 +244,11 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
 
     motion_.updateAcceleration
     (
+        firstIter,
         ramp*(f.forceEff() + motion_.mass()*g_),
         ramp*(f.momentEff() + motion_.mass()*(motion_.momentArm() ^ g_)),
-        t.deltaTValue()
+        t.deltaTValue(),
+        t.deltaT0Value()
     );
 
     Field<vector>::operator=

src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -152,13 +152,15 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     const Time& t = mesh.time();
 
     // Store the motion state at the beginning of the time-step
+    bool firstIter = false;
     if (curTimeIndex_ != t.timeIndex())
     {
         motion_.newTime();
         curTimeIndex_ = t.timeIndex();
+        firstIter = true;
     }
 
-    motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
+    motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
 
     vector gravity = vector::zero;
 
@@ -173,9 +175,11 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     // Do not modify the accelerations, except with gravity, where available
     motion_.updateAcceleration
     (
+        firstIter,
         gravity*motion_.mass(),
         vector::zero,
-        t.deltaTValue()
+        t.deltaTValue(),
+        t.deltaT0Value()
     );
 
     Field<vector>::operator=

src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C

@@ -259,22 +259,38 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
 
 void Foam::sixDoFRigidBodyMotion::updatePosition
 (
+    bool firstIter,
     scalar deltaT,
     scalar deltaT0
 )
 {
-    // First leapfrog velocity adjust and motion part, required before
-    // force calculation
-
     if (Pstream::master())
     {
-        v() = tConstraints_ & (v0() + aDamp_*0.5*deltaT0*a());
+        if (firstIter)
-        pi() = rConstraints_ & (pi0() + aDamp_*0.5*deltaT0*tau());
+        {
+            // First simplectic step:
+            //     Half-step for linear and angular velocities
+            //     Update position and orientation
 
-        // Leapfrog move part
+            v() = tConstraints_ & (v0() + aDamp_*0.5*deltaT0*a());
-        centreOfRotation() = centreOfRotation0() + deltaT*v();
+            pi() = rConstraints_ & (pi0() + aDamp_*0.5*deltaT0*tau());
 
-        // Leapfrog orientation adjustment
+            centreOfRotation() = centreOfRotation0() + deltaT*v();
+        }
+        else
+        {
+            // For subsequent iterations use Crank-Nicolson
+
+            v() = tConstraints_
+              & (v0() + aDamp_*0.5*deltaT*(a() + motionState0_.a()));
+            pi() = rConstraints_
+              & (pi0() + aDamp_*0.5*deltaT*(tau() + motionState0_.tau()));
+
+            centreOfRotation() =
+                centreOfRotation0() + 0.5*deltaT*(v() + motionState0_.v());
+        }
+
+        // Correct orientation
         Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
         Q() = Qpi.first();
         pi() = rConstraints_ & Qpi.second();
@@ -286,16 +302,15 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
 
 void Foam::sixDoFRigidBodyMotion::updateAcceleration
 (
+    bool firstIter,
     const vector& fGlobal,
     const vector& tauGlobal,
-    scalar deltaT
+    scalar deltaT,
+    scalar deltaT0
 )
 {
     static bool first = false;
 
-    // Second leapfrog velocity adjust part, required after motion and
-    // acceleration calculation
-
     if (Pstream::master())
     {
         // Save the previous iteration accelerations for relaxation
@@ -315,9 +330,23 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
         }
         first = false;
 
-        // Correct velocities
+        if (firstIter)
-        v() += tConstraints_ & aDamp_*0.5*deltaT*a();
+        {
-        pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
+            // Second simplectic step:
+            //     Complete update of linear and angular velocities
+
+            v() += tConstraints_ & aDamp_*0.5*deltaT*a();
+            pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
+        }
+        else
+        {
+            // For subsequent iterations use Crank-Nicolson
+
+            v() = tConstraints_
+              & (v0() + aDamp_*0.5*deltaT*(a() + motionState0_.a()));
+            pi() = rConstraints_
+              & (pi0() + aDamp_*0.5*deltaT*(tau() + motionState0_.tau()));
+        }
 
         if (report_)
         {
@@ -329,6 +358,7 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
 }
 
 
+
 void Foam::sixDoFRigidBodyMotion::status() const
 {
     Info<< "6-DoF rigid body motion" << nl

src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,6 +42,9 @@ Description
     url = {http://link.aip.org/link/?JCP/107/5840/1},
     doi = {10.1063/1.474310}
 
+    Changes to Crank-Nicolson integration if subsequent iterations to handle
+    implicit coupling between pressure and motion.
+
     Can add restraints (e.g. a spring)
     and constraints (e.g. motion may only be on a plane).
 
@@ -298,18 +301,22 @@ public:
 
         // Update state
 
-            //- First leapfrog velocity adjust and motion part, required
+            //- First symplectic velocity adjust and motion part, required
             //  before force calculation.  Takes old timestep for variable
             //  timestep cases.
-            void updatePosition(scalar deltaT, scalar deltaT0);
+            //  Changes to Crank-Nicolson integration for subsequent iterations.
+            void updatePosition(bool firstIter, scalar deltaT, scalar deltaT0);
 
-            //- Second leapfrog velocity adjust part
+            //- Second symplectic velocity adjust part
-            //  required after motion and force calculation
+            //  required after motion and force calculation.
+            //  Changes to Crank-Nicolson integration for subsequent iterations.
             void updateAcceleration
             (
+                bool firstIter,
                 const vector& fGlobal,
                 const vector& tauGlobal,
-                scalar deltaT
+                scalar deltaT,
+                scalar deltaT0
             );
 
             //- Report the status of the motion

src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,6 +85,7 @@ Foam::sixDoFRigidBodyMotionSolver::sixDoFRigidBodyMotionSolver
     patchSet_(mesh.boundaryMesh().patchSet(patches_)),
     di_(readScalar(coeffDict().lookup("innerDistance"))),
     do_(readScalar(coeffDict().lookup("outerDistance"))),
+    test_(coeffDict().lookupOrDefault<Switch>("test", false)),
     rhoInf_(1.0),
     rhoName_(coeffDict().lookupOrDefault<word>("rhoName", "rho")),
     scale_
@@ -179,31 +180,15 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
             << " points." << exit(FatalError);
     }
 
-    // Store the motion state at the beginning of the time-step
+    // Store the motion state at the beginning of the time-stepbool
+    bool firstIter = false;
     if (curTimeIndex_ != this->db().time().timeIndex())
     {
         motion_.newTime();
         curTimeIndex_ = this->db().time().timeIndex();
+        firstIter = true;
     }
 
-    // Patch force data is valid for the current positions, so
-    // calculate the forces on the motion object from this data, then
-    // update the positions
-
-    motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
-
-    dictionary forcesDict;
-
-    forcesDict.add("type", forces::typeName);
-    forcesDict.add("patches", patches_);
-    forcesDict.add("rhoInf", rhoInf_);
-    forcesDict.add("rhoName", rhoName_);
-    forcesDict.add("CofR", motion_.centreOfRotation());
-
-    forces f("forces", db(), forcesDict);
-
-    f.calcForcesMoment();
-
     dimensionedVector g("g", dimAcceleration, vector::zero);
 
     if (db().foundObject<uniformDimensionedVectorField>("g"))
@@ -218,12 +203,51 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
     // scalar ramp = min(max((this->db().time().value() - 5)/10, 0), 1);
     scalar ramp = 1.0;
 
-    motion_.updateAcceleration
+    if (test_)
-    (
+    {
-        ramp*(f.forceEff() + motion_.mass()*g.value()),
+        motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
-        ramp*(f.momentEff() + motion_.mass()*(motion_.momentArm() ^ g.value())),
+
-        t.deltaTValue()
+        motion_.updateAcceleration
-    );
+        (
+            firstIter,
+            ramp*(motion_.mass()*g.value()),
+            ramp*(motion_.mass()*(motion_.momentArm() ^ g.value())),
+            t.deltaTValue(),
+            t.deltaT0Value()
+        );
+    }
+    else
+    {
+        dictionary forcesDict;
+
+        forcesDict.add("type", forces::typeName);
+        forcesDict.add("patches", patches_);
+        forcesDict.add("rhoInf", rhoInf_);
+        forcesDict.add("rhoName", rhoName_);
+        forcesDict.add("CofR", motion_.centreOfRotation());
+
+        forces f("forces", db(), forcesDict);
+
+        f.calcForcesMoment();
+
+        // Patch force data is valid for the current positions, so
+        // calculate the forces on the motion object from this data, then
+        // update the positions
+        motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
+
+        motion_.updateAcceleration
+        (
+            firstIter,
+            ramp*(f.forceEff() + motion_.mass()*g.value()),
+            ramp
+           *(
+               f.momentEff()
+             + motion_.mass()*(motion_.momentArm() ^ g.value())
+            ),
+            t.deltaTValue(),
+            t.deltaT0Value()
+        );
+    }
 
     // Update the displacements
     pointDisplacement_.internalField() =

src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.H

@@ -70,6 +70,10 @@ class sixDoFRigidBodyMotionSolver
         //- Outer morphing distance (limit of linear interpolation region)
         const scalar do_;
 
+        //- Switch for test-mode in which only the
+        //  gravitational body-force is applied
+        Switch test_;
+
         //- Reference density required by the forces object for
         //  incompressible calculations, required if rhoName == rhoInf
         scalar rhoInf_;