液冷系统中均热板气液相变的热质传输模拟

2022-06-07李维平李隆键崔文智郭梦婷

哈尔滨工业大学学报 2022年7期

李维平，李隆键，崔文智，郭梦婷

(1.重庆大学能源与动力工程学院，重庆400030；2.低品位能源利用技术与系统教育部重点实验室(重庆大学)，重庆400030)

液体冷却方式是电池热管理的主要散热方式，具有换热能力高等优点[1]，可以有效降低电池表面的最高温度，使之工作在合适的温度范围内，并同时保持较好的温度均匀性，提高电池在使用过程中的安全性，延长使用寿命[2]。然而，冷却液在液冷换热系统中会产生较大的温差，需要提高其换热表面的均温性。常见的液冷流道结构采用蛇形或盘管形结构[3-5]，该设计较为复杂，会导致流阻较高、冷却液循环泵功率较大等。为了改善液冷系统换热表面均温性的要求，可以将液冷系统和高效换热器进行结合。

均热板是一种利用相变工质在封闭腔室中相变传热的高效换热元件，具有较高导热能力和装配紧凑的特点[6]，常应用于电子电力设备的散热系统中，如芯片、LED生热元件[7]，可以在较小热源面积下改善生热表面的换热性能[8]，强化竖直方向上相变工质的传热与传质能力[9-11]。在实际电池热管理系统应用中，电池表面的换热面积较大[12]，常利用平板热管、热管等高导热换热元件[13-16]，将产生的热量及时地导入散热系统中，但未着重研究均热板内部的气液传热传质过程与其表面均温性之间的联系。对于液冷散热系统而言，冷却剂温度会在沿程换热过程中不断提高，这会对换热表面的均温性造成影响，影响元器件的寿命和工作效率。因此，需要研究一种均热板与液冷板的复合冷却方式，来改善换热表面的均温性能，探究均热板内部相变工质在水平方向上传热传质对于其换热表面均温性的影响关系。

本文研究了一种铝槽式均热板和直流式液冷板组合而成的复合液冷系统，建立了耦合有直流式液冷板对流换热与均热板相变传热的三维数值模型，采用Volume-of-fluid(VOF)多相流模型和蒸发冷凝模型，模拟了相变工质在均热板蒸汽腔中输运与传热的过程，并和实验结果进行了对比和验证，揭示了均热板内部相变工质在水平方向上的传热传质过程可以使加热表面获得较好的温度均匀性，为均热板在电池热管理液冷系统中的应用提供了支撑和建议。

1 实验对象及方法

本文研究了一种用于电子电力设备热管理的复合液冷系统，由铝槽式均热板(170 mm×75 mm×10 mm)和直流式液冷板组合而成。图1为均热板的俯视结构示意图，腔室内部尺寸为(140 mm×45 mm×4 mm)，由矩形支撑翅柱(30 mm×1 mm×4 mm)将腔室均匀分隔成横向和纵向通道，其中，纵向通道(140 mm×4 mm×4 mm)的方向与直流式液冷流道平行，用于引导相变工质在相变过程中沿长度方向流动，并起支撑作用。均热板底面布置有加热表面(150 mm×50 mm)，施加有恒定的热流。均热板顶部布置有直流式液冷板，如图2所示，液冷板内部由12个180 mm×2 mm×2 mm的直流式单流道组成，并通过肋片均匀分隔，冷却水沿X方向流动，在对流换热的影响下冷却水温度沿程增加。均热板的加热表面与冷却表面面积之比近似为1。另外，在均热板的上下表面，分别开有测温沟槽，沿冷却水流动方向依次布置有5个T型热电偶(Tc1/Th5～Tc5/Th5)，用以研究均热板在长度方向上的温度分布与温差。

图1 均热板内部结构的俯视示意图(mm)

图2 液冷板内部结构的俯视示意图(mm)

图3为复合液冷系统的实验装置与系统示意图，主要分为热源加热系统、冷却散热系统和测试传感系统。均热板布置在陶瓷加热片和液冷板之间，为了减小接触热阻，其接触表面铺有一层厚度为0.3 mm的导热垫片。直流电源为陶瓷加热片提供恒定40 W(±6.04 %)的输入功率。恒温水箱提供进口温度为25 ℃的冷却水，并在浮子流量计(±2.5 %)和阀门的控制下，将冷却水进口流量控制在0.05 L/min。另外，在冷却水的进出口处分别布置有标定后的T型热电偶(±0.2 ℃)，通过数据采集仪监测冷却水的进出口温度，以及均热板表面的温度分布。

