張 揚，白俊強，華 俊，徐晶磊

（1．西北工業(yè)大學(xué)航空學(xué)院，陜西 西安 710072；2．中國航空研究院，北京 100012；3．北京航空航天大學(xué) 能源與動力工程學(xué)院，北京 100191）

連續(xù)尺度變換方程在穩(wěn)態(tài)及非穩(wěn)態(tài)流動模擬中的應(yīng)用

張 揚1，白俊強1，華 俊2，徐晶磊3

（1．西北工業(yè)大學(xué)航空學(xué)院，陜西 西安 710072；2．中國航空研究院，北京 100012；3．北京航空航天大學(xué) 能源與動力工程學(xué)院，北京 100191）

為了提高渦粘性假設(shè)的湍流模型對于非穩(wěn)態(tài)流動的求解精度，同時兼顧其對于穩(wěn)態(tài)流動的求解性能，將雷諾應(yīng)力項與連續(xù)變換方程（CSSE）結(jié)合而形成新的應(yīng)力項，使其根據(jù)流場尺度、網(wǎng)格尺度及Kolmogorov尺度來自動調(diào)節(jié)當(dāng)?shù)氐膽?yīng)力雷諾應(yīng)力模化水平，避免網(wǎng)格因素在流場模擬中產(chǎn)生不利影響，改正了混合RANS／LES方法的速度型偏離對數(shù)率問題；同時，該方法并未引入顯式亞格子模型（SGS），因此回避了亞格子系數(shù)確定對于流場模擬精度產(chǎn)生的影響，改善了湍流模型對于流動不穩(wěn)定性的辨識精度。在湍流平板算例中，CSSE方法計算的邊界層速度型精度與雷諾平均方法（RANS）相當(dāng)，而對于圓柱尾跡的模擬則證明了CSSE方法具有混合RANS／LES方法的優(yōu)點，即能夠準(zhǔn)確模擬流動的不穩(wěn)定性特征。

湍流模型；剪切應(yīng)力；雷諾方程；混合RANS／LES方法

0 引 言

自然界中的流動大部分以湍流形態(tài)存在，湍流的精確求解仍然是一個開放的科學(xué)問題。雖然流動控制方程早已被Navier和Stocks提出，但是人們一直以來無法得到一種適合于各種流體的解析答案。湍流的Kolmogorov尺度伴隨雷諾數(shù)的增長迅速減小，Pope［1］認(rèn)為最小的空間尺度應(yīng)是Kolmogorov尺度的2倍才能完全解析出流場的全部脈動，這引起直接數(shù)值模擬方法（DNS）所需要的網(wǎng)格數(shù)量巨大；同時，為了求解流場中小尺度脈動，所需時間步必須調(diào)至極其微小才能兼顧空間推進(jìn)和時間推進(jìn)。這些方面注定了DNS方法現(xiàn)在只能運用于低雷諾數(shù)簡單流動，如槽道流動等［2］。DNS方法求解效率低下原因在于對于能譜耗散區(qū)域的直接模擬，假設(shè)流場中耗散尺度都為各項同性，可以將這些尺度對于流場的影響量以一種解析表達(dá)式給定出來，流場中其余尺度脈動采用直接求解，這種方法被稱為大渦模擬（LES），該方法雖然在一定程度上降低了DNS方法的網(wǎng)格需求，在一些方面也有了一定的應(yīng)用［3］。但隨著雷諾數(shù)的提高，所謂的“大渦”也會減小為較小尺度的渦旋結(jié)構(gòu)，因此該方法在高雷諾數(shù)的應(yīng)用仍然存在一定的障礙。

雷諾平均方法（RANS）是求解N-S方程中最常用的統(tǒng)計平均方法，該方法通過?；鲌龅姆嵌ǔＰ?yīng)來達(dá)到降低計算代價的目的，由于該方法流場脈動的空間尺度隨時間的變化被平均運算抹平，因而面對附著流動等時間伴隨性并不十分明顯的流場求解可以獲得滿意的精度。從統(tǒng)計平均方面來說，流場脈動速度的統(tǒng)計值為零，但是在流場不穩(wěn)定因素增多時，如大分離流動，脈動項會對速度場的平均項產(chǎn)生很大的影響，繼而影響流場的動量輸運作用，因此RANS方法在伴隨分離的流動求解中精度較差。

