零概述有限元法编写二维四边形二阶用户自定义单元，采用abaqus提供的uel二次开发平台，包括单元的刚度矩阵K、质量矩阵M和阻尼矩阵C，uel需要向abaqus主程序输出AMARTIX、RHS等矩阵。编写的uel子程序适用于静动力计算、频域计算等等，帖子分为两个部分->理论部分CPE8/CPS8UEL用户自定义单元开发（1）结果部分CPE8/CPS8UEL用户自定义单元开发（2）壹形函数/位移插值函数在经典的三角形单元中，位移的插值是直接在子单元中进行，但是随着单元形状的复杂，直接在子单元进行计算回带来计算的困难的公式推导的困难，因此等参单元法应运而生。在等参元中，单元形状映射中的参数和位移场插值的参数数目相等，因此得名等参元。采用等参单元法计算二维八节点单元，其中，母单元为子单元，即我们绘制的网格真实形状为假设的坐标映射函数为将母单元中八个节点的坐标导入上述式子，求解出16个系数a0~a15，然后求得八节点等参元的形函数这是形状映射中x坐标的映射，对于y方向得映射，将上述所有x替换为y即可。写成求和形式写成矩阵形式为其中，对于边角点有对于边中点有上述公式，对于位移插值同样适用。贰雅各比矩阵等参元中的数值积分中涉及到自然坐标和局部坐标的转换，这种转换通过雅各比矩阵实现。根据链式求导法则有写成矩阵形式雅各比矩阵为将式代入雅各比矩阵有形函数对母单元坐标的导数，可以直接求导，因此，等参元的雅各比矩阵可以求得。则有因此，雅各比矩阵的逆矩阵也可求得。叁B矩阵有限元中，单元的刚度矩阵为D矩阵是弹性矩阵，在线弹性问题中，D矩阵为常数矩阵，根据物理方程可以推导出。B矩阵是应变-位移矩阵，用作将位移转化为应变，下面重点求解B矩阵。应变与位移有如下关系其中算子L为单元的位移u为则应变与位移的关系可表示为则B矩阵为B=LN，具体表达式为将B矩阵写成一个矩阵为其中Bi为在求B矩阵的时候，我们的目标是求得形函数对xy的倒数，这时候我们已经求得了形函数对局部坐标的倒数，因此，通过雅各比矩阵我们就可有求得形函数对xy的倒数，首先，引入以下矩阵通过雅各比矩阵，上述矩阵可以表示为上面等式右面的雅各比逆矩阵，我们已经求得，形函数对st坐标的导数也已经求得，因此，上面引入的矩阵就可以求出来，然后再按照元素对应的规律将上述矩阵以此存放到B矩阵中。在计算质量矩阵的时候，首先计算单元的额协调质量矩阵，然后将行质量集中到对角元，让后将非对角元素置为零，具体可参考帖子：质量矩阵的求解有限元先生，公众号：有限元先生C3D8BBARUEL用户自定义单元二次开发（1）单元的阻尼矩阵可按照比例阻尼求解，即刚度矩阵和质量矩阵的线性组合。求解刚度矩阵的主程序为subroutineKKmartix(KK,coords,DMatx,mcrd,nnode,jelem)INCLUDE'ABA_PARAM.INC'doubleprecisioncoords(mcrd,nnode),dmatx(3,3)doubleprecisionINTEGERPOINT(3),wgt(3)doubleprecisionKK(2*nnode,2*nnode)doubleprecisionB(3,2*NNODE),J(2,2)doubleprecisionDETJ,G,H,COFFdoubleprecisionNPXY(2,NNODE)INTEGERPOINT(1)=0.D0INTEGERPOINT(2)=0.7745966692D0INTEGERPOINT(3)=-0.7745966692D0wgt(1)=0.8888888889D0wgt(2)=0.5555555556D0wgt(3)=0.5555555556D0B=0.d0J=0.d0DETJ=0.d0G=0.d0H=0.d0COFF=0.d0NPXY=0.d0doI=1,3doII=1,3G=INTEGERPOINT(I)H=INTEGERPOINT(II)callCalBJ(B,J,DETJ,G,H,MCRD,NNODE,JELEM,COORDS)COFF=DETJ*WGT(I)*WGT(II)KK=KK+MATMUL(MATMUL(TRANSPOSE(B),DMatx),B)*COFFenddoenddoreturnend求解质量矩阵的主程序为subroutineMMmartix(MM,COORDS,PROPS,MCRD,NNODE,JELEM)INCLUDE'ABA_PARAM.INC'INTEGERMCRD,NNODE,JELEMDOUBLEPRECISIONMM(2*NNODE,2*NNODE),COORDS(MCRD,NNODE)DOUBLEPRECISIONNN(2,16)DOUBLEPRECISIONINTEGERPOINT(3),WGT(3)DOUBLEPRECISIONG,H,DETJ,DENSITY,COFFDOUBLEPRECISIONJ(2,2),PROPS(*)MM=0.D0wgt=0.D0JJ=0.D0NN=0.D0G=0.D0H=0.D0j=0.D0DETJ=0.D0Coff=0.D0DENSITY=0.D0INTEGERPOINT(1)=0.D0INTEGERPOINT(2)=0.7745966692D0INTEGERPOINT(3)=-0.7745966692D0wgt(1)=0.8888888889D0wgt(2)=0.5555555556D0wgt(3)=0.5555555556D0DENSITY=PROPS(3)DOI=1,3DOII=1,3G=INTEGERPOINT(I)H=INTEGERPOINT(II)CALLCalDETJ(J,DETJ,G,H,MCRD,NNODE,JELEM,COORDS)CALLCalNN(NN,G,H)COFF=DETJ*WGT(I)*WGT(II)*DENSITYMM=MM+(MATMUL(TRANSPOSE(NN),NN))*COFFENDDOENDDOreturnend至此，单元的刚度矩阵、质量矩阵和阻尼矩阵已经求得，uel程序的核心部分已经完成，下一个帖子会更新自编CPE8/CPS8uel程序的计算结果。来源：有限元先生
