分享：基于序参量梯度的改进相场模型及对大尺度体系马氏体相变的模拟研究
魏铖
预测和调控材料的微观组织进而设计和制备不同工程应用所需的材料已成为材料科学与工程研究的一个热点。马氏体相变作为多种材料中一种重要的结构转变,对其微观组织的实验研究已获得许多进展[1,2,3,4,5]。但由于马氏体相变速率快,实验难以准确捕捉随时间而变化的组织学和动力学信息;加之马氏体相变受多种复杂因素的影响,包括温度场、应力场和磁场等,这对实验研究(如原位观察)提出了严苛的要求,且实验难以准确反映单因素的影响。目前计算机模拟已成为一种有效的替代手段。构建合理高效的数值模型,对马氏体相变组织学和动力学进行模拟研究,是对实验研究的有效补充,具有重要的理论意义和应用价值。
相场法是目前模拟微观组织演化最有效的数值分析方法之一,它能够耦合温度场、应力场、磁场以及其它外部因素对微观组织的影响,且无需直接对具有复杂形貌的界面进行追踪[6]。Wang和Khachaturyan[7]基于微观弹性理论[8],开创性地建立了适用于马氏体相变的相场模型,成功预测了马氏体形核和长大过程。随后,运用相场模型对马氏体相变的模拟研究逐步发展起来,包括建立了超弹性[9,10]、弹塑性[11,12]的相场模型,分析了惯性力[13]、自协调效应[14]、孔隙[15]、非均匀形核[16]等因素对马氏体相变的影响。然而,目前对马氏体相变的相场模拟研究仍局限于较小尺度(模拟体系格点尺度为纳米级),主要原因是相场模型中与实际相变体系中界面能相关联的界面宽度往往较小;同时,模拟过程需要在界面区域划分足够的格点数以保证相场界面的扩散性,使得相场模型对大尺度马氏体相变的模拟需要极大的计算量。已有的增大相场模型模拟尺度的方法主要基于实际体系中的界面能密度[17],通过调整相场模型中的体积化学能系数与梯度能系数,进而增大界面宽度并扩大相场模拟尺度。由于在马氏体相变相场模型中体积化学能系数与梯度能系数由体系属性所决定,难以直接通过调整系数的方式增大模型界面宽度。相场模型的重要特征是界面具有扩散特性,通过序参量在数值上的梯度进行实现;因此,从序参量梯度入手是一种值得尝试的方法。
本工作基于序参量梯度值在界面区域远大于非界面区域的特性,构造出与序参量梯度相关的全局修正函数,通过调整界面区域的体积化学能系数与梯度能系数,实现对大尺度体系下的马氏体相变的高效率模拟。具体研究内容主要包括：(1) 分析采用全局修正函数后马氏体新相的生长速率;(2) 研究全局修正函数对变体间位向关系形成的影响;(3) 探讨体系大小与全局修正函数对NiTi合金B2-R相变微观组织形貌的影响。
马氏体相变相场模型中决定相变演化路径的能量包括体积化学自由能、梯度能和弹性应变能。其中,体积化学自由能促进马氏体相变;梯度能与体积化学自由能在界面区域的变化量共同组成界面能,控制界面移动;弹性应变能影响马氏体的形核取向及生长方向。
体积化学自由能密度φchem通常由Landau多项式[18]表示：
式中,Δf为母相与马氏体相的体积自由能密度差;ηp (p=1, 2,
体系的梯度能密度φgrad通常用梯度项表示：
式中,kη为序参量场的梯度能系数,由界面能密度确定。
界面能由体积化学能的变化量和体系梯度能组成,界面能密度γ可表示为[20]：
式中,l1和l2分别为界面两侧的位置;l为界面的法向坐标;Δφchem为界面区域体积化学能密度的变化量,其值为体积化学能密度φchem与马氏体相具有的体积化学能密度之差;因此,Landau多项式θ (η1, η2,
当界面宽度趋于稳定时,梯度能密度与体积化学能密度的变化量相等,界面能密度可采用下式[17]计算：
式中,ap和bp分别为界面两侧序参量ηp的值,通常取0或1。参数ap、bp、Δf和θ (η1, η2,
在马氏体相变中,弹性应变能Gelast主要包括新相和母相间由于晶格畸变产生的弹性应变能及马氏体变体之间自协调而产生的能量,可表示为：
式中,V为所研究体系的体积;Cijkl表示弹性模量张量(i, j, k, l为空间向量);
式中,应变张量
将马氏体相变中的体积化学能、梯度能和弹性应变能构成的Ginzburg-Landau函数代入非守恒序参量场的Allen-Cahn方程[21]?,得到以下控制方程：
式中,t为体系组织演化的时间;L为相场动力学系数;Δ为Laplace算子;
式中,
相场模型中的界面宽度λ、母相与马氏体相的体积自由能密度差Δf及梯度能系数kη关系可表示为[17]：
其中,b和a为界面两侧序参量值,通常分别取1和0;Δf和kη均由体系属性所决定,难以直接调整[22]。界面宽度λ内通常需要3~5个格点以保证模拟的准确性。若采用的格点间距d大于λ,则界面宽度将因格点限制被增大至λ′,致使原有的界面平衡被破坏。根据式(4),若λ′=Nλ时,积分长度相应地也增大N倍,体积化学能密度变化量的积分
由此可见,当格点间距d大于λ时会高估界面能密度,影响相变演化路径。由于在大尺度体系下加密格点需要极大的计算量,只有选取合适的方法增大界面宽度,才能实现高效模拟。
相场模型的模拟区域可分为界面区域与非界面区域(马氏体相或母相区域)。非界面区域序参量为0或者1,其梯度为0;而在界面区域中序参量沿界面法向有显著变化,其梯度较大。因此,根据序参量的梯度值,能够很好区分界面区域与非界面区域。基于此特性构造序参量梯度的全局修正函数gp;在不改变界面能密度的前提下,仅针对界面区域调整参数Δf和kη,在增大界面宽度进而扩大模拟体系尺度的同时,尽量避免了参数调整影响到基于界面能密度的体系属性。全局修正函数gp的表达式为：
其中,λ为理论界面宽度,λm为调整后的界面宽度;νp为等效梯度：
式中,
图1为λm=10、λ=1和λm=8、λ=1时,gp分别随νp变化的曲线。由图1和式(16)可见,当νp=0时,gp=1,表示在非界面区域不进行调整;当νp≈1时,gp≈λm/λ(如图1所示,当νp为1时,点A和点B对应的gp分别为10和8),表示在界面区域对参数Δf和kη进行相应调整可以增大界面宽度。
图1 全局修正函数gp随等效梯度νp的变化曲线(λm/λ=8和λm/λ=10)
Fig.1 Global modification function (gp) vs equivalent gradient (νp) at different width ratios (λm/λ), in which gp does not adjust any coefficients in the non-interface region (νp=0, gp=1), but modifies the coefficients in the interfacial region (see points A and B)
利用gp对Δf和kη进行调整：
式中,Δfp为第p个变体调整后的体积自由能密度差;kp为第p个变体调整后的梯度能系数。根据式(12)和(17),同时调整Δf和kη后,界面能密度γ保持不变,界面宽度λ则增至gpλ。需要指出的是,若调整后的Δfp /kp过小,由于数值计算的不稳定性,可能无法形成界面,界面宽度的增大在不同实际体系中有相应的限度。
选取Ti-50.5Ni (原子分数,%)作为模拟研究对象,研究其中B2到R相的相变过程。B2相属于高对称立方晶系;R相属于菱方晶系,有4个变体;在B2相参考坐标系下4个变体的相变应变张量分别为[23]：
模拟采用的温度为245 K;假设B2相与R相具有相同的弹性系数,C11=162 GPa,C12=129 GPa,C44=34 GPa[24]。相场方程中的各参数选取如下[25]：L=5.4774 m2/(Ns);A=0.14,B=12.42,C=12.28;Δf=1.3×107 J/m3;kη=1.1575×10-11 N。未采用全局修正函数所计算得到的相场模型界面宽度λ=0.67 nm。
采用3种相同网格数(800
图2和3分别为修正前和修正后新相在不同尺度体系中的演化过程及新相的等效半径随相变时间变化(等效半径r =
图2 新相在不同尺度体系中的演化过程和等效半径变化
Fig.2 Morphological evolution of new phase (red color) simulated in systems with scales of 134 nm×134 nm (a1~a3), 1072 nm×1072 nm (b1~b3) and 4288 nm×4288 nm (c1~c3) at time of t=0.00 s (a1~c1), t=0.15 s (a2~c2) and t=0.45 s (a3~c3), and change of the equivalent radius with time (In which each simulation system owns the same initial square nucleus size of 10.72 nm×10.72 nm) (d)
图3 引入全局修正函数后新相在不同尺度体系中的演化过程和等效半径变化
Fig.3 Morphological evolution of new phase (red color) simulated with employing the global modification function in systems with scales of 1072 nm×1072 nm (a1~a3) and 4288 nm×4288 nm (b1~b3) at time of t=0.00 s (a1, b1), t=0.15 s (a2, b2) and t=0.45 s (a3, b3), and change of the equivalent radius with time (In which each simulation system owns the same initial square nucleus size of 10.72 nm×10.72 nm) (c)
采用3种具有相同网格数(600
图4 2个R相变体在不同尺度体系中的位向关系及应变能密度变化
Fig.4 Morphological evolution of two R-phase variants in their initial states (variant 1 in blue and variant 2 in red) (a) simulated in systems with different scales of 101 nm×67 nm (b), 2010 nm×1340 nm (c) and 20100 nm×13400 nm (d), and change of strain energy density with time (e)
从图5所示引入全局修正函数(λm分别为13.4和134 nm)后2个大尺度体系的模拟结果可见,采用式(17)同时调整梯度能系数kη和体积自由能密度差Δf后,变体间界面区域的体积化学能与梯度能重新处于平衡状态,界面能不再限制变体间的界面移动,弹性应变能得以在演化过程中持续降低,促使变体协调配置,弹性应变能最终接近于0,体系呈现合理的位向关系。
图5 修正函数后2个R相变体在不同尺度体系中的位向关系及应变能密度变化
Fig.5 Morphologicl evolution of two R-phase variants (variant 1 in blue and variant 2 in red) simulated with employing the global modification function in systems with scales of 2010 nm×1340 nm (a) and 20100 nm×13400 nm (b), and change of strain energy density with time (c)
采用具有相同网格数(600
理论上,R相4个变体在B2相所在参考坐标系的二维模拟中所形成的组织[27]具有如图6所示特征。
图6 R相4个变体分别在(001)面上形成2种界面示意图
Fig.6 Schematic configuration of four R-phase variants in the (001) plane with variants 1 and 4 (a), and variants 1 and 2 (b)
图7和8分别是引入全局修正函数前后B2-R相变在不同尺度体系下的演化过程,其中蓝色代表母相B2相,浅蓝色代表R相变体1,绿色代表R相变体2,黄色代表R相变体3,红色代表R相变体4。由图7a1~a3和d可见,在随机噪声的形核阶段R相各变体都以均匀方式形核,由于变体相互挤压,体系的弹性应变能增大;随相变进行,小尺度体系中R相变体以逐渐清晰且有序的结构形成“鱼骨”状的二维组织,体系的弹性应变能降低,模拟结果与实验观察结果[27]符合得很好。当体系尺度增大后,如图7b1~b3和c1~c3所示,由于界面宽度被强制增大,致使体积化学能与梯度能在界面区域分布不合理,阻碍弹性应变能促进变体协调配置,结果R相变体无法进一步形成具有清晰位向关系的自协调组织。
图7 B2-R相变在不同尺寸体系中的演化过程及应变能密度变化
Fig.7 Morphological evolution of B2-R phase transformation simulated in systems with scales of 101 nm×101 nm (a1~a3), 2010 nm×2010 nm (b1~b3) and 20100 nm×20100 nm (c1~c3) at time of t=0.05 s (a1~c1), t=0.50 s (a2~c2), t=40.10 s (a3), t=16.52 s (b3) and t=16.13 s (c3), and change of strain energy density with time (d) (In which the blue region stands for B2-phase, the light blue region stands for variant 1, the green region stands for variant 2, the yellow region stands for variant 3, and the red region stands for variant 4)
由图8所示2种较大尺度体系中引入全局修正函数(λm分别取13.4和134 nm)后的模拟结果可见,采用式(17)同时调整梯度能系数kη和体积自由能密度差Δf对界面区域的体积化学能与梯度能进行修正后,R相得以形成具有清晰位向关系的结构,与实验观察到的组织形貌特征符合较好。
图8 修正后B2-R相变在不同尺寸体系中的演化过程及应变能密度变化
Fig.8 Morphological evolution of B2-R phase transformation simulated with employing the global modification function in systems with scales of 2010 nm×2010 nm (a1~a3) and 20100 nm×20100 nm (b1~b3) at time of t=0.05 s (a1, b1), t=0.50 s (a2, b2), t=67.71 s (a3), t=83.76 s (b3), and change of strain energy density with time (c) (In which the blue region stands for B2-phase, the light blue region stands for variant 1, the green region stands for variant 2, the yellow region stands for variant 3, and the red region stands for variant 4)
(1) 根据相场模型中界面与非界面区域序参量变化的差异,构造了能有效区分2种区域的全局修正函数,并以此改进了马氏体相变的相场模型;改进的模型可在界面区域调整体积自由能密度差Δf和梯度能系数kη,避免了影响体系整体属性的同时相当于增大了原相场模型中的界面宽度,进而扩大了马氏体相变的模拟尺度。
(2) 相比原有的马氏体相变相场模型,运用改进后的模型在相同的网格数下能够准确分析大尺度体系的界面能密度,可避免相同物理属性的新相随模拟尺度的增大而生长速率加快的问题。
(3) 应用改进的相场模型在不同尺度下对B2-R相变进行模拟,R相变体都得以形成合理的位向关系并呈现清晰有序的“鱼骨状”微观组织,与实验结果相符。
1 相场模型
1.1 模型构造
1.2 模型改进
2 模拟结果与分析
2.1 生长速率
2.2 位向关系
2.3 B2-R相变的微观组织形貌
3 结论
来源--金属学报
沪公网安备31011202020290号
