戴余良，王长湖，苗 海，刘祖源

（1海军工程大学 船舶与动力学院，武汉 430033；2青岛92196部队装备部，山东 青岛 266011；3武汉理工大学 交通学院，武汉 430063）

潜艇水下运动稳定性非线性分析研究

戴余良1，王长湖2，苗 海1，刘祖源3

（1海军工程大学 船舶与动力学院，武汉 430033；2青岛92196部队装备部，山东 青岛 266011；3武汉理工大学 交通学院，武汉 430063）

应用非线性系统运动稳定性理论和同伦延拓数值计算方法，分析了潜艇水下运动的稳态响应及运动稳定性，潜艇应急上浮运动稳定性实例分析结果通过其大攻角六自由度运动方程数值积分的动态仿真进行了验证，表明稳定性分析结果是正确的。从而说明了采用同伦延拓法分析潜艇水下运动稳定性是有效的。

潜艇；运动稳定性；非线性动力学；同伦延拓法

1 引 言

潜艇是主要在水下活动的航行体。为了保证潜艇水下航行安全，特别是在舱室破损进水、舵卡等事故状况下的安全，清楚地了解潜艇的动态响应特性和稳定性是至关重要的。

上述稳定性评价指标基于如下假设：在小扰动基准直线运动中，水平面运动与垂直面运动之间的耦合相对较弱，可以忽略。但是，对于高速航行的潜艇，在航行的极限状态或应急状态下，上述关于水平面与垂直面之间运动不耦合的假设不再成立，此时在所有六个自由度都可能发生大幅运动，各种运动模态之间的非线性相互作用变得非常显著。因此，简单地使用指标Gv和Gh可能导致不正确的结论。此时必须考虑水平面与垂直面运动之间相互耦合的运动特性。

潜艇的水下运动，尤其是应急情况下的操纵运动通常是非线性的。随着非线性理论的日益成熟，运用非线性动力学理论对潜艇的运动微分方程进行分析，能够更准确地对潜艇的运动规律作出描述，从而得到更准确的结论。本文运用非线性系统运动稳定性理论，对潜艇水下运动稳定性的分析方法进行了深入研究。

2 非线性系统运动稳定性

运动稳定性理论自Lyapunov于十九世纪九十年代创建以来，在物理科学和工程技术等各个领域都获得了广泛的应用。

Lyapunov意义下的运动稳定性理论是研究干扰因素对物体运动的影响。众所周知，微小的干扰因素对物体运动的影响，对于不同的运动是不一样的。对于一些运动，这种影响可能并不显著，因而受扰运动与未扰运动相差很小；反之，对于另外一些运动，干扰的影响可能很显著，以至于干扰力无论多么小，受扰运动与未扰运动随着时间的推移而可能相差得很大。简单地说，属于前者的运动称为稳定的，而属于后一类型的运动则称为不稳定的。由于在现实世界中，干扰因素总是不可避免地存在着，所以运动稳定性的问题有很重要的理论意义和实际意义，这也正是稳定性理论蓬勃发展的重要原因。

2.1 运动稳定性判定方法

在研究工程问题时，系统的渐近稳定性具有特别重要的意义。那么如何判断所考察的系统是稳定的还是不稳定的？Lyapunov直接法和Lyapunov一次近似法是两种最基本的分析方法。由于Lyapunov直接法的成功与否取决于构造合适的Lyapunov函数。Lyapunov函数的构造需要一定的经验和技巧，目前尚无一般规律。因此，Lyapunov一次近似法被广泛应用于非线性系统的稳定性分析。

Lyapunov一次近似法的基本思想[2-3]是：将非线性系统线性化，仅保留一次项，研究该线性系统的稳定性。（1）如果一次线性近似系统的所有特征值都具有负实部，则原非线性系统的零解（平衡点）是渐近稳定的；（2）如果一次线性近似系统的特征值中只要有一个具有正实部，则原非线性系统的零解是不稳定的；（3）如果一次线性近似系统存在实部为零的特征值，而其余所有特征值都具有负实部，此时非线性系统零解的稳定性取决于系统的非线性项，不能用Lyapunov一次近似法来研究，它的判定方法十分复杂和困难，本文试图利用奇异性和分叉理论来分析。

2.2 运动稳定度

物体的运动稳定性通常均可通过分析系统方程的特征值λ的性质来确定，当特征值的实部Re(λ)＜0时，物体的扰动运动逐渐衰减，且Re(λ)越负，衰减越快，系统越稳定；反之，当Re(λ )＞0时，物体的扰动运动越来越强，且Re(λ)越大，系统越不稳定；如果存在一对共轭复特征值，则物体运动会发生振荡。因此，系统特征值λ实部的正负决定了系统的运动稳定性，而λ实部的大小决定了物体扰动运动衰减的快慢，如果λ存在一对共轭复数，则物体扰动运动周期性振荡。

