冷氣式氣體發(fā)生器爆破片極限破壞載荷的計算

單津暉，魏學(xué)哲，王 承

（1. 同濟(jì)大學(xué)汽車學(xué)院，上海 201804；2. 均勝汽車安全系統(tǒng)研發(fā)（上海）有限公司CAE部，上海 201707）

簾幕式安全氣囊可以在汽車發(fā)生側(cè)面碰撞或翻滾時和安全帶一起挽救司乘人員的生命，是汽車被動安全系統(tǒng)的重要組成部分。冷氣式氣體發(fā)生器將高壓冷氣存儲在發(fā)生器腔體內(nèi)，點火后火藥爆炸將起閥門作用的爆破片炸開，高壓冷氣快速釋放并充滿氣袋。充氣過程中沒有對人體有害的氣體產(chǎn)生，且高壓冷氣充滿氣袋后吸熱膨脹，對氣袋內(nèi)部氣體壓力保持十分有利，因此對氣體保壓時間有要求的防翻滾簾幕氣囊多選用冷氣式氣體發(fā)生器。

冷氣式氣體發(fā)生器內(nèi)存儲的是高壓縮比的混合氣體，爆破片是氣體唯一的泄放裝置，在未點火時爆破片必須有足夠的強(qiáng)度，而在點火時爆破片必須在設(shè)計破壞壓力載荷下按預(yù)設(shè)的破壞模式發(fā)生破壞，其破壞載荷及破壞方式都有嚴(yán)格的設(shè)計要求［1-2］。

在設(shè)計開發(fā)過程中，為了達(dá)到盡可能高的預(yù)設(shè)壓力并加快驗證速度，通常會采用水壓試驗獲取爆破片的極限載荷，即通過灌注高壓水直至爆破片破壞，試驗過程中記錄的最高水壓值即為極限壓力破壞載荷。根據(jù)實驗室數(shù)據(jù)，準(zhǔn)靜態(tài)水壓和發(fā)生器動態(tài)點爆加載條件下的爆破片破壞載荷基本一致，因此在設(shè)計開發(fā)階段對水壓試驗爆破片極限載荷的準(zhǔn)確預(yù)測是發(fā)生器結(jié)構(gòu)件設(shè)計的關(guān)鍵。

1 極限載荷的分析方法

確定極限載荷的方法主要分為試驗法、理論計算法和有限元法三大類，試驗法只能用于后期驗證，在設(shè)計早期主要采用理論計算法和有限元法。當(dāng)前針對爆破片的極限壓力破壞載荷相關(guān)計算理論并未統(tǒng)一，理論方法對材料非線性及幾何非線性的考慮較少，且無法考慮零件發(fā)生破壞前的局部失穩(wěn)，因此理論計算結(jié)果精度尚不能滿足工程開發(fā)需要。目前基于非線性彈塑性理論的有限元計算正逐漸成為獲取極限載荷的主要手段［3-5］。

在有限元分析中，破壞載荷極值的判定方法較多，如塑性功法、2 倍斜率法、雙切線相交法、零曲率法、弧長法等［6-11］。其他方法相比零曲率法、弧長法精度不足，在實際中應(yīng)用也不易掌握，本文也不再討論。

1.1 零曲率法

極限載荷分析屬于塑性力學(xué)問題，假定材料為理想彈塑性，結(jié)構(gòu)處于小變形狀態(tài)，服從Mises 屈服條件及相關(guān)流動準(zhǔn)則。在上述假定基礎(chǔ)上，在有限元分析中，通過多個載荷步、小的載荷增量逐步對結(jié)構(gòu)施加載荷。當(dāng)載荷增量逐漸趨近于零、即載荷對應(yīng)位移曲線的曲率趨近于零時，此時結(jié)構(gòu)的承載能力到達(dá)極限，再增加載荷會引起極限載荷下的結(jié)構(gòu)失穩(wěn)，即發(fā)生塑性垮塌。這種將載荷位移曲線上零曲率點對應(yīng)的載荷值定義為極限載荷的方法稱為零曲率法。

從力學(xué)分析角度看，極限載荷問題求解，就是通過不斷求解計入幾何非線性和材料非線性的結(jié)構(gòu)平衡方程尋求結(jié)構(gòu)極限荷載的過程。非線性結(jié)構(gòu)的增量平衡方程一般可以表示為［8-9］

