ナビエ・ストークス方程式

ナビエ・ストークス方程式
春日悠
2016 年 2 月 23 日
目次
1
レイノルズの輸送定理
1
2
質量保存の法則
2
3
運動方程式
3
4
角運動量保存の法則
4
5
ニュートン流体の構成方程式
5
6
ナビエ・ストークス方程式
7
1
レイノルズの輸送定理
物体 B があるとする．時刻 t に B の占める領域を Ω とする．空間の位置ベクトルを
x として， B 上のある量を ϕ( x, t) としよう．ϕ はスカラーでもベクトルでもテンソル
でもよい．ϕ の体積積分の物質時間導関数を，次式で定義する．
(∫
)
∫
∫
D
1
ϕ( x, t)dv = lim
ϕ( x, t + ∆t)dv −
ϕ( x, t)dv
(1)
Dt Ω
∆t→0 ∆t
Ω′
Ω
図 1: 物体の運動
1
ここで Ω′ は時刻 t + ∆t において B が占める領域を表す (図 1)．∆Ω = Ω′ − Ω とする
と，Ω′ = Ω + ∆Ω なので
D
Dt
∫
ϕ( x, t)dv
)
)}
{ (∫
(∫
∫
1
1
ϕ( x, t + ∆t)dv −
ϕ( x, t)dv +
ϕ( x, t + ∆t)dv
= lim
∆t
∆t→0 ∆t
Ω
Ω
∆Ω
(2)
Ω
右辺の最後の項は， B の境界面 Γ 上の各点が ∆t の間に通過した体積に ϕ をかけたも
のの総和に等しいと考えられる． B の速度を u とし，Γ の外向き法線ベクトルを n と
すると，Γ 上の微小面積 ds が通過した体積は dv = u · ndtds と書けるので
D
Dt
∫
Ω
∫
ϕdv =
Ω
∂ϕ
dv +
∂t
∫
Γ
ϕu · nds
右辺第 2 項にガウスの発散定理を適用して
}
∫ {
∫
D
∂ϕ
ϕdv =
+ ∇ · (ϕu) dv
Dt Ω
∂t
Ω
(3)
(4)
上式をレイノルズの輸送定理という．
2
質量保存の法則
物体 B の密度を ρ( x, t) とする． B の質量を m とすると
∫
m=
Ω
ρdv
(5)
質量保存の法則は ṁ = 0 と表される．レイノルズの輸送定理 (式 (4)) を用いて
ṁ =
D
Dt
∫
Ω
∫ {
=
Ω
ρdv
}
∂ρ
+ ∇ · (ρu) dv
∂t
(6)
=0
これより，次式が成り立つ．
∂ρ
+ ∇ · (ρu) = 0
∂t
上式を連続の式という．これは次のようにも書ける．
(7)
Dρ
+ ρ∇ · u = 0
(8)
Dt
非圧縮性流体の場合，流体の密度は時間的に変化しない (Dρ/Dt = 0) ため，連続の式
は次のように表される．
∇·u = 0
(9)
非圧縮性流体だからといって ∂ρ/∂t = 0 とは限らない (密度分布が一様である必要は
ない) ことに注意。
2
運動方程式
3
物体 B にかかる外力 f は，表面力を t ，物体力を b とすると，次式のように書ける．
∫
f =
∫
Γ
tds +
bdv
Ω
(10)
応力テンソルを σ とし，B の境界面 S の外向き法線ベクトルを n とすると，コーシー
の公式
t = σ·n
(11)
より
∫
f =
=
∫Γ
σ · nds +
Ω
∫
Ω
bdv
(∇ · σ + b) dv
(12)
ガウスの発散定理を用いた．
運動方程式は，運動量を p とすると
ṗ = f
と表される．運動量 p は
∫
p=
Ω
ρudv
(13)
(14)
であるので，レイノルズの輸送定理 (式 (4)) を用いて
D
Dt
∫
ρudv
Ω
{
}
∫
∂
=
(ρu) + ∇ · (ρuu) dv = f
Ω ∂t
ṗ =
(15)
式 (13) に式 (12) と (15) を代入すると，次式を得る．
∂
(ρu) + ∇ · (ρuu) = ∇ · σ + b
∂t
(16)
上式の左辺は次のように変形できる．
∂ρ
∂u
∂
(ρu) + ∇ · (ρuu) =
u+ρ
+ u∇ · (ρu) + ρu · ∇u
∂t
∂t
∂t
{
}
{
}
∂ρ
∂u
=
+ ∇ · (ρu) u + ρ
+ u · ∇u
∂t
∂t
Du
=ρ
Dt
(17)
連続の式および以下の関係を用いた．
∇ · (uu) = u∇ · u + u · ∇u
3
(18)
Du
∂u
=
+ u · ∇u
Dt
∂t
以上より，運動方程式は次のようにも書ける．
ρ
4
(19)
Du
= ∇·σ+b
Dt
(20)
角運動量保存の法則
座標系の原点に関する物体 B の角運動量保存の法則は，次式のように表される．
D
Dt
∫
Ω
x × ρudv =
∫
Γ
x × tds +
∫
Ω
x × bdv
(21)
ベクトル a, b の外積 a × b は，ϵijk をエディントンのイプシロン
ϵijk



1


= −1



0
(i, j, k が 1, 2, 3 の偶置換)
(i, j, k が 1, 2, 3 の奇置換)
(22)
(i, j, k のうちどれか 2 つが等しい)
とすると，総和規約を用いて a × b = ϵijk a j bk ei と表される．ここで ei は座標系の基
底ベクトルである．式 (21) の右辺第 1 項は
∫
Γ
x × tds =
∫
∫Γ
=
Γ
x × σ · nds
ϵijk x j σkl nl ei ds
= ϵijk
∫
∫
Γ
x j σkl nl dsei
∂
( x j σkl )dvei
Ω ∂xl
)
(
∫
∂σkl
= ϵijk
δjl σkl + x j
dvei
∂xl
Ω
(
)
∫
∂σkl
=
e dv
ϵijk σkj ei + ϵijk x j
∂xl i
Ω
∫ (
)
ϵijk σkj ei + x × ∇ · σ dv
=
= ϵijk
(23)
Ω
ここで δij はクロネッカーのデルタ
δij =

1
(i = j )
0
(i ̸ = j )
である．また，ガウスの発散定理を用いた．
4
(24)
したがって式 (21) は，レイノルズの輸送定理を用いて
}
∫ {
∫
∫
∂
x×
(ρu) + ∇ · (ρuu) dv = x ×
ϵijk σkj ei dv
(∇ · σ + b) dv +
Ω ∂t
Ω
Ω
}
∫ {
∫
∂
(ρu) + ∇ · (ρuu) − ∇ · σ − b dv =
ϵijk σkj ei dv
x×
Ω ∂t
Ω
(25)
左辺は運動方程式によって 0 なので，次式が得られる．
∫
ϵijk σkj ei dv = 0
(26)
ϵijk σjk = 0
(27)
ϵijk σjk = ϵikj σkj = −ϵijk σkj = ϵijk σkj = 0
(28)
ϵijk σjk = ϵijk σkj
(29)
Ω
これから
上式より
したがって
上式が成り立つには σij = σji でなければならない．つまり，応力テンソル σ は対称で
ある．
5
ニュートン流体の構成方程式
せん断応力が変形速度に比例すると想定できる流体を，ニュートン流体という．こ
こでは，ニュートン流体の構成方程式を導く．
変形速度テンソル D は，速度勾配テンソル L の対称成分を表す．
L = ∇u
D=
1
( L + LT )
2
(30)
せん断応力テンソル τ と変形速度テンソル D の関係は，次式で表される．
τij = Cijkl Dkl
(31)
静止流体の場合，圧力を p とすると，静水圧の等方性より応力テンソル σ は次式のよ
うに書ける．
σij = − pδij
(32)
運動する流体の場合は，上式にせん断応力が加わるから
σij = − pδij + τij
= − pδij + Cijkl Dkl
5
(33)
流体が等方であるとしよう．そのとき，テンソル C は等方テンソルである．4 階の
等方テンソル C の成分 Cijkl は，次式のように表すことができる．
Cijkl = λδij δkl + µ(δik δjl + δil δjk )
(34)
ここで，スカラー µ, λ はそれぞれ粘性係数，第 2 粘性係数と呼ばれる．上式を用いて，
応力と変形速度の関係は次式のようになる．
{
}
σij = − pδij + λδij δkl + µ(δik δjl + δil δjk ) Dkl
= − pδij + λδij δkl Dkl + µ(δik δjl Dkl + δil δjk Dkl )
= − pδij + λDkk δij + µ( Dij + D ji )
(35)
= − pδij + λDkk δij + 2µDij
テンソル形式で表示すると
σ = − pI + λ(trD ) I + 2µD
(36)
ここで I は単位テンソル，trD = Dii は D のトレースである．上式を速度で表すと
σ = − pI + λ(trL) I + µ( L + L T )
{
}
= − pI + λ∇ · uI + µ ∇u + (∇u) T
(37)
ここで trL = ∇ · u を用いた．
平均垂直応力 σkk /3 が体積ひずみ速度 Dkk と無関係であると仮定すると
σkk
1
= − p + (3λ + 2µ) Dkk
3
3
(38)
3λ + 2µ = 0
(39)
より
これをストークスの仮説という．これより次式を得る．
2
λ=− µ
3
これを式 (37) に代入すると
{
}
2
σ = − pI − µ∇ · uI + µ ∇u + (∇u) T
3
(40)
(41)
ところで，偏差応力 σ′ は
1
σij′ = σij − σkk δij
3
1
= − pδij + λDkk δij + 2µDij + pδij − (3λ + 2µ) Dkk δij
3
2
= − µDkk δij + 2µDij
3
}
{
2
′
σ = − µ∇ · uI + µ ∇u + (∇u) T
3
(42)
つまり，ストークスの仮説は，せん断応力が偏差応力で表されると仮定するというこ
とである．
6
6
ナビエ・ストークス方程式
運動方程式 (16) に構成方程式 (41) を代入すると，次式が得られる．
(
)
[ {
}]
∂
2
T
(ρu) + ∇ · (ρuu) = −∇ p + ∇ · µ ∇u + (∇u)
−∇
µ∇ · u + b
∂t
3
(43)
上式はニュートン流体の支配方程式であり，ナビエ・ストークス方程式と呼ばれる．
粘性係数 µ が場所に依存しないとすると，∇ · σ の成分は
{ (
)}
(
)
∂σij
∂u j
∂
2 ∂
∂uk
∂ui
∂
=−
( pδij ) −
µ
+
µ
δ +
∂x j
∂x j
3 ∂x j
∂xk ij
∂x j
∂x j
∂xi
(
)
(
)
∂u j
∂p
2 ∂
∂uk
∂
∂ui
=−
− µ
+µ
+
∂xi
3 ∂xi ∂xk
∂x j ∂x j
∂xi
(
)
(
)
∂u j
∂
2 ∂
∂2 u i
∂p
∂uk
− µ
+µ 2 +µ
=−
∂xi
3 ∂xi ∂xk
∂xi ∂x j
∂x j
(
)
∂p
1 ∂
∂uk
∂2 u
=−
+ µ
+ µ 2i
∂xi
3 ∂xi ∂xk
∂x j
これより
1
∇ · σ = −∇ p + µ∇(∇ · u) + µ∇2 u
3
(44)
(45)
これを運動方程式 (16) に代入すると
1
∂
(ρu) + ∇ · (ρuu) = −∇ p + µ∇2 u + µ∇(∇ · u) + b
∂t
3
(46)
非圧縮性流体の場合，連続の式 (9) を式 (43) に適用し，両辺を密度で割ると，次式
が得られる．
[ {
}] 1
∂u
1
+ ∇ · (uu) = − ∇ p + ∇ · ν ∇u + (∇u) T + b
∂t
ρ
ρ
(47)
ここで ν = µ/ρ を動粘性係数という．粘性係数 µ が定数なら，式 (46) より
∂u
1
1
+ ∇ · (uu) = − ∇ p + ν∇2 u + b
∂t
ρ
ρ
参考文献
[1] Y.C. ファン：連続体の力学入門改訂版，培風館 (1980)．
[2] 久田俊明：非線形有限要素法のためのテンソル解析の基礎，丸善 (1992)．
[3] 石原繁：テンソル -科学技術のために-，裳華房 (1991)．
7
(48)

Download Report