R 型彎管流動的數(shù)值模擬分析

馬 影, 孫樂萌

(上海工程技術大學 航空運輸學院, 上海201620)

0 引 言

彎管流道廣泛存在于工業(yè)設備中,例如液壓系統(tǒng)、閥門等,流體經過這些彎管時會發(fā)生流動分離,引起漩渦,從而引起壓降和能量損失,在回流區(qū)也可能會產生雜質沉降現(xiàn)象；此外,在拐角處往往會發(fā)生沖蝕現(xiàn)象,影響管道使用壽命。 因此對于彎管流動的準確預測有助于對管道的優(yōu)化設計和質量把控。

一般可采取實驗或者數(shù)值模擬的方法對彎管流場進行分析,主要涉及到流場的速度分布、壓力分布、回流區(qū)大小和分離點等特征。 隨著計算機的高速發(fā)展,計算流體動力學(CFD)方法的高效、準確和低成本較實驗手段有巨大優(yōu)勢［1］。 本文將以CFD 方法對彎管流道進行流動分析,以期獲取可信的流場信息及流動規(guī)律。

1 方法

1.1 結構建模

如圖1 為一典型的R 型彎曲流道,左側為其三維視圖,右側橫截面均為全等的R 型彎曲流道。 流體經由入口流入,經過彎道后從出口流出,由于垂直于流動方向的管道長度相較于流動路徑很長,因此該三維問題可以簡化成二維流動問題。 只需對圖1中右側所示平面幾何進行流動分析,其結果可以代表三維流場中的主要信息。

由于流動雷諾數(shù)的變化會導致出口流動狀態(tài)的改變,為了保證在出口處使流動得到充分發(fā)展,隨著雷諾數(shù)的增大,適當加大拐角后的流道長度,具體幾何流場尺寸見圖2。 可見雷諾數(shù)越大,賦予V2 的值越大,其余尺寸不作改變。

圖1 二維剖面圖Fig. 1 Two-dimensional cross-sectional view

圖2 二維流場Fig. 2 Two-dimensional flow field

1.2 網(wǎng)格劃分

在流體力學問題中,流場信息可以通過求解關于流體微團的控制方程得到,對于不可壓縮粘性流體建立了N-S 方程。 但對于復雜的外形,直接求解該方程是困難的,因此需要對該方程進行離散處理,尋求數(shù)值解。 一般通過對空間的離散處理,即以網(wǎng)格的形式對穩(wěn)態(tài)流場進行求解。 目前的離散方式主要有：有限體積法、有限差分法和有限元法［2］。 目前,ANSYS CFX、ANSYS FLUENT 和OpenFOAM 等都是采用的有限體積法。 本文采取ICEM 對計算域進行網(wǎng)格劃分,由于是二維問題且流場區(qū)域較規(guī)范,故全流域均以四邊形網(wǎng)格進行填充。 圖3 展示了Re=100 對應幾何的以0.2m 作為最大單元尺寸劃分的局部流場網(wǎng)格,網(wǎng)格單元總數(shù)為1090,以長細比和扭曲度為評判標準的網(wǎng)格質量分別在0.483 和0.527 以上,網(wǎng)格質量較高。

圖3 最大單元尺寸0.2m 處的域網(wǎng)格劃分Fig. 3 Domain meshing at the maximum element size of 0.2m

1.3 模型設置

本文解決的是不可壓縮粘性流體的數(shù)值求解,該類問題的控制方程(1)～(3)如下：

其中：p 為流場的無量綱壓力,Re 為雷諾數(shù),ν為運動粘度。 式(1)連續(xù)性方程,反映了流體質量守恒,式(2)和(3)為動量方程,反映了流體動量守恒。 雷諾數(shù)Re 定義為公式(4)：

