数学物理学报  2017, Vol. 37 Issue (1): 199-216   PDF    
扩展功能
加入收藏夹
复制引文信息
加入引用管理器
Email Alert
RSS
本文作者相关文章
王贺元
平面不可压缩磁流体动力学五模类Lorenz方程组的动力学行为及其数值仿真
王贺元     
辽宁工业大学理学院 辽宁锦州 121001
摘要:该文研究了平面正方形区域上不可压缩的磁流体动力学方程组五模截断所得到的十维模型的动力学行为问题.首先,利用模式截断方法推导了十模系统,讨论了该方程组定常解及其稳定性,其次,发现了Hopf分叉和混沌,证明了该方程组吸引子的存在性和全局稳定性,最后,给出了系统从分叉到混沌整个过程所呈现的动力学行为演变的详细数值模拟结果,分析了磁性对系统动力学行为的影响.基于分岔图、Lyapunov指数谱和庞加莱截面图,返回映射和功率谱等数值模拟结果揭示了这个低维系统的动力学行为特征.这个新混沌系统通过周期倍分岔过渡到混沌(费根鲍姆途径).
关键词Hopf分歧     李雅普诺夫函数     奇怪吸引子     混沌    
The Dynamical Behaviors and the Numerical Simulation of a Five-Mode Lorenz-Like System of the MHD Equations for a Two-Dimensional Incompressible Fluid on a Torus
Wang Heyuan     
College of Sciences, Liaoning University of Technology, Liaoning Jinzhou 121001
Abstract: In this paper, we investigate a ten-dimensional model of plane flow, obtained by a five-mode truncation of the magnetic hydrodynamic (MHD) equations for a two-dimensional incompressible fluid on a torus. First, we derive a model system by taking five modes in the Fouriers expansion, and then discuss the stability of stationary solutions of the five-mode system. Second, we find phenomena of Hopf bifurcation and chaos, and prove the existence of its attractor and global stability of the five-mode system. Finally, we present a detailed numerical result of the whole process from bifurcation to chaos, and analyze the influence of magnetism on the dynamical behavior of the system. Based on numerical simulation results of bifurcation diagram, Lyapunov exponent spectrum, Poincare section, power spectrum and return map of the system, some basic dynamical behavior of this low-dimensional model are revealed, the new chaos system exhibits a route to chaos via a period doubling cascade (Feigenbaum's Scenario).
Key words: Hopf bifurcation     Lyapunov function     Strange Attractors     Chaos    
1 引言

对不可压缩Navier-Stokes方程的研究引起了物理学家和数学家们极大的兴趣. 目前流动现象相当的进展归功于计算机的出现和快速发展. 事实上,正是由于有了计算机,才使数值模拟非线性系统成为可能,即使支配流体流动的偏微分方程作傅立叶展开后截取有限项("模式") 得到的常微分方程组的规模巨大,也是可以用计算机来进行仿真的. 此领域第一个重要的贡献是Lorenz系统[1-2],它是由Saltzman提出的对流板之间的微分方程作三模式截断得到的方程. 十五年后,Navier-Stokes方程模式截断方法再次被考虑,成为更普遍的一种研究方法. Boldrighini和Franceschini开始了对具有周期边界条件的二维Navier-Stokes方程的模式截断的系列详细的研究,他们将平面正方形区域$T^2=[0,2\pi]\times [0,2\pi]$上不可压缩的Navier-Stokes方程的速度$u$,外力$f$和压力$p$ 进行傅立叶展开并截取其中的有限项,获得了五模和七模的类Lorenz方程组[3-5]. 对具有周期性边界条件的Navier-Stokes方程,文献[6-7]也作了一些相应的研究. 后来研究工作扩展到三维空间,1988年 Franceschini V,Inglese G 和 Tebaldi C在Commun Mech Phys上发表了三维空间上的有关Navier-Stokes方程五模截断的文章[8]; 1991年Franceschini V 和 Zanasi R在三维空间对Navier-Stokes方程傅立叶展开,进行七模截断后得到十四个非线性微分方程组成的方程组,并对其进行了详细的分析和讨论[9].

本文我们考虑磁流体动力学(简称MHD)方程的模式截断问题,把平面正方形区域上的磁流体动力学方程展成傅立叶级数形式,获得了全新的五模类Lorenz方程组. 理论分析并数值模拟了当参数值变化时这个新混沌系统的动力学行为,包括基本的动力学性质,分岔,和通向混沌的道路等. 我们的意图就是探讨流体的磁性对系统动力学行为的影响. 而且,对新混沌系统吸引子的存在性和全局稳定性进行了分析和讨论,这些理论可以应用于其它类似的模型. 考虑$T^2=[0,2\pi]\times [0,2\pi]$上的MHD方程[10]

$\left\{\begin{array}{l} \frac{\partial u}{\partial t}+(u\cdot \nabla) u-\lambda \Delta u+\nabla p +s\nabla (\frac{1}{2}B^2)-s(B\cdot \nabla)B=f,\\[3mm] \frac{\partial B}{\partial t}+(u\cdot \nabla) B-(B\cdot \nabla)u+\lambda_m \mbox{rot}(\mbox{rot} B)=0,\\[1mm] \mbox{div} u=0,\mbox{div} B=0,\\ \int_{T^2} u{\rm d}x=0, \int_{T^2} B{\rm d}x=0,\end{array}\right. $ (1.1)

其中$u$为速度场函数; $B$为磁场磁感应强度; $p$为液体压力; $\lambda =\frac{1}{ Re}$,$Re$为流动雷诺数;$\lambda_m =\frac{1}{R_m}$,$R_m$为磁雷诺数; $s=\frac{M^2}{{ Re}\cdot R_m}=\lambda\cdot \lambda_m\cdot M^2$,$M$ 为Hartman数,$f$ 代表流体的(周期)外力.

2 新五模类Lorenz方程组

我们把$u,B,p,f$ 展成如下傅立叶级数形式

$u(X,t)=\sum_{K\neq 0}{\rm e}^{{\rm i}K\cdot X}r_K\frac{K^\perp}{\mid K\mid},$ (2.1)
$B(X,t)=\sum_{K\neq 0}{\rm e}^{{\rm i}K\cdot X}b_K\frac{K^\perp}{\mid K\mid}, $ (2.2)
$f(X,t)=\sum_{K\neq 0}{\rm e}^{{\rm i}K\cdot X}f_K\frac{K^\perp}{\mid K\mid},$ (2.3)
$p(X,t)=\sum_{K\neq 0}{\rm e}^{{\rm i}K\cdot X}p_K(t),$ (2.4)

其中 $K=(h_1,h_2)$ 是具有整数分量的波向量,$K^\perp=(h_2,-h_1)$,$r_K=r_K(t),b_K=b_K(t),f_K=f_K(t),p_K=p_K(t)$ 均为时间 $t$ 的函数. 将(2.1)-(2.4)式代入到方程组(1.1),经过一系列运算得到如下MHD方程组的截断方程组为

$\left\{ \begin{array}{rl} \dot{r}_K=& -\lambda\mid K\mid^2r_K-{\rm i}\sum_{K\in L,K_1+K_2+K= 0}\frac{K_1^\perp\cdot K_2(K^2_2-K^2_1)}{2\mid K\mid\mid K_1\mid\mid K_2\mid}\overline{r}_{K_1}\overline{r}_{K_2}\\[4mm] & +s\cdot {\rm i}\sum_{K\in L,K_1+K_2+K= 0}\frac{K_1^\perp\cdot K_2(K^2_2-K^2_1)}{2\mid K\mid\mid K_1\mid\mid K_2\mid}\overline{b}_{K_1}\overline{b}_{K_2}+f_K,\\[4mm] \dot{b}_K=& -\lambda_m\mid K\mid^2b_K-{\rm i}\sum_{K\in L,K_1+K_2+K= 0}\frac{K_1^\perp\cdot K_2\mid K\mid}{2\mid K_1\mid\mid K_2\mid}(\overline{r}_{K_1}\overline{b}_{K_2}-\overline{b}_{K_1}\overline{r}_{K_2}),\end{array} \right.%$ (2.5)

其中 $L$ 为波向量集合,并且满足:若$K\in L$,则 $-K\in L$. 本文取

$L=\{\pm K_1,\pm K_2,\pm K_3,\pm K_4,\pm K_5\},$

其中 $K_1=(1,1),K_2=(2,1),K_3=(1,0),K_4=(3,-1) ,K_5=(2,-2)$,并且作下列变换

$ \begin{array}{l} r_{K_1}=2\sqrt{10}x_1,r_{K_2}=-2\sqrt{10}{\rm i}x_2,r_{K_3}=2\sqrt{10}x_3,r_{K_4}=-2\sqrt{10}{\rm i}x_4,r_{K_5}=2\sqrt{10}x_5,\\ b_{K_1}=2\sqrt{10}x_6,b_{K_2}=2\sqrt{10}{\rm i}x_7,b_{K_3}=2\sqrt{10}x_8,b_{K_4}=2\sqrt{10}{\rm i}x_9,b_{K_5}=2\sqrt{10}x_{10}. \end{array} $ (2.6)

对得到的方程组只施加外力$f_5$,令

$\lambda=\lambda_m=s=1,r=\mid f_{K5}\mid,$

代入到方程组(2.5)经大量计算得到如下类Lorenz方程组为

$\left\{ \begin{array}{ll} \dot{x}_1=-2x_1+4x_2x_3+4x_7x_8+2x_4x_5+2x_9x_{10},& (1)\\ \dot{x}_2=-5x_2-x_1x_3+x_6x_8 ,& (2)\\ \dot{x}_3=-x_3-3x_1x_2-3x_6x_7 ,& (3)\\ \dot{x}_4=-10x_4+6x_1x_5-6x_6x_{10},& (4)\\ \dot{x}_5=-8x_5-8x_1x_4-8x_6x_9+r ,& (5)\\ \dot{x}_6=-2x_6-2x_2x_8-2x_3x_7-2x_4x_{10}-2x_5x_9,~~&(6)\\ \dot{x}_7=-5x_7-5x_1x_8+5x_3x_6 ,& (7)\\ \dot{x}_8=-x_8+x_1x_7+x_2x_6 ,& (8)\\ \dot{x}_9=-10x_9-10x_1x_{10}+10x_5x_6 ,& (9)\\ \dot{x}_{10}=-8x_{10}+8x_1x_9+8x_4x_6 ,&(10) \end{array} \right.$ (2.7)

这里$x_i=x_i(t)\ (i=1,2,\cdots,10)$ 为谱展开系数. 与文献[3]中的五维模型比较,系统(2.7)的方程(6)-(10)源于磁场,并且 (2.7) 式具有下列对称性质

$(x_1,-x_2,-x_3,-x_4,-x_5,x_6,-x_7,-x_8,-x_9,-x_{10})\Leftrightarrow (x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10}).$
3 定常解及其稳定性

$F(x,r)=F(x_1,\cdots,x_{10},r)=\left( \begin{array}{c} -2x_1+4x_2x_3+4x_7x_8+2x_4x_5+2x_9x_{10} \\ -5x_2-x_1x_3+x_6x_8 \\ -x_3-3x_1x_2-3x_6x_7 \\ -10x_4+6x_1x_5-6x_6x_{10} \\ -8x_5-8x_1x_4-8x_6x_9+r \\ -2x_6-2x_2x_8-2x_3x_7-2x_4x_{10}-2x_5x_9 \\ -5x_7-5x_1x_8+5x_3x_6 \\ -x_8+x_1x_7+x_2x_6 \\ -10x_9-10x_1x_{10}+10x_5x_6 \\ -8x_{10}+8x_1x_9+8x_4x_6 \end{array} \right ),$ (3.1)

经计算我们获得$F(x,r)$的下列导数

$ D_xF(x,r)= \left[\begin{array}{cccccccccccc} -2&~ 4x_3~~& 4x_2&~ 2x_5~~& 2x_4& 0& 4x_8& ~~4x_7~~& 2x_{10}& 2x_9\\ -x_3& -5& -x_1& 0& 0& x_8& 0& x_6& 0& 0\\ -3x_2& -3x_1& -1& 0& 0&~ -3x_7~~& -3x_6& 0& 0& 0\\ 6x_5& 0& 0& -10& 6x_1& -6x_{10}& 0& 0& 0& -6x_6\\ -8x_4& 0& 0& -8x_1& -8& -8x_9& 0& 0& -8x_6& 0\\ 0& -2x_8& -2x_7& -2x_{10}& -2x_9& -2& -2x_3& -2x_2& -2x_5& -2x_4\\ -5x_8& 0& 5x_6& 0& 0& 5x_3& -5& -5x_1& 0& 0\\ x_7& x_6& 0& 0& 0& x_2& x_1& -1& 0& 0& \\ -10x_{10}& 0& 0& 0& 10x_6& 10x_5& 0& 0& -10& -10x_1&\\ 8x_9& 0& 0& 8x_6& 0& 8x_4& 0& 0& 8x_1& -8\\ \end{array} \right].$ (3.2)

$F(x,r)=0$,经大量运算可以求出(2.7)式的定常解(平衡点),下面给出定常解,并讨论其稳定性.

1) 对于 $0\leq r\leq R_1=8\sqrt{\frac{5}{3}}$,系统(2.7)只有唯一的定常解

$ {\rm (i)}X_1=(0,0,0,0,\frac{r}{8},0,0,0,0,0,). $

将其代入(3.2)式,我们得到此平衡点处的雅可比矩阵

$D_xF(x,r)\mid_{X_1}=\left[\begin{array}{cccccccccccc} -2&~ 0~~& 0& ~~\frac{2r}{8}~~& 0&~ 0~~& 0&~ 0~~& 0& 0\\ 0& -5& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& -1& 0& 0& 0& 0& 0& 0& 0\\ \frac{6r}{8}& 0& 0& -10& 0& 0& 0& 0& 0& 0 \\ 0& 0& 0& 0& -8& 0& 0& 0& 0& 0 \\ 0& 0& 0& 0& 0& -2& 0& 0& -\frac{2r}{8}& 0 \\ 0& 0& 0& 0& 0& 0& -5& 0& 0& 0 \\ 0& 0& 0& 0& 0& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& 0& \frac{10r}{8}& 0& 0& -10& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& -8 \end{array} \right]$ (3.3)

对应的特征方程为

$(\lambda+1)^2(\lambda+5)^2(\lambda+8)^2(\lambda^2+12\lambda+20-\frac{3r^2}{16}) (\lambda^2+12\lambda+20+\frac{5r^2}{16})=0.$ (3.4)

$0\leq r<R_1=8\sqrt{\frac{5}{3}}$时,矩阵(3.3)所有特征值均为负实部,所以平衡点$X_1$稳定,数值计算表明当$r$充分小时,它是稳定的全局吸引子.

2) 当 $R_1<r\leq R_2=16\sqrt{\frac{5}{3}}$时,系统(2.7)有3个平衡点: 此时 $X_1$ 变得不稳定,在临界点$R_1$$X_1$经树枝分岔产生如下2个新平衡点:

$ {\rm (ii)}X_{2,3}=(\pm \sqrt{\frac{5}{3}}\sqrt{\frac{r}{8\sqrt{\frac{5}{3}}}-1},0,0,\pm \sqrt{\frac{r}{8\sqrt{\frac{5}{3}}}-1},\sqrt{\frac{5}{3}},0,0,0,0,0). $

$ X_{2,3}$ 矩阵(3.2) 特征方程

$(\lambda^2+6\lambda+10-\frac{5r}{8\sqrt{\frac{5}{3}}}) (\lambda^3+20\lambda^2+\frac{12r}{\sqrt{\frac{5}{3}}}\lambda+640(\frac{r}{8\sqrt{\frac{5}{3}}}-1))$ (3.5)
$\times (\lambda^2+6\lambda+\frac{25r}{24\sqrt{\frac{5}{3}}}-\frac{10}{3}) (\lambda^3+20\lambda^2+\frac{56r}{3\sqrt{\frac{5}{3}}}\lambda+\frac{1600}{3}+\frac{320r}{24\sqrt{\frac{5}{3}}})=0.$

由于在此两点处雅可比矩阵所有特征值均为负实部,所以 $ X_{2,3}$ 是稳定的,数值结果表明任意选定的初值都被它们吸引,所以它们是全局吸引子. 当$r>R_2=16\sqrt{\frac{5}{3}}$时,他们变得不稳定.

3) 当 $r>R_2$,系统(2.7)有7个平衡点: 原来的 $X_1$,$X_{2,3}$,和如下的 $X_{4-7}$

