基于沖蝕耦合動網(wǎng)格方法的搖臂式噴灌噴頭沖蝕變形數(shù)值模擬研究

申恒峰，黃劍峰，周傳亮

(北部灣大學 建筑工程學院，廣西 欽州 535011)

1 研究背景

在使用渾水進行噴灌的過程中，噴頭由于持續(xù)地被沖蝕磨損會發(fā)生沖蝕變形，進而影響噴頭的水力性能。目前對于噴頭使用過程中各個階段的沖蝕變形情況，以及沖蝕變形對噴頭水力性能的具體影響程度的研究還不夠全面。因此，本文通過探究噴灌噴頭使用160 h內(nèi)出口直徑與出口流量之間的變化關系、出口直徑與噴頭射程之間的變化關系以及出口直徑、出口斷面形狀系數(shù)與末端水滴之間的變化關系開展數(shù)值模擬研究。

我國對噴灌噴頭的研究起步相對較晚[1]，諸多研究主要集中在噴頭自身結(jié)構(gòu)的優(yōu)化[2-3]和組合模式創(chuàng)新[4]方面，近年來，噴頭沖蝕變形問題也逐漸得到關注。為了探究泥沙對噴頭的沖蝕磨損規(guī)律，李英能等[5]通過開展4組不同含沙量水流對噴頭沖蝕磨損試驗研究，分析了噴頭流量、噴嘴磨損量和噴嘴尺寸與含沙量及沖蝕時間之間的關系；李紅等[6]通過研究泥沙磨損對噴頭可靠性的影響，得到了渾水質(zhì)量濃度上升會使得噴頭諸多水力性能變化率增加的結(jié)論，并認為提高噴嘴的耐磨性會使噴頭的使用壽命延長；翟恩昱等[7]通過對固體顆粒做受力分析建立了運動方程，確定了固體顆粒對塑性材料噴頭內(nèi)表面沖蝕磨損時的速度。

鑒于開展沖蝕變形的物理試驗研究較為困難，部分學者進行了相關的數(shù)值模擬分析。李文穎[8]、王子君[9]開展了8034D型噴頭的沖蝕磨損試驗研究和數(shù)值模擬，并進行了結(jié)果對比，但數(shù)值模擬并未考慮噴頭的全流道流場，也未考慮噴頭的沖蝕變形；Xu等[10]通過編寫網(wǎng)格重塑程序，利用Fluent計算的沖蝕率數(shù)據(jù)更新壁面網(wǎng)格坐標，探究了噴頭的沖蝕變形速率，發(fā)現(xiàn)時間步長對計算結(jié)果有較大影響。然而編寫網(wǎng)格重塑程序一直是沖蝕變形研究的難點，于是石油管道領域的學者利用沖蝕耦合動網(wǎng)格方法進行構(gòu)件沖蝕變形研究。Adedeji等[11]基于Oka沖蝕模型，結(jié)合沖蝕耦合動網(wǎng)格方法進行了彎管沖蝕變形研究，與試驗結(jié)果對比后發(fā)現(xiàn)該方法可以準確預測構(gòu)件的沖蝕變形；Wang等[12]通過將沖蝕耦合動網(wǎng)格方法應用到兩個工程案例中，發(fā)現(xiàn)該方法在針對構(gòu)件優(yōu)化和預測構(gòu)件壽命上具有可行性；孫巖等[13]利用該方法對繞絲篩管進行了沖蝕變形研究，分析了繞絲沖蝕變形發(fā)生的位置，并得到了不同沖蝕歷時下繞絲的沖蝕變形情況。

綜上所述，沖蝕耦合動網(wǎng)格方法的提出為開展噴頭沖蝕變形的數(shù)值模擬研究提供了一種新思路，利用該方法研究噴頭沖蝕變形規(guī)律及變形后對其水力性能的影響，可以為預測噴頭使用壽命和噴頭結(jié)構(gòu)優(yōu)化等提供參考依據(jù)。

2 研究方法

2.1 沖蝕耦合動網(wǎng)格方法

