简介面向需要开展稀薄气体动力学仿真的高校学生与科研人员这份基于DSMC算法的Matlab工具包提供了一种微观尺度的气体模拟方案。压缩包体积仅5KB共包含3个文件Matlab主程序.m用于实现DSMC碰撞与运动逻辑说明文档.md介绍使用方法与算法思路许可证文件则明确代码使用条款整体结构简洁清晰。代码采用参数化编程关键参数可灵活调整并配有详细注释便于初学者理解算法流程也方便进阶用户根据具体工况进行二次开发案例数据可直接运行帮助读者快速获得模拟结果。目前已有79人学习下载适配Matlab 2014a、2019b、2024b等多个版本适合作为课程设计、期末大作业或毕业设计的参考资料。1. 稀薄气体模拟为什么非 DSMC 算法不可当气体密度低到分子平均自由程和流动特征尺度可比时连续性假设就崩塌了。你拿 Navier-Stokes 去算高空中飞行器绕流算出来激波位置不对壁面热流也可能差一倍很多时候不是网格不够细而是流体微元本身不存在。DSMCDirect Simulation Monte Carlo是工程上处理这类稀薄气体问题的主流手法把成千上万模拟分子按真实运动规则推进用随机数处理碰撞再统计平均还原宏观流场。标题里那个.zip多半是项目包或课程文件但真正值钱的是算法流程、参数标定和验证手段。下面的内容会从理论依据讲到能在笔记本上跑完的最小实现并把几个容易让结果崩坏的参数坑点摊开说。2. DSMC 算法原理与适用范围从 Boltzmann 方程到粒子碰撞模型2.1 DSMC 求解的对象不是流体微元DSMC 的物理起点是 Boltzmann 方程而不是 Navier-Stokes 方程。Boltzmann 方程描述单粒子分布函数 f(x, v, t) 在相空间中的演化碰撞积分是求解难点。DSMC 把碰撞积分拆成“在一个网格内随机挑选粒子对、以一定概率执行碰撞”的蒙特卡洛过程。这个拆法的前提是分子混沌假设同一网格内任意两个分子的速度分布互不相关可以随机配对。在稀薄气体里这个假设基本成立在稠密区网格内分子间存在速度关联DSMC 的碰撞频率会偏大宏观量容易出现系统性偏差。2.2 Knudsen 数什么时候选 DSMC 而不是 CFD判断稀薄程度的标准是 Knudsen 数Kn λ/L其中 λ 是平均自由程L 是流动特征长度。平均自由程随密度增大而变小因此同一个几何在海拔高和气压低的场景下稀薄程度完全不同。工程上一般这样划分范围控制方程常用方法Kn 0.01Navier-Stokes 无滑移边界CFD 有限体积/有限元0.01 Kn 0.1Navier-Stokes 滑移边界CFD Maxwell 滑移模型0.1 Kn 10Boltzmann 方程过渡区DSMCKn 10自由分子流DSMC 关闭碰撞或射线法注意DSMC 并不是只能用在 0.1 到 10 这个区间。Kn 很小时把每个时间步内的碰撞次数加大DSMC 也能给出和 CFD 一致的结果但计算成本成倍增加没人会这么干。反过来Kn 很大时碰撞频率趋近于零DSMC 退化成纯粒子跟踪依然可以处理“几乎没有碰撞”的稀薄羽流。正是这种从自由分子流到滑移区都能覆盖的能力让 DSMC 在高空飞行器、真空设备、微尺度喷管等 CFD 的空白区里不可替代。2.3 碰撞模型硬球与 VHS 的差别碰撞模型决定碰撞截面和速度分布。最简单的硬球模型把分子看成固定半径的弹性球碰撞截面是常数程序好写但无法复现粘度随温度的变化。工程代码里更常用 VHSVariable Hard Sphere模型碰撞截面按相对速度的幂次变化σ_T(g) σ_T,g_ref × (g / g_ref)^(1 - 2ω)这里 g 是分子对相对速度ω 是粘度温度指数硬球对应 ω0.5氩气约 0.80氧气约 0.73。DSMC 实现时 VHS 不改变碰撞搜索流程只把配对数公式里的碰撞截面换成这个速度相关值。别小看这个改动在强温度梯度算例里硬球模型给出的热流和温度剖面会有明显偏差而 VHS 在相同网格和粒子数下就能对上实验数据。2.4 宏观量统计从离散粒子到连续场DSMC 的最终输出是宏观量。给定一个网格单元体积为 V_cell权重系数为 W每个模拟粒子代表 W 个真实分子单元内粒子数为 N_cell数密度就是 n W·N_cell/V_cell。宏观速度是粒子速度的平均温度来自速度方差。在二维模拟里温度与平均平动动能的关系是 T m·⟨c²⟩/(2kB)其中 c 是粒子相对宏观速度的热运动速度三维则是 m·⟨c²⟩/(3kB)。初学者常把自由度数漏掉直接套三维公式结果温度偏高 50%。这里也埋着 DSMC 统计噪声的根源温度是速度方差方差估计需要大量样本单元内粒子数和采样时间步数共同决定噪声水平后面参数设置章节会专门展开。3. 搭建可运行的 DSMC 模拟程序从主循环到最小代码3.1 DSMC 主循环的五个阶段DSMC 一个时间步分成五件事粒子运动、边界处理、碰撞对选择、碰撞执行、宏观量取样。时间步长 dt 要小于平均碰撞时间这样每个时间步内“先移动再碰撞”的拆解误差才可接受。伪代码如下for step in range(steps): for particle in particles: particle.x particle.vx * dt apply_boundary_conditions(particles) rebuild_cell_list(particles) for cell in mesh.cells: select_collision_pairs(cell) exec_collisions(cell) sample_macroscopic_fields(particles)运动与碰撞分开处理是 DSMC 能绕开复杂碰撞积分的关键。真实情况下分子边飞边撞但在稀薄气体里平均自由程大于网格尺度拆分后每个时间步内的轨迹误差可以接受。如果密度升高到 DSMC 适用范围的边缘这个拆分误差会变大表现在宏观量上就是输运系数偏大。3.2 粒子与网格的最小数据结构工程代码里我不会给每个粒子写一个类尤其当粒子数超过百万时。常见做法是把坐标和速度放在连续数组里再单独存一个网格编号x np.empty(n_particles) y np.empty(n_particles) vx np.empty(n_particles) vy np.empty(n_particles) cell_id np.empty(n_particles, dtypenp.int64)碰撞时按cell_id排序粒子索引把同一个网格里的粒子连续放在一起遍历缓存友好。用 C 实现时也建议用结构体数组SoA而不是数组结构体AoS同样的算法SoA 在百万粒子规模下经常比 AoS 快 30% 以上。3.3 一段可直接运行的最小松弛算例代码下面是一段二维方盒子温度松弛模拟。盒子左半边初始温度高右半边低粒子通过 DSMC 碰撞逐步均匀化。代码去掉了真实气体常数用无量纲温度和碰撞系数展示主流程但碰撞搜索、反射边界和统计部分都是完整可用拿回去改参数即可运行。import numpy as np from numpy.random import default_rng def reflect(pos, vel, L): 反射边界: 粒子越界后速度反向, 位置折叠回盒内 mask pos 0 pos[mask] -pos[mask] vel[mask] -vel[mask] mask pos L pos[mask] 2.0 * L - pos[mask] vel[mask] -vel[mask] def collide_cell(vel, vr_max, rng, relax_coeff): 对一个网格单元内的粒子执行 DSMC 碰撞 n len(vel) if n 2 or vr_max 1e-12: return vel # 配对数上限, 真实代码里替换成 NTC 公式即可 npair_max int(0.5 * n * (n - 1) * relax_coeff) for _ in range(npair_max): i, j rng.choice(n, size2, replaceFalse) vi, vj vel[i], vel[j] vr vi - vj vr_mag np.hypot(vr[0], vr[1]) # 接受概率正比于相对速度 if vr_mag / vr_max rng.random(): continue theta rng.uniform(0.0, 2.0 * np.pi) vr_new vr_mag * np.array([np.cos(theta), np.sin(theta)]) cm 0.5 * (vi vj) vel[i] cm 0.5 * vr_new vel[j] cm - 0.5 * vr_new return vel def dsmc_relax(n_particles12000, nx20, ny20, L1.0, dt2e-3, steps300, seed42): rng default_rng(seed) x rng.random(n_particles) * L y rng.random(n_particles) * L # 左半边初始温度 4, 右半边初始温度 1 (以 kBm1 的约化单位计) vx np.where(x L/2, rng.normal(0.0, 2.0, n_particles), rng.normal(0.0, 1.0, n_particles)) vy np.where(x L/2, rng.normal(0.0, 2.0, n_particles), rng.normal(0.0, 1.0, n_particles)) dx, dy L / nx, L / ny vr_max 6.0 for step in range(steps): x vx * dt y vy * dt reflect(x, vx, L) reflect(y, vy, L) cx np.floor(x / dx).astype(int) cy np.floor(y / dy).astype(int) # 逐网格收集粒子并执行碰撞 for ix in range(nx): for iy in range(ny): idx np.nonzero((cx ix) (cy iy))[0] if idx.size 1: vel np.column_stack([vx[idx], vy[idx]]) vel collide_cell(vel, vr_max, rng, relax_coeff0.15) vx[idx] vel[:, 0] vy[idx] vel[:, 1] if step % 50 0: # 二维温度: T mean(0.5 * c^2), 已在公式中消去质量和玻尔兹曼常数 T_left np.mean(0.5 * (vx[x L/2]**2 vy[x L/2]**2)) T_right np.mean(0.5 * (vx[x L/2]**2 vy[x L/2]**2)) print(fstep{step:4d} T_left{T_left:.3f} T_right{T_right:.3f}) if __name__ __main__: dsmc_relax()逻辑说明四个连续数组承载全部粒子状态collide_cell在一小段连续内存上完成配对和速度更新。代码里的relax_coeff0.15只是教学参数它把碰撞截面、数密度和时间步长打包成了一个浮点数真正的 NTC 公式会让配对数随粒子数密度、VHS 截面和 dt 变化。vr_max在这里取常数 6.0是该算例初始速度分布的最大相对速度估算值。如果初始速度从单侧 Mach 数很大的激波出发vr_max会更大写死这个值会让碰撞接受概率偏低松弛偏慢这时需要在每个单元内动态扫描一遍相对速度。3.4 调试时重点看什么调试 DSMC 程序第一步不是盯宏观量而是检查守恒。记录每个时间步的总粒子数和总动量。边界反射写错时粒子数可能不变但总动量会漂移外在表现是温度缓慢下降或出现整体移动。第二步检查碰撞频率频率过高或粒子在不同网格间来回穿越宏观量振荡剧烈先缩dt再调配对数系数。第三步看统计噪声DSMC 的宏观量天然带涨落单元内样本不够时温度曲线看起来像噪声这时的正解是增加粒子数而不是加密网格。4. DSMC 算法参数怎么设时间步长、网格尺度、碰撞频率与统计噪声4.1 时间步长和网格尺度的两个约束DSMC 的稳定性由两个物理尺度控制。时间步长必须远小于平均碰撞时间工程上一般取 dt ≈ τ_coll/5 到 τ_coll/10其中 τ_coll λ/v_thv_th 是热运动速度。网格尺度必须小于平均自由程通常取 Δx λ/3。如果网格比平均自由程大同一个网格内的分子会在进入网格的一瞬间就被强制可能与任何邻居发生碰撞空间关联被抹平输运系数虚高。如果 dt 太大分子一个时间步穿过多个网格粒子“从哪个单元出发”变得模糊碰撞对不同空间位置的颗粒被误判为同室。4.2 NTC 配对数公式和参数速查表NTCNo Time Counter是主流 DSMC 碰撞搜索算法。每个网格里的配对数上限为npair_max 0.5·N_cell·(N_cell−1)·(σ_T·g_max·W·dt)/V_cell其中 W 是单粒子权重g_max 是当前单元最大相对速度σ_T 是 VHS 碰撞截面。这里的g_max必须是单元内最大值而不是整个流场。实际编码时我会在碰撞前先做一次轻遍历取该网格内相对速度的最大值或者用上一时间步的值加上 20% 裕量。用全局均值顶替会让接受概率失真最终松弛曲线变慢。参数调整速查表参数过大过小参考范围dt碰撞步内空间关联强计算量浪费τ_coll/10Δx碰撞对空间关联失真每单元样本不足λ/3 ~ λ每网格粒子数碰撞配对数涨大宏观量噪声大15 ~ 40权重 W粒子数少噪声大计算量不必要地大按真实数密度/粒子预算配对数上限同对反复碰撞碰撞被截断NTC 公式严格计算每网格 1540 个粒子是经验值。少于 15宏观温度噪声肉眼可见多于 40单步碰撞配对计算量按 Nc² 增长但统计精度提升有限。权重 W 和粒子数 N_particles 要一起定先根据几何尺度和数密度算真实分子数再除以你愿意承担的模拟粒子数得到 W。4.3 碰撞子程序的常见数值陷阱最容易翻车的是同对粒子在一个时间步里被选中两次。DSMC 假设一个时间步内每对分子最多碰撞一次配对循环次数超过 N_cell(N_cell−1)/2 后重复碰撞会让松弛速度虚高。解决办法是控制配对上限并且在循环里加小标记表记录本步已碰撞的粒子对。另一个陷阱是二维模拟直接用三维碰撞散射公式本代码里 2D 散射只旋转相对速度方向到一个圆周角如果你把三维球面散射搬进二维温度松弛结果不会错太多但碰撞后速度分布各向异性会留下偏差宏观上出现假切应力。4.4 边界条件漫反射和热适应系数静止等温壁面不能只把粒子速度反射回来否则壁面会“冷却”或“加热”气体。常见做法是漫反射反射速度按壁面温度对应的 Maxwell 分布采样切向速度各向同性。这个模型隐含热适应系数为 1即壁面完全让分子忘记入射速度。工程上更精确的是 CLLCercignani-Lampis-Lord反射模型它对法向和切向分别给能量适应系数 α_n、α_t默认 1 就是漫反射改小可以模拟镜面反射占比。相邻壁面如果温度不同边界处理对整体热流影响非常大务必单独写测试用例别和壁面几何混在一起调试。5. DSMC 算法验证技巧、并行化思路与工程打包习惯5.1 低成本验证温度松弛曲线对照写完 DSMC 程序先跑一遍第 3 章的温度松弛算例。初始左右温度 T1、T2 各占一半质量平衡温度由能量守恒直接算出。在二维约化单位里T_eq (T1T2)/2如果模拟取均值时权重相同可以用这个数对照最终温度。误差 1% 内基本说明碰撞和统计逻辑正确。下一步再跑一维激波管或平板 Couette 流把温度剖面、速度剖面和文献数据对比。真正考验算法的是近连续区Kn 约 0.050.1这个区间壁面滑移和热流对参数最敏感。5.2 并行化空间分解比粒子分解可靠DSMC 天然适合空间分解。把计算域切成若干块每个 MPI 进程只负责自己的块进程间只在块边界交换跨块粒子。粒子分解需要全局通信来匹配碰撞对几乎没人用。负载均衡是个难点稀薄气体在几何阻塞或激波附近密度差异可能超过一个数量级常见做法是每几百步按当前粒子数重新分配分区边界。并行版本里权重 W 和 VHS 参数在全局必须保持一致否则不同进程统计的宏观量合不起温度平均后出现非物理跳变。5.3 把工程包归档成 zip 目录时的好习惯回到标题里的.zip。一个能长期使用的 DSMC 工程包我会把目录分成src/、input/、output/、verification/四层。输入文件里显式写出dt、dx、dy、每网格目标粒子数、权重因子、VHS 直径和粘度指数绝不把参数埋在代码里。单位系统单列一个配置文件统一 SI 或约化单位混用是 DSMC 代码最常见的“跑起来全是数但全都不对”的原因。再把第 5.1 节提到的松弛测试脚本放进verification/每次改完碰撞模型或并行分区自动读 output 目录里的时间序列和理论平衡温度比较超限即非零退出码。把这些文件和验证脚本一起打进 zip后续维护时能省掉大量重复试错。本文还有配套的精品资源点击获取