為了獲得求解效率與精度的平衡，學(xué)者們試圖將RANS方法和LES方法進(jìn)行混合，Spalart［4］等人發(fā)展了一種脫體渦模擬方法（DES），以一方程Spalart-Allmaras湍流模型［5］為基礎(chǔ)，采用網(wǎng)格距離對原方程中的壁面距離進(jìn)行變換，使得方程毀滅源項迅速增大，從而迅速降低當(dāng)?shù)氐睦字Z應(yīng)力水平，將更多的脈動尺度釋放出來。該方法已經(jīng)得到了廣泛的應(yīng)用，Squres［6］在2002年的文獻(xiàn)中，展示了對于C-130，F(xiàn)-16及F-15E的全機DES模擬；Saqib等人采用DES方法對渦流發(fā)生器進(jìn)行了計算［7］；Shivaji等人采用DDES方法預(yù)測了三維的翼型失速過程［8］；國內(nèi)的吳晶峰等人［9］采用DES方法對后臺階流動進(jìn)行了計算，反映出該方法對于回流運動的精確捕捉；陳江濤等基于DES開展了綜合應(yīng)用分析［10］。

由于大分離流動中流動充分失穩(wěn)，附面層內(nèi)部的求解精度對于整體結(jié)果精度影響較小，但是在諸如翼型的抖振邊界等小尺度分離流場中，DES方法所得結(jié)果通常過早失速，這主要由于尺度轉(zhuǎn)換方程過于簡單，引起網(wǎng)格加密過程中附面層內(nèi)被LES侵入，當(dāng)?shù)孛}動未被充分激發(fā)而導(dǎo)致雷諾應(yīng)力迅速下降，引起當(dāng)?shù)厮俣刃团c壁面對數(shù)率不匹配，甚至引發(fā)當(dāng)?shù)爻霈F(xiàn)非物理分離現(xiàn)象［11］，所以在穩(wěn)態(tài)流動如平板流、槽道流動中無法得到滿意的求解精度。

為了克服該缺點，本文以較常用的Menterk-ω SST［12］方程為基礎(chǔ)進(jìn)行改造，使其在保留穩(wěn)態(tài)流動求解精度的同時，能夠提高對于非定常流動的求解精度，形成一種對于流動失穩(wěn)強度自適應(yīng)的湍流模型。

1 湍流控制方程

目前在工程中主要以線性兩方程湍流模型為主，兩方程湍流模型主要面向的是兩個基本量，即湍動能和耗散率，基于此發(fā)展出來了常用的k-ω［13］和k-ε［14］模型，而k-ε模型在附面層附近的模擬精度較差，因為湍流耗散率在附面層內(nèi)數(shù)值極小，對網(wǎng)格要求比較苛刻；而Wilcox［13］提出的k-ω方程雖然改善了近壁模擬精度，但是對于自由來流條件極為敏感。Menter［12］提出基于k-ω方程的剪切應(yīng)力輸運模型，將兩者的優(yōu)點通過混合函數(shù)結(jié)合起來，附面層附近采用k-ω方法，遠(yuǎn)壁面位置采用k-ε方法，保留了兩種方法各自的優(yōu)點。其方程形式為：

其中Pk和Pω分別為湍流生成項，表達(dá)式為：

模型中的Cω、σk、σω、β、F1、F2均為常數(shù)和特定的函數(shù)表達(dá)式，其數(shù)值和處理方法采用文獻(xiàn)［12］中的做法。最終，在剪切應(yīng)變率模型中渦粘性系數(shù)表達(dá)式為：

其中a1＝0.3，Ω為渦量，其張量形式為：Ω＝0.5（Ui，j－Uj，i），Ui，j和Uj，i為速度導(dǎo)數(shù)的張量。

為了與DES方法對比，本文在算例中同時引入了SST-DES方法進(jìn)行驗證，上述SST湍流模型中湍流尺度lk-ω表達(dá)式為：

將湍流尺度lk-ω替換為min（lk-ω，CDESΔ），其中Δ為網(wǎng)格三個方向的最大邊長，常數(shù)CDES＝0．78，這樣便形成了SST-DES［15］方法，具體做法可以參考文獻(xiàn)［15］。

