白金汉大学-写人作文600字

第29卷 第3 期
2 0 0 8 年6月


文章编号：0258-0926(2008)03-0089-05
核 动 力 工 程
Nuclear Power Engineering



Vol. 29. No.3
Jun. 2 0 0 8
球体表面强制对流膜态沸腾换热的数值模拟


袁明豪，杨燕华，李天舒，胡志华

（上海交通大学核科学与工程学院，上海，200240）

摘要：使用改进后的流 体体积法(VOF)追踪汽液自由界面，开发了用于计算带相变的自由界面问题的基
于适体坐标的计算流 体力学(CFD)程序，对球体表面强制对流膜态沸腾换热进行了数值模拟。将数值模拟结
果与实验关联 式进行了对比。结果表明，数值计算较好地模拟了球体表面强制对流膜态沸腾的物理过程。
关键词：自由界面；VOF；膜态沸腾；相变
中图分类号：TK124 文献标识码：A

1 引 言
在轻水堆的堆芯熔化事故或快中子增殖堆
的堆芯解体事故中，高温难挥发的熔融燃料与低
温易挥发的冷却剂之间相互作用。该相互作用可
分为以下几个阶段：高温熔融物团块被释放到冷
却剂中的破裂阶段、粗混合阶段、高温熔融物液
滴的细粒化阶段、压力波传播阶段以及膨胀做功
阶段等。
在粗混合阶段，高温熔融物液滴与 冷却剂形
成一个混合区域，混合区域的状态与范围决定后
续爆炸性膨胀的初始条件与可能转换成 破坏性机
械能的总的热能
[1]
。在该阶段，高温熔融物液滴
的周围被一层蒸 汽膜覆盖，冷却剂液体和熔融物
液滴不能直接接触。高温熔融物液滴形状一般保
持为球形或近似 球形(以下简称高温颗粒)，且在
重力的作用下，在冷却剂液体中向下运动。
Bromley最 早对水平圆柱表面的强迫对流膜态沸
腾进行了理论研究
[2]
，文献[3]使用Bro mley的方
法分析了饱和液体内球体表面的强迫对流膜态沸
腾。文献[4，5]对高温金属球 在水中快速运动时
从膜态沸腾过渡到核态沸腾的过程进行了实验研
究。文献[6，7]通过实验 得到了膜态沸腾条件下
球体的瞬态温度及传热系数，并在此基础上建立
了理论分析模型。文献[ 8]提出了一个新的计算模
型，并用文献[7]的实验数据进行了验证。
随着直接模拟自由界 面的数值方法的发展，
将其用于研究球体表面的强制对流膜态沸腾换
热，有助于揭示一些在实验 中难以观察的现象，
并可对理论模型中的各种假设进行计算验证，进
而更好地理解这一复杂的流 动传热过程。本文使
用流体体积法(VOF)
[9~12]
追踪汽液自由界面，开< br>发了基于适体坐标的计算流体力学(CFD)程序，
对球体表面的强制对流膜态沸腾换热进行数值 模
拟，并与实验关联式进行了对比。

2 物理模型
将液体和蒸汽视为不可压缩牛顿流体，计算
平面中的动量方程为：
∂
(ρ
u
)
+
1
∂
(
ρ
Uu
)< br>+
1
∂
(
ρ
Vu
)
=−
1
∂
tJ
∂
ξ
J
∂
η
J
1
∂
⎡
µ
⎤
−
α
u
β
u
ξη
⎥< br>J
∂
ξ
⎢
⎣
J
⎦
1
∂
⎡< br>µ
⎤

+−
β
u
+
γ
u
ξη
⎥
+2
σκδ
s
n
x
(1)
J
∂
η
⎢
J
⎣⎦
∂
(
ρ
v
)
+
1
∂
(
ρ
Uv
)
+
1
∂
(
ρ
Vv
)
=−
1
∂
tJ
∂
ξ
J
∂
η
J

1
∂
⎡
µ
⎤
−

−
p< br>ξ
x
η
+
p
η
x
ξ
+
ρ< br> g
y
+
α
v
β
v
ξη
⎥
J
∂
ξ
⎢
⎣
J
⎦
1
∂
⎡
µ
⎤

+

−
β
v
+
γ
v
ξη
⎥
+2
σκδ
s
n
y

(2)
⎢
J
∂
η
⎣
J
⎦
⎧
⎪< br>U
=
uy
η
−
vx
η


(3)
⎨
=−
Vvxuy
⎪
ξξ
⎩
p
ξ
y
η
−
p
η
y
ξ
+
ρ
g
x
+

