地源热泵垂直地埋管传热三维数值模型及其建立方法与流程

本发明涉及地源热泵技术领域，具体涉及一种地源热泵垂直地埋管换热器传热数值模型建立方法，适用于在各种不同工况现场热响应测试的岩土体热物性参数计算与垂直地埋管换热器设计。

背景技术：

地源热泵空调系统采用循环液在地埋管中的流动来实现与地下岩土蓄热体的换热。地下岩土体相对空气温度更加稳定，储热与传热能力较强，使地源热泵系统具有高效节能、经济环保的特点。地源垂直埋管方式由于具有相对占用地表面积较少的特点，得到较广泛的应用。在地源热泵系统工程建造中，垂直地埋管换热器的设计是关键部分，而合理设计地埋管换热器的前提是准确测量地下岩土体的换热能力。我国《地源热泵系统工程技术规范《gb50366-2005》(2009年版)新增了岩土热响应试验部分内容，明确了地源热泵系统的岩土热响应试验相关要求，推荐“恒热流法”作为试验实施方法。

目前工程上常用的现场测试的方法有两种：恒定热流法与恒定进水温度法。恒定热流测试法通过对岩土体施加恒定的功率，采集测试过程中各时刻地埋管的进回水温度。利用恒热流的线热源模型或柱热源模型反演计算岩土体的综合热物性参数。求取参数过程简便，但恒定热流法在地埋管换热能力测试方面略有不足，难以获取单位延米换热量；在恒定热流测试过程中，恒定热流的加载功率是随机设定的，不能完全反映实际地源热泵运行状态；另外求取热物性参数时，由于没有使用流量参数与管道口径参数，对流体的质量流量和流速变化产生的换热能力变化不能反映准确。恒定进水温度热响应测试方法为通过控制地埋管的进水温度和水流流量恒定，对岩土体施加冷却或加热的负荷。测试过程中采集地埋管出口温度数据，可以直观地得到一定时间时接近稳定的换热量。由于测试过程加载热流非恒定，获取岩土热物性参数较复杂，并且目前各种求解方法尚不统一，难以按照规范中的计算方法设计地埋管总量。

恒定热流法与恒定进水温度法各有侧重特征，在工程上均有较广泛的应用与研究。专利号cn200810201626、cn200810238160、cn201110223663对恒热流测试方法进行模型修正，提高了求取热物性参数的精度，但仍未定义加载功率与介质流量对测试结果的影响。专利号cn201110224680通过对不同加载负荷、介质流速的进行多次试验，以得到不同结果，进行多因素分析，全面了解测试地埋管换热孔的换热能力。但是测试时间较长，每次测试均受到上次测试时的岩土体温度变化影响。专利号cn201120152222针对地埋管换热器的传热复杂性与传热求解非线性，提出一种面向设计参数应用的孔群总热阻测试方法，简化了岩土热物性参数到设计应用的衔接步骤，但该方法仍使用基本线性热源模型，不能辨识孔内传热的非稳态过程及垂向不同深度传热特征的区别。专利号cn201210026702、cn201210207869、cn201210414520、cn201410479673通过建立测试孔旁侧的监测孔与分层测试处理的方法，利用不同位置处的温度传感器，得到不同深度的温度变化，可以识别不同深度处的岩土体传热特征差异。不足之处为监测孔与测试孔的孔斜角度不易控制，随深度增加二者的距离误差增大，影响数据精度；在整个测试孔的加载功率恒定情况下，不同深度各层的加载功率并非在测试时间上恒定。这些模型中各层的加载功率在测试时间中视为恒定，使误差进一步增大。专利号cn201510181300建立一种单孔和群孔的数值传热计算模型，不但可以计算单个测试孔，还可以通过模拟建筑负荷逐时变化，计算地源换热器孔群的传热响应。该方法使用进回水平均温度概念，对钻孔内的传热视为稳态化，忽略竖向变换及周向温度变化，简化的群孔边界条件也会产生一定误差。

在地埋管换热器的测试与工程应用中，由于其加载负荷的非稳定性与传热非稳态化，逐时数值模拟是最接近实际的传热计算方法。作为理论研究，一些商业软件可以做到地下换热器的数值逐时模拟，对管道间距、回填材料、介质流速等孔内要素也可模拟计算。但作为工程应用则过于复杂、成本较高、计算时间长，缺乏部分应用参数，如地温梯度、分层参数、换热管道形式等。

