周陳龍，丁保庚，王 穎，李興華

（核工業(yè)理化工程研究院，天津 300180）

長期以來，國內(nèi)外在研究高速旋轉(zhuǎn)流場時，一直以等溫剛體旋轉(zhuǎn)狀態(tài)為基礎(chǔ)，對納維－斯托克斯方程（N－S方程）進(jìn)行線形化求解。大量的研究已證明，旋轉(zhuǎn)圓柱筒內(nèi)流場不會發(fā)生不穩(wěn)定情況[1]，但現(xiàn)在考慮的是高速短圓柱筒轉(zhuǎn)速很大，其內(nèi)部流場的最高馬赫數(shù)達(dá)到了7左右，這時流場平衡時的狀態(tài)，是一個非常重要而在一些工程中又急需解決的問題。

關(guān)于高速旋轉(zhuǎn)流場的試驗(yàn)研究一直是我們的一個薄弱環(huán)節(jié)。由于圓柱筒高速旋轉(zhuǎn)，其內(nèi)部流場徑向跨越了分子流、過渡流、粘性流區(qū)域，目前還未找到一種合適的、不干擾流場的、具有高精度的測量高速強(qiáng)旋轉(zhuǎn)流場的方法。因此，采用數(shù)值解法進(jìn)行模擬顯得十分必要。

本工作利用數(shù)值分析法，采用二階迎風(fēng)格式、總能守恒模型，對封閉狀態(tài)下高速旋轉(zhuǎn)的短圓柱筒內(nèi)的三維流場進(jìn)行數(shù)值研究，并與理論計(jì)算結(jié)果進(jìn)行對比分析。

1 等溫剛體旋轉(zhuǎn)狀態(tài)下的平衡解[2]

在等溫剛體旋轉(zhuǎn)狀態(tài)下，圓柱筒內(nèi)氣體角速度與圓柱筒角速度相等，溫度是均勻的，都等于初始狀態(tài)溫度T0，故有：

式中：vr為氣體速度徑向分量；vθ為氣體速度角向分量；Ω為圓柱筒角速度；r為半徑；vz為氣體速度軸向分量；T為氣體溫度。

如果用p和ρ來表示這種狀態(tài)下的壓強(qiáng)和密度，給定的邊界條件是圓柱筒側(cè)壁處的壓強(qiáng)pw和密度ρw，則p和ρ的表達(dá)式為：

將相關(guān)參數(shù)代入式（1）和（3）中，得到等溫剛體旋轉(zhuǎn)狀態(tài)下的速度和壓強(qiáng)分布（圖1、2）。

圖1 等溫剛體速度分布Fig.1 Velocity distribution in isothermal rigid body rotational state

圖2 等溫剛體徑向壓強(qiáng)和－ln（p／pw）沿徑向分布Fig.2 Radial pressure distribution and－ln（p／pw）radial distribution in isothermal rigid body rotational state

2 物理模型

目前研究的短圓柱筒是高速旋轉(zhuǎn)設(shè)備，內(nèi)部絕大多數(shù)氣體均集中在靠近邊壁的區(qū)域。徑向上，流場跨越了分子流、過渡流、粘性流等流動區(qū)域。采用封閉模型，選取高度z與半徑r比為1∶2的短圓柱筒，得到的計(jì)算模型如圖3所示。

圖3 計(jì)算模型Fig.3 Computing model

3 數(shù)學(xué)模型

3.1 基本假設(shè)

基本假設(shè)如下：1）圓柱筒內(nèi)的流體假設(shè)全部為連續(xù)介質(zhì)，流體運(yùn)動遵循納維－斯托克斯方程；2）流體所受的重力與離心力相比可忽略；3）忽略由于輻射或其他原因傳到流體上的熱量；4）流體的粘性系數(shù)、熱傳導(dǎo)系數(shù)和等壓比熱均為常數(shù)；5）對于可壓縮流體的膨脹粘性系數(shù)采用斯托克斯假設(shè)，令膨脹粘性系數(shù)等于零；6）流體的狀態(tài)方程由完全氣體定律來描述。

