アイスヘルの数値計算 ～センチュリースープを求めて～

はじめに

ジャンプが誇る本格的グルメ漫画「トリコ」に登場するセンチュリースープというメニューをご存じでしょうか。

かつて大昔のグルメ家たちが己のフルコースの食材を保存するために用いたアイスヘルという氷の大陸(美食の冷蔵庫)があります。アイスヘルで、その食材たちを氷漬けにして保存している氷山を「グルメショーウィンドー」といいます。
ここでは、およそ100年に一度、メタンハイドレートが大量発生することで氷が融け、その中のうまみがスープとなって流れ出ます。これこそ、「センチュリースープ」なのです。

このアイスヘルなのですが、流体力学的に気になる記述があります。

[トリコ　9巻より]

グルメショーウィンドーに熱源があるために局地風が生じているというのです。こういうのを見ると、このような流れ場ができるのか、数値シミュレーションで確かめたくなってしまいます。

そこで今回は、「思い立ったが吉日、それ以降は凶日」ということで、単純な流体モデルを用いてこのアイスヘルの大気の流れを数値計算で再現することを目指します。アイスヘルでの大気の流れ場が理解されることで、センチュリースープを捜索する大いなる助けになるでしょう。したがって、このことは、グルメ時代において非常に有用なものとなると思われます。(とはいえ、トリコ10巻の段階でグルメショーウィンドーのうまみは既に枯渇しているのですが)

数値モデル

ここでは、次のようなモデルで考えます。

• 非圧縮
• 2次元 (軸対称を想定)
• ブシネスク近似

トリコの世界といえども地球すなわち回転座標系なので、コリオリ力を無視しても良いのかとも思いますが、ここでは簡単のため一旦無視します。密度変化もわりと大きいかもしれませんが、一旦ブシネスク近似でいきます(密度変化をちゃんと考えるとポアソン方程式がめんどくさい)。

そもそも、トリコの世界の地球は我々のそれに比べて非常に大きいことや、アイスヘルの緯度がどのあたりなのかよくわからない(極域ではないのは確実、詳しくはコミックスを読んでください)ので、回転の効果を考えるといろいろ余計なことを考えないといけないのです(面倒)。

支配方程式

ナビエストークス方程式、連続の式、温度の移流拡散方程式を考えます。

