201311183 MATSURA - 1

流体の第二粘性が密度ベース
圧縮性流れ数値計算の解に及ぼす影響
THE EFFECT OF SECOND VISCOSITY ON CALCULATION OF
COMPRESSIBLE FLOW BY DENSITY BASED METHOD松浦大志

Taishi MATSURA（指導教員 山本 亨輔）
Abstract – Numerical method for compressible fluid is classified into density-based
method and pressure-based one depending on calculation order of density andpressure. Density-based method calculates density at first. The density error causes
the solution divergence, especially for models of low Mach number flow. This researchevaluates the effect of the second viscosity on numerical solution by calculating several
models of fluid. The models are Couette flow and natural convection in a cavity, andtheir result are compared with previous studies. The comparison suggests that thesecond viscosity affects the stability of calculation.
そのため、密度ベースの解法を低マッハ数の流れに

1 はじめに数値流体力学において、遅い流れ、一般的にはマッ
ハ数が 0.3 以下の流れでは、圧力による流体の密度変化を無視して非圧縮性流体として計算される[1]。また、
熱流体として解く場合、温度による流体の密度変化に対しては、しばしばブシネスク近似[2]が適用される。
これは密度変化が流体の運動に及ぼす影響を、温度差に比例した浮力としてのみ考慮するという近似である。しかし、温度差が大きい流れに対しては適用が難
しく [3]、このような流れは低マッハ数流れでも流体の密度変化を無視できない。そのため、高低マッハ数ど
ちらの流れに対しても適用可能な数値計算アルゴリ
ズムが求められている。
圧縮性流体の数値計算方法は圧力と密度の計算順
序の違いによって密度ベースの解法と圧力ベースの解法に分けられる。圧力を先に求める圧力ベースの解法は非圧縮性流体の解法から研究がすすめられ、
MAC（Marker And Cell）法[4]や SMAC 法[5][6]などがこれに当てはまる。一方、密度ベースの解法では、まず密度と内部エネルギーを支配方程式から計算し、それから圧力を密度と内部エネルギーから求める。しかし、マッハ数が小さい時、微小な密度の変化は非常に大きな圧力の変化を引き起こすため、密度ベースの解法は密度の誤差に対して非常に敏感になり、正しい解が得られないことがある[7]。

対しても適用できるようにするため、様々な研究がなされてきた。Weiss ら[8]は基礎方程式に前処理を施すことで、密度ベースの解法を低マッハ数流れに対しても適用可能にした。しかし、この方法では、基礎方程式の特性を変えてしまうため、非定常問題を解くときは二重時間ステップを用いる必要がある。また、青木ら[9]は保存形で表された基礎方程式系をコロケート格
子上で、有限体積法を用いて離散化することにより、質量保存則を高い精度で満足する計算アルゴリズムを提案した。
上に述べたような既往の研究では、流体の体積粘性をゼロとし、第二粘性係数を消去するストークスの仮説[10]に基づいた基礎方程式を用いている。第二粘性係数は流体の膨張圧縮に対する抵抗性を表すパラメー
タである。通常、ストークスの仮説に基づく流体の数
値計算では、計算の安定性はレイノルズ数をはじめとする様々な無次元数に支配される。しかし、これらの無次元数では第二粘性係数を評価対象としない。したがって、この係数が密度ベースの解法において、計算の安定性を支配するパラメータの一つである場合、その安定条件を満たすような要素や離散化法を用いる必要がある。
そこで本研究では、低マッハ数の圧縮性流れを密度ベースで数値的に解くとき、その数値安定性に対して、流体の第二粘性係数が及ぼす影響を明らかにする。基201311183 MATSURA - 2
礎方程式には連続の式、ナビエ・ストークス方程式、を代入すると、以下のように運動方程式を得る。

内部エネルギー方程式、理想気体の状態方程式を用い𝜌る。形状の自由度や高次精度化を考慮して、空間の離散化には有限要素法を用いる。例題として、定常クエ

𝜕𝜎𝑖𝑗
𝜕𝑢𝑖
𝜕𝑢𝑖
+ 𝜌𝑢𝑗
= 𝜌𝑏𝑖 +
𝜕𝑡
𝜕𝑥𝑗
𝜕𝑥𝑗

（7）

上式の両辺と𝑢𝑖 の内積を求め、𝐾 = (1/2)𝜌𝑢𝑖2 とおい

ット流およびキャビティ内自然対流を取り上げた。そて、整理すると、単位体積当たりの力学的エネルギー

れらの理論解や既往の研究結果と比較を行い、第二粘の保存を表す次式を得る。ただし、整理においては連性係数が流れの数値解に及ぼす影響を検討する。

続の式を用いた。
𝜕𝜎𝑖𝑗 𝑢𝑖
𝜕𝐾 𝜕𝐾𝑢𝑗𝜕𝑢𝑖+= 𝜌𝑢𝑖 𝑏𝑖 +
− 𝜎𝑖𝑗
𝜕𝑡
𝜕𝑥𝑗
𝜕𝑥𝑗
𝜕𝑥𝑗

