在极地复杂环境下，进行科学考察的船舶和极地保障平台会面临结冰问题。过冷水滴冻结累积形成的冰层将对外部设备造成危害。现有的船舶结冰研究多集中于水滴冻结或船舶上层建筑结冰量计算^[1-5]，多忽略了冰的微观生长过程，而冰晶生长行为的研究有助于深入揭示结冰微观机理及过程。六角冰晶也称枝晶，是凝固过程中最常见的一种晶体生长形式，六角冰晶粒相互冻结在一起，最终形成表面光滑的冰壳^[6]。在枝晶生长问题中，一大难题是对固液相变界面的处理，界面的移动往往很难追踪^[7]，而相场法通过引入一个宽度非常小的过渡层来代替跳跃间断或者尖锐界面，避免了复杂的界面追踪，不仅可以通过耦合相场与外场方程将微观与宏观结合起来，还能模拟分析晶体生长过程中的物性参数对冰晶生长的影响。相场法最先应用于金属模拟领域，而后经过发展，已广泛应用于材料科学、固体物理、电化学等领域。Kobayashi^[8]通过数值求解了考虑各向异性的二维相场模型，并得到了过冷纯金属凝固的二维复杂枝晶生长形貌。Wheeler等^[9]提出了模拟二元合金等温凝固的相场模型(WBM模型)。Wheeler等^[10]模拟计算了过冷镍的凝固过程，结果与Ivantsov理论相一致。1999年，Kim等^[11]建立了可用于合金模拟的相场模型(KKS模型)，消除了WBM模型中界面厚度限制。2005年，Qin等^[12]建立了多元合金凝固模型，与实验结果基本相符。

近年来，国内学者开始将相场法运用到水的凝固。陈梅英等在食品存储领域的冰晶生长有较多的研究^[13-16]；李方方等^[17]实现了细胞尺度下的冰晶数值模拟；邹阳^[18]模拟了制冷剂中冰晶生长过程，结果与实验观察相一致。但这些都没有与真实物理场建立较好的联系。徐立等^[19]讨论了无量纲过冷度对冰晶生长的影响，但没有对Wheeler模型的准确性及各主要参数意义进行讨论分析。

本文将基于Wheeler模型，采用海水物理参数来实现海水冰晶的数值模拟，为过冷水滴冻结的再辉阶段(即枝晶生长过程)提供理论基础，进一步揭示结冰机理，并为模拟船舶结冰过程及评估极地船舶安全提供技术参考。

1 Wheeler相场模型 1.1 自由能函数形式及控制方程

Wheeler相场模型能够将相场参数与实际物理参数联系起来，使模拟结果更加真实^[20]，这也是本文选取该模型的原因。该模型以金兹堡-朗道理论为基础，引入相场变量ϕ，ϕ=0表示固相，ϕ=1表示液相，其采用的自由能函数形式为：

$\begin{aligned} f(\phi, T)=& W \int_{0}^{\phi} p(\phi)[1-p(\phi)] \cdot \\ &\left[p(\phi)-\frac{1}{2}-\beta(T)\right] \mathrm{d} p(\phi) \end{aligned} $

(1)

式中：W为常数，表示单位体积的能量；β(T)是关于温度的函数，表示热力学驱动力，β(T) < 1/2；p(ϕ)为固相分数，1-p(ϕ)为液相分数，p(ϕ)=ϕ³(10-15ϕ+6ϕ²)。f(ϕ, T)为双稳势函数。

相场及温度场的最终控制方程为：

