src/rigidBodyDynamics/rigidBodyMotion: Added support for acceleration relaxation

src/rigidBodyDynamics/Make/files

@@ -4,6 +4,7 @@ bodies/jointBody/jointBody.C
 bodies/compositeBody/compositeBody.C
 bodies/subBody/subBody.C
 bodies/sphere/sphere.C
+bodies/cuboid/cuboid.C
 
 joints/joint/joint.C
 joints/null/nullJoint.C

src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.C

@@ -125,6 +125,19 @@ Foam::spatialTransform Foam::RBD::rigidBodyMotion::X00
 }
 
 
+void Foam::RBD::rigidBodyMotion::forwardDynamics
+(
+    rigidBodyModelState& state,
+    const scalarField& tau,
+    const Field<spatialVector>& fx
+) const
+{
+    scalarField qDdotPrev = state.qDdot();
+    rigidBodyModel::forwardDynamics(state, tau, fx);
+    state.qDdot() = aDamp_*(aRelax_*state.qDdot() + (1 - aRelax_)*qDdotPrev);
+}
+
+
 void Foam::RBD::rigidBodyMotion::solve
 (
     scalar deltaT,

src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H

@@ -158,6 +158,16 @@ public:
 
         // Update state
 
+            //- Calculate and optionally relax the joint acceleration qDdot from
+            //  the joint state q, velocity qDot, internal force tau (in the
+            //  joint frame) and external force fx (in the global frame)
+            void forwardDynamics
+            (
+                rigidBodyModelState& state,
+                const scalarField& tau,
+                const Field<spatialVector>& fx
+            ) const;
+
             //- Integrate velocities, orientation and position
             //  for the given time-step
             void solve

src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C

@@ -77,7 +77,7 @@ void Foam::RBD::rigidBodySolvers::CrankNicolson::solve
     model_.forwardDynamics(state(), tau, rfx);
 
     // Correct velocity
-    qDot() = qDot0() + aDamp()*deltaT()*(aoc_*qDdot() + (1 - aoc_)*qDdot0());
+    qDot() = qDot0() + deltaT()*(aoc_*qDdot() + (1 - aoc_)*qDdot0());
 
     // Correct position
     q() = q0() + deltaT()*(voc_*qDot() + (1 - voc_)*qDot0());

src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C

@@ -85,12 +85,12 @@ void Foam::RBD::rigidBodySolvers::Newmark::solve
 
     // Correct velocity
     qDot() = qDot0()
-      + aDamp()*deltaT()*(gamma_*qDdot() + (1 - gamma_)*qDdot0());
+      + deltaT()*(gamma_*qDdot() + (1 - gamma_)*qDdot0());
 
     // Correct position
     q() = q0()
       + deltaT()*qDot0()
-      + sqr(deltaT())*beta_*qDdot() + sqr(deltaT())*(0.5 - beta_)*qDdot0();
+      + sqr(deltaT())*(beta_*qDdot() + (0.5 - beta_)*qDdot0());
 
     correctQuaternionJoints();
 

src/rigidBodyDynamics/rigidBodySolvers/rigidBodySolver/rigidBodySolver.H

@@ -94,10 +94,6 @@ protected:
         //- Return the previous time-step
         inline scalar deltaT0() const;
 
-
-        //- Acceleration damping coefficient (for steady-state simulations)
-        inline scalar aDamp() const;
-
         //- Correct the quaternion joints based on the current change in q
         void correctQuaternionJoints();
 

src/rigidBodyDynamics/rigidBodySolvers/rigidBodySolver/rigidBodySolverI.H

@@ -85,10 +85,4 @@ inline Foam::scalar Foam::RBD::rigidBodySolver::deltaT0() const
 }
 
 
-inline Foam::scalar Foam::RBD::rigidBodySolver::aDamp() const
-{
-    return model_.aDamp_;
-}
-
-
 // ************************************************************************* //

src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C

@@ -70,7 +70,7 @@ void Foam::RBD::rigidBodySolvers::symplectic::solve
     // First simplectic step:
     //     Half-step for linear and angular velocities
     //     Update position and orientation
-    qDot() = qDot0() + aDamp()*0.5*deltaT0()*qDdot();
+    qDot() = qDot0() + 0.5*deltaT0()*qDdot();
     q() = q0() + deltaT()*qDot();
 
     correctQuaternionJoints();
@@ -88,7 +88,7 @@ void Foam::RBD::rigidBodySolvers::symplectic::solve
 
     // Second simplectic step:
     //     Complete update of linear and angular velocities
-    qDot() += aDamp()*0.5*deltaT()*qDdot();
+    qDot() += 0.5*deltaT()*qDdot();
 }
 
 

src/rigidBodyMeshMotion/rigidBodyMeshMotion.C

@@ -250,19 +250,16 @@ void Foam::rigidBodyMeshMotion::solve()
             const label bodyID = bodyMeshes_[bi].bodyID_;
 
             dictionary forcesDict;
-
             forcesDict.add("type", forces::typeName);
             forcesDict.add("patches", bodyMeshes_[bi].patches_);
             forcesDict.add("rhoInf", rhoInf_);
             forcesDict.add("rhoName", rhoName_);
-            forcesDict.add("CofR", model_.X0(bodyID).r());
+            forcesDict.add("CofR", vector::zero);
 
             forces f("forces", db(), forcesDict);
-
             f.calcForcesMoment();
 
-            fx[bodyID].l() = f.forceEff();
-            fx[bodyID].w() = f.momentEff();
+            fx[bodyID] = spatialVector(f.momentEff(), f.forceEff());
         }
 
         model_.solve

src/sixDoFRigidBodyMotion/sixDoFSolvers/Newmark/Newmark.C

@@ -96,7 +96,7 @@ void Foam::sixDoFSolvers::Newmark::solve
             tConstraints()
           & (
                 deltaT*v0()
-              + sqr(deltaT)*beta_*a() + sqr(deltaT)*(0.5 - beta_)*a0()
+              + aDamp()*sqr(deltaT)*(beta_*a() + (0.5 - beta_)*a0())
             )
         );
 
@@ -105,7 +105,7 @@ void Foam::sixDoFSolvers::Newmark::solve
         rConstraints()
       & (
             deltaT*pi0()
-          + sqr(deltaT)*beta_*tau() + sqr(deltaT)*(0.5 - beta_)*tau0()
+          + aDamp()*sqr(deltaT)*(beta_*tau() + (0.5 - beta_)*tau0())
         );
     Tuple2<tensor, vector> Qpi = rotate(Q0(), piDeltaT, 1);
     Q() = Qpi.first();