界面张力梯度驱动对流向湍流转捩的研究

郭子漪 李 凯 康 琦 段 俐 胡文瑞

中国科学院力学研究所微重力重点实验室,北京100190

中国科学院大学工程科学学院,北京100049

1 引 言

自然对流的流动规律是流体力学的首要研究问题之一.自然对流是指流体在无外界直接驱动力 (如压力) 情况下自发产生的流动.浮力对流是体积力引起的一类重要的自然对流形式,其内在驱动力来源于重力效应下流体内部温度差或组分浓度差等形成的密度梯度,流体力学中通常用Rayleigh 数来表征温度变化产生的浮力效应与流体黏性及热扩散效应之比,Rayleigh 数定义为Ra=gβ∆Td3/νκ,其中g为重力加速度,β为热膨胀系数,∆T为温度差,d为特征长度,ν为运动黏度,κ为热扩散系数; 另一类重要的自然对流形式是流体界面力引起的流动,流体界面(自由面)给出流体所占物理区域的边界,包括液/固界面、液/气界面及分割不混溶液体的液/液界面.这一过程中流体界面张力及其梯度对界面位形及其形变起着至关重要的作用.流体界面张力大小与界面上温度和组分浓度等相关.对于大多数流体界面,界面张力随温度升高而减小.当流体界面上存在温度梯度时,界面热点处界面张力较小,冷点处界面张力较大,由此在界面上形成了从热点到冷点的切向力,驱动流体沿界面运动;流体内部则在流体黏性的作用下形成回流,由此引起的流体整体流动称为界面张力梯度驱动对流(胡文瑞,徐硕昌1999,Walter 1987).由温度梯度引起的界面张力梯度驱动对流又被称为热毛细对流,其中常用的无量纲参数为Marangoni 数,反应了热毛细对流中对流输运与热扩散之比,定义为Ma=∆Td/ρνκ,其中为界面张力温度系数,ρ为流体密度.具有流体界面的流体体系中,浮力对流和界面张力梯度驱动对流常常是耦合的,流体力学中通常用无量纲静态或动态Bond 数(Bo或Bd) 来表征重力(或浮力) 效应与界面张力效应的相对重要性,动态 Bond 数定义为Bd= ∆ρgβd2.当动态Bond 数较大时,重力效应在体系内占优,体系内浮力对流占主导; 当Bond 数较小时,即对于小尺度(特征尺度d小)或微重力环境下(g小)的体系,界面张力起主导作用.对于具有流体界面的流体体系,地面环境下,界面张力梯度驱动对流通常作为次级作用为浮力对流所掩盖;而在微重力环境下,重力效应被极大地减弱,此时界面张力梯度驱动对流成为自然对流热质输运的基本形式(Johns and Narayanan 2002).例如,前期的空间晶体生长探索实验中,曾预期空间微重力环境下,重力效应极大减弱,生长过程为纯扩散控制,可以生长高品质大尺寸晶体.然而实验结果发现,微重力环境下界面张力梯度驱动对流占主导地位,其流动强度甚至可以比拟地面生长过程中的浮力对流.界面张力梯度驱动对流同样影响空间生长晶体的质量,当其驱动力大小超过临界值时,对流变为周期性振荡的非定常对流; 随着驱动力的不断增大,对流最终发展为湍流.对于半导体硅这样低Prandtl 数(Pr=ν/κ) 的熔体,非定常对流极易发生(Croll et al.1994),非定常的对流热输运过程引起温度场的脉动,使得生长速率的脉动与对流过程间接耦合起来,产生与地面生长过程中的浮力对流类似的杂质条纹.

混沌是自然界所固有的现象.Poincare (1899) 于19 世纪末研究三体问题时指出确定的动力学方程的某些解存在不可预见性,从而开创了混沌理论.Lorentz (1963)于1963 年根据前期研究结果开创了用数值实验方法研究混沌现象的方法.自此混沌现象研究因其在理论和实际应用中的重要性受到了持之以恒的关注并获得了突飞猛进的发展.在确定性系统中,转捩过程是连接“确定” 与“随机” 之间的桥梁,是一个强非线性过程.流体运动时空演化过程中的转捩现象最早由Reynolds通过实验观测发现(唐登斌2015): 圆管中的水流流动在超过某个临界状态时会从层流转变为湍流.流动从层流向湍流的转捩过程主要分为两个阶段: 首先,随着流动驱动力增大,流动中的小扰动线性发展,超过临界状态时,稳态层流失稳发展成为非稳态层流,这个过程可以视为临界转捩;随着流动驱动力继续增大,流动中的非线性作用不断增强,大量谐波不断产生,扰动成为随机,导致非稳态层流进一步发展成为湍流,这个过程可以视为超临界转捩.边界层理论(Prandtl 1904)表明层流和完全湍流间存在过渡区域(转捩区).与直接研究极其复杂湍流的内部结构相比,从转捩过程的角度去探索层流如何转化为湍流为研究复杂的湍流问题提供了“钥匙” 和新视角,因而成为相关流体耗散系统中流体力学的前沿研究.

