楊鯤， 盧倪斌， 隋海琛， 王磊峰， 李曄

(1.交通運(yùn)輸部 天津水運(yùn)工程科學(xué)研究所，天津 300456； 2.天津水運(yùn)工程勘察設(shè)計院 天津市水運(yùn)工程測繪技術(shù)重點實驗室，天津 300456； 3.哈爾濱工程大學(xué) 水下機(jī)器人技術(shù)重點實驗室，黑龍江 哈爾濱 150001)

波浪滑翔器作為新型無人海洋探測平臺，從波浪能和太陽能中分別獲得前進(jìn)的動力和電力，克服了大范圍、長航時的難題，具有續(xù)航能力強(qiáng)、自主控制、綠色環(huán)保、價格經(jīng)濟(jì)等突出特點[1]。它能長期、自主地執(zhí)行監(jiān)測環(huán)境、調(diào)查水文、預(yù)報氣象、追蹤生物、預(yù)警危害、中繼通訊等任務(wù)。這些特點也使得波浪滑翔器在軍事和民用范圍內(nèi)皆具有廣泛的應(yīng)用前景[2]。波浪滑翔器最先由美國的Roger Hine于2005年研制[3]。2年后，他成立了Liquid Robotics公司并進(jìn)行更加深入的研究[4]。后來，因為其優(yōu)點多，并且應(yīng)用需求與發(fā)展空間大，波浪滑翔器受到世界眾多學(xué)者的關(guān)注。然而，人們對其的研究主要集中于工程應(yīng)用，對其動力學(xué)模型的理論研究仍然很少。對于任何海洋航行器而言，操縱性能研究是必不可少且極其重要的。建立合適的動力學(xué)模型是開展運(yùn)動機(jī)理分析、運(yùn)動預(yù)報、運(yùn)動控制等研究的基礎(chǔ)，而一個優(yōu)秀的動力學(xué)模型往往是簡化研究難度的關(guān)鍵。

Kraus等[5]給出了波浪滑翔器沿前進(jìn)方向和升沉方向的二維模型。齊占峰等[6]利用Kane方程建立了波浪滑翔器在垂直面的多剛體動力學(xué)模型。在此基礎(chǔ)上，周春琳等[7]將浮體的垂直運(yùn)動加入到動力學(xué)模型中，考慮推力和阻力與波浪運(yùn)動的耦合。田寶強(qiáng)等[8-9]分別基于牛頓-歐拉方程和拉格朗日方程建立了波浪滑翔器縱剖面動力學(xué)模型。上述研究主要考慮波浪滑翔器在豎直平面上的運(yùn)動，而不能很好地描述波浪滑翔器在水平平面上的運(yùn)動響應(yīng)。Smith等[10]對影響波浪滑翔器前進(jìn)速度的因素的研究。雖然該研究可預(yù)測波浪滑翔器的前進(jìn)速度，但是未考慮波浪滑翔器其他自由度的運(yùn)動。

現(xiàn)有的海洋機(jī)器人動力學(xué)模型通常是基于單體的動力學(xué)分析建立起來的，而波浪滑翔器是一個多體結(jié)構(gòu)，先前的模型適用性不足，不能清楚地表示波浪滑翔器浮體與滑翔體之間的強(qiáng)耦合關(guān)系。D-H方法最早用于描述工業(yè)機(jī)械臂，能夠較好表示多個機(jī)械臂之間的關(guān)系，且理論發(fā)展較成熟，于是Caiti等[11]提出可以借鑒D-H方法用于表示浮體與滑翔體之間的耦合關(guān)系。田寶強(qiáng)等[12]在其基礎(chǔ)上進(jìn)一步完善，但是僅僅考慮到二維情況。

本文在國內(nèi)外波浪滑翔器的動力學(xué)模型研究基礎(chǔ)上，結(jié)合Caiti A的理論基礎(chǔ)，利于D-H方法建立三維的波浪滑翔器動力學(xué)方程，并進(jìn)行運(yùn)動仿真。通過對模型進(jìn)行動力學(xué)仿真，對波浪滑翔器進(jìn)行動力學(xué)分析與操縱性預(yù)報。

