模擬退火算法在巖土熱物性參數確定中的應用

張長興，王德水，劉玉峰，孫始財，彭冬根

（1 山東科技大學山東省土木工程防災減災重點實驗室，山東 青島 266590；2 青島大學基建處，山東 青島 266071；3 南昌大學建筑工程學院，江西 南昌 330031）

引 言

在土壤源熱泵系統(tǒng)設計時，獲取建設地點的巖土熱物性參數是實施地埋管換熱器（borehole heat exchanger，BHE）設計的先決條件。因此，中國2009年修訂的《地源熱泵系統(tǒng)工程技術規(guī)范》（以下簡稱《規(guī)范》）對熱物性參數的測試和確定方法進行了明確規(guī)定，該測試又稱為巖土的熱響應試驗(thermal response test,TRT)[1]。通過TRT 試驗確定的熱物性參數主要是巖土的熱導率和容積比熱容，同時由于施工工藝、填充材料、U 形管的布置形式等客觀因素的影響，確定現場地埋管換熱器的實際熱阻值也是非常必要的[2]。

根據地埋管換熱器的傳熱原理，通過TRT 試驗的方法確定巖土熱物性參數是典型的熱傳導反分析問題，該理論由Shumakov 在1957年提出，很多研究者對簡單形體的熱傳導反問題的數值解進行了研究[3]。隨著計算機技術的發(fā)展，熱傳導反問題的求解也日趨成熟，其求解過程的本質是最優(yōu)化問題，它包含兩方面的內容：其一是建立數學模型，即用數學語言描述最優(yōu)化問題，模型中的數學關系式反映了最優(yōu)化問題所要達到的目標和各種約束條件；其二是數學求解，數學模型建好以后選擇合理的最優(yōu)化方法進行求解。最優(yōu)化方法主要包括基于迭代計算的數值方法和智能算法兩類，數值方法主要有牛頓法、最速下降法、共軛梯度法等[4]，智能算法主要有遺傳算法、混沌優(yōu)化算法、模擬退火算法、混合算法等[5]。本研究通過建立TRT 試驗系統(tǒng)模型，利用模擬退火算法(simulated annealing algorithm,SAA)確定變熱流條件下巖土的熱導率、容積比熱容和地埋管換熱器的實際熱阻，為土壤源熱泵系統(tǒng)的設計和應用提供重要的基礎數據。

1 巖土熱物性參數識別的數學模型

要獲得準確的巖土熱物性參數，參數識別的優(yōu)化模型必須可靠、準確。這就對TRT 試驗系統(tǒng)的數學模型提出了更高的要求，同時要求實施優(yōu)化的目標函數對于待識別參數的變化有足夠的敏感性，通過目標函數獲得最小值來最終確定最優(yōu)的熱物性參數組合。

1.1 TRT 試驗系統(tǒng)的數學模型

系統(tǒng)的數學模型是實施最優(yōu)化技術的關鍵，其準確性直接影響參數識別的精度。對于TRT 試驗系統(tǒng)而言，其主要的3 個組成部分分別為地埋管換熱器、電加熱器和循環(huán)水泵，土壤與循環(huán)水通過地埋管換熱器完成熱交換，電加熱器用以控制換熱器加熱負荷的強度，循環(huán)水泵控制換熱器熱交換的循環(huán)水量，如圖1所示。其中，地埋管換熱器傳熱過程的數學建模是最復雜的[6]。根據地埋管換熱器的物理模型，其傳熱分析一般是以鉆孔壁為界分成兩個計算區(qū)域，孔內區(qū)域按穩(wěn)態(tài)傳熱計算，孔外區(qū)域按非穩(wěn)態(tài)傳熱計算。

圖1 巖土熱響應測試系統(tǒng)示意圖Fig.1 Schematic diagram of rock-soil thermal response test system

