引言CO2作为拥有优良环境性能和换热性能的自然工质，近年来在汽车空调、热泵、复叠式制冷等领域被广泛应用。与常规换热器相比，CO2微通道平行流蒸发器具有高效、耐压、制冷剂充注量小等优点。Kim[1]等研究CO2汽车空调蒸发器的换热性能，得到空气侧和制冷剂侧的换热系数计算关联式。石冬冬[2]等研究CO2微通道平行流蒸发器在不同结构参数和工况下的换热规律。为了解决空气流量或迎风面积较小时空气出口计算温度异常的问题，改进了空气侧的计算方法，并对运行参数状况和换热特性进行仿真分析。1建立蒸发器分布参数模型1.1CO2侧换热系数CO2在蒸发器中的两相沸腾过程经历环状流、间歇流、干涸区和雾状流等阶段。制冷剂侧换热系数计算采用Cheng[3]的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.67Wev0.17Frv0.348ρl/ρv0.25q/qcrit0.7 (2)式中：We——韦伯数；Fr——弗劳德数；下标v——气态相关参数；下标l——液态相关参数；qcrit——临界热流密度。qcrit=0.131ρv0.5γgρl-ρvδ0.25 (3)式中：γ——流体表面张力；g——重力加速度，取9.8 m/s2；δ——液膜厚度，m。（3）干涸前两相区CO2侧换热系数为：hnb=131P/Pcrit-0.006 3-lgP/Pcrit-0.55m-0.5q0.58 (4)hcb=0.01334Gr(1-x)δμl(1-ε)0.69Prl0.4klδ (5)S=1                                                                                                                                   ,   xxIA1-1.14D/Dref21-δ/δIA2.2,   x≥xIA (6)htp=Shnb3+hcb31/3 (7)式中：hnb——池沸腾换热系数，W/(m2·K)；Pcrit——CO2的临界压力，MPa；m——制冷剂的摩尔质量，g/mol；hcb——核态沸腾换热系数，W/(m2·K)；Gr——制冷剂的质流密度，kg/(m2·s)；ε——横截面蒸汽空隙率；r1——液相状态对应的特朗普数；kl——液相对应导热系数；Dref——0.007 53；δIA——间歇流与环状流分界点处的液膜厚度，m；htp——两相区总换热系数，W/(m2·K)；S——强化因子。（4）干涸后两相区CO2换热系数为：hmist=0.011 7ReH0.79Prv1.06Y-1.83kv/D (8)Y=1-0.1[(ρl/ρv-1)(1-x)]0.4 (9)hdryout=htpxDI-x-xDIxDE-xDI[htpxDI-hmistxDE] (10)式中：xDI——干涸区起始点；xDE——干涸区结束点，取0.999；k——导热系数，W/(m·K)。为了确保压缩机处于干工况，蒸发器出口需要具有一定的过热度。过热区CO2的湍流和层流换热分别采用Gnielinski公式和Sieder-Tate公式[4]计算。Re≥2 300时，采用Gnielinski公式计算CO2湍流区域努塞尔数：Nuf=fD/8Ref-1 000Prf1+12.7fD/80.5Prf2/3-1 (11)式中：fD——管内湍流流动的Darcy阻力系数，按照Filonenko公式计算。fD=1.82lgRe-1.64-2 (12)Re2 300时，采用Sieder-Tate公式计算：Nuf=1.86RefPrfl/D1/3ηfηw0.14 (13)式中：l——管长，m；D——当量直径，m；ηw——壁面温度对应的动力黏滞系数。1.2空气侧换热系数蒸发压力能够决定蒸发温度，从而影响壁温。壁面温度高于空气露点温度时，空气侧为干工况，反之为湿工况，需考虑传质与潜热。Chen[5]拟合百叶窗翅片干工况下的换热和压降关联式；Kim[6]等对干、湿工况下的换热因子(j)分别进行拟合。本模型采用Kim[6]等的关联式计算换热因子。空气侧百叶窗翅片结构如图1所示。10.3969/j.issn.1004-7948.2023.01.009.F001图1空气侧百叶窗翅片结构空气侧为干工况时：j=ReLp-0.457θ900.257FpLp-0.13FlLp-0.29 TdLp-0.235LlLp0.65TpLp-0.279δfLp-0.05 (14)空气侧为湿工况时：j=ReLp-0.512θ900.25FpLp-0.171FlLp-0.29 TdLp-0.248LlLp0.68TpLp-0.275δfLp-0.05 (15)求得j后，干、湿工况下的显热换热系数αa为：αa=Gacpaj/Pr2/3 (16)湿工况下需要考虑传质，由契尔顿-柯本尔热质交换公式可以得到湿工况下的传质系数：hmαa=le2/3/cpa (17)式中：αa——对流换热系数，W/(m2·s)；hm——传质系数，kg/(m·s)；le——刘易斯数，针对温度为-40 ℃到常温的空气，取1。1.3数学模型以单层扁管与其上下二分之一的百叶窗翅片作为一个换热单元，流量分配均匀的情况下，平行流蒸发器可以被看作由35个相同的换热单元组成，将换热单元沿制冷剂流动方向划分。单个控制单元结构如图2所示。10.3969/j.issn.1004-7948.2023.01.009.F002图2单个控制单元结构空气侧质量守恒方程为：∂ρa∂t+∂∂zρaua=sQcondVa (18)空气侧能量守恒方程为：∂∂tρaha+∂∂zρauaha=-Q1Va (19)空气侧冷凝水质量守恒：∂∂tρad+∂∂zρauad=-QcondVa (20)式中：s——开关量，不发生减湿时为0，发生减湿时为1；Va——微元空气侧体积，m3；Qcond——质交换量。考虑空气与管壁的热平衡与传质，空气侧总换热量和质交换量为：Qcond=hmAaηa(dam-dswm) (21)Q1=αaAaηa(Tam-Tpm)+rQcond (22)式中：下标swm——壁面温度对应的饱和参数；下标am——空气在该微元段内的平均参数；r——壁面温度下水的潜热，J/kg；Tpm——壁面温度，K。空气与管壁之间的热湿交换以焓差驱动表示为：Q1=hmAaηa(ham-hswm) (23)因为模型为稳态仿真，非稳态项为0。沿空气流动方向的自由流通面积不变，能量守恒方程和冷凝水质量守恒方程可简化为：GaAy(hao-hai)=-Q1 (24)GaAy(dao-dai)=-Qcond (25)式中：Ga——空气的质流密度，kg/(m2·s)；Ay——空气侧迎风面积，m2；下标ai——空气入口参数；下标ao——空气出口参数。制冷剂侧质量守恒方程为：∂∂tρr+∂∂yρrur=0 (26)制冷剂侧能量守恒方程为：∂∂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算法与实现在微通道平行流蒸发器中，制冷剂和空气不接触，通过壁面和翅片进行耦合传热。为了准确地模拟换热过程，必须确定壁温。文中采用二分法，壁温假设值的初始上限为空气入口温度，初始下限为制冷剂入口温度，壁温假设值为初始上限和初始下限的算数平均数。判断壁温假设值和空气露点温度的大小关系，分别代入干、湿工况进行运算。干工况下含湿量不发生变化，仅需使微元内空气内能变化和空气与壁温之间的换热达到平衡。沿空气流动方向（扁管宽度）尺寸较小，可认为微元内该方向壁温为一定值。空气与壁面之间为变温差传热，且空气侧的进出口温差很大。文献[2]采用算术平均温度得到的壁面与空气间换热量偏大。文中采用对数平均温差计算壁面与空气之间的换热。Qa=αaAηaηa(Tai-Tao)/lnTai-TpmTao-Tpm (30)湿工况下空气侧的含湿量发生变化，湿空气的状态由焓值和含湿量共同确定。因此，需要分别对焓值和含湿量进行迭代。湿空气的焓值与干球温度和含湿量有关，含湿量的变化会引起焓值变化。根据空气侧凝结水质量守恒方程对含湿量进行迭代，收敛后再判断是否满足热平衡关系。空气出口焓值假设值的上限为空气入口焓，下限为壁温对应的饱和湿空气焓，利用二分法对出口焓进行调整，直至同时满足热湿平衡。根据空气侧和制冷剂侧热平衡可以计算制冷剂侧的热流密度，按照Cheng关联式可得制冷剂侧换热系数，从而求得壁温，与假设壁温比较迭代，直至满足误差要求。壁面温度计算流程如图3所示。10.3969/j.issn.1004-7948.2023.01.009.F003图3壁面温度计算流程蒸发器的运行工况分为两种情况，分别为出口有过热度和出口无过热度，分别对应全程只有两相区和两相过热共存的工况。为了提高运算效率，对这两部分程序进行整合。初始默认蒸发器全程均为两相区，沿制冷剂流动方向将管长均分为N个微元。在已知制冷剂、空气入口参数的情况下，调用求壁温函数，可以得到各个微元的壁温和该微元出口的空气状态和制冷剂焓值。判断微元出口制冷剂是否处于过热状态，如果处于过热状态，说明蒸发器为两相过热共存的工况，跳出循环进入两相过热共存程序，反之继续计算下一微元，直至N个微元计算结束，第N个微元出口参数即为蒸发器出口状态。制冷剂由两相区转变为过热区时，传热急剧恶化，壁面温度和制冷剂温度均迅速上升，迭代容易不收敛。为了保证程序的稳定性，先确定两相区和过热区的分界点，再分别计算。假设两相区的换热长度，将该长度划分为N个微元，调用壁温函数计算该假设长度对应的出口焓是否等于两相区过热区的临界焓，如不满足误差要求，重新假设两相区长度迭代至收敛，从而得到两相区长度和两相区各个微元的进出口参数。两相区出口制冷剂状态即为过热区入口状态。平行流-百叶窗蒸发器为叉流换热器，每一个微元空气侧进口状态参数均相同。将过热区划分为M个微元，调用壁温函数进行计算，得到蒸发器沿制冷剂流动方向各个微元的换热状况。两相和过热共存时微元划分方式如图4所示。10.3969/j.issn.1004-7948.2023.01.009.F004图4两相和过热共存时微元划分方式2数值模拟结果及分析2.1平行流微通道蒸发器结构参数计算采用的微通道换热器为单流程，由35层平行微通道扁管组成，每一扁管内有18个当量直径为1.096 mm的微通道[2]。空气侧为百叶窗翅片。扁管及百叶窗翅片结构参数如表1所示。10.3969/j.issn.1004-7948.2023.01.009.T001表1扁管及百叶窗翅片结构参数参数数值扁管高度Dm/mm2微通道管道长度Ly/mm810扁管间距Tp/mm9.62扁管宽度Td/mm25百叶窗间断长度Lp/mm1百叶窗开缝长度Ll/mm6.35百叶窗翅片倾斜角度θ/(°)30翅片高度Fl/mm7.62翅片间距Fp/mm0.7翅片宽度Fd/mm25.4翅片厚度δf/mm0.1中间微通道宽度pw/mm0.94中间微通道高度ph/mm1.3两端微通道半径pr/mm0.65两端微通道宽度pwo/mm1.0052.2平行流微通道蒸发器模型验证为了减少试验中由重力作用造成的分液不均现象，选取第18排扁管为研究对象，共布置8个测点，测量迎面风速、空气温度湿度、壁面温度等参数。蒸发器迎风面积较大，表面各个位置的风速不均匀，取各个测点风速的平均值(1.8 m/s)为迎面风速。CO2质量流量15.67 g/s，入口干度0.276 7，蒸发压力3.2 MPa，空气入口温度23 ℃，入口相对湿度25%时，第18层扁管试验测量值如表2所示。10.3969/j.issn.1004-7948.2023.01.009.T002表2第18层扁管试验测量值项目测点910111213141516风速/(m/s)1.451.911.611.431.811.821.961.87出口温度/℃22.422.121.721.119.08.26.55.9出口相对湿度/%26.5326.8627.6928.9032.7266.1560.0761.12壁面温度/℃22.7722.1221.5320.5117.126.00-1.70-1.90将试验测得的CO2温度、壁面温度、空气进出口温度与数值模拟结果进行对比。温度参数数值解与试验数据对比如图5所示。10.3969/j.issn.1004-7948.2023.01.009.F005图5温度参数数值解与试验数据对比由图5可知，管长0.20 m处为两相区与过热区的分界点。在两相区，制冷剂侧对流换热系数很大，壁温接近蒸发温度；CO2进入过热区后传热急剧恶化，以显热的方式进行热交换，CO2温度与壁温迅速上升，与空气换热温差减小，换热量减小，出口处壁温与CO2温度均接近空气入口温度。因为微通道平行流蒸发器为叉流换热器，空气入口状态相同，两相区空气与壁面温差大，在两相区的换热量远大于过热区，两相区的长度占总管长的23.53%，但换热量占总换热量的83.3%。CO2侧对流换热系数数值解与试验数据对比如图6所示。沿管长的累计换热量与试验数据对比如图7所示。10.3969/j.issn.1004-7948.2023.01.009.F006图6CO2侧对流换热系数数值解与试验数据对比10.3969/j.issn.1004-7948.2023.01.009.F007图7沿管长的累计换热量与试验数据对比2.3两相区换热特性分析在两相区，CO2侧对流换热系数高，空气与制冷剂温差大，换热效率高。因此，工程实践中要求CO2平行流蒸发器的换热段以两相区为主，仅在出口处适当设置过热区。两相区的换热特性决定了蒸发器的换热性能。2.3.1蒸发压力对两相区对流换热系数的影响平行流微通道蒸发器的蒸发压力为3.0~4.5 MPa，CO2质流密度为400 kg/(m2·s)，入口干度为0.2，空气侧入口温度为25 ℃，相对湿度为30%，迎面风速为3 m/s时，CO2侧对流换热系数随蒸发压力的变化如图8所示。10.3969/j.issn.1004-7948.2023.01.009.F008图8CO2侧对流换热系数随蒸发压力的变化由图8可知，随着蒸发压力的增大，干涸点逐渐后移。因为蒸发压力越大，制冷剂蒸发温度越高，与空气温差小，换热量小，从而使CO2干度上升缓慢，不易发生干涸。蒸发压力较低时，蒸发温度低，换热更强，但发生干涸后会传热恶化。工程实践中，应综合考虑二者的影响，确定合适的蒸发压力。2.3.2CO2质流密度对两相区对流换热系数的影响CO2质流密度为250~400 kg/m2s，蒸发器的蒸发压力为3 MPa，入口干度为0.2，空气侧入口温度为30 ℃，相对湿度为30%，迎面风速为3 m/s时，CO2侧对流换热系数随CO2质流密度的变化如图9所示。10.3969/j.issn.1004-7948.2023.01.009.F009图9CO2侧对流换热系数随CO2质流密度的变化由图9可知，CO2质流密度较小时，微通道内易发生干涸，且干涸后换热状况急剧恶化。随着质流密度的增大，干涸点的位置逐渐后移，因为质流密度的提高会增加二次流和液滴沉积，对干涸起推迟效应。干涸发生前，随着质流密度的增大，CO2侧的换热系数不断增大，因为提高质流密度使微通道内湍流发展更充分，更有利于换热。因此，增大CO2质流密度是提高微通道蒸发器换热性能的有效途径。2.3.3CO2入口干度对两相区对流换热系数的影响微通道平行流蒸发器的CO2入口干度为0.20~0.35，蒸发压力为3 MPa，质流密度为400 kg/(m2·s)，空气侧入口温度为25 ℃，相对湿度为30%，迎面风速为3 m/s时，CO2侧对流换热系数随CO2入口干度的变化如图10所示。10.3969/j.issn.1004-7948.2023.01.009.F010图10CO2侧对流换热系数随CO2入口干度的变化由图10可知，在较低蒸汽干度下，扁管从入口到出口均为两相区，换热方式以核态沸腾为主，传热强度高。随着CO2入口干度的增大，管内干涸点的位置前移，处于干涸状态的区间增多，传热恶化。因此，在入口流量一定的情况下，高入口干度不利于换热。3结语文中提出一种适用于CO2的微通道平行流蒸发器的稳态仿真算法，在空气侧对出口焓和出口含湿量进行双重迭代，考虑换热过程中空气物性参数的变化，提高了对空气侧换热特性预测的准确性。将两相、过热区共存工况下的数值模拟结果与试验数据进行对比。在两相区，CO2侧换热系数和壁温的数值解与试验值基本一致，空气出口温度的数值模拟结果与试验值具有一定偏差，可能与空气侧沿管长方向风速不均有关，利用平均风速代替不同测点的实际风速时，换热量数值解与试验值吻合良好。文中还探究了两相区CO2侧换热系数随质流密度、蒸发压力和入口干度的变化规律，可用于指导生产实践和优化设计。