趙志宏,劉桂宏,王佳鋮,徐浩然

(清華大學(xué) 土木工程系,北京 100084)

如期實現(xiàn)碳達(dá)峰、碳中和,是推動我國經(jīng)濟高質(zhì)量發(fā)展和生態(tài)環(huán)境高水平保護(hù)的內(nèi)在要求,而建設(shè)清潔低碳安全高效的能源體系,是實現(xiàn)碳達(dá)峰、碳中和的必由之路?！吨腥A人民共和國國民經(jīng)濟和社會發(fā)展第十四個五年規(guī)劃和2035年遠(yuǎn)景目標(biāo)綱要》指出要因地制宜開發(fā)利用地?zé)崮?。地?zé)崮苁且环N清潔低碳、分布廣泛、資源豐富、安全穩(wěn)定的優(yōu)質(zhì)可再生能源[1-4]。根據(jù)埋藏深度可將地?zé)崮芊譃?類:200 m以淺為淺層地?zé)崮?200～3 000 m為中深層地?zé)崮?3 000 m以下為深層地?zé)崮?。目?已實現(xiàn)規(guī)模化開發(fā)利用的水熱型地?zé)豳Y源在中深層和深層熱儲中均有廣泛分布。主要針對中深層和深層水熱型地?zé)豳Y源開展研究,為描述簡便,以下統(tǒng)稱為深層地?zé)崮堋?/p>

地?zé)崮荛_采過程涉及復(fù)雜的溫度場、滲流場、應(yīng)力場、化學(xué)場耦合效應(yīng)[5-7]。溫度場與滲流場的耦合是地?zé)崮荛_采的核心過程,主要是回灌的尾水在熱儲層中滲流換熱重新恢復(fù)高溫。當(dāng)需要對熱儲層進(jìn)行增產(chǎn)穩(wěn)產(chǎn)改造時,或當(dāng)考慮儲層變形、回灌堵塞等問題時,則需要在考慮溫度場與滲流場耦合的基礎(chǔ)上,進(jìn)一步考慮應(yīng)力場、化學(xué)場的耦合作用。數(shù)值模擬已成為研究評價深層地?zé)崮荛_采多場耦合效應(yīng)的重要方法,比如TOUGH-FLAC、OpenGeoSys、COMSOL等多場耦合模擬軟件已廣泛應(yīng)用于深層資源評價與開采方案優(yōu)化設(shè)計[8-26]。雖然多場耦合數(shù)值模型可以精細(xì)化地預(yù)測熱儲溫度場、滲流場、應(yīng)力場、化學(xué)場的動態(tài)響應(yīng)過程,但計算效率通常不高,特別是對于復(fù)雜的井-儲模型。為此,將三維地?zé)峋喕癁橐痪S線單元,實現(xiàn)了地?zé)崽锍叨葟?fù)雜群井系統(tǒng)多場耦合效應(yīng)的高效模擬[7,10,23,27-28]。盡管考慮多場耦合效應(yīng)的深層熱儲數(shù)值模型與算法已日臻成熟,但在深層地?zé)崮荛_采過程中的多場耦合災(zāi)變效應(yīng)仍未受到足夠的重視,尚未形成統(tǒng)一的標(biāo)準(zhǔn)體系對地?zé)衢_采多場耦合災(zāi)變效應(yīng)進(jìn)行評價與調(diào)控。

本文總結(jié)筆者團隊近年來針對水熱型地?zé)豳Y源開發(fā)利用多場耦合效應(yīng)開展的數(shù)值模擬研究成果,提出地?zé)崮芸沙掷m(xù)開采的多指標(biāo)評價體系,充分發(fā)揮地?zé)峋?儲系統(tǒng)多場耦合數(shù)值模型的優(yōu)勢,建立面向城市深層地?zé)崮芸沙掷m(xù)開采的資源量評價與采灌方案優(yōu)化設(shè)計方法,旨在促進(jìn)城市深層地?zé)崮芨笠?guī)模和更高質(zhì)量的開發(fā)利用。

1 城市深層地?zé)崮芸沙掷m(xù)開采的多指標(biāo)評價體系

