小展弦比飛翼標模非定常流動及自由搖滾特性

王方劍，解克，劉金，宋玉輝，秦漢，陳蘭

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

飛翼是一種高度融合翼身組合體布局，全機沒有平尾、垂尾等安定面，也沒有傳統意義上的機身，在外形上體現出平滑過渡、高度融合的幾何特征，使戰斗機在各種高度、姿態下的隱身性和機動性都得到很好的兼顧，是實現隱身/氣動一體化設計的理想選擇。國內外在以戰斗機為背景的飛翼布局研究中，主要集中在中小展弦比飛翼布局，其中典型的有4 個研究外形：NASA 的ICE 外 形［1］、國 際 合 作 計 劃 中 的1303 外 形［2］、NATO/RTO AVT-161 項 目 組 采 用 的SACCON 外形［3］以及中國“十二五”期間自主研發的小展弦比飛翼標模（FDM）［4］。

在飛翼布局靜態流動特性方面，Tormalm等［5］分別采用SA、SARAL、SST 湍流模型以及DDES 方法對SACCON 外形靜態流場進行計算并對比分析，分析表明在小攻角區域，SARAL、SST 湍流模型都能夠得到較好的模擬效果，而在較大攻角區域，對背風側分離區的模擬建議采用DES 方 法 模 擬 復 雜 的 流 動 結 構。Frink［6］分 析 了SACCON 隨攻角增加，背風側的流動發展非常復雜，存在頭部渦、厚度渦、翼尖渦，攻角增加之后厚度渦與翼尖渦融合成為一個大渦，體現出飛翼布局復雜的流動特性。蘇繼川等［4］通過數值模擬研究了小展弦比飛翼標模的跨聲速典型氣動特性，主要以渦升力為出發點，分析了渦的產生、發展直至破裂的整個過程。張耀冰等［7］采用數值模擬方法研究了雷諾數對小展弦比飛翼標模流場氣動特性的影響問題。結果表明大攻角時，雷諾數增大，除引起摩阻減小外，壓阻也有明顯增大。

對于無尾布局飛行器非定常流動特性研究主要集中在三角翼模型上，Gursul［8］在試驗中發現三角翼旋渦的破裂點會沿著渦軸進行振蕩，并研究了前緣后掠角對振蕩頻率的影響。Menke等［9］通過采用風洞試驗方法，將三角翼的流動頻率分為4 種模態，分別是渦破裂點振蕩頻率、渦脫落頻率、螺旋波的失穩頻率、剪切層的失穩頻率。當模型攻角逐漸增大，先出現的是旋渦破裂（一般是螺旋式破裂），渦破裂點出現非定常振蕩，并且其頻率特性與后掠角相關。試驗發現［10］旋渦破裂的頻率由準周期的頻率以及高頻率低幅值的振蕩組成，與此同時還伴隨著非對稱渦的位置左右切換，隨著攻角增加，渦破裂點沿渦軸前移，左右旋渦干擾加劇。祝立國和呂志詠［11］也在水洞中發現了此現象，Ericsson［12］在之后的研究中詳細分析了旋渦破裂的流態及其頻率特性。張明祿等［13］對75°后掠三角翼進行脈動壓力測量，發現通過對信號的濾波處理，可知整個三角翼上翼面的壓力脈動以螺旋波的影響為最大，破裂渦流中螺旋波的準周期壓力脈動造成了機翼的抖振。Li［14］和Lin［15］等采用數 值模擬方 法，對雙 三角翼的大攻角渦破裂特性進行分析，發現了破裂形態的轉換、破裂點的振蕩和破裂后螺旋渦結構的脫落，氣動力頻率特征決定于螺旋渦結構的脫落頻率。

飛行器大攻角飛行時，由于背風側分離流動極為復雜，易引發運動失穩，其中較為常見的為以繞體軸滾轉振蕩為主的機翼搖滾運動。Ross等［16］對大量戰斗機非指令運動進行統計研究，將機翼搖滾現象分為兩類，一種為跨聲速較小攻角下的機翼搖滾運動［10，17］，其運動失穩的主要原因為激波誘導的流動失速；另一種主要發生在低、亞聲速，在較大攻角出現的機翼搖滾運動主要由前體渦、機翼渦相互干擾帶來。 Ericsson等［12，18-21］在 機 翼 搖 滾機理方面做了多方面的工作，在理論方面做出了基礎性的貢獻。唐敏中［22］、鄧學鎣［23］、王方劍［24］、楊曉峰［25］等采用風洞自由滾轉試驗、流動顯示的方法，對細長體機翼搖滾的機理做出了很有意義的研究，劉偉［26-27］、楊云軍［28］等采用數值模擬方法，對細長三角翼在搖滾過程中的流動變化進行進一步的深入分析，并采用非線性動力學方法對其運動歷程進行建模，取得很好的成效。

