涂宏茂，孫志禮，姬廣振，劉 勤

(1．東北大學機械工程與自動化學院，遼寧 沈陽 110819；2．中國兵器科學研究院，北京 100089)

引信MEMS微結構可靠性仿真及實現(xiàn)方法

涂宏茂1,2，孫志禮1，姬廣振2，劉 勤1,2

(1．東北大學機械工程與自動化學院，遼寧沈陽110819；2．中國兵器科學研究院，北京100089)

針對目前引信MEMS微結構設計中缺乏高效的可靠性量化仿真方法的問題，提出了基于結構性能仿真的可靠性仿真分析及其實現(xiàn)方法。該方法以ANSYS軟件為例給出CAE程序的集成方法和實現(xiàn)步驟，解決可靠性分析程序對有限元軟件工具的自動調用問題；采用試驗設計和近似建模相結合的方法，有效降低可靠性分析程序對有限元仿真的調用次數(shù)，在此基礎上，給出結構可靠度計算和靈敏度分析的基本原理；提出可靠性仿真實現(xiàn)的總體程序框架，給出各主要模塊的功能及其主要函數(shù)的偽代碼。實際驗證結果表明，該仿真方法具有可行性和實用性。

引信可靠性；MEMS微結構；可靠性仿真；有限元方法；可靠度

0 引言

目前，有限元等數(shù)值仿真方法越來越多地應用到引信MEMS的設計分析中。尤其是對于MEMS微結構(如引信用的MEMS懸臂梁、彈簧等，結構特征尺寸在1 μm至1 mm范圍內)而言，它們常具有相對的運動或相互的表面作用，利用有限元等仿真方法可以在設計階段很好地揭示其運動規(guī)律和力學特性[1-5]，并發(fā)現(xiàn)潛在的失效模式和失效原因，為設計改進與優(yōu)化提供了有益的參考。如果能在此基礎上，進一步開展量化的可靠性分析，一方面給出結構設計方案的可靠性水平，另一方面明確影響結構可靠性的主要因素，對于實現(xiàn)設計階段的可靠性預測和改進，具有重要的工程意義。從國外的一些研究文獻[6-8]來看，解決這一工程需求的基本思路是：將結構可靠性理論與有限元仿真方法相結合，進行MEMS微結構的可靠性仿真分析與優(yōu)化設計。但在實踐過程中，仍然存在兩個方面的問題：一是如何在可靠性分析程序中方便地集成CAE仿真程序或工具(如ANSYS、COMSOL Multiphysics等，后續(xù)簡稱為CAE程序)，從而提高計算效率；二是如何減少可靠性分析程序對有限元軟件的調用次數(shù)，從而降低計算成本。本文結合上述的基本思路，并針對如何在可靠性分析中高效集成CAE程序且有效減少調用次數(shù)問題，給出引信MEMS微結構可靠性仿真和實現(xiàn)方法。

1 主要問題及其解決思路

根據結構可靠度理論[9]，引信MEMS微結構失效概率模型可以寫為：

Pf=P(g(x)<0)

(1)

式(1)中，P(*)表示事件*發(fā)生的概率；g(x)是表征微結構狀態(tài)的函數(shù)，常稱為功能函數(shù)，g(x)<0表示微結構失效，g(x)>0表示微結構安全，g(x)=0表示極限狀態(tài)；x(x1,…,xn)T為隨機變量向量。對應的可靠度模型為：

R=1-Pf=P(g(x)≥0)

(2)

在實際應用中，功能函數(shù)g(x)可以進一步改寫為

g(x)=s(x)-y(x)

(3)

式中，s(x)為廣義強度，可以是微結構材料強度、允許最大變形量、允許最大應力值等；y(x)為廣義應力，可以是微結構應力響應、應變響應等；二者都可以是隨機變量向量x的函數(shù)。

將式(3)代入式(1)可以得到

Pf=P(s(x)-y(x)<0)

(4)

