考慮二次本構(gòu)關(guān)系的湍流模型對(duì)翼身組合體阻力預(yù)測(cè)的影響分析

趙 輝, 張耀冰, 陳江濤, 馬明生

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 綿陽 621000)

0 引 言

飛行器氣動(dòng)力和力矩的精細(xì)預(yù)測(cè)是其設(shè)計(jì)過程中的關(guān)鍵技術(shù)之一。阻力量化和減阻設(shè)計(jì)對(duì)燃料和成本的節(jié)省十分關(guān)鍵。Gilyard等[1]指出，對(duì)于一架現(xiàn)代寬體遠(yuǎn)程運(yùn)輸機(jī)，在最大起飛重量的情況下，減阻1%每年可以帶來$4 000 000的收益，包括燃油的節(jié)省和載貨量的增加。不同外形之間阻力增量的準(zhǔn)確預(yù)測(cè)是減阻嘗試的基礎(chǔ)。近些年來計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)在網(wǎng)格生成技術(shù)、流場(chǎng)求解和高性能計(jì)算機(jī)等方面取得了很大的進(jìn)步，這使得通過數(shù)值模擬來評(píng)估三維真實(shí)運(yùn)輸機(jī)外形的氣動(dòng)性能成為可能，CFD預(yù)測(cè)工業(yè)相關(guān)外形氣動(dòng)特性的能力得到了廣泛的證明[2-5]。然而，對(duì)于一些非線性效應(yīng)主導(dǎo)的流動(dòng)問題，比如二次流、激波誘導(dǎo)分離等，CFD的預(yù)測(cè)精度仍然不足[6]。

與此同時(shí)，CFD的驗(yàn)證和確認(rèn)在CFD業(yè)界受到廣泛的關(guān)注。為了評(píng)估現(xiàn)有的CFD技術(shù)預(yù)測(cè)工業(yè)相關(guān)外形氣動(dòng)特性的能力，AIAA應(yīng)用空氣動(dòng)力學(xué)委員會(huì)在2001年發(fā)起了AIAA阻力預(yù)測(cè)會(huì)議(Drag Prediction Workshop, DPW)，到2016年6月已經(jīng)成功舉辦了六屆。會(huì)議期間，參會(huì)者深入討論了氣動(dòng)力和力矩預(yù)測(cè)的可靠性。有關(guān)會(huì)議的細(xì)節(jié)可以在參考文獻(xiàn)[7-11]中查閱。

正如DPW-V的總結(jié)報(bào)告里展示的一樣，當(dāng)流動(dòng)是附著流或者分離可以忽略的時(shí)候，大部分的數(shù)值模擬都可以得到與試驗(yàn)比較一致的結(jié)果，滿足工業(yè)要求。但是CFD預(yù)測(cè)精度在大迎角區(qū)域仍顯不足，比如在翼身組合體外形機(jī)翼和機(jī)身結(jié)合處分離流動(dòng)的預(yù)測(cè)上，數(shù)值模擬的結(jié)果跟網(wǎng)格拓?fù)浜褪杳芊植肌⑼牧髂Ｐ汀ば皂?xiàng)處理等密切相關(guān)，分散度很大。很多學(xué)者研究發(fā)現(xiàn)，基于Bousinessq線性渦黏假設(shè)的湍流模型會(huì)高估分離區(qū)的大小，導(dǎo)致預(yù)測(cè)的氣動(dòng)力和俯仰力矩與試驗(yàn)值偏差較遠(yuǎn)。究其原因是線性渦黏模型不能準(zhǔn)確預(yù)測(cè)拐角流動(dòng)中雷諾應(yīng)力的各向異性引起的二次流動(dòng)。在湍流模型中考慮二次本構(gòu)關(guān)系(Quadratic Constitutive Relation, QCR)后，分離區(qū)大小和氣動(dòng)特性的預(yù)測(cè)與試驗(yàn)結(jié)果更加吻合。不過以往的研究主要集中在QCR對(duì)氣動(dòng)特性預(yù)測(cè)影響的外在表現(xiàn)上，沒有深入分析產(chǎn)生影響的內(nèi)在機(jī)理，而這也是本文的研究重點(diǎn)。

