张 琦，苑中显，郭占全

（北京工业大学 环境与能源工程学院，北京 100124）

0 引言

太阳能是一种安全、清洁且分布广泛的可再生能源。为了充分地利用太阳能，在没有热量需求时，对太阳能进行收集并加以储存，从而提高太阳能的全年利用率。目前，太阳能长周期跨季节储热技术为研究热点[1]～[3]。

带有储热子系统的太阳能供热系统如图1所示。

图1 带有储热子系统的太阳能供热系统Fig.1 Trans-seasonal solar heating system with a heat storage tank

该系统中的储热水箱埋于地下。在非供暖季，该系统将太阳能集热器所获得的热量以显热的形式存储于储热水箱内；在供暖季，将储热水箱中的热水送至建筑末端，从而实现跨季节供暖。储热水箱的侧壁斜度会影响自身的热量储存效率，这是由于贴近壁面处的水受壁面冷却，密度随之增大，当水在水箱侧壁自上而下流动时，若水箱侧壁为竖直面，则水流不会受到壁面的阻碍，流速较大，产生的热量损失也较大；若水箱侧壁为倾斜面，则水流会受到倾斜壁面的阻碍，流速相对减小，产生的热量损失也较小，因此倾斜的壁面有利于提高储热水箱的热量储存效率。

目前，鲜有研究储热水箱壁面倾斜度对自身自然对流热损失的影响。本文主要研究埋于地下的漏斗状储热水箱（以下简称为储热水箱）侧壁斜度对其自然对流换热所产生的热损失的影响，并将储热水箱简化为由土壤包裹的封闭腔体，以便进行模拟研究。

1 物理模型

储热水箱周围的土壤温度远低于其内部的热水温度。储热水箱中的热水会与土壤发生热交换，导致储热水箱内的水温随之降低，自身温度场也变得不均匀，使得水箱内的水产生流动，形成自然对流。当储热水箱侧壁存在一定的倾斜度时，该壁面会阻碍流体的流动，壁面倾斜度不同，此阻碍作用也不同。改变侧壁与地面之间的夹角（水箱侧壁倾角）θ，可得到不同的自然对流换热强度。

储热水箱的模型图见图2。

图2 储热水箱模型图Fig.2 Numerical domain of the water tank of funnel shape

由图2可知，本文的储热水箱呈四棱台结构。令该水箱上、下底面的边长分别为L1，L2，高为H；土壤的长、宽、高分别为a，b，h。本文将储热水箱的下底面几何尺寸和高设为定值，通过改变θ，构造出不同侧壁倾角的水箱。

模拟过程中，储热水箱的相关参数见表1。

表1 储热水箱的相关参数Table 1 Related parameters of hot water storage tank

计算区域中，水和土壤的物性参数随温度的变化很小，故视为定值。其中，水的密度为979.8 kg/m3，定压比热为 4 214 J/（kg·K），导热系数为0.682 W/（m·K），动力粘度η为0.000 3 Pa·s，容积膨胀系数为0.000 5；土壤的密度为2 000 kg/m3，定压比热为 1 010 J/（kg·K），导热系数为 1.16 W/（m·K）[4]。

2 数学模型

式中：ρ0为土壤密度；为土壤被加热后的温度；TS为土壤初始温度；β为土壤的体积膨胀系数。

2.1 控制方程

在上述假设的基础上，整个计算区域内的流体流动以及传热过程满足以下微分方程。

连续性方程为

式中：u为流体流速；t为时间；ρ为流体密度。

动量方程为

式中：p为流体压力；T为流体温度；T0为参考温度；μ为流体运动粘度；β为流体容积膨胀系数；μt为流体湍动粘度。

能量方程为

式中：σT为湍流普朗特数，取0.9。

湍动能方程为

式中：k为湍流动能；σk为湍流动能k相应的普朗特数，取1；ε为湍流耗散率；Gk为流体平均速度梯度所引起的湍流动能的产生项。

式中：C1ε为小涡拉伸项，取 1.44；C2ε为黏性破坏项，取1.92；σε为ε相应的普朗特数，取1.3。

计算区域内流体湍动粘度μt的计算式为

对于计算区域内流体平均速度梯度所引起的湍动能k，其产生项Gk的表达式为

2.2 初始条件

储热水箱中水的初始温度T0=368K，外围土壤的初始温度TS=295 K，水与土壤之间耦合换热。

2.3 边界条件

远离水箱处的土壤温度不会受到水箱影响，故设置土壤模型的边界为绝热边界。

3 数值求解

3.1 网格划分

本文用CFD前处理软件Gambit建立储热水箱和外围土壤的三维模型并生成三维实体计算网格。整个计算区域划分为水箱和土壤两部分：采用非均分网格在水体和土壤交界处进行局部细化，远离水箱的土壤受水箱影响很小，因此，该处的网格较为稀疏。由于控制了水箱的底部边长和高，各倾角水箱的体积不同，生成的网格数也不同。