2 数値解法の概要
2.1 流体の支配方程式
流体の支配方程式は質量、運動量、エネルギーの 3
つの保存則と、理想気体の状態方程式を用いる。いま、空間中に固定された閉曲面𝛤によって囲まれた領域Ωを考える。𝛺内で質量の発生がないとき、質量保存則
は以下のように表される。
𝜕𝜌𝑑𝛺 + ∫ 𝜌𝑢𝑗 𝑛𝑗 𝑑𝛤 = 0
𝛺 𝜕𝑡
𝛤

∫

つぎに、力学的エネルギーに熱力学的エネルギーを加えた単位体積当たりの全エネルギー𝐸(= 𝐾 + 𝜌𝑒)の保存則は以下のように表される。ここで𝑒は単位質量あたりの内部エネルギーである。
𝜕𝐸𝑑𝛺 + ∫ 𝐸𝑢𝑗 𝑛𝑗 𝑑𝛤
𝛺 𝜕𝑡
𝛤

∫

= ∫ 𝜌𝑢𝑗 𝑏𝑗 𝑑𝛺 + ∫ 𝑢𝑖 𝜎𝑖𝑗 𝑛𝑗 𝑑𝛤（1）𝛺

流体の密度、𝑢𝑗 は速度、𝑛𝑗 は𝛤の外向き法線方向単位ベ
クトルの𝑥𝑗 成分である。
（1）式は𝛺内の質量の増減が単位時間当たりに境界面を通して流入、流出する質量と等しいことを表す。
（1）式の左辺第二項に発散定理を適用し、質量保存則が任意の𝛺で成り立つとすると、

𝛤

ここで、𝑞𝑗 はxj 方向の熱流束、𝐽は単位体積あたりの発熱量である。発散定理により面積積分を体積積分で表し、任意の𝛺で保存則が成り立つとすると以下の式が得られる。

（2）

以下の内部エネルギー方程式を得る。
𝜌（3）

つぎに、𝛺内の流体に作用する単位質量あたりの体積力を𝑏𝑖 （𝑖 = 1,2,3）、応力テンソルを𝜎𝑗𝑖 とすると、流体に働く力Fi は
𝐹𝑖 = ∫ 𝜌𝑏𝑖 𝑑𝑉 + ∫ 𝜎𝑗𝑖 𝑛𝑗 𝑑𝑆
𝛺

（4）

𝛤

とかける。発散定理を適用すると、上式は
𝜕𝜎𝑖𝑗
𝐹𝑖 = ∫ (𝜌𝑏𝑖 +
) 𝑑𝛺
𝜕𝑥𝑗
𝛺

𝜕𝑒
𝜕𝑒
𝜕𝑢𝑖 𝜕𝑞𝑗
+ 𝜌𝑢𝑗
= 𝜎𝑖𝑗
−
+𝐽
𝜕𝑡
𝜕𝑥𝑗
𝜕𝑥𝑗 𝜕𝑥𝑗（11）

次に、本研究で扱う流体はニュートン流体と仮定す
ると、応力の構成式は以下のように表される[10]。
𝜎𝑖𝑗 = −𝑝𝛿𝑖𝑗 + 2𝜇𝜖̇𝑖𝑗 + 𝜆𝜖̇𝑘𝑘 𝛿𝑖𝑗（12）

ここで、𝑝は圧力、𝜇は粘性係数、𝜆は第二粘性係数、
𝛿𝑖𝑗 はクロネッカーのデルタ、𝜖̇𝑖𝑗 は次式で定義される
ひずみ速度である。

（5）とかける。このとき、𝛺内の運動量保存則は以下のように表される。
𝜕𝜌𝑢𝑖∫𝑑𝛺 + ∫ 𝜌𝑢𝑖 𝑢𝑗 𝑛𝑗 𝑑𝛤 = 𝐹𝑖
𝜕𝑡
𝛺
𝛤

（10）

式（10）に式（2）、式（8）を代入して整理することで式（2）は、次のように書くこともできる。
𝜕𝑢𝑗
𝜕𝜌
𝜕𝜌
+ 𝑢𝑗
+𝜌
=0
∂𝑡
∂𝑥𝑗
∂𝑥𝑗

𝛺

𝜕𝜎𝑖𝑗 𝑢𝑖 𝜕𝑞𝑗
𝜕𝐸 𝜕𝐸𝑢𝑗+= 𝜌𝑢𝑖 𝑏𝑖 +
−
+𝐽
𝜕𝑡
𝜕𝑥𝑗
𝜕𝑥𝑗
𝜕𝑥𝑗

以下のような連続の式を得る。

（9）𝛤− ∫ 𝑞𝑗 𝑛𝑗 𝑑𝛤 + ∫ 𝐽𝑑𝛺

