袁顯寶, 劉 超, 譚 偉, 張永紅*, 張彬航, 程子陽
(1.三峽大學機械與動力學院, 宜昌 443002; 2.湖北省水電機械設備設計與維護重點實驗室,宜昌 443002)
隨著溫室效應、氣候變化等環保問題的日益突出,碳排放的有效控制是當今社會在發展中需要考慮的重要問題。利用碳捕獲、利用與封存(carbon capture,utilization and storage,CCUS)技術將CO2資源化的技術受到各國重視[1-5]。成熟的CCUS技術可將含有雜質的CO2氣體經過分離、凈化和壓縮處理為超臨界CO2,并通過管道運輸。當輸送流質中混雜有游離水時,在溫壓產生波動條件下易形成CO2固態水合物,其物理性質與冰類似。在輸送過程中, CO2水合物會與管壁發生碰撞將對管壁造成沖蝕破壞,影響管道傳輸安全性;沖蝕坑為水合物提供生長附著點,水合物沿管壁堆積將影響輸送效率,甚至可能堵塞整個管道[6-7]。彎管是維系管道輸送系統結構完整的必要結構,而據大量研究表明,含有固體顆粒的兩相流對于彎管的侵蝕危害程度相對于直管段更加嚴重[8-11]。
目前,根據應用方向的差異,對管道壁面沖蝕研究主要從混雜有固體顆粒物的水、油、氣三種材料進行開展,且相關研究成果較多。但對于超臨界CO2輸送管道內水合物對管道沖蝕的研究較少。王博等[12]研究了氣固混合流集輸管道的彎管段受顆粒沖蝕影響的因素;劉丹[13]利用實驗方法研究了超臨界CO2流動狀態對碳鋼表面的腐蝕效應的影響;李樂等[14]利用實驗方法研究了在流量變化對高壓環境下CO2水合物生成過程的影響;巨熔冰[15]利用FLUENT軟件研究了超臨界CO2流體產生的應力和湍流動能對管壁的影響;劉洪亮等[16]使用FLUENT軟件求解了壓裂工況下T形三通沖蝕影響因素。超臨界CO2輸送管道中水合物對管道的沖蝕效應未見更詳細的研究報道。因此,根據超臨界CO2輸送管線的輸送工藝和固液兩相物質的屬性,采用數值模擬方法,選取輸送管路中的彎管段的沖蝕情況進行計算,對于輸送管道的防護和安全性評估有重要參考價值,有利于CCUS技術的發展。現基于COMSOL多物理場耦合軟件研究超臨界CO2輸送彎管的沖蝕效應,求解不同流速、顆粒直徑和不同彎角對彎管沖蝕效應的影響。
使用COMSOL Multiphysics軟件的流體流動顆粒追蹤模塊與湍流模塊耦合求解流場中微粒的運動軌跡。通過對每個分量建立相應的常微分方程,計算每個時間步內,粒子受連續相的作用力與粒子間相互作用力的合力,求解粒子運動狀態。離散相顆粒的受力與慣性力保持平衡[17],流體對粒子的曳力FD占主導作用,粒子的受力方程可表示為
(1)
(2)

(3)

(4)
式中:Vp為離散相粒子速度,m/s;dp為粒子直徑,m;ρp為粒子密度,kg/m3;Rep為粒子雷諾數;CD為阻力系數;FD為流體拖曳力;Vf、ρf分別為流體速度、密度大小,m/s1、kg/m3;μf為流體的動力黏度,Pa·s;a1、a2、a3為經驗常數,與粒子雷諾數和顆粒形狀相關。
除受到曳力外,粒子在相載流體中受到浮力FB、壓力梯度FP、附加質量力FA、旋轉力FR,各力的表達式為

(5)

(6)

