Commit ebf00110 authored by Henry Weller's avatar Henry Weller
Browse files

src/rigidBodyDynamics/rigidBodyMotion: Added support for acceleration relaxation

parent 51bb5093
......@@ -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
......
......@@ -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,
......
......@@ -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
......
......@@ -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());
......
......@@ -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();
......
......@@ -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();
......
......@@ -85,10 +85,4 @@ inline Foam::scalar Foam::RBD::rigidBodySolver::deltaT0() const
}
inline Foam::scalar Foam::RBD::rigidBodySolver::aDamp() const
{
return model_.aDamp_;
}
// ************************************************************************* //
......@@ -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();
}
......
......@@ -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
......
......@@ -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();
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment