Tresca本构模型VUMAT-3D（源代码免费）

二维Tresca本构模型VUMAT，请参考（完全免费）：https://www.jishulink.com/post/1867025

一、ABAQUS自带Tresca本构与VUMAT-3D对比

图1 模型示意图

图2  ABAQUS 自带计算应力云图

图3  VUMAT-Tresca-3D 计算应力云图

图4  Tresca应力对比（Path1）

图5  S11应力对比（Path1）

表5 Job监控数据对比

二、Tresca本构模型介绍和算法实现

本人已在二维模型中写明，在此不作赘述，请参考：https://www.jishulink.com/post/1867025

三、Tresca-3 VUMAT 源代码（完整）

      module mymod 

      integer,parameter::Q=8

      real(Q),parameter::zero=0.0D0,one=1.D0,two=2.0D0,three=3.0D0, 

     &                         four=4.0D0,ref=3.D-6,half=0.5D0        

      contains

      subroutine ElasticMatrix(E,v,M)

           real(Q),intent(out),dimension(6,6)::M

           real(Q),intent(in)::E,v

           M(1,1:6)=(/one-v,v,v,zero,zero,zero/)

           M(2,1:6)=(/v,one-v,v,zero,zero,zero/)

           M(3,1:6)=(/v,v,one-v,zero,zero,zero/)

           M(4,1:6)=(/zero,zero,zero,half*(one-two*v),zero,zero/)

           M(5,1:6)=(/zero,zero,zero,zero,half*(one-two*v),zero/)

           M(6,1:6)=(/zero,zero,zero,zero,zero,half*(one-two*v)/)

           M=M*E/((one+v)*(one-two*v))

      end subroutine ElasticMatrix

      subroutine Matrix_inverse(M,Inve)

           real(Q),intent(in),dimension(3,3)::M

           real(Q),intent(out),dimension(3,3)::Inve

           real(Q),dimension(3,3)::temp

           integer::i,j

           real(Q)::Det

           temp(1,1)=M(2,2)*M(3,3)-M(2,3)*M(3,2)

           temp(2,1)=M(2,3)*M(3,1)-M(2,1)*M(3,3)

           temp(3,1)=M(2,1)*M(3,2)-M(2,2)*M(3,1)

           temp(3,2)=M(1,2)*M(3,1)-M(1,1)*M(3,2)

           temp(2,2)=M(1,1)*M(3,3)-M(1,3)*M(3,1)

           temp(3,3)=M(1,1)*M(2,2)-M(1,2)*M(2,1)

           temp(1,2)=M(1,3)*M(3,2)-M(1,2)*M(3,3)

           temp(1,3)=M(1,2)*M(2,3)-M(1,3)*M(2,2)

           temp(2,3)=M(1,3)*M(2,1)-M(1,1)*M(2,3)

           Det=M(1,1)*temp(1,1)+M(1,2)*temp(2,1)+M(1,3)*temp(3,1)

           Inve=temp/Det

      end subroutine Matrix_inverse      

      subroutine returntoplane(p,s)

           real(Q),intent(inout)::p

           real(Q),intent(in)::s

           real(Q)::a,b,c,d,t

           dimension::p(3),s(4)

           t=-(s(1)*p(1)+s(2)*p(2)+s(3)*p(3)+s(4))/two!!!!(a**2+b**2+c**2)

           p(1)=s(1)*t+p(1)

           p(2)=s(2)*t+p(2)

           p(3)=s(3)*t+p(3)

      end subroutine returntoplane

      subroutine returntoline(p,l)

           real(Q),intent(inout)::p

           real(Q),intent(in)::l

           real(Q)::t

           dimension::p(3),l(3)

           t=(p(1)+p(2)+p(3)-l(1)-l(2)-l(3))/three

           p(1:3)=(/t+l(1),t+l(2),t+l(3)/)

      end subroutine returntoline      

      

      

      subroutine  modi_pstress(s,c)

           integer::rp

           real(Q),intent(inout),dimension(3)::s

           real(Q),intent(in)::c

           real(Q)::plane,line,f1,f2,f3,f4,f5,f6,two3,four3

           dimension ::plane(4),line(6)

           two3=c*two/three

           four3=c*four/three

           f1=s(1)-s(2)-c*2.0D0

           f2=s(2)-s(1)-c*2.0D0

           f3=s(2)-s(3)-c*2.0D0

           f4=s(3)-s(2)-c*2.0D0

           f5=s(3)-s(1)-c*2.0D0

           f6=s(1)-s(3)-c*2.0D0

           if ((f1>=f4*2.D0).and.(f1>=f6*2.D0).and.(f1>=0.D0)) then

                plane(1:4)=(/one,-one,zero,-c*two/)

                call returntoplane(s,plane)

           else if((f2>=f3*2.D0).and.(f2>=f5*2.D0).and.(f2>=0.D0))then

                plane(1:4)=(/-one,one,zero,-c*two/)

                call returntoplane(s,plane)

           else if((f3>=f2*2.D0).and.(f3>=f6*2.D0).and.(f3>=0.D0))then

                plane(1:4)=(/zero,one,-one,-c*two/)

                call returntoplane(s,plane)

           else if((f4>=f1*2.D0).and.(f4>=f5*2.D0).and.(f4>=0.D0))then

                plane(1:4)=(/zero,-one,one,-c*two/)

                call returntoplane(s,plane)

           else if((f5>=f2*2.D0).and.(f5>=f4*2.D0).and.(f5>=0.D0))then

                plane(1:4)=(/-one,zero,one,-c*two/)

                call returntoplane(s,plane)

           else if((f6>=f3*2.D0).and.(f6>=f1*2.D0).and.(f6>=0.D0))then

                plane(1:4)=(/one,zero,-one,-c*two/)

                call returntoplane(s,plane)

           else if((f1<f4*2.D0).and.(f4<f1*2.D0).and.(f1>=0.D0))then

                line(1:3)=(/two3,-four3,two3/)

                call returntoline(s,line)

           else if((f4<f5*2.D0).and.(f5<f4*2.D0).and.(f4>=0.D0))then

                line(1:3)=(/-two3,-two3,four3/)

                call returntoline(s,line)           

           else if((f2<f5*2.D0).and.(f5<f2*2.D0).and.(f2>=0.D0))then

                line(1:3)=(/-four3,two3,two3/)

                call returntoline(s,line)

           else if((f2<f3*2.D0).and.(f3<f2*2.D0).and.(f3>=0.D0))then

                line(1:3)=(/-two3,four3,-two3/)

                call returntoline(s,line)

          else if((f3<f6*2.D0).and.(f6<f3*2.D0).and.(f6>=0.D0))then

                line(1:3)=(/two3,two3,-four3/)

                call returntoline(s,line)

          else if((f6<f1*2.D0).and.(f1<f6*2.D0).and.(f6>=0.D0))then

                line(1:3)=(/four3,-two3,-two3/)

                call returntoline(s,line)

          end if

      end subroutine modi_pstress

      

      subroutine getnewstress(Ematrix,Pstress,newstress)

           real(Q),intent(out),dimension(6)::newstress           

           real(Q)::Pstress,Ematrix,inve,temp

           dimension::Ematrix(3,3),Pstress(3),inve(3,3),temp(3,3)

            call Matrix_inverse(Ematrix,Inve)   

            

            Ematrix(1:3,1)=Ematrix(1:3,1)*Pstress(1)

            Ematrix(1:3,2)=Ematrix(1:3,2)*Pstress(2)

            Ematrix(1:3,3)=Ematrix(1:3,3)*Pstress(3)

            temp=matmul(Ematrix,Inve)

            newstress(1:3)=(/temp(1,1),temp(2,2),temp(3,3)/)

            newstress(4:6)=(/temp(1,2),temp(2,3),temp(1,3)/)

      end subroutine getnewstress

     

      

          end module mymod

