基于离散相模型的相变微胶囊流体传热特性数值模拟

吴兴辉，杨震，陈颖，段远源

（1 清华大学热科学与动力工程教育部重点实验室，北京市二氧化碳利用与减排技术重点实验室，北京100084；2广东工业大学广东省功能软凝聚态物质重点实验室，广东广州510006）

引 言

相变材料在物态变化时吸收/释放大量热量，可用于储能及传热强化，已被应用于建筑[1-2]、能源[3-4]、微电子[5-6]等领域。石蜡作为相变材料由于潜热较大、价格便宜、过冷度小等优点受到广泛关注[7]。然而由于相变材料热导率低、融化后不定形，限制了其工程应用。将相变材料包裹在微胶囊中，可使其在融化后保持所需形状，并可提高比表面积来强化传热[8]。将相变微胶囊颗粒添加至单相热流体中形成相变微胶囊悬浮液（micro-encapsulated phase change material slurry, MPCS）（图1）用作传热流体，在强化传热方面具有发展潜力[9]。

微胶囊包裹技术运用成膜材料将固体、液体或气体包覆成具有核壳结构颗粒，所得颗粒称为微胶囊。微胶囊粒径通常在2～2000 μm, 外壳厚度0.5～150 μm，多为球形。1953 年，美国NCR 公司将微胶囊技术应用于无碳复写纸，标志着高聚物外壁微胶囊开始工业运用，微胶囊技术迅猛发展[10]，已广泛应用于储能[11-12]、化工[13-14]、医药[15-16]、农业[17]等领域。

MPCS可用于传热流体[18]，相比于传统纯流体具备诸多优势：（1）与基液相比，悬浮液温度在传热过程几乎不变，保持较大的传热温差；（2）相同传热温差下，石蜡胶囊悬浮液比热容随石蜡质量分数增大而增大，当质量分数为40%、60%时分别比基液增大约20、30倍；（3）相同传热量下需要的工质质量流量低，泵功减少；（4）添加相变微胶囊颗粒可以增大基液的热导率，强化传热。

学者们已对MPCS 流动和传热特性开展了实验研究，Goel 等[19]实验研究了具有恒定热通量的圆管中MPCS 的层流受迫对流换热，分析了微胶囊体积分数、Stefan 数等对传热性能的影响，结果表明与水相比，采用MPCS 后壁温降低了约50%；Chen 等[20]开展了恒定热通量的圆管中MPCS 的层流受迫对流换热实验研究，发现当传热速率为750 W，微胶囊质量分数为15.8%时，MPCS比纯水功耗降低67.5%；Choi等[21]对MPCS 湍流中对流传热特性开展了实验研究，提出了预测沿管温度分布的三区融化模型，研究结果表明MPCS 传热系数在温度低于熔点时增加，在温度高于熔点后快速下降。

对MPCS 传热的数值模拟研究，Charunyakorn等[22]利用隐式差分法模拟恒热流下圆管内流动，发现MPCS 比纯流体对流换热更好；Zhang 等[23]考虑了过冷度，发现对流换热强度随相变温度范围增加而下降；Roy 等[24]的研究表明特定条件下Stefan 数是唯一决定对流换热强度的因素；Hu 等[25]的结论与Charunyakorn 等[22]的类似，并进一步探讨了过冷度、相变温度和微胶囊尺寸对对流换热的影响；Sabbah等[26]利用二维对称模型的模拟结果表明，MPCS传热系数与融化界面位置相关。

上述模拟研究均基于等效均质模型，即忽略固相颗粒的影响，将两相流等效为均相流体，其物性参数如密度、热导率、黏度等均采用等效公式计算，未考虑颗粒分布的影响[27]，也未考虑微胶囊与基液间的相互作用，而无疑微胶囊颗粒的增多会增大颗粒之间的碰撞和微对流效应[28]，而颗粒浓度和颗粒行为对MPCS 传热特性又有重要影响。显而易见的，高颗粒浓度流体黏滞阻力大、颗粒碰撞频繁，牛顿流体假设不再成立，对流传热机理也更加复杂。有研究认为层流下相变微胶囊悬浮液对流传热系数随颗粒浓度升高而单调增加[29]，但刘丽等[30]制备了以丙醇/水溶液作为载流体，芯材为正二十六烷的高颗粒浓度相变微胶囊悬浮液，发现热导率随颗粒浓度先增加后减小。

