蔡曉偉，譚俊杰，王園丁，任登鳳，WANG Huasheng，石 清

（1.南京理工大學(xué) 能源與動力工程學(xué)院，江蘇 南京 210094；2.倫敦瑪麗皇后大學(xué) 材料科學(xué)與工程學(xué)院，英國 倫敦 E14NS；3.中國空氣動力研究與發(fā)展中心，四川 綿陽 621000）

0 引 言

隨著計算機計算能力的不斷提高，計算流體力學(xué)（CFD）在實際工程應(yīng)用中扮演的角色也越來越重要，同時也面臨著越來越復(fù)雜的問題。高質(zhì)量網(wǎng)格（計算精度高、計算效率高、收斂性好）的生成一直是耗費人力和機時的難題。無網(wǎng)格方法［1］摒棄了傳統(tǒng)CFD方法的網(wǎng)格概念，其計算區(qū)域只涉及離散點的布置，因而越來越受歡迎。

在模擬湍流甚至是層流流動時，由于近壁面附近高度各向異性點云的存在，傳統(tǒng)無網(wǎng)格法很難準(zhǔn)確地擬合空間導(dǎo)數(shù)。Munikrishna在近壁區(qū)采用格林－高斯公式代替最小二乘法，耦合BL（Baldwin－Lomax）湍流模型，模擬了包括兩段翼型在內(nèi)的一系列流動問題［2］。王剛等采用相似的方法，對高雷諾數(shù)下的繞流進(jìn)行了模擬，取得了較好的結(jié)果［3］。點云重構(gòu)技術(shù)是在局部計算區(qū)域?qū)Ω叨雀飨虍愋缘狞c云進(jìn)行轉(zhuǎn)化重構(gòu)，從而形成適合于移動最小二乘無網(wǎng)格法的近各向同性的點云結(jié)構(gòu)，已在層流的數(shù)值模擬中取得了較好的效果。

在工程湍流問題中得到廣泛應(yīng)用的是渦粘性模型。BL模型對大多數(shù)的附體流動和弱分離流動都具有較高的準(zhǔn)確性和可靠性，但對于大多數(shù)分離流動卻很難準(zhǔn)確模擬。k－ω模型是最為人所知和應(yīng)用最為廣泛的兩方程渦粘性模型之一，該模型主要求解湍動能k及其比耗散率ω的對流輸運方程，其主要貢獻(xiàn)來自于 Wilcox［4］。且k－ω 型模型對自由剪切湍流、附著邊界層湍流和適度分離湍流都有著較高的精度［5］。

本文借鑒k－ω型湍流模型在網(wǎng)格方法中的成功應(yīng)用，通過對局部計算區(qū)域的高度各向異性點云進(jìn)行重構(gòu)，提出了一套基于無網(wǎng)格法的湍流計算方案。在此基礎(chǔ)上，通過自主編程對廣泛應(yīng)用的兩種k－ω型湍流模型（k－ω SST［6］和 k－ω TNT［7］）方程與雷諾平均N－S方程進(jìn)行解耦計算，比較分析了兩種模型在無網(wǎng)格方法中的湍流預(yù)測能力。

1 控制方程

使用自由來流參數(shù)（如：密度，聲速，粘性系數(shù)）及特征長度（如：翼型弦長）對控制方程無量綱化，則二維雷諾平均N－S方程可以表示為：

式中，W 為守恒變量，f和g為無粘通量，F(xiàn)和G為粘性通量。其具體表達(dá)形式如下：

在粘性通量表達(dá)式中，剪切應(yīng)力和熱通量可以表示為：

本文數(shù)值計算時，取γ＝1.4；Pr＝0.72；PrT＝0.90。層流粘性系數(shù)μ利用Sutherland公式計算得到，湍流粘性系數(shù)μt由湍流模型提供。

2 湍流模型

