強吸氣旋轉圓筒壁面湍流邊界層建模及計算1)

Ievgen Mochalin 林靜雯 蔡建程,2) Volodymyr Brazhenko 鄂世舉

*(浙江師范大學工學院,浙江金華 321004)

?(浙江省城市軌道交通智能運維技術與裝備重點實驗室,浙江金華 321005)

引言

旋轉可滲透圓柱外的流動是在Taylor-Couette 流動的基礎上疊加徑向穿透流動[1-2],它具有很多應用,旋轉過濾是一個典型例子[3-5].含雜質溶液通過旋轉圓筒膜可以大為提高過濾效率,并且旋轉膜表面的固相沉積率低不容易阻塞,因此有學者專門分析旋轉過濾圓柱表面的流體運動[6-10].

旋轉柱體的冷卻是Taylor-Couette 流動的另一個研究方向,F′enot 等[11]綜述了轉子-靜子系統中的同心圓筒間隙內的不可壓縮流動與換熱問題研究.Mochalin 等[12]研究了外筒靜止、內筒旋轉工況下的間隙流動,發現內筒做成完全多孔或是條狀多孔可以使換熱效率提高1.5 倍.因此,研究旋轉多孔圓筒表面附近的流動很有意義.

旋轉圓筒外的流動由于受到離心不穩定的作用,會導致二次渦流即Taylor 渦的產生.對于同心不可滲透圓筒間的流動,在較低的旋轉速度時就能發生流動不穩定性,體現為軸對稱Taylor 渦的產生,隨著旋轉速度的提高,流動模式會改變[13].對于經典Taylor-Couette 流動,相比于由離心不穩定弓起的Taylor 渦產生,兩圓筒間的流動由層流到湍流的轉捩發生較遲.通過多孔圓筒弓入徑向穿透流,可以使間隙流動變得穩定,從而Taylor 渦的產生大為延后[14-15].Mochalin 和Khalatov[16]的研究表明強吸氣穿透流可以防止Taylor 渦的產生直到非常大的旋轉速度,并且發現旋轉吸氣圓筒壁面弓起的旋轉流體運動集中在其壁面附近,而外面靜止圓筒對里邊旋轉圓筒附近的流動影響很小,所以旋轉流動可以視為旋轉多孔圓筒表面的邊界層流動.邊界層厚度依賴于旋轉速度和徑向穿透流強度.

研究表明,多孔壁面徑向流可以使邊界層從層流轉捩到湍流發生在離心不穩定(Taylor 渦) 之前,這就為邊界層湍流強度控制提供了新思路.旋轉壁面流動具有很多應用價值,比如相分離、表面冷卻等[12,16].所以,有必要研究旋轉多孔圓筒表面的湍流邊界層流動.

Mochalin 和Khalatov[16]通過求解Reynolds 平均Navier-Stokes 方程和Reynolds 應力湍流模型,研究了同心圓筒間的流動.梁田等[17]則利用Spalart-Allmaras 湍流模型計算邊界層抽氣,進行葉柵角區分離流動控制的研究.在這些研究中,近壁面流動使用非常稠密的低Reynolds 數計算網格,計算量大.本文的研究目標是對于強吸氣旋轉壁面流動使用簡單而可靠的預測方法,并提供計算步驟.其研究思路是使用恰當的代數湍流模型,基于詳細數值模擬結果進行代數湍流模型參數的調整和校準,使其適用于一定工況范圍.

1 強吸氣壁面旋轉流動速度分布的解析表達式

1.1 基本假設與無量綱參數化

本文所研究的同心圓筒間流動為半徑R外的多孔圓柱壁外的穩定流動(在離心不穩定發生之前).旋轉流體運動主要集中在厚度為δ 的邊界層內.多孔壁面吸氣,假設均勻吸氣速度分布Vs.旋轉流動認為在軸向不變化,即Vz≡0,?/?z≡0,且軸對稱即?/?φ ≡0(φ 為方位角).

