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

==概述==

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

（1）   基础理论

（2）   商软操作

（3）   自编程序

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

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

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

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

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

==第24篇：显式求解 Step By Step==

动力学问题是将力的方程和运动学方程耦合在一起的理论，在实际问题中，运动过程无时不刻的与力相关，因此数学上两个方程联立求解是最理想的。但数值计算中却难以实现，也没必要联立求解，只需先求一个方程，然后再求另一个方程，只要时间步长足够短，那么精度依然可以保证。有限元亦是如此解耦，有限元两个时刻点之间会认为力的状态不变，那么可以分两步求解。

（1）   在每个时刻点只求力学方程，得到运行学的初始条件。

（2）   在两个时刻点之间由于力的状态不变，那么可以只按运动学来求的运动学量。

本文首先研究商用有限元中最常用的显式求解算法中心差分法的理论，并给出了一个弹簧显式动力学分析的Step by Step例子，通过这个例子猜测了Abaqus中采用中心差分法求解显式动力学问题的过程，然后在自编有限元程序iSolver采用同样的算法，验证iSolver的结果和Abaqus完全一致，从而证明Abaqus的内部算法和我们设想的一致。具体验证过程也可以参考我们的演示录像：

https://www.jishulink.com/college/video/c12884  4分析篇.3-弹簧显示动力学分析

1.1 中心差分法的理论

中心差分法的标准理论可查看相应的论文，由于和Abaqus的中心差分法比较接近，所以在此不累述。

1.2 Abaqus中心差分法的理论

• 注：本节公式均摘自《Abaqus Theory Manual 2.4.5 Explicit dynamic analysis》

1.2.1 差分公式

取i-0.5时刻的速度和                                               时刻的加速度，则有下式，

        

式中，

表示速度，

表示加速度，

表示时间步大小。

带入i时刻的位移，则有，

                                  

其中，i时刻的加速度可根据牛顿定律计算得出，即，

                                                 

式中，M表示集中质量阵，F表示外力，I表示内力。

针对初始化、处理某些约束条件以及在后处理过程中，Abaqus对速度有特殊的处理，如下式，

                             

1.2.2 初始化

在初始时，即（t=0），除非用户自定义，一般情况下速度和加速度都设为0。假定有下式，

                                         

带入式（1），则有

                                         

将i=-1带入式（4）中，也可得到式（6），故前后是一致的。

初始化流程如下，

1.2.3 整体流程

以某一时刻t为例，整体流程如下：

• 其中，内力和时间步是同时求的，这个我们会在后续的VUEL文章中讲到。需要特别注意初始化时内力和时间步对应时刻点和整体流程的差别，因此我们调整了求时间步和内力的顺序以和Abaqus一致。

1.3 Abaqus的实现验证

1.3.1 模型例子：弹簧的显式动力学分析

Example 1：弹簧的显式动力学分析

（模型详见附件Job-ExplicitSpring-AddF-NoBC.inp）

创建一根线弹簧，在弹簧两侧各加两个点质量，无约束，右端X方向拉力。

参数如下：

尺寸：X方向长度L=1；

质量：各向同性，大小为100；

刚度：弹簧刚度为1000；

力：分析时间内恒定大小，为800；

时间设置：总时间为1，时间增量固定为0.2。

1.3.2 增量步零（0s）

1.3.2.1 加速度对比

1.3.2.2 速度对比

1.3.2.3 位移对比

1.3.3 增量步一（0.2s）

1.3.3.1 加速度对比

1.3.3.2 速度对比

1.3.3.3 位移对比

1.3.4 增量步二、三、四

与第一个增量步类似，不再累述。

1.3.5 增量步五（1.0s）

1.3.5.1 加速度对比

1.3.5.2 速度对比

1.3.5.3 位移对比

1.3.6 验证结果

可以发现iSolver和Abaqus完全一致。

==总结==

本文概要性地介绍了Abaqus中心差分法的理论以及算法实现的整体流程，并通过简单弹簧显式动力学分析算例与Abaqus计算结果进行对比，验证了算法和整体流程的正确性。后续文章，我们会逐步深入显式动力学的一些细节，敬请期待。

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

snowwave02 From www.jishulink.com

email: snowwave02@qq.com

以往的系列文章：

1.5.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.5.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.5.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