縱掠平板速度和溫度邊界層湍流轉(zhuǎn)捩區(qū)的積分方法

趙 波， 劉 建， 李開勇

(四川大學機械工程學院 空天動力燃燒與冷卻教育部工程研究中心， 成都 610065)

1 引 言

邊界層由簡單穩(wěn)定狀態(tài)的層流向具有強烈非線性現(xiàn)象、多尺度效應(yīng)的湍流過渡區(qū)域稱為轉(zhuǎn)捩區(qū). 長期以來，邊界層轉(zhuǎn)捩一直是流體力學中最重要的前沿問題和難點之一[1-4]. 邊界層轉(zhuǎn)捩伴隨著壁面摩擦和傳熱特性的急劇變化，減小飛行器表面摩擦阻力和熱流、燃氣渦輪發(fā)動機葉片冷卻和燃燒時燃料摻混效率等問題均與轉(zhuǎn)捩特性密切相關(guān)[5-6]；轉(zhuǎn)捩對揭示湍流真正形成機理也具有重要價值[7]. Reynolds最早通過圓管實驗發(fā)現(xiàn)有層流和湍流兩種不同流態(tài)，這被視為研究邊界層轉(zhuǎn)捩問題的肇始. Emmons[8]發(fā)現(xiàn)轉(zhuǎn)捩區(qū)間歇出現(xiàn)的湍斑，它們沿流動方向不斷變大增多，最終在完全湍流區(qū)布滿，這表明邊界層轉(zhuǎn)捩區(qū)以間歇出現(xiàn)的湍斑為主要特征. 為此，Dhawan和Narasimha[9]提出間歇因子的概念以描述轉(zhuǎn)捩區(qū)湍斑的分布特征：間歇因子為0時流動為層流，為1時則代表完全湍流，他們以間歇因子為加權(quán)系數(shù), 可將任意位置的流場視為層流與湍流的線性組合[10]，這種疊加思想能較好地預(yù)測零壓和順壓梯度的轉(zhuǎn)捩過程[4]，但并未考慮層流與湍流間的相互作用. 針對這個不足，Libby[11]采用條件平均法，通過建立含間歇因子的輸運方程加以解決. 張若凌等[12]采用熱力學平衡系統(tǒng)連續(xù)相變方法，類比研究了管內(nèi)轉(zhuǎn)捩區(qū)的流動和振蕩特性，視轉(zhuǎn)捩流動為層流和湍流的線性組合，所提出的合成系數(shù)考慮了振蕩特性. Langtry等[13-14]發(fā)展了完全由局部變量構(gòu)造的新型模式,使轉(zhuǎn)捩計算與流場結(jié)構(gòu)直接相關(guān)，與Coupland[15]的空氣外掠平板轉(zhuǎn)捩實驗的結(jié)果對比表明準確性較高. 董平等[16-17]在上述模型基礎(chǔ)上，分別對外掠平板邊界層轉(zhuǎn)捩區(qū)的速度和溫度進行了數(shù)值計算，并與Coupland[15]試驗結(jié)果進行了比較. 針對完全湍流流動，Prandtl和Taylor沿邊界層厚度方向?qū)⑦吔鐚觿澐譃閷恿鞯讓雍屯牧骱诵膮^(qū)，提出分別描述速度和溫度邊界層的Prandtl-Taylor兩層模型；von Kármán進一步在層流底層和湍流核心區(qū)中間增加了緩沖層，建立了三層模型. Khademi等[18-19]采用完全湍流邊界層的兩層模型，通過積分方法分別建立了動量方程和能量方程，獲得的速度場和溫度場分布規(guī)律與以往理論結(jié)果一致性較好. 最近，Lozano-Durán等[20]利用非線性拋物化理論預(yù)測了轉(zhuǎn)捩的起始位置，預(yù)測結(jié)果與大渦模擬(LES)和直接數(shù)值模擬(DNS)計算吻合度較高. 由上可見，目前對邊界層轉(zhuǎn)捩區(qū)的研究多以數(shù)值模擬方法為主，通常計算量較大(如LES和DNS相當于Re9/4[3])，由于問題的復雜性，采用理論研究的方法相對較少. 就作者所知，目前為止尚未發(fā)現(xiàn)采用積分方法預(yù)測邊界層自然轉(zhuǎn)捩區(qū)的速度和溫度場分布方面的理論研究. 為此，本文在以往研究基礎(chǔ)上[21-22]，將以空氣外掠平板邊界層為研究對象，采用完全湍流區(qū)的兩層模型思想，首先將轉(zhuǎn)捩區(qū)沿邊界層厚度方向劃分為層流底層和準湍流層，然后在準湍流層利用間歇因子定量描述該區(qū)域的湍流時間份額，最后利用積分方法，分別建立轉(zhuǎn)捩區(qū)的動量方程和能量方程，以獲得邊界層轉(zhuǎn)捩區(qū)速度場和溫度場的理論分布，同時將之與Coupland實驗結(jié)果和數(shù)值仿真模型相比較，以證明理論模型的正確性. 邊界層轉(zhuǎn)捩最關(guān)注的因素有三：轉(zhuǎn)捩的起始位置、轉(zhuǎn)捩區(qū)長度和轉(zhuǎn)捩區(qū)內(nèi)的時均流動特性[23]. 本文將僅限于轉(zhuǎn)捩區(qū)的時均流動特性和溫度場分布規(guī)律的討論，研究對象限于零壓力梯度條件下外掠平板的自然轉(zhuǎn)捩特性.

