2.3 多自由度结构的地震反应

2.3.1 多自由度结构地震运动方程的建立

前面研究的都是一个自由度的体系，质点的位置由质点的水平位移唯一地决定。多层建筑、地基、土坝可看作为由许多质点组成的多质点系，或多自由度体系。下面以两质点体系为例，推导得出地震作用下的运动方程式，然后将它推广到多质点体系的一般情况。

为了计算简便，我们先研究自由振动两自由度体系，如图2 10（a）中的两层钢筋混凝土框架结构，当框架横梁的线刚度视作无限大时，可简化成下图所示的框架计算简图，如图2 10（b）、2 10（c）所示。

假定质点在振动过程中只发生水平向的移动而无转动，此图中m₂、h₂、k₂及δ₂为顶层框架的质量、层高度、柱的剪切刚度及水平位移；m₁、h₁、k₁ 及δ₁ 为底层框架的质32

图2 10 两层框架的计算简图

量、层高度、柱的剪切刚度及水平位移。现分别取m₁、m₂为隔离体，见图2 11，根据达朗贝尔原理，可写出微分方程：

{-m₂δ··

-m₁δδ（δ12+-kδ21（δ）₂=-0δ₁）=0

（2 38）

即

··

{m₂δm₁δ+-k（k2δ11++kk22δ）δ21=-0k₂δ₂=0

（2 39）

令k₁₁=k₁+k₂，k₁₂=-k₂，k₂₁=-k₂，k₂₂=k₂。式（2 39）变为

{m₂δ··

m₁δ++kk1211δδ11++kk1222δδ22==00

（2 40）

式中 k₁₁———点2保持不动，为使点1产生单位位移而在点1

图2 11 质点m₁、m₂

处所需施加的水平力；

为隔离体

k₂₂———点1保持不动，为使点2产生单位位移而在点2

处所需施加的水平力；

k₁₂———点1保持不动，由于点2产生单位位移而在点1处所引起的弹性反力；k₂₁———点2保持不动，由于点1产生单位位移而在点2处所引起的弹性反力。

观察式（2 40）很容易得出公式中包含了弹性恢复力及质点受到的惯性力。它们可以分别表示为

{ss12==kk1211δδ11++kk1222δδ22

（2 41）

若将其写成矩阵形式：

{S}=[k]{δ}

（2 42）

在两自由度体系中[k]=「㊣||kk1211 kk1222㊣」||；k_ij反映结构刚度的大小，称为刚度系数，表示j质点产生单位位移在i质点上产生的弹性恢复力。

两自由度体系在自由振动过程中只受弹性恢复力和惯性力，因此振动方程为（以矩阵形式表示）：

[m]{δ}+[k]{δ}=0

（2 43）

若体系在外力作用下产生振动，并考虑阻尼的影响，参照式（2 14）我们可以写出含有阻尼系数和外力的动力方程：[m]{δ}+[c]{δ}+[k]{δ}={F}

（2 44）

在两自由度体系中，c_ij为阻尼系数，它表示j质点发生单位速度时，在i质点上引起阻尼力。在两自由度体系中，矩阵形式为[R]=[c]{δ}

（2 45）

式中[c]———阻尼矩阵，[c]=㊣㊣|cc1211 cc1222㊣㊣|。

在两自由度体系中，惯性力为

{F₂=m₂{（δ··

··

F₁=m₁ δ

g+δ+δ）}

（2 46）

写成矩阵的形式：

[F]=㊣㊣|m01 m02㊣㊣|{δ}+㊣㊣|m01 m02㊣㊣|{δ}

（2 47）

同理，根据达朗贝尔原理可以写出两自由度体系强迫振动方程：

{FF12++RR12++SS12==00

将式（242）、式（245）及式（247）代入上式可得

{mm12δδ++cc1211δδ++cc1222δδ++kk1211δδ11++kk1222δδ22==--mm12δδ

若写成矩阵形式：

[m]{δ}+[c]{δ}+[k]{δ}=-[m]{I}δ

（2 48）

式中

[m]———质量矩阵；[c]———阻尼矩阵；[k]———刚度矩阵；

{δ}、{δ}、{δ}———简化质点的相对位移、速度、加速度向量：{δ}={δδ…δ}^T{δ

}={δδ…δ}^T

{δ}={δ₁δ₂…δ_n}^T

式（2 48）不仅适用于两自由度体系的强迫振动方程，而且适用于多自由体系的强迫振动方程。

2.3.2 多自由度体系的运动自振特性

研究多自由度体系自由振动的主要目的之一是求解该结构体系的自振频率，以便可以应用地震反应谱来衡量该结构的地震反应。由前面学到的知识我们可容易写出多自由度体系的自由振动方程和强迫振动方程。当阻尼值不大时有阻尼体系和相应的无阻尼体系的自振频率数值相当接近，因此我们在求解频率过程忽略阻尼项，为了简便以两自由度体系的自由振动为例。

现考虑m₁及m₂作同频率及同相位的简谐振动，则具有下列的解答：

δδ12（（tt））==XX12ssiinn（（ωωtt++φφ））}

（2 49）

将式（249）代入式（243）中得

k（k211X1-1+m₁（kω222）-Xm1+2ωk212）XX22==00}

（2 50）

式（2 50）为振幅方程，由于振幅X₁、X₂不可能同时为0，因此奇次线性方程组有非0解，所以

k₁₁-m₁ω²

k₁₂

k₂₁

k₂₂-m₂ω² =0

（2 51）

所以（-m₁ω²+k₁₁）（-m₂ω²+k₂₂）-k₁₂k₂₁=0就是振动频率方程。

此方程是关于ω²的一元二次方程，则ω²的两个正实根为

ω²_1.2=2m11m₂[（m₁k₂₂+m₂k₁₁）±

（m₁k₂₂+m₂k₁₁ ）²-4m₁m₂（k₁₁k₂₂-k₁₂k₂₁）]

（2 52）

其中数值较小者为ω₁，称第一自振频率；数值较大者为ω₂，称第二自振频率。相对应的自振周期为：T₁=ω2π1，T₂=ω2π2。

推广到一般的多自由度体系振幅方程变为

（k₁₁-mk₂11ωkXs2，s1）-X+1X1（k+s-2k121+-2Xm（k22s…ω=，s20-k）_nXm，ns2-ω…+1X2k）n2X-3X1s++3k=（ks0，sn+，n1X-_sm+1nω=²0）X_n=0㊣㊣㊣||||||||

（2 53）

写成矩阵的形式

[k]-[m]ω² {X}=0

这是关于{x}的齐次方程，{x}=0说明体系没产生振动，无意义，所以由齐次方程有非零解可得。频率方程为

[k]-ω² [m]=0

（2 54）

式中[k]———刚体矩阵；

[m]———质量矩阵，[m]=㊣㊣|||m︙01…⋱…m0︙n㊣㊣|||_m×n。对于将横梁刚度视作无限大的剪切变形结构，即

[k]=「㊣||||||||||kk01211kkk123222kkks…，s2333-1kkkn…，sn3，4s-1kks0n，s，+n1㊣」||||||||||_n×n

式（2 54）展开后是关于频率的n次方程，求解后可以求得n个正实根。2.3.3 多自由度体系主振型及其正交性

2.3.3.1 主振型

对于两自由度体系而言，由频率方程可求出ω₁、ω₂，然后将其分别代入振幅方程，不能求出相应的位移幅值X₁、X₂，因两个式并非是独立的，只能由其中任一公式求出振幅比。

{k（k211X1-1+m₁（kω222）-Xm1+2ωk212）XX22==00

即

XX21=m₁ωk21-2k₁₁ XX21=m₂ωk221-k₂₂

（2 55）

当ω=ω₁时，振幅比

XX21=m₁ωk211-2k₁₁=m₂ωk2121-k₂₂

当ω=ω₂时，振幅比

XX21=m₁ωk221-2k₁₁=m₂ωk2221-k₂₂

（2 56）

因为体系的ω、m、k、δ为定值，所以振幅比是一个常数，与时间无关。说明在振动过程中，各质点的比值始终保持不变，对应某一个自振频率就有一个振幅比，体系就按某一弹性曲线形状发生振动，这种振动形式通常称为主振型简称振型。对应于第一自振频率ω₁称为第一振型或基本振型；对应于第二自振频率ω₂称为第二振型。

例如若第一振型、第二振型用Φ表示，以两自由度体系的阵型为例：

Φ=[Φ₁ Φ₂]=㊣㊣|xx1211 xx1222㊣㊣|=㊣㊣|12-21㊣㊣|

我们可以按比例绘出该体系的振动图形（图2 12）。

因为主振型只取决于各集中质点的振幅之间的比值大小，而与这些振幅的绝对值大小无关，所以有时为计算或绘制主振型曲线方便，可令某一集中点振幅为1.0，其余集中质点的振幅，按相应比值关系作相应的增大或缩小，以保持原来的弹性曲线的形状。

在一般初始条件下，每一个上质点的振动都是由各振型的简谐振动叠加而成的复合

图2 12 振型

振动，而本身不再是简谐振动。质点之间位移的比值也不再是常数，而是随时间变化的。反过来说体系的振动曲线常可分解为若干个主振型的叠加。因此多自由度弹性系的运动方程的解写成

n

δ_i（t）=∑

X_jisin（ω_jt+φ_j）

j=1

其中每项只是其特解，将其展开数学表达式为

{δδ12（（tt））==XX1112ssiinn（（ωω11tt++φφ11））++XX2212ssiinn（（ωω22tt++φφ22））

（2 57）

一般情况下第一振型曲线偏在表示原始位置的竖直线的一个侧面而不与其相交，第二振型曲线有一个穿越点；第三振型的曲线有两个穿越点，以此类推。主振型这个概念很有用，在研究多自由度体系的强迫振动时将用到它。

2.3.3.2 主振型的正交性

从两自由度体系的两个振型状态下的质点受力和变位分析，主振型的弹性曲线可看作是体系按某一频率振动时体系上相应惯性力所引起的静力弹性曲线形状。如图2 13所示，在ω₁振型下弹性曲线可以视作由集中质量m₁、m₂ 上分别作用惯性力m₁ω²₁X₁₁、m₂ω²₁X₁₂产生的静力位移曲线，在ω₂振型下弹性曲线可以视作由集中质量m₁、m₂上分别作用惯性力m₁ω²₂X₂₁、m₂ω²₂X₂₂产生的静力位移曲线。

图2 13 两自由度体系的主振型图

根据虚功原理：“第一个力在其加力点由于第二个力引起的位移上所做的功等于第二个力在其加力点由第一个力引起的位移所做的功”，则图2 13（a）中惯性力对此虚位移做的静力虚功为

T₁₂=m₁ω²₁X₁₁X₂₁+m₂ω²₁X₁₂X₂₂

同理

T₂₁=m₁ω²₂X₂₁X₁₁+m₂ω²₂X₂₂X₁₂

T₁₂=T₂₁

所以

m₁ω²₁X₁₁X₂₁+m₂ω²₁X₁₂X₂₂=m₁ω²₂X₂₁X₁₁+m₂ω²₂X₂₂X₁₂

（ω²₁-ω²₂）（m₁X₁₁X₂₁+m₂X₂₂X₁₂ ）=0

因为

ω₁≠ω₂

所以

m₁X₁₁X₂₁+m₂X₂₂X₁₂=0

（2 58）

式（2 58）即自由振动两个主振型具有相互正交的特性。对于两个以上的多自由度体系，任意两个振型j与k之间（j≠k）都有上述正交的特性，可以表示为

m₁X_j1X_k1+m₂X_j2X_k2+…+m_nX_jnX_kn=0

n

或

∑

m_iX_jiX_ki=0

（2 59）

i=1

用矩阵表示为

{x}_j^T [m]{x}_k=0（j≠k）

（2 60）

式中{x}_j^T={X_j1，X_j2，X_j3，…，X_jn}；

{x}_k=㊣㊣㊣||||||XXXX…kkkk123n㊣㊣㊣||||||；

㊣

m₁

㊣

m₂

[m]=

|||||

⋱

|||||

。

㊣

m_n

㊣

式（2 60）表示多自由度体系任意两个振型对质量矩阵的正交性。因为

{[k]-ω²_k [m]}{X}_k=0[k]{X}_k=ω²_k [m][X]_k

两面同乘{X}_j^T，得

{X}^t_j[k]{X}_k=ω²_k[X]_j^T [m][X]_k

由于

[X]_j^T[m][X]_k=0

所以

[X]_j^T [k][X]_k=0（j≠k）

（2 61）

上式表示多自由度体系任意两个振型对刚度矩阵的正交性。式（2 60）和式（2 61）都是振型的重要特性，对下节振型分解法公式的推导和演化起关键的作用。

2.3.4 多自由度体系地震反应的振型分解法

由前面学到的知识可知多自由度弹性体系的振动方程如下：38[m]{δ}+[c]{δ}+[k]{δ}=-[m]{I}δ

此式中质点的位移δ_i（t）是以几何坐标来描述的。用这种坐标系表示的多质点体系的振动方程是耦连的，自由度数较多时会给求解方程带来困难，运用手工运算直接求解几乎是不可能的。若采用坐标变换的方法，利用振型的正交性将方程组解耦，使其变为一组各自独立的方程，而每个方程中只含有一个未知量从而可单独求解，最终会使计算量大大简化。这里提到的坐标变换是指采用振型为基底的广义坐标来描述质点的位移。

若采用的振型为基底的广义坐标来描述质点的位移，就可以利用振型的正交性使多自由度体系的振动方程组内的方程式互相解耦而各自独立，从而使多自由度体系的地震反应计算大大地简化。

2.3.4.1 阻尼矩阵的假定

振动方程组的耦连关系表现在阻尼矩阵[c]与刚度矩阵[k]中，而质量矩阵和刚度矩阵具有正交性，而阻尼矩阵不具有正交性，为消除阻尼矩阵各质点的耦连作用，通常可将阻尼矩阵[c]表示为[m]和刚度矩阵[k]的线性组合，即

[c]=a[m]+b[k]（a、b为比例常数）

（2 62）

2.3.4.2 广义坐标

以图2 14所示的三个自由度的悬臂柱为例，在任一时刻t各质点的位移分别为δ₁（t）、δ₂（t）、δ₃（t），可以用三个振型的线性组合来表示：

δδδ123（（（ttt）））===XXX111123qqq111（（（ttt）））+++XXX222123qqq222（（（ttt）））+++XXX333123qqq333（（（ttt）））㊣㊣㊣||||

（2 63）

式中 δ_i（t）———质点在几何坐标的位移，它是时间t的函数；

X_ji———第j振型的质点i处的幅值，以它作为广义坐标的基底；q_j（t）———第j个振型的广义坐标（j=1，2，3，…，n）。

图2143个自由度悬臂柱的振型曲线

对任意一时刻t来说，若式（2 63）中δ₁（t）、δ₂（t）、δ₃（t）已知，就可以解出q_j（t），反之若q_j（t）已知，则可求出δ₁（t）、δ₂（t）、δ₃（t）。因此几何坐标描述的位移与广义坐标描述的位移是可以互为线性变换的。在多自由体系弹性动力反应计算中，通过式（2 63）的线性变换求解几何坐标的未知位移δ₁ （t）的问题转化为求解广义坐标q_j（t）的问题。式（2 63）写成矩阵：

3

{δ}={X₁}q₁+{X₂}q₂+{X₃}q₃=∑

{X_j}q_j

（2 64）

j=1

上述为三个自由度的几何坐标与广义坐标的关系，推广到n个自由度体系时得

n

{δ}={X₁}q₁+{X₂}q₂+…+{X_n}q_n=∑

{X_j}q_j

（2 65）

j=1

或

{δ}=[X]{q}

（2 65a）

将上式对时间t求导{δ}=[X]{}

（2 66）{δ}=[X]{}

（2 67）

2.3.4.3 振型参与系数

为了利用阵型的正交性，将式（2 65a）左右两边分别左乘{X_j}^T[m]，即得

{X_j}^T [m]{δ}={X_j}^T [m][X]{q}

（2 68）

因振型对质量矩阵的正交性得：上式除了{X_j}^T [m]{X_j}≠0外，其余各项均为0。于是可得

{X_j}^T [m][X]{q}={X_j}^T [m]{X_j}q_j

（2 69）

所以根据式（268）和式（269），得到

{X_j}^T [m]{δ}={X_j}^T [m]{X_j}q_j

q_j={X_j}^T [m]{δ}

（2 70）

{X_j}^T [m]{X_j}

若δ₁=δ₂=δ₃=…=δ_n=1，且用r_j代替q_j，则

r_j={X_j}^T [m][I]

{X_j}^T [m]{X_j}

n

或

r_j=∑

m_iX_ji

i=1

n

（2 71）

∑

m_iX²_ji

i=1

r_j称地震反应中第j振型参与系数，它实际上是各质点位移均等于1的第j振型的广义坐标值。由振型参数公式可容易得出∑

n

r_iX_ji=1。

j=1

2.3.4.4 广义坐标微分方程及其解

㊣㊣㊣||||||{{{[cδδδ}]}}====a{[XX[[Xm1]}]]q{q+{+}}[k{]X₂}q₂+…+{X_n}q_n=[X]{q}

代入[m]{δ}+[c]{δ}+[k]{δ}=-[m]{1}δ，并左乘以{X_j}^T得{X_j}^T{[m][X]{}+a[m][X]{}+b[k][X]{}+[k][X]{q}}=-{X_j}^T [m]{1}δ

（272）

根据振型的正交性可知

{X_j}^T [m][X]{}={X_j}^T [m][X_j ]

{X_j}^T [m][X]{}={X_j}^T [m][X_j ]

{X_j}^T [k][X]{}={X_j}^T [k][X_j ]=ω²_j{X_j}^T [m][X_j ]

{X_j}^T [k][X]{q}={X_j}^T [k][X_j ]q_j=ω²_j{X_j}^T [m][X_j ]q_j

将上述关系式代入式（2 72）中得

{X_j}^T [m]{X_j}+（a+bω²_j ）{X_j}^T [m]{X_j}+ω²_j{X_j}^T [m]{X_j}q_j=-{X_j}^T [m]{1}δ

将等式两侧除以{X_j}^T [m]{X_j}，并引入振型参数r_j代入后则可写成

+2ξ_jω+ω²_jq_j=-r_jδ

（2 73）

2ξ_jω_j=a+bω²_jξ_j=12㊣㊣|ωaj+bω_j㊣㊣|

式中 ξ_j———第j个振型的阻尼比。

根据实测表明，结构各振型的阻尼比数量级相同，数值上高振型略大一些，但为简单起见在抗震分析中常采用一个阻尼比，即假定ξ_j=ξ。比较几何坐标系下的振动方程和广义坐标系下的振动方程

+2ξ_jω+ω²_jq_j=-r_jδ

mδ（t）+cδ（t）+kδ（t）=-mδ（t）

可得

q_j=-ωrjj∫t0δ（τ）e^-ξjωj（t-τ）sinω_j（t-τ）dτ

（2 74）

q_j=r_jΔ_j（t）（j=1，2，3，…，n）

Δ_j=ω-_j1∫t0δ（τ）e^-ξjωj（t-τ）sinω_j（t-τ）dτ

（2 75）

至此多自由度体系的运动方程组通过振型分解后，可按单自由度进行计算，其步骤如下：

（1）求各振型的自振频率ω_j，振型幅值x_ji，阻尼比ξ_j，通常ξ_j=ξ，有阻尼的自振频率ω′_j取ω′_j=ω_j。

（2）求出各振型参与系数r_j。

（3）按式（2 73）求出各振型等效单自由度位移地震反应Δ_j（t）。（4）最后按式（2 76）求原坐标下的质点i位移地震反应。

n

n

δ_i（t）=∑

X_jiq_i（t）=∑

r_jX_jiΔ_i（t）（j=1，2，3，…，n）

（2 76）

j=1

j=1

理论上分析时采用全部振型，可实际工程中通常只采用频率较小的（因为前几个振型的影响最大），控制了体系最大反应的前几个振型进行叠加就可以了。

2.3.5 多自由度体系地震反应的逐步积分法

振型叠加法在理论上是简明的，但在具体应用时仍较繁杂，且不适用于非线性动力分析。应用逐步积分法可以不求体系的自振频率，而直接解出各时刻的位移和应力。

2.3.5.1 线性加速度法

线性加速度方法与计算步骤如下。

（1）将输入地震加速度的计算时间分成许多微小时间间隔Δt（例如Δt=0.01s），求出各时刻的输入地震加速度。

图2 15 Δt间隔内的加速度

线性变化

（2）假定在Δt时间间隔内，地震加速度与反应加速度均呈线性变化，如图2 15所示。因此，在Δt内（时刻t-Δt到时刻t之间），时刻τ的加速度{δ}_τ变化规律为{δ}_τ={δ}_t-Δt+（{δ}_t-{δ}_t-Δt）τΔt（0≤τ≤Δt）

（2 77）

将上式对τ积分两次，分别可得

{δ}_τ={δ}_t-Δtτ+（{δ}_t-{δ}_t-Δt ）2τ²Δt+C

（2 78）{δ}_τ={δ}_t-Δtτ²2+（{δ}_t-{δ}_t-Δt ）6τ³Δt+Cτ+D

（2 79）

利用下列初始条件：

τ=0时，{δ}_τ={δ}_t-Δtτ=0时，{δ}_τ={δ}_t-Δt

可求出积分常数为C={δ}_t-Δt，D={δ}_t-Δt

将C、D代入式（278）和式（279），并令τ=Δt，便求得t时刻的速度和位移为{δ}_t={δ}_t-Δt+Δ2t{δ}_t-Δt+Δ2t{δ}_t

（2 80）{δ}_t={δ}_t-Δt+{δ}_t-ΔtΔt+（2{δ}_t-Δt+{δ}_t）Δt²

（2 81）

6

由此得到

{δ}_t=Δ6t²{δ}_t-{A}_t

（2 82）

{δ}_t=3

{δ}_t-{B}_t

（2 83）

Δt{A}_t=2{δ}_t-Δt+Δ6t{δ}_t-Δt+Δ6t²{δ}_t-Δt

（2 84）

{B}_t=Δ2t{δ}_t-Δt+2{δ}_t-Δt+Δ3t{δ}_t-Δt

（2 85）

由动力方程式，在时刻t有

[K]{δ}_t+[C]{δ}_t+[M]{δ}_t={F（t）}

将式（282）、式（283）、式（284）、式（285）代入上式，得

{δ}_t [K]={F}_t

（2 86）

式中

[K]=[K]+6Δ[tM2]+3Δ[Ct]

{F}_t={F}_t+[M]{A}_t+[C]{B}_t

式（2 86）中的{δ}_t可用一般的三角分解法求解。

（3）在{δ}_t求得后，代入式（2 83）求速度{δ}_t，代入式（2 82）求加速度{δ}_t。（4）根据{δ}_t求出应变{ε}_t和应力{σ}_t。

（5）按照（2）～（4）的步骤，计算下一时刻的位移、速度、加速度、应变、应力，一个一个时段地逐步进行，直至地震结束。

应当指出，采用该法求解动力方程时，Δt要分得足够小，以保证解答的稳定性。由于Δt很小，这样就会花费很长的计算时间。一些文献建议采用其他形式的差分方案以改进解答的稳定性。

计算实践表明，在各种不同的差分方案中，外点法（即威尔逊θ法）是一个较有效的方法。该法的要点是：用上述普通差分方法求时间增量为τ=θt的加速度列阵{δ}_{（t+Δt）+τ}（其中θ=1.37～1.40），然后再由它求出{δ}_t、{δ}_t、{δ}_t。容易看出，当θ=1时，外点法就等于上述普通方案。当θ=2时，外点法就是著名的双步长法。

2.3.5.2 纽马克（Newmark）β法

纽马克β法是将线性加速度法加以修改，引入了γ、β两个参数。取结构体系的速度和位移表达式为{{{δδ}}tt++ΔΔtt=={{δδ}}tt++{δ[·（1}t-Δtγ+）{[δ㊣㊣|12}t+-γβ{㊣㊣|δ{·δ}t}+tΔ+tβ]Δt{δ}_t+Δt]（Δt）²

（2 87）

将式（2 87）代入多自由度体系的运动方程即可求出{δ}_t+Δt，然后可求出相应的{δ}_t+Δt、{δ}_t+Δt，式中的γ、β参数数值的选取是否适合对解法的稳定性和精度有直接的影响。根据经验和一些研究结果表明取γ=1/2不宜产生人为阻尼。若取β=1/4，式（2 87）就应用中点加速度法，即时间步长中间点处的加速度作为此时间步长内不变的加速度进行计算，其值为（{δ}_t+{δ}_t+Δt ）/2。