國家能源行業(yè)標(biāo)準(zhǔn)《地?zé)崮苄g(shù)語》(NB/T 10097—2018)中,有3個術(shù)語的定義均提到了“可持續(xù)”的標(biāo)準(zhǔn)[29]:① “可開采量”,在地?zé)崽锟辈?、開采和監(jiān)測的基礎(chǔ)上,考慮到可持續(xù)開發(fā),經(jīng)擬合計算允許每年合理開采的地?zé)崃黧w量和熱量;② “動態(tài)監(jiān)測”,地?zé)豳Y源在勘探、開采及停采階段,連續(xù)記錄水位、井口溫度、井口壓力、開采量、回灌量和蒸汽比例等,并定時分析地?zé)崃黧w化學(xué)組分和同位素的過程,基于此來判斷熱儲溫度、壓力、流體化學(xué)組份含量及資源量的動態(tài)變化,為地?zé)豳Y源的可持續(xù)利用與管理提供依據(jù);③ “地?zé)崮軆?yōu)化開采”,采用最優(yōu)化的方法開采熱儲,在可持續(xù)且不會帶來環(huán)境危害的基礎(chǔ)上獲得最大采熱量?？梢?可持續(xù)已成為地?zé)豳Y源勘探、評價、開發(fā)的重要前提,然而對于可持續(xù)的評價指標(biāo)及其標(biāo)準(zhǔn)仍未形成共識。以下從溫度、水位、變形3個方面給出深層地?zé)崮芸沙掷m(xù)開采的評價標(biāo)準(zhǔn):

(1)熱突破時間。開采井井口溫度下降不超過允許值,本文暫定溫降2 ℃對應(yīng)的時間為熱突破時間,本文暫定最短熱突破時間為100 a。

(2)水位降深。開采井水位降深不超過允許值,本文暫定允許水位降深為-150 m。

(3)地表變形。長期采灌引起的地面變形不超過允許值,本文暫定地表變形允許值為20 mm。

開采井需同時滿足上述3個環(huán)境約束條件才可運行,否則應(yīng)予以關(guān)停,即地?zé)峋畨勖?tlife)由上述3個約束條件共同決定。因此,可采地?zé)崮蹺prod的計算式為

Eprod=qprotlifeρfCp,f(Tpro-Tin)

(1)

其中,qpro為開采流量;ρf為流體密度;Cp,f為流體比熱容;Tpro為開采溫度;Tin為回灌溫度。前人的研究中,tlife在30～100 a間變化,本文按tlife= 100 a考慮,是對深層地?zé)崮芸沙掷m(xù)開采較為嚴(yán)格的限制。因此,需因地制宜選擇地?zé)衢_采模式,并對采灌井布局、井深、采灌量、回灌溫度以及回灌方式等參數(shù)進(jìn)行優(yōu)化設(shè)計,在可持續(xù)開采條件下最大限度地發(fā)揮地?zé)峋淖匀划a(chǎn)能。

需要特別說明的是,深層地?zé)崮荛_發(fā)利用與地表沉降的相關(guān)性研究仍在進(jìn)行當(dāng)中,本文將地表變形作為限制性條件之一,而且允許值極為嚴(yán)苛,主要目的是完全杜絕深層地?zé)崮荛_采對地表變形的影響。根據(jù)干涉合成孔徑雷達(dá)(InSAR)數(shù)據(jù),雄安新區(qū)在2015—2019年的地表沉降速率為37.9～104.8 mm/a,有學(xué)者推斷地表沉降的原因主要是地?zé)豳Y源的開發(fā)利用[30]?；诙嗫捉橘|(zhì)有效應(yīng)力原理,有學(xué)者認(rèn)為雄安新區(qū)地表沉降與第四系地下水的超采有關(guān),其中地下水超采是地面沉降的主要因素,砂巖熱儲層的無序且沒有采灌均衡的開發(fā)成為部分地區(qū)誘發(fā)沉降的次要因素[31]。根據(jù)水準(zhǔn)測量數(shù)據(jù),2015—2019年山東德城區(qū)500 m以淺地層的地表沉降速率為37.9 mm/a,500 m以深地層的地表沉降速率為11.1 mm/a,其中500～800 m層段土體壓縮量約為8.2 mm/a,800 m以深地層壓縮量約為2.9 mm/a。因此,在多場耦合作用下,深層地?zé)衢_采確有可能導(dǎo)致儲層變形,但對地表沉降的影響有限[7],后續(xù)應(yīng)結(jié)合實測數(shù)據(jù),并建立數(shù)值模型精細(xì)評估深層地?zé)衢_采對地表沉降的影響規(guī)律。

2 城市深層地?zé)峋?儲系統(tǒng)多場耦合模擬方法

深層地?zé)嵯到y(tǒng)主要包括兩部分,一是蓋層、儲層、斷層等地質(zhì)結(jié)構(gòu),二是回灌井、開采井等工程結(jié)構(gòu)。地面換熱系統(tǒng)不在本文考慮范圍之內(nèi),相關(guān)內(nèi)容見文獻(xiàn)[31]。以下簡要介紹深層地?zé)峋?儲系統(tǒng)多場耦合模型及非均質(zhì)儲層構(gòu)造方法[7,26]。

2.1 井-儲系統(tǒng)多場耦合數(shù)值模型

圖1給出了井-儲系統(tǒng)多場耦合模型原理。將熱儲層按飽和孔隙介質(zhì)考慮,采用多孔介質(zhì)熱彈性理論描述熱儲層應(yīng)力、滲流、傳熱三場耦合過程。蓋層與儲層類似,但不考慮其中的滲流過程。將斷層按飽和的光滑平板模型來考慮,分別采用簡化的線性模型描述其法向、切向變形過程,同時考慮剪脹對裂隙變形的影響;采用立方定律描述裂隙滲流過程;傳熱過程包括熱對流、熱傳導(dǎo)、熱交換等。為了提高計算效率而又不降低精度,引入一維線單元對三維井筒結(jié)構(gòu)進(jìn)行簡化,考慮沿井筒軸向的對流傳熱過程,井筒內(nèi)流體與圍巖的熱交換通過等效換熱系數(shù)來近似考慮。具體的控制方程可見文獻(xiàn)[7,26],本文不再贅述。熱儲層和熱蓋層采用四面體單元進(jìn)行剖分,地?zé)峋捎镁€單元進(jìn)行剖分,為兼顧計算精度與效率,地?zé)峋車爸饕獰醿拥木W(wǎng)格尺寸應(yīng)相對細(xì)化,而在遠(yuǎn)離地?zé)峋吧w層區(qū)域則可使用較粗化網(wǎng)格,有助于減少網(wǎng)格數(shù)量,節(jié)省計算時間。建議地?zé)峋車木W(wǎng)格尺寸不大于井徑的10倍。

圖1 井-儲系統(tǒng)多場耦合數(shù)值模型原理[26]Fig.1 Schematic diagram of coupled multi-field numerical model of well-reservoir system[26]

2.2 非均質(zhì)儲層構(gòu)建方法

熱儲層的滲透率、比熱容等參數(shù)具有較強的非均質(zhì)性,其空間非均質(zhì)性不僅決定了流體的流動過程,還影響了傳熱和變形過程[21,24,32]。轉(zhuǎn)向帶法[33-34]、傅里葉積分法[35]等各類地質(zhì)統(tǒng)計學(xué)方法已被應(yīng)用于生成非均質(zhì)熱儲模型。轉(zhuǎn)向帶法的主要特點是將高維的隨機過程簡化為一維情況,而傅里葉積分法可基于多種協(xié)方差函數(shù)構(gòu)建非均質(zhì)模型。無論采用何種方法生成隨機模型,其結(jié)果均具有較好的一致性。決定熱儲參數(shù)非均勻分布特征的控制參數(shù)主要有平均值、方差、相關(guān)長度等。生成非均質(zhì)熱儲模型的通常做法是(圖2中,σ為標(biāo)準(zhǔn)差;x為隨機變量;μ為數(shù)學(xué)期望):① 選定合適的熱儲參數(shù)非均勻分布形式(正態(tài)分布、對數(shù)正態(tài)分布等)及其控制參數(shù)(平均值、方差、相關(guān)長度等);② 采用地質(zhì)統(tǒng)計學(xué)方法生成隨機參數(shù)場;③ 將熱儲層離散為一定數(shù)量的代表性體積單元,賦予隨機參數(shù)場;④ 基于熱儲參數(shù)間的物理關(guān)系,生成其他熱儲參數(shù)的隨機分布場。

圖2 非均質(zhì)熱儲模型生成過程示意Fig.2 Schematic diagram of generating heterogeneous geothermal reservoir models

3 基于數(shù)值模擬的可采資源量評價

3.1 可采資源量評價方法

國家標(biāo)準(zhǔn)《地?zé)豳Y源地質(zhì)勘查規(guī)范》(GB/T 11615—2010)[36]中給出了地?zé)釤崃髁糠ā醿Ψ?、解析模型法、統(tǒng)計分析法、數(shù)值模型法、比擬法等地?zé)豳Y源/儲量計算方法。在地?zé)崽锏目辈槌潭缺容^高,并且具有一定時期的開采歷史和比較齊全的監(jiān)測資料時,尤其是城市地區(qū)的地?zé)崽?應(yīng)建立地?zé)崽锏臄?shù)值模擬模型,用以計算/評價地?zé)醿α?并作為地?zé)崽锕芾淼墓ぞ?。因?筆者基于深層地?zé)峋?儲系統(tǒng)多場耦合模擬方法,計算滿足可持續(xù)開采標(biāo)準(zhǔn)的地?zé)豳Y源可采量,具體流程如下:

(1)地質(zhì)模型構(gòu)建。根據(jù)地?zé)崽镞吔缗c評價區(qū)域確定研究區(qū)范圍,逐層導(dǎo)入各地層頂?shù)装濉鄬用鎺缀螖?shù)據(jù),構(gòu)建所有地層界面、斷層面,并通過地層界面、斷層面、模型實體間的相互分割,刪除多余部分,形成研究區(qū)三維地質(zhì)模型。

