CN106202675A - 预测钛合金等温成形与动态再结晶演化耦合响应的方法 - Google Patents

Abstract

一种预测钛合金等温成形与动态再结晶演化耦合响应的方法，将得到的钛合金等温成形过程中晶粒的应力响应、晶粒间不均匀变形及其导致的不均匀位错密度作为变量传递给元胞自动机，得到该晶粒尺度不均匀变形条件下的动态再结晶演化，获得动态再结晶形核、长大的组织形态及其导致的晶界演化和更新的位错密度。将获得的动态再结晶形核及长大等晶粒与晶界信息及位错密度返回晶体塑性有限元方法，更新每个晶粒单元的位错滑移抗力，从而影响钛合金的后续变形，并通过本构关系计算晶粒尺度的应力响应。本发明实现了钛合金等温成形晶粒尺度不均匀变形、动态再结晶组织形态演变、再结晶晶粒尺寸演化、再结晶动力学、变形体和晶粒的流变应力的同步预测。

Description

预测钛合金等温成形与动态再结晶演化耦合响应的方法

技术领域

本发明涉及塑性加工领域，具体是一种预测钛合金等温成形与动态再结晶演化耦合响应的方法

背景技术

等温成形技术是实现钛合金复杂构件成形成性一体化精确制造的有效途径。然而，钛合金等温成形是一个多参数、多场耦合作用下的不均匀变形过程，强烈的不均匀变形使得钛合金等温成形过程中动态再结晶形核与长大规律十分复杂，而动态再结晶又显著地影响钛合金的力学响应，进而对后续变形产生显著影响，加剧变形的不均匀性。等温成形与动态再结晶演化的耦合作用使得钛合金变形与微观组织的预测进而实施控制面临巨大挑战。因此，发展钛合金等温成形与动态再结晶耦合响应的精确预测方法，对发展钛合金等温成形与组织演化全耦合预测理论与方法，进而实现钛合金构件等温成形和组织性能的有效控制的具有重要的理论意义和工业应用价值。

在ZL200910012127.4的发明中公布了采用元胞自动机法对板带钢再结晶组织演变进行预测的方法，获得了板带钢高温变形过程中的动态再结晶形核与长大及组织形态演化规律，并基于平均位错密度结合Kocks-Mecking(K-M)方程计算了材料的应力响应。但是，这些元胞自动机模型是以晶粒均匀变形为前提条件，未考虑实际成形过程中晶粒尺度的不均匀变形的作用；同时，基于平均位错密度计算的材料的应力响应不能反映材料实际成形过程中位错的不均匀分布对应力响应的影响；另外，没有考虑动态再结晶组织演化及材料力学响应变化对后续高温变形过程的影响。从而使得预测结果的可靠性差。针对这些存在的问题，西北工业大学的武川等(Wu Chuan,Yang He,Li Hong Wei.Modeling ofdiscontinuous dynamic recrystallization of a near-a titanium alloyIMI834during isothermal hot compression by combining a cellular automatonmodel with a crystal plasticity finite element method.Computational MaterialsScience.2013,79:944～959.)将晶体塑性有限元与元胞自动机通过网格映射和变形信息单向传递的方法结合起来，以描述晶粒尺度不均匀变形对动态再结晶形核与长大的影响，并应用于钛合金等温成形的不连续动态再结晶预测，获得了钛合金等温成形的晶粒变形形态和动态再结晶动力学规律。但该模型依然采用平均位错计算材料的应力响 应，且该模型中晶体塑性有限元与元胞自动机采用相互独立的网格模型和坐标系，且两者的网格密度不一致，这使得数据传递比较困难，且最主要的是仅实现了单向传递，导致材料变形与动态再结晶演化被割裂成两个系统，仍然没有考虑动态再结晶对材料后续变形的影响。加拿大学者Popova等(E.Popova,Y.Staraselski,A.Brahme,R.K.Mishra,K.Inal.Coupled crystalplasticity-probabilistic cellular automata approach to model dynamicrecrystallization in magnesium alloys.International Journal ofPlasticity.2014,66:85～102)在武川模型的基础上，通过采用相同网格模型，改善了武川模型中数据单向传递的精度，但没有从根本解决上述问题。

另外，上述工作均采用二维元胞自动机模型，难以描述实际材料的真实微观特征，如三岔晶界、四岔晶界等，显著的制约和限制了材料成形与动态再结晶规律的预测精度。

发明内容

为克服现有技术中存在的未考虑组织演化对后续变形影响的不足，本发明提出了一种预测钛合金等温成形与动态再结晶演化耦合响应的方法

本发明的具体实施步骤是：

步骤1，生成α+β两相钛合金等轴晶初始组织：

所述两相钛合金等轴晶初始组织中，α相体积分数为22％，晶粒总数为40的等轴晶初始组织。

步骤2，建立包含晶粒和晶界的钛合金等温成形三维有限元模型：

将步骤1生成的α+β两相钛合金等轴晶初始组织信息导入ABAQUS软件，生成包含α+β两相钛合金等轴晶组织的变形体。采用立方体单元对所述变形体划分网格，从而建立有限元几何模型。得到的限元几何模型为晶体塑性有限元方法和元胞自动机共用，即每个有限元单元也就是一个元胞。

对该有限元几何模型中的α相单元和β相单元分别赋予材料参数。该有限元几何模型采用周期性边界条件，加载方式为单轴压缩。从而建立了包含晶粒和晶界的钛合金等温成形三维有限元模型。

步骤3，读取相关数据：

读取各有限单元上一增量步的变形梯度、拉伸张量U^*和晶粒滑移系位错密度ρ^a、晶粒半径、晶粒取向Q、相信息参数Ph、晶界参数Bdy、格林应变的共轭应力T^*、滑移系a上滑移方向矢量m^a、滑移系a上滑移面单位矢量n^a以及当前增量步的变形梯度。所述晶粒滑移系位错密度ρ^a的初始位错密度取10⁹m^-2；所述滑移系滑移方向m^a、滑移面矢量n^a均位于晶体轴坐标系中。

若读取的相信息参数Ph为1则采用α相参数进行后续步骤的计算；若为读取的相信息参数Ph为0则采用β相参数进行后续步骤的计算。若读取的晶界参数Bdy为0表示非晶界元胞，若读取的晶界参数Bdy为1表示晶界元胞。

步骤4，确定晶粒间不均匀变形、单元的应力响应和晶粒取向演化：

所述不均匀变形指应变在晶粒间和晶粒内的不均匀分布。

确定晶粒间不均匀变形、单元的应力响应和晶粒取向演化的具体过程是：

第一步、通过分解剪应力τ^a和滑移系抗力s^a确定剪切应变速率

其中，γ₀为参考剪切应变速率，χ为率敏感因子。

第二步、计算变形梯度的塑性部分、非塑性部分和格林应变。

第三步、通过弹性本构关系计算单元的柯西应力。

第四步、计算晶粒取向的转动。

至此，确定了晶粒间的不均匀变形、应力响应和晶粒取向演化。

所述计算分解剪应力τ^a和滑移系抗力s^a的具体过程是：

第一步、根据公式确定滑移系分解剪应力τ^a，其中T^*为格林应变的共轭应力，其中：为Schmid张量，τ^a是滑移系分解剪应力。

第二步、根据位错密度计算滑移系抗力s^a

其中，s₀为初始滑移抗力，s_d为位错交互作用对滑移抗力的影响，s_HP为霍尔佩奇效应对滑移抗力的影响，λ为晶界软化因子。初始滑移抗力表示为变形温度的函数：

