各向异性弹塑性VUMAT

HXD900910

本人参照帮助文档的例子写了一个各向异性的弹塑性VUMAT子程序，刚计算就出现单元扭曲，停止计算。存在以下几个疑问：1、子程序是否完整运行了一次？（运行过后，算法和参数问题导致出现单元扭曲？）2、刚调用子程序就出错跳出，子程序并没有完整运行？刚接触子程序，好多问题都不太理解，还望各位前辈能帮忙指点一下，谢谢！

HXD900910

这是我写的程序
! 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
      parameter ( zero = 0., one = 1., two = 2., three = 3.,
     *       half = 0.5, op5 = 4./ 3., a11 = 0.0033, a22 = 1,
     *       a33 = 1, a12 = -0.015, a13 = -0.015, a23 = -0.65,
     *       a44 = 6.9, a55 = 9.1, a66 = 6.9)
C
C The state variable is stored as:
C
C     STATE(*,1) = equivalent plastic strain
C       
! 输入材料的弹性常数
      E1   =   props(1)   
      E2   =   props(2)
      E3   =   props(3)
      v12  =   props(4)
      v13  =   props(5)
      v23  =   props(6)
      G12  =   props(7)
      G13  =   props(8)
      G23  =   props(9)
C Elastic constant
C 初始刚度常数
      v21   =   v12*E2/E1
      v31   =   v13*E3/E1
      v32   =   v23*E3/E2
      delta =   (1-v12*v21-v23*v32-v31*v13-2*v21*v32*v13)/(E1*E2*E3)
      D1    =   (1-v23*v32)/(E2*E3*delta)
      D12   =   (v21+v31*v23)/(E2*E3*delta)
      D21   =   (v12+v32*v13)/(E1*E3*delta)
      D13   =   (v31+v21*v32)/(E2*E3*delta)
      D31   =   (v13+v12*v23)/(E1*E2*delta)
      D2    =   (1-v13*v31)/(E1*E3*delta)
      D23   =   (v32+v12*v31)/(E1*E3*delta)
      D32   =   (v23+v21*v13)/(E1*E2*delta)
      D3    =   (1-v12*v21)/(E1*E2*delta)
      D4    =   two*G12
      D5    =   two*G23  
      D6    =   two*G13
      nvalue = nprops/2-1
C
      if ( stepTime .eq. zero ) then
        do k = 1, nblock
          stressNew(k,1) = stressOld(k,1)+ D1 * strainInc(k,1)
     *            + D12 * strainInc(k,2)+ D13 * strainInc(k,3)
          stressNew(k,2) = stressOld(k,2)+ D21 * strainInc(k,1)
     *            + D2 * strainInc(k,2)+ D23 * strainInc(k,3)
          stressNew(k,3) = stressOld(k,3)+ D31 * strainInc(k,1)
     *            + D32 * strainInc(k,2)+ D3 * strainInc(k,3)
          stressNew(k,4)=stressOld(k,4) + D4 * strainInc(k,4)
          stressNew(k,5)=stressOld(k,5) + D5 * strainInc(k,5)
          stressNew(k,6)=stressOld(k,6) + D6 * strainInc(k,6)
        end do
      else
        do k = 1, nblock
          peeqOld=stateOld(k,1)
          call vuhard(yieldOld, hard, peeqOld, props(10), nvalue)
          s11 = stressOld(k,1)+ D1 * strainInc(k,1)
     *            + D12 * strainInc(k,2)+ D13 * strainInc(k,3)
          s22 = stressOld(k,2)+ D21 * strainInc(k,1)
     *            + D2 * strainInc(k,2)+ D23 * strainInc(k,3)
          s33 = stressOld(k,3)+ D31 * strainInc(k,1)
     *            + D32 * strainInc(k,2)+ D3 * strainInc(k,3)
          s12 = stressOld(k,4) + D4 * strainInc(k,4)
          s23 = stressOld(k,5) + D5 * strainInc(k,5)
          s31 = stressOld(k,6) + D6 * strainInc(k,6)
