非等溫線彈性UMAT子程序及Abaqus內(nèi)部實(shí)現(xiàn)方法

力學(xué)阿東

2021年5月26日 10:25

瀏覽：5970 評(píng)論：21 收藏：47

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

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

其中，l和m是拉梅常數(shù)，a是線膨脹系數(shù)，T是溫度。對(duì)方程兩邊取材料時(shí)間導(dǎo)數(shù)，得到Jaumann率形式：

寫(xiě)成增量形式為：

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

簡(jiǎn)單介紹了一下理論，下面是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)

INCLUDE 'ABA_PARAM.INC'

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 彈性應(yīng)變（不含熱應(yīng)變）

C ETHERM - THERMAL STRAINS 熱應(yīng)變

C DTHERM - INCREMENTAL THERMAL STRAINS 增量熱應(yīng)變

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

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

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

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 線膨脹系數(shù)和參考溫度

C ELASTIC PROPERTIES AT START OF INCREMENT

C 技術(shù)鄰@snowwave02 大神的方法，配合他自研的isolver求解器可以實(shí)現(xiàn)斷點(diǎn)調(diào)試，這一步可以定位到第*次進(jìn)入U(xiǎn)MAT程序。

integer,save::number=0

number=number+1

if(number.eq.3) then

number1=0

endif

C 初始溫度下的拉梅常數(shù)

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 ELASTIC PROPERTIES AT END OF INCREMENT

C 增量溫度下的拉梅常數(shù)

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 求解雅克比矩陣和雅克比矩陣的增量（簡(jiǎn)單說(shuō)是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 CALCULATE THERMAL EXPANSION

C 計(jì)算熱膨脹應(yīng)變和熱膨脹應(yīng)變的增量（熱應(yīng)變僅有主應(yīng)變，沒(méi)有切應(yīng)變）

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 CALCULATE STRESS, ELASTIC STRAIN AND THERMAL STRAIN

C 按應(yīng)力的增量公式計(jì)算應(yīng)力增量，進(jìn)而計(jì)算應(yīng)力。

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 STORE ELASTIC AND THERMAL STRAINS IN STATE VARIABLE ARRAY

C 存儲(chǔ)彈性應(yīng)變至狀態(tài)變量SDV1-SDV12

DO K1=1, NTENS

STATEV(K1)=EELAS(K1)

STATEV(K1+NTENS)=ETHERM(K1)

END DO

RETURN

END

建立單個(gè)單元的熱力耦合有限元模型。

有限元材料參數(shù)為（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。另外設(shè)置密度、比熱容、熱導(dǎo)率。為了排除傳熱因素的影響，我在后續(xù)載荷步里采用直接指定溫度的方法，因此上述參數(shù)可任意設(shè)置。

Q8FF]Q9{HIH%8M~%Y}KQ{Z1.png

作為對(duì)照，采用Abaqus自帶的材料設(shè)置Model-2：

$Q35Y%NAGKM97ZNXWSY{7F_Q.png$ M7675UQABCAT@EAP(9OLU1B.png

幾何模型為1m*1m*1m的立方體，采用C3D8T繪制1個(gè)單元。 GSZLGHGS5({XV$PHN$JX7F1.png 設(shè)置兩個(gè)分析步：初始狀態(tài)設(shè)置域溫度為293K；分析步1約束X負(fù)方向上四個(gè)節(jié)點(diǎn)X方向位移，X正方向上四個(gè)節(jié)點(diǎn)加載X方向位移0.1m；分析步2設(shè)置整個(gè)域溫度為393K。

$4B$@VY{6R_6Q537{$8`I2]T.png$

經(jīng)過(guò)對(duì)比，Model-1和Model-2的計(jì)算結(jié)果完全一致。下面僅展示Model-1的計(jì)算結(jié)果：前兩個(gè)圖是分析步1的X方向應(yīng)力和總應(yīng)變，后圖為分析步2的X方向應(yīng)力和彈性應(yīng)變（采用位移加載，分析步2的總應(yīng)變不變）。

G3YSO]GUWN@A8$%9`8{X(L8.png Y%ZT8NNGX~~A6~J0APXK%G3.png X8IO]U5V}7`)N)}VEMA54FA.png M0NA6H)K(21V{W(LT@1GAU4.png

下面對(duì)計(jì)算結(jié)果進(jìn)行討論。筆者經(jīng)過(guò)有限元計(jì)算，如果不考慮熱膨脹，分析步2計(jì)算完成后應(yīng)力為3e10Pa。由于自由熱膨脹不產(chǎn)生應(yīng)力，因此線膨脹系數(shù)的引入造成了應(yīng)力的減小。

特別注意的是，在有限元建模過(guò)程中，自定義的材料并未使用Abaqus自帶的熱膨脹項(xiàng)，完全是通過(guò)編寫(xiě)的UMAT來(lái)計(jì)算材料的熱膨脹貢獻(xiàn)。那么，如果我們額外增加Abaqus自帶的熱膨脹項(xiàng)（材料屬性里的expansion項(xiàng)），那么彈性應(yīng)變會(huì)變?yōu)?.098，應(yīng)力為變?yōu)?.94e10MPa。也就是說(shuō)我們計(jì)算了兩次熱膨脹系數(shù)。筆者做出了如下猜想：一旦設(shè)置了熱膨脹屬性，進(jìn)入Abaqus的應(yīng)變就從總應(yīng)變張量變?yōu)榱藦椥詰?yīng)變張量。因此，如果需要編寫(xiě)含溫度參數(shù)的UMAT子程序，我們要么使用Abaqus自帶的expansion項(xiàng)，要么在UMAT中將熱膨脹貢獻(xiàn)考慮，切勿重復(fù)。

NN8G@OQ624(R9I9I%~]@AVQ.png 上述討論僅針對(duì)了應(yīng)變和應(yīng)變?cè)隽浚S多本構(gòu)模型采用的是變形梯度描述的，后續(xù)有時(shí)間筆者會(huì)對(duì)此再做一次驗(yàn)證。