基于LBM-CA 的六角冰晶非平衡生长模拟研究

蓝碧兰 郭 晴 全晓军

(上海交通大学机械与动力工程学院 上海 200240)

1 引言

结冰现象在自然界和工业生产中广泛存在,包括飞机机翼结冰、冷冻食品工艺及冰模板法制备多孔材料等。冰晶的生长分为非平衡阶段及平衡阶段。六角枝晶也称枝晶,是非平衡凝固过程中最常见的一种晶体生长形式。枝晶的生长提供平衡阶段的初始生长条件,深入探究枝晶生长的机理有利于加强对冰晶形貌的调控。Zhang 等[1]基于冰模板法制备了具有可调孔径和结构的纳米纤维素气凝胶,通过改变冷冻温度、温度梯度证明了冷冻温度和温度梯度会对冰晶的形貌和生长速率产生影响。Fukasawa 等[2]发现宏观温度梯度场会改变冰晶的形貌。研究表明冷冻工艺条件对冰晶的尺寸、形貌具有较大影响,但实验操作过程中往往由于结晶的随机性以及生长时间较短等因素,导致直接对冰晶进行定量观测和微观机理研究存在局限性。

为了进一步探究冰晶生长的微观机理,打破实验探究的局限性,一系列数值模拟方法被广泛用于结晶过程的研究,包括蒙特卡洛法、顶点法、相场法,元胞自动机法等。蒙特卡洛法简洁方便、易于编程,但存在模拟时间与实际时间无法对应的问题。顶点法高效性只体现于顶点的计算,使用受限。相场法(Phase Field,简称PF)和元胞自动机法是目前在晶体生长过程中应用较为广泛的两种方法。韩端锋等[3]采用经典Kobayashi 相场模型讨论了界面宽度均值、融化温度等多个参数对冰晶生长及形貌的影响,结果表明冰晶形貌对参数的依赖性较高。白旭等[4]采用Wheeler 相场模型,分析不同过冷度对海水冰晶生长的影响,结果表明冰晶生长速率随过冷度增加而增加。然而该模拟受计算精度及收敛性要求的限制,网格数量需求较多,对算力要求相对较高。相场法普遍存在界面厚度较难确定、结果对参数较为敏感的局限性。并且相场法计算量较大,而当前计算机数据处理能力有限,因此相场法仅适用于微小区域的模拟。

与相场法相比,元胞自动机(cellular automaton,CA)方法计算速度较快,能够适应较大范围内的凝固模拟。LBM(Lattice Boltzmann Method,LBM)是一种区别于传统传输模型的模拟手段。一方面LBM 避免了求解非线性方程组,简化了计算过程;另一方面,LBM 方法具有先天并行性、物理意义清晰的优点。耦合CA 方法与LBM 方法的模型具备并行性高、计算速度快、计算范围广的优点。周靖超等[5]基于LBM-CA 方法模拟强制对流作用下单相合金的非对称生长行为,模拟结果表明对流速度越大、初始过冷度越小、冷却速率越小会使枝晶非对称生长行为趋向更强。张勇等[6]基于LBM 方法研究了热障陶瓷涂层凝固生长过程,等离子喷涂涂层在正温度梯度方向呈柱状晶生长,温度梯度越大,柱状晶越发达。上述单相合金固溶体及等离子喷涂涂层晶体均呈四角晶形貌,而过冷水结冰多为六角冰晶。与四角晶相比,六角冰晶存在网格离散导致的各向异性的问题,是模拟中需要重点解决的问题。Guo 等[7]基于LBM-CA 方法研究自然对流及强制对流对冰晶的生长速率及形貌的影响。模型引入减弱因子削弱了各向异性,但该文章并未研究温度梯度及冷却速率等冷冻工艺对冰晶生长的影响。

目前采用LBM-CA 方法模拟冰晶生长的研究仍处于起步阶段。鲜少有学者研究温度梯度与冷却速率对枝晶生长的影响。因此,本文基于GPU 加速,建立LBM-CA 模型,探究过冷度、温度梯度及冷却速率等因素对枝晶形貌及生长速率的影响,并引入减弱因子和优化网格差分方法以削弱网格各向异性及离散各向异性对结果的影响[7]。

