一种大跨度桥梁地震响应时程分析方法与流程

本发明涉及一种大跨度桥梁地震响应时程分析方法，具体涉及一种计算大跨度桥梁在地震作用下响应的无条件稳定、高阶精度、具有最佳数值阻尼的动力时程分析方法。

背景技术：

桥梁是交通运输系统的枢纽工程,是生命线工程的重要组成部分,保证桥梁在地震中的安全和正常使用,对城市和地区的抗震防灾工作和震后恢复重建工作中都具有重要的意义。在设计阶段，常常采用数值模拟手段，对桥梁进行地震作用下的分析，以保证桥梁结构在突发强震作用下的安全。为便于计算，人们提出了一些例如底部剪力法、反应谱法等方法。然而随着土木工程结构形式的不断发展和对动力分析可靠度越来越高的要求，我们需要对结构在动力荷载下的行为进行更为精确的预测。

动力时程分析法是随着强震记录的增多和计算机技术的广泛应用而发展起来的,是公认的精细分析方法。国内外大多数工程抗震设计规范中都指出,对于大跨度桥梁的地震反应分析,要采用动态时程分析法。也就是说,要采用多节点自由度的结构有限元模型,把地震强迫振动的激励一地震加速度时程直接输入,对结构进行地震时程反应分析。动力时程分析法可以精确的考虑地基和结构的相互作用,地震波相位差及不同地震波多分量多点输入,结构的各种复杂非线性因素以及分块阻尼等问题。

目前现有技术针对地震作用下大跨桥梁结构进行动力时程分析主要采用的是newmark法，wilson法和hht法等，但目前这些方法存在或多或少的问题。对于大跨度桥梁结构，由于结构的最长周期和最短周期相差悬殊，因此该问题属于刚性问题，对于刚性问题的动力时程分析，尤其是对于加速度响应的计算，容易产生寄生振荡，这种震荡是由于有限元离散造成的，对结构的实际响应的影响是负面的，需要通过数值阻尼过滤掉。常用的时程分析方法由于没有足够大的数值阻尼，无法过滤掉这种虚假震荡，甚至造成方法不稳定，导致计算不收敛。

为解决此问题，bathe发展了一种有效的方法，对于刚性问题的求解非常有效，目前已经嵌入到adina软件当中。然而bathe方法的计算效率较低，每个时间步上需要进行两次矩阵分解，且方法仅仅具有2阶精度。在使用这些具有2阶精度的分析方法进行动力时程分析时，为保证精度，必须采用较短的时间步长。因为大跨桥梁模型自由度数量庞大，时间步长选的过短，地震输入时长不变，总的时间步数量过多，也就造成了分析过程时间耗费巨大。

技术实现要素：

针对传统动力时程分析方法在求解大跨度桥梁动力时程分析方面的缺陷，本发明的目的在于提供一种计算大跨度桥梁在地震作用下响应的无条件稳定、高阶精度、具有最佳数值阻尼的动力时程分析方法。

本发明的技术方案如下：

一种大跨度桥梁地震响应时程分析方法，其步骤如下：

第一步，对大跨度桥梁结构进行空间有限元离散，建立该结构的有限元模型离散系统，主梁、斜塔、斜拉索、主缆等构件均采用空间伯努利欧拉梁单元，并采用rayleigh阻尼建立单元阻尼矩阵，由单元刚度矩阵、单元质量矩阵和单元阻尼矩阵集成整体刚度矩阵、整体质量矩阵和整体阻尼矩阵，并由hamilton原理导出离散系统的运动方程组：

第二步，选取时间步长，时间步长取为加速度记录时间间隔的n倍，其中n为整数，由于方法具有高阶精度，n可选取为较大的值；

第三步，逐时间步计算，对于第i个时间步，已知ti-1时刻的位移ui-1和速度vi-1，由以下几个子步骤，计算ti时刻的位移ui、速度vi；

a.由方程d11ui＝r1求解位移ui，其中，

r1＝p1+g11ui-1+δtg12vi-1

b.由方程d11vi＝r2求解位移vi，其中，

r2＝(p2+g21ui-1+δtg22vi-1-d21ui)/δt

c.由方程d11ei＝δr1求解位移的误差ei，其中，

d.由方程d11εi＝δr2求解速度的误差εi，其中，

e.由下式修正ui和vi

ui＝ui+eivi＝vi+εi

