CN102779217A - 一种结霜工况下的制冷系统计算机仿真性能计算方法 - Google Patents

Abstract

本发明公开了一种结霜工况下的制冷系统计算机仿真性能计算方法，包括以下步骤：分别建立湿空气物性计算模块、冷媒物性计算模块、冷凝器结构计算模块、蒸发器结构计算模块和蒸发器结霜计算模块；建立制冷循环系统计算模块。本发明在模拟了蒸发器霜层变化的基础上，蒸发温度的变化会导致整个制冷循环工作状况的变化。本发明可以有效地指导设计人员避免由于结霜而产生的系统非正常工作地风险。同时也为制冷系统的除霜控制提供有效地理论数据，更有效的提高系统的制冷量及能效比。本发明首先将计算时间进行均匀步长分割，利用计算机辅助工具对程序内的结构与换热参数进行复杂的迭代计算。本发明实现了整个系统在结霜工况下的动态仿真。

Description

一种结霜工况下的制冷系统计算机仿真性能计算方法

技术领域

本发明涉及一种制冷与低温技术中的计算机辅助设计技术，特别是一种结霜工况下的制冷系统计算机仿真性能计算方法。

背景技术

利用计算机辅助技术进行制冷设备的仿真与模拟已经成为各个制冷产品开发的主要手段。随着计算机技术的蓬勃发展，计算机的计算能力也显著提高，越来越多的企业把关注度投入到仿真模拟中来，同时也利用计算机辅助技术使设计开发工作达到了一个前所未有的高度。

结霜现象普遍存在于制冷、低温、空气调节、航天航空、化工机械等领域。当湿空气流经温度比其露点温度低的表面时，会发生水蒸汽在冷表面凝结现象，若冷表面的温度低于冰点温度时，结霜现象就会发生。

在很多领域，结霜是有害的，在制冷系统、低温系统、热泵等工业设备中常出现结霜而影响设备的正常运行。在换热领域，霜层的存在会恶化传热性能。在结霜状况下的制冷系统换热能力下降35%以上。此外，对于很多内部有流体流动的设备，霜层的存在还会阻碍流体在设备内正常流动，使得流动阻力大大增加。

对结霜现象深入分析，可以做到扬长避短，充分利用其有利方面，并尽可能的削弱其负面影响。对于目前很多工程应用而言，预测换热设备在结霜工况下的运行特性是设计的基础。设备运行中的除霜工作是保证正常运行的必须措施，合理的制定除霜周期，可以保证设备的正常运行，最大限度的降低能耗，提高设备运行性能。

目前对制冷系统的仿真模拟计算技术主要针对两个领域开展研究：

一个领域是针对干工况情况下的蒸发换热情况，吴业正著《小型制冷装置设计指导》（第一版，北京：机械工业出版社，1998）介绍了一种制冷循环仿真计算方法，但这种方法的缺陷在于：并不涉及低温环境中不可避免的结霜问题。对于冷冻冷藏中出现结霜的情况，该方法与实际情况差异极大，所以这种方法并不能如实反映实际问题。

另一个领域是针对结霜问题，《一种结霜工况下翅片管蒸发器的性能计算方法》介绍了具体的单一针对蒸发器的仿真计算方法。但该方法将蒸发器独立于制冷系统之外，对蒸发器由于结霜对整个制冷系统的影响没有进行考虑，并不能真实有效地反映实际工作中各设备运转状况，不能够达到设计人员的要求。

发明内容

为解决现有技术存在的上述问题，本发明设计了一种既考虑结霜工况又能反映制冷系统随结霜变化情况的结霜工况下的制冷系统计算机仿真性能计算方法。

为了实现上述目的，本发明的技术方案如下：一种结霜工况下的制冷系统计算机仿真性能计算方法，包括以下步骤：

A、分别建立湿空气物性计算模块、冷媒物性计算模块、冷凝器结构计算模块、蒸发器结构计算模块和蒸发器结霜计算模块

A1、建立湿空气物性计算模块

利用Visual C++程序工具按美国采暖制冷与空调工程师协会提供的湿空气计算公式，建立湿空气物性计算模块；所述的湿空气物性包括湿空气的干球温度、湿球温度、含湿量、露点温度、焓、密度、比热和粘度；

A2、建立冷媒物性计算模块

利用美国制冷标准协会提供的制冷剂物性参数，建立DLL动态链接库；通过DLL文件与Visual C++程序的链接来建立冷媒物性计算模块；所述的冷媒物性包括冷媒的温度、压力、焓、熵、比热、密度、比容、粘度和表面张力；

A3、建立冷凝器结构计算模块

根据冷凝器结构参数建立冷凝器结构计算模块；所述的冷凝器结构参数见表1；

表1采用翅片管式冷凝器的冷凝器结构参数

翅片管外径(mm)	翅片型式：平片
		翅片管列数(列)	翅片间距(mm)
翅片管排数(排)	翅片厚度(mm)
		翅片管分流支路数(个)	风扇电机数量(个)
翅片管有效长度(mm)	冷凝器风扇风量

翅片管纵向间距(mm)	翅片管横向间距(mm)
		风扇电机输入功率(W/台)

A4、建立蒸发器结构计算模块

根据蒸发器结构参数建立蒸发器结构参数模块，所述的蒸发器结构参数见表2；

表2采用翅片管式冷风机的蒸发器结构参数

翅片管外径(mm)	翅片型式：平片
		翅片管列数(列)	翅片间距(mm)
翅片管排数(排)	翅片厚度(mm)
		翅片管分流支路数(个)	风扇电机数量(个)
翅片管有效长度(mm)	蒸发器风扇电机曲线
		翅片管纵向间距(mm)	翅片管横向间距(mm)
风扇电机输入功率(W/台)

表2所述的数据作为已知条件由用户输入确定；利用Visual C++语言将表2所述的数据编译为可调用的蒸发器结构计算模块；

A5、建立蒸发器结霜计算模块

A501、设置蒸发器初始工况

所述的蒸发器初始工况包括蒸发器的结霜周期、计算时间步长、初始的霜层密度、初始的霜层厚度和迭代收敛条件；

结霜周期：即整个蒸发器从开始工作至开始进行除霜的总时间，单位为h；

