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

(1. 重慶大學 煤礦災害動力學與控制國家重點實驗室，重慶 400044； 2. 仲愷農業(yè)工程學院，廣州 510225)

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

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

1 均勻七次B樣條

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

(1)

(2)

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

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

(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樣條基函數對t的導數階次。

2 基本原理及算法實現

2.1 動力平衡方程遞推格式

單自由度系統結構動力學微分方程可寫為

(11)

取0=t0<…

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

(12)

式中：

(14)

為便于推導，令

(15)

(16)

(17)

由式(12)有

(18)

式(18)可寫為矩陣形式

(19)

式中：

(20)

(21)

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

(22)

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

(23)

令

(24)

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

(25)

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

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

(26)

將上式運算簡化為

(27)

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

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)為計算x(3)(t0)的直接初始化方法。

(30)

(31)

2.3 計算步驟

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

3 穩(wěn)定性分析

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

圖2 本文算法不同參數r對應譜半徑

本文算法與傳統Newmark方法、Wilson 方法譜半徑見圖3。由圖3曲線對比知，本文算法對低頻段具有更好的保持作用，能有效濾掉高頻分量?？紤]工程應用中有限元網格計算的高頻分量往往不準確，且對結構相應貢獻較小，因而可通過增大參數r2,r3能快速有效濾掉高頻分量。

4 數值算例與分析

4.1 單自由度強迫振動

表1 不同方法所求位移及相對誤差

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

表2 不同方法所求速度相對誤差

表3 不同方法所求加速度相對誤差

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

表4 不同方法計算時間對比

4.2 簡支梁彎曲振動

圖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, 阻尼系數c=0，橫向載荷q(x,t)=F0sin(ω0t)δ(x-L/2)，F0=1 kN,ω0=4 Hz。該載荷作用下，簡支梁振動位移理論解為

圖4 簡支梁彎曲振動示意圖

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

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

圖5 給定時刻不同方法數值解

圖6 中點位置不同方法數值解

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

表5 不同方法計算時間對比

5 結 論

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

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

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

[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樣條函數在動力響應分析中的應用[J]. 清華大學學報(自然科學版)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.