图3 实验装置及系统示意图

实验前，先对均热板内部进行真空排气，然后注入体积分数为50 %的丙酮进行测试。实验采用丙酮作为均热板内部的相变工质，其相变温度在动力电池的工作温度范围内，与铝金属不发生化学反应，是一种合适的相变工质。实验时设定好直流电源的电流和电压，以及冷却水从恒温水箱输入液冷板的进口水温，并调节阀门将冷却水进口流量设定到实验工况，在该实验工况稳定后，测试均热板表面的温度分布和冷却水的进出口温度。冷却水温度沿流动方向依次增高，根据加热表面沿流道方向上布置测温点，定义加热表面的最大温差为

ΔTh=Max(Th1-5)-Min(Th1-5)

(1)

2 物理模型及数值模拟方法

2.1 几何模型

根据实验系统中槽式均热板和直流式液冷板布置结构的周期性和对称性特点，如图4所示，选取均热板内部腔室的纵向通道(140 mm×4 mm×4 mm)和4个矩形支撑翅柱(30 mm×1 mm×4 mm)在宽度上的一半，以及所对应的均热板盖板和基底，作为均热板的模拟对象(170 mm×2.5 mm×10 mm)，对应投影面积相同的液冷板区域(180 mm×2.5 mm×6 mm)作为散热单元，宽度方向上的横切面均处理为对称面。热量从均热板底部的加热表面传递给腔室内的相变工质，相变工质受热蒸发，在温度较低的冷却侧冷凝，再通过导热的形式将热量传递给上方的液冷板，与直流式流道中的冷却水换热，最后热量随冷却水沿X方向流动带走。该简化的三维计算物理模型，综合考虑了矩形支撑翅柱对于腔室内气液输运和传热的影响，以及实际冷却流道肋片的导热影响，研究了均热板加热表面沿长度方向(X方向)上的温度分布和均温性。

图4 计算物理模型示意图

2.2 多相流动与相变传热模型

气液两相流动过程的模拟，常采用3种多相流模型：VOF模型、混合模型(Mixture model)、欧拉模型(Eulerian model)。其中VOF模型可利用界面重构方式(Re-construction)观察腔室内部的气液两相流动，即气泡的生成与脱离，实时追踪相变过程，但其求解速度较慢，计算时间长；Mixture模型将气液多相流动过程简化为一种混合相流动形式，考虑了气液两相流动间的相互渗透和扩散，在处理可压缩气液多相流动的复杂问题上具有一定的优势；Eulerian模型则计算精度更高，但需要求解多个耦合方程式，这也使得该方法具有计算成本较高、计算稳定较差的不足。本文为了分析均热板内部气液工质在相变过程中清晰的流动过程，以及获得均热板加热表面的温度分布，因而选用VOF模型对均热板内部的气液流动过程进行数值模拟，并与实验结果进行了对比和验证，讨论了气液相变的传热与传质对于均热板加热表面均温性的影响。

在VOF模型中，均热板内部气液工质的流动过程，可以通过控制单元内各自i相的体积分数来表示，如气相体积分数αvap、液相体积分数αliq，其取值范围为0～1，在同一单元内气液相体积分数之和为1，并满足连续性方程：

(2)

式中：ρi和ui为i相的密度和速度，Sm代表相变过程中气液相之间传递的质量源项。

气液工质流动过程的速度场分布，可以通过动量方程求得：

(3)

ρ=αvapρvap+αliqρliq

(4)

μ=αvapμvap+αliqμliq

(5)

能量方程可以表示为

(6)

式中:k和E分别为气液工质的平均导热系数和能量，p和T分别为控制单元内气液相的平均压力和温度，能量源项SE为气液工质在流动与相变过程中的能量传递。

基于Lee[17]和Schepper等[18]的研究，均热板内部工质相变过程的模型采用了蒸发-冷凝模型(Evaporation-condensation model)，在能量方程和质量方程中分别引入质量源项和能量源项，可以实现蒸发冷凝过程中的质量传递和能量传递[19]。具体的质量源项Sm和能量源项SE如表1所示：

表1 蒸发-冷凝模型中质量和能量源项

表1中，T为气液相平均温度；Tsat为丙酮的饱和温度，当平均温度高于饱和温度，丙酮液体开始蒸发，当平均温度小于饱和温度，丙酮蒸汽开始冷凝；βe、βc分别为蒸发松弛时间因子和冷凝松弛时间因子，单位为s-1。有研究表明：蒸发/冷凝松弛时间因子之比设定为相变工质气/液相密度之比[20]，对于计算的收敛和稳定性有利，根据丙酮的物性参数，本文的蒸发和冷凝松弛时间因子分别设定为0.1和105。

