微分求積法在結(jié)構(gòu)動力分析中的應(yīng)用1

任爭爭 梅雨辰 李鴻晶

（南京工業(yè)大學(xué)，土木工程學(xué)院，南京 211816）

引言

地震等災(zāi)害環(huán)境的作用使結(jié)構(gòu)產(chǎn)生復(fù)雜的動態(tài)響應(yīng)，可能導(dǎo)致結(jié)構(gòu)失效、破壞甚至倒塌，從而形成災(zāi)害、造成損失。對新建結(jié)構(gòu)進(jìn)行抗災(zāi)設(shè)計及對已建結(jié)構(gòu)進(jìn)行抗災(zāi)加固是防御和減輕工程災(zāi)害及其損失的有效途徑，這就要求認(rèn)識結(jié)構(gòu)在這些災(zāi)害作用下的性態(tài)和響應(yīng)行為，而結(jié)構(gòu)動力分析則是實現(xiàn)這一目標(biāo)的基本手段。

結(jié)構(gòu)動力分析的本質(zhì)是實現(xiàn)對動力荷載激勵下的結(jié)構(gòu)運動微分方程（組）的求解，屬于常微分方程的初值問題。目前用于分析結(jié)構(gòu)動力的方法大致可以分為2類：一類是變換方法，如振型疊加法利用振型的正交性和完備性將結(jié)構(gòu)動態(tài)響應(yīng)向各階振型分解，再通過疊加各階振型響應(yīng)以獲得結(jié)構(gòu)動態(tài)響應(yīng)的結(jié)果，這類方法采用疊加原理，一般只用于線彈性結(jié)構(gòu)動力分析，且要求結(jié)構(gòu)具有經(jīng)典阻尼特性；另一類為直接方法，即直接對結(jié)構(gòu)運動微分方程進(jìn)行求解而不必引入任何假定，這類方法既可用于線彈性結(jié)構(gòu)，也可用于非線性結(jié)構(gòu)的分析，以Newmark-β法（Newmark，1959）等逐步積分法為代表，一般采用數(shù)值求解手段。

由于現(xiàn)代結(jié)構(gòu)不斷向大型化、復(fù)雜化發(fā)展，加之結(jié)構(gòu)精細(xì)化模型的采用，導(dǎo)致結(jié)構(gòu)動力分析的計算需求呈爆發(fā)式增長，因而需要尋求高效率的分析方法。微分求積法（DifferentialQuadrature Method，DQM）是由Bellman等（1971，1972）發(fā)展起來的1種求解微分方程的數(shù)值方法，可通過較小的計算工作量獲得較高的計算精度。DQM要求求解域比較規(guī)則，在工程分析中一般用于時間無關(guān)問題的求解。如在結(jié)構(gòu)力學(xué)領(lǐng)域，多用來求解靜力問題和固有振動問題等（Bert等，1988，1993，1996，1997；Jang等，1989；Kang等，1995，1996；Liew等，1996a，1996b；Sherbourne等，1991；Striz等，1988；Wang等，1993，1994；Zeng等，2001）。這類問題都屬于邊值問題，即在DQM中計算的是解函數(shù)對空間坐標(biāo)的導(dǎo)數(shù)。而結(jié)構(gòu)動力分析屬于初值問題，使用DQM對其求解的研究工作相對較少，F(xiàn)ung（2001a，2001b）、Liu等（2008）、李鴻晶等（2011a，2011b）和廖旭等（2013）開展過相關(guān)的研究工作。本文在此基礎(chǔ)上發(fā)展1種基于DQM的結(jié)構(gòu)動力分析的高精度方法。不同于靜力邊值問題，實際結(jié)構(gòu)的動力響應(yīng)問題有其特殊性，許多情況下時間跨度很長，像邊值問題一樣一次性對所有時間區(qū)域進(jìn)行離散求解，將出現(xiàn)病態(tài)問題而產(chǎn)生錯誤。本文借鑒單元法的思想，以期提高結(jié)構(gòu)動力分析的計算效率。

1 微分求積法基本原理

微分求積法是1種用于求解微分方程的數(shù)值方法。它的實質(zhì)是將函數(shù)在某一離散節(jié)點處的各階導(dǎo)數(shù)值，近似表示成計算域內(nèi)所有節(jié)點處離散函數(shù)值的線性加權(quán)和，從而將復(fù)雜的微分方程化為關(guān)于離散點的線性方程（組）。由于本文討論的是結(jié)構(gòu)動力反應(yīng)的計算，求解的運動微分方程僅是關(guān)于時間t的常微分方程，因此僅介紹一維區(qū)域內(nèi)的微分求積原理。

