diff --git a/projects/hydrotech-beam/studies/01_doe_landscape/README.md b/projects/hydrotech-beam/studies/01_doe_landscape/README.md
index 662a6e13..307926eb 100644
--- a/projects/hydrotech-beam/studies/01_doe_landscape/README.md
+++ b/projects/hydrotech-beam/studies/01_doe_landscape/README.md
@@ -4,7 +4,7 @@
 
 ## Purpose
 
-Map the design space of the Hydrotech sandwich I-beam to identify feasible regions, characterize variable sensitivities, and converge on a minimum-mass design that satisfies displacement and stress constraints.
+Map the design space of the Hydrotech sandwich I-beam to identify feasible regions, characterize variable sensitivities, and prepare data for Phase 2 (TPE optimization). This is Phase 1 of the two-phase strategy (DEC-HB-002).
 
 ## Quick Facts
 
@@ -13,59 +13,175 @@ Map the design space of the Hydrotech sandwich I-beam to identify feasible regio
 | **Objective** | Minimize mass (kg) |
 | **Constraints** | Tip displacement ≤ 10 mm, Von Mises stress ≤ 130 MPa |
 | **Design variables** | 4 (3 continuous + 1 integer) |
-| **Algorithm** | Phase 1: LHS DoE (50 trials) → Phase 2: TPE (60-100 trials) |
-| **Total budget** | 114–156 NX evaluations |
-| **Estimated compute** | ~4–5 hours |
-| **Status** | DRAFT — Awaiting pre-flight checks and review |
+| **Algorithm** | Phase 1: LHS DoE (50 trials + 1 baseline) |
+| **Total budget** | 51 evaluations |
+| **Constraint handling** | Deb's feasibility rules (Optuna constraint interface) |
+| **Status** | Code complete — ready for execution |
 
-## Design Variables
+## Design Variables (confirmed via NX binary introspection)
 
-| ID | Name | Range | Type | Baseline |
-|----|------|-------|------|----------|
-| DV1 | `beam_half_core_thickness` | 10–40 mm | Continuous | 20 mm |
-| DV2 | `beam_face_thickness` | 10–40 mm | Continuous | 20 mm |
-| DV3 | `holes_diameter` | 150–450 mm | Continuous | 300 mm |
-| DV4 | `hole_count` | 5–15 | Integer | 10 |
+| ID | NX Expression | Range | Type | Baseline | Unit |
+|----|--------------|-------|------|----------|------|
+| DV1 | `beam_half_core_thickness` | 10–40 | Continuous | 25.162 | mm |
+| DV2 | `beam_face_thickness` | 10–40 | Continuous | 21.504 | mm |
+| DV3 | `holes_diameter` | 150–450 | Continuous | 300 | mm |
+| DV4 | `hole_count` | 5–15 | Integer | 10 | — |
 
-## Baseline Performance
+## Usage
 
-| Metric | Value | Constraint | Status |
-|--------|-------|------------|--------|
-| Mass | ~974 kg | minimize | Overbuilt |
-| Tip displacement | ~22 mm | ≤ 10 mm | ❌ FAILS |
-| VM stress | unknown | ≤ 130 MPa | ⚠️ TBD |
+### Prerequisites
 
-## Key Decisions
+```bash
+pip install -r requirements.txt
+```
 
-- **Single-objective** formulation (DEC-HB-001)
-- **Two-phase** algorithm: LHS → TPE (DEC-HB-002)
-- **True integer** handling for hole_count (DEC-HB-003)
-- **Deb's feasibility rules** for constraint handling (infeasible baseline)
-- **Baseline enqueued** as Trial 0 (LAC lesson)
+Requires: Python 3.10+, optuna, scipy, numpy, pandas.
+
+### Development Run (stub solver)
+
+```bash
+# From the study directory:
+cd projects/hydrotech-beam/studies/01_doe_landscape/
+
+# Run with synthetic results (for pipeline testing):
+python run_doe.py --backend stub
+
+# With verbose logging:
+python run_doe.py --backend stub -v
+
+# Custom study name:
+python run_doe.py --backend stub --study-name my_test_run
+```
+
+### Production Run (NXOpen on Windows)
+
+```bash
+# On dalidou (Windows node with NX):
+python run_doe.py --backend nxopen --model-dir "C:/path/to/syncthing/Beam"
+```
+
+### Resume Interrupted Study
+
+```bash
+python run_doe.py --backend stub --resume --study-name hydrotech_beam_doe_phase1
+```
+
+### CLI Options
+
+| Flag | Default | Description |
+|------|---------|-------------|
+| `--backend` | `stub` | `stub` (testing) or `nxopen` (real NX) |
+| `--model-dir` | — | Path to NX model files (required for `nxopen`) |
+| `--study-name` | `hydrotech_beam_doe_phase1` | Optuna study name |
+| `--n-samples` | 50 | Number of LHS sample points |
+| `--seed` | 42 | Random seed for reproducibility |
+| `--results-dir` | `results/` | Output directory |
+| `--resume` | false | Resume existing study |
+| `-v` | false | Verbose (DEBUG) logging |
 
 ## Files
 
 | File | Description |
 |------|-------------|
-| `OPTIMIZATION_STRATEGY.md` | Full strategy document (problem formulation, algorithm selection, risk mitigation) |
-| `atomizer_spec_draft.json` | AtomizerSpec configuration skeleton (DRAFT — open items must be resolved) |
-| `README.md` | This file |
+| `run_doe.py` | Main entry point — orchestrates the DoE study |
+| `sampling.py` | LHS generation with stratified integer sampling |
+| `geometric_checks.py` | Pre-flight geometric feasibility filter |
+| `nx_interface.py` | NX automation module (stub + NXOpen template) |
+| `requirements.txt` | Python dependencies |
+| `OPTIMIZATION_STRATEGY.md` | Full strategy document |
+| `results/` | Output directory (CSV, JSON, Optuna DB) |
 
-## Open Items Before Execution
+## Output Files
 
-1. Beam web length needed (geometric feasibility of hole patterns)
-2. Displacement extraction method (sensor vs .op2 node parsing)
-3. Stress extraction scope (whole model vs element group)
-4. Baseline stress measurement
-5. SOL 101 runtime benchmark
-6. Corner test validation (16 bound combinations)
+After a study run, `results/` will contain:
+
+| File | Description |
+|------|-------------|
+| `doe_results.csv` | All trials — DVs, objectives, constraints, status |
+| `doe_summary.json` | Study metadata, statistics, best feasible design |
+| `optuna_study.db` | SQLite database (Optuna persistence, resume support) |
+
+## Architecture
+
+```
+run_doe.py
+  ├── sampling.py          Generate 50 LHS + 1 baseline
+  │     └── scipy.stats.qmc.LatinHypercube
+  ├── geometric_checks.py  Pre-flight feasibility filter
+  │     ├── Hole overlap:  span/(n-1) - d ≥ 30mm
+  │     └── Web clearance: 500 - 2·face - d > 0
+  ├── nx_interface.py      NX solver (stub or NXOpen)
+  │     ├── NXStubSolver   → synthetic results (development)
+  │     └── NXOpenSolver   → real NX Nastran SOL 101 (production)
+  └── Optuna study
+        ├── SQLite storage (resume support)
+        ├── Enqueued trials (deterministic LHS)
+        └── Deb's feasibility rules (constraint interface)
+```
+
+## Pipeline per Trial
+
+```
+Trial N
+  │
+  ├── 1. Suggest DVs (from enqueued LHS point)
+  │
+  ├── 2. Geometric pre-check
+  │     ├── FAIL → record as infeasible, skip NX
+  │     └── PASS ↓
+  │
+  ├── 3. NX evaluation
+  │     ├── Update expressions (exact NX names)
+  │     ├── Rebuild model
+  │     ├── Solve SOL 101
+  │     └── Extract: mass (p173), displacement, stress
+  │
+  ├── 4. Record results + constraint violations
+  │
+  └── 5. Log to CSV + Optuna DB
+```
+
+## Phase 1 → Phase 2 Gate Criteria
+
+Before proceeding to Phase 2 (TPE optimization), these checks must pass:
+
+| Check | Threshold | Action if FAIL |
+|-------|-----------|----------------|
+| Feasible points found | ≥ 5 | Expand bounds or relax constraints (escalate to CEO) |
+| NX solve success rate | ≥ 80% | Investigate failures, fix model, re-run |
+| No systematic crashes at bounds | Visual check | Tighten bounds away from failure region |
+
+## Key Design Decisions
+
+- **Baseline as Trial 0** — LAC lesson: always validate baseline first
+- **Pre-flight geometric filter** — catches infeasible geometry before NX (saves compute, avoids crashes)
+- **Stratified integer sampling** — ensures all 11 hole_count levels (5-15) are covered
+- **Deb's feasibility rules** — no penalty weight tuning; feasible always beats infeasible
+- **SQLite persistence** — study can be interrupted and resumed
+- **No surrogates** — LAC lesson: direct FEA via TPE beats surrogate + L-BFGS
+
+## NX Integration Notes
+
+The `nx_interface.py` module provides:
+- **`NXStubSolver`** — synthetic results from simplified beam mechanics (for pipeline testing)
+- **`NXOpenSolver`** — template for real NXOpen Python API integration (to be completed on Windows)
+
+Expression names are exact from binary introspection. Critical: `beam_lenght` has a typo in NX (no 'h') — use exact spelling.
+
+### Outputs to extract from NX:
+| Output | Source | Unit | Notes |
+|--------|--------|------|-------|
+| Mass | Expression `p173` | kg | Direct read |
+| Tip displacement | SOL 101 results | mm | TBD: sensor or .op2 parsing |
+| Max VM stress | SOL 101 results | MPa | ⚠️ pyNastran returns kPa — divide by 1000 |
 
 ## References
 