本文對(duì)DPW-VI外形進(jìn)行了數(shù)值模擬，著重研究考慮二次本構(gòu)關(guān)系的湍流模型對(duì)模擬結(jié)果的影響，并深入分析拐角附近的局部流動(dòng)。重點(diǎn)分析了考慮二次本構(gòu)關(guān)系的湍流模型提高模擬精度的原因和機(jī)理，合理地利用該湍流模型能夠得到更加合理的壓力分布和流動(dòng)分離情況，對(duì)氣動(dòng)特性的預(yù)測(cè)尤其是阻力預(yù)測(cè)具有一定的指導(dǎo)意義。

1 計(jì)算外形和網(wǎng)格

DPW-VI的外形如圖1所示，是NASA和阻力會(huì)議組委會(huì)聯(lián)合設(shè)計(jì)發(fā)展的CRM模型。該模型代表了典型現(xiàn)代跨聲速運(yùn)輸機(jī)外形，設(shè)計(jì)巡航狀態(tài)是Ma=0.85，CL=0.5。在NASA NTF風(fēng)洞、NASA Ames 11-ft風(fēng)洞和歐洲ETW風(fēng)洞進(jìn)行了測(cè)試，有大量數(shù)據(jù)可供CFD確認(rèn)使用。該模型也是DPW-IV和DPW-V的標(biāo)模。在DPW-IV[10]和 DPW-V[11]的總結(jié)報(bào)告中指出，CFD和風(fēng)洞數(shù)據(jù)在俯仰力矩曲線上有很大的差異。可能的原因之一是風(fēng)洞試驗(yàn)?zāi)Ｐ驮跍y(cè)試中因?yàn)轱L(fēng)載荷發(fā)生變形，而CFD計(jì)算的外形是沒有考慮靜氣動(dòng)彈性變形的。因此為了更好地評(píng)估CFD水平，組委會(huì)將在ETW風(fēng)洞測(cè)量得到的模型靜氣動(dòng)彈性變形量考慮了進(jìn)來，提供了在七個(gè)迎角下的變形外形。

(a) 翼身組合體外形

(b) 翼尖的靜氣動(dòng)彈性變形

近些年來，非結(jié)構(gòu)網(wǎng)格因?yàn)槠浔容^容易處理復(fù)雜外形的特點(diǎn)得到了廣泛的應(yīng)用。對(duì)于工程外形來說，生成非結(jié)構(gòu)混合網(wǎng)格的時(shí)間比生成多塊對(duì)接結(jié)構(gòu)網(wǎng)格少很多，只用較少的人為干預(yù)，網(wǎng)格生成的效率得到極大提高。非結(jié)構(gòu)網(wǎng)格另一個(gè)吸引人的方面是可以方便地使用基于流場(chǎng)的網(wǎng)格自適應(yīng)技術(shù)[12]。目前超過一半的阻力會(huì)議參與者使用非結(jié)構(gòu)網(wǎng)格進(jìn)行計(jì)算。

本文計(jì)算使用的是自行生成的非結(jié)構(gòu)網(wǎng)格，生成過程嚴(yán)格按照組委會(huì)提供的網(wǎng)格生成指南進(jìn)行。四套逐漸加密的網(wǎng)格分別記做Grid T、Grid C、Grid M、Grid F，網(wǎng)格的信息列在表1中。最粗的一套網(wǎng)格的表面網(wǎng)格如圖2所示。

表1 網(wǎng)格量信息Table 1 Information about generated grids

圖2 粗網(wǎng)格的表面網(wǎng)格示意圖Fig.2 Surface mesh of tiny grid in the grid series

首先為了驗(yàn)證計(jì)算程序和網(wǎng)格，使用逐漸加密的四套網(wǎng)格完成定升力CL=0.5狀態(tài)下的網(wǎng)格收斂性研究。計(jì)算的外形是迎角為2.75°下的靜氣動(dòng)彈性變形模型。第二個(gè)算例是使用中等規(guī)模的網(wǎng)格完成靜氣動(dòng)彈性影響的預(yù)測(cè)。迎角從2.5°到4°，以0.25°為間隔。使用的外形是在相應(yīng)迎角下的靜氣動(dòng)彈性變形模型。計(jì)算狀態(tài)為馬赫數(shù)0.85，基于平均氣動(dòng)弦長(zhǎng)的雷諾數(shù)為5×106。

2 數(shù)值方法

