地震碰撞分析中接触单元模型的精细积分解法

张瑞杰, 李青宁, 尹俊红

(西安建筑科技大学 土木工程学院，西安　710055)



地震碰撞分析中接触单元模型的精细积分解法

张瑞杰, 李青宁, 尹俊红

(西安建筑科技大学 土木工程学院，西安710055)

摘要：地震碰撞数值分析是结构抗震研究的重要内容，将精细积分法引入到地震碰撞数值分析中。结合5种常用接触单元模型的特点，将碰撞力计算式表达成统一的形式，推导了碰撞阶段的精细积分公式。将地震碰撞过程划分为碰撞阶段和非碰撞阶段，对不同的阶段采用不同的积分步长，并且提出一种基于线性加速度假定的碰撞时刻搜索法。用两小球自由弹性碰撞模型对精细积分法程序进行验证。并采用本文方法对Chau振动台碰撞试验进行了数值仿真。结论为：数值分析方法既保证了模拟地震碰撞的精度又具有较高的计算效率；在地震碰撞分析中应合理选择接触单元模型和模型参数；对于近似弹性碰撞情况Hertz模型和Hertz-damp模型数值仿真结果与试验结果吻合较好。

关键词：地震;结构碰撞;数值分析;精细积分法;接触单元

震害调查表明，地震时相邻结构或构件由于动力特性的差异导致的碰撞会引起局部破坏、落梁，甚至倒塌[1-3]。要揭示地震中桥梁结构的碰撞机理、进行碰撞研究及抗震设计，需要精确地模拟桥梁在地震中的碰撞行为。国内外众多学者已相继开展了地震碰撞试验及数值分析研究[4-6]。目前，地震碰撞数值分析大多采用接触单元法[7-8]，该方法是在普通有限元时程积分基础上发展起来的，通过在可能碰撞的结构体之间设置接触单元，当碰撞发生时激活接触单元，实现结构碰撞行为的模拟。迄今，已发展了多种接触单元模型[9-12]。

接触单元法本质上是将接触碰撞过程看作碰撞力的作用过程，在结构的动力平衡方程中加入碰撞力，得到考虑碰撞作用的动力平衡方程。动力平衡方程在时域内的数值求解有多种方法[13-14]。传统的显式积分法如中心差分法、隐式积分方法如Newmark-β法、Wilson-θ法都是将微分化成为差分算子，使得这些方法仅具有二阶精度或三阶精度，而且对时间步长非常敏感[14]。在碰撞动力平衡方程求解中，碰撞力是在当前阶段的位移和速度响应基础上计算得到，然后将其作为下一时间点的已知力进行求解。这造成了碰撞力施加时间的滞后。因而在地震碰撞模拟中，只能采用较小的时间步长来减轻这种不良影响，但这大幅增加了计算量。

钟万勰[15]提出了一种求解结构动力响应的精细积分法。鲍四元等[16]将其应用于弹性撞击问题的动态响应求解，并指出精细积分法具有受时间步长限制小、精度高和无条件稳定的优点。但目前鲜有文献报道将精细积分法应用到地震碰撞过程的模拟中。本文将精细积分法引入到地震碰撞过程模拟中，结合5种常用接触单元模型的特点，将其碰撞力表达成统一形式，推导了精细积分公式。利用精细积分法对步长不敏感的特点，本文提出了一种变步长的措施，在没有接触碰撞时用较大步长，提高计算效率，在接触碰撞阶段用较小步长，以提高模拟地震碰撞的精度。编制了计算程序，通过对Chau振动台试验进行数值仿真验证了本文方法的有效性和优越性。

1接触单元模型

考察几种常用的接触单元模型，其模型假定和碰撞力特点见图1。线弹簧模型、线弹簧阻尼(Kelvin-Voigt)模型、非线性弹簧(Hertz)模型、非线性弹簧阻尼(Hertz-damp)模型在碰撞阶段的碰撞力表达式均与相对压入量及速度差有关，可写成如下统一形式：

(1)

当u1(t)-u2(t)-gp≥0

图1　接触单元模型Fig.1 Model of contact elements