-- [BREAKDOWN.md](../../BREAKDOWN.md) — Tech Lead's technical analysis
-- [DECISIONS.md](../../DECISIONS.md) — Decision log
-- [CONTEXT.md](../../CONTEXT.md) — Project context
+- [OPTIMIZATION_STRATEGY.md](./OPTIMIZATION_STRATEGY.md) — Full strategy document
+- [../../BREAKDOWN.md](../../BREAKDOWN.md) — Tech Lead's technical analysis
+- [../../DECISIONS.md](../../DECISIONS.md) — Decision log
+- [../../CONTEXT.md](../../CONTEXT.md) — Project context & expression map
 
 ---
 
-*Created by ⚡ Optimizer Agent | 2026-02-09*
+*Code: 🏗️ Study Builder (Technical Lead) | Strategy: ⚡ Optimizer Agent | 2026-02-10*
diff --git a/projects/hydrotech-beam/studies/01_doe_landscape/geometric_checks.py b/projects/hydrotech-beam/studies/01_doe_landscape/geometric_checks.py
new file mode 100644
index 00000000..70910b29
--- /dev/null
+++ b/projects/hydrotech-beam/studies/01_doe_landscape/geometric_checks.py
@@ -0,0 +1,205 @@
+"""Geometric feasibility pre-flight checks for Hydrotech Beam.
+
+Validates trial design variable combinations against physical geometry
+constraints BEFORE sending to NX. Catches infeasible combos that would
+cause NX rebuild failures or physically impossible geometries.
+
+References:
+    OPTIMIZATION_STRATEGY.md §4.2 — Hole overlap analysis
+    Auditor review 2026-02-10 — Corrected spacing formula
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import NamedTuple
+
+
+# ---------------------------------------------------------------------------
+# Fixed geometry constants (from NX introspection — CONTEXT.md)
+# ---------------------------------------------------------------------------
+BEAM_HALF_HEIGHT: float = 250.0  # mm — fixed, expression `beam_half_height`
+HOLE_SPAN: float = 4000.0  # mm — expression `p6`, total distribution length
+MIN_LIGAMENT: float = 30.0  # mm — minimum material between holes (mesh + structural)
+
+
+class FeasibilityResult(NamedTuple):
+    """Result of a geometric feasibility check."""
+
+    feasible: bool
+    reason: str
+    ligament: float  # mm — material between adjacent holes
+    web_clearance: float  # mm — clearance between hole edge and faces
+
+
+@dataclass(frozen=True)
+class DesignPoint:
+    """A single design variable combination."""
+
+    beam_half_core_thickness: float  # mm — DV1
+    beam_face_thickness: float  # mm — DV2
+    holes_diameter: float  # mm — DV3
+    hole_count: int  # — DV4
+
+
+def compute_ligament(holes_diameter: float, hole_count: int) -> float:
+    """Compute ligament width (material between adjacent holes).
+
+    The NX pattern places `n` holes across `hole_span` mm using `n-1`
+    intervals (holes at both endpoints of the span).
+
+    Formula (corrected per auditor review):
+        spacing = hole_span / (hole_count - 1)
+        ligament = spacing - holes_diameter
+
+    Args:
+        holes_diameter: Hole diameter in mm.
+        hole_count: Number of holes (integer ≥ 2).
+
+    Returns:
+        Ligament width in mm. Negative means overlap.
+    """
+    if hole_count < 2:
+        # Single hole — no overlap possible, return large ligament
+        return HOLE_SPAN
+    spacing = HOLE_SPAN / (hole_count - 1)
+    return spacing - holes_diameter
+
+
+def compute_web_clearance(
+    beam_face_thickness: float, holes_diameter: float
+) -> float:
+    """Compute clearance between hole edge and face sheets.
+
+    Total beam height = 2 × beam_half_height = 500 mm (fixed).
+    Web clear height = total_height - 2 × face_thickness.
+    Clearance = web_clear_height - holes_diameter.
+
+    Args:
+        beam_face_thickness: Face sheet thickness in mm.
+        holes_diameter: Hole diameter in mm.
+
+    Returns:
+        Web clearance in mm. ≤ 0 means hole doesn't fit.
+    """
+    total_height = 2.0 * BEAM_HALF_HEIGHT  # 500 mm
+    web_clear_height = total_height - 2.0 * beam_face_thickness
+    return web_clear_height - holes_diameter
+
+
+def check_feasibility(point: DesignPoint) -> FeasibilityResult:
+    """Run all geometric feasibility checks on a design point.
+
+    Checks (in order):
+        1. Hole overlap — ligament between adjacent holes ≥ MIN_LIGAMENT
+        2. Web clearance — hole fits within the web (between face sheets)
+
+    Args:
+        point: Design variable combination to check.
+
+    Returns:
+        FeasibilityResult with feasibility status and diagnostic info.
+    """
+    ligament = compute_ligament(point.holes_diameter, point.hole_count)
+    web_clearance = compute_web_clearance(
+        point.beam_face_thickness, point.holes_diameter
+    )
+
+    # Check 1: Hole overlap
+    if ligament < MIN_LIGAMENT:
+        return FeasibilityResult(
+            feasible=False,
+            reason=(
+                f"Hole overlap: ligament={ligament:.1f}mm < "
+                f"{MIN_LIGAMENT}mm minimum "
+                f"(d={point.holes_diameter}mm, n={point.hole_count})"
+            ),
+            ligament=ligament,
+            web_clearance=web_clearance,
+        )
+
+    # Check 2: Web clearance
+    if web_clearance <= 0:
+        return FeasibilityResult(
+            feasible=False,
+            reason=(
+                f"Hole exceeds web: clearance={web_clearance:.1f}mm ≤ 0 "
+                f"(face={point.beam_face_thickness}mm, "
+                f"d={point.holes_diameter}mm)"
+            ),
+            ligament=ligament,
+            web_clearance=web_clearance,
+        )
+
+    return FeasibilityResult(
+        feasible=True,
+        reason="OK",
+        ligament=ligament,
+        web_clearance=web_clearance,
+    )
+
+
+def check_feasibility_from_values(
+    beam_half_core_thickness: float,
+    beam_face_thickness: float,
+    holes_diameter: float,
+    hole_count: int,
+) -> FeasibilityResult:
+    """Convenience wrapper — check feasibility from raw DV values.
+
+    Args:
+        beam_half_core_thickness: Core half-thickness in mm (DV1).
+        beam_face_thickness: Face sheet thickness in mm (DV2).
+        holes_diameter: Hole diameter in mm (DV3).
+        hole_count: Number of holes (DV4).
+
+    Returns:
+        FeasibilityResult with feasibility status and diagnostic info.
+    """
+    point = DesignPoint(
+        beam_half_core_thickness=beam_half_core_thickness,
+        beam_face_thickness=beam_face_thickness,
+        holes_diameter=holes_diameter,
+        hole_count=hole_count,
+    )
+    return check_feasibility(point)
+
+
+# ---------------------------------------------------------------------------
+# Quick self-test
+# ---------------------------------------------------------------------------
+if __name__ == "__main__":
+    # Baseline: should be feasible
+    baseline = DesignPoint(
+        beam_half_core_thickness=25.162,
+        beam_face_thickness=21.504,
+        holes_diameter=300.0,
+        hole_count=10,
+    )
+    result = check_feasibility(baseline)
+    print(f"Baseline: {result}")
+    assert result.feasible, f"Baseline should be feasible! Got: {result}"
+
+    # Worst case: n=15, d=450 — should be infeasible (overlap)
+    worst_overlap = DesignPoint(
+        beam_half_core_thickness=25.0,
+        beam_face_thickness=20.0,
+        holes_diameter=450.0,
+        hole_count=15,
+    )
+    result = check_feasibility(worst_overlap)
+    print(f"Worst overlap: {result}")
+    assert not result.feasible, "n=15, d=450 should be infeasible"
+
+    # Web clearance fail: face=40, d=450
+    web_fail = DesignPoint(
+        beam_half_core_thickness=25.0,
+        beam_face_thickness=40.0,
+        holes_diameter=450.0,
+        hole_count=5,
+    )
+    result = check_feasibility(web_fail)
+    print(f"Web clearance fail: {result}")
+    assert not result.feasible, "face=40, d=450 should fail web clearance"
+
+    print("\nAll self-tests passed ✓")
diff --git a/projects/hydrotech-beam/studies/01_doe_landscape/nx_interface.py b/projects/hydrotech-beam/studies/01_doe_landscape/nx_interface.py
new file mode 100644
index 00000000..0b483e01
--- /dev/null
+++ b/projects/hydrotech-beam/studies/01_doe_landscape/nx_interface.py
@@ -0,0 +1,398 @@
+"""NX automation interface for Hydrotech Beam optimization.
+
+Stub/template module for NXOpen Python API integration. The actual NX
+automation runs on Windows (dalidou node) via Syncthing-synced model files.
+
+This module defines the interface contract. The NXOpen-specific implementation
+will be filled in when running on the Windows side.
+
+NX Expression Names (confirmed via binary introspection — CONTEXT.md):
+    Design Variables:
+        - beam_half_core_thickness  (mm, continuous)
+        - beam_face_thickness       (mm, continuous)
+        - holes_diameter            (mm, continuous)
+        - hole_count                (integer, links to Pattern_p7)
+    Outputs:
+        - p173                      (mass in kg, body_property147.mass)
+    Fixed:
+        - beam_lenght               (⚠️ TYPO in NX — no 'h', 5000 mm)
+        - beam_half_height          (250 mm)
+        - beam_half_width           (150 mm)
+
+References:
+    CONTEXT.md — Full expression map
+    OPTIMIZATION_STRATEGY.md §8.2 — Extractor requirements
+"""
+
+from __future__ import annotations
+
+import logging
+from dataclasses import dataclass
+from typing import Protocol
+
+logger = logging.getLogger(__name__)
+
+
+# ---------------------------------------------------------------------------
+# Data types
+# ---------------------------------------------------------------------------
+@dataclass(frozen=True)
+class TrialInput:
+    """Design variable values for a single trial."""
+
+    beam_half_core_thickness: float  # mm — DV1
+    beam_face_thickness: float  # mm — DV2
+    holes_diameter: float  # mm — DV3
+    hole_count: int  # — DV4
+
+
+@dataclass
+class TrialResult:
+    """Results extracted from NX after a trial solve.
+
+    All values populated after a successful SOL 101 solve.
+    On failure, success=False and error_message explains the failure.
+    """
+
+    success: bool
+    mass: float = float("nan")  # kg — from expression `p173`
+    tip_displacement: float = float("nan")  # mm — from SOL 101 results
+    max_von_mises: float = float("nan")  # MPa — from SOL 101 results
+    error_message: str = ""
+
+
+# ---------------------------------------------------------------------------
+# NX expression name constants
+# ---------------------------------------------------------------------------
+# ⚠️ These are EXACT NX expression names from binary introspection.
+# Do NOT change spelling — `beam_lenght` has a typo (no 'h') in NX.
+EXPR_HALF_CORE_THICKNESS = "beam_half_core_thickness"
+EXPR_FACE_THICKNESS = "beam_face_thickness"
+EXPR_HOLES_DIAMETER = "holes_diameter"
+EXPR_HOLE_COUNT = "hole_count"
+EXPR_MASS = "p173"  # body_property147.mass, kg
+EXPR_BEAM_LENGTH = "beam_lenght"  # ⚠️ TYPO IN NX — intentional
+
+
+# ---------------------------------------------------------------------------
+# Interface protocol
+# ---------------------------------------------------------------------------
+class NXSolverInterface(Protocol):
+    """Protocol for NX solver backends.
+
+    Implementors must provide the full pipeline:
+    1. Update expressions → 2. Rebuild model → 3. Solve SOL 101 → 4. Extract results
+    """
+
+    def evaluate(self, trial_input: TrialInput) -> TrialResult:
+        """Run a full NX evaluation cycle for one trial.
+
+        Args:
+            trial_input: Design variable values.
+
+        Returns:
+            TrialResult with extracted outputs or failure info.
+        """
+        ...
+
+    def close(self) -> None:
+        """Clean up NX session resources.
+
+        ⚠️ LAC CRITICAL: NEVER kill NX processes directly.
+        Use NXSessionManager.close_nx_if_allowed() only.
+        """
+        ...
+
+
+# ---------------------------------------------------------------------------
+# Stub implementation (for development/testing without NX)
+# ---------------------------------------------------------------------------
+class NXStubSolver:
+    """Stub NX solver for development and testing.
+
+    Returns synthetic results based on simple analytical approximations
+    of the beam behavior. NOT physically accurate — use only for
+    testing the optimization pipeline.
+
+    The stub uses rough scaling relationships:
+        - Mass ∝ (core + face) and inversely with hole area
+        - Displacement ∝ 1/I where I depends on core and face thickness
+        - Stress ∝ M*y/I (bending stress approximation)
+    """
+
+    def __init__(self) -> None:
+        """Initialize stub solver."""
+        logger.warning(
+            "Using NX STUB solver — results are synthetic approximations. "
+            "Replace with NXOpenSolver for real evaluations."
+        )
+
+    def evaluate(self, trial_input: TrialInput) -> TrialResult:
+        """Return synthetic results based on simplified beam mechanics.
+
+        Args:
+            trial_input: Design variable values.
+
+        Returns:
+            TrialResult with approximate values.
+        """
+        try:
+            return self._compute_approximate(trial_input)
+        except Exception as e:
+            logger.error("Stub evaluation failed: %s", e)
+            return TrialResult(
+                success=False,
+                error_message=f"Stub evaluation error: {e}",
+            )
+
+    def _compute_approximate(self, inp: TrialInput) -> TrialResult:
+        """Simple analytical approximation of beam response.
+
+        This is a ROUGH approximation for pipeline testing only.
+        Real physics requires NX Nastran SOL 101.
+        """
+        import math
+
+        # Geometry
+        L = 5000.0  # mm — beam length
+        b = 300.0  # mm — beam width (2 × beam_half_width)
+        h_core = inp.beam_half_core_thickness  # mm — half core
+        t_face = inp.beam_face_thickness  # mm — face thickness
+        d_hole = inp.holes_diameter  # mm
+        n_holes = inp.hole_count
+
+        # Total height and section properties (simplified I-beam)
+        h_total = 500.0  # mm — 2 × beam_half_height (fixed)
+
+        # Approximate second moment of area (sandwich beam)
+        # I ≈ b*h_total^3/12 - b*(h_total-2*t_face)^3/12 + web contribution
+        h_inner = h_total - 2.0 * t_face
+        I_section = (b * h_total**3 / 12.0) - (b * max(h_inner, 0.0) ** 3 / 12.0)
+
+        # Add core contribution
+        I_section += 2.0 * h_core * h_total**2 / 4.0  # approximate
+
+        # Hole area reduction (mass)
+        hole_area = n_holes * math.pi * (d_hole / 2.0) ** 2  # mm²
+
+        # Approximate mass (steel: 7.3 g/cm³ = 7.3e-6 kg/mm³)
+        rho = 7.3e-6  # kg/mm³
+        # Gross cross-section area (very simplified)
+        A_gross = 2.0 * b * t_face + 2.0 * h_core * h_total
+        # Remove holes from web
+        web_thickness = 2.0 * h_core  # approximate web thickness
+        A_holes = n_holes * math.pi * (d_hole / 2.0) ** 2
+        V_solid = A_gross * L
+        V_holes = A_holes * web_thickness
+        mass = rho * (V_solid - min(V_holes, V_solid * 0.8))
+
+        # Approximate tip displacement (cantilever, point load)
+        # δ = PL³/(3EI)
+        P = 10000.0 * 9.81  # 10,000 kgf → N
+        E = 200000.0  # MPa (steel)
+        if I_section > 0:
+            displacement = P * L**3 / (3.0 * E * I_section)
+        else:
+            displacement = 9999.0
+
+        # Approximate max bending stress
+        # σ = M*y/I where M = P*L, y = h_total/2
+        M_max = P * L  # N·mm
+        y_max = h_total / 2.0
+        if I_section > 0:
+            stress = M_max * y_max / I_section  # MPa
+        else:
+            stress = 9999.0
+
+        return TrialResult(
+            success=True,
+            mass=mass,
+            tip_displacement=displacement,
+            max_von_mises=stress,
+        )
+
+    def close(self) -> None:
+        """No-op for stub solver."""
+        logger.info("Stub solver closed.")
+
+
+# ---------------------------------------------------------------------------
+# NXOpen implementation template (to be completed on Windows/dalidou)
+# ---------------------------------------------------------------------------
+class NXOpenSolver:
+    """Real NXOpen-based solver — TEMPLATE, not yet functional.
+
+    This class provides the correct structure for NXOpen integration.
+    Expression update code uses the exact names from binary introspection.
+
+    To complete:
+        1. Set NX_MODEL_DIR to the Syncthing-synced model directory
+        2. Implement _open_session() with NXOpen.Session
+        3. Implement _solve() to trigger SOL 101
+        4. Implement _extract_displacement() and _extract_stress()
+           from .op2 results or NX result sensors
+    """
+
+    def __init__(self, model_dir: str) -> None:
+        """Initialize NXOpen solver.
+
+        Args:
+            model_dir: Path to directory containing Beam.prt, etc.
+        """
+        self.model_dir = model_dir
+        self._session = None  # NXOpen.Session
+        self._part = None  # NXOpen.Part
+        logger.info("NXOpenSolver initialized with model_dir=%s", model_dir)
+
+    def evaluate(self, trial_input: TrialInput) -> TrialResult:
+        """Full NX evaluation pipeline.
+
+        Pipeline:
+            1. Update expressions (beam_half_core_thickness, etc.)
+            2. Rebuild model (triggers re-mesh of idealized part)
+            3. Solve SOL 101
+            4. Extract mass (p173), displacement, stress
+        """
+        raise NotImplementedError(
+            "NXOpenSolver.evaluate() is a template — implement on Windows "
+            "with NXOpen Python API. See docstring for pipeline."
+        )
+
+    def _update_expressions(self, trial_input: TrialInput) -> None:
+        """Update NX expressions for a trial.
+
+        ⚠️ Expression names are EXACT from binary introspection.
+        ⚠️ `beam_lenght` has a typo (no 'h') — do NOT correct it.
+
+        This is the correct NXOpen code pattern (to be run on Windows):
+
+        ```python
+        import NXOpen
+
+        session = NXOpen.Session.GetSession()
+        part = session.Parts.Work
+
+        # Update design variables
+        expressions = {
+            "beam_half_core_thickness": trial_input.beam_half_core_thickness,
+            "beam_face_thickness": trial_input.beam_face_thickness,
+            "holes_diameter": trial_input.holes_diameter,
+            "hole_count": trial_input.hole_count,
+        }
+
+        for expr_name, value in expressions.items():
+            expr = part.Expressions.FindObject(expr_name)
+            unit = expr.Units
+            part.Expressions.EditWithUnits(expr, unit, str(value))
+
+        # Rebuild (update model to reflect new expressions)
+        session.UpdateManager.DoUpdate(
+            session.SetUndoMark(NXOpen.Session.MarkVisibility.Invisible, "Update")
+        )
+        ```
+        """
+        raise NotImplementedError("Template — implement with NXOpen")
+
+    def _solve(self) -> bool:
+        """Trigger NX Nastran SOL 101 solve.
+
+        ```python
+        # Open the sim file
+        sim_part = session.Parts.OpenDisplay(
+            os.path.join(self.model_dir, "Beam_sim1.sim"), None
+        )
+
+        # Get the solution and solve
+        sim_simulation = sim_part.Simulation
+        solution = sim_simulation.Solutions[0]  # First solution
+        solution.Solve()
+
+        # Check solve status
+        return solution.SolveStatus == NXOpen.CAE.Solution.Status.Solved
+        ```
+        """
+        raise NotImplementedError("Template — implement with NXOpen")
+
+    def _extract_mass(self) -> float:
+        """Extract mass from NX expression p173.
+
+        ```python
+        mass_expr = part.Expressions.FindObject("p173")
+        return mass_expr.Value  # kg
+        ```
+        """
+        raise NotImplementedError("Template — implement with NXOpen")
+
+    def _extract_displacement(self) -> float:
+        """Extract tip displacement from SOL 101 results.
+
+        Options (TBD — need to determine best approach):
+        1. NX result sensor at beam tip → read value directly
+        2. Parse .op2 file with pyNastran → find max displacement at tip nodes
+        3. NX post-processing API → query displacement field
+
+        Returns:
+            Tip displacement in mm.
+        """
+        raise NotImplementedError(
+            "Template — extraction method TBD. "
+            "Options: result sensor, .op2 parsing, or NX post-processing API."
+        )
+
+    def _extract_stress(self) -> float:
+        """Extract max von Mises stress from SOL 101 results.
+
+        ⚠️ LAC LESSON: pyNastran returns stress in kPa for NX kg-mm-s
+        unit system. Divide by 1000 for MPa.
+
+        Options (TBD):
+        1. NX result sensor for max VM stress → read value directly
+        2. Parse .op2 with pyNastran → max elemental nodal VM stress
+        3. NX post-processing API → query stress field
+
+        Returns:
+            Max von Mises stress in MPa.
+        """
+        raise NotImplementedError(
+            "Template — extraction method TBD. "
+            "Remember: pyNastran stress is in kPa → divide by 1000 for MPa."
+        )
+
+    def close(self) -> None:
+        """Close NX session gracefully.
+
+        ⚠️ LAC CRITICAL: NEVER kill NX processes directly.
+        Use NXSessionManager.close_nx_if_allowed() only.
+        If NX hangs, implement a timeout (10 min per trial) and let
+        NX time out gracefully.
+        """
+        raise NotImplementedError("Template — implement graceful shutdown")
+
+
+# ---------------------------------------------------------------------------
+# Factory
+# ---------------------------------------------------------------------------
+def create_solver(
+    backend: str = "stub",
+    model_dir: str = "",
+) -> NXStubSolver | NXOpenSolver:
+    """Create an NX solver instance.
+
+    Args:
+        backend: "stub" for development, "nxopen" for real NX (Windows only).
+        model_dir: Path to NX model directory (required for nxopen backend).
+
+    Returns:
+        Solver instance implementing the NXSolverInterface protocol.
+
+    Raises:
+        ValueError: If backend is unknown.
+    """
+    if backend == "stub":
+        return NXStubSolver()
+    elif backend == "nxopen":
+        if not model_dir:
+            raise ValueError("model_dir required for nxopen backend")
+        return NXOpenSolver(model_dir=model_dir)
+    else:
+        raise ValueError(f"Unknown backend: {backend!r}. Use 'stub' or 'nxopen'.")
diff --git a/projects/hydrotech-beam/studies/01_doe_landscape/requirements.txt b/projects/hydrotech-beam/studies/01_doe_landscape/requirements.txt
new file mode 100644
index 00000000..962559ff
--- /dev/null
+++ b/projects/hydrotech-beam/studies/01_doe_landscape/requirements.txt
@@ -0,0 +1,7 @@
+# Hydrotech Beam — Phase 1 LHS DoE Study
+# Python 3.10+
+
+optuna>=3.5.0
+scipy>=1.10.0
+numpy>=1.24.0
+pandas>=2.0.0
diff --git a/projects/hydrotech-beam/studies/01_doe_landscape/run_doe.py b/projects/hydrotech-beam/studies/01_doe_landscape/run_doe.py
new file mode 100644
index 00000000..cda9e7df
--- /dev/null
+++ b/projects/hydrotech-beam/studies/01_doe_landscape/run_doe.py
@@ -0,0 +1,659 @@
+#!/usr/bin/env python3
+"""Phase 1 LHS DoE — Hydrotech Beam Optimization.
+
+Main entry point for the Design of Experiments landscape mapping study.
+Generates 50 LHS sample points + 1 baseline (Trial 0), runs pre-flight
+geometric checks, evaluates via NX (or stub), and manages the Optuna study.
+
+Usage:
+    # Dry run with stub solver (development/testing):
+    python run_doe.py --backend stub --study-name hydrotech_doe_dev
+
+    # Real run with NXOpen (on Windows/dalidou):
+    python run_doe.py --backend nxopen --model-dir /path/to/nx/models
+
+    # Resume interrupted study:
+    python run_doe.py --backend stub --study-name hydrotech_doe_dev --resume
+
+References:
+    OPTIMIZATION_STRATEGY.md — Full strategy document
+    CONTEXT.md — Confirmed expression names and values
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import logging
+import os
+import sys
+import time
+from datetime import datetime, timezone
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+import optuna
+
+from geometric_checks import (
+    DesignPoint,
+    FeasibilityResult,
+    check_feasibility,
+)
+from nx_interface import TrialInput, TrialResult, create_solver
+from sampling import DV_DEFINITIONS, generate_lhs_samples, points_to_dicts
+
+# ---------------------------------------------------------------------------
+# Configuration
+# ---------------------------------------------------------------------------
+DEFAULT_STUDY_NAME = "hydrotech_beam_doe_phase1"
+DEFAULT_N_SAMPLES = 50
+DEFAULT_SEED = 42
+DEFAULT_RESULTS_DIR = "results"
+DEFAULT_DB_PATH = "results/optuna_study.db"
+
+# Constraint limits (hard limits — OPTIMIZATION_STRATEGY.md §1.2)
+DISPLACEMENT_LIMIT = 10.0  # mm
+STRESS_LIMIT = 130.0  # MPa
+
+# Infeasible placeholder values (for geometric pre-check failures)
+INFEASIBLE_MASS = 99999.0  # kg
+INFEASIBLE_DISPLACEMENT = 99999.0  # mm
+INFEASIBLE_STRESS = 99999.0  # MPa
+
+logger = logging.getLogger(__name__)
+
+
+# ---------------------------------------------------------------------------
+# Optuna constraint callback
+# ---------------------------------------------------------------------------
+def constraints_func(trial: optuna.trial.FrozenTrial) -> list[float]:
+    """Compute constraint violations for Deb's feasibility rules.
+
+    Returns list of constraint values where:
+        value ≤ 0 → feasible (constraint satisfied)
+        value > 0 → infeasible (constraint violated)
+
+    This implements Deb's feasibility rules (Deb 2000):
+        1. Feasible vs feasible → lower objective wins
+        2. Feasible vs infeasible → feasible wins
+        3. Infeasible vs infeasible → lower total violation wins
+
+    References:
+        OPTIMIZATION_STRATEGY.md §3.2
+
+    Args:
+        trial: Optuna frozen trial with user attributes set.
+
+    Returns:
+        [displacement_violation, stress_violation]
+    """
+    disp = trial.user_attrs.get("tip_displacement", INFEASIBLE_DISPLACEMENT)
+    stress = trial.user_attrs.get("max_von_mises", INFEASIBLE_STRESS)
+    return [
+        disp - DISPLACEMENT_LIMIT,  # ≤ 0 means displacement ≤ 10 mm
+        stress - STRESS_LIMIT,  # ≤ 0 means stress ≤ 130 MPa
+    ]
+
+
+# ---------------------------------------------------------------------------
+# Trial evaluation
+# ---------------------------------------------------------------------------
+def evaluate_trial(
+    trial: optuna.Trial,
+    solver: Any,
+) -> float:
+    """Evaluate a single trial: geometric check → NX solve → extract.
+
+    Args:
+        trial: Optuna trial (with parameters already suggested/enqueued).
+        solver: NX solver instance (stub or real).
+
+    Returns:
+        Objective value (mass in kg). Returns INFEASIBLE_MASS for
+        geometrically infeasible or failed evaluations.
+    """
+    # Extract design variables from trial
+    dv1 = trial.suggest_float(
+        "beam_half_core_thickness", 10.0, 40.0,
+    )
+    dv2 = trial.suggest_float(
+        "beam_face_thickness", 10.0, 40.0,
+    )
+    dv3 = trial.suggest_float(
+        "holes_diameter", 150.0, 450.0,
+    )
+    dv4 = trial.suggest_int(
+        "hole_count", 5, 15,
+    )
+
+    trial_num = trial.number
+    logger.info(
+        "Trial %d: DV1=%.2f, DV2=%.2f, DV3=%.1f, DV4=%d",
+        trial_num, dv1, dv2, dv3, dv4,
+    )
+
+    # Store DVs in user attributes for logging
+    trial.set_user_attr("beam_half_core_thickness", dv1)
+    trial.set_user_attr("beam_face_thickness", dv2)
+    trial.set_user_attr("holes_diameter", dv3)
+    trial.set_user_attr("hole_count", dv4)
+
+    # Pre-flight geometric check
+    point = DesignPoint(
+        beam_half_core_thickness=dv1,
+        beam_face_thickness=dv2,
+        holes_diameter=dv3,
+        hole_count=dv4,
+    )
+    geo_result: FeasibilityResult = check_feasibility(point)
+    trial.set_user_attr("geo_feasible", geo_result.feasible)
+    trial.set_user_attr("ligament", geo_result.ligament)
+    trial.set_user_attr("web_clearance", geo_result.web_clearance)
+
+    if not geo_result.feasible:
+        logger.warning(
+            "Trial %d: GEOMETRICALLY INFEASIBLE — %s",
+            trial_num, geo_result.reason,
+        )
+        trial.set_user_attr("status", "geo_infeasible")
+        trial.set_user_attr("geo_reason", geo_result.reason)
+        trial.set_user_attr("tip_displacement", INFEASIBLE_DISPLACEMENT)
+        trial.set_user_attr("max_von_mises", INFEASIBLE_STRESS)
+        trial.set_user_attr("mass", INFEASIBLE_MASS)
+        return INFEASIBLE_MASS
+
+    # NX evaluation
+    trial_input = TrialInput(
+        beam_half_core_thickness=dv1,
+        beam_face_thickness=dv2,
+        holes_diameter=dv3,
+        hole_count=dv4,
+    )
+
+    t_start = time.monotonic()
+    nx_result: TrialResult = solver.evaluate(trial_input)
+    t_elapsed = time.monotonic() - t_start
+
+    trial.set_user_attr("solve_time_s", round(t_elapsed, 2))
+
+    if not nx_result.success:
+        logger.error(
+            "Trial %d: NX SOLVE FAILED — %s (%.1fs)",
+            trial_num, nx_result.error_message, t_elapsed,
+        )
+        trial.set_user_attr("status", "solve_failed")
+        trial.set_user_attr("error_message", nx_result.error_message)
+        trial.set_user_attr("tip_displacement", INFEASIBLE_DISPLACEMENT)
+        trial.set_user_attr("max_von_mises", INFEASIBLE_STRESS)
+        trial.set_user_attr("mass", INFEASIBLE_MASS)
+        return INFEASIBLE_MASS
+
+    # Record successful results
+    trial.set_user_attr("status", "solved")
+    trial.set_user_attr("mass", nx_result.mass)
+    trial.set_user_attr("tip_displacement", nx_result.tip_displacement)
+    trial.set_user_attr("max_von_mises", nx_result.max_von_mises)
+
+    # Check constraint feasibility
+    disp_ok = nx_result.tip_displacement <= DISPLACEMENT_LIMIT
+    stress_ok = nx_result.max_von_mises <= STRESS_LIMIT
+    trial.set_user_attr("displacement_feasible", disp_ok)
+    trial.set_user_attr("stress_feasible", stress_ok)
+    trial.set_user_attr("fully_feasible", disp_ok and stress_ok)
+
+    logger.info(
+        "Trial %d: mass=%.2f kg, disp=%.2f mm%s, stress=%.1f MPa%s (%.1fs)",
+        trial_num,
+        nx_result.mass,
+        nx_result.tip_displacement,
+        " ✓" if disp_ok else " ✗",
+        nx_result.max_von_mises,
+        " ✓" if stress_ok else " ✗",
+        t_elapsed,
+    )
+
+    return nx_result.mass
+
+
+# ---------------------------------------------------------------------------
+# Results export
+# ---------------------------------------------------------------------------
+def export_csv(study: optuna.Study, output_path: str) -> None:
+    """Export all trial results to CSV.
+
+    Columns: trial_number, status, beam_half_core_thickness, beam_face_thickness,
+    holes_diameter, hole_count, mass, tip_displacement, max_von_mises,
+    geo_feasible, displacement_feasible, stress_feasible, fully_feasible,
+    ligament, web_clearance, solve_time_s
+
+    Args:
+        study: Completed Optuna study.
+        output_path: Path to output CSV file.
+    """
+    fieldnames = [
+        "trial_number",
+        "status",
+        "beam_half_core_thickness",
+        "beam_face_thickness",
+        "holes_diameter",
+        "hole_count",
+        "mass_kg",
+        "tip_displacement_mm",
+        "max_von_mises_MPa",
+        "geo_feasible",
+        "displacement_feasible",
+        "stress_feasible",
+        "fully_feasible",
+        "ligament_mm",
+        "web_clearance_mm",
+        "solve_time_s",
+    ]
+
+    with open(output_path, "w", newline="") as f:
+        writer = csv.DictWriter(f, fieldnames=fieldnames)
+        writer.writeheader()
+
+        for trial in study.trials:
+            ua = trial.user_attrs
+            writer.writerow({
+                "trial_number": trial.number,
+                "status": ua.get("status", "unknown"),
+                "beam_half_core_thickness": ua.get("beam_half_core_thickness", ""),
+                "beam_face_thickness": ua.get("beam_face_thickness", ""),
+                "holes_diameter": ua.get("holes_diameter", ""),
+                "hole_count": ua.get("hole_count", ""),
+                "mass_kg": ua.get("mass", ""),
+                "tip_displacement_mm": ua.get("tip_displacement", ""),
+                "max_von_mises_MPa": ua.get("max_von_mises", ""),
+                "geo_feasible": ua.get("geo_feasible", ""),
+                "displacement_feasible": ua.get("displacement_feasible", ""),
+                "stress_feasible": ua.get("stress_feasible", ""),
+                "fully_feasible": ua.get("fully_feasible", ""),
+                "ligament_mm": ua.get("ligament", ""),
+                "web_clearance_mm": ua.get("web_clearance", ""),
+                "solve_time_s": ua.get("solve_time_s", ""),
+            })
+
+    logger.info("CSV results exported to %s (%d trials)", output_path, len(study.trials))
+
+
+def export_summary(study: optuna.Study, output_path: str) -> None:
+    """Export study summary as JSON.
+
+    Includes metadata, statistics, best feasible design, and
+    constraint satisfaction rates.
+
+    Args:
+        study: Completed Optuna study.
+        output_path: Path to output JSON file.
+    """
+    trials = study.trials
+    n_total = len(trials)
+
+    # Count by status
+    n_solved = sum(1 for t in trials if t.user_attrs.get("status") == "solved")
+    n_geo_infeasible = sum(
+        1 for t in trials if t.user_attrs.get("status") == "geo_infeasible"
+    )
+    n_failed = sum(
+        1 for t in trials if t.user_attrs.get("status") == "solve_failed"
+    )
+    n_feasible = sum(
+        1 for t in trials if t.user_attrs.get("fully_feasible", False)
+    )
+
+    # Best feasible trial
+    best_feasible = None
+    best_mass = float("inf")
+    for t in trials:
+        if t.user_attrs.get("fully_feasible", False):
+            mass = t.user_attrs.get("mass", float("inf"))
+            if mass < best_mass:
+                best_mass = mass
+                best_feasible = t
+
+    best_info = None
+    if best_feasible is not None:
+        ua = best_feasible.user_attrs
+        best_info = {
+            "trial_number": best_feasible.number,
+            "mass_kg": ua.get("mass"),
+            "tip_displacement_mm": ua.get("tip_displacement"),
+            "max_von_mises_MPa": ua.get("max_von_mises"),
+            "design_variables": {
+                "beam_half_core_thickness": ua.get("beam_half_core_thickness"),
+                "beam_face_thickness": ua.get("beam_face_thickness"),
+                "holes_diameter": ua.get("holes_diameter"),
+                "hole_count": ua.get("hole_count"),
+            },
+        }
+
+    summary = {
+        "study_name": study.study_name,
+        "phase": "Phase 1 — LHS DoE",
+        "project": "Hydrotech Beam Structural Optimization",
+        "timestamp": datetime.now(timezone.utc).isoformat(),
+        "configuration": {
+            "n_lhs_samples": DEFAULT_N_SAMPLES,
+            "seed": DEFAULT_SEED,
+            "baseline_included": True,
+            "algorithm": "Latin Hypercube Sampling (scipy.stats.qmc)",
+            "constraint_handling": "Deb's feasibility rules",
+            "displacement_limit_mm": DISPLACEMENT_LIMIT,
+            "stress_limit_MPa": STRESS_LIMIT,
+        },
+        "design_variables": [
+            {
+                "name": dv.name,
+                "nx_expression": dv.nx_expression,
+                "lower": dv.lower,
+                "upper": dv.upper,
+                "baseline": dv.baseline,
+                "type": "integer" if dv.is_integer else "continuous",
+            }
+            for dv in DV_DEFINITIONS
+        ],
+        "results": {
+            "total_trials": n_total,
+            "solved": n_solved,
+            "geo_infeasible": n_geo_infeasible,
+            "solve_failed": n_failed,
+            "fully_feasible": n_feasible,
+            "solve_success_rate": round(n_solved / max(n_total, 1) * 100, 1),
+            "feasibility_rate": round(n_feasible / max(n_solved, 1) * 100, 1),
+        },
+        "best_feasible": best_info,
+        "phase1_gate_check": {
+            "min_feasible_5": n_feasible >= 5,
+            "solve_success_80pct": (n_solved / max(n_total, 1)) >= 0.80,
+            "gate_passed": n_feasible >= 5 and (n_solved / max(n_total, 1)) >= 0.80,
+        },
+    }
+
+    with open(output_path, "w") as f:
+        json.dump(summary, f, indent=2)
+
+    logger.info("Summary exported to %s", output_path)
+
+
+# ---------------------------------------------------------------------------
+# Main study runner
+# ---------------------------------------------------------------------------
+def run_study(args: argparse.Namespace) -> None:
+    """Execute the Phase 1 LHS DoE study.
+
+    Steps:
+        1. Generate LHS sample points + baseline
+        2. Create/load Optuna study with SQLite storage
+        3. Enqueue all trials (baseline first, then LHS)
+        4. Run optimization (all trials evaluated via objective fn)
+        5. Export results to CSV and JSON
+
+    Args:
+        args: Parsed command-line arguments.
+    """
+    results_dir = Path(args.results_dir)
+    results_dir.mkdir(parents=True, exist_ok=True)
+
+    # -----------------------------------------------------------------------
+    # 1. Generate sample points
+    # -----------------------------------------------------------------------
+    logger.info("=" * 70)
+    logger.info("HYDROTECH BEAM — Phase 1 LHS DoE")
+    logger.info("=" * 70)
+
+    points = generate_lhs_samples(
+        n_samples=args.n_samples,
+        seed=args.seed,
+        include_baseline=True,
+    )
+    n_trials = len(points)
+    logger.info("Generated %d trial points (1 baseline + %d LHS)", n_trials, n_trials - 1)
+
+    # -----------------------------------------------------------------------
+    # 2. Create Optuna study
+    # -----------------------------------------------------------------------
+    db_path = results_dir / "optuna_study.db"
+    storage = f"sqlite:///{db_path}"
+
+    if args.resume:
+        logger.info("Resuming existing study: %s", args.study_name)
+        study = optuna.load_study(
+            study_name=args.study_name,
+            storage=storage,
+        )
+        logger.info("Loaded study with %d existing trials", len(study.trials))
+    else:
+        study = optuna.create_study(
+            study_name=args.study_name,
+            storage=storage,
+            direction="minimize",  # minimize mass
+            load_if_exists=False,
+            sampler=optuna.samplers.TPESampler(seed=args.seed),
+        )
+        logger.info("Created new study: %s (storage: %s)", args.study_name, db_path)
+
+        # Enqueue all trial points
+        trial_dicts = points_to_dicts(points)
+        for i, td in enumerate(trial_dicts):
+            study.enqueue_trial(td)
+            if i == 0:
+                logger.info("Enqueued Trial 0 (baseline): %s", td)
+
+        logger.info("Enqueued %d trials (1 baseline + %d LHS)", n_trials, n_trials - 1)
+
+    # -----------------------------------------------------------------------
+    # 3. Create solver
+    # -----------------------------------------------------------------------
+    solver = create_solver(
+        backend=args.backend,
+        model_dir=args.model_dir,
+    )
+
+    # -----------------------------------------------------------------------
+    # 4. Run all trials
+    # -----------------------------------------------------------------------
+    logger.info("-" * 70)
+    logger.info("Starting trial evaluations...")
+    logger.info("-" * 70)
+
+    t_study_start = time.monotonic()
+
+    # Suppress Optuna's verbose logging during trials
+    optuna.logging.set_verbosity(optuna.logging.WARNING)
+
+    study.optimize(
+        lambda trial: evaluate_trial(trial, solver),
+        n_trials=n_trials,
+        callbacks=[_progress_callback],
+    )
+
+    t_study_elapsed = time.monotonic() - t_study_start
+
+    # -----------------------------------------------------------------------
+    # 5. Export results
+    # -----------------------------------------------------------------------
+    logger.info("-" * 70)
+    logger.info("Study complete. Exporting results...")
+    logger.info("-" * 70)
+
+    csv_path = str(results_dir / "doe_results.csv")
+    json_path = str(results_dir / "doe_summary.json")
+
+    export_csv(study, csv_path)
+    export_summary(study, json_path)
+
+    # -----------------------------------------------------------------------
+    # 6. Print summary
+    # -----------------------------------------------------------------------
+    _print_summary(study, t_study_elapsed)
+
+    # Cleanup
+    solver.close()
+
+
+def _progress_callback(study: optuna.Study, trial: optuna.trial.FrozenTrial) -> None:
+    """Log progress after each trial."""
+    n_complete = len(study.trials)
+    status = trial.user_attrs.get("status", "unknown")
+    mass = trial.user_attrs.get("mass", "N/A")
+    feasible = trial.user_attrs.get("fully_feasible", False)
+    logger.info(
+        "  [%d/%d] Trial %d: status=%s, mass=%s, feasible=%s",
+        n_complete,
+        DEFAULT_N_SAMPLES + 1,
+        trial.number,
+        status,
+        f"{mass:.2f} kg" if isinstance(mass, (int, float)) else mass,
+        feasible,
+    )
+
+
+def _print_summary(study: optuna.Study, elapsed: float) -> None:
+    """Print a human-readable summary to stdout."""
+    trials = study.trials
+    n_total = len(trials)
+    n_solved = sum(1 for t in trials if t.user_attrs.get("status") == "solved")
+    n_geo_inf = sum(1 for t in trials if t.user_attrs.get("status") == "geo_infeasible")
+    n_failed = sum(1 for t in trials if t.user_attrs.get("status") == "solve_failed")
+    n_feasible = sum(1 for t in trials if t.user_attrs.get("fully_feasible", False))
+
+    print("\n" + "=" * 70)
+    print("PHASE 1 DoE — RESULTS SUMMARY")
+    print("=" * 70)
+    print(f"  Total trials:           {n_total}")
+    print(f"  Solved:                 {n_solved}")
+    print(f"  Geometrically infeasible: {n_geo_inf}")
+    print(f"  Solve failures:         {n_failed}")
+    print(f"  Fully feasible:         {n_feasible}")
+    print(f"  Solve success rate:     {n_solved/max(n_total,1)*100:.1f}%")
+    print(f"  Feasibility rate:       {n_feasible/max(n_solved,1)*100:.1f}%")
+    print(f"  Total time:             {elapsed:.1f}s ({elapsed/60:.1f} min)")
+
+    # Phase 1 → Phase 2 gate check
+    print("\n  GATE CHECK (Phase 1 → Phase 2):")
+    gate_feasible = n_feasible >= 5
+    gate_solve = (n_solved / max(n_total, 1)) >= 0.80
+    print(f"    ≥5 feasible points:     {'✓ PASS' if gate_feasible else '✗ FAIL'} ({n_feasible})")
+    print(f"    ≥80% solve success:     {'✓ PASS' if gate_solve else '✗ FAIL'} ({n_solved/max(n_total,1)*100:.0f}%)")
+    print(f"    GATE:                   {'✓ PASSED' if gate_feasible and gate_solve else '✗ BLOCKED'}")
+
+    # Best feasible
+    best_mass = float("inf")
+    best_trial = None
+    for t in trials:
+        if t.user_attrs.get("fully_feasible", False):
+            m = t.user_attrs.get("mass", float("inf"))
+            if m < best_mass:
+                best_mass = m
+                best_trial = t
+
+    if best_trial is not None:
+        ua = best_trial.user_attrs
+        print(f"\n  BEST FEASIBLE DESIGN (Trial {best_trial.number}):")
+        print(f"    Mass:           {ua['mass']:.2f} kg")
+        print(f"    Displacement:   {ua['tip_displacement']:.2f} mm (limit: {DISPLACEMENT_LIMIT} mm)")
+        print(f"    VM Stress:      {ua['max_von_mises']:.1f} MPa (limit: {STRESS_LIMIT} MPa)")
+        print(f"    Core thickness: {ua['beam_half_core_thickness']:.3f} mm")
+        print(f"    Face thickness: {ua['beam_face_thickness']:.3f} mm")
+        print(f"    Hole diameter:  {ua['holes_diameter']:.1f} mm")
+        print(f"    Hole count:     {ua['hole_count']}")
+    else:
+        print("\n  ⚠️ NO FEASIBLE DESIGN FOUND — see OPTIMIZATION_STRATEGY.md §7.1")
+
+    print("=" * 70)
+
+
+# ---------------------------------------------------------------------------
+# CLI
+# ---------------------------------------------------------------------------
+def parse_args() -> argparse.Namespace:
+    """Parse command-line arguments."""
+    parser = argparse.ArgumentParser(
+        description="Hydrotech Beam — Phase 1 LHS DoE Study",
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+        epilog=(
+            "Examples:\n"
+            "  # Development (stub solver):\n"
+            "  python run_doe.py --backend stub\n\n"
+            "  # Real NX evaluation:\n"
+            "  python run_doe.py --backend nxopen --model-dir /path/to/models\n\n"
+            "  # Resume interrupted study:\n"
+            "  python run_doe.py --backend stub --resume\n"
+        ),
+    )
+    parser.add_argument(
+        "--backend",
+        choices=["stub", "nxopen"],
+        default="stub",
+        help="NX solver backend: 'stub' for testing, 'nxopen' for real (default: stub)",
+    )
+    parser.add_argument(
+        "--model-dir",
+        default="",
+        help="Path to NX model directory (required for --backend nxopen)",
+    )
+    parser.add_argument(
+        "--study-name",
+        default=DEFAULT_STUDY_NAME,
+        help=f"Optuna study name (default: {DEFAULT_STUDY_NAME})",
+    )
+    parser.add_argument(
+        "--n-samples",
+        type=int,
+        default=DEFAULT_N_SAMPLES,
+        help=f"Number of LHS sample points (default: {DEFAULT_N_SAMPLES})",
+    )
+    parser.add_argument(
+        "--seed",
+        type=int,
+        default=DEFAULT_SEED,
+        help=f"Random seed for LHS (default: {DEFAULT_SEED})",
+    )
+    parser.add_argument(
+        "--results-dir",
+        default=DEFAULT_RESULTS_DIR,
+        help=f"Output directory for results (default: {DEFAULT_RESULTS_DIR})",
+    )
+    parser.add_argument(
+        "--resume",
+        action="store_true",
+        help="Resume an existing study instead of creating a new one",
+    )
+    parser.add_argument(
+        "--verbose", "-v",
+        action="store_true",
+        help="Enable verbose (DEBUG) logging",
+    )
+
+    return parser.parse_args()
+
+
+def main() -> None:
+    """Entry point."""
+    args = parse_args()
+
+    # Configure logging
+    log_level = logging.DEBUG if args.verbose else logging.INFO
+    logging.basicConfig(
+        level=log_level,
+        format="%(asctime)s [%(levelname)-7s] %(name)s: %(message)s",
+        datefmt="%Y-%m-%d %H:%M:%S",
+    )
+
+    # Run
+    try:
+        run_study(args)
+    except KeyboardInterrupt:
+        logger.warning("Study interrupted by user. Progress saved to Optuna DB.")
+        logger.info("Resume with: python run_doe.py --resume --study-name %s", args.study_name)
+        sys.exit(1)
+    except Exception:
+        logger.exception("Study failed with unexpected error")
+        sys.exit(2)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/projects/hydrotech-beam/studies/01_doe_landscape/sampling.py b/projects/hydrotech-beam/studies/01_doe_landscape/sampling.py
new file mode 100644
index 00000000..6e88ca86
--- /dev/null
+++ b/projects/hydrotech-beam/studies/01_doe_landscape/sampling.py
@@ -0,0 +1,264 @@
+"""Latin Hypercube Sampling for Hydrotech Beam DoE.
+
+Generates LHS sample points for the 4 design variables with:
+- Maximin LHS for space-filling coverage
+- Integer rounding for hole_count (DV4)
+- Stratified integer sampling to ensure all 11 hole_count levels are covered
+- Baseline (Trial 0) always included
+
+References:
+    OPTIMIZATION_STRATEGY.md §2.3 — Phase 1 configuration
+    OPTIMIZATION_STRATEGY.md §1.4 — Integer handling
+"""
+
+from __future__ import annotations
+
+import logging
+from dataclasses import dataclass
+from typing import Any
+
+import numpy as np
+from scipy.stats.qmc import LatinHypercube
+
+from geometric_checks import DesignPoint
+
+logger = logging.getLogger(__name__)
+
+
+# ---------------------------------------------------------------------------
+# Design variable bounds
+# ---------------------------------------------------------------------------
+@dataclass(frozen=True)
+class DVBounds:
+    """Design variable bounds and metadata."""
+
+    name: str
+    nx_expression: str
+    lower: float
+    upper: float
+    baseline: float
+    is_integer: bool = False
+
+
+DV_DEFINITIONS: list[DVBounds] = [
+    DVBounds("beam_half_core_thickness", "beam_half_core_thickness", 10.0, 40.0, 25.162),
+    DVBounds("beam_face_thickness", "beam_face_thickness", 10.0, 40.0, 21.504),
+    DVBounds("holes_diameter", "holes_diameter", 150.0, 450.0, 300.0),
+    DVBounds("hole_count", "hole_count", 5.0, 15.0, 10.0, is_integer=True),
+]
+
+N_DVS = len(DV_DEFINITIONS)
+BASELINE_VALUES = [dv.baseline for dv in DV_DEFINITIONS]
+
+
+def get_baseline_point() -> DesignPoint:
+    """Return the baseline design point (Trial 0).
+
+    Returns:
+        DesignPoint with confirmed baseline values from NX introspection.
+    """
+    return DesignPoint(
+        beam_half_core_thickness=25.162,
+        beam_face_thickness=21.504,
+        holes_diameter=300.0,
+        hole_count=10,
+    )
+
+
+def generate_lhs_samples(
+    n_samples: int = 50,
+    seed: int = 42,
+    include_baseline: bool = True,
+) -> list[DesignPoint]:
+    """Generate LHS sample points for the DoE study.
+
+    Strategy for integer coverage (hole_count = 5..15, 11 levels):
+        1. Generate `n_samples` LHS points with continuous DV4
+        2. Round DV4 to nearest integer
+        3. Check coverage — if any of the 11 levels are missing, replace
+           the closest duplicate trial with the missing level
+        This ensures all integer levels are represented while maintaining
+        the space-filling property of LHS for the continuous variables.
+
+    Args:
+        n_samples: Number of LHS sample points (default: 50).
+        seed: Random seed for reproducibility.
+        include_baseline: If True, prepend baseline as Trial 0.
+
+    Returns:
+        List of DesignPoint instances. If include_baseline is True,
+        the first element is the baseline (Trial 0).
+    """
+    logger.info(
+        "Generating %d LHS samples (seed=%d, baseline=%s)",
+        n_samples, seed, include_baseline,
+    )
+
+    # Generate unit hypercube LHS samples
+    sampler = LatinHypercube(d=N_DVS, seed=seed, optimization="random-cd")
+    unit_samples = sampler.random(n=n_samples)  # shape: (n_samples, 4)
+
+    # Scale to design variable bounds
+    samples = _scale_samples(unit_samples)
+
+    # Round hole_count to nearest integer
+    samples[:, 3] = np.round(samples[:, 3]).astype(int)
+    samples[:, 3] = np.clip(samples[:, 3], 5, 15)
+
+    # Ensure full integer coverage for hole_count
+    samples = _ensure_integer_coverage(samples, rng=np.random.default_rng(seed))
+
+    # Convert to DesignPoint list
+    points: list[DesignPoint] = []
+    if include_baseline:
+        points.append(get_baseline_point())
+
+    for row in samples:
+        points.append(
+            DesignPoint(
+                beam_half_core_thickness=float(row[0]),
+                beam_face_thickness=float(row[1]),
+                holes_diameter=float(row[2]),
+                hole_count=int(row[3]),
+            )
+        )
+
+    logger.info(
+        "Generated %d total points (%d LHS + %s baseline)",
+        len(points),
+        n_samples,
+        "1" if include_baseline else "0",
+    )
+    _log_coverage(points)
+
+    return points
+
+
+def _scale_samples(unit_samples: np.ndarray) -> np.ndarray:
+    """Scale unit hypercube [0,1]^d samples to design variable bounds.
+
+    Args:
+        unit_samples: Array of shape (n, 4) with values in [0, 1].
+
+    Returns:
+        Scaled array with values in [lower, upper] for each DV.
+    """
+    lower = np.array([dv.lower for dv in DV_DEFINITIONS])
+    upper = np.array([dv.upper for dv in DV_DEFINITIONS])
+    return lower + unit_samples * (upper - lower)
+
+
+def _ensure_integer_coverage(
+    samples: np.ndarray,
+    rng: np.random.Generator,
+) -> np.ndarray:
+    """Ensure all 11 hole_count levels (5-15) are represented.
+
+    If any integer level is missing, replace a duplicate from the most
+    over-represented level with a sample at the missing level.
+    Continuous DVs for replacement samples are drawn randomly within bounds.
+
+    Args:
+        samples: Array of shape (n, 4) with rounded hole_count in col 3.
+        rng: NumPy random generator for reproducibility.
+
+    Returns:
+        Modified samples array with full integer coverage.
+    """
+    all_levels = set(range(5, 16))  # {5, 6, 7, ..., 15}
+    present_levels = set(int(x) for x in samples[:, 3])
+    missing_levels = all_levels - present_levels
+
+    if not missing_levels:
+        logger.info("All 11 hole_count levels represented ✓")
+        return samples
+
+    logger.warning(
+        "Missing hole_count levels: %s — patching with replacements",
+        sorted(missing_levels),
+    )
+
+    for missing_level in sorted(missing_levels):
+        # Find the most over-represented level
+        unique, counts = np.unique(samples[:, 3].astype(int), return_counts=True)
+        most_common_idx = np.argmax(counts)
+        most_common_level = unique[most_common_idx]
+
+        # Find indices with the most common level
+        candidates = np.where(samples[:, 3].astype(int) == most_common_level)[0]
+        replace_idx = rng.choice(candidates)
+
+        # Replace: keep continuous DVs random within bounds, set hole_count
+        for j, dv in enumerate(DV_DEFINITIONS):
+            if dv.is_integer:
+                samples[replace_idx, j] = missing_level
+            else:
+                samples[replace_idx, j] = rng.uniform(dv.lower, dv.upper)
+
+        logger.info(
+            "  Replaced trial at idx %d: hole_count %d → %d",
+            replace_idx, most_common_level, missing_level,
+        )
+
+    return samples
+
+
+def _log_coverage(points: list[DesignPoint]) -> None:
+    """Log hole_count coverage statistics."""
+    counts: dict[int, int] = {}
+    for p in points:
+        counts[p.hole_count] = counts.get(p.hole_count, 0) + 1
+
+    logger.info("Hole count coverage:")
+    for level in range(5, 16):
+        n = counts.get(level, 0)
+        logger.info("  hole_count=%2d: %d trials", level, n)
+
+
+def points_to_dicts(points: list[DesignPoint]) -> list[dict[str, Any]]:
+    """Convert DesignPoint list to list of dicts (for Optuna enqueue).
+
+    Args:
+        points: List of DesignPoint instances.
+
+    Returns:
+        List of dicts with DV names as keys.
+    """
+    return [
+        {
+            "beam_half_core_thickness": p.beam_half_core_thickness,
+            "beam_face_thickness": p.beam_face_thickness,
+            "holes_diameter": p.holes_diameter,
+            "hole_count": p.hole_count,
+        }
+        for p in points
+    ]
+
+
+# ---------------------------------------------------------------------------
+# Quick self-test
+# ---------------------------------------------------------------------------
+if __name__ == "__main__":
+    logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
+
+    points = generate_lhs_samples(n_samples=50, seed=42)
+    print(f"\nGenerated {len(points)} total points")
+    print(f"  Trial 0 (baseline): {points[0]}")
+    print(f"  Trial 1 (first LHS): {points[1]}")
+    print(f"  Trial {len(points)-1} (last LHS): {points[-1]}")
+
+    # Verify coverage
+    hole_counts = {p.hole_count for p in points}
+    expected = set(range(5, 16))
+    assert hole_counts == expected, (
+        f"Missing hole_count levels: {expected - hole_counts}"
+    )
+    print("\nAll 11 hole_count levels covered ✓")
+
+    # Verify bounds
+    for p in points[1:]:  # skip baseline
+        assert 10.0 <= p.beam_half_core_thickness <= 40.0
+        assert 10.0 <= p.beam_face_thickness <= 40.0
+        assert 150.0 <= p.holes_diameter <= 450.0
+        assert 5 <= p.hole_count <= 15
+    print("All samples within bounds ✓")