叢 濱,董龍昌,蔡慶港,張厚堯,楊 衡,李陳峰,3
(1. 中國(guó)人民解放軍92941部隊(duì),遼寧 葫蘆島 125000;2. 哈爾濱工程大學(xué) 船舶工程學(xué)院,黑龍江 哈爾濱150001;3. 哈爾濱工程大學(xué) 青島船舶科技有限公司,山東 青島 266400)
我國(guó)部分水面艦船由于服役海區(qū)海況條件較好,實(shí)際航行時(shí)間也低于設(shè)計(jì)要求,同時(shí)加上維護(hù)保養(yǎng)得當(dāng),艦船平臺(tái)狀態(tài)良好,有很大的潛力超期服役[1]。此類超役艦船一部分會(huì)被改裝為實(shí)體靶船,用以部隊(duì)訓(xùn)練、試驗(yàn)等科研活動(dòng)。極限強(qiáng)度是艦船結(jié)構(gòu)強(qiáng)度的重要方面,也是實(shí)體靶船生命力的重要保障,關(guān)系到靶船航行安全及試訓(xùn)活動(dòng)組織。目前,國(guó)軍標(biāo)等相關(guān)水面艦艇規(guī)范及標(biāo)準(zhǔn)缺少實(shí)體靶船極限強(qiáng)度評(píng)估的明確方法[2]。
船體的極限承載能力一般指總縱彎曲作用下船體梁的極限彎矩,常用的計(jì)算方法有解析法、簡(jiǎn)化逐步破壞法和非線性有限元法等[3]。其中,非線性有限元法通過(guò)建立結(jié)構(gòu)有限元模型,模擬結(jié)構(gòu)的初始缺陷,并考慮材料非線性和幾何非線性,計(jì)算精度高,被國(guó)際船舶與海洋工程結(jié)構(gòu)大會(huì)(ISSC)作為極限承載能力預(yù)報(bào)的衡準(zhǔn)方法。如果考慮材料機(jī)械性能、結(jié)構(gòu)尺寸和腐蝕等參數(shù)的不確定性,極限承載能力為符合一定概率分布的隨機(jī)變量[4]。響應(yīng)面法結(jié)合蒙特卡羅法隨機(jī)抽樣是一種求解隨機(jī)變量統(tǒng)計(jì)特征值的較好方法,該方法先通過(guò)有限元數(shù)值計(jì)算并采用最小二乘法等方法擬合一個(gè)響應(yīng)面來(lái)替代未知的、真實(shí)的極限狀態(tài)方程,再根據(jù)隨機(jī)變量的概率分布運(yùn)用蒙特卡羅法抽樣獲得目標(biāo)變量的概率分布,二者的結(jié)合使結(jié)構(gòu)分析的工作量大大減少,克服了傳統(tǒng)蒙特卡羅法龐大的有限元計(jì)算工作量和傳統(tǒng)響應(yīng)面因?yàn)榫€性化帶來(lái)誤差的缺點(diǎn)[5]。
為了合理地評(píng)估實(shí)體靶船的承載能力,考慮材料機(jī)械性能和結(jié)構(gòu)腐蝕等參數(shù)的不確定性,本文將采用響應(yīng)面-蒙特卡羅法結(jié)合非線性有限元法研究建立一種實(shí)體靶船極限強(qiáng)度的計(jì)算方法,分析響應(yīng)面設(shè)計(jì)點(diǎn)的選取方案,給出實(shí)體靶船極限承載能力的概率分布。
工程結(jié)構(gòu)的失效概率可以表示為[6]:
式中:G(X)是極限狀態(tài)方程,即響應(yīng)面方程,當(dāng)G(X)<0,結(jié)構(gòu)失效;Df是與G(X)相對(duì)應(yīng)的失效區(qū)域;f(X)=f(x1,x2,···,xn)是 基本隨機(jī)變量X的聯(lián)合密度函數(shù),當(dāng)X為一組相互獨(dú)立的隨機(jī)變量時(shí),有f(X)=
用蒙特卡羅法計(jì)算結(jié)構(gòu)可靠性時(shí),式(1)可表達(dá)為:
式中:N為抽樣模擬總數(shù);當(dāng)用相對(duì)誤差來(lái)表示蒙特卡羅法的抽樣誤差有:
式中: 是待定系數(shù);和 分別是基本隨機(jī)變量。
為使響應(yīng)面函數(shù)更好地逼近真實(shí)狀態(tài),設(shè)計(jì)點(diǎn)選取非常重要。一方面盡量要少,以降低計(jì)算工作量,另一方面還要盡量包含較多極限狀態(tài)信息,以保證響應(yīng)面的擬合精度。
響應(yīng)面-蒙特卡羅法[7]首先根據(jù)隨機(jī)變量的概率分布選取確定功能函數(shù)值的設(shè)計(jì)點(diǎn)。設(shè)計(jì)點(diǎn)一般以隨機(jī)變量的均值點(diǎn)為中心,根據(jù) 3 σ 原則在 μi±kσi范圍內(nèi)選取。進(jìn)而,根據(jù)設(shè)計(jì)點(diǎn)方案運(yùn)用有限元數(shù)值計(jì)算得到設(shè)計(jì)點(diǎn)對(duì)應(yīng)的功能函數(shù)值。根據(jù)有限元計(jì)算獲得的一系列功能函數(shù)值,采用最小二乘法等方法擬合響應(yīng)面,來(lái)替代真實(shí)的極限狀態(tài)方程。在此基礎(chǔ)上,運(yùn)用蒙特卡羅法對(duì)隨機(jī)變量進(jìn)行抽樣,采用響應(yīng)面進(jìn)行N次數(shù)值計(jì)算,得到N個(gè)功能函數(shù)值Zj(j=1,2,···,N),進(jìn)而確定功能函數(shù)的概率分布,或根據(jù)統(tǒng)計(jì)Zj<0的個(gè)數(shù),計(jì)算失效概率或可靠度指標(biāo)。因此,響應(yīng)面與蒙特卡羅法結(jié)合使結(jié)構(gòu)分析的工作量大大減少,提高效率。
基于非線性有限元法的船體極限強(qiáng)度計(jì)算精度主要取決于模型構(gòu)建的準(zhǔn)確性和非線性求解器的選擇與參數(shù)設(shè)置[8]。
為了避免邊界條件的影響,模型范圍一般選擇艙段模型。板和主要縱向構(gòu)件應(yīng)采用4節(jié)點(diǎn)單元,且盡量保證網(wǎng)格形狀接近正方形。為了準(zhǔn)確模擬極限狀態(tài)結(jié)構(gòu)的塑性大變形和失效模式演化,骨材腹板一般要求至少布置3~4個(gè)單元,并以此為基準(zhǔn)劃分整個(gè)艙段模型的網(wǎng)格往往可以取得較好的效果。同時(shí),極限強(qiáng)度計(jì)算需要考慮結(jié)構(gòu)的初始撓度[9],包括:
板架的整體變形
板的局部變形
骨材的側(cè)傾變形
式中:a和b分別為橫向和縱向構(gòu)件的間距;S為縱桁間距;hw為 骨材腹板高度;m為 半波數(shù),取滿足a/b≤的 最小正整數(shù);A0,B0和C0為初始撓度的幅值,分別取
由于極限強(qiáng)度計(jì)算需要考慮材料的非線性,在缺少船用鋼材料拉伸曲線的情況下,一般采用理想彈塑性或理想塑性強(qiáng)化模型以考慮材料的非線性效應(yīng)。
目前,極限強(qiáng)度分析的非線性求解方法主要有Newton-Raphson迭代法、Riks法和顯式動(dòng)態(tài)算法等。其中,Newton-Raphson法以載荷增量為控制量,根據(jù)節(jié)點(diǎn)力平衡判斷收斂,但結(jié)構(gòu)大變形引起的節(jié)點(diǎn)力不平衡容易造成不收斂。Riks法采用弧長(zhǎng)作為控制變量,改進(jìn)了Newton-Raphson法的迭代策略,計(jì)算效率和收斂性主要取決于弧長(zhǎng)相關(guān)參數(shù)的設(shè)置。顯式動(dòng)態(tài)算法采用中心差分法,用上一步和當(dāng)前步的結(jié)果計(jì)算下一步的計(jì)算結(jié)果,無(wú)需迭代運(yùn)算,收斂性好,計(jì)算效率較高,目前被廣泛地應(yīng)用于結(jié)構(gòu)靜態(tài)和準(zhǔn)靜態(tài)的非線性響應(yīng)分析。
以某型超役水面艦艇為例,該船已服役近30年,將改裝成為實(shí)體靶船??紤]材料屈服限和結(jié)構(gòu)腐蝕的隨機(jī)變化,采用響應(yīng)面-蒙特卡羅法結(jié)合非線性有限元法,開(kāi)展了極限強(qiáng)度可靠性分析。
選取船中艙段為分析對(duì)象,采用Abaqus通用有限元軟件建立有限元模型,網(wǎng)格尺寸整體上取50 mm×50 mm,同時(shí)縱骨腹板劃分3個(gè)單元,單元類型選用4節(jié)點(diǎn)薄殼單元S4R,如圖1所示。該艦材料屈服限為235 MPa,楊氏模量為206 GPa,泊松比0.3
圖1 艙段有限元模型Fig. 1 FE model of target cabin
由于是實(shí)體靶船,因此結(jié)構(gòu)腐蝕是重點(diǎn)考慮的一個(gè)隨機(jī)變量,同時(shí)由于材料屈服限對(duì)于艦船極限強(qiáng)度影響較大,因此主要考慮腐蝕量和材料屈服限的隨機(jī)變化。材料屈服限服從正態(tài)分布,變異系數(shù)取0.05,即σy~N(235,138.06)[4]。
國(guó)內(nèi)外學(xué)者通過(guò)對(duì)結(jié)構(gòu)腐蝕統(tǒng)計(jì)分析和試驗(yàn)研究,提出了一些腐蝕模型[10]。秦圣平等[11]在現(xiàn)有腐蝕模型比較分析的基礎(chǔ)上,提出了一種適用于船體結(jié)構(gòu)時(shí)變可靠性分析的非線性腐蝕模型,該模型可以較好模擬鋼結(jié)構(gòu)在腐蝕環(huán)境下的腐蝕損傷過(guò)程。表達(dá)式為:式中:d∞,β,η,Tst為待定參數(shù),均服從正態(tài)分布,其統(tǒng)計(jì)特征值見(jiàn)表1。
表1 腐蝕參數(shù)的統(tǒng)計(jì)特征值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
采用顯式動(dòng)態(tài)算法分析該艦中垂?fàn)顟B(tài)的極限強(qiáng)度,圖3為艙段的彎矩-轉(zhuǎn)角曲線,圖4為極限狀態(tài)結(jié)構(gòu)應(yīng)力云圖??梢园l(fā)現(xiàn),該艦的中垂極限彎矩為36.21 MN·m,極限狀態(tài)該艦強(qiáng)橫梁間的甲板板架發(fā)生了屈曲破壞。
根據(jù)材料屈服限和結(jié)構(gòu)腐蝕量的概率分布及其統(tǒng)計(jì)特征值,采用 3 σ原則確定49個(gè)響應(yīng)面設(shè)計(jì)點(diǎn),采用非線性有限元法計(jì)算獲得了不同腐蝕量和材料屈服限組合下的極限強(qiáng)度,如表2所示。
圖3 彎矩-轉(zhuǎn)角曲線Fig. 3 Bending moment vs. angle of rotation relationship
圖4 極限狀態(tài)艙段應(yīng)力云圖Fig. 4 Stress patterns of cabin at ultimate strength level
表2 不同材料屈服限和腐蝕量下的艙段極限強(qiáng)度(MN·m)Tab. 2 Ultimate strength of cabin under different yield strength and corrosion amount
從艙段極限強(qiáng)度的計(jì)算結(jié)果可以發(fā)現(xiàn):隨著材料屈服限的提高,極限強(qiáng)度逐步增加;隨著結(jié)構(gòu)腐蝕量的增加,極限強(qiáng)度不斷降低??傮w而言,材料屈服限對(duì)極限強(qiáng)度的影響較腐蝕量的影響要大。
在此基礎(chǔ)上,采用三次多項(xiàng)式表達(dá)式對(duì)響應(yīng)面進(jìn)行擬合,設(shè)計(jì)了10種不同設(shè)計(jì)點(diǎn)方案,分別包括49,25,21,21,17,17,17,13,13和13個(gè)設(shè)計(jì)點(diǎn),如圖5所示。表3為相關(guān)設(shè)計(jì)點(diǎn)方案的響應(yīng)面擬合結(jié)果,圖6為部分設(shè)計(jì)點(diǎn)方案擬合得到的響應(yīng)面。為了評(píng)價(jià)響應(yīng)面擬合效果,引入了局部擬合優(yōu)度和整體擬合優(yōu)度兩個(gè)評(píng)價(jià)指標(biāo)。其中,局部擬合優(yōu)度反映當(dāng)前設(shè)計(jì)點(diǎn)方案n個(gè)設(shè)計(jì)點(diǎn)與該方案下擬合得到的響應(yīng)面之間的擬合優(yōu)度;總體擬合優(yōu)度指全部49個(gè)設(shè)計(jì)點(diǎn)與n個(gè)設(shè)計(jì)點(diǎn)方案擬合得到的響應(yīng)面之間的擬合優(yōu)度,表達(dá)式為:
圖5 響應(yīng)面設(shè)計(jì)點(diǎn)方案Fig. 5 Design point schemes for response surface
表3 響應(yīng)面擬合結(jié)果Tab. 3 Fitting results of response surfaces
圖6 不同設(shè)計(jì)點(diǎn)方案的響應(yīng)面擬合結(jié)果Fig. 6 Response surfaces of different design point scheme
式中:zi為設(shè)計(jì)點(diǎn)的有限元計(jì)算結(jié)果,即原始數(shù)據(jù);為根據(jù)響應(yīng)面函數(shù)得到設(shè)計(jì)點(diǎn)的極限強(qiáng)度;表示原始數(shù)據(jù)的均值,即49個(gè)有限元計(jì)算結(jié)果的均值。
從響應(yīng)面的擬合結(jié)果可以發(fā)現(xiàn):前9個(gè)方案的局部擬合效果均大于0.99;同時(shí),前8個(gè)方案的整體擬合優(yōu)度均大于0.99,且方案1、方案4、方案6和方案8達(dá)到了0.999 8,對(duì)應(yīng)的設(shè)計(jì)點(diǎn)數(shù)量分別為49,21,17和13,說(shuō)明響應(yīng)面的整體擬合優(yōu)度與設(shè)計(jì)點(diǎn)的選取息息相關(guān),并不是設(shè)計(jì)點(diǎn)數(shù)量越多越好。方案10的設(shè)計(jì)點(diǎn)選取沒(méi)有考慮材料屈服限和腐蝕量的相關(guān)性,擬合得到的響應(yīng)面已失真,如圖6(d)所示,局部擬合優(yōu)度僅0.482。方案8雖然只選取了13個(gè)設(shè)計(jì)點(diǎn),但其整體擬合優(yōu)度與49個(gè)設(shè)計(jì)點(diǎn)一致,達(dá)到了0.999 8,局部擬合優(yōu)度甚至達(dá)到了1.0,擬合優(yōu)度比方案2(25個(gè)設(shè)計(jì)點(diǎn))和方案3(21個(gè)設(shè)計(jì)點(diǎn))都要好。
因此,響應(yīng)面設(shè)計(jì)點(diǎn)的合理選取,不僅可以有效降低計(jì)算工作量,同時(shí)可以保證響應(yīng)面擬合精度。
獲得極限強(qiáng)度的響應(yīng)面后,結(jié)合相關(guān)參數(shù)的概率分布,采用蒙特卡羅法隨機(jī)抽樣,進(jìn)一步分析實(shí)體靶船極限強(qiáng)度的概率分布。
首先,開(kāi)展了蒙特卡羅法的收斂性分析,抽樣次數(shù)分別取1 000,10 000,100 000和1 000 000次,樣本的統(tǒng)計(jì)特征值見(jiàn)表5??梢园l(fā)現(xiàn),當(dāng)樣本數(shù)量達(dá)到10 000次,極限強(qiáng)度的均值與標(biāo)準(zhǔn)差已趨于穩(wěn)定;當(dāng)樣本數(shù)量增加到100 000次或更高,極限強(qiáng)度的統(tǒng)計(jì)特征值變化不大。在正態(tài)分布的基礎(chǔ)上,進(jìn)一步采用對(duì)數(shù)正態(tài)分布和Weibull分布擬合極限強(qiáng)度的概率分布,擬合結(jié)果如圖7所示。其中,正態(tài)分布和對(duì)數(shù)正態(tài)分布的分布函數(shù) 分 別 為和 l nWeibull分布的形狀參數(shù)為37.16,尺度參數(shù)為21.47,分布函數(shù)為
表4 樣本的均值與標(biāo)準(zhǔn)差Tab. 4 Mean and standard deviation of samples
圖7 極限強(qiáng)度的概率分布Fig. 7 Probability distribution of ultimate strength
可以發(fā)現(xiàn):極限強(qiáng)度的隨機(jī)分布與Weibull分布偏差較大。因此,實(shí)體靶船極限強(qiáng)度的概率分布更趨近于正態(tài)分布或?qū)?shù)正態(tài)分布。
本文采用蒙特卡羅法與響應(yīng)面法相結(jié)合的可靠性分析方法,開(kāi)展考慮材料屈服限和腐蝕影響的艙段極限強(qiáng)度研究,研究表明:在保證擬合優(yōu)度的前提下,響應(yīng)面設(shè)計(jì)點(diǎn)的合理選取,可以有效降低計(jì)算工作量;考慮材料屈服限和腐蝕聯(lián)合影響的實(shí)體靶船極限強(qiáng)度概率分布更趨近于正態(tài)分布或?qū)?shù)正態(tài)分布。本文研究建立的基于響應(yīng)面-蒙特卡羅法的可靠性分析方法,對(duì)于實(shí)體靶船極限強(qiáng)度評(píng)估具有指導(dǎo)意義。