便携式水下航行器近水面航行波浪力计算分析与预报

潘光，刘亚楠，杜晓旭，Saeed Akram Malik

(西北工业大学 航海学院，陕西 西安710129)

0 引言

20世纪90年代以来，无人水下航行器以其能够提供最为直接、有效的海洋信息，且开发成本较低等独特优势，得到迅速的发展，被用于海洋环境监测与生物资源调查等方面。无人水下航行器在海洋中发射和航行，其运动不可避免地受到波浪的影响，尤其是水下航行器在近水面航行跟踪目标或获取信号时，波浪运动的压力场使水下航行器受到附加的压力和惯性力的影响，会引起水下航行器的颠簸与摇摆，从而影响到水下航行器的操纵性。便携式水下航行器具有体积小、质量轻、附体较大等特点，为了准确地预报其在近水面航行或作业时受到的波浪力以提高它的操控性，需要考虑附体(尤其是鳍舵)的影响。

当结构物长度远远大于波浪的幅值时，结构物受到的波浪力主要是惯性力，粘性力很小可以忽略，因此在水动力学领域，广泛基于势流理论计算结构物在波浪中运动的载荷响应。近年来，潜体在波浪中运动的理论和计算方法从二维切片理论发展至三维边界元法、从零航速发展至有航速，从频域法发展至时域法，并取得了一系列重要的进展，这些方法都是应用势流理论，并对自由面条件进行了线性化近似。

切片理论［1］和它的改良方法，用于预报水下航行器以中低航速航行时受到波浪力的线性预报。切片理论因具有易于理解和计算快捷等优点而得到广泛地应用，但却未能考虑切片间的流体干扰作用，并且对于模型的建立过于简化，不能很好地考虑附体的影响。

自从Hess 等［2］开创性地提出以分布源和分布偶极的三维面元来求解其运动理论以来，各国学者纷纷对其进行了研究使其逐渐发展成为势流问题求解中应用最广的数值方法。边界元法引入格林函数，将流体域内的体积积分转换到流体域边界上的面积分。目前，边界元法主要有常值面元法、高阶面元法［3-5］和无网格法［6］。频域无航速无限水深自由面格林函数的数值计算方法［7］己经比较成熟，对于频域有航速问题，很多学者进行了大量的研究［8-10］。为了避免自由面格林函数计算的复杂性，简单格林函数在水动力计算中，尤其是时域问题计算中迅速发展起来。Cummins［11］提出脉冲响应的概念，把任一结构物运动的时间历经看做一系列瞬时的脉冲运动所组成，同时波浪力也看作是一系列脉冲响应的线性迭加，从而建立起运动微分方程。

基于势流理论，文中运用三维频域格林函数法计算便携式水下航行器近水面航行作业时受到的波浪力，从而得到便携式水下航行器在不同情况下的传递函数、附加质量和阻尼系数，在此基础上对便携式水下航行器进行时域分析。建立与文献［12］中试验模型吻合的数学模型并进行数值计算，将数值计算结果与试验结果进行对比;计算入射波不同参数和航行器不同姿态组合下的航行器受到的波浪力。分析不同因素对波浪力和力矩的影响，分析发现可以采用实时监测的航行器纵向位置来代替有航速情况下的遭遇频率;考虑各因素影响并对波浪力的数值计算结果进行函数拟合，将拟合函数与数值计算结果进行比较，二者差别较小，为能够准确地预报水下航行器的运动轨迹带来很大的方便。

1 理论基础

基于势流理论［13-17］，假设流体为理想流体，不计流体自由表面上的张力，流动是无旋的。Froude-Krylov 力是无扰动的入射波浪引起的压力，衍射力是静止结构的存在影响了波浪密度分布由压差引起的压力，辐射力是结构的运动或振动激起的波浪产生的波浪力。航行器表面的法向指向航行器内部，假设波浪是微幅的。用φ(p，t)表示时刻t 点p 处的流体总速度势，由定常和非定常两部分叠加组成。定常部分为结构物在静水中行进达到定常状态后的速度势φS，不考虑波浪影响，其余的非定常部分有入射波势φI和扰动势φB，其中扰动势φB包括衍射波势和辐射波势。

式中:φB(p，t)=为航行器k 态运动的辐射势，φ7(p，t)是衍射势。取大地坐标系Oxyz，xy 平面为静水面，z 轴垂直向上。衍射势和辐射势的定解条件为

此外还有远方条件和底部条件，其中ηj(j =1，2，3，4，5，6)表示结构物j 态位移，w 为诱导速度，s(t)为航行器表面。使用源分布求解扰动势。考虑二阶Stokes 波为入射波，则入射波势为

式中:ω2=gktanh(kH)，H 为水深;A 为波幅;ω 为入射波频率;θ 为入射波角度;k 为波数。将坐标原点下移到航行器内部的参考点，用随航行器运动的平移坐标系Oxyz.在平移坐标系中，航行器表面是固定不变的，与时间t 无关，用s0表示。

航行器受到的作用力为

航行器受到的力矩为

2 数值计算的试验验证

为了验证数值计算的准确性，将数值计算结果与文献［12］的试验数据进行对比。建立与试验模型相同的数学模型进行数值计算，在零航速和低航速两种工况下对比单位波幅下力的幅值随入射波频率变化曲线;模型的主要参数:长1 825 mm，直径为127 mm，质心离水面距离为z =0.379 m，有航速情况下航速为v =0.489 m/s.图1～图7中‘×’，‘·’，‘+’为模型在逐渐增长的波幅下的试验数据，实线为数值计算结果。波幅:A1= 0.014 m，A2=0.028 m，A3=0.037 m，有航速情况下采用两种波幅进行试验。对比结果表明:数值计算结果和试验数据在趋势上是一致性，但在具体数值上有一定的差距。差距一部分是由于数值计算模型造成的，一部分是由于试验造成的，差距在可以容许的范围内，且不影响整体趋势的判断。

图1 迎浪零航速时单位波幅纵向力幅值Fig.1 Longitudinal forces in unit amplitude at zero velocity in head sea

图2 迎浪零航速时单位波幅垂向力幅值Fig.2 Vertical forces in unit amplitude at zero velocity in head sea

图3 迎浪零航速时单位波幅俯仰力矩幅值Fig.3 Pitching moments in unit amplitude at zero velocity in head sea

3 模型和结果分析

3.1 水动力模型的建立与网格划分

依照便携式水下航行器实体在UG 中建立模型，模型主要参数:长1 850 mm，直径为200 mm，质量为50 kg.将建立好的模型导入ANSYS-apdl，选择SHELL181 单元建立面网格［18］，壳法向方向指向模型外部。由宏命令生成读入文件，并对读入文件进行参数修改。

图4 迎浪低航速时单位波幅纵向力幅值Fig.4 Longitudinal forces in unit amplitude at low velocity in head sea

图5 迎浪低航速时单位波幅垂向力幅值Fig.5 Vertical forces in unit amplitude at low velocity in head sea

图6 迎浪低航速时单位波幅俯仰力矩幅值Fig.6 Pitching moments in unit amplitude at low velocity in head sea

对水下航行器进行粗、细、超细3 种网格划分，为保证精度要求，3 种网格需满足一个波长至少覆盖7 个单元，相邻单元面积比不得小于1/3.将3 种网格计算的水下航行器波浪力进行比较，发现3 种网格计算的结果基本无差别，为保证计算精度，之后的计算均采用细网格。图7为最终选择的网格划分示意图。

图7 无人水下航行器网格划分Fig.7 The mesh of UUV

3.2 水动力分析

应用AQWA-LINE 对水下航行器进行频域水动力分析，采用线性波作为入射波，研究在不同频率和不同深度下便携式水下航行器受到的波浪力。

由图8～图10可以看出，3 个自由度上的波浪力随着入射波频率的变大趋势均是先变大达到极值后变小;波浪力随着航深增加迅速减小，使波浪力达到极值的入射波频率随航深增加逐渐变小，侧向力较纵向力和垂向力较小，所以在时域分析和波浪力的函数拟合阶段忽略对侧向力的分析。

图8 不同水深下的纵向力Fig.8 Longitudinal forces in different depths

由图11～图13可以看出，横滚和偏航力矩相比俯仰力矩小很多，所以在时域分析和波浪力矩的函数拟合阶段忽略对横滚力矩和偏航力矩的分析。

运用AQWA-NAUT 对水下航行器进行时域水动力分析。采用二阶Stokes 波作为入射波，假设水下航行器有良好的操控性，研究水下航行器定深恒速航行时受到的波浪力。水下航行器受到的波浪力为正弦函数、余弦函数或它们的线性组合，不同工况下的时域结果为水下航行器波浪力拟合提供数据。

图9 不同水深下的侧向力Fig.9 Yawing forces in different depths

图10 不同水深下的垂向力Fig.10 Vertical forces in different depths

图11 不同水深下的横滚力矩Fig.11 Rolling moments in different depths

3.3 参数影响分析与函数拟合

通过对水下航行器的频域和时域分析结果研究，开始研究不同因素对于便携式水下航行器波浪力的影响。对于固定于水下某点的物体所受的波浪力为正弦曲线，设为

式中:A*为纵向力幅值;δ 为初相位。将静止时与定深恒速时受到的波浪力进行对比，得到:

图12 不同水深下的俯仰力矩Fig.12 Pitching moments in different depths

图13 不同水深下的偏航力矩Fig.13 Yawing moments in different depths

式中:ωe=ω +ω2u/g;t' =t +x/c，c =ω/k.对于有航速情况下波浪力和力矩的拟合，采用实时的航行器纵向位置来代替遭遇频域，即纵向位置为x 的航行器在t 时刻受到的波浪力等同于初始位置的航行器在t'时刻受到的力。遭遇频率不适用速度变化的情况，但对于航行器而言，在近水面航行时速度不可能会保持常值，所以用实时航行器的位置信息可以解决变航速问题。

图14中给出了ω 分别为0.5 rad/s、1.0 rad/s、1.3 rad/s，水下深度h 为1 ～15 m 时所对应的力的幅值A*.由最小二乘拟合得到3 条曲线的函数表达式为

由波频对应的波数k 分别为0.025 5、0.101 6、0.172 3，可得水下深度对幅值的影响为指数的。当ω=1.0 rad/s时，A*=56.48Ae-kh.

图15中给出，ω 为0.3 ～1.6 rad/s，h =4 m 时所对应的A*/Ae-kh，拟合得到A*/Ae-kh=56.48ω1.92.由于水下航行器的对称性，认为正负侧滑角对波浪力的影响相同，数值结果也证明了这一点，考虑到水下航行器的常用姿态，表1给出了h =2 m 时侧滑角β 为0° ～15°所对应的A*/Ae-khω1.92.

图14 水下深度对于力的幅值A* 影响Fig.14 Impact of depth on force amplitude

图15 频率对力的幅值影响Fig.15 Impact of frequency on force amplitude

图16 俯仰角对初相位的影响Fig.16 Impact of pitch angle on initial phase

由图16可以看出俯仰角θ 对初相位的影响为线性的，由最小二乘拟合得

由表1计算发现带β 角度侧滑角的幅值与零侧滑角幅值关系为

考虑水下航行器的常用姿态，研究俯仰角的范围为-12° ～12°.表2给出h =2 m 时，θ 为-12° ～12°所对应的A*/Ae-khω1.92.

表1 侧滑角对于幅值的影响Tab.1 Impact of sideslip angle on force amplitude

表2 俯仰角对于幅值的影响Tab.2 Impact of pitch angle on force amplitude

在研究俯仰角对于力的幅值变化影响时发现，当水下航行器抬头时，幅值随着仰角的增大而增大;当水下航行器低头时，幅值随着俯仰角的俯角变大先变小一直增大，在θ = ±2°左右时幅值达到最小。

根据对不同工况下水下航行器波浪力的时域计算，对于重要因素进行函数拟合，得到水下航行器纵向波浪力、侧向波浪力、垂向浪力和俯仰力矩力矩的拟合函数。

纵向力的拟合函数为

式中:K = e-khω1.92cos β;δ(θ)=1.014θ;α1(θ)=89.22θ2+ 3.291θ + 56.48;α2(θ)= 86.31θ2-1.882θ+56.48.

垂向力的拟合函数为

式中:a0= 5.529a2.014ω4.007e-2kh+ 1.55;a1= A·e-khω1.988(0.935 1θ2+0.960 9θ -111.5);b1=(A·e-khω4.035+0.08β1.45)(-21.89θ-0.192 3).

俯仰力矩的拟合函数为

式中:a2=Ae-khω1.964(2.684θ2+5.399θ -4.736);b2=Ae-khω3.99cos β(5.381θ2-4.648θ-2.727).

3.4 拟合函数的准确性验证

选取拟合函数自变量参数定义域内任意两种工况，验证拟合函数的精确定。图17～图19分别为两种不同工况下航行器受到的纵向力、垂向力和俯仰力矩的数值解与函数拟合值的对比。

工况1:A =0.5 m，ω =1.1 rad/s，β =8°，θ =8°，h=1.5 m，k=0.123 4;

工况2:A =0.4 m，ω =1.2 rad/s，β =8°，θ =-5°，h=7 m，k=0.146 8，u0=5 m/s.

图17 两种工况下拟合函数Fx校验Fig.17 Verification of surge between numerical simulation and fitting function results in two conditions

图18 两种工况下拟合函数Fz校验Fig.18 Verification of heave between numerical simulation and fitting function results in two conditions

图19 两种工况下拟合函数My校验Fig.19 Verification of pitch between numerical simulation and fitting function results in two conditions

表3为两种工况下纵向波浪力、垂向波浪力、俯仰力矩拟合函数100 步(5 s)的均方误差(MSE)，步长为0.05 s.

表3 拟合函数的均方误差Tab.3 Mean square error of fitting functions

4 结论