2 数值模型

2.1 CA 模型

过冷度ΔT是结晶过程的动力,主要包括热过冷(ΔTt)、曲率过冷(ΔTr)及成分过冷(ΔTc),本文是纯水结晶过程的模拟,因此ΔTc=0。

式中:Γ0为平均吉布斯-汤姆逊系数,εΓ为由各向异性界面张力引起的各向异性系数,K(x)为固液界面处的界面曲率,即:

不同的网格格点可以根据固相率分为固相元胞(S)、界面元胞(S/L)、液相元胞(L)。

对于尚未完全凝固的界面元胞,固相率将继续增加直至为1,成为固相网格,并且按照元胞自动机Moore 类型进行界面元胞的捕获。

界面处的法向角可以用式(6)计算得到:

根据尖锐界面模型,枝晶尖端生长速率是过冷度(ΔT)的函数:

式中:μk(θ)为各向异性动力学系数,可以通过式(8)计算:

式中:εμ为关于生长速率波动的函数,θp为优先生长方向。

2.2 格子玻尔兹曼模型(LBM)

温度场碰撞迁移格式为:

式中:τΓ为温度场的松弛函数,τΓ=3α+0.5,α为热扩散系数;φ(x)为温度场的源项,主要受凝固过程放热的影响:

温度场的平衡分布函数表达式为:

2.3 各向异性缩减处理

2.3.1 网格各向异性

为了降低网格各向异性,引入一个假想存在的各向同性标量场[8],在该标量场中标量的定义如下:

扩散方程为:

式中:D=0.24。

引入缩减因子bred(x):

经过比较,使用虚拟各向同性场ϕ(x)能够在界面处得到一个更为准确且稳定的法向角,因此法向角的计算公式调整为式(15):

2.3.2 离散各向异性

离散各向异性对界面曲率的大小有较大的影响,通过调整固相元胞周围8 个相邻网格的权重以及添加缩减因子可以使得界面曲率的取值进一步贴近实际情况。其中:

同时在界面曲率K(x)的分母项中添加各向异性削减因子,表示如下:

2.4 CA 模型与LBM 模型的耦合

CA 模型主要进行凝固过程冰晶生长过程的模拟,而LBM 模型主要进行温度场的模拟。通过CA模型计算的尖端生长速度可以计算出每一次迭代的固相率增量:

为了保证LBM 模型与CA 模型的相关性,LBM模型中无量纲扩散系数及动力粘度的选取如式(19):

式中:LBM 中δx及δt均为1,CA 模型中Δt及Δx均为7.02 ×10-7。

上述各式中部分常数取值见表1。

表1 水的特性参数Table 1 Physical parameters of water

3 过冷度对枝晶生长的影响

3.1 过冷度对枝晶生长形貌的影响

过冷度是凝固的主要驱动力,在枝晶生长中起重要作用。不同初始过冷度下,枝晶形貌显著不同。如图2 所示,枝晶的尖锐程度和枝状分化程度随过冷度增加而增加,主要原因是初始过冷度增大会导致生长驱动力增大、热扩散层减薄,从而提高枝晶生长速率。如图3 所示,低过冷度下枝晶尖端被更厚的热扩散层包围,热量扩散受阻碍。高过冷度下表现为界面层减薄,热量扩散速率增加,枝晶分化更明显。

图1 Moore 型格点捕获规则示意图Fig.1 Schematic of the Moore type capture rules

图2 t=1.20 s 不同过冷度条件下枝晶生长形貌图Fig.2 Dendrite growth morphology under different subcooling degree in same time

图3 不同过冷度枝晶温度云图对比图Fig.3 Contrast diagram of dendrite temperature cloud of different undercooling degrees

如图4a 所示,过冷度越大,枝晶臂越发达,尖端越深入过冷水。图4b 显示枝晶臂间出现不同程度的“缩颈”效应。过冷度越大,“缩颈”效应越明显。主要原因是过冷度越大,尖端生长越快,从而潜热释放加快。但枝晶臂间空间狭小,热量扩散困难。因此过冷度越大,臂间热积聚效应越明显,尖端生长受更大程度地抑制。

