收稿日期：2021-11-29

基金項目：國家自然科學(xué)基金（41402228）；天津市教委科研計劃（2017KJ050）；中國能源建設(shè)股份有限公司重大科技專項課題

（CEEC2020-KJZX-10-01）

通信作者：馬玖辰（1980—），男，博士、副教授，主要從事可再生能源利用方面的研究。thermaltju@163.com

DOI：10.19912/j.0254-0096.tynxb.2021-1450 文章編號：0254-0096（2023）05-0030-10

摘 要：基于有限長移動線熱源（MFLS）傳熱模型，根據(jù)時空疊加原理推導(dǎo)出含水層非穩(wěn)態(tài)過余溫度解析解ΔTMFLS；在格子單松弛模型（LBGK）的演化方程中引入離散力源項，建立格子Boltzmann法（LBM）與ΔTMFLS的耦合計算模型。通過熱響應(yīng)實驗，驗證耦合計算模型與求解方法的正確性。研究表明，在不同計算工況下含水層各區(qū)域的滲流速度均具有啟動—下降—回升—穩(wěn)定的4個連續(xù)階段。隨著含水層孔隙率的降低，虛擬流體粒子動能損失增大，滲流速度降幅增大，回升過程緩慢。然而隨著進水流速的提高，孔隙率的變化對于含水層水動力場演化過程的影響程度減弱；含水層熱量運移過程的方向性顯著增強。

關(guān)鍵詞：地?zé)崮?；傳熱；多孔介質(zhì)；數(shù)值模型；含水層；格子Boltzmann方法

中圖分類號：TK529 文獻標(biāo)志碼：A

0 引 言

作為一種可再生能源，地?zé)崮艿拈_發(fā)利用對于建成“近零排放”的能源體系，實現(xiàn)碳中和遠(yuǎn)景目標(biāo)具有深遠(yuǎn)的戰(zhàn)略意義［1］。預(yù)計到2025年，中國地?zé)崮芄┡ㄖ评洌┟娣e將比2020年增加50%［2］；淺層地?zé)崮芸蓪崿F(xiàn)建筑供暖（制冷）面積達320億m2［3］。由于地源熱泵系統(tǒng)采用“取熱不取水”的運行方式，當(dāng)前已成為中國淺層地?zé)崮艿闹饕媚Ｊ剑?］。研究表明，通過改變地下水滲流過程，有助于緩解土壤“冷熱堆積”現(xiàn)象，強化地埋管傳熱過程［5］。

由于含水層孔隙結(jié)構(gòu)的復(fù)雜性、多樣性與隨機性，利用有限元與有限差分法開展數(shù)值計算，然而將滲流溶液作為宏觀連續(xù)介質(zhì)，計算結(jié)果受限于網(wǎng)格劃分方式，且存在運算時間長、收斂困難等問題［6-7］，因此國內(nèi)外學(xué)者嘗試采用格子Boltzmann法（lattice Boltzmann method，LBM）從多孔介質(zhì)的均質(zhì)性、滲透性以及滲流溶液飽和度入手對含水層中滲流過程進行研究［8-10］。文獻［11］提出表征單元體積尺度級聯(lián)格子Boltzmann法，探究孔隙率和比熱容對多孔介質(zhì)內(nèi)混合對流及自然對流的影響；文獻［12］通過LBM與X射線成像技術(shù)，研究了孔隙結(jié)構(gòu)對多孔介質(zhì)堵塞的影響程度；文獻［13］采用雙粒子分布函數(shù)LBM，分析了滲流溶液的流動模式與孔隙結(jié)構(gòu)非均質(zhì)性的相關(guān)性。

綜上，將LBM引入含水層熱-流運移模型，從介觀尺度探究含水層滲流過程與傳熱機制變化的研究鮮見報道。本文以有限長移動線熱源模型（MFLS）與格子單松弛模型為基礎(chǔ)，建立含水層非穩(wěn)態(tài)過余溫度解析解ΔTMFLS與LBM耦合計算模型。通過搭建地埋管井群模擬實驗系統(tǒng)，將理論計算與熱響應(yīng)實驗相結(jié)合，預(yù)測含水層水動力場與溫度場的演化過程；從介觀尺度分析孔隙結(jié)構(gòu)與滲流速度對于含水層傳熱性能的影響。以期為探究地下水滲流對于改變地埋管井群傳熱特性與熱泵機組運行效率提供理論依據(jù)與技術(shù)支持。

1 理論基礎(chǔ)與控制方程

1.1 含水層滲流-傳熱控制方程

