考虑晶体滑移面分解正应力的细观损伤模型1)

赵伯宇 胡伟平 孟庆春

∗(北京航空航天大学航空科学与工程学院,北京 100191)

†(中国工程物理研究院机械制造工艺研究所,四川绵阳 621900)

引言

材料在外载荷的作用下,其内部结构会发生破坏,生产微小的裂纹或孔洞等细观缺陷,这些细观缺陷会导致材料宏观力学性能的劣化,这种现象即称为损伤[1-2].经典的损伤力学是在连续介质力学的框架下,并考虑材料变形的热力学过程而建立起来的,故又被称为连续损伤力学,Sumio[3]对该理论做了详尽的研究.在工程应用领域中,将连续损伤力学理论与有限元方法相结合能够很好地分析构件的疲劳破坏与延性破坏等问题,其中国内学者詹志新等[4-7]在该方面做了深入地研究工作.

连续损伤力学采用宏观唯象的研究方法,研究宏观尺度上代表性体积单元(representative volume element,RVE) 的损伤与其力学性能的关系,没有深入到材料内部,从细观尺度上研究损伤机理.考虑到宏观研究尺度的局限性,研究者们提出了一系列的细观损伤研究方法,例如考虑非椭球夹杂的复合材料细观损伤模型[8]、基于相场理论的非局部宏-微观损伤模型[9]、固体损伤破坏的统一相场模型[10]等.

在众多的细观尺度研究方法中,晶体塑性理论将滑移系作为研究的基本对象,已经发展为一套成熟的细观研究方法[11-15].基于晶体塑性理论建立细观损伤模型也成为近年来研究的热点[16-19].晶体尺度的细观损伤模型的研究工作主要集中在3 个方面:损伤变量的定义、被劣化对象的选取以及损伤演化方程的建立.赵红彬[20]将损伤定义为承载位错滑移的滑移面劣化率,以标量损伤度来折减弹性模量与剪切应变率,使用剪切应变来控制损伤的演化过程.Boudifa 等[21]将晶体由宏观到细观划分为3 个尺度:RVE 尺度、晶粒尺度以及滑移系尺度,分别定义了各级尺度的损伤变量,并逐级耦合建立了各个尺度之间的损伤关系,实现了从细观到宏观的延性损伤计算.Hfaiedh 等[22]为每个滑移系定义了一个损伤标量,将每个滑移系的损伤标量进行线性加权叠加得到RVE 的损伤标量,实现对RVE 整体刚度的折减.Zghal 等[23]研究了晶体塑性理论的高周循环疲劳模型,将每个滑移系的损伤标量耦合为一个四阶损伤张量,并用于刻画晶体损伤的各向异性.

对于材料损伤的描述,除了使用连续损伤力学的思想外,还有其他的定义方法.Forest 等[24-26]在位移场中引入了一个新的自由度,该附加自由度被称为微损伤变量.他们通过严格的数学推导,建立了微损伤理论的力学模型,并将该模型与晶体塑性理论相结合,很好地模拟了宏观力学性能的退化.本文采用微损伤理论的思想方法,推导了有限变形框架下考虑晶体解理型损伤的本构方程、塑性流动方程以及损伤演化方程,建立了相应的数值计算方法,完成了材料的细观参数标定,研究了单轴拉伸模拟过程中的塑性流动与损伤演化规律,以期为多晶金属材料的延性损伤、疲劳损伤、蠕变损伤等问题提供新的分析思路与方法.

1 基于滑移系分析的细观损伤模型

为了描述晶体滑移系在受力和变形过程中产生的解理型微观裂纹对晶体后续变形的影响,本文在经典晶体塑性理论的基础上引入损伤变形梯度张量Fd,即认为晶体的解理型损伤会导致附加的变形,进而推导了考虑晶体损伤的变形运动学方程、本构方程、塑性流动方程以及损伤演化方程.

1.1 损伤耦合的变形运动学