本文定义系统所有特征值的最大实部为系统的稳定度 （degree of stability），记为e，即e=max Re λi()｛｝i=1,2,…,n，其中λi为系统的特征值，则e越负系统越稳定，e越正系统越不稳定。

3 同伦延拓非线性数值分析方法

同伦延拓方法是求解非线性问题的一种非常有效的方法，它克服了传统迭代法局部收敛的弱点，对初值的选取没有严格限制，能够全局收敛，可以求出问题的全部解。它虽然于二十世纪70年代才开始发展起来，但在短短的30多年里，同伦延拓方法已经获得了广泛应用[4-9]。传统的求解方法如迭代法、消去法等存在的最大问题是依赖于初值的选取，初值选择不当常常导致迭代过程不收敛，而且很难求出全部解。因此，同伦延拓方法被誉为二十世纪数学研究中一项突破性的新成果[4]。

3.1 同伦延拓法的基本思想

同伦延拓法的基本思想[4-9]是借助同伦函数H，从辅助映射g的零点集出发，跟踪同伦方程组H(x, λ )=0的解曲线，走到目标映射f的零点集。

为了求解非线性方程组f(x)=0，其中光滑映射 f：Rn→Rn，x=｛x1,x2,…,xn｝T，选取具有已知零点的光滑辅助映射 g：Rn→Rn，将其与映射f构成一线性同伦函数H：Rn×Rn→Rn，使得H( x, 0 )=g(x)和H( x, 1 )=f(x)，典型的同伦函数为：

对于同伦方程组H( x, λ )=0，当λ=0时，其解与辅助方程组g(x)=0的解相同；让同伦参数λ从0逐渐变化到1，并跟踪同伦方程组的解，则当λ=1时，同伦方程组的解便是待解方程组f(x)=0的解。

3.2 同伦延拓算法

根据同伦延拓法的基本思想，利用同伦延拓法求解非线性方程组f(x)=0首先应该选择合适的辅助函数g(x)。g(x)的选取方法很多，本文选择g(x)=f(x)-f( x( 0 ))。那么同伦函数为：

因此，H( x, 0 )=f(x)-f( x( 0 ))和H( x,1)=f(x)，其中x(0)是系统f(x)=0对应于λ=0时的解，x(1)是系统f(x)=0对应于λ=1时的解，x(λ)是同伦方程组H( x( λ ),λ )=0的唯一解。

利用同伦延拓法，非线性方程组f(x)=0的求解问题转化成了求解微分方程：

其中x(0)是系统f(x)=0的初始已知解。

微分方程（1）的求解可以采用任一种数值计算方法，本文采用4阶Runge-Kutta法求解，则4阶Runge-Kutta法的迭代公式为：

图1 同伦延拓算法流程图Fig.1 Flow of homotopy and continuation algorithm

式中：h为计算步长。同伦参数λ从0逐渐变化到1，则当λ=1时的解x(1)即是系统f(x)=0的近似解。

为了避免由于Jacobi矩阵J(x)不可逆时上述Runge-Kutta算法失效，将计算ki的等式（2）改写成求解如下线性方程：

具体算法流程见图1。

4 潜艇水下运动非线性方程稳态解

4.1 潜艇水下运动稳态方程

为了研究潜艇运动的非线性特性，必须采用大攻角六自由度非线性运动方程，本文采用文献[10-11]的潜艇运动方程，可化成一般形式：

式中状态变量 x=[u,v,w,p,q,r,φ,θ,ψ ]，控制变量 u= [δr,δb,δs,δB,xGB,yGB,zGB,n ]，其中 δB=B-W，xGB=xG-xB，yGB=yG-yB，zGB=zG-zB，n 为螺旋桨转速。

4.1.1 稳态条件

潜艇水下运动的稳态条件是：处于定常运动状态（定常直航或定常回转状态），即

4.1.2 稳态运动方程

将稳态条件（6）式代入方程（5）便可得到潜艇水下稳态运动方程：

4.2 潜艇水下运动方程稳态解

由于稳态方程（7）是高度非线性的，其解可能出现分叉和多解现象，因此本文使用同伦延拓法求解稳态方程（7）的可能解。

潜艇在水下运动过程中，潜浮速度Vζ、横倾角φ、纵倾角θ和攻角α是影响潜艇安全的主要参数，所以重点分析稳态的潜浮速度Vζ、横倾角φ、纵倾角θ和攻角α随控制参数u= [δr,δb,δs,δB,xGB,yGB,zGB,n ]的变化情况。

潜艇水下运动稳态方程求解及运动稳定性分析过程：

