复杂几何边界下微流体器件中气体滑移流动的计算方法
【专利摘要】本发明涉及流体力学领域，是一种复杂几何边界下微流体器件中气体滑移流动的计算方法，其步骤：1)基于有限体积方法，采用网格划分软件，对计算域和计算边界进行离散；2)基于Maxwell分子模型和Knudsen层内质量、动量、能量守恒，得到适用于复杂几何表面的速度滑移模型；3)利用该模型，结合流体计算商业软件ANSYS FLUENT，编写用户自定义函数，计算稀薄气体在复杂几何表面的滑移速度。本发明几何边界适应性强，能与现有的商业软件结合，使用方便。
【专利说明】复杂几何边界下微流体器件中气体滑移流动的计算方法
[【技术领域】]
[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方程等。下面逐一分析各个方法的优劣。DDSMC-种统计方法，它等效于求解 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] Us = 2入^\、, (〇?"0
[0015] 式中，Us是壁面处的滑移速度，〇 t是切向动量适应系数（TMC)，X代表气体的平 均分子自由程，如/也匕是壁面处法向速度梯度。由于在模型的推导过程中，假设气体来流 平行于壁面，因此该模型适用于光滑平面。对于任意几何表面，需要适应性更广的滑移模 型。
[0016] 在有限体积法的离散过程中，复杂的几何表面被离散成了一定数量的小光滑表 面，设原几何表面为S，每个小光滑表面记为S k,则S = U sk。对于每个sk，设它的单位法向 量为《*=(&，％，&)。入射的气体分子与壁面发生碰撞，根据Maxwell模型，一部分分子（〇 t) 被固体表面吸收后重新喷出，喷出的角度是完全随机的，该部分称为散射分子；而另一部分 分子（l-〇t)则以镜面反射的形式从壁面返回气体流动区域，该部分称为折射分子。因此， 在靠近壁面的气体薄层内（Knudsen层，厚度与平均自由程同级），存在三种分子：入射分 子、散射分子、折射分子。这三种分子向固体表面传递动量，而气体的滑移速度可以通过传 递动量的守恒特性推导得出。
[0017] 设整个计算区域的笛卡尔坐标系为（X，y，Z)，为了方便推导滑移速度，在小光滑表 面S k处建立局部坐标系（xk，yk，zk)，其中，单位法向量&作为y k轴，另外两轴根据右手定则 任意选取。设气体在总体坐标系中的速度记为（u，v，w)，在局部坐标系中则为（u k，vk，wk)。 在局部坐标系中，Knudsen层内分子向壁面传递动量的总通量为
[0018]
【权利要求】
1. 一种复杂几何边界下微流体器件中气体滑移流动的计算方法，包括以下步骤： 1) 对计算区域进行离散，其中计算区域可以是任意形状的，在有限体积法的离散过程 中，复杂的几何表面被离散成了一定数量的小光滑表面，设原几何表面为S，每个小光滑表 面记为S k,则S = U Sk,对于每个Sk,设它的单位法向量为m =(4,/?,?)，入射的气体分子 与壁面发生碰撞，根据Maxwell模型：
式中，Us是壁面处的滑移速度，〇t是切向动量适应系数（TMC)，λ代表气体的平均 分子自由程，dU/dz |w是壁面处法向速度梯度；一部分分子（〇t)被固体表面吸收后重新 喷出，喷出的角度是完全随机的，该部分称为散射分子；而另一部分分子（l-〇 t)则以镜面 反射的形式从壁面返回气体流动区域，该部分称为折射分子，因此，在靠近壁面的气体薄层 内（Knudsen层，厚度与平均自由程同级），存在三种分子：入射分子、散射分子、折射分子， 这三种分子向固体表面传递动量，而气体的滑移速度可以通过传递动量的守恒特性推导得 出； 2) 设整个计算区域的笛卡尔坐标系为（x，y，z)，在小光滑表面Sk处建立局部坐标系 (xk，y k，zk)，其中，单位法向量k作为yk轴，另外两轴根据右手定则任意选取，设气体在总体 坐标系中的速度记为（u，v，w)，在局部坐标系中则为（u k, vk, wk)，在局部坐标系中，Knudsen 层内分子向壁面传递动量的总通量为：
其中，（Uk，Vk，Wk)是分子热运动的速度，#?ΛΑ，％)代表单个分子传递的动量， f (Uk, Vk, Wk)是分子的分布函数，函数供可以表示为供=+%)、P = W(?+匕）或者 p = 式中m代表单个分子的质量，（Usk，vsk，w sk)是局部坐标系中的滑移速度，总通 量P由三部分组成： P = Pi+(l-〇t) Pr+OtPe (0. 3) 式中，Pi、P1?和Pe分别是由入射分子、折射分子以及散射分子携带的动量，气体分子在 散射过程中，散射角度是完全随机的，因此散射分子携带的动量能互相抵消，Pe = 0 ;而在折 射过程中，切向动量不变，法向方向相反，因此有已=-Pi，公式（1.3)可以简化为： P= OtPi (〇· 4) Pi可以表示为： -00 -00 -00
公式（1. 2)和（1. 5)中的分布函数f满足Boltzmann方程，这里采用由Chapman和 Cowling给出的f关于Kn数的一阶近表达式：
其中，η是分子的数密度，k是Boltzmann常数，k。是热传导系数，β = /^?2ΜΤ，R是气 体常数，（C1, C2, C3) = (Uk，Vk，Wk)，（X1, x2, x3) = (xk，yk，zk)，（C1, c2, c3) = (uk，vk，wk)，Cit5C j = CiCrC2 δ ,j/3 ； 3)把公式（I. 6)带入（I. 2)和（I. 5)，并对Uk、Vk和Wk进行积分，可得：
把公式（1. 7)、（1. 8)带入（1. 4)可得：
根据上式，进行简单的代数运算，即可得到滑移速度：
公式中，与Maxwell滑移模型相比，分子上的Svt 也可以根据壁面处的切向剪切应力 得到，分母上多出的项代表yk方向的剪切应力与静压强的比值，可以被忽略，公式（1. 10) 给出的滑移速度Usk是局部坐标系下的结果，总体坐标系下的滑移速度（us，vs，w s)可以通过 坐标变换的原理得到：
公式中，5/也=以/0Χ+?^咐ν + ?0/δζ，Vk = lkU+mkV+nkW，对于等温流体，公式中与温度 相关的热蠕动项可以忽略， 在有限体积法中，以公式（1. 11)代替无滑移速度边界条件，即得到了适用于复杂几何 边界的气体滑移流动的计算方法。
【文档编号】G06F17/50GK104376183SQ201410719871
【公开日】2015年2月25日 申请日期:2014年12月1日 优先权日:2014年12月1日
【发明者】闫寒, 张文明 申请人:闫寒, 张文明