基于横摇衰减数据的横摇阻尼估算方法对比研究

章 东，王文涛，卜淑霞，刘 伟

（1.中国船舶科学研究中心，江苏无锡 214082；2.水动力学重点实验室，江苏无锡 214082）

0 引 言

横摇是船舶在波浪中6个自由度运动中最难理解和显示的复杂现象[1]，其运动阻尼是船舶耐波性与波浪稳性预测精度的决定性因素。一般来说，船舶运动阻尼主要由船体辐射出的重力波所致，然而相比于其他自由度，船舶横摇受流体粘性影响大。横摇阻尼是由波浪、旋涡、升力、摩擦等因素共同产生，这导致基于势流理论的边界元方法（BEM）预测横摇运动精度过低。针对此难题，ITTC（国际拖曳水池会议）的波浪稳性委员会与耐波性委员会开展了大量工作[2-5]。目前横摇阻尼估算方法主要包括三类：模型试验、CFD计算、经验公式。中国船舶科学研究中心顾民等[6-8]将模型试验与CFD计算结合，对一系列标准模型做了深入研究。模型试验与CFD 计算结果都需要通过数据处理以获得阻尼值，其中自由衰减试验具有试验难度小、成本低的特点，该方法被广泛采用。

从横摇衰减曲线估算横摇阻尼的方法众多，自1872 年Froude[9]的工作开始，人们在该领域做了大量工作。目前除经典的线性方法（消灭曲线法）、Froude能量法[1]外，意大利的Bulian 等[10]、德国的Wassermann 等[11]在对数衰减法方面做了深入研究，英国的Roberts[12]提出的能量方法、国际海事组织推荐的最小二乘法[13]、上海交通大学曾智华[14]的粒子群方法、香港城市大学Chan 等[15]的渐近方法、韩国的Kim 等[16]的Hilbert 变换方法、中国船舶及海洋工程设计研究院封培元等[17]的改进型消灭曲线法、中国海洋大学孙金伟等[18]的Prony-SS 方法、随机减量法[19]等方法均成功用于横摇阻尼估算。其中随机减量法是从波浪中船舶横摇运动数据中提取横摇衰减运动的方法，该方法能从非规则波和规则波船舶横摇运动中提取横摇衰减曲线。

船舶横摇运动特别是大角度横摇运动非线性强，基于此特点，从横摇衰减曲线估算横摇阻尼的方法主要分为两类：一是基于分段线性假设，将横摇衰减曲线分成多段，分别处理各小段拟线性运动的阻尼关系；二是从横摇运动微分方程出发，先假设阻尼模型及复原力臂模型（未知时），后通过渐近方法、数值方法等手段经计算机迭代识别出阻尼模型参数，该类方法称为参数识别技术（parameter identification technique，PIT）。本文的目的是针对这两类方法进行分析，以获取高精度阻尼预报方法。

1 横摇运动数学模型

横摇运动时，粘性项与波浪项产生的阻尼量级相当，直接求解船体周围流体速度场以预报船舶横摇运动难度大且精度低，现阶段一般通过弹簧阻尼滑块模型来模化横摇运动。解耦并归一化单自由度横摇运动方程[3]如下：

表1 横摇阻尼估算方法分类Tab.1 Classification of roll damping estimation methods

2 基于分段线性假设的横摇阻尼估算方法

该类方法的基本思想是分段研究每个周期内横摇，假设每单个周期内的运动均为正弦（或余弦）运动，再按如图1所示流程估算横摇阻尼。除消灭曲线法与改进型消灭曲线法外，基于分段线性假设的其他方法，都要引入一个重要概念：等效线性阻尼系数Beq及其归一化形式μeq，某一时间段内式（1）中的阻尼项d(φ̇ )表示如下[13]：

图1 基于分段线性假设的横摇阻尼估算流程图Fig.1 Procedure of roll damping estimation based on piecewise linear assumption

式中，μeq为该时间段内等效阻尼系数，φa为该周期内横摇角幅值，ωE为该周期内横摇圆频率。获得等效阻尼μeq后，由能量关系可推导出非线性阻尼模型（如式（3））系数与等效阻尼的关系[10]：

2.1 消灭曲线法

假定横摇复原力臂——GZ( )φ为线性的情况下，横摇衰减运动由下式控制：

需要在时历曲线的三个不同位置分别取值，获得三个三元一次方程以求解式（15）中的三个非线性阻尼模型系数b1、b2、b3。

2.2 改进型消灭曲线法

封培元等[17]提出了复原力臂——GZ( )φ为3阶模型、阻尼模型为2项情况下的改进型消灭曲线法：

