数值模拟的重现性

曾泳春

       我在这几年间，陆续收到过几次国内外的邮件，询问关于聚合物纤维运动模型的建模和计算问题。邮件应该都来自学生，学生承接了这样的课题，查文献，看着不算难，于是埋头开始做，做着做着就觉得难了。数值模拟就是这样，模型看懂了，照着人家的模型做，可就是计算不出来。几年前我有篇paper刚上线，就有个记不得哪个国家的学生写来邮件，说太好了，他正想做这个，能不能把计算程序给他？我直接就晕菜了，不知教授们收到这样的请求会怎样处理？按理说模型有了，程序也必然产生了，只是你得去写代码，对吧。写代码调程序这样的活，相信很多博士生心有余悸，那是要熬心血的。我去NC state访学时，认识了一个北航去的博后，我想着北航博士嘛，流体力学一定比我强得多，就赶紧抓着人家虚心请教一些流体力学的问题。那博后抹着头上的汗珠子说：整个博士阶段光写代码了，都没去深究那些物理意义。然后恶狠狠地说，博后再也不做模型了，做实验，用死他们的仪器设备猛发paper。说回那个要计算程序的邮件，我想了想，给他回复说：代码是学生写的，学生不同意给。

    写这篇博文是因为最近又有一个学生写邮件来求教静电纺丝射流不稳定运动的模拟。他说采用的是Reneker的模型（晓明同学一定知道Reneker），建立了动量、质量和电荷守恒方程，但在计算中发现库仑力远远大于其它力，包括重力、电场力、粘弹力和表面张力，它们不在一个数量级上，于是他算不下去了。而我大约在4、5年前也做了这个工作并发了paper，于是他就来问了。

    说到这个工作，要追溯到我的博士论文工作。我博士论文研究的纤维运动，是指纤维成纱过程中的运动。我当时借鉴了聚合物动力学中的大分子模型，从最简单的哑铃模型到珠—簧模型进行数值模拟，早就经历了这个学生所经历的问题。当我把模型的所有方程写好，踌躇满志地开始计算时，发现那些力们不在同一个数量级上。这真是直接崩溃，我博士期间流的眼泪，就是为这些力们不在一个数量级上而流的。所以这个学生来求教这个问题，我深深感知他崩溃的心情。而我后来转到聚合物纤维成丝的领域，是因为发现Reneker的模型和我的纤维模型很相似，只是将珠子之间的连接从弹性变成粘弹性。而且Reneker的模型解决了困扰了我很久的一个缺陷，即我是用挠度来计算纤维的弯曲变形，无法用于大的形变。

    说回这个来信求教的同学，他说采用Reneker的模型，只是参数有些不同，但无法重现模拟结果。这个问题相当有趣，有些像这两天大家在热议的“日本美女科学家”的实验工作，实验不能重现，经常受到质疑，而数值模型不能重现，好像还没听说过受到质疑，毕竟，实验研究远远多于数值模拟研究。数值模型一旦建立了方程，应该说不会出现不能重现方程的解的问题，但实际中却屡屡出现。我从几个方面来分析：（1）模型非常依赖于参数，参数的变化会引起解的剧烈变化。而Reneker的模型正是如此，事实上他的文章中也给出了两组参数所得到的差别很大的纤维运动轨迹。（2）有量纲和无量纲的计算，这是我特别迷惑的一个问题，数值计算都必须采用无量纲吗？有量纲的计算会产生什么问题？（3）编程中的处理。一般paper只会写出数理方程和计算步骤，但在编程中，往往会有一些处理。比如我用珠—簧模型计算时，因为实际纤维在无约束的运动中不可能伸长，所以我在每一step之后，都会把伸长的纤维往回拉，这里面有些小技巧，是不会写到paper中去的，也因此如果有人去重复我的模拟结果的话，我相信是重现不了的。

    小议一下，给有兴趣的人讨论讨论。

另：终于写了一篇基于科学网主题的博文来点缀一下博客，学习王晓明同学，也松一口气，闭关喝茶！同梦们玩微信的玩微信，玩博客的博客，老文在微博天天给自己灌心灵鸡汤，大家使劲玩！

蔡淳佳，《等一个晴天》。

等一个晴天 我们会再相见  
你说了 风吹我就听见  
笑着说再见 就一定会再见  
心晴朗 就看得到永远