基于主體參數(shù)化分析的潛水器多學科優(yōu)化

劉 峰 趙彥凱 姚競爭 王 賀

哈爾濱工程大學船舶工程學院，哈爾濱，150001

0 引言

潛水器在海洋開發(fā)和利用中占有重要的地位，隨著人類海洋活動的增加和技術的進步，對潛水器綜合性能、設計效率和設計周期等提出了更高的要求。傳統(tǒng)潛水器設計采用串行模式，對學科間的耦合效應處理不夠充分[1]，矛盾協(xié)調(diào)和計算量大,設計效率低，且學科間的相互影響得不到充分利用，導致最佳設計方案獲取困難。多學科設計優(yōu)化在實現(xiàn)了復雜系統(tǒng)的學科集成、有效解耦、設計過程的有效組織和管理等基礎上，通過充分利用學科間的協(xié)同效應，最終獲得系統(tǒng)最優(yōu)解[2]，可在一定程度上克服傳統(tǒng)設計模式的不足，已在包括潛水器在內(nèi)的眾多領域得到了應用，并取得了明顯的優(yōu)化效果。

潛水器多學科優(yōu)化過程涉及多個學科，建立能夠?qū)W科準確表達的計算模型，對優(yōu)化方案準確度和設計效率提升等具有重要意義。目前，潛水器多學科優(yōu)化計算模型多采用經(jīng)驗公式[3-7]，盡管易于實現(xiàn)，但精度不高。在眾多學科中，艇型和結構構成了潛水器主體，對其總體性能影響最大,因此，建立準確的艇型、結構計算模型就成為潛水器多學科優(yōu)化研究的重點。應用CAE軟件進行艇型和結構分析具有精度高、計算成本低等優(yōu)點，但需要反復修改和分析，設計效率低[8]。參數(shù)化可實現(xiàn)模型的重建和自動分析，已在潛水器阻力[8-9]、結構[10-11]分析與優(yōu)化中得到了應用，但潛水器多學科設計優(yōu)化過程需要對多個方案進行分析，計算成本依然很高。近似模型采用數(shù)學模型對變量和響應進行逼近，可在保證模型精度的前提下，提高設計效率，因此，針對潛水器主體中的艇型、結構選擇樣本點進行參數(shù)化分析，利用近似模型建立計算模型可平衡計算精度和設計效率之間的矛盾。

本文進行了一型潛水器的學科分解，設計了主體參數(shù)化分析流程。選擇樣本點進行了艇型和結構的參數(shù)化分析，建立了擬合精度滿足要求的艇型、結構近似模型，并完成了其他學科計算模型的建立，對設計變量與目標函數(shù)之間的關系進行分析，建立了基于多學科可行優(yōu)化方法的多學科設計優(yōu)化模型，進行了優(yōu)化求解，為潛水器多學科設計優(yōu)化提供了參考。

1 潛水器學科分解與設計參數(shù)

1.1 學科分解

目標潛水器最大潛深700 m，巡航速度2.5節(jié)，最大航速4節(jié)，主要由結構、推進、導航控制、能源、通信、壓載、觀通、姿態(tài)調(diào)節(jié)、載荷等組成。潛水器總體設計模型見圖1。

圖1 潛水器總體設計模型Fig.1 Model of overall design of submersible

圖1中，導航、運動控制、能源、通信、壓載、觀通、姿態(tài)調(diào)節(jié)、均衡系統(tǒng)等的指標，如用電量、質(zhì)量、浮容積等為常量，將這些系統(tǒng)納入機電設備學科；直航阻力是艇型學科主要研究的內(nèi)容，艇型出于總布置、水動力等方面的考慮是可變的，而艇型變化會對結構、推進等產(chǎn)生影響，推進、機電設備直接影響能源。最終將潛水器劃分為艇型、結構、推進、能源、機電設備共5個學科。

1.2 設計參數(shù)的確定

