關于ATAM應用時注意的幾個問題

溫彥良,李 彬,唐小微

(1.遼寧科技大學 礦業工程學院，遼寧 鞍山 114051；2.泰山學院 機械與建筑工程學院，山東 泰安 271000；3.大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024)

如何有效地提高計算效率，一直是一件非常重要和有意義的事情，眾多學者在這方面的努力一直沒有停止過。例如，王林等[2]研究基于改進MP的稀疏表示快速算法，該算法能準確提取滾動軸承故障特征且提高了計算效率；為減少諧波合成法中功率譜矩陣分解的計算量,祝志文等[3]提出了譜解矩陣雙軸插值算法和遞歸插值算法，兩種方法均使風場模擬計算效率大幅度提高。

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

ATAM在應用時，我們不可能為了一個具體的算例重新編寫程序，那樣就失去了提高計算效率的意義。目前已經有很多成熟的數值計算程序，將ATAM嵌入到現有程序中，提高原程序的計算效率和計算穩定性是我們開發這個算法的最終目的。

在《與荷載同步變化的時間步自動調整方法》公式的推導過程中，選擇了Newmark-β法的解做為近似解。當采用中心差分法的解做為近似解時[15]，會得到與前者相似的結果，但略有不同。本文將比較兩者的不同，從而得出ATAM在實際應用中需要注意的問題。

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

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

(1)

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

(2)

(3)

(4)

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

(5)

(6)

將式(2)、(5)和(6)代入式(4)中，整理后得

(7)

(8)

由于ATAM只調整時間步長，不改變原程序的計算方法，所以自適應后程序的收斂完全依賴于原計算方法。

當采用Newmark-β法的解做為近似解時

(9)

兩者關系為

(10)

兩者在時間步長確定公式上的細小差別，在實際應用中卻代表著兩種不同的類型。類型A：原程序在時域離散時采用的方法，與自適應時間步長推導過程中采用近似解的方法完全相同，包括相關的參數設置也相同。例如在《與荷載同步變化的時間步自動調整方法》中，原程序采用Newmark-β法對時域進行了離散，自適應時間步的確定采用了Newmark-β法解做為近似解，且控制參數相同β=0.25。類型B：方法不相同。例如在自適應時間步長的確定中采用了中心差分解為近似解，而原程序采用Newmark-β法。在實際應用中，我們遇到的情況更多是后者。為了方便比較，本文選用了與《與荷載同步變化的時間步自動調整方法》中同一算例的同一結點進行了分析。

2 自適應分析

Bernoulli-Euler梁：一個等截面均質的簡支梁，初始位移與初始速度均為0，有一個集中力荷載從梁最左端向右端勻速運動，其數值模型如圖1所示。

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

(11)

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

(12)

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

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

(13)

2.1 以中心差分解為近似解的自適應分析

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

相對誤差隨時間變化的結果，見圖3。相對誤差在0～0.005 6的范圍內變動，取得最大時所對應的計算精度就是Newmark-β法在固定步長為0.02時所達到的計算精度。

對比自適應前后的兩個局部放大圖(圖2(b)與圖5(b))，可以看出自適應前，近似解等時間步間距不規則的分布在解析解附近。自適應后，近似解與解析解的相對誤差基本保持一致，時間步間距不相同。自適應后的點明顯少于自適應前的點，而點的多少代表計算次數的大小，點越少計算次數越小，計算時間就越少。

相比自適應前，自適應后的近似解明顯偏離解析解，但是相對誤差基本保持一致。這是因為相對誤差與自適應前整個計算過程中的最大相對誤差一致造成的，因此，自適應前后的計算精度保持一致。由此可見，自適應過程提高原程序計算效率的效果是明顯的。

(a)時間步為0.02

(b)局部放大

圖3 Newmark-β法的相對誤差Fig.3 Relative error of Newmark-β method

圖4 自適應后的相對誤差Fig.4 Relative error after time adaptive

(a)自適應時間步

(b)局部放大

2.2 兩種類型的比較與分析

2.2.1 相同點

在提高計算效率方面，兩者都可以顯著的提高原程序的計算效率；在荷載與自適應時間步長的關系上，兩者也是一致的，可以解決在沖擊荷載、振動荷載、爆炸荷載或地震荷載的動力數值計算中，因時間步長的設置不當導致計算結果發散甚至計算中斷的問題，在一定程度上提高了原程序的計算穩定性，相關內容可參考《與荷載同步變化的時間步自動調整方法》。下面就兩種類型在應用時的另一個共同優點做一個簡單介紹。