含水層熱-流運移過程滿足下述假設(shè)條件：1）含水層的孔隙結(jié)構(gòu)不隨滲流與傳熱過程發(fā)生變化；2）含水層為各向異性、均質(zhì)多孔介質(zhì)，水文地質(zhì)與熱物性參數(shù)不隨時間變化；3）地下水沿同一方向勻速滲流，忽略含水層垂向滲流過程以及垂向溫度梯度；4）含水層中固、液相基質(zhì)在傳熱過程中瞬間達到局部熱平衡狀態(tài)，具有相同溫度。

根據(jù)假設(shè)條件，基于達西定律與多孔介質(zhì)傳熱、傳質(zhì)理論［14］，建立含水層三維非穩(wěn)態(tài)質(zhì)量守恒、動量守恒以及能量守恒控制方程，如式（1）～式（3）所示。含水層傳熱過程包括：1）含水層中液相基質(zhì)與固相基質(zhì)間的熱傳導(dǎo)；2）地下水滲流引起的對流傳熱；3）由于含水層孔隙結(jié)構(gòu)差異所引起的地下水流動方向與強度的變化，造成熱量平均化，從而產(chǎn)生的熱彌散效應(yīng)。

[[εsγf+（1-εs）γs]?h?t+εs?v?x=0] （1）

[v=-kρfgεsμf??h?x] （2）

[ρece?T?t+ρfcfεsv?T?x=λx?2T?x2+λy?2T?y2+λz?2T?z2+QT] （3）

在能量守恒控制方程中，含水層有效比熱容ρece通過含水層固相基質(zhì)與液相基質(zhì)的體積比熱容加權(quán)平均（式（4））確定。含水層各方向的熱傳導(dǎo)系數(shù)由式（5）確定，其中有效導(dǎo)熱率λe為固相基質(zhì)與液相基質(zhì)導(dǎo)熱率的加權(quán)平均值；熱機械彌散系數(shù)張量采用地下水滲流速度一次方模型，由式（6）確定。

[ρece=εsρfcf+（1-εs）ρscs] （4）

[λx=εsλf+（1-εs）λs+λdispx=λe+λdispxλy=εsλf+（1-εs）λs=λeλz=εsλf+（1-εs）λs=λe] （5）

[λdisp=ρfcfεsαTvδxy+（αL-αT）vxvyv] （6）

1.2 含水層過余溫度解析解

在含水層儲能以及地源熱泵系統(tǒng)中，抽水-回灌井與地埋管換熱器垂向長度通常為井孔直徑的400倍以上，因此熱（冷）源在含水層中的傳熱過程可概化為半無限大多孔介質(zhì)中移動線熱源（MLS）的傳熱過程?；诜欠€(wěn)態(tài)點熱源格林函數(shù)形式［15］，通過Laplace變化對含水層能量守恒偏微分方程（式（3））中t、x進行轉(zhuǎn)化；通過Fourier變化對y、z進行轉(zhuǎn)化，進而聯(lián)立求解得到無限長移動線熱源模式（MILS）下含水層某點（x，y，z）過余溫度瞬態(tài)解析解ΔTMILS，如式（7）所示。

[ΔTMILS（x，y，z，t）=QT8ρeceλxλyλzρeceπt32?0tt-ξ-3/2exp-ρecex-x-εsvρfcft-ξ24ρeceλxt-ξ2-ρecey-y24λyt-ξ2-ρecez-z24λzt-ξ2dξ] （7）

在含水層能量守恒方程中源（匯）項通常由多組熱（冷）源構(gòu)成，且放（?。峁β手饡r變化。因此，在單一MILS模式解析解ΔTMILS（式（7））的基礎(chǔ)上，通過虛擬熱源鏡像方法［16］得到有限長單一線熱源在時間步長Δt內(nèi)引起含水層瞬態(tài)溫度響應(yīng)模型。根據(jù)時空疊加原理，針對多個非穩(wěn)態(tài)有限長移動線熱源（MFLS）運行模式，推導(dǎo)出含水層過余溫度解析解ΔTMFLS（式（8）），解析解中積分函數(shù)f （x， y， z， Δt）為指數(shù)互補誤差函數(shù)形式（式（9））。

[ΔTMFLS（x，y，z，t）=k=1mj=1nlqk（j?Δt）-lqk（j-1）?Δt2πλxλe?expρfcfεsv（x-x）2λx0lf（x，y，z，Δt）dz--l0f（x，y，z，Δt）dz] （8）

[f（x，y，z，Δt）=14（x-x′）2+（y-y′）2+（z-z′）2·" " " " " " " " " " " exp-ρfcfεsv（x-x′）2λxλx+（y-y′）2λxλe+（z-z′）2λxλe2erfc（x-x′）2λxλx+（y-y′）2λxλe+（z-z′）2λxλeρece-ρfcfεsvΔtλx2ρeceλe-1Δt+" " " " " " " " " " " expρfcfεsv（x-x′）2λxλx+（y-y′）2λxλe+（z-z′）2λxλe2erfc（x-x′）2λxλx+（y-y′）2λxλe+（z-z′）2λxλeρece+ρfcfεsvΔtλx2ρeceλe-1Δt] （9）

2 格子Boltzmann方法

LBM模擬計算流體運動過程中，通過引入虛擬粒子對流體速度進行離散，利用格子Boltzmann演化方程確定各網(wǎng)格節(jié)點處粒子的遷移和碰撞過程。本文基于格子單松弛模型（lattice Bhatnagar-Gross-Krook，LBGK）標(biāo)準(zhǔn)速度場的演化方程［17］，同時考慮含水介質(zhì)對于滲流溶液的外部體積力，在演化方程碰撞項引入離散力源項Fi，從而得到具有外力項的LBGK速度演化方程（式（10））。演化方程中fi（x， t）為t時刻x處速度矢量為Ci的速度分布函數(shù)；平衡態(tài)分布函數(shù)fieq（x， t）與離散力源項Fi如式（11）、 式（12）所示。

[fi（x+CiΔtL，t+ΔtL）-fi（x，t）=-1τ[fi（x，t）-feqi（x，t）]+FiΔtL]

（10）

[feqi（x，t）=wiρf（x，t）1+Ci?vC2s+（Ci?v）22C4s-v22C2s] （11）

[Fi=wi1-12τCi-vC2s+（Ci?v）C4sCi?Ff] （12）

Ff為含水層多孔介質(zhì)對滲流溶液產(chǎn)生的外部體積力，如式（13）所示；含水層滲透率k如式（14）［18］所示。

[Ff=-υfkv-1.75150εskvv+ρfg] （13）

[k=dp2ε3s1801-εs2] （14）

滲流溶液局部密度ρf（x， t）、局部流速v（x， t） 分別由速度分布函數(shù)fi（x， t）的第0和第1階矩確定，如式（15）、式（16）所示；通過建立運動黏滯系數(shù)υf與格子參數(shù)、無量綱松弛時間的關(guān)系式（式（17）），從而將微觀與宏觀流體運動過程相聯(lián)系。

[ρf（x，t）=ifi（x，t）] （15）

[ρf（x，t）v（x，t）=iCifi（x，t）+ΔtL2Ff] （16）

[υf=C23（τ-0.5）ΔtL] （17）

3 實驗研究與模型驗證

3.1 實驗系統(tǒng)與方案

本文利用實驗室已建成使用的地源熱泵實驗系統(tǒng)開展可控物理模擬實驗。將天津市已建成的地源熱泵系統(tǒng)中部分地埋管井群作為工程原型，設(shè)計、搭建一套實驗室尺度下便于操作、組裝的實驗系統(tǒng)，如圖1所示［19］。該實驗系統(tǒng)包括可變水頭滲流砂箱以及多熱源地埋管換熱系統(tǒng)、滲流砂箱水動力及熱源系統(tǒng)、熱泵機組與空調(diào)末端、系統(tǒng)測試與控制平臺等?；跓嶙枧c熱容優(yōu)化模型（TRCM）［20］，對含水層以及井孔內(nèi)、外的傳熱模型進行優(yōu)化；根據(jù)方程分析法，設(shè)計滲流砂箱及地埋管換熱器［21］。為保證滲流砂箱裝置可再現(xiàn)地埋管在含水層中的傳熱過程，要求工程原型與實驗系統(tǒng)的地埋管換熱器中的普朗時數(shù)Pr、雷諾數(shù)Re、傅里葉數(shù)Fo和努賽爾數(shù)Nu以及含水層中Fo、Nu和佩克萊數(shù)Pe等相似準(zhǔn)則數(shù)相同。

為確保砂箱含水層的滲流速度均勻穩(wěn)定，將滲流砂箱設(shè)計為長、寬、高為1.2 m×0.8 m× 1.0 m的長方體，其中滲流區(qū)域長為1 m，兩端分別對稱設(shè)置長為0.1 m的進（排）水區(qū)，在進（排）水端表面沿中軸線開設(shè)6個間隔為0.1 m的出水孔（Φ20 mm），如圖2所示。實驗過程中，通過開啟進（排）水區(qū)域不同高度的出水孔橡膠塞，控制進水與排水端的水頭差值，滿足不同滲流速度要求。選用GDH0510系列恒溫水浴作為冷熱源設(shè)備，確保砂箱填充介質(zhì)的初始溫度與上游邊界溫度保持恒定。滲流砂箱、供水箱、回收水箱以及連接管道均貼有10 mm厚的橡塑保溫材料。

結(jié)合實驗方案，首先將兩根長1 m、內(nèi)徑8 mm、外徑10 mm、線圈直徑80 mm的不銹鋼螺旋管式地埋管換熱器固定在砂箱內(nèi)部。將取自地下含水層的天然原砂經(jīng)過除雜、干燥后，利用濕填法進行等重度填充，使得砂箱內(nèi)含水介質(zhì)重度達到（1.6±0.05） kg/L，確保與地下含水層重度相近。砂箱中含水層的填充高度為0.8 m，在其上部分別鋪設(shè)厚度為0.1 m的黏土層與礫石層，將含水介質(zhì)與外部環(huán)境隔離。在距離砂箱底板0.5 m處，埋設(shè)13組精度為±0.1 ℃的K型熱電偶，地埋管換熱器與熱電偶測點分布如圖3所示。

以0.1 m水頭差由滲流砂箱上游連續(xù)通入15 ℃去離子水，對含水層進行飽水排氣處理。當(dāng)砂箱整體溫度穩(wěn)定為15 ℃且流出溶液的體積穩(wěn)定，測壓管內(nèi)無氣泡與固體顆粒出現(xiàn)時，認(rèn)為砂箱飽水排氣完成。通過實驗系統(tǒng)控制平臺將熱泵機組切換至制冷模式，針對存在地下水滲流模式，開展螺旋管換熱器連續(xù)排熱實驗。實驗運行周期為400 min，每6秒輸出一組螺旋管換熱器進出口水溫以及含水層各觀測點溫度。實驗過程中，確保砂箱進出溶液的水頭差值保持0.2 m，每10分鐘采用秒表-量筒法測試砂箱的出水流量。通過反演計算，得到在實驗階段砂箱含水層的滲流速度在（6～8）×10-6 m/s之間。

3.2 求解方法與模型驗證

將地源熱泵實驗系統(tǒng)中可控水頭滲流砂箱作為物理模型，設(shè)定計算區(qū)域水平面為1×0.8 m2，垂直方向計算范圍0～0.8 m。計算過程中含水層被劃分為8組水平斷層，每組斷層厚度為0.1 m。將計算區(qū)域上部黏土層與礫石層、底部以及左右兩個側(cè)面定義為隔水、絕熱邊界；沿地下水滲流方向前后兩個側(cè)面均作為定水頭邊界；兩根螺旋管式地埋管換熱器作為變熱流密度有限長線熱源處理。

基于所建立的演化方程與解析模型，在計算軟件Matlab上開發(fā)LBM-ΔTMFLS耦合求解程序。本文基于LBGK模型提出包含外力項的速度場演化方程（式（10）），因此根據(jù)實驗系統(tǒng)含水層中的實際滲流過程，選擇三維空間具有18個離散速度的D3Q19模型對演化方程進行離散求解。每個方向的離散速度和權(quán)重因子分別為：

