計算熱流体講義ノート (初版)

埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
計算熱流体講義ノート
(初版)
埼玉工業大学工学部機械工学科
小西克享
1
1/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
2/66
はじめに
本書は，工学部の機械系の大学院に所属する学部生を対象として，半期 15 回程度の授業で，
物性値の推定方法，特性曲線法による管路内流体解析，エネルギ保存式の数値解法などを教育す
るために必要な事項を講義ノートとしてまとめたものである．特に特性曲線法に関しては，その
基礎から油圧回路に対する適用例を詳述するとともに，サンプルプログラムを記載しており，初
学者にとって貴重な資料となるものと期待している．
内容は今後とも加筆修正の予定である．内容に関して不明な点やお気付きの点があれば著者ま
でご連絡いただきたい．
平成 21 年 12 月 1 日
埼玉工業大学 工学部 機械工学科 小西克享
2
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
3/66
目次
第1章
文献
･･････････････
5
第2章
物性値の推定式
2.1.1 密度
･･････････････
6
2.1.2 粘度
･･････････････
6
2.1.3 比熱
･･････････････
6
2.1.4 蒸気圧
･･････････････
6
2.1.5 沸点
･･････････････
7
2.1.6 蒸発潜熱
･･････････････
8
2.2.1 密度
･･････････････
8
2.2.2 粘度
･･････････････
8
2.2.3 比熱
･･････････････
9
2.3.1 燃焼ガス成分のモル定圧比熱
･･････････････
9
2.3.2 燃焼ガス成分のモルエンタルピ
･･････････････
9
･･････････････
10
2.1 水
2.2 空気
2.3 燃焼ガスの熱力学的性質
2.4 石油
2.5 補間によるデータの算出
2.5.1 補間法
･･････････････ 11
2.5.2 水の動粘性係数
･･････････････ 11
第3章
特性曲線法による油圧管路動特性解析
3.1 仮定
･･････････････
15
3.2 運動方程式と連続の式
･･････････････ 15
3.3 特性曲線法による変換
･･････････････ 15
3.4 数値計算法
･･････････････ 17
3.5 油圧要素
3.5.1 基礎式
･･････････････ 18
3.5.2 閉塞端
･･････････････ 18
3.5.3 タンク
･･････････････ 19
3.5.4 ポンプ
･･････････････ 19
3.5.5 オリフィス
･･････････････ 20
3.6
Water hammering （水撃現象）の計算
3.6.1 基礎式
･･････････････ 21
3.6.2 初期条件
･･････････････ 22
3.6.3 境界条件
･･････････････ 22
3.6.4 物性値
･･････････････ 22
3.6.5 計算準備
･･････････････ 23
3.6.6 プログラム例“Water hammering”
･･････････････ 23
3
埼玉工業大学（小西克享）
第4章
計算熱流体講義ノート（初版）
常微分方程式の数値解法＝ルンゲクッタ法
4.1
1 階常微分方程式
･･････････････ 26
4.2
2 階常微分方程式
･･････････････ 26
4.3 プランジャポンプとノズルの計算
4.3.1 基礎式
･･････････････ 27
4.3.2 初期条件
･･････････････ 28
4.3.3 境界条件
･･････････････ 28
4.3.4 物性値
･･････････････ 29
4.3.5 計算準備
･･････････････ 29
4.3.6 プログラム例“Nozzle”
･･････････････ 30
4.4 プランジャポンプ＆油圧ジャッキの計算
4.4.1 基礎式
･･････････････ 34
4.4.2 初期条件
･･････････････ 35
4.4.3 境界条件
･･････････････ 35
4.4.4 物性値
･･････････････ 35
4.4.5 計算準備
･･････････････ 36
4.4.6 プログラム例“Cylinder”
･･････････････ 36
4.5 プランジャポンプ＆燃料弁の計算
4.5.1 基礎式
･･････････････ 42
4.5.2 初期条件
･･････････････ 43
4.5.3 境界条件
･･････････････ 43
4.5.4 物性値
･･････････････ 44
4.5.5 計算準備
･･････････････ 44
4.5.6 プログラム例“Fuel Valve”
･･････････････ 44
第5章
１次元圧縮性流れ
5.1 質点の運動
･･････････････ 50
5.2 剛体の運動
･･････････････ 50
5.3 ガス体の運動
5.3.1 流速に関する式
･･････････････ 50
5.3.2 圧力に関する式
･･････････････ 51
第6章
熱の流れ
6.1 偏微分方程式
･･････････････ 58
6.2 差分近似
6.2.1 導関数の近似法
･･････････････ 58
6.2.2 陽解法
･･････････････ 59
6.2.3
Crank-Nicolson の陰解法
･･････････････ 59
6.2.4 偏微分方程式の無次元化
･･････････････ 61
第7章
数値計算上の注意事項
･･････････････
4
12
4/66
埼玉工業大学（小西克享）
第1章
計算熱流体講義ノート（初版）
第1章
文献
5/66
文献
熱流体に関する数値計算には任意の温度圧力における物性値が必要
方法
１．物性値に関する推定式を利用
推定式の適用できる温度，圧力の範囲を無視しないこと
最大誤差に注意すること．
２．推定式が見あたらない場合は実測値より推定式を作成
最小２乗法
補間多項式を利用
問題点
温度の対する値は計算可．多くの場合，圧力に対する値が不明
熱物性値に関する文献
物性値データの記述がある文献
日本熱物性学会編「熱物性ハンドブック」養賢堂，¥9,270
日本熱物性学会編「熱物性資料集
断熱材編」養賢堂，¥4,635
日本化学会編「改訂３版
化学便覧
日本機械学会「技術資料
流体の熱物性値集」丸善，¥48,410
日本機械学会「伝熱工学資料
基礎編」丸善，¥20,085
改訂第４版」丸善，¥12,360
日本機械学会「燃焼のデータベース
DACOS」¥60,000
Chase, M. W. 他５名 JANAF Thermochemical Tables, 3rd Edition, (1985), NIST
物性値推算の記述がある文献
東京天文台編纂「理科年表」丸善
佐藤一雄「物性定数推算法
第８版」丸善，¥2,900
化学工学協会編「改訂５版
化学工学便覧」丸善
Reid, R. C. Prausnitz, J. M. et al.: The Properties of Gases and Liquids, 3rd ed., McGRAW-HILL(1977)
平田光穂監訳「気体・液体の物性値推算ハンドブック」マグロウヒル(1985)
Warren J. Lyman 他２名, Handbook of Chemical Property Estimation Methods, American Chemical
Society. ISBN0-8412-1761-0
Robert C, Weast, CRC Handbook of Chemistry and Physics, 68th Edition, CRC Press. Inc.
ISBN 0-8493-0468-7
John A. Dean, LANGE’S HANDBOOK OF CHEMISTRY, 11th ed. McGRAW-HILL BOOK COMPANY
水谷幸夫「燃焼工学 第２版」森北出版
西山哲男「流体力学(Ⅰ)」日刊工業新聞社，¥1,400
千輝淳二「伝熱計算法」工学図書，¥3,800
5
埼玉工業大学（小西克享）
第2章
計算熱流体講義ノート（初版）
第2章
物性値の推定式
物性値の推定式
2.1 水
2.1.1
密度
推定式 1
ρ = (999.83952 + 16.945176t − 7.987041 × 10 −3 t 2 − 46.170461 × 10 −6 t 3
+ 105.56302 × 10 − 9 t 4 − 280.54253 × 10 −12 t 5 ) /(1 + 16.879850 × 10 − 3 t )
ρ ；1 atm における密度, kg/m3
t；温度（C）
文献：CRC Handbook of Chemistry and Physics
推定式 2
ρt =
ρ4
1+ α t − 4
（g/cm3）
ρ 4 ；4C における密度=0.99997g/cm3
α ；膨張係数=0.00043
t；温度（C）
文献
2.1.2
西山哲男「流体力学(Ⅰ)」
粘度
μ=
文献
2.1.3
183
. × 10 −4
1 + 0.036t + 0.000185t 2
kgw ⋅ s / m 2
西山哲男「流体力学(Ⅰ)」日刊工業新聞社
比熱
推定式 1（文献：The Properties of Gases and Liquids, 3rd ed.）
C p = A + BT + CT 2 + DT 3
A = 7.701, B = 4.595 × 10 − 4 , C = 2.521 × 10 − 6 , D = −0.859 × 10 − 9
C p ；cal/gram-mole/K
T ；K
推定式 2（文献：
佐藤一雄「物性定数推算法
C p = A + BT + CT 2 ,
第 8 版」）
A = 7.256, B = 2.298 × 10 −3 , C = 0.283 × 10 −6
C p ；cal/gram-mole/K
T ；K (298～1500K)
2.1.4
蒸気圧
6
6/66
埼玉工業大学（小西克享）
①
計算熱流体講義ノート（初版）
第2章
物性値の推定式
Antoine（アントワン）の式
例1
（文献：改訂 5 版
化学工学便覧）
B
,
C +θ
log10 Psat [kPa ] = A -
A = 7.07406,
B = 1657.46, C = 227.02
θ ；温度，C (適用範囲 10～168C）
例2
（文献：LANGE’S HANDBOOK OF CHEMISTRY）
B
C+t
.
, B = 1750.286, C = 235.0
t = 0 − 60C
A = 810765
t = 60 − 150C A = 7.96681, B = 1668.21, C = 228.0
log 10 P[ mmHg] = A -
②
Wexler-Hywood の式
（文献：伝熱工学資料（改訂第 4 版）p.355）
温度 T, K における飽和水蒸気圧 PS , MPa は
ln PS =
a
+ b + cT + dT 2 + eT 3 + fT 4 + g ln T
T
ただし，
水と接する場合
T = 27316
. − 47315
. K
氷と接する場合
T = 17315
. − 27316
. K
× 10 3
a = −58002206
.
b = 13914993
.
a = −5.6745359 × 10 3
b = 6.3925247
c = −4.8640239 × 10 −2
c = −9.6778430 × 10 −3
d = 4.1764768 × 10 −5
d = 6.2215701 × 10 −7
× 10 −8
e = −14452093
.
e = 2.0747825 × 10 −9
f =0
g = 6.5459673
f = −9.4840240 × 10 −13
g = 4.1635019
2.1.5
沸点
推定式 1
（文献：理科年表）
TBP = 100.00 + 0.0367( p − 760) − 0.000023( p − 760) 2
TBP ；沸点，C
p ；圧力，mmHg (適用範囲 680～790 mmHg)
推定式 2
Antoine の式
log P = A -
B
T+C
より
7
7/66
埼玉工業大学（小西克享）
T[ K] =
計算熱流体講義ノート（初版）
第2章
物性値の推定式
B
−C
A - log P[ Pa ]
誤差が最小となるように係数を決定すると
A = 10.3527
B = 1885.08
標準沸点，臨界点(374.15℃)間で
・・・広い温度範囲には不向き
C = −20.603
最大誤差　1768%
.
範囲を限定すると
3.92MPa，14.7MPa 間で
A = 10.6435
B = 224117
.
・・・良好
C = 310569
.
最大誤差　0.09077%
2.1.6
蒸発潜熱
沸点 t℃における水の蒸発潜熱(J/g)は
⎧⎪
⎛ t ⎞
⎟
⎝ 100 ⎠
1.15
λ = ⎨287.20 + 3.4173⎜
⎪⎩
tc；
⎛ t ⎞
− 5.7569 × 10 −3 ⎜
⎟
⎝ 100 ⎠
臨界点温度(374.15℃)
適用範囲；0C～臨界点
文献
物性定数推算法
第8版
2.2 空気
2.2.1
密度
ρ=
0.0012932 H
1 + 0.00367t 760
ρ；密度, g/cm3
H：水銀柱 mmHg
t；温度（C）
文献：
改訂 3 版
化学便覧
基礎編
CRC Handbook of Chemistry and Physics
理科年表
2.2.2
粘度
(1) 文献：西山哲男「流体力学(Ⅰ)」日刊工業新聞社
(
μ = 1.75 × 10 −6 1 + 0.0028t − 0.000015t 2
(2) 文献：佐藤一雄「物性定数推算法
)
kgws/m2
第８版」
8
6.5
⎛ t ⎞
− 8.3736 × 10 −17 ⎜
⎟
⎝ 100 ⎠
30 ⎫
⎪
0.365
⎬(t c − t )
⎪⎭
8/66
埼玉工業大学（小西克享）
⎛ T ⎞
μ = 170.9⎜
⎟
⎝ 273⎠
計算熱流体講義ノート（初版）
第2章
物性値の推定式
9/66
0.768
*
10 −6 poise
T ；K (範囲 0～100C)
2.2.3
比熱
C p = A + BT + CT 2 ,
A = 6.27, B = 2.09 × 10 −3 , C = −0.459 × 10 −6
C p ；cal/gram-mole/K
T ；K (273～2073K)
文献：佐藤一雄「物性定数推算法
2.3
第８版」
燃焼ガスの熱力学的性質
2.3.1
燃焼ガス成分のモル定圧比熱
燃焼ガス成分 i のモル定圧比熱［J/(mol K)］は
T ：ガス温度 K
T
θ=
1000
とおくと
θ ≥ 12
. ( T ≥ 1200K) において cp ,i = a1 +
a2
θ
+
a3
θ
2
+
a4
θ
3
+
a5
θ4
θ < 12
. (T < 1200K) において cp ,i = a6 + a7θ + a8θ 2 + a9θ 3 + a10θ 4
係数 a1 ~ a10 は下表のとおり．
2.3.2
燃焼ガス成分のモルエンタルピ
モルエンタルピ［J/mol］を標準反応熱 Δ f H ( T0 ) 基準で
0
H 0 ( T0 ) = ∫ cp ( T )dT + Δ f H 0 ( T0 )
T
0
T0
と定義した場合，
⎛
⎝
θ ≥ 12
. (T ≥ 1200K) において Hi = a11 + 1000⎜ a1θ + a2 ln θ −
⎛
⎝
θ < 12
. (T < 1200K) において Hi = a13 + 1000⎜ a6θ +
a3
θ
−
a4
a ⎞
− 53 ⎟
2
2θ
θ ⎠
a7θ 2 a8θ 3 a9θ 4 a10θ 5 ⎞
+
+
+
⎟
2
3
4
5 ⎠
係数 a1 ~ a15 は下表のとおり．
a1
a2
CO 2
H 2O
O2
N2
CO
H2
O
OH
H
NO
70.909
64.939
45.141
40.721
40.478
53.784
27.214
42.667
20.786
40.217
-54.152
-28.323
-3.882
-20.754 -17.729 -99.861 -46.960 -25.865 0
9
-4.559
埼玉工業大学（小西克享）
a3
a4
a5
a6
a7
a8
a9
a10
a11
a12
a13
a14
a15
132.189
2.999
計算熱流体講義ノート（初版）
-72.113 47.006
物性値の推定式
29.760
-171.147 -8.880
137.533 -63.980 -50.554 -255.638 -147.515 -43.464 0
-8.487
78.519
11.286
-74.290 30.495
23.588
113.023 63.015
22.061
0
17.438
18.522
33.662
29.984
31.796
25.758
31.340
20.786
34.292
82.929
-7.756
-14.183 -17.082 -20.360 19.486
0
-3.052
-80.827
30.428
54.295
4.681
0
77.848
43.865
-19.928
-50.853 -24.725 -35.597 36.750
-17.990 3.401
0
-2.721
-10.181
4.807
15.631
5.555
9.065
0
17.620
-358382
-278517
-50363
5870.1
-108970 98116
310446 38863
211791 88610
243.60
205.60
230.99
217.06
224.70
169.45
139.77
-402094
-251745
-8690.6 -8914.0 -119629 -8250.3 242318 29944
211792 80968
217.25
230.60
243.57
139.77
53051.9
44398.2
38459.8 37026.6 37363.7 35049.5 25921.1 36239.0 24943.2 38447.6
37.444
233.28
221.536 126.439 36.002
10/66
0
31.509
37.968
第2章
48.281
24.773
-15.513 -6.615
-40.469 24.716
-11.315 4.928
119.90
240.26
157.45
194.62
-2.141
202.87
223.43
240.27
259.08
単位： a11 , a13 , a15 は J/mol，それ以外は J/(mol K)
文献
水谷幸夫著「燃焼工学 第２版」森北出版，ｐｐ．６３－６５
2.4 石油
石油類は温度（℃）が T1 から T2 に変化するとき，
T2の比重 = T1の比重 − 比重補正係数×(T2 − T1 )
の関係がある．
測定温度にお １℃当たりの 測定温度にお １℃当たりの 測定温度にお
１℃当たりの
ける測定比重 比重補正係数 ける測定比重 比重補正係数 ける測定比重
比重補正係数
3
3
3
3
kg/m
kg/(m ･℃)
kg/m
kg/(m ･℃)
kg/m
kg/(m3･℃)
640～643
0.95
742～747
0.80
871～890
0.65
644～648
0.94
748～753
0.79
891～970
0.64
649～654
0.93
754～759
0.78
971～1000
0.63
655～661
0.92
760～765
0.77
662～667
0.91
766～771
0.76
668～674
0.90
772～777
0.75
675～681
0.89
778～783
0.74
682～688
0.88
784～790
0.73
689～696
0.87
791～799
0.72
697～703
0.86
800～808
0.71
704～711
0.85
809～818
0.70
712～719
0.84
819～828
0.69
720～726
0.83
829～838
0.68
727～734
0.82
839～852
0.67
735～741
0.81
853～870
0.66
10
3
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
文献：西山善忠「燃料油・潤滑油」海文堂
2.5
第2章
物性値の推定式
11/66
pp.183～185
補間によるデータの算出
推定式が見つからない場合でも，物性値表があれば補間法によって値を推定することができる．
2.5.1
補間法
補間法にはいくつかの種類があるが，一例としてラグランジ補間を解説する．
n 組のデータ ( xi , yi )
(i = 1,2,3,L, n ) が与えられるとき， x における y の補間値は次式で求めら
れる．
n
y=∑
k =1
Pk ( x )
yk
Pk ( x k )
ここで，
n
Pk ( x ) =
特徴；
∏( x − xi )
i =1
x − xk
基点が等間隔でなくても良い
基点数が多くなると計算時間が著しく増加
2.5.2
水の動粘性係数
補間法で水の動粘性係数を求める．
フォームに貼り付けて実行(Visual Basic Ver.2)
Dim TotalData%, Xdata(100), Ydata(100)
‘メインルーチン
Sub Form_Click()
call main(41) ‘温度をセットし，main をコール
End
Sub main(x0)
Call DataSet 'データ設定
Call lagrange(x0, f0) '補間サブルーチンコール
Print f0
End Sub
Sub DataSet ()
'水の動粘度設定
TotalData% = 14
'C
m2/s
11
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第2章
物性値の推定式
12/66
Xdata(1) = 0:Ydata(1) = .000001792
Xdata(2) = 5:Ydata(2) = .00000152
Xdata(3) = 10:Ydata(3) = .000001307
Xdata(4) = 15:Ydata(4) = .000001139
Xdata(5) = 20:Ydata(5) = .0000010038
Xdata(6) = 25:Ydata(6) = .000000893
Xdata(7) = 30:Ydata(7) = .000000801
Xdata(8) = 40:Ydata(8) = .000000658
Xdata(9) = 50:Ydata(9) = .000000658
Xdata(10) = 60:Ydata(10) = .000000475
Xdata(11) = 70:Ydata(11) = .000000413
Xdata(12) = 80:Ydata(12) = .000000365
Xdata(13) = 90:Ydata(13) = .000000326
Xdata(14) = 100:Ydata(14) = .000000295
End Sub
Sub lagrange (x0, f0)
'lagrangian interpolation ラグランジ補間
'TotalData%= ﾃﾞｰﾀ数 (maximum 100)
'Xdata(i%) = x(100)
'Ydata(i%) = f(100)
Static p(100)
Dim i%, j%
For j% = 1 To TotalData%
p(j%) = 1
For i% = 1 To TotalData%
If i% <> j% Then
p(j%) = p(j%) * (x0 - Xdata(i%)) / (Xdata(j%) - Xdata(i%))
End If
Next
Next
f0 = 0
For i% = 1 To TotalData%
f0 = f0 + p(i%) * Ydata(i%)
Next
End Sub
例題
P0 = 760mmHg, t 0 = 30° C, 相対湿度φ 0 = 70%, 体積1m 3 の湿り空気を 10°C まで冷却す
ると空気の相対湿度は何%となるか．また，水蒸気凝縮量（ドレン量）はいくらか
12
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第2章
物性値の推定式
熱量 Q
湿り空気
P0 , t 0 , φ 0
P1 , t 1 , φ 1
t1
冷却器
水滴
飽和蒸気圧は LANGE’S HANDBOOK OF CHEMISTRY より
log10 P[mmHg] = A -
B
C +t
変形すると
P = 10
A-
B
C +t
係数は，t=0℃～60℃の範囲において
A = 8.10765,B = 1750.286,C = 235.0
飽和水蒸気圧は
t=30℃において 31.827mmHg
t=10℃において 9.1966mmHg
相対湿度は
φ=
( = 水蒸気分圧 )
PH 2 O
PS
( = 飽和水蒸気圧 )
× 100 %
であるから，除湿前の湿り空気の水蒸気分圧は
PH 2 O = φPS =
70
× 31827
= 22.2789 mmHg > 9.1966mmHg
.
100
となって，除湿後の温度における飽和水蒸気圧より高くなる．
すなわち，除湿後の空気中の水蒸気は飽和状態となり，相対湿度は 100%となる．
除湿される水蒸気は
22.2789-9.1966=13.0823 mmHg
分である．
湿り空気の体積が 1m3 とすると，除湿される蒸気の体積は
1×
13.0823
= 0.01721 m 3
760
標準状態では
0.01721 ×
27315
.
= 0.01551 m 3
27315
. + 30
質量にすると
13
13/66
埼玉工業大学（小西克享）
0.015510 ( m 3 ) ×
計算熱流体講義ノート（初版）
第2章
物性値の推定式
18.02 (g)
= 12.5 g
22.4 × 10 −3 ( m 3 )
レポート；富士山頂では何Ｃで水が沸騰するか
ヒント；蒸気圧に関するアントワンの式を利用する．
解答
理科年表によれば富士山頂の年平均気圧は 638.3mb (1941-1970 平均)
638.3 mb = 638.3 ×
760
= 478.74 mmHg
1013.3
LANGE’S HANDBOOK OF CHEMISTRY より
B
C+t
t = 0 − 60C
A = 810765
.
, B = 1750.286, C = 235.0
t = 60 − 150C A = 7.96681, B = 1668.21, C = 228.0
log 10 P[ mmHg] = A -
変形し，t=60-150C と仮定して，数値を代入すると
t=
1668.21
B
−C =
− 228.0 = 87.5 C
7.96681 − log 10 478.74
A − log 10 P[ mmHg]
仮定した範囲内なので，解とする．
14
14/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
第3章
特性曲線法による油圧管路動特性解析
3.1 仮定
管壁の弾性を無視
管内の温度変化はなく，流体の密度は一定
１次元流れ
層流
3.2
運動方程式と連続の式
ナビエストークスの方程式において微小項を無視すると運動方程式は
∂p
∂u
+ρ
+ ρRu = 0
∂x
∂t
(1)
連続の方程式は
∂u 1 ∂ p
+
=0
∂ x K ∂t
(2)
ここで
p：
圧力
x：
管路に沿う座標
ρ： 油の密度
u：
流速
t：
時間
R = 8ν
r2
ν： 油の動粘度
r：
管内半径
K：
油の体積弾性係数
圧力波の伝播速度を c とおくと
c=
3.3
K ρ よりK = ρ c 2
(3)
特性曲線法による変換
運動方程式と連続の式を変形し L1, L 2 とおくと
1 ∂ p ∂u
+
+ Ru = 0
ρ ∂ x ∂t
∂u ∂ p
L2 = ρ c 2
+
=0
∂ x ∂t
L1 =
(4)
(5)
未定定数 λ を用いて線形結合する．
15
15/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
L = L1 + λL2 =
=
16/66
⎛
∂u ∂ p⎞
1 ∂ p ∂u
⎟
+
+ Ru + λ ⎜⎜ ρ c 2
+
ρ ∂ x ∂t
∂ x ∂ t ⎟⎠
⎝
∂ p ∂u
∂u
1∂p
+λ
+
+ λρ c 2
+ Ru
ρ ∂x
∂t ∂t
∂x
(6)
⎛ 1 ∂ p ∂ p⎞ ⎛
∂u ∂u⎞
⎟⎟ + ⎜⎜ λρ c 2
⎟ + Ru = 0
= λ ⎜⎜
+
+
∂ x ∂ t ⎟⎠
⎝ λρ ∂ x ∂ t ⎠ ⎝
圧力および流速の解を p = p( x , t ), u = u( x , t ) とすると，微分して
dp ∂ p dx ∂ p
=
+
dt ∂ x dt ∂ t
(7)
du ∂ u dx ∂ u
=
+
dt ∂ x dt ∂ t
つぎに
dx
1
=
= λρ c 2
dt λρ
(8)
であるとき
λ
⎛ ∂ p dx ∂ p ⎞ ∂ u dx ∂ u
dp du
⎟⎟ +
+
+ Ru = λ ⎜⎜
+
+
+ Ru
dt dt
⎝ ∂ x dt ∂ t ⎠ ∂ x dt ∂ t
⎛ 1 ∂ p ∂ p⎞ ⎛
∂u ∂u⎞
⎟⎟ + ⎜⎜ λρ c 2
⎟ + Ru = 0
= λ ⎜⎜
+
+
∂ x ∂ t ⎟⎠
⎝ λρ ∂ x ∂ t ⎠ ⎝
式(8)より
λ=±
1
ρc
(10)
また
dx
= ±c
dt
(11)
となる．＋は圧力の進行波，－は圧力の後退波となる．
式(10)を式(9)に代入すると
±
1 dp du
+
+ Ru = 0
ρ c dt dt
よって
dp ± ρ cdu ± ρRucdt = 0
(12)
となる．＋は圧力の進行波，－は圧力の後退波である．
進行波では式(11)より
dx = cdt
式(12)より
dp + ρ cdu + ρRucdt = dp + ρ cdu + ρRudx = 0
同様に後退波では
dx = −cdt
16
(9)
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
dp − ρ cdu − ρRucdt = dp − ρ cdu + ρRudx = 0
以上，まとめると
3.4
⎧dx − cdt = 0
進行波 ⎨
⎩dp + ρ cdu + ρRudx = 0
(13)
⎧dx + cdt = 0
後退波 ⎨
⎩dp − ρ cdu + ρRudx = 0
(14)
数値計算法
管路を Δx = cΔt で分割する． x − t 平面上での特性線は下図のとおりである．
P
t
Δt
進行波
A
後退波
Δx
Δx
B
x
i-1
i
i+1
式(13)より
dp + ρ cdu + ρ cRudt = 0
差分近似すると
pi − pi −1 + ρ c(ui − ui −1 ) + ρ cui −1 Δt = 0
よって，進行波では
pi = pi −1 − ρ cui + ρ cui −1 (1 − RΔt )
(15)
同様に，後退波では
pi = pi +1 − ρ cui + ρ cui +1 (1 − RΔt )
(16)
式(15)(16)より圧力と流速について解けば
ρc
1
( pi −1 + pi +1 ) + (1 − RΔt )(ui −1 − ui +1 )
2
2
1
1
ui =
( pi −1 − pi +1 ) + (1 − RΔt )(ui −1 + ui +1 )
2ρ c
2
pi =
次に，
C1,i = ui − C2 pi − ui RΔt
C2 = 1 ρ c
C3,i = ui + C2 pi − ui RΔt
17
(17)
(18)
17/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
とおくと，式(17)(18)は
C3,i −1 − C1,i +1
pi =
ui =
3.5
(19)
2 C2
C1,i +1 + C3,i −1
(20)
2
油圧要素
3.5.1
基礎式
境界における計算式は
管左端；
後退波の式(16)を適用
管右端；
進行波の式(15)を適用
方程式１つに未知数２つ．もう一つ方程式が必要．
残りの方程式は管路に接続された油圧要素の特性式を適用する．
管路左端 i = 1 において
p1 =
u1 =
C3,0 − C1, 2
(21)
2 C2
C1, 2 + C3, 0
(22)
2
i = 0 の要素は存在しないから，消去する必要がある．式(21)より
C3, 0 = C1, 2 + 2C2 p1
式(22)に代入すると
u1 =
C1, 2 + C3, 0
=
C1, 2 + C1, 2 + 2C2 p1
2
同様に，管路右端 i = n では
pn =
un =
2
= C1, 2 + C2 p1
(23)
C3,n −1 − C1,n +1
2 C2
C1,n +1 + C3,n −1
2
i = n + 1 の要素は存在しないから，消去する必要がある．
C1,n +1 = C3,n −1 − 2C2 pn
したがって
un =
3.5.2
C1,n +1 + C3,n −1
2
=
C3,n −1 − 2C2 pn + C3,n −1
2
閉塞端
左端が閉じている場合
18
= C3,n −1 − C2 pn
(24)
18/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
u1 = 0
したがって，式(23)より
p1 = −
C1, 2
C2
右端が閉じている場合
un = 0
したがって
pn =
C3,n −1
C2
左端閉
i=1
i=n
右端閉
3.5.3
タンク
左端がタンクに接続されている場合
p1 = const
流速は式(23)が適応でき
u1 = C1, 2 + C2 p1
右端がタンクに接続されている場合
pn = const
流速は式(24)が適応でき
un = C3,n −1 − C2 pn
p1 = const
左タンク
i=n
i=1
pn = const
右タンク
3.5.4
ポンプ
左端がポンプに接続されている場合
u1 = const
19
19/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
20/66
圧力は式(23)より
p1 =
u1 − C1, 2
C2
右端がポンプに接続されている場合
un = const
圧力は式(24)より
pn =
− un + C3,n −1
C2
u1 = const
左ポンプ
i=1
i=n
un = const
右ポンプ
3.5.5
オリフィス
オリフィスを挟んで２つの管路が接続されている場合，オリフィスの左側の管路を添字 L，右
側の管路を R で表すこととする．
uL, n = uR ,1
管路Ｌ
管路Ｒ
流量連続の関係から
uL, n = uR ,1
オリフィスの流量特性は
Q = u L, n A = Cd Aori
2
ρ
(p
L, n
− pR ,1
)
ここで
A
Cd
管路断面積
オリフィス流量係数
20
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
21/66
オリフィス開口面積
Aori
圧力は式(23),(24)より
pR ,1 =
pL, n =
3.6
( )
u L , n − C1, 2
R
C2
− u L , n + ( C 3 , n −1 ) L
C2
Water hammering （水撃現象）の計算
長さ L=10 m，半径 r=10mm の管内を水が一様に流れている．時刻 t=0 から１秒間で管路右端を
完全に閉めきるとき，管内の圧力変化を求めよ．
p0
ρg
u0
x
1
3.6.1
2
i
基礎式
基礎式は式(17)(18)より
pi =
ui =
C3,i −1 − C1,i +1
2 C2
C1,i +1 + C3,i −1
2
ここで
C1,i = ui − C2 pi − ui RΔt
C2 = 1 ρ c
C3,i = ui + C2 pi − ui RΔt
R = 8ν
r2
管路長 L を m 分割すると，
Δx =
L
m
計算刻みは
21
n
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
Δt =
3.6.2
Δx
c
初期条件
初期値は管路内のすべての場所において
p = p0 = 1 × 10 6 Pa , u = u0 = 100mm / s
3.6.3
①
境界条件
タンク接続部
タンク接続部には油圧要素 5.3 を適用
p1 = p0
u1 = C1,1 + C2 p1
②
バルブ側
a. t < dt ：バルブ閉動作時
右端は図のように流速が減少するものとする．
u0
0
dt
t⎞
⎛
un = u0 ⎜ 1 − ⎟
⎝ dt ⎠
pn =
− un + C3,n −1
C2
b. t ≥ dt ：バルブ閉後
un = 0
C
pn = 3,n −1
C2
3.6.4
物性値
c:水中の音速=1.500x106 mm/s (@23-27C) 理科年表
ν:水の動粘性係数=0.893 mm2/s (@25C) 理科年表
ρ:水の密度 0.99704x10-3 g/mm3 理科年表
22
22/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
3.6.5
計算準備
m=10
L 1000
=
= 100 mm
m
10
Δx
100
2
Δt =
=
= × 10 −5 s
6
c
3
. × 10
15
8ν 8 × 0.893
R= 2 =
1/ s
r
10 3
Δx =
3.6.6
プログラム例“Water hammering”
使用言語
Visual Basic
Dim C1i(200), C3i(200)
Dim p(200), u(200)
Dim c, nu, rho, L, j%, m%, n%, deltaX, deltaT, R, C2, p0, u0, dt
Private Sub Form_Click()
Call main
End Sub
Private Sub main()
Call DataSet
For j% = 1 To 2000 '時間刻み 2000 まで
t = deltaT * j%
'時刻
For i% = 1 To n%
C1i(i%) = C1(i%)
C3i(i%) = C3(i%)
Next
'境界（管路端）における計算
'タンク接続部
p(1) = p0
u(1) = C1i(2) + C2 * p(1)
'バルブ端
If t < dt Then
u(n%) = u0 * (1 - t / dt)
p(n%) = (-u(n%) + C3i(n% - 1)) / C2
Else
u(n%) = 0
p(n%) = C3i(n% - 1) / C2
End If
'境界内部の計算
23
23/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
For i% = 2 To n% - 1
p(i%) = (C3i(i% - 1) - C1i(i% + 1)) / 2 / C2
u(i%) = (C1i(i% + 1) + C3i(i% - 1)) / 2
Next
label1 = j%
'
Call plot1
Call plot2
Next
Print "end"
End Sub
Private Sub DataSet()
'mm/s 水中の音速
c = 1.5 * 1000000#
'mm2/s 動粘性係数
nu = 0.893
rho = 0.99704 * 0.001 'g/mm3 密度
'mm 管路長
L = 10 * 1000
'mm 管半径
Rp = 10
m% = 50
'管路分割数
n% = m% + 1
'境界端番号
'mm 管路刻み
deltaX = L / m%
deltaT = deltaX / c
's 時間刻み
R = 8 * nu / (Rp * Rp)
C2 = 1 / rho / c
p0 = 0.1 * 1000000#
'Pa 初期圧力
u0 = 100#
dt = 0.04
'mm/s 初期流速
's バルブ閉時間
For i% = 1 To n%
p(i%) = p0
'管内圧力初期値設定
u(i%) = u0
'管内流速初期値設定
Next
End Sub
Private Function C1(i%)
C1 = (1 - R * deltaT) * u(i%) - C2 * p(i%)
End Function
Private Function C3(i%)
C3 = (1 + R * deltaT) * u(i%) + C2 * p(i%)
End Function
Public Sub plot1()
'バルブ側の圧力履歴表示
24
24/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第3章
特性曲線法による油圧管
路動特性解析
If j% Mod 5 = 0 Then
Scale (1, -0.4 * 1000000#)-(2000, 0.2 * 1000000#)
Line (1, 0)-(2000, 0), QBColor(7)
PSet (j%, -p(n%)), QBColor(7)
DoEvents
End If
End Sub
Public Sub plot2()
'管内の圧力分布と速度分布表示
If j% Mod 5 = 0 Then
Cls
'圧力分布
Scale (1, -0.2 * 1000000#)-(51, 0.2 * 1000000#)
Line (1, 0)-(51, 0), QBColor(7)
For i% = 1 To n% - 1
Line (i%, -p(i%))-(i% + 1, -p(i% + 1)), QBColor(7)
Next
'速度分布
Scale (1, -200#)-(51, 200#)
For i% = 1 To n% - 1
Line (i%, -u(i%))-(i% + 1, -u(i% + 1)), QBColor(1)
Next
DoEvents
End If
End Sub
25
25/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
第4章
4.1
26/66
常微分方程式の数値解法＝ルンゲクッタ法
1 階常微分方程式
１階常微分方程式
dy
= f ( x, y)
dx
を初期条件 x = x 0 , y = y 0 のもとで解くとき， x = x 0 + h における解は
y = y0 + k
で与えられる．ただし，
k=
1
( k1 + 2 k 2 + 2 k 3 + k 4 )
6
ここで
k 1 = hf ( x 0 , y 0 )
k ⎞
h
⎛
k 2 = hf ⎜ x 0 + , y 0 + 1 ⎟
⎝
2
2⎠
k ⎞
h
⎛
k 3 = hf ⎜ x 0 + , y 0 + 2 ⎟
⎝
2
2⎠
k 4 = hf ( x 0 + h, y 0 + k 3 )
4.2
2 階常微分方程式
一般に，2 階常微分方程式は 2 つの常微分方程式に変換し，連立して解くことで求められる．
例えば，2 階常微分方程式
d2y
=y
dx 2
は１階連立常微分方程式
dz
= y,
dx
dy
=z
dx
に変換できる．
例，重力加速度に関する式
d2x
=g
dt 2
は
dv
= g,
dt
dy
=v
dt
と変換できる．
１階連立常微分方程式
dy
dz
= f ( x, y, z ),
= g ( x, y , z )
dx
dx
を初期値 x = x 0 , y = y 0 , z = z 0 のもとで解くとき， x = x 0 + h における解は
26
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
27/66
y = y0 + k
z = z0 + l
で与えられる．ただし，
1
( k1 + 2 k 2 + 2 k 3 + k 4 )
6
1
l = (l1 + 2l2 + 2l3 + l4 )
6
k=
ここで
⎪⎧k 1 = hf ( x 0 , y 0 , z 0 )
⎨
⎪⎩l1 = hg( x 0 , y 0 , z 0 )
k1
l1 ⎞
⎧
h
⎛
⎪k 2 = hf ⎜⎝ x 0 + 2 , y 0 + 2 , z 0 + 2 ⎟⎠
⎪
⎨
⎪l = hg ⎛⎜ x + h , y + k1 , z + l1 ⎞⎟
⎪⎩ 2
⎝ 0 2 0 2 0 2⎠
k2
l2 ⎞
⎧
h
⎛
⎪k 3 = hf ⎜⎝ x 0 + 2 , y 0 + 2 , z 0 + 2 ⎟⎠
⎪
⎨
⎪l = hg ⎛⎜ x + h , y + k 2 , z + l2 ⎞⎟
⎪⎩ 3
⎝ 0 2 0 2 0 2⎠
⎧⎪k 4 = hf ( x 0 + h, y 0 + k 3 , z 0 + l3 )
⎨
⎪⎩l4 = hg( x 0 + h, y 0 + k 3 , z 0 + l3 )
4.3
プランジャポンプとノズルの計算
l
d
y plunger ,
u plunger =
dy plunger
Vplunger, Pp lunger
ymax
dt
Aplunger
dt
プランジャーの動き
4.3.1
基礎式
基礎式は式(17)(18)より
27
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
pi =
C3,i −1 − C1,i +1
2 C2
ここで
C1,i = ui − C2 pi − ui RΔt
C2 = 1 ρ c
C3,i = ui + C2 pi − ui RΔt
R = 8ν
r2
管路長 L を m 分割すると，
Δx =
L
m
計算刻みは
Δt =
4.3.2
Δx
c
初期条件
初期値は管路内のすべての場所において
pi = p0 = 1 × 10 6 Pa
ui =Cd
4.3.3
(1)
2
ρ
( p0 − pair )
Anozzle
A pipe
境界条件
プランジャポンプ側
フックの法則より，流体の圧力変化は体積歪みに比例し，
dp = K
− dv
v
と表される.ここで，K は体積弾性率．
質量一定なら，体積が減少すると(dv<0)，圧力が増加する(dp>0)．
プランジャポンプ室では
dv = − A plunger u plunger + A pipe u1
プランジャ上昇により体積減少
流出による質量減少＝残りは体積増加
28
28/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
したがって,プランジャポンプ室圧力変化は次式を解いてもとめる．
dPplunger
dt
K
=
V plunger
( A plunger u plunger − A pipe u1 )
接続部の圧力,流速は
p1 = p plunger
u1 = C1,1 + C2 p1
(2)
ノズル側
ノズルから噴出する水の速度は
2
unozzle = Cd
ρ
( pn − pair )
接続部の圧力,流速は
un = unozzle
pn =
4.3.4
Anozzle
A pipe
− un + C3,n −1
C2
物性値
c:水中の音速=1.500x106 mm/s (@23-27C) 理科年表
ν:水の動粘性係数=0.893 mm2/s (@25C) 理科年表
ρ:水の密度 0.99704x10-3 g/mm3 理科年表
4.3.5
計算準備
m=50
L 1000
=
= 20 mm
m
50
20
2
Δx
Δt =
=
=
× 10 −5 s
6
15
c
. × 10
15
Δx =
29
29/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
R=
8ν 8 × 0.893
=
1/ s
10 3
r2
4.3.6
プログラム例“Nozzle”
使用言語
Visual Basic
Dim C1i(200), C3i(200)
Dim p(200), u(200)
Dim c, nu, rho, Koil
Dim L, Rpipe, Apipe, m%, n%, deltaX
Dim Cd, Rnozzle, Anozzle
Dim Vmin, Ymax, dt, Rplunger, Aplunger
Dim j%, deltaT, R, C2, p0, u0
Dim Y1, DefEq1
Private Function C1(i%)
C1 = (1 - R * deltaT) * u(i%) - C2 * p(i%)
End Function
Private Function C3(i%)
C3 = (1 + R * deltaT) * u(i%) + C2 * p(i%)
End Function
'データ設定
Private Sub DataSet()
'定数
Const pi = 3.14159
'流体
c = 1.5 * 1000000#
nu = 0.893
'mm/s 流体中の音速
'mm2/s 動粘性係数
rho = 0.99704 * 0.001 'g/mm3 密度
Koil = rho * c * c
'g/mm/s2 体積弾性係数
'大気
Pair = 0.1 * 1000000# 'Pa 大気圧力
'管
L = 1 * 1000
'mm 管路長
Rpipe = 10
'mm 管半径
Apipe = pi * Rpipe * Rpipe 'mm2 管断面積
m% = 50
'管路分割数
n% = m% + 1
'境界端番号
deltaX = L / m%
'mm 管路刻み
'ノズル
30
30/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
'- 速度係数
Cd = 1
'mm ノズル半径
Rnozzle = 0.5
Anozzle = pi * Rnozzle * Rnozzle
'mm2 ノズル断面積
'プランジャポンプ
Vmin = 200
'mm3 ポンプ室最小容積
Ymax = 50
'mm ストローク
's 移動時間
dt = 0.1 '0.04
'mm プランジャ半径
Rplunger = 20
Aplunger = pi * Rplunger * Rplunger 'mm2 プランジャ断面積
deltaT = deltaX / c
's 時間刻み
R = 8 * nu / (Rpipe * Rpipe) '粘性抵抗
C2 = 1 / rho / c
p0 = 0.1 * 1000000#
'Pa 初期圧力
u0 = Cd * Sqr(2 / rho * (p0 - Pair)) * Anozzle / Apipe
'mm/s 初期流速
'u0 = 100#
For i% = 1 To n%
p(i%) = p0
'管内圧力初期値設定
u(i%) = u0
'管内流速初期値設定
Next
End Sub
'微分方程式
Sub DifEq1(t, y)
Dim Uplunger
Vplunger = Vmin + Ymax * (1 - t / dt) * Aplunger
If t < dt Then
Vplunger = Vmin + Ymax * (1 - t / dt) * Aplunger
Uplunger = Ymax / dt
Else
Vplunger = Vmin
Uplunger = 0
End If
DefEq1 = Koil / Vplunger * (Aplunger * Uplunger - Apipe * u(1))
End Sub
'メインルーチン
Private Sub main()
Call DataSet
Y1 = p0
'R.K.M.初期値
For j% = 1 To 2000 '時間刻み 2000 まで
t = deltaT * j%
'時刻
31
31/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
For i% = 1 To n%
C1i(i%) = C1(i%)
C3i(i%) = C3(i%)
Next
'境界（管路端）における計算
'プランジャポンプ
Call RK4M1(t, Y1)
If Y1 < 0 Then Y1 = 0
p(1) = Y1
u(1) = C1i(2) + C2 * p(1)
'ノズル
If p(n%) - Pair < 0 Then
u(n%) = 0
Else
u(n%) = Cd * Sqr(2 / rho * (p(n%) - Pair)) * Anozzle / Apipe
End If
p(n%) = (-u(n%) + C3i(n% - 1)) / C2
If p(n%) < 0 Then p(n%) = 0
'境界内部の計算
For i% = 2 To n% - 1
p(i%) = (C3i(i% - 1) - C1i(i% + 1)) / 2 / C2
u(i%) = (C1i(i% + 1) + C3i(i% - 1)) / 2
Next
label1 = j%
Call plot1
'
Call plot2
Next
Print "end"
End Sub
'圧力履歴表示
Public Sub plot1()
If j% Mod 1 = 0 Then
Scale (1, -100)-(2000, 100)
Line (1, 0)-(2000, 0), QBColor(7)
PSet (j%, -p(1) / 1000000#), QBColor(7)
DoEvents
End If
End Sub
'管内の圧力分布と速度分布表示
Public Sub plot2()
32
32/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
If j% Mod 1 = 0 Then
Cls
'圧力分布
Scale (1, -20)-(51, 20)
Line (1, 0)-(101, 0), QBColor(7)
For i% = 1 To n% - 1
Line (i%, -p(i%) / 1000000#)-(i% + 1, -p(i% + 1) / 1000000#), QBColor(7)
Next
'速度分布
Scale (1, -4000#)-(51, 4000#)
For i% = 1 To n% - 1
Line (i%, -u(i%))-(i% + 1, -u(i% + 1)), QBColor(1)
Next
DoEvents
End If
End Sub
'ルンゲクッタ４次解法
Sub RK4M1(x, y)
Static JR%
Static K1(4) As Single
Static KC(2), HC(3)
X0 = x
Y0 = y
Call DifEq1(X0, Y0)
K1(1) = deltaT * DefEq1
Call DifEq1(X0 + deltaT / 2, Y0 + K1(1) / 2)
K1(2) = deltaT * DefEq1
Call DifEq1(X0 + deltaT / 2, Y0 + K1(2) / 2)
K1(3) = deltaT * DefEq1
Call DifEq1(X0 + deltaT, Y0 + K1(3))
K1(4) = deltaT * DefEq1
Y1 = Y0 + (K1(1) + 2 * (K1(2) + K1(3)) + K1(4)) / 6
End Sub
33
33/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
4.4
34/66
プランジャポンプ＆油圧ジャッキの計算
Acyl
y w , uw =
W
dy w
dt
l
Vcyl, Pcyl
d
y plunger ,
Vplunger, Pplunger
dy plunger
u plunger =
dt
Aplunger
ymax
ymax
dt
dt
プランジャーの動き
4.4.1
おもりの動き
基礎式
基礎式は式(17)(18)より
pi =
ui =
C3,i −1 − C1,i +1
2 C2
C1,i +1 + C3,i −1
2
ここで
C1,i = ui − C2 pi − ui RΔt
C2 = 1 ρ c
C3,i = ui + C2 pi − ui RΔt
R = 8ν
r2
34
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
管路長 L を m 分割すると，
L
m
Δx =
計算刻みは
Δx
c
Δt =
4.4.2
初期条件
初期値は管路内のすべての場所において
pi = p0 = 1 × 10 6 MPa , ui = u0 = 0m / s
4.4.3
(1)
境界条件
プランジャポンプ側
プランジャポンプ室圧力変化は次式を解いてもとめる．
dPplunger
dt
=
K
V plunger
( A plunger u plunger − A pipe u1 )
接続部の圧力,流速は
p1 = p plunger
u1 = C1,1 + C2 p1
(2)
シリンダ側
おもりの運動方程式は
{ (
)
}
duw
1
Acyl Pcyl − Pair − M w g
=
dt
Mw
dy w
= uw
dt
シリンダ内の圧力変化
dPcyl
dt
=
K
( A pipe un − Acyl uw )
Vcyl
接続部の圧力,流速は
pn = pcylinder
u n = C 3 , n −1 − C 2 p n
4.4.4
物性値
c:水中の音速=1.500x103 m/s (@23-27C) 理科年表
ν:水の動粘性係数=0.893x10-6 m2/s (@25C) 理科年表
ρ:水の密度 0.99704x103 kg/m3 理科年表
35
35/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
g:重力加速度=9.8m/s
pair:大気圧=0.1x106Pa
4.4.5
計算準備
m=200
L
1
=
= 0.005 m
m 200
Δx
0.005
1
Δt =
=
= × 10 −5 s
3
c
3
15
. × 10
8ν 8 × 0.893
R= 2 =
1/ s
r
0.010 2
Δx =
4.4.6
プログラム例“Cylinder”
使用言語
Visual Basic
Dim C1i(201), C3i(201)
Dim p(201), u(201)
Dim grav, Pair
Dim c, nu, rho, Koil
Dim L, Rpipe, Apipe, m%, n%, deltaX
Dim Vmin, Ymax, dt, Rplunger, Aplunger, Yplunger
Dim j%, deltaT, R, C2, p0, u0
Dim y1, y2, y3, y4
Dim DefEq1, DefEq2, DefEq3, DefEq4
Dim Mw, Vcylmin, Rcyl, Acyl
Private Sub Form_Click()
Call main
End Sub
Private Function C1(i%)
C1 = (1 - R * deltaT) * u(i%) - C2 * p(i%)
End Function
Private Function C3(i%)
C3 = (1 + R * deltaT) * u(i%) + C2 * p(i%)
End Function
'データ設定
Private Sub DataSet()
'定数
Const pi = 3.14159
36
36/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
grav = 9.81
'm/s2 重力加速度
'流体
'm/s 流体中の音速
c = 1.5 * 1000#
nu = 0.893 / 1000000# 'm2/s 動粘性係数
rho = 0.99704 * 1000# 'kg/m3 密度
'kg/m/s2 体積弾性係数
Koil = rho * c * c
'大気
Pair = 0.1 * 1000000# 'Pa 大気圧力
'管
'm 管路長
L=1
'm 管半径
Rpipe = 10 / 1000
Apipe = pi * Rpipe * Rpipe 'm2 管断面積
m% = 200
'管路分割数
n% = m% + 1
'境界端番号
'm 管路刻み
deltaX = L / m%
'油圧シリンダ
'kg 重り質量
Mw = 1000
Vcylmin = 20 / 1000000
'm3 最小容積
Rcyl = 20 / 1000
'm シリンダ半径
Acyl = pi * Rcyl * Rcyl 'm2 シリンダ断面積
'プランジャポンプ
Vmin = 20 / 1000000
'm3 ポンプ室最小容積
Ymax = 5 / 1000
'm ストローク
dt = 0.02 '0.04
's 移動時間
Rplunger = 20 / 1000
'm プランジャ半径
Aplunger = pi * Rplunger * Rplunger 'm2 プランジャ断面積
deltaT = deltaX / c
's 時間刻み
'dt/deltaT
R = 8 * nu / (Rpipe * Rpipe) '粘性抵抗
C2 = 1 / rho / c
p0 = 0.1 * 1000000#
'Pa 初期圧力
u0 = 0
For i% = 1 To n%
p(i%) = p0
'管内圧力初期値設定
u(i%) = u0
'管内流速初期値設定
Next
37
37/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
y1 = p0
y2 = 0
y3 = 0
y4 = p0
End Sub
'微分方程式
Sub DifEq1(t, y1, y2, y3, y4)
Dim Uplunger
'ポンプ室圧力
If t < dt Then
Yplunger = Ymax * t / dt
Vplunger = Vmin + (Ymax - Yplunger) * Aplunger
'Vplunger = Vmin + Ymax * (1 - t / dt) * Aplunger
Uplunger = Ymax / dt
Else
Yplunger = Ymax
Vplunger = Vmin
Uplunger = 0
End If
DefEq1 = Koil / Vplunger * (Aplunger * Uplunger - Apipe * u(1))
'シリンダ側
'おもり運動方程式
DefEq2 = (Acyl * (y4 - Pair) / Mw - grav)
DefEq3 = y2
'圧力
Vcyl = Vcylmin + y3 * Acyl
DefEq4 = Koil / Vcyl * (Apipe * u(n%) - Acyl * y2)
End Sub
'メインルーチン
Private Sub main()
Call DataSet
For j% = 1 To 20000 '時間刻み 20000 まで
t = deltaT * j%
'時刻
For i% = 1 To n%
C1i(i%) = C1(i%)
C3i(i%) = C3(i%)
Next
'境界（管路端）における計算
Call RK4M1(t, y1, y2, y3, y4)
38
38/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
'プランジャポンプ
If y1 < 0 Then y1 = 0
p(1) = y1
u(1) = C1i(2) + C2 * p(1)
'油圧シリンダ
'If y4 < 0 Then y4 = 0
p(n%) = y4
u(n%) = C3i(n% - 1) - C2 * p(n%)
'境界内部の計算
For i% = 2 To n% - 1
p(i%) = (C3i(i% - 1) - C1i(i% + 1)) / 2 / C2
u(i%) = (C1i(i% + 1) + C3i(i% - 1)) / 2
If p(i%) <= 0 Then p(i%) = 0
Next
label1 = j%
Call plot1
'
Call plot2
Next
Print "end"
End Sub
'圧力履歴表示
Public Sub plot1()
If j% Mod 10 = 0 Then
Scale (1, -100)-(2000, 0)
Line (1, 0)-(2000, 0), QBColor(7)
PSet (j% / 10, -p(n%) / 1000000#), QBColor(7)
Scale (1, -0.006)-(2000, 0.006)
Line (1, 0)-(2000, 0), QBColor(7)
PSet (j% / 10, -y3), QBColor(7)
PSet (j% / 10, -Yplunger), QBColor(1)
DoEvents
End If
End Sub
'管内の圧力分布と速度分布表示
Public Sub plot2()
If j% Mod 10 = 0 Then
Cls
'圧力分布
Scale (1, -20)-(201, 20)
39
39/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
Line (1, 0)-(201, 0), QBColor(7)
For i% = 1 To n% - 1
Line (i%, -p(i%) / 1000000#)-(i% + 1, -p(i% + 1) / 1000000#), QBColor(7)
Next
'速度分布
Scale (1, -4#)-(201, 4#)
For i% = 1 To n% - 1
Line (i%, -u(i%))-(i% + 1, -u(i% + 1)), QBColor(1)
Next
DoEvents
End If
End Sub
'ルンゲクッタ４次解法
Sub RK4M1(x, y1, y2, y3, y4)
Static JR%
Static K1(4) As Single
Static K2(4) As Single
Static K3(4) As Single
Static K4(4) As Single
X0 = x
Y10 = y1
Y20 = y2
y30 = y3
y40 = y4
If y40 < 0 Then y40 = 0
Call DifEq1(X0, Y10, Y20, y30, y40)
K1(1) = deltaT * DefEq1
K2(1) = deltaT * DefEq2
K3(1) = deltaT * DefEq3
K4(1) = deltaT * DefEq4
y10n = Y10 + K1(1) / 2
y20n = Y20 + K2(1) / 2
y30n = y30 + K3(1) / 2
y40n = y40 + K4(1) / 2
Call DifEq1(X0 + deltaT / 2, y10n, y20n, y30n, y40n)
K1(2) = deltaT * DefEq1
K2(2) = deltaT * DefEq2
K3(2) = deltaT * DefEq3
40
40/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
K4(2) = deltaT * DefEq4
y10n = Y10 + K1(2) / 2
y20n = Y20 + K2(2) / 2
y30n = y30 + K3(2) / 2
y40n = y40 + K4(2) / 2
Call DifEq1(X0 + deltaT / 2, y10n, y20n, y30n, y40n)
K1(3) = deltaT * DefEq1
K2(3) = deltaT * DefEq2
K3(3) = deltaT * DefEq3
K4(3) = deltaT * DefEq4
y10n = Y10 + K1(3)
y20n = Y20 + K2(3)
y30n = y30 + K3(3)
y40n = y40 + K4(3)
Call DifEq1(X0 + deltaT, y10n, y20n, y30n, y40n)
K1(4) = deltaT * DefEq1
K2(4) = deltaT * DefEq2
K3(4) = deltaT * DefEq3
K4(4) = deltaT * DefEq4
y1 = Y10 + (K1(1) + 2 * (K1(2) + K1(3)) + K1(4)) / 6
y2 = Y20 + (K2(1) + 2 * (K2(2) + K2(3)) + K2(4)) / 6
y3 = y30 + (K3(1) + 2 * (K3(2) + K3(3)) + K3(4)) / 6
y4 = y40 + (K4(1) + 2 * (K4(2) + K4(3)) + K4(4)) / 6
If y4 < 0 Then
y2 = 0 '負圧になったとき，おもりの運動はいったん停止
y4 = 0
End If
End Sub
41
41/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
4.5
プランジャポンプ＆燃料弁の計算
y needle ,
F
l
uneedle =
d
y plunger ,
Aneedle
Ahole
dt
Pcham
Aplunger
ymax
ymax
dt
dt
プランジャーの動き
4.5.1
ニードルの動き
基礎式
基礎式は式(17)(18)より
ui =
C3,i −1 − C1,i +1
2 C2
C1,i +1 + C3,i −1
2
ここで
C1,i = ui − C2 pi − ui RΔt
C2 = 1 ρ c
C3,i = ui + C2 pi − ui RΔt
R = 8ν
r2
管路長 L を m 分割すると，
42
dy needle
dt
Vcav, Pcav
mneedle
Vplunger, Pplunger
dy plunger
u plunger =
pi =
42/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
L
m
Δx =
計算刻みは
Δx
c
Δt =
4.5.2
初期条件
初期値は管路内のすべての場所において
pi = p0 = 1 × 10 6 MPa , ui = u0 = 0m / s
4.5.3
(1)
境界条件
プランジャポンプ側
プランジャポンプ室圧力変化は次式を解いてもとめる．
dPplunger
dt
=
K
V plunger
( A plunger u plunger − A pipe u1 )
接続部の圧力,流速は
p1 = p plunger
u1 = C1,1 + C2 p1
(2)
燃料弁側
ニードルの運動方程式は
duneedle
1
=
dt
mneedle
dy needle
= uneedle
dt
{( A
needle
− Ahole ) Pcav + Ahole Pcham − F − mneedle g}
ここで，ばねがニードルを押す力は
F = F0 + k s y needle
噴孔から流出する流体の速度
u jet =Cd
2
ρ
( pcav − pcham )
キャビティー内の圧力変化
dPcav
K
=
( A pipe un − Aneedle uneedle − Ahole u jet )
dt
Vcav
接続部の圧力,流速は
pn = pcav
u n = C 3 , n −1 − C 2 p n
43
43/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
4.5.4
物性値
c:水中の音速=1.500x103 m/s (@23-27C) 理科年表
ν:水の動粘性係数=0.893x10-6 m2/s (@25C) 理科年表
ρ:水の密度 0.99704x103 kg/m3 理科年表
g:重力加速度=9.8m/s
pcham:容器圧力=0.1x106Pa
4.5.5
計算準備
m=400
L
1
=
= 0.0025 m
m 400
0.0025 5
Δx
Δt =
=
= × 10 − 6 s
c
15
. × 10 3 3
8ν 8 × 0.893
R= 2 =
1/ s
r
0.012
Δx =
4.5.6
プログラム例“Fuel Valve”
使用言語
Visual Basic
Dim C1i(1001), C3i(1001),p(1001), u(1001)
Dim grav, c, nu, rho, Koil
Dim L, Rpipe, Apipe, m%, n%, deltaX
Dim Vmin, Ymax, dt, Rplunger, Aplunger, Yplunger
Dim F0, Ks, Rhole, Ahole, Mneedle, Rneedle, Aneedle, Vcav, Ylift, Cd
Dim Pcham
Dim j%, deltaT, R, C2, p0, u0
Dim y1, y2, y3, y4, DefEq1, DefEq2, DefEq3, DefEq4
Dim Vcylmin, Rcyl, Acyl
Private Function C1(i%)
C1 = (1 - R * deltaT) * u(i%) - C2 * p(i%)
End Function
Private Function C3(i%)
C3 = (1 + R * deltaT) * u(i%) + C2 * p(i%)
End Function
'データ設定
Private Sub DataSet()
'定数
Const pi = 3.14159
grav = 9.81
'm/s2 重力加速度
44
44/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
'流体
'm/s 流体中の音速
c = 1.5 * 1000#
nu = 0.893 / 1000000# 'm2/s 動粘性係数
rho = 0.99704 * 1000# 'kg/m3 密度
'kg/m/s2 体積弾性係数
Koil = rho * c * c
'管
'm 管路長
L=1
'm 管半径
Rpipe = 10 / 1000
Apipe = pi * Rpipe * Rpipe 'm2 管断面積
m% = 400
'管路分割数
n% = m% + 1
'境界端番号
'm 管路刻み
deltaX = L / m%
'油圧シリンダ
Vcylmin = 20 / 1000000
'm3 最小容積
Rcyl = 20 / 1000
'm シリンダ半径
Acyl = pi * Rcyl * Rcyl 'm2 シリンダ断面積
'プランジャポンプ
Vmin = 20 / 1000000
'm3 ポンプ室最小容積
Ymax = 5 / 1000
'm ストローク
dt = 0.01 '0.04
's 移動時間
Rplunger = 20 / 1000
'm プランジャ半径
Aplunger = pi * Rplunger * Rplunger 'm2 断面積
'燃料弁
F0 = 200
'N 初期締付力
Ks = 300
'N/m ばね定数
Rhole = 1 / 1000 'm 噴孔半径
Ahole = pi * Rhole * Rhole 'm2 噴孔断面積
Mneedle = 100 / 1000
'kg ニードル質量
Rneedle = 5 / 1000
'm ニードル半径
Aneedle = pi * Rneedle * Rneedle 'm2 ニードル断面積
'm3 キャビティ容積
Vcav = 2 / 1000000#
'm ニードルリフト量
Ylift = 1 / 1000
Cd = 1
'流量係数
'容器
Pcham = 0.1 * 1000000#
deltaT = deltaX / c
'Pa 容器圧力
's 時間刻み
R = 8 * nu / (Rpipe * Rpipe) '粘性抵抗
C2 = 1 / rho / c
p0 = 0.2 * 1000000#
'Pa 初期圧力
u0 = 0
45
45/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
For i% = 1 To n%
p(i%) = p0
'管内圧力初期値設定
u(i%) = u0
'管内流速初期値設定
Next
y1 = p0
y2 = 0
y3 = 0
y4 = p0
End Sub
'微分方程式
Sub DifEq1(t, y1, y2, y3, y4)
Dim Uplunger
'ポンプ室圧力
If t < dt Then
Yplunger = Ymax * t / dt
Vplunger = Vmin + (Ymax - Yplunger) * Aplunger
'Vplunger = Vmin + Ymax * (1 - t / dt) * Aplunger
Uplunger = Ymax / dt
Else
Yplunger = Ymax
Vplunger = Vmin
Uplunger = 0
End If
DefEq1 = Koil / Vplunger * (Aplunger * Uplunger - Apipe * u(1))
'燃料弁側
'ニードル運動方程式
F = F0 + Ks * y3
DefEq2 = ((Aneedle - Ahole) * y4 + Ahole * Pcham - F - Mneedle * grav) / Mneedle
DefEq3 = y2
'圧力
If y3 > 0 Then
Ujet = Cd * Sqr(2 / rho * (y4 - Pcham))
Else
Ujet = 0
End If
DefEq4 = Koil / Vcav * (Apipe * u(n%) - Aneedle * y2 - Ahole * Ujet)
End Sub
'メインルーチン
Private Sub main()
46
46/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
Call DataSet
For j% = 1 To 20000 '時間刻み 20000 まで
t = deltaT * j%
'時刻
For i% = 1 To n%
C1i(i%) = C1(i%)
C3i(i%) = C3(i%)
Next
'境界（管路端）における計算
Call RK4M1(t, y1, y2, y3, y4)
'プランジャポンプ
If y1 < 0 Then y1 = 0
p(1) = y1
u(1) = C1i(2) + C2 * p(1)
'油圧シリンダ
p(n%) = y4
u(n%) = C3i(n% - 1) - C2 * p(n%)
'境界内部の計算
For i% = 2 To n% - 1
p(i%) = (C3i(i% - 1) - C1i(i% + 1)) / 2 / C2
u(i%) = (C1i(i% + 1) + C3i(i% - 1)) / 2
If p(i%) <= 0 Then p(i%) = 0
Next
label1 = j%
Call plot1
Next
Print "end"
End Sub
'表示
Public Sub plot1()
If j% Mod 100 = 0 Then
'燃料弁入り口圧力
Scale (1, -50)-(4000, 5)
DrawWidth = 2
PSet (j% / 10, -p(n%) / 1000000#), QBColor(12)
'リフト
Scale (1, -0.006)-(4000, 0.006)
PSet (j% / 10, -y3), QBColor(9)
PSet (j% / 10, -Yplunger), QBColor(14)
DoEvents
End If
47
47/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
End Sub
'ルンゲクッタ４次解法
Sub RK4M1(x, y1, y2, y3, y4)
Static JR%
Static K1(4) As Single
Static K2(4) As Single
Static K3(4) As Single
Static K4(4) As Single
X0 = x
Y10 = y1
Y20 = y2
y30 = y3
y40 = y4
Call DifEq1(X0, Y10, Y20, y30, y40)
K1(1) = deltaT * DefEq1
K2(1) = deltaT * DefEq2
K3(1) = deltaT * DefEq3
K4(1) = deltaT * DefEq4
y10n = Y10 + K1(1) / 2
y20n = Y20 + K2(1) / 2
y30n = y30 + K3(1) / 2
y40n = y40 + K4(1) / 2
Call DifEq1(X0 + deltaT / 2, y10n, y20n, y30n, y40n)
K1(2) = deltaT * DefEq1
K2(2) = deltaT * DefEq2
K3(2) = deltaT * DefEq3
K4(2) = deltaT * DefEq4
y10n = Y10 + K1(2) / 2
y20n = Y20 + K2(2) / 2
y30n = y30 + K3(2) / 2
y40n = y40 + K4(2) / 2
Call DifEq1(X0 + deltaT / 2, y10n, y20n, y30n, y40n)
K1(3) = deltaT * DefEq1
K2(3) = deltaT * DefEq2
K3(3) = deltaT * DefEq3
K4(3) = deltaT * DefEq4
48
48/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第4章
常微分方程式の数値解法
＝ルンゲクッタ法
y10n = Y10 + K1(3)
y20n = Y20 + K2(3)
y30n = y30 + K3(3)
y40n = y40 + K4(3)
Call DifEq1(X0 + deltaT, y10n, y20n, y30n, y40n)
K1(4) = deltaT * DefEq1
K2(4) = deltaT * DefEq2
K3(4) = deltaT * DefEq3
K4(4) = deltaT * DefEq4
y1 = Y10 + (K1(1) + 2 * (K1(2) + K1(3)) + K1(4)) / 6
y2 = Y20 + (K2(1) + 2 * (K2(2) + K2(3)) + K2(4)) / 6
y3 = y30 + (K3(1) + 2 * (K3(2) + K3(3)) + K3(4)) / 6
y4 = y40 + (K4(1) + 2 * (K4(2) + K4(3)) + K4(4)) / 6
If y3 < 0 Then
y2 = 0
y3 = 0
ElseIf y3 > Ylift Then
y2 = 0
y3 = Ylift
End If
If y4 < 0 Then
y4 = 0
End If
End Sub
49
49/66
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
第5章
5.1
１次元圧縮性流れ
50/66
１次元圧縮性流れ
質点の運動
図に示す質点の運動を考える．質点とは，質量を持つものの大きさが無視できる物体のことで
ある．
x
m
F=
d2x
dt 2
質点の運動方程式は
m
d2x
= F = Fg (1)
dt 2
もしくは
du
= F = Fg
dt
dx
=u
dt
m
(2)
と表わすことができる．ここで
m：質点の質量
u：質点の速度
F：質点に作用する力．例，重力 Fg=mg
数値解法（ルンゲ・クッタ法）を用いて(2)式を解けば，質点の位置と速度が求まる．
5.2
剛体の運動
周囲との摩擦が存在する場合
m
d2x
= F = Fg − FR
dt 2
ここで
Fg:重力
FR：摩擦抵抗
5.3
ガス体の運動
5.3.1
流速に関する式
剛体の運動方程式をガス体の運動に適用することにより解を求める．ただし，剛体の代わりに
流体内部に微小な要素（コントロールボリューム）を考え，この要素の運動を考えることになる．
50
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
F
１次元圧縮性流れ
51/66
F’
mp
P
P’
FR
⇓ 等価
F
mp
F’
FR
コントロールボリュームに関して，運動方程式は
mp
d2x
= FP − FR
dt 2
mp
du
= FP − FR
dt
(3)
もしくは
dx
=u
dt
(4)
ただし，
mp：要素（コントロールボリューム）内の質量＝質点に置き換える
u：流速＝質点の速度
FR：粘性抵抗＝管摩擦抵抗
ここで
F P= F − F '
(5)
はコントロールボリューム前後の圧力差により発生する力（＝質点に作用する力に相当）を表わ
す．
計算格子は固定されるため，質点の位置 x を解くことは無意味．
mp
du
= FP − FR
dt
(6)
のみを解けばよい．
5.3.2
圧力に関する式
圧力はエネルギ式から求める．
(1)
エネルギ保存の式
コントロールボリュームには，左側から流体が速度 u& で流入し，右側から速度 u& ′ で流出するも
のとし，外部から熱量 Q が流入する．
51
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
Q
u&
u& ′ (H’)
(H)
ρ, C P , T
ρ′ , C ′P , T ′
コントロールボリュームに関するエネルギ関係式は
dH + dQ = dU + dL
(7)
ただし，
dH：要素のエンタルピ変化＝流出入する流体のエンタルピ差
dH = H − H '
dQ：要素へのエネルギ流入（出）量
＝要素の加熱 or 冷却
dU：要素の内部エネルギ変化
dU =
1
(VdP + κPdV )
κ −1
dL：要素のする膨張仕事
dL = PdV
コントロールボリュームの体積を一定とすると dV=0
1
VdP
(VdP + κPdV ) + PdV =
κ −1
κ −1
κ −1
∴ dP =
(dH + dQ)
V
dH + dQ =
(8)
単位時間当たりでは
dP κ − 1 ⎛ dH dQ ⎞
=
+
⎜
⎟
dt
V ⎝ dt
dt ⎠
(9)
従って，1 次元圧縮性流れに関して解くべき基礎式は
質点の運動に関して
m
du
= FP − FR
dt
(10)
圧力に関して
dP κ − 1 ⎛ dH dQ ⎞
=
+
⎜
⎟
dt
V ⎝ dt
dt ⎠
(2)
(11)
計算格子
52
１次元圧縮性流れ
52/66
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
１次元圧縮性流れ
53/66
圧力と流速の格子点をずらす．（スタッガード格子）
格子点は圧力に関し１から n まで．速度に関し，１から n-1 までとする．
圧力に関するｺﾝﾄﾛｰﾙﾎﾞﾘｭｰﾑ
Pi-1
Pi
Pi+1
ui-1
ui
P i+2
ui+1
速度に関するｺﾝﾄﾛｰﾙﾎﾞﾘｭｰﾑ
圧力格子点間（圧力に関する i 番目の格子内）のガス質量を質量 mp,i の剛体に置き換える．剛体の
速度は流速 ui に相当する．ただし，下線付き添え字は圧力の格子からずれた（スタッガード）格
子点であることを示す．
速度の決定に関する格子
Pi-1
mp,i-1
Pi
ui-1
Pi+1
mp,i
ui
圧力差により i 番目の要素に作用する力は
F P,i = Ai Pi − Ai +1 Pi +1
(12)
ここで
Ai：管断面積
運動方程式は
m p ,i
du i
dt
= FP,i − FR,i
(13)
= Ai Pi − Ai +1 Pi +1 − FR,i
ただし，
m p ,i = ρiVi ,R + ρi +1Vi +1,L
(14)
Vi,R：圧力の格子点 i より右半分の体積
Vi,Ｌ：圧力の格子点 i より左半分の体積
密度は
53
mp,i+1 ui+1
P i+2
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
ρ=
P
RT
１次元圧縮性流れ
(15)
より，(13)式は
Pi
P
Vi ,R + i +1 Vi +1,L
RTi
RTi +1
m p ,i =
(16)
⎞
P
1⎛P
= ⎜ i Vi ,R + i +1 Vi +1,L ⎟
R ⎝ Ti
Ti +1
⎠
したがって，(13)式に代入すると，
⎞ du i
P
1 ⎛ Pi
= Ai Pi − Ai +1 Pi +1 − FR ,i
⎜ Vi ,R + i +1 Vi +1,L ⎟
R ⎝ Ti
Ti +1
⎠ dt
整理すると，
du i
dt
=R
Ai Pi − Ai +1 Pi +1 − FR ,i
Pi
P
Vi ,R + i +1 Vi +1,L
Ti
Ti +1
(17)
圧力の決定に関する格子
Pi
Ti
mi
Pi-1
ui-1
Pi+1
Ti+1
mi+1
P i+2
ui
ui+1
エンタルピの増加は
dHi
= A i −1 ρ i −1u i −1C P ,i −1T i −1 − A i ρ i u i C P ,i T i
dt
P i −1
Pi
= A i −1
u i −1C P ,i −1T i −1 − A i
u i C P ,i T i
RT i −1
RT i
=
(
1
A i −1 P i −1u i −1C P ,i −1 − A i P i u i C P ,i
R
(18)
)
(11)式より
dPi κ − 1 ⎛ dHi dQi ⎞
=
+
⎟
⎜
dt
Vi ⎝ dt
dt ⎠
=
κ −1⎧ 1
(
⎨ A i −1 P i −1u i −1C p , i −1 − A i P i u i C p , i
Vi ⎩ R
圧力は平均値として
54
)
dQ ⎫
+ i⎬
dt ⎭
(19)
54/66
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
Pi −1 + Pi
2
P + Pi +1
Pi = i
2
１次元圧縮性流れ
55/66
P i −1 =
(20)
を用いることにすると，
dPi κ − 1 ⎧ 1 ⎛
Pi −1 + Pi
P + Pi +1
⎞ dQ ⎫
=
u i −1C p , i −1 − A i i
u i C p, i ⎟ + i ⎬
⎨ ⎜ A i −1
⎠ dt ⎭
2
2
dt
Vi ⎩ R ⎝
(21)
各要素の質量変化は
dmi
= ρ i −1u i −1 − ρ i u i
dt
(22)
各要素の温度は
Ti =
PV
i i
mi R
(24)
まとめ
1 次元圧縮性流れに関して解くべき式は
質点の運動に関して
du i
dt
=R
Ai Pi − Ai +1 Pi +1 − FR ,i
Pi
P
Vi ,R + i +1 Vi +1,L
Ti
Ti +1
(23)
圧力に関して
Pi −1 + Pi
P + Pi +1
dPi κ − 1 ⎧ 1 ⎛
⎞ dQ ⎫
u i −1C p , i −1 − A i i
u i C p, i ⎟ + i ⎬
=
⎨ ⎜ A i −1
⎠ dt ⎭
dt
Vi ⎩ R ⎝
2
2
質量に関して
dmi
= ρ i −1u i −1 − ρ i u i
dt
(25)
温度に関して
Ti =
(3)
PV
i i
mi R
(26)
境界条件
Ａ．管路両端の境界値を圧力値とする場合
管路上の格子点は下図の通り
55
(24)
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
i=1
i=2
u1
i=1
P n-1
Pi
u2
ui
i=i
i=2
56/66
i=n
i=i
P2
P1
１次元圧縮性流れ
Pn
un-1
i=n-1
境界条件 P1 と Pn を計算する．
流速は(23)式より u1 から un-1 まで計算する．
圧力は(24)式より P2 から Pn-1 まで計算する．
①
管路端開放
Pn = const = P0
(27)
とおく．
②
管路端閉
エンタルピの増加は
P n −1
dH n
1
u n −1C P ,n −1T n −1 = A n −1 P n −1u n −1C P ,n −1
= A n −1 ρ n −1u n −1C P ,n −1T n −1 = A n −1
dt
RT n −1
R
(28)
(11)式より
dPn κ − 1 ⎛ dH n dQn ⎞
=
+
⎜
⎟
dt
Vn ⎝ dt
dt ⎠
=
κ −1⎧ 1
Vn
dQn ⎫
⎨ A n−1 P n−1u n−1C p , n−1 +
⎬
dt ⎭
⎩R
(29)
圧力は平均値として
P n−1 =
Pn−1 + Pn
2
(30)
を用いることにすると，
Pn−1 + Pn
dQ ⎫
dPn κ − 1 ⎧ 1
u n−1C p , n−1 + n ⎬
=
⎨ A n−1
dt ⎭
dt
Vn ⎩ R
2
を解けばよい．
Ｂ．管路の１端の境界値を流速とする場合
管路上の格子点は下図の通り
56
(31)
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
i=1
i=2
i =i
P1
P2
Pi
Pn-1
Pn
u1
u2
ui
un-1
un
i =1
i=2
i=i
i=n-1
i=n
流速は(23)式より u1 から un-1 まで計算する．
圧力は(24)式より P2 から Pn まで計算する．
流速一定
un = const = u 0
②
57/66
i=n
境界条件 P1 と un を計算する．
①
１次元圧縮性流れ
差圧から計算
un =
2
ρ
( pn − pback )
57
埼玉工業大学（小西克享）
第6章
計算熱流体講義ノート（初版）
第6章
熱の流れ
熱の流れ
熱伝導方程式（偏微分方程式）を解いて温度分布を求める．
6.1
偏微分方程式
代表的な偏微分方程式は次の 3 種類
6.2
∂ 2 f ( x, y) ∂ 2 f ( x, y)
+
= ρ( x , y )
∂x 2
∂y 2
Poisson 方程式
ρ = 0 のとき
Laplace 方程式
∂ f ( x, y)
∂ 2 f ( x, y)
=a
∂t
∂x 2
熱伝導方程式
（放物型偏微分方程式）
2
∂ 2 f ( x, y)
2 ∂ f ( x, y)
=
a
∂t 2
∂x 2
波動方程式
（双曲型偏微分方程式）
（楕円型偏微分方程式）
差分近似
6.2.1
導関数の近似法
テイラー展開を利用し，導関数を近似する．代表的なものは次の 3 種類．
u(x)
u(x+h)
u(x)
u(x-h)
x-h
x
x+h
x
中心差分近似
1 階導関数
du
1
≈
{u( x + h) − u( x − h)}
dx 2h
誤差 −
1 2
h f ′′′
6
誤差 −
1 2
h f
12
2 階導関数
d 2u
1
≈ 2 {u( x + h) − 2u( x ) + u( x − h)}
2
dx
h
前進差分近似
1 階導関数
du 1
≈ {u ( x + h) − u ( x)}
dx h
誤差 −
1
hf ′′
2
2 階導関数
d 2u
1
≈ 2 {u( x + 2h) − 2u( x + h) + u( x )}
2
dx
h
58
誤差 − hf ′′′
Ⅳ
58/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第6章
熱の流れ
後退差分近似
1 階導関数
du 1
≈ {u ( x) − u ( x − h)}
dx h
誤差 −
1
hf ′′
2
2 階導関数
d 2u
1
≈ 2 {u( x ) − 2u( x − h) + u( x − 2h)}
2
dx
h
6.2.2
誤差 − hf ′′′
陽解法
１ステップ前の計算時刻における値のみを用いて，次の時刻における値を計算する．
t
ui,j+1
j+1
j
ui-1,j ui,j
k
ui+1,j
h
i-1
i
i+1
x
方程式
∂ u ∂2u
=
∂t ∂ x2
は
ui , j +1 − ui , j
k
=
ui +1, j − 2ui , j + ui −1, j
h2
空間微分に対する中心差分近似
時間微分に対する前進差分近似
と近似することができる．
6.2.3
Crank-Nicolson の陰解法
１ステップ前の計算時刻における値のみを用いて，次の時刻における値を計算する．
方程式
∂ u ∂2u
=
∂t ∂ x2
は
ui , j +1 − ui , j
k
=
1 ⎧ ui +1, j +1 − 2ui , j +1 + ui −1, j +1 ui +1, j − 2ui , j + ui −1, j ⎫
+
⎨
⎬
2⎩
h2
h2
⎭
59
59/66
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第6章
熱の流れ
60/66
と近似することができる．
t
ui-1,j+1
ui,j+1 ui+1,j+1
j+1
j
ui-1,j ui,j
k
ui+1,j
h
i-1
i
i+1
x
例題（陽解法）
次の Fourier の偏微分方程式で表される１次元温度場の平行平面板内の温度変化を求めよ．ただし，
刻み幅は h=1, k=1 とする．
∂φ
∂ 2φ
1
=a 2 a=
∂t
2
∂x
境界条件 φ ( 0, t ) = φ ( 4, t ) = 20
初期条件 φ ( x , t ) = 5x( x − 4) + 20
解答
差分化式は
φi , j +1 − φi , j
k
=
1 φi +1, j − 2φi , j + φi −1, j
2
h2
h=1, k=1 なので
φi , j +1 − φi , j =
φi , j +1 =
(
φi , j +1
(
(
1
φi +1, j − 2φi , j + φi −1, j
2
)
)
1
φi +1, j − 2φi , j + φi −1, j + φi , j
2
1
= φi +1, j + φi −1, j
2
)
境界条件と初期条件から
φ0,0 = φ4 ,0 = 20, φ1,0 = 5, φ2 ,0 = 0, φ3,0 = 5
t=1 では
1
1
φ2 ,0 + φ0,0 ) = (0 + 20) = 10
(
2
2
1
1
= (φ3,0 + φ1,0 ) = (5 + 5) = 5
2
2
1
1
= (φ4 ,0 + φ2 ,0 ) = (20 + 0) = 10
2
2
φ1,1 =
φ2 ,1
φ3,1
t=２では
60
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第6章
熱の流れ
61/66
1
1
φ2 ,1 + φ0,1 ) = (5 + 20) = 12.5
(
2
2
1
1
= (φ3,1 + φ1,1 ) = (10 + 10) = 10
2
2
1
1
= (φ4 ,1 + φ2 ,1 ) = (20 + 5) = 12.5
2
2
φ1,2 =
φ2 ,2
φ 3, 2
t
7
6
5
4
3
2
1
0
x
6.2.4
20
20
20
20
20
20
20
20
0
18.75
18.125
17.5
16.25
15
12.5
10
5
1
18.125
17.5
16.25
15
12.5
10
5
0
2
18.75
18.125
17.5
16.25
15
12.5
10
5
3
20
20
20
20
20
20
20
20
4
偏微分方程式の無次元化
熱伝導方程式
∂T
∂2T
=a 2
∂t
∂x
(1)
において各記号の意味と次元は次の通りである．
T: 温度(K)
t: 時間(h)
(2)
x: 距離(m)
a: 温度伝導率(m2/h)
従って，(1)式は(K/h)の次元を持つ式となっている．
(1)式を無次元化するには，まず，各変数を無次元化する．無次元化のためには各変数と同じ次
元を持つ基準を設定して，各々除した変数を作ればよい．無次元変数を，￣の記号をつけて表記
する，(2)式を無次元化した変数は次のようになる．
T − T0
(T: T0 → T1 , T: 0 → 1)
T1 − T0
t
t=
(t:0 → t 0 , t: 0 → 1)
t0
x
x=
( x:0 → x0 , x: 0 → 1)
x0
a
a=
a0
T=
(3)
(4)
(5)
(6)
ここで，
61
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
第6章
熱の流れ
62/66
T1: 例，板の初期温度(K)
T0: 例，周囲温度(K)
t0: 例，t=t0 となるまでの時間(h)
(7)
x0: 例，板厚(m)
2
a0 =
x0
: 温度伝導率(m2/h)
t0
次に，微分定理を用いて(4-1)式の偏微分を無次元変数の微分に変換する．
∂
∂ ∂t
∂ 1
1 ∂
=
= ⋅ =
∂t ∂t ∂t ∂t t 0 t 0 ∂t
(8)
∂
∂ ∂x
∂ 1
1 ∂
=
=
⋅
=
∂x ∂ x ∂x ∂ x x0 x0 ∂ x
(9)
∂2
∂ ⎛ ∂⎞
∂ ⎛ ∂ ⎞ ∂x
∂ ⎛ 1 ∂⎞ 1
1 ∂ ⎛ ∂⎞
1 ∂
=
= 2
=
⎜
⎟⋅
⎜ ⎟= 2
⎜ ⎟⋅
⎜ ⎟=
2
∂x ⎝ ∂x ⎠ ∂ x ⎝ ∂x ⎠ ∂x ∂ x ⎝ x 0 ∂ x ⎠ x 0 x 0 ∂ x ⎝ ∂ x ⎠ x 0 ∂ x
∂x
(10)
(3)から(10)式を用いて(1)式を変換すると
1 ∂
1 ∂2
{(T1 − T0 )T + T0 } = a x 2 2
t 0 ∂t
0 ∂x
{(T
1
− T0 )T + T0 }
∂T
a ∂2 T
=
∂t ⎛ x 0 2 ⎞ ∂ x
⎜
⎟
⎝ t0 ⎠
∂T
∂2 T
=a
∂t
∂x
(11)
となり，(1)式を無次元化した方程式が得られる．
なお，(11)式の差分化において時間と空間の刻み幅は，実時間，実空間での刻みを k と h で表すと，
それぞれ
k
t0
h
h=
x0
k=
(12)
(13)
となる．
無次元化方程式を解くことの意味
(1)式を解くことはある板厚，温度におけるただ 1 つの解を得ることを意味しているのに対して，
(11)式を解くことは相似となるすべての条件について解くことを意味している．
たとえば板の材質が同じ場合，次の例 1 と例 2 は無次元化した(11)式ではまったく同じとなる．
例 1．計算時間 t0=100sec，時間刻み k=1sec，空間刻み h=1mm
62
埼玉工業大学（小西克享）
計算熱流体講義ノート（初版）
T1 = 110 C
（一様）
T0 = 10 C
x 0 = 10 mm
例 2．計算時間 t0=400sec，時間刻み k=4sec，空間刻み h=2mm
T1 = 100 C
（一様）
T0 = 0 C
x 0 = 20 mm
.
例 3．計算時間 t 0 = 1 ，時間刻み k = 0.01 ，空間刻み h = 01
T1 = 1
（一様）
T0 = 0
x0 = 1
63
第6章
熱の流れ
63/66
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
第7章
１次元圧縮性流れ
数値計算上の注意事項
打切り誤差に注意が必要
例 1. 01
.4
+4
01
.2
+44
L+3
01
. ≠1
1
10
コンピュータで計算すると
01
.4
+4
01
.2
+44
L+3
01
. − 1 → −1110
.
× 10 −16
1
10
小数の計算には注意が必要．
となる．
20 = 1
2 −1 = 0.5
2 − 2 = 0.25
.
2 − 3 = 0125
2 − 4 = 0.0625
2 −5 = 0.03125
2 − 6 = 0.015625
2 − 7 = 0.0078125
2 −8 = 0.00390625
２進数で 0.1 を表そうとすると
01
. → 2 −4 + 2 −5 + 2 −8 +L
２進数では，すべての小数を正確に表すことができない．
for i = 0 to 1 step 0.1
:
next
とすると，このループは無限ループとなる可能性がある．
例 2. 01
. * 10 ≠ 1
コンピュータで計算すると
01
. * 10 − 1 → 5551
.
× 10 −17
また
01
.4
+4
01
.2
+44
L+3
01
. ≠ 01
. * 10
1
10
であることに注意が必要．
例 3． 10./1. = 10. but
1./0.1 ≠ 10.
コンピュータで計算すると
64
打切り誤差が必ず発生
64/66
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
10./1.−10. → 0.
× 10 −16
.
1./0.1 − 10. → −5551
. ≠ 01
. × 01
.
例 4． 01
コンピュータで計算すると
2
01
. 2 − 01
. × 01
. → 8.326 × 10 −19
累乗の計算には級数展開式を使用．
例 5．整数型の演算
A, B が整数型，C が実数型の変数のとき
A=1
B=3
C=A/B
とすると C の値は？
答
0.
A, B, C とも実数型の変数のときは 0.333333.....
65
１次元圧縮性流れ
65/66
埼玉工業大学（小西克享） 計算熱流体講義ノート（初版） 第 5 章
１次元圧縮性流れ
計算熱流体講義ノート(初版)
計算熱流体講義ノート(初版)
平成 21 年 12 月 1 日 第 2 版
著者承認検印省略
無断転載を禁ず
個人的な学習の目的以外での使用，転載，配布等はできません．
著者
印刷・製本
非売品
小西克享
Copyright ⓒ 2009 小西克享, All Rights Reserved.
66
66/66