基于ANSYS二次开发的高阶连续梁数值模拟★

殷玉沉 朱彦涛 孙玉周

(1.中原工学院，河南 郑州 450007； 2.西藏职业技术学院，西藏 拉萨 850000)

0 引言

高阶连续理论能够反映出特殊本构反应，例如在偶应力理论、应变梯度连续理论中引入尺度因子描述高阶连续结构的局部弯曲效应[1-3]。在研究高阶连续结构过程中，尺寸效应的影响日益被重视，但是传统有限元法构造高阶形函数需要大量的工作[4]，而且没有引入尺度因子不易实现数值模拟。采用传统有限元法分析计算，会遇到网格划分困难以及近似函数不连续等问题，但应用适当的基函数和权函数，移动最小二乘法(MLS)可以比较容易地构造出高阶形函数[5,6]。近年来，一些学者对相关问题进行了研究，李雷等[7]将无网格法和高阶理论相结合，对岩土工程中的一些问题进行数值模拟；李莉等[8]利用虚功原理建立各向异性修正偶应力理论并引入不同的材料细观参数,建立适用于层合板/夹层板的偶应力理论模型；孙玉周等[9，10]使用无网格法研究基于高阶连续理论下的碳纳米管的多尺度效应以及弯曲梁等问题。上述的研究主要集中在通过无网格法与高阶理论相结合，以及考虑结构尺度效应影响等方面，但是将高阶连续理论应用在有限元分析中以及将其在有限元软件中实现数值模拟的研究较少。

ANSYS作为有限元领域的大型通用软件，被世界各工业领域广泛接受。为适应具体工程需求，ANSYS还提供了APDL,UPFs,UIDL(用户界面设计语言)及TckTk(工具命令语言/图形开发工具箱)4种二次开发工具[11]。针对本文提出的问题，由于ANSYS在网格划分以及高阶连续结构分析上存在局限性，需要对ANSYS软件进行二次开发。本文采用移动最小二乘法可以方便的构造高阶连续形函数的优点，建立高阶连续梁模型，通过Fortran语言编写用户单元子程序，进行ANSYS的二次开发，实现高阶连续梁在有限元软件中的数值模拟，并引入尺寸因子对高阶连续梁进行分析。

1 基于MLS的计算框架

在传统弯曲梁(如图1所示)的基础上，高阶连续梁模型的应力和应变分别定义为：

应变：

(1)

应力：

σxx=Eεxx+lxεxxx
τxxx=g2Eεxxx+lxEεxx
τyxx=g2Eεyxx

(2)

其中,E为弹性模量；w为梁的挠度；g为体本征尺度因子；lx为面本征尺度因子。

梁在应变能及外力的作用下，根据虚功原理可得:

(3)

(4)

令m=4,则一维三次基函数为:

aT(x)=(1,x,x2,x3)

(5)

(6)

其中，wI(x)=wI(x-xI)为节点xI处的权函数，当权函数wI(x)在节点支撑域范围内时，wI(x)>0，否则，wI(x)=0。

由式(3)得：

j=1,2,3,…，m

(7)

由此得:

(8)

即:

A(x)b(x)=B(x)u

(9)

令：

B(x)=[wI(x)a(xI)w2(x)a(x2)…wN(x)a(xN)]

(10)

则形函数为：

(11)

权函数ωI(x)在节点xI处的值最大，而且支撑域边界附近趋于0，即当r=(x-xi)/dmI>1时，w(r)=0，dmI为节点支撑域的直径。

采用三次样条权函数，综上可求得形函数φ(x)的导数为:

φI,x(x)=H,x(x)B(x)+H(x)B,x(x)l
φI,xx(x)=[H,xx(x)B(x)+2H,x(x)B,x(x)+H(x)B,xx(x)]l
φI,xxx(x)=[H,xxx(x)B(x)+3H,x(x)B,xx(x)+3H(x),xxB,x(x)+H(x)B,xxx(x)]l

(12)

在无网格法的计算框架下，挠度的二阶和三阶导数近似为:

(13)

由式(12),式(13)代入式(3)可得高阶连续梁的总刚度矩阵表达式：

(14)

由式(14)可知，高阶连续问题中出现高阶位移导数，传统有限元法处理高阶形函数的工作量比较大，利用移动最小二乘近似(MLS)构造高阶形函数，位移高阶导数可以直接通过节点位移插值得到。

2 ANSYS二次开发

基于ANSYS为用户提供的二次开发功能，通过UPFs的特性实现ANSYS软件模拟一维高阶连续梁。由于ANSYS中传统弯曲梁缺少高阶连续理论，下面将含有的高阶项从高阶连续梁中分离出来即K1。通过添加附加项K1到用户自定义程序UELMatx.F中，实现高阶连续理论在有限元软件中的应用。

从式(14)得：

(15)

ANSYS的用户子程序UELMatx用于访问单元矩阵和载荷向量，主要目的是修改或监视单元刚度矩阵和荷载向量。程序编译过程中需要解决的关键问题主要在以下方面：

1)ANSYS二次开发编译链接环境建立；

2)利用APDL(参数化设计语言)创建节点支撑域覆盖范围内的单元；

3)将附加项K1通过编译用户自定义程序UELMatx实现与有限元软件结合。

ANSYS源程序是基于FORTRAN语言编写，为方便程序编译链接，采用FORTRAN语言对程序进行编译开发。本文基于ANSYS15.0进行二次开发，程序编译链接过程如下:

1)将编译完成的子程序UELMatx.F放置在ANSYS Incv150ansyscustomuserwinx64文件夹中，使用该文件夹下的ANSCUST.BAT批处理文件，编译连接成功后生成用户自定义ANSYS.exe文件；

