秦宣云 文金侶 朱世存

(1.中南大學數(shù)學與統(tǒng)計學院,長沙,410083?2.長沙智能駕駛研究院,長沙,410208)

1 引言

除塵用的風機屬于葉輪機械的一種,在研究除塵的風機內(nèi)流場過程中,一方面,實驗工作有助于理解流動復雜性,另一方面,由于近幾十年來,數(shù)值模擬的興起,數(shù)值模擬方法已被證明是對實驗的有效補充,計算流體力學(CFD)在流場分析中起到越來越重要的作用.

OpenFOAM是用C++開發(fā)的一款嚴謹和高質(zhì)量的開源庫[1],它實際上是面向?qū)ο蟮钠⒎址匠虜?shù)值模擬的工具,并且擁有活躍的開源社區(qū).本文基于OpenFOAM開源框架做風機的數(shù)值模擬,通過開發(fā)適用于除塵風機的數(shù)值模擬求解器PimpleLPTDyFoam,提供了對除塵風機數(shù)值模擬的一種探索.

OpenFOAM開源框架對除塵風機做數(shù)值模擬主要涉及到旋轉(zhuǎn)機械幾何模型、動網(wǎng)格技術、湍流模型、并行計算和拉格朗日粒子追蹤等模塊.在旋轉(zhuǎn)機械的研究方面,Wilhelm D[14]使用OpenFOAM研究旋轉(zhuǎn)流動問題的數(shù)值模擬,介紹了MRF等技術在OpenFOAM的應用?在動網(wǎng)格技術方面,有GGI網(wǎng)格接口技術,用于處理兩套網(wǎng)格連接表面不相匹配的情況[15],P.E.Farrell和J.R.Maddison提出了一個能運用在不連接但相鄰的網(wǎng)格區(qū)域上的AMI[16]網(wǎng)格接口技術,這種任意網(wǎng)格接觸面的算法可應用于非結構化網(wǎng)格的幾何邊界上?在旋轉(zhuǎn)機械下的拉格朗日粒子追蹤方面,Adel Ghenaiet研究了葉輪風機內(nèi)粒子的運動軌跡和沖蝕[17]?在使用OpenFOAM做風機葉片的研究上,Rickard E Bensow[18]使用OpenFOAM研究了非穩(wěn)態(tài)的轉(zhuǎn)子葉片的負載?Tuomas Turunen等[19]用OpenFOAM研究了螺旋葉輪的Open-Water數(shù)值模擬?O.Petit等[20]采用OpenFOAM研究了旋流產(chǎn)生器產(chǎn)生的非穩(wěn)態(tài)流動模擬?另外,Miao Wang等[21]用OpenFOAM對各種圖劃分算法的性能做了分析.

本文基于OpenFOAM開源框架對除塵風機的氣相做數(shù)值模擬主要涉及到旋轉(zhuǎn)機械幾何模型、動網(wǎng)格技術、湍流模型和并行計算等模塊,通過風機內(nèi)氣相流場的Navier-Stokes輸運方程和k?ω?SST湍流方程模型和數(shù)值離散方法,研究在特定邊界條件和初始條件下,氣相流場場量的分布規(guī)律.另外,通過氣象耦合的常微分方程對粉塵運動進行建模,并在特定初始條件下,對粉塵顆粒在除塵風機內(nèi)的運動進行數(shù)值模擬,計算得到粉塵在風機內(nèi)的運行軌跡和速度大小.

2 基本理論和方法

2.1 不可壓縮的Navier-Stokes控制方程

風機內(nèi)的氣相流場由下面不可壓縮的Navier-Stokes方程控制[11]:

其中v,p,ρf,μf分別代表氣相速度,壓力,密度和氣相流體的粘度.

2.2 粉塵顆粒運動軌跡方程

粉塵顆粒p,可以由它的中心位置xp,它的直徑Dp,它的速度vp和它的密度ρp來確定.粉塵顆粒的質(zhì)量在拉格朗日框架下,每個粉塵顆粒的位置xp滿足下面的方程式:

并且粉塵粒子的運動軌跡受牛頓第二定律控制:

