鞏 亮 韓東旭 陳 崢 汪道兵 焦開拓 張 旭 宇 波

1.中國石油大學（華東）新能源學院 2. 北京石油化工學院機械工程學院 3. 西安交通大學動力工程多相流國家重點實驗室

0 引言

地熱能因具有清潔性高、運行穩(wěn)定和空間分布廣泛等優(yōu)點，逐漸成為世界各國重點研究和開發(fā)的新型清潔能源。國家能源局《關于促進地熱能開發(fā)利用的若干意見》[1]明確指出，到2025年全國地熱能供暖（制冷）面積比2020年增加50%，地熱能發(fā)電裝機容量比2020年翻一番；到2035年，地熱能供暖（制冷）面積及地熱能發(fā)電裝機容量力爭比2025年翻一番。地熱資源根據(jù)地質構造特征、熱流傳輸方式、溫度范圍以及開發(fā)利用方式等因素可分為淺層地熱能、水熱型地熱（地下熱水）和干熱巖3種類型[2]。其中干熱巖是指不含或僅含少量流體，溫度高于180℃，其熱能在當前技術經濟條件下可以利用的巖體[3]。經初步測算，我國地下3～10 km范圍內干熱巖資源折合標準煤860×1012t[4]，利用其中的2%，相當于我國2020年能源總消耗量的約3 400倍。我國青海共和地區(qū)已鉆探發(fā)現(xiàn)溫度高達236 ℃的高品質干熱巖[5]，干熱巖開發(fā)在我國具有廣闊的應用前景。

1 增強型地熱系統(tǒng)研究簡況

增強型地熱系統(tǒng)（Enhanced Geothermal System，EGS）是開發(fā)干熱巖型地熱資源的主要手段。EGS形成過程為：鉆注入井，經注入井向儲層內注入高壓水，迫使巖石開裂從而形成具有滲透性的多尺度人工縫網結構（圖1-a），之后依據(jù)人工縫網走向確定采出井位置，形成一注一采、一注多采或多注多采等開發(fā)模式。EGS采熱過程為：經注入井將水、CO2等低溫工質注入具有人工縫網的干熱巖儲層，與高溫巖體充分換熱后，經采出井采出，由地面發(fā)電系統(tǒng)將熱能轉化為電能（圖1-b）。據(jù)統(tǒng)計，全球在建及投運的EGS工程已達30個，其中14個實現(xiàn)了運行發(fā)電，目前尚在運行的有5處，總裝機容量為12.2 MW。整體來看EGS尚未實現(xiàn)規(guī)?；?、商業(yè)化運行[6]，究其原因為一系列瓶頸有待突破，主要有：①EGS開發(fā)過程中，缺乏有效的干熱巖人工造縫調控技術，導致人工裂隙單一且尺寸較大、流體短路循環(huán)，造成過早熱突破、采熱效率低；②EGS形成和采熱過程受滲流、傳熱、介質變形、水巖反應等多種因素的影響，地熱儲層內多尺度多場耦合發(fā)展規(guī)律和機理仍不明確；③EGS采出井舉升過程中高壓降引起的流體閃蒸問題，影響井內的流動換熱特性，制約井內熱流體的高效提??；④EGS地上發(fā)電模式眾多，但具體適用場景并不明晰，仍存在熱電轉換效率低的問題。

圖1 增強型地熱系統(tǒng)示意圖

因此，圍繞上述重大需求和瓶頸問題，筆者針對EGS開采干熱巖型地熱資源涉及的4項關鍵技術進行了綜述，包括干熱巖儲層人工壓裂技術、干熱巖開采數(shù)值模擬技術、井筒熱流體高效提取技術以及干熱巖地熱發(fā)電技術，以期為我國干熱巖高效開發(fā)相關研究的發(fā)展提供一定的指導。

2 干熱巖儲層人工壓裂技術

干熱巖儲層較為致密，滲透率極低，為保證注入流體能夠與儲層巖石進行充分熱交換并順利經采出井采出，一般需要通過儲層改造形成連通良好、導流能力高的流動通道，提高儲層滲透性。但在儲層改造過程中，易形成單一的高滲裂縫，造成流體流動短路、過早產生熱突破，影響EGS的可持續(xù)開發(fā)利用。因此，明晰人工縫網形成機制、精準預測壓裂裂縫以及有效調控縫網結構，是形成復雜人工縫網和提高EGS采熱效率的一大關鍵。下面將從干熱巖儲層改造方法、人工縫網形成機制和裂縫擴展預測模型3個部分總結干熱巖儲層人工壓裂技術方面的研究進展和發(fā)展趨勢。

2.1 干熱巖儲層改造方法

干熱巖儲層改造方法主要包括[7]：水力改造、化學改造和熱改造，表1對3種方法進行了簡要總結。水力改造又稱水力壓裂，指將低溫高壓的流體注入地層，在井底憋壓產生人工裂縫，從而提高地熱儲層滲透性的方法[8-10]，是最主要的儲層改造方法。注入流體一般為清水，也可以是高黏壓裂液，一般垂直于最小應力方向，巖石最易發(fā)生開裂和擴展。油氣開發(fā)中壓裂易造成巖石拉伸破壞，而EGS主要為水力剪切破壞，能夠依靠裂縫粗糙度使裂縫維持張開狀態(tài)，裂縫較寬但開度較小，有利于流體與儲層的充分換熱[11]，但這種自支撐機制并不總是有效，依然需要注入支撐劑，保證裂縫張開[12]。

表1 干熱巖儲層改造的主要方法表

化學改造，包括酸化和堿性化學刺激等手段，通過向井中注入酸性液體/堿性液體，溶解儲層的礦物質[13]，達到改善裂縫連通性和滲透率的目的。酸化過程一般需加緩蝕劑，減緩酸液對完井管柱的腐蝕。酸化過程通常包括前置酸、主體酸和后置酸3個階段[14]，前置酸一般用HCl溶解碳酸鹽類巖石礦物，主體酸用土酸（HCl/HF混合物）來溶解硅酸鹽類巖石礦物，后置酸一般用清水將管柱中的酸液頂替進入地層中。室內實驗結果表明[15]，12%HCl+5%HF的土酸可顯著改善花崗巖的滲透性，在最優(yōu)反應時間下，滲透率可提升4個數(shù)量級。微觀結構和酸巖反應離子探測結果表明[15]，土酸與長石和云母礦物反應較快，而與石英反應相對不活躍。因此，石英顆粒可充當壓裂支撐劑，使裂縫保持張開。Watanabe等[16]研制了一種新的化學改造方法，主要使用螯合劑來持續(xù)溶解特定的礦物，進而形成流動通道，研究結果表明該方法可使得儲層滲透率增大近6倍。Xu等[17]對比研究了高溫高壓條件下酸性液體和堿性液體對干熱巖樣品的溶蝕情況，實驗結果發(fā)現(xiàn)，酸蝕后巖石滲透率提高了1.62倍，而與堿性溶性反應后，巖心滲透率可提高2.45倍。環(huán)境掃描電鏡微結構分析結果表明，干熱巖與酸液和堿性液體反應后，巖心內部均產生了二次礦物，包括綠泥石、球形二氧化硅和蒙脫石等。但現(xiàn)階段化學改造主要用于地層近井附近，溶解巖石或裂縫里的充填物[18-19]，作用范圍有限，是進行儲層改造的次要手段。

熱改造，指在低于巖石破裂壓力時將冷流體注入高溫地層，隨著注入時間延長，地層的注水能力逐步得到提高[20-21]。熱改造機理包括以下2個方面：①注入的冷流體會引起巖石冷卻收縮，使裂縫張開；②地熱儲層冷卻引起的熱應力會促進裂縫擴展[22-23]，縫面剝落碎片可起到支撐裂縫的作用，使裂縫保持永久性張開，進而提升注水能力。武治盛等[24]利用超聲波檢測技術，實驗研究了裂隙充填干熱巖體熱沖擊后縱波波速演化規(guī)律，0 ℃水冷卻沖擊對巖石的損傷比室溫冷卻沖擊更強；對比了不同熱沖擊實驗條件下的波速降幅變化規(guī)律，由于膠結面附近的礦物顆粒尺度差異最大，熱破裂程度較強，含膠結面充填花崗巖的波速降幅比母巖、充填體大。Baujard等[25]研究了法國Rittershoffen EGS工程中花崗巖的熱性質和水力性質，由于初始生產指數(shù)較低，通過低排量注入冷水、局部化學注入和水力刺激改造，可極大地改善井間水力連通性，從而提高EGS取熱效果。目前，熱改造方法的研究與應用較少，一般作為儲層改造的輔助手段，與水力壓裂和化學改造方法聯(lián)合使用。

