接上一章内容，继续かくしゅう：

下一个：

代码片段 2：

C*****Make B the Matrix           DO I=1,2             DO K=1,NNODE               BL(I,K)=0.0D0             ENDDO           ENDDO
           DO I=1,2             DO J=1,NNODE               DO K=1,2                  BL(I,J)=BL(I,J)+BR(I,K)*BN(K,J)               ENDDO             ENDDO           ENDDO
      write(6,*) 'BR Jacobian Inverse Matrix'      DO K1=1,2       write(6,*)(BR(K1,K2),K2=1,2)      ENDDO         write(6,*) 'BN'      DO K1=1,2       write(6,*)(BN(K1,K2),K2=1,4)      ENDDO         write(6,*) 'BL'      DO K1=1,2       write(6,*)(BL(K1,K2),K2=1,4)      ENDDO  
C*****Make B the Matrix                    DO I=1,3             DO K=1,8               B(I,K)=0.0D0             ENDDO           ENDDO
           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)           ENDDO  
      write(6,*) 'B'      DO K1=1,2       write(6,*)(B(K1,K2),K2=1,8)??????ENDDO??

--------------1。B应变矩阵的组成意义------------

B矩阵是一个3x8的应变矩阵。可以看出，B矩阵由以下四部分B1、B2、B3、B4组成。这是因为任何Bi都对应Ni,s和Ni,t（即局部坐标偏导的形函数）刚度矩阵计算公式，i代表节点号，在4节点平面上为1,2,3,4压力单位。特别地，下式B中的s和t分别是高斯积分点的局部坐标，这是因为组成每个Bi的Ni,s和Ni,t都是s,t的函数（即局部坐标高斯积分点函数），由于有4个高斯积分点，对应代码段1中Kintk=1,4的外层大循环。

--------------2.B应变矩阵的重要性--------------

B矩阵是3x8的应变矩阵，D是3x3的本构弹性矩阵。如下图所示，可以得到单元刚度k矩阵。因此，B应变矩阵是计算刚度矩阵不可或缺的重要矩阵。由于直接对下式进行积分得到刚度矩阵非常困难刚度矩阵计算公式，因此一般采用四个高斯积分点的形式来等效代替刚度矩阵的计算。

使用4个高斯积分点的刚度矩阵的计算公式如下，取决于4个高斯积分点的局部坐标。

------3。B应变矩阵的求解过程——BL矩阵求解------

从代码段2可以看出，B应变矩阵的求解主要与代码段中的BL矩阵（2*4）、BR（2*2）矩阵、BN（2*4）矩阵有关2、BR矩阵(2*2)是雅可比矩阵BJ矩阵(2*2)的逆矩阵，然后BL矩阵(2*4)由BR(2*2)之间的相关运算组成矩阵和BN(2*4)矩阵。其中，这个BL矩阵（2*4）表示由4个部分组成的Bi，每个部分Bi只有两个行列式，也就是2*4维度的意思。如下式所示，以B1(3*2)为例，决定B1的因素只有a(N1,s)-b(N1,t)和c(N1,t)-d( N1,s) 。

在决定B1的两个因素中，a,-b,c,-d是构成雅可比矩阵的逆矩阵BR(2*2)的四个参数，如下式所示。这里比较一下雅可比矩阵和雅可比逆矩阵，以免混淆。

同时(N1,s)和(N1,t)是前面提到的BN(2*4)矩阵，即形函数对局部坐标s,t的导数。注意BN(2*4)矩阵是4个高斯积分点的函数，取值取决于4个高斯积分点。

其中，第二维4表示4个节点，即Ni,t和Ni,s中的i；第一维2表示对s和t的偏导数的两个差异，即在Ni,t和Ni,ss和t中，而s的偏导数离开t或eta(η)，t的偏导数离开s 或 xi (ξ)。这里说明一下，代码中的s局部坐标对应有限元课本上的希腊字母xi(ξ)，代码中的t局部坐标对应有限元课本上的希腊字母eta(η)。如下所示。

4. B应变矩阵的求解过程-B应变矩阵（装配）

解完BL矩阵，即得到4个部分和每个部分中Bi的两个决定元素，可以根据下面的代码和两个决定元素在B应变矩阵中的对应位置进行集成拼装，就可以得到4个部分，每个部分的Bi。

注意B矩阵是一个3*8的矩阵，3表示对应三个应变分量，即εx，εy，γxy，8表示对应决定每个Bi的两个决定元素，一共of 4 Bi，即 4X2=8 ，即8个自由度（4个节点，每个节点x,y两个自由度）。下图说明了张量剪应变 εxy 和工程剪应变 γxy 之间的关系。工程剪应变 γ 是总夹角变化，是张量剪应变 εxy 的两倍。

------------------5.附录--------------------

以下是上一节中介绍的 UEL 子例程：

代码片段 1：

C--------------------------------------------------------------------------------------------------CC  USER DEFINED ELEMENT:2D 4NODES ISOPARAMETER LINEAR-ELASTIC MATERIA C  THIS CODE IS PROGRAMMED IN TWO purposed C  (1) : TO GET FAMILAR WITH THE PROCEDURE OF ABAQUS C  (2) : TO GET TO KNOW THE THE ORDINARY ELEMENT DEFINITON AND FE C--------------------------------------------------------------------------------------------------C  C       C--------------------------------------------------------------------------------------------------CC  MOST GENERAL  2D STRUCTURAL ELEMENT HAS TWO MATERIAL PROPERTIES: E  THICK V C STATE VARIABLES ARE STRESS AND STRAIN IN EACH POI        C  NSVARS:4(INTEGRATION POINT)X(7)=28,THEN SVARS   C  NDOFEL:THE DEGREE OF FREEDOM= 4*2    C--------------------------------------------------------------------------------------------------C      SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,      1     PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,     2     KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,     3     NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,     4     PERIOD)
      INCLUDE 'ABA_PARAM.INC'C      DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),     1     SVARS(*),ENERGY(7),PROPS(3),COORDS(MCRD,NNODE),     2     U(NDOFEL),DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),     3     PARAMS(*),JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),     4     DDLMAG(MDLOAD,*),PREDEF(2,NPREDF,NNODE),LFLAGS(4),     5     JPROPS(*)CC--------------------------------------------------------------------------------------------------CC     DEFINE  PARAMETERS AND M   c       I-N开头的变量名会被fortran识别为整型，切记！C----------------------------------           PARAMETER (HALF=0.5D0,ZERO= 0.D0,ONE = 1.D0,TWO=2.D0)             DIMENSION Gauss(NNODE,2)                                  !高斯积分点坐标      DIMENSION BN(2,NNODE)                                     !形函数对xi和eta的偏导      DIMENSION BJ(2,2),BL(2,NNODE),BR(2,2)                     !雅克比矩阵和过渡矩阵         DIMENSION B(3,NDOFEL),D(3,3),S(NDOFEL,3)                  !B阵、D阵、S阵      DIMENSION SSTRAIN(3),SSTRESS(3),Psresult(3)               !应力应变矩阵      DIMENSION CenCoord(1,3),GaussCoord(4,2)                   !单元中心点坐标，高斯积分点坐标C      C--------------------------------------------------------------------------------------------------CC    INITAIL  RHS、AMATRX、Gauss     C----------------------------------           DO K1 = 1, NDOFEL                                 DO KRHS = 1, NRHS             RHS(K1,KRHS) = ZERO           ENDDO           DO K2 = 1, NDOFEL            AMATRX(K2,K1) = ZERO           ENDDO      ENDDO
        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

      DO Kintk=1,4           xi=Gauss(Kintk,1)           eta=Gauss(Kintk,2)   CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  C          Make B Matrix  in Natural Coordinate System                  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC  C                  DO I=1,2             DO J=1,NNODE               BN(I,J)=ZERO             ENDDO           ENDDOC                   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)C                 c******BJ  Jacobian Matrix           DO I=1,2             DO K=1,2                BJ(I,K)=0.0D0             ENDDO           ENDDO
           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)           ENDDO        WRITE (6,99) BJ(1,1),BJ(1,2),BJ(2,1),BJ(2,2)99    FORMAT(//,'BJ(1,1),BJ(1,2),BJ(2,1) BJ(2,2) are equal to',3X,D14.7,     1  3X,D14.7,3X, D14.7,3X, D14.7) C******BR Jacobian Inverse Matrix              DO I=1,2              DO K=1,2                   BR(I,K)=0.0D0              ENDDO           ENDDOC                      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
       WRITE (6,100) BH,BR(1,1),BR(1,2),BR(2,1),BR(2,2) 100    FORMAT(//,'BH,BR(1,1),BR(1,2),BR(2,1), BR(2,2) are equal to',D14.7,?????1???D14.7,?D14.7,?D14.7,D14.7)?

课程推荐

本文对应的课程链接请点击文末原文链接。课程还在更新中，待定！

UEL 自定义单元子例程动手案例研究（语言）

过去的推荐

三月底，同济樱花盛开