設(shè)函數(shù)f(x)為在區(qū)間［a，b］上k階連續(xù)可微，將區(qū)間［a，b］劃分為m段，共（m＋1）個互不相同的節(jié)點，分別記為x0，x1，……，xm-1，xm，其中x0＝a，xm＝b。

根據(jù)計算數(shù)學(xué)的函數(shù)逼近理論，函數(shù)f（x）可做如下逼近：

其中，qj(x)為函數(shù)空間中各線性無關(guān)的基函數(shù)。

對式（1）求k階導(dǎo)數(shù)，然后將所有節(jié)點代入，得到：

式（3）即為一維區(qū)域微分求積的基本公式。

常用的函數(shù)空間是（m＋1）維多項式空間，選擇該空間中的所有基函數(shù)都能得到相同的權(quán)系數(shù)。最常見的多項式基函數(shù)是冪指數(shù)插值函數(shù)和拉格朗日插值函數(shù)2種，分別如式（4）、（5）所示：

選用式（4）的冪指數(shù)函數(shù)計算權(quán)系數(shù)，得到權(quán)系數(shù)的隱式表達(dá)式，需要求解范德蒙矩陣的逆矩陣，但當(dāng)m較大時，不但計算量大，而且矩陣將容易出現(xiàn)病態(tài)。為了解決這些問題，目前大多采用式（5）給出的拉格朗日插值基函數(shù)，因為拉格朗日插值的各項系數(shù)為各節(jié)點的函數(shù)值，形式與式（1）完全相同，直接求k階導(dǎo)數(shù)，便可得到相應(yīng)的k階權(quán)系數(shù)的顯示表達(dá)式，計算效率大幅度提高。

實際計算權(quán)系數(shù)時往往只需計算1階導(dǎo)數(shù)的權(quán)系數(shù)，其它各階導(dǎo)數(shù)的權(quán)系數(shù)可由1階權(quán)系數(shù)以矩陣的形式方便地表示：

其中，A表示1階權(quán)系數(shù)矩陣，A（k）表示k階權(quán)系數(shù)矩陣。1階權(quán)系數(shù)可通過選定拉格朗日基函數(shù)，直接得到如下的表達(dá)式：

應(yīng)用微分求積法時還需確定采樣網(wǎng)格節(jié)點的位置，大致分為均勻網(wǎng)格點和非均勻網(wǎng)格點2大類。雖然均勻網(wǎng)格點的精度總體上沒有非均勻點高，但是其在處理離散荷載，如地震荷載或風(fēng)荷載時，可直接使用原始采樣點作為節(jié)點，不需要對荷載進(jìn)行額外的插值，計算效率比不均勻網(wǎng)格點更高。

2 結(jié)構(gòu)動力反應(yīng)微分求積分析方法

用上述微分求積法求解結(jié)構(gòu)的動力反應(yīng)，并以線彈性單自由度體系為例進(jìn)行討論。雖然實際結(jié)構(gòu)大都為多自由度體系，且為非線性體系，但其動力反應(yīng)都可以通過線性迭代和振型分解轉(zhuǎn)化為線彈性單自由度體系。因而，討論線彈性單自由度體系的微分求積法更具普遍意義，且簡單直觀、易于理解。

線彈性單自由度體系的動力反應(yīng)運動方程為：

其中，ω、ξ分別表示體系的自振頻率和阻尼比；u表示體系的位移；分別表示體系的速度和加速度；p（t）為結(jié)構(gòu)所受到的隨時間變化的荷載。

在實際工程中，動荷載隨類型的不同，作用時間差異很大，短則瞬間（如沖擊荷載），長則幾分鐘甚至幾十分鐘（如風(fēng)荷載）。對于一些長時間作用的荷載，在整個荷載作用時域內(nèi)實施微分求積法，要保證計算結(jié)果的精確性，需要成千上萬的時間節(jié)點，這將導(dǎo)致求解時系數(shù)矩陣的階數(shù)過于龐大，引起矩陣的條件數(shù)大，病態(tài)效果嚴(yán)重，無法得到可靠的結(jié)果。為了解決此問題，可以借鑒有限單元法的思想，將整個荷載持時劃分為多個等間距的時間單元，稱為1個時步，在每個時步內(nèi)，使用微分求積法求解。