干熱巖一般埋藏較深且地應力各向異性較強，僅依靠上述常規(guī)儲層改造技術一般難以形成復雜的裂縫網絡，注水采熱時易產生快速熱突破現(xiàn)象，導致采熱效率較低。近年來，有學者提出暫堵轉向壓裂技術，實現(xiàn)人工縫網準靶向調控，是未來提高干熱巖采熱能力的一種新途徑[26]。該技術一般使用可降解纖維或變形顆粒作為暫堵材料，通過暫堵增加縫內凈壓力，從而逼迫裂縫轉向形成多條裂縫，以提高人工縫網的形成效率（圖2）。初步研究結果表明[26]，高溫條件下，干熱巖人工裂縫暫堵過程中會出現(xiàn)明顯的壓力波動，溫度越高，壓力波動次數(shù)越多。同時，高壓作用下縫內突破封堵的次數(shù)與壓力波動存在明顯的關系，可以通過提升排量、增加縫寬來提高暫堵劑對裂縫的封堵效果。

圖2 暫堵轉向壓裂示意圖

2.2 人工縫網形成機制

合理的人工縫網形成不僅依賴于相關調控技術，還需要明晰干熱巖所處環(huán)境下的巖石力學特性以及人工縫網形成機制。

干熱巖由于高溫高壓作用，巖石塑性增強，壓裂呈現(xiàn)剪切破壞模式，實驗研究結果表明[26]：溫度增加40 ℃，可使干熱巖的斷裂韌性減小30%。高溫將使得干熱巖的峰值強度、彈性模量等巖石力學特征參數(shù)出現(xiàn)顯著的劣化，從而加劇干熱巖內部微裂隙的產生。此外，干熱巖的破裂特征一定程度上受熱流過程的影響[27]，在干熱巖熱采過程中，注入冷水可使裂縫面冷卻，導致縫寬增加，這將影響人工裂縫擴展，或者誘導二次裂縫的產生（圖3）。熱沖擊形成的裂縫具有自相似特征，且熱沖擊形成新裂縫面所需的能量與特征長度、地應力等因素有關。實驗研究結果表明，熱破裂的閾值溫度為260 ℃，此溫度下發(fā)育大量裂隙或形成貫通網絡。熱應力誘導的裂縫開度是增大近井區(qū)域導流能力的主要原因。在高溫條件下，干熱巖將出現(xiàn)軸向劈裂、雙縫共軛剪切、多縫平行剪切、單縫剪切4種破壞模式[28]。Liu等[29]通過熱流固耦合三軸壓縮實驗系統(tǒng)，研究了干熱巖在不同破壞模式下水力傳導能力的演化特征，干熱巖縫面粗糙度和特征應力隨著圍壓增加而增加，復合破壞模式下干熱巖滲透率最大（11.4 μm2），“Y”形剪切破壞模式下滲透率次之（9.7 μm2），最小滲透率情況為單純剪切破壞裂縫（7.2 μm2），增加圍壓大小可減小滲透率。Liu等[30]利用真三軸水力壓裂實驗系統(tǒng)，研究了循環(huán)注入條件下干熱巖的裂縫起裂與擴展機制，水力裂縫的起裂可分為2個階段：僅2～4條簡單裂縫階段和大于4條的復雜裂縫網絡階段，水力裂縫的起裂與擴展主要受流體循環(huán)、循環(huán)注入時間和注入速率等因素影響。在低排量注入條件下，高頻循環(huán)下（循環(huán)時間為10 s、20 s）水力裂縫的起裂與擴展主要受注入壓力影響。干熱巖可壓性，即儲層形成復雜縫網的能力，Wang等[31]通過高溫高壓三軸試驗壓機，以聲發(fā)射為監(jiān)測手段，開展不同溫度條件下干熱巖可壓性評價的實驗研究，基于應力應變曲線獲得不同應力加載水平下的體積應變，據(jù)此建立了考慮熱效應的脆性指數(shù)計算模型，該脆性指數(shù)不僅可以描述初始加載階段的熱損傷微裂紋密度，還可以描述不同應力水平下的裂紋起裂、擴展和成核特征。在上述脆性指數(shù)模型基礎上，發(fā)展了考慮熱效應的干熱巖可壓性評價新方法，該方法可適用于單軸和三軸應力狀態(tài)情形。

圖3 注冷水誘導熱應力產生的二次裂縫圖

此外，干熱巖水力壓裂過程中需注入成千上萬立方米的壓裂液以形成復雜的裂縫網絡，流體與干熱巖間相互作用對干熱巖縫網擴展的影響研究較為缺乏。為揭示干熱巖在不同飽和流體、不同加載應力水平下的損傷斷裂特征，國內學者利用高溫高壓巖石力學三軸壓機，通過監(jiān)測干熱巖的超聲波波速和聲發(fā)射響應特征，開展了蒸餾水、納流乳液等不同流體飽和干熱巖的巖石力學性質實驗研究。干熱巖處于飽和流體狀態(tài)后，縱波波速比干燥干熱巖升高40%，并且裂縫損傷應力比、斷裂韌性和黏聚力等力學參數(shù)明顯降低。飽和流體干熱巖具有更高的聲發(fā)射率，累積聲發(fā)射曲線呈現(xiàn)多個臺階狀上升特征。實驗結果表明[32]：飽和流體可以增加干熱巖的孔隙壓力，促進干熱巖水力壓裂過程中復雜縫網的形成，且納米乳液比蒸餾水形成復雜縫網的能力更強。

近年來，國內學者借鑒物體“熱脹冷縮”物理思想，對干熱巖實施交變溫載處理，增強壓裂改造過程中干熱巖儲層人工縫網復雜度，并研究了交變溫載作用下干熱巖人工縫網擴展規(guī)律[33]。如圖4所示，交變溫載后形成縫網的分形維數(shù)較單純熱處理后明顯提升，隨溫度升高，裂縫分形維數(shù)也逐漸增大，在500 ℃時分形維數(shù)最大提高了20%。交變溫載后，水力壓裂過程中干熱巖的破裂壓力和裂縫擴展時間顯著降低，在500～600 ℃時，破裂壓力降幅最大為51%，說明交變溫載作用能夠降低巖心破裂強度，并加速巖心破裂，形成宏觀貫通裂紋。

圖4 交變溫載與單純加熱處理干熱巖樣品水力壓裂裂縫形態(tài)對比圖[33]

2.3 裂縫擴展預測模型

干熱巖中人工縫網的形成機制復雜，要實現(xiàn)對儲層縫網的準確預測，需要構建綜合考慮熱應力、裂縫內流體流動以及巖石應力等因素影響的裂縫擴展預測模型。目前，干熱巖水力壓裂模型主要包括二維模型、擬三維模型和全三維模型，各種模型進展情況詳述如下：

二維模型假定裂縫高度恒定，一般可得到縫長、縫寬、壓力等參數(shù)的解析解，較為著名的模型有PKN模型、KGD模型和Penny模型[34-39]。PKN模型[9-10]適用于縫長遠大于縫高的情形，它假設每個垂向截面是相互獨立的，相對于縫長方向的壓力來說，縫高方向上的壓力起主導作用，而KGD模型[36-37]適用于縫高遠大于縫長的情況，模型假設平面應變在水平方向上，即所有水平截面相互獨立，裂縫寬度在垂向上比水平方向變化的更加緩慢。實際上，由于巖石具有滲透性，水力壓裂過程中流體沿著縫面濾失，因此PKN和KGD模型需要考慮縫內液體濾失和流體存儲效應。1957年，Carter提出了流體濾失方程，隨后Nordgren將濾失效應和流體存儲效應引入經典的PKN和KGD模型中，彌補了原有模型的不足。Penny模型[38-39]也稱徑向裂縫模型，適用于水力壓裂形成水平縫的情形。