1 波浪滑翔器

1.1 基本結(jié)構(gòu)

波浪滑翔器總體上分為水面浮體、水下滑翔體和系索3大部分，如圖1、2所示。

圖1 工作中的波浪滑翔器Fig.1 Wave glider at work

圖2 波浪滑翔器三維模型Fig.2 3-D model of wave glide

水面浮體由太陽能電池板、浮體材料包圍的密封艙體、控制系統(tǒng)、各類傳感器負(fù)載以及可充電電池等組成；滑翔體由可轉(zhuǎn)動的翼板、翼板支撐框架及舵機(jī)組成；系索主要起到連接浮體和滑翔體作用，并且傳遞信息。

1.2 推進(jìn)原理

如圖3所示，當(dāng)波浪抬升水面浮體時，在系索的拉力下，水下滑翔體做上升運(yùn)動。受到水動力的影響，水翼板尾部向下偏轉(zhuǎn)。當(dāng)攻角在一定范圍內(nèi)時，水翼板產(chǎn)生升力，其分力可分解到水平和垂直2個方向。其中，水平分力推動水下滑翔體向前運(yùn)動，繼而拉動水面浮體，促使其前進(jìn)。

圖3 波浪滑翔器上升和下降運(yùn)動Fig.3 Wave glider rises and falls

當(dāng)水面浮體越過波峰時，受到重力的影響，整個系統(tǒng)將向下運(yùn)動，由于水動力作用，這時水翼板尾部向上翻轉(zhuǎn)。與上升過程一樣，由于水動力作用會有升力產(chǎn)生，水平分力方向向前，使整個系統(tǒng)向前運(yùn)動[13]。

2 動力學(xué)模型

2.1 假設(shè)

為了簡化波浪滑翔器動力學(xué)模型建模的復(fù)雜程度，引用文獻(xiàn)[11-12]，對波浪滑翔器以及它的環(huán)境做出如下假設(shè)：

1)浮體和滑翔體被剛性連接，并且質(zhì)量恒定不變。

2)系索應(yīng)總是處于緊張狀態(tài)，而且不會在水動力阻力作用下彎曲。

3)系索節(jié)點與浮體重心、滑翔體重心之間的距離忽略不計，因為它們相對于系索的長度都是十分短的。

2.2 D-H方法

D-H方法是由Denavit等[14]在1955年最先提出，用4×4齊次變換矩陣來描述機(jī)械臂的2個相鄰連桿之間的位置和姿態(tài)。目前，它已成為一種通用方法，廣泛應(yīng)用于機(jī)械臂動力學(xué)領(lǐng)域[15]。在每個連桿上固定一個坐標(biāo)系，來描述其相對位置和姿態(tài)關(guān)系。所以，每一個連桿可以用ai-1、αi-1、di和θi這4個參數(shù)來描述，其中ai-1和αi-1用來描述連桿i-1自身，di和θi用來描述2個相鄰連桿之間的相對位置關(guān)系。從連桿坐標(biāo)系i-1變換到坐標(biāo)系i的坐標(biāo)轉(zhuǎn)換矩陣可以通過以下步驟獲得：繞xi-1軸旋轉(zhuǎn)一個角度αi-1，使得zi軸與zi-1軸互相平行；沿xi-1軸平移一段距離ai-1，使得zi軸與zi-1軸重合；繞zi軸旋轉(zhuǎn)一個角度θi，使得xi-軸與xi-1軸互相平行；沿zi軸平移一段距離di，使得xi-軸與xi-1軸重合。

根據(jù)以上步驟，可得知αi-1為兩坐標(biāo)系z軸之間的夾角，ai-1為兩坐標(biāo)系z軸之間的距離，θi為兩坐標(biāo)系x軸之間的夾角，di為兩坐標(biāo)系x軸之間的距離。

