diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index f6e57df771d897333bcc161ab381f66f4767f4b5..f0e9dab4036dd0cde2e72f377e37290158bd008f 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -6,7 +6,8 @@
     phic = min(interface.cAlpha()*phic, max(phic));
     surfaceScalarField phir(phic*interface.nHatf());
 
-    if (pimple.firstIter() && MULESCorr)
+    //***HGW if (pimple.firstIter() && MULESCorr)
+    if (MULESCorr)
     {
         fvScalarMatrix alpha1Eqn
         (
diff --git a/src/OpenFOAM/primitives/Tuple2/Tuple2.H b/src/OpenFOAM/primitives/Tuple2/Tuple2.H
index 26ca000575bba9c812ea540b3cf6b5fd76592ca4..4b7b45a1f84cb4a2eb673b91ea12411aa4cb69cc 100644
--- a/src/OpenFOAM/primitives/Tuple2/Tuple2.H
+++ b/src/OpenFOAM/primitives/Tuple2/Tuple2.H
@@ -25,7 +25,7 @@ Class
     Foam::Tuple2
 
 Description
-    A 2-tuple.
+    A 2-tuple for storing two objects of different types.
 
 SeeAlso
     Foam::Pair for storing two objects of identical types.
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
index 59ff974315322ef3ac6bcad5e39631ee047aa9e0..df4fe4a3ae05c45ba3cbf64e5b74a5cb38a66fdc 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -53,7 +53,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_("rho"),
     lookupGravity_(-1),
     g_(vector::zero),
-    relaxationFactor_(1)
+    curTimeIndex_(-1)
 {}
 
 
@@ -71,7 +71,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_(dict.lookupOrDefault<word>("rhoName", "rho")),
     lookupGravity_(-1),
     g_(vector::zero),
-    relaxationFactor_(dict.lookupOrDefault<scalar>("relaxationFactor", 1))
+    curTimeIndex_(-1)
 {
     if (rhoName_ == "rhoInf")
     {
@@ -115,7 +115,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_(ptf.rhoName_),
     lookupGravity_(ptf.lookupGravity_),
     g_(ptf.g_),
-    relaxationFactor_(ptf.relaxationFactor_)
+    curTimeIndex_(-1)
 {}
 
 
@@ -133,7 +133,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_(ptf.rhoName_),
     lookupGravity_(ptf.lookupGravity_),
     g_(ptf.g_),
-    relaxationFactor_(ptf.relaxationFactor_)
+    curTimeIndex_(-1)
 {}
 
 
@@ -203,6 +203,13 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     const Time& t = mesh.time();
     const pointPatch& ptPatch = this->patch();
 
+    // Store the motion state at the beginning of the time-step
+    if (curTimeIndex_ != this->db().time().timeIndex())
+    {
+        motion_.newTime();
+        curTimeIndex_ = this->db().time().timeIndex();
+    }
+
     // Patch force data is valid for the current positions, so
     // calculate the forces on the motion object from this data, then
     // update the positions
@@ -231,7 +238,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
         g_ = g.value();
     }
 
-    motion_.updateForce
+    motion_.updateAcceleration
     (
         f.forceEff() + g_*motion_.mass(),
         f.momentEff(),
@@ -240,11 +247,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
 
     Field<vector>::operator=
     (
-        relaxationFactor_
-       *(
-           motion_.currentPosition(initialPoints_)
-         - initialPoints_
-       )
+        motion_.currentPosition(initialPoints_) - initialPoints_
     );
 
     fixedValuePointPatchField<vector>::updateCoeffs();
@@ -267,9 +270,6 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::write(Ostream& os) const
         os.writeKeyword("g") << g_ << token::END_STATEMENT << nl;
     }
 
-    os.writeKeyword("relaxationFactor")
-        << relaxationFactor_ << token::END_STATEMENT << nl;
-
     motion_.write(os);
 
     initialPoints_.writeEntry("initialPoints", os);
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H
index db17ef3fae7444467023b5d9a7f055de18e451b4..355b44a99463bfbf05589b2fad2354770deabfc4 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H
@@ -84,8 +84,8 @@ class sixDoFRigidBodyDisplacementPointPatchVectorField
         //- Gravity vector to store when not available from the db
         vector g_;
 
-        //- Optional under-relaxation factor for the motion
-        scalar relaxationFactor_;
+        //- Current time index (used for updating)
+        label curTimeIndex_;
 
 
 public:
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
index cb6d4cbe63e20616c8017fd7e3f59e9095bd2576..4e218fc0d316944711bffbe75cf7cdb2663fdec7 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
@@ -164,6 +164,7 @@ void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
 Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
 :
     motionState_(),
+    motionState0_(),
     restraints_(),
     restraintNames_(),
     constraints_(),
@@ -173,8 +174,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
     initialQ_(I),
     momentOfInertia_(diagTensor::one*VSMALL),
     mass_(VSMALL),
-    cDamp_(0.0),
-    aLim_(VGREAT),
+    aRelax_(1.0),
     report_(false)
 {}
 
@@ -191,8 +191,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     const point& initialCentreOfMass,
     const tensor& initialQ,
     const diagTensor& momentOfInertia,
-    scalar cDamp,
-    scalar aLim,
+    scalar aRelax,
     bool report
 )
 :
@@ -205,6 +204,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
         pi,
         tau
     ),
+    motionState0_(motionState_),
     restraints_(),
     restraintNames_(),
     constraints_(),
@@ -214,8 +214,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     initialQ_(initialQ),
     momentOfInertia_(momentOfInertia),
     mass_(mass),
-    cDamp_(cDamp),
-    aLim_(aLim),
+    aRelax_(aRelax),
     report_(report)
 {}
 
@@ -223,6 +222,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
 Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
 :
     motionState_(dict),
+    motionState0_(motionState_),
     restraints_(),
     restraintNames_(),
     constraints_(),
@@ -238,8 +238,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
     ),
     momentOfInertia_(dict.lookup("momentOfInertia")),
     mass_(readScalar(dict.lookup("mass"))),
