孟德虹， 孫巖， 王運(yùn)濤,*， 李偉

1.中國空氣動力研究與發(fā)展中心 計(jì)算空氣動力研究所， 綿陽　621000 2.中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室， 綿陽　621000

?

戰(zhàn)斗機(jī)垂尾脈動壓力數(shù)值模擬

孟德虹1， 孫巖2， 王運(yùn)濤1,*， 李偉1

1.中國空氣動力研究與發(fā)展中心 計(jì)算空氣動力研究所， 綿陽621000 2.中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室， 綿陽621000

在亞跨超計(jì)算流體力學(xué)(CFD)軟件平臺(TRIP)上開發(fā)了基于RANS/LES混合思路的IDDES流動模擬技術(shù)，并通過NACA0021翼型60° 大迎角分離流動與串列圓柱繞流模擬對RANS/LES混合方法的精確度進(jìn)行了驗(yàn)證，針對某戰(zhàn)斗機(jī)外形的垂尾脈動壓力開展了數(shù)值模擬研究。戰(zhàn)斗機(jī)來流馬赫數(shù)為0.1，基于全機(jī)長度的雷諾數(shù)為2×106，模型迎角為20°、30°和40°。分別通過脈動壓力系數(shù)、脈動壓力功率譜密度、空間流動結(jié)構(gòu)以及側(cè)向力響應(yīng)曲線等對戰(zhàn)斗機(jī)的垂尾脈動壓力進(jìn)行了分析。脈動壓力模擬結(jié)果表明：當(dāng)垂尾完全沉浸在邊條翼脫體渦破碎后的寬頻湍流脈動氣流中時(shí)，垂尾翼梢位置的脈動壓力會發(fā)生明顯的增大。

戰(zhàn)斗機(jī)； 垂尾； 脈動壓力； RANS/LES； 計(jì)算流體力學(xué)(CFD)

為了提高機(jī)動性和操縱性，現(xiàn)代高性能戰(zhàn)斗機(jī)普遍采用邊條翼/垂直尾翼的氣動布局。大后掠角的邊條翼在前緣位置產(chǎn)生高強(qiáng)度的脫體渦，為主機(jī)翼上方的氣流提供能量，延遲主機(jī)翼上的流動分離，提高主機(jī)翼的升力。但隨著戰(zhàn)斗機(jī)迎角的增加，邊條翼產(chǎn)生的脫體渦會在垂尾前方發(fā)生破裂，使戰(zhàn)斗機(jī)垂尾沉浸在寬頻的湍流脈動中，從而誘導(dǎo)垂尾發(fā)生抖振[1-2]。垂尾抖振會引起結(jié)構(gòu)疲勞破壞，甚至導(dǎo)致飛行事故發(fā)生。戰(zhàn)斗機(jī)設(shè)計(jì)前，需要對垂尾抖振邊界與抖振響應(yīng)進(jìn)行準(zhǔn)確預(yù)測，避免戰(zhàn)斗機(jī)在飛行包線內(nèi)長時(shí)間處于抖振區(qū)域飛行，提高飛機(jī)的安全性能和使用壽命。大量的研究結(jié)果表明，垂尾抖振是由邊條渦破裂導(dǎo)致的垂尾表面脈動壓力引起的。因此，垂尾表面脈動壓力的預(yù)測是判斷抖振起始邊界、分析抖振響應(yīng)特性及研究抖振作用機(jī)理的關(guān)鍵[3-4]。

目前，垂尾表面脈動壓力預(yù)測主要有風(fēng)洞試驗(yàn)和數(shù)值計(jì)算兩種方法。風(fēng)洞試驗(yàn)方法[5]通過在垂尾模型表面安裝測壓孔，獲取垂尾表面的脈動壓力響應(yīng)。該方法能夠獲得物理真實(shí)的脈動壓力響應(yīng)特性，但脈動壓力測量點(diǎn)數(shù)量受限于垂尾模型結(jié)構(gòu)，且空間脈動壓力測量困難；另一方面，開展垂尾脈動壓力測量風(fēng)洞試驗(yàn)的時(shí)間周期較長、成本較高。數(shù)值計(jì)算方法通過模擬飛行器繞流下的流場，獲取垂尾表面的脈動壓力分布。當(dāng)前，工程應(yīng)用中通常采用離散求解雷諾平均Navier-Stokes方程RANS(Reynolds Average Navier-Stokes)的方式模擬飛行器的繞流流場，該方法具有計(jì)算量小、程序魯棒性好的特點(diǎn)，但該方法的湍流渦黏性值較大，會快速耗散小尺度高頻渦，不利于模擬具有多重尺度的寬頻流場結(jié)構(gòu)。因此，傳統(tǒng)RANS方法模擬戰(zhàn)斗機(jī)垂尾表面脈動壓力存在一定的局限性。大渦模擬(Large Eddy Simulation， LES)技術(shù)具有較好的湍流分辨能力，但計(jì)算量過大、壁面邊界條件處理困難，尚無法在工程中取得廣泛應(yīng)用。綜合RANS和LES優(yōu)缺點(diǎn)的分離渦模擬(Detached Eddy Simulation， DES)技術(shù)成為近十年流動數(shù)值模擬領(lǐng)域的研究熱點(diǎn)[6-11]，并成功應(yīng)用于大范圍分離流動數(shù)值模擬[12-15]。