沖蝕耦合動網(wǎng)格是一種基于準穩(wěn)態(tài)運行，并實現(xiàn)固-液-沖蝕完全耦合的構(gòu)件變形計算方法。Fluent提供的沖蝕模型能夠計算初始沖蝕率，目前已有諸多基于此模型開展的沖蝕磨損研究[14-15]。但沖蝕率實質(zhì)上是瞬態(tài)的，因為沖蝕磨損是一個持續(xù)的動態(tài)過程，所以預測較長時間后的沖蝕率需要考慮構(gòu)件的沖蝕變形。沖蝕耦合動網(wǎng)格方法是通過將總沖蝕時間拆分，并在每一段時間上計算出壁面沖蝕率，然后據(jù)此更新壁面網(wǎng)格坐標，而壁面網(wǎng)格坐標則通過動態(tài)網(wǎng)格子模型基于物理時間步長實現(xiàn)更新。

耦合動網(wǎng)格方法首先需要在穩(wěn)態(tài)條件下計算得到噴頭的壁面沖蝕率等數(shù)據(jù)，在開始動態(tài)網(wǎng)格沖蝕計算后，網(wǎng)格的坐標位置會根據(jù)設定的時間步長內(nèi)噴頭壁面沖蝕率進行更新，其計算流程如圖1所示，單個單元面的變形量由以下公式計算：

圖1 噴頭變形量計算流程

(1)

式中：Δxface網(wǎng)格單元變形量，m；Δtmm為網(wǎng)格運動時間步長，s；Eface為特定時間步長的壁面沖蝕率，kg/(m2·s)；Δtmm為計算時段，s；ρwm為被沖蝕靶材的密度，kg/m3。

2.2 湍流模型及沖蝕模型

利用標準k-ε湍流模型對噴頭不同工況進行模擬，該模型適用性較好[16]，由湍流脈動能k和擴散率ε組成，在實際工程中應用廣泛，其準確性得到了較好的驗證。選用DPM(deformable part model)作為離散相模型[17]，其提供了Generic、Mclaury、Oka和Dnv等沖蝕模型，選用其中適用范圍較廣的Generic模型，該模型由Huser等[18]提出，沖蝕率被定義為：

(2)

式中：mp為顆粒質(zhì)量流量，kg/s；C(dp)為顆粒粒徑函數(shù)；α為顆粒路徑與壁面的沖擊角，(°)；f(α)為沖擊角函數(shù)；v為顆粒相對速度，m/s；b(v)為顆粒相對速度函數(shù)；Aface為壁面面積，m2。

2.3 求解方法與邊界條件

研究基于有限體積法對控制方程進行離散，利用SIMPLEC算法對速度和壓力進行耦合求解，擴散項采用二階中心差分格式，對流項采用二階迎風格式。分別設置壓力入口和出口，入口壓力設為350 kPa。設置連續(xù)相為清水，其密度為998 kg/m3，顆粒密度為2 650 kg/m3。

2.4 幾何模型與網(wǎng)格劃分

首先對PY1-20sh型噴頭進行全流道三維建模，模型尺寸比例為1 ∶1，具體模型幾何參數(shù)如圖2(a)所示。由于DPM模型的局限性，無法計算壁面邊界處的沖蝕率，使得正常情況下無法預測出口直徑的變化，因此需考慮噴嘴約為0.2 mm的實際壁面厚度，且仍以原始出口斷面面積計算出口當量直徑。

不考慮噴頭內(nèi)螺紋等細部結(jié)構(gòu)，假設各部分連續(xù)光滑，利用ICEM軟件對幾何模型進行六面體非結(jié)構(gòu)化網(wǎng)格劃分，如圖2(b)所示。通過多次網(wǎng)格劃分并以相同初始條件計算出口流量進行網(wǎng)格密度分析，計算結(jié)果見表1，考慮在保證精度的前提下節(jié)約計算資源，最終取網(wǎng)格數(shù)為91 808。

圖2 噴頭模型幾何參數(shù)及三維模型網(wǎng)格劃分(單位：mm)

表1 網(wǎng)格密度分析

3 沖蝕條件設置與理論計算公式

3.1 顆粒屬性與時間步長設置

參考林管[19]對PY1-20sh型噴頭開展的沖蝕磨損試驗研究，本文設置與其一致的初始條件，并通過結(jié)果對比驗證本文數(shù)值模擬結(jié)果的可靠性，表2為具體的顆粒注入?yún)?shù)。

