The physics of Mechanisms in the space

TL;DR

Some algebra, calculus, physics put together to create.

Intro

The kind of thing that you put together in few afternoons nowadays.

You thought that your 9-5 is killing you?

Is your 5pm-9pm giving you some life?

Have you ever thought about whats Ikigai for you?

Ops, wrong post to talk about this.

Let me go back: 3D mechanisms! YES!

Here it goes all that unseen talent that the world was eager to consume.

About MBSD

Code is law, specially for multibody system dynamics.

#git clone https://github.com/JAlcocerT/mbsd

NEW MBSD

Source Code of a Python fwk to simulate mechanisms

The next sections will be theoretical, but they will be a base to simulate cool thingies.

2D MBSD is cool but

With just 2D kinematics in place you can create very nice 3D renders.

For example, this geneva drive:

cd ./3Design/mbsd-to-render/geneva-drive
#make help
make all
rsync -avP jalcocert@192.168.1.2:/home/jalcocert/3Design/mbsd-to-render/geneva-drive/render/geneva_drive.mp4 .
mpv geneva_drive.mp4

Or this Scotch…the pantograph…you name it:

cd ./3Design/mbsd-to-render/scotch-yoke
make all
rsync -avP jalcocert@192.168.1.2:/home/jalcocert/3Design/mbsd-to-render/scotch-yoke/render/scotch_yoke.mp4 .
mpv scotch_yoke.mp4
cd ./3Design/mbsd-to-render/pantograph
make all
rsync -avP jalcocert@192.168.1.2:/home/jalcocert/3Design/mbsd-to-render/pantograph/render/pantograph.mp4 .
mpv scotch_yoke.mp4

cd /home/jalcocert/Desktop/3Design/mbsd-to-render
printf "file '%s'\n" *.mp4 > concat.txt
ffmpeg -f concat -safe 0 -i concat.txt -c copy all_mechanisms.mp4
mpv all_mechanisms.mp4

3D Kinematics

But if you want to do 3D mechanics, you need to get 3D kinematics right first.

There are many interesting effects in 3D that simply dont exist in 2D.

3D Dynamics

CAD x Blender

I was using a simple slider-crank to test that its possible to go from a ~mbsd to STL’s and blender renders.

Call it blender as a code or how to bring mechanisms to life.

git clone https://github.com/JAlcocerT/3Design
cd ./3Design/z-cadquery
make check #make help
#make scene-ui #this starts blender UI
#make all
tmux new-session -d -s cad "make all" #if you will be leaving this for the night
#tmux attach-session -t cad #to see hows going
#rsync -avP jalcocert@192.168.1.2:/home/jalcocert/3Design/z-cadquery/render/slider_crank.mp4 .
#mpv slider_crank.mp4

Even though…bringing them to life would be to manufacture them, right?

#git clone git clone https://github.com/JAlcocerT/mbsd
cd ./mbsd/z-cad-render

Well…give me some time for that :)

Real Time Bike Simulator in Python

A 9-DOF multibody dynamics simulator for a bicycle, ported from MATLAB to Python.

cd ./mbsd/bike-real-time-simulator/Python_version
# Install tkinter (system dependency for matplotlib GUI)
sudo apt install python3-tk
# Install Python dependencies
#pip install -r requirements.txt

Yea…this is also 3D :)

Just that it does not uses the reference coordinate system

Why?

Because this has to be solved in real time

For that, i needed to think and optimize the model:

#uv init --no-readme .

# Add packages from requirements.txt one-shot
#uv add numpy matplotlib pynput
uv pip install -r requirements.txt
#uv sync

# Basic run — bicycle rolls forward at ~10 m/s with 3D animation
#python3 main.py

# With keyboard steering
#python3 main.py --keyboard
uv run main.py --keyboard

# Without animation (faster, plots results at the end)
#python3 main.py --no-animate

Keyboard Controls (with --keyboard)

Project Structure

Execution flow:

main.py  →  run_simulation()
              │
              ├── MBodySystem()          [parameters.py]
              ├── runge_kutta_2()        [main.py]
              │     └── equations_of_motion()  at each step
              │           ├── mass_matrix()    [forces.py]
              │           ├── gravity_forces() [forces.py]
              │           ├── coriolis_forces()[forces.py]
              │           ├── aero_forces()    [forces.py]
              │           ├── contact_forces() [forces.py]
              │           │     └── A2, A5, H2, H5, H2p, H5p  [kinematics.py]
              │           └── G2g, G4g         [kinematics.py]
              └── BicycleAnimator.update()  [visualization.py] (every 2 steps)

