薄壁加筋殼強度變差系數(shù)仿真預(yù)估方法研究

徐衛(wèi)秀，楊 帆，蔣亮亮，吳會強，郝 鵬

(1.北京宇航系統(tǒng)工程研究所，北京 100076；2.大連理工大學(xué)工程力學(xué)系工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，大連 116024)

0 引言

強度變差系數(shù)是結(jié)構(gòu)可靠性設(shè)計與分析的重要參數(shù)，是開展可靠性安全系數(shù)量化設(shè)計、結(jié)構(gòu)可靠性分析評估的重要依據(jù)。對運載火箭這類大型結(jié)構(gòu)，受成本限制，通過試驗獲得整個部段的強度數(shù)據(jù)很少，一般只有1～3個，很難獲得其強度變差系數(shù)?，F(xiàn)階段運載火箭結(jié)構(gòu)可靠性設(shè)計、評估參考的仍為20世紀(jì)六七十年代的數(shù)據(jù)。半個多世紀(jì)以來，隨著材料體系、制造加工工藝、產(chǎn)品測量、設(shè)計分析等技術(shù)的發(fā)展，影響運載火箭結(jié)構(gòu)強度散差因素的不確定性和獲取能力都發(fā)生了相應(yīng)變化，非常有必要基于運載火箭結(jié)構(gòu)研制設(shè)計、制造加工工藝現(xiàn)狀應(yīng)用新技術(shù)研究預(yù)估強度變差系數(shù)新方法。

薄壁結(jié)構(gòu)理論強度由材料強度、剛度性能，結(jié)構(gòu)尺寸參數(shù)確定，實際工程中，受制造、加工過程材料性能、結(jié)構(gòu)尺寸、幾何形貌隨機不確定性影響，同一結(jié)構(gòu)不同產(chǎn)品的強度有一定的散布。材料性能、結(jié)構(gòu)尺寸的不確定性較容易通過真實測量數(shù)據(jù)統(tǒng)計表征，對金屬材料強度、彈性模量數(shù)據(jù)及加工的結(jié)構(gòu)尺寸一般能較好地符合正態(tài)或近似正態(tài)分布[1-2]。幾何形貌的偏差又常稱為初始幾何缺陷，是引起受壓薄壁結(jié)構(gòu)屈曲失效極限強度分散性的關(guān)鍵因素[3-4]，一直是相關(guān)行業(yè)研究的熱點[5-9]。幾何缺陷測量技術(shù)最早可追溯到20世紀(jì)50年代[10]，經(jīng)歷了從接觸式到非接觸式測量方法的發(fā)展，目前應(yīng)用廣泛的為非接觸式光學(xué)形貌測量技術(shù)[11-13]，表征圓柱殼缺陷應(yīng)用最廣泛的為縱向半波、環(huán)向整波雙重傅立葉級數(shù)分解法[11，14，20]。2007年起，NASA[15-18]歷時十余年開展了殼體屈曲折減系數(shù)(SBKF)項目研究，基于考慮初始幾何缺陷的現(xiàn)代精確仿真分析技術(shù)，研究開發(fā)了現(xiàn)代運載火箭結(jié)構(gòu)屈曲折減系數(shù)和設(shè)計分析技術(shù)，并以8 in(1 in=25.4 mm)縮比和27.5 in全尺寸貯箱網(wǎng)格加筋結(jié)構(gòu)為例，開展了相關(guān)試驗驗證研究。我國大連理工大學(xué)王博教授團隊聯(lián)合北京宇航系統(tǒng)工程研究所也開展了相應(yīng)技術(shù)研究[19-20]，對1 m直徑整體機加成形光筒殼的初始幾何缺陷進行測量，并引入到理想形貌有限元模型，實現(xiàn)了考慮初始幾何缺陷的精確仿真分析技術(shù)，取得了很好效果。

在實現(xiàn)結(jié)構(gòu)強度精確仿真基礎(chǔ)上，預(yù)估強度變差系數(shù)的關(guān)鍵，是實現(xiàn)結(jié)構(gòu)極限強度的模擬打靶計算。對存在明顯非線性的薄壁結(jié)構(gòu)后屈曲強度問題，直接采用有限元分析方法開展極限強度模擬打靶效率很低。對于復(fù)雜系統(tǒng)的不確定性問題，為簡化計算，常采用基于概率分析的隨機擴展法建立響應(yīng)和變量不確定性間的代理模型，如混沌多項式展開法(PCE)、Karhumen-Loeve(KL)變換法等[21]，來實現(xiàn)基于代理模型的快速仿真打靶。

