離心泵內纖維懸浮流的PIV實驗及CFD-DEM模擬研究

張啟華 張為棟 閆召旭 謝之昂 康 順

（江蘇大學流體機械工程技術研究中心，江蘇鎮江，212013）

采用離心泵輸送纖維懸浮液在制漿、造紙流程及廢水處理過程中應用廣泛。尤其近年來，隨著紙張消費量的上升，作為紙漿循環的動力設備，紙漿泵所產生的能耗約占工藝流程總能耗的15%。同時，泵的穩定運行對造紙終端產品的質量也至關重要。因此，深入開展紙漿泵內部流動的研究，對了解長纖維懸浮液的流動特性及推進紙漿泵節能降耗均有重要的工程意義。

纖維懸浮液的運動規律以及纖維的取向分布狀態造成的流變特征，是纖維懸浮液研究的兩個主要方面。可通過引入黏度作為剪切速率的函數來表征紙漿纖維懸浮液的非牛頓行為。Cotas等［1］提出了一種擬均勻方法來研究管道中桉樹纖維懸浮液的湍流。目前，大多數泵內兩相流模擬采用基于Mixture的模型，即視兩相均為具有連續速度場的流動介質。Liu等［2］基于Mixture模型模擬了粒徑1～12 mm的顆粒在不同工況下的運動狀態并總結了不同粒徑粒子的分布規律。Zhang等［3］基于Mixture模型研究了顆粒濃度、體積分數及粒徑對泵性能的影響。但Mixture模型忽略了顆粒間的碰撞以及顆粒與流體的相對運動，存在一定的模擬誤差。且顆粒體積的變大或顆粒濃度的提高，顆粒的拓撲結構、自身運動以及顆粒間相互作用不可忽略。隨著待解決問題的復雜程度越來越高，具有較高計算效率的離散元（DEM）法在泵兩相流工程中得到了推廣應用。基于豐富的顆粒表征、更復雜的顆粒初始化，以及顆粒-顆粒與顆粒-部件的接觸模型，DEM法能更直觀地展現各類復雜流程中顆粒的運動。Varga等［4］使用流體力學（CFD）-DEM方法模擬管內固液兩相流動并與實驗中不同材料顆粒對管道的磨損情況結果對比，發現材料屬性對管壁磨損程度影響較大。Huang等［5］采用CFD-DEM方法具體分析泵內顆粒的運動情況，發現顆粒趨向于在葉片工作面聚集。以往泵工程的CFD-DEM研究主要以球顆粒為主，而對于長顆粒方面的研究仍鮮見報道。

目前，兩相流可視化中顆粒相的追蹤主要采用高速攝影［6］。近年來，粒子圖像測速技術（PIV）以其非接觸式測量方式，通過瞬態記錄大量空間點上的速度分布信息，從而能夠提供豐富的流場空間結構以及流動細節。通過PIV測試技術，科研人員已在離心泵的清水流場［7-8］以及基于玻璃球顆粒固液兩相流實驗測量方面［9-11］取得較大的進展。宋天雄［12］采用PIV及MATLAB軟件后處理技術，針對紙漿纖維懸浮液在磁力攪拌器流場中的流變特性開展實驗研究，但僅就懸浮液的渦量進行了分析。在紙漿泵方面主要集中在性能實驗研究，進而探討紙漿濃度及含氣率對泵性能的影響［13-14］。Vaezi等［15］研究了輸送麥秸和玉米秸顆粒對離心泵性能的影響，發現了顆粒獨特的摩擦損失行為（即減阻效應）而引起泵效率提高的現象。葉道星等［16-18］對中濃漿泵及湍流發生器開展了一系列性能研究，這些研究都由于漿液不透明等因素，無法進一步通過PIV實驗對內流場進行可視化研究。在紙漿流場研究中，由于紙漿到達一定濃度時會纏結成絮狀及紙漿纖維膨脹的特性，會導致激光無法穿透，給PIV實驗的追蹤測量增加難度。隨著時間推移，紙漿會產生發酵、變質等現象，不易儲存。并且若要控制紙漿纖維尺寸相對均勻，需對其進行機械或化學處理，工作量大且會增加新的化學元素。自然材料中，頭發纖維的性質很接近紙漿纖維，其具有較好的懸浮性且對激光不敏感，因此本研究采用頭發纖維作為實驗材料以便進行后續的PIV可視化研究。