(2)

式中：m1、m2分别为碰撞体的质量；ζ为与恢复系数e有关的碰撞阻尼常数，当e=0时表示两碰撞体之间发生完全塑性碰撞，e=1时表示两碰撞体之间发生完全弹性碰撞。

线弹簧阻尼接触单元假设在碰撞期间均存在能量耗散。由于阻尼系数ck为一常数，因而在接触开始时刻会出现突加的阻尼力，而在恢复过程中，阻尼力可能超过弹性恢复力，引起碰撞力为负值的情况，与实际情况不符。

(3)

(4)

(5)

2精细积分法求解碰撞动力方程

对于图2所示的相邻的两个单自由度系统碰撞模型，其动力平衡方程表达式为：

(6)

图2　碰撞结构模型Fig.2 Model of colliding structures

方程(6)分两种情况进行讨论，其一在非接触阶段，Fc(t)=0，即为两个独立的动力平衡方程；其二，在接触碰撞阶段Fc(t)采用公式(1)计算。下面详细推导接触碰撞阶段的精细积分法公式。

将公式(1)代入方程(6)，将碰撞力项中的弹性恢复力在左端进行分解并与结构刚度矩阵合并，而阻尼力仍作为非齐次项处理。可写成如下形式：

(7)

对方程(7)，引入变量P，令

(8)

则有，

(9)

(10)

方程(7)可转化为Hamilton体系下的一阶常微分方程。

(11)

其中系统的转换矩阵H为：

(12)

式(11)为结构碰撞阶段动态响应的状态方程。其一般解为：

(13)

假定接触开始的时间点是t0时刻，Vt0是非接触结束时刻的状态向量。将接触碰撞阶段分成时间步长为Δt的若干时间间隔， 则tk=t0+kΔt(k=0,1,2,…)，而tk+1=tk+Δt。可以认为Δt时段内荷载按线性变化，满足：F(τ)=F(tk)+[F(tk+1)-F(tk)](τ-tk)/Δt=F0+F1(τ-tk)，根据式(11)的解答可得V(tk)、V(tk+1)之间的转换关系。

V(tk+1)=T[V(tk)+H-1(F0+H-1F1)]-

H-1[F0+H-1F1+F1Δt]

(14)

式中T为指数矩阵：

(15)

Ta矩阵采用泰勒级数展开求解：

(16)

式中：l为截断阶数，在碰撞问题求解中可取l=5,N=20[14-15]。求出Ta矩阵后，由于Ta矩阵中的元素非常小，若直接代入式(15)求矩阵T，则由于计算机的舍入误差会使得指数矩阵T的精度丧失殆尽。利用指数函数的加法原理，指数矩阵T的求解最终采用下式进行。

T=I+Ta

(17)

对于非接触碰撞阶段的动力方程，同样可以通过“增元降阶”的方式转化成如(11)式的状态方程。此阶段，系统的转换矩阵H为：

非齐次项：

(19)

利用式(15)~式(17)可求得非接触阶段的指数矩阵T，代入式 (14)可递推任意非碰撞时刻的状态向量，进而求解速度和加速度。值得注意的是，除时刻为0时初始状态为已知条件，其余非接触阶段的初始状态是接触阶段结束时的状态。

综上所述，将接触单元模型的弹性恢复力作为未知量处理，通过分解合并得到接触碰撞阶段的动力平衡方程。之所以将阻尼力作为已知量处理，采用前一时刻的值近似代替，是因为在Kelvin-Voigt模型和Jan-Hertz-damp模型中阻尼力是突变的，将其作为未知量在左端分解，会引起数值计算的不稳定。对于线弹簧模型，本文解法为精确解；对于Kelvin-Voigt模型、Hertz模型、Hertz-damp模型、Jan-Hertz-damp模型，碰撞弹簧刚度以及阻尼力是时间相关函数，其tk+1时刻值由tk时刻的值近似代替。由于Δt很小，这种近似引起的误差也很小。精细积分法的求解分为两个阶段进行。对于不同的阶段可以采用不同的积分步长。为进一步提高计算精度，需要精准的找到开始碰撞的时刻，为此，采用了线性加速度的假定。

采用积分步长为Δt，假定ti时刻u1-u2-gp<0(没有接触)，而ti+1时刻u1-u2-gp>0(已经接触)，说明接触开始时刻ti+τ介于ti时刻与ti+1时刻之间。通过精细积分法已解出ti时刻和ti+1时刻的位移、速度、加速度。假定ti时刻到ti+1时刻加速度按线性变化，则位移变化为三次曲线形式。ti+τ时刻的位移表达式为：

(20)

假定ti+τ开始接触，即满足u1-u2-gp=0条件：

(21)

3算例验证

3.1理论验证

通过对两小球自由弹性碰撞的数值模拟，验证精细积分法在求解碰撞问题时的精确性和高效率。两小球质量均为2 kg，刚度均为210.25 N/m，左侧小球初始位移0.1 m，释放后自左向右运动与右侧小球发生碰撞，不考虑系统的能量损失，碰撞模型[5]见图3(a)。该模型满足动量守恒的理论位移时程见图3(b)。在用精细积分法求解时，接触单元采用线性弹簧模型，碰撞弹簧刚度参数取1.547 4×106N/m。

图4~图6给出了Matlab平台上的精细积分法、商业软件Ansys 平台上的Newmark法、Matlab平台上Willson-θ法计算的前2 s位移响应时程和加速度响应时程，三种算法的分析过程及结果对比列于表1。

图3　两质点自由弹性碰撞模型及其理论解[5]Fig.3 An illustration of two point masses on free vibration and its theoretical solution

图4　精细积分法求解位移时程及加速度时程Fig.4 Time history of displacement and acceleration with PIM

图5 Newmark法求解位移时程及加速度时程Fig.5 Time history of displacement and acceleration with Ansys method

图6　Willson-θ法求解位移时程及加速度时程Fig.6 Time history of displacement and acceleration with Willson-θ method

算法步长/s碰撞非碰撞计算时间/s加速度峰值/(m·s-2)质点2位移峰值/m质点2精细积分0.00010.0011.2638.40.100Newmark0.00010.001313637.40.100Willson-θ0.000010.000011201709.30.106

与理论解相比，精细积分法位移时程的周期，幅值均与理论解一致，但是碰撞的位置与理论解有微小的差异，并且随着时间推移差异有增加的趋势。这是因为理论解基于动量守恒定律，碰撞是在瞬间完成的，而采用线性弹簧接触模型的碰撞力是有作用过程的，本例接触过程约为0.002 5 s，这一过程随着接触弹簧刚度的增大而减小，其位移解也更接近图5理论解。

总体上，精细积分法和Ansys 软件Newmark法计算的位移时程响应与满足动量守恒的理论值基本吻合，表明精细积分法和Ansys软件Newmark法均具有较高的计算精度。Wilson-θ法计算的前2 s位移峰值与理论解相差6%，加速度峰值相差11%，且位移峰值和加速度峰值有随时间增加的趋势，说明该方法计算精度相对较差。

精细积分法在非碰撞阶段的步长为0.001 s，碰撞阶段的步长为0.000 1 s，2 s内总积分点2145个；Ansys 软件Newmark法的积分步长与精细积分法相同，2 s内积分点数同样为2145个，但精细积分法总计算时间仅1.2 s，而在Ansys软件Newmark法需313 s，可知在保持相同计算精度的情况下，采用精细积分法效率得到很大提高。Wilson-θ法计算程序由于未采用变步长算法，当积分步长为0.000 01 s时，2 s内积分点数达到20 000个，计算时间为1 201 s，若积分步长取0.000 1 s，得到的位移时程及加速度时程是发散的，说明Wilson-θ法对积分步长非常敏感，不适用于接触过程短暂的结构碰撞分析。

3.2试验验证

Chau进行了地震激励下两钢塔的振动台碰撞试验[6]，文献[13]应用基于FENAP平台的隐式积分功能对该试验进行了数值仿真并对比了不同接触单元计算结果。本文基于精细积分法在Matlab平台上编制了地震碰撞分析程序，并针对Chau试验进行数值仿真。

