応力特異点とメッシュ収束 — FEMで応力が発散する原理と実務対策

バグじゃないんだ。弾性力学の理論解そのものが「応力無限大」を予測しているケースがある。90度のリエントラントコーナー（L字の内角270°の角）やき裂先端がその典型だよ。

そう、まさにそこがポイント。FEMは弾性理論の近似解を求めているから、メッシュを細かくするほど「理論上の無限大」に近づいていく。つまり $\sigma_{\max}$ はメッシュサイズ $h \to 0$ で際限なく増加する。通常のメッシュ収束性チェック——3水準のメッシュで応力が一定値に収束するか確認する——が成立しないんだ。

まさにその通り。FEM初心者が最大von Mises応力で許容応力と比較する——これが特異点を含むモデルでの典型的な過設計パターンだ。メッシュを粗くした人は「OK」、細かくした人は「NG」になるから、結果がメッシュ担当者のさじ加減で変わってしまう。これはV&V的に完全にアウトだよ。

Williams（1952）が導いた古典的な固有値解析がある。くさび形領域の頂点からの距離を $r$、角度を $\theta$ とすると、応力場は次のように表される：

ここで $K$ は応力強度に関する係数、$f_{ij}(\theta)$ は角度依存の関数、そして $\lambda$ がWilliams固有値だ。この $\lambda$ はくさびの開き角 $2\alpha$ に依存して決まる。ポイントは、$\lambda < 1$ のとき $r \to 0$ で $\sigma \to \infty$ になるということ。

Williamsの固有値方程式を解くと、内角270°（リエントラントコーナー）では $\lambda \approx 0.545$ だ。つまり応力は $\sigma \propto r^{-0.455}$ のオーダーで発散する。き裂先端（内角360°）ではおなじみの $\lambda = 0.5$ で $\sigma \propto r^{-1/2}$ になる。

その通り。だから実務ではフィレットR（丸み）を付けて角をなくすことで特異性を回避する。ただし、CADモデルにフィレットがあってもメッシュが角を十分に表現できないほど粗いと、数値的に特異性が残ることもあるので注意が必要だ。

正確に言うと、特異点から離れた位置の応力はちゃんと収束する。問題は特異点直近の最大応力だ。標準的なFEMの誤差評価式を見てみよう：

ここで $k$ は要素の多項式次数、$\lambda$ はWilliams固有値だ。滑らかな解なら $\lambda \geq k$ で $h^k$ の収束率が得られる。しかし特異点があると $\lambda < 1 < k$ なので、収束率は $h^\lambda$ に制限される。2次要素を使っても $h^2$ ではなく $h^{0.545}$ 程度にしかならない。

グローバルなエネルギーノルム誤差に対しては改善が遅いだけで意味がないわけではない。ただし、特異点での応力値そのものは要素次数を上げてもメッシュを細かくしても収束しない——これは本質的な限界だ。FEMのバグではなく、物理モデル（線形弾性体の鋭角コーナー）が持つ本質的な性質なんだよ。

その通り。だから「どこの応力で評価するか」が特異点問題の本質なんだ。最大応力のコンター図を見て「赤いところが危ない！」と反射的に判断するのは、特異点があるモデルでは完全に間違いだよ。

代表的なものを挙げるね：

そう。ボルト締結の座面やベアリングの接触面を「1節点に集中荷重」で代替すると、その節点近傍の応力は物理的に意味がない。少なくとも荷重を面に分布させるか、荷重作用点近傍の応力は評価対象から外すのが鉄則だ。

いい質問だね。破壊力学では「応力値」ではなく「応力拡大係数 $K$」や「$J$ 積分」で評価するから、特異点での応力値そのものは使わない。ただし、FEMで $K$ を精度よく求めるための工夫としてquarter-point要素（1/4点要素）がある。

Barsoum（1976）が考案した手法だ。8節点四辺形の2次要素で、き裂先端側の中間節点を辺の1/4の位置にずらす。すると要素内の変位場に $\sqrt{r}$ の項が自然に現れて、き裂先端の $r^{-1/2}$ 特異性を要素レベルで正確に捕捉できる。

AbaqusやAnsysではき裂先端周りにこのquarter-point要素を配置する機能がビルトインされている。これと後述のJ積分やcontour integral法を組み合わせることで、$K_I, K_{II}, K_{III}$ を高精度に求められるんだ。

hp-refinementが有効だ。特異点に向かって幾何学的級数（geometric grading）でメッシュを細分化し、同時に要素の多項式次数を上げる。Babuskaらの研究によると、この組み合わせにより指数関数的な収束が達成できる：

ここで $N$ は自由度数。通常のh-refinementでは代数的（$N^{-\lambda/d}$）にしか収束しないのに対し、hp-refinementは指数的に収束するから圧倒的に効率がいい。ただし実装は複雑で、対応しているソルバーは限られる。StressCheckやp-FEMベースのソルバーが代表例だね。

主に3つの方法がある：

その通り。だから特異点直上の応力値は使わず、やや離れた複数コンターの平均値で $J$ を求める。Abaqusでは *CONTOUR INTEGRAL で5〜10コンターを自動計算してくれる。外側のコンターほどメッシュの影響を受けにくくて安定するんだ。

