軌控噴流干擾氣動
--運動一體化非定常數值模擬

謝信輝，柳煜瑋，劉耀峰

(中國航天空氣動力技術研究院，北京 100074)

0 引言

大推力噴流控制技術可實現大機動、強突防[1]，在多種飛行器上都有廣泛的應用，如低速狀態下橫流中煙囪羽流[2]和高速狀態下飛行器攔截系統[3]。特別是在攔截任務的最后階段，要求攔截飛行器具有高度的機動性，提供直接作用力、快速響應的側向噴流控制推進器，這是對傳統氣動控制面的有效補充。

噴流產生的三維流場是高度復雜的，包含復雜的波系、渦系結構及其相互干擾，呈現高度非定常、非線性、強湍流特性[4]，現階段對運動條件下非定常效應的研究很少。早在1964年Zukoski等[5]就開展了不同馬赫數下聲速噴管欠膨脹噴流噴入超聲速外流的實驗研究，指出了弓形激波、分離激波、分離區等主要特征。隨著計算流體力學的發展，CFD技術在噴流干擾流場研究領域得到廣泛應用。1991年張涵信等[6]利用NND格式，完成了鈍錐完全氣體側噴干擾流場的計算。2005年徐敏等[7]研究了脈沖發動機噴流的非定常性對飛行器控制的關系。

當噴流開啟與關閉時，存在強烈的非定常特性，同時伴隨飛行器運動與姿態變化，存在氣動特性與運動特性的強耦合，為了有效評估飛行性能，須要考慮氣動運動耦合帶來的非定常問題。關于空氣動力學與飛行力學耦合方法的研究，國外近年來發展和應用了大規模模擬的最新數值方法[8-11]。Sahu等[12]在超級計算機上對飛行器的“虛擬飛行”進行了研究，使用耦合的CFD/RBD技術對飛行器的實際飛行軌跡和所有相關的非定常自由飛行空氣動力學進行數值預測。Silton等[10]通過與自由飛行實驗數據進行比較，完成了對虛擬飛行模擬計算結果的驗證。劉耀峰等[1]應用非定常方法數值模擬了側向噴流干擾流場建立與消退過程。陳堅強等[13]模擬了舵面運動與側向噴流間耦合干擾問題。以上研究具有豐富的價值，本文的重點是將氣動--運動一體化方法應用到噴流干擾問題中，對運動條件對噴流干擾的影響做進一步的探索研究，為工程實踐的應用提供參考。

本文基于雷諾平均N-S方程和剛體六自由度運動方程，介紹了將飛行力學方程耦合進CFD解算器的方法，飛行器運動時用到的剛性動網格方法，采用松耦合實現飛行力學方程與流場控制方程的耦合計算，建立非定常氣動--運動一體化方法，以及計算方法的驗證等。對一種采用軌控噴流直接力控制的錐柱裙外形進行了數值模擬，給出了錐柱裙外形在噴流作用下做軌控運動的非定常計算結果，研究了干擾流場結構、力與力矩和飛行器狀態隨時間變化特性，并與定常結果進行對比，分析了外流參數變化對氣動運動特性的影響。結果表明軌控運動非定常狀態流場結構和氣動力與定常狀態存在差異，工程實踐中須要考慮運動非定常效應的影響。

1 數值模擬方法

1.1 流體力學方程

三維可壓縮流動N-S方程的量綱為1化形式為

(1)

采用有限體積法，無黏通量使用二階精度Roe格式[14]進行離散；黏性通量使用中心差分格式；時間推進使用雙時間步方法[15]；湍流模型使用一方程S-A湍流模型[16]。為加速收斂、提高計算效率，使用了基于信息傳遞接口(MPI)的網格分塊并行技術。邊界條件為壁面無滑移邊界，遠場考慮網格動態運動邊界。

1.2 飛行力學方程

飛行力學平動方程定義在慣性系，它是固定系，如下

(2)

式中，V為線速度，m為質量，F為力。

飛行力學的轉動方程定義在體軸系，它是動坐標系，與慣性系相互轉換定義姿態角(γ，θ，ψ)，如下

(3)

式中，ω為角速度，I為轉動慣量，M為力矩。

對于飛行力學方程，使用顯式四階龍格-庫塔方法求解。

