数值分析底部剪力法与振型分解反应谱法对比分析

摘要：建筑结构抗震设计是建筑结构设计中必不可少，也是非常重要的一部分。结构抗震在建筑结构的总成本中占有相当大的比例。建筑抗震设计规范中有关于结构抗震计算的方法以及适用范围，水平地震力的计算方法主要是底部剪力法和振型分解反应谱法，底部剪力法适用于质量和刚度沿高度分布比较均匀的结构，而振型分解反应谱法能反应结构的真实情况，对一般结构都适用。本文通过对五层、八层、十层，质量和刚度分布均匀和不均匀框架结构的各层剪力计算，来比较两种方法的计算结果，验证底部剪力法的适用范围以及有效性。本文对结构特征周期的计算是用广义Jacobi方法，通过Fortran语言编程实现的。

关键词：底部剪力法；振型分解反应谱法； Jacobi方法；Fortran语言

Comparative Analysis between Equivalent Base Shear Method and Modal Analysis Method

Abstract: Seismic design plays an essential and important part in the structure design. It also makes up a significant proportion of the total cost. About the horizontal seismic force, the code has detailed specification of the calculation principle and applicable scope. The calculation method for horizontal seismic force mainly has the equivalent base shear method and modal analysis method. The equivalent base shear method is suitable for mass and stiffness along the height of structure with uniform distribution, and the modal analysis method reflects the true action of the structure and has a wide usage. By calculating the shear of five-story, eight-story and ten-story framework with mass uniform or non-uniform distribution, this paper verified the scope and the effectiveness of the equivalent base shear method. The eigenperiod of the structure is calculated by generalized Jacobi method though the Fortran language programming.

Key words: Equivalent Base Shear Method; Modal Analysis Method; Jacobi Method; Fortran Language

引言

实际的建筑结构其质量一般均是连续分布的，因此，严格的说，其动力自由度均是无限的。但采用无限自由度模型一方面计算过于复杂；另一方面也没有这种必要，因为选用有限自由度模型的计算结果已能充分满足一般情况下工程设计的精度要求[1]。因此，对于多高层房屋，一般均采用多自由度模型。对于建筑结构而言，一般每层楼面及屋面可作为一个质点，而楼面与楼面（屋面）之间墙、柱的质量则分别向上、向下集结到楼面及屋面质点处。这种多自由度模型即为工程上的层间模型。图1所示为这种层间模型的计算简图，这种简图即所谓的“具有n个质量的悬臂柱”，本文所讨论的即“悬臂柱具有n个质量，比较基底剪力法和振型分解反应谱法的结果”。

关于建筑结构水平地震力的弹性计算方法，建筑抗震设计规范中[2]的5.1.2条规定，各类建筑结构的抗震计算，应采用的方法中：

（1） 高度不超过40m、以剪切变形为主且质量和刚度沿高度分布比较均匀的结构，以及近似于单质点体系的结构，可采用底部剪力法等简化方法。

（2） 除（1）款外的建筑结构，宜采用振型分解反应谱法。

那么，规范为什么要如此规定呢，具体有哪些方面的考虑。本文将通过对框架结构的两种方法计算的结果进行对比分析，说明规范规定的原因及考虑，这里只考虑规则结构，也就是不计结构扭转耦联。

1 底部剪力法

对于一般的建筑结构，应采用振型分解反应谱法计算其地震作用效应，但当房屋结构满足下述条件时，可采用更为简便的底部剪力法计算其地震作用效应[1]：

（1）结构的质量和刚度沿高度分布比较均匀；

（2）房屋的总高度不超过40m；

（3）建筑结构在地震作用下的变形以剪切应变为主；

（4）建筑结构在地震作用时的扭转效应可忽略不计。 

满足上述条件的结构在地震作用下其反应通常以第一振型为主，且近似为直线。

图2 结构水平地震作用计算简图

建筑抗震设计规范文献5.2.1条规定：采用底部剪力法时，各楼层可仅取一个自由度，结构的水平地震作用标准值，应取下列公式确定（图2）：

   （1）（i=1,2,3…n） （2）

 （3）

