李 彬, 唐小微(大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024)

隨著科技的不斷進(jìn)步，計(jì)算機(jī)硬件的高速發(fā)展非常有效地提高了計(jì)算效率，但是不論硬件如何提高，通過(guò)對(duì)計(jì)算方法的改進(jìn)與對(duì)新計(jì)算方法的開發(fā)來(lái)提高計(jì)算效率，依然是一件非常重要和有意義的事情。

提高計(jì)算效率的方法有很多種，眾多學(xué)者在這方面的努力一直沒有停止過(guò)。例如，王林等[1]研究基于改進(jìn)MP的稀疏表示快速算法，該算法能準(zhǔn)確提取滾動(dòng)軸承故障特征且提高了計(jì)算效率；為減少諧波合成法中功率譜矩陣分解的計(jì)算量, 祝志文等[2]提出了譜解矩陣雙軸插值算法和遞歸插值算法，兩種方法均使風(fēng)場(chǎng)模擬計(jì)算效率大幅度提高；張斌等[3]提出一種將有限元法和非線性接觸理論相結(jié)合的交叉迭代數(shù)值改進(jìn)算法，極大地提高了動(dòng)力學(xué)方程數(shù)值計(jì)算效率。張弛等[4]提出一種改進(jìn)的混合免疫算法，提高了種群多樣性和收斂性,減少了時(shí)間復(fù)雜度,提高了計(jì)算效率。

在眾多方法中，時(shí)間自適應(yīng)分析方法受到越來(lái)越多的關(guān)注。它將誤差控制在某一允許范圍內(nèi)，同時(shí)盡可能地減少計(jì)算時(shí)間，這樣就在保證計(jì)算精度的同時(shí)節(jié)省了計(jì)算時(shí)間，最終提高了計(jì)算效率。在時(shí)間自適應(yīng)研究領(lǐng)域方面，許多國(guó)內(nèi)外學(xué)者做出了努力。在巖土工程固結(jié)分析中，Sloan等[5]應(yīng)用了時(shí)間自適應(yīng)算法。Zienkiewicz等[6]提出了一種簡(jiǎn)單的時(shí)間自適應(yīng)方法，并成功應(yīng)用到動(dòng)力分析中。Zeng等[7-8]對(duì)Zienkiewicz提出的方法進(jìn)行了改進(jìn)。Tang等[9]通過(guò)控制全局誤差，將時(shí)間自適應(yīng)方法應(yīng)用到地下水流動(dòng)分析中。Kavetski等[10]采用了啟發(fā)式和誤差評(píng)估的方法進(jìn)行了時(shí)間自適應(yīng)。Kuo等[11]結(jié)合時(shí)間元素具有大時(shí)間步幅及動(dòng)量平衡具對(duì)不連續(xù)載重的平滑化的優(yōu)點(diǎn)，提出了一套加權(quán)式動(dòng)量時(shí)間元素的逐步時(shí)間積分方法。王開加等[12]在海上浮基風(fēng)電平臺(tái)繞流數(shù)值模擬分析中，采用了時(shí)間自適應(yīng)技術(shù)。毛衛(wèi)男等[13]提出了一種適應(yīng)時(shí)間的方法，并應(yīng)用于土壤凍結(jié)的水熱耦合模型中。Tang等[14-15]開發(fā)了一套完整的時(shí)間自適應(yīng)方法，并在地震液化模擬中進(jìn)行了應(yīng)用。Naveed等[16]將高階變分離散時(shí)間的自適應(yīng)時(shí)間控制應(yīng)用到了對(duì)流-擴(kuò)散反應(yīng)方程中。Vahid等[17]應(yīng)用時(shí)間適應(yīng)方法，在裂紋擴(kuò)展的相場(chǎng)公式中進(jìn)行了大規(guī)模的并行處理。Marc等[18]對(duì)非線性和時(shí)間相關(guān)的誤差估計(jì)和自適應(yīng)問題作了一個(gè)簡(jiǎn)單的概括，并指出了目前的時(shí)間自適應(yīng)方法為后驗(yàn)式方法。

