卧式制冰机制冰的仿真研究

摘要：本文根据卧式海水制冰机的工作原理和结构建立了数学物理模型，对其制冰过程进行了数值模拟，并将计算结果与试验数据进行了对比分析，二者基本吻合，并研究了冰筒转速、海水温度和蒸发温度对成品冰的厚度产生的影响。

1 前言

当前片冰机的运用比较广泛，主要有立式和卧式两种。由于结构问题，立式片冰机存在布水不均匀、成品冰较湿、冰温较高和厚度薄等诟病。卧式片冰机则较好地解决了上述问题，其成品冰的质量好，冰温更低、冰片更厚，更适合要求较高的冷藏保鲜，并且它既适合淡水冰又适合海水冰的生产，应用更为广泛。国内外对这方面的研究不多。本文根据卧式海水制冰机的结构和工作原理对其结冰过程进行了热力学分析和数值模拟，并结合样机的试验数据，验证了数值模拟的可靠性，并讨论了影响其结冰性能的主要因素。

2 卧式制冰机结冰过程分析

卧式制冰机的结冰工作过程如图1所示：

图1 制冰机的工作原理

首先，低温的筒状蒸发器(冰筒)在水箱内低淋水区速旋转。在冰筒浸入水中的部分(如图中①处)，水直接在其表面冻结成冰，并且冰层的厚度随着冰筒的旋转逐渐增加。当结冰部分旋转到离开水箱内的水面之后，由喷淋机构的洒水使结冰处冰层继续增厚(图中②处)，直至转过喷淋点(图中③处)，冰层开始过冷，旋转到冰刀位置处，在冰刀的作用下剥离制冰筒表面，沿着导向板在重力的作用下滑入贮冰室。

如图1所示，按照冰层在制冰筒上的位置不同，将制冰筒的整个工作区分为以下三个：

1) 浸水区

在浸水区中，制冰筒与海水直接接触，其表面附近的海水逐渐被强制冷却降温，直至降到海水的凝固温度开始发生相变，液态的海水冻结成冰，并且冰层的厚度随着时间的推移而增加。浸水区是结冰的积累区，也是对制冰性能影响最大的一个区。同时海水结冰过程中排挤出来了盐份可以直接被冰表面流动的海水带走。

2) 淋水区

淋水区的主要功能是进一步增加冰层的厚度，需要注意的是必须要控制淋水的水温和流量。水的温度过高或流量过大都会导致冰层厚度的减少，同时也会带来过多的冷量损失，降低制冰效率。所以一般对于喷淋的水采取预冷处理。同时淋水机构的位置也对制冰性能有很大影响。

3) 过冷区

经过前面两个区的冻结工作区，冰层厚度已经达到了最大值，在过冷区内继续降温，直至在冰刀的作用下剥离制冰筒表面，得到最终的海水冰。

3 数学物理模型的建立

3.1 模型的简化

在进行数学建模之前，首先进行以下假设：

由于制冰筒直径远远大于冰层的厚度，所以在计算中忽略圆柱曲面的影响，将筒状制冰筒展开，认为是一个无限大平面状冷壁面;整个制冰过程为一维换热;制冰筒表面温度恒定且均匀，为T1;制冰过程中水箱中的水温恒定均匀，为T2;由于制冰筒转动较慢，在浸水区中忽略流动对其影响，没有对流换热，只考虑热传导，最后对水的传热量进行折算修正;结冰表面即发生相变处的温度恒为海水凝固温度Tm;冻结过程中，固态冰和液态水之间存在一个相变分界区域，并随着冻结时间的推移而移动。海水结冰的固—液分界面的移动速度对温度的扩散不产生影响;冻结前的海水和冻结后冰的密度、比热、导热系数等热物性参数恒定不变。

由于海水在冻结成冰的过程中，液态的海水凝固变成固态的海水冰发生相变，其对应的热物性参数也在相变前后发生很大的变化。为了方便分析，假设相变在一个平面上进行，这个平面定义为相变界面，通过相变界面把整个传热分析区划分为冻结区、非冻结区和相变区，如图2所示，其中相变界面是随着冻结时间的推移而移动的。

图2 从液态到固态的相变过程

3.2 导热微分方程

根据不同的状态区，建立如下不同的导热微分方程。

符号说明

Tk：固态冰的温度(℃);

Tq：液态水的温度(℃);

Tm：水的凝固温度(℃);

Cq，Ck：海水，海水冰的比热(J/kg.℃);

aq，ak：海水，海水冰的热扩散率(m/s2);

λq，λk：海水，海水冰的导热系数(w/m.℃);

ρq，ρk：海水，海水冰的密度(kg/m3);

H：水凝固潜热(J/kg);

W：冰筒转速(min/r);

S(τ)：相变界面的位移(m)。

1) 冻结区

在冻结区中，海水是以固态冰的形式存在的，假设与制冰筒接触的冰温等于制冰筒的蒸发温度T1，随着相变界面的推移冰层的厚度逐渐增加。在冻结区中的导热微分方程为：

a) ， （1）

此处的边界条件为当x=0时，温度恒为制冰筒的蒸发温度T1，其数学描述为：

b) （2）

2) 相变区