本文主要討論基于有限元仿真獲得廣義應力y(x)的情況。針對式(4)所示的失效概率模型的求解，需要解決兩個方面問題：一是如何在計算中集成CAE程序，方便地將可靠性分析的過程數(shù)據(如變量樣本值)傳遞給CAE程序作為輸入數(shù)據，又可以及時獲得CAE程序的求解結果，作為可靠性分析下一步驟的輸入數(shù)據；二是如何用盡可能少的CAE程序調用次數(shù)，完成預定的可靠性分析功能。

針對第一個問題，本文將以ANSYS軟件為例說明CAE集成方法，包括命令流文件更新、軟件調用、求解結果獲取等。針對第二個問題，本文將構建y(x)的近似模型，再利用蒙特卡洛法模擬法或迭代方法(如一次二階矩法及其改進方法等)進行可靠性求解。其中，近似建模所需要的樣本數(shù)據由試驗設計方法獲得。

2 可靠性仿真方法

2.1 仿真流程

結合上述問題及其相應的解決思路，可以得到基于有限元仿真的引信MEMS微結構可靠性仿真流程，如圖1所示。

圖1 引信MEMS微結構可靠性仿真流程Fig.1 Reliability simulation procedure for micro-structure in fuze

圖1中，需要結合引信MEMS微結構失效判據，以及應力-強度干涉的基本理論，得到功能函數(shù)g(x)的基本表達形式。常見的功能函數(shù)的基本形式如表1所示。

表1 常見的功能函數(shù)的基本形式

2.2CAE集成方法

在試驗設計中，采用全因子設計或部分因子設計等方法[10-11]獲得x的樣本數(shù)據，針對每一組樣本數(shù)據，調用一次CAE程序進行求解，獲得y(x)的一組輸出數(shù)據。對所有的x樣本數(shù)據進行同樣的操作，可以得到(x,y(x))的樣本數(shù)據。以ANSYS為例，給出各組樣本計算的流程，如圖2所示。

圖2 集成ANSYS的樣本計算流程Fig.2 Samples calculation procedure based on ANSYS

可以看出，在每一次的ANSYS調用過程中，需要完成三個核心的步驟：

1)依據i組樣本數(shù)據xi修改ANSYS輸入文件中的變量值，即將命令流中的“*set,M,5”與“*set,N,10”分別修改為“*set,M,5.225”和“*set,N,10.125”。在具體實現(xiàn)時，可以通過搜索關鍵字符串(如“*set,M”或“*set,N”)的方式，獲取待處理的數(shù)據文本在命令流文件中的位置，然后用新的樣本數(shù)據替換之，形成新的命令流文件。

2)以批處理方式(Batch Mode)調用ANSYS進行求解，這是采用后臺調用的方式，避免ANSYS軟件界面的反復閃現(xiàn)問題。在具體實現(xiàn)時，若采用Windows系統(tǒng)，那么可以將調用命令寫入后綴名為bat的文件，然后執(zhí)行該文件，即可實現(xiàn)ANSYS的后臺調用。

3)獲取本次ANSYS計算結果，如最大應力值或應變值等。在具體實現(xiàn)時，與命令流文件更新方法類似，可以通過搜索關鍵字符串(如“ITEM=MAX VALUE”)，獲取待提取數(shù)據文本所謂的位置，然后依據該位置信息獲取相應的結果數(shù)據。

其他的CAE軟件，如ADAMS、HyperWorks等，盡管各自的輸入、輸出文件格式，以及后臺調用命令的都不盡相同，但是都可以采用這一方式進行集成。

2.3 近似建模

為了便于后續(xù)的可靠性分析，在獲得(x,y(x))樣本數(shù)據的基礎上，需要構造y(x)的近似模型。常用的近似建模方法有曲線(面)擬合法和插值法，其中擬合法包括最小二乘法、正交多項式法等，插值法包括Kriging法[11-13]、徑向基函數(shù)插值等?？紤]到后續(xù)算例將采用Kriging法建立近似模型，下面簡要介紹該方法的基本原理，其他方法可以參考有關的文獻。

假設系統(tǒng)的響應值與自變量之間的真實關系可以表示成如下的形式：

y(x)=fT(x)b+z(x)

(5)

