管壳式相变储热单元偏心布置强化换热数值研究

Numerical Study on Enhanced Heat Transfer of Shell-and-Tube Latent Heat Storage Units with Eccentric Arrangement

摘要: 目前，相变储热技术得益于储热密度高，温度变化稳定以及安全可靠等优势，展现出广阔的应用前景。然而，该技术仍具有传热速率较慢的问题，严重制约了其工程化应用与性能提升。为实现强化换热的目标，提升相变储热单元的传热速率，本文通过等效热容法模拟相变储热单元熔化过程，首先设置规则肋片以探究对于相变材料熔化过程的具体影响。其次再进行偏心布置对于相变储热单元熔化过程的具体影响。结果表明，肋片可以有效提升相变储热单元熔化速率，偏心布置将进一步强化换热。这些研究内容进一步提升相变储热系统的综合性能，为相变储热单元的高效利用与结构优化设计奠定了理论基础。

Abstract: At present, latent heat thermal energy storage (LHTES) technology has demonstrated broad application prospects due to its advantages of high energy storage density, stable operating temperature, and reliable safety performance. However, the relatively slow heat transfer rate remains a major limitation, severely restricting its engineering applications and further performance improvements. To achieve enhanced heat transfer and improve the thermal response rate of LHTES units, this study employs the effective heat capacity method to numerically simulate the melting process within a shell-and-tube latent heat storage unit. Initially, conventional fin configurations are introduced to investigate their specific effects on the melting characteristics of phase change materials (PCMs). Subsequently, the influence of eccentric tube arrangements on the melting process is further analyzed. The results indicate that fins can significantly accelerate the melting rate, and eccentric arrangements further enhance the overall heat transfer performance. These findings not only improve the comprehensive performance of LHTES systems but also provide a theoretical foundation for the efficient utilization and structural optimization of latent heat storage units.

1. 引言

能源是经济社会发展的核心驱动力，其稳定供应对全球工业、交通、建筑以及居民生活具有至关重要的意义。然而，随着全球能源需求的持续攀升，传统化石能源的利用面临严峻挑战。为应对能源危机和环境问题，各国纷纷加快能源结构调整，推动可再生能源的发展。由于可再生能源具有间歇性与波动性特性，其大规模应用仍然受到储热技术的限制[1]。其中，相变储热凭借较高的能量存储密度、稳定的温度调节能力和较长的循环寿命，在可再生能源存储、工业余热回收和建筑节能等领域表现出显著优势[2]。但相变材料(PCM)的导热系数较低，导致热量存储和释放过程缓慢，从而制约了其工程应用[3]，使储热和放热过程受到阻碍。因此，提升相变储热单元的传热性能已成为当前研究关注的核心问题。而增设肋片可以有效强化储热单元换热能力。本质上，肋片的作用在于增加换热面积，进而增强单元内部热量传递速率。

Guo等人[4]提出一种新型肋片优化方法，即在非均匀排列条件下结合不同结构的肋片。研究表明，凹面肋片在初始熔化阶段强化对流换热效果好，凸面肋片最差，直肋片在熔化后期维持对流换热能力较强。而均匀排列时，凹面肋片位置对熔化过程影响小，凸面肋片位置影响大，组合肋片相比直肋片无明显优势。可知肋片的布置与结构存在一定要求，并不是任意肋片都会强化传热。Zhu等人[5]提出河内塔形肋片分布方式，通过数值模拟研究其对相变储热单元熔化过程的影响，而河内塔形肋片分布的间距呈等差数列递增，肋片长度按固定比例增加。研究通过数值模拟比较了无肋片、均匀肋片分布和河内塔形肋片分布对熔化过程的影响。结果表明，河内塔形肋片分布显著提升了相变材料的熔化速率和均匀性，总熔化时间缩短了18.04%。优化后，底部间距17.82 mm、间距公差3.38 mm、长度比1.2的肋片参数组合使熔化时间进一步缩短。这种分布方式在保持热储存能力的同时，有效缩短了相变时间，提高了温度均匀性。Ao等人[6]设计19种不同V形肋片结构，研究肋片结构参数对PCM熔化性能的影响，为提高管壳式相变储热单元性能提供了理论依据和设计参考。结果表明，研究发现，自然对流对PCM熔化速率有显著影响，底部熔化时间占总熔化时间的62.5%。优化后的肋片结构参数为：肋片角度55˚、长度9 mm、数量5、分支角度55˚、分支长度3 mm。与无肋片相比，优化后的肋片结构使热储存效率提高了2.57倍，热储存时间缩短了61.05%。

