增強(qiáng)型地?zé)嵯到y(tǒng)的多區(qū)域多物理場耦合三維數(shù)值模擬*

丁軍鋒，王世民

(中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院， 北京 100049; 中國科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室， 北京100049)

地?zé)崮苁莵碓捶€(wěn)定的可再生清潔能源，儲(chǔ)量巨大、無污染。開發(fā)利用地?zé)崮芗饶軡M足人類的能量需求，又能保護(hù)地球環(huán)境免遭破壞，因而對(duì)國民經(jīng)濟(jì)的可持續(xù)發(fā)展意義重大。地?zé)崮荛_發(fā)的重中之重是地?zé)岚l(fā)電,其能源利用系數(shù)遠(yuǎn)高于水力、風(fēng)力、太陽能發(fā)電和地?zé)嶂苯永肹1-3]。由于地?zé)犭娬究刹皇芗竟?jié)、氣候、晝夜等自然條件的影響而不間斷運(yùn)行，地?zé)岚l(fā)電可以作為國家電網(wǎng)的基礎(chǔ)載荷，也易于調(diào)峰和實(shí)施熱電聯(lián)供[4-5]。目前地?zé)崮艿拈_發(fā)主要是中低溫地?zé)崮艿闹苯永?，地?zé)岚l(fā)電沒有得到應(yīng)有的發(fā)展[6]。

1 研究背景

地?zé)豳Y源按其成因和產(chǎn)出條件可分為水熱型和干熱巖型[7]。水熱型地?zé)豳Y源賦存于高滲透性的孔隙或裂隙介質(zhì)中，與年輕火山活動(dòng)或高熱流背景相伴生形成高溫水熱系統(tǒng)，而處于正?；蚱蜔崃鞅尘跋碌牡叵滤h(huán)通常形成中低溫水熱系統(tǒng)，主要用于地?zé)崮艿闹苯永?。干熱巖型地?zé)豳Y源賦存于地下較深處的高溫但低滲透性巖體中，原則上不受地區(qū)分布限制，但其開采需要借助人工壓裂進(jìn)行儲(chǔ)層改造(reservoir stimulation)，進(jìn)而通過注水循環(huán)形成增強(qiáng)型地?zé)嵯到y(tǒng)(enhanced geothermal system，EGS)[8]。

近年來EGS的研究和開發(fā)廣受關(guān)注，在美國、英國、法國、德國、瑞士、日本、澳大利亞等國家已經(jīng)進(jìn)行了EGS工程試驗(yàn)[5-6, 9-11]，但商業(yè)運(yùn)行EGS電站的成功案例還很少[1]。EGS是一個(gè)有廣泛應(yīng)用前景，值得深入研究的問題。EGS發(fā)電需要滿足兩個(gè)主要條件：一個(gè)是要從儲(chǔ)層中提取出足夠高溫度的能量，另一個(gè)是采熱過程中要有足夠的流體流量供給[12]。一般來講，EGS選用地下3～10 km深處較高溫度的干熱巖作為儲(chǔ)層。但因干熱巖通常具有低孔隙度、低滲透率的特點(diǎn)，需要對(duì)儲(chǔ)層進(jìn)行人工壓裂形成流體運(yùn)移通道，以滿足EGS發(fā)電的流體流量要求，并保證流體在儲(chǔ)層內(nèi)的運(yùn)移過程中能夠和圍巖進(jìn)行充分的熱交換。