[Ci=0，0，" i=0±1，0，0C，0，±1，0C，0，0，±1C，" i=1～6±1，±1，0C，0，±1，±1C，±1，0，±1C，" i=7～18] （18）

[wi=1/3，" "i=01/18，" "i=1～61/36，" " i=7～18] （19）

由于計算模型為1 m×0.8 m×0.8 m的長方體，根據(jù)所確定的格子距離，將計算區(qū)域離散為4096×1010個節(jié)點（40000×32000×32000）。將含水層上下、前后邊緣均設(shè)置為封閉的實體邊界，采用二階精度的修正反彈邊界格式處理；滲流溶液從左側(cè)進入，右側(cè)流出，因此左右側(cè)面采用非平衡態(tài)外推格式處理。根據(jù)文獻［20］與實驗結(jié)果，確定含水層基本物性參數(shù)與LBM基本設(shè)置參數(shù)，如表1所示。由于LBM與ΔTMFLS解析解均選用顯式循環(huán)算法，且具有各自的時間步長，因此在求解過程中將ΔTMFLS解析解時間步長Δt縮短以確保Δt/ΔtL為整數(shù)，進而可將LBM循環(huán)作為ΔTMFLS的子循環(huán)從而確保耦合計算時間的統(tǒng)一。

將螺旋管換熱器均等效為有限長線熱源，把實驗所得螺旋管逐時換熱量作為邊界條件導(dǎo)入LBM-ΔTMFLS耦合求解模型中，計算得到觀測點5、7、9的逐時過余溫度。同時選擇過余溫度的均方根誤差（root mean square error，RMSE）作為計算值與實驗結(jié)果的相似度判定指標(biāo)。如圖4所示，耦合模型計算結(jié)果在運算周期中的動態(tài)變化與實驗數(shù)據(jù)基本一致，3組觀測點的RMSE（ΔT）均小于10%，最大標(biāo)準(zhǔn)誤差為0.55。因此認(rèn)為LBM-ΔTMFLS耦合求解方法可準(zhǔn)確預(yù)測在地下水滲流作用下含水層中的傳熱過程。

4 模擬計算與分析

4.1 含水層流場變化規(guī)律

將圖2所示滲流砂箱作為研究對象，應(yīng)用所建LBM-ΔTMFLS耦合求解模型，分別針對滲流速度為3×10-6、8×10-6 m/s，孔隙率為0.3、0.5這4類工況開展模擬計算。在這4類工況中均將地埋管換熱器等效為定熱流密度熱源，每延米排熱功率為120 W/m，運算周期設(shè)置為400 min。

圖5、圖6為在不同時刻下滲流砂箱埋深0.5 m處水平斷面水動力場分布的演化過程。由于在LBM的計算方法中首先需對流場與溫度場初始化，將流體粒子節(jié)點均標(biāo)記為初始流速，因此在初始狀態(tài)（0 min）下含水層水平斷面流體粒子的滲流速度與入口邊界處賦值相同。當(dāng)耦合求解程序啟動后在4種計算工況中，水平斷面中下游滲流速度均出現(xiàn)大幅降低；隨著模擬計算的進行，含水層滲流速度逐漸回升，至運算期結(jié)束（400 min）含水層水動力場的整體分布趨于均勻一致。

通過對比發(fā)現(xiàn)，隨著含水層孔隙率與進水側(cè)邊界流速的降低，含水層中后部滲流速度的變化幅度更加明顯。通過分析認(rèn)為，滲流速度出現(xiàn)以上的演化過程是由構(gòu)成含水層的多孔介質(zhì)幾何結(jié)構(gòu)的復(fù)雜性與多樣性所造成的。當(dāng)滲流溶液由邊界層進入到多孔介質(zhì)區(qū)域，由于流動空間的急劇減少，流體的流線方向發(fā)生改變，然而虛擬流體粒子在慣性作用下會保持原有的運動方向，引起粒子之間以及粒子與多孔介質(zhì)固相基質(zhì)之間發(fā)生碰撞，從而造成動能的損失，導(dǎo)致滲流速度下降。與此同時，含水層中孔隙結(jié)構(gòu)大小不一也會使不同孔隙間沿軸向的最大流速存在差異；在孔隙發(fā)育程度較差的不連通或半連通區(qū)域內(nèi)，會引起流體粒子的停滯，導(dǎo)致局部滲流速度驟降。

由于計算過程中假設(shè)含水層滲透率、孔隙率等水文地質(zhì)參數(shù)均為定值，因此所確定的LBGK速度演化方程（式（10））中外部體積力Ff（式（12））不變，隨著滲流溶液不斷流入砂箱，含水層中下游滲流速度逐漸回升。與采用空間平均方法從宏觀尺度描述滲流過程的有限元法不同，LBM計算過程中可根據(jù)邊界條件模擬流體粒子遷移與碰撞過程，從而確定格子新的分布函數(shù)，得到在不同位置、不同時間步長下的流體粒子局部滲流速度。因此，在含水層水動力場中可直觀反映出流體粒子所在滲流通道的演化過程與孔隙結(jié)構(gòu)的基本形貌特征。

分別選擇距離砂箱進水端300、600與900 mm的斷面作為研究對象，通過計算確定各斷面全部流體粒子的平均流速，從而擬合生成動態(tài)變化曲線。如圖7、圖8所示，在4種工況下滲流速度擬合曲線均出現(xiàn)啟動—下降—回升—穩(wěn)定的4個連續(xù)階段。以進水流速3×10-6 m/s為例，當(dāng)孔隙率為0.3時，3組斷面平均流速最大降幅分別達到40%、43%、50%，在計算周期結(jié)束時穩(wěn)定在2.7×10-6、2.5×10-6、2.4×10-6 m/s；當(dāng)孔隙率提高到0.5時，3組斷面平均流速最大降幅依次為27%、37%、43%，且僅在砂箱末端900 mm斷面處的平均流速穩(wěn)定在2.9×10-6 m/s，前兩組斷面的平均流速均回升到3×10-6 m/s。

對比圖7與圖8可發(fā)現(xiàn)，在較低的滲流速度下，隨著含水層孔隙率的降低，各斷面平均流速的降幅度增大，回升過程緩慢。計算表明，含水層孔隙結(jié)構(gòu)的變化直接決定了流體粒子的運動空間，進而影響流體粒子的運動軌跡。隨著孔隙率的降低，流體粒子運動空間減少，繞過多孔介質(zhì)的運動路程增大，粒子之間以及粒子與多孔介質(zhì)的碰撞幾率增大，動能損失增大。然而隨著進水流速的提高，進水側(cè)邊界動能增大，導(dǎo)致孔隙率的變化對于含水層水動力場演化過程的影響程度減弱。

4.2 含水層溫度場變化規(guī)律

通過LBM-ΔTMFLS耦合求解方法，計算得到孔隙率為0.3，運算周期結(jié)束時砂箱埋深0.5 m含水層中20、25、30、35 ℃這4組等溫線分布。為了準(zhǔn)確比較兩類滲流速度下含水層溫度場的分布差異，將含水層初始溫度升高5 ℃的區(qū)域定義為熱擴散范圍。如圖9所示，當(dāng)砂箱進水側(cè)流速為3×10-6 m/s時，含水層熱擴散范圍以線熱源為中心向四周呈對稱輻射狀擴散；當(dāng)進水側(cè)流速增至8×10-6 m/s時，4組等溫線沿地下水滲流方向的偏移程度增大，位于線熱源下游的熱擴散半徑增大0.11 m，線熱源上游的熱擴散半徑則減少0.06 m。隨著滲流速度的增大，含水層熱量運移過程的方向性顯著增強。

以觀測點5、7、9為例，計算得到各觀測點溫度響應(yīng)的動態(tài)變化曲線，如圖10所示。隨著含水層滲流速度的提高，位于熱源中心區(qū)域的觀測點7的過余溫度下降1.7 ℃。由于增強地下水滲流過程可及時遷移熱源周圍堆積的熱量，加快沿滲流方向的溫度鋒面運移速度，因此位于下游觀測點9的過余溫度升高8 ℃；然而位于上游的觀測點5過余溫度下降1 ℃?；诤畬觽鳠崮Ｐ团c過余溫度解析解ΔTMFLS，當(dāng)水文地質(zhì)參數(shù)與熱物性參數(shù)不發(fā)生變化時，對流傳熱強度與熱彌散效應(yīng)均隨滲流速度的上升而增強，從而向下游區(qū)域的熱擴散能力增強。

5 結(jié) 論

1）基于格子單松弛模型，引入離散力源項，得到具有外力項的LBGK速度演化方程；根據(jù)有限長移動線熱源模型，針對多個非穩(wěn)態(tài)有限長移動線熱源運行模式，利用時空疊加原理推導(dǎo)出含水層過余溫度ΔTMFLS的解析解，建立含水層LBM-ΔTMFLS耦合計算模型。通過滲流砂箱的熱響應(yīng)實驗，將測試數(shù)據(jù)與計算結(jié)果相比較，3組觀測點過余溫度的RMSE（ΔT）均小于10%，因此認(rèn)為LBM-ΔTMFLS耦合求解方法可準(zhǔn)確反映在地下水滲流作用下含水層中的傳熱過程。

2）通過對滲流砂箱開展模擬計算，在不同工況下含水層滲流速度均具有啟動—下降—回升—穩(wěn)定4個連續(xù)過程。當(dāng)進水流速為3×10-6 m/s，孔隙率由0.5降至0.3時，3組斷面平均流速的最大降幅分別由27%、37%、43%升至40%、43%、50%。研究表明，隨著含水層孔隙率的降低，流體粒子運動空間減少，粒子之間以及粒子與多孔介質(zhì)的碰撞幾率增加，動能損失增大，從而導(dǎo)致滲流速度的降幅增大，回升過程緩慢。然而隨著進水流速的上升，孔隙率的變化對于含水層水動力場演化過程的影響程度減弱。

3）LBM-ΔTMFLS耦合計算過程中可根據(jù)邊界條件模擬流體粒子遷移與碰撞過程，確定格子新的分布函數(shù)，從而直觀反映流體粒子的運動軌跡與所在含水層孔隙結(jié)構(gòu)的基本形貌特征。當(dāng)含水層孔隙結(jié)構(gòu)不變時，隨著進水流速的上升，流體粒子初始動能增加，含水層不同區(qū)域的滲流速度回升幅度增大，最終趨于進水流速值。因此含水層中的對流傳熱強度與熱彌散效應(yīng)均隨之增強，向下游區(qū)域的熱擴散能力增強。

符 號

Ci i方向格子速度矢量

Cs 格子聲速，Cs=C/30.5

c 定壓比熱容，J/（kg·K）含水介質(zhì)顆粒平均直徑，mm

Ff 含水介質(zhì)對滲流溶液的外部體積力，N

g 重力加速度，m/s2

h 含水層水頭，m

k 含水層滲透率，m2

l 線熱源長度，m

Q 地埋管換熱器換熱功率，W

QT 熱源（匯）項排（?。釓姸?，W

q 地埋管換熱器單位埋深換熱功率，W/m

T 溫度，℃

t 時間，s

v 滲流速度，m/s

wi i方向的權(quán)重因子

x、y、z 直角坐標(biāo)方向，m

x′、y′、z′ 熱源所在坐標(biāo)點，m

λ 導(dǎo)熱率，W/（m·K）

ε 孔隙率

ρ 密度，kg/m3

γ 壓縮系數(shù)，（m·s）/ kg

μ 動力黏滯系數(shù)，Pa·s

υ 運動黏滯系數(shù)，m2/s

αT、αL 熱儲層橫向、縱向熱彌散度，m

ξ 時間積分變量，s

τ 無量綱松弛時間

δxy 克羅內(nèi)克張量

η 地埋管換熱器能效系數(shù)

下標(biāo)

k 線熱源數(shù)，1～m

j 時間步長數(shù)，1～n

in 滲流砂箱進水側(cè)

s 含水層固相基質(zhì)

f 含水層滲流溶液

［參考文獻］

［1］ 中國電子信息產(chǎn)業(yè)發(fā)展研究院. 碳中和愿景下儲能產(chǎn)業(yè)發(fā)展白皮書［M］. 北京： 中國電子信息產(chǎn)業(yè)發(fā)展研究院出版社， 2021.

China Institute of Electronic Information Industry Development. White paper on energy storage industry development under the vision of carbon neutralization ［M］. Beijing： China Institute of Electronic Information Industry Development Publishing House， 2021.

［2］ 國家發(fā)展和改革委員會能源局， 等. 關(guān)于促進地?zé)崮荛_發(fā)利用的若干意見［M］. 北京： 中國計劃出版社， 2021.