本研究以頭發纖維懸浮液作為研究對象，針對離心泵內細長纖維懸浮液內流場，首次采用CFD-DEM耦合計算及PIV可視化實驗相結合的方法，研究了三葉片開式葉輪離心泵內不同工況、不同軸截面以及不同質量濃度的長纖維懸浮液的內流場流動特性。通過流場的對比分析，探討長纖維對流態畸變的作用；同時，也對泵外特性進行了測量和分析。

1 耦合仿真分析

CFD-DEM數值模型的耦合計算包括兩個步驟，首先，通過EDEM軟件最大程度地反映顆粒特征，模擬工程中纖維顆粒真實形態，并生成研究所需的顆粒群；然后，通過與FLUENT耦合對泵內流場進行CFD模擬，建立纖維顆粒與流體間的信息傳遞。其中，EDEM軟件負責纖維顆粒的運動計算，流場通過FLU?ENT軟件進行耦合求解。

1.1 泵基本參數

本研究采用的是一臺扭曲葉片開式葉輪紙漿泵，葉輪葉片數為3，模型泵的主要幾何參數為：葉輪進口直徑D1=70 mm，出口直徑D2=184 mm，葉輪出口寬度b2=17.4 mm，葉片出口角β2=43°，蝸殼基圓直徑D3=191 mm，蝸殼出口直徑D4=70 mm，其三維模型如圖1所示。本研究主要探究離心泵葉輪流道內纖維懸浮液的運動情況，從而在模擬中忽略間隙的影響。

1.2 顆粒模型的建立

顆粒形狀與物料屬性可通過EDEM軟件建模了解。為盡可能模擬頭發纖維，采用多個圓球顆粒相連成鏈狀結構，如圖2所示。單根頭發纖維的長度為6 mm，質量為4.60×10-7kg。

圖2 物料顆粒示意圖Fig.2 Schematic diagramof solid material particles

1.3 模型控制方程

在本研究的流體模擬計算中，所有計算區域的流體相均在絕對坐標系內。流體的控制方程如下：

式中，ρf、u、t、P、g、μf分別表示流體的密度、流體速度、時間、壓力、重力加速度和動力黏度，α表示流體相的體積分數，F表示其他附加力的合力。

本研究中，EDEM軟件模擬采用基于歐拉-拉格朗日法追蹤顆粒的受力與運動情況，顆粒體積分數均小于10%。顆粒受力情況求解的表達式如下：

1.4 參數設置

在耦合模擬前，需對EDEM軟件和FLUENT軟件做相關參數設置。本研究利用FLUENT求解流體運動，采用歐拉-拉格朗日法對泵內固液兩相流進行數值模擬。考慮到流體與顆粒間的相互作用，采用雙向耦合的瞬態模擬。在離心泵入口處定義進口速度為2.03 m s以模擬額定工況下的流體運動情況，入口處湍流強度和湍流黏性比分別為5%和10。出口設置為壓力出口，其靜壓值設置為一個標準大氣壓。在葉輪入口前和蝸殼出口之后各延伸一段長度為3倍管徑的直管以保證顆粒與流體充分混合運動。計算中假設壁面為無滑移壁面，不考慮粗糙度對顆粒的影響。考慮到渦流及旋轉的因素，流體計算模型選擇RNG k-ε湍流模型，選取SIMPLEC算法來耦合泵內流場的速度和壓力，而動量、壓力、湍流動能和耗散率的離散化求解均采用二階迎風格式。

本研究利用EDEM軟件求解顆粒相的受力情況，顆粒間的接觸選用Hertz-Mindlin（無滑移）連接模型。在EDEM中，泵表面均被定義為有機玻璃，其密度為1200 kg/m3，泊松比為0.4066，剪切模量為1.41 GPa，楊氏模量為3.97 GPa，用來模擬頭發纖維的顆粒密度為680 kg/m3，泊松比為0.3。葉輪轉速設置為1450 r/min，與FLUENT設置保持一致。在直管進口處設置顆粒工廠，用于計算的顆粒均勻且隨機地從直管入口處跟隨流體一起運動。定義每秒入射20239、30344及40478個顆粒以表示顆粒質量濃度分別為0.1%、0.15%和0.2%。

