CN102750425B - 一种焊接过程热影响区组织演变的模拟方法 - Google Patents

Abstract

一种焊接过程热影响区组织演变的模拟方法，它涉及一种焊接接头微观组织模拟方法，以解决目前焊缝微观组织演变的定量化，主要基于经验或者半经验的确定性模型或者解析计算，只进行组织含量的计算，而不能动态地反映组织形态、尺寸和分布的问题。本发明的方法是通过以下步骤实现的：步骤一：计算热影响区温度场；步骤二：根据温度场分布，计算不同位置晶粒在β相区以上的晶粒长大过程；步骤三：根据温度场计算获得的冷却速度和β相晶粒分布计算结果，计算连续冷却固态相变；本发明用于焊接接头微观组织模拟。

Description

一种焊接过程热影响区组织演变的模拟方法

技术领域

本发明涉及一种焊接接头微观组织模拟方法，具体涉及一种焊接过程热影响区组织演变的模拟方法。

背景技术

焊接热影响区的晶粒长大、淬硬脆化、时效软化或回火软化、脆性相析出等问题，一直都是导致焊接裂纹或接头性能下降的薄弱环节。因此，焊接热影响区组织演化的研究具有非常现实和重要的意义，其中，数值模拟是最有效的研究手段之一。

焊接接头组织模拟是一项涉及空间尺度、时间尺度、多物理场的巨大研究工程，难度很大，以往的焊接微观组织模拟工作侧重在凝固组织和晶粒长大模拟方面，研究成果相对成熟，对于固态相变尤其是焊接热影响区组织动态演变的模拟较少，且有一些问题始终未解决，如：焊缝微观组织演变的定量化研究，目前主要基于经验或者半经验的确定性模型或者解析计算，只进行组织含量的计算，不能动态地反映组织形态、尺寸和分布；数学模型尚不完善，特别是形核模型，一方面模型中的一些参数还取决于经验公式，另一方面没有考虑焊接过程中不同冷却速度引起的形核方式的改变；由于焊接过程具有高温过热和快速冷却的特点，热影响区的相变过程偏离平衡态，然而大多研究只局限于平衡相变或者固态相变一个部分，即组织动态演变和组织历史相关性的模拟罕有报道；目前关于相变的模拟没有考虑相变晶体学特征即惯习长大现象；对于特定的合金如钛合金，尚无工程用三元或多元合金平衡相图、等温转变曲线(TTT曲线)、连续冷却转变曲线(CCT曲线)，热力学数据的困乏限制了对组织演变过程的精确描述；目前，还没有工程可用的焊接接头组织模拟可操作软件。

发明内容

本发明的目的是为了解决目前焊缝微观组织演变的定量化，主要基于经验或者半经验的确定性模型或者解析计算，只进行组织含量的计算，而不能动态地反映组织形态、尺寸和分布的问题，提供一种焊接过程热影响区组织演变的模拟方法。

本发明的焊接过程热影响区组织演变的模拟方法是通过以下步骤实现的：

步骤一：计算热影响区温度场：

通过求解导热微分方程，并且附上初始条件和换热边界条件，得到焊接过程的温度场分布：

公式一： $\mathrm{\ρc}\frac{\∂T}{\∂t}=\frac{\∂}{\∂x}\left(\mathrm{\λ}\frac{\∂T}{\∂t}\right)+\frac{\∂}{\∂y}\left(\mathrm{\λ}\frac{\∂T}{\∂t}\right)+\frac{\∂}{\∂z}\left(\mathrm{\λ}\frac{\∂T}{\∂z}\right)+\stackrel{\‾}{Q}$

其中，T为热力学温度，ρ为材料的密度，λ为材料的导热系数、c为材料的比热，为内热源强度，

将公式一计算的热影响区温度场分布，通过插值方法赋给晶粒长大的计算域，提取热影响区不同位置的热循环曲线，获得热影响区不同位置在加热阶段从低温α到高温β相变的实际温度、在冷却阶段从高温β到低温α相变的实际温度、峰值温度、晶粒长大温度区间的加热速度及冷却速度，其中，加热速度根据加热阶段的热循环曲线在晶粒长大区间的切线斜率获得，冷却速度根据冷却阶段的热循环曲线在晶粒长大区间的切线斜率获得，峰值温度即焊接粗晶区达到的最高温度；晶粒长大完毕后，开始发生高温β到低温α相变，此时冷却速度根据冷却阶段的热循环曲线在相变区间的切线斜率获得；

步骤二：根据温度场分布，计算不同位置晶粒在β相区以上的晶粒长大过程：

⑴、定义初始退火态β相组织，以及初始参数包括焊接工艺参数，热源参数，晶粒长大热力学参数；

⑵、计算晶粒长大速度，包括晶界曲率和界面能：晶界迁移速度V表示为：

公式二：V＝M_gF

其中，M_g为晶界迁移率，是衡量晶界迁移快慢的参数，F为晶粒长大驱动力，晶粒长大驱动力F由界面能和晶界的弯曲程度求解：

公式三： $F=\frac{2\mathrm{\γ}}{r}$

其中，γ为晶界的界面能，r为晶界的曲率半径，晶界长大驱动力亦即晶粒长大前后的自由能之差，

小角度晶界的界面能由相邻两晶粒间的取向差决定：

公式四： $\mathrm{\γ}={\mathrm{\γ}}_m\frac{\mathrm{\θ}}{{\mathrm{\θ}}_m}(1-1n\frac{\mathrm{\θ}}{{\mathrm{\θ}}_m})$

其中，θ为任意晶界的取向差，γ_m为界面能的最大值，θ_m为界面能对应的晶界取向差，取θ_m＝10°，

基于晶粒长大的原子跃迁机制，求得晶界迁移的速度：

公式五： $V=\frac{A_2{N}_1{\mathrm{\ν}}_1{V}_{\mathrm{\β}}^2}{N_a\mathrm{RT}}{e}^{-\frac{{\mathrm{\ΔG}}_A}{\mathrm{RT}}}\frac{2\mathrm{\γ}}{r}$

其中，晶界迁移率A₂为晶粒2接纳原子的概率，N₁为单位面积上处于能跳动的有利位置的平均原子数目，υ1为晶粒1中原子跳跃的频率，V_β为β相的摩尔体积，N_a为阿伏伽德罗常数，R为气体常数，T为热力学温度，为原子能量高于激活能的概率，△G_A为原子从晶粒1分离出来靠热激活获得的激活能势垒；

⑶、建立晶粒长大动力学条件：假设原子内能的分布符合麦克斯威尔-波尔兹曼统计规律，原子内能的计算公式：

公式六：G_T＝-RTln(RAND)

其中，RAND()表示大于0并且小于等于1的随机数，G_T为原子的内能，ln()为自然对数函数，

晶界元胞跨过晶界即长大转变的动力学能量条件为：

公式七：G_T,_i≥△G_A-∑△G_i

其中，G_T,_i为晶界原子i的内能，∑△G_i为晶界原子i向临近晶粒转变前后的自由能之差，

在上述方法中，只要存在晶界迁移驱动力，晶粒就可以长大，即晶粒长大的条件为△G>0，或者是晶界原子转变前后的自由能与内能变化之和大于0，即热力学能量条件△G+△G_T≥0，通过公式一至公式七建立了晶粒长大的动力学能量条件，即当晶界原子转变前后的内能与自由能变化，大于该原子从当前所属晶粒分离出来的跳跃激活能势垒时，晶粒才会长大，体现了原子跳跃是一个热激活过程；

⑷、计算晶粒长大概率：为了体现晶粒长大过程的不可逆性，只有达到最小能量条件的晶粒长大过程才是最稳定和最该发生的；为体现随机性，以最大晶粒长大速度为参照计算晶界原子的转变概率：