计算时间步长：由用户将结霜周期平均划分成若干份，其中每一份时间规定为一个时间步长；同一时间步长中的蒸发温度、冷凝温度、冷媒流量、霜层厚度和霜层密度为稳态参数；

初始的霜层密度：即开始结霜时刻霜层的密度大小，单位为kg/m³；

初始的霜层厚度：即开始结霜时刻霜层的厚度大小，单位为mm；

迭代收敛条件：即假设值与实际计算参数值之间的差异；

A502、根据用户输入的蒸发器结构参数，调用蒸发器结构计算模块，并进行蒸发器计算区间的划分；所述计算区间的划分方法分两种情况：沿冷媒流动方向，按等管长均匀划分计算区间；沿空气流动方向，按照管路行间距与列间距将翅片划分为相对独立的若干个计算区间；

A503、按照冷媒流动的管路的顺序，用户在人机交互程序界面中输入一个管路编号的顺序数组，所述的管路编号通过Visual C++界面设定程序自动生成；

利用Visual C++文件读入读出程序将所述的顺序数组保存为可读取的固定数据结构，通过读取所述的数据结构，在程序内部确定换热管连接关系，并按照所述的数据结构确定计算区间的计算顺序；

A504、假设迭代计算前的焓差和空气进出口温度差

由于单一计算区间焓差和空气进出口温度差无法直接计算求出，需要通过迭代计算的方式进行求解；因此，首先设定所述焓差和空气进出口温度差的假设值，通过迭代计算得到所述焓差和空气进出口温度差的计算值，当所述焓差和空气进出口温度差的假设值与计算值的差异小于收敛迭代条件时，则将所述焓差和空气进出口温度差的计算值作为最终计算结果；

同时根据所述焓差和空气进出口温度差的假设值计算出每个单一计算区间空气的含湿量的差；

A505、根据蒸发器结构计算模块和霜层厚度的初始设定值确定蒸发器的静压损失，同时利用风扇静压特性曲线计算此时流经蒸发器的风量，计算公式如下：

$\mathrm{\Δ}{p}_a=5.88\×{10}^{-4}{N}_{\mathrm{rd}}\×\frac{\frac{2}{p_t}[S_2-\frac{\mathrm{\π}{(d_0+2{\mathrm{\δ}}_{\mathrm{fr}})}^2}{{4S}_1}]+\frac{\mathrm{\π}(d_0+2{\mathrm{\δ}}_{\mathrm{fr}})}{S_1}}{{S_2}^{0.3}}$

$\×{\left[\frac{p_t\·S_1}{[p_t-(t_f+2{\mathrm{\δ}}_{\mathrm{fr}})][S_1-(d_0+2{\mathrm{\δ}}_{\mathrm{fr}})]}\right]}^3\×{W_F}^{1.7}$

Δp_a：空气侧静压损失，单位Pa；

N_ri：管排数；

d₀：管外径；

δ_fr：霜层厚度，单位m；

S₁：管列间距，单位m；

S₂：管行间距，单位m；

P_t：翅片间距，单位m；

t_f：翅片厚度，单位m；

W_F：迎面风速，单位m/s；

调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算；

A506、根据冷媒的假设的焓差，计算冷媒侧单一计算区间的换热量；

$Q_{\mathrm{iref}}=\stackrel{\·}{m}({\mathrm{hi}}_{\mathrm{out}}-{\mathrm{hi}}_{\mathrm{in}})$

Q_iref：冷媒侧单一计算区间的换热量，单位kW；

hi_out：单一计算区间的出口焓值，单位kJ/kg；

hi_in：单一计算区间的入口焓值，单位kJ/kg；

m：质量流量，单位kg/s；

A507、进行冷媒侧局部换热系数计算；

冷媒在蒸发器中流动，是一个蒸发过程，冷媒在蒸发器中具有三种不同的相态：冷媒在蒸发器中为液态时，所在的相区为液相区；冷媒在蒸发器中为气液混合物时，所在的相区为两相区；冷媒在蒸发器中为气态时，所在的相区为气相区；液相区与气相区由于都为单一相态，故称为单相区；

在蒸发过程中，冷媒的变化为液体—气液混合物—气体的过程；在不同的相态下其计算方法完全不同；

单相区的蒸发局部换热系数计算公式如下：

$a_{\mathrm{ir}}=\frac{{\mathrm{\λ}}_i}{d_i}\·0.023\·{{\mathrm{Re}}_g}^{0.8}\·{{\mathrm{Pr}}_g}^{0.4}$

Re_g：气相雷诺数；

Pr_g：气相普朗特数；

λ_i：液相导热系数；

d_i：管内径；

两相区的蒸发局部换热系数计算公式如下：

a_ir=Fa_cv+Sa_pb

a_cv为强制对流项；

F：强制对流系数；

a_pb为沸腾项；

S：沸腾系数；

调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算；

A508、由导热方程计算管壁温度：

$Q_{\mathrm{iref}}=a_{\mathrm{ir}}{A}_i(T_{\mathrm{iw}}-\frac{T_{\mathrm{irin}}+T_{\mathrm{irout}}}{2})$

Q_iref：冷媒侧单一计算区间的换热量，单位kW；

a_ir：空气侧换热系数；

A_i：管外换热面积；

T_iw：单一计算区间内管壁温度，单位℃；

T_irin：单一计算区间入口冷媒温度，单位℃；

T_irout：单一计算区间出口冷媒温度，单位℃；

调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算；

A509、由结霜判断条件判断单一计算区间是否结霜；

所述的结霜判断条件为同时满足以下两个条件时：计算所得的管壁温度低于0℃和单一计算区间的湿空气的露点温度高于管壁温度；

如果单一计算区间满足所述的结霜判断条件，则为结霜工况，转入步骤A510进行计算；

如果单一计算区间不满足所述的结霜判断条件，则为不结霜工况，确定霜层厚度变化值为0，空气侧含湿量差为0；转入步骤A514进行计算；

A510、计算霜层表面温度的初值

$T_{\mathrm{ifr}}=T_{\mathrm{iw}}+\frac{2\·Q_{\mathrm{iref}}\·m_{\mathrm{ir}}\·i_{\mathrm{sv}}}{2\·A_i\·l_{\mathrm{ifr}}}$