k－ωSST模型在近壁面采用k－ωWilcox模型，在邊界層邊緣和自由剪切層采用標(biāo)準(zhǔn)k－ε模型，其間通過一個混合函數(shù)來轉(zhuǎn)換。k－ωTNT模型通過假定湍流區(qū)與非湍流區(qū)的界面存在一個弱解，從而推導(dǎo)出一套不同于k－ωWilcox模型的參數(shù)。兩種湍流模型都是基于k－ω湍流模型的改進(jìn)，在保證對近壁面湍流預(yù)測的準(zhǔn)確性的同時，通過不同的途徑解決了k－ω模型對于自由來流條件的依賴性。為方便比較異同，本文列出了兩種模型的表達(dá)式：

β＊為模型系數(shù)，取β＊＝0.09。CD為交叉擴散項，兩種模型處理形式不同，將在下文中詳細(xì)表述。Pk、Pω分別為湍動能k和比耗散率ω的生成項，其表達(dá)式為：

2.1 k－ωSST模型

k－ωSST模型與一般的兩方程模型最大的不同之處在于混合函數(shù)的應(yīng)用。交叉擴散項CD以及湍流粘性系數(shù)μt的表達(dá)式如下：

混合函數(shù)的構(gòu)成以及其他參數(shù)值詳見文獻(xiàn)［6］。

2.2 k－ω TNT模型

k－ωTNT模型的交叉擴散項CD為：

k－ωTNT模型中的具體參數(shù)詳見文獻(xiàn)［7］。

通過對兩種湍流表達(dá)式的比較，可以看出不同的處理思維：k－ωSST利用到壁面距離的差異，通過混合函數(shù)實現(xiàn)兩種湍流模型（k－ω Wilcox與k－εStandard）的過渡，原理簡單，但形式表達(dá)復(fù)雜，計算量大。k－ωTNT采用不同于原始k－ω模型的參數(shù)，原理較為復(fù)雜，但表達(dá)式更為簡單，計算效率高。

3 數(shù)值方法

無網(wǎng)格算法根據(jù)離散控制方程的形式可以大致分為積分型和非積分型。積分型無網(wǎng)格法用數(shù)值積分來求解控制微分方程的弱解形式，例如：無網(wǎng)格伽遼金法（EFG）、擴散單元法（DEM）以及hp－clouds法等。非積分型無網(wǎng)格法根據(jù)節(jié)點信息進(jìn)行偏導(dǎo)數(shù)的擬合求解。移動最小二乘法擬合空間導(dǎo)數(shù)是目前在計算流體力學(xué)中得到廣泛應(yīng)用的無網(wǎng)格方法之一［8］。

3.1 移動最小二乘空間導(dǎo)數(shù)擬合

在計算區(qū)域內(nèi)，根據(jù)需要布置有限個離散點，并將其標(biāo)記為i＝1，2，…，np。對于每一個離散點，可以搜索出它的衛(wèi)星點，記為Ci，如圖1所示。與網(wǎng)格方法不同的是，無網(wǎng)格的計算模板沒有特殊的空間拓?fù)浣Y(jié)構(gòu)，只存在點與點之間的聯(lián)系。使用最小二乘法擬合中心點i的空間導(dǎo)數(shù)，需要滿足以下的兩個條件：

（1）模板中衛(wèi)星點的數(shù)目需大于空間維數(shù)；

（2）模板中所有衛(wèi)星點不能共處一條直線上。

假設(shè)u（x，y）為計算區(qū)域內(nèi)任一空間離散點的流場值，則可以基于泰勒級數(shù)展開：

空間導(dǎo)數(shù)由移動最小二乘法則決定：

式（22）的解可以寫成下面的矩陣形式：

對矩陣進(jìn)行顯式求解，并寫成下面的形式：

系數(shù)αij、βij具體表達(dá)由文獻(xiàn)［9］給出。

圖1 各向同性點云結(jié)構(gòu)示意圖Fig.1 Isotropic structure of cloud of points