技术实现要素：

为了综合两种测试工况的数据处理方法，在选取任何一种工况测试后均可求取岩土热物性参数与换热能力，使用同一种计算方法，便于两种工况相互对比；提供地埋管换热的长期逐时模拟，得到准确的温度响应数据。本发明提供一种适当简化但参数齐全的数值模型，该数值模型作为地源热泵地埋管换热器的专用非稳态传热计算工具，要求计算流程明确，多种工况输入，单孔、群孔均可适用。

本发明为解决其技术问题所采用的技术方案为：地源热泵垂直地埋管传热三维数值模型，对模型中的计算单元应用能量守恒原理建立离散方程，每个计算单元具有固定热容，计算单元之间用相应的热阻连接。模型分两个结构层次：水平结构与垂直结构；水平结构部分模拟计算以地埋管为中心的单元径向与圆周方向的传热，垂直结构部分模拟计算地埋管内部流体的轴向温度变化；

水平计算单元为分解的圈层结构，2根或4根换热管保持固定间距，每根管周边的计算单元呈辐射状、层次扩展；

在水平结构上，基本计算公式如下：

每个计算单元单位温度变化需要的热量：

qv＝v0*ρc*k

qv——计算单元的体积热容量；

v0——计算单元的体积；

ρc——计算单元的体积热容；

k——边界系数，在近似模拟群孔状态时，以地埋管换热孔间距为参照，水平模型的外部圈层相对与原有的当量体积减少的比率，取值范围：0＜k≤1；

每个计算单元在不同方向的热阻：

(径向)

(圆周方向与不规则单元的近似处理)

λs＝λ0+a*t

ri——与计算网格单元相邻的单元的热阻；

r0——计算单元的热阻；

m——辐射状分区数目，当圈层单元在距离换热管1-1.5m以外，m＝1；

λs——岩土体地层或管道在一定温度下的热导率；

d0——计算单元弧段所在的圆的外径；

di——计算单元弧段所在的圆的内径；

l——计算单元长度或弧段长度；

s——计算单元截面积；

λ0——岩土体地层或管道在基准温度下的热导率；

a——岩土体地层或管道在不同温度下的热导率变化系数；

t——某时刻的网格单元温度；

r0的计算方法与ri相同；

管道对流换热热阻rf由管道内径、流速、流体温度与流体参数控制，使用经验公式拟合求出，计算需要时叠加到ri之中；

在一个时间步长内，逐个计算所有网格单元的温度响应，而后转入下一个时间步长重新赋值计算：

t＇——下一步长时需计算的网格单元温度；

ti——与计算网格单元相邻的单元在某时刻的温度；

js——每一个计算步长的时间；

n——与计算网格单元相邻的单元的数量；

在垂直结构上，基本计算方法如下：

垂直结构部分，将水平网格结构分成多组叠加构造，用以代表不同深度岩土体地层，换热管内流体在实际中是连续运动的，利用微分原理，将流体切割成一系列长柱状体，将连续流动视为断续运动，在每一个计算步长内，通过每一步的断续推移，模拟管内温度传递状态，管内流体每一步的推移距离与流速、管径相关，与分层模型垂直长度不一致，需要对每一步长的流体管柱重新切割，划分为中间段、进入段、引出段。在每一个计算步长中，将引出段温度计算得出结果后，增加一定温差(恒热流模式)或重新设定温度(恒温度模式)，赋值给予进入段，进入下一步长的计算。

首先求取转换系数kjⁱ：

当(xi-yj)*(yj+1-xi+1)＜0时

当(xi-yj)*(yj+1-xi+1)≥0时

将得到的数值判别是否大于1，大于1的数取值为1，此时以下公式成立：

以tⁱ表示上一时刻各段温度，以tj表示下一时刻各段温度，则

上式中设置地层分层为n层，使用下标变量表示一定流量下水流柱体的长度，上部x1—x2为进入段，x2—x3为第一层，x3—x4为第二层，依次排序，至换热管底部后层号递减，x2n+1—x2n+2为第一层。下部y1—y2为第一层，y2—y3为第二层，依次排序，至换热管底部后层号递减，y2n—y2n+1为第一层，y2n+1—y2n+2为引出段。

