1 引言

　　冰蓄冷空调系统是在空调负荷很低时制冷储存，而在新空调负荷高峰时化冰取冷，因而可以全部或部分转移制冷设备的运行时间，从而缩小供电网络的峰从电力负荷差。因此，冰蓄冷空调系统在美国、日本、澳大利亚等国家发展非常迅速，已广泛使用十几年之久。
　　冰蓄冷空调系统由蓄冰设备和空调系统组成。根据ARI《蓄冷设备热性能指南》，蓄冰设备可分为冰盘管、冰球、冰片滑落、冰晶等几种形式^[1]。其中应用最广泛的是冰盘管和冰球 ^[2] 。
　　近年来许多学者对冰盘管式蓄冰设备进行了研究^[3]，生产厂家也开发出许多系列的冰盘管产品。但由于冰盘管管道较细，流动阻力较大，技术要求较高，目前国内还没有一家生产冰盘管的厂家。而冰球式蓄冰设备由于结构简单，已成为蓄冰系统的一个重要发展方向^[2]。
　　由于密封球体内的结冰和融冰过程是一个伴随着相变的导热与自然对流换热的复杂过程，掠过冰球的载冷剂的流动和换热过程也很复杂，因此目前很难看到有关冰球蓄冷和取冷方面的文章和厂家资料。国内较早研究冰球特性的是清华大学的赵庆珠教授，她在1994年发表两篇关于单个冰球蓄冰和融冰过程的文章^[4，5]。而实际的冰球蓄冰设备是将许多冰球堆集在一个圆形断面蓄冰罐中，因此了解蓄冰罐的蓄冰和融冰特性对冰球式蓄冰罐的生产、设计和使用具有重要的意义。

2 数学模型的建立和求解

　　载冷剂在冰球式蓄冰罐中的流动是一个很复杂的过程，涉及载冷剂掠过冰球的换热以及冰球内的要变换热过程。对此作如下简化：①冰球在罐内均匀分布，载冷剂和冰球中介质的热工参数只沿流程发生变化；②冰球的传热性能只与冰球中的冰量有多少有关^[4，5]。根据这两个简化条件即可得到冰球式蓄冰罐的数学模型。

　　2．1 数学模型
　　将冰罐中的冰球沿载冷剂流程分成若干排，每排中的冰球个数为N，每个冰球的表面积为A，每排冰球占据的长度为E，冰罐中载冷剂的有效流通面积为F。在每一排中，载冷剂的温度相同，冰球的温度相同。载冷剂侧的能量方程为

　　　　 　　　　　　　　　　　（1）

　　 　　　　　　　　　　　　　　　　　　　　（2）

　　　　 　　　　　　　　　　　　　　　　　　（3）

　　　　　　　　　　　　　　　　　　　　　　（4）

　　式（1）～（4）中，T_f为载冷剂温度，τ为时间，u _f为载冷剂流速，a _f为载冷剂的导温系数，q _ht为单位体积载冷剂与冰球间的换热量，q_fls为单位体积载冷剂通过罐壁的热损失，ρ_f，c _f分别为载冷剂的密度和比热容，K _c为冰球的温度，D为冰罐直径，T_amb为周围环境温度，R_w为冰罐壁热阻，a_w为载冷剂与冰球外表面的换热系数，a_n为冰球中的介质与冰球内表面的换热系数，R₀为冰球壳热阻。
　　其中a_w由下式确定[6]：

　　 　　　　　　　　　　　　（5）

　　 　　　　　　　　　　　　　　　　　　（6）

　　式（5）（6）中，d_p为冰球直径，λ_f为载冷剂的导热系数，Pr为载冷剂的普朗特数，Re_p为载冷剂的雷诺数，v₀为罐体中载冷剂的表现速度，ν_f为载冷剂的运动粘度。a_n由文献[7]回归得到。冰球侧的能量方程为

　　　　 　　　　　　　　　　　　（7）

　　冰球中无冰时： 　　　　　　　　　（8）
　　冰球中冰水共存时： 　　　　　　　 （9）
　　冰球中全冰时： 　　　　　　　　　 （10）
　　初如条件：T_b =T₀ ，T_f =T₀，ε=ε₀　　　　　　　　　(11)

　　边界条件：T_f | _x=0 =T_fin　　　　　　　　　　　　　 (12)

　　式（7）～（12）中，Q_ht为单位冰球与载冷剂间的传热量，m为单个冰球中介质的质量，c_w 为水的比热容，ε为冰球中的结冰率，q_lt为冰的常融解潜热，c_i为冰的比热容，T₀为初初始温度，ε₀为初始结冰率，T_fin为冰罐中载冷剂的入口温度。

　　2．2 方程的离散和求解
　　采用隐式格式对上述方程进行离散，用追赶法求解离散后的方程组。具体计算过程如图1所示。
　　
　　　　　　　　　　　　　　　　　　　
　　　　　　　　　　　　　　　　　　　　　　　　图１ 框图STAIT程序的计算
　　
　　根据图1所示的计算框图，笔者编制了计算冰球式蓄冰罐的计算程序STAIT(Software of Time-dependent Analysis for Ice Tube).