From 577ae885a0c6c7f3b0e0babbc052fe4db45ac3ff Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Wed, 8 Mar 2017 11:08:20 +0100 Subject: [PATCH] Updates to sixDOFODE classes 2nd order accurate updates of force and moment during the ODE solution. --- .../sixDOF/geometricSixDOF/geometricSixDOF.C | 4 ++-- .../quaternionSixDOF/quaternionSixDOF.C | 4 ++-- src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C | 11 +++++++++- src/ODE/sixDOF/sixDOFODE/sixDOFODE.C | 22 ++++++++++++++++++- src/ODE/sixDOF/sixDOFODE/sixDOFODE.H | 21 +++++++++++++++--- src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H | 21 ++++++++++++++++-- 6 files changed, 72 insertions(+), 11 deletions(-) diff --git a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C b/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C index 7484f831d..fdf234216 100644 --- a/src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C +++ b/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(); diff --git a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C b/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C index 11599388c..d7162c705 100644 --- a/src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C +++ b/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(); diff --git a/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C b/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C index 3fa462aea..c7bdb1b31 100644 --- a/src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C +++ b/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() ); diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C index df3202d87..f50fd140d 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C +++ b/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_() {} diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H index e9fe063c2..1187ae262 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H +++ b/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 vector& externalMoment + const dimensionedVector& externalForce, + const dimensionedVector& externalMoment ); diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H index ed391a329..3a5c9618b 100644 --- a/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H +++ b/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 vector& externalMoment + const dimensionedVector& externalForce, + const dimensionedVector& externalMoment ) { // Initialise the state before setting external forces and moments