単位荷重法とは？仮想仕事によるたわみ計算の導出と実装

橋の中央がどれくらい下がるか、飛行機の主翼の先端がどれくらいしなるか、衛星に載せた望遠鏡の取付架台が打上げ荷重でどれだけ傾くか。構造設計の現場では「応力が許容値以内か」と同じくらい「特定の1点が、特定の方向にどれだけ動くか」が問われます。強度は足りていても、たわみが大きすぎれば橋は歩けたものではありませんし、架台が0.01度傾いただけで望遠鏡は狙った天体を外します。

ところが、この「1点だけの変位」を求めるのは意外と面倒です。素直にやるなら弾性曲線方程式 $EI\,v” = M(x)$ を4回積分し、境界条件を4本立てて連立させ、たわみ曲線 $v(x)$ を全区間で求めてから、最後に知りたい点の座標を代入する——知りたいのは1個の数字なのに、途中で関数を1本まるごと求めている。荷重が複数あって区間が分かれていれば、積分定数は区間の数だけ倍増します。トラスに至っては、そもそも「たわみ曲線」という概念すらありません。

単位荷重法（unit load method）は、この回り道を一発で飛び越えます。やることはたった1つ、知りたい点に、知りたい方向へ、大きさ1の架空の荷重を置くだけ。すると、その1点の変位だけが積分の形でぽろりと取り出せます。式にすれば

$$ 1 \cdot \delta = \int \frac{M m}{EI}\,dx $$

これだけです。$M$ は本物の荷重による曲げモーメント、$m$ は「大きさ1の架空の荷重」による曲げモーメント。この2つを掛けて $EI$ で割って積分すれば、答えが出ます。積分定数も境界条件の連立も要りません。

左の実状態では、荷重が構造全体に載っていてたわみは無限個の点・無限個の方向に分布しており、そのうち点Cの1点だけを取り出す手立てがありません。右の仮想状態では本物の荷重をすべて外し、知りたい点・知りたい向きに大きさ1の荷重をただ1本だけ置きます。この2つを掛け合わせると、荷重が1本しかない仮想状態のおかげで外力仕事の項も1つしか生まれず、点Cの変位だけが残る——これが単位荷重法のすべてです。

この手法は、単に手計算のテクニックとして便利なだけではありません。応用先を2つ挙げておきます。

• 不静定構造の解析：釣り合い式だけでは解けない構造（両端固定はり、連続はり、ラーメン、多くの実橋梁）は、「余分な支点反力を未知数 $X$ に置き換え、その点の変位がゼロになるように $X$ を決める」という手順で解きます。この「変位がゼロ」を書き下す道具がまさに単位荷重法で、得られる式 $\delta_{10} + X\delta_{11} = 0$ を整合条件（compatibility condition）と呼びます。構造力学の応力法（力法）はこれ一本で成り立っています。
• 有限要素法の柔性行列と感度解析：単位荷重法から自然に出てくる $\delta_{ij} = \int m_i m_j / EI\,dx$ は、そのまま構造の柔性行列（flexibility matrix）の成分です。この行列が対称であること（$\delta_{ij}=\delta_{ji}$、Maxwellの相反定理）は、剛性行列が対称であることの裏返しで、FEMソルバの数値検証にも使えます。さらに「断面 $I$ を少し変えたら変位がどう変わるか」という設計感度も、この積分表現から直接微分するだけで得られます。

この記事では、単位荷重法を「公式の暗記」ではなく「なぜその形になるのか」から丁寧に組み立てます。仮想仕事の原理から出発して $\delta=\int Mm/EI\,dx$ を導き、軸力・せん断・ねじりまで含めた一般形を示し、細長いはりでなぜ曲げ項だけ残るのかをオーダー評価で確かめます。そのうえで、トラスへの拡張、相反定理、1次不静定問題までを一本の筋道で通し、最後にsympyとnumpyで全部を数値的に検証します。

本記事の内容

• 「知りたい点に単位荷重を置くと、その点の変位だけが取り出せる」という核心の直感
• 仮想仕事の原理からMaxwell–Mohrの式 $\delta=\int Mm/EI\,dx$ を省略なく導出
• 軸力・せん断・ねじりを含む一般形と、細長いはりで曲げ項が支配的になる理由
• トラスへの拡張 $\delta=\sum N_i n_i L_i/(E_iA_i)$ とMaxwellの相反定理
• 1次不静定はりを整合条件 $\delta_{10}+X\delta_{11}=0$ で解く
• sympy・numpyによる公式の記号的検証、FEMとの数値照合、BMDの再分配の可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解がスムーズです。

• 内力図（SFD・BMD）の描き方 — $M(x)$ と $m(x)$ を書き下せることが単位荷重法の前提です
• 梁のたわみ（弾性曲線方程式） — 単位荷重法の答え合わせに使います
• 仮想仕事の原理と単位荷重法による構造変位の理論と実装 — 本記事の理論的な土台
• カスティリアーノの定理とエネルギー法を基礎から徹底解説 — 単位荷重法と表裏一体の関係にあります

単位荷重法とは — 知りたい成分だけを取り出す「フィルタ」

いきなり式に入る前に、単位荷重法が何をしている道具なのかを掴んでおきましょう。

ベクトルの内積を思い出してください。3次元ベクトル $\bm{u}=(u_1,u_2,u_3)$ の第2成分だけが知りたいとき、私たちは $\bm{e}_2=(0,1,0)$ という単位ベクトルとの内積 $\bm{e}_2 \cdot \bm{u} = u_2$ を取ります。$\bm{u}$ の中身を全部調べなくても、「知りたい成分に1を立てたベクトル」を掛けるだけで、その成分だけがすっと抜き出せる。他の成分には0が掛かるので、勝手に消えてくれます。

単位荷重法がやっていることは、これとまったく同じです。構造物の変位は、無限に多くの点・無限に多くの方向に分布しています。そのうち「点Cの鉛直方向」だけを知りたい。そこで点Cに鉛直方向の大きさ1の荷重を置くと、仕事という掛け算を通じて、点Cの鉛直変位だけに1が掛かり、それ以外の変位には（その架空の荷重が仕事をしないので）0が掛かって消えます。結果として、外力側の仕事は $1 \times \delta_C$ というただ1項になります。

左のグラフでは、ベクトル $\bm{u}$ の3成分に単位ベクトル $\bm{e}_2$ を掛けると、第1・第3成分には0が掛かって消え、第2成分だけがそのまま残っています。右の構造でもまったく同じことが起きていて、単位荷重を置いていない点の変位には「0」という係数が掛かるので、仕事の総和には現れません。両者の違いは「成分が有限個か、無限個の連続体か」だけで、抜き出しの原理は同一です。

もう少し噛み砕きます。仕事とは「力 × その力の方向の変位」です。構造物にいくつ荷重が載っていようと、架空の荷重を1つしか置かなければ、外力仕事の項も1つしか出てきません。しかも大きさを1にしておけば、係数も1。だから外力仕事はそのまま「知りたい変位そのもの」になります。あとはエネルギー保存（正確には仮想仕事の原理）で「外部の仕事 = 内部の仕事」と等号を結べば、知りたい変位が内部の積分で表される、という寸法です。

用語も整理しておきます。単位荷重法では、2つの状態を同時に考えます。

この2つは別々の問題であり、混ぜてはいけません。仮想状態は変形を計算するための「測定用のプローブ」であって、実際に構造に載る荷重ではありません。逆に実状態の内力は「変形の情報源」であって、仮想状態には一切影響しません。この分離こそが単位荷重法の要で、初学者が最も混乱するポイントでもあります。以下では、記号を大文字（$M,N,V,T$）と小文字（$m,n,v,t$）で徹底的に区別します。

なぜこの2つを自由に組み合わせてよいのか、そしてなぜ「外部の仕事 = 内部の仕事」が成り立つのか。次節でその根拠となる仮想仕事の原理を確認します。

仮想仕事の原理 — 釣り合う力系と、適合する変形は独立に選べる

仮想仕事の原理には2つの顔があります。よく教科書の最初に出てくるのは仮想変位の原理で、「釣り合っている実際の力系に、架空の微小変位を与えると、外力の仕事と内力の仕事が釣り合う」というものです。単位荷重法で使うのはその双対にあたる仮想力の原理（相補仮想仕事の原理）で、力と変位の役割が入れ替わります。

言葉で述べると、こうです。

ある弾性体について、釣り合いを満たす任意の力系（外力 $\bar{\bm{P}}$ と内部応力 $\bar{\bm{\sigma}}$ の組）と、適合条件を満たす任意の変形場（変位 $\bm{u}$ とひずみ $\bm{\varepsilon}$ の組）を、互いに無関係に選んでよい。このとき次が恒等的に成り立つ。

$$ \begin{equation} \sum_i \bar{P}_i\, u_i = \int_V \bar{\bm{\sigma}} : \bm{\varepsilon}\, dV \end{equation} $$

左辺は「架空の外力が実際の変位に対してする仕事」、右辺は「架空の応力が実際のひずみに対してする仕事」です。ここで強調したいのは、$\bar{\bm{\sigma}}$ と $\bm{\varepsilon}$ の間に構成則が要らないことです。$\bar{\bm{\sigma}}$ は釣り合っていればよく、$\bm{\varepsilon}$ は適合していればよい。両者は別の問題の答えでも構いません。

