（水力学及河流动力学专业论文）基于fluent的冰塞数值模拟.pdf

摘要 冰凌是冬季寒冷地区河流上普遍存在的一种水文现象。这些地区河流 冬季常形成冰盖、冰塞或冰坝，它们会带来各种冰凌灾害。现在随着计算 机技术的飞速发展，人们在研究冰水两相流问题时，越来越多的采用了数 学模型的方法。数学模型以其方便高效，费用低廉的优势，已经成为研究 冰水两相流运动规律的重要手段之一。依据传热理论、数值传热的方法， 使用f l u e n t 软件中凝固与融化解算器，对河流冰塞的形成及消融进行模 拟，研究了不同流速、温度、风速条件下的断面温度、冰盖厚度等流场参 数的变化情况。模拟结果显示：减小入口流速，降低空气温度，增大风速 等，形成冰盖的厚度都会增大；增大入口流速，增大空气温度，减小风速 等，会使冰盖的融化加快。并将实验室u 型弯道冰盖下冰颗粒堆积的实验 所得冰盖厚度与数值模拟所得冰盖厚度进行比较，比较结果显示：u 型弯 道时，0 。断面上冰盖厚度最大，随着断面角度的增大，弯道内凸岸的冰 盖厚度逐渐增大，而凹岸的冰盖厚度逐渐减小。 关键词：冰塞；数值模拟；f l u e n t ；冰塞形成；消融 a bs t r a c t i c ei sac o m m o np h e n o m e n o ni nt h er i v e r so fc o l dr e g i o n s r i v e r s m t h e s er e g i o n so f t e nf o r mt h ei c ec o v e r ，t h ei c eja mo rt h ei c ed a m i nw m t e r ， a n d t h e y w i l lb r i n gv a r i o u s i c ed i s a s t e r s n o w a d a y sw i t h t h e r a p i d d e v e l o p m e n to ft h ec o m p u t e rt e c h n o l o g y ，t h em a t h e m a t i cm o d e l 1 sa d o p t e d m o r ea n dm o r ew h e ni tc o m e st os t u d yi n i c e 。w a t e rt w o 。p h a s ef l o w lh e m a t h e m a r l c a lm o d e lh a sb e c o m ea ni m p o r t a n tm e t h o dw h e ns t u d y i n g1 n t h e i c e - w a t e rt w o p h a s ef l o wd u et o i t sc h a r a c t e r i s t i c so fc o n v e n i e n c e ，h l g n e f f e c t i v e n e s sa n dl o wc o s t a c c o r d i n gt h eh e a t t r a n s m i s s i o nt h e o r ya n dt h e m e t h o do fn u m e r i c a lh e a tt r a n s f e r ，b yu s i n gt h e s o l i d i f i c a t i o n m e l t i n g c a l c u l a t i o nm o d e li nf l u e n t ，s i m u l a t i o n so fr i v e ri c ej a mf o r m a t i o na n d m e l t i n ga r er e s e a r c h e d o np a r a m e t e rc h a n g e so fd i f f e r e n tt e m p e r a t u r em s e c t i o na n di c et h i c k n e s su n d e r i n gc o n d i t i o n so fd if f e r e n t f l o wv e l o c i t y 、 t e m d e r a t u r ea n dw i n ds p e e d t h er e s u l t so fc a l c u l a t i o nr e v e a lt h a t ：t h i c k n e s s o fi c es h e e tf o r m a t i o nc a ni n c r e a s eb yr e d u c i n g e n t r a n c ef l o wv e l o c i t y ， r e d u c i n ga i rt e m p e r a t u r ea n di n c r e a s i n gw i n ds p e e d ，e t c i n c r e a s i n g e n t r a n c e f l o wv e l o c i t y ，i n c r e a s i n ga i rt e m p e r a t u r ea n dr e d u c i n gw i n ds p e e d ，e t c w 1 儿 s p e e du pt h em e l t i n gi c es h e e t t h i sp a p e rm a k e sc o m p a r a t i v ei n v e s t i g a t l o n o n i c et h i c k n e s sb e t w e e na c c u m u l a t i o ne x p e r i m e n t u n d e rt h ei c es h e e t i n l a b o r a t o r yu c u r v ea n dn u m e r i c a ls i m u l a t i o n t h e r e s u l to fc o m p a r i s o n s h o wt h a t ：t h em a x i m u mi c e t h i c k n e s so c c u r r e d o no 。c r o s ss e c t l o n t h i c k n e s so fi c es h e e ti n c u r v ec o n v e xb a n k i n c r e a s e sg r a d u a l l yw l t n c r o s s s e c t i o na n g l ei n c r e a s i n g ，w h i l er e d u c i n gg r a d u a l l y i nc o n c a v eb a n k i n c r e a s i n gt h es p e e do fw a t e rf l o w u n d e rt h e i c es h e e tw i l l r e d u c et h e t h i c k n e s so fi c es h e e t k e yw o r d s ：i c ej a m ；n u m e r i c a l s i m u l a t i o n ；f l u e n t ；r i v e r i c ej a m f o r m a t i o n ；m e l t i n g 插图清单 2 1 基本结构程序1 3 3 1 a 不同入口流速条件下断面的温度分布3 3 3 1 b 水流速0 1 m s 时各断面温度3 5 3 1c 不同入口流速条件下的沿程冰厚3 5 3 2 a 不同空气温度条件下断面的温度分布。3 7 3 2 b 大气温度2 6 3 1 5 k 时各断面温度3 8 3 2 c 不同大气温度条件下的沿程冰厚3 8 3 3 a 不同空气温度条件下断面的温度分布3 8 3 3 b 风速1m s 各断面温度4 l 3 3 c 不同风速条件下的沿程冰厚4 2 4 1 不同水流速时的下壁面温度4 4 4 2 不同流速时各参考断面温度4 6 4 3 a 水流速o 4 m s 时冰塞的融化状态。4 6 4 3 b 水流速o 8 m s 时冰塞的融化状态。4 6 4 3 c 水流速1 2 m s 时冰塞的融化状态4 6 4 3 d 水流速1 6 m s 时冰塞的融化状态4 6 4 4 不同大气温度时的下壁面温度4 7 4 5 不同大气温度时的各断面温度4 9 4 6 a 大气温度2 7 1 1 5 k 时冰塞融化状态4 9 4 6 b 大气温度2 7 2 1 5 k 时冰塞融化状态4 9 4 6 c 大气温度2 7 3 1 5 k 时冰塞融化状态4 9 4 6 d 大气温度2 7 4 1 5 k 时冰塞融化状态4 9 4 7 不同风速时下壁面温度。5 0 4 8 不同风速下各断面温度5 2 4 9 a 风速1 m s 时冰塞融化状态。5 2 4 9 b 风速4 m s 时冰塞融化状态。5 2 4 9 c 风速8 m s 时冰塞融化状态5 2 4 9 d 风速1 2 m s 时冰塞融化状态5 2 5 1 实验计算模型图5 2 5 2 实验室试验所得冰厚与模拟所得冰厚的对比。5 5图图图幽图图图图图图图图图图图幽图图图图图图图图图图图图图图 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的研究 成果。据我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人已 经发表或撰写过的研究成果，也不包含为获得 佥巴王些厶堂 或其他教育机 构的学位或证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均 已在论文中作了明确的说明并表示谢意。 学位敝作者签名：乃 ， 7 角签字眺砷年午月毋 学位论文版权使用授权书 本学位论文作者完全了解金蟹_ i 些太堂有关保留、使用学位论文的规定，有 权保留并向国家有关部门或机构送交论文的复印件和磁盘，允许论文被查阅和借阅。 本人授权金巴工些太堂可以将学位论文的全部或部分内容编入有关数据库进行检 索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后适用本授权书) 学位论文作者躲乃磊 签字日期咿7 年￥月) 7 日 学位论文作者毕业后去向： 工作单位： 通讯地址： 导师签名： 签字日期：7 7 年争月筇日 电话：| 弓苫6 弓参器纪 邮编： 致谢 本文是在导师王军教授悉心指导下完成的。两年多来，导师严谨的治 学作风、深厚的学术造诣、渊博的科学知识、朴实的生活作风使我受益匪 浅。正是导师的言传身教、专业上的耐心指导、为人处世方面的谆谆教诲 和生活上无微不至的关怀，才使我得以顺利完成硕士研究生学业。论文凝 聚了导师的心血，值此论文完成之际，谨向恩师及其家人表示衷心的感i 勇十 和最诚挚的敬意! 同时感陈胖胖，盈尚轩在生活和学习上给予的多方面的帮助和支持， 在此向他们表示我心中最诚挚的感激和诩 意! 在此特别感谢我的父母，多年来正是他们的理解和无私奉献，为我创 造了一片宁静的学习环境，是他们在精神上物质上的极大支持才使我得以 安心完成学业。无以为报，谨以此文献给他们。 最后，再次向曾经关心、鼓励和培养我成长的老师、同学、朋友和家 人表示我最诚挚的祝福和最真挚的感谢! 作者：石磊 2 0 0 9 年4 月2 0 日 第一章绪论 1 1 冰情的概述 河流冰凌问题是冬季寒冷地区江河中普遍发生的一种水文现象，它的 产生、发展及消失都是一个极其复杂的过程。与冰凌形成变化有关的因素 包括：热力因素、水流过程的变化、河道形态、地貌状况、水流的水力学 特性和人类活动等，如修建水库、引水工程及其它水工建筑物等，在不同 程度上改变了冰情的自然变化规律。 我国是世界上冰凌灾害较为严重的国家之一。每年的冬季，大约四分 之三以上面积的河流、湖泊不同程度地存在凌情。有冰情的江河，在流冰 封冻、解冻时，往往会发生各种冰凌灾害【1 2 3 】。 1 2 冰凌灾害的几种型式 ( 1 ) 冰坝洪水 冰坝洪水是天然河流中最严重的冰凌灾害。冰坝是由大量的冰块在河 道中堆积而成的。冰坝形成以后，堵塞河道，大大减少过水断面面积、增 加水流阻力。流冰开始堆积的原因是由于河道障碍物、弯道、卡口、河道 束窄、断面变化等。冰坝使水位上涨，在某些河段流水漫堤，造成凌洪灾 由 3 】 口0 ( 2 ) 冰花堵塞 悬浮的冰花遇到过冷的固体时则贴附在外表，层层冻结，逐渐加厚， 减少甚至完全堵塞过水断面，如电站进水口拦污栅，使电站不能运行，同 时电站上游会因水位壅高漫出河堤形成凌洪灾害。 ( 3 ) 影响航运和建筑物安全 流动的冰块会产生很大的动冰压力和撞击力，碰撞船舶和其他建筑 物，使河流冬季无法通航，水工建筑物也会遭到破坏。 ( 4 ) 损坏岸坡和水工建筑物 在温度升高时，冰盖膨胀产生巨大的静冰压力使河岸护坡和水工建筑 物( 如进水塔、桥墩和胸墙等) 遭到破坏。静冰压力与气温条件关系密切， 在华北地区一般为l o t m 2 ，在东北、西北地区为1 5 2 0 t i m 2 ，个别寒冷地 区可达3 0 t m 2 。 1 3 国内外冰情的主要研究方法 目前，国内外主要从实地观测，实验研究，数值模拟等方面对冰水力 学的若干问题进行研究： 在原型观测方面，国内外学者进行了很多研究，如b e l t a o s 【4 5 1 、马 喜祥等【6 】、s u i ，j ve ta l 7 1 ，原型观测对了解特定场所下的河冰运动及规 律并用于相应区域的防灾减灾有积极的作用，但其缺陷在于：自然条件千 变万化，很难控制，难以理解各因子的作用和机理，通常得出的是经验关 系，一般难以移用其它场所。 试验研究的主要目的是希望借助于试验水槽定性或定量研究冰水力 学的若干问题，以探讨机理、寻求一般性的规律，意图建立相关的理论体 系。随着科技和经济的发展，河冰物模试验的手段和方法也日渐先进， 模型冰就很能说明问题，如前所述，最早开始是使用石蜡，其实较早的研 究还有使用木块等作为模型冰材料，现今较多的是用人工合成材料模拟 冰，例如，聚乙烯，泡沫聚苯乙烯，天然冰及掺了添加剂的真冰等。 河冰物模试验分为两大类：一是针对工程原型的比尺试验，另一个是 并不针对某一具体特定原型的物理过程模型试验。显然，两者的目的不同， 比尺试验面向具体原型，要模拟解决特定条件下的工程实际问题，而物理 过程模型试验是一般机理、一般规律或某些特定过程的研究。关于河冰模 型试验的相似关系已有很多学者进行研究，如u r r o z ，g e 和e t t a m a ，l 引、 孙肇初【13 1 、陈建国等【14 1 、s a a d e ，r a a f a t 和g s a r r a f ，s 1 5 】、d h e a l y 和f h i c k s 1 6 1 7 】、e t t a m a ，r 【1 8 】、王军 1 9 - 2 2 1 等通过水槽试验研究了冰塞形 成与演变的一般机理和相应的一些规律，单纯的水槽试验研究近年来已较 少。依据目前国内外对弯槽段冰塞试验研究的可见文献，可以得到这些研 究的共同特性是：基本是探索性的定性描述，系统性不够，试验场次很少， 并不足以较为完全的定性描述弯槽段冰塞形成的规律。 在数值模拟方面，主要运用数学模型和计算机技术进行研究，并结合 水槽试验结果或原型观测结果加以验证。随着计算机运算速度的提高，数 值模拟技术有了长足的发展，开辟了许多研究冰塞的新方法1 2 弘2 7 j 。 z u f e l te t a l 2 8 】建立了一维输冰耦合模型，该模型的基本方程分为冰、 水质量守恒方程和动量守恒方程，利用该模型，z u f e lte t a 研究了冰的动 量和水的非恒定性对冰塞厚度值大小及其分布的影响，得出了基于 b e lt a o s 2 9 】一维静态平衡冰塞厚度计算公式得出的厚度值明显偏小的结 论。 在一维模型应用的基础上，s h e n 和l u 3 0 】较为完整的推出了河冰数值 模拟的两维模型，用以模拟冰塞体的溃决。模型分为几个模块：水动力模 块、连续性方程、动量方程、冰块输运方程、冰质量守恒方程、冰面积浓 度守恒方程、本构定律等。对水动力模块，s h e n 和l u 采用有限元方法计 算，而对输冰方程则采用d p m ( d is c r e t ep a r c e lm e t h o d ) 求解方法。运 用模型分别对一假设的矩形断面明渠和n i a g a r ar i v e r 上游河段进行了模 拟，获得了一些规律性的认识。 2 近年来，国内河冰数值模拟亦取得了长足的进展。如杨开林p l j 等根据 冰塞形成发展的机理提出了冰塞形成的发展方程；吴剑疆等p - 2 j 对敞露河段 内水内冰花的体积分数以及水温的沿程分布进行了模拟研究，所得规律与 理论分析相符；茅泽育等【3 3 】针对天然河道弯曲复杂的特点，建立了适体坐 标下的二维河冰数值模型，经实地验证取得较好的效果等。国内基于神经 网络对河流冰情研究也取得了可喜的成果，如陈守煜等【34 。、王涛等【) 引、王 军等p 6 j ：另外，李亚伟等【37 j 尝试了s v r 方法。 现在许多学者正在运用计算流体力学( 简称c f d ) 的方法进行研究，其。 基本思想是：把原来在时间域及空间域上连续的物理量的场，如速度场， 压力场，用一系列有限个离散点上变量值的集合来代替，通过一定的原则 和方式建立起关于这些离散点上场变量之间的关系的代数方程组，然后求 解代数方程组获得场变量的近似值【3 8 】，本文应用的计算软件f l u e n t 就是 建立在这种思想基础上的，但数值模拟也有其自身的缺陷，譬如模拟时考 虑主要的影响因素而忽略次要的影响因素，这就会影响模拟的精度等，但 数值模拟的方法仍是我们研究的一种重要方法。 1 4 数值模拟中数值方法的分类 当前数值模拟的方法有很多种，但应用最为广泛的是有限差分法、有 限元法、有限分析法和有限体积法。 1 4 1 有限差分方法( f i n i t ed i f f e r e n c em e t h o d ，简称f d m ) 有限差分法( f d m ) 是计算机数值模拟最早采用的方法，至今仍广泛 运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求 解域【3 9 】。f d m 以t a y l o r 级数展开等方法，把控制方程的微商用差商代替进 行离散，从而建立代数方程组来求解。该方法数学概念直观，表达简单， 其解的存在性、收敛性和稳定性早已有较完善的研究成果，是比较成熟的 数值方法，目前应用最广。 由于实际应用中采用的时间和空间的差分形式不同，差分法又可分为 显式、隐式及显式一隐式交替等方法。显式差分格式是指任一网格节点上 的的待求因变量在新的时间层的值可以通过已知的时间层上的变量值显 式解出。显式差分法应用较早、简单，可避免试算，但为了保证其稳定性， 须严格遵守柯朗条件( c a t a x l ，c 为小扰动波的波速，c = 办+ u ) ， 时间步长和空间步长受到限制。隐式差分格式是指未知网格节点上的待求 变量不能由已知的时间层的函数值直接求出，还需同一时间层相邻节点函 数值( 未知) 作为信息，通过联立求解方程组才能得到未知量的解。从理 论上讲，隐式格式是无条件稳定，但在实际应用中，由于空间、时间步长 为有限量，其时间步长也有一定的限制。隐式差分格式的优点是时间步长 可以取得较大，稳定性能好，但在计算中需迭代，计算量较大。 交替方向隐格式法( a d i ) 4 0 】是由d o u q l a c e 和r a c h f o r d 等1 9 5 5 年提出 的，后来被l e e n d e r t s e 结合交替网格建立起来并首次用于计算平面二维 流场。a d i 方法是一种显一隐格式交替使用的有限差分格式，该方法同时 具有显式和隐式两种差分格式的优点，与完全隐格式相比较，它不必每一 个时间步骤都要求求解一个大型代数方程组，因而所需内存少，计算量也 相应减少。同时a d i 方法有较好的计算稳定性和计算精度，目前已广泛应 用于流体计算中。另外，在2 0 世纪7 0 年代初，前苏联学者y a n e n k o 等提 出破丌算子法【4 ，按维数或按照方程的性质对控制方程进行算子的破开。 但目前许多学者对破开算子法提出疑义，所以该方法应用不是很广。 在二维模型计算中，由于使用t a y l o r 级数展开，f d m 一般只适用于矩 形或正交曲线网格，在计算域概化和数值解精度方面，存在着根本性的困 难【4 2 1 。 1 4 2 有限元方法( f i n i t ee 1 e m e n tm e t h o d ，简称f e m ) 有限元法【4 3 】( f e m ) 产生于2 0 世纪5 0 年代，最先应用于固体力学， 6 0 年代开始在流体力学中有所应用。有限元法的基础是极值原理和剖分插 值，它吸收了有限差分法中离散处理思想，同时采用了变分计算中选择逼 近函数及对任意形状( 三角形或四边形) 的许多微小单元进行积分处理的 合理方法，因而具有很广泛的适应性，特别适合于几何、物理条件比较复 杂的问题。该方法具有较强的适应性，计算精度较高，但存在计算格式复 杂、计算量及贮存量较大，大型的系数矩阵较难求解等缺点。 常见的有限元计算方法有直接法、变分法、加权余量法及能量平衡法 等。其中变分法类的里兹法、加权余量法类的g a l e r k i n 法和最小二乘法 常用于河流的数值模拟。采用常规的有限元方法时，对于对流效应比较强 的情况，常常由于有限元网格不恰当而造成数值解的失真或震荡。在有限 元方法中通常时附加人工粘性和采用迎风格式，但这样使方程失去了加权 余量的意义，从而使近似方程不满足相容性，即在保证稳定的条件下，失 去了精度。后来又提出流线迎风p e t r o v g a le r k i n 有限元法，简称s u p g ， 该方法在稳定性、收敛性及精确度等方面已有很大的提高。 在多维数学模型计算中，因有限元方法贮存量比较大，且大型的系数 矩阵求解较困难等，直接影响着计算速度，因而在非恒定问题及其他对计 算速度要求较高的问题中应用不是很多【40 1 。 4 1 4 3 有限分析法( f i n i t ea n a l y t i cm o t h o d ，简称f a m ) 有限分析法是在有限元基础上的一种巧妙的改动，其核心是：在离散 单元上的解是采用求解方程本身在局部线性化后单元上的解析解，不像有 限差分法只限于网格点上的值，也不像有限元那样用插值函数式来表达。 该方法是2 0 世纪8 0 年代美籍华人( c h i n g j e nc h e n ) 教授等提出的一 种新的数值方法【3 9 】。 f a m 有较高的计算精度，并具有自动迎风特性，计算稳定性好，收敛 较快，在我国有了较多的研究和应用。但由于有限分析系数中含有交错级 数，这给实际计算与理论分析都带来了一些困难1 4 引。 1 4 4 有限体积法( f i n i t ev o l u m em e t h o d ，简称f v m ) 有限体积法( f v m ) 又称有限控制容积法，该方法与有限元方法一样， 把计算区域离散分为若干点，以这些点为中心，把整个计算区域分为若干 互相连接但不重叠的控制体。在有限体积法计算中，对每一控制体分别进 行能量和动量平衡计算，便得到一组以控制体特征量平均的物理量为未知 数的代数方程组，同时沿坐标方向对方程组进行离散，形成的离散方程与 有限差分法有些相似。因为跨控制体间界面运输的通量，对相邻控制体来 说，大小相等，方向相反，故对于整个计算域，沿所有内部边界的通量相 互抵消，对由一个或多个控制体组成的区域，以及整个计算区域，都严格 满足物理守恒定律，不存在守恒量的误差。 有限体积法中，若采用相邻控制体形心处通量平均，在矩形网格的计 算中，便相当于二阶中心的有限差分法，而在三角形或四边形网格中，若 物理量定义在网格顶点，则又与线性三角形和双线性四边形单元的 g a l e r k i nf e m 等价。若采用特征逆风格式计算通量，有限体积法适于处理 对流占优的输运问题，且在矩形网格上相当于守恒逆风有限差分格式。 因为有限体积法从物理规律出发，每一个离散方程都是有限大小体积 上某物理量的守恒表达式，在推导过程中物理概念清晰，并可以保证离散 方程的守恒特性，同时该方法能像有限元一样适用于不规则网格和复杂边 界情况，且处理效率与有限差分法相似，远高于有限元方法，在数值模拟 中有着很大的发展潜力【3 2 】。本文应用的计算流体力学软件f l u e n t 就是基 于有限体积法来模拟计算的。 1 5 流体流动的控制方程 水流的运动要受物理守恒定律的支配，基本的守恒规律包括质量守恒 定律、动量守恒定律、能量守恒定律。由于水流作为不可压缩流体，没有 密度变化就没有能量方程和质量方程及动量方程之间的联系，往往只用质 量守恒方程及动量方程求解流场。下面就控制流体运动的基本方程【3 8 1 予以 介绍。 1 5 1 连续性方程 连续性方程【3 9 1 是质量守恒原理在流动问题中的表述，对于不可压缩流 式中，甜( 或“，v ，w ) 分别为( 或x ，y ，z ) 三个方向的流速分量。 1 5 2 动量方程 动量方程是牛顿第二定律或动量守恒定律在流体流动问题上的表述， 分析微小控制体受力和运动情况，可得到不可压缩牛顿流体在z ，y ，z 方向应 满足的方程，其张量形式和微分形式为： 一a u i + 墼：只一土鱼+ 刁旦( 1 - 2 a ) a ta ) c i i pa x , i 瓠j xi 丝+ 望+ 塑+ 塑：x 一土至+ 7 7 v z 甜( 1 2 b ) 挑礅 动 a z p o x 堡+ 塑+ 竺+ 塑：y 一! 望+ 7 7 v z l ，( 1 2 c ) 硗瓠 融 a z pa y 塑+ 塑+ 塑+ 笙：z 一土鱼+ v z w( 1 2 d ) d t礅 谚现p o z 式中，f = ( x ，】，z ) 为单位质量力沿x ，y ，z - - 个方向的分量；p 为压强； 运动粘度叩= p 。( 1 2 a 、b 、c 、d ) 式以1 9 世纪两位独立推导此方程科 学家名字命名，即著名的n a v i e r s t o k e s 方程( n s 方程) 。 1 5 3 紊流模型 ( 1 1 ) 一( 1 2 a 、b 、c 、d ) 式是实际液体运动的微分方程式，既适用于 层流也适用于瞬时紊流。但因紊流的瞬时运动要素有脉动现象，要对高度 复杂的紊流运动进行直接的数值计算，必须采用很小的时间与空间步长， 才能分辨出紊流中详细的空间结构及变化剧烈的特性。现代计算机的存储 能力和运算速度尚不足以求解任何一个实际的紊流问题，因此研究紊流运 动，可求出紊流的时均运动方程式及连续性方程式【4 4 | 。 6 按照r e y n o l d 平均法，任一变量矽的时间平均值的定义为： 万2 挚删矽o ) a t ( i - 3 ) 物理量的瞬时值矽，时均值歹及脉动值之间有下列关系： 矽= 歹+ ( 卜4 ) 对于n a v i e r s t o k e s 方程作时均演算，在直角坐标系下导出时均的 控制方程为： 连续性方程： 翌型2 + 型+ 塑型2 ：堡 i l - 堡+ 翌+ 堡+ 翌+ 堡：o( 1 5 ) 一+ 一+ 一= 一一十一- t 一一r 一1 一一= u 、1 u ， 缸 砂砂苏砂出 出 砂 毖 显然有： i o h - + 娶+ 罢：o 眺铆 o z a 万万细 + + = u 瓠 却 跳 动量方程 以z 方向的动量方程为例，作类似于上面的处理，有： 堡+ 旦盟+ 旦显+ 旦婴+ 垫塑+ 垫塑 o to x a y o z 砂 o z ( 1 7 ) = x 一吉罢州窘+ 雾+ 窘， 把上式左端脉分量乘积的时均值项移到等号的右端，得： 堡+ 塑+ 塑+ 塑 一a t 净o x 童嚷新特u , v - - - - i 张一丽 8 ) = x 一吉罢+ 去 叩芸一两 + 专号一 + 丢 刁鲁一丽 其他两个方向的方程也可以作类似的推导，且采用b o u s s i n e s q 关于 紊流粘性系数的假设，即把湍流应力表示成湍流粘性系数的函数，可得到 紊流时均运动的动量方程，即r e y n o ld 方程，写成张量形式： 7 1 一一 孕铷孚警 互+ a ( u , u j ) ：f 一三至+ 仇盟 ( 1 9 ) 现 瓠j j p 瓠i 。瓠j 0 x j 上式中仍为紊流粘性系数。引入b o u s s i n e s q 假设以后，计算紊流流动 的关键就在于如何确定仇。所谓紊流模型也就是把仍与紊流时均参数联系 起来的关系式。根据确定绣方法的不同，可以将紊流模型区分为所谓的 零方程模型、一方程模型及两方程模型等。 ( 1 ) 零方程模型【3 9 】： 所谓零方程模型，是指不需要微分方程而是用代数关系式把湍流粘性 系数与时均值连系起来的模型。确定砚的方法有两种，一个是仍为常数； 二是采用p r a n d t l 混合长度模型。 p r a n d t l 按照分子运动论自由长度的概念，提出紊流中涡团混合长度 的概念，认为紊动粘性系数可以写为： 功：焉剖 ( 1 - l o ) j 哕i 其中，“为主流的时均速度；y 是主流方向相垂直的坐标；乙为混合 长度，是这种模型需要确定的参数，而该值的确定比较复杂，常与流动的 类型有关。混合长度理论是基于紊流局部平衡的假定基础上的，在物理概 念上有一定的不足。 ( 2 ) 一方程模型 在混合长度理论中，仇的确定仅与几何位置及时均速度场有关，而与 湍流的特性参数无关。湍流粘性系数应当与湍流本身的特性量有关。 p r a n d t l 从上述考虑出发，提出了计算纭的下列表达式： 研= c 。，。肚2 ， ( 卜1 1 ) 其中c ，。为经验系数，是湍流脉动的长度标尺，一般的它不等于混合 长度乙。采用此式来确定仇时，关键在于确定流场中各点的脉动动能及长 度标尺。 紊动能的方程为 夕i o k + p u j 婺：l ( 叩+ 要) 婺l + 铳挈冬+ 要) 一p 竽 ( 1 也) 夕百 瓦2 瓦1 研+ 毒虿j 托茜嚼+ 司嘞p 丁 u 1 动 这里吼称为脉动动能的p r a n d t l 数，其值在1 o 左右。系数c d 在不同 的文献中取值不定。对于，湍流脉动的长度标尺，常用的方法是采用类似 于混合长度理论中的，计算式。 在一方程模型中，湍流粘性系数与能表征湍流流动特性的脉动动能联 系了起来，这无疑优于混合长度理论。但是在一方程模型中仍要用经验的 方法规定长度标尺的计算公式，这是一方程模型的主要缺点。 ( 3 ) 两方程模型 在一方程模型中，湍流的长度标尺，是由经验公式给出的。其实长度 标尺也是一个变量，可以通过求解偏微分方程而得出。由于对靠近壁面地 区的计算而言，s 方程的计算最为方便。( f 值为紊流中单位质量流体脉 动动能的耗散率，即各向同性的小尺度涡的机械能转化为热能的速度。) 因而，在湍流的工程计算中，k s 的两方程模型应用最广。 紊动能k 和紊动能耗散率占的方程依次为 警+ 叶考2 毒ic 刁+ 毒，考卜最吖 n 1 3 ， 害栅专2 苦ic 刁+ 善，考j + 妻c g 最一g s ， c 卜m ， e = 一咖：导 ( 卜1 5 ) a x j r l , = c k 2 f ( 1 - 1 6 ) 在k 一占模型中引用了5 个模型常数，其值依次为c 。= 0 0 9 ，c 1 = 1 4 4 ， c ：= 1 9 2 ，吼= 1 0 ，仃。= 1 3 。上述模型以各向同性紊动粘性为基础，又 称为标准的k s 模型。上述模型又称为标准的k 一占模型，是以各向同性紊 流黏性为基础的，因此对各向异性紊流模拟的预报能力较差。两方程的紊 流模型是所有紊流模型中最简单的完全模型。 ( 4 ) r e y n o l d s 应力模型 上面介绍的各种紊流模型中都采用了各向同性的湍流动力粘度来计 算湍流应力，这些模型难以考虑流动方向表面曲面变化的影响。我国著名 科学家周培源早在2 0 世纪4 0 年代就导出世界上第一个计算湍流应力的l7 方程模型【45 1 ，直接对湍流脉动应力一p u ；u , 9 t 行r 计算。但限于当时计算资 源的影响，并未得到广泛应用。就目前而言，r e y n o l d s 应力模型已开始用 于工程数值计算的是二阶矩模型以及在此基础上经简化而得出的代数应 力模型，不少研究者认为这是目前最有发展前途的湍流模型 4 6 】。 1 6 课题的来源、目的及意义 本课题是国家自然科学基金项目一弯槽段冰塞演变机理的试验研究 与数值模拟( n o 1 0 3 7 2 0 2 8 ) 研究中有关天然河道冬季水位、冰盖冰塞厚 度部分的研究。 冰塞是冰和水在河道运动的产物，由于其错综的复杂性使得冰塞水位 9 的研究仍处于初步阶段。研究冰塞形成、演变及消融规律的最终目的，是 要利用这些规律来预测冰塞河段的壅水幅度，减轻或避免冰塞的危害。 本文主要应用计算流体力学软件f l u e n t 模拟冰的堆积及凝结过程， 分析一般规律，论文共分六部分。 第一部分是绪论，主要对本文要做的工作进行简要的介绍； 第二部分简要介绍了计算流体力学软件f l u e n t 的基本理论及应用， 对它的网格适应、可以计算的物理问题类型、边界条件、湍流模型、求解 器加以说明： 第三部分应用f l u e n t 软件中的凝固与融化解算器模拟模型中的冰盖 形成状况，并对模拟结果进行了比较分析； 第四部分应用f l u e n t 软件中的凝固与融化解算器模拟模型中的冰盖 消融状况，并对模拟结果进行了比较分析； 第五部分将实验室实验所得冰塞厚度的数据与数值模拟所得数据进 行比较，分析两种方法的异同； 第六部分是结论与展望，对本文研究的问题进行综述，并对所研究问 题的研究前景做了展望。 1 0 第二章f l u e n t 基本理论及应用介绍 计算流体软件f l u e n t 是用于模拟具有复杂外形的流体流动以及热传 导的计算机程序。它提供了完全的网格灵活性，可以使用非结构网格，例 如二维三角形或四边形网格、三维四面体六面体金字塔形网格来解决具 有复杂外形的流动，甚至可以用混合型非结构网格。允许根据解的具体情 况对网格进行修改( 细化粗化) 。同时也可以纳入p a t r a n ，a n s y s ，i d e a s 及i c e m c f d 等专门生成网格的软件。可以计算的物理问题类型有：定常与 非定常流动，不可压缩与可压缩流动，含有粒子，液滴的蒸发，燃烧的过 程，多组份介质的化学反应过程等。1 9 9 8 年f l u e n t 公司推出了自己研制 的新的前处理网格生成软件g a m b i t ，并且将f l u e n t u n s 与r a m p a n t 合并 为f l u e n t 5 ，本文应用的是f l u e n t 6 。 2 1f l u e n t 软件的简介 2 1 1 程序结构 f l u e n t 是用c 语言编写的，因此具有很大的灵活性。f l u e n t 可以应 用多种形式的网格，我们可以用g a m b i t 直接产生所需的几何结构以及网 格，然后将网格导入f l u e n t 中进行计算，如果所研究问题的几何结构比 较复杂，直接在g a m b i t 中产生较为困难，也可以在c a d c a e 等绘图软件 中产生几何结构，然后导入g a m b i t 中进行网格划分。也可能用其他软件 产生f l u e n t 所需要的网格，比如a n s y s ，i d e a s ( s d r c ) 或 m s c a r i e s 。m s c p a t r a n 以及m s c n a s t r a n 。一旦网格被读入f l u e n t ，剩下 的任务就是使用解算器进行计算了。其中包括边界条件的设定，流体物性 的设定，解的执行，网格的优化，结果的查看与后处理，如图2 1 所示。 2 1 2f l u e n t 的模拟能力 f l u e n t 解算器应用于流体力学里有如下模拟能力： 用非结构白适应网格模拟2 d 或者3 d 流场，它所使用的非结构网格主 要有三角形四面体、四边形六面体，或者混合网格，其中混合网格有棱 柱形和金字塔形；不可压或可压流动；定常状态或者过渡分析；无粘，层 流和湍流；牛顿流或者非牛顿流；两相流，包括气穴现象；复杂外形的自 由表面流动等。 对于模拟复杂流场结构的不可压缩可压缩流动来说，f l u e n t 是很理 想的软件。 2 1 3 解决问题的步骤 应用f l u e n t 进行数值模拟时，确定所解决问题的特征之后，需要以 下几个基本的步骤来解决问题： ( 1 ) 创建网格 ( 2 ) 运行合适的解算器：2 d ( 二维单精度) 、3 d ( 三维单精度) 、2 d d p ( 二 维双精度) 、3 d d p ( 三维双精度) ( 3 ) 输入网格 ( 4 ) 检查网格 ( 5 ) 选择解的格式 ( 6 ) 选择需要解的基本方程：层流还是湍流( 无粘) ( 7 ) 指定材料物理性质 ( 8 ) 指定边界条件 ( 9 ) 调节解的控制参数 ( 1 0 ) 初始化流场 ( 1 1 ) 计算解 ( 1 2 ) 检查结果 ( 1 3 ) 保存结果 ( 1 4 ) 必要的话，细化网格，改变数值和物理模型 图2 - 1 基本程序结构 1 2 2 2f l u e n t 中的网格 2 2 1 网格类型 f l u e n t 是非结构解法器，它使用内部数据结构来为单元和表面网格点 分配顺序，以保持临近网格的接触。f l u e n t 采用非结构网格以缩短产生网 格所需要的时间，简化了几何外形的模拟以及网格产生过程。和传统的多 块结构网格相比，它可以模拟具有更为复杂几何结构的流场，并且具有使 网格适应流场的特点。f l u e n t 也能够使用适体网格，块结构网格( 比如 f l u e n t4 和许多其它的c f d 解算器的网格) 。f l u e n t 可以在2 d 流动中处 理三角形网格和四边形网格，在3 d 流动中可以处理四面体网格，六面体 网格，金字塔网格以及楔形网格。f l u e n t 可以接受单块和多块网格，以及 二维混合网格和三维混合网格。另外还接受f l u e n t 有悬挂节点的网格( 即 并不是所有单元都共有边和面的顶点) 。这种灵活处理网格的特点使我们 在选择网格类型时，可以确定最适合特定应用的网格拓扑结构。 2 2 2 网格类型的选择 网格的选择依赖于具体的问题，在选择网格的时候，需要考虑下列问 题： 1 初始化的时间 很多实际问题是具有复杂几何外形的，对于这些问题采用结构网格或 块结构网格可能要花费大量的时间，甚至根本无法得到结构网格。复杂几 何外形初始化时间的限制刺激了人们在非结构网格中使用三角形网格和 四面体网格。 2 计算花费 当几何外形太复杂或者流动的长度尺度太大时，三角形网格和四面体 网格所生成的单元会比等量的包含四边形网格和六面体网格的单元少得 多。四边形和六面体单元的一个特点就是它们在某些情况下可以允许比三 角形四面体单元更大的比率。三角形四面体单元的大比率总会影响单元 的歪斜。因此，在相对简单的几何外形，而且流动和几何外形很符合的情 况下，比如本论文主要涉及的u 型弯道，就可以使用大比率六面体单元。 这种网格可能会比三角形四面体网格少很多单元。 3 数值耗散 多维条件下主要的误差来源就是数值耗散又被称为虚假耗散。 关于数值耗散有如下几点：( 1 ) 真实耗散很小时，即对流占主导地位 时，数值耗散是显而易见的。( 2 ) 所有的解决流体问题的数值格式都会有 数值耗散。( 3 ) f l u e n t 中所用的二阶离散格式可以帮助减少解的数值耗散 的影响。( 4 ) 数值耗散量的大小与网格的分辨率成反比。因此解决数值耗 散问题的一个方法就是精化网格。( 5 ) 当流动和网格成一条直线时数值耗 散最小。( 6 ) 使用三角形四面体网格流动永远不会和网格成一条直线，而 如果几何外形不是很复杂时，四边形网格和六面体网格可能就会实现流动 和网格成一条线。 2 2 3 网格质量 网格质量对计算精度和稳定性有很大的影响。网格质量包括：节点分 布，光滑性，以及歪斜的角度( s k e w n e s s ) 。 ( 1 ) 节点分布 边界层解( 即网格近壁面间距) 在计算壁面剪切应力和热传导系数的 精度时有重要意义。这一结论在层流流动中尤其准确，网格接近壁面需要 满足： f 少p 、丝1 ( 2 1 ) ye x 其中： j ，。一从临近单元中心到壁