式（19）即扩展的改进型消灭曲线关系，其系数ε1、ε2、b1、b2、b3均可为0，与文献[17]相比能适应更多复原力臂模型与阻尼模型。

2.3 对数衰减法

对数衰减法（Bulian）[10]将相邻峰值时间段[ti+1-ti]内横摇衰减运动类比为二阶线性系统振动过程：

根据线性系统振动理论易得

若式（25）～（26）中无阻尼自然横摇频率ωx未知，可利用式（27）代入横稳性高-- --GM与船宽BWL估算，

2.4 Hilbert变换法

Hilbert变换相当于90°移相[16]，记横摇衰减信号φ(t)的Hilbert变换为φ̂(t)，则可以构造复平面内解析信号:

式中，μeq(t)为随时间变化的等效阻尼项，ω2x(t)为随时间变化的横摇自然频率项。将式（28）代入式（30）可得以下关系：

根据式（31）与式（32），只要从横摇衰减曲线中拟合出A(t)与ωE( )t，即可获得等效阻尼随时间的变化关系，进而得到每个横摇周期的等效阻尼μeq。

2.5 Froude能量法

根据系统做功能量关系，相邻横摇幅值间阻尼做功应等于系统势能的变化量。当从横摇幅值φi运动到φi+2时,系统势能变化量与阻尼耗散的能量分别为

令式（33）与式（34）相等，可得

归一化的等效阻尼系数μeq的估算同式（26）。

2.6 Roberts能量法

系统机械能EC的耗散率应等于阻尼做功的功率。

横摇运动到幅值φi时，角速度为0，式（36）右侧第一项为0。

根据式（37）可求得幅值φi处的EC值，拟合EC与时间关系曲线，再求导并结合式（22）可得

归一化的等效阻尼系数μeq的估算同式（26）。

3 基于参数识别的横摇阻尼估算方法

该类方法需预设未知项的模型，一般根据经验选择合适的阻尼和复原力臂模型，得到有待定参数的二阶常微分方程，再通过数学手段寻找最接近实验数据的解的模型参数。

3.1 最小二乘法与粒子群方法

以采用式（9）所示阻尼模型为例。

步骤一：先猜测一组初始参数如b1、b2、b3。

步骤二：利用数值方法(如Runge-Kutta法)或解析方法求解横摇运动方程，生成解曲线。

步骤三：定义优化目标函数，如式（39），利用最优化方法（非线性最小二乘优化中的Levenberg-Marquardt优化、粒子群优化[14]等方法均被成功应用于横摇阻尼识别）寻找使得该式最小的一组参数：

式中，φi为计算所得横摇角数据，φexp,i为试验横摇角数据，既可采用整条衰减曲线数据也可只选用幅值数据。为减少计算量，一般采用幅值数据。

3.2 Prony-SS方法与渐近方法

Prony-SS方法最先应用于信号处理领域，对于给定横摇衰减曲线，利用Prony-SS方法重构出衰减曲线的近似表达式[18]，后通过矩阵运算获得阻尼参数值。Prony-SS 方法能提取出横摇信号中的主要成分，具有内在的噪声抑制机制。

渐近方法是求解非线性微分方程的有力手段，对于利用衰减曲线估算横摇阻尼，该方法[15]先利用摄动方法求解出横摇衰减运动方程的渐近解表达式，再通过曲线拟合获得表达式系数，利用消灭曲线反复迭代优化阻尼系数值。

4 横摇阻尼估算-研究结果与讨论

4.1 研究对象及数据简介

为验证并比较上述方法的特点，本文采用标模DTC[20]的横摇衰减试验数据估算横摇阻尼，DTC 标模是德国杜伊斯堡埃森大学按照超巴拿马型集装箱船体以1∶59.407 的比例设计建造的，船体剖面如图2 所示。船模数据如表2 所示，采用压载（ballast）、满载（full condition）两种工况下的无航速自由横摇衰减数据（如图3 所示）作为样本，其复原力臂曲线如图4 所示，从图中可以看出，在试验横摇数据范围内，复原力臂曲线拟合效果良好。采用Wassermann[11]文中推荐的Froude 能量法复现其阻尼估算结果后，利用本文中所述方法处理该数据并进行结果对比。

图2 DTC标模剖面图[20]Fig.2 Hull sections of DTC[20]

图3 横摇衰减曲线Fig.3 Roll decay curve[20]

图4 采用5阶级数拟合的——GZ曲线[11]Fig.4 Fitting curve of——GZ based on 5th oder series[11]

表2 DTC标模参数[20]Tab.2 Calculation parameters of the DTC model[20]

