2.2 微观组织数值模拟的元胞自动机模型

2.2.1 形核模型

根据形核理论，在均质形核与非均质形核两种基本的形核机理中，对实际生产有实际意义的是整体液相中通过外来质点或基底所发生的非均质形核过程。这种质点可以是金属液中原有的，也可以是人工添加的。在微观组织的模拟过程中，对非均质形核的处理有两种方法，即瞬时形核和连续形核。

1.瞬时形核模型

瞬时形核(Instaneous Nucleation)模型，假设所有晶核都是在同一温度即形核温度下形成的。在一定过冷度ΔT时的形核率可表示为

式中，k1是正比于熔体原子与形核质点碰撞频率的常数；n0是初始形核基底密度； n1是下一时刻晶粒和形核质点密度；T是温度；k2是与晶核、形核质点和液体间界面能相关的常数。

该模型依据经典凝固理论，形核率与过冷度成指数关系，表面存在形核临界过冷度。该模型已用于等轴共晶生长模型，便于计算凝固过程中的固相率，但不能解释凝固过程中其他条件对最终晶粒大小、形态的影响，不能准确预测晶粒度。对于具有很窄结晶区的合金，推荐使用瞬时形核模型，因为这种模型节省运算时间。

2.基于正态分布的形核模型

采用基于正态分布的形核模型，来描述铝带铸轧过程中的晶核数量、位置及其变化。在这个形核模型中，假设形核现象发生在一系列不同的形核位置上，而这些形核位置可由连续的而非离散的分布函数dn/d(ΔT)来描述。在某一过冷度ΔT时所形成的晶粒密度n(ΔT)即可由对该分布曲线的积分求得，即

式中，ΔT是过冷度；ΔTN是平均形核过冷度；ΔTσ是标准曲率过冷度；nmax是异质形核基底的数目。

考虑到型壁和液相内部异质形核基底的差异很大，因此采用两组不同参数的正态分布来进行控制。正态分布的参数一般只能通过参照和估计来确定，在具体的模拟中可根据不同工艺条件对形核的影响进行调整。这一形核模型的具体实现一般是预先按照分布参数把形核位置和形核过冷随机分配到元胞中去，这样在计算开始之后，只需要比较当前元胞的实际过冷和核心元胞的形核过冷就可以判断是否形核。

2.2.2 生长模型

晶体的生长从原子的尺度上说，意味着原子在固液界面上不断地堆积。根据固液界面的微观结构不同，可以将其划分为粗糙界面和光滑界面。在此基础上，晶体可能以连续生长、二维台阶生长、螺型位错生长等机制长大。另一方面晶体的生长则表现为宏观上的固液界面向前推进。这一过程首先受到热扩散的影响，因为结晶释放的潜热需要排走。对于合金而言，晶体的生长还受到溶质扩散的影响，因为固、液相的成分不同造成了界面前沿的溶质再分配。此外，晶体的生长还与界面能作用和界面动力学效应有关。

根据固液界面前沿的温度梯度和成分梯度的不同，界面可能呈现不同形态。当界面前沿的温度梯度为负时，界面处的凸起将处于过冷熔体中，因而其生长得到加强，最终将会破坏界面的稳定形成枝晶组织。当界面前沿的温度梯度为正时，对于纯金属而言，界面处的凸起将处于过热状态而熔化，生长界面将保持为稳定的平界面；但对于合金而言，界面前沿可能会由于溶质富集(对于平衡分配系数小于1的合金)而形成成分过冷，使稳定的平界面生长遭到破坏而产生胞状组织或枝晶组织。在实际的凝固过程中，枝晶组织最为常见。

晶粒的生长不仅涉及尖端的稳态生长行为，而且涉及非稳态分枝结构不断产生和演化。下面不仅对晶粒生长的有关理论问题进行了阐述，并且对相应的理论模型进行了推导。在此基础上建立晶粒生长的数学模型，并提出相应的数值计算处理方法。对于合金而言，在常规铸轧条件下凝固时，晶粒的生长主要伴随着溶质扩散和表面能效应这两个物理过程。其中热扩散和溶质扩散将促使界面失稳，而表面能则作为一种稳定性力量使界面维持平衡，从而使晶粒保持稳态生长。若以固液界面为原点建立移动坐标系，则当界面以速度V沿界面法向n方向移动时，界面前沿的温度场和溶质场可以由以下两式给出，即

但对上述两式的稳态求解并不容易，因为界面的形状并不是已知条件，而需要作为求解的一部分，同时还必须通过一种自恰的方式使界面形状在生长中保持不变。严格的形状恒定求解非常复杂，求解结果表明它非常接近一个旋转抛物面。Invanstov首先对这一旋转抛物面模型进行了严格的数学求解，得到

