基于元胞自动机-格子玻尔兹曼模型的枝晶碰撞行为模拟*

张士杰 王颖明 王琦 李晨宇 李日

(河北工业大学材料科学与工程学院,天津 300401)

合金凝固过程中,游离枝晶在熔体中的运动行为是研究合金凝固组织形成过程的关键问题之一.元胞自动机-格子玻尔兹曼耦合模型是近年来进行凝固微观组织数值模拟的主要数值模型.本文改进了模拟枝晶生长的元胞自动机和格子玻尔兹曼模型,使其能够模拟过冷熔体中等轴晶的运动行为.改进的模型采用伽利略不变的动量交换法计算流体力,通过求解质心运动方程计算枝晶的运动位移,并通过动网格技术实现枝晶的运动,运用硬球模型处理枝晶的碰撞.采用该模型模拟了Al-4.7%Cu 合金过冷熔体中单枝晶的沉降、牛顿流体中两圆形粒子的沉降和两枝晶的弹性碰撞.模拟结果表明,本模型在计算枝晶生长运动过程时可以很好地维持枝晶的形貌.采用本模型计算等轴枝晶的碰撞过程表明,枝晶的运动会扰动其周围的熔体,造成周围熔体浓度显著变化,进而影响枝晶的生长,加剧枝晶生长的不对称性.

1 引言

合金凝固过程中,铸型壁面的表层细晶粒可能会脱落,已经凝固的枝晶臂(特别是二次枝晶臂)也可能发生熔断现象,再加上在熔体中随机出现的数量可观的自由等轴晶,所以在液相区会出现大量游离的等轴晶,这些游离的等轴晶在重力和对流的推动下在合金熔体中运动(平移和旋转),且常常会发生碰撞现象[1].熔体中游离枝晶的运动对铸件的组织形貌和成分分布会有很大的影响[2].

当前,对合金凝固过程中枝晶生长的数值模拟研究已经比较成熟,可以有效地计算出纯扩散和对流作用下枝晶的生长过程[3-5],并已发展到了三维(3 dimension,3D)[6-8].而关于枝晶运动的模拟,由于在计算枝晶运动的同时还要同时计算枝晶生长,所以有一定难度.目前该方面的研究比较少,已发表的枝晶运动的文献中大多采用相场(phase field,PF)方法[9-13].采用PF 计算出的界面是连续扩散界面,对运动界面的处理有天然的优势,但是计算量非常大,计算效率较低,计算规模小.而具有计算简单快速、天然并行等优点的元胞自动机(cellular automata,CA)模型在模拟微观组织方面的应用,大多数都只模拟了静止枝晶的生长[14-17],对于枝晶运动的模拟研究很少,且不能很好地维持枝晶的形貌[18,19].本文旨在改进CA 模型,使其在模拟枝晶运动时能够较好地维持枝晶的形貌和成分,在提升效率的同时提高模拟精度;同时,在模型中还加入了碰撞处理机制,使该模型能够处理枝晶运动引起的碰撞现象.

2 模拟方法

采用格子玻尔兹曼(lattice Boltzmann,LB)模型计算流场、温度场和浓度场,采用CA 方法模拟微观枝晶生长,通过求解运动方程来计算枝晶的运动位移(平动和转动),采用硬球模型计算枝晶碰撞(假设枝晶碰撞为完全弹性碰撞),耦合上述方法建立起枝晶运动生长的数值模型以研究枝晶在熔体中的运动行为.

2.1 LB 模型

采用单松弛时间的D2Q9 双分布函数模型[20],流场、温度场和浓度场的玻尔兹曼方程为

其中ωi是权重函数;c是格子速度;u是流体宏观速度,温度场和浓度场的格式和流场一样.流体的流场、温度场和浓度场的松弛时间τf,τt,τc分别和动力学黏度ν、热扩散系数α、溶质扩散系数D有关:

流场、温度场和浓度场的源项为

其中 Δfs是固相率增量;L是潜热;Cp是比热容;Cl是液相浓度;k是平衡分配系数;F是粒子受到的合力,由Boussinesq 近似[21]给出:

其中g为重力加速度;ρ0,T0,C0分别为初始密度、温度和浓度;βT和βC分别为温度膨胀系数和溶质膨胀系数.

流体宏观密度ρ,速度u,温度T和浓度C分别为

2.2 CA 模型

溶质扩散模型最早由Rappaz和Thévoz[22]于1987 年提出的,之后经过多次改进.2007 年,Zhu和Stefanescu[23]提出了用于定量化计算合金凝固过程中树枝晶生长的虚拟前沿跟踪模型.该模型认为合金枝晶生长是由局部液体平衡成分和局部液体实际成分的差异所驱动的.根据局部平衡的热力学概念,界面平衡温度T*定义为

