13-2.靜力隱式-顯式和動力顯式有限元列式

彈塑性大變形有限元方法胡平§大變形彈塑性本構(gòu)方程最簡單的材料大變形彈塑性本構(gòu)方程是不考慮加載歷史的形變理論的全量本構(gòu)方程。這類本構(gòu)方程基于加載歷史是比例加載的基本假定。對于許多金屬沖壓成形問題，比例加載的基本假定是近似正確的，因此，基于這類理論的所謂onestepinversealgorithm，由于其高效高速的計算效率，近年來被廣泛應(yīng)用于汽車概念設(shè)計階段的沖壓工藝性快速校核。關(guān)于全量理論的One-stepinversealgorithm的理論知識將另行介紹（參見理論手冊）。然而，更精細(xì)和準(zhǔn)確的成形問題數(shù)值分析應(yīng)該是能夠考慮加載歷史的增量本構(gòu)理論。為了反映與加載歷史的相關(guān)性，大變形彈塑性問題需要采用速率型（增量型）的本構(gòu)方程。首先，需要在大變形條件下，按照符合客觀性的要求，建立增量型本構(gòu)方程。由（7.26）式可知，變形率張量［d］是客觀張量，但由（7.36）式可知Cauchy應(yīng)力張量的物質(zhì)導(dǎo)數(shù)匚］不是客觀張量，所以若用變形率張量和Cauchy應(yīng)力張量來建立本構(gòu)方程,則將不滿足本構(gòu)方程的客觀性條件（7.42）式。為此，需要另外定義Cauchy應(yīng)力張量的導(dǎo)數(shù)1．Cauchy應(yīng)力張量的Jaumann導(dǎo)數(shù)（率）而Cauchy應(yīng)力的Jaumann導(dǎo)數(shù)定義為(1)(1)推導(dǎo)后得到&—lim丄j推導(dǎo)后得到&—lim丄jdt—>0dt-aiJl+dtt2）3）4）CV=Q-QO-QO

jjikkjkjki其矩陣形式為第一Piola應(yīng)力（名義應(yīng)力）張量的本構(gòu)導(dǎo)數(shù)t(t)=G+Ol-Qlijijijkkkjik用Cauchy應(yīng)力Jaumann導(dǎo)數(shù)表示的第一Piola應(yīng)力的本構(gòu)導(dǎo)數(shù)為t(t)=Qv—Qd—Qd+Ql+Qlijijikkjkjkiikjkijkk其矩陣形式為[t(t)]=[QV]—LqJldJ-LdLqJ+LqJI/片+【Q」tLlJr第二Piola應(yīng)力張量的本構(gòu)導(dǎo)數(shù)仿照第一Piola應(yīng)力本構(gòu)導(dǎo)數(shù)的推導(dǎo)，有T（t）=Q+Ql-Ql-Ql（、ijijijkkkjikikjk（6）把（7.179）式代入（7.182）式，便得用Cauchy應(yīng)力Jaumann導(dǎo)數(shù)表示的第二Piola應(yīng)力的本構(gòu)導(dǎo)數(shù)為T（t）=QV-Qd-Qd+Ql”ijijikkjkjkiijkk（7）其矩陣形式為88）)4）單元應(yīng)力本構(gòu)導(dǎo)數(shù)單元應(yīng)力本構(gòu)導(dǎo)數(shù)為{V}=rD1[b]{v}eepL5）單元速率剛度方程由于fcvSddv二d}tbvldvijij=￡V}e)J也A[D][BIe'ePJQ5(dd)dv=J{§d}tb]tgvkzzldvije]dv{v}eLepL丿e4o)ikkid

e(5v}e)J[B]t[q]Ib]dv{v}eLdL丿e41)Jq5(ll)dv二J{5l}r[o]{l}dv2ijkjkiV{v}e42)=(5v}e)ff[B]t[q{v}e42)vvv丿e其中「[q][o][o]][q]=[o][q][0]

v[o][o][q]5353)[o]=ooo111231ooo122223ooo312333ii工f(k)8v(k)=v}e)44)45)k=1各式代入fp8vdvii()8

?f]T[Nf)T卜v({?(f8{2?p)8?f3ef~p8vda=iiaa(40)式，并考慮到CV}e)的任意性，便得單元速率剛度方程,?8vb)f[血I~p|da48)其中eprC]+其中eprC]+[k]){v}e=0o49)505152k]=f[Bb「D1[B]dv0LepL[k]e=6A[o][B]-[B}[o][b])dvovvvLdLe{}{}p=f[N]rpdveep>=f[N]r<丿ao6）總體速率剛度方程仿照線性彈性有限元法，可由單元速率剛度方程組集成總體速率剛度方程，即([k方程，即([k]+[k]){v}=0b{}54)其中{v}='V1V112{v}='V1V112

[k]=FGa])b][a]

