CN107140110B - 一种船舶大幅横摇运动非线性阻尼系数识别方法 - Google Patents

Abstract

本发明公开了一种船舶大幅横摇运动非线性阻尼系数识别方法，包括如下步骤：确定船舶非线性横摇运动方程，确定船舶横摇角幅值随时间变化解析表达式，推导相邻半个周期的横摇角幅值变化解析式，绘制船舶横摇消灭曲线，并进行曲线拟合，计算非线性阻尼系数c₁和c₂的估算值，进行阻尼系数修正，计算阻尼系数最终识别值，本方法能够克服现有船舶横摇非线性阻尼系数识别技术的不足，具有较高的识别精度。

Description

一种船舶大幅横摇运动非线性阻尼系数识别方法

技术领域

本发明涉及一种船舶大幅横摇运动非线性阻尼系数识别方法。

背景技术

船舶作为航运业重要的交通工具，其在波浪中的运动性能，特别是横摇运动性能备受关注。船舶横摇阻尼准确估算是正确预报波浪中横摇运动的前提。小角度横摇时，船舶横摇运动可以用线性方程描述，但是当横摇幅值增大时，非线性效应必须考虑。横摇非线性主要表现在回复力矩的非线性和阻尼的非线性，船舶横摇非线性回复力矩可由流体静力学计算得到。但是对于船舶非线性横摇阻尼项，目前尚未有完备的理论计算方法。因此，开展船舶大幅横摇非线性阻尼系数识别，实现船舶在波浪中的横摇性能的准确预报，对于指导船舶设计、降低船舶横摇幅值，进而保障船舶航行安全具有重要意义。

对于船舶横摇非线性阻尼系数识别，现有船舶耐波性评估的势流理论不能预报横摇运动中由于摩擦、旋涡和流动分离等因素产生的粘性阻尼；基于模型试验获得的经验公式在实际使用中较为方便，但仅对特定船型适用。因此在工程实践中常用衰减试验的方法来确定船舶横摇非线性阻尼系数。

国内外学者对基于船舶衰减试验的横摇非线性阻尼系数识别方法开展了研究，并取得了一定成果。Mathisen和Price提出了估算衰减横摇和强迫横摇试验阻尼系数的方法，研究发现线性加平方项阻尼模型比线性加立方项阻尼模型具有更好的优越性。Roberts提出了一种能量法，由模型试验测量的横摇衰减曲线来确定船舶横摇阻尼系数。首先给出了反映船舶横摇能量平均衰减情况的能量损失函数，结合某假定的横摇阻尼力矩形式，假设在一个周期内横摇角等幅变化，将能量损失函数表示成未知横摇阻尼系数的解析形式，然后根据模型试验横摇衰减曲线的一系列峰值点，通过最小二乘法确定横摇阻尼力矩系数，该种方法可考虑大幅横摇时船舶横稳心高的变化，但该方法数值实现时需要较多的横摇周期数，否则数值精度会受到一定影响。Bass和Haddara对横摇阻尼力矩的各种形式适用性进行了研究，引进了两种数值方法来确定阻尼力矩表达式中的系数，其中能量法基于横摇机械能守恒的原理利用横摇衰减曲线连续时间历程来获得阻尼系数，适合于衰减过程较快和初始横摇幅值较大的横摇过程。Spouge对于不同学者提出的各类不同阻尼系数识别方法和结果进行了总结比较，并探讨了各种方法的识别精度。Chan等将徐兆提出的一种新的渐进法应用于船舶非线性横摇阻尼系数估算中，该方法适用于大角度横摇以及强非线性回复力矩。李红霞等根据衡量耗散的观点利用横摇试验衰减曲线，提出了一种非线性阻尼识别方法，由实际能量衰减函数与估算能量衰减函数之间方差最小来确定非线性阻尼系数。

综合国内外研究现状可知，不同学者提出了不同的基于自由衰减数据的船舶横摇非线性阻尼系数识别方法。其中，能量法和渐进法是主流的研究方法。该两种方法能实现大幅横摇运动和强非线性回复力矩条件下非线性阻尼系数的识别；但是阻尼系数识别结果受初始横摇幅值影响较大，特别是在大横摇幅值时，阻尼系数识别存在相对较大的误差。因此大横摇幅值和强非线性回复力矩下船舶自由横摇非线性阻尼系数的精确识别仍有待提高。

