機械(工学の)学習に必要な数学をやり直してみた話

ちまたでは、機械学習がブームのようです。
が、、まったく時代についていけていません。

機械学習は、数学が大事と言われています。

そこで、機械(工学の)学習に必要な数学をやり直してみようか
という、まさかの左斜め下方向からのアプローチで、すこしタイムスリップしてみました。

※注意:　流行りの機械学習の話は一切出てきません

1.機械工学とは

機械工学(Mechanical Engineering)とは、機械の設計／製作／運用の全てを対象とする工学の分野です。

次の4つの力学(通称:よんりき)を基礎として、あらゆる産業の基幹を支える重工業をはじめ、航空機／鉄道／自動車／船舶／ロボット／ロケット／人工衛星などありとあらゆる分野に応用されている工学分野です。

1. 機械力学
2. 熱力学
3. 流体力学
4. 材料力学

機械工学が取り扱う学問は、広範囲にわたるので、このうち熱力学と流体力学を取り上げます。

■熱力学とは

熱力学とは、熱や物質の輸送現象やそれに伴う力学的な仕事を解明する学問です。熱力学には次の基本法則があります。

1. 第0法則（熱平衡）
2. 第1法則（エネルギー保存則）
3. 第2法則（エントロピー増大則）
4. 第3法則（ネルンストの定理）

特に、エネルギーを消費することなく継続的に仕事をする永久機関は存在しない！という第1法則（エネルギー保存則）は、有名です。また、第一種永久機関と呼ばれる「外部から何も受け取ることなく、仕事を外部に取り出すことができる機関」は、ある機関が仕事をするためには外部から熱を受け取るか、外部から仕事をなされるかが必要なので、第一種永久機関は何もエネルギー源の無いところからエネルギーを発生させているので、保存則に反している、と言えます。

熱力学ではエントロピーと呼ばれる、運動状態の混沌性や不規則性の程度を表す重要な状態量があります。これは一般では系の乱雑さを表す量と言われていて、自然界では常にエントロピーが小→大に方向に進むというのが、エントロピー増大則です。つまり、閉じた系では、必ず秩序から無秩序へ向かう！というわけです。私は、方向性のある現象は、一方向には進むものの、逆方向には戻らないものとざっくり理解しています。熱力学におけるエントロピーは、この不可逆性の度合いを、数値にしたものです。

■流体力学とは

流体力学とは、水や空気などの流体(Fluid)が流れる現象を、力学的に解析する学問です。流体に対してさまざまな力が作用することで流れ場がつくられます。流体にかかる力は、圧力／重力／慣性力／せん断応力などがあります。これらの力の収支を取ると状態方程式が導出されます。質量の収支式は連続の式、運動量の収支式はナビエ・ストークスの運動方程式と呼ばれています。

流体には、層流と乱流があり、一般的にはレイノルズ数という流体の慣性力と粘性力の比を表す無次元数がありますこの、がある一定の数値より大きくなると層流から乱流に遷移します。蛇口から出る水の量が少ないときは、規則正しい流れができますがこれが層流です。水量をだんだん増やすと、あるところで不規則な流れに変わります。これが乱流(Turbulence)です。

乱流は、流速／圧力／温度などの物理量が時間的にだけでなく、空間的にも不規則に変動します。そして、大小さまざまなスケールの渦運動が発生します。渦運動には、流体でありながら剛体と同じように、速度が渦中心からの半径に比例して変化する強制渦領域(Forced vortex)と、速度が渦中心からの半径に反比例して変化する自由渦領域(Free vortex)があります。流れ場に渦が多数存在すると、フローパターンはものすごく複雑になります。よって、乱流はいまだ解明できてないことが数多くあります。

だから、おもしろい。というわけでCFDの登場です。

2.数値流体力学(CFD)とは

流体力学や熱力学で導出される偏微分方程式の多くは、解析解を求めることができないといわれています。そこで、境界条件(Boundary Condition)と初期条件(Initial Condition)を決め、支配方程式の近似解を求めるため、計算機でシミュレーションして、実際の工業の分野に応用させていきます。熱エネルギーの移動も含める場合は、数値熱流体力学と呼びます。

数値流体力学は、流体力学の支配方程式の理解に加えて、非線形偏微分方程式を計算機上で解析できるようにする離散化という処理やプログラミング言語やアルゴリズムの知識が必要になってきます。
あと、仮定した理論や計算結果の正当性を検証するため、実験を行い結果を評価したりもします。