公式八： $P_i=\left\{\begin{array}{c}0,G_{T,i}<{\mathrm{\&Delta;G}}_A-{\mathrm{\&Sigma;\&Delta;G}}_i\\ \frac{M_i}{M_{\max }}\frac{F_i}{F_{\max }},G_{T,i}\&GreaterEqual;{\mathrm{\&Delta;G}}_A-{\mathrm{\&Sigma;\&Delta;G}}_i\end{array}\right.$

其中，M_i为原子i所在晶粒的晶界迁移率，F_i为原子i的晶粒长大驱动力，M_max为最大晶界迁移率，F_max为最大晶粒长大驱动力，当晶界原子可以向一个以上临近晶粒转变时，实际上只会向转变概率最大的晶粒迁移；

⑸、显示温度场分布以及焊接接头晶粒形态和大小拓扑形貌；

步骤三：根据温度场计算获得的冷却速度和β相晶粒分布计算结果，计算连续冷却固态相变：

⑴′、定义初始参数，包括溶质扩散系数、界面迁移率、冷却速度、过冷度、溶质原子在β和α相中的平衡溶解度，设置焊接热影响区晶粒长大计算结果为初始组织：

利用Thermo-Calc相图计算软件进行热力学计算，获得相变特征温度，即β和α相共存温度区间，以及溶质原子在该温度区间的溶解度，由于焊接过程快速加热和冷却，偏离热力学平衡状态，模拟β相向α相转变需考虑过冷度与特定温度、压力下的稳态平衡相比，溶质在两相中的实际溶解度按照相变平衡温度和相变过冷度计算，采用热模拟试验测量合金的CCT曲线，通过CCT曲线，获得不同冷却速度下连续冷却非平衡相变的起始温度和终止温度，将不同冷却速度的过冷度和相变温度区间、时间作为相变计算的输入条件，设置焊接热影响区晶粒长大计算结果为初始组织；

⑵′、计算相变的形核率，设置形核：

晶界形核稳态形核率I_稳为：

公式九：

其中，N_v为单位体积核胚的数量，k为波尔兹曼常数，h为普朗克常量，Q为溶质原子在β相中的热扩散激活能，△G^*为形核驱动力，

当过冷度较大或者冷却速度较大时，认为形核在相变初期瞬间完成且形核位置从晶界向晶内均匀形核转移；

⑶′、计算生长态元胞的α相长大相变驱动力、速度，计算α相前沿生长态元胞的α相转变分数、α相溶质浓度、β相溶质浓度及整个计算域的溶质场分布：

相界面移动速度V_P：

公式十：V_P＝M_PF_P/V_β

其中，M_P是相界面迁移率，FP为相变驱动力，

当以理想稀溶液处理时，相变驱动力F_P为：

公式十一： $F_P\&ap;\frac{\mathrm{RT}}{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{(X_i-X_{\mathrm{Al}}^{\mathrm{\&beta;}})}^2$

其中，X_i是界面处β相中溶质原子铝的摩尔百分数，X_Al ^β为β/α相平衡时β相的溶质浓度；

对于界面或混合控制型相变，由相长大速度获得界面处母相中铝原子的摩尔百分数：

公式十二： $X_i=X_{\mathrm{Al}}^{\mathrm{\&beta;}}+{\left(\frac{{2V}_P\&CenterDot;X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}\right)}^{\frac{1}{2}}$

β相溶质扩散方程为：

公式十三： $\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;x}\left(D_{\mathrm{\&beta;}}\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(D_{\mathrm{\&beta;}}\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;y}\right)$

α相溶质扩散方程为：

公式十四： $\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;x}\left(D_{\mathrm{\&alpha;}}\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(D_{\mathrm{\&alpha;}}\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;y}\right)$

其中，D_β为β相溶质扩散系数，D_α为α相溶质扩散系数，C_β为β相溶质浓度，C_α为α相溶质浓度，表示x方向求偏导数、表示y方向求偏导数；

根据界面处溶质守恒，获得：

公式十五：VP(XAlα-Xi)＝ $V_P(X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_i)=D(n\&CenterDot;\&dtri;C{|}_{\mathrm{interface}})$

其中，X_Al ^α为β/α相平衡态α相的溶质浓度，V_P(X_Al ^α-X_i)表示单位面积单位时间相界面处从β相脱溶析出的溶质原子量，表示界面处β相中溶质梯度引起的扩散流量，为界面处β相的溶质浓度梯度，n为β/α相界面的法向矢量；

因此，对于混合型相变，相界面移动速度V_P在x轴方向的差分表达式为：

公式十六：

$\begin{array}{c}{V}_{P,x}(i,j)=\frac{D_{\mathrm{\&beta;}}}{\mathrm{\&Delta;x}[X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M_P\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}{V}_P\right)}^{1/2}]}\left[\begin{array}{c}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M_P\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}{V}_{P,x}\right)}^{1/2}+C_{\mathrm{\&beta;}}(i-1,j)f_{\mathrm{\&beta;}}(i-1,j)\\ +X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}\right)}^{1/2}+C_{\mathrm{\&beta;}}(i+1,j)f_{\mathrm{\&beta;}}(i+1,j)\end{array}\right]\\ -\frac{D_{\mathrm{\&alpha;}}}{\mathrm{\&Delta;x}[X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}]}\left[\begin{array}{c}(-X_{\mathrm{Al}}^{\mathrm{\&alpha;}}+C_{\mathrm{\&alpha;}}(i-1,j))f_{\mathrm{\&alpha;}}(i-1,j)\\ +(-X_{\mathrm{Al}}^{\mathrm{\&alpha;}}+C_{\mathrm{\&alpha;}}(i+1,j))f_{\mathrm{\&alpha;}}(i+1,j)\end{array}\right]\end{array}$

其中V_P,_x(i,j)为计算域中坐标为(i，j)的节点在x方向的α相长大速度，△x为空间步长，V_P,_x为上一时间步长α相尖端在x方向的长大速度，C_β(i-1,j)为坐标(i-1，j)的节点的β相溶质浓度，f_β(i-1,j)为坐标(i-1，j)的节点的β相体积分数，C_β(i+1,j)为坐标(i+1，j)的节点的β相溶质浓度，f_β(i+1,j)为坐标(i+1，j)的节点的β相体积分数，C_α(i-1,j)为坐标(i-1，j)的节点的α相溶质浓度，f_α(i-1,j)为坐标(i-1，j)的节点的α相体积分数，C_α(i+1,j)为坐标(i+1，j)的节点的α相溶质浓度，f_α(i+1,j)为坐标(i+1，j)的节点的α相体积分数；

⑷′、计算α相尖端曲率，修正热力学平衡β相溶质浓度：根据Thermo-Calc软件计算可知，当热力学温度一定时，α相的平衡溶解度X_Al ^α确定，此时，当α相具有无限大曲率半径时，溶质原子在β相中的热力学平衡溶解度为X_Alβ；当α相尖端曲率为时，还需考虑β/α相界面曲率效应，即溶质原子在β相中的热力学平衡溶解度应为即将X_Al ^β修正为 $X_{\mathrm{Al}}^{\mathrm{\&beta;}}\stackrel{\&OverBar;}{\left(k\right)},$ 对于稀溶液应满足 $\ln \frac{X_{\mathrm{Al}}^{\mathrm{\&beta;}}\stackrel{\&OverBar;}{\left(k\right)}}{X_{\mathrm{Al}}^{\mathrm{\&beta;}}}=\mathrm{\&Gamma;}\stackrel{\&OverBar;}{k},$ 由此：

公式十七： $X_{\mathrm{Al}}^{\mathrm{\&beta;}}\stackrel{\&OverBar;}{\left(k\right)}=X_{\mathrm{Al}}^{\mathrm{\&beta;}}(1+\mathrm{\&Gamma;}\stackrel{\&OverBar;}{k})$

其中，为α相尖端曲率，Γ为Gibbs-Thomson系数，

α相尖端曲率：

公式十八： $\stackrel{\&OverBar;}{k}=\frac{1}{\mathrm{\&Delta;x}}(1-2\frac{f_{\mathrm{\&alpha;}}+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_{\mathrm{\&alpha;}}\left(k\right)}{N+1})$

其中，f_α生长态元胞的α相分数，N为元胞自动机的邻居元胞数目，对于冯·诺依曼型邻居关系N＝4，对于摩尔型邻居关系型邻居关系N＝8；

Gibbs-Thomson系数：

