費 陽，解紅雨，張治宇，馮 翔

（中國人民解放軍96901 部隊，北京，100094）

0 引 言

裝藥燃面計算是固體火箭發(fā)動機內(nèi)彈道性能預(yù)示的重要內(nèi)容，其計算精度直接影響發(fā)動機優(yōu)化設(shè)計效果和性能預(yù)示精度。傳統(tǒng)的藥柱燃面計算方法包括作圖法、解析公式法和通用坐標(biāo)法；對于三維復(fù)雜藥型，可用CAD 實體造型法和數(shù)值算法等。CAD 實體造型法目前在發(fā)動機工業(yè)設(shè)計部門使用較普遍，需要對UG、SolidWorks 和ProE 等商業(yè)軟件進行二次開發(fā)，燃面設(shè)計時，需要專業(yè)人員干預(yù)操作；有學(xué)者利用最小符號距離算法和惠更斯原理，實現(xiàn)了復(fù)雜三維藥柱燃面計算，但是該方法在非并行計算條件下計算效率較低?；隗w素離散化的藥柱燃面計算方法，利用MC 界面重構(gòu)算法實現(xiàn)了復(fù)雜構(gòu)型裝藥燃面計算，但MC 算法的體素表存在二義性問題，在計算網(wǎng)格稍粗的情況下可能影響計算精度。Level-set 界面追蹤法通用性較強，但是該方法初始燃面定義十分復(fù)雜，需要手動進行代碼定義，因此藥型定義精度受限；另外，該算法追蹤燃面時，需要燃面附近多層“窄帶”網(wǎng)格點信息的迭代計算（重新初始化計算），此特性不僅降低計算效率，還導(dǎo)致燃面過早下降。

針對現(xiàn)有燃面計算方法的缺陷，本文利用SolidWorks 商業(yè)軟件對三維復(fù)雜構(gòu)型藥柱芯模進行建模，將芯模表面（初始燃面）生成STL 數(shù)據(jù)文件，將初始燃面三維笛卡爾坐標(biāo)系中進行離散描述，通過空間距離場算法計算空間離散點到初始燃面的最小符號距離，然后利用空間離散點最小符號距離信息快速計算燃面。此方法充分利用已有實體建模軟件，初始燃面定義快捷精確，燃面計算過程中無需迭代計算；針對最小符號距離計算量較大的問題，可將三維笛卡爾離散空間劃分成多個計算區(qū)域，進行多線程并行計算，極大提高計算效率。

1 燃面計算方法

1.1 總體思路

本文燃面計算方法總體思路如圖1 所示。

圖1 燃面計算總體思路Fig.1 The Outline for Burning Surface Calculation

首先利用SolidWorks或ProE等商業(yè)軟件對推進劑藥柱芯模進行實體造型，然后將芯模實體構(gòu)成初始燃面的表面生成STL 文件。構(gòu)造一個恰好能包絡(luò)藥柱、以{,,}為單位正交基底的三維笛卡爾空間區(qū)域，將空間區(qū)域在、、方向上離散成N ×N ×N個空間點，將方向設(shè)置為藥柱軸向，且空間坐標(biāo)軸原點與STL 文件坐標(biāo)軸原點重合。通過空間幾何關(guān)系分析，得到各空間點到初始燃面的最小符號距離Φ , 。然后，利用面積積分算法對燃面進行計算。假設(shè)燃燒肉厚為Δ，將各空間點最小符號距離Φ, 進行更新Φ , ,Φ , Δ，重復(fù)上述步驟，即可得到不同燃燒肉厚的燃面。

1.2 裝藥芯模初始燃面STL 文件

STL 文件格式用三角面片表達實體表面數(shù)據(jù)信息。ASCII 碼STL 文件格式如表1 所示，ASCII 碼格式STL 文件逐行給出三角面片的幾何信息，每行以1個或2 個關(guān)鍵字開頭。在STL 文件中，由多個三角形面片facet 信息構(gòu)成，每1 個facet 由7 行數(shù)據(jù)組成，第1 行facet normal 為三角形面片的單位法向量，在本文中，單位法向量由燃面指向藥柱內(nèi)部。第3 至第5行為三角形的3 個頂點坐標(biāo)，且單位法向量方向與3個頂點的依次順序遵循右手法則關(guān)系。

表1 初始燃面STL 文件格式Tab.1 The Data Format of the Initial Burning Surface STL File

1.3 空間最小符號距離計算