f.由下式计算加速度ai

公式中所涉及的系数矩阵和向量分别为

优选的，所述的第一步当中，单元刚度矩阵为：

其中为空间杆件在局部坐标系中的单元刚度矩阵，是12×12的对称矩阵，对进行坐标变换，即可得到整体坐标系下的单元刚度矩阵刚度矩阵，即

优选的，所述的第一步当中，单元质量矩阵为：

采用hrz法对进行对角化，得到集中质量矩阵，即为

对进行坐标变换，即可得到整体坐标系下的单元刚度矩阵刚度矩阵，即

优选的，所述的第一步当中，采用rayleigh阻尼，建立单元阻尼矩阵

c^e＝a0m^e+a1k^e

其中

ωi和ωj一般分别取结构的第1阶和第3阶频率，ζ为阻尼比，一般为0.05。

优选的，所述的第一步当中，对于f的计算，需对每个时刻的地震波加速度记录采用线性插值得到üg(t)，由下式计算

f＝-müg(t)

优选的，所述的第二步当中，n取为4-10。

优选的，所述的第三步当中，方程右边的项可采用单元层级并行方法，以r1的计算为例进行说明，首先计算各个单元的此处可采用各个单元并行方法，然后由集成得到r1，其他方程右端项的计算采用同样方式处理；

本发明与现有技术相比，优点在于：

1)本发明的动力时程分析方法具有4阶精度，与newmark法、bathe方法等传统二阶动力时程分析方法相比，精度高2阶。

2)本发明的动力时程分析方法中，每个时间步仅需对维数neq(质点总自由度数目)的矩阵d11求一次逆，计算中能够保持空间有限元离散矩阵的带状稀疏性质，也可并行计算，计算量和newmark法等传统二阶动力时程分析方法相当，低于bathe方法的计算量，因此，对于大跨度桥梁结构，本发明的动力时程分析方法可以采用较长的时间步长，也即较少的时间步，得到与传统动力时程分析方法精度相当的结果，可大幅提升求解效率。

3)本发明的动力时程分析方法与bathe方法均具有渐进消去特性，为l-稳定方法，性能远远超越newmark法，可滤掉由于有限元离散引起的虚假振荡，保证计算稳定。

附图说明

图1某大跨悬索斜拉组合桥结构示意图；

图2有限元模型示意图；

图3空间梁单元示意图；

图4el-centro波示意图；

图5不同步长下bathe法和本发明动力时程分析方法计算的塔顶位移示意图；

图6不同步长下bathe法和本发明动力时程分析方法计算塔顶速度示意图；

图7不同步长下bathe法和本发明动力时程分析方法计算塔顶加速度示意图。

具体实施方式

下面结合具体实施例来对本发明进行进一步说明，但并不将本发明局限于这些具体实施方式。本领域技术人员应该认识到，本发明涵盖了权利要求书范围内所可能包括的所有备选方案、改进方案和等效方案。

下面结合附图对本发明的结构原理和工作原理作具体的描述：

以一个大跨桥梁结构为实例，具体阐述本发明动力时程分析方法，某大跨悬索斜拉组合桥如图1所示，此桥为自锚式斜拉—悬吊协作体系桥。不考虑引桥部分，该桥的跨径为：132m+400m+132m＝664m，具体结构构造为：(1)主梁：采用整体现浇单箱六室钢筋混凝土箱梁断面。主梁高2.6m，为主跨的1/155，顶板厚18-26cm，腹板40cm，底板25-26cm，人行道和钢筋锚固区布置在箱梁外侧。主梁端部加粗加厚，内部设置自锚锚碇。主梁中部120m为钢箱梁布置。(2)索塔：为a字形混凝土结构，桥面以上塔高70m，索塔从上向下逐渐加粗，塔顶面积30m²，塔根处为56.5m²，截面线性变大。塔顶做成圆形头，塔墩固结，塔梁分离。(3)斜拉索：索塔上采用扇形双索面，每个索面由17对拉索组成，非对称布置。边墩和辅助墩之间的斜拉索间距6m，其余部分间距8m，全桥4个索面，加上主塔处的0号索，一共为140根索。所有斜拉索均锚固在主梁内。斜拉索处主梁设置横隔梁。(4)主缆：因为主缆代替了边索的位置，所以顺桥向斜拉索布置并非完全对称。主缆锚固在边墩之上的主梁加粗处。(5)桥墩基础。初步材料表明，此处桥址物理状况比较好，可以采用端承桩，承台尺寸22.75、16.5m，承台为一整体，连接段尺寸为22.75×16.5m，采用的灌注桩，每承台24根桩。本申请的一种计算大跨度桥梁在地震作用下响应的无条件稳定、高阶精度、具有最佳数值阻尼的动力时程分析方法包含如下步骤：

