固體火箭發動機裝藥設計的一種燃面推移方法研究 ①

王 棟，劉登航，卓長飛

(南京理工大學 機械工程學院，南京 210094)

0 引言

燃面變化規律是固體火箭發動機裝藥設計的一個核心問題，直接決定著發動機推力曲線走向[1-2]。目前已有的燃面推移計算方法包括解析法、作圖法、通用坐標法、實體造型法、網格法、level-set法、最小距離函數法等。其中，解析法、作圖法與通用坐標法存在對三維復雜藥型計算困難的局限性[1-6]。實體造型法可利用大型通用CAD軟件，但不易實現與流場耦合計算，對于不同藥型需要進行不同的推移構造和幾何尺寸定義，不具有通用性[7-13]。網格法通用性好，但對復雜裝藥在推移過程中出現部分型面的交匯、分離和變形的情況時誤差較大[14-19]。level-set法不需要顯式的追蹤活動界面，易于處理復雜界面及其在移動過程中發生的拓撲變形，但對不同藥型的初始型面，需要手工給出，或者需要將初始型面用非結構網格離散輸入，裝藥初始型面設置過程非常復雜[20-28]。最小距離函數法具有較強的通用性，但計算過程復雜，計算速度偏慢[29-35]。

目前，基于CAD軟件的燃面算法多只是計算燃面變化規律，少有與一維內彈道相結合揭示整個燃燒過程內流場變化情況的。本文采用離散網格思想，針對任意型面非槽平端的單通道裝藥或者實心裝藥，研究一種通用的燃面推移計算方法，該計算方法與一維內彈道相結合，其優點在于操作簡單，計算量小，基于C++編寫的控制臺程序計算速度快，具有良好的通用性，可以計算多種非槽的單通道裝藥或者實心裝藥，能夠較為準確捕捉燃面拓撲結構變化，可得到裝藥在燃燒過程中的一維內流場參數變化情況，與裝藥實際燃燒過程吻合度較高，對于固體火箭發動機設計具有良好的輔助預示效果。

1 裝藥初始表面數據獲取

1.1 裝藥建模

利用NX 12.0軟件的二次開發工具Block UI Styler模塊以及SNAP編程接口，編寫參數驅動的裝藥設計程序。該程序對于常見的非槽平端單通道裝藥或者實心裝藥給出參數驅動的建模方式，包括端面為平面的星孔裝藥、車輪孔裝藥、錐柱孔裝藥、梅花裝藥、圓柱孔裝藥、多臂形裝藥、十字型裝藥等；對于用戶自定義的非槽平端單通道裝藥或者實心裝藥模型也給出了模型導入接口，如圖1所示。

圖1 自定義裝藥模型導入界面

對于常見裝藥，用戶只需要給出參數，程序就可以自動建立三維模型。例如，要建立如圖2所示的車輪孔裝藥模型，在圖3所示的常見裝藥生成模型界面中選擇“車輪孔”，并輸入尺寸輪臂數5、裝藥外直徑40 mm、特征長度10 mm、裝藥長度500 mm、輪邊夾角70°、角分數0.55、輪臂高度4 mm以及3個過渡圓半徑均為0.5 mm。

圖2 車輪孔裝藥模型

圖3 常見裝藥模型生成界面

1.2 表面離散

進行裝藥表面離散是為了獲取初始表面數據。離散點越密集，越能夠反映藥真實物理形狀。本文的離散方法是在前人已有的射線旋轉相交法[18-19]的基礎上改進得到的射線搖擺旋轉相交法。射線旋轉相交法是從裝藥軸線上引出一條射線，與裝藥的內外表面相交求離散點，射線旋轉一周得到在該軸線位置處的一圈裝藥離散點，然后挪動射線到軸線上的下一位置，繼續旋轉求交，直到離散完整個裝藥內外表面。新方法射線搖擺旋轉相交法增加了相交射線的搖擺回旋操作，使其能夠處理裝藥離散面與射線有多個交點的情況。