通過建立耦合接口傳遞顆粒與流體的計算數據，FLUENT計算時間步長設為5×10-5s，即葉輪旋轉0.435°計算1次，最大迭代步數60，總步數20000。而EDEM中時間步長應小于FLUENT中時間步長且與FLUENT中時間步長成倍數關系，因此，EDEM中時間步長設置為1×10-6s。葉輪每旋轉87°，數據保存1次，總時間0.8 s，即葉輪旋轉19圈。

1.5 無關性驗證

計算網格通過ICEM劃分工具創建，考慮到非結構網格對于耦合計算適應性較好等因素，本研究采用四面體非結構網格對計算區域進行網格劃分。為驗證網格數對計算結果的影響，采用3組不同網格數方案進行計算，結果如表1所示。通過比較，且考慮到計算精度和計算資源的情況下，選用網格B進行后續數值計算。

表1 網格無關性模擬結果Table 1 Simulation results of grid independence

2 PIV實驗研究

2.1 實驗模型

為驗證上述仿真模型的正確性并進一步分析流場，開展了相應的PIV實驗研究。本研究采用的實驗泵是一臺典型的開式葉輪離心式紙漿泵，其主要部件為葉輪、蝸殼、耐磨擋板、軸和密封系統等。為便于進行PIV拍攝，采用半螺旋吸水室結構，而電機軸穿過吸水室與葉輪直連。蝸殼與半螺旋吸入室之間通過4個夾具固定。實驗泵整體結構如圖3所示，泵主要參數見1.1節。

考慮到有機玻璃透光性好，拋光處理后粗糙度較低等特點，模型泵中蝸殼和葉輪均由有機玻璃制造而成。前擋板使用不銹鋼制造，為減小其對PIV激光造成的反射，將前擋板進行涂黑處理。實物圖如圖4所示。

2.2 實驗裝置及實驗介質

為開展相關測試工作，搭建了一座開式實驗臺，主要由模型泵、進出口閥門、罐體和管道等組成。采用敞口塑料制罐，罐內采用攪拌器連續攪拌，以便于顆粒均勻混合。管路、連接彎頭等全部采用DN80不銹鋼管制作。

圖3 模型泵結構圖及零部件Fig.3 Structural drawingand part drawingof themodel pump

圖4 模型泵蝸殼、葉輪和前擋板Fig.4 The volute,impeller，and front plate of the model pump

水泵實驗臺的測試精度通常用水泵效率測量的不確定度來評估。泵效率測量不確定度（eη）計算公式如下：

式中，eQ為流量測量的不確定度，eH為揚程測量的不確定度，eP為驅動機輸入功率測量的不確定度。

實驗中使用的DN80流量計經過法定計量部門的檢定，符合0.2級基本精度要求，流量測量的不確定度為：

實驗中使用的差壓變送器經過法定計量部門的檢定，則揚程H測量的不確定度為：

實驗中使用的電流互感器精度等級為0.2級，即線性誤差E=0.2%，則驅動機輸入功率測量的不確定度為：

因此，泵效率測量的不確定度為：

綜上，該實驗臺精度能夠滿足GB/T 3216—2016標準規定的2級精度要求。

PIV實驗裝置由PIV測試系統及外觸發同步控制系統兩部分構成。實驗中采用美國TSI公司的粒子圖像測速系統。該測試系統主要包括：YAG200-NWL型脈沖激光器、610015-SOL型光臂及其片光源透鏡系統、610035型同步器、630059POWERVIEW型CCD相機。采用Insight 3G和Tecplot軟件對PIV實驗圖像進行采集處理。

PIV實驗采用了具有較好跟隨性和散射型的二氧化硅空心玻璃球作為示蹤粒子，該粒子直徑為20μm，材質密度為1.05 g/cm3。PIV實驗采取相平均的方法進行，即在每個相位點拍攝多組圖片，使用Insight 3G軟件將所拍攝的圖片進行平均處理，得到對應的流場信息。本研究通過軸編碼器和外觸發同步控制系統控制相機在葉輪的葉片掠過蝸殼隔舌位置時進行拍攝。在頭發纖維懸浮液PIV實驗中，示蹤粒子有一定幾率被頭發顆粒遮擋，繼而影響后處理中示蹤粒子的識別，因此實驗中每個相位點拍攝200張圖片以確保實驗結果的準確性。

