Atomizer/.claude/skills/modules/DYNAMIC_RESPONSE_MASTER_PLAN.md

Dynamic Response Processor - Master Implementation Plan

Atomizer Feature: Modal Superposition Engine for Dynamic Optimization

Document Version: 1.0
Created: 2025-12-22
Author: Antoine + Claude
Status: Implementation Ready
Priority: High - Differentiating Feature

---

Executive Summary

This document provides the complete implementation roadmap for adding analytical dynamic response processing to Atomizer. This feature enables optimization of structures under random vibration, sine sweep, and shock environments — without running expensive frequency response FEA solutions for each design iteration.

Value Proposition

Traditional Workflow	Atomizer Dynamic Workflow
Run SOL 103 (modal)	Run SOL 103 (modal)
Run SOL 111 (frequency response)	Skip - compute analytically
Run SOL 112 (random)	Skip - compute analytically
~15-30 min per iteration	~30 sec per iteration
100 iterations = 25-50 hours	100 iterations = 50 minutes

Key Differentiator: SATK is a post-processor. Atomizer becomes an optimization engine with built-in dynamic response intelligence.

---

Table of Contents

1. Technical Foundation
2. Architecture Design
3. Implementation Phases
4. Module Specifications
5. Integration Points
6. Validation Strategy
7. Performance Targets
8. Risk Mitigation
9. Future Extensions

---

1. Technical Foundation

1.1 The Physics: Why This Works

The fundamental insight is that dynamic response can be computed from modal data alone.

For a structure with N degrees of freedom, the equation of motion is:

[M]{ü} + [C]{u̇} + [K]{u} = {F(t)}

Modal analysis (SOL 103) decomposes this into N independent single-DOF equations:

q̈ᵢ + 2ζᵢωᵢq̇ᵢ + ωᵢ²qᵢ = Lᵢ·a_base / mᵢ

Where:

• qᵢ = modal coordinate (how much mode i is excited)
• ωᵢ = natural frequency of mode i (rad/s)
• ζᵢ = damping ratio of mode i
• Lᵢ = modal participation factor (how mode i couples to base motion)
• mᵢ = generalized (modal) mass

The key: Once we extract (ωᵢ, ζᵢ, Lᵢ, φᵢ) from SOL 103, we can compute the response to ANY excitation profile analytically — no more FEA needed.

1.2 Transfer Function Mathematics

Single Mode Transfer Function

For base acceleration input, the displacement transfer function of mode i is:

def H_displacement(f, f_n, zeta):
    """
    Displacement per unit base acceleration.
    
    H(f) = 1 / [(ωₙ² - ω²) + j·2·ζ·ωₙ·ω]
    
    In dimensionless form with r = f/fₙ:
    H(r) = 1 / [(1 - r²) + j·2·ζ·r]
    """
    r = f / f_n
    return 1 / ((1 - r**2) + 2j * zeta * r)

Acceleration Transmissibility

For base excitation, the ratio of response acceleration to input acceleration:

def transmissibility(f, f_n, zeta):
    """
    T(f) = (1 + j·2·ζ·r) / [(1 - r²) + j·2·ζ·r]
    
    |T| = 1 at low frequency (rigid body)
    |T| = Q at resonance (amplification)
    |T| → 0 at high frequency (isolation)
    """
    r = f / f_n
    return (1 + 2j * zeta * r) / ((1 - r**2) + 2j * zeta * r)

Multi-Mode Combination

Total response sums all modal contributions:

def total_response(frequencies, modes, participation_factors):
    """
    H_total(f) = Σᵢ Lᵢ · φᵢ · Hᵢ(f) / mᵢ
    
    Where:
    - Lᵢ = participation factor
    - φᵢ = mode shape at point of interest
    - Hᵢ = modal transfer function
    - mᵢ = generalized mass
    """
    H_total = np.zeros(len(frequencies), dtype=complex)
    for mode in modes:
        H_mode = H_displacement(frequencies, mode.frequency, mode.damping)
        H_total += mode.participation * mode.shape * H_mode / mode.gen_mass
    return H_total

1.3 Random Vibration Response

Random vibration is characterized by Power Spectral Density (PSD), typically in G²/Hz.

Response PSD

S_response(f) = |H(f)|² × S_input(f)

RMS (Root Mean Square)

σ² = ∫ S_response(f) df  # Mean square
RMS = √σ²                 # Root mean square

Peak Response

For Gaussian random processes:

Peak = k × RMS

Where k (crest factor) depends on:

• k = 3.0: Standard 3-sigma (99.7% probability not exceeded)
• k = 3.5-4.0: Von Mises stress (non-Gaussian distribution)
• k = 4.0-4.5: Fatigue-critical applications

Miles' Equation (SDOF Approximation)

For quick estimates with white noise PSD:

Grms = √(π/2 × fₙ × Q × PSD_level)

Where Q = 1/(2ζ) is the quality factor.

1.4 Stress Recovery from Modal Data

Element stresses can also be computed via modal superposition:

σ(t) = Σᵢ qᵢ(t) × σᵢ_modal

Where σᵢ_modal is the stress pattern when mode i has unit amplitude.

For random vibration, stress PSD:

S_σ(f) = |Σᵢ Lᵢ × σᵢ_modal × Hᵢ(f)|² × S_input(f)

1.5 What We Extract from SOL 103

Data	Source	Used For
Natural frequencies (fₙ)	OP2: eigenvalues	Transfer function poles
Mode shapes (φ)	OP2: eigenvectors	Response at specific nodes
Generalized mass (m)	OP2/F06: eigenvalues	Modal scaling
Modal effective mass	F06: MEFFMASS	Participation estimation
Participation factors (L)	F06: MEFFMASS	Base excitation coupling
Modal stress shapes	OP2: modal stress	Stress recovery

---

2. Architecture Design

2.1 Module Structure