潛水器的質(zhì)量mz對其造價、母船配套、維護使用等均有重要影響；儲備浮力Fm將直接影響潛水器所能攜帶的載荷、任務模塊，以及升級改造等；直航阻力Fx在潛水器能源消耗中占比最高，對其快速性、續(xù)航能力等均有重要影響。因此，確定Fx、mz、Fm為目標函數(shù)，見表1。狀態(tài)變量分為系統(tǒng)級和學科級，具體見表2。潛水器艇型、結構構成的主體見圖2。結合圖2，確定設計變量見表3。潛水器的強度和穩(wěn)定性、儲備浮力、舷間距離、總質(zhì)量等需滿足表4的約束條件。總體設計模型中的數(shù)據(jù)傳遞關系見圖3。

表1 目標函數(shù)列表Tab.1 List of objective functions

表2 狀態(tài)變量列表Tab.2 List of state variables

圖2 潛水器主體及參數(shù)Fig.2 Main body and parameters of submersible

表3 設計變量列表

表4 約束列表Tab.4 List of constraints

圖3 總體設計模型中的數(shù)據(jù)傳遞關系Fig.3 Data exchange relationship of general design model

2 主體參數(shù)化分析

2.1 試驗設計

拉丁超立方設計(Latin hypercube design，LHD)具有空間填充能力有效、可擬合非線性響應的優(yōu)點，但存在不可重復性、分布不均等不足。最優(yōu)拉丁超立方方法(Opt LHD)是LHD的改進，可有效改善設計空間的均勻性，使構建的近似模型更加準確，擬合效果更好。本文采用Opt LHD方法進行分析。

2.2 艇型參數(shù)化分析

艇型學科主要研究直航阻力，利用RANS法求解。RANS法在時均化N-S方程中，對瞬態(tài)脈動量以某種模型方式體現(xiàn)。雷諾時均化后的不可壓連續(xù)性方程為

(1)

RANS方程為

(2)

(3)

(4)

式中，Gk、Gω為湍流動能，由平均速度梯度引起；Yk、Yω分別為與k和ω相關湍流耗散項；Γk、Γω分別為k和ω的有效擴散項；Dω為正交發(fā)散項。

利用STAR-CCM+軟件進行直航阻力計算，以模型特征長度向四周拉伸相同距離，建立計算域(尾部縱向拉伸特征長度的2倍)。采用對稱平面邊界條件將模型以中縱剖面分成兩部分，網(wǎng)格如圖4所示。

圖4 阻力計算網(wǎng)格Fig.4 Resistance calculation grid

圖4中，網(wǎng)格在艏、艉設置外、中、內(nèi)三層加密網(wǎng)格，采用切割體網(wǎng)格生成器與棱柱層網(wǎng)格生成器，啟用表面重構和自動表面修復功能，網(wǎng)格基礎尺寸0.3 m，最小表面尺寸6%，面網(wǎng)格增長率1.3，棱柱層6層、延伸1.3、總厚度為15%。

利用STAR-CCM+軟件Java宏錄制功能，對Java宏文件進行二次編譯，通過*.bat批處理文件直接運行，調(diào)用cmd.exe并運行批處理文件中的各個命令，艇型參數(shù)化分析流程見圖5。

圖5 艇型參數(shù)化分析流程Fig.5 Parametric analysis process of boat shape

圖5中，通過編寫 * .bat文件調(diào)用STAR-

CCM+，初始設置以*.sim格式保留在文件中，Java宏文件進行外形參數(shù)設置的保留。分析過程涉及的控制參數(shù)、計算結果等利用iSight軟件進行后處理，將計算生成的結果文件存放于相應的工作目錄下。

2.3 結構參數(shù)化分析

根據(jù)規(guī)范[13]計算載荷pjs，按每個標準大氣壓(10 m水深)條件下的壓力0.0098 MPa換算成

pjs=0.0098hjs

(5)

hjs=Khjx

(6)

hjx=hgz/(0.85～0.90)

(7)

式中，hjs為計算深度；K為安全系數(shù)，取1.5；hjx為極限下潛深度；hgz為工作下潛深度。