其中，是结构总水平地震作用标准值；是相应于结构基本自振周期的水平地震影响系数，根据规范按图3进行计算；是结构等效重力荷载，单质点取总重力荷载代表值，多质点可取总重力荷载代表值的85%；是质点i的水平地震作用标准值；是集中于i质点的重力荷载代表值；是i质点的计算高度；是顶部附加地震作用系数按表1；是顶部附加水平地震作用。

图3 地震影响系数曲线[2]

2 振型分解反应谱法

适用范围：除上述可以采用底部剪力法以外的建筑结构，宜采用振型分解反应谱法。

建筑抗震设计规范5.2.2条规定：采用振型分解反应谱法时，不进行扭转耦联计算的结构，应按以下规定计算其地震作用和作用效应。

表1 顶部附加作用影响系数[2]

注：为结构基本自振周期。

（1）结构第j振型i质点的水平地震作用标准值，应按下列公式计算：

 （i=1,2,…n,j=1,2,…m）（4）

 （5）

其中，是j振型i质点的水平地震作用标准值；是相应于第j自振周期的地震影响系数；是j振型i质点的水平相对位移；是j振型的参与系数。

（2）水平地震作用效应（弯矩、剪力、轴向力和变形），当相邻振型的周期比小于0.85时，可按下式确定：

其中，是水平地震作用标准值的效应；是j振型水平地震作用标准值的效应，可只取前2~3个振型，当基本自振周期大于1.5s或是房屋高宽比大于5时，振型个数应当适当增加。

3 用Jacobi方法求解基本振型

3.1 经典Jacobi方法的计算基本步骤

经典Jacobi方法就是用一系列平面转换矩阵，从实对称矩阵A出发，逐次作相似变换而使之对角化，则此对角阵的对角元素就是矩阵A的全部特征值，而且所作正交变换的正交矩阵乘积的各个列向量即为矩阵A的相应的特征向量[3]。

经典Jacobi方法的基本步骤如下： 

（1）选定Ak-1的非对角元素的绝对值最大的元素；

（2）计算sinθ和cosθ的值

                      （6）

             （7）

           （8）

                     （9）

       （10）

于是平面旋转变换矩阵为：

（3）计算变换之后的Ak

由于对A进行平面旋转变换时，只是改变A的p行p列和q行q列，所以只需计算Ak的p行p列和q行q列元素即可。

                                 （11）                                 （12）                                 （13）                                 （14）                                 （15）（4）计算特征向量Xk

 （16）

                                （17）（18）

3.2 广义Jacobi方法的计算基本步骤

经典Jacobi方法只适用于求解A为实对称矩阵的标准特征值问题的特征值和特征向量。而广义Jacobi方法可以用来求解广义特征值问题的全部特征值和特征向量。在广义特征值问题中，A为实对称矩阵，B为正定近似对角矩阵[4]。

本例中，需要求解质量分布不均匀悬臂柱的特征值和特征向量，所以必须采用广义Jacobi方法来做。广义Jacobi方法的基本原理和经典Jacobi方法类似：

（1）按次序（1，2），（1，3），…，（1，n），（2，3），…，（n-1，n）的顺序选取（i，j）的值。用广义Jacobi矩阵代替经典Jacobi方法中的平面转换矩阵作合同的矩阵序列。

，，

（2）令， 

，

，

。

则当时，我们取，

。

否则，我们取，。

（3）计算变换之后的，

矩阵中各项计算公式：

  （19）

 （20）

 （21）

 （22）

 （23）

对于矩阵中的各项元素，只要把上式中的换成即可。

（4）进行迭代收敛判断 

其中，。

迭代收敛后，所求特征值为，特征向量为矩阵中对应的列向量[5]。否则，从第一步重新开始计算。

3.3广义Jacobi方法的Fortran程序[6]

程序内容

！CALCULATE EIGENVALUES BY JACOBI METHOD

SUBROUTINE JACO(K,M,U,E,N)

REAL K,M

INTEGER P,Q

DIMENSION M(N,N),K(N,N),U(N,N),E(N)

IF(N.EQ.1)THEN

E(1)=K(1,1)/M(1,1)

U(1,1)=1

RETURN

ENDIF

N1=N-1

FMK=0.0

FMM=0.0

DO 10 I=1,N

DO 10 J=1,N

FMK=AMAX1(FMK,ABS(K(I,J)))

