Constrained Dynamics (II)

This post is a continuation to the first entry in this series:

• Part I: Constrained Dynamics (I) - Unconstrained dynamics.

Constrained dynamics

We left off with a single unconstrained particle in free-fall.

In order to build more interesting setups we need to define rules that constrain the motion of the particle in a certain desired way. For example, a rigid-arm pendulum could be modeled with a particle mass plus a constraint to keep the particle at all times at exactly the same distance from the origin.

Example: Distance constraint

This will be our goal in this post: to simulate the motion of a particle that is subject to gravity, but is also constrained to keep a given distance to the origin at all times. Such type of constraint is called a distance constraint.

In all the examples here we will make the particle start resting at arm’s length to the right side of the origin. So in each of the animations below, we will be expecting a yellow circle arc plot like this one. Any yellow plot that deviates from a circle arc will reveal that the distance constraint is not being successfully satisfied.

Springs as constraints (don’t!)

We might attempt to model such constraint with a sufficiently strong (mass-less) spring connecting the particle to the origin. The spring would allow the particle to orbit freely, but forbid it from abandoning the spring-long orbit.

While this setup should work in principle, the spring would require a very large stiffness constant in order to not feel “goopy”. The spring would apply a corrective force to the particle as soon as it drifts off orbit. But this would eventually overshoot, and then rectify back, and so on, ending up in oscillatory motion. The stiffer the spring, the more acute the micro-oscillations. This is a recipe for disaster, specially as soon as the setup becomes more complex and more masses and springs are involved.

Springs

The type of coiled metal spring that we’re all familiar with, remains at some rest length until you try to compress or extend it.

If you try to compress the spring, it will resist your force, trying to extend back to its rest length. If you try to extend the spring, it will also resist your force, trying to compress back to its rest length.

Hooke’s law states that the force needed to compress/extend an ideal spring by some distance \(\Delta{L}=(L_t-L)\) above/below its rest length \(L\) scales linearly with respect to said length increment:

The source code below extends the previous post in order to incorporate a spring connecting the particle to the origin alongside the force of gravity.

struct spring_mass_t
{
  explicit spring_mass_t
  
      ( const f64_t rest_length = .5,
        const f64_t stiffness   = 1e3,
        const f64_t damping     = 1e3 )

      : m( 100 ), p( rest_length, 0 ), v( 0 ),

        L( rest_length ), ks( stiffness ), kd( damping ) {}

  void step( const f64_t dt, const vec2d_c& gravity )
  {
    // Current spring state.

    const vec2d_c s_D  =   p;             // Spring dir. Anchored at (0,0).

    const   f64_t s_L  =   s_D.length();  // Spring len.

    const vec2d_c s_d  = ( s_D / s_L );   // Spring dir (normalized).

    const vec2d_c s_Fs = ( s_d * ( ( L - s_L ) * ks ) );  // Hooke's Law.
    const vec2d_c s_Fd = ( s_d * ( ( v * s_d ) * kd ) );  // Damping.

    // Forces sum => Acceleration.

    const vec2d_c a = ( gravity + ( ( s_Fs - s_Fd ) / m ) );

    // Semi-implicit Euler.

    v += ( a * dt );
    p += ( v * dt );
  }

  // Particle state.

    f64_t m;
  vec2d_c p;
  vec2d_c v;

  // Spring.

    f64_t L;   // Rest length.
    f64_t ks;  // Stiffness.
    f64_t kd;  // Damping.
};

Here’s a goopy spring (\(k_s=1e3\)):

Here’s a stiff spring (\(k_s=1e6\)):

Numerical error accumulation aside, these systems neither gain nor lose energy. So pendular motion would never stop should the simulation run forever.

Note that, apparently, the stiff spring successfully satisfies the distance constraint: the particle draws a seemingly perfect semicircle, which is exactly what we want. But later we’ll see how easy-to-break this apparent stability is as soon as we complicate the system a little bit.

Damping

