什么是UMAT(User-defined material)用戶定義材料?

2016-11-16  by:CAE仿真在線  來源:互聯(lián)網(wǎng)

1、什么時候用用戶定義材料(User-defined material, UMAT)?


很簡單,當(dāng)ABAQUS沒有提供我們需要的材料模型時。所以,在決定自己定義一種新的材料模型之前,最好對ABAQUS已經(jīng)提供的模型心中有數(shù),并且盡量使用現(xiàn)有的模型,因為這些模型已經(jīng)經(jīng)過詳細的驗證,并被廣泛接受。


2、好學(xué)嗎?需要哪些基礎(chǔ)知識?


先看一下ABAQUS手冊(ABAQUS Analysis User's Manual)里的一段話:

Warning: The use of this option generally requires considerable expertise. The user is cautioned that the implementation of any realistic constitutive model requires extensive development and testing. Initial testing on a single element model with prescribed traction loading is strongly recommended.


但這并不意味著非力學(xué)專業(yè),或者力學(xué)基礎(chǔ)知識不很豐富者就只能望洋興嘆,因為我們的任務(wù)不是開發(fā)一套完整的有限元軟件,而只是提供一個描述材料力學(xué)性能的本構(gòu)方程(Constitutive equation)而已。當(dāng)然,最基本的一些概念和知識還是要具備的,比如

應(yīng)力(stress),應(yīng)變(strain)及其分量; volumetric part和deviatoric part;模量(modulus)、泊松比(Poissons ratio)、拉美常數(shù)(Lame constant);矩陣的加減乘除甚至求逆;還有一些高等數(shù)學(xué)知識如積分、微分等。


3、UMAT的基本任務(wù)?

我們知道,有限元計算(增量方法)的基本問題是:

已知第n步的結(jié)果(應(yīng)力,應(yīng)變等) ,; 然后給出一個應(yīng)變增量, 計算新的應(yīng)力 。 UMAT要完成這一計算,并要計算Jacobian矩陣DDSDDE(I,J) =。是應(yīng)力增量矩陣(張量或許更合適), 是應(yīng)變增量矩陣。DDSDDE(I,J) 定義了第J個應(yīng)變分量的微小變化對第I 個應(yīng)力分量帶來的變化。該矩陣只影響收斂速度,不影響計算結(jié)果的準確性(當(dāng)然,不收斂自然得不到結(jié)果)。


4、怎樣建立自己的材料模型?

本構(gòu)方程就是描述材料應(yīng)力應(yīng)變(增量)關(guān)系的數(shù)學(xué)公式,不是憑空想象出來的,而是根據(jù)實驗結(jié)果作出的合理歸納。比如對彈性材料,實驗發(fā)現(xiàn)應(yīng)力和應(yīng)變同步線性增長,所以用一個簡單的數(shù)學(xué)公式描述。為了解釋彈塑性材料的實驗現(xiàn)象,又提出了一些彈塑性模型,并用數(shù)學(xué)公式表示出來。

對各向同性材料(Isotropic material),經(jīng)常采用的辦法是先研究材料單向應(yīng)力-應(yīng)變規(guī)律(如單向拉伸、壓縮試驗),并用一數(shù)學(xué)公式加以描述,然后把講該規(guī)律推廣到各應(yīng)力分量。這叫做“泛化“(generalization)。


5、一個完整的例子及解釋

下面這個UMAT取自ABAQUS手冊,是一個用于大變形下的彈塑性材料模型。希望我的注釋能幫助初學(xué)者理解。需要了解J2理論。


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

1 DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,

2 CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,

3 PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)


STRESS--應(yīng)力矩陣,在增量步的開始,保存并作為已知量傳入UMAT ;在增量步的結(jié)束應(yīng)該保存更新的應(yīng)力;

STRAN--當(dāng)前應(yīng)變,已知

DSTRAN應(yīng)變增量,已知。