從飛翼布局流動特性研究現狀可以看到，相比于三角翼，飛翼布局復雜的外形設計會對流動特性產生較大影響，并且對于飛翼布局流動特性的研究主要集中在定常流動范疇，對于飛翼布局的非定常流動以及運動失穩特性研究較少。本文采用延遲脫體渦模擬（DDES）方法模擬飛翼標模單自由度自由搖滾運動，并采用動力學模態分解（DMD）、物面脈動壓力分析等方法分析飛翼標模非定常特性、自由搖滾機理，為中國飛翼布局飛行器的研制提供理論支撐。

1 DDES 方法及算例驗證

控制方程為三維可壓縮Navier-Stokes 方程，廣義坐標表示為

式中：F、G、H為對流通量；Fv、Gv、Hv為黏性通量。

空間離散采用基于MUSCL 插值方法的FDS-Roe 格式，黏性項采用二階中心差分離散，對流和壓力項使用三階迎風偏置格式。時間推進格式為隱式LUSGS 方法，保證了較高的時間計算精度。非定常計算采用雙時間步長方法，同時采用多重網格算法加快子迭代的收斂速度。在用湍流模型進行流動模擬時，控制方程中應加入雷諾應力相關項。本文采用的湍流模型是基于Menter 的兩方程剪切應力輸運湍流模型（SST），表達式為

為了具有DES 性能，DES 方法的一般思想是重新定義SST 湍流模型k方程的耗散項，其形式為

為了解決MSD（Modeled Stress Depletion）與GIS（Grid Induced Separation）問題，Spalart等［27］提出了延遲函數fd，將SST 湍流模型的湍流能量輸運方程的耗散項重新定義為與網格尺度和湍流黏度有關的DES 長度尺度。延時函數的表達式為

式 中：Ui，j為速度梯度；υt為湍流黏度系數；κ=0.41 為卡門常數；d為網格到物面的距離。結合延遲函數fd的長度尺度dDES被重新定義為

因此，dDES不僅與網格尺度有關，而且與湍流黏度有關，可以防止湍流模擬過早地切換到LES模式。

采用基于三階迎風偏置格式的DDES 方法來開展研究，主要基于以下考慮：

1） DDES 方法對網格依賴性較好，能夠增強RANS 與LES 分界面的識別和過渡能力，從而減少人工計算參數設置帶來的不確定性。雖然DES 類方法存在固有的“灰區”問題，在局部位置對小尺度湍流脈動解析會引入數值模擬誤差，但小尺度旋渦能量較低，誘導物面產生的附加氣動力較小，對整體流場結構及大尺度旋渦演化規律影響較小，所以可認為該方法是適用的。

2） 三階迎風格式與中心格式相比耗散相對較大，這會影響對小尺度旋渦的解析。但由于本文重點在于對機翼搖滾及其非定常現象的研究，小尺度旋渦特征頻率比飛行器搖滾頻率高兩個量級，小尺度、高頻率的湍流脈動將難以對飛行器的整體運動構成明顯影響。因此，這里選用迎風型格式來開展研究理論上是可行的。

測試算例采用的計算外形為NACA0021 翼型，翼型弦長為125 mm，展長為500 mm，翼型兩端采用周期邊界條件。計算網格采用結構網格進行求解，總數為500 萬，圖1 展示了橫截面網格分布，為了保證對翼型后緣微小分離的模擬精度以及翼型背風側的大尺度分離流動的模擬效果，在進行網格劃分時，分別在翼型后緣與模型的背風側法向進行了加密。由于算例計算風速為低速，在計算時引入預處理，使得計算魯棒性與精度都得到了相應的提高。計算參數為：Ma=0.1，攻角為α=60°，Re=0.27×106，物理時間步為1×10-4s，共計算了10 000 個物理時間步，計算物理時間為1 s。圖2 為NACA0021 橫截面的瞬時渦量，從圖中可以看到，通過采用DDES 計算方法，能夠捕捉到NACA0021 翼型背風側的復雜且豐富的旋渦結構。

圖1 橫截面網格分布Fig. 1 Mesh of cross section

圖2 NACA0021 橫截面瞬時渦量Fig. 2 Transient vorticity of NACA0021 cross section

