基于ABAQUS的固體火箭發(fā)動機復合材料殼體快速化建模方法及驗證分析①

祖 磊，葛 慶，李德寶，張桂明，吳喬國，張 騫

(1.合肥工業(yè)大學 機械工程學院，合肥 230009；2.合肥工業(yè)大學 土木與水利工程學院，合肥 230009)

0 引言

碳纖維增強樹脂基復合材料(CFRP)具有高比模量、高比強度、優(yōu)異的耐疲勞耐腐蝕性等諸多優(yōu)異特性。其中，纖維纏繞壓力容器作為重要的復合材料能源貯存裝置，在固體火箭、飛行器、海洋船舶、汽車等各個領域有著廣泛的應用。國內(nèi)外針對纖維纏繞壓力容器研究已有較完備的分析設計流程，包括線型設計、失效分析、疲勞分析和爆壓預測等。

近年來，海內(nèi)外大量學者針對復合材料壓力容器進行了大量的研究，JUAN通過纖維纏繞壓力容器插件WCM對復合材料壓力容器的纏繞層進行了精確建模，并對設計壓力為70 MPa的高壓儲氫Ⅳ型氣瓶進行有限元仿真和爆破實驗驗證。JOHANSEN等設計了纖維纏繞分析程序，通過CAD/CAE/CAM一體化設計方法，實現(xiàn)了對任意的軸對稱回轉體及其組合體進行纏繞分析。矯維成基于封頭筒身的纏繞紗帶總體積相等的原理，提出了一種基于三次樣條函數(shù)預測復合材料壓力容器封頭厚度的方法。許家忠等通過多項式逼近算法實現(xiàn)壓力容器纏繞層封頭的厚度分布預測，解決了傳統(tǒng)厚度預測公式在封頭極孔壁厚預測值和實際值存在較大誤差的問題。ZU等通過微分理論推導出了壓力容器的非測地線方程，研究極孔切點數(shù)和紗線寬度對非測地線纏繞方式的影響。王華畢等基于網(wǎng)格理論設計纏繞層，通過三維掃描儀測量壓力容器輪廓驗證厚度預測，然后進行爆壓預測和失效分析。

目前，針對復合材料壓力容器的鋪層設計方法主要包括網(wǎng)格理論和有限元方法。其中，建立準確的壓力容器模型是有限元分析準確的前提。由于碳纖維復合材料壓力容器在封頭處的纏繞層厚度為連續(xù)變厚度、變角度。導致壓力容器有限元模型建模難度大、建模工作繁瑣。針對纏繞層數(shù)較多的模型，考慮時間成本只能對模型簡化處理，影響最終的有限元計算精度。

國內(nèi)外學者在纖維纏繞壓力容器的設計及建模分析等方面進行了大量的研究，但針對殼體的快速化建模方法研究較少。綜上，本文在前人研究的基礎上，提出了一種固體火箭發(fā)動機復合材料殼體快速建模的方法。實現(xiàn)對商業(yè)纏繞軟件信息文件進行數(shù)據(jù)處理，導入到ABAQUS中參數(shù)化生成有限元模型的功能。通過復合材料殼體的水壓爆破實驗驗證，與生成的有限元模型計算結果進行對比，驗證了建模方法的準確性。這種復合材料殼體快速化建模方法提高了建模效率，能較快速地實現(xiàn)從纏繞設計軟件導入生成有限元幾何模型，同時可以便捷地生成纏繞層封頭補強，縮短了設計周期，該研究成果對殼體的設計分析有著重要的意義。

1 復合材料殼體快速建模方法

1.1 快速化建模流程

復合材料殼體有限元模型的準確性主要是由封頭上纖維厚度和纏繞線型所決定的。傳統(tǒng)的殼體設計及分析流程中，需要手動在ABAQUS中建立有限元模型。因殼體在封頭處的角度和厚度是連續(xù)變化的，導致實際建模的過程中封頭處的建模難度大，需要花費大量時間。考慮時間成本問題，人工對封頭進行切分賦予角度時，無法保證鋪層角度的連續(xù)性，影響模型計算精度。

