1.1. 弹性动力学有限元基本解法
结构系统的通用运动学方程为：
+C U +KU =R （1） M U t
求解该动力学振动响应主要有三类方法：（1）时域法（2）频域法（3）响应谱法
时域法又可分为：（1）直接积分法，（2）模态叠加法。直接积分法又可分为中心差分法（显式），Wilson θ（隐式）法以及Newmark （隐式）法等。
本文介绍中心差分法（显式）与Newmark （隐式）法。
1 中心差分法（显式）
假定0，t 1，t 2，…，t n 时刻的节点位移，速度与加速度均为已知，现求解t n (t +∆t ) 时刻的 结构响应。中心差分法对加速度，速度的导数采用中心差分代替，即为：
=1(U -2U +U ) U t t -∆t t t +∆t ∆t 2
1 U t =(U t +∆t -U t -∆t ) (2) 2∆t
将（2）式代入（1）式后整理得到
ˆU t +∆t =R ˆt （3） M
式（3）中
ˆ=1M +1C M ∆t 22∆t
ˆt =R t -(K -2M ) U t -(1M -1C ) U t -∆t R ∆t 2∆t 22∆t
分别称为有效质量矩阵，有效载荷矢量。R ，M ，C ，K 为结构载荷，质量，阻尼，刚度矩阵。
求解线性方程组（3），即可获得t +∆t 时刻的节点位移向量U t +∆t ，将U t +∆t 代回几何方程与物理方程，可得t +∆t 时刻的单元应力和应变。
ˆ，有效载荷矢量R ˆ，即可中心差分法在求解t +∆t 瞬时的位移U t +∆t 时，只需t +∆t 时刻以前的状态变量U t 和U t -∆t ，然后计算出有效质量矩阵M t
求出U t +∆t ，故称此解法为显式算法。
中心差分法，在开始计算时，需要仔细处理。t =0时，要计算U ∆t ，需要知道U -∆t 的值。因此应该有一个起始技术，因而该算法不是自动起步
，U 是已知的，由t =0时的（2）式可知： 的。由于U ，U 000
U -∆t 2∆t + =U 0-∆t U U 002
中心差分法中时间步长∆t 的选择涉及两个方面的约束：数值算法的稳定性和计算时间。中心差分法的实质是用差分代替微分，并且对位移和加速度的导数采用线性外插，这限制了∆t 的取值不可过大，否则结果可能失真过大。
可以证明：中心差分法是条件稳定的。即当时间步长∆t 必须小于由该问题求解方程性质所决定的一个时间步长的临界值。LS-DYNA 中，采用“变时间步长法”，即每一时刻的步长∆t 由当前结构的稳定性条件来控制。具体算法为：计算每一个单元的极限时间步长∆t ei ，i =1,2…，取∆t =min(∆t ei ) 为下一个时刻的时间步长。各种单元的∆t ei 计算方法如下。
（1）1D 杆，梁单元
L ∆t e =α C
其中α为时间步长因子，系统默认为0.9。L 为杆，梁单元的长度。C =
（2）2D 板，壳单元
L ∆t e =αmin C
其中α为时间步长因子，系统默认为0.9。L min 为 壳单元的最小单元边长度。C =
（3）3D 单元
L e ∆t e =α Q +(Q 2+C 2) 1/2
其中
kk (ε kk
C 0和C 1为无量纲常数，默认C 0=1.5，C 1=0.06. E ρ材料声速。 E (1-ν) ρ2为材料的声速。
⎧V e (对8节点体单元）⎪ L e 为单元等效长度， V e 为单元体积，A e max 为单元最大侧面积。 L e =⎨A e max ⎪L （对4节点体单元）⎩min
E （1-）为材料声速。 (1+ν（) 1-2ν）ρ
时间步因子α可由用户设置，减小α相当于减少时间步长。设置时间步长因子α的关键字为*CONTROL_TIMESTEP，控制参数为TSSFAC 。 另外，质量缩放可以人为控制时间步长。即调整单元密度ρ来改变时间步长。以壳单元为例说明质量缩放改变时间步长。 C =
L min （假如α=1） C
E C =2(1-ν) ρ∆t e =α
(∆t specified ) 2E (1-ν2) ρi 由(得到ρi =2 ) =2L min E L min (1-ν)
LS-DYNA 中有2种质量缩放方案，修改*CONTROL_TIMESTEP中的参数。
方案1：DT2MS 为正的时间步
通过调整单元密度，使所有单元时间步相同，只用于惯性效应不重要的情况。
方案2：DT2MS 为负的时间步 质量缩放只用于小于指定时间步长DT 的单元。惯性效应应该通过变形体的动能与内能的比例进行衡量（一般应该小于10%等）。 ∆t specified 2
2 Newmark法（隐式）
Newmark 假定在时间间隔[t , t +∆t ]内，加速度线性变化，即采用如下的加速度，速度公式：
+[(1-δ) U +δU ]∆t U =U t +∆t t t t +∆t
∆t +[(1-α) U +αU ]∆t 2 (4) U t +∆t =U t +U t t t +∆t 2
式中α，δ为按积分的精度和稳定性要求可以调整的参数。
和U 用U t +∆t ，U ，U 表示的表达式，代入（1）式中整理得到 根据（4）式可给出U t +∆t t +∆t t t
ˆU =R ˆ （5） K t +∆t t +∆t
其中
1δM +C +K α∆t 2α∆t
+(1-1) U ]+C [δU +(δ-1) U +(δ-1) ∆t U ] 称之为有效刚度矩阵和有效载荷矢量。由上式可以看出求解ˆt +∆t =R t +∆t +M [1U t +1U R t t t t t 2α∆t α∆t 2αα∆t α2α
当前U t +∆t ，需要用到当前时刻的R t +∆t ，因此该算法为隐式算法。 ˆ=K
当载荷历史全部已知时，F t +∆t 为已知量，求解需要迭代实现。
可以证明，当参数δ≥0. 5，α≥0. 25(0. 5+δ) 2时，Newmark 法是无条件稳定的，即∆t 的大小不影响数值稳定性。此时时间步长∆t 的选择主要根据解得精度确定。一般，Newmark 法可以比中心差分法的时间步长大得多。
3．结论
比较两种算法，显式中心差分法非常适合研究波的传播问题，如碰撞、高速冲击、爆炸等。分析式（3）发现，显式中心差分法的M 与C 矩阵是对角阵，如给定某些有限元节点以初始扰动，在经过一个时间步长后，和它相关的节点进入运动，即U 中这些节点对应的分量成为非零量，此特点正好和波的传播特点相一致。另一方面，研究波传播的过程需要微小的时间步长，这也正是中心差分法的特点。
而Newmark 法更加适合于计算低频占主导的动力问题，从计算精度考虑，允许采用较大的时间步长以节省计算时间，同时较大的时间步长还可以过滤掉高阶不精确特征值对系统响应的影响。隐式方法要转置刚度矩阵，增量迭代，通过一系列线性逼近（Newton-Raphson ）来求解。正因为隐式算法要对刚度矩阵求逆，所以计算时要求整体刚度矩阵不能奇异，对于一些接触高度非线性问题，有时无法保证收敛。
1.1. 弹性动力学有限元基本解法
结构系统的通用运动学方程为：
+C U +KU =R （1） M U t
求解该动力学振动响应主要有三类方法：（1）时域法（2）频域法（3）响应谱法
时域法又可分为：（1）直接积分法，（2）模态叠加法。直接积分法又可分为中心差分法（显式），Wilson θ（隐式）法以及Newmark （隐式）法等。
本文介绍中心差分法（显式）与Newmark （隐式）法。
1 中心差分法（显式）
假定0，t 1，t 2，…，t n 时刻的节点位移，速度与加速度均为已知，现求解t n (t +∆t ) 时刻的 结构响应。中心差分法对加速度，速度的导数采用中心差分代替，即为：
=1(U -2U +U ) U t t -∆t t t +∆t ∆t 2
1 U t =(U t +∆t -U t -∆t ) (2) 2∆t
将（2）式代入（1）式后整理得到
ˆU t +∆t =R ˆt （3） M
式（3）中
ˆ=1M +1C M ∆t 22∆t
ˆt =R t -(K -2M ) U t -(1M -1C ) U t -∆t R ∆t 2∆t 22∆t
分别称为有效质量矩阵，有效载荷矢量。R ，M ，C ，K 为结构载荷，质量，阻尼，刚度矩阵。
求解线性方程组（3），即可获得t +∆t 时刻的节点位移向量U t +∆t ，将U t +∆t 代回几何方程与物理方程，可得t +∆t 时刻的单元应力和应变。
ˆ，有效载荷矢量R ˆ，即可中心差分法在求解t +∆t 瞬时的位移U t +∆t 时，只需t +∆t 时刻以前的状态变量U t 和U t -∆t ，然后计算出有效质量矩阵M t
求出U t +∆t ，故称此解法为显式算法。
中心差分法，在开始计算时，需要仔细处理。t =0时，要计算U ∆t ，需要知道U -∆t 的值。因此应该有一个起始技术，因而该算法不是自动起步
，U 是已知的，由t =0时的（2）式可知： 的。由于U ，U 000
U -∆t 2∆t + =U 0-∆t U U 002
中心差分法中时间步长∆t 的选择涉及两个方面的约束：数值算法的稳定性和计算时间。中心差分法的实质是用差分代替微分，并且对位移和加速度的导数采用线性外插，这限制了∆t 的取值不可过大，否则结果可能失真过大。
可以证明：中心差分法是条件稳定的。即当时间步长∆t 必须小于由该问题求解方程性质所决定的一个时间步长的临界值。LS-DYNA 中，采用“变时间步长法”，即每一时刻的步长∆t 由当前结构的稳定性条件来控制。具体算法为：计算每一个单元的极限时间步长∆t ei ，i =1,2…，取∆t =min(∆t ei ) 为下一个时刻的时间步长。各种单元的∆t ei 计算方法如下。
（1）1D 杆，梁单元
L ∆t e =α C
其中α为时间步长因子，系统默认为0.9。L 为杆，梁单元的长度。C =
（2）2D 板，壳单元
L ∆t e =αmin C
其中α为时间步长因子，系统默认为0.9。L min 为 壳单元的最小单元边长度。C =
（3）3D 单元
L e ∆t e =α Q +(Q 2+C 2) 1/2
其中
kk (ε kk
C 0和C 1为无量纲常数，默认C 0=1.5，C 1=0.06. E ρ材料声速。 E (1-ν) ρ2为材料的声速。
⎧V e (对8节点体单元）⎪ L e 为单元等效长度， V e 为单元体积，A e max 为单元最大侧面积。 L e =⎨A e max ⎪L （对4节点体单元）⎩min
E （1-）为材料声速。 (1+ν（) 1-2ν）ρ
时间步因子α可由用户设置，减小α相当于减少时间步长。设置时间步长因子α的关键字为*CONTROL_TIMESTEP，控制参数为TSSFAC 。 另外，质量缩放可以人为控制时间步长。即调整单元密度ρ来改变时间步长。以壳单元为例说明质量缩放改变时间步长。 C =
L min （假如α=1） C
E C =2(1-ν) ρ∆t e =α
(∆t specified ) 2E (1-ν2) ρi 由(得到ρi =2 ) =2L min E L min (1-ν)
LS-DYNA 中有2种质量缩放方案，修改*CONTROL_TIMESTEP中的参数。
方案1：DT2MS 为正的时间步
通过调整单元密度，使所有单元时间步相同，只用于惯性效应不重要的情况。
方案2：DT2MS 为负的时间步 质量缩放只用于小于指定时间步长DT 的单元。惯性效应应该通过变形体的动能与内能的比例进行衡量（一般应该小于10%等）。 ∆t specified 2
2 Newmark法（隐式）
Newmark 假定在时间间隔[t , t +∆t ]内，加速度线性变化，即采用如下的加速度，速度公式：
+[(1-δ) U +δU ]∆t U =U t +∆t t t t +∆t
∆t +[(1-α) U +αU ]∆t 2 (4) U t +∆t =U t +U t t t +∆t 2
式中α，δ为按积分的精度和稳定性要求可以调整的参数。
和U 用U t +∆t ，U ，U 表示的表达式，代入（1）式中整理得到 根据（4）式可给出U t +∆t t +∆t t t
ˆU =R ˆ （5） K t +∆t t +∆t
其中
1δM +C +K α∆t 2α∆t
+(1-1) U ]+C [δU +(δ-1) U +(δ-1) ∆t U ] 称之为有效刚度矩阵和有效载荷矢量。由上式可以看出求解ˆt +∆t =R t +∆t +M [1U t +1U R t t t t t 2α∆t α∆t 2αα∆t α2α
当前U t +∆t ，需要用到当前时刻的R t +∆t ，因此该算法为隐式算法。 ˆ=K
当载荷历史全部已知时，F t +∆t 为已知量，求解需要迭代实现。
可以证明，当参数δ≥0. 5，α≥0. 25(0. 5+δ) 2时，Newmark 法是无条件稳定的，即∆t 的大小不影响数值稳定性。此时时间步长∆t 的选择主要根据解得精度确定。一般，Newmark 法可以比中心差分法的时间步长大得多。
3．结论
比较两种算法，显式中心差分法非常适合研究波的传播问题，如碰撞、高速冲击、爆炸等。分析式（3）发现，显式中心差分法的M 与C 矩阵是对角阵，如给定某些有限元节点以初始扰动，在经过一个时间步长后，和它相关的节点进入运动，即U 中这些节点对应的分量成为非零量，此特点正好和波的传播特点相一致。另一方面，研究波传播的过程需要微小的时间步长，这也正是中心差分法的特点。
而Newmark 法更加适合于计算低频占主导的动力问题，从计算精度考虑，允许采用较大的时间步长以节省计算时间，同时较大的时间步长还可以过滤掉高阶不精确特征值对系统响应的影响。隐式方法要转置刚度矩阵，增量迭代，通过一系列线性逼近（Newton-Raphson ）来求解。正因为隐式算法要对刚度矩阵求逆，所以计算时要求整体刚度矩阵不能奇异，对于一些接触高度非线性问题，有时无法保证收敛。