基于D-H方法建立波浪滑翔器的參考坐標(biāo)系，如圖4所示。為了描述其多體系統(tǒng)，假設(shè)基本的參考系xbybzb和水下機(jī)器人通常用于導(dǎo)航的坐標(biāo)系一樣是北-東-下(NED)坐標(biāo)系，而xb1yb1zb1是一個過渡坐標(biāo)系，實現(xiàn)從xbybzb到x1y1z1的坐標(biāo)系統(tǒng)轉(zhuǎn)換。

圖4 波浪滑翔器坐標(biāo)系系統(tǒng)Fig.4 Coordinate systems of wave glider

根據(jù)D-H方法的定義，波浪滑翔器對應(yīng)的D-H參數(shù)在表1中總結(jié)。

其中d1、d2和d3分別是浮體在x、y、z方向上的位移，θ1是浮體的航向角，θ2為浮體的縱傾角，θ3為浮體船首方向與系索的夾角，θ4為系索與滑翔體船首方向的夾角，θ5為浮體與滑翔體航向角的夾角，a1為系索長度。

表1 波浪滑翔器D-H參數(shù)表Table 1 Wave glider D-H parameters

所以從滑翔體坐標(biāo)系8變換到慣性坐標(biāo)系b的坐標(biāo)轉(zhuǎn)換矩陣可以由式(1)計算所得：

(1)

2.3 運(yùn)動傳遞

根據(jù)剛體動力學(xué)，2個相鄰的連桿具有關(guān)系：

對于旋轉(zhuǎn)關(guān)節(jié)，有：

(2)

對于移動關(guān)節(jié)，有：

(3)

于是，每個連桿的質(zhì)心速度可表示為：

ivci=ivi+iωi×ipci

(4)

并且，質(zhì)心速度相對于慣性坐標(biāo)系可表示為：

(5)

式中：ivi和iωi分別是在連桿坐標(biāo)系i下線速度vi與角速度ωi；i+1zi+1是在坐標(biāo)系i+1下z軸方向的單位向量；ipi+1是在坐標(biāo)系i下，坐標(biāo)系i+1的位置向量。

2.4 動能與勢能

考慮波浪滑翔器的運(yùn)動穩(wěn)定性，水下滑翔體可以通過結(jié)構(gòu)設(shè)計和平衡特性來實現(xiàn)平移?？梢缘玫綆缀侮P(guān)系：

θ2+θ3+θ4=0

(6)

即浮體的縱傾角、浮體船首方向與系索的夾角、系索與滑翔體船首方向的夾角，三者之和為0。

為了方便地計算能量，關(guān)節(jié)變量可以被統(tǒng)一替代：

(7)

如果在慣性坐標(biāo)系下，連桿質(zhì)量為mi，平移速度為vci，質(zhì)心角速度為ωi，并且相對應(yīng)質(zhì)心的慣性張量為Ii。這個連桿的動能可表示為：

(8)

波浪滑翔器的動能可表示為：

(9)

考慮到波浪滑翔器在水中航行，這將導(dǎo)致水加速并產(chǎn)生附加質(zhì)量效應(yīng)。所以，周圍流體的動能為Tf。所以，整體的動能應(yīng)為波浪滑翔器的動能與周圍流體的動能之和:

(10)

式中M是廣義慣性矩陣，M=Mw+Mf。

如果將慣性坐標(biāo)系原點的水平設(shè)為零勢能面，系統(tǒng)的總勢能可以寫成：

(11)

由于固定坐標(biāo)系x1y1z1、x2y2z2、x3y3z3、x4y4z4的連桿1、2、3、4都是虛擬的，因此它們都沒有質(zhì)量，即m1=m2=m3=m4=0。而m5、m6、m7分別等于浮體、系索、滑翔體的質(zhì)量。

2.5 動力學(xué)方程

首先定義拉格朗日函數(shù)為：

L=T-U

(12)

于是，拉格朗日方程可表示為：

(13)