T_ifr：单一计算区间霜层表面温度，单位℃；

T_iw：单一计算区间内管壁温度，单位℃；

Q_iref：冷媒侧单一计算区间的换热量，单位kW；

m_ir：单一计算区间冷媒质量流量，单位kg/s；

i_ir：冷媒汽化潜热，单位kJ/kg；

A_i：管外换热面积，单位m²；

l_ifr：单一计算区间的管长，单位m；

式中i_sv为水和冰的相变潜热，A_T为计算区域的换热面积，m_r为导致霜层密度变化的质量，l_ifr为单一计算区间内的霜层的导热系数；

A511、计算空气侧的换热量，所述的换热量包括显热换热量和潜热换热量；

空气侧的显热换热量计算公式为

Q_isen=m_a·Cp_a(T_ain-T_aout)

Q_isen：单一计算区间空气侧显热换热量，单位kW；

m_a：单一计算区间空气流量，单位kg/s；

Cp_a：单一计算区间空气比热，单位kJ/kg K；

T_ain：单一计算区间空气入口温度，单位℃；

T_aout：单一计算区间空气出口温度，单位℃；

空气侧的潜热换热量计算公式为

Q_ilat=m_a·(d_in-d_out)

Q_ilat：单一计算区间空气侧潜热换热量，单位kW；

d_in：单一计算区间空气入口含湿量，单位kg；

d_out：单一计算区间空气出口含湿量，单位kg；

调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算；

A512、蒸发器单一计算区间换热量的迭代计算

由于空气侧换热量与冷媒侧换热量必须相等，比较两个换热量大小，如果其差异小于计算收敛精度10^-4时，则认为步骤A504假设的焓差和单一计算区间空气进出口的温度差为准确值，这样就完成了单一计算区间换热量的计算，否则回到步骤A504重新调整所述焓差和空气进出口温度差的假设值进行计算，直到满足收敛迭代条件为止；

A513、蒸发器整体计算区间计算

完成了一个单一计算区间的计算后，在冷媒流动方向上，将上一个计算区间的出口冷媒的温度、焓差和压力损失作为下一个计算区间的入口条件；在空气流动方向上，将上一个计算区间的空气出口温度、含湿量作为下一个计算区间的入口条件；

按照步骤A504-A512的步骤进行计算；直到将步骤A502中划分的所有计算区间计算值全部求得；即得到整个蒸发器每一个计算区域的冷媒的焓差、压力损失、空气侧的温差、含湿量差以及霜层的密度和厚度的变化；

A514、将上一个时间步长的计算的霜层厚度作为下一个时间步长的初始霜层厚度，由于每个单一区间的霜层厚度计算结果不同，将对不同的计算区间分别赋予相应的初始霜层厚度；

A515、将计算所得的霜层厚度、霜层密度、蒸发器换热量与压力损失、空气侧温度差与含湿量差作为已知数据进行保存，供制冷循环计算模块调用；

B、建立制冷循环系统计算模块

B1、设置蒸发器的结霜周期、计算时间步长、初始的霜层密度和迭代收敛条件，所述的迭代收敛条件的默认值为10^-4；

所述的结霜周期为蒸发器从开始工作至开始进行除霜的总时间，单位为h；

计算时间步长：由用户将结霜周期平均划分成若干份，其中每一份时间规定为一个时间步长；同一时间步长中的蒸发温度、冷凝温度、冷媒流量、霜层厚度和霜层密度为稳态参数；

所述的时间步长由用户将结霜周期进行平均划分，其中每一份时间规定为一个时间步长；同一时间步长中的蒸发温度、冷凝温度、冷媒流量、霜层厚度和霜层密度为稳态参数；

所述的初始的霜层密度为开始结霜时刻霜层的密度大小，单位为kg/m³；

所述的迭代收敛条件为假设参数值与实际计算参数值之间的差异；

B2、设定运行参数，根据用户针对不同的工作条件设定蒸发器过热度、冷凝器过冷度、吸气压力损失和排气压力损失；同时设定室内侧空气温度与湿度、室外侧空气温度与湿度；

蒸发器过热度：冷媒在蒸发器内吸收热量；冷媒从液态至两相，在两相状态下的温度为蒸发温度，在两相状态下冷媒温度不变；又从两相变为气相，温度升高。从蒸发温度至蒸发器出口时的气相温度差为蒸发器过热度；蒸发器过热度作为已知条件由用户输入确定，单位为℃；

冷凝器过冷度：冷媒在冷凝器内吸收热量；冷媒从气态至两相，在两相状态下的温度为冷凝温度，在两相状态下冷媒温度不变；又从两相变为液相，温度下降；从冷凝温度至冷凝器出口时的液相温度差为冷凝器过冷度；冷凝器过冷度作为已知条件由用户输入确定，单位为℃；

吸气压力损失：制冷设备从蒸发器至压缩机吸气口会有一段的管路连接，在管路内由于摩擦力作用，压力会有部分损失，导致冷媒压力下降，温度下降；另外，由于管路与外界的热交换也会产生一部分损失；这些损失统称为吸气压力损失；吸气压力损失用温度来表征，作为已知条件由用户输入确定，单位为℃；

排气压力损失：制冷设备从压缩机排气口至冷凝器会有一段的管路连接，在管路内由于摩擦力作用，压力会有部分损失，导致冷媒压力下降，温度下降；另外，由于管路与外界的热交换也会产生一部分损失；这些损失统称为排气压力损失；该排气压力损失用温度来表征，作为已知条件由用户输入确定，单位为℃；

B3、首先假定蒸发温度及冷凝温度，计算出冷凝器热负荷，通过冷凝器热负荷确定冷凝器所需的换热面积A₁；同时调用冷凝器结构计算模块以及由用户输入冷凝器侧空气入口干球温度、含湿量和冷凝器风扇风量确定实际冷凝器换热面积A₂；冷凝器所需的换热面积A₁与实际冷凝器换热面积A₂进行迭代循环后确定冷凝器实际换热面积；再利用冷凝器热负荷、室外侧空气入口干球温度、含湿量和冷凝器风扇风量，确定冷凝器冷凝温度T_c；同时确定压缩机的制冷能力、蒸发器换热量Q₂、输入功率以及电流；