如图1 所示,在单晶体均匀滑移模型[11]的基础上,引入损伤变形梯度张量Fd,进而将单晶体的损伤耦合变形运动过程用4 个构型来描述:未发生变形的参考构型B0、损伤构型B1、塑性构型B2、变形结束后的当前构型B.

图1 损伤耦合的变形梯度分解Fig.1 Deformation gradient decomposition considering damage

根据Shanthraja 等[27]提出的损伤耦合变形运动分解方法,将变形梯度张量F作如下分解

式中,Fd为损伤变形梯度张量,用于描述由参考构型B0到损伤构型B1的变形过程中因晶体解理而产生张开型裂纹的晶格缺陷效应;Fp为塑性变形梯度张量,用于描述由损伤构型B1到塑性构型B2的变形过程中晶体的塑性位错滑移;Fe为弹性变形梯度张量,用于描述由塑性构型B2到当前构型B 的变形过程中晶体的弹性伸长与刚性转动.

需要指出,本文模型所考虑的损伤类型主要为晶体解理型微观损伤,因此引入的损伤变形梯度张量Fd主要用于描述解理型微观裂纹的产生对晶体沿滑移面法线方向变形的影响,并以此来反映材料力学性能的劣化.同时,出于对数值实现的考虑,本文在式(1)中采用了损伤、塑性、弹性的变形运动学分解顺序.

文献[24,27] 采用面心立方晶体的滑移系晶格矢量来描述晶体内部的3 种微裂纹形式,并将损伤引起的附加变形定义为3 种微观裂纹引起的附加变形的叠加,即认为其中的解理断裂发生在滑移面上.对于本文所考虑的与晶体解理相关的张开型微观裂纹损伤模式,可以通过滑移面法向张开变形来描述其附加影响,于是将损伤速度梯度张量Ld定义为

需要指出,本文建立的损伤模型没有直接考虑另外两种微观裂纹形式对损伤附加变形的贡献,但在下文式(11)所定义的损伤自由能中将滑移面滑移的影响考虑进来,并通过其中的系数β 来反映晶体滑移对损伤能量耗散的影响程度.

根据速度梯度张量L可定义变形率张量D与自旋率张量W,并对其作如下分解

式中,De,Dp,Dd分别为弹性、塑性、损伤变形率张量,We,Wp,Wd分别为弹性、塑性、损伤自旋率张量.不难发现,tr(Dp)=0,这说明塑性应变为偏斜张量,其不会引起体积变化;tr(Dd)0,这表示由微观裂纹引起的体积变化;Wd=0,这说明损伤变形不会引起旋转效应.对比经典晶体塑性理论中的相关方程,损伤变形梯度张量Fd的引入没有改变弹性、塑性下的变形率张量和自旋率张量的数学形式,这也符合图1中对损伤变形梯度张量Fd的假设.

1.2 损伤耦合的本构方程

材料经历不可逆的变形过程,其内部结构会发生变化,进而产生塑性应变、微观裂纹等,因此损伤材料中Helmholtz 自由能密度函数ρΨ 的值取决于应变状态、损伤程度和位错硬化结构等[3].于是,对于考虑损伤效应的等温过程,Helmholtz 自由能密度函数ρΨ 可分为弹性、塑性、损伤3 个部分[26],即

式中,ρ 为材料密度,εe为弹性应变张量,γ 为累积剪切应变,d为累积损伤变量.

对于式(8) 所示的自由能密度函数表达形式需要进行说明.在经典的连续损伤力学中通常需要考虑损伤对弹性模量的折减作用.但是在细观损伤模型中,可以通过附加损伤应变来描述损伤的效果,当损伤应变显著大于弹性应变时可以忽略损伤对弹性模量的影响[24].在本文后续的计算结果也反映了tr(εd) ≫tr(εe) 的现象,因此本文忽略了损伤对弹性模量的影响.

损伤与塑性行为之间的耦合关系则较为复杂,众多学者提出各种耦合模型.Al-Rub 等[28-29]在建立自由能时独立考虑了损伤与塑性的贡献部分,并在推导的状态变量方程中体现了损伤与塑性的耦合关系.Balieu 等[30-31]在建立的损伤演化方程与塑性屈服方程中考虑了损伤与塑性对其的影响.本文采用了法向附加变形来考虑解理型微观裂纹的影响,因此假设损伤主要由法向分解正应力引起,而该法向应力与控制滑移的切向分解剪应力同时产生,于是认为损伤必然伴随着塑性行为,即没有微观塑性就没有微观损伤,故而在损伤自由能中引入了塑性滑移的影响.同时从变形运动学考虑,损伤附加变形的方向垂直于塑性滑移变形的方向,因此假设损伤并不影响滑移,进而在塑性自由能中不考虑损伤的影响.

弹性自由能密度函数ρΨe可表示为弹性应变张量εe的二次型,即

式中,C 为四阶弹性模量张量.

塑性自由能密度函数ρΨp由累积剪切应变γ 所决定,并采用二次型的数学形式给出[26],即

式中,R0为初始滑移阈值,h(γ) 为非线性硬化模量,kh,h0,hs分别为h(γ)的初始变化斜率、初始值、饱和值.

损伤自由能密度函数ρΨd由累积损伤变量d、累积剪切应变γ 所决定,同样地以二次型给出[26],即

式中,Y0为初始损伤阈值,β 为损伤与塑性之间的半耦合系数,H(d)为非线性软化模量,kH,H0,Hs分别为H(d)的初始变化斜率、初始值、饱和值.

对自由能密度函数ρΨ 求偏导数,可得到状态变量(εe,γ,d)相对应的热力学广义力(σ,R,Y),即

式中,R为临界分解剪应力,Y为损伤阈值,其中β 体现了损伤与塑性之间的半耦合关系,即R和Y均为关于γ 和d的函数.

在实际计算过程中采用Hill 等[12]基于塑性构型B2提出的Jaumann 率形式本构方程,并且考虑到临界分解剪应力R与损伤阈值Y均大于零,于是采用式(13)替换式(12),即

1.3 塑性流动方程与损伤演化方程

下面建立屈服准则与损伤准则.如图2 所示,分解剪应力τ 使滑移系产生面内滑移运动,分解正应力σn使滑移面产生法向张开运动.于是根据Schmid 定理,第α 个滑移系的分解剪应力τα、第β 个滑移面的分解正应力分别定义为

图2 分解剪应力与分解正应力Fig.2 Resolved shear stress and resolved normal stress

从能量耗散的角度来建立塑性流动方程和损伤演化方程.热力学第二定律揭示了不可逆热力学过程的熵增现象,于是对于等温绝热的有限变形过程,可采用Clausius-Duhem 不等式给出其定量描述[3],即

将式(8) 定义的自由能密度函数ρΨ(εe,γ,d) 对时间求导可得

2 数值计算方法

通用有限元分析软件ABAQUS 中允许用户通过UserMATerial(UMAT)子程序引入复杂的材料本构关系.本文借助UMAT 子程序中实现损伤耦合本构方程的数值计算过程,并主要实现两个功能:更新应力与状态变量、提供Jacobian 矩阵∂Δσ/∂Δε.

2.1 应力与状态变量的更新

式中,xn为增量步开始时刻tn的状态变量,xn+1为增量步结束时刻的状态变量,Δt为时间增量.

在式(28) 的迭代过程中,其他变量取增量步开始tn时刻的值.当迭代结束后,其他变量可根据其增量更新,且增量可由Δγα和Δdβ线性表出.考虑到更新形式是固定的,这里仅给出弹性应变与应力的更新式,即

2.2 Jacobian 矩阵∂Δσ/∂Δε 的计算

UMAT 子程序在完成应力与状态变量的更新后,还需提供Jacobian 矩阵∂Δσ/∂Δε,ABAQUS 利用每一个积分点的Jacobian 矩阵生成结构整体的刚度矩阵,用于判断结构整体是否达到平衡状态.

将损伤耦合的本构方程改写为增量形式可得

上式对Δε 求偏导并采用ABAQUS 中Voigt 记法可得

式中,C为二阶弹性模量张量,σ为一阶应力张量.进一步将上述列阵存储在矩阵中,即

式中,∂tr(Δεd)/∂Δε 等价于将∂Δd/∂Δε 沿列求和得到的向量.在变形过程中近似认为‖mβ‖ ≈1,于是有tr(Mβ)≈1.不难发现,{σ}⊗{∂tr(Δε)/∂Δε}为非对称矩阵,于是Jacobian 矩阵∂Δσ/∂Δε 也为非对称矩阵.

在式(33) 中∂Δγ/∂Δε 与∂Δd/∂Δε 仍是未知的,因此还需得到这两个矩阵的具体表达式.采用文献[34]中的线性内插法可得

式中,δij为Kronecker 符号,即i=j时δij=1,i≠j时δij=0,且计算过程中各个状态变量均采用增量步结束时的值.

综上,根据式(36) 计算得到矩阵∂Δγ/∂Δε 与∂Δd/∂Δε,然后将其代入式(33),便可得Jacobian 矩阵∂Δσ/∂Δε,在这个过程中需要求解6 个(N+M)阶线性方程组.

3 材料参数标定

在上述建立的细观损伤模型中包含若干材料细观参数,其中一些参数可以通过单晶材料的常规试验获得,但是大部分参数并不能通过试验直接获取,需要采用间接的方法进行标定.常用的一种方法是通过比对有限元模拟和试验得到的应力应变曲线来给出合适的参数值.于是,将待标定的材料参数作为设计变量,并建立一个目标函数用于定量地描述两条应力应变曲线之间的差距,那么材料参数的标定就变成了对于目标函数求极小值的最优化问题,这就实现了参数的逆向识别[35].

有限元模拟得到的应力应变曲线来自于单晶铜RVE 的单轴拉伸过程,试验得到的应力应变曲线来自于文献[36]中单晶铜[100]取向的单轴拉伸过程.为了保证模拟和试验中晶体取向的一致性,使有限元计算中晶体局部坐标系与全局坐标系重合,即两个坐标系不发生相对转动.

如图3 所示,有限元模拟中的RVE 为0.1 mm×0.1 mm×0.1 mm 的正方体,并沿其3 个边长方向等分为5 份,共计得到125 个C3D8 网格单元.位移载荷使整个加载端面沿着z轴正方向平移,且uz=0.02 mm,边界条件为

图3 z 正方向拉伸的代表性体积单元Fig.3 RVE subjected to z-axis tension

本文采用文献[34] 中的单晶铜材料,并给出细观损伤模型的基本材料参数,弹性常数C11,C12,C44分别为168.4 GPa,121.4 GPa,75.4 GPa,塑性指数np为10,损伤指数nd为10.

单晶弹性常数可以通过试验测得,塑性流动指数np与损伤演化指数nd作为幂指数被引入到式(24)、式(25)中,用于描述模型的率相关特性,且这两个参数在一定范围内影响较小.将这些量作为基本材料参数可以减少参数标定过程的计算量.

除基本材料参数外,将待标定的11 个材料参数列出,即

为了标定式(37) 中的11 个参数,采用最优化的思想,使用粒子群优化(particle swarm optimization,PSO)算法进行求解,并最终给出了标定结果,即

为了验证该标定方法所给出的材料参数的可靠性,将上式中的标定结果用于有限元计算中,图4 给出了有限元计算结果与试验结果[36]的对比.

图4 材料参数标定结果的验证Fig.4 Verification of material parameters calibration results

结果表明,根据标定结果给出的计算曲线与试验曲线一致,说明了标定方法的可行性和适用性.

4 RVE 单轴拉伸的计算分析

本节以单晶铜材料为例,其材料参数按式(38)取值,通过RVE 单轴拉伸有限元模拟,从细观尺度进行了塑性与损伤分析,以说明本文理论模型的合理性.

下面选取[001]加载取向,研究分析单晶铜RVE在单轴拉伸载荷下的塑性流动与损伤演化规律.有限元计算所采用图3 中的模型,位移载荷uz=0.06 mm.

