一种地埋管换热器传热耦合模拟降阶方法

摘要

一种地埋管换热器传热耦合模拟降阶方法，包括以下步骤：步骤1，建立控制方程；步骤2，边界确定；步骤3，得到地埋管换热器传热问题数学模型的状态空间形式的原始系统；步骤4，对步骤3得到的原始系统，采用模型降阶方法进行降阶，得到降价系统式；步骤5，对降阶系统式进行数值求解，得到降阶系统的解；步骤6，通过步骤4求得的降阶系统系数矩阵和步骤5求得的降阶系统状态变量，计算得到输出变量，即为原始系统输出变量的近似解，即为控制方程式在相应的边界条件下输出变量的数值解。本发明降阶系统的规模远小于原始系统，只需求解一个规模较小的系统，从而节约了计算资源，大大缩短计算程序的运行时间。

权利要求

1.一种地埋管换热器传热耦合模拟降阶方法，其特征在于，包括以下步骤：步骤1，分析地埋管换热器与岩土之间传热的特征，建立控制方程：地埋管换热器径向和埋深方向的二维计算模型，具体过程为：利用等效直径法将U型管地埋换热器简化为一直管段，并利用径向的对称性，将控制方程表达为径向r向和埋深方向z向的二维计算模型；

基于换热器内的流体和岩土固体的不同传热机理，控制方程包括地埋管换热器内流体对流换热的控制方程和地埋管换热器外岩土导热的控制方程两部分：地埋管管内流体对流换热控制方程如下：

岩土导热控制方程如下：

3

ρ：密度，kg/m ；cp：比热容，J/(kg·℃)；w：轴向流速，m/s；T：温度，℃；λ：导热系数，W/(m·℃)；λt：湍流导热系数，W/(m·℃)；t：时间，s；

在层流状态下，湍流导热系数λt为0，轴向流速w由基于流体截面平均流速wm的解析解获得；

在湍流状态下，通过求解动量控制方程并应用混合长度理论求解湍流粘性系数得到速度分布，控制方程如下：η：动力粘性系数，kg/(m·s)；

ηt：湍流粘性系数，kg/(m·s)；

p：压力，Pa；

先假定一个速度场，根据经验公式计算混合长度lm，进而计算出湍流粘性系数ηt，然后根据湍流粘性系数ηt和动力粘性系数η求解动量方程，获得新的速度场，反复上述过程，直到收敛，获得速度场及ηt；再求得湍流导热系数λt，就进行温度控制方程的求解；

步骤2，地埋管管内流体进口温度、岩土远边界处温度、上下边界的温度和岩土的初始温度的确定，具体过程为：地埋管管内流体进口温度根据所模拟的工况确定；岩土远边界处温度与岩土的原始温度相同；上下边界的温度为此深度下岩土的原始温度，根据地质勘测得到的地温数据确定；岩土的初始温度应根据模拟工况开始时岩土的温度确定；

步骤3，采用有限容积法，对步骤1建立的控制方程式(1)和(2)进行空间离散，结合步骤

2确定的边界条件，得到地埋管换热器传热问题的状态空间形式的原始系统式，如下：式中， 为系统的状态变量，即为各离散节点的温度,其中n为空间离散节点数；

为系统输入变量，在此模型中为岩土的边界温度及地埋管的进口流体温度，其中p为输入变量数； 为系统输出变量，为地埋管的出口流体温度，m是输出变量数；

为系数矩阵，其中A和B主要与流体和岩土的物性参数、节点的位置有关，C和D则由所需要的输出变量决定；

步骤4，对步骤3得到的原始系统式，采用模型降阶方法进行降阶，得到r阶降阶系统如下，使得r＜＜n：式中， 为降阶系统相应的系数矩阵，

V为通过模型降阶方法获得的投影矩阵； 为降阶系统的状态变量， 为降阶系统的输出变量；

式(5)中各系数矩阵 的求解过程如下：

1)输入系数矩阵A、B、C、D及降阶的阶数r；

‑1

