4. Reading CAE simulation results

In an automated injury-report pipeline, the simulation solver is only half the story. The other half is the reader: the code that opens whatever the solver wrote, understands its structure, and pulls out the handful of channels that actually matter for injury assessment. CAE post-processing rarely produces a single tidy table. A typical crash or occupant-safety study emits CSV exports for nodal time histories, JSON metadata for model and load-case definitions, HDF5 containers for large multi-channel result sets, and raw binary dumps for solver-internal arrays. Each format has its own access pattern, its own failure modes, and its own cost. This chapter builds a practical reader layer that normalizes all of them into one in-memory structure: a time vector plus named channels, ready for injury metrics such as HIC, NIC, chest deflection, and femur load.

What "injury-relevant" actually means for the reader

Before writing a single parser, decide what the downstream metrics require. Injury criteria are defined over specific signals and specific windows:

• HIC needs the resultant head acceleration time history, sampled fast enough (typically 10 kHz or higher) to resolve the 36 ms window.
• NIC needs upper-neck and lower-neck x-axis accelerations, plus the relative displacement derived from them.
• Chest deflection needs the deflection channel in millimetres, with a clear sign convention (positive = compression).
• Femur load needs the axial force channel in kN, with the compression direction defined.

This list drives the reader's contract: it must preserve sample rate, units, sign convention, and channel identity. A parser that silently drops the time base or rescales units will produce plausible-looking but wrong HIC values. Treat metadata as first-class data, not decoration.

Format-by-format access patterns

CSV is the simplest and the most fragile. A typical export has a header row of channel names and a first column of time. Use pandas.read_csv with explicit dtype and na_values, and validate that the time column is strictly monotonic:

import pandas as pd
import numpy as np

def read_csv_history(path, time_col="time"):
    df = pd.read_csv(path, na_values=["", "NaN", "-9999"])
    t = df[time_col].to_numpy(dtype=float)
    if not np.all(np.diff(t) > 0):
        raise ValueError(f"Non-monotonic time base in {path}")
    return t, {c: df[c].to_numpy(dtype=float) for c in df.columns if c != time_col}

The -9999 sentinel is common in solver exports and must be mapped to NaN before any metric computation, otherwise it will dominate a peak search. Also check for a units row: some exporters write a second header line such as s, m/s^2, mm, kN. If present, strip it and record the units alongside each channel.

JSON usually carries metadata rather than bulk arrays: model version, load case, sensor locations, units, and sometimes short time histories. Parse it with the standard library and validate required keys explicitly:

import json

def read_json_meta(path):
    with open(path, "r", encoding="utf-8") as f:
        meta = json.load(f)
    required = {"load_case", "units", "channels"}
    missing = required - meta.keys()
    if missing:
        raise KeyError(f"Missing metadata keys: {missing}")
    return meta

Do not trust JSON numbers to be floats — some exporters write integers for what should be millimetres. Coerce with float() at the boundary and record the original type if it matters.

HDF5 is where large result sets live. Use h5py and read datasets lazily; never load a 2 GB acceleration array just to inspect its shape:

import h5py

def read_hdf5_history(path, group="/results", time_key="time"):
    with h5py.File(path, "r") as f:
        g = f[group]
        t = g[time_key][...]
        channels = {k: g[k][...] for k in g.keys() if k != time_key}
    return t, channels

HDF5 files are hierarchical, so the same channel name can appear under different groups (e.g. /results/head and /results/neck). Always resolve the full path, not just the leaf name, and keep the group path as part of the channel identifier. Also check dataset.attrs for units and sample rate — HDF5 is the one format where this metadata is reliably present, so use it.

Binary dumps are the least standardized. You need the layout: endianness, element type, record length, and channel order. A common pattern is a header of int32 counts followed by interleaved float64 records. Read with numpy.fromfile and reshape explicitly:

import numpy as np

def read_binary_history(path, n_channels, dtype="<f8"):
    raw = np.fromfile(path, dtype=dtype)
    if raw.size % n_channels != 0:
        raise ValueError("Record length mismatch")
    data = raw.reshape(-1, n_channels)
    t = data[:, 0]
    return t, {f"ch{i}": data[:, i] for i in range(1, n_channels)}

Binary readers are where silent corruption happens. Validate the record count against the header, check that the time column is monotonic, and compare a checksum or a known first/last value against the solver log. If the file was written on a big-endian machine, <f8 will produce garbage that still looks like floats — always confirm endianness from the export documentation.

Normalizing into a single time-history object

Once each format is parsed, convert everything into one structure so downstream injury metrics do not care where the data came from. A small dataclass works well:

from dataclasses import dataclass, field

@dataclass
class TimeHistory:
    time: np.ndarray
    channels: dict
    units: dict = field(default_factory=dict)
    source: str = ""
    sample_rate: float = 0.0

    def resample(self, rate):
        # linear interpolation onto a uniform grid
        t_new = np.arange(self.time[0], self.time[-1], 1.0 / rate)
        out = {k: np.interp(t_new, self.time, v) for k, v in self.channels.items()}
        return TimeHistory(t_new, out, self.units, self.source, rate)

Resampling matters because HIC and NIC are defined over fixed windows and assume a uniform time step. If the solver wrote variable-step output, interpolate onto a uniform grid before computing any metric, and record the original sample rate in the object so the report can state it.

Extracting the injury channels

With a normalized TimeHistory, extraction becomes a lookup with validation. Define a channel map per load case, then assert that every required channel exists and has the expected unit:

REQUIRED = {
    "head":  ("accel_resultant", "m/s^2"),
    "neck":  ("accel_x_upper", "m/s^2"),
    "chest": ("deflection", "mm"),
    "femur": ("force_axial", "kN"),
}

def extract_injury_channels(th, required=REQUIRED):
    out = {}
    for part, (name, unit) in required.items():
        if name not in th.channels:
            raise KeyError(f"Missing channel {name} for {part}")
        if th.units.get(name) != unit:
            raise ValueError(f"Unit mismatch for {name}: {th.units.get(name)} != {unit}")
        out[part] = th.channels[name]
    return out

The unit check is not pedantic. A chest deflection exported in metres instead of millimetres will pass every shape check and fail every injury threshold. Fail loudly at extraction time, not silently at report time.

Reading the result: what to verify before trusting a time history

After extraction, run a short sanity pass on each channel before it reaches the metric code. Check that the time vector starts near zero and ends at or beyond the event duration, that the peak lies inside the expected window (e.g. head acceleration peak within the first 150 ms), and that no channel is entirely NaN or constant. A constant channel usually means a failed sensor or a wrong column index; a channel peaking at the very first sample usually means a time-offset error. Log these checks per file so the automated report can flag suspect runs instead of publishing them.

A reproducible check on a small example

Suppose a load case exports head_accel.csv with columns time, accel_resultant at 10 kHz, and a companion case_meta.json declaring "units": {"accel_resultant": "m/s^2"}. Run the CSV reader, confirm the time vector is monotonic and spans 0–0.2 s, then confirm the peak of accel_resultant falls between 0.02 s and 0.15 s. If the peak sits at the first sample, suspect a time-offset or a column misalignment; if the channel is constant, suspect a failed sensor or a wrong column index. These are the same checks the automated report will run, so exercising them manually on one file before scaling to a batch is the fastest way to catch a reader bug.

The reader layer is now complete enough to feed the injury-metric stage: it opens four formats, normalizes them into one object, validates units and time bases, and extracts the channels that HIC, NIC, chest deflection, and femur load require. The next step is computing those metrics on the normalized histories and assembling them into the report structure.