3. Python environment for CAE

A parametric crash model is, underneath the geometry, a numerical pipeline. Section coordinates, thickness fields, material curves, and barrier positions all become arrays; the solver consumes those arrays and returns time histories that must be reduced, filtered, and compared. If the Python environment is assembled carelessly, the parametric layer becomes fragile long before the first LS-DYNA or Radioss deck is written. This section establishes a reproducible environment for parametric CAE development, covering the four pillars: NumPy for array geometry and vectorized math, SciPy for interpolation, spatial queries, and signal processing, Pandas for parameter tables and result aggregation, and a mesh library for reading and writing solver meshes.

Choosing the interpreter and isolating the environment

CAE toolchains are sensitive to interpreter version and binary compatibility. A practical baseline is Python 3.11 or 3.12, 64-bit, managed through conda (or mamba) because mesh and solver-adjacent packages frequently ship compiled extensions that are easier to resolve through conda-forge than through pure pip wheels. Create a dedicated environment rather than installing into the base interpreter:

conda create -n cae-param python=3.11 -c conda-forge
conda activate cae-param
conda install -c conda-forge numpy scipy pandas matplotlib pyvista meshio h5py
pip install lxml pyyaml pytest

Pin the environment so a colleague or a CI runner reproduces it exactly. Export after the first successful import test:

conda env export --no-builds > environment.yml
python -c "import numpy, scipy, pandas, meshio; print(numpy.__version__, scipy.__version__, pandas.__version__, meshio.__version__)"

If the printed versions differ from the pinned file, the environment is not reproducible and any parametric result is suspect. Treat the version string as part of the model definition, not as incidental metadata.

NumPy: the coordinate and parameter substrate

In a parametric crash model, the node table is the central data structure. A node array of shape (n_nodes, 3) with float64 dtype supports vectorized translation, mirroring, and scaling without Python loops. The following example builds a symmetric rail section and mirrors it about the vehicle centerline, which is the Y = 0 plane in the standard SAE vehicle coordinate system (X forward, Y left, Z up):

import numpy as np

# Half-section profile in the Y-Z plane, X = station coordinate
y = np.array([0.0, 40.0, 40.0, 0.0])      # mm
z = np.array([-30.0, -30.0, 30.0, 30.0])  # mm
x = np.full_like(y, 250.0)                # station X = 250 mm

half = np.column_stack((x, y, z))
mirror = half.copy()
mirror[:, 1] *= -1.0                      # reflect across Y = 0
section = np.vstack((half, mirror[::-1]))

print(section.shape, section.dtype)       # (8, 3) float64

Two conventions matter here and must stay consistent through the whole project. First, units: keep everything in millimetres and tonnes-millimetre-seconds, matching the solver deck, and never mix metres into a mesh array. Second, handedness: mirroring across Y = 0 flips the winding of shell elements, so the element connectivity must be reordered when the mirrored section is meshed. Recording these conventions in a module-level docstring prevents the classic failure where a mirrored rail has inverted normals and the solver reports negative volumes.

SciPy: interpolation, spatial queries, and filtering

Parametric studies rarely reuse the same mesh; they reuse the same topology with different dimensions. SciPy provides the numerical glue. scipy.interpolate.interp1d or PchipInterpolator maps a design variable such as rail gauge to a section width; scipy.spatial.cKDTree finds the nearest node to a barrier reference point; scipy.signal.butter and sosfiltfilt apply the SAE J211 CFC 60 filter to acceleration histories.

from scipy.interpolate import PchipInterpolator
from scipy.signal import butter, sosfiltfilt

# Design table: gauge (mm) -> rail width (mm)
gauge = np.array([600.0, 700.0, 800.0, 900.0])
width = np.array([120.0, 140.0, 165.0, 190.0])
width_of = PchipInterpolator(gauge, width)
print(width_of(750.0))                    # ~152.5 mm, monotone, no overshoot

# SAE J211 CFC 60 low-pass filter at 10 kHz sampling
fs, fc = 10000.0, 60.0
sos = butter(2, fc / (fs / 2.0), btype="low", output="sos")
a_filtered = sosfiltfilt(sos, a_raw)