-    cDamp_(dict.lookupOrDefault<scalar>("accelerationDampingCoeff", 0.0)),
-    aLim_(dict.lookupOrDefault<scalar>("accelerationLimit", VGREAT)),
+    aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
     report_(dict.lookupOrDefault<Switch>("report", false))
 {
     addRestraints(dict);
@@ -254,6 +253,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
 )
 :
     motionState_(sDoFRBM.motionState_),
+    motionState0_(sDoFRBM.motionState0_),
     restraints_(sDoFRBM.restraints_),
     restraintNames_(sDoFRBM.restraintNames_),
     constraints_(sDoFRBM.constraints_),
@@ -263,8 +263,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     initialQ_(sDoFRBM.initialQ_),
     momentOfInertia_(sDoFRBM.momentOfInertia_),
     mass_(sDoFRBM.mass_),
-    cDamp_(sDoFRBM.cDamp_),
-    aLim_(sDoFRBM.aLim_),
+    aRelax_(sDoFRBM.aRelax_),
     report_(sDoFRBM.report_)
 {}
 
@@ -376,91 +375,79 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
 
     if (Pstream::master())
     {
-        vector aClip = a();
-        scalar aMag = mag(aClip);
-
-        if (aMag > SMALL)
-        {
-            aClip /= aMag;
-        }
-
-        if (aMag > aLim_)
-        {
-            WarningIn
-            (
-                "void Foam::sixDoFRigidBodyMotion::updatePosition"
-                "("
-                    "scalar, "
-                    "scalar"
-                ")"
-            )
-                << "Limited acceleration " << a()
-                << " to " << aClip*aLim_
-                << endl;
-        }
-
-        v() += 0.5*(1 - cDamp_)*deltaT0*aClip*min(aMag, aLim_);
-
-        pi() += 0.5*(1 - cDamp_)*deltaT0*tau();
+        v() = v0() + 0.5*deltaT0*a();
+        pi() = pi0() + 0.5*deltaT0*tau();
 
         // Leapfrog move part
-        centreOfMass() += deltaT*v();
+        centreOfMass() = centreOfMass0() + deltaT*v();
 
         // Leapfrog orientation adjustment
-        rotate(Q(), pi(), deltaT);
+        Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
+        Q() = Qpi.first();
+        pi() = Qpi.second();
     }
 
     Pstream::scatter(motionState_);
 }
 
 
