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

ナビエ・ストークス方程式
春日 悠
2014 年 7 月 5 日
目次
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
= lim
ϕ( x, t + ∆t)dv −
ϕ( x, t)dv +
ϕ( x, t + ∆t)dv
∆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
Dt
(8)
非圧縮性流体の場合,密度が一定なので,連続の式は次のように表される.
∇·u = 0
2
(9)
運動方程式
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)