2有限元求解二维介质辐射——相变耦合传热

(6)
ρL
对(7)式进行有限元离散，得到
(7)
V
∫ ρ L ∂t Wdxdy = ∆t (∫ ρ L( f (t + ∆t ) − f (t ))Wdxdy)
l l V
∂f l
1
(8)

令 W = Γ i ，且有 fl (t + ∆t ) − fl (t ) = ∑ ( f l (t + ∆t ) − f l (t )) j Γ j ，由此可得：
kl / ks = 0.6 ， kmz / ks = 0.76 ， cl / cs = 1.2 ， cmz / cs = 1.12 ， L = 1000 J kg -1 ， cs = 1 J kg −1 K −1 ， ρ = 1000 kg m −3 ， ∆t = 0.01 s 。其中下标“ s ”代表固相区，“ mz ”代表 两相区，“ l ”代表液相区，方腔被划分为 894 个三角形单元。本文假设温度和相变不引 起介质折射率变化。
j =1
N
V
∫ ρL
N ∂f l ρL ( ∫ (∑ ( f l (t + ∆t ) − f l (t )) j Γ j )Γ i dxdy ) Wdxdy = ∂t ∆t V j =1
(9)
于是辐射—相变耦合换热能量方程(5)的有限元离散格式为：
− ∫ k ∇(∑ Ti Γ i )g ∇Γ j dV −
T * = (T − T0 ) /(Ti − T0 ) ， Ts* = 0.6 ， Tl * = 0.8 ， St = cs (Ti − T0 ) / L ， 无量纲时间 t * = α s t / Y 2 ， Y 为特征长度。下标 l 表示液相区参数， s 表示固相区参数， mz 表示两相区参数； Ti 表
中国工程热物理学会 学术会议论文
传热传质学 编号：093523
有限元求解二维介质辐射—相变耦合传热
易红亮 1，何凯 1，甄仌 1,2，谈和平 1
1 2
哈尔滨工业大学能源科学与工程学院，哈尔滨 150001 Tel: 0451-86412674, E-mail: yihongliang@
V i =1
N
上式可简写为：
([ M ] +
其中，
{P}t [C ] [C ] ){T }t = {N }t + {T }t −∆t + ∆t ∆t ∆t M ij = − k ∫ ∇Γ j g ∇Γi dV
V
(12)
N i = κ a ∫ (4σ T *4 − G )Γ i dV
V
Cij = ρ c p ∫ Γi Γ j dxdy
1
控制方程的有限元离散
半透明介质热辐射-相变耦合传热能量方程（焓方程）的微分形式为：
∂( ρ H ) = ∇g(k ∇T − qr ) ∂t
其中，
(1)
H = c pT + f l L
其中， H 为总焓， fl 为液相率， L 为潜热。液相率 fl 确定如下：
(2)
0 fl = ( H − H s ) /( H l − H s ) 1
V
Pij = ρ L( fl (t + ∆t ) − f l (t ))i ∫ Γ i Γ j dxdy
V
2
相变传热有限元离散模型的验证
本文结果与文献[22]给出的精确值进行比较以验证相变传热有限元模型的正确性。
可供比较的精确值对应的模型为一维模型，因此本文所建立的二维模型需假定方腔长宽 比为 10：1。参数设定为： kl / k s = 0.6 ， kmz / k s = 0.76 ， cl / cs = 1.2 ， cmz / cs = 1.12 ，
将(2) 式代入(1)式得：
H < Hs H s ≤ H ≤ Hl H > Hl
(3)
∂ ( ρ c pT ) ∂t
= ∇g( k ∇T − qr ) −
∂ ( ρ fl L) ∂t
(4)
假设 ρ 在整个区域内保持不变，且不随时间变化。在每一个时刻内， c p 保持上一个时刻 的值不变，于是忽略 c p 对时间的求导，则(4)式可以变为：
V i =1
N
N 1 ( ∫ ρ c p (∑ Ti Γi )t Γ j dxdy ) ∆t V i =1 N 1 ( ∫ ρ c p (∑ Ti Γi )t −∆t Γ j dxdy ) ∆t V i =1
= κ a ∫ (4σ T *4 − G )Γ j dV −
V
(10)
+
N j
ρL
∆t
( ∫ (∑ ( f l (t + ∆t ) − f l (t )) j Γ j )Γ i dxdy )
ρcp
∂f ∂T = ∇g(k ∇T ) − ρ L l − ∇gqr ∂t ∂t ∂f ∂T = ∇g(k ∇T ) − (∇gqr + ρ L l ) ∂t ∂t ∂f l f (t + ∆t ) − fl (t ) = ρL l ∂t ∆t
(5)
将相变潜热作为内热源项处理，
ρcp
时间采用向后差分，
y/Y
0.4 0.2 0 0
y/Y
0.4 0.2 0 0
0.2
0.4
x/X
0.6
0.8
1
0.2
0.4
x/X
0.6
0.8
1
(a) 无辐射 图3
(b)
n=2
t = 2 s 时热辐射对相变传热介质温度场的影响
图 2 描述了相变过程中热辐射所起的作用，通过观察可以发现，热辐射会影响相变 过程中各个时刻的温度分布，而且随着折射率的增大，这种影响逐渐加大。在相同的时 刻， 考虑热辐射作用的相变过程的温度低于纯相变过程的温度， 而且纯相变过程在第 919 步达到稳态，而折射率为 1 的辐射—相变耦合传热过程在第 906 步达到稳态，折射率为 2 的辐射—相变耦合传热过程在第 778 步达到稳态，由此可见热辐射加快了相变传热。 观察稳态温度分布发现，折射率为 2 的辐射—相变耦合传热过程的稳态温度高于其他两 种情况，说明辐射会影响相变过程的稳态温度，随着辐射作用的增大，稳态温度升高。 图 3 绘出了在 t = 2 s 时方腔内的等温线，可以看出在 t = 2 s 时，不考虑辐射的相变
哈尔滨商业大学土木与制冷学院，黑龙江 哈尔滨 150028
摘要：许多工程领域都涉及到半透明材料的高温相变传热问题，这类问题的实质是辐射-导热-相变瞬 态耦合传热问题。本文采用焓方程考虑导热、辐射及相变传热的耦合作用，将相变潜热项视为方程的 内热源项，采用有限元法求解导热、辐射及相变耦合传热的焓方程，并将不考虑辐射的相变传热结果 和文献提供的精确解进行比较，验证本文采用有限元法建立的导热-相变耦合模型的正确性。最后分析 了热辐射对半透明介质相变传热的影响。 关键词：半透明材料；辐射－相变耦合传热；焓方程；有限元
1.0 0.9 t=0.1s 0.8
T/Tr t=0.5s
t=1s
1.00 0.96 0.92
T/Tr
t=0.1s no radiation n=1 n=2 t=0.5s t=1s steady state
0.88 0.84 0.80 0.76 0.72 0.0 0.2 0.4
y/Y
0.7 0.6 0.5 0.0 0.2 n=no radiation n=1 steady state n=2 0.4
基金项目：国家自然科学基金 (No. 50806018)
Furmański [21]采用有限体积法求解焓方程和辐射传递方程，研究了半透明介质的二维瞬 态凝固过程。 本文将焓方程中的相变潜热项作为源项处理， 采用有限元法对瞬态辐射-导热耦合能 量方程进行离散，研究二维方腔半透明介质的凝固传热行为，并将纯凝固传热计算结果 和文献及精确解进行了比较。
0
前言
许多工程领域重要的物理过程和现象都涉及到半透明材料的高温相变（专指熔化和 凝固）问题。其典型工业应用包括：光学晶体生长[1-3]、半透明材料的加工、等离子镀 膜、重返大气层运载工具的烧蚀冷却、发电厂煤渣的凝固、考虑核反应堆安全时铀燃料 的熔化和凝固、高温热能存贮系统中相变材料的熔化和凝固、太阳能选择性涂层的制备 [4]以及高性能半导体薄膜材料的生产[5]。与不透明材料相比，半透明材料的相变机理和 行为更加复杂，因此大大增加了数值求解的难度。在分析中如果忽略热辐射，预测液相 份额和温度分布时，将导致严重偏差。 研究热辐射对半透明材料相变过程的影响始于上世纪 70 年代[6]。文献[6]假定存在 纯液态和纯固态两个相， 由平直界面隔开。 上世纪 80 年代， Chan [7]、 Hsu [8]、 Oruma [9]、 Webb [10]等学者研究了一维半透明材料的相变问题；文献[7]和[8]假设在纯液相和纯固 相区域之间存在等温“糊状区” 。1996 年至 2000 年，Yimer [11-13]采用焓法和球谐函数 法求解了相变—辐射的轴对称问题， 数值分析了半透明介质相变过程中辐射对温度分布、 固/液界面位置和能量存储能力的影响。2001 年，Rubtsov 和 Savvinova [14,15]研究了一 维经典 Stefan 问题（移动界面问题）中各向同性和各向异性散射对透明材料熔化和凝固 的影响[14]，还研究了各向异性散射对“糊状区”边界位置的影响[15]。2002 年，Yao 等人[16]拓展了早期 Chan 等人的工作[7,8]，采用焓法和离散坐标法对纯半透明材料受到 辐射和对流冷却时的一维 “糊状区” 平衡凝固过程进行了研究。 最近三年， Mishra 等[17,20] 连续发表多篇论文，对半透明材料的一维[17,18]、二维[19]及三维[20]凝固问题进行了研 究。假设半透明材料的凝固在某个温度范围内发生，因此相变发生时，具有明确的液相、 糊状和固相区。采用格子-波尔兹曼法[17,20]求解焓能量方程，采用离散传递法[17]、离 散坐标法[18]及有限体积法[19, 20]求解能量方程中的辐射信息。分析了衰减系数、散射 反照率、折射率等参数对介质内温度和液相或固相份额分布的影响。2008 年，Łapka 和