この記事は格子法を説明した流体力学のプログラムを作りたい!その1記事の補足になります。
ナビエ・ストークス方程式の難しさの一つに、次元が増えるごとに速度情報とその偏微分成分が加速度的に増えることが挙げられる。
二次元で考えれば速度情報はx方向速度、y方向速度はもちろん、その2つのx軸方向偏微分、y軸方向偏微分で2*2=4つの成分を考えることになる。それが三次元ではx速度、y速度、z速度がさらにx,y,z軸方向偏微分で3*3=9の成分を考えなくてはいけなくなる。
普通に頭がこんがらがる。

だが避けては通れないので自分用のメモとしても、ここに整理しておこうと思う。

スカラー場とベクトル場

ベクトル場

ベクトル場で分かりやすい例として天気予報等ででてくる風速図
vec
ある地点の風速がベクトルで表示されている。これがベクトル場。

ベクトル表記ようにx方向成分とy方向成分に分離することができる。
vec

x速度成分

vecx

y速度成分

vecy

スカラー場

一方スカラー場は、ある地点の情報が1つの変数で表せられるもの。例えばモノクロ色画像
lena

1ピクセルあたり濃淡の1変数しか情報を持たない。
さきほどの分離した成分
vecx
これもスカラー場と言える。

でここからがややこしいところだが、このスカラー場に∇を作用させるとベクトル場になったりベクトル場に∇を作用させるとスカラー場になったりする。


∇演算子の作用について

作用3パターン

∇(ナブラ)は表記の仕方によって3種類にわかれる。
勾配発散回転
表記∇f∇⋅A∇×A
∇の右に
くるもの
スカラー場ベクトル場ベクトル場
計算結果ベクトル場スカラー場ベクトル場
https://physics-school.com/grad-div-rot/
まずこれを知っていないと理解できない。
ナビエ・ストークス方程式を理解するのに必要なのは左二つの∇fと∇⋅A

∇f


早速∇fという操作をやってみる。

lenasiki0
スカラー場 fは例によってlena画像を使う。

lenasiki1

ここでは、全ピクセルの色情報についてx方向、y方向に差をとっている。具体的にはピクセル位置(x+1,y)から(x,y)の色を引いた値がx方向偏微分の値、(x,y+1)から(x,y)の色を引いた値がy方向偏微分の値ということになる。

そして、結果をみればわかるよう∇fの結果はベクトル場になった。(スカラー場を作用させたらベクトル場になった=情報量が増えている)

∇・A

次に∇⋅Aという操作をしてみる。ここでAがベクトル場でないといけないので
lenasiki2

として計算してみる。
lenasiki3

lenasiki4

ということで結果がスカラー場になった。色の加減算では、灰色(r,g,b=127,127,127)が±0として、白がプラス、黒がマイナスの値として計算している。

改めて2次元平面での∇演算子の定義をみてみると
teigi

とあり、これはベクトルの形である、ということのようだ。
確かに上の表にあるよう ベクトル*スカラー=ベクトルになるし、ベクトル・ベクトル=スカラー になる。


2Dナビエ・ストークス式の移流項をひも解く

ns

自分がしばらく理解できなかったのはこのうち
nabiekao
移流項の部分。


まずカッコ内のところからひも解いてみる。
nabiekaokakko

これは速度場vが∇の右ではなく左にある形だ。なので∇は作用せず、単にベクトルの内積という計算をすればよい。
ここで速度場vのイメージは
v=vecimg2=(vecimg2x,vecimg2y)

こんな感じ。これをまずはv=(vx,vy)とし式を書くと

nabie1
こうなる。もちろんこのvx,vyはスカラー場だ。

そして次にカッコの右側のvを作用させると
nabie2

こんな感じのイメージになる。

ここで速度場の偏微分がでてきた。まずはx方向の偏微分をみてみると
nabie3

こうなる。
ベクトル場の偏微分なので答えもベクトル場になる。
ここでは「速度場のx速度x方向偏微分」と「速度場のy速度x方向偏微分」を要素に持つベクトルを求めているわけだ。

同じように
nabie4

ここでは「速度場のx速度y方向偏微分」と「速度場のy速度y方向偏微分」を要素に持つベクトルを求めている。

くどいようだが重要なのでここを間違いないようにしたい。
    • x速度x方向偏微分
    • y速度x方向偏微分
    • x速度y方向偏微分
    • y速度y方向偏微分
というように4パターンのスカラー場がでてくるのが分かると思う。
これが3次元だと9パターン・・・結構大変ではある。

それでは話を進め、これら4パターンの実際の偏微分の値が求まったとする。

nabie5

これを上の式に代入すると

nabie6


これを成分ごとにまとめ整理すると


nabie7
こうなる。
ここでvx,vyは思い出してもらうと

vecimg2x,vecimg2y

こんな感じのスカラー場であったので、これも結局は1変数の数値。
矢印の長さを数値化すると
vx,vy
こうなる。


したがって上の式は最終的に

nabie8


という計算をすることで求まることが分かった!


一次風上差分のおさらい

流体力学のプログラムを作りたい!その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文でごちゃごちゃしているが、偏微分の計算のところは単に差をとっており、あとはさっきの画像の通りになっていることが分かるはずだ。