同心圓柱體間旋轉(zhuǎn)氣隙內(nèi)對(duì)流換熱的數(shù)值模擬研究

孫玉昕，趙旭峰，吳 波

（1 武漢工程大學(xué) 計(jì)算機(jī)科學(xué)與工程學(xué)院，湖北 武漢 430074；2 華中科技大學(xué) 能源與動(dòng)力工程學(xué)院，湖北 武漢 430074；3 中國(guó)航空動(dòng)力機(jī)械研究所，湖南 株洲 412000）

0 前言

發(fā)電機(jī)內(nèi)部溫度場(chǎng)的分布是保證發(fā)電機(jī)安全穩(wěn)定持續(xù)運(yùn)行的關(guān)鍵因素。如果溫升過高，會(huì)導(dǎo)致絕緣的損壞，使絕緣的電氣性能下降并引發(fā)電機(jī)內(nèi)部的各種放電、短路故障，造成汽輪發(fā)電機(jī)結(jié)構(gòu)上的損壞，最終導(dǎo)致汽輪發(fā)電機(jī)的故障，直接影響電站的可靠、經(jīng)濟(jì)運(yùn)行[1]。定子和轉(zhuǎn)子間環(huán)形氣隙作為大型汽輪發(fā)電機(jī)運(yùn)行過程中通風(fēng)換熱的重要組成部分，其內(nèi)部流場(chǎng)的流動(dòng)換熱情況需要進(jìn)行深入研究。

轉(zhuǎn)子外壁和定子內(nèi)壁分別形成了環(huán)形氣隙的內(nèi)、外壁面。由于發(fā)電機(jī)定子線圈和轉(zhuǎn)子線圈的存在，環(huán)形氣隙的內(nèi)外壁表面并不是光滑的，而是形成了與線圈位置相對(duì)應(yīng)的軸向凹槽。壁面上凹槽的存在使得環(huán)形氣隙中的流動(dòng)換熱情況異常復(fù)雜。

流體在旋轉(zhuǎn)同軸圓筒間的流動(dòng)稱為泰勒-庫(kù)特流。泰勒-庫(kù)特流因其獨(dú)特的流動(dòng)特性得到了越來越廣泛的應(yīng)用。例如，可以利用湍流泰勒渦流抑制膜污染和濃差極化現(xiàn)象來強(qiáng)化旋轉(zhuǎn)管式微濾膜的分離過程；也可利用其徑向混合性能來研制均混反應(yīng)器；還可以深入認(rèn)識(shí)高速機(jī)械的潤(rùn)滑液流動(dòng)狀態(tài)。除此之外，泰勒-庫(kù)特流的徑向流動(dòng)特性可以強(qiáng)化流場(chǎng)的傳熱傳質(zhì)過程,減小壁面附近流體與主流體之間的溫度和濃度差異[2]。因此，為研究汽輪發(fā)電機(jī)內(nèi)環(huán)形氣隙的流動(dòng)換熱特性，可以首先對(duì)定子、轉(zhuǎn)子間環(huán)形氣隙內(nèi)流動(dòng)的簡(jiǎn)化模型—泰勒-庫(kù)特流模型進(jìn)行研究，通過數(shù)值模擬的方法研究其流動(dòng)換熱規(guī)律，了解流場(chǎng)的分布對(duì)流動(dòng)換熱的影響，以增進(jìn)對(duì)汽輪發(fā)電機(jī)冷卻通流部分的了解，為其通流結(jié)構(gòu)的優(yōu)化設(shè)計(jì)和可靠運(yùn)行提供依據(jù)。