圖2 通量求解示意圖Fig.2 Illustration for flux calculation

對于式（1）的對流項，可以采用多種成熟的迎風(fēng)格式計算通量。鑒于本文算例的計算范圍是從不可壓流動到可壓縮跨聲速流動，AUSM＋－up［10］格式除保留原始AUSM格式間斷、激波、粘性分辨率好的優(yōu)勢外，通過改進(jìn)原始AUSM格式的壓力項和數(shù)值聲速項，使其能夠適應(yīng)更廣馬赫數(shù)范圍的模擬。因此，本文的數(shù)值通量（圖2中的fij與gij）用AUSM＋－up格式計算。式（1）中的粘性通量中包含了速度和溫度的梯度，這兩種梯度可以通過式（24）直接計算出。求點ij的粘性通量時，為避免產(chǎn)生奇偶失聯(lián)，本文對代數(shù)平均的梯度值進(jìn)行修正。

利用最小二乘無網(wǎng)格法空間離散后，式（1）變成一個常微分方程。此時，時間推進(jìn)求解格式可以借鑒網(wǎng)格算法中成熟的應(yīng)用。本文選擇三階SSP（Strong Stability Preserving）型Runge－Kutta格式，并使用局部時間步長推進(jìn)求解。

3.2 點云重構(gòu)方法

對于一般性湍流流動，物面法向的流場梯度變化要比切向的更為劇烈。此時，為減少計算代價，一般在物面法向上的離散點間的距離比切向要小很多，甚至是10－3以上量級的差距。這樣，就會導(dǎo)致在局部計算區(qū)域出現(xiàn)高度各向異性的點云結(jié)構(gòu)，而式（23）的條件系數(shù)會達(dá)到105量級。這時的矩陣系統(tǒng)是病態(tài)的，很難準(zhǔn)確地進(jìn)行空間導(dǎo)數(shù)擬合。

只有大幅度地減少矩陣條件數(shù)，才能使得流場得以穩(wěn)定地、準(zhǔn)確地迭代到收斂值。點云重構(gòu)技術(shù)利用了無網(wǎng)格方法布點靈活的特點，對各向異性點云結(jié)構(gòu)中的衛(wèi)星點進(jìn)行增刪，從而使得矩陣條件數(shù)接近1（健康矩陣系統(tǒng)的條件數(shù)）。如圖3所示，點i為計算點，衛(wèi)星點1、2、3、4、5、6、7、8形成了各向異性的點云結(jié)構(gòu)（衛(wèi)星點2、7到計算點i的距離為其余衛(wèi)星點的10－3倍）。假設(shè)衛(wèi)星點離計算點i最短距離為R，以R＋δ（δ為極小正數(shù)）為半徑，以離散點i為原點，構(gòu)建緊支域，并與衛(wèi)星點和計算點的連線相交，從而形成新點1′、3′、4′、5′、6′、8′。加上原先的兩個點2（2′）和7（7′），這樣，一套新的近各向同性的點云結(jié)構(gòu)便形成了。

圖3 點云重構(gòu)示意圖Fig.3 Illustration for cloud of points reconstruction

采用點云重構(gòu)技術(shù)生成的新點，需要插值計算其流場值，才能參與離散點i的空間導(dǎo)數(shù)擬合。本文采用工程中常用的拉格朗日插值方法插值計算。如算例4.1及算例4.2中的離散點布置，在翼型近壁區(qū)和尾跡區(qū)的計算點會出現(xiàn)高度各向異性的點云（點云內(nèi)部最長、最短距離之比達(dá)到4×103），經(jīng)過點云重構(gòu)后，會形成健康的最小二乘矩陣系統(tǒng)。經(jīng)過數(shù)值試驗表明，此種方法不僅不會產(chǎn)生數(shù)值振蕩，而且不影響流場的迭代收斂。

4 算例分析