${\rm (iii)} X_{4-7}=(\varepsilon \sqrt{\frac{5}{3}},\sigma\sqrt{\frac{\frac{3}{5}(\frac{r}{16})^2-1}{6}},-\sqrt{15}\varepsilon \sigma\sqrt{\frac{\frac{3}{5}(\frac{r}{16})^2-1}{6}},\varepsilon \frac{\sqrt{\frac{3}{5}r}}{16},\frac{r}{16},0,0,0,0,0),\\ \varepsilon =\pm 1,\sigma=\pm 1.$

此时前三个平衡点总是不稳定的. 当$r<R_3=65.951$时四个新平衡点是稳定的,数值计算表明它们吸引任意初值. 当 $r\geq R_3$时,由于一对复共轭特征值穿过虚轴四个新平衡点开始不稳定,在临界点 $R_3$ 处环绕平衡点$X_{4-7}$经由Hopf分岔产生4个稳定的周期轨道,在第5段我们将给出详细的讨论.

4 吸引子的存在性和全局稳定性分析
4.1 吸引子的存在性

下面我们证明系统(2.7)吸引子的存在性,设 $H=R^{10},u(t)=(x_1,\cdots,x_{10})$,作下列运算:

$(1)\times x_1+(2)\times x_2+(3)\times x_3+(4)\times x_4+(5)\times x_5+(6)\times x_6+(7)\times x_7+(8)\times x_8+(9)\times x_9+(10)\times x_{10}$

