仿生减摇鳍升力模型的数值模拟及分析

李冬松 金鸿章

1哈尔滨工程大学 自动化学院，黑龙江 哈尔滨150001

2哈尔滨工业大学 数学系，黑龙江 哈尔滨150001

仿生减摇鳍升力模型的数值模拟及分析

李冬松1，2金鸿章1

1哈尔滨工程大学 自动化学院，黑龙江 哈尔滨150001

2哈尔滨工业大学 数学系，黑龙江 哈尔滨150001

船舶在系泊或低航速状态下，发动机主机停止工作，船舶失去自控方向能力，船体随波浪左右摇晃，比航行时摇摆更为剧烈。因此，研究零航速减摇鳍非常重要。文章研究零航速时仿生减摇鳍产生升力的模型，其基础为Weis－Fogh机构理论。首先讨论如何旋转才能产生给定的升力。依据吕卡提方程理论，给出模型周期解的存在性和稳定性条件；采用单步Runge－Kutta方法求出模型的数值解，给出保证数值解稳定的条件。

减摇鳍；零航速；吕卡提方程；周期解；稳定性；仿生学

1 引言

常规减摇鳍是目前最常用且应用最成功的船舶主动式减摇装置，减摇效果可达90%以上。然而，只有船舶的航速较高时，减摇鳍才可以有效地减摇，船舶在低航速或零航速情况下，减摇鳍几乎不能进行减摇。这是由其减摇机理决定的，即减摇鳍升力的产生需要高速水流流过鳍表面。当速度很小时，鳍上的升力也变得很小，在零航速时升力也同时消失了。对于在低航速或系泊状态下仍需要减摇的船舶来说，传统的减摇鳍就不再适用了。而减摇水舱减摇效果有限，有时甚至出现增摇现象，这就要求发展零航速减摇鳍技术，此技术在减摇领域一直是个空白。本文依据Weis－Fogh机构理论设计零航速减摇鳍。Weis－Fogh机构是英国的生物学家发明的一种仿生机构，这种机构能在无来流速度的流体中产生升力，这一特性在零航速减摇鳍中得到了应用［1－7］。图1为Weis－Fogh机构的简化模型，它由两个平板翼组成，当Weis－Fogh机构张开时，两个平板翼的翼根点紧紧靠在一起，前缘点B分开，两翼绕A点转动张开，张开角为2απ，角速度为Ω，两翼间形成空隙，迫使周围流体充填该空隙。由于流体不可压缩，当张开角很小时，流体充填空隙的速度非常高，造成两翼上下表面很大的压力差。当Weis－Fogh机构闭合时，两翼对空隙间的流体产生压力，迫使流体流出空隙。与张开时类似，当张开角很小时，流体流出空隙的速度非常高，造成两翼上下表面很大的压力差。当Weis－Fogh机构张开或闭合时，机构运动是对称的，这种对称性可以在一定程度上简化理论分析。

图1 Weis－Fogh机构物理模型

采用Weis－Fogh机构来设计零航速减摇鳍，首先要对它的升力特性进行分析。由于Weis－Fogh机构的对称性，只需对单翼进行升力分析，由对称性即可得到两翼共同作用的升力特性。由伯努利方程求出物面压强p的表达式，将压强沿物面积分，经复杂的推导得出流体对翼的作用力。由文献［8］得，作用于机翼上的力可以表示为角加速度和来流加速度的一次函数的线性组合；角速度和来流速度的二次函数的线性组合，即升力：

式中，u为来流速度；u′为来流加速度；Ω为翼的旋转角速度；Ω′为翼的旋转角加速度。

当讨论零航速时，来流速度u＝0，来流加速度u′＝0。如令T＝f1（t），Ω＝y（t），Ω′＝y′（t），从而式（1）可以简化成：

本文研究翼如何旋转才能产生船舶减摇所需要的指定升力f1（t），此时问题转化为求解下面微分方程：

由文献［7］，m1＜a（t）＜M1＜0。本文仅考虑升力控制，由于翼回摆时产生的反向力涉及到翼的变形，我们将再另文研究。在此简化条件下，模型中a（t）成为周期函数，周期为2π。

2 模型数值解的分析及在减摇鳍控制中的应用

首先介绍模型精确解的性质，当所需升力f（t）是周期为2π的连续函数时，模型为吕卡提方程，吕卡提方程近年来得到广泛的研究［9－10］。

模型的示性代数方程为

由文献［8］，a（t）＜0，升力f（t）＞0，故a（t）f（t）＜0。当a（t），f（t）为周期为2π的连续函数时，依据吕卡提方程性质，已知方程（3）存在两个周期为2π的连续解y1（t），y2（t），且y＝y1（t）为稳定的周期解，y＝y2（t）为不稳定的周期解。

由于常微分方程的精确解难以求出，故需要采用数值方法求方程的近似解，本文采用单步Runge－Kutta方法，它具有精度高，稳定性好的特点。对升力模型：

采用Runge－Kutta方法离散得：

式中，f∈C∞，f（t+2π）=f（t），y0是已知数，A= aij）s×s，b=（b1， …，bs），c=（c1， …，cs） 由具体Runge－Kutta方法确定。

下面首先进行几个数值试验，探讨离散模型的性质。

试验1 y＝－2y2＋（2－sin（t）），y（0）＝5。

使用Heun’s公式得：图2。

图2 试验1中数值解在t=5时呈现周期性

图3 试验2中数值解在t=2时呈现周期性

试验2 y′=-2y2+（2-sin t），y（0）＝0.1。仍使用Heun’s公式，计算结果见图3。

试验3 y′=-10y2+（2-cos（2t）），y（0）＝0.1，h＝0.1。

当方程为y′＝f（t，y）时，Euler公式为yn+1=yn+f（tn，yn），应用于模型得：

