非齊次動力學(xué)方程的一種精細(xì)積分單步方法

王 永，馬 駿，李靖翔，陳汝斯，許 楊，郝躍東，胡 鵬

(1.國網(wǎng)上海市電力公司特高壓換流站分公司，上海 201413; 2.國網(wǎng)湖北省電力有限公司，武漢 430070;3.南方電網(wǎng)超高壓公司廣州局，廣州 510000; 4.國網(wǎng)湖北省電力有限公司電力科學(xué)研究院，武漢 430070)

1 引 言

目前，針對非齊次動力學(xué)方程數(shù)值解的研究成果十分豐富。常用的數(shù)值積分方法，如中心差分法、Wilson類方法、Houbolt法、Newmark類方法以及四階經(jīng)典Runge-Kutta法等，雖然應(yīng)用十分廣泛，但是在長時間的積分過程中會表現(xiàn)出能量耗散和相位誤差[1]。

鐘萬勰[2]提出的精細(xì)積分法，由于其具有高精度和強(qiáng)穩(wěn)定性的特點(diǎn)，為非齊次動力方程的時程積分提供了一種新的技術(shù)途徑。對于齊次動力方程，精細(xì)積分法可以得到計(jì)算機(jī)精度的解；而對于非齊次方程，由于需要求解Duhamel積分項(xiàng)的數(shù)值解從而會引入數(shù)值誤差[3]。針對Duhamel積分項(xiàng)的數(shù)值計(jì)算問題，研究人員提出了不同的技術(shù)途徑[4]。(1) 近似處理，采用一般多項(xiàng)式、正余弦函數(shù)、指數(shù)函數(shù)及其組合形式近似代替非齊次項(xiàng)，并且推導(dǎo)出相應(yīng)精細(xì)積分的遞推格式； (2) 增維處理，將非齊次項(xiàng)也視為結(jié)構(gòu)動力方程的狀態(tài)變量，化非齊次為齊次，從而在源頭上避免Duhamel積分項(xiàng)的數(shù)值誤差； (3) 數(shù)值積分，采用常見的數(shù)值積分公式近似計(jì)算Duhamel積分項(xiàng)，從而衍生出多種精細(xì)積分格式，其中比較典型的有高斯精細(xì)法[5]、精細(xì)庫塔法[6]以及高精度直接積分法[7]； (4)響應(yīng)矩陣法，將Duhamel積分項(xiàng)轉(zhuǎn)化為一系列響應(yīng)矩陣的精細(xì)積分，無需進(jìn)行矩陣求逆，也適用于非線性問題[8]。

近似處理法要求非齊次項(xiàng)本身具有特定的形式，且需要對系統(tǒng)矩陣進(jìn)行求逆運(yùn)算，因此具有很大的局限性和數(shù)值不穩(wěn)定風(fēng)險[3]。增維精細(xì)積分法[9]對于當(dāng)Duhamel積分項(xiàng)為非線性和時變時的計(jì)算精度不高[10]。高斯精細(xì)法雖然計(jì)算精度較高，但是需要計(jì)算每個積分節(jié)點(diǎn)上的矩陣指數(shù)[8]，計(jì)算效率不高。響應(yīng)矩陣法針對Duhamel積分項(xiàng)為非多項(xiàng)式和非正余弦函數(shù)也存在精度不高的問題[10]。此外，富明慧等[11]提出的廣義精細(xì)積分能夠得到Duhamel積分項(xiàng)計(jì)算機(jī)上的精確解，但是對于其形式有限制且對非線性問題無法直接使用。為此，王海波等[12]采用拉格朗日插值公式提出了對于非線性動力學(xué)分析的廣義精細(xì)積分，但是計(jì)算中需要進(jìn)行預(yù)估-校正；為長時間保持系統(tǒng)辛結(jié)構(gòu)，學(xué)者將辛算法與精細(xì)積分法相結(jié)合提出了辛精細(xì)積分法[13]，但是該算法額外地增加了大量矩陣求逆運(yùn)算。為減少矩陣求逆計(jì)算量，都琳等[14]提出將非齊次方程按湊微分的思想近似化為齊次形式，從而減小了矩陣求逆運(yùn)算量。