0nV1???VnVnVn3123e=1[k]=蘭Ga])[kMa]bb丿p55)56575859)由于剛度矩陣［k］、［k］與時間t無關(guān)，所以求解時可把（54）式寫成以節(jié)0b點位移增量為基本未知量的增量形式，即Qk］+［k］）{au}={ap}（60）0b顯然，需用增量法求解，而把第i-1增量步所得［b］和［b］作為第i增量步Vd的初應(yīng)力，按（7.236）式求得該增量步的節(jié)點位移增量{au}，把i-1步所得單元形狀作為計算單元剛度矩陣和單元體、面力節(jié)點載荷增量的體積和面積。這樣求解，就是把第i增量步起始時刻（t時刻）的構(gòu)形作為參考構(gòu)形，把第i增量步終了時刻C+At時刻）的構(gòu)形作為現(xiàn)時構(gòu)形。所以，看起來

似乎用的是Lagrange描述方法，只不過是逐步修改的Lagrange方法，即U.L法。大變形彈塑性動力顯式有限元法動力顯式有限元算法的理論基礎(chǔ)是連續(xù)介質(zhì)振動理論中的動力學(xué)方程和中心差分格式算法。首先，我們先簡單地回顧一下最基本的單自由度動力學(xué)方程及其解法。kmAAAA/AAAA/P,x,v,a圖7－7.1線性彈簧系統(tǒng)1單自由度系統(tǒng)對于只有一個自由度的線性彈簧系統(tǒng)，如圖7-7.1所示，該系統(tǒng)的平衡方程為：ma+cv+kx=Pnnnn(/?23/丿其中，n代表時間增量步的第n瞬時。

⑴?無阻尼系統(tǒng)的自由振動如果(7.237)式中的阻尼系數(shù)等于零，而且右端載荷項也為零，則(7.237)式可以寫成無阻尼的自由振動方程，人k令w2二mma+kx=0nn(7.238)方程(7.238)改寫為X+w2x=0nn(7.239)其解為xn=Acoswt+Asinwt12(7.240)式中，A］和A2為待定常數(shù)，由系統(tǒng)的初始條件確定。設(shè)t=0時，x=x0，x=x0，則由(7.240)式得￡二x°，A2二/式(7.240)還可以變換為相角形式的簡諧函數(shù)，(7.241)x=Acos(①t—申)(7.241)n式中，A=：A2+A2,申二arctg212A1簡諧函數(shù)的性質(zhì)可用圖(7-7.2)來表示。圖7-7.2簡諧振蕩圖7-7.2簡諧振蕩（2）．有阻尼系統(tǒng)的自由振動如果（7.237）式中僅右端載荷項為零，則（7.237）式可以寫成有阻尼的自由振動方程，(7.242)ma+cv+kx=0(7.242)nnn引入臨界阻尼系數(shù)C,CC=2jmk=2m?(7.243)⑴欠阻尼（C<Cc）（圖C=2jmk=2m?(7.243)⑴欠阻尼（C<Cc）（圖7-7.3）方程（7.242）的解為x=Ae-ecos（\：⑵臨界阻尼（C=Cc）（圖7-7.4）方程（7.242）的解為(7.244)x=(x+(xx)t)e—wtn000⑶過阻尼(C>Cc)(圖7-7.5)方程(7.242)的解為(7.245)x=Aeg/+Aes2tn12(7.246)其中，a=xo+(G+￡2-l)?xo,a=12叭匕2-12—x—(^+l,'^2—1)?x

=002叭U2—1S1,2=(%土山2—1)?0.8-八O.S-\1\…0.-4-1\1\/\f'、0.2-\\/\r\Ob)11-1；11/1/1\1./1\=0.2--□.斗--O.S-\jU:2、動力學(xué)平衡方程的中心差分求解格式為了求解動力平衡方程（7.237）,ma+cv+kx=Pnnnn將時間歷程分成無數(shù)個離散的時間點，每兩個時間點之間稱為時間步長，這樣，通過中心差分法就可以計算出每個時間點處的位移、速度、加速度、應(yīng)力、應(yīng)變…等物理量。中心差分方法的具體推導(dǎo)過程如下（見圖8.6）：設(shè)t時刻的狀態(tài)為n，t及t時刻之前的物理量已知，定義:Atn-1Atn(n-1)(n)(n+1)(n-訂(n+訂n-1t—+At2n-1tt+iAtt+At2nnt-At設(shè)：設(shè)：圖7-7.6離散的時間點t-At：n-1狀態(tài)；n-1t—+At:n—+狀態(tài);2n-12t++At:n++狀態(tài)；2n2t+At：n+1狀態(tài)。n再設(shè)：卩=巴’并將（7-237）式n—1ma+cv+kx=Pnnnn中的速度蔦和加速度"P分別用差分格式改寫為：1-v+vn+丄1+Dn一丄22(V—v)(1+卩)Atn—1叱n—2t+At（n+1狀態(tài)）時刻的總位移可以累加得到nx=x+vAtn+1nn+亠n2將（7.247）、（8.248）兩式帶入（7.237）式，并令尼）得7.247)7.248)7.249)c=a-m(比例阻2m(1+0)Atn-1(v-v)+am(n+十20v1+0n+2P-kxnn7.250)由上式解出：對于多自由度系統(tǒng)，(1+0)At+n-1—(2+a0At)mnn

n-12-aAtv2+a0Atn-1只要將(7?247)~(7.251)式中的v,a,x,m,k,P分別-(P-kx)7.251)用矢量或矩陣表示，{v}=n7.252){a}=n2(TT7

{x}二n+1l\m\({v}-{v})At、n+1n-Jn-1{x}+{v}Atnn+in27.253)7.254)(1+0)At二{p}-HUnn({v}-{v}n+12)+an-127.255)如果［m］可以寫成對角矩陸則上式又可以寫成分量的形式，從而得到速度的顯式表達(dá)格式，2-aAtvi=n-1vin+十2+a0Atn-12n-12其中，i表示第i個自由度。(1+0)Atn-1+?(Pi—kxi)(2+a0At)iminn(7.256)n-1為了保證系統(tǒng)計算的穩(wěn)定性，中心差分法對允許使用的時間步長有一個限制（參見振動分析的數(shù)值方法），即，At<Atcr7.257)式中，Atcr表示離散系統(tǒng)的臨界時間步長。對于單自由度系統(tǒng)而言，臨界時間步長按下式計算7.258)At7.258)crAtcr對于多自由度系統(tǒng)，Atcr7.259)其中，3表示系統(tǒng)的最高頻率分量。m動力顯式算法(1)動力顯式有限元方程fCyP,u,uxxx\P,u,uxxx\P,yIu,uyyTxy卓TyEdx左dxdxc+ydyy6y取出一個微小的正平行六面體，它在x和y方向的尺寸分別為dx,dy,它在Z方向的尺寸取單位長度。在X軸方向，有