综上所述，增设肋片是强化换热研究领域的重点之一，本文将研究偏心布置对管壳式相变储热单元熔化过程的影响，以此总结规律为后续研究提供经验。

2. 模型建立

2.1. 物理模型

建立二维相变传热单元如图1所示，布置8个相同的铜金属肋片等角度分布。其中储热单元外径30 mm，内径7.5 mm，内管壁厚1 mm。规则肋片的长度为16 mm，肋片宽2 mm。进行相变传热模拟过程中，模型的初始温度为20℃，热源边界80℃。

Figure 1. Diagram of basic model

图1. 基础模型图

2.2. 数学模型及材料物性参数

为了进行固–液相变传热模拟，在COMSOL Multiphysics采用等效热容法计算模拟相变传热过程[7]，进行以下假设：

(1) 肋片与PCM材料为各向同性材料；

(2) PCM考虑自然对流，密度采取Boussinesq假设[8]，在动量方程中添加体积力，引起由密度变化而造成的自然对流；

(3) PCM熔化过程被定义为低速、非稳态过程和不可压缩流动过程，流体为层流状态。

本文计算相变传热利用等效热容法，将相变潜热定义为相变过程比热容的变化，与焓法不同的是，等效热容法待求变量只有温度。先给出计算的相关方程。

连续性方程：

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ (1)

式中，u与v分别为x方向和y方向的速度分量。

动量守恒方程：

x方向：

$\frac{\partial u}{\partial t}+\left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\mu \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)-S_t\left(T\right)\cdot u$ (2)

y方向：

$\frac{\partial v}{\partial t}+\left(u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\right)=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\mu \left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)+F_B-S_t\left(T\right)\cdot v$ (3)

式中，μ表示动力粘度；单位为kg/(m·s)。

考虑自然对流由密度差引起，引用Boussinesq近似，在y方向动量方程引入浮力项F_B，表达式：

$F_B={\rho }_l\alpha g\left(T-T_m\right)$ (4)

式中：α为热膨胀系数，是描述浮力相的关键参数，ρ_l为液相PCM密度，单位为kg/m³。

$S_t\left(T\right)=A_m\frac{{\left(1-\phi \left(T\right)\right)}^2}{\phi {\left(T\right)}^3+\epsilon }$ (5)

将阻力项S_t(T)引入到动量方程，S_t(T)为Kozeny-Carman [9]源项，用于区分动量方程中PCM内的固液相。固体区域中的阻力较大，液体区域阻力较小，从而确保固相和液相之间的平滑过渡。其中ε是一个非常小的数字，以避免被零除，而A_m是糊状区域的常数，通常取值范围为10⁴~10⁷之间，本文取10⁷作数值研究，单位为kg/(m³·s)。当PCM为固体，φ(T)取0值，使得S_t(T)数值较大，而当PCM为液体，φ(T)取1值，S_t(T)数值较小，当阻力阻碍了流动，可视为固体存在。

区分PCM是处于固相或是液相，判断流动条件：