この式が成り立つ理由は、実は純粋に数学的です。釣り合い式 $\nabla\cdot\bar{\bm{\sigma}} + \bar{\bm{b}} = \bm{0}$ に実変位 $\bm{u}$ を掛けて体積積分し、発散定理を使って部分積分すると、表面項が左辺（外力の仕事）に、体積項が右辺（内部の仕事）になります。物理法則というより、釣り合い式を変位で重み付けしただけの恒等式だと思ってください。だからこそ、材料が線形か非線形かにも依存しません。

左の箱には「釣り合いさえ満たしていればよい力系」、右の箱には「適合条件さえ満たしていればよい変形場」が入っており、この2つは別々の問題の答えでも構いません。両者を掛け合わせた下段の等式が仮想力の原理で、材料則を一度も経由していないことが図の構造から読み取れます。右の変形場の原因が荷重でも温度でも製作誤差でもよい、という後の一般化（温度項）は、この「原因を問わない」性質から直接出てきます。

この「独立に選んでよい」という性質が決定的です。私たちは仮想状態として「知りたい点に単位荷重を置いた、静定で簡単な釣り合い系」を勝手にでっち上げ、実状態として「本物の荷重による変形」を使う。この2つを式(1)に放り込むだけで、目的の変位が現れます。

ここまでで道具は揃いました。次は式(1)の右辺を、はりの理論の言葉（曲げモーメントと $EI$）に書き下していきます。

単位荷重法の導出 — Maxwell–Mohrの式

ゴールを先に宣言します。式(1)の右辺を、はりの断面力 $M$ と $m$ で表し、$\delta = \int Mm/(EI)\,dx$ を導くのが目標です。

実状態のひずみを書く

はりの曲げ理論（Euler–Bernoulli梁）では、断面内の軸方向応力は中立軸からの距離 $y$ に比例して分布します。実状態の曲げモーメントを $M(x)$、断面二次モーメントを $I$ とすると

$$ \sigma_{xx}(x,y) = \frac{M(x)\, y}{I} $$

材料がフックの法則に従うなら、対応する軸ひずみは応力をヤング率 $E$ で割って

$$ \varepsilon_{xx}(x,y) = \frac{\sigma_{xx}}{E} = \frac{M(x)\, y}{E I} $$

となります。これが式(1)の右辺に入る「実の変形」です。

仮想状態の応力を書く

仮想状態も同じ形をしています。単位荷重によって生じる曲げモーメントを $m(x)$ とすれば

$$ \bar{\sigma}_{xx}(x,y) = \frac{m(x)\, y}{I} $$

こちらは釣り合いさえ満たしていればよいので、材料則を経由する必要はありません。単位荷重を載せた静定はりのBMDを描けば、それで $m(x)$ は決まります。

体積積分を実行する

式(1)の右辺に代入します。曲げだけを考えるので $\bar{\bm{\sigma}}:\bm{\varepsilon}$ は $\bar\sigma_{xx}\varepsilon_{xx}$ の1項だけです。体積積分を「断面内の積分 $dA$」と「軸方向の積分 $dx$」に分けて書くと

$$ \int_V \bar{\sigma}_{xx}\varepsilon_{xx}\,dV = \int_0^L \!\!\int_A \frac{m(x)y}{I}\cdot\frac{M(x)y}{EI}\,dA\,dx $$

ここで、被積分関数のうち $x$ だけの関数（$m$, $M$, $I$, $E$）を断面積分の外に出すと、断面内に残るのは $y^2$ だけです:

$$ = \int_0^L \frac{m(x)M(x)}{EI^2}\left(\int_A y^2\,dA\right)dx $$

括弧の中は断面二次モーメントの定義そのもの $\int_A y^2 dA = I$ です。これを代入すると $I^2$ の片方が約分され

$$ = \int_0^L \frac{M(x)\,m(x)}{EI}\,dx $$

きれいに1次元の積分になりました。$y$ 方向の情報がすべて $I$ に吸収されるのが、この導出の気持ちのよいところです。

左の実状態も中央の仮想状態も、応力は中立軸からの距離 $y$ に比例した同じ形の分布をしています。この2つを掛けると符号によらず $y^2$ という正の重みだけが残り、右のグラフのようにその断面積分が $\int_A y^2 dA = I$ になります。$I^2$ の片方がここで約分されるので、断面という2次元の情報が完全に $EI$ の中に畳み込まれ、残るのは軸方向の1次元積分だけになるわけです。

左辺と結ぶ

一方、式(1)の左辺は「仮想の外力が実の変位に対してする仕事」でした。仮想状態の外力は、知りたい点Cに置いた大きさ1の荷重ただ1つです（支点反力も生じますが、支点では実変位がゼロなので仕事をしません）。したがって

$$ \sum_i \bar{P}_i u_i = 1\cdot \delta_C $$

の1項だけ。両辺を結んで、Maxwell–Mohrの式が得られます。

$$ \begin{equation} 1 \cdot \delta_C = \int_0^L \frac{M(x)\,m(x)}{EI}\,dx \end{equation} $$

符号のルール

式(2)を使うときのルールは2つだけです。

1. $M$ と $m$ は同じ符号規約で描く（例えば下面引張を正）。両方をひっくり返しても積 $Mm$ は変わらないので、規約そのものは何でも構いません。大事なのは「揃えること」です。
2. 答え $\delta_C$ の符号は、単位荷重の向きを正とする。$\delta_C>0$ なら単位荷重と同じ向きに変位し、$\delta_C<0$ なら逆向きです。

もう1つ実用上重要なのが、たわみ角を求めたいときは単位「モーメント」を置くことです。仕事は「力 × 変位」だけでなく「モーメント × 回転角」の形でも定義されるので、単位モーメント $1$ を置けば左辺は $1\cdot\theta_C$ となり、まったく同じ式でたわみ角が求まります。求めたい量に対応する（共役な）荷重を単位量だけ置く、というのが一般則です。

4つのパネルはどれも「求めたい量に共役な荷重を単位量だけ置く」という同じルールの現れです。たわみなら力、回転角ならモーメント、2点の相対量なら逆向きの対——置き方さえ変えれば、右辺の積分式 $\int Mm/EI\,dx$ は一切変わりません。左下の逆向きの単位力の対では、外力仕事が $1\cdot\delta_B – 1\cdot\delta_A$ となって相対変位がそのまま抜き出される点に注目してください。

ここまでで曲げによる変位は求まるようになりました。しかし実際の部材には軸力もせん断力もねじりも働きます。次はそれらを式(2)に足し込みます。

一般形 — 軸力・せん断・ねじりまで含める

式(2)の導出をよく見ると、やったことは「応力とひずみの掛け算を断面で積分し、断面定数にまとめた」だけでした。同じ操作を他の断面力についても繰り返せば、対応する項が同じ形で出てきます。

軸力：軸応力は断面に一様に分布するので $\sigma = N/A$、ひずみは $\varepsilon = N/(EA)$。仮想側は $\bar\sigma = n/A$ です。断面積分すると

$$ \int_A \frac{n}{A}\cdot\frac{N}{EA}\,dA = \frac{nN}{EA} $$

となり、軸方向に積分して $\int nN/(EA)\,dx$ の項が加わります。

ねじり：円形断面のねじりでは、せん断応力は中心からの距離 $r$ に比例して $\tau = Tr/J$（$J$ は断面二次極モーメント）、せん断ひずみは $\gamma = \tau/G$。曲げのときとまったく同じ流れで $\int_A r^2 dA = J$ が現れ、$\int Tt/(GJ)\,dx$ の項が出ます。

せん断力：これだけは少し厄介です。せん断応力は断面内で一様ではなく（矩形断面なら中立軸で最大の放物線分布）、$\tau = VQ(y)/(I b(y))$ という形をしています。素直に断面積分すると

$$ \int_A \frac{vQ}{Ib}\cdot\frac{VQ}{GIb}\,dA = \frac{vV}{GA}\underbrace{\left[\frac{A}{I^2}\int_A \left(\frac{Q}{b}\right)^2 dA\right]}_{\textstyle \kappa} $$

括弧の中は断面形状だけで決まる無次元量で、せん断補正係数（形状係数）$\kappa$ と呼ばれます。矩形断面なら $\kappa = 6/5 = 1.2$、円形断面なら $\kappa = 10/9 \approx 1.11$、I形断面ではおおむね「全断面積 ÷ ウェブ面積」に近い値になります。これを使うと $\int \kappa Vv/(GA)\,dx$ の項が加わります。

以上をまとめると、単位荷重法の一般形は次のようになります。

$$ \begin{equation} 1\cdot\delta = \int \frac{Mm}{EI}\,dx + \int \frac{Nn}{EA}\,dx + \int \frac{\kappa Vv}{GA}\,dx + \int \frac{Tt}{GJ}\,dx \end{equation} $$

4つの項はすべて「実の内力 × 仮想の内力 ÷ 対応する剛性」という同じ顔つきをしています。曲げなら曲げ剛性 $EI$、軸なら伸び剛性 $EA$、せん断ならせん断剛性 $GA$、ねじりならねじり剛性 $GJ$。この統一感は偶然ではなく、どの項も「応力 × ひずみ」を断面で潰した結果だからです。

さらに、温度変化や製作誤差といった「荷重によらない変形」も同じ枠組みで足せます。温度勾配 $\Delta T$ による曲率変化を $\alpha\Delta T/h$、一様温度上昇による伸びひずみを $\alpha T_0$ とすれば

$$ \delta_{\text{温度}} = \int m\,\frac{\alpha\,\Delta T}{h}\,dx + \int n\,\alpha T_0\,dx $$

