VUMAT状态变量

yinbin_12345

为什么我的子程序不能运行，求帮助！！！

! VUMAT头文件
      subroutine vumat(
! 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,
! Write only (modifiable) variables -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
!
      include 'vaba_param.inc'

      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)
!
      character*80 cmname
      parameter(zero=0.d0,half=0.5,two=2.d0,one=1.d0)
      dimension CFULL(6,6)
      dimension ss(6)

! 输入材料强度参数，用于失效判断
      xt=props(1)                 
      xc=props(2)                  
      yt=props(3)                 
      yc=props(4)
      zt=props(5)
      zc=props(6)               
      s12=props(7)               
      s23=props(8)                 
      s13=props(9)                 
! 输入材料的弹性常数（下标顺序均按照沈观林的教材规定）
! props(13)=emu21,为主泊松比
! props(14)=emu32,props(15)=emu31为次泊松比
      e1=props(10)
      e2=props(11)
      e3=props(12)
      emu12=props(13)
      emu23=props(14)
      emu13=props(15)
      g12=props(16)
      g23=props(17)
      g13=props(18)
! ******Standard strain****
      sn1t=xt/e1
      sn1c=xc/e1
      sn2t=yt/e2
      sn2c=yc/e2
      sn3t=zt/e3
      sn12=s12/g12
      sn23=s23/g23
      sn13=s13/g13
      E11=e1
      E22=e2
      E33=e3
      GG12=g12
      GG23=g23
      GG13=g13
   
! 在初始步给状态参数赋初值，状态参数为0说明材料完好
        NTENS=NDIR+NSHR
       do 100 i = 1,nblock
           
      if(totalTime .eq.0.02)then
       do j=1,6
          stateOld(i,j)=0.0
          strainInc(i,j)=0.0
       end do
      endif

       do j=1,6
          stateNew(i,j)=stateOld(i,j)+strainInc(i,j)
          ss(j)=stateOld(i,j)+strainInc(i,j)
      end do
         
      
        damL=0.0
        damT=0.0
        damZ=0.0
        damTZ=0.0
        damLZ=0.0
        damLT=0.0
      term1=(1/sn1t-1/sn1c)*ss(1)+1/sn1t/sn1c*ss(1)*ss(1)
      term2=(1/sn2t-1/sn2c)*ss(2)+1/sn2t/sn2c*ss(2)*ss(2)
      term3=(1/sn3t-1/sn3c)*ss(3)+1/sn3t/sn3c*ss(3)*ss(3)
      term4=(2*ss(5)/sn23)**2
      term5=(2*ss(6)/sn13)**2
      term6=(2*ss(4)/sn12)**2
      term7=1/(sqrt(sn1t*sn1c*sn2t*sn2c))*ss(1)*ss(2)
     *     +1/(sqrt(sn2t*sn2c*sn3t*sn3c))*ss(2)*ss(3)
     *     +1/(sqrt(sn3t*sn3c*sn1t*sn1c))*ss(1)*ss(3)
      TWD=term1+term2+term3+term4+term5+term6-term7
      
      if(TWD .gt.one) then
         
           if(term1 .ge. term2) then
          term8=term1
          damL=1.0
          damT=0.0
          damZ=0.0
          damTZ=0.0
          damLZ=0.0
          damLT=0.0
          else
          term8=term2
          damL=0.0
          damT=1.0
          damZ=0.0
          damTZ=0.0
          damLZ=0.0
          damLT=0.0
          end if
         
           if(term8 .ge. term3) then
          else
          term8=term3
          damL=0.0
          damT=0.0
          damZ=1.0
          damTZ=0.0
          damLZ=0.0
          damLT=0.0
          end if
         
           if(term8 .ge. term4) then
          else
          term8=term4
          damL=0.0
          damT=0.0
          damZ=0.0
          damTZ=1.0
          damLZ=0.0
          damLT=0.0
          end if
         
           if(term8 .GE. term5) then
          else
          term8=term5
          damL=0.0
          damT=0.0
          damZ=0.0
          damTZ=0.0
          damLZ=1.0
          damLT=0.0
          end if
         
           if(term8 .GE. term6) then
          else
          term8=term6
          damL=0.0
          damT=0.0
          damZ=0.0
          damTZ=0.0
          damLZ=0.0
          damLT=1.0
          end if
           
          else
        end if
      
      stateNew(i,7)=max(stateOld(i,7),damL)  
      stateNew(i,8)=max(stateOld(i,8),damT)  
      stateNew(i,9)=max(stateOld(i,9),damZ)  
      stateNew(i,10)=max(stateOld(i,10),damLT)  
      stateNew(i,11)=max(stateOld(i,11),damTZ)  
      stateNew(i,12)=max(stateOld(i,12),damLZ)  
      
      
              
      e11=e1
      e22=e2
      e33=e3
      xemu12=emu12
      xemu23=emu23
      xemu13= emu13
      gg12=g12
      gg23=g23
      gg13=g13
      
      
      if(stateNew(i,7).EQ.1.0)then
          e11=0.01*e1
          gg13=0.2*g13
          gg12=0.2*g12
       end if
      
       if(stateNew(i,8).EQ.1.0)then
          e22=0.01*e2
          gg23=0.2*g23
          gg12=0.2*g12
       end if
      
       if(stateNew(i,9).EQ.1.0)then
          e33=0.01*e3
          gg23=0.2*g23
          gg13=0.2*g13
       end if
       if(stateNew(i,10).EQ.1.0)then
          e22=0.01*e2
          e33=0.01*e3
          gg23=0.2*g23
          gg12=0.2*g12
          gg13=0.2*g13
       end if
       if(stateNew(i,11).EQ.1.0)then
          e11=0.8*e1
          e33=0.01*e3
          gg23=0.2*g23
          gg13=0.2*g13
       end if
      
       if(stateNew(i,12).EQ.1.0)then
          e11=0.8*e1
          e22=0.01*e2
          gg23=0.2*g23
          gg12=0.2*g12
       end if
      
      
