ansys-umat.for

ansys或ABAQUS二次材料属性的开发,实现材料非线性.

*****************************************************************top***
**: <Fully implicit integration of the Combined OW-AF model
**:  for 3D/axisymmetric solid and plane strain elements>
************************************************************************
**:
**:  Notice!
**:
**:  Set '*depvar'=ntens*(mah+1)+1 in input file.
**:
************************************************************************
**: <ABAQUS user subroutine, UMAT>
************************************************************************
**:
**: <UMAT heading>
**:
      subroutine umat(ss,statev,ddsdde,sse,spd,scd,
     & rpl,ddsddt,drplde,drpldt,
     & sn,dsn,time,dtime,temp,dtemp,predef,dpred,cmname,
     & ndi,nshr,ntens,nstatv,props,nprops,coords,drot,pnewdt,
     & celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)

      include 'aba_param.inc'

      character*8 cmname
      dimension ss(ntens),statev(nstatv),
     & ddsdde(ntens,ntens),ddsddt(ntens),drplde(ntens),
     & sn(ntens),dsn(ntens),time(2),predef(1),dpred(1),
     & props(nprops),coords(3),drot(3,3),dfgrd0(3,3),dfgrd1(3,3)
**:
************************************************************************
**: <Local variable definitions>
************************************************************************
**:
**: <elastic constants>
      double precision yg,pn,sg
      parameter (yg=1.65d5,pn=0.3d0,sg=0.5d0*yg/(1.0d0+pn))
**:
**: <yield stress>
      double precision sy
      parameter (sy=2.40d2)
**:
**: <kinematic hardening parameters>
      integer mah
      parameter (mah=4)
      double precision zi(mah),ri(mah)
      data zi(1),ri(1)/2.0d3,5.0d1/
      data zi(2),ri(2)/5.0d2,5.0d1/
      data zi(3),ri(3)/2.0d2,4.0d1/
      data zi(4),ri(4)/5.0d1,1.0d2/
      double precision owaf
      parameter (owaf=0.0d0)
      save zi,ri
**:
      double precision em(6,6),vem(6,6)
      double precision snp(6),dsnp(6),dp,pacc
      double precision bhi(mah,6),tah(6),tbhi(mah,6),thi(mah),fbi(mah)
      double precision tvs(6),vsa(6),fy
      double precision pm(6,6),khm(6,6)

      integer ite,itemax
      parameter (itemax=25)
      integer fakn
      double precision odp,xakn(3)

      integer i,j,k
      double precision xx
**:
***************************************************main routine start***
**: <Labeling on state variables>
************************************************************************
**:
      k=1
      pacc=statev(k)
      k=k+1
      do 10 j=1,ntens
         snp(j)=statev(k)
         k=k+1
   10 continue
      do 30 i=1,mah
         do 20 j=1,ntens
            bhi(i,j)=statev(k)
            k=k+1
   20    continue
   30 continue
**:
************************************************************************
**: <Elastic stiffness>
************************************************************************
**:
      call kmkem( em,sg,pn,ndi,nshr )
      call kmkvem( vem,sg,ntens,ndi )
**:
************************************************************************
**: <Elastic predictor>
************************************************************************
**:
**: <trial deviatoric stress>
      do 120 i=1,ntens
         tvs(i)=0.0d0
         do 110 j=1,ntens
            tvs(i)=tvs(i)+vem(i,j)*(sn(j)+dsn(j)-snp(j))
  110    continue
  120 continue
**:
**: <trial back stress>
      do 140 j=1,ntens
         tah(j)=0.0d0
         do 130 i=1,mah
            tah(j)=tah(j)+ri(i)*bhi(i,j)
  130    continue
  140 continue
**:
************************************************************************
**: <Yield function>
************************************************************************
**:
      xx=0.0d0
      do 210 j=1,ndi
         xx=xx+(tvs(j)-tah(j))**2
  210 continue
      do 220 j=ndi+1,ndi+nshr,1
         xx=xx+0.5d0*(tvs(j)-tah(j))**2
  220 continue
      fy=sqrt( 1.5d0*xx )-1.0001d0*sy
**:
************************************************************************
**: <Plastic corrector>
************************************************************************
**:
      if (fy .gt. 0.0d0) then

         ite=1
         fakn=0
         odp=0.0d0
         do 310 i=1,mah
            thi(i)=1.0d0
  310    continue