式中：K為結(jié)構(gòu)整體切線剛度；ΔX為結(jié)構(gòu)位移增量向量；Δλ為載荷比例系數(shù)；F為載荷向量；R為不平衡力向量。

在極限載荷狀態(tài)下，載荷比例系數(shù)Δλ趨近于零，當(dāng)計算嘗試給Δλ一個系統(tǒng)設(shè)定的最小增量時，整體剛度矩陣K出現(xiàn)奇異，剛度矩陣變成“病態(tài)”矩陣，即系數(shù)矩陣的微小誤差將導(dǎo)致計算結(jié)果的巨大誤差，引起計算不收斂而導(dǎo)致計算終止。在有限元計算中，零曲率法將計算最后一步能收斂的載荷定義為極限載荷。

零曲率法基于隱式算法求解，計算位移增量ΔX需要計算結(jié)構(gòu)剛度K的逆矩陣，因此如果結(jié)構(gòu)非線性程度較高，尤其是含復(fù)雜接觸的大載荷極限工況，該方法收斂難度大。另外零曲率法的計算終止是以載荷增量步為零作為標(biāo)志，屬于計算非正常終止，這要求分析人員有較高的能力對結(jié)果做出判斷。

1.2 弧長法

弧長法是一種隱式非線性求解的迭代控制方法，不需要事先設(shè)定失穩(wěn)準(zhǔn)則或斷裂準(zhǔn)則，可以在荷載和位移增量均不確定的情況下，追蹤結(jié)構(gòu)加載路徑，避免傳統(tǒng)隱式算法的收斂問題［12-13］?；￠L增量法以載荷比例系數(shù)和位移增量平方和的開方作為增量參數(shù)進(jìn)行收斂計算分析，其增量參數(shù)為

在收斂計算過程中，弧長法以前一步增量計算得到的平衡點作為圓心，基于超平面搜索、求解半徑為ΔS的超曲面S與結(jié)構(gòu)平衡路徑的交點。在達(dá)到極限載荷后，增量弧長法在超曲面S內(nèi)可沿載荷增量Δλ負(fù)曲率方向搜索平衡狀態(tài)，通過產(chǎn)生負(fù)的載荷增量Δλ以匹配結(jié)構(gòu)失穩(wěn)后負(fù)定的剛度矩陣K，繼續(xù)追蹤結(jié)構(gòu)位移ΔX增加對應(yīng)載荷F的變化。

弧長法將載荷增量為正的最后一個迭代步對應(yīng)的載荷定義為極限載荷。其極限載荷求解精度與零曲率法基本相當(dāng)，但相比零曲率法，弧長法可以繼續(xù)求解結(jié)構(gòu)在失穩(wěn)狀態(tài)下剛度負(fù)定、載荷增量為負(fù)的力學(xué)平衡方程。這對研究整個結(jié)構(gòu)在達(dá)到極限載荷后應(yīng)力、應(yīng)變變化規(guī)律具有重要意義。

對每一個計算迭代步，弧長法都要對結(jié)構(gòu)剛度矩陣正定或負(fù)定進(jìn)行判定，并對載荷增量變化在正負(fù)曲率變化方向進(jìn)行搜索，因此相比只考慮載荷增量正曲率變化的零曲率法，弧長法的計算增量步長較小，每一迭代步的計算時間更長。

1.3 針對顯式算法的應(yīng)力矩陣狀態(tài)分析法

基于隱式算法的零曲率法和弧長法對計算收斂的要求較高，對需要考慮接觸的復(fù)雜結(jié)構(gòu)，在極限大載荷加載狀態(tài)下計算通常都很難收斂。對于瞬態(tài)高度非線性工況，通常采用顯式算法求解。t時刻結(jié)構(gòu)的動力學(xué)微分方程可表示為［14］

式中：M為質(zhì)量矩陣；C為阻尼矩陣；U?t為加速度向量；U?t為速度向量；Ut為位移向量。

若t時刻節(jié)點位移、速度加速度均為已知，使用中心差分法對加速度、速度的導(dǎo)數(shù)采用中心差分替代，式（3）可轉(zhuǎn)化為［15］

