劉峰，姚競爭，黃潯，2

(1.哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001；2.上海航空電器有限公司，上海 201100)

球殼是載人潛器重要的耐壓結(jié)構(gòu)形式之一，隨著載人潛器的輕量化、低造價的技術(shù)發(fā)展要求[1]，對于耐壓球殼的設(shè)計提出了更高的要求。耐壓球殼的設(shè)計目標是在保證其強度與穩(wěn)定性的前提下，通過尺寸的優(yōu)化達到減小其重量、降低建造成本的目標。然而，載人潛器耐壓球殼是典型的多開孔結(jié)構(gòu)，其中，人員進出口、觀察窗等較大的開孔不僅會嚴重削弱耐壓殼體的強度和剛度，還提高了其設(shè)計的難度。圍繞著多開孔耐壓球殼的設(shè)計，眾多學(xué)者開展了研究工作，孫善萍等[2]針對開孔耐壓球殼的極限強度進行了非線性有限元分析，計算了開孔的大小及加強圍壁參數(shù)的變化對殼體極限強度的影響；劉峰[3]采用第二代非支配排序遺傳算法，針對載人艙出入孔加強結(jié)構(gòu)進行了質(zhì)量和極限強度的多目標優(yōu)化。余俊等[4]針對開孔耐壓球殼，運用多島遺傳算法和Hook-Jeeves法進行了優(yōu)化。然而，載人潛器耐壓球殼需要設(shè)置觀察窗，觀察窗所采用材料的有機玻璃為粘彈性材料，恒壓作用下一定時間內(nèi)容易產(chǎn)生蠕變變形，變形過大很可能造成窗玻璃擠出，進一步對于載人潛器的安全產(chǎn)生不利影響。楊青松等[5]對于半球形觀察窗接觸、蠕變等問題進行了試驗研究和數(shù)值模擬；Guo等[6]建立了球形觀察窗的廣義 Maxwell 模型，利用 ABAQUS 軟件進行了計算；Du等[7]采用數(shù)值法建立了接觸有限元模型，確定了觀察窗應(yīng)力破壞準則。Pranesh等[8-9]開展了觀察窗的優(yōu)化研究。Wang等[10]針對全海深載人潛水器觀察窗完整窗口模型，進行了準靜態(tài)加載-卸載循環(huán)試驗，所得出的研究結(jié)論為全海深觀察窗設(shè)計提供了依據(jù)和參考；黃潯等[11]利用有限元分析方法對于觀察窗蠕變特性進行了研究，總結(jié)了觀察窗厚度-直徑比、傾角和接觸面摩擦系數(shù)等設(shè)計參數(shù)對于觀察窗軸向位移和應(yīng)力的影響。目前針對耐壓球殼的優(yōu)化設(shè)計中，未有將觀察窗蠕變變形納入優(yōu)化范圍，僅對球殼質(zhì)量或極限強度的優(yōu)化可能忽略了觀察窗的大變形，從而帶來一系列安全隱患。

本文在考慮觀察窗蠕變特性的基礎(chǔ)上，設(shè)計了耐壓球殼參數(shù)化分析流程，利用該分析流程進行了樣本點的選取及分析，在建立了基于蠕變效應(yīng)的耐壓球殼近似模型的基礎(chǔ)上，進行了耐壓球殼的優(yōu)化的求解。

1 考慮蠕變特性的耐壓球殼分析模型

1.1 觀察窗結(jié)構(gòu)形式

某載人潛器設(shè)有2個側(cè)觀察窗，1個主觀察窗。由于各個觀察窗之間球面距離滿足《潛器規(guī)范》[12]要求，獨立且互不影響，因此，僅針對主觀察窗進行研究。耐壓球殼觀察窗和窗座結(jié)構(gòu)形式及主要設(shè)計參數(shù)見圖1。

圖1 觀察窗結(jié)構(gòu)形式Fig.1 Structure of observation window

圖1中，初始設(shè)計參數(shù)為：D2=470 mm，D1=410 mm，d=200 mm，d1=220 mm，t=105 mm，α=β=45°。窗座與殼體均采用高強度鋼材，觀察窗為有機玻璃，材料參數(shù)見表1。

表1 材料物理屬性Table 1 Physical properties of materials

1.2 觀察窗蠕變位移求解

根據(jù)試驗擬合出常溫有機玻璃回歸公式為[13]：