**:
**: <return mapping>
 1000    call krmap( dp,dsnp,vsa,tvs,bhi,thi,zi,ri,sg,sy,
     &               mah,ndi,nshr,ntens )
         call kaccel( dp,dsnp,xakn,fakn,ntens )
**:
**: <hardening parameters correction>
         call kkinh( thi,tbhi,fbi,bhi,dsnp,dp,owaf,zi,
     &               mah,ntens,ndi,nshr )
         if ((abs( dp-odp ) .gt. 1.0d-5*dp) 
     &       .and. (ite .lt. itemax)) then
            ite=ite+1
            odp=dp
            goto 1000
         endif

      endif
**:
************************************************************************
**: <internal variable update>
************************************************************************
**:
      if (fy .gt. 0.0d0) then

         do 410 j=1,ntens
            snp(j)=snp(j)+dsnp(j)
  410    continue
         pacc=pacc+dp
         do 430 j=1,ntens
            do 420 i=1,mah
               bhi(i,j)=thi(i)*tbhi(i,j)
  420       continue
  430    continue

      endif
**:
************************************************************************
**: <Stress update>
************************************************************************
**:
      do 520 i=1,ntens
         ss(i)=0.0d0
         do 510 j=1,ntens
            ss(i)=ss(i)+em(i,j)*(sn(j)+dsn(j)-snp(j))
  510    continue
  520 continue
**:
************************************************************************
**: <Tangent operator for global Newton iteration>
************************************************************************
**:
      if (fy .gt. 0.0d0) then
         call kpm( pm,vsa,sy,dp,ntens,ndi,nshr )
         call kkhm( khm,owaf,zi,ri,thi,bhi,fbi,vsa,sy,
     &              ntens,ndi,nshr,mah )
         call kctop( em,pm,khm,vem,ntens )
      endif

      do 620 i=1,ntens
         do 610 j=1,ntens
            ddsdde(i,j)=em(i,j)
  610    continue
  620 continue
**:
************************************************************************
**: <State variables update>
************************************************************************
**:
      k=1
      statev(k)=pacc
      k=k+1
      do 710 j=1,ntens
         statev(k)=snp(j)
         k=k+1
  710 continue
      do 730 i=1,mah
         do 720 j=1,ntens
            statev(k)=bhi(i,j)
            k=k+1
  720    continue
  730 continue
**:
*****************************************************main routine end***
**: <End of UMAT>
************************************************************************
**:
      return
      end

************************************************************************
**: <User-defined local functions and subroutines>
************************************************************************
**:
**: <Return mapping>
**:
      subroutine krmap( dp,dsnp,vsa,tvs,bhi,thi,zi,ri,sg,sy,
     &                  mah,ndi,nshr,ntens )
      include 'aba_param.inc'
      integer mah,ndi,nshr,ntens
      double precision dp,dsnp(6),vsa(6),tvs(6)
      double precision bhi(mah,6),thi(mah),zi(mah),ri(mah)
      double precision sg,sy

      integer i,j
      double precision tvsa(6),ahn(6),tsef,hh

      do 20 j=1,ntens
         ahn(j)=0.0d0
         do 10 i=1,mah
            ahn(j)=ahn(j)+thi(i)*ri(i)*bhi(i,j)
   10    continue
   20 continue

      hh=0.0d0
      do 30 i=1,mah
         hh=hh+thi(i)*ri(i)*zi(i)
   30 continue

      tsef=0.0d0
      do 40 j=1,ndi
         tvsa(j)=tvs(j)-ahn(j)
         tsef=tsef+tvsa(j)**2
   40 continue
      do 50 j=ndi+1,ndi+nshr,1
         tvsa(j)=tvs(j)-ahn(j)
         tsef=tsef+0.5d0*tvsa(j)**2
   50 continue
      tsef=sqrt( 1.5d0*tsef )

      dp=(tsef-sy)/(3.0d0*sg+hh)

      do 60 j=1,ntens
         vsa(j)=sy*tvsa(j)/tsef
         dsnp(j)=1.5d0*dp*tvsa(j)/tsef
   60 continue

      end

