View Issue Details

IDProjectCategoryView StatusLast Update
0002297OpenFOAMBugpublic2016-10-19 15:13
Reportersose Assigned Tohenry  
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version15.04
Fixed in Versiondev 
Summary0002297: 6DOF solvers do not capture the nonlinearity in the rotational DOFs correctly
DescriptionWe are simulating a torque free motion where an analytical solution exists. A rectangular body is floating freely in vacuum free of any external force/moment. Initially the angular velocity is set to omega(x,y,z) = {0, 10, 0.05} (body frame). Due to the nonlinearity in the rotational DOF the small rotation in z-axis will induce a periodic flip motion. These are the results we g
et using the available solvers in vers. 3 and 4

- sixDofRigidBodyMotion (symplectic): the nonlinear motion is captured correctly
- sixDofRigidBodyMotion (CrankNicolson): the flip motion is missing.
- rigidBodyModel (symplectic/CrankNicolson): there are artificial oscillations in the angular velocity and the amplitude is increasing steadily
Steps To Reproduceunpack the attached tar-file
execute Allclean/Allrun

run the attached matlab/octave script (main.m) to show the results
TagsNo tags attached.



2016-10-18 11:04


testing6DOF.tgz (133,231 bytes)


2016-10-18 11:22

manager   ~0007020

Have you tested the Newmark scheme for the rigidBodyModel?


2016-10-18 11:37

reporter   ~0007021

With the default parameters Newmark and CrankNicolson are exactly the same scheme. I have run the rigidBodyModel with the Newmark scheme. I see the same results as for the CN.


2016-10-18 12:01

manager   ~0007022

- sixDofRigidBodyMotion (symplectic): the nonlinear motion is captured correctly

What are the problems with this method? I would expect it to be most appropriate for your case.


2016-10-18 12:54

reporter   ~0007023

Thanks for your effort on this. Yes, the symplectic scheme when applied with sixDofRigidBodyMotion can capture the nonlinear motion correctly in this test case. However, I believe you must have a reason for introducing the CN/Newmark scheme in the package. If you think there are bugs in these solvers, we should fix them. Here the symplectic scheme works OK for the sixDofRigidBodyMotion but not for the rigidBodyModel.


2016-10-18 13:01

manager   ~0007024

I do not think there are bugs in the CN/Newmark schemes but if you believe there are and can provide patches to fix them I would be happy to apply them.

I am not sure why the symplectic scheme in the rigidBodyModel does not behave the same as for the sixDofRigidBodyMotion, could you investigate and provide any changes necessary?


2016-10-18 16:44

reporter   ~0007029

I agree that the problem with CN/Newmark schemes can not be fixed straight away. The symplectic scheme in the sixDofRigidMotion has problems with stability when the coupling between the flow solution and the motion is strong. This tiny modification (pls. see the patch attached to this note) should help to improve the stability of the symplectic scheme.
symplectic.patch (2,086 bytes)   
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
index a91f17e..2bc1cde 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
@@ -269,6 +269,9 @@ public:
             //- Store the motion state at the beginning of the time-step
             inline void newTime();
+            //- Reset the motion state after the first iteration
+            inline void resetTime();
             //- Return non-const access to the centre of rotation
             inline point& centreOfRotation();
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
index ffe8981..5d68c6d 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
@@ -274,6 +274,10 @@ inline void Foam::sixDoFRigidBodyMotion::newTime()
     motionState0_ = motionState_;
+inline void Foam::sixDoFRigidBodyMotion::resetTime()
+    motionState_ = motionState0_;
 inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfRotation()
diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C b/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C
index 61d94ab..b34e1e5 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C
@@ -69,10 +69,12 @@ void Foam::sixDoFSolvers::symplectic::solve
     if (!firstIter)
