単純支持梁の等分布荷重 — FEM検証ベンチマーク
理論と物理
概要
単純支持梁のベンチマークって何をチェックするんですか？
FEMの最も基本的な検証問題だよ。等分布荷重 $w$ を受ける単純支持梁には、Euler-Bernoulli梁理論から厳密なたわみ解が得られる。中央の最大たわみは
この理論解と数値解を比較して、要素タイプの精度を確認するのがこのベンチマークの目的だ。特に重要なのは、1次要素（線形要素）を使うとシア・ロッキングが発生して、梁が実際よりも硬く振る舞い、たわみが過小評価されること。この現象を定量的に把握し、低減積分やB-bar法での改善効果を検証するのが実務的なポイントになる。
え、たわみの式がひとつあるだけで、そこまでいろんなことがわかるんですか？
そう。シンプルだからこそ使い勝手がいいんだ。NAFEMSのベンチマーク集にも収録されているし、ASME V&V 10-2006のCode Verification（コードの数学的正確性の確認）の第一歩としても推奨されている。新しいソルバーを導入したとき、まずこの問題で理論解と0.1%以内の一致を確認するのが業界の鉄板手順だ。たとえば自動車メーカーの構造解析部門でも、ソルバーのバージョンアップのたびにこのベンチマークを回帰テストとして走らせている。
支配方程式
じゃあ、理論式の全体像を教えてください。たわみだけじゃなくて応力も出るんですよね？
Euler-Bernoulli梁の支配方程式から始めよう。等分布荷重 $w$（N/m）を受ける長さ $L$ の単純支持梁のたわみ曲線は
ここで $x$ は左端からの距離、$E$ はヤング率、$I$ は断面二次モーメントだ。中央 $x = L/2$ に代入すると最大たわみが得られる。
次に曲げモーメント分布。
中央で最大となり $M_{\max} = wL^2/8$ だ。ここから最大曲げ応力は
$Z = I/c$ は断面係数、$c$ は中立軸から最外縁までの距離だ。矩形断面 $b \times h$ なら $I = bh^3/12$、$c = h/2$、$Z = bh^2/6$ になる。
せん断力の分布もありますよね？
せん断力は
支点で最大 $V_{\max} = wL/2$、中央でゼロだ。Euler-Bernoulli理論では断面が平面を保ちかつ法線方向を維持する（せん断変形を無視する）仮定だから、スパン/せい比 $L/h$ が大きい薄い梁で精度が高い。
Timoshenko梁との比較
$L/h$ が小さいときはどうなるんですか？ Timoshenko理論との差が気になります。
Timoshenko梁理論ではせん断変形を追加で考慮する。中央たわみは
第2項がせん断変形によるたわみ増分だ。$\kappa$ はせん断補正係数（矩形断面で $5/6$）、$G$ はせん断弾性係数 $G = E/\{2(1+\nu)\}$ だ。
| $L/h$ | 曲げたわみ比率 | せん断たわみ比率 | Euler-Bernoulli誤差 |
|---|---|---|---|
| 50 | 99.7% | 0.3% | 0.3% |
| 20 | 98.1% | 1.9% | 1.9% |
| 10 | 92.9% | 7.1% | 7.1% |
| 5 | 78.1% | 21.9% | 21.9% |
$L/h \geq 20$ ならEuler-Bernoulliで十分だ。$L/h < 10$ の深い梁ではTimoshenko理論か3Dソリッド要素が必要になる。ベンチマークでは通常 $L/h = 50$ 程度を使って、Euler-Bernoulli理論が厳密に成立する条件で検証する。
ベンチマーク検証データ
具体的な数値で確認したいんですが、どんなパラメータを使うのが標準ですか？
標準的なベンチマーク設定はこうだ。
| パラメータ | 記号 | 値 |
|---|---|---|
| スパン長 | $L$ | 1000 mm |
| 断面幅 | $b$ | 100 mm |
| 断面高さ | $h$ | 20 mm |
| 等分布荷重 | $w$ | 1.0 N/mm |
| ヤング率 | $E$ | 200,000 MPa |
| ポアソン比 | $\nu$ | 0.3 |
$L/h = 50$ だからEuler-Bernoulli理論が十分成立する。断面二次モーメント $I = bh^3/12 = 100 \times 20^3/12 = 66{,}667$ mm$^4$。
ここで $Z = bh^2/6 = 100 \times 20^2/6 = 6{,}667$ mm$^3$ だ。
要素タイプごとの結果はどうなりますか？
上記パラメータでの比較結果だ。
| 要素タイプ | メッシュ | DOF | $\delta_{\max}$ [mm] | 変位誤差 [%] | $\sigma_{\max}$ [MPa] | 応力誤差 [%] |
|---|---|---|---|---|---|---|
| BEAM2（Euler-Bernoulli） | 20要素 | 126 | 0.9766 | 0.00 | 18.75 | 0.00 |
| BEAM3（Timoshenko） | 20要素 | 126 | 0.9795 | 0.30 | 18.75 | 0.00 |
| QUAD8（二次シェル） | 20×2 | 2,520 | 0.9764 | 0.02 | 18.72 | 0.16 |
| HEX8 完全積分 | 50×10×4 | 19,800 | 0.8112 | 16.93 | 17.92 | 4.43 |
| HEX8 低減積分 | 50×10×4 | 19,800 | 0.9738 | 0.29 | 18.70 | 0.27 |
| HEX20（二次ソリッド） | 20×4×2 | 15,120 | 0.9762 | 0.04 | 18.73 | 0.11 |
| TET4（線形四面体） | 自動 | ~30,000 | 0.6834 | 30.02 | 15.21 | 18.88 |
| TET10（二次四面体） | 自動 | ~25,000 | 0.9741 | 0.26 | 18.68 | 0.37 |
HEX8の完全積分で約17%、TET4に至っては30%もの変位誤差が出ている。これがシア・ロッキングの威力だ。一方、低減積分のHEX8やHEX20なら誤差0.3%以内。この差をはっきり見せるのが単純支持梁ベンチマークの意義なんだ。
数値解法と実装
シア・ロッキングの物理
さっきの表でHEX8やTET4がひどい結果だったのは「シア・ロッキング」のせいですよね。これって具体的に何が起きてるんですか？
本質を理解するには、1次要素の形状関数の限界を考えるといい。
1次六面体要素（HEX8）の変位場は各方向に線形だ。純曲げでは断面が回転するから、上面と下面で逆方向の水平変位が生じる。理論的にはこのときせん断ひずみ $\gamma_{xy} = 0$ のはずだが、線形の変位場ではこの条件を満たせない。
つまり本来存在しないはずの寄生的なせん断ひずみ（parasitic shear）が生じてしまう。このひずみに対応するせん断応力が余計な剛性として作用し、梁を実際よりも硬く、つまりたわみを小さく見積もる。これがシア・ロッキングだ。
なるほど…要素の形状関数が曲げの変形パターンを表現できないから、「嘘のせん断」が出てくるわけですね。
そうだ。メッシュを細かくすれば個々の要素のアスペクト比が改善されて緩和はするが、収束が非常に遅い。厚さ方向に2要素程度では全然足りなくて、8〜16要素くらい必要になる。計算コストが膨大になるから、実務的にはロッキング対策技法を使うほうがはるかに効率的だ。
ロッキング対策技法
対策の方法にはどんなものがありますか？
主なアプローチは3つある。
1. 低減積分（Reduced Integration）
完全積分（HEX8なら2×2×2=8点）の代わりに、低減積分（1×1×1=1点）を使う。積分点を減らすことで、寄生的なせん断ひずみのエネルギー寄与を排除できる。ただし代償としてアワーグラスモード（hourglass mode）という零エネルギーの変形モードが生じる。砂時計状にぐにゃっと変形してもエネルギーがゼロと評価されてしまう。
2. B-bar法（選択的低減積分）
体積変化（ダイレーション）成分だけを低減積分で計算し、偏差成分は完全積分を維持する。体積ロッキングとせん断ロッキングの両方に効果がある。Abaqusの C3D8R + hourglass control や Nastranの CHEXA reduced がこの発想に基づいている。
3. EAS法（Enhanced Assumed Strain）
通常の変位場に加えて、要素内部で追加のひずみモードを仮定する。内部自由度として追加のひずみパラメータを導入し、静的凝縮で外部には見えないようにする。AnsysのSOLID185（keyopt(2)=2）がこれに該当する。
| 手法 | 利点 | 注意点 |
|---|---|---|
| 低減積分 | 実装が簡単、計算コスト低い | アワーグラス制御が必須 |
| B-bar法 | 体積・せん断ロッキング両対応 | 要素技術の理解が必要 |
| EAS法 | 高精度、アワーグラスなし | 計算コストやや増、非線形で不安定な場合あり |
| 2次要素 | 根本解決、曲げモードを正確表現 | DOF数が約3倍、メッシュ生成に手間 |
アワーグラスモードの制御って、具体的にどうやるんですか？
アワーグラスモードを抑制するために人工的な剛性（hourglass stiffness）を加える。アワーグラスエネルギーが全体の内部エネルギーの5〜10%を超えたら危険サインだ。Abaqusなら *SECTION CONTROLS, HOURGLASS=ENHANCED、LS-DYNAなら IHQ=6（Belytschko-Bindeman）が実績ある設定だ。このベンチマーク問題でアワーグラスエネルギー比率を確認して、1%以下であれば合格と判断できる。
メッシュ収束性の理論
メッシュを細かくしていくと、どのくらいの速さで理論解に近づくんですか？
FEMのアプリオリ誤差評価（Ceaの定理）によれば、$p$ 次要素でエネルギーノルム誤差は $O(h^p)$、変位の $L^2$ ノルム誤差は $O(h^{p+1})$ で減少する。
つまり線形要素（$p=1$）で要素サイズ $h$ を半分にすると変位誤差は約 $1/4$ に、二次要素（$p=2$）なら約 $1/8$ になる。この理論的収束率が実際のメッシュ収束テストで再現されるかを確認するのがCode Verificationの核心だ。
具体的には3つ以上の異なるメッシュ密度で解析し、両対数プロットで傾きが理論値 $p+1$ に近いかを確認する。傾きが理論値から大きく外れる場合は、バグ、境界条件の誤り、あるいは特異点の影響を疑うことになる。
実践ガイド
解析の手順
実際にソルバーでこのベンチマークを走らせるとき、どんな手順でやればいいですか？
基本的な流れはこうだ。
Step 1: ジオメトリ作成
$L \times b \times h = 1000 \times 100 \times 20$ mm の直方体。対称性を利用して半分モデル（$L/2$）にしてもよい。
Step 2: メッシュ生成
最初は粗いメッシュ（スパン方向10分割、厚さ方向2分割）で始めて、順次細分化する。メッシュ収束テストの推奨はスパン方向 10→20→40→80 の4段階。
Step 3: 材料定義
線形弾性。$E = 200{,}000$ MPa、$\nu = 0.3$。
Step 4: 境界条件設定
左端：$u_y = 0$（鉛直方向拘束）+ $u_z = 0$（面外変位拘束）
右端：$u_y = 0$ + $u_z = 0$
片方の端で $u_x = 0$（軸方向剛体移動を拘束）
Step 5: 荷重適用
上面に面圧 $p = w/b = 1.0/100 = 0.01$ MPa。または線荷重 $w = 1.0$ N/mm。
Step 6: 求解・後処理
中央断面の最大たわみと最大応力を理論解と比較。
境界条件の注意点
境界条件で失敗しやすいポイントはありますか？
よくある落とし穴が3つある。
1. 過拘束
「単純支持」なのに両端で $u_x = 0$ を設定してしまうケース。軸方向に熱膨張のような伸びが拘束されて、寄生的な軸力が発生する。片方だけ $u_x = 0$ にすること。
2. 点拘束によるキンクの発生
ソリッド要素で支点を1点（単一節点）だけで拘束すると、そこに応力集中が生じてしまう。支点側の辺または面に沿って節点群を拘束する、あるいはMPC（多点拘束）で支持線を定義するのが正しい。
3. 面外方向の拘束忘れ
3D モデルでは $z$ 方向の剛体回転を防ぐために $u_z = 0$ を面外方向に設定する必要がある。忘れると特異剛性マトリクスになりソルバーがエラーを吐く。
品質チェックリスト
解析が終わったあと、何をチェックすればOKと言えますか？
ベンチマーク合格のためのチェックリストだ。
| 項目 | 合格基準 |
|---|---|
| 中央たわみ $\delta_{\max}$ | 理論解との誤差 ±1% 以内 |
| 最大曲げ応力 $\sigma_{\max}$ | 理論解との誤差 ±2% 以内 |
| 支点反力の合計 | $wL$ と一致（力の釣り合い確認） |
| メッシュ収束性 | 2段階以上の細分化で単調収束 |
| 収束率の傾き | 両対数プロットで理論値 $p+1$ の ±0.2 以内 |
| アワーグラスエネルギー比率（低減積分時） | 全内部エネルギーの1%以下 |
| 対称性の確認 | 中央に対して左右対称のたわみ分布 |
ソフトウェア比較
ソルバー別の要素選択
ソルバーごとにどの要素を使えばいいか教えてください。
主要ソルバーでの推奨要素と設定をまとめよう。
| ソルバー | 推奨要素 | ロッキング対策 | 入力カード/設定 |
|---|---|---|---|
| Abaqus | C3D20R（二次六面体低減積分） | 内蔵 | *ELEMENT, TYPE=C3D20R |
| Abaqus | C3D8R（一次六面体低減積分） | アワーグラス制御 | *SECTION CONTROLS, HOURGLASS=ENHANCED |
| Ansys Mechanical | SOLID186（二次六面体） | 内蔵 | デフォルト設定でOK |
| Ansys Mechanical | SOLID185（一次六面体） | EAS法 | KEYOPT(2)=2 (Enhanced Strain) |
| MSC Nastran | CHEXA（20節点） | 内蔵 | PSOLID, FCTN=SMECH |
| CalculiX | C3D20R | 内蔵 | Abaqus互換入力 |
| Code_Aster | HEXA20 | 内蔵 | MODELISATION='3D' |
クロスバリデーション結果
ソルバー間で結果に差は出ますか？
二次六面体（HEX20相当）で統一した場合のクロスバリデーション結果だ。
| ソルバー | 要素 | $\delta_{\max}$ [mm] | 誤差 [%] | $\sigma_{\max}$ [MPa] |
|---|---|---|---|---|
| Abaqus 2024 | C3D20R | 0.9762 | 0.04 | 18.73 |
| Ansys 2024R2 | SOLID186 | 0.9763 | 0.03 | 18.74 |
| MSC Nastran 2024 | CHEXA-20 | 0.9761 | 0.05 | 18.72 |
| CalculiX 2.22 | C3D20R | 0.9762 | 0.04 | 18.73 |
| 理論解 | — | 0.9766 | — | 18.75 |
全ソルバーで誤差0.05%以内に収まっている。二次要素を正しく使えば、商用・オープンソースを問わずほぼ同じ結果が得られる。逆にソルバー間で結果が乖離する場合は、要素定式化か境界条件の設定ミスを疑うべきだ。
先端技術
非線形への拡張
この問題を非線形に拡張するとどうなりますか？
2つの方向で拡張できる。
幾何学的非線形（大変形）
たわみが梁のせい $h$ と同程度以上になると、曲率変化が無視できなくなりEuler-Bernoulli理論の線形仮定が破綻する。大変形解析では、梁の中立軸の伸びと回転を考慮したvon Karmanのひずみ-変位関係を使う。
たわみが $\delta/h > 0.5$ を超えると非線形効果が5%以上になるから、その場合はこちらの定式化に切り替える必要がある。
材料非線形（弾塑性）
応力が降伏応力を超える場合、von Mises降伏条件＋等方硬化則などの弾塑性構成則が必要になる。このベンチマークの線形弾性版で検証が通ったら、次のステップとして弾塑性解析に進み、降伏荷重 $w_y$ や崩壊荷重 $w_c$ を評価する問題に拡張する。
Richardson外挿
厳密解がわからない問題では、どうやって「真の解」を推定するんですか？
Richardson外挿を使う。3つ以上のメッシュ密度の結果から、メッシュサイズ $h \to 0$ の極限値を推定する方法だ。
ここで $f_1, f_2$ は細かい順のメッシュの解、$r$ はメッシュ細分化比（通常2）、$p$ は観測された収束次数だ。この問題では厳密解がわかっているから、Richardson外挿の推定値が厳密解と一致するかどうかを検証できる。実務で厳密解のない複雑な問題に移行したとき、この技法が生きてくるわけだ。
トラブルシューティング
理論値と合わない場合
解析結果が理論値と合わないとき、何から疑えばいいですか？
切り分けの優先順位はこうだ。
1. たわみが小さすぎる（剛性が高い）場合
- シア・ロッキング → 要素タイプを確認。完全積分のHEX8/TET4なら低減積分に切替、または二次要素に変更
- 過拘束 → 両端で $u_x = 0$ にしていないか確認
- 材料定数のミス → $E$ の単位を確認（GPa と MPa の混同が多い）
2. たわみが大きすぎる場合
- 荷重の単位ミス → $w$ の単位が N/mm なのか N/m なのか確認
- 断面寸法の入力ミス → $h$ を mm で入れるべきところに m で入れると $I$ が $10^{12}$ 倍ずれる
- アワーグラスモード暴走 → 低減積分でアワーグラス制御が無効になっていないか確認
3. 応力値だけがずれる場合
- 応力の評価位置 → 積分点の値と節点外挿値で差が出る。理論値との比較では積分点が最も信頼性が高い
- ポアソン比の影響 → 3Dソリッドでは $\nu$ によるポアソン効果で梁理論と微小な差が生じる（正常）
応力の評価位置
応力は中央断面のどこを見ればいいんですか？ 上面？ 下面？ 積分点？
理論上、$\sigma_{\max}$ は中央断面 ($x = L/2$) の最外縁（上面または下面）に生じる。FEMでは以下の点に注意が必要だ。
積分点（Gauss点）の応力が最も直接的な計算値で、外挿や平均化の誤差を含まない。ただし積分点は要素内部にあるため、表面の値を直接見ていない。
節点外挿値は積分点から外挿して節点に出した値で、最外縁の応力を見るにはこちらを使う。隣接要素間で値の不連続が見られる場合は、平均化（nodal averaging）の方法も確認する必要がある。
実務的には、メッシュが十分に細かければ積分点値と節点外挿値の差は1%以内に収まる。差が大きい場合はメッシュ不足のサインだ。このベンチマーク問題でその差を定量的に確認しておくと、実際の設計解析でも応力値の信頼度を判断する力がつくよ。
すごく整理されました。シンプルな問題ほどきちんとやると、複雑な問題の基礎体力になるってことですね。
その通り。単純支持梁ベンチマークは「FEMの体温計」のようなものだ。ここで要素の特性、ロッキングの影響、収束性の感覚を掴んでおけば、実務のどんな問題に遭遇しても「この結果は信頼できるのか？」という判断軸ができる。V&Vの第一歩として、ぜひ自分の手でソルバーを動かしてみてくれ。
なった
詳しく
報告