射流閥門開/閉過程無網格法數值模擬

王園丁,郭曼麗,譚俊杰,任登鳳

(1.上海空間推進研究所,上海 201112;2.上海空間發動機工程技術研究中心,上海 201112;3.南京理工大學 能源與動力工程學院,江蘇 南京 210094)

姿軌控發動機內流場存在著非常復雜的非定常流動現象,尤其當各種射流閥門開/閉時其流場更加復雜。射流閥開/閉過程與射流控制技術密不可分。射流控制理論由美國Harry Diamond實驗室于20世紀五十年代提出[1],并已經在工業工程中得到了較為廣泛的應用。例如,在石油工業中用于鉆采的射流式沖擊錘[2]、兵器工業中用于改變導彈姿態的雙穩換向閥[3]以及用于化工工業中電解鋁溫控系統的射流控制系統等[4]。射流閥通過控制端射流與主流的相互作用改變主流的流動特性,使主流流動方向發生偏轉,并在某一個方向產生一定的推力,以改變受控體的姿態和軌跡,實現對受控體的姿態控制。

國內外在射流控制方面也做了大量研究工作,大部分工作都是圍繞較基礎的數值模擬展開的。國外約翰霍普金斯大學應用物理實驗室的Roger等[5]、韓國航空航天大學的Heo[6]對射流閥進行了細致研究。北京理工大學是國內最早開展燃氣射流閥研制的高校,并成功研制了兵器系統某遠程火箭彈修偏用燃氣射流閥。文獻[7-9]通過數值模擬手段對啟動射流閥的舵機壓力特性、驅動機構的瞬態響應、噴口遮擋的影響等一系列問題進行了研究。文獻[3]對某遠程多管火箭彈姿控發動機雙穩閥超音速射流流場進行數值模擬,基于A100遠程多管火箭彈BK183-5姿控發動機50閥模型以及改進的42.5閥模型,采用隱式算法以及先進的V2F湍流模型對射流閥流場進行求解,通過與已有實驗數據的對比,合理解釋了實驗中發生的現象,驗證了利用數值模擬手段求解超音速射流閥的可靠性和準確性,取得了較為滿意的結果。

與傳統的網格方法相比,無網格方法不需要通過離散節點組成的網格單元進行求解區域的劃分,而是通過“點云”結構代替常規的網格單元離散計算域,布點較為靈活,因而國內外開展了大量的研究工作[10-17],主要涉及基礎算法研究,如將迎風格式引入無網格算法、求解各項異性點云結構的無網格算法等,但鮮有應用于姿軌控發動機射流閥流場的數值模擬。本文采用耦合Wilcoxk-ω湍流模型的最小二乘無網格,對含拓撲結構改變的姿軌控發動機超音速射流閥流場進行數值模擬,成功模擬出附壁流動、主流的切換等復雜的非定常流動狀態。

1 射流閥基本原理

圖1為射流換向閥結構示意圖[5],主要由主流入口、控制通道、尖劈和排氣通道組成,主流通道結構造型為拉瓦爾噴管。圖中,b為主流通道寬度,d為噴管厚度,φ為擴張角。

射流換向閥的基本工作過程可以根據圖1描述如下:主流從主流入口處進入漏斗形噴管后,經垂直段通道到達喉部,若打開左側控制端的電磁閥,則控制流從左側控制通道進入噴管,并在喉部對主流流動產生向右的作用力,主流流動將發生向右偏轉,在擴張段形成右側附壁流動,并經尖劈分流后從右側排氣通道流出。圖1中尖劈結構起到限制分流的作用。若要使得主流發生向左偏轉,則只需要關閉左側控制端電磁閥,并打開右側電磁閥使控制流從右控制端流入噴管,則控制流會在喉部對主流產生向左的作用力而在擴張段形成左側附壁流,并沿左側排氣通道流出。

圖1 射流換向閥結構示意圖

2 射流閥開閉過程數值模擬

采用耦合Wilcoxk-ω湍流模型的最小二乘無網格,對含拓撲結構改變的姿軌控發動機超音速射流閥流場進行數值模擬。

2.1 物理模型

圖2為不帶閥芯裝置的計算模型,射流口高度為1.8 mm,喉部主流入口寬度為1.8 mm,喉部出口寬度(擴張段入口)為2.1 mm,擴張半角為12°,z方向為5.0 mm(將z方向中間位置定義為z=0坐標起點)。對此模型進行數值模擬,可以對控制流與主流的相互作用進行初步研究,以觀測在控制流作用下的主流切換效果。

