CA-LBM模型模拟自然对流作用下的枝晶生长
- 格式：pdf
- 大小：828.56 KB
- 文档页数：12
CA-LBM模型模拟自然对流作用下的枝晶生长1杨朝蓉，孙东科，潘诗琰，戴挺，朱鸣芳东南大学，江苏省先进金属材料高技术研究重点实验室，南京（211189）E-mail: zhumf@摘要：本文建立了一个二维的元胞自动机－格子玻尔兹曼方法(cellular automaton-lattice Boltzmann method，CA-LBM）的耦合模型，对自然对流作用下枝晶的生长行为进行模拟研究。
本模型采用CA方法模拟枝晶的生长，采用LBM对自然对流及由对流和扩散控制的溶质和热传输进行数值计算。
通过计算方腔自然对流问题对模拟自然对流、溶质和热传输的LB模型进行了验证。
应用所建立的CA-LBM耦合模型模拟研究了合金中单枝晶和多枝晶在自然对流作用下的生长规律，并将单枝晶的上游尖端的稳态生长模拟数据与解析模型的预测结果进行了比较。
结果表明，模拟结果和理论预测值吻合良好，自然对流的存在会对枝晶的生长产生重要的影响。
关键词：数值模拟；枝晶生长；自然对流；元胞自动机；格子玻尔兹曼方法中图分类号：TG111.41.引言枝晶是一种典型的凝固组织。
众所周知，凝固过程中液态金属除强制对流外，还存在由温度梯度和浓度梯度引起的自然对流，这种流动在凝固过程中是始终存在的，它将直接影响到材料凝固后的组织，成分偏析以及凝固疏松和孔洞等缺陷的分布。
近年来，国内外学者们针对纯金属和合金，应用相场(phase field, PF)模型、前沿跟踪(front tracking, FT) 模型或元胞自动机 (cellular automaton, CA)模型耦合基于Navier-Stokes(NS) 方程的流场数值计算方法研究了对流作用下枝晶的生长规律[1-8]。
但由于传统流体力学计算模型是建立在宏观连续介质的假设条件下的，比较适合于计算单相流体流动，而在凝固过程中存在固、液两相，当固相分数比较大的时候，流场计算不易收敛，影响了计算效率。
此外，利用NS方程方法计算自然对流格式复杂，较难实现。
格子玻尔兹曼方法[9](Lattice Boltzmann Method, LBM)是上世纪八十年代出现并迅速发展起来的一种新的流体计算方法。
与传统的流体力学计算方法不同，LBM不是从宏观方程出发，而是针对流体系统建立一种简化的微观或介观动力学模型。
与传统的计算流体力学方法相比，LBM具有计算效率高、本质并行、容易实现等显著优点，这些优点吸引了众多领域的学者对其理论和应用进行研究。
目前，LBM在热流动、多相流、多孔介质流等众多领域中都得到了比较成功的应用[10]。
近年来，LBM在凝固领域也逐渐得到了应用。
国外一些学者将LBM与PF相耦合模拟了纯物质在对流作用下的枝晶生长规律[11-13]。
国内有学者开展了将LBM应用于凝固的宏观流动模拟研究[14-15]。
但至今国内关于将LBM应用于凝固微观组织的模拟研究方面的报道还较少。
此外，目前国外关于LBM在对流枝晶模拟方面的研究大多是针对纯物质，而应用LBM模拟合金在自然对流作用下枝晶生长规律方面的报道还较少。
CA模型能模拟实验所观测到的各种凝固微观组织特征并且具有较高的计算效率。
因而近年来在微观组织模拟方面备受关注并得到了较大的发展[16-18]。
本文将计算枝晶生长的CA方法与计算流场、浓度场和温度场的LBM进行耦合，建立了一个二维的CA-LBM模型，1本课题得到国家自然科学基金50671025、江苏省自然科学基金BK2006105和教育部博士点基金的资助。
对合金凝固过程中在自然对流作用下单枝晶和多枝晶的生长规律进行模拟研究。
2. 模型和算法2.1 模拟自然对流的格子玻尔兹曼模型本文采用基于单步松弛时间的2维9速LB 模型，即LBGK(BGK: Bhatnagar ，Gross ，Krook)-D2Q9模型[19-20]计算二维的流动问题，其分布函数演化方程为：i eq i i i i i F t f t f t f t t t f +−−=−∆+∆+−)],(),([),(),(1x x x c x τ i= 0, 1… 8, (1) 式中，()t f i ,x 为粒子分布函数，()t f eqi ,x 为相应的平衡分布函数，τ为无量纲的单步松弛时间，i F 为外力分量,i c 为粒子在i 方向的运动速度，定义为{}0,00=c ，()[]()[]{}2/1sin ,2/1cos 41ππ−−=−i i c c , ()[]()[]{}4/92sin ,4/92cos 285ππ−−=−i i c c其中，格子速度t x c ∆∆=,x ∆为空间步长，t ∆为时间步长。
式(1)中平衡分布函数()t f eq i ,x 由密度ρ和速度u 决定]2)(:1[422ss i i s i i eqic c c f I c c uu u c −+⋅+=ρω (2) 其中权重系数ω0=4/9，ω1-4=1/9，ω5-8=1/36，格子声速3/c c s =。
宏观量密度和速度可由分布函数求出：),(80t f i i x ∑==ρ，F u 280tf c i i i ∆+=∑=ρ (3) 无量纲单步松弛时间τ 与流体运动黏度相关50)(2.t c s +∆=ντ，计算中取ν=6.0×10-9 m 2/s 。
采用Boussinesq 近似，即假设除浮力项外，流体的密度ρ是常数，而浮力项中的ρ是局部温度和浓度的线性函数。
在此假设下，粘性耗散可以忽略，并且重力项中的密度可以假设为温度梯度和浓度梯度的线性函数[9]，即：)]()(1[000C C T T c T −−−−=ββρρ (4)其中0ρ、0T 、0C 分别为初始的密度、温度和浓度，T β、C β为温度和浓度膨胀系数，C 为溶质浓度、T 为溶质温度。
则产生自然对流的浮力项为)()(0000C C T T C T −+−=βρβρg g F 其中g 为重力加速度。
LB 演化方程式(1)中的外力分量为：F c u c u c ⋅∆⋅+−−=t c c F i si s i i i ][)211(22ωτ (5) 由温度梯度和浓度梯度产生的自然对流的流动特征可以用4个无量纲数描述：温度Rayleigh 数ναβ/3TL g Ra T T ∆=、浓度Rayleigh 数D CL g Ra C C νβ/3∆=、Prandtl 数ν=Pr 和Schmidt 数D Sch ν=。
其中，α为热扩散系数、D 为溶质扩散系数，L 为计算区域的长度，T ∆和C ∆分别为液相中的最大温度差和最大浓度差。
2.2 模拟溶质和热传输的格子玻尔兹曼方程本文采用被动标量法的LBM 计算浓度场和温度场，其演化方程为：i eq i i D i i i G t g t g t g t t t g +−−=−∆+∆+−)],(),([),(),(1x x x c x τ i = 0,1…8, (6) i eq i i i i i H t h t h h t t t h +−−=−∆+∆+−)],(),([)t ,(),(1x x x c x ατ i =0,1…8, (7) 式中，),(t g i x 和),(t h i x 分别表示浓度和温度的分布函数。
D τ和ατ表示浓度场和温度场的无量纲单步松弛时间，5.0/)5.0(+−=Sch D ττ，5.0/)5.0(+−=Pr ττα。
i G 为溶质源项，i H 为温度源项，分别由下式求得：s i i k C G ϕω∆−=)1(，P s i i C H H /∆∆=ϕω (8)其中，C p 为热容，k 为溶质分配系数，s ϕ∆为一个时间步长中界面网格的固相分数增量。
计算中分别取：k=0.103，3105×=∆ΗJ/m 3， C p =1937.5J/(Kg.K)。
与),(t g i x 和),(t h i x 相对应的平衡分布函()t g eq i ,x 和 ()t h eqi ,x 为：()]2)(:1[,g 422s s i i s i i eqic c c C t I c c uu u c x −+⋅+=ω (9)()]2)(:1[,422s s i i si i eqic c c T t h I c c uu u c x −+⋅+=ω (10) 其中，ωi 权重系数与式(2)中相同。
浓度C 和温度T 可由其相应的分布函数求出：∑==80),(i i t g C x ，∑==80),(i it h T x (11)这样，就得到了一个计算由温度梯度和浓度梯度产生自然对流的LBM 模型。
2.3 模拟枝晶生长的元胞自动机方法在CA 方法中，枝晶的生长由局部过冷驱动。
在t n 时刻，局部过冷度()n t T ∆由下式计算()()()()()()n n n l n t K C t C m t T T t T θΓ−−⋅+−=∆0 (12)式中，T l 、m 、C 0 和 Γ(θ)分别是液相线温度，液相线斜率，初始浓度，Gibbs-Thomson 系数。
其中，m 取-2.16k/(mol.%)，溶质分布C (t n )与温度分布T (t n )通过前面所介绍的计算动量、质量、能量传输的LBM 模型求出。
界面曲率K (t n )由下式求出：⎥⎦⎤⎢⎢⎣⎡∂∂⎟⎟⎠⎞⎜⎜⎝⎛∂∂−∂∂⋅⎟⎠⎞⎜⎝⎛∂∂−∂∂∂⋅∂∂⋅∂∂⋅⎥⎥⎦⎤⎢⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛∂∂+⎟⎠⎞⎜⎝⎛∂∂=−22222222/3222)(x y y x y x y x y x t K s s s s s s s s s n ϕϕϕϕϕϕϕϕϕ (13) 本模型采用传统的尖锐界面模型计算界面生长速度和过冷度之间的关系)(n k g t T V ∆⋅=µ (14)式中，µk 为界面动力学系数，通过一种反推法进行确定。
ZS (Zhu and Stefanescu) 提出了一个基于界面平衡浓度和实际浓度之差确定界面生长动力学的模型。
因此，该模型不需界面动力学系数，并已经过验证，具有良好的定量模拟能力[21]。
将CA-LBM 和ZS 模型在完全相同的纯扩散条件下进行模拟计算，比较两个模型计算的枝晶尖端速度，以此确定在本文的物性参数计算条件下取 µk ＝0.008m/(sk)。
为考虑界面能的各向异性因素，Gibbs-Thomson 系数()θΓ用下式计算()[]{}04cos 151)(θθεθ−−Γ=Γ (15)⎟⎟⎠⎞⎜⎜⎝⎛∂∂∂∂=x y arctan s s ϕϕθ (16) 式中Γ、ε分别是平均Gibbs-Thomson 系数和表面能各向异性强度系数，取810626−×=Γ.mk ，0266.0=ε。