が加わるだけです。式(1)の右辺は「仮想の力 × 実の変形」であって、その変形の原因を問わないからです。宇宙構造物のように日陰・日向で数百Kの温度差がつく場面では、この項が力学荷重より支配的になることさえあります。

こうして4項が出揃いましたが、実務ではたいてい第1項しか使いません。なぜ他を無視してよいのか、オーダーで確かめておきましょう。

なぜ細長いはりでは曲げ項だけでよいのか

「せん断変形は無視できる」と教科書はさらっと書きますが、無視できるかどうかははりの細長さ次第です。定量的に見積もってみます。

長さ $L$ の片持ちはりの先端に集中荷重 $P$ を載せ、先端たわみを2つの項で比べます。先端から測った座標を $x$ とすると、実状態は $M = -Px$、$V = P$、仮想状態（先端に単位荷重）は $m=-x$、$v=1$ です。

曲げ項は

$$ \delta_b = \int_0^L \frac{(-Px)(-x)}{EI}dx = \frac{PL^3}{3EI} $$

せん断項は $V$ も $v$ も一定なので、そのまま積分できて

$$ \delta_s = \int_0^L \frac{\kappa P \cdot 1}{GA}dx = \frac{\kappa P L}{GA} $$

比を取ります。$P$ と $L$ が一部約分されて

$$ \frac{\delta_s}{\delta_b} = \frac{\kappa P L/(GA)}{PL^3/(3EI)} = \frac{3\kappa E}{G}\cdot\frac{I}{A}\cdot\frac{1}{L^2} $$

ここで $I/A$ は断面回転半径の2乗 $i^2$ そのもので、長さの2乗の次元を持ちます。つまりこの比は

$$ \frac{\delta_s}{\delta_b} = 3\kappa\,\frac{E}{G}\left(\frac{i}{L}\right)^2 $$

と書け、細長比の逆数の2乗でスケールします。$L$ が2倍になれば比は4分の1になる、というのが本質です。

具体的な数字を入れましょう。矩形断面（幅 $b$、せい $h$）なら $I/A = h^2/12$、$\kappa=6/5$。等方材でポアソン比 $\nu=0.3$ なら $E/G = 2(1+\nu) = 2.6$。これらを代入すると

$$ \frac{\delta_s}{\delta_b} = 3\cdot\frac{6}{5}\cdot 2.6\cdot\frac{1}{12}\left(\frac{h}{L}\right)^2 = 0.78\left(\frac{h}{L}\right)^2 $$

という覚えやすい式になります。スパン／せい比 $L/h$ ごとに数値を並べると次の通りです。

両対数グラフでプロットすると、せん断項と曲げ項の比はきれいな傾き $-2$ の直線になり、実際に個別積分した実測値（緑）とオーダー評価 $0.78(h/L)^2$（黒破線）が完全に重なります。左のピンク帯（$L/h<5$）ではせん断寄与が3%以上あり、$L/h=2$ では19.5%にも達するのでTimoshenko梁が必要です。一方、右の緑帯（$L/h\ge10$）では1%の目安線を下回っており、一般的なはりでは曲げ項だけで実用上十分だと一目でわかります。

一般的な鋼はりの $L/h$ はおよそ15〜25、木造の梁でも10〜20程度ですから、せん断項は1%を大きく下回ります。材料特性のばらつきや支持条件の理想化による誤差のほうがはるかに大きく、わざわざ計算する意味がありません。これが「はりでは曲げ項だけ」の正体です。

逆に言えば、$L/h$ が5を切るような深いはり、サンドイッチパネルのようにコア材のせん断剛性 $G$ が極端に低い構造、コンクリートの短スパン梁などでは、せん断項が効いてきます。この場合はEuler–Bernoulli梁ではなくTimoshenko梁の枠組みに切り替えるのが素直です。単位荷重法の一般形(3)は、そのまま第3項を残すだけで対応できます。

理論の骨格が固まったので、次は実際に手を動かして古典的な公式を再現してみましょう。

例1：片持ちはり先端の集中荷重

長さ $L$、曲げ剛性 $EI$ 一定の片持ちはり。先端に下向き $P$ を載せたときの先端たわみ $\delta$ を求めます。

手順1：座標を選ぶ。 積分が簡単になるよう、自由端を原点にして $x$ を固定端に向かって測ります。こうすると $x=0$ でモーメントがゼロになり、積分定数のような煩わしさが消えます。

手順2：実状態のBMD。 位置 $x$ の断面で切って自由端側の釣り合いを取ると、その断面には距離 $x$ だけ離れたところに $P$ が載っています。上面引張（hogging）を負とする規約で

$$ M(x) = -P x \qquad (0 \le x \le L) $$

手順3：仮想状態のBMD。 求めたいのは先端の鉛直たわみなので、先端に下向きの単位荷重1を置きます。形は実状態とまったく同じで、$P$ が $1$ に変わるだけです。

$$ m(x) = -1 \cdot x = -x $$

手順4：積分する。 式(2)に代入します。マイナス同士が掛かって正になることに注意してください。

$$ \delta = \int_0^L \frac{(-Px)(-x)}{EI}dx = \frac{P}{EI}\int_0^L x^2 dx = \frac{P}{EI}\cdot\frac{L^3}{3} = \frac{PL^3}{3EI} $$

上段の $M=-Px$ と中段の $m=-x$ はどちらも同じ直線で、大きさが $P$ から $1$ に変わっただけであることが見て取れます。下段はその積を $EI$ で割った $Mm/EI$ で、負×負が正になるため全区間で正の放物線になり、この面積がそのまま先端たわみです。$L=3$ m、$P=10$ kN、$EI=4\times10^6$ N·m$^2$ で数値を入れると 22.500 mm となり、公式 $PL^3/3EI$ の値と完全に一致しています。

見慣れた公式が出ました。弾性曲線方程式を2回積分して境界条件を2本使う手順と比べると、積分1回で終わっているのがわかります。しかも $\delta>0$ なので、先端は単位荷重の向き（下向き）に変位します。物理的にも当然です。

同じ要領で、等分布荷重 $w$ の場合も一瞬です。実状態は $M=-wx^2/2$、仮想状態は同じ $m=-x$ なので

$$ \delta = \int_0^L \frac{(-wx^2/2)(-x)}{EI}dx = \frac{w}{2EI}\cdot\frac{L^4}{4} = \frac{wL^4}{8EI} $$

こちらも教科書通りです。実状態が変わっても仮想状態は使い回せる——これが単位荷重法のもう1つの利点です。同じ点の変位を、荷重パターンを変えて何通りも調べたいとき、$m(x)$ は一度描けば済みます。

片持ちはりは区間が1つしかない易しい例でした。次は区間が分かれるケースを見ます。

例2：単純支持ばりの等分布荷重、中央たわみ

スパン $L$ の単純支持ばり（左端ピン、右端ローラ）に等分布荷重 $w$ が載っています。中央 $C$ の鉛直たわみを求めます。

実状態のBMD。 左端からの座標を $x$ とします。両端の反力はそれぞれ $wL/2$ なので

$$ M(x) = \frac{wL}{2}x – \frac{w x^2}{2} \qquad (0\le x\le L) $$

これは全区間で1本の式（放物線）で書けます。中央で最大値 $wL^2/8$ を取ります。

仮想状態のBMD。 中央に下向きの単位荷重1を置くと、両端反力は $1/2$ ずつ。荷重点で式が切り替わるので、区間を2つに分ける必要があります。

$$ m(x) = \begin{cases} \dfrac{x}{2} & (0\le x\le L/2)\\[8pt] \dfrac{L-x}{2} & (L/2\le x\le L) \end{cases} $$

三角形のBMDで、中央で最大値 $L/4$ です。

積分する。 実状態 $M$ は左右対称、仮想状態 $m$ も左右対称なので、積 $Mm$ も対称です。したがって左半分だけ積分して2倍すれば済みます。

$$ \delta_C = 2\int_0^{L/2}\frac{1}{EI}\left(\frac{wL}{2}x-\frac{wx^2}{2}\right)\frac{x}{2}\,dx $$

括弧を展開して $x$ の多項式にまとめると

$$ = \frac{1}{EI}\int_0^{L/2}\left(\frac{wL x^2}{2}-\frac{wx^3}{2}\right)dx = \frac{w}{2EI}\left[\frac{Lx^3}{3}-\frac{x^4}{4}\right]_0^{L/2} $$

上限 $x=L/2$ を代入します。$x^3 = L^3/8$、$x^4=L^4/16$ なので

$$ = \frac{w}{2EI}\left(\frac{L^4}{24}-\frac{L^4}{64}\right) $$

括弧の中を通分します。24と64の最小公倍数は192で、$1/24 = 8/192$、$1/64=3/192$ ですから差は $5/192$。したがって

$$ \delta_C = \frac{w}{2EI}\cdot\frac{5L^4}{192} = \frac{5wL^4}{384EI} $$

上段の実状態は全区間で1本の放物線ですが、中段の仮想状態は荷重点である中央でぴたりと折れており、ここが区間を2つに分けなければならない理由です。下段の積 $Mm/EI$ も中央で折れ点を持ち、左右対称なので左半分の面積を2倍すれば済むことが図から確認できます。面積の合計 1.318 mm は、公式 $5wL^4/384EI$ に $L=3$ m、$w=5$ kN/m、$EI=4\times10^6$ N·m$^2$ を入れた値と一致します。

構造力学で最も有名な公式の1つが出ました。分母の384という半端な数字が、$2\times 192$ から来ていたことがわかります。

