Acta Metallurgica Sinica 2016 , 52 (10): 1297-1310 https://doi.org/10.11900/0412.1961.2016.00361

Orginal Article

合金凝固过程中显微组织演化的元胞自动机模拟^*

朱鸣芳¹^,, 汤倩玉¹, 张庆宇¹, 潘诗琰², 孙东科³

1 东南大学江苏省先进金属材料高技术研究重点实验室, 南京 211189
2 南京理工大学工程训练中心, 南京 210094
3 上海交通大学上海市先进高温材料及其精密成形重点实验室, 上海 200240

CELLULAR AUTOMATON MODELING OF MICRO-STRUCTURE EVOLUTION DURING ALLOY SOLIDIFICATION

ZHU Mingfang¹^,, TANG Qianyu¹, ZHANG Qingyu¹, PAN Shiyan², SUN Dongke³

1 Jiangsu Key Laboratory for Advanced Metallic Materials, Southeast University, Nanjing 211189, China
2 Engineering Training Center, Nanjing University of Science and Technology, Nanjing 210094, China
3 Shanghai Key Laboratory of Advanced High-Temperature Materials and Precision Forming, Shanghai Jiao Tong University, Shanghai 200240, China

通讯作者: Correspondent: ZHU Mingfang, professor, Tel:(025)83793355, E-mail: zhumf@seu.edu.cn

责任编辑: ZHU MingfangTANG QianyuZHANG QingyuPAN ShiyanSUN Dongke

收稿日期: 2016-08-5

网络出版日期: 2016-10-27

基金资助: * 国家自然科学基金项目51371051和51501091, 中央高校基本科研业务费专项资金项目2242016K40008, 以及东南大学优秀博士论文培育基金YBJJ1627资助

作者简介:

作者简介: 朱鸣芳, 女, 1957年生, 教授

展开

摘要

元胞自动机(cellular automaton, CA)方法能够有效地描述凝固过程中显微组织形貌的复杂演化过程, 且计算效率较高, 展现出良好的实际应用潜力. 近20年来, CA模型取得了很大的进展. 本文简要综述几种模拟凝固组织的CA模型, 包括纯扩散和对流作用下的枝晶生长、共晶凝固、多元合金中的热力学和动力学耦合、枝晶耦合凝固气孔的生长、以及多尺度的耦合模拟. 最后, 对今后CA模型的发展提出作者的几点思考.

关键词： 合金 ; 凝固 ; 显微组织 ; 数值模拟 ; 元胞自动机

Abstract

Microstructure evolution during solidification is a complex process controlled by the interplay of heat, solute, capillary, thermodynamics and kinetics. Computational modeling can provide detailed information about the interactions between transport phenomena and phase transformation. Thus, it has emerged as an important and indispensable tool in studying the underlying physics of microstructural formation in solidification. During the last two decades, extensive efforts have been dedicated to explore the numerical models based on the methods of phase field (PF), cellular automaton (CA), front tracking (FT), and level set (LS), for the simulation of solidification microstructures. The CA approach can reproduce various realistic microstructure features with an acceptable computational efficiency, indicating the considerable potential for practical applications. It has, therefore, drawn great interest in academia and achieved remarkable advances in the simulation of microstructures. This paper gives an overview of CA based models, spanning from the meso-scale to the micro-scale, for the prediction of microstruc ture evolution during alloy solidification. The governing equations and numerical algorithms of CA based models and derived coupling models are summarized, including the calculations of nucleation, growth kinetics, interface curvature, surface tension anisotropy and crystallographic orientation, thermal and solutal transport, melt convection utilizing the lattice Boltzmann method (LBM), the coupling of CA with control volume (CV) method, the coupling of CA with CALPHAD approach for multi-component alloy systems, as well as the approaches for eliminating the artificial anisotropy caused by the CA square cells. The main achievements in this field are addressed by presenting examples encompassing a wide variety of problems involving dendritic growth in pure diffusion and with melt convection, eutectic solidification, microstructure formation in multi-component alloys, dendritic growth with gas pore formation, and multi-scale simulation. Finally, the future prospects and challenges for the CA modeling of solidification microstructures are discussed.

Keywords： alloy ; solidification ; microstructure ; numerical modeling ; cellular automaton

PDF (7078KB) 元数据多维度评价相关文章收藏文章

本文引用格式导出 EndNote Ris Bibtex

朱鸣芳, 汤倩玉, 张庆宇, 潘诗琰, 孙东科. 合金凝固过程中显微组织演化的元胞自动机模拟^*[J]. , 2016, 52(10): 1297-1310 https://doi.org/10.11900/0412.1961.2016.00361

ZHU Mingfang, TANG Qianyu, ZHANG Qingyu, PAN Shiyan, SUN Dongke. CELLULAR AUTOMATON MODELING OF MICRO-STRUCTURE EVOLUTION DURING ALLOY SOLIDIFICATION[J]. Acta Metallurgica Sinica, 2016, 52(10): 1297-1310 https://doi.org/10.11900/0412.1961.2016.00361

显微组织是材料科学和工程关注的核心问题之一, 是材料成分、制备加工过程和性能之间最重要的连接点. 在材料制备过程中, 几乎所有的金属材料都经历了一次或多次的凝固过程. 金属材料的凝固显微组织不仅影响了后续的热加工工艺, 也直接影响了金属制品的最终力学性能.

随着计算机硬件和数值计算技术的迅速发展, 当今数值模拟已成为与实验技术和理论研究并行发展的第三种科学研究方法. 它对材料科学和技术的发展起到了日益重要并且不可或缺的作用. 显微组织的数值模拟可以将相变热力学、动力学、界面能、晶体生长的各向异性、温度场、浓度场、流场、应力-应变场等诸多因素有效地综合起来, 能够可视化地再现合金在相变过程中显微组织和溶质偏析的演化. 通过模拟研究, 可以帮助人们深入理解凝固过程中显微组织形成的内在机制, 为材料成分设计、改进现有工艺以及发展新工艺提供科学依据^[1,2].

目前, 凝固显微组织的模拟方法主要有相场(phase field, PF), 元胞自动机(cellular automaton, CA), 前沿跟踪(front tracking, FT)和水平集(level set, LS)等方法. CA 是一种尖锐界面模型, 能够有效地描述相变时显微组织形貌的复杂演化过程, 具有模型和算法简单、计算效率高等优势. 因此, 在显微组织模拟领域得到了广泛应用.

CA方法最初被应用于模拟形变再结晶过程中的晶粒生长. 1993年, Rappaz和Gandin ^[3]率先将CA方法应用于模拟介观尺度的凝固晶粒组织的演变. 他们所提出的CA模型包含了非均质形核、生长动力学、晶粒的择优取向以及晶粒间的竞争生长等凝固过程的物理现象. 可以模拟出等轴晶和柱状晶凝固组织随时间的演化, 以及柱状晶向等轴晶的转变(CET). 此后, 研究者们应用介观尺度的CA模型模拟研究了不同铸造工艺过程中的晶粒组织形貌^[4-11]. 然而, 介观尺度的 CA模型只考虑固/液界面处的温度过冷, 而未考虑凝固过程中的溶质再分配和界面曲率的影响. 因此, 只能模拟预测晶粒组织形貌, 而不能描述晶粒内的枝晶形貌、微观成分偏析以及第二相(共晶)组织的形成.

