rust-1d-premix-laminar-flame

実験的プロジェクト。 Rust で 1D 予混合層流火炎ソルバーを実装する個人実験です。

燃料・酸化剤混合気の未燃条件（温度・圧力・当量比）を与えると、火炎全域の温度・ 化学種質量分率の空間プロファイルおよび層流燃焼速度（固有値）を計算するソルバー。

アプローチ

PREMIX / Cantera の手法に準拠：

• 支配方程式: 1D 非一様格子上の化学種保存式 + エネルギー方程式
• 化学反応: Arrhenius、Troe フォールオフ、PLOG 圧力依存反応を含む詳細反応機構
• 輸送: 混合平均輸送（Chapman-Enskog 粘性・拡散、Mason-Monchick 熱伝導率、Wilke 混合則）
• 衝突積分: Monchick-Mason (1961) 2D テーブル（T* × δ*、極性分子対応）
• 熱力学: NASA 7 係数多項式
• 離散化: 有限差分（対流：風上差分、拡散：中心差分）
• ソルバー: 減衰ニュートン法（帯行列 LU + 枠付きシュア補完）+ 擬似過渡継続法 + 適応格子精細化（GRAD/CURV）
• 入力: Cantera YAML 反応機構ファイル、TOML 設定ファイル
• 出力: z, T, M, u, ρ, HRR, Xk, Yk の CSV プロファイル

参考実装: Kee et al., "PREMIX" (Sandia, 1985) および Cantera FreeFlame。

検証結果

H₂/air φ=1、300 K、1 atm（data/h2o2.yaml、Cantera 3.1.0 参照）:

残差 0.1% は風上差分の O(Δz) 打ち切り誤差によるもの

H₂/air 当量比スイープ（data/h2o2.yaml、300 K、1 atm）

ビルド・実行

cargo build
cargo run -- examples/hydrogen_air.toml

設定ファイル例（TOML）:

[mechanism]
file = "data/h2o2.yaml"
format = "cantera_yaml"

[flame]
pressure = 101325.0
fuel = { H2 = 1.0 }
oxidizer = { O2 = 0.21, N2 = 0.79 }
equivalence_ratio = 1.0
t_unburned = 300.0
domain_length = 0.02

[grid]
initial_points = 100
max_points = 1000
grad = 0.05
curv = 0.10

[solver]
atol = 1.0e-9
rtol = 1.0e-6
max_newton_iter = 50

可視化

scripts/plot_flame.html — ブラウザで動作するインタラクティブな火炎プロファイルビューア（Plotly.js 使用）。

• シングルモード: 4 パネル（T+HRR、主要種、ラジカル種対数、流速+密度）
• 比較モード: 最大 3 ファイルを重ねて比較（6 パネル + Su/T_max サマリーテーブル）
• 全トレースで lines+markers 表示（格子点位置を可視化）

テスト

cargo test           # 全テスト（ユニット + 統合）
cargo test --lib     # ユニットテストのみ

現在 83 テスト、全 PASS（統合テスト 1 件含む）:

参考文献

• Kee, Grcar, Smooke, Miller (1985). PREMIX. Sandia Report SAND85-8240.
• Kee, Rupley, Miller (1989). Chemkin-II. Sandia Report SAND89-8009.
• Monchick, Mason (1961). Transport Properties of Polar Gases. J. Chem. Phys. 35, 1676.
• Mason, Monchick (1962). Heat Conductivity of Polyatomic and Polar Gases. J. Chem. Phys. 36, 1622.
• Neufeld, Janzen, Aziz (1972). Empirical Equations to Calculate 16 of the Transport Collision Integrals. J. Chem. Phys. 57, 1100.
• Wilke (1950). A Viscosity Equation for Gas Mixtures. J. Chem. Phys. 18, 517.
• Cantera: https://cantera.org

Getting Started

必要なもの

• Rust 1.75 以上（rustup でインストール推奨）
• Git

確認コマンド:

rustc --version   # rustc 1.75.0 以上
cargo --version

インストール

git clone https://github.com/oktomoya/rust-1d-premix-laminar-flame.git
cd rust-1d-premix-laminar-flame
cargo build --release

初回ビルドは依存クレートのコンパイルで数分かかります。

最初の計算を走らせる

同梱の H₂/air φ=1.0 設定ファイルで計算してみましょう:

cargo run --release -- examples/hydrogen_air.toml

正常に収束すると、こんなサマリーが表示されます:

---------------------------------------------------
  Laminar flame speed Su = 2.3327 m/s
  Max temperature       = 2356.9 K
  Grid points           = 221
  Mass flux M           = 2.8294e-1 kg/(m²·s)
---------------------------------------------------
Solution written to /tmp/h2air_phi_1.0.csv

出力ファイルを可視化する

scripts/plot_flame.html をブラウザで開き、生成された CSV ファイルをドラッグ＆ドロップします。

または Python スクリプト:

pip install pandas matplotlib numpy
python scripts/plot_flame.py /tmp/h2air_phi_1.0.csv

設定ファイルの書き方

examples/hydrogen_air.toml を複製して編集するのが最も簡単です。

[mechanism]
file = "data/h2o2.yaml"       # Cantera YAML 反応機構ファイル
format = "cantera_yaml"

[flame]
pressure = 101325.0           # 圧力 [Pa]
fuel     = { H2 = 1.0 }       # 燃料組成（モル分率）
oxidizer = { O2 = 0.21, N2 = 0.79 }
equivalence_ratio = 1.0       # 当量比 φ
t_unburned = 300.0            # 未燃温度 [K]
domain_length = 0.02          # 計算領域長さ [m]