当储热水箱侧壁倾角θ为60°时，储热水箱和外围土壤三维模型的网格划分如图3所示。

图3 当θ为60°时，储热水箱和外围土壤三维模型的网格划分情况Fig.3 When θ=60 °，meshing of 3D model of hot water storage tank and surrounding soil

3.1.1 网格无关性验证

当 θ为 90°，网格数分别为 28 061，36 286，50 419，75 993时，通过计算得到储热水箱的平均温度，而后根据计算结果进行网格无关性验证。网格无关性验证结果见图4。

图4 网格无关性验证结果图Fig.4 The average tank temperature changing with the grid number

由图4可知，当网格数由36 286逐渐增加至50 419时，储热水箱平均温度的变化量小于0.01数量级，因此，可认为当网格数为50 419时，已达到网格无关。

本文还分别对 θ为 15，30，45，60，75，90 °条件下的模型进行网格无关性验证。验证结果表明，当 θ分别为 15，30，45，60，75，90 °时，本文模型的网格数分别应取 1 113 027，3 649 938，416 826，207 819，98 221，50 419。

3.1.2 时间步长无关性验证

通过计算发现，当时间步长分别为1，0.5，0.1，0.05s时，储热水箱的平均温度分别为367.407，367.526，367.655，367.657 K。由此可见，当时间步长为由0.1 s逐渐减小至0.05 s时，储热水箱的平均温度变化量在0.01数量级，因此，本文取时间步长为0.1 s。

3.2 求解器设置

本文采用的数值方法为基于原始变量的SIMPLEC方法，压力项采用PRESTO!方案进行离散，动量方程和能量方程采用二阶迎风格式进行离散。

土壤与水体采取耦合换热，这样能够避免将壁面设置为恒温或恒热流边界条件时，土壤对流换热模拟结果的误差。

4 数值模拟结果与分析

图5为当θ为60°时，储热水箱及其周围土壤的温度分布情况。由图5可知，土壤的温度低于储热水箱，导致储热水箱的壁面温度较低，储热水箱壁面处的水被冷却，密度变大并沿壁面向下流动，远离壁面处水的温度较高，密度也较小，并填补壁面处向下流动的水，所以水箱中存在自然对流换热。被冷却的水沉积于底部，使储热水箱出现温度分层现象，水箱底部水的温度最低。当储热水箱底部的水逐渐远离该处，温度逐渐升高。综上可知，本文的储热水箱可以实现分层取水，具体的取水方式：需要低温水时，从储热水箱的下部区域取水；需要高温水时，从储热水箱的上部区域取水。

图5 当θ为60°时，储热水箱及其周围土壤的温度分布情况Fig.5 Temperature distribution of the water tank and the surrounding soil as θ＝60 °

4.1 倾角对流场的影响

流体沿储热水箱壁面的自然对流换热强度会受到自身运动高度的影响。储热水箱侧壁倾角会影响流体沿壁面的运动高度。建立储热水箱剖面的坐标系如图6所示。

图6 储热水箱剖面的坐标系Fig.6 The coordinate system of the water tank for the numerical simulation

图7为当储热水箱侧壁倾角θ为15～90°时，该水箱内部流体的流动情况。由图7可知，当θ较小（θ=15，30 °）时，接近壁面处出现多个涡旋（当θ=15°时，水箱两侧各存在3个涡旋；当θ=30°时，水箱两侧各存在2个涡旋）。这是由于当θ较小时，水箱壁面斜度较平缓，当流体沿壁面向下流动时，受到倾斜壁面的阻碍作用较大，因此生成涡旋较多。涡旋越多，涡旋附近水的流速越小，导致水箱壁面的对流换热强度减弱；随着θ逐渐增大（θ=45～90°），水箱两侧的涡旋各减少至1个，且涡旋逐渐向壁面靠近。由图6，7可知，当θ=45，60，75，90 °时，储热水箱右侧涡旋中心的坐标分别为 （9.48，0.34），（8，-0.26），（6.85，-0.83），（5.9，-0.71），该涡旋与储热水箱右侧壁面之间的垂直距离分别为2，1.5，1.1，1.1 m。这是由于随着θ逐渐增大，水箱壁面对流动的阻碍作用逐渐减弱，导致涡旋与储热水箱壁面之间的距离大体上逐渐减小。由图7还可看出，涡旋距离壁面越近，壁面处水的速度梯度和流动速度越大，从而导致水箱壁面的对流换热强度越大。综上可知，随着θ逐渐增大，水箱壁面的对流换热强度逐渐增大。

