橫向過載下考慮顆粒聚合行為的顆粒阻尼特性研究①

楊文婧,吳學婷,劉貝貝,劉佩進

(1.西北工業大學 航天學院,燃燒、熱結構與內流場重點實驗室,西安 710072;2.中國空空導彈研究院,洛陽 471009)

0 引言

近年來由于導彈突防能力需求的不斷提升,高機動性能下的固體火箭發動機需要承受復雜的過載狀況,因此需要提升發動機工作時的安全性和可靠性。通常為了提高發動機的比沖、抑制燃燒不穩定性,同時增加裝藥的能量,需要采用含鋁復合推進劑,燃燒后產生的大量Al2O3凝相顆粒,形成復雜的發動機兩相流系統。發動機中凝相顆粒在過載力作用下,會發生偏移聚集行為,增加了顆粒之間的碰撞概率,更容易產生大粒徑的顆粒。固體火箭發動機中存在多種增益與阻尼因素,主要的增益因素包括推進劑的燃燒增益、聲渦耦合,主要的阻尼因素包括粒子阻尼、噴管阻尼等,是否會出現不穩定現象并維持壓強振蕩是各種因素綜合作用的結果。

王國輝等[1]對有過載和無過載多種工況下發動機內流場進行了計算,得出了氣相及顆粒相在過載條件下的流動特性,氣相受橫向過載的影響不大,而凝相顆粒的粒徑越大,受過載的影響越大。SABNIS等[2]認為粒徑在50 μm以下的顆粒受氣相流場的影響大于橫向過載的影響,50 μm以上受過載的影響較大,在過載力作用下其運動軌跡發生改變,向發動機承載一側偏移,加劇凝相粒子之間的碰撞以及凝相粒子與壁面之間的碰撞。田維平等[3]總結了粒子間碰撞與粒壁碰撞的規律,其中,不同粒徑的液滴相互碰撞的結果可能為反彈、聚合及破碎等,導致液滴粒徑分布發生變化;而液滴碰撞壁面機制一般依據韋泊數(We)大小,分為吸附、反彈、展成液膜及飛濺等四種結果,在燃燒室中產生的Al2O3液滴與壁面之間主要表現為吸附與反彈共存;而在噴管收斂段的主要表現為反彈。武淵等[4]認為由于過載對凝相顆粒分布的影響導致燃燒室內的溫度場發生畸變,溫度最高點發生在藥柱后段與絕熱層連接部位,該處流場的溫度達到了4000 K以上。許團委等[5]采用歐拉-拉格朗日模型對過載下戰術發動機的內流場進行了數值仿真,研究表明凝相顆粒受橫向過載的影響易在發動機承載面的絕熱層表面沿著流場方向形成一條粒子濃度緩慢增大的聚集帶。苗琳[6]計算了多種過載形式下發動機內的兩相流場,其中,橫向過載以及軸向過載易導致凝相顆粒在發動機承載面以及后封頭處沉積,當過載較小時,適當的轉速可減少承載面和后封頭處的粒子濃度。劉長猛等[7]的研究表明隨著軸向過載的增大,凝相粒子在發動機內部的駐留時間會減小。劉中兵等[8]認為在橫向長時間中小過載下,在凸起較大的人脫根部和后接頭拐彎處因凝相粒子流聚集會出現較大的燒蝕坑,燒蝕極不均勻,易發生燒穿故障,應在結構型面設計中避免較大的臺階出現。WANG 等[9]的研究表明隨著橫向過載的增大,噴管喉部處的顆粒濃度幾乎呈線性增加,導致噴管喉道的異常燒蝕急劇惡化。王立武等[10]對典型喉襯材料抗高橫向過載燒蝕性能進行了研究,研究認為因顆粒沖刷引起的機械侵蝕的加劇是導致喉部偏燒蝕的主要原因,在高橫向過載下,軸編C/C抗過載燒蝕性能最好。張翔宇等[11]通過建立固體發動機火箭橇地面過載模擬試驗方法,復現了飛行過載誘發的發動機不穩定燃燒現象。之后,張翔宇等[12]通過對故障發動機開展試驗研究,驗證了導彈飛行過載是引起發動機不穩定燃燒的最主要原因。橫向過載會使顆粒的濃度發生改變,承載一側聚集大量的大粒徑顆粒,而另一側只有一些隨流性好的小顆粒。檀葉等[13]利用不穩定燃燒線性理論分析認為過載引起的顆粒濃度分布變化是引起不穩定燃燒的關鍵因素。游艷峰等[14]通過離散元法得出在橫向過載條件下,顆粒的非均勻分布式燃燒會使粒子阻尼減小,對發動機的穩定性十分不利。在發動機中,高濃度區域粒子阻尼增大,而低濃度區域粒子阻尼減小。30 μm以內的粒子在過載下粒子阻尼基本不變;較大的粒子在大橫向過載下,極短時間(5 ms)后粒子阻尼開始減小,且粒徑越大的粒子其阻尼減小得越早、越快。劉佩進等[15]總結了固體火箭發動機中燃燒不穩定現象的研究進展,提出了橫向過載會改變發動機內凝相顆粒的分布,影響顆粒的阻尼效應,同時顆粒空間分布的改變也會使燃燒室內鋁粒子燃燒釋熱區域發生變化,易導致熱聲不穩定。

