多点地震动激励下的高效反应谱方法＊

王君杰郭进

1） 中国上海 200092 同济大学土木工程学院

2） 中国石家庄 050043 石家庄铁道大学土木工程学院

引言

地震动是空间变化的，由于波列传播速度的有限性和相干性损失以及局部场地地质条件不同，结构各支承点的地震激励出现显著差异，这种差异可能对平面大尺度结构（例如桥梁、水坝）的地震反应产生重要影响.欧洲规范 8 （European Committee for Standardization，2004）和城市轨道交通结构抗震设计规范GB 50909—2014 （中华人民共和国住房和城乡建设部，2014）也提到了地震动空间变化对结构的影响.地震动空间变化可以通过以下三种方法来分析：① 时间历程方法；② 随机振动方法；③ 反应谱方法.后两种方法的优点是应用了地震动的统计特征对结构响应进行分析.

对于随机振动分析，地震动输入由自功率谱密度矩阵（auto-power spectral density，缩写为APSD）确定；对于反应谱分析，地震动输入由地震反应谱确定，模态组合系数通常由结构随机振动理论推导得出（王君杰，1992；Berrah，Kausel，1992；Kiureghian，Neuenhofer，1992；王君杰等，1995）.一般情况下，随机振动法的数值积分运算量很大，基于随机振动理论的反应谱组合系数的计算也非常耗时.为减少积分时间，Harichandran和Vanmarcke （1986）提出了行波激励的解析公式，Loh和Ku （1995）基于简化的自功率谱密度函数提出了一种有效的计算方法，但采用了简单的空间相干函数.孙建梅等（2003）对Kiureghianh和Neuenhofer （1992）提出的多点输入反应谱方法进行了简化，即在相关系数的计算中忽略频率比介于0.7—1.2范围以外的部分，减少了计算量，但损失了精度.

为提高多点地震动激励下谱分析方法中相关系数的计算效率，本文拟基于对自功率谱和相关系数积分表达式的近似处理，提出相关系数的近似解析公式.

1 结构地震响应

假定在空间某点r（r＝ 1，2，···，m）处的地震激励是零均值平稳过程，其位移、速度和加速 度 函 数 分 别 记 为u（rt），u˙ （rt）和u¨（rt）， 结 构 平 稳 响 应 的 方 差 为 （王 君 杰 ，1992；Kiureghian，Neuenhofer，1992）

式中：σk为第k个结构反应量的均方差；σgr和σgs分别为第r个空间点和第s个空间点的地面运动位移的均方差；σir和σjs分别为第r个地面运动激励下结构第i振型反应的均方差和第s个地面运动激励下结构第j振型反应的均方差；参数a和b均为与结构有关的系数；m为多点地震动的数目；n为参与计算的结构的振型数目；ρgrgs，ρgrjs和ρirjs定义为

其中

式中，H（iω）为结构自由振动第i阶振型的频响函数，）为第r空间点地震加速度与第s空间点地震加速度的互功率谱密度函数，ω为圆频率.根据式（1），多点地震激励下结构反应的组合方法为

式中，SD（rωi，ζi）和SD（sωj，ζj）分别代表第r个地震激励下i阶振型最大位移和第s个地震激励下j阶振型最大位移，即平均位移反应谱.同时，ur，max＝SD（r0，ζi），us，max＝SD（s0，ζj）.

对于式（3），根据复相干函数，互功率谱密度函数Sr（sω）可以写为

式中，γ为相关因子， ρr（sω，Δrrs）为地震动的空间相干函数， θr（sω，Δrrs）为地震动的相位差，Δrrs为空间距离的向量.

国内外研究人员已经提出了多个 ρr（sω，Δrrs）表达式，例如Harichandran和Vanmarcke（1986）及Qu等（1996）提出的表达式等.Qu等（1996）提出的地震动空间相干函数模型如图1所示，其经验表达式为

图1 空间相干函数模型（Qu et al，1996）（a） ρrs-空间距离关系；（b） ρrs-频率关系Fig.1 Spatial coherence function model （Qu et al，1996）（a） ρrs-spatial distance relationship；（b） ρrs-frequency relationship

式中：a1＝1.678×10-5，a2＝1.219×10-3；b1＝-5.5×10-3，b2＝0.7674.

