基于動(dòng)網(wǎng)格技術(shù)的機(jī)翼動(dòng)態(tài)失速仿真分析

錢 宇， 蔣 皓

(中國(guó)民用航空飛行學(xué)院飛行技術(shù)學(xué)院， 廣漢 618307)

失速時(shí)機(jī)翼表面氣流發(fā)生嚴(yán)重分離，導(dǎo)致升力驟降、阻力驟升，進(jìn)而造成航空器的失控，嚴(yán)重影響運(yùn)行安全。有別于航空器設(shè)計(jì)時(shí)考慮的靜失速，動(dòng)態(tài)失速是航空器在實(shí)際環(huán)境下的運(yùn)行狀態(tài)，關(guān)注其演化機(jī)理有助于預(yù)防和改出失速。因此，對(duì)于動(dòng)態(tài)失速的研究不容忽視。

動(dòng)態(tài)失速計(jì)算通常采用雷諾時(shí)均，直接數(shù)值模擬需要大量的計(jì)算資源，大渦模擬雖然在一定程度上減小了計(jì)算量，但對(duì)于工程應(yīng)用來說其計(jì)算時(shí)間和計(jì)算量依舊是不可以被接受的。在計(jì)算精度可接受的情況下，雷諾時(shí)均法具有計(jì)算時(shí)間較短的優(yōu)勢(shì)。張彥軍等[1]將雷諾時(shí)均與剪切應(yīng)力傳輸(shear stress transfer, SST)k-ω湍流模型結(jié)合，計(jì)算了不同雷諾數(shù)工況下邊界層變化情況及動(dòng)態(tài)失速相關(guān)特性。Li等[2]用動(dòng)網(wǎng)格技術(shù)計(jì)算了對(duì)翼型在不同振蕩運(yùn)動(dòng)狀態(tài)下的動(dòng)態(tài)失速特性。Karbasian 等[3]用SSTk-ω湍流模型對(duì)S809機(jī)翼在不同頻率下的動(dòng)態(tài)失速進(jìn)行了數(shù)值計(jì)算，得出在穩(wěn)態(tài)和非穩(wěn)態(tài)情況下，所得計(jì)算結(jié)果與實(shí)驗(yàn)均相似。張衛(wèi)國(guó)等[4]采用等離子體流動(dòng)控制對(duì)風(fēng)力機(jī)翼的翼型的動(dòng)態(tài)失速特性進(jìn)行了數(shù)值研究。Surekha等[5]對(duì)振蕩的NACA0012翼型的深動(dòng)態(tài)失速現(xiàn)象進(jìn)行了數(shù)值模擬，并與實(shí)驗(yàn)對(duì)比以研究深動(dòng)態(tài)失速狀態(tài)的空氣動(dòng)力學(xué)特性。

研究針對(duì)機(jī)翼動(dòng)態(tài)失速問題，用O-Block法生成結(jié)構(gòu)網(wǎng)格，用于構(gòu)建旋轉(zhuǎn)區(qū)域。區(qū)域的整體旋轉(zhuǎn)有助于網(wǎng)格更新。通過翼型動(dòng)態(tài)失速計(jì)算進(jìn)行了驗(yàn)證。用Thin-cut 技術(shù)生成機(jī)翼的非結(jié)構(gòu)網(wǎng)格，以改善尖錐薄壁面網(wǎng)格質(zhì)量。將SSTk-ω湍流模型與Navier-Stokes(N-S)方程相結(jié)合對(duì)翼型和機(jī)翼的動(dòng)態(tài)失速過程進(jìn)行了數(shù)值求解，進(jìn)而分析機(jī)翼動(dòng)態(tài)失速演化特性。

1 飛機(jī)俯仰運(yùn)動(dòng)方程

根據(jù)動(dòng)態(tài)失速情況下機(jī)翼的運(yùn)動(dòng)特性，以模型重心為擺動(dòng)中心，俯仰擺動(dòng)以模擬動(dòng)態(tài)失速過程中迎角的變化情況，如圖1所示。

圖1 機(jī)翼運(yùn)動(dòng)過程示意圖Fig.1 Schematic diagram of wing motion process

俯仰運(yùn)動(dòng)方程[6]為

α(t)=α0+αlsin(ωt)

(1)

(2)

式中：α為運(yùn)動(dòng)過程的平均迎角；α0為運(yùn)動(dòng)起始的迎角；ω為運(yùn)動(dòng)過程的角速度；αl為運(yùn)動(dòng)過程的擺動(dòng)振幅；ω2為模型繞z軸旋轉(zhuǎn)的角速度。