橫向過載會改變固體火箭發動機內凝相粒子原有的運動軌跡,增大顆粒與顆粒,顆粒與壁面的碰撞概率,改變凝相粒子的粒徑和濃度分布,使顆粒的阻尼特性發生變化,造成發動機的燃燒不穩定。故而,研究橫向過載下Al2O3的聚合行為以及其對顆粒阻尼的影響具有很大的意義。

1 數值模型

本文考慮徹體力作用下燃燒室凝相顆粒運動過程,分析并獲得徹體力對凝相顆粒運動規律、相互作用規律和空間動態分布的影響,因此選擇燃燒室作為計算區域。考慮顆粒數量會很大程度地影響計算量,采用一定比例縮小的圓柱形幾何模型來代替發動機燃燒室模型,模型長度330 mm,直徑55 mm,如圖1所示。

圖1 計算幾何模型圖

1.1 數值方法

1.1.1 氣相控制方程

氣相控制方程是基于連續介質流場進行求解計算,主要的控制方程是連續方程和動量方程,見式(1)、式(2)所示。

連續方程:

(1)

動量方程:

(2)

1.1.2 顆粒相控制方程

以單個顆粒為研究對象進行分析,則其運動方式一般情況下可以分為平動和轉動。離散元模型中,單個顆粒的運動主要由牛頓第二定律控制進行的,見式(3)、式(4)所示。

平動方程:

(3)

轉動方程:

(4)

1.1.3 氣固耦合模型

在固體發動機中,凝相顆粒將在燃氣的拖曳作用下運動,此時氣相也會受到顆粒相的反作用力,兩者之間將會出現動量傳遞,兩相流中顆粒相與氣相的雙向耦合主要通過動量傳遞實現。氣固兩相流耦合計算流程示意圖如圖2所示。

圖2 氣固兩相流耦合計算流程示意圖

本文所采用的CFD-DEM氣固耦合計算流程如下:

(1)給定一個初始的氣相流場,待氣相流場收斂后,給定初始顆粒位置、粒徑等信息,然后在一個時間步長內,計算顆粒與顆粒、顆粒與壁面的作用力、流體曳力以及徹體力等作用力,基于DEM方法計算出顆粒新的位置、速度等信息;

(2)將DEM方法所得顆粒位置、粒徑及顆粒數等信息應用于流體網格內計算出流體網格的孔隙率以及流體的曳力,在孔隙率和流體曳力已知情況下,再次利用N-S方程計算新的流場信息,待流場收斂后重新迭代計算顆粒相信息;

(3)交替求解顆粒相與氣相的控制方程,求得兩相流中顆粒相的運動過程。

1.2 粒子碰撞聚合模型