[grid]
initial_points = 100          # 初期格子点数
max_points = 1000             # 最大格子点数
grad = 0.05                   # GRAD 精細化基準（小さいほど細かい）
curv = 0.10                   # CURV 精細化基準

[solver]
atol = 1.0e-9
rtol = 1.0e-6
max_newton_iter = 50
su_initial_guess = 2.3        # 火炎速度初期推定 [m/s]（省略可）

[output]
file = "/tmp/h2air_phi_1.0.csv"

当量比を変えるには

equivalence_ratio を変えるだけです。収束が遅い場合は、隣の φ で計算した CSV を initial_profile に指定すると擬似過渡継続フェーズをスキップして高速に収束します:

[solver]
su_initial_guess = 1.6
initial_profile = "/tmp/h2air_phi_0.9.csv"   # 前の計算結果を初期値に使用

テストを実行する

cargo test                    # 全 83 テスト（ユニット + 統合）
cargo test --lib              # ユニットテストのみ（統合テストを除く）
cargo test -- --nocapture     # 標準出力を表示しながら実行

統合テスト (tests/h2air_validation.rs) はフルソルバーを走らせるため、 数秒かかります。

よくある問題

cargo run が遅い

--release フラグを忘れていませんか？デバッグビルドは最適化なしなので 化学反応計算が 10 倍以上遅くなります。

Newton が収束しない

• su_initial_guess を実際の火炎速度に近い値にしてみてください
• initial_profile に収束済み CSV を指定して PT フェーズをスキップしてみてください
• grad / curv を少し大きく（0.1 / 0.2）して格子を粗くしてみてください

出力 CSV が空 / 途中で止まる

ターミナルの出力に Newton: stalled が多発している場合は PT フェーズで詰まっています。 time_steps を増やすか（例: 1000）、dt_initial を小さく（例: 1e-8）してみてください。

chemistry モジュール

1D 火炎ソルバーの化学反応層。化学種データの保持・熱力学計算・反応速度計算・ 機構ファイルのパースを担当する。

モジュール構成

chemistry/
├── mod.rs            — サブモジュールの公開宣言
├── species.rs        — 化学種データ構造 (Species, NasaPoly, TransportParams)
├── mechanism.rs      — 機構全体のデータ構造 (Mechanism, Reaction, RateType)
├── thermo.rs         — NASA7 多項式による熱力学関数
├── kinetics.rs       — 反応速度計算・生成速度
└── parser/
    ├── mod.rs
    ├── cantera_yaml.rs — Cantera YAML パーサー ✅ 実装済み
    └── chemkin.rs      — CHEMKIN-II パーサー   ✅ 実装済み

species.rs — 化学種データ構造

Species

1 化学種が持つ全データをまとめた構造体。

nasa_poly(t) メソッドが温度 t に応じて nasa_low / nasa_high を自動選択する (切り替え温度 = nasa_low.t_high = T_mid)。

NasaPoly

NASA 7係数多項式 1レンジ分。有効温度範囲 [t_low, t_high] と係数配列 coeffs[7]。

cp/R = a[0] + a[1]·T + a[2]·T² + a[3]·T³ + a[4]·T⁴
h/RT = a[0] + a[1]/2·T + a[2]/3·T² + a[3]/4·T³ + a[4]/5·T⁴ + a[5]/T
s/R  = a[0]·ln(T) + a[1]·T + a[2]/2·T² + a[3]/3·T³ + a[4]/4·T⁴ + a[6]

TransportParams

Lennard-Jones パラメータと分子形状 (輸送係数モジュールで使用)。

mechanism.rs — 機構データ構造

RateType

反応速度定数の種別を enum で表現する。

#![allow(unused)]
fn main() {
pub enum RateType {
    Arrhenius { a: f64, b: f64, ea: f64 },
    Falloff { high: Arrhenius, low: Arrhenius, troe: Option<TroeParams> },
    Plog { rates: Vec<(f64, Arrhenius)> },
}
}

• Arrhenius: 通常の Arrhenius 式。三体反応も含む (third_body フィールドが Some のとき)。
  • ea の単位は J/mol (パーサーが変換済み)
  • a の単位はパース元ファイルのまま (Cantera YAML は CGS: cm³/mol/s 等)
• Falloff: 高圧限界 / 低圧限界の Arrhenius に加え、オプションの Troe 補正係数。 Lindemann 形式は troe = None。
• Plog: 圧力依存反応。(圧力 [Pa], Arrhenius) のリストで定義される。

TroeParams

Troe 補正のパラメータ。Cantera YAML の Troe: ブロックに対応。

Fcent = (1-A)·exp(-T/T3) + A·exp(-T/T1) + exp(-T2/T)

フィールド: a, t3, t1, t2: Option<f64> (T2 は省略可)。

ThirdBodySpec

三体衝突パートナーの拡張効率を保持する。

#![allow(unused)]
fn main() {
pub struct ThirdBodySpec {
    pub efficiencies: Vec<(usize, f64)>,  // (species_index, eff)
}
}

デフォルト効率は 1.0 (暗黙)。非デフォルト種のみ格納する。 [M]_eff = Σ[Ci] + Σ(eff_k - 1)·[Ck] で有効第三体濃度を計算する。

Reaction