公式十九： $\mathrm{\&Gamma;}=\frac{{2\mathrm{\&sigma;V}}_{\mathrm{\&beta;}}}{\mathrm{RT}(X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}})}$

其中，σ为β/α的相界界面能；

生长态元胞α相分数：

公式二十： ${\mathrm{\&delta;f}}_{\mathrm{\&alpha;}}=\frac{\mathrm{\&delta;t}}{\mathrm{\&Delta;x}}(V_x+V_y-V_x{V}_y\frac{\mathrm{\&delta;t}}{\mathrm{\&Delta;x}})$ 并且f_α ⁿ＝f_α ⁰+δf_α

其中，δt为计算时间步长，Vx为生长态元胞在x方向的相界面移动速度，V_y为生长态元胞在y方向的相界面移动速度，f_α ⁰为上一时间步长的生长态元胞的α相分数，f_α ⁿ分别为当前时间步长的生长态元胞的α相分数，δf_α为生长态元胞的α相分数增量；

⑸′、计算各向异性相界面迁移率：钛合金从β相区连续冷却时，α相与β相基体之间存在严格的伯格斯(Burgers)晶体学取向关系，以界面迁移率的各向异性表征，在择优生长方向和与其垂直方向的界面迁移率分别记为M₁和M₂，根据相界迁移的微观物理机制，获得界面迁移率的表达式：

公式二十一：M＝M₀exp(-Q_m/RT)

公式二十二： $M_0=\frac{1}{16}{(V_{\mathrm{\&alpha;}}^{\frac{1}{3}}+V_{\mathrm{\&beta;}}^{\frac{1}{3}})}^4{V}_D{\left(\mathrm{kT}\right)}^{-1}$

其中，M₀为指数前因子，V_α为α相的摩尔体积、ν_D≈5×10¹³/s为德拜频率；由于M₁/M₂的大小直接决定析出新相的长细比，取M₁/M₂＝10、100或1000，分别进行析出α相形态的模拟，并且通过扫描电镜观测相同焊接工艺条件下热影响区的组织形态，将析出α相形态的模拟结果与实验结果对比，两者最吻合时采用的M₁/M₂最佳；

⑹′、计算晶内不同择优取向的α相长大形貌：新相长大方向与母相晶粒的晶体位向有关，当母相晶体取向为优先取向时，按照与β相晶粒的晶体取向一致方向长大，当β相晶粒的晶体取向不是优先取向时，则与β相晶粒的晶体取向呈一特定角度长大，由于对于钛合金而言，β-α相变存在多个惯习面和惯习方向，所以每个β相晶粒内α相长大的位相取向并不唯一，因此，α相生长方向与β相的晶界呈一定角度，最终经历连续冷却非平衡相变后，在焊接热影响区形成网篮组织；

⑺′、更新元胞状态，显示组织形态及溶质场分布、不同位置不同时间的α相变分数，并以文件形式保存。

本发明具有以下有益效果：

一、本发明的方法是利用计算机在虚拟条件下进行的。采用本发明的焊接热影响区晶粒长大和固态相变组织演变元胞自动机模型，实现了焊接热影响区不同位置尤其是粗晶区的组织演变模拟，能够动态再现焊接全过程的组织形态，定量预测焊接接头最终晶粒尺寸、相含量、组成相的特征尺寸，如片状α相的宽度，α相集束的宽度等，从而建立了工艺-组织与力学性能之间的映射关系，为制定钛合金的焊接工艺，改善焊接接头的显微组织和力学性能提供理论基础，因此，本发明解决了目前焊缝微观组织演变的定量化，主要基于经验或者半经验的确定性模型或者解析计算，只进行组织含量的计算，而不能动态地反映组织形态、尺寸和分布的问题。

二、由于本发明采用元胞自动机法能较精确地模拟材料的介观和微观组织，从而实现固态相变元胞自动机法模拟与介观尺度上晶粒长大元胞自动机法模拟的“无缝连接”，并以模拟晶粒长大为桥梁建立与宏观温度场之间的联系，实现宏观-介观-微观多尺度耦合模拟，揭示焊接过程组织遗传规律，即焊接热影响区晶粒过热粗大对后续接头组织演化规律的影响。

三、本发明通过数值模拟、相变理论和材料热力学多学科结合，在不同尺度和不同位置上再现焊接接头组织演化的形貌，定量描述晶体学、热力学和动力学机制，并且可预测TTT曲线，填充了航空新材料钛合金TTT/CCT图的空白。

四、组织场、宏观应力应变场和温度场三者相互作用，共同影响焊接接头性能，通过组织模拟能定量计算溶质偏析，并能建立微观局部应力场和组织场的耦合模型，进而结合断裂力学和损伤力学，分析产生焊接缺陷如冷裂纹产生的冶金条件和力学条件。

附图说明

图1为具体实施方式一中焊接热影响区组织演变的计算方法的流程图；图2为焊接热影响区温度场分布的模拟结果图；图3为焊接热影响区距离焊缝不同位置的热循环曲线图；图4为具体实施方式一中步骤二的流程图；图5为具体实施方式一中步骤三的流程图；图6为TA15钛合金在焊接热影响区晶粒长大和固态相变组织形貌模拟结果图；图7为TA15钛合金在焊接热影响区不同位置随时间变化的最终相变分数计算结果图；图8为TA15钛合金在焊接热影响区不同大小的晶粒内部不同角度的α相长大形貌图。

具体实施方式

具体实施方式一：结合图1～图5说明本实施方式，本发明的焊接过程热影响区组织演变的模拟方法是通过以下步骤实现的：

步骤一：计算热影响区温度场：

通过求解导热微分方程，并且附上初始条件和换热边界条件，得到焊接过程的温度场分布：

