abaqus dload子程序实现复杂移动荷载

外太空土豆儿

4小时前浏览7

零

概述

abaqus给用于提供了dload子程序用以满足复杂的荷载需求，帖子采用该子程序模拟了四轮小车在平板上移动的运动过程。

壹

几何模型信息

模型为四轮车，每个轮子与地面的接触区域为0.1x0.3（不考虑实际情况），车辆在平板居中放置，从右边界启动，向左移动

采用C3D8单元划分网格

网格横断面

车轮经过的区域网格设置加密，车辆中间区域与外部区域网格尺寸设置渐变。其中，车轮经过区域的网格

贰

程序

!  模拟四个车轮经过桥面，动荷载            SUBROUTINE DLOAD(F,KSTEP,KINC,TIME,NOEL,NPT,LAYER,KSPT,     1 COORDS,JLTYP,SNAME)C      INCLUDE 'ABA_PARAM.INC'C      DIMENSION TIME(2), COORDS (3)      CHARACTER*80 SNAME! FIXED INCRANTMENT TIME IS 0.01S            DOUBLE PRECISION INC_TIME      DOUBLE PRECISION F,DISP            INC_TIM=0.01
!  10s，增量步长为0.01s，总增量步数为1000      !  荷载表达式：F=sin(2*3.14/10)!  桥梁尺寸：桥长40，桥宽6      X=COORDS(1)      Y=COORDS(2)      Z=COORDS(3)
      DISP=TIME(2)*(40/10)            IF((Z>=2.2.AND.Z<=2.5).OR.(Z>=3.5.AND.Z<=3.8))THEN              IF((X>DISP-0.1).AND.X<(DISP))THEN          F=1        ELSEIF((X>DISP-0.1-2).AND.(X<DISP-2))THEN          F=1        ELSE          F=0        ENDIF      ELSE        F=0      ENDIF
      RETURN      END

来源：有限元先生

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 来源：有限元先生