View Issue Details

IDProjectCategoryView StatusLast Update
0001675OpenFOAMBugpublic2015-04-30 19:14
Reporteruser135Assigned Tohenry  
PrioritynormalSeveritymajorReproducibilityalways
Status closedResolutionunable to reproduce 
Summary0001675: sixDoFRigidBodyMotion: accelerationDamping is used inconsistently (well, actually wrong)
DescriptionIn six-DOF calculations, "accelerationDamping" should reduce the calculated accelerations by a certain factor. Instead, it reduces both velocity and acceleration.

In sixDoFRigidBodyMotion::updatePosition we read:

 v() = tConstraints_ & aDamp_*(v0() + 0.5*deltaT0*a());
pi() = rConstraints_ & aDamp_*(pi0() + 0.5*deltaT0*tau());

where it should be:

 v() = tConstraints_ & (v0() + aDamp_*0.5*deltaT0*a());
pi() = rConstraints_ & (pi0() + aDamp_*0.5*deltaT0*tau());

In the second step in sixDoFRigidBodyMotion::updateAcceleration it is correct:

v() += tConstraints_ & aDamp_*0.5*deltaT*a();
pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();

TagsNo tags attached.

Activities

henry

2015-04-28 13:17

manager   ~0004670

I tried several formulations to improve stability and convergence. Have you found that damping only the acceleration is preferable to damping both velocity and acceleration? For which cases and what was the benefit in terms of stability and convergence?

michele

2015-04-28 14:03

reporter   ~0004671

Last edited: 2015-04-28 14:06

Just my two cents, for a floating breakwater study, I found most useful to filter the forces instead of motions.
Here attached the modified sources from my 2.3.x installation.

[forceRelax.tgz]

michele

2015-04-28 14:04

reporter  

forceRelax.tgz (6,285 bytes)

henry

2015-04-28 14:13

manager   ~0004672

What is the difference in behavior between filtering the forces rather than the accelerations?

After a lot of messing around I found that relaxing the acceleration provided better behavior than damping it.

Have either of you tested the changes you are proposing on the tutorial cases:

multiphase/LTSInterFoam/DTCHull
multiphase/interDyMFoam/ras/DTCHull
multiphase/interDyMFoam/ras/floatingObject

user135

2015-04-28 14:16

  ~0004673

To be honest I didn't check if damping the acceleration is better than damping the velocities. But if it is called "accelerationDamping", it should damp the acceleration, not the velocity :). If needed, introduce "velocityDamping".

Currently any value for "accelerationDamping" smaller than 1 reduces the velocities further in each timestep . This is not what you would expect.

By the way: I think that an "accelerationLimiter" would be very useful to improve stability in critical cases.

user135

2015-04-28 14:25

  ~0004674

Just a comment:

Relaxing the accelerations is much better than damping, because the result is still physically correct for transient movements (if nOuterCorrectors is big enough). But for stationary systems, acceleration damping can be helpful to reach the steady state.

henry

2015-04-28 14:45

manager   ~0004675

I agree that the naming is not consistent with the action and I will remove the damping of the velocity unless the arises in which case I will add "velocityDamping".

How would you apply an "accelerationLimiter"? How would the limit be set?

user135

2015-04-28 15:09

  ~0004676

"accelerationLimiter" could be a scalar that is specified in dynamicMeshDict:sixDoFRigidBodyMotionCoeffs

In pseudo-code (sorry, don't know how to write this in OF-syntax), this might look like:

 if (magnitude(a()) > accelerationLimiter)
   a() *= accelerationLimiter/magnitude(a());

One could think about using a vector for accelerationLimiter and limit the components separately (indices are directions in space):

 if (abs(a(0)) > accelerationLimiter(0))
   a(0) *= accelerationLimiter/abs(a(0));
 if (abs(a(1)) > accelerationLimiter(1))
   a(1) *= accelerationLimiter/abs(a(1));
 if (abs(a(2)) > accelerationLimiter(2))
   a(2) *= accelerationLimiter/abs(a(2));

I think something similar would be useful for tau.

About setting the limit: In my usecases I have a good idea about the range of accelerations I can expect from a physical point of view. But that depends on the user ...

henry

2015-04-28 15:16

manager   ~0004677

The "accelerationDamping" issue is fixed by
commit fd0b7d87c61ef1261cf6f1eae58cc79349d01398

I am not convinced that an acceleration limiter in which the user must supply a vector of the maximum acceleration allowed would be generally useful. This would need some testing on a range of cases and more input from potential users of this functionality.

user135

2015-04-28 15:43

  ~0004678

Well, maybe you're right.

If I find some time, I'll try it on my cases ...

michele

2015-04-28 16:19

reporter   ~0004679

Based on my observations, the problem relies on pressure force' "spikes" that may appear during dynamic motions (so I'm doubtful of the usefulness of an acceleration limiter).