在結(jié)構(gòu)動(dòng)力時(shí)程分析中，本文針對(duì)由時(shí)間步長(zhǎng)引起的誤差，建立了相對(duì)誤差與時(shí)間步長(zhǎng)之間的關(guān)系式。并將設(shè)定的目標(biāo)誤差限代入此關(guān)系式中，反算出時(shí)間步長(zhǎng)。這樣，得到的每一步時(shí)間步長(zhǎng)所對(duì)應(yīng)的相對(duì)誤差與給定目標(biāo)誤差限相等，由此保證了計(jì)算精度。同時(shí)，極大的節(jié)省計(jì)算時(shí)間，有效的提高了計(jì)算效率。

1 自適應(yīng)原理與自適應(yīng)時(shí)間步長(zhǎng)確定方法

1.1 原理

時(shí)間步長(zhǎng)固定的時(shí)間離散動(dòng)力分析方法中，一般來(lái)說(shuō)，時(shí)間步長(zhǎng)越小，所產(chǎn)生的誤差越小，計(jì)算精度越高。固定的時(shí)間步長(zhǎng)在計(jì)算過(guò)程中產(chǎn)生的誤差有大有小，而一個(gè)算法所能達(dá)到的計(jì)算精度是由計(jì)算結(jié)果最差情況，也就是產(chǎn)生的最大誤差所決定的。當(dāng)取某一固定時(shí)間步長(zhǎng)進(jìn)行計(jì)算時(shí)，每一時(shí)間步所產(chǎn)生的誤差不一定相同，誤差中較小部分對(duì)應(yīng)的固定時(shí)間步長(zhǎng)無(wú)疑占用了較多的時(shí)間資源，浪費(fèi)了時(shí)間。本文研究的方法出發(fā)點(diǎn)就是通過(guò)增加誤差中較小部分對(duì)應(yīng)的時(shí)間步長(zhǎng)，使得增加后的時(shí)間步長(zhǎng)所產(chǎn)生的誤差剛好達(dá)到所有固定步長(zhǎng)產(chǎn)生誤差中的最大誤差，這樣，就能在保證計(jì)算精度的同時(shí)，大大的節(jié)省計(jì)算時(shí)間。

在計(jì)算過(guò)程中，先設(shè)定某一相對(duì)誤差限，這個(gè)相對(duì)誤差限就對(duì)應(yīng)于所有固定步長(zhǎng)產(chǎn)生誤差中的那個(gè)最大誤差。理論上，由該相對(duì)誤差限計(jì)算得到的每一個(gè)時(shí)間步長(zhǎng)所對(duì)應(yīng)的相對(duì)誤差是一致的(與設(shè)定的相對(duì)誤差限相等)，這也使得計(jì)算精度得到了保證。如果時(shí)間步長(zhǎng)再增加一點(diǎn)點(diǎn)，所對(duì)應(yīng)的相對(duì)誤差都會(huì)超過(guò)相對(duì)誤差限。所以，計(jì)算得到的每一個(gè)時(shí)間步長(zhǎng)已經(jīng)達(dá)到最大。理論上，該方法在保證計(jì)算精度的同時(shí)，最大化地節(jié)省了計(jì)算時(shí)間，提高了計(jì)算效率。

1.2 自適應(yīng)時(shí)間步長(zhǎng)的確定

結(jié)構(gòu)動(dòng)力數(shù)值分析方法中，動(dòng)力控制微分方程通常具有如下形式

(1)

由Newmark-β法，上述動(dòng)力微分方程tn+1時(shí)刻位移近似解ui+1(t)可以表達(dá)為

(2)

(3)

將式(3)代入(2)得

(4)

Ο(Δt3)=

(5)

(6)

(7)

計(jì)算相對(duì)誤差時(shí)需要精確解，然而精確解通常是難以獲得的，這里采用泰勒級(jí)數(shù)的展開法獲得ti+1時(shí)刻位移的解來(lái)替代精確解進(jìn)行誤差評(píng)估，其表達(dá)式為

(8)

(9)

將式(4)、(5)、(8)和(9)代入式(7)中，整理后得

(10)

(11)

(12)

(13)

ax3+bx2+cx+d=0

(14)

根據(jù)盛金公式

A=b2-3ac,B=bc-9ad,C=c2-3bd,

Δ=B2-4AC

(15)

當(dāng)A=B=0,式 (13) 的解為

(16)

當(dāng)Δ=B2-4AC>0,式(13)的解為

(17)

Δt2,3=x2,3=

