水轮混沌旋转的力学机理与能量演化研究*

王贺元, 肖胜中, 梅鹏飞, 张 熙

(1.广东科技学院 通识教育学院, 广东 东莞 523083;2.沈阳师范大学 数学与系统科学学院, 沈阳 110034;3.广东农工商职业技术学院, 广州 510507)

0 引 言

混沌研究的历史最早可追溯到Poincare对三体问题的研究.1963年美国的气象学家Lorenz在研究局部区域小气候的数值实验时发现了混沌现象,开启了混沌研究的先河,引发了后续的大量研究,相关文献非常丰富[1-12].受数值结果的启发,Lorenz设计了与我国古代的水车相似的混沌水轮实验装置,其转动的示意图如图1所示.其主要组成部分为存在摩擦阻力的水平轴和竖直的与其相连的轮,水轮顶端连有水管,可以将速度可调节的水流注入到挂在轮边缘的水杯中,每只水杯的底部均有一个能恒定漏水的小孔.如果注入水流的速度很慢,顶部杯中水量较少,因而不能克服轮轴摩擦力,水轮静止不动.如果水流加大,杯中的水来不及漏出,随着顶部杯中水的增多,产生力矩驱动水轮转动.如果保持水流速度恒定,水轮开始匀速旋转,随着水流继续加大,水会不断地注入其他水杯中,这些水杯来不及装满就转到下侧,下侧的水杯来不及漏空就转到顶端,形成了反向力矩,从而导致转速减小,甚至发生逆转.旋转逐渐变得无序混乱,发生非匀速旋转、倒转、周期反转等现象,最终呈混沌状态,转动的方向和速度会因系统内在的非线性而显现出非常复杂的运动特征,可以观察到十分丰富的混沌旋转现象.1972年,Malkus改进了Lorenz混沌水轮实验装置,完成了通过实验演示Lorenz方程混沌行为的这一挑战性工作,引起了许多学者的广泛关注,开展了分析和解释水轮混沌旋转现象的一系列研究工作,具体见文献[6,8,10].上述研究工作都是分析和解释实验中观察到的水轮的旋转状态及其相互演化的过程等,而对水轮为什么会混沌旋转、混沌旋转的生成机理、其能量转换等问题目前还没有文献涉及,本文将在这方面展开研究与探索.

图1 混沌水轮示意图

探讨混沌生成的力学机理及其能量演化方面的研究已取得了一些进展.文献[13]利用Kolmogorov系统描述了具有Hamilton函数的不同强迫动力系统、流体动力系统等.文献[14]对Lorenz系统进行了研究,给出了统一的Kolmogorov和Lorenz系统.文献[15]利用Casimir函数分析了Chen系统的动力学行为和能量转换,揭示了能量转换和轨道与平衡点间的距离,估计了混沌吸引子的边界.文献[16]分析了Chen系统的动力学机理与能量转换.文献[17]研究了Qi四翼混沌系统的力学机理与能量转换,通过力矩的5种耦合模式分析了四翼混沌系统产生混沌的动力学机制.借助扩展的Kolmogorov系统,文献[18]研究了Lorenz系统的能量演化.文献[15-17]中的混沌系统均为通过数值计算发现的混沌系统,没有明确的物理意义.文献[11-12]探讨了同轴圆筒间旋转流动的动力学机理,研究了具有明确物理意义的C-T流问题,分析了C-T流的动力机制、物理意义和能量转换.

基于这些研究工作,本文探讨Malkus水轮混沌旋转的内在机理.将Malkus水轮系统转化为Kolmogorov系统, 分析和仿真了不同动力学模式下Malkus水轮模型的旋转行为, 进而解释了Malkus水轮混沌旋转的力学机理.

因为Malkus水轮系统的轨道和平衡点之间的距离随参数和时间变化,所以理论分析和数值模拟Malkus水轮系统的混沌吸引子是非常困难的.由于Casimir函数与距离息息相关,本文借鉴文献[18]关于Lorenz吸引子的研究思想,引入Casimir函数研究了Malkus水轮系统全局吸引子的边界和能量演化问题.

1 Malkus水轮系统及其Kolmogorov系统

文献[6]利用质量守恒方程、动量矩方程结合Fourier展开等数学方法推导出 Malkus水轮的数学模型为如下三维非线性微分方程组(以下简称为Malkus水轮系统):

(1)

这里x,y,z是Fourier系数,它们是时间t的函数,x是与水轮旋转角速度有关的物理量,y和z是与水杯中水的质心坐标有关的物理量,σ是正的参数,r为Rayleigh数,它们是与漏水率和注水率有关的变量.由于系统(1)是著名的Lorenz系统的特例(b=1),其混沌现象是显著的.经数值仿真得当σ=5,5.45≤r≤210.51,212.07≤r≤301.53,348.13≤r≤510.53时,系统(1)存在混沌吸引子.其分岔图和最大Lyapunov指数如图2所示.

