ナビエの方程式の導出｜変位と応力の微分方程式とは？

これまでの材料力学の問題では変形を伸び、曲げ、ねじりに分類してそれぞれに対して変形量を個別に解いて行きました。

しかし、一般の変形はこれらの変形が組み合わさったものであり、きっちりと分類して計算できるものではありません。

今回は材料の変形について一般的に成立する方程式、ナビエの方程式を導出していきます。

※ ナビエ・ストークス方程式とは異なります。

平衡方程式の導出

材料力学で扱う物体は伸びたり縮んだりします。

物理学の専門用語ではこのような物体を弾性体と呼びます。

力学や剛体の力学で扱う物体は変形しないことが前提となりますが、材料力学では弾性体を扱い、物体の変形を認めることが特徴になります。

材料力学のもう一つのポイントは、材料力学が静力学に分類されることです。（→静力学・動力学とは？）

したがって、材料力学で扱う弾性体は静止状態にあることが前提になります。

ということで、まずは弾性体の釣り合いについて考えていきましょう。

図のような微小な立方体を考えます。

この立法体の一辺の長さをそれぞれ$\diff x, \diff y, \diff z$とします。

さて、立方体の各面には垂直応力やせん断応力が図のように働いているとします。

ここで、応力の表し方について次のように約束します。

例えば、$x$軸に対し垂直な面に対し$x$方向に働く力は$\sigma_{xx}$と表すことになります。

また、二つの添え字が同じ応力は垂直応力であり、それ以外の応力はせん断応力であることが図より分かります。（→応力・垂直応力・せん断応力とは？）

準備ができたので、立方体の釣り合いについて考えましょう。

今、この微小要素は静止しているため、合力とモーメントの和がそれぞれ$0$となっているはずです。

力の釣り合い

まず、合力の釣り合いから考えましょう。

$x$方向に関する釣り合いを考えると、次のような式が立てられます。

ただし、$F_x$は体積力（重力・電磁気力・浮力など）を表します。

これを整理すると次のようになります。

その他の方向に関しても同様に整理すると、釣り合いの方程式は次のように整理できます。

モーメントの釣り合い

次に、重心を通る$x$軸回りのモーメントの釣り合いについて考えます。（→モーメントとは？）

ただし、体積力は重心に対して対称に分布しているため、体積力によるモーメントは考えなくても良いことに注意してください。

具体的に計算すると次のようになります。

高次の微小項に関する量を無視すると、次のような関係が成立することが分かります。

\begin{split}
\sigma_{yz} = \sigma_{zy}
\end{split}

$y$軸、$z$軸回りのモーメントの釣り合いについても考えると、以下の関係が成り立つことが分かります。

\begin{split}
\sigma_{xy} = \sigma_{yx} \EE
\sigma_{zx} = \sigma_{xz} \EE
\end{split}

平衡方程式と未知数の数

さて、最初に計算した釣り合いの式のことを釣合い方程式、または平衡方程式と呼びます。

平衡方程式には応力に関して9個の未知数が登場します。

しかしながら、式は3つしかないため、方程式を解くことはできません。

そのため、変数の数を減らすことを考えなければなりません。

ここで、モーメントの釣り合いの計算結果を利用します。

モーメントの釣り合い計算の結果より、せん断力に関する変数の中で3つは従属変数であることが分かります。

したがって、未知数（独立変数）は9個から6個に減るのです。

平衡方程式を行列を用いて表記すると、次のようになります。

ここで突然、$\sigma_{ij}$なる表記が登場しましたが、$i, j$は$1,2,3$の数を取る変数とします。

そして$1,2,3$の順で$x,y,z$の対応関係を持つと約束します。

これより、$\sigma_{11} = \sigma_{xx}, \sigma_{23} = \sigma_{yz}$という関係があると言えます。

複数の変数が登場すると、表記が煩雑になるためこのように簡略化して表記する場面が増えます。

これを発展させると、アインシュタインの総和規約としてまとめられます。

変位とひずみの関係式

次に変位とひずみの関係式について導出しましょう。