由于二維模型假設縫高恒定，難以適用于多層水力壓裂縫高垂向擴展的情況，因此擬三維模型應運而生。該模型假設縫長方向的延伸大于縫高方向的延伸，并將縫內流體流動簡化為一維或二維流動，水力裂縫可以三維延伸，各個截面上變形為二維平面應變狀態(tài)且相互獨立，因此計算量要比全三維模型低[37-40]。基于平衡高度假設，Simonson等[40]推導了垂向含有3個小層、對稱分布的擬三維壓裂模型；隨后，F(xiàn)ung等[41]得到了垂向含有多個非對稱小層的擬三維壓裂模型；Zhang等[42]考慮了層狀巖石材料屬性和地應力的差異性，構建了一種新的擬三維壓裂模型，垂向平面裂縫沿著橫向劃分為一定數(shù)目的位移不連續(xù)單元，每個單元內橫截面變形為平面應變，流體壓力可沿著垂向和水平向2個方向發(fā)生變化。

二維模型和擬三維水力壓裂模型僅能考慮平面裂縫擴展，無法考慮非平面裂縫轉向擴展。因此，國內外學者發(fā)展了全三維水力壓裂模型，該模型可以真實反映水力壓裂在三維空間的裂縫擴展情況，裂縫擴展可以同時在3個方向上進行且可發(fā)生空間扭曲，流體可以在2個方向流動，因此全三維模型的計算量也最大[43-44]。Settari等[45]發(fā)展了一種全三維水力壓裂模型，該模型耦合了垂直平面內三維水力裂縫內的兩相流、支撐劑運移和熱傳遞過程，此外還可模擬停泵后的裂縫閉合過程，能預測壓力的變化趨勢。Yan等[46]建立了一種新的全三維水力壓裂模型，該模型中裂縫內流體流動通過裂隙滲流來表征，基質中的流體流動通過孔隙滲流來表征，可以模擬任意復雜裂縫網絡的擴展。Wang等[47]發(fā)展了基于生死彈簧單元的暫堵本構模型，構建了干熱巖變形、基巖滲流、裂隙流完全耦合的三維計算模型，該模型可自適應建立復雜人工縫網結構。彈簧元在加入暫堵劑時可以自動激活，停止注入暫堵劑時彈簧元自動刪除，可模擬暫堵劑加入過程中縫內凈壓力的增加過程。

3 干熱巖開采數(shù)值模擬技術

干熱巖儲層含有多儲滲結構，包括孔隙和裂縫，其中裂縫尺寸從厘米級跨越到米級（甚至千米級），具有多尺度、強非均質特征；另外，低溫流體注入高溫儲層，破壞儲層內溫度場、壓力場、應力場和化學場，產生強烈擾動，熱流固化多場耦合效應顯著。因此，干熱巖開采屬于典型的多尺度多場耦合問題。筆者將從巖心尺度多場耦合模型、儲層尺度多場耦合模型和尺度升級方法3個部分總結干熱巖開采數(shù)值模擬技術方面的研究進展和發(fā)展趨勢。

3.1 巖心尺度多場耦合模型

在巖心尺度下可利用計算機圖像處理技術或者分形建模技術對巖石的孔隙結構與礦物組成特征進行精細刻畫，進而對巖石內部的流動傳熱和受力變形特征進行精確描述，方便定量研究各種微觀因素對巖石物理特性的影響。此外，巖心尺度數(shù)值模擬也可從微觀角度揭示宏觀現(xiàn)象，彌補傳統(tǒng)巖石物理實驗的諸多不足，揭示物理實驗無法展現(xiàn)的機理[48-49]，常見的巖心尺度多場耦合模型如表2所示。

表2 巖心尺度多場耦合模型表

巖心尺度中的力場求解涉及多孔介質內部應力較小時的連續(xù)變形過程，和內部應力較大時由完整材料到起裂以及裂縫擴展的非連續(xù)變形過程，大致可分為連續(xù)介質力學模型和離散力學模型。連續(xù)介質力學模型下，除裂縫以外其他區(qū)域中孔隙和固體骨架被假設為連續(xù)變形的塊體，通過有限元方法（Finite Element Method，F(xiàn)EM）對基巖中連續(xù)變形區(qū)域的靜力學平衡方程進行離散求解。對存在裂縫的非連續(xù)變形區(qū)域，其周圍應變梯度往往較大，通常需結合塑性變形理論和損傷力學理論求解[50-51]。但當巖心內存在復雜裂縫網絡時，由于裂縫周圍區(qū)域具有較強的非連續(xù)性，連續(xù)介質力學模型的模擬效果往往不佳[52]。相較連續(xù)介質力學模型而言，離散力學模型將材料視為離散粒子或塊體的集合（圖5），采用拉格朗日方法追蹤每個粒子的運動，更適于計算非連續(xù)性大應變。常見的離散力學模型有離散元法（Discrete Element Method，DEM）[53]和非連續(xù)變形分析法（Discontinuous Deformation Analysis，DDA）[54]。離散力學模型中時間步長較小，計算負荷較大，適于巖心尺度下的力學分析，隨著計算機性能的提升也可進行一些大尺度的破裂模擬[55-56]。隨著模擬技術的發(fā)展，連續(xù)介質力學模型和離散力學模型之間的界限已不再清晰，連續(xù)介質力學模型在非連續(xù)變形處與離散力學模型結合能夠形成優(yōu)勢互補。例如有限元/離散元耦合方法（Finite-Discrete Element Method，F(xiàn)DEM）[57-58]，即在FEM中的黏結單元失效后，兩側單元由黏結關系轉變?yōu)榻佑|關系，接觸力的計算采用DEM中的方法。

根據(jù)孔隙的描述方法，巖心尺度下的流動與熱質傳遞模型可分為以下2種類型：①孔隙重構模型，即在詳細刻畫巖石孔隙結構的基礎上，求解巖石內部裂縫、孔隙中的流動，再根據(jù)流體速度場求解孔隙內部的傳熱和傳質方程（圖6-a）；②孔隙網絡模型[61]（Pore Network Model，PNM），將巖石內部裂縫、孔隙的流動通道簡化為球形孔隙和喉道的流通關系，流體存儲在孔隙中，孔隙之間流體的質量交換由一維的喉道描述，流體與固體之間的相互作用如對流換熱、流固作用力等也需經喉道計算得到（圖6-b）。

圖6 巖心尺度下裂縫及孔隙的描述方法圖

在第一種類型中，可采用傳統(tǒng)的CFD方法如有限體積法（Finite Volume Method，F(xiàn)VM）、有限差分方法（Finite Difference Method，F(xiàn)DM）對孔隙空間中的Naiver-Stokes方程、傳熱和傳質的對流擴散方程進行求解，具有較高的計算效率和數(shù)值穩(wěn)定性，也適用于巖心多相流動中大密度比、黏度比的模擬。另一方面，近年來快速發(fā)展的格子玻爾茲曼方法（Lattice Boltzmann Method，LBM），由于在網格上近似求解時間、空間、速度離散化后的玻爾茲曼方程而具有粒子動力學內在特征，在處理復雜形狀的流固界面時具有顯著優(yōu)勢，LBM中多分布函數(shù)模型的發(fā)展也使其適合多場耦合方面的模擬[62]。在無網格方法中，光滑粒子流體力學（Smoothed Particle Hydrodynamics，SPH）方法近年來也廣泛應用于多場耦合模擬[63]，其核心是根據(jù)粒子之間的距離計算粒子的未知量場[64]。值得注意的是，LBM、SPH在多場耦合中的數(shù)值穩(wěn)定性較差，容易產生數(shù)值振蕩。對于孔隙網絡模型而言，真實裂縫和孔隙的復雜幾何結構被簡化，因此具有很高的計算效率，是研究多孔介質中熱質傳輸規(guī)律的有效辦法。因為孔隙和喉道的尺度大小和空間分布存在多種選擇方式，且與多孔介質的結構特征存在密切關系，如何從巖心的CT數(shù)字掃描中提取出合理的孔隙網絡模型是該類問題研究的難點和熱點之一[65]。具體到流動模擬時，一維喉道中的流動通常采用達西或泊肅葉穩(wěn)態(tài)流動進行近似計算，孔隙之間的流動總體滿足質量守恒，并且孔隙壓力波動由孔隙中流體體積變化引起[66-67]。另外，在計算喉道的傳熱和傳質方程時，為描述流固界面的信息交換，不可避免地需要引入簡化后的經驗公式。