-void Foam::sixDoFRigidBodyMotion::updateForce
+void Foam::sixDoFRigidBodyMotion::updateAcceleration
 (
     const vector& fGlobal,
     const vector& tauGlobal,
     scalar deltaT
 )
 {
+    static bool first = true;
+
     // Second leapfrog velocity adjust part, required after motion and
-    // force calculation
+    // acceleration calculation
 
     if (Pstream::master())
     {
-        a() = fGlobal/mass_;
+        // Save the previous iteration accelerations for relaxation
+        vector aPrevIter = a();
+        vector tauPrevIter = tau();
 
+        // Calculate new accelerations
+        a() = fGlobal/mass_;
         tau() = (Q().T() & tauGlobal);
-
         applyRestraints();
 
+        // Relax accelerations on all but first iteration
+        if (!first)
+        {
+            a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
+            tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
+        }
+        first = false;
+
+        // Apply constraints after relaxation
         applyConstraints(deltaT);
 
-        vector aClip = a();
-        scalar aMag = mag(aClip);
+        // Correct velocities
+        v() += 0.5*deltaT*a();
+        pi() += 0.5*deltaT*tau();
 
-        if (aMag > SMALL)
+        if (report_)
         {
-            aClip /= aMag;
+            status();
         }
+    }
+
+    Pstream::scatter(motionState_);
+}
 
-        if (aMag > aLim_)
-        {
-            WarningIn
-            (
-                "void Foam::sixDoFRigidBodyMotion::updateForce"
-                "("
-                    "const vector&, "
-                    "const vector&, "
-                    "scalar"
-                ")"
-            )
-                << "Limited acceleration " << a()
-                << " to " << aClip*aLim_
-                << endl;
-        }
 
-        v() += 0.5*(1 - cDamp_)*deltaT*aClip*min(aMag, aLim_);
+void Foam::sixDoFRigidBodyMotion::updateVelocity(scalar deltaT)
+{
+    // Second leapfrog velocity adjust part, required after motion and
+    // acceleration calculation
 
-        pi() += 0.5*(1 - cDamp_)*deltaT*tau();
+    if (Pstream::master())
+    {
+        v() += 0.5*deltaT*a();
+        pi() += 0.5*deltaT*tau();
 
         if (report_)
         {
@@ -472,7 +459,7 @@ void Foam::sixDoFRigidBodyMotion::updateForce
 }
 
 
-void Foam::sixDoFRigidBodyMotion::updateForce
+void Foam::sixDoFRigidBodyMotion::updateAcceleration
 (
     const pointField& positions,
     const vectorField& forces,
@@ -493,7 +480,7 @@ void Foam::sixDoFRigidBodyMotion::updateForce
         }
     }
 
-    updateForce(fGlobal, tauGlobal, deltaT);
+    updateAcceleration(fGlobal, tauGlobal, deltaT);
 }
 
 
