王沖霄，劉忠樂，文無敵，張志強(qiáng)，趙 苗

(海軍工程大學(xué) 兵器工程學(xué)院， 武漢 430033)

水下拖曳系統(tǒng)，是一種廣泛應(yīng)用于海洋監(jiān)測、海洋研究以及軍事等領(lǐng)域的水下探測裝置，在開發(fā)海洋的先進(jìn)技術(shù)手段中，拖曳系統(tǒng)裝備具有極其重大的意義[1]。線列陣是水下拖曳系統(tǒng)的一個(gè)關(guān)鍵組成部分，其水動(dòng)力特性直接影響整個(gè)水下航行器系統(tǒng)的快速性、操縱性及穩(wěn)定性，開展水下拖纜的水動(dòng)力特性研究具有重要的理論意義和工程實(shí)用價(jià)值。近年來，基于響應(yīng)面的優(yōu)化方法已廣泛應(yīng)用在穩(wěn)健設(shè)計(jì)和多目標(biāo)與多學(xué)科優(yōu)化設(shè)計(jì)的代理模型里[2]。這種近似模型技術(shù)是在初始數(shù)據(jù)集合基礎(chǔ)上構(gòu)造逼近目標(biāo)函數(shù)和約束條件的方法，同時(shí)也為快速優(yōu)化和敏感性分析提供了一種高效的解決方法[3]。本文針對水下拖纜，在王飛[4-5]的穩(wěn)態(tài)運(yùn)動(dòng)求解和分析的基礎(chǔ)上，引入多目標(biāo)優(yōu)化設(shè)計(jì)理論，將復(fù)雜求解方法所得結(jié)果進(jìn)行回歸處理，建立起拖纜尾端拖曳深度和首端張力的二次響應(yīng)面模型，給出多目標(biāo)優(yōu)化算法下的Pareto最優(yōu)解集，并分析了拖纜參數(shù)對尾端拖曳深度和首端張力的影響。

1 拖纜運(yùn)動(dòng)模型及數(shù)值計(jì)算

為簡化討論過程，本文僅討論穩(wěn)定海流下自由端拖纜的運(yùn)動(dòng)模型，其拖曳系統(tǒng)如圖1。

圖1 拖曳系統(tǒng)示意圖

其他的拖曳形式，例如尾端加載拖船、尾繩、水下拖體及其他的纜載設(shè)備，相比自由端僅在邊界條件上存在區(qū)別，此處不詳細(xì)討論。將拖纜視為理想的柔性圓形纜繩，由水下航行器(AUV)搭載，建立空間固定的慣性坐標(biāo)系O-XYZ，單位矢量定義為(i,j,k)，附著在拖纜上的局部坐標(biāo)系btn，單位矢量定義為(b,t,n)。軸t表示拖纜切向，方向?yàn)槔|長s的增長方向；軸n表示拖纜的法向，處在軸t和軸t在OXY平面內(nèi)的投影所組成的平面內(nèi)，并垂直于軸t；軸b與軸n和軸t共同組成右手笛卡爾坐標(biāo)系。定義歐拉角θ，φ為拖纜微元相對慣性坐標(biāo)系的姿態(tài)角，θ為Otn平面偏離OX軸的角度，φ為軸t偏離OXY平面的角度，θ∈(-180°,180°],φ∈(-90°，90°]，兩個(gè)歐拉角均以逆時(shí)針方向?yàn)檎较?。慣性坐標(biāo)系和局部坐標(biāo)系通過姿態(tài)角相關(guān)聯(lián)，轉(zhuǎn)換關(guān)系如下：

(1)

(2)

對于每個(gè)拖纜微元ds，在穩(wěn)定直航的狀態(tài)下，其重力、浮力、流體阻力的合力處于平衡狀態(tài)，有以下平衡方程：

(3)

式中:T代表拖纜張力，始終指向拖纜的切向，B和G分別代表拖纜單位長度的浮力和重力，D代表流體阻力。

單位矢量(b,t,n)對拖纜長度s微分用“′”表示，則有[6]：

(4)

則可得：

(5)

所以張力在局部坐標(biāo)系可展開為下式：

(6)

