1、剛塑性有限元法及其在軋制中的應用剛塑性有限元法及其在軋制中的應用軋制技術及連軋自動化國家重點實驗室1.1.學習目的和要求學習目的和要求l了解現代軋鋼生產和軋制技術的發展概況；l了解現代軋制理論研究的基本任務；l掌握剛塑性有限元的基本概念；l掌握剛塑性有限元的基本理論；l掌握剛塑性有限元的基本方法；2.2.學習的主要內容學習的主要內容l剛塑性有限元的基本概念和基本理論；l剛塑性有限元相關技術問題的處理方法；l求解軋制過程的剛塑性有限元程序。3.3.本課程的基礎和相關知識本課程的基礎和相關知識l 現代塑性加工力學現代塑性加工力學 基本方程、變分原理、有限元基礎知識；l 工程數學工程數學 矩陣分析、

2、優化方法、數值分析；l計算機基礎知識計算機基礎知識 操作系統、FORTRAN語言和FORTRAN4.0編程軟件。4.4.講課和學習方法講課和學習方法l 課堂講授課堂講授 基本概念、基本理論、基本方法 程序剖析；l 課外自學課外自學 消化理解、閱讀程序；l上機實踐上機實踐 調試程序1.1.緒緒 論論 1.1 1.1 現代軋制理論研究的發展概況現代軋制理論研究的發展概況1.1.1 1.1.1 現代軋鋼生產的發展現代軋鋼生產的發展1.1.2 1.1.2 軋制技術的發展軋制技術的發展1.1.3 1.1.3 現代軋制理論研究的基本任務現代軋制理論研究的基本任務1.2 1.2 軋制理論數值方法軋制理論數值

3、方法1.2.1 1.2.1 初等理論中的數值方法初等理論中的數值方法1.2.2 1.2.2 滑移線理論及其數值解法滑移線理論及其數值解法1.2.3 1.2.3 能量法及其數值解法能量法及其數值解法1.2.4 1.2.4 彈塑性有限元法彈塑性有限元法1.1 1.1 現代軋制理論研究的發展概況現代軋制理論研究的發展概況 20世紀世紀60年代前年代前，軋鋼生產過程手工手工操作和使用單體單體設備。 軋制理論主要解決問題軋制力、力矩、功率、寬展和前滑等參數的近似近似計算。 主要進展提出卡爾曼和奧羅萬方程，采用一些假設條件推導出軋制力和寬展等公式，逐步形成了以工程法為核心的傳統軋制理論體系。 20世紀世紀

4、60年代以后年代以后，隨著軋鋼生產軋鋼生產和軋制技術軋制技術的飛躍發展和用用戶對產品質量戶對產品質量要求的日益提高，以計算機計算機為工具，以現代數值現代數值分析方法分析方法的為特征的現代軋制理論現代軋制理論得到了迅速發展。1.1.1 1.1.1 現代軋鋼生產的發展現代軋鋼生產的發展l 20世紀世紀5070年代年代發展趨勢是大型化、高速化和連續化發展趨勢是大型化、高速化和連續化 1960年前建立的熱帶鋼軋機，輥身范圍11202490mm，年生產能力100200萬噸，帶鋼卷重614噸，最大精軋速度為1012m/s，技術進步是將AGC應用于精軋機； 20世紀6070年代，軋機向現代化技術方面發展，同

5、時連鑄技術發展成熟。大型連鑄坯、步進式加熱爐、大型化的粗軋機、7機架精軋機組、AGC、升速軋制、層流冷卻技術以及軋制過程計算機控制的全面應用。l 20世紀世紀80年代以后年代以后軋鋼生產主要向提高產品質量、降低消軋鋼生產主要向提高產品質量、降低消耗、優化軋制過程、開發新鋼材和新品種方向發展。耗、優化軋制過程、開發新鋼材和新品種方向發展。 板形、厚度及超級鋼我國軋鋼生產的發展我國軋鋼生產的發展 1957年鞍鋼第一套2800/1700mm半連續式板帶鋼軋機 到目前為止，輥身長度在1422mm以上的熱軋寬帶鋼軋機8套、薄板坯連鑄連軋帶鋼軋機10余套。 武鋼、本鋼1700mm3/4連續式熱帶鋼軋機各一

6、套 寶鋼2050mm3/4連續式熱帶鋼軋機 攀鋼1450mm半連續式熱帶鋼軋機 太鋼1549 mm半連續式熱帶鋼軋機 梅鋼1422mm全連續式熱帶鋼軋機 寶鋼1580、鞍鋼1780mm半連續式熱帶鋼軋機各一套 珠鋼1500、邯鋼1900和包鋼1750薄板坯連鑄連軋機各一套。1.1.2 1.1.2 軋制技術的發展軋制技術的發展軋鋼生產的發展促進了軋制技術的進步l 連鑄技術l 連鑄直接軋制技術(CC-DR)l 連鑄熱裝直接軋制技術(CC-HCR)l AGC、AFC、ATCl SFR及無頭軋制技術l ISP及CSP薄板坯連鑄連軋技術 1.1.3 1.1.3 現代軋制理論研究的基本任務現代軋制理論研究

7、的基本任務求解軋制變形區各種分布量，如應力場、應變場、速度場和溫度場等，為板形板厚控制和型鋼孔型設計提供理論基礎。對軋制過程中工具及工件的溫度溫度與變形變形進行綜合研究，為鋼的高精度軋制及軋機的高精度控制服務。對軋件不均勻變形及軋件頭尾不穩定變形過程的理論研究，為提高產品質量和成材率、進一步優化軋制規程服務。提高軋制過程參數的理論解析精度，建立和完善建立和完善控制軋制過程的數學模型。開展軋制過程熱力學及冶金學參數的綜合研究，對軋制過程的變形溫度、變形程度、金屬的微觀組織及產品的最終性能進行綜合模擬，實現實現根據產品使用進行鋼材成份及軋制過程的預設計。1.2.1 1.2.1 初等理論中的數值方法