图1 相变微胶囊悬浮液(MPCS)示意图Fig.1 Schematic diagram of MPCS

尽管目前对于MPCS 对流换热特性已有一定研究，但对于微胶囊颗粒行为对流动和换热的影响规律和细致机理仍不明晰。本研究因此发展了MPCS的离散相模型[31]更真实地模拟两相流动，深入认识颗粒行为对于传热的影响规律。本研究采用离散相模型对微胶囊颗粒质量分数、粒径大小和相变潜热对壁温控制和传热效果的影响，微胶囊颗粒分布对温度分布的影响展开了模拟研究，对各参数影响大小进行了比较分析，同时针对不同质量分数下离散相模型与传统单相流模型的模拟结果进行了比较分析。

1 数学模型

1.1 基本假设

研究中模拟区域为长1000 mm、半径5 mm 的恒热流水平圆管。采用较为简单、网格数较少的二维模型。为简化计算，进行了下列假设：（1）忽略重力以便采用对称轴模型；（2）不可压稳态层流，0＜Re＜2000；（3）除比热容外，流体物性参数不受温度影响；（4）微胶囊颗粒扩散运动可以忽略；（5）微胶囊颗粒为均质相变材料，忽略轴向换热和黏度耗散；（6）综合考虑颗粒与颗粒、颗粒与流体、颗粒与壁面相互作用产生的微对流效应，通过有效热导率考虑微对流效应的影响。

1.2 控制方程

图2 展示了二维模型的数学计算域，相变微胶囊流体从管道入口以一定的温度和速度流入，在沿管流动过程中，温度逐渐升高直至达到相变材料的融化温度，胶囊内的相变材料开始融化，受之影响微胶囊流体和壁面的温升会放缓。

对连续相，计算控制方程如下：连续方程

动量方程

能量方程

悬浮液的总焓由表观比热容积分得到

对离散相，其运动方程参见式（5）

式（5）为单位质量颗粒的运动方程，其中右边第一项为曳力项，第二项为重力项，第三项为附加的加速项（单位颗粒质量所受的力），主要包括压力梯度力和可以忽略的Saffman升力，自旋运动产生的Magnus 力，温度梯度产生的热泳力等；本研究中离散相与连续相密度几乎相同，故颗粒动量方程中第二项重力项较小可以忽略。式（6）中Re 为相对Reynolds 数，参见式（7）。CD为曳力系数，其表达式参见式（8），其中a1、a2、a3均为常数，取自Morsi 等[32]的研究结果。

微胶囊颗粒与流体间传热表达式参见式（9），忽略颗粒内部热阻和辐射换热[33]

在本研究中，对流传热系数由式（10）得到

1.3 边界条件

为便于说明，进行无量纲化，见式（11）～式（16）。

在圆管入口处需满足入口温度和入口流速恒定

在管道出口处，流场和温度场充分发展即轴向速度不再沿轴变化

在管道中心处，径向速度、热流和切向剪应力均为零

管道内壁满足无滑移和恒热流条件

1.4 相变微胶囊颗粒

相变微胶囊颗粒由内部相变材料和囊壁材料两部分组成，视为均质颗粒时其整体参数主要由核心材料和囊壁材料所占的体积分数决定，其中相变材料为石蜡，囊壁材料为聚甲基丙烯酸甲酯（polymethyl methacrylate, PMMA），石蜡占整体的体积分数取60%，石蜡和PMMA 的物理性质如表1所示。

表1 PMMA和石蜡的物理性质Table 1 Physical properties of PMMA and paraffin

微胶囊颗粒的整体物性参数可由下列关于相变材料体积分数的关系式得到。微胶囊颗粒密度参见式（21）[26]

相变材料质量分数可由式（22）得到

相变微胶囊颗粒比热容由式（23）得到[26]

相变微胶囊颗粒热导率由式（24）得到[34]

石蜡的相变潜热为190 kJ/kg，融化温度范围为

307～310 K。

本研究中在进行数值模拟时，相变材料的融化效应均采用等效比热容法来进行替代，等效比热容法即将比热容视为关于温度的单值函数，见图3。选取正弦函数来拟合石蜡的相变融化曲线，具体表达式见式（25），通过添加UDF来进行实现[35]。

1.5 初始条件

