3. Surrogate modeling for crashworthiness
Full-frontal and offset crash simulations routinely consume 8–24 wall-clock hours per design point on 32–128 cores. A design-of-experiments (DOE) campaign with 200 samples therefore costs weeks of cluster time, which makes direct optimization of intrusion, acceleration, and force metrics impractical. Surrogate modeling replaces the finite element (FE) solver with a cheap analytical approximation trained on a limited set of high-fidelity runs. This chapter examines four families — response surface methods (RSM), Gaussian process regression (GPR), polynomial chaos expansion (PCE), and deep learning (DL) surrogates — and how each maps to the crashworthiness responses that matter most: intrusion (toe pan, dash, B-pillar), acceleration (peak and 3 ms clip of the occupant compartment), and force (section force, peak force, and force–displacement curve integrals).
3.1 What the Surrogate Must Reproduce
Before choosing a model, define the response vector. A typical frontal crashworthiness study uses:
- Intrusion metrics: maximum toe-pan intrusion
d_toe(mm), dash panel intrusiond_dash, B-pillar intrusion at beltlined_bpillar. - Acceleration metrics: peak resultant acceleration
a_peak(g), 3 ms clipa_3ms, and occupant compartment deceleration pulse. - Force metrics: peak section force
F_peak(kN), mean crush forceF_mean, and energy absorbedE_abs(kJ).
These responses are non-linear, non-monotonic, and often discontinuous in the design variables (e.g., gauge thickness, material grade, rib layout, trigger geometry). A surrogate that fits d_toe smoothly may fail badly on a_peak, because acceleration is dominated by local buckling events that produce sharp peaks. This is the central engineering constraint: one surrogate per response family, not one surrogate for all.
3.2 Response Surface Methods (RSM)
RSM fits a low-order polynomial — usually quadratic — to the DOE data. For k variables, the full quadratic model has 1 + k + k(k+1)/2 coefficients. With k = 5 (e.g., three thicknesses, one material grade, one trigger radius), that is 21 coefficients, requiring at least 21–30 samples for a stable fit.
# Quadratic RSM for toe-pan intrusion (example, synthetic data)
# Variables: t1, t2, t3 (thickness mm), mat (grade index), r_trig (mm)
# Model: d_toe = b0 + sum(bi*xi) + sum(bii*xi^2) + sum(bij*xi*xj)
import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
X = np.array([[1.2, 1.5, 1.8, 2, 8],
[1.4, 1.5, 1.8, 2, 8],
[1.2, 1.7, 1.8, 2, 8],
[1.2, 1.5, 2.0, 2, 8],
[1.2, 1.5, 1.8, 3, 8],
[1.2, 1.5, 1.8, 2, 12]])
y = np.array([42.1, 39.8, 40.5, 38.9, 36.2, 44.7]) # d_toe mm
poly = PolynomialFeatures(degree=2, include_bias=True)
Xp = poly.fit_transform(X)
model = LinearRegression().fit(Xp, y)
print("R2 =", model.score(Xp, y))
RSM is fast, interpretable, and adequate when the response is smooth and the design space is small. Its weakness is clear: it cannot capture buckling-induced discontinuities, and extrapolation beyond the DOE box is unreliable. Use RSM for early screening and for responses such as F_mean that vary smoothly with gauge thickness.
3.3 Gaussian Process Regression (GPR)
GPR places a prior over functions and updates it with observed data. The prediction at a new point x* is:
mu(x*) = k(x*, X) [K(X, X) + sigma_n^2 I]^-1 y
var(x*) = k(x*, x*) - k(x*, X) [K(X, X) + sigma_n^2 I]^-1 k(X, x*)
The kernel k encodes smoothness assumptions. For crash responses, a Matérn 5/2 kernel is a common choice because it allows rougher functions than the squared-exponential kernel:
k(r) = (1 + sqrt(5) r / l + 5 r^2 / (3 l^2)) exp(-sqrt(5) r / l)
GPR returns both a mean prediction and a variance, which is invaluable for adaptive DOE: the next FE run is placed where the predicted variance or expected improvement is highest. For a_peak, GPR with Matérn 5/2 and a noise term typically achieves R² ≈ 0.85–0.95 on 60–100 samples, whereas quadratic RSM often stalls near R² ≈ 0.6–0.75. The cost is O(n³) training and O(n²) prediction, which becomes noticeable beyond ~500 samples.
a_peak, consider a classification layer that first predicts whether a buckling mode is triggered, then regresses within each mode.3.4 Polynomial Chaos Expansion (PCE)
PCE expands the response in orthogonal polynomials of the input random variables:
Y(x) = sum_{alpha in A} c_alpha * Psi_alpha(xi)
where Psi_alpha are multivariate orthogonal polynomials (Legendre for uniform inputs, Hermite for Gaussian inputs) and c_alpha are coefficients computed by projection or regression. PCE is particularly strong when inputs have known probability distributions and the response is smooth in the stochastic sense. It yields analytical moments (mean, variance, Sobol indices) directly from the coefficients, which makes it ideal for uncertainty quantification of intrusion and force metrics under manufacturing scatter (e.g., thickness tolerance ±0.1 mm, material property scatter ±5%).
For crashworthiness, PCE with total degree 3–4 and sparse regression (e.g., LARS) can handle 5–10 variables with 50–150 samples. The limitation is the same as RSM: Gibbs-type oscillations near discontinuities. A hybrid approach — PCE for smooth global trends plus GPR for local residuals — often outperforms either alone.
3.5 Deep Learning Surrogates
When the DOE exceeds several hundred samples and the response is highly non-linear or high-dimensional (e.g., full force–displacement curves, 3D intrusion fields), neural networks become competitive. A typical architecture for scalar crash metrics is a fully connected network with 3–5 hidden layers, 64–256 units per layer, ReLU or Swish activation, and dropout for regularization. For curve prediction, an encoder–decoder or a temporal convolutional network (TCN) maps design variables to the full F(s) curve.
# Example: PyTorch MLP for a_peak prediction
import torch, torch.nn as nn
model = nn.Sequential(
nn.Linear(5, 128), nn.SiLU(),
nn.Linear(128, 128), nn.SiLU(),
nn.Linear(128, 64), nn.SiLU(),
nn.Linear(64, 1)
)
# Input: [t1, t2, t3, mat, r_trig]
# Output: a_peak (g)
DL surrogates require careful normalization, learning-rate scheduling, and early stopping. They do not provide native uncertainty, so ensemble methods (5–10 models) or Monte Carlo dropout are used to estimate prediction variance. In practice, a well-tuned ensemble reaches R² ≈ 0.90–0.97 for d_toe and F_peak with 300–500 samples, but can overfit badly below 150 samples.
3.6 Choosing and Validating the Surrogate
The selection depends on sample budget, response smoothness, and whether uncertainty is needed. The table below summarizes the trade-offs for the three crash metrics.
| Method | Typical samples | Best-suited response | Uncertainty | Key risk |
|---|---|---|---|---|
| Quadratic RSM | 20–40 | F_mean, E_abs |
No | Poor at discontinuities |
| GPR (Matérn 5/2) | 60–150 | d_toe, d_dash |
Yes | O(n³) scaling |
| PCE (sparse, degree 3–4) | 50–150 | d_toe, F_peak under scatter |
Analytical | Oscillations near jumps |
| DL ensemble | 300–500 | a_peak, full F(s) curve |
Approximate | Overfitting, data hunger |
Validation must use a held-out set that was never used for training or hyperparameter tuning. Report R², root-mean-square error (RMSE), and — critically — the maximum absolute error at the design point that matters most (e.g., the optimum). A surrogate with R² = 0.95 but 8 mm error on d_toe at the candidate optimum is not acceptable if the target is 50 mm with a 5 mm margin. Always run a confirmation FE simulation at the final optimum.
3.7 Workflow Summary
The practical sequence is: (1) define the response vector and target tolerances; (2) run a space-filling DOE (Latin hypercube or Sobol sequence) with 30–60 initial samples; (3) fit a fast RSM to screen variables; (4) fit GPR or PCE on the reduced variable set; (5) use GPR variance or PCE Sobol indices to place 20–40 adaptive samples; (6) if a_peak or curve prediction is required, train a DL ensemble on the accumulated 300+ samples; (7) validate on a held-out set and confirm the optimum with one FE run. The surrogate is not a replacement for the FE model — it is a filter that decides where the expensive FE runs should go.
The map above shows the intended pairing: intrusion and force metrics are usually well served by GPR or PCE, acceleration peaks demand a DL ensemble or a mode-aware hybrid, and full force–displacement curves require a sequence model. The validation gate is not optional — it is the only safeguard against a surrogate that looks accurate on average but fails at the design point you actually care about.