具体的な数値も入れておきましょう。$L=3$ m、$w=5$ kN/m、$E=200$ GPa、$I=2000$ cm$^4$（$EI = 4\times10^6$ N·m$^2$）とすると

$$ \delta_C = \frac{5\times 5000\times 3^4}{384\times 4\times 10^6} = 1.318\times 10^{-3}\ \text{m} = 1.318\ \text{mm} $$

スパンの $1/2276$ です。一般的なたわみ制限（$L/300$ など）に対して十分小さく、この断面で問題ないと判断できます。

区間分割さえ間違えなければ、あとは機械的な積分です。ただし毎回この積分を手でやるのは面倒なので、実務では次の「表」を使います。

積分表（Vereshchaginの図式積分）

$M$ と $m$ はどちらもBMD、つまり図形です。そして単位荷重法の積分は「2つの図形の積の積分」に過ぎません。片方が直線なら、この積分は非常に簡単な形になります。

$m$ が区間 $[0,L]$ で直線だとします。すると

$$ \int_0^L M m\,dx = \int_0^L M(x)\,m(x)\,dx = m(x_G)\int_0^L M\,dx = A_M \cdot m(x_G) $$

ここで $A_M$ は $M$ 図の面積、$x_G$ はその図心の位置です。2行目の変形は、$m$ が直線なら $\int Mm\,dx$ が「$M$ 図の1次モーメント」に比例することから従います。つまり

一方の図が直線なら、積分は「もう一方の図の面積 × 直線図の図心位置での値」になる。

これがVereshchaginの定理（図式積分法）です。よく使う組み合わせを表にしておきます（$M_0$、$m_0$ はそれぞれの図の代表値）。

例1の片持ちはり（等分布）で確かめます。$M$ 図は自由端で0、固定端で $M_0=wL^2/2$ の2次放物線（頂点が自由端側）、$m$ 図は自由端0から固定端 $m_0=L$ の三角形。表の該当セルは $LM_0m_0/4$ なので

$$ \delta = \frac{1}{EI}\cdot\frac{L\cdot(wL^2/2)\cdot L}{4} = \frac{wL^4}{8EI} $$

先ほどの積分結果と一致しました。

上段の $M$ 図は放物線なので図心は自由端から $3L/4=2.250$ m の位置にあり、下段の $m$ 図はそこで $-2.250$ m の値を取ります。この2つを掛けた「面積 × 図心位置の値」が $50.6250$ kN·m$^3$ で、素直に $\int Mm\,dx$ を数値積分した値と小数点以下4桁まで一致しました。曲がっている側は面積と図心だけ、直線側は1点の値だけあればよいので、BMDを描いた時点で積分が終わっているのがVereshchaginの威力です。

なお、$m$ が全区間で直線でない場合は区間を分けます。例2の単純ばり中央たわみは $m$ が中央で折れているので、左右に分けて足す必要があります。この場合の答えは表からは直接読めず、$5LM_0m_0/12$（$M_0=wL^2/8$、$m_0=L/4$）という別のセルになります。実際 $5L(wL^2/8)(L/4)/12 = 5wL^4/384$ で、確かに一致します。「表を引く前に $m$ が折れていないか確かめる」のが鉄則です。

はりの話が固まったので、次はまったく形の違う構造——トラス——に同じ式を適用します。

トラスへの拡張

トラスの部材は、両端ピンで曲げを伝えないため、軸力しか働きません。しかも各部材内で軸力は一定です。ということは、式(3)のうち生き残るのは軸力の項だけで、しかも積分が単なる掛け算になります。

部材 $i$ について、実状態の軸力を $N_i$（引張正）、仮想状態の軸力を $n_i$、部材長を $L_i$、断面剛性を $E_iA_i$ とすると

$$ \int_{\text{部材}i} \frac{N_i n_i}{E_iA_i}dx = \frac{N_i n_i L_i}{E_i A_i} $$

全部材について足せば、トラスの単位荷重法の公式が得られます。

$$ \begin{equation} 1\cdot\delta = \sum_{i=1}^{M} \frac{N_i\, n_i\, L_i}{E_i A_i} \end{equation} $$

積分記号すら消えて、ただの掛け算と足し算になりました。手順も明快です。

1. 実荷重で全部材の軸力 $N_i$ を求める（節点法・断面法、または直接剛性法）
2. すべての実荷重を外し、知りたい点・方向に単位荷重1だけを置いて $n_i$ を求める
3. $N_i n_i L_i/(E_iA_i)$ を表にして縦に足す

ここで注意したいのが、$N_i n_i L_i/(E_iA_i)$ という各項の意味です。$N_iL_i/(E_iA_i)$ は部材 $i$ の実際の伸び $e_i$ ですから、式(4)は

$$ \delta = \sum_i n_i\, e_i $$

と書き直せます。つまり「各部材の伸びが、知りたい変位にどれだけ寄与するか」を $n_i$ が重み付けしている。$n_i$ が大きい部材ほど、その伸びが着目点の変位に直結する——という設計上の読み方ができます。逆に $n_i=0$ の部材は、いくら伸びても着目点は動きません。トラスの部材断面を増やして変位を減らしたいとき、どの部材を太くすれば効くかが $n_i$ を見れば一目でわかるわけです。これは感度解析そのものです。

注意点として、$n_i$ は $N_i$ に比例するとは限りません。荷重が1つだけでその位置・方向に単位荷重を置いた場合に限り $n_i = N_i/P$ になりますが、荷重が複数あったり、荷重の載っていない点の変位を求めたりすると、$n_i$ は $N_i$ とまったく違うパターンになります。後半のPython実装では、あえて荷重を2つ載せてこの点を確認します。

トラスの式を眺めていると、ある美しい性質に気づきます。次節でそれを取り出しましょう。

Maxwellの相反定理

式(2)や式(4)で、実状態の荷重も単位荷重にしてみるとどうなるでしょうか。

点 $j$ に単位荷重を置いたときの曲げモーメントを $m_j(x)$、点 $i$ に単位荷重を置いたときを $m_i(x)$ とします。「点 $j$ の単位荷重によって点 $i$ に生じる変位」を $\delta_{ij}$ と書くと、式(2)で実状態を $M = m_j$、仮想状態を $m = m_i$ として

$$ \begin{equation} \delta_{ij} = \int \frac{m_i(x)\,m_j(x)}{EI}\,dx \end{equation} $$

この右辺を見てください。$i$ と $j$ について完全に対称です。掛け算の順序を入れ替えても値は変わりません。したがって

$$ \begin{equation} \delta_{ij} = \delta_{ji} \end{equation} $$

これがMaxwellの相反定理です。言葉にすると

点 $j$ に単位荷重を載せたときの点 $i$ の変位は、点 $i$ に単位荷重を載せたときの点 $j$ の変位に等しい。

直感的には驚くべき主張です。橋の1/4点に人が立ったときの中央の沈下量と、中央に同じ人が立ったときの1/4点の沈下量が、まったく同じだというのですから。構造が非対称でも、断面が場所ごとに違っても、支持条件が複雑でも成り立ちます。線形弾性でありさえすればよいのです。

なぜこんなことが起きるのか。式(5)の対称性がすべてです。より深いところでは、線形弾性体のひずみエネルギーが変位の2次形式 $U=\frac{1}{2}\bm{u}^\top\bm{K}\bm{u}$ で書け、その2階微分である剛性行列 $\bm{K}$ が対称だからです。柔性行列 $\bm{F}=\bm{K}^{-1}$ の成分がまさに $\delta_{ij}$ で、対称行列の逆行列は対称なので $\delta_{ij}=\delta_{ji}$。数学的には当たり前ですが、物理的には非自明な結論です。

一般化した形（Bettiの相反定理）もあります。2つの荷重系 $\{P_i\}$ と $\{Q_j\}$ について、

$$ \sum_i P_i\, u_i^{(Q)} = \sum_j Q_j\, u_j^{(P)} $$

「荷重系Pが、荷重系Qによる変位に対してする仕事」＝「荷重系Qが、荷重系Pによる変位に対してする仕事」。単位荷重法はこの特別な場合（両方が単位荷重1つずつ）に過ぎません。

実務での使いどころは3つあります。第一に計算量の削減：知りたい変位が計算しにくい配置なら、荷重と測定点を入れ替えてしまえばよい。第二に影響線の作成：橋梁で「荷重がどこにあるとき中央のたわみが最大か」を知りたいとき、中央に単位荷重を置いた1回の解析で影響線が丸ごと得られます。第三に数値解析の検証：FEMソルバが正しく実装されていれば $\delta_{ij}=\delta_{ji}$ が機械精度で成り立つはずで、成り立たなければどこかにバグがあります。

相反定理は「同じ積分の対称性」から出てきました。次は、その積分を未知数を含む形で書いて、不静定問題を解きます。

1次不静定はりを解く — 整合条件

ここまでは静定構造、つまり釣り合い式だけで内力が決まる構造を扱ってきました。しかし実際の構造の多くは不静定です。例えば「一端固定・他端単純支持」のはり（propped cantilever）は、未知反力が4つ（固定端の反力・モーメント、支点反力、水平反力）あるのに釣り合い式は3本しかなく、1つ足りません。この「足りない本数」を不静定次数と呼びます。

不静定を解く発想は次の通りです。

1. 余分な拘束（冗長力）を1つ選んで外し、静定な「基本系」を作る
2. 外した拘束の代わりに、未知の力 $X$ を作用させる
3. 本来その点で満たされていた変位条件（例：支点だから鉛直変位ゼロ）を書き下し、$X$ について解く

3番目のステップが整合条件です。そしてこの変位を計算する道具が、まさに単位荷重法です。

