模拟退火算法在岩土热物性参数确定中的应用

张长兴，王德水，刘玉峰，孙始财，彭冬根

（1 山东科技大学山东省土木工程防灾减灾重点实验室，山东 青岛 266590；2 青岛大学基建处，山东 青岛 266071；3 南昌大学建筑工程学院，江西 南昌 330031）

引 言

在土壤源热泵系统设计时，获取建设地点的岩土热物性参数是实施地埋管换热器（borehole heat exchanger，BHE）设计的先决条件。因此，中国2009年修订的《地源热泵系统工程技术规范》（以下简称《规范》）对热物性参数的测试和确定方法进行了明确规定，该测试又称为岩土的热响应试验(thermal response test,TRT)[1]。通过TRT 试验确定的热物性参数主要是岩土的热导率和容积比热容，同时由于施工工艺、填充材料、U 形管的布置形式等客观因素的影响，确定现场地埋管换热器的实际热阻值也是非常必要的[2]。

根据地埋管换热器的传热原理，通过TRT 试验的方法确定岩土热物性参数是典型的热传导反分析问题，该理论由Shumakov 在1957年提出，很多研究者对简单形体的热传导反问题的数值解进行了研究[3]。随着计算机技术的发展，热传导反问题的求解也日趋成熟，其求解过程的本质是最优化问题，它包含两方面的内容：其一是建立数学模型，即用数学语言描述最优化问题，模型中的数学关系式反映了最优化问题所要达到的目标和各种约束条件；其二是数学求解，数学模型建好以后选择合理的最优化方法进行求解。最优化方法主要包括基于迭代计算的数值方法和智能算法两类，数值方法主要有牛顿法、最速下降法、共轭梯度法等[4]，智能算法主要有遗传算法、混沌优化算法、模拟退火算法、混合算法等[5]。本研究通过建立TRT 试验系统模型，利用模拟退火算法(simulated annealing algorithm,SAA)确定变热流条件下岩土的热导率、容积比热容和地埋管换热器的实际热阻，为土壤源热泵系统的设计和应用提供重要的基础数据。

1 岩土热物性参数识别的数学模型

要获得准确的岩土热物性参数，参数识别的优化模型必须可靠、准确。这就对TRT 试验系统的数学模型提出了更高的要求，同时要求实施优化的目标函数对于待识别参数的变化有足够的敏感性，通过目标函数获得最小值来最终确定最优的热物性参数组合。

1.1 TRT 试验系统的数学模型

系统的数学模型是实施最优化技术的关键，其准确性直接影响参数识别的精度。对于TRT 试验系统而言，其主要的3 个组成部分分别为地埋管换热器、电加热器和循环水泵，土壤与循环水通过地埋管换热器完成热交换，电加热器用以控制换热器加热负荷的强度，循环水泵控制换热器热交换的循环水量，如图1所示。其中，地埋管换热器传热过程的数学建模是最复杂的[6]。根据地埋管换热器的物理模型，其传热分析一般是以钻孔壁为界分成两个计算区域，孔内区域按稳态传热计算，孔外区域按非稳态传热计算。

图1 岩土热响应测试系统示意图Fig.1 Schematic diagram of rock-soil thermal response test system

1.1.1 钻孔外传热的数学模型 地埋管换热器钻孔外传热的数学模型主要有解析解模型和数值解模型两类。在进行TRT 试验系统建模时，解析解模型中的无限长线热源模型（line source model，LSM）、有限长线热源模型（finite line source model，FLSM）和柱热源模型（cylinder source model，CSM）均得到了广泛的应用[2,6]，同时一维、二维和三维的地埋管换热器数值模型在进行热物性参数确定时亦显示了一定的优势[7-10]。然而，考虑到现场岩土TRT 试验的复杂条件和测试现场的电压波动等因素的影响，测试过程中经常出现《规范》要求的恒热流和定循环水流量难以保证的情况，如重新测试，需等待10～14 d 的时间[11]，会对测试的成本和工程进度造成影响，这对地埋管换热器的数学模型提出了更高的要求。针对解析解模型在变热流条件下计算速度慢[12-13]而数值模型难以准确确定换热器热阻的缺点[10]，本研究以应用格林函数法求解的柱热源模型作为钻孔外传热的数学模型，在变热流条件下应用SAA 优化算法实施热物性参数的识别。

在变热流条件下，由于柱热源模型在进行Laplace 变换后的解析解在确定钻孔壁温时需实施负荷叠加才能确定即时的钻孔壁温，会影响热物性参数识别的计算速度。本研究应用格林函数变换求解钻孔壁温的方法，通过计算方式的改变提高参数识别的速度[14]。通过引入格林函数，在数学计算方式的处理上可将积分计算进行离散求和的转换，得到钻孔壁温的计算公式。

在时刻：

在+ 时刻：

