基于轉捩模型的位移邊界層研究

鄧家鈺，王成恩

(上海交通大學機械與動力工程學院，上海 201100)

1 引言

在航空航天領域中，邊界層的層流到湍流的轉捩過程的研究是計算流體力學(CFD)領域的一個關鍵且困難的問題。工程上對轉捩流動的模擬大多是采用由實驗總結出來的經驗關系式，國外研究人員提出的一些轉捩經典模型，如Abu-Ghanam&Ahaw(AGS)模型、Steelant & Dik模型、Suzen & Huang(S-H)模型、Menter & Langry(M-L)[1-5]模型等。其中M-L轉捩模型對之前提出的轉捩模型進行修正，在基于SST k-ω[6]湍流模型的基礎上，融入基于流場當地變量的轉捩模型。該模型結合了轉捩經驗關系式和低雷諾數湍流模型的優勢，提高了模型的通用性和準確性。

轉捩現象與邊界層外的流動密切相關，對于平板邊界層來說，當自由流湍流強度低于1%發生自然轉捩，即邊界層失去穩定導致層流瓦解；當湍流強度大于1%的，邊界層對于自由流脈動的感受而導致bypass轉捩[7]。Jacobs&Du-bin[8]以及Brandt[9]、Schaltter[10]關于自由流湍流影響下的bypass轉捩的數值分析研究做了詳細的工作。

關于平板邊界層的研究，G.B.Schubauer[11]和A.M.Sacill[12]進行了平板邊界層轉捩實驗。國內學者董平[13]對經典的T3系列[12]平板邊界層實驗進行實驗和仿真對比，顧金生[14]提出了不依賴經驗系數的湍流模型，郭玉波等人[15]進行了粘性流場的數值模擬實驗。

在平板邊界層轉捩的研究中，何克敏[16]所研究的低湍流度風洞中湍流度對平板邊界層轉捩影響的實驗、徐國亮[17]和符松[18]等人對平板的bypass轉捩進行的實驗研究。J.H.M.Fransson等人[19]研究的受自由流湍流影響的平板轉捩實驗、Schrader[20]對自由流為湍流的橫掠平板轉捩研究的數值仿真、陳奕[21]對平板邊界層轉捩的模擬仿真，引發了本文的研究點：對于來流為自由流湍流的平板的位移邊界層厚度隨著自由流湍流強度改變有著什么樣的變化規律。

本文利用商用CFD軟件FLUENT的SST tran-sition模型，模擬了經典的平板邊界層轉捩實驗T3A和T3B[13]，以及Schrader[20]所做的平板邊界層數值模擬實驗。分析模擬結果得出以下結論：在自由流湍流強度大于1%的前提下，隨著湍流強度的減小，位移邊界層厚度沿著流向方向的變化曲線在轉捩區域出現下凹現象，且隨著湍流強度減少，現象越來越明顯；在自由流湍流強度小于1%的提前下，在湍流強度接近1%，位移邊界層厚度沿著流向方向的變化曲線在轉捩區域出現下凹現象，但隨著湍流強度減少，該現象逐漸消失。

2 理論方法

SST transition模型是SST k-ω傳輸方程和其它兩個傳輸方程耦合的結果，其中一個運輸方程用于控制間歇，另外一個方程為控制轉捩起始條件，用動量厚度雷諾數表示。

2.1 流體一般控制方程

質量守恒方程，即連續性方程

(1)

動量守恒方程，即Navier-Stokes方程

(2)

能量守恒方程

(3)

(4)

其中，Φ是粘性導致的能量耗散函數。

2.2 Transition SST模型的傳輸方程

傳輸方程的間歇系數γ定義如下

(5)

傳輸源項定義為：

PJ1=2FlengthρS[γFlength]CJ3

(6)

(7)

其中，S為應變速率大小，Flength是一個控制轉捩區域長度經驗相關的因子。消減項/再層流化項被定義為

PJ2=(2CJ1)ρΩγFturb

