固体表面玻璃化冷冻方法降温结晶特性数值研究

朱文欣 赵芷慧 黄永华 程锦生 陈 威 李 铮

(1 上海交通大学机械与动力工程学院 上海 200240;2 上海东富龙医疗装备有限公司 上海 201108;3 上海交通大学附属第一人民医院泌尿中心男科 上海 200080)

细胞疗法被认为是未来医学的三大支柱之一,生物样本将成为重要战略资源[1]。由于生物的活性在低温环境中受到抑制,代谢等化学反应减慢近乎停滞[2-3],使细胞等生物长期冻存成为可能[4]。然而,冷冻过程温跨大,造成细胞等生物在冷冻过程中产生不可逆损伤,如氧化性损伤、渗透性损伤和机械损伤等[5-6]。

为减少冷冻损伤,多种冷冻方法在细胞等生物样本冻存领域获得应用[7],如直接接触液氮冷却法[8]、环网状冷冻法[9]、微流控冷冻法[10]等。然而上述方法均易造成样品受热不均、回收困难并存在样品污染的风险。而固体表面玻璃化冷冻法(solid surface vitrification,SSV)由于降温速率快、降温速率控制性好、样品不直接与液氮接触、易回收、无污染等优点,被越来越多地应用于低温生物冻存[11-12]。目前,针对SSV方法的研究多数为实验或经验性报道[13-14],主要聚焦于降温速率的大小,对于SSV方法在冻存过程中造成的冷冻样品结晶机理与结晶度变化规律鲜有涉及。此外,由于冻存过程时间短,通过实验手段以高时间分辨率获得全冷冻过程样品全域温度、结晶度变化仍较为困难,如采用热电偶测温不仅侵入性大、易结露、漏热严重,而且无法测量结晶度;低温显微系统价格昂贵且无法定量分析冷冻过程的结晶变化[15]等。

鉴于此,本文采用数值计算方法探究SSV方法实现快速降温的核心要素,并通过建立冻存过程的非等温结晶预测模型,研究预冷温度、样品体积及形状对冷冻结晶效果的影响,为冷冻方法优化、冷冻对象体积形状等设计提供理论指导。

1 基于非等温结晶模型的SSV方法建模

SSV方法采用一定质量的高导热金属(铜、银或铝合金)块或片浸入液氮或其他低温介质预冷,将样品滴在预冷后的低温表面,通过瞬间接触导热实现快速降温冷冻,物理过程如图1所示。

图1 SSV方法工作原理

聚焦液滴碰触金属表面之后的传热过程,建立如图2所示的数学物理模型。鉴于轴对称性,该模型采用二维轴对称形式,取一半,其中左上角凸出部分为液滴,其余长方形区域为金属基底载冷台。以5 μL液体冷冻为例,真实滴落液滴的几何尺寸可以近似为半径1.82 mm,中心高度为0.89 mm的球冠,液滴的底部与金属块的上表面直接接触。此处考虑无氧铜材质的基底载冷台,其导热系数在低温下达400 W/(m·K)。为了消除边缘效应,取铜块的宽度为50 mm,高度为20 mm,此时相对液滴而言,可以认为铜块无限大。

图2 数学物理模型和边界条件

基于傅里叶导热定律,金属块内的导热可表示为:

(1)

式中:T为温度,K;ρ为密度,kg/m3;λ为导热系数,W/(m·K);cp为比定压热容,J/(kg·K);t为时间,s。

(2)

(3)

式中:L为熔融潜热,kJ/mol;X为结晶度。忽略结晶界面处由于分子扩散对冰晶前沿生长速率的影响,根据半经验非等温结晶模型[16-20]:

(4)

式中:Tm为熔融温度,K;k1为结晶常数,与液体种类、质量浓度相关;Q为活化能,kJ/mol;R为气体常数8.314 J/(mol·K)。

由于实际采用的抛光铜表面粗糙度Ra为0.8,视作光滑表面条件,且该铜表面与液滴直接接触,因此在进行建模时忽略了金属块与液滴样品间的接触热阻:

TS=TC

(5)

(6)

在边界条件方面,将金属块与液氮接触边界设为第一类边界条件,由于金属块与液氮蒸气对流换热,因此其他边界设为第三类边界条件,对流换热表面传热系数h取为100 W/(m2·K)[21]。

在热物性方面,铜块选用常物性参数,如表1所示。

表1 固体材料热物性参数

需要说明的是,该模型侧重于了解被冷冻液滴所经历的传热过程,并实现与系统的传质过程解耦,其分析尺度大于冰晶生长的微结构,且不包含溶液传质浓度变化对冰晶微结构的影响。