do(◎+xdx)dyx1-odyx1+(t+牛dy)dxx1-tdxx1+(P-pu一cU)xdxdyx1=0yxTOC\o"1-5"\h\zxdxxyxdyyx等式兩端同除以dxdy得dodt匸+匹+P-pU一cu=0dxdyxxx同理，沿y軸方向有，dodt匕+江+P-pu-cu=0dydxyyy采用動力顯式算法，綜合上面兩個方程，得板料的運動學(xué)微分方程為，do..(7.260)j+p-Pu-cU=0dxiii(7.260)j式中，P是材料的質(zhì)量密度，c是阻尼系數(shù)，u和u分別是材料內(nèi)任一點的ii速度和加速度，p.是作用在該點上的體積力，o是該點處的Cauchy應(yīng)力。iij根據(jù)散度定理j§udVq6uddV―J..Jvijij式中，5u.是虛速度，壯是對應(yīng)于Cauchy應(yīng)力o的虛應(yīng)變速率。iijij把物體離散成m個單元，對于任一單元，有°個結(jié)點，取其形函數(shù)為N°，單元內(nèi)任意點的位移分量u、速度分量u和加速度分量u分別為iii

u=u=Nauai<u=Nauaiiu=Nauaii(7.262)由幾何方程可得<=Baua(7.263)ijji式中，ua、ua和u，a分別為結(jié)點a的位移分量、速度分量和加速度分iii量，Baj為第五章給出的應(yīng)變張量(具體求解時，采用線性化處理，可以按j照小變形的應(yīng)變矩陣考慮，但要進(jìn)行構(gòu)形更新)。將(7.262)式和(7.263)式代入(7.261)式，可得JpNauaN卩§u卩dV+fcNauaN卩6u卩dVTOC\o"1-5"\h\ziiiieVe=JpN卩6u卩dV+JqN卩6u卩dT-JcB卩§u卩dViiiiijjiVTVeee(7.264)式中6u0是結(jié)點p的虛速度，上式還可以寫成iJpNaN卩dV〃a+JcNaN卩dVuaiiVeVe(7.265)(7.266)=JpN0dV+JqN0dT—JcB0dV(7.265)(7.266)iiijjTOC\o"1-5"\h\zVTVeee寫成矩陣形式為fpNtNdVU+JcNtNdVU=JNtpdV+JNtqdT—JBtgdVVTV.eee將單元方程集合，即得整體有限元方程