(8)

Er3=Cr3

(9)

其中，Ω是渦量大小。轉捩的開始由以下的系數控制

(10)

(11)

(12)

(13)

(14)

Fonset=max(Fonset2-Fonset3，0)

(15)

(16)

CJ1=0.03；CJ1=50；CJ3=0.5；σJ=1.0

2.3 轉捩修正

對分離引起的轉捩的修正如下

(17)

(18)

γeff=max(γ，γsep)

(19)

(20)

源項定義為

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

Reθt=f(Tu，λ)

(29)

(30)

(31)

第一個經驗相關式子為當地湍流強度方程，Thawaite壓力梯度因子λθ定義為

λθ=(θ2/v)dU/dS

(32)

2.4 耦合轉捩模型和SST傳輸方程

給出SST傳輸方程后，為了計算更為準確，需要將SST傳輸方程和轉捩模型耦合

(33)

Pk=γeffPkP

(34)

(35)

(36)

(37)

Ft=max(Florig，F3)

(38)

其中，Dk和Pk為SST模型的原始生成項和原始銷毀項，F1orig為原始SST混合方程。

2.5 指定入口湍流等級

在入口給定的湍流強度根據粘性比不斷衰減。粘性比越大，則湍流衰減率越小。湍動能的衰減由以下公式計算

(39)

對于在自由流中SST湍流模型的常數為

β=0.09，β*=0.0828

時間尺度由下面定義

(40)

x為進口流向距離，V為平均對流速度。渦流粘度被定義

(41)

湍動能衰減方程可由入口湍流強度和渦流粘性比重新定義為

(42)

2.6 位移邊界層的定義

位移邊界層厚度δ1的計算公式為

(43)

其中，u為流體速度，u∞為無窮遠處自由流速度。

3 零壓力梯度平板邊界層仿真

使用商業CFD軟件FLUENT的不可壓縮非結構隱式RANS求解器對流動進行仿真，基本原理是數值求解流體力學基本方程Navier-Stokes(N-S)方程和SST k-ω湍流模型的控制方程、γ-Reθt轉捩模型的控制方程。通過對經典的平板邊界層轉捩實驗T3A和T3B[13]仿真，以及對Schrader[20]的數值仿真進一步深入探究，簡化其幾何模型，設置不同的入口湍流強度進行計算，分析零壓力梯度平板邊界層流動的位移邊界層變化規律。

3.1 T3A和T3B仿真

通過對董平[13]T3A和T3B的實驗，該實驗的示意圖如圖1所示，進口自由流的湍流強度、進口自由流的速度以及粘性比如表1所示，將其看作二維流動。幾何模型大小為1700mm×200mm，網格劃分為1700×600。為了更為準確的預測轉捩過程以及得到邊界層內的流動速度分布，故對轉捩發生的區域和壁面網格都進行了加密，距壁面的第一個網格點y+≤1。選擇的物理模型為SST transition四方程轉捩模型，差分格式均為二階迎風格式，所有物理參考量均為入口參數。入口參數設置與實驗進口條件相同，底面為No-Slip條件，頂部為Symmetry條件，出口設置為Outflow條件。

圖1 平板轉捩實驗示意圖

表1 實驗的進口參數

3.2 仿真結果

對于T3A和T3B實驗的計算結果，首先分析其壁面摩擦系數(Cf)能夠很好的表征轉捩發生的位置，通常認為Cf最小處為轉捩的起始點，最大處為轉捩的終止點)，與實驗所測的值進行比較如圖2和圖3所示。對于自由流湍流強度為3%的T3A的計算結果，計算得到壁面摩擦系數與實驗測量得到的值幾乎相同；對于湍流強度為6%的T3B的計算結果，計算得到的壁面摩擦系數與實驗值也基本一致。故可以分析得出使用SST transition模型模擬自由流橫掠平板的邊界層實驗較為準確，能夠比較準確捕捉轉捩點的起始位置和終止位置。

