槽流的亚临界转捩与局地湍流
2020-03-13陶建军
陶建军
(北京大学 工学院 力学与工程科学系, 湍流与复杂系统国家重点实验室, 北京 100871)
0 引 言
转捩是指流动从层流态转变为湍流态的过程，是生产生活中的常见现象。由于层流和湍流的动量、能量、质量传输特性有着很大不同，在工程设计中为了减阻、隔热、降噪需要流场尽量保持在层流状态，而在增混、促燃和强化换热时又需要促进转捩过程。因此，航空航天航海以及能源化工设备的优化设计同层流-湍流转捩的高效控制是密不可分的。转捩研究可追溯至雷诺的管流实验，而转捩的起始与流场的失稳密不可分[1]。亥姆霍兹、凯尔文、瑞利等人的早期稳定性研究专注于流场对无穷小扰动的响应，通过线性化的控制方程求解不稳定的简正模态(normal mode),从而确定临界控制参数。黏性剪切流最重要的控制参数即为雷诺数。对于小幅值扰动，一些流动如瑞利-贝纳尔对流随着控制参数的增加可经历一系列的超临界分岔过程，分岔的临界参数是可以精确确定的[2]。但对于边界层流、黏性剪切驱动的槽流或压差驱动的管流和槽流，其转捩行为可明显不同。它们的层流速度型在中等雷诺数时对小幅值扰动(线性扰动)是稳定的[3-4]，但扰动幅值足够大时可变为湍流，即发生亚临界转捩[5-6]。以管流为例，对于不同形式和强度的扰动，其转捩临界雷诺数的变化范围可很大，如1700～100 000[7-8]。剪切流的亚临界转捩过程对外部扰动的敏感性给转捩预测带来极大困难。
在边界层流中，人们发现一对斜Tollmien-Schlichting(TS)波的非线性演化对边界层流的转捩至关重要[9-10]；基于对转捩过程中湍流间歇性的认识，构造了基于间歇因子输运方程和当地流动变量的转捩模式[11]；实验研究了壁湍流转捩中的三维非线性波和二次涡环等结构[12]。对于光滑壁的平行剪切流，如管流，人们发现初始人为引入的大幅值流向涡扰动并不会完全衰减掉[13]，随后的实验研究显示局地湍流(puff)[14-18]的瞬态增长特性和分裂对管流的亚临界转捩过程有重要意义。Eckhardt等和Mullin分别于2007年[19]和2011年[20]在流体力学年鉴(Annu. Rev. Fluid Mech.)中对管流亚临界转捩的研究进展进行了综述，不断涌现的成果相继发表在了Science[21-23]和Nature[24-25]，引起广泛关注。本文将重点讨论平行平板间压差驱动的平板泊肃叶流 (plane-Poiseuille flow，PPF)的亚临界转捩过程，同时简述剪切驱动的平板库埃特流(plane-Couette flow，PCF)和压差驱动的圆管流动(Hagen-Poiseuille flow，管流)的相关研究进展。
1 槽流的亚临界转捩现象
能量稳定性分析或能量法可以给出基本流失稳的必要条件。对于槽流这类黏性剪切流而言，能量法所定义的临界阈值远远低于实际的转捩雷诺数[26-28]。由于线性稳定性分析所得特征函数并不一定是两两正交的，初始扰动有可能经历瞬态增长[3,29]，因此人们寄希望于经瞬态增长机制将小扰动放大至有限幅值后的非线性阶段，也研究了槽流的最大瞬态增长模态和三维扰动瞬态增长的最大能量放大倍数[30-32]。需指出的是，由于扰动能量的增长机制是线性的，能量法和针对扰动方程线性算子的瞬态增长分析均与扰动的初始幅值无关。采用同伦分析的方法[33-36]、在小尺度周期盒子中用迭代法[37-38]或通过二分的边界追踪(edge-tracking)[39-40]等数值方法，人们可求得空间上有严格周期性的非线性解和有局地特性的边界态(localized edge state)。尽管这些Navier-Stokes方程的解在中低雷诺数时会衰减掉，其流动结构、动力学和统计学特性对人们认识低雷诺数湍流是有益的[41]。近年来的一个重要研究进展是发现倾斜的大尺度的湍流带同槽流的亚临界转捩进程直接相关。
对于平板泊肃叶流(PPF)，常用的雷诺数定义为Re=3Umh/(2ν)，其中h为半槽高，Um为截面平均流速，ν为运动学黏性系数。1928年，Davies和White实验测量了不同横截面宽高比的槽道流的壁面摩擦阻力系数与雷诺数的关系[42]。实验结果显示，摩阻系数在小雷诺数时同层流解相符，在Re高于某临界值时开始高于层流解。更大的宽高比意味着更小的侧壁影响，数据显示临界雷诺数随着宽高比的增大而增大，当宽高比足够大后则不再增加，维持在Re=667.5附近。这一临界值在此前一直未受重视，因为随后实验发现的转捩临界值要高得多。1969年，Patel和Head[43]的槽流实验发现当Re<1035时，流动维持在层流状态；当1035
Carlson等人通过流动显示[45]发现，雷诺数在1000左右时壁面小孔的脉冲射流可激发出湍斑结构(见图1a)，在往下游的发展过程中湍斑会发生分裂，如图1(b)所示。当Re大于1200时，无需人工激励即可由背景扰动随机地形成湍斑。在Re<840时，激励触发的湍斑无法充分发展，最终衰减掉。由于这些实验结果显示可自维持的湍流要求雷诺数高于1000，随后大量的实验、理论和数值模拟工作主要集中在了Re>1000的湍斑研究上[46-50]。Alavyoon等人在实验中也采用壁面小孔射流来激发湍斑[51]，他们发现在Re=1100附近可以观察到未发展完全的湍斑；当Re<1100时，任意大幅值的扰动都无法产生湍斑；直到Re>1200时，才能观察到完整的湍斑，即在展向有对称的两个分支的湍斑。随着计算能力的提高，人们在Re≥1070的数值模拟中发现存在规则排列的被层流区间隔的湍流带(turbulent stripe或turbulent band)结构[52,53]，随后的实验也证实了湍流带结构的存在[54](见图2)。实验显示在Re<1275时湍流带会断裂；在Re=1200时仅发现一两个湍斑；在Re=975时流场为稳定的层流[55]。以Re≥1567的层流场叠加扰动涡对作为初始场，分别在保持净压力梯度不变[56]和保持流量不变[57]的条件下进行的数值模拟发现，初始扰动诱发的湍斑均会扩展、演化。
(a) x/h=64, Re=1000
(b) x/h=260, Re=1000
(a) 数值模拟所得的流向扰动速度等值面, δ为半槽高, Re=1440
(b) 流动显示，Re=1500
Tao等人重新分析了Davies和White的实验数据[42]并将其分为三个区(见图3)[58]。I区对应于层流，II区和III区的斜率明显不同。II区随着槽道高宽比的减小而减小，而III区的斜率并不随高宽比变化，反映了在该区间不受高宽比影响的内在的动力学特性。II、III区的分界线外插至层流态的临界值估计为Red=1127或Re=845。采用大计算域的初步的数值模拟也显示在雷诺数低至900时存在增长的湍斑，因此该阈值被解释为可发展湍斑的临界值[58]。这一阈值的提出提示人们在Re<1000的情况下，PPF有可能存在着可自维持的湍流态。
图3 对Davies和White实验数据[42]的分析[58]Fig.3 Analysis of the experimental data[42,58]
Tuckerman对PPF用倾斜的狭小计算域进行了数值模拟[59]，发现在雷诺数大于850时存在湍流带。需注意的是，采用狭小计算域虽可大为节省计算资源并定性模拟局地湍流结构，但人为设定的计算域倾角和展向的周期边界条件对转捩过程及临界阈值的确定有显著影响[60]。为减少周期边界条件对局地湍流结构的影响，2013年Tao等采用流向和展向尺度均很大的计算域数值模拟了PPF的亚临界转捩过程，发现存在湍流带的雷诺数可低至720[61]；得到了孤立的湍流带(见图4)，即在计算域内湍流带周围为层流，没有其他的湍流结构。由于孤立的湍流带远离其他局地湍流结构从而大大减弱了同它们的相互影响，其存在为研究湍流带自身的时空演化特性提供了可能。
根据数值模拟结果，存在统计定常态的转捩阶段被定义为平衡局地湍流态(equilibrium localized turbulence)[61]。根据全流场的耗散率特性定义的平衡局地湍流态的临界值为Re=859，该阈值同之前分析Davies和White实验结果所得阈值845[58]基本一致。数值模拟还发现在Re>950时湍流带发生分裂的频次增加从而大为促进湍流的扩展，由层流区间隔开的湍流带结构最终会铺满整个流场；在Re>1300时层流区不再明显，此时流态被定义为均匀湍流(uniform turbulence)[61]，如图5所示。随后，人们开始重视PPF在雷诺数低于1000的转捩特性[62-63]。
(a)
(b)
(c)
(d)
通过对PPF在大计算域(流向和展向计算域尺度均为102h量级)进行系统的直接数值模拟，Xiong等人[64]发现，局地扰动诱发的湍斑要么演化为湍流带，要么衰减掉，因此其自维持时间远远小于湍流带的自维持时间；在中等雷诺数时，孤立的湍流带会倾斜伸长(见图6)；当Re<660时湍流带扰动动能经瞬态增长后会衰减掉，见图6(a)(b)。另外，湍流带周围存在大尺度环流(见图6c)，这些环流之间的相互作用会导致湍流带的断裂和衰减，而在适当的条件下断裂的湍流带还会发生重连现象[65]。倾斜狭长的计算域无法模拟湍流带的伸长现象，而小的计算域由边界条件带来的较强的环流相互作用会导致湍流带的断裂和衰减，从而无法准确获得可自维持湍流带的下临界雷诺数。比如近期人们发现用小展向尺度的计算域进行的数值模拟只能得到最终衰减的类似于管流中puff的局地湍流结构[66]。综上所述，在中低雷诺数发现的湍流带结构使人们对槽流亚临界转捩的认识进入了一个新阶段。
(a)
(b)
(c)
(d)
(a)
(b)
(c)
2 亚临界转捩的理论模型
为了理解亚临界转捩中局地湍流结构的生成与演化机理，人们建立了多种动力学模型。Barkley等借用化学中的反应扩散模型和生命科学中的神经元模型来定性描述管流从puff的瞬态增长到最后全场湍流的整个亚临界转捩过程[67-70]。金融领域的极值理论(Extreme Value Theory)和生态捕食模型(predator-prey model)也分别被用来分析平板库艾特流亚临界转捩阈值问题[71]和管流转捩现象[72]。Pomeau在1986年将转捩过程与定向逾渗(Directed Percolation，DP)的相变类型联系了起来[73]。在定向逾渗理论中，一个活跃状态可以通过相邻的节点以一定概率p进行传播。随后人们对管流和平板库艾特流转捩过程的实验和数值模拟数据用DP模型进行了拟合分析[74-76]。
对于PPF，2016年Sano和Tamai通过增加多孔板的方式提高试验段入口流场的湍流度，并利用流动显示技术观测下游流动结构的变化，计算了流场中湍流占比与Reynolds数的关系[77]，如图7所示，发现湍流占比(湍流区域所占的面积百分比)与参数ε= (Re-Rec)/Rec呈幂次分布，同DP模型相符，并通过外插法定义了PPF亚临界转捩的临界雷诺数为830。但从图7可以看出，在Re<830时湍流占比并不是零，即转捩已经开始，这是DP模型无法解释的。实际上，用DP幂次关系外插实验数据所得的转捩阈值同之前依据早期实验数据外插所得的阈值845[58]和由数值模拟数据定义的平衡局地湍流态的临界值859[61]是符合的，反映的是存在统计定常态的局地湍流的临界性质。在雷诺数低于830时，之前的数值模拟已经表明湍流带是可以增长演化的(如图6)。近期在Re= 660～800范围内的数值模拟进一步表明，大计算域中的周期湍流带是可以自维持的，这一阶段每一雷诺数对应的湍流占比不唯一但存在一个统计上的最大值，即PPF在转捩之初存在一种稀疏湍流态[65]。
图7 实验所得PPF湍流占比随雷诺数的变化[77]Fig.7 Turbulent fraction obtained experimentally for PPF at different Reynolds numbers[77]
3 槽流转捩的相似性
无论对管流还是槽流，如平板库艾特流和平板泊肃叶流，局地湍流结构的产生与时空演化特征显示亚临界转捩过程可包括不同的阶段(见图8)，而不同黏性剪切流的转捩阶段之间的相似性似乎预示它们遵循着相同的转捩路径[78-79]。但对这种相似性的量化分析是困难的，比如不同流动发生转捩的雷诺数可相差巨大，因为雷诺数的定义中选取的特征长度和特征速度因流动不同而不同。那么不同的黏性剪切流是否可定义统一的稳定性控制参数呢？对这一问题的研究已有六十年的历史。
图8 平行剪切流的亚临界转捩阶段[79]Fig.8 Transition stages of the subcritical transition of parallel shear flows[79]
人们分别尝试依据平均流特性[80-81]和流场的局部特性来定义有普适性的稳定性参数。1959年，Ryan 和 Johnson[82]依据线性化的二维能量方程为管流提出了一个稳定性参数Zmax，是“流体微团的能量输入与能量耗散之比的函数”。Zmax代表着横截面内该能量比值的最大值。管流实验数据显示牛顿流体和几种假塑性流体的转捩临界Zmax值是一致的，这在一定程度上证实了一个观念：流态是由具有最大能量输入与能量耗散之比的流体微团或局地稳定性参数Zmax决定的。由于Zmax不依赖于边界条件，仅取决于流场的局地特性，因此相较于传统的雷诺数定义Zmax更能反映不同流场转捩过程的共性。随后Hanks[83]基于Navier-Stokes (NS)方程将这一观念扩展到了槽流(PPF)，并给出了另一个局地控制参数K，其对压力驱动管流和槽流的定义分别为：
(1)
普适的转捩或稳定性控制参数对不同的流动应具有统一的定义式。依据NS方程定义的局地雷诺数Rem反映了局地的非线性对流与黏性扩散效应之比[58]。与之前的稳定性参数定义中的分母(dp/dx或μ2U)仅由基本流决定的情况不同，Rem定义中的基于物理模型所表征的黏性扩散项为包含了扰动成分。η是受扰流线的局地斜率，l是扰动的流向特征长度尺度。对于舒缓的扰动，η/l≪1，扩散项的第二项相较于第一项(如在泊肃叶流中)可忽略；对于平板库艾特流，其线性速度剖面使得扩散项中第一项为零，则包含η/l的第二项成为主导项。 受Ryan和Johnson工作的启发，Tao等[79]依据动能方程定义了另一局地雷诺数ReL，用来表征受扰流体微团从基本流所获能量与其能量耗散之比在截面内的最大值，其定义式是：
(2)
Rh是水力半径(截面积与湿周之比),τw是壁面切应力，x2是壁面法向坐标， “max” 代表截面内最大值。传统雷诺数的定义式包含的是反映流场平均特征的量，如截面平均流速和槽高，但这些平均量难以反映具有不同边界条件的流场特性(如PCF和PPF)。尽管不同形式的流动有不同的边界条件，但流体微团的动力学行为是统一的，其动量与能量控制方程是相同的，因此反映流体微团能量输入与耗散特性的局地雷诺数ReL相较于传统的雷诺数定义具有更好的普适性。另外，不同于之前的局地控制参数 Zmax和K，ReL既适用于压差驱动的泊肃叶流，又适用于黏性剪切驱动的库艾特流，而且具有统一的定义形式。将管流、平板库艾特流和平板泊肃叶流的基本流代入公式(2)所得的各类流动的ReL表达式同Rem的表达式是一致的；管流和平板泊肃叶流的ReL表达式同Hanks定义的K也是一致的[79]。
通过分析管流、PPF和PCF以往的实验和数值模拟结果[79]并整理近期结果[77,86]，三种不同剪切流亚临界转捩的不同阶段的ReL阈值基本相符，分别为：255～290(稀疏湍流)、320～340(平衡局地湍流)、380～394(空间间歇性分布但可自维持湍流)、500～520(均匀湍流)。采用传统雷诺数定义的阈值则相差巨大，比如对于管流和平板库艾特流，其平衡局地湍流的阈值分别为1750和325左右。显然，采用基于流场局地信息定义的稳定性参数是可以定量表征不同剪切流亚临界转捩过程的相似性的。当然，需要提及的是不同流动的转捩过程也会有其特殊性。比如在转捩的最初阶段，管流中的puff结构是瞬态的，平板泊肃叶流中的湍流带结构在较狭窄的流道或计算域中也是瞬态的[66]，但当这些局地湍流结构相距很远时是可以维持倾斜伸长的。因此，亚临界转捩的最初阶段被称为稀疏湍流态。
4 结论与展望
剪切流的亚临界转捩是一个剧烈的过程，但人们通过研究局地湍流结构的产生、演化(如湍流带的断裂、伸长、分裂、重连)，认识到该过程仍可细分为不同的阶段，而不同流动的转捩阶段是存在相似性的。这种内在的相似性来自于流体微团受控于相同的动力学规律(方程)这一事实，而基于局地流场特性所定义的稳定性参数揭示了不同流动在不同的转捩阶段具有相符合的临界阈值，从而定量地反映了黏性剪切流亚临界转捩的统一性。尽管人们已经取得了许多可喜的进展，但对于剪切流的亚临界转捩过程仍有许多重要问题有待解决。首先，对于槽道流亚临界转捩的实验研究仍是不充分的，目前的认识多是来自于大规模的并行计算，局地湍流结构在各转捩阶段的临界阈值仍需实验确认，局地结构的演化特征(如湍流带的倾斜伸长、分裂、重连)仍需相应的实验研究[87]。其次，局地湍流结构的产生和演化机理仍需更为深入的探索。目前的动力学模型分析多是从类比、形似的定性角度来解释现象，仍需直接根植于流动控制方程的机理研究[88]。另外，目前的亚临界转捩研究多是由强初始扰动触发的，由复杂边界触发的亚临界转捩过程的相关研究仍有待深入。相信随着实验技术、计算能力、数值方法和稳定性理论的发展，这些挑战性的问题在不久的将来会被一一攻克，从而进一步优化人们的转捩控制方法、加强人们在工程实践中的流动控制能力。
