·今日铸造 Today Foundry

晶体塑性理论进展及晶体塑性有限元模型的应用

孙 进,张同先,王寅杰

(科华控股股份有限公司,江苏常州 213001)

要:晶体塑性作为突破宏观塑性均匀化局限性的重要理论手段,有助于关联微观变形与宏观性能,实现材料力学行为的精准预测。本文综述了晶体塑性变形运动学与晶体塑性变形基础,对比了唯象型晶体塑性本构模型与机理型晶体塑性本构模型的优缺点。根据单晶体塑性模型与多晶体塑性模型的本质及其联系引出均匀化模型,对全场晶体塑性模型与平均场晶体塑性模型进行介绍,列举了晶体塑性有限元模型在工程应用中的代表性实例。阐述了该模型应用的突破性进展,提出多晶体塑性模型的不足与晶体塑性理论的局限性,并对其未来的发展趋势进行了展望。

关键词:晶体塑性;塑性变形;本构模型;均匀化模型;有限元模型

在金属材料塑性成形过程中,微观组织演化、各向异性、循环加载、马氏体相变、变形孪晶与再结晶、损伤断裂疲劳等行为一直是塑性力学研究与工程应用的焦点。金属材料属于晶体材料,描述晶体材料的塑性理论可分为数学塑性理论与物理塑性理论[1]。虽然数学塑性理论较为简化,但其中的大多数参数是经验性的,无法通过实验获取。因此,当需要研究晶粒尺寸与大小,或者需要描述大变形塑性过程中的各种组织演变与材料力学响应行为时,物理塑性理论通常更为适用。

晶体塑性理论是物理塑性理论的分支。随着计算机硬件技术与微观组织建模软件的迅速发展,基于晶体塑性理论的晶体塑性模拟已广泛应用于材料、机械等工程应用领域,特别是结合有限元等偏微分求解器开展的晶体塑性有限元法(crystal plasticity finite element method,CPFEM)。这种方法可以模拟晶体材料在各种复杂的物理边界条件下的宏微观力响应以及宏观边界条件与微观组织演化之间的交互作用。大量科研人员近40年来对CPFEM进行了深度研究,扩展了其应用领域,并随着杜塞尔多夫高级材料模拟工具包(düsseldorf advanced material simulation kit,DAMASK)与粘塑性自洽模型(visco-plasitic self consistant,VPSC)等开源程序的发展[2-3],CPFEM对晶体材料塑性成形过程中的各种行为描述愈发准确。时至今日,晶体塑性理论已趋于成熟。本文旨在回顾前人对晶体塑性的研究,综述其发展历程,并对晶体塑性未来发展趋势进行展望。

1 晶体塑性变形理论

晶体塑性领域的开拓者当属Taylor[4]、Schmid[5]、Orwan[6]和Elam[7]等。1934年,Taylor发现可以用晶体薄层间的剪切变形来描述单晶的拉伸与压缩行为,于是将单晶体的塑性变形定义为位错在滑移系上运动引起的晶体切变。随后,Schmid提出Schmid定律,即当外力在某一滑移系中的分切应力达到一定临界值时,该滑移系首先发生滑移,这个分切应力称为滑移的临界分切应力,且只有当晶体满足上述条件时,才会发生塑性变形。Orwan等陆续提出了金属塑性变形的位错运动观点。在大量实验基础上,Taylor认为晶体内至少需要5个独立的滑移系才能实现其任意应变增量,首次提出了率无关本构关系,即材料的力学性能与载荷速率及变形速率都无关。基于这一理论,Hill与Rice[8]对塑性变形的几何学及运动学进行了详细的数学描述。在随后的几十年内,Needleman等[9]对晶体塑性理论的发展做出了较大贡献,其中当属发展了一套基于晶体塑性的率相关有限变形理论框架。率相关理论考虑了加载率的影响,并且很好地解决了滑移系判定问题,易于数值实现。下文主要介绍晶体塑性变形基础及其变形的运动本质。

1.1 晶体塑性变形基础

传统塑性理论认为,塑性变形是由位错在滑移系上运动导致。因此,塑性变形的本质是各向异性的,即与晶粒的取向密切相关。各种微观结构之间的相互作用将会对材料的流变行为和加工硬化产生影响[10],这些相互作用包括位错的增殖与湮灭、位错与晶界、亚结构、增强相或强化颗粒。塑性本构理论的核心任务是建立适当的数学模型来描述上述微观物理现象,从而表现出材料宏观的塑性流动行为与加工硬化。

面心立方结构A1(FCC)、体心立方结构A2(BCC)与密排六方结构A3(HCP)为金属材料最为常见的3种晶体结构。FCC晶体滑移系为{111}<110>,即FCC晶体具有固定的原子排列最密集的晶面与晶向。BCC晶体的原子密排程相对较低,滑移面可能为{110}、{112}与{123}晶面,具体受材料、温度等因素影响,并不具有最突出的原子排列和最密集的晶面,但其滑移方向是固定的,为<111>晶向,于是有些学者将BCC金属的滑移机制称为Pencil-Glide[11]。至于HCP晶体,其滑移方向通常为<11-20>,滑移面与其轴比(c/a)有关。当c/a≥1.633时,滑移面通常为{0001}晶面;当c/a<1.633时,原子密排面并不唯一,滑移可能发生在{0001}、{10-11}或{10-10}等晶面[12]。由于具有足够多的滑移系,满足了von Mises提出的变形协调准则[13],即宏观塑性变形至少需要5个激活的相互独立的滑移系,因此立方金属具有较好的塑性变形能力。对于HCP金属而言,在室温下只有基面<a>滑移或柱面<a>滑移容易激活,并且都只具有两个独立的滑移系,无法满足变形协调准则,因此HCP金属的室温塑性变形能力较差。当然,通过提高变形温度可以降低非基面滑移系的临界分切应力(critical resolved shear stress,CRSS),从而提高高温塑性变形能力。当材料处于高的静水压力或复杂的应力状态时,临界分切应力较高的非基面滑移更易被激发。常见的金属结构及其滑移系如表1所示。

表1 常见晶体结构及滑移系
Tab.1 Common crystal structure and slip system

BCC{112}<111>12 α-Fe、Zr Crystal Slip Slip Slip Common structure plane direction system metal FCC{111}<110>12 Cu、Al、Au、Ag{110}12{123}24{0001}<11-20>3(c/a≈1.633) Cd、Mg{10-10}<11-20>3(c/a<1.633) Ti HCP{10-11}<11-20>6(c/a>1.633) Mg、Zn{11-22}<112-3>6(c/a≈1.633) Ti、Zr、Cd

变形孪晶(机械孪晶)是HCP金属或低层错能的立方金属除位错滑移外另一种常见的塑性变形模式。绝大部分HCP金属都具有4种孪晶模式,每种模式都有6种对应的孪晶系[14]。孪生虽然常作为滑移不易进行时的补充,但其变形本质上却是一种位移型的快速切变,具有很快的形核和长大速度。由于产生的切变γtw取决于孪晶面法向、孪晶方向及c/a比值的几何常数,因此γtw也被称作本征切应变[15]。由于孪生能够改变晶体位向,晶体变形部分与未变形部分呈镜面对称关系,因此变形孪晶具有极性与方向性[16]。若要依据孪晶剪切变形对晶粒在<c>轴方向变形产生的特点将HCP金属的变形孪晶进行细分,则晶粒沿<c>轴方向被拉长或压缩对应的是拉伸孪晶与压缩孪晶。由于基面或柱面滑移均不能协调沿<c>轴方向的分切变形,故变形孪晶不仅为HCP金属提供更多的变形模式来满足von Mises所提出的变形协调准则,而且还提供了必要的沿<c>轴方向的分切变形[17]。孪晶形核所需克服的CRSS较大且通常在晶界等滑移受阻或应力集中的地方形核,因此孪晶形核受到层错能、晶粒大小与取向等方面的影响,即形核所需驱动力较大。然后,一旦孪晶形核完成,后续长大过程所需的驱动力将大大降低。即使孪生对塑性变形总量的贡献量通常较小,但变形孪晶会将某些不利于滑移的取向转变成有利于滑移的取向从而促进滑移的产生,因此孪生对于塑性变形的贡献依旧突出。

就目前而言,对于变形孪晶的微观物理机理研究仍存在较大的争议[18-19],现有的晶体塑性理论仍采用伪滑移的方式来考虑孪生。其最大的优势在于,可以允许在连续介质力学框架下,使用成熟的晶体塑性理论来统一描述滑移和孪生这两种不同的塑性变形机制,便于处理滑移与孪生之间的相互作用。当然,其最为突出的劣势在于,孪生具有方向性,在完整的伪滑移模型中,需要通过本构关系来引入这种方向性。除了考虑孪生的方向性,常用的流动模型及硬化模型几乎与滑移所用模型类似,认为其依旧遵循Schmid定理,甚至将孪生的CRSS用常数来表示[19-24]

当持续施加外力且晶体不发生滑移与孪生等塑性变形时,此时晶体将采用其他塑性变形方式。以镉单晶纵向压缩变形为例,如果外力恰好与该HCP单晶的底面(0001)晶面(六方最密堆积面)平行,此时取向因子为0,滑移系上的分切应力亦为0,晶体无法发生滑移。假设在这种情况下,孪生过程受到阻碍,若继续增加外力,当外力超过某一阈值时,为了适应外力,晶体将产生局部弯曲,形成扭折带,该区域晶体取向将发生不对称性变化,这种变形方式称为扭折。

1.2 晶体塑性变形运动学

大变形晶体塑性理论与宏观塑性理论类似,认为材料的变形可分解为弹性变形部分与塑性变形部分[25]。假设该材料上任意质点在tn时刻的位置为x,在tn+1时刻的位置为y,将其变形梯度设为F,则有:

材料由初始构型经过变形梯度转变为当前构型。晶体的变形梯度可以分为弹性变形梯度与塑性变形梯度,将该张量F进行乘法分解可得:

式中,Fe是晶格由外部载荷和位移(以及刚体旋转)的可逆响应引起的弹性变形梯度,也可以认为弹性变形梯度是由晶格畸变或晶体刚性转动引起的。Fp是位错在滑移系上运动从而对晶体产生的剪切所引起的塑性变形梯度[26],如图1所示。Fp是一种不可逆的永久变形,当产生变形的所有外力和位移都被消除时,这种变形仍然存在。从这个意义上说,Fp对初始构型的转换导致一个自由的不受外界应力影响的中间构型,并且该构型通常被认为可以保持完美的晶格[27]

图1 变形梯度示意图:(a)纯位错变形实例;(b)纯晶格变形实例[26]
Fig.1 Schematic diagram of the deformation gradient:(a) an instance of pure dislocation deformation;(b) an example of pure lattice deformation[26]

如图2所示,初始构型通过Fp转变为中间构型,中间构型又通过Fe转变为当前构型。

图2 构型转变方式[27]
Fig.2 Configuration transformation mode[27]

