《纯物质枝晶凝固的元胞自动机模型》由会员分享，可在线阅读，更多相关《纯物质枝晶凝固的元胞自动机模型（7页珍藏版）》请在金锄头文库上搜索。
1、纯物质枝晶凝固的元胞自动机模型*单博炜 林 鑫魏 雷 黄卫东( 西北工业大学凝固技术国家重点实验室, 西安 710072)( 2008 年1 月 2日收到; 2008 年4 月 6 日收到修改稿)采用元胞自动机模型模拟了纯物质凝固微观组织演化行为. 与以往采用 KGT 模型简化考虑不同, 采用界面能量守恒方程来自洽地确定液固界面的生长动力学, 考察了模型的网格尺寸相关性和时间相关性, 并对模型的稳定性进行了验证. 采用这个模型模拟了纯丁二腈在过冷熔体中的自由枝晶生长, 并与 Glicksman的等温枝晶生长实验( isothermal dendritic growth experiment,
2、IDGE) 和LM -K 临界稳定性原理进行了对比. 模拟的形貌再现了实验形貌的特征, 并且特征参数与理论结果及实验结果符合得很好.关键词: 元胞自动机, 界面, 凝固, 临界稳定性原理PACC: 6150, 6470* 国家自然科学基金( 批准号: 50471065, 50771083) 和新世纪优秀人才支持计划( 批准号: NCET-06 -0879) 资助的课题. 通讯联系人. E -mail: xlin nwpu. edu. cn11 引言采用适当的模型对材料在凝固过程中的微观组织的形成和演化进行完整描述一直是材料研究中所关注的重要课题. 但是目前的分析模型 1无法描述复杂的枝晶演化过
3、程, 存在较多的局限性. 近年来,随着计算机模拟技术的迅速发展, 深入直观的研究凝固过程中的材料微观组织, 例如枝晶的形成和长大, 已经成为可能. 而在众多的数值模拟方法中, 从20世纪 90 年代发展起来的元胞自动机法( cellularautomaton, CA) 以其简洁快速的优点而受到关注.1993 年, Rappaz 和 Gandin 2首先使用 CA 模型模拟了铝硅合金在均匀温度场下凝固的柱状晶 P 等轴晶转变过程. 不过, 由于其模型较为简单, 因此仅能模拟介观尺度下的晶粒结构. 2001 年, Zhu 和Hong 等 3, 4在 经 典 CA 的 基 础 上 提 出 了 MCA
4、(modified cellular automaton) 模型, 考虑了界面溶质再分配和曲率的影响因素, 模拟了枝晶的完整形貌演化过程. 他们模型的共同点在于都是采用KGT 5模型来描述生长动力学过程. 这中间存在的问题是,KGT 模型是仅适用于枝晶尖端的近似的分析描述,用它来模拟整个枝晶的生长在动力学上是不准确的.最近, Nastac 6, Stefanescu 等 7和 Lee 等8也分别使用 CA 模型进行了合金枝晶的模拟. 包括自由生长和 CET 转变等等. 他们的共同点在于采用更为严格的界面通量守恒条件来确定枝晶的生长动力学行为. 然而, 由于他们模拟的都是金属系合金, 实验中准确
5、测量枝晶的参数尖端半径和生长速度非常困难, 很难将模型与实验结果进行验证. 2004 年以来,李强等 9采用 CA 模型重点针对合金枝晶的自由生长进行了模拟. 郭大勇 10, 11等则进一步考察了多组元和外加流场的影响, 采用 CA 模型分别针对 A- l Si合金和Fe -C -Si 三元合金的凝固过程进行了模拟, 模型不但再现了枝晶的生长、 粗化以及微观偏析, 而且还考虑了外加流场对于枝晶相貌的影响. 但目前的大多数 CA 模型仍然缺乏将枝晶的特征参数, 如尖端半径和生长速度等, 同实验和理论模型进行充分的比较.相比于金属合金体系, 类金属的透明有机物由于可以对其在凝固过程中的组织演化行为
6、进行实时直接观察, 多年以来一直为凝固研究所关注和采用.丁二腈就是这样一种透明有机物, 凝固理论的许多重大进展都是首先以丁二腈凝固的实验结果作为依据而得到学术界认可, 并在其后从多方面由包括金第 58 卷 第 2期 2009 年 2月 1000 -3290P 2009 P58(02)P1132 -07物 理 学 报 ACTA PHYSICA SINICAVol. 58, No. 2, February, 2009 n 2009 Chin. Phys. Soc.属在内的其他合金的实验进一步证实 12) 15. 由于纯丁二腈的各项物性参数准确、 全面, 同时, 可以对其枝晶的尖端半径和尖端速度进行
7、精确的实时测量, 便于对数值模型进行验证.因此, 本文基于Nastac, Stefanescu 和 Lee 等的工作, 发展了一种描述过冷熔体中纯物质枝晶自由生长的 CA 模型, 并将模型与 Glicksman 16, 17等前期采用丁二腈进行过冷熔体自由生长的实验结果以及现有的理论模型进行了比较, 进而对模型进行了验证.21 模型描述2111 宏观物理场对于纯物质过冷熔体枝晶生长, 宏观的物理场只有温度场. 考虑到存在一层液固界面层, 二维情况下的温度场控制方程为9T 9t= A92T9x2+92T9y2+L Q Cp9fs 9t,( 1)其中 T 是温度, Q是密度, t 是时间, x,
8、y 是二维坐标, L 是凝固潜热, Cp是比热, A是热传导系数. 方程右边最后一项只在界面处考虑. 在区域边界处热传导的边界条件为T = T0,( 2)其中 T0是区域边界处的恒定温度.2121 生长动力学生长动力学的核心是确定生长速度. 在许多其他学者的工作中, 对生长动力学的描述都是采用一个描述枝晶尖端的分析模型进行近似处理, 例如Rappaz 等采用的 KGT 模型, Zhu 等采用的 KGT 模型和LGK 模型等, 这使得在处理非尖端的界面行为时容易产生较大的偏差. 为此, 本文应用动力学上严格正确的界面通量守恒条件v =1 Lks9Ts 9z- kl9Tl 9z,( 3)其中9Ts
9、 9z是固相一侧的温度梯度,9Tl 9z是液相一侧的温度梯度, ks和 kl分别是固相和液相的导热系数, v 是界面的法向生长速度.至于界面的温度, 则可以根据过冷度推导出来,如图 1所示. 图中 Tm是纯物质熔点, T*是固液界面温度, T是无穷远处熔体温度, $TR是曲率过冷度, $TT是热过冷度, $T 是总过冷度.图 1 枝晶尖端的温度场对于纯物质枝晶, 当忽略界面原子沉积动力学 过冷度时, 界面的温度可表示为T*= Tm- #K ,( 4)其中 是 Gibbs -Thomson 系数, K 是界面曲率.2131 曲率算法本文确定曲率计算方法采用数元胞算法, 算法的具体形式为 6J =
10、1 a1- 2fs+ENi= 1fs( i )N + 1,( 5)其中, J是平均曲率, N 是邻居元胞的数目, fs是元胞的固相分数, 对于固相元胞, fs是 1, 对于液相元胞fs是 0, 对于界面元胞,fs在 0到 1 之间.2141 捕获规则纯物质熔体的温度当低于熔点时, 晶体一旦形核就会自发生长. 本文不考虑具体的形核过程, 假定在过冷熔体中都能形成晶核. 考虑在熔体中心生 成一个单一晶核, 在模拟时选定熔体中心点的一个元胞设为形核点. 一旦此元胞形核后其状态设为固态, 它的邻居元胞的状态设为界面, 它将会沿着优先 生长方向( 此时是30014方向) 进行生长, 而其他所有的元胞的状
11、态都设为液态. 元胞自动机在每一个离散的时间步长 $t 内进行演化. 一个时间步长 $t 内的固相分数的增量 $ fs为$fs=$t $sVx+ Vy- VxVy$t $s,( 6)其中 Vx, Vy分别是界面法向速度v 在x, y 方向上的分量; $s 是元胞大小. 一旦生长的元胞的fs1 则发生状态转变. 它就由界面变为固态. 当一个元胞 转变为固态之后, 它会捕获周围的液态邻居元胞使它们变为界面元胞, 并使它们和自己具有相同的晶11332 期单博炜等: 纯物质枝晶凝固的元胞自动机模型向. 在捕获时使用 von Neumann 邻居定义( 如图 2所示) , 进行邻居的捕获.图 2 von
12、 Neumann邻居2151 数值模拟采用二维空间的元胞自动机. 在我们的模型中, 整个模拟区域是一个二维平面, 划分成了具有相同大小的一组元胞, 每一个元胞都被一个状态( 例如固态、 液态或界面) 和一组值所标示( 例如温度、 固相分数、 晶向等) . 凝固过程开始以后, 宏观的温度场被首先求解, 从而确定出每一个元胞的宏观参数, 然 后根据这些参数值以及 CA 转变规则, 来确定每一个元胞状态的转变, 而元胞状态的改变又反过来影响宏观的温度场, 这是一个耦合的过程. 整个过程耦合交替进行, 使得模拟出的枝晶不断的长大.热传导方程使用交替方向的隐式中心差分算法来求解. 隐式算法本身保证了网格
13、划分的大小不会影响算法的稳定性. 进行模拟区域是一个正方形区域, 区域大小为 500 L m 500 L m. 整个区域采用正方形网格进行离散, 为了考察网格相关性, 比较了网格大小在 012) 3 L m 范围变化时模拟结果的变化 情况.尖端生长速度的确定方法如下: 记录凝固时间t, 然后测量出生长的格点数 n, 由下式即可计算出来瞬时尖端生长速度 vv =n$s t,( 7)式中 $s 是单个网格的大小. 而尖端半径的确定则采用抛物线拟合的方式. 图 3 所示的是在过冷度$T= 115K 下, 模拟的枝晶尖端点及其抛物线的拟合曲线.图3 在过冷度$T = 1 15K 下, 模拟的纯丁二腈枝
14、晶尖端采样点及抛物线拟合曲线如果拟合曲线的方程为y = ax2+ bx + c,( 8) 则尖端半径按照下式计算R =1 2a.( 9)模拟对象材料高纯丁二腈的热物性参数如下 13: Tm= 58109e , ks= 5136 10- 4cal#cm- 1#s- 1#K- 1, k l=5132 10- 4cal#cm- 1#s- 1#K- 1, A s= 1114 10- 3cm2#s- 1, A l= 1116 10- 3cm2#s- 1, Q s= 11016 g#cm- 3, Q l= 01970 g#cm- 3, C = 2114 cal#cm- 2, $h f= 88511cal#
15、mol- 1, $S= 2168cal#mol- 1#K- 1, c p= 38125 cal#mol- 1#K- 1.31 结果与讨论3111 模型稳定性图 4 显示了在给定过冷度 $T = 115 K 下瞬时凝固速度随时间的变化关系. 在凝固开始时, 尖端生长速度很大, 随着凝固过程的继续, 尖端速度逐渐 变小, 最后收敛到一个固定值. 因此在所有的模拟中, 都是在速度达到稳态之后才确定速度的数值.在模拟的过程中, 网格大小的选择会导致网格 的各向异性强弱不同, 从而会影响到模拟的尖端速度的稳定性. 图 5 显示了网格尺度对模型计算的影响. 将过冷度固定在 $T= 115 K, 整个计算区
16、域的大小固定不变, 将网格大小 $s 从 3 L m 变化到 012 L m, 考察稳态尖端速率随网格的变化. 可以看出, 随着网格的逐渐变小, 速度也逐渐减小, 但是最后收敛到一个稳定的值. 总体上, 当网格大小 $s 1 L m1134物 理 学 报58卷图 4 模拟枝晶尖端生长速率随时间的变化后, 网格对尖端速度的选择已经没有明显影响, 因此在我们所有的模拟计算中都将 $s 取为 1 L m. 虽然初始设置的元胞尺寸并不一定满足全场方程和真实 过程, 上述两种情况都趋于稳定表明模型最后是收敛的.图 5 模拟枝晶稳态尖端生长速率随网格大小的变化3121 凝固形貌的演化图6 显示了 115 K 的过冷度下, 不同凝固时间 下枝晶形貌的演化以及温度场的分布. 可以看出,随着凝固时间的推进, 模拟出的枝晶经历了一个从仅有微小的一次枝晶臂, 到逐渐出现微小的二次臂,再到出现发达的二次臂的发展过程. 在图 6( c) , ( d) 所显示的温度场分布上可以看出, 所有固相的温度几乎相等; 在靠
