初めに
以前の投稿の内容を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$1には，Lilly1966が与えた理論値$C_s \simeq 0.17$では減衰が強すぎることが知られているので2，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$の流れをシミュレートしたときの，中央部断面での物理量を可視化する．
| Results | |
|---|---|
| $p$ |  |
| $| \mathbf{u} |$ |  |
| $\mathbf{u}$ |  |
| $\boldsymbol{\omega}$ |  |
Re = 10,000
$\mathrm{Re} = 1,000$では，意外と静かな流れになった．（醍醐味の？）3次元的な渦の広がりを見てみたいので，もう少しReynolds数を大きくしてみる．
| Results | |
|---|---|
| $p$ |  |
| $| \mathbf{u} |$ |  |
| $\mathbf{u}$ |  |
| $\boldsymbol{\omega}$ |  |
Re = 100,000
既に解像できない領域が大きすぎるので最早遊びなのだが（LESでも乱流エネルギーの80%程度は解像するのが望ましい，とされている），もう少しだけ見た目の迫力がある動画を見たくなった（$100,000$という数字を選んだのは，後述の津波のReynolds数を意識してのことである．もちろん，Reynolds数を考えるべきかFroude数を考えるべきかとか，そもそも流れの性状が全く違うとか，話が違うところは数えきれないくらい沢山あるのだが）．
| Results | |
|---|---|
| $p$ |  |
| $| \mathbf{u} |$ |  |
| $\mathbf{u}$ |  |
| $\boldsymbol{\omega}$ |  |
終わりに
2次元から3次元に拡張しただけで，本質的な新しさは無いが，2次元ではどうにかなったシミュレーションも，3次元では計算コストの観点からかなり難しいことを実感した．
周辺の話題として，
- Ku+1987が，同Reynolds数の流れを2次元と3次元で計算し，両者の差分を議論している
- E&Liu1997が，渦度-ベクトルポテンシャルの定式化を与え，原始変数を用いたMAC法との比較を行っている
ことなどが個人的には面白かった．
Appendix A: どの程度のReynolds数まで考えるか？
ところで，津波3の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}$程度になるのだろう．実験室スケールでさえ，きちんとシミュレートするのは非常に難しいことなのだと感じる．