(1)

式中：ε為應(yīng)變；εc為臨界斷裂應(yīng)變；tc為臨界斷裂時間，εc和tc與應(yīng)力水平相關(guān)：

(2)

(3)

式中：σ為應(yīng)力；E為彈性模量。

將式(2)、(3)代入式(1)，得到觀察窗有機玻璃本構(gòu)模型為：

logε=0.036 8(logt+0.021 94σ)2+

0.278 1(logt+0.021 94σ)-0.028 95

(4)

1.3 有限元分析模型

觀察窗與窗座之間接觸定義采用接觸對算法，設(shè)定軟質(zhì)材料觀察窗為從面，高強度鋼材窗座為主面，接觸面摩擦系數(shù)為0.1。遵循從面網(wǎng)格尺寸小于主面網(wǎng)格，考慮計算量，觀察窗全局種子取10，殼體取20，網(wǎng)格劃分見圖2。分析步設(shè)為靜力、粘性，分析時間為潛器標準作業(yè)時間6 h。

圖2 模型網(wǎng)格劃分Fig.2 Model mesh generation

為消除耐壓球殼剛體位移，采用3點約束球殼6個方向自由度，具體見圖3。

圖3 三點約束模式Fig.3 Three points constraint mode

1.4 觀察窗蠕變位移與耐壓球殼強度求解

目標潛器工作水深1 500 m，根據(jù)《潛器規(guī)范》[12]，安全系數(shù)取1.5，則計算壓力Pj為24.5 MPa。耐壓球殼表面施加24.5 MPa載荷進行蠕變分析，結(jié)果見圖4、圖5。計算得到的耐壓球殼最大Mises應(yīng)力為606.63 MPa，滿足規(guī)范強度要求，觀察窗最大蠕變位移為2.914 mm。

圖4 耐壓球殼應(yīng)力云圖Fig.4 Stress nephogram of pressure spherical shell

圖5 觀察窗蠕變位移云圖Fig.5 Creep displacement nephogram of observation window

1.5 耐壓球殼穩(wěn)定性分析

基于上述有限元模型，采用弧長法進行耐壓球殼極限載荷Pcr的計算，弧長法能很好對于位移-載荷曲線下降段進行捕捉，在結(jié)構(gòu)的后屈曲分析中已經(jīng)得到了廣泛應(yīng)用。計算中的初始缺陷通過模型特征值屈曲模態(tài)的提取得到，并對于幾何以及材料非線性進行考慮。將材料的名義應(yīng)力應(yīng)變進行轉(zhuǎn)化、擬合得到真實應(yīng)力應(yīng)變曲線。替換模型分析步為屈曲和線性攝動，進行單位載荷施加，分析特征值屈曲。提取得到前六階模態(tài)，部分模態(tài)如圖6所示。

圖6 部分模態(tài)Fig.6 Partial modal diagram

圖6中，各階模態(tài)特征值相差不大，可見，模型對初始缺陷比較敏感，說明在添加初始缺陷時提取全部六階模態(tài)的節(jié)點位移來使計算結(jié)果偏于保守，模態(tài)比例因子通常取殼體厚度的1%。再次替換模型分析步為靜力、弧長，載荷根據(jù)第六階模態(tài)特征值取109.85 MPa，進行后屈曲分析，球殼失穩(wěn)破壞形式如圖7，弧長-載荷比例因子曲線見圖8。

圖7 球殼失穩(wěn)破壞Fig.7 Instability failure diagram of spherical shell

圖7、圖8計算得到的耐壓球殼極限載荷Pcr=33.933 MPa，滿足規(guī)范穩(wěn)定性要求。

圖8 弧長-載荷比例因子曲線Fig.8 Arc length load scale factor curve

2 考慮蠕變的耐壓球殼參數(shù)化分析

2.1 基于Python語言的Abaqus二次開發(fā)

命令執(zhí)行前借助Python解釋器生成.rpy格式文件提交計算將生成inp格式文件，通過分析求解生成便于后處理的odb文件。

Abaqus創(chuàng)建的模型包括根(root)對象：進行調(diào)整模型視角等的視圖模塊session；存儲模型樹中包含的部件、材料、裝配體、分析步、載荷等數(shù)據(jù)模型的數(shù)據(jù)模塊mdb和odb；存儲所有模型信息、輸出數(shù)據(jù)等的數(shù)據(jù)輸入輸出模塊odb。每個根對象下面又有很多命令分支，通過編寫腳本命令可以準確地對不同對象進行訪問和編輯。

