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

==概述==

本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过

（1）   基础理论

（2）   商软操作

（3）   自编程序

三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。

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

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

自主结构有限元求解器iSolver介绍视频：

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

==第25篇：显式分析的稳定时间增量==

相对隐式分析，显式分析在单元函数中的计算要相对简单，只需计算应力的更新，而少了单元刚度矩阵的求解，但显式在增量步的自动步长、网格畸变的控制、质量缩放等方面又比隐式要多做许多工作。其中求解过程中的稳定时间增量是每个显式分析都会遇到的问题，由于没有迭代，显式分析都能得到某个结果，不存在收敛问题，而结果的正确性和稳定时间增量密切相关，很多人做显式分析都一般直接交给商软去做自动步长，当商软显式结果发散或者结果差异很大时也不清楚怎么修改时间步长。同时，对于显式自定义单元VUEL等子程序，Abaqus要求用户自己计算稳定时间增量dtimestable，Abaqus会根据用户计算的稳定时间增量来限制它的增量步长。本文将简单介绍一下稳定时间增量的概念和理想及工程应用上的两种计算方式，并用Abaqus中一个简单的算例来验证工程上的稳定时间增量的计算公式，便于你对稳定时间增量的理解和自己编程实现。

1.1.1 稳定性的含义

在本系列13篇：显式和隐式的区别中提到，显式分析都是条件稳定的，譬如下面求解一个微分方程：

                                                                        y'(x) = -y+x+1

其中y(0) = 1。显式分析可得到下面的增量表达式：

当h<2时，y收敛，但h>=2时，y将发散。当h=2.1时，可发现显式分析的蓝色线将和理论值越来越远。

上面针对的是一个具体的数学函数，对一个任意的有限元动力学系统，只要采用显式分析，也同样存在稳定性问题。动力学问题如果把质量导致的惯性力考虑在内，也是在解下面的平衡方程：

如果上式所有量都是t或者所有量都是t+dt，那么上式肯定是平衡的，没有一点问题。但在有限元求解过程中，只能通过已知时刻点t来求未知时刻点t+dt的值，而某些值在t+dt的值是不知道的，譬如刚度矩阵K，所以，只能退而求其次，上式中部分是时刻t的，譬如刚度矩阵K，而部分是时刻t+dt的值，譬如载荷F和质量阵M。

这么取值后上式左边就不再=0了，产生了人为导致的误差。而这种人为的误差，如果是隐式分析，那么可以通过迭代来解决，迭代多次后最后的时刻的Tol可以达到预定的小量，方程依然是平衡的。但对显式分析，流程如下所示，不会再次迭代，只会进行下一次的t+2*dt的求解，这样，就可能导致两种结果：

（1）   t时刻的误差将在t+dt时刻累积，此时系统就是不稳定的。譬如上图中的蓝色曲线。

（2）   t时刻的误差在t+dt时刻不累积，此时系统就是稳定的。譬如上图中的如果h<2。

1.1.2 稳定时间增量的理想计算方式

总结显式分析中增量步之间的关系，可以发现增量步和上一个增量步有简单的关系：

那么当前时刻和第一个时刻的一个简单关系：

而这个C类似于一个比例系数，与dt的取值有关，譬如上述函数例子中C=1-dt。对这种只有一个自由度的问题，可以发现|C|<1时，系统稳定，否则系统不稳定。在数学上可以证明C与系统的特征值是密切相关的。也就是说dt与系统的特征值有关，对一个有限元系统来说，数学上已经证明dt如果满足下式就是稳定的：

其中

为有限元系统的最大特征值。我们称右端为稳定时间增量的理想计算方式。

1.1.3 稳定时间增量简化的必要条件

在稳定时间增量dt_ideal的理想计算方式中，需要计算有限元系统的最大特征值。而如果每次显式分析前都要计算一遍最大特征值，且对几何变化大的问题，在每个增量步都要更新一次，那么计算量将非常大。更何况对非对称的刚度等系统，还没有办法得到特征值。在实际有限元的工程计算中，会找另一种简单的稳定时间增量dt_engeer的工程计算方式，这个dt_engeer必须满足两个条件：

（1）   条件1：这种计算方式并不一定严格求出上面的dt_ideal，但求出的稳定时间一定比理想的稳定时间增量更小。即当有限元的时间增量步长dt<dt_engeer，必然有dt<dt_ideal，此时系统依然满足稳定性要求。

（2）   条件2：dt_engeer比dt_ideal不能小太多，应该是同一量级，否则增量步将大幅增加，导致系统耗时过大。

1.1.4 一般有限元的稳定时间增量的工程计算方法

本着这个原则，一般有限元上做了两次计算简化工作：

（1）第一步：是将整个系统的最大特征值取为所有单个单元的特征值,由于系统的约束会压缩总体频率，导致dt_element<dt_ideal。也就是：

（2）单个单元依然有特征值求解问题，因此，再进一步简化，实际取的稳定时间增量不再需要计算特征值，而是估算为：

其中min表示对所有单元遍历后的最小值，Le为单元特征长度，Cd为材料疏密波速度，对杆件有Cd=sqrt(E/d), E、d分别为杆长度、杨氏模量和密度。按此公式计算将非常简单，只需遍历所有单元，将得到的参数代入上式即可。