公式一： $\mathrm{\&rho;c}\frac{\&PartialD;T}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;x}\left(\mathrm{\&lambda;}\frac{\&PartialD;T}{\&PartialD;t}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(\mathrm{\&lambda;}\frac{\&PartialD;T}{\&PartialD;t}\right)+\frac{\&PartialD;}{\&PartialD;z}\left(\mathrm{\&lambda;}\frac{\&PartialD;T}{\&PartialD;z}\right)+\stackrel{\&OverBar;}{Q}$

其中，T为热力学温度，ρ为材料的密度，λ为材料的导热系数、c为材料的比热，为内热源强度，

将公式一计算的热影响区温度场分布，如图2所示，通过插值方法赋给晶粒长大的计算域，提取热影响区不同位置的热循环曲线，如图3所示，获得热影响区不同位置在加热阶段从低温α到高温β相变的实际温度、在冷却阶段从高温β到低温α相变的实际温度、峰值温度、晶粒长大温度区间的加热速度及冷却速度，其中，加热速度根据加热阶段的热循环曲线在晶粒长大区间的切线斜率获得，冷却速度根据冷却阶段的热循环曲线在晶粒长大区间的切线斜率获得，峰值温度即焊接粗晶区达到的最高温度；晶粒长大完毕后，开始发生高温β到低温α相变，此时冷却速度根据冷却阶段的热循环曲线在相变区间的切线斜率获得；

步骤二：根据温度场分布，计算不同位置晶粒在β相区以上的晶粒长大过程：如图4所示，

⑴、定义初始退火态β相组织，以及初始参数包括焊接工艺参数，热源参数，晶粒长大热力学参数；

⑵、计算晶粒长大速度，包括晶界曲率和界面能：晶界迁移速度V表示为：

公式二：V＝M_gF

其中，M_g为晶界迁移率，是衡量晶界迁移快慢的参数，F为晶粒长大驱动力，

晶粒长大驱动力F由界面能和晶界的弯曲程度求解：

公式三： $F=\frac{2\mathrm{\&gamma;}}{r}$

其中，γ为晶界的界面能，r为晶界的曲率半径，晶界长大驱动力亦即晶粒长大前后的自由能之差，

小角度晶界的界面能由相邻两晶粒间的取向差决定：

公式四： $\mathrm{\&gamma;}={\mathrm{\&gamma;}}_m\frac{\mathrm{\&theta;}}{{\mathrm{\&theta;}}_m}(1-1n\frac{\mathrm{\&theta;}}{{\mathrm{\&theta;}}_m})$

其中，θ为任意晶界的取向差，γ_m为界面能的最大值，θ_m为界面能对应的晶界取向差，取θ_m＝10°，

基于晶粒长大的原子跃迁机制，求得晶界迁移的速度：

公式五： $V=\frac{A_2{N}_1{\mathrm{\&nu;}}_1{V}_{\mathrm{\&beta;}}^2}{N_a\mathrm{RT}}{e}^{-\frac{{\mathrm{\&Delta;G}}_A}{\mathrm{RT}}}\frac{2\mathrm{\&gamma;}}{r}$

其中，晶界迁移率A₂为晶粒2接纳原子的概率，N₁为单位面积上处于能跳动的有利位置的平均原子数目，υ₁为晶粒1中原子跳跃的频率，V_β为β相的摩尔体积，N_a为阿伏伽德罗常数，R为气体常数，T为热力学温度，为原子能量高于激活能的概率，△G_A为原子从晶粒1分离出来靠热激活获得的激活能势垒；

⑶、建立晶粒长大条件：假设原子内能的分布符合麦克斯威尔-波尔兹曼统计规律，原子内能的计算公式：

公式六：G_T＝-RTln(RAND)

其中，RAND()表示大于0并且小于等于1的随机数，G_T为原子的内能，ln()为自然对数函数，

晶界元胞跨过晶界即长大转变的动力学能量条件为：

公式七：G_T,_i≥△G_A-∑△G_i

其中，G_T,_i为晶界原子i的内能，∑△G_i为晶界原子i向临近晶粒转变前后的自由能之差，

在上述方法中，只要存在晶界迁移驱动力，晶粒就可以长大，即晶粒长大的条件为△G>0，或者是晶界原子转变前后的自由能与内能变化之和大于0，即热力学能量条件△G+△G_T≥0，通过公式一至公式七建立了晶粒长大的动力学能量条件，即当晶界原子转变前后的内能与自由能变化，大于该原子从当前所属晶粒分离出来的跳跃激活能势垒时，晶粒才会长大，体现了原子跳跃是一个热激活过程；

⑷、计算晶粒长大概率：为了体现晶粒长大过程的不可逆性，只有达到最小能量条件的晶粒长大过程才是最稳定和最该发生的；为体现随机性，以最大晶粒长大速度为参照计算晶界原子的转变概率：

公式八： $P_i=\left\{\begin{array}{c}0,G_{T,i}<{\mathrm{\&Delta;G}}_A-{\mathrm{\&Sigma;\&Delta;G}}_i\\ \frac{M_i}{M_{\max }}\frac{F_i}{F_{\max }},G_{T,i}\&GreaterEqual;{\mathrm{\&Delta;G}}_A-{\mathrm{\&Sigma;\&Delta;G}}_i\end{array}\right.$

其中，M_i为原子i所在晶粒的晶界迁移率，F_i为原子i的晶粒长大驱动力，M_max为最大晶界迁移率，F_max为最大晶粒长大驱动力，当晶界原子可以向一个以上临近晶粒转变时，实际上只会向转变概率最大的晶粒迁移；

⑸、显示温度场分布以及焊接接头晶粒形态和大小拓扑形貌；

步骤三：根据温度场计算获得的冷却速度和β相晶粒分布计算结果，计算连续冷却固态相变：如图5所示，

⑴′、定义初始参数，包括溶质扩散系数、界面迁移率、冷却速度、过冷度、溶质原子在β和α相中的平衡溶解度，设置焊接热影响区晶粒长大计算结果为初始组织：

利用Thermo-Calc相图计算软件进行热力学计算，获得相变特征温度，即β和α相共存温度区间，以及溶质原子在该温度区间的溶解度，由于焊接过程快速加热和冷却，偏离热力学平衡状态，模拟β相向α相转变需考虑过冷度与特定温度、压力下的稳态平衡相比，溶质在两相中的实际溶解度按照相变平衡温度和相变过冷度计算，采用热模拟试验测量合金的CCT曲线，通过CCT曲线，获得不同冷却速度下连续冷却非平衡相变的起始温度和终止温度，将不同冷却速度的过冷度和相变温度区间、时间作为相变计算的输入条件，设置焊接热影响区晶粒长大计算结果为初始组织；

⑵′、计算相变的形核率，设置形核：

晶界形核稳态形核率I_稳为：

公式九：

其中，N_v为单位体积核胚的数量，k为波尔兹曼常数，h为普朗克常量，Q为溶质原子在β相中的热扩散激活能，△G^*为形核驱动力，

当过冷度较大或者冷却速度较大时，认为形核在相变初期瞬间完成且形核位置从晶界向晶内均匀形核转移；

⑶′、计算生长态元胞的α相长大相变驱动力、速度，计算α相前沿生长态元胞的α相转变分数、α相溶质浓度、β相溶质浓度及整个计算域的溶质场分布：

相界面移动速度V_P：

公式十：V_P＝M_PF_P/V_β

其中，M_P是相界面迁移率，FP为相变驱动力，

当以理想稀溶液处理时，相变驱动力F_P为：

公式十一： $F_P\&ap;\frac{\mathrm{RT}}{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{(X_i-X_{\mathrm{Al}}^{\mathrm{\&beta;}})}^2$

其中，X_i是界面处β相中溶质原子铝的摩尔百分数，X_Al ^β为β/α相平衡时β相的溶质浓度；

对于界面或混合控制型相变，由相长大速度获得界面处母相中铝原子的摩尔百分数：

公式十二： $X_i=X_{\mathrm{Al}}^{\mathrm{\&beta;}}+{\left(\frac{{2V}_P\&CenterDot;X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}\right)}^{\frac{1}{2}}$

β相溶质扩散方程为：

公式十三： $\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;x}\left(D_{\mathrm{\&beta;}}\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(D_{\mathrm{\&beta;}}\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;y}\right)$

α相溶质扩散方程为：

公式十四： $\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;x}\left(D_{\mathrm{\&alpha;}}\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(D_{\mathrm{\&alpha;}}\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;y}\right)$

其中，D_β为β相溶质扩散系数，D_α为α相溶质扩散系数，C_β为β相溶质浓度，C_α为α相溶质浓度，表示x方向求偏导数、表示y方向求偏导数；

根据界面处溶质守恒，获得：

公式十五：VP(XAlα-Xi)＝ $V_P(X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_i)=D(n\&CenterDot;\&dtri;C{|}_{\mathrm{interface}})$

其中，X_Al ^α为β/α相平衡态α相的溶质浓度，V_P(X_Al ^α-X_i)表示单位面积单位时间相界面处从β相脱溶析出的溶质原子量，表示界面处β相中溶质梯度引起的扩散流量，为界面处β相的溶质浓度梯度，n为β/α相界面的法向矢量；

因此，对于混合型相变，相界面移动速度V_P在x轴方向的差分表达式为：

公式十六：

$\begin{array}{c}{V}_{P,x}(i,j)=\frac{D_{\mathrm{\&beta;}}}{\mathrm{\&Delta;x}[X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M_P\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}{V}_P\right)}^{1/2}]}\left[\begin{array}{c}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M_P\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}{V}_{P,x}\right)}^{1/2}+C_{\mathrm{\&beta;}}(i-1,j)f_{\mathrm{\&beta;}}(i-1,j)\\ +X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}\right)}^{1/2}+C_{\mathrm{\&beta;}}(i+1,j)f_{\mathrm{\&beta;}}(i+1,j)\end{array}\right]\\ -\frac{D_{\mathrm{\&alpha;}}}{\mathrm{\&Delta;x}[X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}]}\left[\begin{array}{c}(-X_{\mathrm{Al}}^{\mathrm{\&alpha;}}+C_{\mathrm{\&alpha;}}(i-1,j))f_{\mathrm{\&alpha;}}(i-1,j)\\ +(-X_{\mathrm{Al}}^{\mathrm{\&alpha;}}+C_{\mathrm{\&alpha;}}(i+1,j))f_{\mathrm{\&alpha;}}(i+1,j)\end{array}\right]\end{array}$