2 基于連續(xù)尺度變換函數(shù)的應(yīng)力方程

以線性渦粘假設(shè)為背景的SST湍流模型對于雷諾平均后產(chǎn)生的不封閉項采用Boussinesq假設(shè)關(guān)系式獲得，但是該過程容易對渦粘性產(chǎn)生過高估計，加之時均效應(yīng)對流場中的中小尺度脈動會造成強烈的壓制，總體上難以反映流場的中高頻脈動變化過程，在雷諾平均動量方程中，需要進(jìn)行建模求解的目標(biāo)量為湍流剪切應(yīng)力τRANSij，該項表征著湍流脈動引起的動量通量輸運，一般通過Boussinesq提出的渦粘性理論進(jìn)行求解，該理論假設(shè)湍流應(yīng)力和平均應(yīng)變率張量中存在一種線性關(guān)系，即：

其中u′i為湍流脈動速度，k為湍動能，δij是Kronecher Delta符號（若i＝j(luò)，則δij＝1；若i≠j，則δij＝0），結(jié)合式（3），可以得到以下結(jié)論：

Batten［16］等人曾經(jīng)提出過一種數(shù)值限制尺度模型（LNS），該方法通過一個調(diào)整因子將湍流模型的生成項降低來削減當(dāng)?shù)氐睦字Z應(yīng)力水平；Spalart等人的DES方法則是將湍流模型的耗散項增大來達(dá)到同樣的目的。這兩種思路是一種湍流求解尺度轉(zhuǎn)換思想，如DES方法中，壁面距離被式（9）代換：

在本文中，為了避免網(wǎng)格誘導(dǎo)分離現(xiàn)象的產(chǎn)生，同時將流場尺度切換函數(shù)變得更為光順，引入連續(xù)型變換因子ZR，令ZR∈（0，1），并將雷諾應(yīng)力改寫為：

圖1 連續(xù)尺度變換方法與DNS和RANS方法關(guān)系Fig.1 Relationship between CSSE，DNS and RANS

區(qū)別于DES方法的是，CSSE方法并未引入顯式的亞格子模型方程，避免了亞格子系數(shù)的不同而導(dǎo)致RANS和LES區(qū)域分割不同。

3 連續(xù)變換因子ZR的推導(dǎo)

在湍流能量譜中，湍動能高波段代表湍流耗散區(qū)，該區(qū)域的湍流尺度通常位于網(wǎng)格分辨尺度以下，對于該區(qū)域在LES方法中通常采用亞格子模型模化，將一定波數(shù)以上的脈動采用關(guān)系式模化以此來反映高頻脈動對于其他頻段脈動影響，其他波數(shù)則直接求解，而將直接求解與模化求解的分界線定義為截斷波數(shù)κc，從圖2中可以看出，κc位于湍動能譜慣性子區(qū)，考慮到Kolmogorov的－5／3次方率，同時湍動能譜函數(shù)、湍流耗散率以及波數(shù)之間有如下估計：

κi和κK分別代表了湍流積分尺度和Kolmogorov尺度對應(yīng)的波數(shù)。

從本質(zhì)上講，除DNS方法外，其他湍流模擬方法仍是一個以截斷波數(shù)為界線的分區(qū)問題，以此為基礎(chǔ)構(gòu)造連續(xù)變換因子ZR，定義：

其中ED和ES分別代表直接求解的湍動能和模化求解的湍動能，將式（12）重新寫為：

為了保證數(shù)值穩(wěn)定性，令κc＝max（κi，min（κK，κc））來保證ZR的上下限。由于波數(shù)κ在求解中不能直接得出，需要將其轉(zhuǎn)換為湍流的尺度量，引入波數(shù)與湍流尺度的關(guān)系式，可以得到：

聯(lián)立式（13）與式（14），將轉(zhuǎn)換因子寫為：

其中：lK＝（ν3／ε）1／4，li＝Cμk3／2／ε，根據(jù)Nyquist采樣定律（如尺度為Δ的網(wǎng)格是無法解析尺度小于2Δ的流場脈動），故定義lc＝2Δ＝2（ΔxΔyΔz）1／3，Δx、Δy和Δz分別代表了當(dāng)?shù)鼐W(wǎng)格在三個方向的法向長度（采用網(wǎng)格體積除以三個方向投影面積獲得）。