旋轉流動可使用Reynolds 平均Navier-Stokes(NS) 方程聯合連續方程進行描述.顯然,使用柱坐標系(r,φ,z)更為方便,在柱坐標系下穩態不可壓縮N-S方程

式中,Vr,Vφ分別為徑向和周向的平均速度,ρ 為流休密度,νe為流體的有效運動黏度,它等于分子黏性ν和湍流黏性νt之和,νe=ν+νt.

連續方程為

積分上述方程,注意到壁面Vr(R)=Vs條件,得到

利用式(2),方程(1)可以轉化成

在本文中,邊界層的計算在徑向劃分為多個微元段Δr,在每一段中νe認為是常數.遵循這一假設,式(4)中最后一項在每個微元段可以忽略.進一步對式(4)中最后一項在整個邊界層內都忽略,這是附加的假設.其有效性將從本文方法的結果與詳細的數值仿真結果相比較看出.利用式(3) 以及上述假設,積分式(4),得到

上式中的積分常數C可由旋轉壁面的邊界條件

確定,其中Ω 為壁面旋轉角速度.式(5)中的速度梯度可用壁面處的摩擦應力

來表示.

聯合式(5)、式(6)和式(8),得到

為了得到無量綱表達,使用圓筒半徑R、壁面旋轉速度W和動壓作為長度、速度和壓力的參考量,通過下式定義無量綱量

其中,ξ 為無量綱徑向坐標,ur,uφ為無量綱徑向和周向速度,q為無量綱壓力.從而得到無量綱化的式(9)為

式中,uw=Vw/W為無量綱摩擦速度,us=Vs/W為無量綱吸氣速度.式(11)兩邊同除以RW,得到

為求解式(12),將旋轉圓柱壁面法向的邊界層分成足夠小間隔Δξi,在每個間隔內湍流黏性看成為常數,式(12)簡化為

積分后得到

積分常數C由下列邊界條件得到

式中ξ0i,uφi分別為無量綱徑向距離和每個徑向間隔Δξi的周向速度;對于第一步,由式(6)、式(10)可知ξ01=1,uφ1=1.把式(14)代入式(13),得到常數C的定義

因此,在每個間隔內ξ ∈[ξ0i,ξ0i+Δξi],周向速度可由式(13)、式(15)得到

1.2 圓筒表面摩擦應力

壁面吸氣使得兩個圓柱面間流動的Taylor 渦產生大為滯后.這意味著徑向Reynolds 數(Rer)足夠大時,兩圓柱面間流動為層流流動,其周向流動速度分布可表示為[16]

上式沿徑向求導,得到

由該速度分布,進而可以得到摩擦系數

表1 給出了所考慮的旋轉圓筒壁面摩擦系數.表中的數值仿真按文獻[16] 中的方法得到,可以看出由式(18) 得到的結果與數值仿真結果非常接近,誤差在1%之內.表1 中的第1、3 行的Reφ,Rer組合,其流態為層流,其余為湍流.基于此,式(18)可以應用于旋轉同心圓筒間的湍流流動.所以,有如下關系式

利用式(19)可以化簡式(16),得到

這里順便指出,已有的對吹吸氣槽道湍流的研究表明吹氣可以使壁面阻力系數下降,而吸氣會使壁面阻力系數上升[18-19],而高頻周期吹氣擾動在狹縫下游產生明顯的減阻效果[20],另外通過在平板壁面施加不同頻率振幅的周期性擾動,能進行湍流邊界層的主動控制減阻[21].表1 的第2 行和第3 行表明,固定旋轉速度(Reφ一定)的前提下,摩擦系數隨著吸氣強度Rer的增加而增加.而隨旋轉速度Reφ的提高,摩擦系數下降.所以當Reφ和Rer同時提高時,摩擦系數有可能下降.

表1 強吸氣旋轉圓筒表面摩擦應力Table 1 Friction coefficient values on the surface of a rotating cylinder with strong suction

2 代數湍流模型

2.1 湍流黏性的基本關系式