thermo.rs — NASA7 熱力学関数

R_UNIVERSAL = 8.314462618 J/(mol·K) を定数として定義。

注意: entropy_species は J/(kg·K) を返す。kinetics.rs::equilibrium_constant では sp.molecular_weight を乗じて J/(mol·K) に戻してから使用している。

kinetics.rs — 反応速度・生成速度

arrhenius(a, b, ea, t)

k(T) = A · T^b · exp(−Ea / (R·T))

ea は J/mol。

production_rates(mech, t, concentrations, pressure)

全化学種のモル生成速度 ωk [mol/(cm³·s)] を返す。concentrations[k] = ρ·Yk/Wk·1e-6 [mol/cm³]。 A 値が CGS 単位 (cm³/mol/s) で格納されているため、濃度も mol/cm³ で渡す必要がある。

各反応について以下の順序で処理する:

1. 速度定数 kf を RateType に応じて計算:

   • Arrhenius → arrhenius(a, b, ea, t)
   • Falloff → falloff_rate(...) (後述)
   • Plog → plog_rate(...) (後述)

2. 三体補正:

   • Falloff の場合: falloff_rate 内で既に [M] を取り込み済みのため 追加乗算なし
   • それ以外で third_body.is_some() の場合: kf × [M]_eff

3. 平衡定数 Kc → 逆反応速度定数 kr = kf / Kc

4. 進行速度 q = kf·Π[Ri]^ν_i − kr·Π[Pi]^ν_i

5. 生成速度 ωk += Σ ν_prod · q − Σ ν_react · q

falloff_rate — Troe/Lindemann フォールオフ

k∞ = Arrhenius(high, T)
k0 = Arrhenius(low,  T)
Pr = k0·[M] / k∞
F  = troe_broadening(Troe, T, Pr)   ← Troe あり
   = 1.0                             ← Lindemann (Troe なし)
kf = k∞ · Pr/(1+Pr) · F

Troe ブロードニング係数 F の計算:

Fcent = (1−A)·exp(−T/T3) + A·exp(−T/T1) + exp(−T2/T)
c = −0.4 − 0.67·log10(Fcent)
n =  0.75 − 1.27·log10(Fcent)
f1 = (log10(Pr) + c) / (n − 0.14·(log10(Pr) + c))
F  = 10^(log10(Fcent) / (1 + f1²))

Cantera の TroeRate::updateTemp と同一の式。

plog_rate — PLOG 圧力依存反応

圧力方向に対数線形補間:

log k = log k1 + (log P − log P1) · (log k2 − log k1) / (log P2 − log P1)

圧力範囲外はクランプ (最低 / 最高圧エントリを使用)。

equilibrium_constant(mech, rxn_idx, t)

ΔG°/RT から Kp を求め、モル数変化 Δν を補正して Kc に変換する。

Kp = exp(−ΔG°/RT) = exp(−Σ νk·(hk − T·sk) / RT)
Kc = Kp · (P_atm / RT)^Δν

parser/cantera_yaml.rs — Cantera YAML パーサー ✅

parse_file(path) または parse_str(content) で Mechanism を返す。

パース手順

1. units: ブロックから活性化エネルギー単位を読む (デフォルト: cal/mol)
2. species: セクション → 各種 MW・NASA7・輸送パラメータ を構築
3. reactions: セクション → 反応タイプ判定 → RateType を構築

活性化エネルギー単位変換 (→ J/mol)

反応タイプ判定ロジック

type: pressure-dependent-Arrhenius  →  Plog
rate-constant: が null              →  Falloff  (high-P / low-P / Troe を読む)
それ以外                             →  Arrhenius

三体反応の判定: type: three-body または式中に (+M) / + M を含む場合、 efficiencies: マップから ThirdBodySpec を構築する。

化学量論係数パース

• 2 OH (スペース区切り係数)
• 2OH (係数が名前に直接付く)
• 係数なし → 1.0
• (+M) / + M / 孤立した M トークンはスキップ

分子量計算

composition: マップのキーを大文字正規化してテーブル参照。 対応元素: H, O, N, C, AR, HE, S, CL, F, BR。

テスト (9 テスト、全 PASS)

Cantera 3.1.0 Python API との数値クロスバリデーション (scripts/cross_validate.py): 164 チェック全 PASS (MW・NASA7係数・全29反応の A/b/Ea・Troe 4パラメータ・三体効率)。

parser/chemkin.rs — CHEMKIN-II パーサー ✅

対応ファイル形式

エントリポイント

#![allow(unused)]
fn main() {
// 文字列から直接パース
parse_mechanism(mech_src: &str, therm_src: Option<&str>, tran_src: Option<&str>) -> Result<Mechanism>

// ファイルパスから読み込み
parse_file(mech_path: &str, therm_path: Option<&str>, tran_path: Option<&str>) -> Result<Mechanism>
}

• therm_src / therm_path が None の場合は mech.inp 内の THERMO ブロックを参照
• tran_src / tran_path が None の場合は輸送データなし（デフォルト値）

パース手順

1. ELEMENTS ブロック: ELEMENTS … END を抽出、大文字正規化
2. SPECIES ブロック: SPECIES … END を抽出、宣言順を保持
3. therm.dat: 1/2/3/4 行インジケータ（列 79）で 4 行セットを認識し、固定列から NASA7 係数・温度範囲・元素組成・MW を読む
4. tran.dat: 名前 geometry ε/k σ μ α Z_rot の空白区切り形式
5. REACTIONS ブロック: 反応式行 + 修飾子行（LOW/, TROE/, REV/, 効率行, DUPLICATE）をレコードに集約して変換

