非平稳随机地震激励在ANSYS中快速模拟

侯中学

(四川省交通勘察设计研究院有限公司，四川 成都 610017)

1 概述

地震动中的非平稳性需要进行着重考虑[1-2]。在几乎所有的时程和随机抗震分析中，地震动的非平稳性用一个均匀调制随机过程来模拟。但是，均匀调制方法得到的随机过程只能表现出激励功率谱与时间变化的对应关系，而不能反映地震前期大量的高频成分以及高频成分迅速地衰减直至只剩低频成分的过程[3-5]，因此，在实现对地震动非平稳性模拟的方法中，用非均匀调制的随机过程来模拟真实地震动在时域和频率内的变化是合适的。

当前本领域研究人员的关注点主要放在寻找提高随机振动计算效率的方法上。在复杂结构的响应分析时考虑多维多点平稳随机地震激励，如果要考虑均匀调制的非平稳随机地震动激励，则在编程建模、计算方法、计算硬件等方面都有很高要求，在这些前提下结构的响应计算已经十分困难，如果再考虑非均匀调制模拟的非平稳随机地震动激励，那么复杂结构的响应分析就会更加费力[6-10]。同时，由于随机振动理论非常复杂，一般情况下，需要计算者自编程序来实现，这样很难被工程师们所接受。不少学者对于随机振动的计算效率和是否需要编程实现等问题进行研究，也取得了不少的成果[11-18]。大跨复杂结构遭受由非均匀调制方法得到的多维多点非平稳随机地震激励下的动力响应结果需要进一步研究，下面将使用随机振动法解决以上问题，通过精细积分法，可以提高计算效率的同时不需编制专门程序(完全通过有限元软件，如ANSYS完成)，此方法对地震动特性有充分的考虑。

为了将非平稳随机振动理论应用到实际工程抗震分析中去，本文在上述研究的基础上，首先，多维多点非平稳虚拟地震动激励法与精细积分法进行结合的计算方法，从计算效果上大大提高非均匀调制的非平稳地震进行激励时结构响应的计算速度和精准程度；给出单自由度体系由非平稳激励计算得到的理论解，同时检验该方法的精准程度；其次，将本文方法在ANSYS中进行快速模拟，节约了自编程序的时间，便于向实际工程应用推广；最后，以实际高墩桥梁为例，对该高墩桥梁在多维多点非均匀调制模拟方法得到的非平稳地震动激励下的结构响应进行分析，完成计算效果的验证。

2 非均匀调制非平稳地震的虚拟激励力的构造

2.1 非平稳激励地震动场模拟

大跨复杂结构m个支撑处受到3个方向激励的多维多点非平稳地震动场功率谱密度函数表示为：

(1)

其中任意子矩阵Skl(iω,t)是两个水平分量(x,y)和一个竖向分量(z)的3×3矩阵,具体表示如下：

(2)

两水平方向的功率谱密度函数被认为是相同的,而竖向与水平向的相干系数取0.6,因此表示如下：

Sklxx(iω,t)=Sklyy(iω,t)=Sklxy(iω,t)=Sklyx(iω,t)

(3)

(4)

互功率谱密度矩阵中的元素可表示为:

(5)

其中,Skkxx(ω,t),Skkxx(iω)均为第k点上x分量的自功率谱密度函数；ρklxx(iω)为k,l两点上x分量的相干函数。

2.2 直接求解虚拟激励法虚拟力构造

将式(1)中互功率谱矩阵进行分解得到：

S0(iω,t)=P*PT=[A(ω,t)][V]*[SCP]{q0}{q0}T[SCP][V]T[A(ω,t)]T

(6)

其中,P为3m×r的虚拟力矩阵,r为矩阵S0(iω,t)的秩,上标“*”和“T”分别为复数的共轭和转置。

SCP为3 m×3 m维矩阵,表示如下：

(7)

其中，SCPmx(ω)为第m个点x分量的功率谱密度函数。