表2 噴頭沖蝕磨損數(shù)值模型顆粒注入?yún)?shù)

文獻[19]的試驗過程中泥沙顆粒不斷發(fā)生碰撞破碎，使顆粒級配始終在發(fā)生變化，表2中的顆粒級配為綜合考慮后的取值。同時設置顆粒由噴頭入口注入，顆粒速度與水流流速一致，水流含沙量設3組數(shù)據(jù)，分別為2.5、5.0和7.5 kg/m3。

選用固定時間步長進行沖蝕變形計算，并與自適應時間步長的計算結(jié)果進行對比。通過分析噴頭累計平均變形量來選擇合理的時間步開展進一步分析，表3為各時間步長計算的累計平均變形量。計算中設置總沖蝕歷時為160 h，以40 h為時間間隔，沖蝕耦合動網(wǎng)格平滑步數(shù)為5步，每個流動迭代次數(shù)為100次。

表3 各時間步長變形量計算結(jié)果10-6m

由表3可以看出，各時間步長計算的0～160 h累計平均變形量結(jié)果比較接近，只有時間步長為2 h時離散度相對較大，時間步長為5 h時離散度相對最小，因此選用時間步長為5 h。

3.2 主噴嘴出口直徑及直徑增加率計算

通過出口斷面面積計算出口當量直徑：

(3)

式中：D為主噴嘴出口當量直徑，m；A為主噴嘴出口斷面面積，m2(該斷面位置的選取見后文圖6)。

直徑增加率計算公式為：

ηD=(Dt-D0)/D0×100%

(4)

式中：ηD為直徑增加率；Dt和D0分別為t時刻與初始時刻主噴嘴出口直徑，m。

3.3 質(zhì)量損失計算

由內(nèi)流道體積變化計算質(zhì)量損失：

Δm=ρ(V1-V0)

(5)

式中：ρ為壁面材料密度，該噴頭壁面材料為鑄鋁，密度為2 719 kg/m3；V0和V1分別為沖蝕前、后的噴頭內(nèi)流場體積，m3。

3.4 流量增加率及主噴嘴流量系數(shù)計算

流量增加率計算公式為：

ηQ=(Qt-Q0)/Q0×100%

(6)

式中：ηQ為流量增加率，%；Qt和Q0分別為t時刻和初始時刻噴嘴出口的總流量，m3/h。

流量系數(shù)計算公式為：

(7)

式中：μ為主噴嘴流量系數(shù)；Qm為主噴嘴流量，m3/h；A為主噴嘴出口斷面面積，m2；h為噴頭壓力水頭，m。

3.5 噴頭射程理論計算

噴頭射程的理論計算模型較多[20-22]，馮傳達[23]基于對PY1-20sh型噴頭開展的試驗研究提出了射程計算公式，因與本文研究噴頭類型一致，所以有較好的適用性，該公式為：

(8)

式中：R為噴頭射程，m；μ為主噴嘴流量系數(shù)；H0為主噴嘴出口壓力水頭，m；θ為噴頭仰角，(°)；D為主噴嘴當量直徑，mm。

3.6 形狀系數(shù)與末端水滴直徑理論計算

噴嘴的形狀常用形狀系數(shù)β進行描述，β值越大說明噴嘴出口斷面越偏離圓形，對末端水滴直徑影響越大[24]，其計算公式為：

(9)

目前，部分學者已進行了噴頭水滴分布的理論研究[25-26]，其中比較全面地考慮各影響因素的李久生[27]公式為：

(10)

式中：dm為末端水滴直徑，m；D為主噴嘴當量直徑，mm；h為噴頭工作壓力水頭，m；β為主噴嘴形狀系數(shù)；μ為主噴嘴流量系數(shù)；R為射程，m。

4 結(jié)果與分析

4.1 噴頭沖蝕磨損部位分析