圖2 T3A 的摩擦系數

圖3 T3B的摩擦系數

圖4記錄了T3A和T3B位移邊界層厚度沿著流向方向的變化。仔細觀察和分析圖4，發現位移邊界層沿著流向方向變化有一個奇怪的現象：對于T3A的仿真結果，位移邊界層厚度沿著流向方向的變化趨勢是先陡峭增大而又有略微趨近平緩增長再陡峭的增大；而對于T3B的仿真結果，位移邊界層厚度沿著流向方向的變化趨勢是一直增大。而且對于T3A來說，從它的壁面摩擦趨勢圖顯示，x=0.5m附近為轉捩的起始位置(壁面摩擦系數最小)，x=0.8m附近(壁面摩擦系數最大)為轉捩的終止位置，恰巧位移邊界層趨于平緩增長的趨勢的位置大致在0.5m到0.8m這個區間。為了更加深入了解自由流橫掠平板的位移邊界層的變化規律，進一步對Schrader[20]所做的自由流湍流橫掠平板的轉捩數值模擬實驗(該數值實驗中，出現了位移邊界層沿著流向方向出現下凹趨勢)進行展開研究。

圖4 T3A和T3B位移邊界層厚度

3.3 Schrader數值仿真

將其實驗簡化為二維幾何模型，模型長度為1200mm，寬度為80mm，選用網格1700×600。為了準確預測轉捩過程和捕捉邊界層內的流動速度分布，故對壁面網格進行局部加密，距壁面的第一個網格點y+≤1。選擇物理模型為SST transition四方程轉捩模型，差分格式均為二階迎風格式，所有物理參考量均為入口參數。底面為No-Slip條件，頂部為Symmetry條件，出口設置為Outflow條件。入口條件設置不同的入口自由流湍流強度，湍流尺度為10mm。

3.3.1 仿真結果

自由流湍流強度大于1%的時候發生bypass轉捩，自由流湍流強度小于1%的時候發生自然轉捩。分別對自由流湍流強度大于1%和湍流強度小于1%計算和討論，分析計算結果和進行比較。

對于自由流的湍流強度Tu=5.06%、Tu=3.73%、Tu=2.95%、Tu=2.53%、Tu=2.11%，為了能夠較為清晰地捕捉轉捩現象的發生和保證仿真的轉捩發生在平板上，故選擇自由流入口速度為V=5m/s。

壁面摩擦系數沿著流向方向的計算結果如圖5所示，發現自由流的湍流強度越大，轉捩位置越靠前。計算其位移邊界層計算結果如圖6所示，發現轉捩位置越靠前，位移邊界層有下凹趨勢或者趨于平緩的位置越靠前，且該位置出現在轉捩區域附近。

圖5 不同湍流強度下的壁面摩擦系數

對于不同的自由流湍流強度，如果入口的自由流速度太小，則在有限長度的平板上可能不發生轉捩；如果入口的自由流速度太大，則在平板 前緣部分可能很早就發生轉捩，且不容易觀測到轉捩。對于自由流湍流強度Tu=2.11%、Tu=1.69%、Tu=1.26%設置來流入口速度V=10m/s，計算得到的壁面摩擦系數沿流向方向的變化如圖5所示，位移邊界層厚度沿流向方向的變化如圖6所示。

圖6 不同湍流強度下的位移邊界層厚度

觀察和分析自由流的湍流強度為Tu=2.11%、Tu=1.69%、Tu=1.26%的計算結果，發現自由流的湍流強度越大，轉捩位置越靠前。隨著自由流的湍流強度減少，相同入口速度下的位移邊界層沿著流向方向的下凹現象越來越明顯。

圖7 不同湍流強度下的壁面摩擦系數

圖8 不同湍流強度下的位移邊界層厚度