\newcommand{\pdv}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdvh}[2]{\partial #1/\partial #2}
\newcommand{\bm}[1]{ \boldsymbol{#1} }
\begin{align}
\displaystyle
\pdv{u}{t} +u \pdv{u}{x} + w \pdv{u}{z} &= -\frac{1}{\rho_0} \pdv{p}{x} + \nu_h \pdv{^2 u}{x^2}+ \nu_v \pdv{^2 u}{z^2} \\
\pdv{w}{t} +u \pdv{w}{x} + w \pdv{w}{z} &= -\frac{1}{\rho_0} \pdv{p}{z} + \nu_h \pdv{^2 w}{x^2}+ \nu_v \pdv{^2 w}{z^2} - \frac{\rho'}{\rho_0}g \\
\pdv{u}{x} + \pdv{w}{z} &= 0 \\
\pdv{T}{t} +u \pdv{T}{x} + w \pdv{T}{z} &= \kappa_h \pdv{^2 T}{x^2}+ \kappa_v \pdv{^2 T}{z^2}
\end{align}

ここで、

• $x$: 水平方向
• $z$: 鉛直方向
• $u$: 水平風速
• $w$: 鉛直風速
• $p$: 圧力
• $T$: 温度
• $\rho_0$: 基準密度
• $\rho'$: 基準密度からのずれ($\rho=\rho_0+\rho'$)
• $\nu_h$: 水平渦粘性係数
• $\nu_v$: 鉛直渦粘性係数
• $\kappa_h$: 水平渦拡散係数
• $\kappa_v$: 鉛直渦拡散係数

としています。ここで、渦度 $\displaystyle\zeta = \pdv{w}{x}-\pdv{u}{z}$として変形を行なって

\begin{align}
\displaystyle
\pdv{\zeta}{t} +u \pdv{\zeta}{x} + w \pdv{\zeta}{z} =  \nu_h \pdv{^2 \zeta}{x^2}+ \nu_v \pdv{^2 \zeta}{z^2} - \frac{g}{\rho_0}\pdv{\rho'}{x}
\end{align}

が得られます。ここで、単純な線型の状態方程式 $\rho=\rho_0(1-\alpha(T-T_0))$を考えて、$\pdvh{\rho'}{x}=-\rho_0 \alpha \pdvh{T}{x}$が得られます。
また、$\nabla\cdot\bm{u}=0$より、流線関数 $\psi: u=-\pdvh{\psi}{z}, w=\pdvh{\psi}{x}$を定義できる。したがって、

\begin{equation}
\displaystyle
\zeta = \pdv{w}{x}-\pdv{u}{z} = \pdv{}{x}\pdv{\psi}{x}+\pdv{}{z}\pdv{\psi}{z}=\Delta \psi \end{equation}

式をまとめると

\begin{align}
\displaystyle
\Delta \psi &= \zeta \\
\pdv{\zeta}{t} -\pdv{\psi}{z} \pdv{\zeta}{x} + \pdv{\psi}{x} \pdv{\zeta}{z} &=  \nu_h \pdv{^2 \zeta}{x^2}+ \nu_v \pdv{^2 \zeta}{z^2} + \alpha g \pdv{T}{x} \\
\pdv{T}{t} -\pdv{\psi}{z} \pdv{T}{x} + \pdv{\psi}{x} \pdv{T}{z} &=  \kappa_h \pdv{^2 T}{x^2}+ \kappa_v \pdv{^2 T}{z^2}
\end{align}

これらの式を数値的に解いていきます。

計算条件

離散化

計算はとりあえず次のようにしました。

• 空間差分化
  • 移流項：1次の風上差分
  • その他：2次の中心差分
• 時間積分：1次のオイラー法
• 物理量の配置：格子の交わるところに配置

境界条件・初期条件

• 滑り壁、すなわち$\zeta=0$ @ すべての壁
• 流線関数は$\psi=0$ @ すべての壁 としておく
• 温度は各壁によって変える
  • 右と上：フラックスなし
  • 下：等温 ($T=T_\mathrm{min}$)
  • 左(上側)：フラックスなし
  • 左(下側)：関数で与える(ガウス型)
• 初期条件は$\psi=0, \zeta=0, T=T_\mathrm{min}$

パラメーター

• $L_x=10000$m, $L_z=4000$m, $H_m=2000$m
• $n_x=100, n_z=40, $ すなわち $\Delta x=\Delta z=100$m
• $g=10$m/s^2, $\alpha=0.003$/K, $T_\mathrm{max}=150$℃, $T_\mathrm{min}=-50$℃(アイスヘルの平均気温は-50℃らしい)
• $f(z)=(T_\mathrm{max}-T_\mathrm{min})\exp(-(z-0.8H_m)^2/L_m^2)+T_\mathrm{min},$ $ L_m=0.1L_z$
• $\nu_h=\nu_v=\kappa_h=\kappa_v=100$m^2/s (計算時間短縮のためとりあえず大きい値)
• $\Delta t= \mathrm{min}(0.5\mathrm{d}x^2/\nu, 0.5\mathrm{d}z^2/\nu, 0.5\mathrm{d}x/\mathrm{max}(u), 0.5\mathrm{d}z/\mathrm{max}(w) )$

水平距離はトリコたちが帰り道に短時間で走り抜けてるのでこれくらいとしておく。

結果

500000ステップで976559秒(大体10日とかそこら)計算を行ないました。計算時間は実時間で5分くらいでした。コードはJuliaで書きました(コードは掲載するかもしれない)。

上から、温度、流線関数、風速、水平風速をプロットしました。(だいたい定常っぽくなっていました)

考察

風速場を見ると、たしかに風が熱源のほうに集まって、上昇流と下降流に分かれている様子がみられました。
温度場では、左の壁の熱源からの拡散よりも、鉛直流による移流が卓越していそうに見えます。上に熱が運ばれてそこで循環が起こり混合しているようにも見えます。

しかし、下降流によって現れる地表付近の風速は10m/sにも満たない程度で、トリコ本編の記述「風速30m/s」には達しませんでした。これに対していくつか考えられるのは

• パラメーターがおかしい
  今回の数値実験ではパラメータを多く設定しているので(熱膨張率や境界の温度、水平スケールなど)、その影響が出ているのかもしれません。また、地表の風速を求めたいという点では、鉛直方向のグリッドの大きさを地表付近でより高解像度にする必要があるかもしれません。
• 支配方程式がおかしい
  温度差が200Kくらいあるのにブシネスク近似が成り立つのか？という問題もあります。また、今回は軸対象として2次元での実験を行ないましたが、現実はやはり3次元なのでその効果を検証する余地も残っています。あとは境界条件も考慮しないといけないかもしれません。
• トリコの世界がおかしい
  そもそも漫画の世界なので、現実の物理が成り立っていないかもしれません。

まだまだ改善の余地がありますね。トリコの世界に倣ってこの学問分野を「グルメ流体力学(GFD)」とでも呼びたいものです。

おわりに

この数値計算は一見ふざけているようにも見えるかもしれませんが、フィクションの物理をシミュレーションで遊んでみるのも楽しいものです。流体現象はとても身近なものなので、こういう検証のタネもたくさんあります。最近はオープンソースのソフトも結構ありますしみなさんも流体シミュレーションで遊んでみましょう。

また、トリコは非常に面白い漫画です。「食」をテーマにあのような漫画が描ける島袋先生は偉大です。みなさんも余裕があればぜひ読んでください。アニメのリメイクも期待しています。

これは単なる趣味の数値計算ですので、詰めが甘い部分もあると思います。何か間違い等ございましたらお伝えください。

参考文献

• 島袋光年, "トリコ", 集英社