（动力工程及工程热物理专业论文）环形双层液体内热对流过程的二维数值模拟.pdf

重庆大学硕士学位论文 中文摘要 i 摘 要 在晶体生长过程中，由表面张力梯度引起的热毛细对流成为影响晶体材料品 质的重要因素。液封 czochralski 生长技术能有效地抑制熔体内的热毛细对流，但 是，目前对该生长技术中环形双层液体系统内热对流过程基本特征的了解非常缺 乏，因此，研究在水平温度梯度作用下环形双层液体内热对流过程的稳定性、发 生转变的物理机制及其基本特性等，对于改进生产过程、提高产品质量具有重要 的理论意义和实用价值。 本文建立了环形双层液体内存在水平温度梯度时热对流过程的物理模型和数 学模型，采用有限容积法进行了二维数值模拟，基于双层流体内的流场和温度场 分布，分析了marangoni (ma)数、浮力以及几何参数等对热对流的影响，得到了各 种不同条件下的临界ma数，探讨了流动失稳的物理机制。 结果表明： （1）当 ma 数较小时，流动较弱，为稳态的二维流动。随着 ma 数 和深宽比的增大，流动强度增强，等温线发生强烈的非线性变形；浮力的引入强 化了熔体层内的流动。 （2）当 ma 数超过临界值后，流动会失去稳定性，在冷壁附 近出现附加的小流胞，流动会转化为非稳定的多胞流动。小的流胞在冷壁附近形 成，并向热壁运动，最后在热壁附近消失。当 ma 数较小时，多胞结构主要集中在 冷壁附近，速度、温度振荡仅在冷壁附近较窄的区域内出现，随着 ma 数和深宽比 的增大，速度、温度振荡增强并向热壁方向扩展，多胞区域在更大的范围内出现。 （3）临界 ma 数随着深宽比的增大而减小，随着半径比的增大而增大。在有自由 表面的环形液池内， 浮力的引入会使临界 ma 数增大， 推迟流动的失稳， 液池越浅， 浮力的影响越明显。相同条件下，上部为固壁时临界 ma 数更小，流动更易失稳。 在上部为固壁的环形腔内，浮力的引入会使临界 ma 数减小，流动更易失稳。 关键词：关键词：数值模拟，双层流体，热对流，自由表面，液-液界面 重庆大学硕士学位论文 英文摘要 ii abstract in the crystal growth process, the thermocapillary flow induced by the interface tension gradient becomes a prominent factor influencing the quality of crystal materials. the liquid encapsulation czochralski (lec) technology can effectively suppress thermocapillary convection in the melt. but there is a lack of the understanding on the basic characteristics of thermocapillary convection in the annular two-layer liquid system. it is of great theoretical significance and practical value to investigate the stability, physical mechanism and characteristics of thermocapillary convection in the two-layer system under the horizontal temperature gradient so as to improve the production processes and quality of materials. the physical and mathematical models of thermal convection in the annular two-layer systems under the horizontal temperature gradient were established. a set of two-dimensional numerical simulations was carried out using the finite-volume method. the distributions of temperature and velocity in the two liquid layers were then obtained and effects of parameters, such as marangoni (ma) number, buoyant and geometric parameter et al, on the thermal convection were analyzed. also, the critical ma number at different conditions was determined and the physical mechanism of the unstable thermal convection was revealed. it is found that: (1) for small ma number, the flow is steady. with the increase of ma number and the aspect ratio, the flow strength becomes much stronger and the deformation of the isotherms increases sharply. in considering the effect of buoyancy, the thermal convection in the melt is enhanced. (2) while ma number exceeds the critical value, an unsteady multicultural structure is developed, and the additional cells appear near the cold lateral wall. during the oscillatory process, the flow cells move to the hot wall and then fade away in the vicinity of the hot wall. when the ma number is small, the multicultural structure and the oscillations of velocity and temperature mainly appear near the cold wall. with the increase of the ma number, the oscillations of velocity and temperature increase and propagate to the hot wall, and the multicultural structure expands to a wider region. (3) the critical ma number decreases with the increase of the aspect ratio and the decrease of the radius ratio. in the annular pool with a free surface, the buoyancy can delay convective instability of the two-layer system and make the critical ma number increase. the effect of buoyancy is much obviously in the shallow pool. the critical ma number decreases and the instability of the system 重庆大学硕士学位论文 英文摘要 iii enhances when the top surface is bounded by the rigid wall, and in this condition the buoyancy makes the system apt to lose stability and the critical ma number decrease. keywords: numerical simulation, two-layer liquid, thermal convection, free surface, interface 重庆大学硕士学位论文 主要符号表 vi 主要符号表 a：无因次流函数振幅 ：密度，kg/m3 g：重力加速度，m/s2 ：表面张力，n/m gr：grashof 数，gr=1g(th-tc)(ro-ri)3/12 ：无因次时间 h：液池深度，m ：流函数，m3/s ma：marangoni 数，ma= t1-2(th-tc)(ro-ri)/(1a1) ：无因次流函数 p：压力，pa p：无因次压力 角标角标 pr：prandtl 数，pr=/a1 1：下层流体 r：径向坐标，m 2：上层流体 r：无因次径向坐标 c：冷壁 t：温度， cri：临界值 t：时间，s h：热壁 u：径向速度，m/s max：最大值 u：无因次径向速度 p：周期 v：轴向速度，m/s v：无因次轴向速度 z：轴向坐标，m z：无因次轴向坐标 希腊字母希腊字母 a：热扩散率，m2/s ：体积膨胀系数，k-1 t：表面张力温度系数，n/(mk) t1-2：界面张力温度系数，n/(mk) ：液池半径比，= ri/ro ：液池深宽比，=h/(ri- ro) ：无因次温度 ：导热系数，w/(mk) ：动力粘性系数，kg/ms ：运动粘性系数，m2/s 重庆大学硕士学位论文 1 绪 论 1 1 绪 论 1.1 引言 晶体材料在工业领域和日常生活中具有重要作用，已成为不可缺少的重要原 料。对高质量晶体的需求日益增加，使得对晶体材料生长过程中流动的物理机制 的研究显得尤为必要。 czochralski（cz）方法是熔体中生长晶体材料的一种最主要的方法，在这种方 法中，在坩埚和晶体间总存在着熔体自由表面，由于温差的存在，在自由表面上 形成表面张力梯度，拽引界面并通过液体的粘性而诱发主流体运动，形成热毛细 对流。在常重力条件下，温差还会引起密度的变化，在熔体内会形成浮力驱动的 自然对流。当坩埚和晶体间温差较小时，流动是稳定的轴对称运动，随着温差的 增大，流动将失去其稳定性转化为振荡对流，不论是由于浮力驱动还是表面张力 驱动，一旦出现了振荡对流，熔体中的流场、温度场以及浓度场都会随时间变化。 温度的不均匀性还会引起晶体内部应力的不均匀，使晶体内出现条纹及各种缺陷， 从而大大影响晶体生长质量。 为了抑制自由界面表面张力梯度引起的热毛细对流，并防止晶体生长过程中 可挥发性成分的挥发、改善结晶过程中的热传递条件，近年来发展起来的液封 czochralski 生长技术1在这方面表现出极大的前景。所谓液封技术即在熔体自由 表面上覆盖一层与熔体不相混溶的、不起化学反应的流体，这样通过自由表面和 液-液界面的相互抑制，在一定的条件下可以使熔体主流区的热毛细对流得到削 弱。 对于液封cz生长技术中的双层液体系统而言，由于两液层界面处的热耦合和 力耦合的复杂性，使得双层流体系统内的热对流过程变得非常复杂，给理论和实 验研究带来了很大的困难。目前，对矩形腔内双层流体内的热毛细对流的研究已 有不少，但对环形双层液体系统内热对流过程的基本特征的认识却非常缺乏。研 究这些对流现象对控制材料生长、得到新兴材料具有重要的指导意义。本课题采 取数值模拟方法研究环形双层液体内存在水平温度梯度时的热对流过程，分析上 层流体对下层流体的抑制效果，确定流动失稳的临界条件，揭示双层流体内流动 过程的基本规律以及流动失稳的物理本质。本课题在理论上可丰富和发展双层流 体内热对流及其稳定性理论，在实践上可指导和改进液封czochralski晶体生长过 程，因此，具有重要的理论意义和实用价值。 重庆大学硕士学位论文 1 绪 论 2 1.2 研究现状 1.2.1 双层流体热对流的实验研究 对双层流体系统热毛细对流的实验研究最早始于 1988 年，villers 和 platten2 测量了矩形腔内两层流体（水和庚醇）热毛细对流时的速度分布，并建立了一个 简单的理论模型对速度分布进行计算。1991 年，prakash 等3通过实验研究了水平 温度梯度下的两不相混合流体系统。在浮力和表面张力的作用下，双层流体内流 动以浮力对流为主，表面张力驱动的对流很弱。1999 年我国首次利用实践五号科 学实验卫星成功完成了两层流体空间实验4，在两种微重力水平下，获得了垂直于 流体界面加热的 marangoni 对流与平行于流体界面加热的热毛细对流两种模型的 流动图形和温度数据。刘秋生等5-7又通过理论分析、数值模拟和空间实验研究了 微重力环境中两种不同加热方式下的双层流体实验系统的对流特征，实验清楚的 观察到了典型的定常 marangoni 对流和热毛细对流现象， 将其实测结果和数值模拟 分析结果进行了比较，二者的对流流动结构一致、速度大小基本吻合，验证了微 重力双层热毛细对流的理论模型和流动特征。2002 年，someya8等用由硅油和氟 液组成两流体层，分别对有无自由表面两种情况进行了实验，他们观察到了水平 温度梯度驱动的腔内流动，并用 piv 技术测量了交界面处的对流流动，实验测得 的流场和数值研究的结果一致。 近年来，人们从实验和理论两方面研究了两层流体系统热毛细对流中的振荡 现象，张嘉锋9等利用实验和数值模拟的方法研究了由 10 号硅油与 fc-70 组成的 双层流体系统，获得了不同厚度条件下两层流体对流流场中心截面的二维速度场 分布以及液-液界面处的二维速度场、温度场分布和涡胞结构，确定了不同厚度条 件下的临界温差，发现上层流体越薄，临界温差越大，实验得到的临界温差与数 值模拟结果基本吻合。2003 年，simanovskii 等10进行空间实验研究了三层系统内 marangonibenard 对流的不稳定性，系统由热扩散率相近的流体组成，首次用实 验方法观察到了微重力条件下的纯热毛细对流现象，结果证实了理论上预测的振 荡不稳定性的存在，并利用有限差分法研究了振荡对流的非线性特点。洪勇等11 实验观察了高温 bi12sio20熔体中表面张力对流从稳态向振荡态的转变过程，发现 稳态热毛细对流的模式由一个主体对流和两个分支组成，而主体对流总是指向坩 埚的低温点；当振荡态出现后，主体对流区域的长度会随着温度的增加而增长， 振荡频率也随之增加。 1.2.2 双（多）层系统的理论研究和数值模拟 目前，已有不少学者对矩形腔内双层流体内的热毛细对流进行了研究。1990 年， villers12等分别在考虑热毛细力和不考虑热毛细力的条件下，通过理论分析研 究了侧壁加热的矩形腔内两不相混流体层内的热对流，获得了各层流体内的速度 重庆大学硕士学位论文 1 绪 论 3 分布，与矩形腔实验测得的速度分布基本一致。1993 年，doi 和 koster13 研究了 上部有自由表面时两不相混溶流体的稳态热毛细对流特征。在无限大水平流体层 内得到了速度、温度分布的分析解，发现随界面张力温度系数与表面张力温度系 数之比 的改变， 存在四种不同的流型， 当 =0.5 时， 对下层流体的抑制效果最佳。 为了研究垂直边壁对双层流体对流的影响，进行了二维数值模拟，结果表明，由 于边壁效应的影响，只有在较低 marangoni（ma）数或毛细雷诺数下，=0.5 时的 抑制效果最好，随着 ma 数或毛细雷诺数的增加，最佳抑制效果时的 值会减小， 上述结论与实验结果基本一致。刘秋生等14给出了一组矩形腔内双层流体热毛细 对流的数值计算结果，模拟了水平加热封闭腔和上部为自由表面时两不混溶流体 内的热毛细对流（微重力），研究了不同流体粘性比和热扩散率比的影响，得到 了水平无限大两层流体速度分布的渐近解15。1994 年，wang 等16对有自由表面 的两流体层的浮力-marangoni 对流进行了数值模拟。在纵横比较大时，得到了流 体流动的四种不同的流型。人们对热毛细对流进行深入研究后，发现在很多液-液 系统中存在非正规的热毛细效应，即界面张力会随温度的升高而增大。braverman 等17研究发现在浮力和这种异常的热毛细力共同作用下，产生了一种特殊形式的 振荡对流。2002 年，boeck 等18通过直接的三维非线性数值模拟，研究了在特定 的液-液双层系统中，浮力和异常的热毛细力共同作用下的对流流动，发现了稳定 状态到振荡不稳定状态的转变， 当 ma 较小时， 流动总是稳定的， 当 ma 足够大时， 随着 grashof（gr）数的增加，稳定流动转变成振荡流动。gr 数较小时，流动为 稳态流动，随着 ma 的增大，流动从稳定的六角形流胞转变为交替的滚胞。 近年来，对温度梯度垂直于界面的情况已经进行了比较多的研究。2002 年， 刘秋生19给出了无限长矩形腔内两层不相混流体的 marangoni 对流的线性稳定性 分析结果， 比较了两层流体的 marangoni 对流与两层流体外加平行于流体界面的温 度梯度所引起的热毛细对流的主要特征，在相同温差条件下，两层流体内的 marangoni 对流通常比热毛细对流弱，且两者间的强度差随 ma 数的增大而增加。 2002 年，nepomnyashchy 等20采用线性稳定性分析和非线性数值模拟研究了 两不相混合流体层系统内的流动，两流体系统上部加热，液-液界面上有热源。分 析了热毛细对流不稳定机制，界面上的温度不均匀引起了热毛细对流的不稳定性。 2003 年，nepomnyashchy 等21又对一个真实的两层流体系统中的振荡流动进行了 线性和非线性模拟。发现浮力使振荡减弱并最终完全抑制振荡，仅在微重力条件 下能观察到振荡的不稳定性，当液层厚度比一定时，振荡不稳定性的临界 ma 数存 在一个最小值。2004 年，nepomnyashchy 等22分别对浮力和热毛细力作用下和界 面有热源时的振荡对流进行了线性和非线性模拟，gr 数较小时，周期性边界条件 和刚性固壁边界条件下都观察到了振荡流动，随着 gr 数的增大，振荡消失。 重庆大学硕士学位论文 1 绪 论 4 在垂直于界面的温度梯度作用下，双层流体系统失去稳定后会产生 rayleigh- marangoni-benard（r-m-b）或 marangoni-benard 对流。tavener 等23阐述了在考 虑界面变形时，垂直界面加热两不相混合流体的 r-m-b 对流。在有限的深宽比范 围内，用数值方法计算了表面张力和浮力驱动的不稳定性，结果发现，临界 ma 数随上、下层体积比增大而减小，随导热系数比增大而增大，随 rayleigh（ra） 数增大而增大。2004 年，周炳红等24对 10 号硅油和 fc70 组成的两层流体系统的 r-m-b 不稳定性进行了数值研究，发现两流体系统的不稳定性主要取决于两层流 体的深度比。 2005 年， 刘秋生等25采用线性稳定分析和二维数值模拟研究了 10cst 硅油和 fc70、2cst 硅油和水等两层系统内的 r-m-b 对流不稳定性，在 10cst 硅油 和 fc70 系统中，随着 ra 数和 ma 之比的减小，振荡增强；对 2cst 硅油和水系统， 随着 ra 数的增大，稳定的对流分叉为振荡对流，他们还提出了二次振荡不稳定性 并成功地解释了实验观察和线性稳定性分析的分歧。2006 年，刘秋生等26对硅油 和水组成的真实的两层系统进行了数值模拟，结果发现，深度比较小时，临界 ra 数随着深度比的增大而增大， 深度比较大时， 临界 ra 数随着深度比的增大而减小。 mcfadden 等27对水平两层流体系统在考虑相变情况下， 进行了线性稳定性分 析，在垂直温度梯度时考虑浮力和热毛细力的作用，结果发现不同相之间熵的差 值是决定系统稳定性的至关重要的因素，当该差值较小时，无论从顶部或底部加 热都会出现不稳定；当差值较大时，只有在底部加热才会出现不稳定。 目前，关于双层液体内水平温度梯度作用下的热毛细对流和浮力-热毛细对流 的研究还比较少。2003 年，madruga 等28采用线性稳定性方法对水平加热无自由 表面的两流体层的热毛细-浮力对流进行了研究。当温度梯度较小时，流动为稳定 的二维流动，随着界面张力系数不同，可能出现四种流型。当温度梯度达到临界 值后，流动失稳变为三维非稳态振荡流动。在不同的相对液层厚度下，可能出现 三种流型：从冷区向热区运动的热波、从热区向冷区运动的热波或稳定的轴向滚 胞。2004 年，madruga 等29又对水平温度梯度作用下氧化硼和砷化镓流体系统的 稳定性进行了研究，结果发现，在热毛细力和浮力的共同作用下会产生稳定的流 动，为此，得到了系统稳定状态下的四种不同的流型和五种温度曲线分布。流动 失稳会引起一个所谓的热流体波，这些热流体波是平行、还是垂直于温度梯度， 则主要由两层流体的相对高度决定。 2006 年，nepomnyashchy 和 simanovskii30研究了两固定平板间具有水平温度 梯度的双层流体的非线性稳定性，考虑了两种类型的边界条件：周期性边界条件 （无限大液层）和侧壁绝热边界条件（封闭腔）。对 5cst 硅油/ht-70 系统进行了 非线性模拟，发现波的传播方向取决于两个因素，液层厚度比和 marangoni 数。 simanovskii 等31研究了矩形腔内的两层流体系统的振荡对流，考虑浮力和热毛细 重庆大学硕士学位论文 1 绪 论 5 力的共同作用，当 gr 数较小时，随着 k(bo 数的倒数)的增大，流动失去稳定转变 为振荡对流，当 gr 数足够大时，流动为稳态流动。 2006 年，gupta等32-33在不考虑重力效应时，对温度梯度平行于液-液界面的 两不相混溶流体的热毛细对流进行了研究，封闭矩形腔的研究结果表明：引入粘 性更大的液封流体会极大地削弱液封内的热毛细对流，液封的引入导致了界面变 形，液封层的流型和界面变形都与液封的厚度和粘性有关；开口矩形腔的结果表 明：上部自由界面削弱了熔体层的热毛细对流，增强了液封层的对流，选择界面 张力对温度更敏感的液封流体，几乎可以完全抑制熔体内的热毛细对流。随着液 封层粘性的增大和液封厚度的减小，液封层内的热毛细对流强度会减弱。 对存在外加磁场情况下热毛细对流的研究也有不少报道。2001 年，rao 等34 对考虑浮力时、均匀磁场对侧壁加热的浅腔内 benard-marangoni 对流的影响进行 了研究。在近似忽略边壁效应的薄液层内，得到了流函数和温度的封闭形式的解， 分析发现，对不同的 ，存在四种不同的流型，通过恰当的选择磁场，可以控制下 层流体的对流。2007 年，ludovisi 等35给出了腔内的两层流体在不同的磁场强度 下的数值计算结果，发现流动除了受到浮力和热毛细力的影响外，还受到磁场的 影响。非均匀磁场的引入可以强化或抵消浮力作用，通过变化磁场强度和梯度可 以改变系统的速度和温度分布，进而控制传热。黄护林等人36采用数值模拟方法 研究了磁场对水平温度梯度作用下的双层流体热对流的影响，结果表明，常重力 条件下，磁场能促进热剪切层内的稳定。2009 年，cha 等人37通过数值模拟研究 了磁场对腔体侧壁加热的热毛细对流的影响，结果发现，磁场对速度的影响可以 是强化、也可以是削弱。 刘秋生等人38对高频振动下的两层流体系统的 r-m-b 对流的热振荡不稳定性 进行了线性稳定性分析，结果表明，高频热振动可以改变两流体系统对流的不稳 定性，水平的高频振动强化对流流动，使系统更易失稳；垂直的高频振动抑制对 流流动，推迟系统的对流不稳定性的发生；与界面成六十度角的振动对系统不稳 定性的影响很小。2008 年，kovskaya39分析了垂直振荡对界面变形的两流体系统 的影响，高频振荡可以使界面变平，并可以抑制不稳定性。2009 年，kovskaya40 进一步对两不相混流体系统中高频渐进振动对热对流的影响进行了研究，发现， 振动的水平分量强化对流流动，导致了系统的失稳；垂直分量抑制对流流动，使 两层流动趋于稳定。 近年来，人们开始研究具有二个流体交界面的多层(三层)流体系统。georis 等41,42从理论上研究了底部加热的对称的三层流体系统的 marangoni-benard 不稳 定性，线性稳定性分析和非线性数值模拟表明，热扩散率比决定了不稳定性的性 质，当两层流体热扩散率相差较大时，流动为稳态流动；当热扩散率相近时，观 重庆大学硕士学位论文 1 绪 论 6 察到了振荡对流。刘秋生等43介绍了自然界和工程技术中的多层不混溶流体系统 及其流动现象，论述了多层流体质量传输和热传递规律的理论和实验研究成果及 进展，利用理论分析、实验研究和数值模拟方法研究两层及三层流体内的自然对 流及热毛细对流，并分析探讨覆盖液体层对被覆盖液体的动力控制特性及其系统 的稳定性。 nepomnyashchy 和 simanovskii 等 44,45分别对两层和三层流体系统中的 marangoni 对流进行了非线性模拟，研究了不同振荡流型间的转变，结果表明，开 始振荡是非常简单的，并且几乎呈正弦形式，流线及温度场满足对称条件，随着 ma 数的增大，对称的振荡流动变得不稳定并且发展为非对称流动。2005 年， shevtsova 等46采用有限差分法模拟了对称的三层系统中的 marangoni-benard 对 流，研究了对流流动的稳定性及不同流型间的转变。随着 ma 数的增大，稳定的流 动发展成非稳定的振荡流动，并且在长腔体中出现了沿温度梯度方向运动的热流 体波。 2007 年，simanovskii 等47,48分别在有无浮力情况下，对水平温度梯度下的三 流体层系统内的热对流进行了数值模拟。在无限大液层中，不考虑浮力且 ma 数较 大时，波在不同液层内的流动方向不同，其余情况下热流体波都是向热壁方向运 动的。simanovskii 等49在封闭腔内模拟了水平温度梯度作用下的三层流体系统从 单胞结构到多胞结构的转变过程。2010 年，他们又在封闭腔内对温度梯度垂直于 界面的三层流体内浮力-热毛细对流进行了数值模拟50，观察到了振荡流动。 1.3 本课题的研究内容 综上所述，目前，对水平温度梯度作用下双层液体内的热毛细对流或热毛细- 浮力对流的基本特征、失稳的临界条件、失稳的物理机制及失稳后的流动结构等 都还不是很清楚，而且研究主要集中在无限大的双液体层或矩形腔内的双液体层。 而在晶体液封 czochralski 生长技术中涉及到的是环形双层液体内存在水平温度梯 度时的热对流过程，这方面的研究仍然是一片空白。本课题首次采用数值模拟方 法研究晶体液封 czochralski 生长技术中涉及到的环形双层液体内存在水平温度梯 度时的热对流过程及稳定性。 本课题的主要研究内容如下： （1）建立常重力及微重力条件下，环形双层液体内存在水平温度梯度时热对 流过程的物理模型和数学模型。采用有限容积法对水平温度梯度作用下双层液体 内的热毛细对流或热毛细-浮力对流进行数值模拟。 （2）获取微重力及常重力条件下稳态流动时速度场和温度场的分布规律，分 析环形双层液体的几何参数及 ma 数等对流动的影响。 重庆大学硕士学位论文 1 绪 论 7 （3）获取失稳后的各种可能的流动结构，确定流动转变的临界条件，寻找达 到最佳抑制效果的各种参数。 （4）揭示双层液体内流动失稳的物理机制，指导晶体材料液层 czochralski 生长过程，提高晶体质量。 重庆大学硕士学位论文 2 物理数学模型及数值求解 8 2 物理数学模型及数值求解 热毛细对流是由热的原因引起的一种流体运动，它不仅在微重力环境下存在， 即使在地球重力场中其产生的效应依然不容忽视，本课题主要研究常重力及微重 力条件下，环形双层液体内存在水平温度梯度时热对流过程及其稳定性。本章的 目的在于通过一定的简化和相关假设来建立物理和数学模型。 2.1 物理模型和相关假设 物理模型如图 2.1 所示， 环形液池内、 外壁分别维持恒定温度 tc和 th (tcth)， 内半径为 ri，外半径为 ro，熔体（下层）高度为 h1，液封（上层）高度为 h2，底 部为固壁。分别考虑了上部为固壁或自由表面两种情况。 在模型中引入如下假设： (1) 熔体和液封均为不可压缩的牛顿流体，满足 boussinesq 近似； (2) 流速较低，流动为轴对称二维层流； (3) 液-液界面和自由表面平整无变形， 在液-液界面和自由表面考虑热毛细力的作 用，在固-液界面满足无滑移条件； (4) 液池顶部和底部边界均绝热； (5) 表面张力是温度的线性函数。 图 2.1 物理模型 fig 2.1 physical model 为简化起见，取 z 轴右侧为研究对象。定义环形池的深宽比为=h1/(rori)，半 径比为=ri/ro。界面张力随温度呈线性变化，即：() ctc tt=，其中， c 是 温度为 c t 时的表面张力， t 为流体的表面张力温度系数，/ t t= , t 可以大 于 0，也可以小于 0。 2.2 数学模型 由上述物理模型及相关假设，在二维坐标系下可以导出两层流体的控制方程 及边界条件。 重庆大学硕士学位论文 2 物理数学模型及数值求解 9 2.2.1 控制方程及边界条件 控制方程 所研究的区域包括熔体和液封流体，研究区域中流体的控制方程包括连续性 方程、动量方程和能量方程。 1) 连续性方程： ()1 0 i i ruv zrr += (2.1) 2) 动量方程： r 方向的动量方程： 22 222 11 iiiiiiii iii i uuupuuuu vu tzrrrrrzr += + (2.2) z 方向的动量方程： () 22 22 11 iiiiiii iiiici i vvvpvvv vugtt tzrzrrrz += + (2.3) 3) 能量方程： 22 22 1 iiiiii iii tttttt vua tzrrrrz +=+ (2.4) 其中，熔体和液封流体分别用下标 i1，2 表示。u，v 分别表示 r，z 方向的速度 分量，p 表示压力，t 为时间，和 a 分别为密度、运动粘度、热膨胀系数和 热扩散率。 边界条件 1) 内壁面，r = ri , 0 z h 1122 0uvuv=, 12c ttt= (2.5a-b) 2) 外壁面，r = ro , 0 z h 1122 0uvuv=, 12h ttt= (2.6a-b) 3) 熔体底部，ri r ro ,z=0 11 0uv=, 1 0 t z = (2.7a-b) 4) 液封顶部，ri r ro ,z=h 自由表面： 2 0v =, 22 2t ut zr = , 2 0 t z = (2.8a-c) 固壁： 22 0uv=, 2 0 t z = (2.9a-b) 5) 液-液界面，ri r ro ,z=h1 12 0vv=， 12 uu=， 121 121 2t uut zzr = ， (2.10a-c) 重庆大学硕士学位论文 2 物理数学模型及数值求解 10 t1=t2， 12 12 tt zz = (2.10d-e) 其中， 和 分别是动力粘度和导热系数， t 和 1 2t 分别为流体的表面和界面张力 温度系数。 2.2.2 控制方程及边界条件的无量纲化 在传热学研究中，为了使某一特定问题的结果适用于满足相似条件的同类问 题，进而将结果推广到一般情况，需要采用无量纲形式的方程。无量纲的采用还 避免了在求解有量纲方程时因单位转换而可能发生的错误，对长度、时间、速度、 温度和压力做无量纲变换如下： 温度： c hc tt tt = ；速度： 1 / oi uu rr = ， 1 / oi vv rr = ；时间： () 2 1 / oi rr t =； 压力： () 2 11 2 / oi pp rr = ；长度： oi r r rr = ， oi z z rr = 。 控制方程的无量纲化 1) 连续性方程： ()1 0 i i ruv zrr += (2.11) 2) 动量方程： 22 1 222 1 1 iiiiiiiii ii i uuupuuuu uv rzrrrrzr += + (2.12) 22 1 22 11 1 iiiiiiiii ii i vvupvvv vugr zrzrrrz += + (2.13) 3) 能量方程： 22 22 1 11 iiiiiii i a vu zrpr arrrz +=+ (2.14) 其中，r和z是无因次坐标，u、v为无因次速度分量，p为无因次压力，为无因 次时间，gr = g1(th- tc) (ro- ri) 3/12为格拉晓夫数，pr=/a1为普朗特数。 边界条件的无量纲化 1) 内壁面，r= ri/(ro- ri), 0zh/(ro- ri) 1122 0uvuv=, 12 0= (2.15a-b) 2) 外壁面，r= ro/(ro- ri), 0zh/(ro- ri) 1122 0uvuv=, 12 1= (2.16a-b) 3) 熔体底部，ri/(ro- ri) r ro/(ro- ri), z=0 11 0uv=, 1 0 z = (2.17a-b) 4) 液封顶部，ri/(ro- ri) r ro/(ro- ri), z=h/(ro-ri) 重庆大学硕士学位论文 2 物理数学模型及数值求解 11 自由表面： 2 0v =, 22 11 2 t t uma zprr = , 2 0 z = (2.18a-c) 固壁： 22 0uv=, 2 0 z = (2.19a-b) 5) 液-液界面，ri/(ro- ri) r ro/(ro- ri), z =h1/(ro-ri) 12 0vv=, 12 uu=, 1221 1 uuma zzprr = , (2.20a-c) 12 =, 12 12 zz = (2.20d-e) 其中，marangoni数定义为ma= t1-2 (th-tc) (ro- ri) /(1a1)。 流函数 ( 3 :/ms) 11 ,uv rzrr = = (2.21) 流函数满足连续性方程，的控制方程： () 22 22 u vrr rzrz += (2.22) 取 1( ) coi rr= 无量纲化得： () 22 22 u vrr rzrz += (2.23) 2.3 数值求解 数值解法是一种离散近似的方法，能得到被研究区域中某些代表性点上物理 量的近似值。对于前面已经建立的环形双层液体内热毛细对流的物理数学模型， 控制方程是复杂的非线性偏微分方程，边界条件也比较复杂，很难获得其精确解， 因此，拟采用数值方法来获取热对流过程的流场和温度场。 本课题采用有限容积法来对控制方程及边界条件进行离散，对流项用二阶迎 风格式，扩散项采用中心差分格式，压力-速度修正采用simple方法。对控制方 程及边界条件的具体离散过程参见文献51。 2.4 计算条件 为了研究影响环形双层液体内热对流的主要因素，本课题分别对有自由表面 的环形池和上部为固壁的环形腔内的流动进行了模拟，模拟流体为5cst硅油 /ht-70及砷化镓/氧化硼流体，物性参数如表2.1所示。 有自由表面的环形池： 5cst 硅油/ht-70双层流体：(a)微重力条件， 半径比=0.2， 深宽比=0.05、0.075、 0.1，(b)常重力条件，半径比=0.2，深宽比=0.05、0.075、0.1； 砷化镓/氧化硼双层流体： 常重力及微重力条件， 半径比=0.5， 深宽比=0.25； 重庆大学硕士学位论文 2 物理数学模型及数值求解 12 上部为固壁的环形腔： 5cst 硅油/ht-70双层流体：(a)微重力条件， 半径比=0.2， 深宽比=0.05、0.075、 0.1，(b)常重力条件，半径比=0.2，深宽比=0.05、0.075、0.1及半径比=0.6，深 宽比=0.1。 在微重力条件下，取gr=0，在常重力条件下，取g=9.86m/s2。 表 2.1 流体的物性参数 table 2.1 physical properties of the fluid 流体种类 ht-70(1) silicon 5cst (2)gaas (1) b2o3 (2) 密度 (kg/m3) 1680 920 5720 1648 热膨胀系数 (k-1) 1.110-3 1.0510-3 1.8710-4 9.010-5 粘性系数 (kg/(ms) 8.410-4 4.610-3 2.7910-3 3.9 热容量 a (j/(kgk) 962 1590 434 1830 导热系数 (w/(mk) 0.07 0.117 17.8 2.0 表面张力系数 t (n/(mk) - 7.210-5 - -3.5710-5 界面张力系数t1-2 (n/(mk) 7.310-5 1.210-3 2.5 程序正确性及网格收敛性验证 2.5.1 上部为自由表面 为了检验程序的正确性和网格的收敛性， 在与文献13相同的条件下取不同的 非均匀网格进行了模拟计算。计算时，取半径比为 =0.99，这样环形液池可近似 看成矩形液池，表中速度和温度均为自由表面中心点处的值，结果如表2.2和图 2.2所示。显然，当网格为121r61z时，自由表面中心点处速度和温度以及最大流 函数的值变化都很小，同时计算结果与文献13的结果非常接近，说明程序是正确 的，网格是收敛的。 2.5.2 上部为固壁 在与文献30完全相同的条件下，当左侧为高温、深宽比=0.0625、并取半径 比 =0.999时，可将环形池近似成矩形池，模拟结果如图2.3所示。在ma超过临 界值的条件下，稳态流动转变成非稳定的振荡流动，流胞向热壁方向流动。随着 流动的进行在冷壁的附近出现的附加流胞与中心处的脉动大流胞结合。 与文献30 对比，在流型结构上基本一致，因此，可以说结果是可靠的，程序是正确的。 重庆大学硕士学位论文 2 物理数学模型及数值求解 13 表 2.2 程序正确性和网格收敛性检验 table 2.2 check for the code validity and the mesh convergence (pr=0.01, ma=8102, 2/1=0.1, 2/1=10, 2/1=1, t2/t1=2.) 网格 项目 结果 文献13 相对误差，% u 1.59 1.666 4.5 t 0.610 0.585 4.2 41r21z (+) 24.592 - - u 1.552 1.612 3.7 t 0.6099 0.602 1.3 81r41z (+) 22.729 - - u 1.5629 1.600 2.3 t 0.6106 0.606 0.76 121r61z (+) 22.764 - - 说明：为便于比较，表中无因次速度的定义与文献13相同 (a) 网格为 41r21z, (+)=4.12, ()=-24.59. (b) 网格为 81r41z, (+)=4.13, ()=-22.73. (c) 网格为 121r61z, (+)=4.24, ()=-22.76. 图