LSSRLI-1与SHAKE2000都是一维土层反应分析程序,以一维剪切波在水平成层介质中的地震反应稳态解为基础编制的。

首先,底层输入频域地震位移为E_Nexp(ik_Nz)的地震波时, 第n层地震波位移的频域稳态一般解为

U_n(z)=E_nexp(ik_nz)+F_nexp(-ik_nz).(1)

式中, k_n=(1-id_n)·ω/c_n, d_n为第n层的阻尼比; c_n=(1+id_n)·(μ_n/ρ_n)^1/2,μ_n为第n层剪切模量。E_n、F_n为波幅系数,需要用传递矩阵进行求解:

[E_n+1

F_n+1]=[(1+α_n)/2exp(ik_nh_n)(1-α_n)/2exp(-ik_nh_n)

(1-α_n)/2exp(ik_nh_n)(1+α_n)/2exp(-ik_nh_n)][E_n

F_n].(2)

式中,α_n=ρ_nc_n/ρ_n+1c_<sub>n+1,为相邻两层的复波阻抗比值。(2)式右边左乘的矩阵可以根据每一层的土层信息计算出来。根据(2)式由E₁、F₁推出各层的E_n、F_n。但E₁、F₁开始是未知的,需要根据已知的波幅系数E_N来得到,关系如下

E₁=F₁=E_N/e_N.(3)

式中,E_N为已知的输入地震波波幅系数; e_N通过(2)式中的传递矩阵得到。由(2)、(3)式可递推出从1层到N层的波幅系数E、F。这样可以通过(1)式得到各层的位移,再进行傅式逆变换就得到了时域位移,进而可得速度、加速度。但是因为实际的土层是非线性系统,因此叠加原理失效,不能进行傅式变换,所以需要引入等效线性化的假设,让系统成为可以使用傅式变换的线性系统(廖振鹏,2003)。

根据建筑抗震设计规范(GB50011—2010)中场地分类的条款,构造了2个Ⅰ类、4个Ⅱ类场地的土层剖面。Ⅰ类场地模型因为比较薄,所以在程序中仅分为4个计算层,前3层为土层,第4层为基岩,土的物理力学特征见表1。Ⅱ类场地模型的计算层厚度一般为每层2～3 m,土层物理力学特征见表2。

表1 Ⅰ类场地土层的物理力学性质
Tab.1 Physical and mechanical properties of class-I sites

表2 Ⅱ类场地土的物理力学性质
Tab.2 Physical and mechanical properties of class-II sites

考虑到输入地震波对两程序计算结果(PGA、反应谱和剪应变)差异程度有影响,输入的地震波需要具有一定的代表性。本文选择5%阻尼比反应谱中长周期成分占优、短周期成分占优以及成分比较均匀的地震波,确定输入地震波为El-Centro地震波、AKTH19地震波和KSRH地震波(图1)。AKTH19与KSRH均为日本Kik-net台站代码,输入地震波信息见表3。从反应谱来看(图2),KSRH地震波短周期成分丰富,AKTH19地震波长周期成分丰富,El-Centro地震波介于之间,以此代表频率成分不均匀的地震波。将地震波都调幅至0.1 g、0.2 g与0.4 g进行输入,分别代表Ⅶ、Ⅷ和Ⅸ度区。这样,Ⅰ类场地的两个剖面、3条地震波和3个幅值,共有18个组合工况,同样Ⅱ类场地共有36个组合工况。

本文使用的非线性曲线是陈红娟(2009)提出的粘土的模量衰减曲线上限与阻尼比增长曲线的下限的组合,代表着土为弱非线性的情况,数值如表4所示。

表3 输入地震波基本信息
Tab.3 Basic information of the inputting waves

