一、ABAQUS自帶Tresca本構與VUMAT對比

二、Tresca本構模型介紹

   以下， 粗體符號表示向量或矩陣，上標“T”表示向量轉置。

    當屈服函數f（σ）的值為零時，材料屈服。這里σ是應力張量（為列矩陣）。如采用相關聯的流動法則，則無窮小的塑性應變增量為

    式中，a是屈服函數的梯度，dλ是一個非負塑性乘子（non-negative plastic multiplier）。

    如果在分析增量步j 和 j+1之間施加一個應變增量Δε，則應力根據下式進行更新

    式中，ε_j 是增量步 j 之后的總應變，D^ep是無窮小彈塑性本構矩陣。

    由于D^ep高度依賴于σ，因此上式（4.2）通常取近似解。應變增量dε由彈性增量dε^e和塑性增量dε^p組成，因此：

    根據胡克定律得到彈性應力增量如下：

    式中D是彈性本構矩陣。

    結合式（4.1）與式（4.4）可得：

    綜上可得下圖中的徑向映射法則（radial return mapping）：

    上圖中，在增量開始時，考慮點處的應力為σ^A，其在，彈性區域（f<0，f 是屈服函數），當然其同樣可以位于屈服面上（f=0）。彈性預測增量為Δσ^e，由此預測應力為σ^B。通過塑性修正增量?Δσ^p將應力返回到屈服面上的最終應力點σ^C。其中Δσ^p由下式得出：

    通常，Δσ^p由下式進行數值計算：

    其中，下標P表示積分路徑上要計算a的點。在后向歐拉算法中，由于P對應于未知的最終應力狀態σ^C，因此需要迭代。對于主應力呈線性的Tresca屈服準則，a沿積分路徑是常數。

    以下介紹Tresca屈服準則在主應力空間中的表示。屈服面采用與線σ₁=σ₂=σ₃對齊的六角棱鏡的形式，可由六個屈服函數定義：

    對于i=1,…,6，方程f_i（σ，k）=0對應于主應力空間中的平面S_i ,這六個屈服面如下圖所示。

三、Tresca VUMAT 算法實現

輸入：

輸出：

1：

2: 轉化至主應力空間

3:

4:

6：

四、Tresca VUMAT 源代碼（帶注釋）

C  ABAQUS 自帶變量

      subroutine vumat(

C  只讀-

     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  只寫-

     7 stressNew,stateNew,enerInternNew,enerInelasNew)

C

      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)

C

      character*80 cmname

*        !

*        ! ************************** 用戶自定義變量 ************* *************

*        !

      integer :: block_count

      dimension DElas (4 ,4) , test(4,4),

     1  stressNewInt (4 ,1) ,  

     2  strainIncInt (4 ,1) ,  

     3  strPr (3 ,1) ,  

     4  a1 (3 ,1) , a2 (3 ,2) ,  

     5  dStressPlas (3 ,1),

     6  stressinc(4,1)

*        ! ************************************************** ****************

C        ! ** 檢查ndir和nshr是否符合2D或axisymmetry（ndir =3 and nshr =1）

      if (( ndir.ne.3) . or . ( nshr . ne . 1)) then

         call xplb_exit

      end if

*        !

C        ! ** 構建彈性剛度矩陣

*        !

C     PROPS(1) - E ,為彈性模量

C     PROPS(2) - v ,為泊松比

C     PROPS(3) - c, 極限剪應力

      E = props (1)

      v = props (2)

      c=props(3)

      zero=0.0D0

      second=0.5D0

      One=1.0D0

*        !

      tempOne = 1.0 D0 - v

      tempTwo = 1.0 D0 - 2.0 D0 *v

*        !

      DElas (1 ,1:4) = (/ tempOne , v , v , zero /)

      DElas (2 ,1:4) = (/ v , tempOne , v , zero /)

      DElas (3 ,1:4) = (/ v , v , tempOne , zero /)

      DElas (4 ,1:4) = (/ zero , zero , zero ,second * tempTwo/)

*        !

  DElas = ( E /((1.0 D0 + v )*( tempTwo )) ) * DElas

*        !

C        ! ************ 開始主循環 ************

*        !

      do 100 block_count = 1, nblock

*        !

C        ! 使用工程剪應變

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

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

      

*        ! Elastic predictor stress 彈性預測應力

      stressinc=matmul(DElas,strainIncInt)

      stressNewInt(1:4,1)=stressOld(block_count,1:4)+stressinc(1:4,1)

*        !

*        ! Principal values of elastic trial stress  彈性預測主應力

        stressMean = 0.5 D0 * ( stressNewInt (1 ,1) +  

     1   stressNewInt (2 , 1) )

        RMohr = sqrt ((0.5 D0 *( stressNewInt (1 ,1) -  

     1   stressNewInt (2 ,1)))**2  

     2   + ( stressNewInt (4 ,1))**2)

*        !

        strPr (1 ,1) = stressMean + RMohr

        strPr (2 ,1) = stressMean - RMohr

        strPr (3 ,1) = stressNewInt (3 ,1)

        twoTheta = atan2 ( stressNewInt (4 ,1) , 0.5 D0  

     1   *( stressNewIn t (1 ,1)  - stressNewInt (2 ,1)) )

*        !