2.2 Abaqus的集成

在iSight軟件進行Abaqus軟件集成，主要包括以下階段：

1)前處理階段(模型參數(shù)化)。對所研究問題的模型進行定義，包括模型的幾何尺寸、材料參數(shù)、載荷邊界等進行參數(shù)化處理，得到Abaqus輸入文件(.inp,.py)。最為簡便直接的是操作Abaqus/CAE建立三維模型，對于模型錄制.rpy格式文件并修改后綴為.py的Python腳本文件，或直接編寫Python的代碼作為Abaqus輸入文件；

2)模擬計算階段。在iSight中使用.bat格式文件調(diào)用Abaqus/Standard或Abaqus/Explict求解器，后臺自動計算得到輸出文件(job.odb、job.dat等)；

3)后處理。在計算結(jié)果文件得到后，進一步通過Abaqus或其他后處理軟件進行結(jié)果數(shù)據(jù)讀取，或者通過iSight讀取輸出文件中輸出的數(shù)據(jù)，進行結(jié)果分析或建立近似模型。

2.3 參數(shù)化分析的實現(xiàn)

iSight軟件提供了與Abaqus軟件的接口，在iSight軟件進行Abaqus軟件集成主要通過文件進行連接，具體步驟及方法為：

1)蠕變計算。輸入.py文件。文件前處理部分在Abaqus軟件用戶界面上操作生成，包含有限元模型和計算設(shè)置信息等。數(shù)據(jù)提取部分主要提取蠕變計算結(jié)果中模型總的最大Mises應(yīng)力、模型質(zhì)量以及目標位置的位移等信息；采用.bat批處理文件調(diào)用Abaqus進行后臺計算；輸出文件為包含模型所有節(jié)點Mises應(yīng)力信息和最大、最小Mises應(yīng)力數(shù)據(jù)的.rpt報告文件；以及包含由.py文件取到的目標位置中的蠕變位移和模型質(zhì)量的.txt文本文件。

2)線性屈曲分析。輸入文件為.py文件。在得到.rpy文件后，修改文件后綴為.py。文件中的前處理部分在前一步蠕變分析的基礎(chǔ)上，對關(guān)鍵詞、替換分析步設(shè)置等進行修改，輸出部分為分析結(jié)果文件中一階模態(tài)的節(jié)點位移信息，作為模型初始缺陷數(shù)據(jù)；采用.bat批處理文件調(diào)用Abaqus進行后臺計算；輸出文件為包含了線性屈曲計算結(jié)果中所得到的各階模態(tài)特征值的.dat文件，特征值與載荷的乘積即為隨后進行的非線性分析中的載荷值。

3)非線性屈曲分析。對于.py、.rpy等格式的文件操作以及修改、設(shè)置等操作與第2步相同。輸出部分主要涉及極限載荷的結(jié)果，輸出模型弧長-載荷因子曲線數(shù)據(jù)信息，提取最大載荷因子；采用.bat批處理文件調(diào)用Abaqus進行后臺計算；輸出包含模型載荷因子-弧長曲線數(shù)據(jù)信息的.rpt報告文件。

將上述輸入、批處理、輸出文件寫入iSight軟件Simcode模塊中，通過iSight對于Abaqus軟件進行驅(qū)動，實現(xiàn)模型參數(shù)化并依次進行強度分析、蠕變分析、線性屈曲分析及非線性分析。圖1所示的耐壓球殼觀察窗結(jié)構(gòu)幾何參數(shù)，即設(shè)計變量有α、β、t、D2、D1、d、d1，輸出為耐壓球殼極限載荷壓力Pcr、觀察窗臨界應(yīng)力σw、球殼臨界應(yīng)力σb、質(zhì)量M以及觀察窗蠕變位移U。上述設(shè)計變量中，D2、D1、d、d1由球殼布置決定，參數(shù)α、t會影響β，也會使U、Pcr和M發(fā)生變化，因此，確定α、t作為設(shè)計參數(shù)進行參數(shù)化，取值范圍分別為：38°≤α≤60°，105 mm≤t≤150 mm。

3 耐壓球殼優(yōu)化求解

3.1 近似模型的建立