National Development and Reform Commission Energy Reform Institute， et al. Several opinions on promoting the development" and" utilization" of" geothermal" energy［M］. Beijing： China Planning Press， 2021.

［3］ 王貴玲， 劉彥廣， 朱喜， 等. 中國地?zé)豳Y源現(xiàn)狀及發(fā)展趨勢［J］. 地學(xué)前緣， 2020， 27（1）： 1-9.

WANG G L， LIU Y G， ZHU X， et al. The status and development trend of geothermal resources in China［J］. Earth science frontiers， 2020， 27（1）： 1-9.

［4］ 清華大學(xué)建筑節(jié)能研究中心. 中國建筑節(jié)能年度發(fā)展研究報告（2020） ［M］. 北京： 中國建筑工業(yè)出版社， 2020.

Tsinghua University Building Energy Conservation Research Center. Annual development report of building energy conservation in China （2020）［M］. Beijing： China Building Industry Press， 2020.

［5］ SMITH D C， ELMORE A C. The observed effects of changes in groundwater flow on a borehole heat exchanger of" a" largescale" ground" coupled" heat" pump" system［J］. Geothermics， 2018， 74： 240-246.

［6］ BADRUDDIN I A， AZEEM K， YUNUS K T M， et al. Heat transfer in porous media： a mini review［J］. Materials today： proceedings， 2020， 24（Part2）： 1318-1321.