考慮動力荷載持時內(nèi)長度Δt的時步［tj，tk］內(nèi)的反應(yīng)，tj、tk一般與荷載p（t）的采樣時刻點重合。在時步內(nèi)定義一局部坐標(biāo)tt∈ [0，Δt]，坐標(biāo)起點為tj，方向同時間坐標(biāo)。為了方便處理，正則化局部坐標(biāo)，作則定義域被正則化為［0，1］，時步內(nèi)關(guān)于τ的運動方程可寫為：

將時步離散為m段，記節(jié)點為τ0，τ1，……，τm，將分別簡記為則在各離散節(jié)點處的方程為：

對這些節(jié)點使用微分求積法：

其中，aij是微分求積的權(quán)系數(shù)，將式（11）、（12）寫成矩陣的形式：

已知時步的初始位移和速度分別為u0和0v，且將式（11）、（12）改寫為：

將式（15）、（16）代入式（10），并將已知量與未知量分離在等式的兩側(cè)，整理后得到如下方程：

式（17）為各時步內(nèi)微分求積的基本方程，求解該方程可得到時步內(nèi)各點的位移反應(yīng)，再運用微分求積原理，可進(jìn)一步求得各節(jié)點速度反應(yīng)：

將本時步末的位移um和速度作為下一時步的初始位移和速度，從初始0時刻開始逐時步進(jìn)行求解，可得到荷載作用的所有時刻的位移和速度反應(yīng)。除了位移和速度外，工程上關(guān)心的加速度可通過運動方程（8）變形得到：

3 算例

通過具體的算例驗證上述微分求積法求解結(jié)構(gòu)動力反應(yīng)計算結(jié)果的精確性和可靠性。選擇3種不同自振周期的單自由度體系，其頻率范圍大致覆蓋低頻、中頻和高頻，阻尼比都取為工程中常見的0.05，在上面施加不同頻率的正弦荷載，體系的基本參數(shù)如表1所示，荷載信息如表2所示。

表1 體系的基本特性Table 1 Basic characteristics of systems

表2 簡諧荷載的信息Table 2 Information of simple harmonic load

用微分求積法求解3種體系在以上簡諧荷載下的動力反應(yīng)。為了方便計算，時步內(nèi)采用均勻網(wǎng)格離散方案，時步的長度Δt取為簡諧荷載的周期。由于1個周期的長度是簡諧荷載的最小重復(fù)單元，這樣取值不僅能方便地分析時步分段數(shù)m對數(shù)值穩(wěn)定性和精度的影響，更能清楚地發(fā)現(xiàn)荷載周期與步長Δt取值的規(guī)律。

首先，計算該時步長度下采用不同m所得到的體系的位移和速度反應(yīng)，然后采用體系在簡諧荷載下的解析解作為精確解進(jìn)行校核。表3—5列出了分段數(shù)m為20內(nèi)的偶數(shù)時，簡諧荷載激勵下3種體系動力反應(yīng)DQM解與精確解的平均相對誤差。

表3 荷載周期1s的平均相對誤差Table 3 Average relative error with load period of 1s

表4 荷載周期0.2s的平均相對誤差Table 4 Average relative error with load period of 0.2s

表5 荷載周期為0.1s的平均相對誤差Table 5 Average relative error with load period of 0.1 seconds

從表3—5可以看出，當(dāng)時步分段數(shù)m很小時，DQM計算結(jié)果嚴(yán)重偏離精確解，但隨著m的增大，除少部分點外，誤差大致呈迅速減小的趨勢，個別數(shù)據(jù)（如m＝14時）反應(yīng)誤差巨大是由于算法在該時步分段下，時步長度取為荷載周期時出現(xiàn)了數(shù)值不穩(wěn)定現(xiàn)象，初始誤差隨著逐時步推進(jìn)不斷放大，最后完全湮滅真實的結(jié)果。當(dāng)m增大到10時，3種體系在不同周期的簡諧荷載下的位移和速度反應(yīng)都減小到了5%以內(nèi)，對于實際工程已足夠精確。為了更形象地描述這種精確程度，圖1—3給出了0—20s內(nèi)荷載周期為1s下，m＝10時DQM的位移計算結(jié)果與精確解，為了更清晰地進(jìn)行對比，對每張圖進(jìn)行了局部的放大。

圖1 簡諧荷載激勵下體系1的位移反應(yīng)Fig.1 Displacement response of system 1 under simple harmonic load

圖2 簡諧荷載激勵下體系2的位移反應(yīng)Fig.2 Displacement response of system 2 under simple harmonic load