$\phi \left(T\right)=\{\begin{array}{ll}0\hfill & \left(T\le T_s\right)\hfill \\ \frac{T-T_s}{T_l-T_s}\hfill & \left(T_s<T<T_l\right)\hfill \\ 1\hfill & \left(T\ge T_l\right)\hfill \end{array}$ (6)

式中：T_l为液相温度，而T_s为固相温度。

能量守恒方程：

$\rho C_p\frac{\partial T}{\partial t}+\rho C_p\left(u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}\right)-\frac{\partial }{\partial x}\left[\lambda \left(T\right)\frac{\partial T}{\partial x}\right]-\frac{\partial }{\partial y}\left[\lambda \left(T\right)\frac{\partial T}{\partial y}\right]=0$ (7)

式中：λ为导热系数，单位为W/(m·K)，C_p为比热容，单位J/(kg·K)。

同时，PCM的热物理性质随液相分数变化而变化，用以区分不同状态下材料的具体参数，密度ρ(T)和热导率λ(T)的表达式如下：

$\rho \left(T\right)={\rho }_s+\left({\rho }_l-{\rho }_s\right)\cdot \phi \left(T\right)$ (8)

$\lambda \left(T\right)={\lambda }_s+\left({\lambda }_l-{\lambda }_s\right)\cdot \phi \left(T\right)$ (9)

由于固态和液相之间的粘度存在显著差异，为了区分固态和液相相变材料，液相材料粘度μ(T)的表达式如下：

$\mu \left(T\right)={\mu }_l+S_t\left(T\right)\cdot {\mu }_l$ (10)

等效热容法用于计算相变材料在相变过程中吸收和释放的热量，同时考虑相变材料本身的潜热。表达式如下：

$C_p\left(T\right)=C_{p,s}+\phi \left(T\right)\cdot \left(C_{p,l}-C_{p,s}\right)+L\cdot D\left(T\right)$ (11)

式中：L为PCM潜热值，单位kJ/kg。

其中，公式(11)使用高斯函数来校正比热容。

$D\left(T\right)=\frac{{\text{e}}^{-\frac{{\left(T-T_m\right)}^2}{{\left(\Delta T/4\right)}^2}}}{\sqrt{\pi \cdot {\left(\Delta T/4\right)}^2}}$ (12)

式中：ΔT为相变区间，e为自然常数。

初始条件与边界条件：

$T\left(x,y\right)=T_0{r}^2\le x^2+y^2\le R^2$ (13)

${T|}_{\Omega }=T_w\Omega :x^2+y^2=r^2$ (14)

$-\lambda {\frac{\partial T}{\partial n}|}_{\zeta }=0\zeta :x^2+y^2=R^2$ (15)

其中T₀为储热单元初始温度；Ω为热源位置，温度T_w为热壁面温度。在熔化过程中，初始温度为20℃，热源温度为80℃，ζ为外壁面，外壁面为绝热条件。

其中，肋片与PCM接触区域传热视为连续，由此表示：

${T_{PCM}|}_{\psi }={T_{fin}|}_{\psi }$ (16)

$-{\lambda }_{PCM}{\frac{\partial T_{PCM}}{\partial n}|}_{\psi }=-{\lambda }_{fin}{\frac{\partial T_{fin}}{\partial n}|}_{\psi }$ (17)

式中：$\psi$ 表示PCM与肋片交界处。

本文所述PCM与铜物性参数如表1所示。

Table 1. Material property parameters

表1. 材料物性参数

参数	相变材料	铜
密度/kg/m^3	771.2	8920
比热容/J/(kg·K)	2176	381
相变潜热/J/kg	232,400
相变温度/℃	55~57
热导率/W/(m·K)	0.089	387.6
热膨胀系数/K^−^1	0.00075
动力粘度/kg/(m·s)	0.00356

3. 网格无关性验证及模型验证

物理模型

