兩相流與溶液化學平衡相結合的地熱井結垢模擬研究*

岑繼文，蔣方明

兩相流與溶液化學平衡相結合的地熱井結垢模擬研究*

岑繼文1,2,3，蔣方明1,2,3?

（1. 中國科學院廣州能源研究所，先進能源系統(tǒng)研究室，廣州 510640；2. 中國科學院可再生能源重點實驗室，廣州 510640； 3. 廣東省新能源和可再生能源研究開發(fā)與應用重點實驗室，廣州 510640）

水熱型地熱系統(tǒng)和干熱巖增強型地熱系統(tǒng)在運行一段時間后都有可能出現井筒和管道設備結垢現象，困擾和阻礙著地熱資源的高效低成本開發(fā)。由于溫度和壓力變化導致CO2從地熱水中逸出，進而使得CaCO3析出結垢，是地熱生產井堵塞的主要原因。將兩相流動換熱計算與水溶液物理化學模擬相結合，針對地熱井內CaCO3垢進行流動換熱與化學組分變化的模擬和計算方法的研究，獲得了特定地熱條件下井筒內流動狀況與化學組分變化情況的詳細分析結果，可為井筒地熱水流動進行結垢分析和預測。

地熱水；地熱井；結垢；計算機模擬；預測

0 引 言

2020年我國首次提出了要在2030年實現碳達峰，2060年實現碳中和的目標。由于傳統(tǒng)石化能源的不可再生性以及污染問題，大力發(fā)展新能源是實現“雙碳”目標的最佳途徑。地熱能作為可再生能源的一種，與太陽能、風能等可再生能源相比，具有資源儲量巨大、穩(wěn)定性相對較高等優(yōu)點。而地熱開發(fā)過程中，尤其是水熱型地熱資源開發(fā)過程中遇到的管道設備結垢問題，困擾和阻礙地熱資源的高效低成本開發(fā)。地熱電站中結垢類型主要為鈣垢、硅垢、金屬氧化垢及地熱流體中含有的微細顆粒物、污泥、微生物、硫化物等。其中鈣垢主要與壓力有關，通常發(fā)生在井筒內壓力降低而導致的閃蒸點附近，閃蒸使得大量二氧化碳析出，鹵水酸性降低，碳酸鈣過飽和而形成鈣垢析出[4]。另外，CaCO3結垢堵塞在地熱尾水回灌過程中是一個普遍的現象[2-3]，就單個井而言，結垢現象在井底、井上及管道設備中都有可能出現，使地熱能的效率大為降低。

劉明言[1]指出，目前在地熱流體的腐蝕與結垢趨勢預測方面的工作還相對較少，針對我國不同地區(qū)地熱流體特性開展系統(tǒng)的腐蝕結垢趨勢分析，以及三維多相傳遞和化學反應數值模擬研究是今后的方向。朱家玲等[5]對地熱水結垢腐蝕趨勢的判斷方法進行了探討，在地熱水的Cl?物質的量分數大于25%時，用雷諾茲指數判斷有局限性，此時，應結合拉伸指數進行水質判斷，才能使判斷結果符合實際。PáTZAY等[6]建立了CaCO3-H2O-CO2體系模擬地熱井結垢方法，采用了Davies和Pitzer活度計算方法，對井內流動壓降和溫度變化采用了簡單的線性化處理方法。ZOLFAGHARROSHAN等[7]用較為詳細的多相流方法模擬了井筒流動換熱情況，但對結垢的模擬僅僅用地面的成分情況進行判斷，并未全程對井筒進行物理化學組分變化的計算跟蹤。張恒等[8]利用美國地質調查局研發(fā)的PHREEQC軟件，采用水文地球化學模擬技術對康定某高溫地熱井結垢問題進行了分析研究。結果表明，該地熱井在開采過程中隨著地熱流體溫度、pH、壓力、氧化還原環(huán)境等條件變化發(fā)生嚴重的結垢現象。此外，除了水熱型地熱井，長期運行的干熱巖系統(tǒng)也可能會在井筒和地面設備中出現礦物結垢現象[9]。

地熱水結垢的基本物理化學機理已被研究得較為清楚，很多文獻已進行過詳細的論述。而在實際應用過程中，要想對結垢現象進行防治，需要在了解機理的基礎上建立較為精確的計算模型對其預測。本文針對地熱井內CaCO3垢形成情況進行詳細的兩相流動換熱模擬，并對地熱流體中化學組分的變化進行化學平衡計算模擬，以期獲得特定地熱條件下結垢情況的細致而精確的預測。

