0 引 言
固体火箭发动机为了达到提高比冲的目的,广泛采用含金属复合推进剂。燃烧后生成大量的固体颗粒,燃气射流通常是由多种气体成分、固体颗粒形成典型的两相流动,这些颗粒既对燃气排导装置产生巨大的冲刷侵蚀作用,又对燃气射流流场产生直接影响,形成复杂的相互作用干扰流场[1,2,3]。
对于高空羽流,一般采用有效仿真稀薄气体流动的DSMC(Direct Simulation Monte Carlo)方法模拟。该方法采用有限数目的仿真分子模拟实际流场中数目巨大的真实分子。通过跟踪流场中仿真分子的运动和分子间的碰撞达到流场模拟的目的。该方法广泛地用于模拟高空火箭或太空船推进器的羽流,在描述这些流动中的气体特性时有很高的精确度[4,5,6,7]。Gallis提出了一种改进的DSMC方法[8],利用Green 函数发展的DSMC 方法适用于求解在任意分子速度分布的气相流场中颗粒所受的力和热,可以模拟包括 稀薄和化学惰性固体颗粒相在内的稀薄流动。Gallis方法只考虑气相对固体颗粒的作用,忽略颗粒相作用于气体的影响。这种假定在固体推进剂火箭羽流中是有一定缺陷的,因为在固体火箭发动机羽流中很多模拟范围内相间动量和能量的传递对气体流动特性有重大的影响。后来经Burt、Boyd等人的发展,建立了一种适用于DSMC方法的双向耦合算法，既考虑气相对固相的力和热的作用,又考虑固相颗粒对气相的作用,能够准确描述固相颗粒在稀薄过渡流中的输运过程[9,10,11,12,13,14,15]。
本文采用双向耦合技术,建立了模拟稀薄两相流的DSMC方法。对相关文献上的算例计算,模拟了气固两相羽流,验证了计算方法,为稀薄条件下的气固两相流动问题提供了一种新的研究手段。 1 数值模拟算法
两相流DSMC模拟运算法则是基于相间动量和能量从微粒特性的暂时变化的传递的解耦。
考虑气相对固体颗粒的作用。假设固体颗粒处于当地自由分子流的状态,不考虑多原子气体的振动激发,在同一个网格里每个DSMC气体仿真分子作用到一个固体颗粒上的力和热流分别为[10]:
式中,Rp为等效颗粒半径；Ng为每个计算分子所模拟的真实气体分子数目；τ为颗粒表面导热适应系数；Vc为网格体积；m为单个气体分子质量；ur为气体分子与相关颗粒的相对速度,cr是ur的绝对值；k为Boltzmann常数；Tp为颗粒温度；Λ为气体转动自由度数；erot为单个分子转动能。
考虑固体颗粒对周围气体的影响,首先要确定在每个时间步长内哪个仿真分子将与颗粒进行碰撞。对Bird的非时间计数方法进行修正,来确定与所选的颗粒可能发生碰撞的计算分子数ns。
式中,Np为一个仿真颗粒所表示的实际固体颗粒的数量；ng为与固体颗粒在同一网格里的气体仿真分子的数量；Δt为时间步长；(cr)max为网格内采样到的分子-颗粒对碰撞前最大相对速度。
一个与这个颗粒发生碰撞的给定的气体仿真分子,要么为以概率等于颗粒热适应系数τ的等温壁漫 反射碰撞,要么为概率为1-τ的镜面反射。如果发生 镜面反射,则相对速度cr在碰撞中不发生改变,碰撞后的相对速度 u *r可通过cr与单位矢量 n 相乘得到。如果发生漫反射,碰撞后相对速度 u *r围绕初始相对速度 u r的方位角ε 在[0,2π]上等概率分布。在漫反射碰撞中,碰撞后相对速度c*r不能假定为等于初始相对速度cr,而是需要通过使用“取舍”法,从如下分布函数来确定c*r的值:
式中,β为气体在颗粒温度处最可几速度的倒数,β=[m/(2kTp)]1/2。
对于漫反射多原子分子气体,碰撞后转动能erot也必须改变。漫反射双原子气体分子的转动能可计算如下:
式中,Rf为(0,1)之间的一个随机数。
整体坐标系下 u r的分量为ur、vr、wr。采用Bird二元弹性碰撞,相对速度 u *r的各分量可以由ur、vr、wr、cr和c*r、角δ和ε计算得到。有：
式中,δ为碰撞偏转角,定义为- u r与碰撞后的相对速度矢量 u *r= u *m- u *p之间的夹角,这里 u *m、 u *p分别为下面碰撞时气体分子和颗粒的绝对速度。
镜面反射碰撞的偏转角分布为:
漫反射碰撞的偏转角分布为:
可以通过“取舍”法得到δ角。
碰撞后气体分子的绝对速度为:
由于增加了大量的固体颗粒,并需要计算气-固两相的相互作用,计算量会大大增加。为了提高计算效率,本文采用了基于MPI的并行计算技术。将计算区域分解为N个子区域,分别分配给N个节点。各节点进行独立且完整的两相DSMC仿真。仿真粒子(包含气相、固相)穿过区域边界时,计算区域之间进行数据交换。 2 算例与讨论 2.1 算例验证
为了验证本文建立两相流DSMC双向耦合算法可靠性,对文献[10]的算例进行了计算。计算区域确定为0.1mm宽和20mm长的矩形区域(图 1),上下为镜面反射边界,在该边界上不会发生能量和动量交换,这样保证流动过程中只在气体分子与固体颗粒之间发生能量和动量交换。入口边界为均匀流动,出口边界为超声速边界。 模拟中气相混合物为H2、CO和N2组成的混合气体,分子数密度分别为2×1023/m3、 1×1023/m3、1×1023/m3；在入口处,气体的速度为2000m/s,温度为1000K。固相有两种铝粒子：直径分别为3×10-6m、6×10-6m,每种粒子的质量流量相等；在入口边界上,粒子的速度为1200m/s,温度为2200K,总的质量流量为13.33kg/s·m2,材料密度为3970kg/m3,热容为765J/kg·K。表面热适应系数设为0.89。计算区域划分为5000个网格,时间步 长为1.5×10-9s。流场稳定后,约有270 000个气体 仿真分子和10 000个固体仿真颗粒。
 |
