ANSYS UPFs 子程序简单介绍及实例

发表于 更新于

转自仿真科技论坛，点击进入原网址，作者师访

首先做一点说明，这里所说的子程序大致分为两类：一类是用户子程序，即可以被用户修改的、多以 User 或 u 开头的、可以在安装文件夹下找到（如下图所示）的、需要重新编译连接的 Fortran 子程序；另一类是一般的子程序，ANSYS 提供给用户方便 ansys 二次开发的，包括子程序（subroutine，无返回值，使用时需要 external 声明语句）和函数（Function，有返回值，使用时不需要 External 声明语句）。

材料本构模型开发用户子程序

  1. usermat：在所有用于本构模型开发的子程序中，usermat 是使用范围最广的子程序，单元的每个积分点在每个载荷子步的每次迭代均会调用 usermat。usermat 支持所有最新技术单元（Current-technology Elements），最新技术单元区别于传统单元。
  2. UserHyper：虽然 usermat 功能强大，但其无法用于开发超弹性材料，若想实现超弹性材料本构的开发，必须利用 UserHyper 用户子程序。
  3. usercreep 及 usercr：ANSYS 提供了两个用户子程序 usercreep 及 usercr 来定义用户蠕变方程，这两个用户子程序可用来模拟材料的时间相关粘塑性 / 蠕变行为。usercreep 利用隐式时间积分算法，而 usercr 利用的是显示时间积分算法。usercreep 用于较长时间段的蠕变分析，如岩石及金属的长期蠕变分析；usercr 用于瞬态蠕变分析。
  4. user_tbelastic：用来定义线弹性材料特征，允许将弹性材料参数定义成温度或坐标的函数，支持大多数最新技术单元。user_tbelastic 仅仅用于定义弹性阶段的材料本构，它可用来与许多非线性材料模型联合使用，例如 BISO、BKIN、MISO、KINH、NLSIO、PLASTIC、CHABOCHE、CAST IRON、CREEP、EDP（Drucker-Prager 模型）、SMA（形状记忆合金）、PRONY 及 GURSON 等。
  5. userpl：如果用户想利用传统单元定义塑性或粘塑性本构关系，则可以使用 userpl。适用于 LINK1、LINK8、PIPE20、BEAM23、BEAM24、PLANE42、SOLID45、PIPE60、SOLID62、SOLID65、PLANE82、SOLID92 及 SOLID95。
  6. usersw：用户膨胀子程序 usersw 用于编写用户自定义的膨胀律，其调用方式是设置 C72=10。
  7. UserVisLaw：定义粘性模型，仅用于 FLUID141 及 FLUID142 单元。
  8. userfic：用于定义用户摩擦模型，用于接触单元：CONTA171、CONTA172、CONTA173、CONTA174、CONTA175、CONTA176、CONTA177、CONTA178。
  9. userfc：定义层合板类材料的破坏准则。

单元开发用户子程序

开发单元有两种方法：用户子程序（1）和（2）联合使用是一种方法，通过 UserElem 用户子程序可以方便的建立新单元，ElemGetMat 为新单元处理材料本构部分，难度适中；（3）和（4）联合使用是另外一种方法，该方法直接访问 ANSYS 单元相关数据库和文件，难度较大。 1. UserElem：从 ANSYS 11.0 开始新增的。使用该子程序创建单元无需直接访问 ANSYS 数据库。 2. ElemGetMat：用于调用 ANSYS 标准结构材料库，与 UserElem 结合使用。 3. uec100-uec105：在通过直接访问单元相关数据库和文件的方式创建单元时，用于定义单元特征参数。 4. uel100-uel105：在通过直接访问单元相关数据库和文件的方式创建单元时，用于计算各种单元矩阵（刚度矩阵、比热矩阵等）、单元载荷向量（力载荷、热流载荷等）以及所有单元输出变量。此外还要完成单元输出过程，计算需要保存的各个变量并储存在结果文件中。