$A_1=\frac{Q_2}{k_0\·{\mathrm{\θ}}_m}$

k₀:冷凝器总换热系数，单位W/(m²·K)；

θ_m:对数平均温差；

${\mathrm{\θ}}_m=\frac{t_{a1}-t_{a2}}{\ln \frac{t_k-t_{a1}}{t_k-t_{a2}}}$

t_a1：空气入口干球温度；

t_a2：出口空气干球温度；

t_k：冷凝温度；

根据冷凝器结构计算模块确定冷凝器实际换热面积A₂；在冷凝器所需的换热面积A₁与冷凝器实际换热面积A₂不相等时，调整假设的冷凝温度使两者相等，此时的冷凝温度为所述的假设的蒸发温度下的真实冷凝温度；但由于蒸发温度目前仍然为假设值，故此冷凝温度尚不是循环系统中的真实冷凝温度；

由结构确定的实际冷凝器换热面积A₂：A₂=(m-1)πD_tL_COREn_dn_r

D_t：翅片管外径

L_CORE：单一计算区间的管程；

n_d：管段数；

n_r：管列数；

m：面积扩大率m=A_f(π·d_ci·P_f)；

d_ci：管内径；

P_f：翅片间距；

B4、根据确定的冷凝温度T_c、设定的冷凝器过冷度和假定的蒸发温度，确定蒸发器入口冷媒温度、冷媒压力、冷媒的流量以及蒸发器热负荷Q₂；

调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算；

B5、调用蒸发器结霜计算模块进行结霜工况下的蒸发器仿真计算；

B6、利用步骤B5得到的蒸发器在结霜工况下实际的蒸发器换热量Q1和步骤B3得到的蒸发器换热量Q2进行迭代计算，从而确定实际蒸发器侧换热量；同时确定蒸发器结霜厚度和结霜密度；最终确定实际的蒸发温度；完成一个时间步长内的计算结果，利用Visual C++数据传递工具，将计算结果输出到Excel表格中；

B7、将上一个时间步长的计算结果作为下一个时间步长的计算初值，并更新整个计算区域；

B8、当满足计算结束的判断条件时，计算结束；

所述的计算结束的判断条件为：计算时间达到设置的结霜周期或计算的霜层厚度大于等于翅片间距的一半；

计算结束后，将各个时间步长的平均结霜厚度、平均结霜密度、冷凝温度、蒸发温度、制冷量、蒸发器热负荷、冷凝器热负荷、能效比COP和冷媒循环量输入到Excel表格，并自动绘制图表。

与现有技术相比，本发明具有以下有益效果：

1、本发明利用计算机辅助技术方法来实现，对蒸发器结霜情况进行数值模拟。首先，蒸发器按照一定的管长分割为若干计算区域，使用有限容积方法建立蒸发器结霜工况的计算模型，并且用C++编程语言编制仿真计算程序对每个计算区域的控制方程进行反复迭代计算求解，模拟蒸发器的霜层增长状况。本发明可以准确的反应霜层在蒸发器不同区域的结霜程度，可以有效地指导蒸发器的结构设计，从而有效地避免局部结霜过多而导致的制冷量下降以及换热器利用率不高的问题。

2、本发明在模拟了蒸发器霜层变化的基础上，蒸发温度的变化会导致整个制冷循环工作状况的变化。在一个计算步长结束后，霜层厚度变化，会导致下一个计算步长中蒸发温度的变化，故对整个制冷系统工作状态都有影响。动态地反映了制冷系统的关键参数如蒸发温度、冷凝温度、冷媒流量、制冷量等，随着霜层变化而产生的相对变化。这种计算方法准确模拟了一定时间下的制冷设备工作状态，可以有效地指导设计人员避免由于结霜而产生的系统非正常工作地风险。同时也为制冷系统的除霜控制提供有效地理论数据，更有效的提高系统的制冷量及能效比。

3、本发明在程序设计方法上首先将计算时间进行均匀步长分割。在每一个计算步长中，建立冷媒物性计算模块、湿空气物性计算模块、冷凝器结构计算模块、蒸发器结构计算模块、蒸发器结霜计算模块和制冷循环计算模块，通过制冷循环计算模块将各个模块相互联系起来。利用计算机辅助工具对程序内的结构与换热参数进行复杂的迭代计算。最终将迭代收敛的计算结果输出在Excel表格中，实现整个系统在结霜工况下的动态仿真。

附图说明

本发明共有附图6张，其中：

图1制冷系统模块关系总图。

图2蒸发器结霜计算模块流程图。

图3制冷循环计算模块流程图。

图4换热管连接关系设置图。

图5蒸发器计算区间划分关系示意图（冷媒流动方向）。

图6蒸发器计算区间划分关系示意图（空气流动方向）。

具体实施方式

为了使研究人员更为清晰准确地了解本发明计算方法。下面结合相关实际运行下的工况对本发明进行进一步的阐述。

图1所示为本发明的流程总图，具体步骤如下：

A、进行湿空气物性计算模块、冷媒物性计算模块、冷凝器结构计算模块、蒸发器结构计算模块和蒸发器结霜计算模块的建立。

A1、建立湿空气物性计算模块

以下是通过湿空气已知条件计算湿空气其他物性参数的程序方法：

已知湿空气的压力为101.325kPa，干球温度为10℃，相对湿度为60%。计算得出其他参数为湿球温度为6.46℃，焓值为21.54kJ/kg，含湿量为4.56g/kg，密度为1.242kg/m³。

例如：

d=CMoistAir::RH2W(doublet,doublerh){return0.621945*rh*Pws(t)/(m_Patm-rh*Pws(t));};

（rh为相对湿度，t为干球温度，d为含湿量）

该语句意思为输入t和rh参数，可得出含湿量d。

公式来源于ASHARE 2009 F01 SI Charper 1 Psychrometrics（美国采暖、制冷与空调工程师学会2009第一章湿度）

通过湿空气物性计算模块，可以根据已知湿空气的状态参数，求得计算所需的其他的湿空气相关物性参数。包括：干球温度、湿球温度、含湿量、露点温度、焓、密度、比热、粘度等。

A2、建立冷媒计算模块

通过调用DLL程序方法，可以根据已知冷媒的状态参数，求得计算所需的其他相关物性参数。包括：温度、压力、焓、熵、密度、比热、粘度、表面张力等。

