结构动力学方程常用数值解法
对于一个实际结构，由有限元法离散化处理后，动力学方程可写为：
Mx?Cx?Kx?F(t)
从数学角度看，这是一个常系数的二阶线性常微分方程组，计算数学领域，常微分数值算法常用的有两大类：－、针对一阶微分方程数值积分法发展的欧拉法，中点法，Rugge-kutta（龙格—库塔）方法。二、直接基于二阶动力学方程发展的方法。
对结构动力学问题的数值求解，常用的有两大类：一是坐标变换法，它是对结构动力方程式，在求解之前，进行模态坐标变换，实际上就是一种Rize变换，即把原物理空间的动力方程变换到模态空间中去求解。现在，普遍使用的方法是模态（振型）迭加法。二是直接积分法，它是对结构动力方程式在求解之前不进行坐标变换，直接进行数值积分计算。这种方法的特点是对时域进行离散，然后将该时刻的加速度和速度用相邻时刻的各位移线性组合而成。通常又称为逐步积分法。
模态迭加方法，比较常用，但如下情况通常使用直接积分方法（即求解之前不进行模态分析）一、非比例阻尼，非线性情况。二、有冲击作用，激起高频模态，力作用持续时间较短，模态迭加计算量太大。
...一 振型迭加法与Duhamel积分数值解
按照有限单元法的一般规则, 经过边界条件的约束处理, 结构在强迫振动时多自由度体系的运动平衡方程可以表示为:
???CU??KU?RMU (1)
其中, M是体系的质量矩阵, C 是体系的阻尼矩阵, 而K 则是刚度矩阵. R 为外荷载向
???量. U、U和U则分别是体系单元节点的位移、速度和加速度向量. 上述动力平衡方程???实质上是与加速度有关的惯性力MU和与速度有关的阻尼力CU及与位移有关的弹性力
KU在时刻t与荷载的静力平衡。
振型叠加法是把多自由度体系的结构的整体振动分解为与振型次数相对应的单自由度体系, 求得各个单自由度体系的动力响应后, 再进行叠加得出结构整体响应. 振型叠加法原理是利用结构无阻尼自由振动的振型矩阵作为变换矩阵, 将结构动力方程式(1)式变换成一组非耦合的微分方程. 逐个地求解这些方程后, 将解叠加即可得到动力方程的解。
将体系单元节点的位移向量表示为如下的变换形式：
U(t)??X(t)
（2）
式中的变换矩阵?是由动力方程对应的无阻尼自由振动方程解出的前m阶振型矩阵.即
??[???m];X(t)是与时间有关的m阶向量, X的各分量称为广义位移。 12...将式（2）代入动力方程（1）并左乘以?，则可得广义位移为未知数的方程：
T??(t)?CX?(t)?KX(t)?R(t)MX
（3） 式中
TM??TM?，C??C?，K??TK?，R??TR
（4）
现在进一步考察式（4） . 考虑到特征向量的正交性, 可得
?TM??I，?TK???
（5）
于是对应于振型的广义位移的平衡方程( 3) 可改写为
??(t)??TC?X?(t)??X(t)??TR(t)X
（6）
其中，?为特征值
??12???2?2????...????2???i??...??2?m????
（7）
将式( 2) 稍加运算可得广义位移用有限元位移表示的形式
X??TMU
（8）
在( 6) 式中, 当忽略了阻尼的影响, 平衡方程为互不耦合的, 可以对每个方程逐个地进行时间积分. 出于相同的考虑, 在对有阻尼的体系进行分析时仍然希望采用相同的计算过程去求解互不耦合的平衡方程式. 问题是式( 6) 中的阻尼阵C 通常不能象体系的质量阵和刚度阵那样由单元的刚度阵和质量阵装配而成. 但当假定阻尼与固有频率成比例, 即假定
?iTC?j?2??ii?ij（9）
式中，i是振型阻尼参数;
??ij??是Kronecker符号( 当i?j时，ij= 1.当i?j时, ij= 0) 。
这时式( 6) 可简化为如下形式的若干个方程式
2??i(t)?2???xiix(t)??ixi(t)?ri(t)
（10）
其中xi(t)的初始条件为下式
xit?0??iTMU0，
?ixt?0??iTMU0
（11）
2?i式（10）表示了一个具有单位质量，刚度为的自由度体系当阻尼比为?i时的运动平衡
控制方程。这个平衡方程的求解可通过计算Duhamel积分求得。
xi(t)?（12）
式中
1?i?t0ri(t)e??i?i?t???sin?i?t???d??e??i?i(?isin?it??icos?it)t
?i??i1??i2 （13）
当利用式( 9) 来考虑阻尼的影响时意味着假设结构的总阻尼是每个振型的阻尼之和， 而每个振型上的阻尼是能够量测的，况且在大多数情况下结构的阻尼比更易于量测。因而便于用来近似地反映结构体系的阻尼特性。 同时在计算上也避免计算阻尼阵而只需计算刚度阵和质量阵。 积分递推公式 对以上方程式( 10) ，考虑某一模态的振动, 并略去下标i可写为
??(t)?2??x?(t)??2x(t)?r(t)x （14）
在初始条件
xt?t0??x0x，
t?t0?0?x （15）
下的定解为
x(t)?e????t?t0???[sin??t?t0??cos??t?t0?]?x0?
1????t?t0?1?esin??t?t0??x0??r(?)e????t???sin??t???d???t0t
（16） 式中，
???1??2，将上式对时间求导，得
?(t)?ex????t?t0????2?2?[????sin??t?t0?]?x0???
[cos??t?t0???e????t?t0????0sin??t?t0?]?x?
?1?t0????t???r(?)e[???sin??t?????cos??t???]d??t （17）
将上式中的t0，t分别代以t，t?2?（其中2?为时间步长），并按抛物线法则计算式中的卷积，有下式：
x(t?2?)?e?2?????[sin2???cos2??]?x(t)?
1e?2??????(t)?sin2???x??2???esin2???r(t)3?
4??????esin???r(t??t)3? （18）
22??????2????(t?2?)?ex[????sin2??]?x(t)???
???(t)?e[cos2???sin2??]?x?
??2??????e(cos2???sin2??)?r(t)3?
?2????4????????e(cos???sin??)?r(t??)??r(t?2?)3?3
（19）
以上两式以矩阵表示为：
?x(t)????x(t)?a15???r(t)??a25???r(t??)???r(t?2?)???
?x(t?2?)??a11?????a?x(t?2?)??21a12a22a13a23a14a24（20） 记
e1?e????，e2?e，
21c0???2?，c1?cos??，s1?sin??，s2?2s1c1，c2?2c1?1
（21）
则上述A矩阵元素为
a11?e2(c2?c0s2)，a12?e2s2/?，a13?a12a25，a14?4e1s1a25/?，a15?0 a21???2a12，a22?e2(c2?c0s2)，a23?a22a25a24?4e1(c1?c0s1)a25，a25??/3 应用上述递推公式， 以前一时刻来求后一时刻的结果。计算不重复。 当后在时域中的步进求解只是一些简单的数组相乘。计算速度很快。
aij求出后，以
二 Newmark类方法
下面讨论直接积分法
?(tn)?x(tn?1)?x(tn)??t?x....121?(tn)??t3???(tn)?O(?t4)?t?xx26 （1）
...12x(tn?1)?x(tn)??t?x(tn)??tx(tn)?O(?t3)2 （2）
x(tn?1)?x(tn)??t?x(tn)?O(?t2)将（3）代入（1），（2）得：
.......... （3）
x(tn?1)?x(tn)x(tn)??O(?t)?t.....
?t2..?t2..x(tn?1)?x(tn)??tx(tn)?x(tn)?x(tn?1)?O(?t4)36
?t..?t..x(tn?1)?x(tn)?x(tn)?x(tn?1)?O(?t3)22
..1、可以直接略去高阶项
2、用变权来调节
?n?[(xn?1?xn??tx1?n????n?1]?t2??)?xx2