高 歌 江立軍 高 琳

(北京航空航天大學(xué) 能源與動(dòng)力工程學(xué)院，北京100191)

翼體角偶流動(dòng)在空氣動(dòng)力學(xué)和水動(dòng)力學(xué)研究應(yīng)用中廣泛存在，例如，飛行器外流，葉輪機(jī)械，潛艇，電子元件冷卻以及河流/橋梁流動(dòng)等.它是一個(gè)強(qiáng)三維湍流邊界層流動(dòng)，有三維邊界層的分離和再附;同時(shí)具有強(qiáng)三維各向異性，伴有多個(gè)馬蹄渦的形成和發(fā)展，具有不同長(zhǎng)度尺度的平板邊界層和翼型邊界層的相互干擾以及非定常特征等.正是因?yàn)檫@種湍流流動(dòng)的復(fù)雜性，針對(duì)它的實(shí)驗(yàn)研究直到20世紀(jì)80年代才開(kāi)始，詳盡的實(shí)驗(yàn)數(shù)據(jù)直到90年代才由文獻(xiàn)［1］給出.三維湍流流動(dòng)的復(fù)雜性，使得以往湍流模式的評(píng)估與校正通常采用相對(duì)簡(jiǎn)單的二維邊界層算例，這就導(dǎo)致將這些湍流模式應(yīng)用于三維湍流流動(dòng)時(shí)，給出的結(jié)果往往差強(qiáng)人意.文獻(xiàn)［2］采用12種湍流模式，其中包括線性、非線性渦粘模式，雷諾應(yīng)力湍流模式等，對(duì)翼體角偶流動(dòng)進(jìn)行了數(shù)值計(jì)算，得到的結(jié)果與實(shí)驗(yàn)對(duì)比相差甚遠(yuǎn)，其中雷諾應(yīng)力模型表現(xiàn)最好，尤其是在尾跡區(qū)的馬蹄渦結(jié)果捕捉較好，但沒(méi)有一個(gè)模式得到與實(shí)驗(yàn)數(shù)據(jù)接近的結(jié)果，無(wú)論是平均流場(chǎng)還是湍流場(chǎng).文獻(xiàn)［3］用RANS/LES混合方法計(jì)算了翼體角偶流動(dòng)，結(jié)果表明DES(Detached Eddy Simulation)預(yù)測(cè)分離過(guò)早，DDES(Delayed-Detached Eddy Simulation)結(jié)果稍有改善，但還是無(wú)法令人滿意.Gao-Yong湍流方程［4］從Navier-Stokes方程出發(fā)，基于側(cè)偏統(tǒng)計(jì)平均方法，在推導(dǎo)過(guò)程中不曾舍棄任何高階湍流脈動(dòng)量，?；^(guò)程中又采用了符合湍流物理實(shí)質(zhì)的動(dòng)量傳輸鏈概念，因而很好的保存了N-S方程的均化非線性特性.其獨(dú)有的正交各向異性假設(shè)，使得該方程在強(qiáng)各向異性流動(dòng)的數(shù)值模擬中表現(xiàn)優(yōu)異.本文將Gao-Yong湍流方程添加進(jìn)開(kāi)源計(jì)算流體力學(xué)軟件OpenFOAM中，并用于翼體角偶流動(dòng)的計(jì)算，給出了令人滿意的結(jié)果.OpenFOAM是一個(gè)用于計(jì)算流體力學(xué)和結(jié)構(gòu)分析的C++類(lèi)庫(kù)，集中了CFD(Computational Fluid Dynamics)中可能用到的全部元素，比如網(wǎng)格類(lèi)、場(chǎng)類(lèi)、時(shí)間類(lèi)、矩陣類(lèi)、物理模型類(lèi)等.正如文獻(xiàn)［5］中介紹的，它采用操作符形式的方法描述偏微分方程的有限體積離散化，支持多面體網(wǎng)格，可以處理復(fù)雜的幾何外形.經(jīng)過(guò)十多年的發(fā)展，現(xiàn)在OpenFOAM已經(jīng)成為開(kāi)源CFD社區(qū)最受歡迎的開(kāi)源CFD軟件，是一個(gè)極好的研究和應(yīng)用平臺(tái)，正受到越來(lái)越廣泛的關(guān)注［6］.由于C++的面向?qū)ο蠛湍K化特性，使得OpenFOAM具有極強(qiáng)的擴(kuò)展性和可維護(hù)性，本文的主要工作即是將Gao-Yong湍流方程加入其湍流模式模塊中，并進(jìn)行了相關(guān)的驗(yàn)證和數(shù)值計(jì)算工作.這一工作將極大地拓展Gao-Yong湍流方程工程應(yīng)用的前景.

