6-有限元之显式算法和隐式算法

流苏kiwi

2020年3月11日 12:17

浏览：6162 评论：13 收藏：45

最近在学有限元，看到不少资料都提到了 EXPLICIT ALGORITHM 和IMPLICIT ALGORITHM，很是费解，查询到一篇文章写的不错，特在此与大家分享，文章出处：公众号 CAE&CFD

显式（explicit）和隐式（implicit）这两个词在有限元分析中大家可能经常看到，特别是涉及到动力学分析时。但其实广义的说他们分别对应着两种不同的算法：显式算法（explicit method）和隐式算法（implicit method），所以不论在动力学或者静力学中都有涉及到。

我们都知道有限元分析FEA在计算微分方程（differential equations）时，由于计算本身的局限，比如计算机储存的位数有限，以及方程本身的复杂性，计算机运用的是数值算法（numerical algorithm）来逼近真实解的。有限元分析中数值算法的基础是欧拉法（Euler method），欧拉法又分为forward Euler method 和backward Euler method，这两种方法被简称为显式法（explicit method）和隐式法（implicit method）。

首先我们来看看这两种算法的区别。

显式算法（explicit method）（forward Euler method）
考虑常微分方程：

初始条件：

设h为每一步的时间步长，在tn时刻

(n=0,1,2,3...)，在T(n+1)时刻有：

所以在显式算法中，T(n+1)时刻的值由T(n)时刻决定，也就是说当前时刻的值由上一时刻的值决定。

隐式算法（implicit method）（backward Euler method）

考虑同一个方程，在T（n+1）时刻有：

所以在隐示算法中，T(n+1)时刻的值不光由T(n)时刻决定，还由当前时刻T(n+1)决定。也就是说当前时刻的值由上一时刻和当前时刻的值共同决定。隐式算法往往需要求解二次方程。

我们来看看一个具体事例：

设常微分方程：

根据上面的方法，对于显示算法有：

得出：

对于隐式算法有：

导出二次方程：

求解得:

所以很明显，在隐式算法中，要求得K+1时刻的值，就需要求解二次方程的根。

关于收敛性

显式算法不存在收敛性的问题（因为不进行收敛计算），从方程中可以看出来，每个时刻的值由上一时刻所确定，所以一步一步进行下去，当时间步取得较大时，就会偏离真实值。

显式算法的过程（蓝色为真实值）

隐式算法是无条件收敛的，在隐式算法中，在求解二次方程的同时，会通过Newton–Raphson method算法对每一步进行迭代收敛，直至收敛到指定的偏差。如下图所示：

时间步长（time integration）的依赖性（时间变量只在动力学中涉及）

显式算法要获得准确的结果，需要取很小的时间步长

隐式算法对时间步长要求不高，由于是绝对收敛的，往往可以取较大的时间步长。

运用上面的方法，我们以方程为例，通过数值算法求得f(u)。当把时间步长取为1时，显式（explicit）和隐式（implicit）的结果如下图所示：

可以看出，隐式算法是绝对收敛的，每一步都没有偏离真实值，而由于时间步长取得很长，所以显式算法的结果远远偏离了真实值。

当把时间步缩小到0.05时，显示算法的结果如下图所示：

可以看出，当把时间步取得很小时，显示算法可以很接近真实值。

上面主要讲了隐式和显式算法的差别，下面我们来看看这两种方法在动力学分析中的运用和差别。

动力学分析（Dynamics Analysis）

静力学（static）分析不考虑质量/阻尼和时间，而动力学分析需要考虑系统阻尼和时间的变化。

首先大家要知道有限元分析FEA的输出是什么，虽然我们可以从仿真后处理中得到很多的结果，如应力，应变，位移等等，但本质上，所有的物理量都是通过先计算出节点处的位移，然后导出应变，再通过应变根据材料力学的理论导出其它物理量的，这一点大家要记住。

在有限元分析中，动力学分析的基本方程是由如下方程导出和决定的：

[M]{a} + [C]{v} + [K]{x} = {F}

其中[M]是质量矩阵，[C]为阻尼矩阵，[K]为刚度矩阵，a为加速度，v为速度，x为位移，{F}表示外力。如果我们把方程写为导数的形式，则有：

所以这里的加速度，速度和位移是彼此关联的，这很有用。

隐式动力学（explicit dynamic）

在动力学分析中，隐式分析直接计算位移x，而要求得位移x，就需要对刚度矩阵K进行求逆（inversion），而计算机在进行矩阵求逆时，需要耗费大量时间和计算机内存。那么我们需要对方程做如下变动：

而

，所以有：

所以对刚度矩阵K求逆是必须的。

显式动力学（explicit dynamic）

在动力学分析中，显式分析先计算加速度a，再通过积分算法分别导出速度v和位移x。比如一旦在时刻n求得了加速度，速度会在n+1/2时刻计算出，然后位移在n+1时刻求出（相当于一个时间步长），然后通过位移导出应变（strain），在通过应变导出应力（stress）。

可能有人会问，和隐式动力学类似，要求得加速度a，不是一样要对质量矩阵M求逆吗，是的，但这里M的导出远远比刚度矩阵简单很多，而且M可以简化为对角矩阵，所以对M不需要直接求逆，通过简单的矩阵乘法就可以得出逆矩阵。

计算的效率

在动力学分析中，显式法和隐式法在计算时间上各有优缺点

显示动力学由于是间接求得位移x的，和前面类似，要取得足够的精度，需要取很短的时间步长，所以需要进行很多步计算，但每一步计算需要的时间很短

隐式动力学由于是直接求得位移x，所以不存在收敛性问题，但由于动力学分析往往涉及到非线性（几何非线性，材料非线性等），可能每一步都要进行刚度的求逆，非常耗时，但时间步长相对显示算法可以取得很大。

适用范围

虽然在计算效率上各有优势，但总的来说，在动力学分析中，主要还是采用显式算法，因为在动力学分析中，如频率响应分析，响应谱分析等等，往往时间很短，要获得较好的结果，需要取很短的时间步长来捕捉瞬时的响应，这时候显示动力学就很有用，而且在计算时间上也具有较大优势。

隐式算法主要还是用于时间周期较长的非线性分析（Nonlinear analysis），属于静力学分析的范畴，当然也包括动力学分析里面特殊的准静力分析（quasi-static）。

静力学分析中是不会用到显式算法的，动力学分析中主要采用显式算法（特别是响应时间很短的问题），当然也可以采用隐式算法，但需要选择较为合适的时间步长，当时间步长取得合适时，抛开计算时间，两种算法的结果是相差不大的，同一个问题大家可以两种方法都试一试，但总的来说，动力学分析大多数问题还是采用显示算法的。

显式算法的例子：