采用ABAQUS軟件進行分析，耐壓結構受到的載荷左右對稱，在進行邊界條件設置時，簡化為固定于結構的剛體位移。兩封頭的端點X軸和Y軸、中縱剖面處的所有節(jié)點在Y軸和Z軸的位移均約束為零，網(wǎng)格采用0.1×0.1進行劃分，見圖6。

圖6 耐壓結構網(wǎng)格劃分Fig.6 Grid generation of pressure structure

穩(wěn)定性分析分為線性和非線性屈曲兩部分，后者采用弧長法[14]：

(8)

k=1,2,…,N

(1)將運動變量變化引起的應變βN和柯西應力σ進行雙點積，積分得到內(nèi)部節(jié)點應力矩陣IN，對IN和位移uM進行偏微分，M為列數(shù)，得到剛度矩陣KNM，即

(9)

(10)

式(10)中的值很小則收斂，反之則需要進一步按照下式求解：

(11)

(12)

(13)

(4)繼續(xù)進行迭代:

(14)

為實現(xiàn)結構參數(shù)化分析，需要針對ABAQUS進行二次開發(fā)，結構建模、屬性設置、載荷和邊界條件設置、網(wǎng)格劃分、分析和后處理等均采用Python語言編寫，后臺執(zhí)行利用ABAQUS批處理文件*.bat實現(xiàn)，生成結果文件后，進行后處理。為實現(xiàn)自動連續(xù)計算，提高分析效率，在ABAQUS二次開發(fā)的基礎上，利用iSight進行ABAQUS的集成，耐壓結構參數(shù)化分析流程見圖7。

圖7 耐壓結構參數(shù)化分析流程Fig.7 Parametric analysis process of pressure structure

圖7的實現(xiàn)過程為：①利用Simcode對輸入文件*.py和*.bat進行解析，在后臺模式下，利用ABAQUS進行強度分析，然后進行模態(tài)分析，得到相應的特征值；②利用Data Exchanger獲取的特征值提供給Simcode，應用ABAQUS進行穩(wěn)定性分析。最終得到極限強度pcr、Mises應力σmax、周向應力σ1、軸向應力σ2。

3 總體設計優(yōu)化模型的建立

3.1 近似模型及精度判斷

響應面模型(response surface methodology，RSM)是一種多項式函數(shù)，RSM可分為一階、二階、三階、四階， RSM的表達式為[15]

(15)

若所有系數(shù)均不為0，則式(15)為四階RSM，擬合所需點數(shù)為(l+1)(l+2)/2+2l；若gl=0，式(15)為三階RSM，擬合所需點數(shù)為(l+1)(l+2)/2+l；若el=0、gl=0，式(15)為二階RSM，擬合所需點數(shù)為(l+1)(l+2)/2。

徑向基神經(jīng)網(wǎng)絡模型(radial basis function, RBF)由輸入層、隱含層和輸出層組成，第q個隱藏層單元輸出的響應為

(16)

式中，cq、σq分別為第q個隱藏層的中心、單元實際的寬度。

得到輸出層中第r個輸出

(17)

式中，wqr為隱藏層q節(jié)點的第r個輸出值所占的權重;f(x)為神經(jīng)網(wǎng)絡的函數(shù)值。

(18)

(19)

(20)

式中，λi為未知的待定加權系數(shù),λi需符合式(19)無偏估計和式(20)的方差結果;γ(xs，xt)為xs和xt之間距離為h的情況下，參數(shù)的半方差值大小；γ(xs，x0)為xs和x0之間距離為h情況下，參數(shù)的半方差值大小。

利用復相關系數(shù)R2判斷近似模型的擬合精度，它與1的接近程度體現(xiàn)了近似模型擬合精度的高低，R2的計算公式為

(21)

3.2 結構學科

結構學科共有4個設計變量，若采用4階響應面模型進行擬合，則所需樣本點至少為23個。采用Opt LHD選擇50個樣本點進行結構參數(shù)化分析，結構質(zhì)量、浮容積與設計變量有固定的表達式，只需計算pcr、σmax、σ1、σ2，部分樣本點見表5。