fs为元胞的固相分数.根据Gibbs-Thomson公式,界面能各向异性函数可以表示为

其中δ为各向异性系数.因此,界面平衡浓度可以计算为

如果当前时刻界面元胞液相浓度小于界面平衡浓度,则该元胞固相增加.固相率增量为

2.3 枝晶运动方程

将运动枝晶看作刚体处理,枝晶在熔体中运动所受到的流体力采用伽利略不变的动量交换法[24]计算,流体对α枝晶的边界点x的作用力为

其中Ω是α枝晶的边界区域;ρl和ρs分别为流体和固体的密度,固体质量

枝晶的速度和角速度为

其中Iα为惯性张量.

2.4 碰撞模型

游离的枝晶在熔体中碰撞的情况很复杂,初步假定枝晶间的碰撞为完全弹性碰撞,因此采用硬球模型来处理碰撞步骤.将枝晶看做刚性粒子处理.对枝晶应用动量定理和动量矩定理:

其中m为枝晶质量;v0和v是碰撞前、后的速度;ω0和ω是碰撞前、后的角速度;S为内力冲量.

枝晶碰撞的瞬间,由于两枝晶碰撞所受内力冲量相等,根据冲量定理和冲量矩定理,可以求得

其中ω10,ω20和ω1,ω2分别为碰撞前、后枝晶1和枝晶2 的角速度;v10,v20和v1,v2分别为碰撞前、后枝晶1和枝晶2 的速度;I1,I2为枝晶1,2 的转动惯量;r1和r2分别为碰撞点与枝晶1和枝晶2 质心的位置矢量.

再引入恢复系数e,并与(26)式联立写成矩阵AX=b形式:

2.5 动网格方案

虽然CA 方法计算效率高,但是在处理移动边界问题上存在困难.由于CA 方法的界面是尖锐界面,所以其界面在空间上的移动是不连续的,然而浓度在空间上的传输必须保证是连续的,这导致了在处理枝晶运动时溶质浓度在枝晶周围积聚,人为地改变了枝晶的环境,影响了枝晶的生长.为了解决这个问题,本文改进了CA-LB 模型,采用局部动网格技术处理枝晶的运动.

每个枝晶都拥有一个局部的可移动网格,动网格会随着枝晶的长大而增大,动网格的角速度和平动速度与内部枝晶一样(如图1所示).在背景网格中求解流场和温度场,计算枝晶受到的力和转矩.在每个动网格和背景网格中同时求解浓度场,将动网格中的溶质浓度覆盖背景网格,动网格中流体节点的速度是相对速度,

图1 动网格方案示意图Fig.1.Schematic diagram of dynamic grid scheme.

其中,M(θ)是旋转矩阵,v1是流体节点在动网格中的速度,vg是流体节点在背景网格中的速度,vd是动网格在该点的合速度,由方程(19)求出.旋转矩阵被定义为

其中,θ是枝晶的转动角度.动网格的边界从全局网格插值得出.枝晶的生长和捕获在动网格中计算,将枝晶的固相分数在每一步插值到背景网格中.

对于全域而言,可能会出现多个动网格部分重叠的情况.对于重叠的部分:如果都是流体节点,方程(3)中的源项为0,方程退化为简单的对流扩散方程,因为动网格中的速度是相对速度,每个动网格中的流动相对于背景网格都是一样的,所以每个动网格中的浓度值都是相等的;如果有界面节点和固体节点,则选定该值为背景网格的值.本方法的优点是将不连续的枝晶移动转化为连续的熔体流动,通过局部网格中流体的相对流动体现枝晶的运动,避免了直接处理枝晶时溶质的积聚现象.

2.6 算法

模型的耦合算法流程如图2所示,算法如下:

图2 算法流程图Fig.2.Algorithm flowchart.

1) 初始化计算域,给计算域中每个元胞的密度、温度和浓度属性赋初值,设置晶核在计算域中的位置和择优生长方向,计算枝晶的质量、质心和转动惯量;

2) 求解方程(1)计算流场,根据方程(18)—方程(21)计算枝晶受到的力和扭矩;

3) 计算全局温度场和溶质场,方程(2)和方程(3);

4) 采用双线性插值方法将背景网格中的浓度、温度和速度插值到局部动态网格中;

5) 计算局部网格中的温度场和浓度场,在局部网格中求解方程(13)—方程(17)计算枝晶的生长和捕获,更新枝晶的质量、质心和转动惯量;

6) 根据方程(22)和方程(23)计算枝晶的位移和转动角,据此移动局部网格;

7) 将局部网格中的固相分数fs、温度和浓度插值到背景网格;

8)枝晶碰撞处理;

9) 重复(2)—(8)步,直至程序结束.

3 模拟结果与讨论

选取Al-4.7%Cu(质量含量)合金为模拟材料,其物性参数如表1所列.

