充液系统中液体晃动的全局分岔和混沌

Global Bifurcation and Chaos of Liquid Sloshing in Liquid-Filled System

摘要: 本文以在各类工程中广泛应用的充液系统为研究对象，建立了柱形储液箱中液体晃动模态的非线性动力学方程组。通过积分，得到了该类系统的同宿轨道的解析表达式，并利用Melnikov方法对此系统的全局动力学行为进行分析，计算出系统在扰动条件下进入混沌的阈值，表明此类系统在受到较大扰动时，将发生全局分岔和混沌现象。同时，文中给出的算例证明了理论分析结果的正确性。

Abstract: The nonlinear governing equations of the liquid sloshing modals in cylindrical storage tank are established. Through integration, the analytical expressions of the homoclinic orbits of this kind of system are got. Using the Melnikov method, the global dynamical behaviors are analyzed. The threshold of the chaotic motion is calculated, which shows that there exists global bifurcation and chaos when the system is subjected to greater disturbance. Meanwhile, numerical simulations are also given, which confirm the analytical results.

1. 引言

充液系统的非线性晃动动力学作为动力学与控制学科及其交叉学科的重要研究领域之一，广泛应用于工业、土木工程、汽车工程、船舶动力学、以及航空航天工业中 [1] [2] [3] 。以航空航天工程为例，液体晃动及其最优控制始终是航空航天飞行器总体设计的核心问题之一，迄今为止国内外航天力学工作者对充液系统的耦合晃动问题进行了广泛的研究 [4] [5] [6] [7] 。由于液体晃动问题可造成灾难性的破坏事故，使得液体晃动动力学的深入研究更加必要，同时，由于此问题存在大量的非线性因素，也使得此问题的理论研究更加复杂。

为建立有效的数学模型以研究充液系统的晃动动力学行为，尹立中 [8] ，马兴瑞等 [9] 将液体的速度势函数在自由液面处展开，利用Lagrange函数，导出了广义模态坐标下的非线性方程组，对系统的局部动力学行为进行理论分析。随着动力学理论的发展，全局动力学和混沌理论被广泛的应用于各种工程系统中 [10] [11] 。因此，将混沌行为的研究应用于充液系统是有必要和有意义的。

本文利用Melnikov方法对液体晃动方程的全局动力学行为 [12] 进行分析。Melnikov方法在解析分析非线性系统的混沌运动方面有着广泛的应用，主要分为单一扰动的情形 [13] ，或同频多扰动 [14] [15] 的情况。本文主要就单一扰动的情况进行分析，同频或非同频多扰动情况应用本文方法同样可以进行类似的分析。

2. 模型的建立

系统模型如图1所示，圆柱形储液箱半径为a，内部液体深度为h，定义OXYZ为绝对坐标系，局部坐标系 $o\theta rz$ 为柱坐标系。设液体燃料无旋，且与储液箱体接触角不变，根据速度势理论可得到液体晃动方程 [8] [9]

$\frac{\partial \Phi }{\partial t}+\frac{1}{2}{\left(\nabla \Phi -\stackrel{\dot{}}{U}-\Omega \times r\right)}^2-\frac{1}{2}{\left(\Omega \times r\right)}^2+\Pi \mp \frac{2\sigma \bar{H}}{\rho }=0$ (1)

其中， $\Phi \left(r,\theta ,z,t\right)$ 为无旋流体的速度势函数， $\stackrel{\dot{}}{U}$ 为储液箱平动速度矢量， $\Omega$ 为储液箱转动角速度， $r$ 为转动半径， $\Pi$ 为势函数， $\sigma$ 为黏滞系数， $\rho$ 为液体密度， $\bar{H}$ 为液面平均主曲率采用接触角不变假设时的常数，t为时间， $\nabla$ 为散度算子。

不考虑储液箱体的转动，对液体的速度势函数进行模态展开 [9] ，其模态示意图如图2所示，得到描述液体晃动规律的方程组

${\ddot{\eta }}_m+2{\varsigma }_m{\stackrel{\dot{}}{\eta }}_m+{\upsilon }_m^2{\eta }_m+{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\sigma }_{ijm}{\eta }_i{\eta }_j}}+{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\displaystyle \underset{k=1}{\overset{N}{\sum }}{\sigma }_{ijkm}}}}{\eta }_i{\eta }_j{\eta }_k+{\tilde{\alpha }}_m\ddot{X}+{\tilde{e}}_m\ddot{Y}=0,$ $m=1,2,3,\cdots$ (2)

Figure 1. System model

图1. 系统模型

Figure 2. Modal of sloshing liquid

图2. 液体晃动模态示意图

