杭州铝材价格联盟

推荐|2018年1期优秀论文选登——基于晶体塑性的难变形材料不均匀变形多尺度建模研究进展

塑性工程学报2018-11-07 07:23:02

作者:詹梅,李宏伟,孙新新,黄栋,武川,王旭,陈姗姗

单位:西北工业大学 材料科学与工程学院


引文格式:詹  梅,李宏伟,孙新新,等.基于晶体塑性的难变形材料不均匀变形多尺度建模研究进展[J].塑性工程学报,2018,251):1-14. ZHAN Mei, LI Hongwei, SUN Xinxinet al. Research progress of the multi-scale modeling of the heterogeneous deformation of hard-to-deform material based on crystal plasticity[J].Journal of Plasticity Engineering2018,251):1-14.


摘要:


分析并综述了不均匀变形预测中的晶体塑性模型与求解算法,论述了考虑多相材料非均质特征与晶界异质性的晶体塑性建模难点与方法;分析并综述了用于微观组织形态演化预测的元胞自动机建模方法,并论述了考虑不均匀变形与组织形态演化交互作用机制的晶体塑性与元胞自动机耦合建模的关键难点及其解决方法。进一步论述了基于宏观有限元、细观晶体塑性有限元与微观元胞自动机的多尺度全耦合建模方法及其在钛合金复杂隔框构件等温局部加载成形多尺度耦合响应规律预测中的成功应用。最后分析了多尺度建模仍存在的问题,并展望了多尺度建模的发展趋势。



引言


我国航空、航天等领域的高端装备正朝着大运力、高机动、长寿命、低能耗的方向发展,这就对其关键构件提出了高性能、轻量化和高可靠的要求。采用轻质高强材料和整体复杂结构,进行精确塑性成形成性一体化制造,是实现该要求的重要途径。同时,也是当前我国航空航天高端装备发展的迫切需求和国际成形领域研究的热点与难点。然而,以钛合金、高强铝合金为代表的轻质高强材料,室温塑性差,常温难变形,高温组织难控制;整体复杂结构中存在大小尺寸极端结合,成形难度极大,成形缺陷难控制。为此,西北工业大学杨合等将局部加载方法与等温成形相融合,发展了等温局部加载成形工艺,解决了钛合金大型整体复杂隔框类构件等的精确成形与组织性能调控难题,实现了该类构件的高性能柔性省力成形制造,为难变形材料复杂构件成形成性一体化制造探索出了一条道路[1]。

然而,一方面,钛合金等难变形材料往往具有各向异性显著的密排六方(HCP)结构,且由性能迥异的两相或多相组成,在受载条件下会产生显著的不均匀变形。这种不均匀变形是由位错滑移、孪生、晶界变形、晶粒旋转等细观变形机制决定的。另一方面,复杂加载条件,尤其是多道次局部加载,涉及多次热力耦合加载,坯料存在多变形区,且各变形区的加载条件随时间变化。这就加剧了材料的不均匀变形,并导致复杂的组织演化。还可导致充填不满、折叠、穿筋、流线紊乱等宏观成形缺陷。因此,复杂构件的精确成形与组织演化的协同调控面临巨大挑战,而采用数值模拟精确预测难变形材料复杂构件成形与组织演化的规律是解决该挑战的关键前提。

其中,有限元建模仿真是预测成形及缺陷产生规律的有效方法。细观变形是由位错的开动、滑移、攀移等实现的,且与加载方向和晶粒取向密切相关。在变形的同时,晶格的弹性扭曲及相邻晶粒的变形协调作用,使得晶粒取向发生变化,从而产生变形织构。基于连续介质力学和位错滑移理论构建的晶粒塑性理论及晶体塑性有限元方法是预测该细观不均匀变形与织构演化的有效方法。在热力加载条件下,除了织构的演化还会产生微观组织的演化,例如晶粒粗化、动态再结晶、相变等,从而改变材料的后续变形性能。元胞自动机是通过对时间和空间的离散而预测组织形态演化的直接方法,与唯象的统计学公式相比,其预测结果更为直观,便于与试验结果进行对比。目前,在上述单个尺度上的建模仿真工作已经开展的比较全面。然而,宏观、细观、微观的变形与组织演化是同步发生又相互作用的,呈现出复杂的多尺度耦合响应机理与规律。这使得单个尺度上的建模仿真难以胜任多尺度耦合响应预测的需求。因此,迫切需要进行多尺度耦合建模,以实现对难变形材料复杂构件成形多尺度耦合响应规律的可靠预测。

本文通过分析晶体塑性、元胞自动机的建模难点,综述了难变形材料不均匀变形预测的晶体塑性建模、组织形态演化预测的元胞自动机建模及其耦合建模的方法。进一步论述了进行复杂构件成形预测的多尺度建模思路及其在钛合金复杂隔框构件等温局部加载成形中的应用。最后分析了目前多尺度建模中存在的问题,展望了多尺度建模的发展趋势。


1 不均匀变形预测的晶体塑性建模


1.1晶体塑性模型及求解算法


晶体塑性理论是基于晶体材料位错滑移机制和细观介质力学理论建立起来一种本构理论。该理论认为晶体材料的塑性变形主要通过位错沿特定滑移系的滑移来实现,并引入塑性剪切应变来描述滑移系上的位错运动,如式(1)所示。该理论运用统计学思想将大量不连续的位错运动看作连续塑性变形过程,从而与连续介质力学和运动学联系在一起[2]。因此,具有广泛的普适性。