Directly filtering the forces is, in my opinion, most appropriate in the sense that the body motions are consistent with the "filtered forces" acting on it.

We may filter directly the acceleration (a=F/m), in this case the motion result should be identical (but if I post-process forces and accelerations, obviously the basic relationship a=F/m is no more preserved, because a is filtered and F is not).

In any case the aDamp_ parameter, as applied in the code, reduces unphysically the motions, so I preferred to adopt the solution as per the sources previously attached.

If I remember well, I made some verification testing in 2012, and find an improved accuracy (with respect to the "aDamp_ solution"). In the last years, I made various studies on moving structures and floating bodies, and find the proposed solution effective.

henry

2015-04-28 17:14

manager   ~0004680

The purpose of the aDamp_ parameter is indeed to "reduce unphysically the motions" to aid convergence towards a steady-state.

In your proposal to "filter the forces", what exactly do you do to the forces and why? Also how is this more physical than damping the acceleration?

michele

2015-04-28 17:52

reporter   ~0004681

As an example, let's consider the case of a floating breakwater.
The dynamic behaviour of the floater in waves is of primary importance (reflection, transmission and damping of waves is directly influenced by the body motions and interaction with waves). The aDamp_ solution, useful for reaching a steady state, is not the right strategy in this case (because alters unphysically the motions).

The filter on the forces is simply a first-order low-pass filter. It is characterised by a time-constant (that is related to the cut-off frequency of the filter).
So this filter only produces a smoothed force signal (with the drawback of a time-delay, that however is under user-control). The direct benefit is the absence of "force spikes", that in my experience, cause the system to became unstable.

Clearly the filter time constant should be a fraction of the dynamics of interest (otherwise it would alter significantly the structural dynamics).

A typical wave range of interest has a minimum period of about 3 seconds, and in my experience a typical filter period of 0.15 s is sufficient to remove spikes. As you can see, the ratio between the filter period and the dynamics of interest is, 0.15s/3s=0.05.
The filter characteristic (i.e. the ratio between the filtered signal amplitude and the original signal) at a frequency ratio of 0.05 is practically 1, so the filter does not interferes with the dynamical behaviour.

I just found these on the internet (but any signal processing text gives enough info)

http://en.wikipedia.org/wiki/Time_constant#Relation_of_time_constant_to_bandwidth

http://en.wikipedia.org/wiki/Time_constant#/media/File:Single-pole_frequency_response.JPG

henry

2015-04-28 18:16

manager   ~0004682

What is the advantage of the force filtering over accelerationRelaxation?

Could you provide a test case which demonstrates the advantage?

michele

2015-04-28 19:15

reporter   ~0004684

May be this example is not the best one, but is the simplest that comes to me.
Let's consider an unrestrained body, of mass m subjected to a constant force F (let's assume a step force variation from 0 at time=0). The body is initially still.
The analytical solution is a constant acceleration, with a=F/m, and v=a*time, for time>0.
Filtering the force guarantees that the solution is close to the analytical one. In particular the acceleration is identical to the analytical, after an initial transient of the filtered force (with a short duration, as per user-defined filter time constant).

Instead, as implemented, the acceleration is always a factor of the real one, i.e. a=aDamp_*F/m

As a consequence, dynamic motions are by definition wrong at any time instant, not only during an initial transient.

Remembering the objective of having a stable dynamic simulation behaviour, it can be concluded that filtering is a solution that has less impact on system's dynamics.

Same considerations may be performed in presence of sinusoidal forces and motions (as per typical cases of floaters in waves)

henry

2015-04-28 19:23

manager   ~0004685

Note that I am discussing accelerationRelaxation NOT accelerationDamping.

accelerationRelaxation does not affect the solution if the motion-pressure coupling is converged within a time-step. So far I have found that with accelerationRelaxation set to 0.3 good convergence behavior is obtained for the tutorial cases I listed earlier.

Can you upload a test-case which demonstrates the advantage of force filtering over accelerationRelaxation?

michele

2015-04-28 20:29

reporter   ~0004686

