用直接模拟蒙特卡罗法求高温混合气体的热传导系数

第 20 卷　第 4 期 1999 年 10月 宇　航　学　报 JOURNAL OF ASTRONAUTICS Vo l. 20 No . 4 Oct. 1999 收稿日期: 1998年 5月 6日, 修回日期: 1998年 11月 25日 * 国家杰出青年科学基金 ( 59725617) 资助项目 用直接模拟蒙特卡罗法求高温混合 气体的热传导系数* 崔国民　谈和平　阮立明　刘林华 (哈尔滨工业大学能源科学与工程学院·哈尔滨·150001) 　　摘　要　本文在局域热平衡条件下, 从传输现象产生的机理出发, 建立了气体热传导系数 数值实验的物理模型。采用 DSMC 法, 对高温混合气体热传导系数的实验过程进行数值模拟。 考虑了真实气体的化学反应、内部自由度的激发等, 对高温 ( 2000-10000K ) 真实气体混合物的 热传导系数的实验过程进行了数值模拟。获得了纯空气组元的真实热传导系数。经与热传导系 数实测数据 ( 473-1473K ) 和国外最新计算结果 ( 2000～10000K ) 的比较 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 明, 本文提出的粒 子模型和概率模型是正确的, 计算 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 有较好的精度。在此基础上, 对热传导系数的非平衡性 进行了分析。 主题词　真实气体　热传导系数　直接模拟蒙特卡罗法 THE STUDY ON THERMAL CONDUCTIVITY OF HIGH TEMPERATURE REAL GAS BY THE DIRECT SIMULATION MONTE-CARLOMETHOD Cui Guom in　Tan Heping　Ruan Liming　Liu Linhua ( S chool of Energy Science an d En gineering, Harbin Inst itute of T echnology·Harbin·150001) 　　Abstract　Under local-equilibr ium Hypot hesis, this paper set up the phy sical model o f numer ical exper iment on therm al conductiv ity of mix ed high-temperatur e gas. By means o f DSMC, t he computer-aided exper iment was conducted by the gener ation mechanism of the conductiv ity . T aking the chemical r eaction of r eal gas, the exciting o f int ernal fr eedom o f particles into consider ation, exper imental pro cess o f the gas-mixed conductiv ity at high temperat ur e ( 2000-10000K ) was sim ulated, and got r eal conductivit y o f pur e air components. The compar ison w ith the r eal-measur ed dat a ( 473-1473K ) and t he lat est da ta abr oad ( 2000～ 10000K ) , show ed that the par ticle model and probability model pr oposed in this paper are r ight and the computation met hod produces fair ly high precision. On the basis o f tho se above, the non-equilibr ium property of conduct ivity w as analyzed. Key words　Real g as　Thermal conduct ivity　Dir ect simulation monte-carlo 1　引言 在飞行器高超声速再入过程的研究中, 高温混合气体的传输特性直接关系到再入过程 的速度场、温度场和浓度场的精确求解。由于极高的再入温度, 使得混合气体产生明显的 真实气体效应。获取高温气体传输数据的最基本方法是实验, 由此获得的数据成为一种比 较的 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 。然而, 由于实验条件的制约, 在混合气体的种类、温度范围等方面受到很大的 限制。基于实验数据和理论分析的数据拟合法, 扩展了实验数据的温度范围, 但在精度及 混合气体的种类上仍受到约束[ 1, 2]。近年来, 随着计算机技术的发展, 直接模拟蒙特卡罗法 ( DSMC) , 得到了重视。该方法对微观气体粒子的行为进行模拟跟踪, 然后由统计结果确定 波尔兹曼方程的解, 进而求得气体的传输系数 [ 3-5]。 本文在局域热平衡条件下, 建立了热传导系数数值实验的物理模型。采用DSMC 法, 对 高温混合气体热传导系数的实验过程进行数值模拟, 其优点是: ( 1) 适应性强, 适合任意 组成的混合气体热传导系数的求解。( 2) 温度范围广, 理论上可达到任意求解温度。( 3) 充 分考虑了真实气体效应。( 4) 求解方便, 给定混合气体原始组成等初始 参数 转速和进给参数表a氧化沟运行参数高温蒸汽处理医疗废物pid参数自整定算法口腔医院集中消毒供应 后即可获得热 传导系数。 2　高温真实气体的能量传输现象 真实气体效应是指气体在高温时产生的振动激发、离解、电离和化学反应对气体特性 的影响。这些现象使气体的物理特性和完全气体产生了明显的差别。真实气体的能量传输 E 可表为 E = ET + E内 + q化学 ( 1) 式中: q化学 =∫T0c pd T + h0 对于分子: E内 = ER + EV + Ee　　　　对于原子: E内 = Ee ( 2) 其中: ET、E内、E R、E V、E e分别为气体粒子的平动能、内能、转动能、振动能、电子激发能; c p、h0 分别为粒子的定压比热和化学生成焓。其热传导系数 k c也由三部分组成: 气体粒子的平动 能贡献、内能贡献和化学反应能贡献。 根据分子运动论理论, 在高温情况下, 气体粒子的平动能E T和转动能ER 总是完全激发 的,分别是 3N kT / 2和 N kT ;但振动能和电子激发能却没有完全激发, 分别表示为: E V = N k� V / [ exp( � V / T ) - 1] ( 3) E e = N k � e g1 g0 e - � e / T / 1 + g1 g0 e - � e / T ( 4) 式中: � V , � e 分别为粒子的振动特征温度和电子特征温度; g0 , g1是简并度。 40 宇航学报 第 20卷　 3　导热系数实验过程的数值模拟 3. 1　物理模型 用直接模拟蒙特卡罗法, 从传输特性产生的机理上, 对其实验过程进行计算机辅助实 验, 建立数值实验的物理模型和数值模拟方法。它不同于真正的实验过程, 既不需要实验 设备, 也避免了苛刻的实验条件的限制, 但在机理上却符合真正的实验, 各种实验条件可 图 1　热传导物理模型 以在理论允许的条件下任意设定。 一个封闭系统长度为 L , 其内为混合气体, 温度分布为( dT / dx ) 0 线性分布,两端 A、B 处的 温度分别为 T A、T B (见图 1)。整个系统划分为相 等的n 个子域, 各子域的温度、压力相等。子域 x 向的尺度应符合直接模拟蒙特卡罗法计算的要 求,即小于粒子的平均自由程, 使在各子域内满 足局域热平衡条件。粒子与系统 A、B 两端面碰 撞后能量既不损失、也不增加。在系统中间,取一截面0-0′, 其温度为T 0。从初始时刻 t 0开始, 对整个系统各子域内的粒子进行跟踪,在单位时间步� t ( t0 + �t < t1 )后更新各子域的参数, 直到时刻 t 1,获得该时刻 0-0′截面处的温度梯度( dT / dx ) 1。在 t 1 - t0足够小,即跟踪时间很 短的情况下, 整个跟踪过程中 0-0′截面的温度梯度可认为是: ( d T / d x ) = 0. 5 �[ ( dT / dx ) 0 + ( dT / dx ) 1] ( 5) 若在此时间间隔内, 由 0-0′截面传输的能量为 E, 则有下式存在: E / ( t 1 - t0 ) = k c �dT / dx ( 6) 即: k c = E / ( t1 - t 0) �dx / dT ( 7) 对于扩散现象和粘性系数同样可以通过这种方法建立物理模型。 3. 2　概率模型 ( 1) 碰撞对选择模型　样本粒子 i与当地流场中其它种类粒子 j ( j = 1, 2, ⋯, n, n为 粒子种类数) 的碰撞频率 Z ij可以表示为: Z ij = n j d ij [ ( 8�kT ) /mij ] 1/ 2 ( 8) 式中: n j 为 j 类粒子的浓度; d ij= ( d i+ d j ) / 2, 为 i、j 类粒子的平均直径; mij= mim j / (mi+ m j ) , 为 i, j 类粒子的折合质量。样本粒子 i 与流场中 j 类粒子的碰撞几率 Pdj为: P dj = Z ij / ∑ n j= 1 Z ij ( 9) 　　选择 0～1区间的碰撞粒子随机数R i, 若: P d, 1 + P d, 2 + ⋯ + P d, i ≤ R i < P d, 1 + P d, 2 + ⋯ + Pd , i + P d, i+ 1 ( 10) 则认为碰撞粒子是 i 类粒子。 ( 2) 粒子所处位置的判断　在各个子控制域内, 认为粒子是均匀分布的。粒子的位置 x 0通过随机数 R l [ 0～1] 来确定, 即: 　　　　x 0= R l·L ( 11) ( 3) 粒子初始状态的判断　由于在两个控制域中, 所有粒子任何时刻都处于随机热运 动之中, 它可能刚刚进行一次碰撞, 也可能即将与其它粒子发生碰撞, 即粒子在它的一个 平均自由程范围内所处的位置是随机的。通过选取 [ 0, 1] 区间的随机数 R t 来确定粒子在 41　第 4期 崔国民等: 用直接模拟蒙特卡罗法求高温混合气体的热传导系数 抽样时刻所处的状态, 即: t s = R t × s-i/ u i ( 12) 其中: t s是粒子距下一次碰撞的时间, s-i 为 i类粒子的平均自由程, u i 为 i粒子的速率。 ( 4) 由于认为粒子在任一子域内处于区域热平衡状态, 因此其速率符合麦克斯韦尔速 率分布规律, 即: f = n ( 2�kT /m) 3/ 2 �exp - v 2 2kT /m ( 13) 式中: n为气体粒子的数密度; m 为粒子的质量。根据这一分布规律, 取无量纲速率 c , 使 得: c = v 2/ ( 2kT /m) ( 14) 在进行速度选择时, 只考虑速率分布在 c= - 3～+ 3范围内的粒子, 因为在此范围以外的粒 子占总粒子数的 0. 0022%以下, 其影响很小。在 [ 0, 1] 上取均匀分布的随机数 Rv , 令: c( R v) = 3 �( - 1 + 2R v)　　c1 = � c( Rv ) - 1�　　f-( c ) = exp( - c21) ( 15) 另取 [ 0, 1] 上均匀分布的随机数 R1 , 如果 R1< f- ( c) , 则 c即是粒子的相对速率, 否则从 新提取 Rv , 重复上述步骤。在平衡的条件下, 粒子处于一种无规律的运动状态, 没有优先 的运动方向。根据文献 [ 6, 7] , 速度方向由余弦定律产生, 天顶角 �、圆周角 �分别为: �= cos- 1( 1 - 2R�)　　　　�= 2�R � ( 16) 式中　R�, R �为 [ 0～1] 区间内均匀分布的随机数。通过上述方法求得的粒子速率和方向, 经验证, 其平均速率和平动温度误差范围在 0. 4%以下[ 8]。 3. 3　高温平衡空气组成 2000-10000K 温度范围内, 空气的主要组元为: O 2、O、N 2、N、NO、NO +、O +、N + 、 e - (不计含量较少的Ar , CO 2, H2O 等) , 它们的成分组成受如下反应制约[ 9] : O 2� O + O NO� NO+ + e- O� O + + e- N 2� N + N NO� N + O N� N + + e- ( 17) 式 ( 17) 的六个反应分别对应六个平衡常数 ki ( i= 1, 2⋯6, 取自文献 [ 9] pp. 311) , 分别 表示为: k 1 = [ O] 2 / [ O 2] k2 = [ N] 2 / [ N 2] k3 = [ N] [ O] / [ NO] k 4 = [ NO + ] [ e- ] / [ NO] k5 = [ O + ] [ e- ] / [ O] k6 = [ N + ] [ e- ] / [ N] ( 18) 此外还受氧、氮以及电子守恒方程的制约, 分别表示为: [ O 2] 0 = 0. 5 �[ O] + 0. 5 �[ NO] + 0. 5 �[ NO + ] + 0. 5 �[ O + ] + [ O 2 ] [ N 2 ] 0 = 0. 5 �[ N] + 0. 5 �[ NO] + 0. 5 �[ NO + ] + 0. 5 �[ N + ] + [ N 2 ] [ e - ] = [ NO + ] + [ O + ] + [ N + ] ( 19) 其中方括号表示粒子的摩尔浓度: [ O 2] 0、[ N 2 ] 0分别为不同温度下氧分子和氮分子的初始 摩尔浓度。 [ O 2] 0 = 9. 363× 10- 6× �/ �0 [ N 2] 0 = 3. 485× 10- 5× �/�0 ( 20) 式中, �、�0 分别为计算温度下和标准温度下的空气密度, �由下式获得: �= ( 1/ Zc) × P/ RT × M 0 ( 21) 42 宇航学报 第 20卷　 其中, Zc 为高温空气的压缩因子 (取自文献 [ 9] ) ; P 为压力; M 0是标准状态下空气的平 均分子量。联立式 ( 18) 和式 ( 19) 九个方程组成方程组, 共有九个未知数, 通过迭代方 法求解。 3. 4　计算步骤 ( 1) 分别求得各子域中粒子在不同温度下的浓度组成 (见 3. 3)。 ( 2) 从时刻 t0开始, 在单位时间步 �t 内( � t要小于粒子的平均碰撞频率) , 对子域内的粒子 进行跟踪。 ( a) 选择样本粒子。所选样本粒子种类由上步的平衡组成抽样决定。 ( b) 判断粒子的初始位置。 ( c) 判断粒子的初始状态并决定粒子的运动方向和速率。 ( d) 粒子运动后的位置。首先根据粒子初始状态判断粒子在时间步 �t 内是否发生碰 撞,当 ts > � t, 认为粒子未与其它粒子发生碰撞;反之则发生碰撞。若发生碰撞,则按二体碰 撞的能量守恒和动量守恒原理,确定碰撞粒子的速率和能量分配。然后重新选择粒子的运 动方向。其最终的位置可表示为: x 1 = x 0 + cos(�) ×∑N c n= 1 un× t n,式中, N c为粒子在 �t时间 内总的碰撞次数; un和 tn分别是每次碰撞后粒子的速率和自由运动时间。 ( e) 对所选取的样本粒子数 Mpar , 重复上述过程分别对各子域内的粒子进行跟踪。 ( f ) 调整各子域的参数, 包括速率、温度、粒子密度。 粒子的平均速度 ( u, v , w ) k = 1 N k �∑ N k i= 1 ( u, v , w ) i ( 22) 温度 T k = 1 N k �∑Nk i= 1 2 3 �[ ( u2i + v 2i + w 2i ) - ( u2k + v 2k + w 2k ) ] ( 23) 粒子浓度 nk = N k / V k ( 24) 其中, i表示 i粒子; u、v、w 分别为粒子三个方向的速度分量; k表示 k子域; N k为 k子域 的粒子个数; V k 代表 k 子域的体积。式中的温度、速度等参量均为无量纲量。 ( g ) 进行下一时间步的跟踪, 直到时刻 t1为止。 ( 3) 统计在时间间隔 t1- t0 内粒子通过 0-0′界面总传输能, 求得热传导系数。 4　算例及讨论 4. 1　低温下纯空气的热传导系数的验证 为了验证本计算方法的可行性, 本文对低温下 ( 473. 15K-1473. 15K) 的纯空气热传导 系数进行了计算, 并与文献 [ 10] 的数据进行了比较。不同样本粒子数 (M par= 10万、100 万、1000万、10000万) 的计算结果如表 1所示。由于本文的目的是计算高温混合体 (主 要为: O 2 , O, N 2 , N, NO, NO + , O + , N + , e- 九种气体组元) 的热传导系数, 所以在 用本方法计算低温下的纯空气热传导系数时, 只考虑了氧气和氮气, 未考虑氩、水蒸汽等 含量相对小的成份。因此从表1中可以看出 , 在温度较低时 , 本文的计算结果比文献 43　第 4期 崔国民等: 用直接模拟蒙特卡罗法求高温混合气体的热传导系数 Administrator 高亮 选择样本粒子有错。 应该以Pd,j求和后 认为是Pi 然后以Pi再选粒子，而且粒子选的时候 文中的应该是选出的事i+1类粒子。 然后 再选和哪一种粒子碰撞。 [ 10] 的数据略小。表中: er1 ( T j ) 为相对误差, er 2为均方根误差: er1( T j ) = { [ �本文计算值 ( T j ) - �文10( T j ) ] /�文10 ( T j ) } × 100　er2 = ∑ N J J= 1 er1( T j ) 2 N J ( 25) 表 1　空气热传导系数的数据比较 T j ( K) 文献 [ 10] 本文计算结果 10万 er 1 ( T j ) 100万 er1 ( T j ) 1000万 er 1 ( T j ) 10000万 er1 ( T j ) 473. 15 0. 0393 0. 03865 - 1. 663 0. 03849 - 2. 068 0. 03816 - 2. 895 0. 03820 - 2. 791 573. 15 0. 0460 0. 04456 - 3. 121 0. 04438 - 3. 530 0. 04402 - 4. 313 0. 04406 - 4. 214 673. 15 0. 0521 0. 05080 - 2. 491 0. 05058 - 2. 914 0. 05019 - 3. 671 0. 05024 - 3. 577 773. 15 0. 0574 0. 05712 - 0. 493 0. 05686 - 0. 936 0. 05644 - 1. 681 0. 05649 - 1. 589 873. 15 0. 0622 0. 06330 1. 767 0. 06301 1. 303 0. 06255 0. 563 0. 06261 0. 653 973. 15 0. 0671 0. 06920 3. 135 0. 06888 2. 655 0. 06839 1. 922 0. 06845 2. 011 1073. 15 0. 0718 0. 07474 4. 094 0. 07439 3. 603 0. 07386 2. 874 0. 07393 2. 962 1173. 15 0. 0763 0. 07985 4. 652 0. 07947 4. 154 0. 07892 3. 428 0. 07898 3. 516 1273. 15 0. 0807 0. 08450 4. 704 0. 08409 4. 203 0. 08351 3. 480 0. 08358 3. 567 1373. 15 0. 0850 0. 08865 4. 296 0. 08823 3. 796 0. 08761 3. 075 0. 08769 3. 162 1473. 15 0. 0915 0. 09230 0. 869 0. 09185 0. 386 0. 09121 - 0. 314 0. 09129 - 0. 230 er2 3. 186 2. 975 2. 844 2. 848 � O 2　○ O　△ O+　+ NO　� e- ★ N 2 ◇ N ●N + � NO+ 图 2　高温空气组成 4. 2　高温真实气体的热传导系数计算 利用前述方法, 本文对 1个大气压情况下 的纯空气的热传导系数进行了计算, 初始计算 参数为: T 0= 273. 15K, P= 1. 01325×105Pa, XN = 79%, XO= 21%, 由计算步骤 ( 1) , 获 得高温下空气组元中各种成分的组成。其结果 见图 2。在计算过程中, 样本粒子数 Mpar的选 择是很重要的。如果选择的数目过小, 将使计 算结果出现较大的涨落; 如果选择数目过大, 又将浪费大量的计算时间。鉴于此, 本文对不 同样本粒子数分别进行了计算, 其结果示于图 3中。从图中可以看出, 在样本粒子数 M par> 100万个时就可以获得较稳定的结果。因此, 本文的计算结果除特殊说明外, Mpar均取 1000 万个。 图 4是本文计算的真实混合气体 (纯空 44 宇航学报 第 20卷　 气) 的热传导系数与文 [ 11] 的比较。分析图 4并参照图 2, 热传导系数可分为六部分: ( 1) 低温区域真实热传导系数与冻结热传导系数基本相等 (见表 1) ; ( 2) 2000-4000K 区域 氧裂解, 热传导系数增大; ( 3) 4000-5000K 氧裂解结束, 氮裂解开始, 热传导系数变化不 大; ( 4) 5000-7000K 氮裂解旺盛区, 热传导系数升高; ( 5) 7000K-9000K氮裂解结束, 电 离反应还不突出, 热传导系数不升反降; ( 6) 9000K 以上电离反应开始并进入旺盛区, 热 传导系数迅速升高。 图 3　高温空气导热系数与 样本粒子数的关系 图 4　真实热传导系数 ( 1. 01325×105Pa) 与文 [ 11] 的比较 从计算过程可看出,用本方法获得的结果, a—化学反应非平衡区　b—转动和振动非平衡区 c—平动非平衡区 图 5　非平衡热传导系数分布示意图 是在气体保持粒子内部自由度平衡和化学反应 平衡状态下求取的, 因此它只适用于平衡气体。 但从某种意义上讲, 非平衡态是一个动态过程, 所涉及的热传导系数数值也是不定的, 其大小 取决于混合气体的非平衡程度, 介于冻结态与 平衡态数值之间。因此, 对于流动处于化学非平 衡态 (即流动特征时间与反应特征时间很接近) 时, 热传导系数可通过定义一个非平衡因子 � ( 0～1) 表示为: �非平衡 = �冻结 + �(�平衡 - �冻结 ) ( 26) 对于转动和振动同样有这样的关系。非平衡热 传导系数区域如图 5所示。 5　结论与讨论 本文在局域热平衡条件下, 建立了高温混合气体热传导系数数值实验的物理模型。采 45　第 4期 崔国民等: 用直接模拟蒙特卡罗法求高温混合气体的热传导系数 用直接模拟蒙特卡罗法 ( DSMC) , 考虑了真实气体的化学反应、内部自由度的激发等, 对 高温 ( 2000-10000K ) 真实气体混合物的热传导系数的实验过程进行了数值模拟。 ( 1) 高温混合气体的热传导系数, 具有明显的真实气体特性, 并且受气体非平衡程度 的影响。 ( 2) 本文的计算方法充分体现了高温混合气体传输特性的机理, 求解高温混合气体热 传导系数方便, 且具有温度范围宽、计算量小等优点。 ( 3) 经与文献 [ 10, 11] 比较表明, 本计算方法可行, 具有一定的精度。 本方法对高温真实气体热传导系数的求解, 有赖于高温气体成分组成。从图2、图4可 以看出, 由于本文求解高温空气组成时, 所用模型适应于电离度小于 1%的情况, 在 10000K 左右, 获得的组份的精度不够高, 所以在 9000K 和 10000K 时热传导系数与文献 [ 11] 有 偏差。因此, 有必要采用更精确的求解高温气体组份的计算模型。一旦此问题能得以更好 的解决, 则本文的方法将有更高的计算精度, 并可适用于求解更高温度下的热传导系数。 参 考 文 献 1　Candler G. Computation of thermo-ch emical nonequilib rium mar tian atmospher ic ent ry flow s . AIAA 90-1695, 1990 2　Bhut ta B A, Lewis C H . A new technique for low -to-h igh alt itude predict ions of ablat ive hyper sonic flowf ields . AIAA 91-1392, 1991 3　Francis J A, Alejandro L G, Bern i J A. A cons istent Bol tzmann algorithm. Ph ysical Review Letters, 1995, 74 ( 26) 4　S antos A, Montanero J M . Monte Carlo simu lation meth od for th e Ensk og equ at ion. Ph ysical Review E, 1996, 54 ( 1) : 438-444 5　Yang Y S , Th om pson C J. Con duct ivity of a lat ti ce gas model M onte-Carlo s imulat ions . Ph ysica A , 1998, 248: 185- 194 6　陈熙 . 动力论及其在传热与流动研究中的应用 . 北京: 清华大学出版社, 1996: 25-30 7　余其铮, 潘迎春, 张东辉, 季建刚, 谈和平 . 蒙特卡罗方法对各向异性介质辐射特性的模拟 . 工程热物理学报, 1996, 17 ( 1) : 96-100 8　崔国民 . 军用目标红外热像理论建模中的蒙特卡罗法和并行计算 . 哈尔滨工业大学, 工学博士学位论文, 1998: 67- 73 9　卞荫贵, 钟家康 . 高温边界层传热 . 北京: 科学出版社, 1986: 320-330 10　杨世铭 . 传热学 (第二版) . 北京: 高等教育出版社, 1987: 442-443 11　Bou los M I, Fauch ais P, Pfender E. Th ermal Plasmas-Fund amental s and Appl ications , Volum e 1. Plenum Press , New York, 1994: 413-417 46 宇航学报 第 20卷