鄧家鈺，王成恩

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

1 引言

在航空航天領(lǐng)域中，邊界層的層流到湍流的轉(zhuǎn)捩過程的研究是計(jì)算流體力學(xué)(CFD)領(lǐng)域的一個(gè)關(guān)鍵且困難的問題。工程上對(duì)轉(zhuǎn)捩流動(dòng)的模擬大多是采用由實(shí)驗(yàn)總結(jié)出來的經(jīng)驗(yàn)關(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)捩模型對(duì)之前提出的轉(zhuǎn)捩模型進(jìn)行修正，在基于SST k-ω[6]湍流模型的基礎(chǔ)上，融入基于流場(chǎng)當(dāng)?shù)刈兞康霓D(zhuǎn)捩模型。該模型結(jié)合了轉(zhuǎn)捩經(jīng)驗(yàn)關(guān)系式和低雷諾數(shù)湍流模型的優(yōu)勢(shì)，提高了模型的通用性和準(zhǔn)確性。

轉(zhuǎn)捩現(xiàn)象與邊界層外的流動(dòng)密切相關(guān)，對(duì)于平板邊界層來說，當(dāng)自由流湍流強(qiáng)度低于1%發(fā)生自然轉(zhuǎn)捩，即邊界層失去穩(wěn)定導(dǎo)致層流瓦解；當(dāng)湍流強(qiáng)度大于1%的，邊界層對(duì)于自由流脈動(dòng)的感受而導(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)捩實(shí)驗(yàn)。國內(nèi)學(xué)者董平[13]對(duì)經(jīng)典的T3系列[12]平板邊界層實(shí)驗(yàn)進(jìn)行實(shí)驗(yàn)和仿真對(duì)比，顧金生[14]提出了不依賴經(jīng)驗(yàn)系數(shù)的湍流模型，郭玉波等人[15]進(jìn)行了粘性流場(chǎng)的數(shù)值模擬實(shí)驗(yàn)。

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

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

2 理論方法

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

2.1 流體一般控制方程

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

(1)

動(dòng)量守恒方程，即Navier-Stokes方程

(2)

能量守恒方程

(3)

(4)

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

2.2 Transition SST模型的傳輸方程

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

(5)

傳輸源項(xiàng)定義為：

PJ1=2FlengthρS[γFlength]CJ3

(6)

(7)

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

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)捩修正

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

(17)

(18)

γeff=max(γ，γsep)

(19)

(20)

源項(xiàng)定義為

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

Reθt=f(Tu，λ)

(29)

(30)

(31)

第一個(gè)經(jīng)驗(yàn)相關(guān)式子為當(dāng)?shù)赝牧鲝?qiáng)度方程，Thawaite壓力梯度因子λθ定義為

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

(32)

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

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

(33)

Pk=γeffPkP

(34)

(35)

(36)

(37)

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

(38)

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

2.5 指定入口湍流等級(jí)

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

(39)

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

β=0.09，β*=0.0828

時(shí)間尺度由下面定義

(40)

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

(41)

湍動(dòng)能衰減方程可由入口湍流強(qiáng)度和渦流粘性比重新定義為

(42)

2.6 位移邊界層的定義

位移邊界層厚度δ1的計(jì)算公式為

(43)

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

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

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

3.1 T3A和T3B仿真

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

圖1 平板轉(zhuǎn)捩實(shí)驗(yàn)示意圖

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

3.2 仿真結(jié)果

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

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

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

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

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

3.3 Schrader數(shù)值仿真

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

3.3.1 仿真結(jié)果

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

對(duì)于自由流的湍流強(qiáng)度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ù)沿著流向方向的計(jì)算結(jié)果如圖5所示，發(fā)現(xiàn)自由流的湍流強(qiáng)度越大，轉(zhuǎn)捩位置越靠前。計(jì)算其位移邊界層計(jì)算結(jié)果如圖6所示，發(fā)現(xiàn)轉(zhuǎn)捩位置越靠前，位移邊界層有下凹趨勢(shì)或者趨于平緩的位置越靠前，且該位置出現(xiàn)在轉(zhuǎn)捩區(qū)域附近。

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

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

圖6 不同湍流強(qiáng)度下的位移邊界層厚度

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

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

圖8 不同湍流強(qiáng)度下的位移邊界層厚度

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

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

圖9 低湍流強(qiáng)度下的位移邊界層厚度

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

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

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

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

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

4 結(jié)論

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

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

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

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

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