本文基于ABAQUS二次開發(fā)功能，實現(xiàn)了從纏繞軟件中快速導入的參數(shù)化建模功能，并可以實現(xiàn)包括纏繞補強、封頭補強、螺旋纏繞等多種鋪層方式的建模。可以將幾天的建模工作縮短至幾分鐘，可以較大程度地縮短設計人員對有限元模型的前處理時間，建模時間成本降低，方便殼體的結構設計及優(yōu)化分析。本插件除了實現(xiàn)標準橢圓封頭結構的殼體建模功能，也可以實現(xiàn)如固體火箭發(fā)動機尾噴管等異形件的非標準結構建模功能。

纏繞軟件設計線型建立模型后，可生成纏繞層信息文件，其中包括纏繞層落紗點坐標、落紗點纏繞角度和落紗點位置紗線厚度。圖1為CADWIND軟件的纏繞層信息文件部分截圖。本文將纖維纏繞層數(shù)據(jù)離散化，通過PYTHON語言對纏繞軟件生成的信息文件進行數(shù)據(jù)處理，在ABAQUS中參數(shù)化建模生成有限元模型，實現(xiàn)復合材料殼體纏繞層的快速化建模。

圖1 數(shù)據(jù)文件Fig.1 Data file

基于纏繞層信息文件和PYTHON語言，編寫了數(shù)據(jù)處理程序。信息文件中的數(shù)據(jù)一般根據(jù)纖維實際纏繞的位置生成對應參數(shù)，纏繞的坐標變化順序是跟著絲嘴軌跡而變化的，故信息文件中纏繞坐標不是線性變化的。數(shù)據(jù)處理程序先讀取數(shù)據(jù)文件中的落紗點坐標，根據(jù)坐標相對位置進行排序，剔除重復數(shù)據(jù)，繪制出纏繞層內(nèi)輪廓。然后，讀取各落紗點位置上的紗線厚度，繪制出纏繞層外輪廓，依次往復完成殼體有限元模型建立。

部分纏繞軟件如CADWIND基于面片法原理，會根據(jù)芯模輪廓打斷成直線進行求解建模，曲線打斷的直線疏密會影響纏繞層信息文件，采樣的數(shù)據(jù)往往會產(chǎn)生噪點，甚至會偏離模型輪廓，如圖2所示。導致生成的有限元模型表面會產(chǎn)生褶皺，不利于模型的分析。這里采取B樣條函數(shù)法對纏繞層外輪廓數(shù)據(jù)進行光滑處理，保證生成的模型表面光滑連續(xù)。

圖2 原始數(shù)據(jù)Fig.2 Original data

由于部分商用纏繞軟件只能生成螺旋層纏繞數(shù)據(jù)信息。為滿足更通用性的建模工作，本方法中增加了環(huán)向層及補強層建模功能。在纏繞層內(nèi)部中出現(xiàn)了環(huán)向層或補強層時，可以根據(jù)相鄰螺旋纏繞層位置判斷生成環(huán)向層和補強層數(shù)據(jù)。如果纏繞層最外層為環(huán)向層或補強層時，信息文件只能獲取螺旋層位置，無法求解環(huán)向層或補強層位置，需要在插件界面輸入最外層鋪層信息，輔助生成對應鋪層的有限元模型。

基于ABAQUS二次開發(fā)方法，編寫了纏繞軟件導入的快速化殼體建模腳本程序。將數(shù)據(jù)處理程序中采集到的數(shù)據(jù)輪廓點導入，建立幾何模型。筒身段，根據(jù)纏繞角進行材料屬性的賦予。封頭段，通過切片法對封頭模型沿母線方向進行切分，將數(shù)據(jù)處理程序中得到的纏繞角信息依次對封頭進行材料屬性的賦予，通過柱坐標系離散化的方法完成封頭模型的材料屬性賦予。然后，根據(jù)需求對殼體模型施加對應的邊界條件、工作載荷以及最終的網(wǎng)格劃分，完成整個殼體模型的前處理工作。

1.2 纏繞軟件導入數(shù)據(jù)的處理方法

為解決數(shù)據(jù)處理程序中波動的外輪廓問題，本文采用了數(shù)值分析中的B樣條方法，對實際的波動數(shù)據(jù)進行處理。樣條線的概念是過一組指定點集而生成的柔性光滑曲線，而B 樣條曲線法的原理即通過控制點集中局部位置，從而達到控制一條曲線形狀的樣條函數(shù)。

B樣條定義方法較多，一般采取截尾冪函數(shù)差商定義，這里采用德布爾-考克斯遞推公式。B樣條曲線的方程為

(1)

其中，為控制頂點，，()為基函數(shù)，而基函數(shù)表示如下：