式中，fT(x)=[f1(x),f2(x),…,fp(x)]T，fi(x)為回歸模型的基函數(shù)，通常取為多項式函數(shù)，b=(b1,…,bp)T，bi為回歸系數(shù)，i=1,…,p；z(x)為一個平穩(wěn)的高斯過程，均值為0，協(xié)方差為：

cov(z(xi),z(xj))=σ2R(xi,Xj)

(6)

式中，i,j=1,2,…,m，m為試驗次數(shù)；R(·,·)為相關函數(shù)；σ2為隨機過程的方差值。

依據文獻[12]，可以得到b和σ2的估計值

(7)

式中，F(xiàn)是由樣本點得到的p×m矩陣，F(xiàn)=[fT(x1) …fT(xm)]；Y是各樣本點對應的輸出，即Y=[y(x1),…,y(xm)]T；R=[Rij]m×m為樣本點的相關性矩陣，Rij表達式如下：

(8)

在利用樣本數(shù)據獲得未知參數(shù)的估計值后，Kriging的近似模型可以表達為：

(9)

2.4 可靠度與靈敏度分析

結合式(3)和式(9)，可以g(x)得到的近似函數(shù)

(10)

(11)

式中，β為可靠度指數(shù)，滿足R=Φ(β)；y為x的等效標準正態(tài)隨機變量向量，一般通過構造Nataf分布，并進行相應變換得到[14]；L為y的相關系數(shù)矩陣經過Choleskey分解得到的下三角矩陣；u是與y相對應的獨立標準正態(tài)隨機變量向量，滿足u=L-1y。

參數(shù)靈敏度的計算表達如下。

(12)

式中，θ為隨機變量x對應的均值或標準差向量；T(x)為u與x的變換函數(shù)[13]；σ=[σij]n×n，當i=j時，σij為第i個隨機變量的標準差，當i≠j時，σij。

由式(11)和式(12)可以看出，靈敏度計算的核心在于偏導數(shù)?β/?u的求解。若采用不同的可靠度計算方法時，則該偏導數(shù)的計算也有所不同。當采用一次二階矩等迭代方法時，由于計算過程含有梯度的求解，因此?β/?u可以直接獲取。當采用蒙特卡洛等抽樣方法時，則需要結合可靠度計算公式作進一步推導，得到?β/?u的計算結果[15-17]。另外，如果要獲得失效概率或可靠度關于隨機變量的靈敏度，則可以通過二者與可靠度指數(shù)的函數(shù)關系，利用復合函數(shù)求導獲得。

3 可靠性仿真實現(xiàn)

根據上述的方法和流程，可以確定可靠性仿真實現(xiàn)的總體程序框架(如圖3所示)，主要包括：仿真工作流、試驗設計、近似建模和可靠性分析四個模塊，分別命名為Graph、DOE、Model和REL模塊。

圖3 可靠性仿真程序框架Fig.3 Program framework of reliability simulation

下面按模塊分析其主要功能及核心函數(shù)的框架代碼：

1)仿真工作流模塊。主要用于定義單次仿真的流程，它包含：輸入變量的定義，主要包括隨機變量的分布類型、分布參數(shù)等；中間過程的定義，包括CAE輸入文件、CAE調用命令、CAE求解結果獲取等；輸出變量的定義，包括廣義應力s(x)和功能函數(shù)g(x)。該模塊核心功能為單次的仿真計算，定義相應的計算函數(shù)為Graph.Execute()，框架代碼如下：

X = GetSample() '獲取隨機變量樣本值

InputFile = GetCAEInputFile() '獲取CAE輸入文件

InputFile = RefreshFile(X) '用X數(shù)據更新文件

StartCAE(InputFile) '調用CAE進行求解

ResultFile = GetCAEOutputFile() '獲取CAE結果文件

S = GetResult(ResultFile) '獲取最大應力等結果數(shù)據

G = Calculate(S) '計算功能函數(shù)值