| 图 1 计算区域与边界条件示意图Fig. 1 Sketch of domain and boundary types |
图 2、图 3是本文对两相流中固体颗粒与气相宏观流动参数模拟结果与文献[10]计算结果的比较情 况。对固体颗粒,在入口边界上速度比气相速度低,会受到气体的推动,速度逐渐增加(图 2)；而初始温 度比气相温度高,热能会逐渐向气相传递,温度逐渐下降。由于固体颗粒的尺寸和质量不一样,由前面介绍的方法可知,在相同条件下,固体颗粒所受到的力与颗粒半径的平方成正比,而颗粒质量与颗粒半径的 立方成正比,综合起来,颗粒的加速度与颗粒的半径 成反比。所以小尺寸的颗粒加速度比较大。 在流场中,不同尺寸的颗粒会逐渐分离,其数密度会逐渐下降。图 2表明颗粒相的速度、数密度、温度流向变化基本上属于线性变化。对气相流动,受较低速度的固相颗粒的阻挡,气体的速度沿流动方向是下降的(图 3),数密度逐渐上升。由于不断从温度较高的固相颗 粒获得能量,气相的温度沿流动方向逐渐升高。相应地,气相流动的参数变化基本上也是线性的。从图 2、图 3显示出,本文对气固两相宏观流动参数模拟结果与文献研究结果吻合很好。
 |
| 图 2 气固两相流中固体颗粒宏观流动参数计算比较Fig. 2 Comparison of parameters of particle in two phase flow |
 |
| 图 3 气固混合物两相流中气体宏观流动参数计算比较Fig. 3 Comparison of parameters of gas in two phase flow |
图 4和图 5绘出流场中不同截面上动量与能量传递速率随位置的变化关系。能量流率和动量流率通过在计算中分别采样统计通过某个截面的不同相的粒子动量和能量得到。
 |
