基于CFD灌溉機井管道液相流場的數(shù)值模擬

王 萍 ，楊 超

（1.天津工業(yè)大學(xué)電氣工程與自動化學(xué)院，天津 300387；2.天津工業(yè)大學(xué)電工電能新技術(shù)天津市重點實驗室，天津 300387）

隨著儀表技術(shù)的發(fā)展，超聲波流量計開始廣泛應(yīng)用于農(nóng)用灌溉領(lǐng)域，其優(yōu)點在于運行穩(wěn)定，計量精確度高，無壓力損失，并且能夠進行非接觸測量[1].當前，超聲波流量計研究領(lǐng)域的重點方向為傳感器探頭制作工藝、信號提取與處理的準確性以及低功耗.邱立存等[2]提出數(shù)字平均技術(shù)可使測量回波信號得到加強，提高測量精度；賀勝等[3]則對超聲波流量計的最優(yōu)聲道位置選取進行研究.但對于流體流場的變化與分布對流量計測量精度的影響研究較少.實際生產(chǎn)中，灌溉機井出水管道屬于非充分發(fā)展流場，這對時差法超聲波流量計的測量精度產(chǎn)生很大干擾.本文從理論上分析機井管道非充分發(fā)展管流情況下的流體速度分布特性，采用計算流體力學(xué)數(shù)值模擬與實流測試相結(jié)合的方法，確定三維流場的數(shù)學(xué)模型，探究流場分布規(guī)律，進而提出平均流速修正公式，通過實驗驗證，提高流量計的測量精確度.

1 時差法流量測量原理

時差法測量的原理是利用超聲波脈沖在流體中順流和逆流傳播過程中的時間差來計算流體流速，再將流速進一步換算成流量，如圖1所示.

圖1 超聲波流量檢測原理圖Fig.1 Ultrasonic flow detection schematics

2個超聲波換能器被斜向?qū)ΨQ安裝在被測量管道兩側(cè)，換能器 A、B 同時作為發(fā)生器（T1，T2）和接收器（R1，R2）交替工作，超聲波聲速為C，水流速度為u.設(shè)超聲波順流傳播時間為t1，逆流傳播時間為t2，流體流速的超聲波束傳播方向分量為ucos θ，W為截面寬度，在本例中即為圓管直徑D.由此導(dǎo)出t1與t2的表達式為：

進而推導(dǎo)出流體流速表達式為：

管道內(nèi)流量等于流體平均流速乘以管道截面積，計算表達式為：

由于W等于直徑D，所以流量進一步表示為：

在天津工業(yè)大學(xué)工程實訓(xùn)中心流量測量實驗室采用既有時差法流量計進行實際流速測量，通過調(diào)節(jié)變頻電機的轉(zhuǎn)速將水流初始速度設(shè)置為2.5 m/s，選取測量管道為90°彎管的機井出水管道模型.表1所示為出水管位置6D處測量速度的3次典型值.

表1 位置6D處水平軸截面的速度數(shù)據(jù)Tab.1 Velocity data of 6D horizontal axis section

由表1可以看出，在實際工況環(huán)境中，穩(wěn)定狀態(tài)下的流體流經(jīng)90°彎管后，在流量計傳感器安裝位置處測得流體流動速度變小.為提高時差法流量計的測量精度，本文建立管道內(nèi)三維流場的數(shù)學(xué)模型，對機井管道非充分發(fā)展管流情況下的流體速度分布特性進行分析.

2 流體力學(xué)求解

2.1 幾何模型的建立

本文選取某水泵公司某型號潛水泵為標準建立模型，其流量為70 m3/h，揚程22 m，輸水管徑為100 mm，經(jīng)計算流速約為2.5 m/s.利用Gambit軟件建立帶90°彎管的機井出水管道模型，模型取豎直方向長度為20倍管徑，水平方向長度為10倍管徑，彎徑曲率為1.5，所建模型如圖2所示.為了得到更好的收斂性與精確度，面網(wǎng)格采用錢幣畫法，在彎道處加密網(wǎng)格，網(wǎng)格總數(shù)量為9.63×105個.