EGS問題是一個(gè)典型的多區(qū)域、多物理場耦合問題。在由注水井、生產(chǎn)井、孔隙儲(chǔ)層、不可滲透圍巖組成的多區(qū)域EGS中，不同區(qū)域服從不同的控制方程，但在區(qū)域邊界上要滿足物理上合理的連接條件。同時(shí)，在EGS模擬中，溫度場和滲流場必須耦合求解，有時(shí)還涉及與固體變形、組分輸運(yùn)的聯(lián)合求解。國內(nèi)外文獻(xiàn)中發(fā)表的代表性研究工作包括：1)Gong 等[13]對(duì)華北油田伴生地?zé)岚l(fā)電項(xiàng)目的注水-采液過程進(jìn)行滲流和傳熱耦合的數(shù)值模擬，定量研究注水速度、注水溫度對(duì)地?zé)崽餃囟入S時(shí)間變化的影響；2)Bataillé等[14]數(shù)值模擬法國Soultz-sous-Forets EGS地?zé)崽镏械淖匀粚?duì)流和強(qiáng)迫對(duì)流，揭示裂隙巖石中自然對(duì)流對(duì)維持地?zé)崽餃囟冗M(jìn)而延長地?zé)犭娬緣勖淖饔茫?)Bl?cher等[15]研究德國一對(duì)地?zé)嶙⑺蜕a(chǎn)井的長期工作狀態(tài)，預(yù)測這對(duì)地?zé)峋畬⒃谑褂?.6 a后開始發(fā)生熱貫通(thermal breakthrough)；4)Zhou和Hou[16]提出一個(gè)模擬水力壓裂過程的數(shù)值模型，研究裂隙的動(dòng)態(tài)擴(kuò)張過程及孔隙流體與裂隙流體之間的質(zhì)量、動(dòng)量交換；5)Jiang等[17-18]以及Chen和Jiang[19]基于Brinkman方程數(shù)值模擬EGS的滲流傳熱過程，考慮滲透率及圍巖與孔隙流體間熱不平衡對(duì)EGS壽命的影響，Cao 等[20]進(jìn)而研究多種地?zé)峋季旨皹?gòu)造應(yīng)力作用下注水流量與采熱效率、EGS壽命之間的關(guān)系；6)針對(duì)以CO2為工質(zhì)的EGS，Luo等[21]研究井中湍流及其與儲(chǔ)層滲流的耦合，Li等[22]用實(shí)驗(yàn)方法研究穩(wěn)態(tài)自然對(duì)流對(duì)采熱的影響，Jiang等[23]比較不同井型(包括水平井)對(duì)采熱的影響；7) Zeng 等[24-26]數(shù)值模擬Desert Peak、羊八井地?zé)崽锏牟蔁徇^程，包括水平井和垂直裂縫情形;8)Pruess[27-28]從可壓性、能量提取率和流體損失等角度對(duì)比基于CO2和水兩種工質(zhì)的EGS;9)Saeid等[29-31]對(duì)低溫地?zé)嵯到y(tǒng)進(jìn)行一維和兩維耦合模擬，考慮層流和湍流熱導(dǎo)率的不同，通過不同的參數(shù)對(duì)比，得出地?zé)崽飰勖蕾囉诳紫抖?、流量、井距、?chǔ)層溫度、注水溫度的定量關(guān)系；10)Huang等[32]研究雙井EGS井中湍流和儲(chǔ)層達(dá)西流的耦合，發(fā)現(xiàn)井中壓力變化主要由靜水壓力主導(dǎo)，進(jìn)而以一維井內(nèi)流動(dòng)假設(shè)為基礎(chǔ)得到流體溫度的解析解。

在典型EGS工作條件下，井內(nèi)流動(dòng)以湍流主導(dǎo)。例如，在一個(gè)貫穿500 m儲(chǔ)層的半徑為0.15 m的注水井或生產(chǎn)井中，如果質(zhì)量流量為50 kg/s，容易證明，只有在井底以上5 m范圍內(nèi)井內(nèi)流動(dòng)的雷諾數(shù)小于2 200, 即圓形管道內(nèi)層流向湍流過渡的臨界雷諾數(shù)，而在99%的井深范圍內(nèi)井內(nèi)流動(dòng)都處于湍流狀態(tài)。為考慮EGS涉及的井內(nèi)湍流作用，采用湍流模型進(jìn)行數(shù)值模擬是一個(gè)可行的途徑，但需要注意不同湍流模型適用的條件。例如，k-ε模型數(shù)值穩(wěn)定性好，尤其適用于大雷諾數(shù)流動(dòng)，但在固體邊界附近區(qū)域表現(xiàn)較差，而k-ω模型能較好刻畫邊界附近湍流，但對(duì)入口湍流條件較為敏感，導(dǎo)致數(shù)值穩(wěn)定性較差[33-34]。結(jié)合k-ε與k-ω兩種湍流模型的優(yōu)點(diǎn)，Menter[35]提出SST(shear stress transport)模型。此外，還有通過引入阻尼函數(shù)在邊界附近區(qū)域?qū)-ε模型修正而建立的小雷諾數(shù)模型[33]。