Harichandran和Vanmarcke（1986）提出的地震动空间相干函数模型如图2所示，其经验表达式为

图2 空间相干函数模型（Harichandran，Vanmarcke，1986）（a） ρrs-空间距离关系；（b） ρrs-频率关系Fig.2 Spatial coherence function model （Harichandran，Vanmarcke，1986）（a） ρrs-spatial distance relationship；（b） ρrs-frequency relationship

式中A， α ，K，θ， ω0和b为经验参数.对于台湾SMART-1台网第24次地震，A＝0.736，α＝0.147，K＝5210， ω0＝6.85 rad/s，b＝2.78.

本文中使用到的两个自功率谱函数经验模型见图3.第一个是Hu模型（胡聿贤，周锡元，1962，1965；王君杰，1992；王君杰，江近仁，1997），即

图3 自功率谱密度函数 S r（rω）的特征（ζg 为场地土的阻尼比，下同）（a） 胡聿贤自功率谱密度模型；（b） 克拉夫-彭津自功率谱密度模型Fig.3 The characteristic for auto-power spectral density function（ζg is the damping ratio of site soil，the same below）（a） Hu’s auto-power spectral density model；（b） Clough-Penzien’s auto-power spectral density model

第二个是克拉夫-彭津模型（Clough，Penzien，1993），即

本文在计算相关系数时，将以上两个模型分别简化为

式中，S0为APSD的强度因子， ωg为 场地土层的频率参数， ζg为场地土层的阻尼参数， ωc为胡聿贤模型在低频段的滤波参数， ωf和 ζf为Clough-Penzien模型在低频段的滤波参数.计算中取ωg＝1.5 Hz， ωc＝0.3 Hz， ωf＝0.25 Hz， ζg＝0.6， ζf＝0.4，S0＝1.

2 近似自功率谱密度函数的精度

分别比较相干系数 ρgrgs， ρgrjs和 ρirjs的近似值［ 式 （10）］ 与 精确值［ 式 （8）］ 之 间的差异，结果如图4—6所示.其中，空间相干函数使用式（6），视波速向量vapp的模取1 000 m/s.

图4 系数 ρ grgs的精确值与近似值的比较Fig.4 The comparison between exact value and approximate value of coefficient ρgrgs

从图4可以看出，式（10）给出的系数 ρgrgs几乎是精确值.图5表明：当 ωj处于低频段时，系数 ρgrjs非常大；随着 ωj的增大及r点与s点之间空间距离的增大， ρgrjs逐渐减小至零，直至可以忽略，说明式（10）可以给出具有高精度的 ρgrjs近似解答.

图5 |Δ r rs|＝300m （a）和1 000 m （b）时系数 ρ grjs 的精确值与近似值的比较（ζj ＝ 0.05）Fig.5 The comparison between exact value and approximate value of coefficient ρg r js （ζj ＝ 0.05）with |Δ r rs|＝300m （a） and 1 000 m （b）

从图6a-d可以看出：当 ωi和 ωj位于低频段时，由式（10）计算得出的相干系数 ρirjs近乎精确值， 而 近 似的 ρirjs在 ωi和 ωj较 大时会 出 现 更 大的误差；在高频段内，地震动自功率谱的值比低频段的值小得多，相比之下重要性不大.同时，系数 ρirjs的值随着 ωi和 ωj的增加而减小.以上结果表明，相干系数 ρirjs的近似值在工程实践中具有很好的准确性.

图6 系数 ρ irjs （ζi = ζj=0.05）精确值与近似值的比较Fig.6 The comparison between exact value and approximate value of coefficient ρi rj s （ζi = ζj=0.05）（a） |Δ r rs|＝300m，ωi＝0.1 Hz；（b） |Δ r rs|＝1000m，ωi＝0.1 Hz；（c） |Δ r rs|＝300m，ωi＝1.0 Hz；（d） |Δ r rs|＝1000m，ωi＝1.0 Hz；（e） |Δ r rs|＝300m ，ωi＝4.0 Hz；（f） |Δ r rs|＝1000m，ωi＝4.0 Hz

3 相关系数的近似公式

3.1 近似空间相干函数 ρr（sω，Δrrs）

如前所述，由于空间相干函数的经验表达式 ρr（sω，Δrrs）的复杂性，需要进行数值积分来计算系数 ρgrgs， ρgrjs和 ρirjs.为了得到这三个系数的近似解析表达式，引入近似的 ρr（sω，Δrrs）.