2)試驗設計模塊。該模塊一方面需要集成各類試驗設計的算法用于生成輸入變量的樣本數(shù)據，另一方面需要與仿真工作流模塊進行交互，將各組樣本數(shù)據發(fā)送給仿真工作流進行求解，然后獲取相應的計算結果，最后得到全部的樣本數(shù)據。該模塊核心功能為輸入輸出變量樣本數(shù)據的獲取，定義相應的計算函數(shù)為DOE.Execute()，框架代碼如下：

Xs = CreateSamples() '依據算法生成樣本數(shù)據

For i=0 to (Xs.Length-1)

'調用工作流進行仿真計算

Y[i]=Graph.Execute(X[i]) '第i組樣本數(shù)據

End

return (X,Y)

3)近似建模模塊。該模塊一方面需要集成各類近似建模的算法，另一方面以試驗設計模塊的樣本數(shù)據為基礎，構建近似模型。該模塊核心功能為建模與預測，其中與建模相對應的計算函數(shù)定義為Model.Execute()，框架代碼如下：

(Xs,Ys)= DOE.Execute () '獲取DOE的樣本數(shù)據

model = SelectModel() '選擇模型的類型

beta = GetBeta(model, Xs, Ys) '依據樣本計算待定系數(shù)

與模型預測相對應的函數(shù)定義為Model.Evaluate()，框架代碼如下：

X=GetSample()'獲取輸入變量的一組樣本數(shù)據

Y=Predict(model, beta, X) '預測

return Y '返回預測結果

4)可靠性分析模塊。該模塊一方面需要集成各類可靠性算法(包括蒙特卡洛法、一次二階矩法等)，另一方面需要與近似建模模塊交互獲得近似模型，并將其用于可靠性分析。該模塊核心功能為可靠性分析，定義相應的計算函數(shù)為REL.Execute()，框架代碼如下：

f=GetFunction()'獲取功能函數(shù)基本形式

m=GetModel()'獲取功能函數(shù)的近似模型

x=GetRandom()'獲取隨機變量向量

R=Solve(f, m, x) '進行可靠度計算

S=GetSensitivity()'進行靈敏度計算

4 算例

某引信MEMS安全保險機構設計方案如圖4所示。彈丸發(fā)射后，當離心轉速w達到30 000 rpm時，在離心力的作用下，離心保險解除。 彈丸繼續(xù)飛行一定距離或時間后，根據有關傳感器信息解除電保險(電保險閉鎖簧片被推開)，滑塊在離心力作用下，克服平面彈簧抗力向卡座位置滑動，直至鎖入卡座。至此，傳爆序列對正，引信處于待發(fā)狀態(tài)。

圖4 引信MEMS安全保險機構示意圖Fig.4 Schematic diagram of micro-fuze safety mechanism

通過仿真分析發(fā)現(xiàn)閉鎖結構的卡頭是該安全保險機構的薄弱環(huán)節(jié)，其與卡座完成閉鎖的瞬間出現(xiàn)最大的彎曲應力，如圖5所示。針對這一薄弱環(huán)節(jié)，綜合考慮外部載荷、材料參數(shù)和尺寸參數(shù)等的影響，進行閉鎖結構的可靠性仿真分析。

圖5 閉鎖結構仿真分析結果Fig.5 Locking structure simulation result

1)仿真工作流建立

建立如圖6所示的單次仿真分析的工作流，為后續(xù)的試驗設計等提供輸入。其中采用HyperMesh軟件進行仿真模型的建立，利用LS-Dyna進行閉鎖過程的仿真分析，再由LS-PrePost讀取仿真分析的結果數(shù)據。

圖6 閉鎖結構仿真工作流Fig.6 Locking structure simulation workflow

圖6中的各節(jié)點的說明見表2，其中沖擊速度Vel由離心轉速w和卡頭質量換算而來；箭頭表示仿真計算的順序。

表2 仿真工作流各節(jié)點說明

2)試驗設計與近似建模

首先，選擇width 、E 、Vel 和thick為試驗設計的輸入變量， maxStress為試驗設計的輸出變量。輸入變量的取值范圍為[μ-3σ,μ+3σ]，其中μ和σ分別為隨機變量的均值和標準差。采用三水平全因子(三水平為：μ-3σ,μ,μ+3σ)試驗設計，在4個輸入變量的情況下，共計需要進行81次仿真試驗，總耗時約36 min。