Sabbah 等[26]的研究表明在恒热流条件下，相比于纯流体，微胶囊流体需要更长的热入口段使流速充分发展，且管道中流体的温度范围与微胶囊颗粒的相变范围重合度越大，融化吸热效应对传热改善越明显，故二维模拟时微胶囊流体入口流速设定为0.0904 m/s（Re=900），以使微胶囊流体在有限长度的管道中尽快达到流速充分发展，微胶囊流体的入口温度为306 K，圆管恒热流为1273 W/m2，以使管内温度范围与相变温度范围尽可能重合。

图3 等效比热容法示意图Fig.3 Schematic diagram of equivalent specific heat method

1.6 计算方法

采用FLUENT中的有限体积法求解流动和传热控制方程。微胶囊颗粒的性质是利用C语言编写的UDF 并编译。计算域网格为矩形网格，综合考虑计算精度和计算时长，网格总数选取为12500个，其中轴向500 个，径向25 个。无量纲控制方程中连续性方程的收敛标准选取10-6，动量和能量方程的收敛标准选取10-7。

2 结果与讨论

2.1 模型验证

为验证采用模型的可靠性，本研究将采用离散相模型的模拟结果与已有圆管内MPCS 流动换热实验研究[36]报道的结果进行比较，计算区域是一个总长1300 mm，内径3.4 mm，外径4.0 mm 的圆管，可分为上游段、测试段、下游段三个部分，其中上游段长700 mm，测试段长400 mm，下游段长200 mm，上游段和下游段均保持绝热边界条件，以保证微胶囊流体达到流速充分发展，测试段为恒热流边界条件，热流为7957.7 W/m2，工况为Re=197，粒径10 μm，颗粒质量分数为2%时的无量纲壁面温度沿管径方向分布的比较结果如图4所示。可以看出模拟结果与实验结果是几乎一致的，由于实验测量时得到的壁温是若干离散点，模拟时得到的是壁温连续分布曲线，相比于实验结果中壁温曲线单调变化，模拟结果在近出口端壁温存在极大值点，这是由于壁面轴向导热引起的[37]。虽然在下面的讨论中本文所采用的管与本节所验证的管道长径比是不同的，但不同管道中MPCS 的流动和传热机理是相似的。因此基于上述数学模型模拟结果与实验结果良好的匹配性，可以认为该模型可适用于MPCS 流动换热特性研究。

图4 MPCS轴向温度分布模拟与实验结果比较Fig.4 Comparisons between numerical results and experimental results:temperature of MPCS along pipe

2.2 颗粒质量分数

为研究微胶囊颗粒质量分数对流动换热的影响，分别选取微胶囊颗粒质量分数为2%、5%和8%，颗粒粒径取10 μm，Re 取900，即入口流速为0.0904 m/s，其中质量分数为5%时的全管段温度场分布，近入口段速度场分布和颗粒浓度分布如图5 所示，由于入口温度接近相变温度的下限，且壁面热通量足够高，进入管路后MPCM 的温度迅速上升到相变阶段，MPCM 的温度在很长一段距离内保持不变，这是由于在这部分管路中能量以潜热的形式被吸收，应该指出的是近壁区温度升高超过相变上限温度，这意味着壁面附近的相变材料已经完全融化。由于本研究中管道长径比过大，全管道的速度场和颗粒浓度场分布展示效果较差，故选取近入口段速度场和颗粒浓度分布进行展示，可以看到随着MPCM的沿程流动，近壁面流速逐渐减小，产生了径向速度梯度，越靠近管道中心，颗粒流速越大，流速大压强小，颗粒会向速度较高区域聚集，符合真实物理规律。

不同质量分数下壁面温度分布对比和壁面Nu分布对比结果如图6 所示，定量分析结果见表2，其中ΔT 表示出口与入口的壁面温差，表中Nu 选取为出口处的局部值，下同。壁面Nu 计算方法见式（26）～式（28）[26]

其中平均温度计算方法如下

图5 Re=900、质量分数5%、颗粒粒径10 μm时的模拟结果Fig.5 Numerical results when Re=900,mass fraction=5%,particle size=10 μm

图6 不同颗粒质量分数下壁面温度和Nu分布Fig.6 Wall temperature and Nusselt number distribution under different particle mass fraction

表2 不同微胶囊颗粒质量分数对比Table 2 Comparison of microcapsules with different mass fraction

