張玉倫，王光學(xué)，孟德虹，王運(yùn)濤

(中國空氣動(dòng)力研究與發(fā)展中心空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，四川 綿陽 621000)

0 引言

隨著飛行器性能越來越高，工程上對氣動(dòng)力數(shù)據(jù)的精度要求也越來越高。由于邊界層轉(zhuǎn)捩對摩擦阻力、流動(dòng)的分離位置等有很大的影響，因此對邊界層轉(zhuǎn)捩位置的準(zhǔn)確模擬就成為一個(gè)非常重要的工程問題。特別在中等雷諾數(shù)范圍，層流區(qū)域和湍流區(qū)域具有相同的量級時(shí)，使用全層流或全湍流計(jì)算方法都會(huì)導(dǎo)致很大的計(jì)算誤差，對轉(zhuǎn)捩位置的模擬就更顯得必要。

轉(zhuǎn)捩的機(jī)制非常復(fù)雜，比如有自然轉(zhuǎn)捩(由T-S波或橫向流失穩(wěn)造成)、旁路轉(zhuǎn)捩(自由流中的高水平湍流度/擾動(dòng)施加在層流邊界層上造成)、分離誘導(dǎo)轉(zhuǎn)捩以及湍流在順壓梯度區(qū)的再層流化等，使得轉(zhuǎn)捩的模擬非常困難。雖然NS方程具有模擬轉(zhuǎn)捩過程的能力，而且隨著計(jì)算機(jī)的發(fā)展，近年來湍流的高級數(shù)值模擬方法，如直接數(shù)值模擬(DNS)和大渦模擬(LES)得到了很大發(fā)展，也取得了許多令人滿意的結(jié)果，但是距離工程實(shí)用還有很長的距離。在工程應(yīng)用中，轉(zhuǎn)捩主要靠經(jīng)驗(yàn)或半經(jīng)驗(yàn)的方法來確定，比如經(jīng)驗(yàn)關(guān)聯(lián)方法、eN方法和低雷諾數(shù)湍流模型方法等。

經(jīng)驗(yàn)關(guān)聯(lián)方法就是把自由流的湍流度和當(dāng)?shù)氐膲毫μ荻鹊韧ㄟ^試驗(yàn)數(shù)據(jù)關(guān)聯(lián)到轉(zhuǎn)捩動(dòng)量厚度雷諾數(shù)。典型的關(guān)聯(lián)如Abu-Ghannam和Shaw的關(guān)聯(lián)。然后，通過比較實(shí)際的和關(guān)聯(lián)的動(dòng)量厚度雷諾數(shù)來確定轉(zhuǎn)捩的起始位置。由于要計(jì)算非局部變量-動(dòng)量厚度，該方法難以與現(xiàn)代CFD方法相匹配，特別對于非結(jié)構(gòu)網(wǎng)格和并行計(jì)算，問題更突出。

eN方法亦稱為線性穩(wěn)定性分析方法，是基于線性穩(wěn)定性理論和試驗(yàn)數(shù)據(jù)的半經(jīng)驗(yàn)方法。該方法假設(shè):轉(zhuǎn)捩過程是由層流邊界層內(nèi)最初的小擾動(dòng)向下游傳播放大達(dá)到eN時(shí)完成。此法歸結(jié)為確定放大因子N，所以N值就成了判斷轉(zhuǎn)捩的準(zhǔn)則。由于eN方法本身沒有解決如何選取N值，也沒有考慮環(huán)境條件對轉(zhuǎn)捩的影響，因而在實(shí)際應(yīng)用中是靠低湍流度試驗(yàn)來確定N的。但是eN方法難以處理旁路轉(zhuǎn)捩和推廣到三維情況。

低雷諾數(shù)湍流模型適用于三維流動(dòng)，也能夠與現(xiàn)代CFD方法相匹配，但是該方法不能模擬旁路轉(zhuǎn)捩，模擬自然轉(zhuǎn)捩的能力也受到質(zhì)疑，這是因?yàn)?，低雷諾數(shù)湍流模型中阻尼函數(shù)的標(biāo)定依據(jù)是再現(xiàn)粘性底層的行為，而不是從層流到湍流的轉(zhuǎn)捩。

