目前，對(duì)于材料力學(xué)行為的研究，ABAQUS UMAT技術(shù)幾乎成了標(biāo)配。只要涉及強(qiáng)度預(yù)測(cè)、失效準(zhǔn)則、蠕變、粘彈性、疲勞、應(yīng)變率效應(yīng)、固化變形等等研究，大家的論文中如果沒有本構(gòu)的討論、UMAT或者VUMAT的內(nèi)容，就會(huì)顯得文章沒有深度。即便是用其他的商用軟件，也會(huì)涉及到自定義本構(gòu)的問題。UMAT之于ABAQUS，就像UDF之于Fluent。

國內(nèi)大規(guī)模搞這些子程序研究有20多年了，我個(gè)人接觸這個(gè)東西也有10年。最初接觸它，總覺得它神秘又高端。后面用的多了，又覺得無非是在現(xiàn)有程序基礎(chǔ)上改改參數(shù)或者加點(diǎn)模型。隨著自己開發(fā)有限元軟件，又發(fā)現(xiàn)如何給用戶開放自定義本構(gòu)權(quán)限，實(shí)際上是軟件架構(gòu)設(shè)計(jì)的一個(gè)難點(diǎn)。本文更想討論的是，UMAT為什么是今天看到的這個(gè)樣子，以及未來有哪些可能性。

• 什么是UMAT
• UAMT程序結(jié)構(gòu)
• 如何編寫UAMT
• 如何調(diào)用UAMT
• 如何調(diào)試UAMT
• 如何配置子程序
• 未來設(shè)想

什么是UMAT

UMAT:User-defined mechanical material。顧名思義，就是用戶自己定義材料力學(xué)行為的子程序。一般而言，就是定義材料的應(yīng)力應(yīng)變關(guān)系，即本構(gòu)。最初接觸它的時(shí)候，有一點(diǎn)畏難心理，感覺這個(gè)不好學(xué)，而且麻煩。原因主要有：

（1）準(zhǔn)備工作復(fù)雜。需要安裝適合版本的Fotran、VS，且有安裝順序。然后和ABAQUS關(guān)聯(lián)。稍有不慎，就得重裝。

（2）編程問題。現(xiàn)在本科大部分都不再教授Fotran，盡管這個(gè)語言并不難學(xué)。

（3）力學(xué)基礎(chǔ)。需要具有彈性力學(xué)、有限元基礎(chǔ)，了解增量法和全量法的編程流程。

  

現(xiàn)在回過頭看，由于以上三點(diǎn)原因，導(dǎo)致我在初學(xué)的時(shí)候一頭扎進(jìn)這些問題里面，在很久以后再意識(shí)到。UMAT本質(zhì)上是子程序，什么是子程序？就是可以被主程序調(diào)用的一個(gè)函數(shù)。比如我們?cè)赑ython中def自定義的函數(shù)，比如在Matlab中定義的function。當(dāng)我明白了這一點(diǎn)，以前學(xué)習(xí)UMAT的困惑也都迎刃而解了。比如ABAQUS的子程序之所以是Fotran語言，除了Fotran本身計(jì)算效率高的原因外，很大程度上是因?yàn)锳BAQUS本身的求解器也是用該語言編寫的。

我們自己在研的求解器和有限元軟件使用Python編寫的，未來我們開放子程序權(quán)限的時(shí)候，也肯定是用Python語言了。

回到本構(gòu)。舉例，最簡(jiǎn)單的一維線彈性（彈簧）本構(gòu)：

F=K·x

我們就可以在UAMT中描述這個(gè)關(guān)系，當(dāng)然是采用矩陣進(jìn)行描述，告訴求解器應(yīng)力和應(yīng)變之間的關(guān)系，這樣求解器在求解的過程中就會(huì)按照我們描述的關(guān)系來進(jìn)行力學(xué)模擬。如果根據(jù)某個(gè)條件，判斷出材料在某個(gè)應(yīng)力狀態(tài)時(shí)，剛度發(fā)生了變化，在程序中加一個(gè)if語句，對(duì)彈性參數(shù)做修改即可。

