【JY】結(jié)構(gòu)動(dòng)力學(xué)之顯隱式

建源之光 - 減隔震

2021年7月20日 15:00

瀏覽：2935 收藏：8

因你精彩

即刻關(guān)注

“通俗講動(dòng)力數(shù)值算法”

隨著科學(xué)研究和工程技術(shù)的不斷發(fā)展，出現(xiàn)了諸如航天飛機(jī)、空間站、海洋石油鉆井平臺(tái)等大型或超大型的復(fù)雜結(jié)構(gòu)，它們不但自由度高，而且還含有非線性、隨機(jī)載荷和復(fù)雜的邊界條件等多種因素。

因此，到20世紀(jì)末開始，結(jié)構(gòu)動(dòng)力學(xué)發(fā)展到一個(gè)新的高度，對(duì)復(fù)雜結(jié)構(gòu)的分析需要借助高階有限元模型進(jìn)行大量計(jì)算。

更由于一些細(xì)觀和微觀結(jié)構(gòu)分析的提出，大型和超大型的計(jì)算是不可避免的。這些問題由于計(jì)算規(guī)模巨大，求解復(fù)雜，在傳統(tǒng)的串行機(jī)上無法得到滿意的解答。因此改變計(jì)算機(jī)系統(tǒng)結(jié)構(gòu)，變串行處理為并行處理，從理論上講，其速度的提高是無止境的。利用并行計(jì)算機(jī)研究和開發(fā)相應(yīng)的并行算法有可能解決串行機(jī)無法解決的存儲(chǔ)量、計(jì)算時(shí)間和效率等多方面的問題。

在談?dòng)?jì)算效率及并串行之前(期待下期講解)，這期先談?wù)撓陆Y(jié)構(gòu)動(dòng)力學(xué)的顯式隱式計(jì)算。

【結(jié)構(gòu)動(dòng)力學(xué)研究對(duì)象】

結(jié)構(gòu)動(dòng)力學(xué)問題通常有兩類研究對(duì)象：

第一類是在運(yùn)動(dòng)狀態(tài)下工作的機(jī)械或結(jié)構(gòu)，例如高速運(yùn)行的車輛、飛行器等，它們承受著本身慣性及與周圍介質(zhì)或結(jié)構(gòu)相互作用的動(dòng)力載荷。

第二類是承受動(dòng)力載荷作用的工程結(jié)構(gòu)，例如建于地面的高層建筑和廠房，石化廠的反應(yīng)塔和管道，它們可能承受強(qiáng)風(fēng)、波浪、地震以及沖擊等各種動(dòng)力載荷的作用。

【結(jié)構(gòu)動(dòng)力學(xué)研究內(nèi)容】

結(jié)構(gòu)動(dòng)力學(xué)的研究內(nèi)容主要有兩個(gè)方面：

一是結(jié)構(gòu)的動(dòng)力響應(yīng)問題，分析結(jié)構(gòu)的動(dòng)力響應(yīng)特性，以計(jì)算結(jié)構(gòu)振動(dòng)時(shí)的動(dòng)應(yīng)力和動(dòng)位移的大小及其變化規(guī)律。

在看本文之前可以看下歷史文章：

【JY】結(jié)構(gòu)動(dòng)力學(xué)初步-單質(zhì)點(diǎn)結(jié)構(gòu)的瞬態(tài)動(dòng)力學(xué)分析

【JY】淺析各動(dòng)力求解算法及其算法數(shù)值阻尼(人工阻尼)

二是模態(tài)(振型)分析問題，尋求（包括模態(tài)計(jì)算、模態(tài)識(shí)別）結(jié)構(gòu)的固有頻率和主振型，從而了解結(jié)構(gòu)固有的振動(dòng)特性，以便更好地利用或減小振動(dòng)。

在看本文之前可以看下歷史文章：

【JY】主成分分析與振型分解

【JY】振型求解之子空間迭代

對(duì)于多自由度體系，一般可采用振型疊加法，但當(dāng)體系存在強(qiáng)非線性行為時(shí)，振型疊加法已不再適用（弱非線性行為可通過FNA，將非線性行為當(dāng)成外力施加，詳情可看：【JY】基于Ramberg-Osgood本構(gòu)模型的雙線性計(jì)算分析）。

因此目前對(duì)結(jié)構(gòu)動(dòng)力響應(yīng)分析常用的有效方法是數(shù)值積分法。數(shù)值積分法是將振動(dòng)平衡方程式中的時(shí)間分割成許多間隔，每個(gè)時(shí)間間隔都非常小以保證計(jì)算精度。針對(duì)每個(gè)時(shí)間間隔點(diǎn)計(jì)算位移、速度及加速度等，利用已經(jīng)求得的當(dāng)前步的分析結(jié)果作為已知條件，通過一定的計(jì)算方法或假定求得未知的下一步的分析結(jié)果，逐步求得結(jié)構(gòu)在地震作用下的響應(yīng)結(jié)果。

特別是結(jié)構(gòu)的恢復(fù)力特性隨結(jié)構(gòu)反應(yīng)的大小而在不斷地變化（材料非線性），因此在每步的分析中必須根據(jù)結(jié)構(gòu)反應(yīng)狀態(tài)確定當(dāng)前的結(jié)構(gòu)恢復(fù)力特性，進(jìn)行下一步計(jì)算。

直接積分法針對(duì)離散時(shí)間點(diǎn)上的值進(jìn)行計(jì)算，十分符合計(jì)算機(jī)存儲(chǔ)的特點(diǎn)，體系的運(yùn)動(dòng)微分方程也不一定要求在全部時(shí)間上都滿足，而僅要求在離散的時(shí)間點(diǎn)上滿足即可。根據(jù)在當(dāng)前時(shí)刻的反應(yīng)值確定下一時(shí)刻反應(yīng)值方法的不同。

因此在函數(shù)表達(dá)式的不同分別形成了不同的數(shù)值積分法，常見的數(shù)值積分法主要包括Runge-Kutta法、Newmark-β法、Wilson-θ法、Duhamel積分法、HHT-α法、中心差分法、分段解析法、精細(xì)積分法等等（看下圖）。

（詳情可下載學(xué)習(xí)：【JY】SDOF計(jì)算教學(xué)軟件開發(fā)應(yīng)用分享）

那么哪種算法好呢？

其實(shí)判斷一種直接積分方法的優(yōu)劣，其標(biāo)準(zhǔn)是多方面的：

a.收斂性

當(dāng)計(jì)算的時(shí)間步長趨于無窮小時(shí)，所得的數(shù)值解是否收斂于精確解。

b.計(jì)算精度