Use PchipInterpolator rather than a cubic spline for design tables: it is shape-preserving, so a monotonically increasing width table cannot produce a non-physical dip between stations. For filtering, always apply zero-phase filtering (sosfiltfilt) and document the CFC class, because a CFC 60 and a CFC 180 result are not comparable.

Pandas: parameter sweeps and result aggregation

A parametric crash study is a table problem. Each row is a design point; each column is an input variable or a response metric. Pandas makes the sweep auditable and the post-processing repeatable.

import pandas as pd

designs = pd.DataFrame({
    "run_id":   ["D01", "D02", "D03"],
    "gauge":    [700.0, 750.0, 800.0],
    "thickness":[1.6, 1.8, 2.0],
    "material": ["DP590", "DP590", "DP780"],
})
designs["rail_width"] = width_of(designs["gauge"].to_numpy())

results = pd.DataFrame({
    "run_id": ["D01", "D02", "D03"],
    "peak_force_kN": [182.4, 196.1, 214.7],
    "intrusion_mm":  [312.0, 288.5, 261.3],
    "mass_kg":       [41.2, 43.8, 46.5],
})
summary = designs.merge(results, on="run_id")
summary["specific_energy_kJ_per_kg"] = (
    summary["peak_force_kN"] * summary["intrusion_mm"] / 1000.0 / summary["mass_kg"]
)
print(summary.sort_values("intrusion_mm"))

The merge key is the contract between the design generator and the post-processor. If a run fails to produce a result row, the merge exposes it as a NaN rather than silently dropping the design point — a missing run must be visible, not invisible. Store the summary as Parquet (summary.to_parquet("sweep.parquet")) so dtypes survive round-trips and large sweeps load quickly.

Mesh libraries: reading and writing solver decks

Parametric geometry must eventually become a mesh the solver accepts. meshio provides a format-agnostic container for points and cells and can read and write Nastran, Abaqus, VTK, and other formats; pyvista builds on VTK for inspection and quality checks. A typical workflow reads a baseline deck, replaces the node coordinates with the parametric section, and writes a new deck:

import meshio

mesh = meshio.read("baseline_rail.nas")
points = mesh.points.copy()
points[:, 1:3] = section[:, 1:3]          # apply parametric Y-Z section
mesh.points = points
meshio.write("rail_D02.nas", mesh, file_format="nastran")

Before writing, verify mesh quality: check for duplicate nodes within a tolerance, zero-area shell elements, and free edges. A minimal check using NumPy is sufficient for a first pass:

from scipy.spatial import cKDTree

tree = cKDTree(mesh.points)
pairs = tree.query_pairs(r=1e-6)
assert not pairs, f"duplicate nodes detected: {len(pairs)} pairs"

If the assertion fires, the parametric section has produced coincident nodes at the mirror plane — a common symptom when the half-section does not start exactly at Y = 0. The fix is to snap the seam nodes to zero before mirroring, not to increase the tolerance.

The diagram above shows how the four libraries divide responsibility. Pandas owns the design table and the result merge; NumPy owns the coordinate arrays; SciPy supplies interpolation, spatial queries, and filtering; the mesh library converts arrays into solver-readable decks. The red dashed path is the feedback loop that closes the parametric study: solver responses return to the design table, where they become inputs to the next iteration. The quality gate sits between mesh generation and the solver, and the conventions box records the coordinate and unit assumptions that every stage must honour.

Reproducibility as an engineering requirement

Parametric crash development is only meaningful if a design point can be regenerated exactly. Three practices make that possible. First, pin library versions in environment.yml and record them alongside each sweep. Second, keep the design table, the mesh-generation script, and the post-processing script under version control, with the run identifier as the join key. Third, validate the environment with the sanity script before every sweep, so that a broken import or a mismatched NumPy ABI is caught before hours of solver time are spent. With this foundation in place, the next stage — building the parametric geometry and mesh generators themselves — can proceed on a stable numerical base.