与分别为初始构型下晶体中第α个滑移系滑移方向的单位向量与滑移面法线方向的单位向量,sαnα分别为当前构型下晶体中第α个滑移系的滑移方向向量与滑移面法线方向向量。经过变形后sαnα已不是单位向量,但依旧保持正交。sαnα可通过下式表示:

变形速度梯度L通过下式表示:

由式(5)可得当前构型的塑性速度梯度Lp与中间构型的塑性速度梯度的表达式为:

源于滑移系(孪晶系)上的位错运动(变形孪晶)等塑性变形。当金属塑性变形的方式只有滑移时,设滑移系个数为N,则由滑移引起的速度梯度表示为:

式中,为第α个滑移系的滑移剪切率。

值得注意的是,当金属的塑性变形方式为滑移与孪生时,通过下式表示:

式中,ftw为孪晶体积分数,为第α个滑移系的滑移剪切应变率与滑移的Schmid张量,与为β个孪晶系的变形孪晶分切应变率与孪晶的Schmid张量。

当滑移系(孪晶系)上的分切应力超过其CRSS时,位错滑移(变形孪晶)将被激活[28]。其分切应力的表达式为:

式中,为滑移系(孪晶系)上的分切应力,Me为Mandel应力张量,表示为:

式中,Te为第二类基尔霍夫应力。对于金属材料,由线弹性胡克定律可得Te为:

式中,Ce为Intermediate构型的四阶弹性张量。Ee为Green-Lagrangian应变张量(物质应变张量),通过下式表达:

式中,I为单位向量。

通过研究晶体塑性变形运动学,可将晶体的塑性变形运动分解为弹性变形部分与塑性变形部分,通过数学方程的形式表达,为后续本构模型的建立提供数学基础。

2 晶体塑性本构模型

晶体塑性本构理论的核心为建立符合晶体塑性变形机理的塑性流动模型、加工硬化模型及内变量(孪晶分数等)演化模型。与传统塑性本构理论不同的是,晶体塑性本构理论需要考虑到滑移系的自硬化、潜硬化作用、滑移-孪生的交互作用、单晶体到多晶体的均匀化过渡方法等。经过多年的探索与总结,最终将单晶体塑性本构模型分为唯象型本构模型与机理型本构模型两类,并衍生出多晶体塑性本构模型。

2.1 唯象型晶体塑性本构模型

2.1.1 率相关流动法则

率相关流动法则是用来描述剪切应力与剪切应变率之间关系的本构模型。在与率相关的唯象型晶体塑性本构模型中,滑移系剪切应变速率通常通过幂指数法则直接与瞬时分解剪切应力相关联,对应的本构模型如下式所示:

式中,为第α个滑移系的剪切应变速率与参考剪切应变速率;为第α个滑移系的分解剪切应力;为第α个滑移系的滑移阻力(CRSS),表明滑移系的当前强度,即临界分切应力;n为应变速率敏感指数的倒数,当1/n→0时,材料的流动行为可以看做近似率无关。考虑到包申格效应,即循环软化效应,需要在晶体塑性本构模型中引入Yaguchi与Takahashi[29]在宏观本构中使用的背应力,如下式所示。

式中,χα为滑移系α的背应力,用来描述非线性运动硬化。背应力可分解为两部分,且各部分与率相关的本构模型,通过下式表示:

式中,C1C2d2为材料常数,d1是用来描述循环软化行为的变量。

为了用数值来实现孪生这一瞬间完成的塑性变形,通常用位错滑移的方式来处理孪生,即伪滑移模型。该模型与滑移模型的区别仅仅是其考虑了孪生的方向性。Doquet[30]和Lebensohn等[31]率先将孪生纳入晶体塑性理论中,并给出对应的孪生系剪切应变速率:

式中,与为第α个孪生系的剪切应变速率与参考剪切应变速率;为第α个孪生系的分解剪切应力;为第α个孪生系的激活阻力。根据目前的研究结果显示,在晶体塑性模型中,更多的关注点集中在孪晶的体积分数上[32],其率形式如下式所示:

式中,表示孪生从开始到结束,取值范围为0~1;γtw为孪生过程中发生的总分切应变。

2.1.2 硬化模型

由于加工硬化的影响,滑移阻力gtαw孪生阻力都随着塑性变形的进行而增大,具体体现在滑移系间的交互作用、滑移-孪生间的交互作用等方面。建立恰当的滑移系强度模型是正确描述加工硬化的前提。

硬化模型具体分为3类,第1类为Peirce-Asaro-Needleman(PAN)[33]硬化模型,第2类为Anand-Kalidindi[34]硬化模型,第3类为Voce[35]硬化模型。

前两类硬化模型通过下式表达:

式中,为滑移阻力;hαβ为滑移系交互强度矩阵;为第β个滑移系的剪切应变速率。

对于滑移系交互强度矩阵,也称之为滑移硬化模量,其表达式分别为:

式中,h0为硬化模量;gsg0为滑移系饱和强度与滑移系初始临界分切应力;α为硬化指数;q为潜硬化系数;δαβ为Kronecker算子。

之后Peirce等[33]改进滑移硬化模量本构模型,使其形式更为简单,其表达式为:

式中,hααhαβ为自硬化模量与潜硬化模量,分别表示滑移系开动对自身造成的硬化影响与滑移系β开动后对滑移系α造成的硬化影响。显然,自硬化模量依靠γ表征,而两种硬化之间通过潜硬化系数进行转化。γ表示所有滑移系累积剪切应变总量,可用下式进行表达:

第3类硬化采用饱和硬化模量进行表征,表达式为:

式中,为初始强度;为饱和强度;为初始硬化模量;为饱和硬化模量。

PAN模型是最经典的唯象硬化模型,但参数过多,硬化矩阵随着滑移系数量的增加而增加,计算效率过低。而Anand-Kalidindi模型克服了PAN模型的缺点,将硬化分为自硬化与潜硬化,引入状态变量来替代复杂的硬化矩阵。因此目前大多将PAN模型归纳为Anand-Kalidindi模型的初始状态,通过PAN模型去理解晶体塑性硬化的基础。Voce模型更多用于自硬化项,与其他模型联合使用,可以很好地描述FCC这类高层错能金属的饱和行为。对于立方金属而言,3类本构模型皆被广泛应用。而对于密排六方金属,为了考虑孪生的影响,通常会对这些本构进行改正。目前应用最广泛的为开源晶体塑性模拟器VPSC,采用的是Voce硬化模型描述位错滑移与孪生的强度,即CRSS。

2.1.3 晶格重取向

塑性变形会使晶体取向发生转变,从而导致织构的产生。大量研究表明,材料的性能20%~50%受织构的影响,即织构会导致弹性模量、泊松比等多种材料物理性能的变化。因此,描述晶体取向的变化是研究晶体塑性变形的基础。

针对滑移导致的晶体取向变化,通常用欧拉角描述晶体在坐标轴中的位置,其表达式[36]为:

式中,φ1ϕφ2为晶体由局部坐标系放入全局坐标系后需要转动的角度,即两坐标系间的夹角。

针对变形孪晶导致的晶格取向、变化,最初是由Van提出解决方式,归纳为晶格重取向(PRT)法,并由Tomé等[37]改进后引入VPSC中。该方法首先记录各个孪晶系的孪晶体积分数,并将其进行求和得到总体孪晶体积分数。当总体孪晶体积分数达到特定阈值时,孪晶的取向转变为具有最大体积分数孪晶的取向。尽管该方法忽略了体积分数较小孪晶的取向,但大量实验结果证实了其在数值模拟中的可行性。

2.2 机理型晶体塑性本构模型

唯象型晶体塑性本构模型虽然应用范围更广,但总体而言是一种经验模型,对塑性变形中微观物理现象的描述依旧不足。与唯象型本构模型相反,基于物理机理的本构模型依赖于内部变量。在塑性变形中,由于位错是塑性变形的载体,因此最重要的微观结构状态变量必定是位错密度。而以位错密度状态变量为根本建立的晶体塑性本构模型则称之为机理型晶体塑性本构模型。值得注意的是,即使位错密度是最重要的内部变量,仍需更多的参数来完全表征微观物理现象,例如晶粒尺寸、晶粒形状、第二相分数及沉淀物形态等[38]。随着DDB-CP模型的发展,材料模型的内容从位错增殖与湮灭引申到位错的扩散与攀移、位错与晶界或第二相粒子的交互作用、位错偶极、位错胞元亚结构[39-47]等复杂模型。

2.2.1 率相关流动法则

位错运动作为一个热激活过程,总会受到各种短程阻力(固溶强化强度阻力和晶格Peierls阻力等)与长程阻力(林位错/静止位错引起的阻力和晶界的弹性应力场引起的阻力等)的影响。通常用MTS模型[48-49]来描述位错运动引起的塑性流动,如下式所示:

式中,b为伯氏矢量大小;为运动位错密度;lα为位错运动平均自由程;wD为Debye频率;Gsp为位错运动激活能;kB为Boltzmann常数;τ0为短程阻力;为长程阻力。

由于长程阻力几乎不受应变速率和温度的影响,故也称之为绝热阻力。当塑性变形温度低于材料的熔点温度,即TTm时,式(28)可以用式(15)代替。通常而言,大多数与式(28)同源的热激活塑性流动模型总是用滑移系上的参考剪切应变速率来代替

2.2.2 硬化模型

在金属的塑性变形过程中,除了滑移系自身导致的自硬化与潜硬化外,不同滑移系上的位错之间也会存在明显的交互作用,导致长程阻力的出现。这种阻力会显著地影响位错的运动。目前,应用最为广泛的与位错密度相关的滑移阻力(硬化)模型当属Taylor模型[50],其表达式为:

式中,为类似于晶界引起的Hall-Petch效应[51]等微观结构引起的长程阻力;c1为数值常数;μ为滑移剪切模量;b为伯氏矢量大小;ρβ为滑移系β上的位错密度;Gαβ为滑移硬化模量系数(滑移系交互强度系数)。Gαβ的取值对于模拟结果的影响极大,目前虽无法准确给出该值确定方式,但已有一系列准确的经验数值可以进行参考。Nicaise等[52]也提出了采用更小尺度的离散位错动力学来确定滑移系交互强度系数。

2.2.3 位错密度演化模型

上文提到在塑性变形过程中,最重要的微观结构状态变量是位错密度,且基于物理机理的材料模型必定需要描述微观结构演化的模型,因此有关位错密度演化的模型一直是研究重点。在众多本构模型中,最基础的模型当属K-M模型[53],即:

式中,yc为位错湮灭的临界距离;lα为位错运动的平均自由程,其表达式为:

式中,d为平均晶粒尺寸;Lαβ为滑移系之间静止位错的贡献系数;c2K为数值常数。

K-M模型作为最基础的有关位错密度演化模型,其所考虑的变形机制自然不完全,并且无法描述温度对位错演化的影响。Peeters等[54]将胞元内静止位错密度、胞元墙内静止位错密度以及具有方向性的移动位错密度与位错密度有关的状态变量进行引入从而建立了可以描述BCC金属流动应力硬化/软化的晶体塑性本构模型。Bertin等[55]也提出了与位错密度相关的可描述温度效应的晶体塑性本构模型,且与VPSC共同被推广,并衍生出多种版本。

虽然与位错密度相关的机理型晶体塑性本构模型可以描述更多的微观组织、结构及变形机制,但科研者们对孪晶生长与变形机制的认知并不够深入,因此现在的孪晶模型大多还是唯象型的孪晶本构模型。例如Guo等[56]对TWIP钢建立的晶体塑性本构模型中,虽然描述滑移的本构模型是与位错密度相关的模型,但描述孪生的本构模型采用的是Kalidinid等[57]所提出的经验模型。

2.3 非局部模型

大多数文献中的本构模型都可以认为是局部模型,其将总体变形梯度分解为弹性变形梯度与塑性变形梯度,整个加载过程均可由该本构模型所描述,且对于多晶体的应力-应变曲线预测与织构预测皆被证明是有效的。然而,当模拟规模变小时,例如纳米压痕模拟[58-60]、微柱压缩模拟[61]及小尺度梁弯曲模拟[62],由于尺寸效应而使局部模型无法准确描述上述行为。

流动应力的晶粒尺寸依赖性最初是由Hall和Petch通过经验方程进行描述[63-64]。许多研究表明,较小尺寸晶粒的强化效果是由于晶界附近发生的非均匀塑性变形导致的。该现象可以解释为移动位错在晶界前面堆积,导致应力集中,增加了滑动阻力与晶界前面的应变梯度。并且这些应变梯度的增加也会导致位错密度的增加,从而增加滑动阻力[65]。除了晶粒尺寸效应,上述的微观尺寸的实验(压痕、压缩、弯曲)也揭示了流动应力的晶粒尺寸依赖性。在这些实验中,通常会发生不均匀塑性变形,这可能导致晶粒取向与应变梯度的变化,这些应变梯度的变化可以与几何必须位错密度(geometrically necessary dislocation,GND)相关联[66]。如何将GND的概念引入晶体塑性本构中或者说如何建立非局部晶体塑性本构模型就成为完善晶体塑性本构理论的重点。

在硬化方程(滑移阻力演化方程)中引入材料内禀尺寸与应变梯度效应是建立非局部晶体塑性本构模型最直接的方式。应变梯度理论的核心就是将位错密度纳入支配材料行为的本构方程中。也就是说,非局部晶体塑性模型既能考虑晶粒取向的影响、材料各向异性的影响与局部变形的影响,又可以考虑应变梯度(位错密度)的变化对硬化的影响。

最初由Arsenlis、Gurtin、Evers等[67-72]将几何必须位错密度纳入晶体塑性理论框架,再由Han等[73-74]将其发展,确定了滑移抗力的演化形式,如下式所示:

式中,第1项与式(20)所述相同,第2项则描述几何必须位错引起的应变硬化。k0为材料常数;b为伯氏矢量大小;G为剪切模量;为无量纲参数,一般取值为1/3;为滑移系上的几何必须位错密度。

为了确定几何必须位错密度的演化模型,需要引入Nye张量,其表达式为:

利用式(33)与式式(5~7)可以将几何必须位错密度的演化模型导出[75],如下式所示:

总的来说,将GND引用到本构模型中就是将几何必须位错投影到静止位错与平行位错中[76]。为了使该投影行为更方便,可以将分解为3组:第1组为切向矢量平行于滑移方向dβ的螺位错,如式(35)所示;另外两组为切向矢量分别平行于与的刃位错,如式(36)和(37)所示:

式中,tβ为刃位错的方向矢量。上述3式满足下式:

近年来,研究者们对位错密度演化模型的应用逐渐增多,其中最具代表性的当属将位错密度参数纳入硬化模型中,从而描述材料循环应变硬化软化及应力松弛行为。人们普遍认为,在随动硬化规则中加入静态恢复项是模拟循环加载下材料随时间变形的必要条件。一方面,具有恒定静态恢复项的传统硬化规则不能再现循环应力松弛行为[77];另一方面,位错间的交互作用以及位错与沉淀相的交互作用是应力松弛中的主要变形机制,这也是修改随动硬化本构模型的基础[78]。2023年,Pin等[79]在A-F经典模型[80]的基础上,提出了一种修正的非线性运动硬化准则,并引入了一个新的静态恢复项,用来预测镍基单晶高温合金在循环载荷下的蠕变、循环应变硬化软化以及应力松弛行为,其表达式为:

式中,C3为建立的静态恢复系数,与林位错密度有关;r0φs为会恢复过程中的材料常数;ρr为参考位错密度,用于控制应力松弛随载荷循环次数衰减速率;ρfor为不同滑移系上林位错的密度。式中引入的衰减函数C3类似于Xuan模型[81]中的静态恢复项,但不同点在于,该式直接将林位错密度与位错数量替代所累积的非弹性应变。当φs=1时,模型退化为经典的A-F模型。

唯象模型大多关注“是什么”这一概念,追求用简单的数学方程高效的复现出宏观实验,在金属成形过程预测方面有较为广泛的应用。机理模型核心在于“为什么”,更注重于将位错动力学等微观机理引入本构关系中,大多应用于需要研究变形物理机制的场景,例如离散位错动力学模拟。非局部模型的优势在于对“尺度效应”的把控,致力于解决局部模型在微纳尺度下的失效问题,应用场景大多是细晶材料的变形。

3 均匀化模型

多晶体是由不同取向的单晶体堆垛组合而成,构成多晶体的单晶体的总和的塑性响应均匀化之后就可以得到多晶体的塑性行为,由此衍生出多晶体均匀化模型。多晶体均匀化模型,作为预测多晶和多相材料弹塑性响应的一种方法,连接了介观尺度与宏观尺度的变形情况。运用不同的均匀化方法,将会得到不同的多晶体均匀化模型。目前所广泛应用的均匀化模型为平均场晶体塑性模型与全场晶体塑性模型两类。

3.1 平均场晶体塑性模型

平均场晶体塑性模型主要分为3类,分别为Taylor模型、Sachs模型以及自洽模型。Taylor模型[82]在1938年被Taylor提出,该模型认为多晶体的宏观应力等于所有单晶体应力总和的平均值,单个晶粒的应变等于多晶体的宏观应变,因此也称为全约束模型(FC模型)。该模型假设多晶体各个晶粒的应变分布是均匀的,且多晶体宏观应变量与其中每个单晶体的应变量高度吻合,因此多晶体内的所有晶粒都满足应变相容性的要求。该模型自动满足晶粒间的变形协调条件,但破坏了晶界处的应力平衡条件。但Taylor模型在描述HCP金属与低层错能的立方金属的织构演化与各向异性预测方面表现不佳,因此出现了许多集中于晶界应力平衡的模型,例如松弛约束模型(RC模型)[83-85],其中最经典的当属“煎饼模型”。RC模型不再执着于晶粒间应变相容性的要求,因此改善了晶粒间的应力相容性,所受认可度较高。时至今日,Taylor模型依旧被大部分人所认可,并应用于晶体塑性模拟中。

Sachs模型是在1928年被提出的一种最简单的多晶体塑性模型。这一模型假设多晶体各晶粒的应力分布是均匀的,且多晶体宏观应力大小与各晶粒应力大小一致,即多晶体内所有晶粒都满足应力平衡条件,但该模型不满足应变协调条件。因此该模型早期只对处理单晶体的塑性变形问题有效,始终无法描述多晶体的塑性变形过程。即使后人经过修正与完善,但其应用程度远不及Taylor与自洽模型。

1954年,Hershey[86]提出自洽模型(self-consistent theory),后由Kronor、Budiansky、Wu、Mori等进一步应用到多晶体塑性模型。自洽晶体塑性理论(SC-CP)认为多晶体是一种均值等效介质(homogeneous equivalent medium,HEM),每个晶粒视为一种椭球形黏性夹杂物,因此需考虑夹杂物形状与局部领域等问题。最初,Kronor认为晶粒作为一种具有材料整体刚度的孤立夹杂物存在于介质中,之后,Mori将晶粒作为黏性夹杂物直接嵌入进介质中,通过格林函数解决夹杂问题[87]的同时,满足应变协调与应力平衡条件,并将平均基体应变看做成整个材料的远场应变。为了将总应力与应变率结合起来,Molinari等[87]采用切线模量法建立了非增量黏塑性自洽模型,Lebenso-hn与Tomé在此基础上,建立了大量黏塑性变形(SC-CP)模型。当前应用最广泛的SC-CP模型当属黏塑性自洽(visco-plastic self consistant,VPSC)晶体塑性模型,大量研究结果表明,VPSC虽未考虑弹性变形,但对于HCP金属与低层错能立方金属的力学行为描述与各向异性演化均精确有效。

3.2 全场晶体塑性模型

平均场晶体塑性模型可以很好的预测多晶体的应力-应变曲线(宏观力学行为)与织构演化行为,但难以预测晶粒的微观力学行为[88]。全场晶体塑性模型借助于偏微分方程求解器,例如有限元法(finite element method,FEM),通过网格将晶粒进行离散,各个晶粒需要一组有限单元进行确定,如图3所示,一种颜色代表一个取向或一个晶粒,每个晶粒被n个网格所离散。

图3 全场多晶体塑性模型
Fig.3 Full-field polycrystalline plastic model

通过这种全场建模方式,单元与晶粒共用边界,采用平衡(动量)方程,依旧可以满足晶粒间变形协调条件与应力平衡条件,且可以考虑到材料的微观组织情况。网格划分越细,模拟结果越精确,计算时间随之增长且对计算机的硬件要求较高。当然,这种对每个晶粒进行网格划分并开展模拟的方式,其计算时间比平均场晶体塑性模拟所需的时间至少高出一个数量级。因此,该方式并不适用于大变形或大尺度规模的有限元模拟。目前最主流的全场晶体塑性方法为全场CPFEM与全场晶体塑性谱方法(crystal plasticity spectral method,CPSM)[89]。近年来,随着建模技术与PC硬件的升级,全场晶体塑性模型的应用越来越广泛。Jialili等[90]提出了一种结合计算均匀化技术的空间分辨CPFE模型,用于模拟挤压ZK60镁合金循环塑性变形的微观机制。Liu等[91]提出了一个适用于具有强冷轧织构的AA3104-H19与AA2024-T3铝合金的多尺度模型,使用快速傅里叶转变的全场晶体塑性谱方法,预测在平面内和平面外的单轴与多轴应力状态下的各向异性力学性能。

均匀化模型作为一种跨尺度模型,可作为桥梁连接材料的微观与宏观状态。该模型认为微观结构是统计均匀的,即一个足够大RVE单元能够代表材料的微观结构特征。当确定材料宏观性能后,可以合理的运用均匀化理论进行逆向设计,寻找最优的微观结构,这是拓扑优化的基石。然而,对于CPFE模型而言,RVE单元需要较为准确的网格划分,尤其是当晶界附近的网格划分出现偏差时,材料整体的力学性能会有明显的变化。过细的网格带来的超强度计算负荷与过粗的网格导致的计算结果失准,两者很难做到平衡。

