D0I:10.13374/i.i8sm1001t53.2010.03.02 第32卷第3期 北京科技大学学报 Vol 32 No 3 2010年3月 Journal of Un iversity of Science and Technolgy Beijing Mar.2010 粗糙微通道内稀薄气体流动与换热的蒙特卡洛直接模 拟 林林)王爱国2)张欣欣” 1)北京科技大学机械工程学院，北京1000832)中国科学院工程热物理研究所，北京100190 3)中国科学院研究生院，北京100190 摘要采用Delaunay三角化方法对计算区域进行网格划分，开发了适合于非结构化网格的蒙特卡洛直接模拟程序，并对程 序的正确性进行了验证·在此基础上模拟分析了粗糙元为三角形的平行平板间微通道内稀薄气体的二维流动与换热.通道 进出口压力固定，上下平板温度恒定·计算分析了粗糙元高度、宽度以及分布密度的影响·结果表明：微通道内粗糙元对流动 与换热有明显的扰动：随着粗糙元的变大，速度跳跃显著，甚至出现漩涡，增加了通道内的压力损失：但粗糙元增强了微通道 壁面与气体之间的换热。 关键词粗糙：微通道：稀薄气体：流动：换热 分类号K124 Prediction of flow and heat transfer in rough m icro channels w ith direct sim ula- tion M onte Carlo m ethod LIN Lin,WANG Ai-guo),ZHANG X inxin) 1)School ofMechanical Engineering University of Science and Technology Beijing Beijing 100083 China 2)Institute of Engineering Themophysics Chinese Acadeny of Sciences Beijing 100190.Chna 3)Gmaduate University of Chnese Acadeny of Sciences Beijing 100190.China ABSTRACT The computational danain wasmeshed by Delaunay triangulation method A M onte Carlo direct smulation program was developed for unstructured gris and its validity was checked On that case the to-din ensional flow and heat transfer of rarefied gas in parallel plate m icro"channels w ith triangular moughness elments were siulated The pressure of the m icmochannels'inlet and outet and the temperature of the plates above and below were constant The influences of the height width and distrbution density of rough- ness elem ents were computed and analyzed The results show that moughness elements in the m icm-channels can create obvious disturb- ance for the flow and heat transfer The greater the roughness elments the more visible the oscillation of velocity even vortexes oc- cur which increases the pressure loss but roughness elements can enhance the heat transfer KEY WORDS roughness m icmo-channeks marefied gas flow:heat transfer 近年来，微电子机械系统(m ic roe lectrom echani 于壁面粗糙度对流动的影响，本文主要研究了微通 cal system,MEMS)、人工超晶格、纳米材料和生物 道内的粗糙元对于通道内气体的流动和换热的影 芯片等技术的诞生和深入发展引起了人们对微细尺 响.由于尺度的微细化及介质的稀薄效应，常规的 度流动特性的关注，由于微电子系统的散热问题， 连续性方程及N一S方程等不再适用.近年来， 微通道内稀薄气体的流动成为人们十分关心的问 DSMC(direct smulation Monte Carlo)方法已经被广 题，微通道流动和换热与常规尺度的流动有较大的 泛的应用于模拟稀薄气体动力学中的问题，并取得 区别，从许多文献来看，一个很重要的原因来自 了非常好的模拟结果[6-).蒙特卡洛直接模拟方法 收稿日期：2009-10-22 作者简介：林林(1968)男，副教授，博士，Email linli@adm in ustb edu cn
第 32卷 第 3期 2010年 3月 北 京 科 技 大 学 学 报 ＪｏｕｒｎａｌｏｆＵｎｉｖｅｒｓｉｔｙｏｆＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏｇｙＢｅｉｊｉｎｇ Ｖｏｌ．32Ｎｏ．3 Ｍａｒ．2010 粗糙微通道内稀薄气体流动与换热的蒙特卡洛直接模 拟 林 林 1） 王爱国 23） 张欣欣 1） 1） 北京科技大学机械工程学院北京 100083 2） 中国科学院工程热物理研究所北京 100190 3） 中国科学院研究生院北京 100190 摘 要 采用 Ｄｅｌａｕｎａｙ三角化方法对计算区域进行网格划分开发了适合于非结构化网格的蒙特卡洛直接模拟程序并对程 序的正确性进行了验证．在此基础上模拟分析了粗糙元为三角形的平行平板间微通道内稀薄气体的二维流动与换热．通道 进出口压力固定上下平板温度恒定．计算分析了粗糙元高度、宽度以及分布密度的影响．结果表明：微通道内粗糙元对流动 与换热有明显的扰动；随着粗糙元的变大速度跳跃显著甚至出现漩涡增加了通道内的压力损失；但粗糙元增强了微通道 壁面与气体之间的换热． 关键词 粗糙；微通道；稀薄气体；流动；换热 分类号 ＴＫ124 Ｐｒｅｄｉｃｔｉｏｎｏｆｆｌｏｗａｎｄｈｅａｔｔｒａｎｓｆｅｒｉｎｒｏｕｇｈｍｉｃｒｏｃｈａｎｎｅｌｓｗｉｔｈｄｉｒｅｃｔｓｉｍｕｌａ- ｔｉｏｎＭｏｎｔｅＣａｒｌｏｍｅｔｈｏｄ ＬＩＮＬｉｎ 1）ＷＡＮＧＡｉ-ｇｕｏ 23）ＺＨＡＮＧＸｉｎ-ｘｉｎ 1） 1） ＳｃｈｏｏｌｏｆＭｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇＵｎｉｖｅｒｓｉｔｙｏｆＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏｇｙＢｅｉｊｉｎｇＢｅｉｊｉｎｇ100083Ｃｈｉｎａ 2） ＩｎｓｔｉｔｕｔｅｏｆＥｎｇｉｎｅｅｒｉｎｇＴｈｅｒｍｏｐｈｙｓｉｃｓＣｈｉｎｅｓｅＡｃａｄｅｍｙｏｆＳｃｉｅｎｃｅｓＢｅｉｊｉｎｇ100190Ｃｈｉｎａ 3） ＧｒａｄｕａｔｅＵｎｉｖｅｒｓｉｔｙｏｆＣｈｉｎｅｓｅＡｃａｄｅｍｙｏｆＳｃｉｅｎｃｅｓＢｅｉｊｉｎｇ100190Ｃｈｉｎａ ＡＢＳＴＲＡＣＴ ＴｈｅｃｏｍｐｕｔａｔｉｏｎａｌｄｏｍａｉｎｗａｓｍｅｓｈｅｄｂｙＤｅｌａｕｎａｙｔｒｉａｎｇｕｌａｔｉｏｎｍｅｔｈｏｄ．ＡＭｏｎｔｅＣａｒｌｏｄｉｒｅｃｔｓｉｍｕｌａｔｉｏｎｐｒｏｇｒａｍｗａｓ ｄｅｖｅｌｏｐｅｄｆｏｒｕｎｓｔｒｕｃｔｕｒｅｄｇｒｉｄｓａｎｄｉｔｓｖａｌｉｄｉｔｙｗａｓｃｈｅｃｋｅｄ．Ｏｎｔｈａｔｃａｓｅｔｈｅｔｗｏ-ｄｉｍｅｎｓｉｏｎａｌｆｌｏｗａｎｄｈｅａｔｔｒａｎｓｆｅｒｏｆｒａｒｅｆｉｅｄｇａｓ ｉｎｐａｒａｌｌｅｌｐｌａｔｅｍｉｃｒｏ-ｃｈａｎｎｅｌｓｗｉｔｈｔｒｉａｎｇｕｌａｒｒｏｕｇｈｎｅｓｓｅｌｅｍｅｎｔｓｗｅｒｅｓｉｍｕｌａｔｅｄ．Ｔｈｅｐｒｅｓｓｕｒｅｏｆｔｈｅｍｉｃｒｏ-ｃｈａｎｎｅｌｓ’ｉｎｌｅｔａｎｄｏｕｔｌｅｔ ａｎｄｔｈｅｔｅｍｐｅｒａｔｕｒｅｏｆｔｈｅｐｌａｔｅｓａｂｏｖｅａｎｄｂｅｌｏｗｗｅｒｅｃｏｎｓｔａｎｔ．Ｔｈｅｉｎｆｌｕｅｎｃｅｓｏｆｔｈｅｈｅｉｇｈｔｗｉｄｔｈａｎｄｄｉｓｔｒｉｂｕｔｉｏｎｄｅｎｓｉｔｙｏｆｒｏｕｇｈ- ｎｅｓｓｅｌｅｍｅｎｔｓｗｅｒｅｃｏｍｐｕｔｅｄａｎｄａｎａｌｙｚｅｄ．Ｔｈｅｒｅｓｕｌｔｓｓｈｏｗｔｈａｔｒｏｕｇｈｎｅｓｓｅｌｅｍｅｎｔｓｉｎｔｈｅｍｉｃｒｏ-ｃｈａｎｎｅｌｓｃａｎｃｒｅａｔｅｏｂｖｉｏｕｓｄｉｓｔｕｒｂ- ａｎｃｅｆｏｒｔｈｅｆｌｏｗａｎｄｈｅａｔｔｒａｎｓｆｅｒ．Ｔｈｅｇｒｅａｔｅｒｔｈｅｒｏｕｇｈｎｅｓｓｅｌｅｍｅｎｔｓｔｈｅｍｏｒｅｖｉｓｉｂｌｅｔｈｅｏｓｃｉｌｌａｔｉｏｎｏｆｖｅｌｏｃｉｔｙｅｖｅｎｖｏｒｔｅｘｅｓｏｃ- ｃｕｒｗｈｉｃｈｉｎｃｒｅａｓｅｓｔｈｅｐｒｅｓｓｕｒｅｌｏｓｓ；ｂｕｔｒｏｕｇｈｎｅｓｓｅｌｅｍｅｎｔｓｃａｎｅｎｈａｎｃｅｔｈｅｈｅａｔｔｒａｎｓｆｅｒ． ＫＥＹＷＯＲＤＳ ｒｏｕｇｈｎｅｓｓ；ｍｉｃｒｏ-ｃｈａｎｎｅｌｓ；ｒａｒｅｆｉｅｄｇａｓ；ｆｌｏｗ；ｈｅａｔｔｒａｎｓｆｅｒ 收稿日期：2009--10--22 作者简介：林 林 （1968— ）男副教授博士Ｅ-ｍａｉｌ：ｌｉｎｌｉｎ＠ａｄｍｉｎ．ｕｓｔｂ．ｅｄｕ．ｃｎ 近年来微电子机械系统 （ｍｉｃｒｏｅｌｅｃｔｒｏｍｅｃｈａｎｉ- ｃａｌｓｙｓｔｅｍＭＥＭＳ）、人工超晶格、纳米材料和生物 芯片等技术的诞生和深入发展引起了人们对微细尺 度流动特性的关注．由于微电子系统的散热问题 微通道内稀薄气体的流动成为人们十分关心的问 题．微通道流动和换热与常规尺度的流动有较大的 区别从许多文献来看 ［1--5］一个很重要的原因来自 于壁面粗糙度对流动的影响．本文主要研究了微通 道内的粗糙元对于通道内气体的流动和换热的影 响．由于尺度的微细化及介质的稀薄效应常规的 连续性方程及 Ｎ--Ｓ方程等不再适用．近年来 ＤＳＭＣ （ｄｉｒｅｃｔｓｉｍｕｌａｔｉｏｎＭｏｎｔｅＣａｒｌｏ）方法已经被广 泛的应用于模拟稀薄气体动力学中的问题并取得 了非常好的模拟结果 ［6--12］．蒙特卡洛直接模拟方法 DOI :10．13374／j．issn1001—053x．2010．03．023
第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数大小范围，可以将气体流动 划分为四个不同区域：Knl0自由分子流动区. 图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·
,386 北京科技大学学报 第32卷 标准Bi程序及文献[15]中的解析解，可见基于非 结构网格的本文模拟结果（图3(b))与结构网格的 4微通道流动换热模拟结果及分析 标准程序（图3(a))的模拟结果吻合良好. 4.1流动特性 4.1.1速度滑移 30405060708090100110120130140V/%m.s 0.5 图5为工况135和7通道中水平方向速度云 图.由图5可以看出，各工况中微通道内的流场整 体趋势大致相同，不同之处在于，壁面处的速度滑 X/uim 移的大小 0.5 30405060708090100110120130140W%m.） 7ms 05 2030405060708090100110120130140150 X/um 图3光滑微通道水平方向速度图·(a)基于结构化网格的Bm 标准程序模拟结果；(b)基于非结构化网格的本文模拟结果 V(m-s) Fig 3 Horizontal velocity charts of sooth m icmo channels (a) 02 -10010 Bim-standand prognm smulation results based on stmuctured mes (b)siulation results based on unstructured mesh 图4所示的是水平距离1.5m处（即通道中 m ■ms 间处)截面上的速度分布，从图中可以看出，本文的 -60-50-40-30-20-100102030405060708090100 模拟结果与Bid标准程序的模拟结果以及文 献[15的解析解基本吻合 44 0.9 V(m-s) 0.8 -60-50-40-30-20-100102030405060708090100 0.7 0.6 ·基于非结构网格下的结算结果 日0.5 04 q基于结构化网格下hbid标准程序 Xum 的计算结果 0.3 一解析解 图5工况1(a,工况3(b).工况5(c)和工况7(d)微通道中 02 水平方向速度图 0.1 Fig5 Horizontal velocity contours ofm icro-flows in Cases 1 (a).3 (b).5 (c)and 7 (d) 20 30 40 5060708090100110 7(m+s 本文将所有工况的壁面滑移速度的平均值计算 图4光滑微通道截面速度分布(x=1.5“m) Fig 4 Canparison of velcity at them ithvay (x=1.5!m)of mooth 出来，列在表3中.从表3可看出，相对光滑微通 m icmo"channe ls 道，粗糙元的存在减小了壁面的平均滑移速度.粗 表3工况1~7中壁面滑移速度的平均值 Table 3 Average slip vebcity of the wall in Cases I to7 工况编号 2 3 6 7 壁面滑移速度平均值/(m·1) 27.667 8.3948 2.7071 2.4386 2.7020 2.2834 2.7343 糙元较小时，壁面滑移速度减小较明显；粗糙元增 4.1.2压力非线性 大，壁面滑移速度变化不明显， 图7给出了各工况微通道内的压力分布图，从 图6所示为各工况局部气体速度矢量图，从图 图中可以看出，粗糙元对压力分布的总体趋势影响 中可以看出，粗糙元的存在使气体产生漩涡的趋势， 不大，但粗糙元会对通道的局部压力产生较为明显 随着粗糙元增大和加密，粗糙元间的漩涡更加明显， 的影响，这主要表现在，在粗糙元的迎风面，局部压 漩涡的出现，必将增加微通道内的压力损失， 力偏大，而在粗糙元的背风面，局部压力偏小，亦即
北 京 科 技 大 学 学 报 第 32卷 标准 Ｂｉｒｄ程序及文献 ［15］中的解析解可见基于非 结构网格的本文模拟结果 （图 3（ｂ））与结构网格的 标准程序 （图 3（ａ））的模拟结果吻合良好． 图 3 光滑微通道水平方向速度图 ∙（ａ）基于结构化网格的 Ｂｉｒｄ 标准程序模拟结果；（ｂ） 基于非结构化网格的本文模拟结果 Ｆｉｇ．3 Ｈｏｒｉｚｏｎｔａｌｖｅｌｏｃｉｔｙｃｈａｒｔｓｏｆｓｍｏｏｔｈｍｉｃｒｏ-ｃｈａｎｎｅｌｓ： （ａ） Ｂｉｒｄ-ｓｔａｎｄａｒｄｐｒｏｇｒａｍ ｓｉｍｕｌａｔｉｏｎｒｅｓｕｌｔｓｂａｓｅｄｏｎｓｔｒｕｃｔｕｒｅｄｍｅｓｈ； （ｂ） ｓｉｍｕｌａｔｉｏｎｒｅｓｕｌｔｓｂａｓｅｄｏｎｕｎｓｔｒｕｃｔｕｒｅｄｍｅｓｈ 图 4所示的是水平距离 1∙5μｍ处 （即通道中 间处 ）截面上的速度分布．从图中可以看出本文的 模拟结 果 与 Ｂｉｒｄ标 准 程 序 的 模 拟 结 果 以 及 文 献 ［15］的解析解基本吻合． 图 4 光滑微通道截面速度分布 （ｘ＝1∙5μｍ） Ｆｉｇ．4 Ｃｏｍｐａｒｉｓｏｎｏｆｖｅｌｏｃｉｔｙａｔｔｈｅｍｉｄｗａｙ（ｘ＝1∙5μｍ）ｏｆｓｍｏｏｔｈ ｍｉｃｒｏ-ｃｈａｎｎｅｌｓ 4 微通道流动换热模拟结果及分析 4∙1 流动特性 4∙1∙1 速度滑移 图 5为工况 1、3、5和 7通道中水平方向速度云 图．由图 5可以看出各工况中微通道内的流场整 体趋势大致相同．不同之处在于壁面处的速度滑 移的大小． 图 5 工况 1（ａ）、工况 3（ｂ）、工况 5（ｃ）和工况 7（ｄ）微通道中 水平方向速度图 Ｆｉｇ．5 Ｈｏｒｉｚｏｎｔａｌｖｅｌｏｃｉｔｙｃｏｎｔｏｕｒｓｏｆｍｉｃｒｏ-ｆｌｏｗｓｉｎＣａｓｅｓ1（ａ）3 （ｂ）5（ｃ） ａｎｄ7（ｄ） 本文将所有工况的壁面滑移速度的平均值计算 出来列在表 3中．从表 3可看出相对光滑微通 道粗糙元的存在减小了壁面的平均滑移速度．粗 表 3 工况 1～7中壁面滑移速度的平均值 Ｔａｂｌｅ3 ＡｖｅｒａｇｅｓｌｉｐｖｅｌｏｃｉｔｙｏｆｔｈｅｗａｌｌｉｎＣａｓｅｓ1ｔｏ7 工况编号 1 2 3 4 5 6 7 壁面滑移速度平均值／（ｍ·ｓ—1） 27∙667 8∙3948 2∙7071 2∙4386 2∙7020 2∙2834 2∙7343 糙元较小时壁面滑移速度减小较明显；粗糙元增 大壁面滑移速度变化不明显． 图 6所示为各工况局部气体速度矢量图．从图 中可以看出粗糙元的存在使气体产生漩涡的趋势 随着粗糙元增大和加密粗糙元间的漩涡更加明显． 漩涡的出现必将增加微通道内的压力损失． 4∙1∙2 压力非线性 图 7给出了各工况微通道内的压力分布图．从 图中可以看出粗糙元对压力分布的总体趋势影响 不大但粗糙元会对通道的局部压力产生较为明显 的影响．这主要表现在在粗糙元的迎风面局部压 力偏大而在粗糙元的背风面局部压力偏小亦即 ·386·
第3期 林林等：粗糙微通道内稀薄气体流动与换热的蒙特卡洛直接模拟 ,387. 0.05 110 100f 一压力 …辅助助线 2 2.2 2.3 90 00 80 2.1 2.2 2.3 704 XAIm 60 0.05 05 1.0 152.02.5 3.0 um 2.1 2.2 23 图8工况1通道内压力分布图 XAuim @ Fig8 Pressure distrbution along the cen terline ofm icro channels in Case l 图6工况2(a),工况4(b)和工况6(c)微通道局部气体速度 矢量图 所致；在出口处温度的下降是由于出口的压力梯度 Fig 6 Local velocity distribution of gas n m icm channels of Cases 2 变大，压力急剧降低所致，随着粗糙元的变大和增 (a),4(b)and6(c) 多，通道内的气体的温度进一步升高，分析其原因， 在粗糙元附近压力的跳跃比较明显，分析其原因， 是因为粗糙元变大和增多以后，气体的壁面滑移速 主要因为气体在粗糙元的迎风面遇到障碍物”，发 度减小，气体在通道内停留时间增长，气体被加热到 生了一定的滞止，引起了压力的升高；而在粗糙元的 了更高的温度 背风面，产生了相应的空隙，由此造成了背风面压力 a■ K 300305310315320325330335340345350 的减小 120 2 110 XAim 100 ■ 国K 300305310315320325330335340345350 0 80 X小m 60 一工况1 (■ K 一工况2 0.5 300305310315320325330335340345350 50 05 1.0 15 2.0 2.5 3.0 X/um XAim 图7工况12通道内粗糙元附近的压力分布图 d 300305310315320325330335340345350 Fig 7 Pressure distrbution near mugh elments abng the flow diree- tion in Case 1 and 2 2 图8给出的是工况1通道内的压力分布图.通 X/um 过图中可以看出，微通道内的压力呈非线性分布，在 图9工况1(a)、工况2(b)、工况4(c)和工况7(d)中微通道 出口处的压力梯度明显高于通道的其他部位，这一 内温度图 结论已经被前人采用实验和数值模拟的方法 Fig9 Temperature contours ofm icmo-flows n Cases 1 (a).2 (b). 4 (c)and 7 (d) 分别得以证实 4.2传热特性 图10所示为各工况中通道中心线处的温度分 4.2.1温度跳跃 布、由图中可以看出，通道中心线处的温度亦是先 图9所示的为工况12、4和7微通道内气体的 升高后降低。沿着气体流动的方向，温度梯度在入 温度图.从图中可以看出，气体的温度都呈先升高 口段较大，这是由于此时气体与壁面温差较大，气体 后降低的趋势，先升高是由于气体被壁面剧烈加热 与壁面换热剧烈；在通道中间段由于气体与壁面温
第 3期 林 林等： 粗糙微通道内稀薄气体流动与换热的蒙特卡洛直接模拟 图 6 工况 2（ａ）、工况 4（ｂ）和工况 6（ｃ）微通道局部气体速度 矢量图 Ｆｉｇ．6 Ｌｏｃａｌｖｅｌｏｃｉｔｙｄｉｓｔｒｉｂｕｔｉｏｎｏｆｇａｓｉｎｍｉｃｒｏ-ｃｈａｎｎｅｌｓｏｆＣａｓｅｓ2 （ａ）4（ｂ） ａｎｄ6（ｃ） 在粗糙元附近压力的跳跃比较明显．分析其原因 主要因为气体在粗糙元的迎风面遇到 “障碍物 ”发 生了一定的滞止引起了压力的升高；而在粗糙元的 背风面产生了相应的空隙由此造成了背风面压力 的减小． 图 7 工况 1、2通道内粗糙元附近的压力分布图 Ｆｉｇ．7 Ｐｒｅｓｓｕｒｅｄｉｓｔｒｉｂｕｔｉｏｎｎｅａｒｒｏｕｇｈｅｌｅｍｅｎｔｓａｌｏｎｇｔｈｅｆｌｏｗｄｉｒｅｃ- ｔｉｏｎｉｎＣａｓｅ1ａｎｄ2 图8给出的是工况1通道内的压力分布图．通 过图中可以看出微通道内的压力呈非线性分布在 出口处的压力梯度明显高于通道的其他部位．这一 结论已经被前人采用实验 ［16］和数值模拟 ［17］的方法 分别得以证实． 4∙2 传热特性 4∙2∙1 温度跳跃 图 9所示的为工况 1、2、4和 7微通道内气体的 温度图．从图中可以看出气体的温度都呈先升高 后降低的趋势．先升高是由于气体被壁面剧烈加热 图 8 工况 1通道内压力分布图 Ｆｉｇ．8 Ｐｒｅｓｓｕｒｅｄｉｓｔｒｉｂｕｔｉｏｎａｌｏｎｇｔｈｅｃｅｎｔｅｒｌｉｎｅｏｆｍｉｃｒｏ-ｃｈａｎｎｅｌｓｉｎ Ｃａｓｅ1 所致；在出口处温度的下降是由于出口的压力梯度 变大压力急剧降低所致．随着粗糙元的变大和增 多通道内的气体的温度进一步升高．分析其原因 是因为粗糙元变大和增多以后气体的壁面滑移速 度减小气体在通道内停留时间增长气体被加热到 了更高的温度． 图 9 工况 1（ａ）、工况 2（ｂ）、工况 4（ｃ）和工况 7（ｄ）中微通道 内温度图 Ｆｉｇ．9 Ｔｅｍｐｅｒａｔｕｒｅｃｏｎｔｏｕｒｓｏｆｍｉｃｒｏ-ｆｌｏｗｓｉｎＣａｓｅｓ1（ａ）2（ｂ） 4（ｃ） ａｎｄ7（ｄ） 图 10所示为各工况中通道中心线处的温度分 布．由图中可以看出通道中心线处的温度亦是先 升高后降低．沿着气体流动的方向温度梯度在入 口段较大这是由于此时气体与壁面温差较大气体 与壁面换热剧烈；在通道中间段由于气体与壁面温 ·387·
,388 北京科技大学学报 第32卷 差减小，气体与壁面间的换热减弱，故而温度变化不 面之间的热流密度在出口处回升 明显；在通道出口处由于压降较大，温度又发生了较 通过与光滑微通道（工况1）比较，还能明显地 明显的下降，粗糙元越大，分布越密，中心线处的温 看出，粗糙微通道壁面热流密度的跳跃十分明显，且 度也就越高，这说明，粗糙元的存在可以使通道内 热流密度小的位置对应于粗糙元的表面，而热流密 的气体达到更高的温度，并且，粗糙元越多、越大， 度较大的位置对应于粗糙元之间的通道壁面， 气体所能达到的温度也就越高，通过图中还可以得 为了方便比较，本文计算了各工况微通道壁面 到，粗糙元较小、较少时，相对于光滑微通道，粗糙元 的热流量和质量流量，见表4 对通道内气体温度跳跃的影响越大（工况2、3），随 表4通道内平均热流量与质量流量 着粗糙元的增大增多，通道内的气体温度虽然进一 Table 4 Average heat transfer rate ofm icmo-fows and mass fhxes 步升高，但温度跳跃量增加不大（工况56和7） 工况 12345 6 7 345 热流量13.1813.0612.7212.3211.3412.0210.77 340 质量流量kg)3.983.573.392.882.792.572.63 335 从表中可以看出：随着粗糙元增大和加密，微通 330 道内的平均热流量和质量流量均呈下降趋势；质量 ◆一工况1 *-工况2 流量的减小，主要是由于粗糙元的阻力效应引起的， 0 +一工况3 --工况4 随着粗糙元的增大和增多，通道的壁面热流虽 315 -。一工况5 然减小，但是由于通道内的质量流量也相应减少，单 310 0-工况6 +-工况7 位质量的气体吸收热量反而增多，这就解释了随着 305 粗糙元的增大和增多通道内的气体温度升高的 3006 0.5 1.0 1.5 2.0 2.5 3.0 原因 X/um 图10各工况微通道中心线处的温度分布 5结论 Fig 10 Change in temperature alng the centerline of m icmchan nels n Cases l to 7 (1)粗糙元对微通道内的流动、换热扰动显著： 随着粗糙元尺寸的增加和个数的增加，扰动更为 4.2.2热流密度 剧烈， 图11所示的为微通道的壁面热流密度，由图 (2)粗糙微通道内的壁面速度滑移小于光滑微 中可以看出，沿着气体流动的方向，壁面的热流密度 通道；随着粗糙元的增大和增多，速度出现更为严重 总体呈下降趋势，在出口处有一小段上升.这是因 的跳跃，甚至出现漩涡，增加了通道内的压力损失， 为，沿着气体流动的方向由于气体一直被壁面加热， (3)微通道的进、出口处温度变化明显：进口处 气体的温度不新升高，与壁面的温差减小，因此气体 温度剧烈升高，出口处温度下降.随着粗糙元尺寸 与壁面之间的热流密度逐渐减小；在出口附近，气体 和个数的增加，通道内气体的温度升高， 的温度下降，与壁面之间的温差变大，因此气体与壁 (4)粗糙元的存在导致微通道的质量流量和壁 工况1 面平均热流密度减小. ·工况2 +工况3 。通道壁面 参考文献 [1]Zhou J Prediction of F low and Heat Transfer in M icm Scale Chan nel w ith Dimet Smulation Monte Carlo Method Dissertation ] ,粗糙元表面 Xi'an:Xi'an Jiaolong University 2004 435158)13美4511年9 (周靖，用直接模拟蒙特卡洛法计算微通道中的流动与换热 [学位论文1西安：西安交通大学，2004) 0.5 1.01520253.0 [2]Wang H L W ang Y.Liu J Pertubation analysis of wall mough- X/um ness effect on flow in m icmo-channel beteen wo parallel plates J 图11工况1、2和6微通道壁面热流密度 Xian Jiaotong Univ 2005 39(1):540 Fig 11 Heat fhx at the surface abng the flw direction in Cases l. (王吴利，王元，刘江平板微通道壁面粗糙度对流场影响的摄 2 and 6 动分析.西安交通大学学报，200539(5)：540)
北 京 科 技 大 学 学 报 第 32卷 差减小气体与壁面间的换热减弱故而温度变化不 明显；在通道出口处由于压降较大温度又发生了较 明显的下降．粗糙元越大分布越密中心线处的温 度也就越高．这说明粗糙元的存在可以使通道内 的气体达到更高的温度．并且粗糙元越多、越大 气体所能达到的温度也就越高．通过图中还可以得 到粗糙元较小、较少时相对于光滑微通道粗糙元 对通道内气体温度跳跃的影响越大 （工况 2、3）随 着粗糙元的增大增多通道内的气体温度虽然进一 步升高但温度跳跃量增加不大 （工况 5、6和 7）． 图 10 各工况微通道中心线处的温度分布 Ｆｉｇ．10 Ｃｈａｎｇｅｉｎｔｅｍｐｅｒａｔｕｒｅａｌｏｎｇｔｈｅｃｅｎｔｅｒｌｉｎｅｏｆｍｉｃｒｏ-ｃｈａｎ- ｎｅｌｓｉｎＣａｓｅｓ1ｔｏ7 4∙2∙2 热流密度 图 11 工况 1、2和 6微通道壁面热流密度 Ｆｉｇ．11 ＨｅａｔｆｌｕｘａｔｔｈｅｓｕｒｆａｃｅａｌｏｎｇｔｈｅｆｌｏｗｄｉｒｅｃｔｉｏｎｉｎＣａｓｅｓ1 2ａｎｄ6 图 11所示的为微通道的壁面热流密度．由图 中可以看出沿着气体流动的方向壁面的热流密度 总体呈下降趋势在出口处有一小段上升．这是因 为沿着气体流动的方向由于气体一直被壁面加热 气体的温度不断升高与壁面的温差减小因此气体 与壁面之间的热流密度逐渐减小；在出口附近气体 的温度下降与壁面之间的温差变大因此气体与壁 面之间的热流密度在出口处回升． 通过与光滑微通道 （工况 1）比较还能明显地 看出粗糙微通道壁面热流密度的跳跃十分明显且 热流密度小的位置对应于粗糙元的表面而热流密 度较大的位置对应于粗糙元之间的通道壁面． 为了方便比较本文计算了各工况微通道壁面 的热流量和质量流量见表 4． 表 4 通道内平均热流量与质量流量 Ｔａｂｌｅ4 Ａｖｅｒａｇｅｈｅａｔｔｒａｎｓｆｅｒｒａｔｅｏｆｍｉｃｒｏ-ｆｌｏｗｓａｎｄｍａｓｓｆｌｕｘｅｓ 工况 1 2 3 4 5 6 7 热流量／Ｗ 13∙1813∙0612∙7212∙3211∙3412∙0210∙77 质量流量／（ｋｇ·ｓ —1 ） 3∙98 3∙57 3∙39 2∙88 2∙79 2∙57 2∙63 从表中可以看出：随着粗糙元增大和加密微通 道内的平均热流量和质量流量均呈下降趋势；质量 流量的减小主要是由于粗糙元的阻力效应引起的． 随着粗糙元的增大和增多通道的壁面热流虽 然减小但是由于通道内的质量流量也相应减少单 位质量的气体吸收热量反而增多这就解释了随着 粗糙元的增大和增多通道内的气体温度升高的 原因． 5 结论 （1）粗糙元对微通道内的流动、换热扰动显著； 随着粗糙元尺寸的增加和个数的增加扰动更为 剧烈． （2）粗糙微通道内的壁面速度滑移小于光滑微 通道；随着粗糙元的增大和增多速度出现更为严重 的跳跃甚至出现漩涡增加了通道内的压力损失． （3）微通道的进、出口处温度变化明显：进口处 温度剧烈升高出口处温度下降．随着粗糙元尺寸 和个数的增加通道内气体的温度升高． （4）粗糙元的存在导致微通道的质量流量和壁 面平均热流密度减小． 参 考 文 献 ［1］ ＺｈｏｕＪ．ＰｒｅｄｉｃｔｉｏｎｏｆＦｌｏｗａｎｄＨｅａｔＴｒａｎｓｆｅｒｉｎＭｉｃｒｏ-ＳｃａｌｅＣｈａｎ- ｎｅｌｗｉｔｈＤｉｒｅｃｔＳｉｍｕｌａｔｉｏｎＭｏｎｔｅＣａｒｌｏＭｅｔｈｏｄ ［Ｄｉｓｓｅｒｔａｔｉｏｎ］． Ｘｉ’ａｎ：Ｘｉ’ａｎＪｉａｏｔｏｎｇＵｎｉｖｅｒｓｉｔｙ2004 （周靖．用直接模拟蒙特卡洛法计算微通道中的流动与换热 ［学位论文 ］．西安：西安交通大学2004） ［2］ ＷａｎｇＨＬＷａｎｇＹＬｉｕＪ．Ｐｅｒｔｕｒｂａｔｉｏｎａｎａｌｙｓｉｓｏｆｗａｌｌｒｏｕｇｈ- ｎｅｓｓｅｆｆｅｃｔｏｎｆｌｏｗｉｎｍｉｃｒｏ-ｃｈａｎｎｅｌｂｅｔｗｅｅｎｔｗｏｐａｒａｌｌｅｌｐｌａｔｅｓ．Ｊ ＸｉʾａｎＪｉａｏｔｏｎｇＵｎｉｖ200539（1）：540 （王昊利王元刘江．平板微通道壁面粗糙度对流场影响的摄 动分析．西安交通大学学报200539（5）：540） ·388·
第3期 林林等：粗糙微通道内稀薄气体流动与换热的蒙特卡洛直接模拟 ,389. [3]Cmoce G.DAgam P.Numerical analysis of mughness effect on Mi Comput Fhids2001.30(6):711 cmotube heat transfer Superla ttices M icmostnict 2004.35(3-6): [11]Qin F H.Sun D J Yn X Y.DSMC for gas flows n a micm- 601 channel Chin J Camput Phys 2001.18(6):507 [4 Hu Y D.Wemer C Li D Q Infhence of three dmensional (秦丰华，孙德军，尹协远，微管道气体流动的蒙特卡洛直接 roughness on pressuredriven flow through m icrichannels Trans 模拟.计算物理，2001,18(6)：507) ASME2003125(2):871 [12]Bin G A.Mokculr Gas Dynam ics and the Direct Sinulation of [5]Andrade J JS Henrique E A A.Heat transport thmugh mugh Gas F kw London Oxfo University Press 1994 channels PhysA 2004.339(34):296 [13]Bin G A.Molecular Gas Dynam ics London Oxford University [6]Le M.Hassan I EsailN The effects of outlet boundary condi Prss 1976 tions on siulating supersonic m icrochannel flws using DSMC [14 Yang X H.Zhu W L The innuence of inner surface roughness Appl Them Eng 2007.27(1):21 on the fiction factor of gas transn ission pipeline Petmo Chan [7]Liu H W,W ang M.Monte Carlo smnulations of gas fow and heat Eqmp200534(1):25 transfer in vacuum packaged MEMS devices Appl Them Eng (杨晓鸿，朱薇玲.粗糙度对输气管道摩阻系数的影响·石油 2007,27(23):323 化工设备，2005,34(1)：25) [8]W ang M R.LiZ X.Smultions for gas flows in m icmogeametries using the direct siulation Monte Carlo method Int J Heat Flid [15]Arkilic E B Schm itM A.Gascous slip flow n long m icmochan- Flw,2004,25(6):975 nels JM icroelectrmech Syst 1997,6(2):167 [9]W ang M R.LiZ X.Gas m ixing in m icmchannels using the direet [16]Ho C M.Tai Y C Micm electricalmechanical systems sinulation Monte Carlo method Int J Heat Mass Transfer 2006, (MEMS)and fhid flows Annu Rev Mech 1998 30.579 49(9h0):1696 [17]Tang G H.Tao W Q.He Y L Lattice Boltaann method for [10]Wu JS Tseng K C Analysis ofm icmo scale gas fows with pres sinulation gas flow in m icmchannels Int J Mod Phys C 2004. sure boundaries using direct simulation Monte Carlo method 15(3):335
第 3期 林 林等： 粗糙微通道内稀薄气体流动与换热的蒙特卡洛直接模拟 ［3］ ＣｒｏｃｅＧＤʾＡｇａｒｏＰ．ＮｕｍｅｒｉｃａｌａｎａｌｙｓｉｓｏｆｒｏｕｇｈｎｅｓｓｅｆｆｅｃｔｏｎＭｉ- ｃｒｏｔｕｂｅｈｅａｔｔｒａｎｓｆｅｒ．ＳｕｐｅｒｌａｔｔｉｃｅｓＭｉｃｒｏｓｔｒｕｃｔ200435（3—6）： 601 ［4］ ＨｕＹ ＤＷｅｒｎｅｒＣＬｉＤ Ｑ．Ｉｎｆｌｕｅｎｃｅｏｆｔｈｒｅｅ-ｄｉｍｅｎｓｉｏｎａｌ ｒｏｕｇｈｎｅｓｓｏｎｐｒｅｓｓｕｒｅ-ｄｒｉｖｅｎｆｌｏｗ ｔｈｒｏｕｇｈｍｉｃｒｉｃｈａｎｎｅｌｓ．Ｔｒａｎｓ ＡＳＭＥ2003125（2）：871 ［5］ ＡｎｄｒａｄｅＪＪＳＨｅｎｒｉｑｕｅＥＡＡ．Ｈｅａｔｔｒａｎｓｐｏｒｔｔｈｒｏｕｇｈｒｏｕｇｈ ｃｈａｎｎｅｌｓ．ＰｈｙｓＡ2004339（3／4）：296 ［6］ ＬｅＭＨａｓｓａｎＩＥｓｍａｉｌＮ．Ｔｈｅｅｆｆｅｃｔｓｏｆｏｕｔｌｅｔｂｏｕｎｄａｒｙｃｏｎｄｉ- ｔｉｏｎｓｏｎｓｉｍｕｌａｔｉｎｇｓｕｐｅｒｓｏｎｉｃｍｉｃｒｏｃｈａｎｎｅｌｆｌｏｗｓｕｓｉｎｇＤＳＭＣ． ＡｐｐｌＴｈｅｒｍＥｎｇ200727（1）：21 ［7］ ＬｉｕＨＷＷａｎｇＭ．ＭｏｎｔｅＣａｒｌｏｓｉｍｕｌａｔｉｏｎｓｏｆｇａｓｆｌｏｗａｎｄｈｅａｔ ｔｒａｎｓｆｅｒｉｎｖａｃｕｕｍ ｐａｃｋａｇｅｄＭＥＭＳｄｅｖｉｃｅｓ．ＡｐｐｌＴｈｅｒｍ Ｅｎｇ 200727（2／3）：323 ［8］ ＷａｎｇＭＲＬｉＺＸ．Ｓｉｍｕｌａｔｉｏｎｓｆｏｒｇａｓｆｌｏｗｓｉｎｍｉｃｒｏｇｅｏｍｅｔｒｉｅｓ ｕｓｉｎｇｔｈｅｄｉｒｅｃｔｓｉｍｕｌａｔｉｏｎＭｏｎｔｅＣａｒｌｏｍｅｔｈｏｄ．ＩｎｔＪＨｅａｔＦｌｕｉｄ Ｆｌｏｗ200425（6）：975 ［9］ ＷａｎｇＭＲＬｉＺＸ．Ｇａｓｍｉｘｉｎｇｉｎｍｉｃｒｏｃｈａｎｎｅｌｓｕｓｉｎｇｔｈｅｄｉｒｅｃｔ ｓｉｍｕｌａｔｉｏｎＭｏｎｔｅＣａｒｌｏｍｅｔｈｏｄ．ＩｎｔＪＨｅａｔＭａｓｓＴｒａｎｓｆｅｒ2006 49（9／10）：1696 ［10］ ＷｕＪＳＴｓｅｎｇＫＣ．Ａｎａｌｙｓｉｓｏｆｍｉｃｒｏ-ｓｃａｌｅｇａｓｆｌｏｗｓｗｉｔｈｐｒｅｓ- ｓｕｒｅｂｏｕｎｄａｒｉｅｓｕｓｉｎｇｄｉｒｅｃｔｓｉｍｕｌａｔｉｏｎＭｏｎｔｅＣａｒｌｏｍｅｔｈｏｄ． ＣｏｍｐｕｔＦｌｕｉｄｓ200130（6）：711 ［11］ ＱｉｎＦＨＳｕｎＤＪＹｉｎＸＹ．ＤＳＭＣｆｏｒｇａｓｆｌｏｗｓｉｎａｍｉｃｒｏ- ｃｈａｎｎｅｌ．ＣｈｉｎＪＣｏｍｐｕｔＰｈｙｓ200118（6）：507 （秦丰华孙德军尹协远．微管道气体流动的蒙特卡洛直接 模拟．计算物理200118（6）：507） ［12］ ＢｉｒｄＧＡ．ＭｏｌｅｃｕｌａｒＧａｓＤｙｎａｍｉｃｓａｎｄｔｈｅＤｉｒｅｃｔＳｉｍｕｌａｔｉｏｎｏｆ ＧａｓＦｌｏｗ．Ｌｏｎｄｏｎ：ＯｘｆｏｒｄＵｎｉｖｅｒｓｉｔｙＰｒｅｓｓ1994 ［13］ ＢｉｒｄＧＡ．ＭｏｌｅｃｕｌａｒＧａｓＤｙｎａｍｉｃｓ．Ｌｏｎｄｏｎ：ＯｘｆｏｒｄＵｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ1976 ［14］ ＹａｎｇＸＨＺｈｕＷ Ｌ．Ｔｈｅｉｎｎｕｅｎｃｅｏｆｉｎｎｅｒｓｕｒｆａｃｅｒｏｕｇｈｎｅｓｓ ｏｎｔｈｅｆｉｃｔｉｏｎｆａｃｔｏｒｏｆｇａｓｔｒａｎｓｍｉｓｓｉｏｎｐｉｐｅＩｉｎｅ．ＰｅｔｒｏＣｈｅｍ Ｅｑｕｉｐ200534（1）：25 （杨晓鸿朱薇玲．粗糙度对输气管道摩阻系数的影响．石油 化工设备200534（1）：25） ［15］ ＡｒｋｉｌｉｃＥＢＳｃｈｍｉｄｔＭＡ．Ｇａｓｅｏｕｓｓｌｉｐｆｌｏｗｉｎｌｏｎｇｍｉｃｒｏｃｈａｎ- ｎｅｌｓ．ＪＭｉｃｒｏｅｌｅｃｔｒｏｍｅｃｈＳｙｓｔ19976（2）：167 ［16］ Ｈｏ Ｃ Ｍ Ｔａｉ Ｙ Ｃ． Ｍｉｃｒｏ-ｅｌｅｃｔｒｉｃａｌ-ｍｅｃｈａｎｉｃａｌ-ｓｙｓｔｅｍｓ （ＭＥＭＳ） ａｎｄｆｌｕｉｄｆｌｏｗｓ．ＡｎｎｕＲｅｖＭｅｃｈ199830：579 ［17］ ＴａｎｇＧＨＴａｏＷ ＱＨｅＹＬ．ＬａｔｔｉｃｅＢｏｌｔｚｍａｎｎｍｅｔｈｏｄｆｏｒ ｓｉｍｕｌａｔｉｏｎｇａｓｆｌｏｗｉｎｍｉｃｒｏｃｈａｎｎｅｌｓ．ＩｎｔＪＭｏｄＰｈｙｓＣ2004 15（3）：335 ·389·
