OpenMBD: An Open-Source Multibody Dynamics Simulator for Biomechanics Research and Education

Implementation

System architecture

OpenMBD is implemented in Python 3 and structured as a modular collection of sixteen source files, each with a clearly defined responsibility. The package has three external dependencies: NumPy¹¹ for all numerical linear algebra; Matplotlib¹² for three-dimensional visualisation and analysis plotting; and Pillow,¹³ used exclusively for GIF animation export. The graphical user interface is built on Tkinter, which is distributed as part of the Python standard library. All remaining imports are standard library modules. A browser-based model editor (openmbd_model_editor.html) is provided as a standalone HTML5/Three.js application requiring no installation.

The codebase is divided into four functional layers: (i) physics - the dynamics engine (physics_engine.py), rigid body representation (physics_rigid_body.py), contact mechanics (physics_constraints.py), joint limit enforcement (physics_joint_limits.py), and utility mathematics (physics_utils.py); (ii) models - the JSON parser (models_parser.py), kinematic multibody tree (models_multibody.py), and simulation configuration (models_config.py); (iii) geometry - ellipsoid geometry (geometry_ellipsoid.py) and mesh generation (geometry_mesh_generator.py); and (iv) user interface - the main window, setup, run, analysis, and controls modules ( Figure 1).

Figure 1. OpenMBD model editor.

Equations of motion

The dynamics of the multibody system are formulated using the Principle of Virtual Power (PVP). For a system of N rigid bodies connected by kinematic joints, the generalised equations of motion take the matrix form ( Equation 1).

Generalised coordinate parameterisation

Joint rotations are parameterised using ZYX Euler angles [α, β, γ] stored in radians in the state vector. The angular velocity mapping from Euler rate vector to body angular velocity is given by the analytic Jacobian in the child body’s local frame ( Equation 3).

Recursive Newton-Euler algorithm for generalised forces

The generalised force vector Q is computed using a Recursive Newton-Euler Algorithm (RNEA),¹⁴ implemented in the rnea() method of PhysicsEngine. The algorithm proceeds in two passes over a breadth-first traversal of the kinematic tree.

Forward pass: For each body, joint velocities and centripetal accelerations are propagated from root to leaves. At each joint, the appropriate kinematic mapping is applied: a scalar axis projection for revolute joints; the ${\mathrm{E}}_{\text{body}}$ ZYX Jacobian for spherical joints; and the identity-plus ${\mathrm{E}}_{\text{world}}$ mapping for the free root joint.

Backward pass: For each body in reverse order, the net force and moment required to produce the bias acceleration under gravity and external contacts is computed, propagated through the tree, and projected onto the joint axes to obtain Q. This formulation is O(N) in the number of bodies.

Rigid body representation

Each rigid body is represented by the RigidBody class (physics_rigid_body.py), which stores mass, a 3 × 3 body-frame inertia tensor, centre of mass (CoM) position in the local body frame, and an ordered list of ellipsoid geometry objects. The inertia tensor is specified as a six-component vector [Ixx, Iyy, Izz, Ixy, Iyz, Ixz] from which the symmetric 3 × 3 matrix is assembled. Positive semi-definiteness is enforced by eigen decomposition with clamping of negative eigenvalues to a minimum of 1 × 10⁻⁶ kg·m². Body orientation is maintained as both a rotation matrix R (for kinematic operations) and a unit quaternion (for output), with periodic Gram-Schmidt re-orthonormalisation to prevent drift over long simulations.

Ellipsoid geometry representation

Body geometry is represented exclusively by ellipsoids, defined by three semi-axis dimensions [a, b, c] and a 4 × 4 local transform matrix specifying the ellipsoid centre position and orientation relative to the body origin frame (geometry_ellipsoid.py). This representation provides smooth, differentiable geometry that avoids the sharp-edge artefacts common in polygonal mesh contact formulations, while remaining computationally efficient for contact detection.

Contact detection and force model

Contact detection is performed at each time step for all ellipsoid pairs satisfying the inter-model condition, and optionally for non-adjacent same-model pairs. For ground contact (z = 0 plane), penetration depth is computed analytically using the quadratic z-extent of each ellipsoid. For body-body contacts, gradient-based outward surface normal are computed at each candidate contact point. A fast bounding-sphere cull eliminates non-overlapping pairs before the more expensive computation. Contact forces ( ${\mathrm{F}}_{\text{total}})$ are computed using a penalty-based nonlinear viscoelastic model ( Equation 4).

Kinematic adjacency exclusion