(7)
(8)
式中:?p為粒子容積范圍內壓強梯度,Pa;ωp為粒子的旋轉角速度,rad/s。
引入Stokes數St來描述粒子運動對流場變化的響應特點,其定義為粒子相應時間與流體傳播時間的比率[18]。其表達式為
(9)
式(9)中:D為外部容器的管徑,m。
參考Menter[19]對于k-ω湍流模型應用于強流線曲率流動準確性的研究工作,采取標準k-ω湍流模型來求解彎管沖蝕過程中的湍流流場。其控制方程為
(10)
(11)
式中:k為湍動動能;ω為湍流動能比耗散率;u為流體微團的速度;i、j分別代表X、Y方向;Γk、Γω為和k、ω的擴散率;Gk為由速度梯度引起的湍流動能生成項;Gω為由比耗散率的生成項;Yk、Yω分別為湍動能、比耗散率的發散項;Sk和Sω分別為湍動能用戶自定義量、比耗散率用戶定義量。
顆粒隨流體運動過程中,與管壁的碰撞伴隨著由摩擦和非彈性效應產生的能量轉移和損失,主要表現在入射粒子和反射粒子在速度分量的改變,通常以碰撞前后速度分量之比定義恢復系數?,F采取Forder等[20]提出的粒子與壁的碰撞回彈模型來描述恢復系數與沖擊角度之間的關系。
(12)
式(12)中:en為法向恢復系數;et為切向恢復系數;vn為法向粒子速度分量;vt為切向粒子速度分量;θ為粒子碰撞角度(入射速度與表面切線之間的角度);下標1、2分別為碰撞前、后的狀態。
沖蝕效應的嚴重程度主要通過沖蝕磨損率表示,其定義為運動的顆粒物的沖擊造成的靶材上單位時間單位面積上的質量損失。選取COMSOL中的E/CRC沖蝕模型來模擬彎管沖蝕過程中的沖蝕效應。該模型由E/CRC沖蝕與腐蝕聯合研究中心(Erosion/ Corrosion Research Center)和Tulsa大學提出[21],沖蝕磨損率EM的計算公式為
(13)
式(13)中:K為侵蝕模型系數,缺省值為2.17×10-9;BH為表面材料的布氏硬度,缺省值為200;V為入射粒子速度,m/s;Vref為入射粒子的參考速度,缺省值為1,m/s;FS為顆粒形狀系數,缺省值為0.2,代表圓形顆粒;n為模型指數,缺省值為2.41;α為粒子入射角,rad。
基于參考文獻[15]所提出的超臨界CO2輸送管道的設計方案和運行參數來建立基準計算模型。建立如圖1所示的幾何模型,其結構由短管段-彎管段-長管段三部分組成,流體由短管段進入,長管段流出。管道直徑D=200 mm,管壁厚度δ=10 mm,曲率半徑R=200 mm,短管段的長度L1=400 mm,為了使彎管后段的湍流充分發展,出口直管段長度設置直管段長度L2=800 mm。選取研究區域彎管體積的1/2作為建模區域,減少所需的計算內存資源和求解時間。

圖1 彎管模型的幾何結構
當超臨界CO2輸送管道外界溫度在10~25 ℃,水體積含量大于21 000 ppm并且壓力大于約8 MPa時,CO2水合物易生成且對彎管沖蝕較為嚴重[15]。管內兩相流中液相為超臨界CO2,在40 ℃時密度為1 100 kg/m3,動力黏度為3×10-5Pa·s;離散相為CO2水合物顆粒,密度為900 kg/m3,粒徑分布在5×10-6~200×10-6m。管道的管壁材料選取實際應用于CO2超臨界與密相傳輸的管材X80管線鋼,其材料屬性數據如表1所示。

表1 材料物性參數
參考李井洋[22]對大口徑輸送管道中粒子沖蝕因素的研究方法,使用COMSOL計算粒子沖蝕問題分兩步進行,使用湍流模塊求解連續相的流場分布,在此基礎上使用流體流動顆粒跟蹤模塊求解顆粒運動和沖蝕效應。
湍流模塊中,以圖1中的面5、面7為速度入口邊界,設置為充分發展的流動,平均流速為2 m/s;面42為壓力10 MPa的出口邊界;施加方向為Z軸負方向的重力場;設置YZ剖面為對稱邊界。
流體流動粒子跟蹤模塊中,為粒子指定密度和直徑,其密度為900 kg/m3,平均粒徑為80 μm,粒子體積濃度設為2.1%;流體與管壁接觸面設定壁反彈邊界條件,同時增設E/CRC沖蝕邊界;設定入口面為面7,速度場繼承自湍流流場計算的結果;出口邊界為面42,壁條件設為消失;流域內的粒子受到斯托克斯曳力作用,湍流彌散模型設為離散隨機游走;設置YZ剖面為對稱邊界。
求解器分步進行研究,使用穩態研究求解湍流流場及對管壁的應力,使用瞬態求解計算2 s內粒子跟隨流體流動狀態,時間步長為10-5s,容差均設為10-3以保證精度;采用內置的MUMPS求解器進行計算。
采取非結構化單元進行網格繪制,網格模型整體概況如圖2所示。邊界層處細化為5層以求解邊界處的湍流流動狀況,生長因子為1.2。以求解得到流場的最大流速為參考標準,設置多組密度的網格進行計算進行網格無關性驗證。計算結果如圖3所示,網格單元數在到達18×104時求解結果基本穩定,最終的網格劃分方案網格單元數為178 587,最小質量0.148,平均質量0.713。