1.3 氣動--運動耦合方法

耦合方法可分為全耦合、松耦合、緊耦合，考慮到松耦合思路簡單，擁有更高的計算效率，本文采用松耦合算法。

松耦合可以使模塊之間求解獨立，示意圖見圖1。在n時刻流場信息都是已知量，CFD將力和力矩傳給RBD，求解飛行器位移和姿態傳給CFD，生成新的網格，求解得到n+1時刻的流場信息，物理時間步向前推進了一步。

圖1 CFD/RBD松耦合Fig.1 CFD and RBD loose coupling

1.4 剛性動網格方法

剛性動網格是最簡單有效的網格變形技術，適用于任意幅度的運動，網格按照飛行器的運動而一起運動，根據平動位移和轉動角度直接給出變形之后的網格坐標，網格速度被指定給每個網格點。

1.5 數值方法驗證

本節驗證了亞聲速和超聲速條件下兩個算例。一是76°后掠平板三角翼模型在Ma=0.3下，繞2/3轉軸做從0°攻角上仰至 90°攻角的俯仰線性運動，攻角變化公式為

α=kt

(4)

式中，k為減縮頻率，取0.024。圖2給出了升力系數隨攻角變化曲線，與文獻[17]結果相比，隨著攻角增大，升力系數先增大后減小，峰值與風洞實驗結果吻合。

圖2 升力系數變化曲線Fig.2 Variation curves of lift coefficient

二是標模Finner在Ma=2.5下，模擬的是無初始角速度情況下，噴流在t=0時開啟，作用時間20 ms，然后關閉的周期振蕩運動。圖3給出了攻角隨時間變化曲線，幅值隨著時間增加而減小，與文獻[18]結果吻合較好。

圖3 攻角變化曲線Fig.3 Variation curves of angle of attack

2 計算模型和計算條件

計算模型是錐柱裙外形模型，如圖4所示。模型總長2 250 mm，坐標原點O位于飛行器頭部頂點，x表示橫向，y表示縱向，軌控噴流中心位于x=1 177 mm 質心處。整體計算網格為C-H型，采用分區對接方式，在物理量變化劇烈處加密，網格總數為146萬。圖5給出了網格整體布局。

圖4 錐柱裙外形模型Fig.4 Cone-cylinder-flare model

圖5 網格布局Fig.5 Grids of model

計算狀態為：馬赫數為6，高度為20 km，攻角為0°，噴口中心位于飛行器下表面的質心處。來流條件、噴流條件和飛行器參數分別見表1、表2和表3。

表1 來流條件Tab.1 Air conditions

表2 噴流條件Tab.2 Jet conditions

表3 飛行器參數Tab.3 Aircraft parameters

3 結果與分析

3.1 軌控運動干擾流場發展歷程

圖6給出了H=20 km噴流使飛行器軌控運動的壓力云圖，圖中反映出了干擾流場細節隨時間變化的特征。t=0 ms時，流場呈現上下對稱狀態，頭部有一道弓形激波，在飛行器尾部形成渦流；t=10 ms時當噴流垂直射入超聲速來流，就像是形成一道障礙，在噴流前方產生一道很強的弓形激波，噴流反作用推動飛行器向上爬升；t=50 ms 時，在噴流產生的干擾俯仰力矩的作用下有約1°抬頭，噴流干擾區由下表面向上表面覆蓋，噴口前高壓影響范圍擴展到了上表面，產生了包裹效應，造成了總體噴流干擾力的降低。

圖7給出了不同攻角下下表面180°子午線上壓力分布比較圖，壓力分布的不同主要集中在3個區域：飛行器前方錐段(區域1)、噴流出口前方高壓區段(區域2)和飛行器后方尾裙段(區域3)。隨著飛行器運動向上抬頭，攻角增大，錐段壓力增大，噴流前方壓力增大，尾裙段壓力減小。

(a)t=0 ms

圖7 不同攻角下下表面180°子午線上壓力分布比較圖Fig.7 Pressure distribution of 180° centerline at different angle of attack

(a)俯仰角速度

3.2 運動狀態和氣動系數隨時間變化特性

圖8分別給出了飛行器俯仰角速度和俯仰角隨時間變化曲線，噴流的持續作用使飛行器各個運動狀態隨時間增大。圖8(a)顯示俯仰角速度跟時間近似成線性關系，圖8(b)顯示在噴流開啟50 ms后飛行器產生了1°左右的抬頭。

圖9給出法向力系數和俯仰力矩系數隨時間變化曲線，包括總氣動力和力矩系數、噴流本身產生的氣動力和力矩系數和固壁表面除噴流其他區域的氣動力和力矩系數。對于法向力系數，在20 km下，隨著飛行器抬頭，總法向力隨時間推進增大，噴流本身推力基本保持不變；對于俯仰力矩系數，噴流開啟時俯仰力矩變化劇烈，而后隨著飛行器抬頭，總俯仰力矩隨著俯仰角的增大而增大，并且增大的越來越快，噴流本身產生的力矩基本保持不變。

(a)法向力系數

3.3 運動非定常與定常結果對比

為了考量氣動--運動耦合效應帶來的影響，從兩大方面進行運動非定常狀態和定常狀態的比較，一是流場結構，包括表面極限流線圖、壓力云圖、分離距離；二是壓力分布和力與力矩對比，包括下表面180°子午線上壓力分布比較圖、法向力放大系數和干擾俯仰力矩。

首先給出運動非定常計算和定常計算表面極限流線圖，圖10和圖11是俯視視角，可以看到運動非定常狀態與定常狀態相比，表面流線結構非常相似。圖12和圖13是側視視角，關注噴口前表面第一道分離線，運動非定常狀態的分離線在飛行器表面所包裹的面積比定常狀態更大。圖14和圖15是噴口處左視視角，運動非定常狀態的下表面高壓區比上表面覆蓋的范圍更大。定義分離距離為噴管中心到前緣分離線的距離，如圖16所示，可以看到運動非定常狀態的分離距離在不同俯仰角下是大于定常狀態的。

圖10 運動非定常θ=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 運動非定常θ=1°表面極限流線側視圖Fig.12 Side view of surface limit streamlines of unsteady state at θ=1°

圖13 定常θ=1°表面極限流線側視圖Fig.13 Side view of surface limit streamlines of steady state at θ=1°

圖14 噴口中心處運動非定常θ=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 運動非定常與定常狀態分離距離對比Fig.16 Separation distance of steady and unsteady state

進一步考察更細致的對稱線上的壓力分布，圖17給出俯仰角為1°時運動非定常狀態下表面180°子午線上壓力分布與定常狀態比較圖，可以看出運動使噴口前分離區極值點前移、壓力減小，使噴口后尾跡區壓力增大，使尾裙上壓力減小。

圖17 運動非定常與定常狀態下θ=1°表面180°子午線上壓力分布對比Fig.17 Pressure distribution of 180° centerline between steady and unsteady state at θ=1°

定義法向力放大系數、干擾俯仰力矩，來分析噴流干擾氣動特性

(5)

圖18給出了定常、不運動非定常和運動非定常狀態計算結果的隨俯仰角變化曲線對比圖。以俯仰角為自變量，定常狀態定點計算0°～1°，取間隔0.25°下定常狀態的法向力放大系數、干擾俯仰力矩，法向力放大系數隨俯仰角的增大而增大，且同一俯仰角下定常狀態比運動非定常狀態更大。干擾俯仰力矩也隨俯仰角的增大而增大，在0.25°時運動非定常和定常狀態差別小，而1°時運動非定常和定常狀態差別大。這是由于運動非定常狀態飛行器的運動速度逐漸增大，非定常遲滯效應增大，流場變化跟不上飛行器運動，形成阻尼作用。為了進一步考察運動特性帶來的影響，增加了一組不運動非定常狀態的計算結果。由圖18中可知，不運動非定常狀態的計算結果收斂到定常狀態，說明運動帶來的額外影響是存在的。

(a)法向力放大系數對比

分析運動非定常與定常狀態差異產生的原因，有飛行器運動帶來的遲滯效應與噴流干擾兩種因素影響。將飛行器分為三段，頭部錐段是遲滯效應主要影響區，噴流柱段和尾裙段是噴流主要的影響區域，如圖19所示。取圖18中θ=1°的一組狀態，表4和表5分別給出了此時運動非定常與定常狀態的法向力和俯仰力矩的差值在飛行器三段的占比。頭部錐段差異占比分別為56.2%和80.0%，是主要部分，反映氣動流場的遲滯效應是引起運動非定常與定常狀態結果差異的主要原因，噴流柱段和尾裙段差異占比更少。

