芯耀
与非AI
133
引言
前文我们主要以中心差分法为例讨论了显式算法的原理、稳定性和精度。本文我们主要讨论结构动力学中时域直接积分法的隐式格式，主要介绍Newmark隐式算法，此法在结构动力学当中应用非常广泛，当前的主流有限元软件如：ansys、abaqus都有应用。
算法原理
中心差分法的核心思想是通过t时刻和t-∆t时刻的位移、速度和加速度递推t+∆t时刻的位移，是两步显式的；而Newmark算法的核心思想是通过t+∆t时刻的加速度递推t+∆t时刻的位移和速度，是单步隐式的，递推方程为：
速度v(t+∆t)=u'(t+∆t)= u'(t)+[(1-δ)*u''(t)+ δ*u''(t+∆t)]∆t
位移u(t+∆t)= u(t)+u'(t)∆t +[(1/2-α)* u''(t)+ α* u''(t+∆t)]∆t2
方程两边同时存在t+∆t时刻的位移、速度和加速度，可以看出是隐式单步递推格式。
在t~ t+∆t时间内的加速度是t时刻和t+∆t时刻加速度的加权和。当δ=1/2，α=1/6时，t~ t+∆t时间内的加速度为平均加速度：
u''(t+τ)=[u''(t)+ u''(t+∆t)]/2 (0<τ<∆t)
即为常值加速度法，此法是无条件稳定的，但是精度较差。为了提高精度，Newmark法采取牺牲稳定性提高精度的做法。
为了将Newmark算法显式化，首先联立速度和位移的递推方程，用t+∆t时刻的位移表示t+∆t时刻的加速度和速度，然后再代入t+∆t时刻的动力学方程：
Mu''(t+∆t)+Cu'(t+∆t)+Ku(t+∆t)=F(t+∆t)
得到位移的单步递推方程：
[1/(α∆t2)*M+δ/(α∆t)*C+K)*u(t+∆t)=F(t+∆t)+
[(1/(α∆t2)*u(t)+ (1/(α∆t)*u'(t)+(1/(2α)-1)*u''(t)]*M+
[δ/(α∆t)*u(t)+( δ/α-1)*u'(t)+ (δ/(2α)-1)*∆t2*u''(t)]*C
Newmark算法的详细过程如下：
(1) 计算初始条件以及参数：
a) 形成M、C、k矩阵；
b)根据初始条件u(t=0)和u'(t=0)，由t=0时刻的动力学方程求出u''(t=0)；
u''(0)=M-1[F(0)-C u'(0)-Ku(0)]
c) 选择合适的积分时间步长∆t和参数δ(δ≥0.5)、α(α≥0.25*(0.5+δ)2)，并计算积分常数；
c0=1/(α∆t2) c1=δ/(α∆t) c2=1/(α∆t) c3=1/(2α)-1
c4=δ/α-1 c5= ∆t/2*(δ/α-2) c6=∆t(1-δ) c7=δ∆t
d) 三角分解等效刚度矩阵Keqv
Keqv= K+c0*M+c1*C=LDLT；
(2) 循环每一时刻计算(t=0, ∆t, 2∆t,…)：
a) 计算t+∆t时刻的等效载荷Feqv(t+∆t)：
Feqv(t+∆t)=F(t+∆t)+
M*[c0u(t)+c2u'(t)+c3u''(t)] +
C*[c1u(t)+c4u'(t)+c5u''(t)]
b) 计算t+∆t时刻的位移u(t+∆t)；
Keqv *u(t+∆t)=LDLT*u(t+∆t)= Feqv(t+∆t)
c)计算t+∆t时刻的加速度u''(t+∆t)和速度u'(t+∆t)；
u''(t+∆t)=c0[u(t+∆t) -u(t)]- c2u'(t) – c3u''(t)
u'(t+∆t)= u'(t)+ c6u''(t)+ c7u''(t+∆t)
关于Newmark法的几点说明：
(1) 隐式算法：由于刚度矩阵K一般不是对角阵，线性分析时，由于时间步长∆t可以只有一个，因此只需要对等效刚度矩阵Keqv进行一次三角分解即可，避免求逆；但是非线性分析时间步长∆t可以是变化的，则需要对Keqv进行多次三角分解；
(2) 无条件稳定：当积分参数满足δ≥0.5，α≥0.25*(0.5+δ)2时，是无条件稳定的，即∆t不会影响稳定性，只需要考虑∆t对数值精度的影响，一般取∆t=(0.1~0.05)Tmin。实际问题中采用的Tmin比系统的最小周期大很多倍，因此相比中心差分法，在保证精度的同时，Newmark法是以对等效刚度矩阵K求逆降低分析效率为代价，换取比条件稳定的中心差分法更大的时间步长∆t，同时较大的∆t可以过滤高频响应，因此，Newmark算法更适合于时间较长的动力学分析。
稳定性
对于稳定性，我们参考前文显式分析的讨论方式，核心在于推导出特征方程。
Newmark算法的递推公式：
[1/(α∆t2)*M+δ/(α∆t)*C+K)*u(t+∆t)=F(t+∆t)+
[(1/(α∆t2)*u(t)+ (1/(α∆t)*u'(t)+(1/(2α)-1)*u''(t)]*M+
[δ/(α∆t)*u(t)+( δ/α-1)*u'(t)+ (δ/(2α)-1)*∆t2*u''(t)]*C
我们令F(t+∆t)=0,C=0，第r阶模态坐标qr，固有频率为ωr，则递推公式可改写为：
(1+α∆t2ωr2)*qr(t+∆t)=qr(t)+∆t*qr'(t)+(1/2-α)∆t2*qr''(t)
以上公式为单步法，缺qr(t-∆t)构成两步法，为此需要引入t时刻的位移速度递推公式：
qr'(t)= qr'(t-∆t)+[(1-δ)*qr''(t-∆t)+ δ*qr''(t)]∆t
qr(t)= qr(t-∆t)+qr'(t-∆t)∆t+[(1/2-α)*qr''(t-∆t)+α*qr''(t)]∆t2
结合模态坐标下的无阻尼自由振动方程：
qr''(t)+ ωr2qr(t)=0
换掉递推方程等号右边的qr'(t)和qr''(t)得到两步递推公式，同时引入p=(∆t*ωr)2，得到特征方程：
λ2(1+αp)+λ[-2+(1/2-2α+δ)p]+[1+(1/2+α-δ)p]=0
两个特征根：
λ1,2=(2-g)/2 ± sqrt[(2-g)2-4(1+h)]/2
g=(1/2+δ)p/(1+αp) h=(1/2-δ)p/(1+αp)
由于实际系统多数为欠阻尼的震动系统，具有震荡特性，因此要求λ为复数，要求(2-g)2-4(1+h)<0，即：
p[4α-(1/2+δ)2]>-4
因为要求无条件稳定，因此我们可以要求:
α≥(1/2+δ)2/4
除此之外，稳定条件下，响应的幅值不应该无限放大，因此要求|λ|≤1，即-1≤h≤0，即：
δ≥1/2
因此Newmark算法的无条件稳定条件为：
δ≥1/2 α≥(1/2+δ)2/4
如果不满足无条件稳定的条件，我们还想得到稳定的解答，就需要:
∆t<∆tcr=Tr/π/sqrt[(1/2+δ)2-4α]
只要 (2-g)2-4(1+h)<0，λ就为复数，且只要δ=1/2，就始终满足|λ|=1，说明幅值保持不变，符合无阻尼自由振动的响应情况。但如果一旦δ>1/2，就会有 |λ|<1，说明幅值会不断衰减，但是物理系统为无阻尼系统，幅值不应该衰减，这种因为参数选取给系统带来额外阻尼效应的现象称为“数值阻尼”。
最后讨论“数值阻尼”对结构高频和低频部分响应的影响，在δ≥0.5时，∆t/T和幅值比|λ|之间的函数关系并画出函数图像：
|λ|=sqrt(1+h)=sqrt[1+(0.5-δ)*4*π2*(∆t/T)2/(1+α*4*π2*(∆t/T)2)]
从图上可以得到几点结论：
a) δ=0.5时，只要α≥0.25，幅值比|λ|=1，说明无数值阻尼；
b) δ>0.5时，就存在“数值阻尼”，α的取值对低频部分基本无影响(|λ|接近1)；但对于高频部分(λ|<1)，α越小，“数值阻尼”效应越明显，幅值衰减越快；
c) 对于一般的非高速冲击的问题，因为载荷持续时间长，且载荷的频率成分主要集中在低频段，从模态叠加的角度上看，高频段对结果贡献较小，选用Newmark算法并不会对结果精度有较大的影响；
d) 对于高速冲击类的问题，一般都伴随有强烈非线性，冲击持续时间较短，载荷的频率成分较宽，高频段的结果不能忽略，Newmark算法会衰减高频响应，造成精度下降，因此选用中心差分法较为合适。
算例
我们仍然用显式算法当中的算例，即一个三自由度MCK系统在阶梯载荷作用下的响应，运动微分方程为：
Mu''(t)+Ku''(t)=F(t) u(0)=0 u'(0)=0
M=[1,0,0;0,3,0;0,0,1] C=0
K=[2,-1,0;-1,4,-2;0,-2,2] F(t)=[0,0,6]
分别列时间出积分步长∆t=Tmin/20、∆t=Tmin/10、∆t=5Tmin情况下的精确解、显式解和隐式解并进行精度对比。其中Tmin=2*π/ωmax
a) 当∆t=Tmin/20和∆t=Tmin/10时，此时(∆t<∆tcr)，隐式算法和显式算法的精度相差不大，且都是开始时间段精度较高，∆t越小，保持精度较高的时间段越长；
b) ∆t=5Tmin 时，此时(∆t>∆tcr)显式算法是不稳定的，但隐式算法依然稳定，但精度较差，因此也可以说明：隐式算法的无条件稳定性，∆t只决定计算精度。
最后
本文对动力学分析中以Newmark法为代表的隐式积分法进行了简单讨论，重点在于算法原理、分析过程、稳定性等方面，并和精确解以及中心差分法进行了精度对比分析，旨在使读者在解决实际工程问题的时候，选择合适的分析方法，限于作者水平有限，观点仅供参考。
