ラメの厚肉円筒解 — 内圧問題のFEM検証ベンチマーク
理論と物理
概要 — なぜラメ解がFEM検証に最適か
ラメの公式って圧力容器の教科書で見ましたけど、FEM検証にも使うんですか？
使うどころか、最良の軸対称ベンチマークと言ってもいい。内圧を受ける厚肉円筒の応力分布——半径方向応力 $\sigma_r$ と周方向応力 $\sigma_\theta$ ——の解析解が閉じた形で既知だから、FEAの数値解と「一発で答え合わせ」ができるんだ。
他のベンチマーク問題、例えば片持ち梁の曲げとかと比べて、何が優れているんですか？
ポイントは3つある：
- 応力勾配が急 — 内面で応力が最大になり外面に向かって急減衰する。これはメッシュの品質と要素次数の差が如実に出る条件だ
- 2次要素なら粗メッシュでも誤差0.1%以内 — 4要素分割のQUAD8軸対称要素で既に0.12%。8分割なら事実上の完全一致になる
- 実務との直結 — ASME Section VIIIの圧力容器設計計算値との整合確認にそのまま使える。つまり「教科書の演習問題」と「実務の設計検証」を一つのモデルで兼ねられる
なるほど、単なる教科書問題じゃなくて、実務のCode Verificationにそのまま使えるってことですね。NAFEMSのベンチマーク集にも入っているんですか？
もちろん。NAFEMSのLE10やLE11が厚肉円筒系の問題だし、ASME V&V 10-2006で推奨されるCode Verification手順でも「閉じた解析解を持つ問題」の代表例として引用される。新しいソルバーを導入したとき、あるいはバージョンアップ後の回帰テストとして、最初に走らせるべきモデルの一つだよ。
支配方程式 — ラメの公式
では具体的な数式を教えてください。内径 $a$、外径 $b$ の円筒に内圧 $p$ がかかっている場合ですよね？
ラメの公式の導出は、軸対称条件下の力の釣り合い方程式と適合条件から出発する。外圧がゼロ（$p_o = 0$）のとき、任意の半径位置 $r$（$a \le r \le b$）における応力は次のようになる：
半径方向応力（常に圧縮 $\le 0$）：
周方向（フープ）応力（常に引張 $\ge 0$）：
2つの式がすごく似てますね。$b^2/r^2$ の符号が違うだけ……ということは、 $\sigma_r + \sigma_\theta$ はrに依存しない定数になりますか？
鋭い！ その通りだ。
これは軸対称弾性の重要な性質で、力の釣り合いの簡易チェックに使える。FEA結果で各節点の $\sigma_r + \sigma_\theta$ を算出して定数にならなければ、何かがおかしい。
内面 $r = a$ と外面 $r = b$ での値も確認しておきたいです。
肉厚比 $k = b/a$ を使うとスッキリ書ける：
| 位置 | $\sigma_r$ | $\sigma_\theta$ |
|---|---|---|
| 内面 $r = a$ | $-p$ | $\displaystyle p\,\frac{k^2 + 1}{k^2 - 1}$ |
| 外面 $r = b$ | $0$ | $\displaystyle \frac{2p}{k^2 - 1}$ |
内面のフープ応力 $\sigma_\theta(a)$ が常に最大であり、これが圧力容器設計の支配因子になる。例えば $k = 2$（$b/a = 2$）のとき $\sigma_\theta(a) = 5p/3 \approx 1.667p$ だ。
変位場と軸方向応力
FEMとの比較では応力だけでなく変位も見たいです。半径方向変位 $u_r$ の解析解もありますか？
もちろん。平面ひずみ条件（$\varepsilon_z = 0$、長い円筒の仮定）での半径方向変位は：
平面ひずみだと軸方向にも応力が出ますよね？ $\sigma_z$ はどうなりますか？
平面ひずみでは $\varepsilon_z = 0$ だから、Hookeの法則より：
$\sigma_z$ は半径位置 $r$ によらない定数だ。これもFEA結果の検証チェックに使える。もし $\sigma_z$ が要素ごとにばらついていたら、拘束条件の設定を見直す必要がある。
平面応力 vs 平面ひずみ — どちらを使うべきか
- 平面ひずみ（$\varepsilon_z = 0$）: 長い円筒（$L/D \gg 1$）に適用。パイプライン、原子炉圧力容器、砲身など。ほとんどの実務ケースはこちら
- 平面応力（$\sigma_z = 0$）: 薄いリング状部材に適用。実務では稀だが、教科書の演習問題として出題されることがある
- 一般化平面ひずみ（$\varepsilon_z = \text{const.} \ne 0$）: 端面に軸力が作用するケース。ASME設計ではclosed-end条件（端面に内圧が作用）を考慮する場合に該当
von Mises相当応力
実務ではvon Mises応力で評価することが多いと思いますが、ラメ解からvon Mises応力も解析的に出せますか？
$\sigma_r$, $\sigma_\theta$, $\sigma_z$ の3つが主応力そのもの（せん断成分はゼロ）だから、von Mises相当応力は直接計算できる：
平面ひずみの場合、$\sigma_r - \sigma_\theta$ が $r$ に依存する唯一の変数だから、$\sigma_\text{eq}$ も内面 $r = a$ で最大値を取る。例えば後で示す標準ベンチマーク（$a = 50$ mm, $b = 100$ mm, $p = 100$ MPa, $\nu = 0.3$）だと、内面の $\sigma_\text{eq}(a) \approx 236.6$ MPa になる。これがFEAの結果と合っているかを確認するわけだ。
ベンチマーク問題設定
具体的な数値で検証してみたいです。標準的な問題設定を教えてもらえますか？
NAFEMSや各ソルバーのマニュアルでよく使われる設定を基に、以下の標準ベンチマークを定義しよう：
| パラメータ | 記号 | 値 |
|---|---|---|
| 内径 | $a$ | 50 mm |
| 外径 | $b$ | 100 mm（肉厚比 $k = 2$） |
| 内圧 | $p$ | 100 MPa |
| ヤング率 | $E$ | 210 GPa |
| ポアソン比 | $\nu$ | 0.3 |
| 条件 | — | 平面ひずみ、線形弾性、等方性 |
この条件でのラメ解析解（参照値）を示す：
| 物理量 | 内面 $r = a$ | 中間 $r = 75$ mm | 外面 $r = b$ |
|---|---|---|---|
| $\sigma_r$ [MPa] | $-100.00$ | $-18.52$ | $0.00$ |
| $\sigma_\theta$ [MPa] | $+166.67$ | $+85.19$ | $+66.67$ |
| $\sigma_z$ [MPa] | $+20.00$ | $+20.00$ | $+20.00$ |
| $u_r$ [mm] | $0.05317$ | $0.04346$ | $0.03810$ |
| $\sigma_\text{eq}$ [MPa] | $236.56$ | $91.21$ | $57.74$ |
理論解 vs FEA — 検証データ比較表
実際にFEAで解くと、メッシュや要素次数でどれくらい結果が変わるんですか？
上の標準ベンチマーク問題をいろいろなメッシュで解いた結果を比較表にまとめた。評価対象は内面フープ応力 $\sigma_\theta(a)$（理論値 = 166.67 MPa）だ：
| 要素タイプ | 半径方向分割 | DOF | $\sigma_\theta(a)$ [MPa] | $\sigma_r(a)$ [MPa] | 誤差 [%] | 判定 |
|---|---|---|---|---|---|---|
| CAX4（4節点軸対称） | 4 | 50 | 162.3 | $-96.5$ | 2.64 | FAIR |
| CAX4（4節点軸対称） | 8 | 100 | 165.5 | $-99.1$ | 0.72 | PASS |
| CAX4（4節点軸対称） | 16 | 200 | 166.4 | $-99.8$ | 0.18 | PASS |
| CAX8（8節点軸対称） | 4 | 100 | 166.5 | $-99.9$ | 0.12 | PASS |
| CAX8（8節点軸対称） | 8 | 200 | 166.7 | $-100.0$ | 0.00 | PASS |
| C3D20（20節点六面体3D） | 8×16×4 | 18,000 | 166.6 | $-99.9$ | 0.06 | PASS |
判定基準: 相対誤差 < 1%: ■ PASS、1〜5%: ■ FAIR、> 5%: ■ FAIL
2次要素（CAX8）の収束がすごいですね……！ たった4分割で0.12%、8分割でほぼ完全一致。1次要素だと4分割で2.64%もずれるのに。
これがラメ問題の教育的価値でもある。「なぜ2次要素が推奨されるのか」を数値で実感できる。応力が $1/r^2$ に比例する分布だから、1次要素の線形補間では捕えきれない曲率がある。2次要素はこの曲率を形状関数で直接表現できるから、少ない要素でも高精度が出るんだ。
数値解法と実装
軸対称モデルの定式化
FEMで厚肉円筒を解くとき、軸対称要素を使うのが定番だと思いますが、3D要素でも検証すべきですか？
Code Verificationとしてはまず軸対称2Dで始める。モデル構築が簡単で、自由度が少なく、結果の解釈も明確だからだ。その後、3Dモデル（1/4対称の90度セクター等）で同じ結果が再現されることを確認する——これが「コードの3Dソルバーが正しく軸対称問題を処理できるか」の検証になる。
軸対称要素では、変位場は $u_r(r, z)$ と $u_z(r, z)$ の2成分で、ひずみベクトルは4成分：
ここで $\varepsilon_\theta = u_r/r$ が軸対称特有の項であり、半径 $r$ が分母に来るため内面近傍で数値的に敏感になる。
要素剛性マトリクスは体積積分に $2\pi r$ が入る：
要素選択と積分スキーム
積分の方法——完全積分と低減積分——でも結果が変わりますか？
ラメ問題では影響は小さいが、知っておくべきポイントがある：
- 完全積分（2×2 for QUAD4, 3×3 for QUAD8）: 剛性過大評価の傾向だが、ラメ問題は薄肉でないので深刻なロッキングは起きない
- 低減積分（1×1 for QUAD4, 2×2 for QUAD8）: 計算効率は上がるが、1次要素ではアワーグラスモードのリスク。ラメ問題では1要素あたりの変位モードが単純なので通常問題なし
- 非適合モード要素（Abaqus CAX4I 等）: 1次要素の精度不足を補正する追加変位モードを持つ。ラメ問題で1次要素を使うなら推奨
| 要素 | ソルバー名称 | 節点数 | 積分点 | ラメ問題への推奨度 |
|---|---|---|---|---|
| 4節点軸対称（完全積分） | CAX4 / PLANE182 / CQUAD4 | 4 | 2×2 | 基本検証向き |
| 4節点軸対称（非適合モード） | CAX4I / PLANE182 w/ K-opt | 4 | 2×2 | 1次要素なら推奨 |
| 8節点軸対称（完全積分） | CAX8 / PLANE183 / CQUAD8 | 8 | 3×3 | 最推奨 |
| 8節点軸対称（低減積分） | CAX8R / PLANE183 R | 8 | 2×2 | 推奨 |
メッシュ収束性とGCI
先の比較表で収束の傾向はわかりましたが、「定量的にどのメッシュで十分か」をどう判断するんですか？
GCI（Grid Convergence Index）を使う。Roache（1998）が提案した方法で、ASME V&V 20-2009でも推奨されている。3水準のメッシュ（粗・中・細）の結果から、Richardson外挿で離散化誤差の95%信頼区間を推定する。
具体的には、メッシュサイズ $h_1 > h_2 > h_3$ に対する解 $f_1, f_2, f_3$ から：
ここで $r = h_1/h_2 = h_2/h_3$（等比のメッシュ細分化）。観測された収束次数 $p_\text{obs}$ が理論値（1次要素: 2、2次要素: 3〜4）に近ければ、メッシュは漸近収束域にあると言える。
ラメ問題で2次要素を使った場合、半径方向4→8→16分割の3水準で$p_\text{obs} \approx 4$が得られ、GCI$_{12}$（8→16分割の信頼区間）が0.01%未満になるのが典型的だ。つまり8分割の時点でもう十分ということを数量的に示せる。
実践ガイド
解析フローと境界条件
実際にFEAソフトでラメ問題を組むとき、最初から最後まで何をすればいいですか？
ステップバイステップで説明しよう：
- 形状作成: $r$-$z$ 平面上に矩形領域（$a \le r \le b$, $0 \le z \le L$）を作成。$L$ は任意（平面ひずみなので結果に影響しない）
- 材料定義: 線形弾性（$E = 210$ GPa, $\nu = 0.3$）
- メッシュ: 半径方向にバイアスメッシュ（内面を密に）。CAX8推奨。最低4分割
- 境界条件:
- $z = 0$ 面: $u_z = 0$（対称条件）
- $z = L$ 面: $u_z = 0$（平面ひずみの拘束）
- $r = a$ 面: 内圧 $p = 100$ MPa（面圧荷重）
- $r = b$ 面: 自由（外圧 = 0）
- 求解: 静的線形解析。直接法（Sparse solver）で十分
- 後処理: $\sigma_r$, $\sigma_\theta$, $\sigma_z$, $u_r$ を半径方向にパスプロットして理論解と重ねる
境界条件の「$z = L$面で $u_z = 0$」を忘れると、平面ひずみにならないってことですか？
そのとおり。両端の$z$拘束を忘れると一般化平面ひずみ状態になり、$\sigma_z$ が変わってくる。結果として $u_r$ も理論値からずれる。FEM検証で一番多い失敗の一つが、この平面ひずみ拘束の設定漏れだ。
商用ソルバー別の設定ポイント
Abaqus、Ansys、Nastranでそれぞれ注意すべき点はありますか？
主要3ソルバーの設定を整理しよう：
| 項目 | Abaqus | Ansys Mechanical | MSC Nastran |
|---|---|---|---|
| 要素タイプ | CAX8R（推奨）/ CAX4I | PLANE183（KEYOPT(3)=1） | CQUAD8 + PSHELL (AXI) |
| 内圧荷重 | *DSLOAD, P1=100 | SFE, PRES=100 | PLOAD4 |
| 平面ひずみ拘束 | *BOUNDARY, ZSYMM | D, UZ, 0 on both ends | SPC1, UZ=0 |
| 応力出力 | S11=$\sigma_r$, S22=$\sigma_\theta$, S33=$\sigma_z$ | SX, SY, SZ（座標系注意） | 円筒座標系でのVonMises |
| 注意事項 | 節点座標出力で$r$取得 | 円筒座標系（LOCAL=11）定義必須 | CORD2C で円筒座標系定義 |
Ansysで応力成分を取り出すとき、座標系を間違えると全然違う値になりそうですね……
よくある失敗だ。Ansysはデフォルトで直交座標系（XYZ）で応力を出力するから、円筒座標系に変換しないと $\sigma_r$ や $\sigma_\theta$ にならない。Abaqusの軸対称要素はS11, S22が最初から $r, \theta$ 成分だから混乱しにくいが、3Dモデルでは同じ問題が起きる。
ASME Section VIIIとの整合確認
ラメの公式ってASMEの圧力容器規格にも出てきますよね。設計計算値とFEAを比べるときの注意点はありますか？
ASME Section VIII Division 2では、厚肉円筒の許容内圧を次式で規定している：
ここで $S$ は許容応力、$k = b/a$。これはラメの公式から $\sigma_\theta(a) = S$ と置いて逆算した式と一致する。
実務で注意すべきは：
- Division 1 vs Division 2: Division 1は簡易式（薄肉近似ベース）、Division 2はラメの厳密解ベース。Division 2のほうがFEAと直接比較しやすい
- 腐食代（CA）: ASMEでは腐食代を差し引いた有効肉厚で計算する。FEA検証では腐食代ゼロ（公称寸法）で比較する
- Closed-end effect: 実際の圧力容器は端板（鏡板）があるため軸方向に $\sigma_z = pa^2/(b^2 - a^2)$ が作用する。平面ひずみの $\sigma_z$ とは異なるので条件を合わせること
ソフトウェア比較
ソルバー別の結果比較
同じ問題を複数のソルバーで解いた場合、結果に差は出ますか？
ラメ問題のような線形弾性の軸対称問題では、同等の要素・メッシュを使えばソルバー間の差はほぼゼロになる。これが逆に有用で、もしソルバー間で差が出たら「設定ミス」を疑える。
| ソルバー | 要素 | 半径方向8分割 | $\sigma_\theta(a)$ [MPa] | 誤差 [%] |
|---|---|---|---|---|
| Abaqus 2024 | CAX8R | 8 | 166.67 | 0.00 |
| Ansys 2024R2 | PLANE183 | 8 | 166.67 | 0.00 |
| MSC Nastran 2023 | CQUAD8 (AXI) | 8 | 166.66 | 0.01 |
| COMSOL 6.2 | 二次三角形 | 8相当 | 166.65 | 0.01 |
| CalculiX 2.21 | CAX8 | 8 | 166.67 | 0.00 |
| Code_Aster 15.8 | QUAD8_AXI | 8 | 166.66 | 0.01 |
オープンソースのCalculiXやCode_Asterでも商用ソルバーと同等の精度が出るんですね。これならソルバー選定の際の信頼性検証にも使えそうです。
まさにその通り。新しいソルバーを評価するとき、「まずラメ問題を解いて6桁一致するか確認する」のは業界の定石だ。逆に言えば、ラメ問題で誤差が出るソルバーは基本的な実装に問題がある可能性が高い。
先端技術
拡張問題 — 熱応力・弾塑性・クリープ
線形弾性のラメ解は理解できました。でも実際の圧力容器って高温で使うこともあるし、降伏もしますよね。そういう拡張問題でもベンチマークとして使えますか？
いいところに気づいた。ラメ問題は以下の拡張ケースでも解析解・半解析解が存在するので、より高度なCode Verificationに使える：
- 熱応力（NAFEMS LE11）: 内外面に温度差を与えた厚肉円筒。$\sigma_\theta$ に熱応力が重畳される。温度の対数分布と応力の$1/r^2$分布が組み合わさり、メッシュ感度が高まる
- 弾完全塑性: 内圧を徐々に上げると、内面から降伏が始まり弾塑性境界が外側に進展する。弾性域と塑性域の境界半径 $c$ の理論式が存在する（Hill解）
- 定常クリープ: Norton則 $\dot{\varepsilon} = A\sigma^n$ に従うクリープ下の応力再分布。定常状態では応力が内面に集中せず均一化に向かう「クリープによる応力緩和」が理論的に示せる
- 自緊処理（Autofrettage）: 高圧力を加えて内面を意図的に塑性変形させた後、除荷すると残留圧縮応力が発生する。疲労寿命向上の効果を定量評価できる
1つのジオメトリで、線形弾性から弾塑性、クリープ、熱連成まで段階的にV&Vを拡張していけるんですね。まさに「最良のベンチマーク」と呼ばれる理由がわかりました。
その通り。特にFEM初心者には、まず線形弾性のラメ解で完璧な一致を確認し、次に弾塑性に拡張するときに何が変わるかを体験することを強く勧める。これが「検証の階段」（Verification Ladder）の考え方だ。
トラブルシューティング
FEM検証で陥りがちなミス
「理論解と合わない！」というとき、よくある原因をまとめて教えてください。
過去に見てきた失敗パターンをランキング形式でまとめよう：
| 順位 | 失敗パターン | 症状 | 対策 |
|---|---|---|---|
| 1 | 平面ひずみ拘束の設定漏れ | $\sigma_z$ がゼロになり、$u_r$ が理論値から5〜15%ずれる | z方向の両端をUZ=0で拘束。ソルバーによっては「Generalized Plane Strain」の設定 |
| 2 | 座標系の取り違え | $\sigma_r$ と $\sigma_\theta$ が入れ替わる、または直交座標成分で比較してしまう | 軸対称要素の応力成分の定義を確認（Abaqus: S11=$r$, S22=$\theta$）。3Dモデルでは円筒座標系を定義して変換 |
| 3 | 内圧の向き・面の指定ミス | 全体の符号が逆、または外面に圧力がかかっている | 荷重の正の向き（法線方向 or 内向き）をマニュアルで確認。変形図で膨張方向を目視チェック |
| 4 | 単位系の不整合 | 応力が桁違いに大きい/小さい | mm-N-MPa系 or m-N-Pa系で統一。$E$ = 210,000 MPa = 2.1×10$^{11}$ Pa |
| 5 | 1次要素の限界を知らない | 「メッシュを細かくしても2%の誤差が消えない」と悩む | CAX4を32分割するよりCAX8を4分割したほうが高精度。まず2次要素で「正解」を確認してから1次要素の限界を調べる |
| 6 | 応力の外挿 vs 積分点値 | 節点外挿値と積分点値で応力が異なる（特に内面） | ラメ問題の検証では積分点値（Element Output）で比較。節点外挿値は平均化で鈍る場合がある |
1位が「平面ひずみ拘束の設定漏れ」って、すごく基本的なところなんですね……。自分もやりそう。
基本的なことほど見落としやすい。対策としては、解析を走らせたらまず $\sigma_z$ の値を確認することを習慣にするといい。理論値は $\sigma_z = 2\nu a^2 p/(b^2 - a^2)$。今回の例だと $\sigma_z = 20.0$ MPa で一定。もし $\sigma_z \approx 0$ なら平面応力になっている、不均一なら拘束に問題がある——とすぐに判断できる。
なるほど。$\sigma_z$ を「カナリアの鳥」として使うわけですね。チェックリストに加えます！
最後にもう一つ。「解析が合わない」と思ったら、以下のデバッグ手順を踏むといい：
- 変形図を見る — 膨張の方向は正しいか？ 不自然な変形モードはないか？
- $\sigma_r + \sigma_\theta$ が一定か確認 — ラメ解の基本性質。ばらつきがあれば設定ミス
- $\sigma_z$ が一定か確認 — 平面ひずみ拘束の検証
- 外面の $\sigma_r = 0$ を確認 — 外圧ゼロの境界条件が正しく反映されているか
- 内面の $\sigma_r = -p$ を確認 — 内圧の符号と大きさが正しいか
この5項目をクリアすれば、ほぼ確実に理論値と一致するはずだ。
なった
詳しく
報告