第3期 林林等：粗糙微通道内稀薄气体流动与换热的蒙特卡洛直接模拟 ,385. 成为解决此类问题的非常有效的方法 2网格划分 1物理问题 由于通道的对称性，为了减少计算量，取通道的 考察无限大平行平板之间构成的微通道，光滑 一半为计算区域.对计算区域采用Delaunay.三角化 通道微流的物理模型如表1所示，由于通道的宽度 方法进行网格划分，并对下壁面及粗糙元附近增加 远远大于高度和长度，因此忽略宽度方向的流动和 了网格点，以对网格进行加密，通道中心线处所设 换热，简化为二维微通道稀薄气体流动问题；微通道 定的网格点的距离为30m,而下壁面粗糙元附近 高H=0.9m,长度L=3.0m介质为氩气，模拟 的网格点间距设定为3m不同工况所取的边界节 分子所采取的作用势为VHS(variable hand sphere) 点数如表2所示. 模型).气流入口温度300K壁面与粗糙元表面 表2工况1一7中所取边界节点的个数 的温度均恒定为350K:通道的入口压力为101325 Tabl 2 Number of boundary nodes in Cases I to 7 Pa出口压力为入口压力的一半，每个模拟分子所 工况编号 1234567 代表的真实分子的个数为8.343×10. 边界节点数1190138515141518176016422018 通道内气体的平均Kn数可以根据下式计 算： 图1所示的是光滑微通道内工况1)前0.6m Kn-1 2xdnl 1 处的网格划分情况，图2所示的是工况5中通道内 (1) 某粗糙元处网格的划分，由图中可以看出，网格的 其中，d为分子的直径；入为分子的平均自由程；L为 “光滑程度和密度较为理想， 通道的特征尺寸，即通道的宽度；n表达式为： 片 (2) 其中P为压力，PK为波尔兹曼常数，小K:T为 反射分子的温度，K 0.2 Kn数是微流动中用来衡量气体稀薄程度的一 个重要参数，根据Kn数大小范围，可以将气体流动 划分为四个不同区域：Kn<0.001,连续流动区； 0.001<Kn<0.1,滑移流动区；0.1<Kn<10过渡 02 0.4 0.6 X/um 流动区：Kn>l0自由分子流动区. 图1工况1中计算区域前0.6“m内网格的划分 由式(1)和式(2)可以计算出通道内的平均Kn Fig 1 Mesh distrbution in Case 1(top 0.6m) 数约为0.22，属于过渡流动区. 为模拟微通道的人工粗糙元，假设通道的上、 0.04E 下表面均匀分布着若干三角形粗糙元，各工况中粗 糙元的高度、宽度以及分布密度见表1表1中相对 0.02 粗糙度定义为粗糙元的面积之和与通道总面积的比 值 表1各工况粗糙元的大小及分布情况 1.45 1.50 Table 1 Size and distrbution of mough elements for various cases X/um 工况 粗糙元粗糙元粗糙元 粗糙元 相对 图2工况5中通道内某粗糙元附近网格的划分 编号 高度/m宽度/m个数 间距/m粗糙% Fig 2 Mesh distribution near a mough elment n Case 5 1 0 0 2 20 15 102 0.278 3模拟方法的验证 3 0 西 450 0.544 4 0 呢 骨 870 1.11 编制了适合于非结构化网格的DSMC程序，在 5 0 0 9 30 2.18 运行5万个时间步长以后视为程序稳定并进行取 6 60 45 25 72 2.50 样.采用如下方法对程序计算结果进行了验证， 7 60 45 49 15 4.90 图3中比较了光滑微通道（工况1计算结果与第 3期 林 林等： 粗糙微通道内稀薄气体流动与换热的蒙特卡洛直接模拟 成为解决此类问题的非常有效的方法． 1 物理问题 考察无限大平行平板之间构成的微通道．光滑 通道微流的物理模型如表 1所示．由于通道的宽度 远远大于高度和长度‚因此忽略宽度方向的流动和 换热‚简化为二维微通道稀薄气体流动问题；微通道 高 Ｈ ＝0∙9μｍ‚长度 Ｌ＝3∙0μｍ；介质为氩气‚模拟 分子所采取的作用势为 ＶＨＳ（ｖａｒｉａｂｌｅｈａｒｄｓｐｈｅｒｅ） 模型 ［13］．气流入口温度 300Ｋ‚壁面与粗糙元表面 的温度均恒定为 350Ｋ；通道的入口压力为 101325 Ｐａ‚出口压力为入口压力的一半．每个模拟分子所 代表的真实分子的个数为 8∙343×10 7． 通道内气体的平均 Ｋｎ数可以根据下式计 算 ［13］： Ｋｎ＝ λ Ｌ ＝ 1 2πｄ 2ｎＬ （1） 其中‚ｄ为分子的直径；λ为分子的平均自由程；Ｌ为 通道的特征尺寸‚即通道的宽度；ｎ表达式为： ｎ＝ Ｐ ＫＴ （2） 其中 Ｐ为压力‚Ｐａ；Ｋ为波尔兹曼常数‚Ｊ·Ｋ —1；Ｔ为 反射分子的温度‚Ｋ． Ｋｎ数是微流动中用来衡量气体稀薄程度的一 个重要参数‚根据 Ｋｎ数大小范围‚可以将气体流动 划分为四个不同区域：Ｋｎ ＜0∙001‚连续流动区； 0∙001＜Ｋｎ＜0∙1‚滑移流动区；0∙1＜Ｋｎ＜10‚过渡 流动区；Ｋｎ＞10‚自由分子流动区． 由式 （1）和式 （2）可以计算出通道内的平均 Ｋｎ 数约为 0∙22‚属于过渡流动区． 为模拟微通道的人工粗糙元‚假设通道的上、 下表面均匀分布着若干三角形粗糙元‚各工况中粗 糙元的高度、宽度以及分布密度见表 1．表 1中相对 粗糙度定义为粗糙元的面积之和与通道总面积的比 值 ［14］． 表 1 各工况粗糙元的大小及分布情况 Ｔａｂｌｅ1 Ｓｉｚｅａｎｄｄｉｓｔｒｉｂｕｔｉｏｎｏｆｒｏｕｇｈｅｌｅｍｅｎｔｓｆｏｒｖａｒｉｏｕｓｃａｓｅｓ 工况 编号 粗糙元 高度／ｎｍ 粗糙元 宽度／ｎｍ 粗糙元 个数 粗糙元 间距／ｎｍ 相对 粗糙／％ 1 2 3 4 5 6 7 — 20 20 40 40 60 60 — 15 15 30 30 45 45 0 25 49 25 49 25 49 — 102 450 870 30 72 15 0 0∙278 0∙544 1∙11 2∙18 2∙50 4∙90 2 网格划分 由于通道的对称性‚为了减少计算量‚取通道的 一半为计算区域．对计算区域采用 Ｄｅｌａｕｎａｙ三角化 方法进行网格划分‚并对下壁面及粗糙元附近增加 了网格点‚以对网格进行加密．通道中心线处所设 定的网格点的距离为 30ｎｍ‚而下壁面粗糙元附近 的网格点间距设定为 3ｎｍ．不同工况所取的边界节 点数如表 2所示． 表 2 工况 1～7中所取边界节点的个数 Ｔａｂｌｅ2 ＮｕｍｂｅｒｏｆｂｏｕｎｄａｒｙｎｏｄｅｓｉｎＣａｓｅｓ1ｔｏ7 工况编号 1 2 3 4 5 6 7 边界节点数 1190 1385 1514 1518 1760 1642 2018 图1所示的是光滑微通道内 （工况1）前0∙6μｍ 处的网格划分情况．图 2所示的是工况 5中通道内 某粗糙元处网格的划分．由图中可以看出‚网格的 “光滑 ”程度和密度较为理想． 图 1 工况 1中计算区域前 0∙6μｍ内网格的划分 Ｆｉｇ．1 ＭｅｓｈｄｉｓｔｒｉｂｕｔｉｏｎｉｎＣａｓｅ1（ｔｏｐ0∙6μｍ） 图 2 工况 5中通道内某粗糙元附近网格的划分 Ｆｉｇ．2 ＭｅｓｈｄｉｓｔｒｉｂｕｔｉｏｎｎｅａｒａｒｏｕｇｈｅｌｅｍｅｎｔｉｎＣａｓｅ5 3 模拟方法的验证 编制了适合于非结构化网格的 ＤＳＭＣ程序‚在 运行 5万个时间步长以后视为程序稳定并进行取 样．采用如下方法对程序计算结果进行了验证． 图 3中比较了光滑微通道 （工况 1）计算结果与 ·385·