本文基于運載火箭結(jié)構(gòu)研制現(xiàn)狀和當(dāng)前設(shè)計手段，提出了基于產(chǎn)品數(shù)據(jù)的強度變差系數(shù)仿真預(yù)估技術(shù)。以當(dāng)前應(yīng)用廣泛、性能穩(wěn)定，從材料體系、制造加工工藝、設(shè)計分析手段等方面都得到改善的運載火箭金屬薄壁網(wǎng)格加筋殼為對象，開展相關(guān)實現(xiàn)技術(shù)研究。提出了實現(xiàn)技術(shù)途徑；結(jié)合試驗產(chǎn)品研制，測量統(tǒng)計材料性能、結(jié)構(gòu)尺寸和幾何形貌的不確定性；建立參數(shù)化有限元模型，實現(xiàn)考慮初始幾何形貌的結(jié)構(gòu)軸壓極限強度精確仿真分析方法；采用多項式混沌展開代理模型實現(xiàn)薄壁結(jié)構(gòu)極限強度的快速仿真打靶和強度總體變差系數(shù)的預(yù)估；并開展了3件1 m直徑網(wǎng)格加筋圓筒殼軸壓強度破壞試驗驗證研究。該技術(shù)研究旨在利用當(dāng)前結(jié)構(gòu)精細化設(shè)計分析能力，拓展運載火箭結(jié)構(gòu)產(chǎn)品可靠性相關(guān)數(shù)據(jù)獲取途徑，實現(xiàn)對結(jié)構(gòu)可靠性安全系數(shù)的精確量化設(shè)計，挖掘結(jié)構(gòu)承載潛力。

1 試驗件設(shè)計

試驗件為圖1所示的1 m直徑正置正交網(wǎng)格加筋殼。結(jié)構(gòu)設(shè)計時，綜合考慮蒙皮厚度加工工藝最小尺寸和加工精度保證能力，為避免尺寸過小引入的蒙皮厚度不確定性過大，首先確定蒙皮厚度為2.0 mm；其他網(wǎng)格尺寸參數(shù)根據(jù)軸壓載荷作用下的蒙皮局部屈曲、整體屈曲和材料屈服3種失效模式的等強度優(yōu)化設(shè)計確定。殼體長度在滿足中長殼要求的基礎(chǔ)上，考慮上下對接需要確定。原材料和成型加工工藝，為2A14T6鋁合金環(huán)軋鍛造成型，五坐標(biāo)機械銑加工。為實現(xiàn)統(tǒng)計驗證，策劃了3件產(chǎn)品，分別用SJ1-SJ3表示。

圖1 正交網(wǎng)格加筋殼幾何尺寸

2 研究內(nèi)容及實現(xiàn)技術(shù)途徑

研究內(nèi)容及實現(xiàn)技術(shù)途徑如圖2所示，主要包括：

圖2 實現(xiàn)技術(shù)流程圖

1)產(chǎn)品本體材料性能、主要結(jié)構(gòu)尺寸及初始幾何形貌測量及不確定統(tǒng)計表征；

2)含初始幾何缺陷參數(shù)化有限元模型的建立，極限強度仿真及其與試驗結(jié)果的對比驗證；

3)參數(shù)敏感性分析，有限元仿真采樣，多項式代理模型系數(shù)求解及回歸擬合；

4)多項式代理模型模擬打靶、強度變差系數(shù)預(yù)估及相關(guān)性分析等。

3 參數(shù)測量及不確定性表征

3.1 材料性能、結(jié)構(gòu)尺寸的測量及表征

材料性能方面，考慮到鍛件三向性能的不同，測量了軸向極限強度Rm、屈服強度Rp0.2、彈性模量E及延伸率δ等參數(shù)；結(jié)構(gòu)尺寸方面，測量了殼體長度L、直徑D，網(wǎng)格蒙皮厚度ts、筋條高度h、縱筋和橫筋寬度twx和twy。強度和彈性模量的測量精度要求為0.01 MPa，延伸率測量精度為0.01%，尺寸測量精度為0.01 mm。

