顯式求解方法和隱式求解方法對(duì)比

采用有限元方法開(kāi)展結(jié)構(gòu)的動(dòng)力學(xué)分析最終歸結(jié)為求解離散后的常微分方程組MU+CU+KU=R。在時(shí)域內(nèi)求解該方程最常用的方法是直接積分法，而t又根據(jù)求解過(guò)程中是否需要迭代求解線(xiàn)性方程組，將直接積分法分為隱式積分方法和顯式積分方法兩類(lèi)。隱式積分法認(rèn)為t+At時(shí)刻系統(tǒng)的狀態(tài)不僅與t時(shí)刻狀態(tài)有關(guān)，且與t+At時(shí)刻某些量有關(guān)。因此隱式算法是根據(jù)tn及tn廠(chǎng)??時(shí)刻體系的物理量值建立關(guān)于以t+1時(shí)刻物理量為未知量的線(xiàn)性方程組，通過(guò)求解方程組確定tn+1時(shí)刻的物理量（常用的方法有線(xiàn)性加速度法、常平均加速度法、Newmark方法、Wilson-0法、Houbolt方法等）。而顯式積分法認(rèn)為t+At時(shí)刻系統(tǒng)的狀態(tài)僅與t時(shí)刻狀態(tài)有關(guān)可，因此可由tn及tn1...時(shí)刻體系的物理量值直接外推tn+1時(shí)刻物理量值（如中心差分法），不需要求解線(xiàn)性方程組，實(shí)現(xiàn)了時(shí)間離散的解耦。兩種算法的比較：（1）隱式算法隱式算法基于虛功原理，要迭代計(jì)算。隱式算法在每一增量步內(nèi)都需要對(duì)靜態(tài)平衡方程進(jìn)行迭代求解，并且每次迭代都需要求解大型的線(xiàn)性方程組，這一過(guò)程需要占用相當(dāng)數(shù)量的計(jì)算資源、磁盤(pán)空間和內(nèi)存。理論上在這個(gè)算法中的增量步可以很大，但是實(shí)際運(yùn)算中上要受到接觸以及摩擦等條件的限制。隨著單元數(shù)目的增加，計(jì)算時(shí)間幾乎呈平方次增加。由于需要矩陣求逆以及精確積分，對(duì)內(nèi)存要求很高。隱式算法的不利方面就是收斂問(wèn)題不容易解決，且在開(kāi)始起皺失穩(wěn)時(shí)，在分叉點(diǎn)處剛度矩陣出現(xiàn)奇異。（2）顯式算法顯示算法基于動(dòng)力學(xué)方程，無(wú)需迭代，包括動(dòng)態(tài)顯式和靜態(tài)顯式算法。動(dòng)態(tài)顯式算法采用動(dòng)力學(xué)方程的中心差分格式，不用直接求解切線(xiàn)剛度，不需要進(jìn)行平衡迭代，計(jì)算速度快，也不存在收斂控制問(wèn)題。該算法需要的內(nèi)存也比隱式算法要少。數(shù)值計(jì)算過(guò)程可以很容易地進(jìn)行并行計(jì)算，程序編制也相對(duì)簡(jiǎn)單。它也有一些不利方面。顯式算法要求質(zhì)量矩陣為對(duì)角矩陣，而且只有在單元級(jí)計(jì)算盡可能少時(shí)速度優(yōu)勢(shì)才能發(fā)揮,因而往往采用減縮積分方法，容易激發(fā)沙漏模式，影響應(yīng)力和應(yīng)變的計(jì)算精度。靜態(tài)顯式法基于率形式的平衡方程組與Euler前插公式，不需要迭代求解。由于平衡方程式僅在率形式上得到滿(mǎn)足，所以得出的結(jié)果會(huì)慢慢偏離正確值。為了減少相關(guān)誤差，必須每步使用很小的增量，通常一個(gè)仿真過(guò)程需要多達(dá)幾千步。由于不需要迭代，所以這種方法穩(wěn)定性好，但效率低。以下是ABAQUS軟件中顯式和隱式常采用的方法：（摘自幫助文檔）結(jié)構(gòu)系統(tǒng)的通用運(yùn)動(dòng)學(xué)方程為：MU+CU+KU二R （1）t中心差分法（顯式）假定0，t1，t2，…，tn時(shí)刻的節(jié)點(diǎn)位移，速度與加速度均為已知，現(xiàn)求解tn（t+At）時(shí)刻的結(jié)構(gòu)響應(yīng)。中心差分法對(duì)加速度，速度的導(dǎo)數(shù)采用中心差分代替，即為：1U=（U-2U+U）tAt2t—At tt+AtTOC\o"1-5"\h\zU=右（U -U） ⑵t 2/\t t+At t—At將（2）式代入（1）式后整理得到MMu 二R （3）t+At t式（3）中11M=——M+——CAt2 2At211R=R—(K———M)U—(——M———C)Utt At2 tAt2 2At t—At分別稱(chēng)為有效質(zhì)量矩陣，有效載荷矢量。R,M,C,K為結(jié)構(gòu)載荷，質(zhì)量，阻尼，剛度矩陣。求解線(xiàn)性方程組(3),即可獲得t+At時(shí)刻的節(jié)點(diǎn)位移向量U，將U代t+At t+At回幾何方程與物理方程，可得t+At時(shí)刻的單元應(yīng)力和應(yīng)變。中心差分法在求解t+At瞬時(shí)的位移U時(shí)，只需t+At時(shí)刻以前的狀態(tài)變量t+Atu和u，然后計(jì)算出有效質(zhì)量矩陣Mi，有效載荷矢量R，即可求出u，故t t—At t t+At稱(chēng)此解法為顯式算法。中心差分法，在開(kāi)始計(jì)算時(shí)，需要仔細(xì)處理。t=0時(shí)，要計(jì)算u，需要知At道U的值。因此應(yīng)該有一個(gè)起始技術(shù)，因而該算法不是自動(dòng)起步的。由于U，—At 0U，U是已知的，由t=0時(shí)的(2)式可知：00At2U=U—AU+——U—At 0 020中心差分法中時(shí)間步長(zhǎng)At的選擇涉及兩個(gè)方面的約束：數(shù)值算法的穩(wěn)定性和計(jì)算時(shí)間。中心差分法的實(shí)質(zhì)是用差分代替微分，并且對(duì)位移和加速度的導(dǎo)數(shù)采用線(xiàn)性外插，這限制了At的取值不可過(guò)大，否則結(jié)果可能失真過(guò)大?？梢宰C明：中心差分法是條件穩(wěn)定的。即當(dāng)時(shí)間步長(zhǎng)At必須小于由該問(wèn)題求解方程性質(zhì)所決定的一個(gè)時(shí)間步長(zhǎng)的臨界值。LS-DYNA中，采用“變時(shí)間步長(zhǎng)法”,即每一時(shí)刻的步長(zhǎng)At由當(dāng)前結(jié)構(gòu)的穩(wěn)定性條件來(lái)控制。具體算法為：計(jì)算每一個(gè)單元的極限時(shí)間步長(zhǎng)At，，=1,2...‘取At=min(At)為下一個(gè)時(shí)刻的時(shí)間步長(zhǎng)。eiei各種單元的At計(jì)算方法如下。e，(1)1D桿，梁?jiǎn)卧狝t=aL為桿，梁?jiǎn)卧拈L(zhǎng)度。L為桿，梁?jiǎn)卧拈L(zhǎng)度。其中a為時(shí)間步長(zhǎng)因子，系統(tǒng)默認(rèn)為0.9。材料聲速。2)2D板，殼單元AtAt=aeL—minC其中a為時(shí)間步長(zhǎng)因子，系統(tǒng)默認(rèn)為0.9°l為殼單元的最小單元邊長(zhǎng)度。 minC=E為材料的聲速。(1—V2)P

