喷管两相流耦合数值模拟与粒子热增量预测

苏鹏辉 艾邦成 潘宏禄

(中国航天空气动力技术研究院， 北京100074)

现代固体火箭发动机为了抑制不稳定燃烧，提高工作效能，常使用含铝粉的复合推进剂，其燃烧产物中含有氧化铝凝聚相粒子，质量分数有时可高达30%，形成典型两相流动［1］.两相流中粒子对喷管热环境的影响主要体现在两个方面:一方面粒子在流场中存在速度和温度滞后，通过与气相间的动量和能量交换引起气相流场变化;另一方面高速运动的高温粒子撞击喷管壁面产生壁面热流增加和粒子冲刷烧蚀等现象，对喷管喉衬热防护提出了更高要求.因此基于燃气-粒子相互作用的喷管两相流研究已成为喷管喉衬热防护研究的重要方面.

喷管两相流研究主要有数值模拟和试验测量这两种手段.在数值模拟方面，目前国内外学者主要基于欧拉-拉格朗日描述的粒子轨道模型对不同状态下喷管两相流的流场参数和粒子运动状态进行数值计算分析.文献［2－4］利用颗粒轨道模型对喷管两相流进行了数值模拟研究，得到了丰富的成果，显示了颗粒轨道模型的应用价值.

在喷管两相流试验，特别是喷管粒子热增量试验方面，文献［5］进行了不同攻角平板的粒子热增量试验，粒子垂直平板速度所对应动能的70%转化为平板表面热能.文献 ［6］进行了超声速球头绕流粒子热增量试验，在球头驻点处粒子动能70%转化为球头驻点热能.

目前已有的喷管两相流数值模拟集中关注于粒子运动以及流场参数变化，而粒子热增量研究则主要集中于试验方面.通过喷管内流场两相流数值计算获得壁面粒子热增量分布并不多见.本文基于颗粒轨道模型，以含有氧化铝粒子的发动机喷管两相流为研究对象，通过有限体积方法和单元内颗粒源 (PSIC，Particle Source in Cell)两相耦合方法，模拟了不同直径粒子对喷管流场及喷管壁面热环境的影响.

1 控制方程及边界条件

1.1 气相控制方程

两相流中气相控制方程组采用守恒形式完全气体三维N-S方程组

其中，F，G和H为对流项通量;Fv，Gv和Hv为粘性及传热耗散项通量.为封闭控制方程组需补充理想气体状态方程p=ρRT，气体动力粘性系数μ通过萨特兰公式得到.右端源项S表示粒子对气相作用，其表达式在下文给出.

1.2 粒子相控制方程

氧化铝粒子的熔点2345K，沸点3250K［7］，当燃烧室温度在3000～3100K时，燃烧室内及喷管喉部前端的氧化铝粒子在表面张力作用下形成球形液滴［8］.使用粒子轨道模型对粒子相建模，粒子相与气体相方程独立求解，利用源项耦合，通过跟踪计算大量粒子的实际运动轨迹及粒子状态随轨迹变化，实现对两相流物理过程的模拟［9］.定常流动中，每个粒子轨道均代表了在流场入口截面上这个粒子所占据的面积微元内所有流经此微元的粒子轨迹及状态变化.此方法节省存储空间，易于并行计算，易于模拟粒子与气相之间的速度和温度滑移及粒子相复杂变化过程.

针对固体火箭发动机喷管两相流特点，对粒子相做出假设:①粒子均为一定直径的球形;②粒子与气相之间无质量交换;③粒子与气相之间存在对流换热;④粒子与气相之间只考虑阻力;⑤粒子与粒子之间没有相互作用;⑥粒子不考虑湍流扩散效应.粒子速度、位移及温度随时间变化的拉格朗日型控制方程组如下:

其中，up，xp，Tp，Ap，mp和 cp分别为粒子速度、位置、温度、表面积、质量和热容;u和T为流体速度和温度，其余系数计算可参考文献［10］.

对有限体积法中每个网格单元，统计穿过此单元的全部粒子在单元内交换的动量和热量，除以该控制体的体积Ω，即得到此控制体内粒子源项［11］，S的表达式如下:

2 数值方法及定解条件

2.1 数值方法

气相流场控制方程采用有限体积方法求解，其中对流项通量使用Roe格式求解，耗散项通量采用高斯方法求解，时间离散采用LU-SGS(Lower-Upper Symmeffic Gauss-Seidel)隐式方法.粒子相控制方程采用欧拉显式格式求解.

两相耦合过程采用PSIC方法求解［12－13］，过程是:首先获取气相场定常解，之后求解粒子轨道方程，统计粒子与气相之间的动量及能量交换，作为源相加入气相流场，重新计算气相流场定常解，反复迭代，直至收敛到预期精度.迭代过程中为减小振荡，加快收敛速度，在更新粒子源项时对粒子源项进行如下低松弛处理:

其中，Sn为第n次迭代中气相控制方程的粒子源项;S*为求解粒子轨道方程得到的新源项;Sn+1为第n+1次迭代中气相控制方程的粒子源项;α为松弛因子，在计算初期取0.5左右.

2.2 定解条件

2.2.1 粒子相入口条件

粒子在喷管入口加入流场，初始速度和初始温度与当地气相参数相同，粒子质量流量在入口截面上满足均匀分布，每条粒子轨道所代表粒子质量流量为mi=m/Nk，其中m为入口截面上总的粒子质量流量，Nk为粒子轨道总数.

2.2.2 粒子相壁面条件

有关试验研究［5－6］表明，粒子撞击壁面时，粒子垂直壁面速度所对应动能的70%转化为壁面热能，设粒子撞击前平行和垂直壁面的速度分别为Vp1、Vn1，撞击之后反弹速度分别为 Vp2、Vn2，不计粒子内能变化，根据能量守恒有:

粒子撞击壁面时，粒子携带的能量部分转化为壁面热量，根据上述试验结论，每条粒子轨道在与壁面碰撞时，产生的壁面热增量为

对全部颗粒轨道的壁面热增量进行累加，可得到总的壁面热增量.

3 结果与分析

选取典型喷管外形和燃烧室工况进行喷管两相流的计算.图1是喷管三维网格示意图，左侧为入口截面，右侧为出口截面.计算所用参数:燃烧室总压5MPa，总温3000 K，喷管壁面为等温壁，壁温300K.粒子由入口处加入流场，粒子质量流量253kg/(m2·s)，粒子热容为cp=880 J/(kg·K)，粒子物质密度为ρp=3960kg/m3.

图1 喷管计算网格示意图

3.1 粒子轨道数量对热增量均匀性影响

粒子轨道模型中需要模拟的轨道数量应根据具体问题的特点而确定，为了确定适合本文算例的轨道数量，进行了轨道数量与粒子热增量均匀性的分析.由于该喷管几何形状以及流场边界条件的轴对称性，在理想情况下同一轴向位置处的周向各面粒子热增量应保持一致，但是离散轨道的不均匀性会导致计算结果轴对称性下降.模拟轨道数量越多，计算结果的均匀性越好，但计算所消耗机时也会成比例增加，因此通过对不同轨道数的模拟结果进行热增量均匀性分析，可以确定本问题需要的轨道数量，在数据均匀性和计算时间之间取得平衡.

图2粒子轨道数量N分别为1万，10万，40万和100万时在每个轴向位置x处的热增量计算值沿周向分布的标准差σ.

图2 不同轨道数量热增量分布标准差

从图2可以看到:当N＜40万时，数据均匀性较差，随着N增加，数据均匀性有明显改善;当N＞40万时，数据均匀性随N增加改善较小.因此权衡数据均匀性和所需计算时间，本研究中轨道数量取为40万条.

3.2 粒子轨道分析

为了考察粒子轨道模型对喷管内燃气-粒子两相流问题求解的有效性，以及观察粒子在喷管两相流中的流动特点，对不同直径粒子的轨道进行了计算和分析.粒子直径决定了粒子质量和受力，随着粒子直径增大，粒子在流场中跟随性变差，分布范围逐渐向轴线靠拢，在扩张段壁面附近出现逐渐增大的粒子稀疏区域.