進行噴頭全流道計算模型的建立可以更真實地反映噴頭內(nèi)部流體流動情況，獲得更準確的顆粒運動軌跡。但由于噴頭整體結(jié)構(gòu)相對較大且沖蝕磨損分布并不均勻，不利于分析沖蝕變形相對嚴重的局部位置，又因為沖蝕變形的計算基于沖蝕率基礎數(shù)據(jù)，因此首先進行沖蝕率分布分析以確定噴頭主要沖蝕磨損部位，從而進一步探究該部位的沖蝕變形情況。圖3為噴頭內(nèi)流場中沿流程方向的DPM顆粒濃度和沖蝕率分布圖。

圖3 噴頭內(nèi)沿流程方向DPM濃度與沖蝕率分布

由圖3可以看到，噴頭內(nèi)流場中DPM濃度隨著噴頭管徑的縮小呈微弱增大的趨勢，在兩個噴嘴處濃度較大。通常情況下DPM濃度與壁面碰撞率呈正相關關系，顆粒碰撞次數(shù)越多則沖蝕率相對越大。但沖蝕磨損的產(chǎn)生是諸多因素共同作用的結(jié)果，可以看到?jīng)_蝕率從副噴嘴附近開始增大，主噴嘴附近沖蝕率最大。為了更準確地解釋上述情況需要結(jié)合噴頭內(nèi)的流場情況進行分析，圖4為噴頭內(nèi)流場中的流速和壓力分布。

由圖4可知，出現(xiàn)上述情況是因為噴頭入口附近速度較小，僅約為2 m/s，水流對泥沙曳力較小，此時泥沙動能較小，所以該段沒有發(fā)生較大沖蝕磨損。主噴嘴附近由于管徑大幅收縮，流速突增，壓力減小，使得水流流態(tài)紊亂，增加了泥沙與壁面碰撞概率，導致該處沖蝕磨損嚴重。

圖4 噴頭內(nèi)的速度與壓力分布云圖

為了分析在水流不同含沙量情況下噴頭沖蝕率分布情況，選擇2.5、5.0、7.5 kg/m33組含沙量工況進行研究，3組含沙量下主噴嘴沖蝕率分布如圖5所示。由圖5可見，主噴嘴的沖蝕程度隨著含沙量的增大而逐步加深，但主噴嘴整體的沖蝕率分布未發(fā)生明顯變化，始終主要分布于主噴嘴收縮段。

圖5 水流不同含沙量時主噴嘴沖蝕率分布

4.2 主噴嘴沖蝕變形分析

通過分析噴頭沖蝕率分布可知沖蝕磨損主要集中在主噴嘴處。因此，為了更清晰地觀察噴頭局部沖蝕變形情況，需分別探究主噴嘴在3組含沙量工況下不同沖蝕歷時的累計平均變形量變化規(guī)律，并結(jié)合不同時間節(jié)點的沖蝕率變化進行對比分析。圖6為噴頭在含沙量7.5 kg/m3下沖蝕160 h的主噴嘴累計變形量分布云圖，圖7為水流不同含沙量時不同沖蝕時段主噴嘴累計平均沖蝕變形總量和增量以及沖蝕率隨歷時變化圖。

由圖6可以看出，主噴嘴網(wǎng)格沖蝕變形程度與沖蝕變形量分布十分吻合，出口收縮段發(fā)生了較大的變形，該變形表現(xiàn)為主噴嘴壁面整體向外不均勻膨脹擴張，使得噴嘴內(nèi)流道體積增大。出現(xiàn)該現(xiàn)象的原因是由于收縮段增大了顆粒與壁面的碰撞角度，使顆粒趨向于正面碰撞而加深了壁面的沖蝕程度，流道體積擴大的根本原因是壁面材料從表面逐漸剝離。

由圖7可知，主噴嘴平均累計變形量隨著水流含沙量的增大而逐漸增大。同時，在160 h的沖蝕過程中，不同時間節(jié)點的平均沖蝕率逐漸減小。這可能是因為沖蝕變形會改變主噴嘴壁面結(jié)構(gòu)，而變形后的壁面結(jié)構(gòu)會避免與一部分泥沙顆粒的碰撞，從而減緩了該區(qū)域后續(xù)變形的發(fā)生。因此主噴嘴每40 h的累計平均變形量表現(xiàn)為逐漸減小，且隨著水流含沙量的增大這種變化規(guī)律愈加明顯。3組含沙量水流沖蝕160 h后主噴嘴累計平均變形量分別約為5.6×10-6、1.0×10-5和1.3×10-5m。

