找回密码
 注册
Simdroid-非首页
查看: 784|回复: 11

[子程序] 能帮忙调试一下johnson cook vumat子程序吗

[复制链接]
发表于 2015-4-20 17:19:17 | 显示全部楼层 |阅读模式 来自 湖北武汉
RT,我是参照庄茁老师的umat,然后改编的vumat程序,但是通不过,不知道是为何,大侠能否出现,帮忙看看啊

      subroutine vumat(
C Read only (unmodifiable)variables -
     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 (modifiable) variables -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
C
      include 'vaba_param.inc'
C
      dimension props(nprops), density(nblock), coordMp(nblock,*),
     1  charLength(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),
     8  defgradNew(nblock,ndir+nshr+nshr),
     9  fieldNew(nblock,nfieldv),
     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock)
C
      character*80 cmname
C
      dimension eelas(6),eplas(6),flow(6),strain(6)
      parameter (one=1.0D0,two=2.0D0,three=3.0D0,six=6.0D0, half =0.5d0,
     1      third=1.d0/3.d0)

      data newton,toler/40,1.D-6/
      con=sqrt(2.d0/3.d0)

C     props(1) - yang's Modulus
C     props(2) - possion radio
C     props(3) - inelastic heat fraction
C     parameters of johnson cook model
C         props(4) - A
C         props(5) - B
C         props(6) - n
C         props(7) - C
C         props(8) - m
C

C state(*,1)=elastic strain component 11
C state(*,2)=elastic strain component 22
C state(*,3)=elastic strain component 33
C state(*,4)=elastic strain component 12
C state(*,5)=elastic strain component 23
C state(*,6)=elastic strain component 31
C
C state(*,7)=plastic strain component 11
C state(*,8)=plastic strain component 22
C state(*,9)=plastic strain component 33
C state(*,10)=plastic strain component 12
C state(*,11)=plastic strain component 23
C state(*,12)=plastic strain component 31
C
C state(*,13)=equivalent plastic strain
C
C state(*,14)=total strain component 11
C state(*,15)=total strain component 22
C state(*,16)=total strain component 33
C state(*,17)=total strain component 12
C state(*,18)=total strain component 23
C state(*,19)=total strain component 31
C
C state(*,20)=equivalent strain

C
C state(*,21)=total stress component 11
C state(*,22)=total stress component 22
C state(*,23)=total stress component 33
C state(*,24)=total stress component 12
C state(*,25)=total stress component 23
C state(*,26)=total stress component 31
C
C state(*,27)==smises stress

      emod=props(1)
      enu=props(2)
      eg=emod/(one+enu)/two
      eg2=eg*two
      eg6=eg*six
      elam=eg2*(emod-eg2)/(eg6-two*emod)  

!      if ( stepTime .eq. zero ) then
!          do i = 1, nblock
!          trace = strainInc(i,1) + strainInc(i,2) + strainInc(i,3)
!          stressNew(i,1) = stressOld(i,1)
!     *       + eg2 * strainInc(i,1) + elam * trace
!          stressNew(i,2) = stressOld(i,2)
!     *       + eg2 * strainInc(i,2) + elam * trace
!          stressNew(i,3) = stressOld(i,3)
!     *       + eg2 * strainInc(i,3) + elam * trace
!          stressNew(i,4)=stressOld(i,4) + eg2 * strainInc(i,4)
!          stressNew(i,5)=stressOld(i,5) + eg2 * strainInc(i,5)
!          stressNew(i,6)=stressOld(i,6) + eg2 * strainInc(i,6)
!         
!          end do
!      else


      do  i=1,nblock
         trace=strainInc(i,1)+strainInc(i,2)+strainInc(i,3)

         sig1=stressOld(i,1)+elam*trace+eg2*strainInc(i,1)
         sig2=stressOld(i,2)+elam*trace+eg2*strainInc(i,2)
         sig3=stressOld(i,3)+elam*trace+eg2*strainInc(i,3)
         sig4=stressOld(i,4) +eg2*strainInc(i,4)
         sig5=stressOld(i,5) +eg2*strainInc(i,5)
         sig6=stressOld(i,6) +eg2*strainInc(i,6)

         smean=third*(sig1+sig2+sig3)
         sig11=sig1-smean
         sig22=sig2-smean
         sig33=sig3-smean
         smisesold=(one/con)*sqrt(sig11**2+sig22**2+sig33**2+two*sig4**2
     1      +two*sig5**2+two*sig6**2)

         do k1=1,6
           eelas(k1)=stateOld(i,k1)+strainInc(i,k1)
           eplas(k1)=stateOld(i,k1+6)
           stateNew(i,k1+13)=stateOld(i,k1+13)+strainInc(i,k1)
         enddo

         eqplas=stateNew(i,13)

         str=(stateNew(i,14)+stateNew(i,15)+stateNew(i,16))/3
         str1=stateNew(i,14)-str
         str2=stateNew(i,15)-str
         str3=stateNew(i,16)-str
         str4=stateNew(i,17)
         str5=stateNew(i,18)
         str6=stateNew(i,19)