! 计算等效应力
          fyield=half*(a11*s11**2+a22*s22**2+a33*s33**2)
     *          +a12*s11*s22+a13*s11*s33+a23*s22*s33  
     *          +a44*s12**2+a55*s23**2+a66*s31**2  
            eqstress=sqrt(three*fyield)
          p11=a11*s11+a12*s22+a13*s33
          p22=a12*s11+a22*s22+a23*s33
                  p33=a13*s11+a23*s22+a33*s33
                  p12=2*a44*s12
                  p23=2*a55*s23
                  p13=2*a66*s31
                  a=D1*p11+D21*p22+D31*p33
                  b=D12*p11+D2*p22+D32*p33
                  c=D13*p11+D23*p22+D33*p33
                  d=D4*p12
                  e=D5*p23
                  f=D6*p13
          u=D1*p11+D12*p22+D13*p33
                  v=D21*p11+D2*p22+D23*p33
                  w=D31*p11+D32*p22+D3*p33
                  x=D4*p12
                  y=D5*p23
                  z=D6*p13
          M=a*p11+b*p22+c*p33+d*p12+e*p23+f*p13
     *       +op5*fyield*hard
          sigdif = eqstress - yieldOld
          facyld = zero
          if ( sigdif .gt. zero ) facyld = one
          deqps = facyld*(a*strainInc(k,1)+b*strainInc(k,2)
     *          +c*strainInc(k,3)+d*strainInc(k,4)
     *          +e*strainInc(k,5)+f*strainInc(k,6))
     *          *(s11*p11+s22*p22+s33*p33+s12*p12+s23*p23
     *          +s31*p13) / ( M * eqstress )     
C
C Update the stress
C
          yieldNew = yieldOld + hard * deqps
          D1_EQ = D1-facyld*a*u/M
                  D12_EQ=D12-facyld*b*u/M
                  D13_EQ=D13-facyld*c*u/M
                  D14_EQ=   -facyld*d*u/M
                  D15_EQ=   -facyld*e*u/M
                  D16_EQ=   -facyld*f*u/M
                  D21_EQ=D21-facyld*a*v/M
                  D2_EQ = D2-facyld*b*v/M
                  D23_EQ=D23-facyld*c*v/M
                  D24_EQ=   -facyld*d*v/M
                  D25_EQ=   -facyld*e*v/M
                  D26_EQ=   -facyld*f*v/M
                  D31_EQ=D31-facyld*a*w/M
                  D32_EQ=D32-facyld*b*w/M
                  D3_EQ = D3-facyld*c*w/M
                  D34_EQ=   -facyld*d*w/M
                  D35_EQ=   -facyld*e*w/M
                  D36_EQ=   -facyld*f*w/M
                  D41_EQ=   -facyld*a*x/M
                  D42_EQ=   -facyld*b*x/M
                  D43_EQ=   -facyld*c*x/M
                  D4_EQ = D4-facyld*d*x/M
                  D45_EQ=   -facyld*e*x/M
                  D46_EQ=   -facyld*f*x/M
                  D51_EQ=   -facyld*a*y/M
                  D52_EQ=   -facyld*b*y/M
                  D53_EQ=   -facyld*c*y/M
                  D54_EQ=   -facyld*d*y/M
                  D5_EQ = D5-facyld*e*y/M
                  D56_EQ=   -facyld*f*y/M
                  D61_EQ=   -facyld*a*z/M
                  D62_EQ=   -facyld*b*z/M
                  D63_EQ=   -facyld*c*z/M
                  D64_EQ=   -facyld*d*z/M
                  D65_EQ=   -facyld*e*z/M
                  D6_EQ = D6-facyld*f*z/M
          stressNew(k,1) = stressOld(k,1)
     *            + D1_EQ * strainInc(k,1)+ D12_EQ * strainInc(k,2)
     *            + D13_EQ * strainInc(k,3)+ D14_EQ * strainInc(k,4)
     *            + D15_EQ * strainInc(k,5)+ D16_EQ * strainInc(k,6)
          stressNew(k,2) = stressOld(k,2)
     *            + D21_EQ * strainInc(k,1)+ D2_EQ * strainInc(k,2)
     *            + D23_EQ * strainInc(k,3)+ D24_EQ * strainInc(k,4)
     *            + D25_EQ * strainInc(k,5)+ D26_EQ * strainInc(k,6)
          stressNew(k,3) = stressOld(k,3)
     *            + D31_EQ * strainInc(k,1)+ D32_EQ * strainInc(k,2)
     *            + D3_EQ * strainInc(k,3)+ D34_EQ * strainInc(k,4)
     *            + D35_EQ * strainInc(k,5)+ D36_EQ * strainInc(k,6)
          stressNew(k,4) = stressOld(k,4)
     *            + D41_EQ * strainInc(k,1)+ D42_EQ * strainInc(k,2)
     *            + D43_EQ * strainInc(k,3)+ D4_EQ * strainInc(k,4)
     *            + D45_EQ * strainInc(k,5)+ D46_EQ * strainInc(k,6)
          stressNew(k,5) = stressOld(k,5)
     *            + D51_EQ * strainInc(k,1)+ D52_EQ * strainInc(k,2)
     *            + D53_EQ * strainInc(k,3)+ D54_EQ * strainInc(k,4)
     *            + D5_EQ * strainInc(k,5)+ D56_EQ * strainInc(k,6)
          stressNew(k,6) = stressOld(k,6)
     *            + D61_EQ * strainInc(k,1)+ D62_EQ * strainInc(k,2)
     *            + D63_EQ * strainInc(k,3)+ D64_EQ * strainInc(k,4)
     *            + D65_EQ * strainInc(k,5)+ D6_EQ * strainInc(k,6)