组合系数 ρgrgs与拟静态响应有关.拟静态响应是由位移贡献的，它由地震动中非常低的频率成分控制.因此，计算系数 ρgrgs的相干函数近似表达式为

式中， ωd，max是地震动中位移的APSD达到最大值时对应的圆频率.

组合系数 ρgrjs与式（1）或式（4）中的第二项耦合响应项有关，此耦合项由地震动中低频成分控制.因此，计算系数 ρgrjs的相干函数近似表达式被定义为

式中， ωgrjs取或（ ωj＋ωd，max）/2.

组合系数 ρirjs与式（1）或式（4）中的第三项动力响应项有关，动力响应项由地震动中靠近以及介于 ωi和 ωj之间的频率成分控制.因此，计算系数 ρirjs的相干函数近似表达式为

从图7—9所示结果可以看出，使用式（6）和式（8）计算得出的上述近似相干函数 ρr（sω，drs）具有很好的准确性.由此可知，本文所提方法在ρr（sω，Δrrs）的近似表达上具有良好的精度.

图7 系数 ρ grgs的精确值与近似值的比较Fig.7 The comparison between exact value and approximate value of coefficient ρgrgs

图8 |Δ rr s|＝300m （a）和1 000 m （b）时系数 ρ grjs （ζj ＝ 0.05）精确值与近似值的比较Fig.8 The comparison between exact value and approximate value of coefficient ρg r js （ζj ＝ 0.05）with |Δ rr s|＝300m （a） and 1 000 m （b）

经过上述简化后，三个相关系数 ρgrgs， ρgrjs和 ρirjs可以简化为

对于胡聿贤自功率谱密度模型有

图9 系数 ρ irjs （ζi ＝ ζj＝0.05）精确值与近似值的比较Fig.9 The comparison between exact value and approximate value of coefficient ρi rj s （ζi ＝ ζj＝0.05）（a） |Δ r rs|＝300m，ωi＝0.1 Hz；（b） |Δ r rs|＝1000m，ωi＝0.1 Hz；（c） |Δ r rs|＝300m，ωi＝1.0 Hz；（d） |Δ r rs|＝1000m，ωi＝1.0 Hz；（e） |Δ r rs|＝300m，ωi＝4.0 Hz；（f） |Δ r rs|＝1000m，ωi＝4.0 Hz

式中： R e 为复数的实部； τ为地震波从r点传播到s点的时间，，当r与s在相同位置时， τ等于零.

对于克拉夫-彭津的自功率谱密度模型则有

3.2 胡聿贤自功率谱密度模型下相关系数的解析表达式

根据复变函数的基本理论，式（16）的积分形式为

式中， R e［s（fzp）］ 为复变函数的留数，zp是函数f（z）的极点，Imzp是zp的虚部.

根据式（18）得到

如果 τ ≥0 ，则

如果 τ ＜0，则

如果 τ ≥0 ，则

如果 τ ＜0，则

3.3 Penzien-Clough 自功率谱密度下相关系数的解析表达式

式（16）的积分形式为

如果 τ ≥0 ，则

如果 τ ＜0 ，则

如果 τ ≥0 ，则

其中，z1＝ωjeiθj，z2＝ - ωje-iθj，z3＝ωfe-iθf，z4＝ - ωfeiθf；z5＝ωfeiθf，z6＝ - ωfe-iθf，z7＝ωie-iθi，z8＝-ωieiθi；这里取下半平面的单极点留数 R e［s（fz3）］ ， R e［s（fz4）］ ， R e［s（fz7）］ 和 R e［s（fz8）］.

如果 τ＜ 0，则

以上诸式中， θi＝arcsinζi， θj＝arcsinζj， θf＝arcsinζf.

4 精度和计算效率的讨论

4.1 相关系数解析表达式的精度

本节以一座斜拉桥为例进行计算，其中跨240 m，边跨117.5 m，其有限元模型示意图如图10所示.假设地震波沿纵向传播，视波速取1 000 m/s.用式（9）或式（11）的自功率谱密度函数以及式（6）的地震动空间相干模型，计算该桥梁在多点地震激励下的反应.

图10 算例桥梁的有限元模型Fig.10 Finite element model of the example bridge