式中:α为滑移系编号;γα为滑移系α的剪切应变率;γα0为参考剪切应变率;τα为滑移系剪切力;为滑移系抗力;m为率敏感系数。

式(1)也可以称为晶体塑性理论的控制方程,其中的1/m次方给该方程组的求解带来了困难。因为,对于金属材料而言,常温变形的率敏感系数m往往在0.01~0.05之间,而高温下m往往在0.1~0.3之间。这使得式(1)成为高阶非线性方程组,尤其在常温下方程的阶数更高,求解十分困难,导致模型求解的数值稳定性差。

对于该难题,学者们从20世纪二三十年代晶体塑性理论被提出开始至今,一直从事着求解算法的研究工作。其中具有代表性是Euler Forward显式算法[3]与麻省理工的Anand L等[4]提出的全隐式模型及Newton-Raphson(N-R)迭代求解方法。

基于Euler Forward显式算法的晶体塑性模型通常可以表示为:



式中:t为时间;Δt为增量步长; F为变形梯度;Fp为塑性变形梯度;F*为非塑性变形梯度;E*为格林应变;C为4阶弹性张量;T*为Kirchhoff应力;0为滑移系Schmidt张量;I为2阶单位张量。由此可以看出,Euler Forward显式算法不存在迭代,是根据上一增量步的解来获得当前增量步解的算法,因此不存在收敛性问题[5-7]。然而,必须有效控制显式算法的增量步长Δt才能满足求解稳定性要求[8],如式(3)所示:

式中:Le为网格特征单元尺寸,其主要由最小的有限元网格尺寸确定;λμ 为材料弹性常数;Θ为材料密度。显式晶体塑性模型主要用于解决高速动态问题。在显式晶体塑性模型的计算中,网格应尽可能形状规则且尺寸均匀,避免局部相对较小尺寸的网格降低稳定增量步长时间,降低求解效率[9]。

相比于显式算法,隐式算法可使用较大的增量步长完成求解。然而,隐式算法往往存在大量的迭代。虽然N-R迭代具有平方收敛速度,但其对迭代初值要求苛刻,距离真实解较远的初值往往导致迭代不收敛,造成计算失败。这个问题使得求解条件被限制:(1)晶体塑性模型往往是在给定变形梯度的条件下进行计算的,难以用于复杂的成形过程计算,尤其是与有限元结合的复杂过程计算;(2)要求增量步之间的变形及应力响应的变化很小,以满足迭代收敛要求。为此,Kalidind S R等[4]通过人为限定的方式约束了应力增量的大小,如式(4)所示。这就使得计算效率显著降低,难以应用于大变形过程的计算。


式中:n为迭代步数;ΔT*为应力增量;η为系数;s0为初始滑移抗力。

针对上述问题,Li H W等[10-11]将同伦延拓算法引入到晶体塑性模型的求解,构建了同伦模型(式(5)),发展了稳健数值算法。


式中:H为同伦函数; l为0~1之间的数;x0为初始值;α(l)为包含参数l并用来保证微分方程f′(x)非奇异的函数。式(5)的具体求解过程如下:

式中:k=1,2,…,N-1,N为迭代次数;h=1/N为迭代步长;δ为常数;ejj分量为1的单位向量。

同伦延拓算法是求解非线性方程组的无条件收敛算法,很好地解决了上述两个方面的限制问题。然而,同伦延拓算法的求解精度较低,因此,Li H W等[11]将其与N-R迭代相结合,以同伦延拓算法求得的解作为N-R迭代的初值,从而既保证了计算收敛,又保证了二次精度。该模型及算法成功应用于材料大变形的有限元计算过程中,可靠地预测了无氧高导铜(Oxygen-Free High-Conductivity Copper, OFHC)圆柱试样压缩过程中的不均匀变形(图1)和应力应变响应(图2)。进一步,应用于晶体塑性有限元模型(Crystal Plasty Finite Element Method, CPFEM),预测了OFHC多晶体试样在拉、压、剪、平面应变压缩条件下真应变达到1.4时的变形织构(图3)[12];嵌入Taylor型多晶体框架预测了100Cr6钢矩形截面环件轧制成形环坯内表面和中心处的变形织构(图4)[13]。预测的织构与试验测得的织构类型完全相符,表明了该模型具有良好的预测精度。

图1 OFHC铜柱镦粗过程的模拟及试验对照

(a)压缩比30%(b)压缩比50%

Fig.1 Comparison between experiment and simulation of OFHC copper cylinder under upset processing

(a) Compression ratio 30% (b) Compression ratio of 50% 

图2 OFHC铜在不同变形条件下的应力-应变响应及试验对照 

(a)单向压缩(b)纯剪切

Fig.2 Comparison of stress-strain responses of OFHC copper from simulations with those from experiments

(a) Uniaxial compression (b) Simple shear



图3 OFHC多晶体试样在拉、压、剪、平面应变压缩条件下真应变1.4时的变形织构预测及试验对照

Fig.3 Comparison between predictions and experiments of the deformation textures of OFHC grain aggregates under conditions of uniaxial tension, uniaxial compression, simple shear, plane strain compression under true strain of 1.4



