作者：蒋迅

这是我的文章“雪花里的数学：四，元胞自动机模型”的后续。本文已发表在《数学文化》杂志2012年第4期 (总第12期) 上。请勿转载。

五，相变的有限元解

Source: Scientific American

还有一个研究以相变 (phase transition) 为出发点。相变是指物质在外部参数（如：温度、压力、磁场等等）连续变化之下，从一种相（态）忽然变成另一种相，最常见的是在一定的条件下，冰变成水和水 变成蒸气等，也有可能是相反的过程。我们把这样的过程称为相变。因为水和冰之间的边界不是褂讪的，所以它形成的热传导方程是一个自由边界问题。对这个自由边界的最简单描述 (或者说，最简单模型) 就是在这个界面上，温度为零摄氏度。让我们考虑一个最简单的一维情况。这个情况和本文讨论的雪花问题不完全一样，但是也许可以帮助读者加深理解。假设在 [0，+∞] 区域上的冰和水。假定在开始时整个区域都是冰。我们从左边提供一个热源，于是冰开始熔化，在 t 时刻，区域 [0，s(t)] 变成了水。忽略边界条件和初始条件，我们得到在 [0，s(t)] 上的热传导方程：

∂u / ∂t = ∂²u / ∂x²       (x, t): 0 < x < s(t), t > 0

这里，u = u(x, t) 是温度，s(t) 是自由边界。在自由边界上温度是零，所以 u(s(t)，t) = 0。注意自由边界 s 是随时间而变的曲线 (或曲面)，我们还应该有一个在 s 上的条件。最常见的就是著名的“史蒂芬条件”(Stefan Condition)：

ds / dt = - ∂u / ∂x(s(t), t)             t > 0

也就是说，自由边界随时间的变化率和温度在自由边界上位移的变化率成正比，方向相反。相应的偏微分方程就是“史蒂芬问题”(Stefan Problem)。为了说明相变的性质，我们再稍微深入一步。物理实验表明，在自由边界上，温度达到零度，但不会立即继续升温。这里有一个积蓄能量的过程，直到增加了 L 单位的能量 (潜热) 后温度才会继续增长。让我们引入一个新的变量 y 来表示水的浓度，并假定 L = 1。我们定义 v = u + Ly = u + y。这里，y 是阶梯函数：y = 1，如果 0 < x < s；y = 0，如果 x > s；0 < y < 1，如果 x = s。

v 决定了相变的热动力。引入函数 h(z) = min (z, 0) + max (z-1, 0)。则上面的热传导方程可以写成：

∂v / ∂t = ∂²v / ∂x²

经过这个变换，上述方程在整个半实数轴 (0，+∞] 上成立，温度u = h(v)。方程本身成为退化的抛物形偏微分方程。这样做的好处在于，人们可以运用变分和有限元的方法在一个褂讪的区域里得到方程的数值解。

上面的讨论不是一个严格的讨论，只是希望帮助读者理解下面要介绍的相变的有限元解法的思想。 从这样的描述看，雪花的形成问题应该是与相变问题紧密相关的，因为在第三节里我们已经看到，雪花是由水珠在一定的温度和饱和度条件下形成的，从这个意义来 说，就是一个结晶的过程。从数学上说，人们需要做的就是研究自由边界表面随时间的变化。一片雪花的形成过程是否也能用相变的数学模型来描述呢？

2012年初，英国伦敦帝国学院的约汉·巴瑞特 (John Barrett) 教授和罗伯特·纽伦伯格 (Robert Nurnberg) 教授与德国雷根斯堡大学的哈罗德·加克 (Harald Garcke) 教授就按照这个思路做出了一些新的工作：“雪晶体生长中分面格式形成的数值计算”(Numerical computations of facetted pattern formation in snow crystal growth)。这是他们在自由边界问题的有限元分析的成果之上对雪花研究方面的一个有意义的新尝试。

让我们先回到在第三节中的“雪花形状与温度、饱和度的关系”那张图。很明显，雪花的形状与温度和饱和度有关。当温度刚刚在冰点之下的时候，如果饱和度比较低的话，出现柱状雪花；如果饱和度比较高的话，就出现树突状雪花。当温度在-5^oC附近时，如果饱和度比较低的话，出现实心柱状雪花；如果饱和度比较高的话，就出现空心柱状雪花和针状雪花。当温度低到-10^oC以下时，如果饱和度比较低的话，出现实心盘状雪片；如果饱和度比较高的话，就出现树突状雪花。当温度到-25^oC 以下时，如果饱和度比较低的话，出现实心盘状雪片；如果饱和度比较高的话，就又出现柱状雪花。从物理意义上说，雪花的形成过程是固体和气体的边界（即自由 边界）变化的过程。而这个自由边界的变化是由于水分子的扩散和附着等过程。在这个过程中满足物质守恒定律，同时表面能量达到极小。另外我们知道，冰是一种 六方晶系的晶体, 基本形态是六角形的片状和柱状。冰晶体的各向异性导致其物理性能（如导热性）随着方向的不同而有所差异。这种六边形的各向异性（hexagonal anisotropy）也必须考虑进去。巴瑞特、纽伦伯格和加克根据上述条件引入了一个与雪花相关的结晶 (Crystallization) 的数学模型 ─ 准静态的扩散问题 (quasi-static diffusion problem)。