表5 結構學科樣本點Tab.5 Sample points of structural discipline

分別采用RBF、Kriging、二階RSM、三階RSM、四階RSM對表5樣本點進行擬合，擬合精度見表6。

表6 結構近似模型擬合精度Tab.6 Fitting accuracy of structure approximate model

表6中，四階RSM的擬合精度最高，最適合用于建立結構學科近似模型。針對四階RSM選取10個樣本點的預測值，與有限元分析得到的計算值進行對比，如圖8所示。

(a)極限載荷pcr (b)最大Mises應力σmax

(c)最大周向應力σ1 (d)最大軸向應力σ2圖8 結構預測值與計算值Fig.8 Predicted values and calculated valuesof structure

圖8中，所有點均十分接近地分布于直線兩側，說明預測值和計算值很接近，進一步說明四階RSM具有較高的擬合精度。模型系數(shù)見表7。

將表7中的系數(shù)代入式(14)可得到σmax、σ1、σ2、pcr，結構的質(zhì)量mjg、浮容積Vjg根據(jù)公式計算得到，則結構學科的相關狀態(tài)變量為

(22)

式中，mother、Vother分別為耐壓殼體外附屬結構質(zhì)量和體積，兩者均為常量。

表7 結構學科模型系數(shù)Tab.7 Model coefficients of structure discipline

3.3 艇型學科

艇型學科共有設計變量5個，若采用4階響應面模型進行擬合，則需樣本點至少為31個，利用Opt LHD選擇49個艇型樣本點進行分析，見表8。

表8 艇型學科樣本點Tab.8 Sample points of boat shape discipline

分別采用RBF、Kriging、二階RSM、三階RSM、四階RSM對表8的樣本點進行擬合，擬合精度見表9。

表9 艇型學科近似模型擬合精度Tab.9 Fitting accuracies of boat disciplineapproximate model

表9中，四階RSM的擬合精度最高，最適合用于建立艇型學科近似模型。選取四階RSM的10個樣本點的預測值與計算值進行對比，見圖9。

圖9 艇型學科預測值與計算值Fig.9 Predicted values and calculated values of boatshape discipline

圖9中，所有列的點均分布于直線兩邊，且十分接近，表明預測值和實際值之間的偏差很小。艇型學科響應面模型系數(shù)見表10。

表10 艇型學科模型系數(shù)Tab.10 Boat shape discipline model coefficients

將表10中的系數(shù)代入式(15)可得到

Fx=Fx(y1,y2,y3,y4,y5)

(23)

3.4 推進學科

根據(jù)Model推進器性能參數(shù)，以推力Fp為設計變量、輸入功率Pp、質(zhì)量mp為響應，采用響應面模型擬合得到Fp與Pp、mp近似模型為

(24)

推進學科Pp、mp、Vp的R2分別為0.984、0.969、0.991，具有較高的擬合精度。預測值與計算值的對比見圖10。

(a)功率PP (b)質(zhì)量mP

(c)浮容積VP圖10 推進學科近似模型預測值與計算值Fig.10 Predicted values and calculated values ofapproximate model of propulsion discipline

圖10進一步驗證了近似模型擬合精度較高。在推力輸出的過程中，推力損失不可避免，推進系統(tǒng)推進系數(shù)Cp為

Cp=Pp/PS=ηhη0ηrηs

(25)

式中，PS為主機功率；ηh為艇身效率，取0.75；η0為敞水效率，取0.7；ηr為相對旋轉(zhuǎn)效率，在1.0～1.3之間，取1.0；ηs為軸系傳遞效率，在0.95～0.99之間，取0.95。

推進器與潛水器匹配時，需考慮推力減額CT，CT取值在1.05～1.15之間，即

1.05≤CT=FP/Fx≤1.15

(26)

取CT=1.1，則式(25)、式(26)聯(lián)立得

Fp=2.2055Fx=2.2055Fx(y1,y2,y3,y4,y5)

(27)

推進系統(tǒng)所需電量為

Qp=Ppt1

(28)

