はりのたわみを求めるとき、弾性曲線方程式 $EIy” = M(x)$ を2回積分して境界条件を適用する方法は確実ですが、知りたいのが特定の1点のたわみだけという場面では少し回り道に感じないでしょうか。「梁の中央のたわみだけ知りたい」「トラスの特定の節点が何mm動くか知りたい」——こうしたピンポイントの問題に対して、積分を1回で済ませ、たわみ曲線全体を求めることなく答えを得られる強力な手法がカスティリアーノの定理（Castigliano’s theorem） です。
カスティリアーノの定理は、構造物の中に蓄えられたひずみエネルギー（strain energy） を荷重で偏微分するだけで変位が得られるという、驚くほどエレガントな定理です。この手法を理解すると、以下のような幅広い場面で活用できます。
- はりのたわみ計算: 集中荷重・分布荷重が作用するはりの任意の点の変位を効率的に求められる
- トラスの変位解析: 多数の部材からなるトラス構造の節点変位を系統的に計算できる
- 不静定構造の解法: 不静定問題に対して適合条件をエネルギー法で表現できる
- 曲がりはり・コイルばねの変形: 曲げ・せん断・軸力が複合する構造の変形を統一的に扱える
本記事の内容
- ひずみエネルギーの概念と各荷重モード（軸力・曲げ・せん断・ねじり）での計算式
- カスティリアーノの第1定理と第2定理の導出
- 仮想荷重法（ダミーロード法）の考え方と使い方
- はりのたわみ計算への応用（手計算 + Python）
- トラス構造の変位計算への応用
- 相反定理（マクスウェル・ベッティの定理）
- Pythonによる数値計算と理論の検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ひずみエネルギーとは
力学的エネルギーの視点
ばねを手で引き伸ばすと、手が加えた仕事はばねの中に蓄えられます。手を離せばばねは元に戻り、蓄えられたエネルギーが運動エネルギーに変換されます。これはエネルギー保存則の具体例です。
弾性体も同様です。外力が構造物に仕事をすると、そのエネルギーは構造物内部にひずみエネルギー（strain energy） として蓄えられます。弾性範囲内であれば、荷重を取り除くとひずみエネルギーはゼロに戻り、構造物は元の形に復元します。大雑把に言えば、ひずみエネルギーとは「構造物という巨大なばねに蓄えられた弾性エネルギー」です。
ひずみエネルギーの一般的定義
微小体積要素 $dV$ における応力 $\sigma$ とひずみ $\varepsilon$ を考えます。この微小要素に蓄えられるひずみエネルギー密度（単位体積あたりのエネルギー）$u_0$ は、応力-ひずみ曲線の下の面積に等しくなります。
$$ u_0 = \int_0^{\varepsilon} \sigma \, d\varepsilon $$
フックの法則が成り立つ線形弾性体では $\sigma = E\varepsilon$ なので、これを代入すると:
$$ u_0 = \int_0^{\varepsilon} E\varepsilon \, d\varepsilon = \frac{1}{2}E\varepsilon^2 = \frac{\sigma^2}{2E} $$
構造物全体のひずみエネルギー $U$ は、この密度を体積全体で積分したものです:
$$ \boxed{U = \int_V u_0 \, dV = \int_V \frac{\sigma^2}{2E} \, dV} $$
この式は一般的なものですが、実際の計算では荷重のモード（軸力・曲げ・せん断・ねじり）ごとに応力分布を代入した具体的な公式を用います。
軸力によるひずみエネルギー
断面積 $A$、長さ $L$ の棒に軸力 $N$ が作用するとき、一様応力 $\sigma = N/A$ を代入します:
$$ U_{\text{axial}} = \int_0^L \frac{(N/A)^2}{2E} A \, dx = \int_0^L \frac{N^2}{2EA} \, dx $$
$N$ と $A$ が一定の場合:
$$ \boxed{U_{\text{axial}} = \frac{N^2 L}{2EA}} $$
これは「ばねのエネルギー $\frac{1}{2}k\delta^2$」と同じ構造です。$k = EA/L$（棒の軸剛性）、$\delta = NL/(EA)$（伸び）と見なせば、$U = \frac{1}{2}(EA/L)(NL/(EA))^2 = N^2L/(2EA)$ と一致します。
曲げによるひずみエネルギー
曲げ応力 $\sigma = My/I$ を代入します。断面上の積分で $\int_A y^2 \, dA = I$（断面二次モーメント）を使うと:
$$ U_{\text{bending}} = \int_0^L \int_A \frac{(My/I)^2}{2E} \, dA \, dx = \int_0^L \frac{M^2}{2EI} \, dx $$
$$ \boxed{U_{\text{bending}} = \int_0^L \frac{M^2}{2EI} \, dx} $$
はりの問題では、この曲げのひずみエネルギーが支配的です。SFDとBMDから $M(x)$ を求め、この積分を実行することで $U$ が計算できます。
せん断力によるひずみエネルギー
せん断力 $V$ による寄与は:
$$ U_{\text{shear}} = \int_0^L \frac{\kappa V^2}{2GA} \, dx $$
ここで $G$ はせん断弾性係数、$\kappa$ は断面形状に依存する形状係数（矩形断面で $6/5$、円形断面で $10/9$）です。通常、はりでは曲げの寄与に比べてせん断の寄与は小さいため、細長いはり（$L/h > 10$ 程度）では無視できます。
ねじりによるひずみエネルギー
ねじりモーメント $T$ が作用する円形断面棒では:
$$ \boxed{U_{\text{torsion}} = \int_0^L \frac{T^2}{2GJ} \, dx} $$
ここで $J$ は断面二次極モーメントです。
ひずみエネルギーの加法性
線形弾性体では、複数の荷重モードが同時に作用する場合、全ひずみエネルギーはそれぞれの寄与の和として計算できます:
$$ U = U_{\text{axial}} + U_{\text{bending}} + U_{\text{shear}} + U_{\text{torsion}} $$
ひずみエネルギーの計算方法がわかったところで、いよいよ本題に入ります。このひずみエネルギーと外力・変位の関係を明確にしたのがカスティリアーノの定理です。
カスティリアーノの第2定理
定理の直感的な意味
イメージとしては、構造物を巨大な非線形ばねと考えてください。ばねに蓄えられたエネルギーが荷重の大きさに依存するとき、「エネルギーを荷重で微分すると変位が出る」——これがカスティリアーノの定理の本質です。
ばね定数 $k$ の線形ばねで考えると、ひずみエネルギーは $U = \frac{1}{2}kx^2$ です。荷重 $F = kx$ を使って $x = F/k$ と書き換えると $U = F^2/(2k)$ となり、$dU/dF = F/k = x$ が得られます。つまりエネルギーを力で微分すると変位になります。カスティリアーノの定理はこの関係を、複数の荷重が作用する一般的な弾性体に拡張したものです。
定理の正確な記述
カスティリアーノの第2定理（Castigliano’s second theorem）: 線形弾性体に $n$ 個の独立な外力 $P_1, P_2, \dots, P_n$ が作用しているとき、ひずみエネルギー $U$ を力 $P_i$ の関数として表すと、$P_i$ の作用点における $P_i$ 方向の変位 $\delta_i$ は:
$$ \boxed{\delta_i = \frac{\partial U}{\partial P_i}} $$
同様に、モーメント $M_i$ に対しては、その作用点での回転角 $\theta_i$ が:
$$ \boxed{\theta_i = \frac{\partial U}{\partial M_i}} $$
で与えられます。ここで $U$ はひずみエネルギーを外力の関数として書いたもの（相補ひずみエネルギー）です。線形弾性体ではひずみエネルギーと相補ひずみエネルギーが等しいため、区別なく $U$ と書けます。
定理の導出
$n$ 個の外力 $P_1, P_2, \dots, P_n$ が作用する線形弾性体を考えます。各力の作用点における力の方向の変位を $\delta_1, \delta_2, \dots, \delta_n$ とします。
外力が行う仕事は、荷重がゼロから最終値まで比例的に増加する場合:
$$ W = \frac{1}{2}\sum_{i=1}^{n} P_i \delta_i $$
エネルギー保存則より、外力の仕事はすべてひずみエネルギーに変換されるので:
$$ U = W = \frac{1}{2}\sum_{i=1}^{n} P_i \delta_i $$
線形弾性体では、重ね合わせの原理により変位は荷重の線形関数です。すなわち柔軟性係数 $f_{ij}$ を用いて:
$$ \delta_i = \sum_{j=1}^{n} f_{ij} P_j $$
と書けます。これを $U$ に代入すると:
$$ U = \frac{1}{2}\sum_{i=1}^{n} P_i \sum_{j=1}^{n} f_{ij} P_j = \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n} f_{ij} P_i P_j $$
$U$ を $P_k$ で偏微分します。$P_i P_j$ を $P_k$ で微分するとき、$i = k$ の項と $j = k$ の項が寄与します:
$$ \frac{\partial U}{\partial P_k} = \frac{1}{2}\sum_{j=1}^{n} f_{kj} P_j + \frac{1}{2}\sum_{i=1}^{n} f_{ik} P_i $$
ここで、後述する相反定理 $f_{ij} = f_{ji}$ を用いると、右辺の2つの和は等しくなります:
$$ \frac{\partial U}{\partial P_k} = \sum_{j=1}^{n} f_{kj} P_j = \delta_k $$
これがカスティリアーノの第2定理です。この導出のポイントは、ひずみエネルギーが荷重の2次形式（二次関数）になること、そして柔軟性行列の対称性（$f_{ij} = f_{ji}$）を使うことです。
カスティリアーノの第1定理
第2定理は「ひずみエネルギーを荷重で微分すると変位」でした。第1定理はその逆の関係を述べます。
カスティリアーノの第1定理: ひずみエネルギーを変位 $\delta_i$ の関数として表したとき、$\delta_i$ で偏微分すると対応する荷重 $P_i$ が得られます:
$$ \boxed{P_i = \frac{\partial U}{\partial \delta_i}} $$
導出は第2定理と同様です。ひずみエネルギーを変位で表すと:
$$ U = \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n} k_{ij} \delta_i \delta_j $$
ここで $k_{ij}$ は剛性係数です。$\delta_k$ で偏微分すると:
$$ \frac{\partial U}{\partial \delta_k} = \sum_{j=1}^{n} k_{kj} \delta_j = P_k $$
第1定理は有限要素法（FEM）の定式化の基盤となる重要な定理ですが、実用的な手計算では第2定理を圧倒的に多く使います。以降、「カスティリアーノの定理」と言えば第2定理を指すことにします。
カスティリアーノの定理の記述と導出がわかったところで、次は「求めたい点に外力が作用していない場合」にどう対処するかを見ていきましょう。
仮想荷重法（ダミーロード法）
問題の設定
カスティリアーノの定理 $\delta_i = \partial U / \partial P_i$ を使うには、変位を求めたい点に外力 $P_i$ が作用している必要があります。しかし実際には、「荷重が作用していない点のたわみを知りたい」という場面が頻繁にあります。
例えば、等分布荷重だけが作用する単純支持はりの中央のたわみを求めたい場合、中央には集中荷重がありません。このままではカスティリアーノの定理を直接適用できません。
ダミーロードのアイデア
解決策は巧妙です。求めたい点・方向に仮想的な荷重（ダミーロード） $Q$ を追加し、ひずみエネルギー $U$ を $Q$ の関数として計算します。カスティリアーノの定理で $\partial U / \partial Q$ を求めた後、最終的に $Q = 0$ とすれば、元の問題の変位が得られます:
$$ \boxed{\delta = \left.\frac{\partial U}{\partial Q}\right|_{Q=0}} $$
この手法を仮想荷重法またはダミーロード法と呼びます。
実用的な計算手順
はりの曲げ問題を例に、具体的な手順を示します。
ステップ1: 変位を求めたい点に、求めたい方向のダミーロード $Q$ を追加する。
ステップ2: $Q$ を含めた曲げモーメント $M(x)$ を求める。$M$ は実荷重と $Q$ の両方の関数になる。
ステップ3: ひずみエネルギーの偏微分を、積分の中で行う:
$$ \delta = \frac{\partial U}{\partial Q} = \frac{\partial}{\partial Q}\int_0^L \frac{M^2}{2EI}dx = \int_0^L \frac{M}{EI}\frac{\partial M}{\partial Q}dx $$
ここで偏微分と積分の順序を交換しました。この形が計算上非常に便利です。
ステップ4: $Q = 0$ を代入して最終結果を得る。
偏微分を積分の中に持ち込むテクニックにより、$M^2$ の積分ではなく $M \cdot (\partial M / \partial Q)$ の積分を計算すればよくなり、計算量が大幅に減ります。
それでは、この手法を実際のはりの問題に適用してみましょう。
はりのたわみ計算への応用
例1: 片持ちはり — 先端集中荷重
固定端を $x = 0$、自由端を $x = L$ とし、自由端に下向きの荷重 $P$ が作用する片持ちはりを考えます。
曲げモーメント: 自由端から測ると単純に $M(x) = P(x – L)$ です（固定端から $x$ の位置で切断し、右側の釣り合いを考える）。ただし、ここでは自由端のたわみを求めたいので、荷重 $P$ そのものについて偏微分します。
$$ M(x) = P(x – L) $$
$$ \frac{\partial M}{\partial P} = x – L $$
カスティリアーノの定理を適用します:
$$ \delta = \int_0^L \frac{M}{EI}\frac{\partial M}{\partial P}dx = \int_0^L \frac{P(x-L)^2}{EI}dx $$
$(x – L)^2$ を展開して積分します:
$$ \delta = \frac{P}{EI}\int_0^L (x-L)^2 dx $$
$t = x – L$ と置換すると $t: -L \to 0$ で:
$$ \delta = \frac{P}{EI}\int_{-L}^{0} t^2 dt = \frac{P}{EI}\left[\frac{t^3}{3}\right]_{-L}^{0} = \frac{P}{EI}\cdot\frac{L^3}{3} $$
$$ \boxed{\delta = \frac{PL^3}{3EI}} $$
これは梁のたわみの弾性曲線方程式から得られる結果 $\delta_{\max} = PL^3/(3EI)$ と完全に一致します。カスティリアーノの定理では2回積分ではなく1回の積分で結果が得られたことに注目してください。
例2: 単純支持はり — 中央集中荷重
長さ $L$ の単純支持はりの中央（$x = L/2$）に下向き荷重 $P$ が作用する場合の、荷重作用点のたわみを求めます。
反力は対称性から $R_A = R_B = P/2$ です。
曲げモーメント: 対称性を利用して $0 \le x \le L/2$ の区間で計算し、結果を2倍します。
$$ M(x) = \frac{P}{2}x \quad (0 \le x \le L/2) $$
$$ \frac{\partial M}{\partial P} = \frac{x}{2} $$
対称性を用いて積分区間を半分にし、2倍します:
$$ \delta = 2\int_0^{L/2} \frac{M}{EI}\frac{\partial M}{\partial P}dx = 2\int_0^{L/2} \frac{1}{EI}\cdot\frac{Px}{2}\cdot\frac{x}{2}dx $$
$$ = \frac{P}{2EI}\int_0^{L/2} x^2 dx = \frac{P}{2EI}\cdot\frac{1}{3}\left(\frac{L}{2}\right)^3 = \frac{P}{2EI}\cdot\frac{L^3}{24} $$
$$ \boxed{\delta = \frac{PL^3}{48EI}} $$
これも教科書の公式と一致します。
例3: 単純支持はり — 等分布荷重（ダミーロード法）
等分布荷重 $w$ [N/m] が作用する単純支持はりの中央たわみを求めます。中央には集中荷重がないので、ダミーロード法を使います。
ステップ1: 中央（$x = L/2$）に下向きのダミーロード $Q$ を追加。
ステップ2: 反力は $R_A = wL/2 + Q/2$（対称性より）。
$0 \le x \le L/2$ の区間での曲げモーメント:
$$ M(x) = \left(\frac{wL}{2} + \frac{Q}{2}\right)x – \frac{wx^2}{2} $$
$Q$ による偏微分:
$$ \frac{\partial M}{\partial Q} = \frac{x}{2} $$
ステップ3: カスティリアーノの定理を適用し、対称性で2倍します:
$$ \delta = 2\int_0^{L/2} \frac{M}{EI}\frac{\partial M}{\partial Q}dx $$
ステップ4: $Q = 0$ を代入すると:
$$ \delta = 2\int_0^{L/2} \frac{1}{EI}\left(\frac{wLx}{2} – \frac{wx^2}{2}\right)\cdot\frac{x}{2}dx $$
$$ = \frac{w}{EI}\int_0^{L/2} \left(\frac{Lx^2}{2} – \frac{x^3}{2}\right)dx $$
各項を積分します:
$$ \int_0^{L/2} \frac{Lx^2}{2}dx = \frac{L}{2}\cdot\frac{1}{3}\left(\frac{L}{2}\right)^3 = \frac{L^4}{48} $$
$$ \int_0^{L/2} \frac{x^3}{2}dx = \frac{1}{2}\cdot\frac{1}{4}\left(\frac{L}{2}\right)^4 = \frac{L^4}{128} $$
$$ \delta = \frac{w}{EI}\left(\frac{L^4}{48} – \frac{L^4}{128}\right) = \frac{w}{EI}\cdot\frac{L^4(128 – 48\cdot 1)}{48 \cdot 128} $$
通分して計算すると:
$$ \frac{L^4}{48} – \frac{L^4}{128} = L^4\left(\frac{8 – 3}{384}\right) = \frac{5L^4}{384} $$
$$ \boxed{\delta = \frac{5wL^4}{384EI}} $$
弾性曲線方程式から求めた公式と完全に一致します。ダミーロード法が正しく機能していることが確認できました。
ここまでで、はりの曲げ問題にカスティリアーノの定理を適用する方法を見てきました。次に、これをPythonで実装して理論値との一致を数値的にも検証しましょう。
Pythonでの実装1: ひずみエネルギーの計算と可視化
まずは、荷重の増加に伴うひずみエネルギーの変化を可視化し、エネルギーの偏微分が変位に等しいことを数値的に確認します。
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ ---
E = 200e9 # ヤング率 [Pa] (鋼)
I = 8.33e-6 # 断面二次モーメント [m^4] (100mm x 100mm 矩形断面)
L = 2.0 # はりの長さ [m]
# --- 片持ちはり先端集中荷重のひずみエネルギー ---
# U = P^2 L^3 / (6EI)
P_arr = np.linspace(0, 50000, 500) # 荷重 [N]
U_arr = P_arr**2 * L**3 / (6 * E * I)
# --- dU/dP の数値微分 ---
dP = P_arr[1] - P_arr[0]
dUdP_numerical = np.gradient(U_arr, dP)
# --- 理論値: δ = PL^3 / (3EI) ---
delta_theory = P_arr * L**3 / (3 * E * I)
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: ひずみエネルギー
axes[0].plot(P_arr / 1000, U_arr, 'b-', linewidth=2)
axes[0].set_xlabel('Load P [kN]', fontsize=12)
axes[0].set_ylabel('Strain Energy U [J]', fontsize=12)
axes[0].set_title('Strain Energy of Cantilever Beam', fontsize=13)
axes[0].grid(True, alpha=0.3)
# 右: dU/dP vs 理論たわみ
axes[1].plot(P_arr / 1000, dUdP_numerical * 1000, 'ro', markersize=2,
label='dU/dP (numerical)')
axes[1].plot(P_arr / 1000, delta_theory * 1000, 'b-', linewidth=2,
label=r'$\delta = PL^3/(3EI)$ (theory)')
axes[1].set_xlabel('Load P [kN]', fontsize=12)
axes[1].set_ylabel('Deflection [mm]', fontsize=12)
axes[1].set_title("Castigliano's Theorem Verification", fontsize=13)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のグラフから、ひずみエネルギーが荷重の2乗に比例して増加する放物線を描くことがわかります。これは $U = P^2L^3/(6EI)$ という理論式と一致しています。
右のグラフでは、ひずみエネルギーの数値微分 $dU/dP$（赤丸）と理論的なたわみ $\delta = PL^3/(3EI)$（青線）が完全に重なっています。これはカスティリアーノの第2定理 $\delta = \partial U / \partial P$ が数値的にも正しいことを示しています。荷重とたわみが線形関係にあること（フックの法則の帰結）も読み取れます。
Pythonでの実装2: はりのたわみをカスティリアーノの定理で解く
次に、カスティリアーノの定理を使ったたわみ計算をPythonで実装します。SymPyの記号計算を利用して、手計算の過程をそのままコード化します。
import sympy as sp
# 記号の定義
x, P, w, L, E_sym, I_sym, Q = sp.symbols('x P w L E I Q', positive=True)
print("=" * 60)
print("例1: 片持ちはり — 先端集中荷重")
print("=" * 60)
# 曲げモーメント: M(x) = P(x - L)
M1 = P * (x - L)
dM1_dP = sp.diff(M1, P)
# カスティリアーノの定理
integrand1 = M1 * dM1_dP / (E_sym * I_sym)
delta1 = sp.integrate(integrand1, (x, 0, L))
delta1_simplified = sp.simplify(delta1)
print(f"M(x) = {M1}")
print(f"dM/dP = {dM1_dP}")
print(f"たわみ δ = {delta1_simplified}")
print()
print("=" * 60)
print("例2: 単純支持はり — 中央集中荷重")
print("=" * 60)
# 0 <= x <= L/2 の区間（対称性で2倍）
M2 = P * x / 2
dM2_dP = sp.diff(M2, P)
integrand2 = M2 * dM2_dP / (E_sym * I_sym)
delta2 = 2 * sp.integrate(integrand2, (x, 0, L / 2))
delta2_simplified = sp.simplify(delta2)
print(f"M(x) = {M2} (0 <= x <= L/2)")
print(f"dM/dP = {dM2_dP}")
print(f"たわみ δ = {delta2_simplified}")
print()
print("=" * 60)
print("例3: 単純支持はり — 等分布荷重（ダミーロード法）")
print("=" * 60)
# 0 <= x <= L/2、ダミーロード Q を中央に追加
M3 = (w * L / 2 + Q / 2) * x - w * x**2 / 2
dM3_dQ = sp.diff(M3, Q)
# Q=0 を代入してから積分
M3_Q0 = M3.subs(Q, 0)
integrand3 = M3_Q0 * dM3_dQ / (E_sym * I_sym)
delta3 = 2 * sp.integrate(integrand3, (x, 0, L / 2))
delta3_simplified = sp.simplify(delta3)
print(f"M(x) = {M3}")
print(f"dM/dQ = {dM3_dQ}")
print(f"M(x)|Q=0 = {M3_Q0}")
print(f"たわみ δ = {delta3_simplified}")
このコードの出力から、3つの例すべてで手計算と一致する結果が得られることが確認できます:
- 例1: $\delta = PL^3/(3EI)$（片持ちはり先端たわみ）
- 例2: $\delta = PL^3/(48EI)$（単純支持はり中央たわみ）
- 例3: $\delta = 5wL^4/(384EI)$（等分布荷重の単純支持はり中央たわみ）
SymPyを使うことで、偏微分と積分を記号的に実行でき、手計算のミスを防ぎながら結果を検証できます。特に例3のダミーロード法では、$Q$ を含む $M(x)$ を定義し、$Q$ で微分した後に $Q = 0$ を代入するという手順がコードで明確に表現されています。
次は、カスティリアーノの定理をトラス構造の変位計算に適用する方法を解説します。トラスでは各部材が軸力のみを受けるため、曲げのひずみエネルギーではなく軸力のひずみエネルギーを使います。
トラス構造の変位計算
トラスにおけるカスティリアーノの定理
トラスの各部材は軸力のみを受けるため、全ひずみエネルギーは:
$$ U = \sum_{i=1}^{m} \frac{N_i^2 L_i}{2E_i A_i} $$
ここで $m$ は部材数、$N_i$, $L_i$, $E_i$, $A_i$ はそれぞれ第 $i$ 部材の軸力、長さ、ヤング率、断面積です。
カスティリアーノの第2定理を適用すると、外力 $P_j$ の方向の節点変位は:
$$ \delta_j = \frac{\partial U}{\partial P_j} = \sum_{i=1}^{m} \frac{N_i}{E_i A_i}\frac{\partial N_i}{\partial P_j}L_i $$
$\partial N_i / \partial P_j$ は、荷重 $P_j$ が単位量変化したときの第 $i$ 部材の軸力変化であり、実質的には $P_j = 1$ としたときの部材力に等しくなります。
計算手順
- 全外力（実荷重 + 必要に応じてダミーロード）を作用させた状態で、各部材の軸力 $N_i$ を節点法や切断法で求める
- ひずみエネルギーの偏微分 $\partial N_i / \partial P_j$ を計算する
- $\delta_j = \sum N_i (\partial N_i / \partial P_j) L_i / (E_i A_i)$ を計算する
- ダミーロードを使った場合は $Q = 0$ を代入する
例: 3部材トラス
3本の部材からなる単純なトラスを考えます。節点Aに水平荷重 $P$ が作用し、節点Aの水平変位を求めます。
三角形トラスを以下のように設定します: – 節点B（左下、固定）、節点C（右下、固定）、節点A（上部、荷重点） – 部材BA（長さ $L_1$、角度 $\alpha$）、部材CA（長さ $L_2$、角度 $\beta$）、部材BC（長さ $L_3$、水平） – 全部材の $EA$ は同一
この問題の具体的な数値計算はPythonで行います。
ここまでで、カスティリアーノの定理のはりとトラスへの応用を解説しました。次に、エネルギー法の重要な帰結である相反定理を紹介します。
相反定理（マクスウェル・ベッティの定理）
直感的な理解
2つの荷重系 A と B を考えます。荷重系Aが引き起こす変位に荷重系Bが行う仕事と、荷重系Bが引き起こす変位に荷重系Aが行う仕事は等しい——これがベッティの相反定理です。
さらに特殊な場合として、点1に単位荷重を加えたときの点2の変位 $f_{21}$ と、点2に単位荷重を加えたときの点1の変位 $f_{12}$ は等しい——これがマクスウェルの相反定理です。
$$ \boxed{f_{12} = f_{21}} $$
数学的な表現
ベッティの相反定理（Betti’s reciprocal theorem）:
線形弾性体に2つの荷重系 $\{P_i^{(A)}\}$ と $\{P_i^{(B)}\}$ を考えます。荷重系Aによる変位を $\{\delta_i^{(A)}\}$、荷重系Bによる変位を $\{\delta_i^{(B)}\}$ とすると:
$$ \boxed{\sum_i P_i^{(A)} \delta_i^{(B)} = \sum_i P_i^{(B)} \delta_i^{(A)}} $$
左辺は「荷重系Bが生じさせた変位に、荷重系Aの力が行う仮想仕事」であり、右辺は「荷重系Aが生じさせた変位に、荷重系Bの力が行う仮想仕事」です。
導出
柔軟性係数を使った表現で示します。変位は:
$$ \delta_i^{(A)} = \sum_j f_{ij} P_j^{(A)}, \quad \delta_i^{(B)} = \sum_j f_{ij} P_j^{(B)} $$
左辺:
$$ \sum_i P_i^{(A)} \delta_i^{(B)} = \sum_i P_i^{(A)} \sum_j f_{ij} P_j^{(B)} = \sum_i \sum_j f_{ij} P_i^{(A)} P_j^{(B)} $$
右辺:
$$ \sum_i P_i^{(B)} \delta_i^{(A)} = \sum_i P_i^{(B)} \sum_j f_{ij} P_j^{(A)} = \sum_i \sum_j f_{ij} P_i^{(B)} P_j^{(A)} $$
添字 $i, j$ を交換して比較すると、$f_{ij} = f_{ji}$（柔軟性行列の対称性）が成り立てば左辺と右辺は等しくなります。
柔軟性行列の対称性は、ひずみエネルギーが経路に依存しないこと（状態量であること）から証明できます。荷重 $P_1$ を先にかけてから $P_2$ をかける場合と、逆の順序でかける場合で最終的なひずみエネルギーは同じでなければなりません:
$$ U = \frac{1}{2}f_{11}P_1^2 + f_{12}P_1 P_2 + \frac{1}{2}f_{22}P_2^2 $$
$$ U = \frac{1}{2}f_{11}P_1^2 + f_{21}P_1 P_2 + \frac{1}{2}f_{22}P_2^2 $$
両者が等しいためには $f_{12} = f_{21}$ が必要です。これをマクスウェルの相反定理と呼びます。
相反定理の物理的意味
マクスウェルの相反定理 $f_{12} = f_{21}$ は、以下のように述べることができます。
点1に単位力を加えたときに点2に生じる変位は、点2に単位力を加えたときに点1に生じる変位に等しい。
この定理は驚くべきことに、2つの点の位置や荷重の方向が全く異なっていても成り立ちます。構造物の形状が複雑であっても、線形弾性であれば常に成立する普遍的な性質です。
相反定理は単に美しいだけでなく、実用上も重要です。例えば、不静定構造の解法において柔軟性マトリクスを構成する際に、$f_{12} = f_{21}$ を利用して計算量を半減できます。また、構造物の影響線（influence line）の導出にも利用されます。
それでは、相反定理の成立をPythonで数値的に検証してみましょう。
Pythonでの実装3: トラスの変位計算
3部材の平面トラスについて、カスティリアーノの定理を適用して節点変位を計算します。
import numpy as np
# --- 3部材トラスの設定 ---
# 節点座標
# B(0, 0) - 固定支持
# C(4, 0) - 固定支持
# A(2, 3) - 荷重作用点
nodes = {
'B': np.array([0.0, 0.0]),
'C': np.array([4.0, 0.0]),
'A': np.array([2.0, 3.0])
}
# 部材の定義 (始点, 終点)
members = [
('B', 'A'), # 部材1
('C', 'A'), # 部材2
('B', 'C'), # 部材3
]
# 材料・断面特性（全部材共通）
E = 200e9 # ヤング率 [Pa]
A = 1e-3 # 断面積 [m^2]
# 部材長さの計算
lengths = []
for (n1, n2) in members:
L_member = np.linalg.norm(nodes[n2] - nodes[n1])
lengths.append(L_member)
print(f"部材 {n1}-{n2}: L = {L_member:.4f} m")
# --- 節点Aに作用する外力 ---
Px = 10000.0 # 水平荷重 [N]
Py = -20000.0 # 鉛直荷重 [N]（下向き負）
# --- 節点法による軸力の計算 ---
# 節点Aの力の釣り合い
# 各部材の方向余弦
def direction_cosines(n1, n2):
"""n1からn2への方向余弦を返す"""
d = nodes[n2] - nodes[n1]
L_val = np.linalg.norm(d)
return d / L_val
# 部材1(B->A), 部材2(C->A) の節点Aでの方向
# 節点Aから見た方向は逆向き（A->B, A->C）
cos_BA = direction_cosines('A', 'B') # A->B方向
cos_CA = direction_cosines('A', 'C') # A->C方向
# 節点Aの釣り合い: N1*(A->B方向) + N2*(A->C方向) + 外力 = 0
# [cos_x(AB) cos_x(AC)] [N1] [-Px]
# [cos_y(AB) cos_y(AC)] [N2] = [-Py]
coeff_matrix = np.array([
[cos_BA[0], cos_CA[0]],
[cos_BA[1], cos_CA[1]]
])
force_vector = np.array([-Px, -Py])
N12 = np.linalg.solve(coeff_matrix, force_vector)
N1, N2 = N12[0], N12[1] # 正: 引張、負: 圧縮
# 部材3(B-C)の軸力: 節点Bの水平釣り合いから
cos_BA_fromB = direction_cosines('B', 'A')
N3 = -(N1 * cos_BA_fromB[0]) # 節点Bの水平方向の釣り合い
print(f"\n--- 軸力 ---")
print(f"N1 (B-A) = {N1:.2f} N")
print(f"N2 (C-A) = {N2:.2f} N")
print(f"N3 (B-C) = {N3:.2f} N")
# --- カスティリアーノの定理による変位計算 ---
# δx = Σ (Ni / (EiAi)) * (∂Ni/∂Px) * Li
# δy = Σ (Ni / (EiAi)) * (∂Ni/∂Py) * Li
# ∂Ni/∂Px, ∂Ni/∂Py を数値微分で計算
def compute_forces(Fx, Fy):
"""節点Aの外力(Fx, Fy)に対する各部材の軸力を返す"""
fv = np.array([-Fx, -Fy])
n12 = np.linalg.solve(coeff_matrix, fv)
n1_val, n2_val = n12[0], n12[1]
n3_val = -(n1_val * cos_BA_fromB[0])
return np.array([n1_val, n2_val, n3_val])
dF = 1.0 # 微小変化 [N]
N_vec = compute_forces(Px, Py)
# Pxに関する偏微分
N_vec_px_plus = compute_forces(Px + dF, Py)
dN_dPx = (N_vec_px_plus - N_vec) / dF
# Pyに関する偏微分
N_vec_py_plus = compute_forces(Px, Py + dF)
dN_dPy = (N_vec_py_plus - N_vec) / dF
# 変位の計算
lengths_arr = np.array(lengths)
delta_x = np.sum(N_vec * dN_dPx * lengths_arr / (E * A))
delta_y = np.sum(N_vec * dN_dPy * lengths_arr / (E * A))
print(f"\n--- 節点Aの変位（カスティリアーノの定理） ---")
print(f"δx = {delta_x * 1000:.4f} mm")
print(f"δy = {delta_y * 1000:.4f} mm")
print(f"|δ| = {np.sqrt(delta_x**2 + delta_y**2) * 1000:.4f} mm")
# --- ひずみエネルギーの確認 ---
U = np.sum(N_vec**2 * lengths_arr / (2 * E * A))
W = 0.5 * (Px * delta_x + Py * delta_y)
print(f"\n--- エネルギー検証 ---")
print(f"ひずみエネルギー U = {U:.6f} J")
print(f"外力の仕事 W = {W:.6f} J")
print(f"差 |U - W| = {abs(U - W):.2e} J")
このコードの出力では、以下の3点が確認できます:
- 各部材の軸力: 節点法により3部材の軸力が計算されます。引張が正、圧縮が負で表示されます。
- 節点Aの変位: カスティリアーノの定理 $\delta_j = \sum N_i (\partial N_i / \partial P_j) L_i / (EA)$ により、水平・鉛直変位がミリメートル単位で得られます。
- エネルギー保存の検証: ひずみエネルギー $U = \sum N_i^2 L_i / (2EA)$ と外力の仕事 $W = \frac{1}{2}(P_x \delta_x + P_y \delta_y)$ が一致していることが確認でき、計算の正しさが裏付けられます。
偏微分 $\partial N_i / \partial P_j$ を数値微分で求めている点は、解析的に求めることもできますが、部材数が増えた場合にはコードの一般性を保つために数値微分が便利です。
Pythonでの実装4: 相反定理の数値検証
最後に、マクスウェル・ベッティの相反定理を単純支持はりで数値的に検証します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
# --- 単純支持はりのパラメータ ---
E = 200e9 # ヤング率 [Pa]
I = 8.33e-6 # 断面二次モーメント [m^4]
L = 3.0 # はりの長さ [m]
EI = E * I
# --- 曲げモーメントの計算 ---
def moment_point_load(x, a, P, L_beam):
"""位置 a に集中荷重 P が作用する単純支持はりの曲げモーメント"""
M = np.where(
x <= a,
P * (L_beam - a) * x / L_beam,
P * a * (L_beam - x) / L_beam
)
return M
# --- たわみの計算（弾性曲線方程式の数値積分）---
def deflection_point_load(x_eval, a, P, L_beam, EI_val):
"""位置 a に集中荷重 P が作用する単純支持はりのたわみ"""
# 解析解を使用
b = L_beam - a
if x_eval <= a:
delta = P * b * x_eval / (6 * L_beam * EI_val) * (
L_beam**2 - b**2 - x_eval**2
)
else:
delta = P * a * (L_beam - x_eval) / (6 * L_beam * EI_val) * (
2 * L_beam * x_eval - x_eval**2 - a**2
)
return delta
# --- マクスウェルの相反定理の検証 ---
# 点1: x = 0.8 m, 点2: x = 2.2 m
x1, x2 = 0.8, 2.2
P_unit = 1.0 # 単位荷重
# f_12: 点1に単位荷重 → 点2のたわみ
f_12 = deflection_point_load(x2, x1, P_unit, L, EI)
# f_21: 点2に単位荷重 → 点1のたわみ
f_21 = deflection_point_load(x1, x2, P_unit, L, EI)
print("=" * 60)
print("マクスウェルの相反定理の検証")
print("=" * 60)
print(f"はり: 単純支持、L = {L} m")
print(f"点1: x = {x1} m")
print(f"点2: x = {x2} m")
print(f"\nf_12 (点1に荷重 → 点2の変位) = {f_12:.10e} m/N")
print(f"f_21 (点2に荷重 → 点1の変位) = {f_21:.10e} m/N")
print(f"差 |f_12 - f_21| = {abs(f_12 - f_21):.2e} m/N")
print(f"相対誤差 = {abs(f_12 - f_21) / abs(f_12):.2e}")
# --- ベッティの相反定理の検証 ---
# 荷重系A: 点1に P_A, 点2にゼロ
# 荷重系B: 点1にゼロ, 点2に P_B
P_A, P_B = 5000.0, 8000.0
# 荷重系Aによる変位
delta_1A = deflection_point_load(x1, x1, P_A, L, EI)
delta_2A = deflection_point_load(x2, x1, P_A, L, EI)
# 荷重系Bによる変位
delta_1B = deflection_point_load(x1, x2, P_B, L, EI)
delta_2B = deflection_point_load(x2, x2, P_B, L, EI)
# ベッティの定理: Σ P_i^(A) * δ_i^(B) = Σ P_i^(B) * δ_i^(A)
lhs = P_A * delta_1B # 荷重系Aの力 × 荷重系Bの変位
rhs = P_B * delta_2A # 荷重系Bの力 × 荷重系Aの変位
print(f"\n{'=' * 60}")
print("ベッティの相反定理の検証")
print("=" * 60)
print(f"荷重系A: P_A = {P_A} N at x = {x1} m")
print(f"荷重系B: P_B = {P_B} N at x = {x2} m")
print(f"\nΣ P^(A) δ^(B) = {lhs:.6f} J")
print(f"Σ P^(B) δ^(A) = {rhs:.6f} J")
print(f"差 = {abs(lhs - rhs):.2e} J")
# --- 複数点での相反定理の可視化 ---
x_points = np.linspace(0.1, L - 0.1, 100)
f_1j = np.array([deflection_point_load(xp, x1, P_unit, L, EI)
for xp in x_points])
f_j1 = np.array([deflection_point_load(x1, xp, P_unit, L, EI)
for xp in x_points])
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
# 上: f_1j と f_j1 の比較
axes[0].plot(x_points, f_1j * 1e6, 'b-', linewidth=2,
label=r'$f_{1j}$: load at $x_1$, deflection at $x_j$')
axes[0].plot(x_points, f_j1 * 1e6, 'r--', linewidth=2,
label=r'$f_{j1}$: load at $x_j$, deflection at $x_1$')
axes[0].axvline(x=x1, color='gray', linestyle=':', alpha=0.5,
label=f'$x_1 = {x1}$ m')
axes[0].set_xlabel('Position $x_j$ [m]', fontsize=12)
axes[0].set_ylabel('Flexibility coefficient [$\\mu$m/N]', fontsize=12)
axes[0].set_title("Maxwell's Reciprocal Theorem Verification", fontsize=13)
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)
# 下: 差のプロット
axes[1].plot(x_points, np.abs(f_1j - f_j1), 'k-', linewidth=1)
axes[1].set_xlabel('Position $x_j$ [m]', fontsize=12)
axes[1].set_ylabel('$|f_{1j} - f_{j1}|$ [m/N]', fontsize=12)
axes[1].set_title('Difference (should be zero)', fontsize=13)
axes[1].grid(True, alpha=0.3)
axes[1].ticklabel_format(axis='y', style='scientific', scilimits=(-15, -15))
plt.tight_layout()
plt.show()
上のグラフでは、$f_{1j}$（点1に荷重を置いて各点 $x_j$ のたわみを測定）と $f_{j1}$（各点 $x_j$ に荷重を置いて点1のたわみを測定）が完全に一致しています。青の実線と赤の破線が重なっていることから、マクスウェルの相反定理 $f_{1j} = f_{j1}$ がすべての点で成立していることが視覚的に確認できます。
下のグラフでは差 $|f_{1j} – f_{j1}|$ をプロットしていますが、機械精度レベル（$10^{-15}$ m/N 程度）でほぼゼロになっています。解析解を用いているため数値誤差は浮動小数点の丸め誤差のみであり、相反定理が厳密に成立することを裏付けています。
ベッティの相反定理についても、出力で $\sum P_i^{(A)} \delta_i^{(B)} = \sum P_i^{(B)} \delta_i^{(A)}$ が数値的に一致していることが確認できます。
エネルギー法の体系的な位置づけ
ここまでカスティリアーノの定理と相反定理を個別に解説してきましたが、これらはエネルギー法という大きな体系の一部です。最後に、この体系を整理しておきましょう。
仮想仕事の原理との関係
仮想仕事の原理は、釣り合い状態にある構造物に仮想的な変位を与えたとき、外力の仮想仕事と内力の仮想仕事が等しいという原理です:
$$ \delta W_{\text{ext}} = \delta W_{\text{int}} $$
カスティリアーノの定理は、この仮想仕事の原理にフックの法則（線形弾性）の条件を加えることで得られます。仮想仕事の原理はより一般的で、非線形材料にも適用できますが、カスティリアーノの定理は線形弾性体に特化した分、計算が簡単になるという利点があります。
最小ポテンシャルエネルギーの原理
カスティリアーノの第1定理は、最小ポテンシャルエネルギーの原理と密接に関連しています。この原理は:
弾性体の変位は、全ポテンシャルエネルギー $\Pi = U – W_{\text{ext}}$ を最小にするものである
と述べます。ここで $U$ はひずみエネルギー、$W_{\text{ext}}$ は外力のポテンシャルです。$\Pi$ を変位で偏微分してゼロとおく条件:
$$ \frac{\partial \Pi}{\partial \delta_i} = \frac{\partial U}{\partial \delta_i} – P_i = 0 $$
から、カスティリアーノの第1定理 $P_i = \partial U / \partial \delta_i$ が直ちに得られます。この原理は有限要素法（FEM）の理論的基盤となっており、工学における数値解析の出発点です。
最小相補エネルギーの原理
一方、カスティリアーノの第2定理は最小相補エネルギーの原理と関連しています。相補エネルギー $U^*$ を力の関数として表し、適合条件を満たす変位場を求めます。線形弾性体では $U = U^*$ であるため、両者の区別を意識する必要はほとんどありません。
各手法の使い分け
| 手法 | 得意な問題 | 適用条件 |
|---|---|---|
| 弾性曲線方程式 $EIy” = M(x)$ | たわみ曲線全体が必要な場合 | はり（曲げ） |
| カスティリアーノの定理 | 特定の点の変位だけ必要な場合 | 線形弾性体 |
| 仮想荷重法（ダミーロード） | 荷重がない点の変位を求めたい場合 | 線形弾性体 |
| 仮想仕事の原理 | 非線形問題、一般的な問題 | 任意の材料 |
| エネルギー法（最小ポテンシャル） | FEM定式化、近似解法 | 保存力系 |
はりのたわみの弾性曲線方程式と比較すると、カスティリアーノの定理は「1点の変位を効率的に求める」という場面で真価を発揮します。一方、たわみ曲線の全体像が必要な場合は弾性曲線方程式の方が適しています。問題に応じて適切な手法を選択することが、効率的な構造解析の鍵です。
不静定構造への応用
カスティリアーノの定理は、不静定問題の解法にも威力を発揮します。
基本的な考え方
不静定構造では、釣り合い条件だけでは反力を決定できません。カスティリアーノの定理を使うと、変位の適合条件をエネルギーの言葉で表現できます。
例えば、両端固定はり（1次不静定）では、余剰反力 $R$ の作用点で「変位がゼロ」という条件を:
$$ \frac{\partial U}{\partial R} = 0 $$
と書けます。これは「余剰反力 $R$ に関するひずみエネルギーの偏微分がゼロ」という条件であり、最小ひずみエネルギーの定理（カスティリアーノの最小仕事の定理） とも呼ばれます。
$$ \boxed{\frac{\partial U}{\partial R_i} = 0 \quad (\text{余剰反力 } R_i \text{ の作用点で変位がゼロのとき})} $$
この条件を使えば、余剰反力を含むひずみエネルギーの式から、代数方程式を解くだけで不静定問題が解けます。
例: 一端固定・一端単純支持はり
長さ $L$ のはりで、左端（$x = 0$）が固定端、右端（$x = L$）がローラー支持とし、等分布荷重 $w$ が作用する場合を考えます。これは1次不静定問題です。
余剰反力を右端のローラー反力 $R_B$ とすると、右端の変位がゼロという条件は:
$$ \frac{\partial U}{\partial R_B} = \int_0^L \frac{M}{EI}\frac{\partial M}{\partial R_B}dx = 0 $$
曲げモーメントは:
$$ M(x) = R_B(L – x) + \frac{w}{2}(L – x)^2 – \frac{w L^2}{2} + \frac{wLx}{1} – \frac{wx^2}{2} $$
この計算はやや複雑になるため、SymPyでの実装に委ねます。重要なのは、不静定問題がカスティリアーノの定理を通じて「偏微分をゼロとおく」という代数問題に帰着されるという点です。
Pythonでの実装5: 不静定はりの解法
カスティリアーノの定理を使って、1次不静定はり（一端固定・他端単純支持）の余剰反力とたわみ曲線を求めます。
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
# --- 記号定義 ---
x, w, L, E_sym, I_sym, R_B = sp.symbols('x w L E I R_B', positive=True)
print("=" * 60)
print("不静定はり: 一端固定・他端単純支持、等分布荷重")
print("=" * 60)
# 固定端 x=0、単純支持端 x=L
# 反力: R_A (固定端鉛直), M_A (固定端モーメント), R_B (単純支持端鉛直)
# 釣り合い: R_A + R_B = wL, M_A = R_A*L - wL^2/2 ではなくモーメントの釣り合いから
# R_Bを余剰反力として扱う
# 曲げモーメント (x=0 から測る)
# M(x) = R_B*(L-x) - w*(L-x)^2/2 ... 右端から見た曲げモーメント
# もう少し系統的に: 右端からの距離 s = L-x として
# 左端からx: M(x) = (wL - R_B)*x - w*x^2/2 - M_A
# ただし M_A は未知なので、別のアプローチ
# 右端(x=L)から距離sの断面でのモーメント (s = L - x)
s = sp.Symbol('s', positive=True)
M_s = R_B * s - w * s**2 / 2
dM_dRB = sp.diff(M_s, R_B)
print(f"M(s) = {M_s} (s: 右端からの距離)")
print(f"∂M/∂R_B = {dM_dRB}")
# カスティリアーノの定理: ∂U/∂R_B = 0 (右端の変位 = 0)
integrand = M_s * dM_dRB / (E_sym * I_sym)
condition = sp.integrate(integrand, (s, 0, L))
condition_simplified = sp.simplify(condition)
print(f"\n∂U/∂R_B = {condition_simplified} = 0")
# R_B について解く
R_B_sol = sp.solve(condition_simplified, R_B)[0]
R_B_simplified = sp.simplify(R_B_sol)
print(f"\nR_B = {R_B_simplified}")
# R_A と M_A
R_A_sol = w * L - R_B_simplified
R_A_simplified = sp.simplify(R_A_sol)
# M_A: x=0でのモーメント = R_B*L - wL^2/2
M_A_sol = R_B_simplified * L - w * L**2 / 2
M_A_simplified = sp.simplify(M_A_sol)
print(f"R_A = {R_A_simplified}")
print(f"M_A = {M_A_simplified}")
# --- 数値計算と可視化 ---
L_val = 3.0 # [m]
w_val = 10000.0 # [N/m]
E_val = 200e9 # [Pa]
I_val = 8.33e-6 # [m^4]
# 数値化
R_B_num = float(R_B_simplified.subs([(w, w_val), (L, L_val)]))
R_A_num = float(R_A_simplified.subs([(w, w_val), (L, L_val)]))
M_A_num = float(M_A_simplified.subs([(w, w_val), (L, L_val)]))
print(f"\n--- 数値結果 (w={w_val} N/m, L={L_val} m) ---")
print(f"R_B = {R_B_num:.2f} N")
print(f"R_A = {R_A_num:.2f} N")
print(f"M_A = {M_A_num:.2f} N·m")
# たわみ曲線の計算（数値積分）
x_arr = np.linspace(0, L_val, 200)
# M(x) = R_A * x + M_A - w*x^2/2 (x: 左端からの距離)
M_arr = R_A_num * x_arr + M_A_num - w_val * x_arr**2 / 2
# 数値的に2回積分してたわみを求める
dx = x_arr[1] - x_arr[0]
slope = np.cumsum(M_arr * dx) / (E_val * I_val) # y'の近似
deflection = np.cumsum(slope * dx) # yの近似
# 境界条件: y(0)=0, y'(0)=0 は自動的に満たされる（累積和の初期値が0）
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
# BMD
axes[0].plot(x_arr, M_arr / 1000, 'b-', linewidth=2)
axes[0].fill_between(x_arr, M_arr / 1000, alpha=0.2)
axes[0].axhline(y=0, color='k', linewidth=0.5)
axes[0].set_xlabel('Position x [m]', fontsize=12)
axes[0].set_ylabel('Bending Moment M [kN·m]', fontsize=12)
axes[0].set_title('Bending Moment Diagram', fontsize=13)
axes[0].grid(True, alpha=0.3)
# たわみ曲線
axes[1].plot(x_arr, deflection * 1000, 'r-', linewidth=2)
axes[1].axhline(y=0, color='k', linewidth=0.5)
axes[1].set_xlabel('Position x [m]', fontsize=12)
axes[1].set_ylabel('Deflection y [mm]', fontsize=12)
axes[1].set_title('Deflection Curve (Statically Indeterminate Beam)', fontsize=13)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"\n最大たわみ: {np.min(deflection)*1000:.4f} mm")
print(f"最大たわみ位置: x = {x_arr[np.argmin(deflection)]:.3f} m")
このコードの出力から、いくつかの重要な点が読み取れます:
- 余剰反力: $R_B = 3wL/8$ が得られます。単純支持はり（$R = wL/2$）に比べて支持端の反力が小さくなっているのは、固定端がモーメントを分担しているためです。
- BMD（曲げモーメント図）: 固定端で負のモーメント（上凸の曲げ）が生じ、スパン途中で正のモーメントに変わります。BMDがゼロを横切る点（反曲点）が存在するのが不静定はりの特徴です。
- たわみ曲線: 固定端で変位・傾きともゼロ、単純支持端で変位ゼロ（傾きは非ゼロ）という境界条件が満たされています。最大たわみは中央よりやや固定端寄りに生じます。
カスティリアーノの定理を使った不静定問題の解法では、余剰反力に関する $\partial U / \partial R = 0$ という1本の方程式を解くだけで全ての反力が決定できます。次数の高い不静定問題でも、余剰反力の数だけ連立方程式を立てれば解くことができ、体系的なアプローチが可能です。
Pythonでの実装6: 各種はりのたわみ比較
最後に、カスティリアーノの定理で求めたたわみの公式を使って、代表的なはりの支持条件ごとのたわみを比較します。
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ ---
E = 200e9 # ヤング率 [Pa]
I = 8.33e-6 # 断面二次モーメント [m^4]
EI = E * I
# はりの長さと荷重
L_values = np.linspace(0.5, 5.0, 100)
P = 10000 # 集中荷重 [N]
w = 5000 # 分布荷重 [N/m]
# --- たわみの公式（カスティリアーノの定理から導出した結果）---
# 1. 片持ちはり・先端集中荷重: δ = PL^3/(3EI)
delta_cantilever_P = P * L_values**3 / (3 * EI)
# 2. 片持ちはり・等分布荷重: δ = wL^4/(8EI)
delta_cantilever_w = w * L_values**4 / (8 * EI)
# 3. 単純支持はり・中央集中荷重: δ = PL^3/(48EI)
delta_simple_P = P * L_values**3 / (48 * EI)
# 4. 単純支持はり・等分布荷重: δ = 5wL^4/(384EI)
delta_simple_w = 5 * w * L_values**4 / (384 * EI)
# 5. 両端固定はり・中央集中荷重: δ = PL^3/(192EI)
delta_fixed_P = P * L_values**3 / (192 * EI)
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左: 集中荷重
axes[0].plot(L_values, delta_cantilever_P * 1000, 'r-', linewidth=2,
label='Cantilever')
axes[0].plot(L_values, delta_simple_P * 1000, 'b-', linewidth=2,
label='Simply supported')
axes[0].plot(L_values, delta_fixed_P * 1000, 'g-', linewidth=2,
label='Fixed-fixed')
axes[0].set_xlabel('Beam Length L [m]', fontsize=12)
axes[0].set_ylabel('Max Deflection [mm]', fontsize=12)
axes[0].set_title(f'Point Load P = {P/1000:.0f} kN', fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(0, None)
# 右: 分布荷重
axes[1].plot(L_values, delta_cantilever_w * 1000, 'r-', linewidth=2,
label='Cantilever')
axes[1].plot(L_values, delta_simple_w * 1000, 'b-', linewidth=2,
label='Simply supported')
axes[1].set_xlabel('Beam Length L [m]', fontsize=12)
axes[1].set_ylabel('Max Deflection [mm]', fontsize=12)
axes[1].set_title(f'Distributed Load w = {w/1000:.0f} kN/m', fontsize=13)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[1].set_ylim(0, None)
plt.suptitle("Beam Deflection Comparison (Castigliano's Theorem)",
fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()
# --- 数値例 ---
L_ex = 3.0
print(f"\n--- L = {L_ex} m での最大たわみ ---")
print(f"片持ち・集中荷重: {P*L_ex**3/(3*EI)*1000:.3f} mm")
print(f"片持ち・等分布荷重: {w*L_ex**4/(8*EI)*1000:.3f} mm")
print(f"単純支持・集中荷重: {P*L_ex**3/(48*EI)*1000:.3f} mm")
print(f"単純支持・等分布荷重: {5*w*L_ex**4/(384*EI)*1000:.3f} mm")
print(f"両端固定・集中荷重: {P*L_ex**3/(192*EI)*1000:.3f} mm")
print(f"\nたわみ比率（片持ち集中荷重 = 1）:")
ref = P * L_ex**3 / (3 * EI)
print(f"片持ち・集中荷重: 1.000")
print(f"単純支持・集中荷重: {P*L_ex**3/(48*EI)/ref:.3f}")
print(f"両端固定・集中荷重: {P*L_ex**3/(192*EI)/ref:.3f}")
このグラフから、支持条件によるたわみの違いが明瞭に読み取れます:
- 片持ちはりが最も大きいたわみを示します。これは一端が完全に自由であるため、構造全体としての剛性が低いことを反映しています。
- 単純支持はりは片持ちの1/16（集中荷重の場合）のたわみとなります。両端で支持されているため、同じ荷重に対する変形が大幅に抑えられます。
- 両端固定はりはさらに1/4（単純支持比）のたわみとなり、固定端のモーメント拘束がいかに構造を硬くするかがわかります。
- たわみはスパンの3乗（集中荷重）または4乗（分布荷重）に比例するため、スパンが長くなるとたわみは急激に増大します。これが長大スパンの構造設計において剛性の確保が重要である理由です。
これらの公式はすべてカスティリアーノの定理から統一的に導出できるものであり、エネルギー法の汎用性を示しています。
まとめ
本記事では、カスティリアーノの定理とエネルギー法について、ひずみエネルギーの基礎から応用まで解説しました。
- ひずみエネルギーは、弾性体に蓄えられるエネルギーであり、軸力・曲げ・せん断・ねじりの各荷重モードで計算式が異なる
- カスティリアーノの第2定理 $\delta_i = \partial U / \partial P_i$ により、ひずみエネルギーを外力で偏微分するだけで変位が求まる
- 仮想荷重法（ダミーロード法） を使えば、荷重が作用していない点の変位も求められる
- トラスの変位計算では、各部材の軸力のひずみエネルギーの総和から節点変位を求められる
- マクスウェル・ベッティの相反定理 $f_{ij} = f_{ji}$ は柔軟性行列の対称性から導かれ、構造解析の計算量削減に役立つ
- 不静定問題にも $\partial U / \partial R = 0$ の条件から余剰反力を求められる
- エネルギー法は有限要素法（FEM） の理論的基盤でもあり、構造力学・材料力学の中核をなす
カスティリアーノの定理は「たわみ曲線全体を求めず、知りたい点の変位だけを効率的に計算する」という実用性と、「エネルギーと力・変位の関係を明示する」という理論的な美しさを兼ね備えた定理です。
次のステップとして、以下の記事も参考にしてください。
- 梁のたわみ（弾性曲線方程式） — たわみ曲線全体を求める方法との比較
- 不静定問題の解法 — エネルギー法を含む不静定構造の解法を体系的に解説
- 応力とひずみ — ひずみエネルギーの出発点である応力-ひずみ関係の復習