4 晶体塑性有限元模型的应用

CPFE模型目前已广泛应用于各种塑性成形过程的模拟、织构预测、各向异性预测与微观组织演化中,并且可以与其他介观组织模拟方法相结合研究变形过程中的组织演变,涵盖了多尺度机械力学行为、材料基础冶金行为及金属成形工艺等领域。这里仅叙述若干具有代表性的有限元模型在工程应用领域方面的实例。

4.1 晶体塑性有限元模型在塑性成形模拟方面的应用

晶体塑性模型的本质是源于晶体的运动,因此晶体的初始取向对材料成形性能的影响极大。通过全场CPFE模型,赋予各个晶粒初始取向,在载荷的作用下,可以直观地体现各个晶粒间的变形协调状况及晶体取向的变化,从而达到预测成形质量的目的。

筒形件拉深制耳是晶体取向影响板材成形质量的重要体现。Raabe等[92]通过织构组分CPFEM研究了不同织构组分及其组和对铝合金筒形件拉深制耳行为的影响。5个离散织构组分Cube、S、Brass、Copper和Goss所引起的制耳的轮廓如下图4所示。模拟结果显示,S、Brass与Copper织构组分所引起的制耳的耳峰出现在与轧制方向呈45°与135°的位置,Cube与Goss织构组分引起的制耳的耳峰出现在轧制方向与其横向附近,即与轧制方向呈0°与90°。

图4 离散织构组分引起的制耳:(a) S、Brass与Copper织构组分引起的制耳;(b) Cube与Goss织构组分引起的制耳[92]
Fig.4 Earing caused by discrete texture components:(a) earing caused by S,Brass and Copper texture components;(b) earing caused by Cube and Goss texture components[92]

研究发现,不同织构组分的组和对拉深制耳行为有极大的影响。针对于如何通过组合不同织构组分控制制耳的大小实现制耳最小化这一问题,Raabe等[92]也进行了大量的模拟实验,如图5所示。结果显示,当S织构组分与Cube织构组分的体积分数比例为1.67∶1.00时,可以得到最小化的制耳。

图5 不同比例的S与Cube织构组分的组和对制耳影响:(a)不同比例的S与Cube织构组分的组和对耳廓的影响;(b)耳高与耳峰的定量分析[92]
Fig.5 Differentproportions of S andCube texture components andtheirinfluence on earing:(a)S withCubeindifferentproportions of the composition of the texture components and the effect on the profiles;(b) quantitative analysis of ear height and ear peak[92]

制耳的优化取决于晶体织构与织构组分的球形散射度。图6a展示了Cube织构组分的初始球形散射度对耳廓的影响;图6b展示了62.75%的S织构组分与37.25%的Cube织构组分组合的初始球形散射度对耳廓的影响。可以看出,制耳的高度随着织构组分初始球形散射度的增加而降低。

图6 织构组分的初始球形散射度对耳廓的影响:(a) Cube织构组分;(b)62.75%的S织构组分与37.25%的Cube织构组分的组合[92]
Fig.6 Effect of initial spherical scattering of texture components on the ear profiles:(a) Cube texture components;(b) combination of 62.75%of the S texture component and 37.25%of the Cube texture component[92]

Li等[93]提出了一种率相关本构模型并应用到VUMAT子程序中,通过显示积分算法,模拟了圆棒料镦粗成形的过程,如图7所示。将3项结果进行对比,证明了该算法在多单元、大步长、大变形的模拟中的稳定性与可靠性。Kalidindi等[94]模拟了等通道挤压成形过程,如图8所示。其中,下部路径处的预测结果是迄今为止等通道挤压成形预测研究中最准确的结果。Kupka等[95]通过CPFE模型,模拟了缺口悬臂梁弯曲变形行为,研究了脆性晶界(Al-Li合金中)的损伤、断裂性能。

图7 基于RDCP与J2塑性的圆棒镦粗实验与模拟结果[93]
Fig.7 Round bar upset process,experiments and simulation results based on RDCP and J2plasticity[93]

图8 等通道挤压成形过程:(a)距模具顶端在模具厚度10%处的上路径对应的变形网格;(b)中路径对应的变形网格;(c)距模具顶端在厚度90%处的下路径对应的变形网格[94]
Fig.8 Equal channel angular pressing process:(a) deformed mesh corresponding to the upper path at 10%of the die thickness from the top;(b) deformed mesh corresponding to the middle path;(c) deformed mesh corresponding to the lower path at 90%of the thickness from the top[94]

4.2 晶体塑性有限元模型在织构预测方面的应用

金属材料在冷热加工或塑性成形过程中,多晶体中的各个晶粒会沿着某些方向有序排列,这种现象称之为织构。当出现织构不均匀现象时,会导致材料在工业加工过程中产生某些不可控因素,从而使材料的性能无法达到预期目标。Tikhovskiy等[96]通过织构组分CPFEM模拟了X6Cr17铁素体不锈钢在冷轧过程中的表层织构、中心织构及样品整体的织构演化,如图9所示。将模拟结果与实验结果进行对比,可以发现,在小应变的情况下有很好的一致性;当下压量达到70%以上时,模拟与实验结果偏差变得明显,实验结果如图10所示。

图9 模拟X6Cr17铁素体不锈钢冷轧过程中的织构演化:(a1~g1)表层织构演化;(a2~g2)中心层织构演化;(a3~g3)整体织构演化[96]
Fig.9 Simulation of texture evolution during cold rolling of X6Cr17ferritic stainless steel:(a1~g1) surface texture evolution;(a2~g2) central layer texture evolution;(a3~g3) overall texture evolution[96]

图10 X6Cr17铁素体不锈钢冷轧过程中的织构演化:(a~c)表层织构演化;(d~f)中心层织构演化[96]
Fig.10 Texture evolution during cold rolling of X6Cr17ferritic stainless steel:(a~c) surface texture evolution;(d~f) central layer texture evolution

Tikhovskiy等[96]将结果进行对比分析,认为模拟与实验的结果偏差主要分为两点:①主要体现在织构的连续性上。在表层与中心层的织构中,实验织构相对模拟织构具有更强的连续性,尤其是在{001}<110>织构组分上更为明显。这种情况可能是由于这种织构组分的再结晶与旋转行为非常缓慢。②体现在织构强度上。模拟结果显示,在大应变的情况下,{111}<112>织构组分强度强于{111}<110>织构组分的强度,但实验揭示了相反的趋势。这一效应可能归因于两种原因:①是由于实验观察到的织构并不够完整;②是由于在建模过程中,相邻织构层之间的织构梯度过于简化。

Anand等[97]模拟了FCC金属高温压缩与室温压缩过程中的织构演化行为。Li等[98]通过一种显示积分算法,实现了对不同变形状态下织构演变的预测,如图11和12所示。Balasubramanian等[99]模拟了多晶纯钛在高温变形下的织构演化。Riyad等[100]将场波动与黏塑性自洽多晶塑性模型相结合的FF-VPSC模型应用到Al-5182-O的变形再结晶织构演化预测中。分别在轧制模拟中使用VPSC模型与FF-VPSC模型,并将得到的Cu织构进行对比(图13),可以看出,使用标准VPSC模型所得到的织构相较FF-VPSC模型缺少了(112)[11-1]方向的强度。

图11 单轴压缩过程中的织构演化:(a1~a3)模拟结果;(b1~b3)实验结果[98]
Fig.11 Texture evolution during uniaxial compression:(a1~a3) simulation results;(b1~b3) experimental results[98]

图12 单轴拉伸过程中的织构演化:(a1~a3)模拟结果;(b1~b3)实验结果[98]
Fig.12 Texture evolution during uniaxial stretching:(a1~a3) simulation results;(b1~b3) experimental results[98]

图13 轧制过程中的Cu织构演变:(a~e)使用VPSC模型所得的结果,下压量分别为20%、40%、60%、80%、95%;(a1~e1)使用FF-VPSC模型所得结果,下压量分别为20%、40%、60%、80%、95%[100]
Fig.13 Cu texture evolution during rolling:(a~e) results obtained using the VPSC model. The pressures are 20%,40%,60%,80%and 95%,respectively;(a1~e1) results obtained using the FF-VPSC model. The pressures are 20%,40%,60%,80%and 95%,respectively[100]

4.3 晶体塑性有限元模型在介观组织模拟方面的应用

为了描述材料在外力或热力作用下的变形与组织演化机理,针对介观状态下的晶体学进行了大量研究。动态再结晶(dynamic recrystallization,DRX)是一种金属在热变形过程中发生的再结晶组织演化现象,通过是否观察到新的等轴晶粒来判断有无发生DRX现象。DRX的起始伴随着新晶核的出现,假定形核部位为晶界处,形核速率则取决于温度与随后的晶粒生长速度。DRX会通过降低位错密度来降低系统能量,并最终影响热加工行为。近些年已进行了大量的实验工作来研究合金热变形过程中的DRX现象。当考虑到DRX对合金流动行为的影响时,通常将DRX与成形温度、应变以及应变速率联系起来[101-102]。为了节约研究成本,大多数情况下都会采用数值模拟的方式进行探索研究。基于实验数据的唯象本构方程通过开发应用到DRX动力学中,这类本构方程也被集成到有限元法中去预测热成形过程中合金的再结晶体积分数与平均晶粒尺寸[103]。显然,这些唯象方程无法捕获晶粒的局部变形信息与直接影响金属材料力学性能的形态和拓扑特征。

Li等[104]针对两相钛合金(α+β)等温热变形过程提出了一种考虑动态再结晶的晶体塑性本构模型。由于新晶粒的生成、新晶粒取向、位错密度与滑移阻力的减小、再结晶晶粒体积分数的演化等因素,在变形拓扑影响DRX的同时,DRX也会影响后续变形,因此该模型将这两方面进行耦合,从而满足同步考虑的要求。模拟与实验结果如图14所示,应力响应的不同揭示了这两个区域的软化机制不同。动态回复与再结晶是钛合金的两种重要软化机制。动态回复会消除位错结构,诱导位错的增殖,并导致位错密度的平衡。因此,动态回复会使应力保持一个稳定状态。在DRX过程中,会产生大量的新晶粒,它们会占据一部分基体,并使所占据部位的位错结构全部消失,所以DRX使等效位错密度降低,从而使应力达到峰值。

图14 α相区与α+β相区实验与预测应力对比[104]
Fig.14 Comparison of the α phase region and α+β phase region experimental and predicted stresses[104]

晶界迁移是位错密度演化的结果。M晶粒(基体晶粒)与R晶粒(再结晶晶粒)的位错密度差异为晶界迁移提供迁移力,Li等[104]针对不同温度与应变速率下M晶粒与R晶粒的位错密度演化进行模拟,如图15所示。对于M晶粒,由于较强的应变硬化与动态回复效应,位错密度在小应变的情况下直接增加到接近饱和值;而对于R晶粒,由于初始位错密度很低,位错密度先开始增加,随后积累到饱和值。位错密度在β相区较低,在α+β相区较高,这是由于β相区温度更高从而导致更强烈的动态回复。然而,无论实在β相区还是在α+β相区,在两个不同的温度下都可以发现几乎相同的位错密度。