式(1)表示迎角隨時(shí)間的變化情況，對(duì)迎角進(jìn)行一次求導(dǎo)可得控制運(yùn)動(dòng)方程[式(2)]。通過用戶自定義函數(shù)對(duì)該方程進(jìn)行編譯，從而可達(dá)到控制網(wǎng)格動(dòng)態(tài)變化的目的。

2 網(wǎng)格劃分

考慮到網(wǎng)格運(yùn)動(dòng)時(shí)形變的問題，若僅使模型轉(zhuǎn)動(dòng)，將帶來巨大的網(wǎng)格形變量，尤其是計(jì)算精度要求很高，需要大量計(jì)算資源與計(jì)算時(shí)間。因此，在保證計(jì)算精度的前提下，可使靠近模型的區(qū)域整體運(yùn)動(dòng)，這樣可以大大減小網(wǎng)格的形變量，減小計(jì)算資源節(jié)約計(jì)算時(shí)間。

計(jì)算流場(chǎng)網(wǎng)格如圖2所示。采用結(jié)構(gòu)網(wǎng)格與非結(jié)構(gòu)網(wǎng)格構(gòu)成運(yùn)動(dòng)內(nèi)域與固定外域。用O-Block技術(shù)生成旋轉(zhuǎn)內(nèi)域結(jié)構(gòu)網(wǎng)格，該區(qū)域隨翼型共同運(yùn)動(dòng)。固定外域采用非結(jié)構(gòu)網(wǎng)格，該區(qū)域的網(wǎng)格采用三角網(wǎng)格劃分，再結(jié)合網(wǎng)格光順與網(wǎng)格重構(gòu)技術(shù)，能很好地處理網(wǎng)格形變問題，將兩個(gè)網(wǎng)格用interface面連接以組成完整的計(jì)算網(wǎng)格。網(wǎng)格質(zhì)量如圖3所示。

圖2 翼型計(jì)算流場(chǎng)網(wǎng)格劃分Fig.2 Calculated flow field meshing of airfoil

圖3 網(wǎng)格質(zhì)量檢查結(jié)果Fig.3 Result of grid quality checking

網(wǎng)格單元質(zhì)量數(shù)值越趨近于1則網(wǎng)格質(zhì)量越高。數(shù)值計(jì)算要求網(wǎng)格質(zhì)量盡量大于0.3且不能存在小于0的負(fù)體積網(wǎng)格。

雖然結(jié)構(gòu)網(wǎng)格的網(wǎng)格質(zhì)量好于非結(jié)構(gòu)網(wǎng)格，但是在模型的網(wǎng)格劃分時(shí)，由于模型存在較多的不規(guī)則面與尖錐面，生成結(jié)構(gòu)網(wǎng)格需要耗費(fèi)大量時(shí)間，非結(jié)構(gòu)網(wǎng)格卻能較好能適應(yīng)這種情況。因此，機(jī)翼網(wǎng)格的生成采用了非結(jié)構(gòu)網(wǎng)格，其計(jì)算網(wǎng)格如圖4所示。旋轉(zhuǎn)內(nèi)域和運(yùn)動(dòng)外域流場(chǎng)網(wǎng)格均采用非結(jié)構(gòu)網(wǎng)格生成，采用interface面進(jìn)行能量傳輸。

圖4 機(jī)翼計(jì)算流場(chǎng)網(wǎng)格Fig.4 Calculated flow field of wing

Thin-cut技術(shù)能使得薄壁尖錐面處的網(wǎng)格附面情況滿足要求，不會(huì)出現(xiàn)鋸齒狀網(wǎng)格丟失。如圖5所示，機(jī)翼后緣是尖錐面，單純縮小網(wǎng)格尺寸雖然可以改善機(jī)翼后緣網(wǎng)格的分布，但是依然無法使網(wǎng)格完全附面，并且大大增加了不必要的網(wǎng)格數(shù)量，如圖5(a)所示。采用Thin-cut技術(shù)后，可以在不增加網(wǎng)格數(shù)量的情況下，使非結(jié)構(gòu)網(wǎng)格在尖錐面生成較為完善的網(wǎng)格，如圖5(b)所示。

圖5 Thin-cut技術(shù)對(duì)非結(jié)構(gòu)網(wǎng)格的影響Fig.5 Influence of thin cut technology on unstructured grid

