（原创，转载请注明出处）

==概述==

本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式。有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。

一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元求解器，通过自研求解器和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。

iSolver介绍视频：

http://www.jishulink.com/college/video/c12884

==第20篇：UEL用户子程序开发步骤==

用户子程序主要是将用户特定的材料本构模型和单元算法等公式编写为计算机语言表示的公式，并实现和商软求解器之间的交互迭代。

常用的商业有限元软件都提供了用户自定义子程序的功能，且一般都是Fortran语言开发，Fortran是上世纪70年代的语言，相对现代化的流行语言编写，格式要求非常严格，编译调试都比较繁琐，使得开发效率低下，而且接口限制较多，除了商软提供的功能外用户基本没法改动，灵活性较差。由于用户子程序很多都涉及复杂的公式编写，用户除了需要扎实的理论基础外，还需要较强的能将公式表达为Fortran语言的编程能力，这对非计算机专业出身的人来说往往在浪费了很多额外精力，使得很多理论高手都对用户子程序望而却步，难以入门。

在实际工作中，很多工程师用Matlab来编写和推导公式，Matlab被认为是市面上最接近草稿纸上推导公式的一款软件了，而且有限元在数值层面上的计算其实就是矩阵运算，所以Matlab这种数据按矩阵来组织非常适合用来开发有限元相关的程序。而现在市面上还没有采用Matlab来开发商软子程序的案例。iSolver是市面上第一款基于Matlab来开发商软用户子程序的软件工具，支持用Matlab编写和调试用户子程序。iSolver子程序的接口完全按照Abaqus的标准实现，而Abaqus的子程序接口在近几年内已经基本不再变化了，同样的，虽然iSolver在不断发展，但iSolver子程序接口将维持不变，所有在iSolver上编写的算法子程序都只要维护自己的算法部分就行，而不是维护整个有限元求解的整个过程。

前面第八、九篇介绍了UMAT用户自定义材料的开发，这里将介绍UEL用户自定义单元的开发，本文首先简单的讨论了UEL的一般含义，并详细的介绍了基于Fortran和Matlab两种方式的UEL的开发步骤，对比发现开发步骤基本相同，同时采用Matlab更加高效和灵活。具体的开发过程可以参考下面我们的Step by Step的录像，包括了整个的有限元基础理论和我们对Abaqus中的单元理解：

https://www.jishulink.com/college/video/c14948

深入浅出有限元：基础理论->Abaqus操作->matlab编程

                                             

1.1 UEL的关键输入输出参数

UEL网上资料很多，大家可以很容易查看，但大部分资料都只是简单提供UEL算例，这里我们列出了UEL接口的关键输入输出参数，如下表所示：

1.2 基于Fortran的Abaqus的UEL的开发步骤

1.2.1 在inp文件中定义UEL

Abaqus中只能通过修改手动inp文件完成用户UEL的定义，通常包含以下关键字及相应属性，如图所示：

• *User element关键字，用于定义单元基本信息，包括以下属性：

  nodes为单元节点数；

  type为单元类型，U表示为用户自定义单元，1001表示单元类型编号，一般从1000开始；

  properties为浮点数属性个数，一般用于输入单元对应的section数据和材料数据；

  coordinates为坐标维数；

  variables为状态量个数，一般求解计算过程中需要在迭代步之间传递的变量；

  第二行开始定义节点激活的自由度；

• *Uel property关键字，与*User element关键字中的properties对应，第二行开始定义具体的数据。

1.2.2 编写

使用任意编辑器编写.for文件，推荐使用Visual Studio Code，微软开源的轻量化代码编辑器，配置灵活高效：

1.2.3 编译（可选）

Abaqus没有自带Fortran编译器，所以用户需要自己去安装Fortran编译器和Visual Studio Build Tools，并配置相应环境。具体配置过程与UMAT一致，可以查看我们关于环境配置的视频：

https://www.jishulink.com/college/video/c13034?chapter=1

在环境配置完成之后，打开命令提示框，输入命令Abaqus make Library=XXX.for，即开始编译，编译过程中的警告和错误都会打印在命令提示框内。

1.2.4 运行

运行有两种方法，第一种就是在命令提示框中输入Abaqus job=XXX user=XXX.for，如下图所示。

第二种就是在Abaqus中创建基于inp文件的任务，然后选择对应的用户子程序for文件，在任务管理器中提交运行，如图所示。

