基于三維模型空間離散的固體發(fā)動機裝藥燃面計算方法

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

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

0 引 言

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

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

1 燃面計算方法

1.1 總體思路

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

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

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

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

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

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

1.3 空間最小符號距離計算

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

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

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

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

1.3.1 點到面

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

1.3.2 點到點

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

1.3.3 點到邊

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

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

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

1.4 燃面積分

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

式中為時刻燃面面積；下標為計算區(qū)域中任意一點(,,) ；() 為當地最小距離函數。式（15）中函數定義為

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

2 算例分析

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

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

2.1 星孔型二維藥柱

星孔型裝藥結構如圖5 所示。

圖5 星孔型裝藥結構示意Fig.5 A Structure Representation of A Star-shaped Propellant Grain

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

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

2.2 翼柱型三維藥柱

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

圖7 翼柱型裝藥結構示意Fig.7 A Simplified Structure of A Back-wing Propellant Grain

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

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

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

3 結束語

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