1 Gao-Yong湍流方程

側(cè)偏平均思想初步形成于20世紀(jì)90年代，經(jīng)過(guò)十幾年的發(fā)展逐漸成熟:依據(jù)一定的原則將一點(diǎn)處湍流脈動(dòng)劃分為二組，每組的統(tǒng)計(jì)平均量均不為0，將各組的概率權(quán)重引入側(cè)偏平均值之后獲得了加權(quán)漂移速度的對(duì)稱性，從而明確了各組側(cè)偏平均量及其方程的對(duì)稱性，為進(jìn)一步引入加權(quán)漂移速度的正交各向異性奠定了基礎(chǔ).采用側(cè)偏平均后獲得了一階統(tǒng)計(jì)平均量——漂移速度的獨(dú)立動(dòng)量方程和機(jī)械能方程，并利用唯象的動(dòng)量傳輸鏈概念解決了封閉問(wèn)題.二維GAO-YONG不可壓湍流方程詳見(jiàn)文獻(xiàn)［4］，其最終完整的三維不可壓湍流方程在文獻(xiàn)［5］中直接給出:

以上方程分別為平均流連續(xù)方程、動(dòng)量方程、漂移流連續(xù)方程、動(dòng)量方程以及機(jī)械能方程(漂移量有波浪上標(biāo)).式中分別為平均流的層流應(yīng)力和漂移流的層流應(yīng)力，遵循廣義牛頓定律所規(guī)定的應(yīng)力和變形率關(guān)系分別為平均流湍流應(yīng)力和漂移流湍流應(yīng)力.詳細(xì)的推導(dǎo)過(guò)程參見(jiàn)文獻(xiàn)［7］，本文不再贅述.

2 計(jì)算模型與計(jì)算方法

本文研究對(duì)象是橢率為3∶2的橢圓頭與NACA0020翼尾在最大厚度處結(jié)合的聯(lián)合體和平板組成的角偶流動(dòng).此突體模型最大厚度為7.17 cm，弦長(zhǎng)30.5 cm，高22.9 cm.自由來(lái)流平均速度27 m/s.流動(dòng)雷諾數(shù) Re0=U0T/ν=115 000，其中U0為自由來(lái)流平均速度，T為模型最大厚度，分別作為參考速度和參考長(zhǎng)度.幾何模型及采用的坐標(biāo)系如圖1所示.在計(jì)算中，取翼型前緣與平板相交處為坐標(biāo)原點(diǎn).x為流動(dòng)方向，y沿展向，z垂直于平板方向.進(jìn)口設(shè)在x/T=－11.2，進(jìn)口邊界條件、速度、湍流強(qiáng)度等均采用ERCOFTAC數(shù)據(jù)庫(kù)的實(shí)驗(yàn)數(shù)據(jù).出口設(shè)在x/T=14.0處，邊界為沿流向零梯度.y方向的計(jì)算區(qū)域?yàn)閥/T=－5到y(tǒng)/T=5，此兩側(cè)邊界均采用可滑移條件.z方向的計(jì)算區(qū)域?yàn)閦/T=0到z/T=3.2.上邊界同樣采用可滑移邊界條件.此外，平板和翼型突體均為固體壁面，采用無(wú)滑移邊界條件.

在x-y平面內(nèi)用191×96的C型網(wǎng)格離散.在垂直平板的z方向，用56個(gè)網(wǎng)格離散.本文利用OpenFOAM現(xiàn)有求解器，采用空間離散的有限體積法、時(shí)間離散的隱式離散方法.因Gao-Yong湍流方程中的漂移流連續(xù)方程和動(dòng)量方程形與平均流基本一致，故數(shù)值計(jì)算方法均采用SIMPLE算法，對(duì)流項(xiàng)離散格式采用二階精度的TVD格式離散，擴(kuò)散項(xiàng)采用中心格式離散.為避免使用非交錯(cuò)網(wǎng)格所引起的非物理壓力波動(dòng)，對(duì)有限體積界面上的速度采用了Rhie-Chow動(dòng)量插值.離散的線性方程組采用預(yù)處理雙共軛梯度法求解.

圖1 計(jì)算幾何模型

3 計(jì)算結(jié)果