得到

$\begin{array}{l} \dot{x}_1x_1+2x^2_1+\dot{x}_2x_2+5x^2_2+\dot{x}_3x_3+5x^2_3+\dot{x}_4x_4+10x^2_4+\dot{x}_5x_5+8x^2_5 +\dot{x}_6x_6+2x^2_6\\ +\dot{x}_7x_7+5x^2_7+\dot{x}_8x_8+x^2_8+\dot{x}_9x_9+10x^2_9+\dot{x}_{10}x_{10}+8x^2_{10}=x_5\cdot r,\\ \end{array} $

$\frac{1}{2}\frac{\rm d}{{\rm d}t}(\sum^{10}_{i=1} x^2_i)+(2x^2_1+5x^2_2+x^2_3+10x^2_4+8x^2_5 +2x^2_6+5x^2_7+x^2_8+10x^2_9+8x^2_{10})=x_5\cdot r.$

假设 $\mid u(t)\mid^2=\sum\limits^{10}_{i=1} x^2_i$,由 Young 不等式[10-11]

$\frac{1}{2}\frac{\rm d}{{\rm d}t}(\sum^{10}_{i=1} x^2_i)+(2x^2_1+5x^2_2+x^2_3+10x^2_4+8x^2_5 +2x^2_6+5x^2_7+x^2_8+10x^2_9+8x^2_{10})\\ =x_5\cdot r\leq \frac{r^2}{4}+x^2_5.$

