4. Surrogate modeling for crashworthiness and impact simulations

A single LS-DYNA crashworthiness run — a full vehicle frontal impact with a detailed finite element mesh — can take anywhere from several hours to more than a day on a compute cluster. When an adaptive DOE needs to explore a design space with ten or twenty variables, running the solver at every candidate point is simply not feasible. Surrogate modeling (also called metamodeling) solves this by fitting a cheap mathematical approximation to a limited set of expensive LS-DYNA results, then using that approximation to predict responses at untested design points. This chapter covers the three surrogate types most commonly paired with LS-DYNA crash data — Kriging, polynomial response surface (PRS), and radial basis functions (RBF) — and how to build, validate, and select among them using real simulation output.

What the Surrogate Must Capture

Crashworthiness responses are rarely smooth. Intrusion, peak acceleration, and energy absorption often exhibit nonlinearity, threshold effects (a rail buckling mode switching from axial to bending), and occasional discontinuities when contact or failure criteria activate. A surrogate is only as good as its ability to represent these features within the sampled region. The three methods differ fundamentally in how they do this:

• Polynomial response surface (PRS): fits a global low-order polynomial (typically quadratic) by least squares. It is smooth, fast, and interpretable, but struggles with highly nonlinear or multimodal crash responses unless the design space is small and well-behaved.
• Kriging (Gaussian process regression): interpolates the sample points exactly and models spatial correlation through a fitted variogram or covariance function. It excels at nonlinear, smooth-but-wiggly responses and provides a prediction variance that quantifies uncertainty — valuable for adaptive sampling.
• Radial basis functions (RBF): interpolates using a weighted sum of kernel functions centered at each sample point. It handles scattered, high-dimensional data well and is robust for noisy or irregular designs, though it lacks Kriging's built-in uncertainty estimate.

Preparing LS-DYNA Data for Fitting

The raw material for any surrogate is a design matrix X (one row per run, one column per design variable) and a response vector y (e.g., maximum intrusion in mm, peak deceleration in g, or mass in kg). Two preparation steps matter enormously:

1. Normalize the inputs. Design variables often span different units and magnitudes (a thickness in mm, an angle in degrees, a material yield stress in MPa). Scale each to [0, 1] or standardize to zero mean and unit variance before fitting; otherwise the correlation model in Kriging and the distance metric in RBF will be dominated by the largest-magnitude variable.
2. Screen the response for run failures. LS-DYNA runs can terminate with negative volume, excessive distortion, or contact errors. A failed run that still writes a plausible-looking peak value will poison the fit. Flag runs by checking the termination status and energy error; exclude or re-run them before fitting.

A typical extraction step pulls the response from LS-DYNA's glstat and nodout files. For example, peak acceleration at a seat-track node can be read from the nodout time history, while total internal energy comes from glstat. The following Python snippet assembles the training set from a directory of completed runs:

import numpy as np, glob, os

def load_run(run_dir):
    # design variables stored in a per-run params file
    params = np.loadtxt(os.path.join(run_dir, "params.txt"))
    # peak acceleration from nodout (column 2), intrusion from a node set
    nodout = np.loadtxt(os.path.join(run_dir, "nodout_hist.txt"))
    peak_acc = nodout[:, 1].max()
    intrusion = nodout[:, 2].max()
    return params, np.array([peak_acc, intrusion])

X, Y = [], []
for d in sorted(glob.glob("runs/run_*")):
    p, y = load_run(d)
    X.append(p); Y.append(y)
X = np.array(X); Y = np.array(Y)
# min-max normalize inputs to [0,1]
X = (X - X.min(0)) / (X.max(0) - X.min(0) + 1e-12)

Fitting the Three Surrogates

With normalized data in hand, each surrogate is fit and validated. A practical workflow uses a held-out test set (say 20% of runs) plus k-fold cross-validation on the remainder, reporting both R² and root-mean-square error (RMSE) in the response's native units so the error is interpretable to the crash engineer.

Polynomial response surface. A full quadratic model in n variables has 1 + n + n(n+1)/2 coefficients. For five variables that is 21 coefficients, requiring at least ~30 runs for a stable least-squares fit. Fit with ordinary least squares and inspect the residuals; a strong pattern in the residuals signals that the quadratic form is too stiff for the response.

Kriging. Fit a constant or linear trend plus a Gaussian correlation function with per-dimension length scales. The length scales are estimated by maximum likelihood. The key output is not just the prediction but the mean squared error (MSE) at each candidate point — this is what drives adaptive infill sampling in later chapters. A well-fit Kriging model on a smooth intrusion response typically reaches R² above 0.95 with 40–60 runs; a response with a buckling-mode switch may need more samples or a non-stationary variant.

Radial basis functions. Choose a kernel (multiquadric, inverse multiquadric, or thin-plate spline) and solve the linear system for the weights. RBF interpolation is exact at training points, so validation must rely on held-out data. It is particularly effective when the DOE is space-filling (Latin hypercube) and the response is irregular.

Comparing and Selecting

The table below summarizes the trade-offs as they apply to LS-DYNA crash data. These are qualitative guidance, not fixed rules — always validate on your own response.

Worked Example: Frontal Rail Intrusion

Consider a simplified frontal impact model with three design variables: front rail thickness t (1.2–2.4 mm), bumper beam gauge g (0.8–1.6 mm), and crash-box length L (120–200 mm). The response is maximum firewall intrusion. A 30-run Latin hypercube DOE is executed in LS-DYNA; 24 runs form the training set and 6 are held out. After normalization, a quadratic PRS yields cross-validated RMSE of 8.1 mm, Kriging 3.4 mm, and RBF 4.7 mm. Kriging is selected. Its prediction at an untested point (t = 1.9 mm, g = 1.1 mm, L = 165 mm) is 42 mm intrusion with an MSE of 2.1 mm², which the engineer can compare against the 50 mm target. The Kriging MSE also flags regions where additional LS-DYNA runs would most reduce uncertainty — the bridge to adaptive sampling.

Practical Guardrails

Three habits separate a usable surrogate from a misleading one. First, always report error in engineering units, not just R², and compare it against the tolerance that matters for the decision. Second, keep a small independent validation set that was never used for fitting or model selection — this catches optimistic bias from tuning length scales or kernel parameters on the same data. Third, version the training set alongside the surrogate: when a design variable range changes or a new LS-DYNA run is added, the model must be refit, and a stale surrogate silently applied to a shifted design space is one of the most common causes of bad optimization results. With these practices in place, the surrogate becomes a reliable, fast stand-in for the solver — the foundation on which the adaptive sampling and optimization stages of this course are built.