材料性能和結(jié)構(gòu)尺寸采用正態(tài)分布表征，先檢驗數(shù)據(jù)是否符合正態(tài)分布，將不符合正態(tài)分布的數(shù)據(jù)轉(zhuǎn)換為近似正態(tài)分布。統(tǒng)計獲得的主要變量分布參數(shù)如表1所示，材料性能中屈服強度變差系數(shù)最大，為0.049；結(jié)構(gòu)尺寸中蒙皮厚度變差系數(shù)最大，為0.025 5。

表1 主要變量分布參數(shù)

3.2 幾何形貌的測量及表征

幾何形貌測量、點云數(shù)據(jù)處理及表征方案如圖3所示。測量采用了機械臂加激光掃描頭方案。掃描前需先建立合理坐標(biāo)系，用于點云數(shù)據(jù)處理和將幾何缺陷正確引入到有限元模型中。掃描獲得百萬點云數(shù)據(jù)后，根據(jù)有限元模型節(jié)點規(guī)模適當(dāng)精簡，將坐標(biāo)轉(zhuǎn)換到有限元模型坐標(biāo)系；經(jīng)濾波降噪處理，與理論外形相比獲得初始幾何缺陷(離面位移)。

圖3 幾何形貌測量和表征方案

已有研究表明，表征圓柱殼初始幾何缺陷最好的方法為軸向采用半波余弦、環(huán)向采用整波的雙重傅立葉級數(shù)方法[11，14，20]，在柱坐標(biāo)系下可表示為

(Amncosnθ+Bmnsinnθ)

(1)

其中，w為離面位移，L為筒殼的名義高度；x為軸向坐標(biāo)，θ為環(huán)向坐標(biāo)(角度)；Amn和Bmn為傅立葉級數(shù)的系數(shù)。

圖4為測量獲得的3件試驗產(chǎn)品離面位移云圖，方向有凹有凸，大小隨空間分布比較隨機。SJ1產(chǎn)品偏差最大，最大值為蒙皮厚度的13/50，SJ2、SJ3產(chǎn)品相對較小，最大值不到蒙皮厚度的1/10。

(a)SJ1

4 考慮初始幾何缺陷的極限強度精確仿真

4.1 含幾何缺陷參數(shù)化有限元模型建立

將結(jié)構(gòu)幾何簡化為殼，應(yīng)用Python腳本建立參數(shù)化有限元模型，采用S4單元模擬。為捕捉到薄壁結(jié)構(gòu)后屈曲失效模式，采用顯式動力學(xué)非線性方法求解[22]；為保證分析結(jié)果的精確性，對網(wǎng)格大小及計算時間進行了收斂性分析。

在Python環(huán)境下，利用Abaqus內(nèi)置函數(shù)得到理想形貌有限元模型各節(jié)點坐標(biāo)，按傅立葉級數(shù)擬合好的離面位移曲面，對各節(jié)點坐標(biāo)進行偏移，實現(xiàn)形貌偏差的引入。包含初始幾何缺陷的有限元模型如圖5所示。

圖5 含幾何缺陷有限元模型(放大60倍)

4.2 結(jié)果對比驗證

仿真分析了每件試驗產(chǎn)品材料性能，結(jié)構(gòu)尺寸取實測數(shù)據(jù)統(tǒng)計平均值，考慮初始幾何缺陷前后的軸壓極限強度，用承載能力T表征。

為獲得試驗件失效過程，試驗采用了DIC數(shù)字圖像測量系統(tǒng)，應(yīng)用數(shù)字圖像技術(shù)拍攝結(jié)構(gòu)在不同載荷作用下的照片，通過對比照片中像素點灰度值的變化，觀測結(jié)構(gòu)形貌。