两钢塔塔高均为2 m，1号钢塔质量m1=204.0 kg，基频f1=2.3 Hz，阻尼比ξ1=1.4%；2号钢塔质量为m1=146.5 kg，基频f1=2.9 Hz，阻尼比ξ1=1.6%；两钢塔间初始间隙为0.011 m。两钢塔固定在振动台上，沿纵向输入1940年El-Centro地震动记录NS分量。该地震动记录持时53.73 s，峰值加速度为0.348 g。分析中将钢塔简化为两个单自由度体系，刚度分别为k1=4 297.5 N/m，k2=4 864.0 N/m，两钢塔的碰撞刚度参数及阻尼参数见表2。

表2　两钢塔接触模型的参数[6,13]

数值分析中，非接触阶段精细积分时间步长为0.01 s，接触碰撞阶段的积分时间步长为0.000 01 s。

图8为Chau试验所测得的1号钢塔速度响应时程。图9(a)~(e)为基于本文精细积分法采用5种接触单元模型计算的1号钢塔速度响应时程。对比图8、图9发现，5种接触单元模型计算的1号钢塔速度响应时程与试验结果趋势吻合，但数值分析结果的速度响应峰值偏小约25%。文献[6,13]数值分析结果也比试验结果明显偏小。这是由于振动台试验虽然输入地震动记录是峰值0.348 g的El-Centro波，但台面实际加速度会产生变异。若在数值分析中采用实测振动台台面加速度时程，数值分析结果与试验结果吻合程度会更高。

图7　两钢塔碰撞振动台试验[6,13]Fig.7 Shaking table test on pounding of two steel towers

图8　试验测得1号钢塔速度响应时程[6]Fig.8 Tested time history of velocity of tower 1

图9　精细积分法计算1号钢塔速度响应时程Fig.9 Calculated time history of velocity of tower 1 by PIM

图10(a)~(e)给出5种接触单元模型计算的相对压入量-碰撞力关系曲线。图10(a)碰撞力与相对压入量呈线性关系，符合线弹簧接触单元模型假定，斜率为接触弹簧刚度；图10(b)碰撞力包含弹性恢复力和阻尼力两项，在接近阶段碰撞力大于弹性恢复力，且在接触瞬间出现突加碰撞力，而恢复阶段碰撞力小于弹性恢复力，且在碰撞结束时出现负值，这符合Kelvin-Voigt模型特点；图10(c)碰撞力与相对压入量关系为1.5次抛物线，接近阶段与恢复阶段碰撞力路径相同，符合Hertz模型假定；图10(d)、(e)分别符合Hertz-damp和Jan-Hertz-damp模型假定。两种模型均是在Hertz模型的基础上引入阻尼项，因而接近阶段与恢复阶段碰撞力路径不再重合。但不同的阻尼项假定所呈现的阻尼力特点有差别，Hertz-damp模型阻尼力随着相对压入量逐渐增加，在接近最大压入量时迅速减小为0，在恢复过程阻尼力为负值；而Jan-Hertz-damp模型阻尼力呈现随着相对压入量快速增大，然后基本保持稳定，在接近最大压入量时逐渐减小为0，在恢复过程不出现阻尼力的特点。图10(a)~(e)计算结果与理论假定吻合良好，说明采用本文精细积分法能精确的模拟符合不同接触单元模型假定的碰撞过程。

计算发现，采用线弹簧、Kelvin-Voigt模型计算的碰撞时间约为0.77 ms，采用Hertz模型、Hertz-damp模型、Jan-Hertz-damp模型碰撞时间为1.89~1.98 ms。对于线弹簧、Kelvin-Voigt模型本文选取的碰撞阶段的积分步长约为碰撞历程的1/77；对于Hertz、Hertz-damp、Jan-Hertz-damp模型积分步长约为碰撞历程的1/200。在碰撞阶段采用很小的积分步长保证了碰撞过程模拟的精度，而在非碰撞阶段采用的积分步长是碰撞阶段积分步长的1 000倍，极大减少了计算工作量，保证了计算效率。

