导语

北京冬奥会在2022年2月4日开幕，开幕式上雪花的形象令人印象深刻。雪花生长过程能否用复杂系统建模软件来实现？本文用NetLogo软件的元胞自动机模型做了尝试。

研究领域：晶体生长，元胞自动机，多主体建模

1. 问题的提出

从在天安门广场中心布置的“精彩冬奥”主题花坛到冬奥会开幕式演绎一朵雪花的故事[1][2]，北京冬奥会“雪花”元素分外引人注目。

图片来源：

图片来源：

图片来源：S.Wolfram.A New Kind of Science.Wolfram Media,2002.

集智俱乐部复杂系统管理学读书会推荐书目《A New Kind of Science》中用元胞自动机按照“从一个六边形的黑色元胞开始，如果一个元胞相邻的元胞有黑色的话，这个元胞就变为黑色。”这一简单的规则得到31步斑图。

1月5日仇玮祎老师在自生成结构读书会交流群中推送了《为啥雪花长这样？》一文[3]。再回溯到2018年9月，集智俱乐部发布了《冰雪奇缘：为什么一粒灰尘可以长成六边形的雪花？》[4]

在《Netlogo多主体建模》课程中，张江教授主张：对复杂科学最好的入门手段就是自己亲手在计算机中搭建一个复杂系统。Netlogo就是一个非常好的入门工具，它可以让你通过简单的设置和代码编写就能搭建出一个超酷炫的多主体仿真（模拟）程序。

基于以上背景，开始追问“模拟雪花生长的奠基性模型是什么？如何在NetLogo平台上实现仿真？”

2. 模型的梳理

Packard[5]模型的元胞状态为 0 或 1，分别代表气相和固相，网格为三角形网格，每个元胞最近邻都为六个，表现了冰晶胞结构的六角特性。Packard 提出一种常用规则为只有一个邻居元胞为冰相时（记为 n＝1）才能生长，而 n≥2 的情形不能生长；另一个常用规则是用邻居中有奇数个为冰相则生长偶数个为冰相则不生长的规则（即：n＝1、3、5），也称为奇偶规则。虽然Packard模型规则极其简单，却能模拟出基本的板状、柱状和简单六角分枝状的雪花。但由于模型过于简单，其物理意义太过粗略。Gravner & Griffeath（2006）[6]对Packard模型进行了数学分析。

Libbrecht（2005）[7]着重强调，冰的表面层是“准液体”层（quasi-liquid layer），这是解开雪晶生长之谜的主要因素。Gravner & Griffeath（2008）[8]提出的二维模型采用了三种基本的物理量：雪花晶体部分，用a(x，t)表示是否属于晶体，c(x，t) 表示晶体的质量密度；准液体（quasi-liquid），用密度量b(x，t)表示；扩散水汽（vapor），用密度量 d(x，t)表示。在模型中这些量的演化过程可分为 5 类：扩散、凝结、附着、熔化、噪声。二维模型中随着初始密度ρ增加，在适当设置其它参数的前提下，雪晶生长形态表现出的不同特征。

图片来源：

图片来源：

Gravner & Griffeath（2008）[9]基于水汽的扩散、各向异性附着和准液层，引入了一个三维的、便于计算的、介观的雪晶生长模型。研究了若干能对大多数观察到的雪晶体形态忠实地模拟复现的案例。特别是，具有刻面和分支特征的物理标本最引人注目，模型解释了这一现象。模型还复制了许多其他观察到的特征，包括脊、肋、夹层板和空心柱，以及动态不稳定性。观测到的现象一致表明，模型中的因素在物理雪晶体的生长中是至关重要的。上图列出了主要的方程。

模型用C语言编写基本的计算引擎，用MATLAB进行绘图和可视化，用光线跟踪软件POV-Ray进行渲染。下图左边是用Matlab生成的图像，右边是在左图基础上用光线跟踪软件POV-Ray渲染后的效果。

图片来源：

3. 程序的设计

Packard模型的元胞状态、三角形网格、六个近邻表现了冰晶胞结构的六角特性，为后续模型所继承。Packard模型奠定了元胞自动机雪花生长模型的基础框架。Wolfram 在《A New Kind of Science》一书中和刘式达，张自国（1989）[10]运用Packard模型进行模拟仿真，为我们编程仿真提供了参考模式。

为构造六个邻居邻域的基本单元，在Netlogo正方形嵌块的基础上用循环嵌套画出有效嵌块组成的网格图。

其基本单元的相对坐标为：

以奇偶规则为例，要实现邻居中有奇数个为冰相则生长，偶数个为冰相则不生长的规则，需区别(t-1)时刻和t时刻状态。为此，在第三象限设置第一象限反对称镜像，按邻域关系在第一象限读取(t-1)时刻周边六个节点嵌块的颜色状态，并计算嵌块颜色为白色的数量作为判据。根据判据在第三象限反对称镜像改变t时刻由海龟坐标确定的当前嵌块的状态。然后用第三象限反对称镜像更新第一象限状态，为下一时刻规则执行做准备。

此外，要注意虚拟空间设置与界面输入必须匹配。

界面

设置

解决了以上三个关键问题后，就可以运用Netlogo四方格相关指令进行编程。程序代码为：

globals [ i j ]

to setup；；初始化

clear-all

create-turtles 1 [ setxy pxcor pycor ]