2 理論模型

2.1 轉(zhuǎn)捩邊界層積分方程組

首先，這里討論的流體指粘性流體，根據(jù)牛頓內(nèi)摩擦定律、層流理論和Prandtl-Taylor湍流邊界層的兩層理論[24]，靠近壁面的區(qū)域因液體粘性的存在，可近似認為是層流流動[10]，故有理由認為對于從層流向湍流過渡的轉(zhuǎn)捩區(qū)也應(yīng)如此，為與完全湍流的“層流底層”稱謂相一致，仍稱轉(zhuǎn)捩區(qū)的近壁面區(qū)域為“層流底層”(laminar sublayer)，如圖1所示. 其次，根據(jù)間歇理論[8-10]，轉(zhuǎn)捩區(qū)內(nèi)的湍斑在流動方向上逐漸增多，直到在完全湍流區(qū)布滿，若遵照完全湍流區(qū)的Prandtl-Taylor兩層理論，布滿湍斑的區(qū)域應(yīng)僅限于湍流核心區(qū)而不包含層流底層；類似地，這里假設(shè)在轉(zhuǎn)捩區(qū)內(nèi)，只有外層區(qū)域(從層流底層外緣到主流區(qū))內(nèi)存在湍斑，為方便，稱之為“準湍流層”(quasi-turbulent layer)，如圖1. 最后，根據(jù)Langtry和Menter[13-14]的觀點：用動量厚度雷諾數(shù)來判斷轉(zhuǎn)捩的開始，用間歇因子來描述轉(zhuǎn)捩的過程[3]，假設(shè)轉(zhuǎn)捩區(qū)初次出現(xiàn)湍流的位置在動量邊界層最大厚度δ*處，如圖1，并假設(shè)：(1)常物性和不可壓縮流體；(2)準湍流層湍斑在主流方向的脈動速度、脈動溫度與完全湍流區(qū)規(guī)律相同[12]；(3)零壓力梯度和常壁溫，忽略耗散熱；(4)經(jīng)計算可知，層流底層和準湍流層的速度邊界層厚度非常接近溫度邊界層的對應(yīng)厚度，為積分方便，設(shè)層流底層和準湍流層的溫度邊界層厚度Δ1和Δ略大于速度邊界層的對應(yīng)厚度δ1和δ；(5)只考慮自然轉(zhuǎn)捩情況，為與Coupland實驗結(jié)果相比較[14-15]，設(shè)轉(zhuǎn)捩起始位置雷諾數(shù)為105，在該位置處的層流底層厚度記為δ*，如圖1.

圖1 空氣外掠平板轉(zhuǎn)捩區(qū)邊界層的兩層模型示意圖