當力場模型和流動與熱質傳遞模型耦合計算時，上述方法兩兩之間結合可以產生多種組合模型，如圖7所示，即流固（Hydro-Mechanical）[55-56]、熱固（Thermal-Mechanical）[52,68]、熱流固（Thermal-Hydro-Mechanical）[50,59,69-70]等巖心尺度模型。化學反應會造成裂縫通道、孔隙結構的變化，其對巖石力學特征影響較大，但缺乏有效數(shù)值手段表征，因此目前有關巖心尺度下力場與化學場耦合的數(shù)值模型如流固化（Hydro-Mechanical-Chemical）、熱流固化（Thermal-Hydro-Mechanical-Chemical）模型的相關研究較少。

圖7 巖心尺度熱—流—固—化多場耦合圖

3.2 儲層尺度多場耦合模型

儲層尺度的模擬區(qū)域更大、時間更長、更符合工程需求，對于整個采熱系統(tǒng)的取熱性能預測和評價具有重要意義。儲層尺度模擬主要涉及裂縫型儲層中裂縫的描述、多場耦合控制方程與求解方法2個方面。

裂縫型儲層的表征方法主要分為3類：等效介質方法、雙重孔隙介質方法、離散裂縫類方法（圖8[71]）。等效介質方法不直接描述巖石中的裂縫和孔隙，而是將裂縫和孔隙中的孔隙度和滲透率等效為巖石的原生孔隙度和滲透率，之后采用連續(xù)介質對應的一系列理論進行建模[72-74]。對于巖石總體近似表現(xiàn)為具有對稱滲透率張量的均質多孔介質，該方法可以得到符合實際的結果，然而對于裂縫主導的流動，該方法忽略了巖石基質和裂縫介質間的流量交換，模擬結果往往誤差較大。雙重孔隙介質方法認為基質原生孔隙與天然裂縫相互連通，同時基質和裂縫即是流體的儲存場所又是流動通道，但模型中基質被分割為正方體、圓柱體、球形、層疊等規(guī)則形狀[75]，不能充分體現(xiàn)裂縫的各向異性和不連續(xù)特征。相較上述2種方法而言，離散裂縫類模型以顯式表征裂縫為核心，再結合連續(xù)介質力學和裂縫滲流力學等理論進行求解，可以提高儲層中裂縫的計算精度。根據(jù)物理情形不同，離散裂縫類模型又可分為離散裂縫模型（Discrete Fracture Model，DFM）[76-77]、離散裂縫網絡模型（Discrete Fracture Network，DFN）[78]、嵌入式離散裂縫模型（Embedded Discrete Fracture Model，EDFM）[79-80]。但離散裂縫模型滿足高精度要求的同時，存在網格劃分困難、計算量大等問題，難以用于具有復雜裂縫結構的三維實際干熱巖多場耦合數(shù)值模擬。

圖8 裂縫巖體表征示意圖（據(jù)李庭宇[71]修改）

基于上述介紹的裂縫描述方法，已開發(fā)出很多儲層尺度下熱流[81-82]、流固[83-84]、熱流固[85-86]、熱流化[87-88]、熱流固化[89-90]耦合模型，此處對熱流固化四場中的經典控制方程和其之間的部分耦合關系進行介紹。

巖石中力場求解多基于線彈性假設的靜力學平衡方程，在孔隙壓力、熱應力的影響下應力應變本構方程為：

式中，σc表示柯西應力張量，Pa；εc表示應變張量；λ和G表示拉梅系數(shù)，Pa；α表示Biot系數(shù)；β表示熱膨脹因子，Pa/℃；I表示單位張量；p0表示參考壓力，Pa；T0表示參考溫度，℃。

基巖和裂縫中的質量守恒方程分別如下式所示：

式中，上標m和cr分別表示基巖和裂縫系統(tǒng)；p表示孔隙壓力，Pa；ρf表示流體密度，kg/m3；φ表示孔隙度；ct表示壓縮系數(shù)，Pa-1；u表示達西速度，m/s；Qm-cri表示裂縫i到基巖的網格體積平均后的流體質量交換，1/s；Qcri-m表示基巖到裂縫i的網格體積平均后的流體質量交換，1/s；Qcri-crj表示裂縫j到裂縫i的網格體積平均后的流體質量交換，1/s；QB表示注入井或抽出井造成的源項，1/s。

基巖和裂縫中的壓力場和滲流場之間滿足達西定律：

式中，μ表示孔隙流體的黏度，(N·s)/m2；K表示滲透率；g表示重力加速度，m2/s。

基巖和裂縫中的能量守恒方程分別如式（5）和式（6）所示，其中基巖內部流體和固體骨架之間被假設為局部熱平衡狀態(tài)。

式中，T表示溫度，℃；cp表示比熱容，J/(kg·℃)；Γ表示熱擴散系數(shù)，m2/s；下標eff表示孔隙水和固體骨架綜合考慮后的等效參數(shù)；E表示能量交換，J/(s·m3)；EB表示注入井或抽出井造成的源項，J/(s·m3)。

基巖和裂縫中的溶質運移方程分別為：

式中，C表示溶質濃度，mol/m3；D表示溶質擴散系數(shù)，m2/s；HB表示注入井或抽出井造成的源項，mol/(m3·s)；γ表示溶質濃度交換，mol/(m3·s)；HC表示化學反應產生的源項，mol/(m3·s)，其與水巖反應速率有關，水巖反應速率可以表示為：

式中，kC表示化學反應速率常數(shù)，mol/(m2·s)；KC表示反應平衡常數(shù)；AS表示反應比表面積，m2/kg；γ表示離子活度積；指數(shù)θ和η與化學反應級數(shù)有關。

化學反應速率常數(shù)采用Arrhenius公式描述：

水巖反應會對基巖和裂縫中的孔隙度產生影響：

式中，φ0表示初始孔隙度；ρs表示固體礦物密度，kg/m3；Ms表示固體礦物摩爾質量，kg/mol。

孔隙度變化造成滲透率變化，采用Kozeny-Carman方程描述基質中孔隙度（φ）和滲透率（K）之間的關系：

式中，K0表示初始滲透率，m2。

孔隙度φ的變化進一步造成反應比表面積AS的變化，它們之間的關系為：

式中，A0表示初始反應比表面積，m2/kg。

上述控制方程的求解方法眾多，常見的傳統(tǒng)方法有FEM、FVM、FDM等，不同方法之間的原理和所得到的離散方程差別較大。近年來，擴展有限元方法（Extended Finite Element Method， XFEM）和擴展有限體積法（Extended Finite Volume Method，XFVM）[91]被提出，其核心思想是在裂縫所在單元的節(jié)點處引入增強函數(shù)來描述單元內部裂縫面處的非連續(xù)變形，同時這些增強函數(shù)不會影響傳統(tǒng)FEM中形函數(shù)的特性。XFEM相較于傳統(tǒng)FEM的主要優(yōu)點是無需在非連續(xù)變形區(qū)域加密網格或者采用大量的非結構化網格。XFEM一經提出便被廣泛關注且快速發(fā)展[92-93]，著名商業(yè)軟件ABAQUS也植入了XFEM方法，可用于進行非交叉裂縫的模擬。XFEM缺點是增強函數(shù)引入后整體離散方程容易接近線性相關從而增加了收斂難度，同時裂縫分布較多時大量裂縫貫穿單元增加了求解自由度使得計算耗時較大。借鑒XFEM的部分思想，Deb等[94-95]提出了基于有限容積法離散非連續(xù)變形問題的XFVM方法，其可以將XFEM中裂縫貫穿單元增加的8個自由度降低至4個，從而適于快速求解。XFVM求解力場并結合FVM求解THC控制方程可以使得流動和力學變量處于相同網格節(jié)點，易于通用式編程和網格量較大時的并行計算[71]。

3.3 尺度升級方法

由于儲層中裂縫和孔隙存在明顯的多尺度特征并且不同尺度的裂縫和孔隙在儲層模擬中表現(xiàn)出的規(guī)律和影響不同，尺度升級方法成為儲層描述的新概念和新發(fā)展方向。截止目前，針對裂縫型儲層尺度升級方法的研究正處于發(fā)展階段且多為HM耦合的水力壓裂模型，按照粗、細尺度之間信息交換的方法主要分為計算均勻化方法、連接尺度方法、多尺度有限元法、多尺度有限體積法等。