由式（4）、（5）可知，基于中心差分法的動力學(xué)微分方程組求解不需要像式（1）那樣需要求解結(jié)構(gòu)剛度K的逆矩陣；另外中心差分法將加速度U?和速度U?表示位移U的某種組合，求解t+Δt時刻的瞬時位移向量Ut+Δt時，只需要根據(jù)t+Δt時刻以前的狀態(tài)變量計算出M?和R?，即可直接計算出Ut+Δt。因此只要設(shè)定好計算時間步長Δt，顯式算法一般不存在收斂性問題，可有效解決各種復(fù)雜非線性問題。

顯式算法基于中心差分動力學(xué)方程（4），不是力學(xué)平衡方程（1）求解，對于極限載荷計算，在結(jié)構(gòu)達(dá)到承載極限后，只要有滿足收斂的時間步長Δt顯式計算就會繼續(xù)，施加的載荷也會按預(yù)先設(shè)定繼續(xù)增加，不會像隱式算法那樣出現(xiàn)計算終止或載荷下降的情況。因此對于基于顯式算法的極限載荷求解，工程上多采用固定約束位置支撐反力曲線是否出現(xiàn)拐點來判定極限載荷：當(dāng)結(jié)構(gòu)達(dá)到承載極限后發(fā)生失穩(wěn)，其固定約束位置的支撐反力不再隨載荷增加而增加，而因結(jié)構(gòu)失去承載能力開始下降，將固定約束位置支撐力反力曲線出現(xiàn)的極值拐點值定義為結(jié)構(gòu)的極限載荷，這種極限載荷定義方法的原理與弧長法類似，計算精度也能滿足工程需要。但對于載荷施加位置靠近固定約束位置的工況，就算結(jié)構(gòu)失穩(wěn)支撐反力也會隨載荷的增加而繼續(xù)增加，這給極限載荷的判定帶來了較大難度，也表明該判定方法有相當(dāng)?shù)木窒扌浴?/p>

結(jié)構(gòu)剛度矩陣K由結(jié)構(gòu)所有單元的剛度組合而成，類似式（1）單元剛度矩陣Ke是表征單元節(jié)點力向量及節(jié)點位移向量的關(guān)系矩陣。在結(jié)構(gòu)幾何確定的情況下，節(jié)點剛度矩陣主要取決于節(jié)點應(yīng)力矩陣σe的狀態(tài)。

對三維結(jié)構(gòu)上任一節(jié)點，其應(yīng)力矩陣σe為一含9個應(yīng)力分量的3×3 對稱矩陣，為了方便分析，通過矩陣線性變換將σe轉(zhuǎn)化為只含有3個應(yīng)力分量的特征值矩陣，轉(zhuǎn)化后的σe可表示為［16］

式中：σ1、σ2和σ3分別為節(jié)點第一、第二及第三主應(yīng)力。3 個主應(yīng)力σ1、σ2和σ3的取值及其方向決定了一點的應(yīng)力狀態(tài)。

應(yīng)力矩陣σe中各應(yīng)力分量大小取決于材料應(yīng)力-應(yīng)變（σ-ε）本構(gòu)關(guān)系。塑性材料σ-ε曲線可以表示為

式中：k為材料應(yīng)變硬化系數(shù)；n為材料應(yīng)變硬化指數(shù)；σT為材料拉伸極限應(yīng)力；εT為拉伸極限應(yīng)變。圖1為爆破片材料應(yīng)力-應(yīng)變曲線，在達(dá)到極限應(yīng)力前，σ-ε按圖中曲線實線段變化，在實線末端應(yīng)力達(dá)到極值、材料剛度最大后，若繼續(xù)增加載荷，對應(yīng)不同算法：①隱式算法（零曲率法）。σ-ε曲線剛度達(dá)到極值，計算因沒有大于極限載荷的收斂解而終止。②隱式算法（弧長法）。載荷增量變?yōu)樨?fù)值，應(yīng)變隨位移增加而增加，σ-ε曲線按圖1 虛線段變化。③顯式算法。載荷繼續(xù)增加，應(yīng)變隨位移增加繼續(xù)增加，σ-ε曲線按圖1虛線段變化。

對比后2種算法，在達(dá)到極限載荷后，由于材料本構(gòu)相同，其應(yīng)力-應(yīng)變對應(yīng)結(jié)構(gòu)位移變化的趨勢接近，關(guān)鍵位置的應(yīng)力矩陣變化趨勢也接近，使用應(yīng)力矩陣狀態(tài)分析法判定顯式算法的極限載荷，即通過分析、對比隱式弧長法和顯式算法2 種算法相同節(jié)點在極限狀態(tài)下的應(yīng)力矩陣變化規(guī)律來判定顯式算法的極限載荷。