目前，并沒有專門的設備來加工纖維顆粒以控制其長度完全一致，所以本實驗通過人工篩選處理以確保頭發纖維長度嚴格控制在2～8 mm，制取難度及工作量較大，因此，本研究對質量濃度分別為0.1%、0.15%、0.2%及0.3%的頭發纖維懸浮液進行研究，暫不研究更高質量濃度的頭發纖維懸浮液。

2.3 實驗方案與測量工況

實驗時，啟動并調節變頻器使得離心泵轉速穩定在1450 r min。為驗證實驗臺和實驗方法的可靠性，需對實驗臺進行重復外特性實驗。通過泵參數測量儀測量參數并求得離心泵的揚程及效率。

使用電子秤測量頭發纖維質量，逐步將對應頭發纖維加入罐中并開啟攪拌器，使頭發纖維在罐中分布均勻，進行離心泵葉輪內部流場的PIV測量，分別對0.6Qd、0.8Qd、1.0Qd、1.2Qd、1.4Qd5個工況點進行PIV實驗，其中，Qd為實驗泵額定工況，大小為33.5 m3/h。圖5為PIV實驗激光截面效果圖。由圖5可知，PIV拍攝效果較好，頭發纖維的存在并沒有遮擋示蹤粒子，因此不需要對固相顆粒圖像和液相示蹤粒子圖像進行分離。

實驗泵采用了空間扭曲葉片，葉輪不同軸向位置將對葉輪流道的流場產生影響。為了研究扭曲葉片對葉輪流道內流場的影響，PIV實驗選取葉輪區域內3個截面進行拍攝（見圖6），這3個截面都垂直于泵軸，分別記為截面Z1、Z2、Z3，其中截面Z2位于葉輪出口寬度中間位置，截面Z1、Z3各距Z2截面6 mm。為表述方便，將隔舌位置處葉片記為葉片1，沿順時針方向的第二、第三個葉片記為葉片2、葉片3。記葉片1左邊流道為L1，流道L2、L3沿順時針方向類推。

圖5 PIV實驗激光截面效果圖Fig.5 Effect drawing of PIV test laser section

圖6 PIV拍攝截面圖Fig.6 Cross-section schematic of PIVmeasurement

3 結果與討論

3.1 離心泵的性能

為得到準確的離心泵外特性結果，對實驗泵的揚程、流量進行了重復性實驗，頭發纖維質量濃度為0.1%時，離心泵的外特性曲線如圖7所示。從圖7可以看出，3組外特性實驗數據變化較小，曲線基本重合，最大誤差在4%以內，證明了實驗臺及實驗結果的可靠性。

圖8為清水、0.1%、0.15%及0.2%頭發纖維質量濃度條件下離心泵性能曲線對比圖。從圖8可以看出，不同頭發纖維質量濃度下，泵揚程變化不明顯，但都高于輸送清水時的揚程。且隨著泵內頭發纖維質量濃度的提高，泵效率逐漸增大，纖維質量濃度達到0.1%及0.15%時已初具減阻趨勢，纖維質量濃度達到0.2%時已完全顯示出泵內因纖維顆粒存在而引起的減阻效應。

圖8 不同纖維質量濃度離心泵性能曲線對比圖Fig.8 Comparison of pumping performance curvesof different fiber concentrations

3.2 泵內流場驗證分析

為了便于分析流體相對于葉輪的運動規律，采用自帶的PIV速度處理軟件將處理后的絕對速度矢量分解，從而得到葉輪流道內流體的相對速度。在進行PIV實驗時，從泵軸左側打入激光，由于泵軸會遮去部分激光，所以泵軸右側沒有激光，因此在后處理時略去這一區域。

