[mAlletto]

Hello all,


I set up a simple test case to test the 6dof solver in OF. The OF version is OF1912 It is a circular cylinder in a fluid at rest attached at a linear spring without damper. The density and viscosity of the fluid is set zero in order that no fluid force acts on the cylinder. After releasing the cylinder from an initial position where the spring is at tension, the amplitude of the oscillations increase instead to keep at a constant value like expected. The increase of the amplitude is smaller for smaller time steps but even for a very small time step of dt = 0.00005 the amplitude increases of 5% after ten cycles.



Any clue if I messed up something. Has anyone similar experiences?



I attached the case files with an Allrun script.

[mAlletto]

I just used the symplectic solver to update the position of the mass connected to the spring and I got the analytical solution for all time step sized used above (see plot).


Hm interesting. Have to check the implementation of the Newmark method.

[mAlletto]

Following wikipedia https://en.wikipedia.org/wiki/Newmark-beta_method the Newmark method is as follows (without external forces):










It is an implicit method. This means that we have three equations for the three unknowns: the position , the velocity and acceleration , at the new time step. Since the equations are linear we can rearrange the equation and find a solution.






I checked the code of the Newmark method implemented in openfoam:

void Foam::sixDoFSolvers::Newmark::solve
(
    bool firstIter,
    const vector& fGlobal,
    const vector& tauGlobal,
    scalar deltaT,
    scalar deltaT0
)
{
    // Update the linear acceleration and torque
    updateAcceleration(fGlobal, tauGlobal);

    // Correct linear velocity
    v() =
        tConstraints()
      & (v0() + aDamp()*deltaT*(gamma_*a() + (1 - gamma_)*a0()));

    // Correct angular momentum
    pi() =
        rConstraints()
      & (pi0() + aDamp()*deltaT*(gamma_*tau() + (1 - gamma_)*tau0()));

    // Correct position
    centreOfRotation() =
        centreOfRotation0()
      + (
            tConstraints()
          & (
                deltaT*v0()
              + aDamp()*sqr(deltaT)*(beta_*a() + (0.5 - beta_)*a0())
            )
        );

    // Correct orientation
    vector piDeltaT =
        rConstraints()
      & (
            deltaT*pi0()
          + aDamp()*sqr(deltaT)*(beta_*tau() + (0.5 - beta_)*tau0())
        );
    Tuple2<tensor, vector> Qpi = rotate(Q0(), piDeltaT, 1);
    Q() = Qpi.first();
}

For me it seems that the method is explicit since the acceleration is calculated before solving for the velocity and the new position. At this time only the values at the previous time step is known. So the system of equations solved in OF look more like this:










So maybe the calculation of the acceleration at the new time step using old velocity and position values lead to divergence.



Or maybe I'm wrong. Are there any comments from the community?



Best


Michael

[mAlletto]

Just wanted to say that I summarized some results in this small tutorial:


https://wiki.openfoam.com/1Dof_sprin...ichael_Alletto

[chennak]

The implementation of the Newmark scheme in Openfoam is incorrect. The implementation is explicit; however, it is not stable for gamma=0.5 and beta=0.25.