其中，c₃为与固溶强化相关的强度系数，c₆为与材料相关的参考温度，θ为单元温度。

位错交互作用对滑移抗力的影响s_d为：

其中，c₄为材料常数，其值在0.1～1之间；b是泊氏矢量；t-△t指上一增量步，相应的t则指的是当前增量步；μ为剪切弹性模量，与温度有关，通过公式(4)确定：

霍尔佩奇效应s_HP代表的晶粒尺寸对变形的阻碍作用表示为：

其中，c₅为霍尔佩奇系数，r为晶粒半径。

钛合金晶界在高温变形中的作用如图3所示，通过公式(6)计算软化因子：

其中，θ_melt为材料熔点；为等效塑性应变；为宏观应变速率。

步骤5，计算晶粒的位错密度演化；

根据晶粒的滑移剪切应变速率计算晶粒的位错密度演化。同时考虑动态再结晶与变形中的动态回复，通过公式(17)计算位错密度：

其中，k_i为由于再结晶软化产生的位错密度修正系数，ρ^a+为可动位错不动化造成的位错密度增加，即加工硬化效应。ρ^a-为位错攀移、异号位错相消造成的位错泯灭。ρ^a+和ρ^a-分别通过公式(18)和(19)计算：

其中，(18)和(19)式中林位错与不动位错密度存在如下关系：

其中，N为滑移系总数；m与n分别为位错滑移方向与滑移面法向方向；为滑移 系b的不动位错密度。

步骤6，模拟动态再结晶演化。

将晶粒位错密度作为变量传递给元胞自动机，计算动态再结晶演化。具体过程是：

第一步、确定元胞状态变量：

所述元胞状态变量包括：

位错密度状态变量：每个元胞的位错密度为以及滑移系位错密度为

相状态变量Ph：用于标定元胞的相信息，Ph＝1为α相，Ph＝0为β相；

再结晶状态变量Drx：用于标定元胞是否为新生成的再结晶相，Drx＝1表示该元胞为再结晶相，Drx＝0表示该元胞为母相；

再结晶潜在形核状态变量Potential：若母相元胞未转换为再结晶潜在形核元胞，则Potential＝0；否则，Potential＝1；

晶界状态变量Bdy：用于标定元胞是否在晶界处，若读取的晶界参数Bdy为0表示非晶界元胞，若读取的晶界参数Bdy为1表示晶界元胞。

晶粒取向状态变量Q：描述晶粒取向的张量；

晶粒尺寸状态变量r：该晶粒的平均晶粒半径；

晶粒编号状态变量Grain：包含母相晶粒与新生成的再结晶晶粒，每个晶粒都有不同的编号，处于同一个晶粒的元胞编号相同，处于不同晶粒的元胞编号不同。

第二步、根据形核规则判断是否形核；

当变形基体内位错密度与临界位错密度满足关系式ρ_i≥ρ_c时，新的无应变的晶核将以形核率在变形严重的区域形核，则该元胞转换为t时刻潜在形核元胞，即有：

其中，Q_int为界面能；l为位错自由程；Q_line为位错线能量；M为晶界迁移率，M为：

其中，δD_0b为特征晶界厚度与绝对0度晶界自扩散系数之积；K_B波尔兹曼常数；Q_diffu自扩散激活能；R为气体常数；b为柏氏矢量。通过公式(23)计算位错线能量Q_line：

形核率由公式(24)计算：

其中，C_N为拟合参数，m为再结晶形核率指数，Q_act为形核激活能，θ为温度。

为保证在变形严重区域优先形核，根据每个元胞的位错密度值，引入形核概率P_nuclei(i)来描述潜在形核元胞i转变为再结晶晶核的概率，P_nuclei(i)计算为：

其中，ρ_(i)为t时刻潜在形核元胞i的位错密度，ρ_max为t时刻潜在再结晶形核元胞中最大位错密度。随机选取N(t)个t时刻潜在再结晶形核元胞并生成随机数ξ，如果ξ≤P_nuclei(i)，则该元胞转换为再结晶晶核元胞，即有：

第三步、再结晶长大；

一旦元胞开始形核，就会以速率v向其临近元胞长大，考虑单独形核和集中形核两种情况，对每一个母相元胞，计算其被相邻再结晶元胞吞噬的概率：

若中心母相元胞i只有一个再结晶元胞邻居j，长大距离s_i,j通过公式(27)：

其中，t₁为上一次s_i,j置零的时间，晶界迁移的速度v为：

v_i,j＝MP (28)

其中，M为晶界迁移率，P为晶界迁移驱动力。晶界迁移驱动力P为：

P＝P_d-P_r (29)

其中，P_d为促进晶界迁移的驱动力，P_r为阻碍晶界迁移的驱动力。P_d为：

其中，ρ_M与ρ_R分别为母相与再结晶相的位错密度。P_r为：

r_eq为晶界等效曲率半径,Q_m大角度晶界能。等效曲率半径为：

其中，N_i为晶粒编号Grain＝i的元胞总数目，V_cell为单个元胞所代表的实际变形组织的体积。大角度晶界能计算为：

其中，O_m为晶粒取向差，v为泊松比，μ为剪切模量。对于元胞i与元胞j，晶粒取向差可以计算为：

其中，Q_i与Q_j分别为元胞i与元胞j的晶粒取向矢量矩阵。

引入球形膨胀概率因子P_i,j来描述中心元胞被再结晶邻居吞噬转换为再结晶元胞的概率，P_i,j计算为：

其中，r₀为参考半径，0.98、0.57与0.17分别表示为以中心元胞的为球心，以1.5r₀为半径的球体被邻居元胞切割的体积与单位元胞的体积的比值。

如果s_i,j≥l，则生成随机变量ζ∈[0,1]。若ζ≤P_i,j，则中心元胞i转换为再结晶元胞，长大距离状态变量s_i,j置零，中心母相元胞i的位错密度状态变量、相状态变量、再结晶状态变量、晶粒取向状态变量以及晶粒编号状态变量均与相邻再结晶元胞j的状态变量相等。

如果中心元胞有多个再结晶邻居，对每个再结晶邻居进行判断，若中心元胞满足能够同时被多个邻居吞噬的条件，则优先被具有较大P_i,j的再结晶元胞吞噬。具体规则如下所示：

若中心元胞具有n个再结晶元胞邻居，则记录该中心的所有再结晶元胞邻居{j₁,j₂,......,j_n}。

计算中心元胞i对所有再结晶元胞邻居{j₁,j₂,......,j_n}的生长概率{P_growth(1),P_growth(2),......,P_growth(n)}。

从再结晶元胞邻居{j₁,j₂,......,j_n}选取最大P_growth对应的邻居元胞j_x，计算长大距离s_i,j，且生成随机变量ζ∈[0,1]。若s_i,j>l且ζ≤P_i,j，则中心元胞被元胞j_x吞噬。若不满足条件，则再结晶元胞邻居集{j₁,j₂,......,j_n}去掉元胞j_x元素，生长概率{P_growth(1),P_growth(2),......,P_growth(n)}去掉P_growth(x)，然后对剩下的n-1个邻居元胞按上述方法继续进行计算。

如果s_i,j<l,则长大距离状态变量s_i,j随时间t继续积累。

第四步、晶界的演化；

计算每次形核和长大发生后的晶界演化。其方法为：

判断中心元胞与邻居元胞的晶粒编号状态变量Grain是否相等，即判断是否满足：

Grain_(j)＝Grain_(i)；j∈{元胞i的所有邻居} (36)

若满足上述条件，则该元胞处于晶粒内部，晶界状态变量Bdy＝0；若不满足则该元胞处于晶界上，即有晶界状态变量Bdy＝1。