ask patches

[

ifelse ( pxcor = 32 and pycor = 32 ) [set pcolor white ] [set pcolor black ]

]

；；从第一象限copy 到第三象限

ask turtles

[

let i1 0

while [ i1 <= ( X / 2 ) ]

[

let j1 0

while [ j1 <= ( Y / 2 ) ]

[

setxy i1 j1

if ( [pcolor] of patch-here = white )

[

ask patch-at ( - 2 * i1 ) ( - 2 * j1 )

[

set pcolor white

]

]

set j1 ( j1 + 1 )

]

set i1 (i1 + 1 )

]

]

end

to go

let n 1；；步数

while [ n <= step ]

[

；；按计算坐标i，j

let j0 0

while [ j0 <= Y / 4 ]

[

set j ( 2 * j0 )

let i0 1

while [ i0 <= X / 4 ]

[

ifelse ( remainder j0 2 = 0 ) [ set i ( 2 * i0 ) ]；；偶

[ set i ( 2 * i0 - 1 ) ] ；；奇

ask turtles[ setxy i j ]

；；按邻域关系计算(t-1)时刻第一象限周边六个节点嵌块白色的数量

ask turtles with [ pcolor = black ]

[

let pcolor_w [pcolor] of patch-at ( - 2 ) 0

let pcolor_e [pcolor] of patch-at 2 0

let pcolor_wn [pcolor] of patch-at ( - 1 ) 2

let pcolor_en [pcolor] of patch-at 1 2

let pcolor_ws [pcolor] of patch-at ( - 1 ) ( - 2 )

let pcolor_es [pcolor] of patch-at 1 ( - 2 )

let sigma 0

if ( pcolor_w = white ) [ set sigma (sigma + 1) ]

if ( pcolor_e = white ) [ set sigma (sigma + 1) ]

if ( pcolor_wn = white )[ set sigma (sigma + 1) ]

if ( pcolor_en = white )[ set sigma (sigma + 1) ]

if ( pcolor_ws = white )[ set sigma (sigma + 1) ]

if ( pcolor_es = white )[ set sigma (sigma + 1) ]

；；在第三象限计算t时刻aij-t

if (pxcor <= X / 2 - 2 and pycor <= Y / 2 - 2)

[ ifelse (sigma = 0) [ ask patch-at (- 2 * i) (- 2 * j ) [ set pcolor black ] ]

[ ifelse (sigma = 1) [ ask patch-at (- 2 * i ) (- 2 * j ) [ set pcolor white ] ]

[ifelse (sigma = 2) [ ask patch-at (- 2 * i ) (- 2 * j ) [ set pcolor black ] ]

[ifelse (sigma = 3) [ ask patch-at (- 2 * i) (- 2 * j ) [ set pcolor white ] ]

[ifelse (sigma = 4) [ ask patch-at (- 2 * i) (- 2 * j ) [ set pcolor black ] ]

[ifelse (sigma = 5) [ ask patch-at (- 2 * i) (- 2 * j ) [ set pcolor white ] ]

[ ifelse (sigma = 6) [ ask patch-at (- 2 * i) (- 2 * j ) [ set pcolor black ] ]

[

]

]

]

]

]

]

]

]

show ( word "n" n "i" i "j" j "sigma:" sigma )

]

set i0 ( i0 + 1 )

]

set j0 ( j0 + 1 )

]

；；完成一步后从第三象限copy到第一象限

ask turtles

[ let i1 0

while [ i1 >= ( - X / 2 ) ]

[

let j1 0

while [ j1 >= ( - Y / 2 ) ]

[

setxy i1 j1

if ( [pcolor] of patch-here = white )

[

ask patch-at ( - 2 * i1 ) ( - 2 * j1 )

[

set pcolor white

]

]

set j1 ( j1 - 1 )

]

set i1 ( i1 - 1 )

]

]

set n (n + 1 )

]

End

4. 仿真的结果

下图为初始化之后运行的1～6步的雪花斑图。

5. 努力的方向

一是扩展现有编程。由于本程序使用的基本邻域单元不是正六边形，以及用正方形嵌块作为基本图素，故所得图形不能保证各向同性。可利用嵌块信息对相对雪花斑图中心的纵向坐标值乘以（

）进行压缩，实现各向同性。二是加强软件联动。现在已出现应用COMSOL仿真平台求解偏微分方程

（PDE）模拟雪花生长的模型[11][12] ，应及时发挥Netlogo与Mathematica等平台链接方便的优势，开展联合编程，拓展求解偏微分方程（PDE）等功能。三是跟踪前沿探索。Gravner J， Griffeath D[8][9]已经开发了更精准的雪花模拟模型，通过运用Netlogo 的pathes own[…]指令增加属性变量可方便地在本文提供的模型基础上进行拓展。Libbrecht[13]相信，借助越来越复杂且健全的计算机模拟，总有一天能从分子到原子甚至是到最为底层的量子力学中构建出一套完善描述体系。我们应及时跟进雪花建模的前沿，探索在Netlogo及其关联平台上的实现方式。

纪洪波 | 作者

邓一雪 | 编辑

商务合作及投稿转载｜swarma@swarma.org

◆ ◆ ◆

搜索公众号：集智俱乐部

加入“没有围墙的研究所”

让苹果砸得更猛烈些吧!