本文地震波 发震时间/年-月-日 台站地点 震级(M)方向 加速度峰值/gal 数据步长El-Centro 1940-05-18 El Centro Imperial Valley M7.1地震 N-S 326 0.02 sAKTH19 2011-03-11 AKTH19,Kik-net “3·11”M9.0大地震 N-S 20.89 0.01 sKSRH 2004-11-29 KSRH09,Kik-net M7.1地震 E-W 30.56 0.005 s

表4 土层非线性曲线
Tab.4 Non-linear curves of profiles

图1 输入El Centro(a)、AKTH19(b)和KSRH(c)不同地震波时的波形
Fig.1 Waveform of the inputting El Centro(a),AKTH19(b)and KSRH(c)earthquakes

图2 输入El Centro(a)、AKTH19(b)和KSRH(c)不同地震波时的加速度反应谱(5%阻尼比)
Fig.2 Acceleration response spectrum of the inputting El Centro(a),AKTH19(b)and KSRH(c)earthquakes(5% damping ratio)

结果表明,两个程序计算出的Ⅰ类场地的地表加速度峰值PGA的相对差变化范围是(0～1.4)%,平均值为0.6%,差异很小。从图3可见,加速度峰值的相对差随埋深变化,而且靠近地表的时候相对差较大; 剖面Ⅰ-2的相对差都比较小,沿深度变化不大。

两个程序计算出的Ⅱ类场地的地表加速度峰值PGA的相对差变化范围是(0.03～30.7)%,平均值为6.6%,可以看出,与Ⅰ类场地相比,差异程度的变化范围变大。由图4可见,不同深度相对差有所不同,而且最大的相对差并不一定出现在地表处。

不同输入地震波对PGA结果的影响见图5。由图可见,输入3种地震波的计算结果具有很相近的分布,相对差都集中在10%以内。这表明,在弱非线性硬场地条件下,对于不同地震波,LSSRLI-1与SHAKE2000的大部分PGA计算结果是相近的。

为了对比反应谱间的差异,得到两个程序计算的地表加速度反应谱的谱值比为

R_S=ln(S_L)-ln(S_S).(4)

其中,S_S是SHAKE2000的反应谱结果(g),S_L是LSSRLI-1的反应谱结果(g)。

图3 Ⅰ类场地加速度峰值相对差随埋深的变化(a)剖面Ⅰ-1工况;(b)剖面Ⅰ-2工况
Fig.3 Relative difference of acceleration peak value varies with depth in class-I sites (a)profile Ⅰ-1 condition;(b)profile Ⅰ-2 condition

图4 Ⅱ类场地加速度峰值相对差随埋深的变化(a)剖面Ⅱ-1工况;(b)剖面Ⅱ-2工况;(c)剖面Ⅱ-3工况;(d)剖面Ⅱ-4工况
Fig.4 Relative deviation of acceleration peak value varies with depth in class-II sites (a)profile Ⅱ-1 condition;(b)profile Ⅱ-2 condition;(c)profile Ⅱ-3 condition;(d)profile Ⅱ-4 condition

图5 不同输入地震波下Ⅰ、Ⅱ类场地加速度峰值相对差统计直方图
Fig.5 Relative deviation histogram of acceleration peak in class-Ⅰ and class-Ⅱsites under different inputting waves

R_S^-=∑ⁿ₁R_Si·ΔT_i〖JB<1*/〗∑ⁿ₁ΔT_i.(5)

其中,R_S^-是平均谱值比; ΔT_i是周期间隔(s); R_Si是每个周期对应的谱值比。

计算表明,Ⅰ类场地平均谱值比均值是0.3%,变化范围是(0.1～0.6)%; Ⅱ类场地的平均谱值比均值是2.3%,变化范围是(0.2～12.3)%。

图6a是Ⅰ类场地计算结果。从各工况结果可见,在周期T=0.03 s附近,R_S=0.6(相当于相对差为45.1%),两程序在这个周期附近的反应谱结果差异较大; 0.1 s之后,反应谱结果差异很小。