圖3 簡諧荷載激勵下體系3的位移反應(yīng)Fig.3 Displacement response of system 3 under simple harmonic load

從圖中可以看出，DQM的計算結(jié)果與精確解幾乎完全重合，充分說明了分段數(shù)m＝10時，DQM計算動力反應(yīng)足夠精確，且數(shù)值穩(wěn)定性能得到很好的保證。以地震荷載為例，中低頻地震荷載的大部分周期一般在0.2s以上，偏于保守取為0.2s，分10段后節(jié)點間距為0.02s。而目前一般的地震地面加速度采樣周期僅為0.005s或0.1s，節(jié)點間距取采樣周期的幾倍以上，都能得到精確的結(jié)果，計算精度較高。

對比表3、4、5可以發(fā)現(xiàn)，當(dāng)m由10繼續(xù)向上增大時（排除m＝14時因失穩(wěn)誤差劇增的情況），誤差基本上呈繼續(xù)減小的趨勢，很多情況下誤差甚至降至1%或1‰以內(nèi)。但這對工程實際已無太大的意義，反而會隨著m的增大，不斷地增大計算量。因此，對于結(jié)構(gòu)在簡諧荷載作用下的反應(yīng)，采用均勻節(jié)點方案，即在1個荷載周期內(nèi)取m＝10，時步內(nèi)插入9個內(nèi)節(jié)點，既能達(dá)到土木工程領(lǐng)域內(nèi)精度的要求，又能最大限度地減少計算工作量，是相對較優(yōu)的時步節(jié)點數(shù)量。而且，以往DQM求解工程結(jié)構(gòu)靜力問題的大量經(jīng)驗表明，m取8—10可得到滿意的計算結(jié)果（王鑫偉，1995），這與本文得出的m＝10的取值較優(yōu)也相符，這雖然是DQM解決靜力問題的經(jīng)驗，但在1個周期內(nèi)的動力問題與靜力問題存在不少的聯(lián)系，從而也可間接說明本文將m＝10作為較優(yōu)參數(shù)的合理性。

由表3—5還可以發(fā)現(xiàn)，對于簡諧荷載激勵下的反應(yīng)，當(dāng)m＝10時，取時步長度Δt與荷載周期相等時，計算誤差已在工程可接受的范圍內(nèi)。因此，計算體系在簡諧荷載激勵下的反應(yīng)，在選定m＝10的前提下，將時步長度和簡諧荷載的周期保持一致即可。

雖然簡諧荷載是1種理想的荷載，但是對計算實際荷載激勵下的反應(yīng)具有重要的意義，這是因為工程中實際的動荷載都是一定持時的暫態(tài)荷載，都可以通過周期延拓表示成一系列簡諧荷載的疊加。實際計算時，可以首先采用諧波分析，檢測出其包含的不同周期的簡諧波的幅值，得到該動力荷載的頻譜；然后以最大幅值所對應(yīng)的周期，即卓越周期為基準(zhǔn)，確定等效周期；最后，取時步的長度為荷載的等效周期，將時步等分成10段進(jìn)行計算，即可得到較精確的計算結(jié)果。

4 結(jié)論

本文將微分求積法引入結(jié)構(gòu)動力反應(yīng)的分析與計算，并通過數(shù)值算例得到如下的結(jié)論：

（1）采用微分求積法求解結(jié)構(gòu)在動力荷載激勵下的反應(yīng)合理可行。在較大的節(jié)點距離下依然能夠得到精確的結(jié)果，計算精度高，且具有普適性，對不同自振周期的結(jié)構(gòu)、不同頻率的動力荷載都適用。

（2）用DQM進(jìn)行動力分析時，能一次性求得多個時刻的反應(yīng)。相比傳統(tǒng)的每次只能求1個時刻的逐步積分法，計算效率得到提高，計算成本也得到了降低。

（3）在時步長度Δt一定時，DQM求解動力反應(yīng)的計算精度和數(shù)值穩(wěn)定性與時步分段數(shù)m有關(guān)。排除一些失穩(wěn)飄移的情況，一般m越大，計算精度越高，但計算量也越大。綜合考慮，對于均勻網(wǎng)格離散方案，實際計算時取m＝10相對較優(yōu)。

（4）使用DQM進(jìn)行實際結(jié)構(gòu)動力反應(yīng)分析時，時間步長Δt可選為動力荷載的等效周期，然后將各時步等分成10段來計算，這樣可獲得較滿意的計算結(jié)果。