所以，在計算和分析的基礎(chǔ)上，動力學(xué)模型可以以向量形式寫成：

(14)

式中M是廣義慣性矩陣。若用mf1、mf2、mf3、If1、If2、If3、If4表示附加慣性質(zhì)量(轉(zhuǎn)動慣量)，則它的每一元素的值為：

M11=m5+m6+m7+mf1
M22=m5+m6+m7+mf2
M33=m5+m6+m7+mf3
M55=I5z+If2，M77=I7z+If4
M44=(0.25m6+m7)(a1cosθ4)2+I5xsin2θ2+
I5ycos2θ2+I6xsin2θ4+I6ycos2θ4+
I7z+If1
M66=(0.25m6+m7)a12+I6z+If3
M14=M41=-(0.5m6+m7)a1sinθ1cosθ4
M16=M61=-(0.5m6+m7)a1cosθ1sinθ4
M24=M42=(0.5m6+m7)a1cosθ1cosθ4
M26=M62=-(0.5m6+m7)a1sinθ1sinθ4
M36=M63=(0.5m6+m7)a1cosθ4
M47=M74=I7z

矩陣M的其他元素全部為零。

G(q)是重力向量，它的每一元素的值為：

G1=G2=G4=G5=G7=0

G3=-(m5+m6+m7)g

G6=-(0.5m6+m7)a1gcosθ4

根據(jù)在虛功原理下的虛位移的廣義坐標(biāo)表達(dá)式，可獲得廣義力τ的表達(dá)式，為：

式中：Fpx和Fpz分別是滑翔體水翼產(chǎn)生的合力在水平和垂直2個方向上(滑翔體的隨體坐標(biāo)系)的分力；Df、Dl和Dg分別是浮體、系索和滑翔體的阻力；Bf、Bl和Bg分別是浮體、系索和滑翔體的浮力；Fsx和Fsy分別是滑翔體轉(zhuǎn)舵裝置在x軸方向和y軸方向上(滑翔體的隨體坐標(biāo)系)的分力，而τs是轉(zhuǎn)舵裝置對滑翔體產(chǎn)生的力矩；Fw和τw分別是波浪對浮體的力和恢復(fù)力矩。τg和τf分別是水對浮體和滑翔體的阻力力矩。

3 數(shù)值仿真與水池試驗

3.1 運(yùn)動仿真系統(tǒng)搭建

將其式(14)中所有加速度項移到方程的左端，所有的非加速度項移到方程的右端，可以得到動力學(xué)方程：

(15)

若已知波浪滑翔器當(dāng)前時刻的運(yùn)動狀態(tài)參數(shù)，就可以得到式(15)等號右邊的合值，對常數(shù)矩陣M求逆，就可以求解該方程，獲得波浪滑翔器的加速度，進(jìn)而確定波浪滑翔器下一時刻的運(yùn)動參數(shù)，如此往復(fù)迭代就可以使得波浪滑翔器在仿真系統(tǒng)中運(yùn)動起來[16]。

3.2 仿真模型參數(shù)

本文所采用的模型基本參數(shù)主要參考文獻(xiàn)[17]。主要參數(shù)如表2。

表2 波浪滑翔器模型基本參數(shù)Table 2 Parameters of wave glider model

由于系索又細(xì)又輕，在計算時可以忽略它的質(zhì)量和轉(zhuǎn)動慣量。在計算轉(zhuǎn)動慣量時，由于實際質(zhì)量分布復(fù)雜，參考一般船型，可采用經(jīng)驗公式來計算：

(16)

波浪滑翔器整個系統(tǒng)的在各個方向上的附加質(zhì)量，在x、y和z方向上分別是總質(zhì)量的0.1、2.1和2.1倍[18-19]。在計算附加慣性矩時，進(jìn)行如下估算：Jxx≈0.1Ixx，Jyy≈Iyy，Jzz≈Iyy，If1≈I7z，If2≈I5z，If3≈I6z=0，If4≈I7z。仿真計算時，采用文獻(xiàn)[17]中計算的不同波幅、不同周期的正弦波浪條件下，水翼產(chǎn)生的驅(qū)動力，見表3。

