地震动力方程求解的几种新方法

地震动力方程求解的几种新方法

结构地震资料计算的主要方法是对多受地震区域的振动分解反应谱法进行分析。这种方法是一种静力分析。将地震剪力等效于平面力，然后通过静力学方法进行分析和计算。这种计算方法同实际地震反应尚有一定的差距,计算精度不够,不一定能够保证地震作用下的结构安全。时程分析法是一种动力分析法,它是将结构物视为一个弹性振动体,将地震时地面运动产生的位移、速度、加速度作用在结构物上,然后用动力学的方法研究其振动情况。显然,时程分析法比振型分解反应谱法能更准确地反映地震是结构物的反应。1区域弹塑性反应数值分析方法结构动力理论是直接通过动力方程求解地震反应,起源于20世纪60年代。由于地震波为复杂的随机振动,对于多自由度体系振动不可能直接得出解析解,只可采用逐步积分法,而这种方法计算工作量大,只有在计算机应用发展的前提下才能实现。多自由度体系地震反应方程为[Μ]{˙u}+[C]{˙u}+[Κ]{u}=-[Μ]{Ι}¨xg.(1)[M]{u˙}+[C]{u˙}+[K]{u}=−[M]{I}x¨g.(1)在式(1)中,地面振动加速度是复杂的随机函数。同时,在弹塑性反应中刚度矩阵与阻尼矩阵亦随时间变化。因此,不可能求出解析解,只能采取数值分析方法求解。常用的地震反应计算数值方法有线性加速度法、Newmark-β法、Wilson-θ法和中心差分法,将式(1)转化为增量方程为[Μ]{Δ˙u}+[C]{Δ˙u}+[Κ]{Δx}=-[Μ]{Δ¨ug}.(2)[M]{Δu˙}+[C]{Δu˙}+[K]{Δx}=−[M]{Δu¨g}.(2)再逐步积分求解,即将时间转化分成一系列微小时间段,在时间内可采取一些假设,从而能对增量式(2)直接积分,得出地震反应增量,以该步的终态值,作为下一时间段的初始值。这样逐步积分,即可得出结构在地震作用下振动反应的全过程。下面简单介绍这几种方法。2线性加速度法2.1基本思想1)假定在时间[t,t+Δt]内,加速度按线性变化。2)结构体系的特征在时间[t,t+Δt]内保持为常量。2.2根据增量方程求解将式(1)在时刻tj和tj+1应满足的方程相减可得到如下增量方程[Μ]{Δ¨u}j+[C]{Δ˙u}j+[Κ]{Δu}j=-[Μ]{Ι}Δxg,j.(3)[M]{Δu¨}j+[C]{Δu˙}j+[K]{Δu}j=−[M]{I}Δxg,j.(3)式中:{Δu}j={u}j+1-{u}j,其余以此类推。由于线性加速度法假定,在时段Δt内,结构的加速度反应是关于时间τ的线性函数。基于这一假定,可以将式(3)化为关于位移增量Δu的线性代数方程。为此,首先将{u}按Taylor级数在tj附近展开{u(tj+τ)}={u}j+{˙u}j1!τ+{¨u}j2!τ2+{⋯u}j3!τ3+⋯.(4){u(tj+τ)}={u}j+{u˙}j1!τ+{u¨}j2!τ2+{u⋯}j3!τ3+⋯.(4)对时间τ求导,可得{˙u(tj+τ)}={˙u}j+{¨u}jτ+12{⋯u}jτ2+⋯.(5){u˙(tj+τ)}={u˙}j+{u¨}jτ+12{u⋯}jτ2+⋯.(5)当τ=Δt时,由于{u(tj+τ)}={u}j+1和{˙u(tj+τ)}={˙u}j+1{u˙(tj+τ)}={u˙}j+1,并结合线性函数的假定,在求解过程中取为增量形式,则式(4)和式(5)可变为{Δu}j=6Δt2(Δu)j-6Δt{˙u}j-3{u}j,(6){Δ˙u}j={¨u}jΔt+12{Δ˙u}jΔt.(7){Δu}j=6Δt2(Δu)j−6Δt{u˙}j−3{u}j,(6){Δu˙}j={u¨}jΔt+12{Δu˙}jΔt.(7)将式(6)代入式(7)可得{Δ¨u}j=3Δt{Δu}j-3{˙u}j-12{u}jΔt.(8){Δu¨}j=3Δt{Δu}j−3{u˙}j−12{u}jΔt.(8)将式(6)和式(8)代入式(3)可得[ˉΚ]j{Δu}j={ΔΡ}j.(9)[K¯¯¯]j{Δu}j={ΔP}j.(9)其中[ˉΚ]j=6Δt2[Μ]+3Δt[C]+[Κ],(10){ΔΡ}j=[Μ](6Δt{˙u}j+3{u}j)+[C](3{˙u}j+Δt2{˙u}j)-[Μ]{Ι}Δxg,j.(11)[K¯¯¯]j=6Δt2[M]+3Δt[C]+[K],(10){ΔP}j=[M](6Δt{u˙}j+3{u}j)+[C](3{u˙}j+Δt2{u˙}j)−[M]{I}Δxg,j.(11)根据微分方程的初始条件和后续计算的过程可知,式(9)可以像静力问题那样求解。通常称式(9)为拟静力增量方程,而[ˉΚK¯¯¯]为拟静力刚度矩阵,{ΔP}为拟静力荷载向量。特别注意到,拟静力刚度矩阵[ˉΚK¯¯¯]不仅与刚度矩阵[K]有关,而且与质量矩阵[M]和阻尼矩阵[C]有关,这一点在后续的弹塑性动力分析中有重要的意义。同时,拟静力荷载向量{ΔP}不仅取决于地震地面运动加速度的增量,而且取决于前一时刻的计算反应值。这使得动力反应计算的误差逐渐积累,严重时甚至导致结果发散。为了尽可能减少这种误差,提出了加速度平衡校正算法,即根据增量动力平衡方程式求得{Δu}j=-{Ι}Δxg,j-[Μ]-1([C]{Δ˙u}j+[Κ]{Δu}j).(12){Δu}j=−{I}Δxg,j−[M]−1([C]{Δu˙}j+[K]{Δu}j).(12)上述推导过程是以增量方程为目的,这样推导出来的结果不仅能用于结构的弹性地震反应分析,而且也能够用于结构的弹塑性地震反应分析。当然,也可以全量方程为目的来推导相应于方程式的代数方程。事实上,式(6)和式(8)可以改写为{u}j+1=6Δt2{u}j+1-6Δt2{u}j-6Δt{˙u}j-2{u}j,(13){˙u}j+1=3Δt{u}j+1-3Δt{u}j-2{˙u}j-12{˙u}jΔt.(14){u}j+1=6Δt2{u}j+1−6Δt2{u}j−6Δt{u˙}j−2{u}j,(13){u˙}j+1=3Δt{u}j+1−3Δt{u}j−2{u˙}j−12{u˙}jΔt.(14)将式(13)和式(14)代入动力平衡方程可得[Κ′]{u}j+1={Ρ}j+1.(15)[K′]{u}j+1={P}j+1.(15)其中[Κ′]=[ˉΚ]j=6Δt2[Μ]+3Δt[C]+[Κ],(16){Ρ}j+1=[Μ](6Δt2{u}j+6Δt{˙u}j+2{u}j)+[C](3Δt{u}j+2{˙u}j+12{¨u}jΔt)-[Μ]{Ι}xg,j+1.(17)称式(15)为拟静力全量方程。对线性加速度算法而言,用增量方程与全量方程求解得到的结果,其计算精度是一样的。2.3收敛临界时间步长线性加速度法在选取时间步长时,应满足Δt<Tmin/a,这里Tmin是有限元离散系统中最小的固有周期。系数a一般取为10,如果Δt取得过大,计算得到的位移值可能会不收敛或者出现其他异常情况。但是,使结果收敛的临界时间步长是很难预先确定的。线性加速度法不光计算量大,而且实际上往往不能保证计算稳定性和精度。3新标记方法3.1控制积分格式的参数NewMark法是一种将线性加速度法普遍化的方法。该法假定位移和速度可表示为{u}j+1={u}j+{˙u}jΔt+(0.5-β){˙u}jΔt2+β{˙u}j+1Δt2,(18){˙u}j+1={˙u}j+(1+δ){u}jΔt+δ{u}j+1Δt.(19)其中,β和δ是控制积分格式计算精度和稳定性的参数。3.2对过阻尼时的情形在式(18)和式(19)中,当β和δ满足条件:δ≥0.5,β≥(0.5+δ)2/4时,NewMark法为无条件稳定的逐步积分格式;当δ=0.5时,NewMark法的计算精度为二阶,否则计算精度为一阶。当δ=0.5且β=1/6时,NewMark法为线性加速度法;当δ=0.5且β=0.25时,NewMark法为平均常加速度法。当δ≠0.5时,可导致系统的过阻尼情形。因此,一般取δ=0.5。对δ=0.5的情况,称为NewMark-β法。通常β的取值为1/6≤β≤1/2。当取δ=0.5时,由式(18)和式(19)可推出{Δu}j=1βΔt2{Δu}j-1βΔt{˙u}j-12β{˙u}j,(20){Δ˙u}j=12βΔt{Δu}j-12β{˙u}j+(1-14β)Δt{˙u}j.(21)将上述两式代入增量方程,可得[ˉΚ]{Δu}j={ΔΡ}j,(22)其中[ˉΚ]=1βΔt2[Μ]+12βΔt[C]+[Κ],(23){ΔΡ}j=[Μ](1βΔt{˙u}j+12β{¨u}j)+[C](12β{˙u}j-(1-14β)Δt{u}j)-[Μ]{Ι}Δ¨xg,j.(24)同样,也可以推导出全量方程的递推格式[Κ′]={u}j+1={Ρ}j+1.(25)其中[Κ′]=[ˉΚ]=1βΔt2[Μ]+12βΔt[C]+[Κ],(26){Ρ}j+1=-[Μ]{Ι}˙xg,j+1+[Μ][1βΔt2{u}j+1βΔt{˙u}j-(1-12β){¨u}j]+[C][12βΔt{u}j-(1-12β){˙u}j-(1-14β)Δt{¨u}j].(27)3.3物理模型的确定NewMark-β法是线性加速度法的推广。当δ≥0.5,β≥(0.5+δ)2/4,NewMark-β法就无条件稳定,即时间步长Δt大小可不影响解的稳定性。此时Δt主要根据解的精度确定,具体说可根据对结构响应有主要贡献若干基本振型的周期来确定。一般说基本振型周期比系统振型的最小振动周期大得多,所以无条件稳定的隐式算法比有条件稳定的显式算法可采用大得多的时间步长。而采用较大的Δt还可以滤掉高阶不精确特征解对系统响应的影响。4行政区域4.1n-加速度法为了改进线性加速度法有条件才能稳定计算的缺陷,得到无条件稳定的线性加速度法,wilson提出了一个简单而有效的wilson-θ法。该方法假定在时段θΔt内加速度随时间呈线性变化,其中θ>1。与线性加速度法的区别在于,线性加速度法在时刻t+Δt使用动力平衡方程,而wilson-θ法则将动力平衡方程应用于更后一点的时刻t+θΔt。4.2求解各运动方程{u(t+θΔt)}={u(t)}+θΔt{˙u(t)}+(θΔt)2{u(t)}/3+(θΔt)2{u(t+θΔt)}/6‚(28){˙u(t+θΔt)}={˙u(t)}+θΔt{u(t)}/2+θΔt{u(t+θΔt)}/2.(29)在t+θΔt时刻的运动方程为[Μ]{u(t+θΔt)}+[C]{˙u(t+θΔt)}+[Κ]{u(t+θΔt)}=-[Μ]{Ι}xg(t+θΔt).(30)由式(28)和式(29)导出的{u(t+θΔt)}和{˙u(t+θΔt)}代入式(30)可得[ˉΚ]{u(t+θΔt)}={Ρ(t+θΔt)}.(31)其中[ˉΚ]=6(θΔt)2[Μ]+3θΔt[C]+[Κ]‚(32){Ρ(t+θΔt)}=[Μ]{Ι}xg(t+θΔt)+[Μ](6(θΔt)2{u(t)}+6θΔt{˙u(t)}+2{˙u(t)})+[C](3θΔt{u(t)}+2{˙u(t)}+12θΔt{u(t)}).(33)将{u(t+θΔt)}代入式(29)求得{˙u(t+θΔt)},则t+θΔt时刻的加速度可按下式内插求得{˙u(t+θΔt)}=(1-1θ){u(t)}+1θ{¨u(t+θΔt)}.(34)4.3计算误差wilson-θ法的实质是线性加速度法推广。当θ≥1.37时,此法是无条件稳定的,但随着θ的增大,计算误差也增大,所以通常取θ=1.4。在地震作用下,对于一般阻尼比5%的钢筋混凝土结构,时间步长Δt≤0.04T(T为地震波的卓越周期)可以取得较好的结果。5分散法的中心5.1初始点函数值的比较中心差分法的基本思想是:在计算函数的中心点差分,并与初始点函数值进行比较,若两者之差足够小则结束,则以中心点差分函数值作为结果;否则,步长减半,并将中心点差分函数值送给初始点函数值,继续迭代,直到满足误差要求为止。5.23unt3+3.将位移增量函数按Taylor级数展开得u(t+Δt)=u(t)+˙u(t)Δt+12¨u(t)Δt2+16⋯u(t)Δt3+⋯,(35)u(t-Δt)=un-1,u(t)=un,u(t+Δt)=un+1.速度和加速度也有同样的关系,由式(35)得前差分式为un+1=un+˙unΔt+12¨unΔt2+16⋯unΔt3+⋯.(36)同样可得后差分公式为un-1=un-˙unΔt+12¨unΔt2-16⋯unΔt3+⋯.(37)把以上两式分别相减或相加可得˙unΔt=12(un+1-un-1)+Ο(Δt3),(38)¨unΔt2=(un+1-2un+un-1)+Ο(Δt4).(39)可用3点(n-1,n,n+1)的位移增量来近似表示t时刻的瞬时速度和加速度˙un=12Δt(un+1-un-1),(40)¨un=1Δt2(un+1-2un+un-1).(41)略去高阶微量O(Δt2),令t时刻(第n个时间步)的位移增量、速度和加速度满足该时刻动态显式有限元方法中的有限元微分方程,则[ˉΜ]¨un=Fextn-Fintn-Fen.(42)式中:M为对角矩阵,Fextn、Fintn和Fen是t时刻节