流体シミュレーション研究ログ No.1
ポテンシャル流れ
『渦度0』の理想的な流れをポテンシャル流れと呼ぶ．このポテンシャル流れを解く為には速度ポテンシャルφと流れ関数ψの２つのスカラ関数が必要になる．
流れ関数
$$ \frac{\partial {\Psi}}{\partial y} = u \tag{1.1} $$
$$ \frac{\partial {\Psi}}{\partial x} = -v\tag{1.2} $$
流れ関数とは
流れ場が\(u(x, y, t), v(x, y, t)\)であるとする．この時ある地点\((x,y)\)について考える． ここに浮きをおいた場合、微小時間後、\((dx, dy)\)だけ移動する． それぞれの移動量は流れ場\(u(x, y, t), v(x, y, t)\)と平行なので微小時間を算出して見ると
$$ \frac{dx}{u} =\frac{dy}{v} $$
つまり $$ -vdx+udy = 0 $$
実際の微少時間はわからないので具体的な移動距離\((dx, dy)\)を知ることはできないが、移動した方向だけ全微分を使うことで推測することが出来る．
先ほどの式を全微分の形に変形していく
$$ d{\Psi} = \frac{{\partial}{\Psi}}{ {\partial}{x} } dx + \frac{{\partial}{\Psi}}{ {\partial}{y} } dy \tag{1.3} $$
つまりこの勾配ベクトルが流速となるスカラー関数\({\Psi(x, y)}\)が流れ関数である .
流れ関数と速度
流れ関数の勾配ベクトルと速度の内積を考える
$$ \boldsymbol{U} {\cdot} {\Delta {\Psi}} = u \frac{\partial {\Psi}}{\partial x} + u \frac{\partial {\Psi}}{\partial y} $$
(1.1), (1.2)より $$ \boldsymbol{U} {\cdot} {\Delta {\Psi}} = -uv + uv \\ = 0 $$
よって流れ関数の勾配方向（等高線にする垂直な方向）と速度ベクトルが垂直、つまり流れ関数の等高線(流線)と速度ベクトルは平行
流線
流れ関数の一定な値を結んだ線(式(1.3)が0となる線上)、また全ての位置での接線は速度ベクトルに一致する
速度ポテンシャル定義
これを満たせば速度ポテンシャルであると言える． $$ \frac{\partial {\Phi}}{\partial x} = u \tag{1.4} $$
$$ \frac{\partial {\Phi}}{\partial y} = v \tag{1.5} $$
速度ポテンシャルと速度
速度ポテンシャルの勾配ベクトルが速度ベクトルに相当するため、速度ポテンシャルの等高線と速度ベクトルは直行する．
速度ベクトルの向きは速度ポテンシャルが増加する方向である．(イメージとは逆なので注意)
また速度ポテンシャルが存在するとき、「回転運動をしない」、「流体中に速度ゼロの所がない」
、「流体全体が動いてる
流れ関数と速度ポテンシャル
上記の導出より、流れ関数と速度ポテンシャルが直行関係にあることがわかる．
複素ポテンシャル
φ,ψは直行しているので、これらの変数で構成されている複素ポテンシャルは以下のようになる．
$$ W = {\Phi} + i{\Psi} $$
一様流れ
一様流れとは速度ベクトルが一定の流れのことを言う．x方向の一定速度を\(U\)とする
$$ U = \frac{d{\phi}} {dx} $$
$$ U = \frac{d{\psi}} {dy} $$
$$ {\phi} = Ux, {\psi} = Uy $$
$$ W = {\phi} + i{\psi} = Ux + iUy = U(x + iy) = Uz $$
ポテンシャル線描画アルゴリズム
目標
すでにグリッド上の点にはポテンシャルの値が算出されている．そのグリッド上でポテンシャルの値が等しい座標を辿って線を引いて行きたい．
到達へのイメージ
四角形の中で出来るだけ滑らかな曲線を引くためには、まず四角形を三角形4つに分解して考える必要がある．
そして『湧き出し』、『吸い込み』などの条件が加わってくるとgrid上のx, y座標それぞれが同じでもポテンシャルの値は変わってくるので一般的に考えると四角形を構成する４点のポテンシャルの値は異なります．
そこで三角形に分解したあとにこの三角形の領域ないでどう線分を引いたら求めたいポテンシャルの値をたどることが出来るかを考える. ①はa,bの比例関係で、②はa,cの比例関係で求まる．これを4箇所で繰り返す感じ．
- まずはじめに表示するポテンシャルの値のmax, minを決める, またポテンシャル線の本数、間隔の値を算出する．
int lineNum = res - 1; float range = 1.0; float minPotential = flowVelocity * -1.0 * range; float maxPotential = flowVelocity * range; //ポテンシャル線一本ごとに変化する値 float dp0 = (maxPotential - minPotential) /(lineNum+1);
- まわり4点のポテンシャルの値をまずはじめに求める．その後ポテンシャルの大小関係を比べて比例式に持ち込んでいる．
for (int k = 0; k < lineNum; k++) { pp = minPotential + (k+1) * dp0; for (int i = 0; i < res-1; i++) { for (int j = 0; j < res-1; j++) { pot[0] = phi[j + i * res ]; pos[0] = ofVec2f(j * rectDelta.x, i * rectDelta.y); pot[1] = phi[j + (i + 1) * res]; pos[1] = ofVec2f(j * rectDelta.x, (i + 1) * rectDelta.y); pot[2] = phi[(j + 1) + (i + 1) * res]; pos[2] = ofVec2f((j + 1)* rectDelta.x, (i + 1) * rectDelta.y); pot[3] = phi[(j + 1) + i * res]; pos[3] = ofVec2f((j + 1)* rectDelta.x, i * rectDelta.y); pot[4] = pot[0]; pos[4] = pos[0]; //center pot[5] = (pot[0] + pot[1] + pot[2] + pot[3]) / 4.0; pos[5] = ofVec2f((pos[0] + pos[1] + pos[2] + pos[3]) / 4.0); for (int m = 0; m < 4; m++) { pos1 = pos2 = ofVec2f(-10, -10); if ((pot[m] <= pp && pp < pot[m + 1]) || (pot[m] > pp && pp >= pot[m + 1])) { pos1 = pos[m] + (pos[m + 1] - pos[m]) * (pp - pot[m]) / (pot[m + 1] - pot[m]); } if ((pot[m] <= pp && pp <= pot[5]) || (pot[m] >= pp && pp >= pot[5])) { if (pos1.x < 0.0) { pos1 = pos[m] + (pos[5] - pos[m]) * (pp - pot[m]) / (pot[5] - pot[m]); } else { pos2 = pos[m] + (pos[5] - pos[m]) * (pp - pot[m]) / (pot[5] - pot[m]); if (type == 0) { ofSetColor(0.0, 0.0, 255.0); }else if (type == 1) { ofSetColor(0.0, 255.0, 0.0); } ofDrawLine(rectStart + pos1, rectStart + pos2); } } if ((pot[m + 1] <= pp && pp <= pot[5]) || (pot[m + 1] >= pp && pp >= pot[5])) { if (pos1.x < 0.0) { pos1 = pos[m + 1] + (pos[5] - pos[m + 1]) * (pp - pot[m + 1]) / (pot[5] - pot[m + 1]); } else { pos2 = pos[m + 1] + (pos[5] - pos[m + 1]) * (pp - pot[m + 1]) / (pot[5] - pot[m + 1]); ofDrawLine(rectStart + pos1, rectStart + pos2); } } }//m }//j }//i }//k
参考
流体方程式の解き方入門 www.amazon.co.jp
『流れ』と『波』のシミュレーション www.amazon.co.jp