第一步，对该大跨桥梁进行空间有限元离散，采用空间伯努利欧拉梁单元分别模拟该桥梁的主梁、索塔、斜拉索、主缆等构件，得到该桥梁的有限元模型，如图2所示，图中箭头方向为地震作用方向。空间梁单元是桥梁结构中最典型的构件，如图3所示。所述空间梁单元的单元位移向量为

所述空间梁单元的单元刚度矩阵为

其中为空间杆件在局部坐标系中的单元刚度矩阵，是12×12的对称矩阵。

对进行坐标变换，即可得到整体坐标系下的单元刚度矩阵，即

质量矩阵采用对角阵，为

对进行坐标变换，即可得到整体坐标系下的单元刚度矩阵，即

采用rayleigh阻尼，建立单元阻尼矩阵

c^e＝a0m^e+a1k^e

其中

对于一般结构。这里ωi和ωj一般分别取结构的第1阶和第3阶频率。ζ为阻尼比，为0.05。以上得到了单元刚度矩阵、质量矩阵和阻尼矩阵之后，集成整体刚度矩阵、质量矩阵和阻尼矩阵，由hamilton原理，可得到运动方程组如下

第二步，选取时间步长，时间步长取为加速度记录时间间隔的n倍，n取10；

第三步，逐时间步计算，对于第i个时间步，已知ti-1时刻的位移ui-1和速度vi-1，由以下几个子步骤，计算ti时刻的位移ui、速度vi；

a.由方程d11ui＝r1求解位移ui，其中，

r1＝p1+g11ui-1+δtg12vi-1

b.由方程d11vi＝r2求解位移vi，其中，

r2＝(p2+g21ui-1+δtg22vi-1-d21ui)/δt

c.由方程d11ei＝δr1求解位移的误差ei，其中，

d.由方程d11εi＝δr2求解速度的误差εi，其中，

e.由下式修正ui和vi

ui＝ui+eivi＝vi+εi

f.由下式计算加速度ai

公式中所涉及的系数矩阵和向量分别为

在基底输入如图4所示的elcentro地震波，该地震波的加速度记录间隔为0.01s，作用方向为垂直于桥梁平面。本发明的动力时程分析方法中，选取分析步长为加速度记录的10倍，即δt＝0.1s，计算塔顶位移响应。阻尼比选取0.05，各个质点的初始位移和初始速度均为0。

为了展示本发明的动力时程分析方法高效性和精确性，首先采用最常用的二阶精度newmark法进行分析，newark法的参数选为α＝1/2，β＝1/4，步长采用δt＝0.01s，分析失稳，计算结果趋于无穷大。进一步采用bathe方法进行分析，分析结果如图5-7所示，耗费时间为65s。仍然采用bathe法，尝试将步长放大原来10倍，即令δt＝0.1s，其结果已经背离δt＝0.01s时的结果，该结果说明，采用bathe法分析，若放大时间步长，会导致分析结果严重背离。而采用本发明的动力时程分析方法，令步长分别为δt＝0.01s和δt＝0.1s，其结果和采用bathe法，步长为δt＝0.01s的结果非常接近(3条响应曲线几乎重合)，如图5-7所示。该比较例非常直观的表明，当δt＝0.1s时，即本发明的动力时程分析方法当步长为bathe法的10倍时，结果仍然非常精确，耗费时间仅仅为20.4s，约为newmark法的1/3。充分验证了本发明大跨桥梁结构的动力时程分析方法的精确性和高效性。

应当理解的是，本发明描述的方法的步骤仅仅是示例性的描述，对其先后进行的时间顺序没有特殊的要求，除非其本身有必然的先后顺序关系。

如上所示，本发明虽然已参照有限的实施例和附图进行了说明，但在本发明所属领域中具备通常知识的人均可以从此记载中进行各种修改和变形。由此，其他实施例及权利要求书与等同物均属于权利要求的保护范围。