液力透平內氣液兩相流動的數值模擬

畢智高 劉甜甜 謝 勰 相玉琳 王金璽

（1.榆林學院化學與化工學院；2.陜西省低變質煤潔凈利用重點實驗室；3.長慶油田分公司規劃計劃處）

泵作液力透平（PAT）回收液體富余壓能因具有結構簡單、安裝方便、占地面積小、可批量生產及運行維護費用低等優點［1，2］，而廣泛應用于油氣儲運、石油化工及煤化工等行業。 這些工業流程中作為能量回收的工質往往呈氣、 液兩相狀態，入口含氣率（IGVF）對液力透平的性能有著直接的影響［3］。

近年來，關于旋轉可逆式水力機械中氣液兩相流動的研究主要集中于泵［4～12］工況，對PAT工況的研究多以純液為工質，而氣液兩相條件下的研究報道相對較少。 文獻［13～16］采用CFD軟件分別分析了IGVF對PAT的外特性、內流規律和所受徑向力的影響，文獻［17］通過瞬態數值模擬，獲得了氣體在PAT內的分布規律和流場的壓力脈動變化規律，但是還不夠深入。 筆者在借鑒泵內氣液兩相流研究成果的基礎上， 考慮入流流型，基于雙流體模型， 以空氣和水作能量回收工質，對一固定徑向導葉式PAT進行定常數值計算， 旨在為氣液兩相條件下PAT的流動特性及其規律提供部分理論參考。

1 數值方法

1.1 Euler-Euler雙流體模型

Euler-Euler多相流模型可分為均相流和非均相流模型。 均相流模型是基于均質平衡流理論，將氣液兩相流視為均勻混合物，不計相間速度滑移，并且忽略了客觀存在的兩相界面、相間相互作用等許多重要因素［18］，導致求解誤差較大；非均相流模型將氣液兩相之間的界面視作移動的邊界，考慮相間速度滑移、質量與動量傳遞等，每相流體都有各自的流場并通過相間傳遞單元進行傳遞， 即每相都有各自的速度場和溫度場，但共用壓力場，通過相間作用力和熱量傳遞使得兩相速度和溫度得到平衡，因而較符合實際。

筆者采用非均相流模型，液相為不可壓連續相，氣相為不可壓分散相，相傳遞單元采用Particle模型，動量傳遞方式采用Schiller Nauman模型，假定氣液兩相入流為泡狀流型條件。

1.2 控制方程

連續性方程：

動量方程：

式中 fk——與葉輪旋轉有關的質量力；

k——代表液相（l）或者氣相（g）；

Mk——相間作用力；

pk——k相壓力；

wk——k相相對（滑移）速度；

ρk——k相密度；

βk——k相體積分數；

μk——k相動力粘度。

氣、液兩相和兩相混合工質滿足：

式中 Qg、Ql、Q——氣相、 液相和兩相混合工質的體積流量，m3/s；

βl——體積含液率，%；

βg——體積含氣率（氣相濃度），%。

1.3 相間作用力

式中 CD——無量綱曳力系數；

dB——氣泡直徑；

ReB——氣泡雷諾數；

wR——滑移速度；

νl——液相運動粘度。

虛擬質量力Mk，VM是由兩相間的相對加速引起的，按文獻［20］中的公式計算如下：

式中 aVM——虛擬質量加速度；

CVM——虛擬質量系數， 通常為氣相含量的函數，球形氣泡取0.5。

2 計算模型與邊界條件

2.1 計算模型及網格劃分

以一徑向導葉式PAT為研究對象， 工況下的性能參數和主要結構參數依據文獻［21］。 采用對復雜幾何邊界適應性強的非結構四面體網格對計算域進行網格劃分，并對網格進行了無關性驗證， 當水力效率波動小于0.45%時即滿足網格無關性假設。 考慮到計算機的運算能力、計算精度和計算要求，最終確定網格數量在120萬，且網格質量均在0.3以上。

2.2 假設條件

氣相為粒徑均勻的球形，其直徑遠小于流道特征尺寸；不考慮兩相重力；忽略系統內部、系統與外界間的質/熱量傳輸， 系統內無化學反應發生，溫度恒定。

2.3 邊界條件及數值求解

液相采用能更好處理高應變率和流線彎曲程度較大流動的RNG k-ε湍流模型， 中等湍流強度，氣相采用湍流零方程模型；光滑壁面，固壁對液相為無滑移作用、氣相為自由滑移，近壁區的壁面函數采用scalable； 熱量傳遞采用等溫模型，設定溫度為25℃； 設定質量流量進口，IGVF分別為5%、15%、25%、40%，壓力出口；葉輪流體域設置為旋轉，其余計算域設置為靜止，導葉-葉輪及葉輪-出水管間的動靜交界面設置為Frozen Rotor模式。 控制方程的離散采用基于有限元的有限體積法，對流項和湍動能相均采用二階迎風格式求解， 計算迭代步數設置為2 000， 求解殘差類型RMS設為10-5。

3 結果分析

3.1 試驗驗證