从上面的结果中可以看出随着微胶囊颗粒质量分数的增加，壁面温度上升越缓慢。这是由于相变微胶囊颗粒质量分数越高，相变时的融化吸热效应越明显；壁面Nu 随颗粒浓度增加而增加，这是由于在恒热流边界条件下，壁温与壁面附近流体之间的温差与Nu 呈反比，即颗粒质量分数越大，壁温与壁面附近流体之间温差越小，微胶囊流体对壁面的冷却效果更充分。本研究的颗粒质量分数影响与Wang 等[29]的结论一致，但在颗粒质量分数过高时微胶囊颗粒会发生团聚沉降等行为破坏MPCS 的稳定性，造成传热效果恶化。故在实际应用需求中，MPCS 颗粒浓度不可过高，应平衡传热效果和稳定性选取颗粒质量分数最优值。

2.3 颗粒粒径大小

为研究微胶囊颗粒粒径对流动换热的影响，分别选取微胶囊颗粒粒径大小为10、50 和100 μm，颗粒质量分数取5%，Re 取900，不同粒径大小壁面温度分布对比和壁面Nu 分布对比结果如图7 所示，定量分析结果见表3。

表3 不同微胶囊颗粒粒径大小对比Table 3 Comparison of microcapsules with different particle size

图7 不同颗粒粒径大小下壁面温度和Nu分布Fig.7 Wall temperature and Nu distribution under different particle size

从上面的结果中可以看出，在测试段管长10%处前，相变微胶囊颗粒未开始融化时，颗粒粒径大小对壁面温度和传热的影响几乎可以忽略，当相变颗粒开始融化时，颗粒粒径大小影响变得显著，整体趋势是颗粒粒径越小，壁面温度增长越缓慢，即壁面温度控制效果越好，这是因为在相同质量流量的情况下，颗粒粒径越小，颗粒比表面积越大，颗粒与流体之间的接触面积增大，换热效果也会得到改善，当颗粒发生融化吸热时效果更显著，壁面温度增长因此放缓；粒径越小，壁面Nu越大，这表明粒径越小，壁温与壁面附近流体的温度差距越小，微胶囊流体对壁面的冷却更好。从定量分析上来看，颗粒尺寸和颗粒质量分数的影响相当。颗粒尺寸影响相应结论与已有研究结论[38]相符合。

2.4 相变潜热

为研究相变潜热对流动换热的影响，分别模拟不发生相变的微胶囊、常规相变的微胶囊和相变潜热为正常值两倍时的相变微胶囊传热情况，常规相变潜热值为190 kJ/kg，微胶囊流体入口Re取900，颗粒质量分数取5%，粒径大小取10 μm，壁面温度分布对比和壁面Nu 分布对比结果如图8 所示，定量分析结果见表4。

表4 不同相变潜热对比Table 4 Comparison of microcapsules with different latent heat

从上面的结果可以看出，不发生相变时潜热视为零，则相变潜热越大，壁面温度增长越缓慢，壁面传热效果越好，在管长达到10%后，壁温曲线增长斜率放缓，表明相变微胶囊颗粒主要在10%管长后开始融化吸热。从不发生相变到常规相变的改善优于从常规相变到相变潜热为正常值两倍时的改善，从定量分析可以看出，相变潜热大小对壁温控制和壁面传热的影响大于颗粒质量分数和颗粒尺寸的影响。

图8 不同相变潜热下壁面温度和Nu分布Fig.8 Wall temperature and Nu distribution under different latent heat

2.5 模型比较

单相流模型与离散相模型最主要的差别，在于单相流模型将相变微胶囊颗粒与基液等效为新的均质流体，只需设定整体参数，且模拟过程中忽略了颗粒与流体间的传热，假设颗粒在管内均匀分布。将微胶囊颗粒与基液等效为新的均质流体时，物性参数主要与微胶囊颗粒的体积分数或质量分数相关，具体表达式如下。

微胶囊流体的密度表达式与式（21）形式相同，只是其中的体积分数含义变为相变微胶囊占整个流体区域的体积分数。

微胶囊流体的比热容计算如式（29）所示[26]