ASME Boiler & Pressure Vessel Code Section VIII Division 2 が定めている応力線形化（stress linearization）が業界標準だ。これは板厚方向の応力分布を「膜応力 $\sigma_m$」と「曲げ応力 $\sigma_b$」に分解する手法だよ。

ここで $t$ は板厚、$x$ は板厚方向の座標だ。膜応力は断面を一様に引っ張る成分、曲げ応力は表裏で正負が反転する成分。特異点によるピーク応力 $\sigma_p$（線形化で残る非線形分布の残差）は一次評価の対象外とされる。

「一次応力評価」では無視する。ただし疲労評価ではピーク応力が重要になる。ASME Div.2 Part 5.5では疲労解析に full-range stress（膜+曲げ+ピーク）を使う。このとき特異点のピーク応力は物理的に意味がないから、フィレットR付きのサブモデルに切り替えて実際の応力集中を評価するのが正攻法だ。

き裂が存在する（または存在を仮定する）場合は、破壊力学パラメータで評価に切り替える。代表的な方法は3つ：

まさにそう。応力は $r^{-1/2}$ で発散するが、$J$ 積分はパス非依存で有限値に収束する。$K$ 値も $J$ から一意に定まる。特異性のある問題ではポイント値（応力）からグローバル値（エネルギー）への転換が本質的な解法なんだ。

サブモデリングは2つの目的で使う：

ただし注意点がある。サブモデルの切断境界（driven boundary）が特異点から十分離れていること、グローバルモデルの応力がその境界位置で十分収束していることが前提条件だ。サブモデルを特異点に近づけすぎると、境界の変位場自体が不正確になってしまう。

完璧なまとめだ。補足すると（0番目として）「フィレットRを付けて特異性自体をなくす」も重要だ。設計段階でR付けが可能なら、そもそも特異点を避けるのがベストだよ。

主要ソルバーの対応状況を比較するとこうなる：

Ansys WorkbenchのStress Linearization Toolは板厚方向のパスを指定するだけで $P_m, P_b, P_m+P_b, P_m+P_b+Q$ を自動計算してくれるから、圧力容器の設計審査では非常に便利だ。Abaqusでも *PATH でパスを定義してPythonスクリプトで線形化計算はできるが、ワンクリックというわけにはいかない。

オープンソースでは、CalculiXはJ積分計算に対応している。Code_AsterはフランスEDFが開発した汎用FEMで、POST_K1_K2_K3コマンドによる応力拡大係数の抽出やXFEMにも対応していて、破壊力学に関しては商用ソルバー並みの機能を持っているよ。

従来のFEMでは、き裂はメッシュの要素境界に沿って表現する必要があった。つまりき裂が進展するたびにリメッシュが必要だった。XFEMでは既存のメッシュに「エンリッチメント関数」を追加することで、メッシュに依存せずにき裂を表現できる。

$H(\mathbf{x})$ はき裂面を横切るHeaviside関数（不連続性の表現）、$F_\ell$ はき裂先端近傍の $\sqrt{r}$ 特異性を表す関数群だ。き裂先端の特異性をメッシュではなく形状関数レベルで取り込むから、粗いメッシュでも高精度な $K$ 値が得られるのが最大の利点だよ。

ただし万能ではない。3Dの複雑な形状でき裂面が分岐・合流するケースや、大きな塑性変形を伴う延性破壊では、まだ実用上の課題がある。Abaqusでは3D XFEMが使えるが、き裂進展は基本的にパリス則ベースの疲労き裂成長に限られる。任意方向への脆性き裂伝播は2Dでは良いが、3Dでは安定性に課題が残っているのが現状だ。

位相場破壊モデル（Phase-field fracture）は、き裂を鋭い不連続面ではなく、連続的なダメージ変数 $d \in [0, 1]$ で表現する手法だ。$d=0$ が健全、$d=1$ が完全破壊。Griffithのエネルギー釣り合いを変分原理で定式化するもので、Bourdinら（2000）が提案した。

第1項が弾性エネルギー（$g(d)$で劣化）、第2項がき裂表面エネルギー。$G_c$ は臨界エネルギー解放率、$\ell$ は正則化長さパラメータだ。この方法の最大の利点はき裂の発生・進展・分岐・合流を自動的に捕捉できること。メッシュにき裂形状を事前定義する必要がない。

2026年現在、商用ソルバーへの実装はまだ限定的だ。Abaqus 2024からUEL/UMATベースの位相場サポートが強化されてきているが、まだ研究段階の色合いが強い。XFEMの方が実務では先を行っている。ただしアカデミアでは位相場法の論文数がここ5年で爆発的に増えていて、今後5〜10年で商用ソルバーにも標準搭載されていくだろう。

よし、実務で見てきた典型的な失敗を5つ挙げるよ。

そう。だからメッシュ収束性のチェックは「最大応力が収束すること」ではなく「評価点の応力が収束すること」を確認するものだと理解してほしい。コンター図の赤い場所が特異点かどうかを判別する目を養うことが、FEMエンジニアとしての重要なスキルだよ。

よし、いいまとめだ。最後にもう一つだけ。社内レビューで「最大応力がメッシュで変わるんですが…」と言われたとき、「それは応力特異点です。Williamsの固有値は $\lambda = 0.545$ で理論的に収束しません。評価方法を応力線形化またはJ積分に切り替えます」と説明できれば、一目置かれるエンジニアになれるよ。頑張れ！