對一確定的顆粒,其所受的主要的力為:

其中FS為離心力和柯氏力,FD為曳力,可以表示為:

τp是粒子響應局部流速變化所需的時間,稱為粒子的松弛時間，按下式計算:

其中,CD是曳力系數(shù)[12],其標準定義為:

而Rep是粒子的雷諾數(shù),定義為:

2.3 除塵風機k?ω?SST湍流控制方程

除塵風機的數(shù)值模擬采取k?ω?SST湍流模型,因為k?ω?SST湍流模型組合了一種高雷諾數(shù)k??湍流模型[10]和k?ω湍流模型[2]的一些優(yōu)勢?k?ω?SST模型如下[3]:

其中

3 氣相流場邊界條件和初始條件

3.1 邊界條件設置

入口邊界條件:v=(0,0,?1),對壓力p應用Neumann邊界條件,即對湍動能k應用Dirichlet邊界條件,即:k=0?intensity值為0.05.

壁面邊界采取noSlip壁面,k=0.

3.2 粒子的初始條件

粒子從入口網(wǎng)格片上注入到風機,采用patchInjection注射模型,每隔0.025秒注射一次,注射3326個粒子,粒子的尺寸分布設為均勻分布,均值為5μm,注射的粒子速度為:v=(10,0,0)m/s.

表1 初始條件設置

表2 kOmegaSST湍流模型的常數(shù)項

4 旋轉(zhuǎn)幾何模型和動網(wǎng)格模型

4.1 除塵風機幾何模型

除塵風機的幾何模型圖如圖1所示.

圖1 除塵風機幾何模型

4.2 除塵風機的cycileAMI動網(wǎng)格

圖2 Impeller,suction,volute幾何的網(wǎng)格

圖3 風機的網(wǎng)格劃分

為了實現(xiàn)動網(wǎng)格的特性,引入AMI網(wǎng)格層[8],考慮到除塵風機的幾何對稱性,這里采用采用滑移接口,即cycileAMI網(wǎng)格[9]:Impeller,suction,volute的AMI網(wǎng)格(圖4).

圖4 風機的AMI網(wǎng)格

5 氣相流場數(shù)值離散方法和并行算法

5.1 控制方程的數(shù)值離散

5.1.1 控制方程時間離散

為了方便表示瞬態(tài)下時間的數(shù)值離散,可以把控制方程表示成:

其中,?(φ)表示空間算子,包括對流項,散度項和源項等.

對(5.1)式積分,得到:

對空間的離散化,可以表示為:

其中,VC表示離散單元的體積,?(?t)表示在t時刻空間離散化的離散算子.

對(5.2)式在時間區(qū)間[t??t/2,t+?t/2]進行積分,得到

化簡得:

即

本文采取的時間離散格式為一階精度的歐拉隱式格式,即使用一階迎風插值,在離散體積中心,ρ?的值由如下關系給出[4]:

由(5.3)和(5.4)式,得

5.1.2 控制方程空間離散

本文采取的空間離散方式為[5]

梯度格式有高斯積分和最小二乘法兩種處理方式[6].通過最小二乘法計算梯度的格式為

另外,曲面法向梯度(??)f的離散為

本文中梯度通過高斯積分并且從單元中心的值到面中心的值的插值方法為線性插值來離散:

5.1.3離散后的線性方程組

由(5.1)式的時間離散和(5.2),(5.5)式等空間離散可以得到線性方程:

對于每個控制體積單元C,都可以得到形如(5.6)式的一個方程式.上述方程式表示?tC依賴于與控制體積單元C與之相鄰的單元的值.所有的這些控制體積單元內(nèi)的線性方程形成了一個線性方程組:

其中[A]是一個稀疏矩陣,系數(shù)αC是其對角元素?[?]是所有控制體積中心處的場量值.壓力校準這里采用代數(shù)幾何多重網(wǎng)格求解(GAMG)[7].

5.2 并行算法