图4 100Cr6钢矩形截面环件轧制成形环坯内表面和中心处的变形织构预测与试验对照

Fig.4 Prediction of deformation textures at center and internal surface of 100Cr6 steel ring with rectangular section and comparison with experiments


1.2 考虑多相材料热变形组织演化影响的晶体塑性建模


对于钛合金、高强钢等两相非均质材料,两相具有不同的晶格类型、滑移系及参数,不均匀变形异常显著。同时,在热加工过程中会产生动态再结晶等组织演化,导致新晶粒的形成、晶粒取向的突变等物理现象,显著影响后续的不均匀变形分布。因此,采用晶体塑性进行建模的关键在于:如何处理由于动态再结晶而大量产生的新晶粒;如何将组织演化引入位错滑移方程,以描述组织演化对变形的作用。

与动态再结晶产生大量新晶粒相似的问题存在于孪生变形。由于孪生使得原来的一个晶粒(即一个晶粒取向)变成了互为镜像的两个(即两个晶粒取向),其集中爆发会使得新晶粒大量涌现。这给晶体塑性建模带来了极大的挑战:因为晶体塑性模型是给每个晶粒赋予取向,并计算其变形响应的,爆发式增长的新晶粒取向会导致晶体塑性模型计算陷入瘫痪。解决孪生过程中大量新晶粒出现问题的有效处理方法有:优势孪晶取向方法(Predominant Twin Reorientation Method,PTR)[14]、体积分数转变方法(Volume Fraction Transfer Method, VFR)[15]等。这些方法的共同点在于将新晶粒与母相晶粒认为成一个晶粒,并引入等效系数或体积分数等以考虑新晶粒的作用,从而避免了晶粒个数剧增给计算造成的麻烦。

受此启发,Li H W等[16]提出了一种动态再结晶随机形核/球形长大的等效建模方法(图5a):动态再结晶只在边界处形核,形核时将材料点晶粒取向更新为新晶粒取向,新晶粒长大将占据其母相晶粒的空间位置。这样既避免了晶粒个数的增加,又较好地处理了新晶粒的长大问题,从而解决了对动态再结晶形核长大的描述难题。该工作中将滑移系剪切应变率γ·引入Ding R等[17]提出的形核率n·方程(式(8)),将位错密度演化引入模型计算,并采用基于位错密度的临界形核条件作为动态再结晶形核判据(式(9)),以描述晶粒不均匀变形对动态再结晶形核的影响;


式中:ρc为临界位错密度;ρSSD为存储位错密度;Qint为单位面积上的晶界能;ε·为宏观应变速率;b为柏氏矢量;l为位错平均自由程;M为晶界迁移率;Qline为位错线能量。

同时,采用基于存储位错密度(Storage Dislocation Density, SSD)与Hall-Petch效应的滑移系硬化方程(式(10)),以描述由于动态再结晶所引起的尺寸效应以及位错密度变化对后续晶粒变形的影响。基于此,建立了考虑钛合金等温变形动态再结晶影响的晶体塑性模型。然后,基于该模型预测了IMI834钛合金在等温成形过程中的应力-应变响应(图5b)及平均晶粒尺寸(图5c),并通过与试验对照证明了该模型的可靠性。

式中:c3为拟合参数;c4为Hall-Petch因子;d为平均晶粒尺寸。


图5 基于晶体塑性有限元法的IMI834钛合金等温成形动态再结晶建模

(a)形核长大的等效建模方法(b)模拟预测的应力-应变响应(c)模拟预测的平均晶粒尺寸

Fig.5 Modeling of dynamic recrystallization of IMI834 titanium alloy in isothermal deformation

(a) Modelingapproach of nucleation and grain growth (b) Prediction of stress-strain response (c) Prediction of average grain size 


值得注意的是,随着试样尺寸减小(如在微观性能表征、微成形过程中[18-19]),尺寸效应显著增长。为更准确描述尺寸对变形的影响,众多学者在模型中加入了对几何位错密度(Geometrically Necessary Dislocations,GND)的考虑[20-22]。


1.3 考虑多相材料晶界变形影响的晶体塑性建模


晶界具有不同于晶粒内部的物理属性。室温下,晶界对位错运动起强烈的阻碍作用,使位错塞积,从而在晶界处产生很高的应力集中。在高温下,一方面晶界由于温度的升高而变得活跃,塞积的位错会被晶界吸收,位错塞积引起的应力集中会被松弛。

图6 基于TA15钛合金真实微观组织的晶体塑性有限元建模仿真 

(a)晶界本构模型(b)几何模型(c)模拟结果(d)试验结果

Fig.6 Crystal plasticity finite element modeling and simulation based on real microstructure of TA15 titanium alloy

(a) Constitutive of grain boundary(b) Geometric model(c) Simulation results(d) Experiment results


针对这一问题,Li H W等[6]在CPFEM中引入缩放因子λ描述晶界中位错密度的变化(式(11)),以描述在高温变形过程中晶界的演化及其对组织演化的影响。

式中:c5、c6为拟合参数;ε·0为参考应变速率;θmelt为材料熔点;εp为材料点的等效塑性应变。

