零概述开发了适用于静力通用计算的三维二十节点（C3D20）的用户自定义单元，在挖孔悬臂梁受剪切荷载算例中，位移计算结果与ABAQUS自带单元保持一致。对比刚度矩阵，与abaqus保持一致。帖子分为两部分->理论部分-C3D20UEL用户自定义单元开发（1）结果部分-C3D20UEL用户自定义单元开发（2）壹二理论与程序设计C3D20单元节点排布示意图形函数为求解刚度矩阵主程序subroutineKKmartix(KK,coords,DMatx,mcrd,nnode,jelem)INCLUDE'ABA_PARAM.INC'doubleprecisioncoords(mcrd,nnode),dmatx(6,6)doubleprecisionintegerpoint(1,3),wgt(1,3)doubleprecisionKK(3*nnode,3*nnode)doubleprecisionB(6,3*NNODE),J(3,3)doubleprecisiondetJ,kesai,eta,zeta,Coffdoubleprecisionnpxy(3,NNODE)!积分点，二阶单元，设置三个积分点integerpoint(1,1)=0.D0integerpoint(1,2)=0.7745966692D0integerpoint(1,3)=-0.7745966692D0!积分权重wgt(1,1)=0.8888888889D0wgt(1,2)=0.5555555556D0wgt(1,3)=0.5555555556D0B=0.d0J=0.d0detJ=0.d0kesai=0.d0eta=0.d0zeta=0.d0Coff=0.d0npxy=0.d0doi=1,3doii=1,3doiii=1,3kesai=integerpoint(1,i)eta=integerpoint(1,ii)zeta=integerpoint(1,iii)callShapefunction(npxy,detj,kesai,eta,zeta,1coords,mcrd,nnode,J)callCalBmartix(B,npxy)write(6,*)"kk-detj=",detjCOFF=detJ*wgt(1,i)*wgt(1,ii)*wgt(1,iii)KK=KK+matmul(matmul(transpose(B),DMatx),B)*COFFenddoenddoenddoreturnend形函数对物理坐标（x、y、z）的导数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),NPGHR(3,nnode)doubleprecisioninvJ(3,3)j=0.d0jj=0.d0invJ=0.d0callCalNpLocalCoord(NPGHR,kesai,eta,zeta)J=matmul(NPGHR,transpose(coords))CallCaldetJ(detJ,J)callCalInvjocabin(invJ,detj,j)npxy=matmul(invJ,NPGHR)returnendB矩阵subroutineCalBmartix(B,npxy)INCLUDE'ABA_PARAM.INC'doubleprecisionnpxy(3,20),B(6,60)B=0.d0doi=1,20B(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)enddoreturnend形函数对局部坐标（G、H、R）的导数SUBROUTINECalNpLocalCoord(NPGHR,KESAI,ETA,ZETA)INCLUDE'ABA_PARAM.INC'DOUBLEPRECISIONNPGHR(3,20)DOUBLEPRECISIONKESAI,ETA,ZETADOUBLEPRECISIONPG,PH,PRDOi=1,20CALLNPG(PG,KESAI,ETA,ZETA,I)CALLNPH(PH,KESAI,ETA,ZETA,I)CALLNPR(PR,KESAI,ETA,ZETA,I)NPGHR(1,i)=PGNPGHR(2,i)=PHNPGHR(3,i)=PRENDDORETURNEND-------------------------!形函数对局部坐标G的导数SUBROUTINENPG(PG,G,H,R,NINDEX)DOUBLEPRECISIONG,H,R,PGINTEGERNINDEXIF(NINDEX.EQ.1)THENPG=-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)THENPG=-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)THENPG=-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)THENPG=-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)THENPG=-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)THENPG=-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)THENPG=-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)THENPG=-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)THENPG=0.25*(-1)*(1+G)*(1-H)*(1-R)1+0.25*(1-G)*(1)*(1-H)*(1-R)ELSEIF(NINDEX.EQ.10)THENPG=0.25*(1-H)*(1+H)*(1)*(1-R)ELSEIF(NINDEX.EQ.11)THENPG=0.25*(-1)*(1+G)*(1+H)*(1-R)1+0.25*(1-G)*(1)*(1+H)*(1-R)ELSEIF(NINDEX.EQ.12)THENPG=0.25*(1-H)*(1+H)*(-1)*(1-R)ELSEIF(NINDEX.EQ.13)THENPG=0.25*(-1)*(1+G)*(1-H)*(1+R)1+0.25*(1-G)*(1)*(1-H)*(1+R)ELSEIF(NINDEX.EQ.14)THENPG=0.25*(1-H)*(1+H)*(1)*(1+R)ELSEIF(NINDEX.EQ.15)THENPG=0.25*(-1)*(1+G)*(1+H)*(1+R)1+0.25*(1-G)*(1)*(1+H)*(1+R)ELSEIF(NINDEX.EQ.16)THENPG=0.25*(1-H)*(1+H)*(-1)*(1+R)ELSEIF(NINDEX.EQ.17)THENPG=0.25*(1-R)*(1+R)*(-1)*(1-H)ELSEIF(NINDEX.EQ.18)THENPG=0.25*(1-R)*(1+R)*(1)*(1-H)ELSEIF(NINDEX.EQ.19)THENPG=0.25*(1-R)*(1+R)*(1)*(1+H)ELSEIF(NINDEX.EQ.20)THENPG=0.25*(1-R)*(1+R)*(-1)*(1+H)ENDIFRETURNEND!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!形函数对局部坐标H的导数SUBROUTINENPH(PH,G,H,R,NINDEX)DOUBLEPRECISIONG,H,R,PHINTEGERNINDEXIF(NINDEX.EQ.1)THENPH=-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)THENPH=-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)THENPH=-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)THENPH=-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)THENPH=-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)THENPH=-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)THENPH=-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)THENPH=-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)THENPH=0.25*(1-G)*(1+G)*(-1)*(1-R)ELSEIF(NINDEX.EQ.10)THENPH=0.25*(-1)*(1+H)*(1+G)*(1-R)1+0.25*(1-H)*(1)*(1+G)*(1-R)ELSEIF(NINDEX.EQ.11)THENPH=0.25*(1-G)*(1+G)*(1)*(1-R)ELSEIF(NINDEX.EQ.12)THENPH=0.25*(-1)*(1+H)*(1-G)*(1-R)1+0.25*(1-H)*(1)*(1-G)*(1-R)ELSEIF(NINDEX.EQ.13)THENPH=0.25*(1-G)*(1+G)*(-1)*(1+R)ELSEIF(NINDEX.EQ.14)THENPH=0.25*(-1)*(1+H)*(1+G)*(1+R)1+0.25*(1-H)*(1)*(1+G)*(1+R)ELSEIF(NINDEX.EQ.15)THENPH=0.25*(1-G)*(1+G)*(1)*(1+R)ELSEIF(NINDEX.EQ.16)THENPH=0.25*(-1)*(1+H)*(1-G)*(1+R)1+0.25*(1-H)*(1)*(1-G)*(1+R)ELSEIF(NINDEX.EQ.17)THENPH=0.25*(1-R)*(1+R)*(1-G)*(-1)ELSEIF(NINDEX.EQ.18)THENPH=0.25*(1-R)*(1+R)*(1+G)*(-1)ELSEIF(NINDEX.EQ.19)THENPH=0.25*(1-R)*(1+R)*(1+G)*(1)ELSEIF(NINDEX.EQ.20)THENPH=0.25*(1-R)*(1+R)*(1-G)*(1)ENDIFRETURNEND!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!形函数对局部坐标R的导数SUBROUTINENPR(PR,G,H,R,NINDEX)DOUBLEPRECISIONG,H,R,PRINTEGERNINDEXIF(NINDEX.EQ.1)THENPR=-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)THENPR=-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)THENPR=-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)THENPR=-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)THENPR=-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)THENPR=-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)THENPR=-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)THENPR=-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)THENPR=0.25*(1-G)*(1+G)*(1-H)*(-1)ELSEIF(NINDEX.EQ.10)THENPR=0.25*(1-H)*(1+H)*(1+G)*(-1)ELSEIF(NINDEX.EQ.11)THENPR=0.25*(1-G)*(1+G)*(1+H)*(-1)ELSEIF(NINDEX.EQ.12)THENPR=0.25*(1-H)*(1+H)*(1-G)*(-1)ELSEIF(NINDEX.EQ.13)THENPR=0.25*(1-G)*(1+G)*(1-H)*(1)ELSEIF(NINDEX.EQ.14)THENPR=0.25*(1-H)*(1+H)*(1+G)*(1)ELSEIF(NINDEX.EQ.15)THENPR=0.25*(1-G)*(1+G)*(1+H)*(1)ELSEIF(NINDEX.EQ.16)THENPR=0.25*(1-H)*(1+H)*(1-G)*(1)ELSEIF(NINDEX.EQ.17)THENPR=0.25*(-1)*(1+R)*(1-G)*(1-H)1+0.25*(1-R)*(1)*(1-G)*(1-H)ELSEIF(NINDEX.EQ.18)THENPR=0.25*(-1)*(1+R)*(1+G)*(1-H)1+0.25*(1-R)*(1)*(1+G)*(1-H)ELSEIF(NINDEX.EQ.19)THENPR=0.25*(-1)*(1+R)*(1+G)*(1+H)1+0.25*(1-R)*(1)*(1+G)*(1+H)ELSEIF(NINDEX.EQ.20)THENPR=0.25*(-1)*(1+R)*(1-G)*(1+H)1+0.25*(1-R)*(1)*(1-G)*(1+H)ENDIFRETURNEND来源：有限元先生
