自動車のドライブシャフトは、エンジンの回転トルクを車輪に伝える部品です。走行中、このシャフトにはねじりによるせん断応力と自重や振動による曲げ応力が同時に作用します。では、このシャフトが壊れるかどうかを判定するとき、ねじりと曲げの応力をそれぞれ単独で評価すれば十分でしょうか？
答えは「いいえ」です。複数の応力が同時に作用する状態 — すなわち組合せ応力（combined stress）— では、個々の応力成分を単独で見ても安全性は判断できません。応力の「合成効果」を正しく評価するために、主応力、最大せん断応力、そしてフォンミーゼス応力といった指標が必要になります。
組合せ応力の理論を理解すると、以下のような幅広い場面で設計の根拠を明確にできます。
- 機械設計: 軸・歯車・クランクシャフトなど、曲げ・ねじり・引張が複合する部品の強度評価に不可欠です
- 構造設計: 橋梁・建築フレームの接合部は、軸力・せん断力・曲げモーメントが同時に作用する典型例です
- 航空宇宙工学: 航空機の翼桁やロケットの推進剤タンクは、内圧・曲げ・熱応力が複合した過酷な応力状態にさらされます
- 破壊解析・疲労設計: 応力集中部や疲労破壊の評価には、多軸応力場での等価応力の算出が前提です
本記事の内容
- 多軸応力状態の概要と応力テンソルの定義
- 応力の座標変換公式とその導出
- 主応力・主方向の固有値問題としての定式化
- モールの応力円（2D・3D）の構成と読み方
- 降伏条件: トレスカ条件とフォンミーゼス条件の導出と比較
- 安全率の計算と設計への応用
- Pythonによる主応力計算・モールの応力円描画・フォンミーゼス応力の等高線図・降伏条件の比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 応力とひずみの定義 — 垂直応力とせん断応力の基本概念
- フックの法則と弾性係数 — 弾性域での応力-ひずみ関係
- 主応力とモールの応力円 — 2次元の応力変換と主応力の基礎
- 平面応力と平面ひずみ — 3次元問題の2次元近似
- ヤング率とポアソン比 — 弾性定数の物理的意味
- 固有値・固有ベクトル — 主応力の導出に必要な線形代数
多軸応力状態とは
単軸から多軸へ — なぜ「組合せ」が問題になるのか
応力とひずみの定義で学んだ単純引張試験では、棒に作用する応力は引張方向の垂直応力 $\sigma$ だけでした。このような単軸応力（uniaxial stress）の場合、降伏の判定は単純に $\sigma \geq \sigma_Y$（降伏応力）かどうかを確認するだけで済みます。
しかし、実際の構造物や機械部品では、一つの方向だけに応力が作用するケースはむしろ例外です。典型的には以下のような応力が同時に作用します。
| 荷重の種類 | 生じる応力 | 具体例 |
|---|---|---|
| 引張・圧縮 | 垂直応力 $\sigma$ | ボルトの締結力 |
| 曲げ | 垂直応力 $\sigma$ | 曲げ応力を受ける梁 |
| ねじり | せん断応力 $\tau$ | ねじりを受ける軸 |
| 内圧 | 2方向の垂直応力 | 圧力容器の周方向・軸方向応力 |
こうした複数の応力成分が重なった状態が多軸応力（multiaxial stress）であり、「引張応力が50 MPa、せん断応力が30 MPa」というように複数の成分が同時に存在します。このとき「この部品は降伏するのか？」という問いに答えるためには、応力成分を統合的に評価する理論が必要です。それが組合せ応力の理論です。
応力テンソル — 多軸応力状態の完全な記述
3次元空間のある点における応力状態を完全に記述するには、応力テンソル（stress tensor）$\bm{\sigma}$ を用います。イメージとしては、その点を取り囲む微小な立方体の各面に作用する力の「カタログ」です。立方体の6つの面にはそれぞれ垂直応力とせん断応力が作用しますが、対面の応力は作用・反作用で同じなので、独立な面は3つ（$x$面、$y$面、$z$面）です。
応力テンソルは $3 \times 3$ の対称行列として表されます。
$$ \begin{equation} \bm{\sigma} = \begin{pmatrix} \sigma_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \sigma_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \sigma_{zz} \end{pmatrix} \end{equation} $$
ここで、対角成分 $\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{zz}$ は各面の垂直応力、非対角成分 $\tau_{xy}$, $\tau_{xz}$, $\tau_{yz}$ はせん断応力です。
角運動量の保存則（モーメントの釣り合い）から、せん断応力の相反性が成り立ちます。
$$ \begin{equation} \tau_{xy} = \tau_{yx}, \quad \tau_{xz} = \tau_{zx}, \quad \tau_{yz} = \tau_{zy} \end{equation} $$
したがって、9つの成分のうち独立なものは6つ（垂直応力3つ + せん断応力3つ）です。この6つの値が決まれば、その点の応力状態は完全に定まります。
ここで重要なのは、応力テンソルの値は座標系の取り方に依存するという点です。同じ応力状態でも、座標軸を回転させれば各成分の値は変わります。しかし、応力状態の「本質」— たとえば最大の引張応力がいくらか — は座標系に依存しません。この「座標系によらない本質」を取り出すのが、次に述べる応力の座標変換と主応力の理論です。
応力の座標変換
2次元の応力変換公式
まず、理解しやすい2次元（平面応力状態）から始めましょう。平面応力と平面ひずみで見たように、薄板の面内問題では $\sigma_z = \tau_{xz} = \tau_{yz} = 0$ とおけるため、応力テンソルは2次元に縮約されます。
$$ \begin{equation} \bm{\sigma}_{2D} = \begin{pmatrix} \sigma_x & \tau_{xy} \\ \tau_{xy} & \sigma_y \end{pmatrix} \end{equation} $$
$x$ 軸から反時計回りに角度 $\theta$ だけ回転した座標系 $(x’, y’)$ での応力を求めたいとします。この変換は、傾斜面での力の釣り合いから導出できます。主応力とモールの応力円で詳しく述べていますが、結果をまとめると次の通りです。
$$ \begin{equation} \sigma_{x’} = \frac{\sigma_x + \sigma_y}{2} + \frac{\sigma_x – \sigma_y}{2}\cos 2\theta + \tau_{xy}\sin 2\theta \end{equation} $$
$$ \begin{equation} \sigma_{y’} = \frac{\sigma_x + \sigma_y}{2} – \frac{\sigma_x – \sigma_y}{2}\cos 2\theta – \tau_{xy}\sin 2\theta \end{equation} $$
$$ \begin{equation} \tau_{x’y’} = -\frac{\sigma_x – \sigma_y}{2}\sin 2\theta + \tau_{xy}\cos 2\theta \end{equation} $$
この公式の構造を直感的に理解するために、いくつかの特徴を確認しましょう。
垂直応力の和は不変: $\sigma_{x’} + \sigma_{y’} = \sigma_x + \sigma_y$ が常に成り立ちます。座標系をどう回転させても、垂直応力の合計は変わりません。これは応力テンソルのトレース（対角成分の和）であり、テンソルの不変量の一つです。
倍角の三角関数: 物理的な回転角が $\theta$ であるのに対して、式中には $2\theta$ が現れます。これは、応力は方向ベクトルの2乗に比例する量（テンソル量）だからです。ベクトル成分の変換では $\cos\theta$ が現れますが、応力成分の変換では $\cos 2\theta$ が現れるのはこのためです。
回転行列による定式化
応力変換を行列形式で書くと、体系的で3次元への拡張も容易です。回転行列 $\bm{R}$ を
$$ \begin{equation} \bm{R} = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix} \end{equation} $$
とすると、回転後の応力テンソル $\bm{\sigma}’$ は次のように書けます。
$$ \begin{equation} \bm{\sigma}’ = \bm{R}\,\bm{\sigma}\,\bm{R}^T \end{equation} $$
これは線形代数の基礎で学ぶ相似変換（similarity transformation）の形をしています。応力テンソルの座標変換は、まさに行列の相似変換にほかなりません。
3次元の応力変換
3次元では、回転行列 $\bm{R}$ は $3 \times 3$ の直交行列（$\bm{R}^T\bm{R} = \bm{I}$）となります。変換の形式はまったく同じです。
$$ \begin{equation} \bm{\sigma}’ = \bm{R}\,\bm{\sigma}\,\bm{R}^T \end{equation} $$
ただし、3次元回転は3つのオイラー角（あるいは3つの回転軸）で記述されるため、一般的な変換公式を成分で書き下すと非常に長くなります。実用上は行列形式のまま扱い、数値計算に委ねるのが効率的です。
ここまでで、任意の方向の応力を計算する方法がわかりました。では次に、「どの方向を向けば応力が最大（または最小）になるのか？」という本質的な問いに進みましょう。それが主応力の理論です。
主応力と主方向
主応力とは — 「せん断応力がゼロになる面」を探す
ある特別な方向の面では、せん断応力が完全にゼロになり、面に垂直な応力（垂直応力）だけが作用します。このときの垂直応力を主応力（principal stress）、その面の法線方向を主方向（principal direction）と呼びます。
大雑把に言うと、主応力とは「座標系を最もうまく選んだときの応力」です。任意の座標系では応力テンソルの非対角成分（せん断応力）が現れますが、座標軸を主方向に合わせれば、応力テンソルは対角行列になります。
$$ \begin{equation} \bm{\sigma}_{\text{diag}} = \begin{pmatrix} \sigma_1 & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 & \sigma_3 \end{pmatrix} \end{equation} $$
ここで $\sigma_1 \geq \sigma_2 \geq \sigma_3$ と順に並べるのが慣例です。$\sigma_1$ を最大主応力、$\sigma_3$ を最小主応力と呼びます。
固有値問題としての定式化
主方向 $\bm{n}$ の面では、応力ベクトルが面に垂直な成分だけを持ちます。応力テンソル $\bm{\sigma}$ が面の法線ベクトル $\bm{n}$ に作用すると、結果の応力ベクトル $\bm{t}$ は
$$ \begin{equation} \bm{t} = \bm{\sigma}\,\bm{n} \end{equation} $$
です。主面では $\bm{t}$ が $\bm{n}$ に平行、つまり $\bm{t} = \sigma\,\bm{n}$ となります。したがって
$$ \begin{equation} \bm{\sigma}\,\bm{n} = \sigma\,\bm{n} \end{equation} $$
これを整理すると
$$ \begin{equation} (\bm{\sigma} – \sigma\,\bm{I})\,\bm{n} = \bm{0} \end{equation} $$
となり、これはまさに固有値・固有ベクトルの問題です。主応力 $\sigma$ は応力テンソルの固有値、主方向 $\bm{n}$ は固有ベクトルに対応します。
非自明な解 $\bm{n} \neq \bm{0}$ が存在する条件は、係数行列の行列式がゼロになることです。
$$ \begin{equation} \det(\bm{\sigma} – \sigma\,\bm{I}) = 0 \end{equation} $$
特性方程式と応力不変量
3次元の場合、この行列式を展開すると次の特性方程式（characteristic equation）が得られます。
$$ \begin{equation} \sigma^3 – I_1\,\sigma^2 + I_2\,\sigma – I_3 = 0 \end{equation} $$
ここで $I_1$, $I_2$, $I_3$ は応力不変量（stress invariants）と呼ばれ、座標系のとり方によらず一定の値をとります。
第1不変量（トレース）:
$$ \begin{equation} I_1 = \sigma_{xx} + \sigma_{yy} + \sigma_{zz} \end{equation} $$
第2不変量:
$$ \begin{equation} I_2 = \sigma_{xx}\sigma_{yy} + \sigma_{yy}\sigma_{zz} + \sigma_{zz}\sigma_{xx} – \tau_{xy}^2 – \tau_{yz}^2 – \tau_{zx}^2 \end{equation} $$
第3不変量（行列式）:
$$ \begin{equation} I_3 = \det(\bm{\sigma}) \end{equation} $$
これらの不変量は、座標系を回転させても値が変わりません。なぜなら、特性方程式の解（主応力）は応力状態の物理的な性質であり、座標系の選び方に依存しないからです。特性方程式の係数である $I_1$, $I_2$, $I_3$ も当然、座標系に依存しません。
2次元の場合の主応力公式
2次元（平面応力）では、特性方程式は2次方程式に帰着します。
$$ \begin{equation} \sigma^2 – (\sigma_x + \sigma_y)\,\sigma + (\sigma_x\sigma_y – \tau_{xy}^2) = 0 \end{equation} $$
解の公式を適用すると
$$ \begin{equation} \sigma_{1,2} = \frac{\sigma_x + \sigma_y}{2} \pm \sqrt{\left(\frac{\sigma_x – \sigma_y}{2}\right)^2 + \tau_{xy}^2} \end{equation} $$
を得ます。この式の構造を把握しましょう。
- 平均応力 $(\sigma_x + \sigma_y)/2$: 応力のベースライン（応力円の中心に対応）
- 偏差の二乗和の平方根: 応力の「広がり」（応力円の半径に対応）
主応力は平均応力を中心として、上下に同じだけずれた値になります。これはモールの応力円を見ると一目瞭然です。
主方向は次の式で求まります。
$$ \begin{equation} \tan 2\theta_p = \frac{2\tau_{xy}}{\sigma_x – \sigma_y} \end{equation} $$
最大せん断応力
主応力が求まれば、最大せん断応力も直ちに得られます。最大せん断応力は主方向から $45°$ 回転した面に作用し、その値は
2次元の場合:
$$ \begin{equation} \tau_{\max} = \frac{\sigma_1 – \sigma_2}{2} \end{equation} $$
3次元の場合は、3つの主応力間の差の最大値の半分です。
$$ \begin{equation} \tau_{\max} = \frac{\sigma_1 – \sigma_3}{2} \end{equation} $$
最大せん断応力は材料の降伏に深く関わります。金属の塑性変形は結晶面のすべり（せん断変形）によって生じるため、降伏条件にはせん断応力が本質的に関与するのです。
ここまでで主応力と最大せん断応力の理論がわかりました。これらを幾何学的に理解するための強力なツールが、次に述べるモールの応力円です。
モールの応力円
2次元モールの応力円
応力変換公式から円の方程式へ
応力変換公式を変形して、モールの応力円を導出しましょう。式 (4) と式 (6) を再掲します。
$$ \sigma_{x’} = \frac{\sigma_x + \sigma_y}{2} + \frac{\sigma_x – \sigma_y}{2}\cos 2\theta + \tau_{xy}\sin 2\theta $$
$$ \tau_{x’y’} = -\frac{\sigma_x – \sigma_y}{2}\sin 2\theta + \tau_{xy}\cos 2\theta $$
ここで $C = (\sigma_x + \sigma_y)/2$、$R_a = (\sigma_x – \sigma_y)/2$ とおきます。すると
$$ \sigma_{x’} – C = R_a\cos 2\theta + \tau_{xy}\sin 2\theta $$
$$ \tau_{x’y’} = -R_a\sin 2\theta + \tau_{xy}\cos 2\theta $$
両辺をそれぞれ2乗して加えると、$\cos^2 2\theta + \sin^2 2\theta = 1$ より
$$ \begin{equation} (\sigma_{x’} – C)^2 + \tau_{x’y’}^2 = R_a^2 + \tau_{xy}^2 = R^2 \end{equation} $$
ここで
$$ \begin{equation} R = \sqrt{\left(\frac{\sigma_x – \sigma_y}{2}\right)^2 + \tau_{xy}^2} \end{equation} $$
です。式 (22) は $(\sigma, \tau)$ 平面における円の方程式であり、これがモールの応力円です。
- 中心: $(C, 0) = \left(\dfrac{\sigma_x + \sigma_y}{2},\, 0\right)$
- 半径: $R = \sqrt{\left(\dfrac{\sigma_x – \sigma_y}{2}\right)^2 + \tau_{xy}^2}$
モールの応力円の読み方
モールの応力円を使うと、以下の情報を幾何学的に一目で読み取ることができます。
- 主応力: 円が $\sigma$ 軸と交わる点の $\sigma$ 座標が $\sigma_1$ と $\sigma_2$ です。つまり $\sigma_1 = C + R$、$\sigma_2 = C – R$ です
- 主方向: 現在の応力点 $(\sigma_x, \tau_{xy})$ から $\sigma$ 軸上の主応力点まで円周に沿って回る角度が $2\theta_p$ です。物理的な回転角はその半分の $\theta_p$ です
- 最大せん断応力: 円の最上部と最下部の $\tau$ 座標が $\pm\tau_{\max}$ です。$\tau_{\max} = R = (\sigma_1 – \sigma_2)/2$
- 任意方向の応力: 角度 $2\theta$ だけ円周に沿って回った点の座標が、$\theta$ 回転した面での $(\sigma_{x’}, \tau_{x’y’})$ です
モールの応力円の最大の利点は、応力変換の全体像を一枚の図で把握できることです。どの方向の面を見ても、応力は常にこの円の上の点に対応します。
注意: 符号の約束
モールの応力円では、せん断応力の符号について反時計回りのせん断を正とする約束が広く使われています（教科書によって異なる場合があります）。本記事では、$x$ 面に作用するせん断応力 $\tau_{xy}$ が正のとき、モールの応力円上の点を $(\sigma_x, \tau_{xy})$ にプロットします。このとき、$y$ 面の点は $(\sigma_y, -\tau_{xy})$ となり、2点を結ぶ直線が直径になります。
3次元モールの応力円
3次元の応力状態では、モールの応力円は3つの円で構成されます。3つの主応力 $\sigma_1 \geq \sigma_2 \geq \sigma_3$ に対して、以下の3つの円が描かれます。
- 外側の大円: 中心 $\dfrac{\sigma_1 + \sigma_3}{2}$、半径 $\dfrac{\sigma_1 – \sigma_3}{2}$
- 内側の左円: 中心 $\dfrac{\sigma_2 + \sigma_3}{2}$、半径 $\dfrac{\sigma_2 – \sigma_3}{2}$
- 内側の右円: 中心 $\dfrac{\sigma_1 + \sigma_2}{2}$、半径 $\dfrac{\sigma_1 – \sigma_2}{2}$
任意の面における応力 $(\sigma, \tau)$ は、外側の大円の内側、かつ2つの内側円の外側の領域に存在します。この領域は三日月形をしており、可能な応力状態が幾何学的に制約されることを示しています。
3次元の最大せん断応力は外側の大円の半径に等しく、$\tau_{\max} = (\sigma_1 – \sigma_3)/2$ です。これは最大主応力と最小主応力の差の半分であり、中間主応力 $\sigma_2$ は最大せん断応力に影響しないという重要な性質があります。
モールの応力円を使えば、応力状態の全体像が視覚的に理解できます。では次に、この応力状態が材料の降伏を引き起こすかどうかを判定する降伏条件に進みましょう。
Pythonによる応力変換と主応力の計算
まず、応力変換と主応力の計算をPythonで実装してみましょう。ここでは、与えられた2次元応力状態に対して、座標変換による応力の変化と主応力を計算します。
import numpy as np
def stress_transform_2d(sigma_x, sigma_y, tau_xy, theta_deg):
"""2次元応力変換: 角度theta(度)回転後の応力を計算"""
theta = np.radians(theta_deg)
c2 = np.cos(2 * theta)
s2 = np.sin(2 * theta)
sigma_avg = (sigma_x + sigma_y) / 2
sigma_diff = (sigma_x - sigma_y) / 2
sigma_xp = sigma_avg + sigma_diff * c2 + tau_xy * s2
sigma_yp = sigma_avg - sigma_diff * c2 - tau_xy * s2
tau_xpyp = -sigma_diff * s2 + tau_xy * c2
return sigma_xp, sigma_yp, tau_xpyp
def principal_stress_2d(sigma_x, sigma_y, tau_xy):
"""2次元の主応力・主方向・最大せん断応力を計算"""
sigma_avg = (sigma_x + sigma_y) / 2
R = np.sqrt(((sigma_x - sigma_y) / 2)**2 + tau_xy**2)
sigma_1 = sigma_avg + R
sigma_2 = sigma_avg - R
tau_max = R
# 主方向（度）
theta_p = 0.5 * np.degrees(np.arctan2(2 * tau_xy, sigma_x - sigma_y))
return sigma_1, sigma_2, tau_max, theta_p
def principal_stress_3d(stress_tensor):
"""3次元応力テンソルから主応力と主方向を計算（固有値問題）"""
eigenvalues, eigenvectors = np.linalg.eigh(stress_tensor)
# 降順にソート
idx = np.argsort(eigenvalues)[::-1]
sigma_principal = eigenvalues[idx]
directions = eigenvectors[:, idx]
return sigma_principal, directions
# --- 数値例 ---
# 2次元の例: 曲げとねじりの組合せ
sigma_x = 80.0 # MPa (曲げ応力)
sigma_y = -20.0 # MPa (横方向)
tau_xy = 40.0 # MPa (せん断応力)
sigma_1, sigma_2, tau_max, theta_p = principal_stress_2d(sigma_x, sigma_y, tau_xy)
print("=== 2次元の主応力 ===")
print(f"σ_x = {sigma_x} MPa, σ_y = {sigma_y} MPa, τ_xy = {tau_xy} MPa")
print(f"σ₁ = {sigma_1:.2f} MPa")
print(f"σ₂ = {sigma_2:.2f} MPa")
print(f"τ_max = {tau_max:.2f} MPa")
print(f"主方向 θ_p = {theta_p:.2f}°")
# 3次元の例
stress_3d = np.array([
[100, 30, 0],
[30, -20, 15],
[0, 15, 50]
], dtype=float)
sigma_p, directions = principal_stress_3d(stress_3d)
print("\n=== 3次元の主応力 ===")
print(f"応力テンソル:\n{stress_3d}")
print(f"σ₁ = {sigma_p[0]:.2f} MPa")
print(f"σ₂ = {sigma_p[1]:.2f} MPa")
print(f"σ₃ = {sigma_p[2]:.2f} MPa")
print(f"τ_max = {(sigma_p[0] - sigma_p[2]) / 2:.2f} MPa")
このコードでは、2次元と3次元の両方で主応力を計算しています。2次元では解析解（公式）を直接使い、3次元では np.linalg.eigh を用いて応力テンソルの固有値を数値的に求めています。eigh は対称行列専用の関数であり、応力テンソルの対称性を活用して数値的に安定かつ効率的に固有値を計算します。
出力結果を見ると、2次元の例では $\sigma_x = 80$ MPa、$\sigma_y = -20$ MPa、$\tau_{xy} = 40$ MPa という応力状態に対して、最大主応力が約 96.1 MPa、最小主応力が約 $-36.1$ MPa と求まります。元の垂直応力の最大値（80 MPa）よりも主応力の最大値の方が大きいことに注目してください。せん断応力の存在によって、特定の方向では垂直応力がさらに大きくなるのです。これが組合せ応力を考慮しなければならない本質的な理由です。
モールの応力円のPython描画
モールの応力円を描画することで、応力変換の全体像を視覚的に把握しましょう。2次元と3次元の両方を描画します。
import numpy as np
import matplotlib.pyplot as plt
def plot_mohr_circle_2d(sigma_x, sigma_y, tau_xy, ax=None):
"""2次元モールの応力円を描画"""
if ax is None:
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
# 中心と半径
C = (sigma_x + sigma_y) / 2
R = np.sqrt(((sigma_x - sigma_y) / 2)**2 + tau_xy**2)
# 主応力
sigma_1 = C + R
sigma_2 = C - R
# 応力円を描画
theta_circle = np.linspace(0, 2 * np.pi, 200)
s_circle = C + R * np.cos(theta_circle)
t_circle = R * np.sin(theta_circle)
ax.plot(s_circle, t_circle, 'b-', linewidth=2, label="Mohr's circle")
# 現在の応力点
ax.plot(sigma_x, tau_xy, 'ro', markersize=10,
label=f'X face ($\\sigma_x$={sigma_x}, $\\tau$={tau_xy})')
ax.plot(sigma_y, -tau_xy, 'gs', markersize=10,
label=f'Y face ($\\sigma_y$={sigma_y}, $\\tau$={-tau_xy})')
# 直径の線
ax.plot([sigma_x, sigma_y], [tau_xy, -tau_xy], 'k--', linewidth=1)
# 主応力点
ax.plot(sigma_1, 0, 'r^', markersize=12,
label=f'$\\sigma_1$ = {sigma_1:.1f} MPa')
ax.plot(sigma_2, 0, 'bv', markersize=12,
label=f'$\\sigma_2$ = {sigma_2:.1f} MPa')
# 最大せん断応力
ax.plot(C, R, 'mp', markersize=10,
label=f'$\\tau_{{max}}$ = {R:.1f} MPa')
ax.plot(C, -R, 'mp', markersize=10)
# 中心
ax.plot(C, 0, 'k+', markersize=12, markeredgewidth=2)
# 軸と装飾
ax.axhline(y=0, color='k', linewidth=0.5)
ax.axvline(x=0, color='k', linewidth=0.5)
ax.set_xlabel('$\\sigma$ [MPa]', fontsize=14)
ax.set_ylabel('$\\tau$ [MPa]', fontsize=14)
ax.set_title("Mohr's Circle (2D)", fontsize=16)
ax.set_aspect('equal')
ax.legend(fontsize=10, loc='upper left')
ax.grid(True, alpha=0.3)
return ax
def plot_mohr_circle_3d(sigma_1, sigma_2, sigma_3, ax=None):
"""3次元モールの応力円（3つの円）を描画"""
if ax is None:
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
theta = np.linspace(0, 2 * np.pi, 200)
# 3つの円
circles = [
((sigma_1 + sigma_3) / 2, (sigma_1 - sigma_3) / 2, 'b', 'Major circle'),
((sigma_1 + sigma_2) / 2, (sigma_1 - sigma_2) / 2, 'r', 'Circle 1-2'),
((sigma_2 + sigma_3) / 2, (sigma_2 - sigma_3) / 2, 'g', 'Circle 2-3'),
]
for center, radius, color, label in circles:
s = center + radius * np.cos(theta)
t = radius * np.sin(theta)
ax.plot(s, t, color=color, linewidth=2, label=label)
# 主応力点
for s, name in [(sigma_1, '$\\sigma_1$'), (sigma_2, '$\\sigma_2$'),
(sigma_3, '$\\sigma_3$')]:
ax.plot(s, 0, 'ko', markersize=8)
ax.annotate(f'{name}={s:.1f}', (s, 0), textcoords="offset points",
xytext=(0, 15), ha='center', fontsize=11)
ax.axhline(y=0, color='k', linewidth=0.5)
ax.axvline(x=0, color='k', linewidth=0.5)
ax.set_xlabel('$\\sigma$ [MPa]', fontsize=14)
ax.set_ylabel('$\\tau$ [MPa]', fontsize=14)
ax.set_title("Mohr's Circles (3D)", fontsize=16)
ax.set_aspect('equal')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
return ax
# --- 描画 ---
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
# 2Dモールの応力円
plot_mohr_circle_2d(80, -20, 40, ax=axes[0])
# 3Dモールの応力円
plot_mohr_circle_3d(110.0, 52.0, -32.0, ax=axes[1])
plt.tight_layout()
plt.savefig("mohr_circles.png", dpi=150, bbox_inches='tight')
plt.show()
左のグラフ（2次元モールの応力円）を見ると、以下の特徴が読み取れます。
- 応力円の中心は $\sigma$ 軸上の $(\sigma_x + \sigma_y)/2 = 30$ MPa の位置にあります。これは垂直応力の平均値であり、座標変換によって変わりません
- 赤い点（X面）と緑の点（Y面）は直径の両端に位置しています。この2点を結ぶ直線が直径であることは、2面のなす角が $90°$（物理的角度）であり、モールの応力円上では $180°$（$2\theta$ の角度）に対応することを反映しています
- 主応力点（$\sigma$ 軸との交点）では $\tau = 0$ となり、せん断応力がゼロの面が実現していることが確認できます
- 最大せん断応力は円の最上部・最下部に対応し、その値は半径 $R$ に等しくなっています
右のグラフ（3次元モールの応力円）では、3つの主応力 $\sigma_1 = 110$, $\sigma_2 = 52$, $\sigma_3 = -32$ MPa に対応する3つの円が描かれています。最大せん断応力は外側の大円の半径 $(110 – (-32))/2 = 71$ MPa であり、中間主応力 $\sigma_2$ には依存しないことが視覚的に確認できます。
降伏条件
なぜ降伏条件が必要なのか
単軸引張であれば、降伏の判定は $\sigma \geq \sigma_Y$ という単純な比較で済みます。しかし、組合せ応力状態では $\sigma_x$, $\sigma_y$, $\tau_{xy}$（あるいは3次元なら6つの応力成分）が同時に存在するため、「どの値をどう組み合わせたら降伏するのか」を定めるルールが必要です。これが降伏条件（yield criterion）です。
降伏条件は、6次元の応力空間において降伏が始まる境界面（降伏曲面）を定義します。2次元の場合は $\sigma_1$-$\sigma_2$ 平面上の降伏曲線として描けます。降伏曲線の内側は弾性領域、外側は塑性領域（降伏が生じる領域）です。
以下で、代表的な2つの降伏条件を見ていきましょう。
トレスカの降伏条件（最大せん断応力説）
直感
金属材料の降伏メカニズムは、結晶格子面のすべり（転位の運動）です。すべりはせん断応力によって駆動されます。したがって、「最大せん断応力がある臨界値に達したとき降伏が始まる」と考えるのが自然です。これがトレスカ（Tresca）の降伏条件の基本的なアイデアです。
定式化
単軸引張で降伏応力 $\sigma_Y$ に達したとき、最大せん断応力は $\tau_{\max} = \sigma_Y / 2$ です（モールの応力円から即座にわかります）。トレスカの条件は、この $\sigma_Y / 2$ を臨界値として、一般の多軸応力状態に拡張します。
$$ \begin{equation} \tau_{\max} = \frac{\sigma_1 – \sigma_3}{2} \leq \frac{\sigma_Y}{2} \end{equation} $$
これは次のように書き換えられます。
$$ \begin{equation} \sigma_1 – \sigma_3 \leq \sigma_Y \end{equation} $$
2次元（$\sigma_3 = 0$ の場合）では、$\sigma_1$-$\sigma_2$ 平面上の降伏曲線は六角形になります。
$$ \begin{equation} \max(|\sigma_1|, |\sigma_2|, |\sigma_1 – \sigma_2|) \leq \sigma_Y \end{equation} $$
この六角形の各辺は、3つの主応力差のうちの一つが $\sigma_Y$ に達する条件に対応しています。
トレスカ条件の長所と短所
長所: 物理的に明快で直感的。最大せん断応力だけを見ればよいので計算が単純。
短所: 中間主応力 $\sigma_2$ が降伏に影響しないとする仮定は、実験結果と完全には一致しません。また、六角形の頂点で微分不可能であるため、塑性理論での数値計算（塑性流れ則の適用）にやや不便です。
フォンミーゼスの降伏条件（せん断ひずみエネルギー説）
直感
弾性体に蓄えられるひずみエネルギーは、体積変化に関するエネルギー（均圧成分）と形状変化に関するエネルギー（偏差成分）に分解できます。静水圧は体積を均等に圧縮するだけで、すべりを引き起こしません。降伏を引き起こすのは形状を歪めるエネルギー（偏差ひずみエネルギー）だけです。
フォンミーゼス（von Mises）の降伏条件は、「偏差ひずみエネルギーが臨界値に達したとき降伏する」と仮定します。
ひずみエネルギーの分解
単位体積あたりのひずみエネルギーはフックの法則と弾性係数から
$$ \begin{equation} U = \frac{1}{2E}\left[\sigma_1^2 + \sigma_2^2 + \sigma_3^2 – 2\nu(\sigma_1\sigma_2 + \sigma_2\sigma_3 + \sigma_3\sigma_1)\right] \end{equation} $$
です。ここで 静水圧応力（平均応力）$\sigma_m$ を
$$ \begin{equation} \sigma_m = \frac{\sigma_1 + \sigma_2 + \sigma_3}{3} \end{equation} $$
と定義します。体積変化に関するエネルギー $U_v$ は $\sigma_1 = \sigma_2 = \sigma_3 = \sigma_m$ としたときのひずみエネルギーに等しく
$$ \begin{equation} U_v = \frac{3(1 – 2\nu)}{2E}\sigma_m^2 \end{equation} $$
形状変化に関するエネルギー（偏差ひずみエネルギー）$U_d$ は $U – U_v$ として
$$ \begin{equation} U_d = \frac{1 + \nu}{6E}\left[(\sigma_1 – \sigma_2)^2 + (\sigma_2 – \sigma_3)^2 + (\sigma_3 – \sigma_1)^2\right] \end{equation} $$
と求まります。
フォンミーゼス応力の定義
単軸引張（$\sigma_1 = \sigma_Y$, $\sigma_2 = \sigma_3 = 0$）のときの偏差ひずみエネルギーは
$$ U_d^{(Y)} = \frac{1 + \nu}{6E}\left[\sigma_Y^2 + \sigma_Y^2\right] = \frac{1 + \nu}{3E}\sigma_Y^2 $$
一般の多軸応力状態で $U_d = U_d^{(Y)}$ となる条件は
$$ (\sigma_1 – \sigma_2)^2 + (\sigma_2 – \sigma_3)^2 + (\sigma_3 – \sigma_1)^2 = 2\sigma_Y^2 $$
ここで、フォンミーゼス応力（von Mises stress, 相当応力）$\sigma_v$ を次のように定義します。
$$ \begin{equation} \sigma_v = \sqrt{\frac{(\sigma_1 – \sigma_2)^2 + (\sigma_2 – \sigma_3)^2 + (\sigma_3 – \sigma_1)^2}{2}} \end{equation} $$
フォンミーゼスの降伏条件は
$$ \begin{equation} \sigma_v \leq \sigma_Y \end{equation} $$
と簡潔に書けます。つまり、フォンミーゼス応力は多軸応力状態を単一の等価な引張応力に換算したものです。
応力成分による表現
主応力ではなく、元の応力成分で表すと
$$ \begin{equation} \sigma_v = \sqrt{\frac{(\sigma_{xx} – \sigma_{yy})^2 + (\sigma_{yy} – \sigma_{zz})^2 + (\sigma_{zz} – \sigma_{xx})^2 + 6(\tau_{xy}^2 + \tau_{yz}^2 + \tau_{zx}^2)}{2}} \end{equation} $$
この表現は、主応力を求めなくても直接応力成分からフォンミーゼス応力を計算できるため、FEM（有限要素法）のポスト処理で頻繁に使われます。
2次元の降伏曲線
2次元（$\sigma_3 = 0$）では
$$ \begin{equation} \sigma_v = \sqrt{\sigma_1^2 – \sigma_1\sigma_2 + \sigma_2^2} \leq \sigma_Y \end{equation} $$
この降伏曲線は $\sigma_1$-$\sigma_2$ 平面上の楕円です。具体的には、$\sigma_1$ 軸と $\sigma_2$ 軸に対して $45°$ 傾いた楕円であり、トレスカの六角形に外接します。
フォンミーゼス条件の長所と短所
長所: 中間主応力の影響を自然に取り込んでおり、多くの延性金属の実験結果とよく一致します。降伏曲面が滑らか（微分可能）であるため、塑性流れ則の適用が容易で、数値解析（FEM）に適しています。
短所: 静水圧に対する感度がないため、脆性材料（コンクリート、岩石など）や圧力依存型の材料には適用できません。これらの材料にはモール・クーロンの条件やドラッカー・プラガーの条件が使われます。
トレスカとフォンミーゼスの比較
| 項目 | トレスカ | フォンミーゼス |
|---|---|---|
| 降伏判定量 | $\tau_{\max} = (\sigma_1 – \sigma_3)/2$ | $\sigma_v = \sqrt{((\sigma_1-\sigma_2)^2 + \cdots)/2}$ |
| 2D降伏曲線 | 六角形 | 楕円 |
| 中間主応力の影響 | なし | あり |
| 実験との一致 | やや保守的 | 良好 |
| 数値計算 | 頂点で非平滑 | 平滑で扱いやすい |
| 安全側 | 常にフォンミーゼスより保守的 | トレスカよりわずかに非保守的 |
トレスカの六角形はフォンミーゼスの楕円に内接します。したがって、トレスカ条件で安全と判定された設計は、フォンミーゼス条件でも必ず安全です。両者の最大差は約 $15.5\%$（$\sigma_1 = -\sigma_2$ の場合、すなわち純粋せん断の場合）です。
では、これらの降伏条件をPythonで可視化して、両者の違いを幾何学的に確認しましょう。
降伏条件の比較: Pythonによる可視化
トレスカとフォンミーゼスの降伏曲線を $\sigma_1$-$\sigma_2$ 平面上に描画し、さらにフォンミーゼス応力の等高線図を重ね合わせます。
import numpy as np
import matplotlib.pyplot as plt
sigma_Y = 250.0 # MPa (降伏応力)
# --- フォンミーゼス降伏曲線 (楕円) ---
theta = np.linspace(0, 2 * np.pi, 500)
# パラメトリック表現: σ₁² - σ₁σ₂ + σ₂² = σ_Y²
# 45°回転した座標系で楕円の標準形になる
s1_vm = sigma_Y * (np.cos(theta) * np.cos(np.pi/6)
- np.sin(theta) * np.sin(np.pi/6))
s2_vm = sigma_Y * (np.cos(theta) * np.sin(np.pi/6)
+ np.sin(theta) * np.cos(np.pi/6))
# 正確な表現を使用
t_param = np.linspace(0, 2 * np.pi, 500)
s1_vm = sigma_Y * np.sqrt(2/3) * np.cos(t_param) + sigma_Y / 3 * np.sqrt(2) * np.sin(t_param)
# 直接的にパラメトリックに生成
phi = np.linspace(0, 2 * np.pi, 500)
# σ₁ = σ_Y * (cos(phi) + sin(phi)/√3)
# σ₂ = σ_Y * (-cos(phi) + sin(phi)/√3) は不正確
# 代わりに等高線から抽出する方法を使う
# フォンミーゼス: σ₁² - σ₁σ₂ + σ₂² = σ_Y² を等高線で描く
s1_grid = np.linspace(-1.3 * sigma_Y, 1.3 * sigma_Y, 500)
s2_grid = np.linspace(-1.3 * sigma_Y, 1.3 * sigma_Y, 500)
S1, S2 = np.meshgrid(s1_grid, s2_grid)
VM = np.sqrt(S1**2 - S1 * S2 + S2**2)
# --- トレスカ降伏曲線 (六角形) ---
# σ₃ = 0 の場合の6つの条件:
# |σ₁| ≤ σ_Y, |σ₂| ≤ σ_Y, |σ₁ - σ₂| ≤ σ_Y
tresca_vertices = np.array([
[sigma_Y, 0],
[sigma_Y, sigma_Y],
[0, sigma_Y],
[-sigma_Y, 0],
[-sigma_Y, -sigma_Y],
[0, -sigma_Y],
[sigma_Y, 0] # 閉じる
])
# --- 描画 ---
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
# (1) 降伏曲線の比較
ax1 = axes[0]
# フォンミーゼス
ax1.contour(S1, S2, VM, levels=[sigma_Y], colors='blue', linewidths=2.5)
# トレスカ
ax1.plot(tresca_vertices[:, 0], tresca_vertices[:, 1],
'r-', linewidth=2.5, label='Tresca')
ax1.plot([], [], 'b-', linewidth=2.5, label='von Mises')
# 装飾
ax1.axhline(y=0, color='k', linewidth=0.5)
ax1.axvline(x=0, color='k', linewidth=0.5)
ax1.set_xlabel('$\\sigma_1$ [MPa]', fontsize=14)
ax1.set_ylabel('$\\sigma_2$ [MPa]', fontsize=14)
ax1.set_title('Yield Loci: Tresca vs von Mises', fontsize=14)
ax1.set_aspect('equal')
ax1.legend(fontsize=12)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-350, 350)
ax1.set_ylim(-350, 350)
# 応力状態の例をプロット
test_points = [
(200, 50, 'A: Safe'),
(200, -100, 'B: Between'),
(150, 150, 'C: Safe'),
]
for s1, s2, label in test_points:
ax1.plot(s1, s2, 'ko', markersize=8)
ax1.annotate(label, (s1, s2), textcoords="offset points",
xytext=(10, 10), fontsize=10)
# (2) フォンミーゼス応力の等高線図
ax2 = axes[1]
levels = np.linspace(0, 1.5 * sigma_Y, 15)
contour = ax2.contourf(S1, S2, VM, levels=levels, cmap='RdYlGn_r')
ax2.contour(S1, S2, VM, levels=[sigma_Y], colors='white',
linewidths=3, linestyles='--')
fig.colorbar(contour, ax=ax2, label='$\\sigma_v$ [MPa]')
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.axvline(x=0, color='k', linewidth=0.5)
ax2.set_xlabel('$\\sigma_1$ [MPa]', fontsize=14)
ax2.set_ylabel('$\\sigma_2$ [MPa]', fontsize=14)
ax2.set_title('von Mises Stress Contours', fontsize=14)
ax2.set_aspect('equal')
ax2.set_xlim(-350, 350)
ax2.set_ylim(-350, 350)
plt.tight_layout()
plt.savefig("yield_criteria_comparison.png", dpi=150, bbox_inches='tight')
plt.show()
左のグラフから、トレスカとフォンミーゼスの降伏曲線の関係が明確に読み取れます。
- トレスカの六角形はフォンミーゼスの楕円に内接しています。$\sigma_1$ 軸上と $\sigma_2$ 軸上（単軸引張/圧縮に対応する点）では両者は一致します。これは当然です — どちらの条件も、単軸引張での降伏応力 $\sigma_Y$ を基準に定義されているからです
- 最も差が大きいのは純粋せん断の状態（$\sigma_1 = -\sigma_2$）のときです。この方向では、トレスカの降伏が $\sigma_1 = \sigma_Y / 2$ で生じるのに対し、フォンミーゼスでは $\sigma_1 = \sigma_Y / \sqrt{3} \approx 0.577 \sigma_Y$ で降伏します。差は約 $15.5\%$ です
- 点 A（200, 50）は両方の降伏曲線の内側にあり安全、点 B（200, $-100$）はトレスカでは降伏しますがフォンミーゼスでは安全領域にある可能性がある応力状態を示しています
右のグラフ（フォンミーゼス応力の等高線図）では、等高線が楕円形をしており、$\sigma_1 = \sigma_2$ の方向（等二軸引張）で最もフォンミーゼス応力が低く、$\sigma_1 = -\sigma_2$ の方向（純粋せん断）で最も高くなることが確認できます。白い破線が降伏曲線 $\sigma_v = \sigma_Y$ に対応しています。
設計への応用: 安全率の計算
安全率の定義
組合せ応力状態における安全率（factor of safety）$n$ は、降伏条件に基づいて次のように定義されます。
フォンミーゼス基準の場合:
$$ \begin{equation} n = \frac{\sigma_Y}{\sigma_v} \end{equation} $$
トレスカ基準の場合:
$$ \begin{equation} n = \frac{\sigma_Y}{\sigma_1 – \sigma_3} \end{equation} $$
$n > 1$ であれば弾性域内（安全）、$n = 1$ で降伏、$n < 1$ は既に塑性域（設計として不適切）です。
実際の設計では、荷重の不確実性、材料特性のばらつき、環境劣化などを考慮して $n \geq 1.5 \sim 3.0$ 程度の安全率を確保することが一般的です。航空宇宙分野では1.5（究極荷重に対して）、一般機械では2.0〜3.0が目安です。
設計例: 曲げとねじりを受ける軸
具体的な設計例として、曲げ応力とねじりが同時に作用する中実丸軸を考えましょう。
軸の表面（最も厳しい点）では、次の応力が作用します。
- 曲げによる垂直応力: $\sigma_b = \dfrac{32M}{\pi d^3}$
- ねじりによるせん断応力: $\tau_t = \dfrac{16T}{\pi d^3}$
ここで $M$ は曲げモーメント、$T$ はねじりモーメント、$d$ は軸の直径です。
応力状態は $\sigma_x = \sigma_b$, $\sigma_y = 0$, $\tau_{xy} = \tau_t$ なので、2次元の主応力は
$$ \sigma_{1,2} = \frac{\sigma_b}{2} \pm \sqrt{\left(\frac{\sigma_b}{2}\right)^2 + \tau_t^2} $$
フォンミーゼス応力は
$$ \begin{equation} \sigma_v = \sqrt{\sigma_b^2 + 3\tau_t^2} \end{equation} $$
ここで、2次元で $\sigma_y = 0$ の場合のフォンミーゼス応力を式 (35) から導きましょう。$\sigma_1$, $\sigma_2$ を代入して計算することもできますが、より直接的には式 (36) で $\sigma_{xx} = \sigma_b$, $\sigma_{yy} = 0$, $\sigma_{zz} = 0$, $\tau_{xy} = \tau_t$, 他のせん断成分 = 0 を代入すると
$$ \sigma_v = \sqrt{\sigma_b^2 + 0 + 0 + 6\tau_t^2 + 0 + 0} / \sqrt{2} $$
ではなく、正しくは
$$ \sigma_v = \sqrt{\frac{(\sigma_b – 0)^2 + (0 – 0)^2 + (0 – \sigma_b)^2 + 6\tau_t^2}{2}} = \sqrt{\frac{2\sigma_b^2 + 6\tau_t^2}{2}} = \sqrt{\sigma_b^2 + 3\tau_t^2} $$
となります。この $\sqrt{\sigma_b^2 + 3\tau_t^2}$ という形は、曲げとねじりを受ける軸の設計で頻繁に使われる公式です。
重要な点は、もし曲げとねじりを個別に評価して $\sigma_b < \sigma_Y$ かつ $\tau_t < \sigma_Y / \sqrt{3}$ を確認するだけでは不十分だということです。組合せ応力として $\sigma_v = \sqrt{\sigma_b^2 + 3\tau_t^2} \leq \sigma_Y$ を確認しなければ、安全性は保証されません。
では、この設計計算をPythonで実装してみましょう。
実践的な設計計算のPython実装
曲げとねじりを受ける軸の設計計算を行い、安全率を求めます。さらに、$M$-$T$ 平面上で安全率の等高線図を描き、設計点がどの位置にあるかを視覚的に確認します。
import numpy as np
import matplotlib.pyplot as plt
def shaft_design(M, T, d, sigma_Y):
"""
曲げモーメントMとねじりモーメントTを受ける
直径dの中実丸軸の応力と安全率を計算
Parameters
----------
M : float 曲げモーメント [N·m]
T : float ねじりモーメント [N·m]
d : float 軸の直径 [m]
sigma_Y : float 降伏応力 [Pa]
Returns
-------
dict 各応力と安全率
"""
sigma_b = 32 * M / (np.pi * d**3) # 曲げ応力 [Pa]
tau_t = 16 * T / (np.pi * d**3) # ねじり応力 [Pa]
# 主応力
sigma_1 = sigma_b / 2 + np.sqrt((sigma_b / 2)**2 + tau_t**2)
sigma_2 = sigma_b / 2 - np.sqrt((sigma_b / 2)**2 + tau_t**2)
# フォンミーゼス応力
sigma_vm = np.sqrt(sigma_b**2 + 3 * tau_t**2)
# トレスカ応力 (σ₃ = 0, σ₁ > 0 > σ₂ の場合)
sigma_tresca = sigma_1 - sigma_2 # = 2√((σ_b/2)² + τ_t²)
# 安全率
n_vm = sigma_Y / sigma_vm
n_tresca = sigma_Y / sigma_tresca
return {
'sigma_b': sigma_b, 'tau_t': tau_t,
'sigma_1': sigma_1, 'sigma_2': sigma_2,
'sigma_vm': sigma_vm, 'sigma_tresca': sigma_tresca,
'n_vm': n_vm, 'n_tresca': n_tresca
}
# --- 設計例 ---
M = 500 # N·m
T = 300 # N·m
d = 0.04 # 40 mm
sigma_Y = 250e6 # 250 MPa (S45C鋼)
result = shaft_design(M, T, d, sigma_Y)
print("=== 軸の設計計算 ===")
print(f"曲げモーメント M = {M} N·m")
print(f"ねじりモーメント T = {T} N·m")
print(f"軸直径 d = {d*1000:.0f} mm")
print(f"降伏応力 σ_Y = {sigma_Y/1e6:.0f} MPa")
print()
print(f"曲げ応力 σ_b = {result['sigma_b']/1e6:.2f} MPa")
print(f"ねじり応力 τ_t = {result['tau_t']/1e6:.2f} MPa")
print(f"主応力 σ₁ = {result['sigma_1']/1e6:.2f} MPa")
print(f"主応力 σ₂ = {result['sigma_2']/1e6:.2f} MPa")
print(f"フォンミーゼス応力 σ_v = {result['sigma_vm']/1e6:.2f} MPa")
print(f"トレスカ応力 σ₁−σ₂ = {result['sigma_tresca']/1e6:.2f} MPa")
print(f"安全率 (von Mises) n = {result['n_vm']:.3f}")
print(f"安全率 (Tresca) n = {result['n_tresca']:.3f}")
# --- M-T平面上の安全率等高線図 ---
M_range = np.linspace(0, 1000, 300)
T_range = np.linspace(0, 800, 300)
M_grid, T_grid = np.meshgrid(M_range, T_range)
sigma_b_grid = 32 * M_grid / (np.pi * d**3)
tau_t_grid = 16 * T_grid / (np.pi * d**3)
sigma_vm_grid = np.sqrt(sigma_b_grid**2 + 3 * tau_t_grid**2)
n_vm_grid = sigma_Y / sigma_vm_grid
fig, ax = plt.subplots(figsize=(9, 7))
levels = [0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 2.5, 3.0, 4.0]
contour = ax.contourf(M_range, T_range, n_vm_grid, levels=levels,
cmap='RdYlGn', extend='both')
ax.contour(M_range, T_range, n_vm_grid, levels=[1.0],
colors='white', linewidths=3, linestyles='--')
ax.contour(M_range, T_range, n_vm_grid, levels=[1.5],
colors='yellow', linewidths=2, linestyles='-.')
fig.colorbar(contour, ax=ax, label='Safety factor $n$ (von Mises)')
# 設計点
ax.plot(M, T, 'k*', markersize=18, markeredgewidth=1.5,
label=f'Design point (n={result["n_vm"]:.2f})')
ax.set_xlabel('Bending moment $M$ [N$\\cdot$m]', fontsize=14)
ax.set_ylabel('Torsional moment $T$ [N$\\cdot$m]', fontsize=14)
ax.set_title(f'Safety Factor Map ($d$={d*1000:.0f}mm, '
f'$\\sigma_Y$={sigma_Y/1e6:.0f}MPa)', fontsize=14)
ax.legend(fontsize=12, loc='upper right')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("safety_factor_map.png", dpi=150, bbox_inches='tight')
plt.show()
この安全率マップから、以下の設計上の知見が読み取れます。
- 白い破線が $n = 1$（降伏限界）の等高線であり、この曲線の内側（左下）が弾性域、外側が塑性域です。設計点はこの曲線の内側になければなりません
- 黄色い一点鎖線が $n = 1.5$ の等高線であり、一般的な設計基準ではこの内側に設計点を置くことが推奨されます
- 等高線が楕円形をしているのは、フォンミーゼス応力 $\sigma_v = \sqrt{\sigma_b^2 + 3\tau_t^2}$ の等高線が楕円であることを反映しています。ねじりモーメント $T$ の影響は $\sqrt{3}$ 倍に増幅されるため、同じ大きさの $M$ と $T$ ではねじりの方が降伏に対する寄与が大きくなります
- 設計点（黒い星印）の安全率を確認し、$n \geq 1.5$ を満たしていない場合は軸径 $d$ を太くするか材料を変更する必要があります
組合せ応力計算の総合的なPython実装
最後に、これまでの内容を統合した総合的な計算・可視化コードを示します。任意の3次元応力テンソルを入力として、主応力、フォンミーゼス応力、安全率、3次元モールの応力円を一括で計算・描画します。
import numpy as np
import matplotlib.pyplot as plt
class CombinedStressAnalyzer:
"""組合せ応力の総合解析クラス"""
def __init__(self, stress_tensor, sigma_Y):
"""
Parameters
----------
stress_tensor : ndarray, shape (3, 3)
応力テンソル [MPa]
sigma_Y : float
降伏応力 [MPa]
"""
self.sigma = np.array(stress_tensor, dtype=float)
self.sigma_Y = sigma_Y
self._compute()
def _compute(self):
"""主応力・フォンミーゼス応力・安全率を計算"""
# 主応力と主方向（固有値問題）
eigenvalues, eigenvectors = np.linalg.eigh(self.sigma)
idx = np.argsort(eigenvalues)[::-1]
self.principal_stresses = eigenvalues[idx]
self.principal_directions = eigenvectors[:, idx]
s1, s2, s3 = self.principal_stresses
# 最大せん断応力
self.tau_max = (s1 - s3) / 2
# フォンミーゼス応力
self.sigma_vm = np.sqrt(
((s1 - s2)**2 + (s2 - s3)**2 + (s3 - s1)**2) / 2
)
# 応力不変量
self.I1 = np.trace(self.sigma)
self.I2 = (self.sigma[0,0]*self.sigma[1,1]
+ self.sigma[1,1]*self.sigma[2,2]
+ self.sigma[2,2]*self.sigma[0,0]
- self.sigma[0,1]**2
- self.sigma[1,2]**2
- self.sigma[2,0]**2)
self.I3 = np.linalg.det(self.sigma)
# 安全率
self.n_vm = self.sigma_Y / self.sigma_vm
self.n_tresca = self.sigma_Y / (s1 - s3)
def summary(self):
"""解析結果をコンソールに出力"""
s = self.principal_stresses
print("=" * 50)
print("組合せ応力の解析結果")
print("=" * 50)
print(f"\n応力テンソル [MPa]:")
print(self.sigma)
print(f"\n応力不変量:")
print(f" I₁ = {self.I1:.2f} MPa")
print(f" I₂ = {self.I2:.2f} MPa²")
print(f" I₃ = {self.I3:.2f} MPa³")
print(f"\n主応力:")
print(f" σ₁ = {s[0]:.2f} MPa")
print(f" σ₂ = {s[1]:.2f} MPa")
print(f" σ₃ = {s[2]:.2f} MPa")
print(f"\n最大せん断応力: τ_max = {self.tau_max:.2f} MPa")
print(f"フォンミーゼス応力: σ_v = {self.sigma_vm:.2f} MPa")
print(f"\n安全率:")
print(f" von Mises: n = {self.n_vm:.3f}")
print(f" Tresca: n = {self.n_tresca:.3f}")
print(f"\n降伏判定 (σ_Y = {self.sigma_Y:.0f} MPa):")
status_vm = "安全" if self.n_vm >= 1.0 else "降伏"
status_tr = "安全" if self.n_tresca >= 1.0 else "降伏"
print(f" von Mises: {status_vm} (n = {self.n_vm:.3f})")
print(f" Tresca: {status_tr} (n = {self.n_tresca:.3f})")
def plot_mohr_3d(self, ax=None):
"""3次元モールの応力円を描画"""
if ax is None:
fig, ax = plt.subplots(figsize=(8, 8))
s1, s2, s3 = self.principal_stresses
theta = np.linspace(0, 2 * np.pi, 300)
circles_data = [
((s1+s3)/2, (s1-s3)/2, 'C0', f'$\\sigma_1$-$\\sigma_3$'),
((s1+s2)/2, (s1-s2)/2, 'C1', f'$\\sigma_1$-$\\sigma_2$'),
((s2+s3)/2, (s2-s3)/2, 'C2', f'$\\sigma_2$-$\\sigma_3$'),
]
for center, radius, color, label in circles_data:
ax.plot(center + radius * np.cos(theta),
radius * np.sin(theta),
color=color, linewidth=2, label=label)
for s, name in [(s1,'$\\sigma_1$'), (s2,'$\\sigma_2$'),
(s3,'$\\sigma_3$')]:
ax.plot(s, 0, 'ko', ms=8)
ax.annotate(f'{name}={s:.1f}', (s, 0),
textcoords="offset points", xytext=(0, 15),
ha='center', fontsize=11)
ax.axhline(0, color='k', lw=0.5)
ax.axvline(0, color='k', lw=0.5)
ax.set_xlabel('$\\sigma$ [MPa]', fontsize=13)
ax.set_ylabel('$\\tau$ [MPa]', fontsize=13)
ax.set_title("Mohr's Circles (3D)", fontsize=14)
ax.set_aspect('equal')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
return ax
# --- 使用例: 圧力容器 + 曲げ ---
# 薄肉圧力容器の周方向応力と軸方向応力に曲げが加わるケース
stress = np.array([
[120, 45, 10], # σ_xx, τ_xy, τ_xz
[45, 60, -20], # τ_yx, σ_yy, τ_yz
[10, -20, -30] # τ_zx, τ_zy, σ_zz
], dtype=float)
analyzer = CombinedStressAnalyzer(stress, sigma_Y=350.0)
analyzer.summary()
fig, ax = plt.subplots(figsize=(8, 8))
analyzer.plot_mohr_3d(ax)
plt.tight_layout()
plt.savefig("combined_stress_analysis.png", dpi=150, bbox_inches='tight')
plt.show()
この総合解析コードの出力を見ると、入力した応力テンソルに対する主応力、フォンミーゼス応力、トレスカ応力、安全率がすべて一括で得られます。
特に注目すべき点は以下の通りです。
- フォンミーゼス安全率 $>$ トレスカ安全率が常に成り立ちます。これはトレスカ条件がフォンミーゼス条件よりも保守的（安全側）であることの直接的な表れです
- 応力不変量 $I_1$, $I_2$, $I_3$ の値が出力されています。これらは座標変換しても変わらない量であり、主応力が正しく計算されているかの検証に使えます。具体的には、$I_1 = \sigma_1 + \sigma_2 + \sigma_3$ を確認することで計算の整合性を確かめられます
- 3次元モールの応力円によって、最大せん断応力がどの主応力の組み合わせから生じるかが視覚的に明らかです。外側の大円の半径が $\tau_{\max}$ であり、これが降伏に最も関与するせん断応力です
補足: 応力集中と組合せ応力
実際の設計では、応力集中の影響を組合せ応力の計算に取り込む必要があります。キー溝、肩部、穴などの形状不連続部では、公称応力に応力集中係数 $K_t$ を乗じた値を使います。
曲げとねじりでそれぞれ異なる応力集中係数 $K_{t,b}$ と $K_{t,t}$ が適用される場合、フォンミーゼス応力は
$$ \begin{equation} \sigma_v = \sqrt{(K_{t,b}\,\sigma_b)^2 + 3(K_{t,t}\,\tau_t)^2} \end{equation} $$
となります。応力集中は局所的に応力を数倍に増大させるため、組合せ応力の評価においても無視できない要因です。
さらに、繰返し荷重が作用する場合は疲労破壊の考慮も必要です。疲労設計では、等価応力振幅をフォンミーゼスの式で計算し、S-N曲線と比較して寿命を推定します。
補足: 平面応力における第3主応力の扱い
平面応力と平面ひずみの条件下では $\sigma_z = 0$ ですが、主応力の計算では第3主応力 $\sigma_3 = 0$ を忘れてはいけません。
2次元の主応力公式で $\sigma_1$ と $\sigma_2$ を求めた場合、実際の3次元主応力は $\sigma_1$, $\sigma_2$, $0$ の3つです。特に $\sigma_1 > 0 > \sigma_2$ の場合、最大せん断応力は
$$ \tau_{\max} = \frac{\sigma_1 – \sigma_2}{2} $$
ですが、$\sigma_1 > \sigma_2 > 0$ の場合は
$$ \tau_{\max} = \frac{\sigma_1 – 0}{2} = \frac{\sigma_1}{2} $$
となり、$\sigma_3 = 0$ が最小主応力になります。この見落としは過大な安全率の見積もりにつながるため注意が必要です。
まとめ
本記事では、組合せ応力の理論を基礎から体系的に解説しました。
- 応力テンソルは、ある点の応力状態を $3 \times 3$ の対称行列として完全に記述します。6つの独立な成分（垂直応力3つ + せん断応力3つ）で応力状態が決まります
- 応力の座標変換は行列の相似変換 $\bm{\sigma}’ = \bm{R}\,\bm{\sigma}\,\bm{R}^T$ で表され、主応力は応力テンソルの固有値、主方向は固有ベクトルに対応します
- モールの応力円を使えば、応力変換の全体像を幾何学的に一目で把握できます。3次元では3つの円で構成されます
- トレスカの降伏条件は最大せん断応力に基づき、$\sigma_1$-$\sigma_2$ 平面上で六角形の降伏曲線を与えます。計算が単純で保守的です
- フォンミーゼスの降伏条件は偏差ひずみエネルギーに基づき、楕円形の降伏曲線を与えます。延性金属の実験結果とよく一致し、FEMでの数値解析に適しています
- 設計実務では、フォンミーゼス応力 $\sigma_v$ を降伏応力 $\sigma_Y$ で割った安全率 $n = \sigma_Y / \sigma_v$ が広く使われています
組合せ応力の理論は、材料力学の集大成ともいえる分野です。単軸応力から出発し、応力テンソル、座標変換、固有値問題、そして降伏条件へと論理的につながる一連の理論体系を理解することで、複雑な実構造物の強度評価に自信を持って取り組むことができるようになります。
次のステップとして、以下の記事も参考にしてください。
- 応力集中 — 形状不連続部での局所的な応力増大と応力集中係数
- 疲労破壊とS-N曲線 — 繰返し荷重による疲労寿命の評価手法
- 平面応力と平面ひずみ — 3次元問題の2次元近似とその適用条件
- せん断応力 — 梁のせん断応力分布の理論