1.1.1 鉆孔外傳熱的數學模型 地埋管換熱器鉆孔外傳熱的數學模型主要有解析解模型和數值解模型兩類。在進行TRT 試驗系統(tǒng)建模時，解析解模型中的無限長線熱源模型（line source model，LSM）、有限長線熱源模型（finite line source model，FLSM）和柱熱源模型（cylinder source model，CSM）均得到了廣泛的應用[2,6]，同時一維、二維和三維的地埋管換熱器數值模型在進行熱物性參數確定時亦顯示了一定的優(yōu)勢[7-10]。然而，考慮到現場巖土TRT 試驗的復雜條件和測試現場的電壓波動等因素的影響，測試過程中經常出現《規(guī)范》要求的恒熱流和定循環(huán)水流量難以保證的情況，如重新測試，需等待10～14 d 的時間[11]，會對測試的成本和工程進度造成影響，這對地埋管換熱器的數學模型提出了更高的要求。針對解析解模型在變熱流條件下計算速度慢[12-13]而數值模型難以準確確定換熱器熱阻的缺點[10]，本研究以應用格林函數法求解的柱熱源模型作為鉆孔外傳熱的數學模型，在變熱流條件下應用SAA 優(yōu)化算法實施熱物性參數的識別。

在變熱流條件下，由于柱熱源模型在進行Laplace 變換后的解析解在確定鉆孔壁溫時需實施負荷疊加才能確定即時的鉆孔壁溫，會影響熱物性參數識別的計算速度。本研究應用格林函數變換求解鉆孔壁溫的方法，通過計算方式的改變提高參數識別的速度[14]。通過引入格林函數，在數學計算方式的處理上可將積分計算進行離散求和的轉換，得到鉆孔壁溫的計算公式。

在時刻：

在+ 時刻：

式中，=αt/，Ts為土壤的原始溫度，Jn和Yn分別為第一類和第二類貝塞爾函數。

從公式可以看出，該卷積算法與傳統(tǒng)的G函數法[12]相比不再需要反復調用動態(tài)熱流和其對應的G函數，極大地縮短了計算時間，因此在土壤源熱泵系統(tǒng)長期動態(tài)運行特性的預測方面顯示了一定的 優(yōu)勢。

1.1.2 鉆孔內傳熱的數學模型 在地埋管換熱器的內部傳熱主要是U 形管內的水與鉆孔壁之間的傳熱，由于鉆孔內的熱阻與鉆孔外比較相對較小，在定熱流TRT 試驗中經常作為穩(wěn)態(tài)傳熱處理，然而這種熱阻的差異會造成在熱物性測試的初始階段換熱器平均水溫實測值與模擬值的差異，因此在實施系統(tǒng)優(yōu)化時通常舍去10 h 內的數據[2,15]，以保證辨識結果的準確性。

在變熱流TRT 試驗條件下，辨識土壤熱物性參數時，鉆孔內的熱阻是否可作為穩(wěn)態(tài)傳熱處理需要考慮熱流變化的時間間隔、熱流的變化幅度、填充材料的熱導率和體積比熱容等因素的綜合影響。確定鉆孔內穩(wěn)態(tài)熱阻時，利用形狀因子法計算鉆孔內填充材料的熱阻時綜合考慮U形管間熱短路現象的影響[16]，同時結合U 形管內的對流換熱、管壁的導熱這兩部分熱阻，鉆孔內穩(wěn)態(tài)熱阻可以表示為[17]

式中，β0、β1為擬合常數，對于本研究所選的單U 形地埋管換熱器形式，β0=17.44，1β=0.6052-；hc,f為管內流體與管內壁之間的對流傳熱系數，hc,f≈0.023Re0.8Pr0.3λf/dinn；Re=υdinn/v，為Reynolds數，與流體的流速υ、管內徑di和運動黏滯系數v有關；Pr=v/α，為Prandtl 數，與流體的運動黏滯系數v和導溫系數α有關。

從式(3)可以看出，實際工程中確定鉆孔內穩(wěn)態(tài)熱阻時，填充材料的熱導率會受到回填密實程度、地下水和土壤含濕量等因素的影響而造成一定的計算誤差，為了降低這一誤差，提高管內流體與管內壁之間的對流傳熱系數的計算精度是非常必要的。因此，在確定對流傳熱系數時，本研究采用循環(huán)水的熱物性隨水溫變化的動態(tài)參數，根據文獻[18]，循環(huán)水的密度、運動黏滯系數和熱導率隨水溫變化的關系式如下

根據地埋管換熱器傳熱原理和能量守恒定律，地埋管換熱器在時刻的出口水溫可以表示為[19]