STATEV--狀態(tài)變量矩陣,用來保存用戶自己定義的一些變量,如累計塑性應(yīng)變,粘彈性應(yīng)變等等。增量步開始時作為已知量傳入,增量步結(jié)束應(yīng)該更新;

DDSDDE=。需要更新

DTIME時間增量dt。已知。

NDI正應(yīng)力、應(yīng)變個數(shù),對三維問題、軸對稱問題自然是3(11,22,33),平面問題是2(11,22);已知。

NSHR剪應(yīng)力、應(yīng)變個數(shù),三維問題時3(12,13,23),軸對稱問題是1(12);已知。

NTENS=NTENS+ NSHR,已知。

PROPS材料常數(shù)矩陣,如模量啊,粘度系數(shù)啊等等;作為已知量傳入,已知。

DROT—對finite strain問題,應(yīng)變應(yīng)該排除旋轉(zhuǎn)部分,該矩陣提供了旋轉(zhuǎn)矩陣,詳見下面的解釋。已知。

PNEWDT可用來控制時間步的變化。如果設(shè)置為小于1的數(shù),則程序放棄當(dāng)前計算,并用新的時間增量DTIME X PNEWDT作為新的時間增量計算;這對時間相關(guān)的材料如聚合物等有用;如果設(shè)為大余1的數(shù),則下一個增量步加大DTIME為DTIME X PNEWDT??梢愿?。

其他變量含義可參看手冊,暫時用不到。

C

INCLUDE 'ABA_PARAM.INC'

定義了一些參數(shù),變量什么的,不用管

C

CHARACTER*8 CMNAME

C

DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),

1 DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),

2 PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3),

3 DFGRD0(3,3),DFGRD1(3,3)

矩陣的尺寸聲明

C

C LOCAL ARRAYS

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

C EELAS - ELASTIC STRAINS

C EPLAS - PLASTIC STRAINS

C FLOW - DIRECTION OF PLASTIC FLOW

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

C

局部變量,用來暫時保存彈性應(yīng)變、塑性應(yīng)變分量以及流動方向

DIMENSION EELAS(6),EPLAS(6),FLOW(6)

C

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

1 ENUMAX=.4999D0,NEWTON=10,TOLER=1.0D-6)

C

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

C UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC MISES PLASTICITY

C CANNOT BE USED FOR PLANE STRESS

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

C PROPS(1) - E

C PROPS(2) - NU

C PROPS(3..) - SYIELD AN HARDENING DATA

C CALLS HARDSUB FOR CURVE OF YIELD STRESS VS. PLASTIC STRAIN

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

C

C ELASTIC PROPERTIES

C

獲取楊氏模量,泊松比,作為已知量由PROPS向量傳入

EMOD=PROPS(1) E

ENU=PROPS(2) ν

EBULK3=EMOD/(ONE-TWO*ENU) 3K

EG2=EMOD/(ONE+ENU) 2G

EG=EG2/TWO G

EG3=THREE*EG 3G

ELAM=(EBULK3-EG2)/THREE λ

DO K1=1,NTENS

DO K2=1,NTENS

DDSDDE(K1,K2)=ZERO

END DO

END DO

彈性部分,Jacobian矩陣很容易計算



注意,在ABAQUS中,剪切應(yīng)變采用工程剪切應(yīng)變的定義,所以剪切部分模量是G而不是2G!


C

C ELASTIC STIFFNESS

C

DO K1=1,NDI

DO K2=1,NDI

DDSDDE(K2,K1)=ELAM

END DO

DDSDDE(K1,K1)=EG2+ELAM

END DO

DO K1=NDI+1,NTENS

DDSDDE(K1,K1)=EG

END DO

C

C RECOVER ELASTIC AND PLASTIC STRAINS AND ROTATE FORWARD

C ALSO RECOVER EQUIVALENT PLASTIC STRAIN

C

