The Application of Differential Quadrature Method in Structural Dynamic Analysis
-
摘要: 微分求积法(DQM)是1种求解微分方程初(边)值问题的数值方法,通常以较小的计算工作量即可获得较高的数值精度。这种方法应用于工程领域时多用来解决梁、板等结构的静力分析或结构特征值分析等问题,即对边值问题的微分方程的求解。结构动力分析属于初值问题,荷载和结构反应都具有特殊性,直接套用DQM求解边值问题并不能获得问题的解。本文尝试利用微分求积原理建立求解结构动力反应的具体方法。借鉴单元法的思想,将荷载持时划分为若干个时步,在每个时步内对动态荷载和结构反应进行离散,然后用DQM对时步逐个进行求解,得到体系在整个时域内的反应过程。通过对3种不同自振周期的线弹性单自由度体系在不同频率简谐激励下反应的计算,阐释了本文方法的可行性以及高精度、高效率的特点,通过数值试验确定了时步内相对较优的节点数,并为时步长度的选取提供了建议。Abstract: The differential quadrature method (DQM) is a numerical technique of solving initial/boundary value problems of differential equations and capable of obtaining a higher numerical accuracy with a smaller calculation workload. This method is often used to solve the problems of structural static analysis of beams and slabs or eigenvalue analysis when it is applied to the engineering fields, which is to solve the differential equation of the boundary value problems. Dynamic analysis of structures is an initial value problem as well as particular loads and structural response. As a result, applying the DQ method of solving the boundary value problem directly cannot obtain solution of problem. The principle of differential quadrature is applied to establish the specific method of solving structural dynamic response in this paper. By analogizing the idea of unit method, the duration of the load is divided into many time steps and dynamic load and structure are discretized in each time step, then the response of the system can be solved in the whole time domain by employing the DQ method step by step. The feasibility of this method and the characteristics of high precision and high efficiency are expounded by calculating the response of three linear elastic single-freedom-degree system of different natural vibration periods excited by simple harmonic loads of different frequencies. By means of numerical experiment, the optimal meshing scheme is determined and the suggestion for the time step is given.
-
引言
唐山地震台内建有跨断裂水准测线和基线测线,主要用于监测唐山断裂的活动状况,40余年来积累了丰富的资料,为研究唐山断裂的活动奠定了基础,近年来唐山地震台又陆续增加了断层CO2测量和地应力测量等手段用于监测唐山断裂活动。其中跨断层形变测量是积累时间最长、观测精度较高、能够直接反映断层活动的形变观测手段。许多学者利用跨断层形变资料进行了研究,并取得了较丰硕的研究成果(李文静等,2009;黄建平等, 2010, 2011;周海涛等,2013)。利用跨断层形变数据可求解断层活动协调比参数,而断层活动协调比又对地震预测研究有一定指示意义(张晶等,2011)。
断层土壤气是地球内部沿着活动板块或块体边界及其它地壳薄弱带向地表迁移释放的气体。监测跨断层土壤气释放浓度的变化,已成为探索地震前兆与地震预测及评价断裂活动性的重要方法(汪成民等,1991;刘菁华等,2006)。断层土壤气CO2是地球内部生成的众多流体组分中最有可能大量迁移至地表并在地表某点集中释放的气体之一,它的异常浓度和通量,可以很好地反映地震活动和断裂带的活动情况(周晓成等,2012;张扬等,2016)。
唐山地震台从2011年7月开始进行断层CO2测量,每日监测1次,台站场地内的断层、基线和水准测线、CO2观测点位如图 1所示。本文通过对唐山地震台跨断层形变和CO2资料进行分析处理,研究唐山断裂的运动特征、断裂活动协调比及CO2变化特征,以期总结出相关规律,为地震监测预测提供有益的借鉴。
1. 断裂形变特征
基线观测可以监测到断层两盘间的水平形变,本文忽略应变部分,将基线变化近似地看作由断层两盘相对水平运动所引起,研究断层水平运动特征。一般的跨断层基线与断层的位置关系见图 2。唐山地震台共布设4条基线测线和4条水准测线,其中2条水准测线跨断层,4条基线测线均跨断层。水平张压、走滑分量采用公式(1)计算,通过4条跨断层基线数据并利用最小二乘法求解断裂的水平向运动量值(周海涛等,2009)。断层垂向运动分量通过2条跨断层水准变化量取平均值得出。
$$ \left. \begin{gathered} \delta {L_1} = x\sin {\alpha _1} + y\cos {\alpha _1} \hfill \\ \delta {L_2} = x\sin {\alpha _2} + y\cos {\alpha _2} \hfill \\ \end{gathered} \right\} $$ (1) 式中:x为张压分量(垂直于断裂走向,张性运动为正);y为走滑分量(沿断裂走向,左旋走滑为正);δL1和δL2分别是基线1、基线2的2期观测值之间变化量(伸长为正);α1、α2分别是基线1、基线2与断裂走向的夹角(由断层走向顺时针转动到基线方向所转过的角度)。
跨断层测量资料为1984年1月1日—2018年3月6日的数据。将计算结果绘制成时序曲线,见图 3。
从图 3(a)可以看出,断层在2003年以前以拉张为主(曲线上升为拉张,下降为压缩,下同),2003—2008年主要表现为压缩,2008—2011年主要表现为拉张,2011年主要表现为压缩,2012年至今张压运动不明显。断层在各时段张压运动互相交替,但总体上来说,断层呈弱张性运动。
从图 3(b)可以看出,断层在1997年以前主要以右旋走滑为主(曲线下降为右旋,上升为左旋,下同);1997—2002年表现为左旋走滑;2003—2012年呈右旋走滑趋势,2012年至今走滑运动不明显。断层在各时段左旋和右旋走滑运动互相交替,但总体上来说,断层呈弱右旋走滑运动。
从图 3(c)可以看出,断层在1997年以前主要表现为正断活动(曲线上升为正断;曲线下降为逆断;下同)。1997年至今正断和逆断活动变化不明显。
唐山断裂运动变化时序曲线表明,唐山断裂近35年来整体上在水平方向呈微弱的右旋张性活动,垂直方向呈正断活动,但不同时段的运动状态有所不同。
2. 断裂活动协调比
张晶等(2011)提出了断层活动协调比参数的概念,认为走滑兼倾滑无障碍蠕动是1种无应变积累的活动状态,其特征可用刚体运动模型来描述。当断层在动平衡系统下处于无障碍自由蠕滑时,可认为是1种相对稳定状态。当作用力的大小和方向在单位时间内不变时,断层的运动增量分别为走滑分量、拉张分量和垂直分量,定义这3个分量中2个分量之比为断层活动协调比参数。
断层活动协调比作为走滑兼倾滑无障碍断层蠕动模型的特征参数,如果接近常数,就可近似认为断层活动基本为无障碍蠕动。当断层活动协调比是变量,偏离正常值很大时,在排除了非构造活动(干扰及人为影响等)的情况下,可以认为断层的常态活动发生了变化,断层蠕动趋向不稳定,可能产生新的断层异常活动,也可能形成新的断层闭锁,造成协调比趋势的偏离,这表明断层或附近一定有新的应力增强或应变积累。当断层失稳发生强震,断层的高应变能释放之后,断层活动又趋于刚性错动模式,协调比又恢复稳定形态(张晶等,2011)。
通过跨断层观测的基线和水准经推导可得到断层活动协调比参数。首先,由基线、水准同测数据可求得断层水平走滑分量(y)、水平张压分量(x)和垂直运动分量(h)。
将断层两盘在台站测量范围内的相对运动近似看成2块刚性块体间的相对运动,并分解成水平走滑分量、水平张压分量和垂直运动分量的相对运动的3个分量。这3个分量即前文得出的断层水平走滑分量(y)、水平张压分量(x)和垂向运动分量(h),将2个分量之比作为断层活动协调比。为了便于分析,把断层活动协调比中超过3倍标准差的数值去掉,然后将所得出的结果绘制协调比散点图,如图 4所示。
图 4(a)中的协调比散点图为断层水平张压分量与水平走滑分量比值的时间序列,从图中可以看出在1984—1986年该协调比处于比较离散的状态,1976年唐山7.8级地震后,断层活动处于调整阶段,局部存在一定的应变积累或释放,因此会与断层的常态活动存在差异;1985—2000年的协调比相对比较稳定;2001—2004年处于相对离散状态,研究发现2004年1月20日在距唐山断裂约40km的滦县-乐亭断裂上发生了4.6级地震,认为可能与这次震前地壳应力场的异常活动有关,引起了断层协调比的异常变化。
图 4(b)中的协调比散点图为断层垂直运动分量与水平走滑分量比值的时间序列,从图中可以看出,此协调比曲线与图 4(a)中的协调比曲线状态几乎一致,可得到与图 4(a)完全一致的解释,同时其幅值更大,说明在本场地内断层垂直运动信息比水平运动信息更显著。
图 4(c)中的协调比散点图为断层垂直运动分量与水平张压分量比值的时间序列,从图中可以看出在1984—1988年该协调比处于比较离散的状态,应是1976年唐山7.8级地震后断层活动调整的反映。1989—2008年的协调比相对比较稳定,但在2008年汶川8.0级地震之后不久,协调比出现1个异常上升现象,而该异常变化的出现主要是由于水平张压量值的减小。该变化是否与汶川8.0级地震后华北地区应力调整有关还有待进一步的研究。
图 4中显示的散点图大部分近似直线,其实是由于比例尺的原因。若放大比例尺,可以看出散点所构成的曲线在中线附近上下波动,且波动幅度较小,在相对较小比例尺下难以分辨。
通过分析唐山地震台所处位置的唐山断裂协调比可以看出,唐山7.8级地震引起的华北地壳应力场调整并没有完全结束,事实上唐山7.8级地震后也陆续发生了一系列的较大余震。而2001—2004年的应力调整之前发生了1998年张北6.2级地震、2004年滦县4.6级地震和2006年文安5.1级地震。与这些地震相应的地壳应力场的异常活动在唐山台断层形变观测中都有所体现。
3. 断层土壤气CO2特征
断裂活动使围岩发生破裂而释放出大量的CO2气体,CO2沿断裂上升并通过附近的土壤层释放出来。因此,在有利于气体释放的部位,如活动断裂带,观测CO2的动态特征,研究CO2的释放规律和异常活动,有可能获得与地壳运动相关的信息,进而分析断裂活动的情况,捕捉到地震前兆异常,为地震预报提供依据。
为了监测唐山断裂土壤气CO2动态变化特征,在唐山地震台设有2处断层土壤气CO2观测点,分别位于唐山断裂的上盘和下盘,采用CO2快速测定管进行测量。CO2快速测定管是长约14cm、直径约0.45cm的玻璃管。管内充填了吸附有百里蓝酚酞的活性氧化铝,测定管本身为蓝色,当CO2进入测定管,即被活性氧化铝吸附,与管内的百里蓝酚酞发生反应,使测定管的颜色变为白色。CO2根据测定管变白色柱的长度,从刻度上即可直接读出CO2的含量数值。
由于断层土壤气CO2与地温有一定的关系,因此在监测断层土壤气CO2同时还在同一地点、同一深度加测地温。另外,降雨量一般也会对断层土壤气CO2产生一定影响,因此本研究中收集了唐山地震台降雨量资料进行综合分析。
对唐山地震台所在位置的唐山断裂土壤气CO2含量(断层下盘CO2含量)、测定孔内温度及台站降雨量进行绘图分析,如图 5所示。由图可见地温变化呈明显的年周期变化,断层土壤气CO2浓度与地温呈明显的正相关关系。而降水量较大时,断层土壤气CO2与降雨量呈一定的负相关关系。
由图 5可以看出,断层气CO2受地温和降水量的影响较大。一般认为降雨的影响是因为CO2在水中溶解度较大,当雨季测孔进水,孔下沿断层上升至测孔的CO2大量溶解于水中,可导致断层气CO2测值比正常值低数倍的变化。对于观测中的降雨干扰,一般根据测点附近的降雨资料定性识别和排除。同时,断层气CO2与地温呈显著的正相关关系。断层气CO2释放量具有明显的夏高冬低年变特点,一般认为这种现象与地下的生物化学作用有关(Sugisaki等,1983;王基华等,1994),在提取断层气CO2前兆异常时,一般可进行地温改正,用以消除地温干扰。
对断层气CO2采用傅立叶滑动去年周期的方法去除年变周期,结果如图 6所示。可以看出,图 6(a)中已没有明显的年变周期现象,2014年之后断层气CO2浓度的离散性比2014年之前大很多;与降水量进行对比研究可以看出两者的负相关性很明显。
为了研究唐山断裂土壤CO2与地震的对应关系,选取2011年7月以来在唐山断裂附近发生的2个4级以上地震进行研究,分别为2012年5月28日发生的4.8级地震和2016年9月10日发生的4.2级地震。研究发现,在2012年5月28日4.8级地震发生之前,3月份断层土壤气CO2浓度处于低值,然后逐渐上升,地震发生之前的5月23日,断层土壤气CO2浓度大幅下降,之后恢复上升状态,并于5日后发生地震,之后保持了一定时期的高值。2016年9月10日4.2级地震发生之前的7月下旬至8月下旬,断层土壤气CO2浓度大幅下降,并剧烈波动,这可能与降水量变大有一定关系,之后在8月底至9月初大幅上升,在地震发生之前小幅波动,后发生地震,地震之后保持了一定时期的高值。由此可见,唐山地震台断层土壤气CO2对唐山断裂附近发生的4级以上地震有一定的显示。
4. 结论
通过对唐山地震台内的跨断层基线和水准、断层土壤气CO2及台站降水量等资料进行分析研究。粗略地给出了唐山地震台内的唐山断裂形变特征、断层活动协调比、CO2变化特征。
(1)通过唐山断裂运动变化时序曲线可以看出,唐山断裂近35年来整体上在水平方向呈微弱的右旋张性活动,垂直方向呈正断活动,但不同时段的运动状态有所不同。
(2)通过分析唐山地震台所处位置的唐山断裂协调比认为,唐山7.8级地震引起的华北地壳应力场调整并没有完全结束;
(3)近年来唐山地震台内的唐山断裂断层土壤气CO2变化与地温、降水量有很大关系,并在台站周边发生地震前有一定的显示。
(4)综合分析认为,唐山地震台形变和流体前兆观测对附近及华北因区域地壳应力场调整引起的中小地震有一定的显示,说明该台具有一定的前兆异常显示能力;随着观测资料和经验的持续积累,有希望在该区域发生更大地震前捕捉到更显著的前兆异常。
致谢: 感谢薄万举研究员对本文的指导和帮助。 -
表 1 体系的基本特性
Table 1. Basic characteristics of systems
体系编号 自振周期Tn/s 阻尼比ξ 1 0.68 0.05 2 0.25 0.05 3 0.08 0.05 表 2 简谐荷载的信息
Table 2. Information of simple harmonic load
荷载函数/N 周期/s 振幅/m 持时/s sin2πt 1 1 40 sin10πt 0.2 1 40 sin20πt 0.1 1 40 表 3 荷载周期1s的平均相对误差
Table 3. Average relative error with load period of 1s
m 体系编号 1 2 3 位移/% 速度/% 位移/% 速度/% 位移/% 速度/% 2 98.734 98.734 98.734 98.734 98.734 98.734 4 88.218 81.849 99.415 389.850 53.532 1293.500 6 42.402 32.547 15.926 4.728 10.979 9.609 8 6.698 8.020 1.543 7.814 2.042 2.741 10 1.623 2.687 0.907 2.338 1.036 0.349 12 0.121 0.138 0.897 3.341 0.819 4.156 14 0.008 0.009 0.935 5.153 0.716 1.290 16 0 0 0.930 3.899 0.635 0.825 18 0 0 1.005 3.292 0.577 2.28 20 0 0 1.268 4.715 0.527 0.899 表 4 荷载周期0.2s的平均相对误差
Table 4. Average relative error with load period of 0.2s
m 体系编号 1 2 3 位移/% 速度/% 位移/% 速度/% 位移/% 速度/% 2 99.749 99.749 99.749 99.749 99.749 99.749 4 799.760 2460.400 403.870 411.780 108.540 225.480 6 90.884 11.229 23.620 19.707 11.436 5.439 8 6.052 15.762 2.102 2.135 1.068 2.520 10 0.264 0.046 0.116 0.140 0.765 1.155 12 0.009 0.022 0.005 0.004 0.807 3.325 14 3681.600 256.240 251.760 2935.800 972.740 1334.400 16 0 0 0 0 0.093 0.167 18 0 0 0 0 0.010 0.013 20 0 0 0 0 0.001 0.001 表 5 荷载周期为0.1s的平均相对误差
Table 5. Average relative error with load period of 0.1 seconds
m 体系编号 1 2 3 位移/% 速度/% 位移/% 速度/% 位移/% 速度/% 2 99.875 99.875 99.875 99.875 99.875 99.875 4 1620.200 1370.200 578.690 459.540 51.868 52.298 6 175.570 9.444 65.345 7.903 6.80×10282 8.63×10282 8 11.571 9.569 4.413 3.254 1.951 2.655 10 0.498 0.046 0.204 0.044 0.142 0.266 12 0.016 0.013 0.007 0.005 0.006 0.007 14 27683 306.960 1322.400 246.430 234.600 3888.600 16 0 0 0 0 0 0 18 0 0 0 0 0 0 20 0 0 0 0 0 0 -
李鸿晶, 王通, 2011a.结构地震反应分析的逐步微分积分方法.力学学报, 43(2):430-435. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=CAS201303040000629782 李鸿晶, 廖旭, 王通, 2011b.结构地震反应DQ解法的两种数值格式.应用基础与工程科学学报, 19(5):758-766. http://d.old.wanfangdata.com.cn/Periodical/yyjcygckxxb201105008 廖旭, 李鸿晶, 孙广俊, 2013.基于DQ原理的结构弹塑性地震反应分析.工程力学, 30(7):161-166. http://www.cnki.com.cn/Article/CJFDTotal-GCLX201307026.htm 王鑫伟, 1995.微分求积法在结构力学中的应用.力学进展, 25(2):232-240. doi: 10.3321/j.issn:1000-0992.1995.02.010 Bellman R., Casti J., 1971. Differential quadrature and long-term integration. Journal of Mathematical Analysis and Applications, 34 (2):235-238. doi: 10.1016/0022-247X(71)90110-7 Bellman R., Kashef B. G., Casti J., 1972. Differential quadrature:a technique for the rapid solution of nonlinear partial differential equations. Journal of Computational Physics, 10 (1):40-52. doi: 10.1016/0021-9991(72)90089-7 Bert C. W., Jang S. K., Striz A. G., 1988. Two new approximate methods for analyzing free vibration of structural components. AIAA Journal, 26 (5):612-618. doi: 10.2514/3.9941 Bert C. W., Wang X. W., Striz A. G., 1993. Differential quadrature for static and free vibration analyses of anisotropic plates. International Journal of Solids and Structures, 30 (13):1737-1744. doi: 10.1016/0020-7683(93)90230-5 Bert C. W., Malik M., 1996. The differential quadrature method for irregular domains and application to plate vibration. International Journal of Mechanical Sciences, 38 (6):589-606. doi: 10.1016/S0020-7403(96)80003-8 Bert C. W., Malik M., 1997. Differential quadrature:a powerful new technique for analysis of composite structures. Composite Structures, 39 (3-4):179-189. doi: 10.1016/S0263-8223(97)00112-8 Fung T. C., 2001a. Solving initial value problems by differential quadrature method——part 1:first-order equations. International Journal for Numerical Methods in Engineering, 50 (6):1411-1427. doi: 10.1002/(ISSN)1097-0207 Fung T. C., 2001b. Solving initial value problems by differential quadrature method——part 2:second-and higher-order equations. International Journal for Numerical Methods in Engineering, 50 (6):1429-1454. doi: 10.1002/(ISSN)1097-0207 Jang S. K., Bert C. W., Striz A. G., 1989. Application of differential quadrature to static analysis of structural components. International Journal for Numerical Methods in Engineering, 28 (3):561-577. doi: 10.1002/(ISSN)1097-0207 Kang K., Bert C. W., Striz A. G., 1995. Vibration analysis of shear deformable circular arches by the differential quadrature method. Journal of Sound and Vibration, 183 (2):353-360. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=cd5d834cd6a527d56839a5876ae0d87a Kang K., Bert C. W., Striz A. G., 1996. Vibration analysis of horizontally curved beams with warping using DQM. Journal of Structural Engineering, 122 (6):657-662. doi: 10.1061/(ASCE)0733-9445(1996)122:6(657) Liew K. M., Han J. B., Xiao Z. M., 1996a. Differential quadrature method for thick symmetric cross-ply laminates with first-order shear flexibility. International Journal of Solids and Structures, 33 (18):2647-2658. doi: 10.1016/0020-7683(95)00174-3 Liew K. M., Han J. B., Xiao Z. M., et al., 1996b. Differential quadrature method for Mindlin plates on Winkler foundations. International Journal of Mechanical Sciences, 38 (4):405-421. doi: 10.1016/0020-7403(95)00062-3 Liu J., Wang X. W., 2008. An assessment of the differential quadrature time integration scheme for nonlinear dynamic equations. Journal of Sound and Vibration, 314 (1-2):246-253. doi: 10.1016/j.jsv.2008.01.004 Newmark N. M., 1959. A method of computation for structural dynamics. Journal of the Engineering Mechanics Division, 85 (3):67-94. http://cn.bing.com/academic/profile?id=04737dd237295d8ef1e115cb61eed609&encoded=0&v=paper_preview&mkt=zh-cn Sherbourne A. N., Pandey M. D., 1991. Differential quadrature method in the buckling analysis of beams and composite plates. Computers & Structures, 40 (4):903-913. doi: 10.1016-0045-7949(91)90320-L/ Striz A. G., Jang S. K., Bert C. W., 1988. Nonlinear bending analysis of thin circular plates by differential quadrature. Thin-Walled Structures, 6 (1):51-62. doi: 10.1016/0263-8231(88)90025-0 Wang X., Bert C. W., Striz A. G., 1993. Differential quadrature analysis of deflection, buckling, and free vibration of beams and rectangular plates. Computers & Structures, 48 (3):473-479. http://cn.bing.com/academic/profile?id=5a95130e65486d67576bb93101054491&encoded=0&v=paper_preview&mkt=zh-cn Wang X., Striz A. G., Bert C. W., 1994. Buckling and vibration analysis of skew plates by the differential quadrature method. AIAA Journal, 32 (4):886-889. doi: 10.2514/3.12071 Zeng H., Bert C. W., 2001. A differential quadrature analysis of vibration for rectangular stiffened plates. Journal of Sound and Vibration, 241 (2):247-252. doi: 10.1006/jsvi.2000.3295 -