图15 R晶粒与M晶粒位错密度的演化[104]
Fig.15 Evolution of the R grain dislocation density and M grain dislocation density[104]

该模型还捕捉到α相的作用。虽然α晶粒没有发生再结晶,但其对β晶粒DRX行为的影响有两点:①α相界可作为形核部位;②也是最重要的一点,α-β相界处存在较大的剪切应变,会促进R晶粒的形核与长大。

在晶体塑性领域中,针对介观组织模拟方面的研究成果显著。Lan等[105]将CA(元胞自动机)法与CPFE结合,模拟了奥氏体向α-Fe转变的过程;Xiao等[106]在此基础上,结合蒙特卡洛法与CPFEM,研究了奥氏体变形后对等温奥氏体向α-Fe转变的影响;Vondrous等[107]将CPFEM与相场法耦合,模拟了低碳钢在热轧-冷轧-退火过程中的静态再结晶行为;Popova等[108]将概率元胞自动机法与CPFEM进行耦合,建立了一种动态再结晶模型,如图16所示。Park等[109]将非均匀变形的CPFEM与模拟微观组织演变的概率元胞自动机方法及动态再结晶模型相结合,建立了一个基于显微组织的多尺度框架,如图17所示,并对AISI 304LN不锈钢在热机械过程中的DRX行为进行了模拟。CPFE模型描述了多晶材料变形时的非均匀晶粒度以及晶体取向,CA模型描述了DRX过程中的组织演变。Shah等[110]提出了一种耦合基于全场位错密度的晶体塑性模型和元胞自动机模型的模拟方法,并预测了热机械加工中的亚动态再结晶行为,如图18所示。其中晶体塑性模型采用的是DAMASK中的位错密度模型,CA模型基于CASIPT三维元胞自动机程序,并将其推广到应用于奥氏体亚动态再结晶过程中晶界迁移的描述。

图16 晶核长大过程[108]
Fig.16 Process of crystal nucleus growth[108]

图17 耦合CPFEM-CA法的多尺度框架[109]
Fig.17 Multiscale framework for the coupled CPFEM-CA approach[109]

图18 热轧过程中的再结晶体积分数演化(层间退火后):(a)下压量30%;(b)下压量44%;(c)下压量60%[110]
Fig.18 Recrystallization volume fraction evolution during hot rolling after interlayer annealing:(a) pressure of 30%;(b) pressure of 44%;(c) pressure of 60%[110]

4.4 晶体塑性有限元模型在宏观-微观模拟方面的应用

应变硬化被认为是最重要的强化机制之一,在评估材料的功能特性方面起着关键作用,并决定了其对工程-结构应用方面的适用性。识别晶体塑性硬化参数最常规的做法是在单晶上进行简单的机械实验,但由于大多材料都是多晶材料,并且对于多晶体中的单个晶粒,其嵌入后边界条件也会随之改变,因此很少被实践。最好的解决方法当属基于最小二乘法或简单试错法的逆分析技术,该技术的主要作用是校准材料参数,直到模拟结果与实验结果相匹配。然而这两种方法都需要初始估计值,并涉及大量的迭代,努力收敛到一个最佳值。与此同时,这些技术反映出参数预测的最大问题:不确定性,这是由于唯象本构中存在许多经验参数或未知参数,这导致一个参数的变化也可通过另一个参数的变化进行补偿。后续,一些研究人员开发的纳米压痕或数字图像技术(digital image correlation,DIC)等方法就是基于上述问题进行的优化算法的逆分析技术,其核心为建立了优化硬化参数值。

然而,在纳米压痕实验下,压头处于复杂应力状态下,很难通过实验获得控制硬化行为的参数,并且由于需要在几个晶粒上进行实验,成本过大。在DIC实验下,需要进行大量的有限元计算。Sahoo等[111]为了在更高效低成本的基础上确定硬化函数中的不确定参数的近似值,提出了一种新的基于幂律关系的宏观应变硬化函数的解析形式,如式(42)所示。不同于其他经验幂律函数,这个函数是唯一的,具有驱动流动应力饱和的饱和应力。

式中,为宏观von-Mises等效应变率;ha为硬化参数;σsat为饱和应力;为材料参考强度,出现在黏塑性行为的应变幂律公式中,其表达式为:

式中,m为应变速率敏感性指数;为参考应变速率。

将式(42)进行分离变量并积分可以得到下列形式的硬化函数:

式中,σ00为材料的屈服强度(时)。除了材料的屈服强度以外,其余3个未确定的值:σsatha都需要使用合适的拟合过程从实验中获取。

Sahoo[111]通过一种最小化技术,大大减小了预测参数的周期与成本,其表达式如下式所示:

式中,N为应力应变曲线上选取的坐标点个数,代表式(44)所得到的应力值与实验所得的应力值。显然,最小化技术的目的是使得f尽量趋于0。实现方式为:给与3个不定参数一个取值区间,进行计算,识别出最小的f值。之后在减小的区间中继续进行迭代,循环若干次后,可达到最小值f所需的精度。该程序计算速度极快,在几秒钟就可以完成收敛。

由于上述所以出的宏观应变硬化函数是由晶体塑性硬化本构模型简化得来的,因此可以用该宏观参数作为晶体塑性模拟的初始值。而宏观与微观硬化参数值之间的差异是由于力学行为尺度的差异,因此必须将这两种尺度联系起来。Sahoo等[111-112]通过Taylor因子将宏观与微观参数进行转换,转换方式如下所示:

式中,τ0是滑移系分切应力;τsat是饱和应力;hcac是硬化参数;M为泰勒因子。

上述研究方式先将微观硬化模型简化成宏观硬化模型,得到的硬化参数可用于预测晶体塑性硬化参数的初始值,并且通过Taylor因子这一桥梁进行宏观-微观的参数转化。

除了针对结构材料的晶体塑性研究外,对功能材料的研究依旧火热。烧结银纳米颗粒(AgNP)作为一种前途无量的电子芯片粘接材料,其宏观与微观的力学性能受到了广泛关注。Long等[113]使用穿透深度为2 000 nm的布氏压头研究了不同应变速率下烧结AgNP的微观结构与宏观本构行为之间的相关性,并基于纳米压痕响应的关键变量进行了无量纲分析,建立了烧结AgNP的宏观本构关系[114-116]。为了获得可靠的长期使用能力,烧结AgNP在各种载荷条件下的断裂行为受到了广泛关注。Chen等[117]发现烧结AgNP的断裂韧性随着试样尺寸的减小而降低。Kim等[118]研究了连接温度对烧结AgNP接头处微观结构与断裂行为的影响。在数值模拟方面,相场断裂模型对模拟烧结AgNP的裂纹萌生与断裂行为由重要意义。Zhu等[119]基于相场模型研究了具有各种孔隙率的烧结AgNP材料的裂纹扩展过程。Su等[120]提出了一种可靠的相场技术来准确描绘孔隙等微观结构并预测了烧结AgNP的断裂演化过程。

细观尺度的变形与宏观尺度的破坏机制间的关系仍有待考证。Long等[121]采用CPFEM模拟了烧结AgNP的晶体变形,提出了一个多晶有限元模型与Voronoi技术,以揭示烧结AgNP宏观与微观变形之间的基本关系,模型采用四节点平面应力单元CPS4离散,如图19所示。烧结AgNP在不同应变速率与温度下的应力应变曲线如图20所示。可以看出,模拟结果与实验结果吻合度较高,说明所建立的CPFE模型能够较好的预测出烧结AgNP多晶材料的应力-应变响应。

图19 具有100个晶粒的FE模型及其边界条件[121]
Fig.19 FE model with 100grains and its boundary conditions[121]

图20 不同温度与应变速率下的应力-应变曲线:(a)应变速率10-5s-1;(b)应变速率10-4s-1;(c)应变速率10-3s-1[121]
Fig.20 Stress-strain curves at different temperatures and strain rates:(a) strain rate of 10-5s-1;(b) strain rate of 10-4s-1;(c) strain rate of 10-3s-1[121]

Choi等[122]通过CPFEM,考虑晶体滑移与变形孪晶的影响,并使用孪晶重取向(predouinant twin reorien tation,PTR)模型,模拟了镁合金的织构演化与宏观性能。实验与模拟结果吻合度较高。Githens等[123]在此基础上,将实验与晶体塑性法相结合,研究了高强度稀镁合金的微观变形机制与宏观拉伸性能,并首次解决了镁合金中单个滑移系上的应变问题。

4.5 晶体塑性有限元模型在材料本征强度方面的应用

由于合金成分-显微组织-相成分之间极为复杂的耦合作用,使得目前人们在合金元素对合金的强化机制以及多相协调变形机理等方面的认知严重不足。基于纳米压痕的细观力学分析测试是研究各相本征强度、实现成分因素与组织因素解耦的有力手段。将纳米压痕实验与相应的CPFE模拟结合,能够有效排除晶体取向的影响,研究各相的本征性能。

Yin等[124-125]针对TiAl-Nb-Mo基TiAl合金的纳米压痕行为进行CPFE模拟,拟合出TiAl合金典型的α2、γ、βo三相的本征强度。Sun等[126]在此基础上,拟合出Nb、Mo含量不同的TiAl-Nb-Mo基合金γ、βo相的本征强度。图21为成分不同的βo/γ两相的纳米压痕载荷位移曲线和模型预测曲线的对比。作者将先前科研者的研究结论与本次实验模拟所得结论进行对比讨论,明确了Al、Nb、Mo元素对两相基体的强化效果。对于γ相基体,确定了Al和Nb对其有一定的强化作用,而Mo的影响较小甚至可能会导致基体软化。而由于βo相的存在必定需要Nb、Mo等β稳定元素的存在,因此难以排除Al元素的影响。对于βo基体而言,仅定性分析出Mo元素的强化作用显著小于Nb元素,βo相的本征强度可能会随着Al含量的降低而显著提升。对比Zhang等[127]拟合出的βo相强度及Chen等[128]拟合出的γ相强度,不难发现,两者结果相差十分明显。因此,合金成分导致的相成分的变化,对各相的本征强度影响极大,从而导致各相变形协调问题的出现,引起合金宏观力学性能的变化。

图21 βo/γ两相的纳米压痕载荷位移曲线和模型预测曲线的对比,其中黄底色为γ相,蓝底色为βo[126]
Fig.21 Comparison of the nanoindentation load-displacement curves of the βo/γ phase and the prediction curves of the model. The yellow background color represents the γ phase,and the blue background color represents the βo phase[126]

