旋转的Rayleigh-Bénard问题Lorenz模型的动力学行为及数值仿真

王贺元, 陈相霆

(沈阳师范大学 数学与系统科学学院, 沈阳 110034)

1963年,Lorenz[1]将截谱方法作用于Rayleigh-Bénard,在流模型上推导出了三维Lorenz方程。后来,国内外的众多学者研究了Lorenz模型的高维数推广、Lorenz吸引子的存在性及性态[2-3]、Couette-Taylor流的力学机理等[4-10]。Bhattacharjee等[11]和张银[12]先后对旋转Rayleigh-Bénard问题的四维Lorenz模型进行了讨论。本文对他们的模型进行了进一步讨论,进行了全局稳定性分析并对系统的动力学行为进行了数值仿真。

1 数学模型

旋转的Rayleigh-Bénard问题无量纲化后的扰动方程可通过如下的偏微分方程来描述[13]:

(1)

边界条件为自由边界∂zu|z=0,1=∂zv|z=0,1=w|z=0,1=θ|z=0,1=0。式(1)中:u=(u,v,w)表示速度的扰动;θ表示温度的扰动;p表示压强;R是Rayleigh数;T是Taylor数,物理意义是旋转的大小;k=(0,0,1);Pr是Prandtl数。

对方程组(1)化简后展开,选取模式后代入式(1)得到方程组(2)[13]:

(2)

方程组(2)中的X与对流的强度成比例,Y与上下层流体间的温度差成比例,Z与垂直温差的非线性强度成比例,G与流体涡旋强度成比例。

2 全局稳定性分析

2.1 全局吸引子的存在性

对于方程(2)作代换xx,yy,zz+r+Pr后作运算取u(t)=(X(t),G(t),Y(t),Z(t)),H=R4,由上述变量表示可知,参数r,Pr,b均为正数。

2.2 全局稳定性分析

下面利用李雅普诺夫函数法对系统(2)进行全局稳定性分析。构造李雅普诺夫函数

V(x,g,y,z)=13x2+13g2+5y2+5(z-63)2=K>0

很明显,当K是常数时,上式表示一个四维椭球面,把这个椭球面所包围的单连通区域记做E,K越大,椭球面越大。在方程组(2)的混沌区(Pr=5,b=8/3,r=50,τ=1)求V的导数:

显然,下式可以视为一个四维椭球面,记为U:

由李雅普诺夫定理[14]的分析可知,E外的轨线都会进入E内。由此可知,E是这个Rayleigh-Bénard系统的捕捉区。虽然该系统的平衡点都不稳定,但是系统依然具有全局稳定性。

3 数值仿真

改变参量r的取值,借由数值仿真给出的7种指标分析系统的混沌现象。取Pr=5,b=8/3,τ=1,在参数r变化时,对系统的动力学行为进行仿真。

图1给出了系统(2)随参数r变化的最大李雅普诺夫指数图像,可以看出系统的分岔过程对应李雅普诺夫指数的变化。图2给出了系统(2)随r变化的分岔图,它展示了系统分岔和混沌演变的全过程。系统发生混沌后,出现了明显的倒分岔现象。

图1 最大李雅普诺夫指数图像Fig.1 Maximum Lyapunov exponent image

图2 分岔图Fig.2 Bifurcationdiagram

如图3和图4所示,当 0≤r≤1 000时,平衡点逐渐不稳定,从全局吸引子开始生成2个不稳定的极限环,并且轨线条数随r的增大而逐渐增多,最终产生了奇怪吸引子。之后吸引子在奇怪吸引子、拟周期吸引子和极限环之间变换,最终变为拟周期吸引子。

图3 r=20时的吸引子图像Fig.3 Attractor image at r=20

图4 r=900时的吸引子图像Fig.4 Attractor image at r=900

以r=50时为例,给出图5～8的混沌指标。参考最大李雅普诺夫指数(图1)和分岔(图2),可知此时系统(2)处于混沌状态。

图5 r=50时的庞加莱截面Fig.5 Poincare section at r=50

图6 r=50时的时间序列Fig.6 Time series at r=50

图7 r=50时的功率谱Fig.7 Power spectrum at r=50

图8 r=50时的返回映射Fig.8 Return mapping when r=50

4 结 语

通过对系统的全局稳定性分析可知,旋转的Rayleigh-Bénard问题的四维Lorenz模型是全局稳定的;通过对系统模型(2)的数值仿真可知,Rayleigh-Bénard系统的动力学行为随着参量r取值而变化。在[0,280],奇怪吸引子、拟周期轨道和极限环并存,系统存在倒分岔过程;在[280,500]吸引子稳定为极限环的形式。这些现象说明系统的稳定性是增加的,与文献[13]中的结论一致。