(2)井-儲模型構(gòu)建。在地質(zhì)模型的基礎(chǔ)上,按照地?zé)峋鴺?biāo)插入一維線單元代表地?zé)峋?線單元的長度即井深。根據(jù)研究區(qū)斷層帶位置和熱儲層的富水性差異,對熱儲參數(shù)進(jìn)行分區(qū),并賦予合適的模型參數(shù)。設(shè)置合適的溫度場、滲流場、地應(yīng)力場初始條件,并根據(jù)研究區(qū)熱儲層地?zé)岬刭|(zhì)構(gòu)造特征、地?zé)嵝纬蓹C理、水文地質(zhì)條件,確定整個模型熱儲系統(tǒng)的熱水力邊界條件。根據(jù)采灌量、回灌溫度等數(shù)據(jù),為每一眼地?zé)峋O(shè)置溫度和流體的流入流出邊界條件。對于數(shù)值模型的校檢,可以對地?zé)崽镞^去某時期的采灌活動進(jìn)行模擬計算,將模擬結(jié)果與現(xiàn)場監(jiān)測數(shù)據(jù)相對比,以此檢驗?zāi)Ｐ徒⒌暮侠硇?。主要通過校正熱儲層參數(shù)使模擬結(jié)果與監(jiān)測數(shù)據(jù)相吻合,并最終確定模型參數(shù)。建模的滲流過程中,可調(diào)整滲透率、孔隙率和儲水系數(shù)使水位模擬值與水位監(jiān)測數(shù)據(jù)相吻合。而在傳熱過程中,可調(diào)整地溫梯度、比熱容及導(dǎo)熱系數(shù)等,使溫度模擬值與測溫數(shù)據(jù)相吻合。

(3)開采方案優(yōu)選。在采灌開始前需先進(jìn)行地應(yīng)力平衡,計算模型在初始熱水力多場耦合作用下的平衡狀態(tài),然后再運行采灌條件下的數(shù)值模型,根據(jù)生產(chǎn)需求確定計算時長,考慮到采灌周期以月為單位,時步設(shè)置應(yīng)不大于1個月。為數(shù)值模型設(shè)置時間周期函數(shù),將每年離散成供暖季和非供暖季2個時間段。計算結(jié)果以開采井水位、溫度和沉降量變化為主。

借鑒地?zé)衢_發(fā)誘發(fā)地震風(fēng)險管控的“紅綠燈系統(tǒng)”[37],針對不同開采方案下開采井的可持續(xù)開采風(fēng)險進(jìn)行分級評價(表1)。

表1 地?zé)峋沙掷m(xù)開采風(fēng)險等級評價標(biāo)準(zhǔn)

3.2 工程案例

以雄安新區(qū)容東安置區(qū)(面積約13 km2)為案例,計算該區(qū)域的地?zé)豳Y源可采量。從容城地?zé)崽锏刭|(zhì)模型中截取邊長為10 km、深為5 km的區(qū)域建立井-儲數(shù)值模型(圖3),具體參數(shù)取值見表2。模型內(nèi)現(xiàn)有12口開采井和10口回灌井,主要位于原容城縣城,另在容東安置區(qū)規(guī)劃建設(shè)20口地?zé)峋?。將該?儲模型共剖分為四面體單元約26萬個,一維線單元2 000個。根據(jù)水位監(jiān)測結(jié)果,以地表作為壓力基準(zhǔn),通過水柱高度換算構(gòu)建初始水壓場。根據(jù)氣象及有關(guān)地?zé)岬刭|(zhì)資料,以恒溫層帶溫度為基準(zhǔn),通過地溫梯度獲得熱儲溫度,并插值得到初始溫度場。聯(lián)合深孔水壓致裂法與壓磁電感法測量熱儲地應(yīng)力,并通過施加邊界條件形成初始地應(yīng)力場。對于滲流過程,側(cè)邊界為恒壓約束,壓力分布與深度呈線性相關(guān);上、下邊界不透水。對于傳熱,側(cè)邊界為開邊界,溫度由地表溫度和插值得到的地溫梯度場根據(jù)深度計算得到,模型的頂部和底部邊界為熱絕緣邊界。對于力學(xué),側(cè)邊界為應(yīng)力約束,大小主應(yīng)力分別為(-0.033 1z+7.529 6)MPa和(-0.020 8z+6.092 3)MPa,其中z為深度;上邊界是自由的,但下邊界是固定的。對于地?zé)峋木?回灌井設(shè)置為回灌流量和回灌溫度,開采井設(shè)置為開采流量和熱流出。對于地?zé)峋阊鄱?回灌井設(shè)置為質(zhì)量流量和熱流出,開采井設(shè)置為質(zhì)量流量和溫度。通過監(jiān)測井水位變化和測溫數(shù)據(jù),對模型進(jìn)行校準(zhǔn)。校準(zhǔn)后的數(shù)值模型能夠較為準(zhǔn)確地反映熱儲中的熱-水-力多場耦合過程,可用于下一步的計算。

