【《船舶非线性运动数值模拟方法》6700字】

船舶非线性运动数值模拟方法目录TOC\o"1-3"\h\u10930船舶非线性运动数值模拟方法 1141301.1.1相对于惯性坐标系的船舶运动 2235071.1.2相对于附体坐标系的船舶运动 2227201.1.3两个坐标系之间的互相转换 215331.2横摇运动数学模型 394491.2.1经典1.5自由度横摇运动方程 3236011.2.2改进的1.5自由度横摇运动方程 388041.3纵荡运动数学模型 5212731.4操纵-耐波统一模型 683821.4.1船舶的操纵性运动模型 770261.4.2船舶的耐波性运动模型 9142361.4.3操纵-耐波性统一模型 10227211.5波浪力计算方法 1113001.6四阶龙格-库塔方法 13160671.7小结 13为了能够更好地描述船舶在航行中的运动，本研究采用描述船舶运动较为常用的两个坐标系：惯性坐标系（即大地坐标系）与附体坐标系统（即随船坐标系），如图1.1所示。其中，表示相对于地球固定的惯性坐标系，是，作为研究过程中的基准坐标系，它用于表示船舶的绝对运动。分别以正北、正东、朝向地心为轴、轴和轴的正方向，单位向量分别记作、和。附体坐标系将船舶的前后左右的中点取为原点，分别以沿船长指向船艏方向、沿船宽指向右舷方向和沿船深指向龙骨方向为为轴、轴和轴正方向，其单位向量分别定义为、、。图1.1两种坐标系统定义Fig.1.1Twocoordinatesystemdefinitions1.1.1相对于惯性坐标系的船舶运动船舶相对于惯性内的运动可以通过船舶位置以及转角来描述，而船舶位置即为附体坐标系原点坐标，船舶转角即为附体坐标系的坐标轴指向。定义为船体在惯性坐标系下的速度矢量，可以将其分解成以下沿三个坐标轴的分量之和： （1.1）1.1.2相对于附体坐标系的船舶运动在附体坐标系下船舶的运动包括平动与转动，平动速度可以分解成沿x，y，z轴的纵荡、横摇以及垂荡方向上的速度，而转动速度可以分解成绕三个坐标轴的横摇、纵摇以及艏摇方向的角速度： （1.2） （1.3）1.1.3两个坐标系之间的互相转换在惯性坐标系中和附体坐标系下的运动变换关系如公式所示ADDINNE.Ref.{4B61B644-3540-438A-933D-5608958A120D}[62]： （1.4） （1.5）式中，皆为变换矩阵，可表示为： （1.6） （1.7）1.2横摇运动数学模型1.2.1经典1.5自由度横摇运动方程船舶的横摇运动是研究船舶倾覆最重要的指标之一。根据单自由度振动系统来建立船舶的横摇运动方程。船舶所受的横摇力矩由以下几种构成：（1）惯性力矩船舶横摇惯性矩可以用下式来表示： （1.8）其中表示自身惯性矩，表示附加惯性矩。（2）阻尼力矩船舶的横摇阻尼力矩受船舶吃水、减摇装置、船舶航速等因素的影响，对船舶横摇运动的模拟起着相当重要的作用。阻尼力矩可以用下式来表示： （1.9）其中和分别为一次项和三次项阻尼系数。（3）回复力矩船舶在做横摇运动时，会产生了复原力矩。复原力矩可以表示为船舶的复原力臂与排水量的乘积，其中与横摇角、波面升高以及时间t关： （1.10）根据GabrieleBulian[63]等人提出的数学模型，船舶横摇回复力臂可以用下式来表示：（1.11）上式中：表示船舶初稳性高，表示初稳性高随时间的变化项，和分别表示横摇回复力矩的三次项和五次项系数，可以由船舶的静稳性曲线拟合得到（4）波浪扰动力矩波浪的扰动力矩主要包括基于Froude-Krylov假设的入射波压力与绕射力。将上述力矩代入到单自由度振动方程中，建立船舶单自由度横摇运动方程： （1.12）1.2.2改进的1.5自由度横摇运动方程为更准确地描述船舶在海洋中的真实运动，需要考虑船舶其它自由度方向的运动与横摇运动的耦合，才能更好地描述横摇恢复力臂在波浪中的变化。以往的1.5自由度运动模型中，常常采用下面两种方法来考虑垂荡、纵摇跟横摇之间的耦合。一是静力平衡方法，即认为当船舶在任何时刻、具有任何横摇角时，其垂荡和纵摇都恰好使得船舶在两个自由度上的静水力和弗雷德-克雷洛夫力跟重力相平衡[64]。但这种方法只是一种假设，并没有反应船舶真实的运动状况。二是垂荡和纵摇直接采用垂荡-横摇-纵摇三自由度线性运动方程的结果来计入[65]。这种方法考虑了垂荡、纵摇对横摇的影响，但未考虑横摇对垂荡、纵摇的影响。使用这两种方法对横摇运动方程进行求解，有时候能得到令人满意的结果，但结果的精确程度受到船型、航行条件等的制约，这是因为垂荡、纵摇跟横摇之间的相互耦合往往是不能忽略的，而上述两种方法均不是真实的耦合方式。本文改进了第二种方法，为了考虑横摇对垂荡、纵摇的影响，首先求解上述三自由度非线性运动方程，其中，非线性体现在两个方面：第一，横摇运动的阻尼采用“一次加三次”非线性阻尼；第二，静水力和弗雷德-克雷洛夫力沿由入射波面和船舶瞬时位置确定的湿表面进行积分。然后用上述方程的解对三自由度线性运动方程解出的垂荡、纵摇进行修正，主要修正两方面：平衡位置和幅值。修正的方法是，平衡位置和幅值直接采用非线性运动方程的结果，而频率和相位不作修改。接下来，使用修正后的垂荡和纵摇计算船舶在波浪中的回复力臂，其中，QUOTEΦ是横摇角，是船舶重心在波浪中的位置，在波峰处为0，波谷处为，并将其拟合成关于、的分离变量形式的表达式[66]。与以往研究[67]不同的是，在斜浪中的表达式中计入了由静水力与弗雷德-克雷洛夫力之和引起的定常横倾力矩，即表达式中关于横摇角Φ的偶数次幂项，这不仅将使拟合结果更好，且更加符合实际。如上所述，本文对垂荡和纵摇运动需要做必要的修正。在某些易于发生非线性横摇的航行条件下，如果船舶发生了非线性大幅横摇运动（例如纯稳性丧失），则垂荡和纵摇的解将与垂荡-横摇-纵摇三自由度线性运动方程解出的垂荡、纵摇的解有显著的差别，这种差别主要体现在运动幅值的改变和平衡位置的移动。另外，垂荡、纵摇运动的频率成分常常也变得不仅仅有遭遇频率，但为了简便起见，本文的研究中不予计入；此外，相位的变化相当小，本文也不予考虑。产生这些差别的原因是非线性横摇运动对垂荡、纵摇运动产生了显著的影响，且不可忽略。下面以垂荡为例，对修正方法进行说明。设由三自由度线性运动方程解出的垂荡为：（1.13）而三自由度非线性运动方程解出的垂荡（稳定后）为：（1.14）则将三自由度线性运动方程解出的垂荡修正为（1.15）由于船舶在波浪上航行时，静水力与波浪入射力之和QUOTEfspkj（j=1,2,…,6）是关于时间的周期函数，且仅由船舶重心在波浪中的相位QUOTEΨ和横摇角QUOTEx4确定。将该力在横摇方向上的分量QUOTEfspk4t,x4记为，则有：（1.16）从而可得波浪中横摇回复力臂的表达式：（1.17）其中，QUOTEΔ是排水量，QUOTESt是由船舶姿态和波浪位置确定的湿表面；QUOTEn'j（j=1,2,…,6）是广义单位法向矢量；QUOTEps是静水压力和入射波一阶动压力之和。基于多项式和傅立叶展开式对此拟合上式可得：（1.18）式中，曲线关于的函数关系可由其一次项、三次项、五次项和七次项叠加拟合而成，K1、K3、K5和K7分别为对应的系数，此外p和q分别为对应的修正系数，会随波浪条件、船舶载况等因素发生变化。以上是GZ在规则波中的拟合形式。在假定不规则波近似为规则波的基础上，提出不规则波GZ的近似数学表达式。当将规则波（波高：H，波周期：T）中的GZ拟合形式为式（1.18）时，不规则波（有义波高：H，谱峰周期：T）中的GZ可以用以下形式表示：（1.19）（1.20）其中，为波浪激励过程的第个子波波高（振幅）。显而易见，式（1.19）比式（1.11）更好地描述横摇恢复力臂在随机波浪中的变化。1.3纵荡运动数学模型为了建立一个描述随浪中船舶纵荡运动的理想模型，需要做一些简化和假设：（1）船舶是左右舷对称的；（2）流体是无粘的、不可压缩的理想流体；（3）基于势流运动；（4）基于微幅波假说。船舶在随机波浪激励下纵荡方向的运动方程为：（1.21）其中：为惯性坐标系下船舶重心的纵向坐标，为纵向坐标的一阶导数，为船舶质量，为纵荡方向的附加质量，为纵荡方向的波浪扰动力，为船舶阻力，为螺旋桨推力。阻力和推力均采用多项式近似表示：（1.22）（1.23）其中：为推力减额，为海水密度，为螺旋桨转速，为螺旋桨直径，为推力系数。（1.24）其中：进速系数将式（2-4）带入式（2-3）可得：（1.25）（1.26）（1.27）（1.28）其中：为伴流系数。1.4操纵-耐波统一模型以上两节均只介绍了船舶的单自由度运动方程，但船舶在航行时，不仅受波浪的作用做多自由度耦合运动，而且一直要保持航向的稳定，即船舶做操纵性运动。船舶的操纵性为艏摇自由度带来了回复力，当船舶在斜浪中航行时，若不考虑船舶的操纵性运动，在周期性外激励中艏摇角会越来越大，模拟结果将不符合实际情况。为准确模拟船舶在波浪中的运动，需要将船舶的耐波性运动和操纵性运动联合起来，建立一个统一的运动模型，以考虑两者之间的耦合作用，保持航向稳定的同时，模拟船舶的运动响应。1.4.1船舶的操纵性运动模型当前研究船舶的操纵性运动可以从两个方面入手—整体型和分离型。整体性操纵模型是将船舶看作一个整体，船体、桨、舵均包含在内，运动模型中的外力即船舶整体受到的各种外力，通过求解整体运动模型得到船舶的运动响应。分离型模型则是分别建立船舶船体、桨、舵的运动模型，各项外力也分开计算，然后再考虑模型之间的耦合作用，进而求解船舶运动响应。由于分离型各项外力的求解更加简单便捷，因而应用更加广泛，具有代表性的分离型模型即为MMG模型，本文建立的操纵性模型即为模型。MMG模型中，假设船舶航行的水域为无限深，忽略浅水的影响，不考虑波浪带来的波浪力变化，不考虑船舶纵摇和垂荡的影响，建立船舶三自由度纵荡、横荡和首摇运动模型，并考虑舵的影响：（1.29）其中：、分别为纵荡、横荡方向的附加质量，为艏摇惯性矩，为艏摇附加惯性矩，为横摇附加惯性矩。、分别为横荡速度、艏摇角速度。、、分别作用在为纵荡、横荡、艏摇方向的船体力，、、为纵荡、横荡、艏摇方向的舵力，、为横荡、艏摇方向的波浪力。考虑到船舶与波浪中运动时，与水质点之间存在较为复杂的相互作用，故做以下简化处理：①模拟运动时忽略船舶水动力系数随船舶航速的变化。②忽略船舶在波浪中的摇荡运动对舵和桨的影响。上述模型中，船体力分为惯性船体力和粘性船体力。其中惯性船体力分为船体惯性力与附加惯性力，各系数的回归估算公式[68]如下：（1.30）（1.31）（1.32）（1.33）（1.34）其中：为船长，为型宽，为平均吃水，为方形系数，为排水量，为初稳性高，为船舶重心高度，为横摇固有周期。而粘性船体力则采用井上模型[69]，其各种无量纲粘性系数采用经验公式估算的方法[70]确定。上述模型中舵力模型如下式表示：（1.35）其中：为舵面积，为舵角（右舵为正），和为舵力作用点的垂向坐标和纵向坐标。为操舵诱导的船体横向力与舵力的比。（1.36）为船体横向力的作用点到船舶重心的距离，（1.37）为舵力系数，（1.38）其中：为舵的展弦比。为舵力减额系数，采用下述公式计算：（1.39）图1.2舵来流和有效攻角Fig.1.2Rudderflowandeffectiveattackangle为舵来流速度，采用下述公式计算：（1.40）（1.41）（1.42）其中：为舵处流速的纵向分量，为舵处流速的横向分量，为整流系数，为舵桨伴流分数比，为舵桨相互作用系数。为有效舵攻角，船舶运动时，随着水流的偏转，舵的有效攻角小于实际舵角。（1.43）1.4.2船舶的耐波性运动模型船舶的耐波性运动采用六自由度运动模型描述，根据基本动力学方程和微幅波理论得到如下耐波性运动模型：（1.44）其中表示船舶的质量或惯性矩，为波频接近无限时船舶的附加质量或附加惯性矩,为船舶第j个自由度的速度或角速度。表示横摇等效线性化粘性阻尼系数，一般可由静水中的横摇自由衰减曲线结果拟合得到。式中和为船舶回复力、力和绕射力，为舵的横摇力矩。模型中还在横摇自由度中计入了横荡船体力诱导的横摇力矩，为横荡船体力等效作用点距船舶重心的纵向距离。由于船舶在水平面上的横荡、纵荡、艏摇并没有回复力，因此船舶的耐波性运动假设船舶在其平衡位置上做微幅的摇荡作用，且不会严重偏离其平衡位置。1.4.3操纵-耐波性统一模型虽然船舶的耐波性运动包含了六个自由度，但仅靠耐波性方程无法很好地描述船舶在波浪中的运动。若采用线性耐波性方程，船舶在斜浪中的运动将无法被模拟：在斜浪中，船舶会发生艏摇运动，且由于没有艏摇回复力，船舶的艏摇运动将做周期运动，且平衡位置非零，这将打破耐波性方程中纵荡自由度中推力与阻力平衡的假定。若再考虑船舶的纯稳性丧失、骑浪/横甩等非线性运动，即使在纵浪情况下，船舶依然会产生横摇、艏摇、横荡运动，此时单靠耐波性方程无法合理地模拟船舶运动。因此需要计入操纵性方程，通过操舵为船舶产生一个持续的艏摇回复力，以此模拟船舶在波浪中的运动。船舶的耐波性方程是一个高频运动，而船舶在波浪中的转向与航速的变化是十分缓慢的，因此操纵性方程可以看作一个低频运动。根据高频与低频运动的独立性，船舶在艏摇自由度的运动可以看作由耐波性方程中的高频振荡与由舵、桨以及船体力引起的操纵性方程中的低频振荡的代数和。另外，在较小的时间区间内，船舶重心在横向、纵向、垂向运动的矢量和可以看作匀速直线运动。根据上述可知，高频的耐波性运动与低频的操纵性运动联系较弱，两者之间的耦合作用仅仅体现在以下三个方面：第一，在每个时间区间内，耐波性运动方程与操纵性运动方程将单独求解，且两者的步长不同，高频运动时间区间短，低频运动时间区间长；第二，耐波性运动的速度矢量与操纵性运动保持一致；第三，操纵性运动每个时间区间的初始艏摇角，除了包括自身运动方程在上个运动区间求解的结果之外，还要加上高频的耐波性运动在上个操纵运动时间区间内累计所得的艏摇角。图1.3操纵-耐波统一模型流程图Fig.1.3FlowChartoftheUnitedModelofmaneuveringandseakeepingmotion上图展示了操纵-耐波统一模型模拟船舶在波浪中的运动过程，船舶耐波性运动的时间区间为dt，操纵性运动为dT，dt可以取dT的1/n，n>1。另外船舶在每个时间区间内的方程系数，例如附加质量矩阵、附加阻尼矩阵均因航速、航向等的变化而发生变化，因此本文通过插值的方式获得，即预先算得船舶模型在多航速、全浪向下的方程系数，然后根据瞬时航速与迎浪角，建立插值函数来求取。1.5波浪力计算方法求解波浪力对于船舶的非线性运动模型至关重要。波浪的绕射力和Froud-Krylov力均需要计算船舶与波浪位置决定的船舶瞬时湿表面积。本文中单自由度运动方程中的波浪力传递函数借助AQWA软件进行计算。根据目标船型的型线资料建立相应的船体模型并划分单元，通过AWQA-LINE软件，依据三维辐射与衍射理论，获得F-K力、绕射力等非线性力的传递函数。操纵-耐波统一模型中的波浪力则是则是依据三维势流理论，沿船舶的瞬时湿表面积进行积分来计算的。求解势流问题的思路如下：首先通过根据边界条件和初始条件求解拉普拉斯方程得到流场的速度分布，接着求解流场的伯努利方程得到其压力分布，然后再沿结构的表面积进行积分，得到流场对结构的作用力ADDINNE.Ref.{CDF28466-7F3D-4164-A836-486E6A92C23A}[71]。求解流场速度势的过程如下。首先假设流体理想流体，并引入速度势函数，流场中的拉普拉斯方程如下： （1.45）流体速度可以表示为速度势相应方向上的偏导： （1.46）拉普拉斯方程作为一个二阶线性微分方程，根据初始条件和边界条件可以求的唯一的速度势。对于较为复杂的边界条件，可以采取数值方法求解。得到速度分布之后，通过求解伯努利方程来得到压力分布。将散射势分为绕射势与辐射势，并将代入到流场的伯努利方程中，求得其动压力分布： （1.47）根据上式可以推得船舶所受到的波浪力：（1.48）或记为： （1.49）其中： （1.50）上式中积分面积为船舶在波浪中的瞬时湿表面积，将代入到式（），得到下式： （1.51） （1.52） （1.53）上式中为F-K力，即一阶波浪动压力，是绕射力，两者共同构成了波浪扰动力。1.6四阶龙格-库塔方法龙格-库塔法在可采用显式迭代或隐式迭代两种方式进行求解。作为求解常微分方程常用的数值手段，龙格-库