()()
()
()()
()

收稿日期：2007-03-09；修回日期：2008-02-02
基金项目： 国家自然科学基金资助项目(50376036)

90
核 动 力 工 程 Vol. 29. No.3. 2008
⎧
J
=
x
ξ
y
η
−
x
η
y
ξ
⎪
22
⎪
⎪
α
=
x
η
+
y
η

⎨
(4)
β
xxyy
=+
ξηξη
⎪
⎪
γ
=
x
2
+
y
2
⎪
ξξ
⎩
式中，t 为时间，
s
；u、v为直角坐标速度分量，
ms
；
ρ
为密度 ，
kgm
3
；
µ
为动力粘度，
kg(m

· s)
；
式中，V为某个区域；S为包含在V内的自由界
面；
ρ
1
和
ρ
2
为相
1
和相
2< br>的密度，
kgm
3
；m为
通过自由界面的质量流量，
kg(m
2
·

s)
。

在导出式
(5)
的过程中，利用了相
1
和相
2
的
不可压缩性。当存在相变时界面上 的能量条件为：

(
q
1
−q
2
)
•n+ Ψ=mh
lg
(9)
q
和
q
为界面 两侧的热流密度矢量，
p为压力，
Pa
；g
x
、g
y
为重力加速度在x、y方向
的分量，
Nkg
；
σ
为表面张力，Nm
；
κ
为当地的
界面曲率，
m
-
1
；
δ
s
为集中在界面
S
上的狄拉克函
数；n
x、n
y
为界面上的法向矢量在x、y方向的分
量；
ξ
、
η
为计算平面的坐标；U、V为逆变速度
分量，
ms
；J、
α
、
β
、
γ
为坐标变换的几何参数。

计算平面中的能量方程为：

∂
∂
t
(
ρ

C
)
1
p
T
+
∂
J
∂
ξ
(
ρ

C
1
p
UT
)
+
∂
J
∂
η
(
ρ

C
p
VT
)
=

1
∂
⎡
k
J
∂
ξ
⎢
⎣
J
(
α
T
ξ
−
β
T
η
)
⎤
⎥
⎦
+


1
∂
⎡
k
J
∂
η
⎢
⎣
J
(
−
β
T
ξ
+
γ

T
η
)
⎤
⎥
⎦
+Φ
(5)
式中 ，T为温度，
K
；C
p
为定压比热，
J(kg·K)
；k为导热系数，
W(m·K)
；
Φ
为界面上的相变潜热
产生的能 量源项，
Wm
3
。公式
(5)
省略了粘性耗散。

计算平面中的连续性方程为：

1
∂
J
∂
ξ
U
+
1
∂
J
∂
η
V
=
Γ
(6)
式中，Γ为在两种密度不同的流体界面上发生传
质而引起的体 积的膨胀或收缩，
s
－
1
。

式
(1)
、 式
(2)
及式
(5)
中的物理性质在同一相中
保持常数，但不同相之 间的物理性质不同。

为了确定计算区域中某个点属于液相还是
气相，引入了一个特征 函数
χ
，在某一相中取为
1
，
而在另一相中取为
0
，则计算平面中描述
χ
的运动
方程为：

∂
∂
t< br>χ
+
U
∂
J
∂
ξ
χ
+
V< br>∂
J
∂
η
χ
=
0 (7)
VOF
方法中的体积份额函数C可以被看作
特征函数
χ
的离散形式。

当界面上存在蒸发或凝结时，Γ和通过自由
界面的质量流量m之间的关系为：

∫∫∫
Γ
d
v
=
⎜
⎜
V
∫∫
m
⎛
11
⎞
ρ
−
⎟
⎟
d
s

S
⎝
1
ρ
(8)
2
⎠
式中，
12
W(m
2

·s)
；n为界面上的单位法向矢量；h
lg
为相
变潜热，
Jkg
；Ψ为 界面上的辐射热流密度，
W(m
2
·

s)
。

在导出式
(6)
时，进行了一些简化假设，如忽
略了动能和粘性力做功。假定 界面温度等于系统
压力下的平衡饱和温度：

T
s
=T
sat
(
P
∞
)
(10)
式中，
T
s
为界面温度，
K
；
T
sat
为饱和温度，
K
；

P
∞
为系统压力，
Pa
。

当球体的温度较高时， 需要考虑热辐射的作
用。由于蒸汽膜很薄，蒸汽被认为不参于辐射换
热。因此，所有的辐射能量 都在汽
-
液交界面上被
吸收，并为液体蒸发提供了能量来源。球体表面
的辐射 热流密度为：

