镜像对称顶盖驱动方腔内流过渡流临界特性研究1)

2022-10-05安博孟欣雨桑为民

力学学报 2022年9期

关键词：顶盖湍流周期性

安博孟欣雨桑为民,2)

* (西北工业大学航空学院,西安 710072)

† (中国空气动力研究与发展中心结冰与防除冰重点实验室,四川绵阳 621000)

** (翼型、叶栅空气动力学国家重点实验室,西安 710072)

引言

过渡流临界特性作为流体力学中的经典问题长期受到学界的关注.Hof 等[1]在Science上撰文强调,过渡流临界特性的研究对揭示流动物理本质起到了至关重要的作用,同时对流场演化模式起到了决定性的作用.他们认为,流场演化是个极为复杂的过程,需要细节化的分类探索来准确把握诸如湍流等复杂流场的物理特性.Avila 等[2]也在Science上撰文阐述,过渡流研究对认识流场的重要作用,对比经典的Landau-Ruelle-Takens 流场演化模式,他们认为混沌的空间扩散是湍流特性的决定性过程和固有本质.同时Graham[3]在Nature上也介绍了过渡流的概念,文中再次强调了过渡流临界特性研究的重要性,同时指出过渡流的相关研究为解释更为复杂的流动现象铺平了道路.此外诸多学者[4-13]在Annual Review of Fluid Mechanics和Journal of Fluid Mechanics上先后强调了过渡流临界特性研究在流体力学研究中的应用价值.虽然他们就不同视角从不同层面探讨了不同的内容,但是学者们一致肯定了过渡流研究的重要性和必要性,在正确认识物理本质的同时不断推进本学科的蓬勃发展,为解决更为复杂的流动问题打下了坚实的基础.

作为过渡流临界特性的核心研究内容,流动分岔点(flow bifurcations)表征的是流场不同流动状态(定常流动、非定常周期性流动、非定常准周期性流动、湍流) 的物理特性临界点.常见的流动分岔点,比如Hopf 流动分岔点,它的出现意味着此时的流动随雷诺数增加已经从定常状态演化至非定常周期性状态,流动会随时间产生周期的循环变化.根本原因是流场稳定性被破坏了,取而代之的是流场周期性特征.随着雷诺数进一步增加,流动的非定常周期特性也逐渐遭到破坏,此时的流动虽然继承了部分周期性流动的物理特性,但逐渐出现了准周期性流动的典型特征,即基于庞加莱映射的双环形结构.此时的流动不再以单一频率振荡,流动表现出长周期和短周期的双频振荡.而Neimark-Sacker 流动分岔点正是非定常周期性流动和非定常准周期性流动的临界点;Period-doubling 流动分岔点出现后,标志着流场从非定常周期性流动跃变至湍流;再比如定常状态下,区别不同拓扑结构流场解的Saddle-node流动分岔点和在研究流场拓扑结构时捕捉到的Pitchfork 流动分岔点.这样的基础研究,对准确认识流场的物理特性意义巨大.

作者在之前的研究工作[14-15]中针对三角腔顶盖驱动内流,经典方腔顶盖驱动内流以及上下边双驱动方腔内流开展了全面的流场稳定性分析研究.较为完整地揭示了以上流动问题的流场过渡流临界特性.同时发现流场拓扑结构的 π 旋转对称性与流场稳定性密切相关,对流场的演化也有较大影响.根据之前的研究经验发现对称驱动条件对于腔体内流的过渡流临界特性影响非常大.当驱动条件布置在相向的两个边界时,会导致流场内的镜像和 π 旋转对称性,而这些流场对称性又与流场的稳定性息息相关,对更好地理解流场物理规律有重要意义.而单边镜像对称驱动内流的相对缺乏,没有形成完整的体系化的研究.本文独特的单边镜像驱动条件使得顶盖的剪应力分布相较于常规边界驱动内流发生了巨大变化,导致了流场具有更强的失稳性,其流动机理也更为复杂.为了进一步明晰镜像对称性与流场稳定性的关系及对流场演化的影响规律,针对镜像对称顶盖驱动内流开展过渡流临界特性研究.