图6b是Ⅱ类场地计算结果。从各工况结果的曲线看,反应谱结果的相对差在3 s以前存在较大的波动,与Ⅰ类场地相比,反应谱的较大差异从短周期段发展到了长周期段。

谱比值反映的是反应谱间的相对差,不同输入地震波下Ⅰ、Ⅱ类场地反应谱相对差统计结果见图7。由图可见,输入3种地震波计算的结果具有相近的分布,而且不超过20%。这表明,在弱非线性硬场地条件下,两个程序在输入不同地震波时,反应谱计算结果是相近的。

LSSRLI-1与SHAKE2000给出了每一计算层内的剪应变峰值,在不同埋深处,剪应变峰值也不同,笔者对剪应变峰值随埋深变化曲线上的最大值进行比较。结果表明,Ⅰ类场地剪应变相对差变化范围是(0.7～9.0)%,均值为4.4%; Ⅱ类场地的剪应变相对差变化范围是(0.7～303.0)%,均值为46.4%。

图6 Ⅰ类(a)、Ⅱ类(b)场地反应谱谱值比汇总
Fig.6 Spectral ratio of response spectrum of class-Ⅰ(a)and class-Ⅱ sites(b)

图7 不同输入地震波下Ⅰ、Ⅱ类场地反应谱相对差统计直方图
Fig.7 Relative deviation histogram of response spectrum in class-Ⅰ and class-Ⅱsites under different inputting waves

按照不同的输入地震波,将剪应变峰值相对差进行统计,结果见图9。由图可见,3种地震波造成的剪应变峰值差异程度具有相近的分布,但剪应变峰值相对差大于100%的情况均来自于KSRH地震波输入。

经过对比分析发现,剪应变、反应谱和PGA相对差间有一定的关联,其中一个较大时,其余两个也较大。图 10a、b给出了一工况中剪应变差异与相应的反应谱差异,从图中可以看出,SHAKE2000计算得到的最大剪应变为17.7×10^-4,而LSSRLI-1为17.6×10^-4,相对差为0.5%。SHAKE2000计算的PGA为0.28 g,LSSRLI-1也为0.28 g,结果相同。反应谱最大相对差为16%。图 10c、d为另一工况的计算结果,由图可知,SHAKE2000计算得到的最大剪应变为83.1×10^-4,而LSSRLI-1的为144×10^-4,相对差为72.6%; SHAKE2000计算的PGA为0.81 g,而LSSRLI-1的为0.61 g,两者之间的相对差为24%。反应谱最大相对差为59%。由此可见,当剪应变差异越大,则PGA和反应谱的差异越大,三者之间存在一定的关联性。

两个程序的基本原理相同,输入的参数相同,计算结果的差异极有可能是程序内部某些参数的计算方法或算法不同造成的。在迭代过程中,会变化的参数是根据等效剪应变来取值的土层模量比、阻尼比。等效剪应变计算都是剪应变峰值乘以折减系数,计算中该系数均为0.65,即等效剪应变计算方法相同。但从前文第2.3节可以看出,两程序的剪应变峰值是不同的,所以两程序迭代计算得到的等效剪应变也是不同的。因此,迭代中剪应变的不同计算方法是导致不同计算结果的原因。此外,迭代中模量比、阻尼比的插值方法也可能不同。有研究表明,土层地震反应结果对模量比与阻尼比很敏感(刘红帅等,2005; 孙锐等,2009),这也能导致两程序的计算结果有差异。

为了表明两程序的计算结果差异来自内部计算方法不同,笔者给出一个算例。采用图 10c、d的工况,即在剖面Ⅱ-1输入0.4 g的El-Centro地震波。在这个工况中,SHAKE2000与LSSRLI-1

图8 Ⅰ类(a)、Ⅱ类(b)场地最大剪应变相对差频率分布直方图
Fig.8 Frequency distribution histogram of maximum shear strain deviation in class-Ⅰ site(a) and class-Ⅱ site(b)

