盛 譽(yù) 路 昕 凌乃陽(yáng) 周陳龍（核工業(yè)理化工程研究院，天津300180）

0 引言

近年來，隨著計(jì)算機(jī)科學(xué)技術(shù)的迅速發(fā)展，計(jì)算流體動(dòng)力學(xué)（CFD）[1]得到了越來越廣泛的應(yīng)用，已成為許多工程領(lǐng)域流動(dòng)問題研究、設(shè)計(jì)的主要手段之一。但是，隨著發(fā)展的逐漸深入，CFD也面臨著越來越多的困難，特別是CFD中的網(wǎng)格效應(yīng)問題，迄今為止很難給出一個(gè)明確的答案，常常只能通過經(jīng)驗(yàn)解決。

CFD是通過數(shù)值方法求解在空間域上離散得到的控制方程，而在空間域上離散控制方程必須使用網(wǎng)格[4]。可見，網(wǎng)格分布將直接影響數(shù)值計(jì)算的精度和穩(wěn)定性。對(duì)于旋轉(zhuǎn)流體中的一般工程問題來說，雷諾數(shù)不大，網(wǎng)格效應(yīng)對(duì)計(jì)算結(jié)果的影響不會(huì)太明顯，在某些條件下甚至可以得到精確解。但是，對(duì)于高速旋轉(zhuǎn)流場(chǎng)，徑向壓強(qiáng)和密度成指數(shù)分布。為了確定流場(chǎng)，數(shù)值求解所需網(wǎng)格的分布應(yīng)與流場(chǎng)的物理規(guī)律相協(xié)調(diào)。在為數(shù)不多的關(guān)于數(shù)值計(jì)算網(wǎng)格效應(yīng)的研究論文中，至今尚未發(fā)現(xiàn)有針對(duì)這種高雷諾數(shù)下旋轉(zhuǎn)流場(chǎng)網(wǎng)格效應(yīng)問題進(jìn)行的研究，而如何確定合適的網(wǎng)格、近壁網(wǎng)格尺寸多少，以便使計(jì)算的精度高、可靠性好，也未有明確的答案。

為解決以上問題，文章采用有限體積法，通過求解柱坐標(biāo)系下的流體動(dòng)力學(xué)方程[2]，對(duì)高速旋轉(zhuǎn)狀態(tài)下圓筒空腔內(nèi)部三維流場(chǎng)的網(wǎng)格尺度效應(yīng)問題進(jìn)行數(shù)值研究。

1 計(jì)算域模型與網(wǎng)格劃分

本文研究的是空筒的旋轉(zhuǎn)問題，采用的模型如圖1所示。計(jì)算模型為半徑ra=0.06 m、高h(yuǎn)=0.01 m。針對(duì)其內(nèi)部高速旋轉(zhuǎn)的流場(chǎng)研究側(cè)壁附近網(wǎng)格對(duì)流場(chǎng)計(jì)算準(zhǔn)確性的影響問題。

圖1 流場(chǎng)計(jì)算域模型

網(wǎng)格分布直接影響到數(shù)值計(jì)算的精度和穩(wěn)定性，采用了結(jié)構(gòu)化的六面體網(wǎng)格，如圖2和圖3所示，根據(jù)引言對(duì)流場(chǎng)的分析可知，流動(dòng)參數(shù)在軸向、角向梯度不大，徑向梯度很大，且大部分氣體都集中在側(cè)壁附近。因此，網(wǎng)格劃分的原則是在近側(cè)壁10%半徑的區(qū)域內(nèi)單獨(dú)分區(qū)，且徑向從中心軸線到側(cè)壁逐漸加密。

圖2 流場(chǎng)計(jì)算域網(wǎng)格

圖3 流場(chǎng)計(jì)算域?qū)ΨQ面網(wǎng)格

對(duì)近壁區(qū)的流動(dòng)，采用壁面函數(shù)[3]的處理方法。壁面函數(shù)法依賴于近壁第一層網(wǎng)格，在劃分網(wǎng)格時(shí)，需要把第一個(gè)網(wǎng)格節(jié)點(diǎn)布置在對(duì)數(shù)分布規(guī)律成立的范圍內(nèi)，通常通過調(diào)節(jié)無量綱長(zhǎng)度參數(shù)y+保證計(jì)算結(jié)果的準(zhǔn)確度[5]。近壁面網(wǎng)格無量綱長(zhǎng)度參數(shù)y+的定義式為：