| 图 4 两相流中动量传递速率随位置变化比较Fig. 4 Comparison of momentum transfer rates |
 |
| 图 5 两相流中能量传递速率随位置变化比较Fig. 5 Comparison of energy transfer rates |
图中显示出，不论是气相还是固相，流动中的动量流率和能量流率基本上呈线性变化过程，通过比较气相和固相的变化梯度，可以看出，两者的梯度大小大致相等，符号相反。这表明在模拟计算中，气相作用在固相的动量和能量，等于固相反作用于气相的动量和能量，即在相间相互作用时，动量和能量是守恒的。在前面介绍的双向耦合算法中，气相对固相的作用和固相对气相的作用过程是解耦的，两个过程独立计算。在计算固相对气相的作用时，其过程是一个随机统计计算过程，不能保证在每一步的计算中严格地遵守动量和能量守恒，但是在大样本数情况下统计平均时，可以保证动量和能量守恒。
通过以上的对比模拟，表明本项目所建立的气- 固两相流动DSMC双向耦合算法用于计算仿真固体 火箭发动机羽流中气相混合物与固体颗粒输运过程是可靠的。 2.2 发动机两相羽流流场计算
采用本文的方法，对某固体发动机两相羽流流场进行了模拟。羽流中气相组元的摩尔分数见表 1。
颗粒数分布按颗粒半径的连续分布密度可以近似地由对数正态分布曲线来拟合[15]。即：
式中,p为颗粒尺寸的平均半径；σ为标准偏差，表示颗粒尺寸分布的分散程度。
本文采用有限分组方法，把Al2O3固体颗粒分成5组不同尺寸和质量的颗粒，作为5种颗粒来处理，5组颗粒数密度相等。每一组颗粒具有相同的尺寸和质量，颗粒直径在20μm~200μm之间。具体尺寸见表 2。
发动机出口马赫数为1.5，出口温度为2960K。出口速度设为均匀流动，其中固体颗粒的质量流量占总流量的30%。计算中，环境设为真空边界条件。图 6是两相羽流的流场数密度分布。从气相流场看，由于喷流出口马赫数较低，只比声速稍高，喷流向真空膨胀十分迅速，数密度和温度均迅速下降。固相颗粒的扩散是气相扩散引起的，图中显示，固相颗粒的扩散是有界的。当气相膨胀到密度很低时，对固相颗粒的影响很小，颗粒会保持自己的运动状态做匀速直线运动。从数密度分布来看，固体颗粒比较集中地分布在核心区域，这表明气流对固体颗粒的扩散作用是有限的。
 |
| 图 6 两相羽流流场数密度分布Fig. 6 Number density distribution of two phase plume |
温度的分布比较复杂(图 7)，沿轴向，随着气相 的迅速膨胀，温度迅速下降，能量从固体颗粒向气相 转化，颗粒的温度也逐渐下降。沿径向，轴线附近的温度比外围的温度稍低，更外围的颗粒温度更低，这是因为5种颗粒与气相相互作用时能量交换不同(参见图 8)，导致颗粒的温度不同；同时动量交换也不相同，颗粒扩散的程度因此会有所差异，这些因素导致 温度的分布比较复杂。不过整个流场温度分布差别 不大，最大差异在300K左右。轴向速度表明，外围颗粒被气相加速较大，这是因为外围主要是小尺寸的颗粒，尺寸颗粒越小，越容易被气相加速。
 |
| 图 7 羽流中固体颗粒温度分布Fig. 7 Temperature distribution of particle in plume |
 |
