UMAT实现Mises本构模型的详细教程

一.  ABAQUS自带mises本构与UMAT对比

mises-path.jpg

无标题.jpg

mises本构模型UMAT(附源代码和详细注释)的图3


二.  如何在分析中使用自定义材料


1.通过ABAQUS/CAE使用自定义材料

在Property模块中,执行【Material】/【Create】命令(或单击相关工具箱区的按钮),弹出Edit Material对话框(如图所示),用户可以通过该对话框选择材料模型、设置材料参数。对于自定义模型,执行对话框的【General】/【User Material】命令,此时在Material Behaviors区域中会出现User Material字样,表明定义的是用户材料。

在Edit Material对话框下方的User Material 区域中的【User material type】下拉列表中有三个选项,分别为Mechanical(力学材料)、 Thermal (热学材料)和Thermomechanical (热力学模型),默认选项为Mechanical。对于一些岩土材料模型,尤其是采用非关联流动法则的模型,雅克比矩阵是不对称的,此时需勾选Use unsymmetric material stiffiiess matrix 复选框。材料的参数在 Data 区域中的【Mechanical Constants】列表中输入,这里的数据会按次序传给UMAT子程序中的PROPS (NPROPS)数组,数据的个数即为NPROPS。

material.jpg

mises本构模型UMAT(附源代码和详细注释)的图5

注:输入材料参数时,在数据列表尾部按回车健会自动增加一行数据.此外,在数据窗口上右击会显示弹出式菜单,可实现数据的拷贝、增加、删除等功能。

如果UMAT子程序用到状态变量,还需设置状态变量的个数。具体操作仍然在Edit Material对话框中进行,执行对话框的【General】/【Depvar】命令,在【Number of solution-dependent state variables】输入框中设罝状态变量的个数。



2.通过inp输入文件使用自定义材料

在inp输入文件中的材料定义选项块(* Material)中使用下列语句即可:

*User material,type=Mechanical,constants= number of constants,Unsymm

关键字User material表明定义的用户材料;type为材料类型,可为Mechanical(力学材料)、 Thermal (热学材料)和Thermomechanical (热力学模型),Mechanical是默认选项:Constants 选项给出的是材料参数的个数;Unsymm选项只有当雅克比矩阵非对称时才采用。该关键字行下的数据行中的数据为材料参数。

状态变量的个数同样在材料定义选项块中定义,相应的关键字行为:

*Depvar

Number of solution-dependent state variables:状态变量个数

Inp.jpg




三.  Mises弹塑性模型

mises本构模型UMAT(附源代码和详细注释)的图8

3.jpg

4.jpg

5.jpg

6.jpg

7.jpg

mises本构模型UMAT(附源代码和详细注释)的图14

4.源程序和注释

需要启用沙漏控制:

mesh沙漏.jpg

     SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,

    1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,

    2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,

    3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,

    4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)

C

     INCLUDE 'ABA_PARAM.INC'

C

     CHARACTER*80 MATERL

     DIMENSION STRESS(NTENS),STATEV(NSTATV),

    1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),

    2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),

    3 PROPS(NPROPS),COORDS(3),DROT(3,3),

    4 DFGRD0(3,3),DFGRD1(3,3)

C    

     DIMENSION EELAS(6),EPLAS(6),FLOW(6)

C     EELAS(1-6)为弹性应变分量,EPLAS(1-6)为塑性应变分量,FLOW(6)为屈服面对应力的导数

     PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)

     DATA NEWTON,TOLER/10,1.D-6/

C

C -----------------------------------------------------------

C     UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC PLASTICITY

C     J2 FLOW THEORY

C     CAN NOT BE USED FOR PLANE STRESS

C -----------------------------------------------------------

C     PROPS(1) - E ,为弹性模量

C     PROPS(2) - NU ,为泊松比

C     PROPS(3) - SYIELD, 为屈服应力

C     CALLS AHARD FOR CURVE OF SYIELD VS. PEEQ

C -----------------------------------------------------------

C 程序不能用于平面应变

     IF (NDI.NE.3) THEN

        WRITE(6,1)

1       FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',

    1          'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')

     ENDIF

