急急急！求大师关于混凝土hjc模型的ABAQUS的vumat二次开发程...

土木嘉年华

最近在做混凝土的损伤研究，用到了hjc模型，现有以下两个问题向大神请教
（1）在调用子程序时，提交成功，但出现这样的错误，该怎么处理？
The job input file "static-123.inp" has been submitted for analysis.
Job static-123: Analysis Input File Processor completed successfully.
ERROR in job messaging system: Error in connection to analysis
Error in job static-123: The executable D:\SIMULIA\Abaqus\6.10-1\exec\package.exe aborted with system error code 5. Please check the .dat, .msg, and .sta files for error messages if the files exist.  If there are no error messages and you cannot resolve the problem, please run the command "abaqus job=support information=support" to report and save your system information.  Use the same command to run Abaqus that you used when the problem occurred.  Please contact your local Abaqus support office and send them the input file, the file support.log which you just created, the executable name, and the error code.
Job static-123 aborted due to errors.
（2）子程序在没有加入输出调试时，计算结果严重偏离实际情况，怀疑子程序有问题，希望有大神可以帮忙解决一下
现附上子程序文件和inp文件，拜托了，急求，多谢了！      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,three=3.d0,
     *  third=1.d0/3.d0,half=0.5d0,op5=1.5d0,toler=10e-3,
     *  twothird=2.d0/3.d0,NEWTON=10,fifty=50,hundred=100)

      dimension strial(ndir+nshr),ssa(ndir+nshr),etrie(ndir+nshr)
     1 ,ssb(ndir+nshr)
      integer k,i,iii
      double precision e,xnu,A,B,C,en,pcrush,evcrush,plock,evlock,tc,fc,
     * k1,k2,k3,SFmax,EFmin,Eodot,dShearModu,d2G,alamda,dBulkModu,
     * d3bulk,ddamage,pn,kav,dmu,ptrial,c1,dmup,pt,dT,mua,ev3,D1,D2,
     * varj2t,qstar,edoteq,c2,sigmay,pstar,phi,kfactor,dgama
      double precision strial,ssa,etrie,ssb,failure
c
c     stateNew(k,1)=总应变
c          stateNew(k,2)
c          stateNew(k,3)
c          stateNew(k,4)
c          stateNew(k,5)
c          stateNew(k,6)
c          stateNew(k,7)=弹性应变
c          stateNew(k,8)
c          stateNew(k,9)
c          stateNew(k,10)
c          stateNew(k,11)
c          stateNew(k,12)
c          stateNew(k,13)=ddamage
c     stateNew(k,14)=dgama
c
      iii=ndir+nshr
      e=props(1)
      xnu=props(2)
      A=props(3)
      B=props(4)
      C=props(5)
      en=props(6)
      pcrush=props(7)
      evcrush=props(8)
      plock=props(9)
      evlock=props(10)
      tc=props(11)
      fc=props(12)
      D1=props(13)
      D2=props(14)
          k1=props(15)
          k2=props(16)
          k3=props(17)
          SFmax=props(18)
          EFmin=props(19)
          Eodot=props(20)
C
          dShearModu=e/(two*(one+xnu))
          d2G=two*ShearModu
          alamda=xnu*d2G/(one-two*xnu)
          dBulkModu=alamda+twothird*dShearModu
          d3bulk=three*dBulkModu
c
          do k=1,nblock
c         
          write(7,110)'k', k
110       format(i6,/)
          ddamage=stateOld(k,13)
          pn=-third*(stressOld(k,1)+stressOld(k,2)+stressOld(k,3))
          write(7,120)'pn', pn
120       format(f10.7,/)         
c        !n平均应力Pn,压为+
          kav=(one-ddamage)*dBulkModu+ddamage*k1
          dmu=-(strainInc(k,1) + strainInc(k,2) + strainInc(k,3))
          write(7,150) 'dmu',dmu
150       format(f10.7,/)
c        ! 体积应变增量,压为+
          ptrial=pn+kav*dmu
          dT=fc*(one-ddamage)
          if((ptrial+dT).lt.zero) then
              pt=-dT      !拉为-,压为+
              dmup=zero
          else if((ptrial .lt. pn) .or. (ptrial .le. pcrush)) then  
c                !卸载 or 第一阶段
              pt=ptrial
              dmup=zero
          else if(ptrial .lt. plock) then    ! 非卸载 第二阶段
              c1=(plock-pcrush)/evlock     ! H=(Pl-Pc)/Ul
              dmup=(ptrial-pn)/(kav+c1)    !塑性体积应变增量
              pt=pn+kav*(dmu-dmup)
          else    ! 第三阶段
              mua=dmu-(stateOld(k,1)+stateOld(k,2)+stateOld(k,3))
c                         !n+1总的体积应变,压为+
              c1=(mua-evlock)/(one+evlock)
              pt=c1*(k1+c1*(k2+c1*k3))
              dmup=zero
          endif
c   体积应变计算结束,得到n+1时刻的平均应力pt,塑性体积应变增量dmup           
c   计算弹性体积应变/3-->ev3,符号仍是压为+
          ev3=pt/d3bulk
c   计算n+1时刻的偏应力.首先计算反对称应变张量产生的偏应力strial(6)
c   n时刻的偏应力ssa(6)
          ssa(1)=stressOld(k,1)+pn
          ssa(2)=stressOld(k,2)+pn
          ssa(3)=stressOld(k,3)+pn
          ssa(4)=stressOld(k,4)
c          if(nshr.gt.one)then
          ssa(5)=stressOld(k,5)
          ssa(6)=stressOld(k,6)
c          endif
c   Sx=2(S12*dW12-S13*dW31)  Sy=2(-S12*dW12+S23*dW23)  Sz=2(-S23*dW23+S31*dW31)
          strial(1)=two*(-relSpinInc(k,3)*ssa(4)
     *+relSpinInc(k,2)*ssa(6))
          strial(2)=two*(relSpinInc(k,3)*ssa(4)
     *-relSpinInc(k,1)*ssa(5))
          strial(3)=two*(relSpinInc(k,1)*ssa(5)
     *-relSpinInc(k,2)*ssa(6))
c   Sxy=(S22-S11)dW12+S31*dW23-S23*dW31
          strial(4)=(ssa(1)-ssa(2))*relSpinInc(k,3)
     1 -ssa(6)*relSpinInc(k,1)+ssa(5)*relSpinInc(k,2)  
c   Syz=-S31*dW12+（S33-S22)*dW23+S12*dW31
          strial(5)=ssa(6)*relSpinInc(k,3)
     1 +(ssa(2)-ssa(3))*relSpinInc(k,1)-ssa(4)*relSpinInc(k,2)
c  Szx=S23*dW12-S12*dw23+(S11-S33)* dW31
          strial(6)=-ssa(5)*relSpinInc(k,3)
     1 +ssa(4)*relSpinInc(k,1)-(ssa(3)-ssa(1))*relSpinInc(k,2)
c   记录Trial弹性偏应变,后面仍用到.（n时刻弹性应变+应变增量）
          c1=third*(stateOld(k,7)+stateOld(k,8)+stateOld(k,9)-dmu)
          etrie(1)=stateOld(k,7)+strainInc(k,1)-c1
          etrie(2)=stateOld(k,8)+strainInc(k,2)-c1
          etrie(3)=stateOld(k,9)+strainInc(k,3)-c1
          etrie(4)=(stateOld(k,10)+strainInc(k,4))/two
          etrie(5)=(stateOld(k,11)+strainInc(k,5))/two
          etrie(6)=(stateOld(k,12)+strainInc(k,6))/two
c  采用弹性Trial偏应变etrie()预估Trial偏应力
          do i=1,iii
              strial(i)=strial(i)+d2G*etrie(i)
          end do
          varj2t=etrie(4)**2+etrie(5)**2+etrie(6)**2
     1+half*(etrie(1)**2+etrie(2)**2+etrie(3)**2)
          qstar=sqrt(three*varj2t)/fc
      write(6,180) 'qstar',qstar
180       format(f10.7,/)
c    计算等效偏应变率
          c1=dmu/three
          ssb(1)=strainInc(k,1)+c1
          ssb(2)=strainInc(k,2)+c1
          ssb(3)=strainInc(k,3)+c1
          ssb(4)=strainInc(k,4)
          ssb(5)=strainInc(k,5)
          ssb(6)=strainInc(k,6)
          if(dt .gt. zero) then
          edoteq=sqrt(two*third*(ssb(1)**2+ssb(2)**2+ssb(3)**2
     1    +two*(ssb(4)**2+ssb(5)**2+ssb(6) **2)))/(dt*Eodot)
          else
              edoteq=zero
          endif
          write(7,*)"strainrate",edoteq
          pstar=pt/fc
          if(pstar .gt. zero) then
              sigmay=(A*(one-ddamage)+B*(pstar**en))*(one+C*log(edoteq))
      write(6,190) 'sigmay',sigmay
190       format(f10.7,/)              
              if((sigmay-SFmax).gt.zero)then
                  sigmay=SFmax
              end if
                else
                    c1=A*(one-ddamage)*(one+C*log(edoteq))
                    c2=tc/fc*(one-ddamage)
              if(abs(c2) .lt. elimit) then
                  sigmay=zero
              else
                        sigmay=c1/c2*abs(c2-abs(pstar))
              end if
          end if
          phi=qstar-sigmay
                if(phi.gt.zero)then
c    塑性屈服
              if(qstar .gt. elimit) then
                  kfactor=sigmay/qstar
                        else
                            kfactor=zero
              end if
              dgama=phi/(three*dShearModu)  !等效塑性应变增量
c              stateNew(k,14)=stateOld(k,14)+dgama
c  计算屈服后的应力  
              stressNew(k,1)=kfactor*strial(1)-pt
                  stressNew(k,2)=kfactor*strial(2)-pt
                  stressNew(k,3)=kfactor*strial(3)-pt
                  stressNew(k,4)=kfactor*strial(4)
              stressNew(k,5)=kfactor*strial(5)
              stressNew(k,6)=kfactor*strial(6)
c   计算弹性应变
              kfactor=kfactor/d2G
                  stateNew(k,7)=kfactor*strial(1)-eve3   !弹性体积应变更新
                  stateNew(k,8)=kfactor*strial(2)-eve3
                  stateNew(k,9)=kfactor*strial(3)-eve3
                  stateNew(k,10)=two*kfactor*strial(4)
              stateNew(k,11)=two*kfactor*strial(5)
              stateNew(k,12)=two*kfactor*strial(6)
          else
c   偏应力为弹性,更新应力
              stressNew(k,1)=strial(1)-pt
                  stressNew(k,2)=strial(2)-pt
                  stressNew(k,3)=strial(3)-pt
                  stressNew(k,4)=strial(4)
              stressNew(k,5)=strial(5)
              stressNew(k,6)=strial(6)
c   偏应力为弹性,更新弹性应变:偏应变为弹性,但体积应变不见得是弹性的
              stateNew(k,7)=etrie(1)-eve3     !弹性体积应变更新
                  stateNew(k,8)=etrie(2)-eve3
                  stateNew(k,9)=etrie(3)-eve3
                  stateNew(k,10)=two*etrie(4)
              stateNew(k,11)=two*etrie(5)
              stateNew(k,12)=two*etrie(6)
          end if
c  结束弹塑性判断,更新总应变
          do i=1,iii
              stateNew(k,i)=stateOld(k,i)+strainInc(k,i)
          end do
c   更新损伤变量   
          c1=D1*((pt+tc)/fc)**D2
          if(c1.lt.EFmin) then
                    c1=EFmin
          end if
          ddamage=(dgama+dmup)/c1
          if(ddamage .gt. one) then
              ddamage=one
          end if
          write(7,*)"D",ddamage
          stateNew(k,13)=ddamage
c
c         
          if((stateNew(k,13)).lt.(0.d85))then
              failure=one
              else
              failure=zero
              end if
          stateNew(k,14)=failure
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
c      plasticWorkInc = half * ( yieldOld + yieldNew ) * dgama
c      enerInelasNew(k) = enerInelasOld(k)
c     * + plasticWorkInc / density(k)
      end do
C  
      return
      end

土木嘉年华

自己顶一个。。。。。。。。。。

土木嘉年华

顶一个、、、、、、、、、、、、、、、、、、、

douxidaa

inp文件呢？

yxlnbu

你好，请问你有相关参考资料么？

孙浩然jun

这是什么材料的本构啊啊

hazard669

你好，请问一下你调试的方法是怎么样呢？

renyoi

可以交流一下吗，我也是做的hjc，计算不收敛，不知道哪不对，可以帮忙看看嘛

努力学仿真

有人会吗

努力学仿真

可以交流一下吗