基于气体动理论的二维Karman涡街数值模拟

吴俊林, 李志辉, 蒋新宇, 彭傲平

(1. 中国空气动力研究与发展中心超高速空气动力研究所, 四川绵阳 621000; 2. 国家计算流体力学实验室, 北京 100191)

引 言

Boltzmann方程[1]是描述气体分子运动状态和运动规律的确定论基本方程, 主要描述单原子气体或简单气体的分子输运特性. 不考虑外力场影响的Boltzmann方程如公式(1)所示

(1)

其中,f为气体分子的速度分布函数；t，r，ξ分别为时间、 位置空间坐标、 分子速度；ξr为碰撞分子的相对运动速度；σ为微分碰撞截面, 与分子之间的作用力模型相关；Ω为两个分子碰撞后相对速度的方向；R3表示ξ1的三维速度空间积分.

Boltzmann方程通过对气体分子速度分布函数矩积分得到流动的宏观参数, 可描述从稀薄流到连续流跨流域气体流动由非平衡态向平衡态演化的过程[2]. 然而, 直接采用理论分析或数值方法研究Boltzmann方程很难实现, 困难在于[3-4]: (1) 几率密度分布函数f具有至少7个自变量, 维度太高; (2) 碰撞项高度非线性, 且具有高维积分特点, 并与分子碰撞模型相关, 极其复杂. 虽然Aqarwal等[5]、 Kolobov等[6]在求解Boltzmann方程或广义Boltzmann方程(generalized Boltzmann equation， GBE)方面取得了一些进展, 但也并未实现对碰撞项的直接求解, 并且由于求解过程复杂, 计算代价很大, 导致很难应用到工程中.

正是由于Boltzmann方程具有高维度、复杂碰撞积分项的特点, 对于跨流域非定常问题, 特别是稀薄过渡区的低速非定常流动, 其数值求解应该考虑基于Boltzmann的运动模型方程[7-12]. 近年, Polikarpov等[8]通过求解非定常Shakhov模型方程得到了气流穿过方形切口的瞬态流动, 并说明定常流场建立的时间与背压压比和气体稀薄度相关. Lihna-ropoulos等[9]通过数值求解时间依赖的BGK运动方程分析了圆柱管道中的气体启动过程. Chigullapalli等[7]求解Boltzmann-ESBGK 模型方程发展了一种三维非定常稀薄流动求解器. 实践表明, 基于Boltzmann模型方程的数值求解方法是模拟稀薄流到连续流跨流域非定常流动的有效手段.

为了求解Boltzmann模型方程以得到跨流域数值模拟方法, 李志辉等在速度空间采用离散速度坐标法, 位置空间应用计算流体力学有限差分方法, 建立起从稀薄流到连续流的跨流域气体动理论统一算法(gas-kinetic unified algarithm, GKUA)[13-14], 并在定常流动研究中得到了较好的应用[15-16]. 实际上, 统一算法的求解过程是通过长时间非定常输运模拟得到最终的定常流动状态, 因此可通过技术改进得到整个非定常气体流动演化过程. 前期已分析了一维、二维跨流域非定常流动特征与稀薄效应对非定常流场影响机制[17], 并用于求解分析稀薄平面喷流扩散进入真空环境的非定常过程[18]. 结果表明, 基于GKUA的非定常流动求解器对于高真空自由分子流、稀薄流、过渡流等流区的超声速非定常流动具有较好适用性. 因此, 本文对连续流区的经典Karman涡街问题进行模拟, 以确认该跨流域非定常流动模拟算法对于连续流区低速流动的适应性.

1 控制方程及计算方法

1.1 控制方程

GKUA能够求解各种模型方程[19], 包括BGK[20], ES-BGK[21]， Shakhov[22]或者其他模型. 对于氮气或者空气, 室温状态下分子的转动能就已完全激发[23], 因此采用考虑转动自由度影响的双原子气体动理论模型方程能够更好地描述氮气流动. 此外, 众多学者已经对高温条件下多原子气体分子的振动能开展了有意义的分子动力学建模分析和数值计算工作[24-27], 后续如果研究高速、高温气体动力学问题时应该加以考虑.

基于对转动自由度松弛变化特性的Rykov模型研究[28-30], 采用转动惯量描述气体分子的自旋运动, 将分子总角动量守恒作为一个碰撞不变量, 在求解Boltzmann模型方程的统一算法理论框架[13-16]下, 构建考虑转动非平衡效应的Boltzmann模型方程

(2)

