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

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

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

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)的計算可以由精細(xì)積分達(dá)到計算機(jī)精度，因此數(shù)值誤差主要來源于Duhamel積分項(xiàng)的數(shù)值計算誤差。本文采用時域微分求積法計算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ù)估計算可以省略，從而可以極大地提高計算效率。

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ì)庫塔法的計算特性相似。

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

以橢圓積分得到該問題的解析解為基準(zhǔn)，分別采用本文方法(DQM4)和精細(xì)庫塔法計算幅角v1，計算步長Δ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)的計算結(jié)果列入表1和表2。圖2是精細(xì)庫塔法和本文方法計算的幅角v1的變化曲線?？梢钥闯觯?xì)庫塔法的數(shù)值解已經(jīng)發(fā)散，而本文方法可以穩(wěn)定地長時間計算。這是因?yàn)橛嬎鉊uhamel積分項(xiàng)的經(jīng)典四階RK方法是條件收斂的，當(dāng)計算誤差積累到一定程度時，積分結(jié)果將嚴(yán)重偏離系統(tǒng)真實(shí)解。相反，由于s階時域微分求積法是無條件穩(wěn)定的，因此可以采用大步長計算Duhamel積分項(xiàng)，而不會發(fā)生數(shù)值發(fā)散現(xiàn)象。

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

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

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

算例3受迫Duffing方程為

圖2 幅角時程曲線

Fig.2 Time history curve of pendulum angle

圖3 長時間計算結(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)。采用本文方法計算該方程的數(shù)值解，繪制出相應(yīng)的相圖(局部)如圖4和圖5所示，計算到400 s為止。

可以看出，振幅f取0.35和0.41時分別對應(yīng)于系統(tǒng)的一周期和四周期響應(yīng)。相應(yīng)的計算結(jié)果與其他方法對比情況列入表3，其中四階顯式 Runge -Kutta 的計算步長為Δt=1.0 ms，本文方法和預(yù)估校正-辛?xí)r間子域法的計算步長為Δt=0.1 s。由表3可知，本文方法的大步長計算結(jié)果與小步長的四階RK方法吻合得非常好，且比預(yù)估校正-辛?xí)r間子域法的計算精度高，說明本文方法具有較高計算精度和良好的數(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計算結(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)，并給出了通用的計算公式。不僅解決了微分求積法求解非線性常微分方程存在的巨大計算量問題，而且為精細(xì)求積法應(yīng)用于求解非線性動力系統(tǒng)提供了一種新的思路。本文方法可以避免狀態(tài)矩陣求逆，無需對非線性項(xiàng)求導(dǎo)數(shù)，計算格式緊致，便于編程實(shí)現(xiàn)。對兩自由度的線性微分動力系統(tǒng)、非線性單擺和Duffing方程的計算分析表明，本文方法具有較好適應(yīng)性，且計算精度比現(xiàn)有的單步法和預(yù)估校正-辛?xí)r間子域法高，即使采用大步長積分，經(jīng)歷長時間的積分過程仍能保持較好的計算精度和數(shù)值收斂性。因此，本文建立的精細(xì)積分-微分求積法是多自由度結(jié)構(gòu)體系的非線性動力學(xué)分析的一種有效方法。