式中， Ω是无量纲过饱和度； C0 是合金原始成分； T∞是熔体温度； T*是界面温度；Iv（P） 是 Invanstov 函数； PC是溶质 Péclet 数； PT是热 Péclet 数； DL是溶质扩散系数； αL是热扩散系数； R 是晶粒尖端半径； V 是晶粒尖端生长速度；是界面平衡液相成分； L 是结晶潜热； CP是恒压热容； k0 是溶质平衡分配系数。

但通过这一关系式并不能求得晶粒尖端生长速度V和尖端半径R，而只能够求得它们的乘积RV，因此，还需要其他的约束条件来进行约束。Langer等根据界面稳定性理论提出将R=λC作为一个约束条件，其中λC是最小临界扰动波长。这一条件可以写成如下形式，即

式中，mL是液相线斜率；GC是界面处的成分梯度；ξC是PC的函数，在PC<1时，ξC≈1；G是界面处的平均温度梯度；σ*是稳定性常数，其理论值为1/（4π2）；Γ是Gibbs表示Tompson 常数。

这一约束条件也称为临界稳定性判据，根据 IMS （Invanstov solution＆Marginal Stability criterion） 得到的晶粒尖端半径、生长速度和过冷度之间的关系与实验数据取得很好的一致，并且实验数据显示，在不同过冷度下， σ*在很窄范围内变动。这说明了确实表现出稳定性常数的特性，与过冷度和熔体流动无关。这一过冷熔体中自由晶的生长模型最早由Lipton等人提出，因此这个模型被称为LGK模型。下面对这个模型进行了详细的介绍。

对于过冷熔体中生长的晶粒，过冷度由以下几部分组成。

曲率过冷度，这是由于晶粒生长界面的曲率引起平衡液相线温度的改变。

式中，ΔTR是曲率过冷度；K是界面曲率；σ是表面能；ΔSf是体积熔化熵。

成分过冷度，这是由于界面前沿富集了大量的溶质，从而引起界面处的温度与前方熔体的液相线温度形成差距（对于平衡液相成分小于1的合金）。

解式（2-6）得到

由于界面温度T*与熔体温度T∞之间还存在着温差，即

解式（2-7）得到

由上述式（2-12）、式（2-14）和式（2-16）可以得到

根据最小临界波长判据即式（2-11），可得到晶粒尖端半径为

利用=k0及式 （2-6） 可以得到

同样的道理考虑到晶粒尖端的热平衡有

晶粒尖端的温度梯度可以表示为

式中，λS是固体的热导率；λL是液体的热导率。若假设固相和液相的热导率相同，并考虑到Invanstov解中晶粒的生长是一个等温面（GS=0），由式（2-21）和式（2-22）可以得到

将式（2-20）、式（2-23）代入式（2-18） 并认为在PC较小的情况下ξC≈1就得到

在给定总的过冷度下，根据式（2-17）、式（2-24） 可以唯一地求得晶粒的尖端生长速度V与晶粒尖端半径R。这个模型被称为LGT模型。

Kurz等对于枝晶在定向条件下的生长进行了研究，由于定向条件下的温度梯度已知，只需要将式（2-20）代入式（2-18）并把R写成2PCDL/V，整理后可以得到一个关于V一元二次方程，即

式中，PC是溶质Péclet数；Iv(PC)是Invanstov函数。

在给定PC后就可以求得速度V,并由PC=RV/2DL就可以得到晶粒尖端半径R。这一过冷熔体中自由晶的生长模型最早由Kurz等提出，因此这个模型被称为KGT模型。

2.2.3 溶质扩散模型

1.溶质扩散模型的建立

要完整地描述枝晶的生长行为，必须综合考虑到热扩散、溶质扩散、流动、界面能以及高速生长时的动力学效应。但由于它们分别在各种尺度和不同的条件下起作用，此处首先建立一个基本模型，在此基础上可以根据情况对基本模型进行扩展。考虑一个最简单的情况：枝晶在过冷熔体中生长，熔体的过冷度ΔT为已知的恒定值。模型的基本假设如下。

1）考虑到热扩散的尺度比溶质扩散高出3～4个数量级，因此忽略枝晶生长中的热扩散，认为在枝晶生长过程中熔体始终是等温的。

2）在过冷度不大的情况下，忽略动力学过冷，也不考虑高速生长时的动力学效应。这一效应包括固、液相线斜率以及分配系数随生长速度的变化。

3）考虑到溶质在液相中的扩散系数比在固相中高出3个数量级，因此忽略溶质在固相中的扩散。

4）不考虑熔体的流动。

5）认为界面始终处于平衡状态，界面两侧的固液相成分服从

2.界面元胞的生长

在假设1）、2）下，熔体总的过冷度由曲率过冷度和成分过冷度组成。在前面计算曲率过冷度时只考虑枝晶尖端的情况，但此处是对枝晶各处生长界面而言，因此在计算曲率过冷度时要考虑界面能各向异性的影响。对于处于面心立方的铝而言，界面能各向异性可以表示为