式中, 下标i表示空间维度方向; 角标t， r分别表示弹性碰撞和非弹性碰撞, 如νt和νr是弹性碰撞频率和非弹性碰撞频率;μ为黏性系数;Z为非弹性碰撞松弛因子；xi为位置空间坐标；ξi为气体分子速度分量；c为气体分子热运动速度；n，T，P分别为气体分子数密度、气体流动温度和压力；k为Boltzmann常数；m为分子质量；δ对于分子间相互作用规律来说是一个常数δ=μt/(mnD),D是自扩散系数.

物理空间宏观流动参数可以由f0和f1在速度空间上积分[13,28]求得

kTr=ε,Pt=nkTt,P=nkT

黏性系数μ和非弹性碰撞松弛因子Z可定义为如下的形式[29-30]

μt=μ(Tt)=μ(T*)t2/3/φ(t),t=Tt/T*

φ(t)=0.767+0.233t-1/6exp[-1.17(t-1)]

对于氮气N2而言,T*=91.5 K,δ=1/1.55.

1.2 无量纲化和速度空间的离散

采用气体动理论数值计算方法可直接捕捉速度分布函数随时间的演化. 空间坐标、时间、数密度、流动速度、温度、黏性的无量纲参考量为[13-16]

而能量(包括平动能和转动能)的无量纲参考量以及分布函数f0和f1的无量纲参考量分别为

mn∞(2RT∞)3/2,n∞(2RT∞)-3/2

mn∞RT∞(2RT∞)-3/2

对于二维流动, 为了对速度空间进行降维处理, 在对速度分布函数方程(2)应用离散速度坐标法之前, 需要引入3个约化速度分布函数, 对模型方程和基于速度空间的宏观流动矩积分进行约化处理, 这样可以减少对计算机内存的需求, 提高计算效率.

得到关于二维约化速度分布函数gi(t,x,y,ξx,ξy),i=1,2,3的方程. 在此基础上, 使用离散速度坐标法对速度分量ξx，ξy进行数值离散.

应用离散速度坐标方法将物理空间的坐标(x,y)转换到计算平面(ζ,η), 模型方程转化成为每个离散速度坐标点上包含非线性源项的双曲型守恒方程

(3)

在离散速度坐标法中, 采用合适的积分规则对分布函数在速度空间进行矩积分, 其中积分所依赖的分布函数仅在某些离散速度坐标点上有确定值. 对分布函数在全速度域上的积分可基于这些离散速度点上的分布函数值得到. 求解方式是采用Gauss-Hermite数值积分法, 将全域积分转变为基于权重因子的求和. 基于此, 构造某些自适应选取的离散速度坐标点(ξxσ,ξyε), 以及对速度空间这些点上的分布函数值进行合理迭代求解至关重要, 这将直接影响数值模拟算法的收敛性和适用性. 值得注意的是, 这些离散速度坐标点确定的速度分布函数值必须保持正定性.

1.3 适于非定常流动计算的差分方法

对于非定常流场的数值计算, 除了差分格式须满足耗散控制、色散控制和激波控制条件以保证计算过程中的稳定性、不产生虚假波动以及良好地捕捉激波外, 还应特别满足保频谱原则, 并且要求计算网格和边界计算方法要与内点计算匹配协调[31-32]. 一般采用固定时间步长的方法进行显式的时间推进, 但时间步长必须满足流场计算的格式稳定性条件.

基于非定常时间分裂方法[4,13-16,19]将控制方程(3)分裂成为包含非线性源项的碰撞松弛方程和计算平面(ζ,η)上的对流运动方程. 采用3阶WENO格式[19,33]离散速度分布函数方程的对流项. 源项的碰撞松弛积分采用3阶Runge-Kutta 方法求解. 考虑到在真实气体流动中对流运动和碰撞松弛具有相似的过程, 在数值计算中需要前向时间步长和后向时间步长的耦合迭代. 基于此, 在时间推进上采用基于非定常时间分裂法的3阶显式Runge-Kutta[34]时间推进, 得到在各个离散速度坐标点处数值求解约化速度分布函数方程的气体运动论耦合迭代数值格式

(4)

其中, 下标s表示源项,ζ,η为计算平面坐标. (4)式表示在Δt时间步长之后的分布函数值通过半时间步长(Δt/2)的两次松弛演化得到, 这样可以提高该差分离散格式的时间精度.

方程(4)的每一项对应于非线性方程(3)的数值计算, 其中Ls,Lη和Lζ分别表示碰撞松弛源项以及两个计算平面坐标方向的对流运动项的差分方程.

引入3阶WENO格式[35]对计算平面进行扫描. 对于时间项, 可以使用Runge-Kutta 法(RK-3)来提高精度.

数值计算过程中的时间步长Δt由格式稳定条件约束

其中, CFL为时间步长调节系数, 一般CFL<1.

2 Karman涡街形成及运动的数值模拟

2.1 问题描述

