摇 2006 年 4 月 第32 卷 第4 期 北 京 航 空 航 天 大 学 学 报 Journal of Beijing University of Aeronautics and Astronautics April摇 2006 Vol.32摇 No郾 4 DSMC 方法中的自适应当地时间步长法 杜永乐摇 摇 阎摇 超 (北京航空航天大学 航空科学与工程学院, 北京 100083) 摇 摇 摇 要: 理论分析表明,Laux 提出的 DSMC(Direct Simulation of Monte Carlo) 方法中 摇 摘摇 的当地时间步长法尽管能够显著缩短流场达到稳定所需的 CPU 计算时间,提高 DSMC 程序的 运行效率,却存在仿真分子运动和碰撞计算复杂,并需要耗费额外计算机内存的缺憾. 对 Laux 的方法中仿真分子的运动处理时机提出了改进,并改变了其碰撞抽样方法,从而简化了仿真分 子的运动和碰撞计算处理,避免了额外的计算机内存消耗. 应用改进后的自适应时间步长法, 对圆柱的稀薄气体绕流进行了采用和未采用改进的自适应当地时间步长法的对比计算. 结果 表明,改进后的自适应当地时间步长法能明显缩短流场达到稳定所需的计算时间,对流场模拟 结果产生的影响却非常小. 关摇 键摇 词: 稀薄气体; DSMC; 自适应当地时间步长 中图分类号: V 211.25 摇 文献标识码: A摇 摇 文 章 编 号: 1001鄄5965(2006)04鄄0387鄄04 摇 Adaptive local time step method for DSMC code Du Yongle摇 Yan Chao (School of Aeronautic Science and Technology, Beijing University of Aeronautics and Astronautics, Beijing 100083, China) Abstract: Theoretic analysis reveals that, ALTS(adaptive local time step) method offered by Laux, al鄄 though greatly reducing the execution CPU time needed for the flow to reach to a steady state and thus leading to high simulation efficiency of DSMC(direct simulation of Monte Carlo) method codes, has some drawbacks in complicating molecules爷 moving and collision calculations and in requiring additional RAM. Modifications were made to time intervals at which molecular moving and collision are computed, and to the collision sam鄄 pling method. These modifications simplified molecules爷 moving and collision computation, and also avoid the consumption of unnecessary RAM. The modified ALTS method was adopted in the DSMC code. Simulations of rarefied gas flow over a cylinder were performed with and without this modified ALTS method, and the result shows that these modifications result in a significant reduction in the execution CPU time, but introduce little difference on the simulation resolution and the flow structure. Key words: rarefied gas; direct simulation of Monte Carlo; adaptive local time step 踪大量仿真分子的运动和碰撞实现对真实流动过 摇 摇 DSMC(Direct Simulation of Monte Carlo) 方法 是模拟稀薄气体流动的一种统计方法,它通过跟 程的模拟[1,2],是目前模拟稀薄气体流动最成功 的方法[3 ~5]. 但 DSMC 方法往往需要消耗极大的 物理内存和 CPU 时间,从而限制了其大规模应 用. 因此,许多学者提出了并行计算[6]、更改计算 机体系结构和数据结构[6]、高效抽样函数[7]、自 摇 收稿日期: 2005鄄04鄄08 摇 作者简介: 杜永乐(1979-),男,湖北恩施人,硕士生, yongle. du@ gmail. com. 适应 当 地 时 间 步 长(ALTS,Adaptive Local Time Step)[8] 等方法来提高 DSMC 方法的运行效率. DSMC 方法将大量仿真分子的运动和碰撞视 为各计算网格单元中发生的 2 个互相独立的过 程,所有网格单元中的仿真分子的运动和碰撞计 算中采用相同的时间步长,并由具有最大粒子运 定[1,2,4]. 因此,在其他单元中这个时间步长可能 动 速 度 和 最 高 粒 子 对 碰 撞 频 率 的 网 格 单 元 决 

摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 2006 年 自网格单元 i 上一次被处理至当前时间步 Nit 已 经历的计算步数;Ei = E(Nit,Ti) 表示从当前时间 步 Nit起,至网格单元 i 下一次被处理还剩余的时 间步数. 显然,如果当前时间步 Nit 满足Fi=0,则 开始计算网格单元 i 内仿真分子的运动和碰撞. 网格单元 i 中仿真分子的下一次运动和碰撞 的计算 时 间 步 长 驻tp 通 常 可 以 简 单 地 表 示 为 驻tp=Ti驻t. 但由于仿真分子在运动中可能进入具 有不同时间区 Tj 的网格单元 j,或者网格单元 i 可 能在下一次被处理前已经产生时间区的更新,此 时 驻tp屹Ti驻t,因此必须记录每个仿真分子下一次 的运动和碰撞计算时间步长 驻tp. 当仿真分子在当前运动计算中由网格单元 i 进入网格单元 j 时,有 驻tp=Ej驻t. 当仿真分子所处网格单元 i 在下一次被处理 前发生 时 间 区 的 调 整 时,驻tp = [ F( Nit, Ti,old) + E(Nit,Ti,new)]驻t,其下一次运动时间步长为时间 区调整以前剩余的运动时间步长 F(Nit,Ti,old)驻t E(Nit,Ti,new)驻t 之和. 1 所示. 和 时 间 区 调 整 以 后 的 剩 余 运 动 时 间 步 长 该算法中仿真分子运动和碰撞计算流程如图 图1摇 Laux 的 ALTS 流程图 北 京 航 空 航 天 大 学 学 报摇 摇 摇 摇 883 显得过小,从而降低了程序的执行效率. Laux 的 ALTS 方法允许各网格单元具有不同 的时间步长,并在计算中随各网格中粒子运动状 态自动调整当地时间步长而获得最佳值,从而使 流场达到稳定所需的 CPU 计算时间大大缩短. 但 Laux 的 ALTS 方法存在额外的内存消耗,其仿真 分子的运动和碰撞处理也更复杂,本文将对这些 不足之处加以改进,然后对圆柱的稀薄气体绕流 进行采用和未采用改进的 ALTS 方法的对比计 算. 1摇 Laux 的 ALTS 方法[8] 最高碰撞频率和最大仿真分子速度的网格单元决 Laux 的 ALTS 方法不再对各网格单元中粒子 的运动和碰撞计算采用统一的时间步长,每个网 格单元被赋予一个时间区(time zone) Ti沂(1,2, …,Tmax),每隔 Ti 倍的基时间步长 驻t( 仍由具有 定)计算一次该单元中仿真分子的运动和碰撞. 当流场达到稳定后开始统计流场参数,真实 流场参数的统计显然应在同一时刻进行,因此将 所有网格单元的时间区 Ti 置为 1,各网格单元再 次按 驻t 计算仿真分子的运动和碰撞. 该方法涉及的主要技术问题如下: 1. 1摇 时间区 Ti 的调整 在计算过程中,每执行 Nadapt 步,统计单元内 最大仿真分子速度 Vi,max 和最大仿真分子对碰撞 概率 Pi,max,据此调整该单元的时间区 Ti,然后将 Vi,max和 Pi,max置零,重新统计该单元内的 Vi,max 和 Pi,max,作为下一次时间区调整的依据. 调整网格单元 i 的时间区 Ti 时需要确定4 个 参数:最大允许的时间区 Tmax,最大允许的仿真分 子对碰撞几率 Pmax,驻t,以及各单元的初始时间区 Ti,0. 通常可以设定 Pmax=1,Ti,0=1. 网格单元 i 的时间区 Ti 的调整方法为: 1) 保存当前的时间区,Ti,old=Ti; 2) 估计新的时间区:如果 Pi>Pmax,说明碰撞 计算时间步长过大,导致了非物理碰撞频率,因此 Ti,new=max(Ti,old-1,1),中止该单元时间区调整, 否则进行下一步; 3) 依 据 该 网 格 单 元 的 Vi,max 和 特 征 尺 寸 Hi,ref,确定其适合时间区 Ti,id =int(Hi,ref /(Vi,max· 驻t)). 若 Ti,idTi,old,则 Ti,new=min(Ti,old+1,Tmax). 1. 2摇 网格单元 i 内仿真分子运动的处理 如 果 定 义: F ( Nit, Ti ) = mod ( Nit, Ti ), E(Nit,Ti)= Ti-F(Nit,Ti),则 Fi = F(Nit,Ti) 表示 