至此，完成动态再结晶演化的模拟。

步骤7，更新晶粒的位错密度；

根据模拟得到的动态再结晶演化更新晶粒的位错密度。

步骤8，确定位错密度修正系数；

步骤9，确定再结晶体积分数；

式中，X_new为新生成的再结晶晶粒，即Drx＝1的元胞所占元胞数目，X为元胞空间总元胞数目。

步骤10，确定当前增量步变形体的流变应力σ；

通过公式(41)计算流变应力并输出：

其中T_i为每一个有限元单元的应力响应，NG为晶粒数目；

步骤11，存储公共变量；

将更新后的相信息Ph、晶界参数Bdy、位错密度ρ_i、晶粒尺寸r_eq、再结晶体积分数f作为公共变量存储，以便于下一增量步采用晶体塑性有限元法计算下一时刻的不均匀变形；

步骤12，当前时刻晶体塑性有限元与元胞自动机的一次耦合。输出动态再结晶的组织形态、再结晶晶粒尺寸、流变应力曲线、再结晶体积分数及其随变形的演化；至此，完成了当前时刻晶体塑性有限元与元胞自动机的一次耦合。

步骤13，返回有限元模型，重复步骤3～12，进行下一增量步的计算，直至完成钛合金等温成形与动态再结晶演化耦合响应的预测。

本发明的特征在于：采用晶粒长大的3D元胞自动机算法生成钛合金三维微观组织模型1，基于ABAQUS/Explicit软件平台建立该钛合金组织模型的有限元网格，同时作为晶体塑性有限元计算变形的网格模型和元胞自动机计算动态再结晶的元胞模型。然后，采用晶体塑性有限元方法计算钛合金等温成形过程中晶粒的应力响应、晶粒间不均匀变形及其导致的不均匀位错密度分布，并将晶粒的位错密度等作为变量传递给元胞自动机，计算在该晶粒尺度不均匀变形条件下的动态再结晶演化，获得动态再结晶形核、长大的组织形态及其导致的晶界演化和更新的位错密度，然后将元胞机 计算获得的动态再结晶形核及长大等晶粒与晶界信息及位错密度返回晶体塑性有限元方法，更新每个晶粒单元的位错滑移抗力，从而影响钛合金的后续变形，并通过本构关系计算晶粒尺度的应力响应。从而反映钛合金等温成形与动态再结晶演化的耦合作用。

本发明基于三维钛合金微观组织的有限元网格模型，将3D元胞自动机算法完全耦合进晶体塑性本构模型，统一采用材料变形的物理时间，从而实现钛合金等温成形与动态再结晶演化耦合响应的同步预测，为研究揭示钛合金大型复杂构件等温成形与组织演化规律，并进而实施有效控制提供依据。

故而，本发明能够获得钛合金等温成形过程中晶粒的不均匀变形分布、不均匀位错密度分布、不均匀动态再结晶形核和长大及其导致的组织形态演化、动态再结晶导致的应力软化和应力软化对后续不均匀变形的作用，从而实现钛合金等温成形与动态再结晶演化耦合响应的统一预测。

本发明通过结合元胞自动机与晶体塑性有限元法，考虑了钛合金材料在热变形过程中不均匀变形、力学响应与动态再结晶之间的强烈耦合作用，实现了钛合金等温成形晶粒尺度不均匀变形、动态再结晶组织形态演变、再结晶晶粒尺寸演化、再结晶动力学、变形体和晶粒的流变应力的同步预测。

基于包含两相晶粒与晶界的初始组织，建立钛合金等温成形三维有限元模型。其中，元胞自动机与晶体塑性有限元共用网格，即每个有限元为一个元胞。通过晶体塑性有限元模型计算晶粒尺度的不均匀变形与力学响应、以及不均匀变形导致的晶粒取向演化与位错密度演化。其中，位错密度与晶粒取向等变形信息被传递给元胞自动机用于计算动态再结晶演化。基于晶体塑性有限元计算的变形信息，通过采用三维26邻居类型以及半概率性元胞转换规则，建立元胞自动机模型以计算钛合金动态再结晶组织演化。随后，动态再结晶导致的晶粒与晶界形态演化、晶粒取向演化、位错密度演化、以及晶粒尺寸演化将被传递给晶体塑性有限元进行下一增量步的变形与力学响应计算，进而实现不均匀变形、力学响应与动态再结晶演化的紧密耦合。同时，由于元胞自动机采用有限元的变形时间，不仅使得组织演化时间具备物理意义，同时也保证了变形、应力响应与组织演化的高度同步性。元胞自动机与有限元的结合，使得该模型不仅可通过有限元分析软件动态输出位错密度、应力应变分布、流动应力等变形信息，分别如附图5、8、9、11所示，从附图5中不同应变量的位错密度分布可以看 出，在不同的变形条件下位错密度分布并不均匀；同时可以可视化输出晶粒尺寸、再结晶动力学等组织演化信息，如附图10、12所示，附图10中可以看出再结晶的形成使平均晶粒尺寸降低，晶粒得到了细化，附图12中不同温度不同应变速率下的再结晶动力学曲线，可以看出再结晶体积分数随着应变的增加而增加，呈现典型的“S”型曲线，这与实际情况相符；此外，附图7中等效塑性应变与再结晶平均半径对比图表明了变形与再结晶组织演化的相互影响。因此，该模型有效的实现了变形与组织演化的统一预测，这对实现钛合金构件等温成形和组织性能的有效控制的具有重要的理论意义和工业指导价值。

附图说明

图1为本发明耦合模型流程图；

图2为耦合模型示意图；

图3为钛合金高温变形时晶界的软化作用；

图4为三维26邻居类型；

图5为位错密度演化，其中，图5a应变为0.02，图5b应变为0.1，图5c应变为0.2，图5d应变为0.5；

图6为再结晶晶粒平均长大速度，其中，图6a应变速率为1s^-1的再结晶晶粒平均长大速度，图6b应变速率为0.01s^-1的再结晶晶粒平均长大速度；

图7为等效塑性应变与再结晶平均半径对比图，其中，图7a是温度在1323K应变速率为0.01s^-1下应变量为0.01的等效塑性应变分布，图7b是温度在1323K应变速率为0.01s^-1下应变量为0.01的平均晶粒半径；

图8为变形条件下的应力分布图，其中，图8a是温度为1323K、应变速率为0.01s^-1、应变为0.05时的MISES应力图8b是温度为1323K、应变速率为0.01s^-1、应变为0.05时的平均晶粒尺寸；

图9为变形条件下的应变分布图，其中，图9a是等效塑性应变分布图，图9b是平均晶粒尺寸分布图；

图10为平均晶粒尺寸变化；

图11为流变应力曲线，其中，图11a是温度为1273K时不同应变速率下的流动应力，图11b是温度为1323K时不同应变速率下的流动应力；

图12为再结晶体积分数变化曲线，其中，图12a是温度为1273K时不同应变速 率下的再结晶体积分数，图12b是温度为1273K时不同应变速率下的再结晶体积分数。

图13是本发明的流程图。

图中：

1.三维微观组织模型；2.晶体塑性模型；3.3D元胞自动机算法；4.第一邻居类型；5.第二邻居类型；6.第三邻居类型；7.元胞邻居；8.中心元胞。

具体实施方式

本实施例具体是一种TA15钛合金等温压缩与动态再结晶演化耦合响应的预测方法。本实施例的实施基于Abaqus软件平台，并开发了用户材料子程序。步骤3～13的均是通过用户材料子程序来实现的，所述步骤3～12为一个时间增量步的具体计算内容，所有时间增量步的计算完成后，则实现了TA15钛合金等温压缩与动态再结晶演化的耦合响应。本实施例所述的TA15钛合金相关拟合参数如表1所示。