圖7 水流不同含沙量時主噴嘴變形情況隨沖蝕歷時的變化

4.3 噴頭總質(zhì)量、出口直徑與出口流量變化分析

沖蝕變形的本質(zhì)是壁面材料逐漸從表面剝離，可以通過計算噴頭質(zhì)量損失進行該過程的分析。同時，為了進一步定量分析主噴嘴的沖蝕變形，需要計算主噴嘴出口直徑增加率，而出口直徑的擴大將對出口流量造成影響，通常情況會使出口流量增大，諸多研究表明大部分型號的噴頭出口流量變化率超過10%便達到報廢的條件。表4為160 h內(nèi)3組含沙量水流不同沖蝕歷時的噴頭出口直徑、質(zhì)量及出口流量的試驗與模擬計算結(jié)果。

表4 3組含沙量水流不同沖蝕歷時的噴頭出口直徑、質(zhì)量及出口流量的試驗與模擬計算結(jié)果

由表4可知，噴頭質(zhì)量損失總量與含沙量大小以及沖蝕時間均呈正相關，每40 h的質(zhì)量損失增量與含沙量大小呈正相關，與沖蝕時間無明顯線性關系，且質(zhì)量損失增量變化隨沖蝕時間的延長而逐漸趨于平穩(wěn)，這是由于噴頭后期的沖蝕磨損強度逐漸減弱所致。在3組含沙量水流沖蝕160 h后的噴頭質(zhì)量損失總量分別約為0.031、0.065和0.096 g，換算后可知其質(zhì)量損失速率分別約為0.19、0.41 與0.60 mg/h。表4的計算結(jié)果表明，噴嘴運行早期出口直徑增量較大，隨沖蝕歷時的延長逐漸趨于平穩(wěn)。林管[19]的試驗結(jié)果可能由于測量精度原因沒有記錄小于0.01 mm的直徑變化，使得部分時間段內(nèi)的直徑增量表現(xiàn)為零。3組含沙量水流沖蝕160 h后直徑增加率分別約為0.62%、0.98%和1.18%，略大于文獻[19]的試驗結(jié)果，但仍在合理范圍內(nèi)。噴嘴部分通常需要設計管徑收縮段以更好地控制水流，但也增加了該段的沖蝕磨損，而由于收縮段沖蝕變形導致的出口直徑增大將會造成噴頭諸多水力性能的改變，其中影響最大的是噴頭出口流量。

數(shù)值計算的噴頭初始流量略大于文獻[19]的試驗結(jié)果，但誤差約1.1%在合理范圍內(nèi)。通過分析表4的出口流量增加率可知，出口流量與出口直徑隨沖蝕歷時的變化趨勢具有一致性，在3組含沙量水流沖蝕160 h后噴頭流量增加率分別約為2.79%、2.91%和3.28%，每40 h的流量增量逐漸減小，數(shù)值計算結(jié)果與文獻[19]的試驗結(jié)果整體較為吻合。在工作壓力不變的情況下出口流量增加說明噴頭對水流的控制能力變差，將改變噴嘴出口斷面的水流平均速度大小，通過分析流量系數(shù)變化可以看到噴頭能量損失增加的過程。圖8為3組含沙量水流條件下主噴嘴流量系數(shù)μ隨沖蝕歷時變化的計算結(jié)果。

如圖8所示，主噴嘴的流量系數(shù)隨沖蝕歷時的增加而逐漸減小，在3組含沙量下運行160 h后流量系數(shù)分別減小約0.9%、1.6%和2.2%，流量系數(shù)綜合反映了出口直徑、工作壓力和出口流量之間的關系，該系數(shù)的減小說明噴頭的能量轉(zhuǎn)化效率降低。

圖8 水流不同含沙量時主噴嘴流量系數(shù)隨沖蝕歷時的變化曲線

4.4 噴頭射程與末端水滴直徑變化理論分析

圖9 水流不同含沙量時主噴嘴的射程、出口斷面濕周和形狀系數(shù)隨沖蝕歷時的變化曲線