1 不可压缩LBM 计算模型

本文的研究目的在于揭示镜像对称顶盖驱动方腔内流的过渡流临界特性,数值模拟工作以低马赫数和低雷诺数为主.所涉及马赫数为Ma=0.173 2,雷诺数为Re≤20 000,因此选取传统格子Boltzmann方法[16-19]中的经典单松弛碰撞迁移模型作为本文数值模拟的计算模型,且使用了目前应用最广泛的LBGK D2Q9 模型[20].其中平衡态分布函数的构造为

式中,ωi和ei分别对应离散时空模型中不同离散方向的权系数和离散速度,即

其中,c=Δx/Δt=1 为格子速度,Δx为网格步长,Δt为时间步长.格子Boltzmann 控制方程为

其中,fi为碰撞迁移前的分布函数,Fi为碰撞迁移前的分布函数,τ 为松弛时间. ρ 和u为流体粒子的宏观密度和速度,即

2 数值模拟背景

2.1 计算域构建及参数设计

均匀直角网格具有网格质量好,鲁棒性高等特点,其网格结构天然契合格子Boltzmann 方法的碰撞迁移理论,因此在LBM 数值模拟中得到了广泛应用.本文使用的均匀直角网格,网格分辨率为1024×1024,网格步长 Δx=1/1024,在之前的研究工作中[15]我们已经开展了相关的网格独立性验证,研究结果表明这个网格尺度足以保证计算结果的精确性和可靠性.尤其对于腔体内流流场在Re≤20 000 的计算状态下,能够保证y+＜0.5,具备捕捉边界层流动信息的能力.如图1 所示,计算域特征长度L=1.0 是方腔的边长.顶盖镜像对称驱动条件为

本文在过渡流临界特性分析研究中设计了两个数值模拟信息采集点(如图1 所示)Plm(x=0.25,y=0.5)和Prm(x=0.75,y=0.5) 用于记录流场中局部速度随时间的变化曲线.同时,为了研究流场拓扑结构的镜像对称性,设计了对称性参数 ξ,定义为

图1 计算域Fig.1 Computational domain

其中,ulm和vlm是点Plm(x=L/4,y=L/2)水平和垂直速度分量,而urm和vrm是对应的点Prm(x=3L/4,y=L/2)的水平和垂直速度分量.

2.2 边界条件处理

本文数值模拟的边界条件均为平直边界(方腔四个边).其中顶盖为驱动边界,其余三边均为物面边界.为此本文采用了经典的非平衡态外推格式[21].该格式的基本思想是,将边界节点上的分布函数分为平衡态和非平衡态两部分.其中,平衡态部分由平衡态分布函数的定义近似获得,而非平衡态部分则用非平衡态外推法求解.

如图2 所示点D,E和F为远场边界点,根据LBM 的演化(碰撞迁移)原理可知,在每次演化之前需要求解每个点的分布函数,对于E点其分布函数可看作两部分: 平衡态分布函数和非平衡态分布函数,即

图2 平直物面边界Fig.2 Straight wall boundary

因此,可以近似求解E点的分布函数

若考虑边界点的碰撞过程,则边界点E的分布函数为

综上,可以确定壁面边界和驱动边界点的分布函数,各边界点宏观物理量构造如下

3 计算结果与分析讨论

3.1 定常流动

根据流动的演化机理,随着雷诺数的增加流动会从定常状态演化为非定常状态[22-23].本文选取Re=1000时的数值模拟结果作为镜像对称顶盖驱动内流定常结果的代表.如图3(a)和图3(b)所示,分别介绍了该雷诺数下的流线图和涡量图.由流场拓扑结构可见,由于对称驱动的作用,流动在此时保持了非常好的对称性,两个对称主涡几乎占据了整个计算域,主涡的周边存在着两对对称次级涡.从主涡和次级涡的演化和分布规律可以观察到腔体内流的典型特征.此外,本文数值模拟结果中其余的所有涡量图都使用了如图3(b)所示的统一色条.

图3(c)描绘了信息采集点Plm和Prm的水平速度分量及对称性参数随时间的变化曲线,其中t是无量纲流场演化步数.从三个曲线中可以看出此时的流动收敛于一个定常状态且进一步从数值方面证明,此时的流动是严格对称的(ξ ≡0).相较于之前的研究工作[14-15],没有发现类似三角腔顶盖驱动内流的Saddle-node 流动分岔点,说明对于定常流动,只存在一种流场解,这与经典的顶盖驱动方腔内流研究结论一致.

图3 Re=1000 时的定常计算结果Fig.3 Steady state atRe=1000

3.2 非定常流动

本文选取了三个计算状态分别揭示周期性流动(Re=1700)、准周期性流动(Re=1735) 和湍流(Re=20 000)的流场特性.

图4 展示了不同雷诺数信息采集点Prm的水平速度分量及通过傅里叶变换之后的速度频谱曲线.可以看出当雷诺数增加至1700 时流动已经演化为非定常周期性流动,此时流场稳定性已被破坏,取而代之的是以频率为f=1/T=0.343 的周期性振荡.当雷诺数进一步增加至1735 时,流动演化为非定常准周期性流动,虽然保留了周期性流动的部分特征,但是此时的流动不在以周期性单一频率振荡演化.从图中可以看到代表准周期性流动典型特征的双频率,其中f1=0.338 是从周期性流动继承而来的基本频率,而f2=0.169 是伴随周期性流动出现的调制频率,其他频峰则是基本频率和调制频率的线性组合.当雷诺数增加至20 000 时,流场已完全被湍流取代,但是从速度频谱图中仍可以观察到之前准周期性流动的双频率,而围绕双频率频峰的基本都是宽带噪音.此时从速度曲线图也可以看到流动演化的无序性和随机性.图5 给出了对应雷诺数的速度相图,其中闭合的单一曲线代表了周期性流动.而准周期性流动的相图不再是单一闭合的曲线,也不是混沌状态的一团乱麻(湍流),而是介于二者之间的过渡状态,其物理特性既继承了周期性流动的部分特征同时预示了可能出现的湍流.

图4 速度频谱图Fig.4 Velocity spectrum

图5 速度相图Fig.5 Velocity phase map

结合图6,展示了周期性流动在一个完整周期内不同时刻的流场拓扑结构.图6(a) 给出了周期性流动完整周期内水平速度随时间的变化曲线以及不同时刻的选取方式.此时,流动仍旧继承了流动定常解拓扑结构的主要特征,可以观察到两个主涡及其附近的次级涡依旧存在,只不过整个流场以f=0.343的频率循环振荡.

图6 完整周期内不同时刻涡量图(Re=1700)Fig.6 Vorticity snapshots at different time steps within a full period T (Re=1700)

图7 描述了准周期性流动(Re=1735)在不同庞加莱交叉点的流场拓扑结构,给出了准周期性流动庞加莱交叉点的选取方式.图中urm代表准周期性解信息采集点Prm处沿x方向的速度分量代表准周期性解调制频率f2对应的长周期在一个整周期内的速度随时间t的曲线.庞加莱交叉点的选取满足urm=-0.003 46 和的条件.图7(b)～图7(d)分别展示了准周期性解在不同时刻(庞加莱交叉点)的瞬时涡量图,虽然整体流动趋势基本一致,但是流动细节有不同的呈现(见图7(d)所示的对应不同时刻的流函数).这是因为此时流场已经演化为准周期性流动,尽管庞加莱交叉点的选取方式一致,但是对应的是准周期性解调制频率f2所对应的长周期在不同时刻的瞬时涡量图.可以想象,如果此时的流动仍为周期性流动,那么这三个时刻对应的计算结果将完全一致.这也进一步证实了流动此时已经演化为准周期性流动的事实.从流场拓扑结构来看,此时的准周期性解与周期性解的差别不是特别明显,需要从流动细节观察区分.因为流动刚从周期性演化至准周期性不久,其准周期性特征不是特别明显,单从流场拓扑结构很难区分.

图7 准周期性解(Re=1735)Fig.7 A quasi-periodic solution (Re=1735)

结合图8 分析了湍流流动状态的流场特性,图8(a)揭示了湍流状态下的能量频谱图,可以看到湍流惯性子区的斜率为 -5/3,跟文献[24-26]的结论一致.伴随能量级串现象的出现,大尺度的涡结构已破碎变成细小的涡结构,而能量也依次传递至小尺度的涡结构直至Kolmogorov 尺度的能量耗散现象.

图8 湍流状态计算结果Fig.8 Results for chaos

图8(b)给出了某一特定时刻的流场拓扑结构,其中实线代表逆时针旋转的涡,虚线代表顺时针旋转的涡,虽然此时流场已演化至湍流,流动已变得随机和无序,但是部分特征仍旧明显,如流场内大尺度的涡结构都发生了破坏,此时主导流场拓扑结构的基本都是小尺度的细碎涡结构.并且每个尺度的涡结构都有对应的振荡频率,如图8(b)所示,频峰f1所对应的频率为主涡的振荡频率,相较于定常结果(见图3),主涡的结构发生了明显破坏,且尺度变小.频峰f2,f3,和f4分别对应不同位置的次级涡的振荡频率.频峰f5和f6对应边角处次级涡的振荡频率.除了这些较大尺度的涡系结构所对应的振荡频率,从能量频谱曲线中还可以观察到其他频峰对应的附着在主涡和次级涡周围的细碎涡的频率.

3.3 Hopf 流动分岔点

随着雷诺数的增加,流动会从定常状态演化至非定常状态,本文针对镜像对称顶盖驱动方腔内流,根据研究不同雷诺数下的扰动衰减系数(Lyapunov指数) ε=ln(urm-U¯),发现流场稳定性最初的破坏伴随Hopf 流动分岔点的出现而发生,这与我们之前研究工作[14-15]中的结论一致.如图9 所示,当雷诺数从1500 增至1692 时,扰动衰减系数的斜率不断增大,从一个负值逐渐趋向于0.当扰动衰减系数的斜率变为0 时,说明此时流动已经演化为非定常周期性流动,意味着Hopf 流动分岔点出现在雷诺数等于1691 和1692 之间.对比经典的顶盖驱动方腔内流[14-15,27-28](ReH=8025±25),我们发现,该流场的稳定性非常差,很容易失稳.说明顶盖镜像驱动对于腔体内流有很强的失稳作用.

图9 扰动衰减系数Fig.9 Perturbation decay rate

3.4 Neimark-Sacker 流动分岔点

随着雷诺数的进一步增加,从1725 增至1735,流动由非定常周期性流动演化为非定常准周期性流动.说明Neimark-Sacker 流动分岔点出现在雷诺数等于1725 和1735 之间.图10(a)和图10(b)分别展示了速度频谱图和相图,其中黑色和红色曲线分别代表雷诺数为1725 和1735 的计算结果.结合图10(a),可以观察到周期性解的振荡频率为f=0.338,其他频峰均为f的整数倍,为基本频率的谐振频率.而准周期性解则有两个频率,分别为基本频率f1=0.333,调制频率f2=0.169,其他频峰则是f1和f2的线性组合,如f3=f1+f2=0.502 ,f4=f1+2f2=0.671.如图10(b)所示,周期性解的相图是一个闭合的单一曲线,而准周期性解的相图则由一组曲线族构成.类似顶盖驱动和顶底双边驱动方腔内流[14-15],研究发现流动非定常周期性的破坏通常都是伴随Neimark-Sacker 流动分岔点的出现而发生.没有发现类似三角腔顶盖驱动和四边驱动方腔内流时所出现的Period-doubling流动分岔点.并且类似其他腔体内流流动,流动的非定常周期性并不能长期保持,随着雷诺数增加,伴随准周期性流动典型双环形结构的出现,流动很快会进一步演化为准周期性流动.

图10 周期性和准周期性计算结果对比Fig.10 Comparison between periodic and quasi-periodic solutions

3.5 湍流始现

