研究の背景

Problems of water/oil boiling
トランスミッションギヤ等の油焼入れ工程を想定した沸騰流れ現象の非定常シミュレーション法の研究開発(ex. quenching) Large latent heat release Large density variation
 Large heat transfer
Model of Heat transfer for water boiling (1atm)

先駆的科学計算に関するフォーラム 2019年4月26日 九州大学
Nucleate boiling

北海道大学大学院工学院 機械宇宙工学部門
大島伸行・WANG Lu

Rate of heat transfer

円柱JIS試験

Film boiling

Temperature on solid wall (℃)

2

研究の背景
焼入れ: 高温に加熱された金属を急冷するプロセス

3

研究の背景
interFoam(OpenFOAM ライブラリーより
改良VOF法 に基づく二相流解析ソルバー)

熱処理後の変形ばらつきと低減対策の必要性
 性能の維持 修正加工の削減

追加

油を用いた焼入れでは以下の3段階で冷却速度が異なる• エネルギー方程式• 相変化項蒸気膜段階・沸騰段階・対流段階

interThermalPhaseChangeFoam [2]この現象が不均一に起こるため、ばらつきが生じる• 相変化に伴う熱輸送を考慮した非定常二相流解析
ソルバー

この不均一性は以下の３つの項目が相互に影響することによって発生するので、これらを同時に解析できる手法が求められている(1) 温度分布(2) 蒸気発生量分布(3) 流速分布と気相/液相分布[1] D. S. MacKenzie, “Heat Treatment of Gears: Control of ResidualStress and Distortion”, Heat Treating Progress, Volume 7, Issue, July
2007

ギア団体焼入れ このソルバーの焼入れ問題への適用が期待される
 本解析では簡略化した壁面温度変化モデルと連成して焼入れ問題における冷却速度の定性的な再現を行った
一般的な焼入れ時の金属温度の時間変化[1][2] Mahdi Nabil, Alexander S. Rattner, ”interThermalPhaseChangeFoam—
A framework for two-phase flow simulations with thermally driven phasechange”, SoftwareX, 2016, Volume 5, 216-226基礎方程式 ～VOF法/ 非圧縮性流れ
連続条件

運動量保存式

界面計算手法 ～VOF法とその改良r∇ ⋅ u = v& pcrrrrr∂ ( ρu )+ ∇ ⋅ ( ρu u ) + ∇ p ′ − ∇ ⋅ τ = f V + f S∂t

速度-圧力カップリング（非圧縮性流れ）

Vapor

α =0r
∂ ( ρh )+ ∇ ⋅ ( ρu h) + ∇ ⋅ (−k eff ∇T ) = −q& pc∂tr∂α
VOF法（気液界面）
位相輸送方程式+ ∇ ⋅ (u α ) = α& pc
カップリング∂t生成項： q& pc =∆t6

解析対象

T = 98, α = 1・・・

LiquidLiquidα =1α = 0.5

改良VOF法における方程式（＝保存型A-C方程式） 00rr
∂α+ ∇ ⋅ (u α ) − ∇ ⋅ (Γnum ∇α ) + ∇ ⋅ (u *α ) = α& pc∂t0.2 0.4 0.3 0.2

風上差分法による数値粘性q& pc  1q&1  ， α& pc = − pc−，v& pc =hLV  ρV ρ L ρhLV

気液相変化モデルの基礎検証

Vaporr∂α+ ∇ ⋅ (u α ) = α& pc
従来の位相輸送方程式∂t

熱輸送方程式

ρC p (T − Tsat )・・・ρ＋改良速度場

00.3 0.601

α = 0.50.7 1
Extended Level-set equation
1 1= −α& pc

Dφ+ ρ V ∇ φ − ρV ( f d + ∇ φ )Dt+ f ' ( f d + ∇ φ ) − ∇ ⋅ d ( f d + ∇φ ) n = 00 改良速度場
0 000.8 0.6 0.4

1

1

10.9

1

11

157

気液相変化モデルの基礎検証
気泡の時間変化

adiabatic

圧力理論値： Pσ = 2σ / r
半径rは初期値と相変化量から算出• チュートリアルのケースを一部変更し3次元に拡張• 重力加速度を0とする
T = 100α=0

温度分布と界面の時間変化
T = 98
α=13D geometry and Initial condition
T = 98 [K], α = 1

 気泡の凝縮が観測できる
 気泡半径が小さくなった段階で不安定化する
ものの、全体的に並列数による違いはない

計算領域中心圧力と理論値の時間変化

8

気液相変化モデルの基礎検証
計算開始から0.024秒における気泡中心軸上の分布9

