一种新思路用于实现ABAQUS用户自定义单元

4小时前浏览23

零

概述

一种新思路用于实现ABAQUS自定义单元，不需要在INP文件中采用“USER ELEMENT”等关键字定义单元，更重要的是，也不需要在UEL接口中进行繁琐的FORTRAN编程，只需要在INP文件中直接导入相应的矩阵即可。目前这种自定义单元适用的计算有：static, frequency extraction, modal dynamic, mode-based steady-state dynamics, complex eigenvalue extraction, and subspace-based steady-state dynamics。更多的功能还在探索之中。

这种自定义单元方式为多种CAE软件协同二次开发提供了一种可能，众所周知，ABAQUS的非线性方程组求解能力是行业翘楚，这时候如果能得到描述方程组的关键矩阵，可直接导入。

如固体力学的运动方程中的KK、MM和CC，直接导入ABAQUS便可以求解，相较于UEL二次开发，该方法节省了向ABAQUS主程序输出关键矩阵的编程工作，至于 KK、MM和CC，可以自己采用高级语言编程生成，也可以从其他软件中导出。

再比如，COMSOL以其多场耦合计算功能著称，采用该方法就可以将COMSOl的多场耦合功能与ABAQUS结合，本质是数学中的方程组求解。

这种技术可用于解决用户自定义单元在ABAUS中无法可视化的问题 。用户自定义单元的可视化，要么是将数据导出到第三方软件进行处理，要么是采用UMAT套一层单元进行可视化，这两种方法都需要大量的编程工作，采用该技术可避免上述工作。

下面以一个悬臂端受切向力作用为例，讲解实现方法：

壹

模型信息说明

模型尺寸为10x10x10，弹性模量1e10，密度2400，泊松比0.24，一端完全固定，另一端受切向力作用，切向力以节点集合的形式添加，力幅为100，边界条件和荷载示意图为：

设置两种工况：

1、ABAQUS的C3D8单元计算。

2、采用新型自定义单元，新型自定义单元的刚度矩阵和质量矩阵采用MATLAB自编程序，然后输出ABAQUS可以识别的外部文件形式，并手动修改INP文件。

贰

INP文件的具体修改说明

1、删除INP中的PART ASSEMBLY等关键词，不采用PART建模方式

2、添加如下关键字

3、将属性修改为极小值

*MATRIX INPUT,NAME=kk,INPUT=kk.txt,TYPE=UNSYMMETRIC

*MATRIX INPUT,NAME=mm,INPUT=mm.txt,TYPE=UNSYMMETRIC

*MATRIX ASSEMBLE,STIFFNESS=kk*MATRIX ASSEMBLE,mass=mm

叁

计算结果

加载向位移云图对比如下

第一主应变云图对比如下

加载向支反力云图对比如下：

来源：有限元先生

C3D20 UEL用户自定义单元开发（1）