(18)

當(dāng)Δ=B2-4AC=0,式(13)的解為

(19)

(20)

當(dāng)Δ=B2-4AC<0,式(13)的解為

(21)

(22)

2 有效性驗(yàn)證

為了驗(yàn)證本文方法在提高計(jì)算效率方面的有效性，應(yīng)選取具有解析解的算例進(jìn)行驗(yàn)證。這里選取Bernoulli-Euler梁作為算例：一個(gè)等截面均質(zhì)的簡(jiǎn)支梁，初始位移與初始速度均為0，有一個(gè)集中力荷載從梁最左端向右端勻速運(yùn)動(dòng)，其數(shù)值模型如圖1所示。

圖1 移動(dòng)荷載下簡(jiǎn)支梁及節(jié)點(diǎn)分布Fig.1 Simply supported beam under moving load and node distribution

移動(dòng)荷載下Bernoulli-Euler梁運(yùn)動(dòng)微分方程[20]

(23)

(24)

將梁平均分為20個(gè)計(jì)算單元，節(jié)點(diǎn)分布如圖1所示，各節(jié)點(diǎn)的初始位移和初始速度均為0，有一個(gè)集中力荷載由梁最左端向右端勻速運(yùn)動(dòng)。計(jì)算參數(shù)均假設(shè)為無(wú)量綱參數(shù)(其中L，l分別為梁長(zhǎng)和單元長(zhǎng))

EI=100,c=0,L=10,l=0.5,

(25)

2.1 計(jì)算精度的對(duì)比

本文選取梁的11號(hào)節(jié)點(diǎn)為考察對(duì)象，采用固定時(shí)間步長(zhǎng)為0.08、0.02的Newmark-β法，在固定荷載從3號(hào)節(jié)點(diǎn)勻速移動(dòng)到19號(hào)節(jié)點(diǎn)的時(shí)間段內(nèi)，給出11號(hào)節(jié)點(diǎn)撓度隨時(shí)間變化的圖形，并與解析解進(jìn)行了對(duì)比，結(jié)果見圖2。

由固定時(shí)間步長(zhǎng)為0.08、0.02時(shí)的近似解和解析解，按照式(7)計(jì)算了11號(hào)節(jié)點(diǎn)撓度相對(duì)誤差隨時(shí)間變化的結(jié)果，見圖3。

(a) 時(shí)間步為0.08

(b) 時(shí)間步為0.02圖2 11號(hào)節(jié)點(diǎn)撓度Fig.2 The node 11 deflection

(a) 時(shí)間步為0.08

(b) 時(shí)間步為0.02圖3 Newmark-β法的相對(duì)誤差Fig.3 Relative error of Newmark-β method

從圖2、圖3中可以看出，取固定時(shí)間步長(zhǎng)為0.08時(shí)撓度的近似解和解析解相差較大，相對(duì)誤差最大為0.225。當(dāng)時(shí)間步長(zhǎng)取為0.02時(shí)，近似解與解析解已經(jīng)非常接近，相對(duì)誤差最大僅為0.005 5。這里，選擇固定步長(zhǎng)為0.02的Newmark-β法的計(jì)算結(jié)果，與本文方法的計(jì)算結(jié)果進(jìn)行對(duì)比。從圖3(b)中可以看出，步長(zhǎng)為0.02時(shí)相對(duì)誤差在0～0.005 5的范圍內(nèi)變動(dòng)，當(dāng)相對(duì)誤差取得最大時(shí)所對(duì)應(yīng)的計(jì)算精度就是Newmark-β法在固定步長(zhǎng)為0.02時(shí)所達(dá)到的計(jì)算精度。

從圖4中可以看出相對(duì)誤差都集中在0.005附近變動(dòng)，這是由于在式(12)的推導(dǎo)過(guò)程中，省略了Ο(Δt3)′、Ο(Δt3)″所致；最大的相對(duì)誤差為0.005 3，與固定步長(zhǎng)為0.02時(shí)Newmark-β法的最大相對(duì)誤差0.005 5基本保持一致。將圖3與圖4比較可以看出，時(shí)間步預(yù)判方法可以把相對(duì)誤差控制在很小的范圍內(nèi)變化，而Newmark-β法無(wú)法實(shí)現(xiàn)這一點(diǎn)。