FMM=AMAX1(FMM,ABS(M(I,J)))

10   CONTINUE

FAC=0.00001

EPSK=FMK*FAC

EPSM=FMM*FAC

KFAC=5

KMAX=KFAC*N*N

AK=1.0E+30

DO 30 I=1,N

DO 20 J=1,N

U(I,J)=0.0

20    CONTINUE

U(I,I)=1.0

30    CONTINUE

DO 100 KAISUU=1,KMAX

P=1

Q=2

IF(MOD(KAISUU,2).EQ.1)GOTO 50

AK=ABS(K(1,2))

DO 40 I=1,N1

I1=I+1

DO 40 J=I1,N

A=ABS(K(I,J))

IF(AK.GE.A)GOTO 40

P=I

Q=J

AK=A

40    CONTINUE

GOTO 70

50 CONTINUE

AM=ABS(M(1,2))

DO 60 I=1,N1

I1=I+1

DO 60 J=I1,N

A=ABS(M(I,J))

IF(AM.GT.A)GOTO 60

P=I

Q=J

AM=A

60    CONTINUE

70    CONTINUE

D1=K(P,P)*M(P,Q)-K(P,Q)*M(P,P)

D2=K(Q,Q)*M(P,Q)-K(P,Q)*M(Q,Q)

D3=K(P,P)*M(Q,Q)-K(Q,Q)*M(P,P)

D=D3*D3+4.0*D1*D2

S=SQRT(D)

X=(D3+S)*0.5

IF(D3.LT.0.0) X=(D3-S)*0.5

A=D2/X

G=-D1/X

DO 80 J=1,N

W=K(P,J)

K(P,J)=K(P,J)+G*K(Q,J)

K(Q,J)=K(Q,J)+A*W

W=M(P,J)

M(P,J)=M(P,J)+G*M(Q,J)

M(Q,J)=M(Q,J)+A*W

80    CONTINUE

DO 90 I=1,N

W=K(I,P)

K(I,P)=K(I,P)+G*K(I,Q)

K(I,Q)=K(I,Q)+A*W

W=M(I,P)

M(I,P)=M(I,P)+G*M(I,Q)

M(I,Q)=M(I,Q)+A*W

W=U(I,P)

U(I,P)=U(I,P)+G*U(I,Q)

U(I,Q)=U(I,Q)+A*W

90    CONTINUE

IF((AK.LT.EPSK).AND.(AM.LT.EPSM))GOTO 120

100    CONTINUE

WRITE(*,'("matrix K or M are wrong")')

120    CONTINUE

DO 130 I=1,N

E(I)=K(I,I)/M(I,I)

130    CONTINUE

DO 180 I=1,N1

I1=I+1

DO 170 J=I1,N

IF(E(I).LE.E(J))GOTO 170

RE=E(I)

E(I)=E(J)

E(J)=RE

DO 160 L=1,N

RU=U(L,I)

U(L,I)=U(L,J)

U(L,J)=RU

160    CONTINUE

170    CONTINUE

180    CONTINUE

RETURN

END

DIMENSION A(8000)

CHARACTER *8 CH1,CH2

WRITE(*,'("Enter the name of input data file1,please")')

READ(*,7)CH1

7    FORMAT(A8)

OPEN(1,FILE=CH1,STATUS='OLD')

WRITE(*,'("Enter the name of output data file,please")')

READ(*,7)CH2

OPEN(2,FILE=CH2,STATUS='NEW')

READ(1,*)N

L=N*N

CALL SUB1(N,A(1),A(1+L),A(1+2*L),A(1+3*L))

END

SUBROUTINE SUB1(N,M,K,U,E)

REAL M,K

DIMENSION M(N,N),K(N,N),U(N,N),E(N)

READ(1,*)((M(I,J),J=1,N),I=1,N)

READ(1,*)((K(I,J),J=1,N),I=1,N)

CALL JACO(K,M,U,E,N)

WRITE(2,'("The eigen value is")')

WRITE(2,*)(E(I),I,1,N)

DO 20 I=1,N

20    E(I)=SQRT(E(I))

WRITE(2,'("The angle frequency w is")')

WRITE(2,*)(E(I),I,1,N)

WRITE(2,'("Number of eigen vector and the eigen vectors are")')