1.2 優(yōu)化模型的建立

根據巖土TRT 試驗系統(tǒng)的數學模型，將試驗中的逐時加熱量和循環(huán)水流量作為輸入條件，結合地埋管換熱器的幾何條件和物性參數可計算出地埋管換熱器的逐時進出水溫度，通過對比地埋管換熱器的進出水溫度平均值與實測值對式(1)中的巖土熱導率和容積比熱容進行辨識，并根據式(3)確定地埋管換熱器的實際熱阻值。本研究采用進出水溫度平均值與實測值的均方根誤差（root-mean-square error，RMSE）作為優(yōu)化模型的目標函數

通過對地埋管換熱器傳熱過程的分析可以看出，地埋管平均出水溫度的計算是非線性的，會增加參數識別工作的難度。因此，本研究利用SAA算法實施熱物性參數的識別，以達到在變熱流工況下快速、方便和有效地確定巖土的熱導率、容積比熱容和地埋管換熱器熱阻的目的。

2 熱物性參數識別的SAA 算法設計

模擬退火算法（simulated annealing algorithm，SAA）是基于Mente Carlo 迭代求解策略的一種隨 機尋優(yōu)算法，其出發(fā)點是基于物理中固體物質的退火過程與一般組合優(yōu)化問題之間的相似性。SAA 算法在某一溫度下，伴隨溫度參數的不斷下降，結合概率突跳特性在解空間中隨機尋找目標函數的全局最優(yōu)解，即在局部優(yōu)解能概率性地跳出并最終趨于全局最優(yōu)。目前，SAA 算法已經廣泛應用在工農業(yè)生產中[20-22]。該算法用Metropolis 算法[23]產生組合優(yōu)化問題解的序列，并由與Metropolis 準則對應的轉移概率P確定是否接受從當前解到新解的轉移。

式(9)中的t∈R+表示控制參數。開始讓t取較大的值（與固體的熔解溫度相對應），在進行足夠多的轉移后緩慢減少t的值，如此重復，直至滿足某個停止準則使算法終止。因此，SAA 算法可視為遞減控制參數值時Metropolis 算法的迭代。圖2描述了SAA 算法的計算流程。

圖2 SAA 算法的流程圖Fig.2 Flow chart of simulated annealing algorithm

表1 單U 形地埋管換熱器相關參數Table 1 Related parameters of single U-pipe BHE

3 應用算例

3.1 TRT 試驗概況

為了準確確定巖土熱物性參數，在青島市經濟技術開發(fā)區(qū)進行了TRT 試驗，實測巖土的原始溫度Ts為15.55℃，測試的地埋管換熱器為單U 形PE 管形式，回填材料為10%膨潤土、90% SiO2砂子的混合物，測孔的幾何尺寸和相關材料的熱物性參數見表1。換熱器進出水溫度的測試時間間隔為5 min，測試時間共計55 h。由于現場用電條件的限制，測試中電壓的波動偏差高于5%，造成實測的循環(huán)水泵流量和電加熱功率出現了不同程度的波動（圖3），實測平均加熱功率為7460 W，標準差為228.3 W，不符合文獻[1,11]中的恒熱流加熱條件。在變熱流TRT 測試中，電加熱功率的最大偏差僅為8%，而且間隔時間較短，運行10 h 后將鉆孔熱阻作為穩(wěn)態(tài)熱阻處理。由于實測的加熱功率和循環(huán)水泵流量均隨電壓的波動發(fā)生變化，必須應用實測的動態(tài)電加熱功率和循環(huán)水流量才能準確確定巖土的熱物性參數和換熱器的熱阻。

圖3 現場TRT 試驗中的加熱功率Fig.3 Dynamic heat power of in-situ thermal response test

3.2 應用SAA 算法確定巖土熱物性參數

利用本研究建立的TRT 試驗系統(tǒng)模型，結合TRT 試驗結果，本研究運用SAA 算法確定了巖土的熱導率和容積比熱容。數學模型中的電加熱器功率、水泵流量和進出水溫度均采用試驗瞬時值，水的熱物性參數采用瞬時溫度對應的計算值[式(4)、式(5)和式(6)]，巖土的原始溫度、填充材料和PE 管的熱導率以及測試孔相關幾何尺寸均作為模型的已 知值。

