C3D20 UEL用户自定义单元开发（1）_Abaqus_通用

C3D20 UEL用户自定义单元开发（1）

零

概述

开发了适用于静力通用计算的三维二十节点（C3D20）的用户自定义单元，在挖孔悬臂梁受剪切荷载算例中，位移计算结果与ABAQUS自带单元保持一致。对比刚度矩阵，与abaqus保持一致。

帖子分为两部分

->理论部分-C3D20 UEL用户自定义单元开发（1）结果部分-C3D20 UEL用户自定义单元开发（2）

壹

二理论与程序设计

C3D20单元节点排布示意图

形函数为

求解刚度矩阵主程序

      subroutine KKmartix(KK,coords,DMatx,mcrd,nnode,jelem)          INCLUDE 'ABA_PARAM.INC'      double precision coords(mcrd,nnode),dmatx(6,6)      double precision integerpoint(1,3),wgt(1,3)      double precision KK(3*nnode,3*nnode)      double precision B(6,3*NNODE),J(3,3)      double precision detJ,kesai,eta,zeta,Coff      double precision npxy(3,NNODE)     !积分点，二阶单元，设置三个积分点            integerpoint(1,1)= 0.D0      integerpoint(1,2)= 0.7745966692D0      integerpoint(1,3)=-0.7745966692D0!积分权重                  wgt(1,1)=0.8888888889D0      wgt(1,2)=0.5555555556D0      wgt(1,3)=0.5555555556D0                  B=0.d0      J=0.d0      detJ=0.d0      kesai=0.d0      eta=0.d0      zeta=0.d0      Coff=0.d0      npxy=0.d0                               do i=1,3        do ii=1,3          do iii=1,3           kesai=integerpoint(1,i)            eta  =integerpoint(1,ii)            zeta =integerpoint(1,iii)            call Shapefunction(npxy,detj,kesai,eta,zeta,     1 coords,mcrd,nnode,J)                      call CalBmartix(B,npxy)             write(6,*)"kk-detj=",detj                         COFF=detJ*wgt(1,i)*wgt(1,ii)*wgt(1,iii)            KK=KK+matmul(matmul(transpose(B),DMatx),B)*COFF          enddo          enddo          enddo               return       end

形函数对物理坐标（x、y、z）的导数

       subroutine Shapefunction(npxy,detj,kesai,eta,zeta,     1 coords,mcrd,nnode,J)       INCLUDE 'ABA_PARAM.INC'       integer mcrd,nnode       double precision kesai,eta,zeta,detJ       double precision coords(mcrd,nnode),npxy(3,nnode)       double precision J(3,3),NPGHR(3,nnode)       double precision invJ(3,3)                            j=0.d0       jj=0.d0       invJ=0.d0                                 call CalNpLocalCoord(NPGHR,kesai,eta,zeta)                J=matmul(NPGHR,transpose(coords))                  Call CaldetJ(detJ,J)              call CalInvjocabin(invJ,detj,j)       npxy=matmul(invJ,NPGHR)                  return       end

B矩阵

      subroutine CalBmartix(B,npxy)      INCLUDE 'ABA_PARAM.INC'      double precision npxy(3,20),B(6,60)              B=0.d0            do i=1,20        B(1,3*i-2)=npxy(1,i)        B(2,3*i-1)=npxy(2,i)        B(3,3*i)=npxy(3,i)                B(4,3*i-2)=npxy(2,i)        B(4,3*i-1)=npxy(1,i)                B(5,3*i-1)=npxy(3,i)        B(5,3*i)=npxy(2,i)                 B(6,3*i-2)=npxy(3,i)        B(6,3*i)=npxy(1,i)                       enddo          return       end