3.2 控制方程

根據(jù)前面所作的基本假設(shè)，圓柱坐標(biāo)系下流體的運(yùn)動基本方程如下。

連續(xù)方程：

動量方程：

能量方程：

狀態(tài)方程：

其中：

式中：e為內(nèi)能；κ為熱傳導(dǎo)系數(shù)；μ為粘性系數(shù)；Φ 為粘性耗散函數(shù)；cp為等壓比熱容；i、j、k分別表示r、θ、z方向的單位矢量。

3.3 初始條件和邊界條件

計(jì)算采用圓柱坐標(biāo)系。初始條件采用等溫剛體旋轉(zhuǎn)條件：vr＝ 0m／s；vθ＝ Ωr；vz＝0m／s；p（r）＝ pwexp｛－ A2[1 － （r／ra）2]｝；T＝300K。

邊界條件：上、下邊界和側(cè)壁以vθ旋轉(zhuǎn)，T＝300K，使用壁面無滑移條件。

4 計(jì)算方法

4.1 計(jì)算網(wǎng)格

網(wǎng)格生成的質(zhì)量直接影響流場數(shù)值求解的精度和穩(wěn)定性，所生成的網(wǎng)格須滿足貼體性、光滑性、合理分布性及正交性的要求。

利用ICEM CFD軟件通過創(chuàng)建獨(dú)立于幾何模型的塊體結(jié)構(gòu)近似幾何模型，將塊拓?fù)渲械木W(wǎng)格自動映射到幾何模型上的方法，得到六面體單元。為了改善近壁局部的邊界層網(wǎng)格精度，采用了O－grid網(wǎng)格方法。由于短圓柱筒內(nèi)徑向壓強(qiáng)成指數(shù)分布，在靠近邊壁處進(jìn)行了高度加密處理，最后得到約11 708個單元（圖4、5）。

圖4 計(jì)算模型整體網(wǎng)格Fig.4 Whole mesh of computing model

圖5 網(wǎng)格剖面圖Fig.5 Section plan of mesh

4.2 計(jì)算方法簡介

本工作采用ANSYS CFX軟件進(jìn)行計(jì)算求解，該軟件采用基于有限元的有限體積法。采用由瞬態(tài)計(jì)算得到穩(wěn)態(tài)解的方法，計(jì)算坐標(biāo)系為柱坐標(biāo)系，按總能守衡進(jìn)行求解，計(jì)算采用的模型為剪切應(yīng)力輸運(yùn)（SST）k－ω模型。

SSTk－ω模型是二方程模型中的一種。與標(biāo)準(zhǔn)二方程k－ω 模型相比，SST k－ω 模型中增加了橫向耗散導(dǎo)數(shù)項(xiàng)，在湍流粘度定義中考慮了湍流剪切應(yīng)力輸運(yùn)過程。因此，該模型提高了強(qiáng)逆壓梯度或強(qiáng)順壓梯度的非平衡湍流邊界層及分離流計(jì)算能力[3]，適用于本工作的計(jì)算。

5 計(jì)算結(jié)果及分析

工作介質(zhì)的相關(guān)參數(shù)采用參考壓強(qiáng)為0.1MPa、溫度為300K時的值。計(jì)算采用4種方案：大時間步長（情況1）、小時間步長（情況2）、計(jì)算域旋轉(zhuǎn)和小時間步長（情況3）、高密度網(wǎng)格下的小時間步長（情況4）。所有計(jì)算方案均在3.3節(jié)的計(jì)算條件下進(jìn)行。

5.1 大時間步長計(jì)算