**:
**: <Aitken's acceleration>
**:
      subroutine kaccel( dp,dsnp,xakn,fakn,ntens )
      include 'aba_param.inc'
      integer fakn,ntens
      double precision dp,dsnp(6),xakn(3)

      integer j

      xakn(1)=xakn(2)
      xakn(2)=xakn(3)
      xakn(3)=dp

      if (fakn .ne. 3) then
         fakn=fakn+1
      elseif (abs( xakn(3)-xakn(2) ) .gt. 1.0d-10*dp) then
         dp=xakn(3)-(xakn(3)-xakn(2))
     &               /(1.0d0-(xakn(2)-xakn(1))/(xakn(3)-xakn(2)))
         if (dp .lt. 0.0d0) then
            dp=xakn(3)
         else
            do 10 j=1,ntens
               dsnp(j)=dsnp(j)*dp/xakn(3)
   10       continue
            xakn(3)=dp
            fakn=2
         endif
      endif

      end

**:
**: <Kinematic hardening>
**:
      subroutine kkinh( thi,tbhi,fbi,bhi,dsnp,dp,owaf,zi,
     &                  mah,ntens,ndi,nshr )
      include 'aba_param.inc'
      integer mah,ntens,ndi,nshr
      double precision thi(mah),tbhi(mah,6),fbi(mah)
      double precision bhi(mah,6),dsnp(6),dp,owaf,zi(mah)

      integer i,j
      double precision ci,bbari

      do 30 i=1,mah

         bbari=0.0d0
         do 10 j=1,ndi
            tbhi(i,j)=bhi(i,j)+zi(i)*dsnp(j)/1.5d0
            bbari=bbari+tbhi(i,j)**2
   10    continue
         do 20 j=ndi+1,ndi+nshr,1
            tbhi(i,j)=bhi(i,j)+zi(i)*dsnp(j)/1.5d0
            bbari=bbari+0.5d0*tbhi(i,j)**2
   20    continue
         bbari=sqrt( 1.5d0*bbari )

         ci=1.0d0/(1.0d0+owaf*zi(i)*dp)
         fbi(i)=ci*bbari-1.0d0

         if (fbi(i) .ge. 0.0d0) then
            thi(i)=1.0d0/bbari
         else
            thi(i)=ci
         endif

   30 continue

      end

**:
**: <Elastic stiffness matrix>
**:
      subroutine kmkem( em,sg,pn,ndi,nshr )
      include 'aba_param.inc'
      integer ndi,nshr
      double precision em(6,6),sg,pn
 
      integer i,j
      double precision xx

      xx=2.0d0*sg/(1.0d0-2.0d0*pn)
      do 20 i=1,ndi
         em(i,i)=(1.0d0-pn)*xx
         do 10 j=i+1,ndi
            em(i,j)=pn*xx
            em(j,i)=em(i,j)
   10    continue
   20 continue
      do 40 i=ndi+1,ndi+nshr,1
         em(i,i)=sg
         do 30 j=i+1,ndi+nshr
            em(i,j)=0.0d0
            em(j,i)=em(i,j)
   30    continue
   40 continue
      do 60 i=1,ndi
         do 50 j=ndi+1,ndi+nshr,1
            em(i,j)=0.0d0
            em(j,i)=em(i,j)
   50    continue
   60 continue

      end

**:
**: <Deviatoric elastic stiffness matrix>
**:
      subroutine kmkvem( vem,sg,ntens,ndi )
      include 'aba_param.inc'
      integer ntens,ndi
      double precision vem(6,6),sg

      integer j,k
      double precision xx

      xx=2.0d0*sg
      do 20 j=1,ntens
         do 10 k=1,ntens
            vem(j,k)=0.0d0
   10    continue
         vem(j,j)=xx
   20 continue
      xx=2.0d0*sg/3.0d0
      do 40 j=1,ndi
         do 30 k=1,ndi
            vem(j,k)=vem(j,k)-xx
   30    continue
   40 continue

      end

**:
**: <Consistent tangent stiffness matrix>
**:
      subroutine kctop( em,pm,khm,vem,ntens )
      include 'aba_param.inc'
      integer ntens
      double precision em(6,6),pm(6,6),khm(6,6),vem(6,6)

      integer j,k
      double precision lm(6,6),lminv(6,6)
      double precision xm(6,6),ym(6,6)

      do 20 j=1,ntens
         do 10 k=1,ntens
            lm(j,k)=vem(j,k)+pm(j,k)+khm(j,k)
   10    continue
   20 continue

      call kminv( lm,lminv,ntens )

      call kmprod( lminv,vem,ym,ntens )
      call kmprod( em,ym,xm,ntens )

      do 40 j=1,ntens
         do 30 k=1,ntens
            em(j,k)=em(j,k)-xm(j,k)
   30    continue
   40 continue

      end

