显式求解方法和隐式求解方法对比

CUKUR。在时域内求解该方程最常用的方法是直接积分法，而程组MUt又根据求解过程中是否需要迭代求解线性方程组，将直接积分法分为隐式积分方法和显式积分方法两类。隐式积分法认为ｔ＋Δｔ时刻系统的状态不仅与ｔ时刻状态有关，且与ｔ＋Δｔ时刻某些量有关。因此隐式算法是根据tn及tn-1...时刻体系的物理量值建立关于以tn+1时刻物理量为未知量的线性方程组，通过求解方程组确定tn+1时刻的物理量（常用的方法有线性加速度法、常平均加速度法、Newmark方法、Wilson-θ法、Houbolt方法等）。而显式积分法认为ｔ＋Δｔ时刻系统的状态仅与ｔ时刻状态有关可，因此可由tn及tn-1...时刻体系的物理量值直接外推tn+1时刻物理量值（如中心差分法），不需要求解线性方程组，实现了时间离散的解耦。

两种算法的比较 ： （1）隐式算法

隐式算法基于虚功原理，要迭代计算。隐式算法在每一增量步内都需要对静态平衡方程进行迭代求解，并且每次迭代都需要求解大型的线性方程组，这一过程需要占用相当数量的计算资源、磁盘空间和内存。理论上在这个算法中的增量步可以很大，但是实际运算中上要受到接触以及摩擦等条件的限制。随着单元数目的增加，计算时间几乎呈平方次增加。由于需要矩阵求逆以及精确积分，对内存要求很高。隐式算法的不利方面就是收敛问题不容易解决，且在开始起皱失稳时，在分叉点处刚度矩阵出现奇异。 （2）显式算法

显示算法基于动力学方程，无需迭代，包括动态显式和静态显式算法。动态显式算法采用动力学方程的中心差分格式，不用直接求解切线刚度，不需要进行平衡迭代，计算速度快，也不存在收敛控制问题。该算法需要的内存也比隐式算法要少。数值计算过程可以很容易地进行并行计算，程序编制也相对简单。它也有一些不利方面。显式算法要求质量矩阵为对角矩阵，而且只有在单元级计算尽可能少时速度优势才能发挥, 因而往往采用减缩积分方法，容易激发沙漏模式，影响应力和应变的计算精度。

静态显式法基于率形式的平衡方程组与Euler前插公式，不需要迭代求解。由于平衡方程式仅在率形式上得到满足，所以得出的结果会慢慢偏离正确值。为了减少相关误差，必须每步使用很小的增量，通常一个仿真过程需要多达几千步。由于不需要迭代，所以这种方法稳定性好，但效率低。

以下是ＡＢＡＱＵＳ软件中显式和隐式常采用的方法：（摘自帮助文档）

CUKUR （1） 结构系统的通用运动学方程为：MUt中心差分法（显式）

假定0，t1，t2，…，tn时刻的节点位移，速度与加速度均为已知，现求解tn(tt)时刻的 结构响应。中心差分法对加速度，速度的导数采用中心差分代替，即为：

1(U2UU) Uttttttt2t1(UttUtt) (2) U2t将（2）式代入（1）式后整理得到 ˆURˆ （3） Mttt式（3）中

ˆ1M1C Mt22tˆtRt(K2M)Ut(1M1C)Utt Rt2t22t分别称为有效质量矩阵，有效载荷矢量。R，M，C，K为结构载荷，质量，阻尼，刚度矩阵。

求解线性方程组（3），即可获得tt时刻的节点位移向量Utt，将Utt代回几何方程与物理方程，可得tt时刻的单元应力和应变。 中心差分法在求解tt瞬时的位移Utt时，只需tt时刻以前的状态变量

ˆ，有效载荷矢量Rˆ，即可求出Utt，故Ut和Utt，然后计算出有效质量矩阵M

t称此解法为显式算法。

