線性動力分析的一種通用積分格式

王海波，何崇檢,賈耀威

(中南大學 土木工程學院，長沙 410075)

結構的動力分析一直是工程研究的熱門領域。傳統的直接積分法，如中心差分法、Newmark法及wilson-θ法等在求解結構動力方程方面有著廣泛的應用。但這些方法在建立積分格式時缺乏相應的變分原理的保證，因此存在能量耗散與相位誤差，不能很好地滿足工程實際的需要。

1994年前后，鐘萬勰[1]放棄了傳統直接積分法中的差分格式，效仿哈密頓體系對偶變量的引入，提出了精細時程積分法，為結構的動力分析開辟了嶄新的途徑。對于線性齊次方程，精細積分法能利用其遞推格式給出計算機意義上的精確解。而對于非齊次動力方程，鐘萬勰采用一次多項式來近似荷載項，其計算精度不是很高，且需要對狀態矩陣求逆。Lin等[2]分別就荷載項為三角函數及傅里葉級數的情況給出了其解析解的表達式。文獻[3-6]提出了增維精細積分法，該方法將非齊次項也看做動力方程的狀態變量，將其納入到求解中去，從而將非齊次動力方程轉化為齊次動力方程求解。該類方法具有很高的計算精度，且不需要對狀態矩陣求逆，但該類算法只適用于荷載項為幾類特殊形式的情況。文獻[7-9]提出了一類利用數值積分法，如辛普森求積公式等來求解非齊次項的動力響應的方法，但該類方法的計算精度略低，無法與通解的精度相匹配,而如若要提高計算精度，則需要增加積分點的個數，其計算量會隨之急劇增加,因此，該類方法的計算精度與效率之間存在難以協調的矛盾。富明慧等[10]結合矩形積分公式與遞推算法，提出了廣義精細積分法。該方法不需要對狀態矩陣求逆，且特解與通解具有相同的計算精度，具有十分優越的算法的特性。儲德文等[11]討論了精細積分法中積分方法選擇的問題，指出了為保證精細算法的高精度，應根據荷載的性質選擇適當的方法。

本文結合泰勒級數展開式和廣義精細積分法，提出了一種避免狀態矩陣求逆的線性動力分析的通用積分格式。該方法可用于任意激勵下結構的動力響應的求解，具有廣泛的適用性。數值算例證明了本文方法的有效性。

1 動力方程的廣義精細積分法

線性定常系統的動力方程為

(1)

式中：M，C和K分別為質量陣、阻尼陣和剛度陣;X為位移向量;F為荷載向量。式(1)的初始條件為

(2)

(3)

式中:v,H及r(t)有兩種表達方式

(4)

或

(5)

式(4)為鐘萬勰效仿哈密頓體系引入偶變量的形式，式(5)為稍后發展出來的另一種形式。式(5)因為能同時反映系統的位移與速度，其應用更廣泛一些。

利用指數矩陣，可將式(3)轉化為同解積分方程

(6)

式中：Δt為時間步長。式(6)中的第一項為非齊次動力方程的通解，可用精細積分法精確得到。第二項為非齊次動力方程的特解，是本文討論的重點。令

(7)

將時間步長離散為m=2N等份，N=20，令τ=Δt/m。由廣義精細積分法可知，當r(t)為多項式時，可采用矩形積分公式與遞推算法精確高效地求出式(7)。以三次多項式為例來說明這一過程，記

r(t)=r0+r1(ti+1-t)+
r2(ti+1-t)2+r3(ti+1-t)3

(8)

將式(8)代入式(7)，可得

(9)

其中,

(10)

式中：k=0,1,2,3,做變量替換s=ti+1-t，并注意到ti+1-ti=Δt，則式(10)可化為

(11)

Sk與所在的區段無關，在整個計算過程中只需計算一次。記Y=eHτ，利用矩形積分公式對式(11)中各項進行數值積分，得

(12)

(13a)

(13b)

(13c)

(13d)

不難驗證遞推公式

(14)

(15a)

(15b)

(15c)

(15d)

不難驗證遞推公式

(16)

式中：Ak為修正系數矩陣，其表達式為

(17)

Ak的選取項數與所需的計算精度有關，選取的項數越多，其修正的精度就越高。

2 特解的通用積分格式

廣義精細積分法的核心思路在于矩形積分公式與遞推算法，本文的筆者曾試圖將其推廣到其它荷載形式的情況，但均以失敗告終。其根本原因在于，只有多項式與指數函數具有這種內在的遞推性質。而將函數采用多項式來擬合，有兩種方法，其一是采用拉格朗日插值多項式來近似，此即傳統數值積分法采用的方式。另一種則是采用泰勒級數展開式來近似荷載函數，此即本文采用的思路。

在式(7)中，做變量替換，令t=ti+1-s，得

(18)

由泰勒中值定理可知，如果函數f(x)在含有x0的某個開區間(a,b)內具有直到n+1階的導數，則對任一x∈(a,b)，有

(19)

式中：Rn(x)為拉格朗日型余項。根據式(19)可將r(ti+1-t)在t=0時刻展開成多項式

(20)

式中：t∈(0,Δt)令