定式化

冗長力 $X$ の作用点・方向における変位を、重ね合わせで書きます。線形なので、実荷重による寄与と $X$ による寄与を足せます。

$$ \delta = \delta_{10} + X\,\delta_{11} $$

• $\delta_{10}$：基本系に実荷重だけを載せたときの、$X$ の作用点・方向の変位
• $\delta_{11}$：基本系に$X=1$ だけを載せたときの、同じ点・方向の変位

支点なら $\delta=0$ なので、整合条件は

$$ \begin{equation} \delta_{10} + X\,\delta_{11} = 0 \quad\Longrightarrow\quad X = -\frac{\delta_{10}}{\delta_{11}} \end{equation} $$

$\delta_{10}$ も $\delta_{11}$ も単位荷重法で計算できます。基本系での実荷重によるBMDを $M_0$、$X=1$ によるBMDを $m_1$ とすれば

$$ \delta_{10} = \int\frac{M_0\,m_1}{EI}dx, \qquad \delta_{11} = \int\frac{m_1^2}{EI}dx $$

ここで見逃せないのが $\delta_{11}$ の形です。被積分関数が2乗なので、$\delta_{11}$ は必ず正（$EI>0$ なら）。つまり分母がゼロになることはなく、$X$ は常に一意に決まります。これは「柔性行列が正定値である」ことの1次元版で、構造が安定であることの数学的な保証になっています。

例：一端固定・他端単純支持のはり

スパン $L$、等分布荷重 $w$、$EI$ 一定。左端Aが固定、右端Bが単純支持です。

基本系の選択。 Bの支点を取り除いて、Aで固定された片持ちはりにします。冗長力は $X$（B点の鉛直反力、上向き正）。

$\delta_{10}$ の計算。 自由端Bを原点に $x$ を取ります。等分布荷重による片持ちはりのBMDは $M_0 = -wx^2/2$。上向き単位荷重によるBMDは $m_1 = +x$（符号規約を揃えていることに注意）。したがって

$$ \delta_{10} = \int_0^L \frac{(-wx^2/2)(x)}{EI}dx = -\frac{w}{2EI}\cdot\frac{L^4}{4} = -\frac{wL^4}{8EI} $$

負符号は「単位荷重（上向き）と逆、つまり下向きに変位する」ことを意味します。片持ちはりの先端が下がるのは当然です。

$\delta_{11}$ の計算。

$$ \delta_{11} = \int_0^L \frac{x^2}{EI}dx = \frac{L^3}{3EI} $$

整合条件を解く。 式(7)に代入すると $EI$ と $L^3$ が約分されて

$$ X = -\frac{-wL^4/(8EI)}{L^3/(3EI)} = \frac{wL^4}{8EI}\cdot\frac{3EI}{L^3} = \frac{3wL}{8} $$

有名な結果 $R_B = 3wL/8$ が出ました。残りの反力は釣り合いから $R_A = wL – 3wL/8 = 5wL/8$ です。

最終的なBMD。 重ね合わせで $M(x) = M_0(x) + X m_1(x)$ を作ります（$x$ はBからの距離）。

$$ M(x) = -\frac{wx^2}{2} + \frac{3wL}{8}x $$

固定端 $x=L$ では

$$ M_A = -\frac{wL^2}{2}+\frac{3wL^2}{8} = -\frac{wL^2}{8} $$

上面引張のモーメント（hogging）が $wL^2/8$ 生じます。一方、最大の正モーメント（sagging）は $dM/dx = -wx + 3wL/8 = 0$ から $x = 3L/8$ で

$$ M_{\max}^{+} = -\frac{w}{2}\left(\frac{3L}{8}\right)^2+\frac{3wL}{8}\cdot\frac{3L}{8} = -\frac{9wL^2}{128}+\frac{9wL^2}{64} = \frac{9wL^2}{128} $$

ここが面白いところです。同じはりを単純支持にすれば最大モーメントは $wL^2/8 = 16wL^2/128$ でしたから、中央付近の正モーメントが $9/16 = 56\%$ に減った代わりに、固定端に $wL^2/8$ の負モーメントが現れています。これがモーメントの再分配で、不静定構造の本質的な利点です。断面を最大モーメントで決めるなら、単純支持と同じ $wL^2/8$ で済むうえに、たわみは大幅に減ります。

たわみ。 重ね合わせで $\delta(x) = \dfrac{w x^2 (L-x)(3L-2x)}{48EI}$（$x$ は固定端Aから測る）となり、最大値は $x = \frac{15-\sqrt{33}}{16}L \approx 0.5785L$ で

$$ \delta_{\max} = \frac{wL^4}{184.6\,EI} \approx \frac{wL^4}{185\,EI} $$

単純支持の $5wL^4/384EI = wL^4/76.8EI$ と比べると、約41.6%にまで減っています。一端を固定するだけでたわみが半分以下になる——不静定化の効果がよくわかります。

上段のBMDでは、単純支持（灰）の山ひとつだった分布が、一端を固定すると正モーメント 3.164 kN·m と固定端の負モーメント $-5.625$ kN·m に分かれています。最大絶対値はどちらも 5.625 kN·m のままで断面設計の要求は変わりませんが、モーメントが構造全体に配り直された点が本質です。下段のたわみでは最大値が 1.3184 mm から 0.5484 mm（41.6%）へ落ち、しかも最大点が中央 1.5 m から固定端と反対側の 1.735 m ($0.5783L$) へずれていることが確認できます。

理論はここまでです。以降はすべてをPythonで検証します。

Pythonでの実装1：sympyで公式を再現する

まずは手計算した公式が正しいか、記号積分で確かめます。ポイントは区間ごとに $M$ と $m$ を立てて、素直に足すだけであること。プログラムがそのまま手順書になります。

import sympy as sp

# 記号の定義（すべて正の実数）
x, L, P, w, E, I = sp.symbols('x L P w E I', positive=True)

def unit_load(segments):
    """segments = [(M式, m式, 下限, 上限), ...] を受け取り δ = Σ∫ Mm/EI dx を返す"""
    return sp.simplify(sum(sp.integrate(M * m / (E * I), (x, a, b))
                           for M, m, a, b in segments))

# --- 例1: 片持ちはり、先端集中荷重 P（xは自由端から測る）---
d1 = unit_load([(-P * x, -x, 0, L)])
print("片持ち・先端集中:", d1)                    # P*L**3/(3*E*I)

# --- 例1b: 片持ちはり、等分布荷重 w ---
d2 = unit_load([(-w * x**2 / 2, -x, 0, L)])
print("片持ち・等分布  :", d2)                    # L**4*w/(8*E*I)

# --- 例2: 単純支持ばり、等分布荷重、中央たわみ（xは左端から）---
M_udl = w * L * x / 2 - w * x**2 / 2
d3 = unit_load([(M_udl, x / 2, 0, L / 2),
                (M_udl, (L - x) / 2, L / 2, L)])
print("単純ばり・中央  :", d3)                    # 5*L**4*w/(384*E*I)

# --- 例3: 単純支持ばり、等分布荷重、左端のたわみ角（単位モーメントを置く）---
d4 = unit_load([(M_udl, 1 - x / L, 0, L)])
print("単純ばり・端回転:", d4)                    # L**3*w/(24*E*I)

出力は順に P*L**3/(3*E*I)、L**4*w/(8*E*I)、5*L**4*w/(384*E*I)、L**3*w/(24*E*I) となり、手計算した4つの公式すべてと一致します。特に注目したいのは最後のたわみ角です。単位「モーメント」を置いた仮想状態のBMDは $m = 1-x/L$ という三角形になり、同じ関数 unit_load に渡すだけで $\theta = wL^3/(24EI)$ が得られました。求めたい量が変位でも回転角でも、コードは1文字も変わりません。これが「共役な単位荷重を置く」という一般則の威力です。

もう1つの読みどころは、例2で区間を2つに分けている点です。$M$ 側は全区間で同じ放物線ですが、$m$ 側が中央で折れるため、タプルを2つ渡しています。単位荷重法で最も間違えやすいのがこの区間分割で、コードにすると「どこで折れるか」が明示されるため、ミスが減ります。

Pythonでの実装2：たわみ曲線を丸ごと描く

単位荷重法は「1点だけ」を求める手法ですが、単位荷重の位置 $a$ をパラメータにして記号のまま解けば、たわみ曲線全体が一度に手に入ります。弾性曲線方程式を解くのとは違う道筋で同じ答えに到達することを確かめましょう。

import sympy as sp
import numpy as np
import matplotlib, matplotlib.pyplot as plt

for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

x, a, L, w, E, I = sp.symbols('x a L w E I', positive=True)

# 実状態：単純ばり＋等分布荷重
M = w * L * x / 2 - w * x**2 / 2
# 仮想状態：位置 a に単位荷重（反力は (L-a)/L と a/L）
m_left, m_right = (L - a) / L * x, a / L * (L - x)

delta_a = sp.simplify(sp.integrate(M * m_left / (E * I), (x, 0, a))
                      + sp.integrate(M * m_right / (E * I), (x, a, L)))
print("δ(a) =", sp.factor(delta_a))
print("a=L/2 :", sp.simplify(delta_a.subs(a, L / 2)))

得られる式は $\delta(a) = \dfrac{w\,a\,(L^3 – 2La^2 + a^3)}{24EI}$ で、$a=L/2$ を代入すると 5*L**4*w/(384*E*I) に戻ります。これは弾性曲線方程式を4回積分して境界条件4本を解いた結果とまったく同じ多項式です。単位荷重の位置を記号のまま残すという一手間だけで、点の答えが曲線の答えに化けました。