圖9（a）顯示了通過CFD-DEM耦合計算，相同質量濃度下葉輪轉到PIV拍攝相同位置時，Z2截面處流體的絕對速度及相對速度分布；圖9（b）給出了通過PIV拍攝得到頭發纖維質量濃度為0.1%時，葉輪區域Z2截面處流體的絕對速度及相對速度分布。從圖9可以看出，靠近葉片工作面處流體的絕對速度高于背面，從葉輪進口到出口，Z2截面處流體的絕對速度逐漸增大，葉輪出口附近Z2截面處流體的絕對速度有所下降。而Z2截面處流體的相對速度會從葉輪進口至出口處逐漸增大，且在葉輪流道中段靠近壓力面附近可觀察到低速區域，在此區域內流體的流速低于周圍區域內的流速，但差距并不明顯。本研究選取相同位置葉片附近流體速度進行分析，并將CFD-DEM耦合計算結果與PIV實驗結果進行對比發現，使用CFD-DEM方法得到的流體速度分布及變化趨勢與PIV實驗結果吻合較好，說明CFD-DEM方法能夠較好地模擬泵內復雜流場下的固液兩相流運動，驗證了采用CFD-DEM方法定性分析泵內流場的可行性，同時反映了PIV實驗的準確性。

由圖9可以看出，由于遮光、反光等因素影響，L1流道的PIV拍攝效果明顯優于L2、L3，所以后續采用L1流道的PIV實驗結果對泵內流場進行分析。關于離心泵內清水流場PIV實驗，這方面已有許多的研究結論，本文不再贅述。

3.3 測試結果及分析

3.3.1 不同流量下葉輪流道內流體的相對速度分布

為更具體地分析葉輪流道內流體的相對速度隨流量的變化情況，分別提取出頭發纖維質量濃度為0.1%時，0.6Qd、1.0Qd和1.4Qd工況下Z2截面L1流道內葉片壓力面、中間流線及吸力面上流體的相對速度分布并繪制曲線圖，結果如圖10～圖12所示。其中，橫坐標L表示距離葉輪進出口的位置，記r0為當前位置半徑，r2為葉輪出口半徑，0＜r0/r2＜0.6為葉輪輪轂區域，0.6＜r0/r2＜1為葉輪流道區域，因此，L=0代表在r0/r2=0.6位置處，L=1代表在r0/r2=1位置處。

由圖10可以看出，在0.6Qd和1.0Qd工況下，葉片壓力面附近流體的相對速度隨著L的增大呈先減小后增大的趨勢，即在流道中段壓力面附近會產生低速區，且隨著流量逐漸增大，低速區強度和面積逐漸減小，流量增加到1.4Qd時，低速區基本消失。由圖11可知，在0.6Qd和1.0Qd工況下，中間流線流體的速度變化規律與壓力面一致，說明葉輪流道中段亦存在低速區，但低速區強度和面積小于葉片壓力面。由圖12可知，在葉片吸力面附近，只有在小流量工況下會形成低速區。葉輪流道內，流體的相對速度總體隨著流量的增大而增大。綜上，當流量較小時，葉輪流道壓力面及中間流道附近流體會形成低速區，隨著流量增大，低速區的面積和強度會逐漸降低，中間流道內流體的低速區會逐漸消失。

3.3.2 額定流量下不同軸截面流體相對速度分布

本研究采用扭曲葉片葉輪，對前文提到的Z1、Z2和Z3截面進行PIV實驗分析。由于不同軸向位置葉片的扭曲程度不同，因此對額定工況下L1流道內頭發纖維質量濃度0.1%時，r0/r2=0.7、0.85及1的3個不同位置處從葉片吸力面到壓力面流體的相對速度分布進行分析，結果如圖13～圖15所示。

圖9 Z2截面處流體的速度分布Fig.9 Velocity distributions of fluids in Z2 cross-section

圖10 不同流量下壓力面流體速度分布Fig.10 Pressuresurfacevelocity distributionsof fluids at different flowrates

圖11 不同流量下中間流線流體速度分布Fig.11 Intermediate flow channel velocity distri?butionsof fluidsat different flowrates

圖12 不同流量下吸力面流體速度分布Fig.12 Suction surface velocity distributions of fluidsat different flowrates

