UEL单元开发（3）—— CST单元
UEL单元开发（3）—— CST单元
有限元编程暂且告一段落，接下来会更新一些有关UEL自定义单元相关的内容，前几期带着大家做了一维杆单元、二维弹簧单元的UEL程序，本期推文继续更新有关常应变三角形单元——CST单元，理论基础可参照栏目《有限元基础教程》。
手撕源代码
IMPLICIT REAL*8(A-H,O-Z)以A-H,O-Z开头的变量是双精度。
REAL INTERGE,JMATRX,JINVER
DIMENSION INTERGE(3,3),DMATRX(3,3),BMATRX(3,8),HH(3),
1 dNdrs(2,3),JMATRX(2,2),JINVER(2,2),CO(2,3),
2 dNdxy(2,3)
PARAMETER (zero=0.0d0,one=1.0d0,two = 2.0D0,THREE=3.0d0,
1 TWICE=0.5d0)定义了程序所需要用的数组、大小、常数值。
YM = PROPS(1)
POIS = PROPS(2)
DO J = 1,3
DO K = 1,3
DMATRX(J,K) = zero
ENDDO
ENDDO
D0 = YM/(one-POIS**two)
DMATRX(1,1) = D0
DMATRX(1,2) = D0*POIS
DMATRX(2,1) = D0*POIS
DMATRX(2,2) = D0
DMATRX(3,3) = D0*(one-POIS)/two定义弹性矩阵
props(1)、props(2)表示单元属性，对应于Inp文件中的*Uel property, elset=Set-1下面的数据行，在该程序中表示的是弹性模量和泊松比。
DO J = 1,2
DO K = 1,3
dNdrs(J,K) = zero
ENDDO
ENDDO
dNdrs(1,1) = one
dNdrs(1,2) = zero
dNdrs(1,3) =-one
dNdrs(2,1) = zero
dNdrs(2,2) = one
dNdrs(2,3) =-one定义了单元形函数的偏导数，等参变换思想将任意三角形规则化，单元形函数
JMATRX = MATMUL(dNdrs,TRANSPOSE(COORDS(1:2,1:NNODE)))
DETJ = JMATRX(1,1)*JMATRX(2,2) - JMATRX(1,2)*JMATRX(2,1)
JINVER(1,1) = JMATRX(2,2)/DETJ
JINVER(2,2) = JMATRX(1,1)/DETJ
JINVER(1,2) =-JMATRX(1,2)/DETJ
JINVER(2,1) =-JMATRX(2,1)/DETJ
dNdxy = MATMUL(JINVER,dNdrs)定义雅可比矩阵、逆矩阵、子单元形函数偏导数，在这里值得注意的地方是计算雅可比矩阵时，我用到了Fortran内置函数Matmul、Transpaose，避免了使用循环从而较少代码量：
DO J = 1,3
DO K = 1,6
BMATRX(J,K) = zero
ENDDO
ENDDO
DO J = 1,3
BMATRX(1,(J-1)*2+1) = dNdxy(1,J)
BMATRX(2,(J-1)*2+2) = dNdxy(2,J)
BMATRX(3,(J-1)*2+1) = dNdxy(2,J)
BMATRX(3,(J-1)*2+2) = dNdxy(1,J)
ENDDO定义了应变矩阵
CO = COORDS
AREA = (CO(1,1)*(CO(2,2)-CO(2,3)) + CO(1,2)*(CO(2,3)-CO(2,1))+
1 CO(1,3)*(CO(2,1)-CO(2,2)))/two
AMATRX = AMATRX+AREA*one*MATMUL(MATMUL(TRANSPOSE(BMATRX),DMATRX),BMATRX)AREA表示的是三角形单元的面积（任意三角形面积公式），AMATRX表示的是单元刚度矩阵，这个也是UEL子程序最为核心的内容，同时也是有限元程序最基本的元素。
RHS残余力矢量RHS计算
方法一：使用内置函数 上策
RHS(1:NDOFEL,1) = RHS(1:NDOFEL,1) - MATMUL(AMATRX,U(1:NDOFEL))方法二：中策
DO K1=1,8
DO K2=1,8
RHS(K1,1)=RHS(K1,1)- AMATRX(K1,K2)*U(K2)
END DO
END DO 方法三：下策
DO K1=1,8
DO K2=1,3
RHS(K1,1)=RHS(K1,1)-B(K2,K1)*SSTRESS(K2)*THICK*BH
END DO
END DO注：在每一个分析步，Abaqus求解平衡方程
式中，P为外力矢量，被Abaqus存放在数组RHS中 在迭代法求解中，令
个人感觉：刚入门阶段先不用管
Abaqus验证
在 Abaqus 建立任意形状的单个单元模型，采用内置的CPS3单元，2、3节点固定，1节点施加 U1 方向位移5。
谢谢你看完木木同学的分享，今日份阅读花费的流量+1M哈哈哈哈哈哈
。
-End-
易木木响叮当
想陪你一起度过短暂且漫长的科研生活