作为空间自然对流热质输运的基本形式,界面张力梯度驱动对流是流动和传热强耦合的复杂非线性过程,同时也是一个多控制参数(包括体系的几何位形、物性参数、界面热边界条件等) 耦合作用的过程,因而呈现出非常多样性的流动时空特征.例如,对于同一个模型问题,在其他参数不变的情况下,随着界面上温度梯度增大,界面张力梯度驱动对流会呈现丰富多样的流型,稳态和非稳态、层流和湍流、周期振荡、流动分岔和混沌等.因此,界面张力梯度驱动对流时空转捩过程、转捩机制、非线性特征及流动向湍流转捩途径等基本规律成为微重力流体物理的重要研究内容和学科前沿.深入系统地开展相关内容的研究将进一步推动小扰动理论、能量法、弱非线性理论等流动稳定性理论的发展,促进混沌理论、耗散结构理论等新兴科学的发展及其与微重力流体力学的融合,极大地丰富微重力流体力学的研究内容.

与此同时,界面张力梯度驱动对流在载人航天、登月计划、火星探索和空间流体利用等诸方面都有较强的应用背景,涉及诸如空间燃料输运过程和空间能源或热管利用(Abe et al.2005) 等空间流体管理问题、空间生命保障系统、空间材料制造、空间焊接和空间生物技术等.因此开展相关流动基本规律的研究对于人类认识、探索和利用空间环境也具有重要的应用价值,同时对于指导和提高地面相关工业过程等都有重要意义.

2 界面张力梯度驱动对流向湍流的转捩

2.1 界面张力梯度驱动对流的临界转捩

液层模型是源自于空间晶体生长实验探索的简化模型,是研究微重力条件下界面张力梯度驱动对流基本流动规律的经典模型.液层模型主要包括矩形液层 (rectangular pool)、环形液层(annular pool)和液桥(liquid bridge)(图1)等.相对而言,液桥是一类特殊的液层模型.液桥中界面张力梯度驱动对流的流动稳定性研究成果极为丰富,相关系统的完整临界转捩描述详见综述(Hu et al.2008).除了极薄液层的情况,在绝大多数界面张力梯度驱动对流的地面实验研究中,界面张力效应和重力效应间存在较强的耦合作用.例如,Riley 和Neitzel (1998) 的地面实验研究结果表明,液层中热流体波的特征,如传播角度和振荡频率,随着液层厚度的不同而改变,表现出与重力效应的强耦合性.这与Chan 和Chen(2010)的研究结论一致.而另一方面,由于空间实验机会昂贵稀缺,界面张力梯度驱动对流的空间实验研究屈指可数.Kamotani 等(1995,1999) 开展了圆形液池中界面张力梯度驱动对流临界转捩的空间实验研究,实验采用的液层厚度相对较厚,高径比约为1.实验研究了层流失稳的临界条件和临界流动振荡模式,讨论了其与加热速率、自由面位形和激光束加热面积的关系,并和相应的数值模拟结果进行了比较.环形液层的特点是周向没有壁面限制,热流体波不会因为壁面影响而衰减(Lappa 2009),可以更好地观测热流体波的基本参数,例如临界波数的变化.Schwabe 等(1999,2002) 开展了外加热环形浅液池中界面张力梯度驱动对流的临界转捩空间实验研究.实验结果表明临界值与液层尺度比无关; 实验中还确认了热流体波的存在,当液层高度很小时,随Marangoni 数的增加,流动由同心多卷圈结构发展为热流体波; 而当液层厚度较大时,失稳后的流动比热流体波更为复杂和不规则;此外,实验中还发现特殊的轴对称临界振荡.但是实验结果分析中没有考虑自由面位形的影响,部分实验结果与二维矩形液层的数值模拟结果有较大差别.近期,Jiang (2017a,2017b)通过实验研究了矩形液层热毛细对流的转捩问题,发现了多种转捩途径,并发现存在表面波动不稳定性.Kang 等(2019a,2019b,2019) 在实践10号返回式科学卫星上开展了内加热环形液池中界面张力梯度驱动对流的临界转捩空间实验研究,实验结果给出了层流失稳的临界条件,临界振荡流动模式及振荡频率等,系统讨论了上述临界转捩特征与液层体积比(自由面位形)关系.这里液层体积比定义为液层实际体积与对应的水平液层的体积之比.此外,Kang 等(2019d,2020) 在天宫二号(TG-2) 空间实验完成了液桥热毛细对流空间实验,讨论液桥高径比– 体积比的几何参数效应、多次转捩、波动模式变换、以及分岔道路的复杂性.

图1

