modular spring damper

Modular Mass-Spring-Damper

By building these components in a modular manner, we can assemble any number of masses, springs, and dampers and solve the system without the need of matrix methods and state-space modeling.

This approach is extremely powerful when designing large, modular simulations.

This code is also in the examples folder of the Ascent repository: modular-spring-damper

Important
When solving modular systems, odeint integration algorithms cannot be used without copying the state vector within the simulation loop. Do not simply swap an Ascent integrator with an odeint integrator, as the odeint integrator will not get the correct solution. Using odeint integration schemes for modular systems is discussed here: Modular Systems with odeint.

The Body class

The asc::state_t is by default a std::vector<double>, it can easily be redefined to std::vector<float> or to another container such as std::deque<double>. However, changing the container from a vector will affect design choices and speed.

The state input to the Body class is the entire state of all parameters to be integrated.

An asc::Param is a single parameter within the state. It contains a reference to an element within the std::vector, which means that the state vector must have enough memory reserved at the time of initializing all asc::Params. If the vector reallocates, then all the asc::Params will become invalid. In most cases it is recommended to use the reserve method to reserve a large amount of memory and therefore avoid reallocation. Use reserve and not resize because resize will allocate memory.

struct Body
{
   Body(asc::state_t& state) : s(state), v(state) {}

   asc::Param s; // position
   asc::Param v; // velocity
   double m{}; // mass
   double f{}; // force

   void operator()(const asc::state_t&, asc::state_t& D, const double)
   {
      s(D) = v;

      if (m > 0.0)
         v(D) = f / m;
      else
         v(D) = 0.0;

      f = 0.0;
   }
};

D is the derivative state vector.

The code s(D) = v; is assigning velocity to the derivative of position.

The code v(D) = f / m is assigning force / mass to the derivative of velocity (i.e. acceleration).

The Damper class

The Damper class takes references to two bodies. Note that the damper becomes invalid if the bodies are destroyed.

Even though the Damper doesn't require state, derivative, or time input to its update (overloaded access operator) method, we still use the same input syntax. This is not required, but it helps when writing generic callers. This simplifies writing interface code Ascent's scripting (via ChaiScript).

struct Damper
{
   Damper(Body& b0, Body& b1) : b0(b0), b1(b1) {}

   Body& b0;
   Body& b1;

   double dv{}; // velocity difference
   double c{}; // damping coefficient
   double f{}; // force

   void operator()(const asc::state_t&, asc::state_t&, const double)
   {
      dv = b0.v - b1.v;
      f = c*dv;

      b0.f -= f;
      b1.f += f;
   }
};

The Spring class

struct Spring
{
   Spring(Body& b0, Body& b1) : b0(b0), b1(b1)
   {
      l0 = b1.s - b0.s;
   }

   Body& b0;
   Body& b1;

   double l0{}; // initial spring length (distance between masses)
   double ds{}; // spring compression/extension
   double k{}; // spring coefficient
   double f{}; // force

   void operator()(const asc::state_t&, asc::state_t&, const double)
   {
      ds = l0 + b0.s - b1.s;
      f = k*ds;

      b0.f -= f;
      b1.f += f;
   }
};

Setup and Run

using namespace asc;

int main()
{
   state_t x;
   x.reserve(100); // We reserve more space than necessary, but Ascent will only allocate what is needed
   double t = 0.0;
   double dt = 0.01;
   double t_end = 1.5;

   Body b0(x);
   Body b1(x);
   b1.m = 1.0;
   b1.s = 1.0;
   b1.v = 40.0;

   Spring spring(b0, b1);
   spring.k = 2000.0;

   Damper damper(b0, b1);
   damper.c = 5.0;

   RK4 integrator;
   Recorder recorder;

   auto system = [&](const asc::state_t& x, asc::state_t& D, const double t)
   {
      // We must run the spring and damper before the body in order to accumulate forces
      spring(x, D, t);
      damper(x, D, t);
      b1(x, D, t);
   };

   while (t < t_end)
   {
      recorder({ t, b1.s });
      integrator(system, x, t, dt);
   }

   recorder.csv("spring-damper", { "t", "b1 position" });

   return 0;
}