中央差分法（陽解法）
理論と物理
陽解法とは
連立方程式を解かない？ それは速そうです。
1ステップあたりの計算は桁違いに速い。ただし安定条件があり、時間刻みが非常に小さい（$\Delta t \propto L_{min} / c$、$c$ は音速）。衝撃のような短時間現象に最適。
中央差分法のアルゴリズム
運動方程式 $[M]\{\ddot{u}\} = \{F\} - [K]\{u\} - [C]\{\dot{u}\}$ に対して：
1. 加速度の計算: $\{\ddot{u}_n\} = [M]^{-1}(\{F_n\} - \{F_{int,n}\})$
2. 速度の更新: $\{\dot{u}_{n+1/2}\} = \{\dot{u}_{n-1/2}\} + \Delta t_n \{\ddot{u}_n\}$
3. 変位の更新: $\{u_{n+1}\} = \{u_n\} + \Delta t_{n+1/2} \{\dot{u}_{n+1/2}\}$
$[M]^{-1}$ が対角行列なら、逆行列は各対角成分の逆数を取るだけ…計算が超高速！
だから陽解法では集中質量マトリクス（対角）が必須。一貫質量（非対角）を使うと逆行列の計算にコストがかかり、陽解法の利点が失われる。
安定条件（CFL条件）
陽解法の安定時間増分：
$$ \Delta t \leq \frac{L_{min}}{c} = \frac{L_{min}}{\sqrt{E/\rho}} $$
陽解法の安定時間増分：
$L_{min}$ は最小要素サイズ、$c$ は材料の音速（弾性波速度）。
鋼（$c \approx 5000$ m/s）で $L_{min} = 5$ mm なら $\Delta t \leq 1 \mu s$…非常に小さい！
だから1秒間の解析に100万ステップ以上必要。衝撃（数ms）なら数千ステップで済むが、準静的（数秒）だと膨大。陽解法は短時間の動的現象に適する。
陽解法 vs. 陰解法
| 特性 | 陽解法 | 陰解法 |
|---|---|---|
| 1ステップの計算 | 非常に軽い | 重い（連立方程式） |
| 時間刻み | 非常に小さい（$\mu s$オーダー） | 大きい（$ms \sim s$） |
| 安定条件 | あり（CFL条件） | なし（無条件安定） |
| 適用場面 | 衝撃、爆発、衝突 | 準静的、振動 |
| 接触 | 容易 | 困難 |
| 非線形 | 自然に扱える | 反復が必要 |
接触が容易なのは大きな利点ですね。
陽解法では接触の検出と力の計算が各ステップで独立に行われ、収束問題がない。車の衝突のように多数の接触が同時に起きる問題では陽解法が圧倒的に有利。
まとめ
中央差分法（陽解法）を整理します。
要点：
- 連立方程式を解かない — $[M]^{-1}$ が対角なら逆数だけ
- $\Delta t \leq L_{min} / c$ — CFL条件。最小要素で決まる
- 短時間の動的現象に最適 — 衝撃、爆発、衝突
- 接触が容易 — 収束問題なし
- 集中質量マトリクスが必須 — 対角行列でないと速くならない
陽解法は「短時間の激しい現象」のためのアルゴリズムなんですね。
LS-DYNAの全解析が陽解法。Abaqus/Explicitも陽解法。自動車の衝突安全はこの手法なしには成り立たない。
中心差分法は1695年の微積分から生まれた
Leibnizが1695年に発表した有限差分の概念を動力学に適用した中心差分法は、加速度をa(t)＝[u(t+Δt)−2u(t)＋u(t−Δt)]/Δt²で近似する。この2次精度スキームはvon Neumann安定性解析から臨界時間刻みΔtcr＝2/ωmaxが導かれ、最小要素の固有振動数を超えないよう制約される。1970年代のLS-DYNAの前身コードDYNA3DでこのアルゴリズムがFEMに初めて大規模実装された。
各項の物理的意味
- 慣性項（質量項）：$\rho \ddot{u}$、つまり「質量×加速度」。急ブレーキで体が前に投げ出された経験はありませんか？ あの「持っていかれる感じ」がまさに慣性力です。重い物体ほど動き出しにくく、動き出したら止まりにくい。地震で建物が揺れるのも、地面が急に動いたのに建物の質量が「置いていかれる」から。静解析ではこの項をゼロにしますが、それは「ゆっくり力をかけるから加速度は無視できる」という仮定です。衝撃荷重や振動問題では絶対に省略できません。
- 剛性項（弾性復元力）：$Ku$ や $\nabla \cdot \sigma$。ばねを引っ張ると「戻ろうとする力」を感じますよね？ あれがフックの法則 $F=kx$ であり、剛性項の本質です。では質問——鉄の棒とゴム紐、同じ力で引っ張るとどちらが伸びるでしょうか？ 当然ゴムです。この「伸びにくさ」がヤング率 $E$ であり、剛性を決めます。よくある勘違い：「剛性が高い＝強い」ではありません。剛性は「変形しにくさ」、強度は「壊れにくさ」で、別の概念です。
- 外力項（荷重項）：体積力 $f_b$（重力など）と表面力 $f_s$（圧力、接触力など）。こう考えてみてください——橋の上のトラックの重さは「中身全体にかかる力」（体積力）、タイヤが路面を押す力は「表面だけにかかる力」（表面力）。風圧、水圧、ボルトの締付力…すべて外力です。ここでありがちな失敗：荷重の方向を間違える。「引張」のつもりが「圧縮」になっていた——笑い話に聞こえますが、3D空間で座標系が回転していると実際に起こります。
- 減衰項：レイリー減衰 $C\dot{u} = (\alpha M + \beta K)\dot{u}$。ギターの弦を弾いてみてください。音は鳴り続けますか？ いいえ、徐々に小さくなりますよね。振動エネルギーが空気抵抗や弦の内部摩擦で熱に変わるからです。車のショックアブソーバーも同じ原理——わざと振動エネルギーを吸収して乗り心地を良くしています。もし減衰がゼロだったら？ 建物は地震の後いつまでも揺れ続けることになります。実際にはそうならないので、適切な減衰の設定が重要です。
仮定条件と適用限界
次元解析と単位系
| 変数 | SI単位 | 注意点・換算メモ |
|---|---|---|
| 変位 $u$ | m（メートル） | mm入力時は荷重・弾性率もMPa/N系に統一すること |
| 応力 $\sigma$ | Pa（パスカル）= N/m² | MPa = 10⁶ Pa。降伏応力との比較時に単位系の不一致に注意 |
| 歪み $\varepsilon$ | 無次元（m/m） | 工学歪みと対数歪みの区別に注意（大変形時） |
| 弾性率 $E$ | Pa | 鋼: 約210 GPa、アルミ: 約70 GPa。温度依存性に注意 |
| 密度 $\rho$ | kg/m³ | mm系ではtonne/mm³（= 10⁻⁹ tonne/mm³ for 鋼） |
| 力 $F$ | N（ニュートン） | mm系ではN、m系ではNで統一 |
数値解法と実装
陽解法の実装
陽解法のFEMコードはどう実装されていますか？
メインループは極めてシンプル：
```
for each time step:
1. 内力の計算: F_int = assemble(sigma B V) $ 要素ループ
2. 外力の計算: F_ext = boundary_conditions()
3. 加速度: a = (F_ext - F_int) / M_diag $ 対角割り算
4. 速度更新: v = v + dt * a
5. 変位更新: u = u + dt * v
6. 接触チェック: contact_detection_and_force()
7. 時間増分更新: dt = CFL * L_min / c
```
ステップ3の「対角割り算」がポイントですね。
そう。ここが陽解法の速さの源泉。陰解法ではステップ3が「$n \times n$ の連立方程式の求解」になり、計算量が桁違い。
安定時間増分の管理
安定時間増分は全要素の中で最も小さい値で決まる：
1つの要素が小さいだけで全体の $\Delta t$ が小さくなる。
1つの悪い要素のせいで全体が遅くなる？
メッシュ品質が計算速度に直結する。アスペクト比の悪い要素やつぶれた要素は $L_e$ が非常に小さく、$\Delta t$ を減少させる。陽解法ではメッシュ品質が計算効率の鍵。
質量スケーリング
小さい要素の質量を人為的に増やして $\Delta t$ を大きくする質量スケーリング：
$\rho$ を増やすと $c$ が下がり $\Delta t$ が大きくなる。準静的問題の陽解法で必須のテクニック。
質量を増やすと慣性効果が変わりませんか？
変わる。だから運動エネルギー/内部エネルギー < 5%を確認する必要がある。この条件が満たされれば「準静的」と見なせる。
ソルバー
| ソルバー | 陽解法 |
|---|---|
| LS-DYNA | 陽解法が主力。衝突安全の業界標準 |
| Abaqus/Explicit | Abaqus StandardのExplicit版 |
| Ansys LS-DYNA | Ansysに統合されたLS-DYNA |
| PAM-CRASH | ESIの衝突解析ソルバー |
| RADIOSS | Altairの陽解法ソルバー |
まとめ
陽解法の実装を整理します。
要点：
- メインループが極めてシンプル — 対角行列の逆数で加速度計算
- 最小要素が$\Delta t$を支配 — メッシュ品質が計算速度の鍵
- 質量スケーリングで準静的にも適用 — 運動エネルギー < 5%で確認
- LS-DYNAが衝突安全の業界標準 — 自動車OEM全社が使用
- Abaqus/ExplicitはAbaqusファミリーの陽解法 — Standard→Explicitの切り替え
質量スケーリングで計算を10倍速く
陽解法の時間刻みは最小要素寸法に支配されるため、数個の微小要素が全体の計算コストを100倍以上高める場合がある。LS-DYNAの質量スケーリング（MASS SCALING）は微小要素に仮想質量を付加し時間刻みを引き上げる技法で、追加質量が全体質量の1〜2%未満であれば解への影響は許容範囲内とされている。ただし慣性力が支配的なフラッシュ衝突では3%超で誤差が顕著になる。
線形要素（1次要素）
節点間を線形補間。計算コストは低いが、応力の精度が低い。せん断ロッキングに注意（低減積分やB-bar法で緩和）。
2次要素（中間節点付き）
曲線的な変形を表現可能。応力精度が大幅に向上するが、自由度は約2〜3倍に増加。推奨：応力評価が重要な場合。
完全積分 vs 低減積分
完全積分：過剰拘束（ロッキング）のリスク。低減積分：アワーグラスモード（零エネルギーモード）のリスク。適材適所で選択。
アダプティブメッシュ
誤差指標（ZZ推定量等）に基づく自動細分化。応力集中部の精度を効率的に向上。h法（要素分割）とp法（次数増加）がある。
ニュートン・ラフソン法
非線形解析の標準的手法。接線剛性マトリクスを毎反復更新。収束半径内で2次収束するが、計算コストが高い。
修正ニュートン・ラフソン法
接線剛性マトリクスを初期値または数反復毎に更新。各反復のコストは低いが、収束速度は線形的。
収束判定基準
力の残差ノルム: $||R|| / ||F_{ext}|| < \epsilon$（一般に $\epsilon = 10^{-3}$〜$10^{-6}$）。変位増分ノルム: $||\Delta u|| / ||u|| < \epsilon$。エネルギーノルム: $\Delta u \cdot R < \epsilon$
荷重増分法
全荷重を一度に負荷せず、小刻みに増加させる。弧長法（Riks法）は荷重-変位関係の極値点を越えて追跡可能。
直接法 vs 反復法のたとえ
直接法は「連立方程式を筆算で正確に解く」方法——確実だが大規模問題では時間がかかりすぎる。反復法は「当て推量を繰り返して正解に近づく」方法——最初は大雑把な答えだが、反復するたびに精度が上がる。辞書で言葉を探すとき、最初のページから順番に探す（直接法）より、見当をつけて開き、前後に調整する（反復法）方が効率的なのと同じ原理。
メッシュの次数と精度の関係
1次要素は「定規で曲線を近似する」——直線の折れ線で表現するため精度に限界がある。2次要素は「フレキシブルカーブ」——曲線的な変化を表現でき、同じメッシュ密度でも格段に精度が向上する。ただし、1要素あたりの計算コストは増えるため、トータルのコスト対効果で判断する。
実践ガイド
陽解法の実務適用
陽解法は実務でどう使いますか？
エネルギーバランスの確認
陽解法の結果の妥当性確認で最も重要なのがエネルギーバランス：
チェック項目：
- 全エネルギー — 保存されているか（数%以内の変動）
- アワーグラスエネルギー — 全エネルギーの5%以下
- 運動エネルギー（準静的の場合） — 全エネルギーの5%以下
- 接触エネルギー — 負にならないか（貫通の兆候）
エネルギーバランスが合っていれば結果は信頼できる？
必要条件だが十分条件ではない。エネルギーが保存されていても、変形パターンが非物理的な場合がある。変形の可視化も必ず行う。
実務チェックリスト
「エネルギーバランスの確認」が陽解法の最重要チェックですね。
エネルギーは「嘘をつかない」。何か問題があれば必ずエネルギーバランスに現れる。全ての陽解法解析でエネルギー出力を確認する習慣をつけよう。
工場での板金プレスにも陽解法が活躍
自動車ボンネットのプレス成形解析では、金型速度を実際の0.5m/sから解析上5m/sに加速する「工程加速法」と質量スケーリングを組み合わせ、Autoformなどのインクリメンタル陽解法コードで1工程あたり15〜30分の計算を実現している。ただし加速比が10倍を超えると動的効果でスプリングバック量が実測と乖離するため、Audi・BMWは加速比5倍を社内基準として規定している。
解析フローのたとえ
解析の流れは、実は料理とそっくりです。まず材料を買い出し（CADモデルの準備）、下ごしらえをして（メッシュ生成）、火にかけて（ソルバー実行）、最後に盛り付ける（後処理で可視化）。ここで大事な問いかけ——料理で一番失敗しやすい工程はどこでしょう？ 実は「下ごしらえ」なんです。メッシュの品質が悪いと、どんなに優秀なソルバーを使っても結果はめちゃくちゃになります。
初心者が陥りやすい落とし穴
あなたはメッシュ収束性を確認していますか？ 「計算が回った＝結果が正しい」と思っていませんか？ これ、実はCAE初心者が最も陥りやすい罠です。ソルバーは与えられたメッシュで「それなりの答え」を必ず返します。でもメッシュが粗すぎれば、その答えは現実から大きくずれている。最低3段階のメッシュ密度で結果が安定することを確認する——これを怠ると「コンピュータが出した答えだから正しいはず」という危険な思い込みに陥ります。
境界条件の考え方
境界条件の設定は、試験の「問題文を書く」のと同じです。問題文が間違っていたら？ どんなに正確に計算しても答えは間違いますよね。「この面は本当に完全固定なのか」「この荷重は本当に一様分布なのか」——現実の拘束条件を正しくモデル化することが、実は解析全体で最も重要なステップだったりします。
ソフトウェア比較
陽解法のツール
陽解法のソルバーを比較してください。
LS-DYNAが圧倒的ですか。
自動車の衝突安全ではLS-DYNAが世界中のOEM/Tier1で標準。EuroNCAP、FMVSS等の規格に対応したモデル（WorldSID、THOR等のダミーモデル）がLS-DYNA向けに提供されている。
選定ガイド
「衝突安全 = LS-DYNA」は動かないんですね。
LS-DYNAの30年以上の実績とダミーモデルの蓄積は他社に真似できない。衝突安全の新規参入には既存のLS-DYNAインフラとの互換性が必須。
陽解法コードの起源はすべてLLNLにあり
現代の主要陽解法コード（LS-DYNA・PAM-CRASH・Radioss）はすべてLawrence Livermore国立研究所のDYNA3D（1976年）を起源に持つ。Radiossは1984年にRegieが欧州向けに独自拡張したもので、後にAltairが2006年に買収。PAM-CRASHはESI Groupが1980年代に自動車向けに特化して開発した。起源が同じだけに基本アルゴリズムは共通だが、接触処理・材料モデルの実装に各社の差別化がある。
先端技術
陽解法の先端研究
陽解法の最前線を教えてください。
GPU加速
陽解法は各要素の計算が独立（大規模並列可能）だから、GPU加速との相性が非常に良い。LS-DYNAのGPU版は特定の要素タイプ（HEX8, SHELL）でCPUの10〜100倍の高速化を報告。
100倍！ 衝突解析が1時間から1分に…。
現実にはI/OやCPU-GPU間のデータ転送がボトルネックで、実効的には5〜30倍程度。それでも設計サイクルの大幅短縮になる。
陽-陰連成
構造の一部は陽解法（衝撃部）、他の部分は陰解法（準静的）で同時に計算する陽-陰連成。計算効率を最大化。LS-DYNAやAbaqusの一部で対応。
IGA（等幾何解析）+ 陽解法
NURBS/T-spline基底のIGA要素を陽解法に適用する研究。LS-DYNAにIGAシェルが実装済み。CAD形状の正確な表現とメッシュ不要の利点。
まとめ
NURBS/T-spline基底のIGA要素を陽解法に適用する研究。LS-DYNAにIGAシェルが実装済み。CAD形状の正確な表現とメッシュ不要の利点。
サブサイクリングで異スケール問題を効率化
流体-構造連成や異種材料接触問題では、構造部と流体部で適切な時間刻みが数十倍異なる場合がある。サブサイクリング（またはマルチタイムステップ）技法は、粗メッシュ領域を大きな刻みで、細メッシュ領域を小さな刻みで計算し、境界で情報交換する。LS-DYNA R12以降のALE-構造連成解析でこの機能が改善され、計算時間を最大60%短縮できると公式ドキュメントに記載されている。
トラブルシューティング
陽解法のトラブル
陽解法でよくあるトラブルを教えてください。
エネルギーが保存されない
全エネルギーが時間とともに増加します。
接触の貫通がエネルギーを生成している。対策：
- 接触のペナルティ剛性を上げる
- 接触面のメッシュを細かくする
- 接触アルゴリズムを変更（NTS→MORTAR等）
アワーグラスモード
低減積分要素（HEX8R, SHELL Q4R）でアワーグラスモードが励起。
対策：
- アワーグラスエネルギーを監視（< 5%）
- アワーグラス制御のタイプを変更（粘性型→剛性型）
- 完全積分要素に切り替え（計算コスト増）
時間増分が極端に小さい
1つの小さい要素が全体の$\Delta t$を制限。
対策：
- メッシュ品質を改善（小さい要素を除去/拡大）
- 質量スケーリング（対象要素のみ）
- 問題の要素をマージ or 削除
計算が飛ぶ（負の体積）
要素が過度に変形して体積がゼロまたは負になる。
対策：
- 要素削除（*MAT_ADD_EROSION等）で過変形要素を除去
- メッシュを細かくして変形を分散
- 材料モデルの破壊基準を追加
まとめ
陽解法のトラブル対処、整理します。
エネルギーを常に監視する。これが陽解法の鉄則ですね。
エネルギーは「宇宙の帳簿」。帳簿が合わなければどこかに問題がある。陽解法では全ての問題がエネルギーバランスに反映される。
負のエネルギーはマイナス剛性要素の警告
陽解法解析中に内部エネルギーが急激に負になる場合、破壊・削除された要素の残留剛性か、接触ペナルティによる過大ひずみが原因のことが多い。LS-DYNAの*CONTROL_ENERGYでSLIDEOPT＝2（接触エネルギー追跡）を有効にし、どの接触ペアでエネルギー収支が崩れているかをRCFORC出力で特定するのが定石的なデバッグ手順である。
関連トピック
なった
詳しく
報告