CFDは、航空や自動車などの空力制御などや、重工業やプラントなどでの管内流れ(非定常流や粘性流体)の解析などのエネルギー分野でも多く使われています。

3.数値流体力学(CFD)のために必要な基礎知識

まず、数値流体解析に必要な知識を大まかにまとめると、次のようになります。

理論

■ 流体力学

流体力学の基本は、やはり力学です。まず、力学を学ぶうえで、質量／長さ／時間などのSI単位や流速／運動量／エネルギーなどの物理量のちゃんとした理解が必要です。高校までの数学や物理は、教科書に公式があって、それを丸暗記してパターン化すれは、テストでよい点をとることはできます。でも、力学の概念的なところを自分の頭で考えることが難しい＆苦手な人は、CFDは難しいと思います。

■ 解析学

CFDでは、ある程度の数学の知識が必要です。特に、解析学です。古典ニュートン力学の基本である微分積分学や複素解析などは基礎的なところを学んでおかないと、CFDで行う計算の意味がまったく理解できないと思います。ただし、理学系の理論数学とは異なり、数値解析の世界では「①ある事象を数式で表現し、②それを解くためのテクニック」が必要です。

たとえば、、、、

有限差分法は、微分方程式を解くため、微分を有限差分近似で置き換えて得られる差分方程式で近似する数値解法で、FDMと呼ばれてます。

有限要素法は、解析的に解くことが難しい微分方程式の近似解を数値的に得る方法の一つです。方程式が定義された領域を要素に分割し、補間関数をつかって近似していきます。ナビエ・ストークスの運動方程式など離散化して行列式に分解するときによく使われます。有限要素法には、メッシュ/節点/要素といった考え方が出てきます。

計算機

■ プログラミング

数値解析の主たるタスクは、離散化した方程式を解くことです。そのため、ニュートン法、ガウスの消去法、最小二乗法、テーラー展開などを理解しておく必要があります。あとは、膨大な計算を行うため、計算機の知識のプログラミング言語が必要です。CFDの場合は、昔はFortranが多く使われてましたが、最近はC++などが多いようです。そして数値解析を行う上では、誤差の概念や、解の収束や発散などの考え方も必要です。

また、良い解を求めるためには、境界条件や初期条件をどう設定するかが重要です。CFDの境界条件の設定は、試行錯誤によるところもありますが「流体力学の深い知識」が有効に効いてくると思います。

■ ハードウエア

大量の演算を行いますので、計算機のマシンリソースを有効に使うための知識が必要です。機械学習と異なり、大量データを扱う知識よりは、ハードウエアに近い低レイヤーの知識（プロセッサやネットワークなど）が重要だと思います。

ビジネス

数値流体解析が使われる分野は、重工業・プラント・自動車・航空機などの輸送機器など多岐にわたります。
シミュレーションは、あくまでもシミュレーションでしかないので、実際にどう活用するかは、それぞれのビジネスに依存します。

というわけで、ぼんやりを思い出したところで、熱力学・流体力学の基礎方程式を導出する簡単な計算ドリルを解きます。

計算ドリル１: 熱伝導の基礎方程式

---

まずは、熱力学の基礎である熱伝導方程式を導出していきます。
熱伝導率をとし、高さ、奥行きのコントロールボリューム内の軸方向のヒートフローを考えます。

微小時間の間に、の位置にあるコントロールボリュームに流入する熱量は

$$
g\left( x\right) \Delta t = -k\dfrac {\partial T} {\partial x}\Delta t
$$

で算出すれば、それより離れたコントロールボリューム外に流出する熱量は

$$
g\left( x +\Delta x\right) \Delta t = \left( -k\dfrac {\partial T} {\partial x}+\dfrac {\partial } {\partial x}\left( -k\dfrac {\partial T} {0x}\right) \Delta x + O\left( \Delta x^{2}\right) \right) \Delta t
$$

となります。

単位時間・単位質量あたりの発熱があれば、微小時間にコントロールボリュームないで内部発熱する熱量は、密度をとすると

$$
\left( \rho \Delta x\right) S_{k}\Delta t
$$

となります。この内部発熱量はコントロールボリューム内の温度上昇として蓄積される熱量と流出する熱伝導量と一致するはずなので、比熱をとして式にすると