1 計算模型

1.1 井筒內CaCO3結垢基本物理化學機理

井筒結垢的形成基本機理是由于地熱水溶液上升、壓力不斷降低。盡管水溶液上升過程中地層溫度下降使得水溶液熱損失溫度有所降低，但由于壓力下降更為迅速對應的飽和溫度下降更為快速，因此在井筒某處會開始出現閃蒸，原先溶解在地熱水溶液中的酸性氣體CO2逸出進入蒸汽中，使得溶液pH升高，進而使CaCO3等溶解度達到飽和，析出固體形成結垢。

1.2 計算模型建立

井筒內CaCO3結垢原理雖已獲公認，但對其進行精確的建模計算具有相當的復雜性，涉及兩相流、地層換熱和溶液物理化學變化，尤其是閃蒸過程。不凝氣體組分和含量對結垢的影響問題涉及井筒的水溶液化學反應平衡和相平衡方面，以及地熱水溶液在管道中的兩相流流阻壓降和地層熱損失換熱等方面，本文參考了文獻[7,10-12]，建立如下計算模型。

1.2.1 兩相流計算模型

質量方程：

其中：為密度，m；為速度，m/s；為時間，s；為井筒長度，m。穩(wěn)態(tài)流動情況下沒有質量累積，因此有：

動量方程：

其中：為壓力，Pa；為摩擦系數，kg/(m?s2)，為井筒內徑，m；為井筒內截面積，m2；為重力加速度，m/s2；為井筒與水平方向的夾角，無量綱。聯(lián)合質量方程，整理可得：

如為兩相流，可將流體當作氣液均相混合物，井為豎直井，上式可變?yōu)椋?/p>

其中：下標mix表示混合參數；為豎直方向長度，m。根據不同的流型計算出不同的摩擦因子，需要用到占空比和混合密度。氣液兩相流可分為bubbly、slug、churn和annular四種流型[10]。

能量方程：

其中：為焓，J/kg；為單位質量流量單位管長內的換熱量，W/(kg/s)，可參考文獻[11]：

式中：k為有效導熱系數，W/(m?K)，其地熱井典型值為4 W/(m?K)。

1.2.2 化學平衡、相平衡模型

擬借鑒電解質相平衡模型[12]進行閃蒸過程計算。電解質溶液相平衡模型如圖1所示，含有一定量化學組分的來流溶液變化到指定的溫度、壓力之后，其化學組分在氣、液、固三相中重新分配。所涉及的方程如下。

物料衡算式：

其中：in,i和out,i分別為電離或離解反應前后各組分的摩爾流量，mol/s；為化學計量系數，無量綱；為反應程度，mol/s；為某個平衡方程；為總的平衡方程數量。

因此，化學反應平衡后更新的各組分組成為：

總物料衡算式為：

其中：為所有相中物料的總和，為液相物料，為氣相物料，單位均為mol/s；為液相組成，為氣相組成，為各相中的物料組成，均為無量綱；為總的物料的種類。

能量衡算式：

其中：L、V和F分別為液相、氣相和所有相總和的摩爾焓，J/mol；為換熱量，W。

相平衡：

其中：為亨利系數，無量綱。上式可理解為根據亨利定律，在一定溫度和平衡狀態(tài)下，氣體在液體中的溶解度（用摩爾分數表示）和該氣體的平衡分壓成正比。

化學反應平衡：

其中：為化學反應平衡常數，單位與具體反應方程式有關；為活度系數，無量綱；為逸度系數，無量綱；為逸度，Pa。液相活度系數的計算采用Pitzer電解質活度系數模型[13]，而氣相逸度系數采用 PR狀態(tài)方程[14]計算。

組分約束：

其中：m為氣相分子種類數，m表示分子。

電中性平衡：

此時的Z表示離子種類所帶電荷的真實值，而不是絕對值，無單位。

物性關聯(lián)式：

1.2.3 結垢判斷

其中：[Me]是二價陽離子；[An]是二價陰離子；結垢指數SI= 0表示平衡，SI > 0表示過飽和產生結垢；sp為溶解度。

對于常見的碳酸鈣結垢，可采用以下關聯(lián)式[7]：

式中：為溫度，℃；為壓力，單位為psi；離子強度為：

其中：C為各組分濃度，mol/kgw（kgw表示1 kg地熱水）。

1.2.4 氣泡點的判斷