圖2給出轉(zhuǎn)捩區(qū)控制體積示意，區(qū)域1-2-3-4為層流底層，3-4-5-6為準湍流層，l為流體厚度，dx為x向微元. 因dx極小，認為層流底層厚度δ1在dx內(nèi)沿x向近似不變.

圖2 轉(zhuǎn)捩區(qū)控制體積微元

根據(jù)動量守恒定理，采用與層流邊界層類似的積分方法[24]，最后獲得到外掠平板轉(zhuǎn)捩區(qū)的動量方程為:

(1)

式中ρ為流體密度(kg/m3)，μ為動力粘度(Pa·s)，δ1、δ分別為層流底層和準湍流層速度邊界層的厚度(m)，見圖1，u1、u2分別為層流底層和準湍流層的時均速度(m/s)，uL和u∞分別為層流底層外緣處和主流區(qū)的速度(m/s). 根據(jù)能量守恒定理，采用完全類似的積分法[24]，得到能量方程為:

(2)

式中c為定壓比熱容[J/(kg·K)]，λ為導熱系數(shù)[W/(m·K)]，Δ1、Δ分別為層流底層和準湍流層溫度邊界層的厚度(m)，T1、T2分別為層流底層和準湍流層的時均溫度(K)，TL和T∞分別為層流底層外緣處和主流區(qū)的溫度(K).

2.2 速度分布函數(shù)

認為外掠平板轉(zhuǎn)捩邊界層具有相似的速度分布，設(shè)在層流底層(0≤y<δ1)和準湍流層(δ1≤y<δ)的時均速度分布具有如下形式:

(3)

式中ai(i=1～5)為待定系數(shù)，m為正的常數(shù).

(4)

γ=1-e-0.412(x-xt)2/χ2

(5)

式中xt為轉(zhuǎn)捩區(qū)起始位置，由臨界雷諾數(shù)確定；χ=x|γ=0.75-x|γ=0.25，由文獻[8]中零壓力梯度條件下外掠平板轉(zhuǎn)捩邊界層的實驗數(shù)據(jù)確定.

由上述邊界條件式并考慮方程(4)和(5)，最終確定方程(3)中的待定系數(shù)為:

a1=a3=0,

a5=u∞δ-1/m.

相應(yīng)地，轉(zhuǎn)捩區(qū)的速度分布式(3)可寫為:

(6)

2.3 溫度分布函數(shù)

同樣地，認為轉(zhuǎn)捩邊界層具有相似的溫度分布，考慮假設(shè)(3)，設(shè)層流底層和準湍流層的過余時均溫度分布為:

(7)

b1=b3=0,

相應(yīng)地，轉(zhuǎn)捩區(qū)的時均溫度分布式(7)可寫為:

(8)

3 有限元模型和常數(shù)m、n的選擇

3.1 有限元模型

以外掠平板的空氣為研究對象，在Fluent軟件中采用如圖3所定義的有限元分析區(qū)域：平板長度取3 600 mm，空氣域高度取220 mm. 所研究的轉(zhuǎn)捩邊界層位于平板近壁面，需精確計算近壁面速度和溫度變化情況，經(jīng)仔細比較[26-28]，最終選取FLUENT軟件中的SSTk-ω模型作為數(shù)值計算的轉(zhuǎn)捩模型.

圖3 空氣外掠平板數(shù)值計算區(qū)域

為了同Coupland實驗結(jié)果[14-15]比較，取自由流速度和溫度分別為u∞=5.3 m/s和T∞=293 K，常壁溫Ts=323 K. 空氣入口位于數(shù)值模型最左端，空氣出口分別位于模型最右端和空氣域上方，兩者皆采用零壓力出口，如圖3. 二維有限元模型采用四邊形面單元，考慮到所關(guān)注的邊界層僅存在于平板近壁面極薄區(qū)域，為節(jié)省計算資源且保證求解精度，取網(wǎng)格全局尺寸10 mm，在平板垂直距離60 mm內(nèi)設(shè)置40層網(wǎng)格，每層網(wǎng)格以1.2倍速度增加，總單元數(shù)19 712，總結(jié)點數(shù)20 121，如圖4所示. 根據(jù)假設(shè)(5), 轉(zhuǎn)捩起始位置處的雷諾數(shù)為105，針對上述給定參數(shù)，數(shù)值和理論計算的結(jié)果表明，轉(zhuǎn)捩邊界層結(jié)束時的雷諾數(shù)為3×105.

