キャビティ流れの有限差分法シミュレーション（3次元）

初めに

以前の投稿の内容を3次元に拡張する．

Methods

支配方程式

非圧縮性のNavier-Stokes方程式は，以下の通り（無次元化している）．

\begin{align}
    \nabla \cdot \mathbf{u} &= 0 \\
    \frac{\partial \mathbf{u}}{\partial t} 
    + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} 
    &= - \nabla p 
    + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u} 
    + \mathbf{f} \
\end{align}

ここで，$\mathbf{u}$は速度，$p$は圧力，$\mathbf{f}$は外力，$\mathrm{Re}$はReynolds数である．

以前と同様の離散化とする:

• 移流項　　: 3次精度風上差分（Kawamura&Kuwahara1984，Kawamura+1986）
• 粘性項　　: 2次精度中心差分
• 圧力勾配項: 2次精度中心差分

時間方向への離散化も以前と同様に，分離型解法の一つであるChorinの射影法（Chorin1968）を用いる．安定性条件（Courant+1967，Neumann&Richtmyer1950）を考慮して時間増分$\Delta t$を選ぶ．

解像度

まず，特徴長さ$L$と特徴流速$U$を$L = 1, \ U = 1$に固定し，以前の投稿と同様に，$\mathrm{Re} = 1,000$の流れを考えることにする．渦が一様で等方的であるとすると，Kolmogorovの長さスケール$\eta$は，

\begin{align}
    \eta 
    \simeq \mathrm{Re}^{-3/4}
    = 1,000^{-3/4}
    \simeq 5.62 \times 10^{-3}
\end{align}

程度であると概算される．特徴長さ$L=1$に対して，DNSを実施するには，$N \simeq 5,832,000$程度，或いはそれ以上の計算点が要求されることになるが，手元のラップトップでは既に現実的な選択肢ではない（2次元では$N \simeq 32,000$程度で済んだのに）．したがって，LESを実施する．$\Delta x = 1 \times 10^{-2}$とし，Smagorinskyの渦粘性モデル（Smagorinsky1963）を用いる（逆カスケードの表現力は無い）．モデル定数$C_s$¹には，Lilly1966が与えた理論値$C_s \simeq 0.17$では減衰が強すぎることが知られているので²，Deardorff1970が示した経験値$C_s = 0.1$を採用する．

Results

先行研究との比較

まずは，Jiang+1994，Wong&Baker2002と中央部断面での$x$方向（壁が動く方向）の速度を比較する．彼らの手法は，

• Jiang+1994: FEM - almost-uniform mesh (50x52x25), trilinear element
• Wong&Baker2002: FEM - solution-adapted nonuniform mesh (48×48×48), trilinear element

とされており，どちらも壁付近に向けて$h$-refinementしているようだ．

今回の計算では，$\mathrm{Re} = 100, 400$程度では，$\Delta x = 1 \times 10^{-2}$で凡そ全ての渦が解像できると考えられるので，渦粘性モデルは用いないことにした．

$\mathrm{Re}$	$100$	$400$	$1,000$
$u (z)$

凡そ合致するが，ピークを捉えきれていない．壁付近での$h$-refinementを無視すれば，現在の方法はJiang+1994，Wong&Baker2002よりも細かい解像度を用いているはずなので，追従しないのは非常に残念なことである．解像度を高めるなどしてFDMが収束する場所を調べるなどしたいが，計算資源の都合から叶わない．

速度は凡そ同じような分布になる，ということと，中央部断面での速度，圧力，渦度の分布が殆ど同じようになる，という結果から，大雑把だが，それらしい計算はできていそうだ，として次に進む．

Re = 1,000

$\mathrm{Re} = 1,000$の流れをシミュレートしたときの，中央部断面での物理量を可視化する．

	$\mathrm{Re} = 1,000$
$p$
$| \mathbf{u} |$
$\mathbf{u}$
$\boldsymbol{\omega}$

Re = 10,000

$\mathrm{Re} = 1,000$では，意外と静かな流れになった．（醍醐味の？）3次元的な渦の広がりを見てみたいので，もう少しReynolds数を大きくしてみる．

	$\mathrm{Re} = 10,000$
$p$
$| \mathbf{u} |$
$\mathbf{u}$
$\boldsymbol{\omega}$

Re = 100,000

既に解像できない領域が大きすぎるので最早遊びなのだが（LESでも乱流エネルギーの80%程度は解像するのが望ましい，とされている），もう少しだけ見た目の迫力がある動画を見たくなった（$100,000$という数字を選んだのは，後述の津波のReynolds数を意識してのことである．もちろん，Reynolds数を考えるべきかFroude数を考えるべきかとか，そもそも流れの性状が全く違うとか，話が違うところは数えきれないくらい沢山あるのだが）．

	$\mathrm{Re} = 100,000$
$p$
$| \mathbf{u} |$
$\mathbf{u}$
$\boldsymbol{\omega}$

