再谈COMSOL中的固液相变
作为“浅谈COMSOL相变”系列文章三部曲的终曲，本文专门对COMSOL中固液相变的模拟方法进行展开介绍。（后面随着本人知识的积累，如有新发现，依然会持续更新相变的相关文章，谢谢大家的支持。）前面的两篇文章《浅谈COMSOL模拟相变的方法》和《再谈COMSOL中的气液相变》分别向同学分析了相变的难点和详细介绍了气液相变的方法，研究相变方向的同学记得一定要翻出来看看，肯定能让你有所收获。
目前有两大类主流的方法可用于处理相变导热的数值问题[1]，分别是界面追踪法（强数值解法）和固定网格法（弱数值解法），如图1所示。其中，COMSOL中自带有界面追踪法的处理功能，其需要与“变形几何”进行配合使用（理论上对于固液相变问题采用界面追踪法进行求解的精度最高，不过其也有一定的局限性，下文进行详细举例说明。）；COMSOL中也自带有固定网格法的处理功能，不过通常只能用于处理同种类的相变，即只能处理气液相变（液液相变）或者固固相变。（有不清楚如何区分“同种类的相变”的同学可详见此文章[2]），目前只有配合一些人为自定义的参量才能较好的在COMSOL中利用固定网格法处理不同种类之间的相变。根据对潜热处理方法的不同，固定网格法可进一步细分为焓法和显热熔法，根据焓法和显热熔法的自定义方式或者使用途径的不同，又进一步衍生出了拟源项法或等效热熔法等别名。
图1 相变数值解法的分类
下面本文选择对“锡熔融前沿”案例下手，为大家详细介绍在COMSOL中五种常用的处理固液相变的方法，按照介绍顺序，方法1为COMSOL自带的界面追踪法，方法2为COMSOL自带的固定网格法（本质应该是属于显热熔法。），方法3为在COMSOL自带的固定网格法基础上加入人为定义的部分，进行改良的显热熔法，方法4为在COMSOL中完全自定义的显热熔法，方法5为在COMSOL中完全自定义的焓法。（请大家理性看待案例，案例更多的是给大家一个启发，本人知识和精力有限无法完全保证案例的准确性。）
方法1COMSOL自带的界面追踪法。用此方法处理的“锡熔融前沿”案例可直接从COMSOL官网进行下载[3]。如图2和图3所示，需要同时设置“变形几何”和“相变界面”边界条件，才能启用COMSOL自带的界面追踪法。其中“变形几何”用于控制网格的变形（界面追踪法的本质就是通过网格的变形来“追踪”每一时刻的固液界面。），网格的变形情况由COMSOL根据用户在传热模块设置的边界条件自动进行计算（非常人性化。）。只有在设置了“变形几何”对应的边界上才能调用“相变界面”边界条件，其用于设置固液相变的相关参数，例如相变温度和潜热。
图2 变形几何
图3 相变界面
如图4所示为计算到10000秒时候的网格变形情况，如图5所示为计算到10000秒时候的固液分布情况，其中蓝色部分为液相白色部分为固相，如图6所示为计算到10000秒时候的液相流动情况，如图7所示为计算到10000秒时候的温度分布。结合图4和图5应该能更形象生动的理解“界面追踪法的本质就是通过网格的变形来‘追踪’每一时刻的固液界面”这句话的涵义，结合图6可以分析出在重力和高低温壁面的共同作用下，液态锡会形成一个顺时针方向的流动，最终导致固液界面形成一个曲面。理论上来说界面追踪法的精度最高、准确性最好，但是界面追踪法也有很明显的缺陷就是不能模拟固液界面的拓扑变化。具体展开说就是：1.不能无中生有，计算域在初始状态必须要同时存在固体和流体，无法处理纯固体融化以及纯流体固化的问题；2.不能模拟固液界面破碎、融合等情况；3.不能模拟完全融化以及完全凝固的情况，即在终止计算时间下也必须存同时存在固体和流体；4.处理固液界面大变形时收敛情况不好，虽然可以用自动重新划分来缓解，但是每次重新划分网格就会出现一次数据的奇点，不方便后处理数据的提取。因为界面追踪法的缺陷导致其只能处理少数的特定工况，因此才有学者提出了适应范围更广的固定网格法
图4 网格变形情况
图5 固液分布情况
图6 液相流动情况
图7 温度分布情况
方法2COMOSL自带的显热熔法。COMSOL自带的显热熔法不能考虑液相的自然对流，如图8-9所示为单纯采用COMSOL自带的显热熔法做出来的“锡熔融前沿”案例在10000秒时固液分布情况和温度分布情况，其各项边界条件与方法1对应的案例一致。对比采用方法1所做出来的案例，当不考虑液相自然对流的时候，固液界面是一条直线，明显不符合实际情况。不过COMSOL自带的显热熔法收敛性能好、设置简单，当遇到“模型计算域足够小”或者“液相粘度非常大”等可以忽略液相流动的情况下，COMSOL自带的显热熔法也是一种不错的选择。
图8 固液分布情况
图9 温度分布情况
方法3改良的显热熔法。当遇到方法1无法处理的情况而又需要考虑到液相自然对流的时候，就可能需要对方法2进行改良，需要在COMSOL自带的显热熔法的基础上增添一些人为定义的参量，以实现考虑液相的自然对流。如图10-12所示，分别为采用方法3所做出的“锡熔融前沿”案例在10000秒时固液分布情况、温度分布和液相流动情况，其各项边界条件与方法1对应的案例一致。从图中可以看出当考虑了液相的自然对流之后，固液界面呈现出曲面，其中曲面上部分液相位置靠前下部分液相位置靠后与采用方法1做出来的案例效果类似。从图11和图12可以看出液相在重力和温差的作用下也是做顺时针的流动，液相流速的最大值为0.02米每秒左右与方法1做出来的结果接近。综上所述，可以判断用方法3做出来的“锡熔融前沿”是符合客观的物理规律的。
图10 固液分布情况
图11 温度分布情况
图12 液相流动情况
方法4完全自定义的显热熔法。如图13-15所示，分别为采用方法4所做出的“锡熔融前沿”案例在10000秒时固液分布情况、温度分布和液相流动情况，其各项边界条件与方法1对应的案例一致。方法4也考虑到了液相的自然对流，与上述几种方法进行对比可以发现，固液界面所形成的曲面和液相的流动符合客观规律。对比方法3，用方法4做出来的固液曲面的效果更接近采用方法1所做出来的结果，从图12可以看到采用方法4所做出来的液相流动速度会更加偏慢一点（速度偏慢可能是由于一些参量没有细调导致，同学可根据实际仿真的结果进行调整。），不过与方法1所得到的速度结果依然属于同一数量级，偏差是在能接受的范围内。在实际做模型的时候发现，方法4的收敛性比方法3的收敛性更好。
图13 固液分布情况
图14 温度分布情况
图15 液相流动情况
方法5完全自定义的焓法。如图16-18所示，分别为采用方法5所做出的“锡熔融前沿”案例在10000秒时固液分布情况、温度分布和液相流动情况，其各项边界条件与方法1对应的案例一致。方法5也考虑到了液相的自然对流，不过所做出的固液界面的曲面形态不太理想（勉强能看到上端液面靠前下端液面靠后的效果。）。从图18可以看出采用方法5所做出来的液相流动效果是最接近采用方法1做出来的结果的，最大液相的流动速度基本一致！在实际做模型的时候，发现方法5的收敛性能也不错，但是方法5涉及的需要用户设置的参数较多，甚至有些参数很难找到其实际的物理意义，因此导致调试非常繁琐，也可能正式因为某些参数没有设置到合理的值才导致固液界面曲面形态不理想。
图16 固液分布情况
图17 温度分布情况
图18 液相流动情况
以上介绍了分别用五种方法分别处理“锡熔融前沿”所得到的结果，其实每种方法各有优缺点，就针对上文所探讨的“锡熔融前沿”这个案例所做出来的效果来说方法1、方法3和方法4都能得出较好的结果，其中理论上来说采用方法1所得出的结果精度最高。但是在现实情况中，锡更可能的是从完全固态的时候开始融化的（这种从完全固态开始融化的状态更符合一般情况。），那么方法1就无法进行求解了，此时只能采用固定网格法进行求解（理论上来说方法2-方法5都可以拿来求解。）。这里本文采取方法4对从完全固态锡开始融化的状态进行了模拟（由于方法2无法考虑到自然对流所以一开始就被排除使用，方法3做了很久发现收敛效果不好也被放弃使用，方法5由于固液界面形态不理想也被排除，最后只能采取方法4。），如图19所示为锡在不同时刻的固液分布情况，如图20为不同时刻液相的流动情况，从图中可以看到模型得到了较为“漂亮”的结果。
图19 不同时刻固液分布情况
图20 不同时刻液相流动情况
[1]党昕,孟多,高慧.焓法与显热容法在建筑相变蓄热技术数值模拟中的应用[J].辽宁工业大学学报(自然科学版),2021,41(03):188-194.DOI:10.15916/j.issn1674-3261.2021.03.013.[3] http://cn.comsol.com/model/tin-melting-front-6234