表3 滑翔體在不同風(fēng)浪等級下產(chǎn)生的驅(qū)動力Table 3 The driving force of the glider in different wave levels

而垂直方向上阻力的作用力被認(rèn)為是相對于波浪力是較小的，在本文中忽略不計。波浪滑翔器航行速度較慢，阻力的主要成分是摩擦阻力。根據(jù)弗勞德假設(shè)，摩擦阻力等于相同速度、相同長度、相同濕表面積的相當(dāng)平板摩擦阻力[19]：

(17)

利用表4中浮體，滑翔體以及系索的阻力系數(shù)，可粗略地求得Df、Dg和Dl。

表4 阻力系數(shù)和濕表面積估算Table 4 Resistance coefficient and wet surface area

本文采用舵的剖面形狀為NACA0012型，取弦長0.2 m，展長0.12 m，舵面積0.024 m2，舵軸安裝于距離前緣0.3倍的弦長處。轉(zhuǎn)舵裝置升力系數(shù)CL，阻力系數(shù)CD分別表示為：

(18)

為了便于計算，將升力系數(shù)CL與舵角δ近似視為正比關(guān)系，取舵角δ為30°時的升力系數(shù)CL=1.0，將升阻比視為一定值，本文將升阻比取為2，可得：

(19)

3.3 水池試驗

為驗證此動力學(xué)模型的正確性和有效性，本文利用哈爾濱工程大學(xué)的“海洋漫步者”號波浪滑翔器的水池試驗數(shù)據(jù)進(jìn)行比對。該波浪滑翔器裝備了多種傳感器，可以測量浮體和滑翔體的速度，艏向角，首搖角速度等信息。“海洋漫步者”號浮體質(zhì)量為55 kg，滑翔體質(zhì)量為40 kg，系索長4 m。

在波高230 mm、波長8 m條件下進(jìn)行直航試驗，仿真試驗和水池試驗的對比結(jié)果如圖5所示。

圖5 浮體縱向速度對比Fig.5 Comparison of longitudinal velocity of the float

圖5顯示了仿真和水池試驗中浮體縱向速度的對比，從圖中可以看出水池試驗值和仿真計算值有較高的吻合度，誤差較小。誤差原因可能是水池試驗中傳感器的采樣頻率和測量精度有限。

在波高230 mm、波長8 m條件下進(jìn)行Z形試驗，仿真舵角輸入和水池舵角相同，結(jié)果如圖6所示。

從圖6(a)和(b)中可以看出，無論是浮體還是滑翔體，當(dāng)舵角變化時，在仿真和水池試驗中的航向響應(yīng)的變化趨勢是相似的。從圖6(c)和(d)中可以看出，浮體與滑翔體的艏向角的變化趨勢在仿真和水池試驗中基本一致。仿真是在理想環(huán)境下進(jìn)行的，而水池試驗存在機(jī)械結(jié)構(gòu)阻力、不規(guī)則波載荷、傳感器噪聲等多種不確定因素，水動力系數(shù)和建模假設(shè)有限的精度給仿真帶來誤差，導(dǎo)致仿真結(jié)果與水池試驗結(jié)果存在差異。

根據(jù)上述仿真與水池試驗的對比，可以判斷基于D-H方法的波浪滑翔器動力學(xué)模型是具備一定正確性和有效性的，用于仿真分析是可行的。

圖6 仿真與水池試驗結(jié)果對比Fig.6 Comparison of simulation and tank test results

3.4 數(shù)值仿真試驗

本文數(shù)值仿真試驗由直航試驗、回轉(zhuǎn)試驗、Z形試驗3個部分組成。

直航試驗是分別在1～3級海況下，舵角為0°，波浪滑翔器按縱向直線航行。通過仿真獲得波浪滑翔器分別在的縱向速度V1、浮體縱傾角θ2、系索夾角θ4以及浮體升沉d3的曲線如圖7所示。

