流动现象是自然界及人类生产科研活动中最为常见的一种物理现象, 流动稳定性是流动现象最为关键的问题.作为流动现象应普遍遵循的Navier-Stokes方程是一种典型的非线性偏微分方程, 刻划着流体的运动规律, 如大气运动、海洋流动、轴承润滑、透平机械内部流动等, 研究它对人们认识和控制湍流至关重要. 1963年美国气象学家E.Lorenz在研究大气对流时, 首次给出了著名的Lorenz方程^[9].所采用的方法是对Navier-Stokes方程和热传导方程进行傅立叶级数展开, 截取级数的前三项, 得到三模的Lorenz系统. 20世纪后期Valter Franceschini又在此方向上进一步扩展, 多次和其他学者合作, 将二维正方形区域 $T^2=[0, 2\pi]\times [0, 2\pi]$上不可压缩的Navier-Stokes方程

$ \begin{eqnarray} &&\frac{\partial u}{\partial t}+(u\cdot \nabla) u=-\nabla p +f +\nu \triangle u, \\ &&{\rm div} u=0, \nonumber\\ &&\int_{T^2} udX=0\nonumber \end{eqnarray} $

(1.1)

(其中 $u$为速度场函数; $p$为液体之间的压力; $f$为外力场函数, $\nu$为动力粘性系数)进行傅立叶展开并截取其中的有限项, 得出五模和七模或者任意模的非线性微分方程组(见文献[1-4]), 讨论当雷诺数变化时方程组解的动力学行为.这种截断后来被扩展到三维空间, 1988年V. Franceschini, Inglese和Tebaldi在Commun. Mech. Phys.上发表了三维空间上的有关Navier-Stokes方程五模截断的文章^[7]; 1991年Franceschini和Zanasi在三维空间上对此方程傅立叶展开, 进行七模截断后得到十四个非线性微分方程组成的方程组, 随后又对这个复杂的方程组进行了详细的讨论^[3].国内王贺元等人选取不同的截断模式, 并把这方面的研究扩展到磁流体, 得到相应类Lorenz方程组并分析了系统的动力学行为^{[10, 11]}. Franceschini在1981年给出的一个七模类Lorenz方程组^[1], 讨论了这个七模模型定常解的线性稳定性, 并对分歧行为进行了数值模拟.本文对此模型的动力学行为进行深入的分析和探讨, 证明了该模型吸引子的存在性, 并讨论了其全局稳定性, 从而在理论上保证了数值模拟的有效性, 并且数值模拟了分歧和混沌吸引子的发生过程.

下面对二维区域 $[0, 2\pi]\times [0, 2\pi]$上Navier-Stokes方程进行傅立叶展开.即对速度函数 $u$, 外力场函数 $f$和流体之间的压力 $p$进行如下傅立叶展开

$ \begin{eqnarray} &&u(X, t)=\sum_{K\neq 0}e^{iK\cdot X}r_K\frac{K^\perp}{\mid K\mid}, \end{eqnarray} $

(2.1)

$ \begin{eqnarray} &&f(X, t)=\sum_{K\neq 0}e^{iK\cdot X}f_K\frac{K^\perp}{\mid K\mid}, \end{eqnarray} $

(2.2)

$ \begin{eqnarray} &&p(X, t)=\sum_{K\neq 0}e^{iK\cdot X}p_K(t), \end{eqnarray} $

(2.3)

其中 $K=(h_1, h_2)$是波向量, $K^\perp=(h_2, -h_1)$, $r_K=r_K(t)$为时间 $t$的函数.将(2.1)-(2.3) 式代入到方程组(1.1), 经过一系列运算得到如下形式的微分方程组^{[1, 2]}

$ {\dot r_K} = - i\sum\limits_{{K_1} + {K_2} + K = 0,{K_1},{K_2} \in L} {\frac{{K_1^ \bot \cdot {K_2}(K_2^2 - K_1^2)}}{{2\mid K\mid \mid {K_1}\mid \mid {K_2}\mid }}} {\bar r_{{K_1}}}{\bar r_{{K_2}}} - \nu \mid K{\mid ^2}{r_K} + {f_K}, $

(2.4)

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

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

其中

$ \begin{array}{l} K_1=(1, 1), K_2=(3, 0), K_3=(2, -1), K_4=(1, 2), \\ K_5=(0, 1), K_6=(1, 0), K_7=(1, -2), \\ -K_1=(-1, -1), -K_2=(-3, 0), -K_3=(-2, 1), -K_4=(-1, -2), \\ -K_5=(0, -1), -K_6=(-1, 0), -K_7=(-1, 2).\\ \end{array} $

(2.5)

在 $\nu=1$时, 分别令 $K$为 $K_1, K_2, K_3, K_4, K_5, K_6, K_7$, 代入到方程组(2.4) 经大量计算, 利用实条件 $r_K=-\overline{r}_{-K}$作代换