3件產(chǎn)品均為結(jié)構(gòu)整體失穩(wěn)破壞。仿真和試驗軸壓強度結(jié)果對比見表2。SJ1、SJ2、SJ3考慮初始幾何缺陷(真實形貌)后的仿真與試驗偏差分別為8.76%、5.21%、4.12%。用DIC數(shù)字圖像測量系統(tǒng)對失效后形貌保持相對完整的SJ1、SJ2進行了測量，如圖6所示，SJ1存在失穩(wěn)波形分布不均、部分區(qū)域未失穩(wěn)的情況，這可能是其試驗強度低，仿真與試驗偏差最大的原因。SJ2、SJ3的偏差接近，圖7為SJ2仿真和試驗軸壓載荷-位移曲線和失效模式對比，曲線斜率和失效模式都吻合。

表2 軸壓極限強度對比

(a)SJ1

圖7 SJ2軸壓載荷-位移曲線及失效模式對比圖

考慮到還存在有限元模型簡化、試驗不確定性等因素，未做修正的仿真結(jié)果相對試驗偏差最大值小于10%，可以認為實現(xiàn)了極限強度的精確仿真。

由3個試驗數(shù)據(jù)獲得強度樣本變差系數(shù)為0.006 3；仿真分析結(jié)果預(yù)估的強度樣本變差系數(shù)為0.016 9，覆蓋了試驗值。同時，根據(jù)3個樣本考慮真實形貌前后分析結(jié)果的折減，預(yù)估3件產(chǎn)品形貌不確定性引起的變差系數(shù)為0.014 3。

5 多項式打靶及強度總體變差系數(shù)預(yù)估

總體變差系數(shù)可通過逐個產(chǎn)品數(shù)據(jù)統(tǒng)計仿真積累方式獲取，但需要足夠的產(chǎn)品數(shù)量和相對較長的時間。這里假設(shè)結(jié)構(gòu)制造加工工藝比較穩(wěn)定，探索一種通過相關(guān)變量統(tǒng)計數(shù)據(jù)仿真打靶預(yù)估方法，主要考慮材料性能、結(jié)構(gòu)尺寸參數(shù)在統(tǒng)計數(shù)據(jù)范圍內(nèi)的不確定，不考慮其隨空間分布的變化，打靶獲得材料、尺寸偏差引起的變差系數(shù)上限。在此基礎(chǔ)上考慮幾個樣本形貌引起的變差系數(shù)，估算強度總體變差系數(shù)上限。

因薄壁結(jié)構(gòu)極限強度有限元分析需要考慮材料非線性、幾何非線性，仿真計算一次需要的時間較長，直接通過仿真打靶效率很低，提出了應(yīng)用混沌多項式代理模型快速模擬打靶計算方案。

混沌多項式擴展法(PCE)是工程上應(yīng)用廣泛的不確定性擬合分析方法。與常規(guī)的多項式擬合不同，混沌多項式擴展法將隨機變量轉(zhuǎn)換為標(biāo)準(zhǔn)正態(tài)分布變量，并采用Hermite正交基函數(shù)，具有良好的均方收斂特性，能更好地估計不確定系統(tǒng)響應(yīng)的變差特性[21]。具體實施流程如圖8所示。

圖8 薄壁結(jié)構(gòu)非線性屈曲強度打靶流程

5.1 變量對強度影響分析

考慮材料彈性模量、屈服極限、極限強度、蒙皮厚度、筋條高度、縱筋和橫筋寬度尺寸共7個參數(shù)的隨機變化特性，為了解各參數(shù)對強度的影響關(guān)系，分析了承載能力隨各參數(shù)在均值和正負3倍標(biāo)準(zhǔn)偏差(μ±3σ)范圍內(nèi)的變化情況(其他參數(shù)取均值)。分析表明，強度隨彈性模量E、屈服強度Rp0.2的變化存在較弱的非線性，其他參數(shù)基本為線性關(guān)系。

5.2 采樣計算

根據(jù)影響分析結(jié)果，確定彈性模量、屈服強度2個參數(shù)采用2階多項式，其他參數(shù)采用線性擬合。

根據(jù)n個隨機變量，p階Hermite多項式展開總項數(shù)P計算公式[21]

(2)

7個隨機變量，2個參數(shù)2階、5個參數(shù)1階的多項式為31項，按2倍樣本計算，需要62個采樣點。應(yīng)用拉丁超立方采樣方法生成7個參數(shù)的62組隨機變量，將每組變量代入理想形貌參數(shù)化有限元模型中，仿真獲得對應(yīng)的極限強度響應(yīng)，實現(xiàn)采樣。