氣泡的起始點可根據各不凝性氣體分壓加和是否超過了當地壓力來進行判斷：

1.3 模型求解算法

將井筒內地熱水流動簡化為一維流動，在長度方向上將井筒分為若干個單元，從底部向上依次逐步推進計算各單元參數，包括換熱、流動、組分變化、結垢情況等。單元計算流程如圖2所示。

圖2 單元計算流程圖

閃蒸算法是最為復雜的部分，可參考化工計算模擬的閃蒸分離過程算法[12]，該方法是對于包含個組分（所有離子和分子）并且發(fā)生個電離或水解等化學反應的兩相體系，根據電離平衡和相平衡關系得到一個+維方程組，并通過Marquardt迭代方法求解非線性方程組實現。

每個單元的計算最后都做一次飽和指數的計算，判斷是否開始結垢，如果結垢發(fā)生，則析出一定的量，使得溶液正好恢復當時條件下的飽和狀態(tài)，從而獲得各組分平衡關系。

1.4 化學反應方程和相變平衡方程

由于各地地熱流體成分的復雜性，不可能做到全部覆蓋，本文目的是闡述CaCO3結垢的計算模擬方法，計算程序中目前僅考慮CaCO3、CO2、H2O、NaCl、CH4這幾種物質的相互作用，在該物質體系下，涉及的化學反應方程和相變平衡見表1[15]。

表1 計算程序所涉及的化學反應方程和相變平衡

注：下角標aq表示在水溶液中，g表示氣相，s表示固相。

2 結果與分析

計算程序采用Python計算機語言，在Visual studio 2019集成開發(fā)環(huán)境中編寫了上述計算模型程序。使用Tkinter庫制作簡單的參數輸入界面，可方便地輸入各個算例參數。

目前常用的結垢模擬軟件有Wellsim、Phreeqc等，均為國外軟件，無法掌握內部計算模型細節(jié)，國內尚未有專門的地熱井結垢計算軟件。本文模擬計算工作的特點是可以考慮各種不凝性氣體的影響，兩相流計算可細化至不同流型的影響等。而且計算過程中使用到的各個計算模型，如兩相流壓降計算、活度系數、逸度系數模型、迭代模型等各方面均可隨時替換或優(yōu)化，所有參數均可根據需要做出調整。

2.1 算例分析

算例輸入參數：井筒長度為1 000 m、直徑為200 mm，管道表面粗糙度0.015，井底溫度70℃，壓力11.2 MPa，流量2.5 kg/s。輸入流體組分Ca2+濃度為4 mg/kgw，HCO3?濃度為12.18 mg/kgw，CO2濃度為50 mg/kgw，CH4濃度為200 mg/kgw。地溫梯度按照地面溫度25℃至井底線性變化計算。

將井筒劃分成500段（可任意設置），經模擬軟件計算后獲得各單元的氣?液組分組成、溫度、壓力、流速、pH、結垢飽和指數等參數，如表2所示。各氣體組分分壓為總壓乘以其氣相摩爾分數及逸度系數；pH可從H+濃度求出；CaCO3(s)累積結垢量為7.0 × 10?5mol/kgw；泡點/開始結垢點位置為118 m。因此，表2中的4個代表深度，分別為井底、井口、泡點前（150 m）和泡點后（60 m）各選取一點作為代表井內各部位參數的變化。

表2 井內地熱水不同深度下各組分濃度

圖3 （a）溫度壓力隨井深變化曲線；（b）泡點之后CO2分壓和pH值隨井深變化曲線；（c）泡點之后氣相組成隨井深變化曲線；（d）泡點后溶液中各組分變化情況

圖3為算例的輸出數據結果繪制的溫度、壓力變化，泡點之后的地熱水pH，氣相、液相組分以及CO2分壓變化曲線圖。在2.5 kg/s流量的地熱條件下，溫度下降較小，壓力呈線性下降。泡點之后由于CO2逸出，地熱水pH逐漸加速升高，CO2分壓呈現加速下降的模式，曲線變化趨勢符合碳酸鈣垢形成的機理。圖3c中氣體的組分變化表明CH4氣體比CO2有更快的增長速度，是由于CH4在水中具有較低的溶解度。圖3d為液相中各組分的變化曲線，可知溶液中除了游離CO2分子、CH4分子、HCO3?、Ca2+外，其他組分濃度較低，小于1 × 10?5mol/kgw。總之，泡點后相應數據逐步變化明顯，而泡點之前除總壓外其他參數幾乎不變。

