基于重要抽樣法的耐壓球殼可靠性計算方法研究

馮士超，萬正權，李艷青

（1.中國船舶科學研究中心，江蘇 無錫 214082；2.深海技術科學太湖實驗室，江蘇 無錫 214082；3.深海載人裝備國家重點實驗室，江蘇 無錫 214082）

深海載人潛水器[1-3]是進入、探測、開發和保護深海的重要技術手段和裝備。載人艙耐壓球殼是保證下潛人員及艙內設備免受外界海水巨大壓力的耐壓結構，是載人潛水器的核心部件[4-5]。目前大深度載人潛水器多選用鈦合金耐壓球殼形式，其具有結構簡單、容重比高、承載能力強、有相同的周向和徑向應力值等優點，其質量約占整個潛水器系統的1/5～1/4[6]。載人艙耐壓球殼的可靠性水平對整個載人潛水器能夠完成設定的任務有決定性影響，同時準確獲知耐壓球殼的可靠性水平將為開展結構優化設計、減輕結構質量提供重要指導，開展耐壓球殼的安全可靠性研究具有重要的意義。

國內學者開展的耐壓球殼可靠性評估結果均表明，在工作載荷下，耐壓球殼失效是極小概率問題[7-8]，且球殼承載能力的計算屬于單失效模式問題。對于失效方程未知的隱式問題，國內外學者在代理模型方法的基礎上，對極小失效概率的可靠性求解問題也開展了相關研究。Echard 等[9]提出了一種結合重要抽樣和自適應Kriging 方法的能用于處理小失效概率問題的方法，該方法首先需要使用改進一次二階矩法求得設計點，然后以設計點為中心，建立自適應Kriging 模型，但在計算設計點時，需要使用失效函數對各個變量的偏微分信息，該要求對于隱式問題往往較為困難。Huang 等[10]提出了一種結合子集模擬法和自適應kriging 方法的能處理小失效概率問題的方法，該方法主要是在生成的Kriging 模型基礎上，采用子集模擬法計算失效概率。Lv 等[11]提出了一種結合方向抽樣和自適應Kriging 方法計算小概率失效問題的方法，該方法是在生成的Kriging 模型基礎上，采用方向抽樣法計算失效概率。Zhang 等[12]提出了一種基于Kriging 的自適應重要抽樣方法，該方法分為2 步，首先通過較稀疏的樣本點識別重要抽樣中心，然后在重要抽樣中心附近建立最終的Kriging 模型。文獻[10-12]在計算極小概率失效問題時，均需要設計較大的樣本空間用于搜索生成Kriging 模型，但部分Kriging 模型位于失效概率較低的空間，實質上低失效概率空間的Kriging 模型對最終計算結果影響不大，這導致此類方法的計算效率不高，且相關文獻中一般均不考慮不同失效條件下的可靠性問題。因此，本文開展了單失效模式下不同失效條件的可靠性問題的計算方法研究。

本文結合自適應Kriging 和重要抽樣法，提出了一種用于求解失效條件為非隨機變量的極小失效概率問題的重要抽樣法。首先，使用自適應Kriging 方法構建較高失效概率下的代理模型，在此基礎上獲得重要方向，在重要方向上計算得到較低失效概率下的設計點，然后以此設計點為中心，使用自適應Kriging方法構建代理模型，并開展可靠性計算。本文首先通過2 個算例驗證了該方法的計算精度與效率，并應用到某耐壓球殼結構可靠性計算問題，在獲得計算載荷可靠性指標的基礎上，計算得到了工作載荷下結構的可靠性指標。

1 自適應Kriging 方法

1.1 Kriging 模型

在可靠性分析中，如果直接采取高精度的仿真模型進行可靠性計算，將產生巨大的計算成本，因此一種以較小計算代價獲得研究對象較為準確狀態函數的方法——代理模型技術應運而生。代理模型技術[13]就是在保證一定擬合精度的前提下，以某種合適的代理模型方法進行數據擬合，近似獲得一個數學函數模型去替代高耗時分析模型進行需要較大樣本量的分析工作。Kriging 模型是一種較為常用的構建代理模型的方法。

Kriging 模型的思想由南非采礦工程師Krige 提出[14]。1989 年，Sacks 等[15]將Kriging 理論進一步推廣并給出了一種較為實用的 Kriging 算法，使得Kriging 模型在地質、氣象、航空航天、船舶等領域得到應用和發展。Kriging 模型一般可表示為一個含參數的線性回歸模型和一個隨機分布之和，即：

Kriging 模型不僅可以對樣本點的函數值進行預測，還能估計預測值的方差，該方差可以用于指導加入新的樣本點，以提高建立的代理模型的精度。

1.2 自適應Kriging

為更高效建立Kriging 模型用于開展可靠性分析，自適應Kriging 方法[16-18]被提出并得到了廣泛的應用。自適應Kriging 的主要思想是使用學習函數篩選對建立Kriging 模型效果最佳的樣本點加入訓練樣本集，并在此基礎上不斷更新代理模型，其主要流程如圖1 所示。