零概述开发了适用于静力通用计算的三维二十节点（C3D20）的用户自定义单元，在挖孔悬臂梁受剪切荷载算例中，位移计算结果与ABAQUS自带单元保持一致。对比刚度矩阵，与abaqus保持一致。帖子分为两部分->理论部分-C3D20 UEL用户自定义单元开发（1）结果部分-C3D20 UEL用户自定义单元开发（2）壹二理论与程序设计 C3D20单元节点排布示意图形函数为求解刚度矩阵主程序 subroutine KKmartix(KK,coords,DMatx,mcrd,nnode,jelem) INCLUDE 'ABA_PARAM.INC' double precision coords(mcrd,nnode),dmatx(6,6) double precision integerpoint(1,3),wgt(1,3) double precision KK(3*nnode,3*nnode) double precision B(6,3*NNODE),J(3,3) double precision detJ,kesai,eta,zeta,Coff double precision npxy(3,NNODE) !积分点，二阶单元，设置三个积分点 integerpoint(1,1)= 0.D0 integerpoint(1,2)= 0.7745966692D0 integerpoint(1,3)=-0.7745966692D0!积分权重 wgt(1,1)=0.8888888889D0 wgt(1,2)=0.5555555556D0 wgt(1,3)=0.5555555556D0 B=0.d0 J=0.d0 detJ=0.d0 kesai=0.d0 eta=0.d0 zeta=0.d0 Coff=0.d0 npxy=0.d0 do i=1,3 do ii=1,3 do iii=1,3 kesai=integerpoint(1,i) eta =integerpoint(1,ii) zeta =integerpoint(1,iii) call Shapefunction(npxy,detj,kesai,eta,zeta, 1 coords,mcrd,nnode,J) call CalBmartix(B,npxy) write(6,*)"kk-detj=",detj COFF=detJ*wgt(1,i)*wgt(1,ii)*wgt(1,iii) KK=KK+matmul(matmul(transpose(B),DMatx),B)*COFF enddo enddo enddo return end 形函数对物理坐标（x、y、z）的导数 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),NPGHR(3,nnode) double precision invJ(3,3) j=0.d0 jj=0.d0 invJ=0.d0 call CalNpLocalCoord(NPGHR,kesai,eta,zeta) J=matmul(NPGHR,transpose(coords)) Call CaldetJ(detJ,J) call CalInvjocabin(invJ,detj,j) npxy=matmul(invJ,NPGHR) return end B矩阵 subroutine CalBmartix(B,npxy) INCLUDE 'ABA_PARAM.INC' double precision npxy(3,20),B(6,60) B=0.d0 do i=1,20 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(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) enddo return end 形函数对局部坐标（G、H、R）的导数 SUBROUTINE CalNpLocalCoord(NPGHR,KESAI,ETA,ZETA) INCLUDE 'ABA_PARAM.INC' DOUBLE PRECISION NPGHR(3,20) DOUBLE PRECISION KESAI,ETA,ZETA DOUBLE PRECISION PG,PH,PR DO i=1,20 CALL NPG(PG,KESAI,ETA,ZETA,I) CALL NPH(PH,KESAI,ETA,ZETA,I) CALL NPR(PR,KESAI,ETA,ZETA,I) NPGHR(1,i)=PG NPGHR(2,i)=PH NPGHR(3,i)=PR ENDDO RETURN END-------------------------!形函数对局部坐标G的导数 SUBROUTINE NPG(PG,G,H,R,NINDEX) DOUBLE PRECISION G,H,R,PG INTEGER NINDEX IF (NINDEX.EQ.1)THEN PG=-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)THEN PG=-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)THEN PG=-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)THEN PG=-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)THEN PG=-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)THEN PG=-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)THEN PG=-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)THEN PG=-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)THEN PG=0.25*( -1)*(1+G)*(1-H)*(1-R) 1 +0.25*(1-G)*( 1)*(1-H)*(1-R) ELSEIF(NINDEX.EQ.10)THEN PG=0.25*(1-H)*(1+H)*( 1)*(1-R) ELSEIF(NINDEX.EQ.11)THEN PG=0.25*( -1)*(1+G)*(1+H)*(1-R) 1 +0.25*(1-G)*( 1)*(1+H)*(1-R) ELSEIF(NINDEX.EQ.12)THEN PG=0.25*(1-H)*(1+H)*(-1)*(1-R) ELSEIF(NINDEX.EQ.13)THEN PG=0.25*( -1)*(1+G)*(1-H)*(1+R) 1 +0.25*(1-G)*( 1)*(1-H)*(1+R) ELSEIF(NINDEX.EQ.14)THEN PG=0.25*(1-H)*(1+H)*( 1)*(1+R) ELSEIF(NINDEX.EQ.15)THEN PG=0.25*( -1)*(1+G)*(1+H)*(1+R) 1 +0.25*(1-G)*( 1)*(1+H)*(1+R) ELSEIF(NINDEX.EQ.16)THEN PG=0.25*(1-H)*(1+H)*(-1)*(1+R) ELSEIF(NINDEX.EQ.17)THEN PG=0.25*(1-R)*(1+R)*(-1)*(1-H) ELSEIF(NINDEX.EQ.18)THEN PG=0.25*(1-R)*(1+R)*( 1)*(1-H) ELSEIF(NINDEX.EQ.19)THEN PG=0.25*(1-R)*(1+R)*( 1)*(1+H) ELSEIF(NINDEX.EQ.20)THEN PG=0.25*(1-R)*(1+R)*(-1)*(1+H) ENDIF RETURN END!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!形函数对局部坐标H的导数 SUBROUTINE NPH(PH,G,H,R,NINDEX) DOUBLE PRECISION G,H,R,PH INTEGER NINDEX IF (NINDEX.EQ.1)THEN PH=-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)THEN PH=-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)THEN PH=-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)THEN PH=-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)THEN PH=-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)THEN PH=-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)THEN PH=-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)THEN PH=-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)THEN PH=0.25*(1-G)*(1+G)*(-1)*(1-R) ELSEIF(NINDEX.EQ.10)THEN PH=0.25*( -1)*(1+H)*(1+G)*(1-R) 1 +0.25*(1-H)*( 1)*(1+G)*(1-R) ELSEIF(NINDEX.EQ.11)THEN PH=0.25*(1-G)*(1+G)*( 1)*(1-R) ELSEIF(NINDEX.EQ.12)THEN PH=0.25*( -1)*(1+H)*(1-G)*(1-R) 1 +0.25*(1-H)*( 1)*(1-G)*(1-R) ELSEIF(NINDEX.EQ.13)THEN PH=0.25*(1-G)*(1+G)*(-1)*(1+R) ELSEIF(NINDEX.EQ.14)THEN PH=0.25*( -1)*(1+H)*(1+G)*(1+R) 1 +0.25*(1-H)*( 1)*(1+G)*(1+R) ELSEIF(NINDEX.EQ.15)THEN PH=0.25*(1-G)*(1+G)*( 1)*(1+R) ELSEIF(NINDEX.EQ.16)THEN PH=0.25*( -1)*(1+H)*(1-G)*(1+R) 1 +0.25*(1-H)*( 1)*(1-G)*(1+R) ELSEIF(NINDEX.EQ.17)THEN PH=0.25*(1-R)*(1+R)*(1-G)*(-1) ELSEIF(NINDEX.EQ.18)THEN PH=0.25*(1-R)*(1+R)*(1+G)*(-1) ELSEIF(NINDEX.EQ.19)THEN PH=0.25*(1-R)*(1+R)*(1+G)*( 1) ELSEIF(NINDEX.EQ.20)THEN PH=0.25*(1-R)*(1+R)*(1-G)*( 1) ENDIF RETURN END!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!形函数对局部坐标R的导数 SUBROUTINE NPR(PR,G,H,R,NINDEX) DOUBLE PRECISION G,H,R,PR INTEGER NINDEX IF (NINDEX.EQ.1)THEN PR=-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)THEN PR=-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)THEN PR=-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)THEN PR=-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)THEN PR=-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)THEN PR=-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)THEN PR=-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)THEN PR=-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)THEN PR=0.25*(1-G)*(1+G)*(1-H)*(-1) ELSEIF(NINDEX.EQ.10)THEN PR=0.25*(1-H)*(1+H)*(1+G)*(-1) ELSEIF(NINDEX.EQ.11)THEN PR=0.25*(1-G)*(1+G)*(1+H)*(-1) ELSEIF(NINDEX.EQ.12)THEN PR=0.25*(1-H)*(1+H)*(1-G)*(-1) ELSEIF(NINDEX.EQ.13)THEN PR=0.25*(1-G)*(1+G)*(1-H)*( 1) ELSEIF(NINDEX.EQ.14)THEN PR=0.25*(1-H)*(1+H)*(1+G)*( 1) ELSEIF(NINDEX.EQ.15)THEN PR=0.25*(1-G)*(1+G)*(1+H)*( 1) ELSEIF(NINDEX.EQ.16)THEN PR=0.25*(1-H)*(1+H)*(1-G)*( 1) ELSEIF(NINDEX.EQ.17)THEN PR=0.25*( -1)*(1+R)*(1-G)*(1-H) 1 +0.25*(1-R)*( 1)*(1-G)*(1-H) ELSEIF(NINDEX.EQ.18)THEN PR=0.25*( -1)*(1+R)*(1+G)*(1-H) 1 +0.25*(1-R)*( 1)*(1+G)*(1-H) ELSEIF(NINDEX.EQ.19)THEN PR=0.25*( -1)*(1+R)*(1+G)*(1+H) 1 +0.25*(1-R)*( 1)*(1+G)*(1+H) ELSEIF(NINDEX.EQ.20)THEN PR=0.25*( -1)*(1+R)*(1-G)*(1+H) 1 +0.25*(1-R)*( 1)*(1-G)*(1+H) ENDIF RETURN END 来源：有限元先生