$ r_{K_1}=x_1, r_{K_2}=-ix_2, r_{K_3}=x_3, r_{K_4}=ix_4, r_{K_5}=x_5, r_{K_6}=ix_6, r_{K_7}=ix_7, $

对得到的方程组只施加外力 $f_3$, 且使雷诺数 ${\rm Re}=\sum\limits^{7}_{k=1}\mid f_K\mid=\mid f_3\mid$, 则截得类Lorenz方程组为

$ \left\{ \begin{array}{l} \dot{x}_1=-2x_1+4\sqrt{5}x_2x_3+4\sqrt{5}x_4x_5, &(1)\\ \dot{x}_2=-9x_2+3\sqrt{5}x_1x_3, &(2)\\ \dot{x}_3=-5x_3-7\sqrt{5}x_1x_2+9x_1x_7+{\rm Re}, &(3)\\ \dot{x}_4=-5x_4-\sqrt{5}x_1x_5, &(4)\\ \dot{x}_5=-x_5-3\sqrt{5}x_1x_4+5x_1x_6, &(5)\\ \dot{x}_6=-x_6-5x_1x_5, &(6)\\ \dot{x}_7=-5x_7-9x_1x_3, &(7)\\ \end{array} \right. $

(2.6)

这里 $x_i=x_i(t)(i=1, 2, \cdots, 7)$为谱展开系数.截得了七模非线性微分方程组的形式和Lorenz方程组相似, 称其为类Lorenz系统.

由于系统(2.6) 在平衡点处Jacobi矩阵与时间无关, 故李雅普诺夫矩阵的特征指数就是Jacobi矩阵的特征值的实部^{[5, 6]}, 它是刻画吸引子性质的重要指标, 尤其对混沌吸引子更为重要.下面对类Lorenz方程组(2.6) 线性化, 然后根据各个平衡点的李雅普诺夫矩阵的特征指数的变化来讨论平衡点的稳定性.令

$ F(X, {\rm Re})=F(x_1, x_2, x_3, x_4, x_5, x_6, x_7, {\rm Re})=\left( \begin{array}{l} -2x_1+4\sqrt{5}x_2x_3+4\sqrt{5}x_4x_5\\[6pt] -9x_2+3\sqrt{5}x_1x_3 \\[6pt] -5x_3-7\sqrt{5}x_1x_2+9x_1x_7+{\rm Re}\\[6pt] -5x_4-\sqrt{5}x_1x_5\\[6pt] -x_5-3\sqrt{5}x_1x_4+5x_1x_6\\[6pt] -x_6-5x_1x_5\\[6pt] -5x_7-9x_1x_3\\[6pt] \end{array} \right ), $

(3.1)

对 $F(X, {\rm Re})$关于 $X$求导数得到如下李雅普诺夫矩阵

$ D_XF(X, {\rm Re})=\left[ \begin{array}{lllllllllrrrrrrrr} -2& 4\sqrt{5}x_3& 4\sqrt{5}x_2& 4\sqrt{5}x_5& 4\sqrt{5}x_4& 0& 0\\ 3\sqrt{5}x_3&-9& 3\sqrt{5}x_1& 0& 0& 0& 0\\ 9x_7-7\sqrt{5}x_2& 7\sqrt{5}x_1&-5& 0& 0& 0& 9x_1\\ -\sqrt{5}x_5& 0& 0&-5& -\sqrt{5}x_1& 0& 0\\ 5x_6-3\sqrt{5}x_4& 0& 0& -3\sqrt{5}x_1& -1& 5x_1& 0\\ -5x_5& 0& 0& 0& -5x_1& -1& 0\\ -9x_3& 0& -9x_1& 0& 0& 0& -5\\ \end{array} \right], $

(3.2)

由 $F(X, {\rm Re})=0$求出(2.6) 式的平衡点, 下面根据Liapunov矩阵的特征指数的变化情况讨论各平衡点的稳定性(具体参见文献[1]).

1) 当 $0\leq {\rm Re}\leq R_1=\sqrt{\frac{15}{2}}$时, 方程组(2.6) 仅有唯一的一个平衡点

$ (0, 0, \frac{\rm Re}{5}, 0, 0, 0, 0), $

(1)

它是稳定的全局吸引子.

2) $R_1 < {\rm Re}\leq R_2=R_1+(\frac{3762}{25\sqrt{30}})\simeq 30.2124$时, 方程组(2.6) 有三个平衡点, 其中包括平衡点(1) 及如下的两点 $P_{\pm}$

