非线性水动力导数的数值计算与研究

第41卷 第1期2017年2月武汉理工大学学报（交通科学与工程版）Journal of W uhan U niversity of Technology(T ransportation Science & Engineering)Vol. 41 No. 1Feb.2016非线性水动力导数的数值计算与研究&赵小仨u徐海祥1>2)(高性能船舶技术教育部重点实验室1:1武汉430063)(武汉理工大学交通学院2)武汉430063)摘要：针对船舶的非线性运动难以界定和非线性运动难以预报的问题，以供应船为研究对象，采用 C F D商用软件F L U E N T，结合动网格技术对大振幅平面运动机构试验进行数值模拟，通过对比不 同工况的流场压力云图，分析得出供应船水动力达到非线性的振幅范围.设计供应船非线性运动的试验方案，分别模拟不同频率时的大漂角斜航运动及大振幅的纯纵荡、纯横荡、纯首摇、组合运 动，拟合得到接近零频率的非线性水动力导数.关键词：非线性水动力；大振幅P M M试验；数值计算；供应船中图法分类号：U661. 33 doi：10. 3963/j. issn. 2095-3844. 2017. 01. 014〇引言船舶操纵性与船舶航行安全紧密相关，是重 要的水动力性能之一.近些年，国际海事组织（in-ternational m aritim e organization, IM O)前后颁布了 A. 751(18)和MSC. 137(76)号决议，针对船 舶操纵性的问题提出了明确的要求，并建议各国 政府机构按要求执行.SIM M A N2008和SIM- M AN 2014 的研讨会，评估T C F D(co m p u tatio n- al fluid dynam ics, C F D)方法预报船舶操纵性的 能力•第 25 届 IT T C(international tow ing tank conference，IT T C)操纵会议对现有的船舶操纵性预报方法做了总结.总之，船舶操纵性能越来越 受到造船界的重视[>3].水动力导数对船舶操纵性的预报至关重要.目前，通过平面运动机构试验（planar motion mechanism test,PM M)确定船舶水动力导数是最可靠的方法之一.从SIM M AN2008发布了针 对三个标准船模进行的P M M试验的实验数据以 来，国内外学者开始对C F D模拟P M M试验进行 探究•T u rnock等[4_12]用C F D软件模拟小振幅P M M试验，求取线性水动力导数；Toxopeus 等〜^建立非线性水动力模型，模拟大振幅P M M试验，求取非线性水动力导数.虽然许多学者对数值模拟P M M试验做了大 量研究工作，但是迄今没有一个定量的标准来判断船模的运动是否达到非线性范畴，相关文章也 较少.评判船舶的运动是否达到非线性，不仅取决 于运动参数，还与船型等因素有关.文中将以供应 船为研究对象[17]，通过数值模拟船模不同运动幅 值的P M M试验，分析出供应船水动力达到非线性的运动幅值范围.在此基础上，设计试验工况，计算零频率附近的非线性水动力导数.1数学模型研究船舶在大振幅下的操纵运动，用线性水 动力模型很难准确的表达船舶所受到的水动力，为了更准确的描述船舶的运动，须考虑运动状态 的非线性项[18].根据经验，在粘性类流体动力和力矩的泰勒级数展开式中保留至三阶项，对描述 船舶在常速域中的运动已足够精确.1)流体惯性力（矩收稿日期：2016-12-14赵小仨（1989—）：女，工学硕士，实验员，主要研究领域为船舶水动力研究国家自然科学基金项目（61301279, 51479158)、中央高校基本科研业务费专项资金项目（163102006)资助第1期赵小仨，等：非线性水动力导数的数值计算与研究• 71 •Xi = X,ii -Y tvr -Y rr2Yj = Y tv + X.ur + Y,f Mj = + d — X A) uv(1)式中：〃，w r 分别为船舶的纵向、横向速度与转首 角速度.2) 粘性类流体动力（矩)X h ，Y h ，M h.X h = X(u) + X mv 2 + X rrr2 +n1 y rr ++ Y rrrr 3 +<Y ^r + Y ^v(2)M h = M ，+ M ，+ M _w 3 + M m.r3 +、M mrv 2 r + M IT V r 2 v3) 水动力导数的无因次化采用M M G 模型系统建议的以为参考面积进行无因次化.2 平面运动机构试验P M M 试验是约束模试验的一种.是通过测 量船模所受到的水动力和力矩，求得计算船舶操 纵运动所需的各种加速度导数、速度导数和耦合 导数.P M M 试验有小振幅和大振幅之分，前者只 能测定线性水动力导数，后者可以测定非线性水 动力导数.本文主要模拟P M M 试验的以下几种 运动形式：定漂角斜航；纯横荡运动；纯纵荡运动； 纯首摇运动；组合运动.3数值计算方法3.1计算模型计算模型为一艘75 m 供应船，缩尺比为 1:20.船模几何参数见表1.三维模型见图1.确定.计算域尺寸船首上游取1. 5倍船长，船尾下 游取3倍船长，船两侧取2倍船长，水深方向取 8. 3倍吃水，见图3.a )纯纵/横荡运动b )纯首摇运动图3供应船的计算域3.4划分网格以及验证收敛性为了保证网格质量，采用分块全结构化网格， 并在首尾部以及呆木处进行适灣加密.以表1中的船模为研究对象，采用3种不同 数量的网格模拟相同工况下的纯首摇运动，进行 网格收敛性的验证.网格数分别为1〇〇万，200万 和300万.图4给出了不同网格数计算得到的纵 向力X 、横向力Y 和转首力矩M 在一个周期内的 曲线.由图4可知，网格数量从100万增至200 万，计算结果变化明显；当网格数纛从200万增至 300万时，计算结果几乎不变.因此选择250万左 右的网格，既能保证收敛性，又能节省计算资源和 计算时间.表1供应船模参数---X (N)-100万网络一只N )-200万网络00万网络——y (N )-100万网络 ~1- Z (N)-200万网络 —八扣-300万网络 ...M (N .m )-100 万网络 +M (N .m )-200万网络 一 M (N .m)-300万网络时间f/s 图4不同网格数的计算结果对比数值计算方法图1供应船模型3.2 坐标系坐标系见图2，0点位于船舯；X 轴指向船首 为正;Y 轴指向左舷为正；Z 轴正向依据右手定则1)边界条件①人流与出流边界条件，人流 面设为速度人口；出流面设为自由流出口，权重为 1|②船体表面，在船体表面施加无滑移壁面条件;③自由面，考虑到供应船舶在定位过程中时低速• 72 •武汉理工大学学报（交通科学与工程版)2017年第41卷航行，忽略自由面兴波的影响，将自由面设为对称面.2)定义动网格(dynamic mesh )编写纯横 荡与纯首摇运动的用户自定义函数（U D F )，并使 之与F L U E N T 相关联；网格更新方法（mesh methods )选择网格光顺方法（smoothing )和动态 层方法（layering ).3)离散格式和求解算法非稳态流动；压强插值格式选用标准格式I 空间离散采用二阶迎风 差分格式；时间积分方案采用一阶隐式;求解算法 采用基于速度-压力耦合的S IM P L E C 算法.4试验设计及计算结果分析4.1供应船非线性模型的振幅范围计算工况见表2〜3.其中，= 和• L /L 7。

分别为纯横荡和纯首摇无因次化的速度幅值.表2纯横荡运动计算工况纵向运动速度队/(m *袭-1》振荡類率_/(ra d •_s—工）横向速度幅值的无因次t/横向位移幅值，/m横荡运动中的最大漂角iSmax/C")10. 20. 311, 30• 30. 4516. 70,430. 6523. 30. 60,40. 60.931. 00,8 1. 238. 71. 0 1. 545. 0表3纯首摇运动计算工况纵向运动速度振荡频率⑴首摇角速度幅竹摇运动中的/:Crad • s -”值的无因次/最大首向角平/〇0. 2480. 53180. 74250. 80. 40. 89301…04351. 3445图5〜6给出了不同运动幅值时的纯横荡与 纯首摇运动所受横向力Y 在一个周期的变化曲 线.图7〜8给出了纯横荡与纯首摇运动横向力Y 高阶量的实际值与E 阶拟合值的对比.从图中可 以得出以下结论.1) 随着横荡幅值和首摇角幅值的增大，横向力幅值明显增大，但横向力的相位几乎没有发生变化.2)当振荡幅度超过某一幅值，横向力已经不冉是一阶正(余）弦的形式，即运动幅值越大，非线 性表现的越强.3) 由图6可知，当初始首向角为8°时，纯首摇运动表现为一阶正余弦的形式；当初始首向角 为18°时，纯首摇运动已经不再是一阶正（余）弦的形式.4)由图7〜8可知，当横荡运动的横向幅值 超过Q . 9 m (运动过程中的最大漂角约为31°)，当 首摇运动的初始首向角超过30°时，横向水动力高 阶量的实际值与三阶拟合值间的误差越来越大.图575 m 供应船纯横荡运动横向力y随横荡幅值的变化时间z7s图6 75 m 供应船纯首摇运动横向力y随初始首向角的变化图7纯横荡运动y 方向高阶量的拟合值与实际值的对比图8纯首摇运动y 方向高阶量的拟合值与实际值的对比因此，对于75 m 供应船模，其非线性水动力模型的振幅范围约为漂角大于10°小于3〇\4.2流场分析图9〜10分别给出了 75 m 供应船首摇角幅 值为8°和25°的纯首摇运动一个周期内的船体表第1期赵小仨，等：非线性水动力导数的数值计算与研究• 73 •方法类似.根据纯横荡的运动规律，非线性水动力 数学模型可变形为：fY == Y.dsin (W ) +Y cl cos (cot) + Yr3 cos (3c 〇t)M = M 3 + Mt,*p +' M 爾齋3 + — X C l)uv =M si sin (cot ) ~h cos (cot ) ~h ^T r 3 cos (3ajt )(3)式中：L = — Y 具⑴，Yfl = Y 義 +1 Y _ ^，Yr 3 =+Y 卿W ，Md = —，Ma = [Mt, +(Y* — U+ — M zw v v 〇, = — M m m v 〇分别计算不同振动频率时船体所受的水动力X ，Y ，iV ，再用最小二乘法把水动力按上式拟合， 分离出系数Y sl ，^，，M sl ，M d ，M u ，进而得到相应的水动力导数.2)直航运动通过模拟0. 4 m /s 〜0.8 m/s 的直航运动，可求取纵向速度导数，X W W ，X _. 图11给出了不同航速时纵向力的变化曲线.用最 小二乘法对图11中的曲线进行拟合，再根据纵向 水动力的三阶非线性表达式，可以直接得出3个纵向水动力导数，结果见表4.图11直航运动纵向力的拟含曲线表4直航运动得到的纵向速度导数无因次量X XL X -数值X-0.42-2* 370:. 213)斜航运动通过模拟10°〜3Cf 不同漂角 下的斜航运动，可求取速度导数X 胃，，N t ,，N _.图12给出了不同漂角时纵向力、横向力和 力矩的变化曲线.与直航运动求取水动力导数的 方法相同，得出斜航运动对应的五个水动力导数， 结果见表5.纵向速f e /(m • s-1)0.10 0.150.20 0.250.30图12斜航运动中纵向力、横向力及首摇力矩的拟合曲线面压力云图.压力变化的总体趋势是:纵向来看， 从船首向船中压力逐渐递减，从船中向船尾压力 逐渐递增；横向来看，随着首向角转动位置与方向 的不同，最大压力的分布左右交替出现.通过对比 图9〜10可知，首摇幅值为25°时的压差比首摇幅 值为8°的压差要大得多.c) 5/8周期d) 7/8周期圈9 75 m 供应船首摇幅值为8°的船体表面压力分布图由图9可知，前半周期船体表面的压力分布 情况与后半周期是关于船体中纵剖面对称的.图 9b )表示首摇运动3/8周期时的压力分布图，它的 左舷压力分布恰好与图9c )中右舷压力分布相 同.出现该种现象的主要原因是:纯首摇运动的前 半周期与后半周期的转首角速度大小相同方向相 反，而且当首摇幅值为8°时，属于小振幅振动，符 合线性假定，因此由运动输出的力在前半周期与 后半周期是左右舷对称的.由图10可知，前半周期船体表面的压力分布 情况与后半周期已不再是关于船体中纵剖面对 称.其主要原因是：当首摇幅值为2^时，属于大 振幅振动，运动的输人与输出已不是线性关系. 4.3非线性水动力导数的计算结果以7S 米供应船为研究对象，根据4, 1节所得 结论，确定非线性水动力的计算工况.进而通过模 拟直航、大漂角斜航，以及频率为〇. 2〜0. 4 md/s 的大振幅纯纵荡、纯横荡、纯首摇和组合运动，计 算各运动模态下的水动力和水动力矩，然后求取 不同频率时非线性模型的水动力导数，最终得到 75 m 供应船在零频率附近的各水动力导数1)数据处理方案以纯横荡运动为例，介绍 本文的数据处理方案.其他运动形式的数据处理2 4 6 8 0N/X -R e 豸_•74•武汉理工大学学报（交通科学与工程版）2017年第41卷表5斜航运动得到的速度导数无因次量M；Y'—数值X10 20. 50-2. 95—1.63-7. 820. 274)纯纵荡运动通过模拟0. 2〜0. 5md/s不同频率下的纯纵荡运动，可求取纵向加速度导数u z°].根据数学模型中纵向力的表达式，用最小二乘方法分离出不同振荡频率时对应的系数xsl，又因为xs l=—x s M。