參數設置如圖4所示，參數M用于將裝藥沿軸線均勻劃分為射線旋轉截面，參數N決定射線在截面處每一次旋轉的角度。

圖4 裝藥面離散參數輸入界面

以車輪孔裝藥內表面的離散為例進行說明，圖5所示為車輪孔裝藥在軸向某一截面的內表面輪廓。前人所使用的旋轉相交法在這種情況下射線直接逆時針從位置1經過位置2、3直接旋轉到位置4，無法處理中間有多個交點的情況；改進后的搖擺旋轉相交法令射線逆時針順時針交替旋轉，在位置1和4之間按照A-C-B-D-F-H-E-G的曲線順序搖擺求相交點。然后，在軸向下一截面位置處的內表面上繼續搖擺旋轉求相交點，直到搖擺旋轉相交完整個車輪孔裝藥內表面，得到該面上所有的離散點。

圖5 車輪孔裝藥內表面離散截面輪廓圖

對裝藥端面進行離散時，將端面看成是特殊的軸向截面，可用射線搖擺旋轉相交法求出端面的邊界離散點，再通過簡單的插值方法插入端面的非邊界點。示例車輪孔裝藥最終離散結果如圖6所示。

圖6 車輪孔裝藥離散結果

1.3 數據結構

為便于管理離散點數據，采用鏈表作為基本數據結構。裝藥側面上的一行鏈表點節是由包絡一圈的離散點集構成的，端面上的一行鏈表則是同一射線上的離散點集構成的。整個裝藥的數據結構如圖7所示，將一個面上的所有鏈表按一定順序編織成網格，再將裝藥的每個面作為節點串成一行面鏈表，所有面的網格組合起來就是整個裝藥的網格。

圖7 離散點數據結構

2 燃面推移

2.1 與一維內彈道耦合推移

本文所采用的推移方法是燃面上的離散點沿燃面法向推移一定距離(肉厚)，該肉厚由該點位置處的燃速乘以指定的推移時間步長得到，而每一推移時刻的各點燃速由一維內彈道計算結果來確定。一維內彈道控制方程組如下：

(1)

當某一推移時刻裝藥的燃面、通氣面積等幾何特性已知，根據推進劑性能、發動機尺寸，使用速度侵蝕燃燒公式將侵蝕效應考慮進去，以四階龍格庫塔法求解一維內彈道微分方程組式(1)，計算出此刻的燃氣質量流率，以及裝藥燃速、溫度、密度、壓力等各物理參數沿發動機軸線方向的分布，可得到實時推力。然后，根據一維內彈道計算結果得到下一時刻每一個離散點應當燒去的肉厚，再進行集合拓撲處理后，又可得到裝藥新的幾何特性，以進行新的一輪一維內彈道計算，直到裝藥燃燒完成。燃面推移和一維內彈道計算的結合，使得每一個時間步長的內彈道計算結果決定下一個推移步的燃面上每一點的推移肉厚，而下一個推移步的燃面推移結果又決定下一個推移步的內彈道計算結果。如此循環迭代，直到裝藥燃燒終止，可得到整個裝藥燃面推移過程中的內流場變化情況。

此外，在推移過程中還要根據設置的離散參數來判斷離散點的疏密情況。在本文研究中，以裝藥初始離散出來的平均點間距為參考，當推移過程中出現同一鏈表行上兩相鄰點距小于參考值某個小數倍時，需要在該行鏈表上進行“點刪除”處理；出現相鄰兩行鏈表間距小于參考值某個倍數時，要進行“鏈表行刪除”處理；出現同一鏈表行上大于參考值的某個倍數時，要用直線擬合插值法在兩點間進行“點插入”處理；出現相鄰兩行鏈表間距大于參考值某個倍數時，同樣要用直線擬合插值法在兩行鏈表間生成新的鏈表行進行“行插入”處理。

2.2 面積計算方法

裝藥的燃面推移過程中需要用到面積計算，而裝藥側面和端面的面積計算方法并不相同。