相比有限的空间实验结果,微重力条件下界面张力梯度驱动对流临界转捩的理论和数值模拟研究的成果更为丰富.Pearson (1958) 重新分析了B´enard (1901) 的实验,通过在边界条件中加入界面张力项,动量方程中忽略重力项,建立了纯界面张力梯度作用下的数学物理模型,由理论分析给出了流动由静止状态到形成胞状对流的临界条件,与B´enard 的实验结果较为符合.Pearson 的研究开启了一系列对于液层界面张力梯度驱动对流的研究.Smith 和Davis(1983)利用线性稳定性方法研究了微重力条件下无限长的液层,其自由面上施加均匀温度梯度,研究发现,在自由面不变形假定下,除了Pearson 指出的定常卷圈结构外还有热流体波的失稳形式; 当考虑自由面变形时,表面波不稳定性的形式是行波.相关临界值、传播方向和临界波数等受流体物性和界面传热的影响.此外,实际液层的有限边界对热流体波传播有明显限制,进而影响其稳定性.Xu 和Davis(1984)研究了轴向施加均匀温度梯度的无限长流体圆柱,同样发现了与轴向成一定角度传播的热流体波.Smith 和Davis (1983,1986) 对热流体波不稳定性的物理机理进行了讨论,指出热流体波是一种与热效应相关的波动,流场主要起对流输运的作用,即使在忽略表面变形时也会产生热流体波.另一种表面波不稳定性则是表面波动变形与内部剪切流动耦合作用的结果,是纯粹的流体动力学效应.Derby 和Brown(1986)最早在微重力科学计划支持下利用环形液层模型开展了提拉法空间晶体生长的模型化研究.Laure 等(1990) 通过线性稳定性理论和局部分岔理论分析了矩形液池界面张力梯度驱动对流的扰动的空间分布等问题.Xu 和Zebib (1998) 对大Prandtl 数界面张力梯度驱动对流进行了二维和三维数值模拟,对于二维模型,研究给出了较为完整的不同Prandtl 数与高径比下流动失稳的临界 Reynolds 数 (Re=γ∆TH/µν,其中γ=−∂σ/∂T,H为特征长度),并从能量角度分析了振荡流产生的机制;对三维模型,研究给出了典型的Prandtl 数下不同高径比时流动失稳的临界Reynolds 数;同时,研究指出侧壁可以起到抑制失稳的作用,而第三个维度方向上的扰动则会促进失稳.Madruga 等(2003,2004) 对两层流体在水平温差作用下的界面张力梯度驱动对流进行了线性稳定性分析,发现了三种不稳定形式: 从冷端向热端传播的行波、从热端向冷端传播的行波,以及纵向卷圈.P´erezgarc´ıa 等(2004) 理论分析了硅油液层厚度对界面张力梯度驱动对流线性稳定性的影响,指出对流会从基本流失稳为热流体波或者纵向涡胞,是由不同的液层厚度所决定的.Li 等(2011a,2011b) 采用渐近分析方法研究了单层和双层环形液层内的定常轴对称界面张力梯度驱动对流.Shi 等(2006,2010) 开展了环形浅液层界面张力梯度驱动对流的线性稳定性研究,分析了浮力效应、旋转等对稳定性的影响,研究表明,随着Marangoni 数增加,流动由单一传播方向和波数的热流体波发展为各种传播方向和波数的热流体波的叠加状态.Sim 等(2003,2004) 采用三维直接数值模拟方法研究了环形液层内界面张力梯度驱动对流从定常轴对称到三维流动的转捩,发现了在不同尺度比 (厚度与半径之比) 下,液层失稳后出现沿周向的行波和驻波,并研究了界面换热的影响.Li 等(2003,2004) 对不同深度外加热环形浅液层中的小Prandtl 数热毛细对流和浮力– 热毛细对流进行了三维数值模拟,考虑了底面绝热和底面恒定热流的情况,研究发现表面温度场呈现的不同振荡形式: 周向传播的热流体波、受到径向扰动调制的热流体波、驻波等.石万元等(2009)发现环形浅液池临界转捩后形成对数螺线形波纹的热流体波,相应传播角为常数并随着驱动力增加而增大.Tang 和Hu (2007) 数值模拟研究微重力环境下矩形浅液池内的热流体波,对热流体波的形成机理进行了分析和讨论.胡文瑞等(2010) 对矩形液池中的界面张力梯度驱动对流起振过程进行了数值模拟,发现微重力条件下液池两侧的温差超过临界温差时,液池中就会出现振荡对流,其涨落值远小于时间平均值,表现为从冷端向热端传播的热流体波.Li 等(2012)对微重力环境下Marangoni 对流和热毛细对流耦合的液层进行了数值模拟,预测了耦合对流的多涡流结构和临界稳定边界,并报道了振荡耦合对流.Ma 和Bothe (2011) 采用基于VOF 方法的直接数值模拟方法研究了动态自由面形变对液层界面张力梯度驱动对流的影响.S´aenz 等(2013) 研究了浅矩形液池内界面张力驱动对流热流体波及其动态自由面形变对不稳定性的影响.

值得指出的是前述对于界面张力梯度驱动对流临界转捩的理论和数值模拟研究大多针对水平自由面模型开展.空间环境下自由面的形状通常是弯曲的,自由面位形对于界面张力梯度驱动对流及其稳定性有重要的影响(Garnier & Normand 2001,Ma & Bothe 2011,Saenz et al.2013).从流体力学基础研究和空间应用两方面来说,开展弯曲自由面液层体积效应对流体流动的临界转捩特征的影响的研究都是非常必要的.

2.2 界面张力梯度驱动对流向湍流转捩的实验研究和直接数值模拟

临界转捩过程处于整个流动向湍流转捩过程的起始阶段.后续由非稳态流动向湍流的超临界转捩过程则表现出强非线性,其流动结构转捩基本规律(分岔途径)的多样性极大丰富了流体力学的内容,同时其复杂性也给相关理论和实验研究带来了巨大的挑战.分岔现象是系统内部本身固有的特性,是指非线性动力系统某一控制参数连续变化达到临界值(分岔点) 时,系统行为突然发生的质变,即在分岔点会发生流动流形拓扑结构的改变.