計(jì)算采用瞬態(tài)計(jì)算的方法，初始步長為1×10－7s，推進(jìn)幾百步后，時間步長改為1×10－4s，收斂效果很好。圖6示出縱切面壓強(qiáng)等值線、軸向中心橫截面徑向壓強(qiáng)、縱切面速度、軸向中心橫截面處徑向速度分布。從圖6a可看出，壓強(qiáng)呈軸對稱分布，從中心區(qū)域沿徑向逐漸增大，在外邊壁處達(dá)到最大。壓強(qiáng)分布近似與軸向坐標(biāo)z無關(guān)，只是徑向坐標(biāo)r的函數(shù)。從圖6b可看出，在0.5z的中心橫切面上，壓強(qiáng)在徑向上的分布與剛體旋轉(zhuǎn)情況下的分布存在很大差異，在中心處壓強(qiáng)約為760Pa，遠(yuǎn)大于剛體假設(shè)時的值，且并未體現(xiàn)出指數(shù)增長的規(guī)律，徑向壓強(qiáng)梯度較小，壓強(qiáng)在邊壁處約為3 260Pa，與剛體旋轉(zhuǎn)時的15kPa的邊壁壓強(qiáng)相差很大。

從圖6c可看出，速度呈軸對稱分布，中心處速度最?。?m／s），沿徑向逐漸增大。在圓柱筒中部速度分布與軸向坐標(biāo)基本無關(guān)，只是徑向坐標(biāo)的函數(shù)，但在靠近上、下邊界處，速度明顯增大。從圖6d可看出，速度在中心區(qū)域基本呈線性分布，當(dāng)r＞0.8ra時，速度急劇增大，在靠近邊壁的區(qū)域內(nèi)表現(xiàn)出很大的速度梯度。速度沿徑向的分布與剛體旋轉(zhuǎn)的速度分布（圖1）存在很大差異。

上述結(jié)果確實(shí)是方程的解，總計(jì)算時間相當(dāng)于圓柱筒旋轉(zhuǎn)約3 200圈，解很穩(wěn)定。但試驗(yàn)數(shù)據(jù)證明這是一非物理解。研究尋找其物理解的方法將是很重要的任務(wù)，具有重要學(xué)術(shù)意義和實(shí)用意義。

5.2 小時間步長的計(jì)算

計(jì)算條件與5.1節(jié)完全相同，但采用小步長1×10－8s，觀察用小步長能否得到物理解。

軸向中心橫截面處圓周速度沿徑向分布和徑向壓強(qiáng)分布示于圖7。由圖7a可看出，速度仍大幅偏離剛體速度分布。由圖7b可看出，中心區(qū)域壓強(qiáng)大幅提高，這與圓周速度大幅低于剛體旋轉(zhuǎn)速度的現(xiàn)象吻合。由于僅計(jì)算10 000步，累計(jì)時間1×10－4s，未達(dá)到穩(wěn)態(tài)，這可由圖7a的曲線形狀反映出，但其趨勢已非常清楚，可以說，在這種情況下，采用小步長的方法仍不可能得到物理解。

5.3 計(jì)算域旋轉(zhuǎn)和小時間步長的計(jì)算

采用3.3節(jié)的計(jì)算條件，計(jì)算域以vθ旋轉(zhuǎn)，計(jì)算坐標(biāo)系為柱坐標(biāo)系，按總能守衡進(jìn)行求解，使用各向異性能考慮哥氏效應(yīng)的雷諾應(yīng)力湍流模型。采用小步長1×10－8s計(jì)算。

圓周速度、軸向速度、徑向速度等值線示于圖8。軸向中心橫截面處圓周速度沿徑向分布示于圖9a。由圖8、9a可看出，圓周速度基本上是剛體轉(zhuǎn)速。由圖8b、c可知，存在有軸向速度和徑向速度，但數(shù)值很小，每秒僅幾米，說明可近似認(rèn)為是剛體旋轉(zhuǎn)速度。

圖6 情況1的縱切面壓強(qiáng)等值線（a）、軸向中心橫截面徑向壓強(qiáng)（b）、縱切面速度（c）、軸向中心橫截面處徑向速度（d）分布Fig.6 Isobars at longitudinal section（a），radial pressure distribution at center transverse section（b），velocity distribution at longitudinal section（c）and radial velocity distribution at center transverse section（d）in case 1

