基于改进相场法的盐度对海水结冰微观特性影响分析

白 旭，杨苏杰，田于逵

（1.江苏科技大学船舶与海洋工程学院，江苏 镇江 212003；
2.中国船舶科学研究中心，江苏 无锡 214082）

海水飞沫是导致船舶结冰的主要原因，96%的船舶结冰与海水飞沫有关[1]。船舶表面结冰会降低船舶的稳性和降低结构稳定性，影响船上设备的使用，对船舶的安全和性能造成了极大影响。研究结冰特性对建立有效的船舶结冰预报方法具有重要意义[2]。

现有的结冰研究多为根据当地气候条件进行的结冰预报，或是利用CFD 软件对整体或构件进行结冰模拟，使用的公式也是根据已有数据推算的经验公式[3]。针对船舶的结冰模型有许多缺陷，如忽略了水与冰的换热、表面粗糙度的变化等[4]。分析及试验表明，冰的生长规律和尖端形貌会影响冰形的预测可靠性。由于成核率和稳定性的限制，所以冰晶生长的试验数据准确性差，而采用数值模拟冰晶生长的研究较少，且均有局限性。

为了掌握结冰的演变机理，开发船舶结冰特性的精细预测方法具有重要意义[5]。基于金兹堡-朗道理论的相场法是模拟凝固组织的一种有效方法，它能够描述复杂的晶体形貌。该方法已成为模拟微观组织演化的一种重要的、通用性极强的方法，成为模拟凝固过程中微观组织演化的有力工具[6]。

1998 年，Kobayashi[7]首先提出了扩散边界概念，并在此基础上建立了第一个相场模型，进行了金属的凝固模拟；
Wheeler 等[8]在模拟中将相场参数和实际物理量相结合，建立了WBM 模型；
Karma 等[9]建立了KKS模型，并消除了界面厚度的影响，扩大了模拟的过冷度范围，加快了计算速度；
Valeria等[10]提出了一个模型用来计算海水结冰过程中的盐分流通现象，描述了高盐度溶液与低盐度冰相分离的时间演化过程，但是没有分析盐度的影响。

国内学者对于相场模型的研究也有一定进展。张玉妥等[11]率先使用相场法模拟纯物质金属的凝固现象和等轴枝晶的生长，通过数值计算显示了等轴枝晶的形态；
邓儒超[6]进行了海水冰晶生长的相场模拟，发现过冷度和各向异性越大，冰晶的生长速度越大；
陈尚海[5]分析了对流情况下的海水冰晶生长凝固过程，枝晶尖端生长速度随着流速、扰动强度和过冷度的增大而增大；
唐杨新[12]给出了适用于海冰演化的多相相场模型的简单推导，应用Alber-Zhu的相场模型来描述不同种类的海冰之间界面的演化，介绍了海冰增长的相场模型和不同海冰在构型力驱动下相变演化模型的新证明方法。

本文基于相场法，建立船舶结冰微观结构特征预测的数值模拟方法，分析不同盐度条件下结冰的微观结构演化特征。海水飞沫引起的结冰过程是一个多晶核生长的复杂竞争过程，它是由振动触发的非均匀形核开始的，但盐度对成核率的影响无从得知，因此在分析时仅考虑单晶核条件下盐度对冰晶生长的影响。

1.1 纯物质相场模型

Wheeler 相场模型使用相场和温度场方程来实现对凝固现象的模拟[13]，对于纯物质的凝固，相场和温度场控制方程如下[8]：

相场方程为

式中：ε(θ)为考虑各向异性影响的相场参数，ε(θ)=εˉ( 1 +γcos(jθ))，γ为各向异性强度，j为各向异性模数，本文中取j= 6，εˉ为无量纲界面厚度；
m为相场迁移率；
α为界面自由能；
Δ为无量纲过冷度。

温度场方程为

式中：p"(ϕ)为固相分数p(ϕ)对ϕ的导数，p"(ϕ)=30ϕ2(1-ϕ2)；
u为无量纲温度，u=(T-TM)/(TM-T0)。

纯物质相场模型无法模拟海水的凝固现象，要进行海水的凝固模拟以及分析盐度的影响，需要引入溶质的影响。

1.2 考虑盐度的海水相场模型

研究盐度对海水凝固的影响时，将海水视作二元溶液。在海水凝固过程中，由于溶质盐在水中的溶解度大于在冰中的溶解度，溶质会在固液界面前沿富集，随着界面推进，固液界面将会向液相排出溶质，进而影响凝固过程。因此，在进行海水凝固模拟时，需要引入溶质再分配的影响。

基于菲克定律，引入溶质场方程，以实现凝固时计算域内的溶质变化：

式中，C取海水中盐的摩尔分数，D为无量纲扩散系数，Ds为固相扩散率，Dl为液相扩散率，k0为平衡分配系数。

计算域内的溶质浓度变化会导致熔点发生变化，因此引入过冷度的实时动态变化来实现盐度在凝固过程的影响。

式中，CP为海水比热容，ΔT为溶液过冷度，mL为液相线斜率，C0为溶液初始浓度，L为单位体积潜热。

添加溶质影响后的相场方程为

温度场方程为

2.1 参数取值

本文所取物性参数均按照浓度为3.50%的海水溶液选取，具体取值如表1所示[13]。

表1 海水溶液物性参数Tab.1 Physical properties of seawater solution

2.2 有限差分法及离散方式

对于相场方程，空间离散采用中心差分方法，时间离散采用显式差分方法。另外，对于方程中出现的拉普拉斯算子，采用九点差分方法进行离散：

计算时为了保证计算结果收敛，计算过程的稳定性由相场方程决定，时间步长需满足以下条件：

