(上海理工大学生物系统热科学研究所 上海 200093)

低温保存是细胞和组织等生物材料长期保存的重要手段，在低温生物学、转化医学等众多学科研究及临床治疗中都具有重要意义[1]。在降温过程中，每种生物体均对应一个最佳冷却速率，若控制不当，极易造成细胞损伤[2]。目前采用的低温保存方法主要分为慢速冷冻、玻璃化和超快速冷冻3种[3-6]。

慢速冷冻可通过液氮程序降温仪实现，通过控制液氮注入量或气相区高度来达到设定降温程序。这种方法使用低浓度低温保护剂，通常降温时间很长，需消耗大量液氮，且仪器昂贵，并会因液氮的使用而受限[7-9]。玻璃化通常使用较高浓度低温保护剂及较快的降温速率来避免冰晶的产生以提高成活率。但其所需的步骤比慢速冷冻复杂，且只能用于小体积样本[10-13]。超快速冷冻是介于慢速冷冻和玻璃化之间的技术，但其应用表现低于前述两种方法。根据阿尔赫尼厄斯(S.Arrhenius)公式，生化反应速率k和绝对温度T之间存在如下关系：

(1)

式中：k为反应速度常数；R为摩尔气体常数，J/(mol·℃)；T为绝对温度，K；Ea为活化能，J/mol；A为Arrhenius因子。对于给定的反应，A为常数。由式(1)可知，温度越低其化学反应的速率越低，保存时间越长。采用上述方法冷冻的样本，由于成本和技术问题，最低只能在液氮温度-196 ℃下长期保存。此外，一些研究表明液氮并非无菌，可能是潜在的污染源。从程序降温装置排出的气体也不是无菌的，使得这些装置不能应用于无菌工作环境中[14-17]。

斯特林制冷机是一种闭式循环的气体机械制冷机，工作温区广，制冷温度可达120 K以下，并且结构紧凑、效率高[18]。将其用于样本速冻保存一定程度上可实现液氮制冷的效果，实现更低的储存温度，装置便携性好且适用于无菌工作环境[19-20]。斯特林制冷机制得的冷量先经冷头传导给冷端适配器，再由导冷板完成对外冷量传导，实现制冷。因此，冷量的快速有效传导是制冷机实现高效制冷的关键[21-23]。由于应用场合和样本实际保存的需求，急需研制出一种高效便携的速冻装置。本文基于斯特林制冷机设计了一种便携式样本速冻装置，并对冷台的形状进行了模拟优化，得到冷台结构设计中对样本降温影响较大的因素，提出了优化方案，为今后冷冻样本的装置设计与实验提供参考。

1 样本速冻平台设计及计算模型

1.1 几何模型

速冻平台几何模型如图1所示。本次模拟选用宁波华斯特林生产的Cryo-s系列80 W斯特林制冷机为冷源，模型主要由导冷平板、冷台、冻存管和空气域组成。选用1.8 mL冻存管，以水为冻存样本。

图1 速冻平台几何模型Fig.1 Geometry model of quick freezing platform

为便于观察平台内部结构，空气域只画出了1/2，如图1(a)所示。各部分材料参数如表1所示。

表1 速冻平台各部分材料参数Tab.1 The material parameters of quick freezing platform

由于速冻平台整体结构对称，可选用1/4体积进行模拟，此外选用两个样本为对照点，如图2所示。

图2 模拟对照点示意图Fig.2 Schematic diagram of the simulation control points

冷冻过程中，制冷机冷头通过导冷平板将冷量传递给冷台及外围空气，然后冷台及空气分别通过热传导及对流与冻存管进行换热，最终将冷量传递给样本。由于设计导冷平板选用紫铜材料，且厚度仅为2 mm，冷量在平板上的传递较为迅速。冷台自身体积、接触热阻及样本相变潜热的释放是影响冷量传递的主要因素。因此，本文主要针对冻存管间距、冷台高度和冷台半径进行模拟分析，具体参数如表2所示。

