搜索什么是UMAT（User-defined material）用戶定義材料？
2016-11-16 by:CAE仿真在線 來(lái)源:互聯(lián)網(wǎng)
1、什么時(shí)候用用戶定義材料(User-defined material, UMAT)?
很簡(jiǎn)單,當(dāng)ABAQUS沒(méi)有提供我們需要的材料模型時(shí)。所以,在決定自己定義一種新的材料模型之前,最好對(duì)ABAQUS已經(jīng)提供的模型心中有數(shù),并且盡量使用現(xiàn)有的模型,因?yàn)檫@些模型已經(jīng)經(jīng)過(guò)詳細(xì)的驗(yàn)證,并被廣泛接受。
2、好學(xué)嗎?需要哪些基礎(chǔ)知識(shí)?
先看一下ABAQUS手冊(cè)(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ǔ)知識(shí)不很豐富者就只能望洋興嘆,因?yàn)槲覀兊娜蝿?wù)不是開(kāi)發(fā)一套完整的有限元軟件,而只是提供一個(gè)描述材料力學(xué)性能的本構(gòu)方程(Constitutive equation)而已。當(dāng)然,最基本的一些概念和知識(shí)還是要具備的,比如
應(yīng)力(stress),應(yīng)變(strain)及其分量; volumetric part和deviatoric part;模量(modulus)、泊松比(Poisson’s ratio)、拉美常數(shù)(Lame constant);矩陣的加減乘除甚至求逆;還有一些高等數(shù)學(xué)知識(shí)如積分、微分等。
3、UMAT的基本任務(wù)?
我們知道,有限元計(jì)算(增量方法)的基本問(wèn)題是:
已知第n步的結(jié)果(應(yīng)力,應(yīng)變等) ,; 然后給出一個(gè)應(yīng)變?cè)隽?/span>, 計(jì)算新的應(yīng)力 。 UMAT要完成這一計(jì)算,并要計(jì)算Jacobian矩陣DDSDDE(I,J) =。是應(yīng)力增量矩陣(張量或許更合適), 是應(yīng)變?cè)隽烤仃嚒?/span>DDSDDE(I,J) 定義了第J個(gè)應(yīng)變分量的微小變化對(duì)第I 個(gè)應(yīng)力分量帶來(lái)的變化。該矩陣只影響收斂速度,不影響計(jì)算結(jié)果的準(zhǔn)確性(當(dāng)然,不收斂自然得不到結(jié)果)。
4、怎樣建立自己的材料模型?
本構(gòu)方程就是描述材料應(yīng)力應(yīng)變(增量)關(guān)系的數(shù)學(xué)公式,不是憑空想象出來(lái)的,而是根據(jù)實(shí)驗(yàn)結(jié)果作出的合理歸納。比如對(duì)彈性材料,實(shí)驗(yàn)發(fā)現(xiàn)應(yīng)力和應(yīng)變同步線性增長(zhǎng),所以用一個(gè)簡(jiǎn)單的數(shù)學(xué)公式描述。為了解釋彈塑性材料的實(shí)驗(yàn)現(xiàn)象,又提出了一些彈塑性模型,并用數(shù)學(xué)公式表示出來(lái)。
對(duì)各向同性材料(Isotropic material),經(jīng)常采用的辦法是先研究材料單向應(yīng)力-應(yīng)變規(guī)律(如單向拉伸、壓縮試驗(yàn)),并用一數(shù)學(xué)公式加以描述,然后把講該規(guī)律推廣到各應(yīng)力分量。這叫做“泛化“(generalization)。
5、一個(gè)完整的例子及解釋
下面這個(gè)UMAT取自ABAQUS手冊(cè),是一個(gè)用于大變形下的彈塑性材料模型。希望我的注釋能幫助初學(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)力矩陣,在增量步的開(kāi)始,保存并作為已知量傳入UMAT ;在增量步的結(jié)束應(yīng)該保存更新的應(yīng)力;
STRAN--當(dāng)前應(yīng)變,已知 。
DSTRAN—應(yīng)變?cè)隽?/span>,已知。
STATEV--狀態(tài)變量矩陣,用來(lái)保存用戶自己定義的一些變量,如累計(jì)塑性應(yīng)變,粘彈性應(yīng)變等等。增量步開(kāi)始時(shí)作為已知量傳入,增量步結(jié)束應(yīng)該更新;
DDSDDE=。需要更新
DTIME—時(shí)間增量dt。已知。
NDI—正應(yīng)力、應(yīng)變個(gè)數(shù),對(duì)三維問(wèn)題、軸對(duì)稱問(wèn)題自然是3(11,22,33),平面問(wèn)題是2(11,22);已知。
NSHR —剪應(yīng)力、應(yīng)變個(gè)數(shù),三維問(wèn)題時(shí)3(12,13,23),軸對(duì)稱問(wèn)題是1(12);已知。
NTENS=NTENS+ NSHR,已知。
PROPS材料常數(shù)矩陣,如模量啊,粘度系數(shù)啊等等;作為已知量傳入,已知。
DROT—對(duì)finite strain問(wèn)題,應(yīng)變應(yīng)該排除旋轉(zhuǎn)部分,該矩陣提供了旋轉(zhuǎn)矩陣,詳見(jiàn)下面的解釋。已知。
PNEWDT—可用來(lái)控制時(shí)間步的變化。如果設(shè)置為小于1的數(shù),則程序放棄當(dāng)前計(jì)算,并用新的時(shí)間增量DTIME X PNEWDT作為新的時(shí)間增量計(jì)算;這對(duì)時(shí)間相關(guān)的材料如聚合物等有用;如果設(shè)為大余1的數(shù),則下一個(gè)增量步加大DTIME為DTIME X PNEWDT。可以更新。
其他變量含義可參看手冊(cè),暫時(shí)用不到。
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
局部變量,用來(lái)暫時(shí)保存彈性應(yīng)變、塑性應(yīng)變分量以及流動(dò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矩陣很容易計(jì)算
注意,在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),分別保存在EELAS和EPLAS中;
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è)沒(méi)有發(fā)生塑性流動(dòng),按完全彈性變形計(jì)算試算應(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計(jì)算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得到對(duì)應(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ā)生,計(jì)算流動(dòng)方向
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è)隽?/span>
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 計(jì)算塑性變形下的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)變量中,并傳到下一個(gè)增量步
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)變,利用插值的方法得到對(duì)應(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
相關(guān)標(biāo)簽搜索：什么是UMAT（User-defined material）用戶定義材料？ abaqus分析培訓(xùn) abaqus技術(shù)教程 abaqus巖土分析 鋼筋混凝土仿真 abaqus分析理論 abaqus軟件下載 abaqus umat用戶子程序編程 Abaqus代做 Abaqus基礎(chǔ)知識(shí) Fluent、CFX流體分析 HFSS電磁分析 Ansys培訓(xùn)
