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