圖2 不帶閥芯裝置的計算模型(z=0截面模型尺寸)

圖3為帶閥芯裝置的計算模型。主流流入射流閥后,控制流在喉部對主流作用,使得主流發生偏轉,并對閥芯產生作用力驅動閥芯運動到另一側,從而起到密封的作用。

圖3 帶閥芯裝置的計算模型

2.2 控制方程

可壓縮雷諾平均N-S控制方程可表示為如下守恒形式:

(1)

式中:Q為守恒變矢量;E,F,G分別為x,y,z坐標方向的無粘通量;Ev,Fv,Gv分別為x,y,z坐標方向的粘性通量。

Q=(ρuρvρwρe)T
E=(ρuρu2+pρuvρuw(ρe+p)u)T
F=(ρvρuvρv2+pρvw(ρe+p)v)T
G=(ρwρuwρvwρw2+p(ρe+p)w)T
Ev=(0τxxτxyτxzuτxx+vτxy+wτxz-qx)T
Fv=(0τxyτyyτyzuτxy+vτyy+wτyz-qy)T
Gv=(0τxzτyzτzuτxx+vτyz+wτzz-qz)T

ρ,p,e分別為流體密度、壓強、單位質量的氣體總內能;u,v,w分別為x,y,z方向的速度分量;γ為比熱比,對于空氣γ=1.4;τ為應力項;q為熱通量項。具體表達式見文獻[12]。

2.3 數值計算方法

2.3.1 空間導數求解

基于無網格方法進行流場計算首要工作便是確定空間導數的擬合方法,這也是無網格算法能否成功的基礎。本文所采用的空間導數求解方法是基于Taylor級數展開而進行的。可以簡單描述為:首先對中心點處的基本變量結合衛星點處的具體變量值進行泰勒展開,利用二次平方極小曲面逼近求解空間導數,進而利用最小二乘求導方法計算數值通量。

假設流場中任意離散點i處的流場值可以表示為函數f(x,y,z),且設函數f(x,y,z)在點云Ci內n階連續可導,對點云Ci內任意一個衛星點s處的流場值f(xs,ys,zs),可以通過Taylor級數展開得到,即f(xs,ys,zs)可表示為

(2)

式中:i表示中心點,s表示以點i為中心點的點云結構中的衛星點,s∈Ci。由點云結構的構造方法可知,衛星點個數大于3,而式(2)中只有3個未知數,所以式(2)是一個矛盾方程組,此時2個未知的偏導數可以通過最小二乘法擬合得到。最小二乘法需要使得擬合曲線盡可能地逼近物理解,即需要滿足監測數據與擬合曲線之間的偏差(亦可理解為式(2)中的高階小量)的平方和最小:

f(O(h2))=

(3)

(4)

顯式求解式(4),可將空間導數表示成如下形式:

(5)

方程組系數αis,βis和γis的具體表達形式為

(6)

式(6)中的系數具體表達式為

(7)

(8)

2.3.2 數值通量求解

本文借鑒文獻[10]的思想,在點云結構中心點與衛星點連線的中點處構造一個虛擬界面J,將MUSCL(monotone upstream-centred schemes for conservation laws,MUSCL)格式引入無網格方法中,用以重構該界面左、右兩側的原始變量值,進而虛擬界面處的無粘通量可通過迎風格式進行求解。本文采用的迎風格式為AUSM+-up格式。

對于粘性通量中的速度和溫度的梯度可以直接應用式(5)求解。為了避免因出現奇偶失聯導致的計算不準或者錯誤的情況,對虛擬中點處的梯度值進行修正:

(9)

通過上述方法對控制方程進行離散之后,控制方程將變成常微分方程,本文采用三階SSP型Runge-Kutta法,并使用當地時間步長進行時間推進求解離散后的控制方程。湍流模型采用Wilcoxk-ω兩方程模型,針對各項異性點云結構,采用經過驗證的點云重構技術[15]進行病態點云結構的修正,以提高模擬粘性流動的準確性。

2.4 計算條件

