diff --git a/src/rigidBodyDynamics/rigidBodyModelState/rigidBodyModelState.C b/src/rigidBodyDynamics/rigidBodyModelState/rigidBodyModelState.C index cad3a5ce04a837e280acba4b37d16c9727dd3f9d..33ea779e7fccd958eeb4fe7fb8130ae5801e22b5 100644 --- a/src/rigidBodyDynamics/rigidBodyModelState/rigidBodyModelState.C +++ b/src/rigidBodyDynamics/rigidBodyModelState/rigidBodyModelState.C @@ -33,7 +33,7 @@ Foam::RBD::rigidBodyModelState::rigidBodyModelState ) : q_(model.nDoF(), Zero), - w_(model.nw(), Zero), + w_(model.nw(), 1), qDot_(model.nDoF(), Zero), qDdot_(model.nDoF(), Zero) {} @@ -46,7 +46,7 @@ Foam::RBD::rigidBodyModelState::rigidBodyModelState ) : q_(dict.lookupOrDefault("q", scalarField(model.nDoF(), Zero))), - w_(dict.lookupOrDefault("w", scalarField(model.nw(), Zero))), + w_(dict.lookupOrDefault("w", scalarField(model.nw(), 1))), qDot_(dict.lookupOrDefault("qDot", scalarField(model.nDoF(), Zero))), qDdot_(dict.lookupOrDefault("qDdot", scalarField(model.nDoF(), Zero))) {} diff --git a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H index 398528a3b6d3edf7c6abadc203d125f11d0dce4a..b2e4562b80f246f6fc6cbdb40d8355323ce89e8a 100644 --- a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H +++ b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H @@ -137,6 +137,9 @@ public: //- Return the motion state inline const rigidBodyModelState& state() const; + //- Return the motion state for modification + inline rigidBodyModelState& state(); + // Edit diff --git a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotionI.H b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotionI.H index e9a430320df7ec4e1f2fe2dc4daaad49a13792ec..0de439f807ebabe43b6a1a956d0f320faabb3ef0 100644 --- a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotionI.H +++ b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotionI.H @@ -38,6 +38,13 @@ Foam::RBD::rigidBodyMotion::state() const } +inline Foam::RBD::rigidBodyModelState& +Foam::RBD::rigidBodyMotion::state() +{ + return motionState_; +} + + inline void Foam::RBD::rigidBodyMotion::newTime() { motionState0_ = motionState_; diff --git a/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C b/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C index 895e9b9e6102daad0249799a4e97ab7b4ac898a4..6d48971be65778e8ce0663b97a2d2f16d1972e69 100644 --- a/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C +++ b/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C @@ -82,6 +82,8 @@ void Foam::RBD::rigidBodySolvers::CrankNicolson::solve // Correct position q() = q0() + deltaT()*(voc_*qDot() + (1 - voc_)*qDot0()); + correctQuaternionJoints(); + // Update the body-state model_.forwardDynamicsCorrection(state()); } diff --git a/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C b/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C index 72e2eecea4cb19253780b287bdd1fd9901d253b1..2d5e5d90eaf3c77d037ebd3b17c53edc0f9548a8 100644 --- a/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C +++ b/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C @@ -92,6 +92,8 @@ void Foam::RBD::rigidBodySolvers::Newmark::solve + deltaT()*qDot0() + sqr(deltaT())*beta_*qDdot() + sqr(deltaT())*(0.5 - beta_)*qDdot0(); + correctQuaternionJoints(); + // Update the body-state model_.forwardDynamicsCorrection(state()); } diff --git a/src/rigidBodyDynamics/rigidBodySolvers/rigidBodySolver/rigidBodySolver.C b/src/rigidBodyDynamics/rigidBodySolvers/rigidBodySolver/rigidBodySolver.C index 669d8c9cc355267dd340d5b7d6fbaa05e6a45292..7cfc931778e0428ce8bdb1dce3ba98d8819a912c 100644 --- a/src/rigidBodyDynamics/rigidBodySolvers/rigidBodySolver/rigidBodySolver.C +++ b/src/rigidBodyDynamics/rigidBodySolvers/rigidBodySolver/rigidBodySolver.C @@ -52,4 +52,40 @@ Foam::RBD::rigidBodySolver::~rigidBodySolver() {} +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::RBD::rigidBodySolver::correctQuaternionJoints() +{ + if (model_.nw()) + { + forAll (model_.joints(), i) + { + const label qi = model_.joints()[i].qIndex(); + + if (model_.joints()[i].nw() == 1) + { + // Calculate the change in the normalized quaternion axis + vector dv((q().block<vector>(qi) - q0().block<vector>(qi))); + scalar magDv = mag(dv); + + if (magDv > SMALL) + { + // Calculate the quaternion corresponding to the change + quaternion dQuat(dv/magDv, cos(magDv), true); + + // Transform the previous time quaternion + quaternion quat + ( + normalize(model_.joints()[i](q0(), w0())*dQuat) + ); + + // Update the joint state + model_.joints()[i](quat, q(), w()); + } + } + } + } +} + + // ************************************************************************* // diff --git a/src/rigidBodyDynamics/rigidBodySolvers/rigidBodySolver/rigidBodySolver.H b/src/rigidBodyDynamics/rigidBodySolvers/rigidBodySolver/rigidBodySolver.H index a1959991b1132104e4c7ccb5ef05317221944a6f..0183a4be7162ab0015a5e21cabd8214d7ec1bd6c 100644 --- a/src/rigidBodyDynamics/rigidBodySolvers/rigidBodySolver/rigidBodySolver.H +++ b/src/rigidBodyDynamics/rigidBodySolvers/rigidBodySolver/rigidBodySolver.H @@ -100,6 +100,9 @@ protected: //- Acceleration damping coefficient (for steady-state simulations) inline scalar aDamp() const; + //- Correct the quaternion joints based on the current change in q + void correctQuaternionJoints(); + public: diff --git a/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C b/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C index 76807a32b22719d461764e05f3e8ea3f3cdd8ef5..18b78ad89f29092f820e4c271fc17cad7ea7e14b 100644 --- a/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C +++ b/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C @@ -73,6 +73,8 @@ void Foam::RBD::rigidBodySolvers::symplectic::solve qDot() = qDot0() + aDamp()*0.5*deltaT0()*qDdot(); q() = q0() + deltaT()*qDot(); + correctQuaternionJoints(); + // Update the body-state prior to the evaluation of the restraints model_.forwardDynamicsCorrection(state());