由于各相间的变形协调问题,使得难以从宏观力学性能上反推材料的本征性能,导致合金元素的固溶强化作用至今都无法定量描述。基于简单的Taylor模型,不难发现,材料的单轴拉伸强度大致为滑移系临界分切应力的3倍。但在上述研究中,针对TNM合金,室温下βo相{110}<111>滑移系的临界分切应力大致为500 MPa,γ相{111}<1-10]滑移系的临界分切应力大致为300MPa,α2相基面<a>滑移系临界分切应力大致为1 000MPa。纳米压痕深度均为150 nm,排除尺寸效应后,各相强度分别为200、120和400MPa。按大致50%的γ相计算,可以粗略估算出合金的强度应为630 MPa左右。实际上,这一估算强度远比不上合金的实际强度。这是因为,合金强度需要考虑晶界强化的作用。目前并没有完全准确的晶界强化晶体塑性本构模型,因此,微观本征强度反演宏观力学性能并不准确。将CPFEM与细观力学测试结合,从微观机理方向出发,研究材料的本征强度,仅能为合金成分设计提供参考,并不能作为理论指导与借鉴。

5 结语及展望

晶体塑性理论的应用及发展有助于更细致地了解材料塑性变形的过程与机理,从而更深层次调控材料的塑性加工过程,完善各类性能指标,满足各种情况与领域的应用。

近年来,大家对于晶体塑性理论的研究主要集中在两个层面:①是完善Taylor模型与FE模型的描述能力并扩展其应用领域,例如FE模型几乎不能预测变形过程中晶粒内部出现的变形带;②是分析FE模型中晶粒划分方式、晶粒形状、网格类型以及密度对模拟过程中产生的影响。

除了继续研究上述内容外,针对现有晶体塑性理论的局限性,需要发展更为丰富的本构理论及模拟技术以实现准确描述每一种塑性变形及塑性变形过程中每一个阶段的目标。这里主要对未来晶体塑性理论发展的3个方向进行展望。

(1)现有大部分晶体塑性本构模型都是针对于滑移这一塑性变形方式所设计,随之通过VPSC、EPSC等开源程序将孪生考虑进去,但除了滑移与孪生外,还有多种塑性变形方式,例如扭折、晶间滑移等,目前尚未有本构模型可以完整准确的描述这类塑性变形过程。

(2)建立多物理场耦合的本构理论,将相变、再结晶、析出、损伤这类热力学或冶金学物理场耦合进现有的本构理论中。发展新的建模与模拟技术,从而将宏观模拟方法与微观纳米尺度的模拟方法联系起来。

(3)晶体塑性有限元计算成本较高。将机器学习与晶体塑性深度耦合,形成“物理模型驱动+数据驱动”的新型计算方式。利用机器学习代替本构关系求解、RVE单元计算等昂贵部分,建立高效的代理模型。

参考文献:

[1]KHAN A S,HUANG S J. Continuum theory of plasticity[M]. New York:John Wiley & Sons,1995.

[2]ROTERS F,DIEHL M,SHANTHRAJ P,EISENLOHR P,REUBER C,WONG S L,MAITI T,EBRAHIMI A,HOCHRAINER T,FABRITIUS H O,NIKOLOV S,FRIÁK M,FUJITA N,GRILLI N,JANSSENS K G F,JIA N,KOK P J J,MA D,MEIER F,WERNERE,STRICKERM,WEYGANDD,RAABED. DAMASK-the düsseldorf advanced material simulation kit for modeling multiphysics crystal plasticity,thermal,and damage phenomena from the single crystal up to the component scale[J]. Computational Materials Science,2019,158:420-478.

[3]MARIN E B,DAWSON P R. On modelling the elasto-viscoplastic response of metals using polycrystal plasticity [J]. Computer Methods in Applied Mechanics and Engineering,1998,165(1):1-21.

[4]TAYLOR G I. The mechanism of plastic deformation of crystals.Part I. Theoretical[J]. Proceedings A,1934,145(855):362-387.

[5]SCHMID E. Yield point of a crystals:Critical shear stress law[A]//Proceedings of the 1st International Congress of Applied Mechanics[C]. Netherlands:Delft,1924.342.

[6]OROWAN E. Zur Kristallplastizität. III [J]. Zeitschrift Für Physik,1934,89:634-659.

[7]TAYLOR G I,ELAM C F. The distortion of an aluminium crystal during a tensile test[J]. Proceedings A,1923,102(719):643-667.

[8]HILL R,RICE J R. Constitutive analysis of elastic-plastic crystals at arbitrary strain [J]. Journal of the Mechanics and Physics of Solids,1972,20(6):401-413.

[9]NEEDLEMANA,TVERGAARDV. Comparison of crystalplasticity and isotropic hardening predictions for metal-matrix composites[J]. Journal of Applied Mechanics,1993,60(1):70-76.

[10]KOCKS U F,MECKING H. Physics and phenomenology of strainhardening:The FCC case[J]. Progress in Materials Science,2003,48(3):171-273.

[11]LOGAN R W,HOSFO R D W F. Upper-bound anisotropic yield locus calculations assuming<111>pencil glide [J]. International Journal of Mechanical Sciences,1980,22(7):419-430.

[12]胡赓祥,蔡珣,戎咏华.材料科学基础[M].上海:上海交通大学出版社,2006.

[13]胡烨之,李文峰,朱冬蓉,张雅琪. Von Mises屈服准则雅可比矩阵的显式表达及收敛性分析[J].江淮水利科技,2023,103(1):15-18.HU Y Z,LI W F,ZHU D R,ZHANG Y Q. Explicit expression and convergence analysis of Jacobian matrix of Von Mises yield criterion[J]. Jianghuai Water Resources Science and Technology,2023,103(1):15-18.

[14]YOO M H. Slip,twinning,and fracture in hexagonal close-packed metals[J]. Metallurgical Transactions A,1981,12(3):409-418.

[15]CHRISTIAN J W,MAHAJAN S. Deformation twinning[J]. Progress in Materials Science,1995,39(1-2):1-157.

[16]SALEM A A,KALIDINDI S R,DOHERTY R D. Strain hardening due to deformation twinning in α-titanium:mechanisms [J].Metallurgical and Materials Transactions A,2006,37(1):259-268.

[17]WASILEWSKI R J. Mechanism of Bcc twinning:Shear or shuffle?[J]. Matallurgical Transactions,1970,1(9):2641-2643.

[18]WANG J,YADAV S K,HIRTH J P,TOMÉ C N,BEYERLEIN I J.Pure-shuffle nucleation of deformation twins in hexagonal-closepacked metals[J]. Materials Research Letters,2013,1(3):126-132.

[19]WANG Y N,HUANG J C. Texture analysis in hexagonal materials[J]. Materials Chemistry and Physics,2003,81(1):11-26.

[20]SALEM A A,KALIDINDI S R,SEMIATIN S L. Strain hardening due to deformation twinningin α-titanium:Constitutive relationsand crystal-plasticitymodeling[J]. Acta Materialia,2005,53(12):3495-3502.

[21]AGNEW S R,BROWN D W,TOMÉ C N. Validating a polycrystal model for the elastoplastic response of magnesium alloy AZ31 using in situ neutron diffraction[J]. Acta Materialia,2006,54(18):4841-4852.

[22]LEBENSOHNR A,TOMÉ CN. A study of the stress state associated with twin nucleation and propagation in anisotropic materials[J].Philosophical Magazine A,1993,67(1):187-206.

[23]KUMAR M A,BEYERLEIN I J,LEBENSOHN R A,TOMÉ C N.Modeling the effect of neighboring grains on twin growth in HCP polycrystals[J]. Modelling and Simulation in Materials Science and Engineering,2017,25(6):064007.

[24]FENG B,BRONKHORST C A,ADDESSIO F L,MORROW B M,CERRETA E K,LOOKMAN T,LEBENSOHN R A,LOW T.Coupled elasticity,plastic slip,and twinning in single crystal titanium loaded by split-Hopkinson pressure bar [J]. Journal of the Mechanics and Physics of Solids,2018,119:274-297.

[25]章海明,徐帅,李倩,刘珈汝,尚晓晴,董湘怀,崔振山.晶体塑性理论及模拟研究进展[J].塑性工程学报,2020,27(5):12-32.ZHANG H M,XU S,LI Q,LIU J R,SHANG X Q,DONG X H,CUI Z S. Progress of crystal plasticity theory and simulations[J].Journal of Plasticity Engineering,2020,27(5):12-32.

[26]韩逢博. Ti-6Al-4V合金纳米压痕变形与高周疲劳行为CPFEM研究[D].西安:西北工业大学,2018.HAN F B. CPFEM study of nanoindentation and high cycle fatigue behavior for Ti-6Al-4V alloy[D]. Xi'an:Northwestern Polytechnical University,2018.

[27]ROTERS F,EISENLOHR P,HANTCHERLI L,TJAHJANTO D D,BIELER T R,RAABE D. Overview of constitutive laws,kinematics,homogenization and multiscale methods in crystal plasticity finite-element modeling:theory,experiments,applications[J]. Acta Materialia,2010;58(4):1152-1211.

[28]刘海军,方刚,曾攀.基于晶体塑性理论的大变形数值模拟技术[J].塑性工程学报,2006,13(2):1-8.LIU H J,FANG G,ZENG P. Numerical simulation technology based on crystal plasticity theory[J]. Journal of Plasticity Engineering,2006,13(2):1-8.

[29]YAGUCHI M,TAKAHASHI Y. Ratchetting of viscoplastic material with cyclic softening,part 2:Application of constitutive models [J]. International Journal of Plasticity,2005,21(4):835-860.

[30]DOQUET V. Twinning and multiaxial cyclic plasticity of a low stacking-fault-energy f.c.c. alloy[J]. Acta Metallurgica et Materialia,1993,41(8):2451-2459.

[31]LEBENSOHN R A,TOMÉ C N. A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals:application to zirconium alloys [J].Acta Metallurgica et Materialia,1993,41(9):2611-2624.

[32]甘元超.镁金属孪晶变形的实验和理论模型研究进[J].高压物理学报,2021,35(4):142-159.GAN Y C. Advances of experimental and theoretical models of magnesium twin deformation[J]. Chinese Journal of High Pressure Physics,2021,35(4):142-159.

[33]PEIRCE D,ASARO R J,NEEDLEMAN A. An analysis of nonuniform and localized deformation in ductile single crystals[J].Acta Metallurgica,1982,30(6):1087-1119.

[34]KALIDINDI S R,BRONKHORST CA,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.

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

[36]杨平.电子背散射衍射技术、几何晶体学与材料科学[J].电子显微学报,2008,27(6):425-431.YANG P. Electron backscatter diffraction technology,geometrical crystallography,and materials science[J]. Journal of Chinese Electron Microscopy Society,2008,27(6):425-431.

[37]TOMÉ C N,MAUDLIN P J,LEBENSOHN R A,KASCHNER G C. Mechanical response of zirconium-I. derivation of a polycrystal constitutive law and finite element analysis [J]. Acta Materialia,2001,49(15):3085-3096.

[38]MA A,ROTERS F. A constitutive model for fcc single crystals based on dislocation densities and its application to uniaxial compression of aluminium single crystals[J]. Acta Materialia,2004;52(12):3603-3612.