反応修飾子

• (+M) または + M を含む方程式は自動的に三体反応と判定
• LOW/ があれば RateType::Falloff、なければ RateType::Arrhenius
• Ea はデフォルト cal/mol として J/mol に変換 (× 4.184)
• PLOG 形式は CHEMKIN-II 非標準のため非対応（Cantera YAML パーサーを使用）

テスト (13 テスト、全 PASS)

テストデータ: data/chemkin_test/{mech.inp, therm.dat, tran.dat}

transport モジュール

1D 火炎ソルバーの輸送係数層。衝突積分・化学種個別輸送係数・混合則を担当する。

モジュール構成

transport/
├── mod.rs                  — サブモジュールの公開宣言
├── collision_integrals.rs  — Neufeld (1972) 多項式 + Monchick-Mason (1961) 2D テーブル
├── species_props.rs        — 化学種単体の μk, λk, Dij
└── mixture.rs              — 混合則: μ_mix, λ_mix, Dkm

collision_integrals.rs — 衝突積分

Chapman-Enskog 輸送理論で使用する換算衝突積分を 2 つの方法で提供する。

Neufeld (1972) 多項式フィット — 非極性化学種向け

引数 t_star = T / (ε/kB) は換算温度（無次元）。非極性化学種（δ* = 0）のみ有効。

omega22(t_star) — Ω*(2,2)

粘性・熱伝導率の計算に使用。

Ω*(2,2) = A/T*^B + C/exp(D·T*) + E/exp(F·T*)

定数: A=1.16145, B=0.14874, C=0.52487, D=0.7732, E=2.16178, F=2.43787

omega11(t_star) — Ω*(1,1)

二成分拡散係数の計算に使用。

Ω*(1,1) = A/T*^B + C/exp(D·T*) + E/exp(F·T*) + G/exp(H·T*)

定数: A=1.06036, B=0.15610, C=0.19300, D=0.47635, E=1.03587, F=1.52996, G=1.76474, H=3.89411

Monchick-Mason (1961) 2D テーブル — 極性化学種対応

極性分子を含む化学種対に対し、換算双極子モーメント δ* の補正を施した衝突積分を テーブル補間で計算する。Cantera の MMCollisionInt と同一データ・同一補間法。

テーブルは T* = 0.1–100（37 点）× δ* = 0–2.5（8 点）の 2D グリッド。

omega22_mm(t_star, delta_star) — Ω*(2,2)(T*, δ*)

δ* = 0 の場合はテーブルの第 0 列から補間（Neufeld 多項式と < 0.3% の差）。

omega11_mm(t_star, delta_star) — Ω*(1,1)(T*, δ*)

A*(T*, δ*) = Ω*(2,2) / Ω*(1,1) テーブルを使用して omega22_mm / astar で計算。

delta_star_reduced(dipole_i, dipole_j, well_depth_ij, diameter_ij) — δ*

δ*_ij = 0.5 × μi × μj / (4πε₀ × εij × σij³)

引数の単位: 双極子モーメント [Debye]、複合井戸深さ [K]、複合径 [Å]。

変換定数 3622.85 は単位変換を吸収する（1 D², 1 K, 1 ų での δ* の値）。

h2o2.yaml では H₂O のみ双極子モーメントを持つ（μ = 1.844 D → δ*_H₂O ≈ 1.22）。

補間法

Cantera の MMCollisionInt::quadInterp と同一:

• δ* 方向: 隣接グリッド間の線形補間
• T* 方向: 3 点二次補間（log T* 座標）

テスト (11 テスト、全 PASS)

species_props.rs — 化学種輸送係数

viscosity(species, t) → [Pa·s]

Chapman-Enskog 理論による純粋化学種の粘性。

μk = 2.6693e-6 × √(Wk[g/mol] × T) / (σk[Å]² × Ω*(2,2)(T*, δ*k))

• δ*k > 0 の場合は omega22_mm を使用（H₂O のみ）
• δ*k = 0 の場合は Neufeld 多項式 omega22 を使用

thermal_conductivity(species, mu_k, cp_k, t, pressure) → [W/(m·K)]

完全な Mason-Monchick (1962) 式による熱伝導率。回転緩和補正（Zrot）を含む。

cv_rot:  0.0 (原子) | 1.0 (線形) | 1.5 (非線形)  [単位: R]
cv_int = cp/R - 2.5 - cv_rot   [振動モード、単位: R]

回転緩和係数 fz(T*):

fz(T*) = 1 + π^1.5/√T* · (0.5 + 1/T*)  + (π²/4 + 2)/T*

補正係数:

A  = 2.5 - f_int
B  = Zrot × fz(298/ε) / fz(T/ε)  +  (2/π) × (5/3 × cv_rot + f_int)
c1 = (2/π) × A / B
f_rot   = f_int × (1 + c1)
f_trans = 2.5 × (1 - c1 × cv_rot / 1.5)
λk = (μk/Wk) × R × (f_trans×1.5 + f_rot×cv_rot + f_int×cv_int)

binary_diffusion(sp_i, sp_j, t, pressure) → [m²/s]

BSL 式による二成分拡散係数。

Dij = 2.6280e-3 × T^1.5 / (P[atm] × σij[Å]² × Ω*(1,1)(T*, δ*ij) × √Wij[g/mol])