To prevent spurious contact forces between anatomically adjacent body segments that overlap by design in the resting pose, OpenMBD maintains a cached set of excluded body-index pairs. Pairs within two kinematic hops are permanently excluded from contact detection. Additionally, initial penetration offsets δ₀ are recorded for all same-model ellipsoid pairs at t = 0, and contact forces are scaled by max(0, δ − δ₀), ensuring that designed overlap generates no force at rest while genuine new contacts during motion are fully resolved.

Physiological joint range-of-motion enforcement

Physiological range-of-motion (ROM) limits are enforced by a penalty spring-damper that activates only when a joint angle violates its prescribed bounds for joint DOF with angle ( $\mathrm{\theta }$ ) limits ( Equation 5)

Passive joint damping

A passive viscous damping torque is applied across all non-root joint DOFs ( Equation 6).

Prescribed joint torque pulses

OpenMBD supports the application of time-limited generalised torques to individual joints, enabling simulation of impulsive joint loading scenarios without requiring a full muscle model. Each prescribed torque is defined by a peak magnitude vector, a start time, and a pulse duration, and is applied as a smooth half-sine envelope over the specified interval.¹⁶ This facility allows researchers to impose controlled joint-level perturbations (e.g., a hip abductor impulse, a lumbar flexion torque, or a shoulder extension pulse) as initial conditions for a subsequent passive dynamics simulation. For each registered torque entry, a scalar envelope function γ(t) modulates the peak torque magnitude over the pulse window [t₀, t₀ + T] ( Equation 7).

Since the prescribed torque is added directly to B in generalised coordinates, it is automatically projected through the mass matrix inversion A⁻¹ B at each timestep, producing physically consistent accelerations that respect the inertial properties and kinematic constraints of the full multibody system. No special handling is required for joint type. The same mechanism applies to spherical (3 DOF), revolute (1 DOF), and free (7 DOF) root joints, with the torque vector sliced to the appropriate DOF count. Multiple torque pulses may be active simultaneously on different joints of the same model or across different loaded models. Each pulse is stored as an independent entry in the prescribed_torques list in PhysicsEngine, keyed by model index and joint name. Entries are populated from ModelConfig.joint_torques at finalise time via _load_prescribed_torques_from_config(), which purges stale entries for the relevant model index before registering new ones, ensuring that repeated calls to rebuild_physics() always produce a clean state.

Mass matrix assembly and numerical integration

At each time step, the symmetric positive-definite generalised mass matrix (M) is assembled from contributions of all rigid bodies using the analytic Jacobians. A regularisation term ε = 1 × 10⁻⁴ is added to the diagonal to ensure numerical non-singularity. The equation of motion is solved using numpy.linalg.solve, with fallback to the pseudoinverse if a singular matrix is encountered. Accelerations are clamped to ±10⁵ rad/s² to prevent blow-up from pathological initial conditions. Numerical integration uses the Symplectic (semi-implicit) Euler scheme ( Equation 9 and 10).

${\dot{\mathrm{q}}}_{\mathrm{n}+1}\text{and}{\mathrm{q}}_{\mathrm{n}+1}$ are the generalised velocity and position for the next time step, respectively. $\mathrm{\Delta t}$ is the integration time step. The default time step is $\mathrm{\Delta t}$ = 0.0001 s (10,000 Hz). ${\mathrm{M}}^{-1}$ is the inverse of the generalised mass matrix. $\mathrm{Q}({\mathrm{q}}_{\mathrm{n}},{\dot{\mathrm{q}}}_{\mathrm{n}},{\mathrm{t}}_{\mathrm{n}})$ is the generalised force vector at the current state.

Model file format

Rigid body systems are defined in a human-readable JSON format (models_parser.py). Each file contains two top-level objects: bodies and joints. Each body entry specifies mass (kg), a six-component inertia vector (kg·m²), CoM position (m), and a list of ellipsoids, each with semi-axis dimensions, a local 4 × 4 transform, and a piecewise-linear force-penetration curve. Each joint entry specifies its type (fixed, revolute, spherical or free), parent and child body names, and two 4 × 4 homogeneous transform matrices T1 and T2 defining the joint anchor location and orientation relative to each body. OpenMBD ships with the model files described below.