圖3中黑色曲線為模型中間橫截面的壓力Cp分布時均值，紅色空心圓點為風洞測壓試驗中相應截面的壓力分布。從圖中可以看到，試驗值［29］與計算值在量值上和趨勢上都保持了較好的一致性，證明SST-DDES 方法在氣動力計算上能夠取得較好的效果。圖4 為計算值與試驗值的頻譜對比情況。橫軸為斯特勞哈爾數St，縱軸為功率譜密度PSD。紅色曲線為通過風洞試驗測得的升力時序數據辨識出的頻譜，黑色曲線為通過CFD 方法得到的升力功率譜。從圖中的對比能夠看到，試驗值與計算值在量級上與趨勢上都能夠得到吻合較好。計算值能夠捕捉到背風側流動的主頻（St=0.2），同時能夠捕捉到2 倍頻（St=0.4），說明SST-DDES 方法在時序過程以及頻率捕捉上也能夠有比較好的計算效果。

圖3 試驗與計算壓力分布對比Fig. 3 Pressure distribution comparison between experiment and CFD

2 動力學模態分解方法

動力學模態分解（DMD）方法只基于流動快照而不受模型限制，通過提取流動中的模態，準確地描述流動結構。對于非線性流動而言，DMD 的模態描述了在數據序列中起支配作用的流動結構。其基本思路是基于線性化假設，假設每一個瞬時的物理量都可以用前一個瞬時線性表示，求解變換矩陣的相似矩陣，對相似矩陣進行特征分解得到對應的流動模態，其優勢在于，每階模態對應都是單一頻率，所以對于多頻復合的復雜流動下，分析單一頻率下的模態特性，能夠對研究復雜物理現象有一定的參考意義。以下為其數學推導過程。

DMD 方法假設每一個瞬時物理量可以用前一個瞬時線性表示，所以定義n個流場快照中，前n-1 個流場快照為X，后n-1 個流場快照為Y，所以可以表示為

式（9）中A表示線性變換矩陣，由于每個流場快照數據量非常龐大，這樣矩陣A的維度就會非常大，不便于求解。所以DMD 方法是通過求解A的相似矩陣來實現。對X進行奇異值分解，得到左奇異矩陣U以及右奇異矩陣V。其中Σ為奇異值矩陣，選取前r階奇異值，從而得到降階的效果。同時將矩陣A向U方向投影得到式（11）。

將式（10）、式（11）代入式（9），得到相似矩陣的表達形式如下：

將相似矩陣進行特征分解，得到主要特征值和特征向量，記第j個特征值和特征向量為μj和wj，則定義第j個模態為

其對應的模態頻率為

最終每個瞬時的流場快照都可以表示為模態、振幅、動力學矩陣乘積。

3 小展弦比飛翼標模及計算網格

采用國內小展弦飛翼布局標模［4］，其三維示意圖與坐標如圖5、圖6所示。模型縮比為1∶25，模型總長L=0.612 m，翼展Span=0.457 m，參考面積S=0.131 6 m2，平均氣動弦長（參考長度）c=0.382 m。計算馬赫數：Ma=0.6，攻角α=0°～30°，側滑角β=0°，Re=8.5×106，力矩參考點為（0.352，0，0.016 5）m，機頭頂點為坐標原點。

圖5 小展弦比飛翼標模外形尺寸Fig. 5 Configuration of low aspect ratio flying-wing

圖6 飛翼標模坐標Fig. 6 Coordinate system of flying-wing

圖7為網格空間拓撲結構，全流場采用全O型結構網格，遠場距物面距離為20 個參考長度，y+＜1。為保證對流動分離細節的捕捉，在飛行器頭部、機翼前緣、翼稍進行適當加密，在模型旋渦相應位置的法向與周向也進行了適當加密。

圖7 網格空間拓撲圖Fig. 7 Mesh space topology

由于針對的是飛翼布局大攻角下非定常流動特性問題，其對旋渦尺度的分辨十分重要，所以需要對網格無關性進行檢查，圖8 為飛翼標模物面網格拓撲結構，網格量分別為750 萬、1 500萬、3 000 萬，3 套網格在流向、周向、法向均進行變化。選取的典型工況為Ma=0.6，α=24.5°，并選取升力CL、阻力CD和俯仰力矩Cm的平均值進行分析。表1 為3 套網格下的升力CL、阻力CD、俯仰力矩Cm計算結果。從表1 可以看到，網格規模對CL、CD、Cm均有一定的影響，隨著網格量的增加，當網格量達到1 500 萬時，其計算結果（CL、CD、Cm）與網格量3 000 萬較為接近，證明網格量1 500 萬對于研究該問題已滿足要求。