對(duì)機(jī)翼網(wǎng)格進(jìn)行質(zhì)量檢測(cè)結(jié)果大于0.3，滿足數(shù)值求解要求。

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

3.1 數(shù)值方法與湍流模型

對(duì)于流體力學(xué)的計(jì)算，采用雷諾時(shí)均，該方法在計(jì)算精度上滿足要求，在計(jì)算時(shí)間和計(jì)算資源上較為合理。雷諾時(shí)均[7]利用N-S方程與湍流模型相結(jié)合，使得方程封閉可用于數(shù)值計(jì)算。N-S方程的矢量形式[8]為

(3)

式(3)中：u為速度；p為壓力；ρ為密度；υ為流體黏度系數(shù)；f為外力；左邊第一項(xiàng)為時(shí)變慣性力，第二項(xiàng)為位變慣性力，右邊第一項(xiàng)為質(zhì)量力，第二項(xiàng)為壓差力，第三項(xiàng)為黏性力。

在湍流模型[9-11]中，一階S-A模型不予考慮，該模型在計(jì)算過程中耗散率大，在處理高轉(zhuǎn)捩區(qū)域時(shí)并不適用，在層流的計(jì)算時(shí)可適當(dāng)運(yùn)用。二階模型有Standardk-ε、Standardk-ω、SSTk-ω等模型，不同二階模型適用范圍不同。

根據(jù)文獻(xiàn)[12]所做仿真計(jì)算與實(shí)驗(yàn)結(jié)果對(duì)比可知，在計(jì)算失速時(shí)，SSTk-ω湍流模型對(duì)于具有高逆壓梯度差轉(zhuǎn)捩區(qū)域的計(jì)算較為友好，其計(jì)算結(jié)果最接近實(shí)驗(yàn)結(jié)果。

3.2 可行性驗(yàn)證

為進(jìn)一步證明SSTk-ω湍流模型的可行性，計(jì)算雷諾數(shù)為2.2×106，來流速度為0.13Ma(Ma為馬赫數(shù))下NACA2412翼型升力系數(shù)穩(wěn)態(tài)變化，并與Seetharam等[13]所做實(shí)驗(yàn)結(jié)果對(duì)比，結(jié)果如圖6所示。

由圖6可知，SSTk-ω湍流模型計(jì)算所得結(jié)果與實(shí)驗(yàn)結(jié)果的相對(duì)誤差為8.63%，在計(jì)算接受的范圍之內(nèi)，因此該模型可用于計(jì)算具有高逆壓梯度差的失速模型。

4 計(jì)算結(jié)果與分析

4.1 翼型動(dòng)態(tài)失速結(jié)果分析

圖7所示為翼型升力系數(shù)隨擺動(dòng)迎角的變化。由圖7可知，翼型擺動(dòng)分為2個(gè)階段：迎角自0°下擺至-25°再上仰至0°；迎角自0°迎角上仰至25°再下擺至0°。

研究著重關(guān)注第2階段。在該階段，當(dāng)翼型上仰擺動(dòng)到18.9°時(shí)升力系數(shù)達(dá)到峰值，其數(shù)值為2.57。當(dāng)迎角達(dá)到動(dòng)態(tài)失速臨界迎角后，渦旋在機(jī)翼上表面的反復(fù)脫離-再生成-附著-再脫離的過程，使得升力系數(shù)呈上下波動(dòng)變化，直至生成的渦旋完全不再附著于機(jī)翼上表面。

渦量能定性定量地分析動(dòng)態(tài)失速情況下渦旋的變化情況。圖8為翼型上仰工況下渦量隨迎角的變化。當(dāng)迎角不斷增加并達(dá)到某一定值，前緣開始生成負(fù)向渦旋并向后緣移動(dòng)，當(dāng)負(fù)向渦旋移動(dòng)到翼型后緣時(shí)，其渦旋力將氣流從翼型下表面吸引至吸引上表面，從而卷起正向渦旋。正向渦旋同樣會(huì)提供額外升力，與負(fù)向渦旋共同為翼型提供額外升力，直至迎角度數(shù)達(dá)到17.8°，此時(shí)機(jī)翼上表面渦旋徹底脫離，再次生成的渦旋也幾乎不附面，僅對(duì)升力有一定的擾動(dòng)。

4.2 機(jī)翼動(dòng)態(tài)失速結(jié)果分析

