1 概述

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

（1）基础理论

（2）商软操作

（3）自编程序

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

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

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

通用结构有限元软件iSolver介绍视频：

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

==第32篇：线性瞬态动力学==

在系列文章第三十篇、三十二篇时，我们介绍了稳态动力学的原理和算法，本章继续介绍瞬态动力学。稳态和瞬态都是在随时间变化的激励作用下产生的运动形式，其中，稳态是时间足够长后，系统慢慢趋于稳定后的运动，譬如当不考虑阻尼时弹簧受瞬间载荷力作用时间足够长后最后会做简谐振动，由于简谐振动的理论公式比较简单，一个简谐振动在时间域上虽然很长，但只要有频率和振幅参数就能完全表示，所以可以跳过前面的不稳定的过程而直接求解最终的稳态形式。但达到稳态前的瞬态过程中间每个时刻点的运动状态是不一样的，无法简单的用几个参数表示，当前时刻的运动状态需要前面运动状态往前推进得到，此时只能用瞬态分析。

瞬态分析有线性和非线性之分，从线性瞬态动力学出发其实更容易理解瞬态动力学的求解方式。本章我们将简单介绍一下瞬态动力学的求解公式，并以一个单摆例子来说明线性瞬态动力学在有限元软件中的内部实现原理。非线性瞬态动力学的原理将在后面章节介绍。

1.1 线性瞬态动力学的理论

1.1.1 线性瞬态动力学的方程

线性瞬态动力学的运动方程和线性静力平衡方程很类似，对于线性静力问题载荷F和位移u是直线关系，可以非常简单的由一个状态类乘以一个系数推到另一个状态。因为静力学可以简单的表示为下方的线性表达形式：

Ku=R

譬如下面的悬臂梁问题，载荷1N的时候计算得到最大位移时1.177mm。

如果载荷是1000N时，不用计算也知道就是最大位移会变为1000倍=1177mm。

但在线性动力学问题中分析复杂一点就是，如果力扩大了1000倍，显然最终的位移不是简单的乘上1000，瞬态动力学和上式原理上很类似，只不过包括了质量阵M相关的惯性力项和阻尼阵C相关的阻尼力项，如下：

一般质量矩阵与时间和位移无关，而K和C如果和位移和时间无关，那么就是线性瞬态问题。为简单起见，我们不考虑阻尼项。得到如下表达式：

位移和加速度与t相关，上述方程就有两个未知数位移和加速度，只有一个方程是无法求出的，所以还有一个加速度和位移的约束关系，实际中必然是加速度先和速度相关，速度再和位移相关。

1.1.2 线性瞬态动力学的有限元求解

在理论上当前时刻的加速度就应该只和当前时刻的位移相关了，但在计算机的数值求解时，无法求出无限个连续时刻的结果，所以只能把整个时间人为的分成多个时间段，方程变为如下：

此时，加速度就必须要通过起码两个时刻的速度或者位移才能计算出来的，可以有不同的取法，取前时刻还是后时刻的位移来决定加速度就有隐式和显式之分（具体可看系列文章13：显式和隐式的区别），任何一种分析理论上都可以得到正确解，只不过显式不用求刚度，隐式需要求解刚度阵，而通过刚度的解释可以更容易理解Abaqus内部的实现原理，在本章中，我们选用隐式瞬态分析。

在隐式方法中，将时间离散，加速度表示为位移的表达式后，得到下方的公式

其中K和F在静力学的增量迭代法中表示切线刚度阵和内力，只不过这边需要加上M项：

其中A是只与时间增量步相关的量，对不同的隐式算法A的值不同，譬如对最简单的Newmark方法中的梯形隐式算法，加速度在增量步内线性变化，下方红线是在输入的载荷转换为加速度的值，但计算中没法处理连续的曲线，所以Abaqus或者iSolver等有限元程序实际上只会认为增量步

内加速度时线性变化的：

此时得到的修正刚度为：

显然当时间增量步固定时，此时的K也是个定值。而F与上一时刻的包括惯性力在内的内力，与位移、速度、加速度相关，可以表示为：

其中B是只与时间增量步相关的量。虽然固定步长K是常量，但右端量还包括了上一个时刻位移、速度等变量，R和u依然为非线性曲线，如下：

由上面的表达式，K和F只与增量步I,j相关，但实际上只要一次迭代左右就能相等平衡了，这个原理就和线性弹性力学的内力F=KU，到时右端在第一个迭代完成后就和外力R一样了，具体可看下面的课程推导：

https://www.jishulink.com/college/video/c14948 深入浅出有限元：基础理论->Abaqus操作->matlab编程

也就是说线性瞬态动力学方程与迭代步j无关，无须迭代，方程可以简单的写成：

而且只要迭代一次，而且内力项在第一次是位移为0，那么右端F就可以简单的写成下方表达式：

1.2 线性瞬态动力学的单摆模型

下面我们将以一个简单的单摆模型例子来说明上述线性瞬态动力学的解法。

1.2.1 实际单摆模型

实际的单摆如下图，线长L，无质量，末端绑定一个质量块m，如果没有空气阻力，在重力加速度g作用下往复运动。我们取初始时刻绳索在水平位置。假定绳索足够硬，在小球作用下变形非常小。显然，这个问题是一个大位移大转动小应变的几何非线性问题（具体可看我们的系列文章18：几何非线性的应变），这个小球位置随时间变化的运动过程就是非线性的瞬态过程。