圖2 各個尺度在能譜中的位置Fig.2 The locations of each scale at energy spectrum

4 算例驗證及分析

4.1 平板邊界層流動模擬

從圖3反映的速度型分布可以看出，由于該網(wǎng)格當(dāng)初是為高雷諾數(shù)的RANS方法模擬設(shè)計，僅僅是為了取得精確結(jié)果而進(jìn)行了一定程度的加密，DES方法的結(jié)果過早上揚，證明了該方法給出的當(dāng)?shù)販u粘性系數(shù)下降過低，說明DES方法在當(dāng)?shù)匾呀?jīng)開始采用LES方法進(jìn)行模擬而削減了當(dāng)?shù)氐睦字Z應(yīng)力水平，而CSSE方法與RANS方法與Spalding理論解十分貼近，證明CSSE方法在附面層內(nèi)確實能夠表現(xiàn)出完全的RANS特征，因此該方法對于穩(wěn)態(tài)流動的魯棒性要高于DES方法。

圖3 平板邊界層速度型Fig.3 Velocity pr of ile of flat plate boundary

4.2 圓柱尾部非穩(wěn)態(tài)流動模擬

圓柱作為經(jīng)典的鈍頭體模型，經(jīng)常被用來驗證非穩(wěn)態(tài)流動的特征，因為圓柱即使在很低的雷諾數(shù)下仍可以產(chǎn)生明顯的流動分離現(xiàn)象和渦旋結(jié)構(gòu)。

本算例采用一個處于亞臨界區(qū)域的圓柱試驗進(jìn)行模擬，實驗雷諾數(shù)Re＝3900。該雷諾數(shù)下是一個層湍共存區(qū)域，由于SST湍流模型并未引入轉(zhuǎn)捩識別模型，所以基于全湍假設(shè)的RANS方法對于該區(qū)域的模擬精度并不理想，采用該算例更能說明RANS、CSSE和DES方法的預(yù)測精度。

試驗數(shù)據(jù)采集自Lourenco［17］等人的試驗，圓柱參考長度為直徑D＝0．01m，計算域如圖。

圖4 圓柱繞流計算域Fig.4 Computation domain of flow over circle cylinder

為了使圓柱展向流動充分發(fā)展，展長給定為πD，保證第一層網(wǎng)格y＋＜1以盡可能的描述小尺度脈動。該網(wǎng)格規(guī)模約為130萬，采用多塊結(jié)構(gòu)網(wǎng)格方便并行分區(qū)計算分區(qū)。計算域入口給定速度型入口，上頂、地板和左右壁面設(shè)定為對稱邊界條件，出口給定壓力出口。圖5為網(wǎng)格示意圖。

圖5 圓柱空間和表面網(wǎng)格Fig.5 Space and surface mesh of cylinder

圖6為圓柱升力和阻力系數(shù)隨時間發(fā)展的變化，CSSE和DES方法都反映出了一種無規(guī)律的振蕩特性，該現(xiàn)象體現(xiàn)了圓柱尾跡中非穩(wěn)態(tài)流動中各尺度脈動相互影響，中小尺度脈動的高頻振蕩使得力系數(shù)的振蕩峰值各不相同。

圖6 力系數(shù)隨時間變化歷程Fig.6 Time history of force coefficient

從圖6觀察可知，兩種方法所得升力系數(shù)基本在0附近振蕩，阻力系數(shù)在1．1附近振蕩，從Lourenco的試驗中采集得到的阻力系數(shù)Cd為0．99，證明CSSE方法在非穩(wěn)態(tài)流動中的求解精度相當(dāng)于DES方法。

圖7和圖8分別為兩種方法的Q值為20時的云圖，Q準(zhǔn)則是Hunt［18］等人為了反映非穩(wěn)態(tài)流場中的瞬時相干結(jié)構(gòu)而發(fā)明的一種可視化準(zhǔn)則，Q準(zhǔn)則可以清晰的反映流場中的渦旋結(jié)構(gòu)。從圖8中可以看出，CSSE方法在釋放流場渦旋結(jié)構(gòu)方面基本等同于DES方法，這反映了該方法的類LES特征，說明連續(xù)尺度變換方程在遠(yuǎn)離附面層后迅速削減了當(dāng)?shù)乩字Z應(yīng)力水平，即保證了附面層內(nèi)的穩(wěn)態(tài)流動求解，也兼顧了非穩(wěn)態(tài)流動的求解精度。

