蹇開林，溫偉斌，，駱少明

(1. 重慶大學(xué) 煤礦災(zāi)害動(dòng)力學(xué)與控制國家重點(diǎn)實(shí)驗(yàn)室，重慶 400044； 2. 仲愷農(nóng)業(yè)工程學(xué)院，廣州 510225)

時(shí)間積分方法為求解結(jié)構(gòu)動(dòng)力響應(yīng)尤其大型結(jié)構(gòu)有限元離散計(jì)算的有效手段，該方法需在每個(gè)時(shí)間步上離散化求解平衡方程；而非線性分析中則需在每個(gè)時(shí)間步上進(jìn)行迭代運(yùn)算并更新系統(tǒng)剛度矩陣，對(duì)自由度數(shù)較大系統(tǒng)計(jì)算量較大。因此，研究時(shí)間積分方法時(shí)不僅考慮算法的穩(wěn)定性質(zhì)，亦應(yīng)考慮計(jì)算精度及效率，期望以較少時(shí)間步計(jì)算獲得高精度數(shù)值解。逐步積分法可分顯式積分法[1-3]與隱式積分法[4-6]兩類。前者優(yōu)點(diǎn)在于無需對(duì)每個(gè)時(shí)間步的平衡方程迭代運(yùn)算，計(jì)算量少，且易于編程實(shí)現(xiàn)[7]。其不足之處在于大多數(shù)方法均為條件穩(wěn)定，計(jì)算結(jié)果連續(xù)性低，需采用較小時(shí)間步計(jì)算才能獲得穩(wěn)定、足夠精度的計(jì)算結(jié)果。由于顯式積分法算法實(shí)現(xiàn)上較隱式積分法簡單，已被用于擬動(dòng)力試驗(yàn)研究[8-10]。分段樣條多項(xiàng)式插值近似已被用于計(jì)算機(jī)輔助設(shè)計(jì)、幾何造型、數(shù)據(jù)擬合及有限元計(jì)算[11-12]。Liu[13]用分段Birkhoff插值多項(xiàng)式求解多自由度系統(tǒng)動(dòng)力響應(yīng)；Inoue等[14]給出基于正交B樣條多項(xiàng)式的逐步積分格式, 相比傳統(tǒng)Newmark法、Wilson-θ法，該方法計(jì)算程序簡單、效率更高。Rostami等[15]用四次B樣條基函數(shù)發(fā)展的具有拋物線加速度的顯式積分法，但其條件為穩(wěn)定的。張琳等[16]在Hamilton型擬變分原理體系下建立時(shí)間子域以三、五次B樣條函數(shù)插值的時(shí)間子域法，但其前處理計(jì)算量較大。

鑒于顯式積分法算法穩(wěn)定性與計(jì)算精度的不足，本文提出參數(shù)控制的逐步積分法。該方法采用均勻七次B樣條插值建立節(jié)點(diǎn)位移、速度、加速度試函數(shù)，給出平衡遞推方程；通過控制參數(shù)可使算法絕對(duì)穩(wěn)定。數(shù)值算例結(jié)果表明，數(shù)值位移、速度、加速度計(jì)算精度較高。

1 均勻七次B樣條

本文采用B樣條基函數(shù)插值離散化求解結(jié)構(gòu)動(dòng)力問題微分方程，B樣條基函數(shù)有多種定義方法，為便于計(jì)算機(jī)編程，用逆推方法定義[12]。設(shè){t0,t1,…,ti,ti+1,…}為單調(diào)不減序列，即，ti≤ti+1，(i=0,1,2，…)。其中ti為B樣條節(jié)點(diǎn)。取Bi,k(t)為第i條k次B樣條曲線，定義為

(1)

(2)

式中：[ti,ti+1)為第i個(gè)B樣條節(jié)點(diǎn)區(qū)間。逆推過程中，需約定0/0=0。

對(duì)均勻B樣條基函數(shù)取ti+1-ti=Δt，τi=(t-ti)/Δt，則由逆推定義可推導(dǎo)t∈[ti,ti+1)內(nèi)，插值所用七次均勻B樣條基函數(shù)及導(dǎo)數(shù)表達(dá)式為