一定条件下的定常来流绕过钝体时, 物体两侧可能会周期性地脱离出旋转方向相反、排列规则的双列线涡, 经过非线性作用后会形成Karman涡街. 以圆柱绕流为例, 低Reynolds数时圆柱后方有一对尾涡; 随着Re的增加, 尾迹出现震荡; 继续增加Re, 近尾迹中涡交替脱落, 进入下游, 形成持续较长距离的稳定涡列, 即Karman涡街[36]. 本文选取Re=1 200 时作为典型计算条件, 对应的计算状态如表1所示, 气体为氮气.

表1 低速Karman涡街算例的典型计算状态Table 1 Classical computing condition for low-speed Karman vortex street case

2.2 圆柱尾迹区的Karman涡街

圆柱直径d=1 m, 数值计算过程中流场网格关于x轴对称(如图1所示), 物理空间的网格量为(301×201), 分子速度空间的离散速度点为(16×16), 关于速度空间的矩积分采用Gauss-Hermite积分法. 圆柱尾迹区Karman涡街的形成、发展、运动及耗散过程如图2所示, 图中给出的是不同时刻流场压力等值线云图. 初始时刻流场完全处于来流状态, 壁面扰动从t=0时刻开始传播. 对于低速流动, 扰动向全流场传播. 在图2(a)所示的t=0.014 s 时刻, 物面边界扰动已经传播到离圆柱较远的距离, 圆柱周围流场基本成形. 此时能够明显看到圆柱头部区域的高压、 圆柱上下两侧的低压以及圆柱后方的扰动区都是上下对称的. 到图2(b)所示的t=0.024 s时刻, 圆柱后端尾迹区出现了上下对称的两个低压旋涡结构. 可以看到, 低压旋涡结构的产生是由逆压梯度导致的流动分离. 分离出的旋涡持续向下游流动, 并逐渐发展、扩大. 到图2(c)所示的t=0.057 s时刻, 尾迹区两个涡在向下游流动的同时相互接触、互相影响, 但此时两个旋涡依然是对称结构. 而在图2(d)的t=0.569 s时刻, 尾迹区两个旋涡结构“对称破缺”, 上下两个涡交替“扩大—收缩”, 并继续向下游流动. 这种对称性被打破的现象来源于流场的非稳态演化, 与外界干扰、网格结构等并无关系. 随着旋涡结构的持续运动、发展、相互作用, 到图2(e)的t=0.687 s时刻形成了“拟序结构”的尾迹区旋涡序列, 并最终形成了图2(f), (g)所示“Karman涡街”非定常流动状态. 其中, 图2(f)和(g)分别是半个周期的 Karman涡街运动耗散过程, 旋涡结构呈现上下交替演化现象, 从图中可以看到Karman涡街脱落、发展、运动、耗散的整个过程.

图1 低速圆柱绕流的物理空间网格Fig. 1 Grid system in physical space for low-speed flow past a cylinder

(a) t=0.014 s

图3给出了圆柱绕流尾迹区在半个周期内的Karman涡街交替脱落现象. 从流场的温度等值线云图中能够更明显地看到圆柱上下顶点稍靠后的位置交替“甩出”涡旋结构, 及其涡脱离的整个过程. 有意思的是, 甩出的涡旋结构为低压、低温, 而另一侧对应被压缩的区域则是高压、高温. 产生的Karman涡街在离圆柱尾部4d距离以后逐渐被耗散掉. 图4 详细绘出了半个周期内的流线图, 展示了圆柱尾迹区下半部甩出旋涡到上半部甩出旋涡的交替演化过程.

图3 圆柱绕流尾迹区半个周期内的Karman涡街交替脱落现象(温度)Fig. 3 Alternately dropping phenomenon of the Karman vortex street in half period at the wake-region of flow past a cylinder(temperature contours)

图4 圆柱绕流尾迹区半个周期内的流线图Fig. 4 Streamlines in half period at the wake-region of flow past a cylinder

图5给出了本文GKUA得到圆柱上下表面的压力分布与连续流区N-S方程数值计算结果的比较, 其中N-S方程计算中分别采用了层流模型和两方程k-ε湍流模型(RANS). 可以看到, GKUA得到的表面压力分布规律与连续流区经典N-S方程基本一致. 对于该低速流动, 迎风面压力值略高于背风区, 且GKUA得到的迎风面压力分布与层流模型和湍流模型都吻合较好. 背风面压力值的差别相对来说要大一些, 但也在可控范围内(最大偏差低于2.3%). 同时可以看出, 上下表面压力分布的不对称性正体现出了流动的非定常特性, 这跟尾迹区涡的交替脱落紧密相关.

(a) Upper-half part of the cylinder surface

