潘光，劉亞楠，杜曉旭，Saeed Akram Malik

(西北工業(yè)大學 航海學院，陜西 西安710129)

0 引言

20世紀90年代以來，無人水下航行器以其能夠提供最為直接、有效的海洋信息，且開發(fā)成本較低等獨特優(yōu)勢，得到迅速的發(fā)展，被用于海洋環(huán)境監(jiān)測與生物資源調查等方面。無人水下航行器在海洋中發(fā)射和航行，其運動不可避免地受到波浪的影響，尤其是水下航行器在近水面航行跟蹤目標或獲取信號時，波浪運動的壓力場使水下航行器受到附加的壓力和慣性力的影響，會引起水下航行器的顛簸與搖擺，從而影響到水下航行器的操縱性。便攜式水下航行器具有體積小、質量輕、附體較大等特點，為了準確地預報其在近水面航行或作業(yè)時受到的波浪力以提高它的操控性，需要考慮附體(尤其是鰭舵)的影響。

當結構物長度遠遠大于波浪的幅值時，結構物受到的波浪力主要是慣性力，粘性力很小可以忽略，因此在水動力學領域，廣泛基于勢流理論計算結構物在波浪中運動的載荷響應。近年來，潛體在波浪中運動的理論和計算方法從二維切片理論發(fā)展至三維邊界元法、從零航速發(fā)展至有航速，從頻域法發(fā)展至時域法，并取得了一系列重要的進展，這些方法都是應用勢流理論，并對自由面條件進行了線性化近似。

切片理論［1］和它的改良方法，用于預報水下航行器以中低航速航行時受到波浪力的線性預報。切片理論因具有易于理解和計算快捷等優(yōu)點而得到廣泛地應用，但卻未能考慮切片間的流體干擾作用，并且對于模型的建立過于簡化，不能很好地考慮附體的影響。

自從Hess 等［2］開創(chuàng)性地提出以分布源和分布偶極的三維面元來求解其運動理論以來，各國學者紛紛對其進行了研究使其逐漸發(fā)展成為勢流問題求解中應用最廣的數(shù)值方法。邊界元法引入格林函數(shù)，將流體域內的體積積分轉換到流體域邊界上的面積分。目前，邊界元法主要有常值面元法、高階面元法［3-5］和無網格法［6］。頻域無航速無限水深自由面格林函數(shù)的數(shù)值計算方法［7］己經比較成熟，對于頻域有航速問題，很多學者進行了大量的研究［8-10］。為了避免自由面格林函數(shù)計算的復雜性，簡單格林函數(shù)在水動力計算中，尤其是時域問題計算中迅速發(fā)展起來。Cummins［11］提出脈沖響應的概念，把任一結構物運動的時間歷經看做一系列瞬時的脈沖運動所組成，同時波浪力也看作是一系列脈沖響應的線性迭加，從而建立起運動微分方程。

基于勢流理論，文中運用三維頻域格林函數(shù)法計算便攜式水下航行器近水面航行作業(yè)時受到的波浪力，從而得到便攜式水下航行器在不同情況下的傳遞函數(shù)、附加質量和阻尼系數(shù)，在此基礎上對便攜式水下航行器進行時域分析。建立與文獻［12］中試驗模型吻合的數(shù)學模型并進行數(shù)值計算，將數(shù)值計算結果與試驗結果進行對比;計算入射波不同參數(shù)和航行器不同姿態(tài)組合下的航行器受到的波浪力。分析不同因素對波浪力和力矩的影響，分析發(fā)現(xiàn)可以采用實時監(jiān)測的航行器縱向位置來代替有航速情況下的遭遇頻率;考慮各因素影響并對波浪力的數(shù)值計算結果進行函數(shù)擬合，將擬合函數(shù)與數(shù)值計算結果進行比較，二者差別較小，為能夠準確地預報水下航行器的運動軌跡帶來很大的方便。

1 理論基礎

基于勢流理論［13-17］，假設流體為理想流體，不計流體自由表面上的張力，流動是無旋的。Froude-Krylov 力是無擾動的入射波浪引起的壓力，衍射力是靜止結構的存在影響了波浪密度分布由壓差引起的壓力，輻射力是結構的運動或振動激起的波浪產生的波浪力。航行器表面的法向指向航行器內部，假設波浪是微幅的。用φ(p，t)表示時刻t 點p 處的流體總速度勢，由定常和非定常兩部分疊加組成。定常部分為結構物在靜水中行進達到定常狀態(tài)后的速度勢φS，不考慮波浪影響，其余的非定常部分有入射波勢φI和擾動勢φB，其中擾動勢φB包括衍射波勢和輻射波勢。

式中:φB(p，t)=為航行器k 態(tài)運動的輻射勢，φ7(p，t)是衍射勢。取大地坐標系Oxyz，xy 平面為靜水面，z 軸垂直向上。衍射勢和輻射勢的定解條件為

此外還有遠方條件和底部條件，其中ηj(j =1，2，3，4，5，6)表示結構物j 態(tài)位移，w 為誘導速度，s(t)為航行器表面。使用源分布求解擾動勢?？紤]二階Stokes 波為入射波，則入射波勢為

式中:ω2=gktanh(kH)，H 為水深;A 為波幅;ω 為入射波頻率;θ 為入射波角度;k 為波數(shù)。將坐標原點下移到航行器內部的參考點，用隨航行器運動的平移坐標系Oxyz.在平移坐標系中，航行器表面是固定不變的，與時間t 無關，用s0表示。

航行器受到的作用力為

航行器受到的力矩為

2 數(shù)值計算的試驗驗證

為了驗證數(shù)值計算的準確性，將數(shù)值計算結果與文獻［12］的試驗數(shù)據(jù)進行對比。建立與試驗模型相同的數(shù)學模型進行數(shù)值計算，在零航速和低航速兩種工況下對比單位波幅下力的幅值隨入射波頻率變化曲線;模型的主要參數(shù):長1 825 mm，直徑為127 mm，質心離水面距離為z =0.379 m，有航速情況下航速為v =0.489 m/s.圖1～圖7中‘×’，‘·’，‘+’為模型在逐漸增長的波幅下的試驗數(shù)據(jù)，實線為數(shù)值計算結果。波幅:A1= 0.014 m，A2=0.028 m，A3=0.037 m，有航速情況下采用兩種波幅進行試驗。對比結果表明:數(shù)值計算結果和試驗數(shù)據(jù)在趨勢上是一致性，但在具體數(shù)值上有一定的差距。差距一部分是由于數(shù)值計算模型造成的，一部分是由于試驗造成的，差距在可以容許的范圍內，且不影響整體趨勢的判斷。

圖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 中建立模型，模型主要參數(shù):長1 850 mm，直徑為200 mm，質量為50 kg.將建立好的模型導入ANSYS-apdl，選擇SHELL181 單元建立面網格［18］，殼法向方向指向模型外部。由宏命令生成讀入文件，并對讀入文件進行參數(shù)修改。

圖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 種網格計算的水下航行器波浪力進行比較，發(fā)現(xiàn)3 種網格計算的結果基本無差別，為保證計算精度，之后的計算均采用細網格。圖7為最終選擇的網格劃分示意圖。

圖7 無人水下航行器網格劃分Fig.7 The mesh of UUV

3.2 水動力分析

應用AQWA-LINE 對水下航行器進行頻域水動力分析，采用線性波作為入射波，研究在不同頻率和不同深度下便攜式水下航行器受到的波浪力。

由圖8～圖10可以看出，3 個自由度上的波浪力隨著入射波頻率的變大趨勢均是先變大達到極值后變小;波浪力隨著航深增加迅速減小，使波浪力達到極值的入射波頻率隨航深增加逐漸變小，側向力較縱向力和垂向力較小，所以在時域分析和波浪力的函數(shù)擬合階段忽略對側向力的分析。

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

由圖11～圖13可以看出，橫滾和偏航力矩相比俯仰力矩小很多，所以在時域分析和波浪力矩的函數(shù)擬合階段忽略對橫滾力矩和偏航力矩的分析。

運用AQWA-NAUT 對水下航行器進行時域水動力分析。采用二階Stokes 波作為入射波，假設水下航行器有良好的操控性，研究水下航行器定深恒速航行時受到的波浪力。水下航行器受到的波浪力為正弦函數(shù)、余弦函數(shù)或它們的線性組合，不同工況下的時域結果為水下航行器波浪力擬合提供數(shù)據(jù)。

圖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 參數(shù)影響分析與函數(shù)擬合

通過對水下航行器的頻域和時域分析結果研究，開始研究不同因素對于便攜式水下航行器波浪力的影響。對于固定于水下某點的物體所受的波浪力為正弦曲線，設為

式中: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 條曲線的函數(shù)表達式為

由波頻對應的波數(shù)k 分別為0.025 5、0.101 6、0.172 3，可得水下深度對幅值的影響為指數(shù)的。當ω=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.由于水下航行器的對稱性，認為正負側滑角對波浪力的影響相同，數(shù)值結果也證明了這一點，考慮到水下航行器的常用姿態(tài)，表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計算發(fā)現(xiàn)帶β 角度側滑角的幅值與零側滑角幅值關系為

考慮水下航行器的常用姿態(tài)，研究俯仰角的范圍為-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

在研究俯仰角對于力的幅值變化影響時發(fā)現(xiàn)，當水下航行器抬頭時，幅值隨著仰角的增大而增大;當水下航行器低頭時，幅值隨著俯仰角的俯角變大先變小一直增大，在θ = ±2°左右時幅值達到最小。

根據(jù)對不同工況下水下航行器波浪力的時域計算，對于重要因素進行函數(shù)擬合，得到水下航行器縱向波浪力、側向波浪力、垂向浪力和俯仰力矩力矩的擬合函數(shù)。

縱向力的擬合函數(shù)為

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

垂向力的擬合函數(shù)為

式中: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).

俯仰力矩的擬合函數(shù)為

式中: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 擬合函數(shù)的準確性驗證

選取擬合函數(shù)自變量參數(shù)定義域內任意兩種工況，驗證擬合函數(shù)的精確定。圖17～圖19分別為兩種不同工況下航行器受到的縱向力、垂向力和俯仰力矩的數(shù)值解與函數(shù)擬合值的對比。

工況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 兩種工況下擬合函數(shù)Fx校驗Fig.17 Verification of surge between numerical simulation and fitting function results in two conditions

圖18 兩種工況下擬合函數(shù)Fz校驗Fig.18 Verification of heave between numerical simulation and fitting function results in two conditions

圖19 兩種工況下擬合函數(shù)My校驗Fig.19 Verification of pitch between numerical simulation and fitting function results in two conditions

表3為兩種工況下縱向波浪力、垂向波浪力、俯仰力矩擬合函數(shù)100 步(5 s)的均方誤差(MSE)，步長為0.05 s.

表3 擬合函數(shù)的均方誤差Tab.3 Mean square error of fitting functions

4 結論

文中對考慮附體(鰭舵)的便攜式水下航行器進行頻域和時域水動力分析?；趧萘骼碚?，采用面元法計算得到了水下航行器近水面航行受到的波浪載荷。分析得到波頻、航行深度、航行位置、航行姿態(tài)對于水下航行器所受波浪力的影響?？紤]各參數(shù)的影響對航行器受到的波浪力和力矩進行函數(shù)擬合。

為了驗證擬合函數(shù)的準確性，采用任意兩種工況對擬合函數(shù)進行校驗。從校驗結果來看，擬合函數(shù)的結果與數(shù)值解相當接近，其中縱向波浪力和俯仰力矩在100 步的均方誤差均小于1%，垂向波浪力的均方誤差也不超過2%.從工況2 的擬合函數(shù)解與數(shù)值解的對比可以看出，可以采用實時的航行器縱向位置來代替有航速情況下的遭遇頻率。采用實時的縱向位置代替遭遇頻率不僅適用于航行器恒速時的波浪力預報，也適用于航行器變速情況下的波浪力預報。任意兩種工況的校驗說明通過此方法得到預報的波浪載荷可以為評估水下航行器在海洋環(huán)境下的操縱性提供具有一定參考價值的數(shù)據(jù)。

References)

［1］ Salvesen N，Tuck E O，F(xiàn)altinsen 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 函數(shù)及其導數(shù)的數(shù)值逼近［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］.北京:國防工業(yè)出版社，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.