圖3 雄安新區(qū)容東安置區(qū)井-儲系統(tǒng)多場耦合數(shù)值模型[26]Fig.3 Coupled multi-field numerical model of well-reservoir system in the Rongdong resettlement area in Xiongan New Area[26]

為評價容東安置區(qū)深層地?zé)崮芸刹闪?擬布置地?zé)峋?0口,井深2 000 m,井距500 m。在供暖季,回灌井和開采井以相同的流量200 m3/h運行,而在非供暖季地?zé)峋贿\行。在采灌均衡條件下,考慮3種布井方案,包括“軌道式”、“棋盤式”和“集中采灌式”(圖4)。

圖4 雄安新區(qū)容東安置區(qū)熱儲初始溫度場、流場與布井方案[26]Fig.4 Initial temperature field,flow field and well layout scheme in the Rongdong resettlement area in Xiongan New Area[26]

考慮到熱儲流體運動方向,開采井布置在上游區(qū)域,防止開采井過早發(fā)生熱突破。

圖5對3種布井方案的水位場、溫度場和位移場進(jìn)行了比較。圖6給出了3種布井方案中10口開采井的開采風(fēng)險等級分布。在“軌道式”布井方案中,開采井0～5號和9號只有熱突破風(fēng)險,開采風(fēng)險等級較低;開采井6～8號同時面臨熱突破和垂直位移風(fēng)險,開采風(fēng)險等級為較高。在“棋盤式”布井方案中,開采井0～2、4、7和9號的開采風(fēng)險等級為低,均滿足地?zé)豳Y源可持續(xù)開發(fā)利用的評價標(biāo)準(zhǔn);開采井3、5、6和8號只有垂直位移風(fēng)險,開采風(fēng)險等級為較低。在“集中采灌式”布井方案中,開采井0、1和4號只有水位風(fēng)險,開采風(fēng)險等級為較低;開采井3、6、8和9號同時面臨熱突破和垂直位移風(fēng)險,開采井7號同時面臨水位和垂直位移風(fēng)險,這5口開采井開采風(fēng)險等級為較高;開采井2和5號完全不滿足地?zé)豳Y源可持續(xù)開發(fā)利用的評價標(biāo)準(zhǔn),開采風(fēng)險等級為高。以上結(jié)果表明:在均質(zhì)模型中,“棋盤式”布井方案開采風(fēng)險最低,其次是“軌道式”布井方案,而“集中采灌式”布井方案的開采風(fēng)險最高。基于地?zé)豳Y源可持續(xù)開發(fā)利用的評價標(biāo)準(zhǔn),可以得到每口開采井的壽命,根據(jù)壽命計算了3種布井方案在100 a內(nèi)的可采地?zé)崮芸偭?如圖7所示。

此外,假設(shè)容東安置區(qū)熱儲的滲透率服從均值為4.32×10-12m2、方差為1的正態(tài)分布,相關(guān)長度為300 m,采用第2.2節(jié)的方法生成非均質(zhì)儲層模型。設(shè)置與均質(zhì)模型相同的初始和邊界條件,對非均質(zhì)熱儲模型“軌道式”、“棋盤式”和“集中采灌式”3種布井方案進(jìn)行對比分析。與均質(zhì)熱儲模型中的近似橢圓形或圓形冷鋒面相比,非均質(zhì)熱儲模型中的冷鋒面由于優(yōu)勢滲流通道而變得不規(guī)則,由回灌引起的垂直位移在非均質(zhì)熱儲模型中不如均質(zhì)熱儲模型顯著(圖5)。與均質(zhì)熱儲模型相比,水位下降到-150 m以下的地?zé)峋當(dāng)?shù)量明顯減少,而熱突破或超過臨界垂直位移的地?zé)峋當(dāng)?shù)量也有減少,這是因為不規(guī)則的流動通道阻礙了熱突破,降低了開采井的熱收縮。因此,部分開采井從由均質(zhì)模型中的高風(fēng)險變?yōu)榉蔷|(zhì)模型中的低風(fēng)險。結(jié)合圖6和圖7可以得出“集中采灌式”是非均質(zhì)熱儲層的最優(yōu)群井開發(fā)模式,在實際工程中應(yīng)充分考慮熱儲非均質(zhì)特性的影響,通過優(yōu)化井位降低開采井風(fēng)險。

圖5 3種布井方案的水位場、溫度場和位移場[26]Fig.5 Water level field,temperature field and displacement field of three well layout schemes[26]