计算中取前50阶振型.表1—4分别列出了桥梁不同构件中某些位置处的绝对位移、相对位移、轴力、剪力和弯矩.表中：Ln（n＝1—6）代表桥梁左侧的反应，Rn代表右侧的反应；n＝1时代表塔顶的绝对位移；n＝2时代表塔顶和塔底之间的相对位移；n＝3时代表边墩顶和梁端之间的相对位移，这项参数在落梁问题中很重要；n＝4时代表塔底的轴力；n＝5时代表塔底的剪力；n＝6时代表塔底的弯矩.

表1 Qu 等（1996）相干模型下的地震位移响应及相对误差Table 1 Seismic displacement response and its relative error under the coherence model of Qu et al （1996）

定义两个相对误差：

式中：A表示基于精确APSD模型和精确空间相干模型的响应的数值积分解；B表示基于精确APSD模型和近似空间相干模型的响应的数值积分解；C为3.2和3.3节中提出的解析解；eBA和eCA是相对误差.

表2 Harichandran 和 Vanmarcke （1986）相干模型下的地震位移响应及相对误差Table 2 Seismic displacement response and its relative error under the coherence model of Harichandran and Vanmarcke （1986）

在方法A和B中，使用了变步长的辛普森积分法.在积分极限从1 000，500，200到100 rad/s的情况下，将基于胡聿贤的功率谱密度函数与屈铁军、王君杰的空间相干函数的系数 ρgrgs， ρgrjs， ρirjs进行了比较，最终方法A采用200 rad/s.当相对积分精度分别为1 0-8，1 0-5，10-4和1 0-3时，分别计算出系数 ρgrgs， ρgrjs， ρirjs，并采用了精度为1 0-4时的积分步长.

表3 Qu 等（1996）相干模型下的地震力或弯矩响应及相对误差Table 3 Seismic force or moment response and its relative error under the coherence model of Qu et al （1996）

表4 Harichandran 和 Vanmarcke （1986）相干模型下的地震力或弯矩响应及相对误差Table 4 Seismic force or moment response and its relative error under the coherence model of Harichandran and Vanmarcke （1986）

从表1—4可以看出，相对误差非常小，这说明近似空间相干函数对结构的反应影响非常小.从3.1节中的图7—9也可以得出这个结论.

除了塔顶的绝对位移（即L1和R1）外，A和C方法的反应相对误差都小于1%.同时在地震分析中，绝对位移并不重要，而相对位移却很重要.这表明解析表达式可以与数值方法很好地吻合.

4.2 计算效率的比较

表5列出了数值方法A和解析方法C对该例桥梁的计算时间，计算中取前50阶振型.表6列出了计算中不同阶振型的计算时间.在本例中，只计算了以下反应：① 四个塔顶的绝对位移；② 四个塔顶与塔底之间的相对位移；③ 两个边墩顶部与梁端之间的相对位移；④ 四个塔底的轴力、剪力和弯矩.在本例中，使用英特尔Visual Fortran 11.1.067编写有限元程序，前50阶的特征值计算耗时50 s.

从表5可以发现，解析法求解的时间消耗只有数值求解的1/52左右.就绝对时间而言，解析法需要93—111 s，这在实际设计中是完全可以接受的.同时，数值解决方案需要5 833—6 273 s，比解析法要多耗时近98%.从表6中也可以得出同样的结论，计算时间随着所取振型阶数的增加而增加.对于大型桥梁，在计算中可能需要取数百阶振型.与用数值法相比，采用解析法计算大型结构响应可大大缩短计算用时.

表5 计算时间比较（单位：s）Table 5 Comparison of time consumption of computation （Unit：s）

表6 计算时间比较（单位：s）Table 6 Comparison of time consumption of computation （Unit：s）

5 结论

本文在提出近似空间相干函数的基础上，建立了相关系数 ρgrgs， ρgrjs和 ρirjs的解析表达式，一方面保证了计算结果的精度，另一方面采用解析法可以更快速地获得三个相关系数的积分值.同时使用了仅需要一至两个地震动参数的APSD近似表达式.通过对数值计算结果的分析，本文所提相关系数 ρgrgs， ρgrjs和 ρirjs的解析式有足够的工程精度和较高的计算效率，可为平面大尺度结构多点激励下的地震响应计算提供高效方法.