例如：已知制冷剂的压力为465kPa，焓为220kJ/kg；

A.Pk=465；A.H=220；(A，B为接口类)

CLoaderMgr::getInstance(R22).getLoader()->getPHFL1(A,B);

ρ=B.RHOMOLL;

T=B.Tk;

计算得出密度为12.63kg/m³；

A3、建立冷凝器结构计算模块

冷凝器结构参数：采用翅片管式冷凝器

计算的相关参数见表3。

表3翅片管式冷凝器结构参数数值表

翅片管外径(mm)：7.2	翅片管内径(mm)：6.6
		翅片管列数(列)：2	翅片型式：平片
翅片管排数(排)：56	翅片间距(mm)：3
		翅片管分流支路数(个)：14	翅片厚度(mm)：0.12
翅片管有效长度(mm)：1125	风扇电机数量(个)：2
		翅片管纵向间距(mm)：21	冷凝器风扇电机曲线
翅片管横向间距(mm)：18.2	风扇电机输入功率：460(W/台)

A4、建立蒸发器结构计算模块

采用翅片管式冷风机的蒸发器结构参数数值见表4。

表4采用翅片管式冷风机的蒸发器结构参数数值表

翅片管外径(mm)：9.52	翅片管内径(mm)：6.6
		翅片管列数(列)：6	翅片型式：平片
翅片管排数(排)：18	翅片间距(mm)：10
		翅片管分流支路数(个)：9	翅片厚度(mm)：0.1115
翅片管有效长度(mm)：1080	风扇电机数量(个)：2
		翅片管纵向间距(mm)：25	蒸发器风扇电机曲线
翅片管横向间距(mm)：21.6	风扇电机输入功率：460(W/台)

A5、建立蒸发器结霜计算模块，其流程如附图2所示。

A501、设置蒸发器的结霜周期，计算时间步长和初始的霜层密度，迭代收敛条件（默认为10^-4）；

蒸发器从工作到开始融霜，该事件总长为120分钟。时间步长设定为5分钟，则每隔5分钟得到1组相关的计算数据，整个计算时间共13组计算数据，同时在这一个时间步长中的各个参数为稳态。如：蒸发温度、冷凝温度、冷媒流量等保持不变。

该算例设定为120min。

时间步长：。该算例设定为5分钟。

初始的霜层密度：即开始结霜时刻，霜层的密度大小。设定为80kg/m³。

在本发明中计算收敛条件默认为10^-4。例如，假设蒸发器壁温与实际计算的蒸发器壁温差异小于10^-4，则认为壁温计算结果为真。

A502、根据用户输入的蒸发器结构参数，调用蒸发器结构计算模块，并进行蒸发器计算区间的划分。

按照图5所示进行冷媒流动方向蒸发器计算区间划分；按照图6所示进行空气流动方向蒸发器计算区间划分。

A503、冷媒管路的流路配置

如图4所示，将蒸发器管路分为7个流路。编号为①—⑦。以第一个流路为例，共占用三排管，管排编号为S1—S3；共占用三列，列编号为R1—R3。按图4所示，用户只要输入S2R1—S1R1—S2R2—S3R2—S3R3—S2R3。程序将读取用户要求的管路连接顺序，计算时将按照所述的管路连接顺序完成冷媒计算。

A504、假设迭代计算前的焓差和空气进出口温度差。

A505、利用风扇静压特性曲线计算此时流经蒸发器的风量。

由发明步骤A505中的静压计算公式，确定蒸发器静压；再通过风扇的静压特性曲线来确定风扇风量。一般该静压特性曲线由蒸发器厂家提供，该曲线主要是体现静压与风量之间的关系。

A506、根据冷媒的假设的焓差，计算冷媒侧单一计算区间的换热量。

A507、进行冷媒侧局部换热系数计算。

A508、计算管壁温度。

A509、判断单一计算区间是否结霜。

A510、计算霜层表面温度的初值。

A511、计算空气侧的换热量，所述的换热量包括显热换热量和潜热换热量。

A512、蒸发器单一计算区间换热量的迭代计算。

A513、蒸发器整体计算区间计算。

A514、确定下一个时间步长的初始霜层厚度。

A515、将计算所得的霜层厚度、霜层密度、蒸发器换热量与压力损失、空气侧温度差与含湿量差作为已知数据进行保存，供制冷循环计算模块调用。

B、按图3所示，建立制冷循环系统计算模块

B1、设置蒸发器的结霜周期、计算时间步长、初始的霜层密度和迭代收敛条件，所述的迭代收敛条件的默认值为10^-4；

B2设定工况参数：

包括蒸发器过热度、冷凝器过冷度、吸气压力损失和排气压力损失，室内侧空气温度与湿度，室外侧空气温度与湿度。

室内侧（蒸发器侧）空气入口温度，湿度：0℃，70%

室外侧（冷凝器侧）空气入口温度，湿度：32℃，60%

蒸发器过热度：5℃

冷凝器过冷度：3℃

吸气压力损失：0.2℃

排气压力损失：0.3

B3、通过假设的蒸发温度与冷凝温度，以及根据循环中所用压缩机的性能曲线，可以确定压缩机的制冷能力，蒸发器换热量Q₂，输入功率以及电流。同时还可以确定冷凝器的热负荷。通过对数平均温差法可求的冷凝器所需的换热面积A₁。

B4、确认该假设蒸发温度下的真实冷凝温度后，同时可以确定该蒸发温度下的蒸发器入口处的冷媒相关参数，蒸发器入口冷媒温度，压力；冷媒流量。

B5、以步骤B4中提供的参数作为蒸发器入口参数，进入蒸发器结霜计算模块。在已知蒸发器入口参数以及蒸发器结构模块确定的蒸发器结构参数条件下，蒸发器结霜仿真计算方法，按照步骤A5的相关计算公式，能够确定以下参数：

结霜条件下蒸发器换热量Q₁、结霜厚度δ₁、结霜密度ρ₁等。由于蒸发温度目前仍然为假设值，故此蒸发器换热量Q₁、结霜厚度、结霜密度尚不是循环系统中的真实值。

B6、判断步骤B5所得的蒸发器换热量Q₁与步骤B3中所得蒸发器换热量Q₂是否相等，如不相等，转到到步骤B5重新设定新的蒸发温度，直到实现Q₁与Q₂相等，此时蒸发温度为制冷循环系统的真实值，在该蒸发温度条件下，步骤B6与步骤B5中计算得到的所有参数都确定为真实值。

B7、完成一个时间步长的迭代计算，将所得结果输入到Excel表格中保存。同时由于结霜厚度与结霜密度的变化，在下一个时间步长开始时刻的相关参数，（如：蒸发器翅片间距、翅片厚度等）都需要进行数据更新。完成更新，进入下一个时间步长计算，计算流程重复步骤B1至步骤B6。

例如：0时刻结霜厚度为零，第一个时间步长结束后结霜厚度增加为1mm，则第二个时间步长初始时刻的结霜厚度为1mm，并由结霜厚度的变化使整个蒸发器结构参数发生相应的变化。

B8、当计算时间结束或者结霜厚度大于翅片间距一半（意味着霜层已经将蒸发器完全堵塞）时，计算完成。将数据保存入Excel表格中。表5为本实施例的最终计算结果。

表5结霜周期120分钟系统仿真模拟计算结果

Claims (1)

1.一种结霜工况下的制冷系统计算机仿真性能计算方法，其特征在于：包括以下步骤：

A、分别建立湿空气物性计算模块、冷媒物性计算模块、冷凝器结构计算模块、蒸发器结构计算模块和蒸发器结霜计算模块

A1、建立湿空气物性计算模块

利用Visual C++程序工具按美国采暖制冷与空调工程师协会提供的湿空气计算公式，建立湿空气物性计算模块；所述的湿空气物性包括湿空气的干球温度、湿球温度、含湿量、露点温度、焓、密度、比热和粘度；

A2、建立冷媒物性计算模块

利用美国制冷标准协会提供的制冷剂物性参数，建立DLL动态链接库；通过DLL文件与Visual C++程序的链接来建立冷媒物性计算模块；所述的冷媒物性包括冷媒的温度、压力、焓、熵、比热、密度、比容、粘度和表面张力；

A3、建立冷凝器结构计算模块

根据冷凝器结构参数建立冷凝器结构计算模块；所述的冷凝器结构参数见表1；

表1采用翅片管式冷凝器的冷凝器结构参数

  翅片管外径(mm)   翅片型式：平片   翅片管列数(列)   翅片间距(mm)   翅片管排数(排)   翅片厚度(mm)   翅片管分流支路数(个)   风扇电机数量(个)   翅片管有效长度(mm)   冷凝器风扇风量   翅片管纵向间距(mm)   翅片管横向间距(mm)   风扇电机输入功率(W/台)

A4、建立蒸发器结构计算模块

根据蒸发器结构参数建立蒸发器结构参数模块，所述的蒸发器结构参数见表2；

表2采用翅片管式冷风机的蒸发器结构参数

  翅片管外径(mm)   翅片型式：平片

  翅片管列数(列)   翅片间距(mm)   翅片管排数(排)   翅片厚度(mm)   翅片管分流支路数(个)   风扇电机数量(个)   翅片管有效长度(mm)   蒸发器风扇电机曲线   翅片管纵向间距(mm)   翅片管横向间距(mm)   风扇电机输入功率(W/台)

表2所述的数据作为已知条件由用户输入确定；利用Visual C++语言将表2所述的数据编译为可调用的蒸发器结构计算模块；

A5、建立蒸发器结霜计算模块

A501、设置蒸发器初始工况

所述的蒸发器初始工况包括蒸发器的结霜周期、计算时间步长、初始的霜层密度、初始的霜层厚度和迭代收敛条件；

结霜周期：即整个蒸发器从开始工作至开始进行除霜的总时间，单位为h；

计算时间步长：由用户将结霜周期平均划分成若干份，其中每一份时间规定为一个时间步长；同一时间步长中的蒸发温度、冷凝温度、冷媒流量、霜层厚度和霜层密度为稳态参数；

初始的霜层密度：即开始结霜时刻霜层的密度大小，单位为kg/m³；

初始的霜层厚度：即开始结霜时刻霜层的厚度大小，单位为mm；

迭代收敛条件：即假设值与实际计算参数值之间的差异；

A502、根据用户输入的蒸发器结构参数，调用蒸发器结构计算模块，并进行蒸发器计算区间的划分；所述计算区间的划分方法分两种情况：沿冷媒流动方向，按等管长均匀划分计算区间；沿空气流动方向，按照管路行间距与列间距将翅片划分为相对独立的若干个计算区间；

A503、按照冷媒流动的管路的顺序，用户在人机交互程序界面中输入一个管路编号的顺序数组，所述的管路编号通过Visual C++界面设定程序自动生成；

利用Visual C++文件读入读出程序将所述的顺序数组保存为可读取的固定数据结构，通过读取所述的数据结构，在程序内部确定换热管连接关系，并按照所述的数据结构确定计算区间的计算顺序；

A504、假设迭代计算前的焓差和空气进出口温度差

由于单一计算区间焓差和空气进出口温度差无法直接计算求出，需要通过迭代计算的方式进行求解；因此，首先设定所述焓差和空气进出口温度差的假设值，通过迭代计算得到所述焓差和空气进出口温度差的计算值，当所述焓差和空气进出口温度差的假设值与计算值的差异小于收敛迭代条件时，则将所述焓差和空气进出口温度差的计算值作为最终计算结果；

同时根据所述焓差和空气进出口温度差的假设值计算出每个单一计算区间空气的含湿量的差；

A505、根据蒸发器结构计算模块和霜层厚度的初始设定值确定蒸发器的静压损失，同时利用风扇静压特性曲线计算此时流经蒸发器的风量，计算公式如下：

$\mathrm{\Δ}{p}_a=5.88\×{10}^{-4}{N}_{\mathrm{rd}}\×\frac{\frac{2}{p_t}[S_2-\frac{\mathrm{\π}{(d_0+2{\mathrm{\δ}}_{\mathrm{fr}})}^2}{{4S}_1}]+\frac{\mathrm{\π}(d_0+2{\mathrm{\δ}}_{\mathrm{fr}})}{S_1}}{{S_2}^{0.3}}$

$\×{\left[\frac{p_t\·S_1}{[p_t-(t_f+2{\mathrm{\δ}}_{\mathrm{fr}})][S_1-(d_0+2{\mathrm{\δ}}_{\mathrm{fr}})]}\right]}^3\×{W_F}^{1.7}$

Δp_a：空气侧静压损失，单位Pa；

N_ri：管排数；

d₀：管外径；

δ_fr：霜层厚度，单位m；

S₁：管列间距，单位m；

S₂：管行间距，单位m；

P_t：翅片间距，单位m；

t_f：翅片厚度，单位m；

W_F：迎面风速，单位m/s；

调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算；

A506、根据冷媒的假设的焓差，计算冷媒侧单一计算区间的换热量；

$Q_{\mathrm{iref}}=\stackrel{\·}{m}({\mathrm{hi}}_{\mathrm{out}}-{\mathrm{hi}}_{\mathrm{in}})$

Q_iref：冷媒侧单一计算区间的换热量，单位kW；

hi_out：单一计算区间的出口焓值，单位kJ/kg；

hi_in：单一计算区间的入口焓值，单位kJ/kg；

m：质量流量，单位kg/s；

A507、进行冷媒侧局部换热系数计算；

冷媒在蒸发器中流动，是一个蒸发过程，冷媒在蒸发器中具有三种不同的相态：冷媒在蒸发器中为液态时，所在的相区为液相区；冷媒在蒸发器中为气液混合物时，所在的相区为两相区；冷媒在蒸发器中为气态时，所在的相区为气相区；液相区与气相区由于都为单一相态，故称为单相区；

在蒸发过程中，冷媒的变化为液体—气液混合物—气体的过程；在不同的相态下其计算方法完全不同；

单相区的蒸发局部换热系数计算公式如下：

$a_{\mathrm{ir}}=\frac{{\mathrm{\λ}}_i}{d_i}\·0.023\·{{\mathrm{Re}}_g}^{0.8}\·{{\mathrm{Pr}}_g}^{0.4}$

Re_g：气相雷诺数；

Pr_g：气相普朗特数；

λ_i：液相导热系数；

d_i：管内径；

两相区的蒸发局部换热系数计算公式如下：

a_ir=Fa_cv+Sa_pb

a_cv为强制对流项；

F：强制对流系数；

a_pb为沸腾项；

S：沸腾系数；

调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算；

A508、由导热方程计算管壁温度：

$Q_{\mathrm{iref}}=a_{\mathrm{ir}}{A}_i(T_{\mathrm{iw}}-\frac{T_{\mathrm{irin}}+T_{\mathrm{irout}}}{2})$

Q_iref：冷媒侧单一计算区间的换热量，单位kW；

a_ir：空气侧换热系数；

A_i：管外换热面积；

T_iw：单一计算区间内管壁温度，单位℃；

T_irin：单一计算区间入口冷媒温度，单位℃；

T_irout：单一计算区间出口冷媒温度，单位℃；

调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算；

A509、由结霜判断条件判断单一计算区间是否结霜；

所述的结霜判断条件为同时满足以下两个条件时：计算所得的管壁温度低于0℃和单一计算区间的湿空气的露点温度高于管壁温度；

如果单一计算区间满足所述的结霜判断条件，则为结霜工况，转入步骤A510进行计算；

如果单一计算区间不满足所述的结霜判断条件，则为不结霜工况，确定霜层厚度变化值为0，空气侧含湿量差为0；转入步骤A514进行计算；

A510、计算霜层表面温度的初值

$T_{\mathrm{ifr}}=T_{\mathrm{iw}}+\frac{2\·Q_{\mathrm{iref}}\·m_{\mathrm{ir}}\·i_{\mathrm{sv}}}{2\·A_i\·l_{\mathrm{ifr}}}$

T_ifr：单一计算区间霜层表面温度，单位℃；

T_iw：单一计算区间内管壁温度，单位℃；

Q_iref：冷媒侧单一计算区间的换热量，单位kW；

m_ir：单一计算区间冷媒质量流量，单位kg/s；

i_ir：冷媒汽化潜热，单位kJ/kg；

A_i：管外换热面积，单位m²；

l_ifr：单一计算区间的管长，单位m；

式中i_sv为水和冰的相变潜热，A_T为计算区域的换热面积，m_r为导致霜层密度变化的质量，l_ifr为单一计算区间内的霜层的导热系数；

A511、计算空气侧的换热量，所述的换热量包括显热换热量和潜热换热量；

空气侧的显热换热量计算公式为

Q_isen=m_a·Cp_a(T_ain-T_aout)

Q_isen：单一计算区间空气侧显热换热量，单位kW；

m_a：单一计算区间空气流量，单位kg/s；

Cp_a：单一计算区间空气比热，单位kJ/kg K；

T_ain：单一计算区间空气入口温度，单位℃；

T_aout：单一计算区间空气出口温度，单位℃；

空气侧的潜热换热量计算公式为

Q_ilat=m_a·(d_in-d_out)

Q_ilat：单一计算区间空气侧潜热换热量，单位kW；

d_in：单一计算区间空气入口含湿量，单位kg；

d_out：单一计算区间空气出口含湿量，单位kg；

调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算；

A512、蒸发器单一计算区间换热量的迭代计算

由于空气侧换热量与冷媒侧换热量必须相等，比较两个换热量大小，如果其差异小于计算收敛精度10^-4时，则认为步骤A504假设的焓差和单一计算区间空气进出口的温度差为准确值，这样就完成了单一计算区间换热量的计算，否则回到步骤A504重新调整所述焓差和空气进出口温度差的假设值进行计算，直到满足收敛迭代条件为止；

A513、蒸发器整体计算区间计算

完成了一个单一计算区间的计算后，在冷媒流动方向上，将上一个计算区间的出口冷媒的温度、焓差和压力损失作为下一个计算区间的入口条件；在空气流动方向上，将上一个计算区间的空气出口温度、含湿量作为下一个计算区间的入口条件；

按照步骤A504-A512的步骤进行计算；直到将步骤A502中划分的所有计算区间计算值全部求得；即得到整个蒸发器每一个计算区域的冷媒的焓差、压力损失、空气侧的温差、含湿量差以及霜层的密度和厚度的变化；