单元开发支持子程序

以下子程序（注意不是用户子程序）用于单元开发： 1. nminfo：给单元特征向量定义单元名称，输入变量为单元特征向量 ielc（见  ANSYS Inc

文件夹下的 echprm.inc），单元参考名称，输出变量为加入了单元参考名称的单元特征向量。 2. svgidx：获得保存的变量的索引。 3. svrget：获得单元保存的变量。 4. mreuse：决定在迭代过程中哪些单元矩阵及向量可被重复使用，哪些单元矩阵及向量需要重新计算。 5. rvrget：获得单元实常数。 6. propev：计算一组材料特性参数。 7. prope1：计算某一个材料特性参数，若要计算多个材料参数，可使用 propev 子程序。 8. pstev1：计算 1 维单元的材料特性参数。 9. plast1：更新单元塑性历史，用于具有一个应力（应变）分量的单元，包括 LINK1、LINK8、BEAM23、BEAM24 和 SOLID65（增强型）单元。 10. plast3：更新单元塑性历史，用于具有 4 或 6 应力（应变）分量的单元，包括 PLANE02、PLANE13、PIPE20、SHELL43、SHELL51、PIPE60、SOLID62、SOLID65、SHELL91、SHELL93、SHELL143、SOLID191。 11. creep1：更新单元蠕变历史，用于 1 维单元，包括 LINK1、LINK8、BEAM23、BEAM24 和 SOLID65 单元。 12. creep3：更新单元蠕变历史，用于 3 维单元，包括 PLANE02、PLANE13、PIPE20、PLANE42、SHELL43、SOLID45、SHELL51、PIPE60、SOLID62、SOLID65、PLANE82、SHELL91、SOLID92、SHELL93、SOLID95、SHELL143 及 SOLID191。 13. swell1：更新单元的膨胀历史，用于 1 维单元，包括 LINK1、LINK8、BEAM23、BEAM24。 14. swell3：更新单元的膨胀历史，用于 3 维单元，包括 PLANE02、PLANE13、PIPE20、PLANE42、SHELL43、SOLID45、SHELL51、PIPE60、SOLID62、PLANE82、SHELL91、SOLID92、SHELL93、SOLID95、SHELL143 和 SOLID191。 15. usereo：在不可累加杂项数据（NMISC）中存储数据。 16. eldwrtL：将单元数据写入文件。 17. eldwrnL：将单元不可累加杂项数据写入结果文件。 18. tmpget：子程序 temget 用于定义当前温度载荷。 19. prsget：定义当前压力载荷。 20. cnvget：定义当前对流载荷。 21. prinst：可用于计算主应力和应力强度。其输入输出输出变量只有一个 s（长度为 11 的数组），s(1)、s(2)、s(3)、s(4)、s(5) 和 s(6) 分别存放σx 、σy 、σz 、σxy 、σyz 和 σzx，经子程序计算后输出主应力σ1、σ2 和σ3 保存在 s(7)、s(8) 和 s(9) 内，应力强度保存在 s(10) 内，s(11) 内输出的是 Mises 等效应力。也可用来计算主应变（ε1 、ε2 和ε3 ）的大小。提醒大家，这是一个很重要的子程序。