2)构造计算空间：将原始系统通过Laplace变换得到传递函数H(s)＝C(sI‑A) B，如果存在一点 使得A‑s0I非奇异，则可构造出原始系统的Krylov子空间其中

3)对得到的两个计算空间应用Arnoldi过程获得计算空间的一组标准正交基底：(1)初始化向量：

(2)计算第二个标准正交向量：

(3)计算第r个标准正交向量：

(4)构建正交投影矩阵V＝[v1,v2,…,vr]；

(5)分别通过关系式 求解降阶系统的系数矩阵，输出降阶系统的系数矩阵：

步骤5，对降阶系统式进行数值求解，得到降阶系统的解；

步骤6，通过步骤4求得的降阶系统系数矩阵 和步骤5求得的降阶系统状态变量由式(4b)计算得到输出变量 即为原始系统式(4)输出变量y(t)的近似解，即为控制方程式(1)、(2)在相应的边界条件下输出变量的数值解。

说明书

一种地埋管换热器传热耦合模拟降阶方法

技术领域

[0001] 本发明属于传热数值计算方法领域，特别涉及一种地埋管换热器传热耦合模拟降阶方法。

背景技术

[0002] 地热能具有储量大、分布广、清洁环保、稳定可靠等特点。实现地热能供暖(制冷)主要是通过地源热泵系统，而提高该系统效率的关键是提高地埋管换热器的换热性能。但地埋管换热器在岩土中换热是一个复杂多变的动态传热过程，涉及的因素很多，如岩土热物性参数、换热器的结构设置、流体运行参数等。利用浅层地热的地埋管换热器埋深为100‑300m，而中深层地埋管换热器的埋深可达2000m，而且不仅仅是单个地埋管换热器的计算，还需要对整个井群进行计算，所计算的空间区域大。此外，由于换热器中流体与土壤传热相互影响，且随着时间而变化，需对地埋管换热器进行全生命周期传热过程的模拟，故计算的时间域长。因此对该类换热问题，数值模拟仿真受到了很大的挑战，急需寻找准确、高效的数值计算方法来指导地埋管换热器的设计，以便能够改善地埋管换热器的效率。

[0003] 模型降阶(MOR)方法是将一个较大复杂系统转化为一个近似的较小系统,从而降低大型复杂系统的理论分析难度，减少数据存储量和运算量，加速系统的模拟计算，因此在大规模集成电路、自动化控制和机械、建筑性能模拟等工程领域得到许多应用。然而，在地埋管换热器传热性能的模拟领域还没有模型降阶方法的应用。

发明内容

[0004] 本发明的目的在于一种地埋管换热器传热耦合模拟降阶方法，以解决上述问题。

[0005] 为实现上述目的，本发明采用以下技术方案：

[0006] 一种地埋管换热器传热耦合模拟降阶方法，包括以下步骤：

[0007] 步骤1，分析地埋管换热器与岩土之间传热的特征，建立控制方程：地埋管换热器径向和埋深方向的二维计算模型；

[0008] 步骤2，地埋管管内流体进口温度、岩土远边界处温度、上下边界的温度和岩土的初始温度的确定；

[0009] 步骤3，采用有限容积法，对步骤1建立的控制方程进行空间离散，结合步骤2确定的边界条件，得到地埋管换热器传热问题数学模型的状态空间形式的原始系统；

[0010] 步骤4，对步骤3得到的原始系统，采用模型降阶方法进行降阶，得到降价系统式；

[0011] 步骤5，对降阶系统式进行数值求解，得到降阶系统的解；

[0012] 步骤6，通过步骤4求得的降阶系统系数矩阵和步骤5求得的降阶系统状态变量，计算得到输出变量，即为原始系统输出变量的近似解，即为控制方程式在相应的边界条件下输出变量的数值解。

[0013] 进一步的，步骤1中具体为：利用等效直径法将U型管地埋换热器简化为一直管段，并利用径向的对称性，将控制方程表达为径向r向和埋深方向z向的二维计算模型；