2 爆破片極限載荷分析

2.1 基于零曲率法計算爆破片極限載荷

使用ABAQUS 軟件對爆破片的極限載荷進(jìn)行計算。爆破片材料牌號為AISI301，材料屈服極限1 100MPa，真實拉伸極限1 482MPa，真實斷裂應(yīng)變0.182，其材料曲線如圖1所示，如圖中實線部分，先按公式σ=1 711ε0.0831對應(yīng)變化，σ達(dá)到拉伸極限后不再增加，σ-ε曲線按圖中虛線段變化。

圖1 AISI301材料應(yīng)力-應(yīng)變曲線Fig.1 Strain versus stress of AISI301

發(fā)生器前端示意圖如圖2 所示，爆破片為正拱型爆破片，爆破片通過激光焊接與擴(kuò)散器連接，擴(kuò)散器和發(fā)生器殼體通過摩擦焊焊接成一個整體。爆破片與擴(kuò)散器焊接后局部截面圖如圖3所示。爆破片厚度0.203 2±0.012 7mm，仿真取厚度下限0.190 5mm。在實際試驗中，除爆破片外，其他零件變形很小，因此為簡化仿真，仿真中只考慮爆破片零件。為保證模型精度，對爆破片模型使用六面體實體網(wǎng)格建模。

圖2 冷氣式氣體發(fā)生器出氣端結(jié)構(gòu)示意Fig.2 Structure of outlet of cold gas inflator

圖3 焊接后爆破片區(qū)域截面示意Fig.3 Section view of burst disk after welding

為了提高仿真計算精度，在極限載荷計算中計入沖壓工藝對爆破片的影響，仿真分2步完成。第1步使用動態(tài)顯式分析模塊將處于自由狀態(tài)的爆破片圓片沖壓成正拱型。沖壓前后爆破片的厚度如圖4所示。沖壓后爆破片厚度沿四周向中心逐漸減薄，中心點位置厚度0.168 4mm，減薄率13.1%。含沖壓仿真產(chǎn)生的預(yù)應(yīng)力、厚度不一的爆破片仿真模型是第2 步仿真的起始條件，與無預(yù)應(yīng)力且厚度均一的爆破片模型相比無疑更符合真實情況。

圖4 沖壓前后爆破片截面形狀對比Fig.4 Section of burst disk before and after stamping

第2步仿真在爆破片周邊焊接位置添加固定約束，并在爆破片沿內(nèi)表面法向施加100Mpa 的壓力（設(shè)計破壞壓力大于82.6Mpa）。使用靜態(tài)隱式非線性分析模塊計算爆破片的極限載荷。在第21 個計算迭代步，載荷為83.5MPa 時仿真計算終止，計算過程中爆破片中心N133號節(jié)點位移始終最高，其位移對應(yīng)壓力載荷曲線如圖5 的虛線所示。由圖5 可知，在載荷較低時，載荷增量增加較快，當(dāng)載荷趨近于83.5Mpa 時，載荷增量逐漸減小并趨近于零。在載荷為83.5Mpa、N133 號節(jié)點位移為0.868mm 時載荷增量為零，此時曲線的斜率為零。根據(jù)零曲率法準(zhǔn)則可以判定爆破片的極限破壞載荷為83.5Mpa。

圖5 壓力載荷對應(yīng)爆破片最大變形量變化曲線Fig.5 Curves of pressure versus maximum displacement of burst disk

根據(jù)試驗室數(shù)據(jù)，15 次水壓試驗的平均破壞載荷為86.0MPa，最大值88.4Mpa，最小值83.1MPa。仿真計算獲得的極限載荷和試驗破壞均值載荷誤差為3.0%，與試驗最小值誤差為0.5%?？紤]到仿真中爆破片厚度取公差下限，這個極限載荷計算結(jié)果精度十分理想。

試驗后爆破片破壞形狀如圖6a 所示，其破壞從中心位置開始；極限載荷作用下爆破片第一主應(yīng)變云圖如圖6b 所示，爆破片中心位置第一主應(yīng)變最大。對比試驗結(jié)果，仿真對應(yīng)變最大位置的預(yù)測是準(zhǔn)確的。

圖6 爆破片試驗及仿真結(jié)果對比Fig.6 Test and simulation results of burst disk

2.2 基于弧長法計算爆破片極限載荷

