基于轉(zhuǎn)捩模型的位移邊界層研究

鄧家鈺，王成恩

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

1 引言

在航空航天領(lǐng)域中，邊界層的層流到湍流的轉(zhuǎn)捩過程的研究是計算流體力學(xué)(CFD)領(lǐng)域的一個關(guān)鍵且困難的問題。工程上對轉(zhuǎn)捩流動的模擬大多是采用由實驗總結(jié)出來的經(jīng)驗關(guān)系式，國外研究人員提出的一些轉(zhuǎn)捩經(jīng)典模型，如Abu-Ghanam&Ahaw(AGS)模型、Steelant & Dik模型、Suzen & Huang(S-H)模型、Menter & Langry(M-L)[1-5]模型等。其中M-L轉(zhuǎn)捩模型對之前提出的轉(zhuǎn)捩模型進(jìn)行修正，在基于SST k-ω[6]湍流模型的基礎(chǔ)上，融入基于流場當(dāng)?shù)刈兞康霓D(zhuǎn)捩模型。該模型結(jié)合了轉(zhuǎn)捩經(jīng)驗關(guān)系式和低雷諾數(shù)湍流模型的優(yōu)勢，提高了模型的通用性和準(zhǔn)確性。

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

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

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

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

2 理論方法

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

2.1 流體一般控制方程

質(zhì)量守恒方程，即連續(xù)性方程

(1)

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

(2)

能量守恒方程

(3)

(4)

其中，Φ是粘性導(dǎo)致的能量耗散函數(shù)。

2.2 Transition SST模型的傳輸方程

傳輸方程的間歇系數(shù)γ定義如下

(5)

傳輸源項定義為：

PJ1=2FlengthρS[γFlength]CJ3

(6)

(7)

其中，S為應(yīng)變速率大小，F(xiàn)length是一個控制轉(zhuǎn)捩區(qū)域長度經(jīng)驗相關(guān)的因子。消減項/再層流化項被定義為

PJ2=(2CJ1)ρΩγFturb

(8)

Er3=Cr3

(9)

其中，Ω是渦量大小。轉(zhuǎn)捩的開始由以下的系數(shù)控制

(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 轉(zhuǎn)捩修正

對分離引起的轉(zhuǎn)捩的修正如下

(17)

(18)

γeff=max(γ，γsep)

(19)

(20)

源項定義為

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

Reθt=f(Tu，λ)

(29)

(30)

(31)

第一個經(jīng)驗相關(guān)式子為當(dāng)?shù)赝牧鲝姸确匠?，Thawaite壓力梯度因子λθ定義為

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

(32)

2.4 耦合轉(zhuǎn)捩模型和SST傳輸方程

給出SST傳輸方程后，為了計算更為準(zhǔn)確，需要將SST傳輸方程和轉(zhuǎn)捩模型耦合

(33)

Pk=γeffPkP

(34)

(35)

(36)

(37)

Ft=max(Florig，F(xiàn)3)

(38)

其中，Dk和Pk為SST模型的原始生成項和原始銷毀項，F(xiàn)1orig為原始SST混合方程。

2.5 指定入口湍流等級

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

(39)

對于在自由流中SST湍流模型的常數(shù)為

β=0.09，β*=0.0828

時間尺度由下面定義

(40)

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

(41)

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

(42)

2.6 位移邊界層的定義

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

(43)

其中，u為流體速度，u∞為無窮遠(yuǎn)處自由流速度。

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

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

3.1 T3A和T3B仿真

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

圖1 平板轉(zhuǎn)捩實驗示意圖

表1 實驗的進(jìn)口參數(shù)

3.2 仿真結(jié)果

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

圖2 T3A 的摩擦系數(shù)

圖3 T3B的摩擦系數(shù)

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

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

3.3 Schrader數(shù)值仿真

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

3.3.1 仿真結(jié)果

自由流湍流強度大于1%的時候發(fā)生bypass轉(zhuǎn)捩，自由流湍流強度小于1%的時候發(fā)生自然轉(zhuǎn)捩。分別對自由流湍流強度大于1%和湍流強度小于1%計算和討論，分析計算結(jié)果和進(jìn)行比較。

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

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

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

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

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

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

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

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

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

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

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

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

3.4 位移邊界層現(xiàn)象分析

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

圖10 邊界層內(nèi)速度分布

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

4 結(jié)論

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

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

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

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

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