其中：L 為參考長度,U 為參考速度,μ = ρ·ν,稱為動力粘度。 求解時,設定流體密度ρ 為1 kg／m3,參考長度L 為1m,入口處速度為1m／s,方向沿x 軸,且入口處壓力設為1 個大氣壓；兩壁面邊界均為無滑移壁面；出口為壓力出口,表壓為零。 由式(4)可知,當雷諾數(shù)Re 為100 時,動力粘度μ 應取為0.01 Pa·s。關于Re = 100 時,計算域的邊界命名及邊界條件設置情況,如圖4 所示。 對于管道內流而言,該雷諾數(shù)處于層流區(qū)間,因此流動模型為Laminar,不考慮流體的傳熱效應,流體特性及流動模型的設置見表1。 采取SIMPLE 算法對壓力和速度進行耦合求解,壓力項為二階格式,動量項為二階迎風格式,求解方法設置如見表2。

圖4 流體域邊界條件設置Fig. 4 Setting of boundary conditions of fluid domain

表1 流體域的材料特性Tab. 1 Material properties of fluid domain

表2 求解方法Tab. 2 Solution method

2 結果和討論

2.1 雷諾數(shù)Re ＝100

由上述所給邊界條件,設置連續(xù)項和速度殘差為10 ～6,計算得到了Re = 100 時管道內的靜壓云圖、速度云圖、速度流線圖、Line1 和Arc1 上靜壓關于縱坐標的梯度變化。

(1)靜壓云圖。 從圖5 可知,最大靜壓出現(xiàn)在入口處。 在水平階段,隨著流動發(fā)展,靜壓逐漸降低；當流動到達Wall1 拐角處時,靜壓降至最低；此后流動繼續(xù)發(fā)展,靜壓又逐步升高,最后達到穩(wěn)定狀態(tài)。 另外,當靜壓梯度由正變?yōu)樨摃r,說明發(fā)生了分離流動,即圖5 中標點所示位置。 Arc1 上靜壓關于縱坐標的梯度曲線如圖6 所示,說明在Arc1 上θ =1.090(62.5°)處發(fā)生分離。

圖5 Re＝100 的靜壓云圖Fig. 5 Static pressure cloud diagram with Re＝100

圖6 靜壓梯度曲線Fig. 6 Static pressure gradient curve

(2)速度云圖與速度流線圖。 如圖7 所示,隨著流動發(fā)展,流體速度先增加再減??；處于管道中部的流體速度高于靠近管道壁面的速度,速度最大值出現(xiàn)在管道1／2 徑向處,這是由于流的粘性效應,使得貼近壁面處的流體速度為零。 因此,隨著與壁面的距離增大,流體速度相應的增大,在距離壁面最遠處即1／2 徑向處,流體速度達到最大；從圖8 可知,在出口附近處的流線已沒有回流現(xiàn)象發(fā)生,結合速度和靜壓云圖,可以確定在出口處的流動已經充分發(fā)展。

圖7 Re＝100 的速度云圖Fig. 7 Speed cloud diagram with Re＝100

圖8 Re＝100 的速度流線圖Fig. 8 Speed streamline diagram of Re＝100

圖9 Line1 上靜壓關于y 軸的梯度變化Fig. 9 The gradient of static pressure on Line1 with respect to the y-axis

(3)靜壓關于縱坐標的梯度變化。 回流點處于靜壓梯度由負變?yōu)檎?即如圖9 表點所示。 回流區(qū)域包含如圖8 的螺旋流線,回流長度Lreatt為分離點處弧長與回流點坐標值的絕對值之和,即公式(5)：

式中,θsepar為分離點對應的弧度,yreatt為回流點對應的縱坐標。 因此,Re = 100 對應的回流長度為1.090＋6.693=6.783 m。

2.2 網(wǎng)格獨立性試驗（僅適用于Re = 100 模型）

