基于GPC的環(huán)肋耐壓圓柱殼結(jié)構(gòu)失穩(wěn)概率分析

張毅博， 孫志禮， 趙中強(qiáng)， 趙經(jīng)武

(1. 東北大學(xué) 機(jī)械工程與自動(dòng)化學(xué)院, 遼寧 沈陽 110819； 2. 中車唐山機(jī)車車輛有限公司， 河北 唐山 063000；3. 中國人民解放軍第93107部隊(duì), 遼寧 沈陽 110141)

隨著對(duì)隱蔽性要求的不斷提高，潛艇下潛深度逐漸增加，環(huán)肋耐壓圓柱殼結(jié)構(gòu)作為潛艇的基本結(jié)構(gòu)，其穩(wěn)定性問題隨靜水壓力的增大而愈發(fā)突出[1].因此，開展深潛環(huán)肋耐壓圓柱殼結(jié)構(gòu)失穩(wěn)概率分析具有重要意義和工程價(jià)值.

可靠性分析的是產(chǎn)品在特定條件下完成特定功能的能力.因此，本文采用可靠性方法評(píng)估環(huán)肋耐壓圓柱殼結(jié)構(gòu)的失穩(wěn)概率.在可靠性領(lǐng)域，結(jié)構(gòu)狀態(tài)通常被分為兩類：安全和失效.近年來，為提高可靠性分析效率，諸多學(xué)者提出了各種基于代理模型(如響應(yīng)面[2]、Kriging模型[3]、支持向量機(jī)SVM[4]、人工神經(jīng)網(wǎng)絡(luò)ANN等)的分析方法.然而，現(xiàn)有可靠性分析方法主要針對(duì)具有連續(xù)響應(yīng)的結(jié)構(gòu)，而對(duì)于具有失穩(wěn)(即屈曲)破壞的拱結(jié)構(gòu)、桁架結(jié)構(gòu)、薄殼結(jié)構(gòu)等的分析研究甚少.盡管可以采用SVM和ANN處理結(jié)構(gòu)狀態(tài)二分類問題，然而SVM和ANN均需要大量已標(biāo)注樣本才能得到較為精確的分類結(jié)果.此外，SVM仍存在核函數(shù)、核參數(shù)以及損失函數(shù)較難選取等問題.

作為高斯過程的一個(gè)分支，高斯過程分類(GPC)不僅具有嚴(yán)格的統(tǒng)計(jì)理論基礎(chǔ)，還能自適應(yīng)獲取超參數(shù)[5].然而，GPC執(zhí)行推理的時(shí)間隨訓(xùn)練樣本數(shù)的增加呈立方增長.為此，學(xué)者們提出了各種自適應(yīng)試驗(yàn)設(shè)計(jì)(DoE)策略，如Kapoor等[6]將最易分類錯(cuò)誤的點(diǎn)(MEMP)作為新訓(xùn)練樣本；Peng等[7]將最可能失效點(diǎn)(MPP)添加到DoE中.然而，現(xiàn)有DoE策略選擇出來的訓(xùn)練樣本容易聚集，即引起不必要的功能函數(shù)評(píng)估.此外，對(duì)于具有小失效概率結(jié)構(gòu)的可靠性評(píng)估，盡管Yang等[8]提出了一種雙循環(huán)且不依賴于MPP的基于Krging和重要抽樣(IS)的可靠性分析方法，然而目前尚未發(fā)現(xiàn)相關(guān)基于GPC和IS的可靠性分析方法.

為此，本文提出了一種單循環(huán)的基于GPC和IS的自適應(yīng)分析方法，既避免了雙循環(huán)中優(yōu)化算法易陷入局部最優(yōu)解的缺陷，又保證了在分類面附近區(qū)域均勻取樣.通過某一分段函數(shù)驗(yàn)證了所提方法的高效性與準(zhǔn)確性，并將其用于評(píng)估某深潛環(huán)肋耐壓圓柱殼結(jié)構(gòu)的失穩(wěn)概率.

1 基本理論

1.1 高斯過程分類

高斯過程分類通常使用類標(biāo)簽+1或-1處理二元分類問題，+1代表結(jié)構(gòu)安全，-1代表結(jié)構(gòu)失效.給定N個(gè)訓(xùn)練樣本X(X=[x1,…,xN])及其類標(biāo)簽Y(Y=[y1,…,yN])，GPC通過建立X和Y間的映射關(guān)系(即潛函數(shù)f(x))來進(jìn)行分類.