2.2.1 使用弧長法計算爆破片極限載荷

使用弧長法計算爆破片的極限載荷，第1 步?jīng)_壓成形仿真和零曲率法一致，第2 步隱式非線性計算仿真除了零曲率算法已經(jīng)設(shè)定的幾何非線性選項外，再增加RIKS（弧長法）算法選項?；￠L法基于小載荷增量步計算，相比零曲率法其迭代步數(shù)增加較多，對于大極限載荷計算，需要預(yù)先設(shè)定一個較大的最大迭代步數(shù)，本算例設(shè)定值為300步。

在弧長法計算過程中，在第183 迭代步載荷達(dá)到峰值83.6Mpa，第184 步后載荷增量為負(fù)值，載荷也隨位移的增加而減小。爆破片中心N133 號節(jié)點位移對應(yīng)壓力載荷曲線如圖5 實線所示。圖5 實線在載荷為83.6MPa 時出現(xiàn)拐點，壓力達(dá)到峰值，此時對應(yīng)節(jié)點位移0.903mm，在此位移之前結(jié)構(gòu)一直加載，曲線曲率為正值，大于此位移后結(jié)構(gòu)失穩(wěn)、可承受的最大載荷也相應(yīng)下降，曲線曲率也變?yōu)樨?fù)值。

將圖5實線的極值83.6Mpa定義為弧長法計算獲得的極限載荷，其值和試驗破壞載荷均值的誤差為2.9%，最小值誤差為0.6%。其極限載荷計算結(jié)果和零曲率法基本一致，只有0.1%的差別。

2.2.2 弧長法爆破片極限載荷計算結(jié)果分析

對比分析零曲率法、弧長法計算結(jié)果，對圖5的2條載荷位移曲線進(jìn)行對比分析：

（1）在達(dá)到極限載荷前，2條曲線走勢基本重合，零曲率法載荷增量步較大，直到接近極限載荷時載荷增量才逐漸減小。

（2）從第1個迭代步開始，弧長法每一個迭代步的載荷增量都很小，其曲線變化趨勢完全體現(xiàn)出材料彈塑性本構(gòu)曲線變化和載荷變化的對應(yīng)關(guān)系。

（3）由圖5實線可知，在極限載荷作用下結(jié)構(gòu)失穩(wěn)后，弧長法可以繼續(xù)計算對應(yīng)位移增加爆破片結(jié)構(gòu)的應(yīng)力-應(yīng)變，這對研究達(dá)到極限載荷后爆破片應(yīng)力-應(yīng)變變化規(guī)律非常有意義。

由式（2）弧長法搜索、計算極限載荷的原理可知，弧長法不依賴于某個強(qiáng)度失效理論，而是基于結(jié)構(gòu)剛度矩陣正定或負(fù)定判斷結(jié)構(gòu)是否失穩(wěn)。因此在進(jìn)行應(yīng)力-應(yīng)變結(jié)果分析時，要對表征剛度的應(yīng)力矩陣σe進(jìn)行分析，即對式（6）中σ1、σ2和σ3及其對應(yīng)應(yīng)變變化規(guī)律進(jìn)行考量。

在仿真過程中，爆破片中心位置變形一直最大，對應(yīng)應(yīng)力、應(yīng)變也最大。對爆破片中心N133號節(jié)點進(jìn)行應(yīng)力-應(yīng)變分析，其受力狀態(tài)為兩拉一壓，沿2個徑向拉伸、沿厚度方向壓縮。

N133號節(jié)點第一、第二主應(yīng)力對應(yīng)位移變化曲線如圖7 所示，在位移為0.6mm 后，第一、第二主應(yīng)力增長趨勢變得非常緩慢，但始終低于材料極限應(yīng)力。N133 號節(jié)點第三主應(yīng)力絕對值對應(yīng)位移變化曲線如圖8 所示，在位移為0.903mm 時第三主應(yīng)力取極大值，后繼計算其值隨位移增加而減小。

圖7 N133號節(jié)點第一、第二主應(yīng)力對應(yīng)的位移變化曲線Fig.7 Displacement versus 1st and 2nd principle stress of Node 133

圖8 N133號節(jié)點第三主應(yīng)力對應(yīng)的位移變化曲線Fig.8 Displacement versus 3rd principle stress of Node 133

