4. Reading CAE simulation output

Once a solver run completes, the raw output is rarely a tidy table of injury metrics. LS-DYNA writes d3plot state databases and nodout/elout ASCII files; PAM-CRASH produces DSY and TH files; Abaqus/Explicit emits .odb plus .fil; Radioss uses Annn animation and T01 time-history files. Each format stores nodal displacement, velocity, acceleration, element stress, and contact forces on its own grid of time steps. The job of this chapter is to turn those heterogeneous files into a clean, physically meaningful set of time histories that feed the injury criteria computed later — HIC, BrIC, chest deflection, Nij, and similar metrics.

The central idea is a two-stage pipeline: discovery (locate and index the result files, understand their time base and entity numbering) and extraction (pull specific channels, resample onto a common time axis, and attach metadata). Getting discovery wrong is the most common source of silent errors — a head acceleration curve sampled at 10 kHz and a chest deflection curve sampled at 1 kHz will produce a meaningless HIC if naively combined.

Understanding the solver time base before touching data

Every explicit solver advances with a stable time step Δt determined by the smallest element and the wave speed. LS-DYNA typically writes d3plot states at a user-controlled interval (e.g., every 0.5 ms) while nodout can be written at a different interval. The binout file is the modern consolidated binary container and is the preferred source for time histories because it preserves the exact solver time vector. Before extraction, always read the header to confirm:

  • the time unit (ms vs. s) — LS-DYNA defaults to ms, Abaqus to s;
  • the number of states and the actual time vector (not assumed uniform);
  • the entity IDs (node sets, part IDs, cross-section IDs) that correspond to your injury-relevant locations.

A practical convention used throughout this course: all internal time histories are stored in seconds, all accelerations in m/s², and all displacements in meters. Conversion happens once, at ingestion, and is recorded in metadata so downstream criteria never guess units.

Mapping injury locations to solver entities

Injury criteria are defined at anatomical landmarks, but the solver only knows node and element IDs. The bridge is a channel map: a small table that links each anatomical target to a solver entity and a component direction. For a frontal impact model this might look like:

ChannelAnatomical targetSolver entityComponentUnit (raw)
head_accHead CGNODE 100234 (rigid head part)resultantmm/ms²
chest_deflUpper sternumNODE 204511 vs. NODE 204590Y (local)mm
pelvis_accPelvis CGNODE 300112Xmm/ms²
neck_forceOccipital condyleXSECTION 5Fz, MyN, N·mm

Note the local coordinate system on the chest deflection channel: the sternum-to-spine direction is not the global Y axis in a rotated occupant model. The channel map must record the reference frame, and the extraction code must apply the same rotation used to define the local system. Mixing global and local components is a classic source of sign-flipped deflection curves.

Extraction with a reproducible Python reader

The following example uses lasso-python (the open-source LS-DYNA reader) to pull channels from a binout file. The same pattern applies to pyNastran for Nastran op2, pyvista for VTK-based solvers, and h5py for solvers that write HDF5.

import lasso
import numpy as np
import pandas as pd

reader = lasso.dyna.Binout("run_042/binout")
t = reader.read("nodout", "time")          # seconds, solver time vector

def channel(node_id, comp):
    data = reader.read("nodout", "node", node_id)
    return data[:, comp]                    # comp: 0=x,1=y,2=z

head = channel(100234, 0)                   # raw mm/ms^2
head_ms2 = head * 1.0                       # mm/ms^2 == m/s^2 numerically
chest = channel(204511, 1) - channel(204590, 1)   # mm

df = pd.DataFrame({"time_s": t,
                   "head_acc_ms2": head_ms2,
                   "chest_defl_m": chest * 1e-3})
df.to_parquet("channels/run_042.parquet")

Two details matter here. First, mm/ms² is numerically identical to m/s² (10⁻³ m / 10⁻⁶ s² = 10³ m/s² per unit, and the raw value already carries the 10⁻³), so the head channel needs no scaling — but this must be verified against the solver's unit system, not assumed. Second, the chest deflection is a relative displacement between two nodes, so it is computed at extraction time rather than stored as two separate channels; this keeps the anatomical meaning attached to the data.

Resampling onto a common time axis

Injury criteria such as HIC integrate over a 15 ms window and require a uniform time step. Solver output is often non-uniform (adaptive time stepping, output at event triggers). Resample every channel onto a common grid — typically 10 kHz for acceleration, 1 kHz for deflection — using linear interpolation, and record the original sample rate in metadata. Never upsample a 1 kHz deflection channel to 10 kHz and treat it as high-fidelity; instead, keep the native rate and let the criterion function handle the mismatch explicitly.

Common failure: reading nodout when the model was run with *DATABASE_BINARY_D3PLOT only. The nodout file may be empty or contain a single state. Check file size and state count before extraction; if empty, switch to binout or enable *DATABASE_ASCII_NODOUT in the deck.

End-to-end data flow

The diagram below shows how solver artifacts flow through discovery, channel mapping, extraction, and resampling into a normalized channel store that later chapters consume for HIC, BrIC, and chest criteria.

Solver artifacts d3plot / binout nodout / elout .odb / .fil T01 / DSY Discovery read header time vector state count entity IDs Channel map anatomy → node component / axis local frame unit (raw) Extraction lasso / pyNastran relative channels unit conversion metadata attach Resample common grid 10 kHz / 1 kHz linear interp Normalized channel store Parquet / HDF5 time_s, value, unit source, frame, rate Injury criteria HIC / BrIC chest deflection Nij Discovery → mapping → extraction → resampling → normalized store; criteria consume the store, never raw solver files.

Verification and typical failure symptoms

After extraction, run three cheap checks before trusting any channel. First, plot the time vector and confirm it is monotonic and spans the event window (e.g., 0–120 ms for a frontal impact). Second, compare the peak of a known channel against the solver's own glstat or nodout summary — a head acceleration peak that differs by more than a few percent usually means the wrong node or a unit error. Third, confirm the channel count matches the channel map; a missing channel should raise an error, not silently produce a zero-filled column.

Typical symptoms and their causes: a flat zero channel usually indicates an empty nodout or a node ID that belongs to a deleted part; a curve that starts at a large nonzero value suggests a coordinate offset or a relative-displacement channel computed against the wrong reference node; a curve with visible stair-stepping indicates the native sample rate is too low for the criterion window and the channel should be flagged rather than upsampled. Recording these checks as assertions in the extraction script turns them from one-off inspections into regression tests that catch deck changes between runs.

Practical rule: store the raw solver time vector alongside every resampled channel. When a criterion later produces an unexpected value, you can re-derive it at the native rate and isolate whether the error came from resampling or from the criterion implementation itself.

The output of this chapter is a normalized channel store — one Parquet or HDF5 file per run, with columns for time, value, unit, source entity, reference frame, and native sample rate. That store is the contract between solver output and every injury criterion computed downstream, and keeping it clean and self-describing is what makes the automated reporting pipeline reproducible across runs and models.