本文介紹Abaqus官方教程的第2個案例——非等溫線彈性UMAT子程序。本文首先簡要介紹非等溫線彈性理論，再采用批注的方式介紹UMAT子程序的實現方法。最后將計算結果與Abaqus自帶材料的進行對比，探討Abaqus熱力耦合的實現方式。筆者水平有限，若有不足之處，煩請指出，不甚感激。

將拉梅常數看做是溫度的函數，則非等溫彈性張量方程為：

其中，l和m是拉梅常數，a是線膨脹系數，T是溫度。對方程兩邊取材料時間導數，得到Jaumann率形式：

寫成增量形式為：

Abaqus Standard采用增量法逐步施加載荷/位移，每步增加的應力即可按上式進行計算。材料的切線剛度矩陣（雅克比矩陣）為應力對應變的偏導數，采用voigt記法寫成6*6的二維矩陣為（引用知乎@Learning STRB）：

簡單介紹了一下理論，下面是UMAT程序：

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,

     1 RPL,DDSDDT,DRPLDE,DRPLDT,

     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,

     3 NDI,NSHR,NTENS,NSTATEV,PROPS,NPROPS,COORDS,DROT,PNEWDT,

     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)

C

      INCLUDE 'ABA_PARAM.INC'

C

      CHARACTER*8 CMNAME

      DIMENSION STRESS(NTENS),STATEV(NSTATEV),

     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),

     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),

     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)

C     LOCAL ARRAYS

C     ----------------------------------------------------------------

C     EELAS - ELASTIC STRAINS 彈性應變（不含熱應變）

C     ETHERM - THERMAL STRAINS 熱應變

C     DTHERM - INCREMENTAL THERMAL STRAINS 增量熱應變

C     DELDSE - CHANGE IN STIFFNESS DUE TO TEMPERATURE CHANGE 雅克比矩陣增量

C     ----------------------------------------------------------------

      DIMENSION EELAS(6), ETHERM(6), DTHERM(6), DELDSE(6,6)

C

      PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, SIX=6.D0)

C     ----------------------------------------------------------------

C     UMAT FOR ISOTROPIC THERMO-ELASTICITY WITH LINEARLY VARYING

C     MODULI - CANNOT BE USED FOR PLANE STRESS

C     ----------------------------------------------------------------

C     PROPS(1) - E(T0)

C     PROPS(2) - NU(T0)

C     PROPS(3) - T0

C     T0溫度下的彈性模量和泊松比

C     PROPS(4) - E(T1)

C     PROPS(5) - NU(T1)

C     PROPS(6) - T1

C     T1溫度下的彈性模量和泊松比

C     PROPS(7) - ALPHA

C     PROPS(8) - T_INITIAL

C     線膨脹系數和參考溫度

C     ELASTIC PROPERTIES AT START OF INCREMENT

C     技術鄰@snowwave02 大神的方法，配合他自研的isolver求解器可以實現斷點調試，這一步可以定位到第*次進入UMAT程序。

      integer,save::number=0

      number=number+1

      if(number.eq.3) then

         number1=0

      endif

C     初始溫度下的拉梅常數

      FAC1=(TEMP-PROPS(3))/(PROPS(6)-PROPS(3))

      IF (FAC1 .LT. ZERO) FAC1=ZERO

      IF (FAC1 .GT. ONE) FAC1=ONE

      FAC0=ONE-FAC1

      EMOD=FAC0*PROPS(1)+FAC1*PROPS(4)

      ENU=FAC0*PROPS(2)+FAC1*PROPS(5)

      EBULK3=EMOD/(ONE-TWO*ENU)

      EG20=EMOD/(ONE+ENU)

      EG0=EG20/TWO

      ELAM0=(EBULK3-EG20)/THREE

C

C     ELASTIC PROPERTIES AT END OF INCREMENT

C     增量溫度下的拉梅常數

      FAC1=(TEMP+DTEMP-PROPS(3))/(PROPS(6)-PROPS(3))

      IF (FAC1 .LT. ZERO) FAC1=ZERO

      IF (FAC1 .GT. ONE) FAC1=ONE

      FAC0=ONE-FAC1

      EMOD=FAC0*PROPS(1)+FAC1*PROPS(4)

      ENU=FAC0*PROPS(2)+FAC1*PROPS(5)

      EBULK3=EMOD/(ONE-TWO*ENU)

      EG2=EMOD/(ONE+ENU)

      EG=EG2/TWO

      ELAM=(EBULK3-EG2)/THREE

     

