與荷載同步變化的時間步自動調(diào)整方法

李 彬, 唐小微(大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024)

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

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

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

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

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

1.1 原理

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

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

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

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

(1)

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

(2)

(3)

將式(3)代入(2)得

(4)

Ο(Δt3)=

(5)

(6)

(7)

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

(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)

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

(16)

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

(17)

Δt2,3=x2,3=

(18)

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

(19)

(20)

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

(21)

(22)

2 有效性驗證

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

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

移動荷載下Bernoulli-Euler梁運動微分方程[20]

(23)

(24)

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

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

(25)

2.1 計算精度的對比

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

由固定時間步長為0.08、0.02時的近似解和解析解，按照式(7)計算了11號節(jié)點撓度相對誤差隨時間變化的結(jié)果，見圖3。

(a) 時間步為0.08

(b) 時間步為0.02圖2 11號節(jié)點撓度Fig.2 The node 11 deflection

(a) 時間步為0.08

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

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

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

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

2.2 計算時間的對比

表1 計算時間對比Tab.1 Comparison of computation time s

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

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

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

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

圖5 簡單模型Fig.5 A simple model

x=t-sint

(26)

(27)

(28)

(29)

圖6 時間步長隨時間的變化Fig.6 Time step change over time

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

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

4 結(jié) 論

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

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

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

參 考 文 獻

[1] 王林, 蔡改改, 高冠琪. 基于改進MP的稀疏表示快速算法及其滾動軸承故障特征提取應(yīng)用[J]. 振動與沖擊, 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] 祝志文, 黃炎. 大跨度橋梁脈動風場模擬的插值算法[J]. 振動與沖擊, 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格式的車輛-軌道耦合迭代過程的改進算法[J]. 中南大學學報(自然科學版), 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] 張弛, 賈麗媛, 王加陽. 改進的混合免疫算法在約束函數(shù)優(yōu)化中的應(yīng)用[J]. 中南大學學報(自然科學版), 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] 王開加, 程建生, 段金輝. 基于時間步長自適應(yīng)技術(shù)的海上浮基風電平臺繞流數(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)時間步長法在土體凍結(jié)水熱耦合模型中的應(yīng)用[J]. 防災(zāi)減災(zāi)工程學報, 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.