圖1為PAT試驗臺，通過試驗測試來驗證數值計算的準確性。

圖1 PAT試驗臺示意圖

用一臺比轉數ns=84.5的單級單吸離心泵作PAT在IGVF為20%時， 性能參數數值計算與試驗測量的結果對比如圖2所示。 由圖2可見，兩者雖有誤差（誤差的主要因素是數值計算過程中忽略了前、后腔內部的流體），但在最高效率時的數值計算與試驗測量的結果誤差小于5%，可見筆者所采用的數值計算方法可以較好地對PAT性能進行預測［17］。

圖2 泵性能參數數值計算與試驗測量的結果對比

3.2 外特性曲線

氣液兩相條件下，PAT的外特性參數透平壓頭H（m）、透平效率η（%）和混合工質流動密度ρf（kg/m3）的計算公式分別為：

式中 g——重力加速度，m/s2；

G——混合工質的質量流量，kg/s；

M——葉輪扭矩，N·m；

pin、pout——進、出口總壓，Pa；

ρg、ρl——氣、液相密度；

ω——葉輪旋轉角速度，s-1。

圖3為不同IGVF下PAT的外特性曲線，表1列出了PAT最高效率時的性能參數。

結合圖3、表1可知，不同IGVF下PAT的流量-效率、流量-壓頭和流量-功率曲線的總體趨勢基本一致，即隨著兩相混合工質流量的增加，PAT的效率先較快增加后緩慢減小，壓頭和功率隨著流量的增加而增加。 隨著IGVF的增加，相同流量下PAT的效率和對外輸出的功率降低，壓頭升高，小流量下的效率降低和壓頭升高較大流量下的明顯，而ρl遠大于ρg導致了大流量下的功率降低較小流量下的明顯。PAT最高效率時的流量發生偏移，與IGVF為5%時相比，最高效率分別下降了1.2%、2.2%、4.5%。

圖3 不同IGVF下PAT外特性曲線

表1 不同IGVF下PAT最高效率時的性能參數

3.3 蝸殼流道氣相分布

圖4為蝸殼流道各截面位置，圖5為不同IGVF下PAT最高效率時蝸殼流道各截面氣相分布情況。 由圖5可見，同一IGVF下，工質從進口沿蝸殼流道（從截面1～8）氣相分布和氣相濃度均勻程度均降低。

圖4 蝸殼流道各截面位置

圖5 蝸殼流道各截面氣相分布

截面1靠近PAT進口端，氣相均勻分布，數值接近初始IGVF， 與泡狀流型的入流假定條件相符，經過截面1后，氣相向流道半徑較小的蝸殼截面中下方匯聚堆積，形成高含氣區，在此可能發生流型轉換和相態分離，這是因為兩相混合工質沿蝸殼流道在做減壓增速運動，而液相的密度較大，所受到的慣性離心力大于氣相，導致大量液相裹挾少量氣泡繞蝸殼外緣壁面運動，而在液相的排擠作用下，大量氣泡被迫向蝸殼內緣出口處轉移，從而形成高含氣區。

另外，截面8處氣相濃度最低，且氣相均布于整個截面；同一截面處的氣相濃度隨IGVF的升高而增大。

3.4 導葉葉輪內氣相分布

圖6給出了不同IGVF下PAT最高效率時導葉葉輪內截面上的氣相分布情況。 由圖6可見，PAT各過流部件中蝸殼內氣相分布相對最均勻，各過流部件的氣相濃度均隨著IGVF的升高而增大，氣相分布的不均勻度增加，流動紊亂度增強，這也是隨IGVF的增加PAT效率下降的原因。 匯聚堆積在蝸殼截面下方的氣相因慣性隨液流順勢就近進入附近的導葉和葉輪流道，因此蝸殼流道內高含氣率截面對應中心角所包含范圍及其附近區域的氣相濃度也相對較高，高低含氣區域呈現出明顯的分界（圖6b）。 高含氣率葉輪流道內氣泡聚向葉片吸力面，且隨著IGVF的增加，氣相聚集程度增強，氣泡聚并成為大的氣團而滯留于葉輪流道，改變了有效過流面積，甚至可能造成“氣堵”。 這是由于在有限葉片數的葉輪旋轉效應和葉形曲率的影響下，液相在進入葉輪后受到較大的慣性離心力和科氏力作用改變了運動狀態而傾向葉片壓力面運動，而氣相受到的慣性離心力和科氏力則相對較小，但在葉片壓力面和吸力面之間壓力梯度的驅動和液相的排擠下，導致了葉輪流道內的氣相向葉片吸力面聚集。 對比圖6c中的液相速度流線，發現在滯留氣團的上游，存在對應漩渦，表明葉輪流道內漩渦形成與氣相的聚并相關。 因此，在對含氣工況下應用的固定徑向導葉式PAT進行優化設計時， 可考慮通過改變導葉開度、 導葉葉輪的葉片數及其出口面積等方案，控制流動狀態和調節壓力梯度以削弱氣相在葉輪流道內聚并與滯留的程度。