Adult male model (openmbd_male.json): The male model (rigid bodies: 21; joints: 21; ellipsoids: 35; height: 175.4 cm; body mass: 90.3 kg¹⁷) is parameterised as follows. Segment masses and principal inertia tensors are derived from the de Leva¹⁸ regression equations. Segment lengths and joint centre locations follow Winter¹⁹ and Clauser et al.²⁰ Hip joint centre positions are computed from the regression equations of Bell et al.,²¹ which predict centre location from inter-ASIS distance and pelvic dimensions. External body proportions draw on the ANSUR II²² anthropometric survey and Pheasant.²³ Piecewise-linear contact force-penetration curves are assigned per segment using published experimental data: thorax stiffness follows Kroell et al.²⁴; head impact properties follow Nahum et al.²⁵; lower-limb soft tissue dynamic properties are from Saraf et al.²⁶; gluteal adipose tissue viscoelasticity draws on Gefen and Haberman²⁷; foot plantar fascia stiffness references Ker et al.²⁸ and Cavanagh and Lafortune²⁹; cortical and trabecular bone stiffness draws on Yamada.³⁰

Adult female model (openmbd_female.json): The female model derives from the male model and uses sex-specific anthropometric data from CDC NHANES 2021–2023¹⁷ for standing height (161.2 cm) and body mass (77.9 kg), with segment mass fractions and radii of gyration from the de Leva¹⁸ female regression tables. Biacromial and biiliac width scaling uses the ANSUR data summarised in Gordon et al.²² All contact force-penetration curves use the same experimental sources as the male model, scaled to female segment geometry and mass.

Car model (openmbd_car.json): The multi-rigid-body car model (total mass: 1,450 kg; height: 158 cm; 12 bodies; 75 ellipsoids) represents a generic mid-size passenger sedan decomposed into physically distinct contact components. Each component is a separate rigid body connected to the central chassis via a fixed joint, preserving single 6-DOF vehicle kinematics. The 13 bodies are: chassis/floor assembly (10 ellipsoids), pillars (6), front bumper (11), bonnet (4), front fenders (6), windscreen (3), roof (3), doors (12), rear screen (2), boot lid (5), rear bumper (5), front wheels (6) and rear wheels (6). Principal inertia tensors were computed analytically using thin-plate approximations for panel components and rectangular-box approximations for the chassis, bumpers and door assembly. Front-end stiffness follows Maki and Kajzer³¹: bonnet surface ellipsoids with post-peak reduction reflecting panel fold with front bumper ellipsoids capture the foam crush plateau followed by beam engagement.³² Windscreen and rear screen ellipsoids reflect laminated glass behaviour.³³ Door and roof panels follow Maltese et al.³⁴ rear bumper follows Untaroiu et al.³⁵ Pillar and structural ellipsoids (chassis rails, sills, cross-members) follows Untaroiu et al.³⁵ Unloading curves are set at 20–25% of loading stiffness, consistent with measured panel coefficients of restitution of 0.10–0.15.³¹

Bicycle model (openmbd_bicycle.json): The bicycle is represented as a simple three-body system (frame, front wheel, rear wheel) with a total mass of 18.3 kg and height of 133.2 cm. Joint definitions and contact geometry are consistent with Maki and Kajzer.^31,36

Operation

System requirements

OpenMBD runs on any platform supporting Python 3.8 or later. The following software environment is required:

Installation

The source code is available from the GitHub repository (see Software Availability). Dependencies are installed with a single command:

No package installation of OpenMBD itself is required. The simulator is launched by executing the entry-point script directly:

Workflow

The standard simulation workflow proceeds through three stages, each corresponding to a tab in the graphical user interface:

Stage 1 — Setup:

Stage 2 — Run:

Stage 3 — Analysis:

Output files

Three output files are automatically generated with a timestamp-based filename on completion of each simulation run:

Model editor

A standalone browser-based model editor (openmbd_model_editor.html) is provided for inspection and modification of JSON model files without running Python. The editor uses Three.js (r128) for interactive 3D visualisation of the ellipsoid geometry and joint hierarchy. Features include a body hierarchy tree view, per-body property editing (mass, inertia, CoM), ellipsoid dimension and transform editing, joint type and anchor editing, interactive force-curve editing, model scaling by height and mass, undo history, and JSON import/export. The editor requires no installation and runs in any modern web browser ( Figure 1).

Use Case 1: Fall simulation

This use case demonstrates the reconstruction of a typical fall experienced by an older adult, a common mechanism of injury in community and residential care settings. Falls among older adults frequently occur during everyday activities such as walking, turning, or losing balance while standing, often resulting in ground impacts involving the hands, knees, hips, or head. The simulation therefore represents a scenario where a human model loses balance and falls under gravity ( Figure 2).

