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

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

·私信获取全部程序·

来源：有限元先生

C3D8 BBAR UEL用户自定义单元二次开发(2)

零概述开发了适用于静力通用、频率分析和动力隐式（固定增量步长和自适应增量步长均可）的三维八节点线性UEL，即ABAQUS自带的C3D8单元，该UEL考虑了B-BAR修正，避免剪切锁死。采用编写的UEL，分别设置了静力通用分析步、频率分析和动力隐式分析步，将计算结果与ABAQUS对比，位移、速度和加速度与ABAQUS均保持一致，说明该UEL复现了一小部分C3D8单元的计算功能。帖子分为两部分理论部分-C3D8BBARUEL用户自定义单元开发（1）->结果部分-C3D8BBARUEL用户自定义单元开发（2）壹模型信息悬臂梁尺寸：10x10x100，弹模1e10，密度200，泊松比0.25。不设置单位，纯验证。每个单元尺寸为：2x2x5贰静力分析边界条件一段固定，另一端受静力作用荷载大小为：1e6，采用固定增量步长，计算总时长为10（静力计算中计算时长无意义，仅为验证设置），增量步长为0.01，总增量步数为1000。总位移云图加载向（U2、Y向）位移云图梁向（U3，z向）位移云图U1，x向位移云图悬臂端角点加载向位移-荷载历程注意这里的时间并没有物理意义，在静力通用分析步中仅仅是增量步的计算意义。叁振动分析边界条件，梁一端固定。计算前100阶频率信息。1~41阶41~56阶55~100阶第1阶振型第2阶振型第3阶振型第20阶振型第50阶振型：第100阶振型叁动力隐式分析给悬臂端节点集合施加简谐荷载，边界条件为简谐荷载曲线INP文件定义*Amplitude,name=Amp-1,definition=PERIODIC1,1.,0.,0.1.,1.荷载幅值为：1e6，对应的INP文件为*Cload,amplitude=Amp-1Part-1-1.forced,2,-1e+06悬臂角点加载向（U2）位移时程曲线悬臂角点轴向（U3）位移时程曲线第3s加载向加速度云图第7s加载向加速度云图第10s加载向加速度云图·私信获取全部程序·来源：有限元先生