另一方面,高温下晶界变形更容易发生,从而成为材料变形损伤断裂的潜在最薄弱的区域,显著影响多晶体材料的宏微观性能。目前,有关晶界变形的研究工作,多关注晶界本身的特性,如粘滞性、率敏感性等。Wu Y Q等[23]针对细观尺度韧性多晶体材料,利用各向同性的Chaboche模型并结合Cohesive单元模拟了高温变形条件下断裂尖端区域率相关的晶界行为。Xu Q等[24]在弹塑性的Cohesive连续型模型中引入粘滞性软化或硬化因子来描述晶界的粘塑性变形,建立了考虑粘弹塑性的晶界变形模型。Chen S S等[25]在Wu Y Q等[23]和Xu Q等[24]的基础上,综合考虑晶界的粘滞性和率敏感性,采用Cohesive单元描述晶界,建立了粘弹塑性晶界变形损伤本构模型,如图6a所示。但是,以上研究均忽略了晶粒的不均匀变形对晶界变形的影响。针对该问题,Li H W等[26]在前期建立的粘弹塑性晶界本构模型的基础上,采用基于真实组织的晶体塑性有限元(CPFEM)描述晶粒的不均匀变形,如图6b所示,综合考虑了晶粒形貌、晶粒取向和相属性对晶界变形损伤的影响。对比图6c和图6d可知,模拟得到的微裂纹与试验结果相吻合,实现了对钛合金高温变形微裂纹萌生与扩展的准确预测。


2组织形态演化预测的晶体塑性与元胞机耦合建模


2.1组织形态演化预测的元胞机建模


元胞自动机法是一种时间、空间、状态都离散的网格动力学模型,利用简单的、局部的规则和离散的方法描述复杂的、全局的、连续系统的高效算法。元胞自动机在生物、交通、气象、材料加工组织演化等众多领域得到了广泛应用。

一个完整的元胞自动机模型系统通常包括元胞类型、元胞状态、邻居类型、元胞空间、边界条件及元胞转变规则6个部分。(1)元胞空间通常为二维和三维。(2)元胞网格类型主要基于元胞空间类型,其通常采用规则形状(如:二维元胞空间的三角形,四边形以及六边形),也可采用不规则形状、随机形状[27-28]。元胞形状可以恒定不变,也可随着元胞自动机的增量步发生变化[6, 29]。(3)每个元胞在每个增量步具有一定的状态(如晶粒位错密度),每个元胞的状态依据于邻居的元胞状态与转换规则进行演化。(4)由于元胞的转化规则只作用于当前元胞及与之相邻的元胞,因此对于一个完整的元胞空间,必须定义元胞的邻居类型。图7所示为二维元胞空间采用的Neumann与Moore邻居类型。(5)在元胞空间中设定边界条件的目的是以有限的计算机模拟区域研究无限的目标领域。元胞空间边界条件通常采用周期性边界条件、反射性边界条件及固定边界条件。(6)元胞转化规则是指由元胞当前状态及其邻居状况确定下一时刻该元胞状态的动力学函数,是元胞自动机法的“动力系统”,决定了整个系统如何演化。在元胞自动机建模过程中,如何基于物理机制制定元胞自动机转换规则是核心问题。如:Wu C等[30]基于溶质原子对界面扩散的拖拽效应,提出了钛合金热处理过程中考虑体扩散与界面扩散多机制控制的晶粒粗化的元胞自动机的建模方法,如图8所示。并基于此开发了数值化程序,在940 ℃保温5 h条件下,模拟获得了与试验结果吻合良好的晶粒粗化组织。

图7 元胞自动机邻居类型

(a)Neumann型(b)Moore型

Fig.7 Neighbor types in cellular automata

 (a) Neumann(b) Moore 


由于元胞不具有特定的物理意义,因此可根据不同的演化机理制定不同的演化规则,实现多种组织演化的模拟研究。德国学者Hesselbarth H W等[31]首次采用元胞自动机法研究了工艺参数对静态再结晶组织演化的影响规律,并预测了再结晶动力学模型,获得了与JMAK预测较吻合的结果。随后,多位学者基于元胞自动机法开展了对静态再结晶的影响的研究[32-33]。Liu Y等[34]基于晶粒的曲率驱动生长机制,首次通过元胞自动机法模拟了晶粒长大过程。同时,多位学者也基于元胞自动机法开展了对材料热加工过程中相转变组织演化的研究[35-37]。

Goetz R L等[38]首次将元胞自动机法引入到动态再结晶组织演化的模拟中,研究了单相合金在等温变形过程中初始晶粒尺寸对应力应变、再结晶形核率、再结晶动力学和加工硬化率的影响。但该研究没有涉及到具体工艺参数,而且假定动态回复是以恒定速率进行。随后,英国学者Ding R 和Guo Z X[39]利用元胞自动机法模拟了OFHC铜高温变形过程中不连续动态再结晶过程。建立了动态再结晶形核率模型、确定了再结晶发生的临界位错密度,并研究了工艺参数(如应变、应变率和温度)对组织演化的影响,但其忽略了材料不均匀变形对组织演化的影响。随后,多名学者基于元胞自动机模拟了金属材料在多步热成形过程中的动态再结晶演化,但仍没有考虑晶粒间的不均匀变形[29, 40-41]。



图8 基于元胞自动机法的TA15钛合金两相区加热多

机制晶粒粗化建模

(a)基于物理机制的数学建模(b)元胞自动机转化规则的制定

(c)940 ℃保温5 h条件下的模拟结果与试验结果对照

Fig.8 Modeling of grain coarsening in thermal two-phase area of 

