某大型固體發動機T800碳纖維殼體封頭結構仿真分析和優化設計①

喻琳峰,任全彬,張愛華,宋學宇

(1.西安航天動力技術研究所,西安 710025;2.航天動力技術研究院,西安 710025)

0 引言

在實際工程應用中,國產T800碳纖維纏繞殼體的力學性能表現往往比不上傳統使用的進口T800碳纖維[1-2],甚至有時候僅僅能與進口T700碳纖維持平。特別是殼體封頭部位往往在壓力未達到設計爆破壓強時就會發生破壞,主要原因在于,封頭處于不同種類材料的交界處,會出現材料不同,剛性不同,對應力的敏感不同的情況。受壓時,由于形成應力的不連續傳遞,導致形成了局部的應力集中區,從而成為承力薄弱點,并首先發生破壞[3]。同時,殼體的赤道位置處由于環向層的加入形成了幾何形狀突變,導致了該位置存在較高的應力。在傳統設計方法以網格理論為基礎的分析計算中,既不考慮基體樹脂的作用,也未考慮纏繞層次的工藝影響,更不可能考慮結構內部缺陷以及纖維與基體界面作用的影響。在實際工程中靠近赤道位置和臨近金屬接頭位置的纖維纏繞層的受力較為復雜,容易發生纖維斷裂或層間剪切的破壞現象。以上原因最終導致了網格理論對靠近赤道位置和臨近金屬接頭破壞位置的設計和分析都不夠準確。因此,應針對T800碳纖維材料的殼體封頭使用試驗或仿真方法,以達到更為理想的研究效果,并進一步開展優化設計研究。

在以往的封頭優化設計研究中,已有較多研究人員對封頭結構進行了仿真分析。崔向斌等[4]使用ANSYS軟件,建立了帶金屬接頭的纖維纏繞殼體模型,并利用殼體解剖數據實現變截面厚度。計算發現在內壓作用下,后封頭在距離赤道附近的區域會因為變形協調產生收縮趨勢,表現出4個屈曲波形態的內凹變形并判斷為殼體封頭薄弱區。程勇等[5]使用ANSYS軟件,將數值應力分析與復合材料成型工藝相結合,采用三維層合實體單元對殼體模型進行有限元計算,發現由于金屬接頭剛度遠大于復合材料殼體剛度,殼體從金屬接頭邊緣附近的應力、應變水平較高;同時由于赤道前后剛度變化,在封頭靠近赤道線處應力、應變水平也較高。PARK等[6]使用非線性有限元分析方法結合ANSYS的APDL語言對殼體的爆破位置進行了預測分析。緱林虎等[7]對橢球封頭模型進行了非線性有限元分析,發現應力應變在金屬接頭附近及接近赤道圓位置變化較大。封頭處的沿纖維方向應變是表現封頭結構可能發生破壞的重要表征指標,特別是封頭靠近赤道位置與金屬接頭位置可能會發生纖維應變的大范圍變化,實際中該位置會因此發生纏繞層的纖維斷裂。由此可知,依據封頭靠近赤道位置或金屬接頭附近的最大沿纖維方向應變得到的優化結構較為可靠。

在補強手段方面,王東等[8]分析了纖維縱向纏繞補強、裁剪條形碳布補強和裁剪環形碳布補強三種補強的效果,通過對這三種補強方法完成的殼體進行水壓爆破試驗,發現采用環形碳布補強對增強殼體后封頭強度的效果較高。關云等[9]利用ABAQUS軟件建立了包括金屬接頭、封頭補強層、彈性層、纏繞殼體層的殼體模型,通過Python對ABAQUS進行二次開發,構建了精細化的復合材料殼體有限元模型,計算結果表明,采用環向補強殼體的封頭內壓承載性能和穩定性能明顯優于縱向補強方法。

目前殼體封頭補強技術主要采用軸向、環向、編織布和封頭帽補強方法[10-11],我國對于傳統的補強手段應用也比較成熟,但對于封頭補強工作研究卻基本依托于經驗,缺乏基礎理論研究作為參考。同時,許多精細化仿真工作僅針對小規模的殼體結構展開,對工程應用缺乏有效的指導作用,基于理論的研究很難與工程實際很好地結合起來。實驗表明,在目前工程中把以往的封頭補強經驗應用于國產碳纖維時,往往沒有辦法達到預期的效果,由于在實際工程中,國產碳纖維的高模量、高強度性能無法很好地體現出來,甚至有時還出現了性能降低的現象,因此有必要對補強結構進行系統化的優化設計,豐富對國產碳纖維纏繞層參數對封頭結構性能影響的認識。而近幾年的文獻顯示,國外對于殼體封頭的結構優化已經傾向于在設計階段就對完整殼體封頭進行優化設計[12-13],這是如今國內研究所欠缺的。因此,有必要系統地基于工程實際和計算基礎,對殼體進行精細化建模,然后進行合理優化設計,以形成一套系統的針對碳纖維纏繞殼體封頭優化設計體系。

