R 型彎管流動的數值模擬分析

馬 影, 孫樂萌

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

0 引 言

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

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

1 方法

1.1 結構建模

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

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

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

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

1.2 網格劃分

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

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

1.3 模型設置

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

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

其中：L 為參考長度,U 為參考速度,μ = ρ·ν,稱為動力粘度。 求解時,設定流體密度ρ 為1 kg／m3,參考長度L 為1m,入口處速度為1m／s,方向沿x 軸,且入口處壓力設為1 個大氣壓；兩壁面邊界均為無滑移壁面；出口為壓力出口,表壓為零。 由式(4)可知,當雷諾數Re 為100 時,動力粘度μ 應取為0.01 Pa·s。關于Re = 100 時,計算域的邊界命名及邊界條件設置情況,如圖4 所示。 對于管道內流而言,該雷諾數處于層流區間,因此流動模型為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 雷諾數Re ＝100

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

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

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

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

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

圖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)靜壓關于縱坐標的梯度變化。 回流點處于靜壓梯度由負變為正,即如圖9 表點所示。 回流區域包含如圖8 的螺旋流線,回流長度Lreatt為分離點處弧長與回流點坐標值的絕對值之和,即公式(5)：

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

2.2 網格獨立性試驗（僅適用于Re = 100 模型）

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

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

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

2.3 回流區尺寸相對于雷諾數的變化

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

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

圖11 Re＝100 情況下θsepar和yreatt的網格獨立性試驗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

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

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

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

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

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

圖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 回流長度隨雷諾數的變化Fig. 17 Variation of backflow length with Reynolds number

2.4 利用CFD 模擬結果導出預測數學模型

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

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

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

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

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

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

3 結束語

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

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

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

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

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

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

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

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

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