基于參數化建模的藥柱傘盤結構形狀優化①

李 磊，段靜波，申志彬，唐國金

(1.總參工程兵科研三所，洛陽 471023;2.國防科技大學 航天與材料工程學院，長沙 410073)

0 引言

藥柱設計是發動機設計的核心部分，藥柱形狀與其體積裝填分數和結構完整性有密切關系，而體積裝填分數和結構完整性是評價藥柱設計優劣的重要指標［1］，所以在藥柱設計過程中，應十分注重藥柱的形狀優化。目前，考慮結構完整性和裝填分數要求的藥形優化還較困難，因為藥柱的幾何形狀往往較復雜，固體推進劑又具有粘彈性特性，準確計算其應力應變狀態具有一定難度，而且在形狀優化過程中分析模型不斷變化，必須不斷地重新生成有限元網格，并進行自適應分析，使得研究難度進一步加大。目前，粘彈性有限元方法和計算軟件的發展，給復雜藥形發動機的結構完整性分析提供了理論與工具支持，計算軟件的二次開發工具又為復雜結構的形狀優化提供了便利條件，利用大型商用有限元軟件的二次開發工具進行各種結構的優化設計成為新興的研究熱點［2－7］。

MSC.Patran軟件是一個并行框架式有限元前后處理及分析仿真系統，可很好地支持 MSC.Nastran和MSC.Marc等有限元分析程序，已廣泛用于航空航天、汽車、造船和國防等各大領域。MSC.Patran具有強大的網格劃分功能，可方便地設定有限元模型的載荷、邊界條件、材料和單元特性，并為用戶提供了便捷的二次開發工具——PCL(Patran Command Language)，使不同領域的用戶可根據自身特點進行專用軟件的二次開發。張國棟［2］利用PCL建立了存在接觸關系的裝配體的參數化模型，并利用遺傳算法進行了結構優化設計。隋允康［4－5］用 PCL 語言對 MSC.Patran/Nastran 軟件進行二次開發，實現了膜結構的截面優化設計，并采用二級控制法對二維連續體進行形狀優化，用PCL語言開發了相應的優化模塊。

錐柱形藥柱中的環向槽通常被稱為傘盤，主要用于調節燃面并改善藥柱的受力狀態。在低溫載荷作用下，傘盤頂部易產生裂紋，是藥形設計關注的重點部位之一。本文主要研究錐柱形藥柱中的傘盤結構，采用PCL語言對有限元軟件MSC.Patran/Marc進行二次開發，建立發動機的參數化有限元模型，同時考慮藥柱的體積裝填分數和結構完整性要求，利用遺傳算法對藥柱傘盤結構進行形狀優化，得到體積裝填分數最大且滿足結構完整性要求的傘盤形狀。

1 發動機模型

1.1 幾何形狀

錐柱形發動機為軸對稱結構，可建立其縱截面模型進行軸對稱分析，以節省計算量。傘盤結構位于發動機頭部，傘盤頂端采用雙圓弧形式，由前后2個90°圓弧組成。發動機縱截面的幾何尺寸有殼體長度Lc、殼體外徑D、殼體頭部大圓半徑R1、殼體頭部過渡圓半徑 R2、殼體厚度 δc、絕熱層厚度 δs、襯層厚度 δd、藥柱前圓柱內徑d、藥柱內徑差Δd、傘盤起始位置Ls、傘盤寬度W、傘盤深度H、傘盤頂端前圓弧半徑r1、傘盤頂端后圓弧半徑r2，其主要幾何尺寸如圖1所示。

圖1 錐柱形發動機縱截面示意圖Fig.1 Sketch of the conocyl SRM axial section

發動機默認的幾何尺寸(單位mm，下同):Lc=4 000，D=1 000，R1=600，R2=200，δc=4，δs=2，δd=1，d=220，Δd=20，Ls=300，W=80，H=400，r1=50，r2=W－r1=30。

1.2 材料參數

發動機的殼體、絕熱層為彈性材料，藥柱、襯層為粘彈性材料。在MSC.Marc軟件中，采用Prony級數模型描述粘彈性材料本構關系。在變溫過程中，可采用WLF方程描述熱粘彈性本構關系。熱粘彈性材料參數主要包括剪切松弛模量的Prony級數、體積松弛模量的Prony級數和WLF方程的方程參數。