Damping is an influence upon an oscillatory system that has the effect of reducing or preventing the oscillation. A damped spring is often modeled like this:

Which means that the force at either end of the spring:

• Resists changes in length, proportionally to a certain stiffness constant \(k_s\).
• Resists velocity along the spring’s direction \(\mathbf{d}_t\), proportionally to a certain damping constant \(k_d\).

Here’s a damped goopy spring (\(k_d=1e3\)):

Damping a spring effectively makes it calmer by reducing its tendency to gain velocity. Which in turn makes it lose energy and stop eventually. Like in real life.

Adding more springs

As explained in the beginning of this post, the oscillatory chaos induced by the use of springs, even when they are stiff and/or damped, becomes quickly uncontrollable as more particles/springs are involved.

Let’s modify our code to connect a second spring-mass to the end of the first one.

The state struct now manages two particle states. When summing forces, note that the particle where both mass-less springs meet is subject to opposing forces coming from both springs. The particle at the end of the contraption is only subject to the force from the second spring.

Remember that our goal is a perfectly semicircular trajectory for the first particle, in yellow.

Things go crazy with a goopy double-spring now:

Damping muffles the craziness a bit, but not quite:

Using stiffer springs seems to work well at first sight… until vibrations crawl up and jerky motion takes over eventually:

Like before, damping helps. But the problem is visibly not fully gone. Also, adding more springs would make the situation spiral out of control anyway:

struct double_spring_mass_t
{
  explicit double_spring_mass_t
  
      ( const f64_t rest_length =  .5,
        const f64_t stiffness   = 1e3,
        const f64_t damping     = 1e3 )

      : m( 100 ), p1( ( rest_length * 1 ), 0 ), v1( 0 ),
                  p2( ( rest_length * 2 ), 0 ), v2( 0 ),

        L( rest_length ), ks( stiffness ), kd( damping ) {}

  void step( const f64_t dt, const vec2d_c& gravity )
  {
    // Current spring states.

    const vec2d_c s1_D  = ( p1      );      // Spring dir.
    const vec2d_c s2_D  = ( p2 - p1 );      //

    const   f64_t s1_L  =   s1_D.length();  // Spring len.
    const   f64_t s2_L  =   s2_D.length();  //

    const vec2d_c s1_d  = ( s1_D / s1_L );  // Spring dir (normalized).
    const vec2d_c s2_d  = ( s2_D / s2_L );  //

    const vec2d_c s1_Fs = ( s1_d * ( ( L  - s1_L ) * ks ) );  // Hooke's Law.
    const vec2d_c s2_Fs = ( s2_d * ( ( L  - s2_L ) * ks ) );  //

    const vec2d_c s1_Fd = ( s1_d * ( ( v1 * s1_d ) * kd ) );  // Damping.
    const vec2d_c s2_Fd = ( s2_d * ( ( v2 * s2_d ) * kd ) );  //

    const vec2d_c s1_F  = ( s1_Fs - s1_Fd );
    const vec2d_c s2_F  = ( s2_Fs - s2_Fd );

    // Forces sum => Acceleration.

    const vec2d_c a1 = ( gravity + ( ( s1_F - s2_F ) / m ) );
    const vec2d_c a2 = ( gravity + ( ( s2_F        ) / m ) );

    // Semi-implicit Euler.

    v1 += ( a1 * dt ); p1 += ( v1 * dt );
    v2 += ( a2 * dt ); p2 += ( v2 * dt );
  }

  // Particle state.

    f64_t m;   // m=m1=m2.

  vec2d_c p1, v1;
  vec2d_c p2, v2;

  // Both springs.

    f64_t L;   // Rest length.
    f64_t ks;  // Stiffness.
    f64_t kd;  // Damping.
};

Representing constraints properly

We would wish for a general method that robustly satisfies rigid constraints without large stiffness constants that induce strong oscillatory forces. One that using only corrective forces (or impulses) gently nudges particles back to valid states constraints-wise.

This is what the next entry in this mini-series will be about.

Stay tuned!