（工程热物理专业论文）多重相变问题和金属熔滴快速凝固过程的双倒易边界元数值模拟.pdf

中国科学技术大学 摘要 摘要 本文研究并发展了二维非定常导热并带有多重相变运动边界问题的双倒易 边界元方法。对一类热储能系统中，相变材料在周期性加热冷却边界条件下出现 的多重相变运动边界问题进行了模拟，得到了系统内部瞬态温度分布及相变界面 随时间推进的变化图象。并分析了一些因素的影响，对这类热储能系统的开发与 设计提供了有一定参考价值的科技资料。 另外，本文发展了同时具有辐射和型逾望昼釜件的袒銮堡垫整堡的轴对称双 倒易边界元方法，并成功地运用于微重力落管和落塔中及喷射成形两种过程中金 属熔滴的快速凝固过程的计算中。得到了过冷度，再辉时间金属温度随时间的 变化以及相变界面随时问的变化，并研究了金属滴大小，环境温度以及辐射率等 因素的影响。 关键词：双倒易边界完务法，多重宿菱，快遗凝固 里型堂垫查奎兰一竺! 塑! l a b s t r a c t g u ow e n ( e n g in e e ri n gt h e r m o p h y s ic s ) d ir e c t e db yl uw e n 一0 i a n g int h isp a p e r ，t h ed u a if e e i p r o e i t yb o u n d a r ye i e m e n tm e t h o d ( d r b e m ) t oc a lc u i a t et w o d i m e n s l o n a im u i t i p i ep h a s e c h a n g e m o v in g b o u n d a r y p r o b i e m s isd a v e i o p e d m e i t in ga n df r e e z i n gp r o b i e m si n v o iv in gm u i t i p i e p h a s e c h a n g em o v in gb o u n d a rie so fp h a s ec h a n g em a t e ria is ( p c m ) int h e r m ai s t o r a g es y s t e m su n d e rp e r i o d i c a iiy h e a t in g c o o ii n gb o u n d a r yc o n d i t i o n a r em o d e i e d a n dt h et r a n s i e n tt e m p e r a t u r ed is tr i b u t i o r sa n dp h a s ef r o n t tim e - m a r c hin gim a g e sinp c ma r eo b t a in e da n dt h ee f f e c t so fs o m ef a c t o r s a r es t u d ie d t h jss h o w st h a td r b e mc a nb e t t e rs oiv et h e s ep r o bie m s a n d t h isn u m e r i c a i a n a i y s is isu s e f u if o rt h eo p t i m iz e dd e s i g no ft h e r m a i s t o r a g es y s t e m s t h ea x is y m m e tr i cd u a i r e c i p r o c i t yb o u n d a r y e i e m e n tm e t h o d is d e v e i o p e d t oc a i c u i a t e p h a s e c h a n g e h e a ttr a n s f e r p r o b i e m s w i t ht h e r a d i a t i v e i ya n dc o n v e c t i v eb o u n d a r yc o n d i t i o n s t h er a p ids o ii d i f i c a t i o n o fm e t a id r o p i e t s int h es p r a yf o r m i n gp r o c e s sa n dm i c r o g r a v i t yf a ili n g t u b eo rt o w e ra r es t u die d t h eu n d e r c o o iin gt e m p e r a t u r e t h etim et o r e c aie s c e n c e 。t e m p e r a t u r eh is t o r yo ft h ed r o pie t p h a s e c h a n g ein t e r f a c e m o v i n gh is t o r ya r ea is oo b t a i n e d t h ee f f e c t so fs o m ef a c t o r ss u c ha st h e s iz eo fd r o p i e t e n v ir o n m e n tt e m p e r a t u r e e m is s i v i t yo nt h isp r o c e s sa r e s t u d i e d k e yw o r d s ：d u a ir e c i p r o c i t yb o u n d a r ye i e m e n tm e t h o d ，m u i t i p i e p h a s e c h a n g e r a p i ds o ii d i f i c a t i o n i i 主垦型堂垫查奎兰堡主堂垡望苎 一量二! ! 一 第一章引言 1 1 研究意义 热储能系统是利用潜热进行能量储存的系统，因而储热密度高，在太阳能领 域、蓄冷空调领域以及航天器件温控方面获得了广泛应用。了解热储能系统中 相变材料的相变传热过程对于此类热储能系统的设计具有重要意义。在一些情况 下，热储能系统中的相变材料可能会出现多个相变运动边界。例如，相变材料被 加热，出现液化继而进一步发生汽化，会出现多个相变运动边界；在微重力环境 下，为了最大限度地利用能源，热储能系统的充放冷周期应与空间站运行于向阳 轨道和地球阴影轨道的交替周期完全相同。1 。在这种情况下，热储能系统中的相 变材料在周期性的加热冷却过程中可能会出现多个相变运动边界。在这些情况 下，相变材料的传热过程与只有单一相变界面时有所不同。了解不同情况下的传 热过程，及影响相变运动界面个数、运动情况的因素，有助于更合理地设计和控 制这类热储能系统。 金属熔滴在快速冷却的条件下可以生成高性能的金属结构，因此已成为材料 加工的一种新型方法。与传统的加工工艺相比，金属熔滴在高冷却率下能够获得 亚稳的非平衡状态，这在传统的金属加工工艺中是无法获得的，理论和实验表明 3 j - s ，小尺度的液滴在快速冷却的条件下。能够得到均匀的晶体结构以及极好的 晶体性质。在快速凝固的方法中，主要有金属雾状微滴的喷射成形( s p r a y f o r m in g ) 和金属熔滴在微重力落管和微重力落塔中下落，它们都是将金属微滴加 热至较高的温度( 超过熔点) ，然后通过与环境高温差的辐射换热和对流换热实 现金属微滴的快速冷却。无论在微重力落管和落塔中的下落过程还是在喷射成形 过程中，由于实验一般只能测量部分温度及再辉的时间，缺乏对金属熔滴温度变 化以及相变过程的精确了解，而数值方法是模拟多种工况的简便且经济实用的途 径，因此在这一领域的研究中，数值方法具有非常重要的意义。 在相变传热问题中，为了准确捕获相变界面，一般用相变过程的物理条件， ， 如s t a f e n 条件来计算相变界面，因此数值解是解决这些问题的一个有效途径。 ! 堕型兰垫查查兰堡主兰堡堡苎一蔓主：! l 有限差分法在执行这类计算时需要反复迭代计算相变界面的贴体坐标，并且需 要在整个区域进行计算，计算效率较差。有限元法也需要在整个区域进行计算， 计算量很大。边界元方法可以直接在运动边界和自由面上离散计算节点，并具有 降维的特点，很适合求解自由面和运动边界问题。但是传统的边界元方法在求解 这类问题时需要附加上特殊的初边界条件，才能将区域积分转化为边界积分 “h “，在一般工程问题中的初边界条件下，不能充分发挥边界元方法的纯边界计 算效益。近年来发展的双倒易边界元方法是一种纯边界积分方法”1 ，可以直接在 边界上离散求解，大大减少了计算量并且已经成功地应用于单一相变运动边界 的问题“”。 1 2 研究概况 i 2 1 多重相变的研究概况 、 对于单一的相变运动边界问题，有很多这方面的研究“。”。 国内外对于多重相变运动边界问题的研究相对很少。 有限差分法( f d m ) ，有限元法( f e m ) ，边界元法( b e m ) 均进行了一维多重相变 问题的计算。z e r r o u k a t 和c h a t w in 使用了针对导热方程的显指数有限差分解， 并使用了根据运动边界速度而变化的网格”。b o n n e r o r 和j a m e t 扩展了适应性 空间一时间有限元方法求解这类多重相变问题。”，他们在相变界面上使用了 曲边三角形元，而在其他部分采用了梯形元。z e r r o u k a t 和w r o b e i 采用了边界 元方法计算这种维多重相变问题“。 对于热储能系统中的相变材料由于受到周期性的加热和冷却而出观的多重 相变问题，v i c k 等人进行了数值模拟工作。v i c k 等人首先针对热储能系统中的 多重相变问题建立了模型并发展了近似解，但是模型假设内部的温度分布均匀， 与边界相同。”。x y u 和n e is o n 对于一维圆柱系统中的多重相变问题进行了模拟 ”。”1 ，v i c k 和n e is o n 又进一步同时考虑了径向和轴向的导热，即考虑了热储能 系统沿着轴向的导热，但是在不同的轴向单元上仍作为一维的径向问题处理”“。 1 2 2 金属熔滴快速凝固过程的研究概况 金属熔滴的快速凝固过程得到了广泛的研究和关注。”。”1 。 中国科学技术大学硕士学位论文 第一章 p o u i i k a k u s 和w a id v o g e l 做了关于喷射成形中输运现象的广泛研究。”，但 是我们对于输运现象的效果的了解仍十分有限，因为大多数的研究集中于材料科 学方面，并主要集中在过程的亚稳特性上。例如，对于单个金属熔滴冷却的传热 分析通常忽略了金属滴内部的热传导。”。3 。b a y a z i t i t o g iu 和c e r n y 对于这种近 似模型进行了修正，他们使用了轴对称的热传导模型，他们的研究表明，对于相 对较低的冷冻率( 低于i 0 4 k s ) ，忽略金属滴内部导热得到的结果足够精确， 但是他们并没有考虑在金属熔滴深过冷过程中出现的非平衡凝固现象”“。 m e g a r i d is 使用了个复杂的同时考虑流体力学和对流换热物理模型，高温 的金属熔滴处于均匀的气流中，但是没有考虑辐射换热，也没有涉及到非平衡凝 固过程的计算1 。 k e h c h j nc h a n g 和c h in m in gc h e n 用有限体积法进行了数值模拟，同时考 虑了外部的散热及金属内部的热传导和非平衡相变过程，并比较了考虑内部传热 和忽略内部传热的不同结果”。 1 2 3 双倒易边界元方法求解相变问题的发展概况 边界元方法可以直接在运动边界和自由面上离散计算节点并具有降维的特 点，很适合求解自由面和运动边界问题。 g r e e n 函数边界元方法对于一维和多维带有相变运动边界的导热问题有过 很多报道“_ 邶3 。 双倒易边界元方法的思想由n a r d in i 和b r e b b i a 于1 9 8 2 年首次提出“，近 年来该方法得到了很大的发展并且已经成功地应用于单相变运动边界的问 题中1 ”。 多年来，我们实验室进一步发展了双倒易边界元方法并延拓应用于相变传热 ”“伽- 1 4 - 1 7 ，热波传播嘲“，非傅立叶导热3 ，多孔介质流动和相变传热“”等问 题。 本文是在实验室上述研究基础上，进步将双倒易边界元方法发展到多重相 变问题的计算，并获得了二维多重相变问题的结果。这些还没有见到报道。 本文还将我们实验室发展的轴对称双倒易边界元方法进一步发展为能计算非线 性辐射对流边界条件的边界元方法。 ! 里型兰垫查奎堂堡主兰竺鲨苎j ：! ! 一 1 3 本文工作 1 研究并发展了求解多重相变问题的双倒易边界元方法。 2 数值模拟了一维及二维多重相变问题。对于热储能系统中相变材料在第 二类及第三类周期性加热冷却条件下出现的多重相变问题，得到了瞬态 温度分布及相变界面随时间推进的图象。 3 研究并发展了非线性边界条件t c h 变传热过程的轴对称双倒易边界元 方法。 4 数值模拟了金属熔滴在快速冷却条件下的深过冷快速凝固过程。分别研 究了在微重力落管和落塔中及喷射成形过程中金属熔滴的快速凝固过 程，得到了过冷度，温度变化及相变界面随时间的变化，并研究了一些 因素的影响。 生里型堂垫查查兰堡主兰垡堡苎一生j 兰一 第二章双倒易边界元方法 2 1 双倒易边界元方法的概述 在传统的边界元方法中，对于诸如泊松方程，非定常扩散方程和带有非线性 项的偏微分方程，在将区域积分转化为边界积分的过程中，总会留有无法转化的 n 个边界点 图2 1 双倒易边界元方法 通常的离散情况 于它可以直接在边界上离散节点 问题。 区域积分项，因而在离散时不可避免地要 进行区域离散处理，这也就丧失了边界元 方法的优越性。如果使用针对某一类方程 的基本解，虽然有可能达到将区域积分全 部转化为边界积分的目的，而且计算精度 也很高，但其通用性又受到了限制。 双倒易边界元方法运用双倒易原理， 将由于方程中非线性项的存在而产生的区 域积分全部转化为边界积分，因此大大节 省了计算时间和对计算机内存的要求。由 非常适合于求解自有面问题和带有相变界面的 2 2 双倒易原理 在传统的边界元方法中，由于方程中的时间项、体积力项、非线性项等的存 在，使得转化过程中残留有区域积分。对于方程 v 2 丁= 6( 2 1 ) 假设己知边界条件是r = 于在厂l 上，g = a t c 3 n = 彳f 2 上，而r l + r 2 = r 是 整个计算域的外边界，”为边界上的单位外法线向量。利用加权残数方法可以写 出： f ( v 2 7 一6 ) r + d q = r g 一虿) r + d r r r f ) g + d r ( 2 2 ) 中国科学技术大学硕士学位论文 第二章 其中q 为计算区域，丁为拉普拉斯方程的基本解，它满足： v 2r + = 一( 2 - 3 ) 上式中的为狄拉克6 函数，根据狄拉克6 函数的性质，即 j t ( v 2 r ) d q = i 、( 2 4 ) n 其中正为未知函数在加载点( 源点) 的值，此处为q 内任意一点的温度值。利用 上述性质，对( 2 2 ) 式左端分部积分两次，得到关于域内及边界上任一点的积分 方程 c 。i p + q d f + f q t d f = t b t + d q ( 2 5 ) r1 1n 上式中，与点的几何性质有关。从以上推导的结果可以看到，上式右端残留 有区域积分项，为把它转化为边界积分，我们可以将原方程的解写成一个齐次方 程的通解和一个非齐次方程的特解的和，即r = f + 于，其中f 满足v 2 f = 0 ， 而于满足v2 于= b ，然而，特解对有些方程来说可能很难找到，或者即使找到了 也很难处理，这里我们将采用双倒易原理来解决这个问题。令： + l b 兰口， ( 2 6 ) ，= 1 其中口，为未知系数，为近似函数，n 为边界节点数，三为内部节点数。双倒 易边界元方法的主要思想就是利用一系列局部特解c 来代替个特解，其总数为 边界节点与内部节点的和，且满足 守t i = 3 ( 2 - 7 ) 对于近似函数，的选择，方法本身并无特别的限制，只是要求由乃构成的 矩阵为可逆矩阵。 将( 2 - 7 ) 代入( 2 - 0 ) 得到： 将( 2 - 8 ) 式代入( 2 1 ) 得到： ( 2 - 8 ) 、ll，一l、星q m 川 j i 6 中国科学技术大学硕士学位论文 第二章 v z ，：艺。，p ：f ) 上式两端乘以基本解丁后在计算域内积分，得到 ( 2 - 9 ) 巾：丁p + 施：芝o 。巾2 j r + 施 ( 2 1 0 ) n j 2 ln 对每一离散节点( 包括边界节点和内点) ，将上式两端分部积分两次，并利 用狄拉克6 函数的性质，我们可以建立关于域内任一点的积分方程： 。z + j g + 册一p + 。订：兰q f q 元+ 卜+ e 刃一巧，订1 ( 2 - i i )c 。z + j g + 册一p + g 订= q 卜元+ 卜+ e 刃一巧，订i rf j 。1 rf ( 2 1 1 ) 式与( 2 - 2 ) 式相比，可以看到，区域积分已全部转化为边界积分。 2 3 二维双倒易边界元方法 在二维问题中，对于l a p ，a c e 方程 v 2 丁：娶+ 婴：0 的基本解为丁：i n 一1 ，g ：娶 苏。却2 r 4 锄 对于，函数的选取使用最多的是距离函数( 或径基函数) ，其形式如下 一= ( 2 1 2 ) 其中吩定义为节点i , j 之间的距离( 如图2 1 所示) 。实际计算中仅取前几项即可。 对于一股的二维问题可以取为厂= 1 + ，文献( 8 ) 用大量的算例说明，此种取法 对工程应用具有足够的精度。 取定乃后，根据( 2 7 ) 式可求得，对于二维情况有 ( 2 - 1 3 ) 相应地，有： 铲等= 等等+ 筹丙o y = 卜+ 筹壤南 c z 州， 品。 = 、 中国科学技术大学硕士学位论文第二章 其中轨、协分别为0 在x 、y 方向上的分量。 对( 2 - 11 ) 式的纯边界方程，将边界离散化为n 个单元，得到： 啊+ 弘1 册一弘k = l 归2 纠j = l 嘱+ 缸k = t 秘舡k = l 叫 2 nr i ln， 鬻：詈凇磷恶 g ( 孝) = m l ( ) g 十+ 巾2 ( 乒) g “lj 、。 其中瓦，瓦+ 。和q 。，q 。分别为第k 个边界单元上两端节点的厂和q 值。f 为 ：熹 ( 2 17 ) l f 2 j。 爱蚓 中：( f ) = 丢( 1 + f ) f 。 伽+ d f = ( 巾，瓦+ a 2 t 。，。) g d f1 、 “， ( 2 - 1 9 ) g 蹦r = 胁”呦。) ，+ d f f “7 r r i 砖= p 。q + d f ，磙= p ：g d f kr g ：= p 。丁+ d f ，赢= p ：丁d f l ( 2 2 0 ) ! 望塑堂兰塞查盔兰堡燮塑一 蔓三兰 则( 2 1 9 ) 式口 写成： j 约+ d f - - h , r , + ；瓦。f 护们幽心f q 。2 ” n j 对( 2 2 0 ) 式中的积分可以用g a u s s 四点积分公式进行数值积分，并考虑到丁 在节点处连续，而9 如果在不光滑边界上的节点处存在角点效应，即角部节点处 q 在相邻单元内有不同的方向，其值可能不连续，因此，g 的系数不做组合， g ，女q t 2 g ：吼+ g * 2 q 。而丁的系数可以组合，令疗m ：碍1 + k ，而 日u 。舟：+ 乩；注意到，疗。= 0 ，将其与c ，组合，并令h 。：疗。+ c ，如，其中 屯= 1 ，i = k ；颤= 0 ，i k ，则方程( 2 15 ) 可以写成离散化的形式： 荟n 风一荟n 瓯吼= 芝j f f i l q 曙n 帆毛一荟2 n 瓯知) 、( 2 2 2 ) 上式可写成矩阵的形式： h u g q = ( h 矗一g q ) 位 ( 2 2 3 ) 将( 2 6 ) 式写成矩阵形式b = f a ，因此有= f - 1 b ，代入( 2 2 3 ) 式有： h u g q = ( 1 硒一g q ) f 一b ( 2 2 4 ) 2 4 轴对称双倒易边界元方法 三维柱坐标下l a p l a c e 方程的基本解为 几击。可磊霜菰丽1i 丽 对于( 2 - 7 f ) 式，可以写为： 嘱+ * t - t * q ) 布= 善q ( c ，o 十弥，观) 嬲) ( 2 - 2 5 ) ( 2 2 6 ) 对于轴对称问题，d s = r j 船，d f ，r 表示所选旋转边界的长度，令： ! 曼型兰垫查查堂堡主兰垡笙奎一望兰：! l 则有 耻舟他， k ，：上旷d p ， 2 r e0 。 蜘卜丁q ) d s = 2 r r j ( t k 一一q k 2 ) r j d r r 进一步得到： 即丽m 卅卿，争等等荨笋脚镑：， + 芒焉鲁鼬，争 妒了丽零菰南菥妒嘉考m ，争 其中，m = 亍妄，k ，争是第一类椭圆积分，e ，争是 x ( r + ) 2 + ( z ，一z ，) 2 22 第二类椭圆积分，n ，n ：为边界上的方向导数。 取，函数为n 砌 ，：1 十三厶e ( m ，要) ( 2 - 2 9 ) c ，正+ 兰j 巩k 。d f 一兰i q r k k ：d f “n ，“r、( 2 3 0 ) ：芝吩fq 元+ 壹亿足，订一妻j q j r k k 2 订) ，= 1 k = l ri = 】l t， 采用常单元进行计算，即当i ，点位于同一个单元上时，取单元的两个端点 1 ，2 的坐标分别为1 、z 。和如、z ：，则点i 、，的坐标分别为 卜三“+ 屹) 卜如) ! 里型堂垫查查堂堡主堂垡堕茎 一翌：三兰一 博：端囊 ( 2 - 3 1 ) 其中孝为相对坐标，在1 点，= 一1 ；h i 点，毒= 0 ；在2 点，孝= 1 。 代入( 2 - 3 0 ) 式，得到 啊+ 善正步墨订一善吼乒k ：订 = 善卜+ 幺窆k = 1 9 k 订一妻k = l 歌磬k z 订j ，= l 、0k 令h h = l k ，i d fg ，= f k 2 r j d f ，( 2 - 3 2 ) 式可以写成矩阵形式 h t c , q = ( h t g q ) 口 ( 2 - 3 2 ) ( 2 - 3 3 ) 将( 2 6 ) 式写成矩阵形式b = f a ，因此有= f b ，代入( 2 3 3 ) 式有： h t g q = ( h 亍一g q ) f 一1 b ( 2 3 4 ) 2 5 非定常导热问题的双倒易边界元解法 对于非定常导热问题 v 2 t ：= 1 a t c a t 令6 = 寿詈，表示为矩阵的腻 b ；f ：一1 塑 c t o t 因为f 可逆，因此有 旺：1 f 一1 翌 c t a t 代入( 2 - 2 4 ) 及( 2 - 3 4 ) 的矩阵公式中， ( 2 - 3 5 ) ( 2 3 6 ) ( 2 - 3 7 ) 中国科学技术大学硕士学位论文 第二章 定义矩阵c 为 h t g q = 壶( 睡g 咖。詈，( 2 - 3 8 ) c = 一寺( 斫一g 蜘 ( 2 3 9 ) 则( 2 - 3 8 ) 式司写为 h t + c 娶：g q ( 2 - 4 0 ) 对上式中的时间导数进行有限差分近似： a i r ：土f t m 一一t m l( 2 4 1 ) 其中，a t 为时间步长。 假设t 和q 在两时i 3 层内线性分布，即 葚1 ：器q 黑e q q ： 融z ， q = ( 一瓯) ”+”。i 其中，吼，0 q 为0 1 之间的常数，通常取最= 0 5 ，0 q = 1 0 。 则( 2 - 4 0 ) 式可以化为： ( 古c + 吼h ) t m + l _ 巳g q 肿1 = 【吉c 一( 1 一吼) h t ”+ ( 1 0 q ) g q ”( 2 - 4 3 ) 上式中，矩阵c ，h ，g 由计算节点的几何位置确定，右端项中的t ”，q ” 是己知值，根据边界条件，可以计算出未知值t ”1 或q ”1 。 中国科学技术大学硕士学位论文 第三章 第三章多重相变问题的数值模拟 3 1 引言 热储能系统是利用潜热进行能量储存的系统，因而储热密度高，在太阳能领 域、蓄冷空调领域以及航天器件温控方面获得了广泛应用”1 。而管式热储能系统 是其中一种常用形式，它一般包含管子和围绕它的相变材料，相变材料与管道内 的流体进行热交换。充冷时，温度低于相变材料凝固点的冷流体通过管道，使管 道外的相变材料不断凝固，从而将冷量以潜热的形式储存在相变材料中；放冷时， 管道中流过温度高于相变材料熔点的热流体，使凝固的相变材料不断融化，将蓄 积的冷量释放出来。了解热储能系统中相变材料的相变传热过程对于此类热储能 系统的设计具有重要意义。 ， 热储能系统中的相变材料可能会出现一个或多个相变运动边界。已经有许多 文章报道了出现单一相变界面的运动边界的问题”h ”3 。但是在一些情况下，相变 材料会出现多重相变运动边界。例如，相变材料被加热，出现液化继而进一步发 生汽化，会出现多个相变运动边界；在微重力环境下，为了最大限度地利用能源， 热储能系统的充放冷周期应与空间站运行于向阳轨道和地球阴影轨道的交替周 期完全相同“3 。在这种情况下，热储能系统中的相变材料在周期性的加热冷却过 程中可能会出现多个相变运动边界。在这些情况下，相变材料的传热过程与只有 单一相变界面时有所不同。了解不同情况下的传热过程，及影响相变运动界面个 数、运动情况的因素，有助于更合理地设计和控制这类热储能系统。 3 2 一维多重相变问题 3 2 1 问题的描述 厚度为口的固体，初始温度给定，一侧( x = 口) 受到恒定热流加热，另一侧 ( x = o ) 为绝热，一段时间后，当固体受热侧的温度达到熔点乇时，将出现第一 个相变界面将整个区域分为固相和液相。当继续加热，相变材料被加热面上的温 里型兰茎查查兰堡主堂垡笙奎 一塑三! ! 一 度达到气化温度兀时，将出现第2 个相变界面，即液气界面。假设气体一出现就 被取走，即热流直接作用在气液界面上如图31 所示。 这个问题在文献“中分别由有限差分法，有限元法和传统边界元法计算 过，本文用双倒易边界元方法进行计算，并将计算结果与有限差分方法”“、有限 元方法“2 2 1 及传统的边界元方法。3 1 的结果进行比较。 图3 1 一维多重相变问题示意图 3 2 2 物理模型 控制方程 在每个区域内，都满足热传导方程 v2 丁：上塑肛so r ， c 。西 其中，c 。= _ 二l ，屯，矶，c 。分别为导热系数，密度和比热。 初始条件 假设整个区域在起始时温度均匀分布 ，= t o ，r = 0 边界条件 加热面上为恒热流条件 ( 3 - 1 ) ( 3 - 2 ) 4 生里型堂垫查查堂堡主兰堡垒壅墨三j ! 一 丑掣：q oo f f 。 d n i 一 婴l ：q 0 o n i ( 3 - 3 ) ( 3 - 4 ) 另- - l 贝r j 边界上的绝热边界条件 以乱= 。 ( 3 _ 5 ) 侧面为绝热边界条件 一以乱! - o ( 3 _ 6 ) 相变界面上为相变边界条件 r ( x l ，y l ，) = l ( 3 7 ) t ( x 2 ，y 2 ，) = 瓦 ( 3 8 ) 班= 一九幽c 3 i ( 。, v , ) 一s 型o t “。， ( 3 9 ) t o z v v = ”_ 孔 ( 3 - 1 0 ) 3 2 3 数值方法 运用上一章2 3 的二维双倒易边界元方法的基本原理及2 5 的针对非定常导 热方程的具体处理方法，可以求解上述的物理模型。 由于在每一时刻，相变界面的位置是未知的，因此要进行相变运动边界的迭 代。本文采用松弛迭代方法，利用对相变速度的迭代来确定下一时刻相变界面的 位置，具体步骤如下： 1用( 2 4 3 ) 式计算固相区域的温度场和边界热流，当右侧边界上的温度达到熔 点时，将出现固液相变界面，将整个区域分为两个计算区域，固相区域1 和液相区域2 。只有一个相变界面时的迭代过程如下： ! 璺型兰垫查查兰堕圭兰垡丝茎 一笙三! ! 一 a 设固液界面上一时刻的位置为x ? 一。假设相变界面的运动速度为( v ? ) 。， 则可以计算出这一时刻相变界面的位置为( x ，) 7 = x + ( v i ，r ) 7 a t b 对于1 区，右侧边界为( x l m ) ，温度为l ，对于2 区，左侧边界为( x ? ) 。， 温度为l ，则根据( 2 4 3 ) 式分别对固相区域和液相区域进行计算，可以 计算出相变界面左右两侧的热流。 c 根据计算出来的温度梯度，可以得到新的速度 c v 1 - 1 = 去c 丑孔一丸乳。， d 两次速度相比较，如果满足条件l ( v ? ) 一( v ? ) 7 i 1 ( v ? ) 7 i s ，则进入下 步，回到步骤a 。否则进行松弛迭代，( v ? ) “1 = ( 】一脚t ) ( v ? ) 7 + f d i ( v ? ) “。， 其中，为松弛因子，0 q 蔓1 ，f 是迭代步数，返回a 继续运算。 2 当液相区域的右边界温度达到汽化温度t 时，将出现第二条相变运动边界。 即液气相变界面。假设气体一出现就被取走，所以只需考虑固相区域和液相 区域。假设固液界面和液气界面上一时刻的位置分别为，则相 变速度的迭代具体过程如下： a 分别假设固液界面和液气界面的运动速度为( v f ) 。和( v ? ) ，则两个相变 界面的位置分别为( z l 卅) 7 = x ? 1 + ( 吖) 7x a t ，( j ? ) 。= j + ( v ? ) 。a t 1 3 对于1 区，右侧边界为( x ? ) 。，温度为l ，对于2 区，左侧边界为( x ? ) 7 ， 温度为l ，右侧边界为( x ? ) 。，温度为t 。根据双倒易边界乖方法，可 以计算出两个相变界面左右两侧的温度梯度。 c 根据计算出来的温度梯度，可以得到新的速度： 广1 = 击c 丑乱，一丑乱 嘲 ) “15 去( q 0 孔) ( 3 _ 1 3 ) v优k h 1 中国科学技术大学硕士学位论文 第三章 d 如果同时满足 则回到步骤 旧) 一( v 圳脉v 圳 s ， a 进行下一步计算， 怫) “1 一( v 圳爪v 圳 s 否则进行松弛迭代。 ( v ，) = ( 1 一出i ) ( v ，) + 甜l ( v f ) 川，( v ? ) “1 = ( 1 一2 ) ( v 罗) + 国2 ( v ? ) “。季e 中， q ，2 为松弛因子，0 。1 ，0 中国科学技术大学硕士学位论文 第三章 3 0 0 0 2 5 0 0 5 0 0 0 图3 4 右边界及中心处的温度随时间变化 0 3 2 5 小结、 由上述算例可以看出，双倒易边界元方法能够很好的处理一维多重相变运动 边界问题。通过在狭长的区域进行计算，双倒易边界元方法得到的二维问题的结 果能够很好地代表一维问题的结果。与有限差分法，有限元法，及传统的边界元 方法相比，双倒易边界元方法是一种纯边界积分方法，因而可以明显地减少计算 时间及对内存的要求。同时也验证了二维双倒易边界元方法求解多重相变问题的 可靠性。 3 3 热储能系统中的二维多重相变问题 3 3 1 问题的描述 ， 研究管式热储能系统中出现的多重相变问题。只考虑热储能系统中的一根管 子及周围的相变材料( 冰) 。相交材料与管道内的流体进行热交换，当管道中流体 的进c i 温度在相变材料的熔点瓦上下周期性变动时，管子外表面的温度也将在 乙上下变动，管道外的相变材料将会从管壁开始交替地出现熔解及凝固，并不 断地向前推进，从而利用相变材料的凝结和融化来储存和释放热能。为了优化设 计这类热储能系统，需要深入了解其内部的多重相变界面运动的情况与相变传热 过程。文献( 2 7 ) ( 2 8 ) 模拟了圆柱系统中的这种多重相变问题，文献( 2 7 ) 模拟了 加 侣 仲 空一芒暑芒adeol 中国科学技术大学硕士学位论文 第三章 径向一维问题，文献( 2 8 ) 同时考虑了径向和轴向的导热，但在每一个轴向部分仍 按照一维径向问题进行计算。本文对这种问题的二维情形进行了模拟，考虑了在 一个加热冷却周期内相变材料的温度分布和多个相变界面的运动情况。 考虑热储能系统处于微重力环境下，或者包含了小密度差的相变材料，或者 周期性加热冷却在小温度差范围，因此其中相变材料的自然对流可以忽略，只考 虑其中的导热。计算区域选择系统的一个轴向截面，分别考虑了两种形状的相变 材料，即在圆管上包围了矩形和圆形的相变材料，数值模拟了相变材料在一个加 热冷却周期内发生的多重相变运动边界问题。 3 3 2 物理模型 相变材料处于如图3 5 所示的区域中初始时为固相相变材料的内表面受到周 期性的加热冷却条件其余边界为绝热条件对称轴上满足对称条件。先受到加 热，当其内表面上的温度达到熔点l 时，将会出现一个相变界面，整个区域被 分为固相和液相，分别用区域1 和区域2 表示，相变界面随着进一步的加热而向 前推进。段时间后，相变材料又受到冷却，当区域2 的边界温度达到凝固点r 。 时，将会出现第二个相变界面，相变界面2 将区域2 分为两相，分别用区域2 和区域3 表示。本文分别计算了在一个加热冷却周期内，相变材料的内表面受到 第二类和第三类加热冷却边界条件的多重相变问题。 b 广竺竺坠b 厂竺竺 jf1 iiij ijj 1 q 2 0 5i q = 0q 2 0 il q 。0 。f ， f。l夕响刚f ii 一一一一、 fff 2 f 1 厂、ii ，、i 5 f f 言，一i t f f 畜i a 图3 5 模型示意图 控制方程假设在每个液相区域和固相区域内，热物性为常数。在区域k 内，控制方程为一般形式的热传导方程 ! 堡堕堕壁叟堕望塑堕主兰垡塑 笙三塞 百o t = g v 2 丁 肛j 2 ，( 3 - 1 4 ) 其中，g = 2 , k i p k c t 为相k 的热扩散系数，丸，n ，q 分别为相k 的导热 系数，密度和比热。 初始条件假设整个区域存起始时温度均匀分布 边界条件相变材料外边界上是绝热边界条件 五i o t ：o ， j ：0 或x = a 或y ：6 d 7 相变材料内边界上受到加热冷却边界条件 1 ) 第二类边界条件：受到恒定热流的加热和冷却 一五娶：q o 砌 ( 3 15 ) ( 3 - 1 6 ) o i x a 2 0 - y - r o ，0 ，s r ( 3 - 17 ) 一以茜一绋耶x s o - y r o ”r ”t ( 3 - 1 8 ) 2 ) 第三类边界条件：由于环境温度周期性的变化而被加热和冷茹 一兄石o t 了= h ( l r ) ，q x 口2 ，o 少s ，o ，“ ( 3 1 9 ) 一以詈= 矗( 一t 一，) ，q x - a 2 , o ，“ ，s “+ t ( 3 2 0 ) 对称轴上满足对称边界条件 ! o ；l = o ，o x d 1 ，y ：o ：口2 x 口y ：o(321)n 。 ”_ ， 相变界面上满足相变边乔条件 r ( x i ，j ，i ，) = 7 j ( 3 2 2 ) t ( x 2 ，y 2 ，f = ? p( 3 2 3 ) 醐v 1 n ) - 如( 荨：k 厂 ( 割o t 。坩 浯2 4 ) 纵v 2 啦也( 荨，i ( x ：, y z ) - - 如( 刳o t 。k 2 5 ) 其中，为相变材料与管壁接触面的半径，三表示相变潜热，v ；，( x 。，儿) 分别表示第k 个相变界面的运动速度和牮标。 2 】 中国科学技术大学硕士学位论文 第三章 3 3 3 数值方法 运用上一章2 3 的二维双倒易边界元方法的基本原理及2 5 的针对非定常导 热方程的具体处理方法，可以求解上述的物理模型。并采用松弛迭代方法，利用 对相变速度的迭代来确定下一时刻相变界面的位置，具体步骤如下： 1 处于加热周期，用( 2 4 3 ) 式计算固相区域1 的温度和边界热流，当区域内边 界上的温度达到熔点l ，时，将出现固液相变界面，将整个区域分为两个计算 区域，固相区域1 和液相区域2 。 a 设固液界面上一时刻的位置为x 。一。假设相变界面的运动速度为( v ? ) 。， 则可以计算出这一时刻相变界面的位置为( x ? ) 。= x r + ( v ? ) 7 a t 。 b 根据( 2 4 3 ) 式分别对固相区域1 和液相区域2 进行计算，则可以计算出 相变界面两侧的热流。 c 根据计算出来的温度梯度，可以得到新的速度 ( ，? l n ) “：喜晦坚喝娶】| ( 3 - 2 6 ) 肚1 d n 珈i ( x l , y j l d 两次速度相比较，如果满足条件 ( v ? n ) “1 一( v ? n ) f ( v ? - n ) 7 f s ，则 进入下一步，回到步骤 。否则进行松弛迭代， ( v ? n ) “1 = ( 1 一q ) ( v f - n ) 7 + l ( v 7 n ) “ 其中，脚。为松弛因子，0 f - o 。1 ，j 是迭代步数，返回a 继续运算。 2 冷却周期开始后，区域2 受到冷却温度下降，但相变界面继续向前推进，当 液相区域2 的内边界温度降低到凝固点咒时，将出现第二条固液相变运动边 界，将区域2 分为两部分，新的固相区域3 和液相区域2 。假设两条相变界 面上一时刻的位置分别为x 卅_ 1 ，x ，则迭代具体过程如下： a 分别假设两条相变界面的运动速度为( v ? ) 7 和( v ；) ，则两个相变界面的 位置分别为( x ? ) 。= x ? 。+ ( v ? ) 7 a t ，( x ；) 7 = x 十( v ：) 。a t 。 b 根据a 计算出的新的边界位置，分别对于1 区，2 区及3 区进行计算， 可以得到两个相变界面两侧的温度梯度。 中国科学技术大学硕士学位论文 第三章 c 根据计算出采的温厦梯厦，口j 以得到耕阴速度： ( v ? “= 瓦1 ，心面o t _ 詈也 ( v ；n ) f + 1 2 瓦1 。瞧面a t t ：i o t 扎。 d 如果同时满足 阳n ) 一( v ? 圳巾v ? n ) 仆s i ( v ；n ) “1 一( v ；n ) v i ( v ；- n ) 仆s 则进行下一步计算，否则进行松弛迭代。 ( v ? - n ) “1 = ( 1 一c 0 1 ) ( v ? n ) + 埘l ( v 7 n ) 卜 ( v 7 - n ) “1 = ( 1 一珊2 ) ( v ；n ) 。+ 2 ( v ；n ) “! ( 3 2 - ) ( 3 2 8 ) 其中，c 0 1 ，( - d 2 为松弛因子，0 q l ，0 甜2 】，是迭代步数，返回 a 继续运算。 、 3 3 4 数值结果与讨论 在3 4 1 提出的物理模型中，以冰( 水) 作为相变材料进行计算。 冰( 水) 的热物性参数如下： l 。= l 2 = 3 3 5 o k j k g ，t 。= 一2 。c ，瓦= 0 = o o c ，c a = 龟= 5 2 7 ( d j ( 堙固， c 2 = 4 2 2 0 0 j ( k g k ) ，p j = 岛= 9 2 0 k g m3 ，岛= 1 0 0 ( k g m 3 ， = = 2 2 6 w ( m 固， = o 5 6 0 v ( m 固。 计算区域为l c m 2 c m 的矩形区域，内边界为半径为o 1 c m 的圆。 a = l c m ，b = 2 c m d 1 = 0 4 c m ，a 2 = 0 6 c m ，f o = o 1 c m 。在上边界上取了1 0 个 相等的边界单元，a x = o 1c ，左右边界上各取了1 0 个边界单元，a y = 0 2 c l ， 在热流边界上取了1 0 个单元，a x = 0 o l c m ，对称轴上取了8 个单元，在区域内 部取了8 1 个内点。采用线性单元进行计算。时间步长a t = 1 o j 。 下面分别讨论相变材料的内边界在第二类边界条件( 即受到恒定热流的加热 和冷却) 和第三类边界条件( 即由于环境温度的周期性改变而被加热和冷却) 下 的数值结果。 中国科学技术大学硕士学位论文 第三章 3 3 4 1 第二类边界条件 内表面的边界条件如下 一k l ( 婴) 。：q 。，o f k t ：擘) ：q o ，m - t “ 一