基于響應面-蒙特卡羅法的實體靶船極限強度分析

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

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

0 引 言

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

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

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

1 基本原理與方法

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

工程結構的失效概率可以表示為[6]：

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

用蒙特卡羅法計算結構可靠性時，式（1）可表達為：

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

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

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

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

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

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

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

板架的整體變形

板的局部變形

骨材的側傾變形

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

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

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

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

2.1 計算模型

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

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

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

2.2 隨機變量的概率分布

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

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

表1 腐蝕參數的統計特征值Tab. 1 Statistical characteristics of corrosion parameters

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

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

2.3 極限彎矩計算與響應面構建

采用顯式動態算法分析該艦中垂狀態的極限強度，圖3為艙段的彎矩-轉角曲線，圖4為極限狀態結構應力云圖?？梢园l現，該艦的中垂極限彎矩為36.21 MN·m，極限狀態該艦強橫梁間的甲板板架發生了屈曲破壞。

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

圖3 彎矩-轉角曲線Fig. 3 Bending moment vs. angle of rotation relationship

圖4 極限狀態艙段應力云圖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

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

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

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

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

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

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

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

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

2.4 極限強度的可靠性分析

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

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

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

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

可以發現：極限強度的隨機分布與Weibull分布偏差較大。因此，實體靶船極限強度的概率分布更趨近于正態分布或對數正態分布。

3 結 語

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