Figure 2. Example OpenMBD reconstruction of a fall.

The articulated multibody representation allows the simulation to resolve segment-level kinematics and contact forces during the impact phase. Sequential contact events are detected between body segments and the ground plane, and the contact model computes forces based on penetration depth, surface stiffness, damping, and friction parameters. These forces propagate through the skeletal chain, influencing joint motion and subsequent impacts. Simulation outputs include the full-body kinematic trajectories, location, and magnitude of ground contact forces. These results enable reconstruction and analysis of impact severity, body segment loading, and potential injury mechanisms, which are particularly relevant for studying fall-related injuries in older adults such as wrist fractures, hip fractures, and head impacts. Such simulations can support research into fall biomechanics, injury risk assessment, and the evaluation of fall prevention strategies, including assistive devices, environmental modifications, or protective equipment.

Use Case 2: Rugby tackle simulation

This use case represents a player-player collision scenario typical of a rugby tackle. Two human models are initialised with velocities, representing the tackler and the ball carrier, and positioned to collide during motion. During simulation, the collision generates multiple contact events between body segments, depending on the initial configuration. Outputs include full-body kinematics, contact force magnitudes, and the location of impact. These simulations can be used to investigate head acceleration events, tackle mechanics, and player safety interventions, making them relevant for sports biomechanics research and injury prevention.

Use Case 3: Vehicle-Pedestrian impact simulation

This use case demonstrates a vehicle-pedestrian collision scenario using a simplified vehicle model interacting with the articulated male model. The vehicle is represented as a rigid body with an initial forward velocity, while the pedestrian model is positioned upright in the vehicle path ( Figure 3). Upon impact, the simulation captures the sequence of contacts between the vehicle and body segments, followed by potential secondary impacts with the ground. The multibody dynamics framework enables analysis of body segment trajectories and ground impact events. Simulation outputs provide detailed kinematic data and contact force histories, supporting analysis relevant to pedestrian injury biomechanics, crash reconstruction, and road safety research.

Figure 3. Example OpenMBD reconstruction of a vehicle-pedestrian impact.

Use Case 4: Custom model definition using the JSON format and model editor

This use case demonstrates the extensibility of OpenMBD by describing the workflow for defining a novel rigid body system using the open JSON model format and the browser-based model editor. Users who wish to simulate systems not included in the distributed model library, for example, a simplified prosthetic limb, an animal model, or an industrial mechanism can define a new model without writing any Python code.

JSON model structure

A minimal OpenMBD model JSON file requires two top-level fields: bodies and joints. The inertia field accepts a six-component vector [Ixx, Iyy, Izz, Ixy, Iyz, Ixz] in kg·m². The force_curve field is a piecewise-linear lookup table mapping penetration depth (m) to elastic contact force (N); the shape of this curve determines the stiffness of the body surface and should be calibrated to the properties of the segment being modelled.

Model editor workflow

The browser-based model editor (openmbd_model_editor.html) provides a graphical interface for all model definition operations. A new model can be assembled through the following steps, entirely without a text editor:

Comparison with existing multibody dynamics tools

The principal differentiating characteristics of OpenMBD relative to the existing tools described in the Introduction are its combination of Python nativity, zero-compilation installation, GUI accessibility, and biomechanics-specific features in a single package. OpenSim remains the most mature and widely validated open-source platform for musculoskeletal systems with a large model repository, extensive documentation, and an active user community.¹ OpenSim is the appropriate tool of choice for studies requiring detailed muscle force estimation, motion-capture-driven inverse kinematics, or integration with the broader OpenSim ecosystem. OpenMBD does not replicate this functionality and is not intended to replace OpenSim in those applications. The relative advantage of OpenMBD lies in scenarios where the analyst requires a self-contained forward dynamics simulation with explicit contact mechanics and requires rapid simulating of novel body geometries. MuJoCo offers superior simulation speed, particularly for reinforcement learning workloads requiring thousands of parallel rollouts and a more sophisticated contact solver based on convex optimisation.³ However, MuJoCo is robotics and control focused and does not provide physiological ROM constraints or a dedicated impact biomechanics GUI. OpenMBD’s penalty-contact approach, while less computationally sophisticated than MuJoCo’s convex solver, is transparently parameterisable through user-defined force-penetration curves, which aligns well with the empirical characterisation of soft tissue material properties from cadaveric impact data. PyBullet and RBDL are appropriate for robotics applications requiring high-speed simulation but lack physiological joint limits, anatomically parameterised body segment inertia, and a biomechanics-oriented user interface.^8,9 PyDy provides a powerful symbolic modelling framework but is intended for users comfortable deriving and manipulating equations of motion manually.⁷ OpenMBD occupies a complementary niche: it abstracts the dynamics formulation entirely and exposes only the model geometry, initial conditions, and physics parameters to the user.