!Calculate stress            
      xemu21 = xemu12 * e22 / e11
      xemu31 = xemu13 * e33 / e11
      xemu32 = xemu23 * e33 / e22
! Initial analysis，assumely pure elastic behaviour

       gg = one / ( one - xemu12*xemu21 - xemu23*xemu32 - xemu31*xemu13
     * - two*xemu21*xemu32*xemu13 )
      CFULL(1,1) = e11 * ( one - xemu23*xemu32 ) * gg
      CFULL(2,2) = e22 * ( one - xemu13*xemu31 ) * gg
      CFULL(3,3) = e33 * ( one - xemu12*xemu21 ) * gg
      CFULL(1,2) = e11 * ( xemu21 + xemu31*xemu23 ) * gg
      CFULL(1,3) = e11 * ( xemu31 + xemu21*xemu32 ) * gg
      CFULL(2,3) = e22 * ( xemu32 + xemu12*xemu31 ) * gg
      CFULL(4,4) = gg12
      CFULL(5,5) = gg23
      CFULL(6,6) = gg13
      CFULL(2,1) = CFULL(1,2)
      CFULL(3,2) = CFULL(2,3)
      CFULL(3,1) = CFULL(1,3)  
      
      
      stressNew(i,1)=CFULL(1,1)*stateNew(i,1)+
     &               CFULL(1,2)*stateNew(i,2)+
     &               CFULL(1,3)*stateNew(i,3)
      stressNew(i,2)=CFULL(2,1)*stateNew(i,1)+
     &               CFULL(2,2)*stateNew(i,2)+
     &               CFULL(2,3)*stateNew(i,3)
      stressNew(i,3)=CFULL(3,1)*stateNew(i,1)+
     &               CFULL(3,2)*stateNew(i,2)+
     &               CFULL(3,3)*stateNew(i,3)
      stressNew(i,4)=CFULL(4,4)*2*stateNew(i,4)
      stressNew(i,5)=CFULL(5,5)*2*stateNew(i,5)
      stressNew(i,6)=CFULL(6,6)*2*stateNew(i,6)
                     
  100 continue
      return
      end
!
为什么删去标红的那几行就能运行了，求帮助

weiliang174

楼主在设置ss(j)的时候有错误，交流qq569250598

weiliang174

你的问题搞定没？我也是做类似的，共同交流，qq569250598

weiliang174

顶一下，楼主问题解决没？

gxp-dhu

楼主问题解决没？共同探讨，QQ569250598

jcz1996

跟着学习

枯藤老树昏鸦

学习一下

枯藤老树昏鸦

学习一下