封闭空间内的丙酮蒸汽可视为可压缩的理想气体。为了更好地模拟相变工质在水平方向上、受温差影响下的传热传质过程，本文根据REFPROP软件中丙酮的物性参数，将丙酮的饱和温度设定成与压力p有关的多项式：

Tsat=276.656 7+8.048 5×10-4×p-2.869 1×10-9×p2

(7)

2.3 边界条件及求解方法

模型的边界条件是根据实验工况而设定的，热量以定热流密度5 333 W/m2的形式施加于均热板底部的加热表面；液冷板进口设定为速度进口边界，冷却液入口流速设定为0.017 5 m/s，出口设置为压力出口边界，内部的流固耦合面为无滑移边界。模拟考虑了表面张力和重力对相变工质的影响，动量方程和能量方程均采用二阶迎风离散格式。VOF模型采用瞬态显式方案，并用PISO算法进行计算，根据界面重构(Re-construction)方式求解气液组分。计算初始，在蒸汽腔中Patch体积分数为50 %的丙酮液体，赋予均热板腔室内部温度为307.32 K，对应饱和压力为45 000 Pa。计算步长采用变时间步长方式，初始时间步长设定为0.000 1 s，当计算物理时间达到60 s，质量、动量残差<10-4，能量残差<10-6时，认为计算结果收敛，达到稳定状态。

2.4 网格无关性验证

为避免网格密度对计算结果的影响，需对物理模型的网格划分进行无关性验证，以此保证数值模拟结果的准确性。然而，三维模型的网格数过多则会大大增加VOF瞬态计算时长，本文综合考虑了以上因素，在计算成本允许的情况下，划分了3种不同网格尺寸的网格，如表2所示，采用网格收敛指数(Grid convergence index)的方式对网格的收敛和可靠性进行验证[21]。

表2 网格无关性验证

网格收敛指数GCI是一种对于网格质量标准化的分析方式，可以用来找到模拟结果对网格分辨率不敏感时最佳的网格数，其公式定义为

GCIfine=Fs|ε|/(rz-1)

(8)

式中：Fs为安全系数，对于3个样本量，其值取1.25；r为网格细化比，本文取值为2；ε为模拟结果f1、f2和f3之间的相对误差，本文选用T4点的温度作为对比参数，其中f1、f2的相对误差可以表示为

ε=(f2-f1)/f1

(9)

z为收敛阶数，可以定义为

z=ln((f3-f2)/(f2-f1))/ln(r)

(10)

计算得到每两种网格数之间的网格收敛系数，即GCI12and GCI23，理想的网格数应当满足以下关系式：

GCI23=rzGCI12

(11)

计算发现，当模型的网格数为238 624时，满足以上条件，且T4点的温度的温度波动范围<1%，可以认为，该模型尺寸满足网格无关性要求。

3 计算结果与分析

本文研究均热板腔室内部的传热传质过程及其与液冷板的耦合传热过程，使用VOF模型对相变工质的多相流动进行瞬态模拟，并与实验结果进行对比，获得均热板对加热表面的温度分布，分析了相变工质在均热板内部的水平流动对加热表面均温性的影响机理。

3.1 实验验证与传热结果

图5为均热板加热表面的温度分布，从图5中可以发现，模拟结果和实验结果趋势一致，高温区域集中在均热板右端。实验的最高温度为312.57 K，对应的VOF模型下模拟值的最高温度为311.58 K，这是由于模型与实际过程存在一定差异，实验中也存在一定的接触热阻，导致计算得到的加热表面温度较实际温度偏低，但该误差均<1%，可认为计算结果可靠。

图5 加热表面温度分布与实验结果对比

表3为模拟得到的加热表面温差与实验结果的对比，可以发现，VOF模型计算得到的加热表面温差与实验结果相近，表明了模拟结果可以较好地反映出实验过程中均热板内部相变过程的传热情况。另外，该工况下直流式液冷板冷却液的进出口温差可以达到11.45 K，实验和模拟结果中加热表面较小的温差也可体现出均热板可以使得加热表面保持较好的均温性。

表3 加热表面温差与实验结果对比