針對一個具體的數值計算方法，在實現時間自適應后，計算效率具體提高了多少？目前的后驗式時間自適應方法無法解決這個問題，因為在沒有解析解的情況下，無法獲得準確的相對誤差，計算精度就無法準確比較，從而計算效率的提高程度也就無法獲得。ATAM可以解決這個問題。

在圖4的自適應過程中，最大自適應步長為0.394 9，最小自適應步長0.020 9(與0.02基本一致)。以最小步長做為自適應前原程序的固定時間步長，相對誤差變化如圖6所示，最大相對誤差0.006 1。

圖6 Newmark-β法的相對誤差Fig.6 Relative error of Newmark-β method

根據前面得出的結果，可從兩方面驗證該方法的可行性。①最小自適應時間步長為0.020 9與固定時間步長0.02基本一致；②圖6的最大相對誤差0.006 1與圖4的0.005 5基本一致。這兩點都可說明自適應前后的計算精度基本保持一致，在此基礎上比較兩者的計算時間，就可實現自適應前后計算精度提高程度的比較。

自適應的過程中，最大時間步長與最小步長之間跨度非常大。在這里最大是最小的18.89倍，如果是沖擊荷載、振動荷載、爆炸荷載或是地震荷載，相信這個跨度會大很多。鑒于此，在程序的實現過程中，應對最大自適應步長加以限制，由于ATAM的推導過程中對時間步長要求Δt<<1，所以建議最大時間步長不超過0.4，所有超過0.4的時間步長都按0.4計算。此建議不影響自適應前后計算效率的比較。

2.2.2 不同點

首先，兩者代表應用時的兩種類型A、B，這在第一部分的自適應時間步長確定方法中已經提到了。由于中心差分法自身的缺陷，很少有數值計算程序選用它來進行時域離散，所以在實際應用中建議采用Newmark-β法的解做為近似解的式(9)來確定自適應時間步長，當遇到原程序采用Newmark-β法對時域進行離散時按A類處理，其他情況按B類處理。

圖7 自適應后的相對誤差Fig.7 Relative error after time adaptive

在本文與《與荷載同步變化的時間步自動調整方法》中，類型A與B之間存在如下的對應關系

(14)

即：類型B設定的相對誤差限是自適應后實際得到的1.333倍。如果ATAM選取Newmark-β法的解做為近似解的式(9)來確定自適應時間步長，而原程序中時域離散的方法是其他方法，這樣的對應關系就難以確定，具體關系需要具體問題具體分析，所以這樣的對應關系只限于此。

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

關于兩者的關系，在《與荷載同步變化的時間步自動調整方法》中已經得出了相關結論，這里補充下由兩者的關系得出的另一個適用條件。

根據《與荷載同步變化的時間步自動調整方法》中的結論，可以得到自適應時間步長隨時間的變化，如圖8所示。

將圖8中的時域擴展，得到結果如圖9所示。

如圖9所示，在位移較小的時候，荷載對時間步長的影響占主導地位，隨著位移的增大，位移對時間步長的影響比重會逐漸增大。位移對時間步長有影響是因為在公式的推導中選取了相對誤差做為衡量標準，在相對誤差保持不變的情況下，位移越大則近似解偏離解析解越遠，則時間步長越大。

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

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

由于ATAM的推導過程中要求時間步長Δt<<1，根據前面的結論，位移越大則時間步長越大，所以ATAM在處理大變形問題時，會受到自適應時間步長最大不超過0.4的限制，使得提高計算效率的程度受到影響。

4 結 論

在實際應用中，ATAM分為兩類情況：類型A與類型B。由于中心差分法自身的缺陷，很少有數值計算程序選用它來進行時域離散，所以在實際應用中建議采用Newmark-β法的解做為近似解的式(9)來確定自適應時間步長，當遇到原程序采用Newmark-β法對時域進行離散時按A類處理，其他情況按B類處理。

兩種類型在計算效率的提高、荷載與自適應時間步長的關系上得出的結論是一致的。兩種類型都可以對自適應前后計算效率的提高程度進行準確比較：以自適應過程中的最小時間步長做為自適應前原程序的固定時間步長，這樣就保證了自適應前后計算精度保持一致，在此基礎上比較兩者的計算時間就可以獲得計算效率的提高程度。

類型A可以設定預期獲得的計算精度，且計算結果與預期相符；類型B不可以。本文中，由于類型A、B存在具體的數值對應關系，所以類型B設定的相對誤差限是自適應后實際得的1.333倍，但這樣的對應關系只限于此。