然后，選擇零階多項式函數(shù)為回歸模型的基函數(shù)、Gauss函數(shù)為相關性函數(shù)的核函數(shù)，進行Kriging模型的建立。取width=0.058 mm、thick=0.25 mm時，maxStress與E、Vel的Kriging曲面如圖7所示(圖中黑點為樣本點)。

同時，增加五組樣本數(shù)據用于檢驗Kriging近似模型，獲得該模型的MSE(Mean Squared Error均方差)為1.790 5×10-4GPa，表明其具有較好的精度。

圖7 maxStress與E、Vel的Kriging曲面Fig.7 Kriging surface of maxStress, E and Vel

3)可靠性分析

將概率表達式P(Strength-maxStress>0)中的maxStress替換為上述的Kriging模型，再利用一次二階矩法進行可靠度和靈敏度的計算，結果分別如表3和表4所示。同時，為驗證可靠度結果的正確性，基于式(4)所示的可靠性模型(不采用近似模型的情況)，采用蒙特卡羅法對進行抽樣計算。為減少標準蒙特卡羅法的抽樣次數(shù)(抽樣次數(shù)106的情況下，總耗時將達305天)，采用文獻[18]的重要抽樣法以較少的抽樣次數(shù)獲得收斂的可靠度結果。

表3 可靠度結果

表4 靈敏度結果

(注：重要性靈敏度和均值靈敏度均已做歸一化處理)

從表3可以看出，采用近似模型的分析方法可以獲得較高精度的可靠度結果。在安全保險機構作用可靠度要求為0.99，且不考慮離心保險、電保險失效的情況下，閉鎖結構可靠度滿足要求。由于不考慮隨機變量之間的相關性，表4所示的重要性靈敏度和均值靈敏度結果是一致的。其中，材料強度，以及閉鎖結構寬度和厚度參數(shù)對于可靠性結果均呈“正”作用，即提高這些參數(shù)的均值，將有利于結構可靠性水平的提高；而材料和速度參數(shù)則呈“反”作用。且各變量的重要性水平可以從量化的數(shù)值大小進行比較。由標準差靈敏度結果可以看出，縮小各變量的標準差均有助于提高結構可靠度。

另外，在具體應用中，若單次有限元計算時間較長，那么采用蒙特卡羅法獲得可靠性結果是不可行的。因此，在利用本文方法進行可靠性分析時，建議從兩個方面保證結果的準確性：一是采用具有較高擬合精度的近似模型；二是采用經過驗證的可靠度和靈敏度計算方法。

5 結 論

本文提出了基于有限元法的引信MEMS微結構可靠性仿真分析及實現(xiàn)方法。該方法可以在結構性能仿真分析的基礎上，進一步開展量化的結構可靠性分析。其中，為實現(xiàn)可靠性分析程序對有限元軟件工具的自動調用，以ANSYS軟件為例給出CAE程序的集成方法和實現(xiàn)步驟；為降低可靠性分析程序對有限元仿真的調用次數(shù)，采用試驗設計和近似建模相結合的方法，構建功能函數(shù)的近似函數(shù)。提出可靠性仿真實現(xiàn)的總體程序框架，給出各主要模塊的功能及其主要函數(shù)的框架代碼。實際算例表明了該可靠性仿真方法的可行性和實用性。

[1]景華, 牛蘭杰, 宋永強. 小口徑榴彈引信微機電安全系統(tǒng)設計與仿真方法[J]. 彈箭與制導學報, 2010, 30(2): 129-132.

[2]李曉杰, 牛蘭杰, 翟蓉, 等. S型MEMS平面微彈簧垂直彈性系數(shù)[J]. 探測與控制學報, 2010, 32(1): 57-60.

[3]田中旺, 趙旭, 宋永強, 等. 改善強度的MEMS隔爆裝置懸臂梁和卡頭[J]. 探測與控制學報, 2011, 33(5): 1-4.

[4]孫誠誠, 牛蘭杰, 趙旭, 等. 引信桿狀彈簧質量塊微機電離心保險裝置[J]. 探測與控制學報, 2014, 36(1): 17-20.

