謝信輝,柳煜瑋,劉耀峰
(中國(guó)航天空氣動(dòng)力技術(shù)研究院,北京 100074)
大推力噴流控制技術(shù)可實(shí)現(xiàn)大機(jī)動(dòng)、強(qiáng)突防[1],在多種飛行器上都有廣泛的應(yīng)用,如低速狀態(tài)下橫流中煙囪羽流[2]和高速狀態(tài)下飛行器攔截系統(tǒng)[3]。特別是在攔截任務(wù)的最后階段,要求攔截飛行器具有高度的機(jī)動(dòng)性,提供直接作用力、快速響應(yīng)的側(cè)向噴流控制推進(jìn)器,這是對(duì)傳統(tǒng)氣動(dòng)控制面的有效補(bǔ)充。
噴流產(chǎn)生的三維流場(chǎng)是高度復(fù)雜的,包含復(fù)雜的波系、渦系結(jié)構(gòu)及其相互干擾,呈現(xiàn)高度非定常、非線性、強(qiáng)湍流特性[4],現(xiàn)階段對(duì)運(yùn)動(dòng)條件下非定常效應(yīng)的研究很少。早在1964年Zukoski等[5]就開展了不同馬赫數(shù)下聲速噴管欠膨脹噴流噴入超聲速外流的實(shí)驗(yàn)研究,指出了弓形激波、分離激波、分離區(qū)等主要特征。隨著計(jì)算流體力學(xué)的發(fā)展,CFD技術(shù)在噴流干擾流場(chǎng)研究領(lǐng)域得到廣泛應(yīng)用。1991年張涵信等[6]利用NND格式,完成了鈍錐完全氣體側(cè)噴干擾流場(chǎng)的計(jì)算。2005年徐敏等[7]研究了脈沖發(fā)動(dòng)機(jī)噴流的非定常性對(duì)飛行器控制的關(guān)系。
當(dāng)噴流開啟與關(guān)閉時(shí),存在強(qiáng)烈的非定常特性,同時(shí)伴隨飛行器運(yùn)動(dòng)與姿態(tài)變化,存在氣動(dòng)特性與運(yùn)動(dòng)特性的強(qiáng)耦合,為了有效評(píng)估飛行性能,須要考慮氣動(dòng)運(yùn)動(dòng)耦合帶來(lái)的非定常問(wèn)題。關(guān)于空氣動(dòng)力學(xué)與飛行力學(xué)耦合方法的研究,國(guó)外近年來(lái)發(fā)展和應(yīng)用了大規(guī)模模擬的最新數(shù)值方法[8-11]。Sahu等[12]在超級(jí)計(jì)算機(jī)上對(duì)飛行器的“虛擬飛行”進(jìn)行了研究,使用耦合的CFD/RBD技術(shù)對(duì)飛行器的實(shí)際飛行軌跡和所有相關(guān)的非定常自由飛行空氣動(dòng)力學(xué)進(jìn)行數(shù)值預(yù)測(cè)。Silton等[10]通過(guò)與自由飛行實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較,完成了對(duì)虛擬飛行模擬計(jì)算結(jié)果的驗(yàn)證。劉耀峰等[1]應(yīng)用非定常方法數(shù)值模擬了側(cè)向噴流干擾流場(chǎng)建立與消退過(guò)程。陳堅(jiān)強(qiáng)等[13]模擬了舵面運(yùn)動(dòng)與側(cè)向噴流間耦合干擾問(wèn)題。以上研究具有豐富的價(jià)值,本文的重點(diǎn)是將氣動(dòng)--運(yùn)動(dòng)一體化方法應(yīng)用到噴流干擾問(wèn)題中,對(duì)運(yùn)動(dòng)條件對(duì)噴流干擾的影響做進(jìn)一步的探索研究,為工程實(shí)踐的應(yīng)用提供參考。
本文基于雷諾平均N-S方程和剛體六自由度運(yùn)動(dòng)方程,介紹了將飛行力學(xué)方程耦合進(jìn)CFD解算器的方法,飛行器運(yùn)動(dòng)時(shí)用到的剛性動(dòng)網(wǎng)格方法,采用松耦合實(shí)現(xiàn)飛行力學(xué)方程與流場(chǎng)控制方程的耦合計(jì)算,建立非定常氣動(dòng)--運(yùn)動(dòng)一體化方法,以及計(jì)算方法的驗(yàn)證等。對(duì)一種采用軌控噴流直接力控制的錐柱裙外形進(jìn)行了數(shù)值模擬,給出了錐柱裙外形在噴流作用下做軌控運(yùn)動(dòng)的非定常計(jì)算結(jié)果,研究了干擾流場(chǎng)結(jié)構(gòu)、力與力矩和飛行器狀態(tài)隨時(shí)間變化特性,并與定常結(jié)果進(jìn)行對(duì)比,分析了外流參數(shù)變化對(duì)氣動(dòng)運(yùn)動(dòng)特性的影響。結(jié)果表明軌控運(yùn)動(dòng)非定常狀態(tài)流場(chǎng)結(jié)構(gòu)和氣動(dòng)力與定常狀態(tài)存在差異,工程實(shí)踐中須要考慮運(yùn)動(dòng)非定常效應(yīng)的影響。
三維可壓縮流動(dòng)N-S方程的量綱為1化形式為
(1)
采用有限體積法,無(wú)黏通量使用二階精度Roe格式[14]進(jìn)行離散;黏性通量使用中心差分格式;時(shí)間推進(jìn)使用雙時(shí)間步方法[15];湍流模型使用一方程S-A湍流模型[16]。為加速收斂、提高計(jì)算效率,使用了基于信息傳遞接口(MPI)的網(wǎng)格分塊并行技術(shù)。邊界條件為壁面無(wú)滑移邊界,遠(yuǎn)場(chǎng)考慮網(wǎng)格動(dòng)態(tài)運(yùn)動(dòng)邊界。
飛行力學(xué)平動(dòng)方程定義在慣性系,它是固定系,如下
(2)
式中,V為線速度,m為質(zhì)量,F(xiàn)為力。
飛行力學(xué)的轉(zhuǎn)動(dòng)方程定義在體軸系,它是動(dòng)坐標(biāo)系,與慣性系相互轉(zhuǎn)換定義姿態(tài)角(γ,θ,ψ),如下
(3)
式中,ω為角速度,I為轉(zhuǎn)動(dòng)慣量,M為力矩。
對(duì)于飛行力學(xué)方程,使用顯式四階龍格-庫(kù)塔方法求解。
耦合方法可分為全耦合、松耦合、緊耦合,考慮到松耦合思路簡(jiǎn)單,擁有更高的計(jì)算效率,本文采用松耦合算法。
松耦合可以使模塊之間求解獨(dú)立,示意圖見圖1。在n時(shí)刻流場(chǎng)信息都是已知量,CFD將力和力矩傳給RBD,求解飛行器位移和姿態(tài)傳給CFD,生成新的網(wǎng)格,求解得到n+1時(shí)刻的流場(chǎng)信息,物理時(shí)間步向前推進(jìn)了一步。
圖1 CFD/RBD松耦合Fig.1 CFD and RBD loose coupling
剛性動(dòng)網(wǎng)格是最簡(jiǎn)單有效的網(wǎng)格變形技術(shù),適用于任意幅度的運(yùn)動(dòng),網(wǎng)格按照飛行器的運(yùn)動(dòng)而一起運(yùn)動(dòng),根據(jù)平動(dòng)位移和轉(zhuǎn)動(dòng)角度直接給出變形之后的網(wǎng)格坐標(biāo),網(wǎng)格速度被指定給每個(gè)網(wǎng)格點(diǎn)。
本節(jié)驗(yàn)證了亞聲速和超聲速條件下兩個(gè)算例。一是76°后掠平板三角翼模型在Ma=0.3下,繞2/3轉(zhuǎn)軸做從0°攻角上仰至 90°攻角的俯仰線性運(yùn)動(dòng),攻角變化公式為
α=kt
(4)
式中,k為減縮頻率,取0.024。圖2給出了升力系數(shù)隨攻角變化曲線,與文獻(xiàn)[17]結(jié)果相比,隨著攻角增大,升力系數(shù)先增大后減小,峰值與風(fēng)洞實(shí)驗(yàn)結(jié)果吻合。
圖2 升力系數(shù)變化曲線Fig.2 Variation curves of lift coefficient
二是標(biāo)模Finner在Ma=2.5下,模擬的是無(wú)初始角速度情況下,噴流在t=0時(shí)開啟,作用時(shí)間20 ms,然后關(guān)閉的周期振蕩運(yùn)動(dòng)。圖3給出了攻角隨時(shí)間變化曲線,幅值隨著時(shí)間增加而減小,與文獻(xiàn)[18]結(jié)果吻合較好。
圖3 攻角變化曲線Fig.3 Variation curves of angle of attack
計(jì)算模型是錐柱裙外形模型,如圖4所示。模型總長(zhǎng)2 250 mm,坐標(biāo)原點(diǎn)O位于飛行器頭部頂點(diǎn),x表示橫向,y表示縱向,軌控噴流中心位于x=1 177 mm 質(zhì)心處。整體計(jì)算網(wǎng)格為C-H型,采用分區(qū)對(duì)接方式,在物理量變化劇烈處加密,網(wǎng)格總數(shù)為146萬(wàn)。圖5給出了網(wǎng)格整體布局。
圖4 錐柱裙外形模型Fig.4 Cone-cylinder-flare model
圖5 網(wǎng)格布局Fig.5 Grids of model
計(jì)算狀態(tài)為:馬赫數(shù)為6,高度為20 km,攻角為0°,噴口中心位于飛行器下表面的質(zhì)心處。來(lái)流條件、噴流條件和飛行器參數(shù)分別見表1、表2和表3。
表1 來(lái)流條件Tab.1 Air conditions
表2 噴流條件Tab.2 Jet conditions
表3 飛行器參數(shù)Tab.3 Aircraft parameters
圖6給出了H=20 km噴流使飛行器軌控運(yùn)動(dòng)的壓力云圖,圖中反映出了干擾流場(chǎng)細(xì)節(jié)隨時(shí)間變化的特征。t=0 ms時(shí),流場(chǎng)呈現(xiàn)上下對(duì)稱狀態(tài),頭部有一道弓形激波,在飛行器尾部形成渦流;t=10 ms時(shí)當(dāng)噴流垂直射入超聲速來(lái)流,就像是形成一道障礙,在噴流前方產(chǎn)生一道很強(qiáng)的弓形激波,噴流反作用推動(dòng)飛行器向上爬升;t=50 ms 時(shí),在噴流產(chǎn)生的干擾俯仰力矩的作用下有約1°抬頭,噴流干擾區(qū)由下表面向上表面覆蓋,噴口前高壓影響范圍擴(kuò)展到了上表面,產(chǎn)生了包裹效應(yīng),造成了總體噴流干擾力的降低。
圖7給出了不同攻角下下表面180°子午線上壓力分布比較圖,壓力分布的不同主要集中在3個(gè)區(qū)域:飛行器前方錐段(區(qū)域1)、噴流出口前方高壓區(qū)段(區(qū)域2)和飛行器后方尾裙段(區(qū)域3)。隨著飛行器運(yùn)動(dòng)向上抬頭,攻角增大,錐段壓力增大,噴流前方壓力增大,尾裙段壓力減小。
(a)t=0 ms
圖7 不同攻角下下表面180°子午線上壓力分布比較圖Fig.7 Pressure distribution of 180° centerline at different angle of attack
(a)俯仰角速度
圖8分別給出了飛行器俯仰角速度和俯仰角隨時(shí)間變化曲線,噴流的持續(xù)作用使飛行器各個(gè)運(yùn)動(dòng)狀態(tài)隨時(shí)間增大。圖8(a)顯示俯仰角速度跟時(shí)間近似成線性關(guān)系,圖8(b)顯示在噴流開啟50 ms后飛行器產(chǎn)生了1°左右的抬頭。
圖9給出法向力系數(shù)和俯仰力矩系數(shù)隨時(shí)間變化曲線,包括總氣動(dòng)力和力矩系數(shù)、噴流本身產(chǎn)生的氣動(dòng)力和力矩系數(shù)和固壁表面除噴流其他區(qū)域的氣動(dòng)力和力矩系數(shù)。對(duì)于法向力系數(shù),在20 km下,隨著飛行器抬頭,總法向力隨時(shí)間推進(jìn)增大,噴流本身推力基本保持不變;對(duì)于俯仰力矩系數(shù),噴流開啟時(shí)俯仰力矩變化劇烈,而后隨著飛行器抬頭,總俯仰力矩隨著俯仰角的增大而增大,并且增大的越來(lái)越快,噴流本身產(chǎn)生的力矩基本保持不變。
(a)法向力系數(shù)
為了考量氣動(dòng)--運(yùn)動(dòng)耦合效應(yīng)帶來(lái)的影響,從兩大方面進(jìn)行運(yùn)動(dòng)非定常狀態(tài)和定常狀態(tài)的比較,一是流場(chǎng)結(jié)構(gòu),包括表面極限流線圖、壓力云圖、分離距離;二是壓力分布和力與力矩對(duì)比,包括下表面180°子午線上壓力分布比較圖、法向力放大系數(shù)和干擾俯仰力矩。
首先給出運(yùn)動(dòng)非定常計(jì)算和定常計(jì)算表面極限流線圖,圖10和圖11是俯視視角,可以看到運(yùn)動(dòng)非定常狀態(tài)與定常狀態(tài)相比,表面流線結(jié)構(gòu)非常相似。圖12和圖13是側(cè)視視角,關(guān)注噴口前表面第一道分離線,運(yùn)動(dòng)非定常狀態(tài)的分離線在飛行器表面所包裹的面積比定常狀態(tài)更大。圖14和圖15是噴口處左視視角,運(yùn)動(dòng)非定常狀態(tài)的下表面高壓區(qū)比上表面覆蓋的范圍更大。定義分離距離為噴管中心到前緣分離線的距離,如圖16所示,可以看到運(yùn)動(dòng)非定常狀態(tài)的分離距離在不同俯仰角下是大于定常狀態(tài)的。
圖10 運(yùn)動(dòng)非定常θ=1°表面極限流線俯視圖Fig.10 Top view of surface limit streamlines of unsteady state at θ=1°
圖11 定常θ=1°表面極限流線俯視圖Fig.11 Top view of surface limit streamlines of steady state at θ=1°
圖12 運(yùn)動(dòng)非定常θ=1°表面極限流線側(cè)視圖Fig.12 Side view of surface limit streamlines of unsteady state at θ=1°
圖13 定常θ=1°表面極限流線側(cè)視圖Fig.13 Side view of surface limit streamlines of steady state at θ=1°
圖14 噴口中心處運(yùn)動(dòng)非定常θ=1°壓力云圖Fig.14 Pressure contours at jet center of unsteady state at θ=1°
圖15 噴口中心處定常θ=1°壓力云圖Fig.15 Pressure contours at jet center of steady state at θ=1°
圖16 運(yùn)動(dòng)非定常與定常狀態(tài)分離距離對(duì)比Fig.16 Separation distance of steady and unsteady state
進(jìn)一步考察更細(xì)致的對(duì)稱線上的壓力分布,圖17給出俯仰角為1°時(shí)運(yùn)動(dòng)非定常狀態(tài)下表面180°子午線上壓力分布與定常狀態(tài)比較圖,可以看出運(yùn)動(dòng)使噴口前分離區(qū)極值點(diǎn)前移、壓力減小,使噴口后尾跡區(qū)壓力增大,使尾裙上壓力減小。
圖17 運(yùn)動(dòng)非定常與定常狀態(tài)下θ=1°表面180°子午線上壓力分布對(duì)比Fig.17 Pressure distribution of 180° centerline between steady and unsteady state at θ=1°
定義法向力放大系數(shù)、干擾俯仰力矩,來(lái)分析噴流干擾氣動(dòng)特性
(5)
圖18給出了定常、不運(yùn)動(dòng)非定常和運(yùn)動(dòng)非定常狀態(tài)計(jì)算結(jié)果的隨俯仰角變化曲線對(duì)比圖。以俯仰角為自變量,定常狀態(tài)定點(diǎn)計(jì)算0°~1°,取間隔0.25°下定常狀態(tài)的法向力放大系數(shù)、干擾俯仰力矩,法向力放大系數(shù)隨俯仰角的增大而增大,且同一俯仰角下定常狀態(tài)比運(yùn)動(dòng)非定常狀態(tài)更大。干擾俯仰力矩也隨俯仰角的增大而增大,在0.25°時(shí)運(yùn)動(dòng)非定常和定常狀態(tài)差別小,而1°時(shí)運(yùn)動(dòng)非定常和定常狀態(tài)差別大。這是由于運(yùn)動(dòng)非定常狀態(tài)飛行器的運(yùn)動(dòng)速度逐漸增大,非定常遲滯效應(yīng)增大,流場(chǎng)變化跟不上飛行器運(yùn)動(dòng),形成阻尼作用。為了進(jìn)一步考察運(yùn)動(dòng)特性帶來(lái)的影響,增加了一組不運(yùn)動(dòng)非定常狀態(tài)的計(jì)算結(jié)果。由圖18中可知,不運(yùn)動(dòng)非定常狀態(tài)的計(jì)算結(jié)果收斂到定常狀態(tài),說(shuō)明運(yùn)動(dòng)帶來(lái)的額外影響是存在的。
(a)法向力放大系數(shù)對(duì)比
分析運(yùn)動(dòng)非定常與定常狀態(tài)差異產(chǎn)生的原因,有飛行器運(yùn)動(dòng)帶來(lái)的遲滯效應(yīng)與噴流干擾兩種因素影響。將飛行器分為三段,頭部錐段是遲滯效應(yīng)主要影響區(qū),噴流柱段和尾裙段是噴流主要的影響區(qū)域,如圖19所示。取圖18中θ=1°的一組狀態(tài),表4和表5分別給出了此時(shí)運(yùn)動(dòng)非定常與定常狀態(tài)的法向力和俯仰力矩的差值在飛行器三段的占比。頭部錐段差異占比分別為56.2%和80.0%,是主要部分,反映氣動(dòng)流場(chǎng)的遲滯效應(yīng)是引起運(yùn)動(dòng)非定常與定常狀態(tài)結(jié)果差異的主要原因,噴流柱段和尾裙段差異占比更少。
表4 運(yùn)動(dòng)非定常與定常狀態(tài)法向力在θ=1°差異占比Tab.4 The difference ratio of the normal force between the unsteady and stady state at θ=1°
表5 運(yùn)動(dòng)非定常與定常狀態(tài)俯仰力矩在θ=1°差異占比Tab.5 The difference ratio of the pitching moment between the unsteady and stady state at θ=1°
圖19 飛行器三段示意圖Fig.19 Three parts of aircraft
為了將遲滯效應(yīng)與噴流干擾隔開對(duì)比,進(jìn)行了一組無(wú)噴強(qiáng)迫俯仰振蕩計(jì)算,運(yùn)動(dòng)方程為
θ(t)=θ0+θmsin(kt)
(6)
其中,θ0為初始俯仰角,設(shè)為0;θm為振幅,取值2°;k為角頻率,取值10π。圖20給出了非定常力與力矩系數(shù)遲滯曲線與定常狀態(tài)對(duì)比,法向力系數(shù)和俯仰力矩系數(shù)始終圍繞靜態(tài)數(shù)據(jù),存在遲滯。向上俯仰時(shí),它不如定常狀態(tài);向下俯仰時(shí)則相反。同樣,表6和表7分別給出了飛行器在θ=1°時(shí)運(yùn)動(dòng)非定常與定常狀態(tài)的法向力和俯仰力矩的差值在飛行器三段的占比,可以看到,無(wú)噴流干擾時(shí),頭部錐段差異占比分別為67.8%和88.4%,比有噴時(shí)占比更大,此時(shí)遲滯效應(yīng)更顯著。
(a)法向力系數(shù)對(duì)比
表6 運(yùn)動(dòng)非定常與定常狀態(tài)法向力在θ=1°差異占比Tab.6 The difference ratio of the normal force between the unsteady and the stationary state at θ=1°
表7 運(yùn)動(dòng)非定常與定常狀態(tài)俯仰力矩在θ=1°差異占比Tab.7 The difference ratio of the pitching moment between the unsteady and the stationary state at θ=1°
對(duì)于遲滯效應(yīng),可以從數(shù)學(xué)分析上進(jìn)行解釋,如圖21所示。來(lái)流速度為U∞,攻角為α,當(dāng)飛行器做俯仰運(yùn)動(dòng)以w的角速度向下旋轉(zhuǎn),前緣處會(huì)有一個(gè)V的向下速度,此時(shí)前緣的實(shí)際速度只有Vreal=U∞sinα-V,實(shí)際攻角被減小,達(dá)到真實(shí)攻角α的時(shí)間被延遲了。
圖21 飛行器俯仰運(yùn)動(dòng)Fig.21 Aircraft pitching motion
本節(jié)關(guān)注不同來(lái)流馬赫數(shù)、飛行高度和初始攻角對(duì)噴流干擾非定常效應(yīng)的影響。以表1的來(lái)流條件作為對(duì)照組,分別改變馬赫數(shù)Ma=6,7,8,高度H=20,25,30 km,初始攻角α0=-2°,0°,2°,對(duì)比分析了不同來(lái)流條件對(duì)運(yùn)動(dòng)狀態(tài)的響應(yīng)。由圖22~圖24可以得出:俯仰角因馬赫數(shù)不同而產(chǎn)生不同的變化趨勢(shì),在Ma=6時(shí)飛行器抬頭,產(chǎn)生正俯仰角;而Ma=8時(shí)飛行器卻低頭,產(chǎn)生負(fù)俯仰角。因高度不同,動(dòng)壓不同,俯仰角增大的幅度不同,高度越高,俯仰角變化越小。俯仰角變化跟初始攻角成非線性關(guān)系,計(jì)算到50 ms 時(shí),α0=2°時(shí)俯仰角變化更大。
圖22 不同馬赫數(shù)下俯仰角比較圖Fig.22 Variation curves of pitching angle at different Mach numbers
圖23 不同高度下俯仰角比較圖Fig.23 Variation curves of pitching angle at different altitudes
圖24 不同初始攻角下俯仰角比較圖Fig.24 Variation curves of pitching angle at different initial attack angles
針對(duì)現(xiàn)階段對(duì)運(yùn)動(dòng)條件下的噴流干擾非定常效應(yīng)研究很少的問(wèn)題開展研究,建立非定常氣動(dòng)--運(yùn)動(dòng)一體化方法,驗(yàn)證了方法的正確性,數(shù)值模擬了錐柱裙外形軌控運(yùn)動(dòng)非定常效應(yīng),研究了干擾流場(chǎng)結(jié)構(gòu)、力與力矩和飛行器狀態(tài)隨時(shí)間變化特性,對(duì)比了運(yùn)動(dòng)非定常計(jì)算和定常計(jì)算的結(jié)果,分析了外流參數(shù)變化對(duì)氣動(dòng)運(yùn)動(dòng)特性的影響。得到如下結(jié)論:
1)飛行器運(yùn)動(dòng)條件下軌控噴流氣動(dòng)干擾存在強(qiáng)烈的非定常和非線性特性。
2)運(yùn)動(dòng)使噴口前分離區(qū)極值點(diǎn)壓力減小,分離距離有所增大。
3)法向力放大系數(shù)和干擾俯仰力矩定常計(jì)算的結(jié)果高于運(yùn)動(dòng)非定常計(jì)算結(jié)果,原因是存在非定常遲滯效應(yīng),流場(chǎng)變化跟不上飛行器運(yùn)動(dòng)。
4)改變來(lái)流馬赫數(shù)會(huì)影響俯仰角變化趨勢(shì),高度越高,俯仰角變化越小,俯仰角變化跟初始攻角成非線性關(guān)系。
本文給出了運(yùn)動(dòng)條件對(duì)軌控噴流控制氣動(dòng)干擾的影響,并得到了一些有意義的結(jié)論。為了滿足工程實(shí)踐中精度要求,使數(shù)值模擬結(jié)果更逼近真實(shí)飛行,須要考慮這種影響,以便反映更真實(shí)的物理現(xiàn)象。