2. Multibody and finite element formulations in MADYMO

MADYMO is a hybrid solver: it lets you build a model from two fundamentally different mathematical descriptions and couple them inside a single simulation. The multibody (MB) formulation represents bodies as rigid bodies connected by kinematic joints, with contact handled through analytically defined surfaces such as ellipsoids. The finite element (FE) formulation represents continua as meshes whose nodes carry displacement degrees of freedom, with stress and strain computed from element formulations. Understanding where each formulation is appropriate — and how they exchange loads — is the single most important modelling decision in an occupant kinematics study.

Rigid bodies and the multibody state

In the MB formulation, each rigid body carries six degrees of freedom: three translational and three rotational. MADYMO expresses the body state through a position vector and an orientation, typically stored internally as a rotation matrix or a set of Euler parameters (quaternions), which avoids the gimbal-lock singularity of Euler angles. The equations of motion are assembled in the form:

M(q) · q̈ + C(q, q̇) = Q_ext + Q_joint + Q_contact

where q is the generalized coordinate vector, M the mass matrix, C the velocity-dependent term (centrifugal and Coriolis contributions), and the right-hand side collects external loads, joint constraint forces, and contact forces. Because bodies are rigid, no element integration is required — the cost per body is low, and stable time steps are large. This is why a full dummy model with dozens of segments runs in minutes rather than hours.

Each rigid body also carries inertial properties: mass, centre of mass location, and the inertia tensor. In MADYMO these are defined in the body's local coordinate system, and the orientation of that local frame relative to the global frame is what the solver integrates. A common modelling error is to define the inertia tensor in the global frame; the result is a body that rotates about the wrong axes and produces non-physical head or torso motion.

Joints as kinematic constraints

Joints remove degrees of freedom between two bodies. MADYMO provides revolute, cylindrical, spherical, translational, universal, and free joints, plus bracket joints that apply a user-defined stiffness. A revolute joint, for example, constrains three translations and two rotations, leaving one rotational degree of freedom. The joint constraint is enforced either by Lagrange multipliers or by a penalty formulation; the penalty approach is more robust for contact-dominated problems but introduces a small compliance that must be stiff enough not to distort kinematics.

Joint stops and friction are defined as part of the joint card. A knee joint with a flexion stop at 120° and a damping coefficient will absorb energy at the stop rather than transmitting an impulsive load to the femur — a detail that materially changes lower-extremity injury metrics in frontal impact. When you review a model, check that every joint's range of motion matches the anatomical range you intend to represent; an unconstrained joint will let the dummy hyperextend and inflate neck injury criteria.

Ellipsoids and contact

Contact in the MB formulation is handled by ellipsoids (and, in some versions, planes, cylinders, and finite-element surfaces). An ellipsoid is defined by three semi-axes and a local frame attached to a rigid body. When two ellipsoids overlap, MADYMO computes the penetration depth along the contact normal and applies a force:

F_n = k · δ^p + c · δ̇   (penetration δ, stiffness k, damping c)

The stiffness and damping are user inputs, and the exponent p is typically 1 to 2. The contact force is applied at the contact point on each body, so it generates both a force and a moment. Ellipsoid contact is cheap and stable, which makes it ideal for dummy-to-airbag, dummy-to-seat, and dummy-to-steering-wheel interactions. Its limitation is geometric fidelity: an ellipsoid cannot represent a sharp edge or a thin shell, so a seat belt modelled as an ellipsoid will not capture the local pressure distribution that drives rib fracture.

Finite element formulation

The FE part of MADYMO uses shell, solid, membrane, and beam elements. Shells are the workhorse for airbags, seat belts, and sheet-metal structures; solids are used for foam and soft tissue. The element formulation may be updated Lagrangian with explicit time integration, which is conditionally stable: the time step is limited by the Courant condition,

Δt ≤ L_min / c,   c = sqrt(E / ρ)

where L_min is the smallest element characteristic length and c the wave speed. A 5 mm shell in steel with E = 210 GPa and ρ = 7800 kg/m³ gives c ≈ 5190 m/s and Δt ≈ 0.96 µs. That is why FE models are expensive: a 100 ms event needs roughly 100,000 steps. Mass scaling can raise the stable step, but it adds artificial inertia and must be checked against the physical response.

Material models matter as much as mesh. Airbags use a fabric model with orthotropic stiffness and a porous leakage term; seat belts use a 1D membrane with a defined force–strain curve; foams use a crushable foam model with a plateau stress. If you assign a linear elastic model to a foam, the computed pelvis deceleration will be far too high because the plateau that absorbs energy is missing.

Coupling MB and FE

Coupling is the mechanism that makes MADYMO a hybrid solver. There are three practical coupling paths:

• MB body to FE surface: an ellipsoid on a rigid dummy segment contacts an FE belt or airbag. The contact force computed from the ellipsoid penetration is distributed to the FE nodes.
• FE node to MB body: an FE node is tied to a rigid body, so the body follows the node. This is used for belt anchor points and for attaching an FE airbag to the steering wheel.
• FE to FE: two FE surfaces contact through a penalty or Lagrange contact algorithm, used for airbag-to-occupant and belt-to-clothing interactions.

The coupling interface must be checked for consistent units and for a stable force transfer. A frequent symptom of a bad coupling is a spike in contact force at the first time step, caused by an initial penetration between the ellipsoid and the FE surface. The fix is to position the ellipsoid with a small initial gap, or to use the solver's initial-penetration removal option.

Choosing the right formulation

The decision is driven by the physics you need to capture. Use MB when the dominant behaviour is gross rigid-body motion and the contact geometry can be approximated: dummy kinematics, restraint system layout studies, parameter sweeps, and real-time or near-real-time applications. Use FE when local deformation, stress, or energy absorption governs the injury metric: airbag deployment, belt loading on the chest, foam crush in the pelvis, and structural intrusion.

A practical rule for occupant safety models is to keep the dummy as MB and the environment as FE where it deforms. This gives you the speed of MB for the occupant and the fidelity of FE for the restraint. The coupling interface then becomes the critical review point: verify that the ellipsoid-to-FE contact stiffness is not so high that it produces a numerical spike, and not so low that the dummy penetrates the belt. A useful check is to plot contact force versus penetration for the first 10 ms; a smooth rise indicates a well-conditioned interface, while an oscillating or spiking curve indicates a stiffness or time-step problem.

Finally, remember that the two formulations share the same global time integration loop. The MB solver advances the rigid-body state, the FE solver advances the nodal state, and the coupling module exchanges forces at each step. If the FE time step is much smaller than the MB step, the solver subcycles the FE part; if the MB step is smaller, the coupling is evaluated more frequently. Monitoring the ratio of coupling force to inertial force over time is a reliable way to confirm that the interface is physically meaningful rather than numerically dominant.