式中，=αt/，Ts为土壤的原始温度，Jn和Yn分别为第一类和第二类贝塞尔函数。

从公式可以看出，该卷积算法与传统的G函数法[12]相比不再需要反复调用动态热流和其对应的G函数，极大地缩短了计算时间，因此在土壤源热泵系统长期动态运行特性的预测方面显示了一定的 优势。

1.1.2 钻孔内传热的数学模型 在地埋管换热器的内部传热主要是U 形管内的水与钻孔壁之间的传热，由于钻孔内的热阻与钻孔外比较相对较小，在定热流TRT 试验中经常作为稳态传热处理，然而这种热阻的差异会造成在热物性测试的初始阶段换热器平均水温实测值与模拟值的差异，因此在实施系统优化时通常舍去10 h 内的数据[2,15]，以保证辨识结果的准确性。

在变热流TRT 试验条件下，辨识土壤热物性参数时，钻孔内的热阻是否可作为稳态传热处理需要考虑热流变化的时间间隔、热流的变化幅度、填充材料的热导率和体积比热容等因素的综合影响。确定钻孔内稳态热阻时，利用形状因子法计算钻孔内填充材料的热阻时综合考虑U形管间热短路现象的影响[16]，同时结合U 形管内的对流换热、管壁的导热这两部分热阻，钻孔内稳态热阻可以表示为[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)可以看出，实际工程中确定钻孔内稳态热阻时，填充材料的热导率会受到回填密实程度、地下水和土壤含湿量等因素的影响而造成一定的计算误差，为了降低这一误差，提高管内流体与管内壁之间的对流传热系数的计算精度是非常必要的。因此，在确定对流传热系数时，本研究采用循环水的热物性随水温变化的动态参数，根据文献[18]，循环水的密度、运动黏滞系数和热导率随水温变化的关系式如下

根据地埋管换热器传热原理和能量守恒定律，地埋管换热器在时刻的出口水温可以表示为[19]

1.2 优化模型的建立

根据岩土TRT 试验系统的数学模型，将试验中的逐时加热量和循环水流量作为输入条件，结合地埋管换热器的几何条件和物性参数可计算出地埋管换热器的逐时进出水温度，通过对比地埋管换热器的进出水温度平均值与实测值对式(1)中的岩土热导率和容积比热容进行辨识，并根据式(3)确定地埋管换热器的实际热阻值。本研究采用进出水温度平均值与实测值的均方根误差（root-mean-square error，RMSE）作为优化模型的目标函数

通过对地埋管换热器传热过程的分析可以看出，地埋管平均出水温度的计算是非线性的，会增加参数识别工作的难度。因此，本研究利用SAA算法实施热物性参数的识别，以达到在变热流工况下快速、方便和有效地确定岩土的热导率、容积比热容和地埋管换热器热阻的目的。

2 热物性参数识别的SAA 算法设计

模拟退火算法（simulated annealing algorithm，SAA）是基于Mente Carlo 迭代求解策略的一种随 机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。SAA 算法在某一温度下，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部优解能概率性地跳出并最终趋于全局最优。目前，SAA 算法已经广泛应用在工农业生产中[20-22]。该算法用Metropolis 算法[23]产生组合优化问题解的序列，并由与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 试验概况

为了准确确定岩土热物性参数，在青岛市经济技术开发区进行了TRT 试验，实测岩土的原始温度Ts为15.55℃，测试的地埋管换热器为单U 形PE 管形式，回填材料为10%膨润土、90% SiO2砂子的混合物，测孔的几何尺寸和相关材料的热物性参数见表1。换热器进出水温度的测试时间间隔为5 min，测试时间共计55 h。由于现场用电条件的限制，测试中电压的波动偏差高于5%，造成实测的循环水泵流量和电加热功率出现了不同程度的波动（图3），实测平均加热功率为7460 W，标准差为228.3 W，不符合文献[1,11]中的恒热流加热条件。在变热流TRT 测试中，电加热功率的最大偏差仅为8%，而且间隔时间较短，运行10 h 后将钻孔热阻作为稳态热阻处理。由于实测的加热功率和循环水泵流量均随电压的波动发生变化，必须应用实测的动态电加热功率和循环水流量才能准确确定岩土的热物性参数和换热器的热阻。

图3 现场TRT 试验中的加热功率Fig.3 Dynamic heat power of in-situ thermal response test

3.2 应用SAA 算法确定岩土热物性参数

利用本研究建立的TRT 试验系统模型，结合TRT 试验结果，本研究运用SAA 算法确定了岩土的热导率和容积比热容。数学模型中的电加热器功率、水泵流量和进出水温度均采用试验瞬时值，水的热物性参数采用瞬时温度对应的计算值[式(4)、式(5)和式(6)]，岩土的原始温度、填充材料和PE 管的热导率以及测试孔相关几何尺寸均作为模型的已 知值。