其中V_P,_x(i,j)为计算域中坐标为(i，j)的节点在x方向的α相长大速度，△x为空间步长，V_P,_x为上一时间步长α相尖端在x方向的长大速度，C_β(i-1,j)为坐标(i-1，j)的节点的β相溶质浓度，f_β(i-1,j)为坐标(i-1，j)的节点的β相体积分数，C_β(i+1,j)为坐标(i+1，j)的节点的β相溶质浓度，f_β(i+1,j)为坐标(i+1，j)的节点的β相体积分数，C_α(i-1,j)为坐标(i-1，j)的节点的α相溶质浓度，f_α(i-1,j)为坐标(i-1，j)的节点的α相体积分数，C_α(i+1,j)为坐标(i+1，j)的节点的α相溶质浓度，f_α(i+1,j)为坐标(i+1，j)的节点的α相体积分数；

⑷′、计算α相尖端曲率，修正热力学平衡β相溶质浓度：根据Thermo-Calc软件计算可知，当热力学温度一定时，α相的平衡溶解度X_Al ^α确定，此时，当α相具有无限大曲率半径时，溶质原子在β相中的热力学平衡溶解度为X_Alβ；当α相尖端曲率为时，还需考虑β/α相界面曲率效应，即溶质原子在β相中的热力学平衡溶解度应为即将X_Al ^β修正为 $X_{\mathrm{Al}}^{\mathrm{\&beta;}}\stackrel{\&OverBar;}{\left(k\right)},$ 对于稀溶液应满足 $\ln \frac{X_{\mathrm{Al}}^{\mathrm{\&beta;}}\stackrel{\&OverBar;}{\left(k\right)}}{X_{\mathrm{Al}}^{\mathrm{\&beta;}}}=\mathrm{\&Gamma;}\stackrel{\&OverBar;}{k},$ 由此：

公式十七： $X_{\mathrm{Al}}^{\mathrm{\&beta;}}\stackrel{\&OverBar;}{\left(k\right)}=X_{\mathrm{Al}}^{\mathrm{\&beta;}}(1+\mathrm{\&Gamma;}\stackrel{\&OverBar;}{k})$

其中，为α相尖端曲率，Γ为Gibbs-Thomson系数，

α相尖端曲率：

公式十八： $\stackrel{\&OverBar;}{k}=\frac{1}{\mathrm{\&Delta;x}}(1-2\frac{f_{\mathrm{\&alpha;}}+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_{\mathrm{\&alpha;}}\left(k\right)}{N+1})$

其中，f_α生长态元胞的α相分数，N为元胞自动机的邻居元胞数目，对于冯·诺依曼型邻居关系N＝4，对于摩尔型邻居关系型邻居关系N＝8；

Gibbs-Thomson系数：

公式十九： $\mathrm{\&Gamma;}=\frac{{2\mathrm{\&sigma;V}}_{\mathrm{\&beta;}}}{\mathrm{RT}(X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}})}$

其中，σ为β/α的相界界面能；

生长态元胞α相分数：

公式二十： ${\mathrm{\&delta;f}}_{\mathrm{\&alpha;}}=\frac{\mathrm{\&delta;t}}{\mathrm{\&Delta;x}}(V_x+V_y-V_x{V}_y\frac{\mathrm{\&delta;t}}{\mathrm{\&Delta;x}})$ 并且f_α ⁿ＝f_α ⁰+δf_α

其中，δt为计算时间步长，V_x为生长态元胞在x方向的相界面移动速度，V_y为生长态元胞在y方向的相界面移动速度，f_α ⁰为上一时间步长的生长态元胞的α相分数，f_α ⁿ分别为当前时间步长的生长态元胞的α相分数，δf_α为生长态元胞的α相分数增量；

⑸′、计算各向异性相界面迁移率：钛合金从β相区连续冷却时，α相与β相基体之间存在严格的伯格斯(Burgers)晶体学取向关系，以界面迁移率的各向异性表征，在择优生长方向和与其垂直方向的界面迁移率分别记为M₁和M₂，根据相界迁移的微观物理机制，获得界面迁移率的表达式：

公式二十一：M＝M₀exp(-Q_m/RT)

公式二十二： $M_0=\frac{1}{16}{(V_{\mathrm{\&alpha;}}^{\frac{1}{3}}+V_{\mathrm{\&beta;}}^{\frac{1}{3}})}^4{V}_D{\left(\mathrm{kT}\right)}^{-1}$

其中，M₀为指数前因子，V_α为α相的摩尔体积、ν_D≈5×10¹³/s为德拜频率；由于M₁/M₂的大小直接决定析出新相的长细比，取M₁/M₂＝10、100或1000，分别进行析出α相形态的模拟，并且通过扫描电镜观测相同焊接工艺条件下热影响区的组织形态，将析出α相形态的模拟结果与实验结果对比，两者最吻合时采用的M₁/M₂最佳；

⑹′、计算晶内不同择优取向的α相长大形貌：新相长大方向与母相晶粒的晶体位向有关，当母相晶体取向为优先取向时，按照与β相晶粒的晶体取向一致方向长大，当β相晶粒的晶体取向不是优先取向时，则与β相晶粒的晶体取向呈一特定角度长大，由于对于钛合金而言，β-α相变存在多个惯习面和惯习方向，所以每个β相晶粒内α相长大的位相取向并不唯一，因此，α相生长方向与β相的晶界呈一定角度，最终经历连续冷却非平衡相变后，在焊接热影响区形成网篮组织；

⑺′、更新元胞状态，显示组织形态及溶质场分布、不同位置不同时间的α相变分数，并以文件形式保存。

不同角度α相的长大形貌的计算按照占小红提出的对角线模拟法则进行，(参见《Cellular automaton simulation of grain growth with different orientation angles duringsolidification process》：X.H.Zhan,Y.H.Wei,Z.B.Dong.Journal of Materials ProcessingTechnology.Vol.12(2007):130)。

本发明建立了在焊接工艺下钛合金的非平衡连续冷却相变模型，预测热影响区晶粒长大和粗大β相晶粒在冷却时经历相变过程的最终相变分数以及显微组织状态，并且定量表征冷却速度、界面能各向异性对组织形貌、新相长大动力学和组织转变分数的影响规律，为制定合理的钛合金焊接工艺提供理论指导。

对于建立的晶粒长大和固态相变计算模型，现以目前应用较广泛的TA15钛合金为例，进行CCT曲线测定实验和相关热力学计算和组织演变模拟，该TA15钛合金的材料成分原子百分含量如下表一所示。

表一TA15的化学成分(mole.％)

其热影响区固态相变数值模拟的方法通过一下步骤实现：

步骤a、针对TA15钛合金，采用最新版本的Thermo-Calc软件开展热力学计算，为简便，以Al当量和Mo当量为参考，将TA15钛合金近似为Ti-Al-Mo三元系合金，采用正规稀溶体近似描述β相和α相，结合材料热力学理论，获得Ti-Al-Mo三元系Al和Mo在母相和新相溶解度随着温度的变化，两相共存温度区间为980℃-774℃，由于焊接过程快速加热和冷却，偏离热力学平衡状态，模拟β相向α相转变需考虑过冷度，即焊接条件下系统处于亚稳平衡态，与特定温度、压力下的稳态平衡相比，平衡成分范围或温度都发生了明显的变化；

步骤b、采用Gleeble1500热模拟试验测量TA15钛合金的CCT曲线，通过CCT曲线，定量的计算连续冷却时相变孕育时间，相变开始时间和温度、结束时间和温度，以及最终转变产物的各相体积分数比，并将热循环试验的热循环曲线作为模拟条件，从而提高模拟结果的可信度。试验棒材的规格为φ6mm×130mm(可接受长度在100mm～150mm之间，粗糙度为Ra3.2，径向正偏差)；热循环工艺如下：加热速度为90℃/s、峰值温度为1360℃，所以β相晶粒粗大，冷却速度依次为18.3℃/s、16.7℃/s、9.9℃/s、5.0℃/s、2.5℃/s、1.5℃/s。整个试验过程在氩气保护或者真空状态下进行，以防止高温下试样被氧化，并由铂金-铂锗热电偶插入预先激光加工的φ1mm×2mm深孔内测量加热区温度变化，测量得到不同冷却速度时试验棒料径向收缩量随温度的变化，通过切线法获得相变起始和终止温度。

