Abaqus子程序之问？代码中变量名写错了有影响吗

4月前浏览1113

Abaqus是知名的通用非线性有限元软件。其擅长处理非线性（材料、几何、接触），是非线性有限元软件的集大成者。在隐式非线性有限元软件中，我个人认为即使不是世界第一，也基本上是世界前三，况且如果还考虑知名度和用户量的因素来看，隐式非线性有限元计算abaqus大概率就是第一。

Abaqus的用户子程序是其重要的功能，其允许用户通过编写子程序实现对Abaqus内核的二次开发。一般采用Fortran编写，但也有比较少见的情况下使用C++.常用的子程序包括以下几种：

. 常用子程序分类

材料类（最常用）

UMAT：Standard 隐式分析，自定义材料本构（应力 - 应变关系），需提供雅可比矩阵（切线刚度）。

VUMAT：Explicit 显式动力学，自定义材料本构，无需雅可比，适合高应变率 / 大变形。

USDFLD：定义场变量（如温度、损伤、浓度），可依赖解状态，用于材料属性随场变量变化。

载荷 / 边界类

DLOAD：自定义非均匀体载荷 / 面载荷（如随空间 / 时间变化的压力、离心力）。

DFLUX：自定义热流 / 热源（如焊接热源、激光加热、内热源）。

UAMP：自定义幅值函数（时间 / 频率 / 空间相关的载荷幅值）。

DISP：自定义位移边界条件（复杂时变 / 空变位移、耦合约束）。

单元 / 接触类

UEL：Standard 中自定义有限元单元（如特殊梁 / 壳 / 实体、用户自定义连接），需定义刚度矩阵与残差。

VUEL：Explicit 中自定义显式单元。

UINTER：自定义接触界面行为（摩擦、粘附、损伤、接触刚度）。

输出与控制类

UVARM：定义用户自定义输出变量（积分点 / 节点结果，如损伤、等效塑性应变）。

URDFIL：读取 / 写入结果文件（.fil），用于后处理二次开发或数据交换。

PNEWDT：控制时间步长（根据求解状态调整步长，提高稳定性 / 效率）。

在编写Abaqus的Fortran子程序时，变量名写错是经常出现的情况。这里一方面是Fortran子程序通常采用了Fortran的隐式声明，即变量不需要声明，直接就可以使用。另一方面则是编写时不可避免地手误导致了变量名写错。

对于使用了隐式声明的Fortran程序，其限定以I,J,K,L,M,N开头为整型数据，其余字母开头为实数型数据。隐式声明的好处是变量随用随写，十分方便，但坏处其实是很多的，编译器无法自动检测出代码前后不一致的变量导致实际上很难发现。比如，把变量stress错写成strees等，编译器并不会直接告诉我们变量名写错了，经常需要重复调试不断检测才能发现这个错误。

本文主要简要介绍两种形式的变量名写错，并分别介绍其是否会对实际计算产生影响。

（1）第一种情况：变量名前后不一致

这种情况通常会有影响，例如导致无法收敛，计算结果错误等。但也可能只是写错了某个中间量最后没用到所以没影响。这里以一个线弹性Umat为例。

以下是一个写错了变量名（变量名前后不一致）的线弹性Umat子程序：

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,     1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,     2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,     3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,     4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)C      INCLUDE 'ABA_PARAM.INC'C      CHARACTER*80 MATERL      DIMENSION STRESS(NTENS),STATEV(NSTATV),     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),     3 PROPS(NPROPS),COORDS(3),DROT(3,3),     4 DFGRD0(3,3),DFGRD1(3,3)C      DIMENSION EELAS(6),EPLAS(6),FLOW(6)      PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)      DATA NEWTON,TOLER/10,1.D-6/CC -----------------------------------------------------------C     UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC PLASTICITYC     J2 FLOW THEORYC     CAN NOT BE USED FOR PLANE STRESSC -----------------------------------------------------------C     PROPS(1) - EC     PROPS(2) - NUC     PROPS(3) - SYIELDC     CALLS AHARD FOR CURVE OF SYIELD VS. PEEQC -----------------------------------------------------------C      IF (NDI.NE.3) THEN         WRITE(6,1) 1       FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',     1          'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')      ENDIFCC     ELASTIC PROPERTIESC      EMOD=PROPS(1)      ENU=PROPS(2)      IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499      EBULK3=EMOD/(ONE-TWO*ENU)      EG2=EMOD/(ONE+ENU)      EG=EG2/TWO      EG3=THREE*EG      ELAM=(EBULK3-EG2)/THREECC     ELASTIC STIFFNESSC      DO 20 K1=1,NTENS        DO 10 K2=1,NTENS           DDSDDE(K2,K1)=0.0 10     CONTINUE 20   CONTINUEC      DO 40 K1=1,NDI        DO 30 K2=1,NDI           DDSDDE(K2,K1)=ELAM 30     CONTINUE        DDSDDE(K1,K1)=EG2+ELAM 40   CONTINUE      DO 50 K1=NDI+1,NTENS        DDSDDE(K1,K1)=EG 50   CONTINUECC    CALCULATE STRESS FROM ELASTIC STRAINSC      DO 70 K1=1,NTENS        DO 60 K2=1,NTENS           STRESS(K2)=STRESS(K2)+ddssde(K2,K1)*DSTRAN(K1) 60     CONTINUE 70   CONTINUECC    RECOVER ELASTIC AND PLASTIC STRAINS      C      RETURN      ENDCC

这里，第70行我们把本应该是DDSDDE的变量写成了ddssde，提交计算，软件报错如下：

