蔡 強，鮑福廷，丁 林，劉 旸

(西北工業(yè)大學 航天學院，西安 710072)

0 引言

固體火箭發(fā)動機翼柱形藥柱具有三維復雜構型，其燃面面積很難用解析公式表達。實體造型方法［1］是目前采用的主要方法，通常對商用CAD軟件進行二次開發(fā)，避免了復雜的數理方程推導，卻容易因特征定義不嚴格，而出現特征消失或自相交，導致幾何模型拓撲畸變而造型失敗，且商用CAD軟件計算大型藥柱的時間較長，不利于方案的驗證與藥柱的優(yōu)化設計。文獻［2］中的方法僅能處理具有1組翼的藥柱，其采用基于造型歷史追溯的燃面算法，容易發(fā)生統(tǒng)計疏漏和重復。

為了解決上述問題，本文結合翼柱形藥柱構型的特點，引入變量化設計方法，在幾何造型內核上自主開發(fā)。通過定義統(tǒng)一的翼形體特征函數，給出翼形體的關鍵設計變量，嚴格控制模型參數的變化，避免了造型失敗;調節(jié)三維圖形變換矩陣中的比例變換因子，縮小幾何模型的尺寸，顯著縮短了大型翼柱形藥柱的仿真時間;本文還研究新的燃面面積算法，避免統(tǒng)計燃面出錯，并針對后翼、前后翼和大小翼等各種類型的翼柱形藥柱，總結了統(tǒng)一的設計流程。

1 翼柱形藥柱變量化仿真建模

1.1 假設條件

(1)整個翼柱形藥柱表面光滑均一，沒有裂紋、尖點等缺陷。

(2)翼柱形藥柱的燃面同時被點燃，燃面各點的條件相同。

(3)燃面退移滿足平行層燃燒定律，即燃面各點都以相同燃速向藥柱內推進。

(4)翼柱形藥柱的各組翼形體周向均勻分布，滿足嚴格的循環(huán)對稱條件。

1.2 建模仿真過程

本文將翼柱形藥柱分為外輪廓、回轉體和翼形體3種特征形體，其變量化設計過程可分為變量化約束草圖繪制、特征形體定義與生成和燃面退移仿真3步，如圖1所示。首先，交互式繪制外輪廓、回轉體和翼形體的二維截面草圖，并對草圖施加約束，定義設計變量與變量方程。然后，根據藥柱相應特征的二維草圖，采用幾何造型算法生成外輪廓、回轉體和翼形體的三維幾何實體，外輪廓依次“減去”回轉體和翼形體，就得到藥柱的初始形狀，并調用統(tǒng)計函數獲得初始燃面面積及質量特性信息。最后，仿真燃面退移時，按照平行層燃燒定律假設，每次退移一定的肉厚，通過改變回轉體和翼形體的特征變量，使得變化后的各形體表面與初始表面等距，外輪廓“減去”變化后的回轉體和翼形體，得到一系列退移后的藥柱形狀，計算燃燒面積和質量特性，得到燃面面積、質心、慣量、形狀的變化規(guī)律。

圖1 翼柱形藥柱變量化設計過程Fig.1 Variational design principle for finocyl grain

2 翼柱形藥柱變量化設計關鍵技術

在翼柱形藥柱的實際應用過程中，采用在藥柱前段增加1組軸向前翼，或在后段周向增加另1組小翼的方法來增大初始燃面，從而為助推段提供更大推力，前者稱為前后翼，后者稱為大小翼。從藥柱二維截面生成三維特征形體以及燃面退移時，為了簡便，通常只處理循環(huán)對稱藥柱的1個對稱單元，然而前后兩組翼通常不在同一對稱平面內，大小2組翼是周向陣列而成，因此需采取措施，確保所有翼形體都處在一份對稱單元中，從而得到合理的布爾運算結果。另外，還要解決翼形體的通用描述方法、穩(wěn)健準確的燃面面積算法和縮短仿真時間等關鍵技術。

2.1 圖形比例變換

本文采用的是支持邊界表達方法的精確幾何造型內核［3］，藥柱的尺寸規(guī)模越大，耗費機時越長。因此，一個直觀的想法就是在保證造型精度的條件下，對外輪廓、回轉體和翼形體的相應二維截面進行圖形比例變換，以縮小尺寸［4］。這樣生成的三維特征形體尺寸較小，能顯著縮短后續(xù)的布爾運算和等距偏移運算所耗費的時間。

式(1)為三維圖形比例變換矩陣:

式中 (x，y，z，1)為原始截面各點的齊次坐標;(x*，y*，z*，1)為縮放后的特征截面相應點的齊次坐標;(xf，yf，zf)為比例變換的參考點坐標;Sx、Sy和 Sz分別為x、y和z坐標分量的比例變換因子，比例變換因子的取值以不損失造型精度的前提下越小越好。