(a)分岔图 (b)最大Lyapunov指数

为探讨水轮混沌旋转的力学机理,我们引进如下的Kolmogorov系统:

(2)

其中X=[x,y,z]T,反对称括号{·,·}表示Hamilton函数H动能部分的代数结构,与cosymplectic矩阵J,或Lie-Poisson结构[19],

{F,G}=Jik∂iF∂kG.

(3)

(4)

系统(4)与原系统(1)平衡点的数量和性质相同,但位置发生了改变.系统(4)的平衡点如下:

下面采用线性稳定性分析方法讨论平衡点的稳定性.

容易求得方程(4)线性化方程在平衡点O处的特征方程为

(λ+1)[λ2+(σ+1)λ+σ(1-r)]=0.

(5)

特征值为

(6)

相应的特征矢分别为

(7)

先看r=r1=1的情形,当rr1时,O变为不稳定的,系统分叉出另外两个稳定的定态P+和P-.故在r=r1=1时出现音叉分岔(图3(a)),具体分析请参考文献[2-5].

(a)r=r1处的音叉分岔 (b)P±点r=rh处的亚临界Hopf分岔

下面讨论定点P+和P-的稳定性.

系统(4)具有反射对称性:当(x,y,z)变为(-x,-y,z)时,方程不变.因此若(x0,y0,z0)是方程的解,则(-x0,-y0,z0)自然也是方程的解,而且此二解具有相同的性质,即P+(x0,y0,z0)和P-(-x0,-y0,z0)性质完全相同,故只需分析其中之一即可,下面只分析P+.

系统(4)的线性化方程在P+处的特征方程为

λ3+(σ+2)λ2+(σ+r)λ+2σ(r-1)=0.

(8)

于是得到其Routh-Hurwitz判别行列式为

(9)

令Δ2=(σ+2)(σ+r)-2σ(r-1)=0时的r为rh,则

(10)

当σ<2时,rh<0,但r(表示无量纲Reynolds数)取负值是没有意义的,因此这种情形应摒弃,而只限于讨论σ>2的情形.由式(10)可知:当r0,Δ3>0;反之,当r>rh时,Δ2<0,Δ3<0.因此根据Routh-Hurwitz判据可知:当rrh时,P+和P-都是不稳定的.

当r=rh时,特征方程化为

(λ+σ+2)[λ2+(rh+σ)]=0.

(11)

定义Hamilton能量H=K+U,其中动能K=(x2+2y2+2z2)/2,势能U=σz.则系统(4)可描述为如下的Kolmogorov系统:

(12)

其中

Λ=diag(σ,1,1),f=(0,0,-r(1+σ))T.

2 系统(4)的动力学机理分析

各种类型的力矩对Malkus水轮旋转都具有不同程度的影响,本节从力矩的不同耦合模式展开讨论,从而阐释了水轮混沌旋转的力学机理.

情形1 系统只包含由动能K产生的惯性力矩:

(13)

显然系统(13)是保守系统,其闭合的周期轨道如图4(a)所示,状态变量z的轨迹如图4(b)所示.由于x是常量,系统(13)是线性系统.

(a)周期轨道的3维视图 (b)状态变量z的轨迹

情形2 系统包含惯性力矩和内力矩,即系统仅包含内能,相应的方程为

(14)

由于没有耗散,系统(14)也是保守系统,其闭合的周期轨道如图5(a)所示,状态z的轨迹如图5(b)所示.比较情形1和情形2,势能U的出现使解更频繁地振荡.由于σ>0,系统是非线性系统,并且势能U的存在使产生的周期解频率比情形1大得多,这意味着内力矩使系统解移动的速度比情形1快得多.频率取决于参数σ,由U释放的内力矩导致了轨道的拉伸和收缩,这对系统不稳定和混沌的生成是潜在有效的.

(a)周期轨道的3维视图 (b)状态变量z的轨迹

情形3 系统包含惯性力矩、内力矩和耗散力矩,即系统包含内能和耗散因素,但不包含驱动因素,此时方程为

(15)

对于σ>0,容易获得

其中V是系统相空间的体积,此时系统是耗散的[20].Hamilton函数的导数为

与文献[14, 18]不同,无法通过Hamilton能量的变化来确定能量耗散.图6(a)给出了Hamilton函数H的时间演化,显示能量是耗散的,由于V收缩从而能量减少,如图6(b)所示.

图6 能量函数的时间演化

情形4 系统在惯性力矩、内力矩和外力矩作用下,此时系统包含内能和驱动因素,但不包含耗散因素,方程为

(16)

(a)螺旋状轨道的3维视图 (b)状态变量z的轨迹

图8 能量的演化

