枝晶生长的傅里叶谱相场模拟

2021-04-19陈慧琴杨斌鑫

太原科技大学学报 2021年2期

陈慧琴,杨斌鑫

(太原科技大学应用科学学院，太原030024)

几个世纪以来，液相凝固过程中复杂微观结构的形成，如雪花的形成和金属铸态微观结构，一直以来都是科学家们所关注的。尤其枝晶在凝固过程中的微观结构尺度的演化决定了金属的许多物理力学性质[1]，是金属实际加工中具有精确尺寸和优良性能的重要依据和基础科学问题，因为几乎所有的金属体系都起源于液态。作为一种免于追踪固-液界面，准确模拟液相凝固过程中枝晶界面形态的重要工具，相场方法(PF)深受国内外广大学者的追捧和广泛使用。相场法是以Ginzburg-Landau理论为物理基础，由van der Waals[2]所提出的扩散界面模型扩展而来，主要优点在于，相场方法的运用使得计算枝晶生长过程中免于追踪复杂的固-液两相界面，通过使用实值0和1，分别代表液相和固相，用0～1之间表示固液界面，简化了运算过程。Wheeler等人建立了合金的第一个相场模型[3-4]，称为WBM模型。Kim等人[5-6]介绍了另一种模型，采用了薄界面极限，被称为KKS模型。

在求解相场方程数值解的过程中，有限差分法是使用最为普遍的求解离散方法。然而，有限差分法的使用会离散出计算极为复杂的九点差分格式，增加了计算运行的时间。傅里叶谱方法作为一种求解偏微分方程的数值计算方法，将偏微分方程的解近似地展开成光滑函数的有限级数展开式，再由所此有限级数展开式与偏微分方程，求解得到有限级数展开式的系数方程组。使用傅里叶谱方法求解相场方程，不需要继续计算复杂的九点差分格式的算子，并且对于与周期性边界条件，傅里叶级数的计算更为便捷，谱方法的精度更为准确，相比有限差分方法，傅里叶谱方法的使用大大加快了求解相场方程的运算速度，提高了数值计算与数值模拟的效率。

1 相场方程的建立

1.1 相场控制方程

下面给出的模型基于Kobayashi[7]最著名的工作，其作为树枝状凝固的最早相位场模型之一。该模型包括两个变量：一个是非守恒相场参数φ(r,t)，它在固体中取值1，在液体中取值0.另一个是温度场T(r,t)，也随着凝固的进展而进化。其中，r是空间位置，t是时间。

在相场方程演化过程中，假设了与时间有关的Ginzburg-Landau或Allen-Cahn方程：

(1)

其中τ是松弛时间，Χ(φ,m)是系统自由能函数。

之后，考虑Ginzburg-Landau型自由能[8]：

(2)

上述积分项中的fgrad(φ)是梯度能量，flocal(φ,m)是局部自由能。

其中局部自由能的表达式为：

flocal(φ,m)=

(3)

flocal(φ,m)对φ的一阶导数为：

(4)

梯度能量的表达式为：

(5)

其中ε是各向异性梯度能量系数，它决定界面层厚度。为了引入界面各向异性，假定在界面处ε依赖于外法向量的方向[9]。

而ε的值可通过以下表达式计算：

(6)

σ(θ)=1+δcos(j(θ-θ0))

(7)

上式中δ是各向异性强度；j是各向异性模数；θ0是初始偏移角，取作常数。这里角度θ被定义为：

(8)

(4)式中的m，假定它取决于过冷度和温度[7]，它的表达式为：

(9)

其中α是一个正常数，Teq是平衡温度。

将(2)，(3)，(4)，(5)代入到(1)中，整理可得：

(10)

1.2 温度场控制方程

由焓守恒定律导出的温度场T(r,t)的演化方程：

(11)