圖4 空氣外掠平板有限元模型網(wǎng)格

3.2 常數(shù)m、n的確定

由第1部分可知，在轉(zhuǎn)捩邊界層的準湍流層速度和溫度分布式(6)和(8)中，冪指數(shù)項1/m和1/n的大小對理論預(yù)測至為關(guān)鍵. 為此，經(jīng)過大量計算、選擇和對比，限于篇幅，僅以m(或n)分別取7、7.5和8為例闡述之. 同時，將本文理論結(jié)果與Coupland平板轉(zhuǎn)捩實驗和上述數(shù)值模型的結(jié)果相對比，以選取合適的m和n常數(shù)值.

圖5 轉(zhuǎn)捩區(qū)內(nèi)壁面的局部摩擦系數(shù)

圖6 轉(zhuǎn)捩區(qū)時均速度沿邊界層厚度方向的分布

圖7給出轉(zhuǎn)捩區(qū)時均溫度沿邊界層厚度方向的變化，通過對比發(fā)現(xiàn)，溫度當n=7時理論解與數(shù)值解(以之為benchmark)最接近. 因此最終確定：選取m=7.5，n=7.

圖7 轉(zhuǎn)捩區(qū)時均溫度沿邊界層厚度方向的分布

相應(yīng)地，速度和溫度分布式(6)(8)變?yōu)?/p>

(9)

(10)

由轉(zhuǎn)捩邊界層的動量厚度?的定義[9, 24]，有:

將式(9)代入上式，若不考慮間歇粘度比?沿y向的變化，轉(zhuǎn)捩區(qū)的無量綱動量厚度可表示為:

由上式可直接寫出動量厚度雷諾數(shù)的顯式表達式.

4 Dhawan和Narasimha的線性疊加思想

為后面對比方便，這里簡要介紹一下Dhawan和Narasimha關(guān)于轉(zhuǎn)捩區(qū)的線性疊加思想. 根據(jù)間歇理論[9-10]，以間歇因子γ為權(quán)重系數(shù)，轉(zhuǎn)捩邊界層的速度和溫度分布可以線性疊加為下式：

u=(1-γ)uL+γuT

θ=(1-γ)θL+γθT

(11)

式中uL、uT分別代表層流速度和完全湍流流動時的時均速度，θL、θT分別是層流過余溫度和完全湍流的過余時均溫度. 對于層流和湍流的速度分布，這里分別采用von Kármán的層流邊界層的積分解和1/7冪律[24, 29]，即

同樣地，由式(11)容易推導出轉(zhuǎn)捩區(qū)努塞爾數(shù)Nu可由層流和完全湍流的努塞爾數(shù)線性疊加得到：

Nu=(1-γ)NuL+γNuT

(12)

其中層流努塞爾數(shù)NuL=0.332Pr1/3Re1/2[24], 湍流努塞爾數(shù)NuT=0.0287Pr0.6Re0.8[30]. 為方便，以下稱方程(11)和(12)為D-N解.

5 結(jié)果討論

5.1 轉(zhuǎn)捩邊界層厚度

圖8給出了通過動量積分方程組(1)獲得的外掠平板動量邊界層中層流底層和準湍流層厚度的理論分布. 由圖8可見，層流底層從層流邊界層的最大厚度δ*(=0.0049 m)處開始突降，直到轉(zhuǎn)捩區(qū)結(jié)束(這里轉(zhuǎn)捩區(qū)結(jié)束時雷諾數(shù)為3×105)后進入完全湍流區(qū)，而準湍流層則從相同的邊界層厚度δ*處開始驟增，直到轉(zhuǎn)捩區(qū)結(jié)束后進入完全湍流區(qū). 圖9給出了通過能量積分方程組(2)獲得的外掠平板熱邊界層中層流底層和準湍流層厚度的理論分布(其中Δ*=0.005 4 m)，主要結(jié)論與動量邊界層類似，不贅述.

圖8 速度邊界層中層流底層和準湍流層的厚度

5.2 時均速度

