《武汉工程大学学报》 2018年06期
673-677
出版日期：2018-12-28
ISSN:1674-2869
CN:42-1779/TQ
基于气体动理学格式的平板绕流数值模拟
描述流体运动有两种方法，第一种是从宏观角度描述流体，例如Euler方程，Navier-Stokes方程等；另一种是从微观角度，它采用粒子速度分布函数描述高维相空间的粒子数密度，如Boltzmann方程。基于介观气体动理论的气体动理学格式（Gas Kinetic Scheme，GKS）从微观角度为气体流动数值计算提供了很好的方法。Xu Kun等通过构造不同的初始气体分布函数以及平衡态分布函数，建立了等价于求解NS方程的GKS-NS格式以及理论上能够求解全流域流动的统一算法（UGKS格式）［1-3］。Jin等［4］在格式中引入速度空间自适应技术，加快了计算速度，减少了多尺度流问题中的内存需求。为高效模拟可压缩流动，学者研究了高精度气体动理学格式［5-7］。另外，在非结构化网格［8］、动网格［9］、浸入边界［10］中应用气体动理学格式进行数值模拟，得到了精确的模拟结果。气体动理学格式在计算过程中表现出稳定性，使得气体动理学格式得到越来越多的应用。绕流是流体力学研究的重要内容之一［11-13］，在工业领域涉及面广、影响面大，有重要的研究价值。目前国内外应用气体动理学格式在绕流方面的研究较少［14］。本文将介绍气体动理学格式在平板绕流中的应用，选取合适的边界条件构造模拟二维平板绕流的物理模型。主要探究平板在不同间距条件下不同的绕流流场特性，揭示气体动理学格式在平板绕流数值模拟中良好的应用前景。1 气体动理学格式的基本原理1.1 控制方程在无外力场作用下，二维气体动理学BGK-Shakhov方程为：[ f t+u f x+v f y=f+-fτ] （1）式（1）中：[f]是气体分布函数；[u]是粒子切向速度；[v]是粒子法向速度。[τ=μp] （2）式（2）中：[τ]是碰撞时间；[μ]是动力粘性系数；p是压力。修正后的分布函数为：[f+=g1+1-Prc·qc2RT-55pRT] （3）式（3）中：Pr是普朗特数；[c]是随机速度矢量；[q]是热流矢量；R是通用气体常数；T是温度。平衡态分布函数为：[g=ρ(λπ)K+22e-λ(u-U2+v-V2+ξ2)] （4）式（4）中：[ρ]是密度；U，V是宏观速度。[ξ2=ξ21+ξ22+?+ξ2K] （5）式（5）中：K是内部自由度。[λ=m2kT=12RT] （6）式（6）中：m是分子质量；k是玻尔兹曼常数。内部自由度K和比热比[γ]之间的关系如下所示：[γ=K+4K+2] （7）宏观守恒量与气体分布函数的关系为：[W=ρρUρVρE=ψfdΞ] （8）式（8）中：[ψ=1uv12(u2+v2+ξ2)]，[dΞ=dudvdξ]由于气体分子在碰撞过程中满足质量、动量和能量守恒，因此有守恒性约束条件要求：[(f+-f)ψdΞ=0] （9）采用有限体积法对公式（1）进行离散。物理空间，时间空间和粒子速度空间被划分为许多小的控制体。在二维空间中，[Ωi,j=?x?y]，[?x=xi+1/2,j-xi-1/2,j]，[?y=yi+1/2,j-yi-1/2,j]。粒子速度空间被划分为矩形网格进行离散，间距为[?u和?v]。（k，l）速度间隔的中心位于（[uk，vl]）=[(k?u，l?v)] 。则[tn]时刻的平均气体分布函数为：[fxi,yj, tn,uk,vl=fni,j,k,l=][1Ωi,j?u?vΩi,j?u?v-∞+∞fx,y,tn,u,vdxdydudv] （10）宏观守恒量的演化更新为：[Wn+1i,j=Wni,j+1?xFi-1/2,j-Fi+1/2,j+][ 1?y(Gi,j-1/2-Gi,j+1/2)] （11）式（11）中：[Fi+12,j=0?tuΨfi+12,j,k,ldΞdt] （12）[Gi,j+1/2=0?tvΨfi,j+1/2,k,ldΞdt] （13）气体分布函数的演化更新为：[fn+1i,j,k,l=(1+?t2τn+1)-1[fni,j,k,l+1Ωi,jtntn+1m?Smumfm,k,ldt+] [ ?t2f+n+1i,j,k,lτn+1i,j+f+ni,j,k,l-fni,j,k,lτni,j]] （14）气体动理学的关键正是利用上述通解得到单元界面上随时间演化的气体分布函数，从而计算出数值通量，得到下一时刻的物理量［15］。1.2 计算域和边界处理计算域为长方形，长L=400 mm，宽W=100 mm。平板布置方式为单平板（布置1）、前后平板（布置2）和上下平板（布置3）3种情况，具体布置方式及计算域如图1所示。布置1、2中平板l=30 mm，布置3中平板长度l=20 mm。求解流动问题需要给出适当的边界条件，边界条件的处理方式对计算结果的精度和计算的稳定性有很大的影响，各种边界条件处理如下：1）滑移壁面边界： 在边界法线上根据镜面反射原则取值，即把和边界节点相邻节点的各流动量矢量值，取相反值赋到计算区域外的虚拟网格节点上。而在边界切线方向上直接取为计算域内与虚拟网格节点相邻节点流动量的值。2）无反射边界：根据零梯度条件，可通过外插的方法由流场的内点得到宏观边界条件。零梯度条件为：[ ? n=0] （15）3） 速度入口边界： 在射流入口处施加速度边界[u=U0]，[v=0]。平板绕流中，入口边界为速度入口，出口边界为无反射边界，上下边界为滑移壁面边界条件。网格数为401[×]101的结构化网格。1.3 算法准确性对比验证选用典型的Poiseuille流即流场无平板时验证程序的准确性。Poiseuille流的解析解为：[uy=-GμL22(yL-y2L2)] （16）式（16）中：[μ]为流体的动力粘性系数；L为管道宽度；[G=-dp/dx]为管道内的压力梯度。当y=L/2时速度取得最大值：[umax=-1μdpdxL28] （17）以最大速度做无量纲处理，则：[uumax=4(yL-y2L2)] （18）计算域尺寸为[0x2，0y1]。采用结构化网格，网格数为128[×]64。上下边界采用滑移壁面边界条件；左边界采用压力进口边界条件，入口压力为1.01 MPa；右边界采用压力出口边界条件，出口压力为1.00 MPa。在流域内取3个截面，section 1，section 2和section 3分别表示与流向垂直的1/4位置截面、中间位置截面和3/4位置截面，截面示意图如图2所示。3个截面流向速度U与解析解的相对误差[δ]如图3所示，从图3中可以看出：管道内流体的速度与解析解最大误差小于1.8%，且误差较大的地方出现在上下边界处，管道中心误差很小，故气体动理学数值计算结果与解析解基本重合，初步验证了程序的正确性。2 平板绕流数值模拟利用气体动理学格式的基本理论，计算二维平板绕流的流体流动情况。假设流体从左端均匀流入。2.1 单个平板绕流对不同雷诺数的单板绕流流场进行模拟。图4给出了布置1在不同雷诺数下的涡线图。由图4可以看出：Re[]100时，流场处于对称尾流区，没有涡的脱落，绕平板分离后的流体在平板后形成对称的漩涡，并且只对平板附近流场有影响。但随着雷诺数的增加，涡不断被拉长，逐渐失去对称性，涡开始在平板两侧周期性的脱落，形成一系列的涡即著名的卡门涡街。2.2 前后平板绕流对雷诺数Re =300，两板间距为10 mm，30 mm和50 mm三种情况下前后平板绕流流场进行了模拟。图5和图6分别给出了布置2在不同板间距下的涡线图和涡量图。由图5和图6可知，对于低雷诺数下前后排列的平板绕流，随着板间距的逐渐增大，其流场结构呈现以下特征：1）当平板间距较小时，两个平板离的很近，此时流场形态类似于一个平板，流场处于对称尾流区，没有涡的脱落，如图5（a）和图6（a）所示；2）当间距稍微增大时，上游平板脱落的剪切层附着于下游平板，在两平板间出现涡量场，但仍没有涡的脱落，如图5（b）和图6（b）所示；3）当间距很大时，上游平板的涡脱落后击打在下游平板上，并在下游出现卡门涡街，如图5（c）和图6（c）所示。当平板距离比较小时没有卡门涡街的形成，当间距较大时则有卡门涡街的形成。2.3 上下平板绕流对雷诺数Re=300，两板间距为10 mm，20 mm和40 mm三种情况，上下板绕流流场进行了模拟。图7和图8分别给出了布置3在不同板间距下的涡线图和涡量图。由图7和图8可知，对于低雷诺数下，上下排列的平板绕流，随着板间距的逐渐增大，其流场结构呈现以下特征：1）当两平板间距较小时，每个平板边缘处有涡的脱落，并在距平板较近的地方涡合并形成共同的涡，出现卡门涡街，如图7（a）和图8（a）所示；2）随着板间距的增大，在上下平板后面附近区域各自形成自己的涡，在离平板比较远的地方两列涡开始逐渐合并形成一列涡，如图7（b）和图8（b）所示；3）当板间距很大时，在两平板下游各自出现卡门涡街，两板之间的影响较小，如图7（c）和图8（c）所示。3 结 语本文基于气体动理学格式对平板绕流进行了数值模拟，主要结论如下：1）对于单平板绕流，当雷诺数较小时，流场处于对称尾流区，没有涡的脱落；当雷诺数较大时，在平板下游会形成卡门涡街。2）对于前后平板绕流，左右平板距离较小时没有形成卡门涡街，但当平板距离较大时就会形成卡门涡街。3）对于上下平板绕流，当平板间距比较小时，相当于一个平板形成的绕流；而当两个平板距离比较大时，在平板下游较近距离处形成各自的涡旋，并且在离平板比较远的地方两列涡开始逐渐合并形成一列涡。