2.2 不凝性氣體含量對結垢量和泡點位置的影響分析

圖4為以上述算例為基準，改變井底溫度和CO2含量兩個變量，獲得的不同地溫條件下CO2含量與最終結垢量的變化曲線。最終結垢量是根據CaCO3飽和指數，將其維持小于零來計算析出量。從圖中可以看出，在70℃地溫條件下，結垢量隨著CO2含量的增加而減少，增加到一定濃度后，結垢量趨近于0。而相同CO2含量下，隨著地溫的升高，結垢量更多，但由于Ca2+數量的限制，結垢量也有一定的最大值。

圖4 不同溫度下結垢量與CO2含量關系

圖5為上述算例條件下，分別改變每千克地熱水中CH4和CO2的含量對泡點位置的影響。從圖中可以看出泡點位置與這兩種不凝性氣體的含量呈現線性增加關系。一般情況下，泡點的位置即為結垢開始出現的位置。隨著不凝性氣體含量的增加，氣體更容易從水溶液中逸出，閃蒸點的位置會更深。

此外，需要說明的是本文算例是根據各類情況綜合考慮，通過改變相應參數條件來考察結垢模擬方法和結垢形成特性問題，并非某個特定井的具體參數。如需要對某個地熱井實際數據進行驗證對比，還需要根據該實際井添加相應的化學組分，進而在程序中關聯(lián)相應的化學反應方程、相平衡方程，結垢也不僅限于CaCO3垢，可能是硅酸鹽、硫酸鹽或多種成分都存在，氣體也不僅有CH4、CO2，還可能有H2S、NH3、N2等，這些都需要在程序中重新做針對性的修改和錄入相關參數。但模型方程和算法與上述一致，是通用的。上述算例中計算出的泡點位置與參考文獻[6]大致數量級相同，因此計算結果在合理范圍內吻合。

3 結 論

將兩相流計算模型與化學反應平衡模型以及化工閃蒸計算模型相結合，對地熱井內的流動進行精細模擬，獲得溫度、壓力、溶液及氣相中各組分的變化情況，分析其pH變化、氣相分壓、沸騰起始點，并進行結垢判斷等。經過試運算，模擬結果的數值較為合理，能夠測算出CO2分壓的指數性減小，pH的變化亦呈指數增大趨勢。同時，經過對比分析還得出以下兩個結論：（1）隨著CO2含量的增加，結垢量逐漸減少直至不結垢，而地熱水溫度越高越容易發(fā)生結垢；（2）隨著不凝性氣體的增加，泡點位置即結垢位置更深。

本文的地熱水化學成分組成僅考慮了Ca2+、HCO3?以及CO2和CH4不凝性氣體組分，真實的地熱井鹵水含有更多更復雜的組分組成，可在此模型框架的基礎上，增加更多物質組分，并盡可能利用更精準的化學反應平衡常數、亨利系數等化工常數，即可準確地模擬真實的地熱井流動和化學變化情況。

[1] 劉明言. 地熱流體的腐蝕與結垢控制現狀[J]. 新能源進展, 2015, 3(1): 38-46. DOI: 10.3969/j.issn.2095-560X. 2015.01.007.

[2] 云智漢, 馬致遠, 周鑫, 等. 碳酸鹽結垢對中低溫地熱流體回灌的影響——以咸陽地熱田為例[J]. 地下水, 2014, 36(2): 31-33. DOI: 10.3969/j.issn.1004-1184.2014. 02.012.

[3] 田濤, 陳玉林, 姚杰. 地熱回灌系統(tǒng)水結垢預測[J]. 地下水, 2011, 33(6): 27-28, 91. DOI: 10.3969/j.issn.1004- 1184.2011.06.012.

[4] 蔡正敏, 李剛, 李源, 等. 肯尼亞地熱電站結垢問題的日常維護[J]. 科技視界, 2018(25): 41-43. DOI: 10.19694/ j.cnki.issn2095-2457.2018.25.018.

[5] 朱家玲, 姚濤. 地熱水腐蝕結垢趨勢的判斷和計算[J]. 工業(yè)用水與廢水, 2004, 35(2): 23-25. DOI: 10.3969/ j.issn.1009-2455.2004.02.007.

[6] PáTZAYG, STáHLG, KáRMáNF H, et al. Modeling of scale formation and corrosion from geothermal water[J]. Electrochimica acta, 1998, 43(1/2): 137-147. DOI: 10.1016/S0013-4686(97)00242-9.