壁面吸氣可以延遲邊界層從層流到湍流的轉捩[22-25].本文研究邊界層除吸氣外還涉及離心力和旋轉對稱等因素影響.文獻[16] 中的數值仿真誠然可以考慮這些因素,但計算量過大.本文提供更為快捷的代數湍流型模型,以求解旋轉邊界層.思路如下:尋找一個合適的代數湍流模型,根據新增加的影響因素進行改進,并依據特定的情況用實驗或計算數值結果進行參數校準.

本文依據Cebeci-Smith(CS)湍流模型[26],校準數據使用文獻[16] 的結果.CS 模型能考慮壓力梯度、壁面吸(吹)氣、可壓縮性和低Reynolds 數效應等因素的影響.對每個影響因素,分別弓入一個因子(經驗參數或是含經驗參數的表達式),然后耦合所有的因子.

平板流動的經典CS 兩層湍流模型可以表述如下

式中νti,νto為湍流邊界層內層與外層的湍流運動黏性系數;ym為νti=νto時的壁面法向距離;Vx,Vy為流向和橫向的速度;lm為Prandtl 混合長度;y為離壁面的法向距離;Vm為邊界層厚度δ 內的最大速度;δ*為邊界層位移厚度.無量綱距離y+定義為

經驗參數

式(22)可從下式

廣義的CS 模型[26]可以考慮如下因素:趨于邊界層外界線時湍流強度的衰減、外層低Reynolds 數對湍流渦的影響、壁面吸(吹)氣、壓力梯度、可壓縮性.本文不考慮最后兩個因素,僅考慮前3 者.

邊界層外層界線附近湍流強度的衰減,由Klebanoff間歇系數來表示

用它乘式(24)右端.

低Reynolds 數的影響通過改變式(24) 右端的Clauser 系數k實現

式中尾跡參數Π0由Schlichting[24]建議為

式中

δ**為動量厚度

廣義CS 模型中通過改變式(23)參數A*以考慮壁面吸氣的影響

CS 模型中并沒有考慮流線曲率的影響,流線曲率會產生離心力.本文利用Bradshaw[27]的思想,這一思想被許多學者所采納[28-30],它把流線曲率的影響與浮力效應相類比,弓入Richardson 數Ri(表示由于離心力和速度梯度產生的湍流能量之比) 到湍流黏性系數中

式中,νt0為不考慮流線曲率的湍流黏性,β,m為經驗常數.β,m的選擇依賴于用于校準的剪切流,這將在下文中簡述.根據Bradshaw[27]的建議,對于具有周向和軸向分量的旋渦流動

所以,在廣義CS 模型中,考慮到趨于邊界層外界線時湍流強度的衰減、外層內低Reynolds 數對湍流渦的影響、流線曲率、壁面吸(吹)氣等方面后,式(22)和式(24)改寫為

式(23)中的A*由式(34)給出,而式(38)中的k由式(30)～式(33)得到而不是式(27)中的常數0.016 8.式(38)中的Klebanoff函數Fkl由式(29)給出.

利用式(3)和式(10),上式可以寫為

對于壁面法向距離y及無量綱壁面距離y+式(26)有

由式(23)、式(34)、式(39)、式(40),邊界層內層湍流黏性式(37)可寫為

經驗參數β,m,C1的選擇在2.3 節討論.

對于式(36)的Richardson 數,考慮到ξ ≈1,有

對所研究的問題Vm=W,所以由式(38)可以得到邊界層外層的黏性系數

由式(29)給出Klebanoff間歇系數Fkl.

由于研究對象為旋轉流動,定義位移厚度應使用相對速度W-Vφ,式(25)可寫為

動量厚度

2.2 外層湍流黏性公式的修正

本文研究的旋轉圓管壁面邊界層與平板邊界層有區別,其軸對稱性使計算需特殊處理.對于平板邊界層,邊界層從平板前緣點開始,向下游(縱向)發展.下游各橫截面上邊界層厚度、位移厚度和動量厚度由上一橫截面上信息計算得到,這是因為平板邊界層沿縱向是拋物型偏微分方程.對于本文中的旋轉軸對稱邊界層,由于初值未知,因此需要使用迭代方法計算,即先設置初始猜想值,通過迭代收斂求解邊界層.

