hello~

    在談數值阻尼之前，咱們先聊聊直接積分法的穩定性問題。

    直接積分法是不對運動方程進行任何變換，直接運動方程進行積分求解。本質上講，直接積分法是基于兩個基本概念：

    1.是將在求解域0≤t≤T內任何時刻t都應滿足運動方程的要求，代之以僅在一定條件下近似地滿足運動方程，例如僅在相隔△t的離散時間點0，△t，2△t，……上滿足運動方程。

    2.是在一定數目的△t區域內，假設位移D、速度V和加速度A的近似函數形式。

即在假定已知時刻t=0的位移D、速度V，并將時間求解域 0 ~ T等分成n個時間段0，△t，2△t，…，n△t（△t=T/n），且已求得時刻0，△t，2△t，…和t的解，現在要求t+△t時刻的解。建立了求解t+△t時刻解的算法后，即可類似地求得所有離散時刻的解，因此該方法也稱逐步積分法。

    詳情可以看下歷史文章【JY】結構動力學初步-單質點結構的瞬態動力學分析。

    逐步積分法是一種近似的求解方法，每一步積分計算都會帶來誤差。解的誤差的來源主要有兩大類：

    第一類是在計算速度和加速度的近似表達式中略去了高階小量，稱為截斷誤差，這類誤差是由方法本身決定的，一般可以作出估計，它隨步長△t增大而增大。

    第二類是由于計算機的字長總是有限位，超過其位數的數必須四舍五入，這類誤差稱為舍入誤差。盡管在每一步中舍入誤差并不大，但這類誤差在求解過程中可能會被不斷放大，以致使計算結果完全失真，這就是直接積分方法的穩定性問題。

    如果無論△t取多大，給定任意初始條件，積分結果都不會無界增大，則稱此方法是無條件穩定的。

    反之，如果△t必須小于某個臨界值時，積分結果才不會無界增大；則稱此方法是條件穩定的。

    在建立直接積分法格式時，我們假定已知時刻0，△t，2△t，…，t-△t和t的解，要求時刻t+△t的解，因此直接積分法利用時刻0，△t，2△t，…，t-△t和t的解遞推求解t+△t時刻的解，即直接積分法的格式可以寫成遞推的形式

    其中 Xt+Δt 和 Xt 為含有位移、速度等結果的向量，矩陣D稱為遞推矩陣或放大矩陣。

    上式右端的第二項與載荷有關，在進行穩定性分析時可令L=0，此時給定初始條件X0后，時刻 t+△t = (n+1)△t 的解 Xt+Δt 可以寫為

    對矩陣D進行譜分解可得

    假定矩陣未有重根(重根情況自行討論)，則可得到

    那么下式即為該遞推矩陣的譜半徑：

    因此可以看出，若是譜半徑≤1時，矩陣是有界的。綜上所述，直接積分法的穩定性可歸結為滿足下式即為穩定：

    如果滿足以上穩定準則，當n→∞時，矩陣D^(n+1)有界。因此如果譜半徑ρ(D)<1，當n→∞時，矩陣D^(n+1)→0，說明該積分方法存在數值阻尼（也稱為人工阻尼）。ρ(D)越小，當n→∞時矩陣D^(n+1)趨于零的速度越快，表明該積分方法的數值阻尼越大。

    因此譜半徑除了度量算法的穩定性外，也度量了算法的數值阻尼。

    對于一般結構動力學問題，系統的響應主要受低階振型控制，高階振型的貢獻很小。

    另外，由于受離散化的影響，有限元法或有限差分法對系統的高階振型的近似程度很差。如果在直接積分中不能有效地濾除這些虛假的高階分量，將會降低結果的精度。

    舉個例子，在【JY】結構動力學初步-單質點結構的瞬態動力學分析，這篇文章中，大型有限元軟件Ansys和Opensees的計算結果出現了較多的虛假高階分量，這篇文章做一個補充說明。

    可以看出Ansys和OpenSees等通用有限元軟件通過計算時，對于絕對加速度計算時會產生虛假高階高頻分量（相對加速度影響不大，因為絕對加速度=相對加速度+地面加速度，會導致誤差顯著化。）

    因此一個好的直接積分方法應在高頻段具有一定的可控的數值阻尼，以有效地濾除虛假的高頻振型對系統響應的影響，同時在低頻段的數值阻尼應盡可能小，以保證結果的精度。除了能濾除虛假的高階振型的影響外，數值阻尼還有助于非線性問題迭代求解的收斂性，并且也有助于接觸等具有約束的問題的求解。可見，在低頻段譜半徑 ρ(D) 應盡可能接近于1，而在高頻段 ρ(D) 應逐步光滑地減小，并趨于某個給定值。

    Wood建議(Wood WL. Practical time-stepping schemes. Oxford: Clarendon,1990)，當△t/ T(結構周期)趨于無窮大時，ρ(D)趨于0.5~0.8較為合適。