最小符號距離場是指在以{,,}為單位正交基底的空間直角坐標(biāo)系中，三維笛卡爾離散空間中（N ×N ×N），每一個空間離散點到燃面的最小符號距離Φ, ，Φ, ,> 0表示當(dāng)前空間點處于藥柱內(nèi)部，Φ , ,=0表示當(dāng)前空間點恰好位于藥柱燃燒表面，Φ, ,< 0表示當(dāng)前空間點位于燃燒室空腔。因此，除了計算最小距離之外，最小距離符號的判斷十分重要。利用STL 文件準(zhǔn)確計算初始最小距離場是實現(xiàn)燃面計算的關(guān)鍵。

燃面計算總體思路如圖2 所示。

圖2 燃面計算總體思路Fig.2 The Outline for Burning Surface Calculation

如圖2 所示，以三角形面片 ΔABC為例，設(shè)為其單位法向量，、為垂直于AB 邊并指向三角形外部的單位向量，、為垂直于BC 邊并指向三角形外部的單位向量，、為垂直于AC 邊并指向三角形外部的單位向量。這些有向線段和三角形面片將平面區(qū)域劃分成7 塊、共3 種情況的子區(qū)域。假設(shè)空間點為，討論不同情況最小符號距離計算方法。

1.3.1 點到面

如圖2 所示，若空間點P 點需滿足式（1），P 到Δ ABC 所在平面的投影在區(qū)域I（ ΔABC內(nèi)部），空間點P 到 ΔABC的最小距離就是P 到 ΔABC的垂直距離，可由式（2）計算，最小距離的符號可由式（2）直接確定。

1.3.2 點到點

如圖2 所示，當(dāng)空間點P 到 ΔABC所在平面的投影在區(qū)域III 時，空間點P 到 ΔABC的最小距離就是P 到Δ ABC 頂點的直線距離。若點P 需滿足式（3），則P點到頂點A 距離最近；若點P 需滿足式（4），則點P到頂點B 距離最近；若點P 需滿足式（5），則點P 到頂點B 距離最近。

多個三角形共點的情況下符號判斷見圖3。

圖3 多個三角形共點的情況下符號判斷Fig.3 Sign Identification when the Adjacent Triangles Share One Point

如圖3a 所示，當(dāng)頂點A 被多個三角形面片公用時，由式（2）確定的符號可能隨三角形的不同而變化，因此，需要構(gòu)造測試點T，T 到所有三角形的最小符號距離符號相同。如圖3b 所示，以A 點為起點，沿著、、、分別構(gòu)造單位向量′、′、′、′，測試點T 在坐標(biāo)系中的位置可由式（6）確定，然后利用式（7）確定符號。

1.3.3 點到邊

如圖2 所示，當(dāng)空間點P 到 ΔABC所在平面的投影在區(qū)域II 時，空間點P 到 ΔABC的最小距離為P 到Δ ABC 某條邊的垂直距離。若點P 需滿足式（8），則點P 到邊AB 距離最近；若點P 需滿足式（9），則點P 到邊BC 距離最近；若點P 需滿足式（10），則點P 到邊AC 距離最近。

三角形共邊情況下符號判斷見圖4。以點到BC邊的最小距離為例，如圖4a 所示，當(dāng)BC 邊被2 個三角形（ ΔABC、 ΔBCD）共享，且兩三角形平面夾角小于90°時（圖4b），若點P 處于1、3 區(qū)域，由 ΔABC、Δ BCD 判斷點P 的符號不一致。如圖4c 所示，在 ΔABC和 ΔBCD的夾角平分面上構(gòu)造測試點T，T 到 ΔABC、Δ BCD 有相同的距離符號。利用式（11）、式（12）在面ABC、面BCD 上構(gòu)造單位向量、，且、垂直于BC。記BC 中點為M，測試點T 坐標(biāo)可由式（13）確定，然后利用式（14）確定符號。

圖4 三角形共邊情況下符號判斷Fig.4 Sign Determination when the Adjacent Triangles Share One Edge

1.4 燃面積分

獲得計算空間內(nèi)各離散點的最小符號距離后，就可利用各種面積計算方法求得燃面。將包覆層等阻燃物質(zhì)轉(zhuǎn)化為數(shù)學(xué)描述的積分限，記燃燒室內(nèi)部區(qū)域為。對于復(fù)雜三維裝藥，直接引入面積計算公式（15），該計算公式被Level-set 計算程序廣泛采用，僅需要燃面附近1~2 層離散點信息：

式中為時刻燃面面積；下標(biāo)為計算區(qū)域中任意一點(,,) ；() 為當(dāng)?shù)刈钚【嚯x函數(shù)。式（15）中函數(shù)定義為