8、初等理論中的數值方法 采用有限差分方法求解卡爾曼或奧羅萬方程。 基本思想：在變形區內取微元體，建立力平衡微分方程，然后在變形區內進行差分網格劃分，在已知邊界條件下，采用差分方法求解微元體上的力平衡方程。 特點：能夠定性地得出變形區中的軋制力和金屬流動規律，但計算精度有待于進一步提高。1.2.2 1.2.2 滑移線理論及其數值解法滑移線理論及其數值解法 滑移線法：把軋制過程變形區劃分為一系列由滑移線族組成的滑移線網絡，每條滑移線均為達到屈服切應力k，根據Henky應力方程可以確定變形區的應力場。 近年來，利用計算機可以形成金屬成型變形區的滑移線網絡，并計算相應的滑移線場。 特點：滑移線法只能處理

9、理想剛塑性體平面變形或軸對稱變形問題，對三維變形問題、溫度和材料性質參數分布不均問題是無能為力的。 1.2.3 1.2.3 能量法及其數值解法能量法及其數值解法 能量法的基礎是剛塑性材料的變分原理。 基本思想：給定邊界條件設定含有待定參數的運動許可速度場或靜力許可應力場建立相應能量泛函使其最小化確定待定參數得到真實的速度場由塑性力學基本關系求出變形及力能參數得到變形區內的應變場。 優點：能量法可以求解三維變形問題，直接得出變形功率、轉矩和由速度場決定的寬展、前滑; 缺點：由于不能直接得出靜水壓力，所以不能直接得出應力分布。此外，能量法也難以處理溫度、變形抗力等不均勻分布的問題。 1.2.4 1

10、.2.4 彈塑性有限元法彈塑性有限元法 彈塑性有限元法分析金屬成型時采用彈塑性材料本構關系，考慮變形的歷史相關性，在求解時需要采用增量加載，在每一個加載步中，只能有少數單元從彈性狀態進入塑性狀態，以便減小計算誤差，因此，所需計算機的容量較大、計算時間長。 優點：不僅可以求解塑性區的擴展、應力、應變分布，而且可以有效地處理卸載問題，計算殘余應力、殘余應變分布。 缺點：存在積累誤差，計算機容量較大，計算時間長。2.2.剛塑性有限元法的基本理論剛塑性有限元法的基本理論 2.1 2.1 有限元法的基本概念有限元法的基本概念 有限元法有限元法：把工件劃分成有限結點相連接的單元，以結點上的速度(位移)作為

11、未知量，利用最小能原理求解相應的方程組確定此未知量，按結點速度與單元內部應變以及單元內部應力之間的關系確定各單元的應力、應變分布。2.2 2.2 剛塑性有限元法及基本思想剛塑性有限元法及基本思想 用有限元方法分析金屬塑性成型過程時，采用剛塑性用有限元方法分析金屬塑性成型過程時，采用剛塑性材料本構模型進行求解材料本構模型進行求解就是剛塑性有限元法。 基本思想基本思想：從剛塑性材料的變分原理變分原理出發，按有限元模式把能耗率泛函能耗率泛函表示為節點速度的非線性函數，利用數學上的最優化理論最優化理論得出滿足極值條件的最優解，即使總能耗率取最小值的運動許可速度場運動許可速度場，根據塑性力學的基本關系和

12、本構方程得出應變速度場、應力場以及變形和力能參數。 2.3 2.3 剛塑性材料模型剛塑性材料模型 金屬成型過程中，材料變形的物理過程非常復雜，為了便于數學處理，必須做出一些假設，把變形中的某些過程理想化。用剛塑性有限元法分析材料變形問題時，材料滿足下列基本假設： (1) 材料均質各向同性； (2) 忽略材料的彈性變形，不計體積力與慣性力； (3) 材料的變形流動服從Levy-Mises流動理論； (4) 材料的體積不變或微可壓縮。 2.3.1 2.3.1 理想剛塑性材料模型理想剛塑性材料模型 理想剛塑性材料模型的基本假設如下： (1) 材料均質各向同性； (2) 忽略材料的彈性變形，不計體積力

13、與慣性力； (3) 材料的變形流動服從Levy-Mises流動理論； (4) 材料的體積不變； (5) 不考慮加工硬化，忽略變形抗力對變形速度的敏感性。2.3.1.1 2.3.1.1 理想剛塑性材料模型的應力應變關系理想剛塑性材料模型的應力應變關系圖2-1 理想剛塑性材料的應力應變關系2.3.1.2 2.3.1.2 理想剛塑性材料模型的特點理想剛塑性材料模型的特點l 只要等效應力達到一恒定數值，材料便發生屈服，而且材料在整個變形過程中屈服應力不再發生變化。l 采用理想剛塑性材料模型進行能量積分時，可以把等效應力做為常數提到積分號之外，從而使積分過程得到簡化。2.3.1.3 2.3.1.3 采用

14、該模型進行采用該模型進行FEMFEM求解應該注意的問題求解應該注意的問題 在軋制變形區中，由于軋件各點的溫度、變形速度和變形程度的不同，屈服應力相差很大，在整個變形區內采用理想剛塑性材料模型必然會給計算結果帶來誤差。因此，用有限元法求解時，把變形區劃分成足夠多的單元，這樣可以認為每個單元內的溫度、變形速度和變形程度相同，在每個單元內采用理想剛塑性材料模型，不同單元采用不同的屈服應力，這樣處理才能得到比較接近實際的結果。2.3.2 2.3.2 剛塑性硬化材料模型剛塑性硬化材料模型 剛塑性硬化材料基本假設如下： (1) 材料均質各向同性； (2) 忽略材料的彈性變形，不計體積力與慣性力； (3)