[0014] 基于换热器内的流体和岩土固体的不同传热机理，控制方程包括地埋管换热器内流体对流换热的控制方程和地埋管换热器外岩土导热的控制方程两部分：

[0015] 地埋管管内流体对流换热控制方程如下：

[0016]

[0017] 岩土导热控制方程如下：

[0018]

[0019] ρ：密度，kg/m3；cp：比热容，J/(kg·℃)；w：轴向流速，m/s；T：温度，℃；λ：导热系数，W/(m·℃)；λt：湍流导热系数，W/(m·℃)；t：时间，s。

[0020] 进一步的，在层流状态下，湍流导热系数λt为0，轴向流速w由基于流体截面平均流速wm的解析解获得；

[0021] 在湍流状态下，需要通过求解动量控制方程并应用混合长度理论求解湍流粘性系数得到速度分布；先假定一个速度场，根据经验公式计算混合长度lm，进而计算出湍流粘性系数ηt，然后根据湍流粘性系数ηt和动力粘性系数η求解动量方程，获得新的速度场，反复上述过程，直到收敛，获得速度场及ηt；再求得湍流导热系数λt，就可以进行温度控制方程的求解。

[0022] 进一步的，步骤2中具体为：地埋管管内流体进口温度根据所模拟的工况确定；岩土远边界处温度与岩土的原始温度相同；上下边界的温度为此深度下岩土的原始温度，根据地质勘测得到的地温数据确定；岩土的初始温度应根据模拟工况开始时岩土的温度确定。

[0023] 进一步的，步骤3中原始系统式如下：

[0024]

[0025] 式中， 为系统的状态变量，即为各空间离散节点的温度,其中n为空间离散节点数； 为系统输入变量，在此模型中为岩土的边界温度及地埋管的进口流体温度，其中p为输入变量数； 为系统输出变量，为地埋管的出口流体温度，m是输出变量数； 为系数矩阵，其中A和B主要与流体和岩土的
物性参数、节点的位置有关，C和D则由所需要的输出变量决定。

[0026] 进一步的，步骤4中具体为：

[0027] 对步骤3得到的n阶原始系统式(3)，采用模型降阶方法进行降阶，得到r阶降阶系统如下，使得r＜＜n；

[0028]

[0029] 式中， 为降阶系统相应的系数矩阵，V为通过模型降阶方法获得的投影矩阵； 为降
阶系统的状态变量， 为降阶系统的输出变量。

[0030] 进一步的，步骤6中具体为：通过步骤4求得的降阶系统系数矩阵 和步骤5求得的降阶系统状态变量 由式(3b)计算得到输出变量 即为原始系统式(3)输出变量y(t)的近似解，即为控制方程式(1)、(2)在相应的边界条件下输出变量的数值解。

[0031] 与现有技术相比，本发明有以下技术效果：

[0032] 本发明中的模型降阶方法是对状态变量在频率域进行级数展开，构建出一个近似的系统，计算此系统的系数矩阵，从而获得降阶系统。降阶系统的规模远小于原始系统，我们只需求解一个规模较小的系统(r阶，r＜＜n)，从而节约了计算资源，大大缩短计算程序的运行时间。

[0033] 本发明在保持所需精度的条件下，减少模拟计算的耗时，为地埋管换热器的设计及全生命周期的模拟提供一种快速的计算方法，从而节约成本和时间，而且所计算的空间与时间域越大该计算方法的高效性月明显。Krylov子空间方法可以通过在不同的点展开从而获得不同计算精度的降阶系统，对不同的模型的计算具有很好的适应性和灵活性，便于在工程设计中应用。

附图说明

[0034] 图1是地埋管换热器与岩土间传热示意图

[0035] 图2是地埋管换热器与岩土间传热的计算模型图

[0036] 图3是地埋管换热器传热模型降阶方法的技术路线图；

[0037] 图4是计算域节点划分示意图；

[0038] 图5是直接求解及实验的进出口水温图；

[0039] 图6是直接求解与KS‑MOR求解的出口温度的相对误差随时间的变化图