無論是外表面還是內表面的面積計算，都可基于側面上兩相鄰行鏈表構成的微元環面來進行，將側面上所有微元環面面積累加，就能得到該側面的面積。

圖8為一車輪孔裝藥內表面微元環面。記微元環面的一閉合車輪邊界輪廓邊長L1，另一閉合車輪邊界輪廓邊長L2，微元環面的母線長l，裝藥設計根數為N，由微元環面上的微元梯形ABCD，可推導出整個微元環面面積計算公式：

S微元環面=N·l·(L1+L2)/2

(2)

該式適用于任意裝藥側面的微元環面面積計算。

圖8 車輪孔裝藥內表面微元環面

計算裝藥的端面面積和軸向截面面積時，只需要用到內外表面上的離散點鏈表。該計算可以利用二維多邊形面積計算公式進行，端面上的離散點僅用作推移過程中的輔助定位推移。如圖9所示的一多邊形Ж，將各頂點按照逆時針方向標記，并設點Pi=(xi，yi)，則其面積計算公式為

(3)

該式即為二維多邊形面積計算公式，式中n表示頂點個數，索引需要與n進行模運算，有xn=x0和y-1=yn-1。

圖9 二維多邊形示意圖

利用式(3)可實現實心裝藥的截面或者端面的面積計算。對于有內孔的單通道裝藥截面，可將其看成外側邊所構成的多邊形與內孔邊所構成的多邊形組合，其截面積為外側多邊形面積減去內孔多邊形的面積。

2.3 推移過程“面交匯”問題的處理

在裝藥燃面推移的過程中，不同面之間會出現奇異交匯的情況，使得裝藥表面拓撲失真。本文針對非槽的實心裝藥和單通道裝藥，可將“面交匯”可歸納為三類：面交叉、面脫落以及面破碎。面交叉和面脫落只發生在鄰面之間，即只發生在側面與端面之間，且由于端面離散點不參與面積計算，故判定標準只需要考慮側面即可；而對于面破碎，通常由于裝藥幾何尺寸特點，則只發生在內外側表面之間，余藥的產生與此相關。

“面交叉”是指裝藥上側面的原邊界脫離了對應端面，且在側面內部形成了與端面新的接壤邊界。以某端面燃燒的圓柱孔裝藥為例進行說明，圖10所示為某一圓柱孔裝藥的半剖圖，當裝藥開始端燃，側面與端面會產生面交叉，比如圖10中紅線圈出位置。圖11為該面交叉情況的簡單示意圖。圖11中，箭頭方向表示端面燃燒推移方向，下同。此時，點A所在的原邊界鏈表節點行脫離了端面。處理方法為在側面找到另一鏈表節點行作為新的邊界，并按照擬合趨勢挪動到端面處以貼合端面，再將側面上多余的超出鏈表節點行刪去。對應在圖11中，則是刪去點A所在的鏈表行，點B所在鏈表行作為新的邊界。

圖10 圓柱孔裝藥半剖圖

圖11 面交叉示意圖

“面脫落”是指裝藥側面與端面的接壤邊界發生脫落，產生了縫隙，需要增加離散點使端面或者側面脫落增加面積，以此來保持兩相鄰面的邊界接壤。舉例如圖12所示，為一內孔燃燒的錐孔裝藥半剖圖，紅線圈出部分會出現面脫落。圖13為處理過程簡易示意圖。圖13中，箭頭表示內表面推移方向，虛線AA1表示將擴展增加的內表面部分。

圖12 錐孔裝藥半剖圖

圖13 面脫落示意圖

