mises本構模型UMAT（附源代碼和詳細注釋）
一. ABAQUS自帶mises本構與UMAT對比
二. 如何在分析中使用自定義材料
1.通過ABAQUS/CAE使用自定義材料
在Property模塊中，執行【Material】/【Create】命令（或單擊相關工具箱區的按鈕），彈出Edit Material對話框（如圖所示），用戶可以通過該對話框選擇材料模型、設置材料參數。對于自定義模型，執行對話框的【General】/【User Material】命令，此時在Material Behaviors區域中會出現User Material字樣，表明定義的是用戶材料。在Edit Material對話框下方的User Material 區域中的【User material type】下拉列表中有三個選項，分別為Mechanical（力學材料）、 Thermal (熱學材料）和Thermomechanical (熱力學模型），默認選項為Mechanical。對于一些巖土材料模型，尤其是采用非關聯流動法則的模型，雅克比矩陣是不對稱的，此時需勾選Use unsymmetric material stiffiiess matrix 復選框。材料的參數在 Data 區域中的【Mechanical Constants】列表中輸入，這里的數據會按次序傳給UMAT子程序中的PROPS (NPROPS)數組，數據的個數即為NPROPS。
注：輸入材料參數時，在數據列表尾部按回車健會自動增加一行數據.此外，在數據窗口上右擊會顯示彈出式菜單，可實現數據的拷貝、增加、刪除等功能。
如果UMAT子程序用到狀態變量,還需設置狀態變量的個數。具體操作仍然在Edit Material對話框中進行，執行對話框的【General】/【Depvar】命令，在【Number of solution-dependent state variables】輸入框中設罝狀態變量的個數。
2.通過inp輸入文件使用自定義材料
在inp輸入文件中的材料定義選項塊（* Material)中使用下列語句即可：
*User material,type=Mechanical,constants= number of constants,Unsymm
關鍵字User material表明定義的用戶材料;type為材料類型，可為Mechanical(力學材料）、 Thermal (熱學材料）和Thermomechanical (熱力學模型），Mechanical是默認選項：Constants 選項給出的是材料參數的個數；Unsymm選項只有當雅克比矩陣非對稱時才采用。該關鍵字行下的數據行中的數據為材料參數。
狀態變量的個數同樣在材料定義選項塊中定義，相應的關鍵字行為：
*Depvar
Number of solution-dependent state variables：狀態變量個數
三. Mises彈塑性模型
4.源程序和注釋
需要啟用沙漏控制：
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,
3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,
4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 MATERL
DIMENSION STRESS(NTENS),STATEV(NSTATV),
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),
4 DFGRD0(3,3),DFGRD1(3,3)
C
DIMENSION EELAS(6),EPLAS(6),FLOW(6)
C EELAS(1-6)為彈性應變分量,EPLAS(1-6)為塑性應變分量,FLOW(6)為屈服面對應力的導數
PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
DATA NEWTON,TOLER/10,1.D-6/
C
C -----------------------------------------------------------
C UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC PLASTICITY
C J2 FLOW THEORY
C CAN NOT BE USED FOR PLANE STRESS
C -----------------------------------------------------------
C PROPS(1) - E ,為彈性模量
C PROPS(2) - NU ,為泊松比
C PROPS(3) - SYIELD, 為屈服應力
C CALLS AHARD FOR CURVE OF SYIELD VS. PEEQ
C -----------------------------------------------------------
C 程序不能用于平面應變
IF (NDI.NE.3) THEN
WRITE(6,1)
1 FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',
1 'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')
ENDIF
C
C ELASTIC PROPERTIES
C
EMOD=PROPS(1)
ENU=PROPS(2)
IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
C 得到拉梅常數λ(ELAM),剪切模量G（EG）
C ELASTIC STIFFNESS
C
DO 20 K1=1,NTENS
DO 10 K2=1,NTENS
DDSDDE(K2,K1)=0.0
10 CONTINUE
20 CONTINUE
C
DO 40 K1=1,NDI
DO 30 K2=1,NDI
DDSDDE(K2,K1)=ELAM
30 CONTINUE
DDSDDE(K1,K1)=EG2+ELAM
40 CONTINUE
DO 50 K1=NDI+1,NTENS
DDSDDE(K1,K1)=EG
50 CONTINUE
C 得到（6-1）中的彈性模量矩陣
C CALCULATE STRESS FROM ELASTIC STRAINS
C 假設應變增量全為彈性，計算應力增量（即第一次應力預測）
DO 70 K1=1,NTENS
DO 60 K2=1,NTENS
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
60 CONTINUE
70 CONTINUE
C
C RECOVER ELASTIC AND PLASTIC STRAINS
C
DO 80 K1=1,NTENS
EELAS(K1)=STATEV(K1)+DSTRAN(K1)
C STATEV(1-6)存儲彈性應變分量，此時假定應變增量全為彈性
EPLAS(K1)=STATEV(K1+NTENS)
C STATEV(6-12)存儲塑性應變分量
80 CONTINUE
EQPLAS=STATEV(1+2*NTENS)
C STATEV(13)為等效塑性應變
C 如果沒有給出屈服應力，材料認為是彈性的
C
IF(NPROPS.GT.2.AND.PROPS(3).GT.0.0) THEN
C
C MISES STRESS
C
SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2)) +
1 (STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3)) +
1 (STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1))
DO 90 K1=NDI+1,NTENS
SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1)
90 CONTINUE
SMISES=SQRT(SMISES/TWO)
C 計算mises應力
C HARDENING CURVE, GET YIELD STRESS
C 等效塑性應變和硬化曲線確定硬化后的屈服應力,初始等效塑性應變為0,
C 對應的屈服應力為初始屈服應力
NVALUE=NPROPS/2-1
CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)
C
C DETERMINE IF ACTIVELY YIELDING
C
IF (SMISES.GT.(1.0+TOLER)*SYIEL0) THEN
C 如果Mises應力超出了屈服應力,判斷發生了屈服,進入了塑性階段
C FLOW DIRECTION
C
SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
C SHYDRO為平均應力
ONESY=ONE/SMISES
DO 110 K1=1,NDI
FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO)
110 CONTINUE
DO 120 K1=NDI+1,NTENS
FLOW(K1)=STRESS(K1)*ONESY
120 CONTINUE
C FLOW(1-6)為 sij/q
C SOLVE FOR EQUIV STRESS, NEWTON ITERATION
C
SYIELD=SYIEL0
DEQPL=0.0
DO 130 KEWTON=1,NEWTON
RHS=SMISES-EG3*DEQPL-SYIELD
C 計算屈服面函數誤差
DEQPL=DEQPL+RHS/(EG3+HARD)
C 按(6-13)計算新的等效塑性應變
CALL AHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(3),NVALUE)
C 由新的等效塑性應變獲得新的屈服應力
IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 140
C 如果滿足收斂要求結束,否則再次迭代
130 CONTINUE
WRITE(6,2) NEWTON
2 FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',
1 'CONVERGE AFTER ',I3,' ITERATIONS')
140 CONTINUE
EFFHRD=EG3*HARD/(EG3+HARD)
C
C
C CALC STRESS AND UPDATE STRAINS
C 計算及更新應力,彈性應變,塑性應變
DO 150 K1=1,NDI
C 更新正應力分量
STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO
C 按(6-16)計算應力
EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO
C 更新塑性應變
EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO
C 更新彈性應變
150 CONTINUE
DO 160 K1=NDI+1,NTENS
C 更新剪應力分量
STRESS(K1)=FLOW(K1)*SYIELD
EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL
EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL
160 CONTINUE
EQPLAS=EQPLAS+DEQPL
SPD=DEQPL*(SYIEL0+SYIELD)/TWO
C
C JACOBIAN
C 計算雅克比矩陣
C
EFFG=EG*SYIELD/SMISES
EFFG2=TWO*EFFG
EFFG3=THREE*EFFG2/TWO
EFFLAM=(EBULK3-EFFG2)/THREE
DO 220 K1=1,NDI
DO 210 K2=1,NDI
DDSDDE(K2,K1)=EFFLAM
210 CONTINUE
DDSDDE(K1,K1)=EFFG2+EFFLAM
220 CONTINUE
DO 230 K1=NDI+1,NTENS
DDSDDE(K1,K1)=EFFG
230 CONTINUE
DO 250 K1=1,NTENS
DO 240 K2=1,NTENS
DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1)
1 *(EFFHRD-EFFG3)
240 CONTINUE
250 CONTINUE
ENDIF
ENDIF
C
C STORE STRAINS IN STATE VARIABLE ARRAY
C
DO 310 K1=1,NTENS
STATEV(K1)=EELAS(K1)
STATEV(K1+NTENS)=EPLAS(K1)
310 CONTINUE
STATEV(1+2*NTENS)=EQPLAS
C
RETURN
END
C
C
SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE)
C 根據導入的硬化曲線及當前的等效塑性應變,求得屈服應力及屈服應力對等效塑性應變的偏導
INCLUDE 'ABA_PARAM.INC'
DIMENSION TABLE(2,NVALUE)
C
C SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
SYIELD=TABLE(1,NVALUE)
HARD=0.0
C
C IF MORE THAN ONE ENTRY, SEARCH TABLE
C
IF(NVALUE.GT.1) THEN
DO 10 K1=1,NVALUE-1
EQPL1=TABLE(2,K1+1)
IF(EQPLAS.LT.EQPL1) THEN
EQPL0=TABLE(2,K1)
IF(EQPL1.LE.EQPL0) THEN
WRITE(6,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 20
ENDIF
10 CONTINUE
20 CONTINUE
ENDIF
RETURN
END
源程序附件：
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP
