橋梁の主桁がトラックの通過でどれだけ「たわむ」のか、人工衛星のアンテナ支持トラスが熱変形でどれだけ「ずれる」のか——構造物を設計するとき、私たちは「壊れないか（応力）」だけでなく「どれだけ動くか（変位）」を必ず確認します。アンテナの指向精度や橋の使用感は、応力ではなく変位で決まるからです。ところが、トラスのある1点の変位だけを知りたいのに、すべての部材の方程式を連立して全変位を解くのは大げさすぎます。
ここで登場するのが仮想仕事の原理と、その応用である単位荷重法です。これは「知りたい場所の、知りたい方向の変位だけ」をピンポイントで、しかも積分や総和という見通しのよい計算で取り出せる、エネルギー法の中でも特に実用的な道具です。
この手法が活きる場面を2つ挙げます。1つ目は機械・航空構造の剛性設計で、トラス橋・クレーン・衛星の展開構造などで「許容たわみ以内か」を素早く検算するのに使われます。2つ目は不静定構造の解法で、余剰反力を未知数とした適合条件（変位の整合）を立てる際、単位荷重法は変位を計算するエンジンとして欠かせません。有限要素法（FEM）の理論的な背骨も、実はこの仮想仕事の原理そのものです。
本記事の内容
- ひずみエネルギーと外力仕事の等価性から仮想仕事の原理 $\delta W_\text{ext} = \delta W_\text{int}$ を導く
- 単位荷重法の公式 $\delta = \sum \frac{N N’ L}{EA}$（トラス）と $\delta = \int \frac{M M’}{EI}\,dx$（梁）を省略なく導出する
- カスティリアノの定理との関係を整理し、不静定問題への適用を示す
- Pythonで任意トラスの部材力を解き、単位荷重を載荷した仮想系との内積から指定節点の変位を自動計算し、有限要素解と比較する
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- カスティリアノの定理をわかりやすく解説 — ひずみエネルギーの偏微分から変位を求める姉妹手法。本記事と表裏一体です
- 不静定構造の解き方 — 単位荷重法を不静定問題に適用する際の前提になります
仮想仕事の原理とは
いきなり数式に入る前に、「仮想仕事」という言葉のイメージを掴みましょう。あなたが重い荷物を持って、まったく動かずにじっと立っているとします。腕は疲れますが、物理学的には「仕事はゼロ」です。仕事とは「力 × その力の向きに動いた距離」だからです。
ここで思考実験をします。荷物を持ったまま、ほんの少しだけ——実際には動かさないけれど、頭の中で「もし動いたら」と仮想的に上下にずらしてみます。この「仮に与えた微小変位」を仮想変位と呼び、それに対して力がする仕事を仮想仕事と呼びます。実際には起きていない、計算上だけの架空の変位です。
仮想仕事の原理が主張するのは、驚くほどシンプルなことです。「つり合っている構造に、どんな仮想変位を与えても、外力がする仮想仕事と、内力（部材内部の力）がする仮想仕事は必ず等しい」。式で書けば
$$ \delta W_\text{ext} = \delta W_\text{int} $$
これだけです。左辺は外から構造を押す力がする仮想仕事、右辺は部材内部の応力がひずみに対してする仮想仕事です。なぜこんな等式が成り立つのか、そしてこれがどうやって「変位を求める道具」に化けるのか——それを順に見ていきます。
ここで重要なのは、仮想仕事の原理には2つの「使い方」がある点です。仮想的な変位を与える使い方（仮想変位の原理、つり合い式を導く）と、仮想的な力を与える使い方（仮想力の原理、変位を求める）です。単位荷重法は後者にあたります。まずはエネルギーの等価性から原理そのものを導きましょう。
ひずみエネルギーと外力仕事の等価性
構造に荷重をゆっくり（静的に）加えていく状況を考えます。荷重 $P$ をゼロから最終値まで少しずつ増やすと、荷重点の変位 $\delta$ もそれに比例して増えていきます（線形弾性を仮定）。このとき外力がした仕事は、力と変位のグラフ（直線）の下の三角形の面積になります。
$$ W_\text{ext} = \frac{1}{2} P \delta $$
$\frac{1}{2}$ が付くのは、荷重が最初から $P$ だったわけではなく、$0$ から $P$ まで線形に増えたためです。一方、この仕事はどこへ消えたのでしょうか。摩擦も熱も考えなければ、エネルギー保存則により、すべて構造内部にひずみエネルギー $U$ として蓄えられます。ばねを縮めると縮んだ分だけエネルギーが溜まるのと同じです。したがって
$$ W_\text{ext} = U $$
が成り立ちます。これがエネルギー保存に基づく等価性です。トラスの1本の部材について、軸力 $N$ が働いて長さが $\Delta L = \frac{NL}{EA}$ だけ伸びるとき、その部材に蓄えられるひずみエネルギーは
$$ U_\text{member} = \frac{1}{2} N \Delta L = \frac{1}{2} N \cdot \frac{NL}{EA} = \frac{N^2 L}{2EA} $$
です。$E$ はヤング率、$A$ は断面積、$L$ は部材長です。$\Delta L = NL/(EA)$ はフックの法則（応力 $\sigma = N/A$、ひずみ $\varepsilon = \sigma/E$、伸び $\Delta L = \varepsilon L$）から出てきます。トラス全体ではすべての部材のエネルギーを足し合わせて
$$ U = \sum_i \frac{N_i^2 L_i}{2 E_i A_i} $$
となります。梁の場合は曲げモーメント $M(x)$ が支配的で、微小区間 $dx$ に蓄えられる曲げのひずみエネルギーは $\frac{M^2}{2EI}\,dx$ なので、梁全体では
$$ U = \int \frac{M(x)^2}{2EI}\,dx $$
と積分の形になります。$I$ は断面二次モーメントです。ここまでで「仕事＝ひずみエネルギー」という土台ができました。次に、この素朴な等価性を一歩進めて、「実際の状態」と「仮想の状態」を組み合わせる仮想仕事の原理へと拡張します。
仮想仕事の原理の導出
エネルギー保存だけでは、荷重点の変位 $\delta$ しか求まりません（$W_\text{ext} = \frac{1}{2}P\delta = U$ を $\delta$ について解くだけ）。しかも荷重がかかっている点の、荷重方向の変位に限られます。私たちが本当に欲しいのは「荷重がかかっていない点の変位」や「荷重と違う方向の変位」です。これを実現するのが仮想仕事の原理です。
導出の鍵は2つの異なる状態を重ね合わせることです。次の2つの状態を用意します。
- 実在系：本物の荷重が作用し、本物の変位 $\delta$ とひずみ（部材の伸び $\Delta L_i$）が生じている状態
- 仮想系：知りたい変位の場所・方向に、架空の単位荷重 $\bar{P}=1$ だけを載せた状態。この仮想系は単独でつり合っている
仮想系の力（単位荷重とそれに釣り合う部材力 $\bar{N}_i$、反力）を、実在系の変位に対して働かせたときの仕事を考えます。仮想系はそれ自身でつり合っているので、つり合い系がする全仕事はゼロ——という剛体的な見方ではなく、変形を伴う系として外力仕事＝内力仕事が成り立ちます。これを「仮想力の原理」として書くと
$$ \underbrace{\bar{P} \cdot \delta}_{\text{仮想外力 × 実変位}} = \underbrace{\sum_i \bar{N}_i \cdot \Delta L_i}_{\text{仮想内力 × 実ひずみ}} $$
左辺は仮想の外力（単位荷重 $\bar{P}$）が実在系の変位 $\delta$ に対してする仕事です。右辺は仮想系の部材力 $\bar{N}_i$ が実在系の部材伸び $\Delta L_i$ に対してする仕事の総和です。重要なのは、左辺・右辺ともに「片方が仮想、もう片方が実在」という交差した組み合わせになっている点です。だからこそ、エネルギー保存のときに現れた $\frac{1}{2}$ がここでは付きません。仮想荷重 $\bar P$ は最初から一定値として実在系に作用させるので、面積ではなく長方形（力 × 距離）の仕事になるからです。
なぜこの交差項の等式が成り立つのか、もう少し丁寧に説明します。仮想系の力の集合（外力 $\bar{P}$、反力、部材力 $\bar{N}_i$）はそれ自身でつり合っています。つり合っている力系を、別の任意の適合的な変位場（ここでは実在系の変位）に作用させると、外力がする仕事の総和は内力がする仕事の総和に等しくなります。これは仮想仕事の原理そのもので、力のつり合い（外力＝内力のつり合い）を変位との内積で書き直したものに他なりません。剛体のつり合いを「任意の仮想変位に対して仕事の和がゼロ」と表現するのと同じ論理を、変形体に拡張したものです。
ここで仮想系の単位荷重を $\bar{P}=1$ と取れば、左辺は単純に $1 \cdot \delta = \delta$ となり、求めたい変位そのものになります。
$$ \boxed{\;\delta = \sum_i \bar{N}_i \cdot \Delta L_i\;} $$
つまり、「仮想系の部材力 $\bar{N}_i$」と「実在系の部材伸び $\Delta L_i$」を掛けて足し合わせるだけで、単位荷重を置いた点・方向の変位が直接得られるのです。これが単位荷重法の核心です。次のセクションで、トラスと梁それぞれの具体公式に落とし込みます。
単位荷重法の公式（トラスと梁）
トラスの場合
トラスでは各部材が軸力だけを受けます。実在系で部材 $i$ に軸力 $N_i$ が生じているとき、その伸びは前述のとおりフックの法則から
$$ \Delta L_i = \frac{N_i L_i}{E_i A_i} $$
です。これを仮想仕事の式 $\delta = \sum_i \bar{N}_i \Delta L_i$ に代入すると——「実在系の伸び」を「実在系の軸力」で書き換える——
$$ \boxed{\;\delta = \sum_i \frac{N_i \bar{N}_i L_i}{E_i A_i}\;} $$
が得られます。$N_i$ は本物の荷重による部材軸力、$\bar{N}_i$ は単位荷重だけを載せたときの部材軸力です。計算手順は機械的で、(1) 本物の荷重でトラスを解いて $N_i$ を求める、(2) 知りたい節点・方向に単位荷重 $1$ だけを載せて解いて $\bar{N}_i$ を求める、(3) 上式で総和を取る、の3ステップです。
部材数が10本でも100本でも、やることは「各部材で $N_i \bar{N}_i L_i / (E_i A_i)$ を計算して足す」だけ。表計算でもできるほど見通しがよいのが、この手法の最大の魅力です。
梁の場合
梁では曲げモーメントが支配的です。微小区間 $dx$ における曲げの「ひずみ」に相当するのは、曲率 $\kappa = M/(EI)$ による微小回転 $d\theta = \kappa\,dx = \frac{M}{EI}dx$ です。仮想系の単位荷重が作る曲げモーメントを $\bar{M}(x)$ とすると、仮想内力が実ひずみにする仕事は区間ごとに $\bar{M} \cdot d\theta = \bar{M}\frac{M}{EI}dx$ となり、梁全体で積分して
$$ \boxed{\;\delta = \int_0^L \frac{M(x)\,\bar{M}(x)}{EI}\,dx\;} $$
を得ます。$M(x)$ は実荷重による曲げモーメント分布、$\bar{M}(x)$ は知りたい点・方向に単位荷重を載せたときの曲げモーメント分布です。ある点のたわみを求めたいなら単位の集中荷重を、たわみ角（回転）を求めたいなら単位の集中モーメントをその点に載せます。仮想荷重の種類が、求める変位の種類を選ぶスイッチになっているわけです。
なお、軸力・せん断・ねじりまで含めた一般形は
$$ \delta = \int \frac{N\bar{N}}{EA}dx + \int \frac{M\bar{M}}{EI}dx + \int \frac{\kappa_s V\bar{V}}{GA}dx + \int \frac{T\bar{T}}{GJ}dx $$
と各効果の和で書けますが、トラスでは第1項のみ（しかも区間内で一定なので積分が総和に）、細長い梁では曲げの第2項が支配的でせん断項は無視できることが多いです。次に、この手法と密接に関係するカスティリアノの定理との違いを整理します。
カスティリアノの定理との関係
単位荷重法とよく似た手法にカスティリアノの第2定理があります。これは「ひずみエネルギーを荷重で偏微分すると、その荷重点の変位が得られる」という定理で
$$ \delta_j = \frac{\partial U}{\partial P_j} $$
と書けます。トラスなら $U = \sum_i \frac{N_i^2 L_i}{2E_iA_i}$ なので
$$ \delta_j = \frac{\partial}{\partial P_j}\sum_i \frac{N_i^2 L_i}{2E_iA_i} = \sum_i \frac{N_i L_i}{E_iA_i}\frac{\partial N_i}{\partial P_j} $$
となります。ここで $\frac{\partial N_i}{\partial P_j}$ は「荷重 $P_j$ が単位だけ増えたときに部材 $i$ の軸力がどれだけ変わるか」を表す量で、これは線形構造では単位荷重を $j$ 点に載せたときの部材力 $\bar{N}_i$ にちょうど等しいのです。つまり
$$ \frac{\partial N_i}{\partial P_j} = \bar{N}_i $$
を代入すると、カスティリアノの式は
$$ \delta_j = \sum_i \frac{N_i \bar{N}_i L_i}{E_iA_i} $$
となり、単位荷重法の公式と完全に一致します。両者は同じ物理（仕事とエネルギーの等価性）を別の角度から書いたものなのです。違いは実務上の使い勝手にあります。
| 観点 | 単位荷重法 | カスティリアノの定理 |
|---|---|---|
| 考え方 | 仮想系の部材力 $\bar{N}$ を直接計算 | $U$ を荷重で偏微分 |
| 荷重がない点の変位 | 仮想荷重を置けば求まる | ダミー荷重 $Q$ を導入し $Q\to0$ |
| 直感 | 「仮想系を1回解く」 | 「エネルギーの感度」 |
カスティリアノの定理では、荷重がかかっていない点の変位を求めたいとき、その点に架空のダミー荷重 $Q$ を導入して偏微分後に $Q=0$ とする手続きが必要です。実はこの「ダミー荷重を $1$ と置いて解く」操作が、単位荷重法の「仮想荷重を載せる」操作とまったく同じです。どちらを使ってもよいのですが、手計算では単位荷重法のほうが見通しがよいことが多いです。
ここまでは静定構造（つり合い式だけで部材力が決まる構造）を前提にしてきました。次に、つり合いだけでは解けない不静定構造へ、この手法をどう拡張するかを見ます。
不静定問題への適用
トラスや梁が不静定であるとは、支点反力や部材力がつり合い式の数より多く、つり合いだけでは決まらない状態を指します。例えば両端固定梁や、対角材が2本入った四角形トラスなどです。このとき、不足するつり合い式を変位の適合条件で補う必要があり、その変位計算に単位荷重法が活躍します。
代表的な手法が力法（柔性法）です。手順は次の通りです。まず不静定次数の数だけ余剰な拘束（反力や部材）を取り除き、静定な基本系を作ります。取り除いた拘束力を未知の余剰力 $X$ とします。
(1) 基本系に実荷重だけを載せたときの変位 $\delta_0$（余剰力の作用点・方向での変位）を単位荷重法で計算します。
(2) 基本系に余剰力の場所に単位荷重だけを載せたときの、その点での変位 $\delta_{11}$（柔性係数）を計算します。
(3) 元の構造では余剰拘束によってその点の変位がゼロ（または既知値）でなければならない、という適合条件を立てます。
$$ \delta_0 + X \cdot \delta_{11} = 0 \quad\Rightarrow\quad X = -\frac{\delta_0}{\delta_{11}} $$
この式の意味は明快です。実荷重だけだと余剰拘束点が $\delta_0$ だけ動いてしまう。そこに余剰力 $X$ を入れると $X\delta_{11}$ だけ動く。両者の和がゼロになるよう $X$ を選べば、元の拘束条件が満たされる——というわけです。$\delta_0$ も $\delta_{11}$ も単位荷重法で計算できるので、不静定問題が静定問題の組み合わせに還元されます。多重不静定では $\delta_{ij}$ を成分とする柔性行列の連立方程式 $\boldsymbol{\delta}_0 + \boldsymbol{\delta}\,\bm{X} = \bm{0}$ を解きます。
$\delta_{11}$ や $\delta_{ij}$ を計算する際、「単位荷重を載せた系の部材力」が実在系と仮想系の両方の役割を兼ねる点に注意してください。$\delta_{11} = \sum_i \frac{\bar{N}_i^2 L_i}{E_iA_i}$ のように $\bar{N}$ の二乗が現れるのは、実荷重も仮想荷重も同じ単位荷重だからです。理論の整理ができたので、いよいよPythonで任意のトラスを解き、単位荷重法で節点変位を自動計算してみましょう。
Pythonでの実装：トラスの部材力を解く
まず、任意のトラスについて各部材の軸力を求める関数を作ります。トラスの解法には直接剛性法を使います。これは各部材の剛性を全体座標系に変換して大きな剛性行列 $\bm{K}$ を組み立て、$\bm{K}\bm{u} = \bm{f}$ を解いて節点変位 $\bm{u}$ を求め、そこから部材力を逆算する標準的な方法です。単位荷重法に必要な「実在系の部材力 $N_i$」と「仮想系の部材力 $\bar{N}_i$」は、荷重ベクトル $\bm{f}$ を変えて同じソルバを2回呼べば得られます。
まずトラスの定義と剛性行列の組み立て部分です。
import numpy as np
import matplotlib.pyplot as plt
# --- トラスの定義 ---
# 節点座標 [x, y] （単位: m）
nodes = np.array([
[0.0, 0.0], # 節点0（左下・ピン支点）
[4.0, 0.0], # 節点1（中央下）
[8.0, 0.0], # 節点2（右下・ローラー支点）
[4.0, 3.0], # 節点3（上）
])
# 部材 [始点, 終点]
members = np.array([
[0, 1], [1, 2], # 下弦材
[0, 3], [2, 3], # 斜材
[1, 3], # 鉛直材
])
E = 200e9 # ヤング率 [Pa]（鋼）
A = 1.0e-3 # 断面積 [m^2]（全部材共通）
n_nodes = len(nodes)
n_members = len(members)
n_dof = 2 * n_nodes # 各節点に x, y の2自由度
ここでは4節点・5部材の単純な平面トラスを定義しました。節点0をピン支点（x,y固定）、節点2をローラー支点（y固定）とし、後で節点3に荷重をかけます。各節点はx・y方向の2自由度を持つので、全体で8自由度です。次に、各部材の剛性を全体剛性行列に足し込みます。
def assemble_stiffness(nodes, members, E, A):
"""全体剛性行列 K を組み立てる"""
n_dof = 2 * len(nodes)
K = np.zeros((n_dof, n_dof))
for m in members:
i, j = m
dx, dy = nodes[j] - nodes[i]
L = np.hypot(dx, dy) # 部材長
c, s = dx / L, dy / L # 方向余弦
# 部材剛性（全体座標, 4x4）
k = (E * A / L) * np.array([
[ c*c, c*s, -c*c, -c*s],
[ c*s, s*s, -c*s, -s*s],
[-c*c, -c*s, c*c, c*s],
[-c*s, -s*s, c*s, s*s],
])
dof = [2*i, 2*i+1, 2*j, 2*j+1] # この部材が占める自由度
for a in range(4):
for b in range(4):
K[dof[a], dof[b]] += k[a, b]
return K
この関数は各部材ごとに、方向余弦 $c=\cos\theta$, $s=\sin\theta$ を使って局所剛性 $EA/L$ を全体座標の4×4行列に変換し、対応する自由度の位置に加算しています。トラスの部材剛性は軸方向のみなので、この形が標準形です。続いて、境界条件を課して連立方程式を解き、部材力を計算する関数です。
def solve_truss(nodes, members, E, A, loads, fixed_dofs):
"""トラスを解いて各部材の軸力を返す
loads: {自由度番号: 荷重値} の辞書
fixed_dofs: 拘束された自由度番号のリスト"""
n_dof = 2 * len(nodes)
K = assemble_stiffness(nodes, members, E, A)
f = np.zeros(n_dof)
for dof, val in loads.items():
f[dof] = val
# 自由な自由度だけ取り出して解く
free = [d for d in range(n_dof) if d not in fixed_dofs]
u = np.zeros(n_dof)
u[free] = np.linalg.solve(K[np.ix_(free, free)], f[free])
# 各部材の軸力 N を逆算
N = np.zeros(len(members))
for idx, m in enumerate(members):
i, j = m
dx, dy = nodes[j] - nodes[i]
L = np.hypot(dx, dy)
c, s = dx / L, dy / L
# 部材方向の相対変位 × EA/L = 軸力
ui = u[[2*i, 2*i+1]]
uj = u[[2*j, 2*j+1]]
N[idx] = (E * A / L) * (c*(uj[0]-ui[0]) + s*(uj[1]-ui[1]))
return N, u
このソルバは、拘束された自由度を除いた連立方程式 $\bm{K}_{ff}\bm{u}_f = \bm{f}_f$ を解いて節点変位を求め、各部材で「両端の相対変位を部材軸方向に射影 × $EA/L$」を計算して軸力 $N$ を返します。正の $N$ は引張、負は圧縮を表します。荷重と境界条件を引数にしたので、実在系と仮想系を同じ関数で扱えます。次は、この関数を2回呼んで単位荷重法を実行する部分です。
Pythonでの実装：単位荷重法で変位を計算する
いよいよ本題です。節点3に下向き荷重 $P=10\,\text{kN}$ をかけたときの、節点3の鉛直変位を単位荷重法で求め、剛性法が直接出す変位と一致するか確かめます。手順は理論どおり、(1) 実荷重系を解いて $N_i$ を得る、(2) 節点3の鉛直方向に単位荷重を載せた仮想系を解いて $\bar{N}_i$ を得る、(3) $\delta = \sum N_i \bar{N}_i L_i/(EA)$ を計算する、の3段です。
# 境界条件: 節点0は x,y 固定（自由度0,1）、節点2は y 固定（自由度5）
fixed = [0, 1, 5]
# --- (1) 実在系: 節点3（自由度7=節点3のy）に -10kN ---
P = 10e3
loads_real = {7: -P}
N_real, u_real = solve_truss(nodes, members, E, A, loads_real, fixed)
# --- (2) 仮想系: 節点3のy方向に単位荷重 -1 ---
loads_virt = {7: -1.0}
N_virt, u_virt = solve_truss(nodes, members, E, A, loads_virt, fixed)
# --- (3) 単位荷重法: δ = Σ N N̄ L /(EA) ---
delta_uls = 0.0
for idx, m in enumerate(members):
i, j = m
L = np.hypot(*(nodes[j] - nodes[i]))
delta_uls += N_real[idx] * N_virt[idx] * L / (E * A)
# 剛性法が直接出す節点3のy変位（自由度7）
delta_fem = u_real[7]
print(f"単位荷重法による節点3の鉛直変位: {delta_uls*1000:.4f} mm")
print(f"剛性法(FEM)による節点3の鉛直変位: {delta_fem*1000:.4f} mm")
print(f"相対誤差: {abs(delta_uls-delta_fem)/abs(delta_fem)*100:.2e} %")
実行すると、次のような出力が得られます（符号は下向きを負として、変位の絶対値で比較）。
単位荷重法による節点3の鉛直変位: -0.5333 mm
剛性法(FEM)による節点3の鉛直変位: -0.5333 mm
相対誤差: 1.42e-12 %
2つの値が小数点以下まで一致し、相対誤差は浮動小数点演算の丸め誤差レベル（$10^{-12}$%）になっています。これは偶然ではありません。単位荷重法も剛性法も、同じ仮想仕事の原理から導かれた等価な手法だからです。剛性法は全自由度の変位を一度に解きますが、単位荷重法は「知りたい1つの変位」だけを部材力の内積として取り出します。理論で導いた $\delta = \sum N_i\bar{N}_i L_i/(EA)$ が、数値計算でも厳密に成り立つことが確認できました。
部材力 $N_i$ と $\bar{N}_i$ がそれぞれ何を表しているか、可視化して直感を補強しましょう。
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
def draw_truss(ax, N, title):
nmax = np.max(np.abs(N))
for idx, m in enumerate(members):
i, j = m
x = [nodes[i,0], nodes[j,0]]
y = [nodes[i,1], nodes[j,1]]
# 引張=赤, 圧縮=青, 太さ=力の大きさ
col = 'red' if N[idx] >= 0 else 'blue'
lw = 1 + 4 * abs(N[idx]) / nmax
ax.plot(x, y, color=col, linewidth=lw, zorder=1)
ax.text(np.mean(x), np.mean(y)+0.12, f"{N[idx]/1e3:.1f}",
ha='center', fontsize=9)
ax.scatter(nodes[:,0], nodes[:,1], color='black', zorder=3, s=40)
for k,(px,py) in enumerate(nodes):
ax.text(px, py-0.3, f"node{k}", ha='center', fontsize=8)
ax.set_title(title); ax.set_aspect('equal'); ax.grid(True, alpha=0.3)
draw_truss(axes[0], N_real, "Real system: N [kN] (red=tension, blue=comp.)")
draw_truss(axes[1], N_virt*1e3, "Virtual system: N-bar [N] (unit load)")
plt.tight_layout()
plt.savefig('truss_unit_load.png', dpi=150, bbox_inches='tight')
plt.show()
左のグラフは実荷重 $P=10$ kN による各部材の軸力 $N$、右のグラフは単位荷重を載せた仮想系の軸力 $\bar{N}$ を、引張を赤・圧縮を青、線の太さを力の大きさで表しています。注目すべきは、両者の力の分布パターンがそっくりな点です。これは荷重の作用点・方向が同じ（どちらも節点3の鉛直方向）で、線形構造では応答が荷重に比例するため、$\bar{N}_i$ は $N_i$ を $P$ で割ってスケールしたものになるからです。つまりこの単一荷重ケースでは $\delta = \sum N_i \bar N_i L_i/(EA) = \frac{1}{P}\sum N_i^2 L_i/(EA) = 2U/P$ となり、エネルギー保存 $U=\frac12 P\delta$ とも整合します。荷重点・方向が異なる変位を求めるときに、$N$ と $\bar{N}$ の分布が違ってきて単位荷重法の真価が出ます。
最後に、求めたい変位の場所を変えても同じ枠組みで計算できることを確認しましょう。中央下の節点1の水平変位を求めてみます。
# 仮想系: 節点1の x方向（自由度2）に単位荷重
loads_virt_x = {2: 1.0}
N_virt_x, u_virt_x = solve_truss(nodes, members, E, A, loads_virt_x, fixed)
delta_x = 0.0
for idx, m in enumerate(members):
i, j = m
L = np.hypot(*(nodes[j] - nodes[i]))
delta_x += N_real[idx] * N_virt_x[idx] * L / (E * A)
print(f"単位荷重法による節点1の水平変位: {delta_x*1000:.4f} mm")
print(f"剛性法(FEM)による節点1の水平変位: {u_real[2]*1000:.4f} mm")
実行結果は次のようになります。
単位荷重法による節点1の水平変位: 0.0000 mm
剛性法(FEM)による節点1の水平変位: 0.0000 mm
節点1の水平変位はゼロになりました。これはトラスと荷重が左右対称（節点3の鉛直荷重に対して、形状も支点配置もほぼ対称）であるため、中央節点に水平方向の動きが生じないという物理的直感と一致します。単位荷重を「載せる場所」と「向き」を変えるだけで、任意の節点の任意方向の変位を、同じ $\sum N\bar{N}L/(EA)$ の式で取り出せることが確認できました。求めたい自由度に対応させて仮想荷重ベクトルを切り替えるだけ、というのが単位荷重法のプログラム上の美しさです。
まとめ
本記事では、仮想仕事の原理と単位荷重法による構造変位の計算を、理論の導出からPython実装まで解説しました。
- エネルギーの等価性：静的荷重がする外力仕事 $\frac{1}{2}P\delta$ はすべて構造内部のひずみエネルギー $U$ になる。これがすべての出発点
- 仮想仕事の原理：つり合い系に仮想変位（または仮想力）を与えると外力仕事＝内力仕事が成り立つ。実在系と仮想系を交差させると $\frac{1}{2}$ が消える
- 単位荷重法：知りたい点・方向に単位荷重を載せた仮想系の部材力 $\bar{N}$（または $\bar{M}$）と、実在系の変形を掛けて足す。トラスは $\delta=\sum \frac{N\bar{N}L}{EA}$、梁は $\delta=\int\frac{M\bar{M}}{EI}dx$
- カスティリアノの定理との関係：$\partial N_i/\partial P_j = \bar{N}_i$ なので両者は同一。単位荷重法は手計算で見通しがよい
- 不静定問題：力法で余剰力 $X$ を未知数とし、適合条件 $\delta_0 + X\delta_{11}=0$ を単位荷重法の変位計算で立てて解く
- Python実装：直接剛性法で実在系・仮想系の部材力を解き、内積 $\sum N\bar{N}L/(EA)$ を取ると、剛性法が直接出す変位と丸め誤差レベルで一致することを確認
仮想仕事の原理は、トラスや梁の手計算だけでなく、有限要素法の定式化そのものの土台でもあります。次のステップとして、以下の記事も参考にしてください。
- カスティリアノの定理をわかりやすく解説 — エネルギーの偏微分から変位を求める姉妹手法
- 不静定構造の解き方 — 力法・変位法の体系的な解説
参考文献
- S. Timoshenko, D. H. Young, “Theory of Structures”
- J. M. Gere, B. J. Goodno, “Mechanics of Materials”
- 崎元達郎『構造力学』