將式(6)與重力、浮力和流體阻力代入平衡方程式(3)中，并在拖纜局部坐標(biāo)系下沿各坐標(biāo)軸方向展開，則平衡方程寫為如下標(biāo)量形式[7]：

(7)

式中:w為拖纜單位長度在水中的質(zhì)量，表示為：w=(u-ρσ)g，其中u為拖纜單位長度質(zhì)量，ρ為流體密度，σ為橫截面積，d為拖纜直徑，ε為拖纜應(yīng)變，Ct和Cn分別為拖纜的切向和法向阻力系數(shù)，ut，ub，un為局部坐標(biāo)系下的速度分量，由系統(tǒng)相對于水流的拖曳速度轉(zhuǎn)換到局部坐標(biāo)系下得到，有：

(8)

式中:v為拖曳速度，J為海流，僅有水平面內(nèi)速度分量而無垂向分量。拖纜在慣性系下的坐標(biāo)由下式給出：

(9)

在均勻海流狀態(tài)下，拖纜自由端即s=0處的張力和歐拉角變化率均為零，拖纜側(cè)向無作用力存在，有θ≡const，此時(shí)可將穩(wěn)態(tài)問題轉(zhuǎn)換到二維空間下求解，在自由端的邊界條件表示如下[4]：

(10)

式中ψ為AUV的航行艏向角，將以上初始值方程和控制方程(7)、(9)聯(lián)立，采用四階龍格庫塔方法進(jìn)行積分計(jì)算，即可求出穩(wěn)態(tài)解。

2 近似模型的建立

正常情況下，拖曳系統(tǒng)均工作于穩(wěn)定狀態(tài)，而穩(wěn)態(tài)下的拖曳深度和拖纜首端張力是設(shè)計(jì)過程中兩個(gè)非常重要的指標(biāo)，拖纜的密度、楊氏模量、阻力系數(shù)以及拖曳速度和流體密度，均會(huì)對此產(chǎn)生一定的影響。王飛[5]對穩(wěn)態(tài)拖曳各參數(shù)在一定的取值范圍下，進(jìn)行了拖曳深度和拖纜首端張力的初步研究，但每個(gè)參數(shù)的變化，是在其余參數(shù)不變的基礎(chǔ)上進(jìn)行的，而且最后僅對各參數(shù)對結(jié)果的影響進(jìn)行了定性分析。以下將采用試驗(yàn)設(shè)計(jì)的方法，在王飛給出的拖纜模型基礎(chǔ)上，構(gòu)建參數(shù)和響應(yīng)之間的近似模型，并得出尾端拖曳深度和首端張力的近似計(jì)算公式。

2.1 試驗(yàn)設(shè)計(jì)

本文選擇拉丁超立方設(shè)計(jì)(Latin hypercube designs)作為試驗(yàn)方法。拉丁超立方設(shè)計(jì)方法是一種優(yōu)秀的強(qiáng)調(diào)樣本點(diǎn)分布均勻的試驗(yàn)設(shè)計(jì)方法[8]。該方法假設(shè)需要n個(gè)實(shí)驗(yàn)設(shè)計(jì)點(diǎn)，則設(shè)計(jì)變量會(huì)被分成n等分，在每等分中選擇一個(gè)參數(shù)作為設(shè)計(jì)點(diǎn)。根據(jù)問題規(guī)模和復(fù)雜程度，樣本點(diǎn)的個(gè)數(shù)也應(yīng)該適當(dāng)增多。通常情況下，對于5～10個(gè)變量的問題，樣本點(diǎn)數(shù)量建議取為1.5×(n+1) ×(n+2)/2個(gè)[9]，因此針對本文的6個(gè)參數(shù)變量，選擇42個(gè)樣本點(diǎn)進(jìn)行試驗(yàn)設(shè)計(jì)。拖纜模型和參數(shù)取值范圍采用文獻(xiàn)[5]給出的模型，參數(shù)如表1所示：

表1 拖纜參數(shù)

尾繩置于引導(dǎo)纜尾端，起到定深和穩(wěn)定的作用。針對引導(dǎo)纜的物理參數(shù)以及流體密度ρ和拖曳速度vy設(shè)定取值范圍如下：