在利用SAA 算法進行熱物性參數識別時，由于TRT 試驗系統(tǒng)模型中地埋管平均出水溫度的計算是非線性的，參數識別過程中一定存在非適定問題[24]，即可能有多個熱物性參數組滿足目標函數RMSE 最小這一條件。為了獲得熱物性參數的真實值，本研究首先設定巖土的熱導率和容積比熱容的識別范圍分別為[1,5]和[1000,5000]，對巖土的熱導率用0.05 的數據間隔在[1,5]范圍內生成81 個值，對容積比熱容用50 的數據間隔在[1000,5000]范圍內生成81 個值，使兩個參數值分別組合成為不同的數組，作為TRT 試驗系統(tǒng)模型的輸入值計算地埋管換熱器的逐時進出口平均水溫，結合TRT 試驗實測的對應值可計算出6561 個RMSE 值，RMSE 計算值的分布如圖4中的等值線所示?？梢钥闯觯赗MSE＜0.2 的范圍內存在6 處滿足RMSE 最小的極值點，為解決參數識別過程中的非適定問題奠定了基礎。根據現場TRT 試驗鉆孔勘察的地質資料，地下巖土以花崗巖、砂巖為主，依照《規(guī)范》提供的相關熱物性參數，并結合圖4的RMSE 值分布狀況，在實施 SAA 算法時設定巖土熱導率為 2 W·m-1·℃-1、容積比熱容為3000 kJ·m-3·℃-1，兩個參數的優(yōu)化范圍分別為[2,3]和[2000,4000]。由于實施SAA 算法需設定算法終止準則，考慮到TRT試驗測溫鉑電阻的測試誤差[25]，本研究設定目標函數RMSE 為0.14℃。

SAA 算法實施過程中，隨著退火溫度t的降低，RMSE 降至0.1361，算法終止，此時對應的巖土熱導率和容積比熱容分別為2.52 W·m-1·℃-1和2956 kJ·m-3·℃-1。如圖4所示，參數識別結果合理可信，結合現場的地質資料有效地解決了傳熱反問題分析中的非適定問題。圖5描述了SAA 計算中退火溫度和目標函數值的變化，可以看出，在退火溫度由100℃降至0.303℃的情況下目標函數值由1.54降至0.1361，達到了識別最優(yōu)熱物性參數的目的。

3.3 參數識別結果的分析

圖4 不同熱物性參數組合對應的RMSE 值分布Fig.4 RMSE distribution corresponding to different parameters’ array

圖5 參數識別過程中退火溫度和目標函數值的變化Fig.5 Annealing temperature and objective function variation in process of parameters identification

圖6 最優(yōu)熱物性參數對應的水溫模擬值與實測值的 溫差平方和動態(tài)變化Fig.6 Dynamic variation of square of temperature difference between simulating water temperature and experiment data corresponding to optimal thermal parameters

圖7 最優(yōu)熱物性參數對應的水溫模擬值與實測值對比Fig.7 Contrast between simulating water temperature and experiment data corresponding to optimal thermal parameters

圖8 地埋管換熱器熱阻隨平均水溫的動態(tài)變化Fig.8 Dynamic variation of borehole heat resistance with mean water temperature

將參數識別結果代入TRT 試驗系統(tǒng)的數學模型，由式(7)可計算出各時刻的地埋管換熱器進出水溫度Tin,sim,i、Tout,sim,i和平均水溫Tav,sim,i。圖6為最 優(yōu)熱物性參數對應地埋管換熱器平均水溫的模擬值和實測值溫差平方和的逐時變化圖，可以看出，在測試10 h 后各數值吻合程度較高，動態(tài)變化趨勢一致，溫差平方和最高值僅為0.15，出現在第41 h，溫差平方和的逐時變化也客觀上驗證了文獻[2,15]的正確性。圖7為10～55 h 各時刻地埋管換熱器平均水溫的實測值和模擬值的對比圖，可以看出，各時刻兩個變量高度吻合，擬合優(yōu)度R2為0.978。