為對(duì)未知樣本x*分類，需由式(1)計(jì)算x*相應(yīng)潛變量f*的后驗(yàn)分布.

(1)

式中：f為訓(xùn)練樣本X對(duì)應(yīng)的潛向量(f=[f(x1),f(x2),…,f(xN)]T)；p(f*|X,x*,f)為f*的條件先驗(yàn)分布；p(f|X,Y)為f的后驗(yàn)分布.

為求上述積分，GPC假定條件分布p(f|X)服從高斯分布，即

p(f|X)=N(0,K)

(2)

式中，K為協(xié)方差矩陣(Kij=k(xi,xj))；k(·)為協(xié)方差函數(shù)，本文選用應(yīng)用最廣泛的平方指數(shù)協(xié)方差函數(shù)(也稱高斯核函數(shù))來計(jì)算K.

(3)

式中：σf用來控制局部相關(guān)程度；指數(shù)項(xiàng)用來表征xi和xj間的距離相關(guān)性，即如果xi和xj的距離相對(duì)于距離尺度l很小，則它們高相關(guān)，否則低相關(guān).

由高斯過程可知，f*和f的聯(lián)合分布也服從高斯分布，即

(4)

式中：Kx*=[k(x*,x1),…,k(x*,xN)]T；kx*=k(x*,x*).

因此，條件先驗(yàn)分布p(f*|X,x*,f)為

(5)

由貝葉斯規(guī)則可得f的后驗(yàn)分布p(f|X,Y)為

p(f|X,Y)=p(Y|f)p(f|X)/p(Y|X)

(6)

將式(5)和式(6)代入式(1)中，則樣本x*的類標(biāo)簽y*為+1的概率為

(7)

然而，式(1)和式(7)沒有解析解，需采用數(shù)值近似法求解上述積分，本文選用Laplace近似法.Laplace采用高斯分布近似p(f*|X,Y,x*)，即

(8)

式中：

因此，式(7)中的分類概率近似為

(9)

式中，σ(·)為響應(yīng)函數(shù)，本文選擇標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)作為響應(yīng)函數(shù).

顯然，分類概率0.5為結(jié)構(gòu)安全與否的界限.

1.2 基于重要抽樣的可靠性分析理論

在可靠性分析中，失效概率定義為

(10)

式中：G(x)為結(jié)構(gòu)的功能函數(shù)；fX(x)為x的聯(lián)合概率密度函數(shù)；x=[x1,…,xM]，M為影響結(jié)構(gòu)性能的不確定因素(即隨機(jī)變量)的個(gè)數(shù).

盡管可以采用各種基于Monte Carlo采樣的方法計(jì)算失效概率，但當(dāng)失效概率很小時(shí)(如Pf<10-4)，需大量的Monte Carlo樣本(>107)才能得到較為精確的結(jié)果.為此，需采用重要抽樣等方差縮減技術(shù)來減少可靠性評(píng)估中隨機(jī)抽樣樣本的數(shù)量.根據(jù)重要抽樣，式(10)可改寫為

(11)

(12)

(13)

2 基于GPC和IS的自適應(yīng)分析方法

2.1 所提自適應(yīng)DoE策略

現(xiàn)有基于GPC的自適應(yīng)DoE策略主要包含兩種：基于最可能失效點(diǎn)(GPC+MPP)和基于最易分類錯(cuò)誤點(diǎn)(GPC+MEMP).

在可靠性分析中，MPP是對(duì)失效概率影響最大的點(diǎn)，且MPP附近區(qū)域高斯過程分類器的擬合精度直接影響著失效概率估計(jì)的準(zhǔn)確性.因此，文獻(xiàn)[7]將MPP作為新訓(xùn)練樣本，其可由式(14)近似獲取.

(14)

由式(9)可知，分類概率越接近0.5，其類標(biāo)簽越容易分類錯(cuò)誤.因此，文獻(xiàn)[6]將滿足式(15)的點(diǎn)，即最易分類錯(cuò)誤的點(diǎn)MEMP作為新訓(xùn)練樣本.

(15)

然而，上述選點(diǎn)策略存在一些缺陷：1)他們的候補(bǔ)樣本集是基于Monte Carlo采樣得到的，難以應(yīng)用于小失效概率評(píng)估；2)選擇出來的訓(xùn)練樣本容易聚集，從而增加不必要的功能函數(shù)評(píng)估次數(shù).為此，本文提出采用馬爾科夫鏈蒙特卡洛法(MCMC)和歐式距離來避免上述缺陷.