計算均勻化方法是建立不同尺度聯(lián)系時最常用的方法，在尺度連接時多基于平均場理論，假設在粗尺度和細尺度下由微元變形造成的虛功是相等的來進行求解[96]。粗尺度下采用FEM方法求解，在有限單元的每個高斯積分點處設置一個可以反映該位置處微觀結構下受力、流動等表征屬性的細尺度單元，通過對細尺度單元界面處的邊界值問題的計算求解獲得單元內部響應，之后進一步獲得積分點處的應力—應變張量、滲透率等參數(shù)用于粗尺度求解。在細尺度單元邊界值處理時，前人研究表明當采用均勻變形、零轉動以及流動傳熱中的Dirichlet邊界條件時會高估表征屬性，而均勻受力、自由轉動以及流動傳熱中的Neumann邊界條件會低估表征屬性，采用周期性邊界條件可以得到處于兩者中間較合理的結果[97-99]。計算均勻化方法對于裂縫、巖性分布均勻的儲層，能很好地融合宏觀和微觀的特征，摒棄宏觀唯象本構假定，根據(jù)微觀結構特性估計其宏觀的等效特性[100-101]。此外，細尺度單元可以采用多種方法離散求解，常見的有FEM方法、DEM方法[102]等（圖9）。

圖9 采用計算均勻化方法的多尺度建模技術圖

連接尺度方法最早是由Wagner等[103-104]在將分子動力學模型和宏觀連續(xù)體模型連接計算時所提出，其將分子動力學的計算區(qū)域設定在局部小部分區(qū)域，而在其他絕大部分區(qū)域中依然采用宏觀連續(xù)體模型計算。該方法將節(jié)點位移未知量分解為相互解耦且正交的粗尺度和細尺度兩個部分，粗尺度的解由FEM形函數(shù)插值表示，而細尺度的解被表示為偏離插值的一部分，由離散力學模型如DEM求解[105]。連接尺度之間的投影算子是該方法的關鍵，其準則是使總動能的粗細尺度部分解耦，即總動能或總動力虛功可分解為粗、細兩種尺度的動能或動力虛功之和，進而解耦求解表示顆粒材料的運動行為。萬柯等[106-107]提出了一種基于DEM-FEM的連接尺度法模型，其在多尺度交界面處提出了高效的準靜態(tài)界面條件，并基于此模型計算了邊坡穩(wěn)定問題[108]，其中在邊坡上方大變形的梯形區(qū)域內采用了DEM粒子計算（圖10）。此外，連接尺度方法也已成功應用于材料動態(tài)斷裂[109]、流體流動[110]的多尺度模擬中。在儲層人工壓裂和干熱巖開采過程中，可在裂縫周圍采用DEM求解來模擬非連續(xù)變形和裂紋擴展，即進行細尺度表征，在熱響應較小的連續(xù)變形區(qū)域采用FEM方法求解，即進行粗尺度表征，在兩者界面處采用連接尺度方法，進而實現(xiàn)高效準確的儲層模擬。

圖10 采用連接尺度方法的多尺度建模技術圖（據(jù)Wan等[108]修改）

多尺度有限元法、多尺度有限體積法的核心思想在于通過對單元內細尺度局部子問題的數(shù)值模擬，得到粗尺度下當?shù)貑卧獌鹊幕瘮?shù)，進一步在粗尺度下基于獲得的基函數(shù)對控制方程離散求解[111]。這些基于細尺度獲得的基函數(shù)不同于傳統(tǒng)方法中采用多項式表示的基函數(shù)，其可以準確描述材料粗單元內部的非均質性，并可以自動地將細尺度下解的信息代入到粗尺度范圍內，得到粗尺度下較為準確的結果。在多尺度有限元法中，若粗尺度單元尺寸與非均質尺寸相近會產生共振效應降低計算精度，此時需采用超樣本技術[112]，即在大于粗尺度單元尺寸的超樣本中求解臨時基函數(shù)。多尺度有限元法可以較好地捕捉小尺度下大變形問題，進而解決傳統(tǒng)FEM模擬裂紋擴展時網格重構和加密問題，適用于研究金屬、巖石等細觀破裂損傷評估[113]。相較而言，多尺度有限體積法通過全局守恒的通量進行計算，求解自由度低，因此被廣泛應用于多孔介質流動換熱、流固耦合問題中[114]。此外，若細觀、宏觀尺度求解控制方程相同時，多尺度有限元法、有限體積法較直接細尺度網格求解均大大降低了計算負荷。

4 井筒熱流體高效提取技術

地熱開采過程中，注入的冷流體在儲層裂隙內流動并與高溫巖石進行充分換熱后，經過采出井傳輸至地面，進行后續(xù)溫泉、供暖、發(fā)電等應用。EGS中，流體與巖石換熱后進入采出井時的溫度較高，同時由于井筒深度可達3～10 km，流體經井筒傳輸至地面的過程中將產生幾十MPa的壓降。在采出井口附近，若壓力降低至流體飽和壓力之下，將引發(fā)流體的閃蒸相變現(xiàn)象，極大地影響井內熱流體的高效提取。因此，明晰地熱井內流體閃蒸過程的流動換熱特性，實現(xiàn)對井內流體閃蒸的精確預測，并進行有效預防是實現(xiàn)地熱井內熱流體高效提取的關鍵。筆者將從地熱井內流體閃蒸的原理，以及井內閃蒸特性研究的實驗和數(shù)值模擬方法3個部分，總結井筒熱流體高效提取技術方面的研究進展和發(fā)展趨勢。

4.1 地熱井內的閃蒸

閃蒸現(xiàn)象主要可以分為以下2種：一種是將過冷流體的壓力迅速降低到初始溫度所對應的飽和壓力之下，使其由過冷狀態(tài)瞬間轉變?yōu)檫^熱狀態(tài)，從而產生劇烈相變的閃蒸現(xiàn)象[115]；另一種閃蒸現(xiàn)象主要發(fā)生在循環(huán)系統(tǒng)中，由于流體在較長的上升段內舉升時，其飽和溫度隨著流體靜壓的減小不斷減小，當流體溫度大于當?shù)仫柡蜏囟葧r將發(fā)生閃蒸[116]，產生閃蒸點以下為單相流體，閃蒸點以上為氣、液兩相流體的流動狀態(tài)。地熱井筒內流體的閃蒸現(xiàn)象與后者類似。

流體中的熱不平衡程度和汽化核心是誘發(fā)閃蒸汽化的充分條件和必要因素[117]，汽化核心可以是流體內部的雜質、不凝性氣體或者從流動壁面產生。在流體舉升過程中，一旦發(fā)生閃蒸，大量的流體顯熱轉化為蒸汽潛熱釋放[118]，汽化核心處的流體產生汽泡形成泡狀流，隨著流體向下游流動過程中的靜壓逐漸減小，閃蒸穩(wěn)定流動狀態(tài)下的流型由上游的泡狀流逐步向下游的彈（帽）狀流和攪混流發(fā)展[117]（圖11）。由于兩相流動狀態(tài)的換熱系數(shù)要大于單相流動狀態(tài)[119]，若井筒沒有良好的保溫措施[120]，井內流體閃蒸的發(fā)生將增大井筒與地層之間的換熱損失，降低熱流體出口溫度，從而導致地熱井的開采熱量減產。此外，閃蒸的發(fā)生會導致井筒內流體的壓力發(fā)生劇烈波動，誘發(fā)系統(tǒng)產生流動震蕩[121]，甚至誘發(fā)水擊現(xiàn)象[122]，威脅系統(tǒng)的安全性。

圖11 閃蒸穩(wěn)定流動狀態(tài)下的兩相流型圖[117]

EGS和水熱型地熱系統(tǒng)均為開式系統(tǒng)，即從地層內部抽采出熱流體，井內流體含有大量Ca2+等礦物質離子及二氧化碳等不凝性氣體。當井內流體發(fā)生閃蒸汽化后，氣相水蒸汽增多，導致二氧化碳分壓降低，溶解在液相的二氧化碳隨之逸出，造成最常見的碳酸鈣結垢[123]（圖12）。目前，美國[124]、匈牙利[125]以及中國西藏羊八井[126]、甘孜[127]、那曲[128]等地的中高溫地熱發(fā)電站地熱井均有井筒結垢問題的相關報道。嚴重的井筒結垢問題將導致地熱井的減產和關閉，制約地熱能的可持續(xù)開發(fā)利用。墨西哥Los Humeros地熱井由于井筒碳酸鈣結垢問題，蒸汽產量不斷下降，由投產時的38 t/h降至4 t/h[129]，羊八井、那曲、川西地區(qū)的地熱井因為嚴重的結垢問題，需要長期依靠機械除垢來維持正常運行，甚至直接停產[128]。目前，國內外學者已經針對地熱采出井筒內流體的閃蒸問題展開了大量研究，但由于采出井高井深的結構特點以及大壓降、閃蒸相變、兩相流動等復雜的物理特征，相關研究存在著實驗研究難以有效展開、數(shù)值模型不夠精確等難點。