108N/m2≤E≤1011N/m2

0.010≤Ct≤0.03

1.2≤Cn≤1.9

0.8 kg/m≤u≤1.2 kg/m

1 m/s≤vy≤3 m/s

1 020 kg/m3≤ρ≤1 030 kg/m3

2.2 二階多項(xiàng)式響應(yīng)面

二階多項(xiàng)式響應(yīng)面的數(shù)學(xué)表達(dá)式為：

(11)

采用復(fù)相關(guān)系數(shù)R2作為響應(yīng)面模型的誤差分析指標(biāo)，定義如下：

(12)

采用多項(xiàng)式回歸技術(shù)對試驗(yàn)設(shè)計(jì)的樣本點(diǎn)和響應(yīng)值進(jìn)行最小二乘擬合，并求出待定系數(shù)，構(gòu)造出近似模型。回歸模型表示如下：

yi=f(xi,θ)+εi，i=1,2,…，n

(13)

Z=33 539.7-258.38x1-2 728.3x2-1.26×10-9x3-

87.61x1x2-1.21×10-11x1x3+

0.19x1x4-27.43x1x5+12.28x1x6-

3.04×10-11x2x3-3.07x2x4+

121.57x2x5-0.47x2x6+1.22×10-12x3x4+

5.44×10-12x3x5-3.12×10-11x3x6-

0.19×10-3x4x5-0.21x4x6-52.96x5x6

(14)

F=-429 504.9-1 752.43x1-84 097.73x2-

9.14×10-9x3-842.3x4+3 016.54x5-

2.03x1x4-241.36x1x5+40.65x1x6+

3.1×10-9x2x3-52.68x2x4-5 144.58x2x5+

47 585.52x2x6+9.16×10-12x3x4+

5.12×10-11x3x5+1.54×10-11x3x6-

1.26x4x5+1.26x4x6-597.65x5x6

(15)

式中變量(x1,x2,x3,x4,x5,x6)分別代表參數(shù)中的Cn，Ct，E，ρ，u，vy。通過計(jì)算復(fù)相關(guān)系數(shù)來驗(yàn)證響應(yīng)面模型的擬合精度。其中，Z的復(fù)相關(guān)系數(shù)為0.996 96，F(xiàn)的復(fù)相關(guān)系數(shù)為0.999 53。該響應(yīng)面模型擬合程度較好。

3 多目標(biāo)遺傳算法

本文所探討的是如何實(shí)現(xiàn)尾端拖曳深度最大、首端張力最小的多目標(biāo)優(yōu)化問題，則本優(yōu)化問題的數(shù)學(xué)模型可以描述為：

minF, maxZ

findCn,Ct,E,ρ,u,vy

s.t. 108N/m2≤E≤1011N/m2

0.010≤Ct≤0.03

1.2≤Cn≤1.9

0.8 kg/m≤u≤1.2 kg/m

1 m/s≤vy≤3 m/s

1 020 kg/m3≤ρ≤1 030 kg/m3

(16)

大多數(shù)情況下，上式的各個(gè)子目標(biāo)往往是互相沖突的，其中一個(gè)子目標(biāo)的優(yōu)化會(huì)帶來其他子目標(biāo)的損失，所以多目標(biāo)優(yōu)化問題的優(yōu)化解是一個(gè)解集，稱為Pareto最優(yōu)解集，解集中的元素稱為Pareto最優(yōu)解，Pareto最優(yōu)解集在目標(biāo)函數(shù)空間中的像稱為Pareto前沿。在目前流行的優(yōu)化算法中，遺傳算法最適合求解多目標(biāo)優(yōu)化問題的Pareto解集。本文采用遺傳算法中的NSGA-Ⅱ算法，基于ISIGHT平臺(tái)，完成整個(gè)優(yōu)化過程。NSGA-Ⅱ算法利用基于Pareto支配的排序方法將個(gè)體進(jìn)行分層排序，并通過計(jì)算擁擠距離的方式對同一層級個(gè)體進(jìn)行具體排序，具有較高的運(yùn)算效率和較好的收斂速度[10]。本優(yōu)化問題的主要流程如下：