為驗證上述無網(wǎng)格算法的可實現(xiàn)性，本文采用該算法模擬了一系列預(yù)測湍流的標(biāo)準(zhǔn)算例。計算算例包括：繞NACA0012翼型和RAE2822翼型的亞聲速流動和跨聲速的激波誘導(dǎo)邊界層分離流動以及大分離流的后向臺階流動。對于翼型繞流，按照“C”型貼體網(wǎng)格的生成方法生成空間離散點（見圖4）。對于上述的算例，分別用耦合k－ωSST模型和k－ωTNT模型的無網(wǎng)格法進(jìn)行計算，將結(jié)果與實驗數(shù)據(jù)進(jìn)行了比對，并分析了不同情況下兩種湍流模型在無網(wǎng)格方法中的適用性。

圖4 NACA0012和RAE2822翼型離散點布置圖Fig.4 Point distributions for NACA0012and RAE2822airfoils

4.1 NACA0012翼型繞流

本算例中共計算了兩種工況：（a）Ma＝0.5，α＝3.51°，Re＝2.93×106；（b）Ma＝0.799，α＝2.26°，Re＝9.0×106。

近壁面第一層離散點距壁面1.0×10－5（以翼型弦長為無量綱參考尺度），翼型表面共布置256個離散點。遠(yuǎn)場邊界離翼型物面大于10倍弦長，整個計算區(qū)域共布置22425個離散點。

圖5 NACA0012翼型表面壓力系數(shù)分布Fig.5 Pressure coefficient distributions around NACA0012airfoil

基于上述的離散點云布置，采用無網(wǎng)格法分別耦合k－ωSST模型與k－ω TNT模型，計算了（a）、（b）兩種工況。圖5顯示了計算的表面壓力系數(shù)與文獻(xiàn)［11］的實驗結(jié)果進(jìn)行的比較?？梢园l(fā)現(xiàn)：最小二乘無網(wǎng)格法在融入點云重構(gòu)技術(shù)后，能夠克服邊界層內(nèi)空間導(dǎo)數(shù)擬合的困難，成功地耦合兩方程湍流模型。對于亞聲速的流動，兩種k－ω型湍流模型都展示出很好的預(yù)測能力，并與實驗結(jié)果吻合得較好。而對于激波誘導(dǎo)邊界分離的跨聲速流動，k－ωSST模型所預(yù)測的激波位置相比于k－ωTNT模型更接近于實驗測量結(jié)果，從而顯示出更好的激波捕捉能力。

4.2 RAE2822翼型繞流

本文的第二個算例是RAE2822翼型。分別計算了兩種工況：（a）Ma＝0.676，α＝1.92°，Re＝5.7×106；（b）Ma＝0.75，α＝2.81°，Re＝6.2×106。

近壁面第一層離散點距壁面1.0×10－5（以翼型弦長為無量綱參考尺度），翼型表面共布置304個離散點。遠(yuǎn)場邊界離翼型物面大于10倍弦長，整個計算區(qū)域共布置23985個離散點。

圖6給出了兩種工況的計算結(jié)果與實驗數(shù)據(jù)［12］及網(wǎng)格算法結(jié)果［13，15］的對比。當(dāng)計算條件為亞聲速時，兩種模型對壁面的壓力分布描述幾乎沒有差別（圖6a），下表面的壓力系數(shù)都和實驗數(shù)據(jù)吻合得很好，而上表面都與實驗數(shù)據(jù)存在一定的差別，前緣吸力峰值都較實驗數(shù)據(jù)略低。當(dāng)發(fā)生激波誘導(dǎo)邊界層分離時，k－ωSST模型捕捉激波位置的能力略強于k－ω TNT模型（k－ωSST計算出的激波位置約在52.17%弦長處，k－ωTNT計算出的位置約在54.66%弦長處，實驗測量位置約在52.32%弦長處）。另外，兩種模型對于激波后的壓力預(yù)測值都小于實驗值。文獻(xiàn)［14］基于網(wǎng)格方法耦合k－ω模型計算時也描述了此種差異。當(dāng)使用更為先進(jìn)的ν2－f四方程模型［15］時，情況依舊沒能得到改觀。作者認(rèn)為：出現(xiàn)此種差異，很可能是因為在文獻(xiàn)［13］的實驗過程中，人為地設(shè)置了轉(zhuǎn)捩點，而本文的模擬中并沒有考慮這點。