圖6 導葉葉輪內氣相分布

3.5 葉片表面氣相分布

圖7給出了不同IGVF下PAT最高效率時葉片表面氣相分布情況。 由圖7可見，不同IGVF下，葉片表面氣相分布并不均勻，IGVF較低時壓力面氣相主要分布在葉片高壓邊側，而葉片低壓邊側有少量氣相聚集。 隨著IGVF的增加，氣相向葉片表面擴散，整個壓力面上葉片高壓邊附近的氣相濃度明顯高于其他區域，且前蓋板側氣相濃度高于后蓋板側。 整個吸力面上葉片低壓邊附近的氣相濃度明顯高于其他區域，這也與氣團在葉輪流道內滯留的位置相符。

圖7 葉片表面氣相分布

對比圖8的氣相速度矢量分布可見， 兩相混合工質進入葉輪流道后，氣相在葉片壓力面的速度方向指向前蓋板側，表明氣相速度場的分布與氣相聚集位置相對應，而葉片吸力面中部的氣相則沿葉片分別向高壓邊和低壓邊的兩側運動，這是因為葉片吸力面中部區域流體的總壓 （機械能）高于兩側高壓邊和低壓邊。 葉片吸力面逆流而上的氣相與從導葉進入的混合工質在葉輪進口區域碰撞摻混，同時還受葉輪、導葉和動/靜葉柵的干涉擾動，導致該區域流動紊亂，產生漩渦；順流而下的氣相匯入氣團而滯留于葉輪出口區域附近。

圖8 葉片表面氣相速度矢量分布（IGVF=25%）

3.6 氣液兩相滑移速度分布

圖9為IGVF=25%時PAT中截面氣液兩相滑移速度分布情況。 由圖9可見，在導葉、葉片頭部及尾部等局部滑移速度為負值，即氣相速度小于液相速度。 總體而言，PAT內氣相速度大于液相速度，這種差異在可能發生“氣堵”的葉輪流道進口區域中更為明顯。 這是由于旋轉葉輪與固定導葉之間劇烈的干涉擾動， 使該區域的流動紊亂，相間作用增強，由式（6）可知，氣液兩相相互作用越強，兩相間的速度差越大，從而導致了氣液兩相滑移速度的增大。 與圖6c對比分析后可知，漩渦區上游的兩相滑移速度大于下游的，而漩渦下游區域也是氣相聚集處，這表明氣相聚結為氣團而發生滯留， 降低了該處氣液兩相的滑移速度，導致了相態分離，同時也說明漩渦上游區域的流動較為紊亂。

圖9 PAT滑移速度分布（IGVF=25%）

3.7 不同流量下PAT氣相分布

圖10為不同流量下IGVF=25%時PAT中截面上的氣相濃度與氣相流線分布情況。 由圖10可見， 不同流量下PAT流道內的氣相分布規律較為相似，蝸殼內氣相分布相對均勻，高、低含氣區域分界明顯，部分葉輪流道靠近葉片吸力面存在氣相聚結區，氣相聚結區存在相態分離。 隨著流量的增加， 葉輪流道內氣相聚結程度和范圍增強（2.2Qd共1個流道，3.0Qd共3個流道，3.4Qd共4個流道）。 這是因為流量越大， 混合工質的動能也越大， 同時氣液兩相各自的體積含量成比例增加，相間的相互作用增強。氣相流線在2.2Qd時較為紊亂，存在較為明顯的漩渦和脫流，加劇了能量耗散，隨著流量的增大，這種情況有所改觀。

圖10 不同工況PAT氣相濃度及氣相流線分布（IGVF=25%）

4 結論

4.1 隨著IGVF的增加，相同流量下PAT的水力效率和輸出功率降低而壓頭升高，最高效率時對應的流量發生偏移，與IGVF為5%時相比，最高效率分別下降了1.2%、2.2%和4.5%。

4.2 相同IGVF下， 沿蝸殼流道截面氣相濃度和氣相分布均勻度均降低，氣相聚集于蝸殼截面中下方；同一蝸殼截面處的氣相濃度隨IGVF的升高而增大。

4.3 PAT各過流部件中，蝸殼內氣相分布相對最均勻，隨IGVF和流量的升高，各過流部件的氣相濃度和分布的不均勻度均增加， 存在明顯的分界。

4.4 葉輪流道內氣相聚向葉片吸力面， 隨著IGVF和流量的增加，聚集程度和范圍增強，氣相聚并成為氣團而滯留于葉輪流道， 發生相態分離，可能造成“氣堵”；氣團附近上游存在對應漩渦， 表明葉輪流道內漩渦的形成與氣相聚并相關；IGVF較低時葉片壓力面氣相主要聚集在高壓邊側，隨著IGVF的增加，氣相在葉片表面擴散，葉片壓力面高壓邊附近靠近前蓋板側的氣相濃度較高，葉片吸力面低壓邊附近的氣相濃度較高。

4.5 PAT內的氣相速度總體上大于液相速度，且在氣團滯留的葉輪流道進口區域較明顯。