圖2 網(wǎng)格劃分圖Fig.2 Diagram of meshing

2.2 設(shè)定邊界條件和流動區(qū)域

選擇入口類型為VELOCITY_INLET，出口設(shè)置為OUT_FLOW，其余管道壁面默認為固壁WALL.設(shè)定流動區(qū)域為機井出水管道的整個幾何實體，類型為“FLUID”.最后，將網(wǎng)格文件保存輸出.

2.3 利用Fluent求解器求解

將網(wǎng)格文件讀取到Fluent軟件中進行三維求解計算.由于流體為不可壓縮粘性流體，故求解器選擇非耦合隱式算法.通過對連續(xù)方程和N-S方程進行時均化處理，得到在笛卡爾坐標下的定常、絕熱、粘性、不可壓縮流體流動的控制方程[4].

連續(xù)性方程：

運動方程：

湍流數(shù)值模擬采用標準k-ε湍流模型來描述流體流動，該模型的湍動能k和耗散率ε的運輸方程為：

式中：Gk為由于平均速度梯度引起的湍動能產(chǎn)生項；Gb為由浮力引起的湍動能產(chǎn)生項；YM為可壓速湍流脈動膨脹對總的耗散率的影響.湍流黏性系數(shù)μt=ρCμk2/ε，模型常數(shù) C1ε=1.44，C2ε=1.92，Cμ=0.09，σk=1.0，σε=1.3[5-6].

本次試驗所選流體為液體，給定初始速度為2.5 m/s，溫度設(shè)定為20℃，密度為998.2 kg/m3，運動黏度為1.003 g/（m·s）.經(jīng)計算，管道雷諾數(shù)大于臨界值2 000，流體流動屬于湍流，由經(jīng)驗公式計算出湍流強度為0.033 8%，圓柱管道水力直徑即為管徑DH=0.1 m.考慮到潛水泵電機在持續(xù)做功克服重力，所以無需考慮重力作用.

2.4 殘差圖分析

初始化流場后，對流場模型進行迭代計算，圖3所示為模擬系統(tǒng)殘差圖.

圖3 系統(tǒng)殘差圖Fig.3 System residual plots

由圖3可以看出，迭代次數(shù)達到195次時模型基本收斂，說明在高雷諾數(shù)條件下網(wǎng)格的收斂速度比較快，計算所需的時間較少，基本達到設(shè)定條件.

2.5 圖形仿真結(jié)果

在模型求解過程中，離散方程組采用SIMPLE算法，選擇二階迎風格式進行求解.由于三維的模擬結(jié)果不便于觀察，所以創(chuàng)建內(nèi)部XY截面圖來查看計算結(jié)果，選擇三點創(chuàng)建等值面的方式來構(gòu)建觀察面.設(shè)定入口水流速度為2.5 m/s時，管道豎直軸截面的模擬云圖如圖4所示.

由圖4可知：

圖4 豎直軸截面的流體運動模擬云圖Fig.4 Fluid motion simulation nephogram of vertical cross-section

（1）在彎管豎直方向流體入口處，由于流體還未受到彎曲段的影響，圖 4（a）所示壓力始終平穩(wěn)，圖 4（d）流線排列整齊，流體速度均勻，由于邊界層摩擦力的存在，圖4（b）靠近管壁處顏色變暗流體速度降低.當流體繼續(xù)向上運動到彎管附近時，管道外側(cè)的流速開始顯著降低，并且越靠近外側(cè)下降幅度越大；在管道內(nèi)壁速度逐漸增加，在彎管最內(nèi)側(cè)最大速度達到4.08 m/s，是初始速度的1.63倍.觀察彎管附近的壓力云圖，彎管外側(cè)壓力明顯大于內(nèi)側(cè)，這是由于流體進入彎管的離心力所致.由于離心力的作用，在彎管橫向段出現(xiàn)二次流，二次流的擾動作用使流體速度提高[7].

