过冷度影响海水结冰形状与速度的相场模拟

doi:10.16183/j.cnki.jsjtu.2020.077

过冷度影响海水结冰形状与速度的相场模拟

白旭^,, 杨苏杰

江苏科技大学,船舶与海洋工程学院,江苏镇江,212003

Simulation of Undercooling Influence on Shape and Velocity of Seawater Freezing Process Using Phase Field Method

BAI Xu^,, YANG Sujie

School of Naval Architecture and Ocean Engineering, Jiangsu University of Science and Technology, Zhenjiang 212003, Jiangsu, China

责任编辑: 陈晓燕

收稿日期: 2020-03-19

基金资助:

国家自然科学基金(51879125)
江苏省高等学校自然科学研究重大项目(18KJA580003)
江苏省“六大人才高峰”高层次人才项目(2018-KTHY-033)
江苏省研究生科研与实践创新计划项目(SJKY19_2664)

Received: 2020-03-19

作者简介 About authors

白旭(1984-),男,江苏省沛县人,副教授,主要从事船舶与海洋结构物安全性分析技术方面的研究.电话(Tel.):0511-84401133;E-mail:baixu@just.edu.cn. E-mail：baixu@just.edu.cn。

摘要

为了研究船舶结构结冰的微观机理,进行了模拟计算.基于Wheeler相场模型,使用有限差分法,对不同过冷度下海水凝固形状进行模拟.在模拟过程中,将海水视为盐和纯水的二元混合物,通过设置冰物理参数和引入晶核方式保证了冰晶生长的真实性.计算时,设置4个晶核模拟冰晶生长,讨论不同过冷度对冰晶生长的影响.结果表明:冰晶生长速度随过冷度增加而增加,且与过冷度变化呈线性关系.当无量纲过冷度小于0.85时,枝晶仅有少量分枝,且枝干较细.当无量纲过冷度大于0.85时,冰晶出现二级以上的多级分枝,且枝干明显变粗.当无量纲过冷度达到1.0时,分枝受到主枝挤压而减少,晶核之间最终被填满.

关键词： 相场法; 海水结冰; 冰晶生长; 有限差分法

Abstract

In order to study the microcosmic mechanism of ice formation in ship structure, a numerical simulation is conducted. Based on the Wheeler phase field model, the finite difference method is used to simulate the solidification shape of seawater at different undercoolings. In the process of simulation, sea water is regarded as the binary mixture of salt and pure water, and the authenticity of ice crystal growth is ensured by setting ice physical parameters and introducing crystal nucleus. In the calculation, four nuclei are set up to simulate the ice crystal growth, and the influence of different undercoolings on the ice crystal growth is discussed. The results show that the growth rate of ice crystal increases with the increase of undercooling, and it is linear with the change of undercooling. When the dimensionless undercooling is less than 0.85, the dendrite has only a few branches which are thin. When the dimensionless undercooling is higher than 0.85, the ice crystal has more than two-stages of branches which are obviously thicker. When the dimensionless undercooling reaches 1.0, the branches are squeezed by the main branches and lessned, and the crystal cores are finally filled.

Keywords： phase field method; seawater icing; ice crystal growth; finite difference approach