图9 不同输入地震波下Ⅰ、Ⅱ类场地剪应变峰值相对差统计直方图
Fig.9 Relative deviation histogram of shear strain peak in class-Ⅰ and class-Ⅱsites under different inputting waves

图 10 剖面Ⅱ-1输入El-Centro地震波的剪应变及反应谱(a)输入0.1 g El-Centro地震波的剪应变;(b)输入0.1 g El-Centro地震波的反应谱;(c)输入0.4 g El-Centro地震波的剪应变;(d)输入0.4 g El-Centro地震波的反应谱
Fig.10 Shear strain and response spectrum of inputting El-Centro wave in profile Ⅱ-1(a)shear strain of inputting 0.1 g El-Centro earthquake wave;(b)response spectrum of inputting 0.1 g El-Centro earthquake wave;(c)shear strain of inputting 0.4 g El-Centro earthquake wave;(d)response spectrum of inputting 0.4 g El-Centro earthquake wave

的计算结果经过若干次迭代已经收敛,现称图 10c、d中的结果为“收敛结果”。在两程序已经收敛的结果上,用SHAKE2000再将其各自迭代一次,称这次结果为“再迭代结果”。如果SHAKE2000与LSSRLI-1内部计算方法相同,那么LSSRLI-1的再迭代结果应该与LSSRLI-1的收敛结果相同或相近。实施这次“额外”迭代的具体做法如下:将图 10c中SHAKE2000与LSSRLI-1的各层剪应变计算结果提取出来,结合前文给出的土非线性曲线(等效剪应变系数取0.65),可以插值获得对应两个程序中各土层的剪切模量与阻尼比,见表5。将这些模量与阻尼比当作初始参数都输入到SHAKE2000中进行1次迭代计算,LSSRLI-1的再迭代计算流程见图 11。因为只有1次迭代,所以剪应变计算结果只由各程序的剪应变计算方法决定,排除了多次迭代带来的其他差异累积,同样的还有反应谱的计算结果。剪应变、反应谱收敛、再迭代结果见图 12。

图 12中的LS代表输入表5中LSSRLI-1收敛结果的G、D列,再用SHAKE2000进行1次计算的再迭代结果,第一个L代表使用LSSRLI-1的

图 11 对LSSRLI-1收敛结果实施再迭代流程图
Fig.11 Flow chart for re-iterate on convergence results of LSSRLI-1

在图 12a中SS与S很接近,即d_S很小,这说明SS仅仅是在S本已收敛的结果基础上又迭代了一次; 而LS与L有大的偏差,与d_S相比,d_L很

大,这说明SHAKE2000程序与LSSRLI-1程序具有不同的剪应变计算方法。而从图 12b中可以看到,SS与S的反应谱很接近,d_S1很小,而LS与L之间的距离则略大,d_L1>d_S1。SS与S的反应谱也仅仅是多迭代一次而已,所以很接近; 而LS与L反应谱有较小的差异,可能是由输入参数上存在

表5 从各层的剪应变计算对应的剪切模量G、阻尼比D
Tab.5 Shear modulus G and damping ratio D calculated by shear strain in every layers

图 12 不同迭代次数下剪应变结果(a)和反应谱结果(b)(剖面Ⅱ-1输入El Centro地震波)
Fig.12 Shear strain(a)and response spectrum results(b)under different iteration number (inputting El Centro wave in profile Ⅱ-1)

的误差导致,也可能是反应谱计算方法存在不同,这些需要进一步研究。

从以上的对比可知,当剪应变不参与计算,即不依据剪应变来插值求得下一次计算所需各层土的剪切模量和阻尼比时,两程序计算结果并无太大差异,如图 12b中L与LS两条线已非常相近。但剪应变参与到计算当中时,两方法所得到的地表反应就相差较大,如图 12b中S与L两条反应谱。因而,可以初步推断,造成两程序间差异的主要原因是剪应变计算方法不同,当然,算法等方面也会有一些差异,但所带来的影响并不大。