图10　碰撞力-相对压入量关系曲线Fig.6 Curve of pounding force varying with penetration displacement

表3给出基于精细积分法的5种接触单元模型数值计算结果与试验结果的对比。从表3分析发现，线弹簧模型和Kelvin-Voigt模型碰撞力峰值比试验结果明显偏高，计算结果不合理；Hertz、Hertz-damp、Jan-Hertz-damp模型碰撞力峰值与试验结果偏差分别为18.7%、2.3%、-10.6%，Hertz-damp模型计算结果最为接近；Hertz模型和Hertz-damp模型计算的碰撞次数与实测结果最为接近；从超过60 kN碰撞次数看，线弹簧、Kelvin-Voigt模型计算明显不合理，而Hertz模型与试验结果最为接近。这是因为Chau所做的振动台碰撞试验接近弹性碰撞，因而与Hertz模型计算结果吻合较好。

以上分析表明，本文提出的精细积分法应用于地震碰撞模拟是可行的，且具有较高的精度。为准确模拟地震碰撞，合理选择接触单元模型并且确定模型的参数是关键。在近似弹性碰撞时Hertz、Hertz-damp模型具有较好的效果。在相同的刚度参数和恢复系数下，Jan-Hertz-damp模型计算的碰撞力最小，说明其能量耗散大于Hertz-damp模型。

表3　数值计算结果与试验结果对比

4结论

本文将5种常用接触单元模型碰撞力表达成统一形式，并与精细积分法相结合，推导了碰撞阶段的精细积分公式。在Matlab环境下编制了计算程序，进行了理论解验证和试验验证，主要结论如下：

(1)精细积分法是一种显式积分方法，对积分时间步长不敏感，在求解时将结构运动过程分为碰撞阶段和非碰撞阶段，在非碰撞阶段求解时采用较大的步长，而在碰撞阶段采用较小的步长，不但保证了碰撞模拟的精度且获得了较高的计算效率，适于结果碰撞过程的模拟。

(2)基于线性加速度假设的碰撞时刻搜索法，可以精确求解碰撞时刻，从而进行碰撞阶段和非碰撞阶段的求解转换。

(3)在地震碰撞模拟中，应合理选择接触单元模型和模型参数。本文采用Hertz模型及Hertz-damp模型计算结果与试验结果较为吻合，说明在近似弹性碰撞时，这两种模型具有较高的模拟精度。

参 考 文 献

[1] Rosenblueth E, Meli R. The 1985 earthquake: causes and effects in Mexico City[C]//Concrete International, 1986, 8: 23-34.

[2] Priestley M J N, Seible F, Calvi G M. Seismic design and retrofit of bridge [M]. USA, New York: John Wiley and Sons INC, 1996: 8-20.

[3] Kasai K, Maison B F. Building pounding damage during the 1989 Loma Prieta earthquake. Engineering Structures 1997;19:195-207.

[4] Jankowski R. Pounding force response spectrum under earthquake excitation.[J]. Engineering Structures, 2006, 28(8): 1149-1161.

[5] Zhu Ping, Abe M, Fujino Y. Modelling three-dimensional non-linear seismic performance of elevated bridges with emphasis on pounding of girders[J].Earthquake Engineering and Structural Dynamics, 2002, 31:1891-1913.

[6] Chau K T, Wei X X, Guo X. Experimental and theoretical simulations of seismic poundings between two adjacent structures [J]. Earthquake Engineering and Structural Dynamics, 2003, 32:537-554.

[7] Chau K T, Wei X X. Pounding of structures modeled as non-linear impacts of two oscillators [J]. Earthquake Engineering and Structural Dynamics, 2001, 30:633-651.

[8] Jankowski R. Non-linear viscoelastic modeling of earthquake-induced structural pounding [J]. Earthquake Engineering and Structural Dynamics, 2005, 34(6):595-611.

[9] Zanardo G, Hao H, Modena C. Seismic response of multi-span simply supported bridges to a spatially varying earthquake ground motion [J]. Earthquake Engineering and Structural Dynamics, 2002, 31: 1325-1345.