圖8 物面網格Fig. 8 Surface mesh

4 飛翼標模靜態流動特性

圖9 為Ma=0.6 下攻角α=0°～30°范圍內的流動變化，從圖中可以看到，攻角α=5°時，飛翼標模在翼稍開始卷起對稱渦流動結構，旋渦起始點在距頭部頂點X/L=0.6左右，隨著攻角持續增大至α=15°，對稱渦起始點逐漸移動至模型頂點，當攻角增大至α=20°時，集中渦變得不穩定，開始出現渦破裂現象，并且模型兩側旋渦開始呈現明顯的非對稱形態，隨著攻角繼續增大至α=30°，旋渦破裂點逐漸向頭部移動，旋渦破裂后呈現出尺度較小的螺旋旋渦結構，流動也體現出了明顯的非定常流動特性。

圖9 飛翼標模靜態流動變化Fig. 9 Static flow of flying-wing

5 飛翼標模非定常流動特性及DMD 分析

飛翼標模在大攻角下具有非常復雜的非定常流動特性，本節以亞聲速為例，對其非定常流動特性進行細致分析，計算Ma=0.6，攻角α=26°，計算物理時間步長為5×10-5s。如圖10 所示，從圖中可以看到，飛翼標模出現了類似三角翼典型的分離渦、渦破裂、螺旋波結構。飛翼前緣誘導分離集中渦，從模型頭部起始，沿流向向后發展，其中渦軸速度在X/L=0.4 位置降至很低，隨即發生旋渦破裂，在旋渦破裂處，渦量由于無法繼續沿渦軸向后發展，會由內而外產生螺旋波流動。為了更好的研究空間流動對物面局部壓力分布的影響，在旋渦影響較大的區域布置監測點，如圖10 所示，其中沿渦軸布置了5 個點，分別監測流動對物面的誘導能力，其在物面X方向位置分別為X/L=0.28，0.41，0.5，0.56，0.7。

圖10 物面監測點布置（Ma=0.6）Fig. 10 Location of monitored points （Ma=0.6）

圖11是飛翼標模空間流動隨時間的變化，從圖中可以看到，隨著時間變化，旋渦飛翼旋渦破裂位置隨著時間也會有一定程度的變化，當時間t=0 s時，旋渦破裂位置在X/L=0.25 的位置，隨著時間推移，其旋渦破裂位置逐漸向下游移動，當t=3.0×10-3s時，旋渦破裂發展到約X/L=0.3的位置，然后旋渦破裂位置繼續向頭部方向移動，從以下圖中也可以看出，即使飛翼標模是靜止的，旋渦破裂位置也會在一定范圍內變化，頻率較低，而螺旋波流動體現出的是一種高頻的流動變化。圖12 是沿旋渦渦軸切面的渦量值隨時間變化，依然能夠較為明顯的看到旋渦破裂位置隨時間的變化。在旋渦破裂以前，渦量處于穩定的狀態，旋渦破裂以后，渦軸的渦量分散成若干高頻小尺度渦向后發展。這也體現出了飛翼標模在大攻角下，背風側復雜非定常流動由不同尺度、不同頻率的流動組合而成。

圖11 非定常旋渦流動隨時間變化（α=26°，Ma=0.6，FDM）Fig. 11 Unsteady vortical flow of flying-wing vs time（α=26°，Ma=0.6，FDM）

圖12 渦軸切面渦量分布（α=26°，Ma=0.6，FDM）Fig. 12 Vorticity distribution on section of vortex axis（α=26°，Ma=0.6，FDM）

圖13是旋渦破裂點在x方向的位置隨時間的變化，從圖中可以清晰的看到，旋渦破裂點隨時間呈現出明顯的非定常性，其變化范圍在X/L=0.29～0.45，波動范圍為模型全長的16%左右。

圖13 渦破裂點位置隨時間變化曲線Fig. 13 Location of vortex breakdown point with time