optimization_engine/
├── processors/
│   └── dynamic_response/
│       ├── __init__.py                 # Public API exports
│       ├── modal_database.py           # ModalDatabase, ModeShape classes
│       ├── transfer_functions.py       # FRF computation engine
│       ├── random_vibration.py         # Random response processor
│       ├── sine_sweep.py               # Harmonic response processor
│       ├── shock_response.py           # SRS computation (Phase 5)
│       ├── stress_recovery.py          # Modal stress combination
│       ├── psd_profiles.py             # Standard PSD definitions
│       └── utils/
│           ├── __init__.py
│           ├── frequency_utils.py      # Frequency array generation
│           ├── integration.py          # Spectral integration methods
│           └── statistics.py           # Peak factors, distributions
│
├── extractors/
│   ├── extract_modal_database.py       # Full modal data extraction
│   ├── extract_random_response.py      # Random vibration objectives
│   ├── extract_sine_response.py        # Sine sweep objectives
│   └── extract_dynamic_stress.py       # Dynamic stress extraction
│
├── hooks/
│   └── dynamic_response/
│       ├── pre_dynamic_analysis.py     # Validate modal data
│       └── post_dynamic_analysis.py    # Cache results, logging
│
└── templates/
    ├── random_vibration_optimization.json
    ├── sine_sweep_optimization.json
    └── multi_environment_optimization.json

2.2 Class Hierarchy

┌─────────────────────────────────────────────────────────────┐
│                    DynamicResponseProcessor                  │
│                        (Abstract Base)                       │
├─────────────────────────────────────────────────────────────┤
│ + modal_db: ModalDatabase                                   │
│ + frf_engine: TransferFunctionEngine                        │
│ + compute_response() → ResponseResult                       │
│ + compute_at_node(node_id) → float                         │
│ + compute_stress(element_ids) → Dict                        │
└─────────────────────────────────────────────────────────────┘
                              │
          ┌───────────────────┼───────────────────┐
          ▼                   ▼                   ▼
┌─────────────────┐  ┌─────────────────┐  ┌─────────────────┐
│ RandomVibration │  │   SineSweep     │  │ ShockResponse   │
│   Processor     │  │   Processor     │  │   Processor     │
├─────────────────┤  ├─────────────────┤  ├─────────────────┤
│ + psd_profile   │  │ + sweep_rate    │  │ + srs_damping   │
│ + compute_grms()│  │ + dwell_freqs   │  │ + compute_srs() │
│ + compute_peak()│  │ + compute_peak()│  │ + compute_pseudo│
└─────────────────┘  └─────────────────┘  └─────────────────┘

2.3 Data Flow

┌──────────────────────────────────────────────────────────────────┐
│                      OPTIMIZATION LOOP                            │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  1. UPDATE DESIGN PARAMETERS                                      │
│     NX Expression Manager → thickness=2.5mm, rib_height=12mm     │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  2. RUN SOL 103 (Modal Analysis)                                  │
│     NX Solver → model.op2, model.f06                             │
│     Time: ~20-60 seconds                                          │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  3. EXTRACT MODAL DATABASE                                        │
│     pyNastran → ModalDatabase object                              │
│     - Frequencies, mode shapes, participation factors             │
│     - Modal stress shapes (if needed)                             │
│     Time: ~0.5-2 seconds                                          │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  4. COMPUTE DYNAMIC RESPONSE (Analytical)                         │
│     RandomVibrationProcessor.compute_response()                   │
│     - Build transfer functions: ~5ms                              │
│     - Compute response PSD: ~10ms                                 │
│     - Integrate for RMS: ~5ms                                     │
│     - Stress recovery: ~50ms                                      │
│     Time: ~50-100 milliseconds                                    │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  5. RETURN OBJECTIVES & CONSTRAINTS                               │
│     {                                                             │
│       'grms': 8.4,                                               │
│       'peak_acceleration': 25.2,                                  │
│       'peak_stress': 145.6,                                       │
│       'first_frequency': 127.3                                    │
│     }                                                             │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  6. OPTUNA SUGGESTS NEXT DESIGN                                   │
│     TPE/CMA-ES/NSGA-II → next trial parameters                   │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
                        [REPEAT 1-6]

2.4 Memory Management

For large models with many modes/elements, we need smart memory handling:

class ModalDatabase:
    """
    Memory-efficient modal data storage.
    
    Strategies:
    1. Lazy loading - only load mode shapes when needed
    2. Sparse storage - only store non-zero eigenvector components
    3. Node filtering - only extract nodes of interest
    4. Disk caching - cache to HDF5 for reuse across trials
    """
    
    def __init__(self, op2_file, lazy_load=True):
        self._op2_path = op2_file
        self._modes_loaded = False
        self._mode_shapes = None  # Load on demand
        
    @property
    def mode_shapes(self):
        if not self._modes_loaded:
            self._load_mode_shapes()
        return self._mode_shapes
    
    def get_modes_for_nodes(self, node_ids: List[int]) -> np.ndarray:
        """Extract mode shapes only for specified nodes."""
        # Efficient partial loading
        pass
    
    def cache_to_hdf5(self, cache_path: Path):
        """Cache modal database for fast reload."""
        pass
    
    @classmethod
    def from_cache(cls, cache_path: Path) -> 'ModalDatabase':
        """Load from HDF5 cache (10x faster than OP2)."""
        pass

---

3. Implementation Phases

Phase 1: Core Infrastructure (Week 1)

Goal: Build the foundational classes and basic functionality.

Tasks

Task	Description	Effort	Deliverable
1.1	Create ModeShape dataclass	2h	modal_database.py
1.2	Create ModalDatabase class	4h	modal_database.py
1.3	Implement OP2 extraction	4h	ModalDatabase.from_op2_f06()
1.4	Implement F06 MEFFMASS parsing	3h	Modal effective mass extraction
1.5	Create TransferFunctionEngine	4h	transfer_functions.py
1.6	Single mode FRF	2h	single_mode_frf()
1.7	Multi-mode combination	2h	base_excitation_frf()
1.8	Unit tests for FRF	3h	tests/test_transfer_functions.py

Validation Checkpoint:

• FRF matches closed-form SDOF solution
• Modal database correctly extracts from sample OP2

Code: ModeShape Dataclass

# optimization_engine/processors/dynamic_response/modal_database.py

from dataclasses import dataclass, field
from typing import Dict, Optional, List
import numpy as np