圖4 時(shí)間步自動(dòng)調(diào)整后的相對(duì)誤差Fig.4 Relative error after time step automatic adjustment

2.2 計(jì)算時(shí)間的對(duì)比

表1 計(jì)算時(shí)間對(duì)比Tab.1 Comparison of computation time s

從上面的分析可以得到：與傳統(tǒng)的時(shí)間離散方法Newmark-β法相比，本文采用的方法在保證計(jì)算精度的同時(shí)，節(jié)省了69.21%的計(jì)算時(shí)間，顯著的提高了計(jì)算效率，并且可以把相對(duì)誤差控制在很小的范圍內(nèi)變化，具有明顯的優(yōu)越性。

3 荷載與自適應(yīng)時(shí)間步長(zhǎng)的關(guān)系

文中方法是取Newmark-β法的解做為近似解推導(dǎo)得來(lái)的，而Newmark-β方法的推導(dǎo)是基于對(duì)加速度的一種假設(shè)，加速度是外荷載引起的(內(nèi)力不計(jì))，那么就有必要對(duì)荷載與時(shí)間步長(zhǎng)的關(guān)系進(jìn)行探討。

Bernoulli-Euler梁算例在這個(gè)問題上，得出的結(jié)果過(guò)于復(fù)雜，很難直接看出荷載與時(shí)間步長(zhǎng)的變化關(guān)系。為了更好的研究這個(gè)問題，本文構(gòu)造了一個(gè)受力體的簡(jiǎn)單模型：物體A，質(zhì)量為m，在光滑水平面上由靜止開始運(yùn)動(dòng)，水平方向上僅受P的作用，如圖5所示。

圖5 簡(jiǎn)單模型Fig.5 A simple model

x=t-sint

(26)

(27)

(28)

(29)

圖6 時(shí)間步長(zhǎng)隨時(shí)間的變化Fig.6 Time step change over time

從圖6、圖7中分析可以看出，自適應(yīng)時(shí)間步長(zhǎng)的變化與荷載的變化密切相關(guān)：荷載變化越快，對(duì)應(yīng)的時(shí)間步長(zhǎng)越?。缓奢d變化越慢，對(duì)應(yīng)的時(shí)間步長(zhǎng)越大。在荷載為0的時(shí)刻，變化最快，對(duì)應(yīng)的時(shí)間步長(zhǎng)取得局部最小；在荷載取得最大值的時(shí)候，變化最慢，對(duì)應(yīng)的時(shí)間步長(zhǎng)取得局部最大。所以本文方法，可以反映自適應(yīng)時(shí)間步長(zhǎng)變化與荷載變化的關(guān)系，做到步長(zhǎng)的調(diào)整與荷載變化同步，能夠反映荷載變化快慢的真實(shí)歷程。

在傳統(tǒng)時(shí)間離散方法中，時(shí)間步長(zhǎng)是固定的，選擇不合適的時(shí)候會(huì)導(dǎo)致計(jì)算結(jié)果發(fā)散或者計(jì)算突然終止，當(dāng)縮小時(shí)間步長(zhǎng)到一定程度后，才能得到滿意的計(jì)算結(jié)果。產(chǎn)生這個(gè)問題的主要原因就在于此：在某一時(shí)間步內(nèi)，荷載變化劇烈，再以原固定步長(zhǎng)進(jìn)行計(jì)算會(huì)產(chǎn)生較大的誤差，這個(gè)較大的誤差可能會(huì)導(dǎo)致后續(xù)計(jì)算結(jié)果的發(fā)散，如果誤差過(guò)大還可能導(dǎo)致計(jì)算的突然終止。而在本文方法中，每一步時(shí)間步長(zhǎng)均是由給定的局部誤差限控制得到的，無(wú)論荷載如何變化，時(shí)間步長(zhǎng)都可以自動(dòng)調(diào)整以滿足局部誤差限的限制，計(jì)算不會(huì)中斷，計(jì)算結(jié)果始終能滿足計(jì)算精度的要求。

4 結(jié) 論