圖14是監測點非定常壓力脈動，中間一列曲線為壓力（P/Pinf： 當地壓力除以來流壓力）隨時間變化，右側為其對應的頻譜。從圖中可以看到，在不同的測壓點由于空間流動不同，誘導出了不同的壓力脈動特性。在X/L=0.28 處，能看到其壓力脈動幅值很小，只是有一些微小的波動，從頻譜中也沒有看到明顯的頻率，此處為穩定的集中渦，所以沒有明顯的脈動特性。在X/L=0.41 處，空間流動處于旋渦剛破裂的區域，由于旋渦破裂點會隨著時間的變化沿渦軸前后移動，所以在壓力脈動上產生了較為明顯的脈動，從頻譜中能夠發現，其中主要頻率St數（St=fc/V∞）在0.12～0.23左右，也揭示了旋渦破裂的主要頻率在這個范圍。在X/L=0.5，0.56 處，旋渦已經發生破裂，并且旋渦破裂點不會下移到此位置，此處流動為渦破裂后的“泡型”膨脹區域，從脈動壓力上能夠看到，此處的脈動壓力值較小，沒有體現出很明顯的主頻，X/L=0.5處，由于受上游渦破裂點移動的影響，在頻譜上輕微的體現出了主頻范圍，與X/L=0.41 處頻率范圍接近，而X/L=0.56 處，渦量已經充分膨脹，沒有體現出明顯主頻。當流動發展至X/L=0.7，渦量膨脹后開始誘導出高頻的螺旋波流動結構，其中誘導的壓力幅值開始變大，從頻譜中可看到，其主頻范圍為St=1.16～2.33 范圍內，體現出了螺旋波的高頻特性。

本文還采用了DMD 方法對Ma=0.6、α=26°工況下大攻角非定常流動進行初步分析，共計算了300 個瞬時空間流動，時間間隔為5×10-4s。由于直接計算全流場的DMD 分解會帶來極大的計算量，所以這里將沿渦軸切面作為表征面來進行分析，將這300個表征面的速度值進行DMD 分析，并提取出流動模態對應的頻率。圖15為不同頻率對應的DMD 模態，可以看到，其中St=0對應模態體現為典型的渦破裂形式，沿渦軸方向速度值在渦破裂前較大，破裂后速度迅速降低，并產生一定回流，該模態類似整體流動的“平均值”，體現了旋渦流動的總體效果。在St=0.25 模態中發現，在渦破裂位置附近體現出了一定的流動非定常變化，同時在下游也體現出了一定的小尺度流動，該模態頻率也在旋渦破裂點振蕩頻率范圍內。在St=0.39 模態中，可以看到，渦破裂點處已經沒有明顯的脈動特征，該頻率下的流動主要體現為下游的小尺度流動，同時有一定的旋渦脫落以及切換特性。在St=1.2 模態中，可以看到，該模態頻率在螺旋波頻率內，從模態中也能夠比較清晰的看到，更多更小尺度的旋渦脫落與切換出現在接近模型尾部的地方，這也是受空間螺旋波流動影響的一種直接體現。

圖15 渦軸切面速度值DMD 模態分析Fig. 15 DMD modal of section of vortex axis

6 飛翼標模自由搖滾特性

為研究飛翼標模在Ma=0.6 下的單自由度橫向失穩運動，采用風洞試驗與數值模擬結合的方法進行研究。風洞試驗是在中國航天空氣動力技術研究院FD-12 風洞中進行，飛翼標模橫向無量綱轉動慣量I*=36.5，無量綱方法為

式中：ρ=1.15 kg/m3為來流密度；A=1.2 為展弦比。

流動/運動耦合策略采用緊耦合方法，并引入了“預估-校正”機制，在雙時間步的內迭代過程中，流動與運動同時收斂，提高系統穩定性。

六自由度運動方程可采用與流動控制方程相同的時間離散方式，根據線性多步法松弛迭代求解，松/緊耦合可用統一的形式描述［30］：

式中：Q6DOF為運動學參數；F6DOF為飛行力學控制方程通量；g、h、ω為松弛因子；a、b、c、d、e為系數。選取不同的系數即可得到不同的時間精度和耦合方式的計算方法，見表2，本文采用的方法為表2 中 的BDF2（2nd-Backward-Difference Formula）方法。

表2 飛行力學離散控制方程系數Table 2 Coefficients of discrete governing equations of flight mechanics

時間步長對自由搖滾準確捕捉也有一定的影響。采用dt=2.25×10-4，4.57×10-4，6.30×10-4，9.00×10-4s 這4種時間步長對Ma=0.6、α=24.5°工況進行測試，圖16 為時間步長對自由搖滾時間歷程的影響，從圖中可以看到，不同的時間步長對搖滾歷程頻率影響很小，對局部滾轉角峰值有一定影響。隨著時間步長不斷減小，當時間步長減小至dt=4.57×10-4s 時，其滾轉角峰值幾乎不隨時間步長變化，證明時間步長dt=4.57×10-4s 已滿足計算要求。

