叢 濱，董龍昌，蔡慶港，張厚堯，楊 衡，李陳峰,3

(1. 中國人民解放軍92941部隊，遼寧 葫蘆島 125000；2. 哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱150001；3. 哈爾濱工程大學 青島船舶科技有限公司，山東 青島 266400)

0 引 言

我國部分水面艦船由于服役海區(qū)海況條件較好，實際航行時間也低于設計要求，同時加上維護保養(yǎng)得當，艦船平臺狀態(tài)良好，有很大的潛力超期服役[1]。此類超役艦船一部分會被改裝為實體靶船，用以部隊訓練、試驗等科研活動。極限強度是艦船結(jié)構(gòu)強度的重要方面，也是實體靶船生命力的重要保障，關(guān)系到靶船航行安全及試訓活動組織。目前，國軍標等相關(guān)水面艦艇規(guī)范及標準缺少實體靶船極限強度評估的明確方法[2]。

船體的極限承載能力一般指總縱彎曲作用下船體梁的極限彎矩，常用的計算方法有解析法、簡化逐步破壞法和非線性有限元法等[3]。其中，非線性有限元法通過建立結(jié)構(gòu)有限元模型，模擬結(jié)構(gòu)的初始缺陷，并考慮材料非線性和幾何非線性，計算精度高，被國際船舶與海洋工程結(jié)構(gòu)大會（ISSC）作為極限承載能力預報的衡準方法。如果考慮材料機械性能、結(jié)構(gòu)尺寸和腐蝕等參數(shù)的不確定性，極限承載能力為符合一定概率分布的隨機變量[4]。響應面法結(jié)合蒙特卡羅法隨機抽樣是一種求解隨機變量統(tǒng)計特征值的較好方法，該方法先通過有限元數(shù)值計算并采用最小二乘法等方法擬合一個響應面來替代未知的、真實的極限狀態(tài)方程，再根據(jù)隨機變量的概率分布運用蒙特卡羅法抽樣獲得目標變量的概率分布，二者的結(jié)合使結(jié)構(gòu)分析的工作量大大減少，克服了傳統(tǒng)蒙特卡羅法龐大的有限元計算工作量和傳統(tǒng)響應面因為線性化帶來誤差的缺點[5]。

為了合理地評估實體靶船的承載能力，考慮材料機械性能和結(jié)構(gòu)腐蝕等參數(shù)的不確定性，本文將采用響應面-蒙特卡羅法結(jié)合非線性有限元法研究建立一種實體靶船極限強度的計算方法，分析響應面設計點的選取方案，給出實體靶船極限承載能力的概率分布。

1 基本原理與方法

1.1 可靠性分析的響應面-蒙特卡羅法

工程結(jié)構(gòu)的失效概率可以表示為[6]：

式中：G(X)是極限狀態(tài)方程，即響應面方程，當G(X)＜0，結(jié)構(gòu)失效；Df是與G(X)相對應的失效區(qū)域；f(X)=f(x1,x2,···,xn)是 基本隨機變量X的聯(lián)合密度函數(shù)，當X為一組相互獨立的隨機變量時，有f(X)=

用蒙特卡羅法計算結(jié)構(gòu)可靠性時，式（1）可表達為：

式中：N為抽樣模擬總數(shù)；當用相對誤差來表示蒙特卡羅法的抽樣誤差有：

式中： 是待定系數(shù)；和 分別是基本隨機變量。

為使響應面函數(shù)更好地逼近真實狀態(tài)，設計點選取非常重要。一方面盡量要少，以降低計算工作量，另一方面還要盡量包含較多極限狀態(tài)信息，以保證響應面的擬合精度。

響應面-蒙特卡羅法[7]首先根據(jù)隨機變量的概率分布選取確定功能函數(shù)值的設計點。設計點一般以隨機變量的均值點為中心，根據(jù) 3 σ 原則在 μi±kσi范圍內(nèi)選取。進而，根據(jù)設計點方案運用有限元數(shù)值計算得到設計點對應的功能函數(shù)值。根據(jù)有限元計算獲得的一系列功能函數(shù)值，采用最小二乘法等方法擬合響應面，來替代真實的極限狀態(tài)方程。在此基礎上，運用蒙特卡羅法對隨機變量進行抽樣，采用響應面進行N次數(shù)值計算，得到N個功能函數(shù)值Zj(j=1,2,···,N)，進而確定功能函數(shù)的概率分布，或根據(jù)統(tǒng)計Zj＜0的個數(shù)，計算失效概率或可靠度指標。因此，響應面與蒙特卡羅法結(jié)合使結(jié)構(gòu)分析的工作量大大減少，提高效率。

1.2 極限承載能力分析的非線性有限元法

基于非線性有限元法的船體極限強度計算精度主要取決于模型構(gòu)建的準確性和非線性求解器的選擇與參數(shù)設置[8]。