式中，σ0是各向同性时的界面能； θ是界面法向与坐标轴正向的夹角；ε是界面能各向异性强度；θ0是择优生长方向与坐标轴正向的夹角。

在考虑了界面能各向异性后，曲率过冷度可以写成

在给定总的过冷度后，界面平衡液相成分可由下式求到，即

由于界面平衡液相成分与远离界面处的平衡液相成分之间存在着成分梯度，因此在这一成分梯度的作用下，会排除多余的溶质并向前推进。通过求解界面速度来求解固相分数增量。有另外的一种方法不用求解界面速度，而是基于界面元胞与液相元胞之间的溶质扩散来获得固相分数增量。考虑到与界面元胞相邻的四个元胞中肯定存在一个液相元胞，从而溶质就会从界面元胞流向液相元胞，流出的溶质为ΔC，则这一溶质平衡关系可以表示为

式中，nL是界面元胞近邻中的液相元胞；CnL是液相元胞nL的成分；Δt是溶质排除的时间间隔，即时间步长；ΔX是元胞尺寸。

固相分数增量与溶质的排除满足如下的关系，即

由于界面元胞在生长过程中，可能还有多余的溶质成分被排入到周围的液相元胞中，因此有必要对液相元胞进行溶质扩散计算。

3.溶质再分配

溶质再分配是整个模型计算中一个非常重要的环节，这一问题以前的解决方法是，在界面元胞内部，可以将其看成由成分为CS的固相fS和成分为CL的液相（1-fS）组成，并且

但需要注意的是：

首先，界面元胞的固相成分CS及液相成分CL并不必须满足

界面元胞的平均成分并不能由下式来表示，即

因为只有对于界面处才满足

而对于所有的尖锐界面模型来说，界面没有厚度，意味着将整个界面元胞的尺寸视为界面厚度，这就违背了尖锐界面模型的基本假设，并且从计算的准确性来说，必须对元胞尺寸进行限制，用这种方法处理溶质的再分配欠合理。

其次， CL并不等于。若假设在每一时间步内，将界面元胞内成分视为，仍然违背了尖锐界面的基本假设。这意味着界面厚度至少是ΔX（1-fS），并且会造成溶质不守恒的矛盾，因为一方面有

而另一方面，从图2-1中可以看出

式中，是界面元胞在上一时刻的液相成分；是界面元胞在上一时刻的固相分数。

第三，多余溶质ΔC不仅可能排到界面元胞剩余液相中，而且有可能排到周围的液相元胞中。若界面元胞在生长过程中需要排走的多余溶质总是排到本胞的剩余液相中，将会使本胞的液相成分不断升高，甚至升高到比当前的界面平衡液相成分还要高。尤其需要注意的是，当本胞剩余液相分数非常少时，会使得1-ΔfS≤0，在数值计算中这会造成剩余液相成分不正常的跳动，严重影响计算结果。Lazaro和Stefanescu 等只是将界面元胞产生的多余溶质排到本胞内部，但并未考虑剩余液相不足的情况。Artemev 等则注意到了这一情况，但他们的处理方法是在界面元胞剩余液相分数大于0.5时，多余溶质排入本胞，而小于0.5时其中一部分多余溶质排入到液相元胞。这一处理方法很难说是合理的，并且在界面元胞的生长速度非常快以至于在某个时间步内的固相转变分数 Δf S仍然超过剩余液相分数1时，显然是无效的。

考虑到以往的研究中出现的这些问题，笔者在保证溶质守恒的前提下，首次提出了一种相对合理的溶质再分配方法，其具体的处理方法如下。

设界面元胞在上一时刻的固相成分为，固相分数为，液相成分为，液相分数为1。当界面元胞在新时刻发生 Δf S的固相转变时，首先判断界面元胞剩余液相是否足以完成 ΔfS的转变。当 ΔfS≥1时，将其修正为 Δf S=1，此时，界面元胞已无剩余液相，令其液相成分 CL=0，而 ΔfS转变需要排出的多余溶质为

这部分溶质被界面元胞以外的液相元胞吸收。同时还需要注意，ΔfS虽然是根据计算出的，但对于界面元胞来说，只存在于一个没有厚度的界面上。因此界面元胞新增的固相只能从界面前成分为的液相中产生，这样才能保证溶质守恒。当界面元胞剩余液相足以完成Δf S的转变，即液相还有剩余时，则首先认为多余溶质完全排到剩余液相内部，使其成分升高到

若此成分低于界面液相平衡成分，则意味着多余溶质可以完全被当前剩余液相吸收。若此成分高于，则将其进行修正为 C L=。这意味着由于多余溶质排入剩余液相后，造成剩余液相成分的升高。但由于界面平衡液相成分最多为，剩余液相己经不可能再吸收更多的溶质，因此需要排到界面元胞以外的溶质为