5.3 多項式擬合及回歸分析

將隨機變量和對應(yīng)的強度響應(yīng)代入公式(3)[21]，即可得擬合多項式系數(shù)B*

B*=(XTX)-1XTY

(3)

式中，X為n×P維的回歸變量矩陣，Y為n×1維的系統(tǒng)響應(yīng)矢量。

進而可獲得擬合響應(yīng)Y*和殘差e

Y*=XB*

(4)

e=Y-Y*

(5)

從而可對回歸方差和殘差進行分析，檢驗代理模型的擬合精確度。

整個過程通過Matlab程序執(zhí)行，多項式的擬合優(yōu)度R2達到0.998 2，調(diào)整的R2達到0.997 9。

5.4 強度打靶及總體變差系數(shù)預(yù)估

應(yīng)用混沌多項式打靶10萬次獲得的極限強度分布如圖9所示，基本仍符合正態(tài)分布，說明結(jié)構(gòu)極限強度隨材料彈性模量、屈服強度參數(shù)的變化無明顯非線性。

圖9 10萬次打靶數(shù)據(jù)分布

預(yù)估的極限強度變差系數(shù)及相關(guān)性分析結(jié)果見表3。材料、尺寸不確定性引起的強度變差系數(shù)為0.042 8。其中與材料部分的相關(guān)性為77.83%，主要由屈服強度引起；尺寸部分的相關(guān)性為22.5%，主要由蒙皮厚度引起。

表3 強度變差與各參數(shù)不確定性的相關(guān)系數(shù)

在考慮材料、尺寸不確定性預(yù)估的強度變差系數(shù)0.042 8基礎(chǔ)上，綜合考慮形貌不確定性引起的強度變差系數(shù)0.014 3(3個樣本預(yù)估值)，可得到強度總體變差系數(shù)預(yù)估值為0.045 1。

5.5 結(jié)果分析

與試驗樣本變差系數(shù)相比，仿真預(yù)估的總體變差系數(shù)達到7倍以上，高于正態(tài)分布變差系數(shù)上限的放大系數(shù)3.12[4]，具有較高的可靠性。但相對歷史統(tǒng)計數(shù)據(jù)0.08～0.12[4]，有較大程度的降低。

仿真預(yù)估的強度總體變差系數(shù)中，因?qū)嶋H結(jié)構(gòu)產(chǎn)品材料、尺寸的不確定性是隨空間分布的場變量，這里簡化為參變量，認為產(chǎn)品各處的材料、參數(shù)尺寸相同，包含整個結(jié)構(gòu)一個或多個變量參數(shù)同時達到上下偏差的極限情況，這比工程實際情況更為惡劣，材料、尺寸不確定性打靶預(yù)估部分相對保守；但形貌部分只考慮了3個樣本數(shù)據(jù)，包絡(luò)不了產(chǎn)品總體情況，需要隨產(chǎn)品數(shù)量增加逐步積累完善。

6 結(jié)論及展望

本文針對運載火箭大型薄壁結(jié)構(gòu)強度變差系數(shù)獲取困難問題，實現(xiàn)了一種基于產(chǎn)品實測數(shù)據(jù)仿真預(yù)估結(jié)構(gòu)強度變差系數(shù)方法。試驗結(jié)果表明，預(yù)估的結(jié)構(gòu)強度變差系數(shù)具有較高可靠性的同時，相對歷史統(tǒng)計數(shù)據(jù)，有一定幅度降低，為降低安全系數(shù)、減小結(jié)構(gòu)質(zhì)量提供了有力數(shù)據(jù)支持。

文中極限強度仿真打靶僅考慮了材料性能、結(jié)構(gòu)尺寸的不確定性，后續(xù)可根據(jù)同類結(jié)構(gòu)統(tǒng)計的初始幾何缺陷幅值及分布特征，開發(fā)產(chǎn)品形貌不確定性隨機抽樣技術(shù)，實現(xiàn)考慮材料性能、結(jié)構(gòu)尺寸和幾何形貌3大類不確定性打靶，進一步提高強度變差系數(shù)仿真預(yù)估精度。