三次元空間内で図のように変形により$\overrightarrow{AB}$がベクトル$\overrightarrow{A’B’}$に移動したとします。

さて、点$\RM{A}$から$\RM{A’}$までのベクトルを$\vec{u}$とし、点$\RM{B}$から$\RM{B’}$までのベクトルを$\vec{u}+\diff \vec{u}$とすると、以下の関係が成立します。

\begin{eqnarray}
\overrightarrow{A’B’} = \overrightarrow{AB}+\diff \vec{u}
\end{eqnarray}

ここで、各ベクトルの成分が$\overrightarrow{AB} = (\diff x, \diff y, \diff z)$、$\overrightarrow{A’B’}=(\diff x’, \diff y’, \diff z’)$、$\vec{u}=(u_x, u_y, u_z)$、$\diff \vec{u} = (\diff u_x, \diff u_y, \diff u_z)$であるとします。

これより、$\overrightarrow{A’B’}$の成分は次のように表せます。

$$
\left\{
\begin{split}
\diff x’ &= \diff x + \diff u_x \EE
\diff y’ &= \diff y + \diff u_y \EE
\diff z’ &= \diff z + \diff u_z \EE
\end{split}
\right.
$$

また、テイラー展開を利用することで、$\vec{u}$と$\vec{u}+\diff \vec{u}$の間に次のような関係を導けます。（一次までの展開で打ち切っています）

以上より、変形前後の長さの変化を次のように計算できます。

これに先程の関係式を代入して整理すると、以下のようにまとめることができます。

このとき、次のようにひずみ$\varepsilon_{ij}$を定義します。

ここで、アインシュタインの総和規約と呼ばれる表記法を導入します。

例えば、$a_i$は$a_1, a_2, a_3$を代表して表し、$a_{ii}$は$a_{11}+a_{22}+a_{33}$の和を表すことになります。

アインシュタインの総和規約を用いて先述のひずみの定義式を整理すると、以下のようにできます。

※ $i, j, m$は$1,2,3$を動き、$\DL{\ff{\del u_m}{\del x_i}\ff{\del u_m}{\del x_j}}$の項に関して和を取ります。

なお、$x_1 = x, x_2 = y, x_3 = z$とします。

以下、特に断りが無い場合はアインシュタインの総和規約を用いていきます。

金属材料では変位は微小であるため、第三項は無視出来て、したがって以下のようにできます。

ここで、$\varepsilon_{ii}$を垂直ひずみ、$\gamma_{ij}$をせん断ひずみと呼ぶことにします。（→良く分かる垂直ひずみとせん断ひずみ）

一般化フックの法則

より一般的な場合でも成立するフックの法則について導いてみましょう。

ただし、材料のどの方向でも同じ弾性率を持つ均一な弾性体に対し、成立するフックの法則を考えます。

なお、このような弾性体を等方弾性体と呼びます。

弾性体に対して図のような$\sigma_{yy}$の応力が働いているとき、物体は$y$方向の変形に加え、実は$x$軸方向に関する変形も起きるのです。

このとき、$x$軸方向に生じるひずみは$-\DL{\nu\ff{\sigma_{yy}}{E}}$と表せるのです。

同様に$\sigma_{zz}$により$x$軸方向に生じるひずみも、$-\DL{\nu\ff{\sigma_{zz}}{E}}E$となります。

なお、$\nu$（ニュー）はポアソン比と呼ばれる定数を表します。

以上より、各方向に働く応力によって$x$軸方向に生じるひずみは次のように表せるのです。

\begin{eqnarray}
\varepsilon_{xx} &=& \ff{1}{E} \big\{ \sigma_{xx} -\nu(\sigma_{yy} + \sigma_{zz}) \big\} \EE
\end{eqnarray}

$y$軸、$z$軸方向のひずみも同様に計算できて、垂直ひずみと垂直応力に関するフックの法則は次のように表せます。

物体に働く応力が判明している場合、これらの式から垂直ひずみを計算できます。

また、横弾性係数$G$とせん断ひずみ$\gamma$の間には、せん断応力を$\sigma_{xy}$とすると、以下の関係が成立します。（→横弾性係数とは？）

\begin{eqnarray}
\sigma_{xy} &=& G\gamma \EE
&=& G\left( 2\varepsilon_{xy} \right) \EE
\therefore\,\,\varepsilon_{xy} &=& \ff{1}{2G}\sigma_{xy}
\end{eqnarray}

$E$と$G$の間には、$E=2(1+\nu)G$の関係があるため、フックの法則は次のようにまとめて表すことができます。

\begin{eqnarray}
\varepsilon_{ij} = \ff{1}{2G}\left( \sigma_{ij}-\ff{\nu}{1+\nu}\sigma_{kk}\delta_{ij} \right)\qquad (i,j = 1,2,3) \tag{2}
\end{eqnarray}

この式が、フックの法則を一般化して表現するためのファーストステップになります。

ただし、$\delta_{ij}$はクロネッカーのデルタを表します。

ラメの定数

さて、$\sigma_{xx}$に関して式(2)を計算すると、次のようにできます。

\begin{eqnarray}
\varepsilon_{xx} &=& \ff{1}{2G}\sigma_{xx} \,- \ff{1}{2G}\ff{\nu}{1+\nu}\sigma_{ii} \EE
\therefore\,\, \sigma_{xx} &=& 2G\varepsilon_{xx} + \ff{\nu}{1+\nu}\sigma_{ii}
\end{eqnarray}

ところで、式(1)の辺々を足すと、

とできます。

ここで、$\varepsilon_{v}$は各方向の垂直ひずみの和を取ったものであり、体積ひずみと呼ばれます。

これより、$\sigma_{ii} = \DL{\ff{E}{1-2\nu}\varepsilon_{v}}$であることが分かるので、これを先程の式に代入すると、

\begin{eqnarray}
\sigma_{xx} &=& 2G\varepsilon_{xx} + \ff{E}{1-2\nu}\ff{\nu}{1+\nu}\varepsilon_{v}
\end{eqnarray}

と整理できます。

$G=\DL{\ff{E}{2(1+\nu)}}$の関係式を利用して、変形を行うと$G$を用いて次のようにきれいに整理できるのです。

\begin{eqnarray}
\sigma_{xx} &=& 2G\varepsilon_{xx} + \ff{2G\nu}{1-2\nu}\varepsilon_{v}
\end{eqnarray}

ここで、ラメの定数なる定数を導入します。

ラメの定数を利用すると、今まで求めた垂直応力とせん断応力の関係式を次のようにすっきりと表すことができるのです。

ナビエの方程式

第２節で求めた応力とひずみの関係式を、先程導出した一般化フックの法則を代入すると次のようにできます。

いよいよクライマックスです。

上の関係式を平衡方程式に代入します。

$x$軸に関して計算を行うと、

となります。

ここで、ラプラシアンを用いて上式を整理すると以下のように表せます。（→ラプラシアンとは？）

$y$軸、$z$軸に関しても同様に計算することができ、アインシュタインの総和規約を用いてこれらをまとめて表すと次のようになります。

\begin{split}
\mu\D u_i + (\lambda + \mu)u_{j, ji} + F_i = 0
\end{split}

とできるのです。

この方程式はナビエの方程式と呼ばれ、力と変位の関係を規定する方程式です。

未知変数の数と、導出した方程式の数を考えましょう。

未知変数は、3個の変位、6個のひずみ、6個の応力があるため、計15個あります。

一方、方程式については３個の平衡方程式、６個のひずみと変位の関係式、そして６個の応力とひずみの関係式（一般化フックの法則）があるので、計15個あります。

これより、未知変数の数と方程式の数が一致するため、完全に変数を決定でき、解くことができるのです。

つまり、境界条件が与えられればナビエの方程式は完全に解くことができるのです。

とは言え、真正面から解くのは大変なので手作業で解くよりも、コンピュータによる数値計算（シミュレーション）で利用されます。