| 图 8 不同固体颗粒沿轴线参数分布Fig. 8 Particle′s parameter distribution along axis |
图 8是喷流轴线上几种粒子的参数分布比较。由于气相膨胀很快，在一定距离之外，气相密度变得很低，对固体颗粒的影响已经很小。本文的计算中，在x=5(以喷流半径为参考无量纲化)以内，气相对固体颗粒有较大的影响，颗粒温度下降较快,x=5以外，温度下降较小，基本上保持在一个温度上。颗粒尺寸越小，温度变化越大。在x=5以内，颗粒的速度快速增加，在 x=5以外，颗粒速度基本上保持在同一个速度上。颗粒尺寸越小，速度变化越大。 3 结束语
本文在DSMC方法的基础上，采用双向耦合技术，建立了适用于DSMC方法的两相流的数值模拟方法，通过相间作用的解耦处理，实现了气固两相动量和能量相互作用的模拟。通过对有关算例的计算表明，本文建立的方法虽然在每一个时间步长内不能保证相间动量和能量交换时守恒，但是在大样本统计平均条件下，能够保证相间动量和能量交换时守恒。同时模拟也表明，固体颗粒尺寸、气体稀薄程度等条件对相间的相互作用有很大的影响。本文的方法为稀薄两相流的计算提供了一种数值模拟手段。
| [1] | Geisler R L. A global view of the use of aluminum fuel in solid rocket motors[R]. AIAA 2002-3748. 38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference. Indianapolis, 2002. |
| [2] | Dettleff G. Plume flow and plume impingement in space technology[J]. Progress in Aerospace Sciences, 1991, 28(1):1-71. |
| [3] | Simmons S. Rocket exhaust plume phenomenology[M]. EI Segundo, California:The Aerospace Press, 2000. |
| [4] | Bird G A. Molecular gas dynamics and the direct simulation of gas flows[M]. London:Oxford Univ. Press, 1994. |
| [5] | Gimelshein S F, Boyd I D, Ivanov M S. Modeling of internal energy transfer in plume flows of polyatomic[R]. AIAA 99-0738, 1999. |
| [6] | Ivanov M S, Khotyanovskyy D V, Kudryavtsev A N, et al. Numerical study of backflow for nozzle plumes expanding into vacuum[R]. AIAA 2004-2687. |
| [7] | Huang Lin, Nie Wansheng, Chen Weifang. Studying of multi-plume interference effects for attitude-control thruster with DSMC method[J]. Acta Aerodynimca Sinica, 2003, 21(1):104-108.(in Chinese)黄琳, 聂万胜, 陈伟芳. 姿控发动机高空羽流流场干扰效应的DSMC方法研究[J]. 空气动力学学报, 2003, 21(1):104-108.[LM] |
| [8] | Gallis M A, Rader D J, Torczynski J R. DSMC simulations of the thermophoretic force on a spherical macroscopic particle[R]. AIAA 2001-2890. |
| [9] | Gallis M A, Torczynski J R, Rader D J. An approach for simulating the transport of spherical particles in a rarefied gas flow via the direct simulation Monte Carlo method[J]. Physics of Fluids, 2001, 13(11):3482-3492. |
| [10] | Burt J M, Boyd I D. Development of a two-way coupled model for two phase rarefied flows[C].42nd AIAA Aerospace Sciences Meeting and Exhibit.Reno, Nevada, 2004. AIAA 2004-1351. |
| [11] | Burt J M, Boyd I D. Monte Carlo simulation of a rarefied multiphase plume flow[C]. 43rd AIAA Aerospace Sciences Meeting and Exhibit. Reno, Nevada, 2005, AIAA 2005-964. |
| [12] | Sergey F Gimelshein, Alina A Alexeenko, Dean C Wadsworth, et al. The influence of particulates on thruster plume/shock layer interaction at high altitudes[C]. 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, 2005. AIAA 2005-766. |
| [13] | Burt J M, Boyd I D. Particle rotation effects in rarefied two-phase plume flows[C]. 24th International Symposium on Rarefied Gas Dynamics. Monopoli, Italy, 2004. |
| [14] | Gosse S, Sarou Kanian V, Veron E, et al. Characterization and morphology of alumina particles in solid propellant subscale rocket motor plumes[C]. 36th AIAA Thermophysics Conference. Orlando, FL, 2003. AIAA 2003-3649. |
| [15] | Li Jie, Ren Bin, Chen Weifang. Modeling and numerical simulation of gas-solid jet in the transitional-rarefied regime[J]. Acta Aerodynamica Sinica, 2005, 23(4):484-489.(in Chinese)李洁, 任兵, 陈伟芳. 稀薄流过渡区气固两相喷流的建模与数值模拟[J]. 空气动力学学报, 2005, 23(4):484-489. |