简介
     所谓显式和隐式，是指求解方法的不同，即数学上的出发点不一样。并不是说显式只能求动力学问题，隐式只能求静力学问题，只是求解策略不通。

     显式求解是对时间进行差分，不存在迭代和收敛问题，最小时间步取决于最小单元的尺寸。过多和过小的时间步往往导致求解时间非常漫长，但总能给出一个计算结果。解题费用非常昂贵。因此在建模划分网格时要非常注意。


  隐式求解和时间无关，采用的是牛顿迭代法（线性问题就直接求解线性代数方程组），因此存在一个迭代收敛问题，不收敛就的不到结果。


        两者求解问题所耗时间的长短理论上无法比较。实际应用中一般感觉来说显式耗时多些。由于两者解题的出发点，所以一般来说显式用于求解和时间相关的动力学问题。隐式用来求解和时间无关的静力学问题。但也不是绝对的。比如，用隐式求解时，为了克服迭代不收敛，改用显式算，但是要多给点时间，这样虽然克服了不收敛的问题，但是求解的时间费用也是相当客观的。另外，隐式也可以求解动力学问题。

        显式求解与隐式在数学上说主要是在求解的递推公式一个是用显式方程表示，一个是用影视方程来表示。比如a(n)=a(n-1)+b(n-1),后一次迭代可以由前一次直接求解，这就是显示方程，如果a(n)=a(n-1)+f[a(n)]，f[a(n))为a(n)的函数，此时a(n) 不能用方程显示表示，及数学上的隐函数，一般很难直接求解，多用迭代试算法间接求解。有限元在求解动力学问题中直接积分法中的中心差分积分就是显示求解，而Newmark积分法则为隐式积分。


隐式积分求解有限元问题
      假设现在一个物体已经被离散成有点个单元。

       显式求解与隐式在数学上说主要是在求解的递推公式一个是用显式方程表示，一个是用影视方程来表示。比如a(n)=a(n-1)+b(n-1),后一次迭代可以由前一次直接求解，这就是显示方程，如果a(n)=a(n-1)+f[a(n)]，f[a(n))为a(n)的函数，此时a(n) 不能用方程显示表示，及数学上的隐函数，一般很难直接求解，多用迭代试算法间接求解。有限元在求解动力学问题中直接积分法中的中心差分积分就是显示求解，而Newmark积分法则为隐式积分。显式求解与隐式在数学上说主要是在求解的递推公式一个是用显式方程表示，一个是用影视方程来表示。比如a(n)=a(n-1)+b(n-1),后一次迭代可以由前一次直接求解，这就是显示方程，如果a(n)=a(n-1)+f[a(n)]，f[a(n))为a(n)的函数，此时a(n)
 不能用方程显示表示，及数学上的隐函数，一般很难直接求解，多用迭代试算法间接求解。
      有限元在求解动力学问题中直接积分法中的中心差分积分就是显示求解，而Newmark积分法则为隐式积分。显式求解与隐式在数学上说主要是在求解的递推公式一个是用显式方程表示，一个是用影视方程来表示。比如a(n)=a(n-1)+b(n-1),后一次迭代可以由前一次直接求解，这就是显示方程，如果a(n)=a(n-1)+f[a(n)]，f[a(n))为a(n)的函数，此时a(n) 不能用方程显示表示，及数学上的隐函数，一般很难直接求解，多用迭代试算法间接求解。有限元在求解动力学问题中直接积分法中的中心差分积分就是显示求解，而Newmark积分法则为隐式积分。


      首先定义几个变量：
      x ： 当前物体位置；则 表示速度，表示加速度。
      为物体的质量矩阵，


对于牛顿第二定律： F=ma，用矩阵的形式表示的话:
其中：表示物体在X位置下，速度为X导的情况下的受力。


首先问题的初始条件是初始位移和初始速度。
目的是求得h时间后的位移和速度。


对于式子：




若用显示欧拉求解，记，，则可得到：


这样，就利用欧拉前向方法在小步长的条件下求得当前状态的位置和速度。


对于欧拉隐式解法，上式就变成了