TA15 titanium alloy based on cellular automata

 (a) Mathematical modeling based on physical mechanism

(b) Establishment of cellular automata switch rules

(c) Comparison between simulation and experiment under condition of 

5 h heat preservation in 940 ℃


2.2 考虑不均匀变形对组织形态演化影响的晶体塑性与元胞机顺序耦合建模


针对元胞自动机法难以考虑材料在热加工过程剧烈不均匀变形对组织演化影响的问题,将晶体塑性有限元法与元胞自动机法进行耦合,是准确预测金属高温变形过程中组织演化的有效途径。

Raabe D等[42]首次提出将晶体塑性有限元与随机性元胞自动机法相结合来模拟铝的静态再结晶过程。随后,Lan Y J等[43]在Raabe D研究基础之上,将耦合晶体塑性有限元的元胞自动机模型用于模拟奥氏体钢的相转变过程。Wu C等[2]将晶体塑性有限元与元胞自动机进行顺序耦合,其主要原理如图9a所示。通过晶体塑性有限元计算晶粒不均匀变形(图9b),并将晶粒的变形以及位错密度等信息传递给元胞自动机,作为输入条件进行组织演化。通过该顺序耦合模型可以预测不均匀变形条件下非均匀再结晶形核和长大等组织演变过程(图9c)。随后,Popova E等[44]在该工作的基础上,建立了镁合金高温成形与组织演化的顺序耦合预测模型,并解决了Wu C等[2]工作中晶体塑性有限元模型与元胞机模型之间需要网格信息对应的问题。


2.3 考虑不均匀变形与组织形态演化交互作用的晶体塑性与元胞机全耦合建模


在金属热加工过程中,不均匀变形与组织演变具有复杂的耦合影响机制,即:不均匀变形会影响组织演化,而复杂的组织演化同时也会影响后续变形。然而,上述晶体塑性有限元法与元胞自动机顺序耦合模型,虽然考虑了不均匀变形对组织演化的影响,但未考虑组织演化对后续变形的影响,变形与组织演化未能实现同步耦合预测。

针对该问题,Li H W等[6]提出了晶体塑性有限元与元胞自动机全耦合建模新思路,如图10a所示。其主要原理为:晶体塑性有限元与元胞机模型共用一个网格模型(即基于初始组织的三维有限元模型),在每个时间增量步中,晶体塑性模型计算晶粒尺度的不均匀变形,并将变形信息(如位错密度演变、晶粒取向)等传递给三维邻居类型的元胞自动机进行组织演化预测;反过来,元胞自动机将组织演化后的信息(如位错密度演变、晶界迁移等)传递给晶体塑性有限元进行下一步变形计算。这样考虑了组织形态演化与不均匀变形的耦合作用机制,最终实现了钛合金等温成形过程中不均匀变形与组织演化的全耦合预测(如图10b与图10c所示)。值得注意的是,已报道的元胞自动机建模工作中,增量步往往不具有物理时间,或通过最大晶界迁移速度与网格尺寸的比值来估算增量步时间[40]。而该模型中,元胞自动机与晶体塑性有限元计算是同步完成的,这使得该模型中元胞自动机的时间具有明确的物理意义,因而具有较高的预测精度。 

图9 基于晶体塑性有限元与元胞自动机的IMI834钛合金等温成形与动态再结晶组织演化顺序耦合建模

(a)建模思路(b)不均匀变形分布(c)动态再结晶不均匀形核与长大

Fig.9 Sequentially coupled modeling of dynamic recrystallization evolution in isothermal processing of IMI834 titanium alloy based on crystal 

plasticity finite element method and cellular automata

 (a) Modeling idea (b) Heterogeneous deformation distribution (c) Non-uniform nucleation and growth of dynamic recrystallization

图10 基于晶体塑性有限元与元胞自动机的TA15钛合金等温成形与组织演化多尺度全耦合建模

(a)建模思路(b)动态再结晶组织演化(c)不均匀变形分布 

Fig.10 Fully coupled modeling of dynamic recrystallization evolution in isothermal processing of TA15 titanium alloy based on crystal plasticity finite 

element method and cellular automata

 (a) Modeling idea (b) Dynamic recrystallization evolution (c) Heterogeneous deformation distribution

3 多尺度模型在钛合金复杂隔框等温成形过程中的应用


3.1 多尺度建模思路


钛合金整体隔框等温成形是多模具、多参数、多场耦合作用下的宏-细-微观多尺度变形过程。宏观模型无法预测大型构件细微观组织性能演化,而细微观模型又由于计算量太大难以负担大型构件成形过程的预测。基于宏观有限元+细观晶体塑性有限元+微观元胞自动机法的全耦合模型,通过采用变形梯度、位错密度等物理量,传递不同尺度间响应信息,可建立高效的难变形材料复杂构件不均匀变形与组织演化预测的机理型多尺度全耦合模型。其主要思路如图11所示:宏观有限元模型用于模拟与预测难变形材料的宏观成形;而以代表性体积单元为几何模型、并以宏观有限元模型中典型截面的变形信息为边界条件的细观晶体塑性有限元模型与微观元胞自动机全耦合模型,不仅可实现对钛合金复杂隔框等温成形过程细观尺度下复杂典型截面中的位错滑移、位错密度演化、织构演化的预测,同时可实现微观组织的动态仿真与预测控制;随后,细微观模型计算的织构等变形信息与微观动态再结晶等组织演化信息被返回给宏观有限元模型,为宏观条件下参数优化调控提供指导,最终实现难变形材料复杂构件成形的宏-细-微多尺度全耦合响应的精确模拟预测。


