結構動力學：理論及其在地震工程中的應用.docx

5章 動力反應的數值計算 如果激勵作用力或地面加速度是隨時間任意變化的，或者體系是非線性的，那么對單自由度體系的運動方程進行解析求解通常是不可能的。這類問題可以通過數值時間步進法對微分方程進行積分來處理。在應用力學廣闊的學科領域中，有關各種類型微分方程數值求解方法的文獻（包括幾部著作中的主要章節(jié)）浩如煙海，這些文獻包括這些方法的數學進展以及它們的精度、收斂性、穩(wěn)定性和計算機實現等問題。 然而，本章僅對在單自由度體系動力反應分析中特別有用的很少幾種方法進行簡要介紹，這些介紹僅提供這些方法的基本概念和計算算法。盡管這些對許多實際問題和應用研究已經足夠了，但是讀者應該明白，有關這個主題存在大量的知識。 5.1 時間步進法 對于一個非彈性體系，欲采用數值求解的運動方程為 或者 (5.1.1)初始條件 假定體系具有線性粘滯阻尼，不過，也可以考慮其他形式的阻尼（包括非線性阻尼），后面會明顯看到這一點。然而由于缺乏阻尼信息因此很少這樣做，特別是在大振幅運動時。作用力由一系列離散值給出: , 到N。時間間隔 （5.1.2）圖5.1.1 時間步進法的記號通常取為常數，盡管這不是必需的。在離散時刻（表示為時刻）確定反應，單自由度體系的位移、速度和加速度分別為、和。假定這些值是已知的，它們在時刻滿足方程 （5.1.3）式中，是時刻的抗力，對于線彈性體系，但是如果體系是非彈性的，那么它會依賴于時刻以前的位移時程和速度。將要介紹的數值方法將使我們能夠確定+1時刻滿足方程(5.1.1)的反應、和，即在+1時刻 （5.1.4）對于=0，1，2，3，連續(xù)使用時間步進法，即可給出=0，l，2，3，所有瞬時所需的反應。已知的初始條件和提供了起動該方法的必要信息。 從時刻到+1時刻的步進一般不是精確的方法，許多在數值上可以實現的近似方法是可能的。對于數值方法，有三個重要的要求：(1)收斂性一隨著時間步長的減少，數值解應逼近精確解；(2)穩(wěn)定性一在存在數值舍入誤差的情況下，數值解應是穩(wěn)定的；(3)精度一數值方法應提供與精確解足夠接近的結果。這些重要的問題在本書中均作簡要的討論，全面的論述可在著重微分方程數值解法的書中找到。 本章介紹三種類型的時間步進法：(1)基于激勵函數插值的方法；(2)基于速度和加速度有限差分表達的方法；(3)基于假設加速度變化的方法。前兩類中各只介紹一種方法，第三類中介紹兩種方法。 5.2 基于激勵插值的方法對于線性體系，通過在每個時間間隔里對激勵進行插值，并利用第4章的方法進行精確求解，能推導出一種非常有效的數值方法。如果時間間隔較短，則線性插值是令人滿意的。圖5. 2.1所示的時間間隔，激勵函數為 （5.2.1a）其中 （5.2.1b）時間變量從0到變化。為數學上簡單起見，我們首先考慮無阻尼體系，后面再將該方法擴展到有阻尼體系。待求解的方程為 （5.2.2）在時間間隔內，反應為三部分之和：(1) =0時刻的初位移和初速度引起的自由振動；(2)零初始條件下對階躍力的反應；(3)零初始條件下對斜坡力的反應。對這三種情況分別采用來自2.1、4.3和4.4中已有的解答，得 （5.2.3a） （5.2.3b）計算時的這些等式，得+1時刻的位移和速度： (5.2.4a) (5.2.4b)將式（5.2.1b）代入后，可將這些等式重寫為如下的遞推公式： (5.2.5a) (5.2.5b)對于欠臨界阻尼體系（即），重復上面的推導，表明式(5.2.5)也適用于有阻尼體系，系數A，B，,的表達式由表5.2.1給出，這些系數取決于體系的參數、k和以及時間間隔。 表5.2.1 遞推公式中的系數（）因為遞推公式是從運動方程的精確解推導出的，因此對時問步長大小的唯一限制條件是，允許它對于激勵函數有一個接近的逼近，并以較密的時間間隔提供反應結果，以使反應峰值不會被漏掉。這種數值方法對于激勵由緊密的時間間隔定義的情況（例如對于地震地面加速度的情況）特別有用，從而使得線性捕值即可得到較完美的結果。如果時間步長是常數，則系數A，B，僅需計算一次。這種數值方法所要求的運動方程精確解答儀對線性體系是可行的。如上所述，這種方法用于單白由度體系較便利，但是對于多自由度體系則是不切實際的，除非它們的反應由振型反應的疊加（第12章和第13章）來獲得。5.3 中心差分法這種方法是基于對位移時間導數（即速度和加速度）的有限差分近似進行的。步長，則時刻的速度和加速度的中心差分表達式為 (5.3.1)將速度和加速度的這些近似表達式代人方程(5.1.3)中，對線彈性體系，得 (5.3.2)在這個方程中，和假定是已知的（來自于前面時間步內方法的執(zhí)行）。將這些已知量移到右側，導得 (5.3.3)或寫成 (5.3.4)其中 （5.3.5） （5.3.6）則未知的由下式給出 （5.3.7）+l時刻的解答是根據時刻的平衡條件即方程(5.1.3)確定的，而不是以時刻+1的平衡條件式(5.1.4)確定的，這種方法稱為顯式方法。觀察（5.3.6），為了計算，需要已知的位移因而，為了確定，需要和。特定的初始位移是已知的，為了確定，我們將式（5.3.1）專門用于=0的情況，得 （5.3.8）從第一個式子解出,然后代入第二個式子,給出 （5.3.9）其中初位移和初速度是已知的,由0時刻()的運動方程 可以得到0時刻的加速度為: （5.3.10）表5.3.1總結了可在計算機上執(zhí)行的上述方法。 表5.3.1 中心差分法如果時間步長選取得不夠短，那么由于數字舍入誤差的存在，中心差分法將會“放大”，而給出無意義的結果。為了穩(wěn)定性，特別要求 （5.3.11）對于單自由體系，上式永遠不會是一個約束，因為為了獲得準確的結果，選擇的時間步長將是非常小的。為了充分地定義反應，通常選擇；在大多地震反應分析中甚至選擇更短時間步長，為了準確地定義地面加速度，通常選取 =0.01到0.02秒。5.4 Newmark法5.4.1 基本方法1959年，NM. Newmark發(fā)展了一類時間步進法，它們基于下面的公式： （5.4.1a） （5.4.1b）參數和定義了時間步內加速度的變化，并決定方法的穩(wěn)定性與精度特征。對于為1/2和1/6l/4的典型選擇，從包括精度的所有觀點來看都是令人滿意的。這兩個等式與時間步結束時的平衡方程(5.1.4)結合，提供了從時刻已知的、和。計算時刻的、和的基礎。執(zhí)行這些計算需要迭代，因為未知的出現在式(5.4.1)的右側。然而，對于線性體系，修正Newmark法的原始公式可以允許使用式(5.4.1)和(5.1.4)求解時不迭代。在闡述這些修正之前，我們先來論證Newmark法的兩種特殊情況，即眾所周知的平均加速度法和線性加速度法。5.4.2 特殊情況對于這兩種情況，表5. 4.1總結了+1時刻的反應、和與時刻相應量之間的關系。式(5.4.2)描述了加速度在步長內變化是常數（等于平均加速度）或線性的假定。對進行積分，給出時間步長內速度的變化，即式(5.4.3)；將代入，得+1時刻的速度，即式(5.4.4)。對進行積分，給出時間步長內的位移的變化，即式(5.4.5)；將代人，得+l時刻的位移，即式(5.4.6)。將式(5.4.4)和式(5.4.6)與式(5.4.1)比較，可見Newmark方程在=1/2和=1/4時與平均常加速度假定推導的那些方程是相同的；= 1/2和=l/6時的方程則與加速度線性變化的假定相符。 表5.4.1 平均加速度法和線性加速度法5.4.3 非迭代表達式我們現在返回到式(5.4.1)，為避免迭代，對其重新進行列式，并使用增量 (5.4.7) (5.4.8)盡管對于線性體系的分析增量形式是不必要的，但是引入它是因為這種形式可以方便地擴展到非線性體系。式(5.4.1)可以重新寫為 (5.4.9)求解第二個等式，可得 (5.4.10)將式(5.4.10)代入式(5.4.9a)中，得 (5.4.11)下面，對具有和的線性體系，由方程(5.1.4)減去(5.1.3)，得到增量運動方程： (5.4.12)將式(5.4.10)和式(5.4.11)代人方程(5.4.12)中，得 (5.4.13)式中 (5.4.14) (5.4.15)由體系的特性,算法參數和,以及時間步長開始的和可知和,則增量位移由下式計算 (5.4.16)一旦求出,則和就可以根據式(5.4.1 1)和式(5.4.10)分別計算出來，從而、和可從式(5.4.7)計算出來。加速度也能從+l時刻的運動方程確定： (5.4.17)這要好于用式(5.4.10)和式(5.4.7)確定的值。式(5.4.17)開始計算需要得到 見式(5.3.10)。在Newmark法中，+l時刻的解由式(5.4.12)確定，這是與使用+1時刻的平衡條件方程(5.1.4)等價的。這種方法稱為隱式法。表5.4.2總結了可在計算機上執(zhí)行的用Newmark法進行時間步進求解的過程。如果時問步長滿足下式，則Newmark法是穩(wěn)定的: (5.4.18)對于 = 1/2，=1/4，這個條件成為 (5.4.19a)這意味著，平均加速度法對任何都是穩(wěn)定的，無論它有多大；然而,如5.3結束時所討論的，僅當足夠小時才是準確的。對于=12和1/6，式(5，4.1 8)表明，如果下式成立，則線性加速度法是穩(wěn)定的 (5.4.19b)然而，像中心差分法的情況一樣，這個條件任單自山度體系分析中沒有什么意義，因為為了獲得激勵和反應的正確表示，必須使用比短得多的時間步長。 表5.4.2 Newmark法:線性體系5.5 穩(wěn)定性與計算誤差5.5.1 穩(wěn)定性如果時間步長比穩(wěn)定性界限短，則稱導致有界解答的數值方法為條件穩(wěn)定方法。不管時間步長的長短，均導致有界結果的方法稱為無條件穩(wěn)定方法。平均加速度法是無條件穩(wěn)定的，線性加速度法在時是穩(wěn)定的，中心差分法在時是穩(wěn)定的。顯然，后兩種方法是條件穩(wěn)定方法。在單自由度體系分析中，穩(wěn)定性準則不是限制性的（即，它們并不規(guī)定時間步長的選擇）。這是因為，為保證數值結果有適當的精度。必須比穩(wěn)定性界限小得多（比方說0.1或更?。?，然而，數值方法的穩(wěn)定性在多自由度體系的分析中則是重要的，經常需要使用無條件穩(wěn)定方法（第1 5章）。5.5.2 計算誤差誤差是運動方程任何數值解答中所固有的，我們不從數學的觀點來討論誤差分析。為了對誤差的本質有一個感性認識，我們先來分析數值解的兩個重要特征，然后介紹一種處理誤差的簡單而有效的方法。考慮自由振動問題 和 理論解為 (5.5.1)采用四種數值方法：中心差分法、平均加速度法、線性加速度法和Wilson法對這個問題進行求解。其申，最后一種方法將在第15章中介紹。取時所得的數值結果與理論結果進行比較，如圖5.5.1所示。這些比較說明有些數值方法可能會預測出位移幅值會隨時間衰減，盡管體系是無阻尼的；并且會預測出固有周期被延長或縮短。 圖5.5.1 采用四種數值方法()和理論解的真由振動解答圖5.5.2示出了四種數值方法中振幅衰減AD和周期延長PE作為函數的圖形，AD和PE分別在圖b和圖c中定義，不過，沒有介紹導出這些數據的數學分析。有三種方法顯示位移幅值無衰減。但是，Wilson法含有幅值的衰減，意味著該方法在體系中引入了數值阻尼，等效粘滯阻尼比示于圖a中。注意，中心差分法的周期誤差在接近方法的穩(wěn)定性界限時迅速增加。中心差分法引入了最大的周期誤差，在這個意義上，它被認為是精度最差的方法。對于比它的穩(wěn)定性界限小的情況，線性加速度法給出最小的周期延長。聯合這個性質和幅值無衰減的特點，使得線性加速度法是單自由度體系已介紹的方法當中最合適的方法。然而，對于多自由度體系，因為穩(wěn)定性需求，我們將得到不同的結論（第15章）。除體系的固有振動周期以外，步長的選擇還與動力激勵隨時間的變化有關。圖5.5.2表明，取會給出相當精確的結果。為了保持激勵函數的最小失真，時間步長也將是足夠短。為了用數字描述地震發(fā)生時所記錄的高度不規(guī)則的地震地面加速度，通常需要非常精細的時間步長，典型的時間步長為 =0. 02秒，對于計算結構反應所選擇的時間步長不應該比這更長。一種雖然不復雜，但非常有用的選擇時間步長的方法是：先用一個看起來合理的時間步長求解問題，然后以稍小的時間步長重復求解，對比結果，持續(xù)這個過程直到連續(xù)的兩個解足夠接近。前述關于穩(wěn)定性和精度的討論嚴格地說適用于線性體系。對于這些問題是怎樣影響非線性反應分析的，讀者應參考其他資料。5.6 非線性反應分析：中心差分法超過線彈性范圍的體泵的動力分析一般不適合用解析方法求解，即使激勵隨時間的變化可以用一個簡單函數來描述。因此，數值方法是非線性體系分析的基本方法?？扇菀椎貙χ行牟罘址ㄟM行修改，以求解時刻的非線性運動方程(5.1.3)。將速度和加速度的中心差分近似式(5.3.1)代人，得到式(5.3.2)，其中由代替，式(5.3.2)可重新寫成求+1時刻反應的表達式 (5.6.1)其中， (5.6.2) (5.6.3)將這些等式與線性體系的那些等式進行比較，可見惟一的差別是對的定義。作此修正后，表5.3.1也適用于非線性體系?？沽雌饋硎秋@式的，因為它僅取決于時刻的反應，而不依賴于未知的+1時刻的反應。因此它容易計算，使得中心差分法對非線性體系來說或許是最簡單的方法。盡管在這方面吸引人，但是由于有更有效的方法可用，因此在實際應用或者研究應用中這種方法并不流行。5.7 非線性反應分析：Newmark法本節(jié)將5.4所述的線性體系的Newmark法擴展到非線性體系。盡管不像中心差分法那樣簡單，但是由于它有非常高的精度，因此這種方法也許足地震反應分析中最流行的方法。式(5.1.3)和式(5.1.4)之差給m增量平衡方程： (5.7.1)增量抗力為 (5.7.2)其中，圖5.7.1中所示的割線剛度不能確定，因為u，是未知的。如果我們做這樣的假定，在一個微小的時間間隔 內，割線剛度可用切線剛度代替，如圖5. 7.1所示，那么，式(5,7.2)可近似地表示為 圖5.7.1 (5.7.3)去掉式(5.7.3)中的下標T，并將該式代入方程(5.7.1)中，得 (5.7.4)這個方程與線性體系的相應方程(5.4.1 2)之間的相似之處啟示我們，前面對線性體系介紹的Newmark法的非迭代公式也可以用于非線性反應分析。所有需要做的是將式(5.4.14)中的用每個時問步開始時求出的切線剛度代替。這個變化意味著表5. 4.2中的步驟1.3應該移到步驟2.1的后面。對于非線性體系，步驟2.5和式(5.4.17)會給出最的不同值，后者給出的值更好一些，因為它滿足+1時刻的平衡方程。這個方法用不變的時間步長會導致不能接受的錯誤結果。重大的誤差起源于兩個原因：(1)用切線剮度代替割線剛度；(2)常數時間步長推遲了力一變形關系中轉折點的發(fā)現。首先，我們考慮誤差的第二個來源，用圖5.7. 2a所示的力變形關系來說明。假設時間步開始時的時刻處位移為，速度是正的（即位移是漸增的），對應圖中的a點。對時間步應用前面描述的數值方法，求得+l時刻的位移和速度，這是圖中的b點。 圖5.7.2如果是負的，那么在時間步內的某點，速度為零并將改變符號，位移開始減少。在數值方法中，如果我們不想麻煩找到點，而是在b點開始下一個時間步繼續(xù)計算，并使用與力一變形圖的卸載分支相關的切線剛度，那么這個方法將在下一個時同步結束時定位于c點，位移為，速度為負的。另一方面，如果能夠確定與點相關的時間瞬時（當速度實際成為零時），那么下一個時間步的計算將從體系在的狀態(tài)開始，給出時間步結束時的位移和速度，記為。不定位，具有超越到b，以及不跟隨力一變形圖上精確路徑的影響。這些與精確路線的偏離將發(fā)生在每一個速度反向處，從而導致數值結果中的誤差。類似的問題也出現在力一變形關系（例如彈塑性體系）的其他尖角處。通過精確定位可以避免這些誤差，這可以南折回去以更小的時間步長(比如說)在到時間間隔內的積分來達到。另外，還可以使用一個迭代過程，用比整個時間步長還小的步長，從時間點重新開始積分，對積分步長進行連續(xù)調整，以使這樣調節(jié)的步長在結束時的速度接近于零?，F在，我們回到誤差的第一種原因上來，它與使用切線剛度代替未知的割線剛度有關，如圖5.7.2b中的力一變形關系所示。在時間步開始時刻i的位移為圖示的口點，在血點使用切線剛度，從時刻到時刻+1的數值積分導出位移，標識為點b。如果我們能沿若正確曲線，那么結果會是點對應的位移。這個偏差經過一系列時間步的累積，會引入非常大的誤差。這些誤差可用迭代方法使其最小化。在Newmark法的每一步中要求解的關鍵方程是式(5.4.13)，對于非線性體系，經改造成為： (5.7.5)式中，由式(5-4.15)給出， （5.7.6）為符號的方便起見，我們將的下標去掉，用T代替，以強調這是切線剛度；還有和的下標也去掉。式(5.7.5)和(5.7.6)成為 (5.7.7) （5.7.8）圖5.7.3a給出了式(5.7.7)的示意圖。關系是非線性的，因為切線剛度取決于位移，因此斜率不是常數。在非線性體系的靜力分析中，=，的非線性特性與的相同。在動力分析中，中的質量和阻尼項的存在減少了非線性特性，因為常數項對于的典型值通常比大許多。 圖5.7.3非線性體系一個時間步內的迭代：(a)修正的Newton - Raphson迭代；(b)Newton - Raphson迭代接下來參考圖5.7. 3a敘述迭代方法。第一個