• σij = (σi + σj) / 2（算術平均）
• εij = √(εi × εj)（幾何平均）
• Wij = 2WiWj/(Wi+Wj)（換算分子量）
• δ*ij = delta_star_reduced(μi, μj, εij, σij)

μi = 0 または μj = 0 の場合（H₂O と非極性種の対）は δ*ij = 0 となり、Neufeld 多項式を使用。

テスト (7 テスト、全 PASS)

mixture.rs — 混合則

mixture_viscosity(mech, mole_fractions, t) → [Pa·s]

Wilke 混合則。

μ_mix = Σk Xk × μk / φk

相互作用係数:

φkj = [1 + (μk/μj)^0.5 × (Wj/Wk)^0.25]² / [√8 × √(1 + Wk/Wj)]
φk  = Σj Xj × φkj

mixture_thermal_conductivity(mech, mole_fractions, mass_fractions, t, pressure) → [W/(m·K)]

算術平均・調和平均の相乗平均（Wassilijewa 混合則）。

λ_mix = 0.5 × (Σk Xk × λk  +  1 / Σk Xk/λk)

mixture_diffusion_coefficients(mech, mole_fractions, mass_fractions, t, pressure) → Vec<f64> [m²/s]

混合平均拡散係数（各化学種 k に対して）。Cantera の getMixDiffCoeffs と一致。

Dkm = (1 - Yk) / Σ_{j≠k} (Xj / Dkj)

分子は Xk（モル分率）ではなく Yk（質量分率）を使用する。Cantera の定式化は (W̄ - Xk·Wk) / (W̄ · Σ Xj/Dkj) = (1 - Yk) / Σ Xj/Dkj。

mass_to_mole_fractions(mech, y) → Vec<f64>

質量分率 → モル分率の変換。W̄ = 1 / Σk(Yk/Wk) を使用。

テスト (3 テスト、全 PASS)

数値的注意事項

• T* は最小値 0.1 にクランプ（超低温での発散を防ぐ）
• sum_jk = 0 の場合（希薄モル分率など）、Dkm = 0 を返す
• Wilke 則でモル分率が < 1e-20 の成分はスキップ
• 二成分拡散係数の計算で Dkj < 1e-30 の場合はクランプして除算を防ぐ

参考文献

• Neufeld, Janzen, Aziz (1972). J. Chem. Phys. 57, 1100.
• Monchick, Mason (1961). J. Chem. Phys. 35, 1676.
• Mason, Monchick (1962). J. Chem. Phys. 36, 1622.
• Bird, Stewart, Lightfoot (2002). Transport Phenomena, 2nd ed. Wiley. §9.3.
• Wilke (1950). J. Chem. Phys. 18, 517.

flame モジュール

1D 火炎ソルバーのコア層。格子・解ベクトル・残差・ヤコビアン・格子精細化・ ソルバー制御を担当する。

モジュール構成

flame/
├── mod.rs            — サブモジュールの公開宣言
├── domain.rs         — 非一様 1D 格子 (Grid)
├── state.rs          — 解ベクトルのレイアウトとインデックス
├── residual.rs       — 残差評価 F(x)
├── jacobian.rs       — 数値ヤコビアン
├── refine.rs         — GRAD/CURV 適応格子精細化
└── solver_driver.rs  — 火炎計算の上位制御

domain.rs — 格子

Grid

非一様 1D 格子。格子点座標を Vec<f64> で保持する。

state.rs — 解ベクトルのレイアウト

解ベクトル構造

x = [T_0, Y_{0,0}, …, Y_{K-1,0},  T_1, Y_{0,1}, …,  …,  T_{J-1}, …, Y_{K-1,J-1},  M]

• 長さ: (K+1) × J + 1
• 各格子点に NATJ = K+1 変数（温度 1 + 化学種 K）
• 末尾の 1 要素は質量流束 M [kg/(m²·s)]（固有値）

FlameState<'a> / FlameStateMut<'a>

解ベクトルへの高レベルアクセサ。

#![allow(unused)]
fn main() {
state.temperature(j)      // T[j]
state.species(k, j)       // Y_k[j]
state.mass_flux()         // M
state.y_slice(j)          // &Y[0..K] at point j
}

initial_guess

シグモイドプロファイルによる初期推定値の生成。

T(z) = T_u + σ((z - z_c)/z_w) × (T_b - T_u)
Y(z) = Y_u + σ((z - z_c)/z_w) × (Y_b - Y_u)

initial_guess_from_csv(path, mech, grid)

Cantera などの外部ソルバーが出力した CSV から初期推定値を読み込む。

• 必須列: z [m], T [K], M [kg/m2/s], Y_<name> (各化学種)
• CSV の z 座標から現在の格子に線形補間し、Y は正規化（和を 1 に）する
• M は CSV 全点の平均値を使用
• solver 設定の initial_profile に CSV パスを指定すると疑似過渡継続（Phase 1）がスキップされ、Newton 法から直接開始できる

residual.rs — 残差評価

支配方程式

内部点 j = 1 … J-2:

化学種方程式（k = 0 … K-2）:

F_yk = M × (Y_k,j - Y_k,j-1)/Δz_m
     + (jk_{j-1/2} - jk_{j-3/2}) / Δz_av
     - ωk × Wk

最終化学種は和の閉鎖式:

F_yK = ΣYk - 1

エネルギー方程式:

F_T = -M × cp × (T_j - T_j-1)/Δz_m
    + (λ_{j+1/2}×(T_{j+1}-T_j)/Δz_p - λ_{j-1/2}×(T_j-T_j-1)/Δz_m) / Δz_av
    - Σk (jk_{j-1/2} + jk_{j+1/2})/2 × cpk × (T_j - T_j-1)/Δz_m
    - Σk ωk × hk

左境界 (j = 0):

F_T  = T_0 - T_u
F_yk = M × (Y_k,0 - Y_k,u) + jk_{1/2}    (inlet flux BC)

右境界 (j = J-1):

F_T  = T_{J-1} - T_{J-2}    (zero gradient)
F_yk = Y_k,J-1 - Y_k,J-2

固有値閉鎖:

F_M = T(j_fix) - T_fix

FlameConfig

eval_residual(x, rhs, mech, grid, config, x_old, rdt)

• 定常計算: x_old = None, rdt = 0.0
• 疑似過渡継続: x_old と rdt = 1/Δt を渡すと内部点に backward-Euler 項を加算:

  T 方程式:  + rdt × ρ × cp × (T - T_old)    [W/m³]
  Yk 方程式: + rdt × ρ × (Yk - Yk_old)        [kg/(m³·s)]
  

  密度・熱容量のスケーリングにより各方程式が次元整合する。

拡散フラックス

混合平均拡散係数によるモル分率基準の拡散フラックスに補正速度を適用して ΣJk = 0 を強制 (PREMIX / Cantera Flow1D と同一の形式):

jk_raw = -ρ × (Wk/W̄) × Dkm × dXk/dz
jk     = jk_raw - Yk × Σj jk_raw,j     (左セル Yk で補正速度を重み付け)

質量分率 Yk ではなくモル分率 Xk の勾配を使うのは PREMIX の定式化に準拠。 W̄ は混合平均分子量。輸送係数はミッドポイント (j, j+1) の平均状態で評価する。

単位系

• 濃度は mol/cm³ で production_rates に渡す（CGS、A値と整合）
• ωk の出力 mol/(cm³·s) × 1e6 → mol/(m³·s)、さらに × Wk で kg/(m³·s)

テスト (2 テスト、全 PASS)

jacobian.rs — 数値ヤコビアン

numerical_jacobian → (BandedMatrix, Vec<f64>, Vec<f64>)

前進差分による数値ヤコビアン J = ∂F/∂x を帯行列 + M 列ベクトルの形で返す。

戻り値:

• J_int: (n-1)×(n-1) 帯行列（内部変数のブロック; M 行・M 列を除く）
• m_col: 各残差の ∂F/∂M（長さ n-1 の列ベクトル）
• jfc: M 残差の ∂F_M/∂x（長さ n-1 の行ベクトル）

これにより bordered_solve_m が 1 回の帯 LU でシュア補完を通じて δM を含む完全な Newton ステップを計算できる（密行列不要、O(n·kl) コスト）。

パラメータ:

• 帯域幅: kl = ku = 2×NATJ（隣接ブロック間結合）
• 摂動幅: h = √(2ε_machine) × max(|x_j|, 1)
• ローカルステンシル eval_residual_range で M 列を除く列を効率的に構築

bordered_solve_m

枠付き帯行列系 [[J_int, m_col], [jfc, jff]] · [δx, δM]ᵀ = [−F_int, −F_M] を解く。

1. J_int·u = m_col を帯 LU で解いて補助ベクトル u を求める
2. シュア補完 S = jff − jfc·u で δM を決定
3. δx_int = J_int⁻¹(−F_int − m_col·δM) を帯後退代入で求める

テスト (2 テスト、全 PASS)

refine.rs — 適応格子精細化

RefineCriteria

TOML で [grid] セクションに記述することで上書き可能。

精細化判定 (Cantera refine.cpp 準拠)

区間 j に中点を挿入する条件（T と全化学種 Yk について評価）:

GRAD:  |φ_{j+1} - φ_j|          > grad × (φ_max - φ_min)
CURV:  |slope[j+1] - slope[j]|  > curv × (slope_max - slope_min)
RATIO: dz[j] > ratio × dz[j-1]  または  dz[j-1] > ratio × dz[j]

• GRAD: 値の変化量がドメイン全体のレンジに対して大きすぎる区間を検出。
• CURV: 一階微分の変化量（slope[j] = dφ/dz）がドメイン全体の傾き範囲に対して大きすぎる区間を検出。両辺とも単位は φ/m で次元整合。
• RATIO: 隣接セル幅の比が ratio を超える場合に挿入。炎帯で局所的に細かくなった格子が隣の粗いセルへ連鎖的に伝播するため、滑らかな格子遷移を保証する。

find_refinement_points(x, mech, grid, criteria) → Vec<usize>

挿入すべき区間インデックスのリストを返す。refine から内部で呼ばれるが単体でも公開。

refine(x, mech, grid, criteria) → Option<(Grid, Vec<f64>)>

精細化が必要な点がなければ None。 挿入は逆順で行いインデックスの安定性を保つ。新しい格子点の解は線形補間で生成する。 M は末尾要素として一旦取り出し、全挿入後に再付加する（インデックスのズレを防ぐ）。

テスト (2 テスト、全 PASS)

solver_driver.rs — ソルバー制御

run_flame(config) の処理フロー