為了避免邊界條件的影響，模型范圍一般選擇艙段模型。板和主要縱向構(gòu)件應采用4節(jié)點單元，且盡量保證網(wǎng)格形狀接近正方形。為了準確模擬極限狀態(tài)結(jié)構(gòu)的塑性大變形和失效模式演化，骨材腹板一般要求至少布置3～4個單元，并以此為基準劃分整個艙段模型的網(wǎng)格往往可以取得較好的效果。同時，極限強度計算需要考慮結(jié)構(gòu)的初始撓度[9]，包括：

板架的整體變形

板的局部變形

骨材的側(cè)傾變形

式中：a和b分別為橫向和縱向構(gòu)件的間距；S為縱桁間距；hw為 骨材腹板高度；m為 半波數(shù)，取滿足a/b≤的 最小正整數(shù)；A0，B0和C0為初始撓度的幅值，分別取

由于極限強度計算需要考慮材料的非線性，在缺少船用鋼材料拉伸曲線的情況下，一般采用理想彈塑性或理想塑性強化模型以考慮材料的非線性效應。

目前，極限強度分析的非線性求解方法主要有Newton-Raphson迭代法、Riks法和顯式動態(tài)算法等。其中，Newton-Raphson法以載荷增量為控制量，根據(jù)節(jié)點力平衡判斷收斂，但結(jié)構(gòu)大變形引起的節(jié)點力不平衡容易造成不收斂。Riks法采用弧長作為控制變量，改進了Newton-Raphson法的迭代策略，計算效率和收斂性主要取決于弧長相關(guān)參數(shù)的設置。顯式動態(tài)算法采用中心差分法，用上一步和當前步的結(jié)果計算下一步的計算結(jié)果，無需迭代運算，收斂性好，計算效率較高，目前被廣泛地應用于結(jié)構(gòu)靜態(tài)和準靜態(tài)的非線性響應分析。

2 實體靶船極限強度可靠性分析

2.1 計算模型

以某型超役水面艦艇為例，該船已服役近30年，將改裝成為實體靶船?？紤]材料屈服限和結(jié)構(gòu)腐蝕的隨機變化，采用響應面-蒙特卡羅法結(jié)合非線性有限元法，開展了極限強度可靠性分析。

選取船中艙段為分析對象，采用Abaqus通用有限元軟件建立有限元模型，網(wǎng)格尺寸整體上取50 mm×50 mm，同時縱骨腹板劃分3個單元，單元類型選用4節(jié)點薄殼單元S4R，如圖1所示。該艦材料屈服限為235 MPa，楊氏模量為206 GPa，泊松比0.3

圖1 艙段有限元模型Fig. 1 FE model of target cabin

2.2 隨機變量的概率分布

由于是實體靶船，因此結(jié)構(gòu)腐蝕是重點考慮的一個隨機變量，同時由于材料屈服限對于艦船極限強度影響較大，因此主要考慮腐蝕量和材料屈服限的隨機變化。材料屈服限服從正態(tài)分布，變異系數(shù)取0.05，即σy～N(235,138.06)[4]。

國內(nèi)外學者通過對結(jié)構(gòu)腐蝕統(tǒng)計分析和試驗研究，提出了一些腐蝕模型[10]。秦圣平等[11]在現(xiàn)有腐蝕模型比較分析的基礎上，提出了一種適用于船體結(jié)構(gòu)時變可靠性分析的非線性腐蝕模型，該模型可以較好模擬鋼結(jié)構(gòu)在腐蝕環(huán)境下的腐蝕損傷過程。表達式為：式中：d∞，β，η，Tst為待定參數(shù)，均服從正態(tài)分布，其統(tǒng)計特征值見表1。

表1 腐蝕參數(shù)的統(tǒng)計特征值Tab. 1 Statistical characteristics of corrosion parameters

將表中相關(guān)參數(shù)代入腐蝕模型，采用蒙特卡羅抽樣，抽樣結(jié)果如圖2所示，獲得了服役30年的腐蝕量服從正態(tài)分布，d(30)～N(1.67,4.36E-3)。

圖2 腐蝕量的概率分布Fig. 2 Probability distribution of structural corrosion

2.3 極限彎矩計算與響應面構(gòu)建

采用顯式動態(tài)算法分析該艦中垂狀態(tài)的極限強度，圖3為艙段的彎矩-轉(zhuǎn)角曲線，圖4為極限狀態(tài)結(jié)構(gòu)應力云圖?？梢园l(fā)現(xiàn)，該艦的中垂極限彎矩為36.21 MN·m，極限狀態(tài)該艦強橫梁間的甲板板架發(fā)生了屈曲破壞。

根據(jù)材料屈服限和結(jié)構(gòu)腐蝕量的概率分布及其統(tǒng)計特征值，采用 3 σ原則確定49個響應面設計點，采用非線性有限元法計算獲得了不同腐蝕量和材料屈服限組合下的極限強度，如表2所示。

圖3 彎矩-轉(zhuǎn)角曲線Fig. 3 Bending moment vs. angle of rotation relationship