本文計(jì)算使用的程序是課題組自行發(fā)展的MFlow[13]。該程序采用基于格心的有限體積法，能夠處理多種網(wǎng)格單元類型(六面體、四面體、三棱柱、金字塔以及使用幾何多重過程中融合產(chǎn)生的多面體)。單元內(nèi)使用線性重構(gòu)達(dá)到二階精度。梯度計(jì)算使用基于格點(diǎn)的格林高斯方法[14]。使用Venkatakrishnan限制器[15]來抑制間斷附近的振蕩。使用Roe格式來計(jì)算無黏通量。計(jì)算使用一階歐拉后差來推進(jìn)到收斂狀態(tài)，使用局部時(shí)間步長(zhǎng)加速收斂過程。Jacobian通量使用一階迎風(fēng)格式得到。使用多重網(wǎng)格來加速收斂。

在方管流動(dòng)中存在著不可忽視的二次渦流，這些二次流動(dòng)從拐角區(qū)域發(fā)展起來，是由拐角附近雷諾應(yīng)力的梯度誘導(dǎo)產(chǎn)生的。Spalart[16-17]指出傳統(tǒng)的Boussinesq線性渦黏假設(shè)無法準(zhǔn)確地預(yù)測(cè)雷諾應(yīng)力的各向異性，即流向、法向和側(cè)向三個(gè)方向雷諾應(yīng)力的區(qū)別，因此線性渦黏模型無法準(zhǔn)確預(yù)測(cè)方管內(nèi)的流動(dòng)。Spalart借鑒了部分代數(shù)雷諾應(yīng)力的形式，提出了雷諾應(yīng)力的二次本構(gòu)關(guān)系。新的雷諾應(yīng)力定義為：

τij,QCR=τij-Cr1[Oikτjk+Ojkτik]-

(1)

其中τij是Boussinesq線性渦黏假設(shè)給出的雷諾應(yīng)力，μt是渦黏系數(shù)，Oik是無量綱化的旋轉(zhuǎn)張量，定義如下：

(2)

Smn是應(yīng)變率張量，定義如下：

(3)

Cr1和Cr2是模型常數(shù)，分別為0.3和2.5。這樣給出的雷諾應(yīng)力τij,QCR并不只是相應(yīng)的應(yīng)變率的函數(shù)，而是同時(shí)考慮了其他方向應(yīng)變的影響，從而引入雷諾應(yīng)力的各項(xiàng)異性。

計(jì)算中使用SA一方程模型[18]，雷諾應(yīng)力分別由線性渦黏假設(shè)(后續(xù)結(jié)果中記為SA)和二次本構(gòu)關(guān)系(記為SA_QCR)給出。本文的計(jì)算從自由來流開始，假定流動(dòng)為全湍流。圖3為本文計(jì)算的殘差收斂和氣動(dòng)力變化曲線，所有計(jì)算密度的殘差都下降了兩個(gè)量級(jí)以上，氣動(dòng)力的振蕩在0.2%以內(nèi)。

圖3 連續(xù)方程的殘差和氣動(dòng)力的收斂曲線Fig.3 Convergence of density residual and force and moment coefficient

3 結(jié)果和討論

3.1 定升力系數(shù)下的網(wǎng)格收斂性研究

本節(jié)中展示了定升力條件下的網(wǎng)格收斂性，結(jié)果如圖4所示。橫坐標(biāo)是N-2/3,N代表網(wǎng)格單元總數(shù)，N-2/3代表網(wǎng)格特征尺度的平方，縱坐標(biāo)中CD代表總的阻力，CDP和CDV分別代表壓差阻力系數(shù)和摩擦阻力系數(shù)，Cm代表俯仰力矩系數(shù)。MFlow在空間為二階精度，因此如果網(wǎng)格位于收斂的漸近區(qū)域，氣動(dòng)力和力矩有線性的收斂特性。除了最粗的網(wǎng)格，升阻力和俯仰力矩有近似線性的網(wǎng)格收斂性，證實(shí)流場(chǎng)解位于漸近解區(qū)域。壓差阻力隨著網(wǎng)格加密減小，摩擦阻力增大，總阻力減小。摩擦阻力隨著網(wǎng)格變化不大，表明網(wǎng)格加密到一定程度后便可以準(zhǔn)確捕捉邊界層。

圖4 定升力CL=0.5條件下網(wǎng)格收斂性曲線Fig.4 Grid-convergence properties under the lift condition of CL=0.5

SA_QCR計(jì)算結(jié)果的規(guī)律跟SA模型的計(jì)算結(jié)果規(guī)律非常一致，在使用同一套網(wǎng)格的情況下，兩者的結(jié)果只是有一個(gè)隨網(wǎng)格變化很小的平移量。SA_QCR得到的總阻力比SA大1.7cnts左右，其中壓差阻力大2.3cnts，摩擦阻力小約0.6cnts(cnt為阻力系數(shù)單位，1cnt=0.0001)。