（采用Hilber-Hughes-Taylor方法，添加數值阻尼）

    注意：結構在進行動力計算時，數值阻尼雖然能有效濾除虛假的高頻振型，但同時會增加結構計算誤差，使得結構輸入輸出總能量不同，多出數值阻尼產生的計算誤差能量。因此，數值阻尼不宜過大，且要施加可控有效的數值阻尼（特別是單自由度或者等效單自由度結構，如：隔震結構），也可為抑制高頻段的振動，增加剛度阻尼(詳見：【JY】結構瑞利阻尼與經濟訂貨模型)，該方法不會使得結構輸入輸出總能量不同，即結構能量是平衡的，但是計算手段視模型而定，不一定可行明顯有效。

   【中心差分法】由于中心差分法所需要的時間步長比較短，實質上會讓該算法的譜半徑的模長等于1，也就是該算法并不能調整數值阻尼。

    順便提一句關于中心差分法的使用關鍵問題：對于n自由度系統，利用直接積分法對其運動方程進行積分，實質上等效于采用相同的時間步長 △t 同時對 n 個解耦后的微分方程進行積分。此時△t的選擇應保證對每個微分方程的積分都是穩定的。

    由于中心差分法的時間步長 △t 必須小于臨界步長△tcr，因此它是條件穩定的。在用有限元法進行動力分析時，系統自由度數很高，結構系統的有效最小周期比較小，因此△t必須選得很小。另外在中心差分法中要求質量矩陣的全部對角元素都大于零，如果有對角元素等于零或接近于零，△tcr也將等于零或接近于零，中心差分法不能使用。

    因此在用有限元法進行動力分析時，如使用中心差分法進行積分，則不宜使個別單元的尺寸比其他單元的尺寸小很多，否則會減小臨界時間步長，因而大幅度增加計算量。

   【Newmark-β法】Newmark-β法在適當的計算參數β、γ的選取下，是可以調整數值阻尼大小，甚至不要數值阻尼(如β=0，γ=1/2，則退化為中心差分法)。

    Newmark-β法使用的關鍵問題：有圖在導出Newmark積分格式時，使用的是t+△t時刻的運動方程，因此Newmark法是隱式方法。總所周知，Newmark法無條件穩定應滿足：

    如果不滿足上述條件，要得到穩定的解，時間步長△t必須小于臨界時間步長。

Newmark-β法的譜半徑為：

    當取為β=0，γ=1/2，此時Newmark法無數值阻尼，并具有二階精度。只有當 γ>1/2 時，Newmark法才具有數值阻尼，其數值阻尼隨著△t的增大而光滑地增大。但是，通過比較Newmark法的放大矩陣與系統的精確放大矩陣可知，此時Newmark法只有一階精度(截斷誤差)，若是 Δt 較長的情況，即便計算也能趨于穩定，但結果誤差也較大。為了使數值阻尼在高頻段最大以有效地濾除系統虛假的高頻響應，應使β取最小值，從而使ρ取最小值。

    【Houbolt法】Houbolt法是隱式計算方法，它的數值阻尼隨 △t/T 的增大而很快增大，因此解的震蕩衰減很快。當△t→∞時，Houbolt法的譜半徑趨于0。

    Houbolt法的使用關鍵問題，在Houbolt法中，不論△t取何值，在計算中Houbolt法的舍入誤差都不會越來越大， Houbolt法是無條件穩定的。當M=C=0時該方法的遞推方程可化為靜力方程，因此Houbolt法可用來求解載荷隨時間變化時的靜態解。但是在中心差分法中當M=C=0時，有效質量陣M=0，因此中心差分法不能用來求解靜態解，好在Houbolt法和中心差分法都是二階精度。

【HHT-α法】這個方法像Newmark和所有隱式方案一樣，此方法的無條件穩定性適用于線性問題，HHT法的數值阻尼大小由α決定。

    當 α = 0.0 時對應于Newmark方法，當滿足以下式子時，HHT法對高頻振動抑制最大。

    其中α可以取0~0.33，且越大時，數值阻尼越大。

（注意：此處表達和OpenSees的α值表達不同，根據OpenSees文檔中，α=1.0 對應Newmark法，和此處表達的內容是一樣的。且對應OpenSees的HHT法為α建議取值再在0.67和1.0之間，且較小的α數值阻尼越大）

    但是HHT方法中，即便是 Δt 較長的情況，計算也能趨于穩定，結果誤差也不會很大。

建源工程俠，下期再見！

~路過別錯過~