[39]ZHANG H M,DONG X H,WANG Q,LI H Z. Micro-bending of metallic crystalline foils by non-local dislocation density based crystalplasticity finite element model[J]. Transactions of Nonferrous Metals Society of China,2013,23(11):3362-3371.

[40]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.

[41]MA A,ROTERS F. A constitutive model for fcc single crystals based on dislocation densities and its application to uniaxial compression of aluminium single crystals [J]. Acta Materialia,2004,52(12):3603-3612.

[42]HUANG M,LI Z,TONG J. The influence of dislocation climb on the mechanical behavior of polycrystals and grain size effect at elevated temperature [J]. International Journal of Plasticity,2014,61:112-127.

[43]NGUYEN T,LUSCHER D J,WILKERSON J W. Adislocationbased crystal plasticity framework for dynamic ductile failure of single crystals[J]. Journal of the Mechanics and Physics of Solids,2017,108:1-29.

[44]PO G,HUANG Y,GHONIEM N. A continuum dislocation-based model of wedge microindentation of single crystals[J]. International Journal of Plasticity,2019,114:72-86.

[45]KAISER T,MENZEL A. A dislocation density tensor-based crystal plasticity framework[J]. Journal of the Mechanics and Physics of Solids,2019,131:276-302.

[46]SUDMANNS M,STRICKER M,WEYGAND D,HOCHRAINER T,SCHULZ K. Dislocation multiplication by cross-slip and glissile reaction in a dislocation based continuum formulation of crystal plasticity[J]. Journal of the Mechanics and Physics of Solids,2019,132:103695.

[47]ZECEVIC M,KNEZEVIC M. A dislocation density based elastoplastic self-consistent model for the prediction of cyclic deformation:Application to AA6022-T4 [J]. International Journal of Plasticity,2015,72:200-217.

[48]FOLLANSBEE P S,KOCKS U F. A constitutive description of the deformation of copper based on the use of the mechanical threshold stress as an internal state variable[J]. Acta Metallurgica,1988,36(1):81-93.

[49]MECKING H,KOCKS U F. Kinetics of flow and strain-hardening[J]. Acta Metallurgica,1981,29(11):1865-1875.

[50]章海明.位错密度晶体塑性模型及在微成形模拟中的应用[D].上海:上海交通大学,2013.ZHANG H M. Dislocation density crystal plasticity model and its application in micro-forming simulation [D]. Shanghai:Shanghai Jiao Tong University,2013.

[51]FRANCIOSI P. The concepts of latent hardening and strain hardening metallic single crystals[J]. Acta Metallurgica,1985,33(9):1601-1612.

[52]NICAISE N,BERBENNI S,WAGNER F,BERVEILLER M,LEMOINE X. Coupled effects of grain size distributions and crystallographic textures on the plastic behaviour of IF steels [J].International Journal of Plasticity,2011,27(2):232-249.

[53]KNEZEVIC M,ZECEVIC M,BEYERLEIN I J,BINGERT J F,MCCABE R J. Strain rate and temperature effects on the selection of primary and secondary slip and twinning systems in HCP Zr[J].Acta Materialia,2015,88:55-73.

[54]PEETERS B,KALIDINDI S R,HOUTTE P V,AERNOUDT E. A crystal plasticity based work-hardening/softening model for b.c.c.metals under changing strain paths[J]. Acta Materialia,2000,48(9):2123-2133.

[55]BERTIN N,TOMÉ C N,BEYERLEIN I J,BARNETT M R,CAPOLUNGO L. On the strength of dislocation interactions and their effect on latent hardening in pure magnesium[J]. International Journal of Plasticity,2014,62:72-92.

[56]GUO X R,MAO N D,KONG T Q,ZHANG J,SHEN J J,WANG C H,SUN C Y,LI P P,XIONG Z P. Strain rate-dependent plastic behavior of TWIP steel investigated by crystal plasticity model[J].Materials Science and Engineering:A,2024,891:145986.

[57]KALIDINDI S R. Incorporation of deformation twinning in crystal plasticity models [J]. Journal of the Mechanics and Physics of Solids,1998,46(2):267-290.

[58]ZAAFARANIN,RAABE D,SINGHRN,ROTERS F,ZAEFFERER S. Three dimensionalinvestigation of the texture and microstructure below a nanoindent in a Cu single crystal using 3D EBSD and crystal plasticity finite element simulations [J]. Acta Materialia,2006,54(1):1863-1876.

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

[60]ZAAFARANI N,RAABE D,ROTERS F,ZAEFFERER S. On the origin of deformation-induced rotation patterns below nanoindents[J]. Acta Materialia,2008,56:31-42.

[61]RAABE D,MA D,ROTERS F. Effects of initial orientation,sample geometry and friction on anisotropy and crystallographic orientation changesin single crystalmicrocompressiondeformation:A crystal plasticity finite element study[J]. Acta Materialia,2007,55:4567-4583.

[62]WEBER F,SCHESTAKOW I,ROTERS F,RAABE D. Texture evolution during bending of a single crystal copper nanowire studied by EBSD and crystal plasticity finite element simulations[J]. Advanced Engineering Materials,2008,10:737-741.

[63]HALL E O. The deformation and ageing of mild steel III:Discussion of results[J]. Procedings of the Physical Society. Section B,1951,64:747.

[64]PETCH N J. The cleavage strength of polycrystals[J]. Journal of Iron and Steel Research International,1953,174:25-28.

[65]EVERS L P,PARKS D M,BREKELMANS W A M,GEERS M G D. Crystalplasticitymodelwithenhancedhardeningbygeometrically necessary dislocation accumulation [J]. Journal of the Mechanics and Physics of Solids,2002,50(11):2403-2424.

[66]ASHBY M F. The deformation of plastically non-homogeneous materials[J]. The Philosophical Magazine,1970,21:399-424.

[67]ARSENLIS A,PARKS D M,BECKER R,BULATOV V V. On the evolution of crystallographic dislocation density in non-homog eneously deforming crystals [J]. Journal of the Mechanics and Physics of Solids,2004,52(6):1213-1246.

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

[69]GURTIN M E. On the plasticity of single crystals:Free energy,microforces,plastic-strain gradients [J]. Journal of the Mechanics and Physics of Solids,2000,48(5):989-1036.

[70]GURTIN M E. A gradient theory of single-crystal viscoplasticity that accounts for geometrically necessary dislocations[J]. Journal of the Mechanics and Physics of Solids,2002,50(1):5-32.

[71]GURTIN M E,ANAND L,LELE S P. Gradient single-crystal plasticity with free energy dependent on dislocation densities [J].Journal of the Mechanics and Physics of Solids,2007,55(9):1853-1878.

[72]EVERS L P. Strain gradient crystal plasticity based on dislocation densities[D]. Eindhoven:TechnischeUniversiteitEindhoven,2003.

[73]HAN C S,GAO H,HUANG Y. Mechanism-based strain gradient crystalplasticity-I. Theory[J]. Journal of the Mechanics and Physics of Solids,2005,53(5):1188-1203.

[74]HAN C S,GAO H J,HUANG Y G,NIX W D. Mechanismbased strain gradient crystal plasticity-II. Analysis [J]. Journal of the Mechanics and Physics of Solids,2005,53(5):1204-1222.

[75]NEMAT-NASSER S,NI L,OKINAKA T. A constitutive model for fcc crystals with application to polycrystalline OFHC copper [J].Mechanics of Materials,1998,30(4):325-341.

[76]BALL J M,JAMES R D. Fine phase mixtures as minimizers of energy[J]. Arch Ration Mech Anal,1987,100(1):13-52.

[77]YAGUCHI M,YAMAMOTO M,OGATA T. A viscoplastic constitutive model for nickel-base superalloy,part 1:Kinematic hardening rule of anisotropic dynamic recovery[J]. International Journal of Plasticity,2002,18(8):1083-1109.

[78]江河,董建新,张麦仓,姚志浩,杨静.服役条件下镍基高温合金应力松弛微观机制[J].金属学报,2019,55(9):1212-1220.JIANG H,DONG J X,ZHANG M C,YAO Z H,YANG J. Stress relaxation mechanism for typical nickel-based superalloys under service condition[J]. Acta Metallurgica Sinica,2019,55(9):1212-1220.

[79]LU P,JIN X C,LI P,SUN Y L,FAN X L. Crystal plasticity constitutive model and thermodynamics informed creep-fatigue life prediction model for Ni-based single crystal superalloy [J].International Journal of Fatigue,2023,176:107892.

[80]CHABOCHE J L. On some modifications of kinematic hardening to improve the description of ratchetting effects[J]. International Journal of Plasticity,1991,7:661-678.

[81]ZHANG S L,XUAN F Z. Interaction of cyclic softening and stress relaxation of 9-12% Cr steel under strain-controlled fatigue-creep condition:Experimental and modeling[J]. International Journal of Plasticity,2017,98:45-64.

[82]段园培. TB8合金热变形组织介观尺度研究与模拟[D].合肥:合肥工业大学,2009.DUAN Y P. Mesoscale study and simulation of hot deformation microstructure of TB8 alloy[D]. Hefei:Hefei University of Technology,2009.

[83]VAN HOUTTE P. A comprehensive mathematical formulation of an extended taylor-bishop-hill model featuring relaxed constraints,the Renouard-Wintenberger theory and a strain rate sensitivity model[J]. Textures and Microstructures,1988,8-9:313-350.

[84]VAN HOUTTE P,DELANNAY L,KALIDINDI S R. Comparison of two grain interaction models for polycrystal plasticity and deformation texture prediction[J]. International Journal of Plasticity,2002,18(3):359-377.

[85]HABRAKEN A M. Modelling the plastic anisotropy of metals[J].Archives of Computational Methods in Engineering,2004,11(1):3-96.

[86]HILL R. The elastic behaviour of a crystalline aggregate [J].Proceedings of the Physical Society. Section A,1952,65(5):349.

[87]MOLINARI A,CANOVA G R,AHZI S. A self consistent approach of the large deformation polycrystal viscoplasticity [J].Acta Metallurgica,1987,35(12):2983-2994.

[88]ERINOSHO T O,COLLINS D M,WILKINSON A J,TODD R I,DUNNE F P E. Assessment of X-ray diffraction and crystal plasticity lattice strain evolutions under biaxial loading [J].International Journal of Plasticity,2016,83:1-18.

[89]ZHANG H M,DIEHL M,ROTERS F,RAABE D. A virtual laboratory using high resolution crystal plasticity simulations to determine the initial yield surface for sheet metal forming operations[J]. International Journal of Plasticity,2016,80:111-138.

[90]JALILI M,SOLTANI B,PAPAZAFEIROPOULOS G,VU Q V.Micromechanical modeling of cyclic plastic deformation of ZK60 Mg alloy using 3D full-field crystal plasticity coupled with computational homogenization[J]. International Journal of Fatigue,2020,140:105794.

