烏 岳，李 卓，魯 榮

(1.內蒙古工業(yè)大學理學院，內蒙古 呼和浩特 010051；2.中國航天科工集團 內蒙古航天紅峽化工有限公司，內蒙古 呼和浩特 010051；3.內蒙古工業(yè)大學化工學院，內蒙古 呼和浩特 010051)

引 言

固體推進劑藥漿兼具黏性體和彈性體的雙重特性，是典型的黏彈性非牛頓流體。固體推進劑低黏度藥漿的澆注通常采用真空花板澆注工藝[1]，推進劑藥漿流動的驅動力為料斗和澆注缸之間的壓差以及料斗內推進劑藥漿自重[2]。其本構方程和流動邊界條件的復雜性，很難采用理論解析方法分析其流動[3]。目前隨著計算機輔助技術的迅速發(fā)展，數(shù)值模擬已成為分析黏彈流體流動的主要手段。

姜愛民等[4]對固體推進劑藥漿在混合釜內的流動過程進行了數(shù)值模擬，結合網格疊加技術和網格加密技術，對混合過程中藥漿流動進行了有限元模擬，詳細研究了捏合機混合過程中混合釜內剪切應力、各點的速度、混合指數(shù)、黏度的分布和極值關系；Connelly C K等[5]利用粒子示蹤法，模擬了行星式混合內玉米糖漿的三維混合流動；Shimada T等[6]利用Polyflow軟件將混有鉛示蹤顆粒的推進劑藥漿澆注到一個透明型腔中，并在周圍架設了兩個方向的X射線透析儀器，通過計算得到藥漿內顆粒的三維跡線，從而模擬得出比較直觀的藥漿內部的流場。

本研究利用黏彈性流體軟件Polyflow提供的廣義牛頓流體的本構模型，對復合推進劑藥漿真空花板澆注流動的全過程進行數(shù)值模擬，針對真空花板澆注系統(tǒng)不同孔徑的藥漿流動速度、擠出脹大效應、藥漿黏度和剪切速率及流平性做了一系列研究，同時通過實驗對模擬結果進行了驗證。

1 數(shù)值模擬

計算過程采用標準的CFD模擬分析流程：模型導入-劃分網格-設置模型及參數(shù)-設置邊界條件-利用CFD-POST對結果進行后處理。

1.1 推進劑澆注流動過程的基本假設

由于Polyflow軟件本身的局限性，本研究依據(jù)實際情況作了以下假設：

(1)流動為低雷諾數(shù)的不可壓縮流體的層流流動[7]，且流動過程不考慮推進劑藥漿與外界環(huán)境的熱交換(實際的澆注工房有循環(huán)水作保溫)；

(2)推進劑藥漿與壁面之間作無滑移邊界條件設置；

(3)由于Polyflow軟件本身的設定，體積分數(shù)云圖中若體積分數(shù)在0.5以上，則認為充滿。

1.2 計算模型參數(shù)及工作原理

真空花板澆注工藝裝置如圖1所示。

工作原理：將混合好的可澆注藥漿定量放入料斗中保溫(溫度50℃)，開啟膠管閥，藥漿按指定下料速度通過膠管閥進入真空度符合要求的燃燒室，將已除氣的藥漿通過真空壓差和自重流入發(fā)動機燃燒室內壁和芯模之間的空隙內[8]。

計算模型在圖1的基礎上作了必要的簡化設置，如圖2所示。

圖2 計算模型Fig.2 Calculation model

選取各個裝置的參數(shù)如下：料斗內徑90mm，花板直徑90mm，孔簡化為36個，孔的直徑依次為3、5、7mm(實際中同一型號的花板孔徑一致，為分析不同孔徑的流速及擠出脹大效應，假設同一花板有不同孔徑的分布)，澆注缸內徑200mm，燃燒室內徑100mm，燃燒室采用后翼、環(huán)型槽結構[9]，12個翼片。藥柱芯孔最大處直徑約40mm，前機口直徑約15mm，后機口直徑40mm，燃燒室總長200mm，直徑100mm。