该状态下圆柱绕流数值计算得到的Strouhal数与文献[37]中的实验值比较如表2所示. 可以看到, 本文数值计算得到的表征Karman涡街非定常流动特征的Strouhal数与实验值和公认理论值0.21都是非常接近的. 考虑到数据量的问题, 在数值计算过程中并没有保存任意时刻的流场结果, 而是隔一定时间保存一次, 导致Karman涡街的周期估算值存在一定偏差. 鉴于此, 本次圆柱绕流尾迹区Karman涡街的数值模拟结果是可信的.

表2 圆柱绕流计算得到的Strouhal数与文献比较Table 2 Comparison of the Strouhal numbers from simulation and experiment

Karman涡街的尺寸比定义为b/a, 该参数可以用来描述涡列, 其中a为同涡列中相邻旋涡之间的距离,b为两列旋涡的间隔(见图6). von Karman[38]给出涡列稳定的必要条件为尺寸比b/a≈0.281. 虽然后人指出Karman涡街稳定的尺寸比并不严格等于 0.281[39], 但不同计算条件下相差不多. 本文数值计算圆柱Karman涡街的尺寸比b/a=0.283, 与 von Karman理论推导的结果非常接近(见表3), 进一步验证了气体动理论统一算法模拟Karman涡街等非定常流动问题是适用的.

图6 Karman涡街尺寸比示意图Fig. 6 Sketch of the length ratios for Karman vortex street

表3 圆柱Karman涡街尺寸比的计算值与理论值比较Table 3 Comparison of the length ratios from simulation and theory analysis

2.3 非对称尖劈产生的涡街

圆柱绕流的Karman涡街是由流动的对称性破缺引起的, 而非对称尖劈绕流则由于流场本身的非对称性, 诱导产生的Karman涡街在空间上是非对称的, 但在时间上呈现出周期变化的规律. 非对称尖劈的物理空间流场网格如图7所示, 规模为58 403个网格点, 分子速度空间采用(32×32)的Gauss-Hermite积分法. 图8给出了一个周期内非对称尖劈诱导涡街的流场压力等值线云图, 其中t=0表示该周期的起始时刻, 并不是初始时刻. 可以看到, 上下尖劈诱导产生的涡系结构在向下游发展时是非常规整的, 形成了相互干扰又彼此独立的拟序结构. 在图8(a)该周期开始时刻, 上部尖劈出现诱导涡分离现象; 到图8(b)Tp/4时刻上部诱导涡发育完全, 脱离尖劈向下游运动; 到图8(c)Tp/2时刻下部尖劈开始产生诱导涡; 到图8(d)3Tp/4时刻下部诱导涡发育完全, 也向下游运动; 到图8(e)t=Tp时刻上部尖劈又产生新的诱导涡. 如此循环往复, 构成了非对称尖劈诱导Karman涡街的非定常流动过程.

图9给出了一个周期内t=Tp/2,t=Tp两个时刻的非对称尖劈绕流下游的涡系结构发展, 其中展示的是流线图. 可以看到, 非对称尖劈下游涡系结构比圆柱绕流下游的拟序结构要复杂得多, 前者上下尖角诱导产生的涡核大小略有差别, 且相隔半个周期的涡系结构并不关于y轴对称. 这些差异是由尖劈的非对称性引起的.

图7 非对称尖劈绕流的物理空间网格Fig. 7 Grid system in physical space for flow past an asymmetric wedge

(a) t=0

3 结论

本文在气体动理论统一算法框架下, 应用离散速度坐标法和3阶WENO格式, 以及3阶Runge-Katta时间推进法, 数值求解考虑转动自由度影响的Rykov模型方程, 得到了能模拟双原子气体跨流域非定常流动问题的一种确定论数值计算方法. 对经典的二维Karman涡街非定常流动现象的模拟, 说明了该跨流域非定常模拟方法能够适用于连续流区的低速流动, 并复现了关于Karman涡街流动现象的一些规律和理论:

(1)低速圆柱绕流的Karman涡街是由于流动演化的对称性破缺而产生周期性变化规律, 与圆柱外形、流场网格等无关;

(2)圆柱绕流的Karman涡街不仅在时间上有周期性变化规律, 而且圆柱上下诱导产生的拟序结构在空间上也具有对称性, 这种空间对称性并不体现在同一时刻, 而是体现在半个周期上的空间对称;

(3)圆柱绕流的诱导涡产生必然是发生在圆柱上下顶点靠后的逆压梯度区;

(4)非对称尖劈诱导产生的Karman涡街并不具有空间对称性, 但依然表现出时间上的周期性变化, 且其涡系结构更加复杂.

致谢本工作得到973计划(2014CB744100), 国家自然科学基金(11902339, 11325212)资助.