有限元学习必须知道：隐式与显式有限元算法

weixin

　　导语：对于有限元的学习，最最关键的其实不是本构方程，也不是屈服准则，而在于对求解算法的理解。本人根据长期学习经验，在美国作访问学者学习经历，纯属兴趣，在这里做一点肤浅的总结，完全原创。

　　隐式与显式有限元最大的区别在于是否迭代，是否所有的物理量在同一时刻获得。采用隐式迭代求解平衡方程(位移、速度和加速度)、而不管是否用隐式与显式的方法(前向或者后向欧拉求解方法)求解本构方程(应力和应变)叫做隐式有限元;用显式时间积分的方法求解叫做显式有限元。

　　首先，对于本构方程的求解，通常分为前向和后向欧拉算法。对于后向欧拉算法求解弹塑性问题，所有的物理量(包括等效塑性应变增量、N+1迭代步的应变和应力以及相关依赖于solution的状态变量)均是同时求解获得，因为涉及到多个物理量，而通常情况下他们是相互依赖、相互成为函数，所以必须通过牛顿迭代同时求解几个方程(如采用试应力方程、屈服函数径向返回算法(对于各向异性，也叫回映算法，最近点的投射算法)联合求解等效塑性应变增量)。对于前向欧拉，直接由N时刻的应力和应变求出N+1时刻的应力和应变，无需迭代。

　　其次，对于平衡方程的求解，通常分为隐式和显式有限元算法。对于隐式有限元算法，由应力平衡方程+边界条件变分之后获得的刚度方程KU=F，隐式求解必须引入雅可比矩阵(二次收敛、只影响计算速率、不影响数值精度;K又称为雅可比)，其是实时更新的，是N+1时刻的应力、应变以及状态变量(如损伤内变量)的函数，隐式求解是很robust的，确保了计算精度，但是不足之处在于计算非常expensive，每次迭代都要计算K的逆矩阵，也容易产生数值收敛性问题，目前解决的方法有弧长法、粘性阻尼法等，个人认为粘性阻尼法效果最好。

　　对于显示算法，采用时间积分，用t+1时刻的积分点应力、应变，获得t+1时刻的节点位移，无需迭代求解，也不需要雅可比矩阵(应力对应变偏导数);如果硬是要有，连续雅可比，基于本构模型而不是刚度方程推导近似的连续雅可比。对于显示算法，单元高斯积分点应力、应变的求解可用前向或者后向欧拉方法，然后通过时间积分求取节点位移。本质上，平衡方程中位移的迭代求解与本构方程中的应力、应变求解没有关联，这点很容易造成误解，很多时候将前、后欧拉算法视为显式和隐式的区别，大大错误。通常应用较广的显示算法纽马克法、威尔逊-sita法，其中改变纽马克法中的两个参数，可以实现隐式与显式求解，其中alpha=0.5和beta=0是中心差分法(二阶精度)。目前一个大的误区认为只有显示算法可以求解动力学问题，隐式只能求解准静态问题(如低速冲击)，alpha=0.5和beta=0.25就是隐式，所有的物理量在t+1时刻同时求解，通常ABAQUS软件中所说的隐式动力学求解采用了斯坦福大学Hilber、HUGHES院士(现在德克萨斯大学奥斯丁分校)和加州大学伯克利分校Taylor院士提出的无条件稳定隐式差分算法，可以求解低速动力学问题，缺点是不适合含阻尼的求解、计算效率不高;alpha=0.5和beta=0时的纽马克法更适合求解动力学问题，主要原因在于比隐式求解计算效率更高，不足之处在于其是条件稳定，时间增量过大位移解容易震荡，根本原因是差分算法的条件稳定导致的，时间增量必须非常小(其值越大，一方面不稳定、另一方面计算误差也更大)，其依赖于波速、弹性模量和最小单元网格尺寸，这是显式算法计算最耗时的地方。相对于隐式算法，显式算法的功能更强大，通常计算依赖于率的变形和应力，也可以求解稳态问题，如alpha=0.5和beta=0时，对于刚度方程中引入阻尼矩阵后，叫做动态松弛法，可以解决静力学问题的一些稳态问题(如重力、预应力引起的初始应力)。此外，一些准静态的剪切自锁问题，本质上有解，但是用牛顿法求解失效，中心差分引入质量矩阵后，可以获得正常的解。需要注意的是，时间积分算法通常采用Lumped集中对角质量矩阵而不是一致质量矩阵，以提高计算效率。总体来说，由于计算效率的问题，隐式时间积分算法ABAQUS-Standard特别适合于低速冲击问题;对于高速冲击问题，由于存在不连续非线性接触的动响应过程，隐式算法解决不好，使用显式时间算法ABAQUS-Explicit更好。此外，对于瞬态和稳态热传导问题，半离散的抛物线方程，中心差分法可以较好获得温度分布。

　　对于依赖于率的粘塑性问题(对于本质上的粘性材料)，与弹塑性材料的根本区别在于，一般来说是一致性条件不满足(排除弹塑性材料在高温下的软化问题，对于这种问题，屈服条件也可以满足)，即屈服条件不满足，N+1时刻的物理量不用回映到N+1时刻的屈服面上，粘塑性模型成为过应力模型，显示和隐式算法都可以求解。对于依赖于率的本构模型，其可解决模拟高速冲击、爆炸、弹道射击问题时存在的动态应变局部化问题(对于动态问题，平衡方程丧失双曲线特性;对于静态问题，平衡方程失去椭圆性)，解决网格尺寸效应，其实质上是引入了适当的阻尼迟滞效应。需要注意的是，对于大变形(又称为有限变形)问题，Cauchy应力率和速度梯度(包括客观和对称的扭曲张量率D、不客观和反对称的spin旋转张量W两个部分)均是不客观的，为解释刚体旋转(如纯剪切变形就包含刚体旋转)，在共旋坐标系下面求解真实应力和应变，应力和应变积分求解的时候应首先求解客观性的Jaumann应力率(相对于真实应力，空间坐标系)，相对于Second Piola-Kirchhoff应力是Truesdell率(材料坐标系)。ABAQUS软件对于大变形问题已经做了旋转。

　　对于一些耦合场问题，由于计算量非常大，同时要求解太多物理量，如热-流-固耦合，要求解位移、压力、温度，采用纯隐式算法或显式算法基本不太获得收敛或准确的解，这时候可采用混合的隐式与显式有限元格式mixed implicit-explicit partitoning方法，将刚度矩阵和阻尼矩阵分成两个部分，在同一区域采取不同算法，提高计算效率和精度、稳定性和收敛性。

　　来源：科学网 刘鹏飞的博客