(3)

(4)

8(2-τi)7-l+28(1-τi)7-l]

(5)

28(2-τi)7-l-56(1-τi)7-l]

(6)

(7)

(8)

(9)

(10)

式中：l為各B樣條基函數(shù)對(duì)t的導(dǎo)數(shù)階次。

2 基本原理及算法實(shí)現(xiàn)

2.1 動(dòng)力平衡方程遞推格式

單自由度系統(tǒng)結(jié)構(gòu)動(dòng)力學(xué)微分方程可寫為

(11)

取0=t0<…

x(l)(t)=B(l)(τi)di

(12)

式中：

(14)

為便于推導(dǎo)，令

(15)

(16)

(17)

由式(12)有

(18)

式(18)可寫為矩陣形式

(19)

式中：

(20)

(21)

將式(19)代入式(12)有

(22)

將式(22)代入式(11)得殘值方程為

(23)

令

(24)

將式(24)代入式(23)有

(25)

取權(quán)函數(shù)W1(τi)=δ(τi-1)，W2(τi)=δ(τi-r1),W3(τi)=δ(τi-r2)，W4(τi)=δ(τi-r3)。即在點(diǎn)τi=1、τi=r1、τi=r2與τi=r3處采用配點(diǎn)法。其中r1,r2,r3滿足r1≠r2≠r3，用以調(diào)節(jié)計(jì)算穩(wěn)定性與精度。

將式(25)利用加權(quán)殘值法，得

(26)

將上式運(yùn)算簡化為

(27)

式(27)表明本文算法構(gòu)造的逐步積分法具有自起步、無中間計(jì)算環(huán)節(jié)特點(diǎn)。

2.2 遞推平衡方程初始化

mx(3)(t)+cx(2)(t)+kx(1)(t)=f(1)(t)

(28)

令式(28)中t=t0,有

x(3)(t0)=m-1[f(1)(t0)-kx(1)(t0)-cx(2)(t0)]

(29)

式(29)為計(jì)算x(3)(t0)的直接初始化方法。

(30)

(31)

2.3 計(jì)算步驟

由計(jì)算步驟知，矩陣A,Fi的計(jì)算會(huì)較大程度影響計(jì)算效率。式(24)中E(τi)為矩陣A,Fi的主要計(jì)算部分?？紤]對(duì)任意等距時(shí)間單元采用相同系數(shù)r1,r2,r3，因而對(duì)不同時(shí)間步無需重復(fù)計(jì)算，尤其自由度較大系統(tǒng)，亦無需對(duì)E(τi)進(jìn)行大量重復(fù)計(jì)算, 只需將式(24)中m,c,k替換為對(duì)應(yīng)矩陣，并引入Kronecker乘積運(yùn)算即可。其它計(jì)算過程與單自由度系統(tǒng)類似。由此分析表明本文算法計(jì)算量不大。

3 穩(wěn)定性分析

以單自由度微分方程(T/2π)2x(2)(t)+x(t)=0考察算法的穩(wěn)定性，其中T為周期。若式(27)中傳遞函數(shù)矩陣A滿足譜半徑ρ(A)≤1，則算法是穩(wěn)定的。式(27)中矩陣A(r1,r2,r3)會(huì)隨r1,r2,r3的取值變化。為便于分析及工程應(yīng)用，本文取r1=0.9。通過數(shù)值計(jì)算繪出ρ(A)隨Δt/T變化曲線見圖2。由圖2(a)看出，0

圖2 本文算法不同參數(shù)r對(duì)應(yīng)譜半徑

本文算法與傳統(tǒng)Newmark方法、Wilson 方法譜半徑見圖3。由圖3曲線對(duì)比知，本文算法對(duì)低頻段具有更好的保持作用，能有效濾掉高頻分量?？紤]工程應(yīng)用中有限元網(wǎng)格計(jì)算的高頻分量往往不準(zhǔn)確，且對(duì)結(jié)構(gòu)相應(yīng)貢獻(xiàn)較小，因而可通過增大參數(shù)r2,r3能快速有效濾掉高頻分量。