V为3 m×3 m维行波因子矩阵,表示如下:

[V]=diag[e-iωT1x,e-iωT1y,e-iωT1z,e-iωT2x,e-iωT2y,
e-iωT2z,…,e-iωTmx,e-iωTmy,e-iωTmz]

(8)

其中，Tmx为空间地面运动x分量地震波传播第m点需要的时间;V为行波效应。

根据空间变化地面运动相干函数的分解方法,可将完全相干、部分相干和完全不相干虚拟力表示如下：

(9)

2.3 基于直接求解虚拟激励法的结构响应计算

(10)

其中:

(11)

将每一个特征值下得到的响应进行叠加就得到了部分相干情况下的响应功率谱密度如下：

(12)

完全相干:

(13)

式中:

(14)

完全不相干:

(15)

式中:

(16)

根据式(11)，式(14)，式(16),可以发现结构的每一个自振周期都对应着一个响应,那么就要求每次达到频率点时,都进行瞬态分析计算,虽然虚拟激励法可以直接求解,相对于使用较多的随机振动理论已经简化了计算过程,但是如果处理时离散的频率点数量太多,那么意味着结构需要计算的瞬态分析也就越多,计算的过程也就会变得复杂得多,采用将多维多点的非平稳虚拟激励法和精细化积分方法结合起来的结构仅需要计算两个瞬态分析就能确定结构任意节点的响应结果,可以有效简化计算,减少计算时间。

3 精细积分理论

将每一个确定频点下的动力方程式写成状态方程形式[10-11]，具体如下：

(17)

(18)

状态方程(17)的一般解为:

(19)

其中，eHt为指数矩阵,对式(19)进行离散积分,设时间步长为Δt,由递推法可知,假如在ti时刻的响应为V(ti)=Vi,那么ti-1时刻的响应V(ti-1)=Vi-1表示如下：

(20)

其中，T=eHΔt为指数矩阵,据文献[10-11]精细算法为:

T=eHΔt=[eHΔt/m]m

(21)

令τ=Δt/m,当取m=2N很大时,Δτ非常小,由泰勒级数展开式有:

(22)

当展开的项数N=20时,将泰勒级数节段后的截断误差与计算机进行计算的舍入误差相比,泰勒级数的截断误差要小很多,由以上对比可知,选取展开级数N=20不会因为截断造成过大数值误差,所以以上T(τ)的计算精度实际上已经比计算机的精确解要高出许多,用泰勒计算展开计算是合理的。

[I+Tai]≡[I+Ta,i-1]2=[I+2Ta,i-1+Ta,i-1×Ta,i-1]
(i=1,2,…,N)

(23)

因此,依次类推：

[I+TaN]≡[I+Ta,N-1]2=…=[I+Ta0]m=T(τ)

(24)

每次运算单元矩阵式时[I]不参与计算,为了避免因为Tai的值过小,以至于在程序中被默认忽略的情况,使用下列递推方法进行求解：

[Tai]=2[Ta,i-1]+
[Ta,i-1][Ta,i-1] (i=1,2,…,N)

(25)

假设在每个时间积分的区间(ti-1,ti)中,荷载值是线性变换的,那么非平稳虚拟力从时间的角度上将F(ω,t)进行等效离散,即分别在t0,t1,t2,…,tk处的随机变量F0,F1,F2…并用函数表示出来，Fk为非平稳虚拟力按照以下公式表示为:

F(ω,t)=r0+r1(τ-ti-1)

(26)

其中,r0,r1均为不变向量。

将式(26)代入式(20)可得：

Vi=T·Vi-1-H-1[H-1(Fi-Fi-1)/Δt+Fi]+
TH-1[H-1(Fi-Fi-1)/Δt+Fi-1](i=1,2,…,k)

(27)

其中H-1的计算可根据式(18)可得:

(28)

根据在每一个固定频率点上,V0=0,则可得到固定频率点各时刻的响应表达式为:

(29)

其中：

S3=TS2+S1

(30)