ここで、𝑡は時間、𝑥𝑗 （𝑗 = 1,2,3）は直交座標系、𝜌は𝜕𝜌 𝜕𝜌𝑢𝑗
+
=0
∂𝑡
∂𝑥𝑗

（8）

（6）

上式の左辺第二項に発散定理を適用し、式（2）、式（5）1 𝜕𝑢𝑖 𝜕𝑢𝑗𝜖̇𝑖𝑗 = (+)
2 𝜕𝑥𝑗 𝜕𝑥𝑖

（13）

式（12）において主成分の平均値を求めると、
𝜎𝑘𝑘
2= −𝑝 + (𝜆 + 𝜇) 𝜖̇𝑘𝑘
3
3

（14）

となり、このとき𝜖̇𝑘𝑘 の係数(𝜆 + 2/3 𝜇)を体積粘性係数と呼ぶ。この体積粘性をゼロとして、平均主応力（応201311183 MATSURA - 3
力の主成分の平均値）が圧力に等しいとする仮定をス
トークスの仮説という。この仮定を導入すると、
2
𝜆=− 𝜇
3

（15）

𝑢𝑖𝑛+1 − 𝑢𝑖𝑛𝜕𝑢𝑖𝑛+1
𝜌𝑛 (
+ 𝑢̅𝑗
)
Δ𝑡
𝜕𝑥𝑗
= 𝜌𝑛 𝑏𝑖 −

𝑛+1
𝑛+1
𝜕𝜖̇𝑖𝑗
𝜕𝑝𝑛
𝜕𝜖̇𝑘𝑘
+ 2𝜇
+𝜆
𝜕𝑥𝑖
𝜕𝑥𝑗
𝜕𝑥𝑖

（22）

となって第二粘性係数は消去される。しかし、第二粘ここで、左辺の移流速度は𝑢𝑗𝑛+1 が未知であるため、初性係数は流体の膨張圧縮に対する抵抗性を表す係数

期条件として 2 次精度の Adams-Bashforth 法[13]を用である。密度の誤差が発散要因になりやすい密度ベー

いて以下のように近似することにより線形化していスの解法においては、第二粘性係数が無視できない可る。

能性がある。したがってその影響を評価することで数

31𝑢̅𝑗 = 𝑢𝑗𝑛 − 𝑢𝑗𝑛−122

値計算の安定性向上が期待できる。そのため本研究で（23）は、第二粘性係数の影響を評価するため、ストークス

式（22）では、密度にnステップ目の値を用いるこの仮説は導入せず、𝜇と𝜆は独立な定数として扱う。

とで𝑢𝑛+1 を求めることができるようになっている。次また、熱流束𝑞𝑗 がフーリエの法則[11]に従うとすると、に、式（22）で求めた𝑢𝑛+1 を用いて式（20）を時間方向に式（24）のように陰的に離散化し、𝑇 n+1 を求める。

熱流束は以下のように表される。
𝜕𝑇
𝑞𝑗 = −𝜅𝜕𝑥𝑗（16）

ここで、𝑇は絶対温度、𝜅は熱伝導率である。
流体に理想気体を仮定すると、内部エネルギーと状態方程式はそれぞれ以下のように表される[12]。
𝑒 = 𝑐𝑣 𝑇（17）

𝑝 = 𝜌𝑅𝑇（18）

ここで、𝑐𝑣 は定積比熱、𝑅は気体定数である。
よって、式（3）、式（7）、式（11）、式（12）、式（13）、式（16）、式（17）より、以下のような式を得る。
𝜕𝑢𝑖
𝜕𝑢𝑖
𝜌(
+ 𝑢𝑗
)
𝜕𝑡
𝜕𝑥𝑗
𝜕𝜖̇𝑖𝑗
𝜕𝑝
𝜕𝜖̇𝑘𝑘
= 𝜌𝑏𝑖 −
+ 2𝜇
+𝜆
𝜕𝑥𝑖
𝜕𝑥𝑗
𝜕𝑥𝑖

𝜌𝑛 𝑐𝑣 (

𝑇 𝑛+1 − 𝑇 𝑛𝜕𝑇 𝑛+1+ 𝑢𝑗𝑛+1
)
Δ𝑡
𝜕𝑥𝑗

=𝜅

𝜕 2 𝑇 𝑛+1
− 𝑝𝑛 𝜖̇𝑘𝑘 + 𝐽 + 𝛷 𝑛+1
𝜕𝑥𝑗2

（24）

次に式（3）において密度を陰的に扱い、時間方向に式（25）のように離散化する。
𝜕𝑢𝑗𝑛+1
𝜌𝑛+1 − 𝜌𝑛𝜕𝜌𝑛+1+ 𝜌𝑛+1+ 𝑢𝑗𝑛+1
（25）
Δ𝑡
𝜕𝑥𝑗
𝜕𝑥𝑗
=0
次に、重みつき残差法の考え方に基づき式（22）、式（24）、式（25）に重み関数をかけて𝛺内で積分する。
運動方程式の圧力項と粘性項、内部エネルギー方程式の熱伝導項に部分積分を施すことで重みつき残差式を得る。得られた重みつき残差式を双一次要素を用い（19）

