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

易公子 2022年8月12日浏览：2252 评论：12 收藏：7

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

有限元编程暂且告一段落，接下来会更新一些有关UEL自定义单元相关的内容，前几期带着大家做了一维杆单元、二维弹簧单元的UEL程序，本期推文继续更新有关常应变三角形单元——CST单元，理论基础可参照栏目《有限元基础教程》。源代码及程序文件可在公众号：易木木响叮当，后台回复：CST_UEL,即可自动获取，让木木带着你“手撕源码”吧！！！

手撕源代码

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

定义弹性矩阵 $D$ ,预先置零再赋值，具体的数值就不讲了吧（弹性力学基本常识）；

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

定义了单元形函数的偏导数，等参变换思想将任意三角形规则化，单元形函数 $N_i$ 为

$\begin{cases} N_1=\xi\\ N_2=\eta\\ N_3=1-\xi -\eta\\ \end{cases}$

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，避免了使用循环从而较少代码量：

$J=\left[ \begin{matrix} \frac{\partial x}{\partial \xi}& \frac{\partial y}{\partial \xi}\\ \frac{\partial x}{\partial \eta}& \frac{\partial y}{\partial \eta}\\ \end{matrix} \right] =\left[ \begin{matrix} \sum_{i=1}^3{\frac{\partial N_i}{\partial \xi}x_i}& \sum_{i=1}^3{\frac{\partial N_i}{\partial \xi}y_i}\\ \sum_{i=1}^3{\frac{\partial N_i}{\partial \eta}x_i}& \sum_{i=1}^3{\frac{\partial N_i}{\partial \eta}y_i}\\ \end{matrix} \right]$ $J^{-1}=\frac{1}{\det J}\left[ \begin{matrix} J_{22}& -J_{12}\\ -J_{21}& J_{11}\\ \end{matrix} \right]$

$\left[ \begin{array}{c} \frac{\partial N}{\partial \xi}\\ \frac{\partial N}{\partial \eta}\\ \end{array} \right] =J\left[ \begin{array}{c} \frac{\partial N}{\partial x}\\ \frac{\partial N}{\partial y}\\ \end{array} \right]$

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

定义了应变矩阵 $B$ 矩阵，平面三节点，每个节点上具有两个自由度，故 $B$ 矩阵的维度是 $3\times 6$ 。以上代码用的是作者对 $B$ 矩阵公式的解读，大家可以将 $B$ 矩阵拆开再乘积，也可以使用面积坐标进行编程，主要是掌握概念，代码风格因人而异。

         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子程序最为核心的内容，同时也是有限元程序最基本的元素。

$k^e=B^TDB\,\,t\,\,A$

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 - F = 0$

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

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

Abaqus验证

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

UEL单元开发（3）—— CST单元的图1

UEL单元开发（3）—— CST单元的图2

谢谢你看完木木同学的分享，今日份阅读花费的流量+1M哈哈哈哈哈哈 UEL单元开发（3）—— CST单元的图3 。

-End-

公众号：易木木响叮当

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

评论 (10条)

默认最新

吗咪吗咪轰

解释了不少接口的含义，还有编程思路，真的有很大帮助！！

2022年8月31日

评论点赞 1

zxkcup

讲解详细👍

2022年8月18日

评论 1 点赞 1

易公子回复大连海事大学

有帮助就好

2022年8月19日

评论点赞

查看更多评论 >

点赞 13 评论 12 收藏 7

易公子

关注

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

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

手撕源代码

Abaqus验证

推荐阅读