对于浮力对流这一经典流体动力学耗散体系,对其流动发展至湍流的超临界转捩过程研究较多,也相对系统化.根据现有研究结果,向湍流转捩的典型分岔途径大致分为:

(1) 倍周期分岔途径(Feigenbaum 途径): 规则运动状态(如周期运动) 通过周期不断加倍方式转捩成为湍流(Feigenbaum 1979);

(2)准周期分岔途径: 运动通过由连续的周期运动变化至准周期运动的Hopf 分岔转捩为湍流,其间需要激发三个不可公约频率(RTN 途径),或者出现具有锁频特征的两个不可公约频率(Ruelle and Takens 1971);

(3) 阵发分岔途径: 运动出现周期运动与准周期运动、周期运动与不同混沌间的反复切换(间歇) (Pomeau and Manneville 1980);

(4) 切分岔途径: 在切分岔点附近运动发生奇数倍周期分岔,而后随着参数继续改变,可能由倍周期或准周期途径通向混沌(Broze & Hussain 1996),此外,切分岔也被认为与某类阵发现象的形成机制有关.

Gollub 和Benson (1980) 通过液层中Rayleigh-B´enard 流动的地面实验确认了若干向湍流转捩的途径 (图2): 随着 Rayleigh 数的增加,稳态流动演化成为周期性的振荡,然后通过倍周期分岔(Feigenbaum 途径)、具有锁频特征的准周期分岔、具有三个不可公约频率的准周期分岔(RTN 途径) 或者阵发途径中的一种转捩成为湍流,具体通过何种途径与液层的几何位形和流体物性有着复杂依赖关系,取决于相关耗散动力学体系控制参数,如Reynolds 数、Prandtl 数和体系几何位形(Mukutmoni & Yang 1993) 等.Braun 等(1998) 通过二维数值模拟研究了流动向湍流转捩的途径,分析了其与Reynolds 数、边界条件和长高比的关系.Saha 等(2000)识别了通过方柱的二维流动经由具有锁频特征的准周期途径向湍流转捩的过程.其他转捩途径,如倍周期分岔途径(Libchaber and Maurer 1978)、RTN 途径 (Ruelle & Takens 1971,Newhouse et al.1978)、阵发途径 (Pomeau &Manneville 1980) 也都在研究中得到确认.

Jiang 等(2017a,2017b) 通过矩形液层中热毛细对流的地面实验,在不同Prandtl 数和高径比下发现了几种流动转捩的切分岔途径,包括三周期分岔、五周期分岔及伴随准周期分岔的切分岔等,还观察到了阵发现象(图3),给出了相应的相轨迹描述(图4); 并利用红外热像仪实时测量了流体表面温度,发现了热毛细对流表面温度分布的波动特性(图5).近期,Kang 等(2019a) 在实践10 号返回式科学卫星上利用内加热环形液池中界面张力梯度驱动对流的空间实验研究了流动的超临界转捩现象,主要研究了体积比和不同加热模式对于稳定性和转捩途径的影响,验证了体积比效应理论; 发现流动失稳从驻波开始,通过竞争转化为行波,然后驻波行波耦合; 发现流动模式以波数4 转3 为主要形式,与体积比有关; 发现了二频准周期和1/2 倍周期这种单一的分岔过程,并与受浮力对流影响的地面实验结果进行了对比.由于空间实验条件所限,流动没有最终转捩进入湍流 (图6∼图9).Kang 等 (2019d,2020) 在 2016 年 —2019 年历时 34 个月,完成了 740 组液桥热毛细对流研究,积累了大量的空间实验数据.空间实验中已经发现了液桥高径比与体积比的综合效应新规律以及“跳变” 的新效应,发现低频模式起振区具有多次转捩; 给出液桥行波和驻波完整的转换图谱,探索了毛细流动体系通向混沌的分岔道路的复杂特征(图10∼图12).

图2

对于界面张力梯度驱动对流这一相对较新的流体动力学耗散体系,在对流动向湍流超临界转捩过程的探索中,无论是实验研究还是理论分析都相当缺乏.由前述可见,在绝大多数界面张力梯度驱动对流的地面实验研究中无法区分存在较强耦合作用的界面张力效应和重力效应;而空间实验的昂贵成本导致界面张力梯度驱动对流的空间实验研究屈指可数;因此,数值研究、理论分析和实验研究一起成为了相关内容科学探索研究的三个相互依存、不可缺少的方法.流动的控制方程组加上边界条件可以看作特定函数空间上的无穷维动力系统,系统的解是时空坐标下的函数,反映了流动状态.流动由层流向湍流转捩的现象可以用系统吸引子随参数改变而相应地发生变化来描述.动力系统以时间为参数的解曲线在相空间中的分布情况构成系统的相图,若在某个参数附近,相图的拓扑结构发生了改变,则称系统在这个参数处发生了分岔.分岔又可分为全局分岔和局部分岔: 较大的不变集与平衡点重叠,导致大范围内相图拓扑结构发生改变的现象,称为全局分岔,比如同宿分岔、异宿分岔等; 而平衡点、闭轨等解的稳定性或者数目发生变化的现象,一般称为局部分岔,研究局部分岔的方法一般是分析局部稳定性.本文重点关注局部分岔.在对流向湍流转捩的数值研究中,期望达到的目标是: (1) 得到给定参数下尽可能多的解,包括定常解、周期解,稳定、不稳定的解等; (2) 得到每个解随参数改变的变化情况,判断不同解之间随参数改变是否有联系,是如何联系的; 正确识别出分岔点,得到分岔的临界参数、分岔发生前后流动状态的变化; (3) 不同参数与其对应的解曲线的拓扑结构之间的关系构成分岔图,希望由上述(1) 和(2) 的内容总结出完整的分岔图,并在此基础上与不同系统的分岔图进行对比以探索更深层次的规律.

