地震動力方程求解的幾種新方法

地震動力方程求解的幾種新方法

結(jié)構(gòu)地震資料計算的主要方法是對多受地震區(qū)域的振動分解反應(yīng)譜法進行分析。這種方法是一種靜力分析。將地震剪力等效于平面力，然后通過靜力學(xué)方法進行分析和計算。這種計算方法同實際地震反應(yīng)尚有一定的差距,計算精度不夠,不一定能夠保證地震作用下的結(jié)構(gòu)安全。時程分析法是一種動力分析法,它是將結(jié)構(gòu)物視為一個彈性振動體,將地震時地面運動產(chǎn)生的位移、速度、加速度作用在結(jié)構(gòu)物上,然后用動力學(xué)的方法研究其振動情況。顯然,時程分析法比振型分解反應(yīng)譜法能更準(zhǔn)確地反映地震是結(jié)構(gòu)物的反應(yīng)。1區(qū)域彈塑性反應(yīng)數(shù)值分析方法結(jié)構(gòu)動力理論是直接通過動力方程求解地震反應(yīng),起源于20世紀(jì)60年代。由于地震波為復(fù)雜的隨機振動,對于多自由度體系振動不可能直接得出解析解,只可采用逐步積分法,而這種方法計算工作量大,只有在計算機應(yīng)用發(fā)展的前提下才能實現(xiàn)。多自由度體系地震反應(yīng)方程為[Μ]{˙u}+[C]{˙u}+[Κ]{u}=-[Μ]{Ι}¨xg.(1)[M]{u˙}+[C]{u˙}+[K]{u}=?[M]{I}x¨g.(1)在式(1)中,地面振動加速度是復(fù)雜的隨機函數(shù)。同時,在彈塑性反應(yīng)中剛度矩陣與阻尼矩陣亦隨時間變化。因此,不可能求出解析解,只能采取數(shù)值分析方法求解。常用的地震反應(yīng)計算數(shù)值方法有線性加速度法、Newmark-β法、Wilson-θ法和中心差分法,將式(1)轉(zhuǎn)化為增量方程為[Μ]{Δ˙u}+[C]{Δ˙u}+[Κ]{Δx}=-[Μ]{Δ¨ug}.(2)[M]{Δu˙}+[C]{Δu˙}+[K]{Δx}=?[M]{Δu¨g}.(2)再逐步積分求解,即將時間轉(zhuǎn)化分成一系列微小時間段,在時間內(nèi)可采取一些假設(shè),從而能對增量式(2)直接積分,得出地震反應(yīng)增量,以該步的終態(tài)值,作為下一時間段的初始值。這樣逐步積分,即可得出結(jié)構(gòu)在地震作用下振動反應(yīng)的全過程。下面簡單介紹這幾種方法。2線性加速度法2.1基本思想1)假定在時間[t,t+Δt]內(nèi),加速度按線性變化。2)結(jié)構(gòu)體系的特征在時間[t,t+Δt]內(nèi)保持為常量。2.2根據(jù)增量方程求解將式(1)在時刻tj和tj+1應(yīng)滿足的方程相減可得到如下增量方程[Μ]{Δ¨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內(nèi),結(jié)構(gòu)的加速度反應(yīng)是關(guān)于時間τ的線性函數(shù)?；谶@一假定,可以將式(3)化為關(guān)于位移增量Δu的線性代數(shù)方程。為此,首先將{u}按Taylor級數(shù)在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)對時間τ求導(dǎo),可得{˙u(tj+τ)}={˙u}j+{¨u}jτ+12{?u}jτ2+?.(5){u˙(tj+τ)}={u˙}j+{u¨}jτ+12{u?}jτ2+?.(5)當(dāng)τ=Δt時,由于{u(tj+τ)}={u}j+1和{˙u(tj+τ)}={˙u}j+1{u˙(tj+τ)}={u˙}j+1,并結(jié)合線性函數(shù)的假定,在求解過程中取為增量形式,則式(4)和式(5)可變?yōu)閧Δ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)根據(jù)微分方程的初始條件和后續(xù)計算的過程可知,式(9)可以像靜力問題那樣求解。通常稱式(9)為擬靜力增量方程,而[ˉΚKˉˉˉ]為擬靜力剛度矩陣,{ΔP}為擬靜力荷載向量。特別注意到,擬靜力剛度矩陣[ˉΚKˉˉˉ]不僅與剛度矩陣[K]有關(guān),而且與質(zhì)量矩陣[M]和阻尼矩陣[C]有關(guān),這一點在后續(xù)的彈塑性動力分析中有重要的意義。同時,擬靜力荷載向量{ΔP}不僅取決于地震地面運動加速度的增量,而且取決于前一時刻的計算反應(yīng)值。這使得動力反應(yīng)計算的誤差逐漸積累,嚴(yán)重時甚至導(dǎo)致結(jié)果發(fā)散。為了盡可能減少這種誤差,提出了加速度平衡校正算法,即根據(jù)增量動力平衡方程式求得{Δ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)上述推導(dǎo)過程是以增量方程為目的,這樣推導(dǎo)出來的結(jié)果不僅能用于結(jié)構(gòu)的彈性地震反應(yīng)分析,而且也能夠用于結(jié)構(gòu)的彈塑性地震反應(yīng)分析。當(dāng)然,也可以全量方程為目的來推導(dǎo)相應(yīng)于方程式的代數(shù)方程。事實上,式(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)為擬靜力全量方程。對線性加速度算法而言,用增量方程與全量方程求解得到的結(jié)果,其計算精度是一樣的。2.3收斂臨界時間步長線性加速度法在選取時間步長時,應(yīng)滿足Δt<Tmin/a,這里Tmin是有限元離散系統(tǒng)中最小的固有周期。系數(shù)a一般取為10,如果Δt取得過大,計算得到的位移值可能會不收斂或者出現(xiàn)其他異常情況。但是,使結(jié)果收斂的臨界時間步長是很難預(yù)先確定的。線性加速度法不光計算量大,而且實際上往往不能保證計算穩(wěn)定性和精度。3新標(biāo)記方法3.1控制積分格式的參數(shù)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)其中,β和δ是控制積分格式計算精度和穩(wěn)定性的參數(shù)。3.2對過阻尼時的情形在式(18)和式(19)中,當(dāng)β和δ滿足條件:δ≥0.5,β≥(0.5+δ)2/4時,NewMark法為無條件穩(wěn)定的逐步積分格式;當(dāng)δ=0.5時,NewMark法的計算精度為二階,否則計算精度為一階。當(dāng)δ=0.5且β=1/6時,NewMark法為線性加速度法;當(dāng)δ=0.5且β=0.25時,NewMark法為平均常加速度法。當(dāng)δ≠0.5時,可導(dǎo)致系統(tǒng)的過阻尼情形。因此,一般取δ=0.5。對δ=0.5的情況,稱為NewMark-β法。通常β的取值為1/6≤β≤1/2。當(dāng)取δ=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)同樣,也可以推導(dǎo)出全量方程的遞推格式[Κ′]={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-β法是線性加速度法的推廣。當(dāng)δ≥0.5,β≥(0.5+δ)2/4,NewMark-β法就無條件穩(wěn)定,即時間步長Δt大小可不影響解的穩(wěn)定性。此時Δt主要根據(jù)解的精度確定,具體說可根據(jù)對結(jié)構(gòu)響應(yīng)有主要貢獻若干基本振型的周期來確定。一般說基本振型周期比系統(tǒng)振型的最小振動周期大得多,所以無條件穩(wěn)定的隱式算法比有條件穩(wěn)定的顯式算法可采用大得多的時間步長。而采用較大的Δt還可以濾掉高階不精確特征解對系統(tǒng)響應(yīng)的影響。4行政區(qū)域4.1n-加速度法為了改進線性加速度法有條件才能穩(wěn)定計算的缺陷,得到無條件穩(wěn)定的線性加速度法,wilson提出了一個簡單而有效的wilson-θ法。該方法假定在時段θΔt內(nèi)加速度隨時間呈線性變化,其中θ>1。與線性加速度法的區(qū)別在于,線性加速度法在時刻t+Δt使用動力平衡方程,而wilson-θ法則將動力平衡方程應(yīng)用于更后一點的時刻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)導(dǎo)出的{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時刻的加速度可按下式內(nèi)插求得{˙u(t+θΔt)}=(1-1θ){u(t)}+1θ{¨u(t+θΔt)}.(34)4.3計算誤差wilson-θ法的實質(zhì)是線性加速度法推廣。當(dāng)θ≥1.37時,此法是無條件穩(wěn)定的,但隨著θ的增大,計算誤差也增大,所以通常取θ=1.4。在地震作用下,對于一般阻尼比5%的鋼筋混凝土結(jié)構(gòu),時間步長Δt≤0.04T(T為地震波的卓越周期)可以取得較好的結(jié)果。5分散法的中心5.1初始點函數(shù)值的比較中心差分法的基本思想是:在計算函數(shù)的中心點差分,并與初始點函數(shù)值進行比較,若兩者之差足夠小則結(jié)束,則以中心點差分函數(shù)值作為結(jié)果;否則,步長減半,并將中心點差分函數(shù)值送給初始點函數(shù)值,繼續(xù)迭代,直到滿足誤差要求為止。5.23unt3+3.將位移增量函數(shù)按Taylor級數(shù)展開得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.速度和加速度也有同樣的關(guān)系,由式(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個時間步)的位移增量、速度和加速度滿足該時刻動態(tài)顯式有限元方法中的有限元微分方程,則[ˉΜ]¨un=Fextn-Fintn-Fen.(42)式中:M為對角矩陣,Fextn、Fintn和Fen是t時刻節(jié)