［7］ ZHANG M Y， ZHAO Q Y， HUANG Z J， et al. Numerical simulation of the drag and heat-transfer characteristics around and through a porous particle based on the lattice Boltzmann method［J］. Particuology， 2021， 58： 99-107.

［8］ BAKHSHIAN S， HOSSEINI S A， SHOKRI N. Pore-scale characteristics of multiphase flow in heterogeneous porous media using the lattice Boltzmann method［J］. Scientific reports， 2019， 9（1）： 464-467.

［9］ MANDZHIEVA R， SUBHANKULOVA R. Practical aspects of absolute permeability finding for the lattice Boltzmann method and pore network modeling［J］. Physica A： statistical mechanics and its applications， 2021， 582： 126249.

［10］ ZHU X F， WANG S， FENG Q H， et al. Pore-scale numerical prediction of three-phase relative permeability in porous media using the lattice Boltzmann method［J］. International communications in heat and mass transfer， 2021， 126： 105403.

［11］ FENG X B， LIU Q， HE Y L. Numerical simulations of convection heat transfer in porous media using a cascaded lattice Boltzmann method［J］. International journal of heat and mass transfer， 2020， 151： 119410.

［12］ PARVAN A， JAFARI S， RAHNAMA M， et al. Insight into particle retention and clogging in porous media： a pore" scale" study" using" lattice" Boltzmann" method［J］. Advances in water resources， 2020， 138： 103530.