图6为不同加热功率下均热板加热表面温度分布的实验结果，可以发现，随着加热功率的增加，均热板加热表面的温度也在相应提高，但仍能保持较好的均温性。当加热功率为20 W时，加热表面温度分布均匀，最低与最高温度分别为304.5、305.4 K;而当加热功率增加到60 W时，加热表面的最低与最高温度分别增加到315.8、319.7 K。与此同时，冷却水进出口的最大温差也会随着加热功率的提高快速增加。从表4中加热表面的温差结果可以发现，当加热功率增加到60 W时，冷却水进出口温差已经达到了17.2 K，此时加热表面温差增加幅度较小，仍能很好地控制在3.9 K以内，说明采用均热板的复合液冷系统可以保证加热表面具有较好的均温性。

图6 不同加热功率下加热表面的温度分布

表4 不同加热功率下加热表面与冷却水温差

3.2 相变过程与流动特性分析

图7为均热板槽道中气液两相体积分数随时间

变化的瞬态计算结果，其中蓝色(深)表示液相，红色(浅)表示气相。随着计算开始，热量通过加热表面传入腔室，底面液体受热升温，当蒸发面的温度高于液体的饱和温度，腔室底部有小气泡生成，并逐渐变大、合并、脱离，最终在气液界面处破裂，进入蒸汽区。与此同时，产生的蒸汽上升至腔室上部，由于上侧的冷凝表面温度低于饱和温度，蒸汽不断冷凝成液体，并形成液滴，最终受重力的作用滴落返回液池。可以发现，随着换热过程的继续(10～60 s)，均热板的沸腾核心区域在不断缩小，最后集中在均热板的右端，如区域3。这表明在实际液冷系统中，均热板冷却表面散热温度不均会对内部的相变过程造成影响。由于液冷板冷却水沿程温升的影响，均热板右端靠近液冷板出口处区域温度较高，因此，该区域相变工质具有较大的过热度，沸腾强度明显强于前端。均热板左端的液体工质受入口处低温冷却水的影响，导致其工作温度未达到相变温度，热量主要以导热形式传递，从而造成该区域沸腾强度较弱。

均热板内部的相变过程在长度方向上的差异性，会影响气液相的传热与传质过程。图8为XY平面上气液两相的速度矢量图，表示了气液两相在长度方向上的流动轨迹，其中红色(上)代表气相，蓝色(下)代表液相，箭头的方向代表流动的方向，箭头的长度代表了速度的大小。由图8可知，均热板腔室上部的丙酮蒸汽沿着均热板长度方向从x=119 mm附近的沸腾核心区向两侧流动，其中区域1、2和3中的丙酮蒸汽主要从右向左流动，区域4中的丙酮蒸汽主要从左往右流动，蒸汽的流速随着流动而逐渐减小。均热板内部蒸汽定向的输运过程，会使得工质吸收的热量可以通过蒸汽的输运从热管右端传递到左端，并在腔室左端的低温区域进行冷凝，冷凝释放的热量可以抑制液冷板进口处的低温冷却水对加热表面均温性的影响。

图7 腔室内气液相的体积分数分布(VOF模型)

图8 XY平面上气液两相速度矢量图(VOF模型)

蒸汽在流动过程中的压强变化与温度变化息息相关。均热板在与直流式液冷板的换热过程中，由于冷却侧的换热温度不均，会导致均热板内部的蒸汽沿X方向存在一定的温差，根据Clausius-Clapeyron方程可知：

(12)

式中Δpvap和ΔTvap表示蒸汽流动中的压强变化和温差变化，Rg为单位质量的气体常数，Tvap和pvap为气体的温度和压力，hfg为相变潜热。因而，蒸汽在长度方向存在的温差，会促使蒸汽产生相应的压差，驱动相变产生的蒸汽从高温区域流向低温区域，及时地将热量从右端输运到左端，从而抑制了长度方向上较大的温差，改善了均热板加热表面的均温性。

图9为YZ平面上气液两相的速度矢量图，可以发现均热板内部的气液工质在不同宽度截面上的流动也存在一定差异性。由于液冷板进口处低温冷却水的影响，气液工质在均热板左端区域相变过程较弱，如截面x=35 mm处，气液工质在宽度方向上的流动较为平缓且界面清晰，呈自然对流流动方式。在均热板右端的沸腾核心区，如截面x=102 mm和x=135 mm处，工质的沸腾现象较为明显，槽道中心位置产生的蒸汽冲破气液交界面，往腔室顶部流动且流速较大。由此可见，气液工质在宽度方向上的流动和相变过程，也与均热板在长度方向上的温差有关。为了保证均热板在长度方向上的均温性，可利用矩形支撑柱对气液工质在宽度方向上的流动进行限制，引导气液工质在长度方向上做定向流动，改善均热板在长度方向上的均温性能。