建物の柱、橋梁のトラス部材、ロケットの胴体フレーム――これらの構造部材は圧縮荷重を受けたとき、材料が壊れるよりもずっと小さな荷重で横方向に急激にたわむことがあります。この座屈現象を理想的な条件で解析したのがオイラーの座屈荷重ですが、現実の柱は理想とはほど遠い環境に置かれています。
荷重が柱の中心軸からわずかにずれている（偏心荷重）、製造上の誤差で柱にわずかな曲がりがある（初期たわみ）、あるいは短い柱では材料が降伏してしまう（非弾性座屈）――こうした現実的な要因を考慮しなければ、安全な構造設計は不可能です。
実際、航空宇宙構造では軽量化のために安全率が小さく設定されるため、偏心や初期不整の影響を定量的に評価することが不可欠です。また、建築鉄骨の設計ではAISCなどの規格がジョンソン式やカラム曲線を採用しており、これらの理論的背景を理解することが設計実務の前提となります。
本記事の内容
- オイラー座屈の復習と、その適用限界
- 偏心荷重を受ける柱の解析とセカント公式の導出
- 初期たわみがある柱の挙動と荷重増幅係数
- 非弾性座屈（接線係数理論・二重係数理論）
- 柱の設計曲線（ジョンソン式・AISC規格）
- 有効座屈長さと境界条件の影響
- Pythonによる4種類の可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
オイラー座屈の復習と限界
オイラーの座屈荷重
まず、前回の記事で導出した結果を簡単に振り返りましょう。両端ピン支持の理想的な柱に軸圧縮力 $P$ が作用するとき、弾性曲線方程式から次の微分方程式が得られました。
$$ EI\frac{d^2v}{dx^2} + Pv = 0 $$
ここで $E$ はヤング率、$I$ は断面二次モーメント、$v(x)$ はたわみです。この方程式の非自明解が存在する条件として、オイラーの座屈荷重が求まります。
$$ P_{cr} = \frac{\pi^2 EI}{L^2} $$
臨界応力 $\sigma_{cr} = P_{cr}/A$ を断面二次半径 $r = \sqrt{I/A}$ を用いて細長比 $\lambda = L/r$ で表すと、次のようになります。
$$ \sigma_{cr} = \frac{\pi^2 E}{\lambda^2} $$
これは細長比 $\lambda$ が大きい（細長い）柱ほど座屈しやすいことを意味しています。
オイラー式の3つの限界
オイラーの座屈公式は非常にエレガントですが、以下の3つの重要な仮定に基づいています。
仮定1: 荷重が完全に中心軸上に作用する。現実には、荷重の作用点が柱断面の重心からずれている（偏心がある）ことがほとんどです。接合部の設計上、完全な中心荷重は実現困難です。
仮定2: 柱が完全にまっすぐである。製造・施工の段階で、柱にはわずかな曲がり（初期たわみ）が必ず存在します。圧延鋼材では一般に $L/1000$ 程度の初期たわみが許容されています。
仮定3: 材料が弾性範囲内にある。オイラー式の $\sigma_{cr}$ が降伏応力 $\sigma_y$ を超える領域（短い柱）では、この式はもはや適用できません。実際には材料の塑性化が先に起こるため、別の理論が必要です。
これらの限界を一つずつ克服していくのが、本記事の目的です。まず最初に、荷重が中心軸からずれている場合を考えましょう。
偏心荷重を受ける柱
問題設定
偏心荷重とは、圧縮力 $P$ が柱断面の重心から距離 $e$ だけずれて作用する状況です。身近な例で考えると、柱にブラケット（持ち出し梁）が取り付けられていて、その先端に荷重がかかる場合がこれにあたります。荷重がずれているということは、圧縮力に加えて曲げモーメント $M = Pe$ が最初から発生しているということです。
両端ピン支持の柱の中央に座標原点を取り、柱の長さを $L$、偏心距離を $e$ とします。端部での曲げモーメントは $M_0 = Pe$ です。
支配方程式の導出
柱の任意の位置 $x$ における曲げモーメントを考えます。ここでは座標原点を柱の左端に取ります。端部での偏心によるモーメント $Pe$ と、たわみ $v(x)$ による追加モーメント $Pv$ を合わせて、曲げ応力の理論から弾性曲線方程式は次のようになります。
$$ EI\frac{d^2v}{dx^2} = -P(v + e) $$
右辺の $-P(v + e)$ は、偏心 $e$ による初期モーメントとたわみ $v$ による追加モーメントの合計です。$k^2 = P/(EI)$ と置くと、次の微分方程式が得られます。
$$ \frac{d^2v}{dx^2} + k^2 v = -k^2 e $$
微分方程式の求解
この方程式は非斉次の2階線形常微分方程式です。一般解は斉次解と特殊解の和で表されます。
斉次方程式 $v” + k^2 v = 0$ の一般解は $A\sin kx + B\cos kx$ であり、特殊解は右辺が定数なので $v_p = -e$ です。したがって一般解は次のようになります。
$$ v(x) = A\sin kx + B\cos kx – e $$
境界条件 $v(0) = 0$ を適用すると、
$$ 0 = A \cdot 0 + B \cdot 1 – e \implies B = e $$
次に境界条件 $v(L) = 0$ を適用すると、
$$ 0 = A\sin kL + e\cos kL – e $$
$A$ について解くと、
$$ A = \frac{e(1 – \cos kL)}{\sin kL} $$
ここで三角関数の半角の公式 $1 – \cos\theta = 2\sin^2(\theta/2)$ と $\sin\theta = 2\sin(\theta/2)\cos(\theta/2)$ を使うと、
$$ A = \frac{2e\sin^2(kL/2)}{2\sin(kL/2)\cos(kL/2)} = e\tan\frac{kL}{2} $$
したがって、たわみの一般解は次の通りです。
$$ v(x) = e\left(\tan\frac{kL}{2}\sin kx + \cos kx – 1\right) $$
最大たわみとセカント公式
柱の中央 $x = L/2$ でたわみが最大となります。$x = L/2$ を代入すると、
$$ v_{max} = v\left(\frac{L}{2}\right) = e\left(\tan\frac{kL}{2}\sin\frac{kL}{2} + \cos\frac{kL}{2} – 1\right) $$
括弧内の $\tan(kL/2)\sin(kL/2) + \cos(kL/2)$ を計算します。$\tan\alpha = \sin\alpha / \cos\alpha$ を代入すると、
$$ \frac{\sin^2\alpha}{\cos\alpha} + \cos\alpha = \frac{\sin^2\alpha + \cos^2\alpha}{\cos\alpha} = \frac{1}{\cos\alpha} $$
ただし $\alpha = kL/2$ です。これを元の式に戻すと、
$$ v_{max} = e\left(\frac{1}{\cos(kL/2)} – 1\right) = e\left(\sec\frac{kL}{2} – 1\right) $$
$k = \sqrt{P/(EI)}$ を代入すると、
$$ \boxed{v_{max} = e\left[\sec\left(\frac{L}{2}\sqrt{\frac{P}{EI}}\right) – 1\right]} $$
この式から、$P \to P_{cr} = \pi^2 EI/L^2$ のとき $\sec$ の引数が $\pi/2$ に近づき、$v_{max} \to \infty$ となることがわかります。つまり、偏心荷重のある柱はオイラー荷重に達する前に大きなたわみが生じるのです。
セカント公式（最大応力の評価）
柱の中央断面における最大圧縮応力は、直接圧縮応力と曲げ応力の重ね合わせで求まります。
$$ \sigma_{max} = \frac{P}{A} + \frac{M_{max} \cdot c}{I} $$
ここで $c$ は断面の中立軸から最外縁までの距離、$M_{max} = P(e + v_{max}) = Pe \cdot \sec(kL/2)$ です。
$I = Ar^2$ の関係を使い、$P/A$ でくくると、
$$ \sigma_{max} = \frac{P}{A}\left(1 + \frac{ec}{r^2}\sec\frac{L}{2r}\sqrt{\frac{P}{EA}}\right) $$
$\sigma_{max} = \sigma_y$（降伏応力）と置くことで、降伏が始まる荷重を求める式が得られます。これがセカント公式（Secant Formula）です。
$$ \boxed{\sigma_y = \frac{P}{A}\left[1 + \frac{ec}{r^2}\sec\left(\frac{L}{2r}\sqrt{\frac{P}{EA}}\right)\right]} $$
この式は $P/A$ について陽に解くことができないため、数値的に解く必要があります。ここで $ec/r^2$ は偏心比（eccentricity ratio）と呼ばれ、荷重の偏心が応力に与える影響の度合いを表す無次元パラメータです。
セカント公式の重要な特徴は以下の通りです。
- 偏心 $e = 0$ のとき、$\sec$ の項が消えて $\sigma_{max} = P/A$ となり、オイラー座屈に帰着します
- $P \to P_{cr}$ のとき、$\sec$ が発散し $\sigma_{max} \to \infty$ となります
- 偏心が大きいほど、より小さい荷重で降伏に達します
偏心荷重の問題を理解したところで、次に「柱がそもそもまっすぐでない場合」を考えましょう。製造上の初期たわみは偏心荷重とは異なるメカニズムですが、結果として似た効果をもたらします。
初期たわみがある柱の挙動
初期たわみの定義
現実の柱は製造や施工の過程で完全にまっすぐにはなりません。この最初から存在する曲がりを初期たわみ（initial imperfection）と呼び、$v_0(x)$ で表します。もっとも単純なモデルとして、初期たわみが正弦波状であると仮定します。
$$ v_0(x) = a_0 \sin\frac{\pi x}{L} $$
ここで $a_0$ は初期たわみの最大振幅です。この仮定は、両端ピン支持の柱において座屈モードの第1次成分に一致するため、最も危険側（保守的）な見積もりを与えます。
支配方程式
初期たわみ $v_0(x)$ がある柱に軸圧縮力 $P$ が作用すると、追加たわみ $v(x)$ が生じます。全たわみは $v_0(x) + v(x)$ であり、曲げモーメントは全たわみに比例します。
$$ EI\frac{d^2v}{dx^2} = -P(v_0 + v) $$
重要なのは、左辺が追加たわみ $v$ のみの2階微分であることです（初期たわみは曲率を持っていますが、それに対応する内力は荷重を加える前からバランスしています）。$v_0 = a_0\sin(\pi x/L)$ と $k^2 = P/(EI)$ を代入すると、
$$ \frac{d^2v}{dx^2} + k^2 v = -k^2 a_0 \sin\frac{\pi x}{L} $$
微分方程式の求解
この方程式の特殊解を $v_p = C\sin(\pi x/L)$ の形で求めます。2回微分すると $v_p” = -C(\pi/L)^2\sin(\pi x/L)$ となるので、元の式に代入すると、
$$ -C\frac{\pi^2}{L^2}\sin\frac{\pi x}{L} + k^2 C\sin\frac{\pi x}{L} = -k^2 a_0 \sin\frac{\pi x}{L} $$
$\sin(\pi x/L)$ で割って $C$ について解きます。
$$ C\left(k^2 – \frac{\pi^2}{L^2}\right) = -k^2 a_0 $$
$$ C = \frac{-k^2 a_0}{k^2 – \pi^2/L^2} = \frac{a_0}{\pi^2/(k^2 L^2) – 1} $$
$k^2 = P/(EI)$ と $P_{cr} = \pi^2 EI/L^2$ の関係から $\pi^2/(k^2 L^2) = P_{cr}/P$ なので、
$$ C = \frac{a_0}{P_{cr}/P – 1} = \frac{a_0 P}{P_{cr} – P} $$
斉次解 $A\sin kx + B\cos kx$ に境界条件 $v(0) = v(L) = 0$ を適用すると $A = B = 0$（初期たわみが既に対称解を含んでいるため）。したがって追加たわみは次のようになります。
$$ v(x) = \frac{a_0 P}{P_{cr} – P}\sin\frac{\pi x}{L} $$
荷重増幅係数
全たわみ $v_{total} = v_0 + v$ を計算すると、
$$ v_{total}(x) = a_0\sin\frac{\pi x}{L} + \frac{a_0 P}{P_{cr} – P}\sin\frac{\pi x}{L} = a_0\sin\frac{\pi x}{L}\left(1 + \frac{P}{P_{cr} – P}\right) $$
括弧内を通分して整理すると、
$$ \boxed{v_{total}(x) = \frac{a_0}{1 – P/P_{cr}}\sin\frac{\pi x}{L}} $$
この結果は非常に重要です。係数 $1/(1 – P/P_{cr})$ は荷重増幅係数（amplification factor）と呼ばれ、初期たわみが圧縮荷重によってどれだけ増幅されるかを表しています。
荷重増幅係数の振る舞いを具体的に見てみましょう。
| $P/P_{cr}$ | 増幅係数 | 意味 |
|---|---|---|
| 0 | 1.0 | 荷重なし。初期たわみのまま |
| 0.5 | 2.0 | たわみが2倍に増幅 |
| 0.8 | 5.0 | たわみが5倍に増幅 |
| 0.9 | 10.0 | たわみが10倍に増幅 |
| 1.0 | $\infty$ | 発散（オイラー荷重に到達） |
$P$ が $P_{cr}$ に近づくと急激にたわみが増大する様子がわかります。この「緩やかに始まり、臨界荷重付近で急激に増大する」という特性は、偏心荷重の場合と定性的に同じです。
最大応力は中央断面で発生し、偏心荷重の場合と同様に直接圧縮と曲げの重ね合わせとなります。
$$ \sigma_{max} = \frac{P}{A} + \frac{P \cdot v_{max} \cdot c}{I} = \frac{P}{A}\left(1 + \frac{a_0 c}{r^2}\cdot\frac{1}{1 – P/P_{cr}}\right) $$
ここで $a_0 c / r^2$ が初期たわみに対する偏心比に相当するパラメータです。
ここまでで、偏心荷重と初期たわみという2つの現実的な不完全性の影響を見てきました。いずれも、オイラー荷重に達する前に柱が大きくたわみ、降伏に至ることを示しています。しかし、これまでの議論は全て材料が弾性範囲にあることを前提としていました。次に、材料の塑性を考慮した非弾性座屈の理論に進みましょう。
Pythonによる可視化 (1): セカント公式による荷重-変位曲線
セカント公式の挙動を視覚的に理解するため、異なる偏心比に対する荷重-たわみ曲線を描きます。
import numpy as np
import matplotlib.pyplot as plt
# --- 材料・断面パラメータ ---
E = 200e9 # ヤング率 [Pa]（鋼）
sigma_y = 250e6 # 降伏応力 [Pa]
L = 2.0 # 柱の長さ [m]
r = 0.05 # 断面二次半径 [m]
c = 0.075 # 中立軸から最外縁 [m]
A = 1e-3 # 断面積 [m^2]
I_val = A * r**2 # 断面二次モーメント [m^4]
P_cr = np.pi**2 * E * I_val / L**2 # オイラー座屈荷重
# --- セカント公式: 偏心比ごとの荷重 vs 最大たわみ ---
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左図: 荷重-たわみ曲線
ec_ratios = [0.0, 0.1, 0.3, 0.6, 1.0]
P_range = np.linspace(0.01 * P_cr, 0.98 * P_cr, 500)
for ec_r2 in ec_ratios:
e_val = ec_r2 * r**2 / c # 偏心距離
v_max = e_val * (1.0 / np.cos(L / (2 * r) * np.sqrt(P_range / (E * A))) - 1)
label = f'$ec/r^2$ = {ec_r2}'
axes[0].plot(v_max * 1000, P_range / P_cr, lw=2, label=label)
axes[0].axhline(1.0, color='k', ls='--', lw=1.5, label='Euler $P_{cr}$')
axes[0].set_xlabel('Max deflection $v_{max}$ [mm]')
axes[0].set_ylabel('$P / P_{cr}$')
axes[0].set_title('Secant formula: Load vs Deflection')
axes[0].legend(fontsize=9)
axes[0].set_xlim(0, 30)
axes[0].set_ylim(0, 1.2)
axes[0].grid(True, alpha=0.3)
# 右図: セカント公式による臨界応力 vs 細長比
lam = np.linspace(20, 200, 300)
axes[1].plot(lam, np.pi**2 * E / lam**2 / 1e6, 'k--', lw=2, label='Euler ($e=0$)')
for ec_r2 in [0.1, 0.3, 0.6, 1.0]:
# 各 lambda に対して sigma_y に達する P/A を数値的に求める
sigma_cr_arr = []
for lam_i in lam:
# セカント公式: sigma_y = (P/A)[1 + (ec/r^2)*sec(lam/(2)*sqrt(P/(EA)))]
# P/A = sigma とおくと: sigma_y = sigma[1 + ec_r2*sec(lam_i/2*sqrt(sigma/E))]
from scipy.optimize import brentq
def func(sigma):
arg = lam_i / 2 * np.sqrt(sigma / E)
if arg >= np.pi / 2 * 0.999:
return -sigma_y
return sigma * (1 + ec_r2 / np.cos(arg)) - sigma_y
try:
sol = brentq(func, 1e3, sigma_y - 1)
sigma_cr_arr.append(sol)
except:
sigma_cr_arr.append(np.nan)
axes[1].plot(lam, np.array(sigma_cr_arr) / 1e6, lw=2,
label=f'$ec/r^2$ = {ec_r2}')
axes[1].axhline(sigma_y / 1e6, color='r', ls=':', lw=1.5, label=f'$\\sigma_y$ = {sigma_y/1e6:.0f} MPa')
axes[1].set_xlabel('Slenderness ratio $\\lambda = L/r$')
axes[1].set_ylabel('Critical stress $\\sigma_{cr}$ [MPa]')
axes[1].set_title('Secant formula: Critical stress curves')
axes[1].legend(fontsize=9)
axes[1].set_ylim(0, 300)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左図の荷重-たわみ曲線から、2つの重要な特徴が読み取れます。
- 偏心比 $ec/r^2 = 0$（理想柱）のとき、荷重はオイラー座屈荷重 $P_{cr}$ まで垂直に上昇し、そこで急に横方向にたわみます。これが典型的な分岐座屈です。
- 偏心比が大きくなるにつれ、荷重-たわみ曲線は滑らかな曲線となり、明確な分岐点が消えます。荷重の増加とともに徐々にたわみが増え、$P_{cr}$ に漸近的に近づきます。これは偏心荷重下では「突然の座屈」ではなく「漸進的な曲げ破壊」が起こることを意味しています。
右図のセカント公式による臨界応力曲線は、オイラー曲線より常に下側に位置しています。偏心が大きいほど許容応力が低下し、特に中間的な細長比の領域でその影響が顕著です。設計上、偏心の影響を無視することの危険性が明確に示されています。
次に、初期たわみの影響をパラメトリックに解析し、さまざまな初期不整の大きさに対する柱の応答を調べましょう。
Pythonによる可視化 (2): 初期たわみの影響のパラメトリック解析
初期たわみの振幅 $a_0$ をパラメータとして、荷重比 $P/P_{cr}$ に対する最大たわみと最大応力の変化を可視化します。
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ ---
E = 200e9; sigma_y = 250e6; L = 2.0; r = 0.05; c = 0.075
A = 1e-3; I_val = A * r**2
P_cr = np.pi**2 * E * I_val / L**2
# --- 初期たわみのパラメトリック解析 ---
fig, axes = plt.subplots(1, 3, figsize=(18, 5.5))
P_ratio = np.linspace(0.01, 0.99, 500) # P/Pcr
a0_values = [L/10000, L/1000, L/500, L/200, L/100] # 初期たわみ振幅
a0_labels = ['L/10000', 'L/1000', 'L/500', 'L/200', 'L/100']
# (a) 荷重比 vs 最大たわみ
for a0, lab in zip(a0_values, a0_labels):
amp = 1.0 / (1.0 - P_ratio) # 増幅係数
v_max = a0 * amp
axes[0].plot(v_max * 1000, P_ratio, lw=2, label=f'$a_0$ = {lab}')
axes[0].axhline(1.0, color='k', ls='--', lw=1.5, alpha=0.7)
axes[0].set_xlabel('Max deflection [mm]')
axes[0].set_ylabel('$P / P_{cr}$')
axes[0].set_title('(a) Load ratio vs Max deflection')
axes[0].legend(fontsize=9); axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 50); axes[0].set_ylim(0, 1.1)
# (b) 増幅係数 vs 荷重比
axes[1].plot(P_ratio, 1.0 / (1.0 - P_ratio), 'b-', lw=2.5)
axes[1].axhline(1.0, color='gray', ls=':', lw=1)
axes[1].fill_between(P_ratio, 1, 1.0 / (1.0 - P_ratio), alpha=0.15, color='blue')
axes[1].set_xlabel('$P / P_{cr}$')
axes[1].set_ylabel('Amplification factor $1/(1-P/P_{cr})$')
axes[1].set_title('(b) Amplification factor')
axes[1].set_ylim(0, 20); axes[1].grid(True, alpha=0.3)
# 注目点を注記
for pr, af_label in [(0.5, '2x'), (0.8, '5x'), (0.9, '10x')]:
af = 1.0 / (1.0 - pr)
axes[1].plot(pr, af, 'ro', ms=6)
axes[1].annotate(af_label, (pr, af), textcoords="offset points",
xytext=(8, 5), fontsize=10, color='red')
# (c) 荷重比 vs 最大応力
for a0, lab in zip(a0_values, a0_labels):
P_vals = P_ratio * P_cr
sigma_max = P_vals / A * (1 + a0 * c / r**2 / (1.0 - P_ratio))
axes[2].plot(P_ratio, sigma_max / 1e6, lw=2, label=f'$a_0$ = {lab}')
axes[2].axhline(sigma_y / 1e6, color='r', ls='--', lw=1.5, label=f'$\\sigma_y$ = {sigma_y/1e6:.0f} MPa')
axes[2].set_xlabel('$P / P_{cr}$')
axes[2].set_ylabel('$\\sigma_{max}$ [MPa]')
axes[2].set_title('(c) Max stress vs Load ratio')
axes[2].legend(fontsize=8); axes[2].grid(True, alpha=0.3)
axes[2].set_ylim(0, 500)
plt.tight_layout()
plt.show()
3つのグラフから、初期たわみの影響について重要な知見が得られます。
- 図(a): 荷重-たわみ曲線 — 初期たわみが大きいほど、低い荷重で大きなたわみが発生します。$a_0 = L/100$ の場合、$P/P_{cr} = 0.5$ の段階で既に数十mmのたわみが生じており、実用上の限界に達する可能性があります。一方、$a_0 = L/10000$ のような微小な初期たわみでは曲線がほぼ垂直に立ち上がり、理想柱に近い振る舞いを示します。
- 図(b): 増幅係数 — $P/P_{cr} < 0.5$ の領域では増幅係数は2以下であり、初期たわみの影響は限定的です。しかし $P/P_{cr} > 0.8$ を超えると急激に増大し、$P/P_{cr} = 0.9$ で10倍、$P/P_{cr} = 0.95$ で20倍に達します。この非線形性が座屈問題の本質的な難しさです。
- 図(c): 最大応力 — 初期たわみが大きいほど、低い荷重比で降伏応力に到達します。これは実際の柱の許容荷重がオイラー荷重よりも大幅に小さくなることを意味しており、設計における安全率の根拠となっています。
初期たわみと偏心荷重の両方が、柱の耐荷力をオイラー理論から低下させることがわかりました。しかし、まだもう一つの重要な問題が残っています。細長比が小さい（短い太い）柱では、オイラー応力が降伏応力を超えてしまいます。この領域での座屈を扱うには、材料の非線形性を考慮する理論が必要です。
非弾性座屈
なぜ非弾性座屈が重要か
オイラーの座屈公式が予測する臨界応力 $\sigma_{cr} = \pi^2 E / \lambda^2$ は、細長比 $\lambda$ が小さくなると降伏応力 $\sigma_y$ を超えます。この境界となる細長比を遷移細長比 $\lambda_c$ と呼びます。
$$ \sigma_y = \frac{\pi^2 E}{\lambda_c^2} \implies \lambda_c = \pi\sqrt{\frac{E}{\sigma_y}} $$
鋼材（$E = 200$ GPa, $\sigma_y = 250$ MPa）の場合、
$$ \lambda_c = \pi\sqrt{\frac{200 \times 10^9}{250 \times 10^6}} \approx 88.9 $$
つまり細長比が約89以下の柱では、座屈前に材料が降伏してしまうため、オイラー式は適用できません。実際の構造部材の多くはこの「中間柱」の領域に属しており、非弾性座屈の理論が設計実務上不可欠なのです。
接線係数理論（Engesser-Shanley理論）
材料が塑性変形を開始した後は、応力-ひずみ関係が非線形になります。フックの法則が成り立たなくなった後の座屈をどう扱えばよいでしょうか？
Engesserは1889年に画期的なアイデアを提唱しました。応力-ひずみ曲線の現在の点における接線の傾き（接線係数 $E_t$）をヤング率 $E$ の代わりに用いるというものです。
$$ E_t = \frac{d\sigma}{d\epsilon} $$
この $E_t$ はひずみレベルに依存し、塑性域では $E_t < E$ です。接線係数理論による臨界応力は、オイラー式の $E$ を $E_t$ で置き換えたものです。
$$ \boxed{\sigma_{cr}^{(T)} = \frac{\pi^2 E_t}{\lambda^2}} $$
しかし、ここで一つ問題が生じます。$E_t$ 自体が応力（したがって $\sigma_{cr}$）に依存するため、この式は陰的な方程式です。具体的な応力-ひずみ関係のモデルが必要になります。
Ramberg-Osgoodモデル
非弾性挙動を表す代表的なモデルとして、Ramberg-Osgood式があります。
$$ \epsilon = \frac{\sigma}{E} + K\left(\frac{\sigma}{\sigma_y}\right)^n $$
ここで $K$ と $n$ は材料定数です。接線係数は逆数の微分から得られます。
$$ \frac{d\epsilon}{d\sigma} = \frac{1}{E} + \frac{Kn}{\sigma_y}\left(\frac{\sigma}{\sigma_y}\right)^{n-1} $$
したがって、
$$ E_t = \frac{1}{\frac{1}{E} + \frac{Kn}{\sigma_y}\left(\frac{\sigma}{\sigma_y}\right)^{n-1}} $$
応力が増加するにつれて $E_t$ は単調に減少し、$\sigma \to 0$ では $E_t \to E$ に復帰します。
二重係数理論（Reduced Modulus Theory）
接線係数理論に対して、Considere と Engesser はもう一つの理論を提唱しました。座屈が始まるとき、柱の断面内で凹側は追加圧縮（接線係数 $E_t$）、凸側は除荷（ヤング率 $E$）が作用するという考え方です。
この場合、断面の曲げ剛性は $E$ でも $E_t$ でもなく、両者を断面形状に応じて加重平均した二重係数（reduced modulus）$E_r$ で表されます。
$$ \sigma_{cr}^{(R)} = \frac{\pi^2 E_r}{\lambda^2} $$
矩形断面の場合、二重係数は次のように表されます。
$$ E_r = \frac{4EE_t}{(\sqrt{E} + \sqrt{E_t})^2} $$
$E_t < E$ のとき、常に $E_t < E_r < E$ が成り立ちます。したがって、
$$ \sigma_{cr}^{(T)} < \sigma_{cr}^{(R)} < \sigma_{cr}^{(Euler)} $$
という関係があります。
Shanleyの解決
1947年、Shanleyは接線係数理論と二重係数理論の論争に決着をつけました。彼の結論は以下の通りです。
- 座屈は接線係数荷重 $P_T$ で始まる
- 座屈後、荷重は増加しながら二重係数荷重 $P_R$ に漸近する
- したがって、設計上は接線係数理論による荷重を用いるのが安全側
この結論は、座屈が完全な分岐現象ではなく、微小な横変位と軸荷重の増加が同時に進む過程であることを意味しています。Shanleyの理論は実験結果ともよく一致し、現代の設計規格の基礎となっています。
非弾性座屈の理論は材料の非線形性を正確にモデル化する必要があり、設計実務では簡便な近似式が用いられます。次のセクションでは、実際の設計で使われるジョンソン式とAISC規格の柱曲線を見ていきましょう。
柱の設計曲線
ジョンソンの放物線公式
エンジニアが日常的に使う設計式として、ジョンソンの放物線公式（Johnson’s parabolic formula）があります。これは非弾性座屈領域をオイラー曲線に滑らかに接続する放物線で近似するものです。
基本的な考え方はシンプルです。$\lambda = 0$（極端に短い柱）で $\sigma_{cr} = \sigma_y$ となり、$\lambda = \lambda_c$（遷移細長比）でオイラー曲線に接する放物線を見つけます。この条件を満たす式は次の通りです。
$$ \boxed{\sigma_{cr}^{(J)} = \sigma_y\left(1 – \frac{\sigma_y}{4\pi^2 E}\lambda^2\right) \quad (\lambda \leq \lambda_c)} $$
この式の導出を確認しましょう。$\sigma_{cr} = a – b\lambda^2$ の形で、$\lambda = 0$ で $\sigma_{cr} = \sigma_y$ なので $a = \sigma_y$。オイラー曲線 $\sigma_{cr} = \pi^2 E / \lambda^2$ との接点条件（値と傾きが一致）から、
$$ \sigma_y – b\lambda_c^2 = \frac{\pi^2 E}{\lambda_c^2} $$
$\lambda_c^2 = \pi^2 E / \sigma_y$ を代入すると、
$$ \sigma_y – b \cdot \frac{\pi^2 E}{\sigma_y} = \sigma_y \cdot \frac{\sigma_y}{\sigma_y} = \sigma_y \cdot \frac{\pi^2 E / \lambda_c^2}{\sigma_y} $$
整理すると $b = \sigma_y^2 / (4\pi^2 E)$ が得られ、ジョンソン式が導かれます。
AISC規格の柱曲線
アメリカ鉄鋼建設協会（AISC: American Institute of Steel Construction）の規格では、安全率を含んだ許容応力の形で柱の設計式が規定されています。AISC-ASD（Allowable Stress Design）の柱曲線は次の通りです。
遷移細長比 $C_c = \sqrt{2\pi^2 E / \sigma_y}$ を境界として、
弾性領域（$\lambda > C_c$）：
$$ \sigma_{allow} = \frac{12\pi^2 E}{23\lambda^2} $$
これはオイラー応力に安全率 $23/12 \approx 1.92$ を適用したものです。
非弾性領域（$\lambda \leq C_c$）：
$$ \sigma_{allow} = \frac{\sigma_y}{FS}\left[1 – \frac{\lambda^2}{2C_c^2}\right] $$
ここで $FS$ は細長比に依存する安全率です。
$$ FS = \frac{5}{3} + \frac{3\lambda}{8C_c} – \frac{\lambda^3}{8C_c^3} $$
AISC規格は、偏心荷重、初期たわみ、残留応力などの現実的な不完全性を経験的に取り込んだ設計曲線であり、ジョンソン式よりも保守的な値を与えます。
設計曲線の理論的背景を理解した上で、Pythonで各曲線を比較しましょう。
Pythonによる可視化 (3): オイラー曲線とジョンソン曲線の比較
import numpy as np
import matplotlib.pyplot as plt
# --- 材料パラメータ ---
E = 200e9 # Pa
sigma_y = 250e6 # Pa
# --- 細長比の範囲 ---
lam = np.linspace(1, 250, 1000)
# --- オイラー曲線 ---
sigma_euler = np.pi**2 * E / lam**2
# --- 遷移細長比 ---
lam_c = np.pi * np.sqrt(E / sigma_y) # ≈ 88.9
C_c = np.sqrt(2 * np.pi**2 * E / sigma_y) # AISC 遷移細長比 ≈ 125.7
# --- ジョンソンの放物線公式 ---
sigma_johnson = np.where(
lam <= lam_c,
sigma_y * (1 - sigma_y / (4 * np.pi**2 * E) * lam**2),
sigma_euler
)
# --- AISC 許容応力 ---
FS = 5/3 + 3*lam/(8*C_c) - lam**3/(8*C_c**3)
sigma_aisc_inelastic = sigma_y / FS * (1 - lam**2 / (2 * C_c**2))
sigma_aisc_elastic = 12 * np.pi**2 * E / (23 * lam**2)
sigma_aisc = np.where(lam <= C_c, sigma_aisc_inelastic, sigma_aisc_elastic)
# --- 接線係数理論（Ramberg-Osgoodモデル）---
n_ro = 20 # 冪指数
K_ro = 0.002 # 0.2%耐力に対応
sigma_range = np.linspace(1e6, sigma_y * 0.999, 500)
Et = 1.0 / (1.0/E + K_ro * n_ro / sigma_y * (sigma_range/sigma_y)**(n_ro-1))
lam_tangent = np.pi * np.sqrt(Et / sigma_range)
# --- プロット ---
fig, ax = plt.subplots(figsize=(10, 7))
ax.plot(lam, sigma_euler / 1e6, 'b--', lw=2, label='Euler curve')
ax.plot(lam, sigma_johnson / 1e6, 'r-', lw=2.5, label="Johnson's parabola")
ax.plot(lam, sigma_aisc / 1e6, 'g-', lw=2.5, label='AISC-ASD allowable')
ax.plot(lam_tangent, sigma_range / 1e6, 'm-', lw=2, label='Tangent modulus (R-O)')
ax.axhline(sigma_y / 1e6, color='k', ls=':', lw=1.5, alpha=0.7, label=f'$\\sigma_y$ = {sigma_y/1e6:.0f} MPa')
ax.axvline(lam_c, color='gray', ls=':', lw=1, alpha=0.5)
ax.axvline(C_c, color='gray', ls='--', lw=1, alpha=0.5)
ax.annotate(f'$\\lambda_c$ = {lam_c:.0f}', xy=(lam_c, 10), fontsize=10, ha='center',
color='gray')
ax.annotate(f'$C_c$ = {C_c:.0f}', xy=(C_c, 10), fontsize=10, ha='center',
color='gray')
# 領域の塗り分け
ax.fill_betweenx([0, sigma_y/1e6], 0, lam_c, alpha=0.05, color='red',
label='Inelastic region')
ax.fill_betweenx([0, sigma_y/1e6], lam_c, 250, alpha=0.05, color='blue',
label='Elastic region')
ax.set_xlabel('Slenderness ratio $\\lambda = L_e/r$', fontsize=12)
ax.set_ylabel('Critical stress $\\sigma_{cr}$ [MPa]', fontsize=12)
ax.set_title('Column design curves comparison', fontsize=14)
ax.legend(fontsize=10, loc='upper right')
ax.set_xlim(0, 250)
ax.set_ylim(0, 400)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このグラフから、柱の設計に関する重要な知見が得られます。
- 弾性領域（$\lambda > \lambda_c \approx 89$）では、オイラー曲線（青破線）が座屈応力の上界を与えます。ジョンソン曲線はこの領域でオイラー曲線と一致します。AISC許容応力（緑線）はオイラー応力に約1.92の安全率を適用しているため、大幅に低い値を示します。
- 非弾性領域（$\lambda < \lambda_c$）では、オイラー曲線が降伏応力を超えて非物理的な値を取るのに対し、ジョンソン曲線は放物線状に $\sigma_y$ へ滑らかに接続します。接線係数理論（マゼンタ線）はRamberg-Osgoodモデルに基づく理論曲線で、ジョンソン曲線と定性的に一致していることが確認できます。
- AISC許容応力は全領域で最も低い値を与えます。これは安全率に加え、残留応力や初期不整などの影響を経験的に含んでいるためです。実際の設計では、この曲線以下で部材を使用することが求められます。
ここまでで、さまざまな理論に基づく柱の設計曲線を比較しました。次に、座屈荷重に大きな影響を与える境界条件（端末条件）について、より詳細に見ていきましょう。
有効座屈長さと境界条件の影響
有効座屈長さの概念
オイラーの座屈の記事でも触れましたが、異なる端末条件での座屈荷重は有効座屈長さ $L_e = KL$ を導入することで統一的に表現できます。
$$ P_{cr} = \frac{\pi^2 EI}{L_e^2} = \frac{\pi^2 EI}{(KL)^2} $$
ここで $K$ は座屈長さ係数（effective length factor）で、境界条件によって異なる値を取ります。直感的には、$K$ は「柱が両端ピン支持の場合の何倍の長さと等価か」を表しています。
各境界条件の座屈モード
4つの代表的な境界条件について、座屈モードの導出を概観します。
Case 1: 両端ピン（$K = 1.0$）
これは基本ケースです。境界条件 $v(0) = v(L) = 0$ から座屈モードは $v(x) = A\sin(\pi x / L)$ となります。
Case 2: 一端固定・一端自由（片持ち柱、$K = 2.0$）
固定端 $x = 0$ で $v(0) = 0$, $v'(0) = 0$。自由端 $x = L$ で $EIv”(L) = 0$（モーメントなし）。自由端では軸力のモーメントと曲げモーメントのバランスから、
$$ EI\frac{d^2v}{dx^2} = -P(v – \delta) $$
ここで $\delta = v(L)$ は自由端のたわみです。$k^2 = P/(EI)$ として解くと、$v(x) = \delta(1 – \cos kx)$ が得られ、$v(L) = \delta$ の条件から $\cos kL = 0$、すなわち $kL = \pi/2$ が最小の座屈条件です。したがって、
$$ P_{cr} = \frac{\pi^2 EI}{4L^2} = \frac{\pi^2 EI}{(2L)^2} $$
有効長さは $L_e = 2L$ で、柱が「2倍の長さの両端ピン柱」と等価であることを意味しています。
Case 3: 両端固定（$K = 0.5$）
両端で $v(0) = v(L) = 0$ かつ $v'(0) = v'(L) = 0$ です。対称性から座屈モードは $v(x) = A(1 – \cos(2\pi x / L))$ の形となり、座屈条件は $kL = 2\pi$ です。
$$ P_{cr} = \frac{4\pi^2 EI}{L^2} = \frac{\pi^2 EI}{(L/2)^2} $$
有効長さは $L_e = L/2$ で、両端を固定することで座屈荷重が4倍に増加します。
Case 4: 一端固定・一端ピン（$K \approx 0.7$）
固定端で $v(0) = 0$, $v'(0) = 0$、ピン端で $v(L) = 0$, $v”(L) = 0$ の条件を適用します。この場合、ピン端に反力モーメントが作用しないため、固定端に水平反力が発生する不静定問題となります。超越方程式 $\tan kL = kL$ を数値的に解くと、最小根は $kL \approx 4.493$ であり、
$$ P_{cr} = \frac{(4.493)^2 EI}{L^2} \approx \frac{\pi^2 EI}{(0.6992L)^2} $$
有効長さは $L_e \approx 0.7L$ です。
有効座屈長さの一覧
| 端末条件 | $K = L_e/L$ | $P_{cr}/({\pi^2EI/L^2})$ | 拘束の強さ |
|---|---|---|---|
| 両端ピン | 1.0 | 1.0 | 基準 |
| 一端固定・一端自由 | 2.0 | 0.25 | 最も弱い |
| 両端固定 | 0.5 | 4.0 | 最も強い |
| 一端固定・一端ピン | 0.7 | 2.04 | やや強い |
端末条件の拘束が強いほど（変位・回転の自由度が少ないほど）、有効座屈長さが短くなり座屈荷重が大きくなります。この直感的な理解は重要です。構造設計では、接合部の剛性を高めて拘束条件を改善することが座屈対策の有効な手段となります。
実際の構造における注意点
理論上の $K$ 値は理想的な境界条件を仮定していますが、実際の構造では以下の点に注意が必要です。
- 完全な固定端は存在しない — 現実のボルト接合や溶接接合にはある程度の回転剛性があり、$K$ は理論値と完全ピンの間の値を取ります
- ブレースの効果 — 中間点にブレース（振れ止め）を設置すると有効座屈長さが短くなり、座屈荷重を大幅に向上できます
- フレーム内の柱 — 建築フレーム内の柱は、梁との接合による拘束と隣接柱の挙動の影響を受けるため、$K$ の評価が複雑になります。この場合、FEMによる固有値解析が有効です
それでは、各境界条件での座屈モードをPythonで可視化して、直感的な理解を深めましょう。
Pythonによる可視化 (4): 各境界条件での座屈モード可視化
4つの代表的な境界条件に対して、第1〜第3座屈モードを可視化します。
import numpy as np
import matplotlib.pyplot as plt
L = 1.0
x = np.linspace(0, L, 500)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
titles = ['(a) Pin-Pin (K=1.0)', '(b) Fixed-Free (K=2.0)',
'(c) Fixed-Fixed (K=0.5)', '(d) Fixed-Pin (K≈0.7)']
colors = ['#1f77b4', '#ff7f0e', '#2ca02c']
# --- (a) 両端ピン ---
ax = axes[0, 0]
for n in [1, 2, 3]:
v = np.sin(n * np.pi * x / L)
ax.plot(v, x / L, lw=2.5, color=colors[n-1], label=f'n={n}')
ax.axvline(0, color='k', lw=0.8)
ax.plot([0, 0], [0, 1], 'k|', ms=15, mew=3) # ピン記号
ax.set_title(titles[0], fontsize=12)
# --- (b) 一端固定・一端自由 ---
ax = axes[0, 1]
for n in [1, 2, 3]:
kn_L = (2*n - 1) * np.pi / 2 # kL = pi/2, 3pi/2, 5pi/2
v = 1 - np.cos(kn_L * x / L)
v = v / np.max(np.abs(v)) # 正規化
ax.plot(v, x / L, lw=2.5, color=colors[n-1], label=f'n={n}')
ax.axvline(0, color='k', lw=0.8)
ax.set_title(titles[1], fontsize=12)
# --- (c) 両端固定 ---
ax = axes[1, 0]
for n in [1, 2, 3]:
kn_L = 2 * n * np.pi # kL = 2pi, 4pi, 6pi
# 座屈モード: v = A(1 - cos(kn*x)) - Bx は対称/反対称
if n == 1:
v = 1 - np.cos(2 * np.pi * x / L)
elif n == 2:
v = np.sin(2 * np.pi * x / L) - 2 * np.pi * x / L + np.pi
v = np.sin(2 * np.pi * x / L) # 反対称モード
else:
v = 1 - np.cos(4 * np.pi * x / L)
v = v / np.max(np.abs(v))
ax.plot(v, x / L, lw=2.5, color=colors[n-1], label=f'n={n}')
ax.axvline(0, color='k', lw=0.8)
ax.set_title(titles[2], fontsize=12)
# --- (d) 一端固定・一端ピン ---
ax = axes[1, 1]
# 超越方程式 tan(kL) = kL の解
kL_roots = [4.4934, 7.7253, 10.9041]
for i, kL in enumerate(kL_roots):
# v = A[sin(kx) - (kx)*sin(kL)/(kL)] ... 簡略化モード形状
v = np.sin(kL * x / L) - (kL * x / L) * np.sin(kL) / kL
v = v / np.max(np.abs(v)) if np.max(np.abs(v)) > 0 else v
ax.plot(v, x / L, lw=2.5, color=colors[i], label=f'n={i+1}')
ax.axvline(0, color='k', lw=0.8)
ax.set_title(titles[3], fontsize=12)
# 共通設定
for ax in axes.flat:
ax.set_xlabel('Normalized deflection $v/v_{max}$')
ax.set_ylabel('Position $x/L$')
ax.legend(fontsize=10, loc='best')
ax.grid(True, alpha=0.3)
ax.set_xlim(-1.3, 1.3)
ax.set_ylim(0, 1)
# 境界マーカー
ax.invert_yaxis()
plt.suptitle('Buckling mode shapes for different boundary conditions',
fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
4つの境界条件における座屈モードから、重要な物理的洞察が得られます。
- 両端ピン（a） — 第1モードは半波の正弦波で、柱全体が一方向に膨らみます。高次モードは $n$ 個の半波を持ち、節点（たわみゼロ）が柱内部に現れます。高次モードの座屈荷重は $n^2$ に比例して増大するため、通常は第1モードのみが実用上重要です。
- 一端固定・一端自由（b） — 固定端で変位・傾きがゼロ、自由端で最大変位となります。第1モードは「四分の一波長」であり、有効座屈長さが物理的な柱長の2倍になることが視覚的に確認できます。実際、この変形形状を対称に延長すると、$2L$ の両端ピン柱の座屈モードと一致します。
- 両端固定（c） — 両端で変位・傾きがゼロのため、第1モードは全波の余弦波（$1 – \cos$型）です。変曲点が $L/4$ と $3L/4$ に現れ、この間の区間が有効座屈長さ $L/2$ に対応しています。拘束が最も強く、座屈荷重は両端ピンの4倍です。
- 一端固定・一端ピン（d） — 非対称な境界条件のため、座屈モードの節点位置が非対称になります。固定端側でより急な変形が生じる一方、ピン端側では緩やかに変形します。
偏心荷重・初期たわみ・非弾性座屈の相互関係
ここまで、偏心荷重、初期たわみ、非弾性座屈を個別に議論してきました。実際の柱ではこれらが同時に作用することが一般的です。ここでは、それぞれの影響がどのように組み合わさるかを整理します。
等価偏心の概念
初期たわみ $a_0$ のある柱と偏心距離 $e$ の柱は、荷重-変位の応答が定性的に似ています。実際、以下の「等価偏心」の概念でこの2つを統一的に扱えます。
初期たわみによる中央断面の最大応力は次の通りでした。
$$ \sigma_{max} = \frac{P}{A}\left(1 + \frac{a_0 c}{r^2}\cdot\frac{1}{1 – P/P_{cr}}\right) $$
一方、セカント公式は次の通りです。
$$ \sigma_{max} = \frac{P}{A}\left(1 + \frac{ec}{r^2}\sec\frac{\pi}{2}\sqrt{\frac{P}{P_{cr}}}\right) $$
$P/P_{cr}$ が小さい領域では $\sec(\pi\sqrt{P/P_{cr}}/2) \approx 1/(1 – P/P_{cr})$ と近似できるため、初期たわみ振幅 $a_0$ は偏心距離 $e$ と等価に扱えます。この近似は $P/P_{cr} < 0.6$ 程度まで良好です。
不完全性に対する感度
応力集中の問題と同様に、座屈問題でも不完全性に対する構造の感度が重要な設計指標です。
偏心比 $ec/r^2$ や初期たわみ比 $a_0/L$ がわずかに変化したとき、臨界荷重がどれだけ変化するかを不完全性感度（imperfection sensitivity）と呼びます。前述の荷重-たわみ曲線から明らかなように、座屈問題は本質的に不完全性に敏感です。特に $P/P_{cr}$ が1に近い領域での急激な応答変化は、安全率設定の根拠となっています。
Perry-Robertsonの公式
イギリスの設計規格で広く使われているPerry-Robertsonの公式は、初期たわみの影響を閉じた形で取り込んだ設計式です。
$$ \sigma_{cr} = \frac{1}{2}\left[(\sigma_y + (1+\eta)\sigma_E) – \sqrt{(\sigma_y + (1+\eta)\sigma_E)^2 – 4\sigma_y\sigma_E}\right] $$
ここで $\sigma_E = \pi^2 E / \lambda^2$（オイラー応力）、$\eta = a_0 c / r^2$（初期たわみに基づくPerry因子）です。
$\eta = 0$（完全柱）のとき $\sigma_{cr} = \min(\sigma_y, \sigma_E)$ に帰着し、$\eta$ が増加するほど臨界応力は低下します。この式はジョンソン式と異なり、初期たわみの影響を物理的に取り込んでいるため、より合理的な設計基準を与えます。
実際の構造設計では、これらの不完全性の影響を定量的に評価し、適切な安全率を設定することが極めて重要です。特に航空宇宙分野では軽量化要求が厳しいため、不完全性の大きさを製造品質管理によって制御し、可能な限り安全率を小さく抑えることが求められます。
次のセクションでは、Perry-Robertson公式とジョンソン式の統合的な比較をPythonで行い、不完全性パラメータの影響を定量的に評価します。
Pythonによる可視化 (5): Perry-Robertson公式とジョンソン式の比較
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ ---
E = 200e9; sigma_y = 250e6
lam = np.linspace(5, 250, 1000)
sigma_E = np.pi**2 * E / lam**2 # オイラー応力
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# --- 左図: Perry-Robertson公式（eta パラメトリック）---
ax = axes[0]
eta_values = [0.0, 0.1, 0.3, 0.5, 1.0]
for eta in eta_values:
a = sigma_y + (1 + eta) * sigma_E
b = a**2 - 4 * sigma_y * sigma_E
b = np.maximum(b, 0) # 数値的安全策
sigma_pr = 0.5 * (a - np.sqrt(b))
sigma_pr = np.minimum(sigma_pr, sigma_y) # 降伏応力で上限
ax.plot(lam, sigma_pr / 1e6, lw=2, label=f'$\\eta$ = {eta}')
ax.plot(lam, sigma_E / 1e6, 'k--', lw=1.5, alpha=0.5, label='Euler')
ax.axhline(sigma_y / 1e6, color='r', ls=':', lw=1.5, alpha=0.5)
ax.set_xlabel('Slenderness ratio $\\lambda$', fontsize=11)
ax.set_ylabel('Critical stress [MPa]', fontsize=11)
ax.set_title('Perry-Robertson formula', fontsize=12)
ax.legend(fontsize=9); ax.grid(True, alpha=0.3)
ax.set_xlim(0, 250); ax.set_ylim(0, 300)
# --- 右図: 各設計曲線の統合比較 ---
ax = axes[1]
# Euler
ax.plot(lam, sigma_E / 1e6, 'b--', lw=2, label='Euler')
# Johnson
lam_c = np.pi * np.sqrt(E / sigma_y)
sigma_j = np.where(lam <= lam_c,
sigma_y * (1 - sigma_y / (4*np.pi**2*E) * lam**2),
sigma_E)
ax.plot(lam, sigma_j / 1e6, 'r-', lw=2.5, label='Johnson')
# Perry-Robertson (eta=0.3 : 典型的)
eta_typ = 0.3
a_pr = sigma_y + (1 + eta_typ) * sigma_E
b_pr = np.maximum(a_pr**2 - 4*sigma_y*sigma_E, 0)
sigma_pr = np.minimum(0.5*(a_pr - np.sqrt(b_pr)), sigma_y)
ax.plot(lam, sigma_pr / 1e6, 'g-', lw=2.5, label=f'Perry-Robertson ($\\eta$={eta_typ})')
# Rankine-Gordon
sigma_rg = sigma_y / (1 + sigma_y * lam**2 / (np.pi**2 * E))
ax.plot(lam, sigma_rg / 1e6, 'm--', lw=2, label='Rankine-Gordon')
ax.axhline(sigma_y / 1e6, color='k', ls=':', lw=1.5, alpha=0.5)
ax.set_xlabel('Slenderness ratio $\\lambda$', fontsize=11)
ax.set_ylabel('Critical stress [MPa]', fontsize=11)
ax.set_title('Design curves comparison', fontsize=12)
ax.legend(fontsize=9); ax.grid(True, alpha=0.3)
ax.set_xlim(0, 250); ax.set_ylim(0, 300)
plt.tight_layout()
plt.show()
この2つのグラフから、設計曲線の特性と使い分けについて重要な知見が得られます。
- 左図（Perry-Robertson公式のパラメトリック解析） — Perry因子 $\eta$ が増加するにつれて、臨界応力曲線が全体的に低下します。$\eta = 0$（完全柱）ではオイラー曲線と降伏応力の小さい方に一致しますが、$\eta = 1.0$ では特に中間細長比（$\lambda \approx 50 \sim 120$）の領域で大幅な低下が見られます。この領域が設計上最も注意を要する部分です。
- 右図（各設計曲線の統合比較） — 4つの代表的な設計曲線を重ねてプロットしています。ジョンソン式とPerry-Robertson式は中間細長比域で異なる値を取りますが、高細長比域ではともにオイラー曲線に収束します。Rankine-Gordonの式（$\sigma_{cr} = \sigma_y / (1 + \sigma_y\lambda^2/(\pi^2 E))$）は最も古典的な経験式で、全領域で滑らかな曲線を与えます。設計者はこれらの曲線の違いと、各規格が想定している不完全性のレベルを理解した上で適切な式を選択する必要があります。
座屈設計のまとめと実務上の指針
本記事で扱った各理論の適用範囲を整理します。
適用範囲の対応表
| 柱の分類 | 細長比の範囲 | 適用理論 | 特徴 |
|---|---|---|---|
| 長柱 | $\lambda > \lambda_c$ | オイラー座屈 | 弾性座屈。材料によらず幾何形状で決まる |
| 中間柱 | $40 < \lambda < \lambda_c$ | 接線係数理論・ジョンソン式 | 非弾性座屈。材料特性が重要 |
| 短柱 | $\lambda < 40$ | 圧縮強度（降伏・圧壊） | 座屈より圧縮破壊が支配的 |
設計上のポイント
-
偏心は避けられない — 荷重の偏心は接合部の設計で最小化できますが、完全にゼロにはできません。設計規格では偏心比 $ec/r^2 = 0.25$ 程度を標準的に仮定しています
-
初期たわみの管理 — 製造公差を厳しくすることで初期たわみを低減できます。航空宇宙分野では $a_0/L < 1/1000$ が一般的な品質基準です
-
残留応力の影響 — 圧延鋼材や溶接構造には製造過程で残留応力が生じます。これは見かけ上の降伏応力を低下させ、非弾性座屈荷重に影響します。AISC規格の柱曲線は残留応力の影響を含んでいます
-
二軸曲げ座屈 — 荷重が2方向に偏心している場合、弱軸と強軸の両方で座屈を検討する必要があります。通常は弱軸（最小断面二次モーメントの軸）まわりの座屈が先に発生します
-
動的座屈 — 衝撃荷重や振動荷重下での座屈は、静的理論とは異なる挙動を示します。応力波の伝播速度と座屈の発展速度の関係が重要になります
Pythonによる可視化 (6): 総合的な座屈解析ダッシュボード
最後に、異なる材料と境界条件の組み合わせに対する柱の座屈特性を一覧できるダッシュボードを作成します。
import numpy as np
import matplotlib.pyplot as plt
# --- 材料データ ---
materials = {
'Steel (S235)': {'E': 200e9, 'sigma_y': 235e6, 'color': '#1f77b4'},
'Steel (S355)': {'E': 200e9, 'sigma_y': 355e6, 'color': '#ff7f0e'},
'Aluminum 6061': {'E': 69e9, 'sigma_y': 276e6, 'color': '#2ca02c'},
'Titanium Ti-6Al':{'E': 114e9, 'sigma_y': 880e6, 'color': '#d62728'},
}
K_vals = {'Pin-Pin (K=1.0)': 1.0, 'Fixed-Free (K=2.0)': 2.0,
'Fixed-Fixed (K=0.5)': 0.5, 'Fixed-Pin (K=0.7)': 0.7}
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
# --- (a) 材料比較（オイラー + ジョンソン）---
ax = axes[0, 0]
lam = np.linspace(5, 250, 500)
for name, mat in materials.items():
E_m, sy = mat['E'], mat['sigma_y']
sigma_e = np.pi**2 * E_m / lam**2
lam_c = np.pi * np.sqrt(E_m / sy)
sigma_j = np.where(lam <= lam_c,
sy * (1 - sy / (4*np.pi**2*E_m) * lam**2), sigma_e)
ax.plot(lam, sigma_j / 1e6, lw=2.5, color=mat['color'], label=name)
ax.axhline(sy / 1e6, color=mat['color'], ls=':', lw=1, alpha=0.4)
ax.set_xlabel('Slenderness ratio $\\lambda$')
ax.set_ylabel('$\\sigma_{cr}$ [MPa]')
ax.set_title('(a) Material comparison (Johnson + Euler)')
ax.legend(fontsize=9); ax.grid(True, alpha=0.3)
ax.set_xlim(0, 250); ax.set_ylim(0, 1000)
# --- (b) 境界条件比較（鋼S235）---
ax = axes[0, 1]
E_s, sy_s = 200e9, 235e6
L_phys = np.linspace(0.5, 10, 500) # 物理長さ [m]
r_s = 0.04 # 断面二次半径 [m]
for label, K in K_vals.items():
lam_eff = K * L_phys / r_s
sigma_e = np.pi**2 * E_s / lam_eff**2
lam_c = np.pi * np.sqrt(E_s / sy_s)
sigma_j = np.where(lam_eff <= lam_c,
sy_s * (1 - sy_s / (4*np.pi**2*E_s) * lam_eff**2), sigma_e)
P_cr = sigma_j * np.pi * r_s**2 * 4 # 仮定: 円形断面 A = 4*pi*r^2 程度
ax.plot(L_phys, sigma_j / 1e6, lw=2.5, label=label)
ax.axhline(sy_s / 1e6, color='k', ls=':', lw=1.5, alpha=0.5)
ax.set_xlabel('Physical column length $L$ [m]')
ax.set_ylabel('$\\sigma_{cr}$ [MPa]')
ax.set_title('(b) Boundary condition comparison (Steel S235, r=40mm)')
ax.legend(fontsize=9); ax.grid(True, alpha=0.3)
ax.set_ylim(0, 300)
# --- (c) 安全率の比較 ---
ax = axes[1, 0]
lam = np.linspace(10, 200, 500)
sigma_e = np.pi**2 * E_s / lam**2
# Johnson (理論)
lam_c = np.pi * np.sqrt(E_s / sy_s)
sigma_j = np.where(lam <= lam_c,
sy_s * (1 - sy_s / (4*np.pi**2*E_s) * lam**2), sigma_e)
# AISC
C_c = np.sqrt(2 * np.pi**2 * E_s / sy_s)
FS_aisc = 5/3 + 3*lam/(8*C_c) - lam**3/(8*C_c**3)
sigma_aisc_in = sy_s / FS_aisc * (1 - lam**2 / (2*C_c**2))
sigma_aisc_el = 12 * np.pi**2 * E_s / (23 * lam**2)
sigma_aisc = np.where(lam <= C_c, sigma_aisc_in, sigma_aisc_el)
# 安全率 = 理論臨界応力 / AISC許容応力
FS_actual = sigma_j / sigma_aisc
ax.plot(lam, FS_actual, 'b-', lw=2.5, label='FS = Johnson / AISC')
ax.axhline(23/12, color='r', ls='--', lw=1.5, label=f'Elastic FS = {23/12:.2f}')
ax.axhline(5/3, color='g', ls='--', lw=1.5, label=f'Min FS = {5/3:.2f}')
ax.set_xlabel('Slenderness ratio $\\lambda$')
ax.set_ylabel('Factor of Safety')
ax.set_title('(c) Implied factor of safety (AISC vs Johnson)')
ax.legend(fontsize=9); ax.grid(True, alpha=0.3)
ax.set_xlim(0, 200); ax.set_ylim(1, 3)
# --- (d) 座屈荷重の正規化比較 ---
ax = axes[1, 1]
P_ratio_range = np.linspace(0.01, 0.99, 300)
# 偏心荷重（セカント公式ベースの正規化たわみ）
ec_r2_vals = [0.1, 0.5]
for ec in ec_r2_vals:
v_secant = ec * (1/np.cos(np.pi/2 * np.sqrt(P_ratio_range)) - 1)
v_secant = np.minimum(v_secant, 20) # 表示上の上限
ax.plot(v_secant, P_ratio_range, lw=2,
label=f'Eccentric ($ec/r^2$={ec})', ls='-')
# 初期たわみ
for a0_r in [0.1, 0.5]:
v_imperfect = a0_r / (1 - P_ratio_range)
v_imperfect = np.minimum(v_imperfect, 20)
ax.plot(v_imperfect, P_ratio_range, lw=2,
label=f'Initial imperf. ($a_0 c/r^2$={a0_r})', ls='--')
ax.axhline(1.0, color='k', ls=':', lw=1.5, alpha=0.5)
ax.set_xlabel('Normalized max deflection')
ax.set_ylabel('$P / P_{cr}$')
ax.set_title('(d) Eccentricity vs Initial imperfection')
ax.legend(fontsize=8); ax.grid(True, alpha=0.3)
ax.set_xlim(0, 10); ax.set_ylim(0, 1.1)
plt.tight_layout()
plt.show()
この総合ダッシュボードから得られる知見を整理します。
- 図(a): 材料比較 — チタン合金（Ti-6Al-4V）は降伏応力が880 MPaと非常に高いため、広い細長比範囲で非弾性座屈が支配的です。一方、アルミ合金（6061-T6）はヤング率が鋼の約1/3であるため、同じ細長比でもオイラー応力が低く、弾性座屈が起きやすくなります。材料選定では強度だけでなく剛性（ヤング率）も座屈性能に直結することが確認できます。
- 図(b): 境界条件比較 — 同じ物理的な柱長でも、端末条件によって臨界応力が大きく異なります。両端固定（$K = 0.5$）と一端固定一端自由（$K = 2.0$）では臨界応力に最大16倍の差が生じます。これは有効座屈長さの比 $(2.0/0.5)^2 = 16$ に対応しています。
- 図(c): 安全率 — AISC規格に暗黙的に含まれる安全率は細長比に依存し、$\lambda = 0$ 付近で約5/3（1.67）、弾性領域で約23/12（1.92）です。中間細長比域ではより大きな安全率が取られており、この領域での不確定性（残留応力、初期不整、偏心）の影響が大きいことを反映しています。
- 図(d): 偏心と初期たわみの比較 — セカント公式（偏心荷重）と増幅係数モデル（初期たわみ）の荷重-たわみ曲線を比較しています。同じパラメータ値（$ec/r^2$ と $a_0 c/r^2$）に対して、偏心荷重の方がわずかに大きなたわみを生じますが、定性的な挙動は非常に類似しています。これが等価偏心の概念の根拠です。
まとめ
本記事では、オイラーの座屈荷重を出発点として、現実の柱が直面する3つの主要な課題 — 偏心荷重、初期たわみ、非弾性座屈 — を解説しました。
- 偏心荷重: セカント公式 $\sigma_y = (P/A)[1 + (ec/r^2)\sec((L/2r)\sqrt{P/(EA)})]$ により、偏心の影響を定量的に評価できます。偏心が大きいほど、より低い荷重で降伏に達します
- 初期たわみ: 荷重増幅係数 $1/(1 – P/P_{cr})$ により、初期たわみが圧縮荷重で指数関数的に増幅されることが示されました。$P/P_{cr} > 0.8$ の領域では特に急激な増幅が起こります
- 非弾性座屈: 接線係数理論 $\sigma_{cr} = \pi^2 E_t / \lambda^2$ と二重係数理論は、短〜中間柱での座屈荷重を予測します。Shanleyの解析により、接線係数荷重が実用的な下界であることが確立されました
- 設計曲線: ジョンソンの放物線公式、Perry-Robertson公式、AISC規格の柱曲線は、これらの理論を設計実務に落とし込んだものです
- 境界条件: 有効座屈長さ $L_e = KL$ により、異なる端末条件を統一的に取り扱えます。拘束が強いほど座屈荷重は大きくなります
座屈問題は構造設計の中核をなすテーマであり、本記事で扱った内容は実務における判断の基礎となります。特に「不完全性に対する感度」と「弾性/非弾性の遷移」の理解は、安全で合理的な設計を行う上で不可欠です。
次のステップとして、以下の記事も参考にしてください。