前人研究工作中，對(duì)多區(qū)域耦合采用簡化近似的方法，或者直接指定區(qū)域界面的速度或壓力[21]，或者以Brinkman方程同時(shí)描述井內(nèi)區(qū)域和儲(chǔ)層區(qū)域(通過孔隙度和滲透率不同取值體現(xiàn)區(qū)域差異，將井內(nèi)區(qū)域看作滲透率趨近無窮大情形)[17-20]。由于Brinkman方程只對(duì)高孔隙度介質(zhì)成立[36]，將其外推應(yīng)用于低孔隙度EGS儲(chǔ)層并沒有足夠的理論或?qū)嶒?yàn)依據(jù)。事實(shí)上，更簡單的經(jīng)典達(dá)西定律能夠更精確地刻畫低孔隙度介質(zhì)內(nèi)的滲流[36-37]。另一方面，對(duì)于井中的自由流體而言，考慮到地?zé)岚l(fā)電要求的流量至少為250 m3/h[11]，對(duì)應(yīng)于井內(nèi)流動(dòng)由湍流主導(dǎo)，而Brinkman方程根本無法與湍流模型結(jié)合模擬井中湍流。前人雖有基于湍流模型模擬EGS井內(nèi)湍流的研究[21]，但僅限于k-ε模型。由于k-ε模型適用于大雷諾數(shù)流動(dòng)，而EGS井內(nèi)湍流屬于中低雷諾數(shù)情形，基于k-ε模型的結(jié)果需要通過與其他湍流模型結(jié)果對(duì)比加以檢驗(yàn)。

本文基于多區(qū)域多物理場耦合的三維有限元模型，系統(tǒng)研究EGS滲流與傳熱過程及其對(duì)EGS電站壽命的影響。通過連接條件實(shí)現(xiàn)不同區(qū)域間溫度場、壓力場和速度場的自然耦合，結(jié)合多種湍流模型模擬井內(nèi)湍流，并探討以二維模型模擬EGS的條件和有效性。

2 計(jì)算模型

本文模擬的EGS分成3個(gè)區(qū)域，即井內(nèi)區(qū)域、滲流儲(chǔ)層區(qū)域和不可滲透圍巖區(qū)域。假設(shè)儲(chǔ)層處于孔隙水飽和狀態(tài)且沒有自然對(duì)流發(fā)生。模擬中，將井內(nèi)和儲(chǔ)層中的流體流動(dòng)設(shè)為定常的，而整個(gè)EGS的溫度場是隨時(shí)間變化的。

井內(nèi)和儲(chǔ)層中的流體速度均滿足不可壓縮連續(xù)性方程

(1)

不同區(qū)域的流動(dòng)滿足不同形式的動(dòng)量方程。儲(chǔ)層中的孔隙流動(dòng)滿足達(dá)西定律

(2)

式中：k,μ,p分別為儲(chǔ)層滲透率、流體黏度和孔隙壓力。井內(nèi)自由流體的湍流則滿足動(dòng)量方程

(3)

式中：I為單位(恒等)張量，ρ為流體密度，μT為湍流黏度。μT的具體形式由湍流模型定義[33-35,38-41]，而層流對(duì)應(yīng)于μT=0情形。

EGS各區(qū)域內(nèi)滿足統(tǒng)一形式的能量方程

(4)

式中：T為溫度，t為時(shí)間，φ為孔隙度，ρs為固體巖石密度，cf、cs分別為流體和固體比熱。通過給定孔隙度取值，方程(4)可以應(yīng)用于包括井內(nèi)區(qū)域(φ=1)和圍巖區(qū)域(φ=0)在內(nèi)的整個(gè)EGS計(jì)算域。