通過對飛行試驗以及地面模擬過載試驗中發動機內顆粒粒度的測量,發現發動機顆粒在過載狀態下粒徑值大于地面靜止環境時所測值,產生大粒徑顆粒的概率顯著增加,因此很有必要考慮過載下顆粒之間相互碰撞引起顆粒聚合行為對發動機中凝相顆粒粒徑改變的影響以及獲取此時顆粒的運動規律。

針對液滴碰撞規律相關學者已開展了大量的理論、實驗及數值模擬研究,然而不同的研究對象其碰撞規律卻不盡相同,尤其高表面張力、高粘度的Al2O3顆粒之間的碰撞,其碰撞規律和水、煤油之間有較大的差別。Al2O3凝相顆粒在固體火箭發動機中處于高溫高壓環境下,因此很難針對此狀態下的顆粒展開顆粒間碰撞聚合實驗研究。西北工業大學夏盛勇等[16]基于VOF方法、利用開源程序Gerris建立了Al2O3顆粒之間的碰撞模型,為本文建立顆粒間碰撞聚合模型提供了數值研究基礎。

本文中顆粒間碰撞等效圖,如圖3所示。

圖3 顆粒碰撞等效圖

本文中顆粒間的碰撞規律主要由韋伯數和碰撞參數控制,見式(5)和式(6)所示。

韋伯數:

(5)

碰撞參數:

(6)

根據文獻[7]可知,隨著韋伯數和碰撞參數的變化,顆粒與顆粒之間碰撞后會發生微小變形后聚合、反彈、大變形后聚合、自反分離以及拉伸分離等運動行為,對應的判別公式見式(7)～式(11)所示。

微小變形后聚合:

(7)

(8)

反彈:

(9)

自反分離:

(10)

拉伸分離:

(11)

大變形后聚合:介于其他區域之間。

由式(7)～式(11)可知,不同的韋伯數、碰撞參數下,顆粒在碰撞后將會有不同的運動行為。上述碰撞結果又可分為兩部分進行研究:第一部分是微小變形后聚合行為,此情況發生在韋伯數很低的情況下,當給定顆粒直徑、顆粒密度以及顆粒表面張力等參數時,韋伯數僅與顆粒間相對速度有關,因此微小變形后聚合行為主要發生在顆粒間相對速度很低的情況;第二部分包括反彈、大變形后聚合、自反分離和拉伸分離等行為,顆粒間相對速度較高。

根據式(7)～式(11),可得本文顆粒碰撞模型圖。由于上述兩種情況對應的韋伯數差別較大,這里分別給出對應的碰撞模型圖,如圖4所示。由圖4可知,Al2O3顆粒碰撞后發生的微小變形后聚合行為主要發生在韋伯數小于0.2的情況下,當顆粒間相對速度很低時,便會出現微小變形后聚合行為;隨著韋伯數的增大,會相繼出現反彈、大變形后聚合、自反分離以及拉伸分離等行為,此時顆粒間相對速度會較大一些。

(a)Small Weber numbers (b)Large Weber numbers

由式(7)～式(11)可知,隨著顆粒碰撞后韋伯數和碰撞參數的不同,顆粒碰撞后的運動行為可分為微小變形后聚合、反彈、大變形后聚合、拉伸分離以及自反分離等行為;由式(5)、式(6)可知,韋伯數和碰撞參數主要與顆粒的密度、顆粒間相對速度、顆粒直徑、顆粒的表面張力以及顆粒碰撞時的角度有關。

綜上可得粒子碰撞模型的嵌入方法:

(1)由氣相程序計算一個初步的流場,凝相顆粒在氣相曳力作用下移動,由其他已知信息可得此狀態下每個顆粒在給定時間步下實時的位置、線速度、旋轉角速度、顆粒粒徑等信息。

(2)由顆粒的位置信息來判斷顆粒間接觸狀態。當兩個顆粒中心間距大于等于兩個顆粒半徑之和時,顆粒跳出碰撞模型計算模塊,保持原來的狀態進入下一步計算模塊;當兩個顆粒中心間距小于兩個顆粒半徑之和時,進入顆粒碰撞模型計算模塊。