圖12 地熱井筒及閥門的碳酸鈣結垢圖[123,127]

4.2 流動系統(tǒng)中的可視化閃蒸實驗

針對3～10 km的地熱井筒開展全尺度的實驗研究顯然是不切實際的，目前以地熱井筒為實驗對象，開展井筒內流體流動換熱特性以及閃蒸現(xiàn)象研究的實驗鮮有報道。值得關注的是，地熱井筒內流體的閃蒸過程與自然循環(huán)系統(tǒng)中的類似，根本區(qū)別在于循環(huán)的驅動力不同，自然循環(huán)依靠流體自身密度差驅動循環(huán)，而地熱井筒內流體的流動主要依靠泵功。

關于自然循環(huán)系統(tǒng)中流體閃蒸實驗的研究目前廣泛存在于核反應堆安全系統(tǒng)的設計中[130]，主要的實驗思路為[117]：冷流體首先進入加熱段，通過電加熱等方式加熱至設定溫度；流體經充分加熱后進入上升管路模擬自然循環(huán)中的流體舉升過程，管路一般采用透明材料，結合高速攝像機可以開展可視化研究，管路沿程布置測溫點、測壓點測量流體的溫度、壓力參數(shù)；流體經上升段進入上部水箱，后經下降段進入加熱段完成循環(huán)，沿程可布置電磁流量計等測量循環(huán)流量（圖13）。

圖13 自然循環(huán)系統(tǒng)中的可視化流體閃蒸實驗平臺圖[117]

國內外學者基于上述實驗過程開展了大量研究，包括對流體閃蒸流型演變的可視化研究、閃蒸的穩(wěn)定性研究、閃蒸的流動換熱特性以及對閃蒸進行控制等。Manera等[131]首次搭建出了垂直管道的流體閃蒸實驗平臺，并對管內靜止狀態(tài)與流動狀態(tài)的流體瞬態(tài)閃蒸過程進行了可視化實驗，采用絲網傳感器實現(xiàn)了對氣泡尺寸分布的測量，探究了不同流速下流體閃蒸流型的變化。閃蒸的發(fā)生將引起系統(tǒng)內部強烈的不穩(wěn)定性流動，大大影響系統(tǒng)內部的阻力特性，郭雪晴等[132]通過可視化實驗分析了自然循環(huán)系統(tǒng)中的瞬態(tài)運行特性和不穩(wěn)定性機理，發(fā)現(xiàn)流體的降壓閃蒸是影響系統(tǒng)穩(wěn)定的流動驅動壓頭的一大根本原因。在流動特性研究的基礎上，通過測量系統(tǒng)的流量和溫度參數(shù)，可以對流動過程中的閃蒸換熱特性進行研究，李魯寧等[133]探究了不同過熱度、壓力、流量等工況下閃蒸換熱系數(shù)的變化規(guī)律，張友森等[118]通過搭建循環(huán)閃蒸實驗臺研究了循環(huán)閃蒸顯熱轉化率受流量、壓力、液膜高度的影響規(guī)律。此外，在了解閃蒸流動換熱特性的同時，預防閃蒸的發(fā)生也是備受關注的方向，研究表明，通過提升系統(tǒng)壓力或減小阻力可以有效減少閃蒸現(xiàn)象的發(fā)生，使系統(tǒng)流動趨于穩(wěn)定[134-135]。程俊等[136]還通過在上升段內設置內插物使系統(tǒng)流型實現(xiàn)了由流動不穩(wěn)定狀態(tài)下的間歇流向穩(wěn)定流型的轉變，有效地減小了系統(tǒng)的震蕩幅度。上述研究可作為地熱井內流體閃蒸流動換熱特性研究以及預防井內流體閃蒸的參考，若要準確把握地熱井內流體閃蒸的特征，仍需進一步開展地熱井內的流體流動閃蒸實驗。

參考自然循環(huán)系統(tǒng)中的閃蒸實驗，要開展地熱井筒內閃蒸過程的實驗需要解決以下2個問題：①改變流體循環(huán)的驅動方式；②實現(xiàn)對高井深、大壓降特點下地熱井筒內流體流動狀態(tài)的等效實驗還原。Ishii等[137]根據(jù)小擾動原理，在單相?；椒ɑA上，應用漂移流模型推導出兩相穩(wěn)態(tài)自然循環(huán)模化準則，意圖構建合理的小比例?；瘻蕜t，在準確反應原型系統(tǒng)流動和傳熱特性的前提下，減小實驗系統(tǒng)的幾何尺寸和運行參數(shù)。郭雪晴等[138]對大尺度開式自然循環(huán)系統(tǒng)內的閃蒸過程進行了模化分析，所得模型系統(tǒng)的運行特性與原型實驗結果吻合較好。由此可以推斷，通過合理的模化準則來進行地熱井筒閃蒸實驗的設備尺寸設計及運行參數(shù)設計，即通過小比例的實驗來實現(xiàn)對大尺度地熱井筒內流體閃蒸過程的模擬，在設計上具有一定的可行性。

4.3 地熱井內流動傳熱的數(shù)值模型

相較于實驗研究，數(shù)值模擬方法在地熱井筒方面的應用更為廣泛，主要可以分為穩(wěn)態(tài)下井內流動的模擬和瞬態(tài)下井內流動的模擬2個方面。

考慮到地熱井筒存在高度方向尺寸與徑向尺寸比例極大的結構特征以及井筒內可能存在的閃蒸相變及復雜兩相流等問題，在進行數(shù)值模擬計算過程中往往需要引入一些假設：①將井筒簡化為沿高度方向的一維物理模型，忽略徑向的參數(shù)變化或將三維方程沿徑向截面積分簡化為一維的守恒方程[139]；②假設井筒內的流動狀態(tài)為穩(wěn)定狀態(tài)[140]，忽略時間項導數(shù)；③假設井內為單相流體，對于兩相的模擬常采用均相模型[141]，將流體看作是一種偽單相流體。對于穩(wěn)態(tài)流動的模型，流量參數(shù)可以通過邊界條件獲得，在整個模擬計算過程中，只需要考慮兩個主要變量即可求解整個模型。求解方法可以是從井口邊界到井底的自上而下的求解（TOP-DOWN），也可以是從井底邊界到井口的自下而上的求解（DOWN-TOP）[142]。壓力、溫度、流動質量分數(shù)或含氣率常作為主要的變量用于穩(wěn)態(tài)模型的求解[143]，其中質量分數(shù)無法用于含逆流的模型求解，因為當質量流量為零時流體質量分數(shù)是無法定義的，此時可選擇采用含氣率來求解相關模型[144]。壓力和流體的焓也是常用的主要變量，當流體不存在相變時，流體的焓可以看作是溫度、壓力的函數(shù)，但地熱井內存在閃蒸相變時，這種方法就不可行了，需要對方程進行離散，根據(jù)邊界條件逐步求解各節(jié)點的壓力、焓、流體質量分數(shù)等參數(shù)[140]?；谏鲜龇椒ǎ瑢W者們采用一維穩(wěn)態(tài)模擬實現(xiàn)了井筒內流體沿程速度、壓力、溫度以及含氣率等參數(shù)的快速計算，完成了包括采出井出口參數(shù)和井底儲層環(huán)境參數(shù)的反演、井內流體閃蒸位置的預測及預防結垢等問題的模擬研究。卜憲標等[145]針對地熱井內的結垢問題，采用DOWN-TOP的求解方法，根據(jù)已知的壓力、溫度、總質量流量、二氧化碳含量等井底參數(shù)對沿程壓力、溫度分布進行了計算，確定了井內流體的閃蒸位置，明確了影響井內閃蒸和結垢的主要參數(shù)。梁海軍等[146]采用DOWN-TOP算法，計算獲得了井筒沿程壓力、溫度分布，反演得到井底流體的溫度參數(shù)，預測了流體閃蒸的起始位置，并提出通過調控井口壓力和流量，可以調節(jié)閃蒸位置，優(yōu)化流量和防垢成本之間的關系。表3給出了主要的地熱井筒模擬器，Bjornsson等[147]開發(fā)出了第一個可以處理多進料層的穩(wěn)態(tài)地熱井筒模擬器HOLA，自此之后，學者們經過對模型的不斷修正和開發(fā)，又相繼提 出 了 GWELL[148]、WELLSIM[149]、MULFEWS[150]、FLOWELL[151]、SWELFLO[152]、PROFILI[153]等 一 系列求解器，可以滿足多井、蒸汽井等不同地熱生產條件下井筒流動傳熱過程的模擬計算。但閃蒸作為一種井內流體相態(tài)突變的現(xiàn)象，采用穩(wěn)態(tài)模型計算不具備足夠的說服力，因此對于井內流體閃蒸流動的模擬需進一步結合瞬態(tài)模型進行。