其中， ${\eta }_m$ ( $m=1,2,3,\cdots$ )为各阶广义模态位移， ${\varsigma }_m$ ( $m=1,2,3,\cdots$ )为各阶模态阻尼， ${\upsilon }_m$ ( $m=1,2,3,\cdots$ )为各阶模态频率， ${\sigma }_{ijm}$ 、 ${\sigma }_{ijkm}$ 、 ${\tilde{\alpha }}_m$ 、 ${\tilde{e}}_m$ ，( $m=1,2,3,\cdots$ )为模态耦合系数。 $\ddot{X}$ 和 $\ddot{Y}$ 分别为X和Y方向上的扰动，上标一点为相应变量的一阶导数，二点为二阶导数。

设系统的扰动形式为 $\ddot{X}=S_1\cos \left({\Omega }_{1}t\right)$ ， $\ddot{Y}=0$ ，取液体晃动方程组中的第一阶主模态作为简化，得到

${\ddot{\eta }}_1+{\upsilon }_1^2{\eta }_1=K_{31}{\stackrel{\dot{}}{\eta }}_1+K_{32}{\eta }_1^2+K_{33}{\eta }_1^3+K_{34}\cos \left({\Omega }_{1}t\right)$ (3)

其中， ${\eta }_1$ 为液体晃动的第一阶主模态， ${\upsilon }_1$ 为一阶模态频率， $S_1$ 为X方向扰动的振幅， ${\Omega }_1$ 为X方向扰动的频率， $K_{31}=-2{\varsigma }_1$ ， $K_{32}=-{\sigma }_{111}/{\mu }_{11}$ ， $K_{33}=-{\sigma }_{1111}/{\mu }_{11}$ 和 $K_{34}={\tilde{\alpha }}_1{S}_1$ 为简化后的方程系数。

3. 全局分岔

为了便于讨论，将(3)改写为

$\ddot{f}+{\lambda }_{1}f+{\lambda }_2{f}^2+{\lambda }_3{f}^3=\epsilon \left(g\cos \omega t-\delta \stackrel{\dot{}}{f}\right)$ (4)

其中， $f={\eta }_1$ ， ${\lambda }_1={\upsilon }_1^2$ ， ${\lambda }_2=-K_{32}$ ， ${\lambda }_3=-K_{33}$ ， $g=\frac{K_{34}}{\epsilon }$ ， $\delta =\frac{-K_{31}}{\epsilon }$ ， $\omega ={\Omega }_1$ 。研究表明，在部分充液条件下，当充液深度相对较小时，上述各参数均为大于零的量 [16] 。

在方程(4)中 $\epsilon$ 为一小量，其对应的无扰动力系统为