綜合自由流的湍流強度大于1%的計算結果，相對比較位移邊界層厚度沿著流向方向的變化趨勢發現，隨著入口自由流的湍流強度減少，位移邊界層沿著流向方向的變化趨勢從沿著流向方向一直陡峭地增大過渡到先陡峭地增大而后平緩地增大再到相對陡峭地增大，最后過渡到出現先陡峭地增大而后減少再增大的趨勢(位移邊界層厚度沿著流向方向的變化曲線出現下凹的現象)。

而對于自由流的湍流強度低于1%的設置，需要進行略微調整。因為在大量數值計算后發現，對于給定低于1%的自由流的湍流強度，橫掠平板的位移邊界層厚度對自由流速度非常敏感，難以捕捉位移邊界層趨于平緩增長或出現下降的變化趨勢。故經過多組速度的仿真計算，得到位移邊界層變化趨勢如圖9所示。

圖9 低湍流強度下的位移邊界層厚度

通過對比數據，發現在湍流強度小于1%的時候，位移邊界層變化規律與自由流湍流強度大于1%的位移邊界層變化規律是不同的，平板邊界層位移厚度沿流向的變化曲線隨著湍流強度變大而容易出現下凹或者趨于平緩的現象。

3.4 位移邊界層現象分析

從圖7和圖8中任意選擇一個湍流強度對應的壁面摩擦系數和位移邊界層厚度的曲線考慮，壁面摩擦系數在層流區域是隨流動發展減小的，在湍流區域也是隨著流動發展減小的，但是在轉捩區域卻是增大的(這里取湍流強度為1.26%作為案例說明，大致在X=1.0m處開始轉捩，X=1.2m處完成轉捩)。位移邊界層厚度δ1產生原因是因為壁面無滑移所導致的速度梯度，進而導致流體流過該區域的質量流量損失(圖10中b的陰影區域)。觀察不同X位置處對應的速度輪廓曲線，如下圖11所示：從速度輪廓曲線中知道質量損失分為兩塊(質量損失為水平速度u在Y方向的積分，在圖上表示為面積如圖10(b)的陰影部分)，大約在Y=0.002以下和Y=0.002以上的兩部分(以水平線Y=0.002分界)。隨流動發展，下半部分的面積逐漸減少，上半部分的面積逐漸增大。在過渡區域，一度出現減少大于增大的情況(X=1.0m～1.2m)，上半部分增大的面積小于下半部分增大的面積，此時相當于總質量流量是減少的，故位移邊界層厚度δ1減少。在湍流區，上半部面積顯著增大，總的陰影面積也增大，此時位移邊界層厚度δ1增大。這就解釋了為什么邊界層厚度δ1在轉捩的過渡區域出現先減少后增大的現象。

圖10 邊界層內速度分布

圖11 基于速度輪廓的原因分析

4 結論

通過對零壓力梯度的平板邊界層實驗T3A、T3B[13]，以及Schrader[20]的仿真，模擬了對于不同自由流湍流度的自由流橫掠平板層流到湍流的轉捩過程，分析其計算結果和對比與之相關的實驗得出以下的結論：

1)在來流自由流的湍流強度大于1%的時候，隨著湍流強度的減少，平板的位移邊界層厚度沿著流向方向先是增大而后稍微有下降而后又增大。并且該現象隨著湍流度的降低，越來越明顯。

2)在來流自由流的湍流強度小于1%的時候，隨著湍流強度的減少，平板的位移邊界層變化的趨勢很難捕捉到，基本上是沿著流向方向一直是增大的趨勢。在湍流強度接近1%的時候，位移邊界層厚度沿著流向方向的變化會出現下降趨勢。

3)湍流強度1%是一個分界點，湍流度大于1%，會發生bypass轉捩，即邊界層對自由流脈動的感應而導致的bypass轉捩；而湍流度小于1%發生的是自然轉捩，即邊界層失穩而導致的層流瓦解。由于轉捩的方式不同，而導致的位移邊界層隨湍流度的變化，而導致趨勢變化不同。

4)位移邊界厚度沿著流向方向出現下降或者趨于平緩的變化與轉捩相關。