The theoretical basis for this is at: https://github.com/JAlcocerT/Bike_dynamic_simulator

Dont forget about the steer to the fall concept!

It computes how the vehicle moves over time under the influence of gravity, aerodynamic drag, tire-ground contact forces, and externally applied steering/pedaling torques.

The Python code is a direct port of the original MATLAB simulator (Matlab_version/).

Conclusions

MBSD

Source Code of a MBSD framework in Python

#git clone https://github.com/JAlcocerT/VideoEditingRemotion
#cd remotion-cc

#git clone https://github.com/JAlcocerT/3Design
cd ./3Design/mbsd-to-render

You see how I write about many different stuff?

find ./content/blog -maxdepth 1 -type f -name "*.md" | wc -l #this post is 350+ and 80+ expected for this year!

find content/blog -name '*.md' -print0 |
xargs -0 awk '
  FNR==1 { 
    post_date=""; 
    printed=0 
  }

  /^date:/ && !printed {
    gsub(/^date:[[:space:]]*/, "", $0)
    post_date = substr($0, 1, 10)
    if (post_date >= "2019-01-01" && post_date <= "2026-03-31") {
      print FILENAME ": " post_date
      printed=1
      count++
    }
  }

  END {
    print "TOTAL:", count+0
  }'

If any of the things I wrote about catches your attention and you want my time:

Consulting Services

Consulting - Tier of Service

DIY via ebooks

Distilled knowledge via web/ooks to enable you to create

Upcoming topics with 3D mechanics:

1. Engine configuration analysis
2. Suspensions: Double Wishbone, macpherson…

https://youtu.be/y-ANdaUthxg

MBSD x Web Dev

Because if mechanism 3D dynamics its kind of trivial now.

So is web development.

#git clone https://github.com/JAlcocerT/Slider-Crank
cd ./Slider-Crank/landing #https://multibodysystemdynamics.pages.dev/

FAQ

The launch strategy: aka, focus strategy

Non for comercial purposes :)

The Tier of Service: DIY (1b - leverages on actual tech stack Ive put together - PaaS x (WP/Ghost or SSG+CMS))

The Tech Stack:

Interactivity

• https://brm.io/matter-js/

Matter.js is a 2D physics engine for the web

• ThreeJS and D3JS

Bike MultiBody Model

The bicycle is modelled as 4 rigid bodies connected by joints: just 9 degrees of freedom remaining!

The simulator uses multibody dynamics in generalized coordinates — a classical analytical mechanics approach where the system state is described by the minimum set of independent coordinates that fully describe its configuration.

What the model does NOT capture

• Road camber, bumps, or surface roughness
• Tire sidewall dynamics or pneumatic pressure effects
• Structural flexibility (frame flex, fork flex)
• Suspension (no spring/damper between wheel and frame)
• Rider body dynamics (rider is rigidly attached to frame, body 3)
• Wind, rain, or temperature effects
• Motor/drivetrain dynamics (pedaling torque is applied directly as a generalized force)

Each rigid body has its own local frame, whose orientation relative to the global frame is given by a rotation matrix Aᵢ(q). Body frames are right-handed and aligned with the principal axes of inertia when at zero lean/pitch/steer.

Angular rotation convention: Euler angles are applied in this sequence for the frame body (body 3).

1. Roll φ (phi) — rotation about the X axis (lean)
2. Yaw ψ (psi) — rotation about the Z axis (heading change)
3. Pitch θ (theta) — rotation about the Y axis (tilt forward/back)

The resulting rotation matrix A3 maps body-frame vectors to global-frame vectors:

r_global = A3(q) · r_body

Each body has:

• A mass mᵢ
• A center of gravity (CoG) position expressed in body-frame coordinates
• A 3×3 inertia tensor Iᵢ in body frame

The system has 9 generalized coordinates q ∈ ℝ⁹ (0-indexed in Python):

The frame reference point (Rx3, Ry3, Rz3) is a point fixed to body 3 near the bottom bracket.