-        FatalErrorIn("sixDoFSolvers::symplectic::solve")
+/*        FatalErrorIn("sixDoFSolvers::symplectic::solve")
             << "The symplectic integrator is explicit "
                "and can only be solved once per time-step"
             << exit(FatalError);
+        body_.resetTime();
     // First simplectic step:
symplectic.patch (2,086 bytes)   


2016-10-18 16:53

manager   ~0007030

> I agree that the problem with CN/Newmark schemes can not be fixed straight away.

I do not see anything wrong with these schemes or that they need to be "fixed" in anyway. I think they are just not appropriate for your case unless you use a very small time-step (this is conjecture).

> This tiny modification (pls. see the patch attached to this note) should help to improve the stability of the symplectic scheme.

Does it help? It appears that you have enabled the use of the symplectic scheme with outer-iterations which I am not sure is consistent with the formulation of the scheme.


2016-10-19 11:39

reporter   ~0007039

From what I have seen in my tests I will never use CN/Newmark in a case where the nonlinear interaction between the rotational dof is expected. For the symplectic scheme I have overlooked the existent of the underrelaxation in the acceleration. Indeed it is not consistent to reset the 6dof solution due to the underrelaxation. However, simply allow the outer-iterations with the resetting will do just fine because the vel. & pos. at the previous time step will not be updated until a new call to newTime(). The following is the algorithm for the symplectic scheme (without the translational dof)

    //1. compute ang.vel. at haft step using acc. at 0
    pi() = pi0() + 0.5*deltaT0*tau0();

    //2. compute orientation at full time step using the haft step vel.
    Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
    pi() = Qpi.second();

    //3. compute acc. at full time step (this is where the outer-iter. makes a difference)
    updateAcceleration(fGlobal, tauGlobal);

    //4. compute vel. at full time step
    pi() += 0.5*deltaT*tau();

In steps 1 & 2 the solution will be recomputed using the previous state (0)
In steps 3 & 4 the outer-iteration will make sure the the acceleration is consistently computed from the flow solution at a full time step.

In my tests this tiny modification is making a huge difference. Attached is the results of one of the cases where the outer-iteration really helps !!!

Btw, why is there aDamp() in the algorithm? I cannot imaging why anyone would want to reduce/amplify the acceleration by a constant amount like this.
symplectic-2iters.png (29,606 bytes)   
symplectic-2iters.png (29,606 bytes)   


2016-10-19 12:15

manager   ~0007042

Thanks for the analysis, I will apply the change.

> Btw, why is there aDamp() in the algorithm?

Useful for stabilizing runs during the early stages of motion development or to aid convergence to a steady-state. However, I would not choose the symplectic scheme for steady-state runs, the Newmark scheme is more robust, stable and convergent.


2016-10-19 12:51

reporter   ~0007043

sorry, i've noticed a critical typos in my previous post. it should have been

... allow the outer-iterations "without" the resetting ...

and thank you for considering pushing this change into the code.


2016-10-19 13:06

manager   ~0007044

The algorithm you propose only updates pi() and is missing the application of the constraints so I cannot apply it. Could you clarify which part of the algorithm you have changes because apart from the missing functionality it appears to be the same.


2016-10-19 13:24

reporter   ~0007045

Please see the patch below.

diff --git a/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C b/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C
--- a/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C
@@ -67,14 +67,6 @@ void Foam::sixDoFSolvers::symplectic::solve
     scalar deltaT0
- if (!firstIter)
- {
- FatalErrorIn("sixDoFSolvers::symplectic::solve")
- << "The symplectic integrator is explicit "
- "and can only be solved once per time-step"
- << exit(FatalError);
- }
     // First simplectic step:
     // Half-step for linear and angular velocities
     // Update position and orientation


2016-10-19 13:40

manager   ~0007046

That just removes the check for the use of symplectic scheme with outer iteration; is that sufficient? Is no other change necessary?


2016-10-19 13:53

reporter   ~0007047

yes, simply remove this check and allow the outer iteration is sufficient c.f. note ~0007039


2016-10-19 15:13

manager   ~0007048

Resolved in OpenFOAM-dev by commit 7ae0766a311fd09377ce97bcc83d73915cd9274c

Issue History

Date Modified Username Field Change
2016-10-18 11:04 sose New Issue
2016-10-18 11:04 sose File Added: testing6DOF.tgz
2016-10-18 11:22 henry Note Added: 0007020
2016-10-18 11:37 sose Note Added: 0007021
2016-10-18 12:01 henry Note Added: 0007022
2016-10-18 12:54 sose Note Added: 0007023
2016-10-18 13:01 henry Note Added: 0007024
2016-10-18 16:44 sose File Added: symplectic.patch
2016-10-18 16:44 sose Note Added: 0007029
2016-10-18 16:53 henry Note Added: 0007030
2016-10-19 11:39 sose File Added: symplectic-2iters.png
2016-10-19 11:39 sose Note Added: 0007039
2016-10-19 12:15 henry Note Added: 0007042
2016-10-19 12:51 sose Note Added: 0007043
2016-10-19 13:06 henry Note Added: 0007044
2016-10-19 13:24 sose Note Added: 0007045
2016-10-19 13:40 henry Note Added: 0007046
2016-10-19 13:53 sose Note Added: 0007047
2016-10-19 15:13 henry Assigned To => henry
2016-10-19 15:13 henry Status new => resolved
2016-10-19 15:13 henry Resolution open => fixed
2016-10-19 15:13 henry Fixed in Version => dev
2016-10-19 15:13 henry Note Added: 0007048