步骤c、计算β-α相变形核率。晶界形核稳态形核率为：

其中，N_v为单位体积核胚(原子团)的数量。

计算表明，当相变温度在750℃～900℃之间时，形核率在0(nuclei/m/s)～0.01(nuclei/m/s)之间，与温度呈抛物线关系变化，成分固定时，在850℃左右形核率达最大值，温度过高或过低都会造成形核率的下降，形核方式即是均匀形核还是晶界形核，以及是瞬间形核还是连续形核，依据过冷度和冷却速度大小而定；当发生连续冷却相变时，与等温相变形核率计算不同，近似认为形核率与冷却速度呈线性关系变化，由于相变过冷度大，认为形核在相变初期瞬间完成。

步骤d、计算相变速度和溶质场分布：

相界面移动速度V_P：V_P＝M_PF_P/V_β

其中，MP是相界面迁移率，FP为相变驱动力。

当以理想稀溶液近似处理时，相变驱动力F_P为： $F_P\&ap;\frac{\mathrm{RT}}{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{(X_i-X_{\mathrm{Al}}^{\mathrm{\&beta;}})}^2$

其中，Xi是界面处母相中铝原子的摩尔百分数，X_Al ^β为母/新相平衡态母相的溶质浓度；

对于界面或混合控制型相变，由相长大速度获得界面处β相中铝原子的摩尔百分数：

$X_i=X_{\mathrm{Al}}^{\mathrm{\&beta;}}+{\left(\frac{{2V}_P\&CenterDot;X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}\right)}^{\frac{1}{2}}$

β相溶质扩散方程： $\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;x}\left(D_{\mathrm{\&beta;}}\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(D_{\mathrm{\&beta;}}\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;y}\right)$

α相溶质扩散方程： $\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;x}\left(D_{\mathrm{\&alpha;}}\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(D_{\mathrm{\&alpha;}}\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;y}\right)$

其中，D_β、D_α分别为β相和α相溶质扩散系数，C_β和C_α为β相和α相的溶质浓度。

根据界面处溶质守恒，获得：V(XAlα-Xi)＝ $V(X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_i)=D(n\&CenterDot;\&dtri;C{|}_{\mathrm{interface}})$

其中，X_Al ^β为β/α相平衡时β相的溶质浓度，等式左侧表示单位面积单位时间相界面处从β相脱溶析出的溶质原子量，为界面处β相中溶质梯度引起的扩散流量，为β相的溶质浓度梯度，n为β/α相界面的法向矢量。

因此，对于混合型相变，速度V_P在x轴方向的差分表达式为：

$\begin{array}{c}{V}_{P,x}(i,j)=\frac{D_{\mathrm{\&beta;}}}{\mathrm{\&Delta;x}[X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M_P\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}{V}_P\right)}^{1/2}]}\left[\begin{array}{c}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M_P\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}{V}_{P,x}\right)}^{1/2}+C_{\mathrm{\&beta;}}(i-1,j)f_{\mathrm{\&beta;}}(i-1,j)\\ +X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}\right)}^{1/2}+C_{\mathrm{\&beta;}}(i+1,j)f_{\mathrm{\&beta;}}(i+1,j)\end{array}\right]\\ -\frac{D_{\mathrm{\&alpha;}}}{\mathrm{\&Delta;x}[X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}]}\left[\begin{array}{c}(-X_{\mathrm{Al}}^{\mathrm{\&alpha;}}+C_{\mathrm{\&alpha;}}(i-1,j))f_{\mathrm{\&alpha;}}(i-1,j)\\ +(-X_{\mathrm{Al}}^{\mathrm{\&alpha;}}+C_{\mathrm{\&alpha;}}(i+1,j))f_{\mathrm{\&alpha;}}(i+1,j)\end{array}\right]\end{array}$

步骤e、计算α相尖端曲率，修正热力学平衡β相溶质浓度：考虑β/α相界面曲率效应后，溶质原子在β相中的热力学平衡溶解度应为即将X_Al ^β修正为对于稀溶液 $X_{\mathrm{Al}}^{\mathrm{\&beta;}}\stackrel{\&OverBar;}{\left(k\right)}=X_{\mathrm{Al}}^{\mathrm{\&beta;}}(1+\mathrm{\&Gamma;}\stackrel{\&OverBar;}{k})$

α相尖端尖端曲率： $\stackrel{\&OverBar;}{k}=\frac{1}{\mathrm{\&Delta;x}}(1-2\frac{f_{\mathrm{\&alpha;}}+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_{\mathrm{\&alpha;}}\left(k\right)}{N+1})$

Gibbs-Thomson系数： $\mathrm{\&Gamma;}=\frac{{2\mathrm{\&sigma;V}}_{\mathrm{\&beta;}}}{\mathrm{RT}(X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}})}$

生长态元胞α相分数： ${\mathrm{\&delta;f}}_{\mathrm{\&alpha;}}=\frac{\mathrm{\&delta;t}}{\mathrm{\&Delta;x}}(V_x+V_y-V_x{V}_y\frac{\mathrm{\&delta;t}}{\mathrm{\&Delta;x}})$ 并且f_α ⁿ＝f_α ⁰+δf_α

步骤f、计算各向异性相界面迁移率：钛合金从β相区连续冷却时，α相与β相基体之间存在严格的伯格斯(Burgers)晶体学取向关系，因此，以界面迁移率的各向异性表征，在择优生长方向和与其垂直方向的界面迁移率分别记为M₁和M₂，根据相界迁移的微观物理机制，获得相界面迁移率的表达式：M_P＝M₀exp(-Q/RT)

其中，为指数前因子，V_α为α相的摩尔体积，ν_D≈5×10¹³/s为德拜频率；由于M₁/M₂的大小直接决定析出新相的长细比，取M₁/M₂＝10、100或1000，分别进行析出α相形态的模拟，并且通过扫描电镜观测相同焊接工艺条件下热影响区的组织形态，将析出α相形态的模拟结果与实验结果对比，两者最吻合时采用的M₁/M₂最佳，为M₁/M₂＝100。随着冷却速度增大，在以上关于界面迁移率的表达式基础上，近似认为相界面迁移率与冷却速度呈线性正比关系。

步骤g、根据TIG焊接工艺下(如表二所示)的温度场和热循环模拟结果，开展焊接热影响区组织演变计算，计算域为200×150μm，空间网格数为400×300和空间步长为0.5μm。焊接热影响区晶粒长大和固态相变组织形貌模拟组织形貌模拟结果如图6所示，焊接热影响区不同位置随时间变化的最终相变分数如图7所示；

表二TIG焊接工艺参数

步骤h、计算晶内不同择优取向的α相长大形貌：固态相变需遵循晶体学位向关系，α相从β相析出长大时，按照与惯习面平行的方向优先选择生长，所以，α相长大方向与β相晶粒的晶体位向有关，当β相晶体取向为优先取向时，按照与β相晶体取向一致方向长大；当β相晶体取向不是优先取向时则与β相晶体取向为一特定角度长大；在模拟时，当设置晶界方向固定为x轴方向，在任意晶粒内的α相长大，即在面心立方体其它晶向的β相内长大，长大方向与x轴呈一定夹角；由于存在多个惯习面和方向，所以每个晶粒内α相长大的取向并不唯一，形成网篮组织，如图8所示。

Claims (1)

1.一种焊接过程热影响区组织演变的模拟方法，其特征在于：所述方法是通过以下步骤实现的：

步骤一：计算热影响区温度场：

通过求解导热微分方程，并且附上初始条件和换热边界条件，得到焊接过程的温度场分布：