還有一類基于間歇因子γ的轉(zhuǎn)捩預(yù)測方法。間歇因子被定義為空間某點(diǎn)的流態(tài)是湍流的概率，是Dhawan＆Narasimha［1］首先引入的。之后出現(xiàn)了大量的間歇因子的轉(zhuǎn)捩預(yù)測方法，例如:Cho＆Chung［2］針對自由剪切流發(fā)展了與k-ε湍流模型聯(lián)合使用的間歇因子輸運(yùn)方程方法，Steelant＆Dick［3］發(fā)展了與條件平均 NS方程聯(lián)合使用的間歇因子輸運(yùn)方程方法，Suzen＆Huang［4］將 Steelant＆ Dick 的模型與 Cho ＆ Chung 的模型相結(jié)合發(fā)展了間歇因子的對流-擴(kuò)散方程。但是所有這些轉(zhuǎn)捩模型都需要積分動(dòng)量厚度，難以與現(xiàn)代CFD方法相匹配。

為了適應(yīng)現(xiàn)代 CFD 計(jì)算，Menter＆Langtry［5］提出了基于當(dāng)?shù)仃P(guān)聯(lián)的γ-Reθ轉(zhuǎn)捩模型(local correlation-based transition modeling，LCTM)。該模型把經(jīng)驗(yàn)關(guān)聯(lián)方法和間歇因子方法有機(jī)地結(jié)合了起來，通過經(jīng)驗(yàn)關(guān)聯(lián)函數(shù)-轉(zhuǎn)捩動(dòng)量厚度雷諾數(shù)來控制邊界層內(nèi)間歇因子的生成，再通過間歇因子來控制湍流模型中湍流的生成。該模型通過把動(dòng)量厚度雷諾數(shù)與當(dāng)?shù)氐淖畲髴?yīng)變率(或渦量)相關(guān)聯(lián)，回避了動(dòng)量厚度的計(jì)算，又通過Reθ輸運(yùn)方程回避了經(jīng)驗(yàn)關(guān)聯(lián)函數(shù)的非當(dāng)?shù)赜?jì)算，從而實(shí)現(xiàn)了模型計(jì)算的當(dāng)?shù)鼗Ｔ撃Ｐ筒⒉辉噲D模擬邊界層內(nèi)轉(zhuǎn)捩的物理過程，只是為把經(jīng)驗(yàn)關(guān)聯(lián)方法融入到現(xiàn)代CFD中去提供一個(gè)有效途徑。該模型的特點(diǎn)是:1)可以使用不同的經(jīng)驗(yàn)關(guān)聯(lián)函數(shù)進(jìn)行標(biāo)定;2)能夠涵蓋不同的轉(zhuǎn)捩機(jī)制;3)不依賴初場，即無論初始邊界層是層流還是湍流，都有相同的解;4)不影響基礎(chǔ)湍流模型在完全湍流區(qū)的行為;5)不依賴坐標(biāo)系的選取;6)適用于三維邊界層流動(dòng)。因此，該模型取得了很大的成功。但是，該方法要用到幾個(gè)經(jīng)驗(yàn)關(guān)聯(lián)函數(shù)，由于模型的標(biāo)定對所使用的CFD軟件平臺(tái)的差異是敏感的，所以，這些關(guān)聯(lián)函數(shù)在不同的CFD軟件平臺(tái)之間缺乏互換性。

本文的主要工作就是補(bǔ)充了相應(yīng)的關(guān)聯(lián)函數(shù)，基于平板的試驗(yàn)數(shù)據(jù)在CARDC的trip軟件平臺(tái)上進(jìn)行了標(biāo)定，并研究了關(guān)聯(lián)函數(shù)的變化對轉(zhuǎn)捩結(jié)果的影響規(guī)律。

1 γ-Reθ轉(zhuǎn)捩模型簡介

兩方程轉(zhuǎn)捩模型的最終目的是求解間歇函數(shù) (即空間某點(diǎn)的流態(tài)是湍流的概率，0≤γ≤1)，并通過它與湍流模型的聯(lián)合來控制轉(zhuǎn)捩的發(fā)生。無量綱的輸運(yùn)方程［5］的守恒形式為:

上式中的S為應(yīng)變率的模，Ω為渦量的模，y為離壁面的最小距離，F(xiàn)length是轉(zhuǎn)捩區(qū)的長度，Reθc是邊界層內(nèi)間歇函數(shù)開始增加位置的臨界動(dòng)量厚度雷諾數(shù)。Flength和Reθc是經(jīng)驗(yàn)關(guān)聯(lián)函數(shù)，它們都是轉(zhuǎn)捩動(dòng)量厚度雷諾數(shù)Reθt的函數(shù)。轉(zhuǎn)捩雷諾數(shù)Reθt是邊界層外自由流湍流度Tu、流向壓力梯度參數(shù)λθ等的函數(shù)。

對于關(guān)聯(lián)函數(shù)Reθt，不同文獻(xiàn)中的形式略有不同，本文采用Langtry［6］建議的形式:

其中U是當(dāng)?shù)氐乃俣?，θ是?dòng)量厚度，s為流線的弧長。F(λθ)代表了壓力梯度的影響。在上式中，Reθt是用當(dāng)?shù)氐耐牧鞫萒u和壓力梯度參數(shù)λθ計(jì)算的，直接用到邊界層內(nèi)顯然不符合實(shí)際，為了解決這個(gè)問題，Menter等專門使用了一個(gè)輸運(yùn)方程，即～Reθt方程。使得邊界層以外的 ～Reθt通過經(jīng)驗(yàn)關(guān)聯(lián) Reθt獲得，邊界層之內(nèi)的～Reθt通過輸運(yùn)方程從邊界層以外的 ～Reθt擴(kuò)散而來?！玆eθt輸運(yùn)方程［5］:

源項(xiàng)Pθt的設(shè)計(jì)目標(biāo)就是使邊界層以外的等于Reθt，所以式子中有個(gè)(Reθt-)?；旌虾瘮?shù) Fθt用來關(guān)掉邊界層內(nèi)的生成項(xiàng)Pθt從而使得從邊界層以外擴(kuò)散得到，所以設(shè)計(jì)Fθt的原則是在邊界層以外(即自由流)Fθt=0而在邊界層內(nèi) Fθt=1.0;Fwake的作用是保證混合函數(shù)Fθt在尾跡區(qū)不被激活。

為了模擬分離誘導(dǎo)轉(zhuǎn)捩，對間歇函數(shù)γ進(jìn)行如下修正［5］:

該轉(zhuǎn)捩模型只是獲得間歇函數(shù)γ，還需要與湍流模型聯(lián)合才能模擬轉(zhuǎn)捩過程。Menter采用的湍流模型為SST模型。具體的聯(lián)合方式就是使用間歇函數(shù)γ來修正k方程的生成項(xiàng)、破壞項(xiàng)和混合函數(shù)［5］，即:

其中:Pk，Dk和Florig是SST模型中原來的生成項(xiàng)、破壞項(xiàng)和混合函數(shù)，Menter推薦的γDkmin=0.1;對ω方程不做修正。

上述方程使用的無量綱化參數(shù)見表1。

表1 參數(shù)無量綱化表Table1 Parameter dimension

2 關(guān)聯(lián)函數(shù)的選擇和標(biāo)定

標(biāo)定工作在CARDC自行開發(fā)的trip軟件平臺(tái)上進(jìn)行。對流項(xiàng)的離散選用MUSCL_roe格式，耗散項(xiàng)的離散采用二階中心格式，離散方程采用LU-SGS方法求解，為了適應(yīng)低速情況和加速收斂，采用了預(yù)處理技術(shù)和多重網(wǎng)格技術(shù)。

標(biāo)定依據(jù)文獻(xiàn)［5］提供的 T3A、T3B、T3A-和 Schubauer＆Klebanof等平板試驗(yàn)數(shù)據(jù)，試驗(yàn)條件見表2。

表2 平板轉(zhuǎn)捩試驗(yàn)的進(jìn)口條件Table2 Inlet condition for the flat plate test cases

