（热能工程专业论文）喷管内气粒两相流场的数值模拟.pdf

哈尔滨工程大学硕士学位论文 第1 章绪论 1 1 喷管内气粒两相流动数值模拟的重要性 所谓相，通常是指某一系统中具有相同成分及相同物理、化学性质的均 匀物质部分。两相流动则通常指的是由两种相组成的流动，如流体相和颗粒 相组成的流动。喷管内的两相流动广泛用于工业、航天等领域，如固体火箭 发动机喷管内的流动，粒子加速喷管内的流动等。以固体火箭发动机喷管为 例，我们来看一看对喷管内气粒两相流动进行数值模拟的重要性。 固体火箭发动机因其结构简单，维护简便，操作使用方便，可靠性高， 长期储存性好，作为一种动力装置，在军事、航空航天等领域内得到了广泛 的应用。它的工作过程为：作为发动机能源的推进剂被点燃后，在燃烧室中 燃烧，经过复杂的物理变化和剧烈的化学反应过程，生成高温( 2 0 0 0 3 0 0 0 k 以上) ，高压( 几个到几十个m p a ) 的燃烧产物，将推进剂所蕴藏的部分化学 能转换为燃烧产物的热能。燃烧产物从燃烧室流入喷管，通过这个特殊形状 的管道膨胀加速，使其流速由亚音速转变为超音速，并从喷管中高速喷出， 将燃烧产物的热能转换为射流的动能，且高速气流从喷管中喷出时会产生直 接反作用力推力，推动飞行器运动，对飞行器做功，并转换为飞行器的飞 行动能。从固体火箭发动机的工作过程我们不难看出燃烧产物在喷管中的流 动是发动机运行时一个重要的能量转换过程。而且现代固体火箭发动机广泛 采用含金属复合推进剂，特别是含铝复舍推进剂，以达到提高比冲的目的， 并在一定程度上抑制高频振荡不稳定燃烧。铝的燃烧产物氧化铝粒子及由于 不完全燃烧而形成的液态铝粒子加入燃气后将形成典型的两相流动，会对发 动机性能产生较大影响。 喷管两相流会造成发动机比冲降低，因为： ( 1 ) 由于凝相粒子的布朗运动比气相要小的多，可以忽略不计，所以它不提 供压强，在流动中不作膨胀功，从而使含有凝相粒子的燃烧产物作的膨 胀功减少。 ( 2 ) 凝相粒子是借助于气体与粒子间的粘性摩擦力，靠气体带动它运动，这 样粒子的速度总小于气相速度，存在速度滞后。这种具有速度滞后的凝 相粒子对气体运动产生阻力，从而使气体的排气速度减小，比冲降低。 将热能传给气体，因此粒子的温度总是高于气体温度，存在温度滞后。 在粒子与气体之间还来不及达到热平衡就从喷管中排出，使一部分热能 没有转换为气体的动能，从而降低了喷管的能量转化效率。 大量的计算和实验证实在含铝推进剂的发动机中，喷管两相流损失约占 喷管中总损失的1 2 到t 3 。且在喷管两相流动中，高速运动的凝相粒子可 能对喷管壁面产生冲刷，一方面会对喷管的传热增大，加剧喷管的烧蚀；另 一方面，由于粒子的冲刷作用会引起壁顽粗糙度增大或者壁面的轻微剥蚀， 前者使喷管摩擦损失增加，后者如果发生在喷管喉部将使发动机的推力一时 间曲线下降。 所以准确掌握固体火箭发动机喷管内两相流动情况，对于减少两相流对 发动机性能的影响，提高发动机的设计水平具有重要意义。故本研究课鼷将 研究颗粒在一种小扩张比喷管内的流动，对其物理参数和物理现象加以预测， 希望能为研究装有此类喷管的火箭发动机性能参数和颗粒加速喷管的特性提 供一定的理论基础。 1 2 国外的发展和研究概况 在计算机技术和计算流体力学迅猛发展势头的推动下，在单相流模拟的 基础上，气粒两相流动的数值模拟也相应地开展起来【3 】。这是因为许多实际 的流体流动是气粒两相流动，例如云和雾、含灰尘的空气、粉尘分离与收集、 液雾喷涂、煤粉燃烧，等等。 以建立在解析研究方法之上的经典两相流体力学为基础。描述气粒流动 的颗粒相模型不断得到改进。最早期的单颗粒动力学模型( s p d ) ，它考虑已知 流场( 常常是均匀的速度场和温度场) 中颗粒的时平均运动或对流运动的轨 道，颗粒速度和温度沿轨道的变化，忽略颗粒对流场的影响。5 0 年代f u c h s 描述的气溶胶力学，6 0 年代航空发动机燃烧室中探讨的液一气两相流浓度场 问题，甚至8 0 年代末期国内外某些探讨选矿及粮食输送的两相流理论，大致 上都基于这类观点。这是一种极端的简化模型。差不多与此同时，t c h e n ，h i n z e 等人又由单颗粒在脉动流场中的行为出发探讨了颗粒湍流扩散问题。到6 0 年代后期及7 0 年代初期提出了单流体模型或无滑移模型( n s ) 该模型假设颗 粒和流体的速度与温度在空闯中处处相等，而颗粒扩散则相当于流体组分的 扩散( 扩散平衡) ，把颗粒与流体作为统一的流体来加以研究，这是另一种极 哈尔滨工程大学硕士学位论文 端的简化模型。由6 0 年代后期开始，s l s o o 开始提出用多连续介质即颗粒 群的小滑移拟流体模型来描述两相流，其中对稀疏悬浮流虽忽略颗粒对流体 的作用，然而对颗粒则由单颗粒的描述过渡到颗粒群拟流体概念。只考虑颗 粒群集体速度与温度，承认颗粒与流体间的滑移，但认为颗粒的滑移与其扩 散漂移是同一现象的两个方面，可称为小滑移模型( s s ) ，该模型基于颗粒连续 介质假设，模型中假设颗粒的运动单纯地由流体运动引起，流体与颗粒间的 速度滑移相对于平均流动来说是小量，可看作是颗粒扩散的结果。小滑移模 型是建立较完整的两相流模型的开端。 由7 0 年代中期以后逐渐发展了更完整的两相流模型，即完整地考虑相间 速度与温度滑移( 包括非起源于扩散的大滑移) ，颗粒扩散，完整地考虑相间 耦合，包括质量，动量和能量间耦合，特别是颗粒对流体的作用，以及颗粒 反应的模型，其中又主要分成颗粒相轨道( 包括不考虑扩散的确定轨道及考 虑扩散的随机轨道) 模型及颗粒相拟流体模型两大类。 关于喷管中两相流的研究比较受到重视，研究的成果比较多实验研究方 面有：h o g l a n d ( 1 9 6 4 ) 首先对气粒喷管流动开展了实验研究，并在实验研究 的基础上对其进行了一维近似求解。l e w i s 和c a r s o n ( 1 9 6 4 ) 通过将含颗粒 的气体由喷管射入气压室的实验，研究了颗粒装填比对气体自由射流结构的 影响，他们发现颗粒装填比越大，马赫盘向上游移动得越明显。对于二维无 粘气粒流动，k l i e g e l 和n i c k e r s o n ( 1 9 6 2 ) 指出在喷管的超音速区域，所要求 解的方程是双曲线型的，可用特征值法进行求解。s i n g l e t o n ( 1 9 6 5 ) 将 m a r b l e ( 1 9 6 3 ) 对含尘流体掠过半无限大平板的不可压层流边界层流动的分析 扩展到可压的情况。f e v e l l ( 1 9 7 4 ) 修改了一维模型。考虑了壁面表面摩擦及 对壁面传热的影响，使用喷管“核心流动”值来计算伴有烧蚀侵蚀喷管内的 两相流动。 数值计算方面有：c h a n g 等人( 1 9 9 6 ) 将欧拉方程和薄层纳维一斯托克 斯方程用于双流体模型，形成了时阀相关数值算法，可用来研究由于粒子的 影响而使喷管产生失效性的问题。对于长尾喷管，m e h t a ( 1 9 9 8 ) 采用双流体 法，将有限容积法和龙格库塔法相结合计算了喷管的粘性两相流流动。 c h a n g ( 1 9 9 6 ) 将矢通量分裂法和双流体模型相结合，计算了喷管跨声速两相无 粘流动。g c a r r i e r ( 1 9 9 0 ) h l 采用极为简化的准稳态逆流模型，计算了星孔 形装药燃烧室和潜入喷管内的气相流场，用拉格朗日法跟踪求解了颗粒相的 运动。 在用欧拉一拉格朗日法计算固体火箭发动机两相流的研究工作中， g o l a f s h a n i d ( 1 9 8 9 ) 【5 】的工作比较有代表性，他采用欧拉一拉格朗日法计算 了j p l 喷管无粘两相流场及两种固体火箭发动机后封头的二维轴对称粘性可 压两相流场。气相采用b e a m - w a r m i n g 的隐式近似因子分解法求解二维轴对称 n - - s 方程，对流项的离散采用了v a nl e e r 的矢通量分裂技术。颗粒相的求 解采用拉格朗日法沿粒子轨迹对运动方程进行数值积分。在g o l a f s h a n i 的研 究中引入了计算粒子的概念，即用有限数量的计算粒子来代表流场中所出现 的数量巨大的真实粒子，具有相同特性( 如位置、速度、温度、尺寸、质量 及组分) 的一组粒子可以用一个计算粒子来代表。气相和颗粒相之间的耦合 计算方法是完全的显式方法，在进行气相计算时，与颗粒相有关的耦合项假 定已知且在计算过程中保持常数。同样，当对颗粒相的动力方程进行积分求 解时，耦合项也处理成常数。当两种计算进行切换时，相闯耦合项就进行更 新，这两种计算的切换重复进行，直到达到稳态为止。为模拟粒子与固体壁 面的相互作用，给出两个恢复系数，分别代表粒子沿法向和切向的动量损失。 恢复系数等于1 时代表完全弹性碰撞，恢复系数等于0 时代表完全非弹性碰 撞( 粒子的动量完全损失并且粘在壁面上) 。为了模拟粒子穿越对称轴而到达 另一半，允许粒子与对称轴发生完全弹性碰撞。针对j p l 喷管进行了全耦合 无粘两相流动的计算。对粒子和壁面的碰撞，模拟了两种极限情况( 即完全 非弹性碰撞和完全弹性碰撞) 。计算中粒子质量百分数取3 0 ，粒子直径分 别为1 、2 0 、5 0 、1 0 0 um 。计算结果表明：1 ) 对于粒子与壁面为完全非弹性 碰撞的情况。1um 的粒子的随流性很好，粒子基本没有与喷管收敛段发生碰 撞，喷管扩张段上也没有出现无粒子区，2 0 um 以上粒子的惯性比较大，一 部分粒子与喷管收敛段发生碰撞而粘在壁面上，无法穿过喷管喉部，因此在 喷管扩张段粒子的轨迹趋向于向对称轴靠拢。这样，在喷管扩张段出现了无 粒子区。随着粒子直径的增大，喷管扩张段的无粒子区逐渐扩大。2 ) 对于完 全弹性碰撞的情况，粒子与喷管壁面之间可以发生完全弹性碰撞，这种碰撞 使得直径为1 0 0 t j m 粒子中的一些穿过中心线，使直径为5 0 i i m 的粒子集中 在粒子的极限流线附近。真实的两相流动，大粒子与壁面碰撞并非是理想的 完全弹性碰撞，也非理想的完全非弹性碰撞，而是非弹性碰撞，所以这种计 算实际上是考虑了两种极端情况。 m a d a b h u s h i ( 1 9 9 1 ) 旧求解了具有潜入喷管的固体火箭发动机后封头的 两相流场，采用可压缩的线性化块隐式格式求解气相流场，湍流模型采用低 4 哈尔滨工程大学硕士学位论文 雷诺数下的k e 两方程模型，用拉格朗日法模拟两相流中凝相粒子的运动轨 迹。他还研究了粒子加入速度在1 至2 5 之间变化时，对粒子运动轨迹的 影响，结果表明粒子的运动特性对粒子初速度的变化比较敏感。 s a b n i s ( 1 9 9 2 ) 啦基用轨道法计算了使用a p h t p b a 1 推进剂的s r m 多 相反应流场。凝相产物的燃烧速率采用了一个经验模型进行计算，连续相的 化学反应按化学平衡过程处理。s a b n i s 也引入了计算粒子，菇且粒子是在计 算坐标下进行追踪计算的。这种方法使得对粒子的追踪及粒子边界条件的处 理变得简化。他的研究中认为凝相粒子垂直于推进剂表面加入，并假定粒子 的加入速度是推进剂表面气相加质速度的l 。 m a w ( 1 9 9 3 ) 【b 】采用f d n s 代码对a s r m 后封头的有化学反应的两相流场 进行了求解，f d n s 采用欧拉一拉格朗日技术，气相控制方程中对流项采用二 阶迎风格式进行离散，粘性项和源项用二阶中心差分离散，对于几组不同粒 径的粒子用一步隐式方法来计算其轨迹，气相和颗粒相进行耦合计算。之后， l i a w 完善了这种模型。在f d n s 代码的基础上加入了颗粒蒸发，燃烧，破碎 模型( 1 9 9 4 ) 咖，并且采用随机轨道法来模拟湍流扩散对颗粒相的影响。 s m it h k e n t ( 1 9 9 3 ) 采用一种基于位势流的两相流模型来预示发动机 的熔渣沉积量。他将气相流场近似为位势流( 无粘，无旋) ，这样气相只求解 位势方程。颗粒相仍然用拉格朗日法计算。在他的模型中颗粒的尺寸，回流 区特性的描述，以及捕获规则等都是基于实验及经验数据。这种模型不考虑 喷管背壁回流区的具体形态，而是根据冷流实验的结果假定一个回流区与主 流的分界线，并且用固体边界来代替这个分界线。其他一些研究者也采用了 这种基于位势流的两相流模型来预示熔渣沉积，】。 位势模型虽然简便且计算省时，但对流场的描述过于简化，而求解湍流 全n s 气相控制方程既复杂又耗时，如果气相与颗粒相进行耦合计算，耗时 非常大。考虑到这种情况s a l i t a ( 1 9 9 4 ) t 4 1 发展了一种基于涡一流函数法的 代码叫做v o r s t r e m ，这种代码计算速度比较快，可以很方便地计算非稳态粘 性不可压流。s a l i t a 在v o r s t r e m 代码的基础上加入用拉格朗日法跟踪粒子 运动的代码t d 2 p ，再加上网格生成代码e g g ，形成了一套比较完整的计算固 体火箭发动机两相流的模型“e v t ”邮1 6 】( e g g - 贴体网格生成代码：v o r s t r e m - 涡流函数法求解不可压层流n s 有限差分代码；t d 2 p - 拉格朗曰法跟踪粒子 代码) 。 9 0 年代以后，由于大型火箭发动机普遍采用含铅复合推进剂及潜入喷 哈尔滨工程大学硕士学位论文 管，潜入喷管背壁熔渣沉积问题越来越受到重视，国外很多学者对这个问题 进行了深入的研究，并且取得了一些研究成果。1 9 9 5 年第3 3 届国际联合推 进会议上专门将熔渣沉积列为一个专题进行研讨，并且在这方面发表了很多 学术论文。这些研究工作基本都是在欧拉一拉格朗日方法的基础之上，进一步 发展和完善已有的模型，并提出了一些用于描述熔渣沉积过程的模型，例如 颗粒的捕获规则等。其中比较有代表性的工作有s a l i t a t ”】，c h a u v o t e ”1 ， l i a 讲- 一1 等人的工作。 s a l i t a ( 1 9 9 5 ) 采用“e v t ”模型计算了8 种s r m 燃烧室的两相流，预测 了熔渣沉积，并且与发动机试车时用r t r 技术观测到的熔渣沉积过程进行了 比较。 c h a u v o t ( 1 9 9 5 ) 计算了三种发动机的熔渣沉积，他对湍流模型，湍流扩 散效应以及捕获规则进行了讨论。 l i a w ( 1 9 9 5 ) 进一步发展了原有的模型，将其用于预测r s r m 发动机的熔 渣沉积。他的模型仍然综合考虑了气相的化学反应，颗粒的燃烧，蒸发，破 碎。计算结果表明由于燃烧，蒸发及破碎的影响，使得a i a 1 。o 。粒子在运动 过程中尺寸会不断的减小。为了模拟处于熔融状态的熔渣的沉积过程，他的 代码中引入了v o f ( v o l u m eo f f l u i d ) 模型，以捕捉熔渣自由面的运动。在 不考虑颗粒的燃烧，蒸发，破碎的情况下，他还计算了3 - dr s r m 的粒子轨迹 及熔渣沉积。 s a l i t a 【2 1 1 ( 1 9 9 5 ) 在他的文章中比较系统全面的分析了现有的s r m 熔渣 沉积模型，指出了现有模型中存在的不足和急待解决的几个关键问题。 1 3 国内的发展和研究概况 国内在固体火箭发动机两相流数值模拟方面的研究起步比较早，并且取 得了一定的研究成果。方丁酉( 1 9 8 8 ) | 2 2 - ”1 对固体火箭发动机聩管两相流的研 究作了大量的工作。周旭( 1 9 9 4 ) 采用欧拉一拉格朗日法计算了燃烧室加质 湍流两相流动，气相采用非交错网格的s i m p l e c 求解二维轴对称不可压粘性 控制方程，湍流模型采用k e 模型，颗粒相采用p i s c 法在拉格朝日坐标下 进行求解，气相和颗粒相进行耦合计算。侯晓嗍( 1 9 9 0 ) 用近似因子法和轨 道模型计算了喷管跨声速两相无粘和粘性流动。曾卓雄伫6 】( 2 0 0 2 ) 用t v d 格 式和双流体模型计算了可压稀相两相流场。严红采用s i m p l e 算法结合分散颗 粒轨道模型，在同位网格上实现了喷管中跨声速两相湍流流场的数值模拟。 6 啥尔滨工程大学硕士学位论文 淡林鹏等吲( 2 0 0 2 ) 采用颗粒轨道模型和有限容积的j a m e s o n 格式计算了长 尾喷管纯气相和两相流流场。 虽然已经有很多研究人员开展了两相流动的研究，但对于较大颗粒在喷 管中的运动规律及运动现象的研究这方面的文献或资料较少。 1 4 本论文主要的研究内容 本论文主要用f o r t r a n 9 0 语言编制计算机程序，为了捕捉到流动中的一些 细节问题采用了自适应非结构化网格，用p o r k 算法求解气相方程组，用拉格 朗日方法求解颗粒相方程组，对小扩张比喷管内的两相流动加以研究，同时也 用商业软件s t a r - c d ，采用标准k 一占模型和颗粒轨道模型中的确定轨道模型， 对同一喷管进行了湍流两相流动的数值模拟。通过理论分析和数值模拟，本 研究课题希望收到以下成果： 1 研究喷管内两相流动中气粒的运动规律，并与已有文献进行比较。验证所 用算法的合理性。 2 研究两相流动中大颗粒的运动规律及出现的物理现象，分析原因得出结 论。重点分析喷喉下游轴线附近流场中的气粒两相流动特性及规律，包括弥 散激波产生的可能性，影响的因素( 粒子负载比，颗粒尺寸等因素) ，区域尺 度，及其存在对喷管下游两相流动特性的影响。 哈尔滨工程大学硕士学位论文 i 第2 章喷管内气粒两相流场的数值模拟 2 1 引言 喷管内的气相流场可以同时兼有亚音速、跨音速和超音速，对应的控制 方程各不相同，因而是十分复杂的，再加入粒子形成两相流流场，使得喷管 内的流场更为复杂，但计算机技术和数值计算方法的飞速发展，使得对其内 流场的数值模拟成为可能。本章对喷管内两相无粘流场进行数值模拟。在这 里先扼要的介绍一下流场模拟问题数值计算过程的几个主要环节，本章将按 照这几个步骤进行喷管内气固两相流流场的数值模拟计算。 流场的数值模拟计算主要包括以下几步： 建立物理与数学模型。首先对所研究的实际问题作出一定简化假设，以确 立其物理模型。饲如，当物理过程中流体的物性变化不大，可作常物性的假定； 物理量在某一方向上的变化相对于其它两个方向很小时可以做二维的假设， 等等。根据所确定的物理模型写出该过程的控制方程及相应的初始条件和边 界条件，这也就是确立数学模型的过程。 划分网格。数值模拟计算中用离散的网格来代替原物理问题中的连续空 间，网格中的节点则是所求解物理量的几何位置。从网格的构造来说，可以 分为结构化网格、块结构化网格及非结构化网格三种。这三种网格各有各的 特点，也就各自适应不同的物理问题，具体的应用将在后文中详细说明。 确定控制方程的离散方法。选择多种离散方法中最适用于所模拟物理问题 的控制方程的一种方法，对控制方程进行离散化处理，得到节点间物理量的 代数方程。在离散方程的过程中，基于一定的物理真实性，需要对所求解的 变量在两个节点之间的变化特性作假设，从而选择不同的离散格式。 求解控制方程离散后形成的代数方程组，然后对数值计算的结果从物理过 程的角度进行解的分析。 2 2 控制方程及计算模型 2 2 1 计算模型 由于本章计算的颗粒直径大于分子平均自由程，不适宜把颗粒相当作“拟 流体”来看待，故采用了颗粒轨道模型( 欧拉拉格朗日模型) ，对颗粒在二 哈尔滨工程大学硕士学位论文 维小扩张比喷管内的流动加以模拟。这种模型仍然把气相当作连续相，但把 颗粒相看成是不连续的离散相，对每个粒子( 或粒子群) 在拉格朗日坐标下 进行跟踪。这样，气相采用欧拉型方程，而颗粒相采用拉格朗日型方程。 本章在计算单相流场时，假定气体为忽略质量力的、具有常比热系数的 完全气体，分为考虑粘性和不考虑粘性两种情况进行计算，在对二维小扩张 比喷管内的两相流动编制f o r t r a n 程序，进行流场模拟时，所作简化及假设 的变化如下： 1 气相成分冻结，除了和颗粒接触处外可不计粘性； 2 颗粒为性质均匀的不可压球形，不考虑颗粒容积分数的影响； 3 不考虑颗粒间的相互作用，认为颗粒是离散的； 4 两相之间的相互作用只考虑颗粒与流体之间的动量传递和能量传递； 5 两相混合物与壁之间没有质量交换、动量交换和能量交换。 2 2 2 基本控制方程 1 气相控制方程 本章在进行两相流动计算前，先进行了单相( 气相) 流场的有粘和无粘 计算，有粘时对应求解层流n a v i e r - - s t o k e s 方程，无粘时对应求解e u l e r 方程。 由于n a v i e r - - s t o k e s 方程中不考虑粘性项即为e u l e r 方程，因此本章虽然在两 相流计算时对气相采用鹊是e u l e r 方程，但具体描述方程组时，包括后面要 讲到的方程的离散，还是以n a v i e r - - s t o k e s 方程为主，其表达形式为： 罢+ 掣+ 掣：至婴+ 旦掣坝( 2 - - 1 ) 式中 9 叫咖阪 l 朋+ 甲j p 妒 卜 斗 2 v e 钆虻咖 uf d 粤肌她 ，。_ i l u 哈尔滨工程大学硕士学位论文 并有 其中 g 3 ( u ) =，g 7 ) = j m y i j 篡a v 蛩2 一f j , , ( p e + p ) y 一甜，吆一岛一岛j = ( z + 2 考+ a 罢 j = 卜二维流动。 j = 1 轴对称流动： e _ _ p 1 + 妊2 + v 2 ) ， 、 p 密度 甜x 方向速度分量 v y 方向速度分量 p 压力 e 内能 k 热传导系数 r 温度 1 0 ( 2 2 ) ( 2 3 ) ( z 一4 ) ( 2 5 ) 坚砂 乐+ o 坚 毽 露卜 勺 + o 屯 塑砂 ， 跏一知 钆一缸 m 耖砌学 印 k q 哈尔滨工程大学硕士学位论文 “动力粘性系数 五第二粘性系数 e ，颗粒和气体间x 方向动量源项 只，颗粒和气体间y 方向动量源项 q p 颗粒和气体间的能量源项 以上参数中的粘性系数是随温度t 的变化而变化的，其取值利用工 程上常用的苏士兰( s u t h e r l a n d ) 渊公式来得到： 趔：f 三r 丝 地l 瓦jr + t 式中；参考温度瓦= 2 7 3 1 6 k ；正为s u t h e r l a n d 常数， 正= 1 1 3 k ；肌为一个大气压、2 7 3 1 6 k 下的空气粘性系数， 1 7 1 6 1 0 - 5 n s m 2 。 第二粘性系数五为： 兄= 一詈 传热系数k 可以根据普朗特数p r 得到，在本文中p r = 0 7 。 p r ：h c p 七 2 无因次化 用激波前的状态 p ，= p l ，只= 只，= 4 p , , o , ， x + = 兰，y + = 兰，广= 圭蚱三 7 三二。 p ：卫，户：卫，r ：乓，f ：垒 p，pfpp ， = 尝，拈一罢3 ， 其中l 参考长度， 上标乖一无因次变量， 下标1 激波前的变量 ( 2 6 ) 对空气取 其值取作 ( 2 7 ) ( 2 8 ) 哈尔滨i - 程大学硕士学位论文 下标r 参考变量 对上述方程进行无因次化，有 了o u * + 三验+ 掣：圭篱骝+ 掣m ( u ) c z 圳 矿。苏+ 却+ r e 础+ 砂+ 。7 式中u = g 。( u + ) = g ( 矿) = s ( u + ) = 艮岵各等 r 【“+ + v 矗+ 万y 耳a * 歹o t *，一1j 1 叫 f ( 扩) = 一j v ( p e + p l y 一群；蜀乏一q ：j 每一q ； p ， p vu p + v 2 + p + p + v + e + + v + p + ( 2 一1 0 ) p z p 十 +。 一 。v f 甜 “ 甜 “ d d p p ，。l i | 、i， u ，i f ， “。v f 加一咖 1 1 f “ 一 甜一：毽 矿+， 嚏够 + = k 哈尔滨工程大学硕士学位论文 川i + 2 + ) 蒡材缸a u _ l 。 ( 2 - - 1 2 ) ， p + = ( r - 0 p 【e 一丢( “2 + v 坨) 】( 2 - - 1 3 ) r e ：= p ，- u , l 社， ( 2 一1 4 ) 本文中所使用的公式均为无因次后的公式，为方便起见，上标在后面 的章节中将不再写出。 3 颗粒相控制方程 颗粒相的控制方程为： 鲁= 警= 咋 堕=一3鬲pdt 4 ( 舭魄f口。觑“、 “o 盟=三矗啪训hidt 4n n “ 亟d t 嘶旦c 去坠p r 叮吲 。d ：d 。 、 ” 式中下标k 第k 束粒子 p 。第七束粒子的直径 k 一 第七束粒子x 方向的滑移速度 f v 一心l 第七柬粒子y 方向的滑移速度 c & 第七柬粒子的滞止系数 c 粒子比容 n u 。第k 柬粒子的努谢尔特数 ( 2 1 5 ) ( 2 一1 6 ) ( 2 1 7 ) ( 2 1 8 ) 粒子的滞止系数可以由标准滞止函数计算得到： ： 蠢( 1 0 + l 6 r e ) 瞅1 0 0 0 ( 2 州) 1 0 4 4r e 1 0 0 0 粒子的努谢尔特数可以由下式得到： n u 。：2 + 0 6 r e ；k ：p ( 2 - - 2 0 ) 粒子的雷诺数和马赫数分别为 瓯：剑翻 弘 m a k ：豳 ( 2 2 1 ) ( 2 2 2 ) 式中 i u - - i 第七束粒子的滑移速度 4 气相与颗粒相之间的耦合 由于本章中，气粒两相间的相互作用只考虑动量传递和能量传递，因此， 颗粒与气相耦合的源项表达形式可写为： x 方向动量源项： 矗。南p ；峨帅训卜u , l c o t ( 2 - - 2 3 ) y 方向动量源项： 2 竞p ；【研m ( p ) l p 叱i 】( 2 - - 2 4 ) 能量源项： q ，= 专鲁；脚t 札即删】( 2 - - 2 5 ) 这里我们假定气相速度和温度在所计算的控制容积内为常数。 上式中 ，所计算的控制容积 哈尔滨工程大学硕士学位论文 l 第k 束粒子中的粒子数 m 与气体速度、喷入时间步长和所计算的网格面积有关，可由下式得到： m = i 6 器瓯( 2 - - 2 6 ) 其中玎粒子装填比 p + 气体的滞止密度 甜+ 气体的滞止速度 厶喷入时间步长 珂喷入的粒子束总数 0 取决于粒子大小和装填比，它可以比计算时间步长垃大一些。并且在每 - - a t , 上喷入的粒子束总数和在网格上粒子束的数目是一定的，但要是碰到 粒子束内的粒子不到个时，则要适当增大粒子喷入的时间步长。 具体的实现过程为：先计算稳态的气相场，然后在此基础上加入颗粒， 得到稳定的流场后，根据颗粒所在的位置寻找最近的气相网格点，利用线性 平均求得对应颗粒所处位置的气相参数值，再对气相进行计算。即对气相和 颗粒相之间采用的是双向耦合，既考虑气相流场对颗粒相流场的影响，又考 虑颗粒相流场对气相流场的影响。可用下式表示 矿1 = 厶学l , ( m ) l c 学q 1 ( 2 - - 2 7 ) 式中f 时间步长 厶气相的计算过程 l 颗粒相的计算过程 啥尔滨工程大学硕士学位论文 2 3 网格划分 网格划分即网格生成技术是计算流体力学( c f d ) 发展的一个重要分支， 它是随着数值计算方法和计算机技术的发展而发展起来的，在数值计算中占 有重要地位。实际上，流动问题数值计算结果的最终精度及计算过程的效率， 主要取决于所生成的网格与所采用的算法，可见网格生成是计算流体力学作 为气体动力学工程应用的有效工具所面临的关键技术之一。 从网格的构造来说，可以分为结构化( s t r u c t u r e d ) 网格、块结构化 ( b i o c k - s t r u c t u r e d ) 网格、以及非结构化( u n s t r u c t u r e d ) 网格三种。结构化 网格的节点以阵列形式排列，与计算机语言自然匹百已，便于矩阵演算与操作a 而对于块结构化网格，严格的说，其实块结构化网格是结构化网格的一种类 型，是将结构化网格应用于复杂外形的一种有效方法。当计算区域比较复杂 时，单是用一种规格的结构网格难以妥善的处理所求解的不规则区域，这时 可将整个求解区域按其形状分为若干个子区域，在每个子区域中独立的生成 适合其外形的结构化网格，然后将各个子区域网格连接起来。按其在连接处 方式的不同大致可分为两种类型：拼片式网格( p a t c h e d 或a b u t t i n g ) 例( 或 称对接式) 和搭接式网格( o v e r l a p p i n gg d d ) ( 或称重叠式) ，前者在块与块 的交界处是不重叠的，而后者则有部分区域重叠。在结构化网格中，每一个 节点及控制容积的信息必须加以储存，但该节点的邻点关系则是可以依据网 格编号的规律而自动得出，因而不必专门储存这一类信息，这正是结构化网 格的一大优点。一般来说，结构化网格在二维问题中为四边形网格，在三维 问题中为六面体网格。 与结构化网格和块结构化网格不同，非结构化网格是另一类型的计算网 格。它对网格的节点没有结构性的限制，它的节点和单元分布是任意的，因 此这种网格具有优越的几何灵活性，对各类复杂的边界都能很好的完成网格 划分，生成的网格质量较高，其随机的数据结构有利子进行网格自适应，所 以能够较好的处理复杂边界问题，本文拟采用这种对不规则区域具有特别适 应性的非结构化网格。 2 3 1 非结构化网格的生成 非结构化网格是处理复杂计算域的一种有效方法，其基本思想给予如下 假设：四面体是三维空间最简单的形状，任何空间区域都可以被四面体单元 所填满，即任何空间都可以被以四面体为单元的网格所划分。所谓“非结构 1 6 哈尔滨工程大学硕士学位论文 化网格”就是在这种网格系统中节点的编号命名无一定规则，甚至是完全随 意的，而且每一个节点的邻点个数也不是固定不变的。由于非结构化网格舍 去了网格节点的结构性限制，易于控制网格单元的大小、形状及网格点的位 置，因此与节点排列有序、每个节点与邻点的关系固定不变、结构严密的结 构化网格相比，非结构化网格表现出一种不规则、无固定结构的特点，这种 特点使得非结构化网格具有更大的灵活性，对复杂外形的适应能力非常强。 对于结构化网格，在计算域内网格线和平面都应保持连续，并正交于物体边 界和相邻的网格线和面：而非结构化网格则无此限制。另外，由于在非结构 化网格系统中，一个节点周围的节点数和单元数都不是固定的，可以根据流 场性质方便的做自适应调整，从而合理地安排网格节点的疏密，以便提高计 算精度，如流场中物理量变化剧烈处。网格节点可安排较密，物理量变化平 缓的地方网格点可安排较稀疏：再如流场中相对小尺寸处，网格可划分较密， 尺寸相对大的地方网格划分可以稀疏一些。从网格本身及与此相关的算法来 看，非结构化网格比结构化网格更接近计算几何学和计算机科学的原则。正 是因为有着一系列的优点，伴随着计算机科学的飞速发展，非结构化网格技 术在二十世纪8 0 年代末和9 0 年代初得到了迅速发展和广泛应用。 生成非结构化网格的方法主要有两大类：前沿推进法( a d v a n c i n gf r o n t m e t h o d ) ；d e l a u n a y 三角形化方法( d e l a u n a yt r i a n g u l a t i o n ) 。 前沿推进法是网格和节点同时生成的非结构化网格生成法，它从边界上 的网格点所形成的一系列线段( 前沿) 出发，逐一与区域内部的点形成三角 元，所形成的三角元的新的边加入到前沿的行列，而生成该三角元的原来的 线段( 前沿) 则从前沿行列中自动消去，如此不断地向区域内部推进，最终 使得前沿的行列为空，而所生成的三角元覆盖整个计算域。 生成非结构化网格的另一种重要的方法是d e l a u n a y 三角形化方法，这种 方法基于计算几何学的严密规则之上，在2 0 世纪的8 0 - 9 0 年代得到了迅速的 发展。此三角形化方法源于d i r i c h l e t 在1 9 世纪5 0 年代提出的一种利用已 知点集将平面划分成凸多边形的理论。用d e l a u n a y 三角形化方法生成的网格 有可能某些网格单元位于区域外或者与边界相交，使流场边界的完整性遭到 破坏，但用这种方法生成的网格单元质量高，而且网格生成效率也相对较高， 所以，d e l a u n a y 三角形化方法目前已经成为生成非结构化网格的种主要方 法洲，本文选用该种网格生成方法。 哈尔滨工程大学硕士学位论文 一、边界节点的生成 为了在任意个计算域上生成二维非结构化网格，首先要做的，就是确 定该计算域上边界节点的空间分布，这样才能既满足所需求解的解的精度要 求，又可以准确表示出实际计算边界的形状。然后再根据边界的形状，创建 计算域内部的网格节点。这样的话可以由边界节点的疏密程度使内部节点规 则地分布在计算域内。本文进一步发展了一种简单的算法( l o1 9 8 5 ) 来生成 内部节点。整个过程描述如下： 1 首先确定所要计算的区域在x 方向上和y 方向上的最大值和最小值( 即 x 。，x 。，y m ，y 。) 。 2 选取合适的长度d ，用来表示节点间的横向平均距离，并且所选的d 值 能准确的表示出所要计算区域的实际边界，同时也满足解的精度要求。 3 使参考点尽量呈等边三角形形状排列覆盖整个计算区域。即每行两个相 邻节点之间距离为d ，行与行之间的距离为车d 。 4 计算域的边界由彼此不连通的封闭环组成。外部边界上的节点按逆时针方 向排列，而内部边界上的节点则按顺时针方向排列，如图2 1 所示，封闭环 s t = 只q ，i = l ，2 ，n b ) 表示外部边界，封闭环s 产 只q f ，i = n l + l ，n 1 + 2 ，n b ) 表示内部边界，8 表示边界点的数目。在封闭环上，线段；酒两个端点的 坐标可以写为( 算( 霉) ，) ，( 只) ) 和( x ( q ，) ，y ( q ) ) ，只、q 之间的关系为： 墨= q w t ，罡= 9 ，曰= q 一。，昂一= q “ y 嘣 ，2 孤 x m m y “o 图2 1 计算域内节点生成示意图 1 8 5 在边界线上需要的地方插入新的边界节点。为此，首先考虑水平边界。如 果满足 m = y 2 ，r x 2 一t 卜1 4 d ( 2 - - 2 8 ) 则需要插入新的节点。对于其它非水平的边界，设想在计算域内存在一些间 隔( 峨) 均匀的水平线，这些水平线把边界线分割成按升序 z = 【石，筋，筋- p z 2 _ 1 排列的均匀的点集，当满足 ( u h 、fv一hk) c ：拭2 ，胎 q 1 式中 ( x ) ，y ) ) ，i = 1 ，2 ，肋边界节点 c 与边界平均节点间距相关的常数，常取为0 7 d 则需要插入新的节点。 。可以由下式计算得到 哪样( 2 - - 3 0 ) 6 在所有的新边界节点都插入后，内部节点的插入只要不断重复步骤5 就可 以实现。设想的每一条水平线上都有一组相交点z 。在该水平线上的参考点 如果在计算域内部，并且满足 甓阮，鼽1 1 f = l ，3 ，2 n i ( 2 - - 3 1 ) 就是网格节点。如果不满足式( 2 - - 3 1 ) ，丽且这些参考点与边界节点相距较 远，且有 ( 。，一x ( 暑) ) 2 + ( 爿：一y ( 日) ) 2 c 2 ，f = 1 , 2 ，n b ( 2 3 2 ) ( k 一石) cf = 1 ，2 ，2 n( 2 3 3 ) 那么，就需要在每条水平线上用步骤5 中的算法，再创建新的内部节点，此 时，有 哈尔滨工程大学硕士学位论文 y = h k ，七= l ，2 ，m ( 2 - 3 4 ) ( 2 3 5 ) ，设想的水平线数 线与线之间的竖直间距 最终可以得到整个计算域的节点( 薯，乃) ，f - 1 ，2 ，n b ，口+ 1 ，p ， 其中，n b 表示边界节点总数，n p 表示内部节点总数。 二、d e l a u n a y 三角形化方法 基于d e l a u n a y 三角形的网格生成方法，利用一个特别简单的几何准则， 获得点的连接而得到合适的不相交的单元。这个几何准则提供了一个连接点 的机制。网格点的生成必须是独立考虑，所以用d e l a u n a y 三角形方法生成网 格包含了两个方面；点的生成和网格点的连接。 d i r i c h l e t 提出了一种方法，对于任何空间中给定的一个区域，可以系 统地分解成一个封闭的凸区域的集合。对( 只) 这个点的集合，空间利用下面 的准则被分割成区域( 巧) 的集合：在区域巧中任意一点到只点的距离比到其 它点的距离更近。这样的间隔称为d i r i c h l e t 分布。一个闭合区域的分布产 生了一个不重叠的凸区域的集合，称为v o r o n o i 区域，它覆盖了整个计算区 域。 从理论上讲，如果点的集合由( 只) 表示，那么v o r o n o i 区域的定义如下： ( k ) = p ：l i p p , l l i i p - 弓l l ，w f ( 2 3 6 ) 也就是说v o r o n o i 区域( k ) 是点p 的一个集合，在此集合中点p 和只之间的 距离比点p 和其它点的距离更近。所有满足这一条件的点p 构成v o r o n o i 区 域。从这个定义可以看出，在二维空间中，构成一个v o r o n o i 多边形的某条 边一定是相邻两点之间连线的垂直平分线。如果用直线连接具有共同边界的 2 0 d 打丁 =胡 中式 哈尔滨工程大学硕士学位论文 一对点，那么就形成了由点集( p ) 组成的三角形。用这种方式形成的三角形 称为d e l a u n a y 三角形。图2 2 给出了这种几何结构的一个例子。 图2 2d e l a u n a y 三角形( 细线) 和v o r o n o i 区域( 粗线) d e l a u n a y 三角形化方法具体步骤如图2 3 所示( 以5 点为例) ： 1 首先构造一个“超级三角形”，这个“超级三角形”包含了所有要插入的 网格点。同时插入网格点l 构成三个网格单元，如图2 3 a 所示； 2 假设已经插入了3 个网格点，现在是插入网格点4 ，如图2 3 b 所示； 3 在所有已生成的d e l a u n a y 三角形中判断它们的外接圆是否包含此点，如 包含此点，那么此单元就归到被打破的序列； 4 把被打破序列中的单元删除，形成一个称为d e l a u n a y 空腔的多边形，如 图2 3 c 中虚线所示： 5 连接插入点与多边形的每一条边，重新形成d e l a u n a y 三角形，如图2 3 d 所示； 6 对网格点5 重复同样的过程，得到图2 3 e 的结果，然后删除所有与最初 “超级三角形”顶点连接的三角形单元，形成了最终单元，如图2 3 f 所 示。 2 1 哈尔滨工程大学硕士学位论文 公公叁 煮丛 ( d )绀 图2 3d e l a u n a y 方法形成三角形过程 由于d e l a u n a y 三角化过程中，只是考虑了所有网格点的坐标。而对于 网格点的具体性质不加考虑，这样在最终生成的网格中，对于凹边界的流场 区域，可能有一些单元会出现在流场外的非结构化网格区域中，这也是 d e l a u n a y 三角化方法的一个缺点，必须要具体处理，保证物理边界和网格边 界的一致性，即物面完整性。为此，可采用根据边界网格点的具体性质，并 利用构成d e l a u n a y 三角形的y o r o n o i 顶点的结构。在构成具体d e l a u n a y 三 角形时，必须要保证三角形的面积为正，即三角形三个顶点的排列顺序是逆 时针方向。其次在离散边界网格点时，利用阵面推进法中的思想，规定边界 线段的方向，要始终保持物理区域在边界线的左边，这样外边界为逆时针方 向，而内边界为顺时针方向。在离散成网格点时，按照边界的方向进行。保 证网格点的序号沿着边界不断增大。对由边界点组成的初始网格单元的三个 顶点，判断它们网格点序列大小关系，从而决定它们的保留与否，一般情况 下若按三角形的顶点排列方向，三个顶点的序列号是增加的话，那么这个单 元就是要保留的单元。 本章采用的是自动插点的d e l a u n a y 方法，它根据边界点的分布情况具 体决定在网格化区域中如何插入网格点，而不需要人为的干预。即直接在三 角形的内部插入网格点，位置在三角形的质心，这样有利于提高网格生成效 率。 哈尔滨工程大学硕士学位论文 2 3 2 网格的优化 网格的质量与计算模拟的结果之间有着十分密切的关系，因此对于非结 构化网格的优化是必不可少的内容。优化的基本目的是提高网格的质量，在 二维网格优化方法中，对利用d e l a u n a y 三角化方法生成的网格基于弹性原理 的方法加以优化，其主要思想就是使得内流场中的网格点移向包围它的多边 形的形心，如图2 4 所示。这样的话，每个三角形能更接近于等边三角形， 有利于提高计算效率和计算精度。移动网格点的新旧坐标间的关系如下公式： 乙【墨，咒j ( 训) = 趔；矿江胭“，n b + 2 ，胛 q 7 式中( 一，m ) 内部节点i 的坐标； n ( i ) 与i 点相连的点数； v t ( i ，j ) i 点周围的临近网格点 b i b 边界节点总数； - 俨内部节点总数。 图2 4 网格的优化 2 3 3 网格自适应方法 在数值计算过程中，流场的变化