diff --git a/docs/HYDROGEN_TANK_SHM.md b/docs/HYDROGEN_TANK_SHM.md index 8ebb254..b7437a4 100644 --- a/docs/HYDROGEN_TANK_SHM.md +++ b/docs/HYDROGEN_TANK_SHM.md @@ -119,11 +119,17 @@ GW グラフ（センサ=ノード）は `build_gw_graph.py` のスキーマを ## 7. 1サンプル FEM 実行計画 (One-Sample Plan) プロジェクト流儀に従い、**まず1サンプルで検証**してからバッチ化する。 +着手済みコード: `src/generate_cryotank_doe.py`（DOE, Abaqus 非依存・テスト済み）、 +`src/generate_cryotank_dataset.py`（Abaqus CAE 雛形, `--dry-run` は非 Abaqus で解析プラン確認可）。 ```bash -# 1) 生成スクリプト雛形（新規, generate_fairing_dataset.py ベース） -# → src/generate_cryotank_dataset.py（healthy 1件） -abaqus cae noGUI=src/generate_cryotank_dataset.py -- --config doe_cryotank_1sample.json +# 0) DOE 生成（healthy 1件, Abaqus 不要） +python src/generate_cryotank_doe.py --healthy_only --n_samples 1 --output doe_cryotank_1sample.json +# ローカルで解析プランを確認（Abaqus 不要） +python src/generate_cryotank_dataset.py --dry-run --defect doe_cryotank_1sample.json + +# 1) FEM 生成（クラスタ, Abaqus）— 雛形の M1/M2 TODO を実装後に実行 +abaqus cae noGUI=src/generate_cryotank_dataset.py -- --job CryoTank_Healthy_0000 --defect doe_cryotank_1sample.json # 2) ODB 抽出（既存を流用） abaqus python src/extract_odb_results.py --odb abaqus_work/Job-CryoTank-Healthy.odb @@ -157,8 +163,9 @@ cd src && python build_graph.py --data_dir ../dataset_cryotank_1sample ## 9. 次アクション +- [x] `src/generate_cryotank_doe.py` / `src/generate_cryotank_dataset.py` 雛形作成（DOE + dry-run はテスト済み） - [ ] 本メモのレビュー・数値の確定（材料・寸法・荷重） -- [ ] `src/generate_cryotank_dataset.py` 雛形作成（healthy 1件） +- [ ] Abaqus 雛形の M1/M2 TODO 実装（sector part・材料・内圧+熱・欠陥）→ クラスタで healthy 1件生成 - [ ] M1 検証 → 本メモに結果追記 - [ ] `ROADMAP.md` / `CLAUDE.md` の Research Lines に「Hydrogen Tank SHM」を正式追加 diff --git a/src/generate_cryotank_dataset.py b/src/generate_cryotank_dataset.py new file mode 100644 index 0000000..f753b26 --- /dev/null +++ b/src/generate_cryotank_dataset.py @@ -0,0 +1,154 @@ +# -*- coding: utf-8 -*- +# generate_cryotank_dataset.py +# Abaqus CAE scaffold — H3 cryogenic LH2 tank (Al-Li) FEM with SHM defects. +# +# Design note: docs/HYDROGEN_TANK_SHM.md +# This is the SCAFFOLD for milestone M1 (memo §7): generate ONE healthy sample, +# eyeball the graph, then batch. Abaqus model-building bodies are marked TODO and +# should reuse the fairing generator logic (generate_fairing_dataset.py: +# sector part, symmetry BCs, mesh control; generate_cohesive_fairing.py / +# generate_czm_sector12.py: cohesive/CZM interfaces for debond & weld flaw). +# +# The Abaqus import is guarded so this module can be imported (and --dry-run +# executed) WITHOUT Abaqus — the cluster runs the real build: +# abaqus cae noGUI=src/generate_cryotank_dataset.py -- --job --defect +# Local plan check (no Abaqus): +# python src/generate_cryotank_dataset.py --dry-run --defect doe_cryotank_1sample.json + +import argparse +import json +import sys + +try: + from abaqus import * # noqa: F401,F403 + from abaqusConstants import * # noqa: F401,F403 + from caeModules import * # noqa: F401,F403 + from driverUtils import executeOnCaeStartup + _HAVE_ABAQUS = True +except Exception: + _HAVE_ABAQUS = False + +# ============================================================================== +# PARAMETERS — Al-Li LH2 tank (design assumption; memo §1-§2, dimensions non-public) +# ============================================================================== +RADIUS = 2600.0 # mm, φ ≈ 5.2 m barrel +H_BARREL = 4000.0 # mm, barrel section modeled +SECTOR_DEG = 30.0 # deg, symmetric sector (1/12) +T_WALL = 4.0 # mm, wall thickness (design assumption) + +# Al-Li material — RT baseline + cryogenic shift (memo §2; DESIGN ASSUMPTION, verify). +MAT_ALLI = { + "name": "AlLi_2219", + "rt": {"E_MPa": 75000.0, "nu": 0.33, "yield_MPa": 400.0, "alpha_per_C": 23.0e-6, + "density_t_mm3": 2.70e-9}, + # cryogenic (20 K) multipliers relative to RT + "cryo_factor": {"E": 1.10, "yield": 1.20, "alpha": 0.55}, +} + +REFERENCE_TEMP_K = 293.0 + + +def material_at(temperature_K): + """Return effective Al-Li properties at a temperature (RT interpolation endpoints).""" + rt = MAT_ALLI["rt"] + if temperature_K <= 20.0: + cf = MAT_ALLI["cryo_factor"] + return {"E_MPa": rt["E_MPa"] * cf["E"], "nu": rt["nu"], + "yield_MPa": rt["yield_MPa"] * cf["yield"], + "alpha_per_C": rt["alpha_per_C"] * cf["alpha"], + "density_t_mm3": rt["density_t_mm3"], "temperature_K": temperature_K} + return {"E_MPa": rt["E_MPa"], "nu": rt["nu"], "yield_MPa": rt["yield_MPa"], + "alpha_per_C": rt["alpha_per_C"], "density_t_mm3": rt["density_t_mm3"], + "temperature_K": temperature_K} + + +def resolve_plan(defect_params): + """Build a human-readable analysis plan from defect_params (no Abaqus needed).""" + dtype = defect_params.get("defect_type", "healthy") + temp = float(defect_params.get("temperature_K", REFERENCE_TEMP_K)) + press = float(defect_params.get("pressure_MPa", 0.3)) + mat = material_at(temp) + plan = { + "geometry": {"radius_mm": RADIUS, "h_barrel_mm": H_BARREL, + "sector_deg": SECTOR_DEG, "wall_mm": T_WALL}, + "material": mat, + "loading": {"internal_pressure_MPa": press, + "thermal": {"reference_K": REFERENCE_TEMP_K, "operating_K": temp, + "delta_K": temp - REFERENCE_TEMP_K}}, + "defect": {"type": dtype}, + } + if dtype != "healthy": + plan["defect"].update({ + "theta_deg": defect_params.get("theta_deg"), + "z_center_mm": defect_params.get("z_center"), + "radius_mm": defect_params.get("radius"), + "model": { + "weld_flaw": "cohesive/element weakening on weld line", + "microcrack": "local stiffness loss + partial contact discontinuity", + "insulation_debond": "interface cohesive degradation (reuse CZM)", + }.get(dtype, "TODO"), + }) + return plan + + +# ============================================================================== +# Abaqus build (cluster only) — bodies TODO, reuse fairing generator logic +# ============================================================================== +def build_model(job_name, defect_params): + if not _HAVE_ABAQUS: + raise RuntimeError("Abaqus not available. Use --dry-run for a local plan check, " + "or run via: abaqus cae noGUI=generate_cryotank_dataset.py -- ...") + executeOnCaeStartup() + plan = resolve_plan(defect_params) + # TODO(M1): create sector shell part (reuse generate_fairing_dataset.py sector logic) + # TODO(M1): assign Al-Li material from plan["material"]; set expansion (alpha) + ref temp + # TODO(M1): sections + assembly + symmetry BCs on circumferential edges + # TODO(M1): Step-1 static: internal pressure (plan["loading"]) + Predefined Field + # temperature (operating_K) for cryogenic thermal stress + # TODO(M2): insert defect per plan["defect"]["model"] + # - insulation_debond / weld_flaw: cohesive interface (generate_cohesive_fairing.py, + # generate_czm_sector12.py) + # - microcrack: local stiffness reduction + seam + # TODO(M1): mesh (continuum/solid shell) + element type + # TODO(M1): create + write Job(job_name) + raise NotImplementedError("Abaqus model build is a scaffold — implement M1 bodies.") + + +def main(): + p = argparse.ArgumentParser(description="Cryogenic tank FEM generator (Abaqus)") + p.add_argument("--job", type=str, default="CryoTank_Healthy_0000") + p.add_argument("--defect", type=str, default=None, + help="JSON: single defect_params, or a DOE {'samples':[...]} (uses sample 0)") + p.add_argument("--dry-run", action="store_true", default=False, + help="Print the resolved analysis plan without Abaqus") + # Abaqus passes script args after '--'; argparse needs them isolated. + argv = sys.argv[1:] + if "--" in argv: + argv = argv[argv.index("--") + 1:] + args = p.parse_args(argv) + + defect_params = {"defect_type": "healthy", "temperature_K": 20.0, "pressure_MPa": 0.3} + if args.defect: + with open(args.defect) as f: + data = json.load(f) + if isinstance(data, dict) and "samples" in data: + defect_params = data["samples"][0]["defect_params"] + elif isinstance(data, dict) and "defect_params" in data: + defect_params = data["defect_params"] + else: + defect_params = data + + if args.dry_run or not _HAVE_ABAQUS: + plan = resolve_plan(defect_params) + print(json.dumps({"job": args.job, "abaqus_available": _HAVE_ABAQUS, "plan": plan}, + indent=2)) + if not args.dry_run and not _HAVE_ABAQUS: + print("\n[note] Abaqus not found — printed plan only. " + "Run on the cluster to build the model.", file=sys.stderr) + return + + build_model(args.job, defect_params) + + +if __name__ == "__main__": + main() diff --git a/src/generate_cryotank_doe.py b/src/generate_cryotank_doe.py new file mode 100644 index 0000000..a151a65 --- /dev/null +++ b/src/generate_cryotank_doe.py @@ -0,0 +1,155 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Cryogenic Hydrogen-Tank DOE Generation — LHS over defect + operating conditions. + +Design note: docs/HYDROGEN_TANK_SHM.md (geometry, materials, loading, defects). +This is the *first* code step of the hydrogen-tank SHM line (memo §7, milestone M1). + +Pure-Python (no Abaqus) — runnable/verifiable outside the cluster. Produces a JSON +DOE consumed by src/generate_cryotank_dataset.py (Abaqus), mirroring the fairing +DOE schema (samples[].defect_params + metadata) so downstream tooling is reusable. + +Defect types (memo §4): + - healthy — no defect (baseline) + - weld_flaw — porosity / lack-of-fusion on a weld line (located ON a weld) + - microcrack — thermal-cycle micro-cracking (local stiffness loss) + - insulation_debond — foam/MLI separation (same interface type as skin-core debond) + +Operating point (memo §2-§3): + - temperature_K : 20 K (LH2 cryo) or 293 K (RT), stratified + - pressure_MPa : operational internal pressure + +Usage: + python src/generate_cryotank_doe.py --n_samples 1 --output doe_cryotank_1sample.json + python src/generate_cryotank_doe.py --n_samples 100 --output doe_cryotank_100.json + python src/generate_cryotank_doe.py --healthy_only --output doe_cryotank_healthy.json +""" + +import argparse +import json + +import numpy as np + +try: + from scipy.stats.qmc import LatinHypercube + _HAVE_QMC = True +except Exception: # scipy missing / too old — fall back to uniform RNG + _HAVE_QMC = False + +# ---- Geometry bounds (design assumption; see memo §1) ------------------------ +# Al-Li cylindrical LH2 tank, barrel section (representative, dimensions non-public). +THETA_RANGE = (5.0, 85.0) # deg, margin from symmetry edges (1/12--1/6 sector) +Z_RANGE = (400.0, 3600.0) # mm, along barrel, avoid dome junctions +WELD_THETA_DEG = [15.0, 45.0, 75.0] # longitudinal weld-line angular positions + +# ---- Defect size tiers (mm): (name, r_min, r_max, fraction) ------------------ +SIZE_TIERS = [ + ("Small", 20.0, 60.0, 0.35), + ("Medium", 60.0, 120.0, 0.40), + ("Large", 120.0, 220.0, 0.25), +] + +# ---- Operating point (memo §2-§3; design assumptions) ------------------------ +TEMPERATURE_K = [20.0, 293.0] # LH2 cryo / room-temperature +PRESSURE_MPA_RANGE = (0.2, 0.5) # operational internal pressure + +DEFECT_TYPES = ["weld_flaw", "microcrack", "insulation_debond"] + + +def _lhs(n, d, seed): + """Latin Hypercube samples in [0,1)^ (n x d); uniform fallback if no scipy.qmc.""" + if _HAVE_QMC: + return LatinHypercube(d=d, seed=seed).random(n) + rng = np.random.default_rng(seed) + return rng.random((n, d)) + + +def _pick_tier(u): + """Map u in [0,1) to a size tier by cumulative fraction, return (name, rmin, rmax).""" + acc = 0.0 + for name, rmin, rmax, frac in SIZE_TIERS: + acc += frac + if u < acc: + return name, rmin, rmax + name, rmin, rmax, _ = SIZE_TIERS[-1] + return name, rmin, rmax + + +def generate(n_samples, seed=42, healthy_only=False, defect_types=None): + """Return a DOE dict: {'samples': [...], 'metadata': {...}}.""" + defect_types = list(defect_types) if defect_types else list(DEFECT_TYPES) + samples = [] + + if healthy_only: + for i in range(n_samples): + t_u = _lhs(n_samples, 2, seed)[i] + samples.append({ + "id": i, + "job_name": "CryoTank_Healthy_%04d" % i, + "defect_params": { + "defect_type": "healthy", + "temperature_K": float(TEMPERATURE_K[int(t_u[0] * len(TEMPERATURE_K)) % len(TEMPERATURE_K)]), + "pressure_MPa": round(float(PRESSURE_MPA_RANGE[0] + t_u[1] * + (PRESSURE_MPA_RANGE[1] - PRESSURE_MPA_RANGE[0])), 4), + }, + }) + return {"samples": samples, "metadata": {"n_samples": n_samples, + "healthy_only": True, "seed": seed}} + + # 5 design dims: [defect_type, theta, z, size, operating(temp,press packed via 2 dims)] + u = _lhs(n_samples, 6, seed) + for i in range(n_samples): + dtype = defect_types[int(u[i, 0] * len(defect_types)) % len(defect_types)] + + # weld_flaw sits ON a weld line; others use continuous theta + if dtype == "weld_flaw": + theta = float(WELD_THETA_DEG[int(u[i, 1] * len(WELD_THETA_DEG)) % len(WELD_THETA_DEG)]) + else: + theta = round(float(THETA_RANGE[0] + u[i, 1] * (THETA_RANGE[1] - THETA_RANGE[0])), 2) + + z = round(float(Z_RANGE[0] + u[i, 2] * (Z_RANGE[1] - Z_RANGE[0])), 1) + tier, rmin, rmax = _pick_tier(u[i, 3]) + radius = round(float(rmin + u[i, 4] * (rmax - rmin)), 1) + temp = float(TEMPERATURE_K[int(u[i, 5] * len(TEMPERATURE_K)) % len(TEMPERATURE_K)]) + press = round(float(PRESSURE_MPA_RANGE[0] + (u[i, 5] * 7.0 % 1.0) * + (PRESSURE_MPA_RANGE[1] - PRESSURE_MPA_RANGE[0])), 4) + + samples.append({ + "id": i, + "job_name": "CryoTank_%s_%04d" % (dtype, i), + "defect_params": { + "defect_type": dtype, + "theta_deg": theta, + "z_center": z, + "radius": radius, + "size_tier": tier, + "temperature_K": temp, + "pressure_MPa": press, + }, + }) + + return {"samples": samples, + "metadata": {"n_samples": n_samples, "seed": seed, + "defect_types": defect_types, + "size_tiers": [t[0] for t in SIZE_TIERS]}} + + +def main(): + p = argparse.ArgumentParser(description="Cryogenic hydrogen-tank DOE generator") + p.add_argument("--n_samples", type=int, default=1) + p.add_argument("--output", type=str, default="doe_cryotank_1sample.json") + p.add_argument("--seed", type=int, default=42) + p.add_argument("--healthy_only", action="store_true", default=False) + p.add_argument("--defect_types", nargs="+", default=None, choices=DEFECT_TYPES) + args = p.parse_args() + + doe = generate(args.n_samples, seed=args.seed, + healthy_only=args.healthy_only, defect_types=args.defect_types) + with open(args.output, "w") as f: + json.dump(doe, f, indent=2) + print("Wrote %d samples -> %s" % (len(doe["samples"]), args.output)) + + +if __name__ == "__main__": + main() diff --git a/tests/test_cryotank_doe.py b/tests/test_cryotank_doe.py new file mode 100644 index 0000000..3698f45 --- /dev/null +++ b/tests/test_cryotank_doe.py @@ -0,0 +1,42 @@ +"""Smoke tests for the hydrogen-tank DOE generator + FEM scaffold (no Abaqus).""" + +import importlib + + +def test_cryotank_doe_generate_one(): + m = importlib.import_module("generate_cryotank_doe") + doe = m.generate(1, seed=42) + assert doe["metadata"]["n_samples"] == 1 + assert len(doe["samples"]) == 1 + dp = doe["samples"][0]["defect_params"] + assert dp["defect_type"] in (["healthy"] + m.DEFECT_TYPES) + assert dp["temperature_K"] in m.TEMPERATURE_K + + +def test_cryotank_doe_batch_and_bounds(): + m = importlib.import_module("generate_cryotank_doe") + doe = m.generate(50, seed=7) + assert len(doe["samples"]) == 50 + for s in doe["samples"]: + dp = s["defect_params"] + assert m.THETA_RANGE[0] <= dp["theta_deg"] <= m.THETA_RANGE[1] or dp["defect_type"] == "weld_flaw" + assert m.Z_RANGE[0] <= dp["z_center"] <= m.Z_RANGE[1] + assert dp["radius"] > 0 + + +def test_cryotank_healthy_only(): + m = importlib.import_module("generate_cryotank_doe") + doe = m.generate(3, healthy_only=True) + assert all(s["defect_params"]["defect_type"] == "healthy" for s in doe["samples"]) + + +def test_cryotank_dataset_dry_run_plan(): + m = importlib.import_module("generate_cryotank_dataset") + # cryogenic operating point → E should be scaled up vs RT + plan = m.resolve_plan({"defect_type": "healthy", "temperature_K": 20.0, "pressure_MPa": 0.3}) + assert plan["material"]["E_MPa"] > m.MAT_ALLI["rt"]["E_MPa"] + assert plan["loading"]["thermal"]["delta_K"] < 0 # cooling from RT reference + rt = m.resolve_plan({"defect_type": "weld_flaw", "temperature_K": 293.0, + "pressure_MPa": 0.3, "theta_deg": 45.0, "z_center": 2000, "radius": 80}) + assert rt["material"]["E_MPa"] == m.MAT_ALLI["rt"]["E_MPa"] + assert rt["defect"]["type"] == "weld_flaw"