至此，基于Fortran的UEL开发流程已经完成，但结果的正确性还需要更加细致的验证，为更方便的查找问题，建议先采用单个单元调试UEL，在确保单个单元正确后再将UEL用于实际问题。

1.2.5 调试（可选）

如果想要知道代码的运行结果是否和预期的一致，一种笨办法是用print打印到log文件中，高效的方法是采用断点调试的方法进行运行中的调试。

Abaqus支持命令行调试，不过命令行反复运行也比较繁琐，用户也可选择用一键调试Abaqus的用户子程序的DUS插件工具。DUS（Debug User Subroutine）是集成在ABAQUS/CAE中的一个插件，能够一键启动用户配置的用户子程序开发平台（如Visual Studio 2008等），并进入对用户定义子程序的单步调试模式。

有兴趣的可到下面网页下载使用。

Abaqus用户子程序调试插件：

https://www.jishulink.com/content/post/424513

1.3 基于Matlab的iSolver的UEL开发步骤

基于Matlab的Abaqus的UEL具体开发步骤和Abaqus类似，只不过某些步骤需要用到自研有限元求解器开发平台iSolver。

1.3.1 在inp文件中定义UEL

与Abaqus相应的操作一致，如图所示：

1.3.2 编写

在Matlab中创建并编写U1001.m的文件，放入Abaqus工作目录下。该文件只包括一个U1001函数，接口和Abaqus的接口参数完全一致，功能也是计算应力应变关系和当前应力状态等，相对Fortran，利用Matlab可以更容易的编写计算公式，同时可以利用Matlab在矩阵计算中各种强大功能和算法库。因为Abaqus的UEL接口和计算功能各个版本相对固定，这个matlab的UEL接口参数也相对固定，不会因为iSolver的版本不同而重新修改接口。

1.3.3 编译（无）

由于matlab是脚本语言，不需要编译。

1.3.4 调试（可选）

在Abaqus菜单栏的Plug-ins里选择iSolver插件的菜单。

点击iSolver->Engine，按照下图所示，在功能项Use Solver中选择iSolver，在Source Type里面选择Matlab，勾选Debug。点击Submit进行调试运行。

程序会自动打开matlab并加载U1001.m文件，手动打上断点

点击在Debug菜单下的Run U1001运行。

程序将在断点处停止，且将鼠标移动到需要调试查看的参数上，能够查看到对应的值。

按F10可以进行单步调试。

1.3.5 运行

在上述步骤的基础上去掉勾选Debug选项，点击Submit运行计算，此时将采用iSolver求解器联合U1001.m进行求解分析，运行完毕点击Result在Abaqus中查看结果。

1.4 算例

具体的壳理论和在iSolver中实现的UEL算例可以参考下面的视频：

https://www.jishulink.com/college/video/c14948

深入浅出有限元：基础理论->Abaqus操作->matlab编程

 

1.5 总结

本文首先简单的讨论了UEL的一般含义，并详细的介绍了基于Fortran和Matlab两种方式的UEL的开发步骤，对比发现开发步骤基本相同，但Matlab更加高效和灵活。同时，由于iSolver基本单元类型和Abaqus算法完全一致，可以发现同一个算例验证两者分析结果完全一致，从而证明基于Matlab的UEL的流程和结果的正确性。

UEL的开发一方面要有扎实的公式推导能力，另一方面需要基础的编程能力和开发工具应用水平，后者不是重点，但往往浪费了大家很多的精力，善用工具方能提高效率，基于Fortran和Matlab两种方式的UEL的开发步骤和开发工具如下表：

 

如果有任何其它疑问或者项目合作意向，也欢迎联系我们：

SnowWave02 From www.jishulink.com

email: snowwave02@qq.com

以往的系列文章：

第一篇：S4壳单元刚度矩阵研究。介绍Abaqus的S4刚度矩阵在普通厚壳理论上的修正。

http://www.jishulink.com/content/post/338859

第二篇：S4壳单元质量矩阵研究。介绍Abaqus的S4和Nastran的Quad4单元的质量矩阵。

http://www.jishulink.com/content/post/343905

第三篇：S4壳单元的剪切自锁和沙漏控制。介绍Abaqus的S4单元如何来消除剪切自锁以及S4R如何来抑制沙漏的。

http://www.jishulink.com/content/post/350865

第四篇：非线性问题的求解。介绍Abaqus在非线性分析中采用的数值计算的求解方法。

http://www.jishulink.com/content/post/360565

