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

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

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

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

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

1 射流閥基本原理

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

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

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

2 射流閥開閉過程數(shù)值模擬

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

2.1 物理模型

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

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

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

圖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 數(shù)值計算方法

2.3.1 空間導數(shù)求解

基于無網(wǎng)格方法進行流場計算首要工作便是確定空間導數(shù)的擬合方法,這也是無網(wǎng)格算法能否成功的基礎。本文所采用的空間導數(shù)求解方法是基于Taylor級數(shù)展開而進行的?？梢院唵蚊枋鰹?首先對中心點處的基本變量結合衛(wèi)星點處的具體變量值進行泰勒展開,利用二次平方極小曲面逼近求解空間導數(shù),進而利用最小二乘求導方法計算數(shù)值通量。

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

(2)

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

f(O(h2))=

(3)

(4)

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

(5)

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

(6)

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

(7)

(8)

2.3.2 數(shù)值通量求解

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

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

(9)

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

2.4 計算條件

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

2.5 計算結果分析

2.5.1 不帶閥芯裝置的數(shù)值模擬

1)穩(wěn)態(tài)計算結果。

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

圖4 穩(wěn)態(tài)計算結果(馬赫數(shù)云圖)

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

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

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

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

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

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

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

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

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

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

2.5.2 帶閥芯裝置的數(shù)值模擬

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

圖7 帶閥芯流場圖

3 結束語

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