Error in job 1elem-static-LOAD-UMAT: Problem during linking - Abaqus/Standard User Subroutines. This error may be due to a mis match in the Abaqus user subroutine arguments. These arguments sometimes change from release to release, so user subroutines used with a previous release of Abaqus may need to be adjusted.

这里只给出了一个变量不匹配的错误提示，但是并没有给出具体哪个量不匹配。

如果是另一种方式的不匹配：

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,     1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,     2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,     3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,     4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)C      INCLUDE 'ABA_PARAM.INC'C      CHARACTER*80 MATERL      DIMENSION STRESS(NTENS),STATEV(NSTATV),     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),     3 PROPS(NPROPS),COORDS(3),DROT(3,3),     4 DFGRD0(3,3),DFGRD1(3,3)C      DIMENSION EELAS(6),EPLAS(6),FLOW(6)      PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)      DATA NEWTON,TOLER/10,1.D-6/CC -----------------------------------------------------------C     UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC PLASTICITYC     J2 FLOW THEORYC     CAN NOT BE USED FOR PLANE STRESSC -----------------------------------------------------------C     PROPS(1) - EC     PROPS(2) - NUC     PROPS(3) - SYIELDC     CALLS AHARD FOR CURVE OF SYIELD VS. PEEQC -----------------------------------------------------------C      IF (NDI.NE.3) THEN         WRITE(6,1) 1       FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',     1          'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')      ENDIFCC     ELASTIC PROPERTIESC      EMOD=PROPS(1)      ENU=PROPS(2)      IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499      EBULK3=EMOD/(ONE-TWO*ENU)      EG2=EMOD/(ONE+ENU)      EG=EG2/TWO      EG3=THREE*EG      ELAM=(EBULK3-eg22)/THREECC     ELASTIC STIFFNESSC      DO 20 K1=1,NTENS        DO 10 K2=1,NTENS           DDSDDE(K2,K1)=0.0 10     CONTINUE 20   CONTINUEC      DO 40 K1=1,NDI        DO 30 K2=1,NDI           DDSDDE(K2,K1)=ELAM 30     CONTINUE        DDSDDE(K1,K1)=EG2+ELAM 40   CONTINUE      DO 50 K1=NDI+1,NTENS        DDSDDE(K1,K1)=EG 50   CONTINUECC    CALCULATE STRESS FROM ELASTIC STRAINSC      DO 70 K1=1,NTENS        DO 60 K2=1,NTENS           STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) 60     CONTINUE 70   CONTINUECC    RECOVER ELASTIC AND PLASTIC STRAINS      C      RETURN      ENDCC

这里，第46行的EG2被错写成了eg22。提交计算，中间没有任何报错提示，直接计算完了。由于是线弹性，我们可以用自带的线弹性材料和umat比较计算结果：

二者虽然趋势一致，但是计算结果数值明显不一致。这种情况相对于前一种来说，有时候可能更难找出具体原因，软件并没有告诉我们是什么原因导致的计算结果不一致。有可能是算法，也可能是代码编写，这对于排查来说就更困难一点了。

（2）第二种情况：表头虚参变量名写错，但整个代码始终统一是错的变量名

虚参指的是子程序表头括号里的变量，如果在表头里这里面的变量写错了某一个变量名，但整个代码都是用的这个错误的变量名，这种情况下，尽管变量名是错的，但其实不影响计算结果。

      SUBROUTINE UMAT(stresss,STATEV,DDSDDE,SSE,SPD,SCD,     1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,     2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,     3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,     4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)C      INCLUDE 'ABA_PARAM.INC'C      CHARACTER*80 MATERL      DIMENSION stresss(NTENS),STATEV(NSTATV),     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),     3 PROPS(NPROPS),COORDS(3),DROT(3,3),     4 DFGRD0(3,3),DFGRD1(3,3)C      DIMENSION EELAS(6),EPLAS(6),FLOW(6)      PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)      DATA NEWTON,TOLER/10,1.D-6/CC -----------------------------------------------------------C     UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC PLASTICITYC     J2 FLOW THEORYC     CAN NOT BE USED FOR PLANE stresssC -----------------------------------------------------------C     PROPS(1) - EC     PROPS(2) - NUC     PROPS(3) - SYIELDC     CALLS AHARD FOR CURVE OF SYIELD VS. PEEQC -----------------------------------------------------------C      IF (NDI.NE.3) THEN         WRITE(6,1) 1       FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',     1          'ELEMENTS WITH THREE DIRECT stresss COMPONENTS')      ENDIFCC     ELASTIC PROPERTIESC      EMOD=PROPS(1)      ENU=PROPS(2)      IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499      EBULK3=EMOD/(ONE-TWO*ENU)      EG2=EMOD/(ONE+ENU)      EG=EG2/TWO      EG3=THREE*EG      ELAM=(EBULK3-EG2)/THREECC     ELASTIC STIFFNESSC      DO 20 K1=1,NTENS        DO 10 K2=1,NTENS           DDSDDE(K2,K1)=0.0 10     CONTINUE 20   CONTINUEC      DO 40 K1=1,NDI        DO 30 K2=1,NDI           DDSDDE(K2,K1)=ELAM 30     CONTINUE        DDSDDE(K1,K1)=EG2+ELAM 40   CONTINUE      DO 50 K1=NDI+1,NTENS        DDSDDE(K1,K1)=EG 50   CONTINUECC    CALCULATE stresss FROM ELASTIC STRAINSC      DO 70 K1=1,NTENS        DO 60 K2=1,NTENS           stresss(K2)=stresss(K2)+DDSDDE(K2,K1)*DSTRAN(K1) 60     CONTINUE 70   CONTINUECC    RECOVER ELASTIC AND PLASTIC STRAINSCC      RETURN      ENDCC