針對(duì)機(jī)翼不同站位的壓力分布情況，一般選取9個(gè)典型的展向站位進(jìn)行分析，如圖5所示，不同站位到對(duì)稱面的距離占機(jī)翼展長(zhǎng)的百分比分別為0.131、0.201、0.283、0.397、0.502、0.603、0.727、0.846、0.950，包含了翼根、翼中間和翼尖三個(gè)不同的區(qū)域。

如圖6和圖7所示，本文分別選取兩個(gè)典型站位展示網(wǎng)格加密和湍流模型對(duì)壓力分布的影響。該迎角下，流動(dòng)基本沒有分離，隨著網(wǎng)格加密，壓力分布變化不大。在激波附近，網(wǎng)格越密，激波越陡。計(jì)算的吸力峰比試驗(yàn)低，很可能是計(jì)算的迎角比試驗(yàn)低的緣故。從η=0.603界面開始預(yù)測(cè)的激波位置偏后。QCR修正使得激波位置稍稍前移，吸力峰更高一點(diǎn)。

可以用Cfx的分布來研究壁面的分離情況，如圖8所示。藍(lán)色的區(qū)域表示負(fù)的Cfx，用來近似表征流動(dòng)分離。QCR修正使得翼身結(jié)合處的分離區(qū)域大大減小，網(wǎng)格加密使得分離變大。但總的來說，分離泡的尺寸很小，流向長(zhǎng)度不大于當(dāng)?shù)叵议L(zhǎng)的4%。

圖8 表面流線和Cfx云圖比較Fig.8 Comparison of surface streamlines and contours of Cfx

3.2 CRM WB抖振特性預(yù)測(cè)

從圖9(a)的升力系數(shù)曲線可以看到，在迎角3.5°處，SA出現(xiàn)了比較明顯的提前的流動(dòng)分離，升力系數(shù)開始轉(zhuǎn)而下降。隨著迎角繼續(xù)增大，升力系數(shù)又開始爬升，在迎角3.75°處升力線拐折，斜率減小。相對(duì)于原始SA模型，QCR在迎角2.5°～3.0°之間對(duì)升力預(yù)測(cè)的影響相對(duì)較小。在更大的迎角下，QCR修正對(duì)計(jì)算結(jié)果影響很大。QCR預(yù)測(cè)的升力系數(shù)跟試驗(yàn)趨勢(shì)一致，只是有一個(gè)隨迎角變化不大的平移量。

圖9(b)為阻力系數(shù)的預(yù)測(cè)情況。SA_QCR計(jì)算得到的阻力系數(shù)跟試驗(yàn)的規(guī)律基本一致，比試驗(yàn)值稍大約30～40cnts。原始的SA模型在迎角3.5°開始出現(xiàn)明顯的斜率變小。在極曲線上，由于升力系數(shù)在迎角3.5°出現(xiàn)的拐折，原始SA模型出現(xiàn)了明顯的下降再上升的過程。SA_QCR結(jié)果跟試驗(yàn)趨勢(shì)一致。

同樣的，SA_QCR預(yù)測(cè)的俯仰力矩系數(shù)跟試驗(yàn)規(guī)律比較一致，如圖9(c)，但是預(yù)測(cè)的力矩過于低頭。有學(xué)者指出，考慮支架干擾會(huì)進(jìn)一步改善俯仰力矩的預(yù)測(cè)。

圖9 湍流模型對(duì)氣動(dòng)特性預(yù)測(cè)的影響Fig.9 Effect of different turbulence models on the prediction of aerodynamic properties

從迎角3.5°開始，SA計(jì)算結(jié)果背離試驗(yàn)的趨勢(shì)，而SA_QCR的結(jié)果一直跟試驗(yàn)趨勢(shì)比較吻合。從圖10(a)的3.5°迎角下的壓力分布可以看出，SA_QCR的預(yù)測(cè)大體上令人滿意。在靠近翼根的截面，原始SA模型計(jì)算出現(xiàn)了激波位置大大靠前、激波后流動(dòng)分離的現(xiàn)象。從圖10(b)整個(gè)機(jī)翼的壓力分布云圖可以清楚的看出，原始模型的結(jié)果，在靠近翼根區(qū)域，激波位置更加靠前。

