時(shí)間積分方法的研究進(jìn)展與挑戰(zhàn)

邢譽(yù)峰 季奕 張慧敏

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

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

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

隨著科學(xué)技術(shù)的發(fā)展,人們對(duì)時(shí)間積分方法的精度、效率、耗散和穩(wěn)定性的要求越來(lái)越高,經(jīng)典方法已經(jīng)難以滿(mǎn)足需求,這推動(dòng)了時(shí)間積分方法在過(guò)去幾十年中的持續(xù)發(fā)展。 激發(fā)人們?cè)O(shè)計(jì)高性能時(shí)間積分方法的主要問(wèn)題可以概括如下:

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

2) 收斂階數(shù)大于2 時(shí),線(xiàn)性多步法只能是條件穩(wěn)定或不穩(wěn)定[3]。

3) 對(duì)線(xiàn)性系統(tǒng)具有無(wú)條件穩(wěn)定性的二階Newmark 辛幾何方法,如TR,對(duì)于簡(jiǎn)單的非線(xiàn)性問(wèn)題可能發(fā)散[4-5]。

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

5) 對(duì)非線(xiàn)性動(dòng)力系統(tǒng)而言,時(shí)間積分方法的失穩(wěn)機(jī)理,二階非線(xiàn)性系統(tǒng)的無(wú)條件穩(wěn)定性方法的設(shè)計(jì)原則。

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

本文著重介紹了針對(duì)上述問(wèn)題提出的先進(jìn)時(shí)間積分方法的發(fā)展?fàn)顩r,并對(duì)時(shí)間積分方法的未來(lái)發(fā)展趨勢(shì)進(jìn)行展望。

1 參數(shù)方法

為了提高Newmark 方法耗散格式的收斂階次,學(xué)者們基于Newmark 方法發(fā)展了具有數(shù)值阻尼的二階參數(shù)方法。 根據(jù)構(gòu)造方式不同,參數(shù)方法可以分成2 類(lèi):一類(lèi)是加權(quán)動(dòng)力學(xué)平衡方程,另一類(lèi)是在差分格式中引入輔助變量。

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

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

由式(2)可以看出,第一類(lèi)參數(shù)方法不滿(mǎn)足離散時(shí)間節(jié)點(diǎn)的平衡方程。 廣義-α方法的參數(shù)可統(tǒng)一寫(xiě)成高頻極限處譜半徑ρ∞的函數(shù),即

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

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

2 高階無(wú)條件穩(wěn)定方法

通過(guò)減小時(shí)間步長(zhǎng)來(lái)提高精度最容易實(shí)現(xiàn),這相當(dāng)于有限元方法中的h收斂性,但不可避免會(huì)導(dǎo)致計(jì)算量增加,而且會(huì)引入舍入誤差。 升高時(shí)間積分方法的收斂階次是另一種提高計(jì)算精度的有效方式,相當(dāng)于有限元方法中的p收斂性。雖然線(xiàn)性多步法的穩(wěn)定性受Dahlquist 定理[3]的約束,即當(dāng)收斂階次大于2 時(shí),算法是條件穩(wěn)定或不穩(wěn)定,不過(guò)學(xué)者們利用多級(jí)升階思想構(gòu)造出了具有無(wú)條件穩(wěn)定性的高階時(shí)間積分方法。

多級(jí)方法的定義是:在每個(gè)時(shí)間步內(nèi),在多個(gè)時(shí)刻或多次在同一時(shí)刻滿(mǎn)足動(dòng)力學(xué)平衡方程的時(shí)間積分方法。 代表性工作有Runge-Kutta 法[2]、微分求積法[13-15]和加權(quán)殘量法[16-18]等。 多級(jí)方法的穩(wěn)定性不受Dahlquist 定理約束,可以同時(shí)實(shí)現(xiàn)高階精度和無(wú)條件穩(wěn)定性,隨著算法級(jí)數(shù)的增加,方法的收斂階次隨之提高,但伴隨著成倍增加的計(jì)算量。

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

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

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

式中:ψj為第j個(gè)時(shí)間點(diǎn)的Lagrange 插值函數(shù);uj、vj和aj分別為第j個(gè)時(shí)間點(diǎn)的位移、速度和加速度。

Lagrange 插值函數(shù)的定義為

式中:tk為第k個(gè)時(shí)間步的初始時(shí)刻;τs和Δt為第k個(gè)時(shí)間步內(nèi)的配點(diǎn)。

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

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

3 保能量方法

著名的TR(即歐拉中點(diǎn)辛差分格式)譜半徑特性表明其是一種無(wú)條件穩(wěn)定方法,但將其應(yīng)用于求解一些簡(jiǎn)單的非線(xiàn)性問(wèn)題時(shí),如剛性單擺[4]和軟彈簧方程[5],發(fā)現(xiàn)若時(shí)間步長(zhǎng)選取的不合適,TR 會(huì)發(fā)散。 雖然數(shù)值阻尼可以增強(qiáng)穩(wěn)定性,但當(dāng)把耗散方法用于非線(xiàn)性問(wèn)題求解時(shí),其結(jié)果仍可能不收斂。

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

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

下面以CEM 為例,簡(jiǎn)述保能量方法的實(shí)現(xiàn)過(guò)程。 考慮如下非線(xiàn)性保守系統(tǒng):

式(17)是所有保能量方法的設(shè)計(jì)準(zhǔn)則。

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

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

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

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

4 線(xiàn)性多步法和復(fù)合方法

為了綜合提高時(shí)間積分方法的數(shù)值性能,包括精度、效率、耗散和穩(wěn)定性,可以采用線(xiàn)性多步法和復(fù)合方法。 線(xiàn)性多步法屬于單步時(shí)間積分方法中的一種,但為了計(jì)算當(dāng)前時(shí)刻的響應(yīng)xt+Δt,需要利用其他自啟動(dòng)方法計(jì)算前幾個(gè)時(shí)刻的響應(yīng),以獲得足夠的啟動(dòng)信息。 與線(xiàn)性多步法不同,復(fù)合方法是將一個(gè)時(shí)間步劃分成若干個(gè)分步,計(jì)算xt+Δt需要每個(gè)分步的信息。

4.1 線(xiàn)性多步法

工程界目前使用較多的線(xiàn)性多步法[2]有Adams 類(lèi)方法、顯式Adams-Bashforth 方法、隱式Adams-Moulton 方法及預(yù)測(cè)校正類(lèi)格式等。 這些方法雖然具有高階精度,但穩(wěn)定性難以提高。 Dahlquist[3]曾證明,所有顯式方法和二階以上的線(xiàn)性多步法均無(wú)法實(shí)現(xiàn)無(wú)條件穩(wěn)定。 因此,在實(shí)際應(yīng)用中,顯式和隱式線(xiàn)性多步方法常作為預(yù)測(cè)和校正方法成對(duì)使用,并結(jié)合自適應(yīng)步長(zhǎng)技術(shù),控制局部截?cái)嗾`差,以避免結(jié)果發(fā)散。 此外,向后差分格式,如Houbolt 方法[59]和Park 方法[1],也屬于經(jīng)典的線(xiàn)性多步法。 它們一般具備二階精度、無(wú)條件穩(wěn)定性及強(qiáng)烈的高頻耗散能力,可以快速濾掉高頻響應(yīng)成分,因此更適用于剛性動(dòng)力學(xué)系統(tǒng)的分析。

為進(jìn)一步提高線(xiàn)性多步法的計(jì)算精度,學(xué)者們構(gòu)造出了更加精確的線(xiàn)性三步法[33]和線(xiàn)性四步法[33]。 另外,為了解決線(xiàn)性多步法無(wú)法自啟動(dòng)的問(wèn)題,通過(guò)引入一些輔助狀態(tài)變量,學(xué)者們還提出了與一些線(xiàn)性多步法具有相同數(shù)值性能的單步法[32-33]。

下面以線(xiàn)性?xún)刹椒╗33](linear two-step method, LTM)為例,介紹該類(lèi)方法的求解流程。 LTM的遞推格式為

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

4.2 復(fù)合方法

復(fù)合方法在一個(gè)時(shí)間步的不同分步內(nèi)使用不同的時(shí)間積分方法,以發(fā)揮不同方法各自的優(yōu)勢(shì),從而獲得優(yōu)于各分步方法的數(shù)值性能。

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

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

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

式中:θ0、θ1和θ2為算法參數(shù)。

聯(lián)立求解t+γΔt和t+Δt時(shí)刻的動(dòng)力學(xué)平衡方程,即可得到一個(gè)時(shí)間步內(nèi)所有狀態(tài)變量。

另外,針對(duì)由TR 和BDF(或BIF)復(fù)合得到的時(shí)間積分方法,邢譽(yù)峰和季奕等[42]研究了如何選擇分步數(shù)、分步步長(zhǎng)、差分點(diǎn)數(shù)、各分步采用的方法,并得到了一些重要的結(jié)論,為復(fù)合方法的構(gòu)造提供指導(dǎo),這些結(jié)論包括:

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

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

3) 當(dāng)分步數(shù)為常數(shù)時(shí),只要采用相同的參數(shù)優(yōu)化方式,就可以構(gòu)造出性能相近的復(fù)合方法。如三分步優(yōu)化Optimal-TR-BDF-BDF 方法中采用TR 的分步的占比大約是三分步優(yōu)化Optimal-TRTR-BDF 方法中的一半,但二者的計(jì)算精度幾乎相同。

5 BN 穩(wěn)定型方法

引入數(shù)值阻尼是改善時(shí)間積分方法在非線(xiàn)性系統(tǒng)中穩(wěn)定性的一種直接且有效的方式。 與非耗散方法相比,一般情況下耗散方法在非線(xiàn)性系統(tǒng)中具有穩(wěn)定性?xún)?yōu)勢(shì),但其精度變低。 為了保證在計(jì)算非線(xiàn)性動(dòng)力學(xué)響應(yīng)時(shí)的穩(wěn)定性,保能量方法是一種選擇。 但保能量方法具有額外的計(jì)算量,并且適用范圍有限。

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

式中:F為包含內(nèi)力和外力的向量。

由式(29)和式(30)可求解出t+c1Δt和t+c2Δt兩個(gè)時(shí)刻的狀態(tài)變量。 時(shí)間步終點(diǎn)t+Δt的信息是由這2 個(gè)時(shí)刻的狀態(tài)變量加權(quán)得到,即

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

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

為了確保ρ∞-TSSBN 方法能夠穩(wěn)定地處理非線(xiàn)性系統(tǒng),即具有BN 穩(wěn)定性,利用BN 穩(wěn)定性理論的充分條件(即代數(shù)穩(wěn)定性)設(shè)計(jì)了該算法參數(shù)。 當(dāng)ρ∞=1 時(shí),算法參數(shù)為

數(shù)值算例表明,在TR 和ρ∞-Bathe 方法失效的非線(xiàn)性問(wèn)題中,ρ∞-TSSBN 方法可以給出穩(wěn)定且精確的結(jié)果。

6 結(jié)束語(yǔ)

本文介紹了用于求解動(dòng)力學(xué)常微分方程的時(shí)間積分方法,并著重介紹了近年來(lái)發(fā)展出的一些先進(jìn)時(shí)間積分方法,包括參數(shù)方法、高階無(wú)條件穩(wěn)定方法、保能量方法、線(xiàn)性多步法、復(fù)合方法和BN穩(wěn)定型方法。 更詳細(xì)的資料讀者可以閱讀相關(guān)的綜述文章[63]和專(zhuān)著[64]。

關(guān)于時(shí)間積分方法的未來(lái)發(fā)展問(wèn)題,可以從其求解的是線(xiàn)性問(wèn)題和非線(xiàn)性問(wèn)題角度進(jìn)行討論。

對(duì)于線(xiàn)性動(dòng)力系統(tǒng),時(shí)間積分方法的設(shè)計(jì)已有成熟的理論基礎(chǔ)。 雖然已經(jīng)存在高精度、快速計(jì)算方法[65-67],但如何更快、更精確地計(jì)算仍然是人們追求的目標(biāo)。 對(duì)于線(xiàn)性單自由度動(dòng)力系統(tǒng),已經(jīng)可以利用時(shí)間積分方法得到精確的結(jié)點(diǎn)位移[68-69],其具有精確的幅值和相位。 對(duì)于線(xiàn)性多自由度系統(tǒng),雖然辛幾何算法[70]具有保幅值或保能量特性,或者說(shuō)從根本上解決了幅值誤差累計(jì)問(wèn)題,但相位誤差累計(jì)問(wèn)題仍然存在。 幅值有界準(zhǔn)則和能量準(zhǔn)則都是針對(duì)幅值衰減和發(fā)散問(wèn)題提出的,但建立一個(gè)針對(duì)相位誤差的約束準(zhǔn)則是值得探索的。

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

現(xiàn)代科學(xué)發(fā)展的特征是多學(xué)科交叉、多物理場(chǎng)耦合和多尺度分析和非線(xiàn)性特征,對(duì)于結(jié)構(gòu)動(dòng)力學(xué)系統(tǒng)時(shí)間積分方法的研究工作也要與現(xiàn)代科學(xué)發(fā)展的步伐協(xié)調(diào)一致。