[5]孫彩云，張康. 機電引信安保機構解除保險過程動態(tài)仿真[J].四川兵工學報, 2014(6):111-114.

[6]Rohit Pathak, Satyadhar Joshi. Optimizing reliability modeling of MEMS devices based on their application [J]. Work Journal ofModeling and Simulation, 2011(2): 139-154.

[7]Adam Martowicz, Tadeusz Uhl. Reliability- and performance-based robust design optimization of MEMS structures considering technological uncertainties [J]. Mechanical Systems and Signal Processing, 2012, 32: 44-58.

[8]Viviana Mulloni, Francesco Solazzi, Giuseppe Resta, etc. RF-MEMS switch design optimization for long-term reliability [J].Analog Integrated Circuits and Signal Processing, 2014, 78(2): 323-332.

[9]秦權, 林道錦, 梅剛. 結構可靠度隨機有限元-理論及工程應用[M]. 北京：清華大學出版社, 2006.

[10]Akira Todoroki，Tetsuya Ishikawa．Design of experiments for stacking sequence optimizations with genetic algorithm using response surface approximation [J]．Composite Structures，2004，64：349-357.

[11]Sei-Ichiro Sakataa, Fumihiro Ashidaa, Masaru Zakob. Approximate structural optimization using kriging method and digital modeling technique considering noise in sampling data [J]. Computers & Structures, 2008, 86(13): 1477-1485.

[12]Soren N Lophaven, Hans Bruun Nielsen, Jacob Sondergaard. DACE A Matlab Kriging Toolbox：Report IMM-REP-2002-13[R] Informatics and Mathematical Modeling， Technical University of Denmark, 2002. http://www.imm.dtu.dk/～hbn/dace/.

[13]Irfan Kaymaz. Application of kriging method to structural reliability problems [J]. Structural Safety, 2005, 27: 133-151.

[14]Ditlevsen O, Madsen HO. Structural Reliability Methods [M].West Sussex, England: John Wiley & Sons Ltd., 2007: 94-101.

[15]Melchers RE, Ahammed M. A fast approximate method for parameter sensitivity estimation in Monte Carlo structural reliability. Computers & Structures, 2004, 82: 55-61.

[16]Wu YT, Sitakanta M. Variable screening and ranking using sampling-based sensitivity measures[J]. Reliability Engineering & System Safety, 2006, 91: 634-47.

[17]Jensen H A, Mayorga F, Papadimitriou C. Reliability sensitivity analysis of stochastic finite element models [J]. Computer Methods in Applied Mechanics and Engineering, 2015, 296: 327-351.

[18]Frank Grooteman. Adaptive radial-based importance sampling method for structural reliability [J]. Structural Safety, 2008, 30: 533-542.

ReliabilitySimulationMethodforMEMSMicro-structureinFuze

TU Hongmao1,2, SUN Zhili1, JI Guangzhen3, LIU Qin1,2

(1. College of Mechanical Engineering and Automation, Northeastern University, Shenyang 110819, China； 2. Ordnance Science and Research Academy of China, Beijing 100089, China)

Reliability simulation method for MEMS micro-structure in fuze and its implementation details for further quantitative reliability assessment based on the corresponding performance simulation were proposed in this paper. By taking ANSYS as an example, the CAE program integration method and its implementation procedure were presented to facilitate the finite element tools utility during reliability simulation. It used design of experiment method and approximation modelling to reduce the number of times of finite element simulations in the whole reliability analysis procedure. Then, basic principles of reliability and sensitivity calculation issues were given. Programme framework for the reliability simulation method and its main code block were also presented. An engineering example demonstrated the feasibility and practicality of the proposed method were feasible and practical.

fuze reliability; MEMS micro-structure; reliability simulation; finite element method; reliability

2017-01-21

:國防技術基礎項目資助(Z092014B001)

:涂宏茂(1981—)，男，福建泉州人，博士研究生，副研究員，研究方向：MEMS可靠性。E-mail：bjthm@126.com。

TB114.3

：A

：1008-1194(2017)04-0005-07