た Galarkin 法に基づく有限要素法[13]によって離散化する。状態方程式は圧力を求めるために用いる。以下

𝜕𝑇
𝜕𝑇
𝜕2𝑇
𝜌c𝑣 ( + 𝑢𝑗
) = 𝜅 2 − 𝑝𝜖̇𝑘𝑘 + 𝐽
𝜕𝑡
𝜕𝑥𝑗
𝜕𝑥𝑗
+𝛷

では簡単のため、2 次元問題について考え、𝑥1 = 𝑥、
（20）x2 = 𝑦とする。{𝑵}を形状関数。{𝒖}、{𝒗}をそれぞれ𝑥方向、𝑦方向の節点における速度を表すベクトルとし、各係数行列を定義すると、以下の離散化式を得る。

ただし、Galarkin 法に基づく有限要素法を用いたため、重み22

1 𝜕𝑢𝑖 𝜕𝑢𝑗𝜕𝑢𝑘
𝛷 = 𝜇(+) +𝜆()
2
𝜕𝑥𝑗 𝜕𝑥𝑖
𝜕𝑥𝑘

（21）

である。
本研究では式（3）、式（18）、式（19）、式（20）を基礎方程式とし、未知量𝑢𝑖 、𝑝、𝜌、𝑇を求める。

2.2 方程式の離散化
いま、n ステップ目の物理量𝜌𝑛 、𝑢𝑖𝑛 、𝑇 𝑛 、𝑝𝑛 が既知の時、式（19）において速度を陰的に扱い、時間方向に次のように離散化する。

関数と未知変数の形状関数は同一のものとなっている。以上をまとめた式（26）、式（27）、式（28）の連
立一次方程式を解くことで、n + 1ステップにおける物理量を得ることができる。ただし、総和記号は要素係数行列の重ねあわせを表す。数値積分にはガウス-ル
ジャンドル法[14]を用い、積分点は 5 とした。連立一次方程式の解法には GMRES 法[15]を用いた。

201311183 MATSURA - 4
𝑛+1[𝐀NS ] {𝒖𝑛+1 } = {𝒃𝑁𝑆 }𝒗（26）[𝐀EN ]{𝑻𝑛+1 } = {𝒃𝐸𝑁 }（27）[𝐀C ]{𝝆𝑛+1 } = {𝒃𝐶 }（28）{𝒑𝒏+𝟏 } = 𝑅{𝝆𝑛+1 } ∘ {𝑻𝑛+1 }（29）

1[𝐀𝐍𝐒 ] = ∑ ( [𝐌U ] + [𝐂U ] − [𝐒])
Δ𝑡

（30）

1
{𝒃𝑁𝑆 } = ∑ ( [𝐌U ] − [𝐆]{𝒑𝑛 } + {𝒇})
Δ𝑡

（31）

1
[𝐀𝐄𝐍 ] = ∑ ( [𝐌𝐄 ] + [𝐂𝐄 ] − [𝐇𝐄 ])
Δ𝑡

（32）

1
[𝐌 ] − [𝐃E ]{𝒑𝑛 } + {𝑾}
Δ𝑡 𝐄

（33）

𝑒

𝜕{𝑵}
{𝑵}𝑇
𝜕𝑥
[𝐆] = ∫
𝑑𝛺
𝜕{𝑵}
𝛺
{𝑵}𝑇
[ 𝜕𝑦
]

（43）

[𝐌E ] = 𝑐𝑣 ∫ {𝑵}{𝑵}𝑇 {𝝆𝑛 }{𝑁}𝑇 𝑑Ω

（44）

𝛺𝑒

[𝐂E ]
= 𝑐𝑣 ∫ {𝑵}{𝑵}𝑇 {𝝆𝑛 }{𝑵}𝑇 [𝐜𝐨𝐧𝐯 𝑛+1 ]𝑑𝛺𝛺𝑒𝑒

𝜕{𝑵} 𝜕{𝑵}𝑇[𝐇E ] = 𝜅 ∫ (
𝜕𝑥 𝜕𝑥
𝛺𝑒

（46）

[𝐃E ] = ∫ {𝑵}𝑑𝑖𝑣 𝑛+1 {𝑵}𝑇 𝑑𝛺

（47）

𝜕{𝑵} 𝜕{𝑵}𝑇
+
) 𝑑𝛺
𝜕𝑦 𝜕𝑦

𝑒

{𝑏𝐸𝑁 } =[𝐀C ] = ∑ (𝑒1[𝐌 ] + [𝐃C ] + [𝐂C ])Δ𝑡 C