图7 当θ为15～90°时，该水箱内部流体的流动情况Fig.7 The streamline pattern of the natural convection flow in the tank as θ varies from 15 °to 90 °

4.2 温度、速度分布情况

本文通过数值模拟计算得到，不同θ条件下，储热水箱中心截面（x=0）处温度场的分布情况如图8所示。

图8 不同θ条件下，储热水箱中心截面处（x=0）的温度场分布情况Fig.8 Distribution of the temperature in the tank at x=0 palne as θ varies

由图8可知，储热水箱顶部、底部区域温度较低，其他区域水的温度大体一致。由于储热水箱顶部和底部区域受土壤冷却作用明显，导致这两个区域中水的温度较低。此外，由于储热水箱底部存在温度分层情况，导致该区域温度梯度较大。

由图8还可看出，当θ不同时，储热水箱底部区域的温度也不同。当θ较大（60～90°）时，储热水箱底部的温度较低，为362.7～363.5 K，则储热水箱底部区域的温度梯度较大，当θ较小（15～45°）时，储热水箱底部的温度较高，为365.7～366.6 K，则储热水箱底部区域的温度梯度较小。

不同θ条件下，储热水箱中心截面（x=0）处速度场的分布情况如图9所示。

图9 不同θ条件下，储热水箱中心截面处（x=0）的速度场分布情况Fig.9 Distribution of the velocity in the tank at x=0 palne as θ varies

由图 9 可知，当 θ较小（15，30 °）时，水箱内水的流动速度自上而下（不包含底部区域）逐渐增大，这是由于当水从储热水箱顶部向下流动时，水箱顶部（y=5 m处）水的流速较小，为 0～0.015 m/s，随后在重力加速度的影响下，水的流动速度逐渐增大。当水流至水箱底部时，由于底面阻挡，速度有所降低。当θ由45°逐渐增大至90°时，水箱内水的流动速度整体上逐渐增大，此外，水箱内水的流动速度受涡旋影响较大，与高度之间的关系不明显。随着θ逐渐增大，储热水箱壁面对流动的阻碍能力逐渐减弱，则流体的流动速度大体上随倾角θ的增大而增大。

由图9还可看出，当θ为75，90°时，储热水箱底部区域 （y=-5 m附近）水的流动速度为0～0.005 m/s；当θ为15～60°时，储热水箱底部水的流动速度为0.006～0.033 m/s。总体来看，对于θ为75，90°的储热水箱，其底部区域水的流速明显小于其他倾角条件下的流速，这是由于当水流到储热水箱底面后，流动方向会发生改变，底面对流动产生堵塞作用，导致水的流速减小，并且这种阻碍能力会随着θ的增大而增强，当θ大于60°时，此影响尤为明显。这一现象也会减弱储热水箱内壁的自然对流换热强度。综上可知，改变θ，对储热水箱内水的流速分布情况影响较大，进而影响储热水箱内的对流换热情况。

4.3 储热水箱侧壁倾角对换热强度的影响

储热水箱的顶部、底部和侧壁的平均Nu数随θ的变化情况如图10所示。

图10 储热水箱的顶部、底部和侧壁的平均Nu数随θ的变化情况Fig.10 Variation of the average Nu number of the top，bottom and side walls of the hot water storage tank with θ

由图10可知，储热水箱顶部的平均Nu数随θ变化不大，说明θ的改变对储热水箱顶部的自然对流换热强度影响较小；当θ由15°逐渐增大至60°时，储热水箱侧壁的平均Nu数逐渐增大，当θ由60°逐渐增大至90°时，储热水箱侧壁的平均Nu数逐渐减小。由上文可知，θ越小，储热水箱侧壁自然对流换热强度越小，则储热水箱侧壁的平均Nu数越小，当其他条件相同时，储热水箱因自然对流产生的热损失就越小。随着θ逐渐增大，储热水箱侧壁的自然对流强度逐渐增强，导致储热水箱侧壁的平均Nu数逐渐增大，储热水箱的散热量也随之增加；当θ大于60°时，热水箱侧壁的平均Nu数有所减小，这是因为在有限空间的自然对流中，当流体流动到水箱底部时，流动方向改变，且θ越大，流向改变得越明显，水箱底面对流动的堵塞作用越大，流体的流速随之减小，导致水箱底面的自然对流换热强度变弱。

储热水箱边界面的平均Nu数随θ的变化情况如图11所示。

将水箱侧壁、底面、顶面的Nu数求平均值，得到储热水箱整体的平均Nu数。由图11可知，当θ由0°逐渐增大至60°时，储热水箱边界面的平均Nu数逐渐增大，这是由于随着θ逐渐增大，水箱侧壁处的自然对流换热强度逐渐增强，导致储热水箱边界面的平均Nu数随之增大；当θ由60°逐渐增大至90°时，由于水箱底面对流动的阻碍能力越发明显，因此水箱底面的自然对流换热强度逐渐减小，导致储热水箱边界面的平均Nu数随之减小。