由圖9(a)可知，噴頭射程隨沖蝕歷時的延長而逐漸縮短，經(jīng)過160 h沖蝕后，水流含沙量為7.5 kg/m3時的射程減幅最大，縮短了1.03 m，3組含沙量射程分別縮短約1.4%、2.6%和3.8%。噴頭射程的縮短會大幅減少原本灌溉范圍內(nèi)距離噴頭安放中心較遠區(qū)域的灑水量，其次由于噴嘴出口斷面流速梯度減小，會增加原本灌溉范圍內(nèi)中間區(qū)域的灑水量，減少噴頭安放中心附近的灑水量，從而影響灌溉質(zhì)量。

由圖9(b)可知，噴嘴出口斷面濕周隨沖蝕歷時的增加而逐漸變大，3組含沙量水流沖蝕160 h后出口斷面濕周分別增大了0.79%，1.35%和1.82%。斷面濕周的增大使得水流與噴嘴壁面的接觸面積增加，而出口平均流速、湍動能等物理量均受邊壁的影響較大，是造成水流能量損耗的主要原因。

由圖9(c)可知，主噴嘴出口形狀系數(shù)變化相對較小且隨沖蝕歷時的延長出現(xiàn)一定的波動，但整體上在逐漸增大，出現(xiàn)該變化是因為主噴嘴的沖蝕變形主要發(fā)生在局部位置，但隨沖蝕歷時的延長局部沖蝕位置發(fā)生了轉(zhuǎn)移，使噴頭出口斷面形狀獲得一定程度的復原。形狀系數(shù)是影響噴頭末端水滴直徑的重要因素，為了探究具體的影響需進行噴頭沖蝕160 h內(nèi)的末端水滴直徑理論計算，圖10為相應的計算結(jié)果。

圖10 水流不同含沙量時主噴嘴末端水滴直徑隨沖蝕歷時的變化曲線

由圖10可看出，噴頭末端水滴直徑隨沖蝕歷時的延長而逐漸減小，但減小幅度相對較小，3組含沙量水流沖蝕160 h后噴頭末端水滴直徑分別減小了0.6%、1.3%和1.9%，說明噴頭的沖蝕變形對其末端水滴直徑的影響十分有限。

由理論計算公式(公式(10))可以看出，噴頭末端水滴直徑受到結(jié)構(gòu)尺寸影響的因素分別為出口當量直徑D和出口斷面形狀系數(shù)β，其中形狀系數(shù)β對末端水滴直徑影響的比重較大，但由于通常情況下噴嘴出口斷面的形狀系數(shù)并不會發(fā)生較大變化，所以水滴直徑并未受到明顯的影響。另外，末端水滴直徑減小會增強噴頭的霧化效果，減小水滴對作物的打擊力度。

5 結(jié) 論

利用沖蝕耦合動網(wǎng)格方法對PY1-20sh型噴灌噴頭進行沖蝕變形計算，實現(xiàn)了固-液-沖蝕-網(wǎng)格運動的完全耦合，通過模擬噴頭在運行160 h內(nèi)的沖蝕變形發(fā)展過程，并利用變形后的網(wǎng)格文件進行噴頭水力性能變化的研究，得出以下主要結(jié)論：

(1) PY1-20sh噴頭在運行過程中，主噴嘴收縮段持續(xù)地受到?jīng)_蝕磨損，造成噴嘴流道體積和出口直徑增大以及出口形狀小幅度偏離圓形，這些變形會在不同程度上對噴頭的水力性能產(chǎn)生影響。其中出口直徑的變化是影響噴頭水力性能的主要原因，其對出口流量的影響十分明顯，也是導致噴頭報廢的主要因素。出口直徑和出口形狀系數(shù)的變化對噴頭射程和末端水滴直徑未造成明顯影響。

(2) 數(shù)值模擬的噴頭出口直徑和流量變化結(jié)果與林管[19]的試驗結(jié)果整體上相吻合，說明沖蝕耦合動網(wǎng)格方法用于計算噴頭沖蝕變形具有可行性，該方法對解決工程實際問題有較大的意義和應用價值。今后可利用該方法進一步對噴頭的抗沖蝕變形結(jié)構(gòu)優(yōu)化進行研究。