[0040] 图7是直接求解与KS‑MOR求解的出口温度的平均相对误差图；

[0041] 图8是在不同阶数下KS‑MOR方法的计算耗时图；

[0042] 图9是不同计算的时间域下KS‑MOR的计算耗时图；

[0043] 图10是不同计算的时间域下的KS‑MOR/直接求解的计算总耗时的比例图；

具体实施方式

[0044] 下面结合附图和具体实施方式对本发明进行详细说明。

[0045] 地埋管换热器传热的降阶求解方法，以Krylov子空间模型降阶方法为例，在频率域上对系统的状态变量展开，然后构造所需降阶矩阵，其技术路线如图3所示，详细步骤如下。

[0046] 步骤一：分析地埋管换热器与岩土之间的传热的过程，建立地埋管换热器传热的物理数学模型。

[0047] 基于换热器内的流体和岩土固体的不同传热机理，控制方程包括地埋管换热器内流体对流换热的控制方程和地埋管换热器外岩土导热的控制方程两部分。

[0048] 地埋管管内流体对流换热控制方程如下：

[0049]

[0050] 岩土导热控制方程如下：

[0051]

[0052] ρ：密度，kg/m3；

[0053] cp：比热容，J/(kg·℃)；

[0054] w：轴向流速，m/s；

[0055] T：温度，℃；

[0056] λ：导热系数，W/(m·℃)；

[0057] λt：湍流导热系数，W/(m·℃)；

[0058] t：时间，s。

[0059] 在层流状态下，湍流导热系数λt为0，轴向流速w由下式确定：

[0060] w/wm＝2[1‑(r/R1)2]  (3)

[0061] R1：管道半径，m。

[0062] wm：截面平均流速，m/s。

[0063] 在湍流状态下，需要通过求解动量控制方程并应用混合长度理论求解湍流粘性系数得到速度分布。

[0064] 动量控制方程：

[0065]

[0066] η：动力粘性系数，kg/(m·s)；

[0067] ηt：湍流粘性系数，kg/(m·s)；

[0068] p：压力，Pa。

[0069] 将动量方程(3)行无量纲化处理：

[0070]

[0071] 其中，ζ＝r/R1，

[0072] 应用混合长度理论求解湍流粘性系数：

[0073]

[0074] 混合长度按下式计算：

[0075] lm＝(DF)V·(lm)N  (7)

[0076] 其中(lm)N按Nikuards公式计算，而阻尼因子(DF)V按下式确定：

[0077] (lm)N＝[0.14‑0.08(r/R1)2‑0.06(r/R1)4]R1  (8)

[0078]

[0079] 其中，τw为壁面切应力，A＝26。

[0080] 湍流导热系数 σT为湍流Prandtl数，通常取0.9～1。

[0081] 因为ηt与速度有关，因此动量方程需要通过迭代进行求解。先假定一个速度场，根据式(7)、(8)及(9)，确定ηt，然后求解动量方程，反复此过程，直到收敛，获得速度场及ηt。再求得λt，就可以进行温度控制方程的求解。

[0082] 步骤二：初始条件和边界条件的确定与建立。

[0083] 为了对控制方程(1)和(2)进行数值求解，需要根据问题确定其初始条件和边界条件。地埋管管内流体进口温度根据所模拟的工况确定。岩土远边界处温度不再受地埋管换热器的影响，因此与岩土的原始温度相同，上下边界的温度也看作是此深度下岩土的原始温度，可根据地质勘测得到的地温数据确定。岩土的初始温度应根据模拟工况开始时岩土的温度确定，根据不同工况确定。在本例中，其具体的边界条件如下：

[0084]

[0085] Tin：进口温度，℃；

[0086] TW：壁面温度，℃。

[0087] 步骤三：地埋管换热器传热控制方程的空间离散，得到求解问题的原始系统。