1. 機構ファイルの読み込み (現在は Cantera YAML のみ; CHEMKIN-II は未接続)
2. 未燃・既燃組成の計算（当量比 φ モード または 直接組成指定モード）
3. 断熱火炎温度 T_ad の推定（エンタルピー収支 Newton 法）
4. ラジカルシード: 既燃組成に H/OH/H₂ を注入（chain-branching 起動）
5. 初期格子の生成（均一、initial_points 点）
6. Phase 1: 疑似過渡継続 (n_steps, dt_initial)
     ※ initial_profile が指定されている場合はスキップ
       → Cantera 収束解を初期値とする場合などに使用
7. Phase 2: Newton 法 + 適応格子精細化のループ（最大 20 パス）
   a. Newton 法を実行し ‖F‖ を計測
   b. 収束判定: ‖F‖ < atol*√n + rtol*‖F‖  または  ‖F‖ が 10% 以上改善
   c. 改善不十分なら PT を追加して再試行
   d. GRAD/CURV/RATIO 基準で格子精細化
   e. 精細化不要 (refine が None) になるまで繰り返す
   f. 各パス後に z_fix を T プロファイルから再計算
8. CSV 出力 + サマリー表示

組成計算: compute_compositions

2 つのモードに対応:

当量比モード (fuel + oxidizer + equivalence_ratio):

n_ox = stoich_O2 / (x_O2_in_ox × φ)
x[k] = (x_fuel[k] + n_ox × x_ox[k]) / (1 + n_ox)

化学量論 O2 量は完全燃焼式 (C + H/4 − O/2 + S) で決定。

直接指定モード (composition モル分率マップ): 指定値を正規化して未燃混合気として使用。

どちらのモードも最後にモル分率 → 質量分率に変換する。

既燃組成: estimate_burned_composition

元素バランスによる完全燃焼計算:

全 C → CO2  (O 不足なら一部 CO)
全 H → H2O  (O 不足なら一部 H2)
残余 O → O2
N → N2, S → SO2, Ar/He はそのまま

リッチ火炎 (φ > 1) では O2 枯渇時に CO/H2 を生成する。あくまで初期推定用。

ラジカルシード: seed_radicals

完全燃焼組成にはラジカルが含まれないため chain-branching が開始できない。 以下の 2 反応で H₂O を部分解離させてラジカルを注入する:

H₂O → H + OH          (f1 = 6%、H₂O 質量を H:OH = 1:17 で分配)
2H₂O → 2H₂ + O₂       (f2 = 2%、質量保存)

H₂/air (φ=1, T_ad ≈ 2500 K) の平衡ラジカル濃度 (x_H ≈ 0.6%, x_OH ≈ 3.8%) に近似。 注入後に ΣYk = 1 に再正規化する。

初期推定パラメータ

su_initial_guess は TOML の [solver] セクションで上書き可能。

収束後に z_fix を動的追跡: find_z_fix

各 Newton パス後に T が t_fix を最初に越える位置を線形補間で特定し、 res_config.z_fix を更新する。格子精細化で格子点位置が変わっても固有値制約が 安定して機能する。

格子精細化の収束挙動

H2/air φ=1 の典型的な実行例（Cantera CSV 初期値、grad=0.05, curv=0.10, ratio=2.0）:

最終 Su = 2.3327 m/s（Cantera 参照値 2.3354 m/s、誤差 0.1%）。

solver モジュール

線形系の求解・Newton 法・疑似過渡継続を担当する。

モジュール構成

solver/
├── mod.rs               — サブモジュールの公開宣言
├── banded.rs            — 帯行列 LU 分解・求解
├── dense.rs             — 密行列 LU 分解・求解 (faer)
├── newton.rs            — 減衰 Newton 法（直線探索付き）
└── pseudo_transient.rs  — 疑似過渡継続 (PT)

banded.rs — 帯行列

BandedMatrix

LAPACK 互換の一般帯行列 (GB) 形式で係数行列を保持する。

格納レイアウト: data[col × ldab + (kl + ku + row − col)]

• ldab = 2×kl + ku + 1（部分ピボッティングの fill-in 分の kl 行を余分に確保）

fill-in の扱い: ピボット行入れ替え後、行 k の非ゼロ範囲は k+ku+kl 列まで 拡張しうる。後退代入でも同じ幅を使用することで fill-in を正確に処理する。

テスト (3 テスト、全 PASS)

dense.rs — 密行列

DenseMatrix

faer の部分ピボッティング LU で n×n 密行列を求解する。

帯行列では帯外に落ちる質量流束 M 列（解ベクトル末尾、全残差に結合）を 正確に扱うために使用する。

テスト (2 テスト、全 PASS)

newton.rs — 減衰 Newton 法

アルゴリズム概要

for iter in 0..max_iter:
    評価: F(x)
    収束判定: ‖F‖ < atol·√n + rtol·‖x‖
    Jacobian: (J_int, m_col, jfc) = numerical_jacobian(x, F)  ← 最大 max_jac_age 反復再利用
    線形系:   枠付き帯行列系を bordered_solve_m で解く
                δx_int = J_int⁻¹(-F_int + F_M / (jfc - m_col·J_int⁻¹·e_M) × m_col)
                δM    = シュア補完で決定
    ステップクリップ: 物理スケール上限を超える成分を縮小
    直線探索: α を 1 から 0.5 刻みで最大 12 回試行
    更新:     x ← x + α·step
    物理クランプ: T, Y, M を有効範囲に制限
    失速検出: best ‖F‖ が 1% 以上改善しない反復が max_stall_iter 回続いたら受理

