金属の棒を手で繰り返し折り曲げると、1回では何も起きなかったのに、何十回か繰り返すとあるとき突然ポキッと折れます。この現象は日常的に経験できますが、工学の世界ではきわめて深刻な破壊モードです。航空機の翼、橋梁の桁、自動車のサスペンション、風力発電のブレード――これらの構造部材は、静的な荷重では安全な応力レベルであっても、繰り返し荷重を何百万回、何千万回と受けることで破壊に至ることがあります。これが疲労破壊（fatigue failure）です。
疲労破壊が恐ろしいのは、静的強度の半分程度の応力でも発生しうること、そして破壊の直前まで目に見える変形がほとんど現れないことです。実際、機械部品の破損原因の約80〜90%は疲労破壊であるとも言われています。1954年のデ・ハビランド コメット連続墜落事故では、客室窓角部の応力集中箇所から疲労き裂が進展し、機体が空中分解しました。この事故は航空工学における疲労設計の重要性を世界に知らしめた歴史的事件です。
疲労破壊の理論を理解すると、以下のような幅広い分野で活用できます。
- 航空宇宙設計: 機体構造や翼の繰り返し荷重に対する寿命予測と検査間隔の設定
- 自動車設計: サスペンション、エンジン部品、車体フレームの耐久性評価
- 土木・建築: 橋梁、鉄道レール、クレーンなどの繰り返し荷重を受ける構造物の安全性評価
- エネルギー産業: 風力タービンのブレードや原子力配管の経年劣化管理
- 機械設計全般: 回転軸、歯車、ばねなどあらゆる繰り返し荷重を受ける部品の寿命設計
本記事の内容
- 疲労破壊の概要と静的強度との本質的な違い
- 疲労き裂の発生・進展メカニズム（ミクロ的視点）
- S-N曲線（ウェーラー曲線）の読み方と数学的表現
- 疲労限度（耐久限度）の概念と設計上の意味
- 平均応力の影響と修正グッドマン線図・ゾーダーベルク線図
- マイナーの累積損傷則（線形損傷則）
- 応力集中と疲労の関係（切欠き係数）
- Python実装: S-N曲線の描画・フィッティング、修正グッドマン線図、マイナー則による累積損傷計算、ランダム荷重下の疲労寿命推定
前提知識
この記事を読む前に、以下の記事に目を通しておくと理解がスムーズです。
- 応力とひずみの定義 — 垂直応力・せん断応力の基本的な定義と単位
- フックの法則と弾性係数 — 弾性域における応力とひずみの線形関係
- 応力集中 — 穴や切欠き周りの応力が局所的に増大する現象
- 対数とべき乗 — S-N曲線の両対数表示に必要な対数の性質
疲労破壊とは — 静的強度では安全なのに壊れる
静的破壊との違い
まず、疲労破壊を正しく位置づけるために、静的破壊と比較してみましょう。
応力とひずみ で学んだように、材料に力を加えると応力が生じます。材料には引張強さ $\sigma_B$（ultimate tensile strength）があり、応力がこの値を超えると破壊に至ります。これが静的破壊です。静的破壊では、荷重を徐々に増加させていくと、弾性変形→塑性変形→くびれ→破断という明確なプロセスを経ます。破壊に至る前には大きな塑性変形が生じるため、「壊れる前兆」がはっきり現れます。
一方、疲労破壊は本質的に異なります。応力が引張強さどころか降伏応力よりも小さい値であっても、その応力が繰り返し加わることで破壊が起こるのです。しかも、破壊面を観察すると、塑性変形はほとんど見られません。つまり、延性材料であっても脆性的に破壊するのが疲労破壊の特徴です。
両者の違いを表にまとめます。
| 特性 | 静的破壊 | 疲労破壊 |
|---|---|---|
| 荷重の種類 | 単調増加荷重 | 繰り返し荷重（変動荷重） |
| 破壊応力 | 引張強さ $\sigma_B$ 以上 | $\sigma_B$ よりはるかに低い応力でも発生 |
| 破壊前の変形 | 大きな塑性変形（延性材料の場合） | ほとんど変形なし（見かけ上、脆性的） |
| 破壊のプロセス | 弾性→塑性→くびれ→破断 | き裂発生→き裂進展→最終破断 |
| 破面の特徴 | カップ＆コーン型（延性材料） | ビーチマーク（貝殻模様）＋最終脆性破面 |
| 予測の困難さ | 材料試験データから比較的容易 | 多くの因子が関与し予測が複雑 |
疲労破壊が「沈黙の殺し屋」と呼ばれる所以は、目に見える前兆がほとんどないまま突然破壊に至ることにあります。だからこそ、疲労のメカニズムを理解し、適切な寿命予測手法を身につけることが工学において決定的に重要なのです。
ここで自然な疑問が浮かびます――繰り返し荷重によって材料内部では何が起きているのでしょうか。次に、疲労き裂の発生と進展のメカニズムを見ていきましょう。
疲労き裂の発生と進展メカニズム
ミクロ的な視点 — すべり帯からき裂へ
疲労破壊がなぜ低い応力で起こるのかを理解するには、材料のミクロ構造レベルまで掘り下げる必要があります。
金属材料は多数の結晶粒（グレイン）から構成されています。各結晶粒には特定の結晶方位があり、外力が作用すると結晶面に沿ったすべり（slip）が起こります。1回の荷重サイクルでは、このすべりの量はごくわずかです。しかし、荷重が繰り返されると、すべりが同じ面に沿って累積していきます。
この累積すべりのプロセスをもう少し詳しく見てみましょう。引張荷重が作用すると、結晶内の特定のすべり面に沿って転位が移動し、表面にわずかな段差（押し出し, extrusion）と溝（引き込み, intrusion）が形成されます。荷重が除荷されてから逆方向に加わると、逆方向のすべりが生じますが、酸化膜の形成や転位の絡み合いにより、すべりは完全には元に戻りません。このような不可逆的なすべりの蓄積が、サイクルを重ねるごとに進行します。
やがて、表面の引き込み部分が深くなり、微小なき裂（マイクロクラック）が形成されます。これが疲労き裂の発生段階（initiation stage）です。
き裂の進展 — 3つの段階
疲労破壊のプロセスは、大きく3つの段階に分けられます。
第I段階（き裂発生・初期成長）: 表面の入り込みから発生した微小き裂は、最初は最大せん断応力方向（荷重方向に対して約45度）に沿って成長します。この段階のき裂は非常に小さく（数結晶粒のスケール）、結晶粒界で一時的に停止することもあります。
第II段階（安定き裂進展）: き裂がある程度の大きさに達すると、成長方向が変化し、最大引張応力に垂直な方向（荷重方向に対して90度）に沿って安定的に進展していきます。この段階では、各荷重サイクルごとにき裂が微小量だけ進展します。破面を電子顕微鏡で観察すると、1サイクルごとの進展量に対応する縞模様（ストライエーション, striation）が見られます。この段階の進展速度は破壊力学入門で解説する応力拡大係数によって定量的に記述できます（Paris則）。
第III段階（最終破断）: き裂が臨界長さに達すると、残った断面積ではもはや荷重を支えきれなくなり、静的な過負荷によって一気に破断します。この最終破断領域は粗い破面を呈し、延性材料ではディンプル模様、脆性材料ではクリベージ模様が観察されます。
破面の特徴 — ビーチマーク
実際の疲労破面を肉眼で観察すると、同心円状の縞模様が見られることがあります。これをビーチマーク（beach mark）または貝殻模様（clamshell mark）と呼びます。ビーチマークは、荷重条件が変化した際（運転条件の切替え、休止期間など）にき裂先端の形状や酸化状態が変わることで形成されます。ビーチマークの中心がき裂の起点であり、そこから同心円状に広がる模様がき裂の進展方向を示しています。
破面解析（フラクトグラフィー）では、ビーチマークの模様と最終破断領域の比率から、破壊に至った荷重の大きさや作用した繰り返し回数を推定することができます。
ここまでで、疲労破壊のミクロ的なメカニズムを理解しました。では、このような疲労破壊を工学的に定量化するにはどうすればよいのでしょうか。19世紀にアウグスト・ウェーラーが発明した画期的なアプローチ、S-N曲線を見ていきましょう。
S-N曲線（ウェーラー曲線）
S-N曲線とは
疲労寿命を定量的に評価する最も基本的なツールがS-N曲線（S-N curve）です。Sは応力振幅（stress amplitude）、Nは破断繰り返し数（number of cycles to failure）を表します。この曲線は、1860年代にドイツの鉄道技師アウグスト・ウェーラー（August Wöhler）が車軸の疲労破壊を調べるために考案したもので、ウェーラー曲線とも呼ばれます。
S-N曲線を得るためには、実験室で以下のような疲労試験を行います。
- 同じ材料・形状の試験片を多数用意する
- 各試験片に異なる応力振幅の繰り返し荷重（通常は完全両振り: 平均応力ゼロ）を与える
- 破断するまでの繰り返し数 $N$ を記録する
- 応力振幅 $S_a$ を縦軸に、$N$ を横軸にプロットする
繰り返し応力の定義
S-N曲線を正しく読むためには、繰り返し応力を特徴づけるパラメータを理解する必要があります。時間とともに変動する応力 $\sigma(t)$ を考えましょう。
最大応力 $\sigma_{\max}$ と最小応力 $\sigma_{\min}$ を用いて、以下の量が定義されます。
応力振幅（stress amplitude）は、応力の変動幅の半分です。
$$ \begin{equation} \sigma_a = \frac{\sigma_{\max} – \sigma_{\min}}{2} \end{equation} $$
平均応力（mean stress）は、最大応力と最小応力の中央値です。
$$ \begin{equation} \sigma_m = \frac{\sigma_{\max} + \sigma_{\min}}{2} \end{equation} $$
応力範囲（stress range）は、最大応力と最小応力の差です。
$$ \begin{equation} \Delta\sigma = \sigma_{\max} – \sigma_{\min} = 2\sigma_a \end{equation} $$
応力比（stress ratio）は、最小応力と最大応力の比です。
$$ \begin{equation} R = \frac{\sigma_{\min}}{\sigma_{\max}} \end{equation} $$
よく使われる荷重条件を整理します。
| 荷重条件 | 応力比 $R$ | 平均応力 $\sigma_m$ | 説明 |
|---|---|---|---|
| 完全両振り | $R = -1$ | $\sigma_m = 0$ | 引張と圧縮が対称。標準的な疲労試験 |
| 片振り（引張） | $R = 0$ | $\sigma_m = \sigma_a$ | 0からある最大値まで変動 |
| 脈動（引張） | $0 < R < 1$ | $\sigma_m > \sigma_a$ | 常に引張側で変動 |
標準的なS-N曲線は $R = -1$（完全両振り）の条件で取得されます。平均応力の影響については後のセクションで扱います。
S-N曲線の数学的表現 — バスキンの式
S-N曲線を対数軸で描くと、多くの金属材料で近似的に直線関係が得られます。これをバスキンの式（Basquin’s equation）として表すと、次のようになります。
$$ \begin{equation} \sigma_a = \sigma_f’ (2N)^b \end{equation} $$
ここで $\sigma_f’$ は疲労強度係数（fatigue strength coefficient）、$b$ は疲労強度指数（fatigue strength exponent、通常 $-0.05$ から $-0.12$ 程度の負の値）です。$2N$ は応力の反転回数（reversal）を表し、1サイクルあたり2回の反転があるため繰り返し数 $N$ の2倍になります。
両辺の対数をとると、次のような線形関係が得られます。
$$ \begin{equation} \log \sigma_a = \log \sigma_f’ + b \log(2N) \end{equation} $$
この式は、両対数グラフ上で S-N 曲線が傾き $b$、切片 $\log \sigma_f’$ の直線になることを意味しています。多くの金属材料について $\sigma_f’$ は引張強さ $\sigma_B$ にほぼ等しく、$b$ は材料によって $-0.05$ 〜 $-0.12$ の範囲の値をとります。
S-N曲線の形状は材料によって異なりますが、鉄鋼材料とアルミニウム合金では大きく異なる特徴があります。この違いを理解する鍵が「疲労限度」の概念です。
疲労限度（耐久限度）
鉄鋼材料の疲労限度
S-N曲線を低応力領域まで延長していくと、鉄鋼材料では非常に興味深い現象が観察されます。ある応力振幅以下になると、繰り返し数を $10^6$ 〜 $10^7$ 回以上に増やしても破断しなくなるのです。このそれ以下では破断が起こらない応力振幅の限界値を疲労限度（fatigue limit）または耐久限度（endurance limit）と呼び、$\sigma_w$ あるいは $S_e$ で表します。
疲労限度が存在する物理的なメカニズムは、転位の固着（ロッキング）と関係しています。鉄鋼材料に含まれる炭素や窒素などの侵入型固溶原子が転位の周囲に偏析してコットレル雰囲気を形成し、転位の運動を阻害します。応力振幅が十分に小さい場合、転位は固着されたまま動くことができず、すべり帯が形成されないため、疲労き裂が発生しません。
実用的な目安として、鉄鋼材料の疲労限度は引張強さのおよそ半分です。
$$ \begin{equation} \sigma_w \approx 0.5 \, \sigma_B \quad (\text{鉄鋼材料、} \sigma_B \lesssim 1400 \, \text{MPa の範囲で}) \end{equation} $$
ただし、引張強さが約 1400 MPa を超える高強度鋼では、この比例関係が崩れ、疲労限度は頭打ちになる傾向があります。これは高強度鋼ほど介在物や表面欠陥に敏感になるためです。
非鉄金属の疲労特性
一方、アルミニウム合金や銅合金などの多くの非鉄金属では、明確な疲労限度が存在しません。S-N曲線はどこまでも右下がりに続き、応力振幅がどんなに小さくても、十分な繰り返し数があれば破断に至る可能性があります。
この場合、設計上は特定の繰り返し数（たとえば $N = 10^7$ や $N = 10^8$）における強度を「実用上の疲労限度」として定義し、時間強度（fatigue strength at $N$ cycles）と呼びます。アルミニウム合金の場合、$N = 5 \times 10^8$ における疲労強度は引張強さの約 $0.3$ 〜 $0.4$ 倍程度です。
| 材料 | 疲労限度の有無 | 目安 |
|---|---|---|
| 炭素鋼・合金鋼 | あり（$10^6$ 〜 $10^7$ 回付近） | $\sigma_w \approx 0.5 \sigma_B$ |
| ステンレス鋼 | あり（ただし腐食環境では消失しうる） | $\sigma_w \approx 0.4 \sigma_B$ |
| アルミニウム合金 | なし | $N = 5\times10^8$ で $\approx 0.35 \sigma_B$ |
| 銅合金 | なし | $N = 10^8$ で $\approx 0.3 \sigma_B$ |
| チタン合金 | あり | $\sigma_w \approx 0.5 \sigma_B$ |
設計の現場では、疲労限度以下の応力で部品を使用する「無限寿命設計」と、有限のサイクル数での使用を前提とする「有限寿命設計」を使い分けます。航空機の構造部材のように検査と交換が可能なものでは有限寿命設計が採用され、定期的な検査と部品の交換が義務付けられています。
ここまではすべて、平均応力がゼロ（$R = -1$、完全両振り）の条件で話を進めてきました。しかし現実の構造部材は、常に引張応力がかかった状態で変動荷重を受けるなど、平均応力がゼロでない場合がほとんどです。平均応力は疲労寿命にどのような影響を与えるのでしょうか。
平均応力の影響
なぜ平均応力が重要なのか
実際の機械部品は、ボルトの初期締付力、回転軸の自重による曲げ、圧力容器の内圧など、定常的な荷重（静荷重）が作用した状態で、追加の繰り返し荷重を受けます。このとき、応力波形の「中心線」がゼロからずれます。すなわち、$\sigma_m \neq 0$ です。
実験事実として、引張側の平均応力（$\sigma_m > 0$）は疲労寿命を減少させ、圧縮側の平均応力（$\sigma_m < 0$）は疲労寿命を増加させます。直感的には、引張の平均応力はき裂を開く方向に作用するため、き裂の発生と進展を促進すると理解できます。
この効果を定量的に表現する代表的な方法が、グッドマン線図とゾーダーベルク線図です。
修正グッドマン線図
修正グッドマン線図（modified Goodman diagram）は、平均応力と応力振幅の安全な組合せを可視化するための図です。
横軸に平均応力 $\sigma_m$、縦軸に応力振幅 $\sigma_a$ をとります。グッドマンの関係式は次のようになります。
$$ \begin{equation} \frac{\sigma_a}{\sigma_w} + \frac{\sigma_m}{\sigma_B} = 1 \end{equation} $$
ここで $\sigma_w$ は疲労限度（$R=-1$ の条件）、$\sigma_B$ は引張強さです。この式は、$\sigma_a$-$\sigma_m$ 平面上で $(\sigma_m, \sigma_a) = (0, \sigma_w)$ と $(\sigma_B, 0)$ を結ぶ直線を表しています。この直線の内側（原点側）が安全領域、外側が危険領域です。
修正グッドマンの式を応力振幅について解くと、平均応力 $\sigma_m$ が存在するときの等価的な疲労限度 $\sigma_{a,\text{allow}}$ が得られます。
$$ \begin{equation} \sigma_{a,\text{allow}} = \sigma_w \left(1 – \frac{\sigma_m}{\sigma_B}\right) \end{equation} $$
この式は明快です。平均応力が引張強さに近づくほど、許容される応力振幅は小さくなり、$\sigma_m = \sigma_B$ で許容振幅はゼロになります。すなわち、静的に破壊する応力に達してしまうということです。
ゾーダーベルク線図
ゾーダーベルク（Soderberg）の関係式は、グッドマンよりも保守的（安全側）な推定を与えます。引張強さ $\sigma_B$ の代わりに降伏応力 $\sigma_Y$ を用いる点が異なります。
$$ \begin{equation} \frac{\sigma_a}{\sigma_w} + \frac{\sigma_m}{\sigma_Y} = 1 \end{equation} $$
降伏応力は引張強さよりも小さいため（$\sigma_Y < \sigma_B$）、ゾーダーベルク線図はグッドマン線図の内側に位置し、より保守的な設計基準を与えます。安全が最優先される構造物の設計では、ゾーダーベルク線図が好まれることがあります。
ガーバー線図
ガーバー（Gerber）の放物線は、実験データにより近い近似を与えることが多く、次のように表されます。
$$ \begin{equation} \frac{\sigma_a}{\sigma_w} + \left(\frac{\sigma_m}{\sigma_B}\right)^2 = 1 \end{equation} $$
グッドマン線図が直線であるのに対し、ガーバー線図は上に凸の放物線です。多くの実験データはガーバーの放物線とグッドマン直線の間に分布するため、グッドマン線図は安全側の推定を与えるのが一般的です。
これら3つの基準を比較すると、保守性の順序は次のようになります。
$$ \text{ゾーダーベルク（最も保守的）} \subset \text{グッドマン} \subset \text{ガーバー（最も非保守的）} $$
設計の現場では、安全率の要求やデータの信頼性に応じて使い分けます。航空宇宙分野ではグッドマン線図が広く使われ、一般機械設計ではゾーダーベルク線図が用いられることが多いです。
ここまでで、一定振幅の繰り返し荷重に対する疲労評価の枠組みを学びました。しかし現実の構造物は、異なる大きさの応力が混在する複雑な荷重履歴を受けます。このような場合に疲労寿命をどう推定するかが次の課題です。
マイナーの累積損傷則
線形累積損傷の考え方
現実の構造物は、単一の応力振幅の繰り返しではなく、さまざまな大きさの応力が不規則に作用します。たとえば、航空機の翼は離着陸時の大きな荷重と巡航中の小さな乱気流荷重が組み合わさった複雑な荷重履歴を受けます。
このような変動振幅荷重下での疲労寿命を推定するための最も基本的な手法が、マイナーの累積損傷則（Miner’s rule）、または線形累積損傷則（linear damage accumulation rule, Palmgren-Miner rule）です。この手法は1924年にパルムグレンが提案し、1945年にマイナーが体系化しました。
マイナー則の基本的なアイデアは、非常にシンプルです。
各応力レベルでの損傷率（damage fraction）は、その応力レベルで実際に受けた繰り返し数と破断寿命の比に等しいと仮定します。応力振幅 $\sigma_{a,i}$ に対するS-N曲線上の破断寿命を $N_i$ とし、実際にその応力レベルで受けた繰り返し数を $n_i$ とすると、その応力レベルでの損傷率 $D_i$ は次のように定義されます。
$$ \begin{equation} D_i = \frac{n_i}{N_i} \end{equation} $$
各応力レベルの損傷率を合算し、その合計が1に達したときに破壊が起こるとします。
$$ \begin{equation} D = \sum_{i=1}^{k} \frac{n_i}{N_i} = 1 \quad \text{（破壊条件）} \end{equation} $$
この式の物理的な意味を具体例で確認しましょう。ある部品のS-N曲線から、応力振幅 200 MPa での破断寿命が $N_1 = 10^5$ 回、150 MPa での破断寿命が $N_2 = 10^6$ 回であるとします。この部品が 200 MPa を $n_1 = 5 \times 10^4$ 回、150 MPa を $n_2 = 5 \times 10^5$ 回受けた場合、累積損傷は次のようになります。
$$ D = \frac{5 \times 10^4}{10^5} + \frac{5 \times 10^5}{10^6} = 0.5 + 0.5 = 1.0 $$
$D = 1$ に達したので、ちょうど破壊に至ると予測されます。
マイナー則の限界
マイナー則は計算が簡単で直感的にもわかりやすいため広く使われていますが、いくつかの重要な限界があります。
荷重の順序効果を無視する: マイナー則では、高応力→低応力の順に荷重が作用しても、低応力→高応力の順でも、累積損傷は同じになります。しかし実験的には、高応力の後に低応力を負荷すると（High-Low loading）、マイナー則の予測より早く破壊する傾向があり、逆に低応力→高応力（Low-High loading）では寿命が長くなる傾向があります。これは、高応力でき裂先端の塑性域が形成され、後続の低応力でのき裂進展が促進されるためです。
応力間の相互作用を無視する: 各応力レベルでの損傷は独立に蓄積されると仮定しており、異なる応力レベル間の相互作用（たとえば過負荷によるき裂先端の鈍化やき裂閉口効果）を考慮しません。
実験値とのずれ: 実験では $D = 0.7$ 〜 $2.2$ 程度のばらつきで破壊が起こり、必ずしも $D = 1$ ちょうどではありません。安全側の設計のために $D < 0.5$ や $D < 0.7$ を破壊条件として用いることもあります。
これらの限界にもかかわらず、マイナー則は「まず第一近似として」広く使われている手法です。より精密な方法としては、非線形累積損傷則（Marco-Starkey則など）や、破壊力学に基づくき裂進展解析がありますが、それらは破壊力学入門で扱う発展的な内容です。
次に、疲労寿命に極めて大きな影響を及ぼす因子である応力集中について見ていきましょう。
応力集中と疲労
応力集中係数と疲労切欠き係数
応力集中 で学んだように、穴、溝、段付き部、キー溝などの形状不連続部（切欠き）が存在すると、局所的に応力が増大します。この応力増大の度合いを表すのが応力集中係数 $K_t$ です。
$$ \begin{equation} K_t = \frac{\sigma_{\max,\text{local}}}{\sigma_{\text{nom}}} \end{equation} $$
ここで $\sigma_{\max,\text{local}}$ は切欠き底の最大応力、$\sigma_{\text{nom}}$ は切欠きがない場合の公称応力です。
静的強度に対する応力集中の影響は材料の延性によって緩和されますが、疲労では応力集中の影響が非常に大きくなります。ただし、疲労における応力集中の影響は $K_t$ そのものではなく、疲労切欠き係数 $K_f$ で表されます。
$$ \begin{equation} K_f = \frac{\text{平滑材の疲労限度}}{\text{切欠き材の疲労限度}} = \frac{\sigma_{w,\text{smooth}}}{\sigma_{w,\text{notch}}} \end{equation} $$
一般に $1 \leq K_f \leq K_t$ であり、$K_f$ は $K_t$ 以下の値をとります。両者の関係は切欠き感度係数 $q$ で表されます。
$$ \begin{equation} q = \frac{K_f – 1}{K_t – 1} \end{equation} $$
$q = 0$ のとき $K_f = 1$ で切欠きの疲労への影響がなく、$q = 1$ のとき $K_f = K_t$ で応力集中の影響がそのまま疲労に反映されます。
切欠き感度係数 $q$ は材料と切欠きの形状（特に切欠き底の曲率半径 $\rho$）に依存します。一般的な傾向として次のことが知られています。
- 高強度材料ほど $q$ が大きい（切欠きに敏感）
- 切欠き底の曲率半径が小さいほど $q$ が大きい
- 鋳鉄のような内部欠陥の多い材料は $q$ が小さい
疲労強度に影響する各種修正係数
実際の部品の疲労限度を推定する際には、標準試験片で得られた疲労限度 $\sigma_w$ にさまざまな修正係数を乗じて、実部品の疲労限度 $\sigma_{w,\text{part}}$ を求めます。
$$ \begin{equation} \sigma_{w,\text{part}} = \frac{C_{\text{surf}} \cdot C_{\text{size}} \cdot C_{\text{rel}} \cdot C_{\text{temp}}}{K_f} \cdot \sigma_w \end{equation} $$
| 修正係数 | 記号 | 影響する因子 |
|---|---|---|
| 表面仕上げ係数 | $C_{\text{surf}}$ | 鏡面研磨→1.0、機械加工→0.7〜0.9、鍛造→0.3〜0.5 |
| 寸法係数 | $C_{\text{size}}$ | 大きい部品ほど小さい（表面積増加で欠陥確率増加） |
| 信頼度係数 | $C_{\text{rel}}$ | 50%信頼度→1.0、99%信頼度→0.814 |
| 温度係数 | $C_{\text{temp}}$ | 高温で疲労限度が低下（クリープと応力緩和 も参照） |
これらの修正係数は、実験データの統計的な分析に基づくものです。疲労データには本質的にばらつきが大きく、正規分布 に基づく統計的処理が重要になります。
以上で疲労破壊に関する理論的な枠組みを一通り学びました。ここからは、Python を使ってこれらの概念を実装し、可視化を通じて理解を深めていきましょう。
PythonによるS-N曲線の描画とフィッティング
まず、疲労試験データからS-N曲線を描画し、バスキンの式でフィッティングしてみましょう。ここでは仮想的な鉄鋼材料（S45C相当）の疲労試験データを用いて、S-N曲線の描き方と対数線形フィッティングの方法を示します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
# --- 仮想的な疲労試験データ（S45C相当の鋼材）---
# 応力振幅 [MPa] と破断繰り返し数
stress_amp = np.array([450, 400, 370, 340, 310, 290, 270, 260, 250, 245])
cycles_to_failure = np.array([
1.2e3, 5.0e3, 1.5e4, 5.0e4, 2.0e5, 6.0e5, 3.0e6, 8.0e6, 1.0e7, 1.0e7
])
# 疲労限度を推定（最後の2点は未破断と仮定）
fatigue_limit = 245 # MPa（推定値）
# バスキンの式でフィッティング: σ_a = σ_f' * (2N)^b
# 対数変換: log(σ_a) = log(σ_f') + b * log(2N)
# 疲労限度以上のデータのみ使用
mask = stress_amp > fatigue_limit
log_reversals = np.log10(2 * cycles_to_failure[mask])
log_stress = np.log10(stress_amp[mask])
# 線形回帰で b と log(σ_f') を求める
coeffs = np.polyfit(log_reversals, log_stress, 1)
b = coeffs[0] # 疲労強度指数
log_sigma_f = coeffs[1] # log(σ_f')
sigma_f = 10**log_sigma_f # 疲労強度係数 [MPa]
print(f"疲労強度係数 σ_f' = {sigma_f:.1f} MPa")
print(f"疲労強度指数 b = {b:.4f}")
# フィッティング曲線の描画
N_fit = np.logspace(2, 8, 200)
stress_fit = sigma_f * (2 * N_fit)**b
# 疲労限度で下限をクリップ
stress_fit = np.maximum(stress_fit, fatigue_limit)
# 描画
fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(cycles_to_failure, stress_amp, c='red', s=80,
zorder=5, label='Experimental data')
ax.plot(N_fit, stress_fit, 'b-', linewidth=2, label="Basquin's fit")
ax.axhline(y=fatigue_limit, color='green', linestyle='--',
linewidth=1.5, label=f'Fatigue limit = {fatigue_limit} MPa')
ax.set_xscale('log')
ax.set_xlabel('Number of cycles to failure, N', fontsize=13)
ax.set_ylabel('Stress amplitude, $\\sigma_a$ [MPa]', fontsize=13)
ax.set_title('S-N Curve (Steel, S45C equivalent)', fontsize=14)
ax.legend(fontsize=11)
ax.set_xlim(1e2, 1e8)
ax.set_ylim(200, 500)
ax.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.show()
上のグラフから、S-N曲線の典型的な特徴がはっきりと読み取れます。
-
両対数軸上での直線性: 応力振幅 250 MPa 以上の領域では、データ点がバスキンのフィッティング直線によく乗っていることが確認できます。これは、バスキンの式 $\sigma_a = \sigma_f'(2N)^b$ が有限寿命領域で良い近似を与えることを意味しています。
-
疲労限度の存在: $10^6$ 〜 $10^7$ 回付近で曲線が水平になり、それ以下の応力振幅（245 MPa 以下）では破断が起こらなくなることが見てとれます。これは鉄鋼材料に特有の疲労限度の存在を示しています。
-
応力と寿命のトレードオフ: 応力振幅を 450 MPa から 270 MPa に40%下げるだけで、寿命は約 $10^3$ 回から $10^6$ 回へと1000倍以上に延びています。疲労設計において応力低減がいかに効果的かがわかります。
フィッティングで得られた疲労強度指数 $b$ の値は、一般的な鋼材の範囲（$-0.05$ 〜 $-0.12$）に入っていることを確認しましょう。
続いて、S-N曲線に加えてアルミニウム合金のデータも重ねて描画し、鉄鋼材料との疲労特性の違いを可視化してみましょう。
import numpy as np
import matplotlib.pyplot as plt
# --- 鉄鋼材料（S45C相当）---
N_steel = np.logspace(2, 8, 300)
sigma_f_steel = 900 # MPa（疲労強度係数）
b_steel = -0.085 # 疲労強度指数
fatigue_limit_steel = 245 # MPa
stress_steel = sigma_f_steel * (2 * N_steel)**b_steel
stress_steel = np.maximum(stress_steel, fatigue_limit_steel)
# --- アルミニウム合金（A2024-T4相当）---
N_al = np.logspace(2, 9, 300)
sigma_f_al = 680 # MPa
b_al = -0.075 # 疲労強度指数（疲労限度なし）
stress_al = sigma_f_al * (2 * N_al)**b_al
# 描画
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(N_steel, stress_steel, 'b-', linewidth=2.5, label='Steel (S45C)')
ax.plot(N_al, stress_al, 'r--', linewidth=2.5, label='Aluminum (A2024-T4)')
ax.axhline(y=fatigue_limit_steel, color='blue', linestyle=':',
alpha=0.5, label=f'Steel fatigue limit = {fatigue_limit_steel} MPa')
ax.set_xscale('log')
ax.set_xlabel('Number of cycles to failure, N', fontsize=13)
ax.set_ylabel('Stress amplitude, $\\sigma_a$ [MPa]', fontsize=13)
ax.set_title('S-N Curves: Steel vs Aluminum', fontsize=14)
ax.legend(fontsize=11)
ax.set_xlim(1e2, 1e9)
ax.set_ylim(50, 600)
ax.grid(True, which='both', alpha=0.3)
# 疲労限度領域をハイライト
ax.fill_between(N_steel, 0, fatigue_limit_steel, alpha=0.1, color='blue')
ax.annotate('Fatigue limit\n(infinite life)',
xy=(1e7, fatigue_limit_steel), xytext=(5e7, 320),
fontsize=11, arrowprops=dict(arrowstyle='->', color='blue'),
color='blue')
ax.annotate('No fatigue limit\n(curve keeps decreasing)',
xy=(1e8, stress_al[np.argmin(np.abs(N_al - 1e8))]),
xytext=(2e6, 120), fontsize=11,
arrowprops=dict(arrowstyle='->', color='red'), color='red')
plt.tight_layout()
plt.show()
この比較グラフからは、鉄鋼材料とアルミニウム合金の疲労特性の根本的な違いが一目瞭然です。
-
鉄鋼材料のS-N曲線は水平部をもつ: 約 $10^6$ 〜 $10^7$ 回を超えると曲線が水平になり、それ以下の応力では疲労破壊が起こりません。この水平部が疲労限度です。
-
アルミニウム合金のS-N曲線は水平部をもたない: 曲線は $10^9$ 回以上でも緩やかに下がり続けます。これは、どんなに小さい応力でも十分な繰り返し数が加われば破壊に至りうることを意味しており、航空機のアルミニウム構造の設計では有限寿命設計が必須である理由がここにあります。
-
低サイクル領域では鉄鋼の方が強い: $10^3$ 回以下の低サイクル領域では鉄鋼材料の方が高い応力振幅に耐えますが、高サイクル領域では両者の差が変化していくことが読み取れます。
この知識を踏まえて、次は平均応力の影響を修正グッドマン線図で可視化してみましょう。
Pythonによる修正グッドマン線図の描画
修正グッドマン線図、ゾーダーベルク線図、ガーバー線図の3つを同一平面上に描画し、それぞれの安全領域の違いを視覚的に比較します。
import numpy as np
import matplotlib.pyplot as plt
# 材料パラメータ（S45C相当の鋼材）
sigma_B = 570 # 引張強さ [MPa]
sigma_Y = 345 # 降伏応力 [MPa]
sigma_w = 245 # 疲労限度（R=-1）[MPa]
# 平均応力の範囲
sigma_m = np.linspace(0, sigma_B, 300)
# 修正グッドマン線: σ_a/σ_w + σ_m/σ_B = 1
goodman = sigma_w * (1 - sigma_m / sigma_B)
# ゾーダーベルク線: σ_a/σ_w + σ_m/σ_Y = 1
sigma_m_sod = np.linspace(0, sigma_Y, 300)
soderberg = sigma_w * (1 - sigma_m_sod / sigma_Y)
# ガーバー放物線: σ_a/σ_w + (σ_m/σ_B)^2 = 1
gerber = sigma_w * (1 - (sigma_m / sigma_B)**2)
# 描画
fig, ax = plt.subplots(figsize=(10, 7))
ax.plot(sigma_m, goodman, 'b-', linewidth=2.5, label='Modified Goodman')
ax.plot(sigma_m_sod, soderberg, 'g--', linewidth=2.5, label='Soderberg')
ax.plot(sigma_m, gerber, 'r-.', linewidth=2.5, label='Gerber')
# 降伏線: σ_a + σ_m = σ_Y
sigma_m_yield = np.linspace(0, sigma_Y, 200)
yield_line = sigma_Y - sigma_m_yield
ax.plot(sigma_m_yield, yield_line, 'k:', linewidth=1.5, label='Yield line')
# 安全領域の塗りつぶし（ゾーダーベルク内側）
ax.fill_between(sigma_m_sod, 0, soderberg, alpha=0.1, color='green',
label='Safe region (Soderberg)')
# 特性値の点をプロット
ax.plot(0, sigma_w, 'ko', markersize=8)
ax.annotate(f'$\\sigma_w$ = {sigma_w} MPa', xy=(0, sigma_w),
xytext=(30, sigma_w + 15), fontsize=11)
ax.plot(sigma_B, 0, 'ko', markersize=8)
ax.annotate(f'$\\sigma_B$ = {sigma_B} MPa', xy=(sigma_B, 0),
xytext=(sigma_B - 80, 20), fontsize=11)
ax.plot(sigma_Y, 0, 'ks', markersize=8)
ax.annotate(f'$\\sigma_Y$ = {sigma_Y} MPa', xy=(sigma_Y, 0),
xytext=(sigma_Y - 80, 30), fontsize=11)
# 設計点の例
sigma_m_design = 150
sigma_a_design = 120
ax.plot(sigma_m_design, sigma_a_design, 'r*', markersize=15, zorder=5)
ax.annotate(f'Design point\n({sigma_m_design}, {sigma_a_design}) MPa',
xy=(sigma_m_design, sigma_a_design),
xytext=(sigma_m_design + 40, sigma_a_design + 30),
fontsize=10, arrowprops=dict(arrowstyle='->'))
ax.set_xlabel('Mean stress, $\\sigma_m$ [MPa]', fontsize=13)
ax.set_ylabel('Stress amplitude, $\\sigma_a$ [MPa]', fontsize=13)
ax.set_title('Fatigue Diagram: Goodman, Soderberg, Gerber', fontsize=14)
ax.legend(fontsize=11, loc='upper right')
ax.set_xlim(-10, sigma_B + 30)
ax.set_ylim(-10, sigma_w + 60)
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
# 各基準での安全率を計算
def safety_factor_goodman(sa, sm, sw, sB):
return 1.0 / (sa / sw + sm / sB)
def safety_factor_soderberg(sa, sm, sw, sY):
return 1.0 / (sa / sw + sm / sY)
def safety_factor_gerber(sa, sm, sw, sB):
# ガーバーの安全率は二次方程式を解く必要がある
# n*σ_a/σ_w + (n*σ_m/σ_B)^2 = 1
a_coeff = (sm / sB)**2
b_coeff = sa / sw
if a_coeff == 0:
return sw / sa
discriminant = b_coeff**2 + 4 * a_coeff
return (-b_coeff + np.sqrt(discriminant)) / (2 * a_coeff) * (sa / sw + (sm/sB)**2) ** 0 # 簡易版
# 簡易的に Goodman と同様の計算
return 1.0 / (sa / sw + (sm / sB)**2)
n_goodman = safety_factor_goodman(sigma_a_design, sigma_m_design, sigma_w, sigma_B)
n_soderberg = safety_factor_soderberg(sigma_a_design, sigma_m_design, sigma_w, sigma_Y)
print(f"\n設計点: σ_a = {sigma_a_design} MPa, σ_m = {sigma_m_design} MPa")
print(f"グッドマン基準の安全率: n = {n_goodman:.2f}")
print(f"ゾーダーベルク基準の安全率: n = {n_soderberg:.2f}")
このグッドマン線図から、以下の重要な設計上の知見が読み取れます。
-
3つの基準の包含関係: ゾーダーベルク線（緑の破線）が最も内側にあり、次にグッドマン線（青の実線）、そしてガーバー放物線（赤の一点鎖線）が最も外側です。安全領域の大きさはゾーダーベルク < グッドマン < ガーバーの順であり、ゾーダーベルクが最も保守的な基準であることが視覚的に確認できます。
-
設計点の評価: 設計点（$\sigma_m = 150$ MPa, $\sigma_a = 120$ MPa）は3つの基準のいずれの内側にもあり、安全と判定されます。グッドマン基準での安全率がゾーダーベルク基準よりも大きいのは、グッドマンの方が許容範囲が広いためです。
-
降伏線との関係: 降伏線 $\sigma_a + \sigma_m = \sigma_Y$ よりも外側の領域では塑性変形が生じます。疲労限度の基準を満たしていても降伏線を超えている場合は、静的な過負荷による塑性変形のリスクがあります。
-
平均応力の影響の定量化: たとえばグッドマン線図から、平均応力が 0 MPa のとき許容される応力振幅は 245 MPa（疲労限度そのもの）ですが、平均応力が 285 MPa（引張強さの半分）になると許容振幅はその半分の約 123 MPa に低下します。平均応力の管理が疲労設計において極めて重要であることがわかります。
次に、マイナーの累積損傷則をPythonで実装し、変動振幅荷重下での疲労寿命を推定してみましょう。
Pythonによるマイナー則の累積損傷計算
一定の応力振幅がブロック状に繰り返されるケースを想定し、マイナーの累積損傷則による疲労寿命の推定を行います。
import numpy as np
import matplotlib.pyplot as plt
# --- S-N曲線のパラメータ（鋼材）---
sigma_f = 900 # 疲労強度係数 [MPa]
b = -0.085 # 疲労強度指数
sigma_w = 245 # 疲労限度 [MPa]
def sn_curve(sigma_a, sigma_f, b, sigma_w):
"""S-N曲線から破断寿命 N を計算する"""
if sigma_a <= sigma_w:
return np.inf # 疲労限度以下は無限寿命
# σ_a = σ_f' * (2N)^b → N = (1/2) * (σ_a / σ_f')^(1/b)
N = 0.5 * (sigma_a / sigma_f) ** (1.0 / b)
return N
# --- 荷重ブロックの定義 ---
# 各ブロック: (応力振幅 [MPa], 繰り返し数)
load_blocks = [
(400, 1000), # 高応力ブロック: 1,000回
(320, 10000), # 中応力ブロック: 10,000回
(270, 100000), # 低応力ブロック: 100,000回
(240, 500000), # 疲労限度以下: 500,000回（損傷寄与なし）
]
print("=" * 65)
print("マイナーの累積損傷則による疲労寿命推定")
print("=" * 65)
print(f"{'応力振幅':>10} {'繰り返し数':>12} {'破断寿命N':>14} {'損傷率D_i':>10}")
print("-" * 65)
D_total = 0.0
damages = []
for sigma_a, n in load_blocks:
N = sn_curve(sigma_a, sigma_f, b, sigma_w)
if N == np.inf:
D_i = 0.0
N_str = "∞ (below limit)"
else:
D_i = n / N
N_str = f"{N:.0f}"
D_total += D_i
damages.append(D_i)
print(f"{sigma_a:>10.0f} MPa {n:>12,} {N_str:>14} {D_i:>10.4f}")
print("-" * 65)
print(f"{'累積損傷 D':>38} = {D_total:>10.4f}")
print()
if D_total >= 1.0:
print("→ 破壊条件 D ≥ 1 を満たします。このブロックで破壊に至ります。")
else:
n_blocks = 1.0 / D_total
print(f"→ D < 1 のため、まだ破壊に至りません。")
print(f"→ 破壊までに同じ荷重ブロックを {n_blocks:.1f} 回繰り返せます。")
# 損傷の内訳を可視化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# 左: 損傷率の棒グラフ
labels = [f'{s} MPa\n({n:,} cycles)' for s, n in load_blocks]
colors = ['#e74c3c', '#e67e22', '#f1c40f', '#2ecc71']
ax1.bar(labels, damages, color=colors, edgecolor='black', linewidth=0.8)
ax1.set_ylabel('Damage fraction $D_i = n_i / N_i$', fontsize=12)
ax1.set_title('Damage Contribution per Load Block', fontsize=13)
ax1.grid(axis='y', alpha=0.3)
for i, d in enumerate(damages):
ax1.text(i, d + 0.005, f'{d:.4f}', ha='center', fontsize=10)
# 右: 累積損傷の推移
cumulative = np.cumsum(damages)
block_labels = [f'Block {i+1}' for i in range(len(load_blocks))]
ax2.step(range(len(cumulative)), cumulative, 'b-o', where='mid',
linewidth=2, markersize=8)
ax2.axhline(y=1.0, color='red', linestyle='--', linewidth=1.5,
label='Failure criterion (D=1)')
ax2.set_xticks(range(len(block_labels)))
ax2.set_xticklabels(labels, fontsize=9)
ax2.set_ylabel('Cumulative damage $D$', fontsize=12)
ax2.set_title('Cumulative Damage Progression', fontsize=13)
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_ylim(0, 1.2)
plt.tight_layout()
plt.show()
この結果から、マイナー則の動作とブロック荷重下の疲労損傷について、いくつかの重要な知見が得られます。
-
高応力ブロックの支配的寄与: 応力振幅 400 MPa のブロックは繰り返し数が最も少ない（1,000 回）にもかかわらず、破断寿命が短いため損傷率への寄与が非常に大きいことがわかります。疲労設計では、まれにしか起こらない大きな応力も無視できないことを意味しています。
-
疲労限度以下のブロック: 応力振幅 240 MPa（疲労限度 245 MPa 以下）のブロックは、繰り返し数が50万回と最大であるにもかかわらず、損傷率はゼロです。マイナー則では、疲労限度以下の応力は寿命に影響しないと仮定されます。ただし実際には、先行する高応力サイクルでき裂が発生している場合、疲労限度以下の応力でもき裂が進展する可能性があります。これはマイナー則の限界の一つです。
-
累積損傷の進展パターン: 右のグラフに示される累積損傷の推移から、高応力ブロックで急激に損傷が蓄積し、低応力ブロックでの損傷蓄積は緩やかであることが視覚的に確認できます。
次に、より現実的なランダム荷重下での疲労寿命推定に進みましょう。
Pythonによるランダム荷重下の疲労寿命推定
実際の構造物は、規則的なブロック荷重ではなく、不規則に変動する荷重を受けます。ランダム荷重から疲労寿命を推定するには、まず荷重履歴から応力サイクルを数え上げる必要があります。この「サイクル計数法」として最も広く使われているのがレインフロー法（rainflow counting method）です。
レインフロー法は、応力-時間履歴から個々の応力サイクル（応力振幅と平均応力の組）を抽出する手法です。名前の由来は、応力-時間履歴のグラフを横に倒して屋根と見なし、雨水が流れ落ちる様子に見立てたことに由来します。ここでは、Python の rainflow ライブラリを使わずに簡易的なレインフロー計数を実装し、マイナー則と組み合わせた疲労寿命推定を行います。
import numpy as np
import matplotlib.pyplot as plt
# --- ランダム荷重履歴の生成 ---
np.random.seed(42)
t = np.linspace(0, 10, 2000)
# 主要な正弦波荷重 + ランダムノイズ（模擬飛行荷重）
stress_history = (
180 # 定常平均応力 [MPa]
+ 80 * np.sin(2 * np.pi * 1.0 * t) # 主荷重成分
+ 40 * np.sin(2 * np.pi * 3.7 * t) # 高周波成分
+ 20 * np.random.randn(len(t)) # ランダムノイズ
)
# 描画: 応力-時間履歴
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(t, stress_history, 'b-', linewidth=0.5, alpha=0.8)
ax.set_xlabel('Time [s]', fontsize=12)
ax.set_ylabel('Stress [MPa]', fontsize=12)
ax.set_title('Random Stress History (Simulated Flight Load)', fontsize=13)
ax.axhline(y=np.mean(stress_history), color='red', linestyle='--',
alpha=0.7, label=f'Mean = {np.mean(stress_history):.1f} MPa')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"応力履歴の統計:")
print(f" 平均応力: {np.mean(stress_history):.1f} MPa")
print(f" 最大応力: {np.max(stress_history):.1f} MPa")
print(f" 最小応力: {np.min(stress_history):.1f} MPa")
print(f" 標準偏差: {np.std(stress_history):.1f} MPa")
上のグラフは、航空機の翼に作用する荷重を模した模擬的な応力-時間履歴です。主荷重成分（1 Hz の正弦波）に高周波成分とランダムノイズが重畳されており、現実の荷重に近い不規則な変動パターンを示しています。平均応力は約 180 MPa であり、ゼロではない定常的な引張応力が存在します。
次に、この荷重履歴から応力サイクルを計数し、マイナー則で累積損傷を計算する簡易レインフロー法を実装します。
import numpy as np
import matplotlib.pyplot as plt
def find_peaks_and_valleys(signal):
"""時系列データからピークと谷を抽出する"""
extrema = []
extrema_idx = []
for i in range(1, len(signal) - 1):
if signal[i] > signal[i-1] and signal[i] > signal[i+1]:
extrema.append(signal[i])
extrema_idx.append(i)
elif signal[i] < signal[i-1] and signal[i] < signal[i+1]:
extrema.append(signal[i])
extrema_idx.append(i)
return np.array(extrema), np.array(extrema_idx)
def simple_rainflow(extrema):
"""簡易レインフロー法（3点法）によるサイクル計数"""
cycles = [] # (応力範囲, 平均応力, 回数)
residue = list(extrema)
i = 0
while i < len(residue) - 2:
s0, s1, s2 = residue[i], residue[i+1], residue[i+2]
r1 = abs(s1 - s0) # 前の振幅
r2 = abs(s2 - s1) # 後の振幅
if r2 >= r1:
# 閉じたサイクルを抽出
stress_range = r1
mean_stress = (s0 + s1) / 2.0
cycles.append((stress_range, mean_stress, 0.5))
residue.pop(i + 1)
residue.pop(i)
i = max(0, i - 1)
else:
i += 1
# 残余をハーフサイクルとして処理
for j in range(len(residue) - 1):
sr = abs(residue[j+1] - residue[j])
sm = (residue[j] + residue[j+1]) / 2.0
cycles.append((sr, sm, 0.5))
return cycles
# --- 前のセクションの応力履歴を再生成 ---
np.random.seed(42)
t = np.linspace(0, 10, 2000)
stress_history = (
180 + 80 * np.sin(2 * np.pi * 1.0 * t)
+ 40 * np.sin(2 * np.pi * 3.7 * t)
+ 20 * np.random.randn(len(t))
)
# ピーク・谷の抽出とレインフロー計数
extrema, _ = find_peaks_and_valleys(stress_history)
cycles = simple_rainflow(extrema)
# サイクルデータの集計
stress_ranges = np.array([c[0] for c in cycles])
mean_stresses = np.array([c[1] for c in cycles])
stress_amplitudes = stress_ranges / 2.0
counts = np.array([c[2] for c in cycles])
print(f"抽出されたサイクル数: {len(cycles)}")
print(f"応力振幅の範囲: {stress_amplitudes.min():.1f} ~ "
f"{stress_amplitudes.max():.1f} MPa")
print(f"平均応力の範囲: {mean_stresses.min():.1f} ~ "
f"{mean_stresses.max():.1f} MPa")
# --- マイナー則による累積損傷計算 ---
sigma_f_prime = 900 # 疲労強度係数 [MPa]
b_exp = -0.085 # 疲労強度指数
sigma_w_base = 245 # 基本疲労限度 [MPa]（R=-1）
sigma_B = 570 # 引張強さ [MPa]
def equivalent_amplitude_goodman(sigma_a, sigma_m, sigma_w, sigma_B):
"""グッドマン線図で等価完全両振り振幅に変換"""
if sigma_m >= sigma_B:
return np.inf
return sigma_a / (1 - sigma_m / sigma_B)
def fatigue_life(sigma_a_eq, sigma_f, b, sigma_w):
"""等価応力振幅からS-N曲線で破断寿命を計算"""
if sigma_a_eq <= sigma_w:
return np.inf
return 0.5 * (sigma_a_eq / sigma_f) ** (1.0 / b)
D_total = 0.0
n_damaging = 0
for sa, sm, cnt in zip(stress_amplitudes, mean_stresses, counts):
# グッドマン補正で等価完全両振り振幅に変換
sa_eq = equivalent_amplitude_goodman(sa, sm, sigma_w_base, sigma_B)
N_f = fatigue_life(sa_eq, sigma_f_prime, b_exp, sigma_w_base)
if N_f < np.inf:
D_total += cnt / N_f
n_damaging += 1
print(f"\n累積損傷計算結果:")
print(f" 損傷に寄与するサイクル数: {n_damaging} / {len(cycles)}")
print(f" この荷重履歴1ブロックの累積損傷: D = {D_total:.6f}")
if D_total > 0:
blocks_to_failure = 1.0 / D_total
print(f" 破壊までのブロック数: {blocks_to_failure:.0f} ブロック")
total_time = blocks_to_failure * t[-1]
print(f" 推定疲労寿命: {total_time:.0f} 秒 "
f"= {total_time/3600:.1f} 時間")
# レインフロー結果の可視化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# 左: 応力振幅のヒストグラム
ax1.hist(stress_amplitudes, bins=30, color='steelblue',
edgecolor='black', alpha=0.8)
ax1.set_xlabel('Stress amplitude [MPa]', fontsize=12)
ax1.set_ylabel('Number of cycles', fontsize=12)
ax1.set_title('Rainflow Cycle Counting: Amplitude Distribution',
fontsize=13)
ax1.axvline(x=sigma_w_base, color='red', linestyle='--',
linewidth=1.5, label=f'Fatigue limit = {sigma_w_base} MPa')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
# 右: 散布図（応力振幅 vs 平均応力）+グッドマン線図
sm_line = np.linspace(0, sigma_B, 200)
goodman_line = sigma_w_base * (1 - sm_line / sigma_B)
ax2.scatter(mean_stresses, stress_amplitudes, c='steelblue',
s=15, alpha=0.5, label='Rainflow cycles')
ax2.plot(sm_line, goodman_line, 'r-', linewidth=2,
label='Modified Goodman')
ax2.set_xlabel('Mean stress $\\sigma_m$ [MPa]', fontsize=12)
ax2.set_ylabel('Stress amplitude $\\sigma_a$ [MPa]', fontsize=12)
ax2.set_title('Rainflow Cycles on Goodman Diagram', fontsize=13)
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, sigma_B + 20)
ax2.set_ylim(0, sigma_w_base + 50)
plt.tight_layout()
plt.show()
この解析結果から、ランダム荷重下での疲労評価に関する重要な知見が得られます。
-
応力振幅分布の特徴: 左のヒストグラムから、レインフロー法で抽出されたサイクルの応力振幅は小さいサイクルが圧倒的に多く、大きな応力振幅のサイクルはごくわずかであることがわかります。このような分布は実際のランダム荷重でも一般的に見られるもので、少数の大振幅サイクルが損傷の大部分を支配します。
-
グッドマン線図上のサイクル分布: 右の散布図では、レインフロー法で得られた個々のサイクルがグッドマン線図上にプロットされています。ほとんどのサイクルがグッドマン線の内側（安全領域）に位置していることが確認できます。しかし、安全領域内であっても各サイクルがわずかずつ損傷を蓄積し、それが何万ブロックにもわたって累積することで最終的に破壊に至ります。
-
寿命推定の実用性: 1ブロック（10秒分の荷重履歴）あたりの累積損傷は非常に小さい値ですが、マイナー則により $D = 1$ に達するまでの繰り返しブロック数を逆算することで、推定疲労寿命が得られます。この推定値は設計の初期段階での目安として有用ですが、前述したマイナー則の限界（荷重順序効果の無視など）を考慮した安全率の設定が実務では不可欠です。
最後に、応力集中の影響を含めた総合的な疲労評価のまとめとして、切欠き係数を考慮したS-N曲線の修正を可視化します。
Pythonによる切欠き係数を考慮したS-N曲線
応力集中がある部品では、S-N曲線を切欠き係数 $K_f$ で修正する必要があります。ここでは、平滑材と切欠き材のS-N曲線を比較し、表面仕上げと寸法効果も含めた修正S-N曲線を描画します。
import numpy as np
import matplotlib.pyplot as plt
# 材料パラメータ（S45C鋼）
sigma_B = 570 # 引張強さ [MPa]
sigma_w = 245 # 平滑試験片の疲労限度 [MPa]（R=-1）
sigma_f = 900 # 疲労強度係数 [MPa]
b = -0.085 # 疲労強度指数
# 修正係数
Kt = 2.5 # 応力集中係数（円孔つき平板）
q = 0.85 # 切欠き感度係数
Kf = 1 + q * (Kt - 1) # 疲労切欠き係数
C_surf = 0.75 # 表面仕上げ係数（機械加工仕上げ）
C_size = 0.85 # 寸法係数（直径 50mm の軸）
# 修正疲労限度
sigma_w_mod = (C_surf * C_size / Kf) * sigma_w
print("疲労強度の修正:")
print(f" 応力集中係数 Kt = {Kt}")
print(f" 切欠き感度係数 q = {q}")
print(f" 疲労切欠き係数 Kf = {Kf:.2f}")
print(f" 表面仕上げ係数 C_surf = {C_surf}")
print(f" 寸法係数 C_size = {C_size}")
print(f" 平滑試験片の疲労限度: σ_w = {sigma_w} MPa")
print(f" 修正後の疲労限度: σ_w,mod = {sigma_w_mod:.1f} MPa")
print(f" 低減率: {sigma_w_mod/sigma_w*100:.1f}%")
# S-N曲線の計算
N = np.logspace(2, 8, 300)
# 平滑材のS-N曲線
stress_smooth = sigma_f * (2 * N)**b
stress_smooth = np.maximum(stress_smooth, sigma_w)
# 修正S-N曲線（切欠き + 表面 + 寸法効果）
# 高サイクル側の修正: 疲労限度を低減
# 低サイクル側(N=1000)では修正の影響は小さいと仮定
sigma_1000_smooth = sigma_f * (2 * 1000)**b # N=1000での応力
sigma_1000_mod = sigma_1000_smooth / Kf # 低サイクル側の修正
# 対数軸で線形補間
log_N = np.log10(N)
log_N_low = np.log10(1000)
log_N_fatigue = 7 # 10^7 回で疲労限度に達する
stress_mod = np.zeros_like(N)
for i, ln in enumerate(log_N):
if ln <= log_N_low:
stress_mod[i] = sigma_1000_mod
elif ln >= log_N_fatigue:
stress_mod[i] = sigma_w_mod
else:
# 対数軸で線形補間
frac = (ln - log_N_low) / (log_N_fatigue - log_N_low)
log_s = np.log10(sigma_1000_mod) * (1 - frac) + \
np.log10(sigma_w_mod) * frac
stress_mod[i] = 10**log_s
# 描画
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(N, stress_smooth, 'b-', linewidth=2.5,
label='Smooth specimen (standard)')
ax.plot(N, stress_mod, 'r--', linewidth=2.5,
label=f'Notched part ($K_f$={Kf:.2f}, $C_{{surf}}$={C_surf}, '
f'$C_{{size}}$={C_size})')
ax.axhline(y=sigma_w, color='blue', linestyle=':', alpha=0.5,
label=f'$\\sigma_w$ (smooth) = {sigma_w} MPa')
ax.axhline(y=sigma_w_mod, color='red', linestyle=':', alpha=0.5,
label=f'$\\sigma_{{w,mod}}$ = {sigma_w_mod:.0f} MPa')
# 修正の効果を矢印で示す
N_show = 1e5
s_smooth_show = sigma_f * (2 * N_show)**b
idx_mod = np.argmin(np.abs(N - N_show))
s_mod_show = stress_mod[idx_mod]
ax.annotate('', xy=(N_show, s_mod_show), xytext=(N_show, s_smooth_show),
arrowprops=dict(arrowstyle='<->', color='green', lw=2))
ax.text(N_show * 1.5, (s_smooth_show + s_mod_show) / 2,
f'Reduction\n~{(1 - s_mod_show/s_smooth_show)*100:.0f}%',
fontsize=11, color='green')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('Number of cycles to failure, N', fontsize=13)
ax.set_ylabel('Stress amplitude, $\\sigma_a$ [MPa]', fontsize=13)
ax.set_title('S-N Curve: Smooth vs Notched Part', fontsize=14)
ax.legend(fontsize=10, loc='lower left')
ax.set_xlim(1e2, 1e8)
ax.set_ylim(50, 800)
ax.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.show()
この比較グラフは、実際の部品設計において考慮すべき疲労強度の低減を明確に示しています。
-
大幅な疲労限度の低下: 平滑試験片の疲労限度 245 MPa に対して、切欠き・表面仕上げ・寸法効果を考慮した修正疲労限度は大幅に低下しています。標準試験片のデータをそのまま使うと、実部品の疲労強度を大きく過大評価してしまう危険性が視覚的にわかります。
-
低サイクル側と高サイクル側の影響の違い: 高サイクル領域では修正の影響が大きく、低サイクル領域では影響が相対的に小さくなっています。これは、高サイクル疲労が弾性域のき裂発生で支配されるため応力集中に敏感であるのに対し、低サイクル疲労では塑性変形が先行するため応力集中の先端が鈍化し、影響が緩和されることに対応しています。
-
設計上の教訓: 疲労設計では、材料データベースの値をそのまま使うのではなく、実際の部品形状・表面仕上げ・寸法に応じた修正を行うことが不可欠です。特に応力集中部の疲労限度は、平滑試験片の値の半分以下になることも珍しくありません。フィレットの曲率半径を大きくする、表面をショットピーニングで処理するなどの対策が重要です。
まとめ
本記事では、疲労破壊の基礎理論からPythonによる実装まで、以下の内容を解説しました。
- 疲労破壊の本質: 静的強度よりはるかに低い応力でも、繰り返し荷重によって破壊に至る。疲労破壊は表面での微小き裂の発生→安定進展→最終破断の3段階で進行し、破壊直前まで目に見える変形がほとんどない
- S-N曲線: 応力振幅と破断繰り返し数の関係を表す基本ツール。両対数軸上ではバスキンの式 $\sigma_a = \sigma_f'(2N)^b$ で近似的に直線になる
- 疲労限度: 鉄鋼材料では $10^6$ 〜 $10^7$ 回付近に疲労限度が存在し、それ以下の応力では無限寿命が得られる。一方、アルミニウム合金では明確な疲労限度が存在しない
- 平均応力の影響: 引張側の平均応力は疲労寿命を低下させる。修正グッドマン線図 $\sigma_a/\sigma_w + \sigma_m/\sigma_B = 1$ やゾーダーベルク線図で安全領域を判定する
- マイナーの累積損傷則: 変動振幅荷重下での疲労寿命を $D = \sum n_i/N_i = 1$ で推定する。計算は簡便だが、荷重順序効果を無視する等の限界がある
- 応力集中と修正係数: 実部品の疲労限度は、切欠き係数 $K_f$、表面仕上げ、寸法効果によって平滑試験片の値から大幅に低下する
疲労破壊の理解は、材料力学・構造設計の基盤となる重要な知識です。本記事で扱った線形累積損傷則（マイナー則）は「第一近似」としては有用ですが、より精密な寿命推定には破壊力学に基づくき裂進展解析が必要になります。
次のステップとして、以下の記事も参考にしてください。
- 破壊力学入門 — 応力拡大係数とき裂進展速度（Paris則）による疲労寿命予測
- クリープと応力緩和 — 高温環境での時間依存型破壊。疲労とクリープの相互作用も重要な設計課題
- 組合せ応力 — 多軸応力下での疲労評価に必要な応力状態の理解
- 主応力とモールの応力円 — 多軸疲労における等価応力の算出に必要な知識
- 材料試験法 — 疲労試験を含む各種材料試験の手法と規格