q
r
=
εσ
b
T
w
4
(11)
式中，
q
r
为球体表面的辐射热流密度，
W(m
2
·

s)
；
ε
为球体表面的发射率；
T
w
为球体表面的温度，
K
；
σ
b
为斯忒藩
-
玻尔兹曼常量，其值为
5.67
×
10
-
8
W(m
2
·

K
4
)
。

根据自由界面的位置和方向，可以将由该式
计算得到的球体表面的辐射热流密度近似分布到
气 液界面上，从而得到式
(9)
中的
Ψ
。


3 数值方法
3.1 算 法
用有限体积法对式
(1)
、式
(2 )
、式
(5)
和式
(6)
进行离散。压力与速度的耦合关系使用SIMPLE
算法求解；其中，式
(6)
用于构造压力修正方程。
在空间 离散时，扩散项采用了二次中心差分，对
流项使用
MUSCL
格式。时间导数的离散使 用全
隐格式的一次向后差分。

用
VOF
方法求解式
(7)
，即使用一个体积份
额函数
C
来表示一个控制体某一相流体
(
如液相
)
所占有的份额。
C
=1
表示全液相控制体；
C
=0
表
示全气相控制体；
0<1<
C

则表示混合物控制体。

袁明豪等：球体表面强制对流膜态沸腾换热的数值模拟
91
为了得到 新时刻
C
的分布，需要知道在一个时间
步长内通过每个控制体表面的某一相流体的体积
通量。
VOF
方法的一个重要特点就是基于控制体
内气液界面的几何信息计算 此体积通量。在已知
的上一时步
C
的分布基础上，通过界面重构，可
以得到一 个控制体内自由界面的方向和位置。

VOF
算法有多种，它们的界面重构方法各不
相同，如
SOLA- VOF
[8]
、
FLAIR
[11]
和
Youngs
VOF
[12]
等。其中，
Youngs VOF
在大多数情况下
性能较好，本文使用
Youngs VOF
进行界面重构。

Youngs VOF
使用由斜率和截距确定的2
维
直线段表示混合物控制体内的自由界面
(PLIC)
，
如图
1
所示。图中的阴影部分表示流体
1
通过混
合物控制体右表面的体积 通量。使用该方法，可
以得到新时刻的
C
的分布。新时刻
C
的分布用 来
确定式
(1)
、式
(2)
和式
(5)
中的物理性 质分布：

⎧
⎪
ρ
=
ρ
l
C
+< br>ρ
g
(
1
−
C
)
⎪
⎨
µ< br>=
µ
l
C
+
µ
g
(
1
−< br>C
)

⎪
k
=
k
l
C
+
k
g
(
1
−
C
)
(12)
⎪
⎩
C
p
=
C
p
l
C
+
C
p
g
(
1
−
C
)
式中，下标< br>l
表示液相，
g
表示气相。



图1 混合物控制体内界面重构及流体1的
体积通量计算
Fig. 1 Interface Reconstruction in a Mixture Cell
and Calculation of Volume Flux for Fluid 1

计算 区域和网格单元在物理平面上是不规
则的，但在计算平面上却是规则区域，网格单元
是规则的矩 形。所以，在计算平面上使用
VOFPLIC
方法求解
(7)
式的过程与直角 坐标下的
求解步骤十分相似。在膜态沸腾的过程中，汽
-
液界面上存在相变，对这一问 题的解决方法详见
文献
[13]
。

3.2 计算步骤
(1)
在上一时刻的流体体积份额
C
的分布和
速度分布的基础上，使用
VOFPLIC
方法取得新
时刻
C
的分布。

(2)在新时刻
C
分布的基础上，使用连续表面
力
(CSF)
模型[14]
计算表面张力。

(3)
由新时刻
C
的分布， 根据式
(12)
确定物理
性质的分布。

(4)
求解能量方程式
(5)
，得到新时刻的温度
场。
< br>(5)
由新时刻
C
的分布计算自由界面的位置
和方向，得到式
(9)
中的辐射热流密度
Ψ
。结合新
时刻的温度场，由式
(9)计算通过界面的相变质量
流量
m
。

(6)
由式
(8)
计算体积的膨胀和收缩，即式
(6)
中的
Γ
。
< br>(7)
用
SIMPLE
算法迭代求解式
(1)
、式
( 2)
和
式
(6)
，得到新时刻的速度场。


