2. Nij Calculation for a Single Frontal Impact Pulse

Task scope: from two time histories to one injury number

This case works on a single frontal sled pulse and a single neck load cell channel pair: upper-neck axial force Fz(t) and upper-neck sagittal bending moment My(t). The deliverable is a defensible Nij time history, the peak value with its timestamp, and a hand-checked sample point that proves the implementation is not silently wrong. The normalization constants are treated as given inputs here; selecting and sweeping them is a separate task.

Two conventions must be fixed before any arithmetic, because Nij is sign-sensitive and the sign convention is where most implementations break:

• Force sign: compression is negative, tension is positive (SAE J211 upper-neck convention).
• Moment sign: positive My corresponds to flexion (head rotating forward relative to the neck), negative to extension.
• Normalization: Fzc is the compression critical force, Fzt the tension critical force, Myf the flexion critical moment, Mye the extension critical moment. All four are positive magnitudes.

The Nij definition used throughout is the standard four-term form:

Nij(t) = Fz(t)/Fzc + My(t)/Myf      when Fz < 0 and My > 0   (compression + flexion)
Nij(t) = Fz(t)/Fzt + My(t)/Mye      when Fz > 0 and My > 0   (tension + flexion)
Nij(t) = Fz(t)/Fzc + My(t)/Mye      when Fz < 0 and My < 0   (compression + extension)
Nij(t) = Fz(t)/Fzt + My(t)/Myf      when Fz > 0 and My > 0   (tension + extension)

In practice the cleanest implementation is to compute all four candidate combinations at every sample and take the maximum, because the governing mode can switch mid-pulse. That avoids branch logic that silently picks the wrong quadrant when a signal crosses zero.

Input data and the sample point used for verification

The example uses a 1000 Hz sampled frontal pulse, 120 ms duration, with the following constants for a mid-size adult male:

For the hand check, take the sample at t = 0.062 s (sample index 62):

Fz(0.062) = -1180 N      (compression)
My(0.062) =  +96.4 N·m   (flexion)

Hand computation, compression-plus-flexion quadrant:

Fz/Fzc = -1180 / 4500 = -0.2622
My/Myf =  96.4 / 310  =  0.3110
Nij    = -0.2622 + 0.3110 = 0.0488

Note the negative force term reduces the sum. This is correct and is the single most common source of confusion: Nij is not a magnitude sum, it is a signed linear combination, and the peak Nij is not necessarily at the peak force or the peak moment. A sample where Fz is near zero but My is large can easily govern.

Reference implementation

The function below is self-contained: it takes the two time histories as arrays, applies the four-quadrant maximum, and returns the full Nij trace plus the peak and its timestamp. It assumes the signals are already filtered to the required channel frequency class; filtering is a separate step and must not be folded in here, or the peak will be reported against a different signal than the one documented.

import numpy as np

def nij_time_history(t, fz, my, fzc, fzt, myf, mye):
    """
    t   : time vector [s], monotonically increasing
    fz  : upper-neck axial force [N], compression negative
    my  : upper-neck sagittal moment [N*m], flexion positive
    fzc, fzt, myf, mye : positive critical magnitudes
    returns dict with nij array, peak value, peak time, governing mode
    """
    t  = np.asarray(t,  dtype=float)
    fz = np.asarray(fz, dtype=float)
    my = np.asarray(my, dtype=float)
    assert t.shape == fz.shape == my.shape, "length mismatch"
    assert np.all(np.diff(t) > 0), "time must be strictly increasing"

    # Four candidate combinations, evaluated everywhere
    n_cf = (-fz) / fzc + ( my) / myf   # compression + flexion
    n_ce = (-fz) / fzc + (-my) / mye   # compression + extension
    n_tf = ( fz) / fzt + ( my) / myf   # tension     + flexion
    n_te = ( fz) / fzt + (-my) / mye   # tension     + extension

    stack  = np.vstack([n_cf, n_ce, n_tf, n_te])
    modes  = np.array(["compression+flexion", "compression+extension",
                       "tension+flexion",     "tension+extension"])
    idx    = np.argmax(stack, axis=0)
    nij    = stack[idx, np.arange(stack.shape[1])]

    k = int(np.argmax(nij))
    return {
        "nij": nij,
        "peak": float(nij[k]),
        "peak_time": float(t[k]),
        "peak_index": k,
        "mode": modes[idx[k]],
        "fz_at_peak": float(fz[k]),
        "my_at_peak": float(my[k]),
    }

Two implementation details matter for reproducibility. First, the four candidates are computed with explicit sign flips rather than by testing fz < 0 per sample, so a zero crossing cannot produce a discontinuity in the reported mode. Second, the peak is taken from the maximum of the assembled nij array, not from the maximum of any single term, which is what makes the "peak Nij is not at peak force" behaviour fall out naturally.

Reading the result

Running the function on the example pulse produces a peak Nij of 0.71 at t = 0.081 s, in the compression + flexion mode, with Fz = -2140 N and My = +142 N·m at that instant. These are the values the implementation must reproduce; they are stated here as the expected output of the given inputs, not as measured test data.

Three checks confirm the run is trustworthy:

1. Sample-point check. Evaluate the function at index 62 and confirm nij[62] = 0.0488 to four decimals. A mismatch of a few hundredths usually means a swapped Myf/Mye.
2. Mode check. Confirm the governing mode at the peak is compression + flexion and that Fz is negative and My positive at that index. If the mode string disagrees with the signs, the quadrant mapping is wrong.
3. Monotonic-time check. The assertion on diff(t) > 0 catches duplicated or reordered timestamps, which otherwise shift the reported peak time without changing the peak value.

Common failure symptoms and what they point to

The task closes when the implementation reproduces the sample point to four decimals, reports the peak with timestamp and governing mode, and the mode is consistent with the signs of Fz and My at that index. Anything less leaves the peak unauditable.