UMAT/VUMAT入门攻略-Davvife留个Simwe的纪念帖

crest

学习了   谢谢

mathew12345

花了好几天，看了论坛大侠的分享帖，终于搞定子程序调试了，接下来就是学习一下用法

liangxiuwen

谢谢楼主啦   先下喽

mathew12345

必须置顶啊

子程序有用的东东太少了

tyrone

感谢感谢，刚要学习

yuyuchen

感谢版主所提供资料！

wangcongkang

在这一经典的vumat好贴中，我弱弱地向各位版主提一个问题。 我用vumat编写的johnson-cook模型，运行时遇到693的错误。不清楚程序到底哪里有误。请daivvife 等版主帮忙看看， 真的感激不尽！

c User subroutine VUMAT
      subroutine vumat (
c Read only -
     1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     2 stepTime, totalTime, dt, cmname, coordMp, charLength,
     3 props, density, strainInc, relSpinInc,
     4 tempOld, stretchOld, defgradOld, fieldOld,
     5 stressOld, stateOld, enerInternOld, enerInelasOld,
     6 tempNew, stretchNew, defgradNew, fieldNew,
c Write only -
     7 stressNew, stateNew, enerInternNew, enerInelasNew )
c
      include 'vaba_param.inc'
c
      dimension coordMp(nblock,*), charLength(nblock), props(nprops),
     1 density(nblock), strainInc(nblock,ndir+nshr),
     2 relSpinInc(nblock,nshr), tempOld(nblock),
     3 stretchOld(nblock,ndir+nshr),
     4 defgradOld(nblock,ndir+nshr+nshr),
     5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     6 stateOld(nblock,nstatev), enerInternOld(nblock),
     7 enerInelasOld(nblock), tempNew(nblock),
     8 stretchNew(nblock,ndir+nshr),
     9 defgradNew(nblock,ndir+nshr+nshr),
     1 fieldNew(nblock,nfieldv),
     2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     3 enerInternNew(nblock), enerInelasNew(nblock)
c
      character*80 cmname
c
      parameter ( zero = 0.d0, one = 1.d0, two = 2.d0,
     1 third = 1.d0 / 3.d0, half = 0.5d0, op5 = 1.5d0 )
c
c For plane strain, axisymmetric, and 3D cases using
c the J2 Mises Plasticity with linear hardening.
c The state variable is stored as:
c STATE(*,1) = equivalent plastic strain
c STATE(*,2) = plastic strain rate
c
c
c User needs to input
c props(1) Young's modulus
c props(2) Poisson's ratio
c props(3) A
c prpps(4) B
c prpps(5) n
c prpps(6) C
c prpps(7) m
c
c
       e = props(1)
       xnu = props(2)
       A = props(3)
       B = props(4)
       EN = props(5)
       C = props(6)
       EM = props(7)
       twomu = e / ( one + xnu )
       alamda = xnu * twomu / ( one - two * xnu )
       thremu = op5 * twomu
     1
      if ( stepTime .eq. zero ) then
        do k = 1, nblock
          trace = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
          stressNew(k,1) = stressOld(k,1)
     1     + twomu * strainInc(k,1) + alamda * trace
          stressNew(k,2) = stressOld(k,2)
     1     + twomu * strainInc(k,2) + alamda * trace
          stressNew(k,3) = stressOld(k,3)
     1     + twomu * strainInc(k,3) + alamda * trace
          stressNew(k,4)=stressOld(k,4) + twomu * strainInc(k,4)
          if ( nshr .gt. 1 ) then
             stressNew(k,5)=stressOld(k,5) + twomu * strainInc(k,5)
             stressNew(k,6)=stressOld(k,6) + twomu * strainInc(k,6)
          end if
        end do
      else
        do k = 1, nblock   
c
c calculate hard and yieldOld
c
          peeqOld=stateOld(k,1)
c         
          if(peeqOld.eq.zero) then
            yieldOld = A
          else
            hard = EN * B * peeqOld ** (EN-1)
            yieldOld = A + B * peeqOld ** EN
          end if
            
          trace = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
          s11 = stressOld(k,1) + twomu * strainInc(k,1) + alamda *
     1      trace
          s22 = stressOld(k,2) + twomu * strainInc(k,2) + alamda *
     1      trace
          s33 = stressOld(k,3) + twomu * strainInc(k,3) + alamda *
     1      trace
          s12 = stressOld(k,4) + twomu * strainInc(k,4)
          if ( nshr .gt. 1 ) then
          s13 = stressOld(k,5) + twomu * strainInc(k,5)
          s23 = stressOld(k,6) + twomu * strainInc(k,6)
          end if
c
          smean = third * ( s11 + s22 + s33 )
          s11 = s11 - smean
          s22 = s22 - smean
          s33 = s33 - smean
            if ( nshr .eq. 1 ) then
              vmises = sqrt( op5*
     1       (s11*s11+s22*s22+s33*s33+two*s12*s12) )
            else
              vmises = sqrt( op5 * ( s11 * s11 + s22 * s22 + s33 * s33 +
     1        two * s12 * s12 + two * s13 * s13 + two * s23 *
     2         s23 ) )
            end if
c
            tvp = C * log(stateOld(k,2))
            tvp1 = tvp + one
            hard1 = hard * tvp1 + yieldOld * C / stateOld(k,2)
            yieldOld = yieldOld * tvp1
            sigdif = vmises - yieldOld
            facyld = zero
            if ( sigdif .gt. zero ) facyld = one
              deqps = facyld * sigdif / (thremu + hard1)
c
c  Update the state variable
              stateNew(k,1) = stateOld(k,1) + deqps
              stateNew(k,2) = deqps/dt
