高 歌 江立軍 高 琳

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

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

1 Gao-Yong湍流方程

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

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

2 計算模型與計算方法

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

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

圖1 計算幾何模型

3 計算結果

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

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

圖2 使用Q準則得到的馬蹄渦

圖3 翼型前對稱面流線圖

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

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

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

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

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

圖7 翼型上游流向速度型

圖8 翼型上游z向速度型

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

4 結論

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

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］.北航學報，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］高歌，熊焰.側偏平均湍流方程研究綜述［J］.中國力學文摘，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)