1{𝒃𝐶 } = ∑ ( [𝐌𝐂 ]{𝝆𝑛 })
Δ𝑡

Ωe

（34）

[𝐌U ] = [[𝐌] = ∫

𝑑𝑖𝑣 𝑛 =

𝜕{𝑵}𝑇 𝑛+1
𝜕{𝑵}𝑇 𝑛+1{𝒖 } +
{𝒗 }
𝜕𝑥
𝜕𝑦

（48）

𝜕{𝑵}𝑇
𝜕{𝑵}𝑇
+ {𝒗𝑛 }
𝜕𝑥
𝜕𝑦

（49）

（35）

𝑒

[𝐌]

[𝐜𝐨𝐧𝐯 𝑛 ] = {𝒖𝑛 }

]
[𝐌]

（36）

{𝑵}{𝑵}𝑇 {𝝆𝑛 }{𝑵}𝑇 𝑑Ω

（45）

（37）

{𝑾}
= ∫ {𝑵}{𝒖𝑛+1 }𝑇 [𝐁]𝑇 [𝐃][𝐁]{𝒖𝑛+1 }𝑑𝛺

（50）

𝛺𝑒

Ω𝑒

[𝐂̅]
[𝐂U ] = [

]
[𝐂̅]

（38）

[𝐌C ] = ∫ {𝑵}{𝑵}𝑇 𝑑Ω

（51）[𝐃C ] = ∫ {𝑵}𝑑𝑖𝑣 𝑛 𝑑Ω（52）

[𝐂C ] = ∫ {𝑵}{𝑵}𝑇 [𝐜𝐨𝐧𝐯 𝑛+1 ]𝑑Ω

（53）

Ωe

Ωe

[𝐂̅] = ∫ {𝑵}{𝑵}𝑇 {𝝆𝑛 }{𝑵}𝑇 [𝐜𝐨𝐧𝐯̅̅̅̅̅̅̅]𝑑Ω

（39）

Ωe

Ωe

[𝐒] = ∫ [𝑩]𝑇 [𝑫][𝑩]𝑑𝛺（40）

本手法のフローチャートを図 1 に示す。

𝛺𝑒

𝜕{𝑵}𝑇
𝜕𝑥
[𝐁] =
𝑇

𝜕{𝑵}
[ 𝜕𝑦

𝜆 + 2𝜇[𝐃] = [ 2𝜇

𝜕{𝑵}𝑇
𝜕𝑦
𝜕{𝑵}𝑇
𝜕𝑥 ]

（41）

2𝜇
𝜆 + 2𝜇]（42）2𝜇図 1 計算のフローチャート

201311183 MATSURA - 5

表 1 モデルパラメータ
計算領域[m] × [m]1×1

時間刻みΔt[s]0.01

計算ステップ数100

気体定数𝑅[𝐽/𝑘𝑔 ∙ 𝐾]

1000

定積比熱𝑐𝑣 [𝐽/𝑘𝑔 ∙ 𝐾]

1000

熱伝導率𝜅[𝑊/𝑚]0.1

初期温度𝑇0 [𝐾]100

初期圧力𝑃0 [𝑃𝑎]105

初期密度𝜌0 [𝑘𝑔/𝑚3 ]1

図 2 モデル概要図
ケース
粘性係数𝜇[𝑃𝑎 ∙ 𝑠]
移動壁速度𝑈[𝑚/𝑠]
レイノルズ数 𝑅𝑒
マッハ数 𝑀

表 2 各ケースのパラメータ
1
2
3
1
0.1
10
0.1
1
1
0.1
10
0.1
2.24×10-4 2.24×10-3 2.24×10-4

4
1
10
10
2.24×10-3

この問題は非圧縮流体の例題であるため、本研究で3 数値計算例

はマッハ数が十分小さいとみなせる範囲で計算を行った。マッハ数は音速を𝑎として理想気体の場合、次式

3.1 定常クエット流れで定義される[1]。

3.1.1 モデル概要

𝑀=

クエット流れは、図 2 のように二枚の平板の間に流体が満たされているとき、一方の平板を固定し、他方の平板を水平方向に移動させると、流体の粘性によって速度分布が発生するという問題である。圧縮性が無視でき、水平方向に圧力勾配が無いとき、その速度分
布は直線分布になることが知られている[16]。非圧縮流体を仮定すると、式（2）において𝜌 = 𝑐𝑜𝑛𝑠𝑡.となり、（𝜕𝑢𝑘 ）/（𝜕𝑥𝑘 ） = 𝜖̇𝑘𝑘 = 0となる。このとき、式（12）の応力の構成式で、体積ひずみ速度𝜖̇𝑘𝑘 がゼロである
とき、流体の運動は第二粘性係数とは無関係に定まる。
よって、マッハ数が十分小さく、流れが非圧縮とみなせるとき、数値解は𝜆によらず、理論解に近づくと考えられる。
モデルのパラメータを表 1 に示す。粘性係数、移動壁速度を変えて表 2 に示す 4 ケースで計算を行った。
これらのケースは、ケース 1 と 3、ケース 2 と 4 は以下の式で定義されるレイノルズ数𝑅𝑒が一致している
ため、定性的に同様の挙動を示すと考えられる。
𝜌𝑈𝐻
𝑅𝑒 =
𝜇