圖1 自適應Kriging 模型建立流程Fig.1 Process of constructing an adaptive Kriging model

ERF（Expected Risk Function）是一種自適應Kriging 常用的學習函數，其基本思想是基于預測值符號（正負）不確定性的思想，選擇候選樣本中符號預測風險最大的樣本用于更新實驗設計[18]，學習函數ERF 一般被定義為：

其中，sign (?)是符號函數，ERF 函數值最大樣本點的符號被誤判的風險最大，候選樣本集中ERF值最大的樣本被識別出來用于更新實驗設計，并重新生成Kriging 模型。樣本的ERF 值越小，樣本符號被誤判的概率就越低，文獻[18]給出收斂條件建議值為

2 重要抽樣法

重要抽樣法通過采用重要抽樣密度函數代替原來的抽樣密度函數，使得樣本落入失效域的概率增加，以此獲得高的抽樣效率和快的收斂速度[19-22]。其基本原理和公式為，通過引入重要抽樣密度函數hx(x)，將求解失效概率的積分式轉換為以hx(x)為密度函數的數學期望形式，即：

在開展重要抽樣計算過程中，由于設計點是失效域內對失效概率貢獻最大的點，重要抽樣的中心一般選擇為設計點。這樣按照重要抽樣密度函數抽取的樣本點就有比較大的概率落在失效域內，從而提高了重要抽樣法的計算效率，使其更快地收斂于真值。

本文結合自適應Kriging 和重要抽樣法提出了一種用于求解失效條件為非隨機變量的極小失效概率問題的重要抽樣法，該方法僅適用于處理單失效模式的問題，且一般僅適用于失效條件為非隨機變量的失效問題的求解，可用于計算較低載荷下耐壓球殼的可靠性問題。該重要抽樣法首先使用自適應Kriging 方法構建較高失效概率下的模型g1Kriging(x)，在此基礎上，使用改進一次二階矩法獲得重要抽樣方向，在重要抽樣方向上仍然使用自適應Kriging 方法計算得到較低失效概率下的設計點MPP2，并以此設計點為中心，使用自適應Kriging 方法構建較低失效概率下的模型g2Kriging(x)，然后驗證g2Kriging(x)對應的設計點MPP3距離g2Kriging(x)樣本空間邊界的距離。若MPP3所有變量距離樣本空間的邊界大于等于隨機變量的3倍標準差，說明建立的Kriging 模型g2Kriging(x)能較好反映設計點附近的失效特征，可在g2Kriging(x)的基礎上開展重要抽樣計算；若小于3 倍的標準差，則重新以MPP3為中心建立Krging 模型，直至滿足判斷標準。計算流程如圖2 所示。其中重要抽樣方向為標準正態空間中從原點指向設計點的方向，該方向的物理意義是功能（失效）函數最速下降方向[23]。

圖2 重要抽樣法流程Fig.2 Process of importance sampling

3 數值算例

3.1 算例1

該算例來源于文獻[24]，其本身是一個非線性程度較高的失效問題，包含2 個標準正態分布隨機變量x1和x2，以其立方和不超過a建立約束，且x1和x2相互獨立，x1和x2的均值為0.5，標準差為0.2。a取1 時，失效概率較高；a取4 時，失效概率較低，該問題的功能函數為：

算例1 失效邊的界構建過程如圖3 所示。在建立的高失效概率Kriging 模型的基礎上，得到設計點MPP1，在均值指向MPP1的方向，通過自適應Kriging 方法計算得到低失效概率下的設計點MPP2，然后以MPP2為中心，建立低失效概率下的Kriging模型經驗證后，為失效函數開展重要抽樣，該算例的可靠性計算結果見表1。本文使用的重要抽樣方法與蒙特卡洛法的計算結果接近，相差2.5%，但調用功能函數次數為39 次，其中建立高失效概率下Kriging 模型為18 次，計算低失效概率下設計點MPP2為6 次，計算低失效概率下kriging模型為15 次，調用功能函數次數遠小于蒙特卡洛法。對于失效概率更低的情況，由于樣本數一般需滿足N=(102～104)/Pf的要求，直接采用蒙特卡洛法進行計算將是工程實際中無法接受的，而本文提出的重要抽樣法可以方便地計算出失效概率。

表1 算例1 計算結果Tab.1 Computed results of example 1

圖3 算例1Kriging 模型構建過程Fig.3 Kriging model construction process of example 1

3.2 算例2

該算例來源于文獻[23]，是一個包含3 個獨立隨機變量懸臂梁結構，其受到水平和豎直方向的載荷X和Y的作用（如圖4 所示），以其自由端位移不超過D0為約束建立功能函數：

圖4 懸臂梁結構Fig.4 Cantilever beam structure

式中：E、ω和t分別為材料的彈性模量、梁的寬度與厚度，ω=2.488 4m,t=3.888 4m、L=100m，E、X和Y為分布參數，見表2。0D取2.2 mm 時，失效概率較高；0D取2.8mm 時，失效概率較低。