圖6 3種布井方案下10口開采井的風(fēng)險等級對比(彩色圓圈代表開采井,白色圓圈代表回灌井)Fig.6 Comparison of risk rating of 10 production wells under three well layout schemes (Colorful circles are production wells,white circles are injection wells)

圖7 容東安置區(qū)不同采灌方案100 a可采地?zé)崮蹻ig.7 Recoverable geothermal energy of 100 years under different exploitation schemes in the Rongdong resettlement area

4 基于代理模型的地?zé)崮荛_采方案優(yōu)化設(shè)計

深層地?zé)豳Y源開發(fā)利用中,開采參數(shù)的優(yōu)化設(shè)計對提高地?zé)岵墒章示哂兄匾绊?但目前還僅停留在基于敏感性的優(yōu)化設(shè)計。通過優(yōu)化開采參數(shù),在最大限度可持續(xù)開采地?zé)崮艿耐瑫r,避免誘發(fā)地下水位下降、熱突破和地面沉降等環(huán)境災(zāi)害。根據(jù)地?zé)衢_發(fā)方案的需求,基于地?zé)峋?儲系統(tǒng)多場耦合數(shù)值模型,結(jié)合實驗設(shè)計、代理模型以及多響應(yīng)優(yōu)化分析確定地?zé)峥刹蓞^(qū)的優(yōu)化對井開發(fā)模式,以實現(xiàn)最大限度地可持續(xù)開發(fā)利用地?zé)豳Y源,但儲層參數(shù)的不確定性對可采地?zé)崮艿挠绊懖辉诒疚难芯康姆秶鷥?nèi)。

4.1 開采參數(shù)優(yōu)化設(shè)計方法

地?zé)衢_采參數(shù)優(yōu)化工作流程如圖8所示。首先,根據(jù)地?zé)衢_發(fā)方案的需求,收集地?zé)崽锏牡責(zé)岬刭|(zhì)資料,在此基礎(chǔ)上建立三維熱儲地質(zhì)模型,并利用監(jiān)測數(shù)據(jù)對模型進(jìn)行校準(zhǔn)。然后,根據(jù)實驗設(shè)計矩陣,設(shè)定人為控制參數(shù)后運行數(shù)值模型,并進(jìn)行統(tǒng)計分析,得到熱儲響應(yīng)對人為控制參數(shù)的敏感性大小。最后,建立熱儲響應(yīng)和人為控制參數(shù)之間的代理模型,通過數(shù)值模型驗證代理模型的準(zhǔn)確性,并利用多響應(yīng)優(yōu)化方法,得到人為控制參數(shù)的最優(yōu)組合來滿足地?zé)衢_發(fā)方案中在地?zé)衢_采區(qū)最大限度地開采地?zé)崮艿男枨?。具體流程如下:

(1)選定擬優(yōu)化的開采參數(shù)并確定其上下限,指定熱儲響應(yīng)。本文選取采灌量、回灌溫度、采灌井間距、井深和回灌井方位等參數(shù)作為主控因子,根據(jù)地質(zhì)模型和地?zé)峋\行狀態(tài)確定每個因子的上下限,將熱突破時間、水位降深、地表變形和可采地?zé)崮茏鳛闊醿憫?yīng)。

(2)基于Plackett-Burman設(shè)計構(gòu)建熱儲響應(yīng)代理模型。Plackett-Burman是二水平的部分試驗設(shè)計,對每個因子分別取高低兩水平,即上限“+”和下限“-”,通過比較各個因子兩水平之間的差異來確定因子的顯著性,最終能夠快速、有效地完成影響因素的篩選。Plackett-Burman設(shè)計矩陣隨機生成,每列所包含的高低水平數(shù)相等,即在N次試驗中,每個因子的高低水平各出現(xiàn)N/2次。使用第2.1節(jié)中的城市深層地?zé)峋?儲系統(tǒng)多場耦合數(shù)值模型單獨運行每次實驗。采用多元回歸方法擬合熱儲響應(yīng)與人為控制參數(shù)之間的多項式代理模型。

(3)開采參數(shù)優(yōu)化。多響應(yīng)優(yōu)化的目的是解決響應(yīng)之間的沖突,并確定聯(lián)合優(yōu)化單個響應(yīng)或多個響應(yīng)的參數(shù)設(shè)置組合。解決多響應(yīng)之間沖突的一種主流策略是降維,該策略將多響應(yīng)模型降為具有單個合意性函數(shù)的模型,然后將其作為單目標(biāo)優(yōu)化問題進(jìn)行求解。合意性函數(shù)法已被廣泛用于解決各種多響應(yīng)優(yōu)化問題,它將每個響應(yīng)轉(zhuǎn)換為單個合意性函數(shù),然后聚合為復(fù)合合意性,其定義為所有單個合意性函數(shù)的加權(quán)幾何均值。本文通過尋找一組人為控制參數(shù)取值使地?zé)衢_采區(qū)的可采地?zé)崮茏畲蠡?而不超過第2.1節(jié)中定義的水位埋深、熱突破時間和垂直位移的臨界值,按這3個環(huán)境約束條件同等重要考慮,即權(quán)重相同,暫未考慮賦予3個環(huán)境約束條件不同的權(quán)重。首先獲得每個響應(yīng)的單個合意性(d),然后將單個合意性加權(quán)以形成復(fù)合合意性(D)。采用帝國主義競爭算法(ICA)使復(fù)合合意性(D)最大化。