式中= 1.5min( Δ,Δ, Δ)。

2 算例分析

利用燃面計算方法對星孔型二維藥柱和翼柱型三維藥柱進行燃面計算，將計算結(jié)果與CAD 實體造型法進行對比，驗證計算方法的正確性和準(zhǔn)確性。計算所用工作站CPU（Xeon E3-1230）核心數(shù)為16，核心主頻為3.5 GHz。

具體過程為：對算例藥柱內(nèi)部空腔（即藥柱芯模）進行SolidWorks 實體建模，選取所有空腔表面（即藥柱初始燃面），以裝藥尾部端面中心為三維笛卡爾坐標(biāo)原點，軸和軸所在平面垂直于藥柱軸向，軸為裝藥軸向由尾部指向頭部，輸出STL 文件（文件格式如表1 所示），燃面計算模塊讀取STL 文件，在相同坐標(biāo)系下，根據(jù)藥柱外徑、長度信息自動構(gòu)建恰好能包裹裝藥的長方體計算空間，根據(jù)芯模最小特征尺度（如翼寬、星尖半徑等）合理設(shè)置計算網(wǎng)格尺度Δ、Δ、Δ，然后進行1.3 節(jié)所述的燃面計算過程。

2.1 星孔型二維藥柱

星孔型裝藥結(jié)構(gòu)如圖5 所示。

圖5 星孔型裝藥結(jié)構(gòu)示意Fig.5 A Structure Representation of A Star-shaped Propellant Grain

算例1 采用圖5 所示星孔型裝藥，兩端及外側(cè)包覆，內(nèi)表面燃燒。星孔型裝藥燃面計算雖有精確的解析公式，但解析公式無法計算燃面達到包覆層后面積迅速下降的過程，為全面檢驗本文方法準(zhǔn)確性，本算例與CAD 實體造型法進行對比。實體造型法被工業(yè)部門廣泛應(yīng)用，精度較高。初始燃面最小符號距離計算過程中，藥柱芯模STL 文件三角面片數(shù)量為1158，計算網(wǎng)格數(shù)量為300×300×500，將包覆層轉(zhuǎn)化為積分上限，若燃面超出包覆層，則不對其積分。計算采用32個并行線程數(shù)，計算時間為25.325 s。圖6 為星孔型裝藥燃面計算結(jié)果對比曲線，二者吻合較好，平均相對誤差為1.25%，最大相對誤差為3.36%。燃去肉厚達50 mm 時，燃面抵達側(cè)面包覆層，燃面迅速下降。

圖6 星孔型裝藥燃面計算結(jié)果對比Fig.6 A Comparison of the Burning Surface Results between this Method and Cad Modeling Method

2.2 翼柱型三維藥柱

為進一步考察燃面計算方法的通用性、準(zhǔn)確性，算例2 采用圖7 所示翼柱型裝藥，兩端及外側(cè)包覆，內(nèi)表面燃燒。利用本文方法對其平行層推移過程進行仿真。三維翼柱型裝藥燃面沒有解析公式，本算例與CAD 實體造型法計算結(jié)果進行對比。初始燃面最小符號距離計算過程中，藥柱芯模STL 文件三角面片數(shù)量為3152，計算網(wǎng)格數(shù)量為350×350×1000。同樣，程序中將包覆層轉(zhuǎn)化為積分限，若燃面超出包覆層，則不對其積分。計算采用32 個并行線程數(shù)，計算時間為127.559 s。

圖7 翼柱型裝藥結(jié)構(gòu)示意Fig.7 A Simplified Structure of A Back-wing Propellant Grain

圖8為翼柱型裝藥燃面計算結(jié)果對比。

圖8 翼柱型裝藥燃面計算結(jié)果對比Fig.8 A Comparison of the Burning Surface Results between this Method and Cad Modeling Method

二者吻合較好，平均相對誤差為1.33%，最大相對誤差為3.07%。圖中A-A 截面為后翼剖視圖。

3 結(jié)束語

本文對固體發(fā)動機藥柱芯模進行CAD 實體建模并生成STL 文件，利用最小符號距離算法，將STL 文件生成初始燃面空間最小符號距離場，通過最小符號距離場信息實現(xiàn)復(fù)雜三維藥柱燃面計算。計算結(jié)果與CAD 實體造型相比，吻合較好。利用多線程并行計算技術(shù)，極大提高了計算效率。該方法無需通過代碼定義藥柱構(gòu)型和迭代計算，可直接用于不同固體發(fā)動機藥柱設(shè)計。