因?yàn)橛羞@樣的特點(diǎn)，所以UMAT是我們研究復(fù)合材料力學(xué)行為最常用的手段之一。比如所謂的損傷折減研究，就是在迭代過程中，根據(jù)材料的應(yīng)力狀態(tài)判斷損傷的發(fā)生，然后對(duì)其剛度性能做出相應(yīng)的折減。因此UMAT是嵌入在求解器迭代中使用的。

對(duì)應(yīng)的有限元求解器基本步驟

1. 組建剛度矩陣K，abaqus自己處理
2. 載荷列陣F，abaqus自己處理
3. x=K^-1*F，所有節(jié)點(diǎn)的位移 abaqus自己處理
4. 根據(jù)x，求每個(gè)單元的應(yīng)變?cè)隽恐?abaqus自己處理
5. 根據(jù)單元應(yīng)變和單元?jiǎng)偠染仃嚕髴?yīng)力。UMAT處理

UMAT在迭代中發(fā)揮作用

F, F/10（總載荷和載荷增量）
Step1 F/10，stran0=0，Dstran1，stran1=stran0+Dstran1→UAMT計(jì)算stress1
Step2 F/10，F(xiàn)/10，stran1，Dstran2，stran2=stran1+Dstran2→UAMT計(jì)算stress2
......

UAMT程序結(jié)構(gòu)

UMAT子程序的結(jié)構(gòu)可以分為三個(gè)部分：

• 頭部。用來定義變量和變量大小，見下面的標(biāo)黃區(qū)域。
• 正文。描述本構(gòu)關(guān)系 。
• 結(jié)束，見下面紅色字體。

子程序基于Fortran語言編寫，形式如下：

      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)      !子程序的變量
C     
      INCLUDE 'ABA_PARAM.INC'                 !include 語句為浮點(diǎn)變量設(shè)置了適當(dāng)?shù)木?
C     
      CHARACTER*80 CMNAME                       !定義材料名稱最多占用的字符位數(shù)
      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),
     4    DSTRES(NTENS)                                 !定義變量的大小
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
.....................................代碼正文，定義應(yīng)力應(yīng)變關(guān)系
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        RETURN
        END

比較重要的兩個(gè)部分是子程序變量、變量大小的定義。有過其他代碼編寫經(jīng)驗(yàn)的話，特別是寫過自定義的函數(shù)代碼可以理解為：

• SUBROUTINE UMAT是函數(shù)名稱，后面括號(hào)跟的是函數(shù)輸入。
• 變量大小的定義是注釋性的設(shè)置，表明每個(gè)變量的數(shù)組大小。

下面一個(gè)重要的問題來了，既然子程序?qū)嵸|(zhì)上是一個(gè)自定義的function，即有著固定的傳入?yún)?shù)和返回參數(shù)。UMAT的返回值是那些呢？我們看官方文檔：

The following quantities must bedefined:Stress(STRESS), SDVs(STATEV), and material Jacobian(DDSDDE) 
The following variables may bedefined: Strain energy(SSE), plastic dissipation(SPD), and “creep” dissipation Suggested new (reduced) time increment(SCD)

從文檔中可以看出，UMAT的返回值分為兩類：

• 必須要計(jì)算的：STRESS、STATEV（如果定義了SDV）、DDSDDE；
• 按需返回的：SSE、SPD、SCD。

也就是說，在程序的正文部分，STRESS、STATEV（如果定義了SDV）、DDSDDE一定要更新。

上面我們看到那么多眼花繚亂的變量，這些變量具體是干啥的呢？看下面的表格，字體傾斜的變量是重要常用的變量。

