利用COMSOL复现SCI二区论文有感
《Numerical study of melting inside concentric and eccentric horizontal annulus》的基本复现——利用COMSOL复现SCI二区论文有感
作者：极度喜欢上课
近段时间本人不是正在复现论文就是在复现论文的路上，在看了大量流体力学、传热、相变相关论文以及尝试复现过一定数量的论文之后，有很多感悟一直想和大家分享，但时间总是太匆忙想分享的内容一转头就忘，所以现在趁此机会，以一篇SCI二区论文《Numerical study of melting inside concentric and eccentric horizontal annulus》[1]为例，好好跟大家分享一下近期我在复现论文过程中的一些心得体会，也希望能给大家带来一定的启发。（想说的话有很多，本人决定采取最朴素的方式，从头到位开始叙述。）在上完第六次COMSOL网络课程[2]之后，感觉本人目前所掌握的利用COMSOL处理相变问题的方法还有一定精进的空间，所以网课之后本人一直都有在找传热相变相关的论文进行阅读并积极尝试复现。（现在本人传授给各位如何查找合适用来复现的论文的经验，总的原则就是“复现别人复现过的论文”。）刚开始查找英文文献的时候本人通过关键词“COMSOL”“PCM”“phase change”在百度学术上进行检索，然后通过论文标题判断是否有下载该篇论文的必要，阅读已经下载好的论文的时候则着重看它是如何验证论文模型有效性的，（较为有效的验证模型有效性的方法是自己做实验对比或者参考前人的数据。），专门看它有没有参考前人的数据，它所参考的前人的数据就是一篇极大概率的“可被复现的论文”！《Numerical study of melting inside concentric and eccentric horizontal annulus》[1]就是通过上述方法被我找到的。一开始本人是下载了《Numerical investigation of heat transfer and melting process in a PCM capsule_ Effects of inner tube position and Stefan number》[3]这篇论文，该篇论文在验证模型准确性时总共参考了两篇论文的数据，其中有一篇就是《Numerical study of melting inside concentric and eccentric horizontal annulus》，然后本人就下载了该篇论文进行阅读，发现其提供的边界条件较为齐全，有较大的进行复现的可能性，所以它才被挑选为本人复现论文道路上的一块炼金石。一篇论文少则七八页多则几十页，每个人看论文的目的都不一样，所以一定要有针对性的看论文，才能尽可能的节省时间和精力。如果我们是以复现论文为目的的，最重要的则是找到论文中所给的建立模型的条件，如图1所示，就是《Numerical study of melting inside concentric and eccentric horizontal annulus》中所给出的主要的条件，包括：模型几何尺寸，材料物性参数以及必要的边界条件。（当然作为一篇采用自定义焓法处理相变问题的论文，它所提供的自定义方法我们也需要看懂，只是该篇论文用的是fluent进行仿真，本人打算用COMSOL进行复现，所以本人还是打算使用在COMSOL里面较为容易实现的自定义显热熔法进行复现。）
论文的复现基本不太可能一次性就能成功的，就这次的复现而言如图2所示，本人总共建立了6个COMSOL模型，耗时约3天才完成了基本复现。
如图3左所示为第一次尝试建立起的COMSOL模型所计算出来的第十分钟时的相变结果（此mph文档为2023.01.02），图3右为参考文献同种工况下第十分钟时的相变结果。通过对比可以发现“2023.01.02”的结果接近参考文献，但是还存在较大差别，主要的差别在于“11点”和“1点”方向的液态石蜡没有触碰到大圆的边缘。看到了差别那么就要开始分析问题了，一通分析下来有本人发现可能是两个原因导致的，1.就是融化速度可能偏慢，一些参数的取值存在一定的范围，可能在建立“2023.01.02”的时候没有取到一个较优值，所以导致相变材料融化偏慢；2.就是几何尺寸有可能有问题，文中提到的中间小圆有一个1.5mm的壁厚，在建立“2023.01.02”的时候可能没有考虑到壁厚导致中间的小圆半径画小了。
既然分析出了潜在的问题，那么就要一步一步验证了，首先验证原因1。如图4右所示为“2023.01.02”的液相率，图4左为参考文献同种工况下的液相率。（由于第一次尝试的模型结果明显不对，所以在计算到十分钟时就停止了计算。参考文献对应的曲线为Array1。）“2023.01.02”在第十分钟时液相率不到0.45，而参考文献的约为0.6，明显偏慢了。
验证原因2。除图1外参考文献其他地方基本没有进一步提及或者解释几何尺寸的地方了，所以如果要验证原因2需要发挥一定的主观能动性。如图5（a）为“2023.01.02”中的几何模型截图，图5（b）为参考文献出现的几何截图。图5（c）为参考文献出现的结果图截图，图5（d）为图5（a）、图5（b）和图5（c）的重合图。在图5（d）中三种图几乎完美的重合在一起，通过对比应该可以确认“2023.01.02”的几何尺寸是没有问题的。（为了更好的进行对比，图5（d）将三幅图的透明度更改为20%。）
上述分析基本确认了较大可能是由于原因1导致“2023.01.02”所计算出的结果和参考文献的结果存在较大差别。既然确认了大概方向，那就开始进一步调试和验证，通过调整“2023.01.02”的有关参数，适当增加石蜡的融化速度。如图6左所示为第二次尝试建立起的COMSOL模型所计算出来的第十分钟时的相变结果（此mph文档为2023.01.02.1），图6右为参考文献同种工况下第十分钟时的相变结果。通过图6可以看出，“2023.01.02.1”比“2023.01.02”更接近参文献的结果，可见上述分析的方向应该是没有错的。
“2023.01.02.1”所得到的结果和参考文献的结果还存在一定差距，所以还可以进一步调试。如图7左所示为第三次尝试建立起的COMSOL模型所计算出来的第十分钟时的相变结果（此mph文档为2023.01.02.2），图7右为参考文献同种工况下第十分钟时的相变结果。从图7中可以看出已经是非常接近参考文献的结果了。
如图8右所示为“2023.01.02.2”的液相率，图8左为参考文献同种工况下的液相率。“2023.01.02.2”的液相率随时间变化的规律与Array1基本一致，液相率曲线在时间为20min节点时都出现了明显拐角，而且液相率等于1的时间节点也大致相同，不过仔细对比一下可以发现“2023.01.02.2”的融化速度比参考文献的稍慢。
综合分析图7和图8所反映的情况，以及考虑到采用的是不同仿真软件以及相变处理方法处理的结果，本人基本认为“2023.01.02.2”是复现成功的。为了进一步验证模型的稳定性以及参数的独立性，本人在“2023.01.02.2”的基础上调整了有关参数、网格大小以及容差，分别额外进行了三次调试，形成了mph文档“2023.01.02.3”、“2023.01.02.4”和“2023.01.02.5”。额外三次调试的结果都表明“2023.01.02.2”是基本正确的，能较好的与参考文献进行对应。（这里还是要呼吁广大同学要理性看待COMSOL案例以及论文复现，由于精力、学识有限，本人无法保证自己所做的案例百分之百正确无误。本文也只是将复现的结果进行客观的呈现。）由于复现过程中做了较多的模型，为了防止日后校对的时候忘记每个模型的具体调试过程，所以本人专门将每次调试的具体内容写在了text文档上，如图9所示。如同学们有需要获取本次复现过程所涉及的案例，可加QQ：2581349746进行联系。
[1]Darzi A , Farhadi M , Sedighi K . Numerical study of melting inside concentric and eccentric horizontal annulus[J]. Applied Mathematical Modelling, 2012, 36(9):4080–4086.[2]https://www.bilibili.com/video/BV1QW4y1v7TN/[3]Bouzennada T , Mechighel F , Filali A , et al. Numerical investigation of heat transfer and melting process in a PCM capsule: Effects of inner tube position and Stefan number[J]. Case Studies in Thermal Engineering, 2021, 27(1–2).