那么 $\frac{\rm d}{{\rm d}t}\mid u\mid^2+2\mid u\mid^2\leq \frac{r^2}{2}$. 利用 Gronwall 不等式[10-11]$\mid u\mid^2\leq \mid u(0)\mid^2{\rm e}^{-2t}+\frac{r^2}{4}(1-{\rm e}^{-2t})$,则 $\lim\limits_{t\rightarrow\infty}\sup\mid u(t)\mid^2\leq\frac{r^2}{4}.$ 因此有

$\lim\limits_{t\rightarrow \infty}\sup\mid u(t)\mid\leq \frac{r}{2}=\rho_0. $

如果 $\rho$充分大,则$B(0,\rho)$是泛函不变集和吸引集,因此系统(2.7)存在全局吸引子[10-11].

4.2 全局稳定性

为讨论系统(2.7)的全局稳定性,构造李雅普诺夫函数为[11-12]

$V(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10})=\sum^{10}_{i=1}x^2_italic>0.$ (4.1)

$V(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10})=k$,很明显,当$k$ 是一正常数时,上式表示$H$上的一球面,记为$E$.求$V$ 的导数,并利用(2.7)式得

$\frac{{\rm d}V}{{\rm d}t}=2\sum^{10}_{i=1}x_i\dot{x}_i \\ =-2(2x^2_1+5x^2_2+x^2_3+10x^2_4+8x^2_5+2x^2_6+5x^2_7+x^2_8+10x^2_9+8x^2_{10}-x_5r) \\ =-2[2x^2_1+5x^2_2+x^2_3+10x^2_4 +8(x_5-\frac{r}{16})^2+2x^2_6+5x^2_7+x^2_8+10x^2_9+8x^2_{10}-\frac{r^2}{32}]. $ (4.2)