摇 摇 摇 摇 摇 第4 期摇 摇 摇 2摇 改进的 ALTS 方法 摇 摇 摇 杜永乐等:DSMC 方法中的自适应当地时间步长法 983 型,物面反射模型为完全漫反射模型,碰撞抽样采 用 Bird 的 NTC 方法. 计算中考虑到转动能松弛, 但不计及振动能的激发. 对比计算在一台配置为 P4 1. 6 GHz,512 M RAM 的微机上进行. 运动和碰撞的模拟共进行 35 000步,流场参数的统计进行 8 000 步,计算每 执行100 步调整一次时间区 Ti. 计算结果如下: 3. 1摇 加 速 比 表1 为采用不同的最大时间区时的计算时间 肄 S),其中 和阻力系数 Cd(定义为 Cd=D/(0.5籽肄 V2 D 为阻力,籽肄 为自由来流密度,V肄 为自由来流速 度,而 S 为以圆柱直径为宽边、圆柱长为长边的矩 形面积). Tmax=1 对应未采用改进的 ALTS 方法的 结果. 表中的实验数据取自文献[9]. Tmax 计算时间/h 加速比 Cd 实验值 Cd 1.494 1 10 1.498 15 1.497 摇 摇 从表1 中可以看出,使用改进后的 ALTS 方 法可显著缩短模拟所需的 CPU 时间,而计算精度 却不会受到明显的影响. 需要指出的是,各网格单元适合的最大时间 区显然与网格单元的尺寸和网格中的最大的仿真 分子运动速度有关. 对不同的流动条件和不同的 网格划分,实际加速比应有所不同. 此外,ALTS 方 用,而在统计流场参数时使用基时间步长. 因此, 实际的加速比还与仿真分子运动计算和统计流场 参数的时间步之比有关. 法只能在仿真分子运动和碰撞计算时起到加速作 表 1摇 采用不同最大时间区时的计算时间 误差/% 3.49 3.23 3.29 102.3 56.5 45.6 1.00 1.81 2.24 1.548 稳定后再进行流场的统计也将影响到结果的准确 此外,适当的物面反射模型和分子模型对准 确模拟稀薄气体流动极为重要,是否在流场达到 性. 因此,采用更真实的气体分子模型,增加统计 流场前粒子运动和碰撞的计算步数,有可能提高 计算精度. 3. 2摇 流场结构对比 图 3 所 示 为 Tmax = 1( 对 应 计 算 中 未 采 用 ALTS 方法的结果)和 Tmax=15 的流场等值线对比 图,各等值线图均按完全相同的参数绘制. 从图中 可 以看出,2种方法下流场等值线几乎完全一致; 为了减少各网格单元内宏观参数的统计涨 网格单元内仿真分子下一次运动和碰撞计算中的 刻对每个网格单元内的仿真分子均进行运动和碰 由于需要记录每个分子下一次运动计算的时间步 落,DSMC 方法要求用相当数目的仿真分子模拟 真实气体的流动,因此仿真分子的状态参数的记 录是内存消耗的最主要来源. Laux 的 ALTS 方法 长,并记录每个网格单元时间区调整前及调整后 的时间区,从而带来更大的内存消耗,并使仿真分 子运动和碰撞处理变得更加复杂. 仔细分析可以发现,只要在时间区调整的时 撞计算,而不论该网格单元是否达到处理时间,各 时间步长就可以简单的表示为 驻tp = Ei驻t,这样 驻tp 就不再是必须记录的参数,同时也不必记录 各网格单元时间区调整前的时间区 Ti,old,从而减 少了不必要的内存消耗. Laux 采用 Baganoff 和 McDonald 的抽样方法 (参见文献[4]) 计算碰撞,因此需要考虑碰撞抽 样几率 Pcol>1 的情况. 本文采用 Bird 的非时间计 数器(NTC)方法,该方法自动保持了 Pcol<1,可以 只统计各网格单元 Vmax 作为时间区 Ti 更新的准 则. 因此,这一改进还简化了时间区更新的方法. 改进的 ALTS 方法计算流程如图2 示. 图2摇 本文改进的 ALTS 流程图 3摇 结果讨论 本文对圆柱的稀薄气体绕流进行了采用和未 采用改进的 ALTS 方法的对比计算. 计算状态为:来流马赫数 Ma=5. 92,克努森 数 Kn=0.016,攻角 琢=0毅. 计算选用 VHS 分子模 

093 图4 所示为 Tmax =1 和 Tmax =15 的物面压力系数 Cp 和传热系数St对比图,该图表明,2 种计算方 北 京 航 空 航 天 大 学 学 报摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 摇 2006 年 法下曲线几乎完全重合. 因此,改进的 ALTS 方法 能极大地缩短模拟时间,对计算结果带来的影响 非常小. 4摇 结摇 论 本文对圆柱的稀薄气体绕流算例采用加入和 本文对 Laux 的 ALTS 方法中仿真分子的运 动和碰撞处理时机进行了改进,并改变了其碰撞 抽样方法,从而减少了不必要内存消耗,并极大地 简化了仿真分子的运动和碰撞计算过程. 未加入改进的 ALTS 方法的 DSMC 程序进行了对 比计算. 计算结果表明:改进后的 ALTS 方法消除 了 Laux 的 方 法 中 的 不 足 之 处, 并 能 显 著 缩 短 DSMC 方法模拟稀薄气体流动时流场达到稳定所 需的 CPU 时间,对计算结果产生的影响却非常 小. 参考文献 (References) [1] Bird G A. Molecular gas dynamics[M]. Oxford: Clarendon Press, 1976 [2] Bird G A. Molecular gas dynamics and the direct simulation of gas flows[M]. Oxford: Clarendon Press, 1994 [3] 沈摇 青. 稀薄气体动力学[M]. 北京: 国防工业出版社, 2003 Shen Qing. Rarefied gas dynamics[M]. Beijing: Press of Na鄄 tional Defense Industry, 2003 (in Chinese) [4] 吴其芬, 陈伟芳. 高温稀薄气体热化学平衡流动的 DSMC 方 法[M]. 长沙: 国防科技大学出版社, 1999 Wu Qifen, Chen Weifang. DSMC method for high鄄temperature rarefied nonequilibrium gas flows[M]. Changsha: Press of Uni鄄 versity of National Defense Science and Technology, 1999(in Chinese) [5] 沈摇 青. DSMC 方法与稀薄气流计算的发展[J]. 力学进展, 1996, 21(1):1 ~12 Shen Qing. Development of DSMC method and rarefied gas flow simulations[J]. Development of Mechnics, 1996, 21(1):1 ~ 12(in Chinese) [6] Matthias Laux. Optimization and parallelization of the DSMC method on unstructured grids[R]. AIAA鄄97鄄2512, 1997 [7] Agbormbai A A. Efficient sampling schemes for DSMC computa鄄 tions [R]. AIAA鄄2002鄄2762, 2002 [8] Matthias Laux. Local time stepping with automatic adaptation for the DSMC method[R]. AIAA鄄98鄄2670, 1998 [9] Maslach G J, Schaaf S A. Cylinder drag in the transition from continuum to free鄄molecule flow[J]. The Physics of Fluids, 1963, 6(3): 315 ~321 图3摇 不同的最大时间区对流场结构影响的比较 a摇 物面压力系数 Cp b摇 物面传热系数 St 图4摇 不同的最大时间区对物面气动参数 分布影响的比较