STRESS（NTENS）	該增量步開始之前的應(yīng)力向量，在增量步結(jié)束之后，必須進(jìn)行更新。如果指定了初始應(yīng)力，則該向量在分析開始始將保持初始應(yīng)力。真實(shí)Cauchy應(yīng)力。	需要定義的變量，在所有分析情況下均適用。
STATEV（NSTATV）	求解過程中的狀態(tài)變量的存貯向量。在該增量步開始時(shí)，用來傳遞狀態(tài)變量，除非進(jìn)行了更新（采用USDEFL或UEXPAN）。在增量步結(jié)束時(shí)，STATEV更新為結(jié)束時(shí)刻的狀態(tài)變量值。STATEV的維數(shù)（NSTATV）由*DEPVAR決定。
DDSDDE（NTENS，NTENS）	Jacobian矩陣，即本構(gòu)關(guān)系矩陣。?σ/?ε。除非聲明采用非對(duì)稱方程求解，否則均為對(duì)稱矩陣DDSDDE(i, j)。
SSE、SPD、SCD	彈性應(yīng)變能、塑性耗能、徐（蠕）變耗能。在該增量步結(jié)束時(shí)進(jìn)行更新，并在下一增量步開始時(shí)進(jìn)行傳遞。這些能量參數(shù)對(duì)于求解結(jié)果不起作用，除非結(jié)果采用能量形式輸出。
RPL	該增量步結(jié)束時(shí)，由于材料的力學(xué)作工而產(chǎn)生的體積發(fā)熱量。	只用于完全耦合的溫度－應(yīng)力分析
DDSDDT（NTENS）	與溫度想對(duì)應(yīng)的應(yīng)力增量的變化量
DRPLDE（NTENS）	與應(yīng)變?cè)隽肯鄬?duì)應(yīng)的體積發(fā)熱量（RPL）的變化量
DRPLDT（NTENS）	與溫度相對(duì)應(yīng)的體積發(fā)熱量（RPL）的變化量
STRAIN（NTENS）	該增量步開始之前的總應(yīng)變向量。如果考慮了熱膨脹效應(yīng)，那么STRAIN僅為力學(xué)應(yīng)變（即已經(jīng)在總應(yīng)變中減去了熱膨脹得到的溫度應(yīng)變）。這些應(yīng)變?cè)谳敵鼋Y(jié)果中以“彈性”應(yīng)變給出。	信息傳遞變量
DSTRAIN（NTENS）	應(yīng)變?cè)隽肯蛄俊Ｈ绻紤]了熱膨脹應(yīng)變，則僅表示力學(xué)應(yīng)變?cè)隽俊?/em>
TIME（1） TIME（2）	當(dāng)前分析步開始時(shí)刻的，時(shí)間步的值。 當(dāng)前分析步開始時(shí)刻的，總時(shí)間的值。
DTIME	該增量步的時(shí)間增量
TEMP	當(dāng)前增量步開始時(shí)刻的溫度
DTEMP	該增量步的溫度增量
PREDEF	在當(dāng)前增量步開始時(shí)刻的，預(yù)定義的場(chǎng)變量（基于讀入的節(jié)點(diǎn)值）的內(nèi)插值向量。
DPRED	預(yù)定義的場(chǎng)變量的增量向量
CMNAME	用戶定義的材料名。由于ABAQUS內(nèi)部的一些給定材料是以“ABQ_”作為材料名，因此應(yīng)盡量不采用“ABQ_”作為CMNAME的名稱。
NDI	該點(diǎn)的直接應(yīng)力分量的個(gè)數(shù)
NSHR	該點(diǎn)的工程剪應(yīng)力分量的個(gè)數(shù)
NTENS	應(yīng)力或應(yīng)變向量的維數(shù)，等于NDI +NSHR。
NSTATV	求解過程中的狀態(tài)變量的個(gè)數(shù)，與材料類型匹配。
PROPS（NPROPS）	用戶定義的材料常數(shù)
NPROPS	用戶定義的材料常數(shù)的個(gè)數(shù)
COORDS	該點(diǎn)的坐標(biāo)向量。如果在當(dāng)前分析步中沒有考慮幾何非線性，COORDS就等于當(dāng)前坐標(biāo)系下的向量。否則，COORDS為最開始的坐標(biāo)向量	信息傳遞變量
DROT（3，3）	轉(zhuǎn)角增量矩陣。代表了剛體的基本坐標(biāo)系中的轉(zhuǎn)角增量（該基本坐標(biāo)系就是應(yīng)力、應(yīng)變向量存儲(chǔ)時(shí)的坐標(biāo)系）。用于用戶定義子程序中的向量或矢量狀態(tài)變量的轉(zhuǎn)角處理，而應(yīng)力及應(yīng)變向量在UMAT調(diào)用之前已經(jīng)進(jìn)行了轉(zhuǎn)角處理。在小位移分析中，該矩陣是一個(gè)單位矩陣；在大位移分析中，如果該材料點(diǎn)的基本坐標(biāo)系隨著材料坐標(biāo)系轉(zhuǎn)動(dòng)（如殼單元或采用了局部轉(zhuǎn)角坐標(biāo)時(shí)），該矩陣亦是一個(gè)單位矩陣。
PNEWDT	建議的新時(shí)間增量與原時(shí)間增量（DTIME）之間的比值大小。該變量允許用戶在ABAQUS/Standard中輸入自動(dòng)時(shí)間增量的計(jì)算法則（如果設(shè)置了自動(dòng)時(shí)間增量）。對(duì)于ABAQUS/Standard的準(zhǔn)靜態(tài)分析中的自動(dòng)時(shí)間增量（基于標(biāo)準(zhǔn)蠕變率積分技術(shù)），不允許在UMAT中控制。 在每一次調(diào)用UMAT前，PNEWDT被設(shè)置為一個(gè)足夠大的值。 如果沒有選擇自動(dòng)時(shí)間增量方法，大于1.0的PNEWDT值將被忽略，而起作用的僅是當(dāng)小于1.0的PNEWDT值時(shí)。	能夠更新的變量
CELENT	特征單元長(zhǎng)度。	信息傳遞變量
DFGRD0（3，3）	該增量步開始時(shí)刻的變形梯度向量。
DFGRD1（3，3）	該增量步結(jié)束時(shí)刻的變形梯度向量。如果在分析步中未考慮幾何非線性，則該向量為零。
NOEL	單元的個(gè)數(shù)
NPT	積分點(diǎn)的個(gè)數(shù)	信息傳遞變量
LAYER	層的個(gè)數(shù)（用于復(fù)合材料的殼單元，及層結(jié)構(gòu)固體單元）
KSPT	當(dāng)前層的截面點(diǎn)的個(gè)數(shù)
KSTEP	分析步的個(gè)數(shù)
KINC	增量步的個(gè)數(shù)