本文建立了時(shí)間步長(zhǎng)與相對(duì)誤差之間的關(guān)系，可依此關(guān)系計(jì)算滿足目標(biāo)誤差限的時(shí)間步長(zhǎng)，作為結(jié)構(gòu)動(dòng)力時(shí)程分析中當(dāng)前計(jì)算步的時(shí)間步長(zhǎng)，實(shí)現(xiàn)了時(shí)間步的自動(dòng)調(diào)整，豐富了時(shí)間自適應(yīng)理論研究；與傳統(tǒng)的時(shí)間離散方法Newmark-β法相比，本文采用的時(shí)間自適應(yīng)方法可以在保證計(jì)算精度的同時(shí)節(jié)省計(jì)算時(shí)間，顯著的提高了計(jì)算效率，并能很好的控制相對(duì)誤差，使其在很小的范圍內(nèi)變動(dòng)。

通過(guò)研究本文方法的自適應(yīng)時(shí)間步長(zhǎng)與荷載之間的關(guān)系，得出以下結(jié)論(內(nèi)力不計(jì))：荷載變化越快，對(duì)應(yīng)的時(shí)間步長(zhǎng)越??；荷載變化越慢，對(duì)應(yīng)的時(shí)間步長(zhǎng)越大。從而證實(shí)了該方法可以很好的反映自適應(yīng)時(shí)間步長(zhǎng)變化與荷載變化的關(guān)系，做到步長(zhǎng)的調(diào)整與荷載變化同步。

在傳統(tǒng)時(shí)間離散方法中，時(shí)間步長(zhǎng)是固定的，選擇不合適的時(shí)候會(huì)導(dǎo)致計(jì)算結(jié)果發(fā)散或者計(jì)算突然終止，當(dāng)縮小時(shí)間步長(zhǎng)到一定程度后，才能得到滿意的計(jì)算結(jié)果。這主要是因?yàn)椋涸谀骋粫r(shí)間步內(nèi)，荷載變化劇烈，再以原固定步長(zhǎng)進(jìn)行計(jì)算會(huì)產(chǎn)生較大的誤差，這個(gè)較大的誤差會(huì)導(dǎo)致后續(xù)計(jì)算結(jié)果的發(fā)散，如果誤差過(guò)大還可能導(dǎo)致計(jì)算的突然終止。而在本文方法中，每一步時(shí)間步長(zhǎng)均是由給定的局部誤差限控制得到的，無(wú)論荷載如何變化，時(shí)間步長(zhǎng)都可以自動(dòng)調(diào)整以滿足局部誤差限的限制，計(jì)算不會(huì)中斷，計(jì)算結(jié)果始終能滿足計(jì)算精度的要求，并且計(jì)算時(shí)間還能大大的節(jié)省，這些都說(shuō)明該方法在數(shù)值計(jì)算中具有重要意義。

參 考 文 獻(xiàn)

[1] 王林, 蔡改改, 高冠琪. 基于改進(jìn)MP的稀疏表示快速算法及其滾動(dòng)軸承故障特征提取應(yīng)用[J]. 振動(dòng)與沖擊, 2017, 36(3): 176-182.

WANG Lin, CAI Gaigai, GAO Guanqi. Fast sparse representation algorithm based on improved MP and its applications in fault feature extraction of rolling bearings[J]. Journal of Vibration and Shock, 2017, 36(3): 176-182.

[2] 祝志文, 黃炎. 大跨度橋梁脈動(dòng)風(fēng)場(chǎng)模擬的插值算法[J]. 振動(dòng)與沖擊, 2017, 36(7): 156-163.

ZHU Zhiwen, HUANG Yan. Interpolation algorithm for fluctuating wind field simulation of long-span bridges[J]. Journal of Vibration and Shock, 2017, 36(7): 156-163.

[3] 張斌, 雷曉燕, 羅雁云. 基于Newmark格式的車輛-軌道耦合迭代過(guò)程的改進(jìn)算法[J]. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2016, 47(1): 298-306.

ZHANG Bin, LEI Xiaoyan, LUO Yanyun. Improved algorithm of iterative process for vehicle-track coupled system based on Newmark formulation[J]. Journal of Central South University (Science and Technology), 2016, 47(1): 298-306.

[4] 張弛, 賈麗媛, 王加陽(yáng). 改進(jìn)的混合免疫算法在約束函數(shù)優(yōu)化中的應(yīng)用[J]. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2016, 47(6): 1940-1946.