本文的除塵風機是3維幾何模型,為了加快計算速度,將整個風機幾何模型劃分為四塊區(qū)域,并選用Scotch算法[13]劃分的網(wǎng)格進行并行計算,不同顏色代表不同處理器處理的網(wǎng)格區(qū)域,采用四個處理器來計算的網(wǎng)格劃分結果如圖5(在圖5中,紅色表示Processor0處理的區(qū)域,綠色表示Processor1處理的區(qū)域,黃色表示Processor2處理的區(qū)域,藍色表示Processor3處理的區(qū)域).

圖5 并行計算的網(wǎng)格區(qū)域劃分

6 計算結果及其分析

6.1 除塵風機的速度和壓力分布結果

采用Paraview軟件對除塵風機的采用瞬時態(tài)計算穩(wěn)定下來的速度和壓力結果進行可視化處理,所得結果如圖6所示.

圖6 瞬態(tài)計算穩(wěn)定下來的速度和壓力分布圖

對風機在初始速度v=(0,0,0)m/s,轉(zhuǎn)速為150rad/s的工況條件下的數(shù)值模擬計算結果,圖7顯示了其在不同時刻繞旋轉(zhuǎn)中心(x,y,z)=(0,0,0)的速度分布?圖8顯示了其在不同時刻繞旋轉(zhuǎn)中心(x,y,z)=(0,0,0)的壓力分布.

從圖7可以看出,在不同時刻葉片轉(zhuǎn)過的位置不同,但到0.016秒時刻以后的速度分布云圖均已穩(wěn)定,說明此時計算已經(jīng)達到穩(wěn)定狀態(tài)?其次,風機內(nèi)大部分非旋轉(zhuǎn)部件顏色為深藍色,對應接近1m/s位置的顏色,說明風機在初始速度v=(0,0,0)m/s,轉(zhuǎn)速為150rad/s的工況條件下非旋轉(zhuǎn)部件的速度比較低,在轉(zhuǎn)動部件和非旋轉(zhuǎn)部件之間的速度分布在10m/s左右,而在旋轉(zhuǎn)區(qū)域內(nèi)的速度分布在16m/s到32m/s之間?這說明除塵風機的旋轉(zhuǎn)機械在運動的過程中大部分流體是直接流過葉片表面的,且在葉片周圍會對來流產(chǎn)生影響.圖7顯示葉尖的速度最大,葉尖對應的最大速度為32.85m/s?從能量守恒的角度分析,來流的能量轉(zhuǎn)換為葉片的動能,所以葉片經(jīng)過的部分其速度減小?在葉尖處相對于周圍流體速度有所升高,是因為葉尖帶動流體向葉尖運動的方向運動,越靠近葉尖的部分速度越大.

圖7 不同時刻繞旋轉(zhuǎn)中心的速度分布圖

圖8 為在不同時刻,風機葉片繞旋轉(zhuǎn)中心(x,y,z)=(0,0,0)的壓力分布云圖.圖中顯示,到0.016秒時刻以后的速度分布云圖均已穩(wěn)定?穩(wěn)定后的計算結果表明風機的入口基本處于淡藍色區(qū)域,即意味著壓力為負壓,且壓力的絕對值比較大,超過150m2/s2?而出口的壓力大部分處于紅色區(qū)域,表示壓力絕對值比較低,接近于0.從圖8中可以看出風機的葉片在旋轉(zhuǎn)的過程中,葉片背面表面呈藍色,說明表面的流體受到的是負壓?同時,葉片正面表面呈深紅色,說明葉片正表面受到的壓力為正,壓力增大.

圖8 不同時刻繞旋轉(zhuǎn)中心的壓力分布圖

6.2 穩(wěn)態(tài)下壓力等值面

圖9 為壓力計算基本到達穩(wěn)態(tài)后,不同壓力值下的壓力等值面分布圖.圖9(a)為壓力每增加10單位值形成的一簇等值面.從圖9(a)的壓力等值面中,可以看出在葉片尾部形成負壓,葉片尖部基本為正壓?并且壓力等值面呈現(xiàn)一定的規(guī)則性.圖9(d)為壓力值為0m2/s2形成的等值面,有四個凸起的部分,對應于四個葉片.等壓面凸起部分沿著葉片旋轉(zhuǎn)方向往前傾一定位置,在凸起部分之間等壓面凹陷,呈現(xiàn)凸凹交替.結合圖9(a)可以發(fā)現(xiàn),壓力值從0m2/s2開始增加時,基本上等值面變化的趨勢并未改變.圖9(c)為壓力值為-180m2/s2對應的等值面,基本上分布在impeller區(qū)域內(nèi).

