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