在固體推進劑的泊松比ν不隨時間變化的前提下，體積松弛模量K(t)和剪切松弛模量G(t)滿足如下關系式:

故而，只需定義推進劑的剪切松弛模量和泊松比，即可確定其體積松弛模量。

發動機各部分的材料參數:

(1)殼體的彈性模量2.06 ×105MPa、泊松比 0.3、線膨脹系數1.1×10－5;

(2)絕熱層的彈性模量 22.0 MPa、泊松比 0.498、線膨脹系數2.2×10－5;

(3)襯層的泊松比0.498、線膨脹系數8.0 ×10－5;

(4)藥柱的泊松比0.498、線膨脹系數8.6 ×10－5。

襯層的松弛模量與藥柱的相同，其剪切松弛模量表達式為

WLF方程參數為

式中 T為熱力學溫度。

1.3 載荷與邊界條件

采用MSC.Marc軟件分析發動機在溫度載荷作用下的結構完整性。復合固體推進劑的固化溫度為50℃，則其零應力溫度為58℃，假設在1 d(86 400 s)時間內，溫度由58℃線性下降到－40℃，且整個發動機為均勻溫度場。發動機殼體尾部約束軸向位移，以使整個發動機不產生剛體位移，殼體頭部和側壁為自由狀態，可沿軸向和徑向自由收縮或膨脹。

1.4 結果輸出

在低溫載荷作用下，藥柱傘盤頂部產生裂紋的主要原因是應變過大，結構完整性遭到破壞。為了避免裂紋的產生，必須保證藥柱的最大Von Mises應變小于其許用應變，故計算完畢后，將藥柱傘盤部位的最大Von Mises應變作為計算結果進行輸出。

另外，在固體發動機藥柱設計過程中，往往希望其體積裝填分數盡量大，以提高發動機的質量比，藥柱的體積裝填分數也需計算輸出。由于錐柱形藥柱為軸對稱結構，將縱截面模型旋轉一定角度，即可得到實體模型，利用PCL語言求出實體模型的體積，即可得到藥柱的體積裝填分數，將其也作為計算結果進行輸出。

2 發動機參數化建模方法

MSC.Patran的建模過程一般是界面操作，大部分界面操作在其命令歷史(Command History)窗口都會顯示對應的命令執行情況，這些命令就是基本的PCL語句，將界面操作對應的PCL語句直接輸入命令輸入行，可起到與界面操作相同的效果。對于每個新建的有限元模型數據庫文件(model_name.db)，都會自動生成記錄其建模歷史命令的日志文件(model_name.db.jou)，從日志文件中可得到建立該有限元模型用到的所有PCL語句，這些語句是進行參數化建模的基本依據。

發動機參數化建模的基本步驟:

(1)根據固體發動機的結構形式和幾何特點，按照分析人員的已有經驗進行手動建模，并通過不同建模方式的比較，選擇一種適用于類似結構形式的通用建模方法，通用建模方法包括幾何模型生成、有限元網格劃分、載荷和邊界條件添加、材料和單元特性設定及分析參數設置等內容;

(2)利用通用建模方法建立發動機的有限元分析模型，進行分析計算并輸出計算結果，在日志文件中，提取所有PCL語句，形成參數化建模的基本命令流;

(3)選擇需要進行參數化的幾何尺寸和有限元網格控制參數(也可是材料參數、載荷數值等)，將其聲明為變量，把基本命令流中與其相關的數值用變量替代，就得到了建模用的PCL程序，程序中所定義的變量即為參數化模型中的參數。

對參數化建模程序中的變量賦值之后，重復調用建模程序，即可建立不同的有限元分析模型，并自動進行分析計算和結果輸出。可采用其他編程語言(如Matlab、C++等)，在 MSC.Patran環境中調用 PCL程序，用以計算發動機的各種響應值。主程序與PCL程序之間的數據傳遞，可采用形狀參數或數據文件的方式進行，將PCL程序作為主程序函數進行調用。

經過多次建模的對比分析，得到錐柱形藥柱發動機的參數化建模程序。圖2所示為發動機縱截面有限元網格，約有3 500個單元和3 800個節點，具體數量根據幾何參數的不同稍有變化。

3 藥柱傘盤結構形狀優化

3.1 優化問題描述

本文對藥柱傘盤結構進行形狀優化的目的是在滿足結構完整性要求的前提下，設計出裝填分數最大的藥柱，該優化設計問題可描述為