加熱床面モデルでの検証
計算モデルおよび初期条件
気相領域

α =0

T = 100 ℃
液相領域

α =1

T = 100 ℃
半径 0.005m の半球+
加熱壁面から厚さ 0.001m の薄膜
温度分布alpha分布 並列数や格子数による計算結果に大きな違いは見られなかった
 格子解像度が細かいと、界面が薄くなることが確認できる

高温壁
初期温度： 150℃

側面は周期境界条件

10

Fluid property

11

計算結果
壁面温度・熱流束・界面・温度分布の時間変化

Non-physical one

HexaDecane

LiquidVapor

LiquidVaporρ [kg/m3 ]

200.0

5.0

777.534.80615

Cp [J/(kg ∙ K)]

400.0

200.0

2374.72820.6

݇ [W/(m ∙ K)]

40.0

2

ν [m /s]

5.0 × 10

1.0−4

1.0 × 100.14080−33.96 × 10

0.029543−6

1.67 × 10

σ [N/m]0.04290.02747hLV [J/kg]

2.0 × 105

3.5685 × 105

Tsat [K]100560−5◎D. Z. Guo, D. L. Sun, Z. Y. Li, and W. Q. Tao,“Phase Change Heat Transfer Simulation for
Boiling Bubbles Arising from a Vapor Film by the VOSET Method”, Numerical Heat Transfer,Part A, 59: 857–881, 2011 膜沸騰、核沸騰、対流状態の再現
 壁面温度変化を定性的に捉えている
1. 蒸気膜段階
2. 沸騰段階
3. 対流段階

1.non-physical fluid case
１２１２ ３３

壁面温度の時間変化

壁面熱流束の時間変化

12

計算結果
計算開始から4.1秒

Dependent on area size
0.078*0.078*0.078
0.156*0.078*0.156

800

I.C. 150K (unsteady)700

Heat Flux (W)600

SMALL AREA

500
400

LARGE AREA

300
200
100
00

123456789Time (s)0.078*0.078*0.078
0.156*0.078*0.156

150

LARGER AREA SIZE
Twall (K)140 液相に素早く温度境界層が発達→飽和温度を超える領域が広い面積で界面と接し、蒸発量が増える→蒸発熱として壁面からの熱が吸収されていくため、熱流束大130120LARGE AREA

SMALL AREA

110

100

2019/4/19012345678913Time (s)Fixed wall temperature
110K

Dependent on area size

⇒ stable
7000

110K
115K
120K

6000

Heat Flux (W)

5000
4000

120K

3000
2000

110K
120K

1000

⇒ unstable00

2468Time (s)SMALL AREA
0.078*0.078*0.078
2019/4/19

120K

14

2019/4/19

LARGE AREA
0.156*0.078*0.156
15

Dependent on area size

Dependent on resolution
1.non-physical fluid case200*100*200
400*200*400

150

Twall (K)

140

130

127.5K120

110

1000246810Time (s)200*100*200 thickness: 0.003m

3500

200*100*200
400*200*400

3000

Heat Flux (W)

2500
2000

686.96W

1500
1000

SMALL AREA 120K
0.078*0.078*0.078

LARGE AREA 120K
0.156*0.078*0.156

500

654.47W0024

2019/4/19

16

2019/4/19

Finer grid solution
1.non-physical fluid case

6Time (s)400*200*400 thickness: 0.0032m

Real Fluid Case
2.Real fluid property Temperature difference ：50K
200

Heat Flux

150

100

50

00.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.4

0.5

0.6

0.7t/s610.5
610.0609.5

Twall

609.0
608.5
608.0
607.5
607.0
606.5
606.0

2019/4/19

2019/4/19

0.0

0.1

0.2

0.3

t/s

810

Cylinder case with physical fluid property

212.Real fluid property Temperature difference ：50K

110K

Initial condition： Liquid：120℃
Silver：800℃

140K
2019/4/19

まとめ•

界面安定化を導入したレベルセット近似モデル（改良ＶＯＦ法＝拡張レベルセット法）を用いた過熱壁面の冷却（焼き入れ）現象のシミュ
レーションにより過熱面の沸騰遷移過程が模擬された．

•

過加熱（過冷却）された実際の（分子スケール）ミクロ界面において，相変化界面の局所温度は連続，かつ，平衡値に近い極小（極大）をとる⇒ ナノスケールの物理的気液界面での実現性検証•

焼き入れ油の実物性値での解析では、膜沸騰段階での浮力流れ場の速度変動が予測された。核沸騰遷移における乱流予測の必要性がある。
⇒ 乱流効果を考慮した高解像度／界面モデルの導入

23