上式中T无量纲化的温度，所以特征冷却温度是零，平衡温度是1[10].κ是一种无量纲潜热，它与潜热成正比，与冷却强度成反比。为了简单起见，无论在固相中还是液相中都将扩散常数设置为相同的。

2 用傅里叶谱方法演化相场模型

2.1 傅里叶谱方法

傅里叶空间中空间导数的一般关系如下[11]：

其中{·}k是括号内量的傅里叶变换，k是k个傅立叶模式的系数。

特别的，当n=1时且在二维空间中有[12]：

其中F(A)=fft(A)，F-1(A)=ifft(A)，fft傅里叶正变换，ifft是傅里叶逆变换。

根据上述表达，相场模型中的空间一阶偏导、时间一阶偏导、拉普拉斯算子的傅里叶变换如下：

空间一阶偏导：

时间一阶偏导：

对于时间的偏导数采用向前差分。

拉普拉斯算子：

F(2φ)=-k2{φ}k，F(2T)=-k2{T}k

若使用有限差分九点格式离散拉普拉斯算子则有：

∇2ψ(i,j,n)=ψxx(i,j,n)+ψyy(i,j,n)=

很直观的就可以观察到，傅里叶谱方法在运算上比有限差分简单的多。

2.2 相场模型的离散

用傅里叶谱方法变换方程(10)两边[13]：

(12)

(12)的半隐式形式是：

(13)

Δt是时间步长n+1和n之间的时间增量[14]，通过重新排列(13)，得：

(14)

对温度场方程(11)进行傅里叶谱变换：

(15)

(15)的半隐式形式是：

(16)

Δt是时间步长n+1和n之间的时间增量，通过重新排列(16)，得：

(17)

3 结果与讨论

为了证明傅里叶谱方法可用于求解此类方程，对求解结果进行了数值模拟。

下面给出处理后方程所需要的相场模型参数的取值：

这些模拟是针对一个方形模拟单元进行的，该模拟单元Nx=Ny=300，网格步长为dx=dy=0.03，时间步长为Δt=0.000 1，时间迭代步数为4 000.

初始结晶核的大小是x2+y2<10.0，如图1所示。

图1 初始晶核

实验所需的其他无量纲化参数在表1中已给出，利用Matlab进行编程，并用Paraview可视化得到下列图像(如图2所示)，分别是在时间迭代步数为400、1 400、2 400、3 400时的图像，其中图2第一行是枝晶凝固的时间演变，图2第二行是温度场的时间演变。在无噪声的情况下，采用粗糙网格进行数值模拟，得到了定性的结果，与杨斌鑫等人的结果非常一致。

表1 模拟所需的无量纲化参数

张晨辉提到随各向异性强度δ增加，枝晶尖端不再分叉，且在主枝上出现了侧枝，主枝变细，优先生长；而平衡温度Teq增大,枝晶侧枝的生长速度变快，侧枝的尖端半径在减小。通过改变模拟时的各向异性强度、平衡温度等参数，得到了如图3、图4所示的实验结果，其中a1-a3(时间迭代步数分别为900、1 900、2 900)以及b1-b3(时间迭代步数分别为600、1 600、2 600)是数值结果，而a4和b4是实验结果，与实验结果对比，可以看到数值结果与实验结果非常一致。

图2 各向异性模数j=6，各向异性强度ε=0.02时，枝晶凝固的时间演变和温度场的演变

图3 各向异性模数j=6，各向异性强度ε=0.06时，所获得的片晶形态

图4 平衡温度Teq=0.9时枝晶凝固的时间演变和实验结果

在图5中(分别是在时间迭代步数为400、1 400、2 400、3 400时的图像)，在界面上加了小噪声，这样更接近现实。采用了随机数字序列arφ(1-φ)表示噪声，其中a为噪声幅值，此处取a=1，x均匀分布在[-0.5,0.5]这个区间上。对比图2第一行发现，加了噪声后的图像主枝稍微粗壮了些，侧枝数量减少，枝晶生长速度减缓。