兩個熱物性參數確定后，試驗系統(tǒng)模型由式(3)可計算出地埋管換熱器熱阻隨地埋管平均水溫的動態(tài)變化值，如圖8所示。隨著TRT 試驗的進行，地埋管換熱器的熱阻逐漸趨于穩(wěn)定，10 h 后熱阻的平均值為0.109 m·℃·W-1，換熱器的實際熱阻值可作為測試地點進行地埋管換熱器組群設計的計算 依據。

為了驗證SAA 算法參數識別結果的準確性，在同一測試地點按照《規(guī)范》的測試方法和數據處理方式確定的兩個參數分別為2.41 W·m-1·℃-1和2995 kJ·m-3·℃-1，與《規(guī)范》方法確定的兩個熱物性參數的相對誤差為4.1%和1.3%。在變熱流TRT 試驗系統(tǒng)模型中應用地埋管換熱器的柱熱源模型，SAA 算法在參數識別過程中具有較強的適用性，而且識別結果準確性高，為地埋管換熱器組群的精確設計提供了可靠的基礎數據。

4 結 論

本研究利用地埋管換熱器的柱熱源模型，通過引入格林函數法建立了變熱流條件下TRT 試驗系統(tǒng)模型，在此基礎上通過應用SAA 算法提出了一種巖土熱物性參數識別的方法，得到以下結論。

（1）在變熱流TRT 試驗系統(tǒng)模型中，應用格林函數法的柱熱源模型顯示了較好的適應性，為巖土熱物性參數的準確識別奠定了理論基礎。

（2）利用RMSE 值分布圖，有效地解決了傳熱反問題參數識別過程中的非適定問題。在此基礎上，結合TRT 試驗數據對巖土熱物性參數進行識別的過程中SAA 算法對應的退火溫度降溫速度快，確定的巖土熱導率和容積比熱容與《規(guī)范》方法對應值的相對誤差分別為4.1%和1.3%，證明了本研究方法的準確性。

（3）在巖土熱物性參數的識別過程中，可計算地埋管換熱器熱阻隨平均水溫的動態(tài)變化，確定測試現場換熱器的有效熱阻值，為地埋管換熱器組群的精確設計提供有效的基礎數據。

符 號 說 明

C——容積比熱容，kJ·m-3·℃-1

c——比熱容，J·kg-1·℃-1

d——直徑，m

G——G函數

H——鉆孔深度，m

h——對流傳熱系數，W·m-2·℃-1

J——第一類貝塞爾函數

M——M函數

m——水泵流量，kg·s-1

n——測試時間，h

P——加熱功率，W

Pr——Prandtl 數

ql——單位鉆孔深度的換熱量，W·m-1

R——熱阻，m·℃·W-1

Re——Reynolds 數

r——半徑，m

T——溫度，℃

t——時間，h

v——速度，m·s-1

Y——第二類貝塞爾函數

α——熱擴散系數，m2·h-1

λ——熱導率，W·m-1·℃-1

υ——運動黏滯系數，m2·s-1

ρ——密度，kg·m-3

下角標

av——平均

b——鉆孔

exp——試驗值

f——流體

g——填充材料

in——進口

inn——內部

o——外部

out——出口

p——管

s——巖土

sim——模擬值

[1]The Ministry of Housing and Urban-Rural Construction of China.Technical code for GCHPs.(GB 50366—2005) 2009 edition.Beijing:Chinese Building Industrial Press,2009

[2]Zhang Changxing (張長興),Guo Zhanjun (郭占軍),Liu Yufeng (劉玉峰),et al.Application of pattern search algorithm for determining heat resistance of ground heat exchanger [J].Transactions of Chinese Society of Agricultural Engineering(農業(yè)工程學報),2013,29 (21):182-187

[3]Yang Chen (楊晨),Ulrich Gross.Estimation of thermal properties based on inverse heat conduction method [J].Journal of Chemical Industry and Engineering(China) (化工學報),2005,56 (12):2415-2420

[4]Zhang Liwei (張立衛(wèi)),Shan Feng (單鋒).Optimization Method (最優(yōu)化方法) [M].Beijing:Science Press,2010:83

[5]Wang Ling (王凌).Intelligent Optimization Algorithm and Its Application (智能優(yōu)化算法及其應用) [M].Beijing:Tsinghua University Press,2001:12

[6]Zhang C,Guo Z,Liu Y,et al.A review on thermal response test of ground-coupled heat pump systems [J].Renewable and Sustainable Energy Reviews,2014,40:851-867