Ah, now I see. When I implemented the force filtering on OpenFoam 2.0, the accelerationRelaxation feature was not available ( https://github.com/OpenFOAM/OpenFOAM-2.0.x/blob/master/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C ).
During the openfoam evolution I just updated my changes to the newer version, as I previously said, my tests were made some time ago, indeed when the accelerationRelaxation was not available.
I just looked at the code of the 2.3.x, and, after a rapid glance, I think that, apart the small inconsistency between the forces (that are unfiltered in the accelerationRelaxation implementation) and the acceleration, the features look more or less equivalent (if time step is constant it is possible to define a correlation between the filter time constant and the relaxation factor used in your implementation).

user135

2015-04-30 12:57

  ~0004696

Sorry for opening again, but I want to continue the discussion thread:

Concerning a use case for an "accelerationLimiter":

We are doing simulations for seagoing vessels, simulating the process of being lifted upwards in a navigation lock. During this process, the sum of all forces on the vessel is near zero (compared to the mass force) and we are trying to compute the forces from the locking process on the vessel (which are very small, thus an accuracy of 0.00001 of the mass force is desired). When doing a restart of a parallel computation (~200 cores, with interDyMFoam), the computed forces on the hull are sometimes waaayyyy of in the first iteration of the first new timestep. Thus, a huge acceleration is computed (e.g. 200 m/s^2) and the whole pressure field jumps up and down for some time. Here, a limiter to e.g. 0.01 m/s^2 would be helpful. Or is it a bug in the restart mechanism?

henry

2015-04-30 13:02

manager   ~0004697

I would need to study the case and the terms for which you would need to arrange for some support.

Are you writing mesh and results in binary? In order to make the restart accurate you will either need to write binary or increase the write precision; binary is preferable.

user135

2015-04-30 13:20

  ~0004698

We tried with both binary and high accuracy ASCII. Further analysis has shown, that the rereading of the fields is not the problem: The force balance is fine directly after rereading. But after the the first p_rgh-computation cycle it is bad. We tried increasing accuracy in the p_rgh-solvers and increase nCorrectors. This makes it better, but doesn't cure it.

Sorry, can't make it a support case: It is relying on "special code", which should not leave the house :(

michele

2015-04-30 13:27

reporter   ~0004699

@carsten : you may find useful the force filtering approach, in the attached files. The restart would be smoothed by filtering directly the force (based on the latest time step).

user135

2015-04-30 13:34

  ~0004700

@Michele: Thanks for the suggestuion, but we already used relaxation. But it is not strong enough. Sometimes the forces are off for several orders of magnitude. In that situation, even the relaxed forces/accelerations are way too big.

michele

2015-04-30 13:48

reporter   ~0004701

There are more subtle differences between approaches, so my suggestion.

However, if I were you, I would better investigate (in detail) the restart process (better on a small case, both serial and parallel).

user135

2015-04-30 13:55

  ~0004703

@Michele: Sigh, I agree.

henry

2015-04-30 18:36

manager   ~0004707

In order to gon any further with this we would need a preferably small test-case which reproduces the problem.

Issue History

Date Modified Username Field Change
2015-04-28 12:58 user135 New Issue
2015-04-28 13:17 henry Note Added: 0004670
2015-04-28 14:03 michele Note Added: 0004671
2015-04-28 14:04 michele File Added: forceRelax.tgz
2015-04-28 14:06 michele Note Edited: 0004671
2015-04-28 14:13 henry Note Added: 0004672
2015-04-28 14:16 user135 Note Added: 0004673
2015-04-28 14:25 user135 Note Added: 0004674
2015-04-28 14:45 henry Note Added: 0004675
2015-04-28 15:09 user135 Note Added: 0004676
2015-04-28 15:16 henry Note Added: 0004677
2015-04-28 15:43 user135 Note Added: 0004678
2015-04-28 16:19 michele Note Added: 0004679
2015-04-28 17:14 henry Note Added: 0004680
2015-04-28 17:52 michele Note Added: 0004681
2015-04-28 18:16 henry Note Added: 0004682
2015-04-28 19:15 michele Note Added: 0004684
2015-04-28 19:23 henry Note Added: 0004685
2015-04-28 20:29 michele Note Added: 0004686
2015-04-28 20:38 henry Status new => resolved
2015-04-28 20:38 henry Resolution open => fixed
2015-04-28 20:38 henry Assigned To => henry
2015-04-30 12:57 user135 Note Added: 0004696
2015-04-30 12:57 user135 Status resolved => feedback
2015-04-30 12:57 user135 Resolution fixed => reopened
2015-04-30 13:02 henry Note Added: 0004697
2015-04-30 13:20 user135 Note Added: 0004698
2015-04-30 13:20 user135 Status feedback => assigned
2015-04-30 13:27 michele Note Added: 0004699
2015-04-30 13:34 user135 Note Added: 0004700
2015-04-30 13:48 michele Note Added: 0004701
2015-04-30 13:55 user135 Note Added: 0004703
2015-04-30 18:36 henry Note Added: 0004707
2015-04-30 18:36 henry Status assigned => closed
2015-04-30 18:37 henry Resolution reopened => unable to reproduce