It is not the CoG of body 3 — the CoG is offset by [XG3, 0, ZG3] = [0.20, 0, 0.10] in the body frame.

The full state vector is y = [q; q̇] ∈ ℝ¹⁸, where q̇ are the generalized velocities.

Why these coordinates?

• The 3 translational DOFs (Rx3, Ry3, Rz3) locate the bicycle in space.
• Roll φ, yaw ψ, pitch θ orient the frame. This ordering produces a rotation matrix without singularities at the operating angles of interest.
• Steering γ is relative to the frame — it is the handlebar rotation about the (tilted) head tube axis.
• Wheel spin angles ε and ν are needed to compute tire slip and gyroscopic forces.

Equations of Motion

The dynamics follow the Newton-Euler formulation in generalized coordinates:

M(q) · q̈ = Q(q, q̇)

Where:

• q — 9-element generalized coordinate vector
• M(q) — 9×9 symmetric positive-definite mass matrix (configuration-dependent)
• Q(q, q̇) — 9-element generalized force vector (sum of all forces)

Q is assembled from four independent contributions:

Q = Q_grav + Q_coriolis + Q_aero + Q_contact + Q_steer + Q_pedal

Mass Matrix

Computed using the composite rigid body algorithm:

M = Σᵢ [ mᵢ · Jvᵢᵀ Jvᵢ  +  Jwᵢᵀ · Iᵢ_global · Jwᵢ ]

Where Jvᵢ and Jwᵢ are the translational and angular Jacobians of body i, and Iᵢ_global = Aᵢ · Iᵢ_body · Aᵢᵀ transforms the body-frame inertia tensor to global frame.

Jacobians are computed numerically via central finite differences on the rotation matrices and CoG position functions (step size ε = 1e-7). This avoids translating 590 lines of symbolic MATLAB code at the cost of ~18 extra mass_matrix() calls per Coriolis evaluation.

Coriolis & Centrifugal Forces

Computed from Christoffel symbols of the first kind:

Qvᵢ = -Σⱼₖ Γᵢⱼₖ · q̇ⱼ · q̇ₖ

Γᵢⱼₖ = ½ (∂Mᵢⱼ/∂qₖ + ∂Mᵢₖ/∂qⱼ − ∂Mⱼₖ/∂qᵢ)

∂M/∂qₖ is evaluated numerically (central finite differences on mass_matrix()).

Tire-Ground Contact

Each wheel uses an elastic contact model:

Normal force — spring-damper with indentation δ:

Fn = Kn · |δ|  +  c_damp · δ̇ · δ       (active when δ < 0)

Tangential force — piecewise linear (bristle model) saturated by a friction ellipse:

Longitudinal:  Fx = μx · Fn · (κ / κc)      if |κ| ≤ κc
               Fx = μx · Fn · sign(κ)         otherwise

Lateral:       Fy = μy · Fn · (α / αc)       if |α| ≤ αc
               Fy = μy · Fn · sign(α)         otherwise

With elliptic saturation: if (Fx/μx)² + (Fy/μy)² ≥ Fn², both components are scaled back onto the friction ellipse.

• κ = longitudinal slip ratio (wheel slip vs. rolling speed)
• α = slip angle (angle between wheel heading and velocity vector)
• κc, αc = critical slip thresholds (both default to 0.12 rad)
• μx = μy = 0.9 — friction coefficients

The contact point geometry accounts for the toroidal tire cross-section (ellipse semi-axes a_n = b_n = 0.02 m), which shifts the effective contact radius with lean angle.

Applied Torques

Steering and pedaling torques are projected into generalized coordinates using the angular Jacobians G4g (handlebar) and G2g (rear wheel):

Q_steer = G4gᵀ · [0, 0, τ_steer]
Q_pedal = G2gᵀ · [0, τ_pedal, 0]

State Vector & Generalized Coordinates

The state vector y has 18 elements: y = [q; q̇]

Generalized Coordinates q (indices 0–8)

Generalized Velocities q̇ (indices 9–17)

Same ordering as coordinates. Initial conditions:

q0 = [0, 0, 0.97, 0, 0, 0, 0, 0, 0]   # starts at 0.97 m height
v0 = [10, 0, 0, 0, 0, 0, 0, 10/0.34, 10/0.34]  # 10 m/s forward, wheels spinning