圖10給出x位置固定時(Re=2.5×105)，外掠平板轉(zhuǎn)捩區(qū)內(nèi)空氣的時均速度沿邊界層厚度方向的分布. 由圖可見，本文理論解與D-N解大體上符合得很好，尤其在準湍流層區(qū)域，而在層流底層內(nèi)二者略有偏差.

圖10 時均速度沿轉(zhuǎn)捩邊界層厚度方向的分布

圖11給出了平板不同x位置處(對應(yīng)不同的γ)的時均速度分布情況. 由圖可見，γ表觀上對近壁面的層流底層速度分布有更大的影響，這與本文理論公式(9)，以及Dhawan和Narasimha[9]的發(fā)現(xiàn)相一致：即γ主要影響x方向的時均速度分布，且在近壁面區(qū)域影響相對明顯. 進一步觀察發(fā)現(xiàn)，γ對遠離壁面的準湍流層速度分布完全沒有影響，這是由于圖中采用無量綱厚度y/δ坐標的緣故，實際上不同的γ會導致不同的δ，具體見公式(9). 注意到圖11中的速度呈相似性分布，γ=0和γ=1曲線分別代表層流和完全湍流流動兩種特殊情況，其中γ=0代表的層流是計算轉(zhuǎn)捩開始的初始條件，轉(zhuǎn)捩開始后，當x(和γ)稍有所增大，速度的變化就十分劇烈，這符合轉(zhuǎn)捩區(qū)速度(和壁面摩擦系數(shù))的突變驟增的特點，注意到速度總體分布與文獻[9]中的圖11結(jié)果幾乎完全相同，只不過后者采用了另一個無量綱厚度y/?作為橫坐標，相當于放大了圖11中的層流底層區(qū)域.

圖11 間歇因子γ對無量綱時均速度的影響

圖12給出了無量綱速度在轉(zhuǎn)捩區(qū)層流底層和準湍流層隨間歇粘度比?參數(shù)變化情況. 由圖12可見，除轉(zhuǎn)捩剛開始時層流底層速度驟升外，其余大部分情形均是速度隨?增加而減小，且兩個區(qū)域下降幅度非常相似. 這是因為：根據(jù)間歇粘度比?定義，y固定時?只與間歇因子γ(或x)有關(guān)，當γ(或x)增大時無量綱厚度y/δ變小，由圖10易知時均速度也相應(yīng)減小.

圖12 間歇粘度比?對無量綱時均速度的影響

5.3 時均溫度

圖13給出x位置固定時(Re=2.5×105)時均溫度沿邊界層厚度方向的分布. 由圖可見，同無量綱速度的分布幾乎完全相同，本文理論解與D-N解大體上符合得很好. 圖14和15給出了平板不同x位置處(對應(yīng)不同的γ，y固定時也對應(yīng)不同的間歇導溫比ψ)的時均溫度分布，主要結(jié)論與前面無量綱速度的分布情形非常相似，不贅述.

圖13 時均溫度沿轉(zhuǎn)捩區(qū)邊界層厚度方向的分布

圖14 間歇因子γ對無量綱時均溫度的影響

5.4 壁面摩擦系數(shù)

根據(jù)Blasius公式[24]，層流時壁面摩擦系數(shù)CfL=0.664Re-0.5；由Prandtl-Schlichting理論公式[24, 29]，完全湍流時的壁面摩擦系數(shù)CfT=0.0587Re-0.2. 基于間歇理論，Emmons[8]以及其他學者[9]給出了轉(zhuǎn)捩區(qū)壁面摩擦系數(shù)Cft公式：Cft=(1-γ)CfL+γCfT，為方便，稱之為Emmons解.