，()的第一個下標表示序號，第二個下標表示階次。如確定第個次B樣條，()，需要確定，+1，…，++1一共+2個節(jié)點。這里把[，++1]定義為，()的支撐區(qū)間。

B樣條光滑曲線原理如圖3所示。

圖3 B樣條數(shù)據(jù)處理Fig.3 B-spline data processing

當兩節(jié)點為同一位置時，遞推公式中組合系數(shù)在特殊情況時會出現(xiàn)分子分母均為0的現(xiàn)象，這種情況直接令系數(shù)為0即可。

該曲線在區(qū)間[，++1]中某一點()是由+1個頂點(=-，-+1，…，)和基函數(shù)決定，與其他頂點無關，如圖4所示。

圖4 B樣條頂點控制原理Fig.4 B-spline vertices approach

由圖4可見，移動第個控制頂點，最多只會對區(qū)間(，++1)上的局部曲線產(chǎn)生影響，而其他區(qū)間的B樣條曲線不會受到干預。因此，階次越高，B樣條函數(shù)會擬合出的數(shù)據(jù)點越光滑，但遞推函數(shù)計算量就越大。

1.3 可視化界面

基于PYTHON語言和ABAQUS二次開發(fā)功能，為建模方法設計了可視化界面。可以在插件中界面輸入相關信息進行模型的前處理工作，方便設計人員操作，界面主要分為兩部分，如圖5所示。

(a)Interface A (b)Interface B圖5 插件界面Fig.5 Plug-in interface

針對任意形狀的構件，可以在圖5(a)中輸入?yún)?shù)實現(xiàn)模型生成。在文件指派模塊中，選擇纏繞軟件生成的纏繞層信息文件，然后在纏繞層參數(shù)指派模塊中，輸入設計爆壓、環(huán)向單層厚度、三維模型旋轉角度以及筒身長度。在設計參數(shù)指派模塊中，輸入工程常數(shù)和纏繞層形式，其中主要纏繞形式為螺旋纏繞、環(huán)向纏繞等。在數(shù)據(jù)處理模塊中，可以將每次使用到的參數(shù)進行保存。再次使用時，相同的參數(shù)選中文件可以直接導入，節(jié)約輸入?yún)?shù)的時間。

針對任意比例的橢圓殼體，可在圖5(b)中實現(xiàn)參數(shù)化建模。首先，在ABAQUS中建立芯模模型；之后，在參數(shù)欄中輸入相關纏繞和建模參數(shù)，實現(xiàn)包括螺旋纏繞、環(huán)向纏繞、擴孔和補強等多種纏繞形式的殼體有限元模型建立。

2 殼體纏繞有限元模型建立

文中插件實現(xiàn)的功能是對已設計完成的殼體結構進行快速化建模功能，側重體現(xiàn)在建模功能，線型設計部分并非此插件功能，故直接取網(wǎng)格理論設計好的殼體鋪層方案，使用合肥工業(yè)大學研發(fā)的纏繞軟件生成對應模型的纏繞層信息文件。使用本插件將信息文件導入到ABAQUS中生成殼體的有限元模型。根據(jù)有限元軟件中仿真結果與對殼體的水壓爆破實驗進行對比，觀察殼體在不同壓力階段的應變、位移和失效形式等。

2.1 復合材料殼體結構

本文使用的殼體芯模主要由橡膠內(nèi)襯、金屬接頭和砂芯組成。芯模的具體尺寸如圖6所示，直徑為150 mm的殼體，左右極孔尺寸分別為30 mm和50 mm。

圖6 芯模尺寸Fig.6 Mandrel dimension

通過網(wǎng)格理論對殼體鋪層進行設計，設計爆壓為30 MPa。設計考慮封頭承壓能力弱，適當增加纏繞層，減少環(huán)向層。由于本插件可以實現(xiàn)補強建模功能，故增加了補強層輔以驗證。其中，纏繞補強是通過將筒身段的復合材料層剪斷，只保留封頭段復合材料層作為補強層，在有限元模型中表現(xiàn)形式如圖7所示，設計的纏繞方案如表1所示。

圖7 補強區(qū)域模型圖Fig.7 Model diagram of the reinforce part

表1 鋪層方案Table 1 Layup scheme

2.2 材料參數(shù)

采用干法纏繞進行實驗，選用T700-12k纖維預浸帶，中溫固化為120 ℃，預浸帶厚度0.15 mm，樹脂含量32%，纖維含量150 g，帶寬6 mm。

