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

任爭(zhēng)爭(zhēng) 梅雨辰 李鴻晶

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

引言

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

結(jié)構(gòu)動(dòng)力分析的本質(zhì)是實(shí)現(xiàn)對(duì)動(dòng)力荷載激勵(lì)下的結(jié)構(gòu)運(yùn)動(dòng)微分方程（組）的求解，屬于常微分方程的初值問題。目前用于分析結(jié)構(gòu)動(dòng)力的方法大致可以分為2類：一類是變換方法，如振型疊加法利用振型的正交性和完備性將結(jié)構(gòu)動(dòng)態(tài)響應(yīng)向各階振型分解，再通過疊加各階振型響應(yīng)以獲得結(jié)構(gòu)動(dòng)態(tài)響應(yīng)的結(jié)果，這類方法采用疊加原理，一般只用于線彈性結(jié)構(gòu)動(dòng)力分析，且要求結(jié)構(gòu)具有經(jīng)典阻尼特性；另一類為直接方法，即直接對(duì)結(jié)構(gòu)運(yùn)動(dòng)微分方程進(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)動(dòng)力分析的計(jì)算需求呈爆發(fā)式增長，因而需要尋求高效率的分析方法。微分求積法（DifferentialQuadrature Method，DQM）是由Bellman等（1971，1972）發(fā)展起來的1種求解微分方程的數(shù)值方法，可通過較小的計(jì)算工作量獲得較高的計(jì)算精度。DQM要求求解域比較規(guī)則，在工程分析中一般用于時(shí)間無關(guān)問題的求解。如在結(jié)構(gòu)力學(xué)領(lǐng)域，多用來求解靜力問題和固有振動(dòng)問題等（Bert等，1988，1993，1996，1997；Jang等，1989；Kang等，1995，1996；Liew等，1996a，1996b；Sherbourne等，1991；Striz等，1988；Wang等，1993，1994；Zeng等，2001）。這類問題都屬于邊值問題，即在DQM中計(jì)算的是解函數(shù)對(duì)空間坐標(biāo)的導(dǎo)數(shù)。而結(jié)構(gòu)動(dòng)力分析屬于初值問題，使用DQM對(duì)其求解的研究工作相對(duì)較少，F(xiàn)ung（2001a，2001b）、Liu等（2008）、李鴻晶等（2011a，2011b）和廖旭等（2013）開展過相關(guān)的研究工作。本文在此基礎(chǔ)上發(fā)展1種基于DQM的結(jié)構(gòu)動(dòng)力分析的高精度方法。不同于靜力邊值問題，實(shí)際結(jié)構(gòu)的動(dòng)力響應(yīng)問題有其特殊性，許多情況下時(shí)間跨度很長，像邊值問題一樣一次性對(duì)所有時(shí)間區(qū)域進(jìn)行離散求解，將出現(xiàn)病態(tài)問題而產(chǎn)生錯(cuò)誤。本文借鑒單元法的思想，以期提高結(jié)構(gòu)動(dòng)力分析的計(jì)算效率。

1 微分求積法基本原理

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

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

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

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

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

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

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

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

實(shí)際計(jì)算權(quán)系數(shù)時(shí)往往只需計(jì)算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)用微分求積法時(shí)還需確定采樣網(wǎng)格節(jié)點(diǎn)的位置，大致分為均勻網(wǎng)格點(diǎn)和非均勻網(wǎng)格點(diǎn)2大類。雖然均勻網(wǎng)格點(diǎn)的精度總體上沒有非均勻點(diǎn)高，但是其在處理離散荷載，如地震荷載或風(fēng)荷載時(shí)，可直接使用原始采樣點(diǎn)作為節(jié)點(diǎn)，不需要對(duì)荷載進(jìn)行額外的插值，計(jì)算效率比不均勻網(wǎng)格點(diǎn)更高。

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

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

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

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

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

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

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

對(duì)這些節(jié)點(diǎn)使用微分求積法：

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

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

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

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

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

3 算例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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