晶体塑性计算及其在细观力学分析的应用

ABAQUS 软件 2003 年度用户论文集 - 1 - 晶体塑性计算及其在细观力学分析的应用* 张克实 通信地址，西北工业大学 350 信箱 邮编 710072 摘要：采用 ABAQUS/STANDARD 用户材料接口 UMAT 作材料本构关系分析，主要工作 有： （1）探讨单晶粘塑性有限变形的增量计算方法，以整体坐标系下的 Cauchy 应力为基本变 量，采用Newton-Raphson方法求解在晶体粘塑性滑移本构关系基础上建立的非线性方程组； （2）考虑材料细观弹性和塑性各向异性，对多晶材料 Voronoi 集合体细观粘塑性变形和细 观应力的分布和统计学特点进行了分析。 关键词：多晶体塑性；有限变形；各向异性；微观应力分布；Newton-Raphson 迭代 1 引言 多晶材料塑性分析的基础是单晶塑性模型和计算方法[1-3]。单晶体的塑性变形依赖于晶体滑移 面和滑移方向（与晶体内原子密排面和密排方向相联系）。晶体变形会造成晶格旋转，因而晶体 滑移面和滑移方向的取向发生变化，同时这种变形又造成材料微元的刚体旋转，使得材料变形主 轴连带转动。晶体塑性变形分析方法主要有两种：当前构形下次弹性形式率本构模型描述的方法 [2]和卸载构形下本构描述的方法[3,6,10,11]。前者可用显式积分过程，后者则多采用隐式积分过程。 本文工作主要包括两个部分。第一部分是介绍作者建议的一种单晶粘塑性有限变形的增量计 算方法：（1）该方法是基于率本构模型的方法；（2）以整体坐标系下的 Cauchy 应力为基本变量 建立晶体粘塑性滑移本构关系的非线性方程组；（3）时间积分是隐－显混合格式；（4）采用 Newton-Raphson 方法求解。本文的第二部分则介绍了：作者编制的 ABAQUS 材料用户程序用 于对多晶材料 Voronoi 集合体细观粘塑性变形和细观应力的分布和统计学特点进行的分析及结 果。 2 变形与应变 按乘法分解理论，弹性速度梯度张量和塑性速度梯度张量表示为 1−⋅= eee FFL & ， 11 −− ⋅⋅⋅= eppep FFFFL & (1) *本文内容部分曾在全国第六届计算力学会议交流 zhanglei Back ABAQUS 软件 2003 年度用户论文集 - 2 - 式中 eF 是总变形梯度， eF 代表晶体中晶格畸变和刚性转动所产生的变形梯度， pF 代表晶体沿 滑移方向的均匀剪切所对应的塑性变形梯度。于是对应于当前构形的速度梯度张量L为 pe LLFFL +=⋅= −1& (2) 塑性变形梯度的演化规律是[4,5] p n p FnmF ⋅   ⊗= Σ= )()()(1 )( αααα γ&& (3) 式中 )(αγ& 是第α滑移系的剪应变率， )(αm 和 )(αn 分别表示第α滑移系滑移方向的单位向量和滑移 面的单位法向量，n是晶体滑移系的数目。材料微元伸长率可表示为： ( ) )( 2 1 LLLD symT =+= (4) 于是可以有伸长率的和分解： pe DDD += , ( ) )( 2 1 eTeee sym LLLD =+= ， ( ) )( 2 1 pTppp sym LLLD =+= (5) 基于极分解，材料微元的增量平均旋转可以用增量刚体旋转 R∆ 表示，它从总的增量变形梯 度 F∆ 中求得： RVF ∆⋅∆=∆ (6) 式中 V∆ 是增量左伸长张量。参考 tt ∆+ 时刻构形定义增量中任意时刻的速度梯度 ( ) )()()()(ˆ tttttt TT ∆+∆⋅+∆⋅⋅+∆⋅∆+∆= RRLRRL ττ (7) 然后可定义对数应变增量[8] ∑ = ⊗∆=∆ 3 1 )ln( i iii nnε λ (8) 这里 in 是 tt ∆+ 时刻当前主伸长 iλ 的方向， iλ∆ 是主伸长增量，且 0>∆ iλ 。 利用(5)式，可以把应变增量张量作和分解 pe εεε ∆+∆=∆ (9) 注意到（3）式，则塑性速度张量可以写为 )(*)(*)( 1 1 )( ααα α γ&nmFLFL ⊗=⋅⋅= Σ= − nepep ， 1−⋅= ppp FFL & (10) 塑性伸长率可以写为 )(*)(*)(*)(*)( 1 )( 2 1 ααααα α γ&mnnmD ⊗+⊗= Σ= n p (11) 式中 )(*)( αα mFm ⋅= e ， 1)(*)( −⋅= eFnn αα (12) )(αm 、 )(αn 更新为 *)(αm 和 *)(αn 是为了考虑晶格的旋转和畸变对滑移的影响。记 )( 2 1 *)(*)(*)(*)(*)( ααααα mnnmP ⊗+⊗= , )( 2 1 *)(*)(*)(*)(*)( ααααα mnnmW ⊗−⊗= (13) 则 zhanglei Back ABAQUS 软件 2003 年度用户论文集 - 3 - ppp ΩDL += , )(*)( 1 αα α γ&PD Σ== n p , )(*)( 1 αα α γ&WΩ Σ== n p (14) 式中 pΩ 称为塑性旋率。 假设在一个从时间 t到时间 tt ∆+ 的时间增量中，总变形梯度的增量是 F∆ ，滑移系α 分解剪 应变增量是 )(αγ∆ 。本文用以下方法计算增量弹性变形梯度和增量晶格旋转： [ ] 121 )( −∆++⋅∆=∆ FFFL ttt , (15) 利用（6）式和（4）式可以给出增量弹性变形梯度的预报 etet n ete FFnmFLF ⋅   ⋅   ∆⊗⋅−∆=∆ − =Σ 1)()()( 1 )(' ααα α γ (16) 然后修正为 ( ) ( ) ( )'')(' 1)()()( 1 eeteet n eete FFFFnmFFLF ∆+⋅   ∆+⋅   ∆⊗⋅∆+−∆=∆ − =Σ αααα γ (17) 然后可以利用极分解计算增量晶格旋转 eR∆ 。 3 当前构形下的应力增量与晶体滑移粘塑性增量本构关系 考虑到有限变形情形在这一增量过程中有一增量刚体旋转，要满足应力描述的客观性，应力 的累积可以按照下式进行[8]： )( DσRσRσ ∆∆+∆⋅⋅∆=∆+ Tttt , (18) 式中σ 是 Cauchy 应力，它与对数应变是功共轭的。增量应力 σ∆ 由材料本构关系确定， D∆ 代 表该时间增量中材料微元的变形增量。 晶轴坐标系与整体固定坐标系标架相差一个旋转 crysT 。在当前时刻 crysettcrystcrystt TRTRT 0⋅=⋅∆= ∆+∆+ (19) 式中 crysT0 初始时刻坐标变换矩阵， ett R∆+ 是累积计算的晶格旋转矩阵。为简便起见，后面将 crystt T∆+ 简记为 crysT 。由于弹性应变很小，晶轴坐标系下对应于当前构形的 Cauchy 应力增量（近 似等于参照 t时刻构形的 Kirchhoff 应力增量）可以依照广义胡克定律求得： ecryscryscrys εCσ ∆=∆ >< : 4 (20) 式中 ><4 crysC 、 ecrysε∆ 分别是晶轴坐标系下的 4 阶各向异性弹性张量和弹性对数应变增量。于是 ( ) TcryscryseTcryscryscrysTttt TTεTCTRσRσ ⋅ ⋅∆⋅⋅+∆⋅⋅∆= ><∆+ :4 (21) 利用增量应变的弹塑性和分解即（9）式，（20）式可以写为 )(: 4 pcryscryscryscrys εεCσ ∆−∆=∆ >< (22) 式中 pcrysε∆ 是晶轴坐标系下的增量塑性应变，由整体坐标系下增量塑性应变通过变换得到 zhanglei Back ABAQUS 软件 2003 年度用户论文集 - 4 - )(*)( 1 )(*)( 1 αα α αα α γγ ∆=⋅   ∆⋅=∆ ΣΣ == crys n crys nTcryspcrys PTPTε (23) 式中 crysTcryscrys TPTP ⋅⋅= *)(*)( αα 。利用线性迭加计算当前增量步分解剪应变增量： tttt ∆γηγηγ∆ α∆αα ] )1[( )()()( && ++−= (24) 式中 10 ≤<η 是常量； )(αγ&t 和 )(αγ&tt ∆+ 分别是滑移系α 分解剪应变率的已知增量初始值和待求的增 量末了值。采用 Hutchinson 建议的粘塑性演化律[9] k g )( )( )( 0 )( )sgn( α α αα ττγγ && = (25) 0γ& 称为参考剪应变率，是待定材料常数；k是反映材料的率敏感性质的材料常数；而 )(αg 是与材 料单元的非弹性滑移的硬化函数，其演化规律为 )()( )()( β β αβ α γγγ && ∑= n hg 及 ∫∑= n dβ βγγ )( (26) 在当前增量步的增量变化可由下式确定 ∑ = ∆=∆ n hg 1 )()( )()( β β αβ α γγγ (27) 4 求解非线性增量应力应变的控制方程与迭代方法 建议按以下方法来求解单晶本构关系。首先将（24）式代入，（23）式右端括号内部分可以写 为：   ∆+  ∆−=∆ ∑∑Σ = ∆+ == n tt n t n tt 1 )(*)( 1 )(*)()(*)( 1 )1( α αα α αααα α γηγηγ && PPP (28) 注意到 Cauchy 应力与滑移系分解剪应力的关系： σP :*)()( αα =τ ，有 ∑Σ = ∆+ ∆+ ∆+∆+ =      = n k tt tt tttt n g τ 1 )( *)( )( 0 *)()(*)( 1 : )sgn( α α α αααα α γγ σPPP && (29) 因而可以得到关于 σtt ∆+ 的非线性方程组 )( σHrσ tttt ∆+∆+ −= (30) 式中 Tt RσRr ∆⋅⋅∆= Tcryscryst nTcryscryscrys t TTPεTCT ⋅        ⋅   ∆−−∆⋅⋅+ Σ= >< )(*)( 1 4 )1(: αα α γη & (31) Tcrysttcryscrystt t TσHCTσH ⋅  ⋅∆= ∆+ ><∆+ )(:)( 1 4 0γη & (32) crys n k tt tt ttTcrystt g τ T σP PTσH ⋅             ⋅= ∑ = ∆+ ∆+ ∆+∆+ 1 )( )*( )()*( 1 : )sgn()( α α α αα (33) （30）∼（33）式中未知量有 tt ∆+ 时刻的 σtt ∆+ 、 *)(αP 、 crysT 和 )(αgtt ∆+ 。本文对（30）∼（33） 式的求解采用迭代结合校正的方法来进行：（1）取 t时刻的 *)(αP 、 crysT 和 )(αg 代入（30）∼（33） 式可以采用 Newton-Raphson 迭代求得 σtt ∆+ ；（2）计算滑移系分解剪应变增量、弹性变形梯度 增量；（3）计算 tt ∆+ 时刻的 *)(αP 、 crysT 和 )(αgtt ∆+ ；（4）计算累积弹性变形梯度。 zhanglei Back ABAQUS 软件 2003 年度用户论文集 - 5 - 5 计算 Jacobi 矩阵的方法 计算 Jacobian 矩阵是为了更新有限元单元的切线刚度矩阵。用本文方法该矩阵可以在晶轴坐 标系下生成再转换至固定整体坐标系下，也可以直接在整体坐标系下生成（计算过程中需将所有 相关张量变换到固定坐标系下）。为节省篇幅，仅介绍在在晶轴坐标系下生成的 Jacobian 矩阵。 以下讨论中，应力和应变张量用向量σr 和εr来表示，四阶本构张量错误！不能通过编辑域代码创 建对象。也采用矩阵错误！不能通过编辑域代码创建对象。来表示。可以将本构关系用矩阵和向 量表示为[9]： Ωε ε σσ +∂ ∂= rr rr dd , (34) 式中 εσ rr ∂∂ / 构成的矩阵称为 Jacobian 矩阵, Ω是不依赖于 εrd 的向量。 将（25）式作对应于 t时刻的 Taylor 展开并略去二阶以上的高阶项，有 )( )( )( )( )( )( )()( β β α β β α αα γττ γγγ g g ttt ∆∂ ∂+∆∂ ∂+=∆+ &&&& , (35) 利用 (22)、(23) 和 Cauchy 应力与滑移系分解剪应力的关系： σP :*)()( αα =τ ，有    ∆−∆=−=∆ ∑ = ∆+ n cryscryscrysTcrysttt 1 )(*)(*)()()()( ˆˆ β ββαααα γτττ PεCP r (36) 即 )(ατ∆ 可以用 )(βγ∆ 来表示。而又由（27）式，知 )(αg∆ 亦可以用 )(βγ∆ 来表示。于是（36）式可 以改写为： γMrτ ∆−=∆ τ1 , cryscrysTcrys εCPr r∆≡ *1 ˆ , [ ] nncrysjcrysTcrysi ×≡ *)(*)( ˆˆ PCPMτ (37) (27)式也可以改写为 γMMγMg ∆=∆=∆ sgn|| hh , [ ]αβhh =M ， [ ]αββ δτ )sgn( )(sgn =M (38) 定义对角矩阵 τA 和 gA ，两者的非对角元素均为零。且 τA 和 gA 的对角元素分别是： )( )( i i tA τ γ ∂ ∂∆= &τii （不对 i求和）， )( )( i i g tA ∂ ∂∆= γ&gii （不对 i求和）， (39) 注意到（24）式，然后可以得到 [ ] ttg ∆+∆+∆=∆ γgAτAγ &τη (40) 将（37）和（38）两式代入（40）式，则可得 ( )[ ] ( )tthg ∆+−+=∆ − γrAMMAMAIγ &11sgn τττ ηη (41) 于是由（22）和（23）式，可以得到晶轴坐标系下 Jacobian 矩阵的近似表达 ( )[ ]{ }crysTcryshgcryscrysepcrys ηη CPAMMAMAIPICC *1sgn* ˆˆ τττ −−+−= (42) 可以用矩阵或四阶张量旋转变换的方法得到固定整体坐标系下的 Jacobian 矩阵 epC 。 6 算例：多晶 Voronoi 集合体应力应变的不均匀性分析 zhanglei Back ABAQUS 软件 2003 年度用户论文集 - 6 - 图 1 600 晶粒 Voronoi 多晶集合体 多数金属材料是多晶材料。由于晶粒构成和取向的随机性，宏观上看来均匀一致的材料细观 上实际是不均匀的，而且在小尺度范围必定是各向异性的。因此，研究这类材料微元的力学行为 应该考虑微元内部构成的随机性和晶粒力学行为的各向异性。 6.1 材料单晶晶粒的力学性质 作为初步分析，本文以文献[2]研究的铜多晶材料为研究对象，其单晶晶粒的材料弹性常数（立 方对称）为： GPaC 227.5111 = 、 GPaC 93.3612 = 、 GPaC 937.2244 = 。铜单晶具有 FCC 八面 体滑移系（见表 1），而粘塑性滑移模型参数为： MPa84.600 =τ 、 MPa51.109s =τ 、 MPah 5.5410 = 、 s/101 30 −×=γ& 、 1=q 表 1 FCC 八面体滑移系的滑移面法向矢量与滑移方向矢量 1 2 3 4 5 6 7 8 9 10 11 12 n )111( )111( )111( )111( )111( )111( )111( )111( )111( )111( )111( )111( m ]110[ ]101[ ]110[ ]101[ ]011[ ]101[ ]101[ ]011[ ]011[ ]110[ ]110[ ]011[ 6.2 多晶集合体几何、有限元划分与边界条件 取一个包含 600 个单晶晶粒的立方体为材料的代表性单元，晶粒是多面体，其几何按 Voronoi 体胞生成法则生成，且各晶粒的取向亦完全随机生成。作为初步分析，模型中每个晶粒均具有相 同的力学性质，不考虑晶界厚度和晶界的力学 特性（即几何 晶界），亦不考虑材料中的微孔洞和微裂纹。 用 Voronoi 方 法生成的多晶集合体[7]如图 1 所示。 多晶集合体的有限元单元划分按晶粒进 行，且保证相 邻晶粒在界面上节点重合。由于晶粒几何非常 复杂，有限元 模型采用四面体单元（C3D4）。模型总共包含 172923 个 单 元，32377 个节点[7]。 对该多晶集合体进行最简单的单轴拉伸模拟分析，加载通过在加载面施加均匀位移来进行， 法线方向与加载方向垂直的面用位移约束方程保证该面垂直位移一致（以满足周期性假设），并 保证没有宏观横向应力（宏观单轴拉伸应力状态）。 6.3 多晶集合体内应变与应力的不均匀性 zhanglei Back ABAQUS 软件 2003 年度用户论文集 - 7 - 多晶集合体内各晶粒由于晶粒弹性性质和粘塑性滑移变形的各向异性，又由于取向的随机 性，在宏观均匀单轴拉伸状态下实际多晶体内是处于复杂应力状态。每个晶粒取向不同，不同晶 粒的变形和应力不相同，每个晶粒内部变形和应力也不均匀。 图 2 显示出多晶集合体表面 Mises 等效应力、三轴拉伸应力和纵向拉伸应变的分布，从该图 可知多晶体中应力非常不均匀，而且 Mises 等效应力、三轴拉伸应力的极大、极小值均多出现在 晶粒交界处（它预示如果晶界不加以强化，材料破坏一般将从晶界开始）。图中也表明多晶体内 应变亦很不均匀，且有呈交叉带状局部剪切的倾向。 导致应力的不均匀的原因单晶晶粒的各向异性和单晶晶体的随机取向；不均匀的程度与单晶 的弹性各向异性、晶格结构、滑移系组成和多晶体集合方式都有关。在本文考虑的条件下，晶粒 和晶向完全随机安排，形成的多晶集合体近似于真实的等轴晶结构。对这样的计算结果可以进行 难以通过试验进行的应变、应力的统计分析，以探讨多晶材料统计意义下的细观力学行为。图 3 给出宏观拉伸应变为 0.039 时按体积平均的纵向应变、应力的分布。很明显，两者都符合高斯分 布的特点。 7 结论 图 2 多晶集合体内的应力应变分布：a)Mises 等效应力；b)三轴拉伸应力；c)纵向拉伸应变 a) b) c) a) b) 图 3 a)纵向应变分量分布；b)纵向应力分量分布 zhanglei Back ABAQUS 软件 2003 年度用户论文集 - 8 - 本文介绍了一种求解单晶粘塑性滑移本构关系的方法。结合相应的 Jacobian 矩阵 计算公式 六西格玛计算公式下载结构力学静力计算公式下载重复性计算公式下载六西格玛计算公式下载年假计算公式 ， 该方法可用于非线性有限变形有限元计算中的单晶本构计算。该方法可以采用 ABAQUS 用户材 料接口子程序（Standard/umat 和 Explicit/vumat 都可以）加以实现。 以下是本文研究的 总结 初级经济法重点总结下载党员个人总结TXt高中句型全总结.doc高中句型全总结.doc理论力学知识点总结pdf ： (1) 本文建议方法的特点是：以 Cauchy 应力为基本未知量，采用 Newton-Raphson 迭代和隐— 显式混合求解（总体是隐式计算但对临界分解剪应力是显示计算）过程。 (2) 单晶单轴拉伸颈缩模拟计算算例证实了方法对晶体粘塑性滑移有限变形计算的有效性。 (3) 可以应用本文方法结合随机方法生成的 Voronoi 多晶集合体模型进行多晶力学行为分析。 (4) 多晶材料中应力应变分布由于晶粒取向的随机性而极不均匀，Mises 等效应力和三轴拉伸应 力的极大、极小值都出现在（三叉）晶界处。 (5) 多晶集合体内应变有呈带状分布倾向。 (6) 多晶集合体中拉伸应力和拉伸应变分量统计意义上可以认为符合高斯分布。 致谢 本文 Voronoi 多晶集合体及有限元网格采用了文献[7]的数据，得到了 Wu MS 教授 （Nanyang Technological University）、Feng R 教授（University of Nebraska-Lincoln, USA）和薛生强硕士的帮助，在此谨表衷心感谢！ 参 考 文 献 [1] Taylor G.L., plasticstrain in metals, J. Inst. Metals, 62 (1938) 307-324. [2] Asaro, R.J., Needleman, A., Texture development and strain hardening in rate dependent polycrystals, Acta Metall. 33 (1985) 923-953. [3] Gorti Sarma, Thomas Zacharia. Integration algorithm for modeling the elastoviscoplastic response of polycrystalline materials. Journal of the Mechanics and Physics of Solids, 47 (1999) 1219-1238. [4] R. Hill and J.R. Rice, Constitutive analysis of elastic-plastic crystal at arbitrary strain, J. Mech. Phys. Solids, 20 (1972) 401-413. [5] Asaro R. J. and Rice J.R., Strain localization in ductile single crystals, J. Mech. Phys. Solids, 25 (1977) 309-338. [6] Surya R. Kalidindi, Modeling anisotropic strain hardening and deformation textures in low stacking fault energy fcc metals, Int. J. Plasticity, 17 (2001) 837-860. zhanglei Back ABAQUS 软件 2003 年度用户论文集 - 9 - [7] Zhang K.S., Wu M.S. and Feng R, Simulation of microplasticity-induced deformation in shocked ceramics by 3-D Voronoi polycrystal modeling, Int. J. Plasticity (待发表). [8] ABAQUS Reference Manuals, V6.2, 2001. Hibbit, Karlsson & Sorenson. [9] Hutchinson, J.W., 1976. Bounds and self-consistent estimates for creep of polycrystalline materials. Proc. R. Soc. Lond. 348(A), 101-127. [10] 张克实：变温情形下单晶体热－粘塑性有限变形的分析方法，机械科学与技术（2001 年增刊）， 第六届全国细观力学会议论文集，2001 年 11 月（西安），第 1-4 页。 [11] 冯露，张克实、张光：延性单晶非均匀变形的三维有限元模拟，应用力学学报，2002 年，第 19 卷，第 4 期，第 5-9 页。 zhanglei Back