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

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

0 引言

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

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

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

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

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

1 理論基礎(chǔ)

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

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

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

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

航行器受到的作用力為

航行器受到的力矩為

2 數(shù)值計(jì)算的試驗(yàn)驗(yàn)證

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

圖1 迎浪零航速時(shí)單位波幅縱向力幅值Fig.1 Longitudinal forces in unit amplitude at zero velocity in head sea

圖2 迎浪零航速時(shí)單位波幅垂向力幅值Fig.2 Vertical forces in unit amplitude at zero velocity in head sea

圖3 迎浪零航速時(shí)單位波幅俯仰力矩幅值Fig.3 Pitching moments in unit amplitude at zero velocity in head sea

3 模型和結(jié)果分析

3.1 水動(dòng)力模型的建立與網(wǎng)格劃分

依照便攜式水下航行器實(shí)體在UG 中建立模型，模型主要參數(shù):長1 850 mm，直徑為200 mm，質(zhì)量為50 kg.將建立好的模型導(dǎo)入ANSYS-apdl，選擇SHELL181 單元建立面網(wǎng)格［18］，殼法向方向指向模型外部。由宏命令生成讀入文件，并對(duì)讀入文件進(jìn)行參數(shù)修改。

圖4 迎浪低航速時(shí)單位波幅縱向力幅值Fig.4 Longitudinal forces in unit amplitude at low velocity in head sea

圖5 迎浪低航速時(shí)單位波幅垂向力幅值Fig.5 Vertical forces in unit amplitude at low velocity in head sea

圖6 迎浪低航速時(shí)單位波幅俯仰力矩幅值Fig.6 Pitching moments in unit amplitude at low velocity in head sea

對(duì)水下航行器進(jìn)行粗、細(xì)、超細(xì)3 種網(wǎng)格劃分，為保證精度要求，3 種網(wǎng)格需滿足一個(gè)波長至少覆蓋7 個(gè)單元，相鄰單元面積比不得小于1/3.將3 種網(wǎng)格計(jì)算的水下航行器波浪力進(jìn)行比較，發(fā)現(xiàn)3 種網(wǎng)格計(jì)算的結(jié)果基本無差別，為保證計(jì)算精度，之后的計(jì)算均采用細(xì)網(wǎng)格。圖7為最終選擇的網(wǎng)格劃分示意圖。

圖7 無人水下航行器網(wǎng)格劃分Fig.7 The mesh of UUV

3.2 水動(dòng)力分析

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

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

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

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

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

圖9 不同水深下的側(cè)向力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ù)擬合

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

式中:A*為縱向力幅值;δ 為初相位。將靜止時(shí)與定深恒速時(shí)受到的波浪力進(jìn)行對(duì)比，得到:

圖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.對(duì)于有航速情況下波浪力和力矩的擬合，采用實(shí)時(shí)的航行器縱向位置來代替遭遇頻域，即縱向位置為x 的航行器在t 時(shí)刻受到的波浪力等同于初始位置的航行器在t'時(shí)刻受到的力。遭遇頻率不適用速度變化的情況，但對(duì)于航行器而言，在近水面航行時(shí)速度不可能會(huì)保持常值，所以用實(shí)時(shí)航行器的位置信息可以解決變航速問題。

圖14中給出了ω 分別為0.5 rad/s、1.0 rad/s、1.3 rad/s，水下深度h 為1 ～15 m 時(shí)所對(duì)應(yīng)的力的幅值A(chǔ)*.由最小二乘擬合得到3 條曲線的函數(shù)表達(dá)式為

由波頻對(duì)應(yīng)的波數(shù)k 分別為0.025 5、0.101 6、0.172 3，可得水下深度對(duì)幅值的影響為指數(shù)的。當(dāng)ω=1.0 rad/s時(shí)，A*=56.48Ae-kh.