在正式計算之前,需要排除網(wǎng)格因素對結果的影響,即要進行網(wǎng)格獨立性分析。 針對Re = 100 時的管道流動,建立了5 種類型的網(wǎng)格,分別是coarse, medium, fine, extra＿fine, super＿fine。 他們的區(qū)別是最大網(wǎng)格尺寸不同,即代表了5 種不同的網(wǎng)格密度；就aspect ratio 和skewness 而言,這5 種網(wǎng)格的質量處于同一水平；網(wǎng)格具體信息見表3。 圖10 中對比了由不同密度的網(wǎng)格計算得到的回流長度,通過左縱坐標可以發(fā)現(xiàn)extra＿fine 與super＿fine網(wǎng)格計算結果幾乎一致,右縱坐標顯示了不同網(wǎng)格的計算迭代步數(shù)變化,super＿fine 網(wǎng)格的迭代步數(shù)幾乎是extra＿fine 的3 倍,為了兼顧計算精確度和計算效率,將采用extra＿fine 類網(wǎng)格進行后續(xù)的計算。 比較了5 種不同網(wǎng)格類型下長寬比(Aspect ratio)和偏度(Skewness)的值,如圖11 所示,隨著網(wǎng)格密度的增加,分離點和回流點的位置均逐漸趨于穩(wěn)定。

表3 Re＝100 情況下五個網(wǎng)格的信息比較Tab. 3 Information comparison of five grids under Re＝100

圖10 Re＝100 情況下Lreatt的迭代次數(shù)和網(wǎng)格獨立性測試Fig. 10 The number of iterations and grid independence test under Re＝100

2.3 回流區(qū)尺寸相對于雷諾數(shù)的變化

討論不同雷諾數(shù)對回流區(qū)域的影響。 首先要保證當雷諾數(shù)增大時,管道幾何要保證在出口處的流動是充分發(fā)展的。 通過觀察不同雷諾數(shù)下管道內的靜壓云圖、流線圖和速度云圖等變化,判斷經過拐角后的流動到出口附近時是否已經得到充分發(fā)展。 如圖12 所示,不同雷諾數(shù)下,壓力分布不同,在同一標準圖例下可以得到以下幾條結論：

(1)隨著雷諾數(shù)的增大入口處的壓力逐漸降低；(2)雷諾數(shù)越大,拐角分離點處的角度越大；(3)分離點隨雷諾數(shù)增加而越來越接近；(4)出口處附近的壓力分布均已平穩(wěn)。

圖11 Re＝100 情況下θsepar和yreatt的網(wǎng)格獨立性試驗Fig. 11 Grid independence test of sum under Re＝100

圖12 Re＝50， 100， 150， 200 情況下分別的壓力云圖Fig. 12 Pressure cloud diagrams when Re＝50， 100， 150， 200

此外,還對比了不同雷諾數(shù)下的速度云圖,如圖13 所示,可以得到以下結論：(1)隨著雷諾數(shù)的增大,速度場中能達到最高速度的區(qū)域越大；(2)出口附近區(qū)域的流動速度分布已經穩(wěn)定。

圖13 Re＝50， 100， 150， 200 情況下分別的速度云圖Fig. 13 Velocity cloud diagrams when Re＝50， 100， 150， 200

最后,對比了不同雷諾數(shù)下的流線圖,如圖14所示,發(fā)生回流的區(qū)域隨雷諾數(shù)增大而增大,可以清晰地看到,在出口附近的流線已均勻分布,達到穩(wěn)定狀態(tài)。 結合壓力云圖和速度云圖,可知本文選區(qū)的隨雷諾數(shù)變化而變化的幾何可以滿足在出口處達到穩(wěn)定流動的要求。

圖14 Re＝50，100，150，200 的流線圖Fig. 14 Streamline diagram of Re＝50，100，150，200

為了比較不同雷諾數(shù)下的回流長度,將四組分離點和四組回流點隨雷諾數(shù)的變化曲線分別放到了同一張圖里。 圖15 顯示的是四組雷諾數(shù)下的分離點對比圖。 可見,當Re = 50 時,分離點所處的角度最小,即最靠近Line1；Re = 100, 150, 200 時,對應的分離點很接近,可以認為幾乎不發(fā)生變化,大約都在64°附近；Re = 150, 200 時,流動分離前,壓力梯度會出現(xiàn)震蕩,說明此處的流動情況較為復雜,如果要得到分離前的流場詳細信息,可以在該處進行網(wǎng)格加密。 圖16 是回流點的對比圖,可知隨著雷諾數(shù)的增大,回流點越靠近出口；且當Re = 100, 150,200 時,在回流發(fā)生之前,壓力梯度均有類似的震蕩現(xiàn)象出現(xiàn),說明在該區(qū)域的流動較為復雜；回流發(fā)生后,壓力梯度均不再有明顯變化。 結合式(5)可得到回流長度隨雷諾數(shù)的變化曲線,如圖17 所示,回流長度隨雷諾數(shù)增加而增加。