Principal strengths

Minimal installation barrier: The three-package dependency (NumPy, Matplotlib, Pillow) and absence of any compilation step represent a substantially lower barrier to adoption than C++ − based alternatives. A researcher with a standard Python installation can be running simulations within minutes of downloading the repository.

Transparent, auditable physics: Every component of the dynamics formulation is implemented in readable Python with inline documentation. This transparency supports pedagogical use and facilitates modification and extension.

Open, human-readable model format: The JSON model format is self-documenting and editable in any text editor or the browser-based model editor.

Human body models: The adult male and female models are parameterised from CDC NHANES 2021–2023 population data,¹⁷ with segment masses and inertia tensors derived from the de Leva¹⁸ regression equations and sex-specific geometric scaling from ANSUR anthropometric surveys.²² This grounding in population-representative anthropometry makes the models appropriate for population-level biomechanics research.

Structured, timestamped output: The automatic export of three output files on completion of every run ensures that all simulation results are immediately available for analysis and that runs are uniquely identified by timestamp, preventing accidental overwriting. The columnar CSV format is compatible with Python, R, MATLAB, and Excel without format conversion.

Limitations

Rigid Body Assumption: All body segments in OpenMBD are treated as perfectly rigid. Biological tissues exhibit viscoelastic deformation under load (e.g., bone flexes, soft tissue deforms, and joint cartilage compresses), all of which are absent from the rigid body model. This approximation is standard in gross-motion multibody dynamics simulations¹⁴ and is generally acceptable for studying segment kinematics and gross contact forces, but it means that OpenMBD is not appropriate for applications requiring internal stress distribution, tissue deformation fields, or bone fracture prediction.

Ellipsoid Contact Geometry: The exclusive use of ellipsoid geometry simplifies contact detection and provides smooth surface normals, but it cannot capture the concave anatomical features of real body segments that influence contact area, pressure distribution, and joint congruency. The current ellipsoid approximation is most appropriate for gross-motion and impact studies where the precise contact geometry is a secondary concern relative to the overall trajectory and force magnitude.

Explicit Euler Integration and Time Step Sensitivity: The Symplectic Euler integration scheme, while exhibiting superior energy conservation properties compared with explicit Euler for oscillatory systems, remains a first-order method and is sensitive to the choice of time step. The default Δt = 0.0001 s was selected to satisfy the explicit-Euler stability criterion for the stiffest joint in the human model and to resolve the fast dynamics of ground impact events.

Single-Threaded Integration: The integration loop runs on a single CPU thread. This architecture is appropriate for the interactive, exploratory use case for which OpenMBD is designed, but it does not exploit multi-core parallelism. For applications requiring large numbers of simulation runs (e.g., Monte Carlo sensitivity analyses, parameter optimisation, or training data generation for machine learning), a batch simulation mode without the GUI overhead would be desirable.

Contact Model Parameterisation: The piecewise-linear force-penetration curves that define contact stiffness for each ellipsoid must be specified by the user and are not automatically calibrated to measured tissue properties. The curves distributed with the models are based on published impact data and provide physically plausible results for impact scenarios, but they have not been formally validated against instrumented impact experiments for the full range of loading conditions. Researchers requiring quantitatively accurate contact force predictions should treat simulation outputs as order-of-magnitude estimates and calibrate force-penetration curves against segment-specific experimental data where possible.

Controls for bias and sources of variability

Several design decisions in OpenMBD were made specifically to control for known sources of bias and variability in multibody dynamics simulation:

Initial overlap suppression: The recording of initial ellipsoid penetration offsets at t = 0 and subtraction of these offsets from contact force calculations prevents the resting geometry of the model from generating spurious forces from the very first timestep.

Adjacency exclusion: The exclusion of contact pairs within two kinematic hops prevents neighbouring body segments from generating contact forces through their designed geometric overlap, analogous to the self-collision exclusion masks used in character animation and robotics simulation.⁸

Inertia tensor regularisation: Eigen decomposition with negative eigenvalue clamping and diagonal mass matrix regularisation prevent singular or indefinite mass matrices that would produce non-physical accelerations.

Future development

Several extensions to OpenMBD are planned or identified as high-priority for future development: