发明名称 一种结霜工况下的制冷系统计算机仿真性能计算方法
摘要 本发明公开了一种结霜工况下的制冷系统计算机仿真性能计算方法,包括以下步骤:分别建立湿空气物性计算模块、冷媒物性计算模块、冷凝器结构计算模块、蒸发器结构计算模块和蒸发器结霜计算模块;建立制冷循环系统计算模块。本发明在模拟了蒸发器霜层变化的基础上,蒸发温度的变化会导致整个制冷循环工作状况的变化。本发明可以有效地指导设计人员避免由于结霜而产生的系统非正常工作地风险。同时也为制冷系统的除霜控制提供有效地理论数据,更有效的提高系统的制冷量及能效比。本发明首先将计算时间进行均匀步长分割,利用计算机辅助工具对程序内的结构与换热参数进行复杂的迭代计算。本发明实现了整个系统在结霜工况下的动态仿真。
申请公布号 CN102779217A 申请公布日期 2012.11.14
申请号 CN201210277558.5 申请日期 2012.08.06
申请人 大连三洋压缩机有限公司 发明人 秦海杰;夏梦心;朱卫英
分类号 G06F17/50(2006.01)I 主分类号 G06F17/50(2006.01)I
代理机构 大连东方专利代理有限责任公司 21212 代理人 李洪福
主权项 一种结霜工况下的制冷系统计算机仿真性能计算方法,其特征在于:包括以下步骤: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/m3;初始的霜层厚度:即开始结霜时刻霜层的厚度大小,单位为mm;迭代收敛条件:即假设值与实际计算参数值之间的差异;A502、根据用户输入的蒸发器结构参数,调用蒸发器结构计算模块,并进行蒸发器计算区间的划分;所述计算区间的划分方法分两种情况:沿冷媒流动方向,按等管长均匀划分计算区间;沿空气流动方向,按照管路行间距与列间距将翅片划分为相对独立的若干个计算区间;A503、按照冷媒流动的管路的顺序,用户在人机交互程序界面中输入一个管路编号的顺序数组,所述的管路编号通过Visual C++界面设定程序自动生成;利用Visual C++文件读入读出程序将所述的顺序数组保存为可读取的固定数据结构,通过读取所述的数据结构,在程序内部确定换热管连接关系,并按照所述的数据结构确定计算区间的计算顺序;A504、假设迭代计算前的焓差和空气进出口温度差由于单一计算区间焓差和空气进出口温度差无法直接计算求出,需要通过 迭代计算的方式进行求解;因此,首先设定所述焓差和空气进出口温度差的假设值,通过迭代计算得到所述焓差和空气进出口温度差的计算值,当所述焓差和空气进出口温度差的假设值与计算值的差异小于收敛迭代条件时,则将所述焓差和空气进出口温度差的计算值作为最终计算结果;同时根据所述焓差和空气进出口温度差的假设值计算出每个单一计算区间空气的含湿量的差;A505、根据蒸发器结构计算模块和霜层厚度的初始设定值确定蒸发器的静压损失,同时利用风扇静压特性曲线计算此时流经蒸发器的风量,计算公式如下: <mrow> <mi>&Delta;</mi> <msub> <mi>p</mi> <mi>a</mi> </msub> <mo>=</mo> <mn>5.88</mn> <mo>&times;</mo> <msup> <mn>10</mn> <mrow> <mo>-</mo> <mn>4</mn> </mrow> </msup> <msub> <mi>N</mi> <mi>rd</mi> </msub> <mo>&times;</mo> <mfrac> <mrow> <mfrac> <mn>2</mn> <msub> <mi>p</mi> <mi>t</mi> </msub> </mfrac> <mo>[</mo> <msub> <mi>S</mi> <mn>2</mn> </msub> <mo>-</mo> <mfrac> <mrow> <mi>&pi;</mi> <msup> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mn>0</mn> </msub> <mo>+</mo> <mn>2</mn> <msub> <mi>&delta;</mi> <mi>fr</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> <msub> <mrow> <mn>4</mn> <mi>S</mi> </mrow> <mn>1</mn> </msub> </mfrac> <mo>]</mo> <mo>+</mo> <mfrac> <mrow> <mi>&pi;</mi> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mn>0</mn> </msub> <mo>+</mo> <mn>2</mn> <msub> <mi>&delta;</mi> <mi>fr</mi> </msub> <mo>)</mo> </mrow> </mrow> <msub> <mi>S</mi> <mn>1</mn> </msub> </mfrac> </mrow> <msup> <msub> <mi>S</mi> <mn>2</mn> </msub> <mn>0.3</mn> </msup> </mfrac> </mrow> <mrow> <mo>&times;</mo> <msup> <mrow> <mo>[</mo> <mfrac> <mrow> <msub> <mi>p</mi> <mi>t</mi> </msub> <mo>&CenterDot;</mo> <msub> <mi>S</mi> <mn>1</mn> </msub> </mrow> <mrow> <mo>[</mo> <msub> <mi>p</mi> <mi>t</mi> </msub> <mo>-</mo> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>f</mi> </msub> <mo>+</mo> <mn>2</mn> <msub> <mi>&delta;</mi> <mi>fr</mi> </msub> <mo>)</mo> </mrow> <mo>]</mo> <mo>[</mo> <msub> <mi>S</mi> <mn>1</mn> </msub> <mo>-</mo> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mn>0</mn> </msub> <mo>+</mo> <mn>2</mn> <msub> <mi>&delta;</mi> <mi>fr</mi> </msub> <mo>)</mo> </mrow> <mo>]</mo> </mrow> </mfrac> <mo>]</mo> </mrow> <mn>3</mn> </msup> <mo>&times;</mo> <msup> <msub> <mi>W</mi> <mi>F</mi> </msub> <mn>1.7</mn> </msup> </mrow>Δpa:空气侧静压损失,单位Pa;Nri:管排数;d0:管外径;δfr:霜层厚度,单位m;S1:管列间距,单位m;S2:管行间距,单位m;Pt:翅片间距,单位m;tf:翅片厚度,单位m;WF:迎面风速,单位m/s;调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算;A506、根据冷媒的假设的焓差,计算冷媒侧单一计算区间的换热量; <mrow> <msub> <mi>Q</mi> <mi>iref</mi> </msub> <mo>=</mo> <mover> <mi>m</mi> <mo>&CenterDot;</mo> </mover> <mrow> <mo>(</mo> <msub> <mi>hi</mi> <mi>out</mi> </msub> <mo>-</mo> <msub> <mi>hi</mi> <mi>in</mi> </msub> <mo>)</mo> </mrow> </mrow>Qiref:冷媒侧单一计算区间的换热量,单位kW;hiout:单一计算区间的出口焓值,单位kJ/kg;hiin:单一计算区间的入口焓值,单位kJ/kg;m:质量流量,单位kg/s;A507、进行冷媒侧局部换热系数计算;冷媒在蒸发器中流动,是一个蒸发过程,冷媒在蒸发器中具有三种不同的相态:冷媒在蒸发器中为液态时,所在的相区为液相区;冷媒在蒸发器中为气液混合物时,所在的相区为两相区;冷媒在蒸发器中为气态时,所在的相区为气相区;液相区与气相区由于都为单一相态,故称为单相区;在蒸发过程中,冷媒的变化为液体—气液混合物—气体的过程;在不同的相态下其计算方法完全不同;单相区的蒸发局部换热系数计算公式如下: <mrow> <msub> <mi>a</mi> <mi>ir</mi> </msub> <mo>=</mo> <mfrac> <msub> <mi>&lambda;</mi> <mi>i</mi> </msub> <msub> <mi>d</mi> <mi>i</mi> </msub> </mfrac> <mo>&CenterDot;</mo> <mn>0.023</mn> <mo>&CenterDot;</mo> <msup> <msub> <mi>Re</mi> <mi>g</mi> </msub> <mn>0.8</mn> </msup> <mo>&CenterDot;</mo> <msup> <msub> <mi>Pr</mi> <mi>g</mi> </msub> <mn>0.4</mn> </msup> </mrow>Reg:气相雷诺数;Prg:气相普朗特数;λi:液相导热系数;di:管内径;两相区的蒸发局部换热系数计算公式如下:air=Facv+Sapbacv为强制对流项;F:强制对流系数;apb为沸腾项;S:沸腾系数;调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算;A508、由导热方程计算管壁温度: <mrow> <msub> <mi>Q</mi> <mi>iref</mi> </msub> <mo>=</mo> <msub> <mi>a</mi> <mi>ir</mi> </msub> <msub> <mi>A</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>T</mi> <mi>iw</mi> </msub> <mo>-</mo> <mfrac> <mrow> <msub> <mi>T</mi> <mi>irin</mi> </msub> <mo>+</mo> <msub> <mi>T</mi> <mi>irout</mi> </msub> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow>Qiref:冷媒侧单一计算区间的换热量,单位kW;air:空气侧换热系数;Ai:管外换热面积;Tiw:单一计算区间内管壁温度,单位℃;Tirin:单一计算区间入口冷媒温度,单位℃;Tirout:单一计算区间出口冷媒温度,单位℃;调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算;A509、由结霜判断条件判断单一计算区间是否结霜;所述的结霜判断条件为同时满足以下两个条件时:计算所得的管壁温度低于0℃和单一计算区间的湿空气的露点温度高于管壁温度;如果单一计算区间满足所述的结霜判断条件,则为结霜工况,转入步骤A510进行计算;如果单一计算区间不满足所述的结霜判断条件,则为不结霜工况,确定霜层厚度变化值为0,空气侧含湿量差为0;转入步骤A514进行计算;A510、计算霜层表面温度的初值 <mrow> <msub> <mi>T</mi> <mi>ifr</mi> </msub> <mo>=</mo> <msub> <mi>T</mi> <mi>iw</mi> </msub> <mo>+</mo> <mfrac> <mrow> <mn>2</mn> <mo>&CenterDot;</mo> <msub> <mi>Q</mi> <mi>iref</mi> </msub> <mo>&CenterDot;</mo> <msub> <mi>m</mi> <mi>ir</mi> </msub> <mo>&CenterDot;</mo> <msub> <mi>i</mi> <mi>sv</mi> </msub> </mrow> <mrow> <mn>2</mn> <mo>&CenterDot;</mo> <msub> <mi>A</mi> <mi>i</mi> </msub> <mo>&CenterDot;</mo> <msub> <mi>l</mi> <mi>ifr</mi> </msub> </mrow> </mfrac> </mrow>Tifr:单一计算区间霜层表面温度,单位℃;Tiw:单一计算区间内管壁温度,单位℃;Qiref:冷媒侧单一计算区间的换热量,单位kW;mir:单一计算区间冷媒质量流量,单位kg/s;iir:冷媒汽化潜热,单位kJ/kg;Ai:管外换热面积,单位m2;lifr:单一计算区间的管长,单位m;式中isv为水和冰的相变潜热,AT为计算区域的换热面积,mr为导致霜层密度变化的质量,lifr为单一计算区间内的霜层的导热系数;A511、计算空气侧的换热量,所述的换热量包括显热换热量和潜热换热量;空气侧的显热换热量计算公式为Qisen=ma·Cpa(Tain‑Taout)Qisen:单一计算区间空气侧显热换热量,单位kW;ma:单一计算区间空气流量,单位kg/s;Cpa:单一计算区间空气比热,单位kJ/kg K;Tain:单一计算区间空气入口温度,单位℃;Taout:单一计算区间空气出口温度,单位℃;空气侧的潜热换热量计算公式为Qilat=ma·(din‑dout)Qilat:单一计算区间空气侧潜热换热量,单位kW;din:单一计算区间空气入口含湿量,单位kg;dout:单一计算区间空气出口含湿量,单位kg;调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算;A512、蒸发器单一计算区间换热量的迭代计算由于空气侧换热量与冷媒侧换热量必须相等,比较两个换热量大小,如果其差异小于计算收敛精度10‑4时,则认为步骤A504假设的焓差和单一计算区间空气进出口的温度差为准确值,这样就完成了单一计算区间换热量的计算,否则回到步骤A504重新调整所述焓差和空气进出口温度差的假设值进行计算,直到满足收敛迭代条件为止;A513、蒸发器整体计算区间计算完成了一个单一计算区间的计算后,在冷媒流动方向上,将上一个计算区间的出口冷媒的温度、焓差和压力损失作为下一个计算区间的入口条件;在空气流动方向上,将上一个计算区间的空气出口温度、含湿量作为下一个计算区间的入口条件;按照步骤A504‑A512的步骤进行计算;直到将步骤A502中划分的所有计算区间计算值全部求得;即得到整个蒸发器每一个计算区域的冷媒的焓差、压力损失、空气侧的温差、含湿量差以及霜层的密度和厚度的变化;A514、将上一个时间步长的计算的霜层厚度作为下一个时间步长的初始霜层厚度,由于每个单一区间的霜层厚度计算结果不同,将对不同的计算区间分别赋予相应的初始霜层厚度;A515、将计算所得的霜层厚度、霜层密度、蒸发器换热量与压力损失、空气侧温度差与含湿量差作为已知数据进行保存,供制冷循环计算模块调用;B、建立制冷循环系统计算模块B1、设置蒸发器的结霜周期、计算时间步长、初始的霜层密度和迭代收敛条件,所述的迭代收敛条件的默认值为10‑4;所述的结霜周期为蒸发器从开始工作至开始进行除霜的总时间,单位为h;计算时间步长:由用户将结霜周期平均划分成若干份,其中每一份时间规定为一个时间步长;同一时间步长中的蒸发温度、冷凝温度、冷媒流量、霜层厚度和霜层密度为稳态参数;所述的时间步长由用户将结霜周期进行平均划分,其中每一份时间规定为一个时间步长;同一时间步长中的蒸发温度、冷凝温度、冷媒流量、霜层厚度和霜层密度为稳态参数;所述的初始的霜层密度为开始结霜时刻霜层的密度大小,单位为kg/m3;所述的迭代收敛条件为假设参数值与实际计算参数值之间的差异;B2、设定运行参数,根据用户针对不同的工作条件设定蒸发器过热度、冷凝器过冷度、吸气压力损失和排气压力损失;同时设定室内侧空气温度与湿度、室外侧空气温度与湿度;蒸发器过热度:冷媒在蒸发器内吸收热量;冷媒从液态至两相,在两相状态下的温度为蒸发温度,在两相状态下冷媒温度不变;又从两相变为气相,温度升高;从蒸发温度至蒸发器出口时的气相温度差为蒸发器过热度;蒸发器过热度作为已知条件由用户输入确定,单位为℃;冷凝器过冷度:冷媒在冷凝器内吸收热量;冷媒从气态至两相,在两相状态下的温度为冷凝温度,在两相状态下冷媒温度不变;又从两相变为液相,温度下降;从冷凝温度至冷凝器出口时的液相温度差为冷凝器过冷度;冷凝器过冷度作为已知条件由用户输入确定,单位为℃;吸气压力损失:制冷设备从蒸发器至压缩机吸气口会有一段的管路连接,在管路内由于摩擦力作用,压力会有部分损失,导致冷媒压力下降,温度下降;另外,由于管路与外界的热交换也会产生一部分损失;这些损失统称为吸气压力损失;吸气压力损失用温度来表征,作为已知条件由用户输入确定,单位为℃;排气压力损失:制冷设备从压缩机排气口至冷凝器会有一段的管路连接,在管路内由于摩擦力作用,压力会有部分损失,导致冷媒压力下降,温度下降;另外,由于管路与外界的热交换也会产生一部分损失;这些损失统称为排气压 力损失;该排气压力损失用温度来表征,作为已知条件由用户输入确定,单位为℃;B3、首先假定蒸发温度及冷凝温度,计算出冷凝器热负荷,通过冷凝器热负荷确定冷凝器所需的换热面积A1;同时调用冷凝器结构计算模块以及由用户输入冷凝器侧空气入口干球温度、含湿量和冷凝器风扇风量确定实际冷凝器换热面积A2;冷凝器所需的换热面积A1与实际冷凝器换热面积A2进行迭代循环后确定冷凝器实际换热面积;再利用冷凝器热负荷、室外侧空气入口干球温度、含湿量和冷凝器风扇风量,确定冷凝器冷凝温度Tc;同时确定压缩机的制冷能力、蒸发器换热量Q2、输入功率以及电流; <mrow> <msub> <mi>A</mi> <mn>1</mn> </msub> <mo>=</mo> <mfrac> <msub> <mi>Q</mi> <mn>2</mn> </msub> <mrow> <msub> <mi>k</mi> <mn>0</mn> </msub> <mo>&CenterDot;</mo> <msub> <mi>&theta;</mi> <mi>m</mi> </msub> </mrow> </mfrac> </mrow>k0:冷凝器总换热系数,单位W/(m2·K);θm:对数平均温差; <mrow> <msub> <mi>&theta;</mi> <mi>m</mi> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>t</mi> <mrow> <mi>a</mi> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>t</mi> <mrow> <mi>a</mi> <mn>2</mn> </mrow> </msub> </mrow> <mrow> <mi>ln</mi> <mfrac> <mrow> <msub> <mi>t</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mi>t</mi> <mrow> <mi>a</mi> <mn>1</mn> </mrow> </msub> </mrow> <mrow> <msub> <mi>t</mi> <mi>k</mi> </msub> <mo>-</mo> <msub> <mi>t</mi> <mrow> <mi>a</mi> <mn>2</mn> </mrow> </msub> </mrow> </mfrac> </mrow> </mfrac> </mrow>ta1:空气入口干球温度;ta2:出口空气干球温度;tk:冷凝温度;根据冷凝器结构计算模块确定冷凝器实际换热面积A2;在冷凝器所需的换热面积A1与冷凝器实际换热面积A2不相等时,调整假设的冷凝温度使两者相等,此时的冷凝温度为所述的假设的蒸发温度下的真实冷凝温度;但由于蒸发温度目前仍然为假设值,故此冷凝温度尚不是循环系统中的真实冷凝温度;由结构确定的实际冷凝器换热面积A2:A2=(m‑1)πDtLCOREndnrDt:翅片管外径LCORE:单一计算区间的管程;nd:管段数;nr:管列数;m:面积扩大率m=Af(π·dci·Pf);dci:管内径;Pf:翅片间距;B4、根据确定的冷凝温度Tc、设定的冷凝器过冷度和假定的蒸发温度,确定蒸发器入口冷媒温度、冷媒压力、冷媒的流量以及蒸发器热负荷Q2;调用步骤A1的湿空气物性计算模块与步骤A2的冷媒物性计算模块进行计算;B5、调用蒸发器结霜计算模块进行结霜工况下的蒸发器仿真计算;B6、利用步骤B5得到的蒸发器在结霜工况下实际的蒸发器换热量Q1和步骤B3得到的蒸发器换热量Q2进行迭代计算,从而确定实际蒸发器侧换热量;同时确定蒸发器结霜厚度和结霜密度;最终确定实际的蒸发温度;完成一个时间步长内的计算结果,利用Visual C++数据传递工具,将计算结果输出到Excel表格中;B7、将上一个时间步长的计算结果作为下一个时间步长的计算初值,并更新整个计算区域;B8、当满足计算结束的判断条件时,计算结束;所述的计算结束的判断条件为:计算时间达到设置的结霜周期或计算的霜层厚度大于等于翅片间距的一半;计算结束后,将各个时间步长的平均结霜厚度、平均结霜密度、冷凝温度、蒸发温度、制冷量、蒸发器热负荷、冷凝器热负荷、能效比COP和冷媒循环量输入到Excel表格,并自动绘制图表。
地址 116600 辽宁省大连市经济技术开发区松岚街8号