Q4: 3D Version: Experiment with different mechanical structures

Q4: 3D Version: Experiment with different mechanical structures#

Model 1 : Extend the double-pendulum to a chain#

Complete code – chian with springs is in ‘chain_spring.ipynb’

Complete code – chain with constraints is in ‘chain_constraint.ipynb’

Chain with spring#

Assume the quantity of each mass: 1, length of the spring:1, the spring constant:1000, still fixed point at the origin.

We only need to change the geometry part. Here we write a for loop to add more masses and springs:

#create for loop
for i in range(N_MASSES):
    
    initial_pos = (SPRING_LENGTH * (i + 1), 0.0, 0.0)
    
    # add new mass
    current_mass_connector = mss.add(Mass(MASS_VAL, initial_pos)) 
    
    # add new spring
    mss.add(Spring(SPRING_LENGTH, 
                   SPRING_STIFFNESS, 
                   (previous_connector, current_mass_connector)))
    
    # refresh
    previous_connector = current_mass_connector

Such that everytime, we only need to change N_MASSES

Test 1: When mass = 5

Test 2: When mass = 10

Chain with constaint#

Similarly, we write a for loop to add mass and constraint.

for i in range(N_MASSES):
    initial_pos = (LINK_LENGTH * (i + 1), 0.0, 0.0)
    
    # add new mass
    current_mass_connector = mss.add(Mass(1, initial_pos)) 
    
    # add new constraint and connet it to the previous mass
    dc = DistanceConstraint(previous_connector, current_mass_connector, LINK_LENGTH)
    mss.addDistanceConstraint(dc)
    
    # refresh
    previous_connector = current_mass_connector

Such that everytime, we only need to change N_MASSES

N_MASSES = 5
LINK_LENGTH = 1.0

Test 1: When mass = 5

Test 2: When mass = 10

Reuslt:

  • With larger masses, the system begins to approximate a continuous body (such as a soft rope or a whip).

  • The end mass resembling the tip of a whip, characterized by the propagation of waves and accelerated whipping motion.

  • The hole system will turn slow, bulk oscillation, and faster, small-scale localized twisting and vibrational patterns emerging throughout the chain.

Model 2 : Crane Structure and Vibration#

Complete code – ‘crane.ipynb’

In this section, we build a complex mechanical structure (a crane) by using the Mass-Spring System combined with Lagrangian Constraints.

  1. Model Construction Strategy

    We implemented a build_crane function in Python. The structure consists of two main types of connections:

    • Rigid Beams: Modeled using DistanceConstraint. We want the frame of the crane to be rigid and not oscillate like jelly.

    • Cables (Elastic): Modeled using Spring. We want to observe vibration when lifting a mass.

    We defined a helper function to easily switch between these two types:

def add_beam(mss, node1, node2, length, stiffness, use_constraint=True):
    if use_constraint:
        # Use Lagrangian constraint for rigid parts (Steel)
        dc = DistanceConstraint(node1, node2, length)
        mss.addDistanceConstraint(dc)
    else:
        # Use Spring for elastic parts (Cables/Rubber)
        mss.add(Spring(length, stiffness, (node1, node2)))

  1. Building the Geometry

The crane consists of a vertical tower (10 floors) and a horizontal arm (6 units). To ensure structural stability we added crossed diagonal supports to every face.

① The tower loop:

for i in range(floors):
    # ... (node creation logic) ...
    # connecting floors
    add_beam(mss, c_curr, c_next, L, 0, use_constraint=True) #vertical union
    add_beam(mss, c_next, c_next_neighbor, L, 0, use_constraint=True)  #horizontal union
    
    diag_len = math.sqrt(2) * L 
    add_beam(mss, c_curr, c_next_neighbor, diag_len, stiffness_beam, use_constraint=False)  #diagonal (X shape) for stability

The stiffness for the structural parts is 0 because their rigidity will be actually being defined by the Lagrange Constraints, while for the elastic springs we fixed a value stiffness_beam= 5000. We coded the arm loop analougsly.

② The Load and Cables:

We attach a heavy mass (\(m=5.0\)) to the tip of the arm using elastic springs (stiffness_cable = 3000.0). This way we introduce the vibration.

# Cables attached to the tip of the arm
mss.add(Spring(dist_h, stiffness_cable, (tip_node_1, load_mass)))
mss.add(Spring(dist_h, stiffness_cable, (tip_node_2, load_mass)))

  1. Simulation and Visualization

To visualize the movement smoothly, we modified the simulation loop in the Jupyter Notebook. Since the code with the Lagrange Constraints is computationally heavy, we decouple the physics steps from the rendering steps using steps_per_frame so the computer do not struggle that much to calculate the physics and draw the crane quickly at the same time.

steps_per_frame = 2  #render 1 frame every 2 physics calculations

for i in range(1000):
    mss.simulate(0.01, 2) 

    if i % steps_per_frame == 0:
        # we just actualise the graphic positions when we plot (1 over 2)
        sleep(0.01)

  1. Result

The simulation shows the crane structure remaining rigid due to the DistanceConstraints (orange lines), while the load and arm bounce realistically due to the Springs (cyan lines).