表1 Al-4.7%Cu(质量含量)合金的热物性参数[26]Table 1.Physical properties of Al-4.7%Cu (weight percent) alloy[26].

3.1 自然对流下单枝晶沉降

将计算区域划分为300× 600 个网格,网格间距为0.5 μm,时间步长为 2.5×10-6s,初始过冷度为7 K,流场边界为非平衡外推边界条件,温度场和浓度场采用无扩散边界条件,固液边界为无滑移边界条件.初始时刻,在计算域(150,450)处设置1 个择优生长角为45°的晶核.图3为自然对流下单枝晶沉降过程中溶质浓度和流体速度演化.在枝晶运动早期,由于枝晶运动速度较小,运动及熔体对流对枝晶生长的影响较小,枝晶生长比较均衡,四个一次枝晶臂形貌比较接近,如图3(a)所示.随着枝晶的下落,枝晶迎流端生长加快,二次枝晶臂相比于顺流端发达,枝晶的下游熔体的溶质浓度高于上游的溶质浓度,如图3(b)和图3(c)所示.

将运动对枝晶生长的影响进行定量分析.图4为图3(b)中O点(150,250)竖直方向上的溶质浓度变化.图中横坐标r为距O点的距离,O点水平线下边为枝晶的迎流端上边为顺流端.从图4可以看出,迎流端的峰值比顺流端的高,这是因为枝晶呈45°生长时,随着枝晶的下落溶质会积聚在枝晶凹形的枝晶臂间;在溶质变化平稳后,顺流端的溶质浓度略高于迎流端,即如图3(b)所示,在枝晶后会形成一段淡淡的拖尾.图5为枝晶尖端生长速度随时间的变化.从图5中可以看出,在枝晶凝固初期,无论是上游尖端还是下游尖端的生长速率基本一致,都处于快速下降,这是由于凝固初期枝晶生长所需的过冷以热过冷为主.随着生长进入稳定阶段,可以看出上游尖端的生长速度明显高于下游尖端,上游尖端的稳态生长速度平均约为2.49 mm/s,下游尖端的稳态生长速度平均约为1.47 mm/s,上游尖端的生长速度比下游尖端的生长速度快了约1.7 倍.从图3可以看出,相比于之前的研究[18,19],采用本文提出的动网格方案计算枝晶移动可以很好地保持枝晶的形貌.对于单个枝晶而言,小的计算域便可以展示枝晶的运动情况,因此设置的计算域较小,应用本模型在Intel Reon E5-2650 v42.20 GHz CPU,内存64 G 的服务器上单线程计算时间为1085 s,而采用openMP加速后,12 线程计算共用时246 s.

图3 自然对流下单枝晶沉降过程中溶质浓度和流体速度的时间演化 (a) t=0.005 s;(b) t=0.015 s;(c) t=0.025 sFig.3.Evolution of solute concentration and fluid velocity during precipitation of single dendrite under natural convection:(a) t=0.005 s;(b) t=0.015 s;(c) t=0.025 s.

图4 t=0.015 s 时,模拟域O 点处竖直方向上的溶质浓度变化Fig.4.When t=0.015 s,the solute concentration change in the vertical direction at point O in the simulation domain.

图5 枝晶尖端生长速度随时间的变化Fig.5.The change of dendrite tip growth rate with time.

3.2 牛顿流体中两圆形粒子的沉降

为了验证碰撞程序的可行性,模拟了牛顿流体中两圆形粒子的沉降.模拟域设置为宽W=2 cm,高H=8 cm,两密度ρp=1.01 g·cm-3,直径D=0.2 cm 粒子位于模拟域水平方向的中间,竖直方向两粒子分别在7.2和6.8 cm 处.域中充满密度ρf=1 g·cm-3,运动黏度νf=0.01 cm2/s 的水.在LB 模型中,将模拟域剖分为200× 800 的网格,松弛时间τ=0.65,上下边界采用非平衡外推格式,左右边界采用无滑移边界条件.这里使用的条件与Feng和Michaelides[27]使用的条件一样.初始时刻,两粒子和流体都处于静止状态,粒子受到重力和浮力的作用开始运动.

图6(a)—(d)分别是t=1.0 s,t=1.5 s,t=2.5 s,t=4.0 s 时两粒子沉降过程的瞬时涡量轮廓.图7(a)和图7(b)分别为粒子质心横、纵坐标随时间变化的图像.多粒子沉降时会出现漂移-接触-翻滚(drafting-kissing-tumblin,DKT)现象,由图6中可以看出,运动早期粒子1 被粒子2 沉降产生的尾涡捕获使其沉降速度加快,到1.5 s 时与粒子1 接触,直至2.5 s 时两粒子发生翻滚.运用本文的模型成功模拟再现了Feng 的结果.在图7中比较了Feng和Michaelides[27]的结果,可以看出粒子下落的轨迹是一致的,存在的一些误差是本文对流体力和碰撞的处理与Feng和Michaelides 的不同,本文计算流体力采用的动量交换法满足伽利略不变性,排除了边界速度的影响,Feng和Michaelides[27]采用的动量交换法不满足伽利略不变性.由于粒子边界速度的影响造成计算流体力时存在误差.