中心差分法，在开始计算时，需要仔细处理。t=0时，要计算Ut，需要知道Ut的值。因此应该有一个起始技术，因而该算法不是自动起步的。由于U0，，U是已知的，由t=0时的（2）式可知： U00tU UtU0tU002中心差分法中时间步长t的选择涉及两个方面的约束：数值算法的稳定性和计算时间。中心差分法的实质是用差分代替微分，并且对位移和加速度的导数采用线性外插，这限制了t的取值不可过大，否则结果可能失真过大。

可以证明：中心差分法是条件稳定的。即当时间步长t必须小于由该问题求解方程性质所决定的一个时间步长的临界值。LS-DYNA中，采用“变时间步长法”，即每一时刻的步长t由当前结构的稳定性条件来控制。具体算法为：计算每一个单元的极限时间步长tei，i=1,2…，取tmin(tei)为下一个时刻的时间步长。

2各种单元的tei计算方法如下。

（1）1D杆，梁单元

Lte

C其中为时间步长因子，系统默认为0.9。L为杆，梁单元的长度。CE材料声速。

（2）2D板，壳单元

Ltemin

CLmin为 壳单元的最小单元边长度。其中为时间步长因子，系统默认为0.9。

CE(1)2为材料的声速。

（3）3D单元

Lete 221/2Q(QC)其中

kk(kk0)C1CC0LeQ

kk0)0(C0和C1为无量纲常数，默认C0=1.5，C1=0.06.

Ve(对8节点体单元）A Le为单元等效长度， Ve为单元体积，LeemaxL（对4节点体单元）minAemax为单元最大侧面积。

E（1-）为材料声速。

(1（)1-2）时间步因子可由用户设置，减小相当于减少时间步长。设置时间步长因子的关键字为*CONTROL_TIMESTEP，控制参数为TSSFAC。

另外，质量缩放可以人为控制时间步长。即调整单元密度来改变时间步长。以壳单元为例说明质量缩放改变时间步长。

Ltemin （假如=1）

CE C2(1)C由(tspecifiedLmin(tspecified)2E(12)i得到i2 )2Lmin(1)E2

Newmark法（隐式）

Newmark假定在时间间隔[t,tt]内，加速度线性变化，即采用如下的加速度，速度公式：

UttUt[(1)UtUtt]t

t[(1)UU]t2 (4) UttUtUtttt2式中，为按积分的精度和稳定性要求可以调整的参数。

和U用U，U，U表示的表达式，代入（1）根据（4）式可给出Utttttttt式中整理得到

ˆURˆ （5） Ktttt其中

ˆ1MCK Kt2t(11)U]C[U(1)U(1)tU]ˆttRttM[1Ut1URttttt2tt2t2 称之为有效刚度矩阵和有效载荷矢量。由上式可以看出求解当前Utt，需要用到当前时刻的Rtt，因此该算法为隐式算法。

当载荷历史全部已知时，Ftt为已知量，求解需要迭代实现。

可以证明，当参数0.5，0.25(0.5)2时，Newmark法是无条件稳定的，即t的大小不影响数值稳定性。此时时间步长t的选择主要根据解得精度确定。一般，Newmark法可以比中心差分法的时间步长大得多。

比较两种算法，显式中心差分法非常适合研究波的传播问题，如碰撞、高速冲击、爆炸等。显式中心差分法的M与C矩阵是对角阵，如给定某些有限元节点以初始扰动，在经过一个时间步长后，和它相关的节点进入运动，即U中这些节点对应的分量成为非零量，此特点正好和波的传播特点相一致。另一方面，研究波传播的过程需要微小的时间步长，这也正是中心差分法的特点。

而Newmark法更加适合于计算低频占主导的动力问题，从计算精度考虑，允许采用较大的时间步长以节省计算时间，同时较大的时间步长还可以过滤掉高阶不精确特征值对系统响应的影响。隐式方法要转置刚度矩阵，增量迭代，通过一系列线性逼近（Newton-Raphson）来求解。正因为隐式算法要对刚度矩阵求逆，所以计算时要求整体刚度矩阵不能奇异，对于一些接触高度非线性问题，有时无法保证收敛。