发明内容

为解决上述技术问题，本发明提供了一种船舶大幅横摇运动非线性阻尼系数识别方法，以达到提高阻尼系数识别精度的目的。

为达到上述目的，本发明的技术方案如下：

一种船舶大幅横摇运动非线性阻尼系数识别方法，包括如下步骤：

步骤1：确定船舶非线性横摇运动方程，

式中，x：横摇角；j：描述船舶回复力矩的幂次数项，正整数；M：幂次数项最大值，正整数；k：船舶静水回复刚度系数；

步骤2：确定船舶横摇角幅值随时间变化的解析表达式，

步骤3：推导相邻半个周期的横摇角幅值变化解析式，

令其中n为正整数，T₀为横摇衰减运动固有周期，则有 为经过半个周期的横摇幅值变化，亦即相邻横摇幅值变化，将该等式带入表达式可得：

式中，a：横摇角幅值；t：时间；p₁₁、q₁₁、p₁₂、q₁₂：拟合系数；相位角；Φ₀：相位角变化率表达式的零阶项；x₀：初始横摇角；

步骤4：绘制船舶横摇消灭曲线，并进行曲线拟合，

根据船舶静水自由横摇衰减数据，采用间隔半个周期的横摇角幅值差与相邻两次的平均横摇角幅值绘制横摇消灭曲线，采用最小二乘法对数据进行拟合，可得到船舶相邻横摇幅值变化表达式：

其中，λ₁和λ₂是拟合系数；

步骤5：令步骤3推导的相邻半个周期的横摇角幅值变化解析式 与步骤4中基于消灭曲线得到的船舶相邻横摇幅值变化表达式中参量a和a²的系数相等，可得到非线性阻尼系数c₁和c₂的估算值：

步骤6：阻尼系数的第一次修正，

将步骤5得到的阻尼系数初次估算值和带入船舶横摇运动方程 中，初始条件不变，采用龙格库塔方法进行数值求解，可仿真生成自由横摇衰减曲线，求得新的横摇消灭曲线：

再次根据步骤5中阻尼系数估算公式，可得阻尼系数的第二次估算值：

前两次阻尼系数估算值的相对误差为：

步骤7：阻尼系数的第二次修正，

重复步骤6，可以求得第三次阻尼系数估算值和以及后两次阻尼系数估算值之间的相对误差和

步骤8：计算阻尼系数最终识别值，

根据步骤5得到阻尼系数估算值和步骤6、步骤7中得到的阻尼系数估算值相对误差η₁、η₂、和得出非线性阻尼系数最终识别值计算公式：

上述方案中，所述步骤1中，船舶非线性横摇运动方程的确定方法为：对强非线性系统，其运动可由下列微分方程描述：

式中，ε为小参数，函数g(x)和满足如下关系：g(-x)＝-g(x)，

根据已有研究，横摇非线性阻尼可采用线性加平方形式或线性加立方形式，以线性加平方阻尼为例，则回复力矩采用奇次多项式函数表示： 于是得到船舶横摇运动方程：

上述方案中，所述步骤2中，确定船舶横摇角幅值随时间的变化解析式具体包括以下步骤：

(1)运动方程广义解，

根据摄动理论，可得运动方程的广义渐进解：

式中，a是振幅，是相位角，均为时间的缓变函数，x₁,x₂,…,x_m-1是相位角的周期函数，其周期为2π；a和由如下微分方程决定：

(2)方程一阶近似解，

为简化计算，取方程一阶近似解代替方程广义解：

式中，

(3)求解A₁(a)和解析表达式，

令ε＝0，对非线性系统运动方程进行积分后，方程两边同乘以并积分，假定运动初始条件为x(0)＝a,经过变换，可得其中，

由方程一阶近似表达式可得运动变量表达式，把函数g(x)展开成的ε幂级数 将和g(x)带入到运动方程中，使等式两边ε¹项的系数相等，可得到如下表达式：

，

方程两边同乘以然后积分，令可以求出A₁表达式 令则可以求出Φ₁表达式，至此方程的一阶近似解可求出；