显然,$2x^2_1+5x^2_2+x^2_3+10x^2_4+8(x_5-\frac{r}{16})^2+2x^2_6+5x^2_7+x^2_8+10x^2_9+8x^2_{10}=\frac{r^2}{32}$ 表示$H$中一椭球面,记此椭球面为$C$,由(4.2)式得

$\frac{{\rm d}V}{{\rm d}t}\left\{ \begin{array}{l} <,~~\mbox{在$C$外},\\ =0,~~\mbox{在$C$上},\\ >0,~~\mbox{在$C$内}. \end{array}\right. $

$k$取得充分大,$E$即可包围$C$. 这样,从上可知,在$C$ 外面,$\frac{{\rm d}V}{{\rm d}t}<0$ ,$V\frac{{\rm d}V}{{\rm d}t}<0$,由李雅普诺夫定理[11-12]得知,$E$外系统(2.7)的解轨线都将进入$E$内. 可见$E$ 就是类Lorenz系统(2.7)的捕捉区. 虽然这时系统平衡点$(X_i,i=1,\cdots,7)$都不稳定,但系统仍具有全局稳定性. 系统最终要收缩到捕捉区内,而区内又无收点,因此系统的轨线只能在区内不停的振荡. 于是轨线最终要在捕捉区内形成一个不变集合,这就是所谓的奇怪吸引子.

4.3 全局指数吸引集和正向不变集

下面我们讨论系统(2.7)的全局指数吸引集和正向不变集. 首先,我们给出一些基本定义,令 $X(t)=(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10})$,并且假设$X(t,t_0,X_0)$ 表示满足初始条件$X(t_0,t_0,X_0)=X_0$系统(2.7)的解,在不会混淆的情况下,该解我们简记为$X(t)$.

定义 4.1 如果存在常数 $L>0$,使得对于 $V(X_0)>L$,$V(X(t))>L$,都有$\lim\limits_{t\rightarrow +\infty} V(X)\leq L$,那么我们称 $\Omega =\{X\mid V(X(t))\leq L\}$ 是系统(2.7)的一个全局吸引集.

如果对任意的 $X_0\in \Omega $ 和任意的 $t>t_0$,都有$X(t,t_0,X_0)\in \Omega $,则称 $\Omega =\{X\mid V(X(t))\leq L\}$ 为正向不变集.

如果存在常数 $L>0$,$M>0$ 对任意的 $X_0\in R^{10}$,使得 $V(X_0)>L$,$V(X(t))>L$,有下列的指数估计不等式

$V(X(t))-L\leq (V (X_0)-L){\rm e}^{-M(t-t_0)},$

那么我们称 $\Omega =\{X\mid V(X(t))\leq L\}$ 是系统(2.7)的一个全局指数吸引集.

定理 4.1

$V(X)=\sum^{10}_{i=1}x^2_i>0,L=\frac{r^2}{15},$

$V(X_0)\geq L$,$V(X(t))\geq L$时,系统(2.7)有全局指数吸引集和正向不变集的估计式 $V(X(t))-L\leq (V(X_0)-L){\rm e}^{-(t-t_0)}$,从而 $\overline{\lim\limits_{t\rightarrow +\infty}} V(X(t))\leq L$. 也就是说 $\Omega =\{X\mid V(X(t))\leq L\}$ 是系统(2.7)的全局指数吸引集和正向不变集.

计算沿系统(2.7)的正半轨线$V(X)$关于时间的导数,有

$\frac{{\rm d}V}{{\rm d}t}=2\sum^{10}_{i=1}x_i\dot{x}_i \\ =-2(2x^2_1+5x^2_2+x^2_3+10x^2_4+8x^2_5+2x^2_6+5x^2_7+x^2_8+10x^2_9+8x^2_{10}-x_5r) \\ =-V+F(X),$ (4.3)

其中 $F(X)=-3x^2_1-9x^2_2-x^2_3-19x^2_4-15x^2_5-3x^2_6-9x^2_7-x^2_8-19x^2_9-19x^2_{10}+2x_5r$. 计算$F(X)$关于$X$的Lagrange极值,因为$F$为二次函数,其局部极大值为全局极大值. 令

$\frac{\partial F}{\partial x_i}=0,~~i=1,\cdots,10$