C     ELASTIC STIFFNESS AT END OF INCREMENT AND STIFFNESS CHANGE

C     求解雅克比矩陣和雅克比矩陣的增量（簡單說是T1溫度和T0溫度的雅克比矩陣差值）

      DO K1=1,NDI

      DO K2=1,NDI

      DDSDDE(K2,K1)=ELAM

      DELDSE(K2,K1)=ELAM-ELAM0

      END DO

      DDSDDE(K1,K1)=EG2+ELAM

      DELDSE(K1,K1)=EG2+ELAM-EG20-ELAM0

      END DO

      DO K1=NDI+1,NTENS

      DDSDDE(K1,K1)=EG

      DELDSE(K1,K1)=EG-EG0

      END DO

C

C     CALCULATE THERMAL EXPANSION

C     計算熱膨脹應變和熱膨脹應變的增量（熱應變僅有主應變，沒有切應變）

      DO K1=1,NDI

      ETHERM(K1)=PROPS(7)*(TEMP-PROPS(8))

      DTHERM(K1)=PROPS(7)*DTEMP

      END DO

       DO K1=NDI+1,NTENS

      ETHERM(K1)=ZERO

      DTHERM(K1)=ZERO

      END DO

C

C     CALCULATE STRESS, ELASTIC STRAIN AND THERMAL STRAIN

C     按應力的增量公式計算應力增量，進而計算應力。

      DO K1=1, NTENS

      DO K2=1, NTENS

      STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*(DSTRAN(K1)-DTHERM(K1))

     1 +DELDSE(K2,K1)*( STRAN(K1)-ETHERM(K1))

      END DO

      ETHERM(K1)=ETHERM(K1)+DTHERM(K1)

      EELAS(K1)=STRAN(K1)+DSTRAN(K1)-ETHERM(K1)

      END DO

C

C     STORE ELASTIC AND THERMAL STRAINS IN STATE VARIABLE ARRAY

C     存儲彈性應變至狀態變量SDV1-SDV12

      DO K1=1, NTENS

      STATEV(K1)=EELAS(K1)

      STATEV(K1+NTENS)=ETHERM(K1)

      END DO

      RETURN

      END

建立單個單元的熱力耦合有限元模型。

有限元材料參數為（Model-1）：E(T0)=200GPa，NU(T0)=0.3，T0=293K；E(T1)=300GPa，NU(T1)=0.3，T1=393K；ALPHA=1e-5，T_INITIAL=293K。另外設置密度、比熱容、熱導率。為了排除傳熱因素的影響，我在后續載荷步里采用直接指定溫度的方法，因此上述參數可任意設置。

作為對照，采用Abaqus自帶的材料設置Model-2：

幾何模型為1m*1m*1m的立方體，采用C3D8T繪制1個單元。設置兩個分析步：初始狀態設置域溫度為293K；分析步1約束X負方向上四個節點X方向位移，X正方向上四個節點加載X方向位移0.1m；分析步2設置整個域溫度為393K。

經過對比，Model-1和Model-2的計算結果完全一致。下面僅展示Model-1的計算結果：前兩個圖是分析步1的X方向應力和總應變，后圖為分析步2的X方向應力和彈性應變（采用位移加載，分析步2的總應變不變）。

下面對計算結果進行討論。筆者經過有限元計算，如果不考慮熱膨脹，分析步2計算完成后應力為3e10Pa。由于自由熱膨脹不產生應力，因此線膨脹系數的引入造成了應力的減小。

特別注意的是，在有限元建模過程中，自定義的材料并未使用Abaqus自帶的熱膨脹項，完全是通過編寫的UMAT來計算材料的熱膨脹貢獻。那么，如果我們額外增加Abaqus自帶的熱膨脹項（材料屬性里的expansion項），那么彈性應變會變為0.098，應力為變為2.94e10MPa。也就是說我們計算了兩次熱膨脹系數。筆者做出了如下猜想：一旦設置了熱膨脹屬性，進入Abaqus的應變就從總應變張量變為了彈性應變張量。因此，如果需要編寫含溫度參數的UMAT子程序，我們要么使用Abaqus自帶的expansion項，要么在UMAT中將熱膨脹貢獻考慮，切勿重復。

上述討論僅針對了應變和應變增量，許多本構模型采用的是變形梯度描述的，后續有時間筆者會對此再做一次驗證。

線彈性變溫本構UMAT.pdf