［13］ ZHANG Y T， JIANG F， TSUJI T. Influence of pore space heterogeneity on mineral dissolution and permeability evolution investigated using lattice Boltzmann method［J］. Chemical engineering science， 2022， 247： 117048.

［14］ MA J C， JIANG Q， ZHANG Q L， et al. Effect of groundwater forced seepage on heat transfer characteristics of" "borehole" "heat" "exchangers［J］." Geothermal" "energy， 2021， 9（1）： doi： 10.1186/S40517-021-00192-1.

［15］ SELCUK E， BERTRAND F. Multilayer analytical model for vertical ground heat exchanger with groundwater flow［J］. Geothermics， 2018， 71： 294-305.

［16］ HUAN J Z. Improved analytical model for vertical borehole ground heat exchanger with multiple-layer substrates and groundwater flow［J］. Applied energy， 2017， 202： 537-549.

［17］ MCCULLOUGH J W S， AMINOSSADATI S M， LEONARDI C R. Transport of particles suspended within a temperature-dependent viscosity fluid using coupled LBM-DEM［J］. International journal of heat and mass transfer， 2020， 149： 119159.

［18］ LIU X C， HUANG H B， LU X Y. Lattice Boltzmann study of effective viscosities of porous particle suspensions［J］. Computers and fluids， 2019， 181： 135-142.

［19］ 馬玖辰， 王文君， 王宇， 等. 基于含水層熱-滲運移機理的地源熱泵實驗系統(tǒng)研發(fā)［J］. 實驗室研究與探索， 2020， 39（8） ： 88-93.

MA J C， WANG W J， WANG Y， et al. Research and development of ground source heat pump experimental system based on the mechanism of aquifer groundwater seepage" " and" " heat" " transferring［J］." " Research" " "and exploration in laboratory， 2020， 39（8）： 88-93.

［20］ 馬玖辰， 王文君， 魏璠， 等. 熱力彌散對地埋管換熱器所在含水層傳熱過程的影響［J］. 太陽能學(xué)報， 2021， 42（3）： 164-170.

MA J C， WANG W J， WEI F， et al. Influence of thermal dispersion on transfer process in aquifers around borehole heat exchangers［J］. Acta energia solaris sinica， 2021， 42（3）： 164-170.

［21］ YANG G C， JING L， KWOK C Y， et al. A comprehensive parametric study of LBM-DEM for immersed granular flows［J］. Computers and geotechnics， 2019， 114： 103100.

NUMERICAL SIMULATION OF FLUID SEEPAGE AND HEAT TRANSFER IN AQUIFER WITH LBM-MFLS COUPLED MODEL

Ma Jiuchen1，2，Lyu Linhai1，Yang Jie1，Cui Afeng1，Wei Fan1，2

（1. College of Energy Safety Engineering， Tianjin Chengjian University， Tianjin 300384， China；

2. Key Laboratory for Efficient Use of Low and Medium Grade Energy， Ministry of Education of China，

Tianjin University， Tianjin 300072， China）

Abstract：Taking the multiple moving finite line heat sources（MFLS） operation mode as the research object， the unsteady state analytical solutions of the excess temperature in aquifer ΔTMFLS are obtained by applying space-time superposition principle， based upon the transient moving finite line heat source model. A calculation model coupled lattice Boltzmann method （LBM） and ΔTMFLS is established， according to introduce the discrete force source term into the Lattice Bhatnagar-Gross-Krook （LBGK） evolution equation of the velocity distribution function. The coupled model and the calculation method are validated by the data determined from the in-situ thermal re-sponse test. The results show that coupled LBM-ΔTMFLS simulation can well reproduce the evolution process of seepage velocity in aquifer which presents four successive stages of starting， declining， rising and stabilizing. With the decreasing of the porosity of aquifer， the movement space of virtual fluid particles reduces， consequently the collision probability of inter-particles and between particles and porous media increases， which leads to the kinetic energy loss enhancing， the descend range of seepage velocity increasing， and the recovery process retarding. While the inflow velocity increases， however， the influence degree of porosity on the hydrodynamic evolution of aquifer weakens. Meanwhile the directivity of the heat transport process in the aquifer enhances significantly.

Keywords：geothermal energy； heat transfer； porous media； numerical models； aquifers； lattice Boltzmann method