圖8 地?zé)衢_采方案優(yōu)化設(shè)計流程[39]Fig.8 Flow chart of optimization design of geothermal exploitation schemes[39]

4.2 工程案例

以北京城市副中心及周邊區(qū)域作為案例,針對地?zé)釋到y(tǒng)開采參數(shù)進(jìn)行優(yōu)化設(shè)計,實現(xiàn)地?zé)豳Y源可開采量的最大化。從北京通州區(qū)三維地質(zhì)模型[38]中截取邊長為10 km、深度為5 km的區(qū)域建立數(shù)值模型,在提高計算效率的同時且有較好的代表性,具體參數(shù)取值見表3。該區(qū)域現(xiàn)有7口開采井和4口回灌井,井深為2 203～3 002 m,開采溫度為45～51 ℃,采灌量為1 603～2 573 m3/d。所建立的三維井-儲數(shù)值模型如圖9所示,共剖分四面單元約20萬個,三角形網(wǎng)格約3萬個,一維線單元約500個。模型的初始溫度場可根據(jù)研究區(qū)內(nèi)各地層的地溫梯度來確定,初始壓力場根據(jù)研究區(qū)水位分布監(jiān)測結(jié)果,以地表作為壓力基準(zhǔn)面,進(jìn)行換算來確定。初始地應(yīng)力通過水力壓裂法對北京地區(qū)的5個鉆孔進(jìn)行原位應(yīng)力測量得到。對于滲流過程,側(cè)邊界為恒壓約束,壓力分布與深度呈線性相關(guān);上、下邊界不透水。對于傳熱,側(cè)邊界為開邊界,溫度場由插值得到,模型的頂部和底部邊界為熱絕緣邊界。對于力學(xué),側(cè)邊界為應(yīng)力約束,大小主應(yīng)力分別為(-0.032 8z+2.5)MPa和(-0.022 1z+2)MPa,上邊界是自由的,但下邊界是固定的。對于地?zé)峋木?回灌井設(shè)置為回灌流量和回灌溫度,開采井設(shè)置為開采流量和熱流出。對于地?zé)峋阊鄱?回灌井設(shè)置為質(zhì)量流量和熱流出,開采井設(shè)置為質(zhì)量流量和溫度。

擬在副中心東北、正南2個位置各布置一對地?zé)峋?需對采灌方案進(jìn)行優(yōu)化設(shè)計。首先確定擬優(yōu)化的開采參數(shù)并確定其上下限,見表4。采用Plackett-Burman設(shè)計給出敏感性計算方案(表5),需運行模型共12次,其中開采參數(shù)在上、下限水平下各運行6次。在回灌開始前先進(jìn)行地應(yīng)力平衡,計算模型在初始熱-水-力多場耦合作用下的平衡狀態(tài),然后將采灌井按照Plackett-Burman設(shè)計中的參數(shù)水平設(shè)置并運行模型,所有模型的初始和邊界條件都相同。計算時長為100 a,并將每年分為供暖季4個月,非供暖季8個月,時步取1個月,在此期間,連續(xù)監(jiān)測儲層模型中的水位、溫度和位移的演變。使用1臺CPU為i7-4790 K、內(nèi)存為32 GB的普通電腦,每次模擬平均需要12 h完成。

表3 通州地?zé)崽锬Ｐ蛥?shù)[39]

圖9 通州地?zé)崽锶S地質(zhì)模型示意[39]Fig.9 Schematic diagram of three-dimensional geological model of Tongzhou geothermal field[39]

表4 北京城市副中心地?zé)釋到y(tǒng)開采參數(shù)[39]

表5 Plackett-Burman設(shè)計方案[39]

根據(jù)對井系統(tǒng)I的計算結(jié)果(表5),通過多元回歸得到熱儲響應(yīng)與開采參數(shù)之間的多項式代理模型:

ER=8.59+0.001 172Q-0.293 5T+0.004 689D+ 0.000 243d-0.012 74α(R2=0.93)

(2)

(3)

(4)

(5)

式中,ER為可采地?zé)崮?LW為水位;DV為垂直位移;TBT為熱突破時間。

