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 CMNAME
      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),DFGRD0(3,3),DFGRD1(3,3)
      ! user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD 
      ! and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT
      RETURN
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) C INCLUDE 'ABA_PARAM.INC' ! 定義了一些相關參數與變量什么，從 ABAQUS 安裝目錄下的子文件夾“… \site”中可找到 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 , 3k = E /(1-2ν ) EG2=EMOD/(ONE+ENU) ! 2G , 2G = E /(1+υ) EG=EG2/TWO ! G , G = E /2(1+υ) EG3=THREE*EG ! 3G ELAM=(EBULK3-EG2)/THREE ! λ , λ = (3k-2G)/3 DO K1=1,NTENS DO K2=1,NTENS DDSDDE(K1,K2)=ZERO END DO END 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,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 ! 讀取彈性應變分量，塑性應變分量，並旋轉（調用了 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 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 根據當前等效塑性應變，調用 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 應力大余屈服應力，屈服發生，計算流動方向 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 根據 J2 理論並應用 Newton-Rampson 方法求得等效塑性應變增量 C SOLVE FOR EQUIVALENT VON MISES STRESS C AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATIO C 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 更新應力σ n+1，應變分量 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... END IF C 將彈性應變，塑性應變分量保存到狀態變量中，並傳到下一個增量步 C STORE ELASTIC AND (EQUIVALENT) PLASTIC STRAINS C IN STATE VARIABLE ARRAYC DO K1=1,NTENS STATEV(K1)=EELAS(K1) STATEV(K1+NTENS)=EPLAS(K1) END DO STATEV(1+2*NTENS)=EQPLAS C RETURN 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 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 END IF CC 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 END IF END DO 10 CONTINUE END IF RETURN END SUBROUTINE HARDSUB