終わりに

2次元から3次元に拡張しただけで，本質的な新しさは無いが，2次元ではどうにかなったシミュレーションも，3次元では計算コストの観点からかなり難しいことを実感した．

周辺の話題として，

• Ku+1987が，同Reynolds数の流れを2次元と3次元で計算し，両者の差分を議論している
• E&Liu1997が，渦度-ベクトルポテンシャルの定式化を与え，原始変数を用いたMAC法との比較を行っている

ことなどが個人的には面白かった．

Appendix

A. どの程度のReynolds数まで考えるか？

ところで，津波³のReynolds数やFroude数はどの程度なのだろう，と思ったので少し調べてみた．Exton+2019によると，2011年東北地方太平洋沖地震による津波のReynolds数やFroude数は，$\mathrm{Re} = 1 \times 10^{7} \sim 2 \times 10^{7}$，$\mathrm{Fr} = 0.5 \sim 1.0$程度と概算されるらしい．また，同論文には，

A large Reynolds number causes substantial mismatch with the corresponding scaled-down laboratory simulations, and the mismatch makes the outcomes from the experiments unreliable. Suppose we use a 1/40 scaled-down model, which is possible in a laboratory water tank. Then, to match the Froude number, the Reynolds number would be $5 \times 10^4 \sim 1 \times 10^5$ which is mismatched by 2 to 3 orders of magnitude.

という記述がある．これは，Iemura+2007が報告している，実験室での津波のReynolds数，$\mathrm{Re} = 1.2 \times 10^{5} \sim 1.8 \times 10^{5}$程度と合致する．また，Imamura+2008は，1771年明和の大津波に関して，実験室での津波遡上のReynolds数が$\mathrm{Re} = 1.0 \times 10^{4}$程度，実際の明和の大津波では$\mathrm{Re} = 1.0 \times 10^{7}$程度であったと推定しており，これもExton+2019の内容と合致する．

Froude数とReynolds数を同時に合致させることは難しいので，どちらを合わせるかを選ぶことになる．津波は重力で駆動されると考えられるので，Froude数を合わせた実験をすることになり，したがって，実験室でのReynolds数は$\mathrm{Re} = 1.0 \times 10^{4} \sim 1.0 \times 10^{5}$程度になるのだろう．実験室スケールでさえ，きちんとシミュレートするのは非常に難しいことなのだと感じる．

B. 高Reynolds数では対称性を失うか？

まず，数値実験において，計算格子は規則正しく，偏りなく配置したことを強調する．

その上で，本文の数値実験では，$\mathrm{Re} = 1,000$では対称性がある流れ場に落ち着いたが，$\mathrm{Re} = 10,000$では対称性が破れた．これは，現象自体の本質的にchaoticな性質（今回の文脈では，乱流）を，数値計算が再現したものと考える．この観測から，対称性が失われるReynolds数は，今回のキャビティ流れの枠組みでは，$1,000 \le \mathrm{Re} \le 10,000$の範囲にあると予想できる．

Jiang+1994は，$\mathrm{Re} \lt 3,200$程度ならば，対称性がある流れ場に落ち着いた，と述べている．$1,000 \lt \mathrm{Re} \lt 10,000$の範囲を少しだけ調べる．全ての場合で，$\Delta x = 1 \times 10^{-2}, \ \Delta t = 1.33 \times 10^{-3}$の解像度で，$C_s = 0.1$のSmagorinskyのモデルを用い，無次元化時間$60$までのシミュレーションを実施した．

$\mathrm{Re}$	$\boldsymbol{\omega}$
$2,000$
$4,000$
$6,000$
$8,000$

手元の計算結果として，$\mathrm{Re} \le 2,000$の間では対称性が残り，$\mathrm{Re} \ge 4,000$では対称性が失われた．即ち，criticalなReynolds数は，$2,000 \le \mathrm{Re} \le 4,000$にあると考えられる．これは概ね，Jiang+1994と合致する．ただし，より高解像度化したり，より長時間のシミュレーションを行ったりすれば，話は変わってくる可能性はある．

どのReynolds数で対称性を失うか，には特別な興味が無い．寧ろtakeawayとして持ち帰りたいのは，対称でなくなる場合がある，ということである．そしてそれは，計算格子を十分規則正しく配置しても起きる，と言うことを強調する．計算資源の都合などから，対称性を加味して1/2モデルや1/4モデルなどを用いる場合があるが，そのようなことができるのは，限られた問題についてのみである，と分かる．

1. $C_s$の定義は，SGSの渦の特徴長さ$\ell_0$とグリッド幅$\Delta$の比である（$C_s = \ell_0 / \Delta$）． ↩

2. Lilly1966の解析は，壁面から十分遠方のfree stream regionが対象である． ↩

3. 本稿とは明らかに違う問題設定だが，突然話題に出したのは，所属研究室で津波が話題に上がることが多いためである． ↩