图11 储热水箱边界面的平均Nu数随θ的变化情况Fig11 Variation of the average Nu number of the boundary surface of the hot water storage tank with θ

4.4 侧壁倾角对储热水箱散热过程的影响

不同的θ条件下，储热水箱平均温度随时间的变化情况如图12所示。

图12 不同的θ条件下，储热水箱平均温度随时间的变化情况Fig.12 Changes of the average temperature of the water tank with time for different θ

由图12可知，当θ为75°时，储热水箱的散热速度最快；当θ为45°时，储热水箱的散热速度最慢。由此可以看出，不同侧壁倾角条件下，储热水箱的散热量与其Nu数的变化规律不完全相符，即散热量不仅与换热强弱有关，可能还与其比表面积相关。

储热水箱的比表面积为其表面积与体积之间的比值。储热水箱比表面积S的计算式为

式中：S′为储热水箱的表面积；V为储热水箱的体积。

储热水箱比表面积S随θ的变化情况如图13所示。

图13 储热水箱比表面积S随θ的变化情况Fig.13 The specific surface area S of the tank changes with θ

由图13可知，随着θ逐渐增大，储热水箱比表面积S呈现出先减小再增大的变化趋势。当θ为60°时，储热水箱的比表面积达到最小值；当θ为15°时，储热水箱的比表面积达到最大值，这与储热水箱各面及整体的平均Nu数的变化趋势恰好相反。

综上可知，储热水箱的比表面和对流换热强度对自身散热情况的影响结果：当θ为15°时，虽然储热水箱内的对流换热强度较弱，但比表面积较大，导致储热水箱单位时间的散热量较大，平均温度随时间的下降速率较大，可见，此时比表面积对储热水箱的散热量影响较大，40 min时，储热水箱的平均温度降至367.1 K；当θ为45°时，储热水箱的比表面积显著减小，并接近于最小值，并且其对流换热强度相对较弱，因此，储热水箱温度下降得比较缓慢，当τ为40 min时，储热水箱平均温度为367.4 K。；当θ为60°时，虽然储热水箱的比表面积较小，但是其对流换热强度较强，所以储热水箱温度下降得比较迅速，此时对流换热强度对储热水箱散热量的影响较大，当τ为40 min时，储热水箱的平均温度降低至367 K。总体来看，在对流换热强度和比表面积的双重作用下，当θ为45°时，储热水箱的温度下降得最慢，散热量也最少。

4.5 比表面积对散热的影响

基于不同θ下，储热水箱温度的变化情况可以看出，比表面积对储热水箱的散热强弱有一定的影响。下面进一步分析储热水箱的比表面积与其平均温度之间的关系。

本文分别以边长为 10，12，15，20，25，30 m 的正方体水箱（θ为90°）为研究对象，对应的比表面积分别为 0.6，0.5，0.4，0.3，0.24，0.2 m-1。水的初始温度T0为368 K，散热时间为20 min。

储热水箱的比表面积与平均温度的关系如图14所示。

图14 储热水箱的比表面积与平均温度的关系Fig.14 Changes of the average temperature of the water tank with time for the cubic shape tank

由图14可知，比表面积越大，储热水箱的平均温度越低，这是由于比表面积越大，单位体积储热水箱与外界之间的接触面积越大，即换热面积越大，因此储热水箱的换热速率越大，最终导致储热水箱的平均温度越低。

5 结论

本文基于数值模拟结果分析了埋于土壤中的储热水箱的侧壁倾角θ对其内部自然对流换热强度的影响，分析结论如下。

①储热水箱侧壁倾角越小，侧壁对水箱内部流体流动的阻碍能力越大，流体流速越小，从而影响水箱内部的自然对流换热强度：当储热水箱侧壁倾角由15°逐渐增大至60°时，该水箱的自然对流换热强度逐渐增强，当储热水箱侧壁倾角由60°逐渐增大至90°时，由于水箱底面对流动的阻碍作用愈发显著，使得该水箱内的对流换热强度开始缓慢减弱。

②当储热水箱侧壁倾角大于60°时，水箱底面对流体流到的阻碍作用较大，从而导致水箱内的换热强度减弱。

③随着比表面积逐渐大，储热水箱的散热强度逐渐增大。

④由自然对流换热强度和比表面积对储热水箱散热量的影响可知，当储热水箱侧壁倾角为45°时，其内部水温的下降速率较慢，散热量较小。

⑤储热水箱侧壁倾角较大时，水箱底部的温度分层比较明显，需要分层取水时，应使得储热水箱的侧壁倾角较大。