圖2 彎管模型的網格模型

圖3 網格數量與對應的最大流速
平均入口流動速度2 m/s,粒徑為80 μm,彎曲角度為90°的彎管,中心對稱面上的超臨界CO2流體速度和壓力分布分別如圖4(a)、圖4(b)所示。
由圖4(a)所示的速度分布云圖可知,上游短管段的湍流流場充分發展,管道中心區域的流速高于近壁區的流速,平均流速約為2 m/s;進入彎管區域后,彎管內側流體獲得較大流速,約為3 m/s,外側區域流速約為1.5 m/s;流體由彎管段向直管段過渡時,發生了一定程度的流動分離,靠近內側管壁附近的流動速度較小而靠近外側區域流速較大。在下游長管段中,高流速區域貼近外側管壁且截面逐漸縮小,管內整體流速降低到2 m/s。根據流場的分布可知,當水合物粒子隨流體流動時,受曳力和粒子自身向心力作用,在長管段外側的壁面容易遭受一次粒子的沖蝕破壞。
由圖4(b)所示的壓力分布云圖可知,彎管附近存在壓力變化,彎管內側由于流體的流速較快,壓力低于外側區域,壓差約為4 000 Pa;超臨界CO2流經彎整個管段過后產生了約8 000 Pa的壓降。

圖4 對稱面上超臨界CO2流場分布
圖5為直角彎管在前面所設條件下壁面受粒子沖蝕影響云圖,由圖5中可知,受到CO2水合物沖蝕較為嚴重的區域集中在彎管段外側管壁與長管段交接處,對應圖4中沿直管壁的高流速區域的起點位置,管道壁面所受到的最大沖蝕速率為1.37×10-9kg/m2/s。

圖5 基準條件下彎管表面沖蝕云圖
為研究入口平均流速大小對彎管沖蝕磨損速度的影響,設置入口平均流速從1 m/s變化到5 m/s,速度每增加0.5 m/s進行一次計算,不改變其他條件。彎管壁面的最大沖蝕率與入口平均流速關系如圖6所示。由圖6可知,對于直角彎管模型,當入口平均流速增大時,彎管最大沖蝕率隨之增大,基本呈現出二次函數增長的關系。

圖6 最大沖蝕速率隨入口速度變化曲線
另外,平均入口速度發生改變時,管壁受到沖蝕影響區域的位置也將改變。圖7為不同的入口平均流速條件下彎管受到沖蝕影響區域的對比。從圖7中可知,流速增大時,沖蝕效應增強,沖蝕中心更加集中并向彎管拐角處偏移。造成這種現象的原因是,水合物粒子在彎管中隨流體運動時,除了受流體曳力影響外還受到自身慣性力作用。當流體的流動速度增加時粒子速度隨之增加,此時慣性力對粒子運動軌跡的影響增強,流體流動方向突然改變時,粒子容易脫離流體運動方向,與管壁發生碰撞,導致沖蝕效應增大且沖蝕中心位置向下方移動。