國(guó)內(nèi)外很多學(xué)者都對(duì)泰勒-庫(kù)特流進(jìn)行過研究。Taylor[3]在研究?jī)瑟?dú)立同心旋轉(zhuǎn)圓柱體間環(huán)形間隙內(nèi)流動(dòng)特性的過程中，用理論預(yù)測(cè)和實(shí)驗(yàn)驗(yàn)證了當(dāng)內(nèi)壁面轉(zhuǎn)速達(dá)到臨界轉(zhuǎn)速時(shí)流場(chǎng)失穩(wěn)，流體形成螺旋軌跡運(yùn)動(dòng)的現(xiàn)象。圓柱壁面轉(zhuǎn)速較低時(shí)，流體沿壁面轉(zhuǎn)動(dòng)方向做簡(jiǎn)單的移動(dòng)；當(dāng)轉(zhuǎn)速不斷增大，達(dá)到臨界轉(zhuǎn)速后，流場(chǎng)失穩(wěn)，流體就會(huì)形成細(xì)胞狀微團(tuán)并在流場(chǎng)中作螺旋軌跡運(yùn)動(dòng)，這種流動(dòng)就稱為泰勒-庫(kù)特流，流場(chǎng)中形成的細(xì)胞狀微團(tuán)就是眾所周知的泰勒渦。Wereley等[4]通過試驗(yàn)的方法研究了內(nèi)圓柱轉(zhuǎn)動(dòng)外圓筒靜止的環(huán)形間隙內(nèi)泰勒庫(kù)特流，Hwang等[5]通過數(shù)值計(jì)算的方法對(duì)文獻(xiàn)[4]中的試驗(yàn)研究情況進(jìn)行了模擬驗(yàn)證。R.Jakoby等[6]和I.Cornet等[7]通過試驗(yàn)研究了旋轉(zhuǎn)圓筒泰勒渦的形成對(duì)內(nèi)壁面換熱的影響，并且得出了平均努賽爾數(shù)與旋轉(zhuǎn)雷諾數(shù)的對(duì)應(yīng)關(guān)系。國(guó)內(nèi)有一些學(xué)者也對(duì)于泰勒-庫(kù)特流的流動(dòng)特性進(jìn)行了研究，文獻(xiàn)[2][8][9]中都是采用數(shù)值模擬的方法對(duì)流動(dòng)特性進(jìn)行的模擬和分析。

本文運(yùn)用商用CFD軟件對(duì)同心圓柱體間環(huán)形間隙內(nèi)流體對(duì)流換熱情況進(jìn)行數(shù)值模擬，對(duì)流場(chǎng)中速度分布、溫度分布、壁面熱流密度分布以及不同泰勒數(shù)下流場(chǎng)的對(duì)流換熱特性進(jìn)行了討論和分析，并對(duì)帶凹槽環(huán)形氣隙內(nèi)的流動(dòng)進(jìn)行了模擬研究。

1 數(shù)值模型

1.1 物理模型及計(jì)算區(qū)域

幾何模型的簡(jiǎn)圖如1所示。其中，內(nèi)圓柱面半徑ri=43.4mm，外圓柱面半徑ro=52.3mm，間隙高度d=ro-ri=8.9mm，內(nèi)外半徑比η=ri/r0=0.83，模型軸向長(zhǎng)度L=356mm，軸向長(zhǎng)度與間隙高度比Γ=L/d=40。

網(wǎng)格模型簡(jiǎn)圖如圖2所示。整個(gè)計(jì)算域采用結(jié)構(gòu)化網(wǎng)格，并對(duì)局部進(jìn)行加密。網(wǎng)格總數(shù)量為1,966,080，其中，徑向、軸向和周向網(wǎng)格數(shù)量分別為 60、256、128。內(nèi)外壁面均設(shè)置邊界層網(wǎng)格，邊界層網(wǎng)格均為10層，起始高度為0.005mm，增長(zhǎng)比例為1.2。

圖1 幾何模型

1.2 控制方程

應(yīng)用Fluent6.3軟件，采用三維不可壓、穩(wěn)態(tài)、層流模型。控制方程[10]：

（1）連續(xù)性方程：

（2）動(dòng)量方程：

（3）能量方程：

對(duì)于湍流問題，還需要引入湍流模型，本文采用k-ε兩方程模型[10]

（4）k-ε兩方程模型：

2 “泰勒-庫(kù)特流”流動(dòng)換熱特性的求解和分析

2.1 物理?xiàng)l件及計(jì)算方法

計(jì)算模型內(nèi)壁面溫度設(shè)為343K，其轉(zhuǎn)速設(shè)置詳見表1，外壁面溫度設(shè)為313K。介質(zhì)為40℃空氣（密度：1.128kg/m3(采用boussinesq假設(shè))；比熱容：1005 J/kg·K；熱導(dǎo)率：0.0276w/m·K；運(yùn)動(dòng)粘性系數(shù)：1.913×10-5m2/s；熱膨脹率(1/K)：0.00319；），沒有軸向流動(dòng)。

模型內(nèi)外壁面間隙內(nèi)的空氣，在常物性和宏觀能量守恒的假設(shè)條件下，其流動(dòng)和能量傳遞滿足連續(xù)性方程、動(dòng)量方程和能量方程。速度壓力耦合方法采用PISO算法，動(dòng)量方程、能量方程以及湍流工況下的湍流耗散方程均采用二階迎風(fēng)格式。