@dataclass
class ModeShape:
    """
    Single mode shape data container.
    
    Stores all information needed for modal superposition
    of one natural mode.
    """
    # Identification
    mode_number: int
    
    # Dynamic properties
    frequency: float  # Hz
    damping_ratio: float = 0.02  # Dimensionless (2% default)
    
    # Modal mass/stiffness
    generalized_mass: float = 1.0  # kg (normalized)
    generalized_stiffness: float = 0.0  # N/m
    
    # Participation factors (for base excitation)
    participation_x: float = 0.0
    participation_y: float = 0.0
    participation_z: float = 0.0
    
    # Modal effective mass (absolute contribution)
    effective_mass_x: float = 0.0
    effective_mass_y: float = 0.0
    effective_mass_z: float = 0.0
    
    # Eigenvector at specific nodes
    # Dict[node_id, np.array([tx, ty, tz, rx, ry, rz])]
    eigenvector: Dict[int, np.ndarray] = field(default_factory=dict)
    
    # Modal stress shapes for stress recovery
    # Dict[element_id, np.array([sxx, syy, szz, sxy, syz, sxz, vm])]
    stress_shapes: Dict[int, np.ndarray] = field(default_factory=dict)
    
    # Derived properties
    @property
    def omega(self) -> float:
        """Angular frequency (rad/s)."""
        return 2 * np.pi * self.frequency
    
    @property
    def omega_squared(self) -> float:
        """ω² for efficiency."""
        return self.omega ** 2
    
    @property
    def quality_factor(self) -> float:
        """Q = 1/(2ζ) - amplification at resonance."""
        if self.damping_ratio > 0:
            return 1 / (2 * self.damping_ratio)
        return np.inf
    
    @property
    def half_power_bandwidth(self) -> float:
        """Δf = fₙ/Q - bandwidth at -3dB."""
        return self.frequency / self.quality_factor
    
    def get_participation(self, direction: str) -> float:
        """Get participation factor for direction ('x', 'y', 'z')."""
        return getattr(self, f'participation_{direction}', 0.0)
    
    def get_effective_mass(self, direction: str) -> float:
        """Get effective mass for direction."""
        return getattr(self, f'effective_mass_{direction}', 0.0)
    
    def get_eigenvector_component(
        self, 
        node_id: int, 
        direction: str
    ) -> float:
        """Get eigenvector component at node in direction."""
        if node_id not in self.eigenvector:
            return 0.0
        dir_map = {'x': 0, 'y': 1, 'z': 2, 'rx': 3, 'ry': 4, 'rz': 5}
        idx = dir_map.get(direction, 0)
        return self.eigenvector[node_id][idx]

Phase 2: Random Vibration Processor (Week 2)

Goal: Full random vibration response calculation with stress recovery.

Tasks

Task	Description	Effort	Deliverable
2.1	Create PSDProfile class	2h	psd_profiles.py
2.2	Standard profiles (GEVS, MIL-STD)	2h	STANDARD_PROFILES dict
2.3	PSD interpolation (log-log)	2h	PSDProfile.interpolate()
2.4	Create RandomVibrationProcessor	4h	random_vibration.py
2.5	Response PSD calculation	3h	compute_response_psd()
2.6	Spectral integration	2h	compute_rms()
2.7	Peak factor calculation	2h	Gaussian + non-Gaussian
2.8	Modal contribution analysis	2h	Identify dominant modes
2.9	Stress recovery	4h	compute_stress_response()
2.10	Result dataclass	2h	RandomVibrationResult
2.11	Integration tests	4h	Miles' equation validation

Validation Checkpoint:

• Grms matches Miles' equation for SDOF
• Multi-mode response reasonable for sample model
• Stress recovery gives correct distribution

Code: PSD Profiles

# optimization_engine/processors/dynamic_response/psd_profiles.py

import numpy as np
from scipy import interpolate
from dataclasses import dataclass
from typing import List, Tuple, Optional
from pathlib import Path


@dataclass
class PSDProfile:
    """
    Power Spectral Density profile for random vibration.
    
    Stores PSD as G²/Hz vs frequency, with log-log interpolation.
    """
    name: str
    frequencies: np.ndarray  # Hz
    psd_values: np.ndarray   # G²/Hz
    description: str = ""
    source: str = ""  # Standard reference
    
    @property
    def freq_min(self) -> float:
        return float(self.frequencies.min())
    
    @property
    def freq_max(self) -> float:
        return float(self.frequencies.max())
    
    @property
    def grms(self) -> float:
        """Overall Grms level."""
        from scipy.integrate import trapezoid
        return np.sqrt(trapezoid(self.psd_values, self.frequencies))
    
    def interpolate(self, frequencies: np.ndarray) -> np.ndarray:
        """
        Interpolate PSD to arbitrary frequencies.
        
        Uses log-log interpolation (standard for PSD).
        Clamps to edge values outside defined range.
        """
        # Safety: ensure positive frequencies
        frequencies = np.maximum(frequencies, 1e-6)
        
        log_f = np.log10(self.frequencies)
        log_psd = np.log10(np.maximum(self.psd_values, 1e-20))
        
        interp_func = interpolate.interp1d(
            log_f, log_psd,
            kind='linear',
            bounds_error=False,
            fill_value=(log_psd[0], log_psd[-1])
        )
        
        return 10 ** interp_func(np.log10(frequencies))
    
    def plot(self, ax=None, **kwargs):
        """Plot PSD profile."""
        import matplotlib.pyplot as plt
        if ax is None:
            fig, ax = plt.subplots()
        ax.loglog(self.frequencies, self.psd_values, **kwargs)
        ax.set_xlabel('Frequency (Hz)')
        ax.set_ylabel('PSD (G²/Hz)')
        ax.set_title(f'{self.name} - {self.grms:.2f} Grms')
        ax.grid(True, which='both', alpha=0.3)
        return ax
    
    @classmethod
    def from_breakpoints(
        cls, 
        breakpoints: List[Tuple[float, float]],
        name: str = "Custom",
        description: str = ""
    ) -> 'PSDProfile':
        """
        Create PSD from breakpoint specification.
        
        Args:
            breakpoints: List of (frequency_hz, psd_g2_hz) tuples
            name: Profile name
            description: Description text
        
        Example:
            psd = PSDProfile.from_breakpoints([
                (20, 0.01),    # 20 Hz: 0.01 G²/Hz
                (80, 0.04),    # 80 Hz: 0.04 G²/Hz (ramp up)
                (500, 0.04),   # 500 Hz: 0.04 G²/Hz (flat)
                (2000, 0.007)  # 2000 Hz: 0.007 G²/Hz (ramp down)
            ], name='MIL-STD-810')
        """
        breakpoints = sorted(breakpoints, key=lambda x: x[0])
        frequencies = np.array([bp[0] for bp in breakpoints])
        psd_values = np.array([bp[1] for bp in breakpoints])
        return cls(
            name=name,
            frequencies=frequencies,
            psd_values=psd_values,
            description=description
        )
    
    @classmethod
    def from_csv(cls, filepath: Path, name: str = None) -> 'PSDProfile':
        """Load PSD from CSV file with freq,psd columns."""
        data = np.loadtxt(filepath, delimiter=',', skiprows=1)
        return cls(
            name=name or filepath.stem,
            frequencies=data[:, 0],
            psd_values=data[:, 1]
        )
    
    @classmethod
    def flat(
        cls,
        level: float,
        freq_min: float = 20,
        freq_max: float = 2000,
        name: str = "Flat"
    ) -> 'PSDProfile':
        """Create flat (white noise) PSD profile."""
        return cls(
            name=name,
            frequencies=np.array([freq_min, freq_max]),
            psd_values=np.array([level, level]),
            description=f"Flat PSD at {level} G²/Hz"
        )