如何編寫UAMT

Fortran基礎(chǔ)

我們?cè)赨AMT中常用的語法有：

• 判斷語句

      IF ((STRESS(1) .GE. 0) .AND. (mod1 .GE. 1)) THEN
          
          STATEV(1)= 1     !纖維拉伸失效
          
      ELSEIF ((STRESS(1) .LT. 0) .AND. (mod2 .GE. 1)) THEN
          
          STATEV(2) = 1    !纖維壓縮失效
          
      END IF

•  循環(huán)語句

      DO I = 1, NTENS
         STRANT(I) = STRAN(I) + DSTRAN(I)
      END DO   

• 比較語法

.GT. 大于(>)
.LT. 小于(<)
.EQ. 等于(==)
還有
.GE. 大于等于(>=)
.LE. 小于等于(<=)
.NE. 不等于(!=)

• 特殊運(yùn)算

a的平方：a**2
自然常數(shù)的2次方：EXP(2)
自然對(duì)數(shù)是 log(x)
十為第的對(duì)數(shù)是 log10(x)

• 注釋

在任意位置使用  !  可以進(jìn)行注釋
在每行的第一個(gè)字符位，用 C 也可以注釋

自定義數(shù)組變量

在頭部的DIMENSION后可以自行定義自己需要的數(shù)組變量名稱和對(duì)應(yīng)的大小

      DIMENSION STRANT(6)
     DIMENSION CFULL(6,6),CDFULL(6,6),A(3,3),V(3,3)

自定義常數(shù)

可以按照如下形式，自己定義一些常數(shù)。

      PARAMETER (ZERO = 0.D0,ONE = 1.D0,TWO = 2.D0, HALF = 0.5D0)

如何調(diào)用UAMT

調(diào)用UAMT需要如下幾個(gè)步驟：

• 參數(shù)定義部分。
• Step設(shè)置部分。
• Job設(shè)置部分。

在參數(shù)定義模塊，>General，定義變量數(shù)目，可以多給但是不能少。>General>User Material，隨便給個(gè)數(shù)，我們?cè)赨MAT里面已經(jīng)定義好參數(shù)值了，此處輸入的值不會(huì)被調(diào)用。（下面的圖只是示意）

 

 