圖7 直航試驗結(jié)果對比Fig.7 Comparison of simulation in the tracking of straight line

從圖7(a)中可以看出，3種海況下，波浪滑翔器穩(wěn)定速度分別為0.11、0.35和0.6 m/s。在一定條件下，波浪滑翔器的速度隨海況增強(qiáng)而增大。從圖7(b)中可以看出，3種海況下，浮體的縱傾角最大振幅分別為1.5°、15°、20°，浮體的縱傾角最大振幅隨海況增強(qiáng)而增大。從圖7(c)中可以看出，系索夾角存在微幅震動，最終基本上在90°附近微幅震動，海況大小對其大小影響不大。從圖7(d)中可以看出，浮體升沉受波浪影響，海況越大波浪振幅和頻率越大，浮體升沉的也振幅和頻率越大。從整個直航試驗可知，可知波浪波幅和頻率影響波浪滑翔器對波浪能的利用率，海況越大，波浪能越高，波浪滑翔器航行越快。

回轉(zhuǎn)試驗時，假定在3級海況下，給定垂直舵角，分別為5°、10°、15°、20°、25°、30°。由波浪滑翔器的回轉(zhuǎn)試驗仿真曲線可得回轉(zhuǎn)直徑，回轉(zhuǎn)直徑與舵角的關(guān)系如圖8所示。

圖8 回轉(zhuǎn)直徑與舵角的關(guān)系曲線Fig.8 Relation between turning diameter and rudder angle

由圖8可知，在3級海況下，波浪滑翔器的回轉(zhuǎn)直徑會隨舵角增大而減小，最終收斂約為12 m，即最小回轉(zhuǎn)直徑，它大約是浮體長度的4倍。

在Z形試驗時，同樣假定在三級海況下，仿真初始舵角為10°，當(dāng)浮體艏向角改變到達(dá)10°時，舵角階躍為-10°，當(dāng)浮體艏向角改變到達(dá)-10°時，舵角階躍為10°，然后按此步驟循環(huán)。仿真結(jié)果如圖9所示。

圖9 舵角，浮體與滑翔體艏向角變化曲線Fig.9 Relation among δ,θ1,(θ1+θ5) and t

從圖9中可以看出，當(dāng)0 s開始操舵時，波浪滑翔器順時針運(yùn)動。在7 s時，舵角從10°階躍到-10°，波浪滑翔器做逆時針運(yùn)動。分析可得其初轉(zhuǎn)期為7 s，超約時間為3 s，超越角為3°，周期為20 s。

通過仿真結(jié)果之間的對比，可以發(fā)現(xiàn)，采用本文方法仿真，可以得到浮體縱向速度v1、浮體橫向速度v2、浮體升沉d3、浮體艏向角θ1、浮體縱傾角θ2、滑翔體艏向角(θ1+θ5)、系索與滑翔體夾角θ4等參數(shù)。其中，浮體升沉d3和浮體縱傾角θ2是一些傳統(tǒng)方法不能計算的。而這2個參數(shù)對波浪滑翔器波浪能利用率的計算起著關(guān)鍵作用。因此，本文所建立的波浪滑翔器動力學(xué)模型，具有較好的仿真結(jié)果。

4 結(jié)論

1)該波浪滑翔器動力學(xué)模型在利用D-H方法的基礎(chǔ)上，更好地表達(dá)浮體與滑翔器之間的強(qiáng)耦合關(guān)系。

2)該模型以三維的形式，提供更多運(yùn)動參數(shù)，方便計算波浪滑翔器波浪能利用效率，且部分參數(shù)與水池試驗結(jié)果相吻合。

3)仿真結(jié)果驗證了不同海況下的波浪滑翔器運(yùn)動特性，包括其速度、角速度和系索夾角等，其結(jié)果與動力學(xué)理論相一致。

然而，水動力計算、模型參數(shù)估算、環(huán)境影響等問題還需要進(jìn)一步研究。