表1拟合参数

c1(α)

c1(β)

c20(α)

c20(β)

c3(α)

c3(β)

c4

c5

c6

8.7×106

9.1×106

2100

1900

5.1×103

2.1×103

0.5

4.2×103

0.5

本实施例的具体实施步骤包括：

步骤1，生成α+β两相钛合金等轴晶初始组织：

首先通过元胞自动机生成规格为80×80×80、α相体积分数为22％，晶粒总数为40的等轴晶初始组织。用不同颜色区分α相与β相的不同晶粒和晶界。

步骤2，建立包含晶粒和晶界的钛合金等温成形三维有限元模型：

将步骤1生成的α+β两相钛合金等轴晶初始组织信息保存为dat格式文件，并将该文件导入ABAQUS软件，生成包含α+β两相钛合金等轴晶组织的变形体。采用立方体单元对所述变形体划分网格，从而建立有限元几何模型。网格尺寸为1.25×1.25×1.25μm³。下述的晶体塑性有限元方法和元胞自动机共用该有限元几何模型，即每个有限元单元也就是一个元胞。对该有限元几何模型中的α相单元和β相单元分别赋予材料参数，所述相关材料参数见表2。该有限元几何模型采用周期性边界条件，加载方式为单轴压缩。从而建立了包含晶粒和晶界的钛合金等温成形三维有限元模型。

表2相关材料参数

步骤3，基于ABAQUS/Explicit的用户材料子程序VUMAT接口，读取相关数据：

基于ABAQUS/Explicit的用户材料子程序VUMAT接口，读取各有限单元上一增量步的变形梯度、拉伸张量U^*和晶粒滑移系位错密度ρ^a、晶粒半径、晶粒取向Q、相信息参数Ph、晶界参数Bdy、格林应变的共轭应力T^*、滑移系a上滑移方向矢量m^a、滑移系a上滑移面单位矢量n^a以及当前增量步的变形梯度。所述晶粒滑移系位错密度ρ^a的初始位错密度取10⁹m^-2；所述滑移系滑移方向m^a、滑移面矢量n^a均位于晶体轴坐标系中。

若读取的相信息参数Ph为1则采用α相参数进行后续步骤的计算；若为读取的相信息参数Ph为0则采用β相参数进行后续步骤的计算。若读取的晶界参数Bdy为0表示非晶界元胞，若读取的晶界参数Bdy为1表示晶界元胞。

步骤4，确定晶粒间不均匀变形、单元的应力响应和晶粒取向演化：

为了描述晶粒间不均匀变形，每一晶粒都被分为若干个有限单元，同一晶粒内的有限单元被赋予相同取向，而表示不同晶粒的有限单元的取向不同；采用晶体塑性本构模型进行计算，由此晶粒间和晶粒内的不均匀变形得以实现。所述不均匀变形指应变在晶粒间和晶粒内的不均匀分布。在确定晶粒间不均匀变形、单元的应力响应和晶粒取向演化之前首先需要计算分解剪应力τ^a和滑移系抗力s^a。

步骤4的具体过程是：

第一步、根据公式确定滑移系分解剪应力τ^a，其中T^*为格林应变的共轭应力，其中：为Schmid张量，τ^a是滑移系分解剪应力。

第二步、根据位错密度计算滑移系抗力s^a

其中，s₀为初始滑移抗力，s_d为位错交互作用对滑移抗力的影响，s_HP为霍尔佩奇效应对滑移抗力的影响，λ为晶界软化因子。初始滑移抗力表示为变形温度的函数：

其中，c₃为与固溶强化相关的强度系数，c₆为与材料相关的参考温度，θ为单元温度。

位错交互作用对滑移抗力的影响s_d为：

其中，c₄为材料常数，其值在0.1～1之间；b是泊氏矢量；t-△t指上一增量步，相应的t则指的是当前增量步；μ为剪切弹性模量，与温度有关，通过公式(4)确定：

霍尔佩奇效应s_HP代表的晶粒尺寸对变形的阻碍作用表示为：

其中，c₅为霍尔佩奇系数，r为晶粒半径。

钛合金晶界在高温变形中的作用如图3所示，通过公式(6)计算软化因子：

其中，θ_melt为材料熔点；为等效塑性应变；为宏观应变速率。

第三步、计算剪切应变速率

其中，γ₀为参考剪切应变速率，χ为率敏感因子。

第四步、计算变形梯度的塑性部分、非塑性部分和格林应变：

总的变形梯度F表示为：

F＝F^*·F^p (8)

式中F^*表示晶格畸变和刚性转动所产生的变形梯度，即非塑性部分；F^p则表示晶体沿着滑移方向的均匀剪切所对应的变形梯度，即塑性部分。塑性应变梯度变化率表示为：

其中，L^p为速度梯度的塑性部分，由下式计算：

由(8)式有：

F^*＝F·(F^p)^-1 (11)

通过公式(12)得到格林应变：

其中，I为二阶identity张量。

第五步、通过弹性本构关系计算单元的柯西应力T：

其中，T^*为格林应变E^*的共轭应力，为四阶弹性张量。

根据下式计算柯西应力T：

第六步、计算晶粒取向的转动：

变形过程中，晶粒除了发生塑性剪切应变，还会产生刚性转动，使得晶粒取向演变。通过公式(15)确定晶粒取向演化：

Q＝R^*Q(R^*)^T, (15)

其中，Q是代表晶粒取向的三个欧拉角表示的旋转张量，R^*代表材料晶格的刚性转动，与非塑性变形梯度有关。所述R^*通过下式计算：

F^*＝R^*U^*, (16)

其中U^*为拉伸张量。

至此，确定了晶粒间的不均匀变形、应力响应和晶粒取向演化。

步骤5，计算晶粒的位错密度演化；

根据晶粒的滑移剪切应变速率计算晶粒的位错密度演化。同时考虑动态再结晶与变形中的动态回复，通过公式(17)计算位错密度：

其中，k_i为由于再结晶软化产生的位错密度修正系数，ρ^a+为可动位错不动化造成的位错密度增加，即加工硬化效应。ρ^a-为位错攀移、异号位错相消造成的位错泯灭。ρ^a+和ρ^a-分别通过公式(18)和(19)计算：

其中，(18)和(19)式中林位错与不动位错密度存在如下关系：

其中，N为滑移系总数；m与n分别为位错滑移方向与滑移面法向方向；为滑移系b的不动位错密度。

步骤6，模拟动态再结晶演化。

将晶粒位错密度作为变量传递给元胞自动机，计算动态再结晶演化。具体过程是：

第一步、确定元胞状态变量：

位错密度状态变量：每个元胞的位错密度为以及滑移系位错密度为

相状态变量Ph：用于标定元胞的相信息，Ph＝1为α相，Ph＝0为β相；

再结晶状态变量Drx：用于标定元胞是否为新生成的再结晶相，Drx＝1表示该元胞为再结晶相，Drx＝0表示该元胞为母相；

再结晶潜在形核状态变量Potential：若母相元胞未转换为再结晶潜在形核元胞，则Potential＝0；否则，Potential＝1；

晶界状态变量Bdy：用于标定元胞是否在晶界处，具体定义规则与步骤3中一致；

晶粒取向状态变量Q：描述晶粒取向的张量；

晶粒尺寸状态变量r：该晶粒的平均晶粒半径；