層流向湍流轉(zhuǎn)捩最顯著的特征是：壁面摩擦系數(shù)取得極小值處轉(zhuǎn)捩開始，在轉(zhuǎn)捩區(qū)內(nèi)突變驟增，轉(zhuǎn)捩結(jié)束位置處達到極大值，且略大于完全湍流時的壁面摩擦系數(shù)CfT，隨后進入完全湍流，壁面摩擦系數(shù)緩慢減小，該特征也是轉(zhuǎn)捩實驗中最易監(jiān)測的變量(見圖16). 將本文理論解和數(shù)值解獲得的局部摩擦系數(shù)，同Coupland實驗結(jié)果、層流和完全湍流理論解以及轉(zhuǎn)捩區(qū)的Emmons解進行了對比，如圖16. 由圖16可見，本文的理論解與試驗結(jié)果符合得非常好，轉(zhuǎn)捩區(qū)內(nèi)最大相對誤差為4.8%，數(shù)值解也與試驗結(jié)果大體符合. 注意到轉(zhuǎn)捩區(qū)內(nèi)Emmons解與試驗值相比偏差相對更大，這表明用積分方法獲得的理論解同以往用層、湍流疊加的方法(即Emmons解)具有較大優(yōu)勢.

圖15 間歇導溫比ψ對無量綱時均溫度的影響

圖16 平板壁面摩擦系數(shù)分布

5.5 努塞爾數(shù)

圖17給出了代表壁面?zhèn)鳠崽匦缘呐麪枖?shù)沿主流方向的分布情況，發(fā)現(xiàn)同壁面摩擦系數(shù)在轉(zhuǎn)捩區(qū)內(nèi)突變驟增的特征不同，努塞爾數(shù)在轉(zhuǎn)捩區(qū)雖有變化，但波動要相對小得多. 圖17將本文理論解和數(shù)值解獲得的局部努塞爾數(shù)，同D-N解進行了對比，結(jié)果表明理論解和D-N解及數(shù)值解結(jié)果符合得非常好，在轉(zhuǎn)捩區(qū)內(nèi)最大相對誤差分別為3.9%和5.3%，再一次證明了理論模型的正確性.

圖17 局部努塞爾數(shù)分布Fig.17 Profiles of local Nusselt numbers

6 結(jié) 論

以常壁溫、零壓力梯度的外掠平板空氣為研究對象，將自然轉(zhuǎn)捩邊界層沿厚度方向劃分為層流底層和準湍流層兩部分，采用積分方法，分別建立了關(guān)于速度和溫度的動量和能量積分方程組，并獲得了速度場和溫度場的顯式表達式，通過四階龍格-庫塔算法獲得了動量和熱邊界層厚度大小，同時與數(shù)值解、D-N解和實驗結(jié)果加以對比表明，在給定條件下，通過理論解獲得的表面摩擦系數(shù)和努塞爾數(shù)與benchmark結(jié)果的最大相對誤差在4.8%以內(nèi)，證明了理論模型的正確性，主要結(jié)論如下：

(1) 利用相似假設(shè)和湍流邊界層厚度分層思想，分別采用三次多項式和1/m(或1/n)冪指數(shù)函數(shù)代表轉(zhuǎn)捩區(qū)層流底層和準湍流層的速度場和溫度場的分布規(guī)律，發(fā)現(xiàn)常數(shù)m=7.5、n=7時具有較好的理論預(yù)測精度.

(2) 速度和溫度理論公式表明，間歇因子不僅對邊界層(準湍流層)厚度有影響，尤其對“層流底層”的厚度、速度和溫度分布具有關(guān)鍵影響，因此，該底層不是真正的層流流動；參數(shù)分析表明，所提出的間歇粘度比?和間歇導溫比ψ是影響轉(zhuǎn)捩區(qū)速度和溫度場分布的關(guān)鍵參數(shù).

(3) 解析模型從理論角度再次證實，在轉(zhuǎn)捩區(qū)表面摩擦系數(shù)具有突變驟增的顯著特征，而且在轉(zhuǎn)捩起始和結(jié)束位置處分別取得極小值和極大值. 而努塞爾數(shù)(對流換熱系數(shù))在轉(zhuǎn)捩區(qū)的變化遠沒有壁面摩擦系數(shù)的變化顯著，因此表面摩擦系數(shù)應(yīng)視為判斷邊界層轉(zhuǎn)捩起始位置和結(jié)束位置的關(guān)鍵指標. 以上結(jié)論同樣有潛力用于外掠多孔壁面轉(zhuǎn)捩邊界層的發(fā)散冷卻等相關(guān)轉(zhuǎn)捩問題研究中.