1.3 網格劃分

本研究中計算域形狀不規(guī)則，為了提高花板網格質量將整個模型分割為38個體，分別對其進行結構化網格劃分。總網格數(shù)為442798個。網格劃分情況如圖3所示。

1.4 本構方程

研究表明[9-10]，固含量較高的推進劑藥漿往往表現(xiàn)出非線性賓漢流體特征：當剪切應力低于屈服值時體系不流動，而高于屈服值時，剪切應力與剪切速率呈非線性冪律關系，且黏度隨剪切速率增大而降低，即屈服假塑性流體[11]，亦稱為Herschel-Bulkley塑性流體。

圖3 網格劃分Fig.3 Mesh division

本研究針對廣義牛頓流體等溫流動模型Herschel-Bulkley進行計算，形式如下[11]：

圖4 廣義賓漢流體的3類塑性流動Fig.4 Three types of plastic flow of generalized Bingham fluids

文獻[11]表明，該固體質量分數(shù)為88%的四組元HTPB推進劑為典型的屈服假塑性流體，當施加的應力超過屈服應力后，體系發(fā)生流動，此時流動符合n<1的冪律定律。

1.5 邊界條件

本研究在數(shù)值計算中涉及的邊界條件有：流動入口邊界條件、流動出口邊界條件和壁面邊界條件。

工程文件任務一的流動入口采用作用法向應力fn和切向應力fs,由于入口邊界大氣壓強為0.1MPa，即fn=1.0×105，fs=0；流動出口邊界采用作用法向速率vn和作用切向速率vs，且均為零；其余為壁面邊界條件；由于選用有限元瞬態(tài)體積分數(shù)VOF(volume of fluid)模型，系統(tǒng)自動生成任務二Fluid Fraction Transport task，該任務邊界條件入口采用Inflow邊界條件，其余為壁面設置。

2 數(shù)值模擬結果與分析

2.1 花板孔徑對藥漿流動速度的影響

當黏彈性流體流經不同直徑的花板孔時，其流動速率會有很大不同。圓截面不同孔徑速度云圖如圖5所示。

圖5 圓截面不同孔徑速度云圖Fig.5 The cloud map with different orifice speeds in the circle section

由圖5可以看出，孔徑越大，流體的流速越大，且最大值出現(xiàn)在圓孔的中心軸線上：7mm孔徑條件下，流速最大值約0.10m/s；5mm孔徑流速的最大值約0.04m/s；3mm孔徑流速最大值約為0.02m/s，這是由于孔徑越大，流動阻力越小，質量流量越大導致流速越大，其次考慮到3種不同孔徑，壓強的相互競爭作用也會導致出現(xiàn)速度差異。

一般認為，與圓孔管壁處接觸的流體層可視為靜止的，因此圓孔管壁處流體流速接近于0[12]。

2.2 花板孔徑擠出脹大效應分析

擠出脹大效應(BARUS效應)是非牛頓流體的典型流變特征[13]。圖6為流體出口擠出截面，其大小用擠出脹大比(B=Dmax/D)表示。

由圖6可以看出，由于花板厚度為5mm，所以3個不同孔徑(3、5、7mm)的長徑比依次為5∶3、5∶5、5∶7。通過Matlab軟件編寫程序擬合圓孔擠出截面直徑，小圓(3mm)、中圓(5mm)、大圓(7mm)擠出截面直徑擬合結果分別為3.1371、6.5546和10.08mm。即擠出脹大比依次為1.05、1.31和1.44。

圖6 圓孔擠出截面Fig.6 Export extrusion cross section

由此可見，花板孔徑的長徑比強烈影響擠出脹大比。擠出脹大比隨孔徑長徑比減小而增大，當孔徑長度一致時，孔徑越大擠出脹大比越大。

2.3 黏度和剪切速率云圖

