3. Design of experiments for crash simulations

Crashworthiness optimization almost always begins with a design of experiments (DOE). Before any surrogate model can be trained, you must decide which finite element (FE) simulations to run. A full-vehicle frontal impact or a thin-walled tube axial crush can take anywhere from 30 minutes to several hours per run, so the DOE is not a formality — it is the single most expensive decision in the whole workflow. This section covers the four DOE families most relevant to crash problems: full factorial, Latin hypercube sampling (LHS), Sobol sequences, and the space-filling criteria that let you judge whether a plan is actually good.

Why crash responses change the DOE problem

Crash responses are not smooth polynomials. Consider a front rail with a trigger bead: below a certain bead depth the tube folds in a stable, progressive mode; above it, the rail may split or fold asymmetrically, and peak force can jump by tens of percent. The response surface therefore contains kinks, plateaus, and even discontinuities. Three consequences follow directly:

• Edges matter. A DOE that clusters points in the interior will miss the mode transition. You need points near the bounds of each variable.
• One-factor-at-a-time is useless. Interactions dominate: gauge thickness and bead depth jointly determine the folding mode, so the DOE must vary variables simultaneously.
• Noise is real. Contact, hourglassing, and mesh sensitivity introduce run-to-run variation of a few percent even for a deterministic solver. Replicates at a few points are worth budgeting.

Full factorial: the reference design

A full factorial with k factors at L levels requires L^k runs. For a 3-factor, 3-level study that is 27 runs — entirely practical, and it gives you every main effect and every interaction, including the curvature that a 2-level design cannot see. For 6 factors at 3 levels it is 729 runs, which is no longer practical for crash. Full factorial is therefore best used as a screening reference on 2–4 factors, or as a small 2-level design (2^k) when you only need to rank variable importance.

A useful compromise is the face-centered central composite design (CCD): a 2^k factorial core, 2k axial points placed on the faces of the cube (α = 1), and center-point replicates. For k = 3 this is 8 + 6 + 3 = 17 runs, and it fits a full quadratic model with curvature — often the right choice when you have three or four thickness variables and want a defensible response surface.

Latin hypercube sampling

LHS partitions each variable's range into N equal-probability intervals and places exactly one sample in each interval, then randomly pairs the levels across variables. The result is a design with perfect one-dimensional stratification: project the points onto any single axis and they are evenly spread. For crash problems this is valuable because it guarantees coverage of the full thickness or bead-depth range without the exponential cost of a grid.

The weakness is that plain LHS says nothing about joint coverage. Two variables can be accidentally correlated, leaving a diagonal void in the design space. The standard fix is to generate many candidate LHS matrices and keep the one that maximizes a space-filling criterion — this is optimized LHS, and it is what most practitioners actually use.

import numpy as np
from scipy.stats import qmc

# 4 design variables: t1, t2, bead_depth, bead_radius
bounds = np.array([[1.2, 2.0],   # t1 [mm]
                   [1.2, 2.0],   # t2 [mm]
                   [3.0, 12.0],  # bead depth [mm]
                   [2.0, 8.0]])  # bead radius [mm]

sampler = qmc.LatinHypercube(d=4, seed=42)
u = sampler.random(n=40)                 # 40 samples in [0,1)^4
X = qmc.scale(u, bounds[:, 0], bounds[:, 1])

# Maximize minimum pairwise distance (maximin) over 200 candidates
best, best_d = None, -np.inf
for s in range(200):
    cand = qmc.scale(qmc.LatinHypercube(d=4, seed=s).random(40),
                     bounds[:, 0], bounds[:, 1])
    d = np.min(np.linalg.norm(cand[:, None] - cand[None, :], axis=2)
               + np.eye(40) * 1e9)
    if d > best_d:
        best, best_d = cand, d
print(f"maximin distance = {best_d:.3f} mm")

Note the units: the maximin distance is reported in the mixed unit space of millimetres, which is only meaningful if the variables are on comparable scales. When one variable spans 0.8 mm and another spans 9 mm, normalize to [0, 1] before computing distances, or the large-range variable will dominate the criterion.

Sobol sequences

Sobol sequences are low-discrepancy sequences: deterministic point sets constructed so that every prefix of the sequence is well distributed. Unlike LHS, they are extensible — you can add samples without discarding the existing ones, which is exactly what you want when a surrogate turns out to be inaccurate in one region and you need to enrich the DOE. Sobol points also fill space more uniformly than random LHS at low sample counts, and their discrepancy converges as O((log N)^d/N).

sampler = qmc.Sobol(d=4, scramble=True, seed=7)
X_sobol = qmc.scale(sampler.random_base2(m=6),   # 64 points
                    bounds[:, 0], bounds[:, 1])

Use powers of two for Sobol sample counts; the balance properties are guaranteed for those prefixes. Scrambling (random digital shifting) preserves the low-discrepancy structure while giving you independent replicates for error estimation.

Space-filling criteria you can actually compute

Three criteria cover most crash DOE decisions:

In practice, generate 100–500 candidate designs, compute all three metrics, and pick a design that is good on maximin and has |φ| < 0.05. A design with excellent maximin distance but a 0.4 correlation between t1 and t2 is worse than a slightly less uniform design with near-zero correlation, because the surrogate will be unable to separate the two effects.

Choosing a plan: a worked comparison

Suppose you have four variables (t1, t2, bead depth, bead radius) and a budget of 40 crash runs. Full factorial at 3 levels would need 81 runs — over budget. A face-centered CCD needs 25 runs for the quadratic model plus replicates, leaving room for 15 extra points. Optimized LHS at 40 points gives the best global coverage but no guaranteed corner coverage. A practical answer: run a 25-run CCD to fit a quadratic screening model, then add 15 Sobol points in the region where the quadratic residual is largest. This is a sequential DOE, and it is the standard industrial pattern.

Practical checklist before you submit the DOE

1. Normalize all variables to [0, 1] before computing any distance or discrepancy metric.
2. Verify |φ| < 0.05 for every variable pair; regenerate if not.
3. Confirm that at least 4–8 points sit on or within 2% of the design-space corners.
4. Add 2–3 replicate runs at the center point to estimate solver noise; if the coefficient of variation exceeds 5%, investigate contact or hourglass control before trusting the surrogate.
5. Record the DOE matrix with variable names, units, and bounds in a CSV alongside the FE results — the surrogate is only as reproducible as its input table.

With a defensible DOE in hand — typically 30–60 optimized LHS or Sobol points plus corner and replicate runs — you have the training data needed for the surrogate models covered in the next section. The quality of that data, not the choice of regression algorithm, will set the ceiling on your optimization accuracy.