時間積分方法的研究進展與挑戰

邢譽峰 季奕 張慧敏

(1. 北京航空航天大學 固體力學研究所, 北京 100083;2. 北京理工大學 宇航學院, 北京 100081; 3. 北京宇航系統工程研究所, 北京 100076)

時間積分方法[1](time integration method,TIM)是求解動力學常微分方程最有力的數值方法之一,也是目前計算機輔助工程(computer aided engineering, CAE)軟件中的主流時域求解器。 時間積分方法的思想[1]是:將待求時間域離散成為一系列時間區間或時間單元,并建立每個時間區間內位移、速度和加速度隨時間的變化規律,進而由前一時刻的已知狀態變量遞推計算當前時刻的未知狀態變量。 由此可見,時間積分方法有2 個特點[1]:①動力學平衡方程僅在離散時間點或廣義時間點上得到滿足,而不是在任意時刻都成立;②需要構造狀態變量在時間區間或時間步內的變化規律。 這些變化規律決定一個方法的精度、耗散、效率和穩定性。 這與空間有限元方法類似,也需要設計單元的容許位移函數。 也就是位移隨著空間坐標的變化規律,并且代數平衡方程也僅在結點上得到滿足。 有限單元的容許位移函數決定了精度、效率和適用性。

級數展開法、Runge-Kutta 法和Newmark 方法是經典的時間積分方法。 基于級數展開思想直接構造的格式皆屬于條件穩定的顯式時間積分方法,代表性工作有Taylor 級數法和Lie 級數法。Taylor 級數法直接對時間坐標求導,而Lie 級數法則對狀態變量求導。 求導方式的不同使得Lie 級數法可以解決Taylor 級數法中高階求導困難和計算量大的問題。 Runge-Kutta 法是基于Taylor 級數展開思想于19 世紀末、20 世紀初提出的高精度方法。 Runge-Kutta 法[2]的求解對象是形如dz/dt=f(z,t)的一階常微分方程(ordinary differential equation, ODE),其利用函數f在若干時間離散點上信息的線性組合來表示f的導數,再利用Taylor級數法確定其中的系數。 與級數展開法和Runge-Kutta 法不同,雖然Newmark 方法設計思想也是基于級數展開,但其求解對象為二階ODE。Newmark 方法包括許多著名的辛格式,如梯形法則(trapezoidal rule, TR)、中心差分方法(central difference method, CDM)、Fox-Goodwin 方法等。

隨著科學技術的發展,人們對時間積分方法的精度、效率、耗散和穩定性的要求越來越高,經典方法已經難以滿足需求,這推動了時間積分方法在過去幾十年中的持續發展。 激發人們設計高性能時間積分方法的主要問題可以概括如下:

1) 具有數值耗散的Newmark 方法[1]只有一階精度。

2) 收斂階數大于2 時,線性多步法只能是條件穩定或不穩定[3]。

3) 對線性系統具有無條件穩定性的二階Newmark 辛幾何方法,如TR,對于簡單的非線性問題可能發散[4-5]。

4) 如何同時提高時間積分方法的精度、效率、耗散性能和穩定性。

5) 對非線性動力系統而言,時間積分方法的失穩機理,二階非線性系統的無條件穩定性方法的設計原則。

針對問題1,通過對動力學平衡方程進行加權處理[6-10],或在差分格式中引入輔助狀態變量[11-12],學者們發展了參數類耗散方法,該類耗散方法具有二階精度。 為了解決問題2,學者們利用微分求積法[13-15]和加權殘量法[16-18]等構造出了具有無條件穩定性的高階時間積分方法。 針對問題3,借助能量有界準則[19],學者們發展了可穩定計算非線性動力響應的保能量方法[20-30]。為了全面提高時間積分方法的數值性能,通過使用更多時刻的已知狀態變量,或不同分步采用不同類型的時間積分方法,學者們提出了線性多步法[31-38]和復合方法[39-50],這2 類方法可有效解決問題4。 問題5 更具有挑戰性。 為了提高穩定性,學者們提出了結構相關時間積分方法[51-54],在該類方法中,算法參數與系統質量、阻尼和剛度矩陣相關,數值結果表明,該類方法能夠有效提高數值積分方法處理軟剛度系統時的穩定性。 此外,學者們還發現可利用幅值有界準則[55]或時變參數[31]解釋在求解非線性問題時已有時間積分方法的失穩機理,并據此幅值有界準則設計了對一階和二階ODE 都是無條件穩定的BN 穩定型方法[56]。

本文著重介紹了針對上述問題提出的先進時間積分方法的發展狀況,并對時間積分方法的未來發展趨勢進行展望。

1 參數方法

為了提高Newmark 方法耗散格式的收斂階次,學者們基于Newmark 方法發展了具有數值阻尼的二階參數方法。 根據構造方式不同,參數方法可以分成2 類:一類是加權動力學平衡方程,另一類是在差分格式中引入輔助變量。

第一類參數方法的代表性工作有廣義-α方法[6-8]、HHT-α方法[9]和WBZ-α方法[10]。 這類方法沿用了Newmark 方法的差分格式,即

式中:M、C、K和R分別為質量矩陣、阻尼矩陣、剛度矩陣和外載荷列向量。

由式(2)可以看出,第一類參數方法不滿足離散時間節點的平衡方程。 廣義-α方法的參數可統一寫成高頻極限處譜半徑ρ∞的函數,即

上述參數可保證廣義-α方法具有二階精度,并且高頻段的數值耗散程度可由ρ∞精確調控。當α=δ=0 時,廣義-α方法退化為Newmark 方法;當α=0 時,廣義-α方法退化為HHT-α方法;當δ=0 時,廣義-α方法退化為WBZ-α方法。 以廣義-α方法為代表的第一類參數方法采用了加權思想,使得動力學平衡方程不能被嚴格滿足,這導致由該類參數方法計算的加速度僅具有一階精度。 為了改善該性能,張慧敏和邢譽峰[11]提出了第二類參數方法,其在時間節點處嚴格滿足動力學平衡方程,即

因此,該類方法被稱為三參數單步法(three parameters single-step method, TPSM)。 引入的輔助變量θ是一個中間變量,其本身不需要滿足任何狀態方程。 因此,在初始時刻,θ可以設置為0以避免額外的啟動程序。 采用相同的構造方式,邢譽峰和張慧敏[12]還發展了四參數保辛單步法(four parameters single-step method, FPTM)。

2 高階無條件穩定方法

通過減小時間步長來提高精度最容易實現,這相當于有限元方法中的h收斂性,但不可避免會導致計算量增加,而且會引入舍入誤差。 升高時間積分方法的收斂階次是另一種提高計算精度的有效方式,相當于有限元方法中的p收斂性。雖然線性多步法的穩定性受Dahlquist 定理[3]的約束,即當收斂階次大于2 時,算法是條件穩定或不穩定,不過學者們利用多級升階思想構造出了具有無條件穩定性的高階時間積分方法。

多級方法的定義是:在每個時間步內,在多個時刻或多次在同一時刻滿足動力學平衡方程的時間積分方法。 代表性工作有Runge-Kutta 法[2]、微分求積法[13-15]和加權殘量法[16-18]等。 多級方法的穩定性不受Dahlquist 定理約束,可以同時實現高階精度和無條件穩定性,隨著算法級數的增加,方法的收斂階次隨之提高,但伴隨著成倍增加的計算量。

基于多級思想建立的Runge-Kutta 法[2]因其精度高而被廣泛用于求解一階ODE。 若用其求解二階ODE,則需將二階ODE 轉化為一階ODE。Runge-Kutta 法按構造格式可以分為顯式、對角隱式和完全隱式3 類。 利用高斯積分,s級完全隱式Runge-Kutta 法可以達到2s階精度。 但是,如此改善精度導致計算成本大幅度增加,因此其實用價值不高。 綜合精度、效率和穩定性,對角隱式Runge-Kutta 法更受歡迎。 Runge-Kutta 法的求解格式為

由式(8)可以看出,對角隱式Runge-Kutta 法的系數矩陣A是一個對角元素不全為0 的下三角矩陣。 通過參數設計和選擇足夠的級數,隱式Runge-Kutta 法可以具有任意高階精度和無條件穩定性。

基于微分求積思想、加權殘量思想、最小二乘法等,Fung[13-14,16,18]構造出一些高階無條件穩定方法。 在Fung 的工作基礎上,一些性能更加優秀的高階方法被提出,如Huang 方法[57]和Kim 方法[58]。 Huang 和Fu[57]將不同格式的Fung 方法混合使用在一個時間步長內,從而有效提高了低頻精度。 在不損失高階精度和無條件穩定性的前提下,Kim 方法[58]優化了一個時間步內時間離散點的位置,使得最后一個時間配點與時間步的終點重合,從而簡化了Fung 方法的執行流程。 下面以Kim 方法為例,簡述該類高階方法的構造流程。 在Kim 方法中,位移、速度和加速度的表達式如下:

式中:ψj為第j個時間點的Lagrange 插值函數;uj、vj和aj分別為第j個時間點的位移、速度和加速度。

Lagrange 插值函數的定義為

式中:tk為第k個時間步的初始時刻;τs和Δt為第k個時間步內的配點。

在式(9)中,位移、速度和加速度是彼此獨立的。 Kim 方法采用加權殘量法建立速度-位移和加速度-速度之間的關系,即定義如下殘量:

結合t+τiΔt(i=1,2,…,n)時刻的動力學平衡方程,即可求解出所有未知的狀態變量。 基于微分求積思想、加權殘量思想、最小二乘法等建立的高階方法在有數值耗散時具有(2n- 1)階精度,在無數值耗散時是2n階精確的。

3 保能量方法

著名的TR(即歐拉中點辛差分格式)譜半徑特性表明其是一種無條件穩定方法,但將其應用于求解一些簡單的非線性問題時,如剛性單擺[4]和軟彈簧方程[5],發現若時間步長選取的不合適,TR 會發散。 雖然數值阻尼可以增強穩定性,但當把耗散方法用于非線性問題求解時,其結果仍可能不收斂。

為了穩定地計算非線性動力學響應,學者們基于能量有界準則[19]建立了保能量方法。 在該類方法中,當前時間步的機械能小于或等于上一時間步的機械能與外力功之和。 據此準則,Hughes 等[22]首先提出了約束能量型方法(constraint energy method, CEM)。 CEM 以TR 為 基礎,通過Lagrange 乘子法施加能量準則約束。 然而,該方法需要額外計算Lagrange 乘子和離散的能量值,在非線性迭代中還會遇到收斂困難問題。不同于CEM,Simo 和Wong[30]提出了另一種保能量方法,即能量-動量法(energy momentum method, EMM)。 EMM 是基于修正的中點法則建立的,在使用Green 應變的非線性系統中可以保守能量,但其難以推廣到一般非線性系統。

針對一般非線性系統, 已有 Gonazlez 方法[26]、Krenk 方法[25]等。 這些方法均基于中點法則,修正后的中點內力能夠滿足能量守恒的要求,但也都需要計算離散點的能量值。 為了解決這個問題,筆者直接從能量準則出發,提出了一種更為通用的保能量方法[21](energy-conserving method,ECM)。 該方法利用Gauss-Legendre 求積法則計算平均內力,數值結果表明,只要使用足夠的節點數,就可以精確地保守系統能量。 更重要的是,該方法僅需已知內力形式即可,不需要計算能量,并且適用于同時包含阻尼和剛度非線性的動力學系統。

下面以CEM 為例,簡述保能量方法的實現過程。 考慮如下非線性保守系統:

式(17)是所有保能量方法的設計準則。

CEM 是基于TR 設計的,對于方程(14),其遞推格式如下:

在CEM 中,定義一個泛函Γ(xt+Δt),其形式如下:

由式(21)可見,式(18)等價于泛函Γ(xt+Δt)的一階變分等于零。 為了滿足能量守恒方程(17),CEM 將Γ(xt+Δt)修正為

通過式(23)和式(24)可聯立求解出xt+Δt和λ,進而利用式(19)計算t+Δt時刻的速度和加速度。 從方程(24)可以看出,每一步CEM 都需要計算能量,與TR 相比計算量略有增加。

4 線性多步法和復合方法

為了綜合提高時間積分方法的數值性能,包括精度、效率、耗散和穩定性,可以采用線性多步法和復合方法。 線性多步法屬于單步時間積分方法中的一種,但為了計算當前時刻的響應xt+Δt,需要利用其他自啟動方法計算前幾個時刻的響應,以獲得足夠的啟動信息。 與線性多步法不同,復合方法是將一個時間步劃分成若干個分步,計算xt+Δt需要每個分步的信息。

4.1 線性多步法

工程界目前使用較多的線性多步法[2]有Adams 類方法、顯式Adams-Bashforth 方法、隱式Adams-Moulton 方法及預測校正類格式等。 這些方法雖然具有高階精度,但穩定性難以提高。 Dahlquist[3]曾證明,所有顯式方法和二階以上的線性多步法均無法實現無條件穩定。 因此,在實際應用中,顯式和隱式線性多步方法常作為預測和校正方法成對使用,并結合自適應步長技術,控制局部截斷誤差,以避免結果發散。 此外,向后差分格式,如Houbolt 方法[59]和Park 方法[1],也屬于經典的線性多步法。 它們一般具備二階精度、無條件穩定性及強烈的高頻耗散能力,可以快速濾掉高頻響應成分,因此更適用于剛性動力學系統的分析。

為進一步提高線性多步法的計算精度,學者們構造出了更加精確的線性三步法[33]和線性四步法[33]。 另外,為了解決線性多步法無法自啟動的問題,通過引入一些輔助狀態變量,學者們還提出了與一些線性多步法具有相同數值性能的單步法[32-33]。

下面以線性兩步法[33](linear two-step method, LTM)為例,介紹該類方法的求解流程。 LTM的遞推格式為

由式(26)可以看出,LTM 需利用其他自啟動方法計算t- Δt時刻的狀態變量,之后才能執行式(26)。

4.2 復合方法

復合方法在一個時間步的不同分步內使用不同的時間積分方法,以發揮不同方法各自的優勢,從而獲得優于各分步方法的數值性能。

Bank 等[60]在1985 年把TR 和后向差分公式(backward difference formula, BDF)組合使用,并將其用于求解一階ODE。 Bathe 和Baig[46]在2005 年將這種復合方法TR-BDF 用于求解二階ODE,并發展了復合方法概念。 隨后,Bathe 等進一步研究了TR-BDF 方法在非線性系統[61]、波傳播[62]等問題中的性能。 TR-BDF 方法優越的低頻精度和高頻耗散能力受到了學術界和工程界的廣泛關注。 在Bathe 方法基礎上,學者們提出了一系列性能更加優越的復合方法,如三分步TR-TRBDF[50]、三分步TR-BDF-Houbolt[48]、四分步Optimal-TR-TR-TR-BDF[42]等。 與 兩 分 步 Bathe 方法[46]相比,這些方法的低頻精度更高。

由TR 和BDF(或Houbolt 方法)組合的復合方法具有L 型穩定性。 為了能夠精確地調控高頻段的數值耗散量,基于標準Bathe 方法[46],Rezaiee-Pajand 和Sarafrazi[37]提 出 了β1/β2-Bathe 方法,該方法的第二分步采用兩點后向差分公式(backward interpolation formula, BIF),通過調整β1和β2的大小和比例,可以控制高頻段的數值耗散量。 之后,Noh 和Bathe[44]進一步優化了β1/β2-Bathe 方法,提出了ρ∞-Bathe 方法,僅使用一個參數ρ∞來光滑和準確地改變數值耗散的大小。 在具有可控數值耗散的前提下,為了進一步提高該類方法的低頻精度,季奕和邢譽峰[43]提出了一種優化三分步方法。

下面以Bathe 方法[46]為例介紹復合方法的求解思路。 在Bathe 方法中,把一個時間步[t,t+Δt]劃分為2 個分步[t,t+γΔt]和[t+γΔt,t+Δt],其中γΔt表示第1 個分步的大小。 TR 被用在第1 個分步中,即

式中:θ0、θ1和θ2為算法參數。

聯立求解t+γΔt和t+Δt時刻的動力學平衡方程,即可得到一個時間步內所有狀態變量。

另外,針對由TR 和BDF(或BIF)復合得到的時間積分方法,邢譽峰和季奕等[42]研究了如何選擇分步數、分步步長、差分點數、各分步采用的方法,并得到了一些重要的結論,為復合方法的構造提供指導,這些結論包括:

1) 當BDF(或BIF)采用的差分點數不變時,隨著分步數的增加,TR 在一個時間步中的占比增加,從而提高了低頻精度。 但是,這種方法不能明顯提高低頻精度。

2) 當分步數固定時,提高BDF 的差分點數(或BIF 中的插值點數),可以有效提高低頻精度。 實際上,隨著差分點數目增加,一個時間步內采用TR 的分步的占比略有降低。 因此,可以認為增加差分點數比增加分步數更有助于提高低頻精度。

3) 當分步數為常數時,只要采用相同的參數優化方式,就可以構造出性能相近的復合方法。如三分步優化Optimal-TR-BDF-BDF 方法中采用TR 的分步的占比大約是三分步優化Optimal-TRTR-BDF 方法中的一半,但二者的計算精度幾乎相同。

5 BN 穩定型方法

引入數值阻尼是改善時間積分方法在非線性系統中穩定性的一種直接且有效的方式。 與非耗散方法相比,一般情況下耗散方法在非線性系統中具有穩定性優勢,但其精度變低。 為了保證在計算非線性動力學響應時的穩定性,保能量方法是一種選擇。 但保能量方法具有額外的計算量,并且適用范圍有限。

下面簡單介紹ρ∞-TSSBN 方法的求解格式和算法參數。 該方法包含2 個分步[t,t+c1Δt] 和[t+c1Δt,t+c2Δt],其中,0 ＜c1＜c2＜1。 第1 個分步的求解格式為

式中:F為包含內力和外力的向量。

由式(29)和式(30)可求解出t+c1Δt和t+c2Δt兩個時刻的狀態變量。 時間步終點t+Δt的信息是由這2 個時刻的狀態變量加權得到,即

從式(29)可以看出,t0+c1Δt時刻狀態變量的求解不需要初始時刻t0的加速度信息,因此,該方法是一種真正的自啟動方法。

從式(31)可以看出,t+ Δt時刻不需要有效剛度矩陣分解和Newton 迭代運算。 因次,從計算量的角度來說,該方法屬于兩分步格式。

為了確保ρ∞-TSSBN 方法能夠穩定地處理非線性系統,即具有BN 穩定性,利用BN 穩定性理論的充分條件(即代數穩定性)設計了該算法參數。 當ρ∞=1 時,算法參數為

數值算例表明,在TR 和ρ∞-Bathe 方法失效的非線性問題中,ρ∞-TSSBN 方法可以給出穩定且精確的結果。

6 結束語

本文介紹了用于求解動力學常微分方程的時間積分方法,并著重介紹了近年來發展出的一些先進時間積分方法,包括參數方法、高階無條件穩定方法、保能量方法、線性多步法、復合方法和BN穩定型方法。 更詳細的資料讀者可以閱讀相關的綜述文章[63]和專著[64]。

關于時間積分方法的未來發展問題,可以從其求解的是線性問題和非線性問題角度進行討論。

對于線性動力系統,時間積分方法的設計已有成熟的理論基礎。 雖然已經存在高精度、快速計算方法[65-67],但如何更快、更精確地計算仍然是人們追求的目標。 對于線性單自由度動力系統,已經可以利用時間積分方法得到精確的結點位移[68-69],其具有精確的幅值和相位。 對于線性多自由度系統,雖然辛幾何算法[70]具有保幅值或保能量特性,或者說從根本上解決了幅值誤差累計問題,但相位誤差累計問題仍然存在。 幅值有界準則和能量準則都是針對幅值衰減和發散問題提出的,但建立一個針對相位誤差的約束準則是值得探索的。

對于復雜非線性動力系統,時間積分方法的精度和穩定性應該是更加本質的問題。 季奕等[56]基于BN 穩定性理論已經設計出兩分步二階無條件穩定時間積分方法,但發展更精確高效的BN 穩定型方法是一項有意義和挑戰的工作。此外,雖然季奕和邢譽峰[31]已經建立了時變參數譜分析理論,并利用此理論設計了對非線性剛度系統是無條件穩定的方法,但尚未將這種方法推廣用于非線性阻尼系統等。 非線性系統的高階時間積分方法的設計及多尺度快速計算問題等也都是值得關注的問題。

現代科學發展的特征是多學科交叉、多物理場耦合和多尺度分析和非線性特征,對于結構動力學系統時間積分方法的研究工作也要與現代科學發展的步伐協調一致。