ZHANG Chi, JIA Liyuan, WANG Jiayang. An improved hybrid immune algorithm for multimodal optimization[J]. Journal of Central South University (Science and Technology), 2016, 47(6): 1940-1946.

[5] SLOAN S W, ABBO A J. Biot consolidation analysis with automatic time stepping and error control, Part 2: applications[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1999, 23(6): 493-529.

[6] ZIENKIEWICZ O C, XIEY M. A simple error estimator and adaptive time stepping procedure for dynamic analysis[J]. Earthquake Engineering & Structural Dynamics, 1991, 20(9): 871-887.

[7] ZENG L F, WIBERG N E, LI X D. A posteriori local error estimation and adaptive time-stepping for newmark integration in dynamic analysis[J]. Earthquake Engineering & Structural Dynamics, 1992, 21(7): 555-571.

[8] WIBERG N E, LI X D. A post-processing technique and an a posteriori error estimate for the newmark method in dynamic analysis[J]. Earthquake Engineering & Structural Dynamics, 1993, 22(6): 465-489.

[9] TANG G, ALSHAWABKEH A N, MAYES M A. Automatic time stepping with global error control for groundwater flow models[J]. Journal of Hydrologic Engineering, 2008, 13(9): 803-810.

[10] KAVETSKI D, BINNING P, SLOAN S W. Adaptive backward Euler time stepping with truncation error control for numerical modelling of unsaturated fluid flow[J]. International Journal for Numerical Methods in Engineering, 2002, 53(6): 1301-1322.

[11] KUO S R, YAU J D, YANG Y B.A robust time-integration algorithm for solving nonlinear dynamic problems with large rotations and displacements[J]. International Journal of Structural Stability and Dynamics, 2012, 12(6):1250051-(1-24).

[12] 王開加, 程建生, 段金輝. 基于時(shí)間步長(zhǎng)自適應(yīng)技術(shù)的海上浮基風(fēng)電平臺(tái)繞流數(shù)值模擬分析[J]. 船舶與海洋工程, 2013(2): 52-57.

WANG Kaijia, CHENG Jiansheng, DUAN Jinhui. Numerical simulation analysis of offshore floating wind power generation platform based on time step adaptive method[J]. Naval Architecture and Ocean Engineering, 2013(2): 52-57.

[13] 毛衛(wèi)南, 劉建坤. 自適應(yīng)時(shí)間步長(zhǎng)法在土體凍結(jié)水熱耦合模型中的應(yīng)用[J]. 防災(zāi)減災(zāi)工程學(xué)報(bào), 2014, 34(4): 510-516.

MAO Weinan, LIU Jiankun. Application of adaptive time step method to water-heat coupling model of soil freezing[J]. Journal of Disaster Prevention and Mitigation Engineering, 2014, 34(4): 510-516.

[14] TANG X W, ZHANG X W, UZUOKA R. Novel adaptive time stepping method and its application to soil seismic liquefaction analysis[J]. Soil Dynamics and Earthquake Engineering, 2015, 71(4): 100-113.

[15] ZHANG X W, TANG X W, UZUOKA R. Numerical simulation of 3D liquefaction disasters using an automatic time stepping method[J]. Natural Hazards, 2015, 77(2): 1275-1287.

[16] NAVEED A, VOLKER J. Adaptive time step control for higher order variational time discretizations applied to convection-diffusion-reaction equations[J]. Comput Methods Appl Mech Engrg, 2015, 285: 83-101.

[17] VAHID Z R, SHEN Y X. Massive parallelization of the phase field formulation for crack propagation with time adaptivity[J]. Comput Methods Appl Mech Engrg, 2016, 312: 224-253.

[18] MARC L, SERGE P, KRISTOFFER G Z. Preface to the special issue on error estimation and adaptivity for nonlinear and time-dependent problems[J]. Comput Methods Appl Mech Engrg, 2015, 288: 1.

[19] BROWN M J M B. The discretization error of Newmark’s methods for numerical intergration in structural dynamics. Earthquake Engineering & Structural Dynamics, 1985, 13: 43-61.

[20] LIN Y H, TRETHEWEY M W. Finite element analysis of elastic beams subjected to moving dynamic loads[J]. Journal of Sound and Vibration, 1990, 136(2): 323-342.

[21] FRYBA L. Vibration of solids and structures under moving loads[M]. London: Thomas Telford, 1999.