Abaqus用戶子程序umat的學習

說明：在文件中，！后面的內容為注釋內容。本文為學習心得，很多注釋是自己摸索得到。如有不正確的地方，敬請指正。
! ——————————————————————————
! 1、為何需要使用用戶材料子程序（ User-Defined Material, UMAT）？
! 很簡單，當 ABAQUS 沒有提供我們需要的材料模型時。所以，在決定自己定義一種新的材料模型之前，最好
! 對 ABAQUS 已經提供的模型心中有數，并且盡量使用現有的模型，因為這些模型已經經過詳細的驗證，并被廣泛接受。
! ——————————————————————————
! UMAT 子程序具有強大的功能，使用 UMAT 子程序：
! (1)可以定義材料的本構關系，使用 ABAQUS 材料庫中沒有包含的材料進行計算，擴充程序功能。
! (2) 幾乎可以用于力學行為分析的任何分析過程，幾乎可以把用戶材料屬性賦予 ABAQUS 中的任何單元。
! (3) 必須在 UMAT 中提供材料本構模型的雅可比（ Jacobian）矩陣，即應力增量對應變增量的變化率。
! (4) 可以和用戶子程序“ USDFLD”聯合使用，通過“ USDFLD”重新定義單元每一物質點上傳遞到 UMAT 中場變量的數值。
! ——————————————————————————
! 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.
! 但這并不意味著非力學專業，或者力學基礎知識不很豐富者就只能望洋興嘆，因為
! 我們的任務不是開發一套完整的有限元軟件，而只是提供一個描述材料力學性能的本構
! 方程（ Constitutive equation）而已。當然，最基本的一些概念和知識還是要具備的，比如：
! 應力(stress),應變（ strain）及其分量； volumetric part 和 deviatoric part；模量（ modul
! us）、泊松比(Poisson’s ratio)、拉梅常數(Lame constant)；矩陣的加減乘除甚至求逆；還
! 有一些高等數學知識如積分、微分等。
! ——————————————————————————
! 3、 UMAT 的基本任務？
! 我們知道，有限元計算（增量方法）的基本問題是：已知第 n 步的結果（應力，應變等）
! σ[n],ε[n],然后給出一個應變增量 dε[n+1]，計算新的應力σ[n+1]。 UMAT 要完成這一
! 計算， 并要計算 Jacobian 矩陣 DDSDDE(I,J) =Δσ/Δε 。 Δσ 是應力增量矩陣（張量或
! 許更合適）， Δε 是應變增量矩陣。 DDSDDE(I,J) 定義了第 J 個應變分量的微小變化對第 I 個應力分量帶來的變化。
! 該矩陣只影響收斂速度，不影響計算結果的準確性（當然，不收斂自然得不到結果）。
! ——————————————————————————
! 4、怎樣建立自己的材料模型？
! 本構方程就是描述材料應力應變（增量）關系的數學公式，不是憑空想象出來的，
! 而是根據實驗結果作出的合理歸納。比如對彈性材料，實驗發現應力和應變同步線性增
! 長，所以用一個簡單的數學公式描述。為了解釋彈塑性材料的實驗現象，又提出了一些
! 彈塑性模型，并用數學公式表示出來。
! 對各向同性材料（ Isotropic material） ,經常采用的辦法是先研究材料單向應力-應變
! 規律（如單向拉伸、壓縮試驗），并用一數學公式加以描述，然后把該規律推廣到各應
! 力分量。這叫做“泛化“(generalization)。
! ——————————————————————————
! 5、一個完整的例子及解釋
! 由于主程序與 UMAT 之間存在數據傳遞，甚至一些公共變量，因此必須遵循有關
! UMAT 的書寫格式， 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,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)      ! 共37個參數INCLUDE 'ABA_PARAM.INC'CHARACTER*80 CMNAMEDIMENSION 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),DFGRD0(3,3),DFGRD1(3,3)! user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD ! and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDTRETURN
END SUBROUTINE UMAT
! ------------------------------------------------------------------------------
! COORDS 當前積分點的坐標
! DDSDDE ( NTENS，NTENS) 大小為 NTENS×NTENS 的 Jacobian 矩陣( Δσ / Δε )，
! DDSDDE (I,J) 定義了第 J 個應變分量的微小變化對第 I 個應力分量帶來的變化。通常 Jacobian 矩陣是一個對稱矩陣，
! 除非在“ *USER MATERIAL”語句中加入了“ UNSYMM”參數；需要更新 DROT
! 對 Finite strain 問題，應變應該排除旋轉部分，該矩陣提供了旋轉矩陣，詳見下面的解釋；已知
! DSTRAN (NTENS) 應變增量 dε[n+1]，已知
! DTIME 增量步的時間增量 dt ；已知
! KSTEP,KINC 傳到用戶子程序當前的 STEP 和 INCREMENT 值   KSTEP為載荷時間步；KINC為增量步
! NDI 直接應力(正應力)、應變個數，對三維問題、軸對稱問題自然是 3( 11,22,33)，平面問題是 2(11,22)；已知
! NOEL,NPT 積分點所在單元的編號和積分點的編號
! NSHR 剪切應力(剪應力)、應變個數，三維問題時 3(12,13,23)，軸對稱問題是 1(12)；已知
! NTENS =NDI+ NSHR，總應力分量的個數；已知
! PNEWDT 可用來控制時間步的變化。如果設置為小于 1 的數，則程序放棄當前計算，并用新的時間增量 DTIME X PNEWDT
! 作為新的時間增量計算；這對時間相關的材料如聚合物等有用；
! 如果設為大余 1 的數，則下一個增量步加大 DTIME 為 DTIME X PNEWDT。可以更新。
! PROPS (NPROPS) 材料常數數組，如模量啊，粘度系數等等；
! 材料參數的個數，等于關鍵詞“ *USER MATERIAL”中“ CONSTANT S”常數設定的值；
! 矩陣中元素的數值對應于關鍵詞“ USER MATERIAL”下面的數據行。作為已知量傳入；已知
! SSE,SPD,SCD 分別定義每一增量步的彈性應變能，塑性耗散和蠕變耗散。它們對計算結果沒有影響，僅僅作為能量輸出
! STATEV (NSTATEV) 狀態變量矩陣，用來保存用戶自己定義的一些變量，如累計塑性應變，粘彈性應變等等。
! 增量步開始時作為已知量傳入，增量步結束應該更新
! STRAN (NTENS) 當前應變數組ε[n]，已知
! STRESS (NTENS) 應力張量數組，對應 NDI 個直接分量和 NSHR 個剪切分量。
! 在增量步的開始，應力張量矩陣σ[n]中的數值通過 UMAT 和主程序之間的接口傳遞到 UMAT 中，
! 在增量步的結束, UMAT 將對應力張量矩陣更新為σ[n+1]。
! 對于包含剛體轉動的有限應變問題，一個增量步調用 UMAT 之前就
! 已經對應力張量進行了剛體轉動，因此 UMAT 中只需處理應力張量的共旋部分。 UMAT 中應力張量的度量為柯
! 西（真實）應力。
! ------------------------------------------------------------------------------
! 下面這個 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)
CINCLUDE 'ABA_PARAM.INC'
! 定義了一些相關參數與變量什么，從 ABAQUS 安裝目錄下的子文件夾“… \site”中可找到
CCHARACTER*8 CMNAME
CDIMENSION 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)
CPARAMETER(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)              !   EENU=PROPS(2)               !   νEBULK3=EMOD/(ONE-TWO*ENU)  !   3K , 3k = E /(1-2ν )EG2=EMOD/(ONE+ENU)         !   2G , 2G = E /(1+υ)EG=EG2/TWO                 !   G  , G = E /2(1+υ)EG3=THREE*EG               !   3GELAM=(EBULK3-EG2)/THREE    !   λ  ,  λ = (3k-2G)/3DO K1=1,NTENSDO K2=1,NTENSDDSDDE(K1,K2)=ZEROEND DOEND DO                            
! 彈性部分，Jacobian矩陣很容易計算
!   |λ+2G   λ      λ          |
!   |λ      λ+2G   λ          |
! J=|λ      λ      λ+2G       | 
!   |                    G      |
!   |                       G   |
!   |                          G|
! 注意，在ABAQUS中，剪切應變采用工程剪切應變的定義γ(ij)=u(ij)+u(ji)，所以剪切部分模量是G而不是2G!
C  
C        ELASTIC  STIFFNESS  C  DO K1=1,NDIDO K2=1,NDIDDSDDE(K2,K1)=ELAMEND DODDSDDE(K1,K1)=EG2+ELAMEND DODO K1=NDI+1,NTENSDDSDDE(K1,K1)=EGEND DO
C
C RECOVER ELASTIC AND PLASTIC STRAINS AND ROTATE FORWARD
C ALSO RECOVER EQUIVALENT PLASTIC STRAIN
C
! 讀取彈性應變分量，塑性應變分量，并旋轉（調用了 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)
! 先假設沒有發生塑性流動，按完全彈性變形計算試算應力Δσ = J.Δε, σ[n+1]=σ[n]+Δσ
CC CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN
CDO K1=1,NTENSDO K2=1,NTENSSTRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)END DOEELAS(K1)=EELAS(K1)+DSTRAN(K1)  ! 彈性應變分量END DO
C 計算 Mises 應力
C CALCULATE EQUIVALENT VON MISES STRESS
CSMISES=(STRESS(1)-STRESS(2))**2+(STRESS(2)-STRESS(3))**21 +(STRESS(3)-STRESS(1))**2DO K1=NDI+1,NTENSSMISES=SMISES+SIX*STRESS(K1)**2END DOSMISES=SQRT(SMISES/TWO)
C 根據當前等效塑性應變，調用 HARDSUB 得到對應的屈服應力
C GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE
CNVALUE=NPROPS/2-1CALL HARDSUB(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)
C
C DETERMINE IF ACTIVELY YIELDING
C 如果 Mises 應力大余屈服應力，屈服發生，計算流動方向IF (SMISES.GT.(ONE+TOLER)*SYIEL0) THEN
C
C ACTIVELY YIELDING
C SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESS
C CALCULATE THE FLOW DIRECTION
CSHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREEDO K1=1,NDIFLOW(K1)=(STRESS(K1)-SHYDRO)/SMISESEND DODO K1=NDI+1,NTENSFLOW(K1)=STRESS(K1)/SMISESEND DO
C 根據 J2 理論并應用 Newton-Rampson 方法求得等效塑性應變增量
C SOLVE FOR EQUIVALENT VON MISES STRESS
C AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATIO
C
CSYIELD=SYIEL0DEQPL=ZERODO KEWTON=1,NEWTONRHS=SMISES-EG3*DEQPL-SYIELDDEQPL=DEQPL+RHS/(EG3+HARD)CALL HARDSUB(SYIELD,HARD,EQPLAS+DEQPL,PROPS(3),NVALUE)IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 10END DO
C
C WRITE WARNING MESSAGE TO THE .MSG FILE
CWRITE(7,2) NEWTON2 FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',1 'CONVERGE AFTER ',I3,' ITERATIONS')
10    CONTINUE
C 更新應力σ n+1，應變分量
C UPDATE STRESS, ELASTIC AND PLASTIC STRAINS AND
C EQUIVALENT PLASTIC STRAIN
CDO K1=1,NDISTRESS(K1)=FLOW(K1)*SYIELD+SHYDROEPLAS(K1)=EPLAS(K1)+THREE/TWO*FLOW(K1)*DEQPLEELAS(K1)=EELAS(K1)-THREE/TWO*FLOW(K1)*DEQPLEND DODO K1=NDI+1,NTENSSTRESS(K1)=FLOW(K1)*SYIELDEPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPLEELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPLEND DOEQPLAS=EQPLAS+DEQPL
C
C CALCULATE PLASTIC DISSIPATION 
CSPD=DEQPL*(SYIEL0+SYIELD)/TWO
C
C 計算塑性變形下的 Jacobian 矩陣FORMULATE THE JACOBIAN (MATERIAL TANGENT)
C FIRST CALCULATE EFFECTIVE MODULI
CEFFG=EG*SYIELD/SMISESEFFG2=TWO*EFFGEFFG3=THREE/TWO*EFFG2EFFLAM=(EBULK3-EFFG2)/THREEEFFHRD=EG3*HARD/(EG3+HARD)-EFFG3
c...if (props(7).lt..001) go to 99
c...DO K1=1,NDIDO K2=1,NDIDDSDDE(K2,K1)=EFFLAMEND DODDSDDE(K1,K1)=EFFG2+EFFLAMEND DODO K1=NDI+1,NTENSDDSDDE(K1,K1)=EFFGEND DODO K1=1,NTENSDO K2=1,NTENSDDSDDE(K2,K1)=DDSDDE(K2,K1)+EFFHRD*FLOW(K2)*FLOW(K1)END DOEND DO
c...
99    continue
c...END IF
C 將彈性應變，塑性應變分量保存到狀態變量中，并傳到下一個增量步
C STORE ELASTIC AND (EQUIVALENT) PLASTIC STRAINS
C IN STATE VARIABLE ARRAYCDO K1=1,NTENSSTATEV(K1)=EELAS(K1)STATEV(K1+NTENS)=EPLAS(K1)END DOSTATEV(1+2*NTENS)=EQPLAS
CRETURN
END SUBROUTINE UMAT
c...
c...子程序，根據等效塑性應變，利用插值的方法得到對應的屈服應力
SUBROUTINE HARDSUB(SYIELD,HARD,EQPLAS,TABLE,NVALUE)INCLUDE 'ABA_PARAM.INC'DIMENSION TABLE(2,NVALUE)PARAMETER(ZERO=0.D0)
C
C SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
CSYIELD=TABLE(1,NVALUE)HARD=ZERO
C IF MORE THAN ONE ENTRY, SEARCH TABLE
CIF(NVALUE.GT.1) THENDO K1=1,NVALUE-1EQPL1=TABLE(2,K1+1)IF(EQPLAS.LT.EQPL1) THENEQPL0=TABLE(2,K1)IF(EQPL1.LE.EQPL0) THENWRITE(7,1)1 FORMAT(//,30X,'***ERROR - PLASTIC STRAIN MUST BE ',1     'ENTERED IN ASCENDING ORDER')CALL XITEND IF
CC CURRENT YIELD STRESS AND HARDENING
CDEQPL=EQPL1-EQPL0SYIEL0=TABLE(1,K1)SYIEL1=TABLE(1,K1+1)DSYIEL=SYIEL1-SYIEL0HARD=DSYIEL/DEQPLSYIELD=SYIEL0+(EQPLAS-EQPL0)*HARDGOTO 10END IFEND DO
10    CONTINUEEND IFRETURN
END SUBROUTINE HARDSUB

轉載于:https://www.cnblogs.com/baowee/p/9556837.html

總結

以上是生活随笔為你收集整理的Abaqus用户子程序umat的学习的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。