工(JpNtNdV)u+工dcNtNdV)U(7.267)(7.268)二口NtpdV+SfNTqdT-SfBtqdV(7.267)(7.268)vrv.eee上式可簡寫為Mu+Cu=P-F式中，M是一致質(zhì)量矩陣，M=UpNTNdV；Vec是阻尼矩陣，CtNdV；Vep是結(jié)點外力向量，P=UNTPdV+UNTqdr；VreF是結(jié)點內(nèi)力向量，F(xiàn)二BrdVoVe通常，動力顯式積分算法采用集中質(zhì)量矩陣，既M是一對角矩陣，并取C=aM。貝（7.268）式表示的聯(lián)立方程組變成了（結(jié)點數(shù)x結(jié)點自由度數(shù)）個相互獨立的方程咖"F咖"Fi(7.269)這里對i不求和，i=1?總自由度數(shù)（2）臨界時間步長的計算由于中心差分格式的算法是條件穩(wěn)定的，為了保證系統(tǒng)計算的穩(wěn)定性，對時間增量步長At的大小必須加以限制。穩(wěn)定性條件通常由系統(tǒng)的最高頻率?決定，滿足穩(wěn)定性條件的時間增量步長為max

2,At<-(J+g2—g)(7.270)max其中g(shù)是具有最高頻率的模型的臨界阻尼比。與工程直覺相反，阻尼的引入實際上降低了系統(tǒng)的臨界穩(wěn)定性條件。系統(tǒng)的最高頻率由網(wǎng)格中最大的單元膨脹模式?jīng)Q定。滿足穩(wěn)定性條件的時間增量步長可以由膨脹波沿網(wǎng)格中任意單元的最小穿越時間近似得到(7.271)其中丫=0.6?0.8，c為膨脹波在材料中的傳播速度7.272)C=(上)7.272)E，Le為第n狀態(tài)(t時刻)單元e的名義長度。nn穩(wěn)定性條件可以保證在一個時間增量步內(nèi)，擾動只傳播網(wǎng)格中的一個單元。如果系統(tǒng)只包括一種材料，則滿足穩(wěn)定性條件的時間增量步長與網(wǎng)格中最小的單元尺寸成正比；如果系統(tǒng)劃分的單元網(wǎng)格尺寸比較均勻，但包括多種不同材料，則具有最高膨脹波速的材料中網(wǎng)格尺寸最小的單元決定系統(tǒng)的穩(wěn)定時間步長。對于一個簡單的桁架單元，在團聚質(zhì)量矩陣的情況下，穩(wěn)定性準(zhǔn)則給出一個臨界時間步長：At<L，式中，c為材料聲速，l為單元長度，At表cAt示膨脹波穿越長度為l的單元所需要的時間。這就是所謂的Courant-Friedrichs-Lewy(CFL)穩(wěn)定性條件。對于三角形單元和四邊形板單元來說，臨界時間步長的選取依賴于單

元名義長度的確定，一般按照圖8.7的原則來確定單元的名義長度lcrit,對于高階單元來說，臨界時間步長遠(yuǎn)比低階單元的臨界時間步長來的小，這一事實使得高階單元對于顯式積分算法相當(dāng)不合適。rl、瓷2鼻l=rl、瓷2鼻l=min(l,l)crit12crit圖8.7單元名義長度12雖然上面給出的穩(wěn)定性準(zhǔn)則嚴(yán)格來說是對線性系統(tǒng)而言的，但對于非線性問題也給出了有用的穩(wěn)定性估計。對線性問題時間步長的80%?90%縮小，對于大多數(shù)非線性問題保持其系統(tǒng)穩(wěn)定性是足夠的。然而，十分重要I的是：在整個計算過程中，要不斷的檢査能量的平衡問題。任何總能量的增加或損失（5%或更多）都將導(dǎo)致失穩(wěn)。Belytschko指出：常增量時間步不能保持解的穩(wěn)定性，即使系統(tǒng)的最高頻率°max不斷減小。3.接觸問題的動力平衡方程對于考慮阻尼的接觸問題的動力平衡方程又可表達(dá)為Mu+Cu+