続いて数値を入れて描画し、後の比較に使う関数を用意します。

Lv, wv, EIv = 3.0, 5e3, 4e6        # L=3m, w=5kN/m, EI=4e6 N·m^2
f_exact = sp.lambdify(a, delta_a.subs({L: Lv, w: wv, E: 1, I: EIv}), 'numpy')

xs = np.linspace(0, Lv, 201)
plt.figure(figsize=(9, 4))
plt.plot(xs, 1e3 * f_exact(xs), lw=2, label='単位荷重法によるたわみ曲線')
plt.plot([Lv / 2], [1e3 * f_exact(Lv / 2)], 'o', ms=9,
         label=f'中央 {1e3*f_exact(Lv/2):.3f} mm （$5wL^4/384EI$）')
plt.gca().invert_yaxis()
plt.xlabel('支点からの距離 $x$ [m]'); plt.ylabel('たわみ（下向き正）[mm]')
plt.title('単純支持ばり・等分布荷重のたわみ曲線（単位荷重の位置を動かして構成）')
plt.legend(); plt.grid(alpha=0.3); plt.tight_layout(); plt.show()

上段では単位荷重の位置 $a$ を 0.75 m・1.50 m・2.25 m と動かしており、そのたびに三角形の頂点が移動して仮想状態のBMDが別物になっていることがわかります。下段はその各 $a$ に対する積分値をつないだもので、$a$ を記号のまま残して1回積分するだけでたわみ曲線が丸ごと得られています。$a=0.75$ m と $a=2.25$ m で同じ 0.939 mm になるのは構造が左右対称だからで、中央の 1.318 mm が $5wL^4/384EI$ と一致します。

描かれる曲線は左右対称で、中央で最大値 1.318 mm を取ります。手計算した $5wL^4/384EI = 1.318$ mm と小数点以下3桁まで一致しました。曲線が両端でゼロ、両端付近で傾きが最大になるという形も、単純支持の物理的な振る舞いと合致しています。「1点の変位を求める手法」を201回繰り返しただけで曲線になるという当たり前の事実が、図として目に見えるのが面白いところです。

Pythonでの実装3：曲げ項とせん断項の競争

次に、オーダー評価で議論した「細長さと項の寄与」を数値で確かめます。矩形断面の片持ちはりについて、$L/h$ を変えながら両項の比を計算します。

import numpy as np
import matplotlib, matplotlib.pyplot as plt

for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

E, nu, kappa = 200e9, 0.3, 6 / 5      # 鋼、矩形断面
G = E / (2 * (1 + nu))
P, b, h = 10e3, 0.10, 0.30            # 荷重10kN、断面 100mm×300mm
A, I = b * h, b * h**3 / 12

ratios = np.linspace(2, 50, 400)      # L/h
Ls = ratios * h
d_b = P * Ls**3 / (3 * E * I)         # 曲げ寄与
d_s = kappa * P * Ls / (G * A)        # せん断寄与

plt.figure(figsize=(9, 4.5))
plt.loglog(ratios, 100 * d_s / d_b, lw=2, label='せん断項 / 曲げ項（実測）')
plt.loglog(ratios, 100 * 0.78 / ratios**2, '--', lw=1.5,
           label=r'オーダー評価 $0.78\,(h/L)^2$')
plt.axhline(1, color='gray', ls=':', label='1 % の目安')
for r in (5, 10, 20):
    plt.plot([r], [100 * 0.78 / r**2], 'o', ms=7)
    plt.annotate(f'L/h={r}\n{100*0.78/r**2:.2f}%', (r, 100 * 0.78 / r**2),
                 textcoords='offset points', xytext=(8, 6), fontsize=9)
plt.xlabel('スパン／せい比 $L/h$'); plt.ylabel('せん断寄与の割合 [%]')
plt.title('細長いはりでせん断変形が無視できる理由')
plt.legend(); plt.grid(alpha=0.3, which='both'); plt.tight_layout(); plt.show()

両対数プロットで直線になり、傾きが $-2$、つまり比が $(L/h)^{-2}$ でスケールすることが確認できます。実測値（曲げ・せん断を個別に積分した値）とオーダー評価 $0.78(h/L)^2$ は完全に重なりました。数値としては $L/h=5$ で 3.12%、$L/h=10$ で 0.78%、$L/h=20$ で 0.195%。一般的なはりの範囲（$L/h \ge 10$）では1%を切っており、材料定数のばらつき以下です。

逆に $L/h=2$ 付近では 19.5% にも達します。深いはりやサンドイッチ構造でTimoshenko梁理論が必要になるのは、この領域です。「無視できる」は絶対的な性質ではなく、細長比が決める相対的な判断だという当たり前のことが、図としてはっきり見えます。

Pythonでの実装4：トラスの節点変位

はりを離れて、トラスに式(4)を適用します。まず検証用に直接剛性法（トラスFEM）のソルバを用意し、それを使って実状態の軸力 $N_i$ と仮想状態の軸力 $n_i$ の両方を求めます。

対象は、スパン4 m・高さ1.5 mのワーレントラスです。節点は $N_0(0,0)$（ピン支点）、$N_1(2,0)$、$N_2(4,0)$（ローラ支点）、$N_3(1,1.5)$、$N_4(3,1.5)$ の5つ。部材は7本、拘束は3つなので $2\times5-3=7$、静定です。荷重はあえて2つ（$N_1$ に10 kN下向き、$N_4$ に6 kN下向き）載せて、$n_i$ が $N_i$ に比例しない状況を作ります。

import numpy as np

nodes = np.array([[0, 0], [2, 0], [4, 0], [1, 1.5], [3, 1.5]], float)
mems = [(0, 1), (1, 2), (3, 4), (0, 3), (3, 1), (1, 4), (4, 2)]
EA = 200e9 * 1e-3      # E=200GPa, A=1000mm^2 → 2e8 N
fixed = [0, 1, 5]      # N0の水平・鉛直、N2の鉛直

def solve_truss(F):
    """節点荷重ベクトルFを受け取り、節点変位uと部材軸力Nを返す"""
    n = len(nodes)
    K = np.zeros((2 * n, 2 * n))
    Ls, Ts = [], []
    for (a, b) in mems:
        d = nodes[b] - nodes[a]
        Lm = np.hypot(*d); c, s = d / Lm
        T = np.array([-c, -s, c, s])        # 軸方向への射影ベクトル
        Ls.append(Lm); Ts.append(T)
        dof = [2 * a, 2 * a + 1, 2 * b, 2 * b + 1]
        K[np.ix_(dof, dof)] += EA / Lm * np.outer(T, T)
    free = [i for i in range(2 * n) if i not in fixed]
    u = np.zeros(2 * n)
    u[free] = np.linalg.solve(K[np.ix_(free, free)], F[free])
    N = np.array([EA / Ls[i] * Ts[i] @ u[[2 * a, 2 * a + 1, 2 * b, 2 * b + 1]]
                  for i, (a, b) in enumerate(mems)])
    return u, N, np.array(Ls)

このソルバは各部材の剛性 $\frac{EA}{L}\bm{T}\bm{T}^\top$ を全体行列に足し込むだけの標準的な直接剛性法です。射影ベクトル $\bm{T}=(-c,-s,c,s)$ が、節点変位から部材の伸びを取り出す役割を果たしています（$e_i = \bm{T}_i \cdot \bm{u}_i$）。これで実状態も仮想状態も同じ関数で解けます。

# 実状態：N1に10kN下向き、N4に6kN下向き
F = np.zeros(10); F[3] = -10e3; F[9] = -6e3
u, N, Ls = solve_truss(F)

# 仮想状態：N1に鉛直下向きの単位荷重（1 N）だけ
F1 = np.zeros(10); F1[3] = -1.0
_, n1, _ = solve_truss(F1)

contrib = N * n1 * Ls / EA                # 各部材の寄与 [m]
print("部材    N[kN]      n1        NnL/EA[mm]")
for i, (a, b) in enumerate(mems):
    print(f"{a}-{b}   {N[i]/1e3:8.3f}  {n1[i]:7.4f}   {1e3*contrib[i]:8.5f}")
print(f"\n単位荷重法  δ = {1e3*contrib.sum():.5f} mm（下向き）")
print(f"直接剛性法  δ = {-1e3*u[3]:.5f} mm（下向き）")

出力を整理すると次の表になります。

単位荷重法の答え 0.26259 mm に対し、直接剛性法の答えも 0.26259 mm。有効数字16桁まで一致します（相対差は $10^{-16}$ オーダー、つまり丸め誤差そのものです）。まったく異なる2つの理屈——エネルギーの等式と、剛性行列の連立方程式——が同じ数値に到達することは、単位荷重法が近似ではなく厳密な等式であることの何よりの証拠です。

左の実状態と中央の仮想状態を見比べると、下弦の2部材はどちらも $n_1=0.333$ で同じなのに実軸力は 4.33 kN と 6.33 kN で違い、$n$ が $N$ に比例しないことがはっきりします。右の棒グラフは各部材の寄与 $N n L/(EA)$ を降順に並べたもので、部材4–2が 0.06185 mm（23.6%）で最大、部材0–1は 0.01444 mm（5.5%）にとどまります。合計 0.26259 mm は直接剛性法の答えと桁まで一致しており、しかも「どの部材を太くすれば効くか」という感度情報が同時に得られています。

