ABAQUS 靜/動力隱式UEL 二次開發

有限元先生

瀏覽：3120 評論：5 收藏：9

概述：結合HHT時程積分法，推導了ABAQUS 靜/動力隱式算法中的UEL關鍵矩陣表達式，并將公式應用到自編CPS4/CPE4、C3D8 BBAR和C3D20用戶自定義單元中，計算結果均與ABAQUS自帶單元保持一致。其中，靜力計算中，關鍵矩陣AMARTX和RHS等可直接按照剛度矩陣和方程右端不平衡力輸出。動力隱式計算中，這兩者的輸出較為復雜，需要結合HHT時程積分法進行推導，將剛度矩陣、質量矩陣和阻尼矩陣依據LFLAGS數組的數值進行組合，RHS同樣需要進行推導計算，并以合適的方法將解相關的狀態量儲存在SVARS中，供后面的增量步調用。

關于UEL的程序設計，只支持靜力通用計算分析步的資料有很多，本帖子內容適不但適用于靜力通用，而且適用于動力隱式、頻率分析等分析步計算，采用模塊化程序設計，所有的矩陣求解均被封裝，調用方便，讀者可以針對自己的需求對相應函數進行改編即可，尤其是涉及到動力隱式計算部分，適用于任何運動方程的動力隱式求解，可直接移植使用。

----------------------------

（）UEL接口基本參數介紹

----------------------------

用戶自定義單元(USER DEFINED ELEEMNT, UEL)適合進階的工程師/學者使用，u、UEL的目的是實現一個單元的力學行為，即力-位移關系，這部分的功能需要的編程工作量和理論功底都比較高，涉及到的知識大致包括：

有限元理論：形函數插值、應力-應變關系，剛度矩陣組裝，質量矩陣求解，阻尼矩陣求解，

數值算法：高斯積分（全積分、縮減積分）、非線性方程組迭代求解流程（增量迭代，常剛度迭代等）、雅各比矩陣求解

彈/塑性力學：應力-應變關系，位移-應變關系

熟練的FORTRAN編程操作是基礎技能。

UEL接口參數眾多，包括自己本身的參數和與其他子程序聯合的參數，下面介紹我比較重要的幾個參數。

<>需要編程定義的參數

RHS（right hand side）：這個命名是從方程組的角度來的，顧名思義，他就是方程的右端量，其本質是：外力-內力，外力部分程序的編寫涉及到與其他子程序的聯合使用，包括DLOAD和UTRACLOAD等等，這部分內容是給用戶自定義單元施加復雜的廣義力，UEL接口為其提供了相應的參數，如JDLTYP、NDLOAD等等，目前沒有做過嘗試。至于內力的求解，在靜力線性計算中，數值上等于-KU（后面有這個公式的推導），即剛度與位移乘積負數。在動力隱式分析中，這個矩陣需要結合HHT時程積分法推導具體的表達式。

AMATRX：這個參數往往被認為是剛度，但并不這樣，他只有在靜力計算的時候才是剛度。他的具體取值依據分析類型的不同而不同，具體的表達形式不唯一。在動力隱式分析中，這個矩陣需要結合HHT時程積分法推導具體的表達式。

SVARS：這個參數官方文檔說是取決于結果的狀態量，他的具體意義和數據由我們自己確定，而我把它理解為一個小倉庫，可以存放我們的數據。需要注意的是里面的數據可以在不同的增量步之間傳輸，就是說，里面的數據會傳到下一個增量步，只要不更新他，他可以一直被使用。這個數據我們自己決定更新與否。

ENERGY：與用戶自定義單元相關的一些動量等能量，也是依據不同的分析類型而不同。

() CPS4/CPE4公式推導

等參單元中，母單元向笛卡爾坐標的轉換為：（注意：這里的箭頭標反了）

要再現位移場，可假設單元內部位移為坐標的函數，即：

將四個節點的位移和坐標代入上式，有：