式中 X為設計變量;N為設計變量數目;?(X)為目標函數，即藥柱的體積裝填分數;ε(X)為藥柱傘盤部位的最大Von Mises應變;εmax為藥柱Von Mises許用應變;XL和XU為設計變量的下限和上限。

圖2 錐柱形發動機縱截面有限元網格Fig.2 Finite element mesh of the conocyl SRM axial section

進行藥柱傘盤結構的形狀優化時，不能將所有的幾何尺寸都作為參數進行優化，而應有所選擇。文獻［8］的研究認為，傘盤寬度與傘盤深度對傘盤頂部的應變有較大影響。文獻［9］的研究指出，傘盤頂部的曲面結構形式也會顯著影響傘盤頂部的應變分布。另外，藥柱前圓柱內徑對藥柱的體積裝填分數和應變都有顯著影響，傘盤起始位置影響傘盤在藥柱中的位置。所以，本文選擇傘盤寬度、深度、起始位置、頂端前圓弧半徑和藥柱前圓柱內徑5個幾何尺寸作為形狀優化問題的設計變量。

藥柱的裝填分數與幾何尺寸的關系復雜，其顯示表達式推導繁瑣，且藥柱的最大Von Mises應變與幾何尺寸沒有明顯的顯式函數關系，其導數信息很難獲取。傳統的優化設計方法對目標函數和約束函數的連續性與可微性有較高要求，一般都要進行求導運算，不宜采用。遺傳算法是一種模擬生物進化過程的優化方法，對目標函數和約束函數無連續性與可微性要求，僅使用適應度函數值，就可確定進一步的搜索方向和范圍，對于各種優化問題具有廣泛的適應性，故本文采用遺傳算法進行藥柱形狀優化。

3.2 遺傳算法設計

采用遺傳算法進行優化設計的基本原理:首先隨機生成初始群體，然后基于自然界優勝劣汰、適者生存的思想，選擇適應性強的個體進行交叉與變異，模擬生物的繁殖與突變過程，得到的新群體適應性更強，整個群體便得到了進化，如此反復，經過多代進化，最后便可得到適應性最強的個體，即優化問題的最優解。遺傳算法設計主要包括編碼方式確定、適應度函數與約束函數處理、遺傳算子與收斂準則設計等內容，根據優化問題的不同類型需要，進行有針對性的算法設計。

藥柱形狀優化問題的設計變量為幾何尺寸，其數據類型為實數，如按照實數類型進行優化，一方面計算量太大，另一方面優化出的實數幾何尺寸在進行加工制造時精度很難達到。為了減少計算量，同時考慮到加工工藝等因素，將設計變量取值限定為整數。根據變量類型，遺傳算法的編碼方式采用整數編碼，算法的搜索空間與設計變量的取值空間一一對應。

適應度函數表征個體適應環境的能力，與目標函數和約束條件直接相關。本文研究的優化問題為藥柱的裝填分數最大化問題，約束條件為不等式約束，可將藥柱的裝填分數直接作為適應度函數，并采用罰函數方法對群體中不滿足約束條件的個體進行懲罰，降低其適應度函數值。適應度函數的計算方法如下:

式中 F(X)為適應度函數;P(X)為罰函數，要求滿足0＜P(X)＜1，本文取 P(X)=εmax/ε(X)。

遺傳算法中的遺傳算子主要有3種:選擇算子、交叉算子和變異算子。選擇算子以適應度函數值為依據，適應度越大，被選擇的概率越大。本文采用線性排序選擇算子，按照個體在群體中的適應度排名分配選擇概率，根據排名由前到后選擇概率線性變化，選擇概率計算公式為

式中 M為群體規模;η+為當前群體最佳個體期望選擇數量;η－為當前群體最差個體期望選擇數量，要求η++η－=2，且 η+≥η－≥0。

為了保證算法的收斂性，采用最優保存策略，以使整個優化過程中的最優個體不會丟失。

交叉算子采用隨機算術交叉，當個體X1={x11，x12，…，x1n}和 X2={x21，x22，…，x2n}進行交叉時，產生的2個新個體為

式中 ri(x，y)為隨機取整函數，表示以等概率選取［x，y］區間內的整數。