HTPB丁羥四組元推進劑為典型的非線性屈服假塑性流體，其主要特征是在某一剪切速率范圍內，黏度隨剪切速率的增大而降低，即剪切稀變流體，從黏度和剪切速率云圖中可以得到證實。

圖7為某一YZ截面的黏度云圖和剪切速率云圖，圖8為某一XY圓截面翼片位置藥漿的黏度云圖和剪切速率云圖。所采用藥漿初始黏度為400Pa·s。

圖7 某一YZ截面藥漿的黏度云圖和剪切速率云圖Fig.7 The viscosity and shear rate clouds for a YZ section slurry

從圖7和圖8可以看出，藥漿以較高的剪切速率流經花板和燃燒室翼片時，黏度降至最小值300Pa·s，當離開噴孔和翼片時，所受的剪切速率迅速衰減。產生以上結果的原因為：

(1) 流體流經花板不同孔徑時，剪切應力和剪切速率大幅度提升，這是因為藥漿從大面積料斗流道進入小截面圓孔時，流動狀況發(fā)生巨大變化，入口區(qū)藥漿會受到拉伸作用，出現(xiàn)明顯的流線收斂現(xiàn)象，進而導致藥漿不同層流之間的摩擦生熱加劇，使得藥漿連續(xù)相分子間熱運動動能增大，從而提高流體的流動性，導致黏度下降[14]；

圖8 某一XY圓截面翼片位置藥漿的黏度云圖、剪切速率云圖和圓截面不同孔徑速度云圖Fig.8 The viscosity and shear rate cloud maps of the fin position for a certain XY circular section

(2) 流體流經燃燒室翼片時，固體推進劑連續(xù)相的黏合劑在翼片流動剪切作用下其高分子鏈段發(fā)生取向，這些高分子鏈時而相互纏繞，時而解開，從而分子之間的平衡力受到破壞，在高剪切速率下這種破壞難以在短時間內恢復，進而破壞了藥漿內部黏滯結構，從而導致層流之間的運動黏滯度下降，黏度變小[15]。

此外，從計算模型的角度分析，所選用的本構方程為H-B方程，當流體超過臨界剪切速率時，流動過程遵從假塑性流體冪律變化規(guī)律。除花板和翼片的其他位置的剪切應變速率很低，進一步證實藥漿經真空花板澆注系統(tǒng)時，燃燒室中的流動是在低剪切應力和低剪切應變速率下的流動。

城市系統(tǒng)的模擬需要嵌入不確定的因素或者用戶期望的因素，從而模擬出不確定性的城市系統(tǒng)或者用戶所預期的城市形態(tài).傳統(tǒng)地理信息系統(tǒng)(GIS)在處理地理現(xiàn)象的時間過程上存在一定的局限性，而許多研究表明[1]，元胞自動機(CA)能更容易模擬各種現(xiàn)象隨時空變化的動態(tài)性，這是因為CA非常適合模擬復雜系統(tǒng).將GIS和CA結合起來，能夠更好的模擬真實城市的發(fā)展，提高模擬的精度.

2.4 藥漿流平性分析

將藥漿離開花板噴口至充滿燃燒室空間并淌平的過程稱為藥漿流平過程。 不同時刻藥漿的體積分數(shù)云圖如圖9所示。

圖9 不同時刻藥漿的體積分數(shù)云圖Fig.9 Volume fraction cloud chart of the slurry at different time

由圖9可知，藥漿在燃燒室均勻緩慢地從中間向四圍展開，在真空壓差和自重的綜合作用力下逐漸充滿燃燒室空間，這一過程與實際流動規(guī)律極為相似。藥漿在4s左右通過花板流入燃燒室，整個過程流平性良好，直到t=46s時，燃燒室芯孔右側出現(xiàn)一微小空洞，但其隨著流動逐漸消失，直到燃燒室充滿，孔洞的體積分數(shù)在0.5以下，可以忽略不計，基本充滿。