他们的巧妙之处在于把雪花的轮廓看成是一个自由边界，然后在自由边界上施加具有物理意义的复杂条件。对相应的偏微分方程建立其变分形式，以便运用有限元方法找到数值解。为了实现不同形状的雪花，在他们通过改变方程组里的饱和参数u_D = u_∂Ω、凝结系数 β 和表面能量六边形的各向异性参数 γ 来实现。 我们同样不准备对他们的方程组进行详细的讨论，而是转到数值计算上来。

uD = 0.004, γ = β = γhex.
t = 0; 5; ...; 50	t = 0; 50; ...; 500
uD = 0.01, γ = β = γhex.
t = 0; 5; ...; 50	t = 0; 50; ...; 200
uD = 0.04, γ = β = γhex.
t = 0; 0.5; ...; 5	t = 0; 5; ...; 40

在建立了一组偏微分方程之后，他们进而用有限元方法进行了数值计算。上图就是他们数值计算的一些结果。为了方便观察雪花界面（即自由边界）随时间的变化，每一行中并列的的两个图分别记录了这个曲线在不同时间段的形状。注意这里的每张图都是多个界 面的叠加图。我们看到，除了参数u_D以外，其他参数都是相同的。这三个雪花分别由 u_D = 0.004，0.01和0.04 生成。从而从数学上解释了雪花形状和饱和度的关系。在他们的论文中有两类雪花：刻面 (facet) 和树突 (dendrite)。这些图形都是在不同的参数选取下得到的。加克对科学美国人记者说，他们“是第一个用能量守恒和热力学理论同时实现这两种生长”的小组。

Source: Harald Garcke

为了在计算机上模拟这组偏微分方程代表的雪花的生成，人们必须准确地描述结晶面 (即自由边界) 随时间的变化。人们通常是把这个曲面用不断加细的三角形来近似，这是有限元法所必须的。但这些三角形经常会退化从而导致模拟失败。他们的办法就是用现在比较时髦的平均曲率 (Mean curvature) 来控制模拟以达到在计算机上实现的目的。他们表示，这个办法可以避免三角刨分蜕化的难点。笔者认为这是他们成果的一个亮点。

Source: Harald Garcke

Source: Harald Garcke

他们也对偏微分方程的数值解做了分析。他们发现，结晶体中表面分子的结合对结晶的生长有很强的影响。他们还发现，雪花尖端的生长速度与空气中的水蒸气的多少成正比。他们认为，雪花的结构源于扩散有限晶体生长在各向异性的表面能量和各向异性吸附动 力。冰晶形态的稳定在很大秤谌上依赖于饱和、晶粒尺寸和温度。他们注意到了尖端速度和饱和度之间有线性关系。他们还得出结论，表面能量的影响尽管很小，却 对雪花的形成有较大的影响。最后一点最为重要，它也许揭示了一个可能最后解决雪花形成问题的新的思路。

他们的模型 ─ 第一次用能量守恒和热力学理论建立的连续模型成功地研究了雪花的生成，这是此模型与格拉夫纳-格里夫耶斯模型的本质区别。由于这个原因，他们的途径自然地被物理学家所欣赏。重要的是，他们开辟了用偏微分方程和有限元方法研究雪花的新方向。麻省理工学院的“科技评论”(Technology Review)、英国的“邮件在线”(Mail Online) 以及“科学美国人”(Scientific American) 都报导了他们的成果。

在这里，笔者需要对“科学美国人”的报导做一点说明。“科学美国人”在文章中也提到了格拉夫纳和格里夫耶斯的工作，但把他们的工作误解成了在分子的层次上的格状自动机模型。其实，格拉夫纳和格里夫耶斯的工作是在介观体系的层次上，也就是说他们只是到了微米的范围。他们在论文标题上就写清楚了这一点。事实上，目前对雪花的形成人们还没有一个完美的解释。从物理上还无法在分子的层次上解释分子的附着和分离的机制。利伯布莱切特认为，可能有某个物理性质还没有被认识到，比如雪晶形状的不稳定性。在这样的情况下，在分子的层次上建立数学模型也就无从提起了。

有一个有意思的巧合是，德国雷根斯堡大学正是第一个研究雪花的天文学家和数学家开普勒逝世的地方。加克教授说，他们可以在办公室的窗前就见到真正的雪花。这是在加州洛杉矶的利伯布莱切特教授和在加州戴维斯的格拉夫纳教授不能享受的优越待遇。

待续。