·1200· 工程科学学报，第39卷，第8期 程，因此需要在屈服准则与弹性预测阶段考虑到损伤 式中：h为塑性硬化模量，h=mK(E。+E)-1;△y为 的影响.与损伤相耦合的米塞斯(Von Mises)屈服准 塑性流动因子；右下标t与t+△：分别表示t与1+△1 则表达式为： 时刻的变量值. fo,D,0,)=G-(1-D)g,=0 (8) 此外，根据塑性应变增量△ε”与试探应力σ"，可 式中：0，为流动应力：G为等效应力，G=√层S5:S 求得t+△时刻的偏应力S,: S+x=S-2G(1-D,)△e°= 为偏应力，S=0-oL.这里，0表示应力，e表示应 S"-2G(1-D.)△yn. (17) 变，1表示单位张量. 联立式(12)与式(17)，可得到关于△y的一元二 与损伤相耦合的广义胡克定律为： 次方程： △o=(1-D)(2GAe+Atrace(△e)I).(9) A△y2+2B△y+C=0. (18) 式中：△e为弹性应变增量；G为剪切模量，G= 式中： (19) E 21M为拉梅常数，A=1+0-2 为将Lemaitre损伤模型嵌入到Abaqus有限元软 B=o-61-n)+是2）， (20) 件中，采用向前欧拉数值积分算法编写VUMAT子程 序，主要通过弹性预测与塑性修正两阶段来更新应力 C=s-√/层1-0) (21) 状态，其中，塑性修正阶段使用了径向回退算法，具体 如下o) 4”号+w+31-2受广门广， 首先假定t时刻的应变增量△ε全部为弹性应变， (22) 即△e=0,根据广义胡克定律，试探应力σ"为： h,=mK(E。+E)-. (23) o"=o+(1-D,)(2GAE+Atrace(AE)I),(10) 根据式(18)可求得△y: $=o-e(). (11) △y=-B±VB-AC (24) A 式中：右下标1表示t时刻的变量值. △y=min(△y:),△y:>0,(i=1,2). (25) 将式(11)得到的试探偏应力S代入到屈服准则 最后，更新t+△1时刻下的各个状态变量： 式(8)中，假如∫≤0，则说明弹性变形假设是正确的. o+=-2G1-D,)△E=o-2G(1-D,)△yn, 然后，更新t+△：时刻的各个状态变量. (26) 假如∫>0，则说明屈服面将要向外扩张，并满足一 致性条件，可表示为： +层4， Ey=8+AE”=8:+ (27) S=R. (12) (28) 式中：R,为t+△时刻的屈服面半径；n为屈服面外 D=D.+△D=D,+√，y, 法线方向上的单位张量： ey=e+△E-△e”, (29) ew=e+△e. (30) (13) 为了验证所编写子程序的可靠性，采用文献[10] 将式(12)代入到屈服准则式(8)中，可得到： 提供的模型参数与边界条件，建立轴对称无模冲孔二 R=1-s加y 维有限元模型，压强p从0MPa开始以线性增长方式 (14) 加载直至冲孔结束.图1给出了压强p分别为12.16 采用Swif形式的硬化法则，屈服应力为o,= MPa与12.32MPa状态时的损伤分布以及最终断面 K(。+8)”,式中：m为应变硬化指数，K为材料常数. 图，对比文献结果（见图2），损伤分布、断后形貌与试 其中： 验结果吻合良好，说明子程序是准确的，可用于塑性加 工有限元模拟中. 1.3 Lemaitre连续损伤模型参数的确定 为了定量衡量塑性变形过程中损伤演化，采用试 (15) 验与数值模拟相结合的方法并根据反复加载-卸载拉 0+=0+h,△E= 伸试验中杨氏模量变化来确定Lemaitre模型参数&o、 +k(+)层 e与D,如图3所示，具体步骤如下) (16) (1)首先，在万能材料试验机上进行反复加载-卸工程科学学报,第 39 卷,第 8 期 程,因此需要在屈服准则与弹性预测阶段考虑到损伤 的影响. 与损伤相耦合的米塞斯(Von Mises) 屈服准 则表达式为: f(滓,D,滓y) = 滓 - (1 - D)滓y = 0. (8) 式中:滓y 为流动应力;滓 为等效应力,滓 = 3 2 S颐 S;S 为偏应力,S = 滓 - 滓H I. 这里,滓 表示应力,着 表示应 变,I 表示单位张量. 与损伤相耦合的广义胡克定律为: 驻滓 = (1 - D)(2G驻着 e + 姿trace(驻着 e )I). (9) 式中: 驻着 e 为 弹 性 应 变 增 量; G 为 剪 切 模 量, G = E 2(1 + 自) ;姿 为拉梅常数,姿 = 自E (1 + 自)(1 - 2自) . 为将 Lemaitre 损伤模型嵌入到 Abaqus 有限元软 件中,采用向前欧拉数值积分算法编写 VUMAT 子程 序,主要通过弹性预测与塑性修正两阶段来更新应力 状态,其中,塑性修正阶段使用了径向回退算法,具体 如下[10] . 首先假定 t 时刻的应变增量 驻着 全部为弹性应变, 即 驻着 p = 0,根据广义胡克定律,试探应力 滓 tr为: 滓 tr = 滓t + (1 - Dt)(2G驻着 + 姿trace(驻着)I), (10) S tr = 滓 tr - 1 3 trace(滓 tr )I. (11) 式中:右下标 t 表示 t 时刻的变量值. 将式(11)得到的试探偏应力 S tr代入到屈服准则 式(8)中,假如 f臆0,则说明弹性变形假设是正确的. 然后,更新 t + 驻t 时刻的各个状态变量. 假如 f > 0,则说明屈服面将要向外扩张,并满足一 致性条件,可表示为: St + 驻t = Rt + 驻tn. (12) 式中:Rt + 驻t为 t + 驻t 时刻的屈服面半径;n 为屈服面外 法线方向上的单位张量: n = St + 驻t | St + 驻t | = S tr | S tr | . (13) 将式(12)代入到屈服准则式(8)中,可得到: Rt + 驻t = 2 3 (1 - Dt + 驻t)滓y,t + 驻t . (14) 采用 Swift 形式的硬化法则,屈服应力为 滓y = K(着0 + 着 p ) m ,式中:m 为应变硬化指数,K 为材料常数. 其中: Dt + 驻t = Dt + Dc 着R - 着 ( D 滓y,t ) K [ 2 2 3 (1 + 自) + 3(1 - 2自) ( 滓H,t 滓 ) t ] 2 2 3 驻酌, (15) 滓y,t + 驻t = 滓y,t + ht驻着 p = 滓y,t + mK (着0 + 着 p t ) m - 1 2 3 驻酌. (16) 式中:h 为塑性硬化模量,h = mK (着0 + 着 p ) m - 1 ;驻酌 为 塑性流动因子;右下标 t 与 t + 驻t 分别表示 t 与 t + 驻t 时刻的变量值. 此外,根据塑性应变增量 驻着 p 与试探应力 滓 tr ,可 求得 t + 驻t 时刻的偏应力 St + 驻t: St + 驻t = S tr - 2G(1 - Dt)驻着 p = S tr - 2G(1 - Dt)驻酌n. (17) 联立式(12)与式(17),可得到关于 驻酌 的一元二 次方程: A驻酌 2 + 2B驻酌 + C = 0. (18) 式中: A = 2 3 2 3 琢tht, (19) B = 1 3 琢t滓y,t - G(1 - Dt) ( 1 + ht 3 ) G , (20) C = S tr 颐 S tr - 2 3 (1 - Dt)滓y,t, (21) 琢t = Dc 着R - 着 [ D 2 3 (1 + 自) +3(1 -2自) ( 滓H,t 滓 ) t ] ( 2 滓y,t ) K 2 , (22) ht = mK (着0 + 着 p t ) m - 1 . (23) 根据式(18)可求得 驻酌: 驻酌 = - B 依 B 2 - AC A , (24) 驻酌 = min (驻酌i),驻酌i > 0, (i = 1,2). (25) 最后,更新 t + 驻t 时刻下的各个状态变量: 滓t + 驻t = 滓 tr - 2G(1 - Dt)驻着 p = 滓 tr - 2G(1 - Dt)驻酌n, (26) 着 p t + 驻t = 着 p t + 驻着 p = 着 p t + 2 3 驻酌, (27) Dt + 驻t = Dt + 驻D = Dt + 2 3 琢t驻酌, (28) 着 e t + 驻t = 着 e t + 驻着 - 驻着 p , (29) 着 p t + 驻t = 着 p t + 驻着 p . (30) 为了验证所编写子程序的可靠性,采用文献[10] 提供的模型参数与边界条件,建立轴对称无模冲孔二 维有限元模型,压强 p 从 0 MPa 开始以线性增长方式 加载直至冲孔结束. 图 1 给出了压强 p 分别为 12郾 16 MPa 与 12郾 32 MPa 状态时的损伤分布以及最终断面 图,对比文献结果(见图 2),损伤分布、断后形貌与试 验结果吻合良好,说明子程序是准确的,可用于塑性加 工有限元模拟中. 1郾 3 Lemaitre 连续损伤模型参数的确定 为了定量衡量塑性变形过程中损伤演化,采用试 验与数值模拟相结合的方法并根据反复加载鄄鄄 卸载拉 伸试验中杨氏模量变化来确定 Lemaitre 模型参数 着D 、 着R 与 Dc,如图 3 所示,具体步骤如下[13] . (1)首先,在万能材料试验机上进行反复加载鄄鄄卸 ·1200·