粒子直径较小时惯性也较小，在流场中的跟随性很好，图3a是直径1 μm粒子的轨道分布，可以观察到粒子轨迹在喷管中分布均匀，在喷管收缩段内粒子轨迹随喷管型线转弯，向喷管中心轴线靠拢，而在喷管扩张段内粒子轨迹随喷管型线一起扩展至整个喷管出口截面，粒子在收缩段内与喷管喉衬表面没有明显碰撞，在扩张段内亦无明显的粒子稀疏区域.图3b～图3d是直径为10～40 μm粒子的轨迹分布，可以观察到，随着粒子直径增加，粒子惯性愈加明显，在喷管收缩段内粒子由于无法及时随喷管型线转向而与喷管壁面发生碰撞.在喷管扩张段，粒子轨道集中分布于轴线周围的圆形区域内，在喷管壁面附近产生了粒子稀疏区域.同一轴向位置处，小直径粒子的分布区域半径较大，粒子稀疏区域较小，而大直径粒子的分布区域半径较小，粒子稀疏区域较大，大直径粒子的轨道更加趋向于喷管轴心.由以上结果可知，使用粒子轨道模型可以得到粒子在两相流中的典型分布特征，粒子轨道模型对喷管两相流问题是适用的.

图3 不同直径粒子轨道分布

3.3 两相流流场分析

含有粒子的两相流流场与纯气相流场相比有较大区别，为了研究粒子对流场的影响，需要采用两相耦合方法对流场进行数值计算.

图4给出了粒子直径为20 μm时两相流与纯气相流流场参数分布对比.图4a～图4c中，上半部分为两相流流场参数分布，下半部分为纯气相流场参数分布.可以看到，两相流中，气相速度和温度受粒子影响较大，而压力分布受粒子影响较小，在喷管扩张段内两相流中气相速度较纯气相流动时降低，两相流中气相温度较纯气相流动时升高，这是由于扩张段内气相速度增加温度降低，而粒子速度和温度都存在滞后，因此粒子对气相产生减速加温作用.图4a和图4b中，喷管扩张段轴线附近，与纯气相流动相比，两相流的速度和温度等值线出现朝向喷管出口方向的偏移，这是由于在轴线附近位置粒子轨道分布较为密集，气相受到粒子的减速加温作用更加明显.由图4d～图4f可见，在扩张段内相同轴向位置处，两相流比纯气相流速度更低、温度更高.由上述结果可知，两相耦合计算可以较为全面地反映流场中粒子与气体之间的相互作用情况，所得结果更符合真实情况.

图4 两相流与纯气相流流场参数对比

3.4 粒子相壁面热增量分析

为了预测粒子撞击喷管壁面时动能转化为热能所产生的壁面热增量，在两相流计算中对粒子轨道与壁面碰撞情况进行统计分析，得到了喷管壁面热增量的分布范围和热流大小.

图5给出了对应于不同粒子直径的壁面热增量分布.可以观察到，粒子热增量集中分布于喉部前端收缩段内，收缩段的前后粒子热增量均迅速减小，喷管扩张段内粒子热增量并不明显.粒子热增量的分布范围随粒子直径增大略有增加，而粒子热增量峰值则随粒子直径的增大迅速增加，这是由于小直径粒子随流性很好，不会与壁面发生剧烈撞击，因此其热增量相对较小，而大直径粒子惯性较大，在喷管收缩段内无法及时随流场转向，因此撞击壁面的现象比小直径粒子更加明显，壁面热增量也就随之增加.

图5 不同直径粒子热增量分布

4 结论

本文采用基于粒子轨道模型和有限体积方法的PSIC两相耦合算法，结合基于粒子热增量试验的粒子壁面反弹模型，对含有氧化铝粒子的固体火箭发动机喷管两相流及其对喷管壁面热环境的影响进行了研究，得到以下结论:①使用粒子轨道模型得到的粒子轨道分布和流场参数分布能够反映出喷管两相流区别于纯气相流动的特点，说明粒子轨道模型对喷管内燃气-粒子流动是适用的;②两相流耦合计算所得到的气相温度和速度分布与纯气相流动相比有较大区别，耦合计算所得结果反映了气相和颗粒相之间的相互作用，对实际流动现象模拟更加准确;③粒子动能向热能转化产生的粒子热增量集中分布于喉部前端收缩段内，收缩段前后热增量迅速减小，喷管扩张段内粒子热增量并不明显，热增量峰值随粒子直径增加明显增大.