針對上述2種計算物理模型,設置主流入口壓力為10 MPa,溫度300 K;左右控制端控制流入口壓力為6.0 MPa,溫度300 K;不帶閥芯裝置的計算模型,出口壓力為3.5 MPa,溫度300 K,帶閥芯裝置下部左右兩側出口壓力為1個大氣壓,溫度300 K。不帶閥芯裝置的計算模型,整個計算域離散點分布為90萬;帶閥芯裝置的整個計算域離散點分布為286萬,其中閥芯周圍離散點為56萬。先打開主流,將流場計算穩定,再打開一側控制端使控制流進入流場,進行流場計算。

2.5 計算結果分析

2.5.1 不帶閥芯裝置的數值模擬

1)穩態計算結果。

圖4給出了不帶閥芯裝置模型穩態計算過程中馬赫數的變化情況。此時左右通道控制端做固壁處理。

圖4 穩態計算結果(馬赫數云圖)

從圖4中可以看出,隨著主流的流入,開始會形成左右對稱的流場結構,如圖4(a)所示,隨著時間的推進,根據射流控制理論及主流自身紊流特性,主流會向右側發生偏轉,如圖4(b)所示,并最終形成穩定的右側附壁流動。此時左右控制端的壓力為4.65 MPa左右。在穩態計算結果的基礎之上,如圖4(c)所示,打開右端控制閥進行非穩態計算。

2)右側控制端打開的計算結果。

圖5為在穩態基礎之上打開右端控制閥后,計算得到的馬赫數隨時間變化的云圖,此時左通道控制端做固壁處理。

圖5 右側控制閥打開后馬赫數變化歷程

從圖5中可以看出,在t=0.4 ms時,主流便發生了明顯偏轉;在t=0.8 ms時,主流已經在下部經過分流尖劈從右出口逐漸轉向左出口;在t=1.2 ms時,大部分的主流已經經由下部左側出口流出;t=1.6 ms時,主流已經實現從右側到左側的切換,只是流場還沒有最終穩定下來;在2.5 ms時氣流實現了從右側附壁流到左側附壁流的完全切換。完全切換后的左通道控制端附近的壓力為5.6 MPa,右通道控制端附近的壓力為6.0 MPa。在此基礎之上,打開左通道控制端,關閉右通道控制端(做固壁處理),進行二次切換的非穩態計算。

3)左側控制端打開的計算結果。

在氣流完全切換至穩定的左側附壁流動后,關閉右側控制端,打開左側控制端,對流場進行計算,所得的馬赫數隨時間變化的云圖如圖6所示。

圖6 左側控制閥打開后馬赫數變化歷程

從圖6中可以看出,在6.0 MPa的控制壓力下,主流也實現了偏轉,但因為此時左右通道壓差較小,所以需要更長的時間實現主流的完全切換。主流實現完全切換的時間為6.5 ms。

對不帶閥芯裝置的模型進行計算,在右通道控制端6.0 MPa壓力的作用下,主流從右側附壁流在2.5 ms內成功切換至左側附壁流。在此基礎之上,關閉右通道控制端,打開左通道控制端讓6.0 MPa控制流對主流進行作用,成功計算出了二次切換的流場結構圖,所需要的時間為6.5 ms。說明在現有計算模型和計算條件下可以實現射流閥的雙向切換。

2.5.2 帶閥芯裝置的數值模擬

針對帶閥芯裝置的模型,設置主流入口壓力為10.0 MPa,將主流控制端閥門打開,關閉左右控制端(做固壁處理),計算流場至穩定狀態,速度云圖如圖7所示。由于閥芯和閥體之間間隙小,在最小的喉道位置處氣體流動速度將為音速。由于閥芯與閥體完全接觸,左側排氣通道被完全關閉,氣流無法通過左側排氣通道流出,只能通過右側排氣通道流出,此時左右通道端部壓力為4.35 MPa,因此要實現左右控制端氣流順利進入通道驅動主流發生偏轉,控制端氣流需大于4.35 MPa。

圖7 帶閥芯流場圖

3 結束語

采用耦合Wilcoxk-ω湍流模型的最小二乘無網格,對含拓撲結構改變的姿軌控發動機超音速射流閥流場進行數值模擬。成功模擬出附壁流動、主流的切換等復雜的非定常流動狀態,捕捉到了較為清晰的流場結構與流動特征,能夠客觀地反應射流閥的基本切換過程。數值模擬結果為姿軌控發動機射流閥門的設計提供了參考,也擴大了耦合兩方程湍流模型的最小二乘無網格方法數值模擬的范圍。