C***********************************************************************      

C***************************  VUMAT  ***********************************

      subroutine vumat(

C  Readonly-

     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  Writeonly-

     7 stressNew,stateNew,enerInternNew,enerInelasNew)

C      

      use mymod

      include 'vaba_param.inc'

      dimension props(nprops),density(nblock),

     1  coordMp(nblock,*),

     2  charLength(*),strainInc(nblock,ndir+nshr),

     3  relSpinInc(*),tempOld(*),

     4  stretchOld(*),defgradOld(*),

     5  fieldOld(*),stressOld(nblock,ndir+nshr),

     6  stateOld(nblock,nstatev),enerInternOld(nblock),

     7  enerInelasOld(nblock),tempNew(*),

     8  stretchNew(*),defgradNew(*),fieldNew(*),

     9  stressNew(nblock,ndir+nshr),stateNew(nblock,nstatev),

     1  enerInternNew(nblock),enerInelasNew(nblock)

      character*80 cmname

C***********************************************************************

      real(Q)::nam,f1,f2,f3,f4,f5,f6,c,

     &         E,DELAS,dstrain,pstrain,yield,

     &           tstress,pstress,pvec,p,Ematrix,newstress,pvecc

   

      integer::bi,ji,ki,li

C      logical::yield

      dimension::DELAS(6,6),stressNewInt (6,1),stressinc(6,1),

     & tstress(nblock,6),pstress(nblock,3),pvec(nblock,3,3),p(3),

     & pstrain(nblock,3),Ematrix(3,3),strainIncInt(6,1),newstress(6),

     & pvecc(3,3) 