得极值点 $P_0=(0,0,0,\frac{r}{15},0,0,0,0,0,0)$,再求$F(X)$$P_0$点的二阶导数,$ \frac{\partial ^2F}{\partial x^2_i}<0,i=1,\cdots,10$,并且 $\frac{\partial ^2F}{\partial x_i\partial x_j}=0,i\neq j$,$\sup\limits_{X\in R^3}F(X)=F(X)\mid_{P_0}=\frac {r^2}{15}=L$. 由(4.3)式,我们有 $\frac{{\rm d}V}{{\rm d}t}\leq -V+L$,从而,当$V(X_0)\geq L$,$V(X(t))\geq L$时,有全局指数估计式 $V(X(t))-L\leq (V(X_0)-L){\rm e}^{-(t-t_0)},$ 对上式两边取上极限,有 $V(X(t))-L\leq (V(X_0)-L){\rm e}^{-(t-t_0)},$ 那么 $\overline{\lim\limits_{t\rightarrow +\infty}} V(X(t))\leq L$. 因此 $\Omega =\{X\mid V(X(t))\leq L\}$ 是系统(2.7)的全局指数吸引集和正向不变集.

5 数值仿真及分析

随着雷诺数$r$的增大,系统(2.7)的稳定性会发生变化,将出现Hopf分岔和混沌等非线性现象. 下面就来详细数值模拟系统(2.7)分岔和混沌发生的全过程.

1) 对于$r<R_3=65.951\cdots$,系统(2.7)定常解 $X_{4-7}$ 是稳定的,数值模拟显示任意选定的初值均被它们吸引,它们为全局吸引子(见图 1,2). 或许是数值仿真的局限性我们仅能画出两个稳定解的轨线图,下面其它的仿真图也是类似的.

图 1 $r=52.85$$X_{4-7}$ 是稳定的

图 2 $r=54.48$$X_{4-7}$ 是稳定的

2) 当 $r=R_3=65.951\cdots$ 时,由于$X_{4-7}$处系统的雅可比矩阵有一对共轭特征值穿过虚轴,使得经先前分岔产生的4个稳定的定常解 $X_{4-7}$ 变得不稳定,它们经由Hopf分岔产生稳定的周期轨道 $\xi^i_0(i=1,2,3,4)$(图 3). $\xi^i_0(i=1,2,3,4)$一直稳定到$r=R_4=68.749\cdots$,数值结果显示他们吸引随机选取的任何点. 正如Hopf定理所预言的,当$r\rightarrow R_3$时,如此轨道的周期趋于$T=\frac{2\pi}{\mid \lambda_0\mid}\approx 0.49247$($\lambda_0$ 是穿越虚轴的特征值的虚部).

图 3 $r=55.90$时周期轨道

3) 当$r=R_4=68.749\cdots$时,雅可比矩阵的特征值在$-1$处穿过单位圆,周期轨道 $\xi^i_0(i=1,2,3,4)$失去稳定性. 正如一般的分歧理论[7]所预言的,一种新的周期轨道出现了,数值结果显示前面四个轨道的每一个分岔进入一个新轨道$\xi^{i}_1(i=1,2,3,4)$,新轨道$\xi^{i}_1(i=1,2,3,4)$在形状上与前面的非常类似,它们仅缠绕两次. 随着参数$r$的增加,轨道数也增加,计算机模拟显示类似的分叉发生在 $r=R_5=68.852\cdots,r=R_6=68.858\cdots,r=R_7=68.876\cdots$,和 $r=R_8=68.882\cdots$. 每一次新轨道环绕不动点两次,因此,周期倍分岔轨道序列 $\xi^i_n(i=1,2,3,4)$出现了(见图 4),系统处于拟周期状态,按照Feigenbaum理论周期倍分岔序列将聚集(见表 1); 在相空间中另一个区域,一对周期轨道产生了 (一条稳定,一条不稳定),滞后现象(吸引子共存)出现了(图 7,8,15,16)).

图 4 $r=55.98$时周期倍分岔轨道 $\xi^i_n$

表 1 周期倍分岔点Ri和相对周期T(Ri)

4) 正如上面所讨论的,当 $r=R_4=68.749\cdots$ 时一种新的周期轨道出现了,前面四个轨道的每一个分岔进入一个新轨道$\xi^{i*}_0\ (i=1,2)$,新轨道缠绕$X_{4-7}$中的两个平衡点,而不是像前面仅绕一个平衡点,更精确的讲,有两个轨道缠绕$X_4$$X_6$,而另外两个轨道缠绕$X_5$$X_7$ (一条稳定,一条不稳定),经周期倍分岔导致的另外的轨道序列$\xi^{i*}_n\ (i=1,2)$出现了(见图 5-8).

图 5 $r=68.879$时的新轨道 $\xi^{i*}_n$

图 6 $r=68.881$时的新轨道 $\xi^{i*}_n$

图 7 $r=68.884$时轨道 $\xi^i_n$$\xi^{i*}_n$

图 8 $r=68.892$时轨道 $\xi^i_n$$\xi^{i*}_n$