图3

图4

图5

图6

图7

需要指出的是,无论实验研究还是数值模拟,都需要最终对获得的速度场或温度场在时空坐标上的值开展非线性分析,用时间序列的频谱特性及相空间重构等方法来定量确定时序混沌.定常解失稳的临界参数及扰动增长后的波动特性可以通过分析速度场及温度场的时空演化图得到.但是随着参数的增大,解变得越来越复杂,流场的分布图不足以定量反应流动特性,也无法确认解是否为混沌解.出于上述考虑,需要进行简化操作,其中一个简化方向是将复杂的时空混沌简化为时序混沌,即固定空间上某个位置,取它的宏观量 (温度或者速度分量) 随时间的变化序列,考虑不同参数下对应的时间序列的特性,这样做实际上是默认了系统宏观量的时间序列可以反应出系统的动力学行为(虽然这样做会损失系统的部分信息).通过对流动的时间序列信息分析判断对流通向混沌的途径已经有了一套较为完整的流程: 首先,周期解及准周期解各自有其对应的频谱特性,进一步,由时频谱可以从整体上反映频率特性随参数改变而变化的过程,通向混沌的各种途径可由时频谱的具体特点来识别.Feigenbaum 途径的特点是原周期吸引子失稳,产生周期加倍的吸引子,重复此倍周期分岔过程,最终产生混沌吸引子,表现在频谱中就是依次出现1/2 基频,1/4基频,······,1/2n基频,······; Ruelle-Takens 途径的特点是经历了至少两次 Hopf 分岔,且第二次Hopf 分岔后通常出现混沌现象,表现在频谱中就是出现两个以上不可公约频率; 阵发途径的特点是原周期吸引子失稳,产生一个新的吸引子,这个新的吸引子时而是周期运动,时而是混乱的“随机”运动,且进入混乱状态的时间点是没有规律的.切分岔途径在频谱中的表现通常是出现1/3 基频或1/5 基频等.一般来说,当频谱呈现出以上几种途径的特性时,就可能预示着随参数的继续增长,最终系统会产生混沌吸引子.但是反之,混沌状态并不全是经由上述途径产生的,在一些研究中还观察到了其他特殊现象(Lappa 2009).其次,对于一个动力系统,一般在相空间中观察解的行为,但是直接数值模拟方法中的分析对象总是一维的时间序列,系统在原相空间中的行为隐藏在一维时间序列的信息中,为了还原这些信息,一般采用相空间重构法,其主要思想是由吸引子的一维时间序列重构出高维吸引子,高维吸引子位于一个重构的相空间中,Takens 定理(Takens 1981)从理论上保证了重构相空间中的吸引子与原系统的吸引子有相似的动力学特性(比如,吸引子的周期、准周期性质、吸引子的维数、是否为混沌吸引子等等),虽然重构出的空间并不是原系统真正的状态空间(因为真正的系统是无穷维的),但是根据以往的经验与研究结果,相空间重构法可以在一定程度上还原出系统主要的动力学特性.相空间重构的算法主要包括延迟时间的选取和嵌入维数的选取两方面,目前常用的方法有 G-P 算法(Grassberger & Procaccia 1983)、互信息法(Packard et al.1980,Fraser & Swinney 1986)、C-C 方法 (Kim et al.1999) 等,其主要思想大致是取使得重构序列在相空间中得到充分展开的延迟时间和嵌入维数为最佳的延迟时间和嵌入维数,以还原尽可能多的动力学信息.此外,在重构相空间中观察吸引子,不但可以进一步验证一维时间序列的周期及准周期特性,还为混沌时间序列的定量判断提供方便.例如,解对初值敏感一般量化为两个初始距离很近的解的距离随时间呈指数发散的程度,这一现象可以用Lyapunov 指数定量刻画,n维系统的Lyapunov 指数含n个指数,分别表示轨道在各个方向上的变化情况.Lyapunov 指数其实就是离散系统在多次迭代中平均每次迭代所引起的指数分离中的指数(吕金虎等2002),对于耗散系统,如果某个吸引子的最大Lyapunov 指数大于0,则可视其为混沌吸引子.在实际计算中,将Lyapunov 指数的前j个指数之和近似为前j个主轴定义的j维立体体积指数增加的长期平均速率,由此可计算所有的n个指数.此外,对耗散系统,更加关心最大Lyapunov 指数,计算时间序列的最大Lyapunov 指数的方法主要有: Wolf 方法(Wolf et al.1985)、Jacobian 方法(Barna &Tsuda 1993)、p- 范数算法 (Barna & Tsuda 1993)、小数据量方法 (Rosenstein et al.1993) 等,它们都是基于相空间重构方法对重构吸引子计算最大Lyapunov 指数,其中Wolf 方法较为常用,它适用于时间序列无噪声,切空间中小向量的演变高度非线性的情况.除此之外,还可以在重构相空间中运用关联维算法(Grassberger & Procaccia 1983)、排列熵算法(Bandt & Pompe 2002)等进一步确认混沌吸引子.

