2.3 多自由度结构的地震反应
2.3.1 多自由度结构地震运动方程的建立
前面研究的都是一个自由度的体系,质点的位置由质点的水平位移唯一地决定。多层建筑、地基、土坝可看作为由许多质点组成的多质点系,或多自由度体系。下面以两质点体系为例,推导得出地震作用下的运动方程式,然后将它推广到多质点体系的一般情况。
为了计算简便,我们先研究自由振动两自由度体系,如图2 10(a)中的两层钢筋混凝土框架结构,当框架横梁的线刚度视作无限大时,可简化成图210所示的框架计算简图。
图210 两层框架的计算简图
假定质点在振动过程中只发生水平向的移动而无转动,此图中m2、h2、k2及δ2为顶层框架的质量、层高度、柱的剪切刚度及水平位移;m1、h1、k1及δ1为底层框架的质量、层高度、柱的剪切刚度及水平位移。现分别取m1、m2为隔离体,见图211,根据达朗贝尔原理,可写出微分方程
··
-m1δ
㊣
㊣㊣
1-k1δ1+k2(δ2-δ1)=0-m2δ
··
(2 38)
2-k2(δ2-δ1)
=0
即
··
m1δ
㊣
㊣㊣
1+(k1+k2)δ1-k2δ2=0 m2δ
··
(2 39)
2-k2δ1+k2δ2
=0
令k11=k1+k2,k12=-k2,k21=-k2,k22=k2。
式(2 39)变为
··
m1δ
㊣
㊣㊣
1+k11δ1+k12δ2=0 m2δ
··
(2 40)
2+k21δ1+k22δ2
=0
式中包含系数k11、k22、k12、k21,都有明确的物理意义,即:
k11———点2保持不动,为使点1产生单位位移而在点1处所需施加的水平力;
k22———点1保持不动,为使点2产生单位位移而在点2处所需施加的水平力;
图2 11 质点m1、m2
k12———点1保持不动,由于点2产生单位位移而在点1处所引起的弹性反力;
为脱离体
k21———点2保持不动,由于点1产生单位位移而在点2处所引起的弹性反力。
观察式(240)很容易得出公式中包含了弹性恢复力及质点受到的惯性力。它们可以分别表示为
㊣㊣
s1=k11δ1+k12δ2
㊣
(2 41)
s2=k21δ1+k22δ
2
若将其写成矩阵形式
{S}=[k]{δ}
(2 42)
在两自由度体系中[k]=[kk1211 kk1222],kij为刚度系数,它反映结构刚度的大小,表示j
质点产生单位位移在i质点上产生的弹性恢复力。
两自由度体系在自由振动过程中只受弹性恢复力和惯性力,因此振动方程为(以矩阵形式表示)
··
[m]{δ}
+[k]{δ}=0
(2 43)
若体系在外力作用下产生振动,并考虑阻尼的影响,参照式(214)我们可以写出含有阻尼系数和外力的动力方程。
··
·
[m]{δ}
+[c]{δ}
+[k]{δ}={F}
(2 44)
在两自由度体系中,cij为阻尼系数,它表示j质点发生单位速度时,在i质点上引起阻尼力。在两自由度体系中,矩阵形式为
·
[R]=[c]{δ}
(2 45)
式中[c]———阻尼矩阵,[c]=㊣㊣cc1211 cc1222㊣㊣。
在两自由度体系中,惯性力为
··
··
F1=m1 δ
(
g+δ
1)
㊣
㊣㊣
··
··
(2 46)
F2=m2 δ
{
g+δ
2
}
写成矩阵的形式
[F]=㊣㊣m01 m02㊣㊣{δ}+㊣㊣m01 m02㊣㊣{δ}
(2 47)
同理,根据达朗贝尔原理可以写出两自由度体系强迫振动方程
{FF12++RR12++SS12==00
将式(242)、式(245)及式(247)代入上式可得
m1δ+c11δ+c12δ+k11δ1+k12δ2=-m1δ
㊣m2δ㊣㊣+c21δ+c22δ
+k21δ1+k22δ2=-m2δg
若写成矩阵形式
··
·
[m]{δ}
+[c]{δ}
+[k]{δ}=-[m]{I}δ
·· g
(2 48)
式中[m]、[c]、[k]同上;
{δ}、{δ}、{δ}———简化质点的相对位移、速度、加速度向量。{δ}={δδ…δ}T{δ}={δδ…δ}T{δ}={δ1δ2…δn}T
式(248)不仅适用于两自由度体系的强迫振动方程,而且适用于多自由体系的强迫振动方程。
2.3.2 多自由度体系的运动自振特性
研究多自由度体系自由振动的主要目的之一是求解该结构体系的自振频率,以便可以应用地震反应谱来衡量该结构的地震反应。由前面学到的知识我们可容易写出多自由度体系的自由振动方程和强迫振动方程。当阻尼值不大时有阻尼体系和相应的无阻尼体系的自振频率数值相当接近,因此我们在求解频率过程忽略阻尼项,为了简便以两自由度体系的自由振动为例。
现考虑m1及m2作同频率及同相位的简谐振动,则具有下列的解答
{δδ12((tt))==XX12ssiinn((ωωtt++φφ))
(2 49)
将式(249)代入式(243)中得
(k11-m1ω2)X1+k12X2=0㊣k21X1+(k22-m2ω2)X2
㊣㊣
(2 50)
=0
式(250)为振幅方程,由于振幅X1、X2不可能同时为0,所以奇次线性方程组有非0解,所以
k11-m1ω2
k12
k21
k22-m2ω2 =0
(2 51)
所以(-m1ω2+k11)(-m2ω2+k22)-k12k21=0就是振动频率方程。
此方程是关于ω2的一元二次方程,则ω2的两个正实根为
1
ω21.2=
2m1m2 [(m1k22+m2k11 )±㊣(m1k22+m2k11 )2-4m1m2(k11k22-k12k21)]
(2 52)
其中数值较小者为ω1,称第一自振频率;数值较大者为ω2,称第二自振频率。相对应的
自振周期为:T1=ω2π1,T2=ω2π2。
推广到一般的多自由度体系振幅方程变为
㊣
(k11-m1ω2)X1+k12X2=0
k21X1+(k22-m2ω2)X2+k23X3=0
…
㊣
(2 53)
ks,s-1Xs-1+(ks,s-msω2)Xs+ks,s+1Xs+1=0
…
㊣
kn,n-1Xn-1+(kn,n-mnω2)Xn
=0
写成矩阵的形式
[k]-[m]ω2 {X}=0
这是关于{x}的齐次方程,{x}=0说明体系没产生振动,无意义,所以由齐次方程有非零解可得。频率方程为
[k]-ω2[m]=0
(2 54)
式中[k]———刚体矩阵;
[m]———质量矩阵,[m]=㊣㊣m…01…㊣…m0…n㊣㊣m×n。
对于将横梁刚度视作无限大的剪切变形结构,有
㊣
k11 k12
0
㊣
k21 k22
k23
k32
k33
k34
[k]=
…
ks,s-1
ks,s
ks,s+1
…
㊣
0
kn,n-1
kn,
n㊣n×n
式(254)展开后是关于频率的n次方程,求解后可以求得n个正实根。2.3.3 多自由度体系主振型及其正交性
2.3.3.1 主振型
对于两自由度体系而言,由频率方程可求出ω1、ω2,然后将其分别代入振幅方程,不能求出相应的位移幅值X1、X2,因两个式并非是独立的,只能由其中任一式子求出振幅比。
(k11-m1ω2)X1+k12X2=0㊣k21X1+(k22-m2ω2)X2
㊣㊣
=0
即
XX21=m1ωk21-2k11 XX21=m2ωk221-k22
(2 55)
当ω=ω1时,振幅比
XX21=m1ωk211-2k11=m2ωk2121-k22
当ω=ω2时,振幅比
XX21=m1ωk221-2k11=m2ωk2221-k22
(2 56)
因为体系的ω、m、k、δ为定值,所以振幅比是一个常数,与时间无关。说明在振动过程中,各质点的比值始终保持不变,对应某一个自振频率就有一个振幅比,体系就按某一弹性曲线形状发生振动,这种振动形式通常称为主振型简称振型。对应于第一自振频率ω1称为第一振型或基本振型;对应于第二自振频率ω2称为第二振型。
例如第一振型、第二振型用Φ表示,以两自由度体系的阵型为例。
Φ=[Φ1 Φ2]=㊣㊣xx1211 xx1222㊣㊣=(12-21)
我们可以按比例绘出该体系的振动图形
(图2 12)。
图212 振型
因为主振型只取决于各集中质点的振幅之间的比值大小,而与这些振幅的绝对值大小无关,所以有时为计算或绘制主振型曲线方便,可令某一集中点振幅为1.0,其余集中
质点的振幅,按相应比值关系作相应的增大或缩小,以保持原来的弹性曲线的形状。
在一般初始条件下,每一个质点上的振动都是由各振型的简谐振动叠加而成的复合振动,而本身不再是简谐振动。质点之间位移的比值也不再是常数,而是随时间变化的。反过来说体系的振动曲线常可分解为若干个主振型的叠加。因此多自由度弹性系的运动方程的解写成
n
δi(t)=Σ
Xjisin(ωjt+φj)
j=1
其中每项只是其特解。将其展开数学表达式为
{δδ12((tt))==XX1112ssiinn((ωω11tt++φφ11))++XX2212ssiinn((ωω22tt++φφ22))
(2 57)
一般情况下第一振型曲线偏在表示原始位置的竖直线的一个侧面而不与其相交,第二振型曲线有一个穿越点;第三振型的曲线有两个穿越点,以此类推。主振型这个概念很有用,在研究多自由度体系的强迫振动时将用到它。
2.3.3.2 主振型的正交性
从两自由度体系的两个振型状态下的质点受力和变位分析,主振型的弹性曲线可看作是体系按某一频率振动时体系上相应惯性力所引起的静力弹性曲线形状。如图213所示,在ω1振型下弹性曲线可以视作由集中质量m1、m2上分别作用惯性力m1ω21X11、m2ω21X12产生的静力位移曲线,在ω2振型下弹性曲线可以视作由集中质量m1、m2上分别
作用惯性力m1ω22X21、m2ω22X22产生的静力位移曲线。
图213 两自由度体系的主振型图
根据虚功原理:“第一个力在其加力点由于第二个力引起的位移上所做的功等于第二个力在其加力点由第一个力引起的位移所做的功”,则图2 13(a)中惯性力对此虚位移做的静力虚功为
T12=m1ω21X11X21+m2ω21X12X22
同理
T21=m1ω22X21X11+m2ω22X22X12
T12=T21
所以
m1ω21X11X21+m2ω21X12X22=m1ω22X21X11+m2ω22X22X12
(ω21-ω22)(m1X11X21+m2X22X12 )=0
因为
ω1≠ω2
所以
m1X11X21+m2X22X12=0
(2 58)
式(258)即自由振动两个主振型具有相互正交的特性。对于两个以上的多自由度体系,任意两个振型j与k之间(j≠k)都有上述正交的特性,可以表示为
m1Xj1Xk1+m2Xj2Xk2+…+mnXjnXkn=0
n
或
Σ
miXjiXki=0
(2 59)
i=1
用矩阵表示为
{x}jT[m]{x}k=0(j≠k)
(2 60)
式中
{x}jT={Xj1,Xj2,Xj3,…,Xjn}
㊣
{x}k=
㊣
Xk1 Xk2 Xk3
㊣
㊣
㊣
…X
kn㊣
㊣
m1
㊣
m2
[m]=
㊣
㊣
m
n㊣
式(260)表示多自由度体系任意两个振型对质量矩阵的正交性。因为
{[k]-ω2k[m]}{X}k=0[k]{X}k=ω2k[m][X]k
两边同乘{X}jT,得
{X}jT[k]{X}k=ω2k[X]jT[m][X]k
由于
[X]jT[m][X]k=0
所以
[X]jT[k][X]k=0(j≠k)
(2 61)
上式表示多自由度体系任意两个振型对刚度矩阵的正交性。式(260)和式(2 61)都是振型的重要特性,对下节振型分解法公式的推导和演化起关键的作用。
2.3.4 多自由度体系地震反应的振型分解法
由前面学到的知识可知多自由度弹性体系的振动方程如下
··
·
[m]{δ}
+[c]{δ}
+[k]{δ}=-[m]{I}δ
·· g
此式中质点的位移δi(t)是以几何坐标来描述的。用这种坐标系表示的多质点体系的振动方程是耦连的,自由度数较多时会给求解方程带来困难,运用手工运算直接求解几乎是不可能的。若采用坐标变换的方法,利用振型的正交性将方程组解耦,使其变为一组各自独立的方程,而每个方程中只含有一个未知量从而可单独求解,最终会使计算量大大简化。这里提到的坐标变换是指采用振型为基底的广义坐标来描述质点的位移。
若采用的振型为基底的广义坐标来描述质点的位移,就可以利用振型的正交性使多自由度体系的振动方程组内的方程式互相解耦而各自独立,从而使多自由度体系的地震反应计算大大地简化。
2.3.4.1 阻尼矩阵的假定
振动方程组的耦连关系表现在阻尼矩阵[c]与刚度矩阵[k]中,而质量矩阵和刚度矩阵具有正交性,而阻尼矩阵不具有正交性,为消除阻尼矩阵各质点的耦连作用,通常可将阻尼矩阵[c]表示为[m]和刚度矩阵[k]的线性组合。
即
[c]=a[m]+b[k](a、b为比例常数)
(2 62)
2.3.4.2 广义坐标
以图214所示的三个自由度的悬臂柱为例,在任一时刻t各质点的位移分别为
δ1(t)、δ2(t)、δ3(t)可以用三个振型的线性组合来表示
图2143个自由度悬臂柱的振型曲线
㊣
㊣
㊣
δ1(t)=X11q1(t)+X21q2(t)+X31q3(t)δ2(t)=X12q1(t)+X22q2(t)+X32q3(t)δ3(t)=X13q1(t)+X23q2(t)+X33q3()
(2 63)
t
式中 δi(t)———质点在几何坐标的位移,它是时间t的函数;
Xji———第j振型的质点i处的幅值,以它作为广义坐标的基底;
qj(t)———第j个振型的广义坐标(j=1,2,3,…,n)。
对任意一时刻t来说,若式(2 63)中δ1(t)、δ2(t)、δ3(t)已知,就可以解出
qj(t),反之若qj(t)已知则可求出δ1(t)、δ2(t)、δ3(t)。因此几何坐标描述的位移与广义
坐标描述的位移是可以互为线性变换的。在多自由体系弹性动力反应计算中,通过式(2 63)的线性变换求解几何坐标的未知位移δ1 (t)的问题转化为求解广义坐标qj(t)的问
题。式(263)写成矩阵
3
{δ}={X1}q1+{X2}q2+{X3}q3=Σ
{Xj}qj
(2 64)
j=1
上述为三个自由度的几何坐标与广义坐标的关系,推广到n个自由度体系时得
n
{δ}={X1}q1+{X2}q2+…+{Xn}qn=Σ
{Xj}qj
(2 65)
j=1
或
{δ}=[X]{q}
(2 65a)
将上式对时间t求导
·
·
{δ}
=[X]{q}
(2 66)
··
··
{δ}
=[X]{q}
(2 67)
2.3.4.3 振型参与系数
为了利用阵型的正交性,将式(2 65a)左右两边分别乘{Xj}T[m],即得
{Xj}T[m]{δ}={Xj}T[m][X]{q}
(2 68)
因振型对质量矩阵的正交性得:上式除了{Xj}T[m]{Xj}≠0外,其余各项均为0。于是可得
{Xj}T[m][X]{q}={Xj}T[m]{Xj}qj
(2 69)
所以根据式(268)和式(269),得到
{Xj}T[m]{δ}={Xj}T[m]{Xj}qj
qj={Xj}T[m]{δ}
(2 70)
{Xj}T[m]{Xj}
若δ1=δ2=δ3=…=δn=1且用rj代替qj,则
rj={Xj}T[m][I]
{Xj}T[m]{Xj}
n
或
rj=Σ
miXji
i=1
n
(2 71)
Σ
miX2ji
i=1
rj称地震反应中第j振型参与系数,它实际上是各质点位移均等于1的第j振型的广义坐标值。由振型参数公式可容易得出Σ
n
riXji=1。
j=1
2.3.4.4 广义坐标微分方程及其解
㊣
[c]=a[m]+b[k]
㊣
{δ}={X1}q1+{X2}q2+…+{Xn}qn=[X]{q}{δ}
·
·
=[X]{q}
··
··
㊣
{δ}
=[X]{q}代入[m]{δ}+[c]{δ}+[k]{δ}=-[m]{1}δ,并左乘以{Xj}T得
{Xj}T{[m][X]{}+a[m][X]{}+b[k][X]{}+[k][X]{q}}=-{Xj}T [m]{1}δ
(272)
根据振型的正交性可知
{Xj}T [m][X]{}={Xj}T [m][Xj]
{Xj}T [m][X]{}={Xj}T [m][Xj]
{Xj}T [k][X]{}={Xj}T [k][Xj]=ω2j{Xj}T [m][Xj]
{Xj}T [k][X]{q}={Xj}T [k][Xj]qj=ω2j{Xj}T [m][Xj]qj
将上述关系式代入式(272)中得
{Xj}T [m]{Xj}+(a+bωj2){Xj}T[m]{Xj}+ω2j{Xj}T[m]{Xj}qj=-{Xj}T[m]{1}δ
将等式两侧除以{Xj}T[m]{Xj},并引入振型参数rj代入后则可写成
+2ξjω+ω2jqj=-rjδ
(2 73)
2ξjωj=a+bω2j,ξj=12(ωaj+bωj)
式中 ξj———第j个振型的阻尼比。
根据实测表明,结构各振型的阻尼比数量级相同,数值上高振型略大一些,但为简单起见在抗震分析中常采用一个阻尼比,即假定ξj=ξ。比较几何坐标系下的振动方程和广义坐标系下的振动方程
+2ξjω+ω2jqj=-rjδ
mδ(t)+cδ(t)+kδ(t)=-mδ(t)
可得
qj=-ωrjj∫t0δ(τ)e-ξjωj(t-τ)sinωj(t-τ)dτ
(2 74)
qj=rjΔj(t)(j=1,2,3,…,n)
Δj=ω-j1∫t0δ(τ)e-ξjωj(t-τ)sinωj(t-τ)dτ
(2 75)
至此多自由度体系的运动方程组通过振型分解后,可按单自由度进行计算,其步骤如下。
(1)求各振型的自振频率ωj,振型幅值xji,阻尼比ξj,通常(ξj=ξ),有阻尼的自振
频率ω′j取(ω′j=ωj)。
(2)求出各振型参与系数rj。
(3)按式(273)求出各振型等效单自由度位移地震反应Δj(t)。(4)最后按式(276)求原坐标下的质点i位移地震反应。
n
n
δi(t)=Σ
Xjiqi(t)=Σ
rjXjiΔ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),求出各时刻的输入地震加速度。
图215 Δt间隔内的加速度
线性变化
(2)假定在Δt时间间隔内,地震加速度与反应加速度均呈线性变化,如图215所示。因此,在Δt内(时刻t-Δt到时刻t之间),时刻τ的加速度{δ}
··
τ变化规律为{δ}τ={δ}t-Δt+({δ}t-{δ}t-Δt)τΔt(0≤τ≤Δt)
(2 77)
将上式对τ积分两次,分别可得
·
··
··
··
τ2
{δ}
τ={δ}
t-Δtτ+{δ}
(t-Δt)
t-{δ}
2Δt+C
(2 78){δ}τ={δ}t-Δtτ22+({δ}t-{δ}t-Δt)6τ3Δ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+Δ2t{δ}t-Δt+Δ2t{δ{δ}t={δ}t
(2 80){δ}t={δ}t-Δt+{δ}t-ΔtΔt+(2{δ}t-Δt+{δ}t)Δ6t2
(2 81)
由此得到{δ}t=Δ6t2{δ}t-{A}t
(2 82){δ}t=Δ3t{δ}t-{B}t
(2 83){A}t=2{δ}t-Δt+Δ6t{δ}t-Δt+Δ6t2{δ}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)~式(285)代入上式,得
{δ}t[K]={F}t
(2 86)
式中
[K]=[K]+6Δ[tM2]+3Δ[Ct]
{F}t={F}t+[M]{A}t+[C]{B}t
式(286)中的{δ}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)β法
纽马克β法是将线性加速度法加以修改,引入了γ、β两个参数。取结构体系的速度和位移表达式为
·
·
··
··
{δ}
t+Δt={δ}
t+(1-γ){δ}
[t+Δt]Δt
t+γ{δ}
㊣{δ}t+Δt={δ}t+{δ㊣㊣}tΔt+[(12-β){δ}t+β{δ}t+Δt](Δt)
(2 87)
2
··
将式(287)代入多自由度体系的运动方程即可求出{δ}
t+Δt,然后可求出相应的
·
t+Δt、{δ}t+Δt,式中的γ、β参数数值的选取是否适合对解法的稳定性和精度有直接的影响。根据经验和一些研究结果表明取γ=1/2不宜产生人为阻尼。若取β=1/4,式(2 87)就应用中点加速度法,即时间步长中间点处的加速度作为此时间步长内不变的加速度进行计算,其值为{δ}
{δ}
··
··
(t+Δt)
t+{δ}
。
2