4. Full Factorial DOE Matrix Generation for LS-DYNA

Why a Full Factorial Matrix Is the Reference Design

A full factorial design crosses every level of every factor with every level of every other factor. For k factors with L₁, L₂, …, L_k levels, the number of runs is the product of the level counts. This is the only design that lets you estimate every main effect and every interaction independently, which makes it the natural baseline before you move to space-filling (LHS) or orthogonal-array (Taguchi) designs. The practical task in this case is not the statistics — it is producing a clean, deterministic run table that a batch LS-DYNA driver can consume without ambiguity.

Assume the factor-level table from the previous step has already been fixed. A representative set for a crash-box study might be:

With 3 × 3 × 3 × 2 = 54 runs, the matrix is small enough to run overnight but large enough that hand-editing decks is error-prone. The generator below produces the run table and, critically, a stable run identifier that survives reordering.

Generating the Matrix Deterministically

The key design decision is ordering. Use a mixed-radix counter so that the last factor varies fastest; this keeps consecutive runs adjacent in the design space and makes partial results interpretable if the batch is interrupted. Each run gets an ID of the form R0001…R0054 and a hash of its factor tuple so that duplicate or missing combinations are detectable.

import csv, hashlib, itertools, json

# Factor table: name -> (keyword_field, [levels])
factors = {
    "t_shell":     ("*SECTION_SHELL/thickness", [1.2, 1.6, 2.0]),
    "rho_scale":   ("*MAT_024/RO_scale",        [0.9, 1.0, 1.1]),
    "yield_scale": ("*MAT_024/SIGY_scale",      [0.85, 1.0, 1.15]),
    "mu_fric":     ("*CONTACT/FS",              [0.10, 0.25]),
}

names  = list(factors.keys())
levels = [factors[n][1] for n in names]

rows = []
for idx, combo in enumerate(itertools.product(*levels), start=1):
    run_id = f"R{idx:04d}"
    key    = json.dumps(dict(zip(names, combo)), sort_keys=True)
    digest = hashlib.sha1(key.encode()).hexdigest()[:8]
    row = {"run_id": run_id, "combo_hash": digest}
    row.update(dict(zip(names, combo)))
    rows.append(row)

# Sanity checks
assert len(rows) == 54, "unexpected run count"
assert len({r["combo_hash"] for r in rows}) == 54, "duplicate combination"

with open("doe_full_factorial.csv", "w", newline="") as f:
    w = csv.DictWriter(f, fieldnames=["run_id", "combo_hash"] + names)
    w.writeheader()
    w.writerows(rows)

print(f"wrote {len(rows)} runs")

Running this produces doe_full_factorial.csv with 54 rows. The first three rows are:

run_id,combo_hash,t_shell,rho_scale,yield_scale,mu_fric
R0001,3f9a1c02,1.2,0.9,0.85,0.1
R0002,7b2e4d11,1.2,0.9,0.85,0.25
R0003,c41a8e77,1.2,0.9,1.0,0.1

Note that mu_fric varies fastest, then yield_scale, and so on. This ordering is intentional: if you later add a fifth factor, the existing 54 runs remain a valid prefix of the new 108-run design, so partial results stay comparable.

Reading the Design Structure

The SVG below shows how the four factors combine into the 54-run matrix. Each column is a factor, each row is a run, and the color bands make the level structure visible at a glance — you can immediately see that mu_fric alternates every row while t_shell changes only every 18 rows.

The block widths encode the level counts directly: t_shell holds each level for 18 consecutive runs, rho_scale for 6, yield_scale for 2, and mu_fric alternates every run. This is the mixed-radix structure that guarantees every combination appears exactly once.

Mapping Runs to LS-DYNA Decks

The run table is only useful if each row can be applied to a baseline deck. The mapping is a simple substitution driven by the keyword field recorded in the factor table. A minimal driver that produces one deck per run:

import csv, re, shutil, pathlib

baseline = pathlib.Path("baseline.k")
outdir   = pathlib.Path("runs"); outdir.mkdir(exist_ok=True)

# field -> (regex, replacement template)
# Each regex captures the numeric token after the keyword line.
patterns = {
    "t_shell":     (r"(\*SECTION_SHELL\s*\n[^\n]*\n[^\n]*\n\s*)([\d.]+)", r"\g<1>{v}"),
    "rho_scale":   (r"(\*MAT_024\s*\n[^\n]*\n[^\n]*\n[^\n]*\n\s*)([\d.Ee+-]+)", r"\g<1>{v}"),
    "yield_scale": (r"(\*MAT_024\s*\n[^\n]*\n[^\n]*\n[^\n]*\n[^\n]*\n[^\n]*\n\s*)([\d.Ee+-]+)", r"\g<1>{v}"),
    "mu_fric":     (r"(\*CONTACT_[A-Z_]+[^\n]*\n[^\n]*\n[^\n]*\n[^\n]*\n\s*)([\d.]+)", r"\g<1>{v}"),
}

with open("doe_full_factorial.csv") as f:
    for row in csv.DictReader(f):
        text = baseline.read_text()
        for name, (pat, repl) in patterns.items():
            text = re.sub(pat, repl.format(v=row[name]), text, count=1)
        (outdir / f"{row['run_id']}.k").write_text(text)

Verification Before Launching the Batch

Three checks catch almost every failure mode at this stage:

• Count and uniqueness. The combo_hash column must have exactly 54 distinct values. A duplicate hash means a level list was accidentally repeated; a short count means a factor was dropped.
• Deck diff. For a sample of runs (first, middle, last), diff against the baseline and confirm that exactly the expected number of lines changed and that each changed line contains the intended value.
• Keyword validity. Run ls-dyna -i runs/R0001.k -n cpu in check-only mode (or a short termination-time run) to confirm the deck parses. A malformed substitution usually surfaces as a keyword error before the solver starts.

With the matrix generated, validated, and mapped to decks, the batch is ready to launch. The next design decision — whether to keep the full 54-run factorial or reduce it — depends on solver budget and on which interactions the engineering question actually requires.