4 數(shù)值算例與分析

4.1 單自由度強(qiáng)迫振動(dòng)

表1 不同方法所求位移及相對(duì)誤差

基于Hamilton擬變分原理建立的B樣條時(shí)間積分法進(jìn)行動(dòng)力學(xué)分析計(jì)算[16]，給出的三次B樣條插值計(jì)算精度僅為傳統(tǒng)方法(Newmark法，Wilson-θ法)的1/17～2/5；而采用五次B樣條計(jì)算時(shí)雖精度有所提高，但前處理復(fù)雜、計(jì)算量大、計(jì)算效率低，與本文算法推導(dǎo)及數(shù)據(jù)對(duì)比可知，本文方法的算法實(shí)現(xiàn)與計(jì)算精度具優(yōu)越性。

表2 不同方法所求速度相對(duì)誤差

表3 不同方法所求加速度相對(duì)誤差

表4為相同MATLAB編程環(huán)境下，取時(shí)間步數(shù)Nt=60 時(shí)各方法計(jì)算時(shí)間。由表4看出，本文算法的初始化時(shí)間遠(yuǎn)高于傳統(tǒng)方法，而主程序計(jì)算時(shí)間少于傳統(tǒng)方法，雖總的計(jì)算時(shí)間多于傳統(tǒng)方法，但計(jì)算精度遠(yuǎn)高于傳統(tǒng)方法，計(jì)算效率仍較高。

表4 不同方法計(jì)算時(shí)間對(duì)比

4.2 簡支梁彎曲振動(dòng)

圖4為兩端簡支等截面梁，長L=6 m, 截面厚h= 2×10-2m, 截面寬b=2×10-2m, 截面面積A=bh，截面慣性矩I=bh3/12, 彈性模量E=210 GPa,泊松比μ=0.3, 密度ρ=4×104kg/m3, 阻尼系數(shù)c=0，橫向載荷q(x,t)=F0sin(ω0t)δ(x-L/2)，F(xiàn)0=1 kN,ω0=4 Hz。該載荷作用下，簡支梁振動(dòng)位移理論解為

圖4 簡支梁彎曲振動(dòng)示意圖

取時(shí)間步長Δt=0.2 s，空間離散采用三次Hermite有限單元，空間單元數(shù)Ns=10。本文算法參數(shù)選r2= 0.85,r3=0.95，其它參數(shù)同前。簡支梁不同時(shí)刻位移、速度、加速度曲線見圖5。由圖5看出，相同時(shí)間步長時(shí)本文算法計(jì)算結(jié)果與理論解吻合較好，明顯優(yōu)于傳統(tǒng)Wilson-θ法、Newmark法。

簡支梁中點(diǎn)(最大變形處)t=iΔt(i=0,1,…,40)時(shí)刻計(jì)算值見圖6。由圖6看出，傳統(tǒng)方法計(jì)算結(jié)果在初始時(shí)間段與理論解基本吻合，但隨時(shí)間的增加逐漸偏離理論曲線，結(jié)果不準(zhǔn)確。而本文算法隨時(shí)間的增加與理論解吻合仍較好。

圖5 給定時(shí)刻不同方法數(shù)值解

圖6 中點(diǎn)位置不同方法數(shù)值解

用不同方法計(jì)算時(shí)間對(duì)比見表5。其中C1為取參數(shù)r2=0.80，r3=0.95；C2為取參數(shù)r2=0.85，r3=0.95；C3為取參數(shù)r2=0.88，r3=0.95；均采用直接初始化方法。由表5知，雖本文算法初始化時(shí)間遠(yuǎn)多于傳統(tǒng)方法，但總計(jì)算時(shí)間遠(yuǎn)少于傳統(tǒng)方法、計(jì)算精度遠(yuǎn)高于傳統(tǒng)方法，表明本文算法計(jì)算效率高，優(yōu)越性明顯。