前向解法和后向解法的区别就是前向的求解只依赖于前一时刻的状态量，而后向求解则不仅依赖于前一时刻的状态，还依赖于当前时刻的状态。
上式是一个非线性方程，不好直接求解，先对 f 进行泰勒级数展开：


代入原式：







用 I 表示单位矩阵，对上式进行移项，有



通过这个式子很容易求得 detaV，然后detaX也可以求出，那么当前时刻的状态就求出来了。


参考
Large Steps in Cloth Simulation

隐式与显式有限元算法


转载于 科学网 刘鹏飞的博客

链接: http://www.360doc.com/showweb/0/0/950200442.aspx.

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

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

介绍

CAE求解方法一般有两种
1. 显式(Explicit)
    第n步结果可以从n-1,n-2,...1步的结果直接推导出来，迭代时每步计算量小，但迭代增量也有个限制，不能太大，否则会出现发散
2. 隐式(Implicit)
    第n步结果不能直接从前面结果推导出来，必须做进一步的求解，这样，迭代时每步计算量大，但迭代增量可以很大
用Abaqus一段时间了，但一直对这两种方法没有一个直观的感念。正好前段时间碰到一个前向欧拉和后向欧拉的问题，顺便研究了一下。


前向欧拉和后向欧拉
前向欧拉和后向欧拉分别是显式和隐式的一个典型方法。
前向欧拉：
　　　　　　　fn+1
(x) = fn
(x) + h*f'n
(x)
后向欧拉： 
              fn+1
(x) = fn
(x) + h*f'n+1
(x)
 
 
例子
现在计算
y'(x) = -y+x+1
初值y(0) = 1;
的数值解。
前向欧拉： 
           yn+1
 = yn
 + h*(-yn
+xn
+1)
这个式子不需要做任何额外运算就能从n推导出n+1,因此每步迭代计算量小，但h存在一个最大限制.
           yn+1 
= (1-h)yn
 + h*(xn
+1)
                   = (1-h)2
yn-1
 + h*(xn
+1) + h*(xn-1
+1)
                   = (1-h)n
+O(s2
)
当h<2时，y收敛，但h>=2时，y将发散。
后向欧拉：
           yn+1
 = yn 
+ h*(-yn+1
+xn+1
+1)
这个式子不能直接得出yn+1,
必须做进一步计算得到
            yn+1
 = (yn 
+ h*(xn+1
+1)) / (1+h)
当h>0时，y无条件收敛。


结果 
h = 1.9: 

 
h = 2.1

 
Matlab代码：
 
clear all;
h = 2.1;
x_End = 30.0;
x_Start = 0.0;
Step = (x_End - x_Start) / h;
x(1) = x_Start;
y_f(1) = 1.0;%forward Euler
y_b(1) = 1.0;%backward Euler
y_theory(1) = 1.0;
for i = 2:Step
  x(i) = x(i-1) + h;
  Dy_f = -y_f(i-1)+x(i-1)+1;
  y_f(i) = y_f(i-1) + h*Dy_f;
  %Dy_b = 2*x(i);
  %y_b(i) = y_b(i-1) + h*Dy_b;
  y_b(i) = (y_b(i-1)+h*(x(i)+1)) / (1+h);
  y_theory(i) = exp(-x(i)) + x(i);
end
x_dense = 0: 0.01: x_End;
y_theory_dense = exp(-x_dense)+x_dense;
figure(1);
hold off;
plot(x, y_f, 'b', 'LineWidth', 2);
hold on;
plot(x,y_b, 'g', 'LineWidth',2);
plot(x_dense, y_theory_dense, 'r', 'LineWidth', 2);
legend('forward Euler', 'backward Euler', 'theory');
Tolerance_f = abs(y_theory - y_f);
Tolerance_b = abs(y_theory - y_b);
figure(2);
hold off;
plot (x, Tolerance_f, 'b', 'LineWidth', 2);
hold on;
plot(x, Tolerance_b, 'g', 'LineWidth', 2);
legend('forward Euler tolerance', 'backward Euler tolerance');

 