[7] ZOLFAGHARROSHAN M, KHAMEHCHI E. A rigorous approach to scale formation and deposition modelling in geothermal wellbores[J]. Geothermics, 2020, 87: 101841. DOI: 10.1016/j.geothermics.2020.101841.

[8] 張恒, 胡亞召, 云智漢, 等. 水文地球化學模擬技術在康定某高溫地熱井結垢研究中的應用[J]. 新能源進展, 2016, 4(2): 111-117. DOI: 10.3969/j.issn.2095-560X.2016. 02.006.

[9] YANAGISAWA N, MATSUNAGA I, SUGITA H, et al. Temperature-dependent scale precipitation in the Hijiori Hot Dry Rock system, Japan[J]. Geothermics, 2008, 37: 1-18. DOI: 10.1016/j.geothermics.2007.08.003.

[10] HASAN A R, KABIR C S. Modeling two-phase fluid and heat flows in geothermal wells[J]. Journal of petroleum science and engineering, 2010, 71(1/2): 77-86. DOI: 10.1016/j.petrol.2010.01.008.

[11] GARG S K, PRITCHETT J W, Alexander J H. A new liquid hold-up correlation for geothermal wells[J]. Geothermics, 2004, 33(6): 795-817. DOI: 10.1016/j.geothermics.2004. 07.002.

[12] 韓莎莎. 含電解質過程模擬研究與開發(fā)[D]. 青島: 青島科技大學, 2019.

[13] 韓莎莎, 鄭俊強, 孫曉巖, 等. Pitzer活度系數模型研究與開發(fā)[J]. 當代化工, 2020, 49(1): 221-227. DOI: 10.3969/j.issn.1671-0460.2020.01.053.

[14] LI J, WEI L L, LI X C. Modeling of CO2-CH4-H2S-brine based on cubic EOS and fugacity-activity approach and their comparisons[J]. Energy procedia, 2014, 63: 3598- 3607. DOI: 10.1016/j.egypro.2014.11.390.

[15] DUAN Z H, LI D D. Coupled phase and aqueous species equilibrium of the H2O–CO2–NaCl–CaCO3system from 0 to 250oC, 1 to 1000 bar with NaCl concentrations up to saturation of halite[J]. Geochimica et cosmochimica acta, 2008, 72(20): 5128-5145. DOI:10.1016/j.gca.2008.07.025.

Geothermal Well Scaling Simulation by Combining Two Phase Flow with Chemical Equilibrium

CEN Ji-wen1,2,3, JIANG Fang-ming1,2,3

(1. Laboratory of Advanced Energy Systems, Guangzhou Institute of Energy Conversion, Chinese Academy of Sciences,Guangzhou 510640, China; 2. CAS Key Laboratory of Renewable Energy, Guangzhou 510640, China;3. Guangdong Provincial Key Laboratory of New and Renewable Energy Research and Development, Guangzhou 510640, China)

Both hydrothermal geothermal system and hot dry rock enhancement geothermal system may experience scaling problem in the wells or facility on the ground after an operation period. Scaling problems inhibit the efficient exploitation of geothermal resources. For carbonate scale in production wells, the temperature and pressure changes in the fluid cause the release of CO2from the geothermal solution, and lead to supersaturation of the solubility of carbonate calcium. In this paper, two-phase flow model and chemical reaction simulation were integrated together to simulate the scaling problems in production wells to reveal detail temperature and pressure of the flow states and the changes of chemical components simultaneously. The model can be used to be a scaling prediction and analytical tool for geothermal production wells.

geothermal water; geothermal well; scaling; simulation; prediction

2095-560X（2022）01-0020-07

TK521

A

10.3969/j.issn.2095-560X.2022.01.004

2021-10-20

2021-11-09

國家重點研究發(fā)展計劃項目（2019YFB1504104）；中國科學院 A 類戰(zhàn)略性先導科技專項項目（XDA21060700）

蔣方明，E-mail：jiangfm@ms.giec.ac.cn

岑繼文（1979-），男，博士，副研究員，碩士生導師，主要從事地熱能開發(fā)利用、制冷熱泵、電子散熱、電動汽車熱管理等方面的研究。

蔣方明（1973-），男，博士，研究員，博士生導師，主要從事電化學能量/動力系統(tǒng)、增強型地熱系統(tǒng)、微熱流體系統(tǒng)、燃料電池水、熱管理及高效節(jié)能技術/產品等研發(fā)工作。