1998年, Dilthey和Pavlik^[12] 报道了用CA法模拟焊接熔池中的枝晶生长形貌, 他们的工作堪称CA方法从介观尺度推广到微观尺度的里程碑. 随后, 学者们建立了各种微观尺度的CA模型, 并模拟了凝固过程中枝晶形貌的演化^[13-26]. 经过近20年的发展, 微观尺度的CA方法已从纯扩散的枝晶模型推广到耦合对流的作用; 从模拟枝晶的单相凝固推广到包括枝晶、共晶以及气孔的多相凝固模拟; 从二元合金推广到多元合金体系; 从单尺度模拟推广到多尺度的耦合模拟.

本文简要综述几种模拟凝固显微组织的CA模型, 包括纯扩散和对流作用下的枝晶生长、共晶组织、多元合金中的热力学和动力学耦合、凝固气孔以及多尺度的耦合模拟. 并对CA模型的发展趋势提出作者的一些看法.

1 纯扩散条件下的枝晶生长

枝晶是一种最常见的凝固组织. 枝晶的形成是一个复杂的物理过程, 涉及到传热、传质、动量传输、曲率效应、界面能各向异性、界面热力学、动力学等各种因素的综合作用. 枝晶的模拟也是其它凝固组织模拟的基础. 本节将对用CA法模拟枝晶组织的主要控制方程和算法进行简要的叙述.

1.1 CA模型的简要描述

CA方法的基本思想是将一个物理系统进行时间和空间上的离散, 即将系统看作由一定数量的基本元胞组成. 每个元胞都赋有自身状态的标记以及相关物理量的数值. 依据与其相邻其它元胞的情况, 按照一定的演化规则来决定该元胞的状态, 从而描述系统整体随时间的复杂演变规律. 应用于凝固组织模拟的CA方法包括: (1) 元胞的几何形状; (2) 元胞的状态(固态、液态或界面); (3) 在每个元胞中相关变量(温度、浓度、枝晶生长择优取向等)的数值; (4) 邻位结构; (5) 元胞状态的转变规则. 一般将模拟区域划分为均匀的正方(二维模型)和立方(三维模型)网格. 但针对六方晶系, 例如镁合金中的枝晶模拟时, 也可采用正六边形的网格体系. 二维CA模型的邻位结构包括只考虑最近邻的4邻位(Von Neumann)构型和考虑最近邻和次近邻的八邻位(Moore)构型. 为了模拟出不同择优取向的枝晶形貌, 一般采用8邻位的构型. 同理, 对于三维CA模型, 一般采用将最近邻、次近邻和第三近邻3类邻位都包括在内的26邻位构型. 元胞状态随时间而变化的转变规则遵循形核和生长动力学原理. 一般认为CA法是属于随机性模型. 其实, CA模型的随机性只适用于处理形核方面, 即形核位置和枝晶生长的择优取向是用随机算法处理的, 而生长动力学计算则采用确定性算法.

1.2 形核

形核条件会对凝固组织特征产生重要影响. 一般认为合金凝固始于非均质形核. 描述非均质形核的模型有2种: 连续和瞬时形核模型. 至今大多数CA模型都采用Thevoz等^[27] 提出的连续形核模型, 图1所示的2组Gauss分布曲线用来描述铸型壁和熔体内的非均质形核. 形核数的增量dn随过冷度增量d(ΔT) 的变化呈Guass分布:

$\begin{matrix} \frac{dn}{d (ΔT)} = \frac{n_{i, \max}}{\sqrt{2 π} Δ T_{i, σ}} \exp [- \frac{1}{2} {(\frac{ΔT - Δ T_{i, \max}}{Δ T_{i, σ}})}^{2}] \end{matrix}$ (1)

式中, ΔT_i_,_σ为标准差; ΔT_i_,max为平均形核过冷度; n_i_,max为最大形核数, n_i_,max值可从0到∞积分Guass分布曲线而得到; i=s和i=b分别表示铸型壁和熔体内的非均质形核. 因此, 在任一过冷度ΔT下的形核数n(ΔT)可由下式计算:

$\begin{matrix} n (Δ T) = \int_{0}^{ΔT} \frac{dn}{d (Δ \overset{′}{T})} d (Δ \overset{′}{T}) \end{matrix}$ (2)

式中, $\begin{matrix} Δ \overset{′}{T} \end{matrix}$ 为积分变量. 目前还不能用解析方法确定形核位置和形核Gauss分布中的3个特征参数(ΔT_i_,_σ, ΔT_i_,max和n_i_,max), 因此, 通常还需要根据实验数据库或经验公式来进行估计.

显示原图| 下载原图ZIP| 生成PPT

图1 在铸型壁和熔体内的形核分布

Fig.1 Nucleation distributions for nuclei formed at the mold wall and in the bulk of the melt (ΔT_s,max and ΔT_b,maxare the mean nucleation undercoolings, ΔT_s,_σ and ΔT_b,_σ are the standard deviations; n_s and n_b are the maximum densities of nuclei; the subscripts s and b indicate the surface of mold wall and the bulk of the melt, respectively)

1.3 生长动力学

微观尺度与介观尺度CA模型的主要不同之处在于生长动力学的计算. 微观尺度的CA模型的生长动力学除了温度过冷, 还要考虑成分过冷和界面曲率过冷. 目前主要有以下几种常用的生长动力学计算方法.

根据界面的Gibbs-Thomson关系, Shin和Hong^[16]提出用下式计算界面生长速度v_n:

$\begin{matrix} v_{n} = μ_{k} (n) (T_{m} - T^{*} + C_{l}^{*} m_{l} - ΓKf (φ, θ_{0})) \end{matrix}$ (3)

式中, $\begin{matrix} μ_{k} (n) \end{matrix}$ 为与动力学各向异性相关的动力学系数, $\begin{matrix} n \end{matrix}$ 为界面的法向方向, $\begin{matrix} T_{m} \end{matrix}$ 为纯金属的熔点, T^*为界面温度, $\begin{matrix} C_{l}^{*} \end{matrix}$ 为界面液相成分, m_l为液相线斜率, Γ为Gibbs-Thomson系数, $\begin{matrix} f (φ, θ_{0}) \end{matrix}$ 为与界面能各向异性相关的函数, $\begin{matrix} φ \end{matrix}$ 为固/液界面的法向方向与水平方向的夹角(生长角度), $\begin{matrix} θ_{0} \end{matrix}$ 为枝晶择优生长方向, K为局部界面曲率. 式(3)适用于合金和纯金属, 可以描述枝晶从非稳态到稳态的生长过程, 且不需引入噪音即可模拟出二次分支. 其不足之处是, 对于定量的模拟计算, 需要合理确定动力学系数.

一些学者也采用基于固/液界面的溶质守恒法计算合金中的枝晶生长动力学, 并用下式表示:

$\begin{matrix} (C_{l}^{*} - C_{s}^{*}) ∙ v_{n} = D \vec{\nabla} C ∙ n |_{+}^{-} \end{matrix}$ (4)

式中, $\begin{matrix} C_{s}^{*} \end{matrix}$ 为界面固相成分, C为溶质浓度, D为溶质扩散系数. 该方法的优点是无需引入动力学系数. 但界面溶质守恒条件只有在枝晶达到稳态生长时才能得到满足. 因此, 该模型不能合理描述枝晶非稳态生长时的界面动力学. 此外, 必须引入噪音才能模拟出二次枝晶臂的分支.

根据式(4)的模型, 单博炜等^[28]提出了基于界面热通量守恒的纯物质枝晶生长的动力学方法.