$\{\begin{array}{c}y=\stackrel{\dot{}}{f}\\ \stackrel{\dot{}}{y}=-{\lambda }_{1}f-{\lambda }_2{f}^2-{\lambda }_3{f}^3\end{array}$ (5)

其对时间的一次积分为

$H=\frac{1}{2}{\stackrel{\dot{}}{f}}^2+\frac{1}{2}{\lambda }_1{f}^2+\frac{1}{3}{\lambda }_2{f}^3+\frac{1}{4}{\lambda }_3{f}^4$ (6)

对于不同的H值，它代表相平面上不同的轨道。

系统(5)有3个不动点，它们分别为： $O\left(0,0\right)$ ， $A\left(\frac{-{\lambda }_2}{2{\lambda }_3}-\frac{\sqrt{{\lambda }_2^2-4{\lambda }_1{\lambda }_3}}{2{\lambda }_3},0\right)$ 和 $B\left(\frac{-{\lambda }_2}{2{\lambda }_3}+\frac{\sqrt{{\lambda }_2^2-4{\lambda }_1{\lambda }_3}}{2{\lambda }_3},0\right)$ ，不难判定O、A点是中心，B点是鞍点 [17] 。

现在来讨论无扰动动力系统的同宿运动，当 $H=0$ ，由(6)式，得到

${\stackrel{\dot{}}{f}}^{*}=\pm \sqrt{-{\lambda }_1{f}^{*2}-\frac{2}{3}{\lambda }_2{f}^{*3}-\frac{1}{2}{\lambda }_3{f}^{*4}}$ (7)

积分即得同宿轨道为：

$-\frac{1}{\sqrt{-{\lambda }_1}}\ln \left|\frac{\sqrt{-{\lambda }_1-\frac{2}{3}{\lambda }_2{f}^{*}-\frac{1}{2}{\lambda }_3{f}^{*2}}+\sqrt{-{\lambda }_1}}{f^{*}}-\frac{{\lambda }_2}{3\sqrt{-{\lambda }_1}}\right|=C_1\pm t$ (8)

或

$\left|\frac{\sqrt{-{\lambda }_1-\frac{2}{3}{\lambda }_2{f}^{*}-\frac{1}{2}{\lambda }_3{f}^{*2}}+\sqrt{-{\lambda }_1}}{f^{*}}-\frac{{\lambda }_2}{3\sqrt{-{\lambda }_1}}\right|=C_2{\text{e}}^{\pm \sqrt{-{\lambda }_{1}t}}$ (9)

其中 $C_1$ 和 $C_2$ 为积分常量。

当 $H<0$ 时，得绕O、A两点的两族闭轨道，代表在稳定平衡位置附近的振动，当H增大时，振动的周期增大，当 $H>0$ 时，代表另一类型的周期性振动。

对动力系统(4)，定义Melnikov函数为：

$M\left(\tau \right)={\displaystyle {\int }_{-\infty }^{+\infty }{\stackrel{\dot{}}{f}}^{*}\varphi \left[f^{*},{\stackrel{\dot{}}{f}}^{*},f\omega \left(t+\tau \right)\right]\text{d}t}$ (10)

其中， $\varphi =\varphi \left(f,\stackrel{\dot{}}{f},\omega t\right)=g\cos \omega t-\delta \stackrel{\dot{}}{f}$ 。

故

$M\left(\tau \right)=-\delta {\displaystyle {\int }_{-\infty }^{+\infty }{\stackrel{\dot{}}{f}}^{*2}}\text{d}t+g{\displaystyle {\int }_{-\infty }^{+\infty }{\stackrel{\dot{}}{f}}^{*}\cos \omega \left(t+\tau \right)\text{d}t}$ (11)

式(11)中的第一个积分为

${\displaystyle {\int }_{-\infty }^{+\infty }{\stackrel{\dot{}}{f}}^{*2}}\text{d}t=\frac{4}{3{\lambda }_3}\left[{\left(-{\lambda }_1\right)}^{\frac{3}{2}}+\frac{{\lambda }_2^2\sqrt{-{\lambda }_1}}{3{\lambda }_3}+{\lambda }_2\frac{2{\lambda }_2^2-9{\lambda }_1{\lambda }_3}{9{\lambda }_3\sqrt{2{\lambda }_3}}\arcsin \frac{\sqrt{2}{\lambda }_2}{\sqrt{2{\lambda }_2^2-9{\lambda }_1{\lambda }_3}}\right]$ (12)

考虑式(11)中的第二个积分

$\begin{array}{c}{\displaystyle {\int }_{-\infty }^{+\infty }{\stackrel{\dot{}}{f}}^{*}\cos \omega \left(t+\tau \right)\text{d}t}={\displaystyle {\int }_{-\infty }^{+\infty }\cos \omega \left(t+\tau \right)\text{d}{f}^{*}}\\ =2C_2\sqrt{{\lambda }_1}\omega {\displaystyle {\int }_{-\infty }^{+\infty }\frac{{\text{e}}^{\sqrt{{\lambda }_1}t}\sin \omega \left(t+\tau \right)}{{\left(C_2{\text{e}}^{\sqrt{{\lambda }_1}t}+\frac{{\lambda }_2}{3\sqrt{{\lambda }_1}}\right)}^2+\frac{1}{2}{\lambda }_3}\text{d}t}\end{array}$ (13)

其中，

${\displaystyle {\int }_{-\infty }^{+\infty }\frac{{\text{e}}^{\sqrt{{\lambda }_1}t}\sin \omega \left(t+\tau \right)}{{\left(C_2{\text{e}}^{\sqrt{{\lambda }_1}t}+\frac{{\lambda }_2}{3\sqrt{{\lambda }_1}}\right)}^2+\frac{1}{2}{\lambda }_3}\text{d}t}=2\text{π}\cdot \frac{\sin \omega \left(\xi +\tau \right)\left(ch\omega {\psi }_2-ch\omega {\psi }_1\right)}{C_2\sqrt{2{\lambda }_1{\lambda }_3}\left(1-ch\frac{2\text{π}\omega }{\sqrt{{\lambda }_1}}\right)}$ (14)

$\{\begin{array}{l}\xi =\frac{1}{\sqrt{{\lambda }_1}}\ln \left[\frac{1}{3C_2}\sqrt{\frac{2{\lambda }_2^2+9{\lambda }_1{\lambda }_3}{2{\lambda }_1}}\right]\\ {\psi }_1=\frac{1}{\sqrt{{\lambda }_1}}\left[\text{π}+\text{arctg}\left(\frac{3}{{\lambda }_2}\sqrt{\frac{{\lambda }_1{\lambda }_3}{2}}\right)\right]\\ {\psi }_2=\frac{1}{\sqrt{{\lambda }_1}}\left[\text{π}-\text{arctg}\left(\frac{3}{{\lambda }_2}\sqrt{\frac{{\lambda }_1{\lambda }_3}{2}}\right)\right]\end{array}$ (15)

将(12)和(13)代入式(11)中，得到Melnikov函数

$\begin{array}{l}M\left(\tau \right)=-\frac{4\delta }{3{\lambda }_3}\left[{\lambda }_1^{\frac{3}{2}}+\frac{{\lambda }_2^2\sqrt{{\lambda }_1}}{3{\lambda }_3}+{\lambda }_2\cdot \frac{2{\lambda }_2^2+9{\lambda }_1{\lambda }_3}{9{\lambda }_3\sqrt{2{\lambda }_3}}\arcsin \frac{\sqrt{2}{\lambda }_2}{\sqrt{2{\lambda }_2^2+9{\lambda }_1{\lambda }_3}}\right]\\ +\frac{4\text{π}g\omega \sin \omega \left(\xi +\tau \right)\left(ch\omega {\psi }_1-ch\omega {\psi }_2\right)}{\sqrt{2{\lambda }_3}\left(ch\frac{2\text{π}\omega }{\sqrt{{\lambda }_1}}-1\right)}\end{array}$ (16)

Melnikov函数如果有限，则稳定流形和不稳定流形相交，Poincare映射有马蹄，从而存在着奇怪的不变集合，对我们研究的耗散系统，就有可能出现混沌，由 $M\left(\tau \right)=0$ 得到

$\sin \omega \left(\xi +\tau \right)=\frac{\delta \left(ch\frac{2\text{π}\omega }{\sqrt{{\lambda }_1}}-1\right)}{6\sqrt{2}g\omega \sqrt{{\lambda }_3}\left(ch\omega {\psi }_1-ch\omega {\psi }_2\right)}\left[{\lambda }_1^{\frac{3}{2}}+\frac{{\lambda }_2^2\sqrt{{\lambda }_1}}{3{\lambda }_3}+{\lambda }_2\cdot \frac{2{\lambda }_2^2+9{\lambda }_1{\lambda }_3}{9{\lambda }_3\sqrt{2{\lambda }_3}}\arcsin \frac{\sqrt{2}{\lambda }_2}{\sqrt{2{\lambda }_2^2+9{\lambda }_1{\lambda }_3}}\right]$ (17)

为使方程(17)有根，必须有 $\sin \omega \left(\xi +\tau \right)\le 1$ ，故此，发生混沌的临界条件为

$\frac{g}{\delta }\ge \frac{\left(ch\frac{2\text{π}\omega }{\sqrt{{\lambda }_1}}-1\right)}{6\sqrt{2}\omega \sqrt{{\lambda }_3}\left(ch\omega {\psi }_1-ch\omega {\psi }_2\right)}\left[{\lambda }_1^{\frac{3}{2}}+\frac{{\lambda }_2^2\sqrt{{\lambda }_1}}{3{\lambda }_3}+{\lambda }_2\cdot \frac{2{\lambda }_2^2+9{\lambda }_1{\lambda }_3}{9{\lambda }_3\sqrt{2{\lambda }_3}}\arcsin \frac{\sqrt{2}{\lambda }_2}{\sqrt{2{\lambda }_2^2+9{\lambda }_1{\lambda }_3}}\right]=R\left(\omega \right)$ (18)

4. 算例

对非线性动力方程(4)进行数值计算，取 ${\lambda }_1=40$ ，， ${\lambda }_3=90$ ， $\omega =\text{π}$ ， $\delta =1$ ，得到关于扰动幅值g的分岔图(图3)和 $g=280$ 时，系统的时程曲线，相图和Poincare截面，如图4所示，表明系统在给定参数条件下，扰动幅值增加，系统将发生同宿分岔并进入混沌。

5. 结论

本文通过对液体晃动模态方程进行解析求解，得到了系统的同宿轨道的解析表达式，并利用Melnikov函数，确定了系统在受扰情况下发生混沌运动的阈值，表明系统在受到一定程度的扰动后，会发生同宿分岔和混沌运动。通过本文的结果，可以看到，全局动力学行为在液体晃动动力学分析中的重要性，局部分析结果往往只能对液体小幅晃动的情况进行预测和控制，当充液系统工作在较恶劣工况下，液体非小幅晃动时，系统的动力学行为更值得关注，本文的研究结果表明，当扰动较大时，系统发生混沌运动

Figure 3. Bifurcation diagram

图3. 分岔图

(a) (b)(c)

Figure 4. Curves of the responses at g = 280 (a) Time history; (b) Phase plane; (c) Poincare mapping

图4. g = 280时的各种曲线 (a) 时程曲线；(b) 相图；(c) Poincare 映射

的机理是同宿分岔，液体晃动虽非小幅，但为有限幅，在工程应用中，这一情况更容易得到控制，而不会对系统造成破坏性的影响，因此，此结果可为工程优化设计和控制提供理论依据。

基金项目

国家自然科学基金青年项目(批准号：11502161)、天津市自然科学基金青年项目(合同号：13JCQNJC04300)。