讀取彈性應(yīng)變分量,塑性應(yīng)變分量,并旋轉(zhuǎn)(調(diào)用了ROTSIG),分別保存在EELASEPLAS中;

CALL ROTSIG(STATEV( 1),DROT,EELAS,2,NDI,NSHR)

CALL ROTSIG(STATEV(NTENS+1),DROT,EPLAS,2,NDI,NSHR)

讀取等效塑性應(yīng)變

EQPLAS=STATEV(1+2*NTENS)

先假設(shè)沒有發(fā)生塑性流動,按完全彈性變形計算試算應(yīng)力


C

C CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN

C

DO K1=1,NTENS

DO K2=1,NTENS

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

END DO

EELAS(K1)=EELAS(K1)+DSTRAN(K1)

END DO

C計算Mises應(yīng)力

C CALCULATE EQUIVALENT VON MISES STRESS

C

SMISES=(STRESS(1)-STRESS(2))**2+(STRESS(2)-STRESS(3))**2

1 +(STRESS(3)-STRESS(1))**2

DO K1=NDI+1,NTENS

SMISES=SMISES+SIX*STRESS(K1)**2

END DO

SMISES=SQRT(SMISES/TWO)

C 根據(jù)當(dāng)前等效塑性應(yīng)變,調(diào)用HARDSUB得到對應(yīng)的屈服應(yīng)力

C GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE

C

NVALUE=NPROPS/2-1

CALL HARDSUB(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)

C

C DETERMINE IF ACTIVELY YIELDING

C 如果Mises應(yīng)力大余屈服應(yīng)力,屈服發(fā)生,計算流動方向

IF (SMISES.GT.(ONE+TOLER)*SYIEL0) THEN

C

C ACTIVELY YIELDING

C SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESS

C CALCULATE THE FLOW DIRECTION

C

SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE

DO K1=1,NDI

FLOW(K1)=(STRESS(K1)-SHYDRO)/SMISES

END DO

DO K1=NDI+1,NTENS

FLOW(K1)=STRESS(K1)/SMISES

END DO

C根據(jù)J2理論并應(yīng)用Newton-Rampson方法求得等效塑性應(yīng)變增量

C SOLVE FOR EQUIVALENT VON MISES STRESS

C AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATION

C

SYIELD=SYIEL0

DEQPL=ZERO

DO KEWTON=1,NEWTON

RHS=SMISES-EG3*DEQPL-SYIELD

DEQPL=DEQPL+RHS/(EG3+HARD)

CALL HARDSUB(SYIELD,HARD,EQPLAS+DEQPL,PROPS(3),NVALUE)

IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 10

END DO

C

C WRITE WARNING MESSAGE TO THE .MSG FILE

C

WRITE(7,2) NEWTON

2 FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',

1 'CONVERGE AFTER ',I3,' ITERATIONS')

10 CONTINUE

C更新應(yīng)力,應(yīng)變分量

C UPDATE STRESS, ELASTIC AND PLASTIC STRAINS AND

C EQUIVALENT PLASTIC STRAIN

C

DO K1=1,NDI

STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO

EPLAS(K1)=EPLAS(K1)+THREE/TWO*FLOW(K1)*DEQPL

EELAS(K1)=EELAS(K1)-THREE/TWO*FLOW(K1)*DEQPL

END DO

DO K1=NDI+1,NTENS

STRESS(K1)=FLOW(K1)*SYIELD

EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL

EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL

END DO

EQPLAS=EQPLAS+DEQPL

C

C CALCULATE PLASTIC DISSIPATION

C

SPD=DEQPL*(SYIEL0+SYIELD)/TWO

C

C 計算塑性變形下的Jacobian矩陣

FORMULATE THE JACOBIAN (MATERIAL TANGENT)

C FIRST CALCULATE EFFECTIVE MODULI

C

EFFG=EG*SYIELD/SMISES

EFFG2=TWO*EFFG