(21)

式中：r為一個常向量；f(ti+1-t)為基本荷載函數，將式(21)代入式(20)可得

(22)

將式(22)代入式(18)中，可得

(23)

式中：Sk,k=0,1,2,…，n即為式(11)中的表達式。而rk,k=0,1,2,…,n的表達式為

r0=f(ti+1)

(24a)

r1=-f′(ti+1)

(24b)

(24c)

(24d)

此時，結合廣義精細積分法即可精確地計算出非齊次項的動力響應。可以發現，式(24)中rk為變量，在計算過程中需要不斷更新，但這只涉及到標量運算，不會顯著增加運算量。下面就幾種常用荷載形式給出其逐步積分公式。值得注意的是，對于荷載項為正、余弦函數的情況，富明慧等將其轉化為復指數函數，也可利用廣義精細積分法求解出結構的動力響應。考慮到復指數函數運算的復雜性與本文方法的完備性，下面仍給出荷載項為正、余弦函數時的逐步積分公式。

(1)非齊次項為sint的情況

(25a)

(2)非齊次項為cost的情況

(25b)

(3)非齊次項為et的情況

(25c)

(4)非齊次項為ln(1+t)的情況

(25d)

(5)非齊次項為(1+t)a的情況

(25e)

(6)非齊次項為tet的情況

(25f)

式(25)給出了高階導數具有通式的荷載項的積分公式，對于其它一般的荷載項，只需求出其前幾階的導數形式代入式(22)即可求出其非齊次項的動力響應。因此，本文提出的算法具有廣泛的適用性。

3 計算過程

綜上所述，線性系統在任意激勵下的動力響應求解可歸納為以下步驟：

步驟2 在各區間[ti,ti+1]上將非齊次項展開成式(22)的形式；

可以發現，步驟1在整個計算過程中只需要計算一次，這也是遞推算法的高效性所在。此外，修正后的系數矩陣具有很高的精度，因此，本文算法的主要誤差來源于非齊次項的擬合誤差。理論上，通過選取冪級數的項數，可以達到任意精度。

4 算例分析

表1 精細數值積分結果與精確解比較Tab.1 Comparison between the result of precise numerical integration and exact solution

算例2考慮線性二自由度動力學方程

取初值v1(0)=2.5，v2(0)=v3(0)=v4(0)=1.0，本文算法在時間步長Δt=0.1 s,Δt=0.2 s,Δt-0.5 s時的前30 s分析結果見表2。

本算例采用五階冪級數來擬合非齊次項。由表2可以看出，當Δt=0.1 s時，數值解基本可以保證七位有效數字；當Δt=0.2 s時，數值解基本可以保證五位有效數字；而當Δt=0.5 s時，數值解可以保證三位有效數字。這是由于數值解的精度主要取決于非齊次項的擬合精度，時間步長越小，其擬合的精度就越高，從而其數值解的精度也就越高。值得注意的是，由于本文算法采用的是利用積分區間端部時刻的冪級數來擬合非齊次項，因此選取的時間步長不宜過大。

表2 v1的數值結果比較Tab.2 Comparisons of numerical result of v1

算例3考慮動力方程

取初值v1(0)=2.5，v2(0)=0，v3(0)=v4(0)=1.0。該動力方程的解析解為

取時間步長Δt=0.5 s，各種方法的計算結果見表3。

本算例選用的時間步長較大，故此處選用六階冪級數來擬合非齊次項。由表3可以看出，當時間步長較大時，HPD-L的線性化假設不能很好地擬合非齊次項，會帶來很大計算的誤差。增維法此時能得到精確解，且由于其計算格式簡單，所以雖然增加了矩陣的維數，其計算時間仍然很短。但值得注意的是，對于非齊次項為其他形式的情況，增維法需要將非齊次項三角函數化，這本身就會帶來較大的誤差，因此其適用范圍有限。同時可以看出，盡管時間步長較大，科茨公式與高斯公式(三階)仍能給出高精度的解，且其計算效率也很高，這說明傳統的數值積分方法能有效地求解非齊次項的動力響應。本文算法由于采用了高階次的冪級數來擬合非齊次項，故其精度很高，但同時，其計算效率也略低。其主要的原因在于計算系數矩陣需要花費較長的時間。但同時也應注意到，系數矩陣只需要計算一次，當時間步數較大時，其簡單的計算格式能彌補這一缺陷。本文算法最大的優勢在于其簡單的計算格式，不必求解積分系數，只需要執行一系列循環語句即可求解出非齊次項的動力響應，同時，通過選取冪級數的項數，可以達到任意精度。

表3 v1,v2的數值結果比較Tab.3 Comparisons of numerical results of v1 and v2

5 結 論

(2)本文的實質是對廣義精細積分法的擴展，通過采用冪級數來擬合非齊次項將其推廣到任意荷載形式的情況。同時，本文算法也可以看作傳統數值積分法的補充，兩者的實質都是采用多項式來擬合非齊次項，而本文算法的優勢在于其簡單的計算格式，不需要求解積分系數，只需要執行一系列循環語句即可求解出非齊次項的動力響應，適宜于編程。