C

C     ELASTIC PROPERTIES

C

     EMOD=PROPS(1)

     ENU=PROPS(2)

     IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499

     EBULK3=EMOD/(ONE-TWO*ENU)

     EG2=EMOD/(ONE+ENU)

     EG=EG2/TWO

     EG3=THREE*EG

     ELAM=(EBULK3-EG2)/THREE

C     得到拉梅常数λ(ELAM),剪切模量G(EG)

C     ELASTIC STIFFNESS

C

     DO 20 K1=1,NTENS

       DO 10 K2=1,NTENS

          DDSDDE(K2,K1)=0.0

10     CONTINUE

20   CONTINUE

C

     DO 40 K1=1,NDI

       DO 30 K2=1,NDI

          DDSDDE(K2,K1)=ELAM

30     CONTINUE

       DDSDDE(K1,K1)=EG2+ELAM

40   CONTINUE

     DO 50 K1=NDI+1,NTENS

       DDSDDE(K1,K1)=EG

50   CONTINUE

C    得到(6-1)中的弹性模量矩阵

C    CALCULATE STRESS FROM ELASTIC STRAINS

C    假设应变增量全为弹性,计算应力增量(即第一次应力预测)

     DO 70 K1=1,NTENS

       DO 60 K2=1,NTENS

          STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)

60     CONTINUE

70   CONTINUE

C

C    RECOVER ELASTIC AND PLASTIC STRAINS

C

     DO 80 K1=1,NTENS

        EELAS(K1)=STATEV(K1)+DSTRAN(K1)

C     STATEV(1-6)存储弹性应变分量,此时假定应变增量全为弹性    

        EPLAS(K1)=STATEV(K1+NTENS)

C     STATEV(6-12)存储塑性应变分量  

80   CONTINUE

     EQPLAS=STATEV(1+2*NTENS)

C     STATEV(13)为等效塑性应变  

C     如果没有给出屈服应力,材料认为是弹性的

C

     IF(NPROPS.GT.2.AND.PROPS(3).GT.0.0) THEN

C

C       MISES STRESS

C

       SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2)) +

    1          (STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3)) +

    1          (STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1))

       DO 90 K1=NDI+1,NTENS

             SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1)

90     CONTINUE

       SMISES=SQRT(SMISES/TWO)

C       计算mises应力

C       HARDENING CURVE, GET YIELD STRESS

C       等效塑性应变和硬化曲线确定硬化后的屈服应力,初始等效塑性应变为0,

C       对应的屈服应力为初始屈服应力

       NVALUE=NPROPS/2-1

       CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)

C

C       DETERMINE IF ACTIVELY YIELDING

C

       IF (SMISES.GT.(1.0+TOLER)*SYIEL0) THEN

C       如果Mises应力超出了屈服应力,判断发生了屈服,进入了塑性阶段

C         FLOW DIRECTION

C

         SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE

C       SHYDRO为平均应力

         ONESY=ONE/SMISES

         DO 110 K1=1,NDI

            FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO)

110      CONTINUE

         DO 120 K1=NDI+1,NTENS

            FLOW(K1)=STRESS(K1)*ONESY

120      CONTINUE

C     FLOW(1-6)为 sij/q

C       SOLVE FOR EQUIV STRESS, NEWTON ITERATION

C

         SYIELD=SYIEL0

         DEQPL=0.0

         DO 130 KEWTON=1,NEWTON

            RHS=SMISES-EG3*DEQPL-SYIELD

C         计算屈服面函数误差

            DEQPL=DEQPL+RHS/(EG3+HARD)

C         按(6-13)计算新的等效塑性应变

            CALL AHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(3),NVALUE)

C         由新的等效塑性应变获得新的屈服应力

            IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 140

C         如果满足收敛要求结束,否则再次迭代

130      CONTINUE

         WRITE(6,2) NEWTON

2        FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',

    1        'CONVERGE AFTER ',I3,' ITERATIONS')

140      CONTINUE

         EFFHRD=EG3*HARD/(EG3+HARD)

C

C

C       CALC STRESS AND UPDATE STRAINS