圖4 極限狀態(tài)艙段應力云圖Fig. 4 Stress patterns of cabin at ultimate strength level

表2 不同材料屈服限和腐蝕量下的艙段極限強度（MN·m）Tab. 2 Ultimate strength of cabin under different yield strength and corrosion amount

從艙段極限強度的計算結(jié)果可以發(fā)現(xiàn)：隨著材料屈服限的提高，極限強度逐步增加；隨著結(jié)構(gòu)腐蝕量的增加，極限強度不斷降低?？傮w而言，材料屈服限對極限強度的影響較腐蝕量的影響要大。

在此基礎上，采用三次多項式表達式對響應面進行擬合，設計了10種不同設計點方案，分別包括49，25，21，21，17，17，17，13，13和13個設計點，如圖5所示。表3為相關(guān)設計點方案的響應面擬合結(jié)果，圖6為部分設計點方案擬合得到的響應面。為了評價響應面擬合效果，引入了局部擬合優(yōu)度和整體擬合優(yōu)度兩個評價指標。其中，局部擬合優(yōu)度反映當前設計點方案n個設計點與該方案下擬合得到的響應面之間的擬合優(yōu)度；總體擬合優(yōu)度指全部49個設計點與n個設計點方案擬合得到的響應面之間的擬合優(yōu)度，表達式為：

圖5 響應面設計點方案Fig. 5 Design point schemes for response surface

表3 響應面擬合結(jié)果Tab. 3 Fitting results of response surfaces

圖6 不同設計點方案的響應面擬合結(jié)果Fig. 6 Response surfaces of different design point scheme

式中：zi為設計點的有限元計算結(jié)果，即原始數(shù)據(jù)；為根據(jù)響應面函數(shù)得到設計點的極限強度；表示原始數(shù)據(jù)的均值，即49個有限元計算結(jié)果的均值。

從響應面的擬合結(jié)果可以發(fā)現(xiàn)：前9個方案的局部擬合效果均大于0.99；同時，前8個方案的整體擬合優(yōu)度均大于0.99，且方案1、方案4、方案6和方案8達到了0.999 8，對應的設計點數(shù)量分別為49，21，17和13，說明響應面的整體擬合優(yōu)度與設計點的選取息息相關(guān)，并不是設計點數(shù)量越多越好。方案10的設計點選取沒有考慮材料屈服限和腐蝕量的相關(guān)性，擬合得到的響應面已失真，如圖6（d）所示，局部擬合優(yōu)度僅0.482。方案8雖然只選取了13個設計點，但其整體擬合優(yōu)度與49個設計點一致，達到了0.999 8，局部擬合優(yōu)度甚至達到了1.0，擬合優(yōu)度比方案2（25個設計點）和方案3（21個設計點）都要好。

因此，響應面設計點的合理選取，不僅可以有效降低計算工作量，同時可以保證響應面擬合精度。

2.4 極限強度的可靠性分析

獲得極限強度的響應面后，結(jié)合相關(guān)參數(shù)的概率分布，采用蒙特卡羅法隨機抽樣，進一步分析實體靶船極限強度的概率分布。

首先，開展了蒙特卡羅法的收斂性分析，抽樣次數(shù)分別取1 000，10 000，100 000和1 000 000次，樣本的統(tǒng)計特征值見表5?？梢园l(fā)現(xiàn)，當樣本數(shù)量達到10 000次，極限強度的均值與標準差已趨于穩(wěn)定；當樣本數(shù)量增加到100 000次或更高，極限強度的統(tǒng)計特征值變化不大。在正態(tài)分布的基礎上，進一步采用對數(shù)正態(tài)分布和Weibull分布擬合極限強度的概率分布，擬合結(jié)果如圖7所示。其中，正態(tài)分布和對數(shù)正態(tài)分布的分布函數(shù) 分 別 為和 l nWeibull分布的形狀參數(shù)為37.16，尺度參數(shù)為21.47，分布函數(shù)為

表4 樣本的均值與標準差Tab. 4 Mean and standard deviation of samples

圖7 極限強度的概率分布Fig. 7 Probability distribution of ultimate strength

可以發(fā)現(xiàn)：極限強度的隨機分布與Weibull分布偏差較大。因此，實體靶船極限強度的概率分布更趨近于正態(tài)分布或?qū)?shù)正態(tài)分布。

3 結(jié) 語

本文采用蒙特卡羅法與響應面法相結(jié)合的可靠性分析方法，開展考慮材料屈服限和腐蝕影響的艙段極限強度研究，研究表明：在保證擬合優(yōu)度的前提下，響應面設計點的合理選取，可以有效降低計算工作量；考慮材料屈服限和腐蝕聯(lián)合影響的實體靶船極限強度概率分布更趨近于正態(tài)分布或?qū)?shù)正態(tài)分布。本文研究建立的基于響應面-蒙特卡羅法的可靠性分析方法，對于實體靶船極限強度評估具有指導意義。