為實(shí)現(xiàn)區(qū)域耦合，不同區(qū)域的界面上需要滿足以下連接條件：1)在注水井、生產(chǎn)井與儲(chǔ)層的界面上井內(nèi)流體壓力與儲(chǔ)層內(nèi)孔隙壓力連續(xù)，井內(nèi)流體速度與儲(chǔ)層內(nèi)達(dá)西速度連續(xù)；2)在注水井和生產(chǎn)井與儲(chǔ)層的界面上、儲(chǔ)層與圍巖的界面上溫度和熱流均連續(xù)。

本文模擬的EGS為一個(gè)500 m×600 m×500 m的長方體區(qū)域?？紤]到問題關(guān)于y=0的對(duì)稱性，僅計(jì)算y>0一側(cè)即可，如圖1所示。整個(gè)計(jì)算域包括注水井、生產(chǎn)井、儲(chǔ)層與圍巖4個(gè)區(qū)域。采用非結(jié)構(gòu)化四面體網(wǎng)格對(duì)計(jì)算域進(jìn)行剖分，共964 137個(gè)單元，其中井內(nèi)406 075個(gè)單元，井中單元最長邊0.63 m，最短邊0.035 mm，儲(chǔ)層中最長邊50 m，最短邊0.63 m。

圖1 有限元模型及計(jì)算域網(wǎng)格剖分Fig.1 Finite element model and computational mesh

有限元模型參數(shù)與物性參數(shù)分別由表1和表2列出。計(jì)算中，在注水井入口指定定常流量(86.8 kg/s)的邊界條件，而在生產(chǎn)井出口指定壓力為零的邊界條件。在模擬的EGS儲(chǔ)層厚度和地溫梯度下，自然對(duì)流不會(huì)發(fā)生[40]。計(jì)算中，溫度場的瞬態(tài)求解取時(shí)間步長0.05 a。

表1 EGS有限元模型參數(shù)Table 1 Parameters of the finite element model for EGS

表2 流體和巖石熱物理性質(zhì)Table 2 Thermo-physical properties of fluid and rock

3 結(jié)果與討論

基于SST湍流模型計(jì)算得到的三維壓力場和流場由圖2給出。圖2(a)和2(b)表明，壓力和流體速度主要在水平方向上變化，而在垂向上基本保持不變，呈現(xiàn)出明顯的分層流動(dòng)特征。圖2(c)顯示流線跨井壁保持連續(xù)。這一結(jié)果清楚地表明，通過施加連接條件能夠成功實(shí)現(xiàn)流場的區(qū)域耦合。圖2(d)反映出井周圍區(qū)域流體速度大，而在離井較遠(yuǎn)區(qū)域流速則較小。兩井之間的流線較為密集，形成相對(duì)優(yōu)勢流體通道。

圖3給出基于SST湍流模型計(jì)算得到的三維溫度場演化，包括EGS工作1、5、10和20 a共4個(gè)不同時(shí)刻的溫度截圖。圖3顯示，由注水井流入的冷水首先冷卻注水井周邊區(qū)域，而后冷卻區(qū)域逐漸向生產(chǎn)井一側(cè)擴(kuò)展。溫度場在垂向上變化很小，同時(shí)在圍巖區(qū)域熱擴(kuò)散很緩慢，說明傳熱過程以水平方向熱對(duì)流主導(dǎo)，而熱傳導(dǎo)貢獻(xiàn)甚微。圖3所示溫度場分布的特征還表明，傳熱的區(qū)域耦合也通過施加連接條件得以成功實(shí)現(xiàn)。