为了分析比较各种力矩对Malkus水轮转动的影响,上面构造的4种力矩缺失模式均为虚拟的状态,下面的全力矩模式5才具有真实的物理意义.

情形5 系统包括全部力矩,此时系统包含内能、耗散因素和驱动因素,方程为

(17)

当r=32时系统(17)的混沌吸引子如图9(a)所示.情形3缺少外力矩,Malkus水轮系统的解趋于平衡点O.情形4的Malkus水轮系统的解由于没有耗散而无限增长.因此,外力和耗散耦合是Malkus水轮系统产生混沌的必要条件.当外力矩与耗散力矩不匹配时,耗散并不能保证Malkus水轮系统的能量衰减.虽然外力和耗散是产生混沌的基本因素,但外力和耗散简单的耦合并不足以保证系统产生混沌,例如,大参数r可以使系统具有周期性或无穷大,如图9(b)所示.图10(a)绘制了当r=376时,系统的准周期吸引子.只有当参数r在一定范围内取值(σ=5,5.45≤r≤210.51,212.07≤r≤301.53,348.13≤r≤510.53),也就是外力与耗散相匹配时,系统才会产生混沌.图10(b)绘制了动能和势能的时间演化(r=30),其中上曲线表示动能,下曲线为势能,当动能接近波峰时,势能到达波谷,说明两种能量相互转化.

(a)r=32 (b)r=34

(a)拟周期轨道的3维视图 (b)动能和势能

3 系统(4)的全局稳定性分析

混沌系统解的边界性质在混沌系统研究中是至关重要的,确定混沌吸引子的边界通常是困难的.本文利用Lagrange乘数法和Casimir函数法分析了Malkus水轮系统混沌吸引子的边界,数值模拟给出了清晰的边界.Casimir函数C由括号(3)的内核定义,即

{C,G}=0, ∀G∈C∞(g*),

这意味着在Lie-Poisson括号下Casimir函数与每个函数交换[19].对于Malkus水轮系统来说,Casimir函数定义为

由方程(3)和文献[19],我们有

{C,G}=-X·(∇C×∇G)=-X·(X×∇G)=0,

其中X=[x,y,z]T.根据方程(4)得

令

定理1 Casimir函数被限制在如下集合内:

Ξ={(x,y,z)|x2+y2+z2≤(r+σ)2}.

证明根据上述分析,C(t)的最大最小值点均位于集合Ξ0内,因此,可由如下条件极值问题来获得C(t)的上界:

(18)

定义Lagrange函数:

其中λ是Lagrange乘子.

经计算得到如下偏导数:

(19)

得到两个解

(x,y,z,λ)=(0,0,0,0) or (0,0,-(r+σ),-1).

所以,C的最大值为(r+σ)2/2,即

定理1给出了Malkus水轮系统混沌吸引子边界的精确估计.Casimir函数的值如图11(a)(r=40)所示,可见函数C(t)以C(t)<(r+σ)2/2为上界振荡.混沌吸引子被包围在矩形边界内,如图11(b)所示.为了显示Casimir函数与平衡点P±之间的关系,定义两种距离:

(a)Casimir能量时间演化 (b)混沌吸引子的边界

D1(t)=|X(t)-P+|,D2(t)=|X(t)-P-|,

(20)

这表示系统(12)的轨道与平衡点P±之间的距离.在图12(a)(r=30)中,实线是Casimir函数,虚线是D1(t)和D2(t).图12(a)表明,当轨线远离平衡点P+或P-时,Casimir函数达到最大值.当轨线同时接近平衡点P±时,Casimir函数达到最小值.图12(b)(r=31)显示了更紧密的关系,在这个图中,虚线是D1(t)与D2(t)之和,它与Casimir函数具有类似的动态.其极值点几乎在同一时刻,而且有相同的上升和下降趋势.

图12 函数与距离D1,D2的关系

当注水率(参数r)增大时,系统(4)的动能增加,如图13所示,外力矩主要增加了动能,导致水轮失稳而出现混沌旋转.图14显示了距离D1和D2的总和与Casimir函数随参数r的变化趋势.

(a)能量 (b)Casimir函数

(a)D1与D2的距离和 (b)Casimir函数

根据以上理论分析和仿真结果,得出如下结论,当r<1时,系统(4)仅存在平衡点O,这种定常的状态表明水轮处于静止状态.当r>1时,系统(4)存在以下3个平衡:O,P+,P-,平衡点O从稳定结点变成鞍结点,新平衡点P±是稳定的,这种稳定的状态代表水轮匀速的旋转,如图1(a)所示.新的平衡点随着r增加逐渐丧失其稳定性,它们从稳定结点发展到稳定焦点,最后变成了鞍结点.同时,系统(4)的轨线趋于P±,在P+和P-之间来回跳跃,这对应水轮的实际状态是出现倒转现象,依次正向或反向旋转如图1(b)、1(c)所示.距离D1与D2之和随r增大而单调递增.动能从最小值逐渐增加,Casimir函数也有类似的增加趋势.随着水轮稳定性的丧失(图1(a))正向或反向旋转相继出现(图1(b)、1(c)),在耗散和驱动的内禀作用下,系统的内能不断增大,从分岔经由暂态混沌最终到达混沌,对应水轮从正向或反向的匀速旋转,然后到不稳定的非匀速旋转,进而到达混沌旋转.表1给出了相关结论的细节.

