ABAQUS UEL二次开发（静力通用C3D20自定义单元 ）
更新于2024年4月1日 10:54浏览：2582 评论：1 收藏：1
概述：开发了适用于静力通用计算的三维二十节点（C3D20）的用户自定义单元，在挖孔悬臂梁受剪切荷载算例中，位移计算结果与ABAQUS自带单元保持一致。对比刚度矩阵，与abaqus保持一致。
（一）模型模型信息
如下图，悬臂梁尺寸10X10X100，设置四个孔洞（孔洞随意画的，具体参数不晓得，详见附件），弹性模量1e6，密度2000，泊松比0.25，荷载为1e10。
给模型施加静荷载，设置计算时长为1，固定增量步长为0.1，总增量步数10，注意静力计算中计算时长无意义，仅为了迭代求解。
（二）位移云图对比
第1增量步加载向位移云图对比如下图：
第2增量步加载向位移云图对比如下图：
第4增量步加载向位移云图对比如下图：
第6增量步加载向位移云图对比如下图：
第8增量步加载向位移云图对比如下图：
第10增量步加载向位移云图对比如下图：
第1增量步梁轴向位移云图对比如下图：
第2增量步梁轴向位移云图对比如下图：
第4增量步梁轴向位移云图对比如下图：
第6增量步梁轴向位移云图对比如下图：
第8增量步梁轴向位移云图对比如下图：
第10增量步梁轴向位移云图对比如下图：
（三）刚度矩阵对比
单元1刚度矩阵部分数据对比如下图（左边为UEL计算结果，右边是C3D20计算结果）
单元2刚度矩阵部分数据对比如下图（左边为UEL计算结果，右边是C3D20计算结果）
单元12刚度矩阵部分数据对比如下图（左边为UEL计算结果，右边是C3D20计算结果）
单元123刚度矩阵部分数据对比如下图（左边为UEL计算结果，右边是C3D20计算结果）
（四）理论与程序设计
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
-------------------------
-------------------------
（五）附件
---------------------
---------------------
ABAQUS
ESTIFFNESS.MTX
JOB-1
MO.INP
RUN.BAT
---------------------
C3D20UEL
C3D20UEL-STD.OBJ
ESTIFFNESS.MTX
JOB-1.INP
MO.INP
RUN.BAT
USER-DEFINED.INP
---------------------
---------------------
技术邻APP
工程师必备
工程师必备
- 项目客服
- 培训客服
- 平台客服
TOP
2
1
1