在利用SAA 算法进行热物性参数识别时，由于TRT 试验系统模型中地埋管平均出水温度的计算是非线性的，参数识别过程中一定存在非适定问题[24]，即可能有多个热物性参数组满足目标函数RMSE 最小这一条件。为了获得热物性参数的真实值，本研究首先设定岩土的热导率和容积比热容的识别范围分别为[1,5]和[1000,5000]，对岩土的热导率用0.05 的数据间隔在[1,5]范围内生成81 个值，对容积比热容用50 的数据间隔在[1000,5000]范围内生成81 个值，使两个参数值分别组合成为不同的数组，作为TRT 试验系统模型的输入值计算地埋管换热器的逐时进出口平均水温，结合TRT 试验实测的对应值可计算出6561 个RMSE 值，RMSE 计算值的分布如图4中的等值线所示。可以看出，在RMSE＜0.2 的范围内存在6 处满足RMSE 最小的极值点，为解决参数识别过程中的非适定问题奠定了基础。根据现场TRT 试验钻孔勘察的地质资料，地下岩土以花岗岩、砂岩为主，依照《规范》提供的相关热物性参数，并结合图4的RMSE 值分布状况，在实施 SAA 算法时设定岩土热导率为 2 W·m-1·℃-1、容积比热容为3000 kJ·m-3·℃-1，两个参数的优化范围分别为[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，达到了识别最优热物性参数的目的。

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 最优热物性参数对应的水温模拟值与实测值的 温差平方和动态变化Fig.6 Dynamic variation of square of temperature difference between simulating water temperature and experiment data corresponding to optimal thermal parameters

图7 最优热物性参数对应的水温模拟值与实测值对比Fig.7 Contrast between simulating water temperature and experiment data corresponding to optimal thermal parameters

图8 地埋管换热器热阻随平均水温的动态变化Fig.8 Dynamic variation of borehole heat resistance with mean water temperature

将参数识别结果代入TRT 试验系统的数学模型，由式(7)可计算出各时刻的地埋管换热器进出水温度Tin,sim,i、Tout,sim,i和平均水温Tav,sim,i。图6为最 优热物性参数对应地埋管换热器平均水温的模拟值和实测值温差平方和的逐时变化图，可以看出，在测试10 h 后各数值吻合程度较高，动态变化趋势一致，温差平方和最高值仅为0.15，出现在第41 h，温差平方和的逐时变化也客观上验证了文献[2,15]的正确性。图7为10～55 h 各时刻地埋管换热器平均水温的实测值和模拟值的对比图，可以看出，各时刻两个变量高度吻合，拟合优度R2为0.978。

两个热物性参数确定后，试验系统模型由式(3)可计算出地埋管换热器热阻随地埋管平均水温的动态变化值，如图8所示。随着TRT 试验的进行，地埋管换热器的热阻逐渐趋于稳定，10 h 后热阻的平均值为0.109 m·℃·W-1，换热器的实际热阻值可作为测试地点进行地埋管换热器组群设计的计算 依据。

为了验证SAA 算法参数识别结果的准确性，在同一测试地点按照《规范》的测试方法和数据处理方式确定的两个参数分别为2.41 W·m-1·℃-1和2995 kJ·m-3·℃-1，与《规范》方法确定的两个热物性参数的相对误差为4.1%和1.3%。在变热流TRT 试验系统模型中应用地埋管换热器的柱热源模型，SAA 算法在参数识别过程中具有较强的适用性，而且识别结果准确性高，为地埋管换热器组群的精确设计提供了可靠的基础数据。

4 结 论

本研究利用地埋管换热器的柱热源模型，通过引入格林函数法建立了变热流条件下TRT 试验系统模型，在此基础上通过应用SAA 算法提出了一种岩土热物性参数识别的方法，得到以下结论。

（1）在变热流TRT 试验系统模型中，应用格林函数法的柱热源模型显示了较好的适应性，为岩土热物性参数的准确识别奠定了理论基础。

（2）利用RMSE 值分布图，有效地解决了传热反问题参数识别过程中的非适定问题。在此基础上，结合TRT 试验数据对岩土热物性参数进行识别的过程中SAA 算法对应的退火温度降温速度快，确定的岩土热导率和容积比热容与《规范》方法对应值的相对误差分别为4.1%和1.3%，证明了本研究方法的准确性。

（3）在岩土热物性参数的识别过程中，可计算地埋管换热器热阻随平均水温的动态变化，确定测试现场换热器的有效热阻值，为地埋管换热器组群的精确设计提供有效的基础数据。

符 号 说 明

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(农业工程学报),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 (张立卫),Shan Feng (单锋).Optimization Method (最优化方法) [M].Beijing:Science Press,2010:83

[5]Wang Ling (王凌).Intelligent Optimization Algorithm and Its Application (智能优化算法及其应用) [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(农业工程学报),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