晶粒编号状态变量Grain：包含母相晶粒与新生成的再结晶晶粒，每个晶粒都有不同的编号，处于同一个晶粒的元胞编号相同，处于不同晶粒的元胞编号不同。

第二步、根据形核规则判断是否形核；

当变形基体内位错密度与临界位错密度满足关系式ρ_i≥ρ_c时，新的无应变的晶核将以形核率在变形严重的区域形核，则该元胞转换为t时刻潜在形核 元胞，即有：

其中，Q_int为界面能；l为位错自由程；Q_line为位错线能量；M为晶界迁移率，M为：

其中，δD_0b为特征晶界厚度与绝对0度晶界自扩散系数之积；K_B波尔兹曼常数；Q_diffu自扩散激活能；R为气体常数；b为柏氏矢量。通过公式(23)计算位错线能量Q_line：

形核率由公式(24)计算：

其中，C_N为拟合参数，m为再结晶形核率指数，Q_act为形核激活能，θ为温度。

为保证在变形严重区域优先形核，根据每个元胞的位错密度值，引入形核概率P_nuclei(i)来描述潜在形核元胞i转变为再结晶晶核的概率，P_nuclei(i)计算为：

其中，ρ_(i)为t时刻潜在形核元胞i的位错密度，ρ_max为t时刻潜在再结晶形核元胞中最大位错密度。随机选取N(t)个t时刻潜在再结晶形核元胞并生成随机数ξ，如果ξ≤P_nuclei(i)，则该元胞转换为再结晶晶核元胞，即有：

第三步、再结晶长大；

一旦元胞开始形核，就会以速率v向其临近元胞长大，考虑单独形核和集中形核两种情况，对与每一个母相元胞，计算其被相邻再结晶元胞吞噬的概率：

(i)若中心母相元胞i只有一个再结晶元胞邻居j，长大距离s_i,j通过公式(27)：

其中，t₁为上一次s_i,j置零的时间，晶界迁移的速度v为：

v_i,j＝MP (28)

其中，M为晶界迁移率，P为晶界迁移驱动力。晶界迁移驱动力P为：

P＝P_d-P_r (29)

其中，P_d为促进晶界迁移的驱动力，P_r为阻碍晶界迁移的驱动力。P_d为：

其中，ρ_M与ρ_R分别为母相与再结晶相的位错密度。P_r为：

r_eq为晶界等效曲率半径,Q_m大角度晶界能。等效曲率半径为：

其中，N_i为晶粒编号Grain＝i的元胞总数目，V_cell为单个元胞所代表的实际变形组织的体积。大角度晶界能可计算为：

其中，O_m为晶粒取向差，v为泊松比，μ为剪切模量。对于元胞i与元胞j，晶粒取向差可以计算为：

其中，Q_i与Q_j分别为元胞i与元胞j的晶粒取向矢量矩阵。

如图4所示，中心元胞易被第一邻居类型4吞噬转换为再结晶元胞，第二邻居类型5次之，最不容易被第三邻居类型6吞噬并转换为再结晶元胞。同时，再结晶元胞的等效曲率半径越大，中心母相元胞则越难被再结晶的元胞邻居7吞噬成为再结晶元胞。因此，引入球形膨胀概率因子P_i,j来描述中心元胞8被再结晶邻居吞噬转换为再结晶元胞的概率，P_i,j计算为：

其中，r₀为参考半径，0.98、0.57与0.17分别表示为以中心元胞的为球心，以1.5r₀为半径的球体被邻居元胞切割的体积与单位元胞的体积的比值。

如果s_i,j≥l，则生成随机变量ζ∈[0,1]。若ζ≤P_i,j，则中心元胞i转换为再结晶元胞，长大距离状态变量s_i,j置零，中心母相元胞i的位错密度状态变量、相状态变量、再结晶状态变量、晶粒取向状态变量以及晶粒编号状态变量均与相邻再结晶元胞j的状态变量相等。

如果s_i,j<l,则长大距离状态变量s_i,j随时间t继续积累。

(ii)中心元胞有多个再结晶邻居

依据(i)中的步骤对每个邻居都进行判断。若中心元胞满足能够同时被多个邻居吞噬的条件，则优先被具有较大P_i,j的再结晶元胞吞噬。具体规则如下所示：

(1)若中心元胞具有n个再结晶元胞邻居，则记录该中心的所有再结晶元胞邻居{j₁,j₂,......,j_n}。

(2)计算中心元胞i对所有再结晶元胞邻居{j₁,j₂,......,j_n}的生长概率{P_growth(1),P_growth(2),......,P_growth(n)}。

(3)从再结晶元胞邻居{j₁,j₂,......,j_n}选取最大P_growth对应的邻居元胞j_x，计算长大距离s_i,j，且生成随机变量ζ∈[0,1]。若s_i,j>l且ζ≤P_i,j，则中心元胞被元胞j_x吞噬。若不满足条件，则再结晶元胞邻居集{j₁,j₂,......,j_n}去掉元胞j_x元素，生长概率{P_growth(1),P_growth(2),......,P_growth(n)}去掉P_growth(x)。

(4)重复(3)步骤。

如果s_i,j<l,则长大距离状态变量s_i,j随时间t继续积累。

第四步、晶界的演化；

元胞的形核及长大都会影响晶界的演化，因此需要计算每次形核和长大发生后的晶界演化。其方法为：

判断中心元胞与邻居元胞的晶粒编号状态变量Grain是否相等，即判断是否满足：

Grain_(j)＝Grain_(i)；j∈{元胞i的所有邻居} (36)

若满足上述条件，则该元胞处于晶粒内部，晶界状态变量Bdy＝0；若不满足则该元胞处于晶界上，即有晶界状态变量Bdy＝1。

至此，完成动态再结晶演化的模拟。

步骤7，更新晶粒的位错密度；

根据模拟得到的动态再结晶演化更新晶粒的位错密度。

步骤8，确定位错密度修正系数；

步骤9，确定再结晶体积分数；

式中，X_new为新生成的再结晶晶粒，即Drx＝1的元胞所占元胞数目，X为元胞空间总元胞数目。

步骤10，确定当前增量步变形体的流变应力σ；

通过公式(41)计算流变应力并输出：

其中T_i为每一个有限元单元的应力响应，NG为晶粒数目；

步骤11，存储公共变量；

将更新后的相信息Ph、晶界参数Bdy、位错密度ρ_i、晶粒尺寸r_eq、再结晶体积分数f作为公共变量存储，以便于下一增量步采用晶体塑性有限元模型2计算下一时刻的不均匀变形；

步骤12，当前时刻晶体塑性有限元与元胞自动机的一次耦合。输出动态再结晶的组织形态、再结晶晶粒尺寸、流变应力曲线、再结晶体积分数及其随变形的演化；至此，完成了当前时刻晶体塑性有限元与元胞自动机的一次耦合。

步骤13，返回有限元模型，重复步骤3～12，进行下一增量步的计算，直至完成钛合金等温成形与动态再结晶演化耦合响应的预测。

在每次重复进行下一增量步计算时，读取各有限元单元前一时刻的变形梯度、当前时刻的变形梯度、拉伸张量U^*、晶粒滑移系位错密度ρ^a、晶粒半径、相信息参数Ph、晶界参数Bdy、格林应变的共轭应力T^*以及晶体轴坐标系中的滑移系滑移方向与滑移面矢量m^a和n^a。

通过本实施例，能够获得如下预测结果：

1、位错密度的演化。