图8

图9

图10

图11

图12

基于以上非线性分析方法,Mukolobwiez 等(1998) 最早于1998 年实验研究了侧壁加热狭长环形管道内界面张力梯度驱动对流超临界热流体波的非线性动力学过程.研究结果表明该流动表现出超临界Eckhaus 不稳定性.后续Garnier 等(2003a,2003b) 研究了侧壁加热狭长环形和封闭管道内界面张力梯度驱动对流超临界热流体波和调制波的非线性动力学过程,并利用Ginzburg-Landau方程组分析了实验结果.Tang 等(1995) 通过地基小尺寸实验和数值模拟对液桥界面张力梯度驱动对流开展了研究,确认了倍周期分岔的转捩途径,并且对倍周期分岔的特性进行了相关分析,得到了与Feigenbaum 常数理论值非常接近的实验常数.Frank 和Schwabe (1997,1999)利用空间实验研究了不同流体介质液桥的界面张力梯度驱动对流的时空流动结构,并从空间结构的角度分析了准周期和倍周期这两种转捩途径.Bucchignani 和Mansutti (2004) 对Prandtl 数为23 的矩形液层中的浮力– 热毛细对流进行了三维数值模拟,发现当Marangoni 数固定,Rayleigh 数增加时,流动的转捩途径与Rayleigh 数的变化方式有关: 观察到了“周期解→二频准周期解→三频准周期解→混沌”这一途径; 并且通过取特定的初值,还得到了一个各阶段均为定常流的转捩过程(图13).Rahal 等 (2007) 以 Marangoni 数为控制参数,实验研究了 B´enard-Marangoni 对流向湍流转捩的过程,在一定参数范围确认了对流经历了周期运动和准周期运动最终发展成为湍流的过程.实验还观察到时间混沌向时空混沌的转捩过程.Chen (2010) 和 Li (2010) 等通过数值模拟方法研究了矩形液池内水平温度梯度和浓度梯度作用下的双扩散Marangoni 对流转捩过程,确认了准周期转捩 (RTN 途径) 途径.Aa 等 (2010) 通过数值模拟研究液桥界面张力梯度驱动对流,发现了倍周期分岔的转捩途径,并且对倍周期分岔各个阶段进行了相空间重构,得到了跟理论值非常接近的Feigenbaum 常数.Li 等(2012)对微重力环境下Marangoni 对流和热毛细对流耦合的液层进行了数值模拟,确定临界转捩的稳定边界,并报道了振荡耦合对流的多涡流结构转捩过程.Zhu 等(2013)对矩形液层中的浮力– 热毛细对流进行了地面实验,研究了Prandtl 数为10,16 和25 时不同长高比液层中对流随Marangoni 数增加发生的转捩过程,发现了“周期状态→准周期状态→混沌”途径及Feigenbaum 途径.Yu (2015) 等通过地面实验研究了环形浅液池界面张力梯度驱动对流热流体波特性,发现随着Re增加,对流依次经历二维稳态流、热流体波和混沌现象.作为国际空间站–中俄合作项目“矩形液池热毛细对流研究”(Zona-K)的预先研究,Li(2016)等通过数值模拟研究了二维液层热毛细对流分岔转捩过程,发现了流动向湍流转捩的具有锁频特征的准周期分岔途径与倍周期分岔途径,研究表明不同的液层长高比可能导致不同的转捩途径(图14).Zhang(2018)对小Prandtl 数为环形液层界面张力梯度驱动对流进行了三维数值模拟,并且考虑了自由面存在热耗散的情况.研究发现界面换热较小 (接近 0) 时,随着Marangoni 数增加,流动依次经历二维轴对称定常流→三维定常流→热流体波与径向行波共存→径向行波→混沌; 当界面换热较大(接近1) 时,随着Marangoni 数增加,对流依次经历了二维轴对称定常流→热流体波→混沌(图15 和图16).

图13

图14

图15

以上相关研究结果表明: 液层界面张力梯度驱动对流的各种工况下确实观察到了流动通过若干典型途径组合通向混沌,证实了这些理论上的转捩途径在实际物理场景中是可能的; 而另一方面,在一些实验中虽然观察到了具有以上若干分岔模式特征的频谱,但是流动最终没有发展进入混沌.这些例外背后可能是实验条件限制的原因,也可能蕴藏着别的机理解释.

2.3 分岔理论及其数值算法在流体力学研究中的应用