表4 運動非定常與定常狀態法向力在θ=1°差異占比Tab.4 The difference ratio of the normal force between the unsteady and stady state at θ=1°

表5 運動非定常與定常狀態俯仰力矩在θ=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

為了將遲滯效應與噴流干擾隔開對比，進行了一組無噴強迫俯仰振蕩計算，運動方程為

θ(t)=θ0+θmsin(kt)

(6)

其中，θ0為初始俯仰角，設為0；θm為振幅，取值2°；k為角頻率，取值10π。圖20給出了非定常力與力矩系數遲滯曲線與定常狀態對比，法向力系數和俯仰力矩系數始終圍繞靜態數據，存在遲滯。向上俯仰時，它不如定常狀態；向下俯仰時則相反。同樣，表6和表7分別給出了飛行器在θ=1°時運動非定常與定常狀態的法向力和俯仰力矩的差值在飛行器三段的占比，可以看到，無噴流干擾時，頭部錐段差異占比分別為67.8%和88.4%，比有噴時占比更大，此時遲滯效應更顯著。

(a)法向力系數對比

表6 運動非定常與定常狀態法向力在θ=1°差異占比Tab.6 The difference ratio of the normal force between the unsteady and the stationary state at θ=1°

表7 運動非定常與定常狀態俯仰力矩在θ=1°差異占比Tab.7 The difference ratio of the pitching moment between the unsteady and the stationary state at θ=1°

對于遲滯效應，可以從數學分析上進行解釋，如圖21所示。來流速度為U∞，攻角為α，當飛行器做俯仰運動以w的角速度向下旋轉，前緣處會有一個V的向下速度，此時前緣的實際速度只有Vreal=U∞sinα-V，實際攻角被減小，達到真實攻角α的時間被延遲了。

圖21 飛行器俯仰運動Fig.21 Aircraft pitching motion

3.4 來流參數的影響

本節關注不同來流馬赫數、飛行高度和初始攻角對噴流干擾非定常效應的影響。以表1的來流條件作為對照組，分別改變馬赫數Ma=6，7，8，高度H=20，25，30 km，初始攻角α0=-2°，0°，2°，對比分析了不同來流條件對運動狀態的響應。由圖22～圖24可以得出：俯仰角因馬赫數不同而產生不同的變化趨勢，在Ma=6時飛行器抬頭，產生正俯仰角；而Ma=8時飛行器卻低頭，產生負俯仰角。因高度不同，動壓不同，俯仰角增大的幅度不同，高度越高，俯仰角變化越小。俯仰角變化跟初始攻角成非線性關系，計算到50 ms 時，α0=2°時俯仰角變化更大。

圖22 不同馬赫數下俯仰角比較圖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

4 結論

針對現階段對運動條件下的噴流干擾非定常效應研究很少的問題開展研究，建立非定常氣動--運動一體化方法，驗證了方法的正確性，數值模擬了錐柱裙外形軌控運動非定常效應，研究了干擾流場結構、力與力矩和飛行器狀態隨時間變化特性，對比了運動非定常計算和定常計算的結果，分析了外流參數變化對氣動運動特性的影響。得到如下結論：

1)飛行器運動條件下軌控噴流氣動干擾存在強烈的非定常和非線性特性。

2)運動使噴口前分離區極值點壓力減小，分離距離有所增大。

3)法向力放大系數和干擾俯仰力矩定常計算的結果高于運動非定常計算結果，原因是存在非定常遲滯效應，流場變化跟不上飛行器運動。

4)改變來流馬赫數會影響俯仰角變化趨勢，高度越高，俯仰角變化越小，俯仰角變化跟初始攻角成非線性關系。

本文給出了運動條件對軌控噴流控制氣動干擾的影響，并得到了一些有意義的結論。為了滿足工程實踐中精度要求，使數值模擬結果更逼近真實飛行，須要考慮這種影響，以便反映更真實的物理現象。