為此，本文在亞跨超CFD計(jì)算軟件(TRIP)的框架內(nèi)發(fā)展了基于RANS/LES混合方法的模擬技術(shù)，通過兩個(gè)標(biāo)準(zhǔn)算例對RANS/LES混合方法進(jìn)行了驗(yàn)證。在此基礎(chǔ)上，利用RANS/LES混合方法研究了某戰(zhàn)斗機(jī)模型的垂尾脈動壓力特性。

1　數(shù)值計(jì)算方法

1.1流動控制方程求解

理想氣體假設(shè)下的非定常N-S方程可以描述為

(1)

式中：Q為守恒變量；Finv為無黏通量；Fvis為黏性通量；Ω為任意封閉域；?Ω為包圍Ω的邊界；n為邊界位置的法向矢量；t為物理時(shí)間。

本文通過有限體積方法離散求解式(1)，并基于TRIP軟件實(shí)現(xiàn)復(fù)雜流動的數(shù)值模擬。TRIP軟件是中國空氣動力研究與發(fā)展中心自主研發(fā)的基于多塊結(jié)構(gòu)對接/拼接網(wǎng)格大型通用計(jì)算流體力學(xué)(ComputationalFluidDynamics，CFD)計(jì)算軟件，廣泛應(yīng)用于多種工程問題的流動模擬，具有較高的計(jì)算精度和數(shù)據(jù)可信度[16-19]。

本文的數(shù)值計(jì)算中，式(1)的無黏項(xiàng)離散采用矩陣耗散的中心型差分格式，二階和四階人工耗散系數(shù)分別取0.000 1和0.01，黏性項(xiàng)離散采用二階中心格式。非定常計(jì)算采用雙時(shí)間方法進(jìn)行求解，偽時(shí)間步采用隱式方法以提高計(jì)算效率。湍流模型包含單方程SA模型[20]與兩方程SST模型[21]。采用低速預(yù)處理技術(shù)、多重網(wǎng)格技術(shù)和大規(guī)模并行技術(shù)提高計(jì)算效率。

1.2RANS/LES混合方法

RANS/LES混合方法的思想是在壁面區(qū)域采用RANS方法對湍流脈動進(jìn)行?；?，而在分離流動區(qū)域利用LES方法分辨多尺度的分離渦。

分離渦模擬方法DES是RANS/LES混合方法的一種典型模式，最早由Spalart[6]等提出，經(jīng)過十余年的發(fā)展，已經(jīng)形成了多種形式的分離渦模擬技術(shù)[7-8,10]。其中，Shur[10]等發(fā)展的IDDES(ImprovedDelayedDES)具有DDES模型與壁面LES模型的雙重優(yōu)點(diǎn)，能夠解決邊界層附近“Log-LayerMismatch”的問題，且能實(shí)現(xiàn)分離區(qū)中RANS到LES的快速轉(zhuǎn)換。本文采用的RANS/LES混合方法是在SST湍流模型基礎(chǔ)上發(fā)展的IDDES流動模擬技術(shù)。

SST湍流模型的控制方程可以表示為

(2)

(3)

式中：常數(shù)β*為0.09。

將式(2)中的湍流長度尺度LRANS替換為大渦模擬的濾波尺度LLES，則式(2)轉(zhuǎn)變?yōu)閮煞匠痰拇鬁u模擬控制方程。通過選擇長度尺度混合方式，可以得到不同的RANS/LES混合方法。IDDES方法的長度尺度LIDDES定義為

(4)

(5)

其中：CDES為DES經(jīng)驗(yàn)參數(shù)，取值為0.65，其他系數(shù)和變量的定義分別為

(6)