$$
\begin{align}
\left( \rho \Delta x\right) S_{k}\Delta t &= \left( p\Delta x\right) c\Delta T +(q\left( x+\Delta x\right) \Delta t -q\left( x\right) \Delta t) \\
&=\left( \rho \Delta x\right) c(\dfrac {\partial T} {\partial t}\Delta t + o\left( \Delta t^{2}\right) )+\left( \dfrac {\partial } {\partial x}(-k\dfrac {\partial T} {\partial x}\right) \Delta x+O\left( \Delta x^{2})\right) \Delta t
\end{align}
$$

となります。

ここで、、とすると、次の式となります。

$$
\rho c\dfrac {\partial T} {\partial t}=\dfrac {\partial } {\partial x}\left( k\dfrac {\partial T} {\partial x}\right) +\rho S_{h}
$$

これは、方向のみなので、コントロールボリュームをの3次元として、同様の式を解いていくと、

$$
\rho c\dfrac {\partial T} {\partial t}=\dfrac {\partial } {\partial x}\left( k\dfrac {\partial T} {\partial x}\right)+\dfrac {\partial } {\partial y}\left( k\dfrac {\partial T} {\partial y}\right) +\dfrac {\partial } {\partial z}\left( k\dfrac {\partial T} {\partial z}\right)+\rho S_{h}
$$

となります。

流体は管内を流れることが多いので、デカルト系ではなく円筒座標系で同じように式を解いていくと

$$
\rho c\dfrac {\partial T} {\partial t}= \dfrac {1} {\gamma } \dfrac {\partial } {\partial r}\left( rk\dfrac {\partial T} {\partial r}\right)+\dfrac {\partial } {\partial x}\left( k\dfrac {\partial T} {\partial x}\right) +\rho S_{h}
$$

となり、一般形として次のように表します。

$$
\rho c\dfrac {\partial T} {\partial t}=\nabla \cdot \left( k\nabla T\right) +\rho S_{n} \tag{1}
$$

これが3次元空間における熱伝導の方程式が導出できます。

計算ドリル２:ナビエ・ストークス運動方程式

---

次は流体の運動量を考えます。
ある流体の微小なコントロールボリュームに作用する力を整理します。

密度を、流体の質量分率をとすると、コントロールボリュームの面より単位時間あたりに流入する質量は

$$
\left( \rho \phi u-\rho\Gamma _{\phi }\dfrac {\partial \phi } {ax}\right) \left( \dfrac {\Delta V} {\Delta x}\right)
$$

となります。ここで、離れた面から流出する質量を考えたとき、テイラー展開の第2項までをとり、高次項を無視すると、

$$
\left( \rho \phi u-\rho\Gamma _{\phi }\dfrac {\partial \phi } {\partial x}\right) \left( \dfrac {\Delta V} {\Delta x}\right) +\Delta x\dfrac {\partial } {\partial x}\left( \rho \phi u-pf_{\phi }\dfrac {\partial \phi } {\partial x}\right) \left( \dfrac {\Delta V} {\Delta x}\right)
$$

したがって、単位時間に方向へ流出する質量は

$$
\dfrac {\partial } {\partial x}\left( \rho \phi u-\rho\Gamma _{\phi }\dfrac {\partial \phi } {\partial \phi } {\partial x}\right) {\Delta V}
$$

となります。

これは、方向のみなので、コントロールボリュームをデカルト系の3次元として、同様の式を解いていくと、

$$
\Biggl( \dfrac {\partial } {\partial x}\left( \rho \phi u-\rho \Gamma _{\phi }\dfrac {\partial \phi } {\partial x}\right) +\dfrac {\partial } {\partial y}\left( \rho \phi v-\rho \Gamma _{\phi }\dfrac {\partial \phi } {\partial y}\right)+\dfrac {\partial } {\partial z}\left( \rho \phi w-\rho \Gamma _{\phi }\dfrac {\partial \phi } {\partial z}\right)\Biggl)\Delta V
$$

となります。単位時間・単位質量あたりの質量をとすると、、とすると、次の式となります。

$$
\rho S\phi =\dfrac {\partial \rho \phi } {\partial t}+\dfrac {\partial } {\partial x}\left( \rho \phi _{u}-\rho\Gamma _{\phi }\dfrac {\partial \phi } {\partial x}\right) +\dfrac {\partial } {\partial y}\left( \rho \phi _{v}-\rho\Gamma _{\phi }\dfrac {\partial \phi } {\partial y}\right) +\dfrac {\partial } {\partial z}\left( \rho \phi _{w}-\rho\Gamma _{\phi }\dfrac {\partial \phi } {\partial z}\right)
$$