PDF (4210KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

白旭, 杨苏杰. 过冷度影响海水结冰形状与速度的相场模拟[J]. 上海交通大学学报（自然版）, 2021, 55(5): 513-520 doi:10.16183/j.cnki.jsjtu.2020.077

BAI Xu, YANG Sujie. Simulation of Undercooling Influence on Shape and Velocity of Seawater Freezing Process Using Phase Field Method[J]. Journal of shanghai Jiaotong University, 2021, 55(5): 513-520 doi:10.16183/j.cnki.jsjtu.2020.077

近年来,随着全球气候变暖,北极地区冰川消融,北极航道得到开发.但是,由于北极地区特殊的地理环境,船舶在北极航行时极易受到结冰的影响^[1].已有的针对船舶结冰的研究大多集中在宏观层面,即不同气候条件下,海洋结构物各部位产生的结冰以及结冰造成的影响,结冰的原理也主要来自于经验公式和观测结果总结,很少有对其产生原因和微观机理的分析^[2,3,4].通过研究海水结冰的微观机理,用实际环境参数来模拟海水结冰过程,可以得到更加真实的冰晶生长情况,为过冷条件下海水飞沫结冰提供微观理论基础.针对海水的微观结冰过程,现在国内相关研究主要集中在单个晶核的凝固现象,而在实际情况下,外界环境所导致的结冰几乎都是多晶核条件下的结冰.因此,本研究不仅会分析过冷度对海水结冰的影响,同时会涉及多晶核条件下的海水凝固情况.研究水滴凝固的微观机理的重点在于追踪冰晶凝固时固液界面的移动,而通过引入相场法,建立温度、相场等多个微分方程,可以将真实环境条件引入模拟计算,模拟固液扩散界面,从而得到冰晶的生成演化过程^[5,6].

相场法最早是在20世纪80年代由Collins等提出^[7].1993年,Kobayashi首先建立了关于相场和温度的相场模型,在考虑各向异性的情况下,模拟了金属凝固时的枝晶生长^[8].Wheeler等^[9]模拟了Ni的凝固过程,并提出了二元合金的相场方程,即WBM模型,并将相场模拟与实际物理量联系起来,定量分析参数对凝固现象的影响.Fan 等^[10]基于平均浓度场和热、质及动量输运现象,模拟了蛋白质溶液在冷冻过程中蛋白质的分离与重分布.Yang等^[11]使用冰模板,考虑石墨烯再分配和冰晶的各向异性,模拟了在水冻结时石墨烯气凝胶的形成过程.

国内学者对于相场模型的研究也有一定进展.杨燕等^[12]利用Wheeler模型模拟了过冷条件下不同过冷度对冰晶生长的影响.朱昌盛等^[13]基于KKS模型模拟了在等温和不等温两种条件下潜热对二元合金凝固的影响.于艳梅等^[14]对单相和二元合金分别进行相场模拟,分析空间尺度、外界噪声及各向异性等参数对过冷熔体凝固过程的影响.韩端锋等^[15,16]分别使用Kobayashi模型和Wheeler模型,模拟了单晶核条件下纯水和海水的凝固过程,定量分析了外界参数对冰晶凝固的影响.

针对船舶结构的宏观结冰与微观机理的研究方面,Ryerson等^[17]利用船舶切割下的冰晶进行研究,将结冰的宏观现象和微观机理相联系,得到了结构结冰的物理性质和结冰部位的微观分析.Singh等^[18]利用分子动力学定量研究了光滑粗糙表面上过冷水滴的行为,分析了在粗糙表面上凝固温度和液滴黏附功与粗糙度之间的联系.Sultana等^[19]利用数值模拟分析了单个盐水滴在撞击极冷表面后的相变情况,探讨了液相率、盐度等物理性质对液滴冻结的影响.

本文基于Wheeler相场模型,运用有限差分法,对不同过冷度下海水结冰情况进行模拟,研究过冷度对海水凝固时冰晶生长形状和速度的影响,为船舶结冰预报提供理论和技术支持.

1 Wheeler相场模型理论

1.1 相场方程

相场法以金兹堡-道朗方程为研究基础^[20],通过引入相场参数ϕ来表示固液界面的移动.当ϕ=1时表示固相,ϕ=0时表示液相,0<ϕ<1时表示固液界面,通过数值的变化来实现对固液界面的跟踪.

金兹堡-道朗自由能方程为

(1)f(ϕ, T)=W

\int_{0}^{ϕ}

p(ϕ)[1-p(ϕ)]×

[p (ϕ) - \frac{1}{2} - β (T)]

dp(ϕ)

式中:f(ϕ, T)和β(T)分别为自由能密度函数和热力学驱动力函数,β(T)<1/2;W为单位体积能量;p(ϕ)为固相分数,1-p(ϕ)为液相分数,p(ϕ)=ϕ³(10-15ϕ+6ϕ²).

相场方程为

(2)

\frac{{\bar{ε}}^{2}}{m} \frac{\partial ϕ}{\partial t}

=ϕ(1-ϕ)

[ϕ - \frac{1}{2} + 30 \bar{ε} αΔuϕ (1 - ϕ)]

\frac{\partial}{\partial x} (ε (θ) ε' (θ) \frac{\partial ϕ}{\partial y})

\frac{\partial}{\partial y} (ε (θ) ε' (θ) \frac{\partial ϕ}{\partial x})

+Δ·(ε²(θ)Δ ϕ)

式中:ε(θ)为考虑各向异性影响的相场参数,其值与晶体和主轴之间夹角有关,ε(θ)= $\bar{ε}$ (1+γcos(jθ)),γ为各向异性强度,j为各向异性模数,本文中取j=6, $\bar{ε}$ 为无量纲界面厚度;m为相场迁移率;t为时间;α为界面自由能;Δ为无量纲过冷度;u为无量纲温度.参数之间满足:

$\bar{ε}$ = $\frac{δ}{ω}$

m= $\frac{μσ T_{M}}{KL}$

α= $\frac{\sqrt[]{2} ω L^{2}}{12 c_{p} σ T_{M}}$

Δ= $\frac{c_{p} ΔT}{L}$

其中:δ为特征长度;ω为长度参考尺寸;μ为界面迁移率;σ为表面能;T_M为海水凝固点;K为热扩散系数;L为单位体积潜热;c_p为比热容;ΔT为溶液过冷度,ΔT=T_M-T₀,T₀为溶液初始温度.

温度场方程为

(3)

\frac{\partial u}{\partial t}

\frac{p' (ϕ)}{Δ} \frac{\partial ϕ}{\partial t}

=Δ²u

式中:p(ϕ)的导数p'(ϕ)=30ϕ²(1-ϕ)²;无量纲温度u=(T-T_M)/(T_M-T₀).

Wheeler相场将实际物理量与相场参数相结合,因此在计算中,应该使用实际物理参数进行计算.但是为了将界面厚度控制在可计算范围内,取ω为2.7×10^-4 cm,将其无量纲化,实际界面厚度等于无量纲界面厚度与长度参考尺寸的乘积,网格间距也按照同样方法进行换算.

1.2 溶质场方程

在模拟海水凝固时,将海水视为水和盐的二元溶液,相比于一元溶液,二元溶液需要考虑溶质的影响.溶质对凝固的影响主要有两种方式,一种是引入溶质对过冷度的影响,另一种是在方程中添加浓度项^[21].本文通过引入过冷度变化,来实现溶质扩散对溶液凝固的影响.溶质场方程为

(4)

\frac{\partial C}{\partial t}

=Δ·D

[Δ C + \frac{(k_{0} - 1) C}{ϕ + k_{0} (1 - ϕ)} Δ ϕ]

(5)D=D_s+(D_l-D_s)

\frac{ϕ}{ϕ + k_{0} (1 - ϕ)}

式中:C为海水中盐的摩尔分数,海水中盐的质量分数约为3.5%,则盐的摩尔分数为0.011;k₀为平衡分配系数;D为无量纲扩散系数;D_s为固相扩散率;D_l为液相扩散率.

模拟过程中,计算域内的过冷度受到溶质浓度影响,实时过冷度为

(6)Δ'=

\frac{C_{p} [Δ T + m_{L} (C - C_{0})]}{L}

式中:m_L为液相线斜率;C₀为溶液初始质量分数.

1.3 参数取值

本文所取物性参数均按照质量分数为3.5%的海水溶液选取^[15],具体取值如表1所示.

表1 海水溶液物性参数

Tab.1 Physical properties of seawater solution

物理量	数值
c_p/(J·K^-1·cm³)	4.212
μ/(m·K^-1·S^-1)	7.4
σ/(J·cm²)	7.654×10^-6
T_M/K	269.75
K/(cm²·s^-1)	1.31×10^-3
L/(J·cm³)	333.5
D_s/(cm²·s^-1)	2.5×10^-7
D_l/(cm²·s^-1)	5.6×10^-6
k₀	0.8
$\bar{ε}$	0.005
m	0.035
α	405.3065

新窗口打开| 下载CSV

2 数值求解与可靠性验证

在得到相场模型后,需要对其进行数值求解.有限差分法原理简单,易于编程,因此本文采用有限差分法来求解控制方程.有限差分法是将要求的函数在时间和空间上离散,用差分代替微分,并用离散的值来表示函数在计算域内的变化,结合初始条件和边界条件,求出函数在计算域内各网格节点的值.

2.1 有限差分法及离散方式

求解相场控制方程时,空间离散采用中心差分方法,时间离散采用显式差分方法.另外,对于方程中出现的拉普拉斯算子,采用九点差分方法进行离散:

$\frac{\partial ϕ}{\partial x}$ = $\frac{ϕ_{i + 1, j}^{n} - ϕ_{i - 1, j}^{n}}{2 Δ x}$

$\frac{\partial ϕ}{\partial y}$ = $\frac{ϕ_{i, j + 1}^{n} - ϕ_{i, j - 1}^{n}}{2 Δ y}$

$\frac{\partial ϕ}{\partial t}$ = $\frac{ϕ_{i, j}^{n + 1} - ϕ_{i, j}^{n}}{Δ t}$

Δ²ϕ=[2(ϕ_i_+1,_j+ϕ_i_-1,_j+ϕ_i_,_j₊₁+ϕ_i_,_j_-1)+ϕ_i_+1,_j₊₁+ϕ_i_-1,_j_-1+ϕ_i_-1,_j₊₁+ϕ_i_+1,_j_-1]/[3(Δx)²]

式中:Δx、Δy为网格间距;Δt为时间步长.对于温度场和溶质场方程,采用与相场方程同样的差分格式.

在计算域内,为了满足计算精度,保证计算结果收敛,需要同时考虑相场、温度场及溶质场.计算域内的时间步长和空间步长需满足以下条件

(7)Δt

[\frac{1}{{(Δ x)}^{2}} + \frac{1}{{(Δ y)}^{2}}]

≤

\frac{1}{2 m}

(8)Δt

[\frac{1}{{(Δ x)}^{2}} + \frac{1}{{(Δ y)}^{2}}]

≤

\frac{1}{2}

2.2 初始条件和边界条件

本文计算域设置为正方形,网格数设置为1 500×1 500,网格间距Δx和Δy在无量纲化后均为0.03,时间步长Δt=0.000 1 s.在计算域内,按照正方形设置4个初始晶核,晶核半径为R.初始晶核距离计算域边界均为50R,晶核之间距离为50R.晶核分布如图1所示.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 初始晶核大小与位置

Fig.1 Initial nuclei size and position

晶核外相场参数ϕ设为1,无量纲温度u为-Δ,质量分数C为C₀,具体晶核位置和条件设置如下式:

(x-500)²+(y-500)²≤R²时,

(9)ϕ=0, u=-Δ, C=C₀

(x-1000)²+(y-500)²≤R²时,

(10)ϕ=0, u=-Δ, C=C₀

(x-500)²+(y-1000)²≤R²时,

(11)ϕ=0, u=-Δ, C=C₀

(x-1000)²+(y-1000)²≤R²时,

(12)ϕ=0, u=-Δ, C=C₀

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

(13)

\frac{\partial ϕ}{\partial n} = \frac{\partial u}{\partial n} = \frac{\partial C}{\partial n} = 0

2.3 模拟冰晶与自然冰晶的对比

为了得到海水结冰的过程及规律,首先对冰晶凝固过程进行模拟,并将其和实际自然条件下生成的冰晶进行对比.使用上文设置参数,网格划分为500×500, 模拟时间为7500Δt,得到结果如图2所示.冰晶在凝固时,先由晶核发展为冰盘,随后主要沿着6个方向发展,形成主枝,主枝之间的夹角约为60°,主枝随后向两侧长出侧枝和多次分枝,最终形成六角冰晶.通过比较模拟的冰晶形状和现实的冰晶形状^[22],发现实验结果吻合良好,可以进行进一步模拟计算.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 模拟冰晶与实际冰晶形状对比

Fig.2 Comparison of simulated ice crystal and actual ice crystal

3 过冷度对晶核生长形状及速度的影响

3.1 过冷度对晶核生长形状的影响

根据凝固理论,过冷度是液体凝固的驱动力,在冰晶的形核和生长中起重要作用.本文通过分析冰晶形貌、生长速度以及枝晶间距来描述过冷度对冰晶生长的影响.

选择0.80、0.85、0.90、0.95及1.00这5个不同无量纲过冷度作为模拟温度,各无量纲过冷度对应实际过冷度分别为 -63.34、-67.30、-71.26、-75.22 及 -79.18 K.总模拟时间为 7500Δt,模拟结果如图3所示.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 冰晶凝固过程

Fig.3 Solidification process of ice crystals

从图中可以看出,当Δ=0.80时,枝晶生长缓慢,仅有少量二次分枝,枝晶直到模拟结束未发生交汇,仍然是晶核各自发展.随着无量纲过冷度的增大,枝晶生长速度加快,主枝变粗,分枝数量增加,但枝晶之间仍未发生交汇.当Δ=0.90时,生长速度进一步增加,出现二级以上多次分枝,枝晶在模拟过程中出现交汇.交汇后主枝停止生长,分枝未受影响,继续向外生长.Δ=1.00时,粗大的主枝占据了极大的空间,使得分枝数量减少.当主枝交汇后,主枝逐渐变粗,分枝继续发展,逐渐填满晶核间空隙,几乎形成冰盘.

图4为各过冷度下分别模拟 10000Δt得到结果.通过比较各过冷度凝固结果可以看出,在同时存在多个晶核时,过冷度越大,枝晶就越早接触,枝晶内部空隙就越小,以至于最终完全接触,内部不再存在空隙.过冷度较小时,横向与纵向主枝有粗细差别,横向主枝较粗,而过冷度较大时,横向与纵向主枝粗细基本相同.

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 10000Δt各无量纲过冷度凝固结果

Fig.4 Solidification results of each dimensionless undercooling at 10000Δt

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

为了定量分析过冷度对海水凝固的影响,使用枝晶尖端生长速度v_tip来描述冰晶生长情况,枝晶尖端生长速度与无量纲时间和无量纲过冷度之间关系如图5所示.

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 枝晶尖端速度变化情况

Fig.5 Changes of dendrite tip velocity

图5(a)为无量纲过冷度为0.8时,尖端速度随无量纲时间变化情况.根据计算结果,枝晶在生长过程中较为稳定,尖端速度在开始时较大,随后慢慢减小,近似呈正弦函数波动.枝晶尖端速度随无量纲过冷度的变化如图5(b)所示,在无量纲过冷度从0.8增至1.0的过程中,枝晶尖端平均速度分别为0.288、0.360、0.432、0.504、0.576 cm/s,枝晶生长速度随着无量纲过冷度的增加而明显增加.无量纲过冷度对冰晶凝结的主要影响在于,无量纲过冷度越大,界面厚度就越薄,凝固时散发的潜热就越多,从而使冰晶凝固速度更快.因此,代入实际物理量后可知,过冷度越低,冰晶凝固越快,且凝固速度与过冷度基本为线性关系.

3.3 过冷度对多晶核生长形状的影响

本文在计算域内设置4个晶核,在晶核初始间距相同的情况下,通过计算不同过冷度下各晶核主枝之间的距离变化来分析过冷度对海水结冰的影响情况.

不同过冷度下枝晶尖端距离D_tip随无量纲时间的变化如图6(a)所示.随着无量纲时间的增加,枝晶尖端间距逐渐减小,且二者呈线性关系.当无量纲过冷度达到0.90时,枝晶间距最终为0,代表不同晶核生长出的枝晶最终接触.根据图中趋势,无量纲过冷度越大,接触所需时间越短.图6(b)为将枝晶尖端距离随时间变化的曲线拟合为一次函数后的曲线斜率K_tip,其代表枝晶间距的变化快慢,即枝晶生长速度的大小.根据曲线所示,斜率与无量纲过冷度之间近似呈线性关系,说明了随着无量纲过冷度的增大,枝晶生长速度的变化也会相应地近似呈线性增大.

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 枝晶尖端距离及枝晶距离斜率变化

Fig.6 Changes of dendrite tip distance and dendrite distance slope

4 结论

基于Wheeler相场模型,使用有限差分法,将海水看作盐-水二元溶液,把实际物理量转换为方程中的无量纲数,提取出枝晶生长速度作为定量分析海水结冰的指标,并通过不同过冷度下各晶核枝晶之间的间距变化来描述过冷度对海水结冰的影响.具体结论如下:

(1) 无量纲过冷度为0.90以下时,枝晶较细,没有或仅有二级分枝,枝晶生长较慢.无量纲过冷度为0.90和0.95时,枝晶较粗,枝晶生长速度较快,出现二级以上分枝.无量纲过冷度为1.00时,主枝与分枝明显变粗,生长速度明显加快,分枝之间几乎无空隙.

(2) 无量纲过冷度为0.80时,枝晶生长慢,形状细,晶核内部空隙大,最终形成各主枝组合成的枝状冰晶.无量纲过冷度为0.90和0.95时,有多次分枝,分枝之间有一定空隙,且空隙随着分枝发展而减少.无量纲过冷度1.00时,枝晶生长速度快,分枝之间空隙小,晶核内部空隙被填满,形成冰盘.

(3) 在各过冷度下,枝晶生长速度较为稳定,尖端最大速度随无量纲时间呈正弦变化.无量纲过冷度从0.80增至1.00时,枝晶尖端生长速度从0.12 cm/s增至0.24 cm/s,生长速度与过冷度之间为线性关系.

(4) 枝晶尖端间距随着时间的增加而减小,且二者呈线性关系,间距曲线斜率与无量纲过冷度之间近似为线性关系,其反映了枝晶生长速度的与无量纲过冷度之间也近似为线性关系,与之前求得的枝晶尖端速度相吻合.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

陆煊, 崔玫, 曹洪波, 等.

船舶防冻除冰技术现状与发展

[J]. 船海工程, 2016, 45(2):37-39.