c        equstrain:equivalent strain
         equstrain=con*(str1**2+str2**2+str3**2+two*str4**2+
     1 two*str5**2+two*str6**2)
         stateNew(i,20)=equstrain

C
C        calculate miser stress
C  

C        call userhard subroutine, get harding rate and yield stress         

         call userhard(syiel0,hard,eqplas,props(4))
C        determine if yielding is actived
         if (smises > (1+toler)*syiel0) then
C          if plastic deformation has occured,  determine the flow direction        
           shydr0=(stressNew(i,1)+stressNew(i,2)+stressNew(i,3))/3.d0
           onesy=one/smises
           do k1=1,3
              flow(k1)=onesy*(stressNew(i,k1)-shydr0)
           enddo
           do k1=4,6
              flow(k1)=onesy*stressNew(i,k1)
           enddo

C          read the johnson cook parameters           
           A=props(4)
           B=props(5)
           EN=props(6)
           C=props(7)
           EM=props(8)

C          begin newton iteration
           syield = syiel0
           deqpl=(smises-syield)/eg3
           dstres=toler*syiel0/eg3
           deqmin=half*dt*exp(1.d0-4/c)

           do knewton=1,newton
              deqpl=max(deqpl,deqmin)
              call userhard(syield,hard,eqplas+deqpl,props(4))
              tvp=c*log(deqpl/dt)
              tvp1=tvp+one
              hard1=hard*tvp1+syield*c/deqpl
              syield=syield*tvp1
              rhs=smises-eg3*deqpl-syield
              deqpl=deqpl+rhs/(eg3+hard1)
              if(abs(rhs/eg3) < dstres) goto 140

           enddo
140       continue

C          update stress and strain
C          whether the following codes is true?divited by 2?
           do k1=1,3
              stressNew(i,k1) = flow(k1)*syield+shydr0
              eelas(k1) = eelas(k1)-three * flow(k1)*deqpl/two
              eplas(k1) = eplas(k1)+three * flow(k1)*deqpl/two
           enddo

           do k1=4,6
              stressNew(i,k1) = flow(k1)*syield
              eelas(k1) = eelas(k1)-three * flow(k1)*deqpl
              eplas(k1) = eplas(k1)+three * flow(k1)*deqpl
           enddo               

           eqplas=eqplas+deqpl
           spd=edqpl*(syiel0+syield)/two
           rpl=props(3)*spd/dt

         endif

         smises=(stressNew(i,1)-stressNew(i,2)) *
     1             (stressNew(i,1)-stressNew(i,2)) +
     2          (stressNew(i,2)-stressNew(i,3)) *
     3             (stressNew(i,2)-stressNew(i,3)) +
     4          (stressNew(i,3)-stressNew(i,1)) *
     5             (stressNew(i,3)-stressNew(i,1))
         do  k1=4,6
              smises=smises+six*stressNew(i,k1)*stressNew(i,k1)
         enddo
         smises=sqrt(smises/two)

         do  k1=1,6
             stateNew(i,k1) = eelas(k1)
             stateNew(i,k1+6) = eplas(k1)
         enddo

         do  k1=1,6
             stateNew(i,20+k1) = stressNew(i,k1)
         enddo
         stateNew(i,13)=eqplas
         stateNew(i,27)=smises


         stressPower = half * (
     *    ( stressOld(i,1) + stressNew(i,1) ) * strainInc(i,1) +
     *    ( stressOld(i,2) + stressNew(i,2) ) * strainInc(i,2) +
     *    ( stressOld(i,3) + stressNew(i,3) ) * strainInc(i,3) ) +
     *    ( stressOld(i,4) + stressNew(i,4) ) * strainInc(i,4) +
     *    ( stressOld(i,5) + stressNew(i,5) ) * strainInc(i,5) +
     *    ( stressOld(i,6) + stressNew(i,6) ) * strainInc(i,6)

         enerInternNew(i) = enerInternOld(i) + stressPower / density(i)