式中：dw為壁面距離；Δmax為網(wǎng)格單元最大尺度；Δwn為壁面網(wǎng)格單元法向高度；Δ為網(wǎng)格長度尺度；常數(shù)Cw為0.15；常數(shù)ct為1.87；常數(shù)cl為5.0；rdt、rdl分別為rd的湍流和層流分量，定義為

(7)

式中：νt、νl分別為湍流與層流運(yùn)動黏性系數(shù)；ui為流動速度分量；xj為空間坐標(biāo)分量；常數(shù)κ為0.41。

2　數(shù)值方法驗(yàn)證

本節(jié)通過NACA0021翼型與串列圓柱兩個(gè)分離流標(biāo)準(zhǔn)算例對本文選擇的RANS/LES混合方法進(jìn)行驗(yàn)證。

2.1NACA0021翼型繞流

NACA0021翼型大攻角繞流計(jì)算是歐盟DES標(biāo)準(zhǔn)算例之一。該算例計(jì)算條件為，來流Ma為0.1，基于弦長的雷諾數(shù)Re為270 000，翼型迎角α為60°。計(jì)算網(wǎng)格采用DESider官方網(wǎng)站提供的網(wǎng)格，圖1給出了展向剖面的網(wǎng)格，單個(gè)剖面網(wǎng)格單元數(shù)量為21 400，展向網(wǎng)格的尺度為0.025c，分別計(jì)算了展向長度Lz為c、2c和4c時(shí)的繞流流場，對應(yīng)的展向網(wǎng)格節(jié)點(diǎn)數(shù)量分別為41、81和161。

圖1　NACA0021翼型流場計(jì)算網(wǎng)格(X-Y截面)Fig.1　Flow field computational grid of NACA0021airfoil (X-Y section)

圖2　NACA0021翼型壓力系數(shù)分布Fig.2　Pressure coefficient distribution of NACA0021 airfoil

圖2給出了不同展向長度下NACA0021翼型的時(shí)均壓力系數(shù)分布曲線與試驗(yàn)結(jié)果的對比。其中，壓力系數(shù)采樣點(diǎn)為中間展向位置的網(wǎng)格點(diǎn)?？梢钥闯?，展向長度值對下翼面壓力系數(shù)分布影響很小，不同展向長度的下翼面壓力系數(shù)與試驗(yàn)值吻合較好，但對于上翼面壓力系數(shù)分布，展向長度的影響明顯。由于上翼面處于大的分離區(qū)內(nèi)，上翼面壓力趨向于一致，表現(xiàn)為一個(gè)明顯的壓力平臺區(qū)。展向長度Lz等于弦長c時(shí)，本文計(jì)算得到壓力平臺區(qū)存在一定的波動，與試驗(yàn)結(jié)果存在明顯的差異；當(dāng)Lz=2c時(shí)，計(jì)算得到的上翼面壓力平臺區(qū)明顯，并與試驗(yàn)結(jié)果比較接近；當(dāng)模型展向長度繼續(xù)增加，如Lz=4c時(shí)，壓力平臺區(qū)也比較明顯，但計(jì)算數(shù)據(jù)與試驗(yàn)結(jié)果之間的差異在逐漸增大。

圖3給出了展向長度Lz=2c時(shí)NACA0021翼型阻力系數(shù)功率譜密度計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)的對比，其中橫坐標(biāo)St為斯特勞哈數(shù)?？梢钥闯?，本文計(jì)算結(jié)果較好地分辨出了前后緣分離渦脫落的主頻率。

圖3　NACA0021翼型阻力系數(shù)功率譜密度曲線(Lz=2c)Fig.3　Power spectrum density curves of drag coefficients of NACA0021 airfoil (Lz=2c)

2.2串列圓柱

串列圓柱大范圍分離流動是歐盟ATAAC項(xiàng)目與美國AIAABANC項(xiàng)目的標(biāo)準(zhǔn)算例之一。該算例計(jì)算條件為，來流速度U為44m/s，基于圓柱直徑D的雷諾數(shù)ReD為166 000，兩個(gè)圓柱軸之間的距離L為圓柱直徑的3.7倍。本文計(jì)算網(wǎng)格采用ATAAC提供的標(biāo)準(zhǔn)網(wǎng)格，圖4給出了該網(wǎng)格的X-Y截面網(wǎng)格分布。X-Y截面網(wǎng)格單元數(shù)量約為80 000，圓柱展向長度Lz取圓柱直徑D的3.0倍，展向網(wǎng)格尺度取0.02D，展向共有151個(gè)網(wǎng)格點(diǎn)。