計(jì)算網(wǎng)格如圖1，平板長度為2.0，來流遠(yuǎn)場邊界距平板前緣0.22，法向遠(yuǎn)場邊界距壁面也是0.22，網(wǎng)格為H型，壁面上第一層網(wǎng)格的法向距離1.0×10-6，平板前緣第一網(wǎng)格的x向距離1.5×10-3。網(wǎng)格規(guī)模為425×129，其中壁面上有392網(wǎng)格點(diǎn)。

圖1 平板的計(jì)算網(wǎng)格Fig.1 computational mesh for flat plate

該轉(zhuǎn)捩模型經(jīng)由關(guān)聯(lián)函數(shù)Reθt強(qiáng)烈依賴當(dāng)?shù)刈杂闪魍牧鞫?，只有正確模擬了自由流湍流度，標(biāo)定工作才有意義。由于在進(jìn)口處指定的湍流度會(huì)迅速衰減，而且進(jìn)口渦粘性比(μt/μ)越小，湍流度衰減越快，可是，如果進(jìn)口渦粘性比給的太大，壁面摩阻又會(huì)偏離層流太多，所以應(yīng)該認(rèn)真指定進(jìn)口湍流度和進(jìn)口渦粘性比，以保證平板前緣自由流湍流度有正確的值和期望的衰減率。在自由流中，湍動(dòng)能的衰減可以根據(jù)下式［7］計(jì)算

其中 β=0.09，β*=0.0828，t=x/V為時(shí)間，x是下游距進(jìn)口的流向距離，V是平均速度。使用湍流度和渦粘性比(μ1/μ)，可以將上述湍動(dòng)能的衰減規(guī)律重新寫成:

圖2給出了各種情況下邊界層外自由流湍流度的衰減與試驗(yàn)的比較，計(jì)算與試驗(yàn)吻合。

關(guān)聯(lián)函數(shù)Flength和Reθc并不是相互獨(dú)立的，形式也不是固定的。為了便于確定他們，本文采取先給定其中的一個(gè)(Reθc)，然后通過平板的試驗(yàn)數(shù)據(jù)來確定另一個(gè)(Flength)。同時(shí)，為了能夠在校驗(yàn)中有靈活而足夠的自由度，F(xiàn)length的形式采用樣條擬合曲線。

圖2 自由流湍流度衰減與試驗(yàn)的比較Fig.2 Comparison of free stream turbulence decay between CFD calculations and wind tunnel data

盡管該模型對轉(zhuǎn)捩過程有一定的模擬能力，但畢竟轉(zhuǎn)捩過程是復(fù)雜的，而且對于大部分雷諾數(shù)比較高的航空問題，能夠正確模擬轉(zhuǎn)捩位置就基本能夠滿足工程需求了，所以本文中關(guān)聯(lián)函數(shù)的標(biāo)定主要依據(jù)轉(zhuǎn)捩位置，適當(dāng)兼顧轉(zhuǎn)捩區(qū)長度。同時(shí)，雖然Flength被稱之為轉(zhuǎn)捩區(qū)的長度，實(shí)際上這種叫法并不準(zhǔn)確。比如，如果固定Reθt和 Reθc，那么增加 Flength，轉(zhuǎn)捩區(qū)確實(shí)會(huì)加長，但此時(shí)轉(zhuǎn)捩位置也會(huì)提前;如果要保持轉(zhuǎn)捩位置不變，那么增加 Flength，就要相應(yīng)增加 Reθc(固定 Reθt)，此時(shí)轉(zhuǎn)捩區(qū)反而會(huì)變短。也就是說，在保持轉(zhuǎn)捩位置不變的前提下，改變 Reθc就會(huì)改變轉(zhuǎn)捩區(qū)長度。因此，在本文中 Reθc的選擇主要依據(jù)轉(zhuǎn)捩區(qū)長度，相反，F(xiàn)length的確定則只依據(jù)轉(zhuǎn)捩位置。

通過分析不同文獻(xiàn)的計(jì)算結(jié)果，我們發(fā)現(xiàn):為了盡可能模擬轉(zhuǎn)捩區(qū)長度，相對于，關(guān)聯(lián)函數(shù) Reθc在低區(qū)要盡可能大(≈1)，而在高 ～Reθt區(qū)則要適當(dāng)?shù)男?。因此，本文選擇如下形式的:

其中Cθc是一個(gè)常數(shù)，用來控制 Reθc的大小。為了考察不同 Reθc的影響，本文分別考察了 Cθc=0.15、0.24、0.33三種情況，對應(yīng)的 Reθc～如圖3所示。標(biāo)定后的關(guān)聯(lián)函數(shù)Flength擬合曲線如圖4所示。如圖5給出了在標(biāo)定點(diǎn)上計(jì)算的壁面摩擦阻力系數(shù)與試驗(yàn)的比較。

圖 3 Reθc ～Fig.3 Reθc ～

圖4 F length～ Fig.4 F length ～

圖5 壁面摩擦阻力系數(shù)與試驗(yàn)的比較Fig.5 Comparisons of skin friction coefficient distribution between CFD calculations and wind tunnel tests

從圖3～圖6示結(jié)果可以看出:對于給定的Reθt，要保持轉(zhuǎn)捩位置不變，Reθc減小，F(xiàn)length也要相應(yīng)減小(見圖3 和圖 4)。同時(shí)，增大 Cθc即減小 Reθc和 Flength，轉(zhuǎn)捩區(qū)會(huì)變寬，轉(zhuǎn)捩區(qū)后的壁面摩阻會(huì)減小(見圖5)，好在對摩阻的這種影響在一段距離后會(huì)消失。另外，Reθc和Flength過大(即Cθc過小)會(huì)導(dǎo)致計(jì)算收斂性變差，如圖6所示。

還可以看出:對于來流湍流度比較低的Schubauer＆ Klebanof和T3A-兩種情況(對應(yīng)于較高的Reθt)，計(jì)算的轉(zhuǎn)捩區(qū)比試驗(yàn)的窄，而轉(zhuǎn)捩區(qū)后的壁面摩阻比試驗(yàn)的小，要增加轉(zhuǎn)捩區(qū)寬度，就要減小Reθc和Flength，這樣就會(huì)導(dǎo)致轉(zhuǎn)捩區(qū)后的壁面摩阻更小，因此對于來流湍流度比較低的情況，Cθc既不能太大也不能太小。對于來流湍流度較高的T3A和T3B兩種情況(對應(yīng)于較低的Reθt)，計(jì)算的轉(zhuǎn)捩區(qū)比試驗(yàn)的還寬些，但是，此時(shí)(～Reθt＜200)，Reθc幾乎已經(jīng)達(dá)到最大(Reθc≈0.99～Reθt)，不可能通過進(jìn)一步增大Reθc來減小轉(zhuǎn)捩區(qū)寬度。同時(shí)對于來流湍流度較高的情況，Cθc的變化對 Reθc幾乎沒有影響。下面的討論就只是針對Cθc=0.24的情況，相應(yīng)Flength樣條曲線的控制點(diǎn)列于表3。

圖6 不同的Reθc收斂情況比較(T3A-FSTI=0.87)Fig.6 Comparisons of convergencefor for different Reθc

圖7給出了間歇函數(shù)γ分布云圖，從圖中可以清楚看出間歇函數(shù)在層流邊界層中增長直到轉(zhuǎn)捩的情況。圖8給出了在轉(zhuǎn)捩前的不同x位置轉(zhuǎn)捩尺度～Reθt在邊界層附近沿法向的分布，可以看出自由流中的Reθt向邊界層內(nèi)部擴(kuò)散影響的情況，邊界層內(nèi)的～Reθt與自由流中的Reθt存在某種滯后。

在標(biāo)定過程中發(fā)現(xiàn)，轉(zhuǎn)捩位置以后的壁面摩阻偏小，說明轉(zhuǎn)捩后邊界層內(nèi)的湍動(dòng)能耗散過大，這主要是對湍動(dòng)能方程進(jìn)行修正時(shí)使用γD1min造成的?？墒?，如果取消γD1min，又會(huì)造成在高自由流湍流度的情況下轉(zhuǎn)捩前的壁面摩阻偏離層流太多。為此，本文把上述常量γDkmin修正為可變的，以消除其對轉(zhuǎn)捩后壁面摩阻的影響，同時(shí)又能保證轉(zhuǎn)捩前壁面摩阻的性質(zhì)。