如果計(jì)算截?cái)嗾`差是相對(duì)于時(shí)間步長N次方的小量，那么就稱該方法具有N階精度。

c.穩(wěn)定性

隨計(jì)算時(shí)間步數(shù)的增大，如果數(shù)值解沒有遠(yuǎn)離精確解，則該方法是穩(wěn)定的。

d.計(jì)算效率

計(jì)算耗時(shí)的多少可直接反映出一種計(jì)算方法的計(jì)算效率。

因此，對(duì)于不同問題下，選擇不同算法是存在優(yōu)劣之分的，而常見的在實(shí)際的動(dòng)力彈塑性分析中還可以根據(jù)是否需要聯(lián)立求解耦聯(lián)方程組，且是否需要在每一步的求解中進(jìn)行剛度矩陣求逆，將數(shù)值積分法分為隱式分析方法和顯式分析方法。

【隱式分析方法】

該方法需迭代求解耦聯(lián)的方程組，計(jì)算工作量大，增加的工作量至少與自由度的平方成正比，如Newmark-β法、HHT-α法，但是積分步長可以取得較大。

聯(lián)立求解耦聯(lián)方程組和迭代都需要求解大型的線性方程組，這一過程需要占用相當(dāng)數(shù)量的計(jì)算資源、磁盤空間和內(nèi)存。該算法中的增量步可以比較大，至少可以比顯式算法大得多，但是實(shí)際運(yùn)算中上要受到迭代次數(shù)及非線性程度的限制，需要取一個(gè)合理值，而且對(duì)于高度非線性問題無法保證收斂。

【顯式分析方法】

該方法直接求解耦聯(lián)的方程組，計(jì)算工作量較小，工作量至少與自由度數(shù)成正比，如中心差分法，但是要求步長很小。

顯式算法基于動(dòng)力學(xué)方程，采用動(dòng)力學(xué)方程的一些差分格式（如廣泛使用的中心差分法、線性加速度法等）對(duì)時(shí)間進(jìn)行差分，不需要進(jìn)行剛度陣、質(zhì)量陣、阻尼陣的組，不需要進(jìn)行矩陣分割或矩陣求解等求解大型方程組，其右端項(xiàng)的形成只需根據(jù)每個(gè)單元對(duì)有效荷載向量的貢獻(xiàn)累加而成，這樣整個(gè)計(jì)算基本上在單元一級(jí)水平上進(jìn)行，計(jì)算量與節(jié)點(diǎn)數(shù)成正比，因此只需很小的高速存儲(chǔ)區(qū)，尤其當(dāng)一系列單元的剛度陣、質(zhì)量陣、阻尼陣相同時(shí)，就不需重復(fù)計(jì)算，效率更高。

另外，顯式算法不需要進(jìn)行平衡迭代，計(jì)算速度快，時(shí)間步長只要取得足夠小（典型的顯式算法的時(shí)間步長是隱式算法的1%~1%），一般不存在收斂性問題。因此需要的內(nèi)存也比隱式算法要少，并且數(shù)值計(jì)算過程可以很容易進(jìn)行并行計(jì)算，程序編制也相對(duì)簡單。但顯式算法要求質(zhì)量矩陣為對(duì)角矩陣，而且只有在單元級(jí)計(jì)算盡可能少時(shí)速度優(yōu)勢才能發(fā)揮，因而往往采用減縮積分方法。顯式計(jì)算的不足是精度不高，必須設(shè)定非常小的時(shí)間步求解以保證穩(wěn)定狀態(tài)，過大和過小的時(shí)間步往往導(dǎo)致求解時(shí)間非常漫長。特別適用于求解需要分成許多的時(shí)間增量來達(dá)到高精度的高速動(dòng)力學(xué)時(shí)間，諸如沖擊、碰撞和爆破等高度非線性問題。

【數(shù)形結(jié)合：顯式算法VS隱式算法】

以下面這個(gè)方程為例，講清楚顯式、隱式的差別：

為了方便閱讀，提前給個(gè)結(jié)論：

a. 顯式求解法：顯式算法的步進(jìn)具有“預(yù)測”的能力，僅需要當(dāng)前步的所有參數(shù)可以預(yù)測下一步的所有參數(shù)。

b. 隱式求解法：隱式算法的步進(jìn)具有自我矯正的能力，利用下一步的導(dǎo)數(shù)矯正當(dāng)前步誤差，進(jìn)行實(shí)時(shí)步步矯正。

1、對(duì)于上述方程的顯式方法的化解求解如下：

上述方程中可得到當(dāng)|1-h|>1時(shí)(即h>2)，整個(gè)方程將開始發(fā)散，需要保證h<2，整個(gè)方程才可收斂。求解函數(shù)繪制如下：

（即工程師們在采用顯式進(jìn)行計(jì)算彈塑性分析時(shí)，模型或步長不當(dāng)時(shí)，計(jì)算會(huì)發(fā)生的發(fā)散同理。）

2、對(duì)于上述方程的隱式方法的化解求解如下：

可以看到隱式計(jì)算需要用到下一步的x對(duì)y進(jìn)行矯正計(jì)算，且僅需要保證h>0 即可保證方程精度穩(wěn)定且收斂。求解函數(shù)繪制如下：

綜上所述，顯式算法具有計(jì)算存儲(chǔ)量小的優(yōu)點(diǎn)，在計(jì)算分析步內(nèi)，不需要對(duì)剛度進(jìn)行重新組裝迭代，但計(jì)算步長需要取得極小才可讓計(jì)算不發(fā)散。

時(shí)間步長的建議可取以下建議：(其中E是彈性模量，L是單元長度，Cd是材料波速)。

隱式算法具有自動(dòng)矯正的能力，且步長可以較大，但剛度矩陣計(jì)算在每一步內(nèi)需要進(jìn)行迭代重組裝，且在迭代過程，剛度重組的過程較容易計(jì)算不收斂，并需要占用計(jì)算存儲(chǔ)大量資源。

最近筆者也開發(fā)小型的計(jì)算內(nèi)核，受曲哲老師的極速牛頓法的啟發(fā)（文獻(xiàn)DOI：https//doi.org/10.1002/nme.6456）這里提供一種思路給大家：

隱式 + PCG + 極速牛頓

解析：隱式算法使得步長拉長，PCG（預(yù)處理共軛梯度法）使得求解超大規(guī)模的稀疏矩陣變得簡便快速，極速牛頓法進(jìn)在迭代時(shí)最多兩步，盡可能避免不收斂問題，集百家之所長的快準(zhǔn)狠算法~

（圖片來源：哲設(shè)計(jì)）

最后，其實(shí)顯式和隱式求解的本質(zhì)在于采用不同的物理學(xué)平衡方程，因此在不同的物理學(xué)問題也有不同的表現(xiàn)。順便說下，在大型通用有限元中，如Abaqus/Standard用的隱式分析方法的為HHT-α法（改進(jìn)Newmark-β法），如Abaqus/Explicit、LS-DYNA顯式分析方法為中心差分法。