如图5 所示,在单轴拉伸过程中应力应变曲线经历了线性、屈服、硬化、下降4 个阶段,这表明本文提出的细观损伤模型能够很好地描述金属材料的宏观单轴拉伸力学响应.

图5 考虑损伤的单轴拉伸应力应变曲线Fig.5 Uniaxial tensile stress-strain curves considering damage

对于单轴拉伸过程中的体积变化,可以用应变张量ε 的迹来描述,并且考虑到式(6)中加法分解形式,不难得到体积变化的分解,即

式中,塑性体积变化tr(εp)=0,这是由塑性变形率张量Dp的定义所决定的.于是,单轴拉伸过程中的体积变化可分为可恢复的弹性体积变化tr(εe) 与微裂纹引起的损伤体积变化tr(εd)两个部分,图6 给出了单轴拉伸过程中体积变化分解.

图6 考虑损伤的体积变化分解Fig.6 Volume fraction decomposition considering damage

在单轴拉伸过程中弹性体积变化tr(εe) 在屈服之前迅速增长,当材料屈服后弹性体积变化tr(εe)以较低的线性增长速率继续增大.随着载荷的进一步加大,在>0.04 后损伤体积变化tr(εd)快速增大,并逐渐占据体积变化的主导地位.同时,有限元计算结果也显示塑性体积变化tr(εp)始终为零.

如图7 所示,可以通过滑移系上的细观状态变量来分析塑性流动与损伤演化规律.

图7(a)给出了第1 个滑移系的分解剪应力τ(1)、临界分解剪应力R以及累积剪切应变γ 的变化过程.由图可知,分解剪应力τ(1)与临界分解剪应力R同步增大,两者的差值几乎保持不变,这就意味着滑移剪切应变率为定值,对于另外11 个滑移系的滑移剪切应变率也可以得出类似的结论,因此累积剪切应变γ 保持线性增长.

图7 滑移系上的细观状态变量Fig.7 Microscopic state variables on slip system

对比图7(a)和图7(b)两图可知,累积损伤变量d在>0.1 后迅速累积增大,考虑到式(13)定义的临界分解剪应力R,累积损伤变量d的增大会导致临界分解剪应力R的减小,因此临界分解剪应力R在

5 结论

基于晶体塑性理论的基本原理,针对晶体解理型损伤的问题,引入了损伤变形梯度张量Fd的概念,提出了考虑滑移面分解正应力的细观损伤模型,建立了相应的数值计算方法,完成了材料的细观参数标定,并采用所建模型进行了单晶RVE 单轴拉伸的模拟与分析.主要结论如下:

(1)本文提出的细观损伤模型中,塑性驱动力控制滑移面内剪切变形产生不可逆塑性变形,损伤驱动力控制晶体解理产生张开型裂纹,并考虑了塑性与损伤之间的半耦合作用关系,能够很好地描述考虑晶体解理型损伤的材料细观失效机理.

(2)基于粒子群优化算法,完成了细观损伤模型中材料参数的标定方法.该方法具有适应性好、无需梯度信息等优点,能够很好地进行多个材料参数的统一标定,并且标定效果良好.

(3)RVE 的单轴拉伸过程表明,本文模型能够很好地描述金属材料单轴拉伸过程的宏观力学响应,即表现为应力应变曲线的线性、屈服、硬化、下降4 个阶段.从细观尺度来看,损伤体积变化tr(εd)随着载荷的增大而逐渐占据主导地位,塑性累积是接近线性增长的,损伤累积近乎是指数上升的,且塑性累积的相对速度是低于损伤累积的相对速度.

本文研究工作主要建立了考虑晶体解理型微观裂纹损伤模式的细观损伤模型,包括理论框架、数值实现和参数标定,并初步验证了本文模型的合理性与有效性.在后续研究中,尚有大量工作需要开展,例如单晶体的原位拉伸试验与分析、多晶体的损伤模型、考虑滑开型与撕开型裂纹的晶体损伤模型等.