其中：Δy為近壁第一個(gè)網(wǎng)格節(jié)點(diǎn)到壁面的距離；ρ為流體密度；uτ為壁面摩擦速度；τw為壁面切應(yīng)力；up、uw分別為近壁第一個(gè)網(wǎng)格節(jié)點(diǎn)處的速度、壁面速度。

表1給出了不同的近壁面Δy所對(duì)應(yīng)的y+值。為了保證不同方案間計(jì)算結(jié)果的可比性，在改變近壁網(wǎng)格距離時(shí)僅對(duì)近壁面單獨(dú)分區(qū)的網(wǎng)格密度進(jìn)行調(diào)整[6]。為了保證數(shù)值計(jì)算的收斂性，一般應(yīng)保證網(wǎng)格的縱橫比（最大邊/最小邊長(zhǎng)度的比率）小于1000，在求解高雷諾數(shù)的邊界層問題時(shí)允許達(dá)到105～106。這里流場(chǎng)的雷諾數(shù)達(dá)到1.2×106，考慮到計(jì)算周期，采用的計(jì)算網(wǎng)格最大縱橫比為2.8×105，網(wǎng)格節(jié)點(diǎn)總數(shù)均為200000。

表1 不同近壁網(wǎng)格距離的y+值

2 計(jì)算方法

2.1 流體動(dòng)力學(xué)基本方程

在高速旋轉(zhuǎn)的空筒內(nèi)部，氣體流動(dòng)所遵守的流體動(dòng)力學(xué)基本方程可以寫成以下的通用形式：

式中，φ為通用變量，代表u、v、w、T等求解變量；Γ為廣義擴(kuò)散系數(shù)；S為廣義源項(xiàng)。對(duì)于不同的方程，Γ、S具有相應(yīng)的形式，表2給出了這些符號(hào)與方程的對(duì)應(yīng)關(guān)系，其中，cp為比熱容；k為導(dǎo)熱系數(shù)；S為黏性耗散項(xiàng)。

表2 通用控制方程中符號(hào)對(duì)應(yīng)的形式

2.2 初始條件與邊界條件

在建立了流體動(dòng)力學(xué)運(yùn)動(dòng)基本方程以后，需要確定運(yùn)動(dòng)的初始狀態(tài)和邊界狀態(tài)才具有唯一的形態(tài)。初始條件為氣體旋轉(zhuǎn)運(yùn)動(dòng)的等溫剛體條件，即對(duì)于空筒內(nèi)部的氣體，氣體溫度T=T0=300 K，旋轉(zhuǎn)速度ω=ω0=3 000 rad/s。空筒的上下端蓋和側(cè)壁定義為壁面邊界，壁面旋轉(zhuǎn)速度和溫度與空筒內(nèi)部氣體的初始狀態(tài)保持一致。

采用Fluent16.0中基于密度的求解器中的隱式方法對(duì)流場(chǎng)求解。根據(jù)流場(chǎng)特性，雷諾數(shù)達(dá)到106量級(jí)，因此，需要采用湍流模型進(jìn)行求解，根據(jù)經(jīng)驗(yàn)選取Transition SST模型。

3 網(wǎng)格尺度對(duì)計(jì)算結(jié)果的影響

針對(duì)近壁面第一層網(wǎng)格在壁面法線方向長(zhǎng)度分別為1.0×10-4m（方案1）、1.0×10-5m（方案2）、1.0×10-6m（方案3）、1.0×10-7m（方案4）時(shí)的旋轉(zhuǎn)流場(chǎng)進(jìn)行了數(shù)值計(jì)算。由于最小網(wǎng)格尺度不同，計(jì)算所需的時(shí)間步長(zhǎng)也有所變化，網(wǎng)格尺度越小[7]，計(jì)算時(shí)間步長(zhǎng)也越小。對(duì)于尺度很小的網(wǎng)格，若強(qiáng)制采用大的時(shí)間步長(zhǎng)，計(jì)算容易發(fā)散，甚至得到非物理解。方案1和方案2的計(jì)算在時(shí)間步長(zhǎng)Δt=1.0×10-7s的情況下很快收斂，方案3和方案4的計(jì)算相當(dāng)緩慢，方案4的Δt甚至小于了1.0×10-8s，經(jīng)過長(zhǎng)時(shí)間的計(jì)算得到了收斂結(jié)果。