经历了Neimark-Sacker 流动分岔点之后,当雷雷诺数增加至1750 时,流场演化逐渐表现出无序性和随机性,说明湍流出现在雷诺数等于1735 和1750之间.图11(a)和图11(b)分别展示了速度频谱图和速度相图,其中黑色和红色曲线分别代表雷诺数为1735和1750 的计算结果.从速度频谱图(图11(a))可以明显观察到准周期性解的两个频率,分别为基本频率f1=0.333,调制频率f2=0.169.当雷诺数增加至1750 时,虽然仍旧可以观察到从准周期性流动中继承来的两个振荡频率,但是被一系列宽频噪音所包围,说明此时流动已然变为湍流.如图11(b)所示,准周期性解的速度相图是闭合的曲线族(图11(b)局部放大图),而湍流的速度相图明显是一个无序、随机的混乱系统,更进一步证实此时流动特性主要表现为湍流.至此,流动的演化路径已基本明晰,这与顶盖和顶底双边驱动内流的研究结论保持一致,即流动先从定常演化为非定常周期性,再演化为准周期性流动,最终演化为湍流.

图11 准周期性和湍流计算结果对比Fig.11 Comparison between quasi-periodic and chaotic solutions

3.6 流场镜像对称性

由于本文的驱动条件为严格的顶盖镜像对称,所以流动在初期表现出了严格的对称性,且对称性参数一直为0(见图3).但是当我们观察对称性参数(图12)时,可以看到当Hopf 流动分岔点出现时,对称性参数不再是0,且随着雷诺数增大而逐渐增加.这就说明当Hopf 流动分岔点出现时,伴随着流场稳定性的破坏,流场镜像对称性也发生了破坏,这与我们之前研究[15]中观察到的结论类似,即流场 π 旋转对称性与流场稳定性同时丧失.同时对比与该流动较为类似的Taylor-Couette 流动,我们发现了相同的结论,镜像对称性的破坏往往伴随着流场稳定性的丧失,流动不可能不经历对称性破坏而直接演化为湍流[29-30].

图12 不同雷诺数的对称性参数Fig.12 Symmetry at different Re

虽然此时对称参数在数量级上非常小,但是足以说明流场对称性在逐渐丧失,这样的微小差别在肉眼观察流场拓扑结构时很难发现.但随着雷诺数进一步增加,流场非对称性就可以直观地显现出来(见图6).图13 展示了流场演化之标准周期性流动时的速度曲线和对称性参数曲线.可以清晰地看到此时的对称参数以0.015 左右的振幅周期性振荡.

图13 水平速度分量及对称性参数(Re=1700)Fig.13 Velocity and symmetry series (Re=1700)

3.7 流动滞后

在本文之前的研究工作中,我们发现对于顶盖镜像对称驱动方腔内流这一特定流场,Hopf 流动分岔点出现在雷诺数等于1691 和1692 之间.这个结论从图14 中也可进一步证实,当雷诺数增加至1700 时,流动从定常状态演化至非定常周期性流动.如图所示,红色符号代表由初始状态计算得到的结果,而黑色符号意味着,首先基于初始状态计算得到周期性结果(Re=1700),然后在此基础上计算雷诺数小于1700 的流动.其中×代表定常结果,△,□ 和 ○ 分别代表非定常周期性流动的最小值、平均值和最大值.

图14 流动滞后现象Fig.14 Flow hysteresis

如图所示,基于Re=1700 的周期性结果将雷诺数分别降至1600,1500 和1400,流动仍然保持周期性特性,并非之前观察到的定常结果(由初始态计算得到),说明存在流动滞后(flow hysteresis)现象.进一步降低雷诺数至1350,此时流动才回落至定常状态,说明流动滞后现象发生在 1350 ＜Re＜1700 这个区间.同时也说明此前捕捉到的Hopf 流动分岔点为亚临界(subcritical)形式.流动滞后现象的出现意味着在 1350 ＜Re＜1700 这个区间,对应的每个雷诺数会有两种解的可能性,一种是定常状态,另一种是周期性状态.并且,根据观察得到,此周期性流动的基本规律与之前基于初始状态计算得到的周期性解特性基本一致.如图15 所示,展示了Re=1600 时的周期性解在一个完整周期内不同时刻的涡量图.