在COMSOL软件中,模型划分自由三角形网格,在样品与预冷表面接触处做细化处理,共划分了1 500、2 500、3 035、4 000四种密度的网格,网格平均质量均优于0.93。经网格无关性验证,取3 035个网格进行计算在速度和精度方面均满足要求。对金属块用固体传热接口求解,对CPA液滴利用导热接口添加广义热源项,通过域常微分和微分代数方程接口定义非等温结晶源项,并耦合求解。添加温度判据以保证结晶在熔融温度下才发生。为保证对极其微小量数值计算正常,将结晶度初始值设为1×10-16(可识别的接近0的最小正值)。采用瞬态求解器,向后差分方法,时间步长不超过0.01 s,相对容差控制为1×10-6。

采用文献[22]的实测数据来验证模型的适用性,针对西林瓶(类薄层式SSV结构)中1.9 mL甘油水溶液置于液氮中冷冻过程温度变化进行对比验证。需要说明的是,验证时模型中固体材料选用玻璃,溶液选用甘油水溶液,与文献[22]的实验保持一致。预测结果和文献数据的对比如图3所示。由图3可知,耦合非等温结晶源项的传热区域模型可以较好地预测冷冻过程,最大误差约为10%,平均误差小于2%。后期温度略高于实验值,可能是由于被测样品溶液的导热系数在低温时有所升高。

图3 模型计算值与文献[22]实验值对比

2 计算结果与分析

模型验证采用的实验样本液滴为生物细胞冷冻操作常用的冷冻保护剂(CPA)丙二醇溶液(质量分数为35%),其玻璃化转变温度Tg为160 K,冷冻结晶的危险温区介于Tg与熔化温度Tm之间(160～255.5 K)。由于CPA物性参数通常缺乏低温区的实测数据,采用常温区参数一般会低估物性参数的变化;而若采用水和冰的物性参数替代低温保护剂物性参数,会高估低温区物性参数的变化。在没有准确合适的物性参数数据时,有必要假设物性参数不随温度变化来参考计算[23]。一般而言,随着温度的降低,CPA导热系数增加,而比热容和密度降低,总体上热扩散系数增加。因此利用常物性参数计算可能会低估冷却速率,即实际降温效果优于计算值。样品液滴(丙二醇溶液)物性参数来源于文献[24-25],如表2所示。利用上述模型和物性参数,开展被冷冻样品滴落后的冷却和结晶特性研究。

表2 冷冻液滴热物性参数

2.1 非等温结晶特性

首先从液滴对称轴的高、中、低三个位点进行局部降温结晶特性分析,如图4所示。三个位点的冷冻降温过程均未呈现温度平台期,无明显的一级相变特征。由于低点最接近预冷表面,传热路径短、热阻小,因此液滴温度在低点处最接近预冷表面温度,降温速率最快,相比于液滴中点及高点的降温速率有数量级优势。不同的降温特性导致液滴起始结晶时间和最终结晶度不同。根据Uhlmann理论,结晶度小于10-6即为玻璃化[26],低点的结晶度几乎不变,趋向初始设定值10-16量级,实现了玻璃化转变。中点降温速度比高点更快,先于高点降温到熔融温度,所以结晶比高点更早产生。同时,由于中点比高点的传热路径更小,传热速率更大,能较快通过结晶危险温区,最终结晶度小于高点。只有以较高的温变速率快速通过结晶危险温区,才能实现被冷冻样品的玻璃化转变。

图4 样品内低、中、高点处温度与结晶度变化规律

对液滴样品选取切面进行整体降温结晶特性分析,如图5所示(注:图中只针对实现玻璃化冷冻及处于危险温区的部分显示云图)。液滴的降温最先从底部发生并实现玻璃化转变。样品中间结晶度增大,并向上扩散结晶,由于液滴边缘薄中间厚,结晶度与等温线均呈上凹形。但最大结晶度未出现在距离冷源最远、温度最高的液滴顶端,而出现于距顶端一定距离的中部。需要说明的是,由于本文结晶模型与传质解耦,因此,液滴中部存在最大结晶度这一现象与冷冻过程中样品界面处存在液气-固气相变无关,今后可进一步分析液滴内部复杂的温度变化情况对样品结晶度的影响。

图5 液滴内部结晶度和温度云图

进一步对冷冻过程中的液滴顶点处与结晶度最大点处开展降温速率和结晶度变化速率特性分析,如图6所示。顶点的降温速率和最大结晶处的降温速率趋势一致,但被冷冻液滴内部的降温速度存在不均匀性,不利于大体积样品的均匀冷冻。最大结晶度变化速率点比最大降温速率点滞后约0.3 s,表明结晶过程存在滞后性。在整个液滴均出现结晶度后,顶点的降温速率始终大于最大结晶点处的降温速率,顶点的液体更快通过结晶危险温区,表明更大的降温速率可以抑制结晶,使被冷冻样品更容易实现玻璃化。

图6 液滴顶点、结晶度最大处的降温速率与结晶度变化速率