C
C Update the state variables
C
         stateNew(k,1) = stateOld(k,1) + deqps
C
C Update the specific internal energy -
C
         stressPower = half * (
     *     ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +
     *     ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +
     *     ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3) ) +
     *     ( stressOld(k,4) + stressNew(k,4) ) * strainInc(k,4) +
     *     ( stressOld(k,5) + stressNew(k,5) ) * strainInc(k,5) +
     *     ( stressOld(k,6) + stressNew(k,6) ) * strainInc(k,6)
         enerInternNew(k) = enerInternOld(k) + stressPower / density(k)
C
C Update the dissipated inelastic specific energy -
C
         plasticWorkInc = eqstress * deqps
         enerInelasNew(k) = enerInelasOld(k)
     *                   + plasticWorkInc / density(k)
        end do
      end if
C
      return
      end




      subroutine vuhard(syield, hard, eqplas, table, nvalue)
      include 'vaba_param.inc'
c
      dimension table(2, nvalue)
c
      parameter(zero=0.d0)
c
c    set yield stress to last value of table, hardening to zero
c
      syield=table(1, nvalue)
      hard=zero
c
c    if more than one entry, search table
c
      if(nvalue.gt.1) then
        do k1=1, nvalue-1
           eqpl1=table(2,k1+1)
           if(eqplas.lt.eqpl1) then
             eqpl0=table(2, k1)
c
c    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 10
           endif
         end do
   10 continue
       endif
       return
       end

HXD900910

程序与公式推导过程已上传，希望各位前辈能够指点一下，非常感谢

zhdming

请教一下，你塑性各向异性的UMAT最后调试出来了吗

cohesive1

最近也在做各向异性弹塑性，好像听别人说要用牛顿迭代求等效塑性应变，但是我搞不定

zzzzzzz11

你好，我想问下如果直接得到了各向异性弹塑性曲线，如何写vumat呢，你写的程序其塑性是直接用公式定义的