（54）

ここで、代表長さは平板間距離を、代表密度は初期条件における密度を用いた。

𝑈𝑈=𝑎 √𝛾𝑅𝑇（55）

𝛾は比熱比であり、理想気体の場合𝛾 = 1 + 𝑅/𝑐𝑣 で表される[12]。
初期条件として理論解の流速を与え、20 ステップ計算を行った。流速の初期値からの変化率を指標として第二粘性の影響を検討した。

3.1.2 結果と考察
ケース 1 から 4 に対してそれぞれ要素数 400 と 100の場合について計算を行った。それぞれの場合における 20 ステップ後の計算領域内における流速の変化率の絶対値の最大値を図 3、図 4 に示す。初期値として理論解の速度分布を与えているため、計算を進める過程で速度が初期値から変化した場合、その変化分を誤差として評価した。計算したすべての場合において、第二粘性係数が大きいとき、流速の初期値からの変化
率が小さくなる傾向が見られた。また、図 3（a）と図3（b）、図 4（a）と図 4（b）を比較すると、レイノルズ数が異なる場合でも、第二粘性係数が大きい程誤差が小さいという傾向は共通して確認できた。一般に、レイノルズ数が 2300 より大きいとき、流れは乱流と201311183 MATSURA - 6（a）𝑅𝑒 =0.1

（a） 𝑅𝑒 =0.1（b） 𝑅𝑒 =1（b） 𝑅𝑒 =1

図 3 流速の変化率（400 要素）

図 4 流速の変化率（100 要素）

なり[16]、計算の不安定化の原因となる。今回のモデル

渉し合って渦が発生する。このモデルは熱対流現象のでは、乱流となるレイノルズ数に対して十分小さいレ

数値計算におけるベンチマーク問題として知られてイノルズ数で計算を行った。このことから、十分に小おり、多くの計算例][17][18][19][20][21]がある。

さいレイノルズ数の流れに対しても、安定性評価の基準としての、第二粘性係数の有用性が示唆される。ケ

3.2.2 第二粘性係数の影響

ース 1 とケース 3、ケース 2 とケース 4 はレイノルズ

第二粘性係数の影響を評価するために、次のような数が等しく、ともに第二粘性係数が大きくなる程誤差計算を行った。まず、図 5 と同じ計算領域において、が小さくなった。このことから、数値計算の安定性は左右の壁に温度差を与えず、𝑦 = 𝐿で基準圧力と等し第二粘性係数の影響を受けると考えられる。また、今くなる静圧分布となるように圧力と密度の初期値を回用いたモデルでは、節点数の違いによる差は第二粘与える。このとき、流体に働く力は領域内の任意の位性係数による差に比べて小さく、クエット流れのよう置で釣り合っており、時刻を進めても流体は初期条件

な単純な流れのモデルに対しては、少ない節点数でもを保ったまま静止し続けるはずである。図 6 のように妥当な解を得られることがわかった。

領域を分割し、同一時間ステップでの反復計算を行わずに 100 ステップ計算した。初期値からの変化率を指3.2 キャビティ内自然対流

標として第二粘性係数の影響を考える。流体の物性値

3.2.1 モデル概要

は空気の値を参照し、第二粘性係数はλ =10k[Pa ∙ s]

つぎに、二次元正方形容器内に発生する自然対流を（1≤ k ≤10）として計算を行った。表 2 にモデルのパ数値的に再現し、既存の研究と比較する。図 5 に示す

ラメータを示す。図 7、図 8、図 9、図 10 に計算領域

ような、一辺が𝐿 =0.04[m]の正方形容器内に流体が満中央における、圧力、温度、密度の初期値からの変化

たされている。鉛直下向きに重力が作用しているとき、率の絶対値の推移と、100 ステップ後における質量の左右の壁面に温度差を与えると温度、密度、圧力が干初期値からの変化率を示す。
図 7 から図 9 より、
λ ≤103

201311183 MATSURA - 7
表 3 モデルパラメータ

図 5 モデル概要図

計算領域[m] × [m]0.04×0.04

時間刻みΔ𝑡[s]0.1気体定数𝑅[𝐽/𝑘𝑔 ∙ 𝐾]

286.8

定積比熱𝑐𝑣 [𝐽/𝑘𝑔 ∙ 𝐾]

717

熱伝導率κ[𝑊/𝑚]0.023

粘性係数𝜇[𝑃𝑎 ∙ 𝑠]

1.60×10-5