本文從數(shù)值積分的角度出發(fā)，采用微分求積法計(jì)算Duhamel積分項(xiàng)的數(shù)值解，并采用與微分求積法同階的顯式 Runge -Kutta 法計(jì)算單步積分中的內(nèi)點(diǎn)近似值，從而將微分求積法顯式化，避免了狀態(tài)矩陣求逆。算例計(jì)算結(jié)果表明，該方法繼承了精細(xì)積分法和微分求積法各自的高精度和強(qiáng)穩(wěn)定性，計(jì)算效率高，比較適用于非線性結(jié)構(gòu)動力學(xué)方程的數(shù)值計(jì)算。

2 精細(xì)積分單步法

結(jié)構(gòu)動力學(xué)方程的一般形式可表達(dá)為

(1)

式中M，C和K分別表示n維方陣，分別代表該系統(tǒng)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣，n表示系統(tǒng)的自由度數(shù)，f(t)表示施加于系統(tǒng)的外力項(xiàng)，n維向量x表示質(zhì)點(diǎn)的位移。

(2)

式中H為n階常數(shù)矩陣，r為非線性廣義外力項(xiàng)。

(3)

在一個積分步長[tk,tk + 1]內(nèi)，非齊次動力方程(2)的解可表達(dá)為同解積分方程(4)

(4)

式中 等式右邊第一項(xiàng)中的指數(shù)矩陣T=eHΔt可以采用加法定理精細(xì)積分得到，第二項(xiàng)積分與系統(tǒng)特征有關(guān)，稱為Duhamel積分項(xiàng)。由于第一項(xiàng)的計(jì)算可以由精細(xì)積分達(dá)到計(jì)算機(jī)精度，因此數(shù)值誤差主要來源于Duhamel積分項(xiàng)的數(shù)值計(jì)算誤差。本文采用時域微分求積法計(jì)算Duhamel積分項(xiàng)。

(i=1，2，…，s) (5)

式中 Δt是積分步長；cj表示網(wǎng)格點(diǎn)；ai j和bj是與網(wǎng)格點(diǎn)相關(guān)的積分系數(shù)[16]。

(6)

式中ti=tk+i×Δt/s。

(7)

式中S1=Hvk+r(tk,vk)

(8)

T1=eH 3Δ t/4,T2=eH Δ t/2,T3=eH Δ t /4

且有T2=T3×T3，T1=T2×T3，T=T1×T3

(9)

注意，對于線性動力學(xué)方程，預(yù)估計(jì)算可以省略，從而可以極大地提高計(jì)算效率。

3 算例分析

線性的二自由度動力學(xué)方程：

(10)

其解析解如下,

式中x1(0)=x2(0)=x3(0)=x4(0)=0。

取仿真步長Δt=0.2 s，積分時程為12 s，分別采用本文方法(s=4，簡記為DQM4)和精細(xì)庫塔法(簡記為RK4)求解，結(jié)果相對于解析解的絕對誤差的絕對值(記為R(·))曲線如圖1所示。

可以看出，對于線性、非剛性的非齊次動力方程，本文方法在精度上略高于精細(xì)庫塔法。盡管經(jīng)典四階RK方法是條件穩(wěn)定的，而s級時域微分求積法是無條件穩(wěn)定的，但是長時間的數(shù)值積分后二者都可以精確刻畫系統(tǒng)的位移變化軌跡。這說明對于非剛性系統(tǒng)，同階的本文方法與精細(xì)庫塔法的計(jì)算特性相似。

式中 初值v1(0)=1.0472,v2(0)=0。

以橢圓積分得到該問題的解析解為基準(zhǔn)，分別采用本文方法(DQM4)和精細(xì)庫塔法計(jì)算幅角v1，計(jì)算步長Δt取0.1 s和0.01 s，并與文獻(xiàn)[1]的預(yù)估校正-辛?xí)r間子域法、文獻(xiàn)[4]的基于精細(xì)積分的Taylor級數(shù)單步法(T-PIM)和文獻(xiàn)[17]的高效單步法做對比。相應(yīng)的計(jì)算結(jié)果列入表1和表2。圖2是精細(xì)庫塔法和本文方法計(jì)算的幅角v1的變化曲線?？梢钥闯觯?xì)庫塔法的數(shù)值解已經(jīng)發(fā)散，而本文方法可以穩(wěn)定地長時間計(jì)算。這是因?yàn)橛?jì)算Duhamel積分項(xiàng)的經(jīng)典四階RK方法是條件收斂的，當(dāng)計(jì)算誤差積累到一定程度時，積分結(jié)果將嚴(yán)重偏離系統(tǒng)真實(shí)解。相反，由于s階時域微分求積法是無條件穩(wěn)定的，因此可以采用大步長計(jì)算Duhamel積分項(xiàng)，而不會發(fā)生數(shù)值發(fā)散現(xiàn)象。