第五篇：单元正确性验证。介绍有限元单元正确性的验证方法，通过多个实例比较自研结构求解器程序iSolver与Abaqus的分析结果，从而说明整个正确性验证的过程和iSolver结果的正确性。

https://www.jishulink.com/content/post/373743

第六篇：General梁单元的刚度矩阵。介绍梁单元的基础理论和Abaqus中General梁单元的刚度矩阵的修正方式，采用这些修正方式可以得到和Abaqus梁单元完全一致的刚度矩阵。

https://www.jishulink.com/content/post/403932

第七篇：C3D8六面体单元的刚度矩阵。介绍六面体单元的基础理论和Abaqus中C3D8R六面体单元的刚度矩阵的修正方式，采用这些修正方式可以得到和Abaqus六面体单元完全一致的刚度矩阵。

https://www.jishulink.com/content/post/430177

第八篇：UMAT用户子程序开发步骤。介绍基于Fortran和Matlab两种方式的Abaqus的UMAT的开发步骤，对比发现开发步骤基本相同，同时采用Matlab更加高效和灵活。

https://www.jishulink.com/content/post/432848

第九篇：编写线性UMAT Step By Step。介绍基于Matlab线性零基础，从零开始Step by Step的UMAT的编写和调试方法，帮助初学者UMAT入门。

http://www.jishulink.com/content/post/440874

第十篇：耦合约束（Coupling constraints）的研究。介绍Abaqus中耦合约束的原理，并使用两个简单算例加以验证。

https://www.jishulink.com/content/post/531029

第十一篇：自主CAE开发实战经验第一阶段总结。介绍了iSolver开发以来的阶段性总结，从整体角度上介绍一下自主CAE的一些实战经验，包括开发时间预估、框架设计、编程语言选择、测试、未来发展方向等。

http://www.jishulink.com/content/post/532475

第十二篇：几何梁单元的刚度矩阵。研究了Abaqus中几何梁的B31单元的刚度矩阵的求解方式，以L梁为例，介绍General梁用到的面积、惯性矩、扭转常数等参数在几何梁中是如何通过几何形状求得的，根据这些参数，可以得到和Abaqus完全一致的刚度矩阵，从而对只有几何梁组成的任意模型一般都能得到Abaqus完全一致的分析结果，并用一个简单的算例验证了该想法。

http://www.jishulink.com/content/post/534362

第十三篇：显式和隐式的区别。介绍了显式和隐式的特点，并给出一个数学算例，分别利用前向欧拉和后向欧拉求解，以求直观表现显式和隐式在求解过程中的差异，以及增量步长对求解结果的影响。

http://www.jishulink.com/content/post/537154

第十四篇：壳的应力方向。简单介绍了一下数学上张量和Abaqus中壳的应力方向，并说明Abaqus这么选取的意义，最后通过自编程序iSolver来验证壳的应力方向的正确性。

https://www.jishulink.com/content/post/1189260

第十五篇：壳的剪切应力。介绍了壳单元中实际的和板壳近似理论中的剪切应力，也简单猜测了一下Abaqus的内部实现流程，最后通过一个算例来验算Abaqus中的真实的剪切应力。

https://www.jishulink.com/content/post/1191641

第十六篇：Part、Instance与Assembly。介绍了Part、Instance与Assembly三者之间的关系，分析了Instance的网格形成原理，并猜测Abaqus的内部组装实现流程，随后针对某手机整机多part算例，通过自编程序iSolver的结果比对验证我们的猜想。

https://www.jishulink.com/content/post/1195061

第十七篇：几何非线性的物理含义。介绍了几何非线性的简单的物理含义，并通过几何非线性的悬臂梁Abaqus和iSolver的小应变情况的结果，从直观上理解几何非线性和线性的差异。

https://www.jishulink.com/content/post/1198459

第十八篇：几何非线性的应变。首先从位移、变形和应变的区别说起，然后通过一维的简单例子具体介绍了几何非线性下的应变的度量方式，并给出了工程应变、 真实应变、Green应变三者一维情况下在数学上的表达方式。

https://www.jishulink.com/content/post/1201375

第十九篇：Abaqus几何非线性的设置和后台。首先介绍了几何非线性一般的分类，然后详细说明了Abaqus中几何非线性的设置方式和常用单元的分类，最后以一个壳单元的简单算例为对象，可以发现应变理论、Abaqus和iSolver三者在线性、小应变几何非线性和大应变几何非线性三种情况下都完全一致，从而验证Abaqus几何非线性后台采用的应变和我们的预想一致。

http://www.jishulink.com/content/post/1203064