[91]LIU W C,CHEN B K,PANG Y,NAJAFZADEH A. A 3D phenomenological yield function with both in and out-of-plane mechanical anisotropy using full-field crystal plasticity spectral method for modelling sheet metal forming of strong textured aluminum alloy[J]. International Journal of Solids and Structures,2020,193:117-133.

[92]ZHAO Z,MAO W,ROTERS F,RAABE D. A texture optimization studyforminimumearinginaluminiumbyuseofatexturecomponentcrystal plasticity finite element method[J]. Acta Materialia,2004,52(4):1003-1012.

[93]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.

[94]KALIDINDI S R,DONOHUE B R,LI S Y. Modeling texture evolution in equal channel angular extrusion using crystal plasticity finite element models[J]. International Journal of Plasticity,2009,25(5):768-779.

[95]KUPKA D,HUBER N,LILLEODDEN E T. A combined experimental-numerical approach for elasto-plastic fracture of individual grain boundaries [J]. Journal of the Mechanics and Physics of Solids,2014,64:455-467.

[96]TIKHOVSKIY I,RAABE D,ROTERS F. Simulation of the deformation texture of a 17% Cr ferritic stainless steel using the texture component crystal plasticity finite element method considering texture gradients [J]. Scripta Materialia,2006,54(8):1537-1542.

[97]KALIDINDI S R,ANANG L. An approximate procedure for predicting the evolution of crystallographic texture in bulk deformation processing of FCC metals[J]. International Journal of Mechanical Sciences,1992,34(4):309-329.

[98]LI H W,YANG H. Prediction of texture evolution under varying deformation states through crystal plasticity finite element method[J]. Transactions of the Nonferrous Metals Society of China,2012,22:222-231.

[99]BALASUBRAMANIAN S,ANAND L. Plasticity of initially textured hexagonal polycrystals at high homologous temperatures:application to titanium[J]. Acta Materialia,2002,50(1):133-148.

[100]RIYAD I A,KNEZEVIC M. Field fluctuations viscoplastic selfconsistent crystal plasticity:Applications to predicting texture evolution during deformation and recrystallization of cubic polycrystalline metals[J]. Acta Materialia,2023,261:119395.

[101]AKHTAR S K,REHAN K,BABAK F. Multiaxial and non-proportional loading responses,anisotropy and modeling of Ti-6Al-4V titanium alloy over wide ranges of strain rates and temperatures[J].International Journal of Plasticity,2007,23:931-950.

[102]SUNG J H,KIM J H,WAGONER R H. A plastic constitutive equation incorporating strain,strain-rate,and temperature [J].International Journal of Plasticity,2010,26:1746-1771.

[103]VAN HOUTTE P,DELANNAY L,KALIDINDIB S R. Comparison of two grain interaction models for polycrystal plasticity and deformationtextureprediction[J]. International Journal of Plasticity,2002,18:359-377.

[104]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.

[105]LAN Y J,XIAO N M,LI D Z,LI Y Y. 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.

[106]XIAO N M,TONG M M,LAN Y J,LI D Z,LI Y Y. Coupled simulation of the influence of austenite deformation on the subsequent isothermal austenite-ferrite transformation [J]. Acta Materialia,2006,54(5):1265-1278.

[107]VONDROUS A,BIENGER P,SCHREIJÄG S,SELZER M,SCHNEIDER D,NESTLER B,HELM D,MÖNIG R. Combined crystal plasticity and phase-field method for recrystallization in a process chain of sheet metal production [J]. Computational Mechanics,2015,55(2):439-452.

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

[109]PARK J,ROUT M,MIN K M,CHEN S F,LEE M G. A fully coupled crystal plasticity-cellular automata model for predicting thermomechanical response with dynamic recrystallization in AISI 304LN stainless steel[J]. Mechanics of Materials,2022,167:104248.

[110]SHAH V,SEDIGHIANI K,DOKKUM J S V,BOS C,ROTERS F,DIEHL M. Coupling crystal plasticity and cellular automaton models to study meta-dynamic recrystallization during hot rolling at high strain rates[J]. Materials Science and Engineering:A,2022,849:143471.

[111]SAHOO S K,DHINWAL S S,VU V Q,TOTH L S. A new macroscopic strain hardening function based on microscale crystal plasticity and its application in polycrystal modeling [J]. Materials Science and Engineering:A,2021,823:141634.

[112]SAHOO S K,TOTH L S,MOLINARI A,LATYPOV M I,BOUAZIZ O. Plastic energy-based analytical approach to predict the mechanical response of two-phase materials with application to dual-phase steels[J]. European Journal of Mechanics/A Solids,2022,91:104414.

[113]LONG X,HU B,FENG Y H,CHANG C,LI M Y. Correlation of microstructure and constitutive behaviour of sintered silver particles via nanoindentation[J]. International Journal of Mechanical Sciences,2019,161:105020.

[114]LONG X,JIA Q P,SHEN Z Y,LIU M,GUAN C. Strain rate shift for constitutive behaviour of sintered silver nanoparticles under nanoindentation[J]. Mechanics of Materials,2021,158:103881.

[115]LONG X,JIA Q P,LI Z,WEN S X. Reverse analysis of constitutive properties of sintered silver particles from nanoindentations [J].International Journal of Solids and Structures,2020,191:351-362.

[116]LONG X,SHEN Z Y,JIA Q P,LI J,DONG R P,SU Y T,YANG X,ZHOU K. Determine the unique constitutive properties of elastoplastic materials from their plastic zone evolution under nanoindentation[J]. Mechanics of Materials,2022,175:104485.

[117]CHEN C T,NAGAO S,SUGANUMA K,JIU J T,SUGAHARA T,ZHANG H,IWASHIGE T,SUGIURA K,TSURUTA K. Macroscale and microscale fracture toughness of microporous sintered Ag for applications in power electronic devices[J]. Acta Materilia,2017,129:41-51.

[118]KIM M S,NISHIKAWA H. Effects of bonding temperature on microstructure,fracturebehavior andjoint strength of Agnanoporous bonding for high temperature die attach[J]. Materials Science and Engineering:A,2015,645:264-272.

[119]SU Y T,FU G C,LIU C Q,ZHANG K,ZHAO L G,LIU C Y,LIU A,SONG J N. Thermo-elasto-plastic phase-field modelling of mechanical behaviours of sintered nano-silver with randomly distributedmicro-pores[J]. Computer Methodsin Applied Mechanics and Engineering,2021,378:113729.

[120]SU Y T,FU G C,LIU C Q,LIU C Y,LONG X. Fatigue crack evolution and effect analysis of Ag sintering die-attachment in SiC power devices under power cycling based on phase-field simulation[J]. Microelectronics Reliability,2021,126:114244.

[121]LONG X,CHONG K,SU Y T,DU L M,ZHANG G Q.Connecting the macroscopic and mesoscopic properties of sintered silver nanoparticles by crystal plasticity finite element method[J].Engineering Fracture Mechanics,2023,281:109137.

[122]CHOI S H,KIM D H,LEE H W,SHIN E J. Simulation of texture evolution and macroscopic properties in Mg alloys using the crystal plasticity finite element method[J]. Materials Science and Engineering:A,2010,527:1151-1159.

[123]GITHENS A,GANESAN S,CHEN Z,ALLISON J,SUNDARA-RAGHAVAN V,DALY S. Characterizing microscale deformation mechanisms and macroscopic tensile properties of a high strength magnesium rare-earth alloy:A combined experimental and crystal plasticity approach[J]. Acta Materialia,2020,186:77-94.

[124]YIN B Q,XUE X Y,ZHANG M Q,DENG T Q,LI J S,TANG B.Parameter identification and pileup behavior of TiAl alloy through nanoindentation and crystal plasticity simulation [J]. Journal of Alloys and Compounds,2023,948:169743.

[125]YIN B Q,XUE X Y,TANG B,WANG W Y,KOU H C,LI J S.Experiments and crystal plasticity simulations for the deformation behavior of nanoindentation:Application to the α2 phase of TiAl alloy [J]. Materials Science and Engineering:A,2022,831:142283.

[126]孙进,程亮,王川云,朱彬,孙凌燕,杨光. TiAl-Nb-Mo基TiAl合金系的细观力学行为研究[J].塑性工程学报,2025,32(7):155-164.SUN J,CHENG L,WANG C Y,ZHU B,SUN L Y,YANG G.Study on mesomechanical behavior of TiAl-Nb-Mo-based TiAl alloy system [J]. Journal of Plasticity Engineering,2025,32(7):155-164.

[127]ZHANG M Q,LI J S,TANG B,KOU H C,FAN J K. Mechanical characterization and strain-rate sensitivity measurement of Ti-7333 alloy based on nanoindentation and crystal plasticity modeling[J].Progress in Natural Science:Materials International,2018,28:718-723.

[128]CHEN L,EDWARDS T E J,GIOACCHINO F D,CLEGG W J,DUNNE F P E,PHAM M S. Crystal plasticity analysis of deformation anisotropy of lamellar TiAl alloy:3D microstructurebased modelling and in-situ microcompression [J]. International Journal of Plasticity,2019,119:344-360.

Advances in Crystal Plasticity Theory and the Application of the Crystal Plasticity Finite Element Model

SUN Jin,ZHANG Tongxian,WANG Yinjie

(Kehua Holdings Co.,Ltd.,Changzhou 213001,China)

Abstract:Crystal plasticity,as an important theoretical approach for overcoming the limitations of macroscopic plastic homogenization,helps to link microscopic deformation with macroscopic properties and enables accurate prediction of material mechanical behavior. The kinematics of crystal plasticity deformation and the basis of crystal plastic deformation are reviewed in this paper. The advantages and disadvantages of the phenomenological crystal plasticity constitutive model and the mechanical crystal plasticity constitutive model were compared. According to the relationship between the single crystal plasticity model and the polycrystalline plastic model as well as their respective natures,the homogenization model is introduced,and the full-field crystal plastic model and the mean-field crystal plastic model are introduced. Representative examples of crystal plasticity finite element models in engineering applications are presented,and the breakthroughs brought by the application of this model are elaborated. The shortcomings of polycrystalline plasticity models and the limitations of the theory of crystal plasticity are proposed,and future development trends are proposed.

Key words:crystal plasticity;plastic deformation;constitutive models;homogenization model;finite element model

DOI:10.16410/j.issn1000-8365.2026.5182

中图分类号:TG146

文献标识码:A

文章编号:1000-8365(2026)01-0077-23

收稿日期:2025-08-20

作者简介:孙 进,2000年生,硕士生.主要从事金属材料铸造研发. Email:sunjin 131213@163.com

通信作者:张同先,1964年生.主要从事铸造合金、铸造工艺及设备研究. Email:570619090@qq.com

引用格式:孙进,张同先,王寅杰.晶体塑性理论进展及晶体塑性有限元模型的应用[J].铸造技术,2026,47(1):77-99.

SUN J,ZHANG T X,WANG Y J. Advances in crystal plasticity theory and the application of the crystal plasticity finite element model[J]. Foundry Technology,2026,47(1):77-99.

(责任编辑:李亚敏)