神經(jīng)網(wǎng)絡(luò)具有強大的處理復(fù)雜非線性問題的能力、高容錯度和無須數(shù)學(xué)假設(shè)等優(yōu)勢，已成為解決許多領(lǐng)域中的復(fù)雜優(yōu)化問題的有力工具。神經(jīng)網(wǎng)絡(luò)由輸入層、中間層和輸出層3部分構(gòu)成，通過非線性變換完成輸入矢量到新空間之間的映射，即中間層空間。中間層空間到輸出層空間則為線性映射，輸出層在新的線性空間中進行線性加權(quán)組合。徑向基(radial basis functions,RBF)網(wǎng)絡(luò)函數(shù)為[14]：

(5)

式中：m為樣本點數(shù)；x為設(shè)計變量向量；λi(i=1,2,…，m)為待定加權(quán)系數(shù)；φ(r0i)為基函數(shù)；r0i=‖x-xi‖為插值函數(shù)中任意一點到第i個插值點的歐幾里得距離；Pi(x)為多項式項；K為多項式項個數(shù)；ci(i=1,2,…，K)為Pi(x)所對應(yīng)系數(shù)。

基函數(shù)一般選用高斯函數(shù)：

(6)

響應(yīng)面模型(response surface methodology，RSM)是一種對于試驗樣本點進行擬合的多項式函數(shù)。RSM的表達式為：

(7)

誤差分析采用工程領(lǐng)域常用的R2分析，R2取值范圍為[0,1]，越接近1則表明模型擬合誤差越小，計算公式為：

(8)

不同設(shè)計變量取值范圍存在差異，量級和量綱的不同往往會降低優(yōu)化效率，甚至影響優(yōu)化結(jié)果。將設(shè)計變量歸一化，使各設(shè)計變量均位于[-1,1]，建立起以原點為中心的立方形設(shè)計空間，便于統(tǒng)一處理，有利于解決量綱差異帶來的麻煩。本文所選取的設(shè)計變量歸一化處理為：

(9)

構(gòu)建近似模型需要通過試驗設(shè)計獲取一定數(shù)量的樣本點，最優(yōu)拉丁超立方設(shè)計(optimal Latin hypercube design,OptLHD)是基于拉丁方設(shè)計(Latin hypercube design,LHD)的改良版，其能使所有試驗點盡可能地均勻分布于采樣空間，良好的填充性和均衡性使得該方法更加精準真實[15]。選取30個樣本點進行計算，分別采用二階響應(yīng)面模型、四階響應(yīng)面模型和徑向基神經(jīng)網(wǎng)絡(luò)模型對樣本點進行擬合，不同近似模型擬合精度見表2。

表2 不同近似模型擬合精度Table 2 Fitting precision of different approximate models

表2中，徑向基神經(jīng)網(wǎng)絡(luò)模型具有很高的精度，滿足工程實際意義，因此，選取徑向基神經(jīng)網(wǎng)絡(luò)模型建立耐壓球殼的近似模型。

3.2 單目標優(yōu)化

以觀察窗傾角和厚度為設(shè)計變量，耐壓球殼質(zhì)量最小、觀察窗蠕變位移最小為目標函數(shù)，殼體臨界應(yīng)力σ小于0.85倍的材料屈服強度σs、極限載荷Pcr大于計算壓力Pj為約束條件，建立耐壓球殼單目標優(yōu)化模型為：

(10)

多島遺傳算法(multi island genetic algorithm，MIGA)是一種全局優(yōu)化算法，基于傳統(tǒng)遺傳算法對并行分布遺傳算法有所改進。MIGA將種群分成多個子群，又稱之為“島”，每個島上子群獨立進化，而非全部種群采用相同的進化機制，一定時間間隔進行“島嶼”間的“遷移”，完成信息交換[16]。MIGA能夠有效提高運算速度，并且若干獨立進化的子群豐富了整個種群的遺傳多樣性，避免落入局部解。圖9為算法流程圖。

圖9 多島遺傳算法流程Fig.9 Multi Island genetic algorithm process

采用MIGA，子群規(guī)模為10，島數(shù)為20，代數(shù)為10，交叉概率1.0，變異概率0.01，計算結(jié)果與有限元計算結(jié)果對比見表3。

表3 優(yōu)化結(jié)果與有限元結(jié)果對比Table 3 Comparison between optimization results and finite element results