圖4　串列圓柱計(jì)算網(wǎng)格(X-Y截面)　　Fig.4　Computational grid of tandem cylinders(X-Y section)

圖5給出了串列圓柱前后圓柱時(shí)均壓力系數(shù)分布曲線與試驗(yàn)結(jié)果的對比，其中橫坐標(biāo)θ為圓柱表面點(diǎn)的方位角。從圖中可以看出，計(jì)算的串列圓柱壓力系數(shù)分布與風(fēng)洞試驗(yàn)結(jié)果基本吻合。

圖6給出了后圓柱給定位置點(diǎn)(θ=45°)的壓力系數(shù)脈動功率譜密度曲線計(jì)算與試驗(yàn)結(jié)果的對比，其中橫坐標(biāo)f為頻率?？梢钥闯觯疚挠?jì)算方法較好地捕捉到了后圓柱分離渦脫落的主頻率。

圖5　串列圓柱壓力系數(shù)分布Fig.5　Pressure coefficient distribution of tandem cylinders

圖6　串列圓柱功率譜密度曲線Fig.6　Power spectrum density curves of tandem cylinders

3　垂尾脈動壓力計(jì)算

利用前面經(jīng)過驗(yàn)證的RANS/LES混合計(jì)算方法，本節(jié)開展了某戰(zhàn)斗機(jī)風(fēng)洞模型垂尾脈動壓力特性的計(jì)算研究。圖7給出了該戰(zhàn)斗機(jī)模型的外形，該模型在機(jī)翼前方有一個(gè)大后掠的邊條，與戰(zhàn)斗機(jī)的進(jìn)氣道融合在一起，機(jī)翼尾部上方有兩個(gè)垂直尾翼。該模型在中國空氣動力研究與發(fā)展中心低速所進(jìn)行了常規(guī)測力試驗(yàn)，模型進(jìn)氣道和尾噴管進(jìn)行了封堵處理，沒有考慮進(jìn)氣道流動對模型繞流的影響。

圖7　戰(zhàn)斗機(jī)模型外形Fig.7　Shape of fighter model

3.1計(jì)算網(wǎng)格與計(jì)算參數(shù)

圖8給出了該戰(zhàn)斗機(jī)模型的流場計(jì)算網(wǎng)格。該網(wǎng)格為多塊結(jié)構(gòu)對接網(wǎng)格，包含310個(gè)網(wǎng)格塊、31 610 240個(gè)網(wǎng)格單元，第一層網(wǎng)格距離為5.0×10-6，對應(yīng)的網(wǎng)格y+約為1.0。來流馬赫數(shù)Ma為0.1，基于模型全機(jī)長度的雷諾數(shù)Re為2.0×106。根據(jù)Spalart的估算，DES采用計(jì)算網(wǎng)格需要達(dá)到107量級，本文計(jì)算網(wǎng)格規(guī)模達(dá)到DES計(jì)算的網(wǎng)格數(shù)量要求。網(wǎng)格計(jì)算域的最大范圍為戰(zhàn)斗機(jī)平均氣動弦長的100倍左右，在戰(zhàn)斗機(jī)背風(fēng)面分離區(qū)域進(jìn)行加密，最小的網(wǎng)格尺度小于平均氣動弦長的0.5%，并且在分離流動區(qū)域，盡量使得網(wǎng)格各向同性(即Δx≈Δy≈Δz)。無量綱非定常計(jì)算時(shí)間步長Δt取0.004，保證每個(gè)流動周期內(nèi)不少于100個(gè)迭代步。模型迎角α取20°、30°和40° 三個(gè)角度。

圖8　戰(zhàn)斗機(jī)模型計(jì)算網(wǎng)格Fig.8　Computational grid of fighter model

3.2垂尾脈動壓力特性

圖9給出了不同機(jī)身迎角下的全機(jī)等Q值圖?？梢钥闯觯S著機(jī)身迎角的增加，邊條翼形成的脫體渦在加速破裂，機(jī)身迎角為20° 時(shí)，邊條翼脫體渦一直延伸到垂尾的弦向中部位置才發(fā)生破裂，垂尾主要受脫體渦的影響；機(jī)身迎角為30° 時(shí)，邊條翼脫體渦在垂尾前附近發(fā)生破裂，整個(gè)垂尾沉浸在破裂后的小渦結(jié)構(gòu)中；機(jī)身迎角為40° 時(shí)，邊條翼脫體渦在主機(jī)翼前緣位置附近就發(fā)生破裂，破裂后的小渦結(jié)構(gòu)沿著流向運(yùn)動，到達(dá)垂尾時(shí)，靠近垂尾翼梢位置的小渦結(jié)構(gòu)被衰減掉許多。由圖9基于等Q值圖顯示的空間流場渦結(jié)構(gòu)可以分析出，機(jī)身迎角為30° 時(shí)，垂尾可能產(chǎn)生較大的抖振響應(yīng)。

