UEL单元开发（3）—— CST单元

有限元编程暂且告一段落，接下来会更新一些有关UEL自定义单元相关的内容，前几期带着大家做了一维杆单元、二维弹簧单元的UEL程序，本期推文继续更新有关常应变三角形单元——CST单元，理论基础可参照栏目《有限元基础教程》。

手撕源代码

IMPLICITREAL*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)
ENDDO
ENDDO

方法三：下策

DO K1=1,8
DO K2=1,3
        RHS(K1,1)=RHS(K1,1)-B(K2,K1)*SSTRESS(K2)*THICK*BH
ENDDO
ENDDO

注：在每一个分析步，Abaqus求解平衡方程 

式中，P为外力矢量，被Abaqus存放在数组RHS中 在迭代法求解中，令     式中：R为残余力矢量，问题转化为求解    .

个人感觉：刚入门阶段先不用管    ，认为它是0即可。

Abaqus验证

在 Abaqus 建立任意形状的单个单元模型，采用内置的CPS3单元，2、3节点固定，1节点施加 U1 方向位移5。

谢谢你看完木木同学的分享，今日份阅读花费的流量+1M哈哈哈哈哈哈。

-End-

易木木响叮当

想陪你一起度过短暂且漫长的科研生活