留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

滞后阻尼体系地震反应的中心差分虚初始条件法

祁文睿 潘旦光 高永涛 付相球

祁文睿, 潘旦光, 高永涛, 付相球. 滞后阻尼体系地震反应的中心差分虚初始条件法[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.10.0596
引用本文: 祁文睿, 潘旦光, 高永涛, 付相球. 滞后阻尼体系地震反应的中心差分虚初始条件法[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.10.0596
Wen-rui Qi, Dan-guang Pan, Yong-tao Gao, Xiang-qiu Fu. A CENTRAL DIFFERENCE VIRTUAL INITIAL CONDITION METHOD FOR SEISMIC RESPONSE ANALYSIS OF SYSTEMS WITH HYSTERETIC DAMPING[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.10.0596
Citation: Wen-rui Qi, Dan-guang Pan, Yong-tao Gao, Xiang-qiu Fu. A CENTRAL DIFFERENCE VIRTUAL INITIAL CONDITION METHOD FOR SEISMIC RESPONSE ANALYSIS OF SYSTEMS WITH HYSTERETIC DAMPING[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.10.0596

滞后阻尼体系地震反应的中心差分虚初始条件法

doi: 10.6052/j.issn.1000-4750.2019.10.0596
基金项目: 土木工程防灾国家重点实验室开放基金项目(SLDRCE15-01)
详细信息
    作者简介:

    祁文睿(1981−),男,江苏人,教授级高工, 博士研究生,主要从事岩土工程研究(E-mail: qwr711@163.com)

    高永涛(1962−),男,山东人,教授, 博士,博导,主要从事岩土工程和采矿工程方面的教学和研究(E-mail: gyt1962@163.com)

    付相球(1994−),男,江西人,博士研究生,主要从事土木工程防灾减灾研究(E-mail: b20180019@xs.ustb.edu.cn)

    通讯作者: 潘旦光(1974−),男,浙江人,教授,博士,博导,主要从事土木工程防灾减灾研究(E-mail: pdg@ustb.edu.cn)
  • 中图分类号: O342

A CENTRAL DIFFERENCE VIRTUAL INITIAL CONDITION METHOD FOR SEISMIC RESPONSE ANALYSIS OF SYSTEMS WITH HYSTERETIC DAMPING

  • 摘要: 针对滞后阻尼体系的直接积分法解发散问题,提出了中心差分虚初始条件法。该方法以实数初始条件为基础,建立相应的虚初始条件,以消除自由振动中的发散项,而使直接积分法的计算结果收敛到稳定解;并以有条件稳定的中心差分数值计算方法为基础,建立中心差分虚初始条件法的直接积分程序;然后,对三个不同自振频率体系在三条地震波作用下进行中心差分虚初始条件法、解析解和频域解的数值计算。计算结果表明,频域分析结果仅包含稳态解,这将导致低自振频率体系在振动的初始阶段将产生显著误差;中心差分虚初始条件法对不同工况的计算结果都是稳定的,且收敛到包括瞬态反应的解析解。
  • 图  1  输入地震波的加速度时程

    Figure  1.  Acceleration time histories of input seismic waves

    图  2  三条地震作用下位移反应(f=10 Hz, u0=v0=0)

    Figure  2.  Displacement responses under three seismic excitations (f=10 Hz, u0=v0=0)

    图  3  天津波作用下的速度和加速度反应(f=10 Hz, u0=v0=0)

    Figure  3.  Velocity and acceleration response under Tianjin wave excitation

    图  4  天津波作用下不同自振频率体系的位移反应(u0=v0=0)

    Figure  4.  Displacement response for various natural vibration frequencies under Tianjin wave excitation

    图  5  天津波作用下不同自振频率体系的位移反应

    Figure  5.  Displacement response of various natural vibration frequencies under Tianjin wave excitation

    图  6  天津波作用下体系的速度和加速度反应(f=10 Hz)

    Figure  6.  Velocity and acceleration response under Tianjin wave excitation

    表  1  任意荷载下中心差分虚初始条件法计算步骤

    Table  1.   Calculation procedure of central differential virtual initial condition method under arbitrary load

    1. 准备工作
    1.1 对加速度时程a(t)根据式(28)和式(29)计算系数
    1.2 计算虚初始条件,形成总初始条件$u_0^t = u_0^{} + u_0^{v1} + u_0^{v{\rm{4}}}$, $v_0^t = v_0^{} + v_0^{v1} + v_0^{v{\rm{4}}}$
    1.3 $\ddot u_0^t = \left[ { {f_0} - k(1 + {\rm{i}}\eta )u_0^t} \right]/m$$u_{ - 1}^t = u_0^t - \Delta t\dot u_0^t - {(\Delta t)^2}\ddot u_0^t{\rm{/2}}$
    1.4 选择Δt
    1.5 $\hat k = m{\rm{/}}\Delta {t^2}$
    2. 每个荷载步计算
    2.1 ${\hat f_{n + 1} } = {f_n} - \left[ {k(1 + {\rm{i}}\eta ) - \dfrac{ {2m} }{ {\Delta {t^2} } } } \right]u_n^t - \dfrac{m}{ {\Delta {t^2} } }u_{n - 1}^t$
    2.2 ${u_{n + 1}} = {\hat f_{n + 1}}/\hat k$
    2.3 ${\ddot u_{n + 1} } = \dfrac{1}{m}[ { {f_{n + 1} } - k(1 + {\rm{i}}\eta ){u_{n + 1} } } ]$${\dot u_{n{\rm{ + 1} } } } = \dfrac{ {\Delta t} }{2}{\ddot u_{n + {\rm{1} } } }{\rm{ + } }\dfrac{ { {u_{n + {\rm{1} } } } - {u_n} } }{ {\Delta t} }$
    2.4 $u_{n{\rm{ + } }1}^{v1} = - {\rm{i}}\dfrac{ { {\rm{Re} } ({ {\dot u}_{n + 1} }) + \omega \alpha {\rm{Re} } ({u_{n + 1} })} }{ {\omega \mu } }$
    2.5 $u_{n + 1}^t = {\rm{Re}} ({u_{n{\rm{ + 1}}}}) + u_{n + 1}^{v1} + u_0^{v{\rm{4}}}$
    3. 用n+1代替n,重复2.1~2.5进行下一荷载步的计算。
    下载: 导出CSV

    表  2  地震波主要参数

    Table  2.   The main parameters of ground motions

    地震波时间地震峰值加速度/(m/s2)
    El Centro1940.5.18California3.417
    迁安波1976.8.31唐山余震0.974
    天津波1976.11.25唐山余震1.458
    下载: 导出CSV

    表  3  本文方法与频域解的峰值相对误差

    Table  3.   Peak relative errors of the proposed method and frequency domain solution

    地震波自振频率/Hz${e_u}$/(%)${e_{\dot u}}$/(%)${e_{\ddot u}}$(%)
    本文方法频域解本文方法频域解本文方法频域解
    El Centro波101.460.00020.2101.370
    10.640.00290.540.020.080.0218
    0.10.6917.130.3637.960.391.74
    迁安波100.0800.00400.270
    11.230.01041.540.0060.110.0033
    0.13.3618.080.932.570.120.10
    天津波100.760.00021.9302.210
    10.190.420.770.420.060.25
    0.10.536.080.320.170.180.07
    下载: 导出CSV

    表  4  不同算法计算时间比较

    Table  4.   Comparison of calculation time of various methods /s

    地震波解析解本文方法频域解
    天津波1.400.0030.0015
    下载: 导出CSV
  • [1] Nashif A D, Jones D I G, Henderson J P. Vibration damping[M]. New York: John Wiley & Sons, 1985.
    [2] 赵密, 杜修力. 时间卷积的局部高阶弹簧-阻尼-质量模型[J]. 工程力学, 2009, 26(5): 8 − 18.

    Zhao Mi, Du Xiuli. High-order spring-dashpot-mass model for localizing time convolution [J]. Engineering Mechanics, 2009, 26(5): 8 − 18. (in Chinese)
    [3] Chopra A K. Dynamics of structures: Theory and applications to earthquake engineering[M]. New Jersey: Englewood Cliffs, Prentice-Hall, 1995.
    [4] Housner G W, Bergman L A, Caughey T K, et al. Structural control: Past, present, and future [J]. Journal of Engineering Mechanics, 1997, 123(9): 897 − 971. doi:  10.1061/(ASCE)0733-9399(1997)123:9(897)
    [5] 隋䶮, 薛建阳, 董金爽, 等. 附设粘滞阻尼器的混凝土仿古建筑梁-柱节点恢复力模型试验研究[J]. 工程力学, 2019, 36(增刊 1): 44 − 53.

    Sui Yan, Xue Jianyang, Dong Jinshuang, et al. The experimental research on the restoring force model of lintel-column joints of concrete archaizing buildings with viscous dampers [J]. Engineering Mechanics, 2019, 36(Suppl 1): 44 − 53. (in Chinese)
    [6] Zhou Li, Su Yisheng. Cyclic loading test on beam-to-column connections connecting SRRAC beams to RACFST columns [J]. International Journal of Civil Engineering, 2018, 16(11): 1533 − 1548. doi:  10.1007/s40999-018-0288-x
    [7] 李瑞山, 陈龙伟, 袁晓铭等. 荷载频率对动模量阻尼比影响的试验研究[J]. 岩土工程学报, 2017, 39(1): 71 − 80.

    Li Ruishan, Chen Longwei, Yuan Xiaoming, et al. Experimental study on influences of different loading frequencies on dynamic modulus and damping ratio [J]. Chinese Journal of Geotechnical Engineering, 2017, 39(1): 71 − 80. (in Chinese)
    [8] 巴振宁, 喻志颖, 梁建文. 山岭隧洞在平面P-SV波入射下的地震动力响应[J]. 地震工程与工程振动, 2018, 38(5): 52 − 60.

    Ba Zhenning, Yu Zhiying, Liang Jianwen. Dynamic response of a mountain tunnel under plane P-SV waves [J]. Earthquake Engineering and Engineering Vibration, 2018, 38(5): 52 − 60. (in Chinese)
    [9] Chen X, Chen HL, Hu XL. Damping predication of sandwich structures by order-reduction-iteration approach [J]. Journal of Sound & Vibration, 1999, 222(5): 803 − 81.
    [10] Sheng Meiping, Guo Zhiwei, Qin Qi, et al. Vibration characteristics of a sandwich plate with viscoelastic periodic cores [J]. Composite Structures, 2018, 206: 54 − 69.
    [11] Idriss I M, Seed H B. Seismic response of horizontal soil layers [J]. Journal of the Soil Mechanics & Foundations Division, 1968, 94: 1003 − 1031.
    [12] Nampally S, Padhy S, Trupti S, et al. Evaluation of site effects on ground motions based on equivalent linear site response analysis and liquefaction potential in Chennai, south India [J]. Journal of Seismology, 2018, 22(4): 1075 − 1093. doi:  10.1007/s10950-018-9751-z
    [13] 王笃国, 赵成刚. 考虑频率参数协调的频率相关等效线性化方法[J]. 工程力学, 2019, 36(9): 169 − 179.

    Wang Duguo, Zhao Chenggang. Frequency-dependent equivalent linear method for seismic site response considering the compatibility of frequency parameters [J]. Engineering Mechanics, 2019, 36(9): 169 − 179. (in Chinese)
    [14] 朱镜清. 关于复阻尼理论的两个基本问题[J]. 固体力学学报, 1992, 13(2): 113 − 118.

    Zhu Jingqing. On the two basic problems of complex damping theory [J]. Acta Mechanica Solida Sinica, 1992, 13(2): 113 − 118. (in Chinese)
    [15] 朱镜清, 朱敏. 复阻尼地震反应谱的计算方法及其它[J]. 地震工程与工程振动, 2000, 20(2): 19 − 23. doi:  10.3969/j.issn.1000-1301.2000.02.004

    Zhu Jingqing, Zhu Min. Calculation of complex damping response spectra from earthquake records [J]. Earthquake Engineering and Engineering vibration, 2000, 20(2): 19 − 23. (in Chinese) doi:  10.3969/j.issn.1000-1301.2000.02.004
    [16] 何钟怡. 复本构理论中的对偶原则[J]. 固体力学学报, 1994, 15(2): 177 − 180.

    He Zhongyi. The dual principle in theory of complex constitutive equations [J]. Acta Mechanica Solida Sinica, 1994, 15(2): 177 − 180. (in Chinese)
    [17] 朱敏, 朱镜清. 逐步积分法求解复阻尼结构运动方程的稳定性问题[J]. 地震工程与工程振动, 2001, 21(4): 59 − 62.

    Zhu Min, Zhu Jingqing. Studies on stability of step -by -step methods under complex damping conditions [J]. Earthquake Engineering and Engineering Vibration, 2001, 21(4): 59 − 62. (in Chinese)
    [18] 张辉东, 王元丰. 复阻尼模型结构地震时程响应研究[J]. 工程力学, 2010, 27(1): 109 − 115.

    Zhang Huidong, Wang Yuanfeng. Study on seismic time-history response of structures with complex damping [J]. Engineering Mechanics, 2010, 27(1): 109 − 115. (in Chinese)
    [19] 潘玉华, 王元丰. 复阻尼结构动力方程的高斯精细时程积分法[J]. 工程力学, 2012, 29(2): 16 − 20.

    Pan Yuhua, Wang Yuanfeng. Gauss precise time-integration of complex damping vibration systems [J]. Engineering Mechanics, 2012, 29(2): 16 − 20. (in Chinese)
    [20] 何钟怡, 廖振鹏, 王小华. 关于复阻尼理论的几点注记[J]. 地震工程与工程振动, 2002, 22(1): 1 − 6. doi:  10.3969/j.issn.1000-1301.2002.01.001

    He Zhongyi, Liao Zhenpeng, Wang Xiaohua. Some notes on theory of complex damping [J]. Earthquake Engineering and Engineering Vibration, 2002, 22(1): 1 − 6. (in Chinese) doi:  10.3969/j.issn.1000-1301.2002.01.001
    [21] Henwood D J. Approximating the hysteretic damping matrix by A viscous matrix for modelling in the time domain [J]. Journal of Sound and Vibration, 2002, 254(3): 575 − 593. doi:  10.1006/jsvi.2001.4136
    [22] 楼梦麟, 潘旦光. 滞后阻尼在土层时域分析中的应用[J]. 同济大学学报, 2004, 32(3): 281 − 285.

    Lou Menglin, Pan Danguang. Hysteretic damping application in time domain analysis of soil layer [J]. Journal of Tongji University, 2004, 32(3): 281 − 285. (in Chinese)
    [23] 孙攀旭, 杨红, 赵雯桐, 刘庆林. 基于复阻尼模型的时域数值计算方法[J]. 地震工程与工程振动, 2019, 39(2): 203 − 211.

    Sun Panxu, Yang Hong, Zhao Wentong, Liu Qinglin. The time-domain numerical calculation method based on complex damping model [J]. Earthquake Engineering and Engineering Vibration, 2019, 39(2): 203 − 211. (in Chinese)
    [24] 孙攀旭, 杨红, 赵雯桐, 等. 基于滞变阻尼模型的时域计算方法[J]. 工程力学, 2019, 36(6): 13 − 20.

    Sun Panxu, Yang Hong, Zhao Wentong, et al. Time domain calculation method based on hysteretic damping model [J]. Engineering Mechanics, 2019, 36(6): 13 − 20. (in Chinese)
    [25] Ribeiro A M R, Maia N M M, Silva J M M. Free and forced vibration with viscous and hysteretic damping: A different perspective[C]. 5th International Conference on Mechanics and Materials in Design, Porto, Portugal, 24-26, July, 2006. Chapter V: 1-10.
    [26] Pan Danguang, Fu Xiangqiu, Qi Wenrui. The direct integration method with virtual initial conditions on the free and forced vibration of a system with hysteretic Damping [J]. Applied Science, 2019, 9(18): 3707. doi:  10.3390/app9183707
    [27] 朱敏, 朱镜清. 复阻尼振动方程频域解法中有关问题的研究[J]. 世界地震工程, 2004, 20(1): 23 − 28. doi:  10.3969/j.issn.1007-6069.2004.01.004

    Zhu Min, Zhu Jingqing. Some problem in frequency domain solution of complex damping system [J]. World Earthquake Engineering, 2004, 20(1): 23 − 28. (in Chinese) doi:  10.3969/j.issn.1007-6069.2004.01.004
  • [1] 申大为, 徐明, 刘鹏飞.  加筋土整体式桥地震反应研究 . 工程力学, doi: 10.6052/j.issn.1000-4750.2017.07.0538
    [2] 杜修力, 李洋, 赵密, 许成顺, 路德春.  下卧刚性基岩条件下场地土-结构体系地震反应分析方法研究 . 工程力学, doi: 10.6052/j.issn.1000-4750.2015.09.0801
    [3] 孙颖, 陈天海, 卓卫东, 谷音, 许智星.  长周期地震动作用下隔震连续梁桥地震反应特性研究 . 工程力学, doi: 10.6052/j.issn.1000-4750.2015.05.S054
    [4] 陆磊, 李敏, 甘宽.  库伦摩擦对二元翼面气动弹性特性的影响 . 工程力学, doi: 10.6052/j.issn.1000-4750.2013.11.1086
    [5] 潘旦光, 高莉莉.  Rayleigh阻尼系数解法比较及对结构地震反应影响 . 工程力学, doi: 10.6052/j.issn.1000-4750.2013.12.1190
    [6] 潘旦光.  地震反应分析中Rayleigh阻尼系数的优化解 . 工程力学, doi: 10.6052/j.issn.1000-4750.2012.05.0366
    [7] 孙广俊, 李鸿晶, 王通.  考虑桥墩及支座影响的梁桥竖向地震反应分析 . 工程力学, doi: 10.6052/j.issn.1000-4750.2011.02.0079
    [8] 霍林生, 李宏男.  调液阻尼器对偏心结构扭转耦联振动控制的研究 . 工程力学,
    [9] 郭长青, 刘红涛, 王晓锋, 张楚汉.  输流管道在分布随从力作用下的振动和稳定性 . 工程力学,
    [10] 桂国庆, 李永华, 李思明.  复杂高层连体结构多点激励地震反应分析 . 工程力学,
    [11] 潘旦光, 楼梦麟.  层状土层随机地震反应分析的近似解法 . 工程力学,
    [12] 邢德进, 李忠献.  应用SMA智能阻尼器的结构模糊控制 . 工程力学,
    [13] 王倩颖, 吴 斌, 欧进萍.  考虑作动器时滞及其补偿的实时子结构实验稳定性分析 . 工程力学,
    [14] 陈兴冲, 高 峰, 吴少海.  冻土层对桥梁地震反应的影响 . 工程力学,
    [15] 郭泽英, 李青宁, 张守军.  结构地震反应分析的一种新精细积分法 . 工程力学,
    [16] 霍林生, 李宏男.  环形调液阻尼器(CTLCD)对结构平移-扭转耦联振动控制的参数研究 . 工程力学,
    [17] 颜桂云, 孙炳楠, 陆鸣.  MR阻尼器对建筑结构地震反应的半主动H_∞控制 . 工程力学,
    [18] 杜修力, 王进廷.  阻尼弹性结构动力计算的显式差分法 . 工程力学,
    [19] 许斌, 李黎.  建筑结构地震反应的模糊开环控制 . 工程力学,
    [20] 栾茂田, 林皋.  场地地震反应一维非线性计算模型 . 工程力学,
  • 加载中
计量
  • 文章访问数:  50
  • HTML全文浏览量:  28
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-10-15
  • 修回日期:  2020-01-13
  • 网络出版日期:  2020-05-27

滞后阻尼体系地震反应的中心差分虚初始条件法

doi: 10.6052/j.issn.1000-4750.2019.10.0596
    基金项目:  土木工程防灾国家重点实验室开放基金项目(SLDRCE15-01)
    作者简介:

    祁文睿(1981−),男,江苏人,教授级高工, 博士研究生,主要从事岩土工程研究(E-mail: qwr711@163.com)

    高永涛(1962−),男,山东人,教授, 博士,博导,主要从事岩土工程和采矿工程方面的教学和研究(E-mail: gyt1962@163.com)

    付相球(1994−),男,江西人,博士研究生,主要从事土木工程防灾减灾研究(E-mail: b20180019@xs.ustb.edu.cn)

    通讯作者: 潘旦光(1974−),男,浙江人,教授,博士,博导,主要从事土木工程防灾减灾研究(E-mail: pdg@ustb.edu.cn)
  • 中图分类号: O342

摘要: 针对滞后阻尼体系的直接积分法解发散问题,提出了中心差分虚初始条件法。该方法以实数初始条件为基础,建立相应的虚初始条件,以消除自由振动中的发散项,而使直接积分法的计算结果收敛到稳定解;并以有条件稳定的中心差分数值计算方法为基础,建立中心差分虚初始条件法的直接积分程序;然后,对三个不同自振频率体系在三条地震波作用下进行中心差分虚初始条件法、解析解和频域解的数值计算。计算结果表明,频域分析结果仅包含稳态解,这将导致低自振频率体系在振动的初始阶段将产生显著误差;中心差分虚初始条件法对不同工况的计算结果都是稳定的,且收敛到包括瞬态反应的解析解。

English Abstract

祁文睿, 潘旦光, 高永涛, 付相球. 滞后阻尼体系地震反应的中心差分虚初始条件法[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.10.0596
引用本文: 祁文睿, 潘旦光, 高永涛, 付相球. 滞后阻尼体系地震反应的中心差分虚初始条件法[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.10.0596
Wen-rui Qi, Dan-guang Pan, Yong-tao Gao, Xiang-qiu Fu. A CENTRAL DIFFERENCE VIRTUAL INITIAL CONDITION METHOD FOR SEISMIC RESPONSE ANALYSIS OF SYSTEMS WITH HYSTERETIC DAMPING[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.10.0596
Citation: Wen-rui Qi, Dan-guang Pan, Yong-tao Gao, Xiang-qiu Fu. A CENTRAL DIFFERENCE VIRTUAL INITIAL CONDITION METHOD FOR SEISMIC RESPONSE ANALYSIS OF SYSTEMS WITH HYSTERETIC DAMPING[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.10.0596
  • 阻尼为系统的固有特征,对结构的动力反应有显著影响。根据材料的耗能特性,可建立相应的阻尼模型[1]。对于土木工程而言,粘滞阻尼和滞后阻尼是两个广泛使用的模型。粘滞阻尼力与速度成正比,具有计算简便的特点而成为广泛使用的模型[2],如Rayleigh阻尼模型和速度相关型的阻尼器[3-5]。但是,粘滞阻尼的耗能和频率相关。已有实验表明,很多材料的耗能与频率无关,如型钢混凝土梁[6]、土[7-8]和粘弹性夹层梁和板[9-10]等,此时采用阻尼力与位移成正比相位差π/2的滞后阻尼模型更符合实验结果。滞后阻尼将导致复刚度运动方程,而常采用频域方法进行求解。频域计算方法以Fourier变换为基础,理论上仅适用于线弹性体系,对于非线性系统,常采用等效线性化进行计算[11-13]。为进行真非线性计算,需要在时域中进行直接积分法进行计算。

    为在时域中进行复本构运动方程的求解,朱镜清[14-15]、何钟怡[16]等根据对偶原则,建立与实荷载对应的对偶虚荷载,完善了滞后阻尼体系的输入理论。但是,采用时域直接积分法计算滞后阻尼体系动力反应,即使是无条件稳定的直接积分方法也易于出现发散的结果[17-19]。微分方程数值解不稳定的原因有两种:一种是数值方法的不稳定;另一种是由于方程本身具有发散解。复阻尼运动方程的逐步积分法的发散解是由后一种原因造成的[17, 20]。为使复阻尼运动方程的逐步积分法计算结果稳定,一种常用的方法是将滞后阻尼等效为近似的粘滞阻尼[21-22],但计算结果的误差较大。

    在直接积分法得到滞后阻尼体系稳定解方面,孙攀旭等分别基于激励插值方法[23]和滞变阻尼时域理论[24],建立了稳定的直接积分计算方法。事实上,滞后阻尼体系的特征值为互为相反数的复数特征对[25],必然有一个特征值的实部是正的,由此导致强迫荷载下补解中的一项是没有物理意义的发散项,在解析解时人工删除发散项而使计算结果稳定。而直接积分法的计算结果是包含补解的,且滞后阻尼体系直接积分法是收敛到包含发散项补解[17]而导致计算结果不稳定。Pan等[26]首先提出虚初始条件的概念,然后基于无条件稳定的Newmark法建立使滞后阻尼直接积分解不出现发散项的计算方法。本文进一步建立了恒载下的虚初始条件,以及基于有条件稳定的中心差分法建立滞后阻尼的中心差分虚初始条件法。在此基础上,通过算例验证验证所提方法的稳定性、计算精度和计算效率。

    • 单自由度滞后阻尼体系的运动方程为[25]

      $$m\ddot u + (1 + {\rm{i}}\eta )ku = 0{\kern 1pt} $$ (1)

      方程的初始条件为:

      $$u(0) = {u_0}, \;\dot u(0) = {v_0}$$ (2)

      式中:mk分别为系统的质量和刚度;η为复阻尼系数;u0v0分别为初始位移和初始速度,是实数。方程(1)的不稳定的通解为:

      $$u = {C_1}{{\rm{e}}^{st}} + {C_2}{{\rm{e}}^{ - st}}$$ (3)

      稳定的通解为[25]

      $$u = C{{\rm{e}}^{ - st}}$$ (4)

      式中:$s \!=\! \omega (\alpha - {\rm{i}}\mu )$$\omega\! =\!\! \sqrt {k/m} $$\alpha \!=\! \sqrt {\dfrac{{ \!-\! 1\! +\!\! \sqrt {1\! +\!\! {\eta ^2}} }}{2}} $$\mu = \sqrt {\dfrac{{1 + \sqrt {1 + {\eta ^2}} }}{2}} $${C_1} = \dfrac{{s{u_0} + {v_0}}}{{2s}}$${C_2} = \dfrac{{s{u_0} - {v_0}}}{{2s}}$$ C =$${u_0} - {\rm{i}}\dfrac{{{v_0} + \omega \alpha {u_0}}}{{\omega \mu }}$。由于s的实部大于零,式(3)中右边第一项${C_1}{{\rm{e}}^{st}}$将产生无物理意义的发散解。在解析解的求解过程中,可以人为去除发散项而得到稳定解,但直接积分法将收敛到式(3)而产生发散解。因此,为使直接积分法稳定的关键是使式(3)和式(4)相同。为此,引入虚初始条件$u_0^{v1}$$v_0^{v{\rm{1}}}$达到这个目的[26]。在引入虚初始条件后,总初始条件为:

      $$u_0^t = u_0^{} + u_0^{v1}, \;\;v_0^t = v_0^{} + v_0^{v1}$$ (5)

      以式(5)为初始条件,并使式(3)等于式(4),可得:

      $$\frac{{s({u_0} + u_0^{v1}) + ({v_0} + v_0^{v1})}}{{2s}} = 0$$ (6)
      $$\frac{{s({u_0} + u_0^{v1}) - ({v_0} + v_0^{v1})}}{{2s}} = C$$ (7)

      求解式(6)和式(7)的代数方程组可得:

      $$u_0^{v1} = - {\rm{i}}\frac{{{v_0} + \omega \alpha {u_0}}}{{\omega \mu }}\quad\quad\quad\quad\quad$$ (8)
      $$ v_0^{v1} = {\rm{i}}\left[\frac{\alpha }{\mu }({v_0} + \omega \alpha {u_0}) + \omega \mu {u_0}\right] $$ (9)

      式(8)和式(9)表明,虚初始条件$u_0^{v1}$$v_0^{v{\rm{1}}}$为纯虚数。因此,对于总初始条件的实部可解释为初始条件中可观察和测量的部分,与此同时还存在与实初始条件相对应的虚部。在这样的复初始条件下,方程(1)的通解为式(4)的稳定解而无需人为干涉。

    • 在简谐荷载$A{{\rm{e}}^{{\rm{i}}\theta t}}$作用下,单自由度体系的运动方程为:

      $$m\ddot u + (1 + {\rm{i}}\eta )ku = A{{\rm{e}}^{{\rm{i}}\theta t}}$$ (10)

      式中:A为简谐荷载的幅值;θ为迫振频率。方程(10)不稳定的通解为:

      $$u = {D_{\rm{1}}}{{\rm{e}}^{st}} + {D_2}{{\rm{e}}^{ - st}} + X{{\rm{e}}^{{\rm{i}}\theta t}}$$ (11)

      稳定的通解为:

      $$u = {C_{\rm{f}}}{{\rm{e}}^{ - st}} + X{{\rm{e}}^{{\rm{i}}\theta t}}$$ (12)

      式中:$X = \dfrac{A}{{ - {\theta ^2}m + (1 + i\eta )k}}$为稳态解的幅值;D1D2Cf为常数,与初始条件有关。当$u(0) = {\rm{0}}$$\dot u(0)\! \!=\!\! {\rm{0}}$时,${D_1}\!\! = \!\!\dfrac{{ \!-\! s\! -\! {\rm{i}}\theta }}{{2s}}X$${D_2} \!\!=\! \!\dfrac{{{\rm{i}}\theta X \!\!-\!\! sX}}{{2s}}$${C_{\rm{f}}}\! =\!- {\rm{Re}} (X)-$${\rm{i}}\dfrac{{\theta {\rm{Im}} (X) - \alpha \omega {\rm{Re}} (X)}}{{\mu \omega }}$,Re表示取实部,Im表示取虚步。直接积分法解收敛到式(11)而导致计算结果发散。为使式(11)等于式(12),引入虚初始位移$u_0^{v2}$和虚初始速度$v_0^{v2}$。并令与虚初始条件$u_0^{v2}$$v_0^{v2}$对应的D1=0和D2=Cf,即:

      $$\frac{{su_0^{v2} + v_0^{v2} - sX - {\rm{i}}\theta X}}{{2s}}{\rm{ = 0}}$$ (13)
      $$\frac{{su_0^{v2} - v_0^{v2} - sX + {\rm{i}}\theta X}}{{2s}}{\rm{ = }}{C_{\rm{f}}}$$ (14)

      由此可得:

      $$u_0^{v2} = {\rm{i}}\frac{{(\mu \omega - \theta ){\rm{Im}}(X) + \alpha \omega {\rm{Re}} (X)}}{{\mu \omega }}\quad\quad\quad\;\;$$ (15)
      $$v_0^{v2} = {\rm{i}}\frac{{\alpha \theta {\rm{Im}} (X) + (\mu \theta - {\alpha ^2}\omega - {\mu ^2}\omega ){\rm{Re}} (X)}}{\mu }$$ (16)

      $u_0^{v2}$$v_0^{v2}$也是纯虚数。式(12)中的${C_{\rm{f}}}{{\rm{e}}^{ - st}}$是由荷载引起的伴生自由振动,与初始条件无关。虚初始条件$u_0^{v2}$$v_0^{v2}$也是与荷载相关的纯虚数,不改变实部可测量部分,但可使收敛到式(11)的直接积分法稳定。

      式(15)和式(16)是关于零初始条件下的虚初始条件。当$u(0) = {u_0}$, $\dot u(0) = {v_0}$时,由叠加原理可得简谐荷载$A{{\rm{e}}^{{\rm{i}}\theta t}}$下稳定的通解为:

      $$u = C{{\rm{e}}^{ - st}} + {C_{\rm{f}}}{{\rm{e}}^{ - st}} + X{{\rm{e}}^{{\rm{i}}\theta t}}$$ (17)

      而不稳定的通解为:

      $$u = {C_1}{{\rm{e}}^{st}} + {C_2}{{\rm{e}}^{ - st}} + {D_{\rm{1}}}{{\rm{e}}^{st}} + {D_{\rm{2}}}{{\rm{e}}^{ - st}} + X{{\rm{e}}^{{\rm{i}}\theta t}}$$ (18)

      同理引入虚初始位移$u_0^v$和虚初始速度$v_0^v$。由式(8)、式(9)、式(15)和式(16)可知,当$u_0^v = u_0^{v1} + u_0^{v{\rm{2}}}$$v_0^v = v_0^{v1} + v_0^{v{\rm{2}}}$时,式(17)和式(18)相同,从而使收敛到式(18)的直接积分法解也稳定。即非零初始条件下,简谐荷载作用的总初始条件为:

      $$u_0^t = u_0^{} + u_0^{v1} + u_0^{v{\rm{2}}}, \;v_0^t = v_0^{} + v_0^{v1} + v_0^{v{\rm{2}}}$$ (19)
    • 在恒载${A_0}(1 + {\rm{i}}\eta )/2$作用下,单自由度体系的运动方程为:

      $$m\ddot u + (1 + {\rm{i}}\eta )ku = {A_0}(1 + {\rm{i}}\eta )/2$$ (20)

      方程(20)不稳定的通解为:

      $$u = {E_1}{{\rm{e}}^{st}} + {E_2}{{\rm{e}}^{ - st}} + {A_0}/2k$$ (21)

      $u(0) = {\rm{0}}$$\dot u(0) = {\rm{0}}$时,可得E1E2为:

      $${E_1} = - {A_0}{\rm{/}}4k, \;{E_2} = - {A_0}{\rm{/}}4k$$ (22)

      零初始条件下稳定的通解为:

      $$u = {C_{{\rm{f}}0}}{{\rm{e}}^{ - st}} + {A_0}/2k$$ (23)

      式中,$ {C_{\rm{f}{0}}} = - {A_0}{\rm{/}}2k + {\rm{i}}\alpha {A_0}{\rm{/}}2\mu k $。为使式(21)等于式(23),引入虚初始位移$u_0^{v3}$和虚初始速度$v_0^{v3}$。在虚初始条件$u_0^{v3}$$v_0^{v3}$下,可得E1E2的解为:

      $$\left\{ {\begin{array}{*{20}{c}} {{E_1} = \dfrac{{su_0^{v3} + v_0^{v3} - s{A_0}/2k}}{{2s}}} \\ {{E_2} = \dfrac{{su_0^{v3} - v_0^{v3} - s{A_0}/2k}}{{2s}}} \end{array}} \right.$$ (24)

      令式(24)中的${E_1} = 0$${E_2} = {C_{{\rm{f}}0}}$可得:

      $$u_0^{v3} = {\rm{i}}\frac{{\alpha {A_0}}}{{2k\mu }}\quad\quad\quad\quad\quad$$ (25)
      $$v_0^{v3} = {\rm{i}}\frac{{( - {\alpha ^2}\omega - {\mu ^2}\omega ){A_0}}}{{2k\mu }}$$ (26)

      如果将恒载看成激振频率为零的简谐荷载,且稳态解的幅值$X = {A_0}/2k$,则式(25)和式(26)可以直接由式(15)和式(16)计算得到。因此,非零初始条件$u(0) = {u_0}$, $\dot u(0) = {v_0}$下,恒载作用的总初始条件的计算方法与式(19)相同。

    • 若已知地面运动加速度为a(t),将a(t)采用Fourier级数展开:

      $$\begin{split} a(t) = &\frac{{{A_0}}}{2} + \sum\limits_{j = 1}^{N/2 - 1} {( {{A_j}\cos {\theta _j}t + {B_j}\sin {\theta _j}t} )} + \\&\frac{{{A_{N/2}}}}{2}\cos {\theta _{N/2}}t \end{split}\qquad$$ (27)
      $${A_j} = \frac{2}{N}\sum\limits_{l = 0}^{N - 1} {a({t_l})} \cos ({\theta _j}{t_l}),j = 0,1,2,\cdots,\frac{N}{{\rm{2}}}\quad\quad$$ (28)
      $${B_j} = \frac{2}{N}\sum\limits_{l = 0}^{N - 1} {a({t_l})} \sin ({\theta _j}{t_l}),j = 1,2,\cdots,\frac{N}{{\rm{2}}} - 1\quad\;\;$$ (29)

      式中,${\theta _j} = 2\pi j{\rm{/}}N\Delta {t_{\rm{s}}}$(j=0,1,2,…, N/2), $\Delta {t_{\rm{s}}}$是加速度时程的采样时间间隔,${t_l} = l\Delta {t_{\rm{s}}}$(l=0,1,2,…,N-1)。根据对偶荷载原则,可构造a(t)的对偶函数b(t)[15]

      $$b(t) = \eta \frac{{{A_0}}}{2} \!+\! \sum\limits_{j = 1}^{N/2 - 1} {( {{A_j}\sin {\theta _j}t \!-\! {B_j}\cos {\theta _j}t} )} \!+\! \frac{{{A_{N/2}}}}{2}\sin {\theta _{N/2}}t$$ (30)

      则地震作用下单自由度体系的运动方程为:

      $$m\ddot u + (1 + {\rm{i}}\eta )ku = f(t)$$ (31)

      式中:$f(t) = - m[a(t) + {\rm{i}}b(t)]$。对于$u(0) = {\rm{0}}$$\dot u(0) = {\rm{0}}$,对f(t)中每一项简谐荷载对应的虚初始条件求和,可得地震作用的虚初始条件为:

      $$u_0^{v4} = \sum\limits_{j = 0}^{N/2} {( {{C_{{\rm{f}}j}} + {X_j}} )} \quad$$ (32)
      $$v_0^{v4} = \sum\limits_{j = 0}^{N/2} {( { - s{C_{{\rm{f}}j}} + {\rm{i}}{\theta _j}{X_j}} )} $$ (33)

      式中,${C_{{\rm{f}}j}} = - {\rm{Re}} ({X_j}) -{\rm{i}}\dfrac{{{\theta _j}{\rm{Im}} ({X_j}) - \alpha \omega {\rm{Re}} ({X_j})}}{{\mu \omega }}$(j=0, 1, 2,···, N/2),${X_0} = \dfrac{{m{A_0}}}{{2k}}$${X_j} = \dfrac{{m({A_j} - {\rm{i}}{B_j})}}{{ - \theta _j^2m + (1 + i\eta )k}}$(j=1, 2, ···, N/2-1),${X_{N/2}} = \dfrac{{m{A_{N/2}}}}{{2( - \theta _{N/2}^2m + (1 + {\rm{i}}\eta )k)}}$

      因此,非零初始条件$u(0) = {u_0}$, $\dot u(0) = {v_0}$下,由叠加原理可得体系的总初始条件为:

      $$u_0^t = u_0^{} + u_0^{v1} + u_0^{v{\rm{4}}}, \;v_0^t = v_0^{} + v_0^{v1} + v_0^{v{\rm{4}}}$$ (34)
    • 若已知tn时刻及以前时刻的反应已知的情况下,根据中心差分法,求解tn+1时刻的位移方程为:

      $$\hat k{u_{n + 1}} = {\hat f_{n + 1}}$$ (35)
      $${\hat f_{n + 1}} = {f_n} - \left[ {k(1 + {\rm{i}}\eta ) - \frac{{2m}}{{\Delta {t^2}}}} \right]{u_n} - \frac{m}{{\Delta {t^2}}}{u_{n - 1}}$$ (36)

      式中:$\Delta t$为计算时间步长;$\hat k = m{\rm{/}}\Delta {t^2}$。虚初始条件计算方法,仅需根据实初始条件计算相应的虚部作为初始条件,其他计算过程与粘滞阻尼体系的中心差分法[3]相同。同时,为避免数值积分过程中舍入误差的影响,对每个荷载步计算完成后,由平衡条件进行速度和加速度的计算,可提高计算精度。即由式(35)得到tn+1时刻的位移后,根据tn+1时刻的平衡条件可得${\ddot u_{n + 1}}$为:

      $${\ddot u_{n + 1}} = \frac{1}{m}\left[ {{f_{n + 1}} - k(1 + {\rm{i}}\eta ){u_{n + 1}}} \right]$$ (37)

      为计算${\dot u_{n + 1}}$,可根据tn+1时刻的差分条件:

      $${\dot u_{n{\rm{ + 1}}}} = \frac{{{u_{n + {\rm{2}}}} - {u_n}}}{{2\Delta t}}\quad\quad\quad\;$$ (38)
      $${\ddot u_{n{\rm{ + 1}}}} = \frac{{{u_{n + {\rm{2}}}} - 2{u_{n{\rm{ + 1}}}} + {u_n}}}{{\Delta {t^2}}}$$ (39)

      将式(38)、式(39)消去${u_{n + {\rm{2}}}}$,可得:

      $${\dot u_{n{\rm{ + 1}}}} = \frac{{\Delta t}}{2}{\ddot u_{n + {\rm{1}}}}{\rm{ + }}\frac{{{u_{n + {\rm{1}}}} - {u_n}}}{{\Delta t}}$$ (40)

      由此可得tn+1时刻的位移虚部修正解

      $$u_{n{\rm{ + }}1}^{v1} = - {\rm{i}}\frac{{{\rm{Re}} ({{\dot u}_{n + 1}}) + \omega \alpha {\rm{Re}} ({u_{n + 1}})}}{{\omega \mu }}$$ (41)

      tn+1时刻的解为:

      $$u_{n + 1}^t = {\rm{Re}} ({u_{n{\rm{ + 1}}}}) + u_{n + 1}^{v1} + u_0^{v4}$$ (42)

      为此,中心差分法具体的计算步骤如表1所示。

      表 1  任意荷载下中心差分虚初始条件法计算步骤

      Table 1.  Calculation procedure of central differential virtual initial condition method under arbitrary load

      1. 准备工作
      1.1 对加速度时程a(t)根据式(28)和式(29)计算系数
      1.2 计算虚初始条件,形成总初始条件$u_0^t = u_0^{} + u_0^{v1} + u_0^{v{\rm{4}}}$, $v_0^t = v_0^{} + v_0^{v1} + v_0^{v{\rm{4}}}$
      1.3 $\ddot u_0^t = \left[ { {f_0} - k(1 + {\rm{i}}\eta )u_0^t} \right]/m$,$u_{ - 1}^t = u_0^t - \Delta t\dot u_0^t - {(\Delta t)^2}\ddot u_0^t{\rm{/2}}$
      1.4 选择Δt
      1.5 $\hat k = m{\rm{/}}\Delta {t^2}$
      2. 每个荷载步计算
      2.1 ${\hat f_{n + 1} } = {f_n} - \left[ {k(1 + {\rm{i}}\eta ) - \dfrac{ {2m} }{ {\Delta {t^2} } } } \right]u_n^t - \dfrac{m}{ {\Delta {t^2} } }u_{n - 1}^t$
      2.2 ${u_{n + 1}} = {\hat f_{n + 1}}/\hat k$
      2.3 ${\ddot u_{n + 1} } = \dfrac{1}{m}[ { {f_{n + 1} } - k(1 + {\rm{i}}\eta ){u_{n + 1} } } ]$,${\dot u_{n{\rm{ + 1} } } } = \dfrac{ {\Delta t} }{2}{\ddot u_{n + {\rm{1} } } }{\rm{ + } }\dfrac{ { {u_{n + {\rm{1} } } } - {u_n} } }{ {\Delta t} }$
      2.4 $u_{n{\rm{ + } }1}^{v1} = - {\rm{i}}\dfrac{ { {\rm{Re} } ({ {\dot u}_{n + 1} }) + \omega \alpha {\rm{Re} } ({u_{n + 1} })} }{ {\omega \mu } }$
      2.5 $u_{n + 1}^t = {\rm{Re}} ({u_{n{\rm{ + 1}}}}) + u_{n + 1}^{v1} + u_0^{v{\rm{4}}}$
      3. 用n+1代替n,重复2.1~2.5进行下一荷载步的计算。
    • 为验证本文算法的有效性,下面以文献[17]分析过的算例进行本文算法的验证。已知体系的滞后阻尼系数η=0.1,对无阻尼自振频率f=0.1 Hz、1 Hz、10 Hz的三个体系进行地震反应计算。以表2中的3条地震波分别作为体系的地震输入。输入地震波的加速度时程如图1所示。在进行中心差分法计算时,根据中心差分法的稳定性条件和输入地震波的采样时间间隔,计算时间步长$\Delta t$取为:

      $$\Delta t \leqslant \min ({T_n}/20,\Delta {t_{\rm{s}}})$$ (43)

      式中:Tn为体系的自振周期;Δts为地震波加速度时程的采样时间间隔。

      表 2  地震波主要参数

      Table 2.  The main parameters of ground motions

      地震波时间地震峰值加速度/(m/s2)
      El Centro1940.5.18California3.417
      迁安波1976.8.31唐山余震0.974
      天津波1976.11.25唐山余震1.458

      图  1  输入地震波的加速度时程

      Figure 1.  Acceleration time histories of input seismic waves

      作为对比,对地震作用形成式(31)的方程后,计算各项简谐荷载的解析解并求和,所得结果作为解析解。当$u(0) = {u_0}$, $\dot u(0) = {v_0}$时,位移、速度和加速度的解析解为:

      $$\tag{44a}u = C{{\rm{e}}^{ - st}} + \sum\limits_{j = 0}^{N/2} {{C_{{\rm{f}}j}}{{\rm{e}}^{ - st}}} + \sum\limits_{j = 0}^{N/2} {{X_j}{{\rm{e}}^{{\rm{i}}{\theta _j}t}}} \qquad\quad\quad$$
      $$\tag{44b}\dot u = - sC{{\rm{e}}^{ - st}} + \sum\limits_{j = 0}^{N/2} {( - s{C_{{\rm{f}}j}}{{\rm{e}}^{ - st}})} + \sum\limits_{j = 0}^{N/2} {({\rm{i}}{\theta _j}{X_j}{{\rm{e}}^{{\rm{i}}{\theta _j}t}})} $$
      $$\tag{44c}\ddot u = {s^2}C{{\rm{e}}^{ - st}} + \sum\limits_{j = 0}^{N/2} {({s^2}{C_{{\rm{f}}j}}{{\rm{e}}^{ - st}})} + \sum\limits_{j = 0}^{N/2} {( - \theta _j^2{X_j}{{\rm{e}}^{{\rm{i}}{\theta _j}t}})} $$

      对式(31)两边进行Fourier变换[27],可得:

      $$ - \theta _j^2mU({\theta _j}) + (1 + {\rm{i}}\eta )kU({\theta _j}) = F({\theta _j})$$

      式中,$F({\theta _0}) = m{A_0}/2$, $F({\theta _j}) = m({A_j} - i{B_j})$(j=1, 2, ···, N/2-1),$F({\theta _{N/2}}) = m{A_{N/2}}/2$,则:

      $$U({\theta _j}) = {X_j}$$ (46)

      利用式(46)进行Fourier逆变换,可得位移、速度和加速度的频域解为:

      $$\tag{47a}u = \sum\limits_{j = 0}^{N/2} {{X_j}{{\rm{e}}^{{\rm{i}}{\theta _j}t}}}\quad\; $$
      $$\tag{47b}\dot u = \sum\limits_{j = 0}^{N/2} {{\rm{i}}{\theta _j}{X_j}{{\rm{e}}^{{\rm{i}}{\theta _j}t}}} \;$$
      $$\tag{47c}\ddot u = - \sum\limits_{j = 0}^{N/2} {\theta _j^2{X_j}{{\rm{e}}^{{\rm{i}}{\theta _j}t}}} $$

      对比式(44)和式(47)可知,解析解包括由初始条件引起的自由振动,由荷载引起的伴生自由振动和荷载所引起的纯强迫振动。而频域解仅包含荷载所引起的纯强迫振动。

      u0=v0=0,f=10 Hz时,三条地震波作用下体系的位移反应如图2所示,f=10 Hz时天津波作用下的速度和加速度反应如图3所示,f=1 Hz和0.1 Hz时天津波作用下的位移反应如图4所示。

      图  2  三条地震作用下位移反应(f=10 Hz, u0=v0=0)

      Figure 2.  Displacement responses under three seismic excitations (f=10 Hz, u0=v0=0)

      图  3  天津波作用下的速度和加速度反应(f=10 Hz, u0=v0=0)

      Figure 3.  Velocity and acceleration response under Tianjin wave excitation

      图  4  天津波作用下不同自振频率体系的位移反应(u0=v0=0)

      Figure 4.  Displacement response for various natural vibration frequencies under Tianjin wave excitation

      图2~图4可知,对于不同的地震输入和自振频率体系,中心差分虚初始条件法计算的位移、速度和加速度都和解析解基本一致,计算结果稳定而不存在临界自振周期的问题[17]。而对于频域解,当f=10 Hz时与解析解几乎重合;当f=0.1 Hz时,频域解与解析解存在明显差别。这是因为频域解得到的是体系的稳态反应,而解析解是包含稳态和伴生自由振动的瞬态反应,自振频率越高,瞬态反应衰减越快,因此,对于自振频率较高的体系,瞬态反应对体系的总反应影响很小而可直接采用稳态解进行描述,但是,对于自振频率低的体系,瞬态反应衰减需要较长的时间,当自振周期和地震波持时为同一量级时,则瞬态反应在整个地震反应过程中都有显著影响,此时采用仅包含稳态解的频域分析方法将产生显著误差。本文所提的中心差分虚初始条件法,计算结果包含瞬态反应和稳态反应,对于不同自振周期体系的反应都和解析解吻合的很好。

      当体系初始条件非零时,在天津波作用下不同自振频率的位移反应如图5所示,f=10 Hz体系的速度与加速度反应如图6所示。初始条件非零情况与零初始条件下的计算结果对比可知,由于初始位移和初始速度瞬态振动的进一步影响,导致频域解在初始阶段误差增大。对于自振频率较高的10 Hz体系,由于瞬态振动很快消失,此时忽略初始条件对体系总反应的影响较小。但是对于自振周期与输入地震波持时为同一量级的0.1 Hz体系,初始位移和初始速度引起的瞬态振动将对总反应产生显著影响,此时忽略初始条件将引起显著误差。

      图  5  天津波作用下不同自振频率体系的位移反应

      Figure 5.  Displacement response of various natural vibration frequencies under Tianjin wave excitation

      图  6  天津波作用下体系的速度和加速度反应(f=10 Hz)

      Figure 6.  Velocity and acceleration response under Tianjin wave excitation

      为定量研究不同方法的精确性,不同方法相对解析解的峰值相对误差为:

      $${e_{\rm{r}}} = \frac{{\left| {{{\left| {r(t)} \right|}_{\max }}{\rm{ - }}{{\left| {{r^*}(t)} \right|}_{\max }}} \right|}}{{{{\left| {{r^*}(t)} \right|}_{\max }}}} \times 100\% $$

      式中:${\left| {r_{}^*(t)} \right|_{\max }}$表示精确解的峰值;${\left| {r(t)} \right|_{\max }}$表示近似解的峰值。

      零初始条件下本文方法及频域解的峰值相对误差如表3所示。由表3可知,当体系的自振频率较高时,频域解的计算精度高;而当体系的自振频率较低时,频域解将产生显著误差;而中心差分虚初始条件法对于不同自振周期体系反应的峰值相对误差都小于5%,显示了良好的精度。

      表 3  本文方法与频域解的峰值相对误差

      Table 3.  Peak relative errors of the proposed method and frequency domain solution

      地震波自振频率/Hz${e_u}$/(%)${e_{\dot u}}$/(%)${e_{\ddot u}}$(%)
      本文方法频域解本文方法频域解本文方法频域解
      El Centro波101.460.00020.2101.370
      10.640.00290.540.020.080.0218
      0.10.6917.130.3637.960.391.74
      迁安波100.0800.00400.270
      11.230.01041.540.0060.110.0033
      0.13.3618.080.932.570.120.10
      天津波100.760.00021.9302.210
      10.190.420.770.420.060.25
      0.10.536.080.320.170.180.07

      表4u0=v0=0,f=1 Hz时,天津波作用下三种方法位移反应的计算时间。解析解为输入地震波的Fourier变换时间和式(44a)各项直接求和的计算时间,中心差分虚初始条件法的计算时间包括输入地震波的Fourier变换时间、计算虚初始条件的时间和每个时刻位移、速度和加速度的递推计算的时间,频域法则为Fourier变换和Fourier逆变换的总时间。在进行解析解计算时,由于自由振动和伴生自由振动的影响,无法采用快速Fourier变换进行计算而采用逐项求和,因此,计算时间最长。频域采用快速Fourier变换,计算时间最短。中心差分虚初始条件法通过前一步的计算结果递推后一步的计算结果,计算效率仅次于频域方法,且计算时间显著小于解析解。结合表3表4的数据分析结果表明,中心差分虚初始条件法兼顾了计算精度和计算效率。

      表 4  不同算法计算时间比较

      Table 4.  Comparison of calculation time of various methods /s

      地震波解析解本文方法频域解
      天津波1.400.0030.0015
    • 针对滞后阻尼体系直接积分法收敛到包含发散项而不稳定的问题,提出了中心差分虚初始条件计算方法。由理论分析和数值计算结果可知:

      (1)对于不同峰值加速度、卓越频率的地震输入,及不同自振频率体系,本文所提的中心差分虚初始条件法对位移、速度和加速度均能得到稳定的计算结果,不存在临界自振周期的问题。实初始条件是可观测和测量部分,虚初始条件是伴随实初始条件而存在的,可无需人为干涉而使滞后阻尼体系计算结果稳定。由此进行直接积分法计算,即使有条件稳定的中心差分法,依然可以得到稳定的计算结果。

      (2)对于自振频率较高的体系,瞬态反应对体系的总反应影响很小而可直接采用稳态解进行描述,而对于自振频率低的体系,如自振周期和地震波持时为同一量级时,瞬态反应在整个地震反应过程中都有显著影响,此时采用仅包含稳态解的频域分析方法将产生显著误差。

      (3)中心差分虚初始条件法的计算结果包括瞬态反应和稳态反应,计算误差和体系的自振频率无关,可适用于不同的自振频率体系。

      (4)中心差分虚初始条件法计算的峰值相对误差小于5%。同时,计算时间显著小于解析解,因此,这种方法兼顾了计算精度和计算效率。

WeChat 关注分享

返回顶部

目录

    /

    返回文章
    返回