圖15 Re＝50，100，150，200 的分離點對比Fig. 15 Comparison of separation points with Re＝50，100，150，200

圖16 Re＝50、100、150、200 的回流點對比Fig. 16 Comparison of reflow points with Re＝50， 100， 150， 200

圖17 回流長度隨雷諾數(shù)的變化Fig. 17 Variation of backflow length with Reynolds number

2.4 利用CFD 模擬結果導出預測數(shù)學模型

四組不同雷諾數(shù)下的回流長度,具體取值如表4 所示。

表4 四組回流長度值Tab. 4 Four groups of return length values

假設這四組值具有三次函數(shù)關系,即以雷諾數(shù)為自變量,以回流長度為因變量,可設該函數(shù)關系為式(6)：

其中,ai(i=0,1,2,3) 為未知系數(shù),x 代表雷諾數(shù),f(x) 代表回流長度。 經過簡單的三次函數(shù)擬合后得到的各系數(shù)為：

可以預測到,當Re = 160 時的回流長度為f(160)= 11.786m。 表5 比較了Re=80, 120, 160時回流長度計算值和擬合值的誤差,誤差均在1.23%左右。 可見該擬合公式能較好地預測雷諾數(shù)與回流長度的關系。

表5 回流長度計算值與擬合值比較Tab. 5 Comparison of calculated value and fitted value of reflux length

3 結束語

通過對三維長流道的簡化得到了該模型的二維特征面,使用前處理工具ICEM 對其進行網(wǎng)格劃分,應用商業(yè)軟件FLUENT 對本案例進行計算,使用后處理軟件TECPLOT 繪制壓力云圖、速度云圖和流線圖。 為了盡可能降低網(wǎng)格對于計算結果的影響,對Re = 100 時的網(wǎng)格進行了獨立性檢驗,結果表明當最大網(wǎng)格尺寸取為0.05 時,能兼顧計算精度和計算效率。 為了保證在出口處的流動已經充分發(fā)展,對于流動雷諾數(shù)較大的情況,適當增加出口流道長度。 通過提取求解結果中的壓力梯度值,分析得到流動的分離點和回流點,進而獲取回流長度。 通過雷諾數(shù)Re = 50,100,150,200 時的回流長度值,擬合了回流長度關于雷諾數(shù)的三次方程公式,并預測了Re=160 時的回流長度。 通過研究,得出點結論：

(1)最大靜壓出現(xiàn)在入口處,在水平階段,隨著流動發(fā)展,靜壓逐漸降低,當流動到達Arc1 拐角處時,靜壓降至最低,此后流動繼續(xù)發(fā)展,靜壓又逐步升高,最后達到穩(wěn)定狀態(tài)；

(2)隨著流動發(fā)展,流體速度先增加再減小,處于管道中部的流體速度高于靠近管道壁面的速度,速度最大值出現(xiàn)在管道1／2 徑向處；

(3)隨著雷諾數(shù)的增大,入口處的壓力逐漸降低；

(4)雷諾數(shù)越大,拐角分離點處的角度越大；

(5)分離點隨雷諾數(shù)增加而越來越接近；

(6) Re=150,200 時,流動分離前,壓力梯度會出現(xiàn)震蕩；

(7) Re = 100, 150, 200 時,在回流發(fā)生之前,壓力梯度均有類似的震蕩現(xiàn)象出現(xiàn)；

(8)擬合得到的公式能較好的預測其他雷諾數(shù)對應的回流長度。