为提高相变传热瞬态模拟的数值稳定性与计算效率，本文在计算过程中采用自适应步长策略，计算步长可根据收敛性自动调整。如图2所示，分别选取0.5 s、0.8 s、1 s作为最大时间步长约束，对相变过程进行模拟。模拟结果表明，不同最大时间步长条件下，液相分数随时间变化趋势基本一致，对熔化过程的影响较小。因此，后续模拟中不再对最大步长施加限制。

为选择合适的网格数来模拟相变传热，本文选择28,000，46,000，64,000的3种网格数来分别进行模拟，如图3所示。结果表明，不同最大时间步长对于熔化过程液相分数变化影响很小，故不考虑最大步长限制。而网格数为46,000与网格数为64,000结果差距很小，为节省计算资源，选择46,000网格数作为计算要求。

Figure 2. Diagram of the mesh count validation

图2. 网格校核图

Figure 3. Diagram of time step validation

图3. 时间步长校核图

为确保相变传热模拟计算的准确性，本文根据Al-Abidi等人[10]实验，重新建立实验中的模型，研究人员设置石蜡的初始温度为27℃，热流体温度为90℃，选取15个实验监测点表征PCM的平均温度变化，本文利用COMSOL Multiphysics模拟相变传热，构建相同的物理模型，使用有限元数值计算方法，求解器采用PARDISO算法，对计算域采用二阶离散格式，在时间离散方面，采用隐式格式与向后差分公式，求解器将根据收敛性自动调整步长。为更好地进行数值模拟，选择三角形网格与网格数为68,644的条件进行仿真。结果如图4所示，可以看出，数值模拟温度变化趋势与实验基本一致，其中最大误差为7.5%，平均误差为2.5%，这主要是因为实验未能完全达到外部绝热条件，影响了温度的变化。综上所述，模拟结果与实验监测温度重合度较好，证明数值模拟的可行性，为后续工作提供依据。

Figure 4. Diagram of numerical methods assessment

图4. 数值方法考核图

4. 结果分析

为全面评估肋片结构与布置方式对相变过程的影响，进一步引入偏心结构，探讨偏心位置对拓扑优化过程中肋片生长特性及熔化行为的作用机理。图5为偏心结构示意图，考虑不同的方向，不同的距离来比较熔化过程。方向考虑竖直向上与竖直向下，偏心距离从2 mm至12 mm，移动步长为2 mm。

Figure 5. Diagram of eccentric arrangement

图5. 偏心布置示意图

4.1. 储热单元向上偏行布置对熔化过程的影响

图6为单管向上偏心单元熔化液相分数云图。熔化初始阶段，热量通过肋片传递到PCM，贴近肋片的PCM率先熔化，熔化进行到500 s时，肋片周围PCM已完全熔化。由于自然对流的作用，相变储热单元上方热量传递进一步加强，肋片的向上布置使得储热单元上半部分熔化速率较快，当整体模型熔化进行到3500 s时，整体PCM已熔化接近90%，但底部区域仍有少量的PCM未熔化。这是由于自然对流的发展使得热量积聚于上半部分，下半部分热量传递效果较差，一部分通过下方肋片热传导传递热量，另一部分通过上半部分热量随着PCM流动缓慢传递到底部。其中，当熔化进行到6500 s时，向上偏心2 mm于4 mm模型熔化过程将要结束，这可以看出在向上偏移的过程中，对于熔化的提升不再明显。

Figure 6. Diagram of the upward eccentric model melting

图6. 向上偏心模型熔化云图

图7为向上偏心模型熔化液相分数随时间变化图。各个模型的熔化时间分别为96.9 min、128.1 min、149.5 min、152.8 min、156.1 min和159.1 min。熔化时间最大减少39.09%。在熔化前期，相变储热单元以热传导为主要传热机制，液相曲线斜率较大。随着熔化过程推进，当熔化进行到70%，基础向下偏心2 mm模型的熔化曲线斜率开始明显下降，熔化速率减慢。这是由于熔化后期底部区域PCM所受自然对流作用有限，换热能力下降，导致局部熔化速率减慢。而其他模型的熔化曲线斜率降低时间点均在70%之前，这意味着底部区域影响更大，整体熔化时间增大。在所有模型中，向上偏心2 mm模型的熔化时间最短。但难以有效改善底部熔化速率缓慢的问题，所以向上偏心布置熔化时间延长。图8为向上偏心模型熔化过程平均速度图。从图中可以看出，向上偏心2 mm模型的平均速度最大为2.1 × 10⁻⁴ m/s，最大提升了60%。且随着向上偏心距离的增加，平均速度大小直线下降，证明向上偏心布置的设置会让熔化过程中PCM流动受阻，使得自然对流作用减弱，这也证明了向上偏心对于熔化过程存在阻碍作用。

4.2. 储热单元向下偏心布置对熔化过程的影响

图9为单管向上偏心模型熔化液相云图。熔化初始阶段，热量通过肋片传递到PCM，贴近肋片的PCM率先熔化，熔化进行到100 s时，肋片周围PCM已完全熔化。当熔化进行到500 s时，熔化产生的液体PCM已初具规模，自然对流开始发挥作用。随着熔化的进行，熔化界面快速变化，沿着肋片的方向延伸与扩张。当熔化进行到1000 s时，可以观察到各个模型熔化过程的明显差异，当向下偏心距离较短时，对于整体熔化界面的影响较小，熔化界面相似。当向下偏心距离增大时，底部肋片延伸至储热单元底部，热量通过底部肋片更快地传递到下半部分单元，熔化界面表现为下半部分PCM已经完全熔化，与单元边界接壤，向下偏心12 mm模型表现最为明显。当熔化进行到最终时刻，由于向下的布置，使得最终熔化区域集中在顶部，而顶部区域受到自然对流作用熔化速率保持稳定。结果表明，向下偏心布置可以有效改善熔化过程，强化对流换热。

Figure 7. Diagram of the upward model melting liquid fraction

图7. 向上模型熔化液相分数图

Figure 8. Diagram of the upward model time-averaged velocity distribution

图8. 向上模型流动平均速度图

Figure 9. Diagram of the download eccentric model melting

图9. 向下熔化液相云图

图10为基础向下偏心模型熔化液相分数随时间变化图。各个模型的熔化时间分别为48 min、34 min、36 min、37 min、37 min和38 min。熔化时间最大减少29%。其中各个模型在熔化进行到50%之前，熔化过程趋于一致，随着熔化的进行，各个熔化曲线逐渐分散，证明各个模型的熔化速率发生变化，这也能证明向下偏心的设置较大程度影响了中后期熔化过程。当熔化进行到90%时，各个模型的熔化差异增大，除去向下偏心2 mm模型，其余模型的熔化曲线较为集中，代表向下偏心布置对于熔化的提升差距不大。

Figure 10. Diagram of the download model melting liquid fraction

图10. 向下模型熔化液相分数图

图11为向下偏心模型熔化过程平均速度图。从图中可以看出，向下偏心12 mm模型的平均速度最大为5.3 × 10⁻⁴ m/s，最大提升了49%。且随着向下偏心距离的增加，平均速度大小直线上升，证明向下偏心布置会让熔化过程中PCM流动得到强化，使得自然对流作用增强。相对于向上偏心模型的速度整体得到提升，两者速度比较最大提升61%，证明向下偏心布置的优势。

Figure 11. Diagram of the download model time-averaged velocity distribution

图11. 向下模型流动平均速度图

5. 结论

本文通过数值模拟方法，对几个偏心模型进行相变传热分析，得出以下结论：向上偏心模型在熔化过程中热性能表现差距较大，随着向上偏心距离的增加，阻碍了底部传热效果，延长了整体的熔化时间，而向下偏心模型相对于向上偏心模型在熔化过程有明显提升，最大提升78%，各个向下偏心模型熔化时间差距很小，表明向下偏心设置会强化换热效果。