3. Compute and Overlay Cumulative Strain Damage Measure (CSDM) Curves

From Principal Strain History to a Cumulative Distribution

Cumulative Strain Damage Measure (CSDM) answers a different question than a peak-strain readout. Instead of reporting the single worst element at the single worst instant, CSDM describes how much tissue volume has experienced a principal strain above a threshold at any point during the event. The result is a monotonically decreasing curve: at a low threshold nearly 100% of the modeled volume qualifies, and as the threshold rises the surviving fraction falls toward zero. Overlaying these curves across simulations lets you compare configurations by the shape of the whole distribution rather than by one scalar.

The definition used here follows the standard Takhounts-style formulation. For each element e and each threshold level ε*, an element is counted as damaged if its maximum principal strain over the entire event exceeds the threshold:

D(e, ε*) = 1  if  max_t ε₁(e, t) > ε*
          = 0  otherwise

CSDM(ε*) = Σ_e [ D(e, ε*) · V(e) ] / Σ_e V(e)

Two implementation details dominate accuracy. First, the maximum must be taken over time before thresholding — CSDM is a cumulative measure, so a brief excursion above threshold counts permanently. Second, the weighting is by element volume, not element count; a mesh with graded refinement will otherwise bias the curve toward the densely meshed region.

Building the Strain-History Inputs

Assume each simulation has been exported as a per-element table of maximum principal strain versus time. A practical working format is one row per (element, time) pair with columns elem_id, time, eps1, and volume. The volume column is constant per element and can be joined from the model's element table. The following Python fragment computes CSDM for one simulation and returns a curve ready for overlay:

import numpy as np
import pandas as pd

def csdm_curve(hist: pd.DataFrame,
               thresholds=np.arange(0.0, 0.51, 0.01)) -> pd.DataFrame:
    """hist columns: elem_id, time, eps1, volume"""
    # 1. Peak principal strain per element over the whole event
    peak = (hist.groupby("elem_id", as_index=False)
                .agg(eps1_max=("eps1", "max"),
                     volume=("volume", "first")))

    total_vol = peak["volume"].sum()
    if total_vol <= 0:
        raise ValueError("Total volume is zero; check the volume column.")

    # 2. Volume fraction exceeding each threshold
    rows = []
    for thr in thresholds:
        damaged = peak.loc[peak["eps1_max"] > thr, "volume"].sum()
        rows.append((thr, damaged / total_vol))
    return pd.DataFrame(rows, columns=["threshold", "csdm"])

Note the strict inequality >. Using >= shifts the curve by one threshold step and makes results disagree with reference implementations; pick one convention and document it. The threshold grid of 0.00–0.50 in 0.01 steps is a common reporting choice for head-injury work, but any monotonic grid is valid as long as every simulation uses the identical grid — otherwise the overlay is comparing different abscissae.

Overlaying Multiple Configurations

Because every curve shares the same threshold vector, the overlay is a simple wide-format pivot. The function below accepts a dictionary of named simulations and produces a single frame with one column per configuration:

def overlay_csdm(sim_histories: dict, thresholds=np.arange(0.0, 0.51, 0.01)):
    curves = {}
    for name, hist in sim_histories.items():
        c = csdm_curve(hist, thresholds).set_index("threshold")["csdm"]
        curves[name] = c
    return pd.DataFrame(curves)   # index = threshold, columns = configs

Plotting the result with threshold on the x-axis and CSDM on the y-axis gives the familiar descending family of curves. A configuration that shifts its curve to the right at a given CSDM level (say 0.10) is exposing less volume to high strain and is the safer of the two at that operating point. Report the crossing threshold explicitly rather than eyeballing it:

def threshold_at_csdm(overlay: pd.DataFrame, level=0.10):
    out = {}
    for col in overlay.columns:
        y = overlay[col].values
        x = overlay.index.values
        # linear interpolation on the descending curve
        out[col] = float(np.interp(-level, -y, x))
    return out

The negation trick works because np.interp requires increasing x; flipping the sign of both axes preserves the interpolation while satisfying that requirement.

Reading the Overlay

The SVG below shows three representative configurations plotted on a common threshold axis. The vertical guide at CSDM = 0.10 marks the comparison point used by threshold_at_csdm; the horizontal offsets between curves at that level are the quantities you report.

In this illustrative overlay, Config C reaches the 10% volume-fraction level at a threshold of roughly 0.19, Config B at 0.27, and Config A at 0.34. The ordering is consistent across the mid-range of the curve, which is the region that usually drives injury-risk decisions; near the tails the curves converge and the ranking becomes less meaningful. Always state the CSDM level at which you compare, because two curves can cross and reverse their ordering at different volume fractions.

Common Failure Modes