1912年4月15日未明、「沈まない船」と称されたタイタニック号は氷山に衝突し、船体が真っ二つに折れて沈没しました。第二次世界大戦中には、アメリカの溶接構造貨物船リバティ船が航行中に突然船体が破断する事故が相次ぎ、約2,700隻のうち400隻以上に深刻なき裂が発生しました。これらの事故には共通点があります。設計上の応力は材料の降伏応力をはるかに下回っていたにもかかわらず、破壊が起きたということです。
従来の強度設計は「部材に生じる応力が材料の許容応力を超えなければ安全」という考え方に基づいていました。しかし、材料内部に微小なき裂（欠陥）が存在する場合、き裂先端には局所的に極めて高い応力が集中し、全体の平均応力がいくら低くても破壊が起こり得ます。この「き裂がある材料はどの条件で壊れるか」を定量的に扱う学問が破壊力学（fracture mechanics）です。
破壊力学を理解すると、以下のような幅広い分野で活用できます。
- 航空宇宙工学: 航空機の機体やエンジン部品に存在する微小欠陥の許容サイズを決定し、検査間隔を合理的に設定する損傷許容設計（damage tolerance design）の基盤です
- 原子力・プラント工学: 圧力容器や配管のき裂進展を予測し、構造健全性を評価するために破壊力学パラメータが使われます
- 土木・建築: 橋梁の鋼部材における疲労き裂の進展速度を予測し、補修・交換時期を判断します
- 材料開発: 新しい合金やセラミックスの開発において、破壊靱性値 $K_{IC}$ は材料選定の最も重要な指標の一つです
本記事の内容
- 従来の強度設計の限界と破壊力学の必要性
- グリフィスのエネルギー基準（表面エネルギーとひずみエネルギー解放率）
- 応力拡大係数 $K$ の定義と3つの破壊モード（モード I / II / III）
- き裂先端の応力場（Westergaard 解）と Python による可視化
- 破壊靱性 $K_{IC}$ と破壊条件
- エネルギー解放率 $G$ と応力拡大係数 $K$ の関係
- 疲労き裂進展とパリス則
- Python による数値シミュレーション（臨界き裂長さ、形状係数、き裂進展予測）
前提知識
この記事を読む前に、以下の記事に目を通しておくと理解がスムーズです。
- 応力とひずみの定義を基礎から徹底解説 — 応力・ひずみの定義と応力テンソルの基礎
- フックの法則と弾性係数 — 弾性体の構成則とヤング率・ポアソン比の関係
- 応力集中 — 切欠きや穴による局所的な応力の増大現象
- 平面応力と平面ひずみ — 2次元弾性問題の定式化
従来の強度設計の限界 — なぜ破壊力学が必要か
許容応力設計の考え方
構造設計の最も基本的なアプローチは、部材に生じる最大応力 $\sigma_{\max}$ が材料の降伏応力 $\sigma_Y$ や引張強さ $\sigma_B$ を安全率 $n$ で割った許容応力以下であることを確認する方法です。
$$ \begin{equation} \sigma_{\max} \leq \frac{\sigma_Y}{n} \end{equation} $$
このアプローチは「材料は均質で欠陥がない」という前提に立っています。しかし、現実の材料には必ず何らかの欠陥が存在します。鋳造品には気孔（ポロシティ）があり、溶接部にはブローホールやスラグ巻き込みがあり、金属材料の内部には製造過程で生じた微小なき裂が存在します。
応力とひずみ で学んだように、応力は「力を面積で割った量」です。しかし、き裂のような鋭い欠陥の先端では、面積が限りなくゼロに近づくため、理論上は応力が無限大に発散します。つまり、き裂が存在する場合、「最大応力を許容値以下に抑える」という発想そのものが破綻するのです。
イングリス（Inglis）の応力集中の解
この問題に初めて定量的に取り組んだのは、イングリスです。1913年、イングリスは無限板中の楕円孔周りの応力解を導き、楕円の長半径 $a$、先端曲率半径 $\rho$ の孔がある場合、先端の最大応力が次のようになることを示しました。
$$ \begin{equation} \sigma_{\max} = \sigma_{\infty}\left(1 + 2\sqrt{\frac{a}{\rho}}\right) \end{equation} $$
ここで $\sigma_{\infty}$ は遠方の一様応力です。応力集中 の記事で詳しく扱ったように、応力集中係数 $K_t = \sigma_{\max}/\sigma_{\infty}$ は欠陥の形状に大きく依存します。
問題は、き裂のように先端が極めて鋭い欠陥（$\rho \to 0$）の場合です。$\rho \to 0$ とすると $\sigma_{\max} \to \infty$ となり、どんなに小さな外力でも先端応力は無限大になってしまいます。これでは破壊条件を応力で記述することができません。
この行き詰まりを打破するために、まったく異なる視点 — エネルギーの視点 — から破壊問題に取り組んだのがグリフィスです。
グリフィスのエネルギー基準
発想の転換 — 応力からエネルギーへ
アラン・アーノルド・グリフィス（A.A. Griffith）は1921年、画期的な論文を発表しました。彼の発想は次の通りです。「き裂先端の応力が無限大に発散するなら、応力ではなくエネルギーで破壊を議論すればよい」。
グリフィスは、き裂の進展を「エネルギーの収支」として捉えました。き裂が微小量 $da$ だけ伸びるとき、以下の2つのエネルギー変化が起こります。
- ひずみエネルギーの解放: き裂が伸びると、き裂面近傍の応力が緩和され、物体に蓄えられていた弾性ひずみエネルギーの一部が解放されます。これはき裂進展を促進する効果です
- 表面エネルギーの増加: き裂が伸びると新しい破面が2面生じ、その表面エネルギーの分だけ系のエネルギーが増加します。これはき裂進展に抵抗する効果です
き裂進展を一種の「エネルギー取引」に例えると、わかりやすいでしょう。弾性ひずみエネルギーの解放が「収入」、新表面の形成が「支出」です。収入が支出を上回るとき、き裂進展は「採算が合う」ため進行し、収入が支出を下回るときは、き裂は進展しません。
エネルギー保存則 で学んだ通り、孤立系のエネルギーは保存されます。グリフィスの議論は、外力が仕事をしない固定変位条件のもとで、全ポテンシャルエネルギー（弾性ひずみエネルギー + 表面エネルギー）が最小となる条件としてき裂の平衡を定義するものです。
数学的定式化
無限板中に長さ $2a$ の中央き裂があり、遠方で一様引張応力 $\sigma$ を受ける問題を考えます（平面応力条件）。イングリスの解から、き裂が存在することで解放されるひずみエネルギーは次のように求められます。
$$ \begin{equation} U_e = -\frac{\pi a^2 \sigma^2}{E} \end{equation} $$
ここで $E$ はヤング率、板の厚さは単位厚さとしています。負号は、き裂の存在によってエネルギーが解放される（減少する）ことを表します。
一方、長さ $2a$ のき裂によって生じる新しい表面は、き裂の上面と下面の2面で、その面積は $2 \times 2a = 4a$ です。単位面積あたりの表面エネルギーを $\gamma_s$ とすると、表面エネルギーは次のようになります。
$$ \begin{equation} U_s = 4a\gamma_s \end{equation} $$
系の全ポテンシャルエネルギーの変化量 $\Delta U$ は、これらの和です。
$$ \begin{equation} \Delta U = U_e + U_s = -\frac{\pi a^2 \sigma^2}{E} + 4a\gamma_s \end{equation} $$
き裂の平衡条件は、全ポテンシャルエネルギーが極値をとる条件、すなわち $\Delta U$ を $a$ で微分してゼロとおくことで得られます。
$$ \begin{equation} \frac{d(\Delta U)}{da} = -\frac{2\pi a \sigma^2}{E} + 4\gamma_s = 0 \end{equation} $$
この式を $\sigma$ について解くと、グリフィスの破壊応力が得られます。
$$ \begin{equation} \sigma_f = \sqrt{\frac{2E\gamma_s}{\pi a}} \end{equation} $$
グリフィスの臨界き裂長さ
上の式を $a$ について解くと、与えられた応力 $\sigma$ のもとで不安定破壊が始まる臨界き裂長さ $a_c$ が得られます。
$$ \begin{equation} a_c = \frac{2E\gamma_s}{\pi \sigma^2} \end{equation} $$
この式は非常に重要な物理的意味を持っています。
- き裂長さ $a < a_c$ のとき: 表面エネルギーの増加がひずみエネルギーの解放を上回るため、き裂は進展しません。き裂は安定状態にあります
- き裂長さ $a = a_c$ のとき: エネルギーの収支がちょうど釣り合い、き裂は臨界状態にあります
- き裂長さ $a > a_c$ のとき: ひずみエネルギーの解放が表面エネルギーの増加を上回り、き裂は不安定に進展します
この結果は定性的には素晴らしいものでしたが、グリフィスの理論がそのまま金属材料に適用できるかというと、問題がありました。ガラスのような脆性材料ではグリフィスの予測は実験とよく一致しましたが、金属材料では予測される破壊応力が実際の値よりもはるかに低くなってしまいました。
その理由は、金属材料ではき裂先端近傍で塑性変形が生じ、そのエネルギー消費 $\gamma_p$ が表面エネルギー $\gamma_s$ をはるかに上回るためです。オロワン（Orowan）とアーウィン（Irwin）は、グリフィスの式の $\gamma_s$ を有効表面エネルギー $\gamma_{\text{eff}} = \gamma_s + \gamma_p$ に置き換えることでこの問題を解決しました。
$$ \begin{equation} \sigma_f = \sqrt{\frac{2E(\gamma_s + \gamma_p)}{\pi a}} \approx \sqrt{\frac{2E\gamma_p}{\pi a}} \end{equation} $$
金属材料では $\gamma_p \gg \gamma_s$ であるため、実質的に塑性エネルギーのみで破壊応力が決まります。
エネルギー基準は破壊力学の基盤ですが、実用上はき裂先端の応力場を直接特徴づけるパラメータがあると便利です。その役割を果たすのが、アーウィンが導入した応力拡大係数です。
ひずみエネルギー解放率 G
G の定義
グリフィスのエネルギー基準をより一般的に定式化するため、アーウィンはひずみエネルギー解放率（strain energy release rate）$G$ を導入しました。$G$ はき裂が単位面積だけ進展するときに解放されるひずみエネルギーとして定義されます。
$$ \begin{equation} G = -\frac{dU_e}{dA} = -\frac{1}{B}\frac{dU_e}{da} \end{equation} $$
ここで $A$ はき裂面積、$B$ は板の厚さ、$a$ はき裂長さです。$G$ の単位は $\text{J}/\text{m}^2$ または $\text{N}/\text{m}$ です。「率」（rate）という名前が付いていますが、時間に対する変化率ではなく、き裂面積に対する変化率であることに注意してください。
中央き裂の場合、先ほど求めたひずみエネルギーの式から $G$ を計算すると次のようになります。
平面応力条件では、
$$ \begin{equation} G = \frac{\pi a \sigma^2}{E} \end{equation} $$
平面応力と平面ひずみ で学んだ平面ひずみ条件では、ヤング率の代わりに $E/(1-\nu^2)$ を使い、
$$ \begin{equation} G = \frac{\pi a \sigma^2 (1-\nu^2)}{E} \end{equation} $$
となります。ここで $\nu$ はポアソン比です。
臨界エネルギー解放率 $G_c$
グリフィスのエネルギー基準を $G$ で書き直すと、き裂進展の条件は次のように表せます。
$$ \begin{equation} G \geq G_c = 2\gamma_{\text{eff}} \end{equation} $$
$G_c$ は臨界エネルギー解放率（critical energy release rate）と呼ばれ、材料固有の定数です。$G < G_c$ であればき裂は安定であり、$G = G_c$ に達した瞬間に不安定破壊が始まります。
$G$ はき裂を「駆動する力」、$G_c$ はき裂の進展に「抵抗する力」と解釈できます。この考え方は、き裂の安定性を評価する上で非常に直感的です。
ここまでエネルギーの観点から破壊を議論してきましたが、き裂先端の応力場に目を向けると、もう一つの強力なパラメータが現れます。それが応力拡大係数 $K$ です。
応力拡大係数 K の定義
3つの破壊モード
き裂先端の変形は、き裂面の開口方向によって3つの基本モードに分類されます。
- モード I（開口モード）: き裂面に垂直な引張力によって、き裂が開口するモード。工学的に最も重要で、破壊事故の大半はモード I で起こります
- モード II（面内せん断モード）: き裂面内でき裂面に平行、かつき裂前縁に垂直な方向にずれるモード。面内せん断によって駆動されます
- モード III（面外せん断モード）: き裂面内でき裂前縁に平行な方向にずれるモード。反平面せん断（ねじり）によって駆動されます
実際の構造物ではこれら3つのモードが混合して現れることがありますが、純粋なモード I が最も危険であり、破壊力学の議論は主にモード I に焦点を当てます。
応力拡大係数の数学的定義
応力拡大係数 $K$ は、き裂先端の応力場の「強さ」を表すパラメータです。モード I の応力拡大係数 $K_I$ は次のように定義されます。
$$ \begin{equation} K_I = \lim_{r \to 0} \sqrt{2\pi r} \, \sigma_{yy}(r, \theta=0) \end{equation} $$
ここで $r$ はき裂先端からの距離、$\theta$ はき裂面からの角度、$\sigma_{yy}$ はき裂面に垂直な方向の応力成分です。この定義は、「き裂先端近傍の応力が $1/\sqrt{r}$ で発散する」という特異性の強さを測っているものです。
無限板中の長さ $2a$ の中央き裂に遠方一様引張応力 $\sigma$ が作用する場合、$K_I$ は次のように求められます。
$$ \begin{equation} K_I = \sigma\sqrt{\pi a} \end{equation} $$
この式は破壊力学で最も基本的かつ重要な式の一つです。$K_I$ の単位は $\text{Pa}\sqrt{\text{m}}$ ですが、工学の実用上は $\text{MPa}\sqrt{\text{m}}$ がよく使われます。
$K_I$ の式をよく見ると、外力の大きさ（$\sigma$）とき裂の大きさ（$a$）が一つのパラメータに統合されていることがわかります。これが応力拡大係数の最大の利点です。き裂のある構造物の危険度を、たった一つの数値で評価できるのです。
同様に、モード II およびモード III の応力拡大係数もそれぞれ定義されます。
$$ \begin{equation} K_{II} = \tau\sqrt{\pi a}, \quad K_{III} = \tau_{\parallel}\sqrt{\pi a} \end{equation} $$
ここで $\tau$ は面内せん断応力、$\tau_{\parallel}$ は面外せん断応力です。
形状係数（補正係数）
上の $K_I = \sigma\sqrt{\pi a}$ は無限板中の中央き裂という理想的な条件に対する解です。有限幅の板、片側き裂、曲げ荷重など現実の構造では、き裂形状や荷重条件の効果を補正する形状係数（geometry factor, correction factor）$Y$ を導入して次のように書きます。
$$ \begin{equation} K_I = Y \sigma \sqrt{\pi a} \end{equation} $$
形状係数 $Y$ はき裂と構造の幾何学的条件によって決まる無次元量です。代表的な値を示します。
| 形状 | 形状係数 $Y$ |
|---|---|
| 無限板中の中央き裂 | $Y = 1$ |
| 無限板中の片側き裂 | $Y = 1.12$ |
| 有限幅 $W$ の板の中央き裂 | $Y = \sqrt{\sec\left(\frac{\pi a}{W}\right)}$ |
| 有限幅 $W$ の板の片側き裂 | $Y = 1.12 – 0.231\left(\frac{a}{W}\right) + 10.55\left(\frac{a}{W}\right)^2 – 21.72\left(\frac{a}{W}\right)^3 + 30.39\left(\frac{a}{W}\right)^4$ |
応力集中 で学んだ応力集中係数 $K_t$ と応力拡大係数 $K_I$ を混同しないよう注意してください。$K_t$ は「応力の比」（無次元量）であるのに対し、$K_I$ はき裂先端の応力場の強さを表す物理量であり、単位を持ちます。
応力拡大係数の定義を押さえたところで、次にき裂先端の応力場がどのような数学的構造を持つのかを詳しく見ていきましょう。
き裂先端の応力場 — Westergaard 解
応力場の一般形
き裂先端近傍の応力場は、線形弾性力学の枠組みで厳密に求めることができます。ウエスターガード（Westergaard）やウィリアムズ（Williams）による解析の結果、モード I のき裂先端近傍の応力場は次のように表されます。
$$ \begin{equation} \sigma_{xx} = \frac{K_I}{\sqrt{2\pi r}}\cos\frac{\theta}{2}\left(1 – \sin\frac{\theta}{2}\sin\frac{3\theta}{2}\right) \end{equation} $$
$$ \begin{equation} \sigma_{yy} = \frac{K_I}{\sqrt{2\pi r}}\cos\frac{\theta}{2}\left(1 + \sin\frac{\theta}{2}\sin\frac{3\theta}{2}\right) \end{equation} $$
$$ \begin{equation} \tau_{xy} = \frac{K_I}{\sqrt{2\pi r}}\sin\frac{\theta}{2}\cos\frac{\theta}{2}\cos\frac{3\theta}{2} \end{equation} $$
ここで $r$ と $\theta$ はき裂先端を原点とする極座標で、$\theta = 0$ がき裂の延長方向（き裂の前方）に対応します。
これらの式から、き裂先端の応力場について極めて重要な特徴が読み取れます。
第一に、全ての応力成分が $1/\sqrt{r}$ に比例して発散します。き裂先端に近づくほど応力は無限大に増大し、$r = 0$ では特異点となります。これがまさに、応力基準による破壊評価が破綻する根本的な理由です。
第二に、応力場の「形」（角度依存性）はき裂の大きさや外力によらず普遍的です。$\cos(\theta/2)$, $\sin(\theta/2)$ などの三角関数で決まるこの角度分布は、全てのモード I き裂に共通です。
第三に、応力場の「強さ」は $K_I$ のみで決まります。すなわち、$K_I$ が同じであれば、き裂の形状や荷重条件が異なっていても、き裂先端の応力場はまったく同一です。これをK 支配の原理（$K$-dominance）と呼びます。
$K$ 支配の原理こそが破壊力学を実用的なものにしている核心です。小さな試験片で測定した $K_I$ の臨界値を、大型の実構造物にそのまま適用できるのは、き裂先端の応力場が $K_I$ のみで特徴づけられるからです。
Westergaard 解の導出の概要
Westergaard 解は複素関数を使った解法で得られます。ここでは導出の骨子を示します。
2次元弾性問題において、応力成分は応力関数（エアリーの応力関数）$\Phi$ を用いて次のように表されます。
$$ \begin{equation} \sigma_{xx} = \frac{\partial^2 \Phi}{\partial y^2}, \quad \sigma_{yy} = \frac{\partial^2 \Phi}{\partial x^2}, \quad \tau_{xy} = -\frac{\partial^2 \Phi}{\partial x \partial y} \end{equation} $$
偏微分の基礎 で学んだ偏微分を使って、この応力関数が重調和方程式 $\nabla^4 \Phi = 0$ を満たすことが弾性体の平衡条件から要求されます。
Westergaard は、モード I 問題に対して複素変数 $z = x + iy$ の解析関数 $Z_I(z)$ を用いて応力関数を構成しました。き裂先端近傍で $Z_I(z)$ を展開すると、$z = a$ の近傍で $Z_I \propto 1/\sqrt{z – a}$ となる特異性が現れ、これが $1/\sqrt{r}$ 応力特異性の起源です。
変数変換 $z – a = re^{i\theta}$ を行い、実部と虚部を分離すると、先ほど示した応力場の式が得られます。
き裂先端の変位場
参考として、き裂先端近傍の変位場も示しておきます。モード I の変位成分は次の通りです。
$$ \begin{equation} u_x = \frac{K_I}{2\mu}\sqrt{\frac{r}{2\pi}}\cos\frac{\theta}{2}\left(\kappa – 1 + 2\sin^2\frac{\theta}{2}\right) \end{equation} $$
$$ \begin{equation} u_y = \frac{K_I}{2\mu}\sqrt{\frac{r}{2\pi}}\sin\frac{\theta}{2}\left(\kappa + 1 – 2\cos^2\frac{\theta}{2}\right) \end{equation} $$
ここで $\mu$ はせん断弾性係数（フックの法則参照）、$\kappa$ は定数で、平面応力では $\kappa = (3-\nu)/(1+\nu)$、平面ひずみでは $\kappa = 3-4\nu$ です。
変位は $\sqrt{r}$ に比例するため、$r \to 0$ で有限値（ゼロ）に収束します。応力が発散するのに変位が有限であるのは、一見矛盾しているように見えますが、ひずみ $\varepsilon \propto du/dr \propto 1/\sqrt{r}$ も発散するため、応力とひずみ の関係（フックの法則）$\sigma = E\varepsilon$ と整合しています。
それでは、き裂先端の応力場を Python で可視化して、数式から読み取った特徴を目で確認しましょう。
Python: き裂先端の応力場の可視化
き裂先端の応力場の空間分布を Python で描画し、$1/\sqrt{r}$ 特異性と角度依存性を視覚的に確認します。
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定
K_I = 1.0 # 応力拡大係数（規格化して 1.0）
# グリッド生成（き裂先端を原点とする）
x = np.linspace(-0.5, 2.0, 500)
y = np.linspace(-1.5, 1.5, 500)
X, Y = np.meshgrid(x, y)
# 極座標への変換
r = np.sqrt(X**2 + Y**2)
theta = np.arctan2(Y, X)
# ゼロ除算を避けるため微小値で置換
r_safe = np.where(r < 1e-10, 1e-10, r)
# モード I のき裂先端応力場
prefactor = K_I / np.sqrt(2 * np.pi * r_safe)
sigma_xx = prefactor * np.cos(theta / 2) * (1 - np.sin(theta / 2) * np.sin(3 * theta / 2))
sigma_yy = prefactor * np.cos(theta / 2) * (1 + np.sin(theta / 2) * np.sin(3 * theta / 2))
tau_xy = prefactor * np.sin(theta / 2) * np.cos(theta / 2) * np.cos(3 * theta / 2)
# き裂先端近傍で応力値が大きくなりすぎるためクリップ
vmax = 3.0
sigma_xx = np.clip(sigma_xx, -vmax, vmax)
sigma_yy = np.clip(sigma_yy, -vmax, vmax)
tau_xy = np.clip(tau_xy, -vmax, vmax)
# 可視化（3成分を並べて表示）
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
titles = [r'$\sigma_{xx}$', r'$\sigma_{yy}$', r'$\tau_{xy}$']
data = [sigma_xx, sigma_yy, tau_xy]
for ax, title, d in zip(axes, titles, data):
im = ax.contourf(X, Y, d, levels=50, cmap='RdBu_r', vmin=-vmax, vmax=vmax)
# き裂を黒線で表示（x < 0 の部分）
ax.plot([-0.5, 0], [0, 0], 'k-', linewidth=3)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title(title, fontsize=16)
ax.set_aspect('equal')
plt.colorbar(im, ax=ax, shrink=0.8)
plt.suptitle('Mode I Crack Tip Stress Field', fontsize=16, y=1.02)
plt.tight_layout()
plt.show()
このグラフから、き裂先端の応力場について以下の重要な特徴が視覚的に確認できます。
-
$\sigma_{yy}$（き裂面に垂直な応力）: き裂の前方（$\theta = 0$）で最も高い応力が集中しており、き裂先端に近づくほど急激に増大しています。これはき裂を開口させようとする駆動力の分布を直接反映しています。き裂面上（$\theta = \pm\pi$）では $\sigma_{yy} = 0$ となり、自由表面の境界条件を満たしていることも確認できます
-
$\sigma_{xx}$（き裂に平行な応力）: き裂前方でも高い引張応力が生じていますが、$\sigma_{yy}$ と比べると角度依存性のパターンが異なります。$\theta = 0$ では $\sigma_{xx} = \sigma_{yy}$ となり、き裂の真正面では等二軸応力状態であることがわかります
-
$\tau_{xy}$（せん断応力）: $\theta = 0$ および $\theta = \pm\pi$ で $\tau_{xy} = 0$ となる反対称のパターンを示しています。これはモード I の対称性（き裂面に関する鏡映対称性）の直接的な帰結です
次に、グリフィスの理論で導いた臨界き裂長さを Python で計算し、材料ごとの違いを確認しましょう。
Python: グリフィスの臨界き裂長さの計算
さまざまな材料について、作用応力に対する臨界き裂長さの関係をプロットします。グリフィスの式に塑性エネルギーを含む修正版を用いることで、脆性材料と延性材料の違いを比較します。
import numpy as np
import matplotlib.pyplot as plt
# 材料パラメータ（代表値）
materials = {
'Glass': {'E': 70e9, 'gamma_eff': 1.0, 'color': '#e74c3c'},
'Al 7075': {'E': 71e9, 'gamma_eff': 20e3, 'color': '#3498db'},
'Steel (mild)': {'E': 200e9, 'gamma_eff': 100e3, 'color': '#2ecc71'},
'Ti-6Al-4V': {'E': 114e9, 'gamma_eff': 50e3, 'color': '#9b59b6'},
}
# 応力範囲
sigma = np.linspace(10e6, 500e6, 500) # 10 MPa ~ 500 MPa
fig, ax = plt.subplots(figsize=(10, 7))
for name, props in materials.items():
E = props['E']
gamma_eff = props['gamma_eff']
# グリフィスの臨界き裂長さ（修正版）
a_c = 2 * E * gamma_eff / (np.pi * sigma**2)
# メートルからミリメートルに変換
a_c_mm = a_c * 1e3
ax.plot(sigma / 1e6, a_c_mm, label=name, color=props['color'], linewidth=2)
ax.set_xlabel('Applied Stress [MPa]', fontsize=13)
ax.set_ylabel('Critical Crack Length $a_c$ [mm]', fontsize=13)
ax.set_title('Griffith Critical Crack Length vs Applied Stress', fontsize=14)
ax.set_yscale('log')
ax.set_xlim([10, 500])
ax.set_ylim([1e-3, 1e6])
ax.legend(fontsize=12)
ax.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.show()
このグラフから、材料特性が臨界き裂長さに与える影響を明確に読み取ることができます。
-
ガラス（Glass）は極めて小さな欠陥でも破壊する: ガラスの有効表面エネルギーが純粋な表面エネルギー（約 1 J/m$^2$）のみであるため、100 MPa の応力下で臨界き裂長さはわずか数マイクロメートル程度です。これがガラスが脆い理由であり、表面の微小な傷が致命的な欠陥となります
-
金属材料は塑性変形によって大幅にき裂に強い: アルミニウム合金や鋼は、き裂先端での塑性変形エネルギーが表面エネルギーを数桁上回るため、同じ応力でもはるかに大きなき裂を許容できます。軟鋼（mild steel）は特に $\gamma_{\text{eff}}$ が大きく、靱性に優れた材料であることがわかります
-
応力の2乗に反比例する関係: 全ての材料で $a_c \propto 1/\sigma^2$ の関係が成り立っており、対数グラフ上で傾き $-2$ の直線となっています。これは応力が2倍になると臨界き裂長さが4分の1になることを意味し、高応力環境では極めて小さな欠陥管理が重要になることを示しています
ここまでエネルギー基準と応力拡大係数を個別に議論してきましたが、実はこの2つのパラメータには深い数学的な関係があります。次にその関係を導きましょう。
エネルギー解放率 G と応力拡大係数 K の関係
G と K の等価性
エネルギー基準のパラメータ $G$ と応力場のパラメータ $K$ は、見た目は異なるアプローチから生まれましたが、以下の関係で結ばれています。
平面応力条件:
$$ \begin{equation} G = \frac{K_I^2}{E} \end{equation} $$
平面ひずみ条件:
$$ \begin{equation} G = \frac{K_I^2(1-\nu^2)}{E} \end{equation} $$
この関係は、アーウィンがき裂の仮想的な閉合エネルギーを計算することで導出しました。具体的には、き裂が微小量 $\delta a$ だけ進展するときに解放されるエネルギーを、き裂面を閉じるのに必要な仕事として計算するき裂閉口積分（crack closure integral）によって得られます。
導出の概要を示しましょう。き裂先端の前方（$r = \delta a – x$, $\theta = 0$）で $\sigma_{yy}$ が作用しており、き裂が $\delta a$ だけ進展した後のき裂面の変位が $u_y$（$\theta = \pi$）であるとします。き裂閉口に必要な仕事は次のように書けます。
$$ \begin{equation} G = \lim_{\delta a \to 0}\frac{1}{\delta a}\int_0^{\delta a}\frac{1}{2}\sigma_{yy}(r, 0) \cdot 2u_y(\delta a – r, \pi)\,dr \end{equation} $$
$\sigma_{yy}$ と $u_y$ にそれぞれき裂先端の漸近式を代入すると、
$$ \begin{equation} \sigma_{yy}(r, 0) = \frac{K_I}{\sqrt{2\pi r}}, \quad u_y(\delta a – r, \pi) = \frac{K_I}{2\mu}\sqrt{\frac{\delta a – r}{2\pi}}(\kappa + 1) \end{equation} $$
これらを積分に代入します。
$$ \begin{equation} G = \lim_{\delta a \to 0}\frac{1}{\delta a}\int_0^{\delta a}\frac{K_I^2}{4\pi\mu}(\kappa+1)\sqrt{\frac{\delta a – r}{r}}\,dr \end{equation} $$
$r = \delta a \cdot t$（$0 \leq t \leq 1$）と変数変換すると、
$$ \begin{equation} G = \frac{K_I^2(\kappa+1)}{4\pi\mu}\int_0^{1}\sqrt{\frac{1 – t}{t}}\,dt \end{equation} $$
この定積分は $\pi/2$ と評価されるため、
$$ \begin{equation} G = \frac{K_I^2(\kappa+1)}{8\mu} \end{equation} $$
平面応力条件では $\kappa = (3-\nu)/(1+\nu)$、$\mu = E/2(1+\nu)$ を代入すると、
$$ \begin{equation} G = \frac{K_I^2}{E} \end{equation} $$
平面ひずみ条件では $\kappa = 3 – 4\nu$ を代入すると、
$$ \begin{equation} G = \frac{K_I^2(1-\nu^2)}{E} \end{equation} $$
が得られます。
この $G$-$K$ の関係は、「エネルギーで考えても応力場で考えても、同じ結論に到達する」ことを保証しています。実用的には $K$ を用いた評価の方が便利なことが多いため、現代の破壊力学では $K$ が主役です。ただし、弾塑性破壊力学では $G$ の一般化である $J$ 積分が重要な役割を果たします。
$G$ と $K$ の等価性を確認したところで、次に破壊力学の最も実用的な概念 — 材料の破壊靱性 — について解説します。
破壊靱性 $K_{IC}$ と破壊条件
破壊靱性とは
破壊靱性（fracture toughness）$K_{IC}$ は、材料がき裂の不安定進展に対してどれだけ耐えられるかを表す材料固有の定数です。添え字の “I” はモード I、”C” は臨界（critical）を意味します。
直感的には、$K_{IC}$ は「材料のき裂に対する抵抗力」を表しています。$K_{IC}$ が大きい材料ほど、大きなき裂やき裂に高い応力が作用しても壊れにくい — つまり靱い（tough）材料です。
破壊条件は極めて簡潔に表現されます。
$$ \begin{equation} K_I \geq K_{IC} \quad \Longrightarrow \quad \text{不安定破壊（急速破壊）} \end{equation} $$
逆に $K_I < K_{IC}$ であれば、き裂は安定であり、急速破壊は起こりません。
代表的な材料の破壊靱性値
代表的な工学材料の $K_{IC}$ の値を示します。
| 材料 | $K_{IC}$ [MPa$\sqrt{\text{m}}$] | 特徴 |
|---|---|---|
| ソーダガラス | 0.7 – 0.8 | 極めて脆い |
| アルミナ (Al$_2$O$_3$) | 3 – 5 | セラミックス |
| Al 2024-T3 | 26 – 37 | 航空機用アルミ合金 |
| Al 7075-T6 | 24 – 29 | 高強度アルミ合金 |
| Ti-6Al-4V | 55 – 75 | チタン合金 |
| 軟鋼 (低炭素鋼) | 140 – 200 | 延性に富む |
| マルエージング鋼 | 80 – 120 | 超高強度鋼 |
この表から、一般に強度と靱性はトレードオフの関係にあることが読み取れます。例えば、Al 7075-T6 は Al 2024-T3 よりも引張強さが高い一方、$K_{IC}$ はやや低くなっています。また、マルエージング鋼は軟鋼よりも引張強さがはるかに高いですが、靱性は劣ります。
このトレードオフは材料選定における重要な設計判断のポイントです。航空機の構造設計では、必要な強度を満たしつつ、損傷許容設計の観点から十分な靱性を確保する材料が求められます。
平面ひずみ条件の重要性
$K_{IC}$ は「材料固有の定数」と述べましたが、厳密には平面ひずみ条件のもとで測定された値のみが真の材料定数として扱えます。薄い板では平面応力に近い状態となり、き裂先端の塑性域が大きくなるため、見かけ上の破壊靱性値が増大します。
ASTM E399 規格では、有効な $K_{IC}$ を得るための最小試験片寸法を次のように規定しています。
$$ \begin{equation} B, a, (W – a) \geq 2.5\left(\frac{K_{IC}}{\sigma_Y}\right)^2 \end{equation} $$
ここで $B$ は板厚、$a$ はき裂長さ、$W$ は試験片幅、$\sigma_Y$ は降伏応力です。この条件により、き裂先端の塑性域が試験片寸法に比べて十分小さく、平面ひずみ条件が成立することが保証されます。
許容き裂長さの設計
破壊力学を設計に適用する際の基本的な考え方は、「構造部材に存在する（あるいは存在し得る）き裂の $K_I$ が $K_{IC}$ を超えないこと」を確認することです。
作用応力 $\sigma$ が既知のとき、許容できる最大き裂長さ $a_{\max}$ は次のように求められます。
$$ \begin{equation} K_I = Y\sigma\sqrt{\pi a_{\max}} = K_{IC} \end{equation} $$
$$ \begin{equation} a_{\max} = \frac{1}{\pi}\left(\frac{K_{IC}}{Y\sigma}\right)^2 \end{equation} $$
非破壊検査で検出できる最小き裂長さが $a_{\max}$ 以下であれば、検査をパスした部材は安全と判断できます。逆に、検出限界が $a_{\max}$ を超える場合は、設計応力の低減、靱性の高い材料への変更、または検査技術の改善が必要になります。
ここまでの議論は全て「一回の荷重で壊れるかどうか」の問題でした。しかし実際の構造物は、繰返し荷重を受けることが圧倒的に多く、一回では壊れない応力レベルでも、荷重の繰返しによってき裂が徐々に進展し、最終的に破壊に至ります。次に、この疲労き裂進展の問題を扱います。
疲労き裂進展とパリス則
疲労破壊のメカニズム
疲労破壊 の記事で解説したように、疲労破壊は構造物の破壊事故の最大の原因です。機械要素や構造部材が設計荷重以下の繰返し荷重を受けて、巨視的な塑性変形を伴わずに突然破断する現象です。
疲労破壊のプロセスは大きく3段階に分けられます。
- き裂発生（crack initiation）: 応力集中部や材料欠陥を起点として微小なき裂が発生します
- き裂進展（crack propagation）: 繰返し荷重によってき裂が徐々に成長します
- 最終破壊（final fracture）: き裂が臨界長さに達し、不安定破壊が起こります
破壊力学は主に第2段階と第3段階を定量的に扱います。特に、「繰返し荷重1サイクルあたりにき裂がどれだけ進展するか」を予測できれば、き裂が臨界長さに達するまでの残存寿命を計算できます。
パリス則
1963年、パリス（Paris）とアードガン（Erdogan）は、疲労き裂進展速度 $da/dN$（$N$ は荷重サイクル数）が応力拡大係数範囲 $\Delta K$ のべき乗で表されることを示しました。
$$ \begin{equation} \frac{da}{dN} = C(\Delta K)^m \end{equation} $$
ここで $\Delta K = K_{\max} – K_{\min}$ は1サイクルの間の応力拡大係数の変動幅、$C$ と $m$ は材料定数です。$\Delta K$ は $\Delta\sigma$（応力範囲）と $a$ を用いて次のように書けます。
$$ \begin{equation} \Delta K = Y \Delta\sigma \sqrt{\pi a} \end{equation} $$
パリス則の物理的な意味を押さえておきましょう。$\Delta K$ が大きいほど（すなわち、応力変動が大きいほど、あるいはき裂が長いほど）、1サイクルあたりのき裂進展量が大きくなります。指数 $m$ は金属材料では典型的に 2 ～ 4 の範囲にあり、$m = 3$ が代表的な値です。$m = 3$ は $\Delta K$ が2倍になるとき裂進展速度が $2^3 = 8$ 倍になることを意味しており、き裂進展が加速的に進むことを示しています。
疲労き裂進展の3つの領域
$\log(da/dN)$ を $\log(\Delta K)$ に対してプロットすると、疲労き裂進展曲線は典型的に3つの領域に分かれます。
領域 I（しきい値近傍）: $\Delta K$ が下限しきい値 $\Delta K_{\text{th}}$ 以下のとき、き裂はほとんど進展しません。$\Delta K_{\text{th}}$ は通常 2 – 10 MPa$\sqrt{\text{m}}$ の範囲にあります。$\Delta K < \Delta K_{\text{th}}$ の条件を満たせば、無限寿命設計が可能です。
領域 II（パリス領域）: 両対数グラフ上で直線的な関係が成り立つ領域で、パリス則が適用されます。実際の寿命予測の大部分はこの領域で行われます。
領域 III（急速破壊近傍）: $K_{\max}$ が $K_{IC}$ に近づくと、き裂進展速度が急激に増大し、最終破壊に至ります。
残存寿命の計算
パリス則を用いて、初期き裂長さ $a_0$ から臨界き裂長さ $a_c$ に至るまでの荷重サイクル数（残存寿命）$N_f$ を計算できます。
パリス則を変数分離して積分します。
$$ \begin{equation} dN = \frac{da}{C(\Delta K)^m} = \frac{da}{C(Y\Delta\sigma\sqrt{\pi a})^m} \end{equation} $$
両辺を積分すると、
$$ \begin{equation} N_f = \int_{a_0}^{a_c}\frac{da}{C(Y\Delta\sigma)^m(\pi a)^{m/2}} \end{equation} $$
形状係数 $Y$ が $a$ に依存しない定数として扱える場合（あるいは定数で近似できる場合）、$m \neq 2$ のとき閉じた形で積分できます。
$$ \begin{equation} N_f = \frac{1}{C(Y\Delta\sigma)^m \pi^{m/2}}\cdot\frac{1}{1 – m/2}\left[a_c^{1-m/2} – a_0^{1-m/2}\right] \end{equation} $$
$m > 2$ の場合（ほとんどの金属材料がこれに該当）、$a_c^{1-m/2} \ll a_0^{1-m/2}$ となるため、残存寿命は主に初期き裂長さ $a_0$ で決まります。これは重要な工学的知見です — き裂が小さいうちの進展速度が遅いため、寿命の大部分は初期の段階で消費されるのです。
それでは、Python を使って応力拡大係数の形状係数や疲労き裂進展シミュレーションを実装し、理論を数値的に検証しましょう。
Python: 応力拡大係数の形状係数
有限幅板のき裂に対する形状係数 $Y$ がき裂長さ/板幅比 $a/W$ にどのように依存するかを可視化します。
import numpy as np
import matplotlib.pyplot as plt
# a/W の範囲
a_over_W = np.linspace(0.01, 0.7, 200)
# 形状係数の計算
# (1) 中央き裂（有限幅板）: Y = sqrt(sec(pi*a/W))
Y_center = np.sqrt(1.0 / np.cos(np.pi * a_over_W / 2))
# (2) 片側き裂（有限幅板）: Tada-Paris-Irwin の多項式近似
Y_edge = (1.12 - 0.231 * a_over_W
+ 10.55 * a_over_W**2
- 21.72 * a_over_W**3
+ 30.39 * a_over_W**4)
# (3) 三点曲げ（SENB）: ASTM E399 の式
S_W = 4.0 # スパン/幅比（典型値）
Y_senb = (3 * S_W * np.sqrt(a_over_W)
/ (2 * (1 + 2 * a_over_W) * (1 - a_over_W)**1.5)
* (1.99 - a_over_W * (1 - a_over_W)
* (2.15 - 3.93 * a_over_W + 2.7 * a_over_W**2)))
fig, ax = plt.subplots(figsize=(10, 7))
ax.plot(a_over_W, Y_center, 'r-', linewidth=2.5, label='Center crack (finite width)')
ax.plot(a_over_W, Y_edge, 'b-', linewidth=2.5, label='Single edge crack')
ax.plot(a_over_W, Y_senb, 'g--', linewidth=2.5, label='SENB (S/W=4)')
ax.axhline(y=1.0, color='gray', linestyle=':', alpha=0.5, label='Infinite plate (Y=1)')
ax.axhline(y=1.12, color='gray', linestyle='--', alpha=0.5, label='Infinite plate edge (Y=1.12)')
ax.set_xlabel('$a/W$ (Crack length / Width)', fontsize=13)
ax.set_ylabel('Geometry Factor $Y$', fontsize=13)
ax.set_title('Stress Intensity Factor Geometry Corrections', fontsize=14)
ax.set_xlim([0, 0.7])
ax.set_ylim([0, 8])
ax.legend(fontsize=11, loc='upper left')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このグラフから、形状係数の振る舞いについて以下のことが読み取れます。
-
$a/W$ が小さい領域では無限板の解に収束する: 中央き裂では $Y \to 1$、片側き裂では $Y \to 1.12$ に近づいています。き裂が板幅に比べて十分小さい場合、有限幅の効果は無視でき、無限板の解が良い近似となります
-
$a/W$ が大きくなると形状係数は急増する: 特に $a/W > 0.5$ では $Y$ が急激に増大しています。これは、き裂が板幅の半分を超えると残りの断面（リガメント）が小さくなり、応力が急激に集中するためです。実際の設計では $a/W < 0.5$ の範囲で使用することが推奨されます
-
き裂形状による差異が顕著: 同じ $a/W$ でも、中央き裂と片側き裂では $Y$ の値が大きく異なります。片側き裂は拘束条件が非対称であるため、常に中央き裂よりも大きな $K_I$ を与えます。構造設計では、き裂の形状に応じた適切な形状係数を選択することが重要です
-
SENB（三点曲げ）試験片は独自の形状係数パターンを示す: SENB は破壊靱性試験で最も一般的に使用される試験片形状であり、曲げ荷重による応力分布の非一様性が形状係数に反映されています
Python: パリス則によるき裂進展シミュレーション
パリス則を用いて、初期き裂が繰返し荷重によって臨界長さまで成長するプロセスを数値的にシミュレーションします。ここではオイラー法による逐次積分を行い、き裂長さと応力拡大係数の時間発展を追跡します。
import numpy as np
import matplotlib.pyplot as plt
# 材料・荷重パラメータ
E = 200e9 # ヤング率 [Pa]（鋼）
K_IC = 60.0 # 破壊靱性 [MPa*sqrt(m)]
sigma_max = 150.0 # 最大応力 [MPa]
sigma_min = 30.0 # 最小応力 [MPa]
Delta_sigma = sigma_max - sigma_min # 応力範囲 [MPa]
Y = 1.12 # 形状係数（片側き裂）
C = 1.0e-11 # パリス則定数 [m/cycle / (MPa*sqrt(m))^m]
m = 3.0 # パリス則指数
# 初期条件と臨界条件
a_0 = 1.0e-3 # 初期き裂長さ [m]（1 mm）
a_c = (1 / np.pi) * (K_IC / (Y * sigma_max))**2 # 臨界き裂長さ [m]
print(f"Critical crack length: a_c = {a_c*1e3:.2f} mm")
# 数値積分（オイラー法）
dN = 100 # 1ステップあたりのサイクル数
N_max = 5_000_000 # 最大サイクル数
N_list = [0]
a_list = [a_0]
K_list = [Y * sigma_max * np.sqrt(np.pi * a_0)]
a = a_0
N = 0
while a < a_c and N < N_max:
# 現在の応力拡大係数範囲
Delta_K = Y * Delta_sigma * np.sqrt(np.pi * a)
# パリス則によるき裂進展速度
da_dN = C * Delta_K**m
# き裂長さの更新
a += da_dN * dN
N += dN
# 記録
N_list.append(N)
a_list.append(a)
K_list.append(Y * sigma_max * np.sqrt(np.pi * a))
N_arr = np.array(N_list)
a_arr = np.array(a_list) * 1e3 # mm に変換
K_arr = np.array(K_list)
# 解析解（定数 Y の場合）
if m != 2:
N_f_analytical = (1 / (C * (Y * Delta_sigma)**m * np.pi**(m/2))
* 1 / (m/2 - 1)
* (a_0**(1 - m/2) - a_c**(1 - m/2)))
print(f"Analytical fatigue life: N_f = {N_f_analytical:.0f} cycles")
print(f"Numerical fatigue life: N_f = {N_arr[-1]:.0f} cycles")
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# 左: き裂長さ vs サイクル数
axes[0].plot(N_arr / 1e6, a_arr, 'b-', linewidth=2)
axes[0].axhline(y=a_c * 1e3, color='r', linestyle='--', linewidth=1.5,
label=f'$a_c$ = {a_c*1e3:.1f} mm')
axes[0].set_xlabel('Number of Cycles $N$ [$\\times 10^6$]', fontsize=13)
axes[0].set_ylabel('Crack Length $a$ [mm]', fontsize=13)
axes[0].set_title('Fatigue Crack Growth', fontsize=14)
axes[0].legend(fontsize=12)
axes[0].grid(True, alpha=0.3)
# 右: K_max vs サイクル数
axes[1].plot(N_arr / 1e6, K_arr, 'r-', linewidth=2)
axes[1].axhline(y=K_IC, color='k', linestyle='--', linewidth=1.5,
label=f'$K_{{IC}}$ = {K_IC:.0f} MPa$\\sqrt{{m}}$')
axes[1].set_xlabel('Number of Cycles $N$ [$\\times 10^6$]', fontsize=13)
axes[1].set_ylabel('$K_{\\max}$ [MPa$\\sqrt{m}$]', fontsize=13)
axes[1].set_title('Stress Intensity Factor Evolution', fontsize=14)
axes[1].legend(fontsize=12)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このシミュレーション結果から、疲労き裂進展の重要な特徴が明瞭に確認できます。
-
き裂進展の加速性: 左のグラフでは、初期段階ではき裂進展がほとんど見えないほど緩やかですが、寿命の終盤で急激にき裂が成長しています。これはパリス則の非線形性（$da/dN \propto (\Delta K)^m \propto a^{m/2}$）に起因します。き裂が長くなるほど $\Delta K$ が増大し、それによってさらにき裂進展速度が加速するという正のフィードバックが働いています
-
寿命の大部分は初期段階で消費される: 例えば全寿命の 80% の時点ではき裂長さはまだ臨界長さの半分にも達していないことがわかります。これは先ほど理論的に予測した通りで、$m > 2$ の場合に残存寿命が初期き裂長さに支配されることを反映しています
-
最終破壊は突然起こる: 右のグラフでは $K_{\max}$ が $K_{IC}$ に達する直前で急峻な上昇を示しています。実構造物において「何の前兆もなく突然壊れた」と報告される疲労破壊の多くは、まさにこの加速的なき裂進展の最終段階で起こっています
-
解析解と数値解の一致: 出力された解析解と数値積分の結果がほぼ一致しており、数値計算の妥当性が確認できます
次に、応力比 $R = \sigma_{\min}/\sigma_{\max}$ の影響を考慮した、より実践的なシミュレーションを実装します。
Python: 応力比の影響と寿命感度分析
応力比 $R$ やパリス則の指数 $m$ が疲労寿命に与える影響を系統的に調べます。
import numpy as np
import matplotlib.pyplot as plt
# 共通パラメータ
E = 200e9
K_IC = 60.0
Y = 1.12
a_0 = 1.0e-3
C = 1.0e-11
m = 3.0
def fatigue_life(sigma_max, R, C, m, Y, K_IC, a_0, dN=100, N_max=20_000_000):
"""パリス則によるき裂進展シミュレーション"""
sigma_min = R * sigma_max
Delta_sigma = sigma_max - sigma_min
a_c = (1 / np.pi) * (K_IC / (Y * sigma_max))**2
if a_0 >= a_c:
return 0, [], []
a = a_0
N = 0
N_hist, a_hist = [0], [a_0]
while a < a_c and N < N_max:
Delta_K = Y * Delta_sigma * np.sqrt(np.pi * a)
da_dN = C * Delta_K**m
a += da_dN * dN
N += dN
N_hist.append(N)
a_hist.append(a)
return N, N_hist, a_hist
# --- (a) 応力比 R の影響 ---
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
sigma_max_fixed = 150.0
R_values = [0.0, 0.1, 0.3, 0.5]
colors_R = ['#e74c3c', '#3498db', '#2ecc71', '#9b59b6']
for R, color in zip(R_values, colors_R):
Nf, N_hist, a_hist = fatigue_life(sigma_max_fixed, R, C, m, Y, K_IC, a_0)
axes[0].plot(np.array(N_hist) / 1e6, np.array(a_hist) * 1e3,
color=color, linewidth=2,
label=f'R = {R:.1f}, $N_f$ = {Nf/1e6:.2f}M')
axes[0].set_xlabel('Number of Cycles [$\\times 10^6$]', fontsize=13)
axes[0].set_ylabel('Crack Length [mm]', fontsize=13)
axes[0].set_title('Effect of Stress Ratio $R$', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# --- (b) パリス則指数 m の影響 ---
m_values = [2.0, 2.5, 3.0, 3.5, 4.0]
colors_m = ['#1abc9c', '#3498db', '#e74c3c', '#9b59b6', '#f39c12']
R_fixed = 0.1
for m_val, color in zip(m_values, colors_m):
Nf, N_hist, a_hist = fatigue_life(sigma_max_fixed, R_fixed, C, m_val, Y, K_IC, a_0)
axes[1].plot(np.array(N_hist) / 1e6, np.array(a_hist) * 1e3,
color=color, linewidth=2,
label=f'm = {m_val:.1f}, $N_f$ = {Nf/1e6:.2f}M')
axes[1].set_xlabel('Number of Cycles [$\\times 10^6$]', fontsize=13)
axes[1].set_ylabel('Crack Length [mm]', fontsize=13)
axes[1].set_title('Effect of Paris Exponent $m$', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
この感度分析の結果から、設計上の重要な知見が得られます。
-
応力比 $R$ の影響: 左のグラフは、$R$ が増加するほど（すなわち平均応力が増加するほど）疲労寿命が短くなることを示しています。$R = 0$（片振り）と $R = 0.5$ では寿命に数倍の差が生じています。これは $R$ が大きいと $\Delta\sigma$ が減少するものの、$K_{\max}$ が大きくなるため、パリス則の基本形では応力範囲の効果が直接反映されています。実際の材料ではウォーカー（Walker）の修正式やフォーマン（Forman）の式で平均応力効果をより精密に扱います
-
パリス則指数 $m$ の影響: 右のグラフでは、$m$ が1増えるだけで寿命が桁違いに変化しています。$m = 2$ では非常に長い寿命を示すのに対し、$m = 4$ では短い寿命で急速に破壊に至ります。$m$ の精密な測定がき裂進展寿命予測の信頼性に直結することがわかります
-
き裂進展曲線の形状の違い: $m$ が大きいほど、き裂進展の加速性が顕著になり、最終破壊前の急激な成長がより急峻になっています。$m = 2$ では比較的穏やかな加速を示すのに対し、$m = 4$ ではほぼ直角に近い急上昇が見られます
Python: $G$-$K$ 関係の検証と破壊条件マップ
最後に、エネルギー解放率 $G$ と応力拡大係数 $K$ の関係を図示し、破壊条件を応力-き裂長さ平面上でマップします。
import numpy as np
import matplotlib.pyplot as plt
# 材料パラメータ（Al 7075-T6）
E = 71e9 # ヤング率 [Pa]
nu = 0.33 # ポアソン比
K_IC = 27.0 # 破壊靱性 [MPa*sqrt(m)]
sigma_Y = 503.0 # 降伏応力 [MPa]
Y = 1.0 # 無限板
# G-K 関係の検証
K_range = np.linspace(0, 40, 200) # [MPa*sqrt(m)]
# 平面応力
G_plane_stress = K_range**2 / (E / 1e6) # [kJ/m^2] (K in MPa*sqrt(m))
# 平面ひずみ
G_plane_strain = K_range**2 * (1 - nu**2) / (E / 1e6)
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# 左: G vs K の関係
axes[0].plot(K_range, G_plane_stress * 1e3, 'b-', linewidth=2.5,
label='Plane stress: $G = K_I^2/E$')
axes[0].plot(K_range, G_plane_strain * 1e3, 'r--', linewidth=2.5,
label=r'Plane strain: $G = K_I^2(1-\nu^2)/E$')
axes[0].axvline(x=K_IC, color='k', linestyle=':', linewidth=1.5,
label=f'$K_{{IC}}$ = {K_IC} MPa$\\sqrt{{m}}$')
# G_IC の計算
G_IC_ps = K_IC**2 / (E / 1e6) * 1e3
G_IC_pe = K_IC**2 * (1 - nu**2) / (E / 1e6) * 1e3
axes[0].plot(K_IC, G_IC_ps, 'bo', markersize=10)
axes[0].plot(K_IC, G_IC_pe, 'rs', markersize=10)
axes[0].set_xlabel('$K_I$ [MPa$\\sqrt{m}$]', fontsize=13)
axes[0].set_ylabel('$G$ [J/m$^2$]', fontsize=13)
axes[0].set_title('Energy Release Rate vs Stress Intensity Factor', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# 右: 破壊条件マップ（応力 vs き裂長さ）
a_range = np.linspace(0.1e-3, 50e-3, 500) # [m]
sigma_crit = K_IC / (Y * np.sqrt(np.pi * a_range)) # [MPa]
# 降伏条件
sigma_yield_line = np.full_like(a_range, sigma_Y)
axes[1].fill_between(a_range * 1e3, sigma_crit, 600,
alpha=0.15, color='red', label='Fracture region ($K_I \\geq K_{IC}$)')
axes[1].fill_between(a_range * 1e3, sigma_yield_line, 600,
alpha=0.1, color='orange', label='Yield region ($\\sigma \\geq \\sigma_Y$)')
axes[1].plot(a_range * 1e3, sigma_crit, 'r-', linewidth=2.5,
label=f'$K_I = K_{{IC}}$ ({K_IC} MPa$\\sqrt{{m}}$)')
axes[1].axhline(y=sigma_Y, color='orange', linestyle='--', linewidth=2,
label=f'$\\sigma_Y$ = {sigma_Y} MPa')
axes[1].fill_between(a_range * 1e3, 0, np.minimum(sigma_crit, sigma_Y),
alpha=0.1, color='green')
axes[1].text(25, 100, 'SAFE\nREGION', fontsize=14, ha='center',
color='green', fontweight='bold')
axes[1].set_xlabel('Crack Length $a$ [mm]', fontsize=13)
axes[1].set_ylabel('Applied Stress $\\sigma$ [MPa]', fontsize=13)
axes[1].set_title('Failure Assessment Diagram (Al 7075-T6)', fontsize=14)
axes[1].set_xlim([0, 50])
axes[1].set_ylim([0, 600])
axes[1].legend(fontsize=10, loc='upper right')
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このグラフから、破壊力学の設計応用に直結する重要な情報が読み取れます。
-
$G$-$K$ 関係は放物線状: 左のグラフでは $G$ が $K_I^2$ に比例する放物線関係が確認できます。平面応力と平面ひずみでは $(1-\nu^2) = 0.89$ の係数差があり、平面ひずみ条件の方が同じ $K_I$ に対して $G$ がやや小さくなります。この差はポアソン比 $\nu$ が大きいほど顕著になります
-
破壊条件マップの実用的な読み方: 右のグラフ（破壊評価線図, Failure Assessment Diagram）は、構造設計者にとって最も重要な図の一つです。横軸がき裂長さ、縦軸が作用応力で、赤い曲線より上の領域が「破壊する条件」を表しています。任意の応力 $\sigma$ に対して許容き裂長さを読み取ったり、既知のき裂長さに対する安全応力を読み取ったりすることができます
-
降伏条件との競合: 短いき裂では破壊条件よりも先に降伏条件に達します。これは、小さなき裂の場合は $K_I$ が $K_{IC}$ に達する前に全面降伏（general yielding）が起こることを意味しています。破壊力学が適用できるのは、き裂先端の塑性域が小さい「小規模降伏」（small-scale yielding）の条件下に限られます
-
Al 7075-T6 の特性: 右のグラフから、例えば $\sigma = 200$ MPa の作用応力下では、約 5.8 mm の片側き裂が臨界長さとなることがわかります。航空機の構造検査では、この値よりも十分小さなき裂を確実に検出できる非破壊検査手法が必要となります
小規模降伏条件と塑性域
き裂先端の塑性域サイズ
ここまでの議論は全て「線形弾性」の枠組みの中で行ってきました。しかし、実際の材料ではき裂先端の応力が降伏応力 $\sigma_Y$ を超えると塑性変形が生じます。応力とひずみ で学んだ応力-ひずみ曲線において、弾性限度を超えた領域に相当します。
$\theta = 0$ 方向（き裂の前方）で $\sigma_{yy} = \sigma_Y$ となる距離 $r_Y$ を求めると、塑性域の大きさの推定が得られます。
$$ \begin{equation} \frac{K_I}{\sqrt{2\pi r_Y}} = \sigma_Y \end{equation} $$
これを $r_Y$ について解くと、
$$ \begin{equation} r_Y = \frac{1}{2\pi}\left(\frac{K_I}{\sigma_Y}\right)^2 \end{equation} $$
平面応力条件ではこの推定値がほぼ正しく、平面ひずみ条件では三軸応力状態の拘束効果により塑性域は約 $1/3$ に縮小します。
$$ \begin{equation} r_Y^{\text{plane strain}} \approx \frac{1}{6\pi}\left(\frac{K_I}{\sigma_Y}\right)^2 \end{equation} $$
小規模降伏の条件
線形弾性破壊力学（LEFM）が有効であるためには、塑性域がき裂長さ $a$ やリガメント（$W – a$）に比べて十分小さい小規模降伏（small-scale yielding, SSY）条件が成立している必要があります。一般的な目安として、
$$ \begin{equation} r_Y \ll a, \quad r_Y \ll (W – a) \end{equation} $$
ASTM 規格では $a \geq 2.5(K_{IC}/\sigma_Y)^2$ を要求しており、これは塑性域がき裂長さの約 $1/25$ 以下であることに対応します。
小規模降伏条件が成立しない場合（大規模降伏）には、弾塑性破壊力学のパラメータ — $J$ 積分（$J$-integral）やき裂先端開口変位 CTOD（Crack Tip Opening Displacement）— を使う必要があります。これらは線形弾性破壊力学の自然な拡張であり、$J$ は $G$ の弾塑性版、CTOD はき裂面の開口変位に基づくパラメータです。
破壊力学の体系を一通り学んだところで、ここまでの内容を整理し、実務への応用の視点も含めてまとめましょう。
破壊力学の実務応用 — 損傷許容設計
損傷許容設計の思想
破壊力学は、航空宇宙分野を中心に「損傷許容設計」（damage tolerance design）の基盤として発展してきました。従来の「安全係数設計」が「欠陥のない理想的な材料」を前提としたのに対し、損傷許容設計は「構造物には必ず欠陥が存在する」という前提に立ちます。
損傷許容設計の基本フレームワークは以下の通りです。
- 初期欠陥の想定: 製造後の構造物に非破壊検査の検出限界に相当するき裂が存在すると仮定します
- き裂進展の予測: パリス則などを用いて、繰返し荷重下でのき裂進展速度を計算します
- 検査間隔の設定: き裂が臨界長さに達するまでの寿命の半分を検査間隔とし、定期的に非破壊検査を実施します
- 残存強度の確認: き裂が存在する状態でも、設計制限荷重に対して十分な残存強度を有することを確認します
主応力とモールの応力円 で学んだ応力解析と破壊力学を組み合わせることで、多軸応力状態にあるき裂の安定性も評価できます。
航空機の損傷許容認証
航空機の構造認証（FAR 25.571, EASA CS 25.571）では、構造部材を以下のカテゴリに分類して管理します。
- フェイルセーフ（fail-safe）構造: 一つの部材が破壊しても、残りの部材で荷重を支えられる冗長設計
- 損傷許容（damage tolerant）構造: き裂が存在しても検査間隔内に臨界サイズに達しないことを実証
いずれの場合も、破壊力学パラメータ（$K_I$, $\Delta K$, $K_{IC}$）に基づく定量的な評価が求められます。
まとめ
本記事では、破壊力学の基礎を体系的に解説しました。
-
従来の強度設計の限界: き裂先端では応力が無限大に発散するため、応力基準だけでは破壊を予測できません。この問題を解決するために破壊力学が生まれました
-
グリフィスのエネルギー基準: き裂進展を「ひずみエネルギーの解放」と「表面エネルギーの増加」の収支として捉える画期的な発想で、臨界き裂長さ $a_c = 2E\gamma_{\text{eff}}/(\pi\sigma^2)$ が得られます
-
応力拡大係数 $K_I$: き裂先端の応力場の強さを表すパラメータで、$K_I = Y\sigma\sqrt{\pi a}$ と定義されます。き裂先端の応力場は $1/\sqrt{r}$ の特異性を持ち、その角度分布は普遍的で、強さは $K_I$ のみで決まります（$K$ 支配の原理）
-
破壊靱性 $K_{IC}$: 材料のき裂に対する抵抗力を表す材料定数で、$K_I \geq K_{IC}$ のとき不安定破壊が生じます
-
$G$ と $K$ の関係: $G = K_I^2/E$（平面応力）または $G = K_I^2(1-\nu^2)/E$（平面ひずみ）という等価関係があり、エネルギー的視点と応力場の視点が統一されます
-
パリス則: 疲労き裂進展速度 $da/dN = C(\Delta K)^m$ により、初期き裂から臨界き裂長さまでの残存寿命を予測できます
破壊力学は、安全な構造設計のための不可欠な学問であり、航空宇宙、原子力、自動車、土木など幅広い分野の根幹を支えています。本記事で扱った線形弾性破壊力学（LEFM）は、大規模降伏条件では弾塑性破壊力学（$J$ 積分, CTOD）に拡張されます。
次のステップとして、以下の記事も参考にしてください。
- 応力集中 — き裂のない欠陥（穴、切欠き）に対する応力評価
- 疲労破壊 — S-N 曲線による疲労寿命評価の基礎
- 平面応力と平面ひずみ — 2次元弾性問題の定式化と破壊力学への応用