From 4de19563e90d384c42e2fb17c47f285f5edd154b Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Fri, 23 Oct 2015 13:43:35 +0200 Subject: [PATCH 1/2] Motion constraints in sixDOFqODE --- .../sixDOF/finiteRotation/finiteRotation.C | 8 +- .../sixDOF/finiteRotation/finiteRotation.H | 8 +- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C | 87 ++++++++++++++++++- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H | 21 +++++ .../basic/sixDOFSolver/sixDOFmotion/0/bodyA | 8 ++ .../basic/sixDOFSolver/sixDOFmotion/0/bodyB | 8 ++ .../springDamper/0/massSpringDamper | 9 ++ 7 files changed, 138 insertions(+), 11 deletions(-) diff --git a/src/ODE/sixDOF/finiteRotation/finiteRotation.C b/src/ODE/sixDOF/finiteRotation/finiteRotation.C index 501687de9..c03d0d037 100644 --- a/src/ODE/sixDOF/finiteRotation/finiteRotation.C +++ b/src/ODE/sixDOF/finiteRotation/finiteRotation.C @@ -1,9 +1,9 @@ /*---------------------------------------------------------------------------*\ ========= | \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 3.2 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of foam-extend. @@ -72,7 +72,7 @@ Foam::scalar Foam::finiteRotation::rotAngle(const tensor& rotT) Foam::vector Foam::finiteRotation::eulerAngles(const tensor& rotT) { - // Define a vector containing euler angles (x = roll, y = pitch, z = yaw) + // Create a vector containing euler angles (x = roll, y = pitch, z = yaw) vector eulerAngles; scalar& rollAngle = eulerAngles.x(); diff --git a/src/ODE/sixDOF/finiteRotation/finiteRotation.H b/src/ODE/sixDOF/finiteRotation/finiteRotation.H index 9b46b0191..2ea2bf884 100644 --- a/src/ODE/sixDOF/finiteRotation/finiteRotation.H +++ b/src/ODE/sixDOF/finiteRotation/finiteRotation.H @@ -1,9 +1,9 @@ /*---------------------------------------------------------------------------*\ ========= | \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 3.2 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright + \\ / O peration | + \\ / A nd | For copyright notice see file Copyright + \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of foam-extend. @@ -69,7 +69,7 @@ class finiteRotation //- Calculate rotation angle from given rotation tensor static scalar rotAngle(const tensor& rotT); - //- Calculate Euler angles (x-y-z (roll-pitch-yaw) convenction) from + //- Calculate Euler angles (x-y-z (roll-pitch-yaw) convention) from // given rotation tensor. Reference: Mike Day, Insomniac Games, // Extracting Euler Angles from a Rotation Matrix. static vector eulerAngles(const tensor& rotT); diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C index bed3ffc8d..5adbf0141 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C @@ -118,11 +118,14 @@ Foam::dimensionedVector Foam::sixDOFqODE::E Foam::sixDOFqODE::sixDOFqODE(const IOobject& io) : IOdictionary(io), + mass_(lookup("mass")), momentOfInertia_(lookup("momentOfInertia")), + Xequilibrium_(lookup("equilibriumPosition")), linSpringCoeffs_(lookup("linearSpring")), linDampingCoeffs_(lookup("linearDamping")), + Xrel_(lookup("Xrel")), U_(lookup("U")), Uaverage_(U_), @@ -134,11 +137,20 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io) omega_(lookup("omega")), omegaAverage_(omega_), omegaAverageAbsolute_(omega_), + force_(lookup("force")), moment_(lookup("moment")), forceRelative_(lookup("forceRelative")), momentRelative_(lookup("momentRelative")), - coeffs_(13, 0.0) + + coeffs_(13, 0.0), + + fixedSurge_(lookup("fixedSurge")), + fixedSway_(lookup("fixedSway")), + fixedHeave_(lookup("fixedHeave")), + fixedRoll_(lookup("fixedRoll")), + fixedPitch_(lookup("fixedPitch")), + fixedYaw_(lookup("fixedYaw")) { setCoeffs(); } @@ -163,11 +175,14 @@ Foam::sixDOFqODE::sixDOFqODE IOobject::NO_WRITE ) ), + mass_(sd.mass_.name(), sd.mass_), momentOfInertia_(sd.momentOfInertia_.name(), sd.momentOfInertia_), + Xequilibrium_(sd.Xequilibrium_.name(), sd.Xequilibrium_), linSpringCoeffs_(sd.linSpringCoeffs_.name(), sd.linSpringCoeffs_), linDampingCoeffs_(sd.linDampingCoeffs_.name(), sd.linDampingCoeffs_), + Xrel_(sd.Xrel_.name(), sd.Xrel_), U_(sd.U_.name(), sd.U_), Uaverage_(sd.Uaverage_.name(), sd.Uaverage_), @@ -179,11 +194,20 @@ Foam::sixDOFqODE::sixDOFqODE sd.omegaAverageAbsolute_.name(), sd.omegaAverageAbsolute_ ), + force_(sd.force_.name(), sd.force_), moment_(sd.moment_.name(), sd.moment_), forceRelative_(sd.forceRelative_.name(), sd.forceRelative_), momentRelative_(sd.momentRelative_.name(), sd.momentRelative_), - coeffs_(sd.coeffs_) + + coeffs_(sd.coeffs_), + + fixedSurge_(sd.fixedSurge_), + fixedSway_(sd.fixedSway_), + fixedHeave_(sd.fixedHeave_), + fixedRoll_(sd.fixedRoll_), + fixedPitch_(sd.fixedPitch_), + fixedYaw_(sd.fixedYaw_) {} @@ -199,9 +223,11 @@ void Foam::sixDOFqODE::setState(const sixDOFqODE& sd) { mass_ = sd.mass_; momentOfInertia_ = sd.momentOfInertia_; + Xequilibrium_ = sd.Xequilibrium_; linSpringCoeffs_ = sd.linSpringCoeffs_; linDampingCoeffs_ = sd.linDampingCoeffs_; + Xrel_ = sd.Xrel_; U_ = sd.U_; Uaverage_ = sd.Uaverage_; @@ -209,6 +235,7 @@ void Foam::sixDOFqODE::setState(const sixDOFqODE& sd) omega_ = sd.omega_; omegaAverage_ = sd.omegaAverage_; omegaAverageAbsolute_ = sd.omegaAverageAbsolute_; + force_ = sd.force_; moment_ = sd.moment_; forceRelative_ = sd.forceRelative_; @@ -217,6 +244,13 @@ void Foam::sixDOFqODE::setState(const sixDOFqODE& sd) // Copy ODE coefficients: this carries actual ODE state // HJ, 23/Mar/2015 coeffs_ = sd.coeffs_; + + fixedSurge_ = sd.fixedSurge_; + fixedSway_ = sd.fixedSway_; + fixedHeave_ = sd.fixedHeave_; + fixedRoll_ = sd.fixedRoll_; + fixedPitch_ = sd.fixedPitch_; + fixedYaw_ = sd.fixedYaw_; } @@ -248,6 +282,23 @@ void Foam::sixDOFqODE::derivatives dydx[4] = accel.y(); dydx[5] = accel.z(); + // Add translational constraints by setting RHS of given components to zero + if (fixedSurge_) + { + dydx[0] = 0; // Surge velocity + dydx[3] = 0; // Surge acceleration + } + if (fixedSway_) + { + dydx[1] = 0; // Sway velocity + dydx[4] = 0; // Sway acceleration + } + if (fixedHeave_) + { + dydx[2] = 0; // Heave velocity + dydx[5] = 0; // Heave acceleration + } + // Set the derivatives for rotation dimensionedVector curOmega ( @@ -258,7 +309,7 @@ void Foam::sixDOFqODE::derivatives const vector omegaDot = OmegaDot(curRotation, curOmega).value(); - dydx[6] = omegaDot.x(); + dydx[6] = omegaDot.x(); dydx[7] = omegaDot.y(); dydx[8] = omegaDot.z(); @@ -266,6 +317,23 @@ void Foam::sixDOFqODE::derivatives dydx[10] = curRotation.eDot(curOmega.value(), 1); dydx[11] = curRotation.eDot(curOmega.value(), 2); dydx[12] = curRotation.eDot(curOmega.value(), 3); + + // Add rotational constraints by setting RHS of given components to zero + if (fixedRoll_) + { + dydx[6] = 0; // Roll angular acceleration + dydx[10] = 0; // Roll angular velocity (roll quaternion evolution RHS) + } + if (fixedPitch_) + { + dydx[7] = 0; // Pitch angular acceleration + dydx[11] = 0; // Pitch angular velocity (pitch quaternion evolution RHS) + } + if (fixedYaw_) + { + dydx[8] = 0; // Yaw angular acceleration + dydx[12] = 0; // Yaw angular velocity (yaw quaternion evolution RHS) + } } @@ -367,6 +435,19 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const sixDOFqODE& sds) os.writeKeyword("momentRelative") << tab << sds.momentRelative() << token::END_STATEMENT << endl; + os.writeKeyword("fixedSurge") << tab << sds.fixedSurge_ << + token::END_STATEMENT << endl; + os.writeKeyword("fixedSway") << tab << sds.fixedSway_ << + token::END_STATEMENT << endl; + os.writeKeyword("fixedHeave") << tab << sds.fixedHeave_ << + token::END_STATEMENT << endl; + os.writeKeyword("fixedRoll") << tab << sds.fixedRoll_ << + token::END_STATEMENT << endl; + os.writeKeyword("fixedPitch") << tab << sds.fixedPitch_ << + token::END_STATEMENT << endl; + os.writeKeyword("fixedYaw") << tab << sds.fixedYaw_ << + token::END_STATEMENT << endl; + return os; } diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H index 89023e608..6b9eb3d4d 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.H @@ -131,6 +131,27 @@ class sixDOFqODE scalarField coeffs_; + //- Motion constraints (given as fixed motion components) + + //- Fixed surge (x-translation) + Switch fixedSurge_; + + //- Fixed sway (y-translation) + Switch fixedSway_; + + //- Fixed heave (z-translation) + Switch fixedHeave_; + + //- Fixed roll (rotation around x) + Switch fixedRoll_; + + //- Fixed pitch (rotation around y) + Switch fixedPitch_; + + //- Fixed yaw (rotation around z) + Switch fixedYaw_; + + // Private Member Functions //- Disallow default bitwise copy construct diff --git a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA index edf2748d2..3331eeae0 100644 --- a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA +++ b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyA @@ -36,4 +36,12 @@ moment m [1 2 -2 0 0 0 0] (0 0 0); forceRelative fRel [1 1 -2 0 0 0 0] (0 0 0); momentRelative mRel [1 2 -2 0 0 0 0] (-0.5773502692 -0.5773502692 -0.5773502692); +// Motion constraints +fixedSurge no; +fixedSway no; +fixedHeave no; +fixedRoll no; +fixedPitch no; +fixedYaw no; + // ************************************************************************* // diff --git a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB index 11e2a69fb..4865643ed 100644 --- a/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB +++ b/tutorials/basic/sixDOFSolver/sixDOFmotion/0/bodyB @@ -35,4 +35,12 @@ moment m [1 2 -2 0 0 0 0] (0 0 0); forceRelative fRel [1 1 -2 0 0 0 0] (0 0 0); momentRelative mRel [1 2 -2 0 0 0 0] (0 0 10); +// Motion constraints +fixedSurge no; +fixedSway no; +fixedHeave no; +fixedRoll no; +fixedPitch no; +fixedYaw no; + // ************************************************************************* // diff --git a/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper b/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper index 89ce8f0c2..2cd557d0c 100644 --- a/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper +++ b/tutorials/basic/sixDOFSolver/springDamper/0/massSpringDamper @@ -35,4 +35,13 @@ moment m [1 2 -2 0 0 0 0] (0 0 0); forceRelative fRel [1 1 -2 0 0 0 0] (0 0 0); momentRelative mRel [1 2 -2 0 0 0 0] (0 0 0); + +// Motion constraints +fixedSurge no; +fixedSway no; +fixedHeave no; +fixedRoll no; +fixedPitch no; +fixedYaw no; + // ************************************************************************* // From d9c65ce85f7968996e3b4b317aa4f07339594ad4 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Thu, 15 Oct 2015 08:01:29 +0200 Subject: [PATCH 2/2] Constraining omegaAverageAbsolute after calculation --- src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C | 46 +++++++++++++++++++++----- 1 file changed, 38 insertions(+), 8 deletions(-) diff --git a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C index 5adbf0141..62757ffb1 100644 --- a/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C +++ b/src/ODE/sixDOF/sixDOFqODE/sixDOFqODE.C @@ -91,12 +91,29 @@ Foam::dimensionedVector Foam::sixDOFqODE::OmegaDot const dimensionedVector& omega ) const { + // Fix the global moment for global rotation constraints + dimensionedVector mAbs = moment(); + vector& mAbsVal = mAbs.value(); + + if (fixedRoll_) + { + mAbsVal.x() = 0; + } + if (fixedPitch_) + { + mAbsVal.y() = 0; + } + if (fixedYaw_) + { + mAbsVal.z() = 0; + } + return inv(momentOfInertia_) & ( E(omega) // To relative - + (rotation.R() & moment()) + + (rotation.R() & mAbs) + momentRelative() ); } @@ -321,18 +338,15 @@ void Foam::sixDOFqODE::derivatives // Add rotational constraints by setting RHS of given components to zero if (fixedRoll_) { - dydx[6] = 0; // Roll angular acceleration - dydx[10] = 0; // Roll angular velocity (roll quaternion evolution RHS) + dydx[10] = 0; // Roll axis (roll quaternion evolution RHS) } if (fixedPitch_) { - dydx[7] = 0; // Pitch angular acceleration - dydx[11] = 0; // Pitch angular velocity (pitch quaternion evolution RHS) + dydx[11] = 0; // Pitch axis (pitch quaternion evolution RHS) } if (fixedYaw_) { - dydx[8] = 0; // Yaw angular acceleration - dydx[12] = 0; // Yaw angular velocity (yaw quaternion evolution RHS) + dydx[12] = 0; // Yaw axis (yaw quaternion evolution RHS) } } @@ -389,7 +403,23 @@ void Foam::sixDOFqODE::update(const scalar delta) ); omegaAverage_.value() = rotation_.omegaAverage(delta); - omegaAverageAbsolute_.value() = rotation_.omegaAverageAbsolute(delta); + + // Calculate and constrain omegaAverageAbsolute appropriately + vector& omegaAverageAbsoluteValue = omegaAverageAbsolute_.value(); + omegaAverageAbsoluteValue = rotation_.omegaAverageAbsolute(delta); + + if (fixedRoll_) + { + omegaAverageAbsoluteValue.x() = 0; + } + if (fixedPitch_) + { + omegaAverageAbsoluteValue.y() = 0; + } + if (fixedYaw_) + { + omegaAverageAbsoluteValue.z() = 0; + } }