2.3 初始条件和边界条件

本文计算域设置为正方形，网格数设置为150×150，网格间距Δx和Δy在无量纲化后均为0.03，时间步长Δt= 0.0001 s。

具体晶核位置和条件设置如公式（11）所示：

相场、温度场和浓度场均采用Neumann边界条件，即

2.4 可靠性验证

对冰晶凝固过程进行模拟，并将其和实际自然条件下生成的冰晶进行对比。使用上文设置参数，网格划分为500×500, 模拟时间为7500Δt，得到结果见图1。通过比较模拟的冰晶形状和现实的冰晶形状[16]，发现实验结果吻合良好，可以进行进一步模拟计算。

世界海水盐度分布大致为3%～4%[14]。盐度对海水凝固的影响，除了会实时影响凝固场内的过冷度，还会直接影响海水的熔点。

取6 个不同的过冷度作为温度条件，在每个温度下，盐度从3.00%到3.90%，每隔0.1%取一个值，计算不同盐度的冰晶生长行为，分析盐度对冰晶生长过程的影响。

3.1 盐度与温度关系

取初始温度T0为240.85 K、243.85 K、246.85 K、249.85 K、252.85 K、255.85 K，在各温度条件下引入盐度引起的熔点变化。

根据Bodnar[15]的研究，用熔点相比零摄氏度的下降温度来描述盐度对熔点的影响，文中计算的盐度条件对应的熔点温度下降如表2所示。

表2 盐度对应的熔点下降Tab.2 Salinities corresponding to melting point depressions

三种不同温度条件下，不同盐度的过冷度参数如表3所示。

表3 初始温度参数Tab.3 Initial temperature parameters

图2 显示了温度为240.85 K 时的结冰图像，在低过冷度条件下，冰晶以六角冰晶的形式生长，但是并未生成分枝。图3 为枝晶尖端速度图像，尖端速度尽管略有波动，但还是逐渐过渡到平稳期，这与Ivantsov枝晶生长理论一致。

3.2 盐度对结冰过程的影响

表4 列举了三个温度条件下枝晶尖端速度的平均值。总体来说，枝晶生长速度随着温度降低而增加。盐度影响下的溶质场如图4 所示。固相中的浓度要低于液相中的浓度，这在固液界面处尤为明显。这是由于溶质在固相中的溶解度低于其在液相中的溶解度，随着固液界面的推进，溶质从固相中被排出并富集在固液界面处，与凝固理论相一致。

表4 各温度尖端速度均值Tab.4 Mean value of tip velocity at different temperatures

图5 展示了在不同温度下，盐度的变化对尖端速度产生的影响，具体表现为：随着盐度的增加，枝晶生长速度降低。这是因为盐度的增长导致了溶液熔点的降低，在外界温度不变的情况下，从而使溶液的过冷度减小，减缓了枝晶的生长。

盐度变化会影响枝晶的生长速度，从图5 可以看出，温度为240.85 K 时，盐度从3.00%增长到3.90%，尖端速度从5.280Δx/s变为4.271Δx/s，降低了19.1%，相比整体速度降低比较明显。

温度在240～246 K 之间时，由于过冷度较大，结冰速度较快，速度变化量明显，分别为1.0085Δx/s、1.0985Δx/s、0.9039Δx/s，而当温度在249～255 K 时，尖端速度的数值较小，此时盐度引起的速度变化量也很小，分别为0.0059Δx/s、0.0041Δx/s、0.0038Δx/s。

温度从240.85 K 升至246 K 时，尖端速度随盐度变化比例分别为19.1%、33.0%和76.9%，比例逐渐上升，这是因为随着温度升高，过冷度降低，盐度引起的过冷度变化占原有过冷度的比例越来越大，导致速度变化越来越明显。而当温度为249.85K 时，此时过冷度进一步降低，尖端速度发生数量级的变化，因此盐度引起的速度改变量与之前相比也显著减少，减少比例分别为27.3%、29.1%和46.3%。

可以看出，结冰速度随盐度变化而变化，且改变量显著，如果以海洋结构物作为考虑对象，那么不管是结冰尺寸还是结冰时间都会显著增大，带来的结冰量变化会更加明显。

为了提高船舶结构结冰特性的预测精度，本文利用Wheeler 相场模型，并用有限差分法对不同温度和盐度下过冷液体晶体生长过程进行了数值模拟，提取枝晶尖端速度作为衡量盐度对冰晶生长影响的标准，得到如下结论：

（1）盐度对海水结冰的影响主要来自于盐度与熔点之间的联系。盐度越大，海水熔点越低，在相同温度下导致过冷度越小。因此盐度越大，海水结冰越慢；

（2）在凝固过程中，溶质从固相析出并被排到液相，被排出的溶质在固液界面前富集，进一步影响了海水的凝固；

（3）盐度对冰晶生长的影响与过冷度有关。过冷度在240～246 K 之间时，结冰尖端速度较大，盐度引起的速度变化量分别为1.0085Δx/s、1.0985Δx/s、0.9039Δx/s，而当过冷度在249～255 K 时，此时尖端速度数值较小，盐度引起的尖端速度变化量随之减小，分别为0.0059Δx/s、0.0041Δx/s、0.0038Δx/s。

本文模拟了盐度对微观尺度海水结冰过程的影响。传统的针对船舶结构结冰的研究尚有许多不足之处，如《结构大气结冰标准》[17]中的结冰计算主要依据经验公式，结冰过程中的许多机理未得到解释。通过进行微观尺度的结冰研究，可以对结冰过程中的机理进行分析。下一步将着重于研究液滴撞击前后的微观结冰现象，探究撞击过程中的能量转化，并将其运用于船舶结构宏观结冰研究中。