C       计算及更新应力,弹性应变,塑性应变

         DO 150 K1=1,NDI

C        更新正应力分量

            STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO

C        按(6-16)计算应力

            EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO

C        更新塑性应变

            EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO

C        更新弹性应变

150      CONTINUE

         DO 160 K1=NDI+1,NTENS

C        更新剪应力分量

            STRESS(K1)=FLOW(K1)*SYIELD

            EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL

            EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL

160      CONTINUE

         EQPLAS=EQPLAS+DEQPL

         SPD=DEQPL*(SYIEL0+SYIELD)/TWO      

C

C       JACOBIAN

C       计算雅克比矩阵

C

         EFFG=EG*SYIELD/SMISES

         EFFG2=TWO*EFFG

         EFFG3=THREE*EFFG2/TWO

         EFFLAM=(EBULK3-EFFG2)/THREE

         DO 220 K1=1,NDI

            DO 210 K2=1,NDI

               DDSDDE(K2,K1)=EFFLAM

210         CONTINUE

            DDSDDE(K1,K1)=EFFG2+EFFLAM

220      CONTINUE

         DO 230 K1=NDI+1,NTENS

            DDSDDE(K1,K1)=EFFG

230      CONTINUE

         DO 250 K1=1,NTENS

            DO 240 K2=1,NTENS

               DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1)

    1                                       *(EFFHRD-EFFG3)

240         CONTINUE

250      CONTINUE

       ENDIF

     ENDIF

C

C STORE STRAINS IN STATE VARIABLE ARRAY

C

     DO 310 K1=1,NTENS

       STATEV(K1)=EELAS(K1)

       STATEV(K1+NTENS)=EPLAS(K1)

310  CONTINUE

     STATEV(1+2*NTENS)=EQPLAS

C

     RETURN

     END

C

C

     SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE)

C     根据导入的硬化曲线及当前的等效塑性应变,求得屈服应力及屈服应力对等效塑性应变的偏导

     INCLUDE 'ABA_PARAM.INC'

     DIMENSION TABLE(2,NVALUE)

C

C    SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO

     SYIELD=TABLE(1,NVALUE)

     HARD=0.0

C

C   IF MORE THAN ONE ENTRY, SEARCH TABLE

C

     IF(NVALUE.GT.1) THEN

       DO 10 K1=1,NVALUE-1

          EQPL1=TABLE(2,K1+1)

          IF(EQPLAS.LT.EQPL1) THEN

            EQPL0=TABLE(2,K1)

            IF(EQPL1.LE.EQPL0) THEN

               WRITE(6,1)

1              FORMAT(//,30X,'***ERROR - PLASTIC STRAIN MUST BE ',

    1                 'ENTERED IN ASCENDING ORDER')

               CALL XIT

             ENDIF

C

C           CURRENT YIELD STRESS AND HARDENING

C

           DEQPL=EQPL1-EQPL0

           SYIEL0=TABLE(1,K1)

           SYIEL1=TABLE(1,K1+1)

           DSYIEL=SYIEL1-SYIEL0

           HARD=DSYIEL/DEQPL

           SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD

           GOTO 20

           ENDIF

10         CONTINUE

20         CONTINUE

     ENDIF

     RETURN

     END

免责声明:本文系网络转载或改编,未找到原创作者,版权归原作者所有。如涉及版权,请联系删

QR Code
微信扫一扫,欢迎咨询~

联系我们
武汉格发信息技术有限公司
湖北省武汉市经开区科技园西路6号103孵化器
电话:155-2731-8020 座机:027-59821821
邮件:tanzw@gofarlic.com
Copyright © 2023 Gofarsoft Co.,Ltd. 保留所有权利
遇到许可问题?该如何解决!?
评估许可证实际采购量? 
不清楚软件许可证使用数据? 
收到软件厂商律师函!?  
想要少购买点许可证,节省费用? 
收到软件厂商侵权通告!?  
有正版license,但许可证不够用,需要新购? 
联系方式 155-2731-8020
预留信息,一起解决您的问题
* 姓名:
* 手机:

* 公司名称:

姓名不为空

手机不正确

公司不为空