修改和监视已存在单元的用户子程序

  1. userou：用于存储用户提供的单元输出，对于需要将单元数据存入 nmisc 记录的所有单元，都要调用该用户子程序。该子程序通过 USRCAL 命令激活：USRCAL,USEROU
  2. useran：单元坐标系控制，修改 shell43、shell63、shell91、shell93、shell99，以及 solid46、solid64、solid191 单元材料性能参数和应力的方向，可通过单元关键选项调用。
  3. userrc：对 COMBIN7 和 COMBIN37 单元的参数进行操作。当设置单元的 keyopt(9) = 1 时调用该子程序，且 C1 或 C3 必须为非零值。
  4. UElMatx：访问单元矩阵和载荷向量，如频域分析的载荷向量等。主要目的是监视（或修改）单元矩阵和载荷向量。
  5. UTHICK：定义 SHELL181、SHELL208、SHELL209 以及 SHELL281 单元的初始厚度。
  6. Us_Surf_Str：获取 PLANE2、PLANE42、PLANE82 以及 SOLID45、SOLID92 和 SOLOD95 单元的单元面应力。
  7. uflex：计算单元 PIPE288 和 PIPE289 的柔性系数。当输入单元的轴向柔性系数为 - 10 时调用该用户子程序。
  8. usflex：计算单元 PIPE16、PIPE17、PIPE18 和 PIPE60 的柔性系数。当输入单元的轴向柔性系数为负数时调用该用户子程序。

载荷用户子程序

以下用户子程序用于定义载荷： 1. usrefl：可将一些标量表示的载荷（如温度、流量、热生成量、含水量等）修改成用户定义的值，该子程序在每个使用标量场载荷的载荷步及载荷子步内、每次平衡迭代过程中均会被调用，在进入该子程序之前，已获得 ANSYS 输入的单元或节点值信息。 2. userpr：用于改变单元压力，该子程序在每个使用单元压力的载荷步及载荷子步内、每次平衡迭代过程中均会被调用，每个单元调用一次。在进入该子程序之前，已获得 ANSYS 输入的单元或节点值信息。 3. usercv：用于改变单元面对流信息，该子程序在每个载荷步及载荷子步内的每次平衡迭代过程中均会被调用，每个单元调用一次。在进入该子程序之前，已获得标准 ANSYS 输入的对流面信息，因此利用这些信息可以随意按照用户意愿修改对流信息。 4. userfx：用于改变单元面热流信息，该子程序在每个载荷步及载荷子步内的每次平衡迭代过程中均会被调用，每个单元调用一次。仅在热流载荷形成阶段调用，热流计算阶段不调用。

计算干预用户子程序

用户子程序 对应的情况
UAnBeg ANSYS 启动时执行
USolBeg 求解开始前执行
ULdBeg 载荷步开始前执行
USsBeg 载荷子步开始前执行
UItBeg 迭代步开始前执行
UItFin 迭代步结束后执行
USsFin 载荷子步结束后执行
ULdFin 载荷步结束后执行
USolFin 求解结束后执行
UAnFin ANSYS 结束前执行

访问 ANSYS 数据库的子程序

这部分内容包括选择或获得节点及单元的子程序、节点信息相关子程序、单元特征相关子程序、耦合及约束相关子程序、节点载荷子程序、单元载荷子程序、结果信息子程序，由于内容较多，参见附件：下载

方便用户开发的子程序

内容包括通用子程序、向量操作子程序和矩阵操作子程序，见附件：下载

简单例子 线弹性材料模型 usermat 开发

材料的本构行为体现为应力 - 应变关系，usermat 子程序的主要任务是定义材料的应力 - 应变关系，具体分为两个方面。首先由给定的应变增量计算得到应力增量 ，从而得到新的应力 ，称做应力更新过程。其次还要求 usermat 给出雅可比矩阵，学名是一致切线算子矩阵。usermat 用户子程序仅用于最新技术单元，不支持传统单元。

usermat fortran 用户子程序的结构

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
*deck,usermat      USERDISTRIB  parallel                                gal
subroutine usermat(
     &                   matId, elemId,kDomIntPt, kLayer, kSectPt,
     &                   ldstep,isubst,keycut,
     &                   nDirect,nShear,ncomp,nStatev,nProp,
     &                   Time,dTime,Temp,dTemp,
     &                   stress,statev,dsdePl,sedEl,sedPl,epseq,
     &                   Strain,dStrain, epsPl, prop, coords,
     &                   var0, defGrad_t, defGrad,
     &                   tsstif, epsZZ,
     &                   var1, var2, var3, var4, var5,
     &                   var6, var7, var8)
c*************************************************************************
#include "impcom.inc"
c
      INTEGER         
     &                 matId, elemId,
     &                 kDomIntPt, kLayer, kSectPt,
     &                 ldstep,isubst,keycut,
     &                 nDirect,nShear,ncomp,nStatev,nProp
      DOUBLE PRECISION
     &                 Time,    dTime,   Temp,    dTemp,
     &                 sedEl,   sedPl,   epseq,   epsZZ
      DOUBLE PRECISION
     &                 stress  (ncomp), statev  (nStatev),
     &                 dsdePl  (ncomp,ncomp),
     &                 Strain  (ncomp), dStrain (ncomp  ),
     &                 epsPl   (ncomp), prop    (nProp  ),
     &                 coords  (3),       rotateM (3,3),
     &                 defGrad (3,3),     defGrad_t(3,3),
     &                 tsstif  (2)
c
c***************** 用户定义部分 *************************************
c*                                                                  *
c*                                                                  *
c***************** 用户定义部分 *************************************
return
end

输入输出变量说明

usermat 输入变量及其数据类型：

  • matId：材料号；整型。

  • elemId：单元号；整型。

  • kDomIntPt：材料积分点号；整型。

  • kLayer：单元层号，仅用于分层单元，用于 usermat 的 18x 族单元中仅 SHELL181 有分层功能；整型。

  • kSectPt：单元截面号；整型。

  • ldstep：载荷步数；整型。

  • isubst：载荷子步数；整型。

  • nDirect：正应力或正应变分量的个数，与单元类型有关；整型。

  • nShear：剪应力或剪应变分量的个数，与单元类型有关；整型。

  • ncomp：应力分量或应变分量个数，ncomp= nDirect + nShear；整型。

  • nStatev：状态变量个数，通过 TB,STATE 命令定义；整型。

  • nProp：材料常数个数，通过 TB,USER 命令定义；整型。

  • Temp：增量步开始时的温度；双精度。

  • dTemp：增量步内的温度增量；双精度。

  • Time：增量步开始时的时间，但并非真正的时间，仅仅是计算过程的标记；双精度。

  • dTime：增量步内的时间增量；双精度。

  • Strain(ncomp)：增量步开始时的应变分量，包括机械应变及温度应变两部分；双精度数组。

  • dStrain(ncomp)：应变增量，仅包含机械应变，温度应变成分已被去除；双精度数组。

  • prop(nProp)：材料常数数组，由 TB,USER 及 TB,DATA 命令定义，可定义不同温度下的材料常数；双精度数组。

  • coords(3)：当前材料积分点的坐标；双精度数组。

  • rotateM(3,3)：大变形分析（NLGEOM,ON）时考虑刚体转动用到的增量型旋转矩阵，小变形分析时为单位矩阵；双精度矩阵。该矩阵老版本中存在，现已被废除（本贴第七部分提供了计算转动矩阵 R 的代码）。

  • defGrad_t(3,3)：增量步开始时的变形梯度矩阵；双精度矩阵。

  • defGrad(3,3)：增量步结束时的变形梯度矩阵；双精度矩阵。

usermat 输出变量及其数据类型：

  • keycut：载荷切分控制参数，0—不切分，1—切分，默认情况下不进行载荷切分，当 ANSYS 主程序检测到收敛困难时进行载荷切分，并输出 keycut=1；整型。

  • epsZZ：平面应力状态时垂直于平面方向的应变，仅在平面单元或壳单元考虑厚度变化时使用；双精度。

  • tsstif (2)：横观剪切刚度，tssif(1)-GXZ，tssif(2)-GYZ；双精度数组。

  • dsdePl(ncomp,ncomp)：一致切线算子矩阵，dsdePl(i,j) 表示增量步结束时由第 j 个应变分量的改变引起的第 i 个应力分量的变化。ANSYS 默认一致切线算子矩阵是对称的，若要使用非对称的一致切线算子矩阵（对于采用非关联流动法则的塑性材料模型，一致切线算子矩阵往往是非对称的），需要使用 NROP,UNSYM 命令打开非对称求解器；双精度矩阵。