實(shí)際上，翼體角偶流動(dòng)是由兩個(gè)邊界層流動(dòng)相互干涉而形成的一類(lèi)復(fù)雜三維流動(dòng).其流動(dòng)機(jī)理如下:平板遠(yuǎn)前方來(lái)流邊界層受到因突體產(chǎn)生的逆壓梯度影響，發(fā)生分離并繞自身卷起，渦團(tuán)產(chǎn)生.當(dāng)這個(gè)渦團(tuán)隨著主流逐漸接近突體時(shí)會(huì)拉伸變形，圍繞突體形成所謂的馬蹄渦.因?yàn)檫@種渦是由剪切層歪斜產(chǎn)生，根據(jù)普朗特的理論這可視為第1類(lèi)二次流.許多研究者證實(shí)這個(gè)馬蹄渦的強(qiáng)度和位置與突體前緣形狀有關(guān)，一般說(shuō)來(lái)，突體越鈍，馬蹄渦強(qiáng)度越大，位置越靠前.圖2所示為使用Q準(zhǔn)則描述的馬蹄渦基本形狀.圖3給出了計(jì)算出的翼型頭部對(duì)稱面(x-z平面)內(nèi)的流線分布，從圖中可清晰地分辨出對(duì)稱面內(nèi)的回流旋渦，其中心大致位于(－0.21T，0.045T)，和實(shí)驗(yàn)測(cè)量結(jié)果(－0.2 T ，0.05 T)很接近.而且文獻(xiàn)［1］中強(qiáng)調(diào)的在平板和翼型前緣相交處的極小區(qū)域內(nèi)出現(xiàn)的二次分離渦，也被精確的捕捉到，見(jiàn)圖3.

圖4展示了平板和翼型表面的壓力分布以及尾跡區(qū)的旋渦.結(jié)合以上兩圖可以得到翼體角偶流動(dòng)的大體印象.從圖4可以看到，由于計(jì)及Gao-Yong湍流方程中機(jī)械能方程的二階項(xiàng)，尾跡區(qū)兩側(cè)旋渦相互干擾引起的非定常流動(dòng)破壞了原本的對(duì)稱性，使計(jì)算在一定程度上具備了模擬非定常流動(dòng)的能力.

圖2 使用Q準(zhǔn)則得到的馬蹄渦

圖3 翼型前對(duì)稱面流線圖

圖4 壁面壓力分布及尾跡區(qū)旋渦

本文采用的實(shí)驗(yàn)數(shù)據(jù)全部出自歐洲流體、湍流及燃燒研究協(xié)會(huì)對(duì)外公開(kāi)的數(shù)據(jù)庫(kù)(ERCOFTAC database:European Research Community ofn Flow，Turbulence and Combustion)，實(shí)驗(yàn)的相關(guān)情況在文獻(xiàn)［1］中已做了詳細(xì)介紹，本文不再贅述.圖5和圖6是不同高度處翼型表面的壓力系數(shù)圖，兩者的位置分別為z/T=0.133和z/T=0.398.由這兩圖結(jié)合圖4可以看出，翼型表面壓力在最大厚度x/T=0.6以前，壓力迅速減小.在這個(gè)位置之后，壓力又漸緩升高.在翼型表面上，大部分區(qū)域是按這個(gè)一維規(guī)律變化的，但在最大厚度位置，受馬蹄渦影響，沿垂直平板方向，越靠近平板，壓力漸緩升高.如圖5和圖6所示，在翼型表面取兩組數(shù)據(jù)與實(shí)驗(yàn)進(jìn)行對(duì)比，兩者定量符合得非常好.

圖5 z/T=0.133處翼型表面壓力系數(shù)分布

圖6 z/T=0.398處翼型表面壓力系數(shù)分布