DO 40 I=1,N

40    WRITE(2,*)I,(U(J,I),J=1,N)

RETURN

END

（a）五层框架结构

（b）八层框架结构

（c）十层框架结构

图4 质量均匀分布框架结构抗震计算简图

（a）五层框架结构

（b）八层框架结构

（c）十层框架结构

图5 质量非均匀分布框架结构抗震计算简图

4 实例分析

某5、8、10层钢筋混凝土框架结构如图3所示，结构质量现在考虑为均匀分布和不均匀分布两种情况，集中于楼盖和屋盖的重力荷载代表值如图4、图5所示。柱的截面尺寸为，采用C35的混凝土，，梁的刚度EI=∞，建筑场地为II类，抗震设防烈度为7度，设计地震分组为第二组，设计基本地震加速度为0.1g，结构阻尼比为ζ=0.05。

表2  质量均匀分布结构特征值和前三阶振型的特征向量

（a）五层框架结构

特征值（

特征值（

特征值（

（a）五层框架结构

特征值（

特征值（

特征值（

4.1 底部剪力法和振型分解反应谱法各层剪力值结果对比

4.1.1 质量均匀分布结构结果比较

把用广义Jacobi方法计算出的结构基本振型代入规范公式中，并查表进行计算，得到如下表所示结构在地震作用下各层剪力值：

表4  质量均匀分布结构底部剪力法和振型分解反应谱法各层剪力值结果对比

（a）五层框架结构

表5  质量非均匀分布结构底部剪力法和振型分解反应谱法各层剪力值结果对比

（a）五层框架结构

5.1 数据分析

图6 质量均匀分布剪力相对误差图

图7 质量非均匀分布剪力相对误差图

将表4和表5中两种方法计算得到的层间剪力相对误差画图（图5和图6）进行分析，可得如下结论：

（1）高度越高，两种方法的误差越大；

（2）结构质量分布越不均匀，两种计算方法误差越大；

（3）四层（20m）以下，两种方法计算误差在10%以内，误差可以接受。而当结构高度超过8层（40m）时，两者误差显著增加；

（4）对于顶层剪力值，两种计算方法的误差最大。

（5）从图6中可明显发现，用基底剪力法进行计算，采用顶部附加地震作用系数对结构顶部地震作用进行修正，可减小两种方法的误差。

5.2 结论

从以上分析可以看出，对于高度小于40m（本例中八层以下）的均匀结构，底部剪力法和振型分解反应谱法计算结果具有良好的吻合性。而当结构高度越高或者质量分布越不均匀时，两种方法的误差将会越大。因此，《建筑抗震设计规范》（GB50011-2010）中规定，对于高度不超过40m、以剪切变形为主且质量和刚度沿高度分布比较均匀的结构，以及近似于单质点体系的结构，可采用底部剪力法等简化方法。并且，因为采用顶部附加地震作用系数对结构顶部地震作用进行修正，可减小两种方法的误差。所以，当结构基本周期较长，场地特征周期较小时，应该按规范规定对结构顶部地震作用进行修正。

之所以认为振型分解反应谱法的结果更精确，是因为底部剪力法是通过对振型分解反应谱法进行化简得到的。但是对于单个锯齿形的反应谱而言，振型分解反应谱法的分析结果与单个波的时程分析，误差可以达到10-30%之间，因此在个别（特殊性）意义上而言，反应谱分析结果是有误差的，因此，规范规定对于复杂的或者高层建筑需要采用时程分析进行补充计算和验证。

参　考　文　献

[1]  吕西林，周德源，李思明，陈以一，陆浩亮.建筑结构抗震设计理论与实例（第三版）[M]. 上海：同济大学出版社, 2011

[2]  中华人民共和国行业标准. GB50011─2011建筑抗震设计规范[S].北京：中国建筑工业出版社，2010

[3]  何汉林.数值分析（第二版）[M]. 北京：科学出版社, 2011

[4]  曹志浩.矩阵特征值问题（第二版）[M]. 上海：上海科学技术出版社, 1983

[5]  沈冯强.广义Jacobi方法的优化算法[J]. 力学学报, 2010, 42(2):319-324

[6]  林晓彤. Fortran90编程基础[M]. 山东：青岛海洋大学出版社, 2006