umatuel编程心得.doc

看到网上好多关于询问UMAT的一致切线模量的问题，以及UEL中的内力与外力残差问题。做子程序也有两年多，这里谈谈自己的一些对ABAQUS编程计算的看法，希望能受到指正批评。本文做静力学分析，所以有关动态和其他的偶合场不概讨论，但是如果掌握静态条件下，相信其他的外场应该比较容易推导。
ABAQUS隐式计算一切得从牛顿算法求解非线性方程组谈起。
首先得查阅数值牛顿算法的书：可以在书中发现这个公式。，这里面全部代表矩阵，可以核对一下。代表我们给定所有变量的一个初值，如果假设这里变量为位移的话，那就是了，代表三个方向的位移初值。为方程式对各个变量求偏导后再代入初值的jocobi矩阵，那么就代表其对应的求解方程式。牛顿解法必须设定一个初值，这样能计算出下一个值，一直循环直到迭代满足精度为止。这就是牛顿迭代算法的全部东西，其实没什么的！~呵呵。

那么来谈谈ABQUS里面的具体计算过程了，这里即适用于UEL二次编程，也适用于UMAT编程。有限元在时刻给定我们要在时刻求解
这里的内变量为个，偶合在各种方程中。有限元的计算关键得求得节点位移增量 ，a这里为节点号。（这里注意：所有的有限元包括线弹性，在ABAQUS刚开始都认为是非线性的，这样必须用增量的方法求得，所以有限单元法这本书中的前面部分并不适用，要看后面的增量有限元法，前面只是作为一个理念问题）。现在我们假设求得n+1时刻的位移增量，那么事情就好办了，单元内部所有的位移变量可以用插值函数插值得到，应变也可以对坐标求偏导得到，其他的量都可以随即求得。
下面这个方程式大家可以查阅有限单元法书籍的一个公式，即平衡方程的等效积分“弱形式”—虚功原理，其中一种为需位移原理，具体公式给出为：
………..（公式中全为矩阵）
如果单元的位移用节点位移增量插值得到，为原坐标的坐标点。那么我们可以得到应力平衡方程的弱形式（忽略体力作用）：

这个公式推导很简单的（自己推一下吧），只要将位移插值函数代入几何方程即可（附加一点：这里插值函数是总体坐标的函数，想要转化为局部坐标，那么还得进行坐标转化，具体不展开讨论。）。这里我们可以看出，其应力和外加载荷都是在时刻的值，为节点位移增量的函数，而其实是常数。其实如果仔细看的话，我们就可以看出来了，这就是所谓牛顿迭代算法的方程式。方程式是节点位移增量的函数，第一项为函数，第二为常数嘛。其实ABAQUS中，或UEL中的内力和外力就指这两项， 为内力，而为外力，这里指平衡下的内力和外力。

有了方程式，我们可以使用刚开始介绍的牛顿迭代方法来计算这个偏微分方程了。现在我们要求的是未知变量是节点位移增量。 ABAQUS会自动给定一个初始节点增量位移值，，编UEL的我们可以知道，为什么编程开始，我们直接可以得到位移增量值，其实这个值是程序自动给定的，我们可以这样认为，这个值为假值，只是牛顿迭代为了得到下一个值而设定的初始值。牛顿迭代公式为：
，这里我把节点号舍去了。

这里R为残差，因为内力现在只是下计算的，所以会出现不平衡，UEL中我们可以知道，ABAQUS需要编程者给定一个内力与外力的残差，就是这里的残差。
下面最难的是求得这个在下的jacobi矩阵。
，这里代表在方程式求偏导后代入其初始值。因为不含位移项了，所以求导为零。所以：
，求到这里，我要特别介绍一下关于一致切线模量与连续切线模量的概念。看过塑性力学的人会知道，我们利用公式可以求得其
，这里的叫连续切线模量，其实换句话说，这里是在内力和外力达到平衡时求得的切线模量。
但是编UMAT的，我们知道，ABAQUS需要我们输入一致切线模量，也就是程序中的DDSDDE矩阵。其实在这里的意思就是，ABAQUS需要计算出任何位移增量下的应力与应变关系。如果有了这个，那么上面的那个公式可以往下进行了，

好了，这样子全部推导完毕了，我们把上面的牛顿迭代整理一下：
，
令为
，代表jacobi矩阵，c代表其位移插值，R就是内力与外力残差。
ABAQUS自动把作为下一次迭代的初始值，直到满足精度要求。关于精度要求，任何ABAQUS非线形算法中都有，一个是满足位移精度，一个是满足力平衡精度。
个人总结：ABAQUS中ULE最重要及最难的点，即定义其内力与外力的残差R，ABAQUS的外力其实程序自动划分总载荷给定了，所以编UEL时，我们不需要定义的外力，只需要定义内力即可以。而UAMT中重点要定义DDSDDE一致切线模量。这些都在上述公式推导过程有所涉及，希望对编程的人有所帮助。