UEL單元開發（2）——一維桿單元（內含大福利！）

易公子

2022年5月20日 16:27

瀏覽：3031 收藏：3

木木做單元開發也有一段日子了，在編寫XFEM的單元時，遇到了瓶頸，打算先告一段落，將自己的一些經驗分享出來給大家，希望幫助大家快速入門，在以后的推文中也會以單元開發為主，就像之前講解Umat一樣，用自己通俗易懂的語言講述下來。

本篇推文由兩部分組成：「其一」，以「一維桿單元」UEL子程序為主，從理論到代碼Step-by-Step向大家講解，并于Abaqus中標準單元做對比，驗證代碼的正確性；「其二」，木木受邀于北鯤云超算平臺做主講人，將在5月24日晚上七點于B站進行大約一個小時的講演匯報，帶著大家熟悉Abaqus的基本操作，基于Abaqus的擴展有限元仿真以及如何正確使用超算平臺，直播過程中還有大量驚喜禮品、豐富算力金進行抽獎~歡迎大家積極參與。

一維桿單元UEL子程序

一維桿單元子程序的代碼應該是UEL中最簡單的了吧，就好像線彈性Umat那樣，囊括了一些UEL的基本思想，可對比之前的二維彈簧單元子程序來學習。

理論講解

彈性力學中利用「最小勢能原理」（或虛功原理）來得到單元剛度方程，具體操作可查看：

?

《有限元基礎教程》——曾攀

《北航彈性力學博后講透彈性力學》——@兵心依舊022（B站Up）

?

單元剛度矩陣：

其中：為楊氏模量，為桿的橫截面積，為桿的長度。UEL需要的「核心」就是求得單元剛度矩陣，「殘余力RHS」可以認為是0（Cohesive單元除外），以后會給大家普及國際上流行的PPR內聚力單元。

代碼講解

篇幅原因，INP文件可在后臺回復“桿單元”，即可獲得，這里僅給出「關鍵行」：

*User element, nodes=2, type=U1001, properties=2, coordinates=2, variables=4
1,2,6
*Element, type=U1001
 1,1,2
*Elset, elset=Set-1
 1
*Uel property, elset=Set-1
 2.0E11,  0.02

程序代碼如下：

      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(*)

      INTEGER I,J
C     獲取材料屬性
      Y_M = PROPS(1)
      AREA = PROPS(2)
C     獲取節點坐標計算單元長度
      LEN_EL = ABS(COORDS(1,2)-COORDS(1,1))
      write(6,*) LEN_EL
C     構造剛度矩陣
      K_EL = Y_M*AREA/LEN_EL
      AMATRX(1,1) = K_EL
      AMATRX(1,2) = -K_EL
      AMATRX(2,1) = -K_EL
      AMATRX(2,2) = K_EL
C     求RHS

      DO I = 1,NDOFEL
          RHS(I,1) = 0
          DO J = 1,NDOFEL
              RHS(I,1) = RHS(I,1) - AMATRX(I,J)*U(J)
          ENDDO
      ENDDO
      
      RETURN
      END