圖16 不同時間步長對搖滾歷程影響Fig. 16 Effects of time step on wing rock

圖17（a）為飛翼標模自由滾轉風洞試驗滾轉時間歷程；圖17（b）為滾轉角平均振幅隨俯仰角的變化，俯仰角運動范圍為22°～24.5°，圖中顯示，當標模在0 s時刻放開后，飛翼標模平衡滾轉角迅速調整至10°左右，滾轉角存在輕微抖動。隨著俯仰角逐漸增加，當俯仰角增加至24°時，平衡滾轉角繼續增大至27°附近，當俯仰角達到24.5°時，標模平衡滾轉角增大至28°，并且其橫向振蕩振幅突然增大，出現了自由搖滾運動。所以飛翼標模出現的搖滾運動為平衡位置偏向一側（平衡滾轉角為±28°左右）極限環運動，這與三角翼搖滾運動有一定的不同。

圖17Ma=0.6 自由搖滾試驗結果Fig. 17 Wind tunnel test of wing rock （Ma=0.6）

通過CFD/RBD 方法，同樣捕捉到了飛翼標模的自由搖滾運動及CFD 方法獲得的飛翼。圖18（a）為滾轉角變化時間歷程計算結果與試驗結果的對比，圖18（b）為二者的頻譜對比。圖中顯示，計算結果與試驗結果獲得的滾轉角時間歷程的平衡位置均為28°左右，計算結果振幅與試驗結果相比偏大，計算獲得自由搖滾平均振幅為12°左右，風洞試驗獲得平均振幅為8°左右。從頻譜分析可以看到，試驗搖滾減縮頻率k（無量綱化k=fc/（2V∞））為0.008 1，計算獲得的搖滾減縮頻率k=0.007 6，計算與試驗獲得搖滾頻率較為接近。

圖18 自由搖滾試驗與計算結果對比（α=26°）Fig. 18 Comparison between wing rock experiment and CFD results （α=26°）

圖19是Ma=0.6 時，在攻角α=22°，24.5°，26°下的單自由度運動特性曲線，左圖是滾轉角隨時間的變化曲線，右圖是相平面曲線，從圖中可以看到，當攻角α=22°時，運動規律是收斂的，隨著時間增加，滾轉逐漸減弱，從相平面中也能夠觀察出這樣的規律，當攻角α=24.5°，26°時，得到了典型的極限環運動，振蕩平衡滾轉角約為28°左右，振幅約為12°，計算獲得的自由搖滾整體規律與風洞試驗較為吻合。

飛翼標模動態失穩運動主要由模型在有滾轉情況下，迎風側與背風側復雜流動相互作用誘導而成。其中迎風側與背風側定義為：當模型在有俯仰角的前提下，滾轉角度為正時，在駕駛艙視角，右側前緣為迎風側，左側前緣為背風側。

圖19Ma=0.6 單自由度自由滾轉運動規律隨攻角變化Fig. 19 Wing rock characteristics with angle of attack（Ma=0.6）

飛翼標模搖滾機理分析選用的典型工況為Ma=0.6、α=24.5°，選用如圖19（b）紅色方框中的一個周期來進行分析，此周期內振幅已較為穩定。圖20 是自由搖滾過程中的空間流動變化，從圖中可以看到，飛翼標模在背風側一直保持著集中渦流動，同時在頭尖處還會誘導一個較小的頭部旋渦，在迎風側則體現為較為復雜的大分離流動，其集中渦在靠近頭部的位置就已經破裂。隨著滾轉角逐漸增加，背風一側集中渦逐漸增強，渦軸逐漸向內側移動，使得旋渦對飛翼標模的影響區域變大，而當滾轉角逐漸變小，渦軸逐漸向外移動，體現出了周期性運動。在迎風側，復雜的大分離流動體現為旋渦破裂后的螺旋波與高頻小尺度渦，在滾轉角較小時，高頻小尺度渦離物面較近，誘導能力較強，隨著滾轉角逐漸增加，小尺度渦逐漸離開物面，旋渦誘導能力減弱。

圖20 飛翼標模搖滾狀態下非定常流動變化Fig. 20 Change of unsteady flow during wing rock of flying-wing