1) 初始化拖纜的二階多項(xiàng)式響應(yīng)面計(jì)算模型，調(diào)入所需參數(shù)，初始化進(jìn)化過程的參數(shù)，并隨機(jī)生成初始種群N；

2) 將初始種群N中的個(gè)體依次賦值給拖纜計(jì)算模型，修改待優(yōu)化參數(shù)并運(yùn)行該模型，輸出對應(yīng)的拖曳深度和張力值，判斷是否滿足約束條件，以懲罰系數(shù)來處理不滿足條件的個(gè)體，使其在進(jìn)化過程中被淘汰；

3) 對種群中所有個(gè)體均以此方法進(jìn)行操作，得到對應(yīng)的輸出值，從而完成種群的初始化工作；

4) 對初始種群N進(jìn)行非支配排序并計(jì)算擁擠距離，得出每個(gè)個(gè)體的優(yōu)劣性指標(biāo)；

5) 以t記錄進(jìn)化次數(shù)，初始值為1，對父種群Nt執(zhí)行交叉、變異等進(jìn)化操作，產(chǎn)生子代種群Na；

6) 對每個(gè)子代種群Na進(jìn)行計(jì)算，得到對應(yīng)目標(biāo)值；

7) 對父代種群Nt和子代種群Na進(jìn)行非支配排序并計(jì)算擁擠距離；

8) 對種群進(jìn)行更新，對上一步驟合并的種群進(jìn)行修正，得到新的種群Nt+1，判斷是否達(dá)到最大代數(shù)，如果是則輸出非支配解集，否則繼續(xù)迭代。

4 計(jì)算結(jié)果

本優(yōu)化過程的種群規(guī)模為24，遺傳代數(shù)為50，交叉概率為0.9，所得深度和張力的Pareto前沿散點(diǎn)分布如圖2，選取3個(gè)優(yōu)化方案，方案1中Z最大，方案2中F最小，方案3位于中間，所得優(yōu)化結(jié)果如表2所示。Pareto前沿為非支配解集，最終方案的選擇依賴于實(shí)際情況和選擇者偏好[11]。

圖2 拖纜優(yōu)化的Pareto前沿散點(diǎn)分布

表2 Pareto解集中的三個(gè)優(yōu)化結(jié)果

本文同時(shí)在ISIGHT平臺(tái)上通過數(shù)據(jù)處理給出各參數(shù)對尾端拖曳深度和首端張力影響的Pareto柱狀圖，如圖3和圖4所示。圖3的第一條、第三條和第六條以及圖4的第四條和第六條表示該參數(shù)與響應(yīng)呈負(fù)相關(guān)關(guān)系，其余表示該參數(shù)與響應(yīng)呈正相關(guān)關(guān)系?？梢钥闯?，拖曳速度vy對尾端拖曳深度和張力的貢獻(xiàn)程度百分比最大，分別達(dá)到了-52.53%和62.56%；其余各參數(shù)中，拖纜密度u和法向阻力系數(shù)Cn對尾端拖曳深度貢獻(xiàn)程度百分比較大，分別為32.47%和-12.24%，切向阻力系數(shù)Ct和拖纜密度u對張力貢獻(xiàn)程度百分比較大，分別為22.18%和12.12%，而楊氏模量E和流體密度ρ對兩者的貢獻(xiàn)程度均可忽略不計(jì)。以上分析結(jié)果基本和文獻(xiàn)[5]所得結(jié)論一致。

圖3 參數(shù)對拖曳深度影響

圖4 參數(shù)對首端張力影響

5 結(jié)論

基于ISIGHT平臺(tái)實(shí)現(xiàn)了參數(shù)的多目標(biāo)優(yōu)化并得到Pareto解集，為設(shè)計(jì)者提供了重要的設(shè)計(jì)信息；通過計(jì)算各參數(shù)對尾端拖曳深度和首端張力的影響，驗(yàn)證了前人的結(jié)論。若拖纜模型和參數(shù)取值范圍有變，本文所得出的尾端拖曳深度和張力的近似計(jì)算公式相應(yīng)會(huì)發(fā)生變化。