圖6 RAE2822翼型表面壓力系數(shù)分布Fig.6 Pressure coefficient distributions around RAE2822airfoil

4.3 后向臺階大分離流動

后向臺階流動作為大分離的典型流動，早在1981年Standford會議上就確定其為評價湍流模型能力的標(biāo)準(zhǔn)算例。此后，諸多學(xué)者運用包含雷諾平均N－S方程方法、大渦模擬方法（LES）和直接數(shù)值模擬（DNS）的多種湍流模型，進(jìn)行了此算例的數(shù)值模擬。本文選用該算例的目的是比較在無網(wǎng)格框架下的k－ω SST和k－ωTNT湍流模型的大分離流預(yù)測能力。算例的計算工況為：Ma＝0.128，α＝0°，Re＝37423。幾何模型如圖7所示。表1具體描述了數(shù)值模擬的條件參數(shù)。與翼型不同的是，該算例以臺階高度（H＝Hs－Hd）為參考長度，上、下表面均采用無滑移邊界條件。沿壁面法向等比例推進(jìn)布點，離散點距離壁面最近距離為1×10－3，固壁面上布點410個，整個計算區(qū)域共布點19700個。

圖7 后向臺階幾何模型Fig.7 Backward－facing step geometric

表1 后臺階模型計算參數(shù)Table1 Flow parameters for backward－facing step

圖8描述了分別采用兩種湍流模型獲得的后向臺階下表面壓力系數(shù)的分布情況，并與實驗數(shù)據(jù)［16］進(jìn)行了對比。結(jié)果表明：在附著區(qū)，兩種模型的預(yù)測能力相當(dāng)，且都與實驗結(jié)果吻合得很好。在恢復(fù)區(qū)，TNT的計算結(jié)果較實驗值偏小，SST的計算結(jié)果與實驗值吻合得較好。圖9描述了回流區(qū)表面摩擦系數(shù)的數(shù)值計算結(jié)果與實驗數(shù)據(jù)的對比。實驗測量得到的回流區(qū)長度為6.26 H，SST模型計算出的回流區(qū)長度為6.12 H，TNT模型計算出的回流區(qū)長度為6.29 H。通過四個不同位置的速度型的對比（圖10）進(jìn)一步驗證了上述的結(jié)論：在附著區(qū)（x／H＝－4.0），兩種模型獲得的速度型幾乎一致；在回流區(qū)（x／H＝2.5），TNT模型獲得的速度型略優(yōu)于SST；在恢復(fù)區(qū)（x／H＝8.0，x／H＝32.0），則反之。

圖8 下表面壓力系數(shù)比較Fig.8 Comparison of pressure coefficient along the bottom wall

圖9 回流區(qū)摩擦系數(shù)比較Fig.9 Comparison of the friction coefficient in the recirculation region

圖10 平均流速度剖面Fig.10 Streamwise mean velocity u－profiles

5 結(jié) 論

本文采用無網(wǎng)格點云重構(gòu)技術(shù)，克服傳統(tǒng)最小二乘法擬合空間導(dǎo)數(shù)的困難，成功耦合k－ωSST和k－ω TNT湍流模型，實現(xiàn)了最小二乘無網(wǎng)格法框架下的湍流流動的數(shù)值模擬。通過對弱逆壓梯度分離、激波誘導(dǎo)分離和大分離流動條件下的數(shù)值模擬，比較分析了兩種模型在無網(wǎng)格框架下的湍流預(yù)測能力。結(jié)論如下：