(3)當兩個顆粒中心間距小于兩顆粒半徑之和時,首先根據兩個顆粒的位置信息和速度信息計算出顆粒間的碰撞參數B,然后由顆粒的密度、相對速度、顆粒直徑以及表面張力等信息計算出顆粒的韋伯數We,接著由式(7)～式(11)計算該碰撞參數B下對應的臨界微小變形后聚合韋伯數、臨界反彈韋伯數、臨界大變形后聚合韋伯數、臨界自反分離韋伯數以及臨界拉伸分離韋伯數等,通過比較該碰撞狀態下的韋伯數與臨界韋伯數確定顆粒碰撞后的運動行為。

(4)當顆粒的韋伯數在圖4(a)所示曲線的左側時,顆粒碰撞后發生微小變形后聚合行為;當顆粒的韋伯數處于圖4(b)所示的1區域時,顆粒碰撞后發生反彈行為;當顆粒的韋伯數處于圖4(b)的2區域時,顆粒碰撞后發生大變形后聚合行為;當顆粒的韋伯數處于圖4(b)的3區域時,顆粒碰撞后發生自反分離行為;當顆粒的韋伯數處于圖4(b)的4區域時,顆粒碰撞后發生拉伸分離行為。由于固體火箭發動機燃燒室中凝相顆粒間的相對速度較低,顆粒的韋伯數也比較低,因此顆粒碰撞后運動行為基本為微小變形后聚合和反彈行為。

(5)當顆粒碰撞后發生微小變形后聚合行為時,本文擬采用“兩個顆粒碰撞后聚合為一個顆粒”的方法。由體積守恒定律、質量守恒定律和動量守恒定律計算出新的顆粒的位置、直徑及速度等信息。

(6)當顆粒碰撞后發生反彈行為時,根據顆粒運動控制方程計算出顆粒碰撞后的物性參數。

(7)循環計算整個計算域內的顆粒,依次判斷每兩個顆粒間的碰撞狀態,根據碰撞時顆粒的韋伯數和碰撞參數確定顆粒碰撞后的運動行為,最后可得該時間步下所有顆粒的碰撞運動結果。

(8)最后迭代計算整個程序,直至要求的時間步下結束程序。

當固體火箭發動機處于高過載狀態時,顆粒會由于過載的作用而發生橫向偏移,從而增加了顆粒的空間濃度,對應的顆粒間碰撞概率會顯著提高,此時不同顆粒碰撞時對應的韋伯數和碰撞參數不同,因此會有不同的碰撞結果。本文針對徹體力作用的固體火箭發動機燃燒室展開顆粒碰撞數值模擬研究,擬得到考慮顆粒碰撞聚合的兩相流模型以及對應的凝相顆粒運動規律。

1.3 粒子阻尼算法

在常規狀態下時,發動機中顆粒處于均勻分布狀態,此時發動機聲場為為小擾動,可利用Culick微粒松弛理論計算粒子阻尼大小,即通過壓力振幅的時間衰減常數αp來表征粒子阻尼大小,見式(12)～式(14)所示。

(12)

(13)

(14)

從式(12)～式(14)可知,影響顆粒阻尼值的因素可分為兩部分:第一部分與粒子空間分布有關,即單位體積中粒子質量分數越大時,顆粒阻尼值越大;第二部分與顆粒相和氣相的物性參數有關,主要包括粒子直徑、粒子密度、聲振蕩頻率及氣相動力粘度等參數。

當固體火箭在作機動飛行時,其發動機會受到飛行過載的作用,此時發動機燃燒室內凝相顆粒會發生偏移聚集行為,一方面改變了凝相顆粒在燃燒室的空間分布狀態,形成局部的低濃度區域及高濃度區域,從而通過顆粒的局部質量分數改變了顆粒的阻尼值;另一方面,當考慮顆粒碰撞聚合行為時,顆粒的粒徑會發生改變,微粒的動力松弛時間和熱松弛時間也都會改變,從而過載通過影響顆粒的粒徑來改變顆粒的阻尼值。因此,過載下不同局部區域顆粒的阻尼值不一致,通過整個發動機內顆粒的平均質量分數及顆粒的初始平均粒徑值來計算顆粒阻尼值是不準確的,且不能表征不同粒徑下顆粒空間分布不均勻性引起的阻尼值的改變。