最終得到推進學科狀態(tài)變量為Qp(y1,y2,y3,y4,y5)、mp(y1,y2,y3,y4,y5)、Vp(y1,y2,y3,y4,y5)。

3.5 機電設備學科

機電設備學科用電量Qjd、質(zhì)量mjd和浮容積Vjd、重心和浮心坐標等通過下式得到：

(29)

潛水器每次執(zhí)行任務設備的工作時間基本相差不大，同時考慮潛水器應急條件下的電力需求，認為機電設備學科Qjd、mjd、Vjd、Pjd均為定值。為提高分析效率，對設備進行布置，設備重心、浮心的縱坐標位置固定、橫坐標為0。

3.6 能源學科

能源分為設備、動力和應急三部分，采用鋰電池，按電池放電90%計算。動力電源主要為推進系統(tǒng)供電。設備用電為機電設備學科供電，應急用電在潛水器緊急情況下，設備用電無法工作時，為潛水器供電。設備用電Qjd和動力用電Qp已經(jīng)得到，應急電源電量為

Qyj=Pjdt3

(30)

則能源學科狀態(tài)變量通過下式得到：

(31)

式中，Cp為推進系統(tǒng)儲備系數(shù)，取值為1.25；ρn為電池能量密度；ρz為電池質(zhì)量密度;mnyfj、Vnyfj分別為能源系統(tǒng)附件質(zhì)量和體積，均為常量。

3.7 系統(tǒng)層

對各學科的質(zhì)量和浮容積求和，得到潛水器初步設計質(zhì)量mcb、浮容積Vcb，同時為獲得較大的儲備浮力Fm，F(xiàn)m需要滿足

(32)

根據(jù)式(32)儲備浮力的要求，需要對潛水器進行壓載，使其滿足重力、浮力平衡條件，需滿足以下關系：

(33)

式中，myz為壓載質(zhì)量；Vyz為壓載浮容積；ρyz為壓載密度。

潛水器的總質(zhì)量mz、重心通過下式得到：

(34)

式中，MX、MY、MZ分別代表質(zhì)量沿x、y、z方向靜矩；xg、yg、zg分別代表重心沿x、y、z方向的坐標。

潛水器總浮容積、浮心通過下式得到：

(35)

潛水器平衡時，還需滿足

(36)

最終Fx、mz、Fm通過下式得到:

(37)

進行設計變量的靈敏度分析，只列出影響程度前10的變量，見圖11。

(a)直航阻力Fx

(b)總質(zhì)量mz

(c)儲備浮力Fm圖11 設計變量靈敏度分析Fig.11 Sensitivity analysis of design variables

圖11中藍色表示正相關，紅色表示負相關。對Fx影響最大的依次是舯段半寬y1、舯段長y4與艉段長y5、y1的相關項；對mz影響最大的依次是耐壓結構殼體厚度x3、耐壓結構半徑x4；對Fm影響最大的是依次是x4、x3。y1與x4之差為常量0.3 m。可見，潛水器的半寬(或耐壓結構半徑)、耐壓結構殼體厚度對其性能影響較大。

4 潛水器多目標優(yōu)化

4.1 優(yōu)化算法

第二代非支配排序遺傳算法(NSGA-Ⅱ)是非支配排序遺傳算法的改進算法，NSGA-Ⅱ個體的交叉運算和變異運算利用模擬二進制交叉方法(SBX)進行操作，生成子個體的交叉運算為

(38)

(39)

根據(jù)該方法生成子個體的變異運算如下:

(40)

由于u是分區(qū)間的，故區(qū)間小擾動δq由下式得到：

δq=

(41)

其中，ηm為常數(shù)，表征變異分布系數(shù)。則

(42)

NSGA-Ⅱ求解過程見圖12。

圖12 NSGA-Ⅱ算法流程圖Fig.12 Flow chart of NSGA-Ⅱ