文中对考虑附体(鳍舵)的便携式水下航行器进行频域和时域水动力分析。基于势流理论，采用面元法计算得到了水下航行器近水面航行受到的波浪载荷。分析得到波频、航行深度、航行位置、航行姿态对于水下航行器所受波浪力的影响。考虑各参数的影响对航行器受到的波浪力和力矩进行函数拟合。

为了验证拟合函数的准确性，采用任意两种工况对拟合函数进行校验。从校验结果来看，拟合函数的结果与数值解相当接近，其中纵向波浪力和俯仰力矩在100 步的均方误差均小于1%，垂向波浪力的均方误差也不超过2%.从工况2 的拟合函数解与数值解的对比可以看出，可以采用实时的航行器纵向位置来代替有航速情况下的遭遇频率。采用实时的纵向位置代替遭遇频率不仅适用于航行器恒速时的波浪力预报，也适用于航行器变速情况下的波浪力预报。任意两种工况的校验说明通过此方法得到预报的波浪载荷可以为评估水下航行器在海洋环境下的操纵性提供具有一定参考价值的数据。

References)

［1］ Salvesen N，Tuck E O，Faltinsen O.Ship motions and sea loads［M］.Oslo:Det norske Veritas，1971.

［2］ Hess J L，Smith A M.Calculation of non-lifting potential flow about arbitrary three bodies，AD0282255［R］.Long Beach:Douglas Aircraft Company，1962.

［3］ Tong K C.A 3D higher order boundary element method for wave forces on offshore structure［C］∥8th Offshore Mechanics ＆ Arctic Engng Conference.London:British Maritime Technology，1989:17 -27.

［4］ Teng B，Taylor R E.New higher-order boundary element methods for wave diffraction/radiation［J］.Applied Ocean Research，1995，17(2):71 -77.

［5］ 刘明日.基于B 样条面元法的浮体二阶水动力计算［D］.哈尔滨:哈尔滨工程大学，2009.LIU Ming-ri.Calculating second order hydrodynamics of floating structure based on B-spline panel method［D］.Harbin:Harbin Engineering University，2009.(in Chinese)

［6］ Qiu W.Apanel-free method for time-domain analysis of floating bodies in waves［D］.Halifax:Dalhousie University，2001.

［7］ 王如森.三维自由面Green 函数及其导数的数值逼近［J］.水动力学研究与进展:A 辑，1992，7(3):277 -286.WANG Ru-sen.3-D free surface green function and numerical approximation of its derivative［J］.Chinese Journal of Hydrodynamics:Ser A，1992，7(3):277 -286.(in Chinese)

［8］ 缪国平，刘应中，杨勤正，等.三维移动脉冲源的Michell 型表达式［J］.中国造船，1995，131(4):1 -11.MIU Guo-ping，LIU Ying-zhong，YANG Qin-zheng，et al.3-D Michell expression of impulse source［J］.Shipbuilding of China，1995，131(4):1 -11.(in Chinese)

［9］ Inglis R B，Price W G.Calculation of the velocity potential of a translating［J］.Transaction，RINA，1981，123:162 -175.

［10］ Fang M C，Chang P E，Luo J H.Wave effects on ascending and descending motions of the autonomous underwater vehicle［J］.Ocean Engineering，2006，33(14):1972 -1999.

［11］ Cummins W E.The impulse response function and ship motions，AD0288277［R］.Washington，DC:David Taylor Model Basin，1962:101 -109.

［12］ Willy C J.Attitude control of an underwater vehicle subjected to waves［D］.Cambridge:Massachusetts Institute of Technology，1994.

［13］ 李远林.波浪理论及波浪载荷［M］.广州:华南理工大学出版社，1994.LI Yuan-lin.Wave theory and wave loads［M］.Guangzhou:South China University of Technology Press，1994.(in Chinese)

［14］ 戴遗山，段文洋.船舶在波浪中运动的势流理论［M］.北京:国防工业出版社，2008.DAI Yi-shan，DUAN Wen-yang.Potential flow theory of ship motions in waves［M］.Beijing:National Defense Industry Press，2008.(in Chinese)

［15］ 竺艳蓉.海洋工程波浪力学［M］.天津:天津大学出版社，1991.ZHU Yan-rong.Wave forces in ocean engineering［M］.Tianjin:Tianjin Univeristy Press，1991.(in Chinese)

［16］ Faltinsen O.Sea loads on ships and offshore structures［M］.Cambridge:Cambridge University Press，1993.

［17］ 梅强中.水波动力学［M］.北京:科学出版社，1984.MEI Qiang-zhong.Wave dynamics［M］.Beijing:Sciences Press，1984.(in Chinese)

［18］ ZHANG Z.Development of adaptive samping power take-off control for a three-body wave energy converter with numerical modeling and validation［D］.Corvallis:Oregon State University，2011.