針對上述兩方面的影響,本文擬通過引入與空間分布相關的權重函數來對粒子阻尼值進行,即粒子阻尼值表征顆粒空間分布不均勻情況以及發動機凝相顆粒聚合引起粒徑變化情況下凝相顆粒的阻尼值。基本思路是:先將計算域劃分為n個小區域,然后計算每個區域內顆粒的阻尼值,最后對每個區域進行加權平均,從而得到粒子阻尼值。

對于給定的網格i,該區域內顆粒數為m,該區域顆粒的阻尼值為αpi,見式(15)所示:

(15)

根據網格大小對每個網格的顆粒阻尼值αpi進行加權平均,從而得到燃燒室內粒子阻尼值αpn,見式(16)所示:

(16)

由此可得考慮顆粒粒徑變化及空間分布不均勻情況下顆粒的粒子阻尼值αpn,見式(17)所示:

(17)

2 結果與分析

2.1 網格無關性驗證

為了檢驗網格數對粒子阻尼值的影響,本文對橫向過載30g的顆粒運動結果進行后處理計算,得到網格數分別為26 480、87 840、206 800、401 200網格數下粒子阻尼值,如圖5所示。由圖5可知,不同網格數下顆粒的粒子阻尼值隨時間的變化趨勢基本一致,因此考慮顆粒空間分布及粒徑變化的粒子阻尼值可表征凝相顆粒在過載下的阻尼值變化情況。

圖5 網格數對粒子阻尼值的影響

2.2 粒子配位數

橫向過載會引起顆粒的橫向偏移,而不同橫向過載對顆粒橫向偏移的影響程度是不同的,所以顆粒的空間位置分布也不一致,因此給出不同橫向過載下粒子配位數圖,如圖6所示。由圖6可知,當橫向過載為10g時,顆粒沿X向偏移較少,其在遠離壁面區域的粒子配位數基本都低于26,而在近壁面區域中,顆粒的粒子配位數增大,達到60左右的顆粒數很少,且集中在中心區域;當橫向過載為20g時,與10g的橫向過載偏移相比,其偏移略高,顆粒在遠離壁面區域的粒子配位數變化不大,均低于26,而在近壁面區域中,粒子配位數為60左右的顆粒數增大,且沿著過載壁面中心區域向兩邊擴散分布;當橫向過載為30g時,顆粒沿X向偏移較多,其在遠離壁面區域的粒子配位數仍在26以下,而在近壁面區域中,粒子配位數為60左右的顆粒數進一步增大,且沿著整個過載面離散分布;當橫向過載為40g時,顆粒沿X向偏移最多,其在遠離壁面區域的粒子配位數仍舊低于26,而在近壁面區域中,粒子配位數為60左右的顆粒數與30g的橫向過載情況相比變化不大,且沿著整個過載面離散分布。

(a)Lateral overload=10g (b)Lateral overload=20g

因此,隨著橫向過載的增大,顆粒的橫向偏移會越來越大,在遠離壁面區域顆粒的粒子配位數基本一致,在近壁面區域顆粒的粒子配位數會顯著增大。故橫向過載是影響近壁面區域顆粒聚集程度的重要因素。

2.3 顆粒間碰撞概率

不同橫向過載下,凝相顆粒所受的徹體力不同,導致顆粒的橫向偏移量有較大差別,顆粒的空間聚集程度也不同,所以對應的顆粒間碰撞概率會有很大的差異性。不同橫向過載下顆粒間碰撞概率隨時間的變化過程圖如圖7所示。