基于近似模型的優(yōu)化結(jié)果與有限元結(jié)果對比，質(zhì)量相對誤差為0.13%，位移相對誤差為2.25%，說明近似模型具有較高的精度。優(yōu)化結(jié)果較初始設(shè)計質(zhì)量減小0.25%，而蠕變位移增大6.42%，說明不考慮觀察窗蠕變效應(yīng)時，質(zhì)量優(yōu)化效果不明顯，而蠕變變形明顯增大，因此，需要將觀察窗蠕變效應(yīng)作為目標函數(shù)，進一步開展多目標優(yōu)化設(shè)計。

3.3 多目標優(yōu)化

以觀察窗傾角和厚度為設(shè)計變量，耐壓球殼質(zhì)量最小、觀察窗蠕變位移最小為目標函數(shù)，殼體臨界應(yīng)力σ小于0.85倍的材料屈服強度σs、極限載荷Pcr大于計算壓力Pj為約束條件，建立耐壓球殼單目標和多目標優(yōu)化模型為：

(11)

第2代非支配排序遺傳算法NSGA-II[15]改進了非支配排序方法，選擇接近Pareto前沿的個體，增強了Pateto的前進能力，擁有良好的探索能力。導(dǎo)入了“擁擠距離”和“擁擠距離排序”的方法，在具有同樣的Pareto順序的層內(nèi)，能夠?qū)€體進行排序，稱為擁擠距離排序。進化過程中，當前父代群體通過交叉和變異生成子群體，將2個群體合并。在目標空間中遵照Pareto最優(yōu)關(guān)系將群體中個體兩兩按其目標函數(shù)向量進行比較，將群體中所有個體分成多個依次控制的前沿層。在屬于不同的Pareto層的情況下，利用評價Pareto優(yōu)越性來評價個體的優(yōu)劣。屬于相同Pareto層的個體，具有更大擁擠距離的個體更優(yōu)秀。交叉和變異運算機制使用SBX方法。根據(jù)該方法生成子個體的交叉運算為：

(12)

(13)

式中：u在[0,1]中隨機取值；ηc為常數(shù)，表征交叉分布系數(shù)，大小直接決定了子代個體與父代個體的接近程度。

根據(jù)該方法生成子個體的變異運算為：

(14)

(15)

(16)

采用NSGA-II算法，種群數(shù)為100、代數(shù)為50、交叉概率0.8、變異概率0.1。計算得出Pareto解集如圖10。圖中，U與球殼M基本呈負相關(guān)，Pareto解集分布均勻，在耐壓球殼觀察窗設(shè)計過程中，設(shè)計者可根據(jù)實際情況進行設(shè)計方案的選擇。選擇若干優(yōu)化方案與初始設(shè)計方案對比，列于表4。

圖10 多目標優(yōu)化Pareto解集Fig.10 Pareto solution set of multiobjective optimization

表4中，各優(yōu)化方案與初始設(shè)計方案相比，球殼質(zhì)量稍有增大，但增幅不大，最大僅為2.07%，而觀察窗蠕變變形有很大的改善，方案2減小最為明顯，達到-18.09%。綜合來看，在載人潛器耐壓球殼優(yōu)化設(shè)計時，考慮觀察窗蠕變效應(yīng)，能大幅減小蠕變變形，減小潛器的安全隱患，使其安全性能顯著提高，而對球殼質(zhì)量增加不大。

表4 優(yōu)化方案與初始設(shè)計方案對比Table 4 Comparison between optimization scheme and initial design scheme

4 結(jié)論

1)通過對于Abaqus軟件進行二次開發(fā)和集成，實現(xiàn)了考慮蠕變特性耐壓球殼的參數(shù)化建模與自動分析，提高了耐壓球殼設(shè)計效率。

2)得到了滿足工程需要的近似模型，統(tǒng)一了設(shè)計變量的量級和量綱，計算表明合理選擇近似模型及處理設(shè)計變量能夠提高耐壓結(jié)構(gòu)設(shè)計和優(yōu)化效率。

3)考慮蠕變效應(yīng)與未考慮時相比，雖然耐壓球殼質(zhì)量有所增加，但觀察窗蠕變變形減小明顯，說明了在載人潛器耐壓球殼設(shè)計優(yōu)化過程中，應(yīng)加以考慮觀察窗蠕變效應(yīng)，有利于改善觀察窗變形不協(xié)調(diào)問題，從而提高載人潛器的安全性。