與荷載同步變化的時間步自動調整方法

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

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

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

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

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

1 自適應原理與自適應時間步長確定方法

1.1 原理

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

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

1.2 自適應時間步長的確定

結構動力數值分析方法中，動力控制微分方程通常具有如下形式

(1)

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

(2)

(3)

將式(3)代入(2)得

(4)

Ο(Δt3)=

(5)

(6)

(7)

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

(8)

(9)

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

(10)

(11)

(12)

(13)

ax3+bx2+cx+d=0

(14)

根據盛金公式

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 有效性驗證

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

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

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

(23)

(24)

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

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

(25)

2.1 計算精度的對比

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

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

(a) 時間步為0.08

(b) 時間步為0.02圖2 11號節點撓度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時，近似解與解析解已經非常接近，相對誤差最大僅為0.005 5。這里，選擇固定步長為0.02的Newmark-β法的計算結果，與本文方法的計算結果進行對比。從圖3(b)中可以看出，步長為0.02時相對誤差在0～0.005 5的范圍內變動，當相對誤差取得最大時所對應的計算精度就是Newmark-β法在固定步長為0.02時所達到的計算精度。

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

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

2.2 計算時間的對比

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

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

3 荷載與自適應時間步長的關系

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

Bernoulli-Euler梁算例在這個問題上，得出的結果過于復雜，很難直接看出荷載與時間步長的變化關系。為了更好的研究這個問題，本文構造了一個受力體的簡單模型：物體A，質量為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中分析可以看出，自適應時間步長的變化與荷載的變化密切相關：荷載變化越快，對應的時間步長越小；荷載變化越慢，對應的時間步長越大。在荷載為0的時刻，變化最快，對應的時間步長取得局部最小；在荷載取得最大值的時候，變化最慢，對應的時間步長取得局部最大。所以本文方法，可以反映自適應時間步長變化與荷載變化的關系，做到步長的調整與荷載變化同步，能夠反映荷載變化快慢的真實歷程。

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

4 結 論

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

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

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

參 考 文 獻

[1] 王林, 蔡改改, 高冠琪. 基于改進MP的稀疏表示快速算法及其滾動軸承故障特征提取應用[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] 張弛, 賈麗媛, 王加陽. 改進的混合免疫算法在約束函數優化中的應用[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] 王開加, 程建生, 段金輝. 基于時間步長自適應技術的海上浮基風電平臺繞流數值模擬分析[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] 毛衛南, 劉建坤. 自適應時間步長法在土體凍結水熱耦合模型中的應用[J]. 防災減災工程學報, 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.