[0088] 如上所示，地埋管换热器传热过程的控制方程(1)、(2)为关于时间和空间的偏微分方程组，采用数值求解时，首先需要进行空间离散，转化为关于时间的常微分方程组。计算区域的节点划分如图3所示，径向k1个节点，轴向k2个节点，共n＝k1·k2个内节点，即n个控制容积。其中(Δr)i表示径向第i个控制容积的径向长度，(δr)i表示第i个节点到第(i+1)个节点的距离；(Δz)j表示轴向第j个控制容积的轴向长度，(δz)j表示第j个节点到第(j+1)个节点的距离。这里采用有限容积法，对地埋管换热器传热偏微分方程(1)、(2)进行空间离散，结合边界条件(10)，得到状态空间形式描述的地埋管换热器传热数学模型的原始系统，如下式(11)。

[0089]

[0090] 式中， 为系统的状态变量，即为各空间离散节点的温度，其中n为空间离散节点数； 为系统输入变量，在此模型中为岩土边界的温度及地埋管换热器的进口流体温度，其中p为输入变量数，p＝4； 为系统输出变量，可以是地埋管换热器的出口流体的平均温度，m是输出变量数，m＝1。

[0091] 为系数矩阵，其中A和B主要与流体和岩土的物性参数、节点的位置有关，C和D则由所需的输出变量确定，

[0092] 为系数矩阵，当以岩土远边界及上下边界的温度及地埋管换热器的进口流体温度为输入变量，以地埋管换热器的出口流体的平均温度为输出变量时，变量T(t)、u(t)、y(t)的表达式和系数矩阵A、B、C、D计算式如下：

[0093]

[0094]

[0095]

[0096]

[0097]

[0098]

[0099] ai，j，P＝0i，j，E+ai，j，W+ai，j，N+ai，js，i＝1，2，…，k1；j＝1，2，…，k2[0100] 步骤四：对原始系统采用Krylov子空间降阶方法进行降阶，得到降阶系统。如下式(12)，使得降阶系统的阶数r远小于原始系统的阶数n。一般r减小，计算误差增大而计算量减少，因此对于每一具体问题，先需实验找到合理的降阶系统的阶数。

[0101]

[0102] 式(12)中各系数矩阵 的求解过程如下：

[0103] 1)输入系数矩阵A、B、C、D及降阶的阶数r；

[0104] 2)构造计算空间：将原始系统通过Laplace变换得到传递函数H(s)＝C(sI‑A)‑1B，如果存在一点 使得A‑s0I非奇异，则可构造出原始系统的Krylov子空间其中

[0105] 3)对得到的两个计算空间应用Arnoldi过程获得计算空间的一组标准正交基底:

[0106] (1)初始化向量：

[0107] (2)计算第二个标准正交向量：

[0108] (3)计算第r个标准正交向量：

[0109] (4)构建正交投影矩阵V＝[v1,v2,…,vr]；

[0110] (5)分别通过关系式 求解降阶系统的系数矩

[0111] 阵，输出降阶系统的系数矩阵：

[0112] 步骤五：对降阶系统式(12a)进行数值求解，得到降阶系统的解

[0113] 步骤六：根据式(12b)，降阶系统中的输出变量 可由步骤四获得的降阶系统系数矩阵 和步骤五获得的降阶系统状态变量 计算得到，即为原始系统式(11)输出变量y(t)的近似解；原始系统式(11)的状态变量T(t)，可根据关系式 及步骤
四求得的投影矩阵V和降阶系统状态变量 计算求得为控制方程式(1)、(2)在边界条件(10)下的数值解。

[0114] 下面结合具体案例，依据上述技术路线的详细步骤，编写并运行计算程序，通过计算结果分析本发明的地埋管换热器传热模型降阶方法的可靠性及有益效果。

[0115] 不失一般性，以某U型管地埋管换热器放热过程为例。采用等效直径的方法将U型管简化为一直管段后各部分的半径为R1＝0.0378m，R2＝0.0406m，R3＝0.08m，R4＝2.5m。地埋管换热器的埋深为L＝100m。岩土的初始温度为18℃，远边界温度为17℃，计算的时间域为48h，流体进口水温随时间变化如图5。流体、管壁、回填土壤及岩土的物性参数如表1所示。