圖21是亞聲速（Ma=0.6，α=24.5°）滾轉力矩隨滾轉角變化曲線，從圖中可以看到，在一個振蕩周期內，呈現出一個雙“8”字環的特征，其中以滾轉角28°為平衡位置，在滾轉角22.5°～35°區間內，遲滯環為順時針，整體效果為負阻尼效果，使得滾轉運動處于發散狀態，在滾轉角16°～22.5°區間與35°～42.5°范圍內，遲滯環為逆時針，整體效果為正阻尼效果，使得滾轉運動處于收斂狀態，順時針遲滯環面積與逆時針遲滯環面積基本相等，使得飛翼標模能夠在16°～42.5°范圍內，呈現出周期性的極限環振蕩運動規律。

圖21 滾轉力矩系數隨滾轉角變化Fig. 21 Roll moment coefficient with roll angle

對自由搖滾的流動機理分析主要可以從兩個方面進行研究，一方面是橫向靜態特性，另一方面是橫向動態效應。橫向靜態特性體現為，在自由搖滾范圍內（圖21），滾轉力矩氣動特征主要呈現靜穩定特性，使得標模在自由搖滾過程中滾轉角夠回到平衡位置，此氣動特性與三角翼搖滾較為一致。為了進一步分析自由搖滾流動機理，先分析靜態特性對失穩運動的影響。

圖22是不同滾轉角下的橫截面流線圖與壓力系數云圖，截面位置為X=400 mm。選取的是在飛翼標模自由搖滾歷程中所經歷的3 個主要滾轉角φ=20.4°，27.4°，38.5°。從圖中也能夠清晰的看到，飛翼標模橫截面主要流動體現為左側的集中渦流動與右側的小尺度旋渦流動，隨著滾轉角的變化，兩種流動也產生了相應的變化。圖23是不同滾轉角下的壓力系數分布。圖中顯示，模型背風面的吸力整體量值較大，迎風面的量值較小，說明對于飛翼標模滾轉力矩的影響主要來自于背風面，其中集中渦流動誘導正滾轉力矩，分離流動誘導負滾轉力矩。當φ=20.4°時，左側集中渦產生的吸力峰值較大，右側的旋渦分離流動吸力相對較小，所以對標模產生正向滾轉力矩。隨著滾轉角逐漸增大，左側集中渦位逐漸向內靠攏，渦吸力逐漸減小，當滾轉達到φ=27.4°時，左右側流動對滾轉力矩的貢獻較為接近，滾轉力矩接近為0。當φ=38.5°時，集中渦誘導的吸力很小，此時分離流動誘導的負滾轉力矩為主導，同時右側迎風面也會出現壓力峰值，所以整體產生負滾轉力矩，通過以上3 個典型滾轉角的分析能看到，飛翼標模體現出的滾轉靜穩定特性，主要為集中渦吸力誘導的正滾轉力矩與分離流誘導的負滾轉力矩隨著滾轉角變化相互“博弈”的結果。

圖22 不同滾轉角下的橫截面流線Fig. 22 Streamlines of cross section at different roll angles

圖23 不同滾轉角下的橫截面壓力分布對比Fig. 23 Pressure distribution of cross section at different roll angles

橫向動態效應體現為，在運動的過程中存在運動的遲滯效應（見圖21 中的雙“8”字遲滯環），使得飛翼標模能夠在平衡位置附近（φ=28°）耗散能量，在雙“8”字環的兩側積累能量，使得飛翼標模能夠產生極限環振蕩運動形式。以下分別選取雙“8”字環中的“負阻尼”（φ=28°附近）與“正阻尼”（φ=38.5°附近）瞬時壓力分布來分析流動對動態效應的貢獻。 圖24 為平衡位置附近（φ=28°）的運動歷程放大圖，其中有上行和下行兩個過程，圖25 的壓力分布分別對應上行與下行運動經過φ=28°時的瞬時壓力分布，圖26 為上行與下行瞬時流動特性，從圖中可以看到，當標模下行經過平衡位置（φ=28°）時，分離流動誘導的吸力較小，集中渦誘導的吸力較大，從而產生正向滾轉力矩，當飛翼標模上行經過平衡位置（φ=28°）時，分離流誘導的吸力較大，同時集中渦誘導的吸力有減小的現象，從而誘導出負向滾轉力矩。

圖24 雙“8”字環局部放大圖（φ=28°附近）Fig. 24 Zoom-in view of double “8” ring （near φ=28°）

圖25 上行與下行運動至φ=28°瞬時壓力分布Fig. 25 Instantaneous pressure distribution when up and down motion to φ=28°