# ============================================================================
# STANDARD PSD PROFILES
# ============================================================================

STANDARD_PROFILES = {
    # NASA GEVS (General Environmental Verification Standard)
    'GEVS': PSDProfile.from_breakpoints([
        (20, 0.026),
        (50, 0.16),
        (800, 0.16),
        (2000, 0.026),
    ], name='NASA GEVS', 
       description='NASA General Environmental Verification Standard'),
    
    # MIL-STD-810G Category 1 (Basic)
    'MIL-STD-810G_CAT1': PSDProfile.from_breakpoints([
        (10, 0.04),
        (40, 0.04),
        (500, 0.04),
        (2000, 0.01),
    ], name='MIL-STD-810G Cat 1',
       description='Basic transportation vibration'),
    
    # MIL-STD-810G Category 4 (Truck)
    'MIL-STD-810G_CAT4': PSDProfile.from_breakpoints([
        (5, 0.001),
        (10, 0.015),
        (40, 0.015),
        (500, 0.015),
        (2000, 0.003),
    ], name='MIL-STD-810G Cat 4',
       description='Truck transportation'),
    
    # MIL-STD-810G Category 7 (Jet Aircraft)
    'MIL-STD-810G_CAT7': PSDProfile.from_breakpoints([
        (15, 0.015),
        (80, 0.015),
        (350, 0.015),
        (2000, 0.006),
    ], name='MIL-STD-810G Cat 7',
       description='Jet aircraft vibration'),
    
    # MIL-STD-810G Category 24 (Helicopter)
    'MIL-STD-810G_CAT24': PSDProfile.from_breakpoints([
        (5, 0.001),
        (20, 0.01),
        (80, 0.04),
        (350, 0.04),
        (2000, 0.007),
    ], name='MIL-STD-810G Cat 24',
       description='Helicopter vibration'),
    
    # NAVMAT P-9492
    'NAVMAT': PSDProfile.from_breakpoints([
        (20, 0.04),
        (80, 0.04),
        (350, 0.04),
        (2000, 0.007),
    ], name='NAVMAT P-9492',
       description='Navy shipboard equipment'),
    
    # ESA PSS-01-401 Qualification
    'ESA_QUAL': PSDProfile.from_breakpoints([
        (20, 0.01),
        (100, 0.08),
        (300, 0.08),
        (2000, 0.008),
    ], name='ESA PSS-01-401',
       description='ESA qualification level'),
    
    # ESA PSS-01-401 Acceptance
    'ESA_ACCEPT': PSDProfile.from_breakpoints([
        (20, 0.005),
        (100, 0.04),
        (300, 0.04),
        (2000, 0.004),
    ], name='ESA PSS-01-401 Accept',
       description='ESA acceptance level'),
    
    # SpaceX Falcon 9 (approximate public data)
    'FALCON9': PSDProfile.from_breakpoints([
        (20, 0.01),
        (50, 0.08),
        (800, 0.08),
        (2000, 0.02),
    ], name='Falcon 9 (approx)',
       description='SpaceX Falcon 9 approximate envelope'),
    
    # Ariane 5 (approximate)
    'ARIANE5': PSDProfile.from_breakpoints([
        (20, 0.015),
        (100, 0.09),
        (400, 0.09),
        (2000, 0.01),
    ], name='Ariane 5 (approx)',
       description='Ariane 5 approximate envelope'),
}


def get_psd_profile(name: str) -> PSDProfile:
    """
    Get a standard PSD profile by name.
    
    Args:
        name: Profile name (case-insensitive, underscores/hyphens flexible)
    
    Returns:
        PSDProfile instance
    
    Raises:
        ValueError: If profile not found
    
    Example:
        >>> psd = get_psd_profile('GEVS')
        >>> psd = get_psd_profile('mil-std-810g-cat7')
    """
    # Normalize name
    normalized = name.upper().replace('-', '_').replace(' ', '_')
    
    if normalized in STANDARD_PROFILES:
        return STANDARD_PROFILES[normalized]
    
    # Try partial match
    for key, profile in STANDARD_PROFILES.items():
        if normalized in key or key in normalized:
            return profile
    
    available = ', '.join(STANDARD_PROFILES.keys())
    raise ValueError(f"Unknown PSD profile: '{name}'. Available: {available}")


def list_available_profiles() -> List[str]:
    """Return list of available standard PSD profile names."""
    return list(STANDARD_PROFILES.keys())

Phase 3: Atomizer Integration (Week 3)

Goal: Connect dynamic response processor to Atomizer's extractor and optimization systems.

Tasks

Task	Description	Effort	Deliverable
3.1	Create extract_modal_database	3h	New extractor
3.2	Create extract_random_response	3h	New extractor
3.3	Create extract_dynamic_stress	3h	New extractor
3.4	Optimization template	2h	random_vibration_optimization.json
3.5	Runner integration	4h	Hook into OptimizationRunner
3.6	Dashboard widgets	4h	Response PSD plot, Grms history
3.7	Result caching	3h	Cache modal DB between similar designs
3.8	Documentation	3h	User guide, API docs
3.9	End-to-end test	4h	Full optimization run