c            
c
c
c
c  Update the stress
              yieldNew = yieldOld + hard1 * deqps
              factor = yieldNew / ( yieldNew + thremu * deqps )
              stressNew(k,1) = s11 * factor + smean
              stressNew(k,2) = s22 * factor + smean
              stressNew(k,3) = s33 * factor + smean
              stressNew(k,4) = s12 * factor
             if ( nshr .gt. 1 ) then
                stressNew(k,5) = s13 * factor
                stressNew(k,6) = s23 * factor
             end if
c
c Update the specific internal energy -
c
             if ( nshr .eq. 1 ) then
              stressPower = half * (
     1         ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +
     2         ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +
     3         ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3) ) +
     4         ( stressOld(k,4) + stressNew(k,4) ) * strainInc(k,4)
             else
              stressPower = half * (
     1         ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +
     2         ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +
     3         ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3) ) +
     4         ( stressOld(k,4) + stressNew(k,4) ) * strainInc(k,4) +
     5         ( stressOld(k,5) + stressNew(k,5) ) * strainInc(k,5) +
     6         ( stressOld(k,6) + stressNew(k,6) ) * strainInc(k,6)
             end if
          enerInternNew(k) = enerInternOld(k) + stressPower / density(k)
c
c Update the dissipated inelastic specific energy -
          plasticWorkInc = yieldNew * deqps
          enerInelasNew(k) = enerInelasOld(k)
     1      + plasticWorkInc / density(k)
        end do
      end if

      return
      end
c   

佐净

貌似都是UMAT哈，最近学习VUMAT，但是貌似相关的东西比UMAT少的多哈！

张屠户

有用VUMAT做GURSON模型的吗？最近写了个这方面的子程序，可是子程序里的迭代怎么都不收敛，有没有高手可以指点一下？

冰川纪的火

支持一下   我也是新手  

喝茶的羊

想弱弱的问一句：vumat怎么样输出自己想要的数据
应该就是write（*，*）
但是这个只能输出到屏幕上
我在vumat中添加了这样的语句，算了半天还没有输出
不知道是语句的问题，还是说输出数据只能一步一步的输
还是说只是最后一步才输出数据呢，反正各种不明白
希望有大侠指导一下我这个小菜鸟吧

billowriver

一直觉得大家没有必要纠结于UMAT或者VUMAT本身上面，其实说白了就是根据输入变量（基本上你想要的都会有），求Cauchy Stress （不是其它任何应力），隐式的话由于迭代的需要还要多算个jacobi，
所以，问题最后都归结成本构和算法的问题，如果要讨论还是讨论本构和算法更好些。LS-Dyna的本构子程序，把输入，输出变量的名称改一下就直接可以用在Vumat上，不会有任何问题，甚至你都可以用Matlab测试你的本构子程序。

billowriver

wangcongkang 发表于 2011-9-29 15:27
在这一经典的vumat好贴中，我弱弱地向各位版主提一个问题。 我用vumat编写的johnson-cook模型，运行时遇到6 ...

一般693错误是被0除了，把变量用print*, xxx 一行一行打印到屏幕上进行调试.
塑性本构是需要迭代求解等效塑性应变增量的，没有仔细看，但是好像没有迭代，你可以查一下目前最常用的半径回归法(radius return method)。

billowriver

敦诚 发表于 2011-7-26 19:01
好吧，支持美女版主，提个问题。vumat中一般会涉及大变形，这个时候应变总量无法通过简单地线性叠加来实现 ...

需要用stretch tensor，U，由于ABAQUS采用的是随材料转动的材料坐标系，所以最好不要用deformation gradient，然后求解对数应变，lnU，求张量的对数需要用到SVD分解，Vumat里面提供的应变增量其实就是对数应变的增量，

billowriver

北鹰南飞 发表于 2011-7-24 10:17
1# Davvife
请教:正好有个历史遗留问题一直没有解决，如图：

Umat是加在积分点上的，如果就一个单元，需要增加积分点的数量，但是如果想在Umat里面区分单元和积分点，需要用到整体变量，由于每个积分点都会调用一次Umat，所以可以根据调用的次数区分积分点和单元

billowriver

Davvife 发表于 2011-7-26 19:55
谢谢敦大的支持啊！我真的还遇到过这个问题啊！嘻嘻嘻
这是一个好问题，很多人遇到过，但是始终没有人详 ...

无论Umat还是Vumat，Abaqus都采用随材料变形一起的旋转坐标系，都是剔除了旋转部分，不同的是到底有多少旋转部分，Umat和Vumat采用不同的率进行计算，一般Vumat采用Green Naghdi Rate，而Umat采用Jaumann rate。没有哪个对哪个错，只能说哪个更适合你的材料本构。
Umat的麻烦在于它的时间增量一般比较大，所以每一个增量步都需要区分初始和结束时的变形梯度等物理量。

billowriver

lyk_302 发表于 2011-8-2 03:25
在看umat帮助时候有这么一段话：
Large volume changes with geometric nonlinearity：
....

很多情况下不需要特别精确的Jacobi，直接给出线弹性的就行，它只影响收敛速度，不影响最后的结果

lungo999

张屠户 发表于 2011-10-28 20:52
有用VUMAT做GURSON模型的吗？最近写了个这方面的子程序，可是子程序里的迭代怎么都不收敛，有没有高手可以 ...

不收敛很正常，一下子收敛了才不正常，Gurson模型是模拟空穴增长的模型吧，模型应该蛮复杂的，算法最重要，还要注意矩阵和张量运算的区别。

喝茶的羊

喝茶的羊 发表于 2011-11-1 17:24
想弱弱的问一句：vumat怎么样输出自己想要的数据
应该就是write（*，*）
但是这个只能输出到屏幕上

啊啊啊，我被忽略鸟

一一沙尘暴一一

个人认为这篇帖子应该长期置顶啊！