4. Overlay Maximum Principal Strain Curves for Brain Tissue

Objective and Data Model

This case builds a strain-overlay workflow for brain tissue in a head-impact CAE model. The task is to extract the maximum principal strain (MPS) time history from a set of brain tissue elements, then overlay those histories across several impact directions so that phase (timing of the peak) and amplitude (peak magnitude) differences become directly comparable. The output is a single overlay plot plus a small summary table of peak strain and time-to-peak per direction.

Maximum principal strain is the largest eigenvalue of the Green–Lagrange (or small) strain tensor at each integration point. In an explicit FE solver such as LS-DYNA, Abaqus/Explicit, or Radioss, the element-level MPS is typically written to a binary or ASCII history file per element (e.g., elout, ELOUT, or a .csv export). For overlay purposes we treat each direction as one curve: the envelope of MPS over the selected element set at each time step, i.e. MPS_env(t) = max over elements of MPS_e(t). This envelope is the quantity most often compared against injury thresholds, and it is what we overlay.

Inputs and Assumptions

• Four impact directions: frontal, lateral_left, lateral_right, oblique_45.
• Each direction has one simulation with a shared element set brain_roi (e.g., 1,200 elements in the cerebrum).
• Time histories are exported at the solver output interval; here we assume a uniform 0.1 ms step for clarity, but the extraction code does not require uniformity.
• Strain is dimensionless; peaks are expected in the 0.05–0.35 range for a moderate-severity impact.
• All curves share the same physical time origin (impact at t = 0). Onset alignment is out of scope here and is handled in a later case.

Extraction: From Element Histories to an Envelope

The extraction step reads one file per direction, each containing columns time and one column per element (or a long-format table with time, element_id, mps). The function below handles the long format, which is the most common export shape, and returns a time vector plus the envelope. It also records the element that produced the peak, which is useful for locating the hotspot.

import numpy as np
import pandas as pd

def extract_envelope(csv_path, time_col="time", elem_col="element_id", val_col="mps"):
    df = pd.read_csv(csv_path)
    # pivot to wide: rows = time, columns = element
    wide = df.pivot_table(index=time_col, columns=elem_col, values=val_col, aggfunc="max")
    wide = wide.sort_index()
    t = wide.index.to_numpy()
    env = wide.max(axis=1).to_numpy()          # envelope over elements
    peak_elem = wide.idxmax(axis=1).to_numpy() # element id at each step
    return t, env, peak_elem

directions = ["frontal", "lateral_left", "lateral_right", "oblique_45"]
curves = {}
for d in directions:
    t, env, peak_elem = extract_envelope(f"mps_{d}.csv")
    curves[d] = {"t": t, "env": env, "peak_elem": peak_elem}

Two practical checks before overlaying: (1) confirm that t is monotonically increasing and starts near zero; (2) confirm that the envelope is not dominated by a single spurious element (e.g., a distorted element at the ROI boundary). A quick sanity check is to compare the envelope peak against the 99th percentile of all element peaks — a large gap usually indicates one bad element.

Overlay and Phase/Amplitude Comparison

With envelopes in hand, the overlay is a straightforward multi-line plot. The value of the overlay comes from reading two quantities per curve: the peak amplitude and the time-to-peak. Plotting them together makes phase differences (a lateral impact peaking later than a frontal one, for example) visible at a glance.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 4.5))
summary = []
for d, c in curves.items():
    ax.plot(c["t"], c["env"], label=d, linewidth=1.6)
    i = int(np.argmax(c["env"]))
    summary.append({"direction": d, "peak_mps": c["env"][i], "t_peak_ms": c["t"][i]})
    ax.plot(c["t"][i], c["env"][i], "o", markersize=5)

ax.set_xlabel("Time (ms)")
ax.set_ylabel("Max principal strain (envelope)")
ax.set_title("Brain tissue MPS envelope by impact direction")
ax.legend(frameon=False)
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig("mps_overlay.png", dpi=150)

print(pd.DataFrame(summary).to_string(index=False))

The printed table is the quantitative companion to the plot. A representative example output (illustrative values, not measured) might look like:

Read the overlay in two passes. First, amplitude: the lateral directions exceed the frontal peak, which is consistent with the larger rotational component in lateral loading. Second, phase: the lateral peaks arrive roughly 1.5 ms later than the frontal peak, indicating a slower strain build-up. The oblique case sits between the two, as expected from its mixed loading path.

Visualizing the Overlay Logic

The diagram below shows how per-element histories collapse into one envelope per direction, and how the four envelopes are then compared on a shared time axis. It also marks the peak point on each envelope, which is the quantity carried into the summary table.

Verification and Common Failure Modes

Before trusting the overlay, verify three things. First, the envelope peak should equal the maximum of the per-element peaks — if not, the pivot or aggregation is wrong. Second, the peak element id should be stable across a few time steps around the peak; a rapidly jumping id suggests numerical noise. Third, the overlay should reproduce the summary table exactly when you re-read the peak markers.

• Flat or zero curves: the value column may be a different strain measure (e.g., effective plastic strain) or the ROI filter returned no elements.
• Spiky envelope: a single distorted element dominates; inspect peak_elem and consider excluding boundary elements.
• Misaligned peaks: if all curves peak at t ≈ 0, the time column may be in solver cycles rather than milliseconds.
• Overlapping curves: identical element sets or a copy-paste error in the direction loop.

Deliverable: Overlay Summary Sheet

The reusable output of this case is a compact summary sheet that pairs the overlay plot with the peak table. Fill one row per direction; the columns are the minimum needed to compare phase and amplitude.

With the envelope extraction, the overlay plot, and this sheet in place, the phase and amplitude differences across impact directions are explicit and reproducible. The next comparison step — aligning onsets when simulations start at different times — can build directly on these envelopes without re-extracting the element histories.