零概述 开发了适用于静力通用、频率分析和动力隐式（固定增量步长和自适应增量步长均可）的三维八节点线性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 ·私信获取全部程序· 来源：有限元先生
