木木带你手撕Abaqus UEL子程序——Q4单元

本次分享的是四边形等参单元 Q4 单元 UEL 子程序相关内容，相信关注木木的同学已经看到了木木做的一些有关有限元编程、Umat、Uel、Fortran语法、Abaqus数值实现讲解，都是基于入门小白阶段的内容，木木想的是要把这些简单的模型，用不同的方式演绎、讲解，才能对以后复杂的模型有进一步理解，而且如果刚上来就接触比较成熟的代码，会很大可能感觉“头疼”，甚至放弃。所以木木想要在目前的这个阶段将这些基本原理、原始模型带着大家过一遍，做一些“入门工作”。

平面四节点单元等参描述

讲到 Q4 单元，那就必须以有限元“等参思想”为基础去理解。关于等参单元已经在往期的有限元编程系列提及过，这里为了照顾第一次看到此文章的童鞋，木木就把之前的描述性话语重新搬出来~是非常重要的理论基础哦  

图 1 等参单元坐标描述

我们可以从上图看出，任意四边形被边长为 1 的正方形所代替，坐标系也相应地发生了改变，改变后的坐标系习惯上称为“母单元坐标系”，母单元的形函数我们是已知的（前辈们推导的），那么我们在参与计算时还需要知道原来坐标系“子单元坐标系”的形函数，这时候出现了一个重要的概念——雅可比(Jacobian)矩阵，本身是一个数学概念，用于有限元中表示两个坐标系的转换，即：

对于平面四节点单元，Jacobian 矩阵可写为：

从以上两个式子对 Jacobian 矩阵的描述，我们可以看出：Jocabian 矩阵本身是表示两个坐标系 Mapping 的关系；可以通过 Jacobian 进行形函数偏导数的转换。

以上就是等参单元的核心思想，接下来的工作就是求解单元刚度矩阵了，采用高斯--勒让德积分，本次讨论的是完全积分，也就是四节点四个积分点。具体细节在下面的 UEL 代码中讲解并实现。

Q4 单元 UEL 子程序解读

PARAMETER (HALF=0.5D0,ZERO= 0.D0,ONE = 1.D0,TWO=2.D0)
DIMENSION Gauss(NNODE,2)  
DIMENSION BN(2,NNODE)      
DIMENSION BJ(2,2),BL(2,NNODE),BR(2,2)                     
DIMENSION B(3,NDOFEL),D(3,3),S(NDOFEL,3)                
DIMENSION SSTRAIN(3),SSTRESS(3),Psresult(3)             
DIMENSION CenCoord(1,3),GaussCoord(4,2)        

定义程序所需要用到的数组与常数：

• • Gauss(NNODE,2)：高斯积分点坐标

• • BN(2,NNODE)：形函数对母单元坐标系的偏导N_r、N_s

• • BJ(2,2),BL(2,NNODE),BR(2,2)：雅克比矩阵和过渡矩阵

• • B(3,NDOFEL),D(3,3),S(NDOFEL,3)：B阵、D阵、S阵

• • SSTRAIN(3),SSTRESS(3),Psresult(3)：应力应变矩阵

• • CenCoord(1,3),GaussCoord(4,2)：单元中心点坐标，高斯积分点坐标

Gauss(1,1)=SQRT(1/3.0D0)*-1.0D0
Gauss(1,2)=SQRT(1/3.0D0)*-1.0D0
Gauss(2,1)=SQRT(1/3.0D0)
Gauss(2,2)=SQRT(1/3.0D0)*-1.0D0
Gauss(3,1)=SQRT(1/3.0D0)
Gauss(3,2)=SQRT(1/3.0D0)
Gauss(4,1)=SQRT(1/3.0D0)*-1.0D0
Gauss(4,2)=SQRT(1/3.0D0)*1.0D0

四个高斯点坐标，相应的积分点，权系数数值分析教材中肯定会涉及，针对两点高斯积分，积分点为    ，权系数为 1。

图 2 高斯积分点坐标

DO Kintk=1,4
...
End Do

对四个高斯点进行循环，

xi=Gauss(Kintk,1)
eta=Gauss(Kintk,2)
DO I=1,2
    DO J=1,NNODE
        BN(I,J)=ZERO
    END DO
END DO
BN(1,1)=0.25D0*(-1+eta)
BN(1,2)=0.25D0*(1-eta)
BN(1,3)=0.25D0*(1+eta)
BN(1,4)=0.25D0*(-1-eta)
BN(2,1)=0.25D0*(-1+xi)
BN(2,2)=0.25D0*(-1-xi)
BN(2,3)=0.25D0*(1+xi)
BN(2,4)=0.25D0*(1-xi)

将高斯点坐标赋予给（xi,eta）四节点等参单元形函数偏导数BN。

DO J=1,NNODE
    BJ(1,1)=BJ(1,1)+BN(1,J)*COORDS(1,J)
    BJ(1,2)=BJ(1,2)+BN(1,J)*COORDS(2,J)
    BJ(2,1)=BJ(2,1)+BN(2,J)*COORDS(1,J)
    BJ(2,2)=BJ(2,2)+BN(2,J)*COORDS(2,J)
END DO
DO I=1,2
    DO K=1,2
        BR(I,K)=0.0D0
    END DO
END DO
BH=BJ(2,2)*BJ(1,1)-BJ(1,2)*BJ(2,1)
BR(1,1)=BJ(2,2)/BH
BR(1,2)=-1*BJ(1,2)/BH
BR(2,1)=-1*BJ(2,1)/BH
BR(2,2)=BJ(1,1)/BH

计算雅可比矩阵BJ，是个      的矩阵，后而计算雅可比逆矩阵BR。Fortran 没有内置的逆矩阵函数，所以计算子程序在计算逆矩阵时可以编写对应的子程序，或者直接写出逆矩阵的直接形式。

call    cal_D(1,D,E,ENU)
·····
subroutine cal_D(Pss,D,E,ENU) 
INCLUDE 'ABA_PARAM.INC'
double precision E ENU D1 D2
INTEGER Pss
double precision, DIMENSION(3,3)::D
PARAMETER (ZERO= 0.D0,ONE = 1.D0,TWO=2.D0)
DO I=1,3
    DO K=1,3
        D(I,K)=ZERO
    END DO
END DO
IF (Pss.EQ.1) THEN
C  平面应力状态
    D1 = E/(ONE-EMU*EMU)
    D(1,1)=ONE*D1
    D(2,2)=ONE*D1
    D(1,2)=ENU*D1
    D(2,1)=ENU*D1
    D(3,3)=(ONE-EMU)/TWO*D1
ELSE IF (Pss.EQ.2) THEN
C   平面应力状态
    D2 =  E*(ONE-ENU)/((ONE+EMU)*(ONE-TWO*ENU))
    D(1,1)=ONE*D2
    D(2,2)=ONE*D2
    D(1,2)=ENU/(ONE-ENU)*D2
    D(2,1)=ENU/(ONE-ENU)*D2
    D(3,3)=(ONE-TWO*EMU)/(TWO*(1-ENU))*D2
END IF
end 

编写弹性矩阵，以上程序是给大家演示如何在 UEL 中调用自己编写的子程序。

DO I=1,2
    DO K=1,NNODE
        BL(I,K)=0.0D0
    END DO
END DO
DO I=1,2
    DO J=1,NNODE
        DO K=1,2
            BL(I,J)=BL(I,J)+BR(I,K)*BN(K,J)
        END DO
    END DO
END DO
DO I=1,3
    DO K=1,8
        B(I,K)=0.0D0
    END DO
END DO
DO K=1,NNODE
    B(1,K*2-1)=BL(1,K)
    B(2,K*2)=BL(2,K)
    B(3,K*2-1)=BL(2,K)
    B(3,K*2)=BL(1,K)
END DO

NNODE，自定义单元节点数。要有预先将矩阵元素清零的习惯，以上程序是计算应变矩阵B，公式比较多，不想敲了，就直接截取曾攀老师《有限元基础教程》里面的公式了~

hr = ONE ! 权因子
hs = ONE       AMATRX=AMATRX+THICK*hr*hs*BH*MATMUL(MATMUL(TRANSPOSE(B),D),B)
RHS(1:NDOFEL,1) = RHS(1:NDOFEL,1) - MATMUL(AMATRX,U(1:NDOFEL))

计算单元刚度矩阵AMATRX，残余力向量RHS，使用到了内置函数，按道理来说，子程序应该按照固定格式来书写，但是Abaqus子程序允许使用一些自由格式的语法，在以后的推文中会专门讲解，敬请期待哦~

Abaqus数值验证

图 3 Abaqus数值验证