15、材料的變形流動服從Levy-Mises流動理論； (4) 材料的體積不變； (5) 考慮加工硬化和變形抗力對變形速度的敏感性。 2.3.2.1 2.3.2.1 剛塑性硬化材料的應力應變關系剛塑性硬化材料的應力應變關系 圖2-2 剛塑性硬化材料的應力應變關系2.3.2.2 2.3.2.2 剛塑性硬化材料的變形抗力剛塑性硬化材料的變形抗力 對于剛塑性硬化材料來說，當材料的化學成份和物理狀態一定時，通常把變形抗力表示成變形溫度、變形速度和變形程度的函數：( , , ) sf T 0nsa()nmsbc(2-1)(2-2)(2-3)2.3.2.3 2.3.2.3 Mises流動法則流動法則 理想剛塑性

16、材料和剛塑性硬化材料都假設材料是不可壓縮的，根據Mises流動法則，變形速度分量與偏差應力分量成正比，即 這種材料的變形速度場與偏差應力場一一對應，但由于靜水壓力是不確定的，所以當以速度為未知量進行求解時，不能直接求得應力場。而體積可壓縮材料模型可巧妙地解決這一問題。 yxyyzxzxzxyzxyyzzxd(2-4)2.3.3 2.3.3 剛塑性可壓縮材料模型剛塑性可壓縮材料模型 剛塑性可壓縮材料的基本假設如下： (1) 材料均質各向同性； (2) 忽略材料的彈性變形，不計體積力與慣性力； (3) 材料的變形流動服從Levy-Mises流動理論； (4) 材料的體積微可壓縮 ； (5) 考慮加

17、工硬化和變形抗力對變形速度的敏感性。(1) (1) 剛塑性可壓縮材料的屈服條件剛塑性可壓縮材料的屈服條件 剛塑性可壓縮材料與剛塑性硬化材料的主要差別是放松了體積不變條件的約束，即假設屈服與靜水壓力有關，屈服條件不僅取決于偏差應力的二次不變量，也取決于應力的一次不變量。122222222223()1 =()()()6()6xyzxyzmxyyzzxxyyzzxxyyzzxxyyzzxJJ (2-5)(2-6)剛塑性可壓縮材料的屈服條件剛塑性可壓縮材料的屈服條件 在主軸條件下，屈服條件：1222221223311()()() 2mg(2-7)122222222122211()()()6()2 39

18、xyyzzxxyyzzxmggJJ(2-8)Misees屈服曲面平面1122 , 23圖2-3 剛塑性可壓縮材料主應力空間的屈服曲面可壓縮參數可壓縮參數對對屈服條件的影響屈服條件的影響 在平面上(1 + 2 + 3 =3m=0 ),橢球體與Mises屈服圓柱相切，即剛塑性可壓縮材料與Mises屈服條件一致。在平面以外(m0)，剛塑性可壓縮材料比理想剛塑性材料容易進入屈服狀態。可壓縮參數g值越小，橢球體長軸延長，越接近于Mises屈服條件。當考慮加工硬化時，橢球體的體積將隨著加工硬化而膨脹。 數量場和矢量場數量場和矢量場l 數量場：對于空間或部分空間的任意一點M，都有一個確定的數值f(M)與之對

19、應，則稱在這個空間或部分空間上確定了一個數量場。該數量場可用數值函數f(M)來確定。例如溫度場就是一個數量場。l 矢量場：對于空間或部分空間的任意一點M，都有一個確定的矢量f(M)與之對應，則稱在這個空間或部分空間上確定了一個矢量場。該矢量場可用矢量函數f(M)來確定。例如速度場、應力場均是矢量場。 勢和勢場勢和勢場l 梯度：設函數u = f ( x, y, z )在空間或部分空間具有一階連續偏導數，則對空間或部分空間的每一個點P ( x, y, z )都可以確定一個矢量： 這個矢量就稱作函數u = f ( x, y, z )在點P ( x, y, z )的梯度。l 利用場的概念，該矢量函數在