表からもう2つ読み取れます。第一に、$n_1$ の列は $N$ の列に比例していません。例えば部材0–1と1–2はどちらも $n_1=0.3333$ ですが、$N$ は 4.333 kN と 6.333 kN で違います。荷重が2つあるためです。第二に、寄与が最も大きいのは部材4–2（0.06185 mm、全体の24%）で、$N$ も $n_1$ も大きいためです。この部材の断面を太くするのが、$N_1$ のたわみを減らす最も効率のよい手だとわかります。逆に部材0–1の寄与は5.5%しかなく、ここを補強してもほとんど効きません。単位荷重法は変位を出すだけでなく、こうした設計の指針まで一緒に返してくれます。

Pythonでの実装5：相反定理の数値検証

相反定理 $\delta_{ij}=\delta_{ji}$ を、上のトラスで直接確かめます。$N_1$（下弦中央）と $N_4$（上弦右）という、まったく性質の違う2点を選びます。

import numpy as np

# 1kN の単位荷重をそれぞれの点に単独で載せる
Fa = np.zeros(10); Fa[3] = -1e3     # N1 に下向き1kN
Fb = np.zeros(10); Fb[9] = -1e3     # N4 に下向き1kN
ua, na, Ls = solve_truss(Fa)
ub, nb, _  = solve_truss(Fb)

d41 = -1e3 * ua[9]                  # N1の荷重 → N4のたわみ [mm]
d14 = -1e3 * ub[3]                  # N4の荷重 → N1のたわみ [mm]
d_int = 1e3 * np.sum(na * nb * Ls / EA) / 1e3   # 積分（総和）表現

print(f"δ41 (N1に載せてN4を測る) = {d41:.15f} mm")
print(f"δ14 (N4に載せてN1を測る) = {d14:.15f} mm")
print(f"Σ n_a n_b L /EA          = {d_int:.15f} mm")
print(f"相対差 = {abs(d41-d14)/abs(d41):.3e}")

出力は3つとも 0.010954467580699 mm、相対差は 0.000e+00（浮動小数点で完全一致）です。$N_1$ に1 kNを載せて $N_4$ の沈下を測っても、$N_4$ に1 kNを載せて $N_1$ の沈下を測っても、10.95 µm で寸分違わない。トラスの形は左右非対称ではないものの、$N_1$ は下弦の節点、$N_4$ は上弦の節点で、支持条件からの距離も部材構成もまったく違います。それでも一致するのが相反定理です。

左右のパネルは変形形状がまったく違います。左はN1が大きく沈み込む形、右は上弦の右端が押し下げられて骨組み全体がねじれるような形で、共通点は見当たりません。それでも「載せた点と測る点を入れ替えた2つの値」は 10.954 µm で小数点以下15桁まで一致しており、変形の見た目とは無関係に $\delta_{ij}=\delta_{ji}$ が成り立つことが数値で確かめられます。

3行目の $\sum n_a n_b L/(EA)$ が同じ値になっていることにも注目してください。これは式(5)のトラス版で、両方を単位荷重にした積分が柔性行列の成分そのものであることを示しています。$n_a$ と $n_b$ の掛け算は順序を入れ替えても同じ、だから $\delta_{ij}=\delta_{ji}$——証明が式の見た目そのままです。FEMを自作したとき、この検証を1行入れておくと、剛性行列の組み立てミス（対称性を壊すバグ）を確実に検出できます。

Pythonでの実装6：1次不静定はりと曲げモーメントの再分配

最後に、整合条件で不静定はりを解き、BMDがどう再分配されるかを可視化します。$L=3$ m、$w=5$ kN/m、$EI=4\times10^6$ N·m$^2$ の一端固定・他端単純支持のはりです。

import sympy as sp
import numpy as np
import matplotlib, matplotlib.pyplot as plt

for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

x, L, w, E, I, X = sp.symbols('x L w E I X', positive=True)

# 基本系＝片持ちはり（xは自由端Bから測る）
M0 = -w * x**2 / 2          # 実荷重によるBMD
m1 = x                      # X=1（上向き）によるBMD

d10 = sp.integrate(M0 * m1 / (E * I), (x, 0, L))   # -wL^4/(8EI)
d11 = sp.integrate(m1 * m1 / (E * I), (x, 0, L))   #  L^3/(3EI)
Xsol = sp.simplify(sp.solve(sp.Eq(d10 + X * d11, 0), X)[0])
print("δ10 =", sp.simplify(d10), " δ11 =", sp.simplify(d11), " X =", Xsol)

M_final = sp.simplify(M0 + Xsol * m1)
print("M(x) =", M_final)
print("固定端 M_A =", sp.simplify(M_final.subs(x, L)))
xstar = sp.solve(sp.diff(M_final, x), x)[0]
print("最大正モーメント位置 x =", xstar, " 値 =", sp.simplify(M_final.subs(x, xstar)))

出力は $\delta_{10}=-wL^4/(8EI)$、$\delta_{11}=L^3/(3EI)$、$X=3wL/8$、$M_A=-wL^2/8$、そして $x=3L/8$ で $M^+_{\max}=9wL^2/128$。手計算の結果と完全に一致しました。整合条件はたった1本の1次方程式で、しかも係数はどちらも単位荷重法の積分で出せる——不静定がこれほど単純に片付くことに驚くはずです。

続いてBMDを比較描画します。

Lv, wv, EIv = 3.0, 5e3, 4e6
sub = {L: Lv, w: wv}
xs = np.linspace(0, Lv, 301)

f_ind = sp.lambdify(x, M_final.subs(sub), 'numpy')          # 不静定（xはBから）
M_ss  = wv * Lv * x / 2 - wv * x**2 / 2                     # 単純支持のBMD
f_ss  = sp.lambdify(x, M_ss, 'numpy')

plt.figure(figsize=(9, 4.5))
plt.plot(xs, f_ss(xs) / 1e3, lw=2, label='単純支持（両端ピン・ローラ）')
plt.plot(Lv - xs, f_ind(xs) / 1e3, lw=2, label='一端固定・他端単純支持')
plt.axhline(0, color='k', lw=0.8)
plt.annotate(f'固定端 {f_ind(Lv)/1e3:.3f} kN·m', (0, f_ind(Lv) / 1e3),
             textcoords='offset points', xytext=(10, -18))
plt.annotate(f'最大正 {9*wv*Lv**2/128/1e3:.3f} kN·m',
             (Lv - 3 * Lv / 8, 9 * wv * Lv**2 / 128 / 1e3),
             textcoords='offset points', xytext=(-30, 12))
plt.xlabel('固定端Aからの距離 [m]'); plt.ylabel('曲げモーメント [kN·m]')
plt.title('不静定化による曲げモーメントの再分配')
plt.legend(); plt.grid(alpha=0.3); plt.tight_layout(); plt.show()

図から3点が読み取れます。第一に、単純支持では中央に $wL^2/8 = 5.625$ kN·m の正モーメントが1つあるだけですが、不静定化すると固定端に $-5.625$ kN·m の負モーメントが現れ、正モーメントは $3.164$ kN·m（$9wL^2/128$）まで下がります。第二に、両者の最大絶対値はどちらも 5.625 kN·m で同じです。断面設計上の要求は変わりません。第三に、モーメントがゼロになる変曲点が固定端Aから $L/4$（＝自由端Bから $3L/4$、この例では0.75 m）に生じます。ここは実質的にヒンジのように振る舞う点で、配筋やジョイント配置を考えるうえで重要な位置です。

たわみも比べておきましょう。

# 重ね合わせによるたわみ（xは固定端Aから）
d_ind = wv * x**2 * (Lv - x) * (3 * Lv - 2 * x) / (48 * EIv)
f_di = sp.lambdify(x, d_ind, 'numpy')
d_ss = wv * x * (Lv**3 - 2 * Lv * x**2 + x**3) / (24 * EIv)
f_ds = sp.lambdify(x, d_ss, 'numpy')

i = np.argmax(f_di(xs))
print(f"不静定 最大たわみ {1e3*f_di(xs)[i]:.4f} mm at x={xs[i]:.4f} m (={xs[i]/Lv:.4f}L)")
print(f"単純支持 中央たわみ {1e3*f_ds(Lv/2):.4f} mm")
print(f"比 = {f_di(xs)[i]/f_ds(Lv/2):.4f}")

plt.figure(figsize=(9, 4))
plt.plot(xs, 1e3 * f_ds(xs), lw=2, label='単純支持（最大 1.318 mm）')
plt.plot(xs, 1e3 * f_di(xs), lw=2, label='一端固定（最大 0.548 mm）')
plt.gca().invert_yaxis(); plt.xlabel('固定端Aからの距離 [m]')
plt.ylabel('たわみ（下向き正）[mm]'); plt.title('不静定化によるたわみの低減')
plt.legend(); plt.grid(alpha=0.3); plt.tight_layout(); plt.show()

不静定はりの最大たわみは 0.5484 mm、位置は301点グリッド上で $x=1.740$ m（$0.5800L$）と出ました。解析的な最大位置 $x=\frac{15-\sqrt{33}}{16}L = 0.5785L$（1.7354 m）とグリッド刻み（0.01 m）の範囲で一致しています。単純支持の 1.3184 mm に対して41.6%まで減りました。手計算で得た $\delta_{\max}\approx wL^4/185EI$ に数値を入れると 0.5473 mm で、約0.2%のずれで一致します（185は $184.6$ を丸めたハンドブック値なので、この程度のずれが出ます）。たわみ曲線の形も非対称になり、最大点が固定端から遠い側へ寄っているのが確認できます。固定端が回転を拘束することで、そちら側のたわみが強く抑えられているためです。