图5是温度为1323K应变速率为0.01s^-1下不同应变量的位错密度分布。在材料变形的初期阶段，晶界的位错密度增殖更快，而基体晶粒内部的位错密度增殖速度较慢，这是因为：晶界的软化作用使得晶界的材料变形更加剧烈，同时晶界两侧的晶粒由于晶粒取向与变形抗力存在差异，晶界的协调作用使晶界的变形量更大，两者的共同作用导致晶界的位错密度增殖更快；同时，由于初期阶段晶界对变形协调的贡献较大，使得基体晶粒内部的变形较小，位错密度增殖速度较慢。随着变形量的增加，基体晶粒的位错密度逐渐增加，部分基体晶粒的位错密度达到临界位错密度发生动态再结晶，从而造成位错密度的降低。随着变形量的继续增加，基体晶粒的位错密度达到饱和值，再结晶晶粒由于变形造成位错密度逐渐增加且再结晶晶粒的位错密度与基体晶粒的位错密度差值逐渐减小。随着变形量进一步增加，再结晶晶粒的位错密度增加并趋近于基体晶粒的位错密度，整体的位错密度趋于稳定。

2、再结晶晶粒平均长大速度。

图6中(a)和(b)所示分别为温度在1273K应变速率为1s^-1与0.01s^-1下的再结晶晶粒平均长大速度。在变形的初期阶段，再结晶平均晶粒尺寸的平均增长速率为零，这是因为基体晶粒的位错密度较小，均未达到临界位错密度，不足以满足再结晶形核条件。随着变形量的增加，母相急剧大量形核，且母相位错密度与新生成的再结晶位错密度差值较大，为再结晶的长大提供了较大的驱动力，再结晶快速长大，再结晶平均晶粒尺寸的增长速率快速增加。随着变形量的增加，再结晶形核数目增长速度逐渐变慢，新生成的再结晶由于变形而导致位错密度逐渐增加，而基体晶粒位错密度逐渐 进入饱和而不发生变化，再结晶晶粒的位错密度与基体晶粒位错密度的差值逐渐降低，再结晶的驱动力逐渐减小，进而导致再结晶平均长大速率呈现降低的趋势。

3、等效塑性应变分布与平均晶粒尺寸分布对比图。

图7所示为温度在1323K应变速率为0.01s^-1下应变量为0.01的等效塑性应变分布(a)与平均晶粒半径对比图(b)。从图中可以看出，再结晶的形核分布并不均匀，再结晶形核位置主要集中在等效塑性应变较大的晶界处。形量大的区域并非一定形核，这主要有两点原因：首先，位错密度的演化速度，如位错增殖与湮灭主要与材料的瞬时变形速率相关，材料的瞬时变形速率越大，位错演化的速率越大，而材料的变形程度只反映材料变形随时间的累积程度；其次，再结晶形核是一个非稳态的热激活过程，其形核过程具有随机性。

4、应力分布图。

图8所示为温度在1323K应变速率为0.01s^-1下应变量为0.05的应力分布图与平均晶粒尺寸分布图。通过对比图8(a)与(b)可以发现，不连续动态再结晶对钛合金变形组织的应力分布有较大的影响，在动态再结晶晶粒的区域呈现出较低的应力，呈现“应力塌陷区”。

5、应变分布图。

图9为温度为1323K应变速率为0.1s^-1下真实应变为0.06的等效塑性应变分布与平均晶粒尺寸分布图。从图中可以看出，不连续动态再结晶对材料的等效塑性应变分布的影响较为明显。在不连续动态再结晶发生的区域，等效塑性应变相对较大，而基体晶粒内部的等效塑性应变相对较小。然而，在动态再结晶的大量聚集形核与长大的区域出现了变形量低于母相的“应变抑制区”。

6、平均晶粒尺寸。

图10为温度在1323k时不同应变速率下平均晶粒尺寸。从图中可以看出，平均晶粒尺寸随着变形量的增加而降低并达到稳定值，晶粒得到细化。应变速率越高，平均晶粒尺寸越小，这是因为：在高应变速率下，材料的变形较剧烈促进再结晶大量形核，然而再结晶晶核没有充足的时间长大，进而形成较小的再结晶晶粒。

7、流动应力预测。

图11(a)和(b)分别为温度在1273k和1323k时、不同应变速率下的流动应力预测，符号点代表实验数据，曲线为模拟值。在变形初始阶段，流动应力迅速增加至 峰值应力。随着变形的增加，流动应力逐渐呈现软化，这是因为变形导致了材料内部温度的升高以及再结晶的软化作用。从图中可以看出，在同一个温度下，应变速率越高，流动应力越高；同一个应变速率，变形温度越高，流动应力越低。

8、再结晶体积分数预测。

图12中(a)和(b)分别为TA15钛合金在温度1273K和1323K时、不同应变速率下的再结晶动力学曲线，图中符号点代表实验数据，曲线为模拟值。从图中可以看出，再结晶体积分数随着应变的增加而增加，呈现“S”型曲线。通过对比同一温度同一应变量下不同应变速率的再结晶体积分数可以看出，应变速率对再结晶体积分数有较大的影响，应变速率越高，再结晶体积分数越。而在相同的应变速率下，温度越高，再结晶体积分数越高。

Claims (3)

1.一种预测钛合金等温成形与动态再结晶演化耦合响应的方法，其特征在于，

具体实施步骤包括：

步骤1，生成α+β两相钛合金等轴晶初始组织：

所述两相钛合金等轴晶初始组织中，α相体积分数为22％，晶粒总数为40的等轴晶初始组织；

步骤2，建立包含晶粒和晶界的钛合金等温成形三维有限元模型：

将步骤1生成的α+β两相钛合金等轴晶初始组织信息导入ABAQUS软件，生成包含α+β两相钛合金等轴晶组织的变形体；采用立方体单元对所述变形体划分网格，从而建立有限元几何模型；得到的限元几何模型为晶体塑性有限元方法和元胞自动机共用，即每个有限元单元也就是一个元胞；

对该有限元几何模型中的α相单元和β相单元分别赋予材料参数；该有限元几何模型采用周期性边界条件，加载方式为单轴压缩；从而建立了包含晶粒和晶界的钛合金等温成形三维有限元模型；

步骤3，读取相关数据：

若读取的相信息参数Ph为1则采用α相参数进行后续步骤的计算；若为读取的相信息参数Ph为0则采用β相参数进行后续步骤的计算；若读取的晶界参数Bdy为0表示非晶界元胞，若读取的晶界参数Bdy为1表示晶界元胞；

步骤4，确定晶粒间不均匀变形、单元的应力响应和晶粒取向演化：

所述不均匀变形指应变在晶粒间和晶粒内的不均匀分布；

确定晶粒间不均匀变形、单元的应力响应和晶粒取向演化的具体过程是：

第一步、通过分解剪应力τ^a和滑移系抗力s^a确定剪切应变速率

${\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}^a={\mathrm{\&gamma;}}_0{\left(\frac{{\mathrm{\&tau;}}^a}{s^a}\right)}^{1/\mathrm{\&chi;}}\mathrm{sgn}\left({\mathrm{\&tau;}}^a\right)---\left(7\right)$

其中，γ₀为参考剪切应变速率，χ为率敏感因子；

第二步、计算变形梯度的塑性部分、非塑性部分和格林应变：

第三步、通过弹性本构关系计算单元的柯西应力：

第四步、计算晶粒取向的转动：

至此，确定了晶粒间的不均匀变形、应力响应和晶粒取向演化；

步骤5，计算晶粒的位错密度演化；

根据晶粒的滑移剪切应变速率计算晶粒的位错密度演化；同时考虑动态再结晶与变形中的动态回复，通过公式(17)计算位错密度：