（1）通过动态方程（5）采用数值积分法求得某一开始稳定状态的状态变量值，以此作为同伦延拓算法的初值；

（2）以某一控制参数作为同伦参数，其它控制参数的值保持不变，计算出各状态变量随同伦参数变化的稳态解曲线及其Jacobi矩阵的特征值；

（3）改变除同伦参数以外的其它控制参数（作为扰动参数）的值，重复上面的计算，得出另一组稳态解曲线。

（4）根据特征值的性质确定平衡解的稳定性。

4.3 潜艇水下运动稳定性

潜艇水下运动方程可写成如下形式：

式中状态变量 x=[u,v,w,p,q,r,φ,θ,ψ ]。

为了分析潜艇水下运动在某个平衡状态x0的稳定性，根据Lyapunov线性化原理，将水下运动方程（8）在平衡点 x0附近线性化，即：

其中A为f(x)在x0处的Jacobi矩阵，即

由运动稳定性理论可知，若矩阵A的所有特征值有负实部，则平衡状态x0是渐近稳定的；只要矩阵A有一个特征值有正实部，则平衡状态x0是不稳定的。因此，当系统的稳定度e＜0时，则平衡状态x0是渐近稳定的；当系统的稳定度e＞0时，则平衡状态x0是不稳定的。

5 潜艇水下运动稳定性分析实例

应急上浮是潜艇在舱室破损进水、舵卡等事故状况下所采取的紧急操纵措施，是潜艇非战时紧急脱险的最有效方法。为了确保采取的操纵方法正确有效，保证潜艇应急上浮安全，清楚地了解潜艇在应急状况下的动态响应特性和稳定性至关重要。潜艇应急上浮通常是一种大攻角强机动，描述潜艇的运动模型是强非线性的。下面针对潜艇应急上浮运动，以xGB为同伦参数分析潜艇运动稳定性为例，说明本文提出的潜艇水下运动稳定性非线性分析方法。潜艇模型参数及水动力系数采用文献[10-11]提供的数据，其中艇重53.4kN，艇长5.3m。

5.1 稳态解及特征值

计算条件：δr=δs=δb=0，n=1 000rpm，yGB=0，zGB=0.061m，以 δB为扰动参数，分别取 δB＝1%W，5%W 和10%W，其中W为艇重，xGB的变化范围为：-10%L～10%L，其中L为艇长。

计算结果如图2～10所示。

5.2 结果与讨论

由图2～10可以看出：

图2 对于不同δB航速U随xGB的变化 Fig.2 Speed versus xGBfor variations in δB

图3 对于不同δB横倾角φ随xGB的变化Fig.3 Roll angle versus xGBfor variations in δB

图4 对于不同δB纵倾角θ随xGB的变化Fig.4 Pitch angle versus xGBfor variations in δB

图5 对于不同δB攻角α随xGB的变化Fig.5 Attack angle versus xGBfor variations in δB

图6 对于不同δB深度速率ζ˙随xGB的变化 Fig.6 Speed depth versus xGBfor variations in δB

图7 对于不同δB稳定度随xGB的变化Fig.7 Degree of stability versus xGBfor variations in δB

图8 δB=1%W时系统特征值随xGB的变化Fig.8 Eigenvalue versus xGBfor δB=1%W

图9 δB=5%W时系统特征值随xGB的变化Fig.9 Eigenvalue versus xGBfor δB=5%W

图10 δB=10%W时系统特征值随xGB的变化Fig.10 Eigenvalue versus xGBfor δB=10%W

（1）潜艇的航速、上浮速度和纵、横倾角和攻角等的稳态响应在正浮力的作用下，在xGB=2%L附近均有一个突变值，由此看到了当δr=0时（即不操方向舵），潜艇在正浮力作用下出现了突倾现象（见图3），即出现了平面外运动的情况。在xGB的其它取值范围内潜艇运动的变化相对平稳。

（2）由图7可知，当xGB=2%L左右时潜艇在正浮力作用下上浮速度、纵、横倾发生突跳现象，主要是由于此时系统的稳定度是一个较大的正数值。同时由图8～10可以看出，系统特征值中出现了幅值较大的共轭虚数，此时运动将发生大幅振荡。

（3）图7系统稳定度表明，潜艇在正浮力作用下当xGB≤0时，稳定度均为负值，其运动是稳定的，同时由图8～10可以看出，系统特征值存在一对幅值较小的共轭复数，此时运动将发生轻微振荡。xGB＞0时稳定度基本上为正值，其运动是不稳定的。

因此，当xGB≤0时潜艇在正浮力作用下，稳定度均为负值，其运动是稳定的，当xGB＞0时，稳定度过零点，由负值变成正值，运动由稳定变为不稳定，出现分叉点，同时由图8～10可以看出，在xGB=0前后系统特征值均为实数，此时为静态分叉。