Validation Checkpoint:

• Complete optimization runs successfully
• Dashboard shows dynamic response metrics
• Results match manual SATK/Nastran verification

Code: Random Response Extractor

# optimization_engine/extractors/extract_random_response.py

"""
Random Vibration Response Extractor for Atomizer Optimization
=============================================================

Extracts random vibration metrics (Grms, peak acceleration, peak stress)
from SOL 103 modal analysis results using analytical modal superposition.

This extractor enables dynamic optimization without running SOL 111.

Example:
    >>> result = extract_random_response(
    ...     op2_file='model.op2',
    ...     psd_profile='GEVS',
    ...     direction='z'
    ... )
    >>> print(f"Grms: {result['grms']:.2f} G")
    >>> print(f"Peak: {result['peak_acceleration']:.2f} G")
"""

from pathlib import Path
from typing import Dict, Any, Optional, Union, List
import logging

logger = logging.getLogger(__name__)


def extract_random_response(
    op2_file: Union[str, Path],
    psd_profile: str = 'GEVS',
    direction: str = 'z',
    damping: float = 0.02,
    node_id: Optional[int] = None,
    peak_sigma: float = 3.0,
    f06_file: Optional[Union[str, Path]] = None,
    n_modes: Optional[int] = None,
    frequency_resolution: int = 500,
) -> Dict[str, Any]:
    """
    Extract random vibration response for optimization objective/constraint.
    
    This is the main interface for dynamic optimization. It reads SOL 103
    results and computes random vibration response analytically.
    
    Args:
        op2_file: Path to OP2 file from SOL 103 analysis
        psd_profile: PSD specification - name ('GEVS', 'MIL-STD-810G_CAT7', etc.)
                     or path to CSV file
        direction: Excitation direction ('x', 'y', 'z')
        damping: Modal damping ratio (default 2%)
        node_id: Specific node for response (None = use effective mass)
        peak_sigma: Sigma level for peak (3.0 = 3-sigma)
        f06_file: Optional F06 for modal effective mass table
        n_modes: Number of modes to include (None = auto 95% mass)
        frequency_resolution: Points in frequency array
    
    Returns:
        Dict containing:
        - 'grms': RMS acceleration (G)
        - 'peak_acceleration': Peak acceleration (G)
        - 'rms_displacement': RMS displacement (mm)
        - 'peak_displacement': Peak displacement (mm)
        - 'dominant_mode': Mode number with highest contribution
        - 'dominant_frequency': Frequency of dominant mode (Hz)
        - 'n_modes_used': Number of modes in calculation
        - 'mass_fraction': Cumulative mass fraction of included modes
        - 'psd_profile': Name of PSD profile used
        - 'direction': Excitation direction
        - 'success': True if extraction succeeded
        - 'error': Error message if failed
    
    Example:
        # For optimization objective
        >>> result = extract_random_response('bracket.op2', psd_profile='GEVS')
        >>> objective = result['peak_acceleration']
        
        # For constraint checking
        >>> result = extract_random_response('bracket.op2', psd_profile='GEVS')
        >>> if result['peak_acceleration'] > 25.0:
        ...     print("Constraint violated!")
    """
    from optimization_engine.processors.dynamic_response import (
        ModalDatabase,
        RandomVibrationProcessor,
        get_psd_profile
    )
    
    op2_path = Path(op2_file)
    f06_path = Path(f06_file) if f06_file else None
    
    # Auto-detect F06 if not provided
    if f06_path is None:
        possible_f06 = op2_path.with_suffix('.f06')
        if possible_f06.exists():
            f06_path = possible_f06
            logger.info(f"Auto-detected F06: {f06_path}")
    
    try:
        # Build modal database
        logger.info(f"Extracting modal database from {op2_path}")
        modal_db = ModalDatabase.from_op2_f06(
            op2_file=op2_path,
            f06_file=f06_path,
            damping=damping,
            node_ids=[node_id] if node_id else None,
        )
        
        logger.info(f"Extracted {modal_db.n_modes} modes, "
                   f"freq range: {modal_db.frequency_range[0]:.1f}-"
                   f"{modal_db.frequency_range[1]:.1f} Hz")
        
        # Load PSD profile
        if Path(psd_profile).exists():
            from optimization_engine.processors.dynamic_response import PSDProfile
            psd = PSDProfile.from_csv(Path(psd_profile))
        else:
            psd = get_psd_profile(psd_profile)
        
        logger.info(f"Using PSD profile: {psd.name} ({psd.grms:.2f} Grms)")
        
        # Create processor
        processor = RandomVibrationProcessor(
            modal_db=modal_db,
            default_damping=damping
        )
        
        # Compute response
        result = processor.compute_response(
            psd_profile=psd,
            direction=direction,
            node_id=node_id,
            n_modes=n_modes,
            frequency_resolution=frequency_resolution,
            peak_sigma=peak_sigma,
            output_type='acceleration'
        )
        
        # Get mass fraction for included modes
        n_modes_used = n_modes or modal_db.modes_for_mass_fraction(0.95, direction)
        mass_fraction = modal_db.get_cumulative_mass_fraction(direction)[n_modes_used - 1]
        
        # Dominant mode info
        dominant_freq = modal_db.modes[result.dominant_mode - 1].frequency
        
        return {
            'grms': result.rms_acceleration,
            'peak_acceleration': result.peak_acceleration,
            'rms_displacement': result.rms_displacement,
            'peak_displacement': result.peak_displacement,
            'dominant_mode': result.dominant_mode,
            'dominant_frequency': dominant_freq,
            'n_modes_used': n_modes_used,
            'mass_fraction': mass_fraction,
            'modal_contributions': result.modal_contributions,
            'psd_profile': psd.name,
            'psd_grms': psd.grms,
            'direction': direction,
            'peak_sigma': peak_sigma,
            'success': True,
            'error': None,
        }
        
    except Exception as e:
        logger.exception(f"Random response extraction failed: {e}")
        return {
            'grms': None,
            'peak_acceleration': None,
            'rms_displacement': None,
            'peak_displacement': None,
            'dominant_mode': None,
            'dominant_frequency': None,
            'success': False,
            'error': str(e),
        }