由上述可知，高斯過程分類器在分類概率為0.5處的擬合精度對(duì)失效概率估計(jì)的準(zhǔn)確性起決定性影響.因此，應(yīng)將分類概率為0.5的樣本作為候補(bǔ)樣本.為此，本文采用MCMC法產(chǎn)生NCS個(gè)滿足式(16)的樣本并將其作為候補(bǔ)樣本集.

|p(y=+1|X,Y,x)-0.5|≤[ε]

(16)

式中，[ε]取為0.01.

為避免訓(xùn)練樣本聚集，本文提出采用歐式距離來保證取樣均勻性，其中歐式距離為

d=|xc-xt|

(17)

式中：xc為候補(bǔ)樣本；xt為現(xiàn)有訓(xùn)練樣本.

因此，本文將滿足式(18)的點(diǎn)作為新訓(xùn)練樣本，

(18)

(19)

式中，dmin為xc距DoE中所有訓(xùn)練樣本的最短距離.

2.2 基于核密度估計(jì)的重要抽樣

重要抽樣的核心是重要抽樣密度函數(shù)h(x)的構(gòu)造，理論上最優(yōu)h(x)為

h(x)=If(x)fX(x)/Pf

(20)

(21)

式中：ω為帶寬；K(·)為核概率密度函數(shù)，本文選取正態(tài)分布的概率密度函數(shù)作為K(·)，其表達(dá)式為

(22)

帶寬ω可由最小化平均積分平方誤差獲取，

(23)

2.3 停止準(zhǔn)則

現(xiàn)有停止準(zhǔn)則主要有兩種，然而它們由于其固有缺陷而難以被工程實(shí)際所采納.

1) 功能函數(shù)調(diào)用次數(shù)達(dá)到某一規(guī)定的閾值[9].顯然，若規(guī)定的閾值過小，則易錯(cuò)估失效概率；反之，引起不必要的功能函數(shù)評(píng)估.

2) 前后兩次迭代失效概率估計(jì)的相對(duì)誤差小于某一閾值[7]，即

(24)

顯然，若選擇的新訓(xùn)練樣本對(duì)高斯過程分類器的精度改善甚微，則提前終止迭代且得到錯(cuò)誤的失效概率估計(jì).

為此，本文基于失效概率估計(jì)的穩(wěn)定性，提出了一種新的、更加精確的停止準(zhǔn)則，即

(25)

(26)

2.4 所提自適應(yīng)分析方法

本文所提自適應(yīng)分析方法的基本步驟如下：

步驟1 通過Nataf變換將非獨(dú)立隨機(jī)變量轉(zhuǎn)換為獨(dú)立標(biāo)準(zhǔn)正態(tài)變量，采用拉丁超立方抽樣隨機(jī)生成N個(gè)初始訓(xùn)練樣本SDoE=[x1,x2,…,xN]并調(diào)用真實(shí)功能函數(shù)獲取類標(biāo)簽YDoE=[y1,y2,…,yN].

步驟3 基于當(dāng)前高斯過程分類器，應(yīng)用MCMC法生成NCS個(gè)滿足式(16)的候補(bǔ)樣本(本文令NCS=3 000).根據(jù)式(18)選擇新訓(xùn)練樣本，調(diào)用功能函數(shù)獲取其對(duì)應(yīng)的真實(shí)類標(biāo)簽，并將新訓(xùn)練樣本添加到DoE中，返回步驟2.

3 算例驗(yàn)證

以一分段函數(shù)為例，驗(yàn)證所提自適應(yīng)分析方法的準(zhǔn)確性及高效性，其表達(dá)式為[10-11]

(27)

式中，x1，x2相互獨(dú)立且均服從標(biāo)準(zhǔn)正態(tài)分布.

