3. Compute HIC15 and HIC36 from Acceleration Time Histories

Why HIC15 and HIC36 Must Be Computed Together

The Head Injury Criterion (HIC) is defined as the maximum of the time-windowed integral of resultant head acceleration:

HIC = max over (t1, t2) of [ (t2 - t1) * ( (1/(t2-t1)) * ∫ a(t) dt from t1 to t2 )^2.5 ]

where a(t) is the resultant translational acceleration at the head center of gravity expressed in g, and the window length t2 - t1 is constrained to 15 ms for HIC15 and 36 ms for HIC36. Because the exponent 2.5 heavily penalizes long, high-magnitude plateaus, HIC15 and HIC36 can rank two impact pulses in opposite order: a short, sharp pulse may dominate HIC15 while a longer, lower-amplitude pulse dominates HIC36. Computing both on the same time base and overlaying the two window-scan curves is therefore the minimum defensible comparison for severity ranking across load cases.

Inputs, Units, and Preconditions

Before any window scan, the acceleration channel must satisfy three conditions that are easy to get wrong and expensive to debug later:

• Resultant, not component. HIC uses the vector magnitude a_res = sqrt(ax² + ay² + az²). If the source file stores three components, compute the resultant first; do not feed a single component into the integral.
• Units in g. LS-DYNA and Radioss typically output acceleration in mm/ms² or m/s². Convert with a_g = a_SI / 9.80665. A factor-of-9.8 error produces a HIC that is off by roughly 9.8^2.5 ≈ 300×, which is unmistakable but only if you check against a reference.
• Uniform time base. The sliding-window integral assumes constant dt. Resample to a fixed rate (commonly 10 kHz, i.e. dt = 0.1 ms) before scanning; otherwise the trapezoidal weights must be recomputed per window.

Sliding-Window Implementation

The naive double loop over all (t1, t2) pairs is O(N²) and, at 10 kHz over a 150 ms pulse (N = 1500), still runs in well under a second — but the vectorized cumulative-integral form is both faster and less error-prone. The key identity is that the mean acceleration over a window equals the difference of the cumulative integral divided by the window length:

import numpy as np

G = 9.80665

def hic(a_g, dt, window_s):
    """Compute HIC over a fixed window length.
    a_g     : resultant acceleration in g, uniformly sampled
    dt      : sample interval in seconds
    window_s: 0.015 for HIC15, 0.036 for HIC36
    returns (hic_value, t1, t2)
    """
    a_g = np.asarray(a_g, dtype=float)
    n = a_g.size
    # cumulative trapezoidal integral of acceleration (units: g*s)
    cum = np.concatenate(([0.0], np.cumsum(0.5 * (a_g[1:] + a_g[:-1]) * dt)))
    t = np.arange(n) * dt
    w = int(round(window_s / dt))          # window length in samples
    if w >= n:
        raise ValueError("window longer than record")
    # mean acceleration over each window = (cum[i+w] - cum[i]) / (w*dt)
    seg = cum[w:] - cum[:-w]
    mean_a = seg / (w * dt)
    hic_vals = (w * dt) * np.power(np.clip(mean_a, 0.0, None), 2.5)
    k = int(np.argmax(hic_vals))
    return float(hic_vals[k]), float(t[k]), float(t[k + w])

Two details matter for correctness. First, np.clip(mean_a, 0, None) guards against negative means (deceleration-only windows) that would otherwise raise a fractional power of a negative number. Second, the returned t1 and t2 are the window endpoints at the maximum, which are needed to annotate the overlay and to cross-check against the reference report.

Validating Against Reference Values

A HIC implementation is only trustworthy after it reproduces a known answer. The standard sanity check is a rectangular pulse: for a constant acceleration A held for exactly the window duration W, the mean equals A and HIC reduces to A^2.5 · W. Use a synthetic pulse to confirm the code path before touching CAE data:

dt = 1e-4
t = np.arange(0, 0.100, dt)
a = np.zeros_like(t)
a[(t >= 0.020) & (t < 0.035)] = 80.0   # 80 g for 15 ms

h15, t1, t2 = hic(a, dt, 0.015)
print(f"HIC15 = {h15:.1f}  window = [{t1*1e3:.1f}, {t2*1e3:.1f}] ms")
# expected: 80**2.5 * 0.015 = 80**2.5 * 0.015 ≈ 858.6

If the printed value deviates from 858.6 by more than rounding, the fault is almost always one of: wrong dt (units of ms vs s), a missing g conversion, or an off-by-one in the window sample count. Once the synthetic case passes, run the same function on the imported CAE channel and compare against the solver's own HIC output or the reference report; agreement within 1–2 % is normal and reflects filter and resampling differences.

Overlaying HIC15 and HIC36 Window-Scan Curves

Plotting only the scalar HIC values hides where in the pulse the maximum occurs. A more informative overlay plots the window-scan value as a function of window start time for both window lengths on shared axes. The two curves then reveal whether the severity is concentrated in a single sharp event (both curves peak at the same t1) or spread across a longer plateau (HIC36 peaks later and higher relative to HIC15).

import matplotlib.pyplot as plt

def scan(a_g, dt, window_s):
    n = a_g.size
    cum = np.concatenate(([0.0], np.cumsum(0.5*(a_g[1:]+a_g[:-1])*dt)))
    w = int(round(window_s/dt))
    seg = cum[w:] - cum[:-w]
    mean_a = np.clip(seg/(w*dt), 0.0, None)
    return np.arange(n-w)*dt, (w*dt)*np.power(mean_a, 2.5)

t15, h15c = scan(a, dt, 0.015)
t36, h36c = scan(a, dt, 0.036)

fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(t15*1e3, h15c, color="#1f77b4", lw=2, label="HIC15 scan")
ax.plot(t36*1e3, h36c, color="#d62728", lw=2, label="HIC36 scan")
ax.axhline(700, ls="--", color="gray", lw=1, label="AIS3 reference (700)")
ax.set_xlabel("Window start time t1 (ms)")
ax.set_ylabel("HIC value")
ax.set_title("HIC15 vs HIC36 window scan — Case 3")
ax.legend(); ax.grid(alpha=0.3)

The dashed reference line at 700 is the commonly cited AIS 3+ threshold for HIC36; keep it as a visual guide only and cite the actual regulatory or internal limit used by the program. When overlaying multiple load cases, keep the same axis limits and reference line so the comparison is honest.

Reading the Overlay and Ranking Load Cases

In the illustrative scan above, HIC15 peaks near 33 ms while HIC36 peaks near 41 ms, and the two peaks differ by roughly 8 ms. That offset is the diagnostic signal: when the HIC36 peak occurs noticeably later than the HIC15 peak, the pulse has a sustained plateau that the longer window captures but the shorter window truncates. When both peaks coincide, the severity is dominated by a single sharp transient and HIC15 is the more discriminating metric. Use this rule to decide which metric to report as primary for a given load case, and always report both.

For multi-case comparison, build a small summary table rather than eyeballing curves:

The values above are illustrative teaching data, not measured results. The point of the table is the ranking logic: LC-B would be judged less severe than LC-A by HIC15 alone but more severe by HIC36, which is exactly the failure mode that computing only one window length produces.

Common Failure Symptoms and Fixes

• HIC in the tens of thousands. Acceleration still in m/s² or mm/ms². Divide by 9.80665 and re-run.
• HIC15 > HIC36 by a large margin on every case. Usually a sign the 36 ms window is being truncated by a short record; check that the record length exceeds the window plus the pulse duration.
• Peak window endpoints land on the record boundary. The pulse is clipped; extend the simulation end time or the extraction range.
• Reference mismatch of 5–10 %. Filter state or resampling differs from the reference. Confirm CFC class and resample rate before assuming a code bug.