通过分析实验或直接数值模拟获得的流动时序数据研究转捩规律的整个过程比较繁琐,需要在大量离散的数据序列中寻找分岔点,而目前对于分岔点的确定还没有更简便的方法.近年来,一种新的思路,即分岔问题的数值算法得到了研究者的关注,这类方法基于分岔理论,构造系统的解及分岔点的代数方程,通过数值解法得到分支关于参数的曲线,在得到完整的分岔图上有其优势.

分岔的数值算法主要包括以下几个内容: (1) 计算系统定常解及其关于参数的曲线 (分支),(2) 确定分岔点,(3) 选取曲线通过分岔点后的分支方向

对于周期解,可用Poincare 截面将问题化为求定常解,所以只需考虑定常解关于参数的曲线,设系统为

图16

式中p为参数u ∈Rn.系统定常解满足方程

由隐函数定理及其相关推广,在系统满足一定光滑性的条件下,存在满足上述方程的状态空间中的整体解曲线u(p),求解这条曲线也就是要在不同参数下解定常解方程,所用的一类方法是延续算法 (continuation),其实本质上就是微分方程的迭代解法,最简单的就是牛顿迭代法,由于迭代初值会影响收敛结果,一般采用预估– 校正算法来得到更精确的迭代初值; 此外,在实际的场景中解的维数通常很高,这时可考虑用拟牛顿法来避免直接对高维矩阵求逆.在曲线的转折点(fold point or turning point) 处,不能再取参数p为曲线参数,此时一般采用拟弧长延续法(pseudoarclength continuation),引入参数s,将参数p也看作变量,建立关于s的迭代方程组.

为了用数值算法确定出分岔点,这里将分岔点定义为满足: 在曲线的分岔点的任意小邻域内,总存在不在曲线上的满足定常解方程的点.此外,由定理,使Jacobi 矩阵Φu(u(p),p) 的行列式符号改变的点就是一个分岔点.基于上述结论,发展出了许多判断分岔点的数值方法,最简单的想法是求解Jacobi 矩阵行列式的零点,可以直接求解,也可以通过计算矩阵特征值来求解,矩阵特征值的求解有一整套完整的理论及许多数值算法,比如Arnoldi 方法等.对于一些特定类型的分岔,也有特定的数值算法: 许多流体力学方程具有特殊的对称性,而破坏这种对称性的分岔自然就成为一个关注重点,Cliffe 等(2000) 总结了其中一种Z2对称破缺分岔的性质及其相应数值算法; Hopf分岔也是一类在理论与实际应用中比较重要的分岔,在Hopf 分岔点处,定常解失稳,产生周期解,体现在Jacobi 矩阵上的变化就是一对共轭特征值的实部由负变正,利用这一条件,Jepson (1981),Griewank 和Reddien (1983) 给出了计算Hopf 分岔的方程,但是转折点也满足他们给出的方程,随后,Werner 和Janovsky (1991)添加了额外的方程,解决了Jepson,Griewank 和Reddien 的方程无法排除转折点的问题.

事实上,直观地来看,分岔点处的有两个不同的非零切方向,沿这两个方向可继续延拓出通过分岔点之后的曲线.数值计算上,可以直接利用分岔点的性质解出两个切方向,作为延续算法下一步的初值,当已知其中一条分支时,还可利用Lyapunov-Schmidt 方法等计算另一条分支,此外,还有扰动方法(perturbed bifurcation) 等,上述方法的详情可见Keller (1987) 的书.

基于分岔理论建立分支曲线和分岔点的方程后,下一步要考虑求解所得到的高维线性方程.数值求解线性方程和特征值的理论庞杂,方法众多.对于实际的流体力学问题,离散后的方程一般维数较高,在运用迭代法解方程时可能会产生算法复杂度过高、无法求解,或收敛速度过慢等问题,需要根据流体系统的具体性质,选取合适的算法.例如,在计算定常解时,一种常采用的迭代算法是Krylov 子空间法,其中,为了降低子空间维数,可以用预处理矩阵或时间积分等方法简化方程(Tuckerman 2020),更为系统的处理方法可以参考Dijkstra 等(2014) 的综述文章.

上述分岔理论与数值算法在浮力对流和界面张力梯度对流的研究中已有相关应用: Salinger等 (2002,2005) 将延续算法推广到分岔点上,计算了二维方形腔内 Rayleigh-B´enard 对流的 Hopf分岔点关于长宽比的曲线,而后又计算了三维立方体腔内两种边界条件下的 Rayleigh-B´enard 对流的转折点、叉形分岔点、Hopf 分岔点,还得到了分岔点关于另一个参数Prandtl 数的变化曲线.Boro´nska 和Tuckerman(2010a,2010b)对圆柱形容器中Rayleigh-B´enard 对流进行直接数值模拟,在给定的Prandtl 数和体积比下、相同的Rayleigh 数处得到了几种不同的定常流和振荡流,确认了流动还与初始状态有关,又通过计算分岔点,得到定常解及周期解通过分岔相互转变的关系,给出了较为完整的分支图(图17).此外,文中对于三维流动的数值处理方法也有一定借鉴意义.Dijkstra和Henk (1992) 数值研究了小长高比的二维矩形腔中的Rayleigh-B´enard-Marangoni 胞状对流结构,将延续算法用于计算分岔点关于参数: 长高比、Biot 数 (Bi=ηH/κ,其中H为容器高度,η为表面传热系数)和Prandtl 数的变化曲线,发现更小的Biot 数或Prandtl 数以及滑移的边界条件会使定常流的斑图结构更加丰富.Chen 等(1997) 研究了自由面不变形的圆柱形液桥中的热毛细对流,通过解分岔方程计算了静态分岔和Hopf 分岔,发现在不同Prandtl 数和高径比下有不同的失稳模式.Cliffe 和Tavener (1998) 数值研究了二维矩形腔中考虑自由面变形的Marangoni-B´enard 对流,发现多种不同类型的分岔现象,其中计算对称破缺分岔点和转折点时用了Werner 和Spence(1980,1984) 文章中的方法.