图6 两个球形粒子沉降过程中,四个时间的瞬时涡量轮廓 (a) 1.0 s;(b) 1.5 s;(c) 2.5 s;(d) 4.0 sFig.6.The instantaneous vorticity profiles of two spherical particles in the process of settling at four times:(a) 1.0 s;(b) 1.5 s;(c) 2.5 s;(d) 4.0 s.

图7 粒子质心坐标 (a) 横坐标;(b) 纵坐标Fig.7.Coordinates of the particles centroid:(a) Transverse coordinates;(b) longitudinal coordinates.

3.3 双枝晶弹性碰撞

计算区域为600× 400 个网格,初始过冷度为7 K,温度场边界为绝热边界条件,溶质场为无扩散边界条件.为了便于研究枝晶的碰撞行为,将初始流场设置为静止状态,在计算域(200,200)和(400,230)的位置放置两个晶核,相对于水平方向的择优生长角分别为0°和45°,组成为Cs=kCl.当枝晶生长一定大小时,枝晶不在生长并直接赋给两个枝晶0.1和—0.1 (lattice unit)的水平速度,忽略重力的影响.

两枝晶运动前的位置和形貌如图8(a),左边是枝晶1 右边是枝晶2,两个枝晶受到熔体的阻力作减速运动直至发生碰撞如图8(b),碰撞前枝晶1,2 的水平速度分别为0.08和—0.082,碰撞后速度分别为—0.029和0.024,碰后角速度为0和0.0036 rad.图8(c)和图8(d)显示两枝晶碰撞后的运动情况.碰撞时,由于碰撞点和枝晶1 的质心在同一高度,所以碰撞产生的冲量对枝晶1 的转矩为0,枝晶1不发生转动.

图8 枝晶偏心碰撞形貌图 (a) 0.008 s;(b) 0.0089 s;(c) 0.009 s;(d) 0.00915 s.箭头表示速度v 的大小和方向,灰度表示流体溶质浓度CFig.8.The morphology of dendrite eccentric collision:(a) 0.008 s;(b) 0.0089 s;(c) 0.009 s;(d) 0.00915 s.The arrow indicates the size and direction of velocity v,and the gray scale indicates the concentration of solute C.

图9为0.008和0.009 s 时,枝晶2 周围熔体的平均溶质浓度Cave,x-轴表示距离枝晶2 质心的元胞个数.Cave的计算规则如图10所示,图10中白色格子表示枝晶的质心,x表示距离质心的元胞数.当x=1 时,Cave等于浅灰色格子中所有液相格点的平均浓度;当x=2 时,Cave等于浅黑色格子中所有液相格点的平均浓度;其他与x值对应的Cave同理计算.从图9中可以看出,在0.008 s时枝晶处于静止生长状态,枝晶2 周围的溶质浓度在界面处最大随后缓慢的减少直至稳定.这是因为枝晶生长时在固液界面处向外排出溶质,排出的溶质向周围均匀的扩散,所以形成了这样的浓度梯度.在0.009 s 时,溶质浓度的变化变得不均衡,在枝晶尖端附近浓度最高.因为枝晶2 在碰撞后发生了转动,对周围熔体进行了扰动,使周围的熔体浓度变化混乱,将溶质甩在了枝晶的尖端附近.通过对双枝晶弹性碰撞的计算可以证明本文所用碰撞模型可以表达不规则等轴枝晶的碰撞过程.

图9 枝晶2 附近熔体平均溶质浓度Fig.9.Average solute concentration of melt near dendrite 2.

图10 熔体平均溶质浓度计算规则Fig.10.Calculation rules of melt average solute concentration.

4 结论

本文建立了一个改进的CA-LB 模型,该模型采用动网格技术实现枝晶的运动,每个枝晶拥有一个可移动的局部网格,在局部网格内通过流体的相对流动来体现枝晶的移动;模型中还加入了处理枝晶的碰撞的机制.应用该模型模拟了自然对流下单枝晶的沉降过程,模拟结果表明,枝晶运动的上游端浓度变化梯度大,生长速度加快,二次枝晶臂发达;同时表明该模型可以很好地模拟枝晶的运动过程,并维持枝晶的形貌.对模型的碰撞机制进行了测试,模拟了牛顿流体中两圆形粒子的沉降,成功再现了DKT 现象,证明了碰撞算法的可行性.计算了两枝晶的弹性碰撞,模拟结果表明该模型可以处理不规则枝晶的偏心碰撞.