天文学中的数值模拟：编外篇
数值模拟这个系列很久没有写新文了，不说别的，这些文章写作的前提是要先挨篇去读原始文献，这个过程并不轻松。这学期倒还没有忙到一点空闲也抽不出的程度，不过没有心思做这件事是真。这回要说的是自己在最近搞的一些小型流体模拟，主要是为了计算天文课程的期末报告准备的。
其实这是本人第二次选修计算天文了，上次是在大四，因为种种原因没有参加期末考试，自然也就不能拿到学分。要说本人的研究生学分也早已满足学校的要求，不过还是对这门课耿耿于怀，于是这学期选了第二次。
课程报告的蓝本也是大四的时候打算做而没有做成的东西：用流体模拟的方法计算喷流的传播过程。由于个人计算机性能的限制，只能搞二维模拟，201*201格点，格点再多的话，系统就会报错说超过了数组允许的最大规模了，至于三维模拟那更是不可能了。有同学的课程报告居然搞出了个三维磁流体的不稳定性，惊讶之余连忙询问之，他的电脑到底是什么配置？答曰不是用的自己的电脑，用的是学校的并行机，边长八十格点，程序运行了整整一天。原来如此。
说来至少对于流体模拟，这二维与三维甚至更高维其实没有什么本质区别，关键还要看计算机能不能扛得住。当然高维的边界条件设置也会比较复杂，数值模拟的功底可能也就表现在这里了。比如这次另有个同学，模拟的是中央有障碍物的二维流场，仅仅添加几个形状规则的障碍就让程序的长度翻倍，那三维要复杂到何种地步就更是不敢想象了。
要问为什么格点数目是201这种不三不四的数字？原因也很简单，模拟区域的边长一般都会设成整数（本人这次设为2），如果是200格点，那么两点的间距等于2/(200-1)，计算起来不就不太方便，也容易引入误差了吗？而如果取上201点的话，间距倒是整齐的0.01了。
先是在均匀背景下传播的高密度喷流，这项计算主要是为了进行程序的测试：
其中白色表示密度最高的区域，颜色越深密度越低；箭头是速度场。除了喷流本身的行为，以声速传播的扰动也是表现得相当明显。至于压力变化就不亮相了吧，实在有点难看，而且说实话，初始条件等压的话，演化过程中压力的变化也不很大。
最初的想法是打算模拟喷流速度与背景密度梯度方向相交时的喷流传播，预计会出现拐折，不过实际计算表明，这个拐折非常的不明显，可能是参数选取的问题，所以动画就不放出来了。然后在接下来的结果亮相之前，先扯一扯牛人的大旗给自己壮壮胆：
图片提供：Weiqun Zhang and Stan Woosley
这是Woosley的小组完成的伽玛射线暴中心能源模拟。根据现在的主流模型，中心能源是坍缩星，恒星中心先坍缩成黑洞，随后形成吸积盘与极向的相对论性喷流。喷流穿过恒星外包层，如上图所示冲出星体表面之后，才能形成爆发并被观测到。由于时间紧迫，这次本人是先对理想流体进行了模拟。相对论性流体遵循的方程和理想流体不同，象Woosley这样的工作更多考虑的也是能量流而非物质流，所以本人当下的程序还不能用来模拟伽玛暴这类极端现象。不过至少可以看一看速度较慢的喷流穿出致密包层时的行为不是？
第3张图中，喷流的初始速度是第2张的3倍（当然时间步长也不相同）。可见不稳定性还是很明显的，也就是由于瑞利—泰勒不稳定性导致的手指结构（喷流前端）。喷流前端两侧的对称的牛眼结构则可以归结为涡流，慢速喷流可能还伴有更强烈的不稳定性。
当然，类似的不稳定性结构并非天文背景的专利，地球上照样存在，只要略加注意就能发现，比如下面这个例子：
其实上图中的牛眼跟本人计算的第一张图没有任何本质区别。既然计算过程中，流体速度已经是用声速作过归一了，只要假设不同的声速，自然可以应用到不同的场合。
个人的感觉是，数值模拟的主要难点并不在于程序的编写，而是调试。差分方程形式很琐碎，稍有不慎某一项输入错误，那结果就差了十万八千里，因此反复的检查是逃不掉的。由于这次本人采取的格式振荡较大，对于人工耗散项也不能偷懒（实际情况是想偷懒未果，还是老老实实地加了上去），话说这耗散系数的调节也是学问哪，过小则无法补偿振荡，过大照样出问题，至于究竟取到多少才合适？只好挨个试验……
至于什么时候才能把这个方法套用到相对论性流体就说不好了，只能看时间安排。虽说对相对论性流体力学方程进行差分并不是太复杂的工作，但整个过程要牵涉到大量的其他事宜，因此，烂尾也是极有可能的事情，走着瞧……