传统的关于微胶囊流体传热的数值模拟研究选取的几乎都是单相流模型，本研究则是选取了离散相模型。为明确不同模型的优缺点和适应场景，改变颗粒质量分数分别采用离散相模型和单相流模型进行模拟，对模拟结果进行比较分析，颗粒粒径取10 μm，Re取900，结果如图9所示，定量分析如表5所示。

从上面的结果中可以看出在其他条件相同情况下，微胶囊颗粒质量分数越大，两种模型模拟的出口和入口的壁面温度差距越大，但壁面Nu无明显差异。为进一步分析造成不同模型差异的主要原因，本研究从模型根本差异出发，通过比较不同颗粒分布对应的温度场分布来分析颗粒分布对温度场分布的影响，结果如图10所示。

表5 单相流模型与离散相模型对比Table 5 Comparison between discrete phase model and single phase model

从模拟结果可以看出，颗粒在管内径向分布并非为均匀分布，颗粒会在圆管中心处发生集聚，且颗粒质量分数越大颗粒集聚程度越高，中间流面径向温度分布越均匀。该结果表明传统单相流模型假设颗粒在管内均匀分布并非真实情况，且颗粒分布会对中间流面温度造成影响，表明单相流模型忽略颗粒与流体间的相互作用也会造成额外误差，进一步降低计算的准确性。因此对于本研究的颗粒质量分数变化范围，单相流模型模拟结果偏差随着颗粒质量分数增加而增大，进而造成与离散相模型模拟结果差异增大。

图9 不同质量分数下两种模型对比Fig.9 Comparison of two models under different mass fractions

图10 不同质量分数下中间流面颗粒浓度和温度分布Fig.10 Concentration and temperature distribution of middle stream surface under different mass fraction

3 结 论

本研究采用离散相模型研究了恒热流圆管中相变微胶囊流体的传热特性，并对微胶囊颗粒质量分数、颗粒尺寸、相变潜热大小对MPCS 换热特性的影响进行了定量分析，同时将不同颗粒质量分数下离散相模型与单相流模型的模拟结果进行比较分析，证实了离散相模型的优越性。本研究主要结论如下。

（1）随着微胶囊颗粒质量分数的增加，壁面温度上升变慢，壁面Nu逐渐增加。

（2）颗粒粒径大小越小，壁面温度增长越慢，壁面Nu越大。

（3）相变潜热越大，壁面温度增长越缓慢，壁面传热效果越好，相变潜热大小对壁温控制和壁面传热的影响大于颗粒质量分数和颗粒尺寸的影响。

（4）颗粒质量分数越大，单相流模型与离散相模型壁温模拟结果差别越大，壁面传热效果则与颗粒质量分数无明显联系。颗粒在管内径向分布并非为均匀分布，颗粒会在圆管中心处发生集聚，且颗粒质量分数越大颗粒集聚程度越高，出流面径向温度分布越均匀，造成单相流模型计算偏差越大。

符 号 说 明

A——换热面积，m2

a1,a2,a3——曳力系数拟合常数

CD——曳力系数

ceff——等效比热容，J/(kg·K)

cp——表观比热容，J/(kg·K)

d——圆管直径，m

dp——颗粒直径，m

E——单位质量流体总能量，J/kg

F——流体所受外力，N

FD——曳力，N

g——重力加速度，m2/s

H——悬浮液总焓，kJ/kg

hj——组分j的焓值，kJ/kg

href——参考温度下焓值，kJ/kg

hsf——固液相变潜热，kJ/kg

kb——整体热导率，W/(m·K)

keff——有效热导率，W/(m·K)

p——压力，Pa

p*——无量纲压力

q″——热通量，W/m2

r——径向位置，m

r*——无量纲径向位置

T——温度，K

Tref——参考温度，K

ΔTmelt——融化温度范围

t——时间，s

u——轴向速度，m/s

uin——圆管入口速度，m/s

up——颗粒轴向速度，m/s

ur——径向速度，m/s

ur*——无量纲径向速度

v——速度，m/s

x——轴向位置，m

x*——无量纲轴向位置

θ——无量纲温度

μ——黏度，Pa·s

μeff——等效黏度，Pa·s

ρ——密度，kg/m3

ρeff——等效密度，kg/m3

ρp——颗粒密度，kg/m3

χ——相变微胶囊体积分数,%

ω——相变微胶囊质量分数,%

下角标

b——体相

eff——等效

f——流体

in——入口

p——颗粒

r——径向

x——轴向