图17

基于流动宏观量的时序数据的直接分析方法,比较繁琐,且数据具有离散性,数据量大,导致得到的分岔图具有一定的离散性和不完整; 分岔问题的数值算法是基于分岔理论,构造系统的解及分岔点的代数方程,通过数值解法得到分支曲线及分岔点,这种方法可以定位可能的分岔点的大致位置,避免了在大量时序数据中查找分岔点的盲目性,所计算出的分支包括稳定解和不稳定解,较之直接数值模拟所得到的稳定解更为全面,而且还可以得到不同分支之间的联系与转变关系的信息,是获得较为完整的流动形态及分岔图的有力工具.但是,不同的流体体系有各自的物理机制,方程及其对应的解可能有特殊的性质(比如对称性),在处理不同对流问题时,需要根据具体条件,采取不同的变换或是简化、近似的手段来得到足够准确的数值结果,总之,对于一个新的体系建立一套合适的分岔计算方法的过程并不比分析时序数据来得简单.由此可见,分岔的数值算法在界面张力梯度驱动对流上的应用仍具有相当的深入研究空间.

3 结论与展望

界面张力梯度驱动对流是空间自然对流热质输运的基本形式,对其时空转捩过程、转捩机制、非线性特征及流动向湍流转捩途径等基本规律的研究,一方面可以丰富非线性动力学的相关理论,另一方面对于人类认识、探索和利用空间环境也具有重要的应用价值,是微重力流体物理的重要研究内容和学科前沿.本文对目前的研究现状进行了总结,重点介绍了研究液层界面张力梯度驱动对流的实验及数值模拟方法,虽然已有的研究已经得到在不同模型和工况下的各种转捩模式,但是在转捩规律上仍需要更深入的探索,可以从以下几个方面着手:

(1)理论分析和数值模拟结果的正确性需要由实验来验证,空间实验可以满足微重力环境、长时间观测的要求,但是空间实验有一定难度且机会来之不易,故而可以考虑进一步发展实验手段,以实现数值模拟中采用的更丰富的工况; 以及加强对实验条件的控制,以降低无关因素的干扰,提高实验精度.

(2) 目前关于液层界面张力梯度驱动对流向湍流的超临界转捩在数值方法上主要有流动时序数据的分析和分岔问题的数值算法.流动宏观量的时序数据来自实验结果或者直接数值模拟,对于后者,需要对于不同参数分别进行数值模拟,再通过时间序列频谱及其混沌特性的定量计算分析流动转捩规律,即在大量的离散的数据序列中寻找分岔点,故此过程比较繁琐.而通过构造分岔方程对分岔进行数值计算的方法虽然可以一步到位,但是在选取分岔方程,解高维线性、非线性方程等过程中均需要根据具体的流动模型进行调整,具有一定难度,且对于更加复杂的流动模式需要更大的计算量,用此算法也无法直接计算得到混沌解.上述两种方法各有优缺点,目前在转捩过程的数值研究中较为常用的仍是在不同参数下进行直接数值模拟,而后对大量数据进行频谱分析,识别分岔点; 在直接对分岔进行数值计算的研究中,也常常需要通过直接数值模拟来验证分岔得到的解的可靠性与准确性,在今后的研究中可考虑进一步将两种方法结合运用,互相补充、验证.

(3) 液层界面张力梯度驱动对流向湍流转捩的过程中会产生丰富的流动模式,转捩过程除了与上文提到的液层模型、无量纲参数(Prandtl 数、高径比、体积比等) 有关,还受到热边界条件(如体系是否绝热)、加热方式及加热速率等因素的影响; 此外,在具体的应用场景中通常有多种流动相互作用,考虑界面张力梯度驱动对流与其他诸如浮力、电磁场、旋转等效应的耦合,对于重新检视已发现的转捩途径以及寻找新的转捩途径均有一定的积极意义.

(4) 目前对于液层界面张力梯度驱动对流向湍流转捩的研究仍不够完善,在对超临界转捩阶段的实验及数值模拟研究中观察到了许多复杂的转捩模式,但大多只是现象上的描述,并未总结出普遍的规律; 对于流动最终能否通向混沌暂无普适的判据,流动通向混沌过程中出现的诸如阵发、锁频等特殊的现象也尚未有更本质的机理上的解释.总之,对于转捩规律的深入理解,需要界面张力梯度驱动对流这一非线性模型在理论上的进一步发展,未来道阻且长.

致谢国家自然科学基金资助项目(11972353,U1738116).