$\begin{array}{*{20}{c}} {\frac{{\overline \varepsilon }}{m}\frac{{\partial \phi }}{{\partial t}} = \phi (1 - \phi )\left[ {\phi - \frac{1}{2} + 30\bar \varepsilon \alpha \Delta u\phi (1 - \phi )} \right] - }\\ {\frac{\partial }{{\partial x}}\left( {\varepsilon (\theta ){\varepsilon ^\prime }(\theta )\frac{{\partial \phi }}{{\partial y}}} \right) + \frac{\partial }{{\partial y}}\left( {\varepsilon (\theta ){\varepsilon ^\prime }(\theta )\frac{{\partial \phi }}{{\partial x}}} \right) + }\\ {\nabla \cdot \left( {{\varepsilon ^2}(\theta )\nabla \phi } \right)} \end{array} $

(2)

$\frac{\partial u}{\partial t}+\frac{p^{\prime}(\phi)}{\Delta} \frac{\partial \phi}{\partial t}=\nabla^{2} u $

(3)

式中：ε为无量纲界面厚度；m为相场迁移率；α为界面自由能；Δ为无量纲过冷度；ε(θ)是关于界面法向与参考坐标轴的夹角θ的函数；ε(θ)=ε(1+γcos(jθ))；γ为各向异性强度；j为各向异性模数；函数ε(θ)的引入能真实反映枝晶生长各向异性。ε、m、α、Δ与真实物理参数的关系为：

$\bar \varepsilon = \frac{\delta }{\omega }, m = \frac{{\mu \sigma {T_M}}}{{\kappa L}}, \mathit{\Delta } = \frac{{c\mathit{\Delta }T}}{L}, \alpha = \frac{{\sqrt 2 \omega {L^2}}}{{12c\sigma {T_M}}} $

(4)

以上各物理参数均为水的物性参数，具体含义及大小见表 1。诸如界面厚度等涉及尺寸的参数本应取真实值以模拟实际枝晶生长，但从计算角度考虑，需要将界面厚度等扩大到可计算范围内，本文引入长度参考尺寸ω(取2.7×10^-4 cm)，无量纲界面厚度ε乘以参考尺寸ω即得真实界面厚度δ，网格间距及计算域大小同理推得真实尺寸。

1.2 溶质场控制方程

海水冰晶可仿照二元合金，视为纯水与盐的混合物，因此其与纯冰晶相比，需要考虑浓度场影响，在模型中引入浓度主要有2种方式：1)在自由能密度函数中包含浓度参数；2)考虑浓度对过冷度的影响^[21]。本文是在纯冰晶Wheeler相场模型的基础上，考虑浓度对过冷度的影响，相场及温度场方程保持不变，基于Fick定律引入溶质场方程为：

$\frac{\partial C}{\partial t}=\nabla \cdot D\left[\nabla C+\frac{\left(k_{0}-1\right) C}{P+k_{0}(1-P)} \nabla P\right] $

(5)

$D=D_{s}+\left(D_{l}-D_{s}\right) \frac{P}{P+k_{0}(1-P)} $

(6)

式中：C为海水中盐的浓度(摩尔分数)；D_s为溶质固相扩散率；D_l为溶质液相扩散率；k₀为溶质平衡分配分数；溶质场方程中浓度C采用海水中盐的摩尔分数，已知海水中盐的质量分数约为3.5%，可求得盐的摩尔分数为0.011。

采用显式差分格式即可求解每个时间步内控制区域的浓度场，并通过式(7)影响区域内过冷度分布：

${\mathit{\Delta }^\prime } = \frac{{{C_p}\left( {\mathit{\Delta }T + {m_L}\left( {C - {C_0}} \right)} \right)}}{L} $

(7)

式中：Δ′为浓度影响下的实时无量纲过冷度，其值在每个时间步迭代完成后更新，进而影响相场控制方程；m_L为液相线斜率，K/(mol%)；ΔT为实际过冷温度。

1.3 噪声强度的引入

在实际的海水凝固结冰的过程中，受到外界的干扰是不可避免的，因此在相场中引入噪声项来模拟固液界面的干扰，其干扰会一直存在，造成界面处发生失稳现象，从而导致二次分枝、三次以及四次分枝的出现^[22]。本文在相场控制方程中加入噪声项，记为A_n(A_r, R_n)：