圖9 壓力等值線分布

6.3 粉塵粒子的數(shù)值模擬計算結果

圖10 為在風機內(nèi)氣相流程的瞬態(tài)計算穩(wěn)定后,對風機注入5μm的粉塵粒子,且粉塵粒子受Lagrange軌跡運動方程控制的氣相流場的壓力分布和粉塵粒子的速度分布圖.從圖中可以看出,粒子沿著葉片徑向運動,受到葉輪旋轉(zhuǎn)的影響,以一定幅度甩向壁面.進入風機impeller inlet的速度大約在0.1m/s左右,隨著葉輪的旋轉(zhuǎn),粉塵在impeller內(nèi)的速度分布在5m/s到20m/s內(nèi),隨著粒子被甩出到volute,在volute內(nèi)的粉塵粒子速度較低,大概分布在0m/s到10m/s內(nèi).從粒子的運動軌跡可以看出,粒子沿著葉輪葉片底部徑向以螺旋線似的方式運動到葉尖,然后部分粉塵粒子甩向壁面.

圖10 5um粉塵粒子在風機內(nèi)的運動情況

為了使圖像更清晰,只顯示四個代表性的粒子的軌跡.圖11顯示的是四個粒子在不同時刻軌跡的變化,圖12顯示的是四個粒子速度大小變化的曲線.圖12中粒子1對應的是圖11中白色粒子,粒子4對應圖11的黃色粒子,粒子1和4的運動軌跡是從impeller內(nèi)到impeller外的suction幾何部分,它們的速度大小在時間幀從6到8時有波動,其他時間幀速度大小單調(diào)下降.粒子2和粒子3的運動軌跡是在impeller內(nèi),它們的速度變化比較復雜.

圖11 粒子的速度大小隨時間的變化情況

圖12 粒子速度大小變化曲線

7 結論

本文基于OpenFOAM開源框架,開發(fā)出適用于除塵風機的數(shù)值模擬求解器Pimple-LPTDyFoam.求解器涉及到旋轉(zhuǎn)機械模塊、動網(wǎng)格模塊、湍流模型和Navier-Stokes輸運方程、拉格朗日粒子追蹤模塊,以及并行計算和數(shù)值離散方法.氣相流場主要受Navier-Stokes偏微分方程控制,而粉塵粒子的運動受二階非線性的常微分方程控制,該常微分方程耦合了氣相流場的速度場量.對這兩個系統(tǒng)通過PimpleLPTDyFoam求解器求解,從而實現(xiàn)對除塵風機內(nèi)的氣相流場和粉塵粒子的數(shù)值模擬,得到關于場量隨時間變化的時間序列,并基于此進行動力學分析.本文主要結論如下:

1.通過動網(wǎng)格方法對風機內(nèi)部的氣相流場的瞬態(tài)變化進行數(shù)值模擬計算,得到了氣相流場的速度、壓力、湍動能和湍耗散率等數(shù)據(jù),并用paraview將計算得到的數(shù)據(jù)可視化,得到了速度、壓力等場量在風機內(nèi)部不同區(qū)域的大小分布情況.結果表明,在風機葉尖處的速度比較大,并且隨著葉片網(wǎng)格的運動,有“尾帶”效應.

2.在氣相流場瞬態(tài)計算相對穩(wěn)定后,通過以一定速率注入粉塵粒子,模擬粉塵粒子在風機內(nèi)的軌跡和速度變化,得到了粒子在風機內(nèi)隨網(wǎng)格運動時,粒子的軌跡和速度的瞬時變化規(guī)律.結果表明,粒子沿著葉片徑向運動,受到葉輪旋轉(zhuǎn)的影響,以一定幅度甩向壁面.