使用Gao－Yong湍流方程計(jì)算翼體角偶流動

高 歌 江立軍 高 琳

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

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

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

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

本文研究對象是橢率為3∶2的橢圓頭與NACA0020翼尾在最大厚度處結(jié)合的聯(lián)合體和平板組成的角偶流動.此突體模型最大厚度為7.17 cm，弦長30.5 cm，高22.9 cm.自由來流平均速度27 m/s.流動雷諾數(shù) Re0=U0T/ν=115 000，其中U0為自由來流平均速度，T為模型最大厚度，分別作為參考速度和參考長度.幾何模型及采用的坐標(biāo)系如圖1所示.在計(jì)算中，取翼型前緣與平板相交處為坐標(biāo)原點(diǎn).x為流動方向，y沿展向，z垂直于平板方向.進(jìn)口設(shè)在x/T=－11.2，進(jìn)口邊界條件、速度、湍流強(qiáng)度等均采用ERCOFTAC數(shù)據(jù)庫的實(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.上邊界同樣采用可滑移邊界條件.此外，平板和翼型突體均為固體壁面，采用無滑移邊界條件.

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

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

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

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

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

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

圖3 翼型前對稱面流線圖

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

本文采用的實(shí)驗(yàn)數(shù)據(jù)全部出自歐洲流體、湍流及燃燒研究協(xié)會對外公開的數(shù)據(jù)庫(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以前，壓力迅速減小.在這個位置之后，壓力又漸緩升高.在翼型表面上，大部分區(qū)域是按這個一維規(guī)律變化的，但在最大厚度位置，受馬蹄渦影響，沿垂直平板方向，越靠近平板，壓力漸緩升高.如圖5和圖6所示，在翼型表面取兩組數(shù)據(jù)與實(shí)驗(yàn)進(jìn)行對比，兩者定量符合得非常好.

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

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

圖7和圖8是翼型上游對稱面上的速度型.圖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)然，這種相似只是定性上的，邊界層的增長被流體沿展向(y方向)的運(yùn)動所限制.分離點(diǎn)下游不遠(yuǎn)處，平均回流只存在于臨近壁面的狹長區(qū)域.計(jì)算結(jié)果和實(shí)驗(yàn)測量符合得非常好.從分離點(diǎn)位置，到回流旋渦的強(qiáng)度，都極好地重現(xiàn)了實(shí)驗(yàn)結(jié)果.在回流區(qū)，回流速度 U最大可達(dá)0.48Uref，計(jì)算結(jié)果基本重現(xiàn)了這一結(jié)論.法向速度W最大強(qiáng)度隨著流動接近翼型突體逐漸增大，最高可達(dá)參考速度的30%，這個趨勢也很好的捕捉到了.如圖8所示，計(jì)算出的W速度型與實(shí)驗(yàn)數(shù)據(jù)基本重合，說明Gao-Yong湍流方程對這種復(fù)雜的三維流動能給出精確結(jié)果.圖9是在x/T=0.76截面的y向速度分布.測量點(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ū)域，流動非常復(fù)雜.從圖9可以看出，計(jì)算結(jié)果與實(shí)驗(yàn)測量結(jié)果基本符合.

圖7 翼型上游流向速度型

圖8 翼型上游z向速度型

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

4 結(jié)論

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