当前,高超声速飞行器是各大国积极研究的热点,由于飞行过程中面临长时间严酷的气动加热,所以热防护系统是高超声速飞行器必须解决的关键问题之一.电弧风洞是开展飞行器热防护系统防热材料、结构热考核试验等必不可少的设备,在高超声速飞行器论证研制阶段发挥着不可替代的作用.采用电弧加热试验气体的地面风洞试验设备中由于电极烧蚀会产生微小金属颗粒, 当这些微小颗粒进入流场后会对试验流场产生污染,进而对设备内流场特性、 试验模型以及喷管壁面等产生一定程度的影响.开展气粒两相流动数值模拟研究有助于更好地认识喷管内部流场特性, 指导和分析电弧风洞试验, 配合地面试验开展高超声速飞行器热防护研究.
国内外研究两相流的数值模拟方法主要有Euler-Euler、Euler-Lagrange、Lagrange-Lagrange等模型.Euler-Euler模型将气相和颗粒离散相全部认为是统计连续,具有在计算上相对易于求解和物理理论方面易于理解的优点,但是对于稠密气粒两相流动,因控制方程是双曲型的,有着非守恒项和刚性源项,给数值模拟带来了极大的困难.Euler-Lagrange模型允许直接给定单个颗粒的物理特性,便于跟踪模拟有蒸发、挥发、燃烧等异相反应情况下的颗粒的经历,可以获得不同尺寸颗粒离散相的颗粒轨迹等具体信息,另外颗粒方程是常微分方程,相对容易实现.但是,由于颗粒被完全考虑成离散体系,颗粒运动方程数多,对颗粒的跟踪定位方法复杂,计算量随颗粒轨迹数目增加而增大,并且当颗粒相的质量分数增大时,颗粒对流体的影响也逐渐增大,在交替迭代中这种影响可能导致收敛困难.Lagrange-Lagrange模型对气相和颗粒离散相均采用瞬态模型计算,对剪切流动中的离散相传输给出了物理描述.目前仅限于流场为大尺度涡控制的二维流场,离散相的传输仅由大尺度结构控制,忽略了颗粒之间的相互作用[1].总的来说,目前对气粒两相流的物理机理认识还很有限,这在一定程度上制约了该学科的发展和实际工程应用,也使得气粒两相流在很多研究领域得到了更加广泛地重视[2-5].
研究人员针对发动机喷流气粒两相流的研究较多[5-7],而针对地面风洞流场开展两相流场数值模拟的研究相对较少,且是针对人为制造的两相流场,如采用添加粉尘研究风工程或采用示踪粒子对流场进行测量等产生的风洞内两相流动等[8-9].电弧风洞喷管内部流场是由于电极烧蚀产生的金属颗粒进入试验气体流场导致的两相流,目前此类颗粒的物理性质还不够明确,尚未发现针对该领域的研究工作.微小颗粒的引入对流场特性的改变甚至对试验段模型热力载荷的影响值得研究,这对结合气相流动控制方程开展数值模拟方法中物理化学模型的适用性研究等有重要意义.针对电弧风洞流场颗粒质量分数相对较小的特点,采用Euler-Lagrange模型建立两相流数值模拟方法.通过研究合适的颗粒运动时间步长选取方法,探索高效准确的颗粒跟踪定位算法,运用颗粒轨道模型建立颗粒运动的求解方法.然后基于动量和能量关系建立喷管内部气粒两相流动模拟方法,在耦合计算时需要将颗粒离散成一系列轨迹,根据轨迹上颗粒的特性求解耦合源项.利用建立的两相流耦合模拟方法对典型状态进行模拟,分析两相流动对喷管内部流场的大致影响规律,为深入研究电弧风洞两相流场打下基础.
直角坐标系中,三维黏性可压缩气相流动控制方程如下:
(1)
其中Q为守恒变量,t为时间,E,F,G和Ev,Fv,Gv分别为x,y,z三个方向的无黏通量和黏性通量,右端最后一项S为气相与粒子相之间的耦合源项.
对上述气相流动控制方程采用基于结构网格的有限体积法离散,无黏通量采用改进的Steger-Warming格式计算[10-11],网格界面上的变量通过MUSCL插值重构以提高精度,黏性通量采用中心格式计算.时间推进采用LU-SGS隐式方法,求解过程中通过局部时间步长加速收敛.在喷管入口给定气相总温总压,出口如果是亚声速流动需要给定环境压力,如果是超声速流动则边界参数都由内场外推,壁面采用无滑移边界.
基于喷管流场中颗粒的特点做如下一些假设[12-13]: 1) 颗粒为球形,且内部温度、密度等性质均匀,与气相之间没有质量交换; 2) 颗粒相非常稀薄,粒子之间没有碰撞和热能交换; 3) 作用在气相控制体上的力只有气体压力和颗粒阻力; 4) 颗粒所占有的体积忽略不计,同时不考虑颗粒自身重力的影响.
采用的颗粒相动量和能量方程如下[14-16]:
(2)
(3)
其中V为速度矢量,h为焓值,cp为定压比热容,下标k代表颗粒相,τrk和τhk分别为颗粒相动量松弛时间和能量松弛时间.通常假设微小粒子在喷管内呈惰性,即不考虑质量的变化也不参与反应,流场中每条轨道上有一个粒子束,该粒子束上所有颗粒尺寸和密度相同.气相与颗粒相之间的动量耦合源项和能量耦合源项分别采用下列表达式:
(4)
(5)
式中,Ω为气相流场网格体积,分别为第j条轨道上粒子数目、单个粒子质量、穿越该网格时单个粒子的阻力因子、所受的气动力以及对流传热率.当颗粒在气相流场计算网格内运动时,颗粒相对某网格单元内气相的贡献应当对经过该网格的所有粒子束求和.
在颗粒运动求解过程中需要对颗粒进行定位和跟踪以获取相应位置的气相流场参数值,定位即判断颗粒是否处于某一个网格单元里,而跟踪则是确定一个时间步长后颗粒运动到哪个网格单元.由于在每个时间步长里每个颗粒都要进行定位计算,所以计算效率就显得非常重要.常见的ZL算法采用如下的判定准则,假定一个网格体的网格面外方向为正,如果一个颗粒处于网格体所有面的负方向,则这个颗粒就处于这个网格体内.该算法所需的CPU时间比较少,但需要通过对很多的网格面进行判断才能找到颗粒,并且有可能出现死循环,最终找不到颗粒所在位置[17].根据所要研究的气相流场采用结构网格的特点,考虑网格特征时间步长τG,即颗粒经过网格所需时间,对颗粒运动时间步长Δtk进行了限制,尽可能满足颗粒在一个时间步长后尽量在原来网格位置相邻的网格中.式(6)中α取略大于1的常数,分别为颗粒距离网格体对应方向两侧边界的距离.这样大幅缩小了颗粒定位时的搜索范围,且能有效避免死循环.
(6)
Δtk=min(τrk,τhk,τG).
(7)
求解时在喷管入口给定颗粒相初始粒径分布、初始速度分布、初始温度、初始质量载荷.到达壁面与壁面碰撞后速度将会发生改变,按镜面反射处理.颗粒相到达对称面边界时会穿过对称面,但在对称方向会有颗粒入射进来,处理方法与镜面反射类似.在处理反射边界时,由颗粒所在空间位置以及入射方向(lr,mr,nr)和边界法线方向(ln,mn,nn)求得入射线与反射线所在平面方程(8),再通过方程(9)求得颗粒入射角φ.然后通过方程组(10)即可求解得到反射方向向量(lx,mx,nx),由颗粒运动速度和运动时间结合反射方向求解获得粒子下一时刻的空间位置:
Ax+By+Cz=0,
(8)
(9)
(10)
基于Euler-Lagrange模型的两相流耦合算法,多采用“计算单元内颗粒源项算法”,即PSIC(particle source in cell)算法[13].颗粒相与气相间的动量和能量交换在单元体内进行,具体如图1所示.
图1 相间的动量能量交换示意图
Fig. 1 Schematic of interphase momentum and energy exchange
PSIC耦合算法计算流程如下:
1) 首先求解没有颗粒的纯气相流场,得到气相流场的速度、密度、温度等参数分布;
2) 在已有流场信息的基础上求解颗粒运动方程和能量方程,直到颗粒运动出计算区域,得到颗粒运动轨迹、速度以及温度等参数;
3) 计算相间耦合源项,将源项代入气相控制方程进行求解,回到步骤2),直到满足收敛条件.
该算例是常见的用于验证两相流数值模拟方法的轴对称喷管外形[5,18],喷管尺寸如图2所示,喷管入口气相和颗粒相温度都为555.6 K,气相总压为1.034 MPa.气相为空气,颗粒为三氧化二铝(Al2O3)颗粒,密度为4 004 kg/m3,比热容为1 380 J/(kg·K),选取的颗粒直径分别为1,2,10,20 μm,入口粒子束数目为10,颗粒质量分数为0.3.
图3比较了喷管轴线和壁面边界上的压力分布,试验数据来自文献[5],其中相对光滑的曲线对应喷管轴线上的参数变化,而出现明显拐弯的曲线对应喷管壁面上的参数变化,横坐标0点为喉道位置,参考长度为喉道半径,参考压力为总压.可以看出,纯气相流场计算结果与文献[5]中的试验值吻合较好,验证了计算方法模拟喷管内单相流场的可靠性.
图4比较了喷管对称面上不同直径颗粒的运动轨迹,可以看出该状态下气粒两相流动在喷管喉道下游扩张段靠近壁面处出现了面积不等的无粒子区,且颗粒直径越小无粒子区越小,规律和文献[5]结果一致.通常采用量化颗粒相动量松弛时间和流场特征时间比值的Stokes数描述颗粒跟随性,Stokes数越小,颗粒跟随性越好.该算例中颗粒密度相同,流场特征时间相同,所以Stokes数随粒子尺寸减小而减小.对直径小的颗粒,其Stokes数小,更容易追随气体改变运动速度和充满流场空间,故无粒子区小;而直径大的颗粒Stokes数大,不容易追随气体改变运动速度和充满流场空间,故无粒子区大.此外尺寸较大的粒子与喷管喉道上游壁面碰撞反弹后远离壁面也是无粒子区变大的原因.
图2 轴对称喷管几何尺寸 图3 纯气相流场边界压力比较
Fig. 2 Geometric dimensions of the axisymmetric nozzle Fig. 3 Comparison of pressures on the boundary of the gas-phase flow field
(a) 直径1 μm的粒子 (b) 直径2 μm的粒子 (a) Particles with a diameter of 1 μm (b) Particles with a diameter of 2 μm
(c) 直径10 μm的粒子 (d) 直径20 μm的粒子 (c) Particles with a diameter of 10 μm (d) Particles with a diameter of 20 μm
图4 不同直径粒子轨迹比较
Fig. 4 Comparison of the particle trajectories with different diameters
图5 轴线上温度分布比较
Fig. 5 Comparison of temperature distributions on the axis
图5对喷管轴线上气相温度分布进行了比较,可以看出在喷管喉道上游区域,温度差异不明显,而在喉道下游区域单相流场温度低于两相流场温度.这是由于颗粒相与气相在上游区域整体速度都较低,初始温度相同,故相互间的热交换不明显,而在下游区域颗粒相与气相相比存在速度滞后,气体在经过喉道后由于膨胀使得温度降低,而颗粒比热容较气体大,由于热惯性温度变化小,处于相同位置的粒子温度高于周围气体温度,颗粒对周围气体有加热作用,所以两相流动的气相温度高于单相流动.基于有限粒子束数目假设,颗粒相在流场中分布不连续,导致两相流场轴线上的温度分布不够光滑,但总的趋势是相同质量分数条件下颗粒尺寸越小,单相流动与两相流动之间差距越大.这是因为在相同的质量分数条件下,颗粒尺寸越小,在单位质量气体中的颗粒数目就越多,颗粒与气体的接触就越充分,相间的相互作用越剧烈,进而热交换就越明显.就该算例而言,由于入口处两相无温度差,且整个流动区域较小,故喷管出口中心位置的温度差异不大.
采用出口Mach数为5的型面喷管,目的是模拟地面电弧风洞设备中的气粒两相流动,并分析随机分布的多尺寸粒子对气相流场均匀性的影响.铜在高温有氧情况下会生成氧化铜(CuO),氧化铜的密度为6 500 kg/m3,其比热尚无公开数据,且在地面风洞中形成的氧化铜颗粒尺寸也尚不清楚.首先对颗粒物性参数做一些假设以验证所建立的方法,假定颗粒比热为500 J/(kg·K),颗粒在喷管入口初始位置采用随机分布,颗粒尺寸采用固体火箭发动机流场中试验数据的质量平均直径[19-21].通过对固体火箭发动机喷管羽流场中的颗粒进行采集和尺寸分析,结果表明,大部分粒子呈球形,表面光滑,其质量平均直径在8~11 μm,该算例模拟过程中采用这个范围的颗粒尺寸.取粒子束总数为50和200两种状态,粒子取4种直径,设定8 μm和11 μm两种尺寸的粒子束数目各占总粒子束数目的20%,9 μm和10 μm两种尺寸的颗粒各占30%,详见表1.喷管入口处颗粒相质量分数为0.1,假设各粒子束上颗粒质量流量相等.试验气体为空气,风洞驻室即喷管入口处气相总压为6 MPa,总温1 500 K,气流方向角为0°.
表1 不同尺寸的粒子束数目
Table 1 Numbers of the particle beams with different sizes
casenumber of the particle beams8 μm9 μm10 μm11 μm110151510240606040
图6给出了喷管对称面上的温度分布比较,在纯气相条件下,流场中参数等值线比较光滑,而在考虑颗粒相的影响后,流场等值线出现了明显的无规则振荡.由于喷管入口边界颗粒初始位置为随机分布,所以两相流场参数分布并非和纯气相流场一样呈严格轴对称分布.就两种粒子束数目条件比较而言,50个粒子束情形比200个粒子束情形振荡严重,而且在喷管出口附近后者的均匀区要比前者大.这是因为在相同颗粒质量流量的条件下,粒子束数目越多,每个粒子束上的质量流量减少,颗粒在整个流场中分布越广泛,颗粒相与气相的动量能量交换能够在更多的网格中进行,导致数值模拟过程中方程耦合源项在流场中分布更加连续,此时颗粒相对气相流场均匀性的影响要小一些.可以看出在此颗粒质量流量条件下,喷管出口位置的气相Mach数和温度并未受到明显影响,也就是说模拟的风洞理想纯气相流动环境并未因为颗粒相的加入而明显改变.
(a) 没有粒子 (a) Without particles
(b) 粒子束数目为50 (b) With 50 particle beams
(c) 粒子束数目为200 (c) With 200 particle beams
图6 Mach 5喷管对称面温度分布比较
Fig. 6 Comparison of temperature distributions on the symmetry plane of the Mach 5 nozzle
注 为了解释图中的颜色,读者可以参考本文的电子网页版本,后同.
在定常流动假设下对经过流场中网格体的粒子进行统计,可以获得喷管流场中的粒子数密度,该数据可以和试验测量值进行比较验证.图7给出了流向不同截面上的粒子数密度分布,实线为其中10个粒子束的运动轨迹,可以看出在该模拟条件下粒子数密度达1013 m-3.采用颗粒轨道模型模拟有限轨道上初始位置随机分布的颗粒运动,不能得到连续严格轴对称的粒子数密度分布.有些区域经过的粒子很多,有的区域甚至没有粒子经过,导致流场中粒子数密度差异很大,即使是增加模拟的粒子束总数,也不能保证能够得到很连续的分布结果.但是粒子数密度分布却可以定性地反映微小烧蚀颗粒在流场中的大致分布趋势,在喷管喉道上游区域,由于运动速度较低,所以粒子数密度很大,大量粒子如果长期在这些位置沉积或冲刷,可能对喷管喉道附近壁面造成破坏,进而影响试验效果.随着粒子束数目的增大,粒子在流场中分布更广泛,导致喷管喉道下游尤其是出口附近区域流场的均匀性更好,两相流场参数与纯气相更接近.
(a) 50个粒子束 (a) With 50 particle beams
(b) 200个粒子束 (b) With 200 particle beams
图7 喷管不同截面粒子数密度分布
Fig. 7 Particle number density distributions at cross sections of the nozzle
基于Euler-Lagrange模型建立了喷管内部气粒两相流数值模拟方法,可用于求解基于三维多块结构网格的气粒两相耦合流场,并可以同时考虑不同尺寸的颗粒运动,可应用于电弧风洞两相流场的模拟分析.研究发现:
1) 采用结构网格对流场进行模拟时,引入网格特征时间步长,可大幅减少颗粒搜索范围,提高计算效率;
2) 在入口位置颗粒分布足够广的情况下,对相同的颗粒质量分数,颗粒尺寸越小,相间相互作用越剧烈,动量能量交换越充分,在喷管相同位置处两相流动的参数与单相流动的流场参数差异越大;
3) 就设定的Mach 5近似电弧风洞喷管而言,在较小的颗粒质量分数条件下,颗粒相对喷管扩展段前段流场均匀性有一定影响,而对喷管出口流场的影响有限,且总的粒子束数目越多,影响越小.
电弧风洞内部两相流动非常复杂,目前开展的数值模拟方法还不足以准确描述其物理现象,只获得了一些定性的结论,需要结合数值模拟和风洞试验对其进一步深入研究.首先电弧风洞电极烧蚀形成的粒子物性参数,如粒子形状、尺寸、质量分数、比热等,需要开展相关试验进行微小粒子收集和特性测量;其次,电弧风洞流场很多情况下处于非平衡状态,对流场进行模拟时需要考虑热化学非平衡效应,尤其当微小颗粒进入流场后是否对化学反应速率造成影响,进而影响试验段模型表面热载荷等是下一步研究方向.
[1] 李洁. 非均相非平衡可压流动的建模与算法研究[D]. 博士学位论文. 长沙: 国防科技大学, 2004.(LI Jie. Study on models and algorithms of inhomogeneous non-equilibrium compressible flows[D]. PhD Thesis. Changsha: National University of Defense Technology, 2004.(in Chinese))
[2] 顾兆林. 风扬粉尘: 近地层湍流与气固两相流[M]. 北京: 科学出版社, 2010.(GU Zhaolin. Wind-Blown Dust: Near-Surface Turbulence and Gas-Solid Two-Phase Flow[M]. Beijing: Science Press, 2010.(in Chinese))
[3] 郑昌浩, 徐旭常. 不同Stokes数和初始条件对炉内固体颗粒弥散的影响[J]. 中国电机工程学报, 2003, 23(1): 171-176.(ZHENG Changhao, XU Xuchang. Analysis of the effect of Stokes number and initial condition on the dispersion character for the solid particles in boiler furnace[J]. Proceedings of the CSEE, 2003, 23(1): 171-176.(in Chinese))
[4] 刘莹, 张伟中, 殷艳飞, 等. 双向流固耦合作用下狭窄左冠状动脉内两相血流分析[J]. 应用数学和力学, 2016, 37(5): 501-509.(LIU Ying, ZHANG Weizhong, YIN Yanfei, et al. 2-phase hemodynamic analysis under bidirectional fluid-structure interaction in the left coronary artery with stenosis[J]. Applied Mathematics and Mechanics, 2016, 37(5): 501-509.(in Chinese))
[5] 于勇, 刘淑艳, 张世军, 等. 固体火箭发动机喷管气固两相流动的数值模拟[J]. 航空动力学报, 2009, 24(4): 931-937.(YU Yong, LIU Shuyan, ZHANG Shijun, et al. Numerical simulation of gas-particle flow in nozzle of solid rocket motor[J]. Journal of Aerospace Power, 2009, 24(4): 931-937.(in Chinese))
[6] ELLIOTT T S, MAJDALANI J. Two-phase flow stability of cylindrically-shaped hybrid and solid rockets with particle entrainment[C]//The 50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. Cleveland, OH, 2014.
[7] LI X T, TIAN H, YU N J, et al. Three-dimensional numerical simulation of two-phase flow in hybrid rocket motor[C]//The 49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. San Jose, CA, 2013.
[8] 张利珍, 李运泽, 王浚. 气固两相流风洞内颗粒运动特性研究[J]. 系统仿真学报, 2007, 19(14): 3200-3202.(ZHANG Lizhen, LI Yunze, WANG Jun. Study on particle flow character in gas/solid two-phase wind tunnel[J]. Journal of System Simulation, 2007, 19(14): 3200-3202.(in Chinese))
[9] 李中华, 李志辉, 蒋新宇, 等. NPLS流场测量技术在高超声速风洞中纳米粒子跟随性数值仿真[J]. 试验流体力学, 2017, 31(1): 73-79.(LI Zhonghua, LI Zhihui, JIANG Xinyu, et al. Numerical simulation of nano-particle following features for NPLS measurement technology used in hypersonic wind tunnel[J]. Journal of Experiments in Fluid Mechanics, 2017, 31(1): 73-79.(in Chinese))
[10] MACCORMACK R W, CANLDER G V. The solution of the Navier-Stokes equations using Gauss-Seidel line relaxation[J]. Computers & Fluids, 1989, 17(1): 135-150.
[11] SCALABRIN L C. Numerical simulation of weakly ionized hypersonic flow over reentry capsules[D]. PhD Thesis. Michigan: University of Michigan, 2007.
[12] JAMES A T, MORRIS M P. Chemically reacting one-dimensional gas-particle flows[R]. NASA-CR-147388, 1975.
[13] 张永飞. 双脉冲固体火箭发动机内流场模拟[D]. 硕士学位论文. 哈尔滨: 哈尔滨工程大学, 2010.(ZHANG Yongfei. Numerical study on dual-pulse solid rocket motor inner floe fields[D]. Master Thesis. Harbin: Harbin Engineering University, 2010.(in Chinese))
[14] HAYASHI A K, MATSUDA M, FUJIWARA T. Numerical simulation of gas-solid two-phase nozzle and jet flows[C]//The 23rd Thermophysics, Plasmadynamics and Lasers Conference. San Antonio, CA, 1988.
[15] 杨顺华. 碳氢燃料超燃冲压发动机数值研究[D]. 博士学位论文. 绵阳: 中国空气动力研究与发展中心, 2006.(YANG Shunhua. Numerical study of hydrocarbon fueled scramjets[D]. PhD Thesis. Mianyang: China Aerodynamics Research and Development Center, 2006.(in Chinese))
[16] 董长银, 栾万里, 周生田, 等. 牛顿流体中的固体颗粒运动模型分析及应用[J]. 中国石油大学学报(自然科学版), 2007, 31(5): 55-59.(DONG Changyin, LUAN Wanli, ZHOU Shengtian, et al. Analysis and application of model for solid particle movement in Newton fluid[J]. Journal of China University of Petroleum(Edition of Natural Sciences), 2007, 31(5): 55-59.(in Chinese))
[17] 胡建新, 夏志勋, 刘君. 颗粒轨道模型中颗粒跟踪与定位算法研究[J]. 弹道学报, 2005, 17(1): 82-87.(HU Jianxin, XIA Zhixun, LIU Jun. Overview of particle locating algorithm for Eulerian-Lagrangian computations of two-phase flow[J]. Journal of Ballistics, 2005, 17(1): 82-87.(in Chinese))
[18] HWANG C J, CHANG G C. Numerical study of gas-particle flow in a solid rocket nozzle[C]//The AIAA/SAE/ASME/ASEE 23rd Joint Propulsion Conference. San Diego, CA, 1987.
[19] HERMSEN R W. Aluminum oxide particle size for solid rocket motor performance prediction[J]. J Spacecraft, 1981, 18(6): AIAA 80-0035R.
[20] 张宏安, 叶定友, 侯晓. 固体火箭发动机凝聚相颗粒分布研究现状[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 Technology, 2000, 23(3): 25-28.(in Chinese))
[21] SAMBAMURTHI J K. Al2O3 collection and sizing from static firing of solid rocket motor plumes[J]. Journal of Propulsion and Power, 1996, 12(3): 598-604.
罗万清(1985—),男,工程师,博士生(E-mail: wtsluo@163.com);
梁剑寒(1972—),男,教授,博士生导师(通讯作者. E-mail: jhleon@vip.sina.com).