NSGA-Ⅱ求解步驟為[16]：①隨機產(chǎn)生初始種群P0，進行秩的賦予，隨后進行種群的非劣排序；隨后進行選擇、交叉與變異等操作完成P0篩選，得到新的種群Q0；②構造新種群Rg，排序得到非劣前端F1，F(xiàn)2，…；③對所有的Fh，按照擁擠比較操作進行排序，將其中最優(yōu)秀的N個個體組成新種群Pg+1；④對Pg+1執(zhí)行選擇、交叉和變異等操作，得到新種群Qg+1；⑤若滿足優(yōu)化條件，則終止計算程序，若不滿足，令g←g+1，并返回至②。

4.2 優(yōu)化求解

多學科可行方法(multi-disciplinary discipline feasible，MDF)只在系統(tǒng)層設置目標函數(shù)，可使用優(yōu)化器直接控制優(yōu)化設計變量、設計約束和目標函數(shù)，易于編程實現(xiàn)，且MDF總能得到系統(tǒng)分析的可行解，當優(yōu)化過程突然中斷或人為中斷時，優(yōu)化結果時刻滿足系統(tǒng)分析要求。目前，MDF已成為求解高精度復雜多學科優(yōu)化問題的重要方法[17]。根據(jù)設計參數(shù)，結合學科分析模型，建立基于MDF的潛水器多學科優(yōu)化模型。系統(tǒng)層優(yōu)化模型為

(43)

式中，DV代表設計變量。

將NSGA-Ⅱ的種群數(shù)量、進化代數(shù)、交叉可能性、交叉分布指數(shù)、變異分布指數(shù)分別設置為200、120、0.9、10、20，得到的Pareto解集如圖13所示。

在圖13的Pareto解集中選擇3個方案與初始方案進行對比，見表11。

表11 優(yōu)化方案與初始方案對比Tab.11 Comparison of optimization scheme andinitial scheme

(a)直航阻力Fx、總質(zhì)量mz

(b)總質(zhì)量mz、儲備浮力Fm

(c)直航阻力Fx、儲備浮力Fm圖13 總體優(yōu)化Pareto解集Fig.13 Pareto solution set of general optimization

優(yōu)化方案1的Fm比初始方案提高了37 132.59 N，提高幅度為69.020%；Fx和mz分別降低了65.633 N、2721.6 kg，降低幅度分別為12.709%和12.504%；優(yōu)化方案2的Fm較初始方案增加了72.41 N，提高幅度為0.135%，F(xiàn)x和mz較初始方案分別降低了158.728 N、4077.59 kg，降低幅度分別為21.669%和18.134%；優(yōu)化方案3的Fm較初始方案降低22 079.13 N，降低幅度為41.040%，F(xiàn)x和mz分別降低了229.44 N、4912.97 kg，降低幅度分別為31.322%和21.849%。上述3個方案的設計變量x4、y1、x3的變化幅度較大，與圖11的結論一致。

5 結論

(1)在對STAR-CCM+軟件和ABAQUS軟件進行了二次開發(fā)，以及實現(xiàn)了STAR-CCM+軟件、ABAQUS軟件與iSight軟件集成的基礎上，設計了艇型和耐壓結構參數(shù)化分析流程，可實現(xiàn)潛水器主體的自動建模與分析，避免了設計過程中頻繁的模型修改，提高了潛水器主體設計和分析效率。

(2)對潛水器的系統(tǒng)組成進行了分析，建立了結構設計矩陣，進行了學科分解，明確了學科之間、學科與系統(tǒng)層之間的關系；利用響應面模型建立了艇型、結構和推進等學科的近似模型，可在保證精度的前提下提高學科分析效率；分析了設計變量對目標函數(shù)的影響，確定了對潛水器性能影響較大的設計變量。

(3)建立了基于MDF的潛水器優(yōu)化模型，進行了優(yōu)化求解，在得到的Pareto解集中選擇了3個方案與初始方案進行了對比。優(yōu)化方案1、優(yōu)化方案2在3個目標函數(shù)上均得到了優(yōu)化，優(yōu)化方案3潛水器質(zhì)量和直航阻力均得到了優(yōu)化，而儲備浮力有所降低，但依然滿足設計任務書的要求，進一步驗證了所得的Pareto解集優(yōu)化效果明顯。