为了验证上述简化满足前面所述两点，我们举一个由杆单元truss组成的简单系统来验证一下dt_engeer的具体值。此时，每个杆单元的频率可以简单的估算为

其中L、E、d分别为杆长度、杨氏模量和密度。注意上式我们用了集中质量。此时的第一步稳定时间增量

显然，dt_engeer和dt_element同一量级，且略小。

由于dt_engeer<dt_ideal，可以发现有限元商软采用的dt_engeer不一定是最优的，如果你能找到一个时间增量步dt，使得dt_engeer<dt<dt_ideal，那么可以比商软默认的耗时更少。后面的例子中也证明了这点。

1.1.5 Abaqus对稳定时间增量的工程计算方法的修正

Abaqus对稳定时间增量的工程计算还做了两处不同的修正：

（1）   考虑Bulk viscosity效应修正工程稳定时间增量,主要是利用Linear Bulk Viscosity系数b1减小dt_engeer。即

b1在Step->Other的Linear Bulk Viscosity可修改。

默认为0.06，此时

（2）   为了避免系统对实数的精度或者截断误差，加了一个Tolerance，使得

显然这个Tolerance是个远小于1的值，Abaqus取为0.01，

• 此值我们没发现在Abaqus界面上怎么修改，如果谁能找到，也希望能交流一下。

1.2 Abaqus的实现验证

我们将在Abaqus中采用一个简单的显式分析算例，来验证两个问题：

（1）   Abaqus采取的稳定时间增量和上述最后的稳定时间增量的计算公式一致。

（2）   对于某些问题，其实Abaqus采取的稳定时间增量是dt_engeer，相对保守，实际上可以取一个更加接近dt_ideal的值，系统依然是稳定的。

1.2.1 模型例子

模型为一根长1m，截面为4平方cm的圆柱，左端5cm由Steel组成，杨氏模量和密度分别是2e11Pa和7800千克每立方米，剩下部分由碎石组成，杨氏模量和密度分别是4.432e6Pa和1560千克每立方米。右端点固定在墙上，左端点面载荷4e6sin(150t)。我们只研究0-0.01s时间内的位移。

在Abaqus中建模如下，我们简单将模型划分为20个单元。采用truss单元。

1.2.2 稳定时间增量的理论值

1.2.2.1 稳定时间增量的理想计算的理论值

理想计算方式需要先计算系统最大模态特征，由于是20个单元，采用truss单元，就相当于只有21个自由度，右端约束后，无约束的自由度为20个，得到的K和M矩阵的秩为20，那么无论用哪种模态计算方法，得到的模态最大为20阶。在Abaqus中计算，结果如下，可得20阶的模态频率为30864Hz。

1.2.2.2 稳定时间增量的工程计算的理论值

最小的工程稳定时间增量显然是左端的Steel单元，此时为：

1.2.3 自动步长

在Abaqus中选择显式分析，dynamic,explicit，同时设置为自动步长。

运行结束后查看.sta文件，Abaqus会在此文件中在第一个增量前记录前十个最小的单元稳定时间增量。可发现如下所示，最小的单元稳定时间增量为第一个单元，且值和理论完全一致：

此时总共增量步为1087次，得到的左端位移随时间的变换曲线如图：

1.2.4 固定步长

Abaqus中改为固定步长：

取固定步长分别为dt2=1e-5和dt3=1.06e-5，即

                                                             dt_engeer<dt2<dt_ideal<dt3

得到的左端位移随时间变化曲线如下：

可发现dt3已经发散，而dt2和自动步长基本一致，但dt2只计算了1000个增量步，比自动步长少了8.7%。从而验证了，对某些问题，为了加速，时间增量可以取一个超过Abaqus默认更加接近dt_ideal的值，系统依然是稳定的。

1.3 联系方式

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

snowwave02 From www.jishulink.com

email: snowwave02@qq.com

以往的系列文章：

1.4.1 ========第一阶段========

第一篇：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

1.4.2 ========第二阶段========

第十一篇：自主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

第二十篇：UEL用户子程序开发步骤。本文首先简单的讨论了UEL的一般含义，并详细的介绍了基于Fortran和Matlab两种方式的UEL的开发步骤，对比发现开发步骤基本相同，但Matlab更加高效和灵活。

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

1.4.3 ========第三阶段========

第二十一篇：自主CAE开发实战经验第二阶段总结。从实战角度介绍自主CAE在推广和工程化应用的过程中的体会，同时说明一个CAE平台最重要的两个特点：可扩展和易维护。

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

第二十二篇：几何非线性的刚度矩阵求解。介绍几何非线性下的刚度矩阵的理论推导和计算机求解方法，最后利用一个简单的算例验证我们对Abaqus几何非线性的刚度矩阵的实现方式的猜测。

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

第二十三篇：编写简单面内拉伸问题UEL Step By Step。通过简单面内拉伸问题UEL的编写，介绍解决有限元问题从理论到算法再到编程实现的一般流程以及面内拉伸问题的基本算法步骤，最后将UEL与Abaqus的S4R单元计算结果进行对比，进行验证。

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

第二十四篇：显式求解Step By Step。本文概要性地介绍了Abaqus中心差分法的理论以及算法实现的整体流程，并通过简单弹簧显式动力学分析算例与iSolver和Abaqus计算结果进行对比，验证了算法和整体流程的正确性。

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