3.2 应用多尺度模型预测的钛合金复杂隔框等温成形的多尺度响应规律


基于上述多尺度仿真预测模型,并结合试验验证,揭示了局部加载和热处理条件对构件材料不均匀变形(应变分布)和组织演化(动态再结晶体积分数\[6]、晶粒尺寸分布[45]、等轴/片层α尺寸[46],等轴α体积分数[47]等)的影响规律,如图12所示。从而发展了通过非对称不等厚坯优化变形区和组合模具[48],进而采用两步热处理调控目标组织的优化成形工艺方案[49],最终实现了钛合金大型复杂构件整体精确成形,并在构件各个部位获得了均匀分布的性能优异的三态组织。


4 多尺度建模存在的问题及其发展趋势


图11 基于宏观有限元、细观晶体塑性有限元与微观元胞机全耦合的多尺度建模思路

Fig.11 Multi-scale modeling ideal based on microscale finite element method, mesoscale crystal plasticity finite element method and microscale cellular automata


图12 多尺度建模在钛合金复杂隔框等温成形过程中的应用

Fig.12 Applications of multiscale modeling on the isothermal forming of titanium alloy complex bulkheads


采用宏观模型无法预测大型构件细微观变形与组织性能演化,而使用细微观模型计算宏观大型构件成形过程将面临极大的计算量等问题。目前的宏-细-微全耦合多尺度模型采用从宏观模型中进行“取样”,减少了计算量,提高了计算效率。然而,“取样”条件尚不能完全反应复杂构件成形过程中典型截面的几何及变形历史,这使得细微观模型的变形条件难以与宏观有限元模型中的材料点变形信息完全吻合。同时,微观组织性能演化及其对本构的影响,尚未在宏观模型中得到反馈。这些问题的存在使得现有多尺度建模不得不割裂构件多尺度机制间的紧密联系。如何将细观模型与宏观模型的变形信息实现动态同步,并将微观组织性能演化信息实时反映到宏观变形响应中,以建立多尺度机制全耦合的统一多尺度模型,是多尺度模型未来发展的方向。采用均匀化方法来映射不同尺度间效应,是实现该发展目标的有效方法。

近几年来,Zhang X和 Oskay C[50]提出一种基于应变特征值的降阶均匀化方法,利用双尺度渐进展开将多晶体塑性分解为宏观和微观问题,有效解决了多晶材料组织演化过程中宏微观尺度同步预测的难题。Traxl R和Lackner R[51]提出了基质材料强度多阶层均匀化方法,用于预测多尺度材料体系的等效屈服行为。Zhou X Y等[52]提出了基于二阶扰动的随机多尺度均匀化方法,计算了复合材料的总体弹性力学性能。Wu L等[53]提出一种梯度增强均匀化方法预测了不同加载条件下碳纤维增强复合材料的损伤过程。Jia J等[54]提出了一个分阶层拓扑优化方法,通过引入宏观设计变量和微观设计变量来实现结构和多相材料的同步优化。由此可见,发展多尺度均匀化理论,建立大型构件成形的分阶层多尺度均匀化模型,实现大型复杂构件成形与组织性能同步耦合预测和协同调控,是未来多尺度建模的发展趋势。


5 结语


本文综述了晶体塑性求解算法、元胞自动机建模方法以及两者分别在材料成形与组织演化预测方面的应用,分析并综述了不均匀变形预测的晶体塑性建模、组织形态演化预测的元胞自动机建模及其耦合建模中的关键难题与解决方法。进而论述了基于宏观有限元、细观晶体塑性有限元与微观元胞自动机的多尺度全耦合建模方法及其在钛合金复杂隔框构件等温局部加载成形多尺度耦合响应规律预测中的成功应用。最后,分析了目前多尺度建模中存在的问题,并展望了多尺度建模未来的发展趋势。


参考文献:


[1]杨合, 孙志超, 詹梅, 等. 局部加载控制不均匀变形与精确塑性成形\[M\]. 北京: 科学出版社, 2014.


[2]WU C, YANG H, LI H W. Modeling of discontinuous dynamic recrystallization of a near-alpha titanium alloy IMI834 during isothermal hot compression by combining a cellular automaton model with a crystal plasticity finite element method[J]. Computational Materials Science, 2013,79: 944-959.


[3]GRUJICIC M, BATCHU S. Crystal plasticity analysis of earing in deep-drawn OFHC copper cups[J]. Journal of Materials Science, 2002,37(4): 753-764.


[4]KALIDINDI S R, BRONKHORST C A, ANAND L. Crystallographic texture evolution in bulk deformation processing of FCC metals[J]. Journal of the Mechanics and Physics of Solids, 1992,40(3): 537-569.


[5]DUMOULIN S, HOPPERSTAD O, BERSTAD T. Investigation of integration algorithms for rate-dependent crystal plasticity using explicit finite element codes[J]. Computational Materials Science, 2009,46(4): 785-799.