C    get soil properties

      E = props(1)

      nam = props(2)

      c = props(3)

      call ElasticMatrix(E,nam,DELAS)

      

C      if(totaltime.gt.0)then

C      read * !输入回车，继续运行

C      end if   

      do 200 bi=1,nblock

      strainIncInt ( 1:3 , 1 ) = strainInc( bi,1:3 )

      strainIncInt ( 4 , 1 ) = 2.0D0*strainInc( bi,4 )

      strainIncInt ( 5 , 1 ) = 2.0D0*strainInc( bi,5 )

      strainIncInt ( 6 , 1 ) = 2.0D0*strainInc( bi,6 )

      stressinc=matmul(DELAS,strainIncInt)

      do 300 ji=1,6

      tstress(bi,ji)=stressOld(bi,ji)+stressinc(ji,1)          

  300 continue

     

  200 continue

      call vsprind(nblock,tstress,pstress,pvec,3,3)

      

C*****************************  Begin the main loop  **************************************         

        do 100 bi=1,nblock

!to check if yield

        p(1:3)=pstress(bi,1:3)

        f1=p(1)-p(2)-c*2.0D0

        f2=p(2)-p(1)-c*2.0D0

        f3=p(2)-p(3)-c*2.0D0

        f4=p(3)-p(2)-c*2.0D0

        f5=p(3)-p(1)-c*2.0D0

        f6=p(1)-p(3)-c*2.0D0

        if((max(f1,f2,f3,f4,f5,f6)<=0.D0).or.(Totaltime<=1.0D0))then

             stressNew(bi,1:6)=tstress(bi,1:6)

        else

       call modi_pstress(p,c)

        pvecc(1,1)=pvec(bi,1,1)

       pvecc(2,1)=pvec(bi,1,2)

       pvecc(3,1)=pvec(bi,1,3)

       pvecc(1,2)=pvec(bi,2,1)

       pvecc(2,2)=pvec(bi,2,2)

       pvecc(3,2)=pvec(bi,2,3)

       pvecc(1,3)=pvec(bi,3,1)

       pvecc(2,3)=pvec(bi,3,2)

       pvecc(3,3)=pvec(bi,3,3)

       call getnewstress(pvecc,p,newstress)

             do ji=1,6

             stressNew(bi,ji)=newstress(ji)

             end do

       end if

  100 continue

        end

以上为完整的Tresca-3D VUMAT源代码

如需包含注释的版本，请参考付费内容，感谢打赏辛苦费，您的支持是我更新的动力。