usermat 输入输出变量及其数据类型：

  • stress(ncomp)：应力分量，增量步开始时由 ANSYS 主程序传入，需要在 usermat 中更新为增量步结束时的值；双精度数组。

  • statev(nStatev)：状态变量，用于存储与材料模型有关的变量，如塑性因子、强化参数、损伤变量等，增量步开始时由 ANSYS 主程序传入，需要在 usermat 中更新为增量步结束时的值；双精度数组。

  • epsPl(ncomp)：塑性应变分量，增量步开始时由 ANSYS 主程序传入，需要在 usermat 中更新为增量步结束时的值；双精度数组。

  • sedEl：弹性功，仅用于后处理，对计算过程无影响；双精度。

  • sedPl：塑性功，仅用于后处理，对计算过程无影响；双精度

  • epsep：等效塑性应变

APDL 命令自定义材料模型

通过 APDL 命令很容易使用自定义材料模型。只需在前处理部分使用下列语句即可： TB, USER, matId, NTEMPS, NPTS 关键字 USER 表示使用用户定义材料模型；关键字 matId 表示材料号；NTEMPS 为温度点的个数，即可定义不同温度下的多组材料常数；NPTS 表示每一温度下材料常数个数

线弹性本构 usermat

\[ \begin{aligned} \bf{\sigma} & = \bf{D}_e \bf{\varepsilon}\\ \bf{\sigma} & = \left[ {\sigma_x \; \sigma_y \; \sigma_z \; \sigma_{xy} \; \sigma_{yz} \; \sigma_{xz}} \right]^T \\ \bf{\varepsilon} & = \left[ {\varepsilon_x \; \varepsilon_y \; \varepsilon_z \; \varepsilon_{xy} \; \varepsilon_{yz} \; \varepsilon_{xz}} \right]^T \\ \\ \bf{D}_e & = \left[ {\begin{array}{*{20}{c}} {\lambda + 2G}&\lambda &\lambda &0&0&0\\ \lambda &{\lambda + 2G}&\lambda &0&0&0\\ \lambda &\lambda &{\lambda + 2G}&0&0&0\\ 0&0&0&G&0&0\\ 0&0&0&0&G&0\\ 0&0&0&0&0&G \end{array}} \right] \end{aligned} \]

以下是线弹性本构的usermat源代码及测试用的命令流。代码格式显示可能有误，仅供参考。下载源文件

usermat.F

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
*deck,usermat      USERDISTRIB  parallel                                gal
      subroutine usermat(
     &                   matId, elemId,kDomIntPt, kLayer, kSectPt,
     &                   ldstep,isubst,keycut,
     &                   nDirect,nShear,ncomp,nStatev,nProp,
     &                   Time,dTime,Temp,dTemp,
     &                   stress,statev,dsdePl,sedEl,sedPl,epseq,
     &                   Strain,dStrain, epsPl, prop, coords, 
     &                   var0, defGrad_t, defGrad,
     &                   tsstif, epsZZ,
     &                   var1, var2, var3, var4, var5,
     &                   var6, var7, var8)
c*************************************************************************

#include "impcom.inc"
c
      INTEGER          
     &                 matId, elemId,
     &                 kDomIntPt, kLayer, kSectPt,
     &                 ldstep,isubst,keycut,
     &                 nDirect,nShear,ncomp,nStatev,nProp
      DOUBLE PRECISION 
     &                 Time,    dTime,   Temp,    dTemp,
     &                 sedEl,   sedPl,   epseq,   epsZZ
      DOUBLE PRECISION 
     &                 stress  (ncomp  ), statev (nStatev),
     &                 dsdePl  (ncomp,ncomp),
     &                 Strain  (ncomp  ), dStrain (ncomp  ), 
     &                 epsPl   (ncomp  ), prop    (nProp  ), 
     &                 coords  (3),       rotateM (3,3),
     &                 defGrad (3,3),     defGrad_t(3,3),
     &                 tsstif  (2)