$A_{n}=-16 \phi^{2}(1-\phi)^{2} A_{r} R_{n} \times 30 \bar{\varepsilon} \alpha \Delta T \phi^{2}(1-\phi)^{2} $

(8)

式中：A_r为噪声强度，描述引入扰动的大小；R_n为(-1, 1)的随机数，符合高斯分布。

2 数值求解 2.1 控制方程离散形式及稳定性条件

有限差分法在均匀网格中进行离散求解，具有原理简单、方程便于建立和易于编程等优点，本文将采用有限差分法求解控制方程，其中时间离散采用一阶向前差分，空间离散采用中心差分，通过九点差分格式离散拉普拉斯算子。

显式差分格式可满足相场及溶质场控制方程计算的稳定性，而温度场方程运用ADI交替显隐式格式求解，其在二维模拟中无条件稳定，因此仅考虑相场及溶质场稳定性条件，本文控制方程中增加的非线性项对时间步长要求更严格^[23]：

$\Delta t \leqslant \frac{(\Delta x)^{2}}{5 m} $

(9)

式中Δx为网格空间尺寸。

2.2 初始条件及边界条件

结晶温度总是低于融化温度，即液体开始结晶时需要一个冷却过程，同时液体中还必须有一定量的微小结晶核，结晶核的存在能够缩短过冷却过程^[24]。晶核设置如图 1所示。本文设置液体区域初始无量纲温度为-Δ，而无量纲融化温度为0，从而形成过冷水区域，同时通过在特定计算区域设定固定半径的圆来引入晶核并设置初始条件：

$x^{2}+y^{2} \leqslant r^{2}: P=0, T=0, C=C_{0} $

(10)

$x^{2}+y^{2}>r^{2}: P=1, T=-\Delta, C=C_{0} $

(11)

相场、温度场及溶质场均采用zero-Neumann边界条件，即: $\partial P / \partial n=\partial T / \partial n=\partial C / \partial n=0$。

2.3 计算优化及参数设置

本文利用Python语言下的Scipy科学计算库求解三对角矩阵，简化了计算代码。鉴于经典的枝晶生长理论都是关于枝晶尖端的局部模型，进行定量分析时主要关注枝晶尖端特征数据，因此本文在计算区域底部中心设置形核并设定枝晶相场值沿生长主轴左右对称，节省了计算时间。本文设定网格区域N_x×N_y为300×300，网格间距Δx及Δy均为0.03，时间步长dt为0.000 1 s，晶核半径0.05，以上长度值已经无量纲化，噪声强度设为3。

3 数值模拟结果及对比 3.1 冰晶生长过程

在没有外部扰动的冷水结冰时，冰晶粒与增长较快的水平轴平行，所以冰晶在水平方向上生长更快，由极小的冰球体发展到薄圆冰盘，冰盘会继续发展成六角星形^[6]，从圆盘到六角星形意味着冰快速生长的开始，六角冰晶彼此冻结在一起后会形成冰壳，厚度逐渐增加后便形成了一定体积的冰，完成整个结冰过程。冰晶生长过程见图 2。

过冷水滴冻结过程通常分为再辉和主冻结2个阶段，再辉阶段主要表现为枝晶生长，原则上对水滴冻结再辉过程模拟即为枝晶生长模拟，由图 3可以观察初始时水滴与壁面接触点形成晶核，晶核迅速生长成枝晶形状而布满整个水滴，形成糊状^[25]。对枝晶生长的准确模拟有助于理解常被忽略的水滴冻结再辉过程。

陶乐仁利用低温显微镜观察了冰晶由晶核生长为类似雪花形状的过程^[26]，实验结果见图 4。本文在给定上述基本参数基础上，分别设置时间迭代步数n_t为5、50、200和2 000。由图 5可见，随着时间迭代步数的增加，初始晶核开始长大，首先由圆形生长为六边形，接着六边形的6个角开始突出生长，进而演化为6个分支，至此形成了六角冰晶的雏形。模拟结果与实验观察结果相一致。