A514、将上一个时间步长的计算的霜层厚度作为下一个时间步长的初始霜层厚度，由于每个单一区间的霜层厚度计算结果不同，将对不同的计算区间分别赋予相应的初始霜层厚度；

A515、将计算所得的霜层厚度、霜层密度、蒸发器换热量与压力损失、空气侧温度差与含湿量差作为已知数据进行保存，供制冷循环计算模块调用；

B、建立制冷循环系统计算模块

B1、设置蒸发器的结霜周期、计算时间步长、初始的霜层密度和迭代收敛条件，所述的迭代收敛条件的默认值为10^-4；

所述的结霜周期为蒸发器从开始工作至开始进行除霜的总时间，单位为h；

计算时间步长：由用户将结霜周期平均划分成若干份，其中每一份时间规定为一个时间步长；同一时间步长中的蒸发温度、冷凝温度、冷媒流量、霜层厚度和霜层密度为稳态参数；

所述的时间步长由用户将结霜周期进行平均划分，其中每一份时间规定为一个时间步长；同一时间步长中的蒸发温度、冷凝温度、冷媒流量、霜层厚度和霜层密度为稳态参数；

所述的初始的霜层密度为开始结霜时刻霜层的密度大小，单位为kg/m³；

所述的迭代收敛条件为假设参数值与实际计算参数值之间的差异；

B2、设定运行参数，根据用户针对不同的工作条件设定蒸发器过热度、冷凝器过冷度、吸气压力损失和排气压力损失；同时设定室内侧空气温度与湿度、室外侧空气温度与湿度；

蒸发器过热度：冷媒在蒸发器内吸收热量；冷媒从液态至两相，在两相状态下的温度为蒸发温度，在两相状态下冷媒温度不变；又从两相变为气相，温度升高；从蒸发温度至蒸发器出口时的气相温度差为蒸发器过热度；蒸发器过热度作为已知条件由用户输入确定，单位为℃；

冷凝器过冷度：冷媒在冷凝器内吸收热量；冷媒从气态至两相，在两相状态下的温度为冷凝温度，在两相状态下冷媒温度不变；又从两相变为液相，温度下降；从冷凝温度至冷凝器出口时的液相温度差为冷凝器过冷度；冷凝器过冷度作为已知条件由用户输入确定，单位为℃；

吸气压力损失：制冷设备从蒸发器至压缩机吸气口会有一段的管路连接，在管路内由于摩擦力作用，压力会有部分损失，导致冷媒压力下降，温度下降；另外，由于管路与外界的热交换也会产生一部分损失；这些损失统称为吸气压力损失；吸气压力损失用温度来表征，作为已知条件由用户输入确定，单位为℃；

排气压力损失：制冷设备从压缩机排气口至冷凝器会有一段的管路连接，在管路内由于摩擦力作用，压力会有部分损失，导致冷媒压力下降，温度下降；另外，由于管路与外界的热交换也会产生一部分损失；这些损失统称为排气压力损失；该排气压力损失用温度来表征，作为已知条件由用户输入确定，单位为℃；

B3、首先假定蒸发温度及冷凝温度，计算出冷凝器热负荷，通过冷凝器热负荷确定冷凝器所需的换热面积A₁；同时调用冷凝器结构计算模块以及由用户输入冷凝器侧空气入口干球温度、含湿量和冷凝器风扇风量确定实际冷凝器换热面积A₂；冷凝器所需的换热面积A₁与实际冷凝器换热面积A₂进行迭代循环后确定冷凝器实际换热面积；再利用冷凝器热负荷、室外侧空气入口干球温度、含湿量和冷凝器风扇风量，确定冷凝器冷凝温度T_c；同时确定压缩机的制冷能力、蒸发器换热量Q₂、输入功率以及电流；

$A_1=\frac{Q_2}{k_0\·{\mathrm{\θ}}_m}$

k₀:冷凝器总换热系数，单位W/(m²·K)；

θ_m:对数平均温差；

${\mathrm{\θ}}_m=\frac{t_{a1}-t_{a2}}{\ln \frac{t_k-t_{a1}}{t_k-t_{a2}}}$

t_a1：空气入口干球温度；

t_a2：出口空气干球温度；

t_k：冷凝温度；

根据冷凝器结构计算模块确定冷凝器实际换热面积A₂；在冷凝器所需的换热面积A₁与冷凝器实际换热面积A₂不相等时，调整假设的冷凝温度使两者相等，此时的冷凝温度为所述的假设的蒸发温度下的真实冷凝温度；但由于蒸发温度目前仍然为假设值，故此冷凝温度尚不是循环系统中的真实冷凝温度；

由结构确定的实际冷凝器换热面积A₂：A₂=(m-1)πD_tL_COREn_dn_r

D_t：翅片管外径

L_CORE：单一计算区间的管程；

n_d：管段数；

n_r：管列数；

m：面积扩大率m=A_f(π·d_ci·P_f)；

d_ci：管内径；

P_f：翅片间距；

B4、根据确定的冷凝温度T_c、设定的冷凝器过冷度和假定的蒸发温度，确定蒸发器入口冷媒温度、冷媒压力、冷媒的流量以及蒸发器热负荷Q₂；

调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算；

B5、调用蒸发器结霜计算模块进行结霜工况下的蒸发器仿真计算；

B6、利用步骤B5得到的蒸发器在结霜工况下实际的蒸发器换热量Q1和步骤B3得到的蒸发器换热量Q2进行迭代计算，从而确定实际蒸发器侧换热量；同时确定蒸发器结霜厚度和结霜密度；最终确定实际的蒸发温度；完成一个时间步长内的计算结果，利用Visual C++数据传递工具，将计算结果输出到Excel表格中；

B7、将上一个时间步长的计算结果作为下一个时间步长的计算初值，并更新整个计算区域；

B8、当满足计算结束的判断条件时，计算结束；

所述的计算结束的判断条件为：计算时间达到设置的结霜周期或计算的霜层厚度大于等于翅片间距的一半；

计算结束后，将各个时间步长的平均结霜厚度、平均结霜密度、冷凝温度、蒸发温度、制冷量、蒸发器热负荷、冷凝器热负荷、能效比COP和冷媒循环量输入到Excel表格，并自动绘制图表。

Images