本文以某T800碳纖維纏繞殼體后封頭為研究對象,在設計階段引入殼體封頭的優化設計。首先通過試驗及仿真手段確定優化目標,然后基于輔助代理模型設計一套后封頭優化設計方法,最后對優化結果進行評估,形成一套完整的殼體設計方法,并論述未來進一步的研究方向,為今后的大型碳纖維纏繞殼體封頭設計工作提供參考。

1 復合材料殼體模型建立與驗證

1.1 復合材料殼體水壓爆破試驗

根據固體火箭發動機工程需求,本文采用某T800碳纖維纏繞殼體后封頭模型,模型由金屬接頭、橡膠墊片組成。封頭型面為橢球比1.7的橢圓封頭,后極孔半徑為390 mm,半筒身長度為790 mm。殼體采用T800環氧樹脂紗帶纏繞成型,紗帶寬度為18 mm,單層厚度為0.36 mm,筒身纏繞角α=18°,由于殼體為非等極孔結構,故采用非測地線纏繞成型,殼體的纏繞角成型公式為[14-15]

(1)

式中R為縱向纏繞層到中心軸線的徑向距離;R0為過渡點到中心軸線的徑向距離;Rtl為封頭和筒身過度段半徑;θtl為筒身段縱向纏繞層纏繞角;δ為測地線和非測地線纏繞角之間的偏差程度。

當δ=0時,式(1)為測地線纏繞角計算式。

非測地線纏繞成型下的縱向纏繞層厚度模型公式如下[14-15]:

(2)

式中ttl為筒身段縱向纏繞層厚度;θr為對應不同平行圓半徑的纏繞角;rtl為赤道圓半徑;r0為縱向纏繞層轉向的半徑;r為封頭上任意一點對應的平行圓的半徑;BW為紗帶寬度。

對復合材料殼體進行纏繞成型和水壓爆破后結果如圖1和表1所示?？梢钥闯?在16 MPa的內壓載荷下,爆破在封頭靠近赤道部位發生,筒身段的環向應變在0.011 8左右。

圖1 未補強殼體水壓爆破后封頭形貌Fig.1 Hydrostatic burst tests of composite cases without reinforcement

表1 未補強殼體水壓爆破試驗結果Table 1 Results of hydrostatic burst tests of composite cases without reinforcement

1.2 復合材料殼體有限元模型建立

依托ABAQUS仿真軟件建立該殼體后封頭模型,模型包含金屬接頭、堵蓋、橡膠墊片以及T800碳纖維纏繞層四部分。金屬接頭材料使用鋁金屬,墊片材料采用橡膠,堵蓋采用鋼材料,金屬接頭、墊片和堵蓋的材料參數如表2所示,復合材料參數如表3所示。

表2 接頭、墊片及堵蓋材料屬性Table 2 Material properties of joint,gasket and cover

表3 T800碳纖維/環氧樹脂復合材料性能參數Table 3 Performance parameters of T800 carbon fiber/ epoxy resin composites

根據上述材料,使用Python對ABAQUS進行二次開發建立精細化殼體模型,如圖2所示,并開展數值仿真,有限模型由堵蓋、金屬接頭、橡膠墊片和纖維纏繞層組成。為實現封頭有限元模型連續變厚度變角度,本文在ABAQUS中通過細化單元,獲得足夠的單元數量,通過編寫Python 腳本,將纏繞角以1°的角度增量進行離散,生成正負角度,并將其賦予封頭不同位置的單元,這樣做既可實現角度連續,又能達到封頭正負鋪層的效果,縱向纏繞層鋪層角度變化如圖3所示。纏繞層分為11層縱向層和20層環向層,各纏繞層的具體工藝參數如表4所示。

圖2 碳纖維纏繞殼體后封頭模型Fig.2 Aft-dome of carbon fiber wound case

圖3 纏繞角變化曲線Fig.3 Winding angle curves

表4 纏繞層工藝參數Table 4 Winding layer process parameters mm