4.2 基于分段线性假设的横摇阻尼估算结果

图5（a）为Wassermann[11]对图3（a）所示的横摇衰减曲线先进行滤波处理后再基于Froude能量法获得的无量纲阻尼系数结果，图5（b）为作者对图3（a）所示的横摇衰减曲线直接使用Froude 能量法后获得的无量纲阻尼系数结果，二者的大振幅结果基本一致，小振幅结果振荡大，这与Wassermann 对数据进行了滤波[11]操作有关。图中无量纲阻尼系数由下式得到：

图5 不同横摇幅值下的无量纲阻尼（压载（D=0.2018 m），无航速）Fig.5 Non-dimensional roll damping with different amplitudes,ballast condition(D=0.2018 m)，without forward speed

按照图1所示流程分别处理图3所示的横摇衰减曲线，阻尼模型采用式（9）。表3及图6～9展示了压载工况下各方法（为方便比较对方法进行了编号）的结果，表4 展示了满载工况下非线性阻尼系数的拟合结果。图6～8中数据点振荡均较大，这与图3（a）的数据尤其是相邻正负幅值的振荡密切相关，说明实验数据质量对基于该类方法的阻尼系数估算有影响。

图6 采用式（15）拟合的消灭曲线与采用式（19）拟合的改进型消灭曲线Fig.6 Fitting curve of extinction curve method based on Eq.(15)and enhanced extinction curve method based on Eq.(19),ballast condition(D=0.2010 m),without forward speed

图7 基于式（7）的等效阻尼拟合，压载（D=0.2018 m），无航速Fig.7 Curve fitting of equivalent linear damping coefficients based on Eq.(7),ballast condition(D=0.2018 m),without forward speed

图8 基于式（7）的等效阻尼拟合，压载（D=0.2018 m），无航速Fig.8 Fitting curve of equivalent linear damping coefficients based on Eq.(7),ballast condition(D=0.2018 m),without forward speed

图9 幅值与频率随时间的关系曲线拟合-Hilbert变换方法，压载（D=0.2018 m），无航速Fig.9 Fitting curve of amplitudes and frequency with time-Hilbert transform method,ballast codition(D=0.2018 m),without forward speed

表3 非线性阻尼系数拟合结果，压载（D=0.2018 m），无航速Tab.3 Results of nonlinear damping coefficients fitting,ballast condition(D=0.2018 m)

表4 非线性阻尼系数拟合结果，满载（D=0.2354 m），无航速Tab.4 Results of nonlinear damping coefficients fitting,full loading condition(D=0.2354 m)

4.3 基于参数识别的横摇阻尼估算方法结果

利用参数识别方法（PIT）估算横摇阻尼需要先预设阻尼模型和复原力模型（未知时），为便于对比分析，本文中复原力采用图4 所示的拟合模型，阻尼模型选用式（9）形式。所得归一化的自由横摇运动微分方程形式如下：

根据图4 拟合结果，式（41）中ε1、ε2已知，ωx使用式（27）估算。根据文献[11]可知，式（27）估算结果偏差不超过1%,则上式中未知参数有b1、b2、b3，本文选取最小二乘法、粒子群算法和Prony-SS 方法（为方便比较对方法进行了编号）识别上述参数。图10和表5展示了部分识别过程及非线性阻尼系数识别结果。

图10 基于Prony-SS方法的横摇衰减曲线重构，压载（D=0.2018 m），无航速Fig.10 Reconstruction of roll decay curves based on Prony-SS method,ballast condition(D=0.2018 m),without forward specd

表5 非线性阻尼系数参数识别结果Tab.5 Results of nonlinear damping coefficients based on PIT

最小二乘法参数初值采用[0；0；0]，粒子群算法的参数采用默认设置。图10（a）展示了使用Prony-SS 方法时对横摇衰减数据构成的Hankel 矩阵进行奇异值分解（singular value decomposition）后获得的奇异值量级分布。根据奇异值量级大小可判断重构的衰减曲线近似表达式的主要项，图10（b）为取前11个奇异值对应的项重构出的近似表达式与实验数据的对比，可见重构精度良好。

表5 中的非线性阻尼系数识别结果为未经调试检验的结果，均与基于分段线性假设的方法识别结果存在较大差异，说明基于参数识别的方法可能存在多解的现象。该类方法结果受识别过程中参数设置的影响，为排除多解和得到足够精度的阻尼系数，既要经验也需要根据每条实验数据反复比对调试。

5 横摇阻尼估算方法特点对比