[10] Anagnostopoulos S A. Equivalent viscous damping for modeling inelastic impacts in earthquake pounding problems [J]. Earthquake Engineering and Structural Dynamics, 2004, 33(8): 897-902.

[11] Van Mier J G M, Pruijssers A F, Reinhardt H W, et al. Load-time response of colliding concrete bodies [J].Journal of Structural Engineering (ASCE), 1991, 117(2):354-374.

[12] Muthukumar S, DesRoches R. A Hertz contact model with non-linear damping for pounding simulation [J].Earthquake Engineering and Structural Dynamics, 2006,35(7): 811-828.

[13] 禚一,李忠献,王菲. 桥梁地震碰撞分析中不同接触单元模型的对比分析[J].工程力学, 2014, 31(3): 11-17.

ZHUO Yi, LI Zhong-xian,WANG Fei. Comparative analysis of different contact element models in seismic pounding analysis of bridges. [J]. Engineering Mechanics, 2014, 31(2): 11-17.

[14] 郭泽英,李青宁,张守军. 结构地震发应分析的一种新精细积分法[J]. 工程力学, 2007, 24(4): 35-40.

GUO Ze-ying,LI Qing-ning,ZHANG Shou-jun. A new precise integration method for structure seismic response analysis.[J]. Engineering Mechanics, 2007,24(4):35-40.

[15] 钟万勰.结构动力方程的精细时程积分法[J].大连理工大学学报,1994, 34 (2): 131-136.

ZHONG Wan-xie. One precise time integration method for structural dynamics[J]. Journal of Dalian University of Technology, 1994, 34 (2): 131-136.

[16] 鲍四元,邓子辰. 结构撞击响应的一种弹性模型及其精细求解[J].工程力学,2008,25 (60):14-17.

BAO Si-yuan, DENG Zi-chen. An elastic model and its precise solution for structure impact response[J]. Engineering Mechanics, 2008,25 (60):14-17.

[17] 邹立华,郭润,黄凯,等.带预应力橡胶支座相邻隔震结构碰撞分析[J].振动与冲击,2014,33(9):131-136.

ZOU Li-hua, GUO Run, HUANG Kai, et al. Pounding of adjacent isolated-structures with prestressed rubber bearings[J]. Journal of Vibration and Shock,2014,33(9):131-136.

Precise integration method for the solution of contact element model in earthquake pounding analysis

ZHANGRui-jie,LIQing-ning,YINJun-hong

(School of Civil Engineering, Xi’an University of Architecture and Technology, Xi’an 710055, China)

Abstract:The numerical analysis of earthquake pounding is one of the important points in bridges’ seismic study. A precise integration method was introduced into the numerical analysis of the structures’ earthquake pounding. In consideration of the characteristics of five kinds of contact element models that are commonly used, a unified expression of pounding force was put forward, and then the formulas for the precise integration method were deduced. The structures’ response process in earthquake was divided into two parts, the one is collision process and the other is none collision process. And different integration time steps were applied in different processes as well as a method to search precise contact moment was developed. The precise integration method program has been tested by using a two lumped masses on free vibration model. And then the numerical simulation of two steel structures’ earthquake pounding experiments implemented by Chau was carried out based on the presented precise integration method. The conclusions are: the precise integration method not only can ensure the precision of the pounding simulation but also has high computing efficiency; the contact element model as well as the model’s parameters should be choosen carefully; the numerical analysis results by the Hertz model and Hertz-damp model are in good agreement with the experiment results in the case of approximate elastic pounding.

Key words:earthquake; structure pounding; numerical analysis; precise integration method; contact element model

中图分类号:TU973;U442.5

文献标志码：A

DOI:10.13465/j.cnki.jvs.2016.03.019

收稿日期：2014-09-10修改稿收到日期：2014-12-19

基金项目：国家自然科学基金(51078306)；高等学校博士学科点专项科研基金(20106120110004)；西安建筑科技大学重大科技项目创新基金(ZX0901)；陕西省重点学科建设专项资金资助项目(E01004)

第一作者 张瑞杰 男，博士生，讲师，1982年生