发明内容
针对现有技术的缺陷,本发明的目的在于提供一种可用于地埋管换热器中断热响应测试的参数反演方法及系统。该方法可在供电恢复后,无需等待地下岩土温度恢复至初始水平,直接恢复热响应测试,并能利用中断前后的测试数据,准确地估计土壤及回填材料热物性等多个参数。
本发明的技术方案如下:
用于地埋管换热器中断热响应测试的参数反演方法,其包括:
S1:获得预埋的U型地埋管换热器的已知物理参数;
S2:通过原位热响应测试测得该U型地埋管换热器中循环流体平均温度即流体温度实测值,及其随时间变化的曲线,即循环流体的温度-时间测试曲线,根据对该测试曲线的数值求导,获得温度对时间的数值导数,即温度导数实测值;
S3:确定待估参数的种类和赋予其初始值,将获得的已知物理参数和待估计参数的初始值输入地埋管换热器传热模型中,获得不同传热时间下的循环流体的温度响应值和温度对时间或时间对数的导数值,即流体温度计算值和温度导数计算值;
S4:基于所述流体温度实测值与所述流体温度计算值,和所述温度导数实测值与所述温度导数计算值,并结合正则化处理,构建单目标优化函数;
S5:基于所述单目标优化函数,通过优化算法获得所述待估计参数的最佳估计结果,即所述待估计参数的反演值。
根据本发明的一些具体实施方式,所述原位热响应测试通过移动热响应测试仪实现。
根据本发明的一些具体实施方式,所述单目标优化函数选自所述流体温度实测值和所述流体温度计算值之差的二范数的平方与所述温度导数实测值和所述温度导数计算值之差的二范数的平方及正则化项,即温度项和导数项及正则化项,或所述温度项和导数项的加权组合。
根据本发明的一些具体实施方式,所述加权组合中,温度项与导数项各自的权重通过所述流体温度实测值的不确定度及所述温度导数实测值的不确定度确定。
根据本发明的一些具体实施方式,所述单目标优化函数设置如下:
其中W
1为所述流体温度实测值的不确定度,W
2为所述温度导数实测值的不确定度,Y为所述流体温度实测值,T
f(β)为所述流体温度计算值,
表示基于流体温度测量值使用中心差分法计算得到的所述温度导数实测值,
表示通过对时间对数偏导得到的所述温度导数计算值,t表示时间,||||
2表示二范数。
或,根据本发明的一些具体实施方式,所述单目标优化函数设置如下:
其中,β表示待估计参数的向量,β0表示其初始值,γ表示正则化参数。
根据本发明的一些具体实施方式,所述地埋管换热器传热模型为基于复合介质的无限长线热源模型。
根据本发明的一些具体实施方式,所述地埋管换热器传热模型为如下的复合介质的无限长线热源模型:
φ=akJn(urb)J’n(aurb)-J’n(urb)Jn(aurb)
ψ=akJn(urb)Y’n(aurb)-J’n(urb)Yn(aurb)
f=akYn(urb)J’n(aurb)-Y’n(urb)Jn(aurb)
g=akYn(urb)Y’n(aurb)-Y’n(urb)Jn(aurb)
其中,G(t)表示从U型管外壁至半无限大介质的非稳态传热过程的热响应函数;kb表示回填材料导热系数;ab表示回填材料热扩散系数;u表示积分变量;t表示时间;i表示求和指标;rA=D-ro和rB=D+ro为以钻孔圆心为坐标原点的两个计算点的径向坐标;J2i表示2i阶贝塞尔函数;D表示U型管圆心的径向坐标;k、a为无量纲参数,k=ks/kb,a=(ab/as)1/2,as表示土壤热扩散系数;Jn和Yn分别表示第一类和第二类n阶Bessel函数;Jn’和Yn’分别表示Jn和Yn的一阶导数;rb表示钻孔半径;Tf为循环流体温度;T0为土壤的初始温度;ql表示单位长度换热器的热流密度(W/m);N表示钻孔内U型管数;Rp表示管子热阻;kp是U型管导热系数;ro和ri分别是U型管的内外径;α表示管内对流传热系数。
根据本发明的一些具体实施方式,所述已知物理参数包括:土壤和回填材料的初始温度T0,地埋管钻孔外径rb、钻孔长度L,地埋管外径ro、内径ri、圆心距2D、管壁导热系数kp,循环流体的体积热容Cf、体积流量Vf,加热器平均加热率q。
根据本发明的一些具体实施方式,所述待估计参数包括以下参数中的一种或多种:岩土导热系数ks,岩土体积热容Cs,回填材料导热系数kb,和回填材料体积热容Cb。
根据本发明的一些具体实施方式,所述S3具体包括:
S31:建立基于复合介质的无限长线热源模型,根据该模型获得循环流体的温度对时间的导数的解析解模型;
S32:通过所述已知物理参数和所述待估计参数的初始值对所述解析解模型中的模型参数进行赋值,根据赋值后的解析解模型获得循环流体的温度-时间响应函数,根据该响应函数获得所述流体温度计算值及所述温度导数计算值。
根据本发明的一些具体实施方式,S31中所述解析解模型包括中断热响应测试前的解析解模型和中断热响应测试后继续测试得到的解析解模型,其中:
所述中断热响应测试前的解析解模型如下:
所述中断热响应测试后继续测试得到的解析解模型如下:
G2(t)=G(t)-G(t-t1)+G(t-t2)
φ=akJn(urb)J’n(aurb)-J’n(urb)Jn(aurb)
ψ=akJn(urb)Y’n(aurb)-J’n(urb)Yn(aurb)
f=akYn(urb)J’n(aurb)-Y’n(urb)Jn(aurb)
g=akYn(urb)Y’n(aurb)-Y’n(urb)Jn(aurb);
根据本发明的一些具体实施方式,所述S2具体包括:
S21根据所述原位热响应测试中获得的不同测试时间下循环流体的温度数据,采用一定的时间间隔,通过多种不同的数值求导公式获得测试时间间隔内循环流体温度对时间的导数,即所述温度导数实测值,并计算各公式获得的温度导数实测值的不确定度。
根据本发明的一些具体实施方式,所述S21具体包括:
S211:根据测试所得的循环流体温度响应的时变特性与不确定度,选取一定的时间间隔用于数值求导;
S212:在选定的时间间隔下,采用3种或以上的具有相同误差精度的数值求导公式,进行温度对时间的求导,获得的数值导数即所述温度导数实测值;
S213:计算不同求导公式获得的数值导数的期望、方差,根据所得期望和方差获得其不确定度,即所述温度导数实测值的不确定度。
本发明进一步提供了一种地埋管换热器多参数协同反演的正则化系统,其包括:一个或多个处理器及用于存储一个或多个程序的存储器,其中,当所述一个或多个程序被所述一个或多个处理器执行时,所述一个或多个处理器可执行上述任一参数反演方法。
本发明具有以下有益效果:
(1)本发明创新性地将地埋管换热器传热模型、温度导数模型与正则化有机结合,能有效解决地埋管换热器多参数反演的不适定性问题,其不仅能充分利用早期高敏感性的热响应测试数据,当被中断的热响应测试的电源恢复后,可直接继续进行测试而无需等待,提高了现场热响应测试的效率和参数估计的可靠性。
(2)在一些具体实施方式中,本发明进一步将温度导数加入到复合介质的无限长线热源模型和Tikhonov正则化的反演算法中来解决反演问题的不适定性,在缩短现场热响应测试时间的同时得到可靠的、稳定的结果。
(3)在一些具体实施方式中,本发明利用短时模型使短时高敏感性数据能够用于参数反演,且正则化方法可以抑制数据噪声的影响,从而获得平滑稳定的近似解。
具体实施方式
以下结合实施例和附图对本发明进行详细描述,但需要理解的是,所述实施例和附图仅用于对本发明进行示例性的描述,而并不能对本发明的保护范围构成任何限制。所有包含在本发明的发明宗旨范围内的合理的变换和组合均落入本发明的保护范围。
实施例1
通过以下步骤进行本发明的用于中断热响应测试的地埋管换热器参数反演的方法:
S1:获取U型地埋管换热器原位热响应测试得到的不同测试时间下的循环流体温度数据及U型地埋管换热器的已知物理参数;其中,已知物理参数包括:土壤和回填材料的初始温度T0、钻孔外径rb、钻孔长度L、U型管外径ro、U型管内径ri、U形管管脚圆心距的一半D、U型管壁导热系数kp、循环流体的体积热容Cf、循环流体的体积流量Vf、电加热器平均加热率q。
S2:获得流体温度实测值及其不确定度、温度导数实测值及其不确定度。
其可具体包括:
S21:根据U型地埋管换热器原位热响应测试中获得的循环流体温度数据,采用一定的测试时间段,通过多种不同的数值求导公式计算出测试时间段内获得的循环流体温度对时间的导数及各种方法下获得的导数的不确定度,即温度导数实测值及其不确定度,并计算该测试时间段内的循环流体的平均温度和不确定度,即流体温度实测值及其不确定度,其可进一步包括:
S211:根据测试所得的循环流体温度响应的时变特性与不确定度,选取合适的时间间隔用于数值求导;
S212:在选定的时间间隔下,采用3种或以上的具有相同误差精度的数值求导公式,对测试时间段内获得的温度数据进行对时间的求导,获得数值导数,即所述温度导数实测值;
S213:计算不同求导公式获得的数值导数的期望、方差,根据所得期望和方差获得温度导数实测值的不确定度;
S214:计算测试时间段内的循环流体的平均温度和不确定度,获得流体温度实测值及其不确定度。
S3将获得的已知物理参数和待估计的物理参数的初始值输入U型地埋管换热器传热模型中,获得不同传热时间下的循环流体的温度响应值和温度对时间的导数值,即流体温度计算值和温度导数计算值;其中,所述U型地埋管换热器传热模型为线热源型传热模型,所述待估计的物理参数包括岩土导热系数ks、回填材料导热系数kb、岩土体积热容Cs和回填材料的体积热容Cb中的一种或多种。
其可具体包括:
S31:建立基于循环流体温度的复合介质无限长线热源模型,并根据该模型得到循环流体平均温度对时间的导数的解析解模型;
具体的,该复合介质无限长线热源模型可建立如下:
其中:
其中,G表示从U型管外壁至半无限大介质的非稳态传热过程,其物理意义为非稳态热阻;kb表示回填材料导热系数;ab表示回填材料热扩散系数,rA=D-ro和rB=D+ro为A、B两点的径向坐标(以钻孔圆心为坐标原点);k、a为无量纲参数,k=ks/kb,a=(ab/as)1/2;Jn和Yn分别表示第一类和第二类n阶Bessel函数;Jn’和Yn’分别表示Jn和Yn的一阶导数;rb表示钻孔半径;Tf为循环流体温度;T0为土壤的初始温度;ql表示单位长度换热器的热流密度(W/m);N表示钻孔内U型管数(本实施例为单U管,N=2);Rp表示管子热阻;kp是U型管导热系数;ro和ri分别是U型管的内外径;α表示管内对流传热系数。
其解析解模型包括:
中断热响应测试前的所述解析解模型:
中断热响应测试后继续测试的解析解模型:
G2(t)=G(t)-G(t-t1)+G(t-t2)
S32:通过已知物理参数和待估计的物理参数的初始值分别对所得循环流体平均温度的解析解模型和所得循环流体平均温度对时间的导数的解析解模型中的模型参数进行赋值,获得赋值后的循环流体温度-时间响应函数和循环流体温度导数-时间响应函数。根据该响应函数获得不同时间下的温度响应值和温度对时间的导数值,即所述流体温度计算值和温度导数计算值。
S4:构建单目标函数。
其中,单目标函数可选择为以下函数中的任一种:流体温度实测值和流体温度计算值之差的二范数的平方及正则项;温度导数实测值和温度导数计算值之差的二范数的平方及正则项;流体温度实测值和所述流体温度计算值之差的二范数的平方与所述温度导数实测值和温度导数计算值之差的二范数的平方的加权组合及正则项。
优选的,其中,加权组合中温度项(即流体温度实测值和流体温度计算值之差的二范数的平方)和导数项(即温度导数实测值和温度导数计算值之差的二范数的平方)各自的权重而根据温度导数实测值的不确定度及流体温度实测值的不确定度确定,如
其中W
1为流体温度测量值的不确定度,W
2为温度导数测量值的不确定度。
S5:通过优化算法,调整待估计物理参数,使S4得到的目标函数最小化,目标函数取最小值时所得的待估计物理参数的值即为对应参数的反演值。
运用优化算法实现地埋管换热器的参数反演,获得其最优的参数估计结果。
实施例2
采用国际标准沙箱实验(Beier,2011)数据对本发明的反演方法进行验证,通过该标准沙箱试验进行的地埋管换热器热响应测试中,加热器和循环水泵的电力中断发生在9至11小时之间,共两个小时,其电功图如附图2所示。该测试中,由于水泵在中断期间是关闭的,因此无法获得有代表性的循环温度数据,但在整个测试期间,包括中断期间,钻孔温度和土壤温度均可以得到,与中断的现场试验的状态相符。
在以上测试情况下,参照附图1及图3~5,采用以下程序过程实现本实施例的参数反演方法:
参照附图1,在主程序下运行,并进行:
步骤1:输入收集的热响应测试数据,其中,热响应测试数据可包括中断前后的两次加热采样的测试数据;其输入方式可如:对第一次加热采样数据,每隔15分钟取一个测试获得的循环流体温度数据输入,对第二次加热采样数据,按对数距离取60个循环流体温度数据输入;
步骤2:根据输入的热响应测试数据,计算循环流体温度的数值导数,即获得所述温度导数实测值;
步骤3:确定需要反演的参数即待估参数的种类和数量;
步骤4:输入地埋管换热器传热模型的模型条件及收集的岩土热物性数据;
步骤5:根据输入的模型条件及岩土热物性数据初始化地埋管换热器传热模型及土壤的相关参数;
步骤6:给定待估参数的初始值;
步骤7:调用循环流体温度计算子程序,计算得到该时间点的循环流体温度;
步骤8:调用循环流体温度对时间的导数计算子程序,计算得到该时间点的循环流体温度的导数;
步骤9:调用参数反演迭代子程序,进行目标函数的最优化;
其中,如图3所示,步骤7中循环流体温度计算子程序的运行过程包括:
步骤701:输入地埋管换热器及岩土参数。其中,地埋管换热器参数包括:钻孔半径、钻孔高度,回填材料导热系数及热扩散系数,U型管导热系数、U型管间距的一半、内径、外径和类型。岩土参数包括:土壤的导热系数、热扩散系数和初始温度;
步骤702:设置参数的初始值,给定参数的上下限,进行初始化;
步骤703:输入传热模型的传热函数即G函数,本实施例采用的传热模型为复合介质无限长线热源模型,这一模型可以有效利用热响应测试的短时间数据,从而达到更好的拟合结果,是目前充分利用热响应测试数据理论上最成熟的一种模型,具体可通过下式得到:
其中:
其中,G表示从U型管外壁至半无限大介质的非稳态传热过程,其物理意义为非稳态热阻;kb表示回填材料导热系数;ab表示回填材料热扩散系数;rA=D-ro和rB=D+ro为A、B两点的径向坐标(以钻孔圆心为坐标原点);k、a为无量纲参数,k=ks/kb,a=(ab/as)1/2,as表示土壤热扩散系数;Jn和Yn分别表示第一类和第二类n阶Bessel函数;Jn’和Yn’分别表示Jn和Yn的一阶导数;rb表示钻孔半径;Tf为循环流体温度;T0为土壤的初始温度;ql表示单位长度换热器的热流密度(W/m);N表示钻孔内U型管数(本实施例为单U管,N=2);Rp表示管子热阻;kp是U型管导热系数;ro和ri分别是U型管的内外径;α表示管内对流传热系数;
步骤704:输入第一次加热采样时间点向量,通过模型计算得到第一次加热采样的循环流体温度;
步骤705:计算第二次加热采样的循环流体的平均温度,其计算模型如下:
设停止加热的时刻为t1,重新开始加热的时刻为t2,根据Duhamel定理,重新开始加热后,有:
G2(t)=G(t)-G(t-t1)+G(t-t2)
步骤706:输入第二次加热采样时间点向量,计算得到第二次加热采样的循环流体温度。
其中,如图4所示,步骤8中所述循环流体温度对时间的导数计算子程序运行如下:
步骤801:根据复合介质无限长线热源模型得到第一次加热采样的循环流体温度对时间的导数的计算模型,如下:
其中
步骤802:输入第一次加热采样时间点向量,根据第一次加热采样的循环流体温度对时间的导数的计算模型计算得到第一次加热采样的循环流体的温度对时间的导数。
步骤803:根据复合介质无限长线热源模型得到第二次加热采样的循环流体温度对时间的导数的计算模型,如下:
其中
步骤804:输入第二次加热采样时间点向量,根据第二次加热采样的循环流体温度对时间的导数的计算模型计算得到第二次加热采样的循环流体温度导数。
其中,如图5所示,步骤9中参数反演迭代子程序的运行包括:
步骤901:确定目标函数,正则化方法的目标函数可认为是修正的最小二乘范数,它是最小二乘范数项与正则项的和,本实施例采用了Tikhonov正则化项,其形式如下:
其中,T
f(β)分别用于表示m个观测时间的模型计算的流体温度,
是用中心差分计算得到的流体温度的数值导数,
表示模型计算的流体温度对对数时间的偏导数,β表示待估计参数的向量,例如,β
T=[k
s,a
s,k
b,a
b];采用的初始值是
y=[Y
1,Y
2,...,Y
m]
T代表在m个观测时间内循环流体平均温度的观测向量;γ是正则化参数,用来控制解的光滑度。
步骤902:输入γ值范围,γ=0时表示正则项为零,此时为普通最小二乘法目标函数;
步骤903:令i为正则化参数的顺序数,初始化循环指针i=1,给定被拟合参数的初始值及上下限值;
步骤904:利用MATLAB系统内置程序非线性拟合命令来实现对目标函数的最小化;
步骤905:令i=i+1;
步骤906:判断i≤M,如果是,则转入步骤904;如果否,则算法结束,其中,M为正则化参数γ的个数。
本实施例中,热响应测试的各项实验参数参考值如表1所示:
表1沙箱热响应测试实验参数
| 参数 |
描述 |
值 |
精度 |
| T0 |
土壤与回填材料初始温度 |
22.06℃ |
±0.03 |
| rb |
钻孔外径 |
6.5cm |
±0.5 |
| L |
钻孔长度 |
18.3m |
±0.5 |
| ro |
U型管外径 |
1.6700cm |
±0.0005 |
| ri |
U型管内径 |
1.3665cm |
±0.0005 |
| D |
U型管圆心距的一半 |
2.6cm |
±0.5 |
| kp |
U型管壁导热系数 |
0.39Wm-1K-1 |
±0.05 |
| ks |
湿沙导热系数 |
2.82Wm-1K-1 |
±0.1 |
| kb |
回填材料导热系数 |
0.73Wm-1K-1 |
±0.04 |
| Cf |
循环流体的体积流量 |
0.197Ls-1 |
±0.004 |
| q |
电加热器平均加热率 |
1056W |
±4 |
| Cs |
湿沙的容积热容估计值 |
3.2E+6Jm-3K-1 |
NA |
| Cb |
回填材料的容积热容估计值 |
3.8E+6Jm-3K-1 |
NA |
a湿沙和回填材料的热容数据为估计值(非独立测量值),这可能导致热扩散系数的估计值存在不确定性。
本实施例要估计的4个参数为:ks、kb、as、ab,其参考值分别为[ks,as,kb,ab]=[2.82,8.8e-7,0.73,1.9e-7]。
根据以上程序过程,表2示出了使用不同的参数初始值和不同的正则化参数的计算结果,如下:
表2使用不同的初始值和不同的正则化参数的估计结果
可以看出,在设置的两组初始值下,当γ=0.7时计算出的结果最好。
其中,对于第一组结果,土壤导热系数ks的相对误差为0.35%,土壤热扩散系数as的相对误差为25.0%,回填材料导热系数kb的相对误差为14.1%,回填材料热扩散系数ab的相对误差为10.5%。
第二组结果表明如果初始值在一定的范围内,正则化参数γ
时还是会有最佳的估计值。
经过计算,对于[k
s,a
s,k
b,a
b]的初始值,在上界为[2.2,4e-7,1,5e-7],下界为[2,1e-7,0.5,1e-7]的范围内,当γ
时会有最佳估计结果。如果初值的选取超过这个范围,则不一定在γ=0.7时得到最佳估计结果,无法确定合适的γ值。
进一步的,图6示出了本实施例使用正则化的四参数估计的结果。从图中可以看出,模型的计算值和沙箱数据测量值之间吻合的非常好,参数估计结果的精度也在可接受的范围内。
图7示出了实施例中测量的循环流体温度数据,可以看到,在画圈部分的数据有轻微的震荡。这部分震荡虽然对测量温度值影响很小,但会让计算出的数值导数产生很大的震荡。如果剔除掉这部分数据,可得到更好的结果。因此,可将整个数据分割成三段。
在三段数据中,如果间隔取得太小,数值导数就会产生很大的震荡,如果间隔取得太大,虽然可以得到更加平滑的导数曲线,但会失掉很多信息,难以充分利用测量数据。因此,为了让早期高敏感数据发挥更大的作用,本实施例进一步取第一段时间间隔为5分钟,第二段时间间隔为15分钟,第三段仍然依照对数等距间隔取60个数据,由此可得到表3及附图8所示的去除掉了震荡点之后的参数估计结果,如下:
表3去除震荡点后的估计结果
可以看出,以上计算结果中,第一组获得土壤导热系数ks的相对误差为7.1%,土壤热扩散系数as的相对误差为20.0%,回填材料导热系数kb的相对误差为17.8%,回填材料热扩散系数ab的相对误差为3.1%。其相比于不去除震荡点的计算结果更佳。
以上实施例仅是本发明的优选实施方式,本发明的保护范围并不仅局限于上述实施例。凡属于本发明思路下的技术方案均属于本发明的保护范围。应该指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下的改进和润饰,这些改进和润饰也应视为本发明的保护范围。