$$
\newcommand{\pdv}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\bm}[1]{ \boldsymbol{#1} }
$$
はじめに
鴨川(高野川・賀茂川)です。
川を渡る飛び石があります。
この飛び石の中には亀や鳥の形をしたものがあって、なんだかかわいらしい。
飛び石とその間の川の流れを眺めると、下流側に砂が積もっていたりして割と面白いです。こんなのを見ると数値シミュレーションで再現して遊んでみたいな～と思ってしまいますよね。
ということで、今回は亀の形をした境界まわりの流れ(とりあえず2次元)のシミュレーションを行なってみます。
(砂とかの河床変動は一旦置いておいて)
境界の表現 (埋め込み境界法)
計算を行ないたいのですが、円柱とかならまだしもこんな亀の形の境界に合わせたグリッドを作るのなんて非常に面倒です。
困っていたら、埋め込み境界法(Immersed Boundary Method, IBM)といういい感じの手法を見つけました。
IBMでは、境界の周りで流速を境界の速度(静止ならゼロ)に近づける強制を与えることで、境界を表現します(Direct Forcing)。
この方法だと、境界の位置が分かれば、その位置の近傍のグリッドに適当な分布関数で得られる外力を与えればよいので、直交格子のままでも計算ができます。
支配方程式
とりあえず非圧縮・2次元のナビエストークス方程式を考えます。(浅水方程式とかでもいいかもしれない)
\begin{align*}
\pdv{\bm{u}}{t} + (\bm{u}\cdot\nabla) \bm{u} &= -\nabla p + \frac{1}{Re}\Delta \bm{u} + \bm{F} \\
\nabla\cdot\bm{u} &= 0
\end{align*}
$\bm{u}=(u,v)$は流速、$p$は圧力、$Re$はレイノルズ数で、方程式は無次元化してあります。$\bm{F}$は外力で、今はダイレクトフォーシングによる埋め込み境界法の強制を考えます。
粘着条件を考えて、
\begin{gather}
\bm{F} = \frac{1}{\Delta t} ( \bm{u}_w - \bm{u} ) \delta_w(\bm{x}_w-\bm{x}).
\end{gather}
ここで $\bm{u}_w$ は移動境界の速度、$\bm{x}_w$は移動境界の境界位置。この外力項は、境界において、境界の速度に緩和させるような項で、静止境界の場合は$\bm{u}_w=0$となります。
$\delta(\bm{x})$はデルタ関数で、境界でのみこの外力がはたらくようにします(数値計算時は分布関数で与える)。
これらの方程式を離散化して、差分法で計算し時間発展させていきます。
離散化
Arakawa C グリッドで離散化を行ないます。
空間差分は、とりあえず移流項は5次のWENO、その他は2次の中心差分で実装しました。
時間発展は3次のルンゲクッタ法(SSP-RK3)を使いました。
フラクショナルステップ法で $u,v,p$の計算を行なって、圧力のポアソン方程式は離散コサイン変換で解きました。
(このあたりの説明は省略します、すみません)
外力項は、次のような離散化を施しました。
\begin{gather*}
\bm{u}_{ib}=\sum_{i,j} \bm{u}_{i,j} \delta(\bm{x}_{i,j}-\bm{x}_w) \\
\bm{f}_{ib}=\frac{1}{\Delta t} ( \bm{u}_w - \bm{u}_{ib} ) \\
\bm{F}_{i,j} = \sum_{n=1}^N (\bm{f}_{ib})_n \delta(\bm{x}_{i,j}-(\bm{x}_w)_n )
\end{gather*}
下付き添え字$ib$は Immersed Boundary を表しており、$u_{ib}$ は壁の移動速度ではなく、境界点での内挿された流速です。$f_{ib}$はその境界点での外力となっています。
IBMでは境界を点の集合として表しますので、$N$はその境界点の個数をあらわします。
すなわち、わかりづらいのですが、$u_{ib},f_{ib},x_{w}$は$N$次元のベクトルとなっており、各成分は境界を構成する点の流速・外力・座標を表します。
ですから、この外力の計算の流れは、以下のようになります。呼び方の都合上、物理量を計算する点($i,j$で表すグリッドにより決められる点)を計算点、境界を表す点($N$個の点)を境界点と呼びます。
ここで、補間関数 $\delta(\bm{x}_{i,j}-\bm{x}_w)$は次のようなものを用いました(Griffith, B. E., & Peskin, C. S. 2005)。
\begin{gather*}
\delta(\bm{x}_{i,j}-\bm{x}_w)= D \left(\frac{x_i-x_w}{dx}\right) D \left(\frac{y_j-y_w}{dy}\right)\\
D (r) \equiv \left\{
\begin{array}{ll}
\frac{1}{8}(3 - 2|r| +\sqrt{1+4|r|-4r^2}) & 0 \leq |r| < 1, \\
\frac{1}{8}(5 - 2|r| -\sqrt{-7+12|r|-4r^2}) & 1 \leq |r| < 2, \\
0 & 2 \leq |r|
\end{array}\right.
\end{gather*}
境界点の取得
さて、埋め込み境界法で計算する準備が整いました。あとは使いたい境界(今回は亀)の点の座標を取得するだけです。
しかし、亀石の形は円や四角形のように単純ではないため、境界点の取得が難しいのです。(直接計るわけにもいかない)
それなら画像から読み取ればいいじゃないかということで、インターネットを探したところ、Level Set 法(Chan-Vese level sets model)による画像のセグメンテーションという方法が使えそうだったので用いてみました。
レベルセット法による画像処理に関してはあまり詳しくないので細かい説明はできませんが、
- レベルセット関数(符号付距離関数)を、輝度値から得られるエネルギー関数を最小化するように反復計算を行う
という感じ(?)のようです。
とりあえず github (https://github.com/kevin-keraudren/chanvese) で公開されているPythonコードをこの亀の画像に用いてみると、次のような結果(レベルセット関数)が得られました。
google map の画質が粗いので細かい構造はわかりませんが、だいたいの輪郭は取得できている気がします。
これでレベルセット関数(符号付距離関数)の空間分布 $\phi(\bm{x})$ が得られたのはいいですが、これを(今考えている)埋め込み境界法に用いるためには境界の点を得る必要があります。
そこまで良い方法が思いつかなかったので、とりあえず愚直に
- 各$x$方向のインデックス $i$ に対して、$y$方向に $\phi=0$ となる点を探索
- 各$y$方向のインデックス $j$ に対して、$x$方向に $\phi=0$ となる点を探索
- 2点が近すぎる場合はそのうち1点を除去
というアルゴリズムで点を取り出しました。
雑なコード(の一部)：
ii,jj=np.shape(phi) # レベルセット関数のサイズ
xc=[] # 境界点のx座標
yc=[] # 境界点のy座標
dd=6 # 探索するインデックスの間隔
for i in range(1,ii,dd):
for j in range(jj-1):
if phi[i,j]*phi[i,j+1]<0:
xc.append(float(i))
yc.append(float(j+abs(phi[i,j])/(abs(phi[i,j])+abs(phi[i,j+1])) ))
for j in range(3,jj,dd):
for i in range(ii-1):
if phi[i,j]*phi[i+1,j]<0:
xc.append(float(i+abs(phi[i,j])/(abs(phi[i,j])+abs(phi[i+1,j]))))
yc.append(float(j))
N=len(xc)
pop_list=set() # 除去する点のインデックス
for i in range(0,N-1):
for j in range(i+1,N):
dis=np.sqrt( (xc[i]-xc[j])**2 + (yc[i]-yc[j])**2 )
if dis <= 4: # とりあえず距離4以下を除去
print(i,j)
pop_list.add(j)
pop_list=list(pop_list)
print(pop_list)
for i in sorted(pop_list, reverse=True):
xc.pop(i)
yc.pop(i)
これで得られた座標を適当に規格化してプロットしたら次のような図が得られました。
まあまあ等間隔で、ある程度亀の形を再現するくらいの点は得られていそうです。もっと細かいのが欲しかったら探索する間隔とか除去する2点間距離を変えればよいと思います。
ということで、境界点が得られたのでいざ数値実験です。
数値実験
境界条件・初期条件・パラメータ
とりあえず適当に決めました。
- 境界条件は以下の図のように設定し、亀は$x=2.5～ \ 5$あたりに置きました (現実の鴨川には亀以外の障害物もありますが一旦無視で)
- 左側境界：流入, $u=1,v=0,\pdv{p}{x}=0$
- 上下境界：滑りなし, $\pdv{p}{x}=0$
- 右側境界：$\pdv{u}{x}=\pdv{v}{x}=\pdv{p}{x}=0$
- 初期条件は、$u=1$、$v$は微小な擾乱：0.05*(rand()-0.5) のように実装
- グリッド数は $n_x=200, n_y=100$、すなわち $\Delta x = \Delta y = 0.1 $
- CFL条件を満たすように、$\Delta t = 0.01$
- わかりやすい現象が見たいので、カルマン渦列が出そうなレイノルズ数：$Re=100$と設定
結果
コードはJuliaで書きました(掲載するかは未定)。$t=100$ まで計算しました。
障害物(亀)の後ろにカルマン渦列ができているのが分かります。渦度の時間変化をプロットすると、以下のようになりました。渦がウズウズしていて面白いです。
問題
2次元非圧縮なので、
u=\pdv{\psi}{y}, \quad v=-\pdv{\psi}{x}
となるような流線関数$\psi$を定義できます。プロットすると、
障害物付近を拡大すると、流線関数が境界を貫通してしまっているのが見えます。
このような流線が境界を貫いてしまう問題は結構議論されているようで、その解決策も多く提案されているようです(詳しくはいろいろ論文を漁らないといけないので割愛)。このあたりも勉強できたら楽しそうです。
おわりに
とりあえず流れ場のみの計算を行なったので、次は温度場とかも入れられたらまた面白いかもしれません。また、モチベーションとなった河川の流れの再現のために、浅水方程式にIBMを導入して、河床変動のモデルを組み合わせても面白そうです。これからの課題です。
これは単なる趣味の数値計算ですので、詰めが甘い部分もあると思います。何か間違い等ございましたらお伝えください。
参考文献
- Griffith, B. E., & Peskin, C. S. (2005). On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems. Journal of Computational Physics, 208(1), 75-105.
- 瀬田剛. (2021). 格子ボルツマン法.
- Chan-Vese level-sets in Python: https://github.com/kevin-keraudren/chanvese