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

# 初めに
[以前の投稿](https://qiita.com/ShotaDeguchi/items/93699a35eada106239d7)の内容を3次元に拡張する．

# Methods
## 支配方程式
非圧縮性のNavier-Stokes方程式は，以下の通り（無次元化している）．
```math
\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](https://doi.org/10.2514/6.1984-340)，[Kawamura+1986](https://doi.org/10.1016/0169-5983(86)90014-6)）
* 粘性項　　: 2次精度中心差分
* 圧力勾配項: 2次精度中心差分

時間方向への離散化も以前と同様に，分離型解法の一つであるChorinの射影法（[Chorin1968](https://doi.org/10.2307/2004575)）を用いる．安定性条件（[Courant+1967](https://doi.org/10.1147/rd.112.0215)，[Neumann&Richtmyer1950](https://doi.org/10.1063/1.1699639)）を考慮して時間増分$\Delta t$を選ぶ．

## 解像度
まず，特徴長さ$L$と特徴流速$U$を$L = 1, \ U = 1$に固定し，以前の投稿と同様に，$\mathrm{Re} = 1,000$の流れを考えることにする．渦が一様で等方的であるとすると，Kolmogorovの長さスケール$\eta$は，
```math
\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](https://doi.org/10.1175/1520-0493(1963)091%3C0099:GCEWTP%3E2.3.CO;2)）を用いる（逆カスケードの表現力は無い）．モデル定数$C_s$[^1]には，[Lilly1966](http://dx.doi.org/10.5065/D62R3PMM)が与えた理論値$C_s \simeq 0.17$では減衰が強すぎることが知られているので[^2]，[Deardorff1970](https://doi.org/10.1017/S0022112070000691)が示した経験値$C_s = 0.1$を採用する．

[^1]: $C_s$の定義は，SGSの渦の特徴長さ$\ell_0$とグリッド幅$\Delta$の比である（$C_s = \ell_0 / \Delta$）．
[^2]: Lilly1966の解析は，壁面から十分遠方のfree stream regionが対象である．

# Results
## 先行研究との比較
まずは，[Jiang+1994](https://doi.org/10.1016/0045-7825(94)90172-4)，[Wong&Baker2002](https://doi.org/10.1002/fld.204)と中央部断面での$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)$ | ![vel_comparison.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/c7b1e971-7f6b-2a7a-4dae-2f0c4fa0bb52.png) | ![vel_comparison.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/82c06e63-cffc-0b48-7e30-f44f38a566fb.png) | ![vel_comparison.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/a3a90ab9-4231-6921-c139-948308165527.png) |

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

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

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

| | $\mathrm{Re} = 1,000$ |
|:-:|:-:|
| $p$                   | ![Re1000_p.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/7425fa9e-b49b-42ad-8cde-6760cd4299cd.gif) |
| $\| \mathbf{u} \|$    | ![Re1000_u.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/227812ac-8d29-b084-fe27-ef3841793f9d.gif) |
| $\mathbf{u}$          | ![Re1000_q.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/1872be5b-41da-1010-c1f3-d389c13fc3aa.gif) |
| $\boldsymbol{\omega}$ | ![Re1000_v.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/2c007963-273d-c5a6-a3e7-011b723b1e35.gif) |

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

| | $\mathrm{Re} = 10,000$ |
|:-:|:-:|
| $p$                   | ![ezgif.com-animated-gif-maker.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/25c2116a-fe6b-6e2d-f263-41debea26216.gif) |
| $\| \mathbf{u} \|$    | ![ezgif.com-animated-gif-maker (2).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/b586835c-e016-8895-8a14-a7146bda26ce.gif) |
| $\mathbf{u}$          | ![ezgif.com-animated-gif-maker (1).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/eb2f8137-210b-b296-b888-4abe843a994c.gif) |
| $\boldsymbol{\omega}$ | ![ezgif.com-animated-gif-maker (3).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/f31c3d46-fa60-8e1f-6e92-83ead3513a19.gif) |

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

| | $\mathrm{Re} = 100,000$ |
|:-:|:-:|
| $p$                   | ![ezgif.com-animated-gif-maker.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/4a7e3d2a-fc7f-6673-5164-dafd4bcfb708.gif) |
| $\| \mathbf{u} \|$    | ![ezgif.com-animated-gif-maker (2).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/adf3ca23-15ec-a48a-5d8f-014a61b81870.gif) |
| $\mathbf{u}$          | ![ezgif.com-animated-gif-maker (1).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/0544273b-3893-c8ae-3f3c-3c4568e88e15.gif) |
| $\boldsymbol{\omega}$ | ![ezgif.com-animated-gif-maker (3).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/029e71a4-efa1-bb6b-e203-43bd3f150166.gif) |

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

周辺の話題として，

* [Ku+1987](https://doi.org/10.1016/0021-9991(87)90190-2)が，同Reynolds数の流れを2次元と3次元で計算し，両者の差分を議論している
* [E&Liu1997](https://doi.org/10.1006/jcph.1997.5815)が，渦度-ベクトルポテンシャルの定式化を与え，原始変数を用いたMAC法との比較を行っている

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

# Appendix
## A. どの程度のReynolds数まで考えるか？
ところで，津波[^3]のReynolds数やFroude数はどの程度なのだろう，と思ったので少し調べてみた．[Exton+2019](https://doi.org/10.1038%2Fs41598-019-47512-x)によると，2011年東北地方太平洋沖地震による津波のReynolds数やFroude数は，$\mathrm{Re} = 1 \times 10^{7} \sim 2 \times 10^{7}$，$\mathrm{Fr} = 0.5 \sim 1.0$程度と概算されるらしい．また，同論文には，
[^3]: 本稿とは明らかに違う問題設定だが，突然話題に出したのは，所属研究室で津波が話題に上がることが多いためである．

> 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](https://doi.org/10.11532/proee2005a.29.902)が報告している，実験室での津波のReynolds数，$\mathrm{Re} = 1.2 \times 10^{5} \sim 1.8 \times 10^{5}$程度と合致する．また，[Imamura+2008](https://doi.org/10.1029/2007JC004170)は，[1771年明和の大津波](https://www.jma-net.go.jp/ishigaki/know/jishin/meiwa.html)に関して，実験室での津波遡上の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](https://doi.org/10.1016/0045-7825(94)90172-4)は，$\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$ | ![cavity_Re2000.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/cbc56716-0b11-a182-7082-15244e1bf407.gif) |
| $4,000$ | ![cavity_Re4000.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/b01f4089-ad50-bac1-bc37-86a8d66a01a9.gif) |
| $6,000$ | ![cavity_Re6000.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/39006d7d-b4af-d830-89c0-fb37303b9464.gif) |
| $8,000$ | ![cavity_Re8000.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/7b67e417-6a30-ac0b-e770-812f6234d3a1.gif) |

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

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