變異算子同時采用非均勻變異和最優變異。非均勻變異以非均勻變異概率對個體的基因進行操作，對個體X={x1，x2，…，xN}中的基因xk進行變異操作時，如果其取值范圍為［xkL，xkU］，則新基因值由下式確定:

式中 ceil(x)為取整函數，函數值為不小于x的整數;r為［0，1］之間的均勻分布隨機數，rdi為［0，1］之間均勻分布的隨機整數;Tt為遺傳代數;Tz為終止代數;b為控制參數，本文取為2。

最優變異按照最優變異概率對最優個體進行變異，取最優個體的一個隨機基因，然后將基因值隨機進行“+1”或者“－1”操作。如果基因值為其取值下限，則只能進行“+1”操作;如果基因值為其取值上限，則只能進行“－1”操作。變異后的個體替換群體中的一個隨機選取的個體。

遺傳算法的收斂準則用于判斷算法是否收斂于最優解，以決定是否終止計算。采用針對最優個體的收斂準則，即全群體最優個體超過Tc代沒有更新就判斷算法收斂。該準則以最優個體的穩定程度作為判斷收斂的標準，必須與最優保存策略同時使用。同時，設定一個終止代數Tz，當遺傳代數Tt達到終止代數時終止計算。遺傳算法的流程見圖3。

圖3 遺傳算法流程圖Fig.3 Flow chart of the GA

3.3 優化實例

將發動機藥柱的傘盤寬度、深度、起始位置、頂端前圓弧半徑和藥柱前圓柱內徑5個幾何尺寸作為設計變量，其他幾何尺寸、材料參數和載荷參數不變，在藥柱傘盤的最大Von Mises應變小于15%的約束條件下，以藥柱的裝填分數最大為優化目標，對傘盤結構進行形狀優化。傘盤寬度W取值范圍為50～150，傘盤深度H取值范圍為300～450，傘盤起始位置Ls取值范圍為200～400，藥柱前圓柱內徑d取值范圍為100～300，前圓弧半徑r1的取值范圍與傘盤寬度有關，r1與W的比值既不能過大，也不能過小，將其限定為0.2～0.8之間。由于采用整數編碼遺傳算法，將r1與W的比值處理為整數較為方便，故定義λ=r1/W×100，將λ作為設計變量，取值范圍為20～80。

文獻［10］中認為，發動機頭部脫粘深度Hd(從軸線到脫粘層根部距離)對傘盤部位的應變影響很大，但脫粘深度與藥柱形狀沒有關系，故脫粘深度不作為設計變量，但選取2種不同脫粘深度的模型進行優化，分別是Hd=200的模型和Hd=450的模型。

(1)Hd=200的模型

遺傳算法的控制參數為群體規模20，交叉概率0.9，非均勻變異概率 0.1，最優變異概率 0.5，收斂準則中的Tc=10，終止代數Tz=100。

算法運行91代收斂，群體適應度函數最大值與平均值隨遺傳代數的變化曲線見圖4。得到的優化結果為{116，444，20，111，353}，對應的幾何尺寸為傘盤寬度W=116，傘盤深度H=444，前圓弧半徑r1=23.2(λ=20)，藥柱前圓柱內徑d=111，傘盤起始位置Ls=353，傘盤部位最大Von Mises應變為14.97%，體積裝填分數0.968 0。形狀最優藥柱傘盤部位的Von Mises應變分布如圖5(a)所示。

從圖4看出，群體適應度函數最大值穩定增長，由于遺傳算法中采用了最優保存策略，保證了最優個體的穩定性，在其保持10代沒有更新的情況下，達到了算法的收斂條件，得到優化結果。圖5(a)中，藥柱傘盤頂部為應變集中部位，這是結構本身的受力特點所決定的，傘盤頂端的應變分布相對均勻，說明頂端曲面的形狀是合理的。

(2)Hd=450的模型

遺傳算法的控制參數與算例(1)完全相同，算法運行90代收斂，得到的最優藥形幾何尺寸為傘盤寬度W=50，傘盤深度 H=300，前圓弧半徑 r1=10.5(λ =21)，藥柱前圓柱內徑d=100，傘盤起始位置Ls=302，傘盤部位最大Von Mises應變為14.97%，體積裝填分數0.980 7。形狀最優藥柱傘盤部位的Von Mises應變分布如圖5(b)所示。圖中所示藥柱頭部的最大Von Mises應變為24.87%，位于脫粘層根部，而傘盤頂部的最大Von Mises應變只有14.97%，脫粘層材料的延伸率一般比推進劑高得多，且本文主要研究傘盤結構。所以，只以傘盤頂部的最大Von Mises應變為約束條件，不考慮脫粘層根部的影響。

