球面座標系におけるNavie - Stoke E.Q.の数値計算

1 SCHEME
Navie Stokes Equation on R
大田黒紘之 2015.5.29
1
Scheme
1.1
Basic
球座標系 (r,θ,ϕ) における非圧縮性流れの速度ベクトルを v(vr , vθ , vϕ ) とする.
連続の式及び運動方程式 (Navie Stokes Equation) は以下のように表現できる.
1 ∂ ( 2 )
1 ∂
1 ∂vϕ
r vr +
(sinθvθ ) +
=0
2
r ∂r
rsinθ ∂θ
rsinθ ∂ϕ
(1)
∂vr
∂vr vθ ∂vr
vϕ ∂vr vθ2 + vϕ2
1 ∂p
+ vr
+
+
−
=−
+R
∂t
∂r
r ∂θ
rsinθ ∂ϕ
r
ρ ∂r
(
)
2vr
2 ∂vθ 2vθ cotθ
2 ∂vϕ
+ ν ∇2 vr − 2 − 2
−
− 2
2
r
r ∂θ
r
r sinθ ∂ϕ
(2)
∂vθ
∂vθ vθ ∂vθ
vϕ ∂vθ vr vθ vp2 cotθ
1 ∂p
+ vr
+
+
+
−
=−
+Θ
∂t
∂r
r ∂θ
rsinθ ∂ϕ
r
r
ρr ∂θ
(
)
2 ∂vr
vθ
2cosθ ∂vϕ
+ ν ∇2 vθ + 2
− 2 2 − 2 2
r ∂θ
r sin θ r sin θ ∂ϕ
(3)
∂vϕ
∂vϕ vθ ∂vϕ
vϕ ∂vϕ vr vϕ vθ vϕ cotθ
1 ∂p
+ vr
+
+
+
+
=−
+Φ
∂t
∂r
r ∂θ
rsinθ ∂ϕ
r
r
ρrsinθ ∂ϕ
(
)
vϕ
2
∂vr
2cosθ ∂vθ
2
+ ν ∇ vϕ − 2 2 + 2 2
+ 2 2
r sin θ r sin θ ∂ϕ
r sin θ ∂ϕ
(4)
ただし,
1 ∂
∇ = 2
r ∂r
2
(
)
(
)
1
∂
∂
1
∂2
2 ∂
r
+ 2
sinθ
+ 2 2
∂r
r sinθ ∂θ
∂θ
r sin θ ∂ϕ2
(5)
1
1
R, Θ, Φ は外力を表す項、P は、圧力を示す項、ν は動粘度であり、ρ は密度を示すパラメータである.
1
1.2
1.2
流速予測子 v ∗ (vr∗ , vθ∗ , vϕ∗ ) の導入
1 SCHEME
流速予測子 v ∗ (vr∗ , vθ∗ , vϕ∗ ) の導入
球面座標系におけるナビエ・ストークス方程式を以下のように変形する.
∂vr
∂vr vθ ∂vr
vϕ ∂vr vθ2 + vϕ2
1 ∂p
= −vr
−
−
+
−
+R
∂t
∂r
r ∂θ
rsinθ ∂ϕ
r
ρ ∂r
(
)
2vr
2 ∂vθ 2vθ cotθ
2 ∂vϕ
2
+ ν ∇ vr − 2 − 2
−
− 2
r
r ∂θ
r2
r sinθ ∂ϕ
(6)
∂vθ
∂vθ vθ ∂vθ
vϕ ∂vθ vr vθ vp2 cotθ
1 ∂p
= −vr
−
−
−
+
−
+Θ
∂t
∂r
r ∂θ
rsinθ ∂ϕ
r
r
ρr ∂θ
(
)
2 ∂vr
vθ
2cosθ ∂vϕ
2
+ ν ∇ vθ + 2
− 2 2 − 2 2
r ∂θ
r sin θ r sin θ ∂ϕ
(7)
∂vϕ
∂vϕ vθ ∂vϕ
vϕ ∂vϕ vr vϕ vθ vϕ cotθ
1 ∂p
= −vr
−
−
−
−
−
+Φ
∂t
∂r
r ∂θ
rsinθ ∂ϕ
r
r
ρrsinθ ∂ϕ
(
)
vϕ
2
∂vr
2cosθ ∂vθ
2
+ ν ∇ vϕ − 2 2 + 2 2
+ 2 2
r sin θ r sin θ ∂ϕ
r sin θ ∂ϕ
(8)
以上の式を用いて、タイムステップ n を考慮した流速予測子 v ∗ (vr∗ , vθ∗ , vϕ∗ ) を定義すると、
(n)
vr∗(n+1)
∗(n+1)
vθ
=
vr(n)
∂vr
+ ∆t
∂t
(9)
(n)
(n)
= vθ + ∆t
∂vθ
∂t
(10)
(n)
∂vϕ
=
+ ∆t
∂t
タイムステップ n の流速から、タイムステップ n+1 の流速予測子が定義できる
∗(n+1)
vϕ
(n)
vϕ
2
(11)
1.3 圧力 P の定義
1.3
1 SCHEME
圧力 P の定義
流速予測子 v ∗ (vr∗ , vθ∗ , vϕ∗ ) を用いて、圧力 P を以下のように計算する. (連続の式を満た
す様な圧力項を探してくる為、物理的な圧力とは若干意味合いが異なる)
∇ · v ∗(n+1)
∆t
(
) (12)
∗(n+1)
)
∂v
1 ∂ (
1
1 ∂ ( 2 ∗(n+1) )
1
ϕ
∗(n+1)
r vr
+
=−
sinθvθ
+
∆t r2 ∂r
rsinθ ∂θ
rsinθ ∂ϕ
∇2 ϕ(n+1) = −
P n+1 = P n − ϕ(n+1)
(13)
また、以上の ϕ と流速予測子 v ∗ を用いて、タイムステップ n+1 の流速は以下のように
更新できる.
∂ϕ(n+1)
∂r
∂ϕ(n+1)
∗(n+1)
= vθ
+ ∆t
∂θ
∂ϕ(n+1)
∗(n+1)
= vϕ
+ ∆t
∂ϕ
vr(n+1) = vr∗(n+1) + ∆t
(n+1)
vθ
(n+1)
vϕ
1.4
(14)
(15)
(16)
まとめ
・球面座標系におけるナビエ・ストークス方程式を時間項について解き、時間方向の離
散化パラメータ ∆t を用いて流速予測子 v ∗ を計算
・流速予測子 v ∗ のダイバージェンスから、ϕ を求め次のタイムステップにおける圧力を
求める
・ϕ から次のタイムステップにおける流速を求める
1.5
Program
#include<iostream>
int main(int argc,char *argv[]){
//Under Constructing
return 0;
}
3

Download Report