この記事は格子法を説明した流体力学のプログラムを作りたい！その１記事の補足になります。
ナビエ・ストークス方程式の難しさの一つに、次元が増えるごとに速度情報とその偏微分成分が加速度的に増えることが挙げられる。
二次元で考えれば速度情報はx方向速度、y方向速度はもちろん、その2つのx軸方向偏微分、y軸方向偏微分で2*2=4つの成分を考えることになる。それが三次元ではx速度、y速度、z速度がさらにx,y,z軸方向偏微分で3*3=9の成分を考えなくてはいけなくなる。
普通に頭がこんがらがる。
だが避けては通れないので自分用のメモとしても、ここに整理しておこうと思う。
ある地点の風速がベクトルで表示されている。これがベクトル場。
ベクトル表記ようにx方向成分とy方向成分に分離することができる。
1ピクセルあたり濃淡の1変数しか情報を持たない。
さきほどの分離した成分
これもスカラー場と言える。
でここからがややこしいところだが、このスカラー場に∇を作用させるとベクトル場になったりベクトル場に∇を作用させるとスカラー場になったりする。
https://physics-school.com/grad-div-rot/
まずこれを知っていないと理解できない。
ナビエ・ストークス方程式を理解するのに必要なのは左二つの∇fと∇⋅A
早速∇fという操作をやってみる。
スカラー場 fは例によってlena画像を使う。
ここでは、全ピクセルの色情報についてx方向、y方向に差をとっている。具体的にはピクセル位置(x+1,y)から(x,y)の色を引いた値がx方向偏微分の値、(x,y+1)から(x,y)の色を引いた値がy方向偏微分の値ということになる。
そして、結果をみればわかるよう∇fの結果はベクトル場になった。(スカラー場を作用させたらベクトル場になった＝情報量が増えている)
として計算してみる。
ということで結果がスカラー場になった。色の加減算では、灰色(r,g,b=127,127,127)が±0として、白がプラス、黒がマイナスの値として計算している。
改めて2次元平面での∇演算子の定義をみてみると
とあり、これはベクトルの形である、ということのようだ。
確かに上の表にあるよう ベクトル*スカラー=ベクトルになるし、ベクトル・ベクトル=スカラー になる。
自分がしばらく理解できなかったのはこのうち
移流項の部分。
まずカッコ内のところからひも解いてみる。
これは速度場vが∇の右ではなく左にある形だ。なので∇は作用せず、単にベクトルの内積という計算をすればよい。
ここで速度場vのイメージは
v=
=(
,
)
こんな感じ。これをまずはv=(vx,vy)とし式を書くと
こうなる。もちろんこのvx,vyはスカラー場だ。
そして次にカッコの右側のvを作用させると
こんな感じのイメージになる。
ここで速度場の偏微分がでてきた。まずはx方向の偏微分をみてみると
こうなる。
ベクトル場の偏微分なので答えもベクトル場になる。
ここでは「速度場のx速度のx方向偏微分」と「速度場のy速度のx方向偏微分」を要素に持つベクトルを求めているわけだ。
同じように
ここでは「速度場のx速度のy方向偏微分」と「速度場のy速度のy方向偏微分」を要素に持つベクトルを求めている。
くどいようだが重要なのでここを間違いないようにしたい。
これが3次元だと9パターン・・・結構大変ではある。
それでは話を進め、これら4パターンの実際の偏微分の値が求まったとする。
これを上の式に代入すると
これを成分ごとにまとめ整理すると
こうなる。
ここでvx,vyは思い出してもらうと
こんな感じのスカラー場であったので、これも結局は1変数の数値。
矢印の長さを数値化すると
,
こうなる。
したがって上の式は最終的に
という計算をすることで求まることが分かった！
で移流項のコードを載せているが、これも上記の計算と同じことをやっていることが分かる。
if文でごちゃごちゃしているが、偏微分の計算のところは単に差をとっており、あとはさっきの画像の通りになっていることが分かるはずだ。
ナビエ・ストークス方程式の難しさの一つに、次元が増えるごとに速度情報とその偏微分成分が加速度的に増えることが挙げられる。
二次元で考えれば速度情報はx方向速度、y方向速度はもちろん、その2つのx軸方向偏微分、y軸方向偏微分で2*2=4つの成分を考えることになる。それが三次元ではx速度、y速度、z速度がさらにx,y,z軸方向偏微分で3*3=9の成分を考えなくてはいけなくなる。
普通に頭がこんがらがる。
だが避けては通れないので自分用のメモとしても、ここに整理しておこうと思う。
スカラー場とベクトル場
ベクトル場
ベクトル場で分かりやすい例として天気予報等ででてくる風速図ある地点の風速がベクトルで表示されている。これがベクトル場。
ベクトル表記ようにx方向成分とy方向成分に分離することができる。
x速度成分
y速度成分
スカラー場
一方スカラー場は、ある地点の情報が1つの変数で表せられるもの。例えばモノクロ色画像1ピクセルあたり濃淡の1変数しか情報を持たない。
さきほどの分離した成分
これもスカラー場と言える。
でここからがややこしいところだが、このスカラー場に∇を作用させるとベクトル場になったりベクトル場に∇を作用させるとスカラー場になったりする。
∇演算子の作用について
作用3パターン
∇(ナブラ)は表記の仕方によって3種類にわかれる。| 勾配 | 発散 | 回転 | |
|---|---|---|---|
| 表記 | ∇f | ∇⋅A | ∇×A |
| ∇の右に くるもの | スカラー場 | ベクトル場 | ベクトル場 |
| 計算結果 | ベクトル場 | スカラー場 | ベクトル場 |
まずこれを知っていないと理解できない。
ナビエ・ストークス方程式を理解するのに必要なのは左二つの∇fと∇⋅A
∇f
早速∇fという操作をやってみる。
スカラー場 fは例によってlena画像を使う。
ここでは、全ピクセルの色情報についてx方向、y方向に差をとっている。具体的にはピクセル位置(x+1,y)から(x,y)の色を引いた値がx方向偏微分の値、(x,y+1)から(x,y)の色を引いた値がy方向偏微分の値ということになる。
そして、結果をみればわかるよう∇fの結果はベクトル場になった。(スカラー場を作用させたらベクトル場になった＝情報量が増えている)
∇・A
次に∇⋅Aという操作をしてみる。ここでAがベクトル場でないといけないのでとして計算してみる。
ということで結果がスカラー場になった。色の加減算では、灰色(r,g,b=127,127,127)が±0として、白がプラス、黒がマイナスの値として計算している。
改めて2次元平面での∇演算子の定義をみてみると
とあり、これはベクトルの形である、ということのようだ。
確かに上の表にあるよう ベクトル*スカラー=ベクトルになるし、ベクトル・ベクトル=スカラー になる。
2Dナビエ・ストークス式の移流項をひも解く
自分がしばらく理解できなかったのはこのうち
移流項の部分。
まずカッコ内のところからひも解いてみる。
これは速度場vが∇の右ではなく左にある形だ。なので∇は作用せず、単にベクトルの内積という計算をすればよい。
ここで速度場vのイメージは
v=
=(,)こんな感じ。これをまずはv=(vx,vy)とし式を書くと
こうなる。もちろんこのvx,vyはスカラー場だ。
そして次にカッコの右側のvを作用させると
こんな感じのイメージになる。
ここで速度場の偏微分がでてきた。まずはx方向の偏微分をみてみると
こうなる。
ベクトル場の偏微分なので答えもベクトル場になる。
ここでは「速度場のx速度のx方向偏微分」と「速度場のy速度のx方向偏微分」を要素に持つベクトルを求めているわけだ。
同じように
ここでは「速度場のx速度のy方向偏微分」と「速度場のy速度のy方向偏微分」を要素に持つベクトルを求めている。
くどいようだが重要なのでここを間違いないようにしたい。
- x速度のx方向偏微分
- y速度のx方向偏微分
- x速度のy方向偏微分
- y速度のy方向偏微分
これが3次元だと9パターン・・・結構大変ではある。
それでは話を進め、これら4パターンの実際の偏微分の値が求まったとする。
これを上の式に代入すると
これを成分ごとにまとめ整理すると
こうなる。
ここでvx,vyは思い出してもらうと
こんな感じのスカラー場であったので、これも結局は1変数の数値。
矢印の長さを数値化すると
,こうなる。
したがって上の式は最終的に
という計算をすることで求まることが分かった！
一次風上差分のおさらい
流体力学のプログラムを作りたい！その１で移流項のコードを載せているが、これも上記の計算と同じことをやっていることが分かる。
vx,vyはそれぞれx方向速度、y方向速度が格納されていて、スタガード格子という前提がある。
vx_after,vy_afterには移流項計算後の値が入る。
vxの更新のコード
u=vx[i][j]
v=( vy[i-1][j] + vy[i][j] + vy[i-1][j+1] + vy[i][j+1] )/4
if (u>=0.0)&(v>=0.0){
vx_after[i][j] = vx[i][j] - u*(vx[i][j] - vx[i-1][j])*Δt - v*(vx[i][j]-vx[i][j-1])*Δt
}
if (u<0 .0)&(v>=0.0){
vx_after[i][j] = vx[i][j] - u*(vx[i+1][j] - vx[i][j])*Δt - v*(vx[i][j]-vx[i][j-1])*Δt
}
if (u>=0.0)&(v<0 .0){
vx_after[i][j] = vx[i][j] - u*(vx[i][j] - vx[i-1][j])*Δt - v*(vx[i][j+1]-vx[i][j])*Δt
}
if (u<0.0)&(v<0.0){
vx_after[i][j] = vx[i][j] - u*(vx[i+1][j] - vx[i][j])*Δt - v*(vx[i][j+1]-vx[i][j])*Δt
}
vy更新のコード
u=( vx[i][j-1] + vx[i+1][j-1] + vx[i+1][j] + vx[i][j] )/4
v=vy[i][j]
if (u>=0.0)&(v>=0.0){
vy_after[i][j] = vy[i][j] - u*(vy[i][j] - vy[i-1][j])*Δt - v*(vy[i][j]-vy[i][j-1])*Δt
}
if (u<0 .0)&(v>=0.0){
vy_after[i][j] = vy[i][j] - u*(vy[i+1][j] - vy[i][j])*Δt - v*(vy[i][j]-vy[i][j-1])*Δt
}
if (u>=0.0)&(v<0 .0){
vy_after[i][j] = vy[i][j] - u*(vy[i][j] - vy[i-1][j])*Δt - v*(vy[i][j+1]-vy[i][j])*Δt
}
if (u<0.0)&(v<0.0){
vy_after[i][j] = vy[i][j] - u*(vy[i+1][j] - vy[i][j])*Δt - v*(vy[i][j+1]-vy[i][j])*Δt
}
if文でごちゃごちゃしているが、偏微分の計算のところは単に差をとっており、あとはさっきの画像の通りになっていることが分かるはずだ。