若用Ai,0,Ai,1,…,Ai,i分别表示F0,F1,…，Fi前面的系数,Ai,i下标的第一个i为ti时刻,第二个i为第i个集中荷载点的位置,则式(30)可简化表示为：

Vi=Ai,0F0+Ai,1F1+…Ai,i-1Fi-1+Ai,iFi(i=1,2,…,k)

(31)

其中,Ai,0,Ai,1,…,Ai,i是一种和结构自身固有属性相关和外界激励无关的变量,但这种固有属性是在结构本身受到外界荷载激励下所体现出来的。Ai,i为通过递推方法进行计算的,用Vi-1的系数Ai-1,0,Ai-1,1,…,Ai-1,i-1来依次推导下一步的系数为：

(32)

在系数矩阵中所有其他元素都是该系数矩阵的前两列元素的重复,同时具有明显的规律性。所以只需要矩阵前两列的元素,便可以确定整个系数矩阵。令F0=1,Fk=0(k>0),进行一个瞬态分析便可以得到系数矩阵的第一列所有元素。同样的,令F0=0,F1=1,Fk=1(k>2),在完成一次瞬态分析以后得到系数矩阵里的第二列元素。所以只需要两个条件——两个频率点的瞬态分析就可以得到完整的系数矩阵,在得到完整的系数矩阵后,结构中任意点在不同频率下的非平稳时程荷载作用后的响应便可以求出。在运用过程中,可以在三个方向按照时间激励先后顺序系分别对结构进行激励,可以求出结构关键点响应的系数矩阵,比如在运用绝对位移直接求解时,分别求出结构关键点响应的系数矩阵A,将施加在结构上的各个方向的激励时程与相对应的系数矩阵A相乘后进行叠加,求出以上固定频率点的响应结果。那么所有频率点都按上述方法代入后,即可得到全部频率下的结构响应。这种方法在一定程度上可以消除对每个频率点都做多维多点的时程分析的必要性,只需通过两个频率点可计算出整个结构的响应结果,缩短了运算过程的步骤和计算时间。

4 工程算例

4.1 计算模型介绍及参数说明

以某高墩桥梁为实际工程算例,图1为该算例中建立的三维有限元模型,采用梁墩固结的形式,墩底受到全部自由度上的约束,不过墩底全部平动方向的自由度在加载时会释放,1号、2号、4号和5号墩的墩顶与梁底是完全耦合的,仅有3号墩的墩顶和梁底在纵向自由度不耦合,桥台处约束了竖向位移和横竖扭转。

本文选用Clough-Penzien功率谱模型,参数谱强度因子S*=0.001 77 m2/s3,场地参数ωg=15.0 rad/s,ξg=0.6,ωf=1.5 rad/s,ξf=0.6,g(t)函数选用文献[14]的形式及其参数。

4.2 结果分析

由于篇幅有限,且本文的目的在于介绍非均匀调制非平稳随机激励的快速模拟计算理论,故这里仅给出2号桥墩受地震作用时的墩顶横向剪力响应时变功率谱和剪力响应均方差如图2，图3所示。

从图3可看出,响应功率谱同时随时间和频率发生变化,说明地震动所具有的非平稳特性,图4从时变均方差的角度也可进行佐证,因此说明该方法能够有效模拟非均匀调制非平稳随机地震激励。

5 结论

通过理论研究和计算验证得出以下结论:

1)本文方法并未专门编制运算程序,节约了编程时间,可高效地分析多维多点非均匀调制非平稳随机地震激励下结构的响应结果。

2)提出了绝对位移直接求解与精细积分相结合的ANSYS快速模拟方法,将N(离散频率点个数)个离散频率点的时程分析转换为2×3×m(m为支撑节点个数)个时程分析,成百上千倍地提高了计算效率。

3)根据本文所述方法,既可以考虑完全非平稳性的地震动,也能直接考虑地震动的场地效应、相干效应和行波效应的影响,与传统意义上的虚拟激励法具有相同的效果。