Updates to sixDOFODE classes

2nd order accurate updates of force and moment during the ODE solution.

src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C

@@ -83,7 +83,7 @@ Foam::dimensionedVector Foam::geometricSixDOF::A
     // solved in the global coordinate system
     const dimensionedVector explicitForcing
     (
-        force() // External force
+        force(t) // External force
       - (linSpringCoeffs() & xR) // Spring force
       - (linDampingCoeffs() & uR) // Damping force
     );
@@ -155,7 +155,7 @@ Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
     const dimensionedVector explicitForcing
     (
         E(omega) // Euler part
-      + (RT & moment()) // External torque
+      + (RT & moment(t)) // External torque
     );
     const vector& efVal = explicitForcing.value();
     const diagTensor& I = momentOfInertia().value();

src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C

@@ -68,7 +68,7 @@ Foam::dimensionedVector Foam::quaternionSixDOF::A
     // solved in the global coordinate system
     const dimensionedVector explicitForcing
     (
-        force() // External force
+        force(t) // External force
       - (linSpringCoeffs() & xR) // Spring force
       - (linDampingCoeffs() & uR) // Damping force
     );
@@ -141,7 +141,7 @@ Foam::dimensionedVector Foam::quaternionSixDOF::OmegaDot
     const dimensionedVector explicitForcing
     (
         E(omega) // Euler part
-      + (R & moment()) // External torque
+      + (R & moment(t)) // External torque
     );
     const vector& efVal = explicitForcing.value();
     const diagTensor& I = momentOfInertia().value();

src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C

@@ -130,10 +130,19 @@ void Foam::sixDOFBodies::solve()
         Info << "Solving 6-DOF for " << names_[bodyI] << " in time"
          << tab << "T = " << runTime_.value() << " s" << endl;
 
+        // Note: set external force and moment needed to initialize the state
+        // of the sixDOFODE to correctly take into account multiple calls per
+        // time step
+        odes_[bodyI].setExternalForceAndMoment
+        (
+            dimensionedVector("zero", dimForce, vector::zero),
+            dimensionedVector("zero", dimForce*dimLength, vector::zero)
+        );
+
         solvers_[bodyI].solve
         (
+            runTime_.value() - runTime_.deltaT().value(),
             runTime_.value(),
-            runTime_.value() + runTime_.deltaT().value(),
             tol,
             runTime_.deltaT().value()
         );

src/ODE/sixDOF/sixDOFODE/sixDOFODE.C

@@ -208,6 +208,26 @@ void Foam::sixDOFODE::setState(const sixDOFODE& sd)
 }
 
 
+Foam::dimensionedVector Foam::sixDOFODE::force(const scalar t) const
+{
+    // Get ODE step fraction
+    const scalar alpha = odeStepFraction(t);
+
+    // Return linearly interpolated external force
+    return (alpha*oldStatePtr_->force() + (1 - alpha)*force());
+}
+
+
+Foam::dimensionedVector Foam::sixDOFODE::moment(const scalar t) const
+{
+    // Get ODE step fraction
+    const scalar alpha = odeStepFraction(t);
+
+    // Return linearly interpolated external moment
+    return (alpha*oldStatePtr_->moment() + (1 - alpha)*moment());
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::sixDOFODE::sixDOFODE(const IOobject& io)
@@ -294,7 +314,7 @@ Foam::sixDOFODE::sixDOFODE(const word& name, const sixDOFODE& sd)
     force_(sd.force_),
     moment_(sd.moment_),
 
-    curTimeIndex_(-1),
+    curTimeIndex_(sd.curTimeIndex_),
     oldStatePtr_()
 {}
 

src/ODE/sixDOF/sixDOFODE/sixDOFODE.H

@@ -49,6 +49,7 @@ SourceFiles
 #include "OutputControlDictionary.H"
 #include "runTimeSelectionTables.H"
 #include "Switch.H"
+#include "foamTime.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -172,7 +173,10 @@ class sixDOFODE
 
             //- Initialise ODE before setting external forces/moments and
             //  solving
-            virtual void initState();
+            void initState();
+
+            //- Calculate current ODE step fraction given time from ODE
+            inline scalar odeStepFraction(const scalar odeTime) const;
 
 
 protected:
@@ -185,6 +189,17 @@ protected:
             virtual void setState(const sixDOFODE&);
 
 
+        // Get external force and moment (used during the solution process)
+
+            //- Return external force given time (linearly interpolated from
+            //  old state and current state)
+            dimensionedVector force(const scalar t) const;
+
+            //- Return external moment given time (linearly interpolated from
+            //  old state and current state)
+            dimensionedVector moment(const scalar t) const;
+
+
 public:
 
     //- Run-time type information
@@ -269,8 +284,8 @@ public:
             //- Set external force and moment
             inline void setExternalForceAndMoment
             (
-                const vector& externalForce,
+                const dimensionedVector& externalForce,
-                const vector& externalMoment
+                const dimensionedVector& externalMoment
             );
 
 

src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H

@@ -35,6 +35,23 @@ Author
 
 \*---------------------------------------------------------------------------*/
 
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
+
+Foam::scalar Foam::sixDOFODE::odeStepFraction(const scalar odeTime) const
+{
+    // Get current global time. Note: assuming that the ODESolver for a given
+    // time step is integrating from t - deltaT to t.
+    const scalar globalTime = dict().time().value();
+    const scalar globalDeltaT = dict().time().deltaT().value();
+
+    // Calculate current ODE step time step size
+    const scalar odeDeltaT = globalDeltaT - (globalTime - odeTime);
+
+    // Calculate and return ODE step size fraction
+    return odeDeltaT/globalDeltaT;
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::dimensionedScalar& Foam::sixDOFODE::mass() const
@@ -99,8 +116,8 @@ const Foam::dimensionedVector& Foam::sixDOFODE::moment() const
 
 void Foam::sixDOFODE::setExternalForceAndMoment
 (
-    const vector& externalForce,
+    const dimensionedVector& externalForce,
-    const vector& externalMoment
+    const dimensionedVector& externalMoment
 )
 {
     // Initialise the state before setting external forces and moments