如此，求解系數，并用節點位移表示單元內部位移：

上式即為單元節點位移與坐標的插值關系。考慮到等參單元，即位移的插值函數與母單元-笛卡爾單元變換采用相同插值函數。則等參變化的插值函數為：

上式的矩陣表達為：

其中：

對形函數求導，得形函數對母單元坐標的導數：

由彈性力學幾何方程，應變矩陣為：

合并為：

JOCABIN矩陣為：

其中：

則JOCABIN的具體表達式為：

剛度矩陣的求解：

質量矩陣分為兩種，集中質量矩陣和協調質量矩陣，下面給出的是一致質量矩陣（又叫協調質量矩陣）矩陣，該種質量矩陣因其采用了與位移插值相同的插值函數而被稱為一致質量矩陣。但是協調質量矩陣為滿陣，不利于數值計算，所以有限元中一般采用集中質量矩陣，即只有對角元素不為零。集中質量矩陣通常由協調質量矩陣處理而來，比如將協調質量矩陣每一行的質量集中到對角元素，然后將非對角元素數值置為零。

阻尼矩陣，采用比例阻尼：

UEL主程序采用常剛度法設計，剛度矩陣和質量矩陣只在第一個增量步計算一次，然后儲存在變量SVARS中，后面一直到分析結束之間的所有增量步不再計算，均直接讀取SVARS數組中的數據并輸出會ABAQUS主程序。該主程序將各個關鍵量的求解分為不同的子程序，這些子程序基本是所有用戶自定義單元都會涉及的，因此該主程序的可移植性很強，在線彈性問題中，即可用常剛度迭代法求解的運動方程中，讀者完全可以按照這個直接編寫自己的程序。

    SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
     1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,
     2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,
     3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD)
C
      INCLUDE 'ABA_PARAM.INC'

C
      DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*),
     1 SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL),
     2 DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*),
     3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),
     4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)
      double precision DMATX(3,3)
      double precision KK(2*nnode,2*nnode),MM(2*nnode,2*nnode)  
      KK=0.D0   
      MM=0.D0       
!      Print work state on screen and ANA_STAT.txt
      CALL PRINTSTATE(JPROPS,LFLAGS,DTIME,JELEM) 
      IF(KSTEP.EQ.1.AND.KINC.EQ.1) THEN           
!        material martix   
        CALL DMATERIAL(DMatx,PROPS)    
!        stiffness martix     
        call KKmartix(KK,coords,DMatx,mcrd,nnode,jelem)  
!        mass martix  
        call MMmartix(MM,coords,props,mcrd,nnode,jelem,dmatx)
!        Store pre-calculated stiff/mass matrices
        CALL STOREMATRICES (SVARS,NSVARS,KK,MM,NDOFEL)
      else
        CALL READMATRICES(SVARS,NSVARS,KK,MM,NDOFEL)
      endif 
!      Output required variables
      CALL OUTPUTVARIABLE(RHS,AMATRX,SVARS,PROPS,ENERGY,U,V,A,
     1 LFLAGS,DTIME,NDOFEL,NRHS,NSVARS,MLVARX,JELEM,PARAMS,KK,MM)
      RETURN
      END

付費文件包含但不限于以下內容：

-------------------------------------------------

（1）CPE4/CPS4單元UEL，適用于靜/動力隱式計算

（2）C3D8 BBAR單元UEL，，適用于靜/動力隱式計算，包含UMAT給應力應變可視化部分

（3）C3D20 單元UEL，適用于靜/動力隱式計算

（4）WORD文檔，推導、講解動力隱式計算中的AMARTX、RHS

贈送內容：

（5）收集的一些UEL-UMAT程序

（6）一對一輔導，包括理論、公式和編程等......

-------------------------------------------------

ps：付款之前咨詢我一下，技術鄰的錢提現太慢了，我孫子輩都花不上！

以下內容為付費內容，請購買后觀看