3. Feature engineering for numerical failure diagnosis

Raw solver output is not directly consumable by a machine-learning model. A crash simulation produces gigabytes of nodal displacements, element stresses, contact forces, and solver log lines, but a classifier or regressor needs a compact, fixed-length, physically meaningful vector. Feature engineering is the bridge between the simulation database and the diagnostic model. In this chapter we build that bridge from three families of descriptors: time-series descriptors extracted from scalar and vector histories, spatial field statistics computed over the mesh, and solver convergence indicators parsed from the nonlinear iteration log. The goal is a feature vector that is stable across mesh sizes, robust to rigid-body motion, and sensitive to the specific failure mode we want to detect.

Why raw fields cannot be fed directly

A single LS-DYNA or Abaqus/Explicit run may contain 2–10 million nodes and 10,000+ time steps. Three problems block direct use:

• Variable dimensionality. Different models have different node and element counts, so a raw field is not a fixed-length vector.
• Mesh dependence. Peak stress at one node depends on local mesh density; the same physical event yields different numbers on a coarse and a fine mesh.
• Redundancy and noise. Most of the field is smooth and uninformative; the diagnostic signal is concentrated in a small region and a short time window.

Feature engineering solves all three by aggregating over space and time into a small set of scalars that carry the physics.

Time-series descriptors

For every scalar history s(t) — kinetic energy, internal energy, contact force, a nodal acceleration component, or a solver residual — we compute a descriptor set. The following Python function is representative of what runs in a preprocessing pipeline:

import numpy as np

def ts_descriptors(t, s):
    s = np.asarray(s, dtype=float)
    dt = np.diff(t)
    v = np.gradient(s, t)          # first derivative
    a = np.gradient(v, t)          # second derivative
    return {
        "peak":      s.max(),
        "t_peak":    t[np.argmax(s)],
        "mean":      s.mean(),
        "std":       s.std(),
        "rms":       np.sqrt(np.mean(s**2)),
        "skew":      float(((s - s.mean())**3).mean() / (s.std()**3 + 1e-12)),
        "kurtosis":  float(((s - s.mean())**4).mean() / (s.std()**4 + 1e-12)),
        "max_slope": np.abs(v).max(),
        "max_curv":  np.abs(a).max(),
        "zero_cross": int(np.sum(np.diff(np.sign(s)) != 0)),
        "decay_rate": _exp_decay_fit(t, s),   # fitted exponential slope
    }

Each descriptor targets a specific failure signature. A negative-energy event shows up as a sign change in internal energy and a spike in zero_cross. Hourglassing produces high-frequency oscillation, visible as elevated kurtosis and max_curv on the hourglass energy history. Mass scaling that is too aggressive inflates kinetic energy; the peak and decay_rate of the kinetic-energy curve expose it. Contact instability appears as a sharp max_slope in the contact-force history with a simultaneous drop in the time step.

Two practical rules apply. First, always normalize time by the total simulated duration so t_peak becomes a dimensionless fraction in [0, 1]; this makes features comparable across models with different end times. Second, compute descriptors on a windowed basis (e.g., 20 equal windows) when the failure is localized in time, then keep both the window-level and the whole-history aggregates.

Spatial field statistics

Time-series descriptors capture when something happened; spatial statistics capture where and how severely. For a nodal field f(x) over the mesh, we compute:

The volume-above-threshold feature deserves emphasis. A single element with 5 GPa stress is usually a numerical artifact; a contiguous region of 200 elements above 1 GPa is a real plastic hinge. By integrating over volume rather than taking the maximum, the feature becomes mesh-insensitive and physically interpretable. Choose the threshold τ from the material yield or failure criterion, not from the data, so the feature has consistent meaning across the dataset.

For vector fields such as displacement or velocity, first remove rigid-body motion by subtracting the mass-weighted mean translation and rotation, then compute the same statistics on the residual field. Without this step, a model that simply translates will show large displacement magnitudes that have nothing to do with deformation.

Solver convergence indicators

Implicit and explicit solvers emit a rich diagnostic stream that is often ignored. For implicit (Abaqus/Standard, Optistruct), the key quantities per increment are:

• Residual norm R and its ratio to the initial residual R₀.
• Displacement correction norm Δu.
• Number of equilibrium iterations n_iter.
• Time-step cutbacks n_cut and the cutback factor.
• Line-search failures and negative eigenvalue warnings.

For explicit (LS-DYNA, Abaqus/Explicit), the critical indicators are the minimum time step Δt_min, the mass-scaling percentage, the added-mass ratio, and the energy ratio (hourglass / internal, kinetic / internal).

These are aggregated into features such as: mean and max n_iter, fraction of increments requiring cutback, slope of log(R/R₀) versus iteration, and the ratio Δt_min / Δt_initial. A convergence failure typically shows a residual that plateaus instead of dropping, or a time step that collapses by more than two orders of magnitude. These patterns are exactly what a classifier can learn.

Putting the three families together

The figure below shows how the three feature families are extracted from a single simulation and concatenated into one diagnostic vector. Each family answers a different question: time-series descriptors answer when, spatial statistics answer where and how much, and convergence indicators answer how hard the solver struggled.

A worked example: detecting hourglassing

Consider a shell-element impact model that exhibits hourglassing in a small patch of elements. The diagnostic features that fire are:

• Time-series: hourglass energy history shows kurtosis > 8 and max_curv an order of magnitude above the baseline; the ratio of hourglass to internal energy peaks at 0.12 (threshold 0.05).
• Spatial: the volume-above-threshold for hourglass energy is small (a few elements) but the gradient magnitude is very high, indicating a localized oscillation rather than a smooth deformation.
• Convergence: the minimum time step drops by 40% in the affected region, and the mass-scaling percentage rises above 1%.

Individually, none of these is conclusive. Together, they form a signature that a Random Forest can separate from mesh distortion, contact penetration, and material instability with high accuracy. The key engineering insight is that the feature set must be designed around the failure modes you expect, not assembled generically. Start from the physics of each failure mode, ask what observable quantity changes, and derive the descriptor from that quantity.

Feature selection and validation

After extraction, the feature matrix typically has 50–200 columns. Not all are useful. Apply three filters in order:

1. Variance filter: drop features with near-zero variance across the dataset (they carry no discriminative information).
2. Correlation filter: drop one of each pair with |Pearson r| > 0.95 to reduce redundancy.
3. Model-based importance: train a gradient-boosted tree and keep the top-k features by gain, or use recursive feature elimination.

Validate the feature set with grouped cross-validation: split by simulation model, not by time step or by sample. Otherwise, features from the same run leak into both training and validation, and reported accuracy will be optimistic. A robust diagnostic model should achieve stable performance when entire models are held out.

The output of this stage is a clean, documented feature table — one row per simulation, one column per descriptor, plus metadata columns for model ID, material, impact velocity, and the ground-truth failure label. That table is the input to the model training and evaluation stages that follow.