材料參數(shù)見表2。

表2 T700-12k預浸帶力學性能參數(shù)Table 2 Mechanical properties of T700-12k prepreg

2.3 有限元模型導入

在纏繞軟件中輸入芯模尺寸參數(shù)和對應的鋪層方案，生成纏繞層模型及導出的纏繞層封頭角度如圖8所示。

圖8 復合材料殼體線型仿真圖Fig.8 Composite case line type simulation

快速化建模功能是導入纏繞軟件生成的數(shù)據(jù)文件進行數(shù)據(jù)處理生成模型。基于本課題組自研的纏繞軟件導入生成的模型，自研纏繞軟件中纖維纏繞層厚度分布采取的三次樣條函數(shù)法進行預測，矯維成等對三次樣條函數(shù)法預測厚度進行了理論推導和實驗驗證，證明了三次樣條函數(shù)厚度預測的準確性，在纖維不滑紗的情況下，實際纖維最大厚度在一個帶寬附近。

為節(jié)約計算時間，建立復合材料殼體的1/72模型。將纏繞層信息文件導入到插件中，殼體纏繞層采用C3D8R單元，生成的殼體模型已經(jīng)完成了材料屬性賦予和網(wǎng)格劃分工作，完成了殼體的前處理工作。本建模流程一般只需幾分鐘即可完成，實現(xiàn)了殼體的快速化建模功能。

因纏繞軟件中只能生成纏繞層信息文件，不包含金屬接頭部分。在ABAQUS中建立金屬接頭模型，將殼體與金屬接頭裝配。對殼體中間部分施加軸向位移約束，避免有限元計算中發(fā)生剛體位移，殼體和接頭模型側邊施加循環(huán)對稱條件。在殼體內(nèi)表面施加均布壓力，從0～30 MPa依次增加。

邊界及載荷施加如圖9所示。

圖9 載荷與邊界條件Fig.9 Load and boundary conditions

2.4 殼體損傷模型

漸進損傷模型(Progressive Damage Model)已被廣泛應用于復合材料結構的漸進破壞和失效分析中。首先，選取適合的失效準則表征材料失效狀態(tài)，根據(jù)不同失效形式的退化模型，對材料進行剛度折減，當損傷到達一定程度后，整體結構失效，完成計算分析。漸進損傷計算流程如圖10所示。

圖10 漸進損傷流程圖Fig.10 Flow chart of progressive damage process

針對殼體結構失效分析，主要考慮的失效模式有纖維拉伸壓縮破壞、基體拉伸壓縮破壞、纖維基體剪切破壞和分層破壞。

纖維拉伸破壞：

(3)

纖維壓縮破壞：

(4)

基體拉伸破壞：

(5)

基體壓縮破壞：

(6)

纖維基體剪切破壞：

(7)

分層破壞：

(8)

式中 下標11、22、33表示正軸坐標系下1、2、3方向的正應力；、、為剪切力；、分別為表單層板1方向的拉伸、壓縮強度；、分別為單層板2方向的拉伸、壓縮強度；、、分別為單層板各向的剪切強度；、分別為3方向的拉伸和壓縮強度；α為經(jīng)驗系數(shù)，=1。

應力滿足破壞條件時，會對材料損傷部分進行剛度折減。本文采用CAMANHO材料退化模型進行漸進失效分析。材料退化參數(shù)如表3所示。

表3 材料性能退化參數(shù)Table 3 Degradation parameters of material properties

3 結果與討論

3.1 實驗結果分析

為殼體設計水壓爆破實驗，如圖11所示，監(jiān)測爆破中殼體的應變位移變化。在殼體上均勻設置應變監(jiān)測點，具體位置如圖12所示。

圖11 實驗測量示意圖Fig.11 Schematic diagram of experimental measurement

圖12 應變測量點分布Fig.12 Distribution of strain measurement points

由于升壓的過程中，殼體變形量較大，殼體伴隨有樹脂開裂，導致有一些應變片損壞，應變值無法顯示，取數(shù)據(jù)完整的測試點1、3、5、7、9、13、19和21的結果，如圖13(a)所示，其中每個應變測點監(jiān)測兩個方向的應變變化，取纖維方向為a，垂直于纖維方向為b，應變單位為με。在小端極孔處的13號測點應變產(chǎn)生了突變。位移變化監(jiān)測設置兩點，與應變監(jiān)測點3、8重合。提取兩點位移結果如圖13(b)所示。兩點均在27 MPa后產(chǎn)生了位移突變。結合實驗中應變位移測試結果可知，此殼體的實際爆破壓強為27 MPa。