下面進行迭代求解的說明.圖1 示意了靜止壁面和旋轉壁面邊界層速度分布.實線下的面積記為S1,S2,虛線下的面積記為如果黏性增加,則將導致層與層之間的動量交換.對于靜止的平板邊界層:.

圖1 靜止壁面(左)和旋轉壁面(右)邊界層:1–νt 小;2–νt 大Fig.1 Typical velocity profiles in the boundary layers on a stationary surface(Left)and a rotating cylinder surface(Right):1–less νt values;2–greater νt values

而對于旋轉邊界層

邊界層位移厚度的計算依賴于式(12)所描述的uφ速度分布.由式(19)=us可知無量綱位移厚度主要依賴于us和Reφ.對R=0.11 m 以空氣為介質和R=0.055 m 以水為介質的情況,在適用于低Reynolds數的計算網格上使用Reynolds 應力湍流模型進行數值仿真,得到不同us,Reφ組合的邊界層厚度(見表2),基于此建立(us,Reφ)的關聯公式.

對表2 中的數據分析擬合,得到無量綱位移厚度

表2 旋轉吸氣圓筒表面無量綱邊界層厚度Table 2 Values of the dimensionless boundary layer thicknesses on the surface of a rotating suction cylinder

無量綱厚度

通過上兩式得到相應邊界層厚度的初始估計值.

式(49)～式(52)作為軸對稱旋轉邊界層迭代求解的初值.

2.3 流線曲率和壁面吸氣所對應的經驗參數選取

由流線曲率以及相應的離心力對湍流影響由Bradshaw 修正式(35) 得到,結合到式(41) 和式(44)式中進行內層和外層湍流黏性修正.對于經驗參數m,取m=1.β 依賴壁面曲率、吸氣率和旋轉速度.從黏性底層到湍流強烈的外層,β 的變化與形狀因子有關,詳細的流場數值模擬表明如下的β 取值較合理

吸氣對邊界層的影響通過式(41)中的經驗參數A*修正.A*指數形式的式(42)中系數C1在平板邊界層的CS 模型中取C1=11.8.對于本文研究的旋轉邊界層,C1與邊界層中的湍流程度有關,而與湍流程度密切相關.對于層流邊界層,=2.湍流黏性為=0,從式(41)和式(42)可知C1→∞;對于本文研究的已經發展的湍流邊界層≈1.17,C1≈11.通過詳細數值模擬的流場數據的校核標定,表明C1的較好近似

這與上述兩個極限情況吻合良好.

注意到系數β 和C1校核與形狀因子的確定是耦合的:一方面β 和C1依賴于形狀因子的值,而反過來β 和C1值的選擇又影響形狀因子的值.基于已驗證詳細的數值模擬結果,計算形狀因子值,建立標定所需數據庫.對于不同Reφ和us的組合,使用所建立的近似方法進行不同β 和C1參數組合的試算,求得相應的值.然后根據梯度下降法改變β 和C1取值,使詳細數值模擬得到的與基于代數湍流模型方法得到的之差在容許誤差之內.有了β 和C1離散形式參數值集后,對其進行擬合,得到最小二乘擬合函數β()、C1(),其結果見式(53)和式(54).

尾跡參數Π0使用經典公式(31),它對壁面曲率和吸氣的依賴通過合理選擇參數β 和C1體現.

3 算法及驗證

3.1 迭代計算過程

邊界層流動方程(20)的近似解以及上述定義代數湍流模型的方程,表明可以用迭代法求解吸氣旋轉圓筒壁面的邊界層流動.主要步驟如下.

(1)確定下列參數的初始值:

由式(20)旋轉Reynolds 數Reφ;無量綱吸氣速度us;徑向空間步長kr;徑向最大的無量綱尺寸=ξmax;無量綱位移厚度計算迭代容差ε*.

(2)第一步計算:

(a)從式(19)得到無量綱摩擦速度uw;