形函数对局部坐标（G、H、R）的导数

   SUBROUTINE CalNpLocalCoord(NPGHR,KESAI,ETA,ZETA)      INCLUDE 'ABA_PARAM.INC'            DOUBLE PRECISION NPGHR(3,20)      DOUBLE PRECISION KESAI,ETA,ZETA      DOUBLE PRECISION PG,PH,PR              DO i=1,20       CALL NPG(PG,KESAI,ETA,ZETA,I)       CALL NPH(PH,KESAI,ETA,ZETA,I)       CALL NPR(PR,KESAI,ETA,ZETA,I)              NPGHR(1,i)=PG       NPGHR(2,i)=PH       NPGHR(3,i)=PR      ENDDO               RETURN      END-------------------------!形函数对局部坐标G的导数      SUBROUTINE NPG(PG,G,H,R,NINDEX)      DOUBLE PRECISION G,H,R,PG      INTEGER NINDEX      IF   (NINDEX.EQ.1)THEN          PG=-0.125*( -1)*(1-H)*(1-R)*(2+G+H+R)     1     -0.125*(1-G)*(1-H)*(1-R)*(  1)          ELSEIF(NINDEX.EQ.2)THEN        PG=-0.125*(  1)*(1-H)*(1-R)*(2-G+H+R)     1     -0.125*(1+G)*(1-H)*(1-R)*( -1)        ELSEIF(NINDEX.EQ.3)THEN        PG=-0.125*(  1)*(1+H)*(1-R)*(2-G-H+R)     1     -0.125*(1+G)*(1+H)*(1-R)*( -1)      ELSEIF(NINDEX.EQ.4)THEN        PG=-0.125*( -1)*(1+H)*(1-R)*(2+G-H+R)     1     -0.125*(1-G)*(1+H)*(1-R)*(  1)      ELSEIF(NINDEX.EQ.5)THEN        PG=-0.125*( -1)*(1-H)*(1+R)*(2+G+H-R)     1     -0.125*(1-G)*(1-H)*(1+R)*(  1)      ELSEIF(NINDEX.EQ.6)THEN        PG=-0.125*(  1)*(1-H)*(1+R)*(2-G+H-R)     1     -0.125*(1+G)*(1-H)*(1+R)*( -1)      ELSEIF(NINDEX.EQ.7)THEN        PG=-0.125*(  1)*(1+H)*(1+R)*(2-G-H-R)     1     -0.125*(1+G)*(1+H)*(1+R)*( -1)      ELSEIF(NINDEX.EQ.8)THEN        PG=-0.125*( -1)*(1+H)*(1+R)*(2+G-H-R)     1     -0.125*(1-G)*(1+H)*(1+R)*(  1)      ELSEIF(NINDEX.EQ.9)THEN            PG=0.25*( -1)*(1+G)*(1-H)*(1-R)     1    +0.25*(1-G)*(  1)*(1-H)*(1-R)      ELSEIF(NINDEX.EQ.10)THEN        PG=0.25*(1-H)*(1+H)*( 1)*(1-R)      ELSEIF(NINDEX.EQ.11)THEN        PG=0.25*( -1)*(1+G)*(1+H)*(1-R)     1    +0.25*(1-G)*(  1)*(1+H)*(1-R)      ELSEIF(NINDEX.EQ.12)THEN        PG=0.25*(1-H)*(1+H)*(-1)*(1-R)      ELSEIF(NINDEX.EQ.13)THEN        PG=0.25*( -1)*(1+G)*(1-H)*(1+R)     1    +0.25*(1-G)*(  1)*(1-H)*(1+R)      ELSEIF(NINDEX.EQ.14)THEN        PG=0.25*(1-H)*(1+H)*( 1)*(1+R)      ELSEIF(NINDEX.EQ.15)THEN        PG=0.25*( -1)*(1+G)*(1+H)*(1+R)     1    +0.25*(1-G)*(  1)*(1+H)*(1+R)      ELSEIF(NINDEX.EQ.16)THEN        PG=0.25*(1-H)*(1+H)*(-1)*(1+R)      ELSEIF(NINDEX.EQ.17)THEN         PG=0.25*(1-R)*(1+R)*(-1)*(1-H)      ELSEIF(NINDEX.EQ.18)THEN        PG=0.25*(1-R)*(1+R)*( 1)*(1-H)      ELSEIF(NINDEX.EQ.19)THEN        PG=0.25*(1-R)*(1+R)*( 1)*(1+H)      ELSEIF(NINDEX.EQ.20)THEN           PG=0.25*(1-R)*(1+R)*(-1)*(1+H)      ENDIF      RETURN      END!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!形函数对局部坐标H的导数      SUBROUTINE NPH(PH,G,H,R,NINDEX)      DOUBLE PRECISION G,H,R,PH      INTEGER NINDEX      IF    (NINDEX.EQ.1)THEN          PH=-0.125*(1-G)*( -1)*(1-R)*(2+G+H+R)     1     -0.125*(1-G)*(1-H)*(1-R)*(    1  )      ELSEIF(NINDEX.EQ.2)THEN        PH=-0.125*(1+G)*( -1)*(1-R)*(2-G+H+R)     1     -0.125*(1+G)*(1-H)*(1-R)*(    1  )      ELSEIF(NINDEX.EQ.3)THEN        PH=-0.125*(1+G)*(  1)*(1-R)*(2-G-H+R)     1     -0.125*(1+G)*(1+H)*(1-R)*(   -1  )      ELSEIF(NINDEX.EQ.4)THEN        PH=-0.125*(1-G)*(  1)*(1-R)*(2+G-H+R)     1     -0.125*(1-G)*(1+H)*(1-R)*(   -1  )      ELSEIF(NINDEX.EQ.5)THEN        PH=-0.125*(1-G)*( -1)*(1+R)*(2+G+H-R)     1     -0.125*(1-G)*(1-H)*(1+R)*(    1  )      ELSEIF(NINDEX.EQ.6)THEN        PH=-0.125*(1+G)*( -1)*(1+R)*(2-G+H-R)     1     -0.125*(1+G)*(1-H)*(1+R)*(    1  )      ELSEIF(NINDEX.EQ.7)THEN        PH=-0.125*(1+G)*(  1)*(1+R)*(2-G-H-R)     1     -0.125*(1+G)*(1+H)*(1+R)*(   -1  )      ELSEIF(NINDEX.EQ.8)THEN        PH=-0.125*(1-G)*(  1)*(1+R)*(2+G-H-R)     1     -0.125*(1-G)*(1+H)*(1+R)*(   -1  )      ELSEIF(NINDEX.EQ.9)THEN            PH=0.25*(1-G)*(1+G)*(-1)*(1-R)      ELSEIF(NINDEX.EQ.10)THEN        PH=0.25*( -1)*(1+H)*(1+G)*(1-R)     1    +0.25*(1-H)*(  1)*(1+G)*(1-R)      ELSEIF(NINDEX.EQ.11)THEN        PH=0.25*(1-G)*(1+G)*( 1)*(1-R)      ELSEIF(NINDEX.EQ.12)THEN        PH=0.25*( -1)*(1+H)*(1-G)*(1-R)     1    +0.25*(1-H)*(  1)*(1-G)*(1-R)      ELSEIF(NINDEX.EQ.13)THEN        PH=0.25*(1-G)*(1+G)*(-1)*(1+R)      ELSEIF(NINDEX.EQ.14)THEN        PH=0.25*( -1)*(1+H)*(1+G)*(1+R)     1    +0.25*(1-H)*(  1)*(1+G)*(1+R)      ELSEIF(NINDEX.EQ.15)THEN        PH=0.25*(1-G)*(1+G)*( 1)*(1+R)      ELSEIF(NINDEX.EQ.16)THEN        PH=0.25*( -1)*(1+H)*(1-G)*(1+R)     1    +0.25*(1-H)*(  1)*(1-G)*(1+R)      ELSEIF(NINDEX.EQ.17)THEN         PH=0.25*(1-R)*(1+R)*(1-G)*(-1)      ELSEIF(NINDEX.EQ.18)THEN        PH=0.25*(1-R)*(1+R)*(1+G)*(-1)      ELSEIF(NINDEX.EQ.19)THEN        PH=0.25*(1-R)*(1+R)*(1+G)*( 1)      ELSEIF(NINDEX.EQ.20)THEN           PH=0.25*(1-R)*(1+R)*(1-G)*( 1)      ENDIF      RETURN      END!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!形函数对局部坐标R的导数      SUBROUTINE NPR(PR,G,H,R,NINDEX)      DOUBLE PRECISION G,H,R,PR      INTEGER NINDEX      IF    (NINDEX.EQ.1)THEN          PR=-0.125*(1-G)*(1-H)*( -1)*(2+G+H+R)     1     -0.125*(1-G)*(1-H)*(1-R)*(      1)      ELSEIF(NINDEX.EQ.2)THEN        PR=-0.125*(1+G)*(1-H)*( -1)*(2-G+H+R)     1     -0.125*(1+G)*(1-H)*(1-R)*(      1)      ELSEIF(NINDEX.EQ.3)THEN        PR=-0.125*(1+G)*(1+H)*( -1)*(2-G-H+R)     1     -0.125*(1+G)*(1+H)*(1-R)*(      1)      ELSEIF(NINDEX.EQ.4)THEN        PR=-0.125*(1-G)*(1+H)*( -1)*(2+G-H+R)     1     -0.125*(1-G)*(1+H)*(1-R)*(      1)      ELSEIF(NINDEX.EQ.5)THEN        PR=-0.125*(1-G)*(1-H)*(  1)*(2+G+H-R)     1     -0.125*(1-G)*(1-H)*(1+R)*(     -1)      ELSEIF(NINDEX.EQ.6)THEN        PR=-0.125*(1+G)*(1-H)*(  1)*(2-G+H-R)     1     -0.125*(1+G)*(1-H)*(1+R)*(     -1)      ELSEIF(NINDEX.EQ.7)THEN        PR=-0.125*(1+G)*(1+H)*(  1)*(2-G-H-R)     1     -0.125*(1+G)*(1+H)*(1+R)*(     -1)      ELSEIF(NINDEX.EQ.8)THEN        PR=-0.125*(1-G)*(1+H)*(  1)*(2+G-H-R)     1     -0.125*(1-G)*(1+H)*(1+R)*(     -1)      ELSEIF(NINDEX.EQ.9)THEN            PR=0.25*(1-G)*(1+G)*(1-H)*(-1)      ELSEIF(NINDEX.EQ.10)THEN        PR=0.25*(1-H)*(1+H)*(1+G)*(-1)      ELSEIF(NINDEX.EQ.11)THEN        PR=0.25*(1-G)*(1+G)*(1+H)*(-1)      ELSEIF(NINDEX.EQ.12)THEN        PR=0.25*(1-H)*(1+H)*(1-G)*(-1)      ELSEIF(NINDEX.EQ.13)THEN        PR=0.25*(1-G)*(1+G)*(1-H)*( 1)      ELSEIF(NINDEX.EQ.14)THEN        PR=0.25*(1-H)*(1+H)*(1+G)*( 1)      ELSEIF(NINDEX.EQ.15)THEN        PR=0.25*(1-G)*(1+G)*(1+H)*( 1)      ELSEIF(NINDEX.EQ.16)THEN        PR=0.25*(1-H)*(1+H)*(1-G)*( 1)      ELSEIF(NINDEX.EQ.17)THEN         PR=0.25*( -1)*(1+R)*(1-G)*(1-H)     1    +0.25*(1-R)*(  1)*(1-G)*(1-H)      ELSEIF(NINDEX.EQ.18)THEN        PR=0.25*( -1)*(1+R)*(1+G)*(1-H)     1    +0.25*(1-R)*(  1)*(1+G)*(1-H)      ELSEIF(NINDEX.EQ.19)THEN        PR=0.25*( -1)*(1+R)*(1+G)*(1+H)     1    +0.25*(1-R)*(  1)*(1+G)*(1+H)      ELSEIF(NINDEX.EQ.20)THEN           PR=0.25*( -1)*(1+R)*(1-G)*(1+H)     1    +0.25*(1-R)*(  1)*(1-G)*(1+H)      ENDIF      RETURN      END

