程素雅， 陈宝明、2， 郭梦雪， 张艳勇， 李佳阳

(1.山东建筑大学热能工程学院，山东济南250101；2.山东省建筑节能技术重点实验室，山东济南250101)

1 概述

相变材料(PCM)以其储能密度高、相变过程温度几乎恒定等优点，在储能领域得到了广泛的应用[1-2]，如太阳能储存[3]、建筑节能[4-5]、工业余热回收[6]、纺织品制造[7]等。根据材料的性质，相变材料可分为金属材料、非金属材料[8]。在实际使用中，与非金属相变材料(如石蜡、熔盐)相比，金属相变材料(如锡、铝)的投资更高。然而导热能力低的特点，限制了非金属相变材料的应用。为此，许多学者开发了增强技术[9-10]来克服这个缺点。其中，金属翅片作为一种简单、可靠、有效的方法，被广泛应用于非金属相变材料的应用中(在容器内部加入翅片，促进腔体内相变材料相变过程更加均匀)[11-15]。

相变材料在带翅片矩形腔内的熔化问题得到了广泛的研究。文献[8]在带有水平翅片的透明矩形腔中进行了相变材料熔化实验。实验结果表明，随着翅片数量的增加，热分层区域变小，总传热速率增大，相变材料完全熔化时间缩短。文献[16]通过搭建固液相变储热实验系统，记录相变材料内实时温度响应数据，研究添加环形翅片蓄热装置的蓄热性能。文献[17]实验研究了基于相变材料的散热器在热管理中的应用，研究了翅片数量、翅片高度、翅片厚度等参数的影响。结果表明，增加翅片的数量和高度，翅片的整体热性能明显提高，而增加翅片厚度的效果并不明显。

与实验方法相比，数值模拟可以提供更详细的传热和相变过程，特别是对于自然对流传热效应。文献[18]建立了数值模型用于模拟带翅片的金属外壳内相变材料熔化行为，分析了翅片数量、翅片长度和厚度、壁面温度等对熔化过程的影响。文献[19]模拟了翅片布置角度对熔化过程的影响，结果表明翅片布置角度影响显著，对于熔化速率的影响程度，翅片向下倾斜优于向上倾斜。文献[20]对翅片长度和间距的影响进行了数值模拟，结果表明，相比于翅片间距，翅片长度对熔化过程的影响有限。针对熔化过程，文献[21]对翅片布置方式进行了数值模拟，发现将翅片布置在蓄热单元下部有利于提高温度分布均匀性。

基于文献调研发现，翅片的长度、数量、厚度、角度、翅片间距等参数是研究的重点，但对于翅片排布方式(主要是长度的排布)的研究比较少，特别是矩形腔内翅片排布方式对相变材料的相变过程和自然对流传热的影响。因此，本文以填充石蜡的矩形腔(分为无翅片矩形腔、带翅片矩形腔)为研究对象，建立数学模型。采用有限元软件COMSOL Multiphysics模拟矩形腔内石蜡的熔化行为，分析不同翅片排布方式对石蜡熔化行为的影响，筛选有利于增强石蜡熔化的翅片排布方式。

2 模型建立

2.1 问题描述

将矩形腔简化为二维几何模型进行数值模拟，见图1，图中尺寸单位为mm，x轴正方向为垂直向上。腔体高度为120 mm，长度为50 mm，3个平行金属翅片垂直腔体右侧排布，长度分别为L1、L2、L3，翅片厚度为3 mm，间距均为27.75 mm。除右侧壁面为恒温受热面(温度恒定为313 K)外，其余壁面均包裹足够厚度的保温材料，视为绝热。相变材料选用石蜡，金属翅片材料为铝。石蜡、铝的物性参数分别见表1、2。

图1 二维几何模型

表1 石蜡物性参数

表2 铝物性参数

2.2 翅片排布方案

在翅片总长为96 mm不变的前提下，设置了4种翅片排布方式，见表3。

表3 翅片排布方式

2.3 数学模型