def extract_random_stress(
    op2_file: Union[str, Path],
    psd_profile: str = 'GEVS',
    direction: str = 'z',
    damping: float = 0.02,
    peak_sigma: float = 3.0,
    element_ids: Optional[List[int]] = None,
    non_gaussian_correction: bool = True,
    f06_file: Optional[Union[str, Path]] = None,
) -> Dict[str, Any]:
    """
    Extract peak Von Mises stress from random vibration.
    
    Args:
        op2_file: Path to OP2 with modal stress results
        psd_profile: PSD specification
        direction: Excitation direction
        damping: Modal damping ratio
        peak_sigma: Sigma level
        element_ids: Specific elements (None = all)
        non_gaussian_correction: Apply VM stress distribution correction
        f06_file: Optional F06 for effective mass
    
    Returns:
        Dict with:
        - 'max_rms_von_mises': Maximum RMS VM stress (MPa)
        - 'max_peak_von_mises': Maximum peak VM stress (MPa)
        - 'critical_element': Element ID with max stress
        - 'margin_of_safety': (allowable - peak) / peak
        - 'n_elements': Number of elements analyzed
    """
    from optimization_engine.processors.dynamic_response import (
        ModalDatabase,
        RandomVibrationProcessor,
        get_psd_profile
    )
    
    op2_path = Path(op2_file)
    f06_path = Path(f06_file) if f06_file else None
    
    try:
        # Build modal database with stress shapes
        modal_db = ModalDatabase.from_op2_f06(
            op2_file=op2_path,
            f06_file=f06_path,
            damping=damping,
            element_ids=element_ids,  # Request stress shapes
        )
        
        # Load PSD
        psd = get_psd_profile(psd_profile) if not Path(psd_profile).exists() \
              else PSDProfile.from_csv(Path(psd_profile))
        
        # Process
        processor = RandomVibrationProcessor(modal_db, damping)
        result = processor.compute_stress_response(
            psd_profile=psd,
            direction=direction,
            element_ids=element_ids,
            peak_sigma=peak_sigma,
            non_gaussian_correction=non_gaussian_correction,
        )
        
        return {
            'max_rms_von_mises': result['max_rms_von_mises'],
            'max_peak_von_mises': result['max_peak_von_mises'],
            'critical_element': result['critical_element'],
            'n_elements': result['n_elements'],
            'success': True,
            'error': None,
        }
        
    except Exception as e:
        logger.exception(f"Stress extraction failed: {e}")
        return {
            'max_rms_von_mises': None,
            'max_peak_von_mises': None,
            'critical_element': None,
            'success': False,
            'error': str(e),
        }


# Convenience functions for optimization objectives

def get_grms(op2_file, psd='GEVS', direction='z', damping=0.02) -> float:
    """Quick Grms extraction for objective function."""
    result = extract_random_response(op2_file, psd, direction, damping)
    return result['grms'] if result['success'] else float('inf')


def get_peak_g(op2_file, psd='GEVS', direction='z', damping=0.02, sigma=3.0) -> float:
    """Quick peak acceleration extraction for constraint."""
    result = extract_random_response(op2_file, psd, direction, damping, peak_sigma=sigma)
    return result['peak_acceleration'] if result['success'] else float('inf')


def get_peak_stress(op2_file, psd='GEVS', direction='z', damping=0.02, sigma=3.0) -> float:
    """Quick peak stress extraction for constraint."""
    result = extract_random_stress(op2_file, psd, direction, damping, sigma)
    return result['max_peak_von_mises'] if result['success'] else float('inf')

Phase 4: Sine Sweep Processor (Week 4)

Goal: Add harmonic/sine sweep response capability.

Tasks

Task	Description	Effort	Deliverable
4.1	Sine sweep transfer function	3h	Magnitude/phase at frequencies
4.2	Resonance search	2h	Find peak frequencies
4.3	Dwell analysis	3h	Response at specific frequencies
4.4	Sweep rate effects	2h	Transient buildup correction
4.5	Create SineSweepProcessor	4h	sine_sweep.py
4.6	Create extractor	2h	extract_sine_response.py
4.7	Optimization template	2h	sine_sweep_optimization.json
4.8	Tests	3h	Validate against analytical

Phase 5: Advanced Features (Future)

Goal: Extended capabilities for complete dynamic qualification.

Shock Response Spectrum (SRS)

• Compute SRS from transient pulse or base motion
• Pseudo-velocity and pseudo-acceleration
• Compare against equipment SRS limits

Multi-Axis Excitation

• SRSS (Square Root Sum of Squares) combination
• CQC (Complete Quadratic Combination) for correlated inputs
• Full 6-DOF base motion

Fatigue from Random Vibration

• Narrow-band fatigue using Miner's rule
• Dirlik method for wide-band random
• Cycle counting from PSD

Neural Surrogate for Modal Database

• Train GNN to predict (frequencies, mode shapes) from geometry
• Skip SOL 103 entirely for similar designs
• 1000x additional speedup potential

---

4. Module Specifications

4.1 ModalDatabase

class ModalDatabase:
    """
    Central storage for all modal analysis data.
    
    Responsibilities:
    - Extract data from OP2/F06 files
    - Store mode shapes, frequencies, participation factors
    - Provide efficient access patterns for dynamic calculations
    - Support caching and serialization
    
    Thread Safety: Read-only after construction
    Memory: O(n_modes × n_nodes × 6) for full eigenvectors
    """
    
    # Class attributes
    SUPPORTED_ELEMENTS = ['CQUAD4', 'CTRIA3', 'CHEXA', 'CPENTA', 'CTETRA']
    
    # Required data
    model_name: str
    total_mass: float
    center_of_gravity: np.ndarray
    modes: List[ModeShape]
    
    # Optional data
    node_coordinates: Dict[int, np.ndarray]
    element_connectivity: Dict[int, List[int]]
    
    # Metadata
    source_op2: Path
    source_f06: Path
    extraction_timestamp: datetime
    nastran_version: str
    
    # Methods
    @classmethod
    def from_op2_f06(cls, op2_file, f06_file, **options) -> 'ModalDatabase'
    
    def get_mode(self, mode_number: int) -> ModeShape
    def get_frequencies(self) -> np.ndarray
    def get_participation_factors(self, direction: str) -> np.ndarray
    def get_cumulative_mass_fraction(self, direction: str) -> np.ndarray
    def modes_for_mass_fraction(self, fraction: float, direction: str) -> int
    
    def to_hdf5(self, filepath: Path)
    @classmethod
    def from_hdf5(cls, filepath: Path) -> 'ModalDatabase'
    
    def to_dict(self) -> Dict
    def summary(self) -> str

4.2 TransferFunctionEngine

class TransferFunctionEngine:
    """
    Computes frequency response functions from modal data.
    
    Responsibilities:
    - Single mode FRF calculation
    - Multi-mode combination
    - Base excitation transmissibility
    - Point excitation mobility (future)
    
    Performance: Vectorized numpy operations
    Accuracy: Double precision complex arithmetic
    """
    
    def __init__(self, modal_db: ModalDatabase)
    
    # Core FRF methods
    def single_mode_frf(
        self,
        frequencies: np.ndarray,
        mode: ModeShape,
        output_type: str = 'displacement'
    ) -> np.ndarray  # Complex FRF
    
    def base_excitation_frf(
        self,
        frequencies: np.ndarray,
        direction: str,
        output_type: str = 'acceleration',
        node_id: int = None,
        n_modes: int = None
    ) -> np.ndarray  # Complex transmissibility
    
    def point_excitation_frf(
        self,
        frequencies: np.ndarray,
        input_node: int,
        input_direction: str,
        output_node: int,
        output_direction: str
    ) -> np.ndarray  # Complex FRF
    
    # Stress FRF
    def stress_frf(
        self,
        frequencies: np.ndarray,
        element_id: int,
        direction: str
    ) -> np.ndarray  # Complex stress transfer function

4.3 RandomVibrationProcessor

class RandomVibrationProcessor:
    """
    Computes random vibration response using modal superposition.
    
    This is the SATK-equivalent processor.
    
    Capabilities:
    - Grms and peak response calculation
    - Response PSD computation
    - Modal contribution analysis
    - Stress recovery with non-Gaussian correction
    - Multiple response locations
    
    Performance Target: <100ms for typical aerospace model
    """
    
    def __init__(self, modal_db: ModalDatabase, default_damping: float = 0.02)
    
    # Main computation methods
    def compute_response(
        self,
        psd_profile: PSDProfile,
        direction: str = 'z',
        node_id: int = None,
        n_modes: int = None,
        frequency_resolution: int = 500,
        peak_sigma: float = 3.0,
        output_type: str = 'acceleration'
    ) -> RandomVibrationResult
    
    def compute_stress_response(
        self,
        psd_profile: PSDProfile,
        direction: str = 'z',
        element_ids: List[int] = None,
        peak_sigma: float = 3.0,
        non_gaussian_correction: bool = True
    ) -> Dict[str, Any]
    
    # Quick objective/constraint methods
    def compute_grms(self, psd_profile, direction, node_id=None) -> float
    def compute_peak_g(self, psd_profile, direction, sigma=3.0) -> float
    def compute_peak_stress(self, psd_profile, direction, sigma=3.0) -> float
    
    # Analysis methods
    def get_modal_contributions(self, psd_profile, direction) -> Dict[int, float]
    def get_response_psd(self, psd_profile, direction) -> Tuple[np.ndarray, np.ndarray]
    def identify_critical_modes(self, psd_profile, direction, threshold=0.1) -> List[int]

---

5. Integration Points

5.1 Extractor Registry

Add to optimization_engine/extractors/__init__.py:

# Dynamic response extractors
from .extract_modal_database import extract_modal_database
from .extract_random_response import (
    extract_random_response,
    extract_random_stress,
    get_grms,
    get_peak_g,
    get_peak_stress
)
from .extract_sine_response import extract_sine_response

__all__ = [
    # ... existing extractors ...
    
    # Dynamic response
    'extract_modal_database',
    'extract_random_response',
    'extract_random_stress',
    'get_grms',
    'get_peak_g', 
    'get_peak_stress',
    'extract_sine_response',
]

5.2 Feature Registry

Add to optimization_engine/feature_registry.json:

{
  "dynamic_response": {
    "random_vibration_processor": {
      "feature_id": "random_vibration_processor",
      "name": "Random Vibration Processor",
      "description": "Modal superposition-based random vibration response calculation",
      "category": "dynamics",
      "implementation": {
        "file_path": "optimization_engine/processors/dynamic_response/random_vibration.py",
        "class_name": "RandomVibrationProcessor"
      },
      "interface": {
        "inputs": [
          {"name": "modal_database", "type": "ModalDatabase"},
          {"name": "psd_profile", "type": "PSDProfile"},
          {"name": "direction", "type": "str"}
        ],
        "outputs": [
          {"name": "grms", "type": "float"},
          {"name": "peak_acceleration", "type": "float"},
          {"name": "response_psd", "type": "np.ndarray"}
        ]
      }
    }
  }
}

5.3 Dashboard Integration

New dashboard components for dynamic response:

// atomizer-dashboard/src/components/DynamicResponse/
├── ResponsePSDPlot.tsx      // Log-log PSD plot with input/output
├── TransmissibilityPlot.tsx // FRF magnitude vs frequency
├── ModalContributionBar.tsx // Bar chart of mode contributions
├── GrmsHistoryChart.tsx     // Grms vs trial number
├── CampbellDiagram.tsx      // For rotating machinery (future)
└── DynamicResponseCard.tsx  // Summary card for dashboard

5.4 CLI Integration

# New CLI commands
atomizer dynamic extract model.op2 --psd GEVS --direction z
atomizer dynamic psd-list
atomizer dynamic psd-plot GEVS --output psd.png
atomizer optimize --template random_vibration_optimization.json

---

6. Validation Strategy

6.1 Unit Tests

# tests/test_dynamic_response/

test_modal_database.py
├── test_mode_shape_creation
├── test_op2_extraction
├── test_f06_parsing
├── test_mass_fraction_calculation
└── test_hdf5_caching

test_transfer_functions.py
├── test_sdof_frf_analytical
├── test_resonance_amplification
├── test_phase_at_resonance
├── test_multi_mode_combination
└── test_base_excitation_transmissibility

test_random_vibration.py
├── test_miles_equation_sdof
├── test_grms_calculation
├── test_peak_factor
├── test_non_gaussian_correction
├── test_modal_contributions
└── test_stress_recovery