圖15中給出，ω 為0.3 ～1.6 rad/s，h =4 m 時(shí)所對(duì)應(yīng)的A*/Ae-kh，擬合得到A*/Ae-kh=56.48ω1.92.由于水下航行器的對(duì)稱性，認(rèn)為正負(fù)側(cè)滑角對(duì)波浪力的影響相同，數(shù)值結(jié)果也證明了這一點(diǎn)，考慮到水下航行器的常用姿態(tài)，表1給出了h =2 m 時(shí)側(cè)滑角β 為0° ～15°所對(duì)應(yīng)的A*/Ae-khω1.92.

圖14 水下深度對(duì)于力的幅值A(chǔ)* 影響Fig.14 Impact of depth on force amplitude

圖15 頻率對(duì)力的幅值影響Fig.15 Impact of frequency on force amplitude

圖16 俯仰角對(duì)初相位的影響Fig.16 Impact of pitch angle on initial phase

由圖16可以看出俯仰角θ 對(duì)初相位的影響為線性的，由最小二乘擬合得

由表1計(jì)算發(fā)現(xiàn)帶β 角度側(cè)滑角的幅值與零側(cè)滑角幅值關(guān)系為

考慮水下航行器的常用姿態(tài)，研究俯仰角的范圍為-12° ～12°.表2給出h =2 m 時(shí)，θ 為-12° ～12°所對(duì)應(yīng)的A*/Ae-khω1.92.

表1 側(cè)滑角對(duì)于幅值的影響Tab.1 Impact of sideslip angle on force amplitude

表2 俯仰角對(duì)于幅值的影響Tab.2 Impact of pitch angle on force amplitude

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

根據(jù)對(duì)不同工況下水下航行器波浪力的時(shí)域計(jì)算，對(duì)于重要因素進(jìn)行函數(shù)擬合，得到水下航行器縱向波浪力、側(cè)向波浪力、垂向浪力和俯仰力矩力矩的擬合函數(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ù)的準(zhǔn)確性驗(yàn)證

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

工況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校驗(yàn)Fig.17 Verification of surge between numerical simulation and fitting function results in two conditions

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

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

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

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

4 結(jié)論

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

為了驗(yàn)證擬合函數(shù)的準(zhǔn)確性，采用任意兩種工況對(duì)擬合函數(shù)進(jìn)行校驗(yàn)。從校驗(yàn)結(jié)果來看，擬合函數(shù)的結(jié)果與數(shù)值解相當(dāng)接近，其中縱向波浪力和俯仰力矩在100 步的均方誤差均小于1%，垂向波浪力的均方誤差也不超過2%.從工況2 的擬合函數(shù)解與數(shù)值解的對(duì)比可以看出，可以采用實(shí)時(shí)的航行器縱向位置來代替有航速情況下的遭遇頻率。采用實(shí)時(shí)的縱向位置代替遭遇頻率不僅適用于航行器恒速時(shí)的波浪力預(yù)報(bào)，也適用于航行器變速情況下的波浪力預(yù)報(bào)。任意兩種工況的校驗(yàn)說明通過此方法得到預(yù)報(bào)的波浪載荷可以為評(píng)估水下航行器在海洋環(huán)境下的操縱性提供具有一定參考價(jià)值的數(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òng)力計(jì)算［D］.哈爾濱:哈爾濱工程大學(xué)，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ù)及其導(dǎo)數(shù)的數(shù)值逼近［J］.水動(dòng)力學(xué)研究與進(jìn)展: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］ 繆國平，劉應(yīng)中，楊勤正，等.三維移動(dòng)脈沖源的Michell 型表達(dá)式［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］ 李遠(yuǎn)林.波浪理論及波浪載荷［M］.廣州:華南理工大學(xué)出版社，1994.LI Yuan-lin.Wave theory and wave loads［M］.Guangzhou:South China University of Technology Press，1994.(in Chinese)

［14］ 戴遺山，段文洋.船舶在波浪中運(yùn)動(dòng)的勢(shì)流理論［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］ 竺艷蓉.海洋工程波浪力學(xué)［M］.天津:天津大學(xué)出版社，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］ 梅強(qiáng)中.水波動(dòng)力學(xué)［M］.北京:科學(xué)出版社，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.