表3 主要的地熱井筒流動模擬器統(tǒng)計表

對于瞬態(tài)下井筒內的流動模擬，Miller[154]在1980年開發(fā)出了第一個瞬態(tài)井筒模擬器WELBORE，后續(xù)相繼有學者發(fā)展出用于求解地熱井內兩相流動的三方程模型，Akbar等[157]采用相關模型實現(xiàn)了高焓值地熱井內閃蒸過程中流體性質的預測，同時考慮了浮力、相變、壓縮性、熱相互作用、壁面摩擦以及相間滑移等因素的影響。近來，Tonkin等[142]總結歸納了大量關于地熱井筒流動的穩(wěn)態(tài)和瞬態(tài)數(shù)學模型，指出了Akbar等一些學者所應用的動量/能量守恒方程中所存在的錯誤，并提出了一個完整的、物理上精確的三方程模型用于模擬地熱井中的不穩(wěn)定流動。此外，石油行業(yè)內井筒內的多相流動與地熱井內的流動具有幾何相似性和物理模型的相似性，目前在石油行業(yè)的油氣井多相流模擬中，已經發(fā)展出瞬態(tài)的五方程模型和六方程模型[158-159]，因此對比石油行業(yè)中的前沿研究成果，發(fā)展用于求解地熱井筒內復雜流動的數(shù)學模型是值得關注的研究方向。

5 干熱巖地熱發(fā)電技術

5.1 地熱發(fā)電基本原理

通常來說，地熱資源的利用方式取決于地熱流體溫度的高低，低溫的地熱資源（小于90 ℃）一般是對其進行直接利用，例如種植業(yè)、溫泉等；對于中溫的地熱資源（90～150 ℃），分布式的制冷和供熱是重要的利用方式；對于溫度較高的地熱資源（大于150 ℃），地熱發(fā)電是實現(xiàn)能量利用的有效途徑[160-161]。

地熱發(fā)電原理是基于朗肯循環(huán)，利用高溫蒸汽驅動膨脹機作功，實現(xiàn)熱電轉換。依據(jù)地熱資源特性，地熱發(fā)電廠的循環(huán)利用方式分為直接發(fā)電（干蒸汽發(fā)電系統(tǒng)、閃蒸發(fā)電系統(tǒng)）和間接發(fā)電（雙工質發(fā)電系統(tǒng)）[162]（圖14）。

針對高溫水熱型、干熱巖型地熱資源中含氣率較高的地熱流體的利用，一般通過降壓閃蒸的方式獲取高溫蒸汽，直接驅動膨脹機作功，即干蒸汽發(fā)電（圖14-a）、單級閃蒸（圖14-b）、雙級閃蒸（圖14-c）。

干蒸汽發(fā)電系統(tǒng)結構簡單、操作維修方便，從地下開采的高焓地熱蒸汽經過必要的凈化和脫酸處理后即可驅動膨脹機作功。但是，由于干蒸汽發(fā)電對地熱條件要求較高，只能在少數(shù)地方使用而未能得到推廣。相較于干蒸汽發(fā)電，閃蒸發(fā)電則具有更強的靈活性，通過降壓閃蒸的方式，將地熱水轉換為飽和蒸汽，進而實現(xiàn)熱電轉換。閃蒸發(fā)電多適用于170 ℃以上的中、高溫地熱流體，目前在全球地熱裝機容量發(fā)電占據(jù)半數(shù)以上[163-164]。我國廣東豐順的鄧屋地熱發(fā)電廠即采用閃蒸發(fā)電模式，裝機容量為300 kW[165]。此外，為了提高閃蒸循環(huán)的效率，在單閃蒸循環(huán)的基礎上進行了雙級閃蒸循環(huán)改進，即對一次閃蒸液相部分進行第二次閃蒸，獲取更高的凈輸出效益，例如冰島東北部的Krafla地熱田[166]。針對閃蒸發(fā)電模式，Jalilinasrabady等[167]對伊朗的Sabalan地熱發(fā)電廠單級閃蒸和雙級閃蒸循環(huán)系統(tǒng)進行熱力學和經濟性的比較，研究表明雙級閃蒸循環(huán)系統(tǒng)的凈輸出功率高于單閃循環(huán)系統(tǒng)，但雙級閃蒸模式的投資較高，回收期較長，綜合考慮熱經濟效益，雙級閃蒸發(fā)電模式具有更加廣闊的應用前景。

閃蒸發(fā)電的循環(huán)方式有結構簡單，熱效率較高的優(yōu)點，但是存在一定的局限性：當?shù)責崃黧w溫度低于130 ℃時，受到水的物性條件限制，閃蒸模式無法獲取滿足地熱發(fā)電所需的高焓蒸汽，采用引入中間介質的雙工質發(fā)電模式（有機朗肯循環(huán)[168]、卡琳娜循環(huán)[169]）是實現(xiàn)對中低溫地熱資源利用的有效途徑。由于有機朗肯循環(huán)（Organic Rankine Cycle，ORC）屬于閉口循環(huán)，純凈的低沸點工質通過換熱器吸收地熱能量，避免了地熱流體對管路的腐蝕，因此ORC在中低溫發(fā)電模式中展現(xiàn)出明顯的優(yōu)越性。在ORC循環(huán)中，根據(jù)地熱流體的特性引入低沸點的工質，使其能夠在較低地熱溫度條件下蒸發(fā)，進入膨脹機推動作功，乏汽進入冷凝器被冷凝，再由泵升壓進入蒸發(fā)器，從而完成如圖14-d所示的循環(huán)過程。隨著深層地熱資源的開發(fā)，EGS技術與ORC相結合的EGS-ORC技術逐漸展現(xiàn)優(yōu)勢。2009年法國蘇爾士地熱田率先建成了1.5 MW地熱電站，采用異丁烷作為工質驅動膨脹機作功，并成功實現(xiàn)了并網運行[170]。與閃蒸發(fā)電模式相似，ORC循環(huán)衍生出的雙級ORC循環(huán)發(fā)電具有更高的作功能力[171-172]，雙級ORC通過對地熱能的分段分品利用，有效降低了循環(huán)不可逆損失，提升了系統(tǒng)性能。在ORC循環(huán)中工質的選擇非常關鍵，基于不同類型的地熱資源，選擇不同沸點的工質與熱源進行匹配，能夠有效提高地熱能的利用率[173]。此外，基于閃蒸和雙工質兩種循環(huán)方式建立的閃蒸—雙工質發(fā)電廠[174]，能夠充分發(fā)揮兩種循環(huán)模式的優(yōu)勢，實現(xiàn)對地熱的溫度分區(qū)利用：高溫區(qū)間的地熱流體用于閃蒸發(fā)電，閃蒸后的中溫地熱流體為ORC循環(huán)供熱，從而完成對地熱能的梯級利用，提高地熱發(fā)電的效率。

根據(jù)地熱資源特性選擇合理的循環(huán)發(fā)電模式是實現(xiàn)提高熱電轉換效率的關鍵：干蒸汽、閃蒸等直接發(fā)電模式多適用于高溫地熱資源，ORC等間接發(fā)電模式則是中低溫地熱資源利用的重要方式。然而，目前針對各種發(fā)電技術的具體適用溫度、壓力及地質條件尚未明晰，僅以溫度為技術選擇標準存在局限性，且地熱采出溫度明顯低于傳統(tǒng)火力電廠，這將導致地熱電效率低下。因此，改善地熱發(fā)電技術、探明地熱發(fā)電技術與地熱資源條件的適配性是未來地熱發(fā)電技術的發(fā)展趨勢。