5) 周期倍分岔轨道$\xi^i_n\ (i=1,2,3,4)$的吸引域随着$r$的增大(在很小的区间范围内)迅速扩张,当$r$位于这个小区间内时,滞后现象(不唯一)开始发生了,一些初值被新的稳定轨道$\xi^{i*}_n\ (i=1,2)$吸引,另外的一些被周期倍分叉轨道$\xi^{i}_n\ (i=1,2,3,4)$吸引(见图 7-16). 很快周期倍分岔轨道$\xi^i_n\ (i=1,2,3,4)$ 的吸引域似乎吞没了新轨道$\xi^{i*}_n\ (i=1,2)$ 所要定位的区域,暂态混沌似乎发生了(见图 15,16,41),滞后现象出现了. 由于磁性的影响,间歇性不同于流体五模系统的间歇性,周期倍分岔轨道$\xi^i_n\ (i=1,2,3,4)$为主导. 然后奇怪吸引子出现了,周期倍分岔过渡到混沌(图 4,9,12,15,16). 图 17-22给出了系统(2.7)在不同雷诺数$r$下的周期轨道,图 23-26给出了混沌区内一些新稳定轨道的相图.

图 9 $r=68.98$时间歇性

图 10 $r=68.99$时间歇性

图 11 $r=68.928$时轨道 $\xi^i_n$

图 12 $r=68.934$时轨道 $\xi^i_n$$\xi^{i*}_n$

图 13 $r=68.96$时轨道 $\xi^i_n$$\xi^{i*}_n$

图 14 $r=68.96$时轨道 $\xi^i_n$$\xi^{i*}_n$

图 15 $r=68.97$时奇怪吸引子

图 16 $r=68.98$时奇怪吸引子

图 17 r=72.54

图 18 $r=72.55$

图 19 $r=72.56$

图 20 $r=72.58$

图 21 $r=73.25$

图 22 $r=73.26$

图 23 $r=78.67$

图 24 $r=78.68$

图 25 $r=78.68$

图 26 $r=78.71$

6) 当$r=69.507\cdots$时轨道 $\xi^i_n\ (i=1,2,3,4)$ 消失,由拟周期分岔产生的轨道$\xi^{i*}_n\ (i=1,2)$$r$的增大继续失去稳定性而连续加倍,最终产生混沌吸引集,图 27-40显示了系统 (2.7)再次进入混沌状态.

图 27 $r=82.42$时极限环

图 28 $r=82.73$

图 29 $r=83.00$

图 30 $r=83.05$

图 31 $r=83.12$

图 32 $r=83.23$

图 33 $r=83.25$

图 34 $r=83.27$

图 35 $r=83.32$

图 36 $r=83.33$

图 37 $r=87.50$

图 38 $r=87.51$

图 39 $r=87.52$时奇怪吸引子

图 40 $r=87.53$时奇怪吸引子

7) 图 41描绘了系统(2.7)的状态变量$x_5$相对于雷诺数 $r$分岔图,图 42是对应的最大李雅普诺夫指数. 从图中可以看出最大李雅普诺夫指数大于零的区域与分岔图显示的混沌区域是一致的. 从分岔图 41我们能看出系统(2.7)转变到混沌整个过程,暂态混沌似乎发生在$r=67.257\cdots$,当 $r=68.753\cdots$ 系统进入周期倍分岔状态,从周期4到周期8,再到周期16,最终到达混沌. 在$r=72.503\cdots$$r=81.875\cdots$时也有类似的情况发生.

图 41 分叉图

图 42 最大Lyapunov指数

8) 图 43-46给出了系统(2.7)在一些雷诺数$r$下的庞加莱截面,图 47-50展示了系统(2.7)在不同雷诺数$r$下的返回映射,图 51,52描绘了系统(2.7)在$r=78.69$$r=84.63$下的功率谱,它们均表明了系统的混沌特征.

图 43 $r=78.68$时Poincare截面

图 44 $r=78.68$时Poincare截面

图 45 $r=84.51$

图 46 $r=84.53$

图 47 $r=78.68$时返回映射

图 48 $r=78.68$时返回映射

图 49 $r=84.51$

图 50 $r=84.53$

图 51 $r=78.69$时功率谱

图 52 $r=84.63$时功率谱

9) 从分岔图 41中我们发现混沌区包含明显的周期窗口,在某些雷诺数奇怪吸引子和极限环交替出现. 通过精细计算我们获得如下细节:在$r=82.04\cdots$ 奇怪吸引子收缩成对称的极限环(图 27-30),当$r=82.72\cdots$时由周期倍分岔产生的周期分岔轨道$\xi^{i*}_n\ (i=1,2)$快速连续演化,最终产生混沌吸引集(图 33-40).

10) 图 53-58给出了在不同雷诺数$r$下系统(2.7)的一些状态轨线.

图 53 $r=78.68$时状态轨线

图 54 $r=78.69$时状态轨线

图 55 $r=78.68$

图 56 $r=78.69$

图 57 $r=84.51$

图 58 $r=84.63$

11) 当$r\geq 155.986$ 系统不在出现稳定的周期轨道,模型展现了湍流行为,图 59-64 显示了在高雷诺数下系统的湍流状态.

图 59 $r=155.98$

图 60 $r=162.58$

图 61 $r=173.56$

图 62 $r=194.24$

图 63 $r=205.36$时奇怪吸引子

图 64 $r=223.37$时奇怪吸引子
6 结论

由于我们的五模磁流体系统动力学现象较为复杂,因而有必要给出系统的总结,首先重新定义$r$的几个临界值,$r'_1=R_3=65.951,r'_2=68.749,r'_3=68.882,r'_4=69.507,r'_5=155.986,$ 总结以上的数值模拟结果得到:

(a) 当 $0<Re\leq r'_1$时,系统(2.7)只有一个平衡点.

(b) 当 $r'_1<r\leq r'_4$时,系统(2.7)通过无穷分岔序列产生四个对称的无穷序列周期倍分岔轨道$\xi^i_n(i=1,2,3,4)$; 每一个都是先前轨道的周期加倍,在$r'_2<r\leq r'_3$内我们发现了几个周期倍分岔的临界值,并且它们满足Feigenbaum原理.

(c) 当 $r'_3<Re\leq r'_4$时,进一步的无限序列分岔引起了两个对称的结构更复杂的无穷序列周期倍分岔轨道$\xi^{i*}_n(i=1,2)$,它们对称放置,但结构更复杂. 然后系统产生滞后现象,周期倍分岔轨道$\xi^i_n$$\xi^{i*}_n$随r增大连续加倍最终产生混沌吸引集,由于磁性的影响,轨道$\xi^i_n$是产生混沌的主要因素,这种情况与文献[3]中流体五模系统的恰好相反. (d) 当$r'_4<r\leq r'_5$时,轨道$\xi^i_n(i=1,2,3,4)$消失,周期轨道$\xi^{i*}_n(i=1,2)$ 变得越来越不稳定,经由周期倍分岔到达混沌,两个对称的"奇怪"吸引子类似于文献[1]中得Lorenz模型中的"湍流" .

(e) 当$r\geq r'_5$时,系统再无周期状态出现,模型展现了湍流行为. 系统的状态完全随机的,并且敏感依赖于初始条件. 由于所有的数值结果执行到$r=8000$,一直显示随机行为,我们认为混沌可能保持到 $r$趋于无穷大,并且在高雷诺数$r$的情形下没有稳定的吸引周期轨道存在,这是与文献[3]中五模系统的显著差别,即由于磁性的影响我们的磁流体五模动力系统没有复制文献[3]中五模系统的本质特征.

为了揭示流动行为的本质特征,低模分析方法是非常有意义的. 本文给出并研究了平面正方形区域上磁流体新五模类Lorenz系统,研究和数值模拟了此新混沌系统的基本动力学性质,分岔,通向混沌的道路等动力学行为,我们发现了拟周期分岔到达混沌,阵发混沌,周期窗口和滞后现象等. 借助于最大李雅普诺夫指数,庞加莱截面,返回映射和功率谱显示了系统的混沌特征. 而且全局吸引子的存在性证明和全局稳定性分析在相关文献中没有涉及,而且我们的讨论方法对其它类似的相关文献中的模型也是适用的.

参考文献
[1] Edward N, Lorenz. Deterministic nonperiodic flow. Journal Atmospheric Sciences, 1963, 20: 130–141. DOI:10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2
[2] Hilborn R C. Chaos and Nonlinear Dynamics. Oxford Unit:Oxford Univ Press, 1994: 10–30.
[3] Boldrighini C, Franceschini V. A five-dimensional truncation of the plane incompressible Navier-Stokes equations. Communications in Mathematical Physics, 1979, 64: 159–170. DOI:10.1007/BF01197511
[4] Franceschini V, Tebaldi C. A seven-modes truncation of the plane incompressible Navier-Stokes equations. Journal of Statistical Physics, 1981, 25(3): 397–417. DOI:10.1007/BF01010796
[5] Franceschini V, Tebaldi C. Breaking and disappearance of tori. Commun Math Phys, 1984, 94: 317–329. DOI:10.1007/BF01224828
[6] 王贺元. Navier-Stokes五模类Lorenz方程组的动力学行为及数值仿真. 应用数学与计算数学学报 , 2010, 24 (2): 13–22.
Wang H Y. The dynamical behavior of a five-modes lorenz equations of the Navier-Stokes equations for a two-dimensional incompressible fluid on a torus. Communication on Applied Mathematics and Computation, 2010, 24(2): 13–22.
[7] Wang H Y. The dynamical behavior and the numerical simulation of a five-mode system of the Navier-Stokes equations for a two-dimensional incompressible fluid on a torus. Journal of Applied Mathematics and Physics, 2016, 4(7): 1245–1253. DOI:10.4236/jamp.2016.47130
[8] Franceschini V, Inglese G, Tebaldi C. A five-mode truncation of the Navier-Stokes equations on a three-dimensional torus. Commun Mech Phys, 1988, 64: 35–40.
[9] Franceschini V, Zanasi R. Three-dimensional Navier-stokes equations trancated on a torus. Nonlinearity, 1992, 5(1): 189–209. DOI:10.1088/0951-7715/5/1/008
[10] Teman R. Infinite Dimensional Dynamic System in Mechanics and Physics. New York: Springer-Verlog, 2000.
[11] 李开泰, 马逸尘. 数理方程HILBERT空间方法(下). 西安: 西安交通大学出版社, 1992.
Li K T, Ma Y C. Hilbert Space Method for Mathematical and Physics Equations. Xi'an: Xi'an Jiaotong University Press, 1992.
[12] 刘秉正, 彭建华. 非线性动力学. 北京: 高等教育出版社, 2004: 406-415.
Liu B Z, Peng J H. Nonlinear Dynamics. Beijing: Advanced Education Press, 2004: 406-415.