（2）觀察彎管出口處的速度云圖和壓力云圖可知，流體速度出現(xiàn)分層現(xiàn)象并由內(nèi)側(cè)到外側(cè)逐漸升高，緊靠管壁內(nèi)側(cè)最低速度降至1.41 m/s，是初始速度的0.56倍.這種現(xiàn)象是由彎管的特殊結(jié)構(gòu)導(dǎo)致的，流動方向的驟變使管壁內(nèi)側(cè)流體動能損耗，速度降低，隨著相近流體的帶動速度逐漸增加.由圖4（b）、（c）、（d）可以看出，內(nèi)側(cè)流速始終小于外側(cè)，呈現(xiàn)出不同流速流體共流的現(xiàn)象[8-9].

由于實際生產(chǎn)中流量計的安裝位置位于彎管末端橫向管道上，所以取橫向管道上彎管90°截面以及2D、4D、6D、8D 4個位置的軸向截面作為觀察對象，為了便于觀察直接采用isometric視角截取上述5個截面的速度云圖進行對比，如圖5所示.

由圖5可知，隨著流體向前流動，二次流現(xiàn)象逐漸減弱，上層速度較高的區(qū)域面積逐漸減小，當流體運動到4D長度時中心顯示速度最接近初始速度2.5 m/s，但是速度分布仍然不均勻，云圖存在較為明顯的漩渦，只是強度減弱.對比后4個截面的速度矢量圖，如圖6所示.

圖5 典型位置截面速度云圖（m/s）Fig.5 Velocity nephogram of section at typical positions（m/s）

圖6 截面速度矢量圖（m/s）Fig.6 Velocity vector of different section（m/s）

由圖6可知，2種不同流速的區(qū)域在相互融合，上層速度逐漸降低，下層速度逐漸升高，并同時向中心區(qū)域的初始速度靠攏.

為了更直觀對比，分別取橫向出水管道豎直軸截面1/4、1/2、3/4處由上及下3條直線上的點繪制出坐標圖顯示速度的變化，如圖7所示.

圖7 豎直截面速度曲線Fig.7 Velocity curves of vertical cross-section

由圖7可以看出：在橫向管道中軸線上，流體在經(jīng)過彎管達到最大速度峰值后呈下降趨勢，在0.4 m附近處達到最小值并逐步上升至初始值2.5 m/s附近.在豎直軸截面上機井出水管道的不同層面以及不同位置的流速差別很大，流體速度分層明顯，數(shù)值紊亂，對實際的流量計測量造成誤差.因此，考慮在機井出水管道的水平軸截面安裝流量計，下面采用isometric視角簡要分析水平軸截面的流體運動狀況，如圖8所示.

圖8 水平軸截面的流體運動模擬云圖Fig.8 Fluid motion simulation nephogram of horizontal cross-section

由圖8可知，流體流經(jīng)彎管，在水平軸截面中心區(qū)域的流速低于兩側(cè)流速，在速度云圖上呈現(xiàn)低速漩渦；隨著流體繼續(xù)向前運動，中心區(qū)域流速逐漸升高，并與兩側(cè)融合接近.分別取水平軸截面由前向后1/4、1/2、3/4處的3條直線描繪速度變化，如圖9所示.

圖9 水平截面速度曲線Fig.9 Velocity curves of horizontal cross-section

由圖9可知，在水平軸截面管道兩側(cè)的流體速度變化趨勢基本重合.中間區(qū)域的流體速度與前文豎直截面分析結(jié)果相同，流體速度在流經(jīng)彎管獲得最大值后下降至最小值繼而逐漸上升與兩側(cè)速度融合，顯然這是受到二次流的影響[8-9].

2.6 數(shù)值模擬計算流速

通過上述分析，相對于豎直軸截面的流體速度分層且速度不對稱分布的情況，水平軸截面的速度分布更適合作為流量計探頭安裝位置進行流量檢測.選取入口速度為2.5 m/s時，在流速相對穩(wěn)定的6D位置處對水平軸截面進行3點數(shù)據(jù)統(tǒng)計.在該截面上取一條與管道夾角成60°的切線，將該切線四等分，測量每一段上的3點速度值.V1方向為水流流動方向，將速度矢量投影到水平和垂直2個方向上，如圖10所示.

圖10 水平軸截面速度矢量俯視圖Fig.10 Velocity vector planform of horizontal cross section

在Fluent軟件中輸入幾何模型6D位置處的切線上各數(shù)值求解點的X軸、Y軸和Z軸的坐標值，模擬出的流速數(shù)據(jù)如表2所示.

表2 水平軸截面的速度數(shù)據(jù)Tab.2 Velocity data of horizontal axis section

表2中的數(shù)據(jù)為經(jīng)Fluent軟件后處理后得出的數(shù)據(jù)，如圖10所示將聲波傳播路徑等分成4段，分別為 ab、bc、cd、de 段，其后再將每一段分別取其 1/4、1/2、3/4點位置計算出其速度值v以及水平方向速度值v1，這12個點的數(shù)據(jù)將作為數(shù)值模擬數(shù)據(jù)，與實驗室實際測量流速值進行比較，并完成流速數(shù)值的校正補償.

3 流量補償計算

在實際測量時，超聲波接收裝置安裝在機井出水管水平軸截面的前后兩側(cè)，通過對不同層面的流量補償系數(shù)的確定計算出最終的補償系數(shù)[10-15].

超聲波順流與逆流傳播時間不同，二者之差Δt為：

由于u cos θ<

表2統(tǒng)計了Fluent模擬條件下來流速度為2.5 m/s時，出水管道水平軸截面上位置6D處測量切線各點水平方向的速度數(shù)據(jù).為了提高測量精度，結(jié)合既有時差法流量計實驗室實流測試數(shù)據(jù)，分別對該切線4段的中點、1/4點以及3/4點處的速度值進行補償計算，進而求得最終管道流量補償值[14].4段中點處的時差之和以及實流測試時傳感器接收脈沖所需時差表達式為：

推導(dǎo)出中點處的補償系數(shù)與速度補償值：

同理，得到表1中每段1/4點處和3/4點處的補償系數(shù)：

由上述3組運算結(jié)果可計算得到平均補償系數(shù)與速度補償值：

經(jīng)過校正補償后管道內(nèi)流體流量為：

通過上述計算可知，當入水速度為2.5 m/s時，實驗室傳感器所測得的流體速度為2.35 m/s，平均流量為66.4 m3/h，約為CFD模擬所得速度的0.951倍，補償速度大約為0.12 m/s，即真實流速需要在傳感器實測值的基礎(chǔ)上增加0.12 m/s.補償后平均流量為69.8 m3/h，經(jīng)修正補償公式計算，測量精確度提升了4.9%.

4 結(jié) 語

本文通過建立農(nóng)用灌溉機井出水管道流體模型，完成了流場分布計算的理論推導(dǎo)，分析得出流體經(jīng)90°彎管產(chǎn)生的二次流使出水管流場分布急劇紊亂，經(jīng)對比，選擇具有分布對稱性的水平橫截面作為超聲波探頭安裝位置.為降低非均勻發(fā)展的彎管二次流流場對時差法超聲波流量計測量精確度的影響，基于本文提出的某型號水泵，結(jié)合CFD仿真數(shù)據(jù)與實驗室實際測量數(shù)據(jù)，提出平均流速修正公式，求解具體情況下的流速與流量補償系數(shù).結(jié)果表明，經(jīng)補償后流量計測量精確度提升4.9%.本研究在一定程度上降低了流場分布不均對流量計測量精確度的影響，對現(xiàn)有各種應(yīng)用在非充分發(fā)展流場的時差法超聲波流量計的程序編程及流量修正具有參考價值.

[1]陳震，范永申，周慶峰，等.常用農(nóng)業(yè)灌溉水表精確性及耐砂耐凍測試評價[J].灌溉排水學(xué)報，2015，34（10）：19-24.CHEN Z，F(xiàn)AN Y S，ZHOU Q F，et al.Evaluating accuracy and duration sediment of wear test and frozen treatment of agriculture common watermeters[J].Journal of Irrigation and Drainage，2015，34（10）：19-24（in Chinese）.