圖1 精細(xì)庫塔法與本文方法的精度比較

Fig.1 Accuracy comparison between precise -Kutta method and this method

從表1和表2可知，本文方法的計(jì)算精度要明顯高于預(yù)估校正-辛?xí)r間子域法、T-PIM和文獻(xiàn)[17]的方法，與橢圓積分解保持高度一致。并且隨著計(jì)算時間的積累，這種現(xiàn)象越明顯。計(jì)算步長取Δt=0.1 s，本文方法長時間計(jì)算結(jié)果的局部與橢圓積分解的對比如圖3所示?？梢钥闯?，本文方法可以長時間跟蹤單擺的運(yùn)動軌跡，說明本文方法具有高精度和良好的穩(wěn)定性。

算例3受迫Duffing方程為

圖2 幅角時程曲線

Fig.2 Time history curve of pendulum angle

圖3 長時間計(jì)算結(jié)果對比

Fig.3 Long -term calculation results comparison

在α=0.4,k=-1,β=1和ω=1.3的條件下，取時間步長Δt=0.1 s，振幅f依次取0.35和 0.41，系統(tǒng)的初始值為[v1,v2]=(2.0,0.0)。采用本文方法計(jì)算該方程的數(shù)值解，繪制出相應(yīng)的相圖(局部)如圖4和圖5所示，計(jì)算到400 s為止。

可以看出，振幅f取0.35和0.41時分別對應(yīng)于系統(tǒng)的一周期和四周期響應(yīng)。相應(yīng)的計(jì)算結(jié)果與其他方法對比情況列入表3，其中四階顯式 Runge -Kutta 的計(jì)算步長為Δt=1.0 ms，本文方法和預(yù)估校正-辛?xí)r間子域法的計(jì)算步長為Δt=0.1 s。由表3可知，本文方法的大步長計(jì)算結(jié)果與小步長的四階RK方法吻合得非常好，且比預(yù)估校正-辛?xí)r間子域法的計(jì)算精度高，說明本文方法具有較高計(jì)算精度和良好的數(shù)值穩(wěn)定性。

表2 幅角v1的數(shù)值結(jié)果對比(Δt=0.01 s)Tab.2 Comparison of pendulum angle (Δt=0.01 s)

圖4f=0.35的周期響應(yīng)

Fig.4 Period response forf=0.35

圖5f=0.41的四周期響應(yīng)

Fig.5 Period response forf=0.41

表3 受迫Duffing方程x計(jì)算結(jié)果對比

Tab.3 Comparison ofxfor Duffing equation

AmlitudeTime/sThispaperRunge-KuttaRef.[1]f=0.351000.56296050.5629610.5629692001.33087891.3308841.330878f=0.411000.37703250.3770340.3772602001.05362701.0536331.053778

4 結(jié) 論

本文將結(jié)構(gòu)動力學(xué)方程轉(zhuǎn)化為非齊次動力方程，通過對未知量的顯式預(yù)估，建立了求解非齊次動力方程精細(xì)積分-微分求積法PIDQM(precise integration-differential quadrature method)，并給出了通用的計(jì)算公式。不僅解決了微分求積法求解非線性常微分方程存在的巨大計(jì)算量問題，而且為精細(xì)求積法應(yīng)用于求解非線性動力系統(tǒng)提供了一種新的思路。本文方法可以避免狀態(tài)矩陣求逆，無需對非線性項(xiàng)求導(dǎo)數(shù)，計(jì)算格式緊致，便于編程實(shí)現(xiàn)。對兩自由度的線性微分動力系統(tǒng)、非線性單擺和Duffing方程的計(jì)算分析表明，本文方法具有較好適應(yīng)性，且計(jì)算精度比現(xiàn)有的單步法和預(yù)估校正-辛?xí)r間子域法高，即使采用大步長積分，經(jīng)歷長時間的積分過程仍能保持較好的計(jì)算精度和數(shù)值收斂性。因此，本文建立的精細(xì)積分-微分求積法是多自由度結(jié)構(gòu)體系的非線性動力學(xué)分析的一種有效方法。