@@ -506,19 +493,14 @@ Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition
 ) const
 {
     vector vTemp = v() + deltaT*(a() + deltaForce/mass_);
-
     vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
-
-    point centreOfMassTemp = centreOfMass() + deltaT*vTemp;
-
-    tensor QTemp = Q();
-
-    rotate(QTemp, piTemp, deltaT);
+    point centreOfMassTemp = centreOfMass0() + deltaT*vTemp;
+    Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
 
     return
     (
         centreOfMassTemp
-      + (QTemp & initialQ_.T() & (pInitial - initialCentreOfMass_))
+      + (QpiTemp.first() & initialQ_.T() & (pInitial - initialCentreOfMass_))
     );
 }
 
@@ -531,12 +513,9 @@ Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
 ) const
 {
     vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
+    Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
 
-    tensor QTemp = Q();
-
-    rotate(QTemp, piTemp, deltaT);
-
-    return (QTemp & initialQ_.T() & vInitial);
+    return (QpiTemp.first() & initialQ_.T() & vInitial);
 }
 
 
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
index 4ead010504bfe5a22ec2bfdee2a4b1b6eb1481a3..dcdc9f441bbd26f8d07ff9515cfb78aa10c5af2f 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.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
@@ -59,7 +59,7 @@ SourceFiles
 #include "pointField.H"
 #include "sixDoFRigidBodyMotionRestraint.H"
 #include "sixDoFRigidBodyMotionConstraint.H"
-
+#include "Tuple2.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -87,6 +87,9 @@ class sixDoFRigidBodyMotion
         //- Motion state data object
         sixDoFRigidBodyMotionState motionState_;
 
+        //- Motion state data object for previous time-step
+        sixDoFRigidBodyMotionState motionState0_;
+
         //- Restraints on the motion
         PtrList<sixDoFRigidBodyMotionRestraint> restraints_;
 
@@ -116,14 +119,8 @@ class sixDoFRigidBodyMotion
         //- Mass of the body
         scalar mass_;
 
-        //- Acceleration damping coefficient.  Modify applied acceleration:
-        //  v1 = v0 + a*dt - cDamp*a*dt
-        //     = v0 + dt*f*(1 - cDamp)/m
-        //  Increases effective mass by 1/(1 - cDamp).
-        scalar cDamp_;
-
-        //- Acceleration magnitude limit - clips large accelerations
-        scalar aLim_;
+        //- Acceleration relaxation coefficient
+        scalar aRelax_;
 
         //- Switch to turn reporting of motion data on and off
         Switch report_;
@@ -143,8 +140,14 @@ class sixDoFRigidBodyMotion
         //  frame z-axis by the given angle
         inline tensor rotationTensorZ(scalar deltaT) const;
 
-        //- Apply rotation tensors to Q for the given torque (pi) and deltaT
-        inline void rotate(tensor& Q, vector& pi, scalar deltaT) const;
+        //- Apply rotation tensors to Q0 for the given torque (pi) and deltaT
+        //  and return the rotated Q and pi as a tuple
+        inline Tuple2<tensor, vector> rotate
+        (
+            const tensor& Q0,
+            const vector& pi,
+            const scalar deltaT
+        ) const;
 
         //- Apply the restraints to the object
         void applyRestraints();
@@ -152,6 +155,7 @@ class sixDoFRigidBodyMotion
         //- Apply the constraints to the object
         void applyConstraints(scalar deltaT);
 
+
         // Access functions retained as private because of the risk of
         // confusion over what is a body local frame vector and what is global
 
@@ -199,6 +203,21 @@ class sixDoFRigidBodyMotion
             //- Return access to torque
             inline const vector& tau() const;
 
+            //- Return access to the orientation at previous time-step
+            inline const tensor& Q0() const;
+
+            //- Return access to velocity at previous time-step
+            inline const vector& v0() const;
+
+            //- Return access to acceleration at previous time-step
+            inline const vector& a0() const;
+
+            //- Return access to angular momentum at previous time-step
+            inline const vector& pi0() const;
+
+            //- Return access to torque at previous time-step
+            inline const vector& tau0() const;
+
 
         // Edit
 