對N133 號節(jié)點的應(yīng)力、應(yīng)變分析表明，達(dá)到極限載荷時，該節(jié)點第一、第二主應(yīng)力接近材料拉伸極限，其位移對應(yīng)第三主應(yīng)力及壓力載荷2 條曲線同時產(chǎn)生拐點，這說明該節(jié)點，即爆破片中心位置沿厚度方向發(fā)生了失穩(wěn)，整個結(jié)構(gòu)的承載能力開始下降，爆破片將從中心位置為起點發(fā)生塑性破壞。

N133號節(jié)點3個方向主應(yīng)變對應(yīng)的位移變化曲線如圖9所示?？梢?，第一、第二主應(yīng)變變化曲線基本重合，第三主應(yīng)變絕對值為第一、第二主應(yīng)變之和；在極限載荷對應(yīng)位移0.903mm位置，3個主應(yīng)變變化都比較平緩，沒有明顯的拐點，這和圖1中應(yīng)變從實線段平緩過渡到虛線段的趨勢相匹配。

圖9 N133 號節(jié)點第一、第二、第三主應(yīng)變對應(yīng)的位移變化曲線Fig.9 Displacement versus 1st, 2nd, and 3rd principle strain of Node 133

2.3 基于應(yīng)力矩陣狀態(tài)分析法計算爆破片極限載荷

對于本算例，通過零曲率法和弧長法都獲得了精度較高的極限載荷值。但受計算方法限制，仿真中未考慮爆破片與下方焊接零件擴(kuò)散器的接觸，由圖3 可知，與擴(kuò)散器的接觸限制了爆破片對應(yīng)位置向下的移動，若忽略此接觸，爆破片的變形結(jié)果精度會有一定失真。

為考察接觸設(shè)定對爆破片極限載荷和變形結(jié)果的影響，進(jìn)一步提高仿真精度，使用DYNA 顯式分析軟件對爆破片的極限載荷進(jìn)行計算。分析過程同樣分為2 步，第1 步為沖壓仿真，第2 步在爆破片內(nèi)表面沿法向施加一個0～100MPa隨時間線性增加的壓力，同時約束爆破片焊接位置并設(shè)置爆破片與擴(kuò)散器零件的接觸。接觸類型為面對面接觸（surfaceto-surface），由于爆破片與擴(kuò)散器都是鋼材，材料剛度接近，對應(yīng)接觸柔度設(shè)定為類型1（Soft設(shè)為1）。

第2步仿真計算在載荷達(dá)到100MPa后按預(yù)先設(shè)定終止，在載荷為88.3MPa 時爆破片中心位置單元因變形過大被刪除，但在整個仿真過程中，爆破片焊接位置附近區(qū)域始終與其下方的擴(kuò)散器接觸并承受不斷增加的法向壓力，導(dǎo)致爆破片約束位置的支撐反力曲線始終保持持續(xù)增加的趨勢，沒有極值拐點出現(xiàn)。若將仿真中爆破片中心位置出現(xiàn)斷裂時對應(yīng)的內(nèi)表面法向壓力值88.3MPa確定為爆破片的極限載荷，其值與試驗結(jié)果上限值相當(dāng)，和試驗下限結(jié)果83.1MPa誤差為6.3%，計算精度低于隱式算法。

使用應(yīng)力矩陣狀態(tài)分析法對顯式算法的極限載荷進(jìn)行判定。與弧長法分析過程相同，對爆破片中心N133號節(jié)點應(yīng)力矩陣狀態(tài)變化進(jìn)行分析，其3個方向主應(yīng)力對應(yīng)位移變化曲線與弧長法對應(yīng)結(jié)果基本一致。限于篇幅，只對出現(xiàn)拐點的第三主應(yīng)力變化曲線進(jìn)行討論。

顯式算法計算獲得的N133 號節(jié)點第三主應(yīng)力絕對值對應(yīng)壓力載荷及位移變化曲線分別如圖10、圖11 所示。由圖10 可知，當(dāng)載荷為84.0MPa 時，N133號節(jié)點第三主應(yīng)力達(dá)到峰值，并隨載荷增加迅速下降。由圖11 可知，當(dāng)N133 節(jié)點第三主應(yīng)力絕對值在節(jié)點位移為0.590mm時達(dá)到峰值，然后隨位移增加其值逐漸減小。

圖10 N133號節(jié)點第三主應(yīng)力的應(yīng)力-位移曲線Fig.10 Pressure versus 3rd principle stress of Node 133