我们使用焓-孔隙率法来模拟石蜡熔化过程，该方法将整个计算域视为一个多孔域，每个单元的孔隙率由液相率表征。为简化计算，对数值模型进行必要的设定：相变材料的物性参数在固、液相中不随温度发生变化，固液共存时随温度线性变化。液相为牛顿不可压缩流体并且符合Boussinesq假设，即忽略压强变化对相变材料密度的影响，仅考虑温度变化引起的相变材料密度变化。液相状态下相变材料流动为非稳态层流流动。忽略恒温受热面的壁厚及热阻。压力设定为标准大气压力。基于以上设定，采用以下控制方程对模型进行计算。

连续性方程：

式中ρ——石蜡的密度，kg/m3

t——时间，s

u——Ox轴方向速度，m/s

x——Ox轴坐标，m

w——Oy轴方向速度，m/s

y——Oy轴坐标，m

动量方程：

x方向：

式中p——压力(绝对压力)，Pa

Su——速度u方向动量方程源项

f——液相率

ε——系数(为防止分母为0)，取10-3

σ——固液相混合物的连续性系数，取5×104

液相率f的计算式为：

f=0，T≤Ts

f=1，TL≤T

式中T——任意时刻石蜡温度，K

Ts——熔化开始温度，K

TL——熔化终止温度，K

y方向：

式中μ——石蜡的动力黏度，Pa·s，本文取4.43×10-3Pa·s

Sw——速度w方向动量方程源项

ρref——石蜡的基准密度(即石蜡的初始密度)，取900 kg/m3

g——重力加速度，m/s2，本文取9.8 m/s2，α——石蜡体膨胀系数，K-1，本文取0.001 K-1

Tref——石蜡的基准温度，即石蜡的初始温度，K，本文取293 K

能量方程：

h=hx+Δh

Δh=fQq

式中h——任意时刻石蜡比焓，J/kg

λ——石蜡的热导率，W/(m·K)

cp——比定压热容，J/(kg·K)

Sh——能量方程源项

hx——显热比焓，J/kg

Δh——潜热比焓，J/kg

href——石蜡基准比焓(即初始比焓)，J/kg

Qq——石蜡相变潜热，J/kg

2.4 初始条件和边界条件

初始条件：石蜡初始温度为293 K，相变开始之前石蜡为固态。受热面温度恒定为313 K。

边界条件：除受热面外，其他面均为绝热面。

2.5 网格无关性检验

采用有限元软件COMSOL Multiphysics通过耦合层流与流固传热物理场，求解上述控制方程。在层流物理场中，设置液态石蜡物性参数，添加浮升力及重力。在流固传热物理场中，设置固态石蜡物性参数及边界条件。最后设置多物理场耦合进行求解。

在数值计算中，网格质量直接影响计算精度。为了能够尽量精确并节省计算资源，在选取合适的模型参数的基础上，首先对网格进行无关性检验，确定计算所需网格数。笔者针对翅片排布方式1，以矩形腔内石蜡的液相面积比(矩形腔内，液相石蜡面积与矩形面积之比)为研究对象，进行网格无关性检验。依次选取网格数为3 890、6 679、10 076、21 289进行网格检验。由受热时间为50、400、800、1 500 s的液相面积比模拟结果可知，在相同受热时间下，网格数量为10 076、21 289的模拟结果基本一致。因此，认为当网格数量达到10 076时，模型计算精度可满足要求。网格划分后的二维模型见图2。

图2 网格划分后的二维模型

3 模拟结果与分析

3.1 无翅片矩形腔内石蜡熔化行为

受热时间为300、1 100、2 700 s时，无翅片矩形腔内石蜡的液相率分布分别见图3～5。文中出现的矩形腔内石蜡液相率分布图与图3共用标值。由图3～5可知，右上角的石蜡最先熔化，然后熔化部分向矩形腔中心扩散，直至矩形腔左下角石蜡完成熔化。

与受热面接触的石蜡受热首先熔化，高温液态石蜡在自然对流传热的作用下向上流动，而矩形腔下部的石蜡温度仍低于熔化开始温度，处于固态。随着受热时间的延长，在导热和自然对流传热的共同作用下，矩形腔底部石蜡得以熔化。因此，自然对流传热导致了矩形腔内石蜡非均匀熔化。

图3 受热时间为300 s时无翅片矩形腔内石蜡的液相率分布

图4 受热时间为1 100 s时无翅片矩形腔内石蜡的液相率分布