(4)推导横摇角随时间变化表达式，将线性加平方阻尼表达式带入到(3)中A₁(a)表达式中，可得：

式中，

P(a)和Q(a)是幅值a的函数，可以用多项式函数和近似表示：

其中，系数p₁₁，p₁₂，q₁₁和q₁₂通过最小二乘法进行数据拟合得到，于是可得横摇角随时间的变化：

通过上述技术方案，本发明提供的船舶大幅横摇运动非线性阻尼系数识别方法能够克服现有船舶非线性阻尼系数识别技术的不足，具有较高的识别精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。

图1为本发明公开的一种船舶大幅横摇运动非线性阻尼系数识别方法的流程图；

图2为本发明的实施例公开的非线性回复力矩曲线图；

图3为本发明的实施例公开的由数值模拟生成的自由横摇衰减曲线图；

图4为本发明的实施例公开的横摇消灭曲线图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。

本发明提供了一种船舶大幅横摇运动非线性阻尼系数识别方法，如图1所示，本实施例中没有使用实际船舶模型的衰减试验数据，而是基于已知的船舶横摇运动方程，通过数值方法模拟生成船舶静水自由横摇衰减曲线，基于该衰减数据对已确定的阻尼系数进行识别，以描述本发明实施例中的技术方案。

步骤1：确定船舶非线性横摇运动方程，

此例中，船舶非线性运动方程已知：(初始条件x(0)＝0.3，)。

其中，回复力矩表达式为g(x)＝2.25x-x³，具有强非线性，如图2所示。船舶的非线性阻尼系数已知：C₁＝0.15，C₂＝0.2。

步骤2：确定船舶横摇角幅值随时间变化表达式，

根据回复力矩g(x)表达式，可得V(x)＝1.125x²-0.25x⁴，继而求得

根据表达式以及P(a)和Q(a)表达式，通过最小二乘法进行数据拟合，得到拟合系数：p₁₁＝0.49572，p₁₂＝0.03359，q₁₁＝0.6381和q₁₂＝0。

步骤3：推导相邻半个周期的横摇角幅值变化，

根据横摇固有周期表达式计算得到T₀＝4.25316s。

步骤4：绘制船舶横摇消灭曲线，并进行曲线拟合，

基于给定的船舶横摇运动方程，通过数值方法模拟生成自由横摇衰减曲线，如图3所示。绘制相应消灭曲线，见图4。根据最小二乘法对消灭曲线数据进行拟合，得到拟合系数λ₁＝0.15573和λ₂＝0.28290。

步骤5：阻尼系数估算值计算，

根据步骤2、3和4计算得到的p₁₁，p₁₂，q₁₁，q₁₂，T₀，λ₁和λ₂值，得到非线性阻尼系数估算值和

步骤6：阻尼系数第一次修正，

将和作为新的阻尼系数带入到横摇运动方程中，重复上述计算过程，得到阻尼系数第二次估算值和阻尼系数两次估算值之间的相对误差为η₁＝-1.45874％和η₂＝0.36716％。

步骤7：阻尼系数第二次修正，

重复步骤6，得到阻尼系数第三次估算值和以及后两次阻尼系数估算值之间的相对误差和

步骤8：计算阻尼系数最终识别值，

根据阻尼系数最终识别值计算公式，可得阻尼系数最终识别值c₁＝0.15000和c₂＝0.20000，与方程中真实阻尼系数μ₁＝0.15和μ₂＝0.2相等，识别精度达到100％。

上述实施例对本发明的技术方案进行了具体展示。为验证本发明的阻尼系数识别方法的精度，基于文献公开发表的阻尼系数识别结果，本发明方法分别与Roberts能量法、Chan渐进法进行了比较，结果分别见表1和表2。从表中可知，本发明方法得到的阻尼系数识别结果要明显优于Roberts能量方法和Chan渐进法，具有较高的识别精度。

表1本发明法和Chan渐进法得到的阻尼系数识别结果比较

(横摇运动方程：初始条件x(0)＝0.866，)

表2本发明方法和Roberts能量法得到的阻尼系数识别结果比较

(横摇运动方程：初始条件x(0)＝0.698，)

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (3)

1.一种船舶大幅横摇运动非线性阻尼系数识别方法，其特征在于，包括如下步骤：