2)激活UPFs：打开ANSYS15.0>Mechanical APDL Product Launcher 15.0，选择用户自定义生成的ANSYS.exe文件路径，也可通过GUI操作进行。该用户子程序在求解之前调用，调用命令如下：USRCAL,UELMATX。

3 算例分析

3.1 简支梁的计算分析

如图1所示，利用简支梁的例子来验证编译开发的用户子程序UELMatx在ANSYS中调用求解时的收敛性及精确度。简支梁受集中力F=1 N作用在跨中，梁的长度为10 mm，梁截面高度h=1 mm，厚度b=0.01 mm，梁的横截面面积A=h×b，弹性模量E=1.0×105MPa。

本文利用Beam3单元建立一维梁模型，为验证高阶连续梁在ANSYS二次开发后的数值计算精度，与基于无网格法的一维高阶连续梁的数值计算结果进行对比。梁上布置不同的节点数，梁中点挠度变化见表1。

表1 不同节点数时梁中点挠度值

假设一维简支梁的面本征尺度和体本征尺度为固定值。由表1的计算结果可以看出，不同节点数的一维高阶连续简支梁通过ANSYS二次开发模块求解的中点挠度值收敛且基本与无网格法框架下的数值计算结果吻合。

从表1的数据可以知道高阶连续理论下的简支梁中点挠度值与ANSYS标准状态下的结果存在极小差异。利用ANSYS二次开发后的计算模块分析尺度因子的变化对梁的挠度影响，相关实常数不变，只改变尺寸因子的数值，计算结果分别见表2,表3。

表2 尺寸因子g变化时对高阶连续梁中点挠度值

表3 尺寸因子lx变化时对高阶连续梁中点挠度值

由表2,表3可知，当尺寸因子g为定值，面本征尺寸因子lx从0.1变化到23，随着lx值的逐渐增大，对高阶连续梁的挠度值影响极小；当面本征尺寸因子lx不变，g从0.5变化到15.0，对梁的挠度变化有较大影响，其影响趋势如图2所示。

3.2 悬臂梁的计算分析

通过导入ANSYS二次开发程序后的计算模块分析当lx保持不变时尺寸因子g对高阶连续悬臂梁的自由端挠度的影响。构建一维悬臂梁，集中力P=10 N作用在悬臂梁自由端，梁长度为5 mm，截面高度h=1 mm，厚度b=0.01 mm，横截面面积A=h×b，弹性模量E=1.05×105MPa(面本征尺寸因子lx=2.5)，见图3。

表4 尺寸因子g变化时对悬臂梁自由端挠度值 mm

g0.51.01.52.55.07.510.012.515.0w4.76344.75894.75194.72904.62464.46094.25064.00833.7478注:g为体本征尺寸因子;w为悬臂梁自由端挠度

由表4数据可以看出，导入ANSYS二次开发的程序后，高阶连续悬臂梁自由端挠度受到体本征尺寸因子g值的变化影响，且悬臂梁自由端挠度值变化趋势随g值的增大而减小(如图4所示)。对比图2与图4可知梁挠度值的变化趋势基本一致。

4 结语

基于ANSYS为用户提供的APDL和UPFs二次开发工具，通过研究高阶连续理论在有限元软件中的应用，得出以下结论：1)利用ANSYS提供的UELMatx.F用户子程序，将高阶连续理论及尺寸因子的概念引入到ANSYS中，实现高阶连续梁在有限元软件中的应用，克服有限元软件对高阶连续结构数值模拟的局限性，并得出尺寸因子的影响趋势。2)计算结果与基于移动最小二乘法框架的Fortran数值计算结果基本一致，通过算例验证了该方法的合理性、可行性。3)证明该方法对移动最小二乘法在有限元软件中应用的可行性，为高阶连续理论在有限元软件中的应用提供新思路。

[1] Xia Z C,Hutchinson J W.Crack tip fields in strain gradient plasticity[J].Journal of the Mechanics & Physics of Solids,1996,44(10):1621-1648.

[2] Fleck N A,Hutchinson J W.Strain Gradient Plasticity[J].Advances in Applied Mechanics,1997(33):295-361.

[3] Toupin R A.Elastic materials with couple-stresses[J].Archive for Rational Mechanics and Analysis,1962,11(1):385-414.

[4] 梁醒培,王 辉.应用有限元分析[M].北京：清华大学出版社,2010.

[5] Belytschko T,Krongauz Y,Organ D,et al.Meshless methods:An overview and recent developments[J].Computer Methods in Applied Mechanics & Engineering,1996,139(1-4):3-47.

[6] Belytschko T,Lu Y Y,Gu L.Element-Free Galerkin methods[J].International Journal for Numerical Methods in Engineering,1994,37(2):229-256.

[7] 李 雷,谢水生,黄国杰.应变梯度塑性理论下超薄梁弯曲中尺度效应的数值研究[J].工程力学,2006,23(3):44-48.

[8] 李 莉,陈万吉,郑 楠.修正偶应力理论层合薄板稳定性模型及尺度效应[J].工程力学,2013,30(5):1-7.

[9] Sun Y,Liew K M.The buckling of single-walled carbon nanotubes upon bending:The higher order gradient continuum and mesh-free method[J].Computer Methods in Applied Mechanics & Engineering,2008,197(33-40):3001-3013.

[10] Sun Y,Liew K M.Application of the higher-order Cauchy-Born rule in mesh-free continuum and multiscale simulation of carbon nanotubes[J].International Journal for Numerical Methods in Engineering,2010,75(10):1238-1258.

[11] 师 访.ANSYS二次开发及应用实例详解[M].北京：中国水利水电出版社,2012.

[12] 张 雄,刘 岩.无网格法(精)[M].北京：清华大学出版社,2004.