有限元动力学分析方程及解法
动力分析中平衡方程组的解法
1前言 描述结构动力学特征的基本力学变量和方程与静力问题类似，但所有的变量都是时间的函数。
基本变量
三大类变量(,)i u t ξ、(,)ij t εξ和(,)ij t σξ是坐标位置(,,)x y z ξ和时间t 的函数，一般将其记为()()()i ij ij u t t t εσ。
基本方程
(1) 平衡方程
利用达朗贝尔原理将惯性力和阻尼力等效到静力平衡方程中，有
,()()()()0ij j i i i t b t u t u t σρν+--= (1)
其中ρ为密度，ν为阻尼系数。
(2) 几何方程
）
,,1()(()())2
ij i j j i t u t u t ε=+ (2)
(3) 物理方程 ()()ij ijkl kl t D t σε= (3)
其中ijkl D 为弹性系数矩阵。
(4) 边界条件
位移边界条件()BC u 为，
()()i i u t u t = 在u S 上 (4)
力的边界条件()BC p 为，
()()ij j i t n p t σ= 在p S 上 (5)
初始条件
0(,0)()i i u t u ξξ== (6)
{ 0(,0)()i i u t u ξξ== (7)
虚功原理
基于上述基本方程，可以写出平衡方程及力边界条件下的等效积分形式，
,()()0p
ij j i i i ij j i S u u b u d n p dA δσρνδσΩ∏=---+Ω+-=?? (8)
对该方程右端第一项进行分部积分，并应用高斯-格林公式，整理得，
()()0p
ijkl ij kl i i i i i i i i S D u u u u d b u d p u dA εδερδνδδδΩΩ-++Ω-Ω+=??? (9) 有限元分析列式
单元的节点位移列阵为，
111222()[(),(),(),(),(),()(),(),()]e t k k k U t u t v t w t u t v t w t u t v t w t = (10)
单元内的插值函数为，
(,)()()e t u t N U t ξξ= (11)
%
其中()N ξ为单元的形状函数矩阵，与相应的静力问题单元的形状函数矩阵完全相同，ξ为单元中的几何位置坐标。
基于上面的几何方程和物理方程及(11)式，将相关的物理量表达为节点位移的关系，有，
(,)[](,)[]()()()()e e t t t u t N U t B U t εξξξξ=?=?= (12)
(,)()()()()e e t t t D DB U t S U t σξεξξ=== (13)
(,)()()e t u t N U t ξξ= (14)
(,)()()e t u t N U t ξξ= (15)
将(12)-(15)供稿到虚功方程(9)中，有，
[()()()()]()0e e e e e e e
T e t t t t t M U t C U t K U t R t U t δδ∏=++-= (16)
由于()e t U t δ具有任意性，消去该项并简写有，
e e e e e
t t t t U C U KU R ++= (17)
其中，
】
e e T M N
Nd ρΩ=
Ω? (18) e
e T C N Nd νΩ=Ω? (19) e
e T K B DBd Ω=Ω? (20)
e M 为单元质量矩阵，e C 为单元阻尼矩阵，e K 为单元刚度矩阵。同样，将单元的各个矩阵进行组装，可形成系统的整体有限元方程，即，
MU CU KU R ++= (21)
其中M ，C 和K 分别是系统的质量、阻尼和刚度矩阵，R 是外荷载向量，U ,U
和U 分别是有限元分割体的加速度、速度和位移向量。方程(21)是通过考虑在时刻t 的静力平衡而推导出来的。
对静力或动力分析的选择(即在分析中是考虑或忽略与速度及加速度有关的力)，一般取决于工程上的判断，其目的在于减少所需要的分析工作量。但是，应该认识到，一个静力分析的假定，应该有理由说明它是正确的，否则，分析的结果就是无意义的。确实，在非线性分析中，采用忽略惯性力和阻尼力的假定，可能严重到难以求得甚至无法求得解答。
在数学上，方程(21)是一个二阶线性微分方程组，原则上可用求解常系数微分方程组的标准过程来求得方程组的解。但是，如果矩阵的阶数很高，则采用求解一般微分方程组的过程可能要付出很高的费用，除非特别利用系数矩阵K ，C 和M 的特殊性质。因此，在实用的有限元分析中，主要对几种有效的方法感兴趣，下面将集中介绍这几种方法。我们所考虑的基本过程，可分为两种求解方法：直接积分法和振型叠加法。初看起来，这两种方法似乎完全不同，但事实上它们有着密切的关系，至于选择这种或那种方法，只取决于它们的数值效果。
2直接积分法
在直接积分中对方程(21)是逐步地进行数值积分的，“直接”的意思是，进行数值积分前没有进行把方程变为另一种形式的变换。实质上，直接积分是基于下面的两个想法，第一个想法是只在相隔t ?的一些离散的时间区间上而不是试图在任一时刻t 上满足方程(21)即包含有惯性力和阻尼力作用的(静力)平衡是在求解区间上的一些离散时刻点上获得的。因此，似乎在静力分析中使用过的所有求解方法，在直接积分法中或许也能有效地使用；第二个想法是假定位移、速
度和加速度在每一时间区间t ?内变化。
下面假设分别用0
00U ,U ,U 来表示初始时刻)t (0=的位移、速度和加速度向量为已知，要求出方程(21)从0=t 到T t =的解。在求解时，把时间全程T 划分为几个相等的时间区间t ?(即n /T t =?)，所用的积分格式是在时刻t ,?0，
T ,,t t ,t ,,t ?+?2上确定方程的近似解。
由于计算下一个时刻的解的算法要考虑到前面各个时刻的解，因此假定在时刻t ,,t , ?0的解为已知，来推导出求时刻t t ?+的解的算法。计算时刻t t ?+的解对于计算自此以后t ?的时刻上的解是有代表意义的，这样就可建立用来计算在所有离散时间点上解的一般算法。
（a ） %
（b ） 中心差分法
若把式(21)的平衡关系看作是一个常系数常微分方程组，便可以用任一有限差分表达式通过位移来近似表示加速度和速度。因此，在理论上，许多不同的有限差分表达式均可使用。但是，我们要求求解格式必须是有效的，这样便只需考虑少数几种计算格式。对某些问题求解是非常有效的一个过程是中心差分法，这个方法假定
{}{}
t t t t t t t t t t t U U t U U U U t U ?+?-?+?-+-?=+-?=21212 (22)
将式(22)代入t 时刻的式(21)，可得
t t t t t t U C t M t
U M t K R U C t M t ?-?+??? ???-?-??? ???--=??? ???+?2112211222 (23)
从式(23)我们可以求出t t U ?+。应该注意，t t U ?+的解是基于利用在时刻t 的平衡条件。因此，该积分过程称为显式积分方法，且这样的积分格式在逐步解法中不需要对(有效)刚度矩阵进行分解。另一方面，以后所考虑的Houbolt ，Wilson θ及Newmark 方法，要利用在t t ?+上的平衡条件，因而称为隐式积分方法。
另外还应注意到，应用中心差分法时，t t U ?+的计算包含有t U 和t t U ?-，因此，
计算在时刻t ?的解，必需用一个具体的起始过程。由于0
00U ,U ,U 都是已知的，由关系式(22)可求
02002
U t U t U U t ?+?-=?- (24) 具体计算步骤为
A ． 初始计算
1．·
2． 形成刚度矩阵K 、质量矩阵M 和阻尼矩阵C 。
3．计算初始值0
00U ,U ,U 。 4．选取时间步长t ?，要求cr t t ?≤?（临界值）。
5．计算系数201t
a ?=，t a ?=211，022a a =，231a a =。 6．计算0300U a U t U U t +?-=?-。
7．形成有效质量矩阵C a M a M ?1
0+=。 8．对M ?作三角分解：T LDL M ?=
B ． 每一时间步长内的计算
1．计算在时刻t 的有效荷载：
()()t
t t t t U C a M a U M a K R R ??-----=102。 2．计算时刻t t ?+的位移：t
t t T R ?U LDL =?+。 3．.
4． 必要时，按照式（）计算时刻t 速度和加速度。
假设所考虑的系统没有物理阻尼，即C 是零矩阵，在这种情形下式(23)可简化为
t t t R ?MU t
=??+21 (25) 其中
()()t
t t t t U C a M a U M a K R R ??-----=102 因此，如果质量矩阵是对角形的，则解方程组（）时就不需要进行矩阵的分解，
即只需进行矩阵相乘便可求得右端项的有效荷载向量t
R ?，从而利用 ???
? ???=?+ii )i (t )i (t t m t R ?U
2 (26)
可得出位移向量的各个分量，其中)i (t t U ?+和)i (t R ?分别表示向量t t U ?+和t
R ?的第i 个分量，而ii m 是质量矩阵的第i 个对角线元素，并且假定0>ii m 。
如果对总刚度矩阵和质量矩阵都不需进行三角分解，也就不必形成总体的K 和M 。此时，求解式(23)可以在单元一级来解决，然后将每个单元的结果累加即可，即
)U U (M t U K R R ?t i
t t i t i i t t 212-?--=∑∑?- (27) 使用式(26)和(27)形式的中心差分法的优点是很明显的，因为它不需要计算总刚度矩阵和总质量矩阵，求解过程基本上是在单元一级上进行，所需要的内存比较少。如果所有相继的单元刚度矩阵和质量矩阵均相同，则该方法就显得更有效，因为这时只需计算或从后备存贮器上连续读出对应于系统中第一个单元的矩阵。
|
至于中心差分法的缺点，必需承认，该过程的效果与对角形质量矩阵中采用和忽略通常依赖于速度的阻尼力有关，若只包合一个对角形阻尼矩阵，则仍然可保持在单元一级上进行求解的优点。从实用上看，只能用于对角形质量矩阵的这个缺点通常是不很严重的，因为可以采用足够精细的有限元离散化来使解有良好的精度。
使用中心差分格式的另一个十分重要的考虑是，该积分方法要求时间步长t ?小于一个临界值cr t ?，可由整个单元分割体的刚度和质量的性质来算出cr t ?。更准确地说，要得到一个有效的解必须
πn
cr T t t =?≤? (28)
其中n T 是分析物体的最小周期，n 是单元系统的阶。
要求使用的时间步长t ?小于临界时间步长cr t ?的差分格式，例如中心差分法，称为条件稳定的。如果使用一个大于cr t ?的时间步长，则积分是不稳定的，这意味着由数值积分或在计算机上的舍入所导致的误差都会增大，并且在许多情形下会使响应的计算失去意义，积分稳定性的概念是非常重要的。
上面我们讨论了中心差分方法的主要缺点：其格式是条件稳定的，许多其它积分方法也会是条件稳定的。由于在结构分析中条件稳定方法的有效使用会受到限制，因此在下面几节我们来考虑一些通常使用的无条件稳定的积分格式。无条
件稳定的积分格式的有效性可从积分能取得较好精度这一事实显示出来，时间步长t ?可以不按式(28)的要求束选取，在许多情形下t ?可以大于式(28)所允许的t ?几个量级。所有下面讨论的积分方法都是隐式的，即在求解时要求对刚度矩阵K ，或者说得更确切些是对有效刚度矩阵进行三角分解。
（c ） Houbolt 法
Houbolt 积分格式与前面讨论的中心差分方法是有联系的，因为它也是根据标准的有限差分表示式用位移的各个分量来近似地表示加速度和速度的各个分量。Houbolt 积分格式采用了下面的有限差分展开式
{}{}t t t t t t t t t t t t t t t t t t U U U U t
U U U U U t U ?-?-?+?+?-?-?+?+-+-?=-+-?=222291811614521 (29) 为了得到在时刻t t ?+的解，现在我们考虑在时刻t t ?+的式(21)(而不是象中心差分法那样考虑在时刻t 的情形)，它给出
t t t t U C t M t R U K C t M t ??
? ???+?+=??? ??+?+??+?+35611222 -
t t t t U C t M t
U C t M t ?-?-??? ???+?+??? ???+?-222311234 (30) 利用上式求解t t U ?+，需要先知道t U ，t t U ?-和t t U ?-2，虽然知道0
00U ,U ,U 对开始Houbolt 积分格式是有用的，但用其它方法计算t t U ,U ??2会更精确，即需要利用特殊的起始方法。
用来积分式(21)以得出t t U ,U ??2，是用一个不同的积分格式并尽可能是用以t ?的几分之一为时间步长，例如用中心差分格式。
Houbolt 法与中心差分法的一个基本差别在于，刚度矩阵K 是作为所求位移t t U ?+的因子出现的。t t KU ?+项的出现是因为在式(21)中取时刻t t ?+的平衡而不像中心差分法取时刻t 的平衡，因此Houbolt 方法是一个隐式积分格式，而中心差分方法则是一个显式过程。关于积分所用的时间步长t ?不存在临界的时间步长限制，一般t ?可以选得比式(28)所给的中心差分法的步长大一些。
一个值得注意的特点是，如果忽略质量和阻尼的影响，基于Houbolt 方法的逐步积分格式可直接简化为静力分析，但中心差分解法则不能这样做。换言之，若0=C 和0=M ，Houbolt 求解方法就得出与时间有关的荷裁作用下的静力解。
（d ） Wilson θ法
Wilson θ法实质上是线性加速度法的推广，线性加速度法假定加速度从时刻t 到时刻t t ?+为线性变化。Wilson θ方程则假定加速度从时刻t 到时刻t t ?+θ为线性变化，其中01.≥θ，当01.=θ时，这个方法就简化为线性加速度格式。要达到无条件稳定，则必须用371.≥θ，通常用401.=θ。
令τ表示时间的增量，其中t ?≤≤θτ0，于是对从t 到t t ?+θ的时间区间，假定
()t t t t t U U t
U U -?+=?++θτθτ (31) 对上式积分得
()()t t t t t t t t t t t t t U U t U U U U U U t U U U -?+++=-?++=?++?++θτ
θτ
θτττθττ6212322 (32) |
当t ?=θτ时
()()t t t t t t t t t t t t t U U t U t U U U U t U U 26
222+?+?+=+?+=?+?+?+?+θθθθθθθ (33) 这样，就可以求出
()()
t t t t t t t t t t t t t t U t U U U t U U U t U U t U 22326622?---?=-?--?=
?+?+?+?+θθθθθθθθ (34) 这样，要得到在时刻t t ?+的位移、速度和加速度的解，就只需考虑在时刻t t ?+θ的平衡方程(21)。然而，因为假定加速度为线性变化，故所用的投影荷载向量是线性变化的，即
()t t t t t t R R R R -+=?+?+θθ (35)
具体计算步骤为
A ． 初始计算
1．形成刚度矩阵K 、质量矩阵M 和阻尼矩阵C 。
2．计算初始值0
00U ,U ,U 。
3．选取时间步长t ?。计算系数()206t a ?=
θ，t a ?=θ31，122a a =，23t a ?=θ，θ04a a =，θ25a a -=，θ316-=a ，2
7t a ?=，628t a ?= 4．；
5． 形成有效刚度矩阵C a M a K K ?1
0++=。 6．对K
?作三角分解：T LDL K ?= B ． 每一时间步长内的计算
1．计算在时刻t t ?+θ的有效荷载：
()()()t
t t t t t t t t t t t U a U U a C U U a U a M R R R R ? 312022++++++-+=?+?+θθ。 2．计算时刻t t ?+θ的位移：t t t t T R ?U LDL ?+?+=θθ。
3．计算时刻t t ?+的位移、速度和加速度：
()t t t t t t t U a U a U U a U 654++-=?+?+θ
()t t t t t t U U a U U ++=?+?+7
()t t t t t t t U U a U t U U 28++?+=?+?+
正如前面所指出的，Wilson θ法也是一个隐式积分方法，因为刚度矩阵K 是未知位移向量的系数矩阵。还可以注意到，该方法不需要特别的初始过程，因为在时刻t t ?+的位移、速度和加速度只是利用在时刻t 的相同的量来表示。
（e ） ￥
（f ） Newmark 法
Newmark 积分格式也可以认为是线性加速度法的推广，它所用的假定如下 ()[]2211t U U t U U U t
U U U U t t t t t t t t t t t t t ???
????+??? ??-+?+=?+-+=?+?+?+?+ ααδδ (36) 其中α和δ是参数，根据积分的精度和稳定性的要求来确定这两个参数。当21
=
δ和6
1=α时，式（）相应于线性加速度法（它也可由Wilson θ法取1=θ得到）。Newmark 最初提出以‘恒定—平均—加速度法’作为无条件稳定的格式，在这种
情形下，21=δ，4
1=α。 要得到在时刻t t ?+的位移、速度和加速度的解，就只需考虑在时刻t t ?+的平衡方程(21)，具体步骤如下：
A ． 初始计算
1．形成刚度矩阵K 、质量矩阵M 和阻尼矩阵C 。
2．计算初始值0
00U ,U ,U 。 3．选取时间步长t ?，参数δ和α。计算积分系数：
50.≥δ，()250250δα+≥..，2
01t a ?=α，t a ?=αδ1，t a ?=α12，1213-=αa ，14-=αδa ，??
? ??-?=225αδt a ，()δ-?=16t a ，t a ?=δ7 4．形成有效刚度矩阵C a M a K K ?1
0++=。 5．`
6． 对K
?作三角分解：T LDL K ?= B ． 每一时间步长内的计算
1．计算在时刻t t ?+的有效荷载：
()()t t t t t t t t t t U a U a U a C U a U a U a M R R ? 541320++++++=?+?+。
2．计算时刻t t ?+的位移：t
t t t T R ?U LDL ?+?+=。 3．计算时刻t t ?+的速度和加速度：
()t t t t t t t U a U a U U a U 320---=?+?+
t t t t t t U a U a U U ?+?+++= 76
应该注意Newmark 法和Wilson θ法在计算机上执行时的密切关系，利用这个关系就有可能在一个简单的计算机程序中方便地使用这两个积分格式。 11.1 振型叠加法
在直接积分法中，质量矩阵是对角形的且无阻尼，则对一个时间步长的运算次数稍多于k nm 2，其中n 和k m 分别是所考虑的刚度矩阵的阶和半带宽。在中心
差分法中，刚度矩阵与位移向量相乘需要k nm 2次运算；而在Houbolt 法、Wilson θ法和Newmark 法中，则在每个时间步长求解方程组时约需要k nm 2次运算。有效刚度矩阵的初始三角分解还要求一些附加的运算，并且，如果在分析中使用一致质量矩阵或考虑阻尼矩阵，对其中任一情形，每一时间步长所需要的附加运算次数正比于k nm 。因此，若略去初始计算的运算次数，则整个积分所要求的总运算次数约为s nm k α，其中α与所用的矩阵的性质有关，2≥α，而s 是时间步长的步数。
^
上面的考虑说明，直接积分所需要的运算次数直接正比于分析中的时间步数。因此，一般说来，当要求较短时间（即n 个时间步效）的响应时，可以预料，使用直接积分法是有效的。但是，如果积分必须对许多时间步数进行，则先把平衡方程(21)变换，使之能以较少的费用进行逐步求解就可能会更有效。具体地说：由于所需要的运算次数直接正比于刚度矩阵的半带宽k m ，因而k m 的减小会按比例地降低逐步解法的费用。
平衡方程组(21)是将有限元插值函数应用在虚功方程的计算中得到的，而最终的矩阵K ，M 和C 的带宽由有限元节点的编号束决定。因此，有限元网格的决定了系统矩阵的阶数和带宽。虽然我们可以通过重新编排节点号达到减小系统矩阵带宽的目的，但是，用这种方式能够得到的最小带宽是有限度的，因而有必要探索另外的途径，这就是振型叠加法。
利用下面的变换，对有限元分析的位移进行变换，即
)t (PX )t (U = (37)
其中P 是个待确定的方阵，)t (X 是与时间有关的n 阶向量，)t (X 的各个分量称为广义位移。将上式代入式(21)，并左乘T P 得到
)t (R ~)t (X K ~)t (X C ~)t (X M ~=++ (38)
其中
MP P M ~T =，CP P C ~T =，KP P K ~T =
实际上，这个变换就是将位移的表达式变为
)t (PX N )t ,z ,y ,x (U )m ()m (= (39)
变换的目的是要得到新的系统刚度、质量和阻尼矩阵M ~,K ~和C ~，使它们的
带宽比原来系统矩阵的小，同时也应适当地选取变换矩阵P 。此外应注意，为了使式(37)所表示的任何向量U 和X 之间有唯一的关系，P 必须是非奇异的。 '
在理论上可以有许多不同的变换矩阵P 都可以减小原来系统矩阵的带宽，但
在实用上，一个有效的变换矩阵是由无阻尼自由振动平衡方程0=+KU U
M 的解来建立的，该方程的可以假定解为
()0t t sin U -=ωφ (40)
其中φ是n 阶向量，t 是时间变量，0t 是时间常数，而ω是表示向量φ的振动频率的常量。把式(40)代入无阻尼自由振动平衡方程，得到广义特征值问题
φωφM K 2= (41)
从中可以确定φ和ω，特征问题(41)式有n 个特征解()121φω,，()222
φω,，……，()
n n ,φω2
，其中特征向量与M 是正交的，即 ??
?≠====)
j i (,)j i (,M j T i 01φφ (42) 而 222210n ωωω≤≤≤≤ (43)
向量i φ称为第i 阶振型向量，i ω是相应的振动频率。应该强调指出，n 个位移解公)n ,,,i )(t t (sin i i 210=-ωφ中的任何一个均满足振动平衡方程。
定义一个矩阵Φ，它的各列为特征向量i φ，并定义一个矩阵2Ω，其对角线上的元素为特征值2i ω，即
[]n φφφ 21=Φ，??????
????????=Ω222212n ωωω (44) ，
则式(41)可以写成
2ΦΩ=ΦM K (45)
由于向量M 是正交的，于是
I M ,K T T =ΦΦΩ=ΦΦ2 (46)
很显然，矩阵Φ就是一个合适的P 。利用
)t (X )t (U Φ= (47)
可得到对应于振型广义位移的平衡方程
)t (R )t (X )t (X C )t (X
T Φ=Ω+ΦΦ+2 (48) 根据正交性，可得初始条件
00MU X T Φ=，0
0U M X T Φ= (49) 方程(48)表明，若在分析中不包含阻尼矩阵，当利用有限元系统的自由振动的振型向量组成变换矩阵P 时，则有限元平衡方程组是不耦合的，由于在许多情形下阻尼矩阵不能明显确定而只能近似地考虑阻尼的影响，采用下述阻尼矩阵是合理的：它包含所有要求的影响同时又能有效地求出平衡方程的解。在许多分析中是完全忽略阻尼影响的。
（a ） ：
（b ） 忽略阻尼的分析
若在分析中不考虑与速度有关的阻尼影响，方程（）就简化为
)t (R )t (X )t (X T Φ=Ω+2 (50) 亦即下面的n 个独立的方程
n ,,,i )t (R )t (r )t (r x )t (x T i i i i i i 212=?
??==+φω (51) 其中应注意到，在(51)中的第i 个典型方程是一个具有单位质量和刚度2i ω的单自由度系统的平衡方程，这个系统运动时的初始条件可从式(49)得出
000U M )(x ,MU )(x T i i T i i φφ== (52) (51)中的每一项积分，都可以用上节的直接积分法求出，也可以利用杜哈密（Duhamel ）积分求得
?++-=t i i i i i i i i t cos t sin d )t (sin )(r )t (x 01
ωβωαττωτω (53)
式中i i ,βα是由初始条件确定的。一般上式的杜哈密积分要采用数值积分来计算。
对于完整的响应，必须算出式(51)中的所有n 个方程的解，然后将每个振型的响应叠加就得到有限元的节点位移，即
。
∑=i
i i )t (x )t (U φ (54)
因此，概括地讲，用振型叠加的响应分析要求：首先求解方程(41)的特征值和特征向量；然后求解不耦合的平衡方程组(51)；最后，如式(54)所示，将每一个特征向量的响应进行叠加。在分析中，特征向量是有限元分割体自由振动的振型。如前面讨论过的，选择振型叠加分析或用上节所述的直接积分法，仅仅取决于其数值效果，而这两种过程所得的解在所用的时间积分格式的数值误差和计算机的舍入误差两方面部是相同的。
对于许多类型的实际荷载条件来说，为了得到系统真实响应的一个好的近似解，仅需要考虑全部不耦合方程中的一部分就够了。通常仅需用到头p 个平衡方程，即为了得到一个满意的近似解，我们在分析中需要考虑式(51)中有关p ,,,i 21=的方程，其中n p <<。这意味着仅需对(41)式解出最小的p 个特征值及相应的特征向量，同时，在式(54)中仅对头p 个振型进行叠加。
事实上，在振型叠加分析中可以仅需考虑少数几个振型，振型叠加过程比直接积分有效得多。但是，由此也可知，振型叠加的效果与分析中必须考虑的振型数有关，通常所考虑的结构和荷载的空间分布以及频率成份决定着所用振型的个数。对于地震荷载，有些情形只需考虑十个最低的振型，虽然系统的阶可能大于1000。另一方面，对于冲击或振动荷载，一般要包含更多的振型，可能要n p 3
2=。 当考虑振型叠加分析中所应包含的振型的个数的选取问题时，应该记住的是，找出的是动力方程(21)的一个近似解。因此，如果所考虑的振型的个数不足，方程(21)的解就不可能足够精确。实际上这意味着计算出的近似响应，并不满足包含惯性力的平衡。以p U 表示考虑p 个振型时由振型叠加而算得的响应，则可以用下面的误差估计p ε，来估计任意时刻t 的分析精度。 []22)t (R )t (KU )t (U
M )t (R )t (p p p +-= ε (55)
如果己得到系统平衡方程(21)一个满意的近似解，则)t (p ε在任何时刻t 都是很小的。但必需注意，)t (U p 的求得，一定要基于所考虑的p 个振型中每一个振型
的响应的精确计算，因为该方法的唯一误差是由于在分析中没有包含足够的振型所致。
应该注意，由上式算出的误差度量决定着包括惯性力的平衡被满足到什么程度，同时它又是惯性和弹性节点力与节点荷载三者平衡的度量.另一方面，因为不是所有振型向量都被用到，所以我们可以说，p ε是没有包括在振型叠加分析中的那些外荷载向量部分的度量，由此可以看出，在直接积分分析中，p ε在时刻 ,t ,t ,??20总是等于零的。总起来说，假定已精确地解出(51)中互不耦合的方程，用n p <的振型叠加分析中的误差是由于没有使用足够的振型所致，而在直接积分分析时，误差的产生是由于使用了过大的时间步长所致。
还要指出一个更为重要的情况，到目前为止仅讨论了有限元系统平衡方程
(21)的精确解。但是，我们真正的目的是要求得所考虑的结构的真实的精确响应的一个良好的近似解。事实上，有限元分析能够求得真实结构的“精确”频率的上界。通常，有限元分析非常好地逼近于最低的精确频率，而可以预料，在逼近高频和相应振型方面就几乎没有精度可言。振型叠加过程在分析中不必包括对应于有限元系统中可能不精确的高频响应，与直接积分相比，这是该方法的特有优点。因而，假定在有限元分析中准确地算出了所有重要的频率，则就不必计算系统中高振型的响应，也不会严重影响解的精度。
（c ） 考虑阻尼的分析
以特征向量)n ,,,i (i 21=φ为基的有限元系统平衡方程的一般形式由式(48)给出。该式表明，假如略去阻尼的影响，则平衡方程就互不耦合，即可以对每一方程逐个地进行时间积分。对不能略去阻尼影响的系统的分析，我们仍然希望基本上能使用相同的计算过程去处理互不耦合的平衡方程(48)。通常，阻尼矩阵C 不能象形成单元分割体的质量矩阵和刚度矩阵那样，由单无阻尼矩阵形成，而它的用处在于近似表示出在系统响应期间的全部能量的消耗。如果能够假定阻尼是成比例的，则振型叠加分析特别有效，该假定为
'
ij i i j T i C δξωφφ2= (56)
其中i ξ是振型阻尼参数，ij δ是Kronecker δ符号。因此，利用上式，假定特征向量)n ,,,i (i 21=φ也是C 正交的，则方程(48)简化为下面形式的n 个方程：
)t (r x )t (x )t (x
i i i i i i i =++22ωξω (57)
同样，这个方程可以用前节的直接积分计算，也可以用Duhamel 积分
[]?++-=---t i i i i t i )t (i i i t cos t sin e d )t (sin e )(r )t (x i i i i 01
ωβωαττωτωωξτξω (58) 其中21i i i ξωω-=，其它同无阻尼一样。
若阻尼满足关系(56)式，在振型叠加分析中就易于考虑阻尼的影响。然而，在真实的阻尼比)p ,,,i (i 21=ξ为已知，则需要用显式算出矩阵C ，把它代人式
(56)中即可得出确定的阻尼比。如果2=p ，可假定Rayleigh 阻尼为如下形式
K M C βα+= (59)
式中α和β是常数。
例如：假设一个多自由度系统，21=ω，32=ω且在该两个振型中我们分别要求有2%和10%的临界阻尼，即要求0201.=ξ和1002.=ξ，试确定Rayleigh 阻尼的常数。
根据Rayleigh 阻尼，有
!
K M C βα+=
而利用关系（），可得
i i i i T i )K M (ξωβωαφβαφ22=+=+
将两队阻尼比和相应的频率代入
60
090804..=+=+βαβα 解得3360.-=α，1040.=β。这样阻尼矩阵为
K .M .C 10403360+-=
由于已给出阻尼矩阵，如果采用Rayleigh 阻尼矩阵，我们就可以确定任何值i ω的阻尼比
i
i i ωβωαξ22+=。 在实际的分析中，很可能已知两个以上频率的阻尼比，在这种情形下，可用两个平均值1ξ和2ξ来计算α和β。
例如：假设对一个多自由度系统的近似阻尼指定如下：
019
0150070030021401000400300205432154321==========ωωωωωξξξξξ,,,,.,.,.,.,. 试选择合适的Rayleigh 阻尼常数。
由于只要两对阻尼比与频率，根据已知数据，选择
1712040302211====ωξωξ,.;,.
代入可得
08
428924016..=+=+βαβα 解得014980.=α，014050.=β。这样，阻尼矩阵为K .M .C 014050014980+=，计算时的阻尼比由下式求出
i
i i ..ωωξ20140500149802+= 对于用两个以上的阻尼比，可用柯西（Caughey ）级数的方式，来建立阻尼矩阵，具体的步骤可参考相关有限元书籍。但是，大多数采用直接积分的实际分析仍采用Rayleigh 阻尼的假定。Rayleigh 阻尼的一个缺点是高振型衰减大于低振型，因为Rayleigh 常数是按低振型选定的。
实际上，对一个具体结构的分析，其合适的Rayleigh 系数常常可以利用一个典型的同类结构的阻尼特性的已有数据来选取，即，把同样的α和β近似地用在同类的结构分析中。Rayleigh 系数的大小，在很大程度上是由结构材料的能量消耗特性来决定的。
在以上的讨论中，不论在振型叠加分析或者在直接积分过程中，我们都假定可以利用比例阻尼来适当地表示结构的阻尼特性。在多数分析中，比例阻尼的假定（即式(56)满足）是满足要求的。然而，在那些材料性质有很大变化的结构的分机中，就可能需要使用不成比例的阻尼，例如，在基础与结构相互影响的问题中，可以观察到基础的阻尼显著大于地面上结构的阻尼。在这种情形下，在形成结构的阻尼矩阵时，按不同的结构部分规定不同的Rayleigh 系数是合理的，这就导致一个不满足关系式(56)的阻尼矩阵。当把阻尼集中在特定的自由度上时(如在结构的支点上)就遇到另一种阻尼不成比例的情形。
具有不成比例阻尼的有限元平衡方程组，可不加修改地利用上节的直接积分算法进行求解，因为阻尼矩阵的性质并不影响到求解过程的推导。另一方面，对用无阻尼自由振动振型作为基向量的振型叠加分析而言，我们看到，式(48)中ΦΦC T 在阻尼不成比例的情形下是一个满矩阵，换言之，以振型向量为基的平衡
方程组不再是互不耦合的。