第二粘性係数λ[𝑃𝑎 ∙ 𝑠]

10k（1≤ 𝑘 ≤10）

初期温度𝑇0 [𝐾]283.15

基準圧力𝑃0 [𝑃𝑎]1.01325×105

基準密度𝜌0 [𝑘𝑔/𝑚3 ]1.2477ではステップが進むと急激に変化率が大きくなり、解1.465[K]の温度差を与えた。この値は、次式で定義さが発散していることがわかる。図 7 から図 9 より、れるレイリー数𝑅𝑎の値が既往の研究[9]で用いられて

λ ≤103 ではステップが進むと急激に変化率が大きくいる1.0 × 104となるように設定した。

なり、解が発散していることがわかる。特に図 7、図𝑅𝑎 =8 と図 9 を比べると、密度の誤差に対して圧力、温度

𝑔𝛽𝛥𝑇𝐿3𝛼𝜈（56）の誤差が大きい。この点は、密度ベースの解法が低マ

ただし、𝛼(= 𝜅/𝜌𝑐)は温度拡散係数、𝛽は体積膨張率、ッハ数流れにおいて、密度の誤差に敏感であるという

𝜈(= 𝜇/𝜌)は動粘性係数であり、初期条件での値を参照

肖らによる文献と一致する。このとき、図 10 より、した。

領域内の質量が減少し、質量保存則が満たされていな
いことがわかる。一方、第二粘性係数が大きいとき、

3.2.4 結果と考察

計算結果はステップ数が増加しても、解は初期条件を

図 12 に本手法による𝑡 =120[𝑠]での水平方向速度、高い精度で保っていることが確認できる。このことか

鉛直方向速度、温度の分布を示す。等高線の間隔はそら、第二粘性係数が計算の安定性に影響を及ぼすことれぞれ 3.190×10-9[𝑚/𝑠]、3.028×10-9[𝑚/𝑠]、0.1465[𝐾]がわかる。次に、図 11 に、計算を始めた最初のステッである。また、既往の研究として、青木らによる結果

プにおける、領域内の体積ひずみ速度の最大値を示す。

[9]を示す。青木らの結果は流れが定常状態に至った後

計算の安定性に対する考察として、式（22）の内部エの分布を表し、等高線の間隔はそれぞれ 1.686×10-

ネルギー方程式の右辺第二項𝑝𝜖̇𝑘𝑘 について考える。圧3[𝑚/𝑠]、2.105×10-3[𝑚/𝑠]、0.1465[𝐾]である。容器内に力がおよそ大気圧に等しいとき、𝑝 ≅105[kg/m3 ]、渦が発生している点で定性的な一致は確認できる。し𝜌 ≅100[𝑃𝑎]であり、密度に対して圧力のオーダーが大かし、速度の大きさや温度分布にはおおきな差がある。

きい。このとき、体積ひずみ速度𝜖̇𝑘𝑘 にわずかな誤差が

本研究の手法では、移流項に対する安定化処理を施し生じると、圧力との乗算によって誤差の影響が大きくていないため、速度分布に数値振動が見られる。また、なり、計算が発散すると考えられる。これに対して、初期状態を静止状態とし、非定常問題として解いてい式（12）の応力の構成式より、第二粘性係数が大きい

るため、𝑡 =120[s]では系が定常状態に至っていない可

とき、すなわち流体の膨張圧縮に対する抵抗が大きい

能性がある。また、渦が形成された後の流れの発達にとき、図 11 のように体積ひずみ速度は小さく、計算

第二粘性係数が影響していることがわかる。

が発散しなかったと考えられる。

最後に、図 13 に本手法による結果における質量の誤差を、図 14 に青木らによる結果を示す。計算を行

3.2.3 自然対流の再現った時間において誤差は十分小さく、質量保存則を高

前節の結果より、第二粘性係数𝜆 = 105 [pa ∙ s]を用い精度で満足していると言える。ただし、本手法によいて、自然対流によりキャビティ内に発生する渦を再って得られた数値解は速度の大きさが青木らの結果現する。物性値は前節と同じ値を使用し、要素数をと異なるため、単純に比較することは困難であると考32 × 32 要 素 と し た 。 左 右 の 壁 面 に 温 度 差 ∆𝑇 =えられる。今後の課題としては長時間の計算と、発達201311183 MATSURA - 8

0.04

0.03

0.02

0.01

0
0

0.01

0.02

0.03

0.04

図 6 メッシュ分割図（20×20）

図 7 圧力の初期値からの変化率

図 8 温度の初期値からの変化率

図 9 密度の初期値からの変化率

図 10 質量の初期値からの変化率

図 11 体積ひずみ速度の最大値

した流れにおける第二粘性係数の影響の検討が挙げう際に、適切な解を得るためには、流体の第二粘性をられる。

考慮して安定性を評価し、適切な要素サイズや時間増分などを設定する必要がある可能性が示唆される。

