第五章 结的强迫振动响应分析

第五章 结构的强迫振动响应分析
§5.1 概述
如果结构已经用有限元方法进行了离散化，当一个结构系统受到外激励作用时，其响应就是一个多自由度系统的强迫振动问题的解。

求解多自由度系统强迫振动响应的方法之一就是直接积分法。

考虑到实际结构的高维数（自由度数很大）而给求解带来的困难，往往在实际求解中采用模态叠加法。

直接积分法和模态叠加法这两种方法都可以得到具有相当精度的振动响应解，并且各有其特点。

§5.2 求解强迫振动响应的直接积分法
对动力学基本方程
)}({}]{[}]{[}]{[t P U K U C U M =++&&& （5－1）
进行直接积分，其含义是指在对方程进行积分之前，不对其进行任何形式的变换，在积分中，实际上是按时间步长逐步积分的。

这样做的实质是基于如下考虑：
（1） 只在相隔t ∆的一些离散时间区间上、而不是在整个时间区间上的任一个
时刻t 上满足方程，即平衡是在求解区间上的一些离散时刻上获得的。

（2） 假定位移、速度、加速度在每一个时间区间t ∆内按一定规律变化，也正
是采用不同的变化形式，决定了各种直接积分解的精度、稳定性和求解速度。

首先，设}{}{}{0
00U U U &&&表示初始时刻（0=t ）的位移、速度和加速度为已知向量，要求出从0=t 到T t =的解，则把时间段T 均分为n 个间隔n T t /=∆，所用的积分是在T t t Λ,2,∆∆上求方程的近似解。

即要在t t t Λ,2,∆∆的解已知的情况下，求解t t ∆+时刻的解。

【中心差分法】
若基本方程式的平衡关系作为一个常系数微分方程组，则可以用任一种差分格式通过位移来表示速度和加速度。

通常采用中心差分格式，这是一个行之有效的求解微分方程的格式。

}){}({21}{}){}{2}({1
}{2
t t t t t
t t t t t t U U t U U U U t U ∆∆∆∆∆∆-++--=+-=
&&& （5－2）
假定}{t U 及前一时刻的位移}{t t U ∆-已经求得，则将}{t U &&}{t U &代入方程
（5－1）得到：
}]){[21][1(}]){[2]([}{}]){[21][1(
22
2t t t t t t U C t
M t U M t K P U C t M t ∆∆∆∆∆∆∆-+----=+ （5－3）
由此式求出}{t t U ∆+
上述格式是一个显式格式。

具体计算时，还有一个步进递推格式启动的问题，即}{}{}{0
00U U U &&&已知时，}{t U ∆-的求解问题。

由}{t U &&}{t U &的差分表达式，可求出：
}{2
}{}{}{0
2
00U t U t U U t &&&∆∆∆+-=- （5－4） 中心差分法的具体步骤为：
1． 用有限元素法形成结构的刚度矩阵][K 、质量矩阵][M 、阻尼矩阵][C
2． 计算初始值}{}{}{0
00U U U &&& 3． 选择步长t ∆，并计算积分常数
2
3021201
,2,21,1a a a a t a t a ====
∆∆ （5－5） 4． 计算
}{}{}{}{0
300U a U t U U t &&&+-=-∆∆ （5－6） 5． 形成
][][]ˆ[1
0C a M a M += （5－7） 6． 分解
T L D L M
]][][[]ˆ[= （5－8） 对每一步长，进行如下计算： 1． 求t 时刻的有效载荷
}]){[][(}]){[]([}{}ˆ{102
t t t t t U C a M a U M a K P P ∆-----= （5－9）
2． 求解在时刻t t ∆+的位移
}ˆ{}{]][][[t
t t T P U L D L =+∆ （5－10） 3． 如果需要，计算t 时刻的加速度和速度：
})
{}{(}{}){}{2}({
}{1
0t
t t
t t
t t t t t t U U a U U U U a U ∆∆∆∆+-+-+-=+-=&&& （5－11）
中心差分格式使用中，一个重要的问题是步长必须小于临界步长n
T t n
cr =
∆，以保证步进递推的数值稳定性，这里，n 为系统的阶数，
n T 为系统最小自然周期，即：
n
T t t n
cr =
≤∆∆ （5－12） 因此，中心差分格式是条件稳定的。

中心差分法作为显式算法的优点是，当质量阵为对角阵，阻尼阵也可以对角化时，可以避免矩阵求逆运算，而矩阵的分解运算非常简单，特别在进行非线性系统的响应分析时，由于需要在每个时间增量步修改刚度矩阵，采用中心差分法可以避免每一增量步对刚度矩阵的分解。

【Houbolt 方法】
Houbolt 方法也是一种差分方法，它是基于拉格朗日插值公式的步进方法，该方法利用向后差分，由位移导出速度和加速度的多步隐式公式。

Houbolt 方法得到的计算结果比较光滑。

其差分的格式为：
}){2}{9}{18}{11(61}{}){}{4}{5}{2(1
}{222t 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 ∆∆∆∆∆∆∆∆--+--+-+-=-+-=
&&& （5－13）
考虑t t ∆+时刻的平衡方程
}{}]{[}]{[}]{[t
t t t t t t t P U K U C U M ∆∆∆∆++++=++&&& （5－14） 从而有：
}]){[31][1})]){[23
][4(}]){[3][5(}{}]){[][161][2(
22222t t t t t t t t t U C t
M t U C t
M t U C t M t P U K C t
M t ∆-∆-∆+∆+∆+∆+∆+∆-∆+∆+=+∆+∆ （5－15）
显然，要求解}{t t U ∆+必须知道}{},{},{2t t t t t U U U ∆∆--
在使用Houbolt 方法时，不是用此格式求初始两个时间步上的位移响应
}{},{2t t U U ∆∆，而是用其它方法如中心差分法，步长取t ∆的几分之一来求得。

Houbolt 方法是一个隐式差分格式，其步长可以取得比中心差分法大一些，不受cr t ∆的限制，但为了保证计算精度，步长也不能取得太大。

Houbolt 方法的求解步骤总结如下：
1．用有限元素法形成结构的刚度矩阵][K 、质量矩阵][M 、阻尼矩阵][C
2．计算初始值}{}{}{0
00U U U &&& 3．取步长t ∆，求
9,2,22,3
,5,611,23
7063504322120a
a a a a a a a t a t a t a t a ==-=-=====
∆∆∆∆ （5－16）
4． 求}{}{2t t U U ∆∆（一般用一个特殊的启动格式如，中心差分法） 5．求解
][][][]ˆ[1
0C a M a K K ++= （5－17） 6．矩阵分解
T L D L K
]][][[]ˆ[= （5－18） 7． 对每一步长，求 （1）
})
{}{}{]([}{}{}{]([}{}ˆ{27532642t t 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 P P ∆∆∆∆∆∆----++++++++= （5－19）
（2） }ˆ{}{]][][[t t t t T P U L D L ∆∆++= （5－20） （3） 若需要，求出
}{}{}{}{}{26420t t t t t t t t t U a U a U a U a U ∆∆∆∆--++---=&& （5－21）
}{}{}{}{}{27531t t t t t t t t t U a U a U a U a U ∆∆∆∆--++---=& （5－22）
【Wilson －θ法】
Wilson －θ法假定加速度从时刻t 到时刻t t ∆+θ为线性变化，所以，可以认为它是线性加速度法的推广。

在37.1≥θ时，它是无条件稳定的。

通常取4.1=θ。

具体方法为：
令τ为时间增量，其中t ∆θτ≤≤0，则在时刻t 到时刻t t ∆θ+的区间，有：
}){}({}{}{t t t t t U U t
U U &&&&&&&&-+=++∆θτ∆θτ
（5－23）
积分上式： }){}({2}{}{}{2
t t t t t t U U t
U U U &&&&&&&&-++=++∆θτ∆θττ （5－24） }){}({61}{21}{}{}{32t t t t t t t U U t
U U U U &&&&&&&-∆+++=∆++θττθττ （5－25）
当t ∆θτ=
}){}({2
}{}{t t t t t t U U t U U &&&&&&++=++∆θ∆θ∆θ （5－26）
})2}({6
}{}{}{2
2t t t t
t t t U U t U t U U &&&&&+++=++∆θ∆θ∆θ∆θ （5－27）
由此解出：
}2}{6}){}({6}{2
2t t t t t t t U U t U U t U &&&&&-∆--∆=∆+∆+θθθθ （5－28） }2
}{2}){}({3}{t t t t t t
t U t U U U t U &&&&∆---∆=∆+∆+θθθθ （5－29） 在t t ∆θ+时刻的平衡方程为：
}{}]{[}]{[}]{[t t t t t t t t P U K U C U M ∆+∆+∆+∆+=++θθθθ&&& （5－30）
其中
}){(}{}{t t t t t t P P P P -+=∆+∆+θθθ （5－31）
为t t ∆θ+时刻的载荷。

Wilson －θ法的具体步骤为：
（1） 用有限元素法形成结构的刚度矩阵][K 、质量矩阵][M 、阻尼矩阵][C
t ∆
{t U &&}t t ∆θ+&&
（2） 计算初始值}{}{}{0
00U U U &&& （3） 取步长t ∆，取4.1=θ，计算：
6
,2,31,,,2,2,3,)(62
8762
50
43
12120t
a t a a a a a a t a a a t a t a ∆∆θθθ
∆θ∆θ∆θ==-=-
==
==== （5－32）
（4）形成
][][][]ˆ[1
0C a M a K K ++= （5－33） （5）矩阵分解
T L D L K
]][][[]ˆ[= （5－34） （6） 对每一步长t ∆，计算时刻t t ∆θ+时刻的有效载荷：
})
{}{2}{]([}){2}{}{]([}){}({}{}ˆ{3
1
20t
t
t
t t t t t t t t t U a U U a C U U a U a M P P P P &&&&&&++++++-+=++∆∆θθ（5－35） （7）求解t t ∆θ+的位移：
}ˆ{}{]][][[t
t t t T P U L D L ∆θ∆θ++= （5－36） （7） 计算在t t ∆+时刻的加速度、速度、位移：
将t t U ∆+θ代入（5-28）解出t
t U ∆+θ&&，再代入（5-23）式，并令t ∆=τ，得到： }{}{}){}({}{654t t t t t t t U a U a U U a U &&&&&++-=∆+∆+θ （5－37）
将（5-23）代入（5-24）、（5-25）式中，令t ∆=τ，得到：
}){}({}{}{7t t t t t t U U a U U &&&&&&++=∆+∆+ （5－38） }){2}({}{}{}{8t t t t t t t U U a U t U U &&&&&+++=++∆∆∆ （5－39）
【Newmark 方法】
Newmark 方法也可以认为是线性加速度方法的推广。

其假定如下：
t U U U U t
t t t t t ∆δδ∆∆}]{}){1[(}{}{+++-+=&&&&& （5－40） 2}]{}){2
1[(}{}{}{t U U t U U U t
t t t t t t ∆αα∆∆∆+++-++=&&&&& （5－41） αδ,为积分常数。

根据积分精度和计算的稳定性来确定。

通常取2)5.0(25.0,5.0δαδ+=≥，可以保证算法是无条件稳定的。

当6
1
,21==αδ时，就是线性加速度法。

其具体方法为，根据t t ∆+时刻的平衡方程：
}{}]{[}]{[}]{[t
t t t t t t t P U K U C U M ∆∆∆∆++++=++&&& （5－42） 由最初的两个假定式（5-40）（5-41），求出用}{t t U ∆+表示的}{t t U ∆+&&和}{t
t U ∆+&的表达式，代入上方程求出}{t t U ∆+，然后回代求出}{t t U ∆+&&和}{t
t U ∆+&。

具体步骤为：
（1） 用有限元素法形成结构的刚度矩阵][K 、质量矩阵][M 、阻尼矩阵][C
（2） 计算初始值}{}{}{0
00U U U &&& （3） 取步长t ∆以及参数δα,（5.0)5.0(25.02≥+=δδα），计算：
t
a t a t a a a t a t a t a ∆δδ∆αδ
∆αδα∆α∆αδ∆α=-=-=-=-====
76543212
0),1(),2(2,1,121
,1,,1 （5－43） （4）形成
][][][]ˆ[1
0C a M a K K ++= （5－44） （5）矩阵分解
T L D L K
]][][[]ˆ[= （5－45） （6） 对每一步长t ∆，计算时刻t t ∆+时刻的有效载荷：
})
{}{}{]([}){}{}{]([}{}ˆ{5
4
1
320t
t
t
t t t t t t t U a U a U a C U a U a U a M P P &&&&&&++++++=∆+∆+ （5－46）
（7）求解t t ∆+的位移：
}ˆ{}{]][][[t
t t t T P U L D L ∆+∆+= （5－47） （8）计算在t t ∆+时刻的速度、加速度：
}{}{}){}({}{320t t t t t t t U a U a U U a U &&&&&---=++∆∆ （5－48） }{}{}{}{76t
t t t t t U a U a U U ∆∆++++=&&&&&& （5－49）
对一个具体问题，究竟采用何种积分格式，要根据各自的经验来选择，为了保证算法的稳定性，一般使用Wilson-θ法和Newmark 法，从使用上看，当用
于求解振动响应时，Newmark 法应用较多一些。

另外，求解微分方程的积分格式，常用的还有四阶龙格－库塔方法，它主要用于求解自由振动响应。

对于非线性系统的振动响应分析，由于严格的讲，不能用模态坐标对方程进行解耦，因而直接积分方法成为一种求解非线性系统振动响应的有效方法。

§5.3 方程解耦和模态响应
【方程解耦】
对于线性振动系统，如果求出了系统的固有振型][Φ，则可以通过模态坐标变换
}]{[}{q U Φ= （5－50）
将方程
}{}]{[}]{[}]{[P U K U C U M =++&&& （5－51）
变为解耦的形式：
}{}]{[}]{[}]{[r r r r P q K diag q C diag q M diag =++&&& （5－52）
这里假定][C 阵可以对角化（或者是比例阻尼，或者是强制对角化）。

从而使一个n 阶耦合微分方程组变成了n 个独立的微分方程，求解物理坐标}{U 的响应问题，变成了求解模态坐标}{q 的响应问题。

实际求解时，只需对其中一个方程进行。

【时域分析】
第r 个模态坐标的振动方程为：
r r r r r r r P q K q C q M =++&&& （5－53）
如果该单自由度系统对于单位脉冲输入的响应是)(t h ，则)}({t P r 在输入的τd 时间内引起的响应为τττd t h P r )()}({-，从而
⎰-=t
r r d t h P t q 0)()()(τττ （5－54）
已知有阻尼系统的单位脉冲响应函数为：
)1sin(11)(2
2
t e M t h r r t r
r r r r ωζωζωζ--=
- （5－55）
第r 个模态坐标的稳态响应为（Duhamel 积分）：
ττζωτωζτωζd t e P M t q r r t t
r r
r r s r r r )](1sin[)(11)(2
)(0
2
---=
--⎰
（5－56）
当考虑初始条件引起的瞬态响应时，
)1sin()(2r r r t r t r t e q t q r r θωζωζ+-=- （5－57）
常数r q ，r θ由出初始条件决定。

【频域分析】
当输入为简谐力或可以分解为简谐力的叠加时（如周期力），则可以采用频域分析方法求解振动响应。

第r 个模态坐标的方程：
t P q K q C q M r r r r r r r Ω=++cos &&& （5－58）
其解为：
)cos(r r r t q q ϕ-Ω= （5－59）
代入方程（5－58）得到：
r
r r r r r r r
r M K C C M K P q 2222tan )()(ΩΩ
ϕΩΩ-=
+-=
（5－60）
即：
)cos()
()()(2
2
2r r r r r
r t C M K P t q ϕ-ΩΩ+Ω-=
（5－61）
§5.4 结构的瞬态响应
瞬态响应是结构振动中常见的现象之一，它反映了结构对冲击激励的动态响应。

实际的结构大多会受到冲击的作用，且对于一般激励，可以将其看成是一系列脉冲激励的叠加。

研究结构的瞬态振动响应也有助于研究结构在一般激励下的响应。

当结构受到一个突加的（有限作用时间）非周期力)(t F 作用时，通常不会产生稳态振动过程，而是产生一个称之为瞬态响应的过程，这种非周期激励
力可以分解为一系列脉冲激励力的叠加，故我们先看脉冲激励和脉冲响应。

【脉冲激励和脉冲响应】
由于多自由度系统的振动方程可以通过模态坐标变换，变换为由若干个模态坐标下的单自由系统的振动方程的叠加，故我们考察单自由度系统在脉冲激励力下的响应。

1. 单位脉冲
冲量定义为力对时间的积分，记为：
⎰
=dt t f F )(ˆ （5－62） 一个幅值为ε/)(ˆt F 、作用时间为ε的冲击力，在0→ε时的力趋于无穷大，
但力的冲量是一个有限值。

当)(ˆt F
为单位冲量时，)(t F 在0→ε的情况下称为单位脉冲，在数学上，单位脉冲表示为δ函数。

δ函数有如下性质：
⎰
⎰∞
∞
=-=-≠=-0
0)
()()(1)(0)(ττδτδτ
τδf dt t f t dt t t t
－63）
2．单位脉冲响应
结构在冲量作用下，动量会发生变化，但位移无变化，故受脉冲作用的结
构（以质量为m 的单自由度系统为例）会产生一个初速度m F
/ˆ，已知单自由度系统在初始条件下的瞬态响应为：
t x t x t x n n n
ωωωcos )0(sin )
0()(+=
& （5－64）
t m F
t x n n
ωωsin ˆ)(= （5－65）
对单位脉冲，则
t m t x n n
ωωsin 1
)(=
（5－66） 冲量F
ˆ作用下，当系统有粘性阻尼时， )1sin()(2ϕωζζω--=-t Ae t x n t n （5－67）
代入初条件
m
F
x x ˆ)0(,0)0(=
=& （5－68） )(ˆ)1sin(1ˆ)(22
t g F
t e m F t x n t n n ⋅=--=
-ωζζωζω （5－69） 其中，)(t g 称为单位脉冲响应。

【任意激励下的响应】
任意激励力可以分解为一系列脉冲激励力的线性组合，因此任意激励力下的响应可以用一系列脉冲激励响应的线性叠加来得到。

考虑一个作用在时刻τ=t 的脉冲，它产生的冲量是
τ∆τ)(ˆf F
= （5－70） 在τ=t 时刻的单位脉冲响应为)(τ-t g ，故对冲量τ∆τ)(ˆf F
=，它在时刻τ=t 产生的响应为)()(ττ∆τ-⋅t g f
根据叠加原理，在时间t 的响应为：
⎰-=t
d t g f t x 0)()()(τττ （5－71）
【响应谱】
研究结构在冲击激励下的脉冲响应，一个主要的指标是响应的最大值，它是结构对脉冲激励的反应严重程度的一个度量。

显然，它与系统动力学特性有关，通常用响应谱来表示，由于冲击响应的反应时间很短，即使系统有阻尼，阻尼机制也来不及产生作用，因此，常常用无阻尼单自由度系统对脉冲激励的响应来研究冲击响应谱，简称响应谱。

它是单自由度系统最大峰值响应与其固有频率间的关系，不同类型的脉冲形成不同类型的响应谱。

在实际分析中，是将响应的最大值画成随T t /1变化的函数曲线来表示。

1
t 为脉冲特征时间如脉冲的间隔、持续时间等，T 为系统的自然周期。

纵坐标也用无量纲量)//(0max k F x 来表示，0F 为脉冲力的幅值。

对矩形脉冲的响应为：
1X m ax
1
⎪⎪⎪
⎩
⎪
⎪⎪⎨⎧>-=---≤≤-=1110101
)2(2sin 2sin 2)
(cos (cos 0)cos 1()(t
t t t t k F t t t k F t t t k F t x n n n n n ωωωωω（5－72） 则其冲击谱如图所示。

§5.5 模态叠加法
除了用直接积分法来求解结构的振动响应外，对于线性系统的振动，工程中也常用模态叠加法来进行振动响应的求解。

它是建立在模态空间中的一种方法。

其理论基础是固有模态正交性和展开定理。

【模态位移法】
物理坐标下结构的振动方程
)}({}]{[}]{[}]{[t P u K u C u M =++&&& （5－73）
初始条件为
}{)}0({}{)}0({00u u u u &&== （5－74）
根据其相应的自由振动方程解得固有模态阵为：
}]{}{}{}[{][321n φφφφΦΛ= （5－75）
利用模态坐标变换（展开定理）
}]{[}{ηΦ=u （5－76）
及固有模态正交性，对方程进行解耦，得到模态坐标系下结构的振动方程：
),2,1()
(n i t P K C M i i i i i i i Λ&&&==++ηηη （5－77） 其中，
)}
({}{)(}
]{[}{}]{[}{}]{[}{t P t P C C K K M M i T
i i i T
i i i T i i i T i i φφφφφφφ==== （5－78）
对初始条件，同样可以利用展开定理，变换成模态坐标下的初始条件:
}]{[}{}]{[}{0000ηΦηΦ&&==u u （5－79）
}
]{[}]{][[][}]{[][}]{[}]{][[][}]{[][000000ηηΦΦΦηηΦΦΦ&&&i T
T
i T T M diag M u M M diag M u M ==== （5－80）
}
]{[}{1}]{[}{1
0000u M M u M M T i i
i T i i
i &&φηφη==
),2,1(n i Λ= （5－81）
利用（5－81）的2n 个初始条件，求解（5－77）的n 个单自由度振动方程（其解有2n 个待定系数），回代到坐标变换式（5－76）中，就可以求得物理坐标系下结构的振动响应。

【模态截断】
从上述过程可以看到，如果结构的自由度数较大，则计算量仍然很大，在实际工程振动问题中发现，并非结构的所有模态都会被外激励激发而对结构响应有所贡献，或者即使有贡献，其成分也相当小，忽略这些次要模态的影响并不会使振动响应的计算精度产生明显降低。

所以，在实际使用模态叠加法进行振动响应计算时，并不需要求出所有的模态响应进行叠加，而只需求出前面若干阶对振动响应有明显贡献的那些模态的响应（主要是前若干阶低阶模态）。

因此在进行模态坐标变换时，就只需要在整个模态序列中截取前面若干阶低阶模态进行坐标变换，这就是所谓的模态截断。

注意在进行模态截断时，一般是从最低阶模态开始，顺序截取前面m 阶模态作为保留模态集，而不是随意选取几阶模态。

至于选取多少个模态，则要根据经验来确定，比如，根据经验，可以取8+=l m 或l m 2=等，l 是认为对响应有明显贡献的模态数。

采用模态截断后，则坐标变换式成为：
∑===m
i i i t u
1}{}ˆ]{ˆ[)}(ˆ{ηφηΦ （5－82） 从而解耦后的振动方程为：
),2,1()
(m i t P K C M i i i i i i i Λ&&&==++ηηη （5－83）
【模态位移法】
我们以一个简单的情况来分析。

假定系统无阻尼，且受到简谐激励力作用，则
}cos {)}({t P t P Ω= （5－84）
模态激励力为：
t P t P t P t P i T i T i i Ω=Ω=Ω=cos cos }{}{}cos {}{)(φφ （5－85）
模态坐标下的振动方程为：
t P K M i i i i i Ω=+cos ηη&& （5－86）
响应为：
t P H t i i Ω⋅Ω=cos )()(η （5－87）
)(ΩH 为频响函数，当输入为t P i Ωcos 时，输出的稳态响应为t Ωcos η，即
响应与激励同频。

代入到方程（5－85）得：
2
2ΩηηΩηi i i
i i
i i i i M K P P K M -=
=+- （5－88）
根据频响函数的定义，
2
1
cos cos )(Ω-=
ΩΩ=
Ωi i i i M K t
P t
H η （5－89）
模态响应为：
t
D t K P t P M K t i i
i i i i i i i ΩΩ-=
ΩΩ-=ΩΩ-=
cos )(1cos ])
(11
[cos 1)(2
2
2ωωη （5－90）
其中，
i
i
i K P D =
（5－91） 称为模态静变形。

考虑模态截断后，系统的近似稳态响应为：
t K P t u
m
i i i i i ΩΩ-=∑=cos ])(11
)[}({)}(ˆ{12
ωφ （5－92） 值得注意的是，在进行模态截断时，还要考虑到外激励的频率成份，如果某阶模态频率刚好与外激励的某个频率相等或接近，则即使该阶模态阶数超过了预先确定的模态截断的最高阶数，该阶模态也不能被截除而必须保留，以保证结果有足够的精度。

因此，在进行模态截断之前，应该进行外激励力的频谱分析和结构的固有频率分析。

在外激励的频率成分较多时，采用模态位移法计算的结果精度不高，要得到较高精度的响应解，就必须采用更多的模态，但随模态数的增加，收敛速度也较慢。

为了解决这个问题，可以采用模态加速度法。

模态加速度法改善了收敛特性，可以用较少的模态截断数目获得精度较好的结果。

【模态加速度法】
为了简化公式表达，采用无阻尼、无刚体位移的多自由度系统，来说明模态加速度法的原理： 由运动方程：
)(}]{[}]{[t P u K u M =+&& （5－93） })]{[}({][}{1u M P K u &&-=- （5－94）
用模态位移法中的}ˆ{u &&近似代替}{u &&，则：
∑=----=-=m
i i i M K P K u M P K u 1
111
}]{[][}{][})ˆ]{[}({][}ˆ{ηφ&&&&（5－95）
由特征值问题
}0{}]){[]([2=-i i M K φω （5－96）
得到
}{1
}]{[][2
1i i
i M K φω
φ=
- （5-97a ）
∑=--=m
i i i i
P K u
1
2
1
}{1
}{][}ˆ{ηφω
&& （5－97b ）
称第一项为“伪静响应项”，第二项称为“模态加速度项”，正是由于模态加速度项中的因子
2
1
i ω，改善了收敛速度，因为随着i 的增加，对应项对总响
应的影响就越小，当01,
2
→∞→i
i ω。

（5－97b ）式中的模态加速度i η&&求解如下：
ττωτωωηωωηηd t P M t t t i
t
i
i
i i i i
i i i )(sin )(1
sin )0(1
cos )0()(0
-+
+
=⎰& （5－98）
根据变上限积分的导数：
⎰⎰+=t t
t t f dt t f dt
d d t f dt d 00),()],([),(τττ （5－99） 有：
τ
τωτωωττωωτττωτd t P t t t P d t P d t P dt d i t
i i i i i i t
i i t
i )(cos )()(sin )()(cos )()(sin )(0
00
-=-+-=-⎰⎰⎰ （5－100）
)
()(sin )()(cos )()(sin )()(cos )(0
00t P d t P t t t P d t P d t P dt
d i i i i
t
i i i i i i t
i i t
i ωττωτωωττωωτττωτ+--=-+--=-⎰⎰⎰ （5－101）
从而求得：
⎰--+--=t
i i i
i i i i i i i i i i d t P M M t P t
t t 0
2
)(sin )()(sin )0(cos )0()(τ
τωτωωηωωηωη&&& （5－102） 由（5－97b ）和（5－102）就组成模态加速度法求响应的公式。

对于简谐激励，模态位移和模态加速度响应为：
t K P t K P i i i i
i i i i ΩΩ-⋅Ω-=ΩΩ-=
cos ])
/(11[cos ])
/(11[22
2ωηωη&& （5－103）
显然，当i ωΩ→时，∞→i η&&，故即使采用模态加速度法，在进行模态截
断时，也必须计及模态频率与激励频率接近或相等的那些模态，才能保证精度。

对有阻尼、无刚体位移的多自由度系统，在阻尼矩阵可以对角化的假设下，其模态方程为：
i
i i i i i i i M t P )
(22=
++ηωηωζη&&& ),2,1(n i Λ= （5－104） 用模态加速度法求得：
i i m
i i i m
i i i
i
t P K t u
ηφωηφωζ&&&}{)1
(
}){2()}({][)}(ˆ{1
2
1
1
∑∑==---=（5－105）
模态位移响应为：
τ
τωζτζωωζηωζηζωωζηητωζωζωζd t e
P M t e
t
e t i i t t
i
i i i i i t
i i i i i i i i t i i i i i i i
i )(1sin )(11
1sin )]0()0([11
1cos )0()(2)
(0
2
22
2---+
-+-+-=----⎰& （5－106）
由上式求得i η&和i η&&，
代入（5－105）即可由模态加速度法求出阻尼系统的响应。

【精细积分法】
将结构强迫振动方程化成状态方程的形式，利用现代控制理论中对状态方程的求解方法，引入状态转移矩阵的概念，进行逐步的递推求解，就形成对动力学响应求解的精细积分方法。

精细积分方法只适合于处理形如：
)}({}]{[}{t r u H u +=&
（5－107）
的微分方程，称为系统的状态方程。

因此在应用精细积分方法之前，要先将结构振动的微分方程化成（5－107）的形式。

1． 自由振动响应的求解
对于自由振动，相应的运动状态方程形式为：
}]{[}{u H u =&
（5－108）
}{u 称为系统的状态向量。

根据矩阵微分方程解的理论，其通解为：
}){]ex p([}{0u t H u = （5－109）
令时间步长为τ，记
])exp([][τH T = （5－110）
称为状态转移矩阵，则在一系列等步长的时刻:
ΛΛτττ
k t t t t k ====20210 （5－111）
系统在τ时刻的状态与初始状态之间有关系：
}]{[}{)}({01u T u u ==τ （5－112）
且有状态递推公式：
}]{[}{1k k u T u =+ （5－113）
因此，精细积分的计算归结于状态转移矩阵][T 的计算。

只要能精细计算出状态转移矩阵，就可以得到各个时刻系统状态的高精度解。

将积分步长τ划分为N m 2=等分，即：
ττ∆N m t -==2/ （5－114）
一般取20=N ，于是ττ∆620102--≈=t ，从而：
m m t H m H H T )]][ex p([)]/][ex p([)]ex p([][∆ττ=== （5－115）
展开成泰勒幂级数形式并取前四项得到：
]
[][2/)12/)]([)3/]([]([)]([][][!4/)]([!3/)]([!2/)]([][][)]ex p([22432a T I t H t H I t H t H I t H t H t H t H I t H +=++⨯++=++++≈∆∆∆∆∆∆∆∆∆ （5－116）
2/)12/)]([)3/]([]([)]([][][22t H t H I t H t H T a ∆∆∆∆++⨯+= （5－117）
从而：
Λ
=+⨯+⨯+⨯+=+⨯+=+=------)
2()
2()
2()
2()
1()
1(2
2
2
22
22])[]([])[]([])[]([])[]([])[]([])[]([])[]([][N N N N N N N
a a a a a a a T I T I T I T I T I T I T I T （5－118）
由于
][][][2][])[]([])[]([a a a a a T T T I T I T I ⨯++=+⨯+ （5－119）
因此令
][][][2][a a a a T T T T ⨯+= （5－120）
迭代计算上式N 次，就可得到：
][][][a T I T += （5－121）
泰勒展开中，略去的高阶项为)(5t O ∆，即)(10530τO -，显然，这样得到的][T 已经具有很高的精度。

2． 强迫振动响应
强迫振动响应的运动方程也可以化成状态方程的形式：
)}({}]{[}{t r u H u +=& （5－122）
假定在时间步),(1+k k t t 内，)}({t r 是线性的，即
)}({}{}]{[}{10k t t r r u H u -++=&， k t t =时，}{}{k u u = （5－123）
令)(k t t -Φ是对应齐次方程的解，即：
]][[][ΦΦH =& ][)]0([I =Φ （5－124）
从而，方程（5－122）的解为：
)})]
({}{][}({][}{][}
{][})][{([}{111
01
11
01k k k t t r r H r H r H r H u t t u -++-++-Φ=---- （5－125）
由于
][)]([)]([1T t t k k ==-+τΦΦ
所以：
)})]
({}{][}({][}{][}
{][][{[}{111
01
11
011τr r H r H r H r H u T u k k ++-++=----+ （5－126）。