表1 工況設(shè)置

2.2 結(jié)果與分析

模擬計(jì)算按泰勒數(shù)Ta不同設(shè)置5個(gè)工況。其中，前四個(gè)工況中流場(chǎng)為層流狀態(tài)，隨著Ta的增大，到第5個(gè)工況時(shí)，流場(chǎng)為湍流狀態(tài)。模型在工況1和工況5下得到的子午面上溫度場(chǎng)分布云圖如圖3所示。

圖3 不同工況下模型子午面溫度云圖

從子午面溫度分布云圖可以看出，由于模型內(nèi)壁面轉(zhuǎn)動(dòng)在流場(chǎng)中形成的泰勒渦的作用，溫度場(chǎng)也形成了類似“泰勒渦”的成對(duì)的渦團(tuán)狀溫度分布。

以工況2（Ta=20000）的模擬結(jié)果為例，在模型中取一條沿軸向方向并經(jīng)過內(nèi)、外壁面間氣隙中心的采樣直線，取該直線上的速度與溫度數(shù)據(jù)，以及內(nèi)壁面上的熱流密度沿軸向的分布數(shù)據(jù)。對(duì)比該采樣直線上速度與溫度沿軸向的分布如圖。由圖可以看出，速度、溫度沿軸向的分布都呈現(xiàn)周期性的上下波動(dòng)，波動(dòng)頻率與流場(chǎng)中泰勒渦的數(shù)量相同。同時(shí)，溫度的波動(dòng)與速度的波動(dòng)趨勢(shì)一致。沿模型軸向方向，每對(duì)渦的中心位置對(duì)應(yīng)圖中速度分布曲線的“波峰”位置，同時(shí)也對(duì)應(yīng)溫度分布曲線的“波峰”位置；相鄰兩渦的交界的位置對(duì)應(yīng)速度分布曲線的“波谷”位置，同時(shí)也對(duì)應(yīng)溫度分布曲線的“波谷”位置。這說明，在每對(duì)渦中，中心位置處的速度最大，溫度最高；在相鄰兩渦的交界的位置，速度最小，溫度也最低。

對(duì)比采樣直線上速度分布與內(nèi)壁面熱流密度沿軸向的分布如圖 4。由圖可以看出，內(nèi)壁面熱流密度沿軸向的分布也呈現(xiàn)周期性的上下波動(dòng)，且波動(dòng)頻率與流場(chǎng)中泰勒渦的數(shù)量相同。不同的是，沿模型軸向方向熱流密度的波動(dòng)與速度波動(dòng)的趨勢(shì)完全相反。速度最大的位置，熱流密度最??；速度最小的位置，熱流密度最大。這是因?yàn)樗俣却蟮奈恢闷錅囟纫蚕鄬?duì)較高，流體與高溫內(nèi)壁面間溫差較小，溫度梯度較小，導(dǎo)致熱流密度較?。幌喾?，速度小的位置溫度較低，流體與高溫內(nèi)壁面間溫差較大，溫度梯度也相對(duì)較大，因而使熱流密度較大。綜上所述，對(duì)于整個(gè)流場(chǎng)，其局部的換熱強(qiáng)度沿軸向也呈現(xiàn)周期性的變化，在每對(duì)渦的中心處最小，沿軸向向兩側(cè)逐漸增加，到相鄰泰勒渦的交界處達(dá)到最大，然后再逐漸減小，到渦的中心處達(dá)到最小。如同速度、溫度的分布一樣，往復(fù)交替變化。

圖4 Ta=20000的工況下速度、熱流密度、溫度軸向分布

由圖5可以看出，各工況下局部對(duì)流換熱系數(shù)的分布曲線與速度、溫度、熱流密度的分布曲線一樣，沿軸向周期性波動(dòng)，其波動(dòng)趨勢(shì)與速度分布和溫度分布相反，與熱流密度的分布相同，即在每對(duì)渦的中心處最小，在相鄰兩對(duì)渦的交界處最大。同時(shí)，從總體上看，泰勒數(shù)Ta越大（內(nèi)壁面轉(zhuǎn)速越高），平均對(duì)流換熱系數(shù)就越大。從圖中也可以看出，隨著泰勒數(shù)Ta的增加，平均努賽爾數(shù)Nu也不斷增加。也就是說，模型內(nèi)壁面的轉(zhuǎn)速越高，壁面與流體的對(duì)流換熱強(qiáng)度越高，流體對(duì)壁面的冷卻效果越好。

圖5 各工況下局部對(duì)流換熱系數(shù)沿軸向的分布曲線

