Rope

Now we’ll try to simulate a rope. Here you’ll see the advantage of PyMbs in that it is fully scriptable using Python. If you want to try it first, or look at the complete source code, see Rope.py.

We start with importing time (you’ll see why in a moment) and PyMbs, and set up our initial frame of reference as before:

import time

from pymbs.input import MbsSystem, diag, pi

world = MbsSystem([0, 0, -1])

To make our script a bit more flexible, we add the length of the rope, the rotational damping, the mass density, the radius and the number of segments used to simulate it as variables. Additionally, we choose which method is to be used to derive the equations of motion:

L = 3.0         # total length
d = 0.2         # rotational damping
rho = 500.0     # 500 kg/m^3
r = 0.1         # radius
segments = 15

# Choose between "Explicit", "OrderN" and "Recursive"
method = "Recursive"

Now we initialize the lists that will hold our bodies and joints. Also, we calculate the length, mass and inertia of each rope segment:

bodies = [None] * segments
joints = [None] * segments

# Parameters
l = L / segments
m = rho * pi * r**2 * l
I = 1 / 12 * m * l**2 * diag([1, 1, 1])

The interesting part is creating all the bodies and joints. We loop over the number of segments, create a body and an additional frame at the end for each segment. We then connect the first segment to the world and the other segments to each other. A simple line is used for visualisation. In the end, a joint load is added to each joint, with its value being calculated from multiplying the angular velocity (s[1]) with the rotational damping d:

for i in range(0, segments):
    # Create Body and Frame
    bodies[i] = world.addBody(m, cg=[l / 2, 0, 0], inertia=I)
    bodies[i].addFrame(name='end', p=[l, 0, 0])

    # Create joints
    if i == 0:
        joints[i] = world.addJoint(world, bodies[i], 'Ry')
    else:
        joints[i] = world.addJoint(bodies[i - 1].end, bodies[i], 'Ry')

    world.addVisualisation.Line(bodies[i], l)

    # Add damping
    s = world.addSensor.Joint(symbol=f'state{i}_', joint=joints[i])
    D = world.addExpression(symbol_str=f'T{i}_', exp=-d * s[1])
    world.addLoad.Joint(joint=joints[i], symbol=D)

Method of Equation Generation

Now we get to the issue of which method to choose for generating the equations of motion. The choice is between Explicit, Recursive and OrderN. The explicit generation scheme is comparatively slow and should not be used for larger systems. On my system, using eight segments as a benchmark, the Explicit scheme took 12.9 seconds to generate the equations, OrderN took 1.06 seconds and the Recursive one 0.05 seconds.

You can test these for yourself using the time module:

t = time.time()
if method == "OrderN":
    world.genEquations.OrderN()
elif method == "Recursive":
    world.genEquations.Recursive()
elif method == "Explicit":
    world.genEquations.Explicit()
print(f'Time needed for generating equations: {(time.time() - t):.2f} s')

Not only is the time to generate the equations different, the complexity of the equations differs as well. Again for eight segments, the python code generated that is used for the simulation has 26.6 million characters (!) when using Explicit, but only 172 thousand when using OrderN and 34 thousand for Recursive. So in this case you should use the Recursive scheme; for other models OrderN might be beneficial.

Now we can take a look at our rope swinging:

world.show('Rope')

This is what it should look like:

If you hit simulate with more than say six segments, the simulation will be quite slow, because the time integration has to evaluate a python function to obtain the accelerations for each time step. However, you can use the Compile Model in Fortran or Compile Model in C buttons to automatically export the equations to Fortran or C, compile them and reimport them for faster integration. This way, even 20 segments appear smooth on my machine.

In the end, the complete source looks like this:

"""
Procedurally generate a planar rope. The rope segments are ideally stiff
and connected by 1D rotational joints.
"""

# Warning: The source code of the examples is quoted in the documentation. If
# you change this file, you'll have to change the corresponding file in the
# documentation (see doc/examples).

import time

from pymbs.input import MbsSystem, diag, pi

world = MbsSystem([0, 0, -1])

L = 3.0         # total length
d = 0.2         # rotational damping
rho = 500.0     # 500 kg/m^3
r = 0.1         # radius
segments = 15

# Choose between "Explicit", "OrderN" and "Recursive"
method = "Recursive"

# Initialize lists
bodies = [None] * segments
joints = [None] * segments

# Parameters
l = L / segments
m = rho * pi * r**2 * l
I = 1 / 12 * m * l**2 * diag([1, 1, 1])

# Create bodies and connect them
for i in range(0, segments):
    # Create Body and Frame
    bodies[i] = world.addBody(m, cg=[l / 2, 0, 0], inertia=I)
    bodies[i].addFrame(name='end', p=[l, 0, 0])

    # Create joints
    if i == 0:
        joints[i] = world.addJoint(world, bodies[i], 'Ry')
    else:
        joints[i] = world.addJoint(bodies[i - 1].end, bodies[i], 'Ry')

    world.addVisualisation.Line(bodies[i], l)

    # Add damping
    s = world.addSensor.Joint(symbol=f'state{i}_', joint=joints[i])
    D = world.addExpression(symbol_str=f'T{i}_', exp=-d * s[1])
    world.addLoad.Joint(joint=joints[i], symbol=D)

t = time.time()
if method == "OrderN":
    world.genEquations.OrderN()
elif method == "Recursive":
    world.genEquations.Recursive()
elif method == "Explicit":
    world.genEquations.Explicit()
print(f'Time needed for generating equations: {(time.time() - t):.2f} s')

world.show('Rope')