（1）對因弱逆壓梯度產(chǎn)生的分離流動，兩種湍流模型的預(yù)測能力相當(dāng)。而TNT模型因為表達(dá)形式簡單，因而較SST模型更具效率性。

（2）對激波誘導(dǎo)分離的流動，SST模型捕捉激波的能力更強。

（3）對于大分離的流動，在附著區(qū)和恢復(fù)區(qū)，SST的計算結(jié)果更接近于實驗值。在回流區(qū)，TNT的計算結(jié)果與實驗值吻合得更好。

綜上，本文拓寬了移動最小二乘無網(wǎng)格法在湍流研究方向的應(yīng)用，并通過對三種不同程度的分離流進(jìn)行了初步的模擬與分析，為模擬不同情況下的湍流選用合適的湍流模型提供了進(jìn)一步的參考。

［1］BATINA J.A gridless Euler／Navier－Stokes solution algorithm for complex two－dimensional application［R］.NASA TM 107631，1992.

［2］MUNIKRISHNA N，BALAKRISHNAN N.Turbulent flow computations on a hybrid Cartesian point distribution using meshless solver LSFD－U［J］.Computers ＆ Fluids，2011，40：118－138.

［3］WANG G，YE Z Y，JIANG C Q，et al.Gridless method for Navier－Stokes equations with high Reynolds number［J］.Chinese Journal of Applied Mechanics，2007，24（3）：348－352.（in Chinese）王剛，葉正寅，蔣超奇，等.高雷諾數(shù)下求解NS方程的無網(wǎng)格算法［J］.應(yīng)用力學(xué)學(xué)報，2007，24（3）：348－352.

［4］WILCOX D C.Turbulence modeling in CFD［M］.Third Edition，DCW Industries，Inc.，La Canada，CA，2006.

［5］閻超.計算流體力學(xué)方法及應(yīng)用［M］.北京航空航天大學(xué)出版社，2006.

［6］MENTER F R.Two－equation eddy－viscosity turbulence models for engineering applications［J］.AIAA Journal，1949，32（8）：269－289.

［7］KOK JC.Resolving the dependence on free－stream values for the k－ωturbulence model［J］.AIAA Journal，2000，38（7）：1292－1295.

［8］CAI X W，TAN J J，MA X J.A 2Dmeshless solver based on AUSM＋and MUSCL scheme［J］.Applied Mechanics and Materials，2012，（105－107）：2140－2143.

［9］SRIDAR D.An upwind finite difference scheme for meshless solvers［J］.Journal of Computational Physics，2003，（189）：1－29.

［10］LIOU M S.A sequel to AUSM，PartⅡ：AUSM＋－up for all speeds［J］.Journal of Computational Physics，2006，214（1）：137－170.

［11］HARRIS CD.Two－dimensional aerodynamic characteristics of the NACA0012airfoil in the Langley 8－foot transonic pressure tunnel［R］.NASA TM 81927，1981.

［12］COOK P H，DONALD M A，F(xiàn)IRMIN M C P.Aerofoil RAE2822pressure distribution，boundary layer and wake measurement［R］.AGARD AR 138，1979.

［13］DELANAYE M，et al.Automatic hybrid－cartesian grid generation for high－Reynolds number flows around complex geometries［R］.AIAA 1999－777，1999.

［14］MORYOSSEF Y，LEVY Y.Unconditionally positive implicit procedure for two－equation turbulence models：application to kωturbulence models［J］.Journal of Computational Physics，2006，220（1）：88－108.

［15］LIEN F S，KALITZIN G.Computations of transonic flow with thev2－fturbulence model［J］.International Journal for Heat and Fluid Flow，2001，22（1）：53－61.

［16］DRIVER D M，SEEGMILLER H L.Features of a reattaching turbulent shear layer in divergent channel flow［J］.AIAA Journal，1985，23（1）：163－171.