[6]LI H W, SUN X X, YANG H. A three-dimensional cellular automata-crystal plasticity finite element model for predicting the multiscale interaction among heterogeneous deformation, DRX microstructural evolution and mechanical responses in titanium alloys[J]. International Journal of Plasticity, 2016, 87: 154-180.


[7]ROSSITER J, BRAHME A, SIMHA M H, et al. A new crystal plasticity scheme for explicit time integration codes to simulate deformation in 3D microstructures: Effects of strain path, strain rate and thermal softening on localized deformation in the aluminum alloy 5754 during simple shear[J]. International Journal of Plasticity, 2010, 26(12): 1702-1725.


[8]LI H W, YANG H. A computational efficient and stable model of rate dependent crystal plasticity[J]. Steel Research International, 2011,SE(1): 814-818.


[9]HAREWOOD F J, MCHUGH P E. Comparison of the implicit and explicit finite element methods using crystal plasticity[J]. Computational Materials Science, 2007,39(2): 481-494.


[10]LI H W, YANG H, SUN Z C. Implicit algorithm and numeralisation of microscopic constitutive relations[J]. Materials Science and Technology, 2006,22(12): 1401-1408.


[11]LI H W, YANG H, SUN Z C. A robust integration algorithm for implementing rate dependent crystal plasticity into explicit finite element method[J]. International Journal of Plasticity, 2008,24(2): 267-288.


[12]LI H W, CHEN S S, YANG H. A new quantitative relation between deformation textures and corresponding mechanical properties of FCC metals[J]. Chinese Science Bulletin, 2014, 59(10): 1049-1056.


[13]LI H W, FENG L, YANG H. Deformation mechanism of cold ring rolling in view of texture evolution predicted by a newly proposed polycrystal plasticity model[J]. Transactions of Nonferrous Metals Society of China, 2013,23(12): 3729-3738.


[14]HOUTTE P V. Simulation of the rolling and shear texture of brass by the Taylor theory adapted for mechanical twinning [J]. Acta Metall, 1978,26(4): 591-604.


[15]LEBENSOHN R A, TOM C N. A self-consistent viscoplastic model: prediction of rolling textures of anisotropic polycrystals[J]. Materials Science and Engineering: A, 1994,175(1): 71-82.


[16]LI H W, WU C, YANG H. Crystal plasticity modeling of the dynamic recrystallization of two-phase titanium alloys during isothermal processing[J]. International Journal of Plasticity, 2013,51: 271-291.


[17]DING R, GUO Z X. Microstructural modelling of dynamic recrystallisation using an extended cellular automaton approach[J]. Computational Materials Science, 2002,23(1-4): 209-218.


[18]WANG Y, RAABE D, KLBER C, et al. Orientation dependence of nanoindentation pile-up patterns and of nanoindentation microtextures in copper single crystals[J]. Acta Materialia, 2004,52(8): 2229-2238.


[19]RAABE D, MA D, ROTERS F. Effects of initial orientation, sample geometry and friction on anisotropy and crystallographic orientation changes in single crystal microcompression deformation: A crystal plasticity finite element study[J]. Acta Materialia, 2007,55(13): 4567-4583.


[20]MA A, ROTERS F, RAABE D. A dislocation density based constitutive model for crystal plasticity FEM including geometrically necessary dislocations[J]. Acta Materialia, 2006, 54(8): 2169-2179.


[21]EVERS L P, BREKELMANS W A M, GEERS M G D. Non-local crystal plasticity model with intrinsic SSD and GND effects[J]. Journal of the Mechanics & Physics of Solids, 2004,52(10): 2379-2401.


[22]ARSENLIS A, PARKS D M. Modeling the evolution of crystallographic dislocation density in crystal plasticity[J]. Journal of the Mechanics & Physics of Solids, 2002,50(9): 1979-2009.


[23]WU Y Q, ZHANG K S. Crack propagation in polycrystalline elastic-viscoplastic materials using cohesive zone models[J]. Applied Mathematics and Mechanics, 2006,27(4): 509-518.


[24]XU Q, LU Z X. An elastic-plastic cohesive zone model for metal-ceramic interfaces at;finite deformations[J]. International Journal of Plasticity, 2013,41(2): 147-164.


[25]CHEN S S, LI H W, YANG H. Elastic-viscoplastic constitutive model of grain boundary deformation and damage[J]. Journal of Plasticity Engineering, 2014, 21(2): 13-19.


[26]LI H W, HUANG D, ZHAN M, et al. High-temperature behaviors of grain boundary in titanium alloy: Modeling and application to microcrack prediction[J]. Computational Materials Science, 2017, 140: 159-170. 


[27]JANSSENS K G F. Random grid, three-dimensional, space-time coupled cellular automata for the simulation of recrystallization and grain growth[J]. Modelling and Simulation in Materials Science and Engineering, 2003, 11(2): 157.


[28]FLACHE A, HEGSELMANN R. Do irregular grids make a difference? Relaxing the spatial regularity assumption in cellular models of social dynamics[J]. Journal of Artificial Societies and Social Simulation, 2001, 4(4):6.


[29]CHEN F, CUI Z, LIU J, et al. Mesoscale simulation of the high-temperature austenitizing and dynamic recrystallization by coupling a cellular automaton with a topology deformation technique[J]. Materials Science and Engineering: A, 2010, 527(21): 5539-5549.