*        ! Evaluate yield functions 計算屈服函數

        tempOne = 2.0 D0 * c

        f1 = strPr (1 ,1) - strPr (2 ,1) - tempOne

        f2 = strPr (2 ,1) - strPr (1 ,1) - tempOne

        f3 = strPr (2 ,1) - strPr (3 ,1) - tempOne

        f4 = strPr (3 ,1) - strPr (2 ,1) - tempOne

        f5 = strPr (3 ,1) - strPr (1 ,1) - tempOne

        f6 = strPr (1 ,1) - strPr (3 ,1) - tempOne

*        !

*        ! 根據屈服函數值進行修正

  if ( ( max ( f1 ,f2 , f3 , f4 ,f5 , f6 ) <= 0.0 D0 ) . OR .  

     1   ( totalTime == 0.0 D0 )) then

*        ! ** ELASTIC 彈性

*        ! 無需修正            

*        !

      else if ( ( f1 >= 2.0 D0 * f4 ).AND.( f1 >= 2.0 D0 * f6 ) ) then

*        ! ** Return to single surface (1) 修正至面（1）

      a1 (1:3 ,1) = (/ -0.5 D0 , 0.5 D0 , 0.0 D0 /)

      strPr = strPr + a1 * f1

*        !

      else if ( ( f4 >= 2.0 D0 * f1 ).AND.( f4 >= 2.0 D0 * f5 ) ) then

*        ! ** Return to single surface (4) 修正至面（4）

      a1 (1:3 ,1) = (/ 0.0 D0 , 0.5 D0 , -0.5 D0 /)

      strPr = strPr + a1 * f4

*        !

      else if ( ( f6 >= 2.0 D0 * f1 ).AND.( f6 >= 2.0 D0 * f3 ) ) then

*        ! ** Return to single surface (6) 修正至面（6）

      a1 (1:3 ,1) = (/ -0.5 D0 , 0.0 D0 , 0.5 D0 /)

      strPr = strPr + a1 * f6

*        !

      else if ( f6 < 2.0 D0 * f3 ) then

*        ! ** Return to line (6/3) 修正至線（6/3）

      a2 (1:3 ,1) =(/ 1.0 D0 /3.0 D0 ,-2.0 D0 /3.0D0,1.0 D0 /3.0 D0 /)

      a2 (1:3 ,2) = (/-2.0 D0 /3.0 D0,1.0 D0 /3.0 D0,1.0 D0 /3.0 D0 /)

      strPr = strPr +matmul ( a2 , reshape ( (/ f3 , f6 /) ,(/2,1 /)))

*        !

      else if (( f1 < 2.0 D0 * f6 ) . AND . ( f6 < 2.0 D0 * f1)) then

*        ! ** Return to line (1/6) 修正至線（1/6）

      a2 (1:3 ,1)=(/-1.0 D0 /3.0 D0,2.0 D0/3.0 D0 ,-1.0 D0 /3.0 D0/)

      a2 (1:3 ,2)=(/ -1.0 D0/3.0 D0,-1.0 D0/3.0 D0 ,2.0 D0 /3.0 D0 /)

      strPr = strPr + matmul ( a2 ,reshape ( (/ f1, f6/) ,(/ 2,1/) ) )

*        !

      else if ( ( f1 < 2.0 D0 * f4 ) . AND . ( f4 < 2.0 D0 * f1)) then

*        ! ** Return to line (4/1) 修正至線（4/1）

      a2 (1:3 ,1) = (/-2.0 D0 /3.0 D0 ,1.0 D0 /3.0 D0,1.0 D0 /3.0D0 /)

      a2 (1:3 ,2) = (/ 1.0 D0 /3.0 D0 , 1.0 D0 /3.0 D0,-2.0D0/3.0D0 /)

      strPr = strPr + matmul ( a2 , reshape ( (/ f1 ,f4/) ,(/2,1/)))

*        !

      else if ( f4 < 2.0 D0 * f5 ) then

*        ! ** Return to line (5/4) 修正至線（5/4）

      a2 (1:3 ,1) = (/ -1.0 D0 /3.0D0,2.0 D0/3.0 D0 ,1.0 D0 /3.0 D0 /)

      a2 (1:3 ,2) = (/2.0 D0 /3.0D0 , -1.0D0 /3.0D0,-1.0 D0 /3.0 D0 /)

      strPr = strPr + matmul(a2,reshape ((/ f4 , f5 /) ,(/ 2, 1 /) ) )

*        !

      else

*        ! ** 錯誤

        call xplb_exit

      end if

*        !

*        ! ** Transform principal stresses back to Cartesian 將主應力轉換回笛卡爾坐標

        twoTheta = - twoTheta

        stressMean = 0.5 D0 *( strPr (1 ,1) + strPr (2 ,1))

        offset = 0.5 D0 *( strPr (1 ,1) - strPr (2 ,1)) * cos (twoTheta)

        stressNewInt (1 ,1) = stressMean + offset

        stressNewInt (2 ,1) = stressMean - offset

        stressNewInt (3 ,1) = strPr (3 ,1)

        stressNewInt (4 ,1) = 0.5 D0 *( strPr (2 ,1) - strPr (1 ,1)) *  

     1  sin( twoTheta )

*        !

*        ! ** Prepare output 輸出，更新應力

        stressNew ( block_count , 1:4 ) = stressNewInt ( 1:4 , 1 )

*        !

100   continue

      end

程序下載：

Tresca-VUMAT-2D-1.zip