圖7和圖8是翼型上游對(duì)稱面上的速度型.圖7為翼型上游不同位置沿流向速度U的分布，圖8為z向速度W的分布.其坐標(biāo)分別為x/T=－0.46，－ 0.4，－ 0.35，－ 0.3，－ 0.25，－ 0.2，－0.14，－0.1，－0.05.分離點(diǎn)上游和其附近的平均速度型與處于逆壓梯度下二維邊界層分離相似.當(dāng)然，這種相似只是定性上的，邊界層的增長(zhǎng)被流體沿展向(y方向)的運(yùn)動(dòng)所限制.分離點(diǎn)下游不遠(yuǎn)處，平均回流只存在于臨近壁面的狹長(zhǎng)區(qū)域.計(jì)算結(jié)果和實(shí)驗(yàn)測(cè)量符合得非常好.從分離點(diǎn)位置，到回流旋渦的強(qiáng)度，都極好地重現(xiàn)了實(shí)驗(yàn)結(jié)果.在回流區(qū)，回流速度 U最大可達(dá)0.48Uref，計(jì)算結(jié)果基本重現(xiàn)了這一結(jié)論.法向速度W最大強(qiáng)度隨著流動(dòng)接近翼型突體逐漸增大，最高可達(dá)參考速度的30%，這個(gè)趨勢(shì)也很好的捕捉到了.如圖8所示，計(jì)算出的W速度型與實(shí)驗(yàn)數(shù)據(jù)基本重合，說(shuō)明Gao-Yong湍流方程對(duì)這種復(fù)雜的三維流動(dòng)能給出精確結(jié)果.圖9是在x/T=0.76截面的y向速度分布.測(cè)量點(diǎn)坐標(biāo)分別為:y/T=0.775，0.85，0.925，1.0，1.075，1.175，1.325，1.525.此處位于翼型最大厚度(x/T=0.6)的下游不遠(yuǎn)處，在靠近翼型突體附近是馬蹄渦充分發(fā)展的區(qū)域，流動(dòng)非常復(fù)雜.從圖9可以看出，計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)量結(jié)果基本符合.

圖7 翼型上游流向速度型

圖8 翼型上游z向速度型

圖9 x/T=0.76截面不同位置y向速度分布

4 結(jié)論

本文將三維不可壓Gao-Yong湍流方程加入開(kāi)源計(jì)算流體力學(xué)軟件OpenFOAM中，并使用其對(duì)三維復(fù)雜邊界層流動(dòng)——翼體角偶流動(dòng)進(jìn)行了數(shù)值計(jì)算.計(jì)算完美捕捉到了此流動(dòng)中的主要流動(dòng)結(jié)構(gòu).無(wú)論是馬蹄渦的位置，還是馬蹄渦強(qiáng)度都得到精確結(jié)果.此外，對(duì)于固體壁面的壓力分布，回流區(qū)及垂直流向截面的速度分布，數(shù)值計(jì)算都與實(shí)驗(yàn)結(jié)果符合得極好，證明Gao-Yong湍流方程有能力精確模擬翼體角偶流動(dòng)這類(lèi)復(fù)雜的三維邊界層流動(dòng).由于整個(gè)計(jì)算基于OpenFOAM這個(gè)強(qiáng)大的開(kāi)發(fā)平臺(tái)，使得網(wǎng)格劃分，結(jié)果分析等前后處理過(guò)程極為方便，這就將Gao-Yong湍流方程向工程實(shí)際應(yīng)用方向推進(jìn)了一大步.

References)

［1］Devenport W J，Simpson R L.Time-dependent and time-averaged turbulence structure near the nose of a wing-body junction ［J］.Journal of Fluid Mechanics，1990，210:23 －55

［2］Apsley D D，Leschziner M A.Investigation of advanced turbulence models for the flow in a generic wing-body junction ［J］.Flow，Turbulence and Combustion，2001，67:25 －55

［3］Fu Song，Xiao Zhixiang，Chen Haixin，et al.Simulation of wingbody junction flows with hybrid RANS/LES methods［J］.International Journal of Heat and Fluid Flow，2007，28:1379 －1390

［4］Gao Ge，Yong Yan.Partial-average-based equations of incompressible turbulent flow ［J］.International Journal of Non-Linear Mechanics，2004，39(9):1407 －1419

［5］Weller H，Tabor G，Jasak H，et al.A tensorial approach to computational continuum mechanics using object-oriented techniques［J］.Computers in Physics，1998，12(6):620 －631

［6］江立軍，高歌，張常賢.三維斜掠后臺(tái)階數(shù)值模擬［J］.北航學(xué)報(bào)，2010，36(1):87 －90 Jiang Lijun，Gao Ge，Zhang Changxian.Numerical simulation of a swept backward-facing step flow［J］.Journal o f Beijing University of Aeronautics and Astronautics，2010，36(1):87 － 90(in Chinese)

［7］高歌，熊焰.側(cè)偏平均湍流方程研究綜述［J］.中國(guó)力學(xué)文摘，2008，22(2):1 －20 Gao Ge，Xong Yan.A review on the research development of Gao-Yong equations of turbulent flows［J］.Chinese Mechanics Abstracts，2008，22(2):1 － 20(in Chinese)