不同井內(nèi)流動(dòng)三維模型預(yù)測的生產(chǎn)井內(nèi)壓力垂向變化與出口溫度時(shí)間變化由圖4給出。對(duì)比的模型包括k-ω，k-ε，小雷諾數(shù)k-ε，SST共4種湍流模型和層流模型。圖4(a)表明，4種湍流模型給出的壓力變化基本一致，而層流模型預(yù)測的井內(nèi)總壓降則只有湍流總壓降的約1/4。這一結(jié)果清楚地表明井內(nèi)湍流因受到比層流更高的摩擦阻力而產(chǎn)生更大的壓降。然而，不管是不足0.02 MPa的井內(nèi)層流總壓降，還是不足0.08 MPa的井內(nèi)湍流總壓降，都要比注水井與生產(chǎn)井之間有必要說明的是，本文采用的幾種湍流模型均基于雷諾時(shí)間平均假設(shè)，即沒有直接計(jì)算湍流對(duì)應(yīng)的速度隨時(shí)間波動(dòng)，而是依據(jù)特定假設(shè)(不同湍流模型采用不同假設(shè))與速度的時(shí)間平均值相聯(lián)系。基于雷諾平均假設(shè)的湍流模型不能刻畫描述湍流的瞬態(tài)細(xì)節(jié)特征，但用于近似計(jì)算湍流壓降等的時(shí)間平均特征已被大量成功實(shí)例所證明是可靠的，事實(shí)上也已成為主流工程分析數(shù)值模擬軟件的“標(biāo)準(zhǔn)配置”。這些相對(duì)簡單的湍流模型能夠滿足本文研究的需要，然而要想深入了解多區(qū)域、多物理場耦合等情況下的井內(nèi)湍流的細(xì)節(jié)，尤其是湍流狀態(tài)下的傳熱機(jī)理與效率，需要借助更精細(xì)的湍流模擬方法，如大渦旋模擬(large eddy simulation, LES)和直接數(shù)值模擬(direct numerical simulation，DNS)。LES和DNS需要求解瞬態(tài)Navier-Stokes方程，為捕捉湍流渦旋形態(tài)還需要保持時(shí)間步長足夠小，其要求的計(jì)算量很大。在EGS數(shù)值模擬中，采用LES或DNS以更精確地模擬井內(nèi)湍流是一個(gè)潛在的研究方向。

圖2 三維SST模型壓力場與速度場結(jié)果Fig.2 Results of the 3D SST model for pressure and velocity fields

圖3 三維SST模型結(jié)果：EGS溫度場演化Fig.3 Results of the 3D SST model: evolution of EGS temperature field

近50 MPa的井間總壓降(圖2(a))小3個(gè)量級(jí)。因此，盡管湍流與層流的傳熱機(jī)理和效率有本質(zhì)上的不同，但就本文模擬的EGS而言，由于井內(nèi)湍流壓降遠(yuǎn)小于井間滲流壓降，導(dǎo)致井內(nèi)湍流對(duì)EGS采熱過程的影響不大。圖4(b)表明，除k-ω模型結(jié)果略有差異(可能由此模型對(duì)入口湍流條件較為敏感所致)，幾種湍流模型模擬井內(nèi)流動(dòng)給出基本一致的生產(chǎn)井出口溫度演化，且與層流模型結(jié)果很接近。這一結(jié)果說明，至少對(duì)以水為工質(zhì)的EGS，在井內(nèi)湍流效應(yīng)影響不顯著。

圖4 不同井內(nèi)流動(dòng)三維模型預(yù)測結(jié)果對(duì)比Fig.4 Comparison of the predicted results among different 3D models for in-well flows

上述EGS三維數(shù)值模擬結(jié)果表明，井內(nèi)湍流因其壓降遠(yuǎn)小于井間滲流壓降而對(duì)EGS工作狀態(tài)影響很小。這一發(fā)現(xiàn)不但驗(yàn)證了前人在EGS三維數(shù)值模擬中直接將井內(nèi)流動(dòng)簡化為層流的合理性，而且揭示了對(duì)適當(dāng)?shù)腅GS問題采用二維模擬的可能性。我們知道，湍流具有三維流動(dòng)的本質(zhì)屬性，一般情況下是不能簡單地用二維模型有效刻畫和描述的。然而，在已經(jīng)清楚了解井內(nèi)湍流對(duì) EGS 運(yùn)行影響可忽略的前提下，如果EGS的結(jié)構(gòu)和物性都沒有隨深度的變化，而且在儲(chǔ)層中沒有自然對(duì)流發(fā)生，則由于井內(nèi)壓降遠(yuǎn)小于儲(chǔ)層內(nèi)的滲流壓降，EGS的流場將以水平流動(dòng)為特征(圖2)。再考慮到EGS中傳熱過程由熱對(duì)流占主導(dǎo)(圖3)，則本文模擬的EGS問題可以用二維水平流動(dòng)和傳熱模型近似模擬。