圖10(c)給出了表面流線和Cfx云圖。原始的SA模型在翼身結(jié)合處出現(xiàn)了比較大的分離，而在QCR的結(jié)果中，此處的分離幾乎不可見。QCR在機(jī)翼中部激波后的分離要大得多。

圖10 在迎角3.5°下，SA與SQ_QCR湍流模型對(duì)流場(chǎng)模擬的比較Fig.10 Comparison of SA and SA_QCR turbulence models for flow field simulations at 3.5° of angle of attack

翼根處的分離同樣可以從η=0.131截面處的馬赫數(shù)云圖和平面流線清楚地看到，如圖11所示。使用原始SA模型時(shí)，流動(dòng)在流向遇到很大的逆壓梯度進(jìn)而分離，形成較大的分離區(qū)，減弱了分離點(diǎn)之前流體微團(tuán)的能量。而使用QCR之后，流動(dòng)一直保持附著流動(dòng)。

(a) SA

(b) SA_QCR

為了深入分析QCR對(duì)翼身結(jié)合處拐角流動(dòng)的影響，取機(jī)翼中弦處等流向截面來研究局部的流動(dòng)現(xiàn)象。圖12給出了機(jī)翼表面的壓力分布沿機(jī)翼展向的變化。在翼身結(jié)合處，原始SA模型得到的壓力比QCR的結(jié)果大很多。為了進(jìn)一步分析SA模型得到更大逆壓梯度的原因，圖13給出了拐角處渦量流向分量ωx的云圖，紅色和綠色分別代表符號(hào)相反的ωx區(qū)域。在原始SA模型的結(jié)果中，在拐角的機(jī)翼一側(cè)ωx為負(fù)。使用QCR后，貼近壁面處ωx為負(fù)，離開壁面一定距離后ωx為正，在兩者相互作用下，當(dāng)?shù)乇诿娴膲毫Σ蝗缭糞A模型的大。Spalart指出QCR會(huì)提高管道流動(dòng)中二次流模擬的精度。在本文算例中，與試驗(yàn)結(jié)果相比，QCR得到了更加合理的壓力分布和流動(dòng)分離情況，也提高了氣動(dòng)特性的預(yù)測(cè)。

圖12 機(jī)翼中弦處壓力分布沿機(jī)翼展向的變化Fig.12 Comparison of Cp along spanwise direction

圖13 機(jī)翼中弦處ωx分布云圖比較Fig.13 Comparison of ωxdistribution at the middle of wing

4 結(jié) 論

本文使用自行發(fā)展的MFlow求解器對(duì)第六屆AIAA阻力預(yù)測(cè)會(huì)議的翼身組合體外形進(jìn)行了數(shù)值模擬研究。

在CL=0.5條件下，數(shù)值模擬得到的氣動(dòng)力和力矩有線性的網(wǎng)格收斂特性，證實(shí)流場(chǎng)解位于漸近解區(qū)域。線性渦黏模型和考慮二次本構(gòu)關(guān)系的湍流模型的結(jié)果差別不大。但是隨著迎角增大，線性渦黏模型的結(jié)果在翼身結(jié)合處一個(gè)很大的分離渦，這也導(dǎo)致了預(yù)測(cè)的升力系數(shù)和阻力系數(shù)曲線出現(xiàn)明顯的拐折。使用考慮二次本構(gòu)關(guān)系的湍流模型計(jì)算，翼身結(jié)合處的分離渦幾乎不可見，這與試驗(yàn)更加吻合，預(yù)測(cè)的氣動(dòng)力和力矩也跟試驗(yàn)趨勢(shì)一致。

對(duì)翼身結(jié)合處局部二次流動(dòng)的詳細(xì)分析展示了考慮二次本構(gòu)關(guān)系的湍流模型提高模擬精度的原因。在使用QCR模型后，非常貼近壁面處和離開壁面一定距離處有相反的ωx，在兩者相互作用下，當(dāng)?shù)乇诿娴膲毫Σ蝗缭糞A模型大。QCR修正對(duì)二次流模擬精度的提高，得到了更加合理的壓力分布和流動(dòng)分離情況，也提高了氣動(dòng)特性的預(yù)測(cè)。

本文僅從局部流動(dòng)分析入手，有效驗(yàn)證了QCR模型在工程應(yīng)用上的潛力，下一步將會(huì)在模型層面對(duì)其進(jìn)行探討，分析規(guī)律。