5.2 地熱發(fā)電市場應用

隨著石化行業(yè)鉆井、壓裂等技術的突破，對高溫深層地熱資源例如高溫水熱型、干熱巖型地熱資源的開發(fā)促進了地熱發(fā)電的發(fā)展。1904年世界上第一座地熱發(fā)電廠在意大利投建并成功運行[175]，地熱發(fā)電的研究自此蓬勃展開，并因其清潔可持續(xù)的特點在世界范圍廣泛應用。除了意大利之外，美國、菲律賓、墨西哥、意大利、冰島等地熱資源豐富的國家的地熱發(fā)電裝機容量近年來迅速增加[163,176-177]，表4、圖15所示為全球典型地熱發(fā)電市場應用。截至2020年底，全球地熱裝機容量已達15 608 MW，其中，美國地熱發(fā)電裝機3 714 MW，占據(jù)世界首位[178]。在我國，高溫的地熱資源集中分布在西南地區(qū)，基于此條件，我國先后在西藏建立了羊八井、羊易、那曲、郎久等地熱發(fā)電站，截至2021年西藏地區(qū)地熱電站裝機達42.18 MW，有效緩解了市區(qū)供電壓力。80 ℃以上中低溫地熱資源的地熱電站則是在廣東豐順、河北后郝窯等地展開，由于資源和技術限制，現(xiàn)僅剩廣東豐順的300 kW機組在間斷運行[179-181]。

表4 全球典型地熱廠統(tǒng)計表[182-186]

圖15 典型地熱發(fā)電廠現(xiàn)場圖

針對干熱巖的地熱發(fā)電利用，目前尚未實現(xiàn)大規(guī)模的商業(yè)應用。20世紀70年代，美國洛斯阿拉莫斯實驗室在新墨西哥州Fenton Hill地區(qū)[191]實施了全球首例干熱巖項目，自此之后英國、日本、德國、澳大利亞、法國等國也相繼開展了干熱巖的項目建設[192]。值得關注的是，法國的Soultz[193]、德國的Landau[194]、澳大利亞的Habanero[195]、美國的Desert Peak[196]和The Geysers[197]等試驗項目目前已初步具備1～5 MW的發(fā)電能力，實現(xiàn)了商業(yè)運營。上述案例也將為后續(xù)干熱巖項目的規(guī)?；l(fā)展提供寶貴經驗，表5對具有代表性的干熱巖發(fā)電項目進行了總結。我國在干熱巖開發(fā)利用方面仍處于資源勘探及EGS技術研究的初步發(fā)展階段，目前已在青海、西藏、福建等地進行了地熱資源勘探工作，其中在青海共和盆地發(fā)現(xiàn)高品質干熱巖體，GR1井底最高溫度236 ℃，孔內3 366 m以下深度平均地溫梯度可達8.8 ℃/100 m。此外，2012年我國啟動了“863”計劃項目“干熱巖熱能開發(fā)與綜合利用關鍵技術研究”，相關研究成果為我國干熱巖資源開發(fā)利用提供了理論依據(jù)和關鍵技術支撐。但距離我國真正實現(xiàn)干熱巖地熱資源的發(fā)電利用，仍需進一步地解決本文所述EGS開發(fā)過程中的一系列關鍵技術問題。

6 結論與展望

干熱巖開采中存在人工壓裂縫網不合理、多尺度多場耦合規(guī)律不明、井筒流體閃蒸造成采熱效率低下、地熱資源熱電轉換效率低等問題，需要解決4個關鍵技術問題。

1）干熱巖儲層人工壓裂技術：干熱巖儲層較為致密，滲透率極低，一般需要通過人工壓裂的手段為流體創(chuàng)造流動換熱通道，但不合理的壓裂縫網易引發(fā)流體短路，造成過早的熱突破，制約地熱系統(tǒng)的可持續(xù)開發(fā)利用。應用暫堵轉向壓裂、液氮壓裂等技術可幫助形成復雜的人工縫網，但仍需明晰縫網的形成機制，開發(fā)精確的裂縫擴展預測模型來指導壓裂工作。對于干熱巖水力壓裂模型方向，將向能考慮干熱巖熱彈塑性變形與破壞、壓裂液與干熱巖相互作用、應力腐蝕引起的亞臨界裂紋擴展等現(xiàn)象的全三維壓裂模型發(fā)展；同時對于干熱巖水力壓裂算法，由于全三維壓裂模型計算成本較高，水力壓裂算法向最佳正交分解（POD）、廣義最佳分解（PGD）、減基（Reduced Basis）等降階算法發(fā)展，從而提高水力壓裂數(shù)值模擬的運算效率。

2）干熱巖開采數(shù)值模擬技術：干熱巖儲層具有多尺度、強非均質特征，開采過程受到熱—流—固—化多場耦合作用的影響，屬于典型的多尺度多場耦合問題，但相關耦合作用機理尚不明確。對于干熱巖巖心尺度多場耦合模型，將向更精細刻畫孔隙結構、更精確計算孔隙內部流動換熱、更精確計算固體骨架局部變形及斷裂擴展的方向發(fā)展，離散力學模型與新型流體力學求解方法如LBM、SPH將更深入的結合和應用，此外巖心尺度THMC耦合數(shù)值模型待需開發(fā)。對于儲層尺度多場耦合模型，精確描述裂縫分布和走向的離散裂縫類模型是開發(fā)熱點，各場控制方程之間會引入更全面的關聯(lián)項和耦合方程進而符合實際物理過程，對于方程求解將大力應用XFEM、XFVM等先進方法以實現(xiàn)高效模擬。另一方面，尺度升級方法可以綜合考慮儲層中裂縫和孔隙多尺度特征，將與多場耦合模型結合，實現(xiàn)宏觀和微觀尺度耦合。

3）井筒熱流體高效提取技術：地熱井筒由于高井深、大壓降的特點，內部流體在舉升過程中易發(fā)生閃蒸，極大地影響井內的流動換熱特性，引發(fā)流動震蕩、井內碳酸鈣結垢等問題，威脅地熱系統(tǒng)的安全性和可持續(xù)開采。通過實驗可以實現(xiàn)對地熱井內流體閃蒸過程的直觀描述，目前相關研究已在核反應堆設計的自然循環(huán)方面廣泛展開，未來可通過建立合理的模化準則開展地熱井筒內流體閃蒸的等效實驗。對于地熱井筒內的數(shù)值模擬，穩(wěn)態(tài)、單相的模型已受到廣泛應用，瞬態(tài)、兩相模型可更準確地描述井內的閃蒸現(xiàn)象，已初步應用于地熱井筒領域。此外油氣開發(fā)領域具有更加成熟和前沿的井內兩相流動模型，可擴展應用于地熱井筒。同時，將地熱井筒全尺度的數(shù)值模擬與局部區(qū)域的可視化實驗研究相結合，可更好地實現(xiàn)對地熱井內流體閃蒸機理的深入探究，為預防井內流體的閃蒸以及提高井筒熱流體提取效率提供理論支持。

4）干熱巖地熱發(fā)電技術：地熱發(fā)電的主要方式包括直接發(fā)電模式和間接發(fā)電模式，適配高溫和中低溫的地熱資源條件。近年來，地熱發(fā)電技術以其清潔可持續(xù)的特點在全球范圍內得到廣泛應用。然而，目前在地熱發(fā)電技術的選擇存在較大局限性，地熱壓力、流體成分、地質條件等也應作為技術選擇的重要參考。針對地熱資源利用效率低的問題，多能互補和能量梯級利用的分布式能源利用將是提高能源利用率的有效途徑。此外，對于地熱發(fā)電的評估不應僅限于電廠的投建運行，以全生命周期的角度將地熱勘探、開采、利用等環(huán)節(jié)的投入與電力輸出目標效益結合，明晰系統(tǒng)運行周期內的經濟收益、“碳”收益，為地熱發(fā)電的產業(yè)化發(fā)展提供指導。

5）EGS是一個技術密集型系統(tǒng)，對其工程開發(fā)需先從機理分析和可行性上著手。然而，由于地下儲層工況的復雜性和地面設備的不穩(wěn)定性，機理研究往往與實際脫節(jié)，需與先導試驗項目密切結合、互相指導，開展產—研結合模式，在不斷往復中提高認知、突破關鍵點，提升EGS的應用價值，為我國實現(xiàn)“雙碳”目標助力。