迭代(iteration)指的是去获取元素的一种方式，一个接一个。当你显式或隐式的使用循环来遍历某个元素集的时候，那就是迭代。

在Python里面，可迭代对象(iterable)和迭代器(iterator)有着特殊的含义。


	
iterable是实现了__iter__()方法的对象，该方法会返回一个iterator对象
	
	

	
iterator是实现了__iter__()和__next__()方法的对象，__iter__()方法返回的是iterator对象本身
	
由此可见，iterable和iterator的本质区别就是后者多了一个__next__()方法。也就是说一个iterator对象必定是一个iterable对象。

当你使用一个for循环或者map，或着一个列表推导，那么会先通过iter()获取相应的迭代器， 然后每次循环自动通过next方法调用这个迭代器(iterator)，从中获取每一个元素，从而完成迭代过程。

在一个iterable对象上执行iter会返回一个iterator对象, 比如iter(obj)

下面一个例子可以非常清晰的解释清楚：


			
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

			
			

			
>>> s = 'cat'      # s is an ITERABLE

>>> t = iter(s)    # t is an ITERATOR
                   # t has state (it starts by pointing at the "c")
                   # t has a next() method and an __iter__() method
>>> next(t)        # the next() function returns the next value and advances the state
'c'
>>> next(t)        # the next() function returns the next value and advances
'a'
>>> next(t)        # the next() function returns the next value and advances
't'
>>> next(t)        # next() raises StopIteration to signal that iteration is complete
Traceback (most recent call last):
...
StopIteration

>>> iter(t) is t   # the iterator is self-iterable

			
		
 

可以使用isinstance()判断一个对象是否是Iterable对象：


			
1
2
3
4
5
6
7
8
9
10
11

			
			

			
>>> from collections import Iterable
>>> isinstance([], Iterable)
True
>>> isinstance({}, Iterable)
True
>>> isinstance('abc', Iterable)
True
>>> isinstance((x for x in range(10)), Iterable)
True
>>> isinstance(100, Iterable)
False

			
		
 

可以使用isinstance()判断一个对象是否是Iterator对象：


			
1
2
3
4
5
6
7
8
9

			
			

			
>>> from collections import Iterator
>>> isinstance((x for x in range(10)), Iterator)
True
>>> isinstance([], Iterator)
False
>>> isinstance({}, Iterator)
False
>>> isinstance('abc', Iterator)
False

			
		
 

生成器都是Iterator对象，但list、dict、str虽然是Iterable，却不是Iterator。

把list、dict、str等Iterable变成Iterator可以使用iter()函数：


			
1
2
3
4

			
			

			
>>> isinstance(iter([]), Iterator)
True
>>> isinstance(iter('abc'), Iterator)
True

			
		
 

你可能会问，为什么list、dict、str等数据类型不是Iterator？

这是因为Python的Iterator对象表示的是一个数据流，Iterator对象可以被next()函数调用并不断返回下一个数据，直到没有数据时抛出StopIteration错误。可以把这个数据流看做是一个有序序列，但我们却不能提前知道序列的长度，只能不断通过next()函数实现按需计算下一个数据，所以Iterator的计算是惰性的，只有在需要返回下一个数据时它才会计算。

Iterator甚至可以表示一个无限大的数据流，例如全体自然数。而使用list是永远不可能存储全体自然数的。

小结

凡是可作用于for循环的对象都是Iterable类型；

凡是可作用于next()函数的对象都是Iterator类型，它们表示一个惰性计算的序列；

集合数据类型如list、dict、str等是Iterable但不是Iterator，不过可以通过iter()函数获得一个Iterator对象。

Python的for循环本质上就是通过不断调用next()函数实现的，例如：


			
1
2

			
			

			
for x in [1, 2, 3, 4, 5]:
    pass

			
		
 

实际上完全等价于：


			
1
2
3
4
5
6
7
8
9
10

			
			

			
# 首先获得Iterator对象:
it = iter([1, 2, 3, 4, 5])
# 循环:
while True:
    try:
        # 获得下一个值:
        x = next(it)
    except StopIteration:
        # 遇到StopIteration就退出循环
        break