**:
**: <Kinematic hardening tangent operator>
**:
      subroutine kkhm( khm,owaf,zi,ri,thi,bhi,fbi,vsa,sy,
     &                 ntens,ndi,nshr,mah )
      include 'aba_param.inc'
      integer ntens,ndi,nshr,mah
      double precision khm(6,6),owaf,zi(mah),ri(mah),thi(mah)
      double precision bhi(mah,6),fbi(mah)
      double precision vsa(6),sy

      integer i,j,k
      double precision xx

      xx=0.0d0
      do 10 i=1,mah
         xx=xx+ri(i)*thi(i)*zi(i)/1.5d0
   10 continue
      do 30 j=1,ntens
         do 20 k=1,ntens
            khm(j,k)=0.0d0
   20    continue
         khm(j,j)=xx
   30 continue

      do 300 i=1,mah
         if (fbi(i) .ge. 0.0d0) then
            xx=ri(i)*thi(i)*zi(i)
            do 130 j=1,ntens
               do 110 k=1,ndi
                  khm(j,k)=khm(j,k)-xx*bhi(i,j)*bhi(i,k)
  110          continue
               do 120 k=ndi+1,ndi+nshr,1
                  khm(j,k)=khm(j,k)-0.5d0*xx*bhi(i,j)*bhi(i,k)
  120          continue
  130       continue
         else
            xx=ri(i)*thi(i)*owaf*zi(i)
            do 230 j=1,ntens
               do 210 k=1,ndi
                  khm(j,k)=khm(j,k)-xx*bhi(i,j)*vsa(k)/sy
  210          continue
               do 220 k=ndi+1,ndi+nshr,1
                  khm(j,k)=khm(j,k)-0.5d0*xx*bhi(i,j)*vsa(k)/sy
  220          continue
  230       continue
         endif
  300 continue

      end

**:
**: <Perfectly plastic tangent operator>
**:
      subroutine kpm( pm,vsa,sy,dp,ntens,ndi,nshr )
      include 'aba_param.inc'
      integer ntens,ndi,nshr
      double precision pm(6,6),vsa(6),sy,dp

      integer j,k

      do 30 j=1,ntens
         do 10 k=1,ndi
            pm(j,k)=-sy/dp*vsa(j)/sy*vsa(k)/sy
   10    continue
         do 20 k=ndi+1,ndi+nshr,1
            pm(j,k)=-0.5d0*sy/dp*vsa(j)/sy*vsa(k)/sy
   20    continue
         pm(j,j)=pm(j,j)+sy/(1.5d0*dp)
   30 continue

      end

**:
**: <Inversion of matricies>
**:
      subroutine kminv( am,invm,nsize )
      include 'aba_param.inc'
      integer nsize
      double precision am(6,6),invm(6,6)

      integer i,j,k
      double precision xx,xm(6,6)

      do 20 i=1,nsize
         do 10 j=1,nsize
            xm(i,j)=am(i,j)
            invm(i,j)=0.0d0
   10    continue
         invm(i,i)=1.0d0
   20 continue

      do 60 k=1,nsize-1
         do 50 j=k+1,nsize
            xx=xm(k,j)/xm(k,k)
            do 30 i=k+1,nsize
               xm(i,j)=xm(i,j)-xm(i,k)*xx
   30       continue
            do 40 i=1,nsize
               invm(i,j)=invm(i,j)-invm(i,k)*xx
   40       continue
   50    continue
   60 continue

      do 90 k=1,nsize
         do 80 j=nsize,1,-1
            xx=invm(k,j)
            do 70 i=j+1,nsize
               xx=xx-xm(i,j)*invm(k,i)
   70       continue
            invm(k,j)=xx/xm(j,j)
   80    continue
   90 continue

      end

**:
**: <Matrix production>
**:
      subroutine kmprod( am,bm,cm,ntens )
      include 'aba_param.inc'
      integer ntens
      double precision am(6,6),bm(6,6),cm(6,6)

      integer i,j,k

      do 30 j=1,ntens
         do 20 k=1,ntens
            cm(j,k)=0.0d0
            do 10 i=1,ntens
               cm(j,k)=cm(j,k)+am(j,i)*bm(i,k)
   10       continue
   20    continue
   30 continue

      end

***************************************************************bottom***