為減少計算量,采用殼體全模型的1/36仿真分析,堵蓋、金屬接頭和橡膠墊片均采用混雜單元,其中堵蓋包含C3D8R縮減積分單元2360個,C3D6實體單元65個;橡膠墊片中C3D8RH單元60個,C3D6單元6個;金屬接頭含C3D8R單元1460個,C3D6單元32個。復合材料層采用C3D8實體單元,個數為232 992個。

殼體載荷約束根據水壓爆破試驗添加,對殼體模型內面施加15 MPa的壓力,由于采用1/60模型,模型截面施加循環對稱約束,并在底面施加對稱約束(圖4(a))。金屬接頭和堵蓋之間、金屬接頭與墊片、墊片與纏繞層間施加綁定約束(圖4(b)),金屬接頭側面與堵蓋側面、接頭與纏繞層頭部之間采用滑動摩擦約束,摩擦系數設為0.1(圖4(c)、(d))。

1.3 有限元模型計算結果分析及驗證

對于未補強優化的殼體模型的計算變形結果如圖5所示。由仿真結果可知,在15 MPa的內壓環境下,封頭處會產生較大的位移變化,殼體的變形主要集中在最外層纏繞層的臨近接頭、墊片與封頭靠近赤道位置附近,如圖6、圖7所示,這與之前的分析結果一致。

圖5 后封頭位移云圖Fig.5 Displacement contour of the aft dome

圖6 金屬接頭附近應變云圖Fig.6 Strain contour of metal joint of outermost winding layer

圖7 赤道位置附近應變云圖Fig.7 Strain contour of equatorial position of outermost winding layer

提取筒身段的環向應變如圖8所示,可見在筒身內壁會達到最大的環向應變,最大應變為0.011 1,與實驗結果差5.9%,且由于仿真設計爆破壓強為15 MPa,因此結果在誤差范圍之中。由此可知,仿真結果是精確可信的。

圖8 筒身段環向應變Fig.8 Toroidal strain curve of barrel

提取各縱向纖維層的剪切應力變化情況,如圖9所示?？梢钥闯?最外層纏繞層在封頭靠近赤道位置和接頭附近可以達到最大的700 MPa和350 MPa,其中又以封頭靠近赤道位置的應力變化情況最為復雜,數值波動也最為劇烈,由此可知,封頭靠近赤道位置很可能發生剪切破壞,有必要對赤道位置進行進一步的分析。

圖9 縱向纏繞層剪切應力Fig.9 Shear stress of longitudinal winding layer

接著提取各縱向纖維層的沿纖維方向應變云圖后可以發現,在封頭赤道附近存在最大的沿纖維方向應變,如圖10所示,且最大應變會第1或第11層縱向纖維層出現,根據應變曲線圖(如圖11、圖12所示)可以看出,排除纏繞層頭部的計算偏差因素后,最大沿纖維方向應變會在封頭靠近赤道處達到0.016 5,而筒身處的纖維方向模量卻只有0.01,封頭靠近赤道位置處的纖維應變相比筒身處增加了65%。在實際工程中,一般情況下殼體發生爆破時,封頭靠近赤道位置的最大沿纖維方向應變在0.013左右。因此,在15 MPa的內壓環境下,該位置極有可能發生纖維斷裂破壞,此處將作為本文接下來的主要研究位置進行研究。

圖10 縱向纏繞層赤道附近纖維應變云圖Fig.10 Fiber strain near equator of helical winding layer

圖11 第1、11層縱向纖維應變Fig.11 Fiber strain in helical layers 1 and 11

圖12 封頭赤道位置纖維應變Fig.12 Fiber strain at equatorial position of dome

根據分析結果可以得出結論,該殼體結構在封頭靠近赤道位置存在較為復雜的變形情況,如圖13所示,因其處于封頭末端,紗帶間堆疊較少,因此纏繞層厚度較薄,同時該位置處在筒身段環向纏繞層起始位置附近,存在厚度突變現象,此位置的受力情況為多種應力耦合的結果,存在復雜的應力應變變化,且數值波動范圍較大。因此,該位置應成為最重點的研究與優化設計位置。各纏繞層中又以最外層縱向纏繞層的應力應變情況最為復雜、多變,因此也最容易首先發生破壞。

圖13 封頭近赤道位置變形情況Fig.13 Deformation near equatorial position of dome

根據祖磊等[16]對殼體封頭結構的損傷失效分析可知,縱向纖維斷裂失效為殼體爆破失效的主要以及最早發生的失效模式,也是封頭靠近赤道位置的主要失效形式。因此,本文的主要設計優化目標為最外纏繞層封頭靠近赤道位置的最大纖維應變。

