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单元的计算功能。帖子分为两部分->理论部分-C3D8BBARUEL用户自定义单元开发（1）结果部分-C3D8BBARUEL用户自定义单元开发（2）壹二理论与程序设计C3D20单元节点排布示意图一阶形函数为求解质量矩阵主程序subroutineMMmartix(MM,coords,props,mcrd,nnode,jelem)INCLUDE'ABA_PARAM.INC'integermcrd,nnode,jelemdoubleprecisionMM(3*nnode,3*nnode),coords(mcrd,nnode)doubleprecisionprops(*),jj(3,8),J(3,3)doubleprecisionintegerpoint(1,2),NN(3,24),NPXY(3,8)doubleprecisionkesai,eta,zeta,densitydoubleprecisiondetj,Coff,ELVOLMM=0.D0jj=0.D0!积分权重均为1，不设置数组储存积分权重integerpoint(1,1)=-0.577350269189626D0integerpoint(1,2)=0.5777350269189626D0NN=0.D0j=0.D0kesai=0.D0eta=0.D0zeta=0.D0density=props(3)detj=0.D0Coff=0.D0doi=1,2doii=1,2doiii=1,2kesai=integerpoint(1,i)eta=integerpoint(1,ii)zeta=integerpoint(1,iii)callCalNpLocalCoord(jj,kesai,eta,zeta)J=matmul(jj,transpose(coords))CallCaldetJ(detJ,J)callCalNN(NN,kesai,eta,zeta)MM=MM+(matmul(transpose(NN),NN))*detJ*densityELVOL=ELVOL+detJenddoenddoenddoCALLConMass(MM,nnode)returnend!ContratemassmartixsubroutineConMass(MM,nnode)INCLUDE'ABA_PARAM.INC'integernnodedoubleprecisionMM(3*nnode,3*nnode)doi=1,3*nnodemm(i,i)=sum(mm(i,1:3*nnode))doii=1,3*nnodeif(i.ne.ii)thenMM(i,ii)=0.D0endifenddoenddoreturnend计算B矩阵subroutineCalBmartix(B,npxy)INCLUDE'ABA_PARAM.INC'doubleprecisionnpxy(3,8),B(6,24)B=0.d0doi=1,8B(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.D0B(5,3*i-2)=0.D0B(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.D0B(6,3*i)=npxy(1,i)enddoreturnend计算形函数subroutineShapefunction(npxy,detj,kesai,eta,zeta,1coords,mcrd,nnode,J)INCLUDE'ABA_PARAM.INC'integermcrd,nnodedoubleprecisionkesai,eta,zeta,detJdoubleprecisioncoords(mcrd,nnode),npxy(3,nnode)doubleprecisionJ(3,3),jj(3,nnode),intpcoord(3,8)doubleprecisioninvJ(3,3)j=0.d0jj=0.d0intpcoord=0.d0invJ=0.d0callCalNpLocalCoord(jj,kesai,eta,zeta)J=matmul(jj,transpose(coords))CallCaldetJ(detJ,J)callCalInvjocabin(invJ,detj,j)npxy=matmul(invJ,jj)returnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutineCalInvjocabin(invJ,detj,j)INCLUDE'ABA_PARAM.INC'doubleprecisioninvJ(3,3),j(3,3)doubleprecisiondetJ,DETInvDETInv=1.D0/detjinvj(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))returnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!系数subroutineCalNpLocalCoord(jj,kesai,eta,zeta)INCLUDE'ABA_PARAM.INC'doubleprecisionjj(3,8)doubleprecisionkesai,eta,zetadoubleprecisionLoCalCoord(8,3)doubleprecisionxp1,xp2,xp3!localcoordCallCaLoCalCoord(LoCalCoord)doi=1,8xp1=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)enddoreturnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutineCaLoCalCoord(LoCalCoord)INCLUDE'ABA_PARAM.INC'doubleprecisionLoCalCoord(8,3)LoCalCoord=0.d0LoCalCoord(1,1)=-1LoCalCoord(2,1)=1LoCalCoord(3,1)=1LoCalCoord(4,1)=-1LoCalCoord(5,1)=-1LoCalCoord(6,1)=1LoCalCoord(7,1)=1LoCalCoord(8,1)=-1LoCalCoord(1,2)=-1LoCalCoord(2,2)=-1LoCalCoord(3,2)=1LoCalCoord(4,2)=1LoCalCoord(5,2)=-1LoCalCoord(6,2)=-1LoCalCoord(7,2)=1LoCalCoord(8,2)=1LoCalCoord(1,3)=-1LoCalCoord(2,3)=-1LoCalCoord(3,3)=-1LoCalCoord(4,3)=-1LoCalCoord(5,3)=1LoCalCoord(6,3)=1LoCalCoord(7,3)=1LoCalCoord(8,3)=1returnend!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutineCaldetJ(detJ,J)INCLUDE'ABA_PARAM.INC'doubleprecisiondetJdoubleprecisionj(3,3)detJ=0.D0detJ=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)returnend·私信获取全部程序·来源：有限元先生