(a)Tendency of strain with pressure (b)Tendency of displacement with pressure圖13 應變與位移曲線Fig.13 Strain and displacement curves

殼體爆破實驗結果如圖14所示。可見，小端極孔接頭的纖維從接頭的肩部斷裂，接頭飛出。因為封頭肩部具有較大的應力集中現(xiàn)象，導致纖維方向上的應力過大。由于封頭肩部軸向承載能力大于環(huán)向承載能力，在均布載荷下，環(huán)向承受應力大于軸向，從而極易導致基體的損傷。隨著環(huán)向應力的增大，基體損傷加劇，而基體損傷也會引起纖維的損傷。直到纖維損傷到達其強度極限后，纖維最終會斷裂，從而導致封頭肩部的纖維全部剪斷，金屬接頭飛出，軸向位移突變。

圖14 實驗結果Fig.14 Experimental result

3.2 有限元結果分析

通過觀察殼體有限元計算結果對比可知，殼體的主要失效原因為纖維失效，失效位置發(fā)生在殼體小端封頭處，這里取殼體有限元計算中纖維拉伸失效結果進行分析，失效部位與未發(fā)生失效部位取高對比度的色階便于觀察。當內(nèi)壓為22.8 MPa時，整個復合材料纏繞層并未發(fā)生失效，如圖15(a)所示。隨著內(nèi)壓增加到23.2 MPa時，小端極孔封頭肩部附近纖維發(fā)生失效，如圖15(b)所示。隨著內(nèi)壓逐漸增加，復合材料層逐層的纖維不斷發(fā)生失效。當內(nèi)壓最終達到26.02 MPa時，金屬接頭肩部附近復合材料層纖維發(fā)生貫穿性失效，最終導致封頭在該位置發(fā)生爆破，如圖15(c)所示。

(a)p=22.8 MPa (b)p=23.2 MPa (c)p=26.02 MPa圖15 纖維破壞云圖Fig.15 Damage contours of fiber

除封頭肩部位置，在筒身赤道處發(fā)生的損傷較多，但沒有形成致命的失效。殼體的軸向位移云圖如圖16所示，提取小端封頭段某節(jié)點的徑向位移和軸向位移，如圖17所示。在內(nèi)壓為26.02 MPa時，殼體小端極孔封頭軸向位移發(fā)生突變，而徑向位移并未發(fā)生突變，表明此時殼體發(fā)生軸向爆破，即外在表現(xiàn)形式為金屬接頭軸向飛出。

圖16 軸向位移云圖Fig.16 Axial displacement contours

圖17 封頭位移曲線Fig.17 Head displacement curves

通過實例分析對比可知，本建模方法下建立的有限元模型可以較準確地對復合材料殼體失效形式進行分析。只需從纏繞軟件中生成線型文件，然后在插件界面中配置相關參數(shù)，即可快速生成殼體的有限元模型。如果計算結果不滿足設計要求，只需更改鋪層參數(shù)等相關信息，或者在纏繞軟件中設計線型，重新計算，直到滿足設計要求即可。在滿足殼體強度剛度的要求下，也可對殼體結構進行優(yōu)化設計，從而達到節(jié)約成本的目的。

4 結論

(1)針對固體火箭發(fā)動機復合材料殼體，提出了一種基于ABAQUS二次開發(fā)的快速化建模方法。通過對纏繞軟件生成的纏繞層信息文件進行數(shù)據(jù)處理，導入到ABAQUS中，實現(xiàn)從纏繞設計軟件到有限元分析軟件的快速化建模方法。

(2)在數(shù)據(jù)處理過程中，針對波動性較大的數(shù)據(jù)，通過B樣條函數(shù)對厚度進行光滑處理，保證生成模型的平滑性。

(3)通過殼體水壓爆破實驗結果的位移和失效形式分析，對比本方法生成的有限元仿真計算結果，對建模方法的準確性進行了驗證。

(4)對上述建模方法設計了可視化界面，方便設計人員使用，大大縮短了設計的建模時間。在實際工程中，固體火箭發(fā)動機復合材料殼體設計與分析有著重要的意義。