4 数值模拟结果的验证
在关于球体表面强制对流膜态沸腾换热的
众多研究中，文献
[7 ]
做了大量不同大小、不同材
料的球体在不同来流速度的水中膜态沸腾的实
验，提出了 计算饱和沸腾换热系数的实验关联式：
⎧
14
⎪
⎪
Nu
=< br>0.5Re
1
l
2
⎛
⎜
µ
l
⎞⎜
⎟
⎛
⎜
R
4
K
⎞
⎟
⎪⎝
µ
g
⎟
⎠
⎜
⎝
Sp
'
⎟< br>⎠
⎪
C
⎪
⎪
Sp'
=
p
g
∆
T
sup
⎨
(
h
lg
+
0.5C
pg
∆
T
sup
)
Pr
g

⎪
⎪
K
=
ρ
(13)
l
⎪ρ
g
⎪
⎪
R
=
⎛
⎜
µ
gρ
g
⎞
⎪
⎩
⎜
⎝
µ
l
ρ⎟
l
⎟
⎠
式中，努塞尔特数
Nu=hdk
g
； 液相雷诺数
Re
l
=
ρ
l
U
∞
dµl
；汽相普朗特数
Pr
g
=C
pg
µ
g
k
g
；
d
为球体
直径，
m
；
h
为表面换热系数，
W(m
2
·

K)
；
U
∞
为
来流速度，
ms
；
∆
T
sup
为球 体过热度，
K
。
对
500

℃的球体在不同来流速度(1
，
1.5
，
2
，
3
，
4 ms)
的饱和冷却剂中膜态沸腾的过程进行数
值模拟，并用实验关联式
(13)
对模 拟结果进行对
比验证。
在计算中，使用了二维柱坐标，网格分辨率
为
120 ×120
，计算区域为
0.3×0.15 m
。由于蒸汽
膜很薄，在球面附近 布置了较多的网格以精确地
捕捉汽膜的位置。计算网格和边界条件见图
2
。
由 于计算中采用了定网格，所以球体保持静止，

92
核 动 力 工 程 Vol. 29. No.3. 2008

图 2 计算网格及边界条件
Fig. 2 Grid for Calculation and Boundary Conditions



t=5.0×10
－
3
s t=9.18×10
－
3
s t=1.61×10
－
2
s
t=3.11×10
－
2
s


t=4.73×10
－
2
s t=5.4×10
－
2
s t=6.2×10
－
2
s
t=6.7×10
－
2
s
图3 汽-液交界面随时间的变化过程
(来流速度=1.0 ms)
Fig. 3 Change Process of Vapor-Liquid Interface
with Time (U
∞
=1.0ms)


流体从下到上流过球体，这与球体在静止的流体
中运动的情形是等效的。

计算中，水的物理性质取为：

ρ
l
= 958.0 kgm
3
；
µ
l
= 2.79
×
10
－
4
kg(m · s)
k
l
= 0.68 W(m · K)
；
C
pl
= 4 217 J(kg · K)

蒸汽物理性质取汽膜平均温度所对应的物
理性质：

ρ
g
= 0.4 kgm
3
；
µ
g
= 1.97
×
10
－
5
kg(m · s)
；

k
g
= 0.04 W(m · K)
；
C
pg
= 2 015 J(kg · K)
其他参数为：表面张力
σ

= 5.8 9
×
10
－
2
Nm
，

饱和温度
T
s

=373

K
，相变潜热
h
lg

=2.257
×

10
6
Jkg
，重力
g=9.8 ms
2
，球体表面发射率
ε

=0.7
。

计算初始化时液面距离球心
0.01 m
。当来流速度
为
1.0 m s
时，计算得到的汽液自由界面随时间的
变化如图
3
所示。当
t=6.7
×
10
－
2
s
时，球体前端
蒸汽膜 及液体内的速度分布见图
4
。图中，蒸膜
内的速度分布呈抛物线形状，与文献
[6
、
7]
中提
出的“由于蒸汽膜内粘性力占主导地位，故速度
分布 沿径向为抛物线形状”相符。在不同来流速
度下，计算得到的球体表面对流换热平均努塞尔
特数 见图
5
。由图
5
可以看出，数值计算结果与
实验关联式
(1 3)
的计算结果较为吻合。



图4 球体前端蒸汽膜及液体内的速度分布
(t=6.7×10
－
2

s，来流速度=1.0 ms)
Fig. 4 Velocity Field at the Front of the Sphere in
Vapor and Liquid when t=6.7×10
-2
s and
U
∞
= 1.0 ms)