對比圖8、圖11弧長法和顯式算法N133號節(jié)點第三主應(yīng)力變化曲線，可知：

圖11 N133號節(jié)點第三主應(yīng)力對應(yīng)其位移變化曲線Fig.11 Displacement versus 3rd principle stress of Node 133

（1）弧長法和顯式算法計算結(jié)果曲線都出現(xiàn)極值拐點，表明在厚度方向爆破片中心位置出現(xiàn)了失穩(wěn)，爆破片結(jié)構(gòu)也將在此位置發(fā)生破壞。

（2）在極值點前，2條曲線變化趨勢基本相同，在極值點后由于載荷繼續(xù)增加顯式算法曲線下降趨勢相比載荷開始下降的弧長法要更快一些。即相比弧長法，顯式算法在厚度方向的失穩(wěn)隨位移增加會更迅速。

將顯式計算過程中爆破片中心位置失穩(wěn)即N133 號節(jié)點第三主應(yīng)力取極值時對應(yīng)的載荷84.0MPa 定義為爆破片的極限破壞載荷，其值與爆破片試驗破壞載荷均值誤差為2.4%，與試驗結(jié)果下限誤差為1.1%。

上述基于隱式弧長法和顯式中心差分算法計算結(jié)果對爆破片中心位置關(guān)鍵節(jié)點應(yīng)力矩陣變化規(guī)律進(jìn)行分析、對比，確定顯式算法極限載荷的方法即為應(yīng)力矩陣狀態(tài)分析法。對比隱式弧長法與針對顯式算法的應(yīng)力矩陣狀態(tài)分析法計算結(jié)果，可知：

（1）使用應(yīng)力矩陣狀態(tài)分析法和弧長法判定的極限載荷分別為84.0MPa和83.6MP，2種方法的匹配度達(dá)到99.5%，表明應(yīng)力矩陣狀態(tài)分析法和弧長法的判定的極限載荷值有極佳的一致性。

（2）顯式計算中考慮爆破片下方零件支撐的影響后，極限載荷作用下爆破片最大變形結(jié)果由0.902mm 減小到0.590mm，相差34.6%，因此若要準(zhǔn)確計算爆破片的變形量，必須使用基于顯式算法的應(yīng)力矩陣狀態(tài)分析法。

3 結(jié)論

首先使用基于隱式算法的零曲率法和弧長法進(jìn)行計算，確定冷氣式氣體發(fā)生器爆破片水壓試驗的極限破壞壓力載荷，考慮了爆破片的沖壓成形工藝后，和試驗結(jié)果相對比，這2種方法對爆破片極限載荷的計算誤差都在3%以內(nèi)。

對弧長法極限載荷計算過程中爆破片中心節(jié)點應(yīng)力矩陣狀態(tài)變化規(guī)律進(jìn)行分析，發(fā)現(xiàn)在極限載荷作用下，中心位置節(jié)點第三主應(yīng)力對應(yīng)位移曲線出現(xiàn)極值拐點，該拐點對應(yīng)載荷步與壓力位移曲線極值拐點對應(yīng)載荷步相同。這表明在達(dá)到極限載荷后，爆破片中心位置在厚度方向發(fā)生了失穩(wěn)。

為了彌補(bǔ)隱式算法接觸非線性難收斂的不足，基于弧長法應(yīng)力-應(yīng)變分析結(jié)果提出針對顯式算法的極限載荷判定方法，即應(yīng)力矩陣狀態(tài)分析法，它將顯式計算過程中爆破片中心在厚度方向發(fā)生失穩(wěn)時對應(yīng)的載荷定義為極限載荷。該方法在利用顯式算法在解決接觸非線性易收斂的優(yōu)勢同時，避免了結(jié)構(gòu)局部微小失穩(wěn)時顯式算法受力狀態(tài)不完全平衡帶來的誤差。最終使用應(yīng)力矩陣狀態(tài)分析法判定的爆破片極限載荷值與弧長法結(jié)果匹配度達(dá)到99.5%，與試驗結(jié)果相比誤差在2.4%以內(nèi)，且相比隱式算法該方法的變形結(jié)果更加準(zhǔn)確。

作者貢獻(xiàn)聲明：

單津暉：理論研究、提煉，論文撰寫。

魏學(xué)哲：論文統(tǒng)籌、規(guī)劃。

王 承：仿真，試驗協(xié)調(diào)。