1、节能基础科学34节能 ENERGY CONSERVATIONNO.01 2023CO2微通道平行流蒸发器换热特性仿真分析梁方 董波*秦妍 陈聪 李维仲（大连理工大学能源与动力学院 海洋能源利用与节能教育部重点实验室，辽宁 大连 116024）摘要：提出一种适用于CO2微通道平行流蒸发器的稳态仿真模型，利用分布参数法分析了百叶窗翅片侧空气与扁管内CO2耦合传热并编制程序。模型考虑了空气侧的干、湿工况与蒸发器出口制冷剂是否有过热度的情况，模型数值解与试验结果吻合，可用于指导优化设计。关键词：二氧化碳；微通道蒸发器；换热特性；百叶窗翅片中图分类号：TK124 文献标识码：A 文章编号：1004-79
2、48（2023）01-0034-05 doi：10.3969/j.issn.1004-7948.2023.01.009引言CO2作为拥有优良环境性能和换热性能的自然工质，近年来在汽车空调、热泵、复叠式制冷等领域被广泛应用。与常规换热器相比，CO2微通道平行流蒸发器具有高效、耐压、制冷剂充注量小等优点。Kim1等研究CO2汽车空调蒸发器的换热性能，得到空气侧和制冷剂侧的换热系数计算关联式。石冬冬2等研究CO2微通道平行流蒸发器在不同结构参数和工况下的换热规律。为了解决空气流量或迎风面积较小时空气出口计算温度异常的问题，改进了空气侧的计算方法，并对运行参数状况和换热特性进行仿真分析。1建立蒸发器分
3、布参数模型1.1CO2侧换热系数CO2在蒸发器中的两相沸腾过程经历环状流、间歇流、干涸区和雾状流等阶段。制冷剂侧换热系数计算采用Cheng3的CO2流动沸腾传热模型，利用环状流和间歇流的分界点、干涸点将整个沸腾过程划分为三大区域。Cheng模型涵盖了CO2沸腾换热的所有流型，预测精度较高。（1）间歇流与环状流分界点干度为：xIA=1.81/0.875(v/l)-1/1.75(l/v)-1/7+1-1(1)式中：v气相密度，kg/m3；l液相密度，kg/m3；l气相动力黏度，N S/m2；v液相动力黏度，N S/m2。（2）干涸点为：xDI=0.58exp0.52-0.67We0.17vFr0.
4、348v()l/v0.25()q/qcrit0.7(2)式中：We韦伯数；Fr弗劳德数；下标v气态相关参数；下标l液态相关参数；qcrit临界热流密度。qcrit=0.1310.5vg(l-v)0.25(3)式中：流体表面张力；g重力加速度，取9.8 m/s2；液膜厚度，m。（3）干涸前两相区CO2侧换热系数为：hnb=131(P/Pcrit)-0.006 3-lg(P/Pcrit)-0.55m-0.5q0.58(4)hcb=0.0133 4Gr(1-x)l(1-)0.69Pr0.4lkl(5)S=1 ,xxIA1-1.14()D/Dref2()1-/IA2.2,xxIA(6)htp=(Shn
5、b)3+h3cb1/3(7)式中：hnb池沸腾换热系数，W/(m2 K)；PcritCO2的临界压力，MPa；m制冷剂的摩尔质量，g/mol；hcb核态沸腾换热系数，W/(m2 K)；Gr制冷剂的质流密度，kg/(m2 s)；横截面蒸汽空隙率；r1液相状态对应的特朗普数；kl液相对应导热系数；Dref0.007 53；IA间歇流与环状流分界点处的液膜厚度，m；htp两相区总换热系数，W/(m2 K)；S强化因子。（4）干涸后两相区CO2换热系数为：hmist=0.011 7Re0.79HPr1.06vY-1.83(kv/D)(8)Y=1-0.1(l/v-1)(1-x)0.4(9)hdryout
6、=htp(xDI)-x-xDIxDE-xDIhtp(xDI)-hmist(xDE)(10)式中：xDI干涸区起始点；xDE干涸区结束点，取0.999；k导热系数，W/(m K)。作者简介：梁方（1997），女，硕士，研究方向为CO2热泵系统。通信作者：董波（1982），女，博士，副教授，研究方向为CO2热泵系统。收稿日期：2022-06-08引用本文：梁方，董波，秦妍，等.CO2微通道平行流蒸发器换热特性仿真分析 J.节能，2023，42（1）：34-38.节能基础科学35NO.01 2023节能 ENERGY CONSERVATION为了确保压缩机处于干工况，蒸发器出口需要具有一定的过热度。
7、过热区CO2的湍流和层流换热分别采用Gnielinski公式和Sieder-Tate公式4计算。Re2 300时，采用Gnielinski公式计算CO2湍流区域努塞尔数：Nuf=()fD/8()Ref-1 000 Prf1+12.7()fD/80.5()Pr2/3f-1(11)式中：fD管内湍流流动的Darcy阻力系数，按照Filonenko公式计算。fD=(1.82lgRe-1.64)-2(12)Re2 300时，采用Sieder-Tate公式计算：Nuf=1.86(RefPrfl/D)1/3(fw)0.14(13)式中：l管长，m；D当量直径，m；w壁面温度对应的动力黏滞系数。1.2空气侧
8、换热系数蒸发压力能够决定蒸发温度，从而影响壁温。壁面温度高于空气露点温度时，空气侧为干工况，反之为湿工况，需考虑传质与潜热。Chen5拟合百叶窗翅片干工况下的换热和压降关联式；Kim6等对干、湿工况下的换热因子(j)分别进行拟合。本模型采用Kim6等的关联式计算换热因子。空气侧百叶窗翅片结构如图1所示。空气侧为干工况时：j=Re-0.457Lp(90)0.257(FpLp)-0.13(FlLp)-0.29(TdLp)-0.235(LlLp)0.65(TpLp)-0.279(fLp)-0.05(14)空气侧为湿工况时：j=Re-0.512Lp(90)0.25(FpLp)-0.171(FlLp)-
9、0.29(TdLp)-0.248(LlLp)0.68(TpLp)-0.275(fLp)-0.05(15)求得j后，干、湿工况下的显热换热系数a为：a=Gacpaj/Pr2/3(16)湿工况下需要考虑传质，由契尔顿-柯本尔热质交换公式可以得到湿工况下的传质系数：hma=le2/3cpa(17)式中：a对流换热系数，W/(m2 s)；hm传质系数，kg/(m s)；le刘易斯数，针对温度为-40 到常温的空气，取1。1.3数学模型以单层扁管与其上下二分之一的百叶窗翅片作为一个换热单元，流量分配均匀的情况下，平行流蒸发器可以被看作由35个相同的换热单元组成，将换热单元沿制冷剂流动方向划分。单个控制单
10、元结构如图2所示。图2单个控制单元结构空气侧质量守恒方程为：at+z(aua)=sQcondVa(18)空气侧能量守恒方程为：t(aha)+z(auaha)=-Q1Va(19)空气侧冷凝水质量守恒：t(ad)+z(auad)=-QcondVa(20)式中：s开关量，不发生减湿时为0，发生减湿时为1；Va微元空气侧体积，m3；Qcond质交换量。考虑空气与管壁的热平衡与传质，空气侧总换热量和质交换量为：Qcond=hmAaa(dam-dswm)(21)Q1=aAaa(Tam-Tpm)+rQcond(22)式中：下标swm壁面温度对应的饱和参数；下标am空气在该微元段内的平均参数；r壁面温度下水的
11、潜热，J/kg；Tpm壁面温度，K。空气与管壁之间的热湿交换以焓差驱动表示为：Q1=hmAaa(ham-hswm)(23)因为模型为稳态仿真，非稳态项为0。沿空气流动方向的自由流通面积不变，能量守恒方程和冷凝水质量守恒方程可简化为：GaAy(hao-hai)=-Q1(24)图1空气侧百叶窗翅片结构节能基础科学36节能 ENERGY CONSERVATIONNO.01 2023GaAy(dao-dai)=-Qcond(25)式中：Ga空气的质流密度，kg/(m2 s)；Ay空气侧迎风面积，m2；下标ai空气入口参数；下标ao空气出口参数。制冷剂侧质量守恒方程为：t(r)+y(rur)=0(26)
12、制冷剂侧能量守恒方程为：t(rhr)+y(rurhr)=QrVr(27)式中：下标r制冷剂相关参数；Qr制冷剂侧微元换热量，J；Vr制冷剂侧微元体积，m3。经化简得：GrArs(hro-hri)=-Qr(28)根据制冷剂与壁面的热平衡，可得：Qr=rAr(Trm-Tpm)(29)式中：Gr制冷剂质流密度，kg/(m2 s)；Ars制冷剂侧自由流通截面积，m2；Ar制冷剂侧微元换热面积，m2。1.4算法与实现在微通道平行流蒸发器中，制冷剂和空气不接触，通过壁面和翅片进行耦合传热。为了准确地模拟换热过程，必须确定壁温。文中采用二分法，壁温假设值的初始上限为空气入口温度，初始下限为制冷剂入口温度，壁
13、温假设值为初始上限和初始下限的算数平均数。判断壁温假设值和空气露点温度的大小关系，分别代入干、湿工况进行运算。干工况下含湿量不发生变化，仅需使微元内空气内能变化和空气与壁温之间的换热达到平衡。沿空气流动方向（扁管宽度）尺寸较小，可认为微元内该方向壁温为一定值。空气与壁面之间为变温差传热，且空气侧的进出口温差很大。文献2采用算术平均温度得到的壁面与空气间换热量偏大。文中采用对数平均温差计算壁面与空气之间的换热。Qa=aAaa(Tai-Tao)lnTai-TpmTao-Tpm(30)湿工况下空气侧的含湿量发生变化，湿空气的状态由焓值和含湿量共同确定。因此，需要分别对焓值和含湿量进行迭代。湿空气的焓
14、值与干球温度和含湿量有关，含湿量的变化会引起焓值变化。根据空气侧凝结水质量守恒方程对含湿量进行迭代，收敛后再判断是否满足热平衡关系。空气出口焓值假设值的上限为空气入口焓，下限为壁温对应的饱和湿空气焓，利用二分法对出口焓进行调整，直至同时满足热湿平衡。根据空气侧和制冷剂侧热平衡可以计算制冷剂侧的热流密度，按照Cheng关联式可得制冷剂侧换热系数，从而求得壁温，与假设壁温比较迭代，直至满足误差要求。壁面温度计算流程如图3所示。蒸发器的运行工况分为两种情况，分别为出口有过热度和出口无过热度，分别对应全程只有两相区和两相过热共存的工况。为了提高运算效率，对这两部分程序进行整合。初始默认蒸发器全程均为两
15、相区，沿制冷剂流动方向将管长均分为N个微元。在已知制冷剂、空气入口参数的情况下，调用求壁温函数，可以得到各个微元的壁温和该微元出口的空气状态和制冷剂焓值。判断微元出口制冷剂是否处于过热状态，如果处于过热状态，说明蒸发器为两相过热共存的工况，跳出循环进入两相过热共存程序，反之继续计算下一微元，直至N个微元计算结束，第N个微元出口参数即为蒸发器出口状态。制冷剂由两相区转变为过热区时，传热急剧恶化，壁面温度和制冷剂温度均迅速上升，迭代容易不收敛。为了保证程序的稳定性，先确定两相区和过热区的分界点，再分别计算。假设两相区的换热长度，将该长度划分为N个微元，调用壁温函数计算该假设长度对应的出口焓是否等于
16、两相区过热区的临界焓，如不满足误差要求，重新假设两相区长度迭代至收敛，从而得到两相区长度和两相区各个微元的进出口参数。两相区出口制冷剂状态即为过热区入口状态。平行流-百叶窗蒸发器为叉流换热器，每一个微元空气侧进口状态参数均相同。将过热区划分为M个微元，调用壁温函数进行计算，得到蒸发器沿制冷剂流动方向各个微元的换热状况。两相和过热共存时微元划分方式如图4所示。图3壁面温度计算流程节能基础科学37NO.01 2023节能 ENERGY CONSERVATION2数值模拟结果及分析2.1平行流微通道蒸发器结构参数计算采用的微通道换热器为单流程，由35层平行微通道扁管组成，每一扁管内有18个当量直径为1.096 mm的微通道2。空气侧为百叶窗翅片。扁管及百叶窗翅片结构参数如表1所示。2.2平行流微通道蒸发器模型验证为了减少试验中由重力作用造成的分液不均现象，选取第18排扁管为研究对象，共布置8个测点，测量迎面风速、空气温度湿度、壁面温度等参数。蒸发器迎风面积较大，表面各个位置的风速不均匀，取各个测点风速的平均值(1.8 m/s)为迎面风速。CO2质量流量15.67 g/s，入口干度0.276