2 結果與討論

結構優化設計是從多個設計方案選擇最優的設計方法,以數學理論為基礎,根據結構設計需要滿足的性能指標,設定約束條件,選擇設計優化變量,建立目標函數,尋找最優解。由于現代大多數的工程實際優化問題都比較復雜,一般仿真軟件對單個模型的計算時間會達到幾分鐘甚至幾天,傳統的仿真優化設計使用有限元軟件對大量樣本進行受內壓變形模擬,再聯合優化算法對結構進行優化設計的方法明顯不可行。

為了降低計算時間,減輕計算負擔,本文采用基于近似模型的優化方法,即在優化前進行實驗設計并構建代理模型[17]。本文的優化設計過程主要分為如下幾個步驟:

(1)確定優化目標和設計參數。本文以最大沿纖維方向應變為優化目標,根據不同纏繞層退移半徑、帶寬、補強位置為設計參數建立如下優化模型:

(3)

式中F(Xi)為最大纖維應力目標函數;S、t、Rb為補強優化參數設計變量,分別為退移半徑、帶寬、補強起始半徑。

(2)初始化種群數量,同時根據不同參數組合為單個種群的染色體。

(3)建立有限元模型樣本,求解有限元模型,輸出封頭靠近赤道位置最大纖維應變。

(4)補強參數尋優,通過基因選擇、雜交、變異,求解適應度函數,并判斷最優解是否滿足目標要求,若不滿足則再次循環新一代,返回基因選擇步驟,同時修改有限元模型的參數變量值;若最優解滿足目標要求,則輸出最優的補強角度和補強層數。

根據國產T800碳纖維殼體的纏繞方式、補強鋪層位置等變化影響殼體力學性能,將設計變量定為退移半徑、帶寬、補強起始半徑三個變量,補強材料選取T800碳纖維補強層,共兩層補強帶前后半徑200 mm,每層補強厚度為0.36 mm,分別設立在第4、第5層縱向纏繞層與第5、第6纏繞層之間,即位于縱向纏繞層中部。各參數取值范圍如表5所示。其中,纖維的退移、紗帶寬度和補強層的起始半徑都是設計階段的重要可設計參數。纏繞層退移會引起封頭型面纏繞角和幾何形狀的變化,帶寬涉及纏繞工藝的變化影響幾何形狀發生改變,補強層的位置會直接影響封頭模量和厚度。以上三種參數都可以對以最大纖維應變為目標函數的優化設計起到至關重要的作用。

表5 優化變量參數取值范圍Table 5 Optimization variable parameter value range

實驗設計采用空間填充實驗設計,采樣方法為Latin Hypercube Sampling(LHS)采樣方法,這樣做的好處是可以以較高的離散度在整個設計空間中進行隨機采樣,由于LHS采樣方法屬于分層采樣,可以使設計空間中的所有區域都會被采樣到,確保采樣點對設計空間的全面覆蓋。而且從方法產生的樣本也非常適用于神經網絡的擬合計算[18]。

使用LHS采樣方法采樣50個采樣點后依據神經網絡算法構建本文的代理模型,相較于其他模型,神經網絡對于復雜非線性問題具有更強的逼近能力[19,20]。構建的神經網絡系統結構為3-3-1,即輸入層3個神經元對應于系統的3個輸入變量,隱藏層設3個神經元,僅有的單個輸出對應于輸出層的一個神經元。訓練結果如圖14、圖15所示?？梢钥闯?網絡輸出與訓練集目標的相關系數為0.999(最佳值為1),網絡綜合性能分析相關系數為0.969。由此可見,該訓練模型較好的實現了對訓練樣本集的擬合。

圖14 網絡輸出-訓練樣本線性回歸分析Fig.14 Network output-training sample linear regression analysis

圖15 綜合網絡性能曲線Fig.15 Integrated network performance curves

在建立代理模型之后,可對三個優化變量進行靈敏度分析。圖16～圖18給出了纖維退移距離、紗帶寬度和補強起始半徑對最大沿纖維方向應變的響應曲線。

圖16 退移層退移距離對目標函數的影響Fig.16 Influence curve of backsliding distance of backsliding layer on objective function

圖17 紗帶寬度對目標函數的影響Fig.17 Influence curve of yarn bandwidth on objective function

圖18 補強層起始半徑對目標函數的影響Fig.18 Influence curve of initial radius of reinforcement layer on objective function