计算结果如图4。

试验4 y=-10y2+2-cos（2t），y（0）＝1，仍使用Euler公式，结果如图5。

图4 试验3中数值解在t=1时呈现周期性

图5 试验4中数值解在t=0.4时呈现周期性

从数值试验中可观察出以下结论，当步长较小时，数值解保持精确解的渐近周期性。但在研究方程y′=-10y2+2-cos（2t）时，我们发现，当初值y（0）较大，如y（0）＝2，h＝0.1，数值解不稳定，发散速度极快。这说明离散控制时，还应考虑步长h，下面讨论步长的选取问题。

3 离散模型的稳定性分析

下面研究Runge－Kutta方法求解方程 （5）的渐近周期性，将R－K法应用于方程（5）可得离散模型：

此处，h＝2π／N，N是正整数。

引理1 设a（t）＜0，f（t）＞0，a（t），f（t）∈C，且a（t）＝a（t＋2π），f（t）＝f（t＋2π），t＞0，令K＝max｝和h0如果0＜h＜h0，则存在正数M，m和N，使得m＜yn＜M，m＜Y＜M，i＝1，2，3，4，n＞N。

证明：设初始值y0＞mF＋MF，则（Y）2+ h 2 f（t0）<0，由式（9）得，Y＜y0。 因为h＜h0＜，我们有

同样道理，

设y0>y1>…>yN-1>yN，yN

式（12）与yN-1＞yN矛盾。使用归纳法即可得，若干步后数值解将在mf和K之间振动。定理证毕。

用类似的方法还可证明

定理1 设a（t）＜0，f（t）＞0，a（t），f（t）∈C，a（t）＝a（t＋2π），f（t）＝f（t＋2π），t＞0，

由以上研究可得在计算升力时，必须注意步长的选取，否则将失去稳定性。

4结论

船舶零航速仿生减摇鳍仍处于研究阶段，有大量的理论和实际问题需要解决。本文讨论Weis－Fogh机构如何旋转才能产生船舶减摇所需要的指定升力及稳定性问题，首先论证了升力模型稳定性，然后给出保证数值解稳定的步长选取条件，为仿生减摇鳍的控制打下理论基础。实际上关于Weis－Fogh机构的理论还有很多工作要做，如可变翼的升力模型的计算，及回摆时力的计算，这些工作将在以后研究。

［1］王宇，金鸿章.船舶零航速减摇鳍升力测量研究［J］.传感器与微系统，2007，26（1）：11－14.

［2］JIN HZ，WANG Y，QI Z G.Study on lift generation of Weis－Fogh flapped fin stabilizer at zero speed［C］／／Korea，Busan：SICE－ICASE International Joint Conference，2006.

［3］JIN H Z，QI Z G，WANG Y.Research on ship roll stabilization at zero speed［C］／／Korea，Busan：SICEICASE International Joint Conference，2006.

［4］DANIEL T，COMBES S A.Flexible wings and fins：bending by inertial or fluid dynamic forces［J］.Integrative Comparative Biology，2002，42（5）：1044－1049.

［5］SANE S P，DICKINSON M H.The aerodynamic effect s of wing rotation and a revised quasi2steady model of flapping flight［J］.Journal of Experimental Biology，2002，205：1087－1096.

［6］王宇，金鸿章，綦志刚，等.船舶零航速减Weis－Fogh减摇鳍升力仿真研究 ［J］.海军工程大学学报，2007，19（3）：21－25

［7］SOHN M H，CHANG J W.Flow visualization and aerodynamic load calculation of three types of clap－fling motions in a weis－fogh mechanism ［J］.Aerospace Science and Technology，2007，11（2-3）：119-129.

［8］章社生，王献孚，吴秀恒.Weis－Fogh机构流体力学［M］.北京：国防工业出版社，1999.

［9］ BALASUBRAMANIAM P，SAMATH J A，KUMARESAN N，et al.Solution of matrix Riccati differential equation for the linear quadratic，singular system using neural networks［J］.Applied Mathematics and Compution，2006，182（2）1832-1839.

［10］ZHU J，LI K.An iterative method for solving stochastic Riccati differential equations for the stochastic LQR problem［J］.Optimization Methods and Software，2003，18（6），721-732.

The Analysis and Numerical Simulation for the Lift Model of Fin Stabilizer

Li Dong－song1，2Jin Hong－zhang1
1 College of Automation，Harbin Engineering Univ.，Harbin 150001，China 2 Dept.of Mathematics，Harbin Institute of Technology，Harbin 150001，China

The engine of ships is laid off at anchor.Ships drift with wave and loss the capacity of controlling navigating direction for self，so the roll is increased and severer than shipping state.It is important to study stabilizer when a ship is at zero speed.In this paper we study the lift model of fin stabilizer that based on potential theory of the Weis－Fogh mechanism when a ship is at zero speed.Firstly，we discussed that Weis－Fogh mechanism how to rotate to generate the lift that has been given.According to the theory of Riccati differential equation，conditions that assure the existence and stability of period solution are given.Then，using Runge－Kutta methods，the numerical solution of the model is obtained.Finally we give the conditions that assure the numerical solution is stable.The results are closely matched with numerical experiments.

antirolling fin；zero speed；Riccati differential equation；periodic solution；stability；bionics

TP212

：A

：1673－3185（2009）01－29－04

2008-10-1 2

国家自然科学基金（50575048）；黑龙江省博士后资助项目（LBH-205052）

李冬松（1963-），男，副教授，博士。研究方向：延迟微分方程精确解和数值解的定性性质。

E-mail：lds_ch@yahoo.com.cn

金鸿章（1946-），男，教授，博士生导师。研究方向：控制理论和控制系统、智能控制、船舶控制系统与减摇装置