${\mathrm{\&rho;}}_i^a\left(t\right)=k_i\&lsqb;{\mathrm{\&rho;}}_i^a(t-\mathrm{\&Delta;}t)+({\mathrm{\&rho;}}^{a+}+{\mathrm{\&rho;}}^{a-})\mathrm{\&Delta;}t\&rsqb;---\left(17\right)$

其中，k_i为由于再结晶软化产生的位错密度修正系数，ρ^a+为可动位错不动化造成的位错密度增加，即加工硬化效应；ρ^a-为位错攀移、异号位错相消造成的位错泯灭；ρ^a+和ρ^a-分别通过公式(18)和(19)计算：

${\mathrm{\&rho;}}^{a+}=c_1\sqrt{{\mathrm{\&rho;}}_F^a}{\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}^a---\left(18\right)$

${\mathrm{\&rho;}}^{a-}=-c_2{\mathrm{\&rho;}}_F^a{\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}^a---\left(19\right)$

其中，(18)和(19)式中林位错与不动位错密度存在如下关系：

${\mathrm{\&rho;}}_F^a=\underset{b=1}{\overset{N}{\&Sigma;}}{\mathrm{\&rho;}}_I^b\left|cos(n^a,m^b)\right|---\left(20\right)$

其中，N为滑移系总数；m与n分别为位错滑移方向与滑移面法向方向；为滑移系b的不动位错密度；

步骤6，模拟动态再结晶演化；

将晶粒位错密度作为变量传递给元胞自动机，计算动态再结晶演化；具体过程是：第一步、确定元胞状态变量：

所述元胞状态变量包括：

位错密度状态变量：每个元胞的位错密度为以及滑移系位错密度为相状态变量Ph：用于标定元胞的相信息，Ph＝1为α相，Ph＝0为β相；

再结晶状态变量Drx：用于标定元胞是否为新生成的再结晶相，Drx＝1表示该元胞为再结晶相，Drx＝0表示该元胞为母相；

再结晶潜在形核状态变量Potential：若母相元胞未转换为再结晶潜在形核元胞，则Potential＝0；否则，Potential＝1；

晶界状态变量Bdy：用于标定元胞是否在晶界处，若读取的晶界参数Bdy为0表示非晶界元胞，若读取的晶界参数Bdy为1表示晶界元胞；

晶粒取向状态变量Q：描述晶粒取向的张量；

晶粒尺寸状态变量r：该晶粒的平均晶粒半径；

晶粒编号状态变量Grain：包含母相晶粒与新生成的再结晶晶粒，每个晶粒都有不同的编号，处于同一个晶粒的元胞编号相同，处于不同晶粒的元胞编号不同；

第二步、根据形核规则判断是否形核；

当变形基体内位错密度与临界位错密度满足关系式ρ_i≥ρ_c时，新的无应变的晶核将以形核率在变形严重的区域形核，则该元胞转换为t时刻潜在形核元胞，即有：

