Kelvin-Helmholtz不安定の(単純な)数値実験例

おはこんばんにちはてれじょんです。

今年こそ記事を色々書きたいと宣言した年末から早くも3ヶ月ちょっとが過ぎてしまいました。皆様いかがお過ごしでしょうか。

私の方はM1になったり、地球連続体力学のTAになったりと、まあ色々やっています。あと最近はお絵描きの練習をしてたりします。一応毎日何かしら書くのを目標にしており、misskeyに模写を投稿したりしてます。たまに行くバーで会った作家さんの勧めで、モルフォ人体デッサンの「箱と円柱で書く」ってやつをやっています。何かおすすめがあれば教えてくれれば嬉しいです。

さて、前置きが長くなりましたが本題です。

いままで色々数値実験をしてきたのですが、実験設定とかをまとめていなかったため、謎の設定でされた計算の結果がサーバーに溜まっている状況になってしまいました。

色々と不便なので、やったそばからブログに残すことにしました。*1今回はKelvin-Helmholtz不安定(以下KH不安定)です。

支配方程式は以下の渦度方程式です。

$\displaystyle\frac{D\zeta}{Dt}=-\nu\nabla^8\zeta$

普通のKH不安定であれば圧縮性もあり、密度成層するような状況を考えるわけですが、ここでは単純な例ということで非圧縮で密度一定の流体を考えています。したがって、バロトロピックな不安定性だけが見られるという設定になっています。また、後述の通りスペクトル法で計算を行ったので、計算の安定のために高階粘性を導入しています。

さて、非圧縮流体なので流線関数 $\psi$ を導入することができます。

$\displaystyle u=-\frac{\partial\psi}{\partial y},~~v=\frac{\partial \psi}{\partial x}$

この流線関数を以下のように関数展開してスペクトル法で計算を行いました。
スペクトル法を用いるため、領域はx-yの二次元とし、両方向に周期境界条件を課しています。
$\displaystyle \psi(x,y)=\sum_{n=-N}^N\sum_{m=-M}^M\hat{s}_{n,m}\exp\left(\frac{2\pi inx}{L_x}\right)\exp\left(\frac{2\pi imy}{L_y}\right)$

その他の実験設定は以下のようになっています。

水平・垂直切断波数 $N,M$	170
水平・垂直格子点数 $J,K$	512
領域の水平長さ $L_x$	4
領域の垂直長さ $L_y$	2
高階粘性係数 $\nu$	$2.0\times 10^{-18}$
時間発展	4段4位Runge-Kutta (線型項は分離)
時間刻み幅	1/50 or 1/100

初期条件には以下のような速度場を用意しました。この速度場になるような流線関数を求め、そこに振幅が $1.0\times10^{-15}$ の乱数をすべての波数に与えました。このように摂動を与えると、発達速度の一番大きな不安定モードが発達し、その波数の波が卓越する様子を観察できます。
$\displaystyle u(x,y,0)=A\tanh\left(B\cos\pi y\right),~~~~v(x,y,0)=0$

さて、実験結果は以下のようになりました。まずは $A=1/3,~~B=3$ の結果です。
youtu.be

t=40くらいから不安定になり、波数3の構造が発達している様子が見えます。
また不安定が十分発達したあと、各々の渦は切離して孤立渦になっています。この孤立渦が相互作用して、乱流状態へと移行する様子も観察できます。

次は $A=1/5,~~B=5$ の結果です。

youtu.be

こちらもt=40くらいから不安定になっていますが、今度は波数5の構造が発達しています。
見える構造が変わった理由ですが、これはシアー層の幅が小さくなったことに起因しています。
詳しくは各自調べてほしい*2のですが、線型安定性解析より、波長の短い波ほど発達率が大きいことがわかります。ただ、波長が短すぎると粘性の効果で散逸されやすくなったり、そもそも空間スケールが小さいためかき混ぜるスケールも小さいなどの理由から、ある程度の長さを持つ波長の波が発達することになります。ところで同じく線型安定性解析より、摂動はシアーの境界からの距離に対して指数関数的に( $e^{-ky}$ のオーダーで)減衰することがわかります。そのため、シアー層の部分が厚くなると、波長の短い波はより減衰しやすくなり、最も発達率の大きな波の波長は更に長くなります。