相對于原結構,最大纖維應變隨退移距離的變大而呈線性下降,下降范圍在0.017 2～0.015 9,即退移縱向纖維層的退移量越大,封頭靠近赤道處的纖維應變就越小。由數據范圍大幅度變化可見,退移距離參數對最大纖維應變的影響是不可忽視的。最大纖維應變隨紗帶寬度和補強層起始半徑增長呈現非線性的先下降后上升趨勢,存在最佳值,但在其余參數沒有變化的情況下,與原未補強結構相比,對最大纖維應變的影響不是特別明顯,但是過大或過小的紗帶寬度、過大的補強起始半徑位置對封頭靠近赤道位置沿纖維方向應變會有不利影響。綜上所述,對封頭靠近赤道處的沿纖維方向應變影響最大且最直接的變量是封頭近赤道位置的型面幾何形狀復雜程度,同時可知,封頭靠近赤道處應有較為復雜的變形和受力情況,有必要對該位置進行優化研究。

應用遺傳算法對所構建的神經網絡代理模型進行優化搜索,獲取目標函數的極小值點。算法基于Matlab軟件平臺,流程包含選擇、雜交和特征突變三個算法[21],如圖19所示。

圖19 遺傳算法流程圖Fig.19 Genetic algorithm flow char

算法最終給出具有最佳適應性個體的各項參數及適應值,如表6所示,可以看出退移量及紗帶寬度均已達到最大設計變量,補強層則處于后封頭中間位置。由于退移量為影響優化目標的主要影響參數,因此其最優取值與靈敏度分析中所得到的結論相同,接近最大取值。而由于退移量的增大,影響了后封頭型面的纏繞角變化規律以及封頭厚度,采取靈敏度分析中紗帶寬度和補強層位置取值可能會導致封頭型面變形產生不利變化,再加上這兩個變量對優化目標的影響要小于退移層退移距離的影響,考慮到三個變量的聯合作用既相互影響,因此最終優化取值與靈敏度分析得出的最佳取值不同。

表6 最佳適應值參數Table 6 Optimum fitness parameters

由于優化算法的函數值基于代理模型得出,這導致所求得的目標函數值只是近似解,因此有必要再次對優化結果所得參數進行仿真計算,以獲取優化設計結構所對應的目標函數真值,從而對優化設計結果進行評估。仿真后第1、21層縱向纏繞層的沿纖維方向應變曲線如圖20所示?？芍?仿真后得到的最大沿纖維方向應變仿真值為0.015 8,與代理模型預測值之間的誤差為0.63%。這說明構建的代理模型的擬合精度較高,結果的可靠性基本可以保證。與未優化結構進行對比可知,優化結構的最大纖維應變減小了4.2%。

圖20 優化模型封頭靠近赤道處纖維應變Fig.20 Fiber strain at the equator of the optimized model dome

3 結論

(1)根據所建的精細化仿真模型,最大沿纖維方向應變出現于最外層縱向纏繞層的封頭靠近赤道處。此外,靠近金屬接頭處特別是與墊片接近處的最內、外層纏繞層也會存在較大的纖維應變變化與剪切應變。對比水壓爆破試驗可知,在最外層的近金屬接頭部位、封頭靠近赤道部位很容易發生纖維斷裂破壞,在最內層的靠近金屬接頭部位出現剪切破壞的可能性較高。

(2)通過構建代理模型簡化了大量的仿真計算時間,LHS的實驗設計方法做到了對材料空間的全面覆蓋效果,擬合后的相關系數可以達到0.99,目標函數預測值與仿真值之間的誤差為0.63%。優化變量中,退移層退移距離對目標函數影響顯著且呈線性關系;紗帶寬度和補強層設置位置對目標函數影響存在最佳值,不當取值會對目標函數起到惡劣影響。

(3)基于遺傳算法對結構進行了優化設計,獲得了以緩解殼體封頭最大纖維應變為導向的最優殼體封頭設計方案,其中砂帶寬度為25.96 mm,退移層退移距離為19.98 mm,補強層的起始半徑為579.38 mm處,優化后殼體封頭的最大沿纖維方向應變為0.015 8,對比原結構減小了4.2%。

(4)對于其他模型,封頭于接頭肩部位置處的最大順纖維應變可考慮與最大纖維應變同時作為優化目標進行多目標優化設計,這樣可以兼顧封頭的其他重點薄弱位置,同時可以考慮以封頭最大應力作為優化對象進行優化設計。該優化設計與損傷失效分析可作為后續的研究內容的指導方向。