翼柱形藥柱特征形體原始截面的各點依次乘以比例變換矩陣，就得到特征形體縮放后截面的各點坐標。

2.2 翼形體特征函數定義

考慮到通用性，本文定義翼形體特征函數如下:

式中 Section為翼形體二維特征截面;w為翼形體厚度;r為翼形體頂端倒角半徑(見圖2);angle為該組翼形體起始角度;symN為藥柱的循環(huán)對稱數目;npc為每個外輪廓對稱單元中該組翼的個數。

圖2 翼形體示意圖Fig.2 Parameters of fin entity

藥柱的循環(huán)對稱數目symN是前后2組翼個數(nf，na)或大小2組翼個數(nb，ns)的最大公約數，則

翼在外輪廓對稱單元中的周向排列如圖3所示。翼形體起始角度angle表示第一個翼距外輪廓對稱單元中心線的角度，其余翼的角度為

圖3 外輪廓對稱單元中多個翼的角度分布Fig.3 Angles of multiple sets of fins defined in one cyclic symmetry cell

每組翼的第1個翼作為幾何造型的樣本，其余通過將樣本的副本繞X軸進行三維圖形旋轉變換得到［4］，即

式中 (x，y，z，1)為樣本點的齊次坐標;(x*，y*，z*，1)為副本點的齊次坐標。

經過如上定義以后，即可確保所有翼形體均處在外輪廓實體對稱單元中，從而得到嚴格旋轉對稱的1/symN份藥柱實體，然后周向陣列即可得到完整的藥柱信息。

2.3 顏色標識的燃面面積計算方法

本文采用邊界表達的幾何模型，其數據結構完整記錄了模型的拓撲、幾何和屬性信息［5］，顏色作為基本屬性，附加在藥柱外輪廓實體和內腔特征實體的幾何模型上，并隨著實體模型之間的布爾運算及圖形變換，保持其合理性，具體體現在:

式中 藥柱實體Sgrain的表面來源于2部分:由外輪廓實體Sout的幾何模型表面分裂而成，記為SFa;由內腔特征實體Score的幾何表面分裂而成，記為SFb。

本文用不同的顏色標記燃面和非燃面:

(1)外輪廓實體所有表面標識為藍色，內腔特征實體表面標識為紅色。因此，所有的SFa為藍色，所有的SFb為紅色。

(2)按照式(6)進行布爾差運算，得到的藥柱實體的表面中，源自外輪廓實體的為藍色，源自內腔特征實體的為紅色。顯然，紅色的面即為翼柱形藥柱的燃面，統(tǒng)計所有紅色面的面積，即為藥柱的初始燃面面積。

(3)內腔特征實體擴張以后，仍保持所有表面為紅色。同理，與不發(fā)生變化的外輪廓實體進行布爾差運算之后，統(tǒng)計得到的擴張后藥柱實體的紅色表面面積，即為退移后的燃面面積。

這樣，只需在初始時進行1次顏色標識，就可保證燃面統(tǒng)計不發(fā)生疏漏，得到精確的計算結果。

3 翼柱形藥柱的變量化設計流程

在解決以上關鍵技術后，整個翼柱形藥柱的變量化設計流程總結見圖4。

圖4 翼柱形藥柱的變量化設計流程Fig.4 Variational design flow chart of finocyl grain

(1)在變量化設計環(huán)境中，交互式繪制外輪廓、回轉形內孔、前后翼或大小翼的約束草圖，得到各個特征實體的原始二維截面。

(2)按照2.1節(jié)的方法對特征實體的截面執(zhí)行比例變換，縮小圖形的尺寸規(guī)模，以節(jié)省求解時間。

(3)由特征截面分別生成外輪廓的對稱單元、回轉形內孔的對稱單元、前后翼或大小翼的三維實體，并按照2.2節(jié)中的處理方法，確保所有翼形體在外輪廓對稱單元中。

(4)按照2.3節(jié)中介紹的顏色標識的燃面面積計算方法，計算初始燃面面積(肉厚Web=0)和退移一定距離(肉厚Web＞0)以后的燃面面積，每退移一步依次計算燃面面積及質量特性。

(5)判斷藥柱是否完全燃燒，即判斷當前藥柱質量是否為0，如果完全燃燒，就輸出燃面面積及質量特性隨燒去肉厚的變化關系;否則，繼續(xù)執(zhí)行退移仿真。

4 實例驗證

本文以LEDAS約束求解器［6］，構建變量化設計環(huán)境，實現了上述針對翼柱形藥柱的處理方法。

4.1 可行性實例驗證及結果分析

為了驗證該方法的可行性，以1個有3組翼形體的藥柱作為示例，截面尺寸如圖5所示。3組翼分為前翼、后大翼、后小翼，其參數見表1。