20、空間或部分空間確定了一個矢量場，即梯度場。它是由數值函數產生的，稱數值函數為這個矢量場的勢，矢量函數則稱為有勢場或勢場。 kzfjyfixfzyxf),( grad(2) (2) 塑性勢和變形速度分量塑性勢和變形速度分量 假設剛塑性可壓縮材料的屈服函數為塑性勢 ：1222222221()()()()6()2ijxyyzzxxyyzzxmFg(2-9)()ijijijFdd根據塑性勢的定義：()ijijijFd則(2-10)dtdddt 兩邊同除以時間增量，令變形速度分量變形速度分量 利用(2-9)式和(2-10)式可直接求出應變速度分量：2(1)92(1)92(1)9222xxmyymzzmx

21、yxyyzyzzxzxgdgdgdddd(2-11) 與與 之間關系之間關系 下面以 為例推導 與 之間關系dd()ijxxFddd1222222221()()()()6()2ijxyyzzxxyzmFg1222222222()()()312()()()6()22223()33 2 ()2 ()3 xyzxmijxxyyzzxxyyzzxmxyzmxxyzmijijxgFgggFF由于2233(1)392 ()2 ()mmxmijijggFF 與與 之間關系之間關系 可得：(2-12)32 ()ijddF dd的簡化形式的簡化形式 根據單位體積的塑性變形功率 ：dxxyyzzxyxyyzyzz

22、xzxdw (2-13)把(2-11)式代入(2-13)式整理可得：212921292129xxmxyxyyymyzyzzzmzxzxgddgddgdd所以 32d (2-14)222222222222222222222222()1()92 2()3192 2()331 333xyzxyyzzxmxyxxyzxyyzzxmxyzxyyzzxmmxydgdgdgd 2222222222222222222222236()921 333()6()231 2222226()231 ()(3zxyyzzxmmxyzxyzxyyzzxmxyzxyyzzxxyyzzxmxygdgdgd 2222222)()

23、6()22 3yzzxxyyzzxmgd(3) (3) 應力應變速度關系應力應變速度關系 把(2-14)式代入(2-11)式21921392219222xxmyymvzzmvmmxyxyyzyzzxzxgdgddgdggddd 剛塑性可壓縮材料剛塑性可壓縮材料 與與 g 和和 關系關系 22333vxyzxyzmmg dgdg(2-15)當 一定時，g 越小，體積變形速度小；當g = 0時，此時剛塑性可壓縮材料的屈服條件變為Mises屈服條件，應力應變關系變為Levy-Mises流動法則，剛塑性可壓縮材料模型變為不可壓縮的理想剛塑性或剛塑性硬化材料模型。從這個意義來講，剛塑性不可壓縮材料可以看

24、作剛塑性可壓縮材料在g=0時的特例，因此更具有普遍性。 根據(2-11)和(2-14)式可得體積變形速率與可壓縮參數和靜水壓力關系 ：mV1222222221()()()6()2xyyzzxxyyzzxmgm(3) (3) 應力應變速度關系應力應變速度關系 2213922 13922 139212 39xxmvxvxxvgggggg0.5/2 3xyxyxyd如果采用張量形式，并引入克羅內克爾符號，可得剛塑性可壓縮材料的應力應變速度關系：212()39ijijijvg(2-16)(4) (4) 等效應變速度等效應變速度 把(2-16) 代入(2-13)式：212()39ijijijvg212

25、=39212 39212 3922 33xxyyzzxyxyyzyzzxzxxxvyyvzzvxyxyyzyxggg 23zxzx 12222222221()()()6()9xyyzzxxyyzzxvg(2-17) 這樣，利用(2-16)式就可以從變形速度分量中直接求出應力場。因此，當從運動許可速度場出發以節點速度為未知量求解塑性變形過程時，能夠簡便地求出變形和力能參數。2222222222222222222124()()()()3932124 ()()3932 3()6(9xyzvxyzvxyzxyyzzxxyzVVxyyzzxxyzVxygg 222222222221)21 ()()()6

26、()9yzzxVxyyzzxxyyzzxVgg(5) (5) 剛塑性可壓縮材料的流動法則剛塑性可壓縮材料的流動法則 根據(2-15)式可得： 將上式代入(2-11)式，可得：v239mmgd(2-18)219221993 2xmmxmxmxmxmmxmgdgdgdd(5) (5) 剛塑性可壓縮材料的流動法則剛塑性可壓縮材料的流動法則 上式表明，剛塑性可壓縮材料的流動法則是偏差應變速度分量與偏差應力分量成正比。與Levy-Mises流動法則一致，只不過是體積變形速率不等于零。32ymxyyzxmzmzxxmymzmxyyzzxd(2-19)(6) (6) 剛塑性可壓縮材料的特點剛塑性可壓縮材料的

27、特點 剛塑性可壓縮材料與剛塑性硬化材料的主要差別是，放松了體積不變條件的約束，即假設屈服與靜水壓力有關。體積變化率取決于靜水壓力，當求出材料的屈服應力、等效應變速率和給定材料的體積可壓縮參數后，可以直接從速度場求得應力場，所得結果與體積不可壓縮條件的解十分接近，而且計算過程得到簡化。(7) (7) 可壓縮參數可壓縮參數g g對體積變形的影響對體積變形的影響 將式(2-15)兩端乘以時間增量dt，可得：vmg(2-20)以單向壓縮應力狀態為例，其應力分量及等效應力為：0 0 0 0 03xyxyyzzxzm (2-21)9/39zmgg所以體積變形為：(2-23)(2-22)9vgg(7)(7)

28、 可壓縮參數可壓縮參數g對體積變形程度的影響對體積變形程度的影響 利用上式計算不同等效應變下，可壓縮參數g對體積變形程度的影響 表2-1 可壓縮參數g對體積變形程度的影響 g10.050.010.0010.00010.13.160.1660.03330.003330.000330.26.320.3320.06660.006670.000670.39.490.4990.09990.010000.00100 從上表可見，在1030%的等效應變條件下，可壓縮參數g取0.010.0001之間時，材料的體積變化不超過0.1%，可以近似滿足體積不可壓縮條件。因此，求解金屬材料的軋制過程，可壓縮參數g可取0

29、.010.0001。2.4 2.4 剛塑性材料的基本方程剛塑性材料的基本方程 剛塑性材料發生塑性變形時，由表面積 S 所圍成的體積 V 中，應力 、速度 和應變速率 應滿足下列基本方程: 1.力平衡微分方程(2-24)ijivij,0000 xyxxzyxyyzij jzyzxzxyzxyzxyz2.幾何方程(2-25),1211()2212yxxxxyyyziji jj iyyzxzzzzxvvvxyxvvvvvxzyvvvxxz2.4 2.4 剛塑性材料的基本方程剛塑性材料的基本方程2.4 2.4 剛塑性材料的基本方程剛塑性材料的基本方程3.本構關系(2-26)33223332223322

30、xxmxyxyijijyymyzyzzzmzxzxs2.4 2.4 剛塑性材料的基本方程剛塑性材料的基本方程 4.屈服準則 5.體積不可壓縮條件 6.邊界條件上在給定速度表面viiSvv ijjipnpS在給定外力表面上(2-27)(2-28)(2-29)(2-30)223ijijs skk0vijijxyz 2.4 2.4 剛塑性材料的基本方程剛塑性材料的基本方程式中： 、 等效應變速率和等效應力； 偏差應力； 屈服剪切應力； 外力表面單位法線矢量的方向余弦。 對于剛塑性材料發生塑性變形時，需要對上述基本方程進行聯立求解來確定變形區的速度、應變速率和應力。由于直接求解這些偏微分方程組非常困難

31、，因此，人們尋求其它求解途徑，即利用變分原理求解。 ijskjn2.5 剛塑性有限元的基本原理剛塑性有限元的基本原理l 剛塑性材料的變分原理是剛塑性有限元法的理論基礎，變分原理通過能量積分把偏微分方程組的求解問題變成了求泛函極值問題，從而為各種實際問題的求解提供了一種新方法。l 材料模型不同，變分原理的形式也不相同。根據附加條件的情況，變分原理可分為一般變分原理、不完全廣義變分原理和完全廣義變分原理。2.5.1 Mapkov原理原理 (第一變分原理第一變分原理) 設剛塑性變形體的體積為V，表面積為S0。在Sp上給定表面力 ，在Sv上給定速度 vi，在滿足變形的幾何條件、速度邊界條件和體積不變條

32、件的一切許可速度場中，真實的速度場使泛函 取極小值。ip23PsijijiiVSdVp v ds (2-31)2.5.2 剛塑性可壓縮材料的變分原理剛塑性可壓縮材料的變分原理 在滿足變形的幾何條件、速度邊界條件和體積微可壓縮條件的一切許可速度場 中，真實的速度場使泛函 取極小值。 ivpiiVSdVp v ds(2-32)2.5.3 剛塑性材料的廣義變分原理剛塑性材料的廣義變分原理 在第一變分原理中，運動許可速度場必須滿足幾何條件、速度邊界條件和體積不可壓縮條件。在處理實際問題時，有些條件比較容易滿足，而有的條件不易滿足，因此引用拉格朗日乘子把這些約束條件全部或部分引入總體能量泛函中，通過泛函

33、變分使這些約束條件得到滿足。引入約束條件后，變分原理的表達式要發生變化，統稱為廣義變分原理。引入部分約束條件的稱為不完全廣義變分原理，全部約束條件同時引入的稱為完全廣義變分原理。 (1) 不完全廣義變分原理不完全廣義變分原理 設定初始速度場時，幾何條件和速度邊界條件容易滿足，而體積不可壓縮條件則不易滿足，所以把體積不變條件通過拉格朗日乘子引入能量泛函中，得到新能量泛函： 在滿足幾何方程和速度邊界條件的一切速度場 中，真實的速度場使上述能量泛函取駐值。 ivpiiijijVSVdVp v dsdV (2-33)(2) 完全廣義變分原理完全廣義變分原理 在一切位移速度、應變速度和應力函數中，使能量

34、泛函 取駐值的vi、 和 必為剛塑性材料的正確解。ijij, ()1 ()2pviiijijiiiVSVSijiji jj iVdVp v dsdVvv dsvvdV(2-34)第一變分原理與廣義變分原理的差別第一變分原理與廣義變分原理的差別 第一變分原理要求能量泛函取最小值，而廣義變分原理僅要求能量泛函取駐值； 第一變分原理要求速度場滿足運動許可條件，靜力許可條件是通過變分過程近似滿足的，而廣義變分原理的速度場不受任何約束，所有方程均由變分過程近似滿足； 第一變分原理比廣義變分原理所得到結果精度高，但前者初始速度場的設定比后者難。 2.6 速率敏感材料的總體能量泛函速率敏感材料的總體能量泛函

35、 求解軋制問題時，總體能量泛函的表達式為： 式中： 塑性變形功率； 接觸表面摩擦功率； 外張力功率，前張力取負號，后張力取正號； 速度不連續面上的剪切功率。pfs(2-35)pfS2.6.1 內部塑性變形功率內部塑性變形功率 對于速率敏感性材料或在高溫下成型的金屬 ,單位體積的塑性變形功率： 速率敏感材料的內部塑性變形功率為： 001() 1nmijEdadm (2-37)1 1pVdVm (2-38)2.6.2 接觸表面上的摩擦功率接觸表面上的摩擦功率 軋制過程中，軋輥與軋件接觸表面上存在中性面，為了避免因摩擦力變向出現第二類奇異點而導致能量泛函不收斂的問題，可采用摩擦應力模型： 摩擦應力是

36、相對滑動速度vg的函數，因此摩擦功率為： 1223gsfgvmvk(2-38)2210()()3gsffggvssvmddsvkk ds (2-39)2.6.3 外張力功率外張力功率 式中： 作用有外張力的表面； 張應力，軋件受拉時取正值，受推時取負值； 相應表面處的位移速度。isv ds(2-40)siv2.6.4 速度不連續面上的剪切功率速度不連續面上的剪切功率 式中： 速度不連續面； 屈服剪應力； 速度不連續面上的速度不連續量。(2-41)ksssvkssssv ds速率敏感材料軋制過程總體能量泛函速率敏感材料軋制過程總體能量泛函 在軋制過程中，邊界上的外力為接觸表面上的摩擦應力、變形區

37、端部的張應力和速度不連續面上的剪切應力，相應的總體能量泛函為： 對于簡單軋制過程來說，總體能量泛函可簡化為： 221 1 ()13sgiksVssssmdVvkk dsvdsv dsm(2-42)221 1 ()13sgVsmdVvkk dsm(2-43)2.7 剛塑性有限元的求解途徑剛塑性有限元的求解途徑l設定運動許可速度場。l建立總體能量泛函，并把泛函表示成速度的函數。l利用數學上的極值理論求解泛函的駐值或最小值。l利用幾何方程和本構方程確定應變速度場和應力場。l通過接觸表面應力積分求出總軋制力和平均單位壓力，利用速度場求出軋件寬展、前滑和軋件側面形狀參數。2.7.1 軋制變形區的有限元離

38、散化軋制變形區的有限元離散化 (1) (1) 研究目標的選擇研究目標的選擇zox根據變形特征確定研究目標。利用對稱面上一個方向速度為零的特點，可使問題簡化。如果變形過程的幾何條件、物理條件在變形區內部及邊界上對稱，那么其真實解也必然對稱。圖2-2 板帶軋制過程的對稱面 (2) 單元類型的選取單元類型的選取 研究目標確定之后，便可著手在所選區域內設置節點、劃分單元。除了速度已知邊界之外，每個節點都有一組速度未知量，對于2維問題節點速度為 ,3維問題為 ，所以設置的節點越多，求解時未知數越多。單元類型的選擇，可根據所處理問題的特點來確定，求解平面變形問題，常用四邊形單元，求解3維變形問題常用立方體

39、單元。 yxvv 、zyxvvv、2.7.2 總能量泛函的離散化總能量泛函的離散化 把軋制變形區劃分為有限個單元之后，總能量泛函便可由單元能量泛函迭加求和來得到。下面以簡單軋制過程為例進行說明。(1) 單元能量泛函單元能量泛函 對于平面變形條件下的簡單軋制過程，設單元的面積為Se，外力已知邊界為Le，則單元能量泛函為： 當單元的速度插值函數設定之后，單元內的 便可確定。在給定材料屈服應力和摩擦邊界條件下，單元能量泛函實際上是單元節點速度的函數： ee221 1 ()13esgSLmdSvkk dLm(2-44)ij 、e112244(,)xyxyxyf vvvvvv(2-45)(2) 總能量泛

40、函總能量泛函 根據(2-44)式對所研究區域全部單元的能量泛函求和，便可得到總能量泛函：式中：nm單元總數； nl 接觸表面單元數目。ee221 1111 ()13enmnmnlsgSLeeemdSvkk dLm(2-46)(3) 接觸表面的速度邊界條件接觸表面的速度邊界條件 軋制過程接觸表面上的節點必須滿足速度邊界條件： tgyxvv (2-47)圖2-3 軋制特征角(4) 總體能量泛函總體能量泛函 由于接觸表面節點的兩個速度分量只有一個是未知量 ，同時對于 x 軸的節點垂直方向的速度已知(為零)，除去節點已知速度分量之后，把節點未知速度分量統一用 表示，總體能量泛函表示成未知速度分量的函數

41、：式中：nx未知量總數。xviv12( ,)inxf v vvv(2-48)2.7.3 總能量泛函的最小化總能量泛函的最小化 用剛塑性有限元法求解時，根據變分原理應求總體能量泛函的極小值。從式(2-48)可知，總體能量泛函是節點速度矢量的多元函數。根據多元函數求極值的條件，則有： 由于變分是任意的，必須有下列方程組成立：1 0 enxeeevvvv(2-49)10 enxeevv(2-50)(1) (1) 線性化處理線性化處理 采用Newton-Raphson方法線性化求解時，假定泛函連續并存在各階導數，則借助于泰勒級數在任一點展開。在第次迭代時，泛函可在第次迭代求得的值上展開，略去二階以上的

42、高階微量，從而得到以為未知量的線性方程組：211111 0 eenmnmnmeeeeTeeekkkkvvvvv211111 eenmnmnmeeeTeeeekkkvvvv (2-51)(2-52)(2) (2) 阻尼因子阻尼因子 采用(2-52)式求解時，先給定一個初始速度場 ，然后用迭代法按上式求出 ，直到相鄰兩次迭代速度近似解的偏差充分小時為止。迭代時，為了防止發散一般采用公式： vv1 kkvvv(2-53)01式中： 阻尼因子， 。2.8 剛塑性有限元的基本公式剛塑性有限元的基本公式 從上節可以看出，剛塑性有限元求解過程的一個基本步驟是計算單元的能量泛函計算單元的能量泛函。為此，需要解

43、決以下幾個方面問題：l設定初始速度場；l設定單元的速度分布函數插值函數；l建立單元變形速度與節點速度的關系；l計算單元能量泛函；l求單元能量泛函的一階偏導數即梯度；l求單元能量泛函的二階偏導數即Hessian矩陣。 2.8.1 單元坐標及速度插值、形狀函數單元坐標及速度插值、形狀函數 單元內任意一點的速度及坐標可以用單元節點的速度及坐標進行插值計算。對速度及坐標的插值使用相同的線性插值函數，該函數只與單元的形狀和節點的配置有關，因此，稱之為形狀函數或形函數，與此相應的單元稱為線性等參單元。(1)(1) 四邊形線形等參單元四邊形線形等參單元 為了便于插值與積分，單元內部采用局部坐標系 ，整體坐標

44、系的任意四邊形單元映射到局部坐標系后就變為正方形單元，單元節點編號及坐標系的映射變換如圖2-4所示。 42y13xO3124O( , ) 圖2-4 單元節點編號及坐標系的映射變換(2)(2) 單元局部坐標與整體坐標的變換單元局部坐標與整體坐標的變換 整體坐標系中4個節點坐標為 ，局部節點坐標 如下： 單元局部坐標 與整體坐標 的變換關系：( , ) , 14kkxyk (, )kkint(1)/ 2)12341231-1-1 11( 1)-1 1-11( 1)kkkk (2-54),(yx( , ) (2-56)41 122334414112233441( , )( , )kkkkkkxN x

45、N xN xN xNxyN yN yN yN yNy (3)(3) 四邊形等參單元的形函數及其性質四邊形等參單元的形函數及其性質 式(2-56)中 稱為單元的形函數，具有以下性質：l形函數 是與坐標插值函數相同形式的雙線性插值函數；l形函數 在節點上， ；l形函數 在節點上， 。 ( , )kN ( , )kN ( , )kN ( , )kN (,)1kkkN (,)0 , kllNlk (4)(4)形形函數表達式的確定函數表達式的確定 以 為例說明其確定過程： 1( , )N 42y13xO3124O上圖局部坐標系中，單元24、34兩條邊方程可表示為： 1-01-0(2-57)(4)(4)形

46、形函數表達式的確定函數表達式的確定 根據形函數的性質(3)， 在節點2、3、4處的值為零，而直線24、34經過這些節點，所以可取： 根據形函數的性質(2)， 在節點1處其值為1，把1點的局部坐標值代入上式得：1( , )N 1( , )(1)(1)NC (2-58)1( , )N 11(1 1)(1 1)4CC(2-59)11111( , )(1)(1)(1)(1)44N (2-60)(5)(5)四邊形等參單元的形函數四邊形等參單元的形函數 22233344411( , )(1)(1)(1)(1)4411( , )(1)(1)(1)(1)4411( , )(1)(1)(1)(1)44NNN 1

47、( , )(1)(1) (1 4)4kkkNk (2-61)(2-62)(6)(6) 單元內任一點的速度與節點速度的關系單元內任一點的速度與節點速度的關系 4141( , )( , )xkxkkykykkvNvvNv (2-63) evNv(2-64) xyvvv 4321432100000000NNNNNNNNN Tyxyxyxyxevvvvvvvvv443322112.8.2 單元變形速度與節點速度關系、單元變形速度與節點速度關系、B矩陣矩陣 (1) 平面變形的幾何方程平面變形的幾何方程xxyyyxxyvxvyvvyx(2-68)(2) 以節點速度表示的幾何方程 4414141411( ,

48、 )( , )()ixxxixixkxkykiyyyiiykykkyiixxyyixixyiNvvxxvNvvNvyyvNvvNNvvvxyyx (2-69)(3) 幾何方程的矩陣形式 11212342123431122334434400000000 xyxxyyxxyyxyvvvbbbbvccccvcbcbcbcbvvv (2-70)4) , 3 , 2 , 1( , , iyNcxNbiiii式中：(3) 幾何方程的矩陣形式 Txyxy4321BBBBB iiiiibccbB00 eBv令 (2-74)(2-71)(2-73)(2-75)BiB式中： 單元應變矩陣或B矩陣； 典型子矩陣。2

49、.8.3 Jacobi矩陣及其逆矩陣和行列式矩陣及其逆矩陣和行列式 單元應變B矩陣中的元素bi及ci是形函數 對x及y的偏導數，為了便于積分計算，形函數 對x及y的偏導數應該用 對 及 的偏導數來表示。根據復合函數的求導法則： ( , )iN ( , )iN ( , )iN TiiiiiTiiiiiNNNNNxyxyxyxyNNNNNxyxyxyxy(2-76)(1) 2維問題的維問題的Jacobi矩陣矩陣 將上式寫成矩陣形式：iiiiNxyNxNNxyy(2-77) xyJxy(2-78)(1) 2維問題的維問題的Jacobi矩陣矩陣 1( , )(1)(1)4( , )1(1) 4 (1

50、4)( , )1(1)4kkkiiiiiiNNiN (2-79)444111444111( , ) ( , )iiiiiiiiiiiiiiiiiiNNxyxyxNxJxyNNyNyxy (2-80)(2) Jacobi矩陣的逆矩陣及行列式矩陣的逆矩陣及行列式 1iiiiiiiiNNxyNNxxJNNNxyNyy(2-81)1 1 yxJyxJ(2-82)xyxyJ(2-83)2.8.4 高斯積分高斯積分 對于2維變形問題，建立單元內任意一點變形速度與節點速度關系之后，在給定初始速度場條件下，即可對單元能量泛函進行積分計算。在單元局部坐標系中，單元能量泛函式可表示成： 剛塑性有限元中的數值積分常

51、用高斯積高斯積分，即按照數學上的規則在單元內選取若干個積分點，用積分點處的函在單元內選取若干個積分點，用積分點處的函數值與求積系數之積的累加結果近似代替原積分數值與求積系數之積的累加結果近似代替原積分。 1112211111 ()13esgfmJ d dvkk Jdm (2-84)(1)單元塑性變形功率單元塑性變形功率高斯積分表示高斯積分表示 式中：n單元積分點個數，Hi、Hj求積系數。 41111111111 111nnepijijiJ d dH HJJmmm (2-85)G13124OG3 G2G4 對于2維問題采用線性單元時，n=2, Hi=Hj =1，此時單元內有4個高斯積分點，積分點

52、坐標為0.57735027，如圖2-5所示。(2)單元摩擦功率的積分單元摩擦功率的積分 單元摩擦功率的積分，可采用1維高斯積分近似計算：12211221122211()3 ()3 =()3esfgfnsigfisgfimvkk JdmHvkk Jmvkk J(2-86) 對于1維問題采用線性單元時，n=2, Hi=1，此時單元內有2個高斯積分點，積分點坐標為0.57735027。12OG1G2(3) 1維問題的維問題的Jacobi矩陣行列式矩陣行列式 式中：x1、x2為摩擦表面單元兩個節點的整體坐標值； N1、N2為1維單元的兩個形函數。212121ifiiNNNxJxxx(2-87)121(

53、1)21(1)2NN(2-88)(4) 單元的能量泛函單元的能量泛函 采用高斯積分后，2維變形簡單軋制過程的單元能量泛函為： 式中： 分別為單元高斯積分點處的等效應力、等效應變及Jacobi 矩陣行列式的值。42221111()13esgfiimJvkk JmfJJ、 、 和(2-90)2.8.5 矩陣分析中的公式矩陣分析中的公式(1) 數值函數對向量變量的偏導數數值函數對向量變量的偏導數TTTTA XAXX AAX 12TnAaaa12TxXxxx()TTf XA XX A設 是給定的向量， 是向量變量，且 ，則有 (2-a)(1) 數值函數對向量變量的偏導數數值函數對向量變量的偏導數 ()

54、Tf XA XAXXijn nAaT12xXxxx()Tf XX AX設 是給定的矩陣， 是向量變量，且 ，則有() 2Tf XAAAXX如果 ，則 (2-b)(2-c)(2) 向量函數對向量變量的偏導數向量函數對向量變量的偏導數 ijn nAaT12nXxxx()f XAX設 是給定的矩陣， 是向量變量，且 ，則有()()TTTf XAXf XAX(2-d)2.8.6 能量泛函的一階偏導數能量泛函的一階偏導數(梯度梯度)(1) 平面變形問題單元的等效應變速度列陣平面變形問題單元的等效應變速度列陣 平面變形問題單元的變形速度列陣 TexyxyB v222411931 xxyyxyvTTTgZC

55、Cg (2-91)(2-92)式中的常數矩陣和向量式中的常數矩陣和向量 4209924 0991003Z 1 1 0C (2-94)(2-93) TZ 222224209924 0991003422419999342241419999393xTxyxyyxyxxyxyxyyxyxxyxyyxyxxyyZ 2xy 1 TTCCg 21 1111 10011TTxxyxyyxyxyxyVCCgggg (1) 2維平面變形條件下的等效應變維平面變形條件下的等效應變 222 4111() 93TexyxyTTxxyyxyvB vZgg 1 TeTTeeTTTeCCvBZ B vvBCC B vg(2-

56、95)(1) 2維平面變形條件下的等效應變維平面變形條件下的等效應變 上式可以進一步寫成更為簡潔的形式： TeTTeZvBZB v412109912141 0991003TggZZCCggg (2-96)(2-97)等效應變不用等效應變不用(2-96)式簡潔形式？式簡潔形式？ 因為采用體積可壓縮方法進行求解時，如果等效應變采用(2-96)式簡潔形式，則要求單元內每一個高斯積分點處的體積變形速度都很小，這樣的約束條件過于苛刻，經常出現迭代計算發散現象，導致有限元數值求解過程無法進行。(1) 2維平面變形條件下的等效應變維平面變形條件下的等效應變 為了排除每一個高斯積分點處體積變形速度都很小這一嚴

57、格約束，單元的體積變形速度取各高斯積分點處體積變形速度的平均值，即：444111111 444eeeVkkkkkkCC BvCBvC B v1 eTTeeTTTevBZ B vvBCC B vg(2-98)(2-100)411 4kkBB(2-99)(2) 等效應變速度對速度向量的一階偏導數等效應變速度對速度向量的一階偏導數 等效應變速度(2-100)式對速度向量的一階偏導數：1 TeTTeePeBZBvBCCBvKvgv1 TTTPKBZBBCC Bg(2-101)(2-102) PK式中：8 8的矩陣。(3) 單元梯度單元梯度 把變形抗力模型代入單元塑性變形功率表達式，并對單元節點速度向量

58、求一階偏導數：44111144114411111111 nmenmPiiemnnmPeeeiiePeiiaJaJmmaJaJmvvvJJKvv(2-103)(4) 單元摩擦功率對速度向量的一階偏導數單元摩擦功率對速度向量的一階偏導數 接觸表面任意一點的速度邊界條件： vxvytgyxvv 接觸表面任意一點的相對滑動速度：secgRxvvv(2-104)單元內任意一點的速度與節點速度的關系單元內任意一點的速度與節點速度的關系 對于2維軋制過程的接觸表面單元而言，單元節點上只有一個未知量，因此，單元內任意一點的速度與節點速度的關系如下： efexffvvNv120.5(1- ) 0.5(1) fT

59、efxxNvvv(2-105)(2-106) 對節點速度向量的一階偏導數對節點速度向量的一階偏導數gv 把(2-105)式代入(2-104)式可得相對滑動速度： secegRffvvNv(2-107)sec gTfefvNv 2222()22 secgggeeffggTfefvkkvvvvvkvNv (2-108)(2-109)單元能量泛函的一階偏導數即梯度單元能量泛函的一階偏導數即梯度 接觸表面單元摩擦功率的一階偏導數： 單元能量泛函的一階偏導數即梯度：22221111() sec 33egfTssfffeeiiffvkkmmJJNvv 42111 sec3eefePefeeTsPffiiv

60、vmJKvJN(2-110)(2-111)2.8.7 能量泛函的二階偏導數能量泛函的二階偏導數(Hessian矩陣矩陣) 1223 (1) (1) (1) nmeTeTnmeTePeTeTePavvmavKvmvvmKv 2.8.7 能量泛函的二階偏導數能量泛函的二階偏導數(Hessian矩陣矩陣)412421324221 (1) (1) eePPeieePPPeeTeTiTePeTTeeePPPPeeTiJKvvJKJKvvvvmKvvKvKvJKmvv 由于 (2-112)2.8.8 總體能量泛函能量泛函的梯度矩陣和總體能量泛函能量泛函的梯度矩陣和Hessian矩陣矩陣 根據迭加原理對所有