3)3D單元L3)3D單元LAt=ae其中\(zhòng)CC+CL|C|(C<0)Q=J1 0ekkkk|0(C>0)kkC和C為無(wú)量綱常數(shù)，默認(rèn)C=1.5，C=0.06.0101(對(duì)8節(jié)點(diǎn)體單元)emaxL(對(duì)4節(jié)點(diǎn)體單元)minQ+(Q2+C2)1/2L為單元等效長(zhǎng)度，V為單元體積，eA為單元最大側(cè)面積。emaxIe(1-v)c二E(1v) 為材料聲速。(1+v)(l-2v)p時(shí)間步因子a可由用戶(hù)設(shè)置，減小a相當(dāng)于減少時(shí)間步長(zhǎng)。設(shè)置時(shí)間步長(zhǎng)因子a的關(guān)鍵字為*CONTROL_TIMESTEP，控制參數(shù)為T(mén)SSFAC。另外，質(zhì)量縮放可以人為控制時(shí)間步長(zhǎng)。即調(diào)整單元密度P來(lái)改變時(shí)間步長(zhǎng)。以殼單元為例說(shuō)明質(zhì)量縮放改變時(shí)間步長(zhǎng)。At=a—min(假如a=1)eC::E\；(1-V2)P,At (1-v2)p/口— (At )2E由(——specified)2二 匸得至HP二 specifiedL E i L2(1-v2)min minNewmark法(隱式)Newmark假定在時(shí)間間隔[t,t+At]內(nèi)，加速度線(xiàn)性變化，即采用如下的加速度，速度公式：U=U+[(1-8)U+bU]Att+At t t t+AtU=U+UAt+[(丄-a)U+aU]At2 (4)t+At tt 2 t t+At式中a，8為按積分的精度和穩(wěn)定性要求可以調(diào)整的參數(shù)。根據(jù)(4)式可給出U和U用U，U，U表示的表達(dá)式，代入(1)t+At t+At t+At tt式中整理得到Ku二R (5)t+At t+At其中K二—M+—C+KaAt2 aAt111888R=R +M[——U+——U+(——-1)/7]+C[ U+(--1)/7+( -l)AtU7]t+At t+At aAt2 t aAt t 2a t aAt tat 2a t稱(chēng)之為有效剛度矩陣和有效載荷矢量。由上式可以看出求解當(dāng)前U ，t+At需要用到當(dāng)前時(shí)刻的R，因此該算法為隱式算法。t+At當(dāng)載荷歷史全部已知時(shí)，F(xiàn)為已知量，求解需要迭代實(shí)現(xiàn)。t+At可以證明，當(dāng)參數(shù)6>0.5,a>0.25(0.5+5)2時(shí)，Newmark法是無(wú)條件穩(wěn)定的，即At的大小不影響數(shù)值穩(wěn)定性。此時(shí)時(shí)間步長(zhǎng)At的選擇主要根據(jù)解得精度確定。一般，Newmark法可以比中心差分法的時(shí)間步長(zhǎng)大得多。比較兩種算法，顯式中心差分法非常適合研究波的傳播問(wèn)題，如碰撞、高速?zèng)_擊、爆炸等。顯式中心差分法的M與C矩陣是對(duì)角陣，如給定某些有限元節(jié)點(diǎn)以初始擾動(dòng)，在經(jīng)過(guò)一個(gè)時(shí)間步長(zhǎng)后，和它相關(guān)的節(jié)點(diǎn)進(jìn)入運(yùn)動(dòng)，即U中這些節(jié)點(diǎn)對(duì)應(yīng)的分量成為非零量，此特點(diǎn)正好和波的傳播特點(diǎn)相一致。另一方面，研究波傳播的過(guò)程需要微小的時(shí)間步長(zhǎng)，這也正是中心差分法的特點(diǎn)。而Newmark法更加適