$ \DeclareMathOperator{\arccosh}{arccosh} \DeclareMathOperator{\arcsinh}{arcsinh} \DeclareMathOperator{\rank}{rank} \DeclareMathOperator{\rot}{rot} \DeclareMathOperator{\grad}{grad} \DeclareMathOperator{\diver}{div} $

Euler-Bernoulli beam equation (モールの定理)

細長い棒が力を受けた場合にどのように変形するか、という問題は、梁などにも適用できるため、建築工学の分野でよく取り上げられるらしいです。
棒の変形が小さい時に近似的に成り立つ方程式として Euler-Bernoulli beam equation[1] というのがあります。
今回はそれを紹介したいと思います。

[1]日本ではモールの定理という名称を使うことが多いらしい

1beam equation (モールの定理)

1.1状況説明

Euler-Bernoulli beam equation
細長い棒が外部から単位長さあたり w の力[1]を受けている状況を想定します。 ただしその外力による棒の変形は微小であるとします。
上の図で y は変位で R は曲率半径です。
F はその部分より左側の棒の部分にかかっている外力 w の合計値、すなわち、
\begin{equation} F = \int w dx \end{equation}
(1)
です。
$M_\sigma$ はその断面にかかる応力が微小部分に与える力のモーメント(トルク)を表します。

1.2モーメントのつり合い

この微小部分に対して力のモーメントのつり合いの条件を考えていきたいと思います。
微小部分にかかっている力は、応力 $\sigma$ と外力 w のみなのでそれぞれについて考えていきましょう。

1.2.1外力のモーメント

まず外力ですが、これは簡単で、微小部分の長さを $\Delta x$ とすれば、
\begin{align*} M_{\text{ex}} & = - \frac{\Delta x}{2} F(x) - \frac{\Delta x}{2} F(x + \Delta x) \\ & = - \frac{\Delta x}{2} \left(F(x) + F(x + \Delta x)\right) \end{align*}
(2)
ですが、$\Delta x \ll 1$ を考えれば $F(x) \approx F(x + \Delta x)$ ですから、
\begin{equation} M_{\text{ex}} = - F(x) \Delta x \end{equation}
(3)
とできます。

1.2.2応力のモーメント

次に応力です。
伸び率と応力の関係
棒の単位長さあたりの伸び h と応力 $\sigma$ の関係は、ヤング率を E として、
\begin{equation} \sigma = E h \end{equation}
(4)
と表すことができるのでした。
そこで棒の伸びを求めてみましょう。
伸びを求める
棒が曲がる時には棒の内側の部分は縮み、棒の外側の部分は伸びます。 棒の中間付近には伸びも縮みもしない場所があり、それを中立面(neutral surface)といいます。 この中立面から r だけ離れた部分での伸びを求めましょう。
$r = r$ での長さは、
\begin{equation} (R + r) \Delta \theta \end{equation}
(5)
であり、中立面での長さは、
\begin{equation} R \Delta \theta \end{equation}
(6)
なので結局、伸びは、
\begin{equation} (R + r) \Delta \theta - R \Delta \theta = r \Delta \theta \end{equation}
(7)
となります。
伸び率を求める
ところがこれは $\Delta x = R \Delta \theta$ あたりの伸びなので、これを単位長さあたりの伸びに直すと、
\begin{equation} \frac{r \Delta \theta}{R \Delta \theta} = \frac{r}{R} \end{equation}
(8)
となります。
従って応力は、
\begin{equation} \sigma = E \times \frac{r}{R} = \frac{Er}{R} \end{equation}
(9)
となります。
曲率半径について
で。
r はいいとして問題は R です。 このままではイケマセンね。 どうにかしましょう。
曲率半径は一般に、
\begin{equation} R = \left( 1 + \left(\frac{\partial y}{\partial x}\right)^2 \right)^{\frac{3}{2}} \left( \frac{\partial^2 y}{\partial x^2} \right)^{-1} \end{equation}
(10)
と表せます。 ところが、いま、変形は微小だとしていたので棒の傾き $\frac{\partial y}{\partial x}$ は 1 に比べるととても小さいはずです。
このことを用いると曲率半径は、
\begin{equation} R = \left( \frac{\partial^2 y}{\partial x^2} \right)^{-1} \end{equation}
(11)
と書くことができます[2]
応力
以上から応力は、
\begin{equation} \sigma = ErR^{-1} = E r \frac{\partial^2 y}{\partial x^2} \end{equation}
(12)
となります。
応力によるモーメント
さて、応力は分かったのでモーメントは積分するだけです。
\begin{align*} M_\sigma (x) & = \int r \sigma \times l(r) dr \\ & = E \frac{\partial^2 y}{\partial x^2} \int r^2 l(r) dr \end{align*}
(13)
となりますが、最後の式の積分部分は棒の形状にのみ依存する量です。 そこでこれを、
\begin{equation} I := \int r^2 l(r) dr \end{equation}
(14)
とおいてしまいます。 この I のことを断面二次モーメント(second moment of area)といいます。
この断面二次モーメントを用いると、応力のモーメントは結局、
\begin{equation} M_\sigma (x) = E I \frac{\partial^2 y}{\partial x^2} \end{equation}
(15)
となります。