圖4顯示了四種方案所得到的氣體旋轉(zhuǎn)速度沿徑向的分布情況。在穩(wěn)定狀態(tài)下，圓筒內(nèi)的氣體運(yùn)動(dòng)理論上遵循剛體運(yùn)動(dòng)的原則，即氣體的旋轉(zhuǎn)速度沿徑向呈線性分布。從圖中可以看出，方案1、2、3所計(jì)算出的氣體旋轉(zhuǎn)速度并不完全遵守氣體沿徑向呈線性分布的規(guī)律，大部分氣體的旋轉(zhuǎn)速度均滯后于相應(yīng)的理論值。但當(dāng)近壁面第一層網(wǎng)格法向長(zhǎng)度減小至1.0×10-7m（即方案4）時(shí)，氣體的旋轉(zhuǎn)速度曲線已基本與理論值吻合。圖5顯示了通過方案4所計(jì)算出的氣體旋轉(zhuǎn)速度云圖。

圖4 旋轉(zhuǎn)速度沿徑向分布

圖5 旋轉(zhuǎn)速度分布云圖

在高速旋轉(zhuǎn)的圓筒內(nèi)部，氣體壓強(qiáng)沿徑向呈指數(shù)分布，越靠近外壁，氣體壓強(qiáng)分布越陡。圖6顯示了四種方案所得到的氣體壓強(qiáng)沿徑向的分布情況。方案1,2,3所得到的氣體壓強(qiáng)分布差別不大，大部分區(qū)域的氣體基本遵守徑向指數(shù)分布的原則，而在徑向上0.6倍半徑的區(qū)域附近，氣體壓強(qiáng)明顯低于理論值。對(duì)比之下，方案4所得到的結(jié)果雖與理論值稍有偏差，但其分布曲線已與理論值曲線十分接近。圖7顯示了通過方案4計(jì)算出的氣體壓強(qiáng)云圖。

圖6 壓強(qiáng)沿徑向分布

圖7 壓強(qiáng)分布云圖

表3給出了不同方案間的轉(zhuǎn)子壁面功耗比較，近壁面網(wǎng)格尺度從10-4遞減至10-7的過程中，計(jì)算出的功耗逐步減小。根據(jù)資料，相同轉(zhuǎn)速與氣體介質(zhì)條件下，半徑0.06 m，高0.1 m長(zhǎng)轉(zhuǎn)子負(fù)載運(yùn)行時(shí)，轉(zhuǎn)子壁面消耗的能量約為0.52 W。本文的研究對(duì)象可視為長(zhǎng)轉(zhuǎn)子軸向的一部分。由于這里圓筒壁面的軸向高度僅為長(zhǎng)轉(zhuǎn)子壁面軸向高度的約10%，可知方案1、2、3所計(jì)算出的壁面功耗偏大，方案4所計(jì)算出的壁面功耗接近實(shí)驗(yàn)值。

表3 不同方案間的功耗比較

4 結(jié)論

本文采用有限體積法，求解了柱坐標(biāo)系下的N-S方程，對(duì)雷諾數(shù)達(dá)到106量級(jí)的高速旋轉(zhuǎn)流場(chǎng)進(jìn)行了數(shù)值研究。針對(duì)不同網(wǎng)格尺度的流場(chǎng)計(jì)算結(jié)果，采用數(shù)值分析、理論、實(shí)驗(yàn)相結(jié)合的方法進(jìn)行了對(duì)比。得到如下結(jié)論：

（1）近壁面網(wǎng)格的法向尺度達(dá)到10-7時(shí)，計(jì)算結(jié)果最接近實(shí)際流動(dòng)狀態(tài)。氣體的旋轉(zhuǎn)速度和氣體壓力在徑向上分別接近線性分布和指數(shù)分布，轉(zhuǎn)子壁面功耗最接近實(shí)驗(yàn)值。

（2）在高雷諾數(shù)的流場(chǎng)計(jì)算時(shí)，應(yīng)根據(jù)實(shí)際流動(dòng)條件來選擇合適的近壁面網(wǎng)格尺度。