この式を、一般形で次のように表します。

$$
\rho S\phi =\dfrac {\partial S\phi } {\partial t}+\nabla \cdot (p\phi u-\rho \Gamma _{\phi }\nabla \phi )　\tag{2}
$$

ここで、質量分率、生成率とすると、流体力学における質量保存則を表した連続の式となります。

$$
\dfrac {\partial \rho } {\partial t}+\nabla \cdot \left( \rho u\right) =0 \tag{3}
$$

デカルト系で書くと

$$
\dfrac {\partial \rho } {\partial t}+　\dfrac {\partial \left( \rho u\right) } {\partial x} +　\dfrac {\partial \left( \rho v\right) } {\partial y} +　\dfrac {\partial \left( \rho w\right) } {\partial z} =0
$$

で、円筒座標系で書くと

$$
\dfrac {\partial \rho } {\partial t}+　\dfrac {\partial \left( \rho u\right) } {\partial x} + \dfrac {1} {r}\dfrac {\partial \left( r\rho v\right) } {\partial r} =0
$$

式(2)と式(3)から、次の式が得られます。

$$
\rho \left( \dfrac {\partial \phi } {\partial x}+\left( u\cdot \nabla \right) \phi \right) =\nabla \cdot \left( \rho \Gamma \phi \nabla \phi \right) +\rho S_{\phi } \tag{4}
$$

ここで、は速度成分、は動粘性係数で、密度が変化しないとみなせる非圧縮性流体を考え、かつ体積力が働かないものとすれば、単位質量あたりの方向の運動量の生成率に対応するは

$$
S_{\phi }=-\dfrac {1} {\rho }\dfrac {\partial p} {\partial x} \tag{5}
$$

です。式(5)を式(4)に代入して、運動量を求めると

$$
\rho \left( \dfrac {\partial u} {\partial t}+u\dfrac {\partial u} {\partial x}+v\dfrac {\partial u} {\partial y}+w\dfrac {\partial u} {\partial z}\right) =-\dfrac {\partial p} {\partial x}+\mu \left( \dfrac {\partial ^{2}u} {\partial x^{2}}+\dfrac {\partial ^{2}u} {\partial y^{2}}+\dfrac {\partial ^{2}u} {\partial z^{2}} \right)
$$

$$
\rho \left( \dfrac {\partial v} {\partial t}+u\dfrac {\partial v} {\partial x}+v\dfrac {\partial v} {\partial y}+w\dfrac {\partial v} {\partial z}\right) =-\dfrac {\partial p} {\partial y}+\mu \left( \dfrac {\partial ^{2}v} {\partial x^{2}}+\dfrac {\partial ^{2}v} {\partial y^{2}}+\dfrac {\partial ^{2}v} {\partial z^{2}} \right)
$$

$$
\rho \left( \dfrac {\partial w} {\partial t}+u\dfrac {\partial w} {\partial x}+v\dfrac {\partial w} {\partial y}+w\dfrac {\partial w} {\partial z}\right) =-\dfrac {\partial p} {\partial z}+\mu \left( \dfrac {\partial ^{2}w} {\partial x^{2}}+\dfrac {\partial ^{2}w} {\partial y^{2}}+\dfrac {\partial ^{2}w} {\partial z^{2}} \right)
$$

となります。

これは流体力学で最も重要な支配方程式であるナビエ・ストークスの運動方程式です。

このナビエ・ストークスの運動方程式をベクトル表示で書くと次の式になります。

$$
\rho \biggl(\dfrac {\partial u} {\partial t}+\left( u\cdot \nabla \right) u\biggl) =-\nabla p-\mu \nabla ^{2}u
$$

上記のようにナビエ・ストークスの運動方程式は非線形偏微分方程式で、解析解を求めるのが難しいといわれています。
そこで、初期条件・境界条件をきめ、系の状態の近似解を解いていくのが、数値流体力学です。

4.まとめ

ずいぶん、昔のことなのできれいさっぱり忘れてましたが、ちょっと手を動かしてみれば、ぼんやり思い出しました。

何年たっても、やり直せるもんだなとも思いました。

流れは人間の目に見えません。
けど、流れを数式で表現すれば、流れがどうなっているのがを知ることができます。
ロマンですね。

おわり