公式一： $\mathrm{\&rho;c}\frac{\&PartialD;T}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;x}\left(\mathrm{\&lambda;}\frac{\&PartialD;T}{\&PartialD;t}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(\mathrm{\&lambda;}\frac{\&PartialD;T}{\&PartialD;t}\right)+\frac{\&PartialD;}{\&PartialD;z}\left(\mathrm{\&lambda;}\frac{\&PartialD;T}{\&PartialD;z}\right)+\stackrel{\&OverBar;}{Q}$

其中，T为热力学温度，ρ为材料的密度，λ为材料的导热系数、c为材料的比热，为内热源强度，

将公式一计算的热影响区温度场分布，通过插值方法赋给晶粒长大的计算域，提取热影响区不同位置的热循环曲线，获得热影响区不同位置在加热阶段从低温α到高温β相变的实际温度、在冷却阶段从高温β到低温α相变的实际温度、峰值温度、晶粒长大温度区间的加热速度及冷却速度，其中，加热速度根据加热阶段的热循环曲线在晶粒长大区间的切线斜率获得，冷却速度根据冷却阶段的热循环曲线在晶粒长大区间的切线斜率获得，峰值温度即焊接粗晶区达到的最高温度；晶粒长大完毕后，开始发生高温β到低温α相变，此时冷却速度根据冷却阶段的热循环曲线在相变区间的切线斜率获得；

步骤二：根据温度场分布，计算不同位置晶粒在β相区以上的晶粒长大过程：

⑴、定义初始退火态β相组织，以及初始参数包括焊接工艺参数，热源参数，晶粒长大热力学参数；

⑵、计算晶粒长大速度，包括晶界曲率和界面能：晶界迁移速度V表示为：

公式二：V＝M_gF

其中，M_g为晶界迁移率，是衡量晶界迁移快慢的参数，F为晶粒长大驱动力，晶粒长大驱动力F由界面能和晶界的弯曲程度求解：

公式三： $F=\frac{2\mathrm{\&gamma;}}{r}$

其中，γ为晶界的界面能，r为晶界的曲率半径，晶界长大驱动力亦即晶粒长大前后的自由能之差，

小角度晶界的界面能由相邻两晶粒间的取向差决定：

公式四： $\mathrm{\&gamma;}={\mathrm{\&gamma;}}_m\frac{\mathrm{\&theta;}}{{\mathrm{\&theta;}}_m}(1-1n\frac{\mathrm{\&theta;}}{{\mathrm{\&theta;}}_m})$

其中，θ为任意晶界的取向差，γ_m为界面能的最大值，θ_m为界面能对应的晶界取向差，取θ_m＝10°，

基于晶粒长大的原子跃迁机制，求得晶界迁移的速度：

公式五： $V=\frac{A_2{N}_1{\mathrm{\&nu;}}_1{V}_{\mathrm{\&beta;}}^2}{N_a\mathrm{RT}}{e}^{-\frac{{\mathrm{\&Delta;G}}_A}{\mathrm{RT}}}\frac{2\mathrm{\&gamma;}}{r}$

其中，晶界迁移率A₂为晶粒2接纳原子的概率，N₁为单位面积上处于能跳动的有利位置的平均原子数目，υ1为晶粒1中原子跳跃的频率，V_β为β相的摩尔体积，N_a为阿伏伽德罗常数，R为气体常数，T为热力学温度，为原子能量高于激活能的概率，△G_A为原子从晶粒1分离出来靠热激活获得的激活能势垒；

⑶、建立晶粒长大动力学条件：假设原子内能的分布符合麦克斯威尔-波尔兹曼统计规律，原子内能的计算公式：

公式六：G_T＝-RTln(RAND)

其中，RAND()表示大于0并且小于等于1的随机数，G_T为原子的内能，ln()为自然对数函数，

晶界元胞跨过晶界即长大转变的动力学能量条件为：

公式七：G_T,_i≥△G_A-∑△G_i

其中，G_T,_i为晶界原子i的内能，∑△G_i为晶界原子i向临近晶粒转变前后的自由能之差，

在上述步骤中，只要存在晶界迁移驱动力，晶粒就可以长大，即晶粒长大的条件为△G>0，或者是晶界原子转变前后的自由能与内能变化之和大于0，即热力学能量条件△G+△G_T≥0，通过公式一至公式七建立了晶粒长大的动力学能量条件，即当晶界原子转变前后的内能与自由能变化，大于该原子从当前所属晶粒分离出来的跳跃激活能势垒时，晶粒才会长大，体现了原子跳跃是一个热激活过程；

⑷、计算晶粒长大概率：为了体现晶粒长大过程的不可逆性，只有达到最小能量条件的晶粒长大过程才是最稳定和最该发生的；为体现随机性，以最大晶粒长大速度为参照计算晶界原子的转变概率：

公式八： $P_i=\left\{\begin{array}{c}0,G_{T,i}<{\mathrm{\&Delta;G}}_A-{\mathrm{\&Sigma;\&Delta;G}}_i\\ \frac{M_i}{M_{\max }}\frac{F_i}{F_{\max }},G_{T,i}\&GreaterEqual;{\mathrm{\&Delta;G}}_A-{\mathrm{\&Sigma;\&Delta;G}}_i\end{array}\right.$

其中，M_i为原子i所在晶粒的晶界迁移率，F_i为原子i的晶粒长大驱动力，M_max为最大晶界迁移率，F_max为最大晶粒长大驱动力，当晶界原子可以向一个以上临近晶粒转变时，实际上只会向转变概率最大的晶粒迁移；

⑸、显示温度场分布以及焊接接头晶粒形态和大小拓扑形貌；

步骤三：根据温度场计算获得的冷却速度和β相晶粒分布计算结果，计算连续冷却固态相变：

⑴′、定义初始参数，包括溶质扩散系数、界面迁移率、冷却速度、过冷度、溶质原子在β和α相中的平衡溶解度，设置焊接热影响区晶粒长大计算结果为初始组织：

利用Thermo-Calc相图计算软件进行热力学计算，获得相变特征温度，即β和α相共存温度区间，以及溶质原子在该温度区间的溶解度，由于焊接过程快速加热和冷却，偏离热力学平衡状态，模拟β相向α相转变需考虑过冷度与特定温度、压力下的稳态平衡相比，溶质在两相中的实际溶解度按照相变平衡温度和相变过冷度计算，采用热模拟试验测量合金的CCT曲线，通过CCT曲线，获得不同冷却速度下连续冷却非平衡相变的起始温度和终止温度，将不同冷却速度的过冷度和相变温度区间、时间作为相变计算的输入条件，设置焊接热影响区晶粒长大计算结果为初始组织；

⑵′、计算相变的形核率，设置形核：

晶界形核稳态形核率I_稳为：

公式九：

其中，N_v为单位体积核胚的数量，k为波尔兹曼常数，h为普朗克常量，Q为溶质原子在β相中的热扩散激活能，△G^*为形核驱动力，

当过冷度较大或者冷却速度较大时，认为形核在相变初期瞬间完成且形核位置从晶界向晶内均匀形核转移；

⑶′、计算生长态元胞的α相长大相变驱动力、速度，计算α相前沿生长态元胞的α相转变分数、α相溶质浓度、β相溶质浓度及整个计算域的溶质场分布：

相界面移动速度V_P：

公式十：V_P＝M_PF_P/V_β

其中，M_P是相界面迁移率，F_P为相变驱动力，

当以理想稀溶液处理时，相变驱动力F_P为：

公式十一： $F_P\&ap;\frac{\mathrm{RT}}{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{(X_i-X_{\mathrm{Al}}^{\mathrm{\&beta;}})}^2$

其中，X_i是界面处β相中溶质原子铝的摩尔百分数，X_Al ^β为β/α相平衡时β相的溶质浓度；

对于界面或混合控制型相变，由相长大速度获得界面处母相中铝原子的摩尔百分数：

公式十二： $X_i=X_{\mathrm{Al}}^{\mathrm{\&beta;}}+{\left(\frac{{2V}_P\&CenterDot;X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}\right)}^{\frac{1}{2}}$

β相溶质扩散方程为：

公式十三： $\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;x}\left(D_{\mathrm{\&beta;}}\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(D_{\mathrm{\&beta;}}\frac{{\&PartialD;C}_{\mathrm{\&beta;}}}{\&PartialD;y}\right)$

α相溶质扩散方程为：

公式十四： $\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;x}\left(D_{\mathrm{\&alpha;}}\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;x}\right)+\frac{\&PartialD;}{\&PartialD;y}\left(D_{\mathrm{\&alpha;}}\frac{{\&PartialD;C}_{\mathrm{\&alpha;}}}{\&PartialD;y}\right)$

其中，D_β为β相溶质扩散系数，D_α为α相溶质扩散系数，C_β为β相溶质浓度，C_α为α相溶质浓度，表示x方向求偏导数、表示y方向求偏导数；

根据界面处溶质守恒，获得：

公式十五： $V_P(X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_i)=D(n\&CenterDot;\&dtri;C{|}_{\mathrm{interface}})$

其中，X_Al ^α为β/α相平衡态α相的溶质浓度，V_P(X_Al ^α-X_i)表示单位面积单位时间相界面处从β相脱溶析出的溶质原子量，表示界面处β相中溶质梯度引起的扩散流量，为界面处β相的溶质浓度梯度，n为β/α相界面的法向矢量；

因此，对于混合型相变，相界面移动速度V_P在x轴方向的差分表达式为：

公式十六：

$\begin{array}{c}{V}_{P,x}(i,j)=\frac{D_{\mathrm{\&beta;}}}{\mathrm{\&Delta;x}[X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M_P\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}{V}_P\right)}^{1/2}]}\left[\begin{array}{c}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M_P\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}{V}_{P,x}\right)}^{1/2}+C_{\mathrm{\&beta;}}(i-1,j)f_{\mathrm{\&beta;}}(i-1,j)\\ +X_{\mathrm{Al}}^{\mathrm{\&beta;}}-{\left(\frac{2X_{\mathrm{Al}}^{\mathrm{\&beta;}}}{M\&CenterDot;(\mathrm{RT}/V_{\mathrm{\&beta;}})}\right)}^{1/2}+C_{\mathrm{\&beta;}}(i+1,j)f_{\mathrm{\&beta;}}(i+1,j)\end{array}\right]\\ -\frac{D_{\mathrm{\&alpha;}}}{\mathrm{\&Delta;x}[X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}}]}\left[\begin{array}{c}(-X_{\mathrm{Al}}^{\mathrm{\&alpha;}}+C_{\mathrm{\&alpha;}}(i-1,j))f_{\mathrm{\&alpha;}}(i-1,j)\\ +(-X_{\mathrm{Al}}^{\mathrm{\&alpha;}}+C_{\mathrm{\&alpha;}}(i+1,j))f_{\mathrm{\&alpha;}}(i+1,j)\end{array}\right]\end{array}$

