3. Physics-aware message passing

A conventional graph neural network treats every edge as a generic channel: it concatenates node features, applies a shared MLP, and aggregates the result. For crash simulation this abstraction is too weak. The quantities that actually govern a transient response — momentum flux across a contact interface, the impulse delivered during a short impact window, and the reaction forces that keep a joint or a rigid body on its kinematic manifold — are all directed, signed, and physically conserved. A physics-aware message-passing layer must therefore be designed so that the edge function computes a physical transfer, the aggregation step respects conservation, and the node update integrates that transfer over a time step. This section shows how to build such a layer concretely.

Reference conventions used throughout

All equations below use a right-handed Cartesian frame. For each node i we store position x_i, velocity v_i, mass m_i, and an accumulated external force f_i^ext. For each directed edge i→j we store a contact normal n_ij pointing from i toward j, a signed gap g_ij (positive when separated), and a friction coefficient μ_ij. Impulses are expressed in N·s and are positive when they push the two bodies apart. Keeping this sign convention fixed across the edge function, the aggregation, and the node update is essential; a sign flip between layers is the single most common source of non-physical drift in learned crash surrogates.

Edge update: from geometry to a contact impulse

The edge function should not predict an arbitrary vector. It should predict the scalar normal impulse λ_n and the tangential impulse vector λ_t, then reconstruct the physical impulse from the contact frame. A practical formulation is a small MLP over physically meaningful inputs, followed by a projection that enforces the complementarity conditions of unilateral contact:

# edge i -> j, per time step dt
g   = dot(x_j - x_i, n_ij) - r_i - r_j        # signed gap, r = contact radius
gdot = dot(v_j - v_i, n_ij)                    # approach velocity
z   = [g, gdot, mu_ij, m_i, m_j, |f_ext_i|, |f_ext_j|]
raw = MLP_edge(z)                              # outputs [lambda_n_hat, lambda_t_hat(2)]
lambda_n = softplus(raw[0]) * (g < 0)          # no tensile impulse
lambda_t = clamp_norm(raw[1:3], mu_ij * lambda_n)   # Coulomb cone
J_ij = lambda_n * n_ij + lambda_t               # impulse on j from i

Three design choices matter here. First, the gap and approach velocity are computed from the current state rather than learned, so the network only has to model the constitutive response (how stiff the contact is, how much energy is dissipated), not the kinematics. Second, the softplus gate and the Coulomb clamp are hard physical priors: they guarantee that contacts never pull and that friction never exceeds the cone, even early in training when the MLP output is noisy. Third, the edge message passed to node j is the impulse itself, J_ij, not a latent embedding — this is what makes the aggregation step interpretable and conservative.

Aggregation: antisymmetry and momentum conservation

Because J_ij is the impulse that i exerts on j, the impulse on i is exactly −J_ij. Aggregation must therefore be an antisymmetric sum, not a symmetric mean or a max-pool:

J_i = sum_{j in N(i)} (-J_ij)      # impulse received by node i
J_j = sum_{i in N(j)} ( J_ij)      # impulse received by node j

With this rule, summing over all nodes gives Σ_i J_i = 0 identically, so total linear momentum is conserved by construction regardless of what the edge MLP outputs. This is the key advantage of physics-aware message passing over a black-box GNN: the network can only redistribute momentum, never create it. For angular momentum, the same argument holds if the edge also carries the lever arm r_ij from the contact point to each node's center of mass and the node update applies the corresponding torque r_ij × J_ij.

Node update: semi-implicit integration with constraint projection

The node update is where the transient dynamics actually advance. A stable and physically consistent choice is a semi-implicit (symplectic Euler) step followed by a constraint projection for joints and rigid links:

v_i^{n+1/2} = v_i^n + dt * (f_ext_i + f_constraint_i) / m_i
v_i^{n+1}   = v_i^{n+1/2} + J_i / m_i
x_i^{n+1}   = x_i^n + dt * v_i^{n+1}
# constraint projection (e.g. ball joint between i and k)
C = x_i - x_k - L0
v_i^{n+1} -= (m_k/(m_i+m_k)) * C / dt
v_k^{n+1} += (m_i/(m_i+m_k)) * C / dt

The constraint force f_constraint can itself be learned — for example, an MLP over the constraint violation C and its rate — but it must be applied through the same mass-weighted projection so that the constraint manifold is restored without injecting energy. In practice, the node update is the only place where dt appears; the edge function is time-step agnostic. This separation lets you train on one dt and deploy at another, provided the impulse magnitudes are rescaled consistently.

Worked example: two-body impact with a learned contact law

Consider a minimal reproducible case: a 2 kg block (node i) moving at 5 m/s toward a 3 kg block (node j) at rest, with dt = 1 ms, μ = 0.3, and a contact radius of 0.02 m each. At the first step where g < 0, suppose the edge MLP outputs λ_n = 4.0 N·s and λ_t = 0.8 N·s along the tangential direction. The aggregation gives J_i = −4.0n − 0.8t and J_j = +4.0n + 0.8t. The node update then yields v_i = 5 − 4.0/2 = 3.0 m/s and v_j = 0 + 4.0/3 ≈ 1.33 m/s along the normal. Total momentum before and after is 10 kg·m/s in both cases, confirming the antisymmetric aggregation did its job. The tangential impulse of 0.8 N·s is within the cone (0.3 × 4.0 = 1.2 N·s), so the clamp is inactive — a useful sanity check that the learned law is not saturating.

What to verify before trusting the layer

Three checks catch most implementation errors. First, run a single edge with a known impulse and confirm that ΣJ_i = 0 to machine precision. Second, drive a head-on elastic impact and confirm that the coefficient of restitution implied by the post-impact velocities matches the dissipation the edge MLP was trained to produce. Third, run a constrained pendulum for several thousand steps and confirm that the constraint violation C stays bounded rather than drifting — a drifting C means the projection is not mass-weighted correctly. Once these pass, the layer can be stacked and trained end-to-end on transient crash data, with the physical priors carrying the burden that a purely data-driven model would otherwise have to learn from scratch.