由式（2-37）或式（2-38）计算出的需要排到液相元胞中的溶质将被这些液相元胞均分。而界面元胞的固相成分和固相分数则按下式进行更新，即

最后，由于在界面元胞的生长过程中，可能有部分多余的溶质被排出到周围的液相元胞中，因此还需要对液相元胞进行溶质扩散计算。由于扩散只发生在液相元胞之间，因此处理相对简单，只需要求解标准的二维非稳态扩散方程，即

在以往的研究中，溶质扩散和再分配问题常常需要处理固相元胞、液相元胞和界面元胞之间复杂的多相交互扩散。而文中提出的这种方法巧妙地避免了这一点，将界面元胞和液相元胞分开处理，在严格保证溶质守恒的条件下，使得整个溶质再分配问题得到了简单而合理的解决。

4.界面曲率的计算

在计算界面平衡液相成分时，需要用到界面曲率。在界面跟踪法中，由于界面是用一系列的锚点来描述的，所以利用锚点的坐标或参数方程在数学上可以很容易地得到界面处的曲率。但对于CA法，关于界面位置和形状的信息只能从元胞的固相分数fS中得到。

Sasikumar等采用了一种基于经验的方法，在界面处画出三倍于网格间距的圆，通过计算圆内平界面和实际界面所包括的固相元胞数目差值，来衡量界面曲率的相对大小。Dilthey等采用了F向量法来计算界面曲率。这种方法也是以界面元胞为圆心，R为半径画圆，F向量的大小将由圆内的固相分数和圆面积比值决定，而其方向则从圆内固相的质心指向当前界面元胞。曲率值由F向量、R值和其他一些经验参数计算得到。Nastac通过对Sasikumar方法的改进提出了用计数法(Counting)计算曲率，并被很多人采用，这个公式是

式中，N是当前元胞的邻胞数，如果取两层元胞则N=24,包括第一层的8个邻胞，第二层的16个邻胞；fS(k)表示第k个邻胞的固相分数；这个公式的曲率值在-1/ΔX～1/ΔX之间。

2.2.4 元胞自动机法的捕获法则

元胞自动机法（CA）是在20世纪50年代初，由计算机创始人、著名数学家Neumann提出的。它的基本思想是一个细胞或系统的基元，依据与其相邻的其他基元的情况，按事先设定的规则来决定自己的状态，从而通过定义局部简单的规则来描述系统整体复杂的演变规律。CA法具有四个要素：①求解区域由具有相同尺寸和几何结构的元胞规则排列而成，在二维情况下，正方或六方点阵是最常见的形式；②元胞具有确定的邻域关系，在二维正方形点阵中，最常采用的是Neumann邻域和Moore邻域，前者由上下左右四个最近邻元胞构成，后者还包括对角上的四个次近邻元胞；③每个元胞用不同的状态值或变量值来标识；④每个元胞自身的状态转变由预先定义的转变规则和邻胞状态决定。

CA法的基本思想简单清晰，在数值计算中处理也很方便，并且很容易和各种物理过程结合起来，因此在晶粒组织和枝晶生长的模拟中都得到了广泛应用。在各种文献中，该方法也称为CA模型或元胞模型，并且还产生了很多改进型。为叙述方便，本书中统称为CA法。在应用CA法模拟凝固组织的过程中，除设定元胞形状、邻域关系、状态值之外，还需要结合形核、生长模型及相应的捕获算法来建立元胞的局部转变规则。

在晶粒组织的模拟中，元胞状态转变可以通过形核或生长两种机制发生。当元胞处于过冷状态且其局部过冷度大于形核过冷度时将发生形核，该元胞的状态被标记为固相，并成为一个核心开始生长，同时该元胞还被赋予一个随机数代表［-π/4，π/4］范围内的择优生长方向。一旦有核心形成，则核心元胞开始沿着择优生长方向生长，并捕获其邻近的液相元胞，使邻胞的状态发生变化，被捕获的元胞将按照核心元胞的生长位向继续捕获其他元胞。当一个元胞的所有邻近元胞都被捕获后，该元胞停止生长。

在处理生长问题时，一般是将枝晶生长理论模型推导的结果进行数学回归处理，将尖端生长速度V表示为界面前沿过冷度ΔT的简单函数，即

这种方法比较简单，但函数的形式和参数的选择会因为合金体系和生长条件的变化而不同，适用性和准确性受到限制。同时需要注意这种方法计算的是枝晶尖端生长速度，在离开尖端的其他各处则无法知道准确的生长速度，因此需要对晶粒的外形轮廓做简化处理，由尖端速度和捕获算法来确定元胞的状态转变规则。