针对式(3)和(4)的不足之处, Zhu和Stefanescu^[29]提出了一种基于界面溶质平衡的枝晶生长动力学的计算方法. 该方法认为, 固/液界面演化的驱动力取决于界面液相平衡成分 $\begin{matrix} C_{l}^{eq} \end{matrix}$ 和实际成分 $\begin{matrix} C_{l}^{*} \end{matrix}$ 之差, 并由下式计算:

$\begin{matrix} Δ f_{s} = \frac{(C_{l}^{eq} - C_{l}^{*})}{(C_{l}^{eq} (1 - k))} \end{matrix}$ (5)

式中, $\begin{matrix} Δ f_{s} \end{matrix}$ 为一个时间步长内的固相分数增量, k为溶质平衡分配系数. 界面液相的实际成分 $\begin{matrix} C_{l}^{*} \end{matrix}$ 由浓度场计算确定, 界面液相平衡成分 $\begin{matrix} C_{l}^{eq} \end{matrix}$ 由下式计算:

$\begin{matrix} C_{l}^{eq} = C_{0} + \frac{T^{*} - T_{l}^{eq}}{m_{l}} + \frac{ΓKf (φ, θ_{0})}{m_{l}} \end{matrix}$ (6)

式中, $\begin{matrix} T_{l}^{eq} \end{matrix}$ 是初始溶质成分为 $\begin{matrix} C_{0} \end{matrix}$ 的平衡液相线温度. 该方法也无需引入动力学系数, 可以合理地描述枝晶从非稳态到稳态的生长过程, 不需引入噪音就可以自然地模拟出枝晶的二次和三次分支. 该模型适用于模拟合金在低Peclet数条件下的枝晶生长.

采用与式(5)类似的方法, Wei等^[30]提出了基于界面温度平衡的纯物质枝晶生长动力学方法:

$\begin{matrix} Δ f_{s} = \frac{(T^{eq} - T^{*}) C_{p}}{L} \end{matrix}$ (7)

式中, $\begin{matrix} T^{eq} \end{matrix}$ 为界面的平衡温度, 实际温度由温度场计算所确定; C_p为比热容; L为凝固潜热.

1.4 界面曲率和各向异性

在微观尺度CA模型中, 需要对固/液界面的曲率进行计算. 一般来说, 计算曲率有以下2种方法. 第一种是根据界面邻位固相分数计算平均曲率K方法^[13]:

$\begin{matrix} K = \frac{1}{Δx} [1 - \frac{2}{N + 1} (f_{s} + \overset{N}{\sum_{j = 1}} f_{s} (j))] \end{matrix}$ (8)

式中, Δx为格子尺寸, $\begin{matrix} f_{s} \end{matrix}$ 为固相分数, N为邻位的个数. 第二种是根据界面固相分数的梯度进行计算^[22]:

$\begin{matrix} K = {[{(\frac{\partial f_{s}}{\partial x})}^{2} + {(\frac{\partial f_{s}}{\partial y})}^{2}]}^{- 3 / 2} \cdot \end{matrix}$

$\begin{matrix} [2 \frac{\partial f_{s}}{\partial x} \frac{\partial f_{s}}{\partial y} \frac{\partial^{2} f_{s}}{\partial x \partial y} - {(\frac{\partial f_{s}}{\partial x})}^{2} \frac{\partial^{2} f_{s}}{\partial y^{2}} - {(\frac{\partial f_{s}}{\partial y})}^{2} \frac{\partial^{2} f_{s}}{\partial x^{2}}] \end{matrix}$ (9)

式(9)的方法具有更明确的物理意义. 但根据作者的经验, 对有些合金系(如Fe-C和Al-Si等), 一般来说, 采用式(9)的方法模拟任意取向的枝晶形貌时,有可能造成枝晶臂的弯曲. 而采用式(8)的方法却可以改善这种现象. 因此, 式(8)仍可以作为曲率计算的一种可行方法.

由于固/液界面能的各向异性, 枝晶生长具有一定的择优取向. 对于立方晶系的合金来说, 界面能呈现四重对称, 并由以下式子表示^[31]:

$\begin{matrix} γ (φ, θ_{0}) = γ_{0} ∙ ψ (φ, θ_{0}) \end{matrix}$ (10)

$\begin{matrix} ψ (φ, θ_{0}) = 1 + εcos [4 (φ - θ_{0})] \end{matrix}$ (11)

式中, $\begin{matrix} γ (φ, θ_{0}) \end{matrix}$ 为固/液界面能, $\begin{matrix} γ_{0} \end{matrix}$ 为平均固/液界面能, $\begin{matrix} ψ (φ, θ_{0}) \end{matrix}$ 为界面能各向异性函数, ε为界面能各向异性强度. 根据Gibbs-Thomson公式, 式(3)和(6)中的各向异性函数 $\begin{matrix} f (φ, θ_{0}) \end{matrix}$ 可由下式计算:

$\begin{matrix} f (φ, θ_{0}) = ψ (φ, θ_{0}) + \frac{\partial^{2}}{\partial φ^{2}} ψ (φ, θ_{0}) = \end{matrix} \begin{matrix} 1 - δcos [4 (φ - θ_{0})] \end{matrix}$ (12)

式中, δ =15ε 为各向异性系数. 生长角度 φ 可根据固/液界面的固相分数梯度用下式计算:

$\begin{matrix} φ = co s^{- 1} [\frac{\frac{\partial f_{s}}{\partial x}}{{({(\frac{\partial f_{s}}{\partial x})}^{2} + {(\frac{\partial f_{s}}{\partial y})}^{2})}^{\frac{1}{2}}}] \end{matrix}$ (13)

由于CA是一种尖锐界面的模型, 用式(12)和(13)模拟的枝晶择优方向有时会受到CA四方网格“各向异性”的影响. Shin和Hong^[16]提出了一种与邻位网格状态相关的形状因子的算法来改善CA网格“各向异性”所造成的影响. 随后, 作者课题组对形状因子G的计算方法做了改进, 并由下式计算:

$\begin{matrix} \{\begin{matrix} G = \min [1, \frac{1}{2} (\overset{4}{\sum_{m = 1}} S_{m}^{I} + \frac{1}{\sqrt{2}} \overset{4}{\sum_{m = 1}} S_{m}^{II})] \\ S^{I}, S^{II} = \{\begin{matrix} 0 (f_{s} < 1) \\ 1 (f_{s} = 1) \end{matrix} \end{matrix} \end{matrix}$ (14)

式中, S^I和S^II分别为最邻近4个网格和次近邻4个网格的状态参数. 几何因子G不仅考虑了邻位固相网格对该界面网格固相分数增长的影响, 同时还考虑了最近邻网格比次近邻网格的更大影响.

也有学者采用一种修正的偏心网格生长算法^[18]或基于MeshTV界面重构的算法^[32]模拟不同的枝晶取向.

将CA模型从二维推广到三维的关键问题是描述三维枝晶生长的不同择优取向. 将耦合界面能各向异性的固/液界面平均曲率定义为权值平均曲率K_wmc, 并表示为:

$\begin{matrix} K_{wmc} = (a (\hat{n}) + \frac{\partial^{2} a (\hat{n})}{\partial θ^{2}}) ∙ K_{1} + (a (\hat{n}) + \frac{\partial^{2} a (\hat{n})}{\partial ϕ^{2}}) ∙ K_{2} \end{matrix}$ (15)

式中, $\begin{matrix} a (\hat{n}) \end{matrix}$ 是界面能各向异性参数, K₁ 和 K₂ 分别对应沿θ和ϕ 2个主方向上的主曲率. 对于使用Euler方法描述三维界面的枝晶生长问题, 仅通过每一时刻三维空间界面上的固相分数信息来求解式(15), 难度较大. 潘诗琰^[33]提出了通过计算Cahn-Hoffman向量 $\begin{matrix} ξ \end{matrix}$ 的曲面散度得到权值平均曲率K_wmc的算法:

$\begin{matrix} K_{wmc} = - \nabla_{s} ∙ ξ (\hat{n}) = (3 ε - 1) ∙ (\partial_{x} n_{x} + \partial_{y} n_{y} + \partial_{z} n_{z}) - \end{matrix}$

$\begin{matrix} 48 ε ∙ [n_{x}^{2} ∙ \partial_{x} n_{x} + n_{y}^{2} ∙ \partial_{y} n_{y} + n_{z}^{2} ∙ \partial_{z} n_{z}] + \end{matrix}$

$\begin{matrix} 12 ε ∙ Q ∙ (\partial_{x} n_{x} + \partial_{y} n_{y} + \partial_{z} n_{z}) + \end{matrix}$

$\begin{matrix} 12 ε ∙ [n_{x} ∙ \partial_{x} Q + n_{y} ∙ \partial_{y} Q + n_{z} ∙ \partial_{z} Q] \end{matrix}$ (16)

式中, $\begin{matrix} n_{x} = \frac{\partial_{x} f_{s}}{|\nabla f_{s}|} \end{matrix}$ , $\begin{matrix} n_{y} = \frac{\partial_{y} f_{s}}{|\nabla f_{s}|} \end{matrix}$ , $\begin{matrix} n_{z} = \frac{\partial_{z} f_{s}}{|\nabla f_{s}|} \end{matrix}$ , 为固/液界面单位法向量的3个分量; $\begin{matrix} Q \equiv n_{x}^{4} + n_{y}^{4} + n_{z}^{4} \end{matrix}$ . 式(16)可以在直角坐标系下, 根据界面的固相分数梯度直接求解, 并且, 可处理二维和三维的权值曲率计算问题.

1.5 温度场和浓度场

在枝晶生长过程中, 根据 $\begin{matrix} C_{s}^{} = k C_{l}^{} \end{matrix}$ 计算固/液界面处的溶质再分配. 溶质扩散方程为:

$\begin{matrix} \frac{\partial C_{i}}{\partial t} = \nabla ∙ (D_{i} \nabla C_{i}) + (C_{l} - C_{s}) \frac{\partial f_{s}}{\partial t} \end{matrix}$ (17)

式中, t为时间, $\begin{matrix} C_{i} \end{matrix}$ 为成分, $\begin{matrix} D_{i} \end{matrix}$ 为溶质扩散系数, 下标 $\begin{matrix} i \end{matrix}$ 表示固、液相, 即 $\begin{matrix} i = s, l \end{matrix}$ . 式(17)右边第二项表示枝晶生长时固/液界面处排出的溶质量.

温度场计算的控制方程为:

$\begin{matrix} ρ C_{p} \frac{\partial T}{\partial t} = \nabla ∙ (λ \nabla T) + ρL \frac{\partial f_{s}}{\partial t} \end{matrix}$ (18)

式中, ρ为密度, λ为热传导系数. 式(18)右边第二项为枝晶生长时固/液界面处所产生的潜热.

一般采用显式有限差分法(finite difference method, FDM) 求解式(17)和(18). 考虑到金属材料的热扩散系数比溶质扩散系数大3~4个数量级, 为了提高计算效率, 可采用两套不同尺度的网格, 对温度场、浓度场和显微组织进行耦合模拟计算. 如图2^[14]所示, 用较大尺度的控制体积(control volume, CV)网格进行温度场计算, 在CV网格内再划分更细的CA元胞进行溶度场和枝晶生长的模拟计算. 通常可认为在同个CV中的所有CA元胞节点的温度是相同的, 也可通过某个CA元胞周围4个相邻CV网格节点的温度进行插值计算. 也有学者采用有限元(finite element, FE)和CA耦合的方法^[4,5]

显示原图| 下载原图ZIP| 生成PPT

图2 CV网格节点和CA元胞的耦合示意图^[14]

Fig.2 Schematic diagram of superimposing cellular automaton (CA) cells with control volume (CV) nodes^[14]

1.6 模拟应用

国内外学者们基于上述CA方法的基本框架, 建立和应用各种CA模型对二维和三维的单枝晶和多枝晶生长进行模拟, 将模拟结果与相关的解析模型和实验结果进行比较验证^[13-46]. 分析研究了等轴和柱状枝晶的竞争生长、一次间距的选择、柱状晶-等轴晶转变、界面能各向异性等物性参数对定向凝固过程中胞晶和枝晶生长形态的影响. 大部分枝晶模拟都是针对立方晶系的合金. 付振南等^[47]、霍亮等^[48]和吴孟武等^[49,50]建立了适合hcp结构特点的CA模型, 对镁合金的单枝晶和多枝晶的生长过程进行模拟研究. 对于hcp结构的CA模型, 需要将计算界面能各向异性函数的式(11)中的四重对称改为六重对称.

图3为用潘诗琰^[33]建立的三维CA模型模拟的Ni-0.4%Cu (质量分数)合金在1.5 K/s冷却速率条件下的等轴多枝晶生长的演化过程. 如图所示, 在凝固的初始阶段, 晶核沿着不同的择优取向生长. 随着凝固的进行, 枝晶主干进一步生长和粗化、并在主干上出现二次分枝. 随后, 枝晶的二次分枝加快生长并逐渐粗化. 在凝固的后期, 枝晶臂互相接触重叠.

显示原图| 下载原图ZIP| 生成PPT

图3 模拟的Ni-0.4%Cu (质量分数)合金在1.5 K/s的冷速下等轴枝晶的演化

Fig.3 Simulated equiaxed dendrite evolution of a Ni-0.4%Cu (mass fraction) alloy with a cooling rate of 1.5 K/s
(a) 1.8 s (b) 3.0 s (c) 4.0 s

2 对流作用下的枝晶生长

由自然对流或强制对流所引起的金属液流动是凝固过程中一种不可避免的现象^[51]. 例如, 在凝固过程中由于凝固潜热的释放和溶质再分配可引起自然对流, 浇注过程或凝固收缩也会引起金属液的流动. 在金属液态成形过程中, 为了促进成分均匀化和细晶组织, 往往通过电磁搅拌、机械搅拌、超声波等方式对金属熔体施加强制对流. 由于流场强烈地改变了凝固过程中浓度场和温度场的分布, 从而对枝晶的生长形貌产生重要影响.

国内外学者们采用CA方法耦合Navier-Stokes (NS) 动量方程和连续性方程的传统流体力学数值计算方法, 对流场作用下枝晶的生长规律进行了模拟研究^[52-61]. 但是, 由于传统的NS求解方法基于连续流体介质假设, 所以, 在应用于对流枝晶的模拟时存在2个问题: (1) 较难处理固/流界面的无滑移边界条件; (2) 对等轴多枝晶进行模拟时, 随固相分数增加, 会使流体计算效率降低, 甚至会造成流体计算发散而使模拟难以继续^[53].

格子Boltzmann方法(LBM)是近20年发展起来的一种新的流体力学计算方法. 与传统的NS求解方法的建模思想不同, LBM将流体抽象为由大量虚拟粒子组成的集合体. 这些粒子按一定的方式在离散格子上进行碰撞和迁移演化, 形成了宏观的流体流动现象. 在LBM中的虚拟粒子用一种分布函数来表示, 在一个满足一定的对称性条件的规则网格中对Boltzmann方程进行离散求解, 以描述虚拟粒子分布函数的演化规律, 然后, 根据分布函数计算流体的密度和速度. 采用基于Bhatnagar-Gross-Krook (BGK)假设的单步松弛时间的二维九速(D2Q9) LB模型^[62] 计算二维的流动问题, 其分布函数的演化方程如下:

$\begin{matrix} f_{i} (x + e_{i} Δt, t + Δt) - f_{i} (x, t) = \end{matrix} \begin{matrix} - [f_{i} (x, t) - f_{i}^{eq} (x, t)] / τ + F_{i} (x, t) \end{matrix}$ (19)

式中, $\begin{matrix} f_{i} (x, t) \end{matrix}$ 为虚拟粒子的分布函数, 表示在位置 $\begin{matrix} x \end{matrix}$ 和时间t出现一个虚拟粒子的概率; Δt 为时间步长; τ为松弛时间; $\begin{matrix} f_{i}^{eq} (x, t) \end{matrix}$ 为相应的平衡分布函数; $\begin{matrix} F_{i} (x, t) \end{matrix}$ 为外力项; $\begin{matrix} e_{i} \end{matrix}$ 为虚拟粒子沿i方向的离散速度. 在二维九速(D2Q9)模型中, 一个虚拟粒子共有9个离散速度(即i=0~8). 平衡分布函数 $\begin{matrix} f_{i}^{eq} (x, t) \end{matrix}$ 和外力项 $\begin{matrix} F_{i} (x, t) \end{matrix}$ 由下式计算:

$\begin{matrix} f_{i}^{eq} (x, t) = w_{i} (\sum_{i}^{=} f_{i} (x, t)) \cdot \end{matrix} \begin{matrix} [1 + 3 \frac{(e_{i} ∙ u)}{c^{2}} + 4.5 \frac{(e_{i} {∙ u)}^{2}}{c^{4}} - 1.5 \frac{u^{2}}{c^{2}}] \end{matrix}$ (20)

$\begin{matrix} F_{i} (x, t) = (1 - \frac{0.5}{τ}) w_{i} [3 \frac{(e_{i} - u)}{c^{2}} + 9 (\frac{e_{i} ∙ u}{c^{4}}) e_{i}] ∙ FΔt \end{matrix}$ (21)

式中, w_i 为权重因子, 其值为w₀=4/9, w_1~4=1/9 和 w_5~8=1/36; c=Δx/Δt 为格子速度, $\begin{matrix} u \end{matrix}$ 为宏观速度, $\begin{matrix} F \end{matrix}$ 为浮力.

LBM也可用于计算由对流和扩散所控制的传热和传质. 相关的LB方程与式(19)和(20)相似, 不同之处在于: 用 $\begin{matrix} g_{i} (x, t) \end{matrix}$ 和 $\begin{matrix} h_{i} (x, t) \end{matrix}$ 替换 $\begin{matrix} f_{i} (x, t) \end{matrix}$ , 用 $\begin{matrix} g_{i}^{eq} (x, t) \end{matrix}$ 和 $\begin{matrix} h_{i}^{eq} (x, t) \end{matrix}$ 替换 $\begin{matrix} f_{i}^{eq} (x, t) \end{matrix}$ , 用τ_D 和τ_α 替换τ, 用 $\begin{matrix} G_{i} (x, t) \end{matrix}$ 和 $\begin{matrix} H_{i} (x, t) \end{matrix}$ 替换 $\begin{matrix} F_{i} (x, t) \end{matrix}$ , 其中 $\begin{matrix} G_{i} (x, t) = \end{matrix} \begin{matrix} w_{i} Δ f_{s} C_{l} (1 - k) \end{matrix}$ 和 $\begin{matrix} H_{i} (x, t) = w_{i} Δ f_{s} L / C_{p} \end{matrix}$ , 为由溶质再分配和凝固潜热所引起的源项. 动量 $\begin{matrix} ρu \end{matrix}$ , 流体密度ρ, 浓度C和温度T, 可通过相关的粒子分布函数计算得到: $\begin{matrix} ρu = \sum_{i}^{=} f_{i} e_{i} + F \frac{Δt}{2} \end{matrix}$ , $\begin{matrix} ρ = \sum_{i}^{=} f_{i} \end{matrix}$ , $\begin{matrix} C = \sum_{i}^{=} g_{i} \end{matrix}$ 和 $\begin{matrix} T = \sum_{i}^{=} h_{i} \end{matrix}$ . 运动学黏度 $\begin{matrix} ν \end{matrix}$ , 溶质和热扩散率D和α, 可分别通过下式与相关的松弛时间τ, τ_D和τ_α相关联: $\begin{matrix} ν = c_{s}^{2} Δt (2 τ - 1) / 6 \end{matrix}$ ; $\begin{matrix} D = c_{s}^{2} Δt (2 τ_{D} - 1) / 6 \end{matrix}$ ; $\begin{matrix} α = c_{s}^{2} Δt (2 τ_{α} - 1) / 6 \end{matrix}$ , 其中 $\begin{matrix} c_{s} = \frac{1}{\sqrt{3}} \end{matrix}$ 为格子声速.

由于LBM中未引入连续流体介质假设, 特别适合于对流枝晶模拟的流场计算. Sun等^[63-66]和Yin^[67]等建立了CA-LBM耦合模型, 用CA方法模拟枝晶生长, 用LBM计算温度场、浓度场和流场, 对强制对流和自然对流作用下的枝晶生长行为进行模拟研究. 与CA-NS模型相比, CA-LBM的计算效率高、且能稳定模拟计算至固相分数达0.9以上. 随后, Sun等^[68-70]将CA-LBM模型与热力学计算引擎PanEngine相耦合, 建立CA-LBM-PanEngine耦合模型, 对三元合金在强制对流作用下的枝晶生长进行了模拟研究.

图4为模拟的Al-2.0%Mg-1.0%Si (质量分数)合金在纯扩散和强制对流下的柱状晶生长形貌. 由图可见, 在凝固初期2种情况下都存在不同择优取向枝晶的竞争生长, 不过, 流动会对后期的柱状晶择优生长产生显著影响. 例如, 在纯扩散条件下, 枝晶A与其左边的2个枝晶, 以及枝晶B, D和E都以相近的速度生长, 而枝晶C受到枝晶B的抑制(图4a). 但在流场作用下, 枝晶A比其左边的2个枝晶生长更快, 枝晶B和D分别被枝晶C和E所抑制. 此外, 在流场作用下, 在枝晶前端有明显的流动漩涡, 影响枝晶两侧的浓度场, 使枝晶主干朝着流入方向倾斜, 并形成了非对称的二次臂生长形貌(图4b).

显示原图| 下载原图ZIP| 生成PPT

图4 模拟的Al-2.0%Mg-1.0%Si (质量分数)合金在纯扩散和强制对流下的柱状晶生长形貌

Fig.4 Simulated columnar dendritic morphologies in pure diffusion (a) and forced convection (b) for an Al-2.0%Mg-1.0%Si (mass fraction) alloy with a cooling rate of 10 K/s and a temperature gradient of 1000 K/m

3 共晶组织

除了枝晶之外, 共晶也是一种常见的凝固组织. 共晶凝固是2个固相从液相中的竞争和协调生长过程, 是合金多相凝固的一种自组织现象. 共晶凝固组织包括规则共晶、非规则共晶和离异共晶. 如包含了共晶反应的Al-Si和Fe-C合金系列, 以其优良的力学性能和良好的铸造工艺性能在工程应用中占有很重要的地位. 这些合金的性能与共晶相的数量和形貌密切相关.

国内外学者们将模拟枝晶组织的CA模型推广到包含一个液相、两个固相的多相系统, 对CBr₄-C₂Cl₆规则共晶、Al-Si非规则共晶的显微组织进行了模拟研究, 探讨了过冷度、生长速率和共晶组织形貌的关系, 其模拟结果不仅和Jackson-Hunt (JH)共晶理论而且与实验结果吻合良好^[71-77].

模拟共晶凝固的CA模型, 需要描述α和β 2个固相从液相中凝固析出时的竞争和协作生长机制. 对于非规则共晶, 还要考虑非小平面和小平面2个固相的不同生长动力学和择优取向关系. 目前, 一般采用类似式(3)的界面Gibbs-Thomson 关系^[72-76]或式(4)的界面溶质守恒法^[77]计算共晶2个固相的生长动力学. 由于α和β相的液相线斜率、Gibbs-Thomson 系数和界面曲率均不相同, 即使在相同的界面温度和界面液相浓度的条件下, α和β相也具有不同的生长驱动力, 从而自然包含了2个相的竞争生长机制. 另一方面, 2个相生长时, 在界面的溶质再分配也不相同: 当α相生长时将排出溶质原子, 增加β相的生长驱动力; 反之, 当β相生长时将吸收溶质原子, 又增加了α相的生长驱动力. 从而实现了2个相的协作生长机制.

球墨铸铁具有良好的力学和铸造性能、并且价格低廉, 在工业上得到了广泛的应用. 球墨铸铁的凝固过程是以离异共晶的凝固方式进行的. 当温度下降到液相线温度时, 初生相(石墨或奥氏体)在液相中形核并生长, 直到共晶温度, 共晶的奥氏体或石墨开始形核和生长. 奥氏体与石墨接触之前, 在液相中以枝晶形式生长. 石墨球与奥氏体发生接触之后, 石墨球迅速被奥氏体壳包围, 并在奥氏体中继续生长. 同时, 奥氏体壳在液相中继续生长, 直到奥氏体壳相互接触, 最终完成凝固过程.

Charbon和Rappaz^[78]最早应用CA方法模拟球墨铸铁的凝固显微组织, 他们基于经典模型用C扩散控制石墨生长. 但其模拟结果显示, 每个共晶晶粒中只包含一个石墨球, 这与实验观测结果不符. 在此模型基础上, Ruxanda等^[79]考虑了初生奥氏体的生长, 对亚共晶球墨铸铁的凝固过程进行模拟, 模拟结果中每个奥氏体晶粒包含多个石墨. 但该模型根据平衡相图确定的界面浓度计算生长动力学, 没有耦合实时的浓度场计算.

上述2个CA模型均未考虑奥氏体枝晶生长的各向异性效应. Burbelko等^[80-82]应用连续形核模型, 考虑了奥氏体和石墨比热容的不同, 将动力学过冷作为固/液界面迁移的驱动力, 用下式来计算一个时间步长中界面网格的新相分数的增量 $\begin{matrix} Δf \end{matrix}$ :

$\begin{matrix} \{\begin{matrix} Δf = vΔt ∙ [Δx {(| cosφ | + | sinφ |)]}^{- 1} \\ v = μΔ T_{μ} \end{matrix} \end{matrix}$ (22)

式中, μ 为动力学系数, φ 为x 轴和固/液界面法向之间的夹角, $\begin{matrix} Δ T_{μ} \end{matrix}$ 为局部动力学过冷度. 他们用该CA模型模拟了球墨铸铁凝固过程中枝晶和石墨的相互作用和不规则球形石墨的形貌, 以及凝固过程冷却曲线的变化. 由于CA网格的各向异性效应, 该模型模拟出的奥氏体枝晶大多沿着45°的方向. 另外, 用式(22)计算生长速度时需要预先确定动力学系数μ.

Zhao等^[83]和张蕾等^[84,85]建立了多相CA模型, 模拟球墨铸铁的离异共晶凝固过程. 在模型中基于类似式(5)的局部溶质平衡法计算石墨球和奥氏体枝晶的生长动力学. 石墨/液相、奥氏体/液相和石墨/奥氏体界面的演化取决于局部的平衡成分和实际成分之差. 在石墨的生长动力学中还考虑了石墨和Fe的密度差. 图5 为模拟的亚共晶球墨铸铁的显微组织和浓度场的演化. 如图所示, 亚共晶球铁的凝固始于奥氏体枝晶的形核与生长. 枝晶的生长释放溶质C, 使区域内C浓度高于初始成分C₀=4.1% (质量分数) (图5a). 当温度降至共晶温度时, 石墨从液相中析出并生长. 当石墨与奥氏体接触时, 石墨迅速被奥氏体所包围, 之后奥氏体从枝晶状转变为团状生长(图5b). 溶质C从液相经过奥氏体壳传输到石墨周围, 使石墨在奥氏体内继续生长(图5c). 凝固末期, 奥氏体壳相互接触, 凝固结束, 每个奥氏体晶粒中包含多个石墨. 之后温度继续下降, 石墨通过吸收奥氏体内的C继续长大, 石墨/奥氏体界面的浓度低于远离石墨区域内浓度(图5d).

显示原图| 下载原图ZIP| 生成PPT

图5 模拟的初始成分C₀=4.1%C (质量分数) 的亚共晶球墨铸铁凝固时的形貌演化

Fig.5 Simulated morphology evolution for a hypoeutectic spheroidal graphite (SG) cast iron with C₀=4.1%C (mass fraction) (Numbers on the figures show the local carbon concentration (mass fraction, %), f_s—total solid fraction, T—temperature)
(a) f_s=7%, T=1163 ℃ (b) f_s=40%, T=1147 ℃ (c) f_s=84%, T=1146 ℃ (d) f_s=100%, T=740 ℃

4 多元合金中的热力学和动力学耦合

在工业中应用的大多数金属材料为多元合金. 研究者们将CA模型从二元系推广到多元系统, 对多元合金的显微组织和微观偏析演化进行模拟研究^[86-93].

如第一节中所述, 溶质分配系数k和液相线斜率m_l是CA模型中的2个重要热力学参数. 在二元合金的CA模型中, 这2个参数可作为常量处理. 但在多元合金系中, 它们会随着温度和成分的不同而发生显著变化. 因此, 建立多元合金系的CA模型必须与热力学相平衡计算相结合, 以获得更准确的热力学参数^[86].

Zhu等^[86]和戴挺等^[87]建立了基于CA和热力学相平衡计算PanEngine引擎相耦合的二维模型, 模拟研究了三元和四元铝基合金的枝晶组织形貌和微观偏析演化. CA模型采用类似式(3)的界面Gibbs-Thomson关系计算生长动力学:

$\begin{matrix} v_{n} = μ_{k} (n) (T_{l}^{eq} - T^{*} - ΓKf (φ, θ_{0})) \end{matrix}$ (23)

根据由浓度场计算所获得的界面液相成分 $\begin{matrix} C_{l}^{*} \end{matrix}$ , 用PanEngine计算出相应的平衡液相线温度 $\begin{matrix} T_{l}^{eq} \end{matrix}$ 和固相成分 $\begin{matrix} C_{s}^{*} \end{matrix}$ , 用式(9)和(12)计算曲率K和各向异性函数 $\begin{matrix} f (φ, θ_{0}) \end{matrix}$ , 据此, 就可用式(23)计算界面的生长速度. 为了提高计算效率, 采用预制表的耦合方式. 即在用 CA模拟显微组织之前, 先用 PanEngine 将液相各组分以l% (原子分数) 的成分间隔分别计算出各成分下的平衡液相线温度及平衡固相成分, 并将计算结果先存储于一个数据表中. 模拟时所需的平衡液相线温度和界面固相组分的热力学数据可由此数据表插值而获得. 石玉峰等^[88]将基于界面溶质守恒的CA模型与PanEngine进行耦合. 用PanEngine获得固/液界面的平衡液相线温度及液相线斜率, 然后用类似式(4)的方法求解枝晶尖端生长速率. 上述模型均未考虑溶质元素之间的互扩散. Zhang等^{[89, 90]}建立了一个耦合热力学、动力学计算的三维CA模型, 考虑了合金元素间相互作用对扩散系数的影响. 他们模拟研究了Al-Cu-Mg合金定向凝固的二次枝晶臂间距, 与实验结果吻合良好. 陈瑞等^[91,92]在统计实验数据的基础上, 通过与热力学、动力学、平衡相图数据库相耦合, 建立了适用于铝合金在较低冷速下凝固时形核密度随最大形核过冷度呈指数性变化的形核模型. 并应用空间坐标变化等算法, 建立了适用于模拟三元铝合金枝晶生长的二维、三维CA模型, 模拟了A356合金在不同凝固条件下的凝固组织, 并与实验结果进行对比, 分析了枝晶三维空间结构的复杂性和多样性.

最近, 王韬涛^[94]将Zhu等^[86]建立的CA-PanEngine耦合模型推广至包含二元和三元共晶转变的多相系统. 图6为应用三元多相CA-PanEngine耦合模型, 对以5 K/s冷却速率凝固的Al-7%Si-1.5%Mg (质量分数)合金中枝晶和共晶组织演化的模拟结果. 图6a~d和e~h分别以各相中的Mg和Si浓度显示. 从图中可以看出, 随着枝晶的生长, 初生相不断向液相中析出Mg和Si溶质, 使2个溶质元素在液相中的浓度升高; 到了共晶凝固阶段, 共晶两相均向液相中析出Mg溶质, 导致Mg在剩余液相中的浓度不断升高. 而对于Si元素来说, 共晶α-Al生长时向液相中析出溶质Si, 而共晶β-Si生长时消耗液相中的Si, 因此, 在共晶凝固阶段, 液相中的Si浓度不会发生显著变化. 在凝固后期发生了三相共晶反应. 即剩余液相转变为共晶α-Al, 共晶β-Si和Mg₂Si相. 图6d和h是凝固结束时的形貌, 其中箭头所指的红色相为Mg₂Si.

显示原图| 下载原图ZIP| 生成PPT

图6 模拟的Al-7%Si-1.5%Mg (质量分数)合金枝晶和共晶组织演化

Fig.6 Simulated morphology evolution of primary dendrites and eutectic structure for an Al-7%Si-1.5%Mg (mass fraction) alloy with a cooling rate of 5 K/s in Mg solute field (a~d) and Si solute field (e~h)
(a, e) f_s=20% (b, f) f_s=55% (c, g) f_s=83% (d, h) f_s=100%

5 凝固气孔

铸件产品的性能不仅与凝固显微组织相关, 还受凝固缺陷的显著影响. 在铝和镁等轻合金铸件中, H气孔是造成其性能降低的主要凝固缺陷之一. 因此, 近几十年来, 凝固过程中显微气孔的模拟预测一直是一个热门的研究课题. 显微气孔模型需要考虑以下因素: (1) 固相(晶粒、枝晶、共晶)的形核与生长; (2) 在固/液界面溶质和H的再分配; (3) 溶质和H的扩散; (4) 显微气孔的形核与生长; (5) 气孔与固相在生长过程中的相互作用.

研究者们也采用CA方法对气孔的形成演化过程进行模拟研究^[95-99]. Lee等^[96]和Dong等^[97]采用CA方法模拟晶粒的形核和生长, 考虑了H在固、液相中的再分配和扩散, 气孔的随机形核和由扩散控制的生长等机制, 模拟了介观尺度的晶粒耦合气孔形成的凝固组织形貌. 李正扬等^[100-102]将介观尺度的气孔CA模型推广到微观尺度, 该模型采用随机形核的方法描述气孔和枝晶的形核, 采用CA模型模拟枝晶和H气孔的生长, 采用FDM计算H和溶质的扩散. 与Lee等^[95,96]建立的模型中由气孔内外的压力差计算气孔生长有所不同, 该模型用局部的H饱和度和H浓度差作为气孔生长驱动力. 应用该模型研究了不同初始H含量、冷却速率等条件下气孔的生长.

亚共晶铝合金的最终凝固组织一般包含初生枝晶和共晶相, 很多铸造铝合金最终的凝固组织中共晶相占到50%以上. 由于共晶反应一般发生在合金凝固过程的最后阶段, 对于显微气孔的形成有重要影响. 最近王韬涛^[94]将先前的CA模型推广至包含气相、液相、枝晶和非规则共晶的多相系统, 模拟Al-Si合金凝固过程中枝晶、共晶和气孔的相互作用. 图7为Al-7%Si (质量分数)合金凝固时枝晶、共晶和气孔耦合生长形貌的演化. 初始H含量为0.7 mol/m³, 冷却速率为10 K/s. 由图7可见, 在凝固初始阶段, 等轴枝晶沿不同择优方向生长(图7a). 当液相中H浓度增加至临界饱和度时, 气孔开始形核和生长. 为了与实验结果^[95]进行比较, 气孔的形核位置预先人为设置. 在与枝晶相遇前气孔呈球状生长(图7b). 与枝晶相遇后, 气孔与枝晶竞争生长, 彼此都受到对方生长的制约, 气孔的形状受周围枝晶的形貌制约变得不规则 (图7c~g). 图7d时达到了共晶温度, 枝晶间的液相达到了共晶成分, 随后发生共晶凝固(图7e~g). 模拟的显微组织和显微气孔形貌(图7g)与实验结果(图7h)^[95]吻合良好.

显示原图| 下载原图ZIP| 生成PPT

图7 模拟的Al-7%Si (质量分数)合金中枝晶生长与H气孔形成的形貌演化

Fig.7 Simulated evolution of dendritic growth with hydrogen pore formation in an Al-7%Si (mass fraction) alloy with a cooling rate of 10 K/s
(a) f_s=10% (b) f_s =15% (c) f_s =30% (d) f_s =55% (e) f_s =65% (f) f_s =80% (g) the final solidified microstructure (h) the experimental micrograph^[95]

然而, 在上述CA模型中, 气孔的位置都是静止不动的. 吴伟等^[103,104]和陈海楠等^[105]将基于多相流的LBM与模拟枝晶生长的CA模型进行耦合, 建立CA-LBM模型, 该模型可以描述在凝固过程中气泡的优先形核位置、气泡的长大、合并、在枝晶间受挤变形以及在液相通道中的迁移等物理现象. 但由于LBM是一个无量纲的模型, 在增加气液两相流的密度比和定量模拟计算方面还有待于进一步的研究. 最近Sun等^[106]将CA-LBM模型推广到三维系统.

6 多尺度耦合

根据模拟尺度, 现有的模拟组织和偏析的模型可大致分为三类: 直接微观模型(10^-6m)、直接宏观模型(10^-3 m)和间接宏观模型(10⁰m)^[107]. 直接微观模型如相场法和微观尺度的CA法等, 可提供枝晶显微组织形貌以及溶质分布. 然而, 由于直接微观模型需要庞大的计算资源, 难以用于模拟大型铸件. 间接宏观模型不能直接模拟相和溶质成分的拓扑分布, 但适合于预测铸锭或大型铸件的温度、宏观偏析、流体流动和缩孔等. 直接宏观模型是将间接宏观模型与模拟晶粒拓扑追踪模型的耦合, 即在较大的网格尺寸下, 采用有限元(FE)或有限差分(FD)等方法计算铸件凝固过程中发生的传热、传质和动量传输等宏观传输过程, 并将铸件内不同部位的降温过程‘保存’下来. 在宏观传输过程计算完成之后, 将计算区域划分为较小的网格, 根据计算‘保存’的冷却曲线, 采用CA方法模拟凝固的形核与长大过程. 这种多尺度耦合方法, 在实际应用研究中具有较大的优势.

1994年, Gandin和Rappaz^[4]最早提出将CA模型与FE方法耦合的思想, 建立了二维CAFE模型, 并将模型推广到三维^[5]. 近年来国内学者也建立了各种二维和三维宏/微观耦合的CAFE和CAFD模型, 对各种凝固工艺过程的温度场和晶粒的形核生长过程进行了模拟并对结果进行了细致的分析^[108-126]. 例如: 在镍基高温合金定向凝固时通过螺旋晶粒选择器, 制得单晶涡轮叶片的过程以及各种因素的影响^[113-118]; 电渣重熔的H13钢锭中的传热、凝固行为和微观组织^[119,120]; 焊接熔池的不同区域中枝晶的竞争生长以及焊缝中枝晶及晶粒组织的形貌演化^[121]; 激光立体成形熔池内固/液界面从平界面失稳到胞/枝晶的非稳态凝固过程^[122]; 高碳钢连铸工艺中溶质扩散控制的凝固过程^[123]; 连续铸造过程中的过热度、铸造速度、二次冷却强度以及电磁搅拌强度等因素对大铸锭凝固组织的影响^[124]; Ti-Al合金胞晶反应的凝固组织^[125]; 微重力条件下的凝固组织和微观偏析^[107,126]等. 为提高计算效率, 上述的一些模型计算中还采用并行计算方法^[110-112]或自适应网格技术^[122].

7 结束语

近20年来, 凝固组织的数值模拟取得了很大的进展. 学者们建立了各种基于CA方法的数值模型, 以模拟研究介观和微观尺度的凝固组织演化. 通过数值模拟, 可获得一些目前的实验手段还难以提供的重要信息, 促使人们更深刻地认识相变过程的内在机制. 实验和模拟的结合是对材料科学问题研究的一个重要和有效途径, 这已成为越来越多材料界学者的共识. 尤其是材料基因组计划 (the materials genome initiative, MGI) 的提出和实施后, 旨在改变传统的材料研发模式, 实现计算技术引领的高效、低成本的材料研发新模式, 为计算材料学的发展带来了新的契机. 作为计算材料学中的一个重要组成部分, 显微组织模拟也面临了一个新的发展机遇. 作者根据在组织模拟领域十多年的研究体会, 对CA法凝固组织模拟在今后几年的发展趋势提出以下几点预测:

(1) 多尺度的耦合: 目前关于凝固组织模拟的多尺度耦合主要侧重于温度场、浓度场和流场的三传宏观模拟和晶粒组织介观CA模拟的耦合. 今后应扩大多尺度耦合的范围, 如将分子动力学(形核、界面能、各向异性以及其它物性参数等方面的模拟计算)、微观尺度CA (显微组织和微观偏析的模拟)、介观尺度CA (晶粒组织的模拟)、FE或FD方法 (传热、传质和动量传输等宏观传输过程的模拟)这些不同尺度的模拟方法进行耦合.

(2) 多场耦合: 目前的CA模型可以耦合温度场、浓度场和流场的三传计算和组织演化的模拟. 近年来发展的许多凝固新技术通过施加磁场、电流脉冲、高压等来改善凝固组织. 为了描述这些凝固过程, 需要在模型中耦合磁场、电流场、外力场等因素对三传和相变的影响.

(3) 发展非平衡凝固的CA模型: 近年来的许多实验研究发现, 通过快速凝固、深过冷、超常凝固等方法, 可使传统材料获得一些特殊的显微组织和优异性能. 但目前的CA模型一般都基于近平衡凝固的假设条件. 因此, 有必要发展能够描述快速凝固和其它极端凝固条件下显微组织演化的非平衡凝固CA模型.

(4) 耦合凝固/熔化效应: 目前模拟凝固组织的CA模型只能描述单纯的凝固过程, 没有耦合熔化效应. 但对于一般的中低速凝固过程, 在糊状区的重熔效应是不可避免的. 如组织粗化和枝晶臂熔断等现象都是在糊状区重熔和再凝共同作用的结果. 此外, 如定向凝固、增材制造、焊接等许多工艺过程都同时包括凝固和熔化效应. 因此, 有必要将目前描述单纯凝固过程的CA模型进行推广, 建立耦合凝固/熔化效应的CA模型, 对糊状区重熔/再凝的显微组织演化和微观偏析进行模拟研究.

(5) 从传统的结构材料拓展到功能材料: 至今的CA组织模拟侧重于传统的结构材料. 近年来各种功能材料发展迅速, 也存在许多与显微组织的预测控制相关的科学问题, 可以借助显微组织的模拟手段进行研究, 为功能材料的设计和研发提供理论指导.

(6) 基础研究和实际应用: 组织模拟本质上属于基础研究的范畴, 但同时也要注重面向实际应用问题. 结合工业实际开展建模与模拟研究.

(7) 组织模拟和性能预测: 组织模拟的最终目标是实现对材料性能的预测和控制, 今后要加强组织模拟和性能预测控制的关联与结合.

(8) 高性能计算: 为了将模拟结果与实际空间和时间尺度的显微组织进行对比, 有必要开展三维大尺度的显微组织模拟. 为了提高三维大尺度模拟的计算效率, 除了目前应用的CPU-MPI技术, 还可将近年来迅速发展的GPU (graphics processing unit)-CUDA (compute unified device architecture)高性能计算技术应用于大尺度显微组织模拟领域.

The authors have declared that no competing interests exist.

〈

〉

合金凝固过程中显微组织演化的元胞自动机模拟*

CELLULAR AUTOMATON MODELING OF MICRO-STRUCTURE EVOLUTION DURING ALLOY SOLIDIFICATION

1 纯扩散条件下的枝晶生长

1.1 CA模型的简要描述

1.2 形核

1.3 生长动力学

1.4 界面曲率和各向异性

1.5 温度场和浓度场

1.6 模拟应用

2 对流作用下的枝晶生长

3 共晶组织

4 多元合金中的热力学和动力学耦合

5 凝固气孔

6 多尺度耦合

7 结束语

合金凝固过程中显微组织演化的元胞自动机模拟^*