決定系數(shù)R2在0.88～0.99,說明所建立的代理模型可以較好地預(yù)測熱儲響應(yīng)與開采參數(shù)之間的對應(yīng)關(guān)系。圖10給出了對井系統(tǒng)I的最優(yōu)開采參數(shù)。提高流量可增大可采地?zé)崮?、垂直位移和水位埋?但會縮短熱突破時間,故最優(yōu)流量(2 181 m3/d)在取值范圍中間。降低回灌溫度均會使4個熱儲響應(yīng)參數(shù)增大,但與可采地?zé)崮芎痛怪蔽灰葡啾?回灌溫度對熱突破時間和水位埋深的影響不太顯著,故最優(yōu)回灌溫度為下限25 ℃。增大井距可增大可采地?zé)崮?、熱突破時間和水位埋深,但會降低垂直位移,故最優(yōu)井距(677 m)在取值范圍中間。增大井深會略微增大可采地?zé)崮?、垂直位移和水位埋?但會縮短熱突破時間?？紤]到增加井深會增加鉆井成本,故最優(yōu)井深為下限2 357 m。增大回灌井旋轉(zhuǎn)角可減小可采地?zé)崮?、熱突破時間和水位埋深,但會增加垂直位移。不過,與可采地?zé)崮芎蜔嵬黄茣r間相比,增大回灌井旋轉(zhuǎn)角對垂直位移和水位埋深的影響較小,故回灌井最優(yōu)方位應(yīng)選在滲流方向的下游。

將上述最優(yōu)開采參數(shù)作為城市副中心井-儲數(shù)值模型的輸入?yún)?shù),采用COMSOL計算熱儲響應(yīng),進(jìn)一步檢驗代理模型的準(zhǔn)確性。結(jié)果表明,由于采灌不均衡,通州副中心出現(xiàn)了大范圍的抽水漏斗,最大水位埋深超過了200 m,無法滿足可持續(xù)開發(fā)利用的要求。而地?zé)釋到y(tǒng)I的水位、熱突破和垂直位移均滿足可持續(xù)開發(fā)利用的評價標(biāo)準(zhǔn)。代理模型與數(shù)值模型的對比結(jié)果表明,4個熱儲響應(yīng)的最大計算誤差為11.6%,說明代理模型能夠較準(zhǔn)確地表征熱儲響應(yīng)和開采參數(shù)之間的關(guān)系。與Plackett-Burman設(shè)計的12次運行結(jié)果相比,采用優(yōu)化后的開采參數(shù)得到了滿足可持續(xù)開采標(biāo)準(zhǔn)的最大可采地?zé)崮?圖11)。

采用類似的流程,得到了對井系統(tǒng)II的代理模型及其最優(yōu)開采參數(shù)(圖10)。相比于對井系統(tǒng)I,對井系統(tǒng)II的最優(yōu)流量和井深較小,而最優(yōu)井距較大?？赡艿脑蚴菍到y(tǒng)II的平均滲透系數(shù)(0.38 m/d)小于對井系統(tǒng)I(0.5 m/d),需采用較小的流量確保水位埋深不超過150 m。因此,由于儲層參數(shù)的非均質(zhì)性和地下流場的各向異性,同一地?zé)崽镏胁煌責(zé)衢_采區(qū)內(nèi)對井系統(tǒng)的最優(yōu)人為控制參數(shù)可能會有明顯的不同。

圖11 可采地?zé)崮芘c開采風(fēng)險等級[39]Fig.11 Recoverable geothermal energy and production risk rating[39]

5 結(jié) 論

(1)從溫度、水位、變形3個方面給出了深層地?zé)崮芸沙掷m(xù)開采的多指標(biāo)評價體系,以及基于數(shù)值模型的可采地?zé)崮苡嬎愎?。但?各評價指標(biāo)的允許值應(yīng)根據(jù)具體情況合理制定。

(2)基于深層地?zé)峋?儲系統(tǒng)多場耦合數(shù)值模型給出了滿足可持續(xù)開采標(biāo)準(zhǔn)的可采地?zé)崮茉u價方法,提出了開采井風(fēng)險評價的“紅綠燈系統(tǒng)”。

(3)基于深層地?zé)峋?儲系統(tǒng)多場耦合數(shù)值模型給出了滿足可持續(xù)開采標(biāo)準(zhǔn)的采灌方案優(yōu)化設(shè)計方法,使用Plackett-Burman設(shè)計構(gòu)建了代理模型,通過多響應(yīng)優(yōu)化實現(xiàn)了參數(shù)的快速優(yōu)化。

致謝研究過程中得到了中國地質(zhì)科學(xué)院水文地質(zhì)環(huán)境地質(zhì)研究所王貴玲研究員和馬峰研究員級高工、北京市地?zé)嵴{(diào)查研究所張進(jìn)平正高工等專家的指導(dǎo)和幫助,在此衷心感謝!