よく「F1マシンの周りの空気の流れをスーパーコンピュータで計算しました」とか「津波を数値解析によってシミュレーションしました」というような話をニュースなどで聞きませんか？ これらの話は，流体（空気や水などの流れ）の動きをコンピュータ上でシミュレーションするCFD（Computational Fluid Dynamics, 数値流体力学）という手法を用いたことを意味しています． 寺本/岡本研究室では，航空宇宙推進系（ロケットエンジンやジェットエンジン）における流体の動きを研究するため，このCFDという手法を実験と同じくらい使用しています． このCFDというものが一体どのようなものなのかイメージがわかないという方も多いのではないでしょうか． ここでは，このCFDというという手法がどのようなものなのか，簡単に（高校生から大学１～２年生くらいの知識を使って）説明します．
まず，質点や剛体の運動を予測するには「運動方程式」，「（角）運動量保存則」，「エネルギー保存則」などの微分方程式を解く必要がありますね． 簡単な例を挙げてみます． 初速 \(u_0 = 0\) m/sかつ質量 \(m = 5\) kgの物体に，外力 \(F = 10t + 5\) [N]（ただし \( t \) は時間[sec]）を加えたときの \(t = 1\) 秒後の物体の速度を求めてみましょう．（ただし \(u\) 方向のみを考える）
この場合，運動方程式
\[ m \frac{du}{dt} = F \]
の微分方程式を解けばいいのは，物理をかじった人であればわかると思います．すなわち，
\[
\begin{align}
& m\frac{du}{dt} = 10t + 5 \ \Leftrightarrow \ \frac{du}{dt} = \frac{10t + 5}{m} \\
∴ \ & u = \int \frac{10t + 5}{m} dt = \frac{5 t^2 + 5t + C}{m} \ (ただしCは積分定数)
\end{align}
\]
\(u_0 = 0\) m/s, \( m = 5 \) kgであるから，
\[
u = t^2 + t
\]
したがって，\( t = 1 \)秒後の質点の速度は，
\[
u = 1^2 + 1 = 2 \ \mathrm{m/s}
\]
と簡単な微分方程式を解くことで求めることができますね．
流体の運動も同様で，いくつかの基礎方程式で記述することが可能です． その方程式は，ナビエ・ストークス方程式と呼ばれており，その方程式を解くことができれば，流体の運動も上で解いた質点の運動と同様に予測することができることになります． そのナビエ・ストークス方程式は，二次元で以下のようになります．[参考: 藤井孝蔵, 流体力学の数値計算法, 東京大学出版会, 1994]
\[
\frac{\partial Q}{\partial t} + \frac{\partial (E – E_\nu)}{\partial \xi} + \frac{\partial (F-F_\nu)}{\partial \eta} = 0, \\
Q = \begin{bmatrix} \rho \\ \rho u \\ \rho v \\ e \end{bmatrix}, \
E = \frac{1}{J} \begin{bmatrix} \rho U \\ \rho u U + \xi_x p \\ \rho v U + \xi_y p \\ (e+p)U \end{bmatrix}, \
F = \frac{1}{J} \begin{bmatrix} \rho U \\ \rho u U + \eta_x p \\ \rho v U + \eta_y p \\ (e + p) U \end{bmatrix}, \\
E_\nu = \frac{1}{J} \begin{bmatrix} 0 \\ \xi_x \tau_{xx} + \xi_y \tau_{xy} \\ \xi_y \tau_{yx} + \xi_y \tau_{yy} \\ \xi_x \beta_x + \xi_y \beta_y \end{bmatrix}, \
F_\nu = \frac{1}{J} \begin{bmatrix} 0 \\ \eta_x \tau_{xx} + \eta_y \tau_{xy} \\ \eta_y \tau_{yx} + \eta_y \tau_{yy} \\ \eta_x \beta_x + \eta_y \beta_y \end{bmatrix}, \\
\tau_{xx} = \frac{2}{3} \mu (2u_x – v_y), \ \tau_{yy} = \frac{2}{3} \mu (2v_y – u_x), \ \tau_{xy} = \tau_{yx} = \mu (u_y + v_x), \\
\beta_x = \tau_{xx} u + \tau_{xy}v + \kappa T_x, \ \beta_y = \tau_{yx} u + \tau_{yy} v + \kappa T_y
\]
それでは，この方程式を質点の運動方程式と同様に解いてみてください（一般解を導いてください）．
たぶん無理だったと思います．「え，解けたけど」っていう人は今すぐ学会に発表する準備をしましょう． ほぼ間違いなく，フィールズ賞がもらえます． 実はこのナビエ・ストークス方程式，数学者でさえまだ誰も解けていない方程式なのです． 解けない理由は色々あるみたいですが，基本的には式自体が非線形（複雑な形）なので，積分が難しい or できないのが原因のようです． じゃあどうやって流体の動きを予測するのかというのを，先ほどの質点の運動方程式の例を使って説明したいと思います．
先ほどの運動方程式から，加速度の時間変化は，
\[
\frac{du}{dt} = \frac{10t + 5}{m} = 2t + 1
\]
と求めることは，積分しなくても導くことができます． それでは，この加速度の式を使って \(t = 1\) 秒後の速度を予測してみましょう． まずは，\(t = 0\) 秒の加速度が \(1\) 秒間続くものだと仮定すると，
\[
\left. \frac{du}{dt} \right|_{t=0} = 1 \ \mathrm{m/s^2}
\]
であるので，\(t = 1\) 秒後の速度は，
\[
u = 1 \ \mathrm{m/s^2} \times 1 \ \mathrm{s} = 1 \ \mathrm{m/s}
\]
と予測できます． しかしこれでは，理論解の \(2\) m/s とかなり離れていますね． そこで今度は，\(t = 0\) と \(t = 0.5\) で加速度を評価して，それぞれの加速度が \(0.5\) 秒続くものだとしてみましょう．
\[
\left. \frac{du}{dt} \right|_{t=0} = 1 \ \mathrm{m/s^2}, \ \left. \frac{du}{dt} \right|_{t=0.5} = 2 \ \mathrm{m/s^2}
\]
であるので，予測される速度は，
\[
u = 1 \ \mathrm{m/s^2} \times 0.5 \ \mathrm{s} + 2 \ \mathrm{m/s^2} \times 0.5 \ \mathrm{s} = 1.5 \ \mathrm{m/s}
\]
となり，先ほどよりも理論解に近づきましたね． 同様にして \(0.1\) 秒ごとに加速度を評価してやると，
\[
u = 1.9 \ \mathrm{m/s}
\]
と予測することができ，かなり理論解に近づきました． さらに細かく \(0.01\) 秒ごとに計算してみると，
\[
u = 1.99 \ \mathrm{m/s}
\]
となります．
このように解析的に積分をしなくても，時間を細かく区切って計算することで，ある程度正確な解を求めることができます． CFDも，同様にしてナビエ・ストークス方程式を時間的に細かく，さらに空間的にも細かく分けて計算することで，この方程式の近似解を求める手法です．
CFDにより，通常では見ることが困難な金属配管内部の流体の流速や圧力分布をCGによって見ることができたり，シミュレーションなので風洞試験などの実験と比べるとかなり安価に結果を得ることができます． しかし，空間を細かく区切って計算するので，エンジン全体のような複雑な流れ場を計算するには，非常に多くの計算コストを必要とするため，現状のCFD技術では再現できないといったデメリットもあります． そこで，寺本/岡本研究室ではこのCFDを使ってどのような現象が起こっているのかを調べることはもちろん，CFDの技術そのものを高めるための手法に関する研究も盛んに行なっています． CFDの結果から得られるCGは綺麗で格好いいものが多いので，興味を持たれた方は是非一度，寺本/岡本研究室を見学に来て下さい．