$\left\{\begin{array}{cc}{\mathrm{\&rho;}}_i\&GreaterEqual;{\mathrm{\&rho;}}_c;& Potential\left(i\right)=1\\ {\mathrm{\&rho;}}_i<{\mathrm{\&rho;}}_c;& Potential\left(i\right)=0\end{array}\right.---\left(21\right)$

其中，Q_int为界面能；l为位错自由程；Q_line为位错线能量；M为晶界迁移率，M为：

$M=\frac{{\mathrm{b\&delta;D}}_{0b}}{K_B\mathrm{\&theta;}}\exp (-\frac{Q_{diffu}}{R\mathrm{\&theta;}})---\left(22\right)$

其中，δD_0b为特征晶界厚度与绝对0度晶界自扩散系数之积；K_B波尔兹曼常数；Q_diffu自扩散激活能；R为气体常数；b为柏氏矢量；通过公式(23)计算位错线能量Q_line：

$Q_{line}=\frac{1}{2}{\mathrm{\&mu;b}}^2---\left(23\right)$

形核率由公式(24)计算：

$\stackrel{\&CenterDot;}{n}=C_N{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}^m\exp (-\frac{Q_{act}}{R\mathrm{\&theta;}})---\left(24\right)$

其中，C_N为拟合参数，m为再结晶形核率指数，Q_act为形核激活能，θ为温度；为保证在变形严重区域优先形核，根据每个元胞的位错密度值，引入形核概率P_nuclei(i)来描述潜在形核元胞i转变为再结晶晶核的概率，P_nuclei(i)计算为：

$P_{nuclei\left(i\right)}=\frac{{\mathrm{\&rho;}}_{\left(i\right)}}{{\mathrm{\&rho;}}_{\max }}---\left(25\right)$

其中，ρ_(i)为t时刻潜在形核元胞i的位错密度，ρ_max为t时刻潜在再结晶形核元胞中最大位错密度；随机选取N(t)个t时刻潜在再结晶形核元胞并生成随机数ξ，如果ξ≤P_nuclei(i)，则该元胞转换为再结晶晶核元胞，即有：

$\left\{\begin{array}{cc}\mathrm{\&xi;}\&le;P_{nuclei\left(i\right)};& DRXhappens\\ \mathrm{\&xi;}>P_{nuclei\left(i\right)};& NoDRX\end{array}\right.---\left(26\right)$

第三步、再结晶长大；

一旦元胞开始形核，就会以速率v向其临近元胞长大，考虑单独形核和集中形核两种情况，对与每一个母相元胞，计算其被相邻再结晶元胞吞噬的概率：

若中心母相元胞i只有一个再结晶元胞邻居j，长大距离s_i,j通过公式(27)：

$s_{i,j}=\underset{t_1}{\overset{t}{\&Integral;}}{v}_{i,j}dt---\left(27\right)$

其中，t₁为上一次s_i,j置零的时间，晶界迁移的速度v为：

v_i,j＝MP (28)

其中，M为晶界迁移率，P为晶界迁移驱动力；晶界迁移驱动力P为：

P＝P_d-P_r (29)

其中，P_d为促进晶界迁移的驱动力，P_r为阻碍晶界迁移的驱动力；P_d为：

$P_d=\frac{1}{2}{\mathrm{\&mu;b}}^2({\mathrm{\&rho;}}_M-{\mathrm{\&rho;}}_R)---\left(30\right)$

其中，ρ_M与ρ_R分别为母相与再结晶相的位错密度；P_r为：

$P_r=\frac{2Q_m}{r_{eq}}---\left(31\right)$

r_eq为晶界等效曲率半径,Q_m大角度晶界能；等效曲率半径为：

$r_{eq}^{\left(i\right)}=\sqrt[3]{\frac{3N_i{V}_{cell}}{4\mathrm{\&pi;}}}---\left(32\right)$

其中，N_i为晶粒编号Grain＝i的元胞总数目，V_cell为单个元胞所代表的实际变形组织的体积；大角度晶界能计算为：

$Q_m=\frac{{\mathrm{\&mu;bO}}_m}{4\mathrm{\&pi;}(1-v)}---\left(33\right)$

其中，O_m为晶粒取向差，v为泊松比，μ为剪切模量；对于元胞i与元胞j，晶粒取向差可以计算为：

$O_m={\cos }^{-1}\&lsqb;\frac{trace\left(Q_i{Q}_j\right)-1}{2}\&rsqb;---\left(34\right)$

其中，Q_i与Q_j分别为元胞i与元胞j的晶粒取向矢量矩阵；

引入球形膨胀概率因子P_i,j来描述中心元胞被再结晶邻居吞噬转换为再结晶元胞的概率，P_i,j计算为：

$P_{i,j}=\left\{\begin{array}{cc}0.97;& j\&Element;\left\{\begin{array}{ccc}neighbor& type& 1\end{array}\right\}\\ 0.46;& j\&Element;\left\{\begin{array}{ccc}neighbor& type& 2\end{array}\right\}\\ 0.17;& j\&Element;\left\{\begin{array}{ccc}neighbor& type& 3\end{array}\right\}\end{array}\right.---\left(35\right)$

其中，r₀为参考半径，0.98、0.57与0.17分别表示为以中心元胞的为球心，以1.5r₀为半径的球体被邻居元胞切割的体积与单位元胞的体积的比值；

如果s_i,j≥l，则生成随机变量ζ∈[0,1]；若ζ≤P_i,j，则中心元胞i转换为再结晶元胞，长大距离状态变量s_i,j置零，中心母相元胞i的位错密度状态变量、相状态变量、再结晶状态变量、晶粒取向状态变量以及晶粒编号状态变量均与相邻再结晶元胞j的状态变量相等；

如果s_i,j<l,则长大距离状态变量s_i,j随时间t继续积累；

中心元胞有多个再结晶邻居：

对每个邻居都进行判断；若中心元胞满足能够同时被多个邻居吞噬的条件，则优先被具有较大P_i,j的再结晶元胞吞噬；具体规则如下所示：

若中心元胞具有n个再结晶元胞邻居，则记录该中心的所有再结晶元胞邻居{j₁,j₂,......,j_n}；

计算中心元胞i对所有再结晶元胞邻居{j₁,j₂,......,j_n}的生长概率{P_growth(1),P_growth(2),......,P_growth(n)}；

从再结晶元胞邻居{j₁,j₂,......,j_n}选取最大P_growth对应的邻居元胞j_x，计算长大距离s_i,j，且生成随机变量ζ∈[0,1]；若s_i,j>l且ζ≤P_i,j，则中心元胞被元胞j_x吞噬；若不满足条件，则再结晶元胞邻居集{j₁,j₂,......,j_n}去掉元胞j_x元素，生长概率{P_growth(1),P_growth(2),......,P_growth(n)}去掉P_growth(x)，按上述方法继续进行计算；

如果s_i,j<l,则长大距离状态变量s_i,j随时间t继续积累；

第四步、晶界的演化；

计算每次形核和长大发生后的晶界演化；其方法为：

判断中心元胞与邻居元胞的晶粒编号状态变量Grain是否相等，即判断是否满足：Grain_(j)＝Grain_(i)；j∈{元胞i的所有邻居} (36)

若满足上述条件，则该元胞处于晶粒内部，晶界状态变量Bdy＝0；若不满足则该元胞处于晶界上，即有晶界状态变量Bdy＝1；

至此，完成动态再结晶演化的模拟；

步骤7，更新晶粒的位错密度；

根据模拟得到的动态再结晶演化更新晶粒的位错密度；

步骤8，确定位错密度修正系数；

${\mathrm{\&rho;}}_i=\left\{\begin{array}{c}{\mathrm{\&rho;}}_0\\ {\mathrm{\&rho;}}_1(t-\mathrm{\&Delta;}t)+\underset{1}{\overset{24}{\&Sigma;}}{\mathrm{\&rho;}}_i^a\end{array}\right.---\left(38\right)$

$k_i=\frac{{\mathrm{\&rho;}}_i}{{\mathrm{\&rho;}}_i(t-\mathrm{\&Delta;}t)}---\left(39\right)$

步骤9，确定再结晶体积分数；

$f=\frac{X_{new}}{X}---\left(40\right)$

式中，X_new为新生成的再结晶晶粒，即Drx＝1的元胞所占元胞数目，X为元胞空间总元胞数目；

步骤10，确定当前增量步变形体的流变应力σ；

通过公式(41)计算流变应力并输出：

$\mathrm{\&sigma;}=\frac{1}{NG}\underset{i=1}{\overset{NG}{\&Sigma;}}{T}_i---\left(41\right)$

其中T_i为每一个有限元单元的应力响应，NG为晶粒数目；

步骤11，存储公共变量；

将更新后的相信息Ph、晶界参数Bdy、位错密度ρ_i、晶粒尺寸r_eq、再结晶体积分数f作为公共变量存储，以便于下一增量步采用晶体塑性有限元法计算下一时刻的不均匀变形；

步骤12，当前时刻晶体塑性有限元与元胞自动机的耦合；输出动态再结晶的组织形态、再结晶晶粒尺寸、流变应力曲线、再结晶体积分数及其随变形的演化；至此，完成了当前时刻晶体塑性有限元与元胞自动机的耦合；

步骤13，返回有限元模型，重复步骤3～12，进行下一增量步的计算，直至完成钛合金等温成形与动态再结晶演化耦合响应的预测。

2.如权利要求1所述预测钛合金等温成形与动态再结晶演化耦合响应的方法，其特征在于，

所述计算分解剪应力τ^a和滑移系抗力s^a的具体过程是：

第一步、根据公式确定滑移系分解剪应力τ^a，其中T^*为格林应变的共轭应力，其中：为Schmid张量，τ^a是滑移系分解剪应力；

第二步、根据位错密度计算滑移系抗力s^a

其中，s₀为初始滑移抗力，s_d为位错交互作用对滑移抗力的影响，s_HP为霍尔佩奇效应对滑移抗力的影响，λ为晶界软化因子；初始滑移抗力表示为变形温度的函数：

$s_0=c_3\exp (-\frac{\mathrm{\&theta;}}{c_6})---\left(2\right)$

其中，c₃为与固溶强化相关的强度系数，c₆为与材料相关的参考温度，θ为单元温度；

位错交互作用对滑移抗力的影响s_d为：

$s_d=c_4\mathrm{\&mu;}b\sqrt{{\mathrm{\&rho;}}_{\left(i\right)}^a(t-\mathrm{\&Delta;}t)}---\left(3\right)$

其中，c₄为材料常数，其值在0.1～1之间；b是泊氏矢量；t-Δt指上一增量步，相应的t则指的是当前增量步；μ为剪切弹性模量，与温度有关，通过公式(4)确定：

$\mathrm{\&mu;}=49.02-\frac{5.821}{\exp (181/\mathrm{\&theta;})-1}GPa---\left(4\right)$

霍尔佩奇效应s_HP代表的晶粒尺寸对变形的阻碍作用表示为：

$s_{HP}=\frac{c_5}{\sqrt{2r}}---\left(5\right)$

其中，c₅为霍尔佩奇系数，r为晶粒半径；

钛合金晶界在高温变形中的作用如图3所示，通过公式(6)计算软化因子：

$\mathrm{\&lambda;}=1-0.087ln\left(\frac{10.0}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}\right)\exp \left(\frac{\mathrm{\&theta;}}{{\mathrm{\&theta;}}_{melt}}\right)\exp (-{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}^p)---\left(6\right)$

其中，θ_melt为材料熔点；为等效塑性应变；为宏观应变速率。

3.如权利要求1所述预测钛合金等温成形与动态再结晶演化耦合响应的方法，其特征在于，读取各有限单元上一增量步的变形梯度、拉伸张量U^*、格林应变的共轭应力T^*和存储的晶粒滑移系位错密度ρ^a、晶粒半径、晶粒取向Q、相信息参数Ph、晶界参数Bdy、、滑移系a上滑移方向矢量m^a、滑移系a上滑移面单位矢量n^a以及当前增量步的变形梯度；所述晶粒滑移系位错密度ρ^a的初始位错密度取10⁹m-²；所述滑移系滑移方向m^a、滑移面矢量n^a均位于晶体轴坐标系中。