依次求出tj后进入水平单元分别计算，计算得到tⁱ的数值结果进入下一循环。

本发明通过多种输入方式实现恒定热流法与恒定进水温度法适用性，并模拟无功测试时流体驱动功耗对岩土体原始温度测试的影响，对孔内回填材料与岩土体的热物性参数准确模拟，对换热管排列等要素进行模拟，并引入地温梯度变化与岩体水平分层参数，选择计算单元体积、每个步长的时间，以体积当量法模拟群孔边界，对群孔传热情况进行长期模拟。

本发明在模型的水平结构上，靠近传热管的传热单元设置为1cm²大小，随距离增大，外围部分传热时的温度梯度逐渐变小，计算单元厚度适当增大，以减少计算单元数量，将局部不规则的单元适当变形处理，试验数据证明对计算精度影响极小，在距离换热管1-1.5m处，将辐射状分层线取消，计算单元变为近似圆形圈层结构，进一步减少计算单元数量，也便于在更靠近外侧的圈层设置孔群参数，计算时间步长设为10-30s。

本发明的有益效果是：

本发明三维数值模型可以对恒热流工况和恒温度工况的现场热响应测试进行数值模拟，计算岩土体热物性参数与地埋管的换热能力，可以使用同一种方法计算两种参数，方便设计应用。

本发明利用数值计算模型，也可以从恒温度模式中导出温度曲线，与实际测试数据拟合，计算岩土体综合参数，这种方法由于引入了流体流量数据，求解结果更加准确。

本发明通过对影响地埋管换热孔传热主要参数的设置，可对实际工程的运行工况进行模拟，如换热孔深度、地温梯度变化、岩土体热导率的温变特征、换热管间距、换热管形式、回填材料传热性能、间歇负荷等，并计算回水或平均进回水温度响应、岩土体空间的温度响应。合理地忽略部分假设条件，使较为复杂的三维数值模型的计算流程相对简洁，在保持较高计算精度同时可占用较少的计算硬件资源，加快运算速度，适用于群孔长期传热计算。

本发明在模拟长期负荷加载时可将换热流体以变流速参数输入，得到不同流速相对应的换热能力，适用于工程设计中调频水泵参数的优化取值研究。

附图说明

图1为单u型数值计算模型水平单元网格局部图。

图2为双u型数值计算模型水平单元网格局部图。

图3为换热管内流体推移计算模型示意图。

图4为岩土体纵剖面温度响应图。

图5为换热管内部流体温度响应图。

图6为岩土体深度分层的单位延米换热负荷变化图。

图7为使用常规线热源方法求取岩土体热导率的误差示意图。

图8为3种流量时的深度-换热能力图。

图9为单个换热孔与完全密集群孔换热能力图。

具体实施方式

本发明建立地埋管换热的三维数值模型，通过多种输入方式实现恒定热流法与恒定进水温度法适用性，对孔内回填材料与岩土体的热物性参数准确模拟，对换热管排列等要素进行模拟，并引入地温梯度变化与岩体体水平分层参数，在保证计算精度情况下，忽略部分次要因素。以岩土体的体积当量法模拟群孔边界，对群孔传热情况进行长期模拟。

本发明三维数值模型计算基本参数：管内流量；进水温度；流动介质在不同温度下的热导率、运动黏度、普朗特数变化；管道热导率；管道外径、内径；管道形式；管道间距；钻孔直径；钻孔深度；恒温层深度；地温梯度；岩土体分层综合热导率、岩土体分层综合热导率在不同温度下的变化率；回填材料热导率、体积热容；回填材料热导率在不同温度下的变化率。

本发明为简化模型，忽略的参数有：流动介质在不同温度下的密度、比热容变化；管道热容视为与回填材料相同；管道与回填材料、回填材料与孔壁的接触热阻；为简化数值计算对体积单元设置的适当变形带来的误差；流动介质在管道内紊流时的温度不均性；孔底u型接头影响；地表气温变化影响；三维模型分层时由地温梯度等因素在垂直方向上温度差异引起的微小热量传递。

本发明数值计算占用计算机资源较大，必须协调计算单元精细程度与计算时间之间的矛盾，选择合理计算单元体积、每个步长的时间。在一般情况下，数值计算模型的网格越精细，越能提高计算精度，反映实际的传热状态。但是过于精细的计算对计算机硬件的要求也较高，使计算时间变长，程序运行调试纠错困难。