“面破碎”指的是裝藥原本不相鄰的內外表面產生了交集，使得這兩個面互相分割破碎。該情況多發生在裝藥燃燒即將結束階段。如圖14所示，是星孔裝藥內孔燃燒推移過程中某一時刻的情形，經過一次推移后，形成了面破碎，形成多塊余藥。圖15為該面破碎過程的截面輪廓圖，由該圖右邊截面輪廓可看出，裝藥內外表面輪廓線部分重疊，這些重疊部分在實際燃燒過程中已經不存在，剩下的未重疊部分形成了若干塊余藥。但面破碎情況不一定都是形成多塊余藥，比如內孔燃燒的圓柱孔裝藥由于侵蝕效應的存在，其裝藥末端會比裝藥頭部先燒完，末端處的裝藥內外表面相互交叉分割的面破碎過程并不會導致多塊余藥產生，只會使裝藥逐漸縮短。此時的處理方法是將實際物理已經不存在的重疊部分的離散點標記，在后續的燃面推移中不參與計算、不進行推移即可。

圖14 星孔裝藥面破碎

圖15 星孔裝藥面破碎輪廓圖

2.4 推移過程“奇異點”的處理

燃面推移過程中同一面上還會出現奇異點，使得計算失真。與面交匯一樣，只需要考慮側面即可。根據奇異點形成的原理不同，將其分為三類：錯位圓弧點、交叉點以及行奇異點。

“錯位圓弧點”是裝藥側面上的圓弧點，這些點經過一次推移后，穿過圓弧的圓心抵達另外一端，形成奇異錯位。舉例如圖16所示，為一內孔燃燒的星孔裝藥，紅線圈出的星根處容易在推移過程中出現本類奇異點。其產生過程如圖17所示，圖17中箭頭指示方向依舊為燃面推移方向，星根圓弧AB經過一次燃面推移到達圓弧A1B1產生了拓撲奇異錯位，這些奇異點即為錯位圓弧點。處理方法為用線段相交函數判斷識別并剔除該類奇異點。所示情況處理后的輪廓線如圖18所示。

圖16 星孔裝藥錯位圓弧點發生處

圖17 錯位圓弧點示意圖

圖18 錯位圓弧點處理后示意圖

“交叉點”主要出現在裝藥內外表面的尖角處，由于尖角的兩側推移法向夾角小于90°，在某一推移時刻經過一次燃面推移后，該處產生輪廓線奇異交叉，這種尖角上的奇異點就是交叉點。如圖19所示，一車輪孔裝藥內孔燃燒，在某一時刻產生了如圖19中紅線圈出部分的尖角。該尖角經過一個推移時間步長、離散點沿著燃燒方向推移一定肉厚后，會出現交叉點，如圖20所示。該類點的處理同樣是通過線段相交函數判斷識別，然后刪去交叉點以及交叉點之間的離散點，將剩下的離散點重新構成鏈表行。

圖19 車輪孔裝藥尖角產生

圖20 交叉點產生示意圖

“行奇異點”一般出現在相對于裝藥軸線有一定拔模傾角的裝藥側面上，且經常以整行鏈表奇異的形式出現。有的裝藥側面在軸線的不同位置有不同的拔模傾角，如果這恰好導致側面上相鄰部分沿燃燒推移法向產生相對推移，那么在推移過程中，側面就會出現一些鏈表行的排列順序與空間順序不一致的情況，將這種鏈表稱為奇異行，奇異行上的點稱為行奇異點。圖21為一內孔燃燒的錐柱孔裝藥，圖中紅線圈出部分即為在推移中會出現奇異行的位置。圖22顯示了該位置奇異行產生的過程，圖中輪廓線的點代表該位置處有一個行鏈表。判斷識別并刪除奇異行，按照原有空間順序編織側面上剩下的鏈表，使其重新組織為面網格，即可消除此類奇異。

圖21 錐柱孔裝藥半剖圖

3 實際應用驗證

3.1 算例

為檢測研究成果對模擬裝藥實際燃燒的計算效果，對某發動機裝藥進行了計算預測。該裝藥長度為730 mm，裝藥外表面為三瓣異型結構，最大外直徑為φ64.3 mm，內孔為三輻車輪狀，三維模型示意如圖23所示，其軸向截面輪廓結構見圖24。