為了能讓損傷結(jié)果順利顯示，需要在Step>Field Output>State>SDV勾選。（下面的圖只是示意）

 

最后的調(diào)用，在Job模塊，調(diào)用子程序。

損傷通過SDV查看。

如何調(diào)試UAMT

提交計(jì)算后，如果子程序本身有語法錯(cuò)誤，軟件會(huì)自動(dòng)終止計(jì)算。并在安裝路徑專門存儲(chǔ)計(jì)算文件的文件夾下面，生成一個(gè)當(dāng)前工況的.log文件，里面會(huì)提升錯(cuò)誤代碼的位置和信息。根據(jù)這些信息我們可以進(jìn)行代碼修改。

除此之外，我們可以在Fortan將關(guān)注的變量值實(shí)時(shí)打印到文件中，幫助我們進(jìn)行問題定位。

如何配置子程序

前面說了，UMAT本質(zhì)是子程序，它實(shí)際上是需要能夠在ABAQUS傳入?yún)?shù)后，然后它完成計(jì)算再將指定變量返回。我們?cè)谝慌_(tái)干凈電腦上是無法直接運(yùn)行C語言代碼、MATLAB代碼或者Python代碼，因?yàn)槲覀儧]有安裝它們的運(yùn)行環(huán)境。這就是UMAT需要做環(huán)境配置的原因。

子程序的使用，需要安裝VS和FORTRAN，并且將ABAQUS與兩者進(jìn)行關(guān)聯(lián)，才能運(yùn)行。不同ABAQUS版本可能需要的VS和FORTRAN版本還有區(qū)別。這個(gè)網(wǎng)上有文件教程可以參考。

https://wenku.baidu.com/view/5766f34dd35abe23482fb4daa58da0116d171f2a.html

嫌麻煩的話，花個(gè)幾十塊找人裝。

未來設(shè)想

前面提到，未來我們會(huì)在自己的有限元軟件中也提供自定義接口。這個(gè)接口大概率是用Python，當(dāng)然這個(gè)不是重點(diǎn)。

重點(diǎn)是AI快速發(fā)展的當(dāng)下，不少軟件都打出了AI賦能的噱頭，具體問AI賦能了什么，又都說不清楚，最基礎(chǔ)的求解器目前仍然需要依賴穩(wěn)定精準(zhǔn)的數(shù)值算法。可以說無論AI如何發(fā)展，數(shù)值求解器仍然是不可替代的，甚至求解器自身就是AI工具的基礎(chǔ)。

那么AI在有限元中比較能發(fā)揮作用的地方在哪呢？我們目前的設(shè)想是，應(yīng)該優(yōu)先在本構(gòu)模型中發(fā)揮作用。即，通過用戶描述材料的力學(xué)行為特征，AI就可以自動(dòng)生成UMAT或者VUMAT子程序，供用戶使用。

  

當(dāng)然，除了AI自身增強(qiáng)對(duì)力學(xué)概念的理解，還需要數(shù)據(jù)支持。無論是各種彈塑性模型還是復(fù)合材料漸進(jìn)失效模型，里面各種系數(shù)實(shí)際上是需要針對(duì)不同的材料做試驗(yàn)測(cè)試才能獲取。尤其是多種多樣的復(fù)合材料。就我所了解，德國就有將各個(gè)試驗(yàn)室和研究單位的材料本構(gòu)數(shù)據(jù)統(tǒng)一管理的計(jì)劃，一旦把數(shù)據(jù)格式統(tǒng)一，隨著時(shí)間積累形成大數(shù)據(jù)庫。在AI技術(shù)的加持下，說不定靠寫材料本構(gòu)就能研究生畢業(yè)的情況，將一去不復(fù)返。與之相對(duì)的，做基礎(chǔ)試驗(yàn)或者尋找材料特性應(yīng)用場(chǎng)景或?qū)⒏又匾?/span>

 

歡迎關(guān)注“靜界有限元”

工作室面向在校學(xué)生、科研院所老師提供結(jié)構(gòu)有限元仿真（含二次開發(fā)）、流體力學(xué)仿真、算法開發(fā)、軟件開發(fā)服務(wù)。

最后，有相關(guān)需求歡迎通過公眾號(hào)聯(lián)系我們。

工眾浩：靜界有限元