表2 算例2 模型參數Tab.2 Model parameters of example 2

該算例的可靠性計算結果見表3。本文使用的重要抽樣方法與蒙特卡洛法的計算結果接近，相差1.4%，但調用功能函數次數為65 次。其中，建立高失效概率下Kriging 模型為31 次，計算低失效概率下設計點MPP2為3 次，計算低失效概率下Kriging 模型為31 次，調用功能函數次數遠小于蒙特卡洛法。對于失效概率更低的情況，由于樣本數一般需滿足N=(102～ 104)/Pf的要求，直接采用蒙特卡洛法進行計算將是工程實際中無法接受的，而本文提出的重要抽樣法可以方便地計算出失效概率。

表3 算例2 計算結果Tab.3 Computed results of example 2

4 耐壓球殼可靠性分析

本文研究的耐壓球殼的工作載荷為115 MPa，材料為鈦合金，安全系數取為1.5，計算載荷為172.5 MPa。主要計算工作載荷和設計載荷下耐壓球殼的可靠性指標，該結果可有助于評估結構的設計方案的安全可靠性以及為耐壓球殼安全使用和結構優化提供支撐。耐壓球殼的變量參數見表4。

表4 耐壓球殼參數Tab.4 Parameters of pressure resistant spherical shell

首先使用自適應 Kriging 方法計算得到的172.5 MPa 下耐壓球殼失效的Kriging 模型，樣本空間范圍為均值加減 5 倍的標準差，候選樣本量為2×105個。在初始樣本的基礎上，采用ERF 學習函數的計算結果不斷更新訓練樣本集。每個樣本點均采用ABAQUS 有限元軟件計算，參數化建模通過宏文件的方式完成，有限元模型如圖5 所示。在球殼承載能力計算過程中，首先使用線性有限元獲得結構的一階屈曲模態，將一階屈曲模態乘以球殼缺陷幅值作為球殼的初始缺陷，然后使用非線性有限元計算結構的承載能力[25]。文獻[18]給出的收斂條件建議值是基于數值算例給出的，數值算例在候選樣本點的響應值為該點響應值的精確解，而本文的響應值是通過有限元軟件給出的，其本質上屬于近似解，且計算過程需按照球殼尺寸重新生成有限元模型，網格等信息也隨之改變。在這些因素的共同作用下，仍選擇文獻[18]給出的建議值將不再合適。經計算驗證，在小于等于1×10–3后繼續加入樣本點，不會有明顯降低，且結構的失效概率也未出現明顯變化，因此收斂條件可取

圖5 耐壓球殼有限元模型Fig.5 Finite element model of pressure resistant spherical shell

在172.5 MPa 的載荷下，選擇10 個樣本點作為初始樣本點，在學習函數的指導下，通過加入14 個樣本點達到收斂指標，此時耐壓球殼的失效概率為4.501×10–2，對應設計點MPP1。然后在重要方向上使用自適應Kriging 方法計算在115 MPa 下的設計點，選擇2 個樣本點作為初始樣本點，在學習函數的指導下，通過加入2 個樣本點達到收斂指標，得到設計點MPP2。以MPP2為中心，仍然使用自適應Kriging 方法建立115 MPa 下耐壓球殼的Kriging 模型，選擇10個樣本點作為初始樣本點，在學習函數的指導下，通過加入11 個樣本點達到收斂指標，115 MPa 下建立的Kriging 模型的對應的設計點為MPP3。3 個設計點的數值見表5。計算可得MPP3設計點各個變量距離115 MPa 的Kriging 模型對應的樣本空間邊界均大于各個變量3 倍的標準差，說明最終建立的Kriging 模型能較好反映設計點附近的失效特征。因此，可以在此基礎上開展可靠性計算，計算可得115 MPa 下耐壓球殼的失效概率為4.094×10–96，失效概率變異系數為0.03。從115 MPa 下耐壓球殼的可靠性計算結果看，115 MPa 的工作載荷下，耐壓球殼的失效概率極低。

表5 設計點數值Tab.5 Value of design points

5 結語

本文提出了一種用于求解失效條件為非隨機變量的極小失效概率失效問題的重要抽樣方法，并將其應用于某耐壓球殼可靠性計算，可得到如下結論：

1）本文提出的重要抽樣方法可用于無失效方程下求解極低失效概率的問題，通過2 個算例的計算結果表明，該方法在保證計算精度的前提下，擁有較高的計算效率，但該方法僅能用于處理單失效模式的問題，且一般僅適用于失效條件為非隨機變量的失效問題的求解。

2）本文提出的重要抽樣法可適用于耐壓球殼的可靠性計算問題，115 MPa 工作載荷下，耐壓球殼的失效概率為4.094×10–96，可以認為該載荷下耐壓球殼失效的概率極低。

本文研究可為耐壓球殼的安全使用提供支撐。