c


      DOUBLE PRECISION var0, var1, var2, var3, var4, var5,
     &                 var6, var7, var8
      
c
c************************用户定义部分*************************************
c     	
      INTEGER          K1, K2
      DOUBLE PRECISION ONE, TWO
	PARAMETER (ONE=1.0D0, TWO=2.0D0)
	DOUBLE PRECISION EMOD, ENU, EG, EG2, ELAM
	DOUBLE PRECISION SesSrn(6,6)
C
C 相关弹性常数
C
	EMOD  =  prop(1)
	ENU   =  prop(2)
	EG    =  EMOD/(TWO*(ONE+ENU))
	EG2   =  TWO*EG
	ELAM=EMOD*ENU/(ONE+ENU)/(ONE-TWO*ENU)
      DO K1=1,6
	    DO K2=1,6
	    SesSrn(K1,K2)=0.0D0
	    END DO
	END DO
	
C
C 获得应力应变矩阵
C
	DO K1=1, 3
	    DO K2=1, 3
	    SesSrn(K2, K1)=ELAM
	    END DO
	    SesSrn(K1, K1)=EG2+ELAM
	END DO
	DO K1=4, 6
	    SesSrn(K1 ,K1)=EG
	END DO
c	
C 更新应力
C
	DO K1=1, 6
	    DO K2=1, 6
	    stress(K2)=stress(K2)+SesSrn(K2, K1)*dStrain(K1)
	    END DO
	END DO
c
c一致切线算子矩阵
c
      DO K1=1,6
	    DO K2=1,6
	    dsdePl(K1,K2)=SesSrn(K1,K2)
	    END DO
	END DO
c
	
      return
      end

uniaxial tension_ANSYS.inp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
!ANSYS自带线弹性材料——单轴拉伸
finish
/clear
/prep7

!单元类型
et,1,45

!设置弹性模量及泊松比
ex,1,180e9
prxy,1,0.3 

!调整视角
/view,1,-1
/ang,1 

!建立几何模型
cyl4,0,0,0.005, , , ,0.1 

!网格划分
lsel, s , , , 1, 8
lesize, all, , ,4
allsel,all
lsel, s , , , 9, 10
lesize, all, , ,20
allsel,all
vsweep, 1, 1, 2, 0           !进行扫略网格划分

!设置边界条件
da, 1, all                        !约束其中一个端面
!施加外载荷
sfa, 2, 1, pres, -200e6      !钢棒另一端面施加拉伸载荷
finish

/solu                             !进入求解器
solve 

/post1                          !进入后处理器
pldisp, 1                     !变形图

uniaxial tension_user.inp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
!用户自定义线弹性材料——单轴拉伸
finish
/clear
/prep7

!单元类型
et,1,185

!选择用户自定义材料类型
tb,user,1,1,2
tbdata,1,180e9,0.3    !给弹性模量及泊松比赋值

!调整视角
/view,1,-1
/ang,1 

!建立几何模型
cyl4,0,0,0.005, , , ,0.1

!网格划分
lsel, s , , , 1, 8
lesize, all, , ,4
allsel,all
lsel, s , , , 9, 10
lesize, all, , ,20
allsel,all
vsweep, 1, 1, 2, 0         !进行扫略网格划分

!设置边界条件
da, 1, all                     !约束其中一个端面

!施加外载荷
sfa, 2, 1, pres, -200e6        !钢棒另一端面施加拉伸载荷
finish

/solu                         !进入求解器
time,1
pred,on                   !打开线性预测
nsubst,10,50,5          !定义载荷子步数
outres,all,all              !输出所有子步结果
solve

/post1                       !进入后处理器
pldisp, 1                    !变形图

usermat 的参考文档

下载