References)

［1］张宏安，叶定友，侯晓.固体火箭发动机凝聚相微粒分布研究现状［J］.固体火箭技术，2000，23(3):25－28

Zhang Hongan，Ye Dingyou，Hou Xiao.Recent status of research on condensed-phase particles distribution in solid rocket motor［J］.Journal of Solid Rocket Technology，2000，23(3):25－28(in Chinese)

［2］李雅娣，陈林泉，蹇泽群.粒子尺寸分级的喷管两相流场计算［J］.固体火箭技术，2003，26(3):32－34

Li Yadi，Chen Linquan，Jian Zequn.Numerical study of twophase nozzle flow with classified particles［J］.Journal of Solid Rocket Technology，2003，26(3):32－34(in Chinese)

［3］李强，李江，刘佩进.模型发动机内凝相颗粒碰撞的数值模拟［J］.推进技术，2008，29(1):18－21

Li Qiang，Li Jiang，Liu Peijin.Numerical simulation of particle collision and impingement for a model SRM ［J］.Journal of Propulsion Technology，2008，29(1):18－21(in Chinese)

［4］Sabnis J S，Jong F J，Gibeling H J.Calculation of particle trajectories in solid-rocket motors with arbitrary acceleration［J］.Journal of Propulsion and Power，1992，8(5):961 －967

［5］姜贵庆，王淑华.喷管内粒子热增量机理研究［J］.固体火箭技术，1990，13(2):18－22

Jiang Guiqing，Wang Shuhua.The investigation of particle impingement heating mechanisms of SRM ［J］.Journal of Solid Rocket Technology，1990，13(2):18－22(in Chinese)

［6］Dunbar L E，Courtney J F，Mcmillen L D.Heating augmentation in particle erosion environments［R］.AIAA-74-607，1974

［7］Patnaik P.Handbook of inorganic chemicals［M］.New York:McGraw-Hill，2002:11 －12

［8］傅德彬，姜毅，袁增凤.燃气射流气固两相流场研究［J］.北京理工大学学报，2004，24(5):387－391

Fu Debin，Jiang Yi，Yuan Zengfeng.Numerical study on gassolid two phase jet flow ［J］.Journal of Beijing Institute of Technology，2004，24(5):387－391(in Chinese)

［9］刘君，黄琳，夏智勋.潜入式摆动喷管两相内流场数值模拟 ［J］.固体火箭技术，2003，26(4):18－21

Liu Jun，Huang Lin，Xia Zhixun.Numerical simulation on two-phase inner flow field of submerged vectored nozzle ［J］.Journal of Solid Rocket Technology，2003，26(4):18 －21(in Chinese)

［10］方丁酉.两相流动力学［M］.长沙:国防科技大学出版社，1988:46－47

Fang Dingyou.Two-phase flow dynamics［M］.Changsha:National University of Defense Technology Press，1988:46 －47(in Chinese)

［11］尚庆.湍流-颗粒反应相互作用的理论分析 ［D］.北京:清华大学航天航空学院，2005

Shang Qing.Studies and application of theoretical model for turbulence-particle reaction interactions［D］. Beijing:School of Aerospace，Tsinghua University，2005(in Chinese)

［12］傅德彬，姜毅.燃气射流气固两相数值模拟与颗粒冲刷分析 ［J］.弹箭与制导学报，2004，24(1):63－66

Fu Debin，Jiang Yi.Numerical study on gas-solid two-phase jet flow and solid particle erosion［J］.Journal of Projectiles Rockets Missiles and Guidance，2004，24(1):63 － 66(in Chinese)

［13］张旭.固体火箭发动机两相流特性研究 ［D］.哈尔滨:哈尔滨工程大学航天与建筑工程学院，2006

Zhang Xu.The investigation of characteristics of two-phase flow in solid rocket motor［D］.Harbin:College of Aerospace and Civil Engineering，Harbin Engineering University，2006(in Chinese)