@@ -244,8 +263,7 @@ public:
             const point& initialCentreOfMass,
             const tensor& initialQ,
             const diagTensor& momentOfInertia,
-            scalar cDamp = 0.0,
-            scalar aLim = VGREAT,
+            scalar aRelax = 1.0,
             bool report = false
         );
 
@@ -279,18 +297,22 @@ public:
             scalar deltaT0
         );
 
-        //- Second leapfrog velocity adjust part, required after motion and
-        // force calculation
-        void updateForce
+        //- Second leapfrog velocity adjust part
+        //  required after motion and force calculation
+        void updateAcceleration
         (
             const vector& fGlobal,
             const vector& tauGlobal,
             scalar deltaT
         );
 
+        //- Second leapfrog velocity adjust part
+        //  required after motion and force calculation
+        void updateVelocity(scalar deltaT);
+
         //- Global forces supplied at locations, calculating net force
         //  and moment
-        void updateForce
+        void updateAcceleration
         (
             const pointField& positions,
             const vectorField& forces,
@@ -354,6 +376,9 @@ public:
             //- Return const access to the centre of mass
             inline const point& centreOfMass() const;
 
+            //- Return const access to the centre of mass at previous time-step
+            inline const point& centreOfMass0() const;
+
             //- Return access to the inertia tensor
             inline const diagTensor& momentOfInertia() const;
 
@@ -366,6 +391,9 @@ public:
 
         // Edit
 
+            //- Store the motion state at the beginning of the time-step
+            inline void newTime();
+
             //- Return non-const access to the centre of mass
             inline point& centreOfMass();
 
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
index 19f6da1711868eb0c20573e63eb836ba718569b4..ec41af0d420ffbbfcab804f3c17acf9bfb271334 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
@@ -61,16 +61,19 @@ Foam::sixDoFRigidBodyMotion::rotationTensorZ(scalar phi) const
 }
 
 
-inline void Foam::sixDoFRigidBodyMotion::rotate
+inline Foam::Tuple2<Foam::tensor, Foam::vector>
+Foam::sixDoFRigidBodyMotion::rotate
 (
-    tensor& Q,
-    vector& pi,
-    scalar deltaT
+    const tensor& Q0,
+    const vector& pi0,
+    const scalar deltaT
 ) const
 {
-    tensor R;
+    Tuple2<tensor, vector> Qpi(Q0, pi0);
+    tensor& Q = Qpi.first();
+    vector& pi = Qpi.second();
 
-    R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
+    tensor R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
     pi = pi & R;
     Q = Q & R;
 
@@ -89,6 +92,8 @@ inline void Foam::sixDoFRigidBodyMotion::rotate
     R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
     pi = pi & R;
     Q = Q & R;
+
+    return Qpi;
 }
 
 
@@ -176,6 +181,36 @@ inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau() const
 }
 
 
+inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q0() const
+{
+    return motionState0_.Q();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::v0() const
+{
+    return motionState0_.v();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::a0() const
+{
+    return motionState0_.a();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::pi0() const
+{
+    return motionState0_.pi();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau0() const
+{
+    return motionState0_.tau();
+}
+
+
 inline Foam::point& Foam::sixDoFRigidBodyMotion::initialCentreOfMass()
 {
     return initialCentreOfMass_;
@@ -281,6 +316,12 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const
 }
 
 
+inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass0() const
+{
+    return motionState0_.centreOfMass();
+}
+
+
 inline const Foam::diagTensor&
 Foam::sixDoFRigidBodyMotion::momentOfInertia() const
 {
@@ -300,6 +341,12 @@ inline bool Foam::sixDoFRigidBodyMotion::report() const
 }
 
 
+inline void Foam::sixDoFRigidBodyMotion::newTime()
+{
+    motionState0_ = motionState_;
+}
+
+
 inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass()
 {
     return motionState_.centreOfMass();
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
index 5cc32dbe5431a83b927e28e9a52f78ed1f82db14..4f6b7b4007ebabdbaeef3162f080871b281a2a8d 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
@@ -40,10 +40,8 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
         << momentOfInertia_ << token::END_STATEMENT << nl;
     os.writeKeyword("mass")
         << mass_ << token::END_STATEMENT << nl;
-    os.writeKeyword("accelerationDampingCoeff")
-        << cDamp_ << token::END_STATEMENT << nl;
-    os.writeKeyword("accelerationLimit")
-        << aLim_ << token::END_STATEMENT << nl;
+    os.writeKeyword("accelerationRelaxation")
+        << aRelax_ << token::END_STATEMENT << nl;
     os.writeKeyword("report")
         << report_ << token::END_STATEMENT << nl;
 
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
index af8a089907a9f5561994c3f7ffe195425f696c65..b24d0a0678999d6dc360a9f9071efdbd509187b9 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -160,7 +160,11 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     }
 
     // Do not modify the accelerations, except with gravity, where available
-    motion_.updateForce(gravity*motion_.mass(), vector::zero, t.deltaTValue());
+    motion_.updateAcceleration
+    (
+        gravity*motion_.mass(),
+        vector::zero, t.deltaTValue()
+    );
 
     Field<vector>::operator=
     (
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
index efed1dd047d74766c45ce412b36d29ca7ddfdd16..2a0a90efa67a634198a39d9b79ec7f9606a5be1a 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
@@ -17,7 +17,7 @@ FoamFile
 
 application     interDyMFoam;
 
-startFrom       latestTime;
+startFrom       startTime;
 
 startTime       0;
 
@@ -29,7 +29,7 @@ deltaT          0.01;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.1;
+writeInterval   0.2;
 
 purgeWrite      0;
 
@@ -47,9 +47,9 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.5;
-maxAlphaCo      0.5;
-maxDeltaT       0.01;
+maxCo           5;
+maxAlphaCo      5;
+maxDeltaT       1;
 
 libs
 (
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
index 852f5e0c74425f84b7f2d62fe4dda0a2ab1fec4a..df7685da654a4a4fa7b159dc88722b7c7a521f5a 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
@@ -27,7 +27,7 @@ gradSchemes
 
 divSchemes
 {
-    div(rho*phi,U)  Gauss vanLeerV;
+    div(rhoPhi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
     div(phi,k)      Gauss upwind;
@@ -55,7 +55,7 @@ fluxRequired
     default         no;
     p_rgh;
     pcorr;
-    alpha;
+    alpha.water;
 }
 
 
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
index d3c114226bcbad9cfd73fbe56db2662662c3a4b9..c539caec767bd1663cc77089af78b605b21446b0 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    cellDisplacement
+    "cellDisplacement.*"
     {
         solver          GAMG;
         tolerance       1e-5;
@@ -29,14 +29,24 @@ solvers
         mergeLevels     1;
     }
 
-    alpha.water
+    "alpha.water.*"
     {
         nAlphaCorr      1;
-        nAlphaSubCycles 3;
+        nAlphaSubCycles 1;
         cAlpha          1;
+
+        alphaOuterCorrectors  yes;
+
+        MULESCorr       yes;
+        nLimiterIter    5;
+
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-8;
+        relTol          0;
     }
 
-    pcorr
+    "pcorr.*"
     {
         solver          PCG;
         preconditioner
@@ -120,9 +130,11 @@ solvers
 PIMPLE
 {
     momentumPredictor no;
-    nCorrectors     2;
+    nOuterCorrectors 6;
+    nCorrectors      1;
     nNonOrthogonalCorrectors 0;
     correctPhi      yes;
+    moveMeshOuterCorrectors yes;
 }
 
 relaxationFactors
@@ -132,7 +144,7 @@ relaxationFactors
     }
     equations
     {
-        "(U|k|epsilon|omega|R|nuTilda).*" 1;
+        ".*" 1;
     }
 }