3.2 模拟冰晶与自然界雪晶对比

加州理工学院Kenneth G. Libbrecht教授在《Field Guide to Snowflakes》一书中展现了其团队拍摄到的形式多样的雪花照片^[27]，呈现了大自然界中雪花(见图 6(d))。本文模型可以模拟雪花的形貌，保持上述其他参数不变，划分网格区域为800×800，无量纲时间为0.65时，由图 6可见，本文模拟的六角冰晶形貌与大自然中真实存在的雪花是相似的，两图中各主支均生长出侧分支，甚至二次分支，有效说明了本文相场模型的真实性。本文结果与雪花的对比只是象征意义上的，相场模型是研究过冷熔体结晶过程，而雪花是水蒸气凝华而成，因此相场模型并不描述雪花形成规律，可作研究参考。

3.3 尖端特征数据的获取

为了更好地描述枝晶生长行为，凝固学提出了一些枝晶生长参数，如尖端生长速度、半径、温度等。本文基于每一时间步的模拟结果，采用数学方法获得上述数据，以便定量分析枝晶生长规律。

定义ϕ(x, y, t)=0.5为固液界面，通过线性插值获得t时刻的尖端位置(x_tip, y_tip)，跟踪每个时间步的尖端位置，即可获得尖端速度和尖端温度，而尖端半径的获取则相对复杂，本文采用拉格朗日方法进行多项式插值，求得尖端半径^[20]：

$R_{\mathrm{tip}}=\frac{\partial \phi / \partial x_{(\mathrm{tip}, 0)}}{\partial^{2} \phi / \partial y_{(\mathrm{tip}, 0)}^{2}} $

(12)

如图 7所示为特征数据随无量纲时间变化曲线，主要计算参数为：γ=0.045，j=4，ε=0.005，Δ=0.6，Δx=Δy=0.01，Δt=0.000 1，网格数300×300。

图 7显示生长过程较稳定，尖端温度接近融点，尖端速度和半径在较小范围内波动，几乎恒定，该现象与Ivantsov枝晶生长理论保持一致，该理论认为枝晶稳定生长时其尖端以恒定速度向过冷熔体生长并且尖端温度等于融点^[28]。而图中初期尖端数据突变是为了迅速过渡至稳定期，表现为尖端温度迅速上升，而尖端半径和速度迅速下降，尖端温度上升是由于界面处潜热的释放，这与过冷水滴冻结第一阶段温度迅速升高现象一致^[29-30]；末期突变是由于枝晶尖端生长至边界；凝固中期生长行为较为稳定，本文将中期的尖端半径、尖端速度、尖端温度分别记为R、V、U，并称为枝晶尖端特征数据。

4 参数敏感性分析

Wheeler模型中无量纲过冷度、各向异性强度、噪声强度等参数会影响冰晶生长(包括尖端速度、尖端半径等)和形貌(侧枝大小、取向等)，本文将对模型主要参数进行敏感性分析，与真实冰晶形貌进行对比，进而提高冰晶生长模拟的准确性。

4.1 无量纲过冷度Δ

过冷度对冰晶的形核、生长过程及最终形貌均有重要影响，本例保持其他参数不变，分别取Δ为0.50、0.60、0.70，对应模拟及分析结果见图 8~9，由图 8可见，随着Δ值增大，枝晶尖端速度变大，尖端温度减小，这是因为在低Δ下，界面热扩散层厚度比较大，会阻碍潜热释放，从而抑制了界面处噪声扰动，导致生长缓慢；而在高过冷度下，界面处温度梯度比较大，热扩散层的厚度小，如此便有利于潜热的释放，从而加快生长。