3.2.5 まとめ
クエット流れ、キャビティ内自然対流の計算結果から、第二粘性係数が密度ベースの圧縮性流体数値計算の安定性に影響を与えることがわかった。このことか
ら、大規模化かつ高精度化が要求される数値計算を行

4 結論
本研究では、密度ベースの圧縮性流れ数値計算の解に対して流体の第二粘性係数が及ぼす影響を調べる
ため、有限要素法による数値計算コードを作成し、定201311183 MATSURA - 9（a）水平速度分布

図 13 質量誤差（本手法）（b）鉛直速度分布（c）温度分布本手法（𝑡 =120[𝑠]）青木ら（定常状態）
図 12 本手法と青木らによる結果の比較
常クエット流れとキャビティ内自然対流を例題とし

図 14 質量誤差（青木ら）[5]

F. Harlow and A. Amsden, Numerical calculation ofて数値計算を行った。計算結果から、第二粘性係数はalmost incompressible flows, J. Comput. Phys., 3,

計算の安定性に影響を及ぼし、十分に遅い流れでは、

80 （1968）

第二粘性係数がある値よりも小さい時、計算が不安定[6]になった。今後の課題として、第二粘性係数が流れのdynamics calculation method for all flow speeds, J.

時間発展に及ぼす影響の評価と、第二粘性係数に関する安定条件の解明が挙げられる。

F. Harlow and A. Amsden, A Numerical fluid
Comput. Phys., 8, 197[7]

肖 鋒, 伊井 仁志, 小野寺 直幸, 計算流体
力学-CIP マルチモーメント法による手法-, コ

参考文献
[1]

松尾

ロナ社, 2009
一 泰 , 圧 縮 性 流体 力学 , オ ー ム 社 ,[8]

Preconditioning Applied to Variable and Constant

1994
[2]
[3]

Helmi MLAOUAH, 辻 俊博, 長野 靖尚, 温
度差の大きい閉空間における熱対流, 日本機械学会論文集（B 編）, vol.62（1996）[4]

Density Flows, AIAA J., 22, 745（1995）

J. S. Turner: Buoyancy Effects in Fluids,（CambridgeUniversity Press, 1973） 9–11

F. Harlow and E. Welch, Numerical calcuration oftime-dependent viscous incompressible flow of
fluid with a free surface, Physics of Fluids, 8, 2182
（1965）

Jonathan M. Weiss and Wayne A. Smith,

[9]

青木

一真, 牛島省, 鳥生

大祐, 圧縮性

流体に関する保存形陰的計算アルゴリズムの提案, 土木学会論文集 A2（応用力学）, Vol. 69,
No. 2, I_51-I_58, 2013[10] 清水昭比古, 連続体力学の話法-流体力学

材料力学の前に-, 森北出版, 2012[11] 中山 司, 流れ解析のための有限要素法入門,
東京大学出版会, 2008[12] 山下 博史, 熱力学, 培風館, 2014

201311183 MATSURA - 10[13] 日本計算工学会流れの有限要素法研究委員会,
続・有限要素法による流れのシミュレーショ
ン, 丸善出版, 2012
善出版,1998
清次, 張

pp.1045-1052[19] 角田 和彦, 川原 祐子, 登坂 宣好, 正方形[14] 日本機械学会, 計算力学ハンドブックⅠ, 丸[15] 藤野

対流, 日本機械学会論文集（B 編）,vol.66（2000）,

キャビティ内自然対流の指数関数型 PetrovGalerkin 有限要素解析, 日本機械学会論文集

紹良, 反復法の数理, 麻倉

書店, 1996（B 編）, vol.59（1993）, pp.2463-2468[20] 姫野武洋, 渡辺紀徳, 低重力環境におけ[16] 日本機械学会, 流体力学, 丸善出版, 2005

る熱流体管理に関する研究, 日本機械学会論[17] Davis, G., Int. J. Numer. Methods Fluids, vol.3

文集（B 編）, vol.69（2003）, pp.266-273（1983）, pp.249-264[18] 柴田 悦太郎, 槙原[21] 白石孝文, 棚橋靖幸, 加藤信介, 石田義洋, 村上孝彦, CIP

周三, 大岡 龍三, サーマルキャビティ内にお有限要素法を用いた正方形キャビティ内自然ける圧縮性自然対流解析, 空気調和・衛生工学会学術講演会講演論文集（1998）, pp.1417-1420

謝辞
本研究を進めるに当たり、日ごろの研究から論文の執筆にわたり、非常に丁寧にご指導いただきました山本
亨輔先生に深謝いたします。また、ゼミのたびに貴重なご指摘、アドバイスをいただきました松島亘志先生、
高橋悠太先輩をはじめとする研究室の先輩同期の皆さん、そして大学での研究生活を支えてくれた家族のおかげでこのような形で論文を執筆できたことに対して、心から感謝いたします。