表5 不同方法計(jì)算時(shí)間對(duì)比

5 結(jié) 論

(1) 本文基于均勻七次B樣條基函數(shù)，通過對(duì)任意局部時(shí)間域位移、速度、加速度插值構(gòu)造，推導(dǎo)出逐步積分法逆推格式與計(jì)算流程。

(2) 本文算法具有自起步、無中間計(jì)算環(huán)節(jié)特點(diǎn)。對(duì)低頻部分具有較好保持作用，可通過選取合適參數(shù)靈活控制高頻段衰減速度。

(3) 由數(shù)值試驗(yàn)給出算法絕對(duì)穩(wěn)定參數(shù)條件。由本文方法所求位移、速度、加速度計(jì)算精度均較高。

[1] Mullen R,Belytchko T. An analysis of an unconditionally stable explicit method[J]. Computers and Structures, 1983,16(6): 691-696.

[2] Gérard R, Soive A, Grolleau V. Comparative study of numerical explicit time integration algorithms[J]. Advances in Engineering Software,2005,36(4):252-265.

[3] Chang S Y. A new family of explicit methods for linear structural dynamics[J]. Computers and Structures,2010, 88(11/12):755-772.

[4] Tamma K K, Zhou X, Sha D. A theory of development and design of generalized integration operators for computational structural dynamics[J]. International Journal of Numerical Methods in Engineering,2001, 50(7): 1619-1664.

[5] Zhou X,Tamma K K. Design, analysis and synthesis of generalized single step solve and optimal algorithms for structural dynamics[J]. International Journal of Numerical Method in Engineering, 2004, 59(5): 597-668.

[6] Leontiev V A. Extension of LMS formulations for l-stable optimal integration methods withu0-v0overshoot properties in structural dynamics: the level-symmetric (LS) integration methods[J].International Journal for Numerical Methods in Engineering, 2007,71(13):1598-1632.

[7] Dokainish M A, Subbaraj K. A survey of direct time integration methods in computational structural dynamics. I. explicit methods[J]. Computers and Structures,1989, 32(6):1371-1386.

[8] Shing P B, Mahin S A. Elimination of spurious higher-mode response in pseudo-dynamic tests[J]. Earthquake Engineering and Structural Dynamics, 1987, 15(4): 425-445.

[9] Chang S Y. Improved numerical dissipation for explicit methods in pseudo-dynamic tests[J].Earthquake Engineering and Structural Dynamics,1997,26(9):917- 929.

[10] Chang S Y. Explicit pseudo-dynamic algorithm with unconditional stability[J].Journal of Engineering Mechanics, 2002, 128(9): 935-947.

[11] Piegl L, Tiller W. The Nurbs book[M]. New York: Springer-Verlag, 1997.

[12] Sevilla R, Fernandez-Mendez S, Huerta A. Nurbs-enhanced finite element method(NEFEM)[J]. Archives of Computational Methods in Engineering, 2011, 18(4): 441-484.

[13] Liu J L. Solution of dynamic response of framed structure using piecewiseBirkhoff polynomial[J]. Journal of Sound and Vibration, 2002, 251(5): 847-857.

[14] Inoue T,Sueoka A. A step-by-step integration scheme utilizing the cardinal B-splines[J]. JSME International Journal, 2002, 45(2): 433-441.

[15] Rostami S, Shojaee S, Moeinadini A. A parabolic acceleration time integration method for structural dynamics using quartic B-spline functions[J]. Applied Mathematical Modelling, 2012, 36(11): 5162-5182.

[16] 張琳,聶孟喜,仝輝. 三次和五次B樣條函數(shù)在動(dòng)力響應(yīng)分析中的應(yīng)用[J]. 清華大學(xué)學(xué)報(bào)(自然科學(xué)版)2006, 46(3): 327-330.

ZHANG Lin, NIE Meng-xi, TONG Hui. Cubic and quintic B-spline functions for dynamic response[J]. Tsinghua Univ(Sci. & Tech.), 2006,46(3):327-330.