費(fèi) 陽，解紅雨，張治宇，馮 翔

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

0 引 言

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

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

1 燃面計算方法

1.1 總體思路

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

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

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

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

STL 文件格式用三角面片表達(dá)實(shí)體表面數(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 個頂點(diǎn)坐標(biāo)，且單位法向量方向與3個頂點(diǎn)的依次順序遵循右手法則關(guān)系。

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

1.3 空間最小符號距離計算

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

1.3.1 點(diǎn)到面

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

1.3.2 點(diǎn)到點(diǎn)

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

多個三角形共點(diǎn)的情況下符號判斷見圖3。

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

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

1.3.3 點(diǎn)到邊

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

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

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

1.4 燃面積分

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

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

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

2 算例分析

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

具體過程為：對算例藥柱內(nèi)部空腔（即藥柱芯模）進(jìn)行SolidWorks 實(shí)體建模，選取所有空腔表面（即藥柱初始燃面），以裝藥尾部端面中心為三維笛卡爾坐標(biāo)原點(diǎn)，軸和軸所在平面垂直于藥柱軸向，軸為裝藥軸向由尾部指向頭部，輸出STL 文件（文件格式如表1 所示），燃面計算模塊讀取STL 文件，在相同坐標(biāo)系下，根據(jù)藥柱外徑、長度信息自動構(gòu)建恰好能包裹裝藥的長方體計算空間，根據(jù)芯模最小特征尺度（如翼寬、星尖半徑等）合理設(shè)置計算網(wǎng)格尺度Δ、Δ、Δ，然后進(jìn)行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)表面燃燒。星孔型裝藥燃面計算雖有精確的解析公式，但解析公式無法計算燃面達(dá)到包覆層后面積迅速下降的過程，為全面檢驗(yàn)本文方法準(zhǔn)確性，本算例與CAD 實(shí)體造型法進(jìn)行對比。實(shí)體造型法被工業(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%。燃去肉厚達(dá)50 mm 時，燃面抵達(dá)側(cè)面包覆層，燃面迅速下降。

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

2.2 翼柱型三維藥柱

為進(jìn)一步考察燃面計算方法的通用性、準(zhǔn)確性，算例2 采用圖7 所示翼柱型裝藥，兩端及外側(cè)包覆，內(nèi)表面燃燒。利用本文方法對其平行層推移過程進(jìn)行仿真。三維翼柱型裝藥燃面沒有解析公式，本算例與CAD 實(shí)體造型法計算結(jié)果進(jìn)行對比。初始燃面最小符號距離計算過程中，藥柱芯模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ā)動機(jī)藥柱芯模進(jìn)行CAD 實(shí)體建模并生成STL 文件，利用最小符號距離算法，將STL 文件生成初始燃面空間最小符號距離場，通過最小符號距離場信息實(shí)現(xiàn)復(fù)雜三維藥柱燃面計算。計算結(jié)果與CAD 實(shí)體造型相比，吻合較好。利用多線程并行計算技術(shù)，極大提高了計算效率。該方法無需通過代碼定義藥柱構(gòu)型和迭代計算，可直接用于不同固體發(fā)動機(jī)藥柱設(shè)計。