谈一谈有限元分析中的显式和隐式解法

   假设时间步i中的所有物理量Xi已知, 如果时间步i+1中的所有物理量Xi+1可以表为Xi的显函数
(1)           Xi+1 = F(Xi)  (注;本文以下标i,i+1等表示时间步)
则该物理量X可以显式计算。而如果
(2)          Xi+1 = F(Xi+1)
则需要隐式解法来计算X。
    上述定义应该没什么歧义,但其应用于有限元时模样大变,似乎不那么好辨认(参考下述文献1,2).本文试图理清一下这一概念。

1.  静力学计算
    以K表示刚度矩阵,U表示位移,F表示外力。时间步i+1的外力已知。则
        KiUi+1=Fi+1     =>   Ui+1 = Ki^-1  Fi+1
为静力学显式方程。由于采用了时间步i的刚度矩阵,在用于计算非线性问题时该方程并不能满足时间步i+1的平衡条件,其使用条件是步长要足够小。而方程
        Ki+1Ui+1=Fi+1
中时间步i+1的刚度矩阵Ki+1为Ui+1的函数(非线性问题),为静力学隐式方程。一般情况下我们用Newton-Raphson法求解此方程。

2.  动力学计算
     在动力学计算中,除了需要考虑刚度矩阵K的线性,非线性,还需要考虑时间的差分格式,时间的差分可以是显式的,也可以是隐式的。
     显式格式; 如Euler前进差分  Ui+1 = Ui + Vi* h + 0.5Ai* h^2+ ... (在这里,我们用h表示时间增分,V,A表示速度和加速度)。在有限元动力学计算中我们常用的是中心差分
         Ai= (Ui-1-2Ui+Ui+1)/h^2
         Vi= (-Ui-1+Ui+1)/(2h)
把i+1时间步的位移Ui+1表为时间步i的速度和加速度的函数。
     隐式格式; 隐式差分格式很多,Euler后退差分属于隐式格式差分。有限元计算中常用Newmark法
         Vi+1= Vi + [(1-a)*Ai +a*Ai+1]h
         Ui+1= Ui + Vi* h + [ (0.5-b)*Ai+ b*Ai+1]*h^2   (a,b为常数)
在这里Ui+1, Vi+1和Ai+1互相相关,不是既知时间步i的物理量的显式函数。
     如果只用时间差分格式来区分,那么我们可以定义采用显式时间差分格式的方法为动力学显式计算,而采用隐式时间差分格式的方法为动力学隐式计算。但是如果我们考虑非线性,情况就有些复杂!比如说下面这个问题:

3.  如何用显式时间差分格式来计算非线性问题?
     如果我们仍用Newton-Raphson法来处理非线性,也就是说时间步i+1中的物理量Xi+1不是Xi的显函数,此时如果将其称为显式动力学计算不是与上面的显式计算的定义(1)相矛盾吗?对这个问题的回答是:
     在显式动力学计算中,我们无视每一计算步长中由于各种非线性带来的非平衡残差,也就是说我们视每一特定计算步中的变形为线性。当然,这种处理的基本前提是,在显式动力学中每一计算步的时间增分足够小。
     也就是说我们在显式动力学计算中不会使用Newton-Raphson法. 至少在作者的认知范围内情况如此。如果存在反例还请告知。
     显式动力学计算的平衡方程为
       Mi Ai + Ci Vi + Ki Ui  = Fi
将时间的中心差分带入可得到
        [ Mi /h^2 + Ci /(2h) ] Ui+1 = Fi -(Ki -2Mi /h^2)Ui - [ Mi /h^2 -Ci /(2h)] Ui-1
可以看到刚度矩阵K出现在方程右端,这一现象带来了该方程的最迷人之处: 如果矩阵M,C为对角阵(比如通过对采用对角化的lumped mass matrix并取衰减阵C比例于M),我们可以不用线性方程组求解器就可以求解这一方程。我们甚至不需要计算刚度矩阵K(因为方程右端的KU代表变形产生的内力,我们只需计算变形产生的内力即可)。所以其计算量,所需内存及并行计算的通信量都大大减少。
     隐式动力学计算,考虑的是时间i+1的平衡方程
         Mi+1 Ai+1 + Ci+1 Vi+1 + Ki+1 Ui+1  = Fi+1
因此采用这种方法计算非线性问题时仍需使用Newton-Raphson法。但是

4.  采用隐式动力学计算格式计算线性问题时,其计算方法属于上面的显式计算(1)还是隐式计算(2)?
     回答是属于显式计算(1). 如此,力学算中的式指的是时间差分格式。也有人(·如参考文献(3))认为考虑时间i的平衡条件的是显, 考虑时间i+1的平衡条件的是隐计算。相对于本文给出的基于显,隐的判别法,这也是不失为一个有效的分类方法。

5.  静力学显式计算和动力学显式计算的比较
      采用动力学显式计算时我们可以不用求解线性方程组,但是静力学显式计算是需要的。因此使用静力学显式计算并没有什么优势。实际上很少见到采用这种方法的报道,有极少的例外如参考文献(4,5),算得上是珍稀品种了。
      有效的静力学显式计算方法应该是采用所谓准静力学显式(动态松弛法)计算,导入假想的质量矩阵M, 使用显式动力学计算公式来进行静力学计算。

结论: 由于存在对于非线性处理和对时间差分的显式和隐式计算方法,给对有限元分析中的显式和隐式解法的理解带来了一些困难。希望本文能澄清这一问题。
   有限元分析中实际存在下述算法:
* 静力学隐式计算: 标配;  所谓静力学显式算法属于珍稀品种。
* 动力学显式计算: 标配;  显式时间差分格式。在计算时间步内不考虑非线性。
* 动力学隐式计算: 标配;  隐式时间差分格式。在计算时间步内考虑非线性,作迭代运算。

1. What is difference between implicit and explicit FEM in non-linear quasi-static problems?
https://www.researchgate.net/post/What_is_difference_between_implicit_and_explicit_FEM_in_non-linear_quasi-static_problems
2.  Implicit and Explicti finite element (http://imechanica.org/node/5396)
3. K.J.Bathe,E.L. Wilson: Numerical Methods in Finite Element Analysis, Prentice-Hall·, 1977
4.Staic-explicit finite element method and its application to drawbead process with spring-back
http://www.sciencedirect.com/science/article/pii/S0924013602004685
5.  http://trialpark.co.jp/en/

VS Code下cmake, c++编译,调试环境的构成步骤

1   下载必须extension      按[Ctrl+Shift+X]打开extension窗口,选择安装"C/C++", "CMake", "CMake Tools" 2   在VSCode下打开作业目录 ...