结构动力学中的常用数值方法

第五章 结构动力学中的常用数值方法5．1．结构动力响应的数值算法....0()(0)(0)M x c x kx F t x a x v ⎧++=⎪⎪=⎨⎪=⎪⎩当c 为比例阻尼、线性问题→模态叠加最常用。

但当C 无法解耦，有非线性存在，有冲击作用（激起高阶模态，此时模态叠加法中的高阶模态不可以忽略）。

此时就要借助数值积分方法，在结构动力学问题中，有一类方法称为直接积分方法最为常用。

所识直接是为模态叠加法相对照来说，模态叠加法在求解之前，需要对原方程进行解耦处理，而本节的方法不用作解耦的处理，直接求解。

（由以力学，工程中的力学问题为主要研究对象的学者发展出来的）中心差分法的解题步骤1. 初始值计算（1） 形成刚度矩阵K ，质量矩阵M 和阻尼矩阵C 。

（2） 定初始值0x ，.0x ，..0x 。

（3） 选择时间步长t ∆，使它满足cr t t ∆<∆，并计算 021()a t =∆，112a t=∆，202a a =（4） 计算...0011122t x x x x a a -∆=-+（5） 形成等效质量阵01M a M a C -=+ （6） 对M -阵进行三角分解T M LDL -= 2．对每一时间步长（1） 计算时刻t 的等效载荷21()()t t t tt Q Q K a M x a Ma C x --∆=---- （2） 求解t t +∆时刻的位移 ()Tt t t L D L x Q -+∆=（3） 如需要计算时刻t 的速度和加速度值，则.1()t t t t t x a x x +∆-∆=-..0(2)t t t t t t x a x x x +∆-∆=-+若系统的质量矩阵和阻尼矩阵为对角阵时，则计算可进一步简化。

纽马克法的解题步骤1.初始值计算（1）形成系统刚度矩阵K ，质量矩阵M 和阻尼矩阵C（2）定初始值0x ，.0x ，..0x 。

（3）选择时间步长t ∆，参数γ、σ。

并计算积分常数0.5γ≥，20.25(0.5)δγ≥+ 021a tδ=∆，1a tγδ=∆，21a tδ=∆3112a δ=-，41a γδ=-，5(2)2t a γδ∆=- 6(1)a t γ=∆-，7a t γ=∆（4）形成等效刚度矩阵K -01K K a M a C -=++（5）K -矩阵进行三角分解 T K L D L -= 2. 对第一时间步长（1）计算t t +∆时刻的等效载荷......623145()()t t t t t t t t t Q Q M a x a x a x C a x a x a x -+∆=++++++ （2）求解t t +∆时刻的位移 ()Tt t t t L D L x Q -+∆+∆=（3）计算t t +∆时刻的加速度和速度.....023()t t t t t t t x a x x a x a x +∆+∆=---......67t t t t t t x x a x a x +∆+∆=++威尔逊-θ法的解题步骤1. 初始值计算（1）形成系统刚度矩阵K ，质量矩阵M 和阻尼矩阵C（2）定初始值0x ，.0x ，..0x 。