图5 受热时间为2 700 s时无翅片矩形腔内石蜡的液相率分布

3.2 带翅片矩形腔内石蜡熔化行为

① 液相率分布

受热时间为80、300、600、1 100 s时，4种翅片排布方式矩形腔内石蜡的液相率分布分别见图6～9。由图6～9可知，在熔化初期，相界面(液相石蜡与固相石蜡的过渡区)几乎与热壁面(受热面和翅片表面)平行。这表明在熔化初期，导热占主导地位。随着受热时间延长，液相石蜡逐渐增多，相界面不再平行于热壁面，说明自然对流传热作用逐渐显现。

由图4、9可知，与无翅片矩形腔内石蜡熔化行为相比，在受热时间为1 100 s时，带翅片矩形腔内石蜡基本完全熔化。因此，矩形腔增加翅片可有效改善石蜡熔化的均匀性，在一定程度上缩短了矩形腔内石蜡的熔化时间。

由图9可知，在受热时间为1 100 s时，采用排布方式2的带翅片矩形腔内石蜡的熔化面积最大。这说明，在4种翅片排布方式中，排布方式2对改善矩形腔内石蜡熔化均匀性的效果最理想，石蜡完成熔化的时间最短。

图6 受热时间为80 s时4种翅片排布方式矩形腔内石蜡的液相率分布

图7 受热时间为300 s时4种翅片排布方式矩形腔内石蜡的液相率分布

图8 受热时间为600 s时4种翅片排布方式矩形腔内石蜡的液相率分布

图9 受热时间为1 100 s时4种翅片排布方式矩形腔内石蜡的液相率分布

② 液相面积比与完全熔化受热时间

矩形腔内石蜡的液相面积比随受热时间的变化见图10。由图10可知，相同受热时间下，带翅片矩形腔内石蜡的液相面积比明显高于无翅片矩形腔。无翅片矩形腔内石蜡完全熔化受热时间为3 522 s，带翅片矩形腔翅片排布方式1～4的石蜡完全熔化受热时间分别为1 874、1 674、2 082、1 910 s。

图10 矩形腔内石蜡的液相面积比随受热时间的变化

③ 熔化增强率

将传统的等长翅片的排布方式1作为基准，采用熔化增强率评价不同受热时间其他3种翅片排布方式对矩形腔内石蜡熔化的增强作用。某受热时间的熔化增强率为正，说明该受热时间下，采用该翅片排布方式时的矩形腔内石蜡熔化面积比采用等长度翅片排布方式更大，熔化效果更好。反之，说明矩形腔内石蜡熔化面积更小，熔化效果逊色。

第i(i=2～4)种排布方式某受热时间的熔化增强率ηi的计算式为：

式中ηi——第i种排布方式某受热时间的熔化增强率

Ai——第i种排布方式该受热时间的液相面积比

A1——排布方式1该受热时间的液相面积比

排布方式2～4不同受热时间的熔化增强率见图11。分析图11可知，翅片排布方式2对矩形腔内石蜡熔化的增强作用明显，增强作用集中在熔化过程的中后期(受热时间400 s至石蜡完全熔化结束)。与排布方式1相比，排布方式3、4起到了相反作用。

图11 排布方式2～4不同受热时间的熔化增强率

4 结论

① 对于无翅片矩形腔，在自然对流传热作用下，右上角的石蜡最先熔化，然后熔化部分向矩形腔中心扩散，直至矩形腔左下角石蜡完全熔化。矩形腔增加翅片可有效改善石蜡熔化的均匀性，缩短了矩形腔内石蜡的熔化时间。

② 排布方式2对改善矩形腔内石蜡熔化均匀性的效果最理想，石蜡完成熔化的时间最短。相同受热时间下，带翅片矩形腔内石蜡的液相面积比(液相石蜡面积与矩形面积之比)明显高于无翅片矩形腔。

③ 无翅片矩形腔内石蜡完全熔化的受热时间为3 522 s，带翅片矩形腔翅片排布方式1～4的石蜡完全熔化的受热时间分别为1 874、1 674、2 082、1 910 s。

④ 翅片排布方式2对矩形腔内石蜡熔化的增强作用明显，增强作用集中在熔化过程的中后期。排布方式3、4起到了相反作用。