图5 不同来流速度下的对流换热平均努塞尔特数
Fig. 5 Average Nusselt Number at Different
Free Flow Velocity

袁明豪等：球体表面强制对流膜态沸腾换热的数值模拟

93
5 结 论
由
CFD
程序模拟的球体表面汽膜内的速度< br>分布与理论分析一致；对流换热努塞尔特数的计
算值与实验关联式符合较好。这表明，数值计算< br>较好地模拟了球体表面强制对流膜态沸腾换热的
物理过程。


参考文献：
[1] Theofanous T G, Yuen W W. The Prediction of Dynamic
Loads from Ex-Vessel Steam Explosions[C]. Proceed-
ings of the International Conference on New Trends in
Nuclear System Thermohydraulics, Pisa, 1994.
[2] Bromley L A, Leroy N R, Robbers J A. Heat Transfer in
Forced Convection Film Boiling[J], Ind Engng Chem,
1953, 45: 2639~2646.
[3] Kobayasi K. Film Boiling Heat Transfer around a Sphere
in Forced Convection[J], Nucl Sci, 1965, 2: 62~67.
[4] Stevens J W, Witte L C. Destabilization of Vapor Film
Boiling around Spheres[J]. Int J Heat Mass Transfer,
1973, 16: 669~679.
[5] Walford F J. Transient Heat Transfer from a Hot Nickel
Sphere Moving through Water [J]. Int J Heat Mass
Transfer, 1969, 12: 1621~1625.
[6] Dhir V K, Purohit G P. Subcooled Film-Boiling Heat
Transfer from Spheres[J]. Nuclear Engineering and De-
sign, 1978, 47: 49~66.



[7] Liu C, Theofanous T G. Film Boiling on Spheres in Sin-
gle and Two-Phase Flows[C]. Part 1: Experimental
Studies ANS Proceedings, Part 2: A Theoretical Study,
National Heat Transfer Conference, Portland, August，
1995.
[8] Kolev N I. Film Boiling on Vertical Plates and Spheres
[J]. Experimental Thermal and Fluid Science, 1998, 18:
97~115.
[9] Hirt C W, Nichols B D, Volume of Fluid(VOF) Method
for Dynamics of Free Boundaries[J]. J Comput Phys,
1981, 39: 201~225.
[10] Lafaurie B, Nardone C, Scardovelli R. Modelling Merg-
ing and Fragmentation in Multiphase Flows with
SURFER[J]. Comput Phys, 1994, 113:134~147.
[11] Ashgriz N, Poo J Y. FLAIR: Flux Line- Segment Model
for Advection and Interface Reconstruction[J]. J Com-
put. Phys, 1991, 93: 449~468.
[12] Youngs D L. Time-Dependent Multi-Material Flow
with Large Fluid Distortion[A]. Morton K W and
Baines M J. Numerical Methods for Fluid Dynamics.
New York: Academic, 1982. 273~285.

[13]
袁明豪, 杨燕华, 李天舒等. 基于VOF方法的带相
变的自由界面的计算[C]. 中国工程热物理学会多相
流学术年会，重庆，2006.
[14] Brackbill J, Kothe D B, Zemach C. A Continuum
Method for Modeling Surface Tension[J]. J Comput
Phys, 1992, 100 : 335~54.
Numerical Simulation of Forced Convection
Film Boiling on a Sphere

YUAN Ming-hao, YANG Yan-hua, LI Tian-shu, HU Zhi-hua

（School of Nuclear Science and Engineering, Shanghai Jiao Tong University, Shanghai, 200240, China）

Abstract
: A CFD code using improved Volume of Fluid (VOF) method to track liquid-vapor interface is
developed to simulate forced convection film boiling on a sphere. The simulation results are compared with
the experimental correlation, and the result show that the numerical method could simulate the physical proc-
ess of forced convection film boiling on a sphere successfully.
Key words:
Free surface, VOF, Film boiling, Phase change

作者简介：
袁明豪(1980— )，男，上海交通大学核科学与工程学院在读博士研究生。主要从事多相流动与传热数值的计算工作。
杨燕华(1962—)，女，教授。1996年毕业于日本东京大学核工程专业，获博士学位。主要从事多相流动 与传热、
反应堆热工水力、核电厂数值仿真及核电站严重事故的研究工作。
李天舒(1973 —)，男，上海交通大学核科学与工程学院在读博士研究生。主要从事核电站严重事故机理研究。

(责任编辑：尚作燕)