圖6 數(shù)值計(jì)算Nu數(shù)與經(jīng)驗(yàn)公式結(jié)果對(duì)比

I.S. Bjorklund和W.M. Kays[11]對(duì)環(huán)形氣隙內(nèi)流體在內(nèi)壁面旋轉(zhuǎn)情況下的對(duì)流換熱特性進(jìn)行了實(shí)驗(yàn)研究，并總結(jié)出計(jì)算該模型努賽爾數(shù)的經(jīng)驗(yàn)公式：

本文所選數(shù)值計(jì)算模型的尺寸及泰勒數(shù)的設(shè)置均滿足該經(jīng)驗(yàn)公式的使用范圍，將數(shù)值計(jì)算得出的努賽爾數(shù)與經(jīng)驗(yàn)公式得出的努賽爾數(shù)對(duì)比，如圖6。

由對(duì)比結(jié)果可以看出，數(shù)值計(jì)算的結(jié)果與文獻(xiàn)中經(jīng)驗(yàn)公式的最大偏差為9.88%，最小偏差為0.46%，平均偏差為4.50%，并且，模擬結(jié)果和經(jīng)驗(yàn)公式得出的努賽爾數(shù)的變化趨勢(shì)一致。這一結(jié)果證明了該數(shù)值計(jì)算的有效性。

3 帶軸向凹槽氣隙模型的流動(dòng)模擬

為使環(huán)形氣隙模型更加接近真實(shí)情況，參照Alstom公司制造的某型號(hào)發(fā)電機(jī)的結(jié)構(gòu)特點(diǎn)，建立帶軸向凹槽氣隙模型，對(duì)氣隙內(nèi)流動(dòng)情況進(jìn)行模擬。為減少計(jì)算時(shí)間，提高計(jì)算效率，根據(jù)環(huán)形氣隙周期性對(duì)稱的結(jié)構(gòu)特點(diǎn)，取整個(gè)環(huán)形氣隙的1/56作為計(jì)算域以減少網(wǎng)格數(shù)量，計(jì)算域網(wǎng)格分布如圖7所示，整個(gè)計(jì)算域同樣采用結(jié)構(gòu)化網(wǎng)格，并對(duì)局部進(jìn)行加密，網(wǎng)格總數(shù)量為 2,662,214。模型內(nèi)壁面轉(zhuǎn)速為 377rad/s（3600r/min），介質(zhì)采用理想氣體。

圖7 帶軸向凹槽氣隙模型網(wǎng)格分布

圖8 帶軸向凹槽氣隙模型子午面速度云圖

計(jì)算得到模型子午面的速度分布如圖8所示?？梢钥闯?，由于環(huán)形氣隙內(nèi)外壁面上凹槽的存在，極大的增加了壁面對(duì)流場(chǎng)的擾動(dòng)作用，氣隙內(nèi)的流動(dòng)更加混亂，使得流場(chǎng)中渦的數(shù)量減少，尺度增大。由光滑壁面模型的流動(dòng)換熱計(jì)算可知，流場(chǎng)中溫度和熱流密度的分布受速度分布的影響很大，可以推測(cè)，壁面凹槽在影響速度分布的同時(shí)，也影響了溫度以及熱流密度的分布。

分別在距離內(nèi)壁面 5mm（R=0.58m）和氣隙中部（R=0.6065m）各取一條沿軸向方向的采樣直線，提取采樣直線上的徑向分速度數(shù)據(jù)，如圖9所示。由該圖可以看出，沿模型軸向方向，流場(chǎng)的徑向分速度呈正負(fù)交替變化的趨勢(shì)。在內(nèi)壁面附近（R=0.58），正向速度的峰值略大于負(fù)向速度的峰值，但二者差距不大；在氣隙中部主流區(qū)域（R=0.6065m），正向速度的峰值要明顯大于負(fù)向速度的峰值。而泰勒渦的形成也正是由于流場(chǎng)中的這種方向相反的徑向流動(dòng)間隔作用的結(jié)果。由于流場(chǎng)中的這種較強(qiáng)的徑向流動(dòng)，使得內(nèi)壁面附近的流體能夠更多的流到外壁面附近，從而強(qiáng)化了流場(chǎng)的傳熱傳質(zhì)過程。

圖9 帶軸向凹槽氣隙模型徑向速度的分布

4 結(jié)論