由圖7可知,當橫向過載為10g時,6 ms時刻后發生顆粒碰撞,大部分時刻下顆粒間的碰撞概率約為 0.000 065,少許時刻下顆粒間碰撞概率為0.000 13左右,顆粒間的碰撞概率基本不隨時間變化;當橫向過載為20g時,3.6 ms時刻后發生顆粒碰撞,大部分時刻下顆粒間的碰撞概率在0.000 065～0.000 073之間,部分時刻下顆粒間的碰撞概率在0.000 13～0.000 15之間,少部分時刻下顆粒間的碰撞概率在0.000 2左右,顆粒間碰撞概率會隨著時間略微增加;當橫向過載為30g時,3.6 ms時刻后發生顆粒碰撞,大部分時刻下顆粒間的碰撞概率在0.000 065～0.000 11之間,部分時刻下顆粒間的碰撞概率在0.000 13～0.000 18之間,少部分時刻下顆粒間的碰撞概率為0.000 2～0.000 3之間,顆粒間碰撞概率隨時間呈現遞增的趨勢;當橫向過載為40g時,2.4 ms時刻后發生顆粒碰撞,大部分時刻下顆粒間的碰撞概率在0.000 065～0.000 11之間,部分時刻下顆粒間的碰撞概率在 0.000 13～0.000 2之間,少部分時刻下顆粒間的碰撞概率在0.000 2～0.000 34之間,且顆粒間碰撞概率隨時間遞增趨勢更加明顯。

(a)Lateral overload=10g (b)Lateral overload=20g

當橫向過載較小時,顆粒群的聚集程度較小,導致顆粒間的碰撞概率較小,且碰撞概率隨時間變化較慢;隨著橫向過載逐漸增大,顆粒群的聚集程度增大,導致顆粒間的碰撞概率明顯增大,且碰撞概率隨時間增加趨勢更加明顯,故橫向過載大小是影響顆粒間碰撞概率的重要因素。

2.4 粒子阻尼

顆粒空間分布狀態的改變及顆粒粒徑的改變皆可以引起顆粒阻尼值的改變,因此給出不同橫向過載下顆粒的粒子阻尼系數及粒子阻尼值,如圖8所示。

由圖8可知,不同過載下初始時刻顆粒阻尼系數為1.0,粒子阻尼值為6.98,阻尼系數和粒子阻尼值都隨著時間逐漸降低,且變化趨勢一致。當橫向過載為10g時,顆粒阻尼系數和粒子阻尼值隨時間變化很小;當橫向過載為20g時,15 ms時刻粒子的阻尼系數為0.9、粒子阻尼值為6.3,阻尼系數和粒子阻尼值變化率都是10%;當橫向過載為30g時,15 ms時刻粒子的阻尼系數為0.78、粒子阻尼值為5.4,阻尼系數和粒子阻尼值變化率都是22%;當橫向過載為40g時,15 ms時刻粒子的阻尼系數為0.63、粒子阻尼值為4.4,阻尼系數和粒子阻尼值變化率都是37%。

圖8 粒子阻尼值

由此可知,在相同橫向過載下,顆粒的阻尼系數和粒子阻尼值隨時間變化規律是一致的,而不同橫向過載下,由于顆粒的空間分布狀態以及粒徑變化規律有很大的差異,所以顆粒的阻尼系數和粒子阻尼值變化規律明顯不同,且隨著橫向過載的增大,阻尼系數和粒子阻尼值降低趨勢更加明顯。因此,向過載是影響顆粒阻尼系數和粒子阻尼值的重要因素。

3 結論

橫向過載下,顆粒橫向偏移顯著增加,顆粒間的碰撞概率增大,并且碰撞概率隨過載作用時間增加大幅度提高。顆粒間因碰撞聚合易生成的大粒徑顆粒,大粒徑顆粒的產生使得顆粒空間分布更具不均勻性,進而導致顆粒的阻尼值明顯降低,最終影響發動機的燃燒穩定性。橫向過載的大小對顆粒行為及發動機性能的影響規律:

(1)橫向過載增大時,會顯著地增加顆粒的橫向偏移,且近壁面區域顆粒聚集程度顯著增加;

(2)橫向過載越大時,顆粒間的碰撞概率越大,且碰撞概率隨時間增加的趨勢更加明顯;

(3)綜合考慮顆粒橫向偏移、顆粒數量以及顆粒粒徑大小等因素,發現橫向過載越大時,顆粒阻尼降低程度越大,發動機燃燒不穩定性越嚴重。