建立的模型为数值离散型，分两个层次：水平结构与垂直结构。

附图1、附图2为水平截面上的计算模型水平单元网格局部结构。附图1为单u型2根换热管道，附图2为双u型4根换热管道。在模型的水平结构上，2根或4根换热管保持固定间距，每根管周边的计算单元呈辐射状、层次扩展，靠近传热管的传热单元设置为1cm²大小，随距离增大，外围部分传热时的温度梯度逐渐变小，计算单元厚度适当增大，以减少计算单元数量。将局部不规则的单元适当变形处理，试验数据证明对计算精度影响极小。在距离换热管1-1.5m处，将辐射状分层线取消，计算单元变为近似圆形圈层结构，进一步减少计算单元数量，也便于在更靠近外侧的圈层设置孔群参数。计算时间步长设为10-30s。该模型特点为在控制误差合理范围内计算单元体积增大，计算单元数量减少，适应于长时间(数月-数年)的传热计算。

数值传热模型基本计算公式如下：

每个计算单元单位温度变化需要的热量：

qv＝v0*ρc*k

qv——计算单元的体积热容量；

v0——计算单元的体积；

ρc——计算单元的体积热容；

k——边界系数，在近似模拟群孔状态时，以地埋管换热孔间距为参照，水平模型的外部圈层相对与原有的当量体积减少的比率，取值范围：0＜k≤1；

每个计算单元在不同方向的热阻：

(径向)

(圆周方向与不规则单元的近似变形处理)

λs＝λ0+a*t

ri——与计算网格单元相邻的单元的热阻；

r0——计算单元的热阻；

m——辐射状分区数目，当圈层单元在距离换热管1-1.5m以外，m＝1；

λs——岩土体地层或管道在一定温度下的热导率；

d0——计算单元弧段所在的圆的外径；

di——计算单元弧段所在的圆的内径；

l——计算单元长度或弧段长度；

s——计算单元截面积；

λ0——岩土体地层或管道在基准温度下的热导率；

a——岩土体地层或管道在不同温度下的热导率变化系数；

t——某时刻的网格单元温度。

r0的计算方法与ri相同。

管道对流换热热阻由管道内径、流速、流体温度与流体参数控制，使用经验公式拟合求出，计算需要时与ri叠加。当换热介质流体为纯水时的对流换热热阻rf如下：

λw＝-1.78*10^-5*tw²+2.84*10^-3*tw+0.55

μ＝5.91*10^-4*tw²-5.1*10^-2*tw+1.79

pr＝5.4*10^-3*tw²-4.4*10^-1*tw+13.67

λw——水的热导率；

tw——水的温度；

μ——水的流动动力黏度；

pr——水的普朗特数。

nμ＝0.027*re^0.8*pr^0.33

re——水的雷诺数；

ρ——水的密度；

v——水的流速。

d——管道内径；

nμ——水的努赛尔数；

η——水的对流换热系数；

rf——水的对流换热热阻。

在一个时间步长内，逐个计算所有单元的温度响应，而后转入下一个时间步长重新赋值计算：

t＇——下一步长时需计算的网格单元温度；

ti——与计算网格单元相邻的单元在某时刻的温度；

js——每一个计算步长的时间；

n——与计算网格单元相邻的单元的数量。

在换热管垂向长度上，不同位置管内流体的温度是有较大差别的，简单的计算往往将这种差别平均化，如《地源热泵工程技术规范》附录内容。为了在数值传热模型中表述这种差别与影响，将附图1、附图2所示网格结构分成多组叠加构造，用以代表不同深度岩土体地层。

在模型的垂直方向，每个计算单元按照约5—10m设定，当换热管内水的流速设为0.4—0.6m/s，在时间步长10-30s时，介质流动距离也在一个垂直计算单元长度左右，垂直模型能使用较少的分层数目来反映出不同深度的传热变化。为简化计算，将实际地温梯度引起的垂向计算单元之间的热传递不予考虑，在较大的轴向长度上，对计算结果的影响较小的，其误差可以忽略。

