超聲速細長飛行器前置小翼展開過程的非定常特性

唐海敏，杜廈，傅建明，李欣益

上海機電工程研究所，上海 201109

現代戰爭中，攔截用飛行器是攔截、攻擊空中移動目標的重要武器。隨著新技術革命浪潮的涌入，空中目標的速度越來越大、機動能力愈來愈高，攔截用飛行器的作戰性能隨之面臨新的挑戰：一方面需要具有較快的飛行速度，另一方面需要其在飛行末端命中目標前很短時間內具有快速改變飛行軌跡和飛行姿態的能力，即在該時段內具有非常強的機動性，以提高攔截成功率。

為應對上述新的性能要求，日前主要采用姿態控制側向直接力方案和軌跡控制側向直接力方案。但這兩種方案因涉及側向力發動機布置使系統和結構更加復雜，側向噴流與來流相互干擾產生的復雜波渦結構[1-6]使飛行器總體氣動特性甚至舵面效率變化隨參數呈現更強的非線性[4]，側向力發動機燃料受限使控制時間較短，增加了側噴控制方案的設計和使用難度。為此，考慮在飛行器頭部附近加裝可展開的前置小翼，在需要時前置小翼快速伸出產生氣動力及控制力矩，實現對飛行器姿態和軌跡的控制。前置小翼伸出之后，飛行器的法向力系數增加，提高了氣動過載，壓心前移，降低了靜穩定度，從而使飛行器的機動性增強。

小翼在展開過程中繞飛行器的流動及飛行器本身所受氣動力都是復雜的非定常問題。翼面展開改變了飛行器的表面形狀，高速條件下會誘導出復雜的波系與渦系結構，如翼前激波、翼后分離渦等，同時其誘導產生的下洗流場影響飛行器彈體及后體尾舵的流動特性和載荷，這些流動機理以及翼面展開過程中飛行器氣動特性的時變規律，亟需開展詳細研究。

本文研究重點是前置小翼在展開過程中波渦等復雜流動結構的形成和發展以及氣動特性隨時間的變化規律。利用適合翼面展開大變形條件的結構動網格技術以及非定常脫體渦模擬(DES)方法，數值研究了大長細比飛行器上前置小翼展開對飛行器氣動特性的影響規律及機理，獲得了詳細的瞬時流場結構特性，分析了法向力、縱向壓心系數隨時間的變化特性。結果表明小翼展開對飛行器超聲速、大迎角條件控制力和力矩產生了較大的有利增益，是一種提高飛行器機動能力的有效途徑。

1 研究外形及條件

圖1給出了本文的研究外形，由彈體、尾舵以及4片前置小翼組成，前置小翼為扇形，展開時繞根弦前緣點旋轉，最大展開角度為30°，展開角速度為660 (°)/s。圖2給出了小翼展開角度φ的定義。計算條件由表1給出，其中滾轉角Φ的定義見圖3。

圖1 模型外形Fig.1 Model geometry

圖2 小翼展開角度定義Fig.2 Definition of wing spread angle

表1 來流條件Table 1 Incoming flow conditions

圖3 滾轉角定義(后視圖)Fig.3 Definition of rolling angle (rear view)

2 數值模擬方法

前置小翼展開將誘導產生復雜的非定常波系與渦系結構，準確模擬這種強非定常現象將對數值方法的模擬精度提出更高要求。傳統雷諾平均Navier-Stokes(RANS)方法難以準確模擬大范圍分離等復雜流動現象，直接數值模擬(DNS)或大渦模擬(LES)方法[7]雖能精確預測但因對計算機硬件能力需求過高仍難以承受，目前可行的途徑是采用RANS/LES混合方法[8-15]。RANS/LES混合方法結合了RANS和LES兩者各自的優勢，在邊界層采用RANS方程，在分離區采用LES，計算量比LES少，計算結果比RANS優，效果與LES相當。在當前計算機能力仍較難滿足DNS和LES的情況下，混合方法無疑是兼顧效果和效率的選擇，在復雜流動問題上，已被證實有很高的可信度[16]。本文采用Spalart等[8]提出的脫體渦模擬，并引入邊界層延遲轉捩函數抑制模型應力損失。

2.1 空間/時間離散方法

采用三維可壓縮Navier-Stokes方程為控制方程，通過有限體積方法進行離散，對流項離散采用二階精度的Roe格式[17]，黏性項離散采用中心差分格式，非定常數值方法時間項離散使用高計算效率的雙時間步方法[18-19]，內迭代采用LU-SGS[20]隱式方法，并采用局部時間步長方法加速內迭代過程收斂。

2.2 DES模型