圖7 情況2下軸向中心橫截面處圓周速度沿徑向（a）和徑向壓強(qiáng)（b）分布Fig.7 Circle velocity radial distribution（a）and radial pressure distribution（b）at center transverse section in case 2

軸向中心橫截面處徑向壓強(qiáng)分布和溫度等值線分別示于圖9b、10。由圖9b、10可見，中心區(qū)域?yàn)?00K，側(cè)壁附近達(dá)324K，不是等溫流動，而且，壓強(qiáng)沿徑向的變化與等溫剛體壓強(qiáng)分布（圖2a）也有很大差別。

圖8 情況3下的圓周速度（a）、軸向速度（b）、徑向速度（c）等值線Fig.8 Isoline of circle velocity（a），axial velocity（b）and radial velocity（c）in case 3

圖9 情況3下軸向中心橫截面處圓周速度沿徑向（a）和徑向壓強(qiáng)（b）分布Fig.9 Circle velocity radial distribution（a）and radial pressure distribution（b）at center transverse section in case 3

5.4 高密度網(wǎng)格下小時間步長的計(jì)算

對圖4的網(wǎng)格在圓周方向和徑向進(jìn)行加密，總的網(wǎng)格數(shù)約為241 600，比原來擴(kuò)大約22倍。用1×10－7s的時間步長計(jì)算，穩(wěn)態(tài)結(jié)果如圖11～13所示。

圖10 情況3下溫度等值線Fig.10 Isotherms of temperature in case 3

從圖11a可見，計(jì)算結(jié)果接近5.3節(jié)的情況。圓周速度基本是線性分布，存在軸向速度和徑向速度，但數(shù)值很小，可近似認(rèn)為是剛體旋轉(zhuǎn)速度。從圖11b可看出，壓強(qiáng)呈指數(shù)分布狀態(tài)，與理論計(jì)算結(jié)果相比，趨勢上基本相同。從數(shù)據(jù)上比較，在中心部位，數(shù)值計(jì)算得到的壓強(qiáng)是幾Pa，而非等溫剛體得到的10－6～10－7Pa，即從實(shí)驗(yàn)結(jié)果分析得出的中心區(qū)域壓強(qiáng)和數(shù)值計(jì)算得到的中心區(qū)域壓強(qiáng)基本吻合。由圖12、13可見，中心區(qū)域的溫度基本保持在300K左右，而邊壁附近的溫度在322K左右，不是等溫流動。

因此，采用高密度網(wǎng)格的小步長計(jì)算能得到物理解，與理論計(jì)算結(jié)果相比，數(shù)值計(jì)算結(jié)果更符合實(shí)際情況。

圖11 情況4下軸向中心橫截面處圓周速度沿徑向（a）和徑向壓強(qiáng)（b）分布Fig.11 Circle velocity radial distribution（a）and radial pressure distribution（b）at center transverse section in case 4

圖12 情況4下的溫度等值線Fig.12 Isotherms of temperature in case 4

圖13 情況4下軸向中心橫截面處徑向溫度分布Fig.13 Radial temperature distribution at center transverse section in case 4

6 結(jié)論

對本工作研究的封閉的短圓柱筒而言，當(dāng)圓柱筒轉(zhuǎn)速很大時，其內(nèi)部流場流動的最高馬赫數(shù)達(dá)到了7左右。研究表明，只有在一定的計(jì)算條件下，數(shù)值計(jì)算才能得到物理解。而在能得到物理解的情況下，高速旋轉(zhuǎn)圓柱筒內(nèi)流場平衡時不是等溫狀態(tài)，但速度和壓強(qiáng)分布基本保持剛體旋轉(zhuǎn)狀態(tài)。

[1]陳懋章.粘性流體動力學(xué)基礎(chǔ)[M].北京：高等教育出版社，2002.

[2]張存鎮(zhèn).離心分離理論[M].北京：原子能出版社，1987.

[3]朱自強(qiáng).應(yīng)用計(jì)算流體力學(xué)[M].北京：北京航空航天大學(xué)出版社，1998.

[4]杜慶華.工程力學(xué)手冊[M].北京：高等教育出版社，1994.