Pythonでの実装7：はり要素FEMとの照合

最後に、単位荷重法の結果をはり要素FEM（Hermite補間のEuler–Bernoulli要素）と突き合わせます。

import numpy as np

Lv, wv, EIv = 3.0, 5e3, 4e6

def beam_fem(ne):
    """単純支持ばり＋等分布荷重を ne 要素で解き、節点たわみを返す"""
    n = ne + 1; le = Lv / ne
    K = np.zeros((2 * n, 2 * n)); F = np.zeros(2 * n)
    ke = EIv / le**3 * np.array([
        [12, 6*le, -12, 6*le], [6*le, 4*le**2, -6*le, 2*le**2],
        [-12, -6*le, 12, -6*le], [6*le, 2*le**2, -6*le, 4*le**2]])
    fe = np.array([wv*le/2, wv*le**2/12, wv*le/2, -wv*le**2/12])   # 等価節点荷重
    for e in range(ne):
        d = [2*e, 2*e+1, 2*e+2, 2*e+3]
        K[np.ix_(d, d)] += ke; F[d] += fe
    fx = [0, 2*n - 2]                                   # 両端のたわみを拘束
    fr = [i for i in range(2 * n) if i not in fx]
    u = np.zeros(2 * n); u[fr] = np.linalg.solve(K[np.ix_(fr, fr)], F[fr])
    return u[::2]

exact = lambda s: wv * s * (Lv**3 - 2*Lv*s**2 + s**3) / (24 * EIv)  # 単位荷重法の解
for ne in (2, 4, 8, 16):
    v = beam_fem(ne); s = np.linspace(0, Lv, ne + 1)
    print(f"要素数 {ne:2d}: 中央 {1e3*v[ne//2]:.9f} mm, "
          f"節点最大誤差 {1e3*np.max(np.abs(v - exact(s))):.3e} mm")

結果は次の通りです。

要素数を増やしても答えが変わらず、わずか2要素で厳密解に一致しています。誤差はすべて丸め誤差のレベルで、しかも要素数が増えるほど（連立方程式が大きくなるぶん）わずかに増えるという、収束とは逆の挙動すら見えます。

これは偶然ではありません。Hermite梁要素の形状関数は3次多項式で、同次方程式 $EIv””=0$ の解空間（3次以下の多項式）を完全に含みます。さらに等分布荷重を等価節点荷重（$wl/2$ と $\pm wl^2/12$）に変換しておけば、節点における変位は分割数によらず厳密になることが理論的に保証されています（要素内部の分布は3次近似なので厳密ではありません）。つまり「単位荷重法の解」と「はり要素FEMの節点値」は、そもそも同じ厳密解を別ルートで表現しているだけです。

この一致は実務でも意味があります。市販FEMを使うとき、単位荷重法で手計算した1点の値と照合すれば、モデル化（境界条件・単位系・断面定数）のミスを短時間で洗い出せます。単位荷重法は、FEM時代においてもなお「電卓で殴れる検算手段」として現役なのです。

よくあるつまずき

最後に、単位荷重法でつまずきやすいポイントをQ&A形式でまとめます。

Q. $M$ と $m$ の符号規約が食い違うとどうなりますか？ 片方だけ逆にすると、積 $Mm$ の符号が全区間で反転し、答えの符号が逆になります。区間ごとに規約がばらつくと、部分的に符号が狂って値そのものが間違います。実用上の対策は単純で、$M$ 図と $m$ 図を同じ紙の上下に並べて描き、同じ側を「正」と決めること。プログラムでは、実状態と仮想状態を同じ関数で解いてしまうのが最も安全です（実装4のトラスがまさにそれです）。

Q. 求めたい方向と逆向きに単位荷重を置いてしまったら？ 答えの符号が反転するだけで、絶対値は正しく出ます。$\delta<0$ が出たら「置いた向きと逆に動いた」と読めばよく、置き直して計算し直す必要はありません。むしろ、変位の向きが直感でわからないときは適当な向きに置いて符号で判定するのが定石です。

Q. 単位荷重の「1」に単位は要りますか？ 式(2)の左辺は $1\cdot\delta$ で、右辺は $\int Mm/EI\,dx$。次元を追うと、$m$ は「単位荷重 × 長さ」の次元なので、右辺全体は「力 × 長さ」＝仕事の次元です。つまり左辺の1も力の次元を持っています。手計算では1 kNや1 Nと決めておき、最後に割り戻すのが安全です。プログラムでは実装4のように 1.0 N を明示的に載せておけば、単位系の混乱が起きません。

Q. カスティリアーノの定理とはどう違うのですか？ 実は同じものです。カスティリアーノの第2定理 $\delta_i = \partial U/\partial P_i$ で、ひずみエネルギー $U=\int M^2/(2EI)dx$ を微分すると $\int (M/EI)(\partial M/\partial P_i)dx$ となり、この $\partial M/\partial P_i$ がまさに $m$ に等しくなります（$P_i$ を1だけ増やしたときのモーメント増分＝単位荷重によるモーメント）。荷重の載っていない点の変位を求めたいときにカスティリアーノでは「ダミー荷重を置いて最後に0にする」という操作が要りますが、単位荷重法ではそれが最初から組み込まれています。手続きが素直なぶん、単位荷重法のほうが使いやすいというのが実感です。詳しくはカスティリアーノの定理の記事を参照してください。

Q. 曲がったはりやラーメンにも使えますか？ 使えます。式(2)の $dx$ を部材軸に沿った弧長 $ds$ に置き換え、部材ごとに局所座標で積分して足すだけです。ラーメンでは柱と梁で別々に $M$、$m$ を立てて合計します。曲率半径が断面せいに比べて十分大きくない曲りばりでは、曲げ理論そのものを補正する必要がありますが、単位荷重法の枠組み自体は変わりません。

Q. 不静定次数が2以上のときは？ 冗長力を $X_1, X_2, \dots$ と複数取り、整合条件を連立させます。

$$ \sum_{j} \delta_{ij} X_j + \delta_{i0} = 0 \quad (i=1,2,\dots,n) $$

係数行列 $[\delta_{ij}]$ は柔性行列で、相反定理から対称です。しかも $\delta_{ii}=\int m_i^2/EI\,dx>0$ で正定値なので、常に一意に解けます。この連立方程式は構造力学ではMaxwell–Mohrの正準方程式と呼ばれ、応力法（力法）の中核をなします。基本系の選び方は自由ですが、非対角項 $\delta_{ij}$ が小さくなるように選ぶと数値的に安定します。

まとめ

本記事では、単位荷重法を仮想仕事の原理から導出し、Pythonで検証しました。

• 核心の直感：知りたい点に大きさ1の仮想荷重を置くと、外部仮想仕事が $1\cdot\delta$ のただ1項になり、その点の変位だけが「フィルタ」で取り出される。ベクトルの成分を単位ベクトルとの内積で抜き出すのと同じ発想
• Maxwell–Mohrの式：仮想仕事の原理 $\sum\bar{P}u=\int\bar{\sigma}\varepsilon\,dV$ に、はりの応力分布 $\sigma=My/I$ を代入して断面積分すると $\int_A y^2dA=I$ が現れ、$1\cdot\delta=\int Mm/EI\,dx$ が得られる
• 一般形：曲げ・軸力・せん断・ねじりの4項はすべて「実の内力 × 仮想の内力 ÷ 剛性」の形。細長いはりではせん断／曲げ比が $0.78(h/L)^2$ でスケールし、$L/h=10$ で0.78%、$L/h=20$ で0.195%と無視できる
• トラス：軸力一定なので積分が総和になり $\delta=\sum N_in_iL_i/(E_iA_i)$。7部材のワーレントラスで単位荷重法と直接剛性法が16桁一致（0.26259 mm）。各項の大小がそのまま「どの部材を太くすべきか」の感度になる
• 相反定理：$\delta_{ij}=\int m_im_j/EI\,dx$ が $i,j$ 対称なので $\delta_{ij}=\delta_{ji}$。トラスで数値検証し、10.954 µm が両方向で完全一致
• 1次不静定：整合条件 $\delta_{10}+X\delta_{11}=0$ の1本で解ける。propped cantileverで $X=3wL/8$、固定端モーメント $-wL^2/8$、最大正モーメントが $wL^2/8 \to 9wL^2/128$（56%）に、たわみが41.6%に低減
• FEMとの関係：Hermite梁要素の節点値は分割数によらず厳密（2要素で1.318359375 mm）。単位荷重法は市販FEMの検算手段として現役

単位荷重法は「古典的な手計算法」に見えますが、その正体は柔性行列の成分を1つずつ取り出す積分表現です。この視点を持つと、応力法（力法）、影響線、感度解析、FEMの検証といった一見バラバラなトピックが、すべて同じ式 $\delta_{ij}=\int m_im_j/EI\,dx$ の別の顔として見えてきます。

次のステップとして、以下の記事も参考にしてください。

• 仮想仕事の原理と単位荷重法による構造変位の理論と実装 — エネルギー法全体の中での位置づけ
• カスティリアーノの定理とエネルギー法を基礎から徹底解説 — 単位荷重法と等価なもう1つのアプローチ
• 梁のたわみ（弾性曲線方程式） — 微分方程式による正攻法との比較
• 内力図（SFD・BMD）の描き方 — $M$ 図・$m$ 図を正確に描くための基礎
• はりの境界条件 — 基本系の選び方を考えるうえでの前提
• オイラーの座屈荷重 — 同じエネルギー的視点から見る安定問題