DES模型最早由Spalart等[8]于1997年提出，在近壁面采用雷諾平均方法，用湍流模型模擬小尺度脈動，在其他區域，通過對傳統湍流模型尺度系數的修正，使模型起到LES亞格子應力模型的作用。DES方法結合了常規RANS方法和LES方法的優點，不僅可以節省時間，而且計算結果準確，尤其是對于大范圍分離流動的模擬。

根據其中RANS采用的湍流模型的不同，目前主要有基于Spalart-Allmaras(S-A)[21]湍流模型以及基于Menter剪切應力輸運(SST)湍流模型的DES方法[8-15]。

基于S-A模型的DES方法中，湍流渦黏性系數定義為

(1)

式中：

(2)

(3)

(4)

g=r+Cw2(r2-r)

(5)

(6)

式中：d為流場中空間點到壁面的最近距離；Δ為當地網格單元中心到臨近單元中心的最大距離；CDES一般取值0.65。

DES模型對網格分布具有較強依賴性，導致其往往會在邊界層內提前開啟LES模式。邊界層內采用LES模式容易導致計算得到的渦黏性不足，造成非物理的網格誘導分離。為解決這一問題，引入了邊界層延遲轉捩函數fd[15]，重新構造了混合方法長度尺度，可以在邊界層內關閉LES模式，這種方法稱為延遲脫體渦模擬(DDES)方法。

邊界層延遲轉捩函數fd的定義為

fd≡1-tanh((8rd)3)

(7)

式中：

(8)

其中：νt為湍流運動渦黏性系數；Ui,j為速度梯度張量；DDES的長度尺度為

(9)

由式(9)可以看出，當邊界層延遲轉捩函數fd=0時，DDES還原為S-A模型，當fd=1時，DDES模型與DES模型相同。由于fd在對數律層為1，到邊界層邊界迅速衰減為零，這樣保證了在邊界層內進行RANS計算時，同時起到了LES區與RANS區光滑過渡的作用。

2.3 網格生成方法

目前動網格技術較多采用基于變形與重構的非結構動態網格生成方法實現，但是非結構網格計算精度較低，同時網格重構又會帶來新的插值誤差。由于前置小翼展開過程外形變化較大，流場拓撲復雜，為了在保持原有網格拓撲結構不變的情況下，實現小翼展開過程的氣動網格變化，避免出現網格扭曲及插值帶來的誤差問題，考慮采用結構化網格，利用基于弧長的超限插值(TFI)方法實現網格跟隨變形邊界的運動，較好地保證了內部網格質量，特別是邊界層內網格質量。同時，采用文獻[22-23]中的體積修正法解決幾何守恒律問題，有效避免了網格變形誘導產生的非物理解。圖4(a)給出了翼面伸出后的局部網格。

在網格生成過程中，首先根據小翼展開的區域及特點，采用分區對接方式劃分初始網格，并預留出變形區域，如圖4(b)虛線框所示。由于翼面展開的軌跡及速度已知，在不同時刻計算出翼面角點及棱線相對上一時刻的變形量，通過二維TFI方法得到翼面網格的變形量，再利用三維TFI方法得到與翼面相關的網格塊內點擾動量，將擾動量與上一時刻網格相加得到本時刻變形后的網格，圖4(c)給出了典型時刻變形后的網格。

圖4 典型時刻結構動網格Fig.4 Moving structure grids at typical times

2.4 邊界條件

1) 入流邊界：來流為超聲速，入流邊界處給定來流的靜壓、靜溫及馬赫數。

2) 出流邊界：出口為超聲速時，下游流場不影響上游流場，將所有參數數值外推。

3) 物面邊界：壁面速度滿足無滑移條件，壁面壓力滿足法向動量方程，壁面密度利用絕熱壁條件和完全氣體狀態方程確定[22]。

2.5 數值模擬方法驗證

根據研究外形及問題特點，采用尖拱形圓柱體模型以及Basic Finner模型分別驗證建立的高精度RANS/LES混合數值模擬方法，驗證模型及試驗條件具體見文獻[24-25]。圖5和表2給出了尖拱形圓柱體模型的外形示意圖與計算條件，其中φ為橫截面子午角。圖6給出了計算得到的尖拱形圓柱體模型中心線壓力分布計算與試驗對比情況，可以看出，計算與試驗結果吻合較好。