6 数值仿真

为了验证前面稳定性分析实例结果的正确性，本节采用稳态方程（7）所对应的文献[10-11]提供的潜艇大攻角六自由度运动方程，潜艇模型参数及水动力系数与前面稳定性分析实例相同，采用文献[10-11]提供的数据，其中艇重53.4kN，艇长5.3m。利用Matlab通过龙格库塔（Runge－Kutta）数值积分进行了仿真，将仿真结果与前述分析结果进行对比。

6.1 仿真结果

仿真条件：δr=δs=δb=0，n=1 000rpm，yGB=0，zGB=0.061m，δB＝1%W， 其中 W 为艇重，xGB=-5%L 和 2%L，其中L为艇长。

仿真结果如图11所示。

6.2 结果讨论

由图11的仿真结果可以看出：

（1）当xGB=-5%L时潜艇纵倾角θ动态响应的稳态解，由图7可知该稳态解是稳定的。

（2）当xGB=2%L时潜艇纵倾角θ动态响应的稳态解，由图7、8可知该稳态解是不稳定的，此时潜艇的纵倾响应是一种周期性的振荡运动。

因此，由上面的分析可知，通过数值积分前面稳定性分析结果的正确性得到了初步验证，它与数值仿真的结果完全一致，说明了前述稳定性分析结果是正确的。

图11 xGB=-5%L和2%L，δB=1%W时潜艇纵倾的动态响应Fig.11 Pitch response for xGB=-5%L and 2%L， δB=1%W

7 结 论

本文对非线性系统运动稳定性理论及其数值分析方法进行了研究，提出了将同伦延拓非线性数值计算用于潜艇水下运动的稳态响应及运动稳定性分析的新方法，潜艇应急上浮运动稳定性实例分析结果通过其大攻角六自由度运动方程数值积分的动态仿真进行了验证，表明稳定性分析结果是正确的。从而说明了本文提出的潜艇水下运动稳定性非线性分析方法是正确有效的。

[1]施生达.潜艇操纵性[M].北京:国防工业出版社,1995.

[2]阮 炯,顾凡及,蔡志杰.神经动力学模型方法和应用[M].北京:科学出版社,2002.

[3]陆启韶.常微分方程的定性方法和分叉[M].北京:北京航空航天大学出版社,1989.

[4]王则柯,高堂安.同伦方法引论[M].重庆:重庆出版社,1990.

[5]Lee Jaewook,Chiang Hsiao-Dong.A singular fixed-point homotopy method to locate the closest unstable equilibrium point for transient stability region estimate[J].IEEE Transactions Circuits and Systems-II:Express Briefs,2004,51(4):185-189.

[6]Burden R L,Douglas Faires J.Numerical Analysis[M].Seventh Edition.Thomson Learning,2000.

[7]陈 永,李 立等.一种基于同伦函数的迭代法─同伦迭代法[J].数值计算与计算机应用,2001(2):149-155.

[8]周 洲,祝小平,张肇炽等.一种新的同伦解曲线跟踪算法[J].西北工业大学学报,1998,16(2):246-250.

[9]黄象鼎,曾钟钢,马亚南.非线性数值分析[M].武汉:武汉大学出版社,2000.

[10]Ibrahim Aydin.Out of plane solutions of submarines in free positive buoyant ascent[R].AD-A271544,1993.

[11]Healey A J,Lienard D.Multivariable sliding mode control for autonomous diving and steering of unmanned underwater vehicles[J].IEEE Journal of Ocean Engineering,1993,18(3):327-339.

Study on nonlinear analysis of motion stability of submarines under water

DAI Yu-liang1,WANG Chang-hu2,MIAO Hai1,LIU Zu-yuan3
(1 Naval Architecture and Power College,Naval University of Engineering,Wuhan 430033,China;2 Equipment Department,Army 92196,Qingdao 266011,China;3 Transportation College,Wuhan University of Technology,Wuhan 430063,China)

The motion stability theory of a nonlinear system and the principles of continuation based on homotopy techniques were applied to stable state response and motion stability analysis of submarines during emergency ascent.The stability analysis results were verified by simulations using direct numerical integrations of the fully nonlinear,coupled six degrees-of-freedom equations of motion for the submarine.It indicates that homotopy and continuation theory techniques could be an effective analysis method for motion stability of submarines under water.

submarine;motion stability;nonlinear dynamics;homotopy and continuation theory techniques

U661.3

A

1007-7294（2011）08-0844-09

2010-10-11 修改日期：2011-03-28

国防预研项目资助（1010501020301）

戴余良（1966-），男，博士，副教授，主要从事船舶操纵运动建模、智能控制与实时仿真研究，E-mail:yuliang_dai@163.com。