[2]邱立存，王汝琳.超聲波氣體流量測量系統(tǒng)的實現(xiàn)[J].傳感器與微系統(tǒng)，2006，25（1）：47-49.QIU L C，WANG R L.Realization of ultrasonic gas flow measurement system[J].Transducer and Microsystem Technologies，2006，25（1）：47-49（in Chinese).

[3]賀勝，彭黎輝，仲里敏.基于CFD的超聲波流量計最優(yōu)聲道位置研究[J].儀器儀表學(xué)報，2009，30（4）：852-856.HE S，PENG L H，ZHONG L M.Computational fluid dynamics based sound path optimization for ultrasonic flow meter[J].Chinese Journal of Scientific Instrument，2009，30（4）：852-856（in Chinese）.

[4]NORLI P，LUNDE P，VESTRHEIM M.Investigation of precision sound velocity measurement methods as reference for ultrasonic gas flow meters[J].Rotterdam，Netherlands:IEEE Ultrasonics Symposium，2005，3（2）：1443-1447.

[5]王萍，萬凱.大口徑矩形彎管氣體流場建模與仿真[J].天津工業(yè)大學(xué)學(xué)報，2013，32（3）：82-84.WANG P，WAN K.Modeling and simulation of rectangular pipe with large diameter gas flow[J].Journal of Tianjin Polytechnic University，2013，32（3）：82-84（in Chinese).

[6]王福軍.計算流體動力學(xué)分析[M].北京:清華大學(xué)出版社，2004：121-125.WANG F J.Computational Fluid Dynamics Analysis[M].Beijing：TsinghuaUniversityPress，2004：121-125（inChinese）.

[7]饒永超，常凱，王樹立，等.90°彎管內(nèi)流體流動特點數(shù)值模擬[J].常州大學(xué)學(xué)報：自然科學(xué)版，2016，28（3）：65-69.RAOYC，CHANK，WANGSL，etal.Numericalsimulationof fluidflowcharacteristicin90°bendpipe[J].JournalofChangzhou University：Natural Science Edition，2016，28（3）：65-69（in Chinese).

[8]ZHAO Z H，F(xiàn)ERNANDO H J S.Numerical simulation of scour around pipelines using an Euler-Euler coupled two-phase model[J].Environmental Fluid Mechanics，2007，7：121-142.

[9]江山，張京偉，吳崇健，等.基于FLUNET的90°圓形彎管內(nèi)部流場分析[J].中國艦船研究，2008，3（1）：37-41.JIANG S，ZHANG J W，WU C J，et al.Numerical simulation of inner flow in 90°bending duct of circular section based on fluent[J].Chinese Journal of Ship Research，2008，3（1）：37-41（in Chinese).

[10]吳春華，鮑敏.超聲波流量計的彎管誤差分析及修正研究[J].機電工程，2015，32（2）：175-179.WU C H，BAO M.Bending pipe error analysis and correction of ultrasonic flowmeter[J].Journal of Mechanical&Electrical Engineering，2015，32（2）：175-179（in Chinese).

[11]BRASSIER P，HOSTEN B，VULOVIC F.High-frequency transducers and correlation method to enhance ultrasonic gas flow metering[J].Flow Measurement&Instrumentation，2001，12（1）：201-211.

[12]方翔.超聲波氣體流量計的研制[D].武漢：華中科技大學(xué)，2004.FANG X.A thesis submitted in partial fulfillment of the requirements for the degreeofmaster of engineering[D].Wuhan：Huazhong University of Science and Technology，2004（in Chinese）.

[13]TABIB M V，ROY S A，JOSHI J B.CFD simulation of bubble column：An analysis of interphase forces and turbulence models[J].Chemical Engineering Journal，2008，139（3）：589-614.

[14]榮雙.大口徑煙道氣體流量測量方法的研究[D].天津:天津工業(yè)大學(xué)，2016.RONG S.Research on large diameter flue gas flow measurement method[D].Tianjin：Tianjin Polytechnic University，2016（in Chinese）.

[15]MERIBOUT M，RAWAHI N Z，NAAMANY A M，et al.A multisensor intelligent device for real-time multiphase flow metering in oil fields[J].IEEE Transactions on Instrumentation&Measurement，2010，59（6）：1507-1519.