test_psd_profiles.py
├── test_interpolation
├── test_grms_integration
├── test_standard_profiles
└── test_csv_loading

6.2 Integration Tests

# tests/integration/

test_dynamic_optimization.py
├── test_full_random_vib_optimization
├── test_multi_objective_dynamic
├── test_constraint_handling
└── test_dashboard_updates

test_nastran_comparison.py
├── test_vs_sol111_response
├── test_vs_sol112_random
└── test_stress_recovery_accuracy

6.3 Validation Cases

Case	Model	Validation Against	Acceptance
SDOF	Single mass-spring	Miles' equation	<1% error
Cantilever beam	Simple beam	SOL 111 FRF	<5% error
Plate	Rectangular plate	SOL 112 random	<5% error
Bracket	Real component	SATK results	<10% error

6.4 Benchmark Against SATK

If possible, run same model through both Atomizer and SATK:

Model: CubeSat bracket
PSD: GEVS
Direction: Z-axis
Damping: 2%

Compare:
- Grms
- Peak acceleration
- Dominant mode
- Response PSD shape
- Peak stress location

---

7. Performance Targets

7.1 Timing Budgets

Operation	Target	Notes
OP2 read (50 modes)	<2s	pyNastran
F06 parse (MEFFMASS)	<0.5s	Regex parsing
Modal DB construction	<1s	Data structuring
FRF calculation (1000 freq)	<20ms	Vectorized
Response PSD	<20ms	Element-wise
Spectral integration	<10ms	Trapezoidal
Stress recovery (1000 elements)	<200ms	Parallel possible
Total per trial	<3s	Excluding SOL 103

7.2 Memory Budgets

Data	Size Estimate	Notes
Modal DB (100 modes, 10k nodes)	~50 MB	Full eigenvectors
Modal DB (100 modes, 100 nodes)	~5 MB	Sparse storage
Response PSD (1000 freq)	<100 KB	Just arrays
Stress shapes (100 modes, 1k elements)	~50 MB	If needed

7.3 Scalability

Tested configurations:

Model Size	Modes	Nodes	Elements	Expected Time
Small	50	5,000	10,000	<2s
Medium	100	50,000	100,000	<5s
Large	200	200,000	500,000	<15s
Very Large	500	1,000,000	2,000,000	<60s

---

8. Risk Mitigation

8.1 Technical Risks

Risk	Probability	Impact	Mitigation
pyNastran can't read NX OP2	Low	High	Test early, have F06 fallback
Modal stress not in OP2	Medium	Medium	Request STRESS(MODAL) output
Poor accuracy vs SOL 111	Low	High	Validate thoroughly, document assumptions
Memory issues large models	Medium	Medium	Lazy loading, sparse storage
Complex modes (damped)	Low	Low	Detect and warn, use magnitude

8.2 Mitigation Actions

1. Early OP2 Testing: Test pyNastran with NX 2412 OP2 files in week 1
2. F06 Fallback: Implement F06 parsing for critical data
3. Validation Suite: Build comprehensive test cases against SOL 111
4. Memory Profiling: Profile with large models, implement streaming if needed
5. Documentation: Clear documentation of assumptions and limitations

---

9. Future Extensions

9.1 Short-Term (3-6 months)

• Sine sweep processor
• Shock response spectrum
• Multi-axis excitation
• Response at multiple nodes
• Automatic mode selection

9.2 Medium-Term (6-12 months)

• Fatigue from random vibration
• Acoustic loading (diffuse field)
• Coupled structural-acoustic
• Time-domain transient from PSD
• Design sensitivity (∂response/∂parameter)

9.3 Long-Term (12+ months)

• Neural surrogate for modal database
• Real-time response visualization
• Uncertainty quantification
• Multi-physics coupling (thermal + dynamic)
• Foundation model for structural dynamics

---

10. Success Metrics

10.1 Technical Metrics

• Grms accuracy <5% vs SOL 111
• Peak stress accuracy <10% vs SOL 112
• Processing time <5s per trial
• Memory <500MB for large models
• 100% test coverage for core modules

10.2 User Metrics

• Complete optimization in <1 hour (100 trials)
• Dashboard shows real-time dynamic metrics
• Documentation enables self-service usage
• Template covers 80% of use cases

10.3 Business Metrics

• Feature demo ready for marketing video
• 3+ real-world optimization case studies
• Positive feedback from beta users
• Competitive differentiation documented

---

Appendix A: Reference Equations

A.1 Modal Superposition

Physical response:
u(x,t) = Σᵢ φᵢ(x) · qᵢ(t)

Modal equation of motion:
q̈ᵢ + 2ζᵢωᵢq̇ᵢ + ωᵢ²qᵢ = Lᵢ·a_base(t) / mᵢ

Frequency response function:
Hᵢ(ω) = 1 / (ωᵢ² - ω² + j·2ζᵢωᵢω)

A.2 Random Vibration

Response PSD:
Sᵧ(f) = |H(f)|² · Sₓ(f)

Mean square:
σ² = ∫₀^∞ Sᵧ(f) df

Peak (Gaussian):
peak = k · σ,  k ≈ 3 for 3-sigma

A.3 Miles' Equation

For SDOF with white noise input:
Grms = √(π/2 · fₙ · Q · PSD)

Where:
- fₙ = natural frequency (Hz)
- Q = 1/(2ζ) = quality factor
- PSD = input PSD level (G²/Hz)

---

Appendix B: File Formats

B.1 OP2 Data Blocks

Table	Contains	Used For
OUGV1	Eigenvectors	Mode shapes
LAMA	Eigenvalues	Frequencies
OES1	Element stress	Modal stress shapes
OGPFB1	Grid point forces	Participation
OQMG1	Modal participation	Effective mass

B.2 F06 Tables

Section	Contains	Parsing Pattern
REAL EIGENVALUES	Frequencies, gen mass	Regex table
MODAL EFFECTIVE MASS	Participation	Regex table
MODAL PARTICIPATION FACTORS	L factors	Regex table

B.3 PSD CSV Format

frequency_hz,psd_g2_hz
20,0.026
50,0.16
800,0.16
2000,0.026

---

Document End - Dynamic Response Processor Master Plan