[0116] 表1流体、管壁、回填土壤及岩土的物性参数

[0117]

[0118] 然后对控制方程进行离散，径向网格的总内节点数为k1＝40，其中，流体的内节点数为kp＝8，轴向网格的总内节点数为k2＝100，总的内节点数n＝4000，即原始系统的阶数为4000阶。时间步长取5s。在利用Krylov子空间模型降阶方法(KS‑MOR)进行降阶处理时，降阶系统的阶数分别取为r＝10,20,40,50,80,100,200。

[0119] 针对该案例的上述已知条件(流体、管道、回填土壤及岩土的参数)，依据图3技术路线的详细步骤，编写并运行计算程序，得到采用KS‑MOR和直接求解方法数值模拟结果。所谓直接求解方法，是直接对步骤三得到的原始系统进行数值求解，即没有上述步骤四(投影矩阵及降阶系统系数矩阵的构造)和步骤五(降阶系统的求解)。计算程序采用Matlab R2016a编写，用于运行程序的工作站主要配置为：CPU E5‑2630v4@2.2GHz，内存64G，显卡NVIDIA Quadro P2000。

[0120] 对于地埋管换热器传热的问题，基于Krylov子空间与直接求解方法的计算结果及其与实验结果的对比如下。

[0121] 图5是直接求解及实验的进出口水温，图6是直接求解与KS‑MOR求解的出口平均温度的相对误差随时间的变化，图7直接求解与KS‑MOR求解的出口平均温度的平均相对误差。

[0122] 相对误差定义及平均相对误差的定义如下：.

[0123]

[0124]

[0125] 其中， 分别是直接求解和KS‑MOR求解的出口平均温度，m为时间节点数。

[0126] 可以看出KS‑MOR方法与直接求解方法的计算结果吻合程度高，两者之间误差很小，在所计算的降阶系统的阶数中，最大的相对误差小于7％，最大的平均相对误差小于0.7％，而且当降阶系统的阶数大于40阶，平均相对误差逐渐趋于稳定且小于0.1％。由此可见，Krylov子空间模型降阶方法的计算准确度高。

[0127] 该KS‑MOR方法的程序运行耗时(即计算耗时)如图7所示，而直接求解方法的计算总耗时、构造矩阵耗时及求解耗时分别为19984.45s、0.14s和19984.31s，可见KS‑MOR方法的计算总耗时及求解耗时远小于直接求解方法的耗时，KS‑MOR方法构造矩阵耗时略大于直接求解方法。

[0128] 综合考虑相对误差和计算耗时，可以看出对于该模型降阶系统的阶数选择40‑50阶比较合适，即为原始系统的1％左右。

[0129] 然后采用KS‑MOR计算获得阶数为40和50阶的降阶系统，并进行了计算时间域分别为6h、12h、24h、36h、48h的数值模拟。图8、9给出了不同的计算时间域下Krylov子空间方法的计算耗时及其与直接求解计算耗时的对比，可以看出KS‑MOR方法的计算总耗时随着计算时间域的增加而增加。而且当计算的时间域是从6h增加到48h时，KS‑MOR方法的计算总耗时从直接求解总耗时的0.182％(r＝40)和0.178％(r＝50)减少到了直接求解耗时的0.043％(r＝40)和0.051％(r＝50)。对于直接求解和在降阶系统的阶数一定时KS‑MOR求解，随着计算时间域的增加，构造矩阵耗时不变，求解耗时与计算的时间域成正比；而且在一个时间节点上，直接求解需要求解n阶的方程而KS‑MOR只需要求解r阶的方程，前者耗时远大于后者，因此随着时间域的增加，直接求解耗时快速增加，而KS‑MOR则相对较慢。

[0130] 由此可见，模型降阶方法不仅计算准确度高，而且能有效缩短程序运行时间，计算效率显著提高。因此，当计算的时间域或者空间增加时，KS‑MOR的计算耗时远远小于直接求解，其高效性也会更加明显。