图 9的模拟结果进一步验证了上述结论：过冷度较小时，生成的枝晶形貌小，枝晶主干较粗，随着过冷度增大，枝晶形貌逐步变大，主干变得更加细长，开始出现二次分枝，当过冷度足够大时甚至出现3次、4次分枝，这符合凝固理论中晶体连续生长机理。

4.2 各向异性强度γ

各向异性强度γ表示界面表面张力、界面厚度及界面动力学各向异性的程度，γ对枝晶生长过程及结果有重要影响。本例保持其他参数不变，分别取γ为0.01、0.02、0.03、0.055、0.065、0.07，对应模拟及分析结果见图 10、11，由图 10可见，随着γ的增大，尖端速度变大，而尖端半径及温度不断减小，这与微观可解性(MSC)理论关于枝晶尖端稳态行为结论一致^[28]。由图 11可见，γ较小时枝晶生长各向异性不明显，各分枝都比较粗大，主支生长优势不明显，趋于各向同性。随着γ增大，晶体沿主支方向生长优势开始凸显，主支逐渐变得细长，γ=0.05开始出现侧向分枝，这是由于γ增大时会使噪声扰动更容易被放大，生长前沿变得更不稳定，有益于枝晶分裂。但同时，随着γ更大，大于0.06时有两大主支被掩盖，仅四主支生长凸显，说明γ值不是越大越好，而是需取折中合理值。

4.3 噪声强度A_r

结冰过程扰动是枝晶侧向分支产生的主要原因，本例保持其他参数不变，分别取A_r为0.0、0.5、1.0、3.0对应结果见图 12~13，图 12给出了不同A_r下尖端速度及温度变化情况，可见随着A_r的增大，尖端数据出现波动现象，但仍在某一稳定值附近波动。由图 13可见，随着A_r增大，枝晶界面由光滑变为出现侧枝，因此A_r能促进生长竞争，进而形成具有发达分枝的枝晶形态^[20]。

5 结论

1) 叙述了结冰原理，解释了晶核—圆盘—六角冰晶—冰壳—厚冰过程，并通过引入晶核方式模拟了海水冰晶的生长，与显微镜实验观察结果及自然界真实存在的雪花进行了对比，验证了本文相场模型的真实性；

2) 借助Python语言中高效可靠的第三方科学计算库对控制方程求解过程进行了优化，有效提高了计算效率和准确性，并从形核生长位置、对称区域等方面节省了大量计算时间；

3) 对无量纲过冷度Δ、各向异性强度γ、噪声强度A_r等参数对枝晶尖端生长的影响进行了定量分析。结果表明，过冷度是海水结冰的驱动力，过冷度直接影响界面热扩散层厚度及潜热释放，Δ越大，冰晶界面由固相向液相生长越快，且能促进侧枝生长，形成分枝发达的枝晶形貌；各向异性强度对枝晶生长行为有重要影响，γ越大，枝晶生长界面处扰动更容易被放大，导致界面前沿不稳定，促进了侧枝生长，枝晶形貌变得更复杂，但γ增大到一定程度后，冰晶会不再对称，形貌开始失真，因此各向异性强度γ应选取适中值；噪声强度A_r的大小表征生长界面处扰动强弱，A_r越大，枝晶界面生长越激烈，越容易出现侧枝，但A_r过大时也会导致界面失稳失真，不宜取过大值；

4) 本文模拟结果与传统枝晶生长结论一致，模拟过程中对不用过冷度下枝晶生长速度进行了预测，并且考虑了界面动力学效应等因素，能为过冷水滴冻结再辉阶段模拟提供依据，进而完善水滴尺度结冰过程研究。

本文模拟了海水六角冰晶的生长，揭示了冰晶形成机理。后续研究将进一步建立相场模拟枝晶生长介观尺度与水滴冻结细观尺度的联系，建立结冰规律，开展多尺度的数理建模与预测手段研究，进而为极地船舶防/除冰设计及安全提供理论基础。