（3）选择时间步长t ∆，并计算积分常数1.4θ= 026()a t θ=∆，13a tθ=∆，212a a =32ta θ∆=，04a a θ=，25a a θ-=631a θ=-，72t a ∆=，286t a ∆=（4）形成等效刚度K -01K K a M a C -=++（5）将等效刚度K -进行三角分解 T K L D L -= 2.对每一个时间步长（1）计算t t +∆时刻的等效载荷...02()(2)t t t t t t t t t R Q Q Q M a x a x x θθ-+∆+∆=+-+++...13(2)t t t C a x x a x +++ （2）求解t t +∆时刻的位移()Tt t t t L D L x R θθ-+∆+∆=（3）计算在t t +∆时刻的加速度、速度和位移.....456()t t t t t t t x a x x a x a x θ+∆+∆=-++5．2 结构动力响应数值算法性能分析对公式（5.1）描述的线性系统结构动力学问题，已经有证明对整个多自由度的积分，等价于将模态分解后对单自由度的积分的结果进行模态叠加，因此可以通过对单自由度问题的分析，来说明算法的特性，其中阻尼均假设为比例阻尼，这样，模态分解后的单自由度结构动力学方程为：)(22t f x x x=++ωξω （5-29） 以下算法的性能分析，均将算法用于这个方程。

5．2．1 算法用于结构动力学方程的有限差分表示将数值计算方法应用于（5-29）, 即分别在相邻的不同时刻应用算法可得如下一般形式k k k L Ay y +=+1 （5-30）A 为放大矩阵或称逼近算子，k L 为载荷逼近算子。

{}{}Tm k k k k Tm k k k k x x x y x x x y -+++--== ,,,,,1111例如将Newmak 方法应用于方程（5-29）有：d t A A A 1-= （5-33）⎢⎢⎣⎡+=2221γωδωh h A t⎥⎥⎦⎤+γξωδξωh h 2122，⎢⎢⎢⎣⎡----=222)1()21(21ωγωδh hA d⎥⎦⎤----ξωγξωδ)1(21])21(1[h h h矩阵A 的特征多项式为02)det(212=+-=--A A I A λλλ （5-34）其中A 1，A 为该矩阵的两个特征向量，分别为矩阵的迹的一半和矩阵的行列式)(212122111A A traceA A +==（5-35） 211222112det A A A A A A -==（5-36）22112,1A A A -±=λ （5-37）对Newmak 方法有：Dv A ])412()12(1[21Ω--+Ω-+=δξγ,DA ])21()22(1[22Ω+-+Ω-+=γδξγ （5-38）其中 h 为时间步长，221,Ω+Ω+==ΩδγξωD h 。

Newmak 方法放大矩阵的规模是二维的，因此特征值也只有两个，可以根据它们进行分析。

有的算法放大矩阵是三维的，例如Wilson-θ方法,在无阻尼情况下放大矩阵为：⎢⎢⎢⎣⎡-+Ω-=Ω-22366)1(13ωθθωD A ωθθθθθΩ+-ΩΩ-Ω+h h 6)6)3(()6(22232⎥⎥⎥⎥⎦⎤-Ω-Ω+Ω--Ω+Ω--Ω+)636)2/336()26(2232222322232θθθθθθθθθh h （5-39）θθ)6(22+Ω=D放大矩阵A 的特征多项式为：2)det(32213=-+-=--A A A I A λλλλ （5－40）其中A 1，A 2，A 3为该矩阵的三个特征向量，分别为矩阵的迹的一半、各阶主子式的和以及矩阵的行列式，对Wilson-θ方法有⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+ΩΩ+Ω+Ω-Ω-+-=+Ω-+Ω-Ω+Ω=+ΩΩ-Ω-Ω-Ω+-=)6(3366)6(1218634)6(233361822232222322223222222222321θθθθθθθθθθθθθθc A A c c A （5-41）此外，在几个不同时刻应用数值算法，然后将方程中的速度和加速度项消去，可得数值算法关于位移的差分方程，例如Newmak 方法，有nn x x ])412()12(1[2)21(212Ω--+Ω-+-Ω+Ω++γδξγδγξ])21()22(1[12=Ω+-+Ω-++-n x γδξγ （5-42）很显然，其特征方程与其放大矩阵A 的特征方程是相同的，使用关于位移的线性多步方式和放大矩阵来说明算法性能是一样的，只不过各有方便之处。

5．2．2 算法的稳定性分析设m i i 2,1,=λ为放大矩阵A 的特征值，则i λρmax =定义为A 的谱半径，若特征值互异，则1≤ρ的算法是稳定的，但若有重特征根，则要求1<ρ。

如果算法的稳定性要求对步长的选取有限制，称算法是有条件稳定的，反之为无条件稳定的。

放大矩阵的谱半径小于等于1成立的充分条件是⎪⎩⎪⎨⎧≥-≥++≥+-0102102122121A A A A A （5-43）对33⨯的放大矩阵⎪⎪⎪⎩⎪⎪⎪⎨⎧≥-+-≥+++≥--+≥+--≥-+-0)2(10210323032302213132321321321321A A A A A A A A A A A A A A A A （5-44）上两式是关于算法自由参数Ω,ξ的不等式，由它可以判断算法是否无条件稳定，若不是，将给出稳定条件。

例5－1分析Newmak 方法、Wilson-θ方法的稳定性 解： 将（5-38）代入（5-43）有)21(22≥-Ω+Ωγξ1)21()2(2≤-Ω-+-Ωξγδγ显然，当2,21γδγ≥≥（5-45）算法无条件稳定。

当2,21γδγ<≥且δγγξδγγξ--+-+-=Ω≤Ω2])21(2[)21(2122c （5-46）算法稳定，但为条件稳定，其中c Ω为临界采样频率。

由于（5-43）式仅仅是充分条件，所以可进一步按照稳定性的定义得到5.1.2节叙述的无条件稳定条件。

对Wilson-θ方法，将（5-41）代入（5-44）得⎪⎪⎪⎩⎪⎪⎪⎨⎧≥-+Ω≥-+Ω-+≥-+-Ω+≥-Ω≥Ω0)312(01224)614(0)661(120)12(606242232222c c c θθθθθ （5-47）容易看出，其中第一，二，五不等式恒成立，对第三，四不等式若希望对任意的Ω均成立，则有：⎪⎩⎪⎨⎧≥+-≥-+-01640661232θθθθ求解上述不等式得37.1231≈+≥θ （5-48）实际使用中通常选取θ=1.4 5.2.3 算法的相容性和收敛性直接积分算法的相容性、收敛性分析同样要使用其位移型的差分方程，或对应的单步多值形式。

在算法（5-30）式中，用精确解代替近似解，即可得到局部截断误差表达式，用符号)(k t e 表示)()()(1k k k k t he L t Ay t y ++=+（5-49）局部截断误差表达式用放大矩阵的特征量以最常用的线性三步法为例可表示为2321/)]2()()(2)([)(hh t x A h t x A t x A h t x t e k k k k k ---+-+= （5-50）其中321,,A A A 分别为对应的33⨯的放大矩阵的三个特征向量，然后将),(h t x k - )2(h t x k -在k t 点进行泰勒展开，然后利用运动平衡方程化简即可。

若局部截断误差表达式为步长的s （s>0）阶小量，则称算法是s 阶相容的。