$ (\sigma 5\sqrt{6}\sqrt{\frac{(2{\rm Re}-\sqrt{30})}{836\sqrt{30}}}, \sigma 5\sqrt{\frac{(2{\rm Re}-\sqrt{30})}{836\sqrt{30}}}, \sqrt{\frac{3}{10}}, 0, 0, 0, -\sigma \frac{27}{\sqrt{5}}\sqrt{\frac{(2{\rm Re}-\sqrt{30})}{836\sqrt{30}}}), $

(2)

其中 $\sigma=\pm1$, 此时平衡点(1) 变得不稳定, 两个新平衡点(2) 是稳定的.

3) 当 ${\rm Re}>R_2$时, 所有平衡点(1), (2) 都是不稳定的.

耗散动力系统的混沌行为是由于存在着一个复杂的吸引子而引起的^[8], 而这个吸引子就是系统的所有轨道当时间趋于无穷时收敛到的集合, 可能是一个分形或康托集或康托集和一个区间的乘积.很自然地这个“吸引子”就成为数学上用来描述观察到的不稳定流的对象, 它的复杂结构就是导致观察到的混沌现象的原因.因此, 研究吸引子的存在性和数值模拟就成为一个重要的问题.下面就来证明系统(2.6) 的吸引子存在性.

取 $H=R^7, u(t)=(x_1, \dots, x_7)$, 对Navier-Stokes方程的七模类Lorenz方程组(2.6) 作如下运算

$ (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 $

得

$ \begin{eqnarray*}&&\dot{x}_1x_1+2x^2_1+\dot{x}_2x_2+9x^2_2+\dot{x}_3x_3+5x^2_3+ \dot{x}_4x_4+5x^2_4\\ &&+\dot{x}_5x_5+x^2_5+\dot{x}_6x_6+x^2_6+\dot{x}_7x_7+5x^2_7=x_3{\rm Re}, \end{eqnarray*} $

因此有

$ \frac{1}{2}\frac{d}{dt}(x^2_1+x^2_2+x^2_3+x^2_4+x^2_5+x^2_6+x^2_7)+(2x^2_1+9x^2_2+5x^2_3+5x^2_4+x^2_5\\ +x^2_6+5x^2_7)=x_3{\rm Re}. $

令 $|u(t)|^2=\sum\limits^7_{i=1} x^2_i$, 利用Young不等式^[10]得

$ \begin{eqnarray*}&&\frac{1}{2}\frac{d}{dt}(x^2_1+x^2_2+x^2_3+x^2_4+x^2_5+x^2_6+x^2_7)+ (2x^2_1+9x^2_2+5x^2_3+5x^2_4+x^2_5+x^2_6+5x^2_7)\\ &=&x_3{\rm Re}\leq \frac{{\rm Re}^2}{4}+x^2_3, \end{eqnarray*} $

所以 $\frac{d}{dt}\mid u\mid^2+2\mid u\mid^2\leq \frac{{\rm Re}^2}{2}$.由Gronwall不等式^[10]得

$ \mid u\mid^2\leq \mid u(0)\mid^2e^{-2t}+\frac{{\rm Re}^2}{4}(1-e^{-2t}), $

因此有

$ \mathop {\lim }\limits_{t \to \infty } \sup \mid u(t){\mid ^2} \le \frac{{{\rm{R}}{{\rm{e}}^2}}}{4}. $

故有

$ \mathop {\lim }\limits_{t \to \infty } \sup \mid u(t)\mid \le \frac{{{\rm{Re}}}}{2} = {\rho _0}. $

记 $\sum =B(0, \rho)$, 其中 $\rho\geq \rho_0$充分大, 则 $\sum =B(0, \rho)$是泛函不变集和吸引集^[12], 因此系统(2.6) 存在全局吸引子^[12].

非线性系统具有全局稳定性时, 其轨线所收敛的单连通闭区域称为系统的捕捉区.只要能证明捕捉区的存在, 不论其中的定常解是否稳定, 系统均具有全局稳定性.而研究系统的全局稳定性主要借助于李雅普诺夫函数方法^{[5, 6]}.李雅普诺夫函数方法的基本思想是构造一个函数, 然后利用它的性质和这个函数沿方程(2.6) 的轨线方向的全导数的性质以确定(2.6) 式平衡点的稳定性, 以确定系统的捕捉区.

对系统(2.6) 构造李雅普诺夫函数为

$ V({x_1},{x_2},{x_3},{x_4},{x_5},{x_6},{x_7}) = \sum\limits_{i = 1}^7 {x_i^2} > 0. $

令 $V(x_1, x_2, x_3, x_4, x_5, x_6, x_7)=K$, 很明显, 当 $K$是一正常数时, 上式表示 $H$上的一球面, 记为 $E$.求 $V$的导数, 并利用(2.6) 式

$ \begin{eqnarray} \frac{dV}{dt}&=&2x_1\dot{x}_1+2x_2\dot{x}_2+2x_3\dot{x}_3+2x_4\dot{x}_4+2x_5\dot{x}_5+2x_6\dot{x}_6+2x_7\dot{x}_7\nonumber\\ &=&-2(2x^2_1+9x^2_2+5x^2_3+5x^2_4+x^2_5+x^2_6+x^2_7-x_3{\rm Re})\nonumber\\ &=&-2[2x^2_1+9x^2_2+5(x_3-\frac{\rm Re}{10})^2+5x^2_4+x^2_5+x^2_6+5x^2_7-\frac{{\rm Re}^2}{100}]. \end{eqnarray} $

(4.1)

显然 $2x^2_1+9x^2_2+5(x_3-\frac{\rm Re}{10})^2+5x^2_4+x^2_5+x^2_6+5x^2_7=\frac{{\rm Re}^2}{100}$表示一椭球面, 记此椭球面为 $C$, 由(4.1) 式得

$ \frac{dV}{dt}\left\{ \begin{array}{l} < 0,\quad \mbox{在$C$ 域以外,}\\ =0,\quad \mbox{在$C$ 上,}\\ >0,\quad \mbox{在$C$ 域内.}\\ \end{array}\right. $

于是若把 $K$取得充分大, $E$即可包围 $C$.这样从式(4.1) 式可知在 $C$外面, $\frac{dV}{dt} < 0$, $V\frac{dV}{dt} < 0$, 由李雅普诺夫定理^[5]的分析得知 $E$外(2.6) 式的解轨线都将进入 $E$内.可见 $E$就是类Lorenz系统(2.6) 的捕捉区.虽然这时系统平衡点(1), (2) 都不稳定, 但系统仍具有全局稳定性:系统最终要收缩到捕捉区内, 而区内又无收点, 因此系统只能在区内不停的振荡.于是轨线最终要在捕捉区内形成一个不变集合, 这就是所谓的吸引子.人们称混沌运动这种具有独特性质和结构的吸引子为奇怪吸引子.它是整体稳定性和局部不稳定性一对矛盾的结合体.其具体形式如何呢？下面就来数值模拟系统(2.6) 的奇怪吸引子.

随着雷诺数的增大, Lorenz方程组(2.6) 的稳定性发生了变化, 出现了Hopf分岔和混沌等非线性现象.下面就来详细数值模拟系统(2.6) 从分岔到混沌的全过程.

1) 当 ${\rm Re} < R_2=30.2123$时, 系统(2.6) 的新平衡点是稳定的, 解轨线为螺旋线(如图 1, 2).

2) 通过数值计算得方程组(2.6) 在 ${\rm Re}=R_2$时平衡点 $p_{\pm}$处的李雅普诺夫矩阵的一对复共轭特征值穿越虚轴, 其实部由负变正, 因而系统(2.6) 发生了Hopf分岔.从不稳定平衡点 $p_{\pm}$分叉出闭轨线.如图 3-6.