圖9　等Q值圖Fig.9　Contour of Q value

圖10給出了左垂尾外翼面(靠近機(jī)翼翼梢方向)脈動壓力系數(shù)CPrms分布云圖，脈動壓力系數(shù)通過脈動壓力均方根值定義得到。從圖10中可以看出，垂尾前緣的壓力脈動要高于后緣位置的壓力脈動值，說明沿著氣流運(yùn)動的方向(X方向)，湍流脈動是在逐步衰減的。當(dāng)模型迎角α為20° 時(shí)，垂尾翼根前緣附近的壓力脈動較大，而翼梢、翼根后緣位置附近的壓力脈動接近零，如圖10(a)，說明此時(shí)邊條翼形成的脫體渦在垂尾前緣還保持大渦狀態(tài)，脫體渦離機(jī)翼上翼面距離較近，脫體渦在垂尾1/2弦線位置附近發(fā)生破裂后僅對翼根的脈動壓力產(chǎn)生了影響；當(dāng)模型迎角α為30° 時(shí)，垂尾翼根和翼梢前緣位置均出現(xiàn)較大的壓力脈動，而翼梢后緣位置的壓力脈動值也較大，如圖10(b)，說明邊條翼形成的脫體渦在機(jī)翼前緣附近發(fā)生了破裂，整個(gè)垂尾沉浸在脫體渦破裂后的寬頻湍流脈動氣流中；當(dāng)模型迎角α為40° 時(shí)，垂尾翼根前緣附近壓力脈動高于其他位置，而翼梢、翼根后緣位置的壓力脈動值均不是很大，如圖10(c)，說明邊條翼形成的脫體渦在垂尾較前位置已經(jīng)發(fā)生了破裂，破裂后的脫體渦對垂尾翼根前緣附近的壓力脈動產(chǎn)生了影響，而垂尾大部分位置，尤其是垂尾翼梢處于充分發(fā)展的寬頻湍流氣流中，脈動壓力值相對機(jī)身迎角為30° 時(shí)有所降低。

圖10　垂尾脈動壓力系數(shù)分布云圖Fig.10　Contour of fluctuating pressure coefficient distribution on vertical tail

為了更好地比較不同機(jī)身迎角下的垂尾壓力脈動值差異，在垂尾外翼表面翼根、翼中和翼梢位置選擇了3條監(jiān)測線，每條監(jiān)測線上均勻分布9個(gè)點(diǎn)，如圖11所示。

圖11　垂尾脈動壓力監(jiān)測點(diǎn)位置Fig.11　Monitoring point location of fluctuating pressure of vertical tail

圖12給出了垂尾不同監(jiān)測點(diǎn)的脈動壓力系數(shù)分布?？梢钥闯?，不同機(jī)身迎角下翼根位置(Line1)和翼中位置(Line2)的脈動壓力系數(shù)分布差異并不明顯，而對于翼梢位置(Line3)，機(jī)身迎角變化對脈動壓力系數(shù)影響很大，機(jī)身迎角30° 時(shí)翼梢位置(Line3)的壓力脈動要明顯大于迎角20° 和40° 時(shí)的壓力脈動。翼梢位置的壓力脈動耦合力臂效應(yīng)，會在垂尾翼根部位形成較大的彎曲力矩，使得垂尾發(fā)生疲勞破壞，因此，需要避免翼梢位置產(chǎn)生較大的脈動壓力。圖12的結(jié)果再次驗(yàn)證了前面關(guān)于不同機(jī)身迎角下垂尾脈動壓力系數(shù)分布的論述。

圖12　監(jiān)測點(diǎn)脈動壓力系數(shù)分布Fig.12　Fluctuating pressure coefficient distribution of monitoring points

