旋轉(zhuǎn)的Rayleigh-Bénard問題Lorenz模型的動(dòng)力學(xué)行為及數(shù)值仿真

王賀元, 陳相霆

(沈陽師范大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院, 沈陽 110034)

1963年,Lorenz[1]將截譜方法作用于Rayleigh-Bénard,在流模型上推導(dǎo)出了三維Lorenz方程。后來,國(guó)內(nèi)外的眾多學(xué)者研究了Lorenz模型的高維數(shù)推廣、Lorenz吸引子的存在性及性態(tài)[2-3]、Couette-Taylor流的力學(xué)機(jī)理等[4-10]。Bhattacharjee等[11]和張銀[12]先后對(duì)旋轉(zhuǎn)Rayleigh-Bénard問題的四維Lorenz模型進(jìn)行了討論。本文對(duì)他們的模型進(jìn)行了進(jìn)一步討論,進(jìn)行了全局穩(wěn)定性分析并對(duì)系統(tǒng)的動(dòng)力學(xué)行為進(jìn)行了數(shù)值仿真。

1 數(shù)學(xué)模型

旋轉(zhuǎn)的Rayleigh-Bénard問題無量綱化后的擾動(dòng)方程可通過如下的偏微分方程來描述[13]:

(1)

邊界條件為自由邊界?zu|z=0,1=?zv|z=0,1=w|z=0,1=θ|z=0,1=0。式(1)中:u=(u,v,w)表示速度的擾動(dòng);θ表示溫度的擾動(dòng);p表示壓強(qiáng);R是Rayleigh數(shù);T是Taylor數(shù),物理意義是旋轉(zhuǎn)的大小;k=(0,0,1);Pr是Prandtl數(shù)。

對(duì)方程組(1)化簡(jiǎn)后展開,選取模式后代入式(1)得到方程組(2)[13]:

(2)

方程組(2)中的X與對(duì)流的強(qiáng)度成比例,Y與上下層流體間的溫度差成比例,Z與垂直溫差的非線性強(qiáng)度成比例,G與流體渦旋強(qiáng)度成比例。

2 全局穩(wěn)定性分析

2.1 全局吸引子的存在性

對(duì)于方程(2)作代換xx,yy,zz+r+Pr后作運(yùn)算取u(t)=(X(t),G(t),Y(t),Z(t)),H=R4,由上述變量表示可知,參數(shù)r,Pr,b均為正數(shù)。

2.2 全局穩(wěn)定性分析

下面利用李雅普諾夫函數(shù)法對(duì)系統(tǒng)(2)進(jìn)行全局穩(wěn)定性分析。構(gòu)造李雅普諾夫函數(shù)

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

很明顯,當(dāng)K是常數(shù)時(shí),上式表示一個(gè)四維橢球面,把這個(gè)橢球面所包圍的單連通區(qū)域記做E,K越大,橢球面越大。在方程組(2)的混沌區(qū)(Pr=5,b=8/3,r=50,τ=1)求V的導(dǎo)數(shù):

顯然,下式可以視為一個(gè)四維橢球面,記為U:

由李雅普諾夫定理[14]的分析可知,E外的軌線都會(huì)進(jìn)入E內(nèi)。由此可知,E是這個(gè)Rayleigh-Bénard系統(tǒng)的捕捉區(qū)。雖然該系統(tǒng)的平衡點(diǎn)都不穩(wěn)定,但是系統(tǒng)依然具有全局穩(wěn)定性。

3 數(shù)值仿真

改變參量r的取值,借由數(shù)值仿真給出的7種指標(biāo)分析系統(tǒng)的混沌現(xiàn)象。取Pr=5,b=8/3,τ=1,在參數(shù)r變化時(shí),對(duì)系統(tǒng)的動(dòng)力學(xué)行為進(jìn)行仿真。

圖1給出了系統(tǒng)(2)隨參數(shù)r變化的最大李雅普諾夫指數(shù)圖像,可以看出系統(tǒng)的分岔過程對(duì)應(yīng)李雅普諾夫指數(shù)的變化。圖2給出了系統(tǒng)(2)隨r變化的分岔圖,它展示了系統(tǒng)分岔和混沌演變的全過程。系統(tǒng)發(fā)生混沌后,出現(xiàn)了明顯的倒分岔現(xiàn)象。

圖1 最大李雅普諾夫指數(shù)圖像Fig.1 Maximum Lyapunov exponent image

圖2 分岔圖Fig.2 Bifurcationdiagram

如圖3和圖4所示,當(dāng) 0≤r≤1 000時(shí),平衡點(diǎn)逐漸不穩(wěn)定,從全局吸引子開始生成2個(gè)不穩(wěn)定的極限環(huán),并且軌線條數(shù)隨r的增大而逐漸增多,最終產(chǎn)生了奇怪吸引子。之后吸引子在奇怪吸引子、擬周期吸引子和極限環(huán)之間變換,最終變?yōu)閿M周期吸引子。

圖3 r=20時(shí)的吸引子圖像Fig.3 Attractor image at r=20

圖4 r=900時(shí)的吸引子圖像Fig.4 Attractor image at r=900

以r=50時(shí)為例,給出圖5～8的混沌指標(biāo)。參考最大李雅普諾夫指數(shù)(圖1)和分岔(圖2),可知此時(shí)系統(tǒng)(2)處于混沌狀態(tài)。

圖5 r=50時(shí)的龐加萊截面Fig.5 Poincare section at r=50

圖6 r=50時(shí)的時(shí)間序列Fig.6 Time series at r=50

圖7 r=50時(shí)的功率譜Fig.7 Power spectrum at r=50

圖8 r=50時(shí)的返回映射Fig.8 Return mapping when r=50

4 結(jié) 語

通過對(duì)系統(tǒng)的全局穩(wěn)定性分析可知,旋轉(zhuǎn)的Rayleigh-Bénard問題的四維Lorenz模型是全局穩(wěn)定的;通過對(duì)系統(tǒng)模型(2)的數(shù)值仿真可知,Rayleigh-Bénard系統(tǒng)的動(dòng)力學(xué)行為隨著參量r取值而變化。在[0,280],奇怪吸引子、擬周期軌道和極限環(huán)并存,系統(tǒng)存在倒分岔過程;在[280,500]吸引子穩(wěn)定為極限環(huán)的形式。這些現(xiàn)象說明系統(tǒng)的穩(wěn)定性是增加的,與文獻(xiàn)[13]中的結(jié)論一致。