3) 当 ${\rm Re}=71.31$时, 平衡点 $p_{\pm}$处分叉出闭轨线开始不稳定, 分叉处环面, 如图 7, 8.

4) 当 ${\rm Re}$进一步增大时出现了滞后现象(各种吸引子共存), 如图 9-13.当 ${\rm Re}=248.23$时系统发生混沌, 出现奇怪吸引子, 如下图 14-16分别给出了不同雷诺数时奇怪吸引子的大体状态.通过数值计算表明系统在高雷诺数下一直是混沌状态, 这与文献[1]的结论是一致的.

5) 图 17, 18分别给出了系统分岔图和最大李雅普诺夫指数, 从分岔图 17表明:当 ${\rm Re} < 71.31$时, 系统是稳定的, 当 ${\rm Re}=71.31$时, 系统开始不稳定, 分叉处环面, 之后系统出现滞后现象, 当 ${\rm Re}=248.23$时系统发生混沌, 出现奇怪吸引子, 算到 ${\rm Re}=1000$系统始终是混沌状态, 在高雷诺数下系统处于湍流状态, 这一点也与lorenz系统有明显的区别. 图 18中给出的最大李雅普诺夫指数与分岔图 17是相符的.

6) 图 19给出了系统的庞加莱截面( ${\rm Re}=252.41$), 图 20给出了系统的功率谱( ${\rm Re}=255.24$), 它们均表明了系统的混沌运动特征.

本文首先给出了七模类Lorenz方程组的推导过程, 对此方程组线性化稳定性分析进行了简单介绍.然后证明了此方程组全局吸引子的存在性, 并对其全局稳定性进行了分析和讨论, 最后数值模拟了雷诺数变化时系统经由不变环面的失稳到达混沌的过程, 运用分岔图、最大李雅普诺夫指数、庞加莱截面和功率谱揭示了系统混沌行为的普适特征.