數字化檔案的文種形式、載體形式、表達形式都具有多樣性，因此數字化檔案信息具有多源異構的特點，要建立統一的檔案信息關聯模型，需要分別采用不同的特征提取方法對多種類型的檔案進行特征提取，其中特征信息可能包括時間、人物、地點、類型、保管期限、著錄機構等。有的特征信息可以通過著錄項獲取，有的特征信息需要使用模式識別等人工智能技術從檔案內容中提取，這些特征信息都是檔案的屬性數據，不同的屬性數據分屬不同的特征視圖，將提取到的特征按照不同的視圖分別構建特征關聯圖，再使用數據融合技術融合多個視圖下的特征關聯圖，最終得到檔案信息關聯模型。為了描述方便，對檔案關聯模型中的基本元素進行如下定義：

圖7和表3分別給出了Basic Finner模型的外形示意圖與計算條件，其中D為直徑。圖8給出了Basic Finner外形動態特性計算與試驗對比情況，并與文獻結果[26]進行了比對，可以看出，采用本文DES方法計算得到的結果與試驗趨勢相同，并優于文獻[26]的結果。

圖5 尖拱形圓柱體模型Fig.5 Model for ogive-circular cylinder body

表2尖拱形圓柱體模型的計算條件

Table2Calculationconditionsofmodelforogive-circularcylinderbody

馬赫數迎角/(°)靜壓/Pa靜溫/K雷諾數1．23041893．31248．44721．40×107

圖6 中心線壓力分布計算與試驗對比Fig.6 Comparison of calculation and experiment results of pressure distribution at central lines (Ma∞=1.2，α=30°)

圖7 Basic Finner 標模外形Fig.7 Model shape for Basic Finner

表3 Basic Finner模型計算條件Table 3 Calculation conditions of Basic Finner model

圖8 Basic Finner外形俯仰阻尼導數計算與試驗對比(Ma∞=1.96,Re=6.233×106)Fig.8 Comparison of calculation and experiment results of pitch damping derivatives of Basic Finner model (Ma∞=1.96,Re=6.233×106)

3 結果與分析

本文采用結構動網格技術以及非定常DES方法，就圖1細長飛行器外形開展了前置小翼展開過程非定常數值仿真研究和分析。

3.1 小翼展開過程定常與非定常流場特性對比

圖9給出了Ma∞=2.0、α=15°條件下4片前置小翼以相同角速度(660 (°)/s)同時展開過程中典型時刻相對初始狀態的壓力差(p-p初值)/p∞云圖，圖中前置小翼展開角度分別為φ=10°、20°、30°。可以看出，小翼展開過程中，翼面附近流場發生顯著變化，翼后近區迎風面高壓沿小翼后緣向側面發展，背風面高壓同時也向側面加強和發展，最后背風面和迎風面的兩片高壓區連通,呈“X”形分布。同時可觀察到，小翼展開對彈體中段和尾段也有較明顯的干擾影響，在載荷計算時需要特別關注。

圖9 小翼展開過程相對初始狀態的表面壓力差云圖Fig.9 Contours of difference in surface pressure relative to the initial state during the spreading of wings (Φ=45°,Ma∞=2.0,α=15°)

為直觀地分析小翼展開過程飛行器表面的非定常流場特性，圖10給出了非定常與定常表面壓力差(p非定常-p定常)/p∞云圖以及翼面和尾舵側面局部放大圖。從圖中可見，非定常效應主要集中在翼面及尾舵附近。一方面，隨著翼面的展開，其附近非定常效應從背風面向迎風面發展，對上翼面產生的法向力和抬頭力矩有所抵消；同時，翼后圓柱體上的高壓區也從側面向迎風面發展，產生附加的法向力和抬頭力矩，兩者互有消長。

圖10 小翼展開過程非定常與定常表面壓力差云圖(Φ=45°,Ma∞=2.0,α=15°)Fig.10 Contours of difference in surface pressure of the unsteady and steady solutions during spreading of wings (Φ=45°,Ma∞=2.0,α=15°)

另一方面，非定常方法真實模擬了小翼對尾舵的洗流影響，尾舵區表現出明顯的非定常效應。由尾舵附近局部放大圖可見，背風面和迎風面云圖結構類似，略有差別。主要特征是上下對稱，在φ=20°時幾乎完全對稱，而φ=10°和φ=30°則有不同，上舵面的高壓產生附加的法向力和低頭力矩，而其根部圓柱體上的高壓則使法向力減小并產生附加的抬頭力矩。非定常效應的集中體現最終取決于小翼展開過程中以上各部分貢獻。