C
C Update the dissipated inelastic specific energy -
C
         plasticWorkInc = half * ( smises + smisesold ) * deqps
         enerInelasNew(i) = enerInelasOld(i)
     *    + plasticWorkInc / density(i)
      enddo

!     endif
      return
      end

C
C
C      

      subroutine userhard(syield,hard,eqplas,table)

      include 'vaba_param.inc'
      dimension table(3)
      A=table(1)
      B=table(2)
      EN=table(3)
      hard=0.0
C
C     calcular current yield stress and hardening rate
C
      if (eqplas .eq. 0.0) then
         syield =A
      else
         hard = EN*B*eqplas**(EN-1)
         syield=A+B*eqplas**EN
      endif
      return
      end

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
 楼主| 发表于 2015-4-20 17:34:38 | 显示全部楼层 来自 湖北武汉
Simdroid开发平台
错误是这样的,八成子程序出问题
Error in job Job-1: Bad Material definition in element number 1 instance PART-1-1: zero or negative initial dilatational modulus caused by bad material data. Please check your material input and any initial conditions if necessary.
Job Job-1: Abaqus/Explicit Packager aborted due to errors.
Error in job Job-1: Abaqus/Explicit Packager exited with an error - Please see the  status file for possible error messages if the file exists.
Job Job-1 aborted due to errors.
回复 不支持

使用道具 举报

发表于 2015-5-18 11:03:03 | 显示全部楼层 来自 大连理工大学
vumat里面是不允许出现迭代计算的,如果确实无法避免的话一定要声明迭代次数,具体你可以参考Writing_User_Subroutines_with_ABAQUS_0这里面关于vumat编译方面的注意事项
回复 不支持

使用道具 举报

 楼主| 发表于 2015-5-21 19:26:57 | 显示全部楼层 来自 湖北武汉
kang.123 发表于 2015-5-18 11:03
vumat里面是不允许出现迭代计算的,如果确实无法避免的话一定要声明迭代次数,具体你可以参考Writing_User_ ...

迭代肯定是允许的,我以前编其他程序的时候,用过,是没问题的,就是我上面的程序里面出现了为零的情况,导致变形巨大,从而不收敛
回复 不支持

使用道具 举报

发表于 2015-7-6 08:51:40 | 显示全部楼层 来自 大连理工大学
楼主你好,不知道你用没用过庄茁老师的JC本构的umat子程序,我现在正在学习,但是不知道怎么回事在使用的时候总是在弹性过渡到塑性时出现不收敛而终止计算。
回复 不支持

使用道具 举报

发表于 2018-10-9 14:20:24 | 显示全部楼层 来自 重庆沙坪坝区
我也遇到这样的问题,请问楼主是怎么解决的?
回复 不支持

使用道具 举报

发表于 2018-10-15 09:14:15 | 显示全部楼层 来自 宁波大学
最近我看庄茁的这个论文,我按照他一模一样的模型,没有做出来,请问有人能指导一下
回复 不支持

使用道具 举报

发表于 2018-10-15 19:05:47 | 显示全部楼层 来自 重庆沙坪坝区
请问楼主问题解决了吗
回复 不支持

使用道具 举报

发表于 2022-4-29 19:18:22 | 显示全部楼层 来自 黑龙江
您好!请问您有JC的切削vumat吗,望指教!
回复 不支持

使用道具 举报

发表于 2022-6-16 19:37:56 | 显示全部楼层 来自 大连理工大学西山生活区
请问楼主这个问题解决了吗,我用vumat也报相同的错误,但是我找不到程序哪里出错了
回复 不支持

使用道具 举报

发表于 2022-12-2 23:25:44 | 显示全部楼层 来自 黑龙江
请问楼主这个问题解决了吗,我最近也在研究这个
回复 不支持

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-3-29 23:26 , Processed in 0.043420 second(s), 12 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表