圖9所示為翼尖截面處渦量分布。由圖9可知，無論機(jī)翼處于任何迎角，其上表面都有渦旋附著。翼尖截面處相對(duì)于機(jī)翼的其他截面較為特殊，上下表面的壓力差對(duì)翼梢處的氣流產(chǎn)生自下而上的吸引，使氣流從下表面向上表面移動(dòng)，從而產(chǎn)生渦旋，即翼尖渦。

圖9 翼尖截面處渦量分布Fig.9 Variation of vorticity at wingtip section

圖10所示為翼根截面處不同迎角下渦量的分布。由圖10可知，當(dāng)迎角小于12.9°時(shí)，在該截面附近的機(jī)翼上表面附著負(fù)向渦旋，隨著迎角增加，負(fù)向渦旋向后緣移動(dòng)，直到13.5°迎角時(shí)負(fù)向渦旋脫落并吸引下表面氣流在上表面產(chǎn)生正向渦旋。隨著迎角度數(shù)的繼續(xù)增大，雖然有正負(fù)向渦旋生成，但均已不依附于此處上表面。

圖10 翼根截面處渦量變化情況Fig.10 Variation of vorticity at wing root section

圖11為距翼根1/4機(jī)翼長(zhǎng)度截面處渦量變化。可知，當(dāng)迎角達(dá)到11.5°時(shí)，該截面處附近的負(fù)向渦旋已經(jīng)開始有脫離的趨勢(shì)，且正向渦旋逐步形成；當(dāng)達(dá)到12.9°時(shí)，截面處附近機(jī)翼表面已完全被正向渦旋覆蓋，并伴著正向渦旋的脫離；之后，隨著迎角度數(shù)的增加，僅有負(fù)向渦旋于機(jī)翼前緣生成且直接脫離，上表面無任何渦旋附著。

圖12為距翼根3/4機(jī)翼長(zhǎng)度截面處渦量變化。

圖12 距翼根3/4機(jī)翼長(zhǎng)度截面處渦量變化情況Fig.12 Variation of vorticity at the section of 3/4 wing length from the root

可知，當(dāng)迎角度數(shù)達(dá)到12.9°時(shí)，該截面處的負(fù)向渦旋開始脫離；之后，隨著迎角繼續(xù)增加，從前緣生成的負(fù)向渦旋脫離機(jī)翼。

由上述分析可知，機(jī)翼在發(fā)生動(dòng)態(tài)失速時(shí)，機(jī)翼不同截面處渦旋脫離機(jī)翼表面的迎角并不相同，即不同截面處有不同的臨界迎角；在負(fù)向渦旋從前緣移動(dòng)到后緣時(shí)，并不是所有截面處都會(huì)吸引下表面氣流向上表面運(yùn)動(dòng)從而生成正向渦旋；翼尖截面處的渦旋同時(shí)受到翼尖渦的影響，其動(dòng)態(tài)失速特性有別于其他截面處。

若要得到更為精確的機(jī)翼動(dòng)態(tài)失速情況，則需在計(jì)算過程中加入更多的擾動(dòng)以模擬航空器在真實(shí)的環(huán)境中的運(yùn)行狀態(tài)，這將成為后續(xù)工作的研究重點(diǎn)。

5 結(jié)論

(1)機(jī)翼處于動(dòng)態(tài)失速狀態(tài)下，不同截面處渦量的變化情況不同。隨著迎角增加，負(fù)向的渦旋從機(jī)翼前緣生成，并向后緣移動(dòng)。當(dāng)負(fù)向渦旋移動(dòng)到后緣時(shí)，其吸引力不僅使機(jī)翼獲得了額外的升力，同時(shí)吸引了下表面氣流向上表面流動(dòng)，從而在上表面生成正向渦旋，這使得動(dòng)態(tài)失速過程中機(jī)翼上表面呈現(xiàn)正負(fù)向渦旋交替出現(xiàn)的現(xiàn)象。

(2)渦旋首先在機(jī)翼中部生成，并向翼梢和翼根發(fā)展，使得機(jī)翼中部首先發(fā)生動(dòng)態(tài)失速，因此機(jī)翼中部動(dòng)態(tài)失速臨界迎角小于翼尖及翼根處臨界迎角。較翼型而言，機(jī)翼具有更為復(fù)雜的氣動(dòng)力特性，同時(shí)受到機(jī)身的影響，導(dǎo)致其各個(gè)截面的動(dòng)態(tài)失速臨界迎角均小于翼型失速迎角。