Model 3 : Replace the Springs of the Double-Pendulum by Distance Constraints#

Complete code – ‘dist_constraints.ipynb’

In the original 3D mass–spring example, the double pendulum is modeled using two springs:

  • one connecting the fixed point to the first mass,

  • one connecting the first mass to the second mass.

In this task, we remove the springs entirely and replace them with holonomic distance constraints enforcing:

\[ \|x_i - x_j\| = L_0. \]

This converts the system from a soft spring model into a rigid “rod-like” double pendulum where the distances remain constant throughout the simulation.

  1. Remove the springs

The original example contained:

mss.add(Spring(1, 200000, (f1, mA)))
mss.add(Spring(1, 100000, (mA, mB)))

These are completely removed.

  1. To avoid external forces at initialization, the rest-lengths are computed directly from the initial distances:

def distance(p, q):
    return ((p[0]-q[0])**2 +
            (p[1]-q[1])**2 +
            (p[2]-q[2])**2 )**0.5

L1 = distance(mss.fixes[f1.nr].pos, mss.masses[mA.nr].pos)
L2 = distance(mss.masses[mA.nr].pos, mss.masses[mB.nr].pos)

  1. Add distance constraints

Using DistanceConstraint springs are removed and replaced by fixed-length distance constraints.

dc1 = DistanceConstraint(f1, mA, L1)
mss.addDistanceConstraint(dc1)

dc2 = DistanceConstraint(mA, mB, L2)
mss.addDistanceConstraint(dc2)

  1. Constraint segments are drawn the same way as springs, but stay at a fixed length:

springpos = []
for c in mss.constraints:
    pA = mss.fixes[c.c1.nr].pos if c.c1.type == 1 else mss.masses[c.c1.nr].pos
    pB = mss.fixes[c.c2.nr].pos if c.c2.type == 1 else mss.masses[c.c2.nr].pos
    springpos.append([pA, pB])

springs.geometry = LineSegmentsGeometry(positions=springpos)

  1. When the Lagrange multiplier becomes very large, the Jacobian matrix used in the Newton solver becomes numerically ill-conditioned.

Error Newton did not converge exist.

We increase the tolerance and max. iteration to 1e-9 and 20. And slow down the video sleep(0.3)

Finally we get:

Obviously, the constraint connection lines maintain a fixed length, regardless of how the mass blocks move. The segments do not stretch or contract noticeably like a soft spring, making the connections appear as metal rods.

Model 4 : Spinning top#

Complete code – ‘bind_speed.cpp’ and ‘spinning.ipynb.ipynb’

  1. First, we need to add initial velocity for the masses. So we can modify bind_speed.cpp

Attention: ‘vel’ has already been defined in hpp file, so we can define the property directly.

//add initial velocity
      .def_property("vel",
                    // getter(return array<double, 3> as output)
                    [](Mass<3> & m) { return py::cast(std::array<double, 3>{m.vel(0), m.vel(1), m.vel(2)}); },
                    // setter (accept array<double, 3> as input)
                    [](Mass<3> & m, std::array<double, 3> v) { m.vel = Vec<3>{v[0], v[1], v[2]}; });

  1. Then we define the geometry in spinning.ipynb

# geometry model
R = 0.5           # radius of the base triangle (m)
H = 1.0           # height (m)
OMEGA = 20.0      # angular velocity (rad/s)
L_BASE = R * np.sqrt(3) # length of the base triangle sides

# fixed point at origin
f1 = mss.add(Fix((0,0,0))) 

# 3 masses as the base triangle vertices
m_mass = 1.0
phi1 = 0
phi2 = 2 * np.pi / 3
phi3 = 4 * np.pi / 3

# initial positions of the 3 masses
mA = mss.add(Mass(m_mass, (R * np.cos(phi1), R * np.sin(phi1), H)))
mB = mss.add(Mass(m_mass, (R * np.cos(phi2), R * np.sin(phi2), H)))
mC = mss.add(Mass(m_mass, (R * np.cos(phi3), R * np.sin(phi3), H)))


# set initial speed
for m in mss.masses:
    # get current position 
    x, y, z = m.pos[0], m.pos[1], m.pos[2] 
    
    # seperate speed
    vx = -OMEGA * y
    vy = OMEGA * x
    vz = 0.0
    
    # input vel
    m.vel = (vx, vy, vz)


# add distance constraints
# fix-mass
L_TIP = np.sqrt(R**2 + H**2)
mss.addDistanceConstraint(DistanceConstraint(f1, mA, L_TIP))
mss.addDistanceConstraint(DistanceConstraint(f1, mB, L_TIP))
mss.addDistanceConstraint(DistanceConstraint(f1, mC, L_TIP))

# inside mass
mss.addDistanceConstraint(DistanceConstraint(mA, mB, L_BASE))
mss.addDistanceConstraint(DistanceConstraint(mB, mC, L_BASE))
mss.addDistanceConstraint(DistanceConstraint(mC, mA, L_BASE))

  1. Finally, we can get the result:

Obviously, at first, the Kreisel rotate fast and with decreasing of energy, the spinning top falls down.

Contents