圖13給出了不同機(jī)身迎角下監(jiān)測點(diǎn)1位置脈動壓力的功率譜密度曲線，監(jiān)測點(diǎn)1的位置如圖11所示，為垂尾翼根靠近前緣的一個(gè)點(diǎn)。從圖13 可以看出，機(jī)身迎角為20° 時(shí)，監(jiān)測點(diǎn)1位置脈動壓力功率譜密度曲線存在一個(gè)峰值，對應(yīng)圖9關(guān)于流場渦結(jié)構(gòu)的分析，該脈動峰值應(yīng)是由邊條翼形成的脫體渦引起的，脈動峰值對應(yīng)的頻率為邊條翼脫體渦脫落的頻率，約為380Hz；而當(dāng)機(jī)身迎角為30°、40° 時(shí)，脫體渦在垂尾前方已經(jīng)發(fā)生了破裂，因此，監(jiān)測點(diǎn)1位置對應(yīng)的脈動壓力功率譜密度曲線表現(xiàn)為寬頻模式，沒有明顯的脈動峰值，如圖13所示。

圖13　監(jiān)測點(diǎn)1脈動壓力的功率譜密度曲線Fig.13　Power spectrum density curves of fluctuating pressure of monitoring point 1

圖14給出了非定常RANS/LES混合方法計(jì)算的垂尾側(cè)向力系數(shù)CC的時(shí)間響應(yīng)曲線。

圖14　垂尾側(cè)向力系數(shù)響應(yīng)曲線Fig.14　Side force coefficient response curves of vertical tail

可以看出，3個(gè)機(jī)身迎角下，RANS/LES混合方法能夠比較好地模擬出垂尾側(cè)向力的響應(yīng)特性。從模擬結(jié)果可以看出，隨著機(jī)身迎角增加，垂尾的側(cè)向力的絕對值在減小，但側(cè)向力的脈動水平在增長。

為了更詳細(xì)地比較側(cè)向力的脈動差異，圖15給出了不同機(jī)身迎角下的垂尾側(cè)向力系數(shù)功率譜密度曲線?？梢钥闯?，在低頻區(qū)域，隨著機(jī)身迎角的增加，垂尾側(cè)向力脈動在增大；在高頻區(qū)域，機(jī)身迎角30° 的垂尾側(cè)向力脈動由于小渦結(jié)構(gòu)在垂尾部分的耗散而與機(jī)身迎角20° 時(shí)的垂尾側(cè)向力脈動相當(dāng)，而機(jī)身迎角為40° 時(shí)，小渦結(jié)構(gòu)在到達(dá)垂尾區(qū)域已經(jīng)被耗散，因此，此時(shí)的垂尾側(cè)向力脈動要低于機(jī)身迎角為20°、30° 時(shí)的垂尾側(cè)向力脈動。

圖15　垂尾側(cè)向力系數(shù)功率譜密度曲線Fig.15　Power spectrum density curves of side force coefficient of vertical tail

由圖15可以看出，機(jī)身迎角30° 和40° 時(shí)，垂尾側(cè)向力的低頻脈動量相當(dāng)，這是由于側(cè)向力是整個(gè)垂尾的積分力，從圖10垂尾脈動壓力系數(shù)分布可以看出，機(jī)身迎角40° 時(shí)翼根位置的脈動壓力要高于30° 時(shí)的脈動壓力，翼梢位置則相反。雖然機(jī)身迎角30° 和40° 時(shí)垂尾側(cè)向力的低頻脈動量相當(dāng)，但翼梢位置的壓力脈動會引起更大的垂尾抖振響應(yīng)。

4　結(jié)　論

在亞跨超CFD軟件平臺(TRIP)的框架內(nèi)，發(fā)展了基于RANS/LES混合方法的流動模擬技術(shù)，通過NACA0021翼型大迎角分離流動和串列雙圓柱繞流算例對RANS/LES混合方法的模擬精度進(jìn)行了驗(yàn)證，并在此基礎(chǔ)研究了某戰(zhàn)斗機(jī)模型在大迎角狀態(tài)下的垂尾脈動壓力。從數(shù)值模擬結(jié)果可以得到以下結(jié)論：

不同迎角下垂尾脈動壓力差異主要體現(xiàn)在翼梢位置，當(dāng)機(jī)身迎角達(dá)到一定值時(shí)，垂尾完全沉浸在脫體渦破裂后的寬頻湍流脈動中，翼梢位置的脈動壓力較明顯，當(dāng)機(jī)身迎角較小時(shí)，脫體渦靠近主機(jī)翼，翼梢位置沒有破碎的脫體渦，而當(dāng)機(jī)身迎角較大時(shí)，翼梢位置的破碎脫體渦被充分衰減。

[1]賈有， 楊智春. 一種飛機(jī)垂尾抖振載荷識別的新方法[J]. 航空學(xué)報(bào)， 2013, 34(10): 2333-2340.