在相变区内，根据假设，相变过程只在相变界面上进行，并且固—液两相为单纯的导热方式，在相变界面上共有三部分热量参与热传导：

①.液相区内靠近相变界面处的海水由于降温而释放的热量q1：

（3）

②.固相区内靠近相变界面处的海水冰由于降温而释放的热量q2：

(4)

③.相变界面处液态海水凝固冻结成冰发生相变而释放的凝固潜热q3：

(5) 根据能量守恒，从固态一侧传出来的热量应等于海水的凝固潜热与从液态一侧传出来的热量和，即： （6） 整理之后的导热方程如下：

， （7） 此处的边界条件为固态冰和液态水在相变界面处的温度都为Tm，数学描述为：

（8） 3) 非冻结区

非冻结区中的导热微分方程和冻结区相似：

， （9） 此处的边界条件：

在x→∞处，液态海水的温度为T2，数学描述为：

（10）

4 数值模拟的方法

利用上面的导热微分方程组，可以进行求解计算来模拟不同参数对制冰性能的影响。由于此类问题属于相变传热问题，相变时伴随相变潜热的释放和吸收，同时相变界面还发生非线性的移动，导致数学模型的方程也呈现出非线性变化;在数值求解过程中，控制方程是导热方程，在相变界面处，连续性方程、能量方程、动量方程和质量守恒方程将发生非线性变化，它们之间需要相互耦合。结合卧式海水制冰机的工作原理和结构特点，本文采用有限差分法对各方程进行离散、求解[2]~[5]，并用Visual Foxpro语言编写出计算程序，对卧式制冰机的结冰过程进行了数值模拟。整个程序的流程图如图3所示。该程序的界面简单方便，只需按照界面提示输入各参数的值，然后点击“开始计算”按钮即可开始计算并得出结果。程序界面如图4所示。计算的结果直接生成“*.dbf”格式文件存入数据库中，可以方便地调出来查阅，并可根据数据绘制出曲线图。

图3 程序流程图

图4 程序操作界面

5 模拟计算与试验测试结果对比分析

本文将模拟出的数据和实际试验所得的数据进行对比、校验，并得出制冰机的各工作参数对成品冰的冰厚的影响.

1) 不同转速下的结冰厚度的对比结果

首先，本文模拟计算了海水制冰机在蒸发温度T1=-40℃，海水温度T2=25℃工况下，不同的冰筒转速下的成品冰的厚度，计算结果与试验结果对比关系曲线如图5所示。

图5 转速与冰厚对比关系

从和图5可知，计算与试验测试结果基本吻合，相对偏差在10%~15%范围之间。且可以得出其他条件不变的情况下，随着冰筒转速的增加，成品冰的厚度越来越小。

2) 不同蒸发温度下的冰厚对比结果

海水制冰机在制冰筒在海水温度T2=25℃，冰筒转速W = 2 min/r，不同的蒸发温度下的成品冰厚度的计算与试验结果对比关系曲线如图5所示。

图6 蒸发温度与冰厚对比关系

从图6可见，计算结果与试验测试结果基本吻合，偏差在3%~15%范围之间，并且可以得出其他条件不变的情况下，随着蒸发温度的增加，成品冰的厚度越来越大。

3) 不同海水温度下的冰厚对比结果

海水制冰机在蒸发温度T1=-40℃，冰筒转速W = 1 min/r，不同海水温度下的成品冰厚度的计算与测试结果对比关系曲线如图6所示。

图7 海水温度与冰厚对比关系

由图7 可知，计算结果与试验测试结果基本吻合，偏差在4%~15%，并可以得出其他条件不变的情况下，随着海水温度的增加，成品冰厚度越来越小。

上面计算结果与试验测试结果表明建立的数学模型和计算方法是可行的。尽管计算与试验存在误差，但是其变化趋势基本一致，偏差在可接受的范围内，为参数的合理性选择提供了理论依据。

6 结论

本文根据海水制冰机的基本工作原理建立了数学物理模型，以微元体作为研究对象，把结冰过程划分为冻结区、相变区和非冻结区三个区域，针对不同区域建立不同的导热微分方程;通过计算机编程求解，并对计算的结果进行了分析;最后将计算结果与试验测试的结果进行对比分析，二者基本吻合，说明了这里建立的数学模型和计算方法是可行的。并且可以得出随着冰筒转速和海水温度的增加、蒸发温度的减少，成品冰的厚度越来越小。这一数学物理模型为卧式海水制冰机的实际设计和性能优化提供了必要的理论依据。

参考文献

[1] 陈天及，李晓辉。 制冰技术研究的新进展。冷藏技术，2005 3，(1)：1～4

[2] 崔小朝等。伴有相变的热传导焓式有限元法。太原重型机械学院学报,16(2):169～173

[3] Koji Kaminishi, Tetsuya Araki, Ryo Shirakashi. Three-Phase model for simulating freezing process of aqueous solutions. International Congress of Refrigeration 2003,Washington,D.C

[4] 任松保，郑正泉。冻结过程中相变界面移动及传热的计算机模拟。华中科技大学学报(自然科学版)，31(2)：61～63

[5] 曾欣,辛明道.耦合求解焓及移动边界位置的多维相变传热问题的数值焓法. 重庆大学学报,15(3):62～67