（1）內(nèi)壁面的轉(zhuǎn)動(dòng)在流場(chǎng)中形成的泰勒渦改變了流體內(nèi)溫度場(chǎng)的分布，使溫度場(chǎng)也形成了類似“泰勒渦”的成對(duì)的渦團(tuán)狀分布；

（2）溫度場(chǎng)沿軸向的周期性波動(dòng)，使內(nèi)壁面的熱流密度沿軸向方向也產(chǎn)生了周期性波動(dòng)，在渦的中心位置，熱流密度最小，對(duì)流換熱系數(shù)最小，該位置對(duì)流換熱的能力最差；在相鄰兩對(duì)渦的交界處，熱流密度最大，對(duì)流換熱系數(shù)最大，該位置對(duì)流換熱的能力最好；

（3）內(nèi)壁面轉(zhuǎn)速越大，平均努賽爾數(shù)Nu越大，流場(chǎng)的對(duì)流換熱性能越好，流體對(duì)壁面的冷卻效果就越好；

（4）環(huán)形氣隙壁面增加凹槽后，壁面對(duì)流場(chǎng)擾動(dòng)作用的增強(qiáng)使得流場(chǎng)中速度的分布發(fā)生了顯著的變化；可以推測(cè)，這種擾動(dòng)也會(huì)使溫度和熱流密度的分布發(fā)生變化，凹槽的存在對(duì)流場(chǎng)換熱性能的影響會(huì)在后續(xù)的工作中進(jìn)行深入研究；

（5）氣隙內(nèi)的泰勒-庫(kù)特流在流場(chǎng)形成了明顯徑向流動(dòng)，這種徑向流動(dòng)強(qiáng)化了發(fā)電機(jī)環(huán)形氣息內(nèi)的傳熱傳質(zhì)過程；

由于大型汽輪發(fā)電機(jī)定子內(nèi)壁面和轉(zhuǎn)子外壁面結(jié)構(gòu)特殊，內(nèi)外壁面間氣隙內(nèi)真實(shí)的流動(dòng)換熱情況更加復(fù)雜，本文的工作僅是對(duì)該問題從數(shù)值模擬的角度開展前期的基礎(chǔ)研究，為了解其實(shí)際的運(yùn)行狀況，提高發(fā)電機(jī)運(yùn)行的安全性和穩(wěn)定性，還需要進(jìn)一步的研究和探索。

[1]丁樹業(yè). 大型發(fā)電機(jī)定子復(fù)雜結(jié)構(gòu)內(nèi)流體流動(dòng)與傳熱特性的研究[D]. 哈爾濱：哈爾濱理工大學(xué)，2008.6.

[2]周先桃，潘家禎，陳理清，等. 湍流泰勒渦特性的數(shù)值模擬[J]. 華東理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 2006， 32（5）:617-622.

[3]Taylor GI. Stability of a viscous liquid contained between two rotating cylinders[J]. Phil Trans R Soc, 1923, A223:289-343.

[4]Steven T. Wereley, Richard M. Lueptow, Velocity field for Taylor-Couette flow with an axial flow[J]. J. PHYSICS OF FLUIDS, 11(12).

[5]Jong-Yeon Hwang, Kyung-Soo Yang, Numerical study of Taylor-Couette flow with an axial flow[J]. J. Computers ＆ Fluids,2004,33: 97-118.

[6]R. Jakoby, S. Kim, S. Wittig. Correlations of the convection heat transfer in annular channels with rotating inner cylinder[J].ASME J. Eng. Gas Turbines Power, 1991, 21: 670-677.

[7]I. Cornet, R. Greif, J.T. Teng, F. Roehler, Mass transfer to rotating robs and plates[J]. International Journal of Heat Mass Transfer ,1980,23: 805-811.

[8]王宇，徐忠鋒. 同軸圓柱間旋轉(zhuǎn)流動(dòng)Taylor-Couette流的數(shù)值模擬[J]. 渤海大學(xué)學(xué)報(bào)(自然科學(xué)版)，2004, 25(3):249-251.

[9]劉難生，董宇紅，陸夕云，等. 旋轉(zhuǎn)同心圓筒間 Couette-Taylor流動(dòng)的數(shù)值模擬[J]. 中國(guó)科學(xué)技術(shù)大學(xué)學(xué)報(bào)，2002,32(1):91-97.

[10]陶文銓. 數(shù)值傳熱學(xué)[M]. 西安: 西安交通大學(xué)出版社, 2001.

[11]I.S. Bjorklund, W.M. Kays, Heat transfer between concentric rotating cylinders[J]. Journal of Heat Transfer, 1959, 81: 175-186.