圖7 DES方法計算的Q云圖Fig.7 Q contour of DES method

圖8 CSSE方法計算的Q云圖Fig.8 Q contour of CSSE method

圖9為平均后的圓柱表面壓力系數(shù)分布，CSSE和DES方法所預(yù)測的壓力系數(shù)與試驗數(shù)據(jù)形態(tài)基本一致，壓力拐點吻合的較好，而RANS方法在壓力拐點后的逆壓梯度遠(yuǎn)大于試驗值，說明RANS方法對于回流區(qū)預(yù)測過長。

圖10分別為X／D＝1．06、1．54和2．02站位的速度型，可以看出，CSSE方法在近壁區(qū)域和遠(yuǎn)壁區(qū)域的預(yù)測精度均與試驗值吻合，具有良好的魯棒性。

圖9 圓柱表面壓力系數(shù)Fig.9 Cpof cycle cylinder surface

圖10 圓柱不同切面速度型Fig.10 Velocity pr of ile at different slice of cylinder

5 結(jié)束語

本文通過一種連續(xù)尺度變換方程，將原始SST方程渦粘性系數(shù)重新構(gòu)造，使其變換為一種流場尺度自分辨方法，通過理論分析和算例驗證，取得了以下結(jié)果：

（1）提出了一種基于連續(xù)變換方程的湍流模擬混合方法，通過理論分析證明該方法可以在DNS和RANS間自動進(jìn)行切換；

（2）構(gòu)造了一種連續(xù)變換因子，除了當(dāng)?shù)鼐W(wǎng)格尺度外，該因子還依賴流場中的慣性尺度以及Kolmogorov尺度，避免了混合RANS／LES方法對于網(wǎng)格的嚴(yán)重依賴性；

（3）基于平板邊界層模擬，證明CSSE方法在穩(wěn)態(tài)流動中可以保持良好的RANS特性，阻止附面層內(nèi)的雷諾應(yīng)力過度損耗，保持正確的附面層內(nèi)速度型分布；

（4）通過圓柱繞流尾部的非穩(wěn)態(tài)流動模擬，證明CSSE方法有類LES方法特征，可以有效的緩解RANS方法對于高波數(shù)段脈動的壓制，其流場仿真精度基本等同于LES方法，說明這是一種新型的混合RANS／LES方法；同時，由于CSSE方法并未引入顯式的亞格子模型，因此擺脫了傳統(tǒng)混合RANS／LES方法中亞格子常數(shù)對于求解精度的影響問題。

［1］POPE S B．Turbulent flows［M］．London：Cambridge University Press，2000：231-232．

［2］JOHN K，PARVIZ M，ROBERT M．Turbulence statistics in fully developed channel flow at low Reynolds number［J］．Journal of Mech．，1987，177：133-166．

［3］PIERRE S．Large eddy simulation for incompressible flows［M］．Berlin：Springer Press，1998：131-132．

［4］SPALART P R，JOU W H，STRELETS M，et al．Comments on the feasibility of LES for wings and on a hybrid RANS／LESapproach［M］．Los Angles：Greyden Press，1997．

［5］SPALART P R，ALLMARAS S R．A one equation turbulence model for aerodynamic flows［J］．RecherchéAerospatiale，1994，1：5-21．

［6］SQUIRES K D，F(xiàn)ORSYTHE J R，MORTON S A．Progress on Detached Eddy Simulation of massively separated flows［R］．AIAA 2002-1021．

［7］SAQIB M，ROLF R．Detached-Eddy Simulation of vortex generator jet using chimera grids［J］．World Academy of Science，Engineering and Technology，2011，79：732-739．

［8］SHIVAJI M，JAMES D B．Numerical investigation of 3-D dynamic stall using delayed detached eddy simulation［R］．AIAA 2012-0099．

［9］WU J F，NING F F．Hybrid RANS／LES method applied to backward facing step［J］．Journal of Beijing University of Aeronautics and Astronautics，2011，6（1）：32-37．（in Chinese）

吳晶峰，寧方飛．后臺階流動的Hybrid RANS／LES模擬［J］．北京航空航天大學(xué)學(xué)報，2011，6（1）：32-37．

