溫彥良,李 彬,唐小微

(1.遼寧科技大學(xué) 礦業(yè)工程學(xué)院，遼寧 鞍山 114051；2.泰山學(xué)院 機(jī)械與建筑工程學(xué)院，山東 泰安 271000；3.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024)

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

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

ATAM在應(yīng)用時(shí)，我們不可能為了一個(gè)具體的算例重新編寫程序，那樣就失去了提高計(jì)算效率的意義。目前已經(jīng)有很多成熟的數(shù)值計(jì)算程序，將ATAM嵌入到現(xiàn)有程序中，提高原程序的計(jì)算效率和計(jì)算穩(wěn)定性是我們開發(fā)這個(gè)算法的最終目的。

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

1 自適應(yīng)時(shí)間步長確定方法

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

(1)

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

(2)

(3)

(4)

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

(5)

(6)

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

(7)

(8)

由于ATAM只調(diào)整時(shí)間步長，不改變原程序的計(jì)算方法，所以自適應(yīng)后程序的收斂完全依賴于原計(jì)算方法。

當(dāng)采用Newmark-β法的解做為近似解時(shí)

(9)

兩者關(guān)系為

(10)

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

2 自適應(yīng)分析

Bernoulli-Euler梁：一個(gè)等截面均質(zhì)的簡支梁，初始位移與初始速度均為0，有一個(gè)集中力荷載從梁最左端向右端勻速運(yùn)動(dòng)，其數(shù)值模型如圖1所示。

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

(11)

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

(12)

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

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

(13)

2.1 以中心差分解為近似解的自適應(yīng)分析

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

相對誤差隨時(shí)間變化的結(jié)果，見圖3。相對誤差在0～0.005 6的范圍內(nèi)變動(dòng)，取得最大時(shí)所對應(yīng)的計(jì)算精度就是Newmark-β法在固定步長為0.02時(shí)所達(dá)到的計(jì)算精度。

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

相比自適應(yīng)前，自適應(yīng)后的近似解明顯偏離解析解，但是相對誤差基本保持一致。這是因?yàn)橄鄬φ`差與自適應(yīng)前整個(gè)計(jì)算過程中的最大相對誤差一致造成的，因此，自適應(yīng)前后的計(jì)算精度保持一致。由此可見，自適應(yīng)過程提高原程序計(jì)算效率的效果是明顯的。

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

(b)局部放大

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

圖4 自適應(yīng)后的相對誤差Fig.4 Relative error after time adaptive

(a)自適應(yīng)時(shí)間步

(b)局部放大

2.2 兩種類型的比較與分析

2.2.1 相同點(diǎn)

在提高計(jì)算效率方面，兩者都可以顯著的提高原程序的計(jì)算效率；在荷載與自適應(yīng)時(shí)間步長的關(guān)系上，兩者也是一致的，可以解決在沖擊荷載、振動(dòng)荷載、爆炸荷載或地震荷載的動(dòng)力數(shù)值計(jì)算中，因時(shí)間步長的設(shè)置不當(dāng)導(dǎo)致計(jì)算結(jié)果發(fā)散甚至計(jì)算中斷的問題，在一定程度上提高了原程序的計(jì)算穩(wěn)定性，相關(guān)內(nèi)容可參考《與荷載同步變化的時(shí)間步自動(dòng)調(diào)整方法》。下面就兩種類型在應(yīng)用時(shí)的另一個(gè)共同優(yōu)點(diǎn)做一個(gè)簡單介紹。

針對一個(gè)具體的數(shù)值計(jì)算方法，在實(shí)現(xiàn)時(shí)間自適應(yīng)后，計(jì)算效率具體提高了多少？目前的后驗(yàn)式時(shí)間自適應(yīng)方法無法解決這個(gè)問題，因?yàn)樵跊]有解析解的情況下，無法獲得準(zhǔn)確的相對誤差，計(jì)算精度就無法準(zhǔn)確比較，從而計(jì)算效率的提高程度也就無法獲得。ATAM可以解決這個(gè)問題。

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

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

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

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

2.2.2 不同點(diǎn)

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

圖7 自適應(yīng)后的相對誤差Fig.7 Relative error after time adaptive

在本文與《與荷載同步變化的時(shí)間步自動(dòng)調(diào)整方法》中，類型A與B之間存在如下的對應(yīng)關(guān)系

(14)

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

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

關(guān)于兩者的關(guān)系，在《與荷載同步變化的時(shí)間步自動(dòng)調(diào)整方法》中已經(jīng)得出了相關(guān)結(jié)論，這里補(bǔ)充下由兩者的關(guān)系得出的另一個(gè)適用條件。

根據(jù)《與荷載同步變化的時(shí)間步自動(dòng)調(diào)整方法》中的結(jié)論，可以得到自適應(yīng)時(shí)間步長隨時(shí)間的變化，如圖8所示。

將圖8中的時(shí)域擴(kuò)展，得到結(jié)果如圖9所示。

如圖9所示，在位移較小的時(shí)候，荷載對時(shí)間步長的影響占主導(dǎo)地位，隨著位移的增大，位移對時(shí)間步長的影響比重會(huì)逐漸增大。位移對時(shí)間步長有影響是因?yàn)樵诠降耐茖?dǎo)中選取了相對誤差做為衡量標(biāo)準(zhǔn)，在相對誤差保持不變的情況下，位移越大則近似解偏離解析解越遠(yuǎn)，則時(shí)間步長越大。

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

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

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

4 結(jié) 論

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

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

類型A可以設(shè)定預(yù)期獲得的計(jì)算精度，且計(jì)算結(jié)果與預(yù)期相符；類型B不可以。本文中，由于類型A、B存在具體的數(shù)值對應(yīng)關(guān)系，所以類型B設(shè)定的相對誤差限是自適應(yīng)后實(shí)際得的1.333倍，但這樣的對應(yīng)關(guān)系只限于此。