（人机与环境工程专业论文）结冰融冰过程的数值模拟.pdf

结冰融冰过程的数值模拟 a b s t r a c t t h ep h e n o m e n ao fi c i n ga n dm e l t i n ga r el a r g er e l a t e dw i t ho u rd a i l yl i f e t h e n u m e r i c a ls i m u l a t i o no ft h et w os i t u a t i o n sh a v e b e e nah o t s p o ti nt h ef i l e do f n u m e r i c a lh e a tt r a n s f e r i nt h i sp a p e r , b a s e do ne n t h a l p ym e t h o di nf i x e dg r i d s ，a n u m e r i c a ls i m u l a t i o no fm e l t i n gp r o c e s so na i r f o i la n dan u m e r i c a ls t u d yo fi c i n g c o n s i d e r i n gn a t u r a lc o n v e c t i o ni nac a v i t yi sp r e s e n t e dr e s p e c t i v e l y f i r s t ，t h ed i s s e r t a t i o ni n t r o d u c e st h ep r i n c i p l eo fi c i n ga n dt h ei m p o r t a n c ea n d a c t u a l i t i e si nr e s e a r c h i n gi c i n ga n d m e l t i n g t h e n ，t h em a i na p p r o a c h e st on u m e r i c a l s i m u l a t et h ep h e n o m e n aa r es u m m a r i z e d t h es e n s i b l eh e a tc a p a c i t ym e t h o dm o d e l a n dt h ee n t h a l p ym e t h o dm o d e l ，w h i c ha r et h et w op r i m a r ym a t h e m a t i c a lm o d e l st o s t u d yp h a s ec h a n g eh e a tt r a n s f e r , w e r eb o t hd e t a i l e da n a l y z e d t h es e l f - c o n t a i n e d m a t h e m a t i c a lm o d e lf o rn a t u r a lc o n v e c t i o np h a s ec h a n g eh e a tt r a n s f e ri se s t a b l i s h e d t h eg r i dg e n e r a t i o nt e c h n i q u e ，d i f f e r e n c es c h e m eo fc o n v e c t i v ea n dd i f f u s i v et e r m s ， p r e s s u r ea n dv e l o c i t yc o r r e c t i o nm e t h o d sa n da r i t h m e t i co fn o n l i n e a re q u a t i o n sa r e d e t e r m i n e d i nn u m e r i c a ls i m u l a t i o no ft h r e ed i f f e r e n tk i n d s h a p eo fi c em e l t i n g p r o c e s so na i r f o i l ，t h en a t u r a lc o n v e c t i o nt e r mi ng o v e r n i n ge q u a t i o ni ss i m p l i f i e d o n n u m e r i c a ls i m u l a t et h ef r e e z i n go fw a t e ri nac a v i t y , t h en a t u r a lc o n v e c t i o ni ss t u d i e d b e s i d e s ，t h ee f f e c to f t h ed e n s i t yi nt h eg r i d sa n dt h es p a c eo ft h e c a l c u l a t i n g t i m ea r e r e s e a r c h e d t h er e s u l t so b t a i n e ds h o wa v e r yg o o da g r e e m e n tb e t w e e nn u m e r i c a la n d e x p e r i m e n t a lr e s u l t s a n a l y s i ss h o w st h a tt h ee f f e c to fn a t u r a lc o n v e c t i o nc a nn o tb e n e g l e c t e d k e yw o r d s ：i c i n g ，m e l t i n g ，n u m e r i c a ls i m u l a t i o n ，p h a s ec h a n g eh e a tt r a n s f e r e n t h a l p ym e t h o d ，n a t u r a lc o n v e c t i o n i i 承诺书 本人郑重声明：所呈交的学位论文，是本人在导师指导下，独立 进行研究工作所取得的成果。尽我所知，除文中已经注明引用的内容 外，本学位论文的研究成果不包含任何他人享有著作权的内容。对本 论文所涉及的研究工作做出贡献的其他个人和集体，均已在文中以明 确方式标明。 本人授权南京航空航天大学可以有权保留送交论文的复印件，允 许论文被查阅和借阅，可以将学位论文的全部或部分内容编入有关数 据库进行检索，可以采用影印、缩印或其他复制手段保存论文。 ( 保密的学位论文在解密后适用本承诺书) 作者签名：丝翌篁 日 期：边丛经瑚 南京航空航天大学硕士学位论文 符号表 英文字母 口热扩散系数 4系数矩阵 巳定压比热 ，液相率 g重力加速度，固相率 h比焓 h焓 ，雅可比因子 三相变潜热 n u努赛尔数 p r普朗特数 r e雷诺数 s t e斯蒂芬数 s源项 丁温度 t时间 “x 方向速度 u 方向速度，逆变速度 v y 方向速度 vt 1 方向速度，逆变速度 p压力 p产生项 q坐标变化参数 m 一 甩 工，y 质量，矩阵 法向向量 矩阵 直角坐标 希腊字母 强 自然对流换热系数 p网格参数 y网格参数 导热系数 p密度 y 运动粘度 毒，玎曲线坐标 求解变量 妒源项函数 甲源项函数 r疏运系数 v 结冰融冰过程的数值模拟 上标 v i 上一时层的迭代收敛值 迭代次数 预估值 东界面 东邻点 南界面固相 南邻点 西界面 西邻点 北界面 北邻点 液相 参考值 初始值 坐标方向，网格节点 +阼 o ” 标 ， 打 ， 雕 。e ，s w 。，阿：萎n 南京航空航天大学硕士学位论文 第一章绪论 在自然界和诸多工业生产、工程领域中，存在着大量的结冰、融冰现象。 例如：地球两极极地冰块的融化和凝固；冰蓄冷空调系统中蓄冰槽内的蓄冰、 融冰；飞机迎风面在飞行中的结冰及除冰；食品的冷冻冷藏：河水的冰封和解 冻等等。 自从1 8 9 1 年，科学家j s t e f a n 在关于极地冰层融化的研究中首次讨论了结 冰、融冰这一课题以来，人们就对此问题开展了研究c l 】【2 j 【3 】。研究的主要内容有： 结冰、融冰的原理；外部环境的变化对结冰、融冰的影响；结冰、融冰空间内 的温度分布；液相中存在自然对流或强制对流对结冰的影响；结冰、融冰过程 中相变界面的移动规律等。 结冰融冰问题的研究对提高节能技术，避免飞机结冰以保证其飞行的安全 性，延长食品的保存期，深入了解地球科学等诸多生产生活中存在的实际问题 都具有重要意义。 1 1 结冰原理 水要结成冰，首先水的温度降到0 以下，成为过冷水。从热力学观点来看， 过冷水处于亚稳态，该状态的解除需要大于临界尺寸的冰核的形成。冰核的形 成存在两种机理：水体内的均匀成核与固液相变界面上的非均匀成核。当过冷 水中出现尺寸大于临界尺寸的冰核时，结冰过程开始，冰核在过冷水中长大， 最终成为宏观意义上的冰。同时，冰核的形成是一种随机现象，即冰核形成时 水的过冷度为一定范围内的概率分布。过冷度的大小对是否发生结冰现象有着 重要的影响。一定过冷度的过冷水是否发生结冰受到许多因素的影响【4 】 ”，其中 可能包括水的流动状态，结冰基体表面几何特性、物理特性、表面面积以及外 加的作用力等等。 1 1 1 结冰过程 由于液体与晶体的结构不同，同一物质的这两种状态在不同温度下的自由 能变化是不同的，如图1 1 所示，因此它们便会在一定的温度下出现一个平衡点， 即理论结晶温度兀。低于理论结晶温度，由于液相的自由能e 。高于固相晶体的 结冰融冰过程的数值模拟 自由能c 。液体向晶体的转变会使能量降低，于是发生结晶。也就是说，要使 液体发结晶，必须造成液体与晶体之间的自由能差( 廿= 一) ，即具有一 定的结晶推动力才行。实际结晶温度正与理论结晶温度兀之间的温差叫过冷度 ( a t = r r o ) 。液体的冷却速度越大，过冷度便越大；而过冷度越大，自由能 差舻便越大，即结晶推动力越 大，结晶倾向越大 6 】。 结晶时同时要发生两方面 的能量变化：一方面时因晶体的。 产生而引起体系自有能的降低：曼 另一方面时因晶体的出现，增加 新的界面而引起的体系自由能 的升高。因此，要使液体真正能 够结晶，不仅要有过冷造成晶体 与液体之间的自有能差以提供 结晶的推动力，而且还必须造成 足够的过冷度以使其自由能的 降低超过表面能的增加，只有这 样才能满足物质自高能态向低 能态变化的规律。这就是结晶必 须要有过冷的原因。 图i 1 液体与晶体在不同的温度下自由能的变化 水在持续降温并最终结冰的过程中，主要经历了以下几个阶段：( 1 ) 液态 水在被冷却到结晶温度o 的过程中，随时不断地产生尺寸较小、大小不一且极 不稳定的晶胚，它也是随后产生的晶核的来源；( 2 ) 当水被过冷至0 以下，某 些具有较大尺寸因而比较稳定的晶胚真正成长为晶核；( 3 ) 晶核渐渐成长，出 现棱角，棱角处的散热条件优于其他部位，因而优先成长，最易把晶间填满。 这种成长方式叫枝晶成长。过冷度越大，枝晶成长的特点便越显著：( 4 ) 大量 晶核成长达到相互接触以后，冰层形成。 1 1 2 冰核形成机理 在各种固体表面上形成的冰核能否长大并最终导致冰的形成，要看冰核自 由能a g 随半径，的变化情况，若冰核自由能a g 随半径，的增大而减小，则冰核 可以继续长大；若冰核自由能a g 随半径，的增大而增大，则冰核的半径会不断 南京航空航天大学硕士学位论文 减小。对于在平面或球面上形成的冰核，存在一临界半径r ，在临界半径时冰 核的自由能最大，即只要冰核的半径大于这一半径，则稳定的冰核就可以形成。 临界半径与过冷度有关，过冷度越大，临界半径越小。 所谓均匀质核生成是在水体中形成冰核。文献 7 中记述了诸多学者对均匀 质核生成进行的理论和实验研究。实验包括将不同体积的纯水的水滴置于疏水 性的比水重的四氯化碳与比水轻的油之间进行冷却，测量水开始结冰时的温度。 实验者还从概率分析入手，从有限的实验数据算出其它情况下过冷水的结冰概 率。结冰概率分析基于这样一个事实：一旦在水体内某点生成了稳定的冰核， 整个水体都将解除过冷状态。 所谓非均匀成核是指在固液界面上形成冰核。s a t o r uc h i r i f u 8 】从单个冰核微 观形成机理入手，考察了固液界面上冰核的形成速度： j ：n z k _ _ f _ te x p f - 丝 e x p f _ 笪1( 1 - - i ) h lk t lk t 式中：月，一单位体积内的水分子数：k 一玻尔兹曼常数；t 一水的温度；a g 一水分子的活化能：a g 一冰核的自由能。 文献【9 对从上面冷却的水体在水体内的结冰概率进行模拟计算，并和传热 表面上的结冰概率进行比较后，证明了以下结论：对于通常的水体，水体内的 匀质成核与在壁面或水体内悬浮粒子表面上的非均匀成核相比，其概率是非常 小的，可以忽略。 1 1 3 过冷度的影响 2 。 褂雠 j 冀、“ 不同的过冷度t 对晶核的形姜藿 成率n ( 晶核形成数目s m m 3 ) 和谴谴 成长率( r a m s ) 的影响如图1 2 所示。 当过冷度等于零时，晶核的形成率 和成长率均为零。随着过冷度的增 加，晶核的形成率和成长率都增大， 并在一定的过冷度时达到最大值。 而后当过冷度进一步增大时，它们 又逐渐减小，直至在很大的过冷度 的情况下，两者又先后各趋于零。 过冷度对晶核的形成率和成长率的 图1 2 晶核的形成率n 和成长率g 与过冷度t 之间的关系 结冰融冰过程的数值模拟 这些影响，主要是因为在结晶过程中有两个相反的因素在同时其作用。其中之 一就是晶体与液体的自由能差 ，它是晶核形成和成长的推动力；另一相反因素 便是液体中原子迁移能力或扩散系数d ，这是晶核的形成和成长的必须条件，因 为原子的扩散系数太小的话，晶核的形成和成长同样也是难以进行的。如图1 3 所示，随着过冷度的增加，液体与晶体的自由能差便愈大，而液体中的原子扩 散系数却迅速减小。由于这两种随过冷度不同而作用相反变化的因素的综合作 用，便使晶核的形成率和成长率 与过冷度的关系上出现一个极大 l o 值。在过冷度较小时，虽然原子 的扩散系数较大，但因作为结晶最 推动力的自由能差较小，以致晶 核的形成率和成长率都较小；在 过冷度度较大时，虽然作为结晶 推动力的自由能差很大，但由于 原子的扩散在此情况下相当困 难，故晶核也难以形成和成长，f ；d 而只有两种因素在中等过冷度情图1 3 液体与晶体的自由能差a f 和扩散系数d 与 况下都不存在明显不利的的影响 过冷度a t 的关系 时，晶核的形成率和成长率才会 达到其极大值。因为过冷度对结冰过程有着重要的影响，国内外学者对影响过 冷度的因素进行了许多研究工作。这些因素包括与水体接触面的材料性质，在 水体中加入的成核添加剂成分与剂量，外加电场超声场等 1 0 。 对于与水体接触的材料表面的影响，有关学者等做过多次实验研究，分别 以机械抛光、电抛光、粗糙、镀金及镀镍铜表面为研究对象，发现过冷水开始 结冰的温度具有统计性质，它在不同温区的分布与表面性质及冷却速度有关。 表面越粗糙，冷却速率越低，则过冷度越小。后来研究的表面扩展至硅表面、 聚四氟乙烯表面及丙稀酸塑料表面。各种表面上的成核能力可以用一个参数t 描述，参数t 的物理意义时在单位表面上有效成核点的个数。其大小与表面粗 糙度有关。若抛光铜表面的t 值为1 ，则聚四氟乙烯表面的t 值为o - 3 ，而粗糙 铜表面的t 值为2 0 0 。 在水体中加入适量的成核添加剂是降低过冷度的有效方法。实验者对黄原 南京航空航天大学硕士学位论文 胶、碘化银、硫化铜、胆淄醇及成核菌等多种添加剂进行了实验研究，其中碘 化银效果最佳，因为它的微观晶体结构与冰晶类似，冰核更易在其表面上形成。 文献 1 2 考察了封装于圆柱单元内水的过冷现象，着重研究了单元内水的体 积及冷却速度对结冰概率的影响，得到以下公式： i n 一l n ( 1 一p ) 】= i n ( k a v 日) + 丑t ( 1 - - 2 ) 式中：p 一结冰概率；足一常数；v 一水体体积；一冷却速度；爿，b 一与水体种类及成核添加剂有关的常数；r 一过冷温度。 电场与超声波等外加因素对过冷水也有重要影响。根据实验研究，若将直 流电场以电极插入的形式施加于过冷水体。则水体在电压低于1 0 0 v 时即可解除 过冷态。直流电场对过冷水的作用为：在电极不光滑或突起处，附近电场强度 很高，极性水分子在电场作用下聚集，很容易形成冰核。超声波对结晶有促进 的作用，这一点早为人知。利用这一点，将超声波用于动态制冰中水过冷状态 的解除，沿着管道布置一系列超声波发生器，在超声波的作用下，水体在过冷 至0 3 c 时即可结冰。与其它动态制冰方式相比，超声波还可以有效防止已形成 的冰晶与管道内壁粘结。 由于结冰过程中水的过冷度大小影响到结冰时间的快慢，因此有必要采取 措施来降低水在结冰前的过冷度。 1 2 结冰融冰问题的研究背景及意义 1 2 i 节能领域的研究意义 能源是人类赖以生存的五大要素之一，是国民经济和社会发展的重要战略 物资。经济、能源与环境的协调发展，是实现中国现代化目标的重要前提。随 着全球工业的迅猛发展，能源问题越来越为人们所关注，能源问题也已经成为 一个全球性战略问题。目前，世界能源结构中所利用的化石能源主要仍然是煤 炭，其次才是石油和天然气，其比例约为6 8 ：1 7 ：1 5 。根据国际上通行的能源预测， 石油将在4 0 年时间内枯竭，天然气将在6 0 年内用光，煤炭也只能用2 2 0 年。 据预测，中国未来能源供需的缺口将越来越大，2 0 1 0 年约缺能8 ，到2 0 5 0 年 将短缺2 4 左右，其中石油缺额可能多达4 4 亿吨标煤。石油进口依存度( 净进 口量与消费量之比) 由1 9 9 5 年的6 6 上升为2 0 0 0 年的2 0 。预计2 0 1 0 年将上 升为2 3 。天然气进口依存度2 0 0 0 年为6 ，2 0 1 0 年为2 0 。另外能源安全性 的问题也将提到议事日程，据专家预测2 0 1 5 年前后将可能发生第三次世界性石 结冰融冰过程的数值模拟 油危机，届时我国油气进口依存度将会达到2 5 左右，在新的世界性石油危机 面前，我国能源供应的安全性将受到威胁。在人类利用能源的历史长河中，石 油、煤炭、天然气等常规能源的使用毕竟是短暂的一瞬间【1 4 】。 面对严峻的现实，一方面必须从现在开始及早寻找新一代的替代能源；另 一方面，在现有的资源配置条件下，必须采用先进技术推进节能，加速可再生 能源开发利用以及依靠市场力量优化资源配置。因此蓄能技术应运而生，蓄能 研究的重要性也日益突出。蓄能可以缓解能量供求双方在时间、强度及地点上 不匹配，是合理利用能源及减轻污染的有效途径，也是广义热能系统优化运行 的重要手段。相变蓄能以蓄能密度高、易与运行系统相匹配、易控制等优点， 日益成为蓄能系统的首选方式。近年来，相变蓄能系统已经越来越广泛地应用 于建筑节能、蓄冷空调和一些余热回收系统。 潜热蓄能之所以具有独特的魅力，主要有以下二个特征：第一，较高的蓄 能密度，蓄能潜力大；第二，相变物质( p c m ) 蓄能时保持恒定的相变温度， 从传热学的角度分析，较小的传热温差意味着有效地减少换热过程不可逆损失， 提高换热效率。 冰蓄冷就是将水制成冰来储存冷量，它是潜热蓄冷的一种方式。水是一种 相变潜热值很大的物质，另外水的比热也很大，而且冰作为蓄冷介质具有性能 稳定，安全无污染，大量易得，价格便宜，使用最为方便等优点。冰蓄冷技术 主要优点在于转移用电高峰，平衡电网负荷，利用用电峰谷的差价以降低制冷 系统的运行费用。因此，冰蓄冷是一种良好的蓄能装置。冰蓄冷技术在发达国 家已普遍应用达十几年之久，现在在我国也有较快的发展。 由于冰的单位体积蓄冷量比水大得多，水的比热容为4 1 8 k j k g ，而相变潜 热为3 3 5 k j k g ，因此，蓄冰槽的容积只有水蓄冷的2 3 0 l 5 ，这大大地减少了 蓄冰槽的体积，也极大地促进了冰蓄冷技术地推广和应用。较为常见的有冰盘 管式、完全冻结式、冰球式等技术方式。 目前冰蓄冷系统中所需解决的主要问题有：对强化制冰技术的研究；高效 低温送风末端设备的研究与开发等。研究蓄冰、融冰对提高冰蓄冷系统的性能， 进一步改进系统的设计方法和优化设备选型，实现系统的经济性运行都有极为 重要的意义吲。 1 2 2 飞机防冰系统中的研究意义 飞机在负温云层飞行，或具有负温表面的飞机在正温云层或无云大气中飞 南京航空航天大学硕士学位论文 行时，飞机部件迎风面会发生结冰现象。飞机升力表面结冰时，会引起飞机的 操纵性和稳定性品质恶化；发动机进气道及动力装置发生结冰时，会导致发动 机的损坏，停车等严重后果；环控系统冷风道结冰时，会导致环控系统的制冷 能力下降，引起机载电子设备的失效：有时，甚至少量的结冰也会导致机毁人 亡的悲剧发生。因此，为保证飞机飞行的安全性，开展飞机防、除冰系统方面 的研究非常重要。 飞机在飞行时与云层中的过冷水滴撞击，水滴撞击到结冰表面后会形成薄 的水膜或由于表面张力形成水滴，水滴有可能在气动力作用下被吹走，也可能 结成冰。在温度范围为- - 4 0 一o 时，所结的冰形状可分为霜状冰( r i m ei c e ) ( 又常被称为雾淞或结晶体冰) 和瘤状冰( g l a z ei c e ) ( 又常被称为明冰或玻璃 冰) ，或两者混合。水滴撞击到机翼表面后立即凝固时，就会产生霜状冰( 如图 1 4 所示) ，霜状冰往往贴附着机翼前沿表面生成，所以对气动性能影响较小。 水滴撞击机翼表面后没有立即完全凝结，而是形成冰水混合物，最后会形成瘤 状冰( 如图1 5 所示) ，瘤状冰会有棱角伸展到机翼的来流中，这会造成很大区 域的流体分离，严重的影响飞机的气动性能。除冰系统工作时，若一次除冰循 环不能把冰全部清除掉，还会产生残余冰( r e s i d u a li c e ) 和间隙冰( i n t e r - c y c l ei c e ) 。 残余冰( r e s i d u a li c e ) 是指除冰系统完成一次除冰操作后，翼面上还残余的 未被除去的冰。间隙冰( i n t e r - c y c l ei c e ) 则是指除冰系统完成一次除冰操作后的 残余冰，再加上除冰系统下一次操作之前新生成的冰。残余冰和间隙冰对飞机 气动性能的影响也较大【l “。当飞机在结冰条件下飞行时，飞机所有部件的迎风 面上都可能结上冰层。当冰在结冰表面形成后，必须采取一定的措施将其除去， 如热融化等。 图1 4 霜状冰图1 5 瘤状冰 飞机防冰系统的主要研究任务是采用物理模拟和数值模拟两类方法研究在 结冰融冰过程的数值模拟 结冰气象条件下飞机迎风面的结冰情况以及结冰后如何除冰。 1 3 认识发展过程及研究现状 人类从十九世纪就开始关注结冰、融冰问题了，并且做了大量的研究工作。 固液两相的热物性不同、液相区的自然或强迫对流、固相区的过冷是处理结冰 和融冰问题时需要面对的问题。而时间和空间维度上的非线性是数值方法处理 相变问题的难点所在。 1 9 8 1 年，s t e f a n 研究了具有移动边界和相变的冰的形成和融解过程。他推 导出潜热近似值，并以此来求解随时间变化的表面温度 1 7 。 二十世纪中期c a r s l a w 、j a e g e r 和c r a n k 等人主要研究了简单初始情况和边 界条件下无限空间和半无限空间下的一维带有相变的导热问题 1 8 。 1 9 8 6 年，h b e e r 和h r i e g e r 给出了圆柱体内固定融化的数值解【2 0 】。 1 9 9 6 年，c h i i ，d o n g 等对半无限大的冰层融化问题用积分法来求解偏微分方 程【2 1 】。 1 9 9 8 年孔祥谦提出“显热容法”的有限元方法用于相变导热分析，并对能 量方程中的相变潜热项提出了两种解决方案：作为附加比热和作为附加源项处 理【2 2 】。 国外早在上个世纪八十年代就开始了对多维结冰过程的数值模拟 2 3 【2 4 1 ，并 考虑了液相对流的影响，而且做了大量实验来验证数值模拟的正确性。而目前 国内的研究对对流换热的影响往往考虑的不够充分，文献 2 5 忡提到了有效导热 系数的概念，即将液相流动的影响通过有效导热系数加以修正，但也没有真正 涉及到液相中的流动问题。就所查阅到的文献来看，对结冰过程中液相区流动 的研究国内尚未见报道。 飞机热气( 电热) 防冰系统工作时，防冰腔内的热气( 电热层的热流) 将 机翼前缘加热结在前缘冰的底层融化，然后在气动力的作用下将冰吹除。文献 2 6 提出了一种改进焓方法，它在处理相变过程中令温度发生一个微小的变化， 从而使温度与焓之间建立起连续性关系。这样在公式中就可以用温度来代替焓 变量，然后用预估校正法进行求解。 南京航空航天大学硕士学位论文 1 4 本文的主要工作 本文采用数值模拟的方法对液相存在自然对流的结冰过程和飞机防冰系统 工作时的融冰过程进行了数值模拟。主要工作有： ( 1 ) 首先介绍了结冰的原理及特点，然后说明了研究结冰、融冰问题的意 义，对目前国内外的研究现状进行了综述。最后简要说明了本文的主要研究工 作。 ( 2 ) 介绍了固液相交传热，详细说明了冰水相变时的特点，并对数值求解 固液相变的数学模型进行了阐述，着重说明了焓法模型和显热容法模型的推导 及实旖过程。 ( 3 ) 用焓方法以机翼电除冰为例进行相变导热问题的计算，分析不同因素 对除冰过程的影响。以方腔内的水的结冰过程为例，进行了含自然对流的结冰 数值模拟，并把本文的计算结果与文献中的计算结果进行了对e e 。 结冰融冰过程的数值模拟 第二章相变传热及其数学模型 从传热学的角度来看，结冰、融冰问题的实质是伴有相变的传热问题，以 导热为主的相变物质的融化和凝固问题在传热学中称为s t e f a n 问题，即伴随有 潜热的移动界面导热问题。 早期对相变导热的研究，由于数学方法还不完备，而且缺少计算机作为必 要的计算工具，因此只能对少数简单情况进行分析求解。到了二十世纪后期， 随着计算机技术突飞猛进的发展，流动与传热问题的数值解法也有了飞速的进 步。人们对结冰、融冰的研究也更为深入。在长期的研究积累中，已经发展形 成了多个求解相变导热问题的较好模型，对结冰、融冰的数值模拟也更加准确， 并且己成功应用到诸多工程领域中。 2 1 固液相变传热 固一液相变问题的特点在于求解域中存在一个位置随时间变化的固一液界 面，对于单纯物质的相变，在确定的相变温度下，该移动界面是比较明晰的。 对于非纯物质，其融熔凝固现象发生在一个温度范围内，因此相变移动界面是 一个的模糊两相区。 s t e f a n 问题在忽略了自然对流和强制对流影响的时，主要特点是：控制方程 是非稳态导热方程，区域内存在着一个随时间移动的两相界面，在该界面上放 出或吸收潜热，因此，又称为“伴有相变潜热的移动界面问题”。但是实际上， 无论是纯物质还是混合物的相变，都会引起由于温度变化而在液相中产生的自 然对流。 在固液相变问题中，影响传热的因素较多：在相变过程中介质变化时，会 产生体积的变化，显著影响传热过程；在相变界面处产生吸热和放热现象，改 变传热系数，使换热过程呈不确定性；固液界面的移动会使问题呈非线性；液 相区的流场以及界面传热率之间的相互作用影响传热面的形状。这些都增加了 问题的复杂性，而实际上这些影响因素之间也是互相耦合的。 2 1 1 结冰、融冰相变问题的特点 与其它纯物质不同，水是一种比较特殊的相变物质，即水存在特有的4 。c 密 南京航空航天大学硕士学位论文 度反转，其密度在4 c 时达到最大值，故冰水相变系统中的自然对流不同于一般 的自然对流。水在相变温度l = 0 时融化，其浮升力b 为： b=g(p。一p一)(2-1) 式中：譬一重力加速度； p 。一为水温瓦时的密度； p 一一为水4 。c 时的最大密度。 浮升力b 的方向d 取决于水温瓦与对应于密度最大值之间的关系，即： d = ( k 一咒) ( t o l ) ( 2 - - 2 ) 如果水温疋4 。则d 0 ，因为p 。p 和b 0 ，这将引起向上的流 动；当0 d 8 ，浮升力b 及流动方向则向下，故0 d 1 2 是自然对流的一个分界条件，对应的温度可认为是流动方向向上流和向下流之 间的逆转温度27 1 。结冰时，在冰表面和周围水间的温度变化区域内存在两个自 然对流区。即紧贴冰层0 。c 表面的向上流和4 c 外的向下流。融冰时，在过程初 期，水膜较薄，传热主要依赖于热传导，随着融化区域的扩大，浮升力驱动的 自然对流开始起主导作用。 实验研究表面，冰水相变系统中自然对流方向的逆转温度大致在4 8 和 5 3 5 。c 之间，当疋处于5 5 c 干n5 6 。c 之间时，努赛尔数( n u ) 最小，当l 处于 4 7 和5 0 c 之间时，流动和传热是非常复杂的，一般只能借助实验的方法来解 决。 自然对流影响下的水的凝固和冰的融化问题，实际上就是相变过程受交界 面移动的冰层内的不稳定导热和冰表面与水间的对流耦合控制。 2 2 数值求解方法 相变传热问题的求解方法可分为分析法和数值法两类。由于这类问题在数 学上是一个强非线性问题( 即使其控制方程可能是线性的，但两相界面的能量 守恒条件是非线性的) ，相变界面的位置是事先无法知道的，而且它还随时间而 异，解的叠加原理不能使用，理论求解的困难即在于此，所以至今只能对很少 的简单情况才能获得分析解，一般情况下都只能用近似方法或数值模拟的方法 来求解。对多维情况，即使用用近似方法求解仍较困难，有时可行，最终结果 也是靠数值计算，所以数值方法是处理这类问题的极重要的手段。数值模拟的 结冰融冰过程的数值模拟 方法主要有：有限差分法，有限容积法，有限元法，有限分析法，边界单元法 等等。 1 )有限差分法：该方法是在每个节点上，将控制方程中每一个导数用相 应的差分表达式来代替，从而在每一个节点上形成了一个代数方程。 求解这些代数方程就获得了所需的数值解。 2 ) 有限容积法：该方法通过将守恒型的控制方程对控制容积做积分来导 出离散方程。在导出过程中，需要对界面上的被求函数本身及其一阶 导数的构成做出假定。用有限容积法导出的离散方程可以保证具有守 恒特性，且离散方程系数的物理意义明确，故目前应用最为广泛。 3 ) 有限元法：在有限元法中把计算区域划分成一系列元体( 在二维情况 下，元体多为三角形或四边形) ，在每个元体上取数个点作为节点， 然后通过对控制方程作积分来获得离散方程。有限元法是对传统变分 计算的归纳和总结，它既吸取了有限差分中离散处理的内核，又继承了 变分计算中选择试探函数并对区域积分的合理方法。 4 ) 有限分析法：用一系列网格线将区域离散，每一个节点与相邻的四个 网格( 二维) 问题组成计算单元，在计算单元中把控制方程中的非线 性项局部线形化，并对该单元上未知函数的变化型线做出假设，把所 选定的型线表达式中的系数和常数项用单元边界节点上未知的变量值 来表示，这样该单元内的被求问题就转化为第一类边界条件下的一个 定解问题，可以找出其分析解，得出该单元中点及边界上8 个邻点上 未知值之间的代数方程，此即为单元中点的离散方程。有限分析法中 的系数不像有限容积法中那样具有明确的物理意义，对不规则区域的 适应性也较差。 2 3 相变导热问题的数学模型 目前，在数值模拟相变导热过程中主要有固定网格技术和追踪界面的动网 格技术【2 。由于固定网格法在解决相变问题时，可以简单地从要求的温度场中 直接得出相变界面，因而被广泛的使用。在使用动网格技术时，由于温度场和 相变界面的位置都是变化的，因此，对相变位置的连续追踪就显得极为重要， 但是执行起来也是相当麻烦的。 使用动网格技术的主要问题是：处理当液相变为固相时，如何解决速度变 南京航空航天大学硕士学位论文 为零的情况。m o r g a n 3 0 运用了一个简单的方法，只要当某个计算网格的平均潜 热日值介于预定值o ( 整个网格均处于固相) 和相变潜热三( 整个网格均处于 液相) 之间时，即认为其速度为零。g a r t l i n g 【3 l 】用了一个更为巧妙的方法，采用 粘度来调整日的作用，当日从三降到0 时，粘度也随之增加了一个显著的量。 相应的，求解相变导热问题的数学模型，也可以主要分为两大类 3 2 】：第一 类：单区域模型( 即统一模型) ，因模型同时应用于固相区、液相区和两相区， 使得这些区域互相耦合，所以只需应用一套固定的网格和一组边界条件就可以 对问题进行求解了，因而统一模型被广泛采用。焓法、显热容法和不动界面法 ( 即坐标变换法) 是目前数值计算中较好的统一模型，且适用于多维问题的数 值求解。第二类：多区域模型，即对各相区分别建立控制方程组，通过在相界 面上建立合适的边界条件来描述固相区、液相区和两相区的耦合关系。 虽然统一模型不用跟踪相界面，但其难以修正界面上的潜热的释放，这影 响了数值解的收敛情况。为此，v o l l e r 和其同仁们提出了通过在迭代过程中恰当 地修正潜热来加快数值解收敛率的方法口引。b e n n o n 和i n c r o p e r a 针对求解统一方 程，就数值方法和差分格式展开了分析和讨论 3 ”。 最初研究相变导热问题一般都是以无限薄的相交界面把相区域d 分成两部 分( 如图2 1 ) ，即固相区域d 。和液相区d ，两部分有各自的温度场和物性参数， 如t ，p ，k ，o ，和乃，p l ，k ，a ，等，分别满足非稳态导热方程及各自的边界 和初始条件。当导热温度场越过相变区间时，物质都会吸收或放出大量的潜热 而且，为了简单，在液相区略去了自然对流或强制对流的作用。 t 0 图2 1 纯物质相变模型 早期的数值求解方法就是遵循上述分区的特点，只是在区域中划分了固定 结冰融冰过程的数值模拟 步长的网格节点，并利用固定的时间步长计算各时刻的温度场 2 ”。由于相变界 面是移动的，因此固定步长的缺点是对于每一个计算时刻，经常发生相变不落 在网格节点上，必须对界面附近的节点进行插值计算。它的改进方法是使用固 定的空间步长和可变的时间步长，使每经过一个时间步长都能保证两相界面正 好移动一个固定的空间步长的距离，这种时间步长都需随时迭代计算来确定。 另一种改进方法使采用自变量变换，将移动边界问题转化为固定边界问题 来处理。例如原来的求解域为： 0 x s ( t )( 移动相变界面) 作无因此变换，令：孝= x s ( t ) ，则求解域变成了0 善 1 ( 固 定相变界面) 为了更适用于多维问题的计算以及在整个区域( 包括固、液和两相区) 建 立一个统一的能量方程而发展了“焓法模型”。焓法模型采用焓和温度同时作为 待求函数。由于焓随时间的变化h f 曲线是连续的，因此用数值方法求解焓分 布时不需跟踪两相界面，从而使液相区和固相区统一处理成为可能。焓场解出 后，温度场就可容易解得。 把物质的相变潜热看作是在一定厚度的相变区域内有一个很大的显热容 量，这就是“显热容法模型”。随着相变过程的逐渐进行，相变潜热不断释出( 或 吸收) ，相变区内的物质温度也随之缓慢下降( 或上升) 。由于整体区域内温度 场是时间的连续函数，因此与焓法模型一样，不需跟踪两相界面，而在整体区 域内求解统一的温度场【3 5 】。 2 3 1 焓方法 2 | 3 1 1 焓方法数学模型 在推导整个区域的统一方程时，假定两相的密度相同而且是均匀的，两相 的其它物性值可以不等，也可以随温度变化。 对任一给定的控制容积q ，其表面积为r 、无内热源的能量方程为： 二i p e d ( ) = 1 2 v t n d f ( 2 - - 3 ) 西g# 式中：e 为单位质量的内能，或称为比内能；为表面的外法线向量。根据 内能e 和焓h 的关系式p e = p h p ，以及压力与时间无关，即： 昙p q 。o ( 2 _ 4 ) 1 4 南京航空航天大学硕士学位论文 所以，式( 2 - - 3 ) 可以改写为： - 未tn j p h d n 2 p a r ( 2 - - 5 ) 上述积分的形式的方程是以焓h 和温度丁一起作为待求函数的，它对整个求 解区域都适用。如果区域内存在相变，相变范围为阢，正】，其中，t = l 一， 互= l + ，只要给出焓与温度的具体关系式h = 仃) ，如 式中，h ；= c 。 t = ，+ 剖；， 0 一z ) c 。， i ，h ，= 。耳+ 。则相变的影响就自动地包含在方程( 2 - - 6 ) 中，因为焓h 中也包含了相变潜热三，式( 2 - - 5 ) 和式( 2 - - 6 ) 即是焓法模 型的基本方程。 为了证明焓法模型的基本方程和描述相变问题的常用方程之间的等价性， 现分别讨论控制容积中不包含两相界面和包含两相界面的两种情况。 1 控制容积中不包含两相界面 由于和九v r 在n 和r 上都是连续的，对式( 2 5 ) 的右端利用散度定理， 并考虑到q 与时间无关，得： 导( 户 p q = v ( x v r ) d n ( 2 7 a ) 即： l 鲁( , o h ) 山( 枷) 卜= o ( 2 _ 7 b ) 由于上式对在固相或液相区内的任一控制容积都适用，因此必然导致 导汹) 一v ( 兄v r ) = 0 ( 2 8 ) 利用d h = c d t 带入，故得： p c i t t l = v ( 胛丁) ( 2 9 ) t t g 这就是单相介质中常用的导热微分方程。 2 控制容积中包含两相界面 如图2 2 所示，在t 时刻，界面s 将控制容积分成固相区q 。和液相区q ，其 表面积分别为l 和l 。经过f 时间后，界面处于新的位置s ，扫过微元体积6 q 。 此时，和九v 丁在穿过移动界面时不再是连续的了，所以为了推导两相界面上 吖q 一 轧蚰 矾 结冰融冰过程的数值模拟 的守恒条件，需分别对固相区和液相区计算方程( 2 3 ) 的两端积分。 r s 6q 图2 2 推导两相界面守恒条件的控制容积 首先考虑将方程( 2 3 ) 左端的积分。在t 时刻： j p h d q = f ( p h ) ，d n + 渤x d q = f ( 肋xd r 2 + f ( 砌) ，d q + f ( 肋i d q ( 2 - - 1 0 ) 最n j 二劬 赢 在t + a t 时刻： f 肪d n = f l o h ) sd q + f 汹) ，d q 矗 o s ；舯o l 劬 = f ) ，d q + f 汹) ，d q + f 油) ，d r 2 ( 2 - - 1 1 ) 疽 扇 n ，二m 将式( 2 - - 1 1 ) 减去式( 2 一l o ) 后除以a t ，取f _ + 0 ，且 f 一6 q ) 斗q ， 得： 鲁肋a q2 导渤上a q + 击汹) ，a q + 1 i mf 堕立盟d q ( 2 - - 1 2 ) “齑 a t 注意到上式右端得最后一项，显然，当a t 斗0 时，d n 出叶v 。d s ( v 。为 两相界面微元面积d s 的局部法向速度，它指向液相区) ，并且6 q 收缩到界面s 。 同时，h ，和h f 分别达到其饱和焓 ：和砰，因此式( 2 - - 1 2 ) 化为： 击肋a q5 击沏xa n + 鲁) ，a q + p g ：一一j ) 1 一a s 2 1 3 其次，将上式中固相区q ，和液相区n ，的积分分别用式( 2 - - 5 ) 代入，并注 意到h ：一h ? = l ，上式化为： 击p 2 。f 。& v t n d f + 。f a v t n d f - p 一如 。1 4 南京航空航天大学硕士学位论文 上式右端的第一个积分中的n 表示q 。的外法线向量，第二个积分中的i 1 表 示q ，的外法线向量，它们的符号刚好相反。于是，由式( 2 - - 1 4 ) 得： 昙p q = p 州r 心警) ，一( a 署) ，钠卜( 2 - - 1 5 ) 此处r = r + f ，将上式减去( 2 - - 5 ) ，即有2 ( 九飘一( 九剽，讪。 q 叫s ， 这就是两相界面上的能量守恒条件。 通过上面两种隋况的讨论，焓法的正确性得到了证明。在采用这种方法时， 可以从如方程( 2 5 ) 的积分形式出发，也可以从与之相等价的微分方程 p 娑= v 伉v r ) 出发进行数值求解。 2 3 1 2 焓方法的应用 为了具体说明焓法的应用，现讨论如 下的二维结冰问题：假设相变是在一个非 常小的温度范围阪一e ，l + 】内发生的， 如图2 3 所示。在相变过程中，焓与温度 在此范围内是成线性关系的，根据温度的 不同，三个状态下总焓值( 包括显热热焓 和潜热热焓) 的表达式分别为： ( 1 ) 冰区： 即当t + 时，有： ( 2 1 7 ) ( 2 一1 8 ) 结冰融冰过程的数值模拟 - h = 岛。c ( r l ) + n ( c l + ) ( 2 1 9 ) 根据上一节推导的二维相变导热的焓方程( 2 1 6 ) ，即： 塑：旦( 马+ 旦( 马( 2 2 0 ) a f缸、缸。乱、却7 求解时根据相态确定所采用的焓的表达式，而相态是根据温度来确定的， 但是温度为待求解的变量，因此需要迭代求解，为了能过提高收敛速度，本文 采用了预估校正法。 所谓预估校正法，即先假定各区所处的状态，然后根据假设选择公式求解。 根据解得的温度确定相态，与假设的相态比较，若相同则开始下一个时