圖4 遺傳算法優化歷程Fig.4 Optimization process of the GA

圖5 形狀最優藥柱Von Mises應變云圖Fig.5 Von Mises strain fringe of the optimal grain

算例(2)的優化結果與算例(1)的優化結果有較大差別，算例(1)得到的最優藥形體積裝填分數為0.968 0，而算例(2)得到的最優藥形體積裝填分數達到0.980 7，這僅因為脫粘深度Hd取值不同。文獻［10］中指出，脫粘深度對傘盤部位的應變影響很大，傘盤深度對傘盤最大Von Mises應變的影響規律隨脫粘深度的不同發生明顯改變。算例(1)中脫粘深度較小，得到的優化結果對應的傘盤深度已接近其取值范圍的上限，因為此時傘盤深度越大，最大Von Mises應變越小;而算例(2)中脫粘深度較大，優化結果對應的傘盤深度為其取值范圍的下限，因為此時傘盤深度越小，最大Von Mises應變越小。

本文2個算例的對比，充分說明脫粘深度對傘盤結構形狀優化的顯著影響。所以，在錐柱形藥柱設計過程中，應充分重視脫粘深度的合理選擇。為了減小低溫載荷作用下傘盤部位的最大Von Mises應變，在傘盤較深時，脫粘深度應取較小值，而在傘盤較淺時，脫粘深度應取較大值。

另外，雖然本文以藥柱的體積裝填分數最大為優化目標，但在藥柱設計過程中，體積裝填分數往往會受到通氣面積的限制。這是因為體積裝填分數的增大，必然導致通氣面積的減小，通氣面積過小，會影響發動機的總體性能指標。所以，將本文優化方法用于實際發動機藥柱設計時，應綜合考慮各方面因素的影響。

4 結論

(1)將遺傳算法與固體發動機結構完整性分析相結合，對錐柱形藥柱的傘盤結構進行了形狀優化設計，得到了體積裝填分數最大且滿足結構完整性要求的藥形。

(2)采用PCL語言對MSC.Patran軟件進行二次開發，確定了發動機參數化建模方法，編制的參數化建模程序能根據輸入參數自動建模計算，并輸出計算結果，保證了優化算法的順利進行。

(3)通過實例表明，對于錐柱形藥柱而言，不同的脫粘深度對應的優化藥形差別較大，在藥柱設計過程中，應充分考慮脫粘深度的影響。

(4)實際應用方面，基于參數化建模的藥形優化方法是可行和有效的，對各種藥形發動機都具有應用前景，具有較高的工程應用價值。

［1］陳汝訓.固體火箭發動機設計與研究［M］.北京:宇航出版社，1991.

［2］張國棟，顧克秋.基于遺傳算法和參數化建模的非線性結構優化［J］.計算力學學報，2003，20(6):764-768.

［3］陳珅艷，袁家軍，黃海.基于MSC.PATRAN/NASTRAN的結構優化程序系統［J］.宇航學報，2005，26(4):514-518.

［4］隋允康，李棟，杜家政，等.用二級控制法對二維連續體進行形狀優化［J］.計算力學學報，2007，24(2):135-141.

［5］陳繼華，隋允康，杜家政.基于滿應力準則的膜結構截面優化［J］.北京工業大學學報，2005，31(2):117-121.

［6］關振群，羅陽軍，王學軍，等.離心式葉輪三維參數化形狀優化設計方法［J］.機械強度，2006，28(6):833-838.

［7］宋鋒，溫衛東，崔海濤.基于改進蟻群算法的結構形狀優化［J］.航空學報，2007，28(5):1110-1115.

［8］蒙上陽，唐國金，雷勇軍.低溫環境下固體火箭發動機藥柱傘盤結構設計［J］.推進技術，2004，25(5):397-400.

［9］李磊，唐國金，雷勇軍，等.固體火箭發動機藥柱傘盤結構應力應變分析［J］.推進技術，2008，29(4):477-480.

［10］李磊，雷勇軍，申志彬，等.溫度載荷下傘盤深度與脫粘深度對藥柱應變的影響［J］.固體火箭技術，2010，33(3):285-288.