横摇阻尼评估是船舶横摇运动预报简化的重要步骤，由于流动的复杂性、船舶各自由度运动的耦合，式（1）形式的运动方程只是对真实横摇运动的近似。好的评估结果应能提取出横摇阻尼的内在特性，即真实的横摇运动应为该特性叠加上各种随机因素的耦合结果。从表3～5 的阻尼参数拟合结果来看，基于分段线性假设的方法一致性较好，基于参数识别的方法结果差别较大，显示横摇阻尼估算结果准确度与使用的方法强相关。

5.1 横摇阻尼评估效果

为评估各方法的阻尼估算效果，参照式（39）定义误差R，表达式如下：

式中，φi为评估得到的阻尼模型参数代入横摇运动微分方程后，使用四阶Runge-Kutta法求得的横摇曲线衰减数据（本文取横摇幅值数据），φexp,i为对应实验数据。各方法的R值见图11，可见R值均不为0。R不为0 有三方面原因：存在测量误差及噪声，阻尼模型不够精确，方法自身引入误差。

图11 DTC标模不同数据处理方法的误差RFig.11 Comparison of roll prediction errors by roll damping of DTC model from different estimation methods

5.2 基于自由横摇衰减曲线的横摇阻尼评估方法误差来源

对于从自由衰减曲线估算横摇阻尼而言，有效横摇周期越多越有利于获取准确的横摇阻尼信息。然而，一般小振幅横摇测量噪声与误差占比大，大幅度横摇运动衰减太快，有效横摇周期少，如图5所示。为尽可能提高横摇阻尼评估精度，需要选择对测量噪声与误差不敏感且评估过程误差引入少的方法。

基于分段线性假设的估算方法涉及的曲线拟合、求导过程都会产生误差，表6所示为各方法曲线拟合与求导状况。对于消灭曲线法与改进型消灭曲线法，分段线性假设后，为了能实现消灭曲线拟合（如式（12）），一般假设所有周期频率都为横摇自然频率，然而，如图9（b）所示，真实的横摇频率会随时间不断变化，即式（15）、式（16）与式（19）都是近似等式，引入了额外误差。

表6 基于分段线性假设的估算方法涉及的曲线拟合及求导次数Tab.6 Times of curve fitting and derivation of estimation methods based on piecewise linear assumption

基于参数识别的估算方法涉及的阻尼模型选择同样会引入误差，根据最优化理论，非线性最小二乘法、粒子群对测量误差与噪声敏感且寻优结果可能为局部最小，即不一定是真实结果，且参数设置对结果影响大（如表5 所示）。Prony-SS 方法重构曲线后，求解阻尼参数的线性方程组为超定方程组，解该类方程组用到线性最小二乘法，存在唯一解，不存在该问题，但只有阻尼模型足够精确其参数识别结果才可信（作者使用自定义阻尼参数生成衰减曲线后，使用Prony-SS方法验证后得出此结论）。

6 结 论

本文以获取高精度阻尼预报方法为目标，对9种（渐近方法未考虑，其阻尼估算过程复杂，自身误差项多）横摇阻尼评估方法的原理进行分析并利用DTC 标模横摇衰减数据进行验证后，得出以下结论：

（1）基于分段线性假设的横摇阻尼估算方法物理意义明确，不需预设阻尼模型，但涉及的曲线拟合与求导过程易引入额外误差；除消灭曲线法外，其他方法均能将非线性横摇力臂纳入考虑；真实的横摇频率随时间不断变化，消灭曲线关系式（15）、改进型消灭曲线关系式（16）和式（19）只是近似成立；Froude能量法对测量噪声与误差不敏感，自身引入误差最少，具有良好性能。

（2）基于参数识别的横摇阻尼估算方法直接从微分方程与实现数据出发，原理清晰，但需预设阻尼模型且识别过程为黑箱子；非线性最小二乘优化、粒子群优化迭代参数对寻优结果影响大即可能为局部最小；Prony-SS方法存在唯一解，但只有阻尼模型足够精确其参数识别结果才可信。

（3）推荐使用对测量噪声与误差不敏感[11]，只有一次曲线拟合、无求导、物理意义明确、识别过程简洁的Froude能量法拟合出等效阻尼曲线，若能找出足够精确的阻尼模型，再利用Prony-SS方法获得更精确的阻尼模型参数。

横摇阻尼估算是一项复杂的工作，本文主要研究了基于衰减曲线估算横摇阻尼的方法，并未对不同影响因素下（航速、船型参数、尺度效应、流动记忆效应等）获得的自由衰减曲线是否能提供足够全的信息、是否反映了横摇的本质进行讨论，因此使用以上方法估算横摇阻尼时应予以考虑。