【文章內(nèi)容簡介】 增量 U 如下： U ＝ Ut － U? ＝ tV?? ＋ ? ?? ?? ?2221 tA t ???? ?? ? /2 (224) 由于全量和增量拉格朗日列式法可以統(tǒng)一寫成： M At ＋ K? U ＝ Rt ＋ stR (225) 將 (224)代入 (225)得： M At ＋ K? tV?? ＋ K? ? ?? ?? ?2221 tA t ???? ?? ? /2＝ Rt ＋ stR (226) 整理得： [M ＋ ? ? ?2t? K? ]At =Rt K? tV?? － K? ? ? A??21? ? ?2t? /2＋ stR (227) 令等效質(zhì)量矩陣 M~ 和等效節(jié)點力向量 Rt~ 為 M~ =M ＋ ? ? ?2t? K? (228) Rt~ =Rt K? tV?? － K? ? ? A??21? ? ?2t? /2 (229) 則有： M~ At ＝ Rt~ ＋ stR (230) 由于 stR 可以表示為 At 的函數(shù)，即 stR ＝ AQtA ＋ AR? ，則： （ M~ － AQ ） At ＝ Rt~ ＋ AR? (231) 當加速度從 (231)得出后便可從方程式 (222)和 (223)中計算出位移和速度。 應當指出，無論實現(xiàn)求出加速 度，還是先求出位移，都會遇到系數(shù)矩陣為變量的情況。這是因為在非線性問題中，剛度矩陣是系統(tǒng)位移和變形的函數(shù)， 12 需要用迭代的方法加以解決。 最終，收斂的條件有不同的表現(xiàn)形式。以能量為判據(jù)的收斂條件表示為： )1()1()1()()( ffeIfIf LRULRU ???? ? (232) 式中 )(IfR 表示方程右邊在第 I 次迭代時的量， )1(fL 表示方程左邊在第 I 次迭代時的量， e? 表示給 定的收斂控制值， U? 表示相鄰兩次迭代所得的位移差： U? ＝ )(IU － )1(?IU (233) 收斂條件也可以用不平衡力為判據(jù)來表示： )1()1()()( ffeIfIf LRLR ??? ? (234) 隱式算法的收斂問題是一個十分重要的問題，但很復雜。對于一個給定的問題，應用 什么樣的收斂判據(jù)，給多大的收斂控制值，都比較難確定。對于沖壓成形這種復雜的非線性問題來說，隱式算法常常無能為力。這就是顯式算法在成形加載過程計算中優(yōu)于隱式算法的一個原因。 板殼理論 在沖壓成形分析中使用的最為普遍及有效的有兩種殼單元，即 HughesLiu和 BelytschkoTsay（ BT）單元?？紤]到研究模擬內(nèi)容是以 BT 殼單元為模擬單元，因為 BT 殼單元是 DYNAFORM 缺省的單元公式，它采用面內(nèi)單點積分，對大變形問題是最穩(wěn)定和有效的公式，另外，采用 Corotational 應力更新，單元坐標 系統(tǒng)置于單元中心，基于平面單元假定，所以對于翹曲的幾何體不適用。因此僅就 BT 殼單元有限元理論加以研究。 HughesLiu 殼單元是退化的殼單元公式，因而可以適應任意復雜外形，并具有較高的計算精度。但是，由于單元公式比較復雜，計算量較大，在求解大型復雜薄板成形問題時需要較長的計算時間，而 BT 殼單元則是用于顯式計算的一種更為有效的單元。 1. 隨動坐標系 對于每一個單元，定義一個局部的隨動坐標系。該坐標系隨單元一起運動，其基底（ 234。1， 234。2， 234。3）由下式定義： 234。3＝ r31r 42/|r31r 42| (235) S1=r21(r21234。 3)234。3 (236) 234。1=S1/|S1|， 234。2=234。3234。 1 (237) 注意， 234。3 總是垂直于參考面，當單元的四個節(jié)點共面時， 234。3 的方向與“厚度”方向 l 重合，當殼體變形時，實際的“厚度”方向與 234。3 之間存在一個角度，在BT 單元中，假定該角度很小，即下式成立： 13 |234。3l－ 1|δ (238) 對大多數(shù)彈塑性工程計算問題，取 σ≈ 102 2. 應變位移關(guān)系 BT 殼單元采用 Mindlin 板殼理論，殼體中任意一點的速度可寫為： V= Vm－ ?234。3 Q (239) 其中 Vm 為參考面上的速度， Q 為角速度， ?為沿厚度方向的距離。 3．有限元方程 由標準的有限元離散，可得有限元方程為： intffvM ext ??? (240) 其中 M 為質(zhì)量矩陣， extf 為外力矢量， intf 為內(nèi)力矢量，單元內(nèi)力為： ???????eeemff (241) 由虛功原理，可得： ????? ev TeITeIeITeIeITeIeITeI dvdfvmQfvmQ ?????? ?????? (242) 其中 Td? ＝ ? ?2313122211 ?2,?2,?2,?,? ddddd (243) ?? ＝ ? ?2313122211 ?,?,?,?,? ????? T (244) eIf ＝ ? ?III fff 321 , T (245) eIm ＝ ? ?III mmm 321 , T (246) 對于四節(jié)點雙線性殼體單元，并采用單點積分，即在一個單元內(nèi)，彎矩 eIm和膜力 eIf 只計算一次，對于彈塑性材料，沿厚度方向可取 5～ 7 個幾分點。 根據(jù)有限元插值，有： ? ??????????????????????IIIImmmxxxNxxx321321,?? （ I＝ 1， 4） (247) 其中（ Ix1 ， Ix2 ， Ix3 ）為節(jié)點坐標， IN 為插值函數(shù)： 1N ＝（ 1－ ? ）（ 1－ ? ） /4 2N ＝（ 1＋ ? ）（ 1－ ? ） /4 3N ＝（ 1＋ ? ）（ 1＋ ? ） /4 4N ＝（ 1－ ? ）（ 1＋ ? ） /4 （ 248) 參考面上的速度和角速度由相同的插值函數(shù)進行插值： ? ? IIm VNV ?? ,? (249) ? ? II QNQ ?? ,? (250) 14 由上述公式，可以求得單元中心處 ? ?0,0 ?? ?? 的速度應變的具體表達式： IIII QBxVBd 2131111 ??? ?? IIII QBxVBd 1232222 ??? ?? )(????2 11223211212 IIIIIIII QBQBxVBVBd ???? IIII QNVBd 23113 ??2 ?? IIII QNVBd 13223 ??2 ?? 其中 11 ?xNB II ??? 22 ?xNB II ??? (251) 單元節(jié)點內(nèi)力的具體表達式為： )( 122111 fBfBAf III ?? )( 121222 fBfBAf III ?? )( 2321313 fBfBkAf III ?? )4( 23121221 fkmBmBAm III ??? )4( 13122112 fkmBmBAm III ??? 03 ?Im (252) 其中 A 為單元面積， k 為剪切修正系數(shù)： ?? 3?? xdf ???? ? ??? 33 ??? xdxm ???? ? ? ?3,2,1, ??? (253) 接觸問題 一般情況下 ,板料的成形實際上是在沖頭運動的過程中完成的。隨著沖頭的運動，沖頭和模具表面因和板料接觸而對板料施加的作用力是板料得以成形的動力。在接觸過程中，板料的變形和接觸邊界的摩擦作用使得部分邊界條件隨加載過程而變化，從而導致了邊界條件的非線性。正確地處理邊界的接觸和 摩擦是得到可信分析結(jié)果的一個關(guān)鍵因素。 接觸力計算 接觸力計算基本算法有兩種，一種是罰函數(shù)法，另一種是拉格朗日乘子法。 在罰函數(shù)法中位于一個接觸面上的接觸點允許穿透與之相接觸的另一個接觸面，接觸力地大小與穿透量成正比，即 : Sfn ??? (254) 15 式中 α 是罰因子， S 是接觸點的法向穿透量，負號表示接觸力與穿透方向相反。罰因子的取值過小會影響精度，過大會降低計算的穩(wěn)定性，在實際計算時要認真選取。 在拉格朗日乘子 法中，接觸力是作為附加自由度來考慮的，其泛函形式除了包含有通常的能量部分外還附加了拉格朗日乘子項： ? ? ? ?DQuFuKuuu TTT 021, ???? ??? (255) 式中 u 是結(jié)點位移向量， K 為剛度矩陣， F 為結(jié)點力向量， ? 是拉格朗日乘子向量， ? ?DQuD 0?? 為接觸點的穿透量向量。 對 能量泛函式變分，建立有限元方程： ???????????????????? DFuQ QK T 00 ? (256) 求解方程即可得到結(jié)點位移和拉格朗日乘子，拉格朗日乘子的分量即為接觸點處的法向接觸力。 拉格朗日乘子法是在能量泛函極小的意義上滿足接觸點互不穿透的邊界條件，它增加了系統(tǒng)的自由度，需要采用迭代算法來求解方程，一般適用于靜態(tài)隱式算法。在顯式算法中，一般采用罰函數(shù)法。這種方法即考慮了接觸力，又不增加系統(tǒng)的自由度，計算效率較高。 摩擦的處理 摩擦力的準確計算對板料成形分 析十分重要 ,目前常用的摩擦定律還是基于經(jīng)典的庫侖摩擦定律 ,只是為了數(shù)值計算的穩(wěn)定性作了一些修正。 由經(jīng)典庫侖摩擦定律可知，當兩接觸物體間的切向摩擦力 ft 小于臨界值 ftc時，兩接觸面間的相對滑移 ut =0，而當 ft = ftc 時，相對滑移是不定的，需由外界載荷和約束條件確定。按照經(jīng)典摩擦定律計算的摩擦力為： tff nt ???? (257) 式中 ? 是摩擦系數(shù)， fn 是接觸點的法向接觸力， t? 是相對滑動方向上的切向單位向量 : rrvvt ??? (258) rv? 是相對滑動速度向量。 由于在板料成形分析中，某些局部的相對速度很小或相對速度方向發(fā)生突變、接觸狀態(tài)由粘著到滑動或相反的變化將導致按式（ 257）計算得到摩擦力大小和方向突變，從而引起計算的不穩(wěn)定。目前一般通過引入光順函數(shù)來修正庫侖摩擦定律，可用的光順函數(shù)有反正切函數(shù)和雙曲正 切函數(shù)，可得到以下修正的庫侖摩擦定律： 16 tvvar c t gffcrnt ?? ?????????? ?? 2 (259) tvvtanhffcrnt ?? ?????????? ? (260) 式中 cv 是一個給定的相對滑動速度，它的大小決定了修正的摩擦模型和原模型的相近程度。太大的 cv 導致有效摩擦力數(shù)值的降低，但使迭代相對容易收斂，而太小的 cv 雖然能夠較好模擬摩擦力的突變，但使求解的穩(wěn)定性下降。 修正的摩擦模型式（ 259）和（ 260）近似如圖 所示。經(jīng)典的庫侖摩擦定律是從最初適用的剛體一般化到變形體，是在遵循“切向力到達某一臨界值時，接觸表面才會在局部產(chǎn)生滑移”這一假設(shè)的前提下應用的。盡管這一假設(shè)在一定的情況下有效，但嚴格地說它是不成立的。實驗發(fā)現(xiàn)，只要有切向力存在，兩接觸表面就會產(chǎn)生滑移，據(jù)此，研究人員又提出了一些非線性的摩擦定律。但是這些摩擦模型有的過于復雜，有的一些系數(shù)很難通過實驗得到，使用的較少。 本章小結(jié) 本 章主要從有限元理論上對板料成形進行描述，通過數(shù)學方法，將板料成形中的三大非線性問題數(shù)值逼近和簡化。在這些理論研究上，應用于板料成形的有限元理論很多，在計算機編程中，得到廣泛引用的主要有兩大方面，即顯式動力算法和隱式算法，前者主要用于沖擊大變形分析，后者主要