Abaqus内置黏弹性本构模型验证

一、Static,General分析步是否会激活Abaqus内置的黏弹性本构1、先说结论Static分析步里默认不会激活Prony时域黏弹性的时变松弛/蠕变行为只按瞬时弹性模量计算。就算材料库定义了 Viscoelastic(Prony)程序直接跳过黏弹性松弛演化只使用瞬时模量得到弹性解没有应力松弛、蠕变、滞后。官方文档原文Viscoelastic material response is always ignored in a static analysis.静力分析会忽略黏弹性响应采用瞬时弹性模量仅耦合场 / 准静态分析Visco才会求解 Prony 黏弹性的时间演化。2、验证算例单轴拉伸杆Prony 黏弹性做蠕变测试CaseA: *STATIC静力通用保持载荷持载 10s位移全程不变没有蠕变 → 黏弹性没激活CaseB*VISCO准静态同样载荷持载 10s位移随时间慢慢变大蠕变→ Prony 黏弹性正常激活。材料设定1 项 Prony简单好理解几何10mm 长杆截面 1mm²C3D8 单元左端全固定右端施加恒定拉力 500N应力 500MPa。分析步设置Static,General:场变量输出​载荷设置网格尺寸设置为1单元类型为C3D8完全积分。节点位移结果如下​更改分析步类型为Visco必须设置 CREEP TOLERANCE此处设为0.001控制黏弹性积分精度否则增量步容易爆:​重新计算节点位移结果​caseA水平线位移不变纯弹性Prony 休眠caseB下降曲线蠕变Prony 正常工作二、验证Abaqus内置黏弹性本构1、想要验证Abaqus内置的黏弹性本构模型有个问题需要先处理Abaqus的帮助文档提供的代码Abaqus - User Subroutines - Abaqus/Standard User Subroutines - UMAT - Example: Simple linear viscoelastic material是不区分分析步类型的即使是在 Static,General 分析步下它也会按照黏弹性本构来进行时间积分这个行为和Abaqus内置的黏弹性本构就不一致。所以在官方提供的UMAT代码里面加一个判断设置为只要分析步类型不是21黏性分析步对应21准静态显式时间积分就按线弹性计算。具体分析步类型对应的编号可以在帮助文档这一页找到Abaqus - Output - File Output Format - Results file output format - Procedure type keys​修改后的UMAT代码SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT, 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC) C INCLUDE ABA_PARAM.INC C CHARACTER*80 CMNAME DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS), 2 DDSDDT(NTENS),DRPLDE(NTENS), 3 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 4 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3), 5 JSTEP(4) DIMENSION DSDEL(NTENS,NTENS),DSDEC(NTENS,NTENS),STREND(NTENS), 1 PSTRAN(3),ASTRAN(3,3),barCp(3,3),Q(6,6), 2 DSTRES(6),D(3,3) C C Elambda PROPS(1) Emu PROPS(2) Elambda1 PROPS(3) Emu1 PROPS(4) Enu1 PROPS(5) C C determine step C IF(JSTEP(2).NE.21) Then ElamInst Elambda1/Enu1 EmuInst Emu1/Enu1 EBULK3 3*ElamInst2*EmuInst EG2 2*EmuInst EG EmuInst EG3 3*EmuInst ELAM ElamInst C DSDEL0.0 Do 10,K11,NDI Do 20,K21,NDI DSDEL(K1,K2)ELAM 20 continue DSDEL(K1,K1)ELAMEG2 10 continue c Do K1NDI1,NTENS DSDEL(K1,K1)EG ENDDO c Do K11,NTENS STREND(K1)STRAN(K1)DSTRAN(K1) ENDDO c DDSDDEDSDEL STRESSMATMUL(DDSDDE,STREND) c ELSE C C EVALUATE NEW STRESS TENSOR C EV 0. DEV 0. DO K11,NDI EV EV STRAN(K1) DEV DEV DSTRAN(K1) END DO C TERM1 .5*DTIME Enu1 TERM1I 1./TERM1 TERM2 (.5*DTIME*ElambdaElambda1)*TERM1I*DEV TERM3 (DTIME*Emu2.*Emu1)*TERM1I C DO K11,NDI DSTRES(K1) TERM2TERM3*DSTRAN(K1) 1 DTIME*TERM1I*(Elambda*EV 2 2.*Emu*STRAN(K1)-STRESS(K1)) STRESS(K1) STRESS(K1) DSTRES(K1) END DO C TERM2 (.5*DTIME*Emu Emu1)*TERM1I I1 NDI DO K11,NSHR I1 I11 DSTRES(I1) TERM2*DSTRAN(I1) 1 DTIME*TERM1I*(Emu*STRAN(I1)-STRESS(I1)) STRESS(I1) STRESS(I1)DSTRES(I1) END DO C C CREATE NEW JACOBIAN C TERM2 (DTIME*(.5*ElambdaEmu)Elambda1 1 2.*Emu1)*TERM1I TERM3 (.5*DTIME*ElambdaElambda1)*TERM1I DO K11,NTENS DO K21,NTENS DDSDDE(K2,K1) 0. END DO END DO C DO K11,NDI DDSDDE(K1,K1) TERM2 END DO C DO K12,NDI N2 K1-1 DO K21,N2 DDSDDE(K2,K1) TERM3 DDSDDE(K1,K2) TERM3 END DO END DO TERM2 (.5*DTIME*EmuEmu1)*TERM1I I1 NDI DO K11,NSHR I1 I11 DDSDDE(I1,I1) TERM2 END DO C C TOTAL CHANGE IN SPECIFIC ENERGY C TDE 0. DO K11,NTENS TDE TDE (STRESS(K1)-.5*DSTRES(K1))*DSTRAN(K1) END DO C C CHANGE IN SPECIFIC ELASTIC STRAIN ENERGY C TERM1 Elambda 2.*Emu DO K11,NDI D(K1,K1) TERM1 END DO DO K12,NDI N2 K1-1 DO K21,N2 D(K1,K2) Elambda D(K2,K1) Elambda END DO END DO DEE 0. DO K11,NDI TERM1 0. TERM2 0. DO K21,NDI TERM1 TERM1 D(K1,K2)*STRAN(K2) TERM2 TERM2 D(K1,K2)*DSTRAN(K2) END DO DEE DEE (TERM1.5*TERM2)*DSTRAN(K1) END DO I1 NDI DO K11,NSHR I1 I11 DEE DEE Emu*(STRAN(I1).5*DSTRAN(I1))*DSTRAN(I1) END DO SSE SSE DEE SCD SCD TDE - DEE ENDIF RETURN END2、验证案例取对应的Abaqus内置材料参数为计算方式见https://zhuanlan.zhihu.com/p/680522083?share_codeGQtpCH2wQA17utm_psn2093420198300852815注意第3部分建模时使用的η实际上是200相关参数根据链接中的计算方式重新计算即可​UMAT参数为3、Abaqus/CAE建模使用尺寸为1*1*1的单个单元模型作拉伸验证。先使用内置材料模型材料参数设置如下分析步设置单元底部完全固定在上表面施加0.1位移网格尺寸设置为1.5单元类型为C3D8。计算结果接下来使用UMAT子程序修改材料计算结果结果内置材料和UMAT材料模型结果一致验证通过。