圖5 大小翼和前后翼藥柱二維截面草圖Fig.5 2-D geometric constraint sketch of rear big fin and rear small fin

表1 3組翼的參數Table 1 Parameters of three sets of fins

藥柱燃面退移過程及燃面面積-肉厚關系曲線如圖6和圖7所示。

圖6 大小翼和前后翼藥柱燃面退移過程仿真及燃面面積-肉厚曲線Fig.6 Grain burnback simulation and profile of burning surface area vs web

這個燃面退移過程分為3段，初始階段燃面不斷增大，直至肉厚Web=30 mm時，后大翼開始燒穿;Web=30～60 mm為減面燃燒，前翼和后小翼逐漸擴張，后大翼逐漸消失;Web=60 mm時，后大翼完全消失，前翼和后小翼開始燒穿，直至Web=85 mm藥柱完全消失，這一階段也為減面燃燒。仿真時，比例變換因子Sx=Sy=Sz=1，肉厚推進步長為1 mm，85步計算共耗時221.63 s，如果令 Sx=Sy=Sz=0.02，即縮小 50 倍，計算時間縮短至151.16 s(CPU主頻2.8 GHz)。

圖7 參考文獻實例的燃面退移仿真過程及燃面面積-肉厚曲線Fig.7 Grain burnback simulation and profile of burning surface area vs web

4.2 計算精度實例驗證及結果分析

為了驗證本文計算方法的精度，對文獻［7］中17#藥柱進行燃面退移仿真，并結合推進劑參數和噴管參數，進行了內彈道計算。整個藥柱肉厚為21.35 mm，燃燒初始階段(Web=0～4 mm)內孔和翼形體同時擴張，為增面燃燒;在Web=4 mm時，后端翼形體開始燒穿，直至Web=12 mm翼形體完全燒完，這一階段為減面燃燒;燃燒最后階段(Web=12～21.35 mm)為近似等面燃燒過程。該例本文設置的比例變換因子為Sx=Sy=Sz=2/3，仿真肉厚步長與文獻［7］相同，都為1 mm，22步仿真時間僅為7.344 s，而文獻［7］中仿真耗費了將近30 min。

本文內彈道計算結果與文獻［7］中壓力-時間進行對比，如圖8所示。由圖8可看出，由于本文內彈道計算未考慮點火段和拖尾段的非定常效應，點火段和拖尾段與試驗曲線差距較大，其余吻合相當好(最大誤差僅3.1%)。

圖8 壓力-時間曲線計算結果與試驗曲線對比Fig.8 Results of pressure vs time compared with experiment data

5 結論

(1)引入翼的起始角度，合理布置大小翼、前后翼的位置，使其在1個外輪廓實體對稱單元中。調節(jié)三維比例變換矩陣的比例變換因子，以縮小藥柱的尺寸規(guī)模，顯著縮短求解時間。

(2)在基于造型歷史追溯的燃面面積統(tǒng)計方法的基礎上，提出改進的顏色標識燃面面積計算方法，該方法具有執(zhí)行效率高、操作簡單、不會統(tǒng)計疏漏等優(yōu)點。

(3)歸納總結了翼柱形藥柱的變量化設計流程，包括特征形體變量化約束草圖的定義、比例變換、燃面退移過程仿真、統(tǒng)計燃面面積等步驟。并通過2個實例，對其有效性及計算精度進行了驗證。

(4)文中研究工作為翼柱形藥柱的設計、參數優(yōu)化和燃面退移仿真提供一套實用的解決方法。同時，本文的研究思路也適用于其他復雜的任意三維藥型。

［1］董新剛，陳林泉，侯曉.基于Pro/E平臺下的固發(fā)裝藥CAD軟件［C］//昆明:2002年中國宇航學會固體推進專業(yè)委員會年會論文集(上)，2002:109-114.

［2］劉旸，鮑福廷，蔡強，等.變量化約束草圖驅動的藥柱三維特征設計研究［J］.固體火箭技術，2010，33(1).

［3］蔡強，鮑福廷.基于ACIS幾何造型平臺的固體火箭發(fā)動機裝藥設計［J］.固體火箭技術，2008，31(3).

［4］孫家廣，等.計算機圖形學(第三版)［M］.北京:清華大學出版社，1998:358-373.

［5］詹海生，李廣鑫，馬志欣.基于ACIS的幾何造型技術與系統(tǒng)開發(fā)［M］.北京:清華大學出版社，2002:30-49.

［6］LEDAS Ltd.LEDAS geometric solver 2D overview［R］.http://ledas.com/products/lgs2d.

［7］Püskülcü G，Ulas A.3-D grain burnback analysis of solid propellant rocket motors:Part 2—modeling and simulations［J］.Aerospace Science and Technology，2008(12):585-591.