圖7 沖蝕區域受入口速度大小變化的影響
為研究顆粒直徑大小對彎管沖蝕磨損速度的影響,設置CO2水合物的顆粒直徑從20 μm變化到250 μm。顆粒的直徑與彎管壁受到的最大沖蝕率之間關系曲線如圖8所示。由圖8可知,隨著粒子粒徑增加,最大沖蝕率隨之增加;當微粒的直徑在20~80 μm變化時,最大沖蝕速率受粒徑變化的影響較小,基本呈二次函數關系,而顆粒直徑在80~250 μm變化時,粒子直徑的增長對沖蝕率的影響較小,呈線性增長關系。圖9為不同粒子粒徑對管壁造成沖蝕區域的對比。由圖9可知,粒徑逐漸增大時,沖蝕區域的位置和形狀發生改變,由直管段外側壁面的彌散分布轉變為彎管外壁的集中橢圓形分布。這是由于,由于粒徑較小時,Stokes數值小,粒子跟隨性強,與管壁發生接觸時位置通常在直管段,沖擊角度與管壁夾角較小,對管壁沖蝕主要以切削作用為主。粒子與管壁碰撞時的動能與粒子速度和粒子質量的增加均相關,因此沖蝕率的增長曲線接近二次函數,但數值變化較??;當粒子直徑增大到一定程度時,Stokes數值較大,粒子跟隨性較弱,在彎管段受到向心力加速的程度弱,粒子容易脫離流體與管壁發生碰撞,彎管段的管壁發生正碰,因此沖蝕程度較為嚴重。此時粒子速度接近初始速度,粒子與管壁碰撞時的動能只受到粒子質量變化的影響,沖蝕增長率與粒徑呈一次函數關系。

圖8 最大沖蝕速率隨粒徑變化曲線

圖9 不同粒徑條件下的沖蝕區域分布
為研究彎管角度(長短直管的夾角)對彎管受到沖蝕磨損速度的影響,設置彎管角度分別為30°、45°、60°、90°、120°、135°、150°。管壁受到的最大沖蝕率與彎管角度的關系如圖10所示。由圖10可知,管壁受最大沖蝕率在彎管角度小于90°時呈上升趨勢,當彎管角度大于90°時呈下降趨勢;角度為90°時,彎管受到的沖蝕效應最大。角度在90°前,最大沖蝕率將隨角度增大而平緩增加,角度超過90°以后,沖蝕效應迅速減弱。

圖10 不同彎管角度下的最大沖蝕率變化曲線
圖11為不同角度彎管受到粒子沖蝕區域位置對比。從圖11中可知,當彎管角度小于90°時,粒子對管壁造成的一次沖擊區域較集中;當彎管角度大于90°時,沖蝕影響區域較分散,且隨著夾角增大,沖蝕區域分散程度增加。彎管角度較小時,管內流體流動方向發生劇烈偏轉,粒子受到流體曳力和向心力不足以抵消粒子自身的慣性力,粒子與彎管的沖蝕區域與直角彎管的沖蝕區域幾乎相同;而流體在夾角大于90°的彎管流道中流動時,方向偏轉小,粒子在與管壁碰撞前將隨流體流動得更遠;與管壁發生碰撞時,主要以切削為主,無法對管壁造成較大損害。

圖11 不同彎管角度條件下沖蝕區域分布
利用COMSOL湍流流動模塊和顆粒跟蹤模塊,耦合求解超臨界CO2輸送管道受到CO2固體水合物沖蝕效應,得出以下結論。
(1)隨著管道入口平均流速的增加,水合物顆粒的速度隨之增加,CO2固體水合物對彎管的沖蝕程度將增大,沖蝕區域由彎管段外側較遠處的彌散分布轉變為彎管外側拐角處集中分布。采取較小的輸送速度能夠有效降低水合物對管壁的沖蝕。
(2)隨著CO2水合物顆粒的粒徑增大,管壁上受到的最大沖蝕率不斷增大,粒徑較小時,沖蝕不明顯;粒徑超過約80 μm時,沖蝕速率的增長與粒徑呈一次線性關系。顆粒粒徑增加將使沖蝕區域由遠離拐角的直管外側彌散分布向拐角靠攏聚集。
(3)固體水合物顆粒隨超臨界CO2流體流經不同角度的彎管時,對管壁的沖蝕程度不同,直角彎管更容易受到粒子沖蝕。當彎管的角度大于90°時,管壁受粒子沖蝕的區域相對分散,且角度增大,沖蝕影響區域會擴散但沖蝕程度會迅速減弱;當彎管角度小于90°時,粒子對管壁的沖蝕效應會隨角度的減小而緩慢降低,但是沖蝕發生的主要區域位置不會發生較大改變。