1.2.3式を立てる

さてこれで全部のモーメントが出そろいましたので、あとは式を立てて計算するだけです。 ただし、符号には十分注意してください。 図をよ〜〜〜く見て、きちんと符号があっていることを確認してくださいねw
\begin{equation} M_\text{ex} + M_\sigma (x + \Delta x) - M_\sigma (x) = 0 \end{equation}
(16)
\begin{equation} M_\sigma (x + \Delta x) - M_\sigma (x) = F \Delta x \end{equation}
(17)
\begin{equation} \frac{M_\sigma (x + \Delta x) - M_\sigma (x)}{\Delta x} = F \end{equation}
(18)
ここらへんで $\Delta x \rightarrow 0$ なる極限をとってみると、
\begin{equation} \frac{\partial M_\sigma}{\partial x} = F \end{equation}
(19)
となります。 F は w の積分でしたので、この両辺を x でさらに偏微分して、
\begin{equation} \frac{\partial^2 M_\sigma}{\partial x^2} = w \end{equation}
(20)
\begin{equation} \frac{\partial^2}{\partial x^2} \left( E I \frac{\partial^2 y}{\partial x^2} \right) = w \end{equation}
(21)
となります。
この式を Euler-Bernoulli beam equation といいます。 beam は「梁」という意味なので、日本語にするとしたら「Euler-Bernoulliの梁方程式」とかでしょうか? ちなみに英語圏では基本的にこの呼び方らしいのですが、日本では主にモールの定理(Mohr's theorem)というらしいです。 モールって誰?
ところで上の式は Wikipedia がこうなっていたからこんな形にしましたけど、断面の形が一定で均質な棒を考えれば E も I も定数ですし、棒の変位 y も普通、位置 x にのみ依存するので、その場合にはより簡単に、
\begin{equation} EI\frac{d^4y}{dx^4} = w \end{equation}
(22)
と書くこともできます。
ちなみにこの EI が大きい物体ほど、同じ力を加えても曲げが小さくなるので、その意味で、
\begin{equation} B := EI \end{equation}
(23)
のことを曲げ剛性率(flexural rigidity)などといいます。

[1]単位面積あたりの力を圧力というのに似てますねw…って Wikipedia に書いてあった
[2]結構この近似はいろいろなところで使うらしいですよ。どうでもいいですが

2具体例

さて、基礎方程式が得られたので具体的な状況を考えて、棒の変形は実際にどんな感じになるのかを計算してみましょう。

2.1境界条件

…といいたいところですが、まずは境界条件について。

2.1.1固定端

まず棒が壁に固定されている場合(固定端(fixed end))には、位置は固定ですから、
\begin{equation} y = 0 \end{equation}
(24)
であり、また壁に接している部分では向きは真横なので、
\begin{equation} \frac{\partial y}{\partial x} = 0 \end{equation}
(25)
となります。

2.1.2自由端

自由端(free end)の場合には、この場所の右 or 左には弾性体がないので応力が働きません。 従って、応力のモーメントも 0 ですから、(15)式は 0 になるので、
\begin{equation} \frac{\partial^2 y}{\partial x^2} = 0 \end{equation}
(26)
となります。

2.1.3ピンによる固定

棒の端がピンにより固定されている場合(pin connection)には、まず位置は固定なので、
\begin{equation} y = 0 \end{equation}
(27)
であり、壁から応力をもらうことはできない[1]ので、
\begin{equation} \frac{\partial^2 y}{\partial x^2} = 0 \end{equation}
(28)
となります。

2.1.4一点に荷重がかかっている場合

ある一点[2]荷重(load) N がかかっている場合には(19)式から、
\begin{equation} N = - \frac{\partial}{\partial x} \left(E I \frac{\partial^2 y}{\partial x^2}\right) \end{equation}
(29)
となります。 負号については、(19)式の F の定義はその地点から左側の棒全体に働く力の合計でしたが、上向きに N の力を加えると壁などから N の力を反対向きに受けてつり合うことから、この壁などから下方向にもらう力をあらわしています。
これはもちろん、 EI が一定であるなどの仮定を採用すれば、
\begin{equation} N = - EI\frac{d^3y}{dx^3} \end{equation}
(30)
とすることができます。

2.2重力による変形

一端を固定した棒全体に重力がかかっている場合にどのように曲がるのかを考えてみましょう。 なお、簡単のために棒の断面や材質は一定とします[3]

2.2.1一般解

$x = 0$ で固定端、$x = L$ で自由端とします。 棒の単位長さあたりの質量、すなわち線密度(line density)$\rho$ とすると、$w = \rho g$ ですから、
\begin{equation} EI\frac{d^4y}{dx^4} = \rho g \end{equation}
(31)
\begin{equation} \frac{d^4y}{dx^4} = \frac{\rho g}{EI} \end{equation}
(32)
です。 右辺が定数なので、これは簡単に解けて、
\begin{equation} y = Cx^4 + ax^3 + bx^2 + cx + d \end{equation}
(33)
という形になります。 ただし、
\begin{equation} C := \frac{1}{4!} \frac{\rho g}{EI} \end{equation}
(34)
とします。

2.2.2境界条件

次に境界条件を考えます。
x = 0 (固定端)
$x = 0$ で固定端ですので、
\begin{equation} y(0) = 0 \end{equation}
(35)
\begin{equation} \therefore d = 0 \end{equation}
(36)
および、
\begin{equation} y'(0) = 0 \end{equation}
(37)
\begin{equation} \therefore c = 0 \end{equation}
(38)
がいえます。
x = L (自由端)
まず、棒の右端には荷重はかかっていませんから、
\begin{equation} y'''(L) = 0 \end{equation}
(39)
が成り立ちます。 「えっ!棒のどこにでも重力がかかっていんだから 0 になるわけないじゃないか!」とツッコミが入りそうなので、一言。
$y'''$ というのは(19)式などから F と同じ値です。 F というのは w の積分値、つまり注目している部分より左側の棒に働く力の合計値のことです。 そうするとこの状況の場合、棒全体が重力から下方向にもらう荷重と壁から上方向にもらう荷重は等しくなるはずです。 なので、棒の右端まで荷重を積分すると 0 になり、$F = 0$ すなわち $y'''(L) = 0$ がいえるわけです。
よろしいでしょうか? これを(33)式に代入して計算すると、
\begin{equation} 24CL + 6a = 0 \end{equation}
(40)
\begin{equation} \therefore a = -4CL \end{equation}
(41)
となります。
また $x = L$ で自由端ですので、
\begin{equation} y''(L) = 0 \end{equation}
(42)
\begin{equation} 12CL^2 + 6aL + 2b = 0 \end{equation}
(43)
\begin{equation} 12CL^2 - 24CL^2 + 2b = 0 \end{equation}
(44)
\begin{equation} \therefore b = 6CL^2 \end{equation}
(45)
がいえます。

2.2.3特殊解

これらを総合すると、棒の変形は、
\begin{align*} L & = Cx^4 - 4CLx^3 + 6CL^2x^2 \\ & = C(x^2 - 4Lx + 6L^2) x^2 \end{align*}
(46)
となります。

2.2.4具体例

材料は鉄で、長さは 1m とします。 形は円形で、半径はとりあえず r とおいておきましょう。
手元の理科年表によれば、鉄のヤング率などの物理特性はこんな感じらしいです。
  • ヤング率 $E = 21.14 \times 10^{10}\ \text{Pa}$
  • 密度 $\mu = 7.874 \times 10^3\ \text{kg/m}^3$
  • 重力加速度[4] $g = 9.798\ \text{m/s}^2$
そうするとまず、線密度 $\rho$ は、
\begin{equation} \rho = \mu \times \pi r^2 \end{equation}
(47)
です。 また断面二次モーメント I を頑張って計算すると[5]
\begin{equation} I = \frac{\pi r^4}{4} \end{equation}
(48)
となるので、C は、
\begin{align*} C & = \frac{1}{24} \frac{\rho g}{E I} \\ & = \frac{1}{24} \frac{\pi\mu r^2 \times g}{E \times \pi r^4 / 4} \\ & = \frac{\mu g}{6 E r^2} \end{align*}
(49)
となります。 棒が重力でお辞儀する量はだいたい半径の二乗に反比例するんですねw
で、理科年表で調べた値を代入しましょう。 とりあえず直径を 1cm($r=0.5 \times 10^{-3}\ \text{m}$) とします。 電卓で計算すると、
\begin{equation} C = 2.433 \times 10^{-3} \end{equation}
(50)
となります。 長さは $L = 1$ としていたので、(46)式は、
\begin{equation} y = 2.433 \times 10^{-3} \times (x^2 - 4x + 6) x^2 \end{equation}
(51)
なので、棒の右端の垂れ具合は $x = 1$ として、
\begin{equation} y = 7.299 \times 10^{-3}\ \text{m} = 7.299\ \text{mm} \end{equation}
(52)
となります。 直径 1cm の棒を横にすると 0.7cm 垂れるらしいです。 本当でしょうか? 私は手元に直径 1cm の鉄パイプがあるような特殊な人材ではないので検証できませんが、やってみたよ、という人がいましたらご一報を。

2.3荷重による変形

次に棒の左端が壁に固定されていて、右端に N ニュートンの荷重がかかってる状況を考えましょう。

2.3.1一般解

まず、棒の右端と左端以外には荷重はかかっていませんので、
\begin{equation} y'''' = w = 0 \end{equation}
(53)
です。 ということは、この微分方程式の一般解は、
\begin{equation} y = ax^3 + bx^2 + cx + d \end{equation}
(54)
という形になります。

2.3.2境界条件

つづいて境界条件を考えましょう。
x = 0 (固定端)
固定端なので、
\begin{equation} y(0) = 0 \end{equation}
(55)
\begin{equation} \therefore d = 0 \end{equation}
(56)
および、
\begin{equation} y'(0) = 0 \end{equation}
(57)
\begin{equation} \therefore c = 0 \end{equation}
(58)
が成り立ちます。
x = L (荷重のかかっている場所)
右端には荷重が N ニュートンかかっているので、
\begin{equation} E I y'''(L) = -N \end{equation}
(59)
\begin{equation} 6a \times EI = -N \end{equation}
(60)
\begin{equation} \therefore a = -\frac{N}{6EI} \end{equation}
(61)
であり、右側から応力をもらえないので、自由端と同じく、
\begin{equation} y''(L) = 0 \end{equation}
(62)
\begin{equation} 6aL + 2b = 0 \end{equation}
(63)
\begin{equation} b = -3aL = \frac{NL}{2EI} \end{equation}
(64)
となります。

2.3.3特殊解

以上から、
\begin{align*} y & = -\frac{N}{6EI}x^3 + \frac{NL}{2EI}x^2 \\ & = \frac{N}{6EI} (3L - x) x^2 \end{align*}
(65)
と分かります。

2.3.4具体例

またまた具体例を考えてみましょう。
素材と形状は前と同じとします。
そうそう、計算する前に確認。 I は半径 r の四乗に比例する値だったので、荷重がかかった場合の変形量は半径の四乗に比例するみたいですね。 重力の場合には二乗だったので、こっちの変形の方が太さが重要になるようです。
それでまた電卓で頑張って計算すると、
\begin{equation} y = 1.606 \times 10^{-3} \times N (3 - x) x^2 \end{equation}
(66)
となります。 右端の荷重は $10 N$ ぐらい[6]として右端($x=1$)での値を求めると、
\begin{equation} y = 3.212 \times 10^{-2}\ \text{m} = 32.12\ \text{mm} \end{equation}
(67)
となります。 約 $3\ \text{cm}$ か…。 そんなに曲がるんですかねぇ…。 こちらも検証できる方がおりましたらお願いいたします m(_ _)m

[1]壁と棒との間に隙間ができてしまうので
[2]端とは限りません
[3]EI が定数の場合を考えるよ、という意味。以下同じ
[4]東京大学理学部化学館地下原点室における値
[5]ただの三角函数の積分だし、できるよね?
[6]だいたい 1 kg ぐらい