表2 模拟扫描参数Tab.2 The parameter for simulation

1.2 数学模型

为方便研究，本文对模拟所需数学模型进行了如下假设[24-25]：1)域内空气为牛顿流体；2)域内空气流动形态为稳态、层流；3)域内空气在固体表面满足无滑移边界条件；4)满足Boussinesq假设[26]；5)空气域外部、导冷平板下部绝热；6)模型无漏热情况发生。

1.3 物理模型及数值方法

1)网格剖分

利用COMSOL Multiphysics 4.4软件对模型进行网格划分。采用结构化网格对冻存管单独扫掠操作，样本域固定单元数为20，管帽和管内空气域固定单元数为5，其余部件采用自由四面体的剖分方法。网格大小按普通物理场超细化处理，最大单元尺寸均不超过1 cm。为保证计算精确度，对冻存管管壁网格局部加密[27]；为防止网格数量对计算结果产生影响，进行网格无关性验证[28]。由于模型的剖分方法相同，因此只需选择一种模型验证即可。选取冷台高度H=30 mm，冷台半径R=22.75 mm，冻存管间距X=15 mm的模型，网格数量分别为106 983，158 031，235 260及351 301。计算得到的样本1平均温度在158 031与106 983间差别为5%，在235 260与158 031间差别为2%，在351 301与235 260间差别为0.4%，即当网格数量达到235 260时，继续增加网格时结果变化不大。而235 260为按照上述方法对模型进行的网格剖分。图3所示为模型网格质量图，质量指标从0～1，趋于0说明网格质量较差，模型平均质量为0.82。综合考虑采用上述扫掠加自由四面体网格剖分模型进行计算。

图3 速冻平台网格质量统计Fig.3 Mesh quality statistics of quick freezing platform

2)边界条件

(1)空气域。冷台外部空气初始温度为25 ℃(298.15 K)，设定空气外壁面绝热，不与外界换热。

(2)冻存管。将管壁设置为薄层边界，管内样本(水)设置为相变材料，相变起始温度0 ℃(273.15 K)，相变区间10 K。

(3)导冷部件。冷台和导冷平板初始温度为25 ℃(298.15 K)，与空气接触部分设置对流传热边界。导冷平板底面绝热。

(4)边界热源。参照Cryo-s 80 W产品参数，计算得出制冷机制冷量Q与冷头温度T关系式为：

Q=0.002 4T2-0.598 3T+40.611

(2)

因此设置热耗率：

(3)

式中：A为制冷机冷头与导冷平板接触面积，m2。

3)数值模拟方法

利用COMSOL Multiphysics 4.4求解器对不同几何参数速冻平台进行数值模拟，模拟冷冻时间300 s，添加传热和层流多物理场耦合。由于为瞬态非等温流模型，瞬态求解器采用计算精度及稳定性均较高的直接式MUMPS求解器。

2 结果分析

2.1 管间距对降温速率的影响

图4 不同管间距下降温300 s模拟结果Fig.4 The simulation results of temperature distribution at 300 s for different tube spacing

模拟选用冷台高度H=30 mm，冷台半径R=27.75 mm，对两组样本间距分别进行模拟，模拟结果如图4和图5所示。由图可知两组外围管样本在降温300 s时均处于相变期，根据图5可知样本1和样本2 在两种管间距下降温速率差别较小，说明样本间距对降温速率影响不显著。此外还可知温度沿导冷平板横向变化较大，可能由于平板厚度较小，无法及时将冷量传递给外围冷台，导致内外冷台温差较大[29]。后期可考虑对导冷平板进行结构优化，底面采用弧形，降低内外冷台传热温差。

图5 不同管间距下样本降温曲线Fig.5 The cooling curve of samples at different tube spacing

2.2 冷台高度对样本降温速率的影响

图6 不同冷台高度下降温300 s模拟结果Fig.6 The simulation results of temperature distribution at 300 s for different cooling stage heights