[30]WU C, YANG H, LI H W, et al. Static coarsening of titanium alloys in single field by cellular automaton model considering solute drag and anisotropic mobility of grain boundaries[J]. Chinese Science Bulletin, 2012,57(13): 1473-1482.


[31]HESSELBARTH H W, GBEL I R. Simulation of recrystallization by cellular automata[J]. Acta Metallurgica Materialia, 1991,39(9): 2135-2143.


[32]GOETZ R L, SEETHARAMAN V. Static recrystallization kinetics with homogeneous and heterogeneous nucleation using a cellular automata model[J]. Metallurgical & Materials Transactions A, 1998, 29(9): 2307-2321.


[33]MARX V, REHER F R, GOTTSTEIN G. Simulation of primary recrystallization using a modified three-dimensional cellular automaton[J]. Acta Materialia, 1999,47(4): 1219-1230.


[34]LIU Y, BAUDIN T, PENELLE R. Simulation of normal grain growth by cellular automata[J]. Scripta Materialia, 1996,34(11): 1679-1683.


[35]SVYETLICHNYY D S. Modeling of grain refinement by cellular automata[J]. Computational Materials Science, 2013,77: 408-416.


[36]ZHENG C W, XIAO N M, HAO L H, et al. Numerical simulation of dynamic strain-induced austenite-ferrite transformation in a low carbon steel[J]. Acta Materialia, 2009,57(10): 2956-2968.


[37]HALDER C, MADEJ L, PIETRZYK M. Discrete micro-scale cellular automata model for modelling phase transformation during heating of dual phase steels[J]. Archives of Civil & Mechanical Engineering, 2014,14(1): 96-103.


[38]GOETZ R L, SEETHARAMAN V. Modeling dynamic recrystallization using cellular automata[J]. Scripta Materialia, 1998,38(3): 405-413.


[39]DING R, GUO Z X. Coupled quantitative simulation of microstructural evolution and plastic flow during dynamic recrystallization[J]. Acta Materialia, 2001,49(16): 3163-3175.


[40]KUGLER G, TURK R. Modeling the dynamic recrystallization under multi-stage hot deformation\[J\]. Acta Materialia, 2004,52(15): 4659-4668.


[41]CHEN F, QI K, CUI Z S, et al. Modeling the dynamic recrystallization in austenitic stainless steel using cellular automaton method[J]. Computational Materials Science, 2014,83: 331-340.


[42]RAABE D. Disrete mesoscale simulation of recrystallization microstructure and texture using a stochastic cellular automation approach[J]. Materials Science Forum, 1998, 273-275(2):169-174.


[43]LAN Y J, XIAO N M, LI D Z, et al. Mesoscale simulation of deformed austenite decomposition into ferrite by coupling a cellular automaton method with a crystal plasticity finite element model[J]. Acta Materialia, 2005,53(4): 991-1003.


[44]POPOVA E, STARASELSKI Y, BRAHME A, et al. Coupled crystal plasticity-probabilistic cellular automata approach to model dynamic recrystallization in magnesium alloys[J]. International Journal of Plasticity, 2015,66: 85-102.


[45]WU C, YANG H, LI H W. Modeling of static coarsening of two-phase titanium alloy in the α+β two-phase region at different temperature by a cellular automata method[J]. Chinese Science Bulletin, 2013,58(24): 3023-3032.


[46]吉喆. 基于建模的TA15钛合金组织—性能定量关系研究 [D]. 西安:西北工业大学, 2015.


[47]JI Z, YANG H, LI H W. Contributions of microstructural features to the integrated hardness of TA15 titanium alloy[J]. Materials Science & Engineering A, 2015, 628:358-365.


[48]ZHANG D W, YANG H, SUN Z C, et al. Deformation behavior of variable-thickness region of billet in rib-web component isothermal local loading process[J]. International Journal of Advanced Manufacturing Technology, 2012,63(1-4): 1-12.


[49]JI Z, YANG H, LI H W. Evolution of two types of α plates in tri-modal microstru-cture of TA15 alloy under varying processing conditions[J]. Rare Metal Materials and Engineering, 2015,44(3): 527-531.


[50]ZHANG X, OSKAY C. Eigenstrain based reduced order homogenization for polycrystalline materials[J]. Computer Methods in Applied Mechanics and Engineering, 2015,297: 408-436.


[51]TRAXL R, LACKNER R. Multi-level homogenization of strength properties of hierarchical-organized matrix-inclusion materials[J]. Mechanics of Materials, 2015,89: 98-118.


[52]ZHOU X Y, GOSLING P D, PEARCE C J, et al. Perturbation-based stochastic multi-scale computational homogenization method for the determination of the effective properties of composite materials with random properties[J]. Computer Methods in Applied Mechanics and Engineering, 2016,300: 84-105.


[53]WU L, NOELS L, ADAM L, et al. A multiscale mean-field homogenization method for fiber-reinforced composites with gradient-enhanced damage models[J]. Computer Methods in Applied Mechanics & Engineering, 2012,s233-236(4): 164-179.


[54]JIA J, CHENG W, LONG K, et al. Hierarchical design of structures and multiphase material cells[J]. Computers & Structures, 2016,165: 136-144.




我们把有用的知识、有趣的话题集中在这里


欢迎留言,提点建议,让我们变得更好 :)

联系我们:

Tel: 010-62912592

E-mail: sxgcxb@263.net

官方微信号:sxgcxb

长按关注