多孔介质冻结过程中温度变化的数值模拟

    多孔介质冻结过程的热质传递问题涉及很多领域,是构成众多自然现象的基本过程。在生物以及工农业生产中有着重要的应用前景,如人体器官保存、低温外科手术、冶金、晶体的成长、食品速冻保鲜、地下建筑等。但是由于多孔介质冻结过程中的固液相变属于多相体系、非连续介质的复杂热物理过程,表现为当在固相和液相的界面处吸收或释性;相变传热涉及多孔介质中已冻结的固相、未冻结的液相和固体骨架以及冻结模糊区(mushyzone)[1-2]。对多孔介质冻结的研究还没有形成一系列完善的理论,多数研究都只是借助于计算机技术建立模型给定数值解。实际上,热质传递是多孔介质冻结过程中重要的基本现象,笔者在借鉴现有数值模型的基础上,考虑多孔介质的内部特点,采用容积平均法建立数学模型,对无限大平板状多孔介质在第三类边界条件下发生的一维冻结过程中温度分布进行分析并加以实验验证。

    1 物理模型及数学描述

    假设将厚度为2δ,初始温度为T0(大于相变温度Tf)的无限大多孔介质平板置于温度低于Tf的对流环境中(对流换热系数为hc,环境温度为T∞)。对此类对称性边界条件可以简化为δ的平板一侧置于温度为T∞的对流环境中,一侧为绝热边界条件。如图1所示。

    系统初始温度为T0,含水量为W0。热量从中心(X=0)向表面(X=δ)传递。冻结从表面开始然后向中心移动。

    为获得数学方程还必须作以下假定:

    ①多孔介质的骨架、固相和液相的热物性是均匀的、各向同性的;

    ②在相变过程中可以不考虑液相区自然对流的影响;

    ③固体骨架与冻结相及液相工质处于局部热力学平衡;

    ④多孔介质骨架和流体混合物是不可压缩的。

    热平衡方程[3]:

    式中,ρ为多孔介质的密度;W为多孔介质中未冻结水的含量;W*为冻结水的含量;Cp为比热容,可由下式求出:

    式中,x为质量含量;L为潜热;下标s,w和i分别代表固体骨架、未冻结水和冰。   

当T < Tf时,

    式中下标tw表示总的含水量;b表示在-40℃时仍不会冻结的含水量。

    式中λc表示连续相的热导率;λd表示弥散相的热导率,在这里连续相为固体骨架,弥散相为水和冰;εd表示弥散相的体积占有率即体积空隙率,可由下式求得:

    边界条件方程:

    在表面:

    计算过程中所用到的物性参数如表1所示。

    2 数值计算与实验分析

    为了简化计算,考虑模型的对称性,如图1所示,将0~δ段用N个节点等分成(N-1)小段,节点1为中心,节点N为表面。分别在各个节点运用(1)~(9)中的方程,并在节点N和节点1处分别满足边界条件方程(10)和(11)。在计算过程中,计算出每一个时间步长时各个节点的温度值Ti和未冻结水Wi,然后根据新的温度值估算热物性参数重新计算。该计算过程使用MATLAB语言编程计算[6]模拟结果(如图2)。