表1 σ=5,r取不同值时,水轮系统(1)的动力学行为与能量演化及其相应的旋转状态

4 能 量 转 换

系统(12)的首项是由Hamilton能量传递的力矩,其似乎与耗散和外力矩无关,但从图13中发现能量与这两项有关.由于每种力矩是耦合的线性或非线性的矢量,其作用于质点上,因此很难研究系统(12)的力矩.由于能量是一个标量,所以通过研究能量来发现力矩特征是很容易的.下面讨论能量与这两种力矩的关系.文献[22]通过Lyapunov函数将Hamilton动力学的Lie-Poisson括号扩展为

(21)

(22)

其中

方程(22)解释为整个混沌系统包含动能、势能、耗散能和外力能,这些能量转换为惯性力矩、内力矩、耗散力矩和外力矩.对于水轮系统,包含注水能量和源自于水重力、摩擦损失以及漏水的耗散能量,作用于水轮上的总力矩产生角加速度.

根据方程(3)和(22),有

(23)

那么

{K,H}+〈K,L〉+〈K,G〉={K,U}+〈K,L〉+〈K,G〉-

σxy-[σx2+2y2+2z2+2(r+σ)z],

(24)

其中{K,U}=-σxy.因此,动能的变化率与势能、耗散能、外能量有一定的关系.同理,也有

(25)

这意味着势能的变化率等于它与动能、耗散能量和外部能量的交换.

-(2L+G)[21].Casimir函数是热力学能,则Casimir函数变化率为耗散与供给能量之间的交换率.类似的有

-[2σ2-2σ+(r+σ)]xy+2σ2x2+2y2+2z2+3(r+σ)z+(r+σ)2.

(26)

然后,水轮系统的能量转换为

(27)

这考虑了在K,U和W交换项中的耗散和力.方程(24)—(27)表明参数r对水轮系统的能量演化有显著影响.

在上面的公式中,当xy<0时,{K,U}是正数,因而有势能转变为动能的净转换.相应地,当xy>0时,{K,U}为负数,有动能变为势能的净转换.

Malkus水轮系统有螺旋式轨道,它的轨道从一个不稳定的平衡点移动到另一个不稳定的平衡点.当轨道远离不稳定的平衡点之一时,D1与D2的距离和因参数r增大而逐渐增大,如图14所示.也就是r增大,使得驱动与耗散相匹配,二者共同作用导致系统内能增大,水轮不稳定最终出现混沌旋转.因此,轨线在两个不稳定的平衡点之间频繁跳跃导致了D1与D2的距离和逐渐增大.

Malkus水轮混沌系统作为Lorenz系统的特殊形式,其能量转换时间演化与文献[18]的结果相似,能量、Casimir函数及D1与D2的距离和与外力矩(Rayleigh数r)间的演化关系是本文的重要发现.

5 结 论

本文研究了Malkus水轮混沌系统的力学机理和能量演化问题,揭示了水轮混沌旋转的力学机制.首先,探讨了Malkus水轮的数学模型作为Kolmogorov系统的物理意义.研究了水轮混沌系统4种力矩耦合的5种情形,分析了导致水轮不稳定和混沌旋转的关键因素.保守情形下,Hamilton能量为常量,对应的方程有周期解.外力矩或耗散力矩加入到保守系统中,Hamilton能量趋于无穷或零,对应的系统不可能产生混沌.只有全力矩模式下,水轮系统才产生混沌,并且动能和势能相互转化.对于水轮混沌系统,内能、耗散和驱动因素并存是产生混沌的必要条件,而且,只有当耗散与驱动相匹配时,系统才能生成混沌.注水为水轮系统提供能量,增加的动能导致水轮不稳定,依次发生正向稳定旋转、反向稳定旋转等经过暂态混沌旋转,最终到达混沌旋转.其次,将水轮混沌系统作为扩展的Kolmogorov系统,引进Casimir函数来分析系统动力学行为和能量转换.经分析获悉:平衡点的距离和与Casimir函数的时间演化是一致的.作为内能的Casimir函数的变化率是耗散与供给能量间的交换率,起着能量转换的作用,平衡点距离的时间演化也是类似的.利用Lagrange乘数和Casimir函数,分析仿真获得了混沌吸引子的边界.