来源：有限元先生

C3D8 BBAR UEL用户自定义单元二次开发（1）

零概述开发了适用于静力通用、频率分析和动力隐式（固定增量步长和自适应增量步长均可）的三维八节点线性UEL，即ABAQUS自带的C3D8单元，该UEL考虑了B-BAR修正，避免剪切锁死。采用编写的UEL，分别设置了静力通用分析步、频率分析和动力隐式分析步，将计算结果与ABAQUS对比，位移、速度和加速度与ABAQUS均保持一致，说明该UEL复现了一部分C3D8单元的计算功能。帖子分为两部分->理论部分-C3D8 BBAR UEL用户自定义单元开发（1）结果部分-C3D8 BBAR UEL用户自定义单元开发（2）壹二理论与程序设计 C3D20单元节点排布示意图一阶形函数为求解质量矩阵主程序 subroutine MMmartix(MM,coords,props,mcrd,nnode,jelem) INCLUDE 'ABA_PARAM.INC' integer mcrd,nnode,jelem double precision MM(3*nnode,3*nnode),coords(mcrd,nnode) double precision props(*),jj(3,8),J(3,3) double precision integerpoint(1,2),NN(3,24),NPXY(3,8) double precision kesai,eta,zeta,density double precision detj,Coff,ELVOL MM=0.D0 jj=0.D0!积分权重均为1，不设置数组储存积分权重 integerpoint(1,1)=-0.577350269189626D0 integerpoint(1,2)= 0.5777350269189626D0 NN=0.D0 j=0.D0 kesai=0.D0 eta=0.D0 zeta=0.D0 density=props(3) detj=0.D0 Coff=0.D0 do i=1,2 do ii=1,2 do iii=1,2 kesai=integerpoint(1,i) eta =integerpoint(1,ii) zeta =integerpoint(1,iii) call CalNpLocalCoord(jj,kesai,eta,zeta) J=matmul(jj,transpose(coords)) Call CaldetJ(detJ,J) call CalNN(NN,kesai,eta,zeta) MM=MM+(matmul(transpose(NN),NN))*detJ*density ELVOL=ELVOL+detJ enddo enddo enddo CALL ConMass(MM,nnode) return end ! Contrate mass martix subroutine ConMass(MM,nnode) INCLUDE 'ABA_PARAM.INC' integer nnode double precision MM(3*nnode,3*nnode) do i=1,3*nnode mm(i,i)=sum(mm(i,1:3*nnode)) do ii=1,3*nnode if (i.ne.ii)then MM(i,ii)=0.D0 endif enddo enddo return end 计算B矩阵 subroutine CalBmartix(B,npxy) INCLUDE 'ABA_PARAM.INC' double precision npxy(3,8),B(6,24) B=0.d0 do i=1,8 B(1,3*i-2)=npxy(1,i) B(2,3*i-1)=npxy(2,i) B(3,3*i )=npxy(3,i) B(4,3*i-2)=npxy(2,i) B(4,3*i-1)=npxy(1,i) B(4,3*i )=0.D0 B(5,3*i-2)=0.D0 B(5,3*i-1)=npxy(3,i) B(5,3*i )=npxy(2,i) B(6,3*i-2)=npxy(3,i) B(6,3*i-1)=0.D0 B(6,3*i )=npxy(1,i) enddo return end 计算形函数 subroutine Shapefunction(npxy,detj,kesai,eta,zeta, 1 coords,mcrd,nnode,J) INCLUDE 'ABA_PARAM.INC' integer mcrd,nnode double precision kesai,eta,zeta,detJ double precision coords(mcrd,nnode),npxy(3,nnode) double precision J(3,3),jj(3,nnode),intpcoord(3,8) double precision invJ(3,3) j=0.d0 jj=0.d0 intpcoord=0.d0 invJ=0.d0 call CalNpLocalCoord(jj,kesai,eta,zeta) J=matmul(jj,transpose(coords)) Call CaldetJ(detJ,J) call CalInvjocabin(invJ,detj,j) npxy=matmul(invJ,jj) return end!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine CalInvjocabin(invJ,detj,j) INCLUDE 'ABA_PARAM.INC' double precision invJ(3,3),j(3,3) double precision detJ,DETInv DETInv=1.D0/detj invj(1,1)=DETInv*(j(2,2)*j(3,3)-j(2,3)*j(3,2)) invj(2,1)=DETInv*(-j(2,1)*j(3,3)+j(2,3)*j(3,1)) invj(3,1)=DETInv*(j(2,1)*j(3,2)-j(2,2)*j(3,1)) invj(1,2)=DETInv*(-j(1,2)*j(3,3)+j(1,3)*j(3,2)) invj(2,2)=DETInv*(j(1,1)*j(3,3)-j(1,3)*j(3,1)) invj(3,2)=DETInv*(-j(1,1)*j(3,2)+j(1,2)*j(3,1)) invj(1,3)=DETInv*(j(1,2)*j(2,3)-j(1,3)*j(2,2)) invj(2,3)=DETInv*(-j(1,1)*j(2,3)+j(1,3)*j(2,1)) invj(3,3)=DETInv*(j(1,1)*j(2,2)-j(1,2)*j(2,1)) return end!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !系数 subroutine CalNpLocalCoord(jj,kesai,eta,zeta) INCLUDE 'ABA_PARAM.INC' double precision jj(3,8) double precision kesai,eta,zeta double precision LoCalCoord(8,3) double precision xp1,xp2,xp3! local coord Call CaLoCalCoord(LoCalCoord) do i=1,8 xp1=LoCalCoord(i,1) xp2=LoCalCoord(i,2) xp3=LoCalCoord(i,3) jj(1,i)=xp1*0.125*(1+xp2*eta )*(1+xp3*zeta) jj(2,i)=xp2*0.125*(1+xp1*kesai)*(1+xp3*zeta) jj(3,i)=xp3*0.125*(1+xp1*kesai)*(1+xp2*eta ) enddo return end!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine CaLoCalCoord(LoCalCoord) INCLUDE 'ABA_PARAM.INC' double precision LoCalCoord(8,3) LoCalCoord=0.d0 LoCalCoord(1,1)=-1 LoCalCoord(2,1)= 1 LoCalCoord(3,1)= 1 LoCalCoord(4,1)=-1 LoCalCoord(5,1)=-1 LoCalCoord(6,1)= 1 LoCalCoord(7,1)= 1 LoCalCoord(8,1)=-1 LoCalCoord(1,2)=-1 LoCalCoord(2,2)=-1 LoCalCoord(3,2)= 1 LoCalCoord(4,2)= 1 LoCalCoord(5,2)=-1 LoCalCoord(6,2)=-1 LoCalCoord(7,2)= 1 LoCalCoord(8,2)= 1 LoCalCoord(1,3)=-1 LoCalCoord(2,3)=-1 LoCalCoord(3,3)=-1 LoCalCoord(4,3)=-1 LoCalCoord(5,3)= 1 LoCalCoord(6,3)= 1 LoCalCoord(7,3)= 1 LoCalCoord(8,3)= 1 return end!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine CaldetJ(detJ,J) INCLUDE 'ABA_PARAM.INC' double precision detJ double precision j(3,3) detJ=0.D0 detJ=j(1,1)*j(2,2)*j(3,3)+j(2,1)*j(3,2)*j(1,3) * +j(1,2)*j(2,3)*j(3,1)-j(1,3)*j(2,2)*j(3,1) * -j(3,2)*j(2,3)*j(1,1)-j(1,2)*j(2,1)*j(3,3) return end ·私信获取全部程序· 来源：有限元先生