JIA Y, YANG Z C. A new approach to identify buffet loads for aircraft vertical tail[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(10): 2333-2340 (in Chinese).

[2]韓冰， 徐敏， 蔡天星， 等. 渦破裂誘導(dǎo)的垂尾抖振數(shù)值模擬[J]. 航空學(xué)報(bào)， 2012, 33(5): 788-795.

HAN B, XU M, CAI T X, et al. Numerical simulation of vertical tail buffeting induced by vortex breakdown[J]. Acta Aeronautica et Astronautica Sinica, 2012, 32(5): 788-795 (in Chinese).

[3]HUANG J T, GAO Z H. Buffeting numerical simulation coupled with aerodynamics and structure based on MDDES[J]. Science China Technological Sciences, 2013, 56: 1550-1560.

[4]AHMED E, OSAMA A K, OKTAY B. F/A-18 twin-tail buffet modeling using non-linear eddy viscosity models： AIAA-2014-2447[R]. Reston: AIAA, 2014.

[5]李勁杰， 楊青， 楊永年， 等. 邊條翼布局雙垂尾抖振表面壓力脈動風(fēng)洞實(shí)驗(yàn)研究[J]. 實(shí)驗(yàn)流體力學(xué)， 2006, 20(3): 29-32.

LI J J, YANG Q, YANG Y N, et al. Wind-tunnel unsteady pressure measurements of twin-vertical tail during buffet of strake-wing configuration[J]. Journal of Experiments in Fluid Mechanics, 2006, 20(3): 29-32 (in Chinese).

[6]SPALART P R, JOU W H, STRETLETS M, et al. Comments on the feasibility of LES for wings and on the hybrid RANS/LES approach[C]//Proceedings of the First AFOSR International Conference on DNS/LES. Colombus: Greyden Press, 1998: 137-147.

[7]DECK S. Zonal-detached-eddy simulation of the flow around a high-lift configuration[J]. AIAA Journal, 2005, 43(11): 2372-2384.

[8]SPALART P R, DECK S, SHUR M, et al. A new version of detached-eddy simulation, resistant to ambiguous grid densities[J]. Theory and Computation Fluid Dynamics, 2006, 20(3): 181-195.

[9]SPALART P R. Detached-eddy simulation[J]. Annual Review of Fluid Mechanics, 2009, 41(1): 181-202.

[10]SHUR MIKHAIL L, SPALART P R, STRELETS M K, et al. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities[J]. International Journal of Heat and Fluid Flow, 2008, 29(6): 1638-1649.

[11]MIKHAIL S G, ANDREY V G, JOCHEN S, et al. Development of DDES and IDDES formulations for thek-ωshear stress tranport model[J]. Flow Turbulence Combust, 2012, 88(3): 431-449.

[12]DECK S. Numerical simulation of transonic buffet over a supercritical airfoil[J]. AIAA Journal, 2005, 43(7): 1556-1566.

[13]吳晶峰， 寧方飛. 后臺階流動的Hybrid RANS/LES模擬[J]. 北京航空航天大學(xué)學(xué)報(bào)， 2011, 37(6): 701-704.

WU J F, NING F F. Hybrid RANS-LES method applied to backward facing step flow[J]. Journal of Beijing University of Aeronautics and Astronautics, 2011, 37(6): 701-704 (in Chinese).

[14]XIAO Z X, LIU J, LUO K Y, et al. Numerical investigations of massively separated flows past rudimentary landing gear using SST-DDES: AIAA-2012-0385[R]. Reston: AIAA, 2012.

[15]HUANG J B, XIAO Z X, LIU J, et al. Simulation of shock wave buffet and its suppression on an OAT15A supercritical airfoil by IDDES[J]. Science China: Physics, Mechanics and Astronomy, 2012, 55(2): 260-271.

[16]王運(yùn)濤， 王光學(xué)， 張玉倫. TRIP2.0軟件的確認(rèn)： DPWII復(fù)雜組合體的數(shù)值模擬[J]. 航空學(xué)報(bào)， 2008, 29(1): 34-40.

WANG Y T, WANG G X, ZHANG Y L. Validation of TRIP2.0: Numerical simulation of DPWII complex configuration[J]. Acta Aeronautica et Astronautica Sinica, 2008, 29(1): 34-40 (in Chinese).

[17]王運(yùn)濤， 王光學(xué)， 張玉倫. 采用TRIP2.0軟件計(jì)算DLR-F6構(gòu)型的阻力[J]. 空氣動力學(xué)學(xué)報(bào)， 2009, 27(1): 108-113.

WANG Y T, WANG G X, ZHANG Y L. Drag prediction of DLR-F6 configuration with TRIP2.0 software[J]. Acta Aerodynamica Sinica, 2009, 27(1): 108-113 (in Chinese).

[18]王運(yùn)濤， 王光學(xué)， 張玉倫. DPWIII機(jī)翼和翼身組合體構(gòu)型數(shù)值模擬[J]. 空氣動力學(xué)學(xué)報(bào)， 2011, 29(3): 264-269.

WANG Y T, WANG G X, ZHANG Y L. Numerical simulation of DPW III wing and wing-body configurations[J]. Acta Aerodynamica Sinica, 2011, 29(3): 264-269 (in Chinese).

[19]王運(yùn)濤， 張書俊， 孟德虹. DPW4翼/身/平尾組合體的數(shù)值模擬[J]. 空氣動力學(xué)學(xué)報(bào)， 2013, 31(6): 739-744.

WANG Y T, ZHANG S J, MENG D H. Numerical simulation and study for DPW4 wing/body/tail[J]. Acta Aerodynamica Sinica, 2013, 31(6): 739-744 (in Chinese).

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

[21]MENTER F R. Two equation eddy viscosity turbulence models for engineering application[J]. AIAA Journal, 1994, 32(8): 1598-1605.

孟德虹男， 碩士， 助理研究員。主要研究方向： 計(jì)算空氣動力學(xué)。

Tel: 0816-2463274

E-mail: mdh157@163.com

孫巖男， 博士研究生， 助理研究員。主要研究方向： 計(jì)算氣動彈性力學(xué)。

Tel: 0816-2463205

E-mail: supersunyan@163.com

王運(yùn)濤男， 博士， 研究員， 博士生導(dǎo)師。主要研究方向： 計(jì)算空氣動力學(xué)。

Tel: 0816-2463015

E-mail: ytwang@skla.cardc.cn

李偉男， 博士研究生， 助理研究員。主要研究方向： 計(jì)算空氣動力學(xué)。

Tel：0816-2463274

E-mail: kuaile06@163.com

Numerical simulation of fluctuating pressure of fighter vertical tail

MENG Dehong1, SUN Yan2, WANG Yuntao1,*, LI Wei1

1. Computational Aerodynamics Institute, China Aerodynamics Research and Development Center,Mianyang621000, China 2. State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center,Mianyang621000, China

IDDES flow simulation technique based on hybrid RANS/LES method is realized on TRIP, a trisonic fluid simulation software platform. Separating flows around airfoil NACA0021 with 60° angle of attack and tandem cylinders are simulated to validate the precision of hybrid RANS/LES method. Then, fluctuating pressure of a fighter vertical tail is calculated with the presented RANS/LES method. Flow Mach number is 0.1, Reynolds number based on length of fighter is 2 million, and angles of attack are 20°, 30° and 40°. Fluctuating pressure of fighter vertical tail is analyzed through fluctuating pressure coefficients, power spectrum density of fluctuating pressure, space flow structure and side force response. Numerical results show that fluctuating pressure on tip of vertical tail increases obviously when vertical tail is immersed in broadband turbulence fluctuating flow behind the breakup of separating vortex from strake wing.

fighter; vertical tail; fluctuating pressure; RANS/LES； computational fluid dynamics (CFD)

2016-01-18； Revised： 2016-01-23； Accepted： 2016-04-11； Published online： 2016-04-2616:02

National Key Research and Development Plan (2016YFB0200700)

. Tel.： 0816-2463015E-mail: ytwang@skla.cardc.cn

2016-01-18； 退修日期： 2016-01-23； 錄用日期： 2016-04-11；

時(shí)間： 2016-04-2616:02

www.cnki.net/kcms/detail/11.1929.V.20160426.1602.002.html

國家重點(diǎn)研發(fā)計(jì)劃 (2016YFB0200700)

.Tel.： 0816-2463015E-mail: ytwang@skla.cardc.cn

10.7527/S1000-6893.2016.0121

V211.3

A

1000-6893(2016)08-2472-09

引用格式： 孟德虹， 孫巖， 王運(yùn)濤， 等． 戰(zhàn)斗機(jī)垂尾脈動壓力數(shù)值模擬[J]． 航空學(xué)報(bào), 2016, 37(8): 2472-2480. MENG D H, SUN Y, WANG Y T, et al. Numerical simulation of fluctuating pressure of fighter vertical tail[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2472-2480.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

URL： www.cnki.net/kcms/detail/11.1929.V.20160426.1602.002.html