3. Design of experiments for LS-DYNA

An LS-DYNA explicit run is expensive: a single 1–2 million element crash or impact model may take 30 minutes to several hours on 32–64 cores, and a parametric study that touches eight or ten inputs can easily consume thousands of core-hours. The design of experiments (DOE) stage is therefore not a formality — it is the decision that determines how much information you extract per simulation. This section covers the three DOE families most often used with LS-DYNA, how to generate a run matrix that respects physical and solver constraints, and how to verify the matrix before committing compute.

Choosing the design family from the response surface you expect

Full factorial evaluates every combination of levels. With k factors at L levels it produces L^k runs: three factors at three levels is 27 runs, five factors at three levels is 243. Its strength is that it resolves all main effects and all interaction terms up to the full order, and it is the only design that lets you compute a complete ANOVA with no aliasing. Its weakness is exponential growth. In LS-DYNA practice, full factorial is appropriate for two to four factors at two or three levels — for example, sheet thickness, yield stress, and friction coefficient in a forming study — or as a screening design at two levels before a refinement stage.

Latin hypercube sampling (LHS) divides each factor's range into N equal-probability strata and places exactly one sample per stratum per factor, then pairs the strata across factors by random permutation. With N runs you get N samples regardless of the number of factors, so LHS scales to ten or twenty inputs where full factorial is impossible. The catch is that a naive LHS can produce poor projections: two factors may end up correlated, or samples may cluster along a diagonal. This is why LHS is almost always paired with an optimization criterion.

Space-filling designs — maximin LHS, orthogonal-array LHS, Halton/Hammersley sequences, and Sobol sequences — explicitly optimize the spread of points in the unit hypercube. Maximin LHS maximizes the minimum pairwise distance; orthogonal-array LHS minimizes the discrepancy between the empirical and uniform distributions in every two-dimensional projection. For LS-DYNA surrogate modeling (kriging, radial basis functions, polynomial chaos), space-filling designs are the default because surrogate accuracy depends on covering the interior of the domain, not on placing points at corners.

Sampling constraints: where the run matrix meets the solver

A raw DOE matrix is a set of points in a hypercube. An LS-DYNA run matrix is a set of *KEYWORD decks that must actually execute. Three classes of constraint intervene between them.

Physical and geometric constraints. Not every combination of inputs is manufacturable or geometrically valid. If thickness t and bend radius R must satisfy R ≥ 2t, a full factorial over independent ranges will generate infeasible decks. The standard remedy is constrained sampling: generate a large candidate pool (e.g. 10,000 LHS points), reject candidates violating the constraint, then select the N runs from the feasible pool using the same space-filling criterion. This preserves stratification approximately while guaranteeing feasibility.

Solver stability constraints. Explicit LS-DYNA time step is governed by the Courant condition Δt = L_c/c, where L_c is the characteristic element length and c the sound speed. If a DOE factor is mesh size, small elements drive the time step down and can push a run from 40 minutes to 6 hours. Similarly, extreme values of strain-rate parameters in *MAT_PLASTIC_KINEMATIC or *MAT_024 can trigger negative volume errors or excessive hourglassing. Practical mitigation: cap the factor range at values you have already validated in single runs, and record the estimated critical time step for each run before submission.

Discrete and categorical factors. Material grade, element formulation (ELFORM = 1, 2, −1, −2), and contact type are categorical. LHS assumes continuous factors; for mixed designs, stratify the continuous factors within each categorical level, or use a sliced LHS so that each categorical slice is itself a valid LHS. A common pattern for a 3-material × 4-continuous-factor study is 3 slices of 20 runs each, giving 60 runs total.

Generating the run matrix

The workflow below uses Python with scipy.stats.qmc and pandas; the same logic applies in LS-OPT, optiSLang, or Dakota. The example studies a simplified frontal impact model with four continuous factors and one categorical factor.

import numpy as np, pandas as pd
from scipy.stats import qmc

# Factor definitions: (name, low, high, unit)
factors = [
    ("thickness_mm",   1.2, 2.4, "mm"),   # rail wall thickness
    ("yield_MPa",    280.0, 420.0, "MPa"),# *MAT_024 SIGY
    ("friction",       0.05, 0.25, "-"),  # *CONTACT friction coeff
    ("strainrate_c",  20.0, 80.0, "1/s"), # Cowper-Symonds C
]
materials = ["DP590", "DP780", "HSLA420"]   # categorical slice

N_PER_SLICE = 20
sampler = qmc.LatinHypercube(d=len(factors), seed=42,
                             optimization="random-cd")  # maximin-ish
rows = []
for mat in materials:
    u = sampler.random(N_PER_SLICE)          # shape (20, 4), in [0,1)
    for i, (name, lo, hi, unit) in enumerate(factors):
        col = lo + u[:, i] * (hi - lo)
        if name == "thickness_mm":
            col = np.round(col, 2)           # manufacturable resolution
        rows.append(pd.DataFrame({
            "run_id": [f"{mat}_{j:02d}" for j in range(N_PER_SLICE)],
            "material": mat, name: col}))
matrix = pd.concat(rows, axis=1)
matrix.to_csv("doe_run_matrix.csv", index=False)

Two details matter here. First, optimization="random-cd" applies a Centered Discrepancy optimization to the LHS, which is a cheap proxy for maximin and typically improves projection uniformity by 20–40% over naive LHS at negligible cost. Second, rounding thickness to 0.01 mm is a real constraint — sheet metal comes in discrete gauges — and rounding after sampling slightly degrades stratification, so verify the rounded matrix rather than the raw one.

Verifying the matrix before submission

Before launching 60 LS-DYNA jobs, run three checks on the CSV. Pairwise correlation: compute the Pearson correlation of every factor pair; absolute values above 0.3 in a 20-run slice indicate a poor permutation and warrant a new seed. Minimum distance: compute the minimum Euclidean distance in normalized [0,1] coordinates; for 20 points in 4-D, values below 0.4 suggest clustering. Constraint audit: evaluate every physical constraint row-wise and confirm zero violations. A short verification script:

from scipy.spatial.distance import pdist
X = matrix[[f[0] for f in factors]].to_numpy()
Xn = (X - X.min(0)) / (X.max(0) - X.min(0))
print("min pairwise dist:", pdist(Xn).min().round(3))
print("max |corr|:", np.abs(np.corrcoef(Xn.T) - np.eye(4)).max().round(3))

For the example above with seed 42, a typical outcome is a minimum distance near 0.55 and maximum off-diagonal correlation near 0.18 — acceptable for surrogate fitting. If the correlation check returns 0.6, do not proceed; regenerate with a different seed or increase the candidate pool.

The output of this stage is a frozen, verified run matrix: one row per LS-DYNA job, one column per factor, plus a categorical slice label and a run identifier. That matrix is the contract between DOE and the simulation campaign, and every subsequent stage — surrogate construction, adaptive infill, and optimization — depends on its integrity.