取12個(gè)拉丁超立方樣本作為初始訓(xùn)練樣本，且每種方法各運(yùn)行10次以消除由不同初始訓(xùn)練樣本及候補(bǔ)樣本的隨機(jī)性所產(chǎn)生的影響，其中GPC+MPP和GPC+MEMP以基于核密度估計(jì)的重要抽樣(KDE-IS)生成的重要抽樣樣本作為候補(bǔ)樣本.不同方法運(yùn)行10次的平均結(jié)果如表1所示.此外，圖1a，1b還分別顯示了由不同方法得到的失效概率估計(jì)和最大相對(duì)穩(wěn)定性ε隨迭代次數(shù)Nit的變化趨勢.顯然，由表1和圖1b可知，所提方法在滿足失效概率分析精度要求的同時(shí)，還需要更少的功能函數(shù)評(píng)估次數(shù).圖2比較了不同方法在調(diào)用真實(shí)功能函數(shù)28次時(shí)DoE中訓(xùn)練樣本的分布及預(yù)測的高斯過程分類邊界.顯然，所提自適應(yīng)分析方法選擇的訓(xùn)練樣本分布更均勻且能更好地收斂到真實(shí)分類邊界，即需要調(diào)用真實(shí)功能函數(shù)次數(shù)更少.

表1 不同方法運(yùn)行10次的平均結(jié)果Table 1 Average results of 10 runs of different methods

4 環(huán)肋耐壓圓柱殼失穩(wěn)概率分析

某深潛環(huán)肋耐壓圓柱殼結(jié)構(gòu)采用矩形肋骨，其截面如圖3所示，其中，r為耐壓殼體內(nèi)徑，t為耐壓殼體厚度，l為肋骨間距，a為肋骨高度，b為肋骨寬度.各參數(shù)及環(huán)肋耐壓圓柱殼所受外部壓力p均服正態(tài)分布，相應(yīng)分布參數(shù)如表2所示[12].

表2 各隨機(jī)變量的分布特征Table 2 Distribution of random variables

為評(píng)估該型環(huán)肋耐壓圓柱殼結(jié)構(gòu)的失穩(wěn)概率，需采用ANSYS中的更符合工程實(shí)際和更具應(yīng)用價(jià)值的非線性屈曲分析來分析其穩(wěn)定性.由于環(huán)肋耐壓圓柱殼結(jié)構(gòu)的對(duì)稱性，僅建立其1/4有限元模型，如圖4所示.模型材料為某高強(qiáng)度鋼，其彈性模量為2.0×105MPa，泊松比為0.3，屈服強(qiáng)度為785 MPa.圓柱殼采用殼單元shell181建模，矩形肋骨采用beam188單元建模.模型一端固定約束，另一端允許軸向位移，1/4截面處施加對(duì)稱約束，殼結(jié)構(gòu)外表面上施加壓力載荷.

由ANSYS中非線性屈曲分析得到的該型環(huán)肋耐壓圓柱殼結(jié)構(gòu)在穩(wěn)定和失穩(wěn)時(shí)的節(jié)點(diǎn)位移云圖分別如圖5、圖6所示.其中，圖6中該型環(huán)肋耐壓圓柱殼結(jié)構(gòu)的失穩(wěn)類型為整體失穩(wěn).

通過建立該型環(huán)肋耐壓圓柱殼結(jié)構(gòu)參數(shù)化有限元模型及MATLAB和ANSYS的聯(lián)合仿真，基于所提分析方法，評(píng)估其失穩(wěn)概率.為此，在標(biāo)準(zhǔn)正態(tài)空間內(nèi)通過拉丁超立方抽樣隨機(jī)選取30個(gè)初始訓(xùn)練樣本并調(diào)用ANSYS仿真得到其相應(yīng)最大節(jié)點(diǎn)位移；構(gòu)造初始高斯過程分類器，根據(jù)2.1節(jié)所提自適應(yīng)DoE策略選擇新訓(xùn)練樣本，調(diào)用ANSYS獲取最大節(jié)點(diǎn)位移并加入到DoE中.經(jīng)過238次循環(huán)仿真后，滿足收斂條件，即迭代停止，所得該型深潛環(huán)肋耐壓圓柱殼結(jié)構(gòu)的失穩(wěn)概率約為8.242×10-5.

5 結(jié) 論

1) 通過某一分段函數(shù)表明，所提方法在滿足失效概率分析精度要求的同時(shí)還極大地提高了分析效率，避免了不必要的功能函數(shù)評(píng)估.

2) 將所提自適應(yīng)分析方法應(yīng)用于評(píng)估某型深潛環(huán)肋耐壓圓柱殼結(jié)構(gòu)的失穩(wěn)概率，結(jié)果表明：所提方法能夠很好地適用于工程實(shí)際，這為預(yù)防潛艇失穩(wěn)、縮短新型潛艇的研制周期及解決潛艇下潛的穩(wěn)定性評(píng)估提供了一種切實(shí)可行的分析方法.