本发明涉及流体力学领域，是一种复杂几何边界下微流体器件中气体滑移流动的计算方法，其步骤：1)基于有限体积方法，采用网格划分软件，对计算域和计算边界进行离散；2)基于Maxwell分子模型和Knudsen层内质量、动量、能量守恒，得到适用于复杂几何表面的速度滑移模型；3)利用该模型，结合流体计算商业软件ANSYS FLUENT，编写用户自定义函数，计算稀薄气体在复杂几何表面的滑移速度。本发明几何边界适应性强，能与现有的商业软件结合，使用方便。
1.一种复杂几何边界下微流体器件中气体滑移流动的计算方法，包括以下步骤：
1)对计算区域进行离散，其中计算区域可以是任意形状的，在有限体积法的离散过程中，复杂的几何表面被离散成了一定数量的小光滑表面，设原几何表面为S，每个小光滑表面记为Sk,则S＝USk，对于每个Sk，设它的单位法向量为 入射的气体分子与壁面发生碰撞，根据Maxwell模型：式中，Us是壁面处的滑移速度，σt是切向动量适应系数(TMAC)，λ代表气体的平均分子自由程，dU/dz|w是壁面处法向速度梯度；一部分分子σt被固体表面吸收后重新喷出，喷出的角度是完全随机的，该部分称为散射分子；而另一部分分子1-σt则以镜面反射的形式从壁面返回气体流动区域，该部分称为折射分子，因此，在靠近壁面的气体薄层内Knudsen层，厚度与平均自由程同级，存在三种分子：入射分子、散射分子、折射分子，这三种分子向固体表面传递动量，而气体的滑移速度可以通过传递动量的守恒特性推导得出；
2)设整个计算区域的笛卡尔坐标系为(x,y,z)，在小光滑表面Sk处建立局部坐标系(xk,yk,zk)，其中，单位法向量 作为yk轴，另外两轴根据右手定则任意选取，设气体在总体坐标系中的速度记为(u,v,w)，在局部坐标系中则为(uk,vk,wk)，在局部坐标系中，Knudsen层内分子向壁面传递动量的总通量为：其中，(Uk,Vk,Wk)是分子热运动的速度， 代表单个分子传递的动量，f(Uk,Vk,Wk)是分子的分布函数，函数 可以表示为 或者式中m代表单个分子的质量，(usk,vsk,wsk)是局部坐标系中的滑移速度，总通量P由三部分组成：P＝Pi+(1-σt)Pr+σtPe (0.3)
式中，Pi、Pr和Pe分别是由入射分子、折射分子以及散射分子携带的动量，气体分子在散射过程中，散射角度是完全随机的，因此散射分子携带的动量能互相抵消，Pe＝0；而在折射过程中，切向动量不变，法向方向相反，因此有Pr＝-Pi，公式(0.3)可以简化为：P＝σtPi (0.4)
公式(0.2)和(0.5)中的分布函数f满足Boltzmann方程，这里采用由Chapman和Cowling给出的f关于Kn数的一阶近表达式：其中，n是分子的数密度，k是Boltzmann常数,kc是热传导系数, R是气体常数，(C1,C2,C3)＝(Uk,Vk,Wk),(x1,x2,x3)＝(xk,yk,zk),(c1,c2,c3)＝(uk,vk,wk)，
3)把公式(0.6)带入(0.2)和(0.5)，并对Uk、Vk和Wk进行积分，可得：把公式(0.7)、(0.8)带入(0.4)可得：
公式中，与Maxwell滑移模型相比，分子上的 也可以根据壁面处的切向剪切应力得到，分母上多出的项代表yk方向的剪切应力与静压强的比值，可以被忽略，公式(0.10)给出的滑移速度usk是局部坐标系下的结果，总体坐标系下的滑移速度(us,vs,ws)可以通过坐标变换的原理得到：公式中， vk＝lku+mkv+nkw，对于等温流体，公式中与温度相
在有限体积法中，以公式(0.11)代替无滑移速度边界条件，即得到了适用于复杂几何边界的气体滑移流动的计算方法。
复杂几何边界下微流体器件中气体滑移流动的计算方法
[技术领域]
[0001] 本发明涉及流体力学的技术领域，具体地说是一种复杂几何边界下微流体器件中气体滑移流动的计算方法。[背景技术]
[0002] 最近几十年来，MEMS(Micro-Electro-Mechanical Systems)技术取得了迅速的发展，MEMS产品已经遍布航天、医疗、汽车制造等众多领域，对人们的生活产生了巨大影响。其中，微流体器件是MEMS器件的重要组成部分，在提取生物、航天、环境监测、化学等领域都有着广泛的应用，例如，微通道可以分离细胞、冷却集成电路，微喷头广泛用于喷墨打印机，微型收缩-扩张喷管用于小型卫星的微型推进系统。因此，微流体的研究渐渐成为MEMS理论中的一个重要分支。由于尺寸的减小，微流体往往不满足连续性假设，因此需要特殊的建模方法。
[0003] 气体是由大量气体分子组成的系统，定义气体的平均分子自由程与特征长度的比值为Kn数，它是表征气体连续性的无量纲参数。不同Kn数下的气体，有不同的建模方法：
[0004] 1)当Kn数小于0.001时，气体处于连续流区。在该流区，气体平均自由程与特征长度相比可以忽略，采用传统的连续模型对气体进行建模。在连续模型中，认为气体是由无穷多的气体质点组成的，气体质点一方面在宏观上足够小，保证连续性；另一方面在微观上足够大，含有大量无序运动的气体分子从而表现出稳定的宏观特性。对于牛顿流体，气体的控制方程是著名的Navier-Stokes(N-S)方程组，它在科学研究和工程实际中取得了巨大成功。
[0005] 2)当Kn数大于0.001而小于0.1时，气体处于滑移流区。气体的平均自由程不再可以被忽略，不过气体分子之间的碰撞与气体-壁面之间的碰撞相比，仍占据主导地位。对于这种气体，一般有两大类建模方法：基于连续模型的方法以及基于粒子模型的方法。基于连续模型的方法认为气体仍然是连续的，只是在接近壁面的部分需要引入速度滑移模型。而基于粒子模型的方法则脱离了气体连续性的假设，从分子运动的角度出发对气体建模。包括直接蒙塔卡洛方法(DSMC)、格子玻尔兹曼方法(LBM)、分子动力学(MD)等。
[0006] 3)当Kn数大于0.1而小于10时，气体处于过渡流区。在过渡流区，气体的连续性假设不再满足，只能采用基于粒子模型的方法，常用的有DSMC、MD。
[0007] 4)当Kn数大于10时，气体处于自由分子流区。可以采用DSMC、MD等方法。
[0008] 在生产实际中，微流体往往处于滑移流区，因此滑移流区的建模就尤为重要。如前文所述，在滑移流区可以采用多种建模方法：DSMC方法、LBM方法、MD方法以及考虑边界滑移的N-S方程等。下面逐一分析各个方法的优劣。1)DSMC一种统计方法，它等效于求解Boltzmann方程。Boltzmann方程有严密的理论背景，它适用于从连续流到自由分子流的所有流区，并且从Boltzmann方程可以推导出N-S方程。但是由于Boltzmann方程中复杂的碰撞项，对它的求解极为困难。DSMC方法取得了巨大的成功，但是DSMC主要应用于高速流动的求解，在流动速度比较低的情况下，DSMC方法存在严重的统计噪声，计算效率较差。而在微尺度流动中，由于粘性效应的影响，在很多情况下气体流动速度较低。2)LBM方法是求解线性Boltzmann方程的一种离散方法，它在边界处理、并行化运算方面有很大优势，在最近十几年得到了很多关注。自2002年以来，人们尝试采用LBM进行微流体的计算，并取得了一定的成功。3)MD是以上方法中能得到最多的微观信息，但是所需要的计算资源也最多。4)考虑边界滑移的N-S方程是对N-S方程的扩展，使N-S方程也能用于滑移流区的建模。其中，速度滑移模型是该方法的关键。另外，该方法所需的计算资源小于前三种方法。
[0009] 在实际应用中，微流体的流动区域一般是三维的，并且流动区域的拓扑通常是不规则的(如表面粗糙度)。因此，选用考虑速度滑移的N-S方程作为建模方法，原因主要有两点：所需计算资源较少，并且得益于发展成熟的针对N-S方程的计算方法，方便处理较为复杂的边界。
[0010] 速度滑移模型是采用N-S方程对滑移流区进行建模的关键，科研人员提出了多种滑移模型描述气体在边界处的滑移，但是这些滑移模型一般是针对光滑平面提出的，对于适用于有曲率的表面甚至粗糙表面的滑移模型，研究较少。
[发明内容]
[0011] 本发明的目的就是要解决上述的不足而提供一种复杂几何边界下微流体器件中气体滑移流动的计算方法。结合求解N-S方程的有限体积法，使得该方法能求解任意几何边界下气体滑移流的流动。
[0012] 为实现上述目的设计一种复杂几何边界下微流体器件中气体滑移流动的计算方法，本方法基于计算流体力学中的有限体积方法，一般情况下，采用有限体积法进行求解的基本步骤如下：1)对计算区域进行离散，其中计算区域可以是任意形状的；2)根据求解问题的不同，采用一定的离散格式(一阶迎风格式、二阶迎风格式等)对变量进行离散，得到离散方程；3)根据特定的边界条件(压力边界、速度边界、周期边界等)，求解离散方程，直至结果收敛。对于传统的宏观问题，在壁面处气体速度等于固壁速度，采用无滑移速度边界条件；
而对于处于滑移流区的气流，气体相对壁面存在滑移速度，需要设定合适的速度边界条件，而合适的速度边界条件依赖于准确的速度滑移模型。
[0013] 著名的Maxwell滑移模型表示为：
[0014]
[0015] 式中，Us是壁面处的滑移速度，σt是切向动量适应系数(TMAC)，λ代表气体的平均分子自由程，dU/dz|w是壁面处法向速度梯度。由于在模型的推导过程中，假设气体来流平行于壁面，因此该模型适用于光滑平面。对于任意几何表面，需要适应性更广的滑移模型。
[0016] 在有限体积法的离散过程中，复杂的几何表面被离散成了一定数量的小光滑表面，设原几何表面为S，每个小光滑表面记为Sk,则S＝USk。对于每个Sk，设它的单位法向量为入射的气体分子与壁面发生碰撞，根据Maxwell模型，一部分分子(σt)被固体
表面吸收后重新喷出，喷出的角度是完全随机的，该部分称为散射分子；而另一部分分子(1-σt)则以镜面反射的形式从壁面返回气体流动区域，该部分称为折射分子。因此，在靠近壁面的气体薄层内(Knudsen层，厚度与平均自由程同级)，存在三种分子：入射分子、散射分子、折射分子。这三种分子向固体表面传递动量，而气体的滑移速度可以通过传递动量的守恒特性推导得出。
[0017] 设整个计算区域的笛卡尔坐标系为(x,y,z)，为了方便推导滑移速度，在小光滑表面Sk处建立局部坐标系(xk，yk，zk)，其中，单位法向量 作为yk轴，另外两轴根据右手定则任意选取。设气体在总体坐标系中的速度记为(u,v,w)，在局部坐标系中则为(uk,vk,wk)。在局部坐标系中，Knudsen层内分子向壁面传递动量的总通量为
[0018]
[0019] 其中，(Uk，Vk，Wk)是分子热运动的速度， 代表单个分子传递的动量，f(Uk，Vk，Wk)是分子的分布函数。函数 可以表示为 或者
式中m代表单个分子的质量，(usk，vsk，wsk)是局部坐标系中的滑移速度。总通
量P由三部分组成：
[0020] P＝Pi+(1-σt)Pr+σtPe (0.3)
[0021] 式中，Pi、Pr和Pe分别是由入射分子、折射分子以及散射分子携带的动量。气体分子在散射过程中，散射角度是完全随机的，因此散射分子携带的动量能互相抵消，Pe＝0；而在折射过程中，切向动量不变，法向方向相反，因此有Pr＝-Pi。公式(0.3)可以简化为：
[0022] P＝σtPi (0.4)
[0023] Pi可以表示为：
[0024]
[0025] 公式(0.2)和(0.5)中的分布函数f满足Boltzmann方程，这里采用由Chapman和Cowling给出的f关于Kn数的一阶近表达式：
[0026]
[0027] 其中，n是分子的数密度，k是Boltzmann常数,kc是热传导系数, R是气体常数，(C1，C2，C3)＝(Uk，Vk，Wk)，(x1，x2，x3)＝(xk，yk，zk)，(c1，c2，c3)＝(uk，vk，wk)，CioCj＝CiCj-C2δij/3。如果气体沿着x轴流动并且平行于光滑平面，只需要考虑 项，可以得到Maxwell滑移模型。然而对于复杂的几何边界，可能会出现垂直于壁面的速度分量，因此，保留公式(0.6)中所有的速度梯度项。把公式(0.6)带入(0.2)和(0.5)，并对Uk、Vk和Wk进行积分，可得：
[0028]
[0029]
[0030] 把公式(0.7)、(0.8)带入(0.4)可得：
[0031]
[0032] 根据上式，进行简单的代数运算，即可得到滑移速度：
[0033]
[0034] 公式中，与Maxwell滑移模型相比，分子上的 也可以根据壁面处的切向剪切应力得到，分母上多出的项代表yk方向的剪切应力与静压强的比值，可以被忽略。公式
(0.10)给出的滑移速度usk是局部坐标系下的结果，总体坐标系下的滑移速度(us，vs，ws)可以通过坐标变换的原理得到：
[0035]
[0036] 公式中， vk＝lku+mkv+nkw。对于等温流体，公式中与温度相关的热蠕动项可以忽略。
[0037] 在有限体积法中，以公式(0.11)代替无滑移速度边界条件，即得到了适用于复杂几何边界的气体滑移流动的计算方法。
[0038] 本发明同现有技术相比，本方法主要有两个优点：
[0039] 1)几何边界适应性强。本方法基于有限体积方法，采用一种扩展的滑移边界代替通常采用的无滑移边界，使得求解结果能够准确描述气体滑移流动的特性。本方法对问题的几何边界没有限制，理论上只要通过有限体积法能够离散的边界，都可以通过该方法进行求解。例如任意截面形状的管道流、考虑随机粗糙度的圆管流等。
[0040] 2)使用方便。N-S方程的有限体积解法经过多年的发展，已经十分成熟，并且有为数众多的基于有限体积法的商业软件，如ANSYS FLUENT14.0等，而本方法可以与这些商业软件联合使用。以ANSYS FLUENT14.0为例，利用FLUENT成熟高效的代码求解N-S方程，并通过用户自定义函数，按照公式(0.11)给出的速度滑移模型计算滑移速度，经过迭代，即可得到最终结果。用户只需要编制简单的计算滑移速度的程序，并根据不同商业软件的特点，将计算的滑移速度代替原本的无滑移边界条件即可。[附图说明]
[0041] 图1三维随机粗糙表面滑移速度的分布图。[具体实施方式]
[0042] 实例1：通过提出的计算方法，计算圆形微通道内压力驱动下的滑移流动。圆形微通道的几何边界不是平面而是圆弧面，而且由于几何形状规则，存在解析解。因此，通过该实例可以验证本方法的正确性。微通道内的流动特性可以用Poiseuille数(Po)表示：
[0043]
[0044] 公式中，F是摩擦因子，Re是雷诺数， 是流动方向的平均压力梯度，Dh是水力直径，是横截面处的平均速度。对于宏观的圆管流，不考虑稀薄效应，Po＝64。而对于圆形微通道，Po与Kn数满足下面的关系：
[0045]
[0046] 取(2-σt)/σt＝1.6，在不同Kn数下，得到的解析解和数值解如下表所示。
[0047] 表1圆形微通道内Poiseuille数的解析结果与数值结果对比
[0048]
[0049] 从表1看出，数值结果与解析结果的吻合度很好，误差在0.5％左右。因此，可以判断本方法是正确的，并且有较高的数值精度。
[0050] 实例2：计算有随机表面粗糙度对圆形微通道内滑移流动的影响。随机粗糙圆形微通道的几何边界是复杂的，常规的滑移边界条件难以处理，而采用本方法，则可以得到较好的结果。图1给出了随机粗糙微通道表面处的滑移速度，可以发现，在相对粗糙度较大的位置(波谷)，滑移速度很小；而在相对粗糙度较小的位置(波峰)，滑移速度很大。科研人员采用不同的方法研究粗糙度轮廓对二维滑移流动的影响时，得到了类似的结果。而本方法可以准确得捕捉三维随机粗糙表面的气体滑移特性。
[0051] 本发明并不受上述实施方式的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。