图7 不同冷台高度下样本降温曲线Fig.7 The cooling curve of samples at different cooling stage heights

图8 不同冷台半径下降温300 s模拟结果Fig.8 The simulation results of temperature distribution at 300 s for different cooling stage radius

图9 不同冷台半径下样本降温曲线Fig.9 The cooling curve of samples at different cooling stage radius

由于管间距对样本的降温速率影响不显著，因此选用样本间距X=15 mm，此外由于冷台高度对样本传热的影响主要在y轴方向，因此冷台半径选用最小的R=22.75 mm，对4组不同冷台高度分别进行模拟。图6所示为模拟结果温度云图，由图6可知同一冷台内的样本温度相差较小，因此每个冷台选取一个样本，对比不同高度冷台的降温曲线如图7所示。可以看出，降温300 s时外围管与中心管温差较大，与冷台接触的样本温度较低，而暴露于空气中的部分温度较高，特别是中心管，H=20 mm时样本两部分温差超过30 K，说明固体接触热传导在样本降温中占主要作用[30]。由图7(a)可知，在降温前期，H=20 mm的冷台降温最快，但后期逐渐变慢。原因是冷台体积较小，前期传热快，但后期样本相变时释放大量潜热，除去制冷机冷量输入外，由于冷台体积受限未能蓄积太多冷量抵消潜热，导致整体降温变慢。H=20 mm的冷台样本上下部分温差过大，样本整体降温效果不佳。其余几组中心管H=35 mm温度最高，外围管样本H=30、35 mm的冷台样本刚开始相变时，H=25 mm的冷台样本冻结已基本结束，说明其外围管的降温最快。

2.3 不同冷台半径对样本降温速率影响

模拟选用冷台高度H=25 mm，样本间距X=15 mm，对3组冷台半径分别进行模拟，结果如图8所示。由图8可知，中心管和外围管R=22.75 mm冷台的样本温度均为最低，随着半径增大，样本的温度逐渐升高。说明冷台半径对样本降温影响较大。此外，观察外围管可知R=32.75 mm的样本刚开始相变时，其余两组已经快结束。图9所示为不同冷台半径下样本降温曲线。由图9可知R=32.75 mm的样本温差更小，原因是冷台体积较大，导致温度波动性较小[31]。后续设计对温度波动较敏感的样本冷台时，可考虑加大冷台半径。

3 结论

对于斯特林制冷机应用研究而言，如何保证机器冷量高效率输出一直是研究的重点[32]，而冷量的高效传导是其中的关键之一。本文基于斯特林制冷机设计了样本速冻平台，得到如下结论：

1)管间距对样本降温无显著影响，后续冷台设计时在不影响样本取放的基础上，可尽量减小管间距以提升冷台利用率及冷量传输效率。此外，制冷机冷头导冷部件对冷量的快速传递至关重要，后期可对导冷部件结构进行优化。

2)冷台高度对样本冷冻过程中降温速率及温度均匀性影响较大。冷台过低时，由于冷冻槽只能接触部分样本，使样本的降温出现明显“断层现象”，且平均降温速率较小，温差很大。当冷台高度超过管内样本高度时，温度均匀性有所提高，但降温速率下降。模拟4组样本效果最好的是H=25 mm，略低于样本高度。

3)冷台半径越小，样本降温速率越快。半径越大管内温差相对越小，但不同半径之间差距较小。后续如对降温速率要求较高的样本可采取缩小冷台半径的方法，而对需要超低温暂存或冷冻转运的样本相应冷台半径可适量增加。

实际应用中为减少冷量耗散，箱体保温的设计也至关重要，如采用导热系数较小的真空绝热板，并提升生产工艺减少“边缘效应”[33]。此外，强化制冷机散热、单批样本冻存数量多的装置采用多制冷机供冷、选用无运动部件的低温制冷机解决冷头震动问题[34]等，这些在今后研究中需要重点关注。