CN116401883B - 用于地埋管换热器中断热响应测试的参数反演方法及系统 - Google Patents
用于地埋管换热器中断热响应测试的参数反演方法及系统 Download PDFInfo
- Publication number
- CN116401883B CN116401883B CN202310396978.3A CN202310396978A CN116401883B CN 116401883 B CN116401883 B CN 116401883B CN 202310396978 A CN202310396978 A CN 202310396978A CN 116401883 B CN116401883 B CN 116401883B
- Authority
- CN
- China
- Prior art keywords
- temperature
- derivative
- value
- representing
- heat exchanger
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000012360 testing method Methods 0.000 title claims abstract description 90
- 230000004044 response Effects 0.000 title claims abstract description 62
- 238000000034 method Methods 0.000 title claims abstract description 42
- 238000011065 in-situ storage Methods 0.000 claims abstract description 9
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 8
- 239000012530 fluid Substances 0.000 claims description 95
- 239000002689 soil Substances 0.000 claims description 37
- 230000006870 function Effects 0.000 claims description 30
- 238000004364 calculation method Methods 0.000 claims description 29
- 239000000463 material Substances 0.000 claims description 29
- 238000012546 transfer Methods 0.000 claims description 25
- 238000010438 heat treatment Methods 0.000 claims description 23
- 238000005259 measurement Methods 0.000 claims description 17
- 238000005553 drilling Methods 0.000 claims description 15
- 238000005457 optimization Methods 0.000 claims description 14
- 239000002131 composite material Substances 0.000 claims description 11
- 238000009795 derivation Methods 0.000 claims description 10
- 238000005316 response function Methods 0.000 claims description 10
- 239000011435 rock Substances 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 6
- 230000004907 flux Effects 0.000 claims description 5
- 238000005070 sampling Methods 0.000 description 12
- 230000010355 oscillation Effects 0.000 description 8
- 238000010586 diagram Methods 0.000 description 6
- 244000035744 Hura crepitans Species 0.000 description 3
- 230000000704 physical effect Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000009529 body temperature measurement Methods 0.000 description 2
- 238000009792 diffusion process Methods 0.000 description 2
- 239000004576 sand Substances 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 241000508725 Elymus repens Species 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000004134 energy conservation Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000004904 shortening Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/17—Mechanical parametric or variational design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E10/00—Energy generation through renewable energy sources
- Y02E10/10—Geothermal energy
Abstract
本发明公开了一种用于地埋管换热器中断热响应测试的参数反演方法及系统。所述反演方法可使用结合了温度响应、温度响应导数和正则化的反演算法，其可在现场热响应测试被迫中断其后恢复供电时，直接继续进行新的热响应测试而无需重新测试。该方法将温度导数用于参数估计，可充分利用测试数据，同时通过正则化可有效抑制数据噪声的影响，从而获得平滑稳定的反演解，不仅能节约大量测试时间与成本，还可以得到可靠、稳定的结果。
Description
技术领域
本发明属于地源热泵技术应用与节能的技术领域，具体涉及获得中断热响应测试的地埋管换热器参数方法的技术领域。
背景技术
地埋管换热器的现场热响应测试是获得当地岩土表观热物性的主要方法，是大规模地源热泵系统设计前必须开展的工作。热响应测试通常要求测试时长≥48小时。原位测试的现场情况复杂，经常发生测试期间电源意外中断而导致测试停止的事故。
若测试中途中断，有两种潜在的应对方案：一是等待土壤温度恢复到测试前的初始温度后，重新进行热响应测试。由于地下岩土温度恢复过程缓慢，重新测试将严重影响了现场施工进度。二是停止测试，仅使用断电前获取的数据估计地下岩土热物性。通过热响应测试计算岩土热物性是一种参数估计问题，属于典型的导热反问题，存在求解过程不适定的问题，即参数反演解的存在性、唯一性和稳定性的问题。热响应测试中断导致的测试数据减少、测试时间缩短，使得反演不适定的问题更加严重，国内尚无合适的参数估计算法解决这一问题。
发明内容
针对现有技术的缺陷，本发明的目的在于提供一种可用于地埋管换热器中断热响应测试的参数反演方法及系统。该方法可在供电恢复后，无需等待地下岩土温度恢复至初始水平，直接恢复热响应测试，并能利用中断前后的测试数据，准确地估计土壤及回填材料热物性等多个参数。
本发明的技术方案如下：
用于地埋管换热器中断热响应测试的参数反演方法，其包括：
S1：获得预埋的U型地埋管换热器的已知物理参数；
S2：通过原位热响应测试测得该U型地埋管换热器中循环流体平均温度即流体温度实测值，及其随时间变化的曲线，即循环流体的温度-时间测试曲线，根据对该测试曲线的数值求导，获得温度对时间的数值导数，即温度导数实测值；
S3：确定待估参数的种类和赋予其初始值，将获得的已知物理参数和待估计参数的初始值输入地埋管换热器传热模型中，获得不同传热时间下的循环流体的温度响应值和温度对时间或时间对数的导数值，即流体温度计算值和温度导数计算值；
S4：基于所述流体温度实测值与所述流体温度计算值，和所述温度导数实测值与所述温度导数计算值，并结合正则化处理，构建单目标优化函数；
S5：基于所述单目标优化函数，通过优化算法获得所述待估计参数的最佳估计结果，即所述待估计参数的反演值。
根据本发明的一些具体实施方式，所述原位热响应测试通过移动热响应测试仪实现。
根据本发明的一些具体实施方式，所述单目标优化函数选自所述流体温度实测值和所述流体温度计算值之差的二范数的平方与所述温度导数实测值和所述温度导数计算值之差的二范数的平方及正则化项，即温度项和导数项及正则化项，或所述温度项和导数项的加权组合。
根据本发明的一些具体实施方式，所述加权组合中，温度项与导数项各自的权重通过所述流体温度实测值的不确定度及所述温度导数实测值的不确定度确定。
根据本发明的一些具体实施方式，所述单目标优化函数设置如下：
其中W1为所述流体温度实测值的不确定度，W2为所述温度导数实测值的不确定度，Y为所述流体温度实测值，Tf(β)为所述流体温度计算值，表示基于流体温度测量值使用中心差分法计算得到的所述温度导数实测值，/>表示通过对时间对数偏导得到的所述温度导数计算值，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为实施例中本发明所述方法的主程序框图。
图2为实施例中所述沙箱实验电功图。
图3为实施例中所述循环流体温度计算框图。
图4为实施例中所述循环流体温度导数计算框图。
图5为实施例中所述参数估计迭代计算框图。
图6为实施例中所述使用正则化的四个参数估计的结果。
图7为实施例中所述测量获得的循环流体温度数据。
图8为实施例中所述去除震荡点后的参数估计结果。
具体实施方式
以下结合实施例和附图对本发明进行详细描述，但需要理解的是，所述实施例和附图仅用于对本发明进行示例性的描述，而并不能对本发明的保护范围构成任何限制。所有包含在本发明的发明宗旨范围内的合理的变换和组合均落入本发明的保护范围。
实施例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：构建单目标函数。
其中，单目标函数可选择为以下函数中的任一种：流体温度实测值和流体温度计算值之差的二范数的平方及正则项；温度导数实测值和温度导数计算值之差的二范数的平方及正则项；流体温度实测值和所述流体温度计算值之差的二范数的平方与所述温度导数实测值和温度导数计算值之差的二范数的平方的加权组合及正则项。
优选的，其中，加权组合中温度项(即流体温度实测值和流体温度计算值之差的二范数的平方)和导数项(即温度导数实测值和温度导数计算值之差的二范数的平方)各自的权重而根据温度导数实测值的不确定度及流体温度实测值的不确定度确定，如其中W1为流体温度测量值的不确定度，W2为温度导数测量值的不确定度。
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正则化项，其形式如下：
其中，Tf(β)分别用于表示m个观测时间的模型计算的流体温度，是用中心差分计算得到的流体温度的数值导数，/>表示模型计算的流体温度对对数时间的偏导数，β表示待估计参数的向量，例如，βT＝[ks，as，kb，ab]；采用的初始值是y＝[Y1，Y2，...，Ym]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％。
第二组结果表明如果初始值在一定的范围内，正则化参数γ时还是会有最佳的估计值。
经过计算，对于[ks，as，kb，ab]的初始值，在上界为[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％。其相比于不去除震荡点的计算结果更佳。
以上实施例仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例。凡属于本发明思路下的技术方案均属于本发明的保护范围。应该指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下的改进和润饰，这些改进和润饰也应视为本发明的保护范围。
Claims (8)
1.用于地埋管换热器中断热响应测试的参数反演方法，其特征在于，其包括：
S1：获得预埋的U型地埋管换热器的已知物理参数；
S2：通过原位热响应测试测得该U型地埋管换热器中循环流体平均温度并作为流体温度的实测值，及其随时间变化的曲线，即循环流体的温度-时间测试曲线，根据对该测试曲线进行数值求导，获得温度对时间的数值导数，即温度导数实测值；
S3：确定待估参数的种类并赋予其初始值，将获得的已知物理参数和待估计参数的初始值输入地埋管换热器传热模型中，获得不同测试时间下的循环流体的温度响应值和温度对时间或时间对数的导数值，即流体温度计算值和温度导数计算值；
S4：基于所述流体温度实测值与所述流体温度计算值，和所述温度导数实测值与所述温度导数计算值，并结合正则化处理，构建单目标优化函数；
S5：基于所述单目标优化函数，通过优化算法获得所述待估计参数的最佳估计结果，即所述待估计参数的反演值；
其中，所述S3具体包括：
S31：建立基于复合介质的无限长线热源模型，根据该模型获得循环流体的温度对时间的导数的解析解模型；
S32：通过所述已知物理参数和所述待估计参数的初始值对所述解析解模型中的模型参数进行赋值，根据赋值后的解析解模型获得循环流体的温度-时间响应函数，根据该响应函数获得所述流体温度计算值及所述温度导数计算值；
其中，所述解析解模型包括中断热响应测试前的解析解模型和中断热响应测试后继续测试的解析解模型，其中：
所述中断热响应测试前的解析解模型如下：
所述中断热响应测试后继续测试得到的解析解模型如下：
其中，G(t)表示从U型管外壁至半无限大介质的非稳态传热过程的热响应函数；kb表示回填材料导热系数；ab表示回填材料热扩散系数；u表示积分变量；t表示时间；i表示求和指标；rA＝D－ro和rB＝D+ro为以钻孔圆心为坐标原点的A、B两点的径向坐标；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型管的内外径；α表示管内对流传热系数；。
2.根据权利要求1所述的用于地埋管换热器中断热响应测试的参数反演方法，其特征在于，所述单目标优化函数选自所述流体温度实测值和所述流体温度计算值之差的二范数的平方与所述温度导数实测值和所述温度导数计算值之差的二范数的平方及正则化项，即温度项和导数项及正则化项，或所述温度项和导数项的加权组合。
3.根据权利要求2所述的用于地埋管换热器中断热响应测试的参数反演方法，其特征在于，所述加权组合中，温度项与导数项各自的权重通过所述流体温度实测值的不确定度及所述温度导数实测值的不确定度确定。
4.根据权利要求3所述的用于地埋管换热器中断热响应测试的参数反演方法，其特征在于，所述单目标优化函数设置如下：
其中W1为所述流体温度实测值的不确定度，W2为所述温度导数实测值的不确定度，Y为所述流体温度实测值，Tf(β)为所述流体温度计算值，表示基于温度实测值Y使用中心差分法计算得到的所述温度导数实测值，/>表示通过对时间对数偏导得到的所述温度导数计算值，t表示时间，||||2表示二范数；
或所述单目标优化函数设置如下：
其中，β表示待估计参数的向量，β0表示其初始值，γ表示正则化参数。
5.根据权利要求1所述的用于地埋管换热器中断热响应测试的参数反演方法，其特征在于，所述地埋管换热器传热模型为基于复合介质的无限长线热源模型，如下：
φ＝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为以钻孔圆心为坐标原点的计算点的径向坐标；D表示U型管圆心的径向坐标；J2i表示2i阶贝塞尔函数；k、a为无量纲参数，k＝ks/kb，a＝(ab/as)1/2，as表示土壤热扩散系数；Jn和Yn分别表示第一类和第二类n阶Bessel函数；Jn’和Yn’分别表示Jn和Yn的一阶导数；rb表示钻孔半径；Tf为循环流体温度；T0为土壤的初始温度；ql表示单位长度换热器的热流密度；N表示钻孔内U型管数；Rp表示管子热阻；kp是U型管导热系数；ro和ri分别是U型管的内外径；α表示管内对流传热系数。
6.根据权利要求1所述的用于地埋管换热器中断热响应测试的参数反演方法，其特征在于，所述已知物理参数包括：土壤和回填材料的初始温度T0，地埋管钻孔外径rb、钻孔长度L，地埋管外径ro、内径ri、U形管间距2D、管壁导热系数kp，循环流体的体积热容Cf、体积流量Vf，加热器平均加热率q；所述待估计参数包括以下参数中的一种或多种：岩土导热系数ks，岩土体积热容Cs，回填材料导热系数kb，和回填材料体积热容Cb。
7.根据权利要求1所述的用于地埋管换热器中断热响应测试的参数反演方法，其特征在于，所述S2具体包括：
根据所述原位热响应测试得到的循环流体温度响应的时变特性与不确定度，选取一定的时间间隔用于数值求导；
在选定的时间间隔下，采用3种或以上的具有相同误差精度的数值求导公式，进行温度对时间的求导，获得的数值导数即所述温度导数实测值；
计算不同求导公式获得的数值导数的期望、方差，根据所得期望和方差获得其不确定度，即温度导数实测值的不确定度。
8.一种地埋管换热器多参数协同反演的正则化系统，其特征在于，所述系统包括：移动热响应测试仪、一个或多个处理器及用于存储一个或多个程序的存储器，其中，当所述一个或多个程序被所述一个或多个处理器执行时，所述一个或多个处理器可执行如权利要求1-7中任一所述的参数反演方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310396978.3A CN116401883B (zh) | 2023-04-14 | 2023-04-14 | 用于地埋管换热器中断热响应测试的参数反演方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310396978.3A CN116401883B (zh) | 2023-04-14 | 2023-04-14 | 用于地埋管换热器中断热响应测试的参数反演方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116401883A CN116401883A (zh) | 2023-07-07 |
| CN116401883B true CN116401883B (zh) | 2023-09-26 |
Family
ID=87007224
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310396978.3A Active CN116401883B (zh) | 2023-04-14 | 2023-04-14 | 用于地埋管换热器中断热响应测试的参数反演方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116401883B (zh) |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101430317A (zh) * | 2008-12-09 | 2009-05-13 | 山东建筑大学 | 非恒定功率条件下岩土热物性现场测试方法 |
| CN102721722A (zh) * | 2012-06-20 | 2012-10-10 | 扬州大学 | 一种地下岩土分层热物性现场热响应测试方法 |
| JP5334221B1 (ja) * | 2012-05-11 | 2013-11-06 | 国立大学法人信州大学 | 熱応答試験および揚水試験の解析方法および解析プログラム |
| CN105021647A (zh) * | 2015-06-30 | 2015-11-04 | 长安大学 | 带有断电过程的土壤源热泵岩土热响应试验数据处理方法 |
| CN113656986A (zh) * | 2021-09-16 | 2021-11-16 | 深能科技(山东)有限公司 | 一种快速计算中深层地热地埋管长期运行换热性能的方法 |
| CN114707367A (zh) * | 2022-06-06 | 2022-07-05 | 浙江陆特能源科技股份有限公司 | 中深层地埋管换热器传热分析模型 |
| CN114818537A (zh) * | 2022-04-19 | 2022-07-29 | 中南大学 | 一种地埋管换热器多参数协同反演的正则化方法及系统 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10288548B2 (en) * | 2015-04-17 | 2019-05-14 | Hamilton Sundstrand Corporation | Wavelet-based analysis for fouling diagnosis of an aircraft heat exchanger |
-
2023
- 2023-04-14 CN CN202310396978.3A patent/CN116401883B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101430317A (zh) * | 2008-12-09 | 2009-05-13 | 山东建筑大学 | 非恒定功率条件下岩土热物性现场测试方法 |
| JP5334221B1 (ja) * | 2012-05-11 | 2013-11-06 | 国立大学法人信州大学 | 熱応答試験および揚水試験の解析方法および解析プログラム |
| CN102721722A (zh) * | 2012-06-20 | 2012-10-10 | 扬州大学 | 一种地下岩土分层热物性现场热响应测试方法 |
| CN105021647A (zh) * | 2015-06-30 | 2015-11-04 | 长安大学 | 带有断电过程的土壤源热泵岩土热响应试验数据处理方法 |
| CN113656986A (zh) * | 2021-09-16 | 2021-11-16 | 深能科技(山东)有限公司 | 一种快速计算中深层地热地埋管长期运行换热性能的方法 |
| CN114818537A (zh) * | 2022-04-19 | 2022-07-29 | 中南大学 | 一种地埋管换热器多参数协同反演的正则化方法及系统 |
| CN114707367A (zh) * | 2022-06-06 | 2022-07-05 | 浙江陆特能源科技股份有限公司 | 中深层地埋管换热器传热分析模型 |
Non-Patent Citations (5)
| Title |
|---|
| A method and case study of thermal response test with unstable heat rate;Pingfang Hu 等;Energy and Buildings;第48卷;第199-205页 * |
| Data analysis and discussion of thermal response test under a power outage and variable heating power;Wei Song等;Case Studies in Thermal Engineering;第20卷;第1-11页 * |
| 地埋管换热器热响应试验的实践与优化;郑湍峰;中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑;第2021卷(第02期);第C038-1590页 * |
| 地温未恢复时热响应试验的线热源叠加法;高宽等;暖通空调;第46卷(第04期);第111-114页 * |
| 基于线热源模型计算岩土热物性参数的研究;杨培志等;暖通空调;第47卷(第03期);第113-116页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116401883A (zh) | 2023-07-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110108328B (zh) | 一种供水管网漏损区域漏水量的获取方法 | |
| Bauer et al. | Transient 3D analysis of borehole heat exchanger modeling | |
| Yu et al. | Thermal response test and numerical analysis based on two models for ground-source heat pump system | |
| Bandos et al. | Improving parameter estimates obtained from thermal response tests: Effect of ambient air temperature variations | |
| CN107907564A (zh) | 一种岩土热物性参数与竖直地埋管换热器热阻的确定方法 | |
| Hu et al. | A method and case study of thermal response test with unstable heat rate | |
| Sass et al. | Improvements on the thermal response test evaluation applying the cylinder source theory | |
| Zhang et al. | Parameter estimation of in-situ thermal response test with unstable heat rate | |
| CN111852463A (zh) | 气井产能评价方法及设备 | |
| Zhang et al. | The coupled two-step parameter estimation procedure for borehole thermal resistance in thermal response test | |
| Yu et al. | A simplified model for measuring thermal properties of deep ground soil | |
| CN108088870B (zh) | 非线性热传导模型试验相似准则的建立方法 | |
| JP5334221B1 (ja) | 熱応答試験および揚水試験の解析方法および解析プログラム | |
| CN109726464A (zh) | 一种土石坝流热耦合模型参数敏感性分析构建方法 | |
| CN116401883B (zh) | 用于地埋管换热器中断热响应测试的参数反演方法及系统 | |
| Christodoulides et al. | A practical method for computing the thermal properties of a Ground Heat Exchanger | |
| CN111721831A (zh) | 基于电刺激的三维层析成像堤坝隐伏渗漏通道扫描方法 | |
| CN114818537A (zh) | 一种地埋管换热器多参数协同反演的正则化方法及系统 | |
| CN101556255A (zh) | 一种围护结构热阻现场检测数据分析方法 | |
| CN110991091A (zh) | 地源热泵岩土热物性网络解析与查询系统 | |
| CN109340583A (zh) | 供热管网泄漏监测系统及方法 | |
| Xu et al. | Two-dimensional coupled model of surface water flow and solute transport for basin fertigation | |
| Corcoran et al. | Calibration of thermal response test (TRT) units with a virtual borehole | |
| CN114813828A (zh) | 一种新型确定含水层热物性参数的微热试验方法 | |
| Bandyopadhyay et al. | A New Approach to Modeling Ground Heat Exchangers in the Initial Phase of Heat-Flux Build Up. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |