基于OpenFOAM的除塵風(fēng)機(jī)的高性能數(shù)值模擬

秦宣云 文金侶 朱世存

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

1 引言

除塵用的風(fēng)機(jī)屬于葉輪機(jī)械的一種,在研究除塵的風(fēng)機(jī)內(nèi)流場(chǎng)過(guò)程中,一方面,實(shí)驗(yàn)工作有助于理解流動(dòng)復(fù)雜性,另一方面,由于近幾十年來(lái),數(shù)值模擬的興起,數(shù)值模擬方法已被證明是對(duì)實(shí)驗(yàn)的有效補(bǔ)充,計(jì)算流體力學(xué)(CFD)在流場(chǎng)分析中起到越來(lái)越重要的作用.

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

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

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

2 基本理論和方法

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

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

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

2.2 粉塵顆粒運(yùn)動(dòng)軌跡方程

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

并且粉塵粒子的運(yùn)動(dòng)軌跡受牛頓第二定律控制:

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

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

τp是粒子響應(yīng)局部流速變化所需的時(shí)間,稱為粒子的松弛時(shí)間，按下式計(jì)算:

其中,CD是曳力系數(shù)[12],其標(biāo)準(zhǔn)定義為:

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

2.3 除塵風(fēng)機(jī)k?ω?SST湍流控制方程

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

其中

3 氣相流場(chǎng)邊界條件和初始條件

3.1 邊界條件設(shè)置

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

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

3.2 粒子的初始條件

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

表1 初始條件設(shè)置

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

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

4.1 除塵風(fēng)機(jī)幾何模型

除塵風(fēng)機(jī)的幾何模型圖如圖1所示.

圖1 除塵風(fēng)機(jī)幾何模型

4.2 除塵風(fēng)機(jī)的cycileAMI動(dòng)網(wǎng)格

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

圖3 風(fēng)機(jī)的網(wǎng)格劃分

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

圖4 風(fēng)機(jī)的AMI網(wǎng)格

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

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

5.1.1 控制方程時(shí)間離散

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

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

對(duì)(5.1)式積分,得到:

對(duì)空間的離散化,可以表示為:

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

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

化簡(jiǎn)得:

即

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

由(5.3)和(5.4)式,得

5.1.2 控制方程空間離散

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

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

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

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

5.1.3離散后的線性方程組

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

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

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

5.2 并行算法

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

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

6 計(jì)算結(jié)果及其分析

6.1 除塵風(fēng)機(jī)的速度和壓力分布結(jié)果

采用Paraview軟件對(duì)除塵風(fēng)機(jī)的采用瞬時(shí)態(tài)計(jì)算穩(wěn)定下來(lái)的速度和壓力結(jié)果進(jìn)行可視化處理,所得結(jié)果如圖6所示.

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

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

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

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

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

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

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

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

圖9 壓力等值線分布

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

圖10 為在風(fēng)機(jī)內(nèi)氣相流程的瞬態(tài)計(jì)算穩(wěn)定后,對(duì)風(fēng)機(jī)注入5μm的粉塵粒子,且粉塵粒子受Lagrange軌跡運(yùn)動(dòng)方程控制的氣相流場(chǎng)的壓力分布和粉塵粒子的速度分布圖.從圖中可以看出,粒子沿著葉片徑向運(yùn)動(dòng),受到葉輪旋轉(zhuǎn)的影響,以一定幅度甩向壁面.進(jìn)入風(fēng)機(jī)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).從粒子的運(yùn)動(dòng)軌跡可以看出,粒子沿著葉輪葉片底部徑向以螺旋線似的方式運(yùn)動(dòng)到葉尖,然后部分粉塵粒子甩向壁面.

圖10 5um粉塵粒子在風(fēng)機(jī)內(nèi)的運(yùn)動(dòng)情況

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

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

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

7 結(jié)論

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

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

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