γ-Reθ轉捩模型的標定研究

張玉倫，王光學，孟德虹，王運濤

(中國空氣動力研究與發展中心空氣動力學國家重點實驗室，四川 綿陽 621000)

0 引言

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

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

經驗關聯方法就是把自由流的湍流度和當地的壓力梯度等通過試驗數據關聯到轉捩動量厚度雷諾數。典型的關聯如Abu-Ghannam和Shaw的關聯。然后，通過比較實際的和關聯的動量厚度雷諾數來確定轉捩的起始位置。由于要計算非局部變量-動量厚度，該方法難以與現代CFD方法相匹配，特別對于非結構網格和并行計算，問題更突出。

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

低雷諾數湍流模型適用于三維流動，也能夠與現代CFD方法相匹配，但是該方法不能模擬旁路轉捩，模擬自然轉捩的能力也受到質疑，這是因為，低雷諾數湍流模型中阻尼函數的標定依據是再現粘性底層的行為，而不是從層流到湍流的轉捩。

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

為了適應現代 CFD 計算，Menter＆Langtry［5］提出了基于當地關聯的γ-Reθ轉捩模型(local correlation-based transition modeling，LCTM)。該模型把經驗關聯方法和間歇因子方法有機地結合了起來，通過經驗關聯函數-轉捩動量厚度雷諾數來控制邊界層內間歇因子的生成，再通過間歇因子來控制湍流模型中湍流的生成。該模型通過把動量厚度雷諾數與當地的最大應變率(或渦量)相關聯，回避了動量厚度的計算，又通過Reθ輸運方程回避了經驗關聯函數的非當地計算，從而實現了模型計算的當地化。該模型并不試圖模擬邊界層內轉捩的物理過程，只是為把經驗關聯方法融入到現代CFD中去提供一個有效途徑。該模型的特點是:1)可以使用不同的經驗關聯函數進行標定;2)能夠涵蓋不同的轉捩機制;3)不依賴初場，即無論初始邊界層是層流還是湍流，都有相同的解;4)不影響基礎湍流模型在完全湍流區的行為;5)不依賴坐標系的選取;6)適用于三維邊界層流動。因此，該模型取得了很大的成功。但是，該方法要用到幾個經驗關聯函數，由于模型的標定對所使用的CFD軟件平臺的差異是敏感的，所以，這些關聯函數在不同的CFD軟件平臺之間缺乏互換性。

本文的主要工作就是補充了相應的關聯函數，基于平板的試驗數據在CARDC的trip軟件平臺上進行了標定，并研究了關聯函數的變化對轉捩結果的影響規律。

1 γ-Reθ轉捩模型簡介

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

上式中的S為應變率的模，Ω為渦量的模，y為離壁面的最小距離，Flength是轉捩區的長度，Reθc是邊界層內間歇函數開始增加位置的臨界動量厚度雷諾數。Flength和Reθc是經驗關聯函數，它們都是轉捩動量厚度雷諾數Reθt的函數。轉捩雷諾數Reθt是邊界層外自由流湍流度Tu、流向壓力梯度參數λθ等的函數。

對于關聯函數Reθt，不同文獻中的形式略有不同，本文采用Langtry［6］建議的形式:

其中U是當地的速度，θ是動量厚度，s為流線的弧長。F(λθ)代表了壓力梯度的影響。在上式中，Reθt是用當地的湍流度Tu和壓力梯度參數λθ計算的，直接用到邊界層內顯然不符合實際，為了解決這個問題，Menter等專門使用了一個輸運方程，即～Reθt方程。使得邊界層以外的 ～Reθt通過經驗關聯 Reθt獲得，邊界層之內的～Reθt通過輸運方程從邊界層以外的 ～Reθt擴散而來。～Reθt輸運方程［5］:

源項Pθt的設計目標就是使邊界層以外的等于Reθt，所以式子中有個(Reθt-)。混合函數 Fθt用來關掉邊界層內的生成項Pθt從而使得從邊界層以外擴散得到，所以設計Fθt的原則是在邊界層以外(即自由流)Fθt=0而在邊界層內 Fθt=1.0;Fwake的作用是保證混合函數Fθt在尾跡區不被激活。

為了模擬分離誘導轉捩，對間歇函數γ進行如下修正［5］:

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

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

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

表1 參數無量綱化表Table1 Parameter dimension

2 關聯函數的選擇和標定

標定工作在CARDC自行開發的trip軟件平臺上進行。對流項的離散選用MUSCL_roe格式，耗散項的離散采用二階中心格式，離散方程采用LU-SGS方法求解，為了適應低速情況和加速收斂，采用了預處理技術和多重網格技術。

標定依據文獻［5］提供的 T3A、T3B、T3A-和 Schubauer＆Klebanof等平板試驗數據，試驗條件見表2。

表2 平板轉捩試驗的進口條件Table2 Inlet condition for the flat plate test cases

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

圖1 平板的計算網格Fig.1 computational mesh for flat plate

該轉捩模型經由關聯函數Reθt強烈依賴當地自由流湍流度，只有正確模擬了自由流湍流度，標定工作才有意義。由于在進口處指定的湍流度會迅速衰減，而且進口渦粘性比(μt/μ)越小，湍流度衰減越快，可是，如果進口渦粘性比給的太大，壁面摩阻又會偏離層流太多，所以應該認真指定進口湍流度和進口渦粘性比，以保證平板前緣自由流湍流度有正確的值和期望的衰減率。在自由流中，湍動能的衰減可以根據下式［7］計算

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

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

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

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

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

通過分析不同文獻的計算結果，我們發現:為了盡可能模擬轉捩區長度，相對于，關聯函數 Reθc在低區要盡可能大(≈1)，而在高 ～Reθt區則要適當的小。因此，本文選擇如下形式的:

其中Cθc是一個常數，用來控制 Reθc的大小。為了考察不同 Reθc的影響，本文分別考察了 Cθc=0.15、0.24、0.33三種情況，對應的 Reθc～如圖3所示。標定后的關聯函數Flength擬合曲線如圖4所示。如圖5給出了在標定點上計算的壁面摩擦阻力系數與試驗的比較。

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

圖4 F length～ Fig.4 F length ～

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

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

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

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

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

在標定過程中發現，轉捩位置以后的壁面摩阻偏小，說明轉捩后邊界層內的湍動能耗散過大，這主要是對湍動能方程進行修正時使用γD1min造成的。可是，如果取消γD1min，又會造成在高自由流湍流度的情況下轉捩前的壁面摩阻偏離層流太多。為此，本文把上述常量γDkmin修正為可變的，以消除其對轉捩后壁面摩阻的影響，同時又能保證轉捩前壁面摩阻的性質。

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

圖7 間歇函數γ云圖Fig.7 Contours of Intermittencyγ

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

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

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

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

圖10 F length～ Fig.10 F length ～

3 結論

關聯函數Flength和Reθc不是相互獨立的，形式也不是固定的。對于給定的 Reθc，要保持轉捩位置不變，Reθc和Flength要同時減小或增大。Reθc和Flength減小，轉捩區會變寬，轉捩位置后的壁面摩阻會減小。

Reθc和 Flength過大會導致收斂性變差。

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

［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.