因你精彩
“通俗讲动力数值算法”
随着科学研究和工程技术的不断发展，出现了诸如航天飞机、空间站、海洋石油钻井平台等大型或超大型的复杂结构，它们不但自由度高，而且还含有非线性、随机载荷和复杂的边界条件等多种因素。
因此，到20世纪末开始，结构动力学发展到一个新的高度，对复杂结构的分析需要借助高阶有限元模型进行大量计算。
更由于一些细观和微观结构分析的提出，大型和超大型的计算是不可避免的。这些问题由于计算规模巨大，求解复杂，在传统的串行机上无法得到满意的解答。因此改变计算机系统结构，变串行处理为并行处理，从理论上讲，其速度的提高是无止境的。利用并行计算机研究和开发相应的并行算法有可能解决串行机无法解决的存储量、计算时间和效率等多方面的问题。
在谈计算效率及并串行之前(期待下期讲解)，这期先谈论下结构动力学的显式隐式计算。
【结构动力学研究对象】
结构动力学问题通常有两类研究对象：
第一类是在运动状态下工作的机械或结构，例如高速运行的车辆、飞行器等，它们承受着本身惯性及与周围介质或结构相互作用的动力载荷。
第二类是承受动力载荷作用的工程结构，例如建于地面的高层建筑和厂房，石化厂的反应塔和管道，它们可能承受强风、波浪、地震以及冲击等各种动力载荷的作用。
【结构动力学研究内容】
结构动力学的研究内容主要有两个方面：
在看本文之前可以看下历史文章：
在看本文之前可以看下历史文章：
对于多自由度体系，一般可采用振型叠加法，但当体系存在强非线性行为时，振型叠加法已不再适用（弱非线性行为可通过FNA，将非线性行为当成外力施加，详情可看：【JY】基于Ramberg-Osgood本构模型的双线性计算分析）。
因此目前对结构动力响应分析常用的有效方法是数值积分法。数值积分法是将振动平衡方程式中的时间分割成许多间隔，每个时间间隔都非常小以保证计算精度。针对每个时间间隔点计算位移、速度及加速度等，利用已经求得的当前步的分析结果作为已知条件，通过一定的计算方法或假定求得未知的下一步的分析结果，逐步求得结构在地震作用下的响应结果。
特别是结构的恢复力特性随结构反应的大小而在不断地变化（材料非线性），因此在每步的分析中必须根据结构反应状态确定当前的结构恢复力特性，进行下一步计算。
直接积分法针对离散时间点上的值进行计算，十分符合计算机存储的特点，体系的运动微分方程也不一定要求在全部时间上都满足，而仅要求在离散的时间点上满足即可。根据在当前时刻的反应值确定下一时刻反应值方法的不同。
因此在函数表达式的不同分别形成了不同的数值积分法，常见的数值积分法主要包括Runge-Kutta法、Newmark-β法、Wilson-θ法、Duhamel积分法、HHT-α法、中心差分法、分段解析法、精细积分法等等（看下图）。
（详情可下载学习：【JY】SDOF计算教学软件开发应用分享）
其实判断一种直接积分方法的优劣，其标准是多方面的：
a.收敛性
当计算的时间步长趋于无穷小时，所得的数值解是否收敛于精确解。
b.计算精度
如果计算截断误差是相对于时间步长N次方的小量，那么就称该方法具有N阶精度。
c.稳定性
随计算时间步数的增大，如果数值解没有远离精确解，则该方法是稳定的。
d.计算效率
计算耗时的多少可直接反映出一种计算方法的计算效率。
因此，对于不同问题下，选择不同算法是存在优劣之分的，而常见的在实际的动力弹塑性分析中还可以根据是否需要联立求解耦联方程组，且是否需要在每一步的求解中进行刚度矩阵求逆，将数值积分法分为隐式分析方法和显式分析方法。
该方法需迭代求解耦联的方程组，计算工作量大，增加的工作量至少与自由度的平方成正比，如Newmark-β法、HHT-α法，但是积分步长可以取得较大。
联立求解耦联方程组和迭代都需要求解大型的线性方程组，这一过程需要占用相当数量的计算资源、磁盘空间和内存。该算法中的增量步可以比较大，至少可以比显式算法大得多，但是实际运算中上要受到迭代次数及非线性程度的限制，需要取一个合理值，而且对于高度非线性问题无法保证收敛。
该方法直接求解耦联的方程组，计算工作量较小，工作量至少与自由度数成正比，如中心差分法，但是要求步长很小。
显式算法基于动力学方程，采用动力学方程的一些差分格式（如广泛使用的中心差分法、线性加速度法等）对时间进行差分，不需要进行刚度阵、质量阵、阻尼阵的组，不需要进行矩阵分割或矩阵求解等求解大型方程组，其右端项的形成只需根据每个单元对有效荷载向量的贡献累加而成，这样整个计算基本上在单元一级水平上进行，计算量与节点数成正比，因此只需很小的高速存储区，尤其当一系列单元的刚度阵、质量阵、阻尼阵相同时，就不需重复计算，效率更高。
另外，显式算法不需要进行平衡迭代，计算速度快，时间步长只要取得足够小（典型的显式算法的时间步长是隐式算法的1%~1%），一般不存在收敛性问题。因此需要的内存也比隐式算法要少，并且数值计算过程可以很容易进行并行计算，程序编制也相对简单。但显式算法要求质量矩阵为对角矩阵，而且只有在单元级计算尽可能少时速度优势才能发挥，因而往往采用减缩积分方法。显式计算的不足是精度不高，必须设定非常小的时间步求解以保证稳定状态，过大和过小的时间步往往导致求解时间非常漫长。特别适用于求解需要分成许多的时间增量来达到高精度的高速动力学时间，诸如冲击、碰撞和爆破等高度非线性问题。
以下面这个方程为例，讲清楚显式、隐式的差别：
为了方便阅读，提前给个结论：
a. 显式求解法：显式算法的步进具有“预测”的能力，仅需要当前步的所有参数可以预测下一步的所有参数。
b. 隐式求解法：隐式算法的步进具有自我矫正的能力，利用下一步的导数矫正当前步误差，进行实时步步矫正。
1、对于上述方程的显式方法的化解求解如下：
上述方程中可得到当|1-h|>1时(即h>2)，整个方程将开始发散， 需要保证h<2，整个方程才可收敛。求解函数绘制如下：
（即工程师们在采用显式进行计算弹塑性分析时，模型或步长不当时，计算会发生的发散同理。）
2、对于上述方程的隐式方法的化解求解如下：
可以看到隐式计算需要用到下一步的x对y进行矫正计算，且仅需要保证h>0 即可保证方程精度稳定且收敛。求解函数绘制如下：
综上所述，显式算法具有计算存储量小的优点，在计算分析步内，不需要对刚度进行重新组装迭代，但计算步长需要取得极小才可让计算不发散。
时间步长的建议可取以下建议：(其中E是弹性模量，L是单元长度，Cd是材料波速)。
隐式算法具有自动矫正的能力，且步长可以较大，但刚度矩阵计算在每一步内需要进行迭代重组装，且在迭代过程，刚度重组的过程较容易计算不收敛，并需要占用计算存储大量资源。
最近笔者也开发小型的计算内核，受曲哲老师的极速牛顿法的启发（文献DOI：https//doi.org/10.1002/nme.6456）这里提供一种思路给大家：
隐式 + PCG + 极速牛顿
解析：隐式算法使得步长拉长，PCG（预处理共轭梯度法）使得求解超大规模的稀疏矩阵变得简便快速，极速牛顿法进在迭代时最多两步，尽可能避免不收敛问题，集百家之所长的快准狠算法~
（图片来源：哲设计）
最后，其实显式和隐式求解的本质在于采用不同的物理学平衡方程，因此在不同的物理学问题也有不同的表现。顺便说下，在大型通用有限元中，如Abaqus/Standard用的隐式分析方法的为HHT-α法（改进Newmark-β法），如Abaqus/Explicit、LS-DYNA显式分析方法为中心差分法。
下期再见！
往期精彩*点击即达
#性能分析
【JY】基于Ramberg-Osgood本构模型的双线性计算分析
#软件讨论
【JY】Abaqus案例—天然橡胶隔震支座竖(轴)向力学性能
欢迎关注微信公众号: 建源学堂