圖23 裝藥三維示意圖

圖24 裝藥截面圖

已知該發動機裝藥主要成分有硝化棉、硝化甘油、氧化鎂、硬脂酸鋁等，其柱剖燃速測試結果見表1。該裝藥推進劑成分類似雙基藥，但由于各成分比例未知，本文選擇用雙芳鎂-3推進劑的性能參數作為參考估計，在表2中給出該裝藥推進劑其他性能參數的經驗估計值，其中Vth為速度侵蝕函數公式中的侵蝕界限流速，Kv為速度侵蝕函數公式中的系數。裝藥不包覆，兩側面兩端面均為燃面，其發動機設計參數裝藥根數為1，燃燒室內徑為φ71 mm，燃燒室長度為740 mm，噴管喉徑為φ27 mm，噴管出口直徑為φ68 mm。

表1 裝藥柱剖燃速測試結果

表2 裝藥推進劑部分性能參數估值

表1中，燃速公式計算所用壓強p的單位為MPa，依照公式計算得到的燃速單位為mm/s。該裝藥模型的離散參數設置周向分布點240個，軸向分布點個數選擇對應的裝藥設計程序的推薦值173，燃面推移計算中的推移時間步長設置為0.006 s，分別演算預示環境溫度為-55、20、60 ℃下的燃面推移結果，并與實驗測得結果進行對比。

3.2 結果分析

圖25～圖30分別為-55、20、60 ℃下的燃燒室頭部位置計算預測結果與實驗測量結果曲線圖。

圖25 -55 ℃下p-t曲線對比

圖26 -55 ℃下F-t曲線對比

圖27 20 ℃下p-t曲線對比

圖28 20 ℃下F-t曲線對比

圖29 60 ℃下p-t曲線對比

圖30 60 ℃下F-t曲線對比

分析上述三種工作環境溫度下的推移計算結果與實驗測量結果對比曲線，發現它們變化趨勢一致，且相比于實驗測量值，計算預計的壓強曲線最大峰值相對誤差不超過3.7%，推力曲線最大峰值相對誤差不超過6.8%，以壓強低于1.5 MPa為工作結束點計算的工作時間誤差不超過6.8%，可認為在該實驗中，燃面推移計算結果在一定精度范圍內能夠較好地預示裝藥實際燃燒情況，表明計算方法是正確的。計算中的離散點數和計算步長選擇能夠滿足計算精度要求，其帶來的誤差相對于圖中的主要誤差來源可忽略其影響。對比圖中計算預測曲線與實驗測得值存在較大誤差的主要來源，一方面是因為計算使用的是一維內彈道，忽略裝藥點火過程，直接認為裝藥初始狀態即為平衡壓強狀態，與實際的起始點不同，峰值有差異，且后效段只是近似處理，不一定符合實際；另一方面，則是因為計算用的推進劑性能參數是估計值。事實上，如果推進劑計算時用的性能參數是準確值，那么燃面推移計算結果將進一步契合裝藥實際工作燃燒情況，提供更好的預示效果。

4 結論

(1)基于離散網格法思想和幾何燃燒定律，結合一維內彈道提出了一種通用的模擬裝藥燃面推移的方法，能夠得到裝藥燃面推移過程中的一維內流場變化情況，且計算速度較快、精確度較好。在此基礎上，發展出了一種固體火箭發動機裝藥的通用設計方式，使裝藥設計過程更加直觀，簡化了設計流程。

(2)該方法對于各種藥形的非槽平端單通道裝藥或者實心裝藥，包括滿足上述條件的不同藥形組合裝藥，都能進行設計計算，覆蓋的裝藥藥形較廣，適用范圍大，通用性較好。

本文研究成果具有較大的工程應用價值，能夠適應多種藥形的復雜三維裝藥，為了解裝藥實際燃燒情況提供指導意義。但本文研究仍存在不足，對于非槽平端單通道裝藥或實心裝藥之外的裝藥，方法則不適用，有待于后續研究改進擴展。