步骤1：确定船舶非线性横摇运动方程，

式中，x：横摇角；j：描述船舶回复力矩的幂次数项，正整数；M：幂次数项最大值，正整数；k：船舶静水回复刚度系数；

步骤2：确定船舶横摇角幅值随时间的变化解析式，

步骤3：推导相邻半个周期的横摇角幅值变化解析式，

令其中n为正整数，T₀为横摇衰减运动固有周期，则有 为经过半个周期的横摇幅值变化，亦即相邻横摇幅值变化，将该等式带入表达式可得：

式中，a：横摇角幅值；t：时间；p₁₁、q₁₁、p₁₂、q₁₂：拟合系数；相位角；Φ₀：相位角变化率表达式的零阶项；x₀：初始横摇角；

步骤4：绘制船舶横摇消灭曲线，并进行曲线拟合，

根据船舶模型静水自由横摇衰减数据，采用间隔半个周期的横摇角幅值差与相邻两次的平均横摇角幅值绘制横摇消灭曲线，对数据进行拟合可得到船舶相邻横摇幅值变化表达式：

其中，λ₁和λ₂是拟合系数；

步骤5：令步骤3推导的相邻半个周期的横摇角幅值变化解析式 与步骤4中基于消灭曲线得到的船舶相邻横摇幅值变化表达式中参量a和a²的系数相等，可得到非线性阻尼系数c₁和c₂的估算值：

步骤6：阻尼系数的第一次修正，

将步骤5得到的阻尼系数初次估算值和带入船舶横摇运动方程 中，初始条件不变，采用龙格库塔方法进行数值求解，可仿真生成自由横摇衰减曲线，求得新的横摇消灭曲线：

再次根据步骤5中阻尼系数估算公式，可得阻尼系数的第二次估算值：

前两次阻尼系数估算值的相对误差为：

步骤7：阻尼系数的第二次修正，

重复步骤6，可以求得第三次阻尼系数估算值和以及后两次阻尼系数估算值之间的相对误差和

步骤8：计算阻尼系数最终识别值，

根据步骤5得到阻尼系数估算值和步骤6、步骤7中得到的阻尼系数估算值相对误差η₁、η₂、和得出非线性阻尼系数最终识别值计算公式：

2.根据权利要求1所述的一种船舶大幅横摇运动非线性阻尼系数识别方法，其特征在于，所述步骤1中，船舶非线性横摇运动方程的确定方法为：对强非线性系统，其运动可由下列微分方程描述：

式中，ε为小参数，函数g(x)和满足如下关系：g(-x)＝-g(x)，

根据已有研究，横摇非线性阻尼采用线性加平方形式，则回复力矩采用奇次多项式函数表示：于是得到船舶横摇运动方程：

3.根据权利要求1所述的一种船舶大幅横摇运动非线性阻尼系数识别方法，其特征在于，所述步骤2中，确定船舶横摇角幅值随时间的变化解析式具体包括以下步骤：

(1)运动方程广义解，

根据摄动理论，可得系统运动方程的广义渐进解：

式中，a是横摇角幅值，是相位角，均为时间的缓变函数，x₁,x₂,…,x_m-1是相位角的周期函数，其周期为2π；a和由如下微分方程决定：

(2)方程一阶近似解，

为简化计算，取方程一阶近似解代替方程广义解：

式中，

(3)求解A₁(a)和解析表达式，

令ε＝0，对非线性系统运动方程进行积分后，方程两边同乘以并积分，假定运动初始条件为x(0)＝a,经过变换，可得其中，

由方程一阶近似表达式可得运动变量表达式，把函数g(x)展开成的ε幂级数 将和g(x)带入到运动方程中，使等式两边ε¹项的系数相等，可得到如下表达式：

，方程两边同乘以然后积分，令可以求出A₁表达式 令则可以求出Φ₁表达式，至此方程的一阶近似解可求出；

(4)推导横摇角随时间变化表达式，将线性加平方阻尼表达式带入到(3)中A₁(a)表达式中，可得：

式中，

P(a)和Q(a)是幅值a的函数，可以用多项式函数和近似表示：

其中，系数p₁₁，p₁₂，q₁₁和q₁₂通过最小二乘法进行数据拟合得到，于是可得横摇角随时间的变化：