(c) 由式(51) 和式(52) 預測無量綱邊界層厚度和動量厚度和;

(8)下面幾步為邊界層外層計算循環:

(a)當前步計算初始湍流黏性由式(30)、式(31)、式(43)、式(45)、式(46)、式(49)得到;

(c)當前步結束時計算導數?uφ/?ξ:

下一步長

(d) 當達到區域邊界時,計算結束,否則回到(a)繼續計算.

3.2 迭代計算方法的驗證

為了驗證上述方法的正確性,對比本文的代數湍流模型和詳細流場數值模擬得到的邊界層周向速度,其結果如圖2 所示.圖中也包括由式(17)描述的層流邊界層.Reφ和Rer的組合覆蓋層流邊界層和離心不穩定發生前的湍流邊界層.圖2 中的Reφ和Rer組合并不是用于校準湍流模型參數的Reφ和Rer組合.

從圖2 可以看出,本文的代數湍流模型所預測的旋轉圓筒周向速度與Reynolds 應力微分形式湍流模型的詳細數值模擬得到的結果幾乎一致.并且它也能預測強吸氣條件下的層流邊界層速度分布,見圖2(f).

圖3 顯示了邊界層內的湍流黏性分布,可以看出代數湍流模型的預測結果與Reynolds 應力模型的數值計算結果比較接近.其中圖3(a)為層流邊界層,理論上湍流黏性為0,本文的代數湍流模型的湍流黏性預測值也非常小.

以上流場結果也證實,對于強吸氣旋轉圓筒Couette 流動,其主要流動為旋轉圓筒壁附近的邊界層流動.隨著旋轉Reynolds 數Reφ的增大,湍流及離心不穩定性增強.吸氣量減弱(即Rer減小),邊界層增厚,邊界層內湍流強度增大.在大Rer情況下,離心不穩定發生大為延后.因此,可以通過合適的Reφ和Rer組合來控制邊界層的厚度和湍流強度,具有很好的創新實用價值.

圖2 旋轉圓筒壁面邊界層的平均速度分布:1–Reynolds 應力湍流模型的數值仿真結果;2–本文代數湍流模型結果;3–層流邊界層,即式(17)Fig.2 Mean velocity profiles in the boundary layer on the surface of a rotating suction cylinder:1–detailed numerical simulation,RSM turbulence model;2–application of the algebraic turbulence model;3–laminar profile,equation(17)

圖3 湍流黏性:1–Reynolds 應力湍流模型的數值仿真結果;2–本文代數湍流模型結果Fig.3 Turbulent viscosity ratio profiles:1–detailed numerical simulation,RSM turbulence model;2–application of the algebraic turbulence model

圖3 湍流黏性:1–Reynolds 應力湍流模型的數值仿真結果;2–本文代數湍流模型結果(續)Fig.3 Turbulent viscosity ratio profiles:1–detailed numerical simulation,RSM turbulence model;2–application of the algebraic turbulence model(continued)

4 結論

本文針對吸氣旋轉壁面的湍流邊界層流動建立了代數湍流模型.基于Cebeci-Smith 兩層代數湍流模型,通過推導分析,對其進行修正和經驗參數的調整,用以考慮離心力場(流線曲率) 和壁面吸氣等因素的影響.利用Reynolds 應力湍流模型的詳細數值模擬得到的流場,校準Cebeci-Smith 代數模型的經驗參數.

建立了壁面吸氣條件下軸對稱旋轉壁面邊界層的迭代算法,它容易用計算機程序實現.該迭代算法能用于離心失穩(Taylor 渦)發生前的層流及湍流邊界層計算.在旋轉多孔圓筒強吸氣條件下,離心不穩定在較大的旋轉速度(Reφ～105) 條件下發生,通過調整吸氣強度可以控制邊界層厚度和湍流強度.這為旋轉壁面的剪切力、動量傳遞、傳熱和傳質的控制提供了新思路.本研究可應用于動態旋轉過濾、轉子軸的熱保護、新型旋轉熱交換器等領域.