［10］CHEN J T，ZHANG P H，ZHOU N C，et al．Application of Detached-Eddy Simulation based on Spalart-Allmaras turbulence model［J］．Journal of Beijing University of Aeronautics and Astronautics，2012，7（38）：2-5（in Chinese）

陳江濤，張培紅，周乃春，等．基于SA湍流模型的DES方法應(yīng)用［J］．北京航空航天大學(xué)學(xué)報，2012，7（38）：2-5．

［11］MENTER F R，EGOROV Y．Re-visting the turbulent scale equation［C］．In：Proc．IUTAM Symp．One Hundred Years of Boundary Layer Research．Springer，Gottingen，2004．

［12］MENTER F R．Two-equation eddy-viscosity turbulence models for engineering applications［J］．AIAA Journal，1994，32（8）：1598-1605．

［13］WILCOX D C．Reassessment of the scale-determining equation for advanced turbulence models［J］．AIAA Journal，1988，26（11）：1299-1310．

［14］JOHNSON D A，KING L S．A mathematically simple turbulence closure model for attached and separated turbulent boundary layer［J］．AIAA Journal，1985，23（11）：1684-1692．

［15］STRELETS M．Detached Eddy simulation of massively separated flows［R］．AIAA 2001-0879．

［16］BATTEN P，GOLDBERG U，CHAKRAVARTHY S．LNS-an approach towards embedded LES［R］．AIAA2002-0427．

［17］LOURENCO L M，KROTHAPALLI A．Application of DPIV to the study of the temporal evolution of the flow past a circular cylinder［J］．Laser Anemometry in Fluid Mech．，1988：161-177．

［18］HUNT J，WRAY A，MOIN P．Proceedings of the 1988 summer program［R］．Stanford：Center for Turbulence Research，1988．

Application of a continuous scale switch equation in steady and unsteady flow simulation

ZHANG Yang1，BAI Junqiang1，HUA Jun2，XU Jinglei3
（1.School of Aeronautics，Northwestern Polytechnical University，Xi′an 710072，China；2.Chinese Aerodynamic Establishment，Beijing 100012，China；3.School of Energy and Power Engineering，Beihang University，Beijing 100191，China）

In order to get the higher precision of the unsteady flow predicted by using liner turbulence model as well as keep the capability of solving the steady flow，the Reynolds-stress term is combined with continuous scale switch equation（CSSE）to form one new stress term which can automatically adjust the local Reynolds-stress level according to the flow field scale，grid scale and Kolmogorov scale．The CSSE method eliminates the disadvantage of grid scale on the flow simulation，the problem of disagreement between the velocity pr of ile obtained by the method of hybrid RANS／LES and the law of the wall is solved．Meanwhile，this method do not use the explicit sub-grid stress（SGS）model which can avoid the disadvantaged impact on the flow simulation caused by SGS coefficient，so the flow instability can be controlled more accurately．The simulations of the boundary flow of the plate and the flow over the cylinder both prove that CSSE method is robust．In the testcase of a turbulent flat-plate，the velocity pr of iles of flat-plate boundary layer solved by CSSE match with that solved by Reynolds averaged Navier-Stocks（RANS）method very closely，the simulation of cylinder wake proves that the CSSE method can precisely simulate the flow instability which is the advantage of hybrid RANS／LES．

turbulence model；shear stress；Reynolds equation；hybrid RANS／LES method

V211

Adoi:10.7638／kqdlxxb-2012.0197

0258-1825（2014）05-0694-06

2012-11-26；

2013-03-28

國家自然科學(xué)基金（11002014）

張揚（1988-），男，陜西西安人，博士研究生，研究方向為工程湍流模擬．E-mail：iamvip2＠163．com

張揚，白俊強，華俊，等．連續(xù)尺度變換方程在穩(wěn)態(tài)及非穩(wěn)態(tài)流動模擬中的應(yīng)用［J］．空氣動力學(xué)學(xué)報，2014，32（5）：694-699．

10．7638／kqdlxxb-2012．0197． ZHANG Y，BAI J Q，HUA J，et al．Application of a continuous scale switch equation in steady and unsteady flow simulation［J］．ACTA Aerodynamica Sinica，2014，32（5）：694-699．