對比圖13～圖15可知，在r0/r2=0.7處，葉片扭曲程度相對較大，因此各個截面上流體的相對速度相差較大且靠近壓力面流體的速度變化幅度較大。而隨著位置向出口方向移動，葉片的扭曲程度逐漸降低，各截面上流體的相對速度差值逐漸減小，到r0/r2=1即葉片出口位置處，葉片形狀可近似于圓柱葉片，各截面流體的速度變化相對平滑，基本歸于一致。

3.3.3 額定流量下不同頭發纖維質量濃度Z2截面處流體的相對速度分布

為進一步研究纖維質量濃度對于紙漿泵內液相流場的影響，限于篇幅影響，本研究分別提取了Z2截面上L1流道內質量濃度0.1%和0.2%頭發纖維懸浮液在0.6Qd、1.0Qd及1.4Qd工況和r0/r2=0.7及r0/r2=1情況下，沿吸力面至壓力面的相對速度分布圖（見圖16和圖17）。

圖13 不同截面上r0/r2=0.7處流體的相對速度分布Fig.13 Relativevelocity distributionsof fluidsat r0/r2=0.7 and different cross-sections

圖14 不同截面上r0/r2=0.85處流體的相對速度分布Fig.14 Relativevelocity distribution of fluidsatr0/r2=0.85 and different cross-sections

圖15 不同截面上r0/r2=1處流體相對速度分布Fig.15 Relative velocity distributions of fluids atr0/r2=1 and different cross-sections

圖16 r0/r2=0.7處不同流量時流體的相對速度分布Fig.16 Relative velocity distributions of fluids at r0/r2=0.75 and different flowretes

由圖16及圖17可知，隨著流量逐漸增大，高質量濃度纖維懸浮液Z2截面內的相對速度會逐漸超過低質量濃度纖維懸浮液、各位置相對速度分布趨勢相似度增大，且隨著位置靠近葉輪出口區域，上述變化趨勢加快。大流量工況下，各質量濃度纖維懸浮液在葉輪流道內從吸力面至壓力面的相對速度總體呈下降趨勢。且隨著質量濃度增大，葉輪出口截面處流場的均勻性在各流量工況均有所改善，葉輪流道內各位置處流體的相對速度變化更趨平滑，有利于降低葉輪出口滑移率，提高葉輪效率。

圖17 r0/r2=1處不同流量時流體的相對速度分布Fig.17 Relative velocity distributions of fluidsat r0/r2=1 and different flowrates

4 結論

本研究使用CFD-DEM耦合方法對離心泵內纖維顆粒懸浮液流場進行模擬和分析，發現其與PIV實驗所得葉輪流道內速度場的分布趨勢基本吻合，從而驗證了CFD-DEM方法初步模擬研究纖維懸浮液內流場的可行性；進一步采用PIV實驗結果對泵內流場進行具體分析，主要結論如下。

（1）隨著泵內纖維質量濃度逐漸增大，泵效率隨之提高，纖維質量濃度達到0.2%時已明顯呈現出減阻效應。

（2）當輸送相同質量濃度纖維時，大流量工況下葉輪流道內流體的流動穩定，未出現漩渦和低速區。而在小流量工況下，沿葉輪流道中段壓力面附近流體會產生低速區，且低速區隨著流量的減小逐漸向中間流道延伸。

（3）在垂直于泵軸的不同截面葉輪流道內，受葉片扭曲程度的影響，各截面處流體的相對速度存在差異。靠近輪轂側，葉片扭曲程度較大時，各位置處流體相對速度分布相差較大。而隨著位置向出口方向移動，葉片的扭曲程度逐漸降低，各截面上流體的相對速度差值逐漸減小，在靠近葉輪出口處各截面處流體相對速度變化相對平滑且基本歸于一致。

（4）頭發纖維質量濃度的變化會影響流道內流體的相對速度分布，當流量逐漸增大時，較高質量濃度纖維懸浮液內液相相對速度大小會逐漸超過低質量濃度纖維懸浮液，在大流量工況下，各質量濃度纖維懸浮液在葉輪流道內從吸力面至壓力面的相對速度總體呈下降趨勢。而隨著纖維懸浮液質量濃度增大，葉輪出口截面流場的均勻性在各流量工況均有所改善，葉輪流道內各位置處流體的相對速度變化更趨平滑。上述結果可對理解長纖維懸浮液減阻效應提供參考。