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

3年前浏览2330

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

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

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

首先，对于本构方程的求解，通常分为前向和后向欧拉算法。对于后向欧拉算法求解弹塑性问题，所有的物理量（包括等效塑性应变增量、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方法，将刚度矩阵和阻尼矩阵分成两个部分，在同一区域采取不同算法，提高计算效率和精度、稳定性和收敛性。

理论科普求解技术结构基础 Abaqus

著作权归作者所有，欢迎分享，未经许可，不得转载

首次发布时间：2021-04-05