图4 t=1.26 s 不同过冷度枝晶形貌对比图Fig.4 Contrast diagram of dendrite morphology at different undercooling degrees in same time

3.2 过冷度对枝晶生长速度的影响

Ivantsov[9]基于热量扩散的假设给出了过冷度与枝晶生长Peclet 数(Pe)的关系:

式中:E1(Pe)表达式为:

Pe表示对流速率与扩散速率之比,表达式为:

Pe既是尖端速率也是尖端半径的函数,仅由式(21)无法求解。

引入Langer 和Muller-Krumbhaar[10-11]提出的临界稳定性理论。LMK 临界稳定性原理指出枝晶尖端稳定性取决于一个无量纲数σ:

定义σ∗为由实验及数值计算确定的稳定性因子,σ∗通常取0.02 或0.025。当σ>σ∗时,将会产生侧向分枝不稳定性使σ减小;当σ>σ∗时,又将导致尖端不稳定性使σ增大。因此,σ趋向于稳定在σ∗附近。因此,根据式(21)—式(24)即可求得生长速率随过冷度的变化曲线。

如图5 所示,初始阶段枝晶生长速率较大,而后生长速率逐渐降低直至趋于稳定。主要原因是枝晶生长初始阶段,枝晶尖端仍保持较高的过冷度,生长速率较大。而后随着界面元胞不断释放潜热,界面前沿温度不断上升直至潜热释放与热扩散达到动态平衡时,温度基本不变,生长速率趋于稳定。此外,图5对比了过冷度为0.65 K 时LMK 理论及本文模拟结果的生长速率变化。LMK 理论中当σ∗给定,枝晶尖端始终保持相同速度生长。而本文模拟结果除在初始阶段受初始过冷度影响而高于LMK 曲线外,整体与LMK 曲线较为一致。

图5 不同过冷度条件下尖端生长速率随时间变化图Fig.5 Diagram of tip growth rate changing with time under different degree of undercooling

如图6 所示,本模拟结果尖端生长速率与初始过冷度呈近似线性关系,主要原因是当Pe较小时,对流换热量与热传导相比处于较低水平,相变释放的潜热主要依靠热传导向周围传递。过冷水与枝晶取相同的导热系数,热传导各向同性。因此,界面元胞过冷度主要受初始过冷度影响。由式(7)可知,当尖端μk(θ)变化较小时,生长速率与初始过冷度呈近似直线关系。

图6 t=0.84 s 时尖端生长速率随过冷度变化关系图Fig.6 Growth rate of tip varies with degree of undercooling at t=0.84 s

令LMK 理论中σ∗分别为0.02 及0.025 计算出不同过冷度下枝晶尖端的生长速率并与本模拟结果进行对比。如图6 所示,二者在小过冷度范围内生长速率均随过冷度增加而增加,趋势一致。当过冷度过大或过小时,受网格尺寸及算力的限制理论分析结果与模拟结果出现一定的偏差。总体上,本模拟结果与LMK 理论分析结果较为吻合。

4 温度梯度

4.1 温度梯度对枝晶形貌的影响

在计算域范围内施加沿X轴负方向的温度梯度。不同的温度梯度下枝晶呈现显著不同的形貌。如图7 所示,负温度梯度方向的枝晶臂均受到不同程度的促进,其中与温度梯度反向平行方向的枝晶臂受到最大程度的促进。这是因为与温度梯度反向平行的方向是温度下降最快的方向,尖端过冷度最大,促进作用最显著。其余两个枝晶臂因为与X轴正方向存在一定夹角,所以促进效应减弱。然而,正温度梯度方向枝晶臂的生长受到了不同程度的抑制,其中与温度梯度平行方向枝晶的生长受到最大程度的抑制。主要原因是与温度梯度平行的方向是温度升高最快的方向,过冷度最小,抑制作用最显著。但由于周围流体处于过冷状态,所以枝晶仍具备生长动力。其余两个枝晶臂因为与X轴负方向存在一定夹角,抑制作用减弱。因此,在宏观温度梯度的作用下,枝晶形貌得到控制,枝晶实现了定向凝固。

图7 t=1.24 s 不同温度梯作用下枝晶形貌图Fig.7 Dendrite morphology under different temperature gradient at t=0.84 s