其中V_P,_x(i,j)为计算域中坐标为(i，j)的节点在x方向的α相长大速度，△x为空间步长，V_P,_x为上一时间步长α相尖端在x方向的长大速度，C_β(i-1,j)为坐标(i-1，j)的节点的β相溶质浓度，f_β(i-1,j)为坐标(i-1，j)的节点的β相体积分数，C_β(i+1,j)为坐标(i+1，j)的节点的β相溶质浓度，f_β(i+1,j)为坐标(i+1，j)的节点的β相体积分数，C_α(i-1,j)为坐标(i-1，j)的节点的α相溶质浓度，f_α(i-1,j)为坐标(i-1，j)的节点的α相体积分数，C_α(i+1,j)为坐标(i+1，j)的节点的α相溶质浓度，f_α(i+1,j)为坐标(i+1，j)的节点的α相体积分数；

⑷′、计算α相尖端曲率，修正热力学平衡β相溶质浓度：根据Thermo-Calc软件计算可知，当热力学温度一定时，α相的平衡溶解度X_Al ^α确定，此时，当α相具有无限大曲率半径时，溶质原子在β相中的热力学平衡溶解度为X_Al ^β；当α相尖端曲率为k时，还需考虑β/α相界面曲率效应，即溶质原子在β相中的热力学平衡溶解度应为即将X_Al ^β修正为 $X_{\mathrm{Al}}^{\mathrm{\&beta;}}\stackrel{\&OverBar;}{\left(k\right)},$ 对于稀溶液应满足 $1n\frac{X_{\mathrm{Al}}^{\mathrm{\&beta;}}\stackrel{\&OverBar;}{\left(k\right)}}{X_{\mathrm{Al}}^{\mathrm{\&beta;}}}=\mathrm{\&Gamma;}\stackrel{\&OverBar;}{k},$ 由此：

公式十七： $X_{\mathrm{Al}}^{\mathrm{\&beta;}}\stackrel{\&OverBar;}{\left(k\right)}=X_{\mathrm{Al}}^{\mathrm{\&beta;}}(1+\mathrm{\&Gamma;}\stackrel{\&OverBar;}{k})$

其中，为α相尖端曲率，Γ为Gibbs-Thomson系数，

α相尖端曲率：

公式十八： $\stackrel{\&OverBar;}{k}=\frac{1}{\mathrm{\&Delta;x}}(1-2\frac{f_{\mathrm{\&alpha;}}+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{f}_{\mathrm{\&alpha;}}\left(k\right)}{N+1})$

其中，f_α生长态元胞的α相分数，N为元胞自动机的邻居元胞数目，对于冯·诺依曼型邻居关系N＝4，对于摩尔型邻居关系型邻居关系N＝8；

Gibbs-Thomson系数：

公式十九： $\mathrm{\&Gamma;}=\frac{{2\mathrm{\&sigma;V}}_{\mathrm{\&beta;}}}{\mathrm{RT}(X_{\mathrm{Al}}^{\mathrm{\&alpha;}}-X_{\mathrm{Al}}^{\mathrm{\&beta;}})}$

其中，σ为β/α的相界界面能；

生长态元胞α相分数：

公式二十： ${\mathrm{\&delta;f}}_{\mathrm{\&alpha;}}=\frac{\mathrm{\&delta;t}}{\mathrm{\&Delta;x}}(V_x+V_y-V_x{V}_y\frac{\mathrm{\&delta;t}}{\mathrm{\&Delta;x}})$ 并且f_α ⁿ＝f_α ⁰+δf_α

其中，δt为计算时间步长，V_x为生长态元胞在x方向的相界面移动速度，V_y为生长态元胞在y方向的相界面移动速度，f_α ⁰为上一时间步长的生长态元胞的α相分数，f_αn分别为当前时间步长的生长态元胞的α相分数，δf_α为生长态元胞的α相分数增量；

⑸′、计算各向异性相界面迁移率：钛合金从β相区连续冷却时，α相与β相基体之间存在严格的伯格斯晶体学取向关系，以界面迁移率的各向异性表征，在择优生长方向和与其垂直方向的界面迁移率分别记为M₁和M₂，根据相界迁移的微观物理机制，获得界面迁移率的表达式：

公式二十一：M＝M₀exp(-Q_m/RT)

公式二十二： $M_0=\frac{1}{16}{(V_{\mathrm{\&alpha;}}^{\frac{1}{3}}+V_{\mathrm{\&beta;}}^{\frac{1}{3}})}^4{V}_D{\left(\mathrm{kT}\right)}^{-1}$

其中，M₀为指数前因子，V_α为α相的摩尔体积、ν_D≈5×10¹³/s为德拜频率；由于M₁/M₂的大小直接决定析出新相的长细比，取M₁/M₂＝10、100或1000，分别进行析出α相形态的模拟，并且通过扫描电镜观测相同焊接工艺条件下热影响区的组织形态，将析出α相形态的模拟结果与实验结果对比，两者最吻合时采用的M₁/M₂最佳；

⑹′、计算晶内不同择优取向的α相长大形貌：新相长大方向与母相晶粒的晶体位向有关，当母相晶体取向为优先取向时，按照与β相晶粒的晶体取向一致方向长大，当β相晶粒的晶体取向不是优先取向时，则与β相晶粒的晶体取向呈一特定角度长大，由于对于钛合金而言，β-α相变存在多个惯习面和惯习方向，所以每个β相晶粒内α相长大的位相取向并不唯一，因此，α相生长方向与β相的晶界呈一定角度，最终经历连续冷却非平衡相变后，在焊接热影响区形成网篮组织；

⑺′、更新元胞状态，显示组织形态及溶质场分布、不同位置不同时间的α相变分数，并以文件形式保存。