從圖1所示三維模型的頂面網(wǎng)格出發(fā)，在注水井壁指定注水流量，而在生產(chǎn)井壁給定零壓力，我們進(jìn)行二維有限元模擬。二維模型得到的壓力、速度、溫度分布(t=20 a)與三維模型結(jié)果高度一致，如圖5所示。

為了更細(xì)致地對(duì)比二維與三維模擬結(jié)果，圖6直接比較二維模型與三維SST模型預(yù)測的生產(chǎn)井出口溫度曲線。兩條曲線總體趨勢一致，但存在明顯的定量差別。為解釋這一差別，我們將二維網(wǎng)格進(jìn)一步加密，即從自由度18 099的2D網(wǎng)格依次增加到自由度357 459的細(xì)網(wǎng)格(fine)和自由度1 415 935的極細(xì)網(wǎng)格(extremely-fine)，并將后兩種網(wǎng)格結(jié)果一并在圖6中比較。很顯然，兩種加密網(wǎng)格給出的二維模擬結(jié)果在圖6中已經(jīng)不可區(qū)分，表明結(jié)果已達(dá)到網(wǎng)格無關(guān)的精確程度，即有限元離散誤差已經(jīng)小到可以忽略不計(jì)的程度，而加密前的二維結(jié)果和三維模擬結(jié)果則包含不容忽略的離散誤差。圖6還清楚地反映出加密前的二維結(jié)果比三維結(jié)果更接近于網(wǎng)格無關(guān)結(jié)果，提示加密前的二維結(jié)果也比三維結(jié)果更精確(因?yàn)槿S結(jié)果還包含在垂直方向上的離散誤差)?？紤]到三維模擬要求的計(jì)算量和計(jì)算時(shí)間遠(yuǎn)大于二維模擬，圖6給出的直接對(duì)比清楚地顯示二維模擬相對(duì)于三維模擬在計(jì)算精度和效率上具有的巨大優(yōu)勢。

圖5 二維模型結(jié)果Fig.5 Results of the 2D model

圖6 三維SST模型預(yù)測的生產(chǎn)井溫度演化與基于3種網(wǎng)格的兩維模型結(jié)果對(duì)比Fig.6 Comparison of production temperature evolution predicted by the 3D SST model with the result predicted by the 2D model employing three meshes

4 結(jié)論

基于有限元模型，對(duì)一個(gè)EGS問題進(jìn)行三維數(shù)值模擬。主要結(jié)論如下：1)通過施加正確的連接條件能夠?qū)崿F(xiàn)EGS在不同區(qū)域之間的多物理場自然耦合；2)4種湍流模型模擬井內(nèi)流動(dòng)給出基本一致的壓力變化，井內(nèi)湍流總壓降約為層流模型的4倍，但比注水井與生產(chǎn)井之間的井間總壓降小3個(gè)量級(jí)；3)井內(nèi)湍流因其壓降遠(yuǎn)小于井間滲流壓降而對(duì)EGS采熱過程總體影響很小，從而在EGS三維數(shù)值模擬中可用簡單層流模型近似描述井內(nèi)流動(dòng)；4)在EGS結(jié)構(gòu)和物性隨深度變化、儲(chǔ)層中自然對(duì)流、井內(nèi)湍流效應(yīng)均可忽略的條件下，EGS中流場以水平方向流動(dòng)為特征，而傳熱過程由水平方向熱對(duì)流占主導(dǎo)，因而可以采用二維模型近似模擬。