EFFG3=THREE/TWO*EFFG2

EFFLAM=(EBULK3-EFFG2)/THREE

EFFHRD=EG3*HARD/(EG3+HARD)-EFFG3

c...

if (props(7).lt..001) go to 99

c...

DO K1=1,NDI

DO K2=1,NDI

DDSDDE(K2,K1)=EFFLAM

END DO

DDSDDE(K1,K1)=EFFG2+EFFLAM

END DO

DO K1=NDI+1,NTENS

DDSDDE(K1,K1)=EFFG

END DO

DO K1=1,NTENS

DO K2=1,NTENS

DDSDDE(K2,K1)=DDSDDE(K2,K1)+EFFHRD*FLOW(K2)*FLOW(K1)

END DO

END DO

c...

99 continue

c...

ENDIF

C將彈性應(yīng)變,塑性應(yīng)變分量保存到狀態(tài)變量中,并傳到下一個增量步

C STORE ELASTIC AND (EQUIVALENT) PLASTIC STRAINS

C IN STATE VARIABLE ARRAY

C

DO K1=1,NTENS

STATEV(K1)=EELAS(K1)

STATEV(K1+NTENS)=EPLAS(K1)

END DO

STATEV(1+2*NTENS)=EQPLAS

C

RETURN

END

c...

c...子程序,根據(jù)等效塑性應(yīng)變,利用插值的方法得到對應(yīng)的屈服應(yīng)力

SUBROUTINE HARDSUB(SYIELD,HARD,EQPLAS,TABLE,NVALUE)

C

INCLUDE 'ABA_PARAM.INC'

C

DIMENSION TABLE(2,NVALUE)

C

PARAMETER(ZERO=0.D0)

C

C SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO

C

SYIELD=TABLE(1,NVALUE)

HARD=ZERO

C IF MORE THAN ONE ENTRY, SEARCH TABLE

C

IF(NVALUE.GT.1) THEN

DO K1=1,NVALUE-1

EQPL1=TABLE(2,K1+1)

IF(EQPLAS.LT.EQPL1) THEN

EQPL0=TABLE(2,K1)

IF(EQPL1.LE.EQPL0) THEN

WRITE(7,1)

1 FORMAT(//,30X,'***ERROR - PLASTIC STRAIN MUST BE `,

1 `ENTERED IN ASCENDING ORDER')

CALL XIT

ENDIF

C

C CURRENT YIELD STRESS AND HARDENING

C

DEQPL=EQPL1-EQPL0

SYIEL0=TABLE(1,K1)

SYIEL1=TABLE(1,K1+1)

DSYIEL=SYIEL1-SYIEL0

HARD=DSYIEL/DEQPL

SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD

GOTO 10

ENDIF

END DO

10 CONTINUE

ENDIF

RETURN

END


開放分享:優(yōu)質(zhì)有限元技術(shù)文章,助你自學(xué)成才

相關(guān)標簽搜索:什么是UMAT(User-defined material)用戶定義材料? abaqus分析培訓(xùn) abaqus技術(shù)教程 abaqus巖土分析 鋼筋混凝土仿真 abaqus分析理論 abaqus軟件下載 abaqus umat用戶子程序編程 Abaqus代做 Abaqus基礎(chǔ)知識 Fluent、CFX流體分析 HFSS電磁分析 Ansys培訓(xùn) 

編輯
在線報名:
  • 客服在線請直接聯(lián)系我們的客服,您也可以通過下面的方式進行在線報名,我們會及時給您回復(fù)電話,謝謝!
驗證碼

全國服務(wù)熱線

1358-032-9919

廣州公司:
廣州市環(huán)市中路306號金鷹大廈3800
電話:13580329919
          135-8032-9919
培訓(xùn)QQ咨詢:點擊咨詢 點擊咨詢
項目QQ咨詢:點擊咨詢
email:kf@1cae.com