帯行列 + 枠付きシュア補完方式: numerical_jacobian が返す (n-1)×(n-1) 帯行列 J_int、 M 列ベクトル m_col、M 行ベクトル jfc を使って bordered_solve_m が 1 回の帯 LU で δM を含む完全な Newton ステップを求める。密 Jacobian 不要で O(n·kl) のコスト。

NewtonConfig

Jacobian の強制再計算: ‖F‖ が前回評価時から 10 倍以上改善した場合、 線形化点が大きく変わったとみなして jac_age をリセットする。

ステップクリップ (step_clip_factor)

物理クランプ (clamp_physical)

テスト (2 テスト、全 PASS)

pseudo_transient.rs — 疑似過渡継続

目的

Newton 法の収束域外（点火前の急勾配・スティッフ化学反応領域）から Newton 域へ 誘導する。TWOPNT スタイルの後退オイラー法を採用。

アルゴリズム概要

for step in 0..n_steps:
    評価: F_steady(x)
    ‖F_steady‖ < 1e7 → Newton 域に達したので終了
    dt < 1e-13 → 収束断念して Newton へ

    rdt = 1/dt
    x_old = x

    内側 Newton (max_inner_iter 回):
        F_pt(x) = F_steady(x) + D/dt·(x − x_old)
        J_pt    = J_steady(x) + D/dt の対角（ρ·cp で T、ρ で Yk、1 で M）
        M 列をゼロクリア（M の疑似時間項は対角のみ）
        step = −J_pt⁻¹·F_pt
        ステップクリップ + 直線探索 + 物理クランプ

    外側受理: ‖F_steady(x_new)‖ ≤ 1.1×‖F_steady(x_old)‖
        成功 → dt *= dt_grow (最大 dt_max まで)
        失敗 → x = x_old, dt *= 0.5; 30 回連続失敗で中断

PseudoTransientConfig

設計上の注意

直線探索の初期化: 内側直線探索の best_norm を ∞ で初期化することで、 J_pt が不定値（D/dt < |λ_chem|）のとき全 α で F_pt が増加しても α = 0 に停滞しない。最小の有限 F_pt を与える α を返す。

外側受理の 10% マージン: 化学反応が非常に剛な点火フェーズでは、 良い内側 Newton ステップでも F_steady が一時的に増加することがある。 1.1 倍の許容マージンにより探索的なステップを通過させ、Newton 域への 到達を可能にする。

D/dt 対角: ρ·cp (T 方程式)、ρ (Yk 方程式)、1 (M) で各変数の 次元が整合した疑似時間スケーリングを実現する（residual.rs の PT 埋め込みと同一）。

テスト (2 テスト、全 PASS)

io モジュール

TOML 設定ファイルの読み込みと CSV 出力を担当する。

モジュール構成

io/
├── mod.rs     — サブモジュールの公開宣言
├── input.rs   — TOML 設定ファイルのデシリアライズ・検証
└── output.rs  — CSV 出力・サマリー表示

input.rs — 設定ファイル

FlameConfig::from_file(path)

TOML ファイルを読み込み、値の検証を行って FlameConfig を返す。 検証エラーは anyhow::Error として返る。

設定構造体

[mechanism]

[flame]

未燃組成は 当量比モード と 直接指定モード のいずれか一方を使用する。 両方同時に指定するとバリデーションエラーになる。

[grid]

[solver]

[output]

バリデーション規則

TOML 例 — 当量比モード

[mechanism]
file = "data/h2o2.yaml"

[flame]
pressure = 101325.0
t_unburned = 300.0
domain_length = 0.02
equivalence_ratio = 1.0

[flame.fuel]
H2 = 1.0

[flame.oxidizer]
O2 = 0.21
N2 = 0.79

[grid]
initial_points = 20
max_points = 500
grad = 0.05
curv = 0.10
ratio = 2.0

[solver]
atol = 1.0e-9
rtol = 1.0e-6
max_newton_iter = 50
time_steps = 100
dt_initial = 1.0e-7
su_initial_guess = 2.0        # 省略可: デフォルト 0.4 m/s
# initial_profile = "data/cantera_h2air_initial.csv"  # 省略時はシグモイド初期推定

[output]
file = "flame_solution.csv"

TOML 例 — 直接組成指定モード

[flame]
t_unburned = 300.0
composition = { H2 = 0.3, O2 = 0.15, N2 = 0.55 }

テスト (5 テスト、全 PASS)

output.rs — CSV 出力

write_csv(path, x, mech, grid, pressure)

収束解ベクトル x を CSV ファイルに書き出す。

列構成 (計 6 + 2K 列、K = 化学種数):

M [kg/m2/s] 列は initial_guess_from_csv が読み込む列と同一のフォーマットであるため、 本ソルバーの出力 CSV をそのまま次の計算（当量比スイープ等）の initial_profile として使用できる。

数値フォーマット: z, M, u, rho, hrr, Xk, Yk は {:.6e}、T は {:.4}。

print_summary(x, mech, grid, pressure)

標準出力にサマリーを表示する。

---------------------------------------------------
  Laminar flame speed Su = X.XXXX m/s
  Max temperature       = XXXX.X K
  Grid points           = XXX
  Mass flux M           = X.XXXXe-X kg/(m²·s)
---------------------------------------------------

• Su = M / ρ(j=0) — 左境界の密度を使って質量流束から火炎速度に換算

テスト (2 テスト、全 PASS)