圖26 上行與下行瞬時流動特性（φ=28°）Fig. 26 Transient flow characteristics of up and down motion （φ=28°）

圖27為雙“8”字環在滾轉角φ=38.5°附近的局部放大圖，圖28 為上行與下行運動至滾轉角φ=38.5°附近的壓力分布對比情況，圖29 為上行與下行瞬時流動特性。從圖中可以看到，在上行與下行運動至φ=38.5°時，左側集中渦變化不大，右側分離流動在靠近邊緣的地方有一定變化，在下行瞬時，分離流動誘導的吸力更大，并且在迎風側也誘導出較大的正壓力，二者均會誘導出負滾轉力矩，所以下行瞬時會比上行瞬時誘導出更大的負滾轉力矩。

圖27 雙“8”字環局部放大圖（φ=38.5°附近）Fig. 27 Zoom-in view of double “8” ring （near φ =38.5°）

圖28 上行與下行運動至φ=38.5°瞬時壓力分布Fig. 28 Instantaneous pressure distribution when up and down motion to φ=38.5°

圖29 上行與下行瞬時流動特性（φ=38.5°）Fig. 29 Transient flow characteristics of up and down motion （φ=38.5°）

7 飛翼標模橫向失穩舵面控制特性

飛翼標模動態失穩是一個復雜的非線性現象，為了嘗試探索飛翼標模自由搖滾的控制方法，采用偏轉上擾流板、外側副翼兩種手段來模擬控制面固定舵偏對自由搖滾現象的控制效果。首先采用了上擾流板偏轉方法，兩側上擾流板向上偏轉30°，其物面網格與局部網格見圖30。

圖30 上擾流板偏轉30°物面網格Fig. 30 Surface mesh of upper spoiler deflects by 30°

圖31 是亞聲速下自由搖滾上擾流板張開與無控狀態下的對比情況，從圖中可以看到，上擾流板張開之后，對于運動規律影響不大，在振幅上略有減小，在相平面中可以看到對于滾轉角速度影響也較小。可見上擾流板在此工況下幾乎無控制效果。圖32 是上擾流板控制與無控狀態下的頻譜對比，圖中顯示上擾流板張開對于幅值有一定影響，但是對于自由搖滾頻率影響微弱，并且依然存在機翼搖滾運動。

圖31 自由搖滾上擾流板控制與無控狀態對比Fig. 31 Upper spoiler control vs uncontrolled state in wing rock

圖32 上擾流板與無控狀態頻譜對比Fig. 32 Spectrum comparison between upper spoiler and uncontrolled state

圖33是外側副翼偏轉物面網格與局部網格，偏轉方式為兩側副翼向下偏轉30°。

圖33 外側副翼向下偏轉30°物面網格Fig. 33 Surface mesh of outer aileron deflects downward 30°

圖34為亞聲速工況下，自由搖滾外側副翼偏轉控制效果，從圖中可以看到，當外側副翼向下偏轉后，自由搖滾運動模式受到了很大的抑制效果，平衡位置從搖滾時的30°減小到25°左右，并且幾乎不存在周期性振蕩運動形式，外側副翼向下偏轉效果明顯。

圖34 自由搖滾外側副翼與無控狀態對比Fig. 34 Outer aileron vs uncontrolled state

8 結 論

采用非定常數值模擬DDES 方法模擬飛翼標模單自由度自由搖滾運動，并采用動力學模態分解（DMD）、物面脈動壓力分析等方法分析飛翼標模非定常特性、自由搖滾機理，得到以下初步結論：

1）飛翼標模大攻角下非定常流動特性主要體現為，頭部發展的集中渦、渦破裂、螺旋波流動結構，其中旋渦破裂點以St=0.12～0.23 的頻率沿渦軸振蕩，螺旋波頻率在St=1.16～2.33 范圍內。

2）飛翼標模亞聲速自由搖滾現象平衡位置在滾轉角28°，滾轉角運動范圍在16°～42.5°之間，在滾轉角22.5°～35°區間內，阻尼整體效果為負阻尼效果，使得滾轉運動處于發散狀態，在滾轉角16°～22.5°區間與35°～42.5°范圍內，整體效果為正阻尼效果，使得滾轉運動處于收斂狀態。

3）自由搖滾現象的流動機理主要是背風側的集中流動與迎風側的具有較多小尺度旋渦的分離流動相互“博弈”。

4）上擾流板打開30°，對亞聲速自由搖滾控制效果不明顯，外側副翼向下打開30°，能夠在一定程度上抑制自由搖滾現象。