高速旋轉圓筒內流場的數值研究

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

（核工業理化工程研究院，天津 300180）

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

關于高速旋轉流場的試驗研究一直是我們的一個薄弱環節。由于圓柱筒高速旋轉，其內部流場徑向跨越了分子流、過渡流、粘性流區域，目前還未找到一種合適的、不干擾流場的、具有高精度的測量高速強旋轉流場的方法。因此，采用數值解法進行模擬顯得十分必要。

本工作利用數值分析法，采用二階迎風格式、總能守恒模型，對封閉狀態下高速旋轉的短圓柱筒內的三維流場進行數值研究，并與理論計算結果進行對比分析。

1 等溫剛體旋轉狀態下的平衡解[2]

在等溫剛體旋轉狀態下，圓柱筒內氣體角速度與圓柱筒角速度相等，溫度是均勻的，都等于初始狀態溫度T0，故有：

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

如果用p和ρ來表示這種狀態下的壓強和密度，給定的邊界條件是圓柱筒側壁處的壓強pw和密度ρw，則p和ρ的表達式為：

將相關參數代入式（1）和（3）中，得到等溫剛體旋轉狀態下的速度和壓強分布（圖1、2）。

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

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

2 物理模型

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

圖3 計算模型Fig.3 Computing model

3 數學模型

3.1 基本假設

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

3.2 控制方程

根據前面所作的基本假設，圓柱坐標系下流體的運動基本方程如下。

連續方程：

動量方程：

能量方程：

狀態方程：

其中：

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

3.3 初始條件和邊界條件

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

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

4 計算方法

4.1 計算網格

網格生成的質量直接影響流場數值求解的精度和穩定性，所生成的網格須滿足貼體性、光滑性、合理分布性及正交性的要求。

利用ICEM CFD軟件通過創建獨立于幾何模型的塊體結構近似幾何模型，將塊拓撲中的網格自動映射到幾何模型上的方法，得到六面體單元。為了改善近壁局部的邊界層網格精度，采用了O－grid網格方法。由于短圓柱筒內徑向壓強成指數分布，在靠近邊壁處進行了高度加密處理，最后得到約11 708個單元（圖4、5）。

圖4 計算模型整體網格Fig.4 Whole mesh of computing model

圖5 網格剖面圖Fig.5 Section plan of mesh

4.2 計算方法簡介

本工作采用ANSYS CFX軟件進行計算求解，該軟件采用基于有限元的有限體積法。采用由瞬態計算得到穩態解的方法，計算坐標系為柱坐標系，按總能守衡進行求解，計算采用的模型為剪切應力輸運（SST）k－ω模型。

SSTk－ω模型是二方程模型中的一種。與標準二方程k－ω 模型相比，SST k－ω 模型中增加了橫向耗散導數項，在湍流粘度定義中考慮了湍流剪切應力輸運過程。因此，該模型提高了強逆壓梯度或強順壓梯度的非平衡湍流邊界層及分離流計算能力[3]，適用于本工作的計算。

5 計算結果及分析

工作介質的相關參數采用參考壓強為0.1MPa、溫度為300K時的值。計算采用4種方案：大時間步長（情況1）、小時間步長（情況2）、計算域旋轉和小時間步長（情況3）、高密度網格下的小時間步長（情況4）。所有計算方案均在3.3節的計算條件下進行。

5.1 大時間步長計算

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

從圖6c可看出，速度呈軸對稱分布，中心處速度最小（為0m／s），沿徑向逐漸增大。在圓柱筒中部速度分布與軸向坐標基本無關，只是徑向坐標的函數，但在靠近上、下邊界處，速度明顯增大。從圖6d可看出，速度在中心區域基本呈線性分布，當r＞0.8ra時，速度急劇增大，在靠近邊壁的區域內表現出很大的速度梯度。速度沿徑向的分布與剛體旋轉的速度分布（圖1）存在很大差異。

上述結果確實是方程的解，總計算時間相當于圓柱筒旋轉約3 200圈，解很穩定。但試驗數據證明這是一非物理解。研究尋找其物理解的方法將是很重要的任務，具有重要學術意義和實用意義。

5.2 小時間步長的計算

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

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

5.3 計算域旋轉和小時間步長的計算

采用3.3節的計算條件，計算域以vθ旋轉，計算坐標系為柱坐標系，按總能守衡進行求解，使用各向異性能考慮哥氏效應的雷諾應力湍流模型。采用小步長1×10－8s計算。

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

圖6 情況1的縱切面壓強等值線（a）、軸向中心橫截面徑向壓強（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）和徑向壓強（b）分布Fig.7 Circle velocity radial distribution（a）and radial pressure distribution（b）at center transverse section in case 2

軸向中心橫截面處徑向壓強分布和溫度等值線分別示于圖9b、10。由圖9b、10可見，中心區域為300K，側壁附近達324K，不是等溫流動，而且，壓強沿徑向的變化與等溫剛體壓強分布（圖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）和徑向壓強（b）分布Fig.9 Circle velocity radial distribution（a）and radial pressure distribution（b）at center transverse section in case 3

5.4 高密度網格下小時間步長的計算

對圖4的網格在圓周方向和徑向進行加密，總的網格數約為241 600，比原來擴大約22倍。用1×10－7s的時間步長計算，穩態結果如圖11～13所示。

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

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

因此，采用高密度網格的小步長計算能得到物理解，與理論計算結果相比，數值計算結果更符合實際情況。

圖11 情況4下軸向中心橫截面處圓周速度沿徑向（a）和徑向壓強（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 結論

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

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

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

[3]朱自強.應用計算流體力學[M].北京：北京航空航天大學出版社，1998.

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