需要说明的是,上文通过引入冷冻样品的结晶度X,更直观地展示了样品在被冷冻过程中的内部状态,可预测样品的被冷冻效果,为设法实现无损冻存的各种尝试提供便捷高效的验证方法。

2.2 固体表面玻璃化冷冻影响因素

针对不同的生物样本,通常需要设定不同的冷却预冷温度。因此,有必要探索预冷温度对样品在固体表面被玻璃化冷冻的影响。液氮是冷却金属载物台最常见的工质,其预冷极限温度Tprec,min为78 K,而35%丙二醇溶液的Tm为255.5 K,Tg为160 K。因此,分别选取Tprec为80、100、120、140、160、200、250 K。图7所示为预冷温度对液滴结晶度和穿越危险温区时降温最慢处的平均降温速率vcr的影响。预冷温度越低,降温速率越慢,结晶度越小,预冷温度与降温速率和结晶度均无线性关系。此外,当预冷温度低于样品玻璃化转变温度时,样品的降温速率均较快,达到103K/min量级,其结晶度在10-5量级;而当预冷温度高于Tg时,降温速率在102K/min量级,结晶度在10-5～10-4量级;当预冷温度越接近熔融温度,降温速率越慢,极易结晶。因此,预冷温度越低,传热温差越大,越有利于实现样品的玻璃化冷冻。

前文非等温结晶模型计算表明,最大结晶处可能在被冷冻样品的中上部。因此,有必要探究样品体积对用SSV方法冷冻的样本的降温速率和结晶度的影响。基于现有实际应用和文献调研情况,样品体积V分别选取1、5、10、20、50、100、500 μL,定义玻璃化比例VR为样品中实现玻璃化冷冻体积与样品总体积之比,计算结果如图8所示。当体积为1 μL时,样品完全玻璃化冷冻,因此,完全玻璃化转变的安全样品体积为1 μL。随着样本体积增加,样本内部导热路径变长,降温速率变慢,结晶显著增加,VR不断下降,且呈非线性关系。在一定范围内,体积越大,结晶度变化越剧烈;而体积过大,结晶度将完全结晶。因此,样本体积是影响SSV方法的关键因素之一。本文提出的耦合热源项的非等温结晶模型可以预测不同体积的样品内部结晶程度和不同条件下玻璃化安全体积,为实际固体表面玻璃化冷冻设计提供一定依据。

图8 不同体积样品中结晶度和玻璃化转变率

此外,根据傅里叶导热定律和牛顿冷却公式,有效传热面积是研究传热问题的关键因素之一。因此,引入传热面积与总体积比值A/V参数,研究样品形状对SSV方法冷冻的影响。样品体积为5 μL,改变样品形状,如图9所示。正球缺即为图 2介绍的模型;反球缺与正球缺尺寸相同倒置在挖空的金属预冷表面;反半球形样品半径为1.34 mm;反圆柱形半径为1.51 mm,厚0.7 mm;反薄层形半径为2.3 mm,厚0.3 mm。

图9 样品形状

不同形状液滴结晶度和玻璃化比例的计算结果如图10所示。由图10可知,相同体积下薄层形样品的A/V更大,玻璃化程度更高,而其他形状样品出现部分结晶。因此,A/V一定程度上可以反映传热拓扑结构的优劣。相同情况下,A/V值越大,结晶度越小,样品越容易实现玻璃化冷冻。若以实现玻璃化冷冻为目标,样品形状应尽量薄层化。说明本文模型可以简便高效地研究不同形状样品在冷冻时内部结晶情况,可为实际冷冻样品形状设计提供依据。

图10 不同形状液滴结晶度和玻璃化比例

3 结论

本文对基于固体表面玻璃化冷冻方法的冷冻液滴建立了耦合非等温结晶源项的传热区域模型,研究了液滴降温特性和冷冻液滴内部的结晶情况,得到如下结论:

1)固体表面玻璃化冷冻方法可以实现快速降温,降温速率可达到103K/min量级。其可获得快速降温能力的核心要素是固体冷台得到了充分预冷并储存足够冷量。

2)降温速率越快,结晶程度越小。结晶易出现在远离预冷表面的样品上部,结晶最严重区域不一定位于样品顶端。

3)预冷温度越低、样品体积越小、样品面积/体积(即薄层化)越大,结晶程度越小,越容易实现玻璃化冷冻。

4)耦合非等温结晶源项的传热区域模型可以通过结晶度X反映冷冻过程中样品内部结构变化,可作为设计开发新型冷冻设备和低温保护剂的辅助手段。

本文受上海市科委创新行动计划(20S31903400);上海交通大学“交大之星”计划医工交叉研究基金(YG2021QN86)项目资助。(The project was supported by Innovation Action Plan of Shanghai Science and Technology Commission (No. 20S31903400), and Shanghai Jiao Tong University Med-Engineering Fund (No. YG2021QN86).)