為了進一步分析小翼展開過程非定常效應對空間流場的影響，圖11給出了小翼附近(X=1 500 mm)、中距(X=2 400 mm)和遠距(X=3 500 mm)的截面橫向流。從圖中可以看出，隨著流場向下游發展，渦的結構尺度逐漸變大，渦核有遠離物面的趨勢，同時渦的形狀從類圓形漸變成類橢圓形。還可看出，隨著小翼的展開，渦的結構開始變得復雜，出現雙渦結構，表明流場中出現了渦的分裂與合并。還可進一步看出，非定常與定常流場結構基本一致，在物面附近卷起典型的渦旋，但渦的尺度、渦核位置有所不同。仔細觀察，可以發現在小翼附近和翼后中距處一致性較高，差異最明顯的是翼后遠距處，說明前置小翼在展開過程中與尾舵之間存在明顯的遠距干擾現象。

3.2 氣動特性隨小翼展開的變化

圖12和圖13分別給出了小翼展開過程中飛行器法向力系數和縱向壓心系數的變化歷程。可以看出，小翼展開過程中展開角度φ<10°時法向力系數變化較為平緩，展開角度φ≥10°時基本上呈線性增加趨勢，縱向壓心系數隨小翼展開角度的增加而減小。對比圖中定常和非定常結果，可以看出在展開角度φ=10°時，定常與非定常結果存在明顯差異。表4詳細給出了初始狀態和小翼完全展開后的法向力和縱向壓心系數對比情況，可以看出小翼完全展開后，法向力增加15.6%，提高了氣動過載；壓心前移4.1%，降低了飛行器的靜穩定性，從而改善了飛行器的機動性能。

圖11 小翼展開過程X=1 500, 2 400,3 500 mm截面橫向流Fig.11 Crossflow at X=1 500, 2 400,3 500 mm cross-section during spreading of wings

圖12 小翼展開過程法向力系數的變化Fig.12 Variation in normal force coefficient during spreading of wings (Ma∞=2,α=15°,Φ=45°)

圖13 小翼展開過程縱向壓心系數的變化Fig.13 Variation in center of pressure coefficient during spreading of wings (Ma∞=2,α=15°,Φ=45°)

表4小翼展開前后法向力和縱向壓心系數結果對比

Table4Comparisonofnormalforceandcenterofpressurecoefficientsbeforeandafterspreadingofwings

系數展開前展開后差量百分比法向力3．0463．5210．47515．6縱向壓心0．523870．48250-0．041374．1

4 結 論

利用結構動網格技術和基于脫體渦(DES)的非定常數值模擬方法，研究了細長飛行器前置小翼展開過程的非定常效應。文中給出了小翼展開過程的詳細動態流場特性，并分析研究了法向力、縱向壓心系數隨小翼展開過程的變化規律。

1) 采用的數值方法計算結果與試驗吻合良好，方法具有較高的可信度與適用性。

2) 超聲速、大迎角條件下前置小翼展開過程非定常效應顯著。

3) 小翼展開過程對小翼附近區域以及尾舵區域產生了強烈的影響。

4) 小翼展開后壓心前移，降低了飛行器的靜穩定性，法向力系數增大，提高了氣動過載，對飛行器機動性能都產生了有利影響。

參 考 文 獻

[1] CHAMBERLAIN R, MCCLURE D, DANG A. CFD analysis of lateral jet interaction phenomena for the THAAD interceptor: AIAA-2000-0963[R]. Reston, VA: AIAA, 2000.

[2] DASH S M, PERRELL E R, ARUNAJATESAN S, et al. Lateral jet aerodynamic interaction simulations for dynamic pressure loads: AIAA-2000-2036[R]. Reston, VA: AIAA, 2000.

[3] EBRAHIMI H B. Numerical investigation of jet interaction in a supersonic freestream: AIAA-2005-4866[R]. Reston, VA: AIAA, 2005.

[4] 劉耀峰. CFD軟件研發和戰術導彈噴流/繞流復雜流動數值模擬[D]. 北京: 北京理工大學, 2005: 56-88.

LIU Y F. Development of CFD software and complex flows numerical simulation for interaction flow field between external flow and jet flow of tactical missiles[D]. Beijing: Beijing Institute of Technology, 2005: 56-88 (in Chinese).

[5] CUI J, AGARWAL R K. 3-D CFD validation of an axisymmetric jet in cross-flow: AIAA-2005-1112[R]. Reston, VA: AIAA, 2005.

[6] RAUSCH J R, ROBERTS B B. Reaction control system plume flow field interaction effects on the space shutter orbiter: AIAA-1974-1104[R]. Reston, VA: AIAA, 1974.

[7] YANG G, YAO Y F, FANG J, et al. Large-eddy simulation of shock-wave/turbulent boundary layer interaction with and without SparkJet control[J]. Chinese Journal of Aeronautics, 2016, 29(3): 617-629.

[8] SPALART P R, JOU W H, STRELETS M, et al. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[C]∥LIU C Q, LIU Z N. Advances in DNS/LES, 1st AFOSR International Conference on DNS/LES. Los Angles: Greyden Press, 1997: 137-147.

[9] STRELETS M. Detached eddy simulation of massively separated flows: AIAA-2001-0879[R]. Reston, VA: AIAA, 2001.

[10] 肖志祥, 符松. 用RANS/LES混合方法研究超音速底部流動[J].計算物理, 2009, 26(2): 221-230.

XIAO Z X, FU S. Study on supersonic base flow using RANS/LES methods[J]. Chinese Journal of Computational Physics, 2009, 26(2): 221-230 (in Chinese).

[11] SOSHI K, KOZO F. Computational study of a supersonic base flow using LES/RANS hybrid methodology: AIAA-2004-0068[R]. Reston, VA: AIAA, 2004.

[12] 趙偉文, 萬德成. 用SST-DES和SST-URANS方法數值模擬亞臨界雷諾數下三維圓柱繞流問題[J]. 水動力學研究與進展, 2016, 31(1): 1-8.

ZHAO W W, WAN D C. Numerical study of 3D flow past a circular cylinder at subcritical Reynolds number using SST-DES and SST-URANS[J]. Chinese Journal of Hydrodynamics, 2016, 31(1): 1-8 (in Chinese).

[13] 孫明波, 梁劍寒, 王振國. 超音速來流橫向狹縫噴流標量輸運的混合RANS/LES模擬[J]. 力學季刊, 2007, 28(3): 395-399.

SUN M B, LIANG J H, WANG Z G. Hybrid RANS/LES simulation for scalar transport of slot injection into a supersonic stream[J]. Chinese Quarterly of Mechanics, 2007, 28(3): 395-399 (in Chinese).

[14] WON S H, JEUNG I S, CHOI J Y. DES study of transverse jet injection into supersonic cross flows: AIAA-2006-1227[R]. Reston, VA: AIAA, 2006.

[15] SHIN J R, CHOI J Y. DES study of base and base-bleed flows with dynamic formulation of DES constant: AIAA-2011-0662[R]. Reston, VA: AIAA, 2011.

[16] 周鑄, 黃江濤, 黃勇, 等. CFD技術在航空工程領域的應用、挑戰與進展[J]. 航空學報, 2017, 38(3): 020891.

ZHOU Z, HUANG J T, HUANG Y, et al. CFD technology in aeronautic engineering field: Applications, challenges and development[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(3): 020891 (in Chinese).

[17] ROE P L. Approximate Riemann solvers, parameter vectors and different schemes[J]. Journal of Computational Physics, 1981, 43(2): 357-372.

[18] 唐海敏. 三維帶操縱面機翼Euler方程非定常氣動力計算[D]. 南京: 南京航空航天大學, 2003: 12-13.

TANG H M. Aerodynamic computation of 3-D unsteady Euler equations for wings with control surface[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2003: 12-13 (in Chinese).

[19] JAMESON A. Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings: AIAA-1991-1596[R]. Reston, VA: AIAA, 1991.

[20] YOON S, JAMESON A. Lower-upper symmetric-Gauss-Sediel method for the Euler and Navier-Stokers equations[J]. AIAA Journal, 1988, 26(9): 1025-1026.

[21] SPALART P R, ALLMARAS S R. A one-equation turbulence model for aerodynamic flows: AIAA-1992-0439[R]. Reston, VA: AIAA, 1992.

[22] 王巍. 有相對運動的多體分離過程非定常數值算法研究及試驗驗證[D]. 長沙: 國防科學技術大學, 2008: 29-31.

WANG W. Numerical simulation technique research and experiment verification for unsteady multi-body flowfield involving relative movement[D]. Changsha: National University of Defense Technology, 2008: 29-31 (in Chinese).

[23] 劉君, 劉瑜, 陳澤棟. 非結構變形網格和離散幾何守恒律[J]. 航空學報, 2016, 37(8): 2395-2407.

LIU J, LIU Y, CHEN Z D. Unstructured deforming mesh and discrete geometric conservation law[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2395-2407 (in Chinese).

[24] BARCHE J. Experimental data base for computer program assessment[R]. London: Technical Editing and Reproduction Ltd, 1979.

[25] USELTON B L, USELTON J C. Test mechanism for measuring pitch damping derivatives of missile configurations at high angles of attack: AEDC-TR-75-43[R]. Tennessee: Arnold Engineering Development Center, 1975.

[26] SCOTT M M. A reduced-frequency approach for calculating dynamic derivatives: AIAA-2005-0840[R]. Reston，VA: AIAA, 2005.