从图8 可看到,温度梯度越大,正温度梯度方向受抑制越严重,负温度梯度方向受促进越显著。主要原因是温度梯度越大,正温度梯度方向温度上升越快,生长速率越小,抑制作用越严重;而负温度梯度方向温度下降越快,生长速率越大,促进作用越显著。因此,合理选择温度梯度是对冰晶形貌进行调控的有效方法。

图8 t=0.84 s 不同温度梯度作用下枝晶形貌对比图Fig.8 Contrast diagram of dendrite morphology under different temperature gradient at t=0.84 s

4.2 温度梯度对枝晶尖端生长速率的影响

如图9 所示,负温度梯度方向枝晶臂生长速率不断上升,而正温度梯度方向枝晶臂生长速率不断下降。该现象与图5 等温枝晶后期生长速率趋于稳定的现象不同。主要原因是负温度梯度方向枝晶尖端可以不断深入温度更低的过冷水中,相变产生的热量能够迅速地向周围过冷液体扩散。相反地,正温度梯度方向枝晶尖端周围过冷水温度提升,凝固推动力减弱,尖端生长速率不断下降。

图9 不同温度梯度下枝晶尖端生长速率随时间变化图Fig.9 Diagram of tip growth rate changing with time under different temperature gradient

图10 是t=0.82 s 时尖端生长速率随温度梯度的变化图。在负温度梯度方向枝晶尖端生长速率随温度梯度的增大而增大;而正温度梯度方向,枝晶尖端生长速率随温度梯度的增加反而下降。原因是在温度梯度作用下枝晶生长偏离对称等温生长,促进作用与抑制作用均随温度梯度的增大而增大。

图10 t=0.82 s 时的尖端生长速率随着温度梯度的变化图Fig.10 Growth rate of tip varies with t temperature gradient at t=0.84 s

值得注意的是,此时生长速率与温度梯度并不为线性关系。两条曲线的斜率均随温度梯度的增大而增大。主要原因是温度梯度的作用具有正反馈特性。温度场的差异导致生长速率的差异,而生长速率的差异又通过影响界面所处位置而作用于生长速率本身。例如负温度梯度方向,温度梯度越大,枝晶生长速率越大,尖端越深入低温过冷水,生长速率将进一步提升。

5 冷却速率

如图11 所示,恒定冷却速率作用下枝晶表现出高度对称性。因为施加恒定冷却速率后冰晶尖端过冷度仍保持相同,6 个枝晶臂的生长速率一致,枝晶对称生长。值得注意的是,枝晶的尺寸随冷却速率的增加而增大。冷却速率越大,枝晶生长速率越快,相同时刻枝晶越大;而冷却速率越小,枝晶生长速率越小,相同时间内枝晶将越小。因此可以通过合理地调控冷却速率达到控制枝晶大小的目的。

图11 不同冷却速率下枝晶生长形貌对比图Fig.11 Contrast diagram of dendrite morphology under different cooling rate

6 结论

基于格子玻尔兹曼方法及元胞自动机方法,利用尖锐界面模型,探究过冷度、温度梯度及冷却速率对枝晶形貌的影响,结论如下:

(1)基于GPU 加速,采用LBM-CA 算法,引入网格各向异性削弱因子,模拟了冰晶在过冷水中的生长,结果与实际情况吻合,各向异性削弱因子表现良好。

(2)过冷度影响枝晶生长速率从而影响枝晶形貌及大小,枝晶尖端生长速率随过冷度的增大而增大,呈现近似直线的关系。模拟结果与LMK 临界稳定性理论吻合,验证了模拟结果的正确性。

(3)通过施加不同温度梯度实现冰晶的定向凝固。与温度梯度相反的方向,枝晶生长受到促进,其中与温度反向平行方向受到最大程度的促进。与温度梯度相同的方向,枝晶生长受到抑制,其中平行方向受到最大程度的抑制。

(4)冷却速率影响枝晶的尺寸。施加恒定的冷却速率可以得到对称的六角冰晶,但枝晶的大小与冷却速率呈现正相关关系。因此,合适的冷却速率是制备优异性能材料过程中需要重要考虑的因素。