ジェットエンジンのタービンブレードは、1000°C を超える高温ガスの中で何万時間も回り続けます。室温では問題ない応力であっても、高温で長時間かかり続けると、材料はじわじわと変形し、最終的には破断に至ります。「一定の荷重なのに変形が進む」——この現象がクリープ（creep）です。
一方、高温のボルト締結では逆の問題が生じます。締め付けたボルトのひずみ（伸び）は変わらないのに、時間とともにボルトの軸力が低下してしまいます。「一定のひずみなのに応力が減る」——これが応力緩和（stress relaxation）です。
クリープと応力緩和は、時間の経過とともに材料の力学的な状態が変わる時間依存変形の二つの代表的な現象です。この二つを理解すると、以下のような幅広い分野で設計や評価に活かすことができます。
- 航空宇宙: ジェットエンジンのタービンブレード、ロケットノズルなど、極限温度環境で長寿命を保証する設計に不可欠です
- 発電プラント: 火力発電のボイラー管や蒸気タービンは数十年にわたりクリープにさらされ、定期的な余寿命評価が求められます
- 半導体実装: はんだ接合部は比較的低温でも融点比（ホモロガス温度）が高いため、クリープが顕著に発生します
- 高分子・ゴム材料: 自動車のシール材やゴムパッキンは常温でも粘弾性的な挙動を示し、応力緩和による密封力の低下が問題となります
本記事の内容
- クリープとは何か——時間依存変形の直感的な理解
- クリープ曲線の3段階（遷移クリープ・定常クリープ・加速クリープ）
- クリープの微視的メカニズム（拡散クリープと転位クリープ）
- 応力緩和の概念と数学的定義
- 粘弾性モデル：マクスウェルモデル、フォークトモデル、標準線形固体モデル
- ラーソン・ミラーパラメータによる高温寿命推定
- 高温設計への応用
- Python による実装：クリープ曲線の描画、粘弾性モデルの応答シミュレーション、応力緩和、寿命推定
前提知識
この記事を読む前に、以下の記事に目を通しておくと理解がスムーズです。
- 応力とひずみの定義 — 応力・ひずみの基本概念と単位
- フックの法則と弾性係数 — 弾性域における応力-ひずみの線形関係
- ヤング率とポアソン比 — 弾性定数の物理的意味
- 微分方程式入門 — 粘弾性モデルの構成方程式で1階常微分方程式を解く
クリープとは——「時間が敵になる」変形
日常の直感から始めよう
クリープを身近に体感するには、スーパーの買い物袋を思い浮かべてください。重い荷物を入れたビニール袋を手に持ったまましばらく立っていると、袋の持ち手が少しずつ伸びてきます。荷物の重さ（力）は変わっていないのに、変形がじわじわと進みます。あるいは、プラスチックの物干し竿に洗濯物をかけっぱなしにしておくと、数年で中央が垂れ下がってきます。これも同じ現象です。
金属材料でも同様のことが起こりますが、室温ではほぼ無視できるほど遅く、温度が高くなると急激に顕著になります。経験的な目安として、材料の絶対温度が融点の約 0.3〜0.4 倍を超えると、クリープが工学的に無視できなくなります。この比率をホモロガス温度 $T_H = T / T_m$ と呼びます（$T$ は使用温度、$T_m$ は融点、ともに絶対温度）。
たとえば、鉄の融点は約 1811 K（1538°C）なので、$0.4 \times 1811 \approx 724$ K、つまり約 450°C 以上でクリープが問題になり始めます。一方、鉛の融点は 600 K（327°C）と低いため、室温（300 K）でも $T_H = 0.5$ であり、室温でクリープが進行します。はんだ（Sn-Pb系で融点約 456 K）も同様の理由で、電子部品の接合部でクリープが重要な設計因子となります。
クリープの数学的定義
クリープ試験は、材料試験片に一定の引張荷重を加え、一定温度の環境下でひずみの時間変化を測定する試験です。応力 $\sigma_0$ と温度 $T$ を一定に保ったとき、ひずみ $\varepsilon$ が時間 $t$ の関数として増加する現象をクリープと定義します。
$$ \begin{equation} \varepsilon = \varepsilon(t) \quad \text{at} \quad \sigma = \sigma_0 = \text{const.}, \quad T = \text{const.} \end{equation} $$
ここで重要なのは、フックの法則が予測する弾性ひずみ $\varepsilon_e = \sigma_0 / E$ は時間に依存しません。クリープひずみとは、この瞬間的な弾性ひずみの上に、時間とともに蓄積される追加のひずみです。全ひずみは次のように分解されます。
$$ \begin{equation} \varepsilon(t) = \varepsilon_e + \varepsilon_c(t) = \frac{\sigma_0}{E} + \varepsilon_c(t) \end{equation} $$
$\varepsilon_e$ は荷重を加えた瞬間に生じる弾性ひずみ、$\varepsilon_c(t)$ が時間依存のクリープひずみです。クリープひずみの時間変化率 $\dot{\varepsilon}_c = d\varepsilon_c / dt$ をクリープ速度（creep rate）と呼びます。
クリープの全体像をつかんだところで、次にクリープ曲線の典型的な形を詳しく見ていきましょう。
クリープ曲線の3段階
典型的なクリープ曲線の形
一定温度・一定応力のクリープ試験で得られるひずみ-時間曲線（クリープ曲線）は、典型的に3つの段階に分けられます。横軸が時間、縦軸がひずみのグラフを思い浮かべてください。
-
第I期：遷移クリープ（transient creep / primary creep） — 荷重を加えた直後、クリープ速度は比較的大きいものの、時間とともに減少していく領域です。材料内部で加工硬化（転位の蓄積や絡み合い）が進むため、変形が徐々に抑制されます。この段階のひずみは、時間のべき乗関数で近似できることが多く、$\varepsilon_c \propto t^{1/3}$ というAndrade のクリープ則が古くから知られています。
-
第II期：定常クリープ（steady-state creep / secondary creep） — クリープ速度がほぼ一定になる領域です。加工硬化と回復（転位の消滅や再配列）がバランスし、動的平衡状態に達しています。この一定のクリープ速度を最小クリープ速度 $\dot{\varepsilon}_s$ と呼びます。工学設計では、この $\dot{\varepsilon}_s$ が最も重要なパラメータです。なぜなら、部材の寿命の大部分は定常クリープ段階で費やされるためです。
-
第III期：加速クリープ（accelerating creep / tertiary creep） — クリープ速度が再び増加し、最終的にクリープ破断（rupture）に至る領域です。ネッキング（断面積の局所的な減少）や内部での微小空孔・亀裂の成長により、実効応力が増加してクリープが加速します。疲労破壊と同様に、最終段階での損傷の蓄積は急速に進行します。
クリープ曲線の経験式
クリープ曲線を数学的に表現するために、いくつかの経験式が提案されています。代表的なものとして、3段階を1つの式で表す次の形があります。
$$ \begin{equation} \varepsilon_c(t) = A t^{1/3} + B t + C \left( e^{\gamma t} – 1 \right) \end{equation} $$
右辺の各項はそれぞれの段階に対応しています。
第1項 $A t^{1/3}$ は遷移クリープを表します。$t$ のべき乗（$1/3$ 乗）であるため、初期にはひずみが急速に増加しますが、時間とともに増加率が鈍化します。
第2項 $Bt$ は定常クリープを表します。$B$ が最小クリープ速度 $\dot{\varepsilon}_s$ に相当し、時間に対して線形に増加します。
第3項 $C(e^{\gamma t} – 1)$ は加速クリープを表します。指数関数的に増加するため、時間の後半でひずみが急激に増大します。$\gamma$ は加速の度合いを決めるパラメータです。
また、定常クリープ速度 $\dot{\varepsilon}_s$ の応力依存性は、Norton の冪乗則（power-law creep）で記述されることが多いです。
$$ \begin{equation} \dot{\varepsilon}_s = A’ \sigma^n \exp\left(-\frac{Q}{RT}\right) \end{equation} $$
ここで $A’$ は材料定数、$n$ は応力指数（通常 3〜8 程度）、$Q$ はクリープの活性化エネルギー [J/mol]、$R$ は気体定数 $8.314 \, \text{J/(mol·K)}$、$T$ は絶対温度 [K] です。この式は、クリープ速度が応力の冪乗に比例し、温度に対してアレニウス型（指数関数型）の依存性を持つことを表しています。
クリープ曲線の定量的な描画は後ほどPythonで実装します。その前に、クリープがなぜ起こるのか、その微視的なメカニズムを掘り下げてみましょう。
クリープのメカニズム——原子スケールで何が起きているのか
クリープは巨視的には「ゆっくり変形する」という単純な現象ですが、その背後では複雑な微視的過程が進行しています。温度と応力の条件に応じて支配的なメカニズムが異なり、大きく分けて拡散クリープと転位クリープに分類されます。
拡散クリープ（低応力・高温）
高温では、原子は格子点に留まっているだけでなく、熱エネルギーを得て格子中を移動します。この原子の拡散によって巨視的な変形が生じるのが拡散クリープです。拡散（フィックの法則）で学んだように、拡散束は濃度勾配に比例しますが、クリープにおいては応力勾配が拡散の駆動力となります。
拡散クリープには2つの主要な経路があります。
ナバロ・ヘリングクリープ（Nabarro-Herring creep） は、原子が結晶粒の内部（格子内）を拡散する機構です。引張方向に垂直な粒界から原子が離れ、引張方向に平行な粒界に向かって拡散します。これは原子の正味の移動が結晶粒を引張方向に伸ばす効果を生みます。クリープ速度は次のように表されます。
$$ \begin{equation} \dot{\varepsilon}_{\text{NH}} = \frac{C_1 D_l \sigma \Omega}{k T d^2} \end{equation} $$
ここで $D_l$ は格子拡散係数、$\Omega$ は原子体積、$k$ はボルツマン定数、$d$ は結晶粒径、$C_1$ は幾何学的定数です。重要なのは、クリープ速度が $\sigma$ の1乗に比例し（$n = 1$）、結晶粒径 $d$ の2乗に反比例することです。つまり、結晶粒を大きくすると拡散クリープは抑制されます。
コブルクリープ（Coble creep） は、原子が結晶粒界に沿って拡散する機構です。粒界は格子内よりも原子が移動しやすい「高速道路」のような経路であり、特に低温側（拡散クリープの範囲内で相対的に低温）では粒界拡散が支配的になります。
$$ \begin{equation} \dot{\varepsilon}_{\text{Co}} = \frac{C_2 D_{gb} \delta_{gb} \sigma \Omega}{k T d^3} \end{equation} $$
ここで $D_{gb}$ は粒界拡散係数、$\delta_{gb}$ は粒界の有効幅です。結晶粒径 $d$ の3乗に反比例するため、ナバロ・ヘリングクリープよりもさらに粒径の影響が大きくなります。微細粒材料ではコブルクリープが支配的になることがあるため、高温用途の材料では粒径の制御が重要な設計指針となります。
転位クリープ（高応力）
比較的高い応力下では、結晶格子中の線欠陥である転位（dislocation）の運動がクリープの主役となります。室温でも塑性変形は転位のすべり運動によって生じますが、高温では転位が障害物を上昇運動（climb）によって乗り越えることができるようになります。
上昇運動とは、転位が滑り面から垂直方向に移動する過程で、空孔の拡散を伴います。室温では空孔の拡散が極めて遅いため上昇運動は起こりませんが、高温ではこれが可能になり、転位が障害物を迂回できるようになります。
転位クリープにおける定常クリープ速度は、先に示した Norton の冪乗則で記述されます。
$$ \begin{equation} \dot{\varepsilon}_s = A’ \sigma^n \exp\left(-\frac{Q}{RT}\right) \end{equation} $$
転位クリープでは応力指数 $n$ は典型的に 3〜8 程度であり、拡散クリープ（$n = 1$）よりも応力依存性が強いことが特徴です。これは直感的に理解できます。応力が高いほど転位の駆動力が大きくなり、転位密度も増加するため、クリープ速度は応力のべき乗で加速します。
変形機構図（デフォメーションマップ）
上で見た複数のクリープ機構は、温度と応力の組み合わせによって支配的なものが切り替わります。この関係を視覚的に整理したものが変形機構図（deformation mechanism map、Ashby map）です。横軸にホモロガス温度 $T/T_m$、縦軸に正規化応力 $\sigma/G$（$G$ はせん断弾性率）を取り、各領域にどの変形機構が支配的かを示します。
| 領域 | 温度 | 応力 | 支配メカニズム |
|---|---|---|---|
| 低温・高応力 | $T/T_m < 0.3$ | 高 | 転位すべり（通常の塑性変形） |
| 中温・高応力 | $0.3 < T/T_m < 0.7$ | 高 | 転位クリープ（上昇運動+すべり） |
| 高温・低応力 | $T/T_m > 0.5$ | 低 | 拡散クリープ（NH / Coble） |
| 極低応力 | 任意 | 極低 | 弾性変形のみ |
この変形機構図は、合金設計や使用環境の評価において非常に有用なツールです。例えば、ジェットエンジンのタービンブレードに用いるニッケル基超合金は、転位クリープ領域で使用されることが多く、$\gamma’$ 析出物による転位のピン止め効果がクリープ抑制に効果的です。
クリープの微視的メカニズムを理解したところで、次は「一定荷重ではなく一定変形を維持したとき」に何が起こるか、すなわち応力緩和の世界に踏み込んでみましょう。
応力緩和——「力がゆるむ」現象
直感的な理解
応力緩和を体感するために、輪ゴムの実験を想像してください。輪ゴムをある長さまで引き伸ばし、その状態で両端を固定します。最初は輪ゴムは強い張力で「元に戻りたい」と引っ張り返していますが、しばらくそのまま放置すると、引っ張る力（張力）が徐々に弱くなっていくのが感じられます。長さ（ひずみ）は固定しているのに、力（応力）が減少する——これが応力緩和です。
金属材料でもまったく同様のことが起きます。ボイラーの高温フランジを締結するボルトは、一定の伸び（ひずみ）を与えて締め付けますが、高温で長期間使用すると軸力が低下します。これが「ボルトの応力緩和」であり、定期的な増し締めが必要となる理由です。
応力緩和の数学的定義
応力緩和試験は、一定ひずみ $\varepsilon_0$ を維持した状態で、応力 $\sigma$ の時間変化を測定する試験です。
$$ \begin{equation} \sigma = \sigma(t) \quad \text{at} \quad \varepsilon = \varepsilon_0 = \text{const.}, \quad T = \text{const.} \end{equation} $$
クリープとの対比を明確にしておきましょう。
| 現象 | 一定に保つもの | 時間変化するもの |
|---|---|---|
| クリープ | 応力 $\sigma_0$（および温度 $T$） | ひずみ $\varepsilon(t)$：増加する |
| 応力緩和 | ひずみ $\varepsilon_0$（および温度 $T$） | 応力 $\sigma(t)$：減少する |
この二つは同じコインの裏表のような関係にあります。どちらも材料の時間依存変形（粘性的な性質）に起因しますが、境界条件（力一定 vs. 変位一定）が異なるだけです。
応力緩和の物理的メカニズムを考えてみましょう。全ひずみ $\varepsilon_0$ は一定ですが、この全ひずみは弾性ひずみ $\varepsilon_e$ とクリープひずみ $\varepsilon_c$ の和で構成されています。
$$ \begin{equation} \varepsilon_0 = \varepsilon_e(t) + \varepsilon_c(t) = \frac{\sigma(t)}{E} + \varepsilon_c(t) = \text{const.} \end{equation} $$
時間が経過すると、クリープひずみ $\varepsilon_c$ が増加します。全ひずみが一定であるためには、弾性ひずみ $\varepsilon_e = \sigma / E$ が対応する分だけ減少しなければなりません。ヤング率とポアソン比で学んだように、弾性ひずみは応力に比例するため、弾性ひずみの減少は応力の減少を意味します。
上式の両辺を時間で微分すると、$\varepsilon_0$ は定数なので左辺はゼロになります。
$$ \begin{equation} 0 = \frac{1}{E}\frac{d\sigma}{dt} + \dot{\varepsilon}_c \end{equation} $$
この式を整理すると、応力の時間変化率が得られます。
$$ \begin{equation} \frac{d\sigma}{dt} = -E \dot{\varepsilon}_c \end{equation} $$
この式は、クリープ速度 $\dot{\varepsilon}_c$ が正である限り、応力は必ず減少することを示しています。マイナス符号が応力の「緩和」を表現しています。クリープ速度が応力に依存する場合（Norton則など）、この微分方程式を数値的に解くことで応力緩和曲線を予測できます。
応力緩和の数学的構造がわかったところで、次はクリープと応力緩和を統一的に記述する枠組み——粘弾性モデルの世界へ進みましょう。
粘弾性モデル——バネとダッシュポットで時間依存挙動を表現する
なぜ粘弾性モデルが必要か
フックの法則で記述される弾性体は、力を加えると瞬時に変形し、力を除くと完全に元に戻ります。一方、蜂蜜やアスファルトのような粘性体（ニュートン流体）は、力を加えると一定の速度で変形し続け、力を除いても変形は元に戻りません。
しかし、現実の多くの材料は、弾性と粘性の両方の性質を併せ持っています。高分子材料は典型的な粘弾性体ですし、高温の金属もクリープ挙動を通じて粘弾性的に振る舞います。このような時間依存的な変形挙動を記述するために、バネ（弾性要素）とダッシュポット（粘性要素）を組み合わせた力学モデルが用いられます。
バネはフックの法則に従い、応力とひずみが比例します。
$$ \begin{equation} \sigma = E \varepsilon \quad \text{（バネ：弾性要素）} \end{equation} $$
ダッシュポットは粘性流体のように、応力がひずみ速度に比例します。
$$ \begin{equation} \sigma = \eta \dot{\varepsilon} \quad \text{（ダッシュポット：粘性要素）} \end{equation} $$
ここで $\eta$ は粘性係数（viscosity）[Pa·s] です。
これら2つの要素の接続のしかた（直列・並列）を変えることで、異なるタイプの粘弾性挙動を表現できます。以下、3つの代表的なモデルを見ていきましょう。
マクスウェルモデル（Maxwell model）——応力緩和を記述する
マクスウェルモデルは、バネとダッシュポットを直列に接続したモデルです。直列接続のポイントは、「両方に同じ応力がかかり、全体のひずみは各要素のひずみの和」であることです。
直列接続における条件を整理しましょう。
$$ \sigma_{\text{spring}} = \sigma_{\text{dashpot}} = \sigma \quad \text{（応力は共通）} $$
$$ \varepsilon = \varepsilon_{\text{spring}} + \varepsilon_{\text{dashpot}} \quad \text{（ひずみは加算）} $$
両辺を時間で微分すると、全ひずみ速度が得られます。
$$ \begin{equation} \dot{\varepsilon} = \dot{\varepsilon}_{\text{spring}} + \dot{\varepsilon}_{\text{dashpot}} = \frac{\dot{\sigma}}{E} + \frac{\sigma}{\eta} \end{equation} $$
これがマクスウェルモデルの構成方程式（微分形式）です。微分方程式入門で学んだ1階線形常微分方程式の形をしています。
応力緩和（$\varepsilon = \varepsilon_0 = \text{const.}$）の解:
$\dot{\varepsilon} = 0$ とすると、
$$ 0 = \frac{\dot{\sigma}}{E} + \frac{\sigma}{\eta} $$
整理すると、
$$ \frac{d\sigma}{dt} = -\frac{E}{\eta} \sigma $$
これは変数分離型の微分方程式です。$\tau_R = \eta / E$ と定義すると（$\tau_R$ は緩和時間と呼ばれます）、
$$ \frac{d\sigma}{\sigma} = -\frac{dt}{\tau_R} $$
両辺を積分すると、
$$ \ln \sigma = -\frac{t}{\tau_R} + C $$
初期条件 $\sigma(0) = \sigma_0 = E \varepsilon_0$ を代入して定数 $C$ を決めると、応力緩和の解が得られます。
$$ \begin{equation} \sigma(t) = \sigma_0 \exp\left(-\frac{t}{\tau_R}\right) = E \varepsilon_0 \exp\left(-\frac{t}{\tau_R}\right) \end{equation} $$
応力は時間とともに指数関数的に減衰し、$t = \tau_R$ で初期値の $1/e \approx 36.8\%$ になります。$t \to \infty$ では $\sigma \to 0$ です。これはマクスウェルモデルの特徴であり、十分な時間が経てば応力が完全にゼロまで緩和することを意味します。
クリープ（$\sigma = \sigma_0 = \text{const.}$）の解:
$\dot{\sigma} = 0$ とすると、構成方程式は $\dot{\varepsilon} = \sigma_0 / \eta$ となり、積分すると、
$$ \begin{equation} \varepsilon(t) = \frac{\sigma_0}{E} + \frac{\sigma_0}{\eta} t \end{equation} $$
これは弾性ひずみ $\sigma_0 / E$ に加えて、ひずみが時間に対して線形に増加する解です。つまり、マクスウェルモデルはクリープの第II期（定常クリープ）の線形増加は表現できますが、第I期（遷移クリープ、ひずみ速度の減少）や第III期（加速クリープ）は表現できません。マクスウェルモデルは応力緩和の記述に適したモデルと言えます。
応力緩和をうまく記述するマクスウェルモデルですが、クリープの遷移挙動をうまく表現できません。次に、クリープの記述に適した別のモデルを見てみましょう。
フォークトモデル（Voigt model / Kelvin-Voigt model）——クリープを記述する
フォークトモデルは、バネとダッシュポットを並列に接続したモデルです。並列接続のポイントは、「両方に同じひずみが生じ、全体の応力は各要素の応力の和」であることです。
並列接続の条件は次の通りです。
$$ \varepsilon_{\text{spring}} = \varepsilon_{\text{dashpot}} = \varepsilon \quad \text{（ひずみは共通）} $$
$$ \sigma = \sigma_{\text{spring}} + \sigma_{\text{dashpot}} \quad \text{（応力は加算）} $$
各要素の構成則を代入すると、フォークトモデルの構成方程式が得られます。
$$ \begin{equation} \sigma = E \varepsilon + \eta \dot{\varepsilon} \end{equation} $$
クリープ（$\sigma = \sigma_0 = \text{const.}$）の解:
構成方程式を $\dot{\varepsilon}$ について整理すると、
$$ \dot{\varepsilon} = \frac{\sigma_0 – E\varepsilon}{\eta} $$
$\tau_c = \eta / E$ を遅延時間（retardation time）と定義し、この微分方程式を初期条件 $\varepsilon(0) = 0$ で解きます。
まず変数を変換します。$u = \sigma_0 / E – \varepsilon$ とおくと、$\dot{u} = -\dot{\varepsilon}$ であり、
$$ \dot{u} = -\frac{Eu}{\eta} = -\frac{u}{\tau_c} $$
これは指数減衰の微分方程式であり、$u(t) = u(0) e^{-t/\tau_c}$ と解けます。$u(0) = \sigma_0/E$ なので、
$$ \frac{\sigma_0}{E} – \varepsilon = \frac{\sigma_0}{E} e^{-t/\tau_c} $$
$\varepsilon$ について解くと、クリープの解が得られます。
$$ \begin{equation} \varepsilon(t) = \frac{\sigma_0}{E} \left(1 – e^{-t/\tau_c}\right) \end{equation} $$
この解は、ひずみが最終値 $\sigma_0 / E$（フックの法則による弾性ひずみと同じ値）に向かって指数関数的に漸近する曲線です。$t = \tau_c$ で最終値の約 $63.2\%$ に到達します。この挙動は遷移クリープ（第I期）の「ひずみ速度が時間とともに減少する」という特徴をよく捉えています。
ただし、フォークトモデルには限界もあります。$t \to \infty$ でひずみが一定値に収束してしまうため、定常クリープ（第II期の線形増加）を表現できません。また、応力緩和を考えると、$\varepsilon = \varepsilon_0$ の瞬間的な変形を加えることは、ダッシュポットに無限大の速度を要求するため物理的に実現できません（$\eta \dot{\varepsilon} \to \infty$）。つまり、フォークトモデルはクリープの記述に適しているが、応力緩和の記述には不向きです。
マクスウェルモデルとフォークトモデルは、それぞれ応力緩和とクリープの一方をうまく記述できますが、他方は苦手です。両方の長所を組み合わせたモデルはないのでしょうか？
標準線形固体モデル（Standard Linear Solid, SLS）
標準線形固体モデル（ゼナーモデルとも呼ばれます）は、マクスウェルモデルとフォークトモデルの限界を克服するために、3つの要素を組み合わせたモデルです。いくつかの等価な構成がありますが、ここではマクスウェル要素（$E_1$ と $\eta$ の直列接続）にバネ $E_2$ を並列に加えた構成を考えます。
この構成では、全体の応力はマクスウェル要素の応力 $\sigma_M$ と並列バネの応力 $\sigma_2$ の和です。
$$ \sigma = \sigma_M + \sigma_2 = \sigma_M + E_2 \varepsilon $$
マクスウェル要素の構成方程式は先ほど導出した通り、
$$ \dot{\varepsilon} = \frac{\dot{\sigma}_M}{E_1} + \frac{\sigma_M}{\eta} $$
$\sigma_M = \sigma – E_2 \varepsilon$ を代入し、整理すると、SLSモデルの構成方程式が得られます。
$$ \begin{equation} \sigma + \frac{\eta}{E_1} \dot{\sigma} = E_2 \varepsilon + \frac{\eta(E_1 + E_2)}{E_1} \dot{\varepsilon} \end{equation} $$
クリープの解:
$\sigma = \sigma_0 = \text{const.}$（$\dot{\sigma} = 0$）として微分方程式を解くと、
$$ \begin{equation} \varepsilon(t) = \sigma_0 \left[\frac{1}{E_1 + E_2} + \frac{1}{E_2} \cdot \frac{E_1}{E_1 + E_2} \left(1 – e^{-t/\tau’}\right)\right] \end{equation} $$
ここで $\tau’ = \eta(E_1 + E_2) / (E_1 E_2)$ です。初期の瞬時弾性ひずみ $\sigma_0 / (E_1 + E_2)$ に、遷移クリープに相当する指数関数的な増加が加わり、最終的に $\sigma_0 / E_2$ に漸近します。
応力緩和の解:
$\varepsilon = \varepsilon_0 = \text{const.}$（$\dot{\varepsilon} = 0$）として微分方程式を解くと、
$$ \begin{equation} \sigma(t) = E_2 \varepsilon_0 + (E_1 + E_2) \varepsilon_0 \cdot \frac{E_1}{E_1 + E_2} \cdot e^{-t E_1/\eta} \end{equation} $$
整理すると、
$$ \begin{equation} \sigma(t) = \varepsilon_0 \left[E_2 + E_1 e^{-t/\tau_R}\right] \end{equation} $$
ここで $\tau_R = \eta / E_1$ です。重要な点は、$t \to \infty$ で $\sigma \to E_2 \varepsilon_0 \neq 0$ であることです。マクスウェルモデルでは応力がゼロまで緩和しましたが、SLSモデルでは並列バネ $E_2$ の存在により、有限の残留応力が保持されます。これは現実の固体の応力緩和挙動（応力が完全にはゼロにならない）をより正しく表現しています。
3つの粘弾性モデルの特徴を表にまとめます。
| モデル | 構成 | クリープ | 応力緩和 | 長所 | 短所 |
|---|---|---|---|---|---|
| マクスウェル | バネ + ダッシュポット（直列） | 線形増加（永久ひずみ） | 指数減衰 → 0 | 応力緩和を良く記述 | クリープが非現実的 |
| フォークト | バネ + ダッシュポット（並列） | 指数的漸近 → 有限値 | 記述困難 | 遷移クリープを良く記述 | 応力緩和・瞬時弾性応答が不可 |
| SLS | マクスウェル + 並列バネ | 瞬時弾性 + 指数的漸近 | 指数減衰 → 有限残留応力 | 両方を定性的に記述 | パラメータが3つ必要 |
粘弾性モデルによる理論的な理解が深まったところで、次は工学的に非常に重要な「寿命推定」の方法であるラーソン・ミラーパラメータについて解説します。
ラーソン・ミラーパラメータ——温度と時間をまとめる魔法の数
短期試験から長期寿命を予測する必要性
実際の発電プラントやジェットエンジンでは、部材は数万〜数十万時間（数年〜数十年）にわたってクリープにさらされます。しかし、そのような長期間のクリープ試験を実施するのは現実的ではありません。そこで、高温・短時間の試験結果から、低温・長時間の寿命を外挿する方法が必要になります。
この目的のために考案されたのが、時間-温度パラメータ法（time-temperature parameter method）です。その中で最も広く使われているのがラーソン・ミラーパラメータ（Larson-Miller Parameter, LMP）です。
ラーソン・ミラーパラメータの導出
ラーソン・ミラーパラメータの理論的基礎は、クリープ速度（や破断時間）がアレニウス型の温度依存性を持つという事実にあります。クリープ破断時間 $t_r$ が応力 $\sigma$ と温度 $T$ に依存するとき、アレニウス型の関係は次のように書けます。
$$ \begin{equation} t_r = t_0 \exp\left(\frac{Q}{RT}\right) \end{equation} $$
ここで $t_0$ は応力に依存する定数、$Q$ はクリープの見かけの活性化エネルギーです。この式の両辺に自然対数を取ります。
$$ \ln t_r = \ln t_0 + \frac{Q}{RT} $$
常用対数に変換すると（$\ln x = 2.303 \log_{10} x$）、
$$ 2.303 \log t_r = 2.303 \log t_0 + \frac{Q}{RT} $$
両辺を $T$ 倍し、$2.303 R$ で割ると、
$$ T (\log t_r – \log t_0) = \frac{Q}{2.303 R} $$
ここで $C = -\log t_0$ とおくと（$C$ は材料定数で、多くの金属では $C \approx 20$ です）、
$$ \begin{equation} \text{LMP} = T (C + \log t_r) \end{equation} $$
この値がラーソン・ミラーパラメータです。$T$ は絶対温度 [K]、$t_r$ はクリープ破断時間 [h]、$C$ は材料定数です。
LMP の優れた点は、異なる温度と時間の組み合わせであっても、同じ応力レベルでは LMP がほぼ一定になることです。つまり、高温・短時間の試験で得たデータを用いて、LMP vs. 応力のマスターカーブを作成すれば、低温・長時間の破断寿命を推定できるのです。
例を挙げましょう。ある合金の $C = 20$ として、以下の2つの条件は同じ LMP を与えます。
- 条件A: $T = 1000$ K, $t_r = 100$ h → LMP = $1000 \times (20 + \log 100) = 1000 \times 22 = 22{,}000$
- 条件B: $T = 900$ K, $t_r = t_r^*$ h → LMP = $900 \times (20 + \log t_r^*) = 22{,}000$
条件Bから $t_r^*$ を求めると、$20 + \log t_r^* = 22{,}000 / 900 = 24.44$ より $\log t_r^* = 4.44$、つまり $t_r^* \approx 27{,}542$ h（約3.1年）となります。温度が 100 K 低下するだけで、同じ応力での寿命が 100 h から約 27,500 h へと大幅に延びることがわかります。
LMP の工学的活用
LMP を用いた寿命推定の手順は次の通りです。
- マスターカーブの作成: 異なる温度・応力条件で複数のクリープ破断試験を実施し、各データ点の LMP と応力 $\sigma$ を計算する
- プロット: 横軸に LMP、縦軸に応力 $\sigma$ を取ってプロットする。データ点は材料固有の1本のマスターカーブに乗る
- 寿命推定: 設計応力と使用温度から LMP を計算し、マスターカーブから破断時間 $t_r$ を外挿する
ただし、LMP法には注意すべき限界もあります。活性化エネルギー $Q$ が温度や応力によって変化する場合（例えばクリープ機構が切り替わる場合）、外挿の精度が低下します。また、$C$ の値は材料や試験条件によって 15〜25 程度の幅があり、$C$ の選択が寿命推定に大きな影響を与えます。
ラーソン・ミラーパラメータの理論を押さえたところで、次は高温設計における実際の応用について見ていきましょう。
高温設計への応用——クリープと応力緩和を考慮した設計指針
許容クリープ応力
高温機器の設計では、室温の降伏応力や引張強さに加えて、クリープに基づく許容応力が設計の支配因子となります。代表的な基準は以下の2つです。
-
10万時間クリープ破断強度: 100,000 時間（約11.4年）でクリープ破断を生じる応力。長期使用を前提とするボイラーや圧力容器の設計基準として広く使われます。
-
10万時間1%クリープひずみ強度: 100,000 時間で累積クリープひずみが 1% に達する応力。変形が制限される部材（薄肉容器のシール面など）の設計に用います。
これらの値は、前述の LMP やクリープ破断試験データから外挿して求めます。設計応力は、これらの値にさらに安全率を掛けて決定します。
クリープ破壊と余寿命評価
長期間運転された高温機器の余寿命を評価するには、以下のようなアプローチが用いられます。
ロビンソンの累積損傷則（Robinson’s life fraction rule）は、疲労破壊のマイナー則と類似した概念です。時間 $t_i$ の間、応力 $\sigma_i$ でクリープが進行し、その条件での破断寿命が $t_{r,i}$ であるとき、損傷分率の総和が1に達したときに破断するとします。
$$ \begin{equation} \sum_i \frac{t_i}{t_{r,i}} = 1 \end{equation} $$
連続的に応力・温度が変化する場合は積分形になります。
$$ \begin{equation} \int_0^{t_f} \frac{dt}{t_r(\sigma, T)} = 1 \end{equation} $$
ここで $t_f$ が予測される破断時間です。
材料選定と組織制御
クリープ強度を高めるための材料設計の指針を整理しておきます。
- 結晶粒の粗大化: 拡散クリープは粒径に反比例するため、粒径を大きくすることでクリープ抵抗が向上します。極限として、タービンブレードには単結晶（粒界なし）が使用されます
- 析出強化: $\gamma’$ 相（ニッケル基超合金の場合）や炭化物などの安定な析出物が転位の運動を阻害します
- 固溶強化: W, Mo, Re などの高融点元素を添加し、格子を歪ませることで転位の運動を抑制します
- 融点の高い材料の選択: ホモロガス温度が下がるため、同じ使用温度でもクリープが抑制されます
また、熱力学第1法則から理解されるように、熱機関の効率を上げるには作動温度を上げることが望ましく、これがより高い耐クリープ性を持つ材料の開発を駆動する背景となっています。
以上で理論的な解説を一通り終えました。ここからは、Python を用いてクリープ曲線、粘弾性モデルの応答、応力緩和、ラーソン・ミラーパラメータの各テーマを実装・可視化し、理論の理解を深めましょう。
Python 実装 (1): クリープ曲線の描画と3段階の可視化
まず、クリープ曲線の3段階を1つのグラフで可視化し、各段階の特徴を目で確認しましょう。先ほど示した経験式 $\varepsilon_c(t) = A t^{1/3} + B t + C(e^{\gamma t} – 1)$ をPythonで実装します。
import numpy as np
import matplotlib.pyplot as plt
# --- クリープ曲線のパラメータ ---
A = 0.008 # 遷移クリープの係数
B = 2.0e-5 # 定常クリープ速度 [1/h]
C = 1.0e-6 # 加速クリープの係数
gamma = 0.0006 # 加速クリープの指数
# 時間 [h]
t = np.linspace(0, 10000, 5000)
# 瞬時弾性ひずみ
eps_elastic = 0.002
# 各段階のひずみ成分
eps_primary = A * t ** (1 / 3) # 第I期: 遷移クリープ
eps_secondary = B * t # 第II期: 定常クリープ
eps_tertiary = C * (np.exp(gamma * t) - 1) # 第III期: 加速クリープ
# 全クリープひずみ
eps_creep = eps_primary + eps_secondary + eps_tertiary
eps_total = eps_elastic + eps_creep
# クリープ速度（数値微分）
dt = t[1] - t[0]
creep_rate = np.gradient(eps_creep, dt)
# --- 描画 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: クリープ曲線
ax = axes[0]
ax.plot(t, eps_total * 100, 'k-', linewidth=2, label='Total strain')
ax.plot(t, (eps_elastic + eps_primary) * 100, 'b--', alpha=0.6, label='Primary only')
ax.plot(t, (eps_elastic + eps_secondary) * 100, 'g--', alpha=0.6, label='Steady-state only')
ax.axhline(y=eps_elastic * 100, color='gray', linestyle=':', alpha=0.5, label='Elastic strain')
# 段階の境界を概略で表示
t_boundary1, t_boundary2 = 1500, 7500
ax.axvline(x=t_boundary1, color='orange', linestyle='--', alpha=0.4)
ax.axvline(x=t_boundary2, color='red', linestyle='--', alpha=0.4)
ax.text(t_boundary1 / 2, ax.get_ylim()[0] + 0.01, 'Stage I', ha='center', fontsize=10, color='blue')
ax.text((t_boundary1 + t_boundary2) / 2, ax.get_ylim()[0] + 0.01, 'Stage II', ha='center', fontsize=10, color='green')
ax.text((t_boundary2 + t[-1]) / 2, ax.get_ylim()[0] + 0.01, 'Stage III', ha='center', fontsize=10, color='red')
ax.set_xlabel('Time [h]', fontsize=12)
ax.set_ylabel('Strain [%]', fontsize=12)
ax.set_title('Creep Curve (3 Stages)', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# 右: クリープ速度
ax2 = axes[1]
ax2.semilogy(t[1:], creep_rate[1:], 'r-', linewidth=2)
ax2.axvline(x=t_boundary1, color='orange', linestyle='--', alpha=0.4)
ax2.axvline(x=t_boundary2, color='red', linestyle='--', alpha=0.4)
ax2.set_xlabel('Time [h]', fontsize=12)
ax2.set_ylabel('Creep rate [1/h]', fontsize=12)
ax2.set_title('Creep Rate vs Time', fontsize=13)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のグラフでは、3段階のクリープ挙動が明確に可視化されています。初期（Stage I）ではひずみの増加率が徐々に小さくなり（上に凸の曲線）、中期（Stage II）ではほぼ線形の増加（定常クリープ）が見られ、後期（Stage III）では指数関数的な加速が始まっています。黒い実線（Total strain）が3つの成分の合成であり、青い破線（Primary only）と緑の破線（Steady-state only）との比較で、各成分がどの時間帯で支配的かが一目でわかります。
右のグラフでは、クリープ速度の時間変化を対数スケールで表示しています。第I期で速度が急速に低下し、第II期でほぼ一定の最小クリープ速度に達した後、第III期で速度が再び上昇する「バスタブ型」の曲線が確認できます。この最小クリープ速度こそが、設計で重要となる定常クリープ速度 $\dot{\varepsilon}_s$ に対応しています。
Python 実装 (2): マクスウェルモデルとフォークトモデルの応答
次に、マクスウェルモデルとフォークトモデルに一定応力（クリープ試験に相当）を加えたときの応答と、一定ひずみ（応力緩和試験に相当）を加えたときの応答を比較します。
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ ---
E = 200e3 # ヤング率 [MPa]
eta = 1e7 # 粘性係数 [MPa·h]
sigma_0 = 100.0 # 一定応力 [MPa]
eps_0 = sigma_0 / E # 一定ひずみ [-]
tau = eta / E # 緩和時間・遅延時間 [h]
print(f"緩和時間 / 遅延時間 τ = η/E = {tau:.1f} h")
t = np.linspace(0, 300, 1000)
# === クリープ応答（一定応力 σ₀） ===
# マクスウェル: ε(t) = σ₀/E + σ₀t/η
eps_maxwell_creep = sigma_0 / E + sigma_0 * t / eta
# フォークト: ε(t) = (σ₀/E)(1 - exp(-t/τ))
eps_voigt_creep = (sigma_0 / E) * (1 - np.exp(-t / tau))
# === 応力緩和応答（一定ひずみ ε₀） ===
# マクスウェル: σ(t) = Eε₀ exp(-t/τ)
sigma_maxwell_relax = E * eps_0 * np.exp(-t / tau)
# フォークト: 応力緩和は本質的に瞬間変形を要求するため記述困難
# σ = Eε₀ + η·δ(0) → 概念上の困難。ここでは σ = Eε₀ (一定) として表示
sigma_voigt_relax = np.full_like(t, E * eps_0)
# --- 描画 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: クリープ応答
ax = axes[0]
ax.plot(t, eps_maxwell_creep * 1000, 'b-', linewidth=2, label='Maxwell')
ax.plot(t, eps_voigt_creep * 1000, 'r-', linewidth=2, label='Voigt (Kelvin)')
ax.axhline(y=eps_0 * 1000, color='gray', linestyle=':', alpha=0.5,
label=f'σ₀/E = {eps_0*1000:.3f} ×10⁻³')
ax.set_xlabel('Time [h]', fontsize=12)
ax.set_ylabel('Strain [×10⁻³]', fontsize=12)
ax.set_title(f'Creep Response (σ₀ = {sigma_0} MPa)', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# 右: 応力緩和応答
ax2 = axes[1]
ax2.plot(t, sigma_maxwell_relax, 'b-', linewidth=2, label='Maxwell')
ax2.plot(t, sigma_voigt_relax, 'r--', linewidth=2, alpha=0.7,
label='Voigt (constant)')
ax2.axhline(y=0, color='gray', linestyle=':', alpha=0.3)
ax2.set_xlabel('Time [h]', fontsize=12)
ax2.set_ylabel('Stress [MPa]', fontsize=12)
ax2.set_title(f'Stress Relaxation (ε₀ = {eps_0*1000:.3f} ×10⁻³)', fontsize=13)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のクリープ応答のグラフから、2つのモデルの本質的な違いがはっきり読み取れます。マクスウェルモデル（青線）は、初期の弾性ひずみの後、一定速度でひずみが線形に増加し続けます。これは定常クリープ（第II期）の挙動を表現していますが、初期の遷移的な減速が見られません。一方、フォークトモデル（赤線）は、ゼロから始まって弾性ひずみ $\sigma_0/E$ に漸近していく指数関数的な曲線です。遷移クリープの「ひずみ速度の減少」は良く捉えていますが、定常クリープ（線形増加）に移行することがなく、有限値で飽和してしまいます。
右の応力緩和のグラフでは、マクスウェルモデル（青線）が応力の指数減衰を示すのに対し、フォークトモデル（赤破線）は応力が一定値を保ったままです。フォークトモデルで瞬間的なひずみを加えるにはダッシュポットに無限の応力が必要となるため、現実的な応力緩和の解が存在しないことを反映しています。この比較から、応力緩和の記述にはマクスウェルモデルが適していることが明確にわかります。
Python 実装 (3): 標準線形固体モデルと応力緩和シミュレーション
次に、標準線形固体（SLS）モデルのクリープ応答と応力緩和応答を実装し、マクスウェルモデルおよびフォークトモデルとの違いを可視化します。加えて、Norton の冪乗則に基づく応力緩和の数値シミュレーションも実装します。
import numpy as np
import matplotlib.pyplot as plt
# === SLS モデルのパラメータ ===
E1 = 150e3 # マクスウェル要素のバネ [MPa]
E2 = 50e3 # 並列バネ [MPa]
eta = 5e6 # ダッシュポットの粘性係数 [MPa·h]
sigma_0 = 100.0 # 一定応力 [MPa]
# 緩和時間
tau_R = eta / E1
# 遅延時間
tau_prime = eta * (E1 + E2) / (E1 * E2)
t = np.linspace(0, 500, 2000)
# === SLS クリープ ===
eps_sls_creep = sigma_0 * (1 / (E1 + E2) + (E1 / (E2 * (E1 + E2))) * (1 - np.exp(-t / tau_prime)))
# === SLS 応力緩和 ===
eps_0_sls = sigma_0 / (E1 + E2) # 瞬時弾性ひずみ
sigma_sls_relax = eps_0_sls * (E2 + E1 * np.exp(-t / tau_R))
# --- 比較用: マクスウェルモデル ---
E_mw = E1 + E2
tau_mw = eta / E_mw
eps_mw_creep = sigma_0 / E_mw + sigma_0 * t / eta
sigma_mw_relax = E_mw * eps_0_sls * np.exp(-t / tau_mw)
# --- 描画 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: クリープ応答の比較
ax = axes[0]
ax.plot(t, eps_sls_creep * 1000, 'b-', linewidth=2, label='SLS model')
ax.plot(t, eps_mw_creep * 1000, 'r--', linewidth=2, alpha=0.7, label='Maxwell model')
ax.axhline(y=sigma_0 / E2 * 1000, color='blue', linestyle=':', alpha=0.4,
label=f'SLS limit: σ₀/E₂ = {sigma_0/E2*1000:.2f}×10⁻³')
ax.set_xlabel('Time [h]', fontsize=12)
ax.set_ylabel('Strain [×10⁻³]', fontsize=12)
ax.set_title('Creep: SLS vs Maxwell', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# 右: 応力緩和の比較
ax2 = axes[1]
ax2.plot(t, sigma_sls_relax, 'b-', linewidth=2, label='SLS model')
ax2.plot(t, sigma_mw_relax, 'r--', linewidth=2, alpha=0.7, label='Maxwell model')
ax2.axhline(y=E2 * eps_0_sls, color='blue', linestyle=':', alpha=0.4,
label=f'SLS residual: E₂ε₀ = {E2*eps_0_sls:.1f} MPa')
ax2.axhline(y=0, color='gray', linestyle=':', alpha=0.3)
ax2.set_xlabel('Time [h]', fontsize=12)
ax2.set_ylabel('Stress [MPa]', fontsize=12)
ax2.set_title('Stress Relaxation: SLS vs Maxwell', fontsize=13)
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のクリープ応答の比較から、SLSモデル（青実線）がマクスウェルモデル（赤破線）とは質的に異なる挙動を示すことがわかります。マクスウェルモデルは無限にひずみが増加し続けるのに対し、SLSモデルは $\sigma_0 / E_2$ の有限値に漸近します。SLSモデルは「初期の弾性ひずみ → 遷移クリープ → 飽和」という流れを表現しており、フォークトモデルの特徴（飽和挙動）をマクスウェルモデルの特徴（初期弾性応答）と組み合わせています。
右の応力緩和の比較では、マクスウェルモデルの応力がゼロに向かって減衰するのに対し、SLSモデルは $E_2 \varepsilon_0$ の有限残留応力に漸近しています。現実の固体材料は応力が完全にはゼロにならない場合が多く、SLSモデルの予測はより実際に近い挙動です。
続いて、Norton の冪乗則に基づく応力緩和をオイラー法で数値的にシミュレーションしてみましょう。
import numpy as np
import matplotlib.pyplot as plt
# === Norton冪乗則に基づく応力緩和シミュレーション ===
# dσ/dt = -E * ε_dot_creep = -E * A' * σ^n * exp(-Q/(RT))
# パラメータ
E = 150e3 # ヤング率 [MPa]
A_prime = 1e-18 # Norton則の係数 [MPa^(-n) / h]
n_values = [3, 5, 7] # 応力指数
Q = 250e3 # 活性化エネルギー [J/mol]
R_gas = 8.314 # 気体定数 [J/(mol·K)]
T = 873 # 温度 [K] (600°C)
sigma_0 = 150.0 # 初期応力 [MPa]
# 時間設定
t_max = 50000 # [h]
dt = 1.0 # 時間刻み [h]
t_arr = np.arange(0, t_max, dt)
N_steps = len(t_arr)
# アレニウス項
arrhenius = np.exp(-Q / (R_gas * T))
fig, ax = plt.subplots(figsize=(10, 6))
for n in n_values:
sigma = np.zeros(N_steps)
sigma[0] = sigma_0
for i in range(N_steps - 1):
# Norton冪乗則によるクリープ速度
eps_dot_creep = A_prime * sigma[i] ** n * arrhenius
# 応力緩和の微分方程式
dsigma_dt = -E * eps_dot_creep
sigma[i + 1] = sigma[i] + dsigma_dt * dt
# 応力が負になるのを防止
if sigma[i + 1] < 0:
sigma[i + 1] = 0
ax.plot(t_arr / 1000, sigma, linewidth=2, label=f'n = {n}')
ax.set_xlabel('Time [×10³ h]', fontsize=12)
ax.set_ylabel('Stress [MPa]', fontsize=12)
ax.set_title(f'Stress Relaxation with Norton Power Law (T = {T} K)', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_ylim(bottom=0)
plt.tight_layout()
plt.show()
このグラフでは、Norton則の応力指数 $n$ の違いが応力緩和の挙動にどう影響するかを比較しています。$n = 3$（拡散クリープ寄り）の場合は緩やかで均一な緩和が見られるのに対し、$n = 7$（転位クリープ寄り）の場合は初期に急激な応力降下が起きた後、低応力領域ではクリープ速度が $\sigma^n$ に比例して急激に減少するため、なかなかゼロに達しません。$n$ が大きいほど応力への依存性が強い（高応力時にクリープが速い）ため、初期の緩和が速くなりますが、応力が低くなるとクリープ速度が急減し、残留応力が長く保持される傾向があります。これは実際のボルト締結部材での応力緩和挙動と定性的に一致しており、$n$ の値が設計上重要なパラメータであることを示しています。
Python 実装 (4): ラーソン・ミラーパラメータによる寿命推定
最後に、ラーソン・ミラーパラメータ（LMP）を用いたクリープ寿命推定を実装します。仮想的な材料のマスターカーブを作成し、異なる温度条件での寿命を推定します。
import numpy as np
import matplotlib.pyplot as plt
# === ラーソン・ミラーパラメータ（LMP）による寿命推定 ===
C_lm = 20 # ラーソン・ミラー定数
def calc_lmp(T, t_r, C=20):
"""LMPを計算する (T: 絶対温度[K], t_r: 破断時間[h])"""
return T * (C + np.log10(t_r))
def calc_rupture_time(T, lmp, C=20):
"""LMPから破断時間を逆算する"""
log_tr = lmp / T - C
return 10 ** log_tr
# --- 仮想的な実験データ（LMP vs 応力のマスターカーブ用） ---
# 異なる温度・時間でのクリープ破断試験データ
test_data = [
# (温度[K], 破断時間[h], 破断応力[MPa])
(1073, 10, 350), # 800°C
(1073, 100, 280),
(1073, 1000, 220),
(1073, 5000, 180),
(1023, 50, 320), # 750°C
(1023, 500, 260),
(1023, 5000, 200),
(1023, 20000, 170),
(973, 100, 340), # 700°C
(973, 1000, 270),
(973, 10000, 210),
(973, 50000, 175),
(923, 500, 330), # 650°C
(923, 5000, 265),
(923, 50000, 205),
]
# LMPを計算
lmp_data = []
sigma_data = []
temp_labels = []
for T, t_r, sigma in test_data:
lmp_val = calc_lmp(T, t_r, C_lm)
lmp_data.append(lmp_val)
sigma_data.append(sigma)
temp_labels.append(T)
lmp_data = np.array(lmp_data)
sigma_data = np.array(sigma_data)
temp_labels = np.array(temp_labels)
# マスターカーブのフィッティング（2次多項式: σ = a₀ + a₁·LMP + a₂·LMP²）
coeffs = np.polyfit(lmp_data, sigma_data, 2)
lmp_fit = np.linspace(min(lmp_data) * 0.95, max(lmp_data) * 1.05, 200)
sigma_fit = np.polyval(coeffs, lmp_fit)
# --- 寿命推定 ---
# 設計条件: σ = 200 MPa での各温度の破断時間を推定
sigma_design = 200 # [MPa]
# マスターカーブから σ_design に対応する LMP を求める（逆算）
# σ = a₂·LMP² + a₁·LMP + a₀ = σ_design を解く
poly_eq = [coeffs[0], coeffs[1], coeffs[2] - sigma_design]
lmp_roots = np.roots(poly_eq)
lmp_target = lmp_roots[lmp_roots > 0].real.max() # 正の実数根を選択
print(f"設計応力 σ = {sigma_design} MPa に対応する LMP = {lmp_target:.0f}")
print(f"\n--- 各温度での推定破断時間 ---")
design_temps = [923, 973, 1023, 1073]
for T in design_temps:
t_r_est = calc_rupture_time(T, lmp_target, C_lm)
years = t_r_est / 8760
print(f" T = {T} K ({T-273}°C): t_r = {t_r_est:,.0f} h ({years:.1f} 年)")
# --- 描画 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: LMP マスターカーブ
ax = axes[0]
unique_temps = sorted(set(temp_labels))
markers = ['o', 's', '^', 'D']
colors_temp = ['#e74c3c', '#e67e22', '#2ecc71', '#3498db']
for i, T in enumerate(unique_temps):
mask = temp_labels == T
ax.scatter(lmp_data[mask] / 1000, sigma_data[mask],
marker=markers[i], color=colors_temp[i], s=80, zorder=5,
label=f'{T} K ({T-273}°C)')
ax.plot(lmp_fit / 1000, sigma_fit, 'k-', linewidth=2, alpha=0.7, label='Master curve')
ax.axhline(y=sigma_design, color='purple', linestyle='--', alpha=0.5,
label=f'σ = {sigma_design} MPa')
ax.axvline(x=lmp_target / 1000, color='purple', linestyle=':', alpha=0.3)
ax.set_xlabel('LMP / 1000', fontsize=12)
ax.set_ylabel('Stress [MPa]', fontsize=12)
ax.set_title('Larson-Miller Parameter Master Curve', fontsize=13)
ax.legend(fontsize=9, loc='upper right')
ax.grid(True, alpha=0.3)
# 右: 温度ごとの破断時間
ax2 = axes[1]
temps_plot = np.linspace(873, 1123, 100)
tr_plot = calc_rupture_time(temps_plot, lmp_target, C_lm)
ax2.semilogy(temps_plot - 273, tr_plot, 'b-', linewidth=2)
ax2.axhline(y=1e5, color='red', linestyle='--', alpha=0.5, label='100,000 h')
for T in design_temps:
t_r_est = calc_rupture_time(T, lmp_target, C_lm)
ax2.plot(T - 273, t_r_est, 'ko', markersize=8)
ax2.annotate(f'{t_r_est:,.0f} h', xy=(T - 273, t_r_est),
xytext=(10, 5), textcoords='offset points', fontsize=8)
ax2.set_xlabel('Temperature [°C]', fontsize=12)
ax2.set_ylabel('Rupture time [h]', fontsize=12)
ax2.set_title(f'Estimated Rupture Life (σ = {sigma_design} MPa)', fontsize=13)
ax2.legend(fontsize=11)
ax2.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.show()
左のグラフは、LMPマスターカーブを示しています。異なる温度（色分け）で得られたクリープ破断データが、LMP を横軸に取ることで1本の曲線上にほぼ集約されていることが確認できます。これがLMP法の最大の利点——温度と時間の効果を1つのパラメータにまとめることで、データの統一的な整理が可能になるということです。紫の水平破線が設計応力 200 MPa を示しており、マスターカーブとの交点から対応する LMP が求まります。
右のグラフは、設計応力 200 MPa における温度ごとの推定破断時間を対数スケールで示しています。温度が上がるにつれて破断時間が劇的に短くなることが一目瞭然です。赤い破線で示した 100,000 時間（約 11.4 年）のラインを安全基準とすると、使用可能な最高温度の目安が読み取れます。わずか 50°C の温度差が破断寿命を1桁以上変化させることから、高温機器における温度管理の重要性が数値的に裏付けられています。
Python 実装 (5): クリープ曲線パラメータのフィッティング
最後に、仮想的な「実験データ」に対してクリープ曲線の経験式をフィッティングし、パラメータを推定する方法を示します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# --- クリープ曲線の経験式 ---
def creep_model(t, A, B, C, gamma):
"""3段階クリープ曲線モデル"""
return A * t ** (1/3) + B * t + C * (np.exp(gamma * t) - 1)
# 仮想実験データの生成（真のパラメータにノイズを加える）
np.random.seed(42)
t_exp = np.linspace(1, 8000, 80) # 測定時刻
# 真のパラメータ
A_true, B_true, C_true, gamma_true = 0.006, 3.0e-5, 5.0e-7, 0.0005
eps_true = creep_model(t_exp, A_true, B_true, C_true, gamma_true)
noise = np.random.normal(0, 0.0003, len(t_exp))
eps_exp = eps_true + noise
# フィッティング
p0 = [0.005, 1e-5, 1e-6, 1e-4] # 初期推定値
bounds = ([0, 0, 0, 0], [0.1, 1e-3, 1e-3, 0.01]) # パラメータの範囲
popt, pcov = curve_fit(creep_model, t_exp, eps_exp, p0=p0, bounds=bounds, maxfev=10000)
perr = np.sqrt(np.diag(pcov)) # 標準誤差
print("--- フィッティング結果 ---")
param_names = ['A', 'B', 'C', 'γ']
true_vals = [A_true, B_true, C_true, gamma_true]
for name, true_val, opt, err in zip(param_names, true_vals, popt, perr):
print(f" {name}: 推定値 = {opt:.4e}, 真値 = {true_val:.4e}, 誤差 = {err:.4e}")
# 定常クリープ速度
print(f"\n定常クリープ速度 (推定): B = {popt[1]:.4e} [1/h]")
print(f"定常クリープ速度 (真値): B = {B_true:.4e} [1/h]")
# --- 描画 ---
t_fine = np.linspace(1, 8000, 1000)
eps_fitted = creep_model(t_fine, *popt)
fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(t_exp, eps_exp * 100, s=20, color='gray', alpha=0.7, label='Experimental data')
ax.plot(t_fine, eps_fitted * 100, 'r-', linewidth=2, label='Fitted curve')
ax.plot(t_fine, creep_model(t_fine, A_true, B_true, C_true, gamma_true) * 100,
'b--', linewidth=1.5, alpha=0.5, label='True curve')
# 残差を小さなサブプロットで表示
ax_inset = ax.inset_axes([0.15, 0.55, 0.35, 0.35])
residual = (eps_exp - creep_model(t_exp, *popt)) * 100
ax_inset.scatter(t_exp, residual, s=10, color='green', alpha=0.6)
ax_inset.axhline(y=0, color='k', linestyle='-', alpha=0.3)
ax_inset.set_xlabel('Time [h]', fontsize=8)
ax_inset.set_ylabel('Residual [%]', fontsize=8)
ax_inset.set_title('Residuals', fontsize=9)
ax_inset.tick_params(labelsize=7)
ax.set_xlabel('Time [h]', fontsize=12)
ax.set_ylabel('Creep strain [%]', fontsize=12)
ax.set_title('Creep Curve Fitting', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
フィッティング結果のグラフでは、仮想的な実験データ（灰色の散布点）に対して、赤い実線のフィッティング曲線が非常によくフィットしていることが確認できます。青い破線の真の曲線ともほぼ一致しており、SciPy の curve_fit による非線形最小二乗法が適切にパラメータを推定できていることがわかります。
インセットの残差プロットでは、残差がゼロ付近にランダムに分布しており、系統的な偏りが見られません。これは、選択した経験式がデータの構造を良く捉えていることを示しています。フィッティングから得られる定常クリープ速度 $B$ は、設計における許容クリープ速度との比較や、余寿命評価のための外挿に直接利用できる重要なパラメータです。
まとめ
本記事では、クリープと応力緩和という時間依存変形の2大テーマについて、基礎概念から粘弾性モデル、工学的応用まで体系的に解説しました。
- クリープは一定応力下でひずみが時間とともに増加する現象であり、典型的に遷移クリープ（第I期）→ 定常クリープ（第II期）→ 加速クリープ（第III期）の3段階を経てクリープ破断に至る
- クリープのメカニズムは、低応力・高温では拡散クリープ（ナバロ・ヘリング、コブル）、高応力では転位クリープ（上昇運動+すべり）が支配的であり、変形機構図で整理される
- 応力緩和は一定ひずみ下で応力が時間とともに減少する現象であり、クリープと表裏一体の関係にある
- マクスウェルモデルはバネとダッシュポットの直列接続で、応力緩和を指数関数的に記述する。フォークトモデルは並列接続で、遷移クリープの飽和挙動を記述する。標準線形固体モデルは両者の長所を併せ持ち、残留応力のある応力緩和と瞬時弾性応答を伴うクリープを表現できる
- ラーソン・ミラーパラメータは温度と時間の効果を1つのパラメータに集約し、高温・短時間の試験データから低温・長時間のクリープ寿命を外挿する強力なツールである
高温機器の設計では、応力とひずみに基づく静的強度評価だけでなく、クリープによる時間依存変形と応力緩和による締結力の低下を考慮することが不可欠です。
次のステップとして、以下の記事も参考にしてください。
- 疲労破壊 — クリープと疲労が同時に作用するクリープ-疲労相互作用は実機の重要な課題です
- ニッケル基超合金 — 高温クリープに対抗するために開発された最先端材料の設計思想を学べます
- 拡散（フィックの法則） — 拡散クリープのメカニズムをより深く理解するための基礎理論です