[7]Gustafsson A M,Westerlund L.Heat extraction thermal response test in groundwater-filled borehole heat exchanger—investigation of the borehole thermal resistance [J].Renewable Energy,2011,36:2388-2394

[8]Raymonda J,Therriena R,Gosselinb L.Borehole temperature evolution during thermal response tests [J].Geothermics,2011,40:69-78

[9]Bozzoli F,Pagliarini G,Rainieri S,Schiavi L.Estimation of soil and grout thermal properties through a TSPEP (two-step parameter estimation procedure) applied to TRT (thermal response test) data [J].Energy,2011,36:839-846

[10]Zhang Changxing (張長興),Hu Songtao (胡松濤),Liu Yufeng (劉玉峰),et al.Determination method for rock-soil thermal properties based on system optimization [J].Journal of Zhejiang University:Engineering Science(浙江大學學報:工學版),2012,46 (12):2237-2242

[11]ASHRAE.ASHRAE Handbook:HVAC Applications [M].Atlanta,GA,USA:ASHRAE,2007

[12]Bernier M A.Ground-coupled heat pump system simulation [J].ASHRAE Transaction.,2001,107 (1):605-616

[13]Bernier M A,Pinel P,Labib R,et al.A multiple load aggregation algorithm for annual hourly simulations of GCHP systems [J].HVAC & Research,2004,10 (4):471-487

[14]Zhang Changxing (張長興),Hu Songtao (胡松濤),Li Xuquan (李緒泉).Application of green function method in calculation on heat conduction of vertical U-tubes heat exchanger [J].Acta Energiae Solaris Sinica(太陽能學報),2010,31 (2):158-162

[15]Richard A Beier,Marvin D Smitha,Jeffrey D Spitler.Reference data sets for vertical borehole ground heat exchanger models and thermal response test analysis [J].Geothermics,2011,40 (1):79-85

[16]Mao Jinfeng (茅靳豐),Li Yong (李勇),Zhang Hua (張華),et al.Thermal short-circuiting and its influence on thermal response in borehole heat exchangers [J].CIESC Journal(化工學報),2013,64 (11):4015-4024

[17]Remund C P.Borehole thermal resistance:laboratory and field studies [J].ASHRAE Transaction,1999,105 (2):439-445

[18]Wang Shuangcheng (王雙成).Physical property calculation of liquid saturated water [J].Chemical Engineering Design(化工設計),1999,9 (6):29-30

[19]Zhang Changxing (張長興),Guo Zhanjun (郭占軍),Liu Yufeng (劉玉峰),et al.A fast forecast method for operation characteristics of ground-coupled heat pump system [J].Transactions of the Chinese Society of Agricultural Engineering(農業(yè)工程學報),2012,28 (24):173-178

[20]Wei Lianwei (魏連偉),Han Wenxiu (韓文秀),Zhang Junyan (張俊艷),et al.Hydrogeological parameter identification on the simulated annealing-genetic algorithm [J].Journal of Tianjin University(天津大學學報),2003,36 (5):618-621

[21]Xu Chuanhua (許傳華),Ren Qingwen (任青文),Zheng Zhi (鄭治),et al.Displacement back analysis of rock mechanic parameters of underground grotto of Suofengying Hydraulic Power Plant [J].Chinese Journal of Geotechnical Engineering(巖土工程學報),2006,28 (11):1981-1985

[22]Zhang Bo (張波),Ye Jiawei (葉家瑋),Hu Yucong (胡郁蔥).Application of optimizing the path by simulated annealing [J].Chinese Journal of Highway and Transport(中國公路學報),2004,17 (1):79-81

[23]Metropolis N,Rosenbluth A W,Rosenbluth M N,Teller A H,Teller E.Equation of state calculation by fast computing machines [J].Journal of Chemical Physics,1953,21 (12):1087-1091

[24]Qian Zhi,Fu Chuli,Xiong Xiangtuan.A modified method for a non-standard inverse heat conduction problem [J].Applied Mathematics and Computation,2006,180 (2):453-468

[25]Wagner V,Bayer P,Kübert M,Blum P.Numerical sensitivity study of thermal response tests [J].Renewable Energy,2012,41:245-253