表3 樣條曲線的控制點(diǎn)列表(Cθc=0.24)Table3 Control point of spline curve(Cθc=0.24)

圖7 間歇函數(shù)γ云圖Fig.7 Contours of Intermittencyγ

圖8 轉(zhuǎn)捩前邊界層附近轉(zhuǎn)捩尺度～Reθt的分布Fig.8 Profiles of the～Reθt for a flat plate

圖9 γDkmin有∕無修正時(shí)壁面摩阻的比較Fig.9 Comparisons of skin friction coefficient distribution between modified and nomodifiedγDkmin

圖9給出了修正前后壁面摩阻的比較，表明修正取得了滿意的效果。

當(dāng)然，使用四個(gè)試驗(yàn)點(diǎn)來標(biāo)定關(guān)聯(lián)函數(shù)實(shí)在太過勉強(qiáng)。比如，圖10給出的兩條Flength～曲線在標(biāo)定點(diǎn)上是重合的，數(shù)值試驗(yàn)也表明這兩個(gè)關(guān)聯(lián)函數(shù)Flength在四個(gè)標(biāo)定點(diǎn)上給出了相同的轉(zhuǎn)捩位置，但是，偏離標(biāo)定點(diǎn)之后，二者給出的轉(zhuǎn)捩位置就有比較大的差異。所以，要使該模型達(dá)到實(shí)用的程度，還需要更多的試驗(yàn)點(diǎn)來完善關(guān)聯(lián)函數(shù)Flength。這也是關(guān)聯(lián)函數(shù)Flength采用樣條形式的原因。

圖10 F length～ Fig.10 F length ～

3 結(jié)論

關(guān)聯(lián)函數(shù)Flength和Reθc不是相互獨(dú)立的，形式也不是固定的。對于給定的 Reθc，要保持轉(zhuǎn)捩位置不變，Reθc和Flength要同時(shí)減小或增大。Reθc和Flength減小，轉(zhuǎn)捩區(qū)會(huì)變寬，轉(zhuǎn)捩位置后的壁面摩阻會(huì)減小。

Reθc和 Flength過大會(huì)導(dǎo)致收斂性變差。

本文對γDkmin的修正有效地改進(jìn)了轉(zhuǎn)捩位置后壁面摩阻偏小的狀況。

［1］DHAWAN S，NARASIMHA R.Some properties of boundary layer flow during transition from laminar to turbulent motion［J］.Journal of Fluid Mechanics，1958，3(4):418-436.

［2］CHO J R，CHUNG M K.A equation turbulence model［J］.Journal of Fluid Mechanics，1992，237:301-322.

［3］STEELANT J，DICK E.Modelling of bypass transition with conditioned Navier-Stokes equations coupled to an intermittency transport equation［J］.International Journal for Numerical Methods in Fluids，1996，23:193-220.

［4］SUZEN Y，HUANG P.Modeling of flow transition using an intermittency transport equation［J］.Journal of Fluids and Engineering，2000，122(2):273-284.

［5］MENTER F R，LANGTRY R B，LIKKI S R，SUZEN Y B，HUANG P G，V?LKER S.A correlation based transition model using local variables:part I-model formulation［R］.ASME-GT 2004-53452.Proceedings of ASME Turbo Expo 2004.Vienna，Austria，2004.

［6］LANGTRY R B.A correlation-based transition model using local variables for unstructured parallelized cfd codes［D］.［Ph.D.Thesis］.Stuttgart，2006.

［7］NIELS N S?RENSEN.CFD modeling of laminar-turbulent transition for airfoils and rotors using the Model［R］.AIAA Paper 2008-7323.

［8］LANGTRY R B，MENTER F R.Transition modeling for general cfd applications in aeronautics［R］.AIAA Paper 2005-0522.

［9］MISAKA T，OBAYASHI S.Application of correlation-based transition models to flows around wings［R］.AIAA Paper 2006-918.

［10］PETTERSSON K，CRIPPA S.Implementation and verification of a correlation based transition prediction method［R］.AIAA Paper 2008-4401.