换热管内流体在实际中是连续运动的，利用微分原理，将流体切割成一系列长柱状体，将连续流动视为断续运动。在每一个计算步长内，通过每一步的断续推移，模拟管内温度传递状态。管内流体每一步的推移距离与流速、管径相关，与分层模型的垂直长度不一致，需要对每一步长的流体管柱重新切割。计算方法见附图3。

附图3上部表示某一时刻水流分段温度情况，下部表示下一时刻温度计算结果。x1、x2、y1、y2等表示一维坐标数值变量，设置地层分层为n层，使用变量表示一定流量下水流柱体的长度，上部x1—x2为进入段，x2—x3为第一层，x3—x4为第二层，依次排序，至换热管底部后层号递减，x2n+1—x2n+2为第一层。下部y1—y2为第一层，y2—y3为第二层，依次排序，至换热管底部后层号递减，y2n—y2n+1为第一层，y2n+1—y2n+2为引出段。将引出段温度计算得出结果后，增加一定温差(恒热流模式)或重新设定温度(恒温度模式)，赋值与进入段，进入下一步长的计算。

首先求取转换系数kjⁱ：

当(xi-yj)*(yj+1-xi+1)＜0时

当(xi-yj)*(yj+1-xi+1)≥0时

将得到kjⁱ的数值判别是否大于1，大于1的数取值为1。此时以下公式成立：

以tⁱ表示上一时刻各段温度，以tj表示下一时刻各段温度，则

依次求出tj后进入附图1所示水平单元分别计算，计算得到tⁱ的数值结果进入下一循环。

实施例

附图4为不同加载时间后的单个地埋管换热孔的岩土体纵剖面温度响应，左图表示在加载时间较短时，换热孔中心的外侧的显示地温梯度垂向层状特征明显；右图为加载时间较长时，地层的原始温度梯度与加载负荷共同影响了岩土体温度变化。

附图5为换热管内部流体温度分布状态，左图为恒热流工况,自下至上4条曲线的时间分别为24h、48h、72h、96h，横坐标深度0m处为换热管水流进出口温度，该处温差在4个时间点上都是相同的；右图为恒温度工况，4条曲线的时间与左图相同，在换热管进口温度不变的情况下，曲线下部逐渐上移，换热管出口温度逐渐变大，随着恒温度负荷加载时间延长，出口温度变化幅度变小。

附图6为在恒热流模式下岩土体不同深度的平均延米换热功率变化情况。在三维数值模型中，设置换热孔深100m，在40-50m处存在热导率较大的异常地层，其他深度均为相对较小的数值。可见，对整个地埋管换热器施加完全恒定的负荷时，地埋管不同深度位置的单位换热功率随时间变化有较大变化，并非与进回水口所测量得到的恒定数据一致，其中约45m深度处的功率变化更加明显。这种换热功率随时间的变化只有在测试时间较长时才会减小(40-50h以后)，所以直接利用恒热流线热源模型来分层计算岩土热导率的方式会产生一定的误差，必须利用参数全面的三维模型予以纠正。

附图7为使用三维模型校正常规线热源模型的计算误差，常规线热源模型在恒热流工况下求取的岩土体热导率参数时，受介质流体流量变化与换热孔深度影响较大。当流量越小时，由于流体质量流量变化与壁面换热系数变化，所使用常规方式求取的岩土体综合热导率偏差越大，并且随换热井深度加大，误差增大，随测试的加载方式不同，夏季工况与冬季工况求取的的岩土体热导率数值偏差反向。

附图8为在恒温度模式下，使用三维模型计算换热孔在不同流体流量参数下受深度变化影响的程度，当换热孔深度加大时，需要提高流体流量，才可得到相应增大的换热能力。受地温梯度影响，夏季工况与冬季工况在深度加大时的换热能力变化有所差异，夏季工况时随着深度加大，其换热能力增加的幅度相对较小。

附图9表示了对单个地埋管换热孔与当量体积法完全密集群孔的换热能力的长时间模拟。可见，完全密集群孔状态下，各种不同参数换热孔的换热能力在10-20d后均逐渐减小，其减小的幅度超过单个换热孔；在较长时间之后，换热能力曲线呈近似线性下降状态，与单个换热孔近似对数曲线下降的形态逐渐分离。在实际工程中，一定数量的群孔的换热能力曲线介于二者之间，其温度响应与换热功率可以使用三维数值模型当量体积法模拟计算。