本算例所采用藥漿的屈服值為47.1Pa，從體積分數(shù)云圖可以看出該藥漿的流平性很好，表明推進劑在流平過程中屈服值越小流平性越好。

3 實驗驗證

本實驗研究對象為某固體質量分數(shù)88%的端羥基聚丁二烯(HPTB)四組元推進劑，配方(質量分數(shù))為：黏合劑，10.5%；高氯酸銨，60%；含能炸藥，10%；金屬粉末，18%；其他功能助劑，1.5%。密度為1.78×103kg/m3。

藥漿在混合出料后約20min內通過落球黏度計進行黏度測定。推進劑藥漿起始黏度值較低，為400Pa·s，模擬中初始黏度的取值與其對應。1h后藥漿黏度升至759.12Pa·s，之后每1h測定一次，藥漿黏度增長平穩(wěn)，說明藥漿整體性能良好。

3.1 模擬實驗器澆注前準備

為盡可能真實地還原真空花板澆注實驗，按照發(fā)動機燃燒室結構尺寸作如圖10的模擬實驗器。內徑100mm，高200mm，體積為0.00157m3。與模擬選用的燃燒室尺寸一致。

圖10 模擬實驗器Fig.10 Simulation experimental device

該模擬實驗器采用6mm厚度普通結構鋼板焊接而成。裝配前，首先用清理劑清理內表面，晾干后在內表面涂敷防黏劑，之后將試驗器各個零件按照先底板后側板的順序進行裝配，零部件之間加裝密封墊。

3.2 模擬實驗器澆注

澆注工裝裝配從下而上，依次為模擬實驗器、套筒、花板、膠管閥和料斗等?；ò宀捎忙?mm×36mm規(guī)格，進行一次除氣。

在澆注過程中，需有效控制膠管閥門開度，保證藥漿下料速度，通過觀察事先安裝的線性速度傳感器，顯示流體經過3mm孔徑時，速度為0.0158m/s，與模擬結果誤差在允許范圍內。記錄四組元配方保溫溫度為50℃，記錄澆注時U型壓力機的讀數(shù)為101000Pa，澆注缸內的余壓為4000Pa，最終澆注歷程總時長為2min。

實際澆注歷程是模擬澆注時長的1.66倍，分析原因認為： (1)藥漿非常黏稠，在流動過程中會在料斗和套筒等形成掛壁現(xiàn)象，掛壁越久越厚，而且在實際澆注過程中藥漿會阻塞花板孔，延遲了總的澆注時間；(2)由于實驗過程和數(shù)值模擬二者所選用花板規(guī)格的差異性，導致總澆注時間的差異。

綜上，誤差在一定的范圍內，可以認定該模擬方法一定程度上正確反映了藥漿的流動規(guī)律。

3.3 藥柱質量檢測

澆注工序結束后，將模擬實驗器進行固化，并立即進行保溫。檢查其外表面無殘藥，并按照裝配過程的逆順序進行脫模。將藥柱分解為Φ20mm×20mm的長條藥柱，如圖11所示。分解過程中，未發(fā)現(xiàn)藥柱存在較大孔洞等質量問題。

圖11 固體藥柱的截面 Fig.11 Cross section of solid grain

4 結 論

(1)花板孔徑越大，藥漿流速越快。擠出脹大比隨孔徑長徑比(L/D)減小而增大，當孔徑長度一致時，孔徑越大，擠出脹大比越大。

(2)推進劑藥漿為典型的屈服假塑性流體，一定剪切速率范圍內，黏度隨剪切速率的增大而降低，且當流經花板和翼片時，流體剪切速率最大，黏度最低。藥漿黏度和剪切速率在流動過程中發(fā)生變化。

(3)目前固體推進劑藥漿的本構模型無法精確擬合，未來需補充實驗數(shù)據(jù)，將本構模型精確化，如引入彈性項，建立線性黏彈性模型或微分形式和積分形式的非線性黏彈性模型，從而更好地描述藥漿的流動狀態(tài)。