行政院國家科學委員會專題研究計畫 成果報告
智慧型金屬板材成形 CAE 軟體之開發與整合研究(3/3)
計畫類別： 整合型計畫
計畫編號： NSC91-2212-E-011-026-
執行期間： 91 年 08 月 01 日至 92 年 07 月 31 日 執行單位： 國立臺灣科技大學機械工程系
計畫主持人： 黃佑民 計畫參與人員： 陳聰嘉
報告類型： 完整報告
處理方式： 本計畫可公開查詢
中 華 民 國 92 年 10 月 23 日
行政院國家科學委員會補助專題研究計畫 成 果 報 告
□期中進度報告
智慧型金屬板材成形 CAE 軟體之開發與整合研究(3/3) ---總計畫
計畫類別：□ 個別型計畫 整合型計畫 計畫編號：NSC 89－2212－E－011－036－
NSC 90－2212－E－011－041－
NSC 91－2212－E－011－026－
執行期間： 89 年 08 月 01 日 至 92 年 07 月 31 日
計畫主持人：黃 佑 民 共同主持人：
專任助理： 陳靜芬 計畫參與人員：陳聰嘉
成果報告類型(依經費核定清單規定繳交)：□精簡報告 完整報告
本成果報告包括以下應繳交之附件：
□赴國外出差或研習心得報告一份
□赴大陸地區出差或研習心得報告一份
□出席國際學術會議心得報告及發表之論文各一份
□國際合作研究計畫國外研究報告書一份
處理方式：除產學合作研究計畫、提升產業技術及人才培育研究計畫、
列管計畫及下列情形者外，得立即公開查詢
□涉及專利或其他智慧財產權，□一年□二年後可公開查詢
執行單位：國立台灣科技大學 機械工程學系
中 華 民 國 九十二 年 十 月 三十一 日
摘 要
本計畫為整合型計畫，分成四項子計畫，以開發金屬板材成形三 維分析 CAE 軟體有關工件變形及應力與應變分佈之研究，提供於塑 性加工時做為預估成形性之依據。在應力與應變之關係式方面擬以 updated Lagrangian formulation (ULF)理論為基礎並結合有限變形的 觀念與rmin技巧來撰寫三維增量型彈塑性大變形有限元素計算程式。
本研究計畫所建立之分析模式，能完整的描述各項變形履歷，並發展 為泛用型之金屬板材成形三維分析CAE 軟體。
關鍵詞：金屬成形，CAE 軟體，彈塑性分析，突緣引伸
Abstract
The main purpose of this proposal is to developing the deformation as well as stress and strain distributions of CAE software of metal forming processes. A methodology for formulating an elasto -plastic finite element CAE software, which is based on an updated Lagrangian formulation (ULF), will be developed to simulate the 3D metal forming process. The CAE computer code will be established to analyze the 3D manufacturing process of sheet metal.
Keyword: metal forming, CAE software, elasto-plastic analysis, flanging drawing.
目 錄
摘 要... I ABSTRACT... II 目 錄...III 圖表索引...VII
一、前 言...1
二、基本理論...7
2.1 基本假設 ...7
2.2 有限變形之應變與應變率 ...7
2.3 有限變形之應力與應力率 ...10
2.4 有限變形之 TOTAL LAGRANGIAN FORMULATION...13
2.5 有限變形之 UPDATE LAGRANGIAN FORMULATION...16
2.6 材料之彈塑性構成關係式 ...17
三、數值分析...23
3.1 簡介...23
3.2 虛速度原理的離散化 ...25
3.3 退化殼元素（DEGENERATED SHELL ELEMENT）之分析...26
3.4 不同積分法則推導退化殼元素之剛性矩陣 ...28
3.5 摩擦處理 ...30
3.6 廣義RMIN法之增量步驟的計算...33
四、金屬板材突緣引伸成形製程之分析...36
4.1 邊界條件 ...36
4.2 材料參數 ...36
4.3 數值模擬與分析 ...37
五、金屬板材圓杯引伸成形製程之分析...47
5.1 邊界條件 ...47
5.2 材料參數與實驗數據 ...48
5.3 數值模擬分析 ...48
5.4 數值分析與實驗結果之比較 ...50
5.4.1 沖頭負荷之比較 ...51
5.4.2 CPU 運算時間之比較 ...52
5.4.3 異向性材料之模擬分析 ...53
5.5 圓杯引伸成形歷程 ...53
5.5.1 圓杯引伸成形歷程之比較 ...54
5.5.2 耳緣成形之分析 ...55
5.5.3 厚度分佈之比較 ...55
5.5.4 應力分佈之比較 ...57
5.5.5 應變分佈之比較 ...58
六、金屬板材深引伸耳緣成形製程之分析...94
6.1 實驗設備 ...94
6.2 實驗原理與步驟 ...94
6.3 邊界條件 ...96
6.4 材料參數 ...97
6.5 數值模擬分析 ...98
6.6 數值分析與實驗結果之比較 ...99
6.6.1 工件之應力分佈 ...99
6.6.2 工件之應變分佈 ...100
6.6.3 工件之厚度分佈 ...100
6.6.4 沖頭負荷與衝程關係之比較 ...101
6.6.5 深引伸耳緣成形製程之成形歷程 ...101
七、金屬板材彎曲成形製程之分析...108
7.1 實驗與材料參數 ...108
7.2 數值分析 ...109
7.3 除荷之設定 ...110
7.4 有限元素模式之實驗驗證 ...111
7.5 數值模擬結果討論 ...111
八、金屬板材方杯引伸成形製程之分析...124
九、平行處理之分析...136
9.1 數值分析 ...139
9.2 結果與討論 ...140
十、參考文獻...147
十一、誌 謝 ...150
十二、計劃成果自評...151
圖表索引
圖1-1 智慧型 CAE 軟體開發流程圖 ...6
圖2-1 連續體在t =t0與t= 時之變化狀況 ...22t 圖2-2 定義應力之連體座標 ...22
圖2-3 平衡時之連體體積與曲面定義 ...22
圖3-1 退化殼元素之座標 ...35
圖3-2 殼元素之自由度 ...35
圖3-3 接觸彈簧元素 ...35
圖4-1 帽型引伸加工之工件與邊界條件...40
圖4-2 帽型引伸加工之沖頭 ...40
圖4-3 帽型引伸加工之模具 ...41
圖4-4 帽型引伸加工之壓料板 ...41
表4-1 帽型引伸成形分析，模具與料片之網格分割數據...42
圖4-5(A) 帽型引伸衝程 20MM厚度圖...42
圖4-5(B) 帽型引伸衝程 40MM厚度圖...43
圖4-6(A) 帽型引伸衝程 20MM應力圖...43
圖4-6(B) 帽型引伸衝程 40MM應力圖...44
圖4-7(A) 帽型引伸衝程 20MM應變圖...44
圖4-7(B) 帽型引伸衝程 40MM應變圖...45
圖4-8(A) 帽型引伸衝程 20MM反力圖...45
圖4-8(B) 帽型引伸衝程 40MM反力圖...46
圖5-1 圓杯引伸成形之胚料與模具系統邊界條件...60
圖5-2 圓杯引伸成形之模具與胚料幾何尺寸圖...61
圖5-3 圓杯引伸成形之胚料與模具系統立體圖...61
圖5-4 圓杯引伸成形實驗的沖頭尺寸 ...62
圖5-5 圓杯引伸成形實驗的下模尺寸 ...62
圖5-6 圓杯引伸成形實驗的壓料板尺寸...63
圖5-7 第一類型圓杯引伸成形模擬之沖頭網格分割...63
圖5-8 第一類型圓杯引伸成形模擬之下模網格分割...64
圖5-9 第一類型圓杯引伸成形模擬之壓料板網格分割...64
圖5-10 第一類型模擬之整體圓形胚料區隔成十二個區域...65
圖5-11 第一類型圓杯引伸成形模擬之整體圓形胚料網格分割...65
圖5-12 第二類型模擬之整體圓形胚料區隔成五個區域...66
圖5-13 第二類型圓杯引伸成形模擬之整體圓形胚料網格分割...66
圖5-14 第三類型圓杯引伸成形模擬之沖頭網格分割...67
圖5-15 第三類型圓杯引伸成形模擬之下模網格分割...67
圖5-16 第三類型圓杯引伸成形模擬之壓料板網格分割...68
圖5-17 第三類型模擬之四分之一圓形胚料...68
圖5-18 第三類型圓杯引伸成形模擬時四分之一圓形胚料之網格分割 ...69 圖5-19 第一類型圓杯引伸成形模擬時二分之一圓形胚料之網格分割 ...69 圖5-20 第一類型圓杯引伸成形模擬時四分之一圓形胚料之網格分割 ...70 圖5-21 第一類型圓杯引伸成形模擬時二分之一圓形胚料之邊界條件 ...70 圖5-22 第一類型圓杯引伸成形模擬時四分之一圓形胚料之邊界條件 ...71 表5-1 第一、二類型模具與胚料有限元素網格分割之相關數據...71 圖5-23 第一類型圓杯引伸整體圓形胚料數值模擬與實驗之沖頭負荷 與衝程關係之比較...72 圖5-24 第一類型圓杯引伸不同網格分割時數值模擬之沖頭負荷與衝 程關係之比較...72 圖5-25 第二類型圓杯引伸整體圓形胚料網格分割時數值模擬與實驗 之沖頭負荷與衝程關係之比較...73 表5-2 第三類型模具與胚料有限元素網格分割之相關數據...73 圖5-26 第三類型圓杯引伸成形模擬時四分之一圓形胚料網格分割之 邊界條件...74 圖5-27 第三類型圓杯引伸四分之一圓形胚料數值模擬與實驗之沖頭 負荷與衝程關係之比較...74
圖5-28 圓杯引伸成形模擬 CPU 運算時間之比較 ...75
圖5-29 圓杯引伸等向性及異向性四分之一圓形胚料數值模擬與實驗 之沖頭負荷與衝程關係之比較...75
圖5-30 圓杯引伸成形模擬時工件之變形歷程...76
圖5-30 圓杯引伸成形模擬時工件之變形歷程（續）...77
圖5-30 圓杯引伸成形模擬時工件之變形歷程（續）...78
圖5-30 圓杯引伸成形模擬時工件之變形歷程（續）...79
圖5-30 圓杯引伸成形模擬時工件之變形歷程（續）...80
圖5-31 等向性與異向性四分之一圓形胚料凸緣輪廓之比較...81
圖5-32 等向性與異向性四分之一圓形胚料耳緣輪廓之比較...81
圖5-33 異向性四分之一圓形胚料數值模擬與實驗厚度分佈之比較82 圖5-34 等向性與異向性四分之一圓形胚料數值模擬厚度分佈之比較 ...82
圖5-35 等向性二分之一與四分之一圓形胚料數值模擬厚度分佈之比 較...83
圖5-36 等向性二分之一圓形胚料最終變形圖...83
圖5-37 圓杯引伸成形模擬時工件之應力分佈...84
圖5-37 圓杯引伸成形模擬時工件之應力分佈（續）...85
圖5-37 圓杯引伸成形模擬時工件之應力分佈（續）...86
圖5-37 圓杯引伸成形模擬時工件之應力分佈（續）...87
圖5-37 圓杯引伸成形模擬時工件之應力分佈（續）...88
圖5-38 圓杯引伸成形模擬時工件之應變分佈...89
圖5-38 圓杯引伸成形模擬時工件之應變分佈（續）...90
圖5-38 圓杯引伸成形模擬時工件之應變分佈（續）...91
圖5-38 圓杯引伸成形模擬時工件之應變分佈（續）...92
圖5-38 圓杯引伸成形模擬時工件之應變分佈（續）...93
圖6-1 實驗設備整體系統配置圖 ...102
照片6-1 實驗模具圖 ...102
圖6-2 深引伸耳緣成形之模具幾何尺寸...103
圖6-3 深引伸耳緣成形製程之料片邊界條件...103
圖6-4 深引伸耳緣成形數值分析之最大主應力分佈圖...104
圖6-5 深引伸耳緣成形數值分析之最大主應變分佈圖...104
圖6-6 衝程達 60.0MM時之工件厚度分佈圖...105
照片6-2 深引伸耳緣成形之工件 ...105
圖6-7 深引伸耳緣成形數值分析與實驗之沖頭負荷與位移關係之比 較...106
圖6-8 深引伸成形模擬時工件之成形歷程圖...107
圖7-1 馬鞍成形製程之模具幾何尺寸 ...114
照片7-1 馬鞍成形製程在不同沖程之成品(BLANK 3) ...114
圖7-2 馬鞍成形製程模擬時沖頭之網格分割...115
圖7-3 馬鞍成形製程模擬時下模之網格分割...115
圖7-4 馬鞍成形製程模擬時板材之網格分割...115
表7-1 模具與板材有限元素網格分割之相關數據...116
圖7-5 馬鞍成形製程數值模擬與實驗之沖頭負荷與沖程關係之比較 ...116
圖7-6 馬鞍成形製程除荷後彎曲角度之計算值與實驗值...117
圖7-7 馬鞍成形製程模擬時板材之變形履歷...117
圖7-8 馬鞍成形製程模擬時不同沖程之應力分佈圖...119
圖7-9 馬鞍成形製程彎曲軸兩端之撓曲情形...119
圖7-10 馬鞍成形製程模擬時之不同沖程主應變分佈圖...121
圖7-11 馬鞍成形製程模擬時不同沖程之反作用力分佈圖...122
圖7-12 馬鞍成形製程模擬時除荷後之厚度分佈圖...123
圖7-13 馬鞍成形製程不同沖程及不同板材寬度時的馬鞍外型...123
表8-1 方杯成形製程之分析數據 ...129
圖8-1 負荷與衝程之關係圖 ...130
表8-2 負荷與衝程之分析數據 ...131
圖8-2 VON MISES 應力分佈圖...132
圖8-3 應變能分佈圖 ...132
圖8-4 變形履歷與等價應力分佈圖（PIII-1G） ...133
圖8-5 變形履歷與應變能分佈圖（PIII-1G） ...134 圖8-6 應變能分佈圖（PIII-1G）...135 表9-1 ALPHASERVER DS20E 6/667 UNIX PACKAGED SYSTEM配備說明
...142 圖9-1 程式未使用平行處理與使用平行處理模式分析深引伸耳緣成
形負荷之比較...143 圖9-2 程式未使用平行處理與使用平行處理模擬深引伸耳緣成形沿
輥軋方向之厚度比較圖...143 圖9-3 程式未使用平行處理與使用平行處理模擬深引伸耳緣成形與
輥軋方向垂直之厚度比較圖...144 圖9-4 深引伸耳緣成形最大主應力分佈圖...144 圖9-5 深引伸耳緣成形最大主應變分佈圖...145 圖9-6 成形工件之數位影像擷取並與數值分析結果套疊之比較....145 圖9-7 程式未使用平行處理與使用平行處理模擬深引伸耳緣成形所
需CPU 時間之比較...146 圖9-8 程式未使用平行處理與使用平行處理模擬深引伸耳緣成形在
不同元素數量下所需CPU 時間之比較 ...146
一、前 言
國內汽車工業、家電工業等產業，其產品製造的過程中，金屬板 材成形的加工製程所佔比率相當高，且隨著消費者導向的時代來臨，
新規範的產品開發、設計越來越多，所開發的產品也朝著少量而多樣 化的方向發展，每一新產品之生命週期越來越短。以汽車為例，每隔 幾年各汽車製造廠為了滿足顧客的需求，必須要將汽車改款，使同一 款式的產品週期越來越短，因此，以往靠經驗與感覺來設計→試作→
試驗→修正設計→再試作→再試驗等步驟之傳統設計方法浪費工 時，降低產品的競爭力。藉由CAE 軟體可將設計→試作→試驗→修 正設計→再試作→再試驗等反覆的步驟，利用電腦支援而縮短開發及 設計時程因而降低成本，亦可提升新產品的品質與設計最佳化。
金屬板金成形之數值模擬為一複雜且困難的工作，因成形之過程 為非線性問題，而造成此非線性的因素大致可分為三個部分：第一為 工件變形的大位移、大旋轉與大應變所造成幾何上的非線性行為；第 二為因變形而產生的材料非線性構成關係；第三為成形製程中，料片 與模具接觸界面之幾何外形的函數之非線性關係。上述的各種非線性 問題，皆增加了金屬板金成形製程模擬的複雜性。本研究之方法乃結 合有限變形理論ULF(updated Lagrangian formulation) 的觀念，採用 Cauchy 應力的 Jaumann rate 關係式，並忽略元素的密度或體積之變化
下，建立一個 incremental updated Lagrangian formulation 之彈塑性大 變形、大應變有限元素解析程式。本程式係以沖頭之增量位移為每一 變形增量步驟之起始增量值，文中採用 Yamada 的 rmin方法來判斷材 料內任一元素的彈塑性狀態變換問題，並將此方法擴展至料片上節點 之接觸、分離的判斷，摩擦方向的轉換及有關最大應變增量、最大旋 轉增量的線性化處理，以決定板金成形時每一增量的位移、應變、應 力與除荷後料片之回彈現象與工件最後形狀。
彎曲成形、拉伸成形及深引伸成形為板金成形之基本加工方法，
而本計畫所探討之各種板金成形加工，其成形加工過程皆是將料片放 置於沖模上，藉由沖頭對沖模之相對位移產生成形，其操作程序簡 單，但在實際的生產過程中則包含了幾個重要的技術問題，如加工除 荷後回彈作用(即 spring-back 或 spring-forward)的模具設計問題、成形 製程形狀精確的控制要求、成形性的極限問題、表面拉伸的破壞之形 成及負荷的預估等。近年來，仍有許多研究學者將有限元素法應用於 分析板金成形問題上，且電腦的快速發展，使得其運算能力與速度大 幅提昇，在處理複雜且高階非線性問題時可縮短其計算的時間，而這 些研究成果可提供工業界作為板金成形實際加工之參考依據。而與本 研究相關文獻有Yamada、Yoshimura 與 Sakurai [1]依據 von Mises 降 伏準則與Prandtl-Reuss 塑流法則，推導出等向性材料之彈塑性應力—
應變關係矩陣，並以rmin法控制元素由彈性進入塑性狀態時，其應力 狀態能先達到降伏曲面上。Hibbitt、Marcal 與 Rice [2]提出參考未變 形前的最初始狀態，以Total Lagrangian Formulation (TLF)推導出適用 於大變形的彈塑性有限元素模式，並建立增量型式的分析方法。Wang [3]計算金屬圓形胚料受沖頭拉伸引伸之大塑性變形，並使用塑性增量 理論分析半球狀沖頭拉伸引伸圓形胚料時之應力與應變分佈情形。
1975 年 McMeeking 與 Rice[4]引用 Hill[5]的變分原理，藉由參考體瞬 時狀態(current state)的 ULF 建立了彈塑性大變形有限元素分析模 式，由於邊界上的處理較為簡單，故較適合分析大變形金屬成形問 題。Pillinger、Hartley、Sturgess 及 Rowe[6]運用線性化的增量方式來 處理近似非線性的大變形旋轉效應，以保持增量之精確性。Kim、Oh 與 Kobayashi [7]以薄膜理論之有限元素法探討半球形沖頭拉伸引伸 成形時的料片厚度變化與相對位置之關係。Liu 等[8]利用退化非線性 殼理論推導合應力退化殼元素，並以圓柱體受到軸向負荷產生挫屈為 例進行數值分析。Keum 等[9]提出可適用於任意剛體模具外形，且僅 需要描述位置資料之模具曲面的新方法，而此方法已由基於 B-曲線 (B-spline)及線性內插法所描述的一般曲面的剛－黏塑性有限元素分 析程式所完成，並驗證此方法較以往所用的三維有限元素分析法更有 效率。Makinouchi 與 Kawka [10]研究有限元素模擬模鍛(stamping)模
具之設計，藉由數值模擬來預測與估計板金成形過程中的缺陷以替代 傳統的試誤法，其理論背景與數值分析是由有限元素程式完成。
Kawka 與 Makinouchi [11]亦提出退化殼元素的構想在靜態顯函有限 元素程式中執行，主要探討各種不同的積分法則比較，分別為完全積 分(FI)、減化積分(RI)、選擇減化積分(SRI)、假設應變場(ASF)與穩定 矩陣(SM)法則等。Lin 與 Tseng[12]為了解決金屬成形在彈塑性區間的 接觸問題，而應用一表面彈簧元素並將其元素之對稱性剛性矩陣疊加 於整體剛性矩陣中，組成有限元素之分析方法，經由實驗結果與其數 值分析結果之驗證，證明其方法可準確地模擬實驗結果。Lange [13]
對板金彎曲系列加工作了廣泛且深入的介紹，Weinmann [14]對熱軋 高強度低合金板金作一系列實驗，探討最大 air bending 負荷與受 coining 作用下的彈性除荷過程受沖頭半徑、模具寬度與板金厚度的 影響與變化，Magnusson 與 Tan [15]以純彎矩彎曲的初等解理論分析 V 型彎曲加工，Huang、Lu 與 Makinouchi [16，17]以彈塑性二維有限 元素增量型式分析平面應變型式的無摩擦 V 型彎曲製程，同時揭露 彎曲除荷時spring-back 與 spring-forward 的回彈現象，並計算由起始 彎曲變形至最後除荷過程的全部彎曲成形形狀與應力、應變分佈，
Makinouchi [18]以 rmin方法處理無摩擦作用的 draw bending 製程的接 觸問題，成功的表達了板材與模具的接觸現象。
本計畫為整合型計畫，其總體目標在於開發一套分析三維金屬板 材成形製程分析之有限元素 CAE 軟體，並於完成後能推廣且應用於 工業界。其各子計畫間之整體相關架構及研究目標，如圖1-1 所示。
本研究之有限元素分析是採用四節點四邊形退化殼元素推導剛性矩 陣，並以商用套裝軟體進行前後處理，將所建立的模具與板材經過網 格分割後轉換成數據檔，輸入增量型彈塑性大變形三維有限元素分析 程式中進行數值解析，來探討各種板材之成形製程，將所模擬出之沖 頭負荷與沖程之結果與實驗結果做比較，並且探討板材厚度分佈，最 後將分析結果輸出至商用套裝軟體顯示變形圖以及應力與應變分佈 圖，以作為模具設計與加工製程時之參考。
智慧型金屬板材成形CAE軟體 之開發與整合研究
金屬板材突緣引伸 成形製程之分析
金屬板材深引伸耳 緣成形製程之分析
金屬板材彎曲成形 製程之分析
金屬板材方杯引伸 成形製程之分析
CAE軟體之後置處 理 CAE軟體之應力與
應變分析 CAE軟體之剛性矩
陣建構 CAE軟體之前置處
理
子計劃一 子計劃二 子計劃三 子計劃四
初步CAE軟體架構
子 計 劃 名 稱
第 一 年
CAE軟體之測試 第
二 金屬板材突緣引伸 年
分析與實驗
金屬板材深引伸耳 緣成形分析與實驗
金屬板材彎曲成形 分析與實驗
金屬板材方杯引伸 成形分析與實驗
智慧型CAE軟體
將程式系統化 試用平行計算改善
程式運算效率 將材料之破裂分析
整合到程式中 將程式簡化為PC版
第 三 年
圖1-1 智慧型 CAE 軟體開發流程圖
二、基本理論
2.1 基本假設
本研究所考慮之材料彈塑變形時，其基本假設條件如下：
（1） 材料為均質性（homogeneous）及平面異向性（planar anisotropic）。
（2） 材料於加工時不考慮溫度影響。
（3） 材料於彈性區域內滿足虎克定律。
（4） 材料遵循 Hill 降伏法則，且於塑性變形後，依循 Prandtl-Reuss 之流動法則。
（5） 材料之構成方程式包含等向性之應變硬化影響。
（6） 材料之應變率可分為彈性與塑性應變率兩部份。
（7） 沖頭、沖模與壓料板皆視為剛體。
2.2 有限變形之應變與應變率
在處理有限變形之彈塑性問題時，因其變形是屬於大變形問題，
須考慮物體之形狀變化與位置變化，故基於不同的座標系統（變形前 和變形後的座標）所描述之物理量將有所差異。在本節及下節中討論 各種不同的應變、應變率、應力與應力率之間的關係，並分別以張量 和矩陣的形式表示。
如圖 2-1 所示，物體某微小線段變形前為AoBo，變形後為AB。
其應變張量(strain tensor)可定義為
) 2(
1
ij j i
ij a
X x X g x
E −
∂
∂
∂
= αβ ∂ α β (2.1)
或
) 2(
1
j i
ij
ij x
X x a X
g
A ∂
∂
∂
− ∂
= αβ α β (2.2)
E 即為 Green 應變張量，係以變形前的座標系統ij X 為準而定義之應i
變。A 則為 Almansi 應變張量，係以變形後的座標系統ij x 為準所定義i 之應變。其中a 為座標系統ij Xi(i =1,2,3)之度量張量(metric tensor)，gij 為座標系統xi(i =1,2,3)之度量張量。座標系統若採用直角卡式座標 (rectangular Cartesian coordinates)或埋入座標(convected coordinates)， 則aij =gij =δij(δ 為ij Kronecker’s symbol)，則Green 應變以張量形式 表示為
) 2(
1 ij
j k i
ij k X
x X
E x ∂∂ −δ
∂∂
= （2.3）
或以矩陣形式表示為 ]) [ ] ][
2([
] 1
[E = J J T − I
其中[J 為 Jacobian 矩陣，] [J]T為[J 之轉置矩陣（transpose matrix），] ]
[I 則為單位矩陣。Green 應變率為 Green 應變對時間的導數，故（2.3）
式微分可得應變率
) 2(
1
j k i k j
k i ij k
X x X x X
x X E x
∂
∂
∂ + ∂
∂
∂
∂
= ∂
•
• •
（2.4）
或
) ] ][
[ ] ][
2([
] 1
[E• = J• J T + J J• T
其中[J 為 •] ] ][
[ ]
[J• = J L （2.5）
其中[L 是以變形後空間固定座標] x,y,z(x1,x2,x3)為基準而定義的物 理量(Lij =∂vi/∂ xj)，稱為速度梯度矩陣（velocity gradient matrix）。
又[J•]T =[L]T[J]T故[E 可改寫成 •]
T T
T
T J L J J J
J L J
E ([ ][ ][ ] [ ][ ] [ ] ) [ ][ ][ ] 2
] 1
[• = + = ε• （2.6）
其中
) ] [ ] 2([
] 1
[ε• = L + L T 或 ( ) 2
1
i j j
ij i
x v x
v
∂ +∂
∂
= ∂
ε• (2.7)
]
[E 稱為 Green 應變率矩陣，係與變形前的座標• X,Y,Z有關的量。[ε•] 是根據物體不斷變形而改變的固定座標x ,,y z為基準而定義的量，稱 為Cauchy 應變率矩陣，為 ULF（updated Lagrangian formulation）所 採用的應變率。
變形前之微量體積dV 與變形後之微量體積 dV 之關係為 0
J X J
x dVdV
j
i = =
∂∂
= det[ ]
0
（2.8）
其中 J 為 ][J 的行列式值。且由質量不變定律可得
ρ ρ0
=
J （2.9）
其中ρo與ρ分別為變形前與變形後之密度。
2.3 有限變形之應力與應力率
在處理有限變形之彈塑性問題時，應力則須給予明確的定義，一 般常使用之應力有 Lagrange 應力、Kirchhoff 應力和 Cauchy 應力三 種。本節對此三種應力之定義加以說明，並推導在 ULF 下應力率之 間的關係。如圖2-2 所示，左側為變形前之狀態，右側為變形後之狀 態。
（1）Lagrange 應力(sij)
} { ] [ } {
) ;
( T o
i o oj L ji
oi s n dS dF dF s dS
dF = = = （2.10）
其中sji為 Lagrange 應力張量，具有非對稱形式，以變形前之面 積dS 為基準，其中o n 為oj dS 之單位法線向量之分量。 o
（2）Kirchhoff 應力(t ij)
} { ] [ ] [ } {
) ;
( i T T o
o oj K ji
oi dF dF J t dS
x dS X
n t
dF = = ∂∂ α =
α （2.11）
其中t 為 Kirchhoff 應力張量。 ji
（3）Cauchy 應力(σ ij)
} { ] [ } {
; dF dS dS
n
dFi = σji j = σ T （2.12）
其中σ 為 Cauchy 應力張量，變形後之力dF，以變形後之面積dSji
為基準所得之應力，係一對稱之應力張量，即σij =σji，n 為 dSj 之單位法線向量之分量。
為了求得各種不同定義應力間之關係，由質量不變定律得知
o i oj
o j
i n dS
x dS X
n ∂
∂ ρρ
= （2.13）
由（2.10）式與（2.12）式可得 Lagrange 應力和 Cauchy 應力間之關 係為
] [ ] [ ]
[s = J G T σ （2.14）
同理，由（2.11）式與（2.12）式可得 Kirchhoff 應力和 Cauchy 應力 間之關係為
] ][
[ ] [ ]
[t = J G T σ G （2.15）
其中[G]= J[ ]−1，即[J][G]=[I]，整理上列關係式後，可獲得各種應 力間之關係為
] [ ] [ ] ][
[ ]
[s = t J = J G T σ （2.16）
（2.16）式對時間的導數，可得到各種應力的變化率 ]
][
[ ] ][
[ ]
[s• = t• J + t J• （2.17）
] [ ] [ ] [ ] [ ] [ ] [ ]
[s• = J• G T σ + J G• T σ + J G T σ• （2.18）
在 ULF 的場合裡，因是以每一時刻下物理性質和幾何形狀為基準，
故以下的關係式成立。
] [ ] [ ]
[J = G = I ] [ ] [J• = L
] [ ]
[G• = − L ] [ ] [E• = ε•
] [ ] [ ]
[s = σ = t
=1 J
ρ ρ•
• =−
J
將上列各式代入（2.17）式及（2.18）式中，得到 ]
][
[ ] [ ]
[s• = t• + σ L （2.19）
] [ ] [ ] [ ] [ ]
[ σ σ
ρ
σ ρ L T
s = − −
• •
• （2.20）
（2.19）式的 ][t 為 Kirchhoff 應力的埋入微分（convected rate），係以• 埋入座標X,Y,Z變形後所產生的新方向來定義，亦即觀察者隨物體一 起變形與轉動所看到的應力變化率。（2.20）式的[σ•]係根據空間固定 座標系x ,,y z所定義，變形後曲面上的作用力、變形後之面積及時間 之關係來定義，稱為Cauchy 應力的固有微分（intrinsic rate），是觀察 者在固定座標系上所看到的應力變化率。另外定義一種應力變化率，
稱為應力的Jaumann 微分，即觀察者隨著物體轉動，但不隨著變形而 觀察所得到的應力變化率。設時間t 時，空間固定座標系X,Y,Z 的 Cauchy 應力為 ][σ ，轉動座標系的 Cauchy 應力為[σω]，兩應力間如
使用一座標轉換矩陣[T ，則得下式成立 ] ]
][
[ ] [ ]
[σω = T T σ T （2.21） 微分（2.21）式，則得Cauchy 應力的Jaumann微分[σ o]
[ ]
[ ][ ]] ][
[ ] [ ] ][
[ ] [ ]
[σo = T• T σ T + T T σ• T + T T σ T• （2.22）
在ULF 時，[T]T =[T]=[I]，[σω]=[σ]其中[T•]之具體形式為
] [ ] [ ] [ 0
0 0
]
[T L
x y
x z
y z
−
=
=
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
−
−
−
= • •
•
•
•
•
•
•
• ω ε
ω ω
ω ω
ω ω
（2.23）
其中
) 2(
1
j i i
k j
x v x
v
∂
− ∂
∂
= ∂
ω• （2.24）
]
[ω• 稱為轉動矩陣。將（2.23）式代入（2.22）式，可得 ]
][
[ ] ][
[ ] [ ]
[σo = σ• + σ ω• − ω• σ （2.25）
結合（2.25）式與（2.20）式，可得
] [ ] ][
[ ] ][
[ ] [ ] [ ] [
] [ ] ][
[ ] [ ] [ ] [ ] [
ρ σ σ ρ
ε σ σ
ε σ
ρ σ ω ρ σ σ
ε σ
• •
•
• •
•
•
− +
−
−
=
−
−
−
=
L s
T T
T
o o
（2.26）
2.4有限變形之 Total Lagrangian Formulation
如圖 2-3 所示，一物體在時間為to之變形前平衡狀態，由座標系 統Xi表示，且其體積為Vo，表面積為So，微小體積為dVo，微小表面
積為dSo，no為dSo之單位法線向量分量。在時間為t之變形平衡狀態，
由座標系統xi表示，其體積為V，表面積為S，微小體積為dV ，微小 表面積為dS，n為dS之單位法線向量分量。若考慮作用在變形後dS面 上之微小表面積力dF，與dV體積上之微小體積力dFb，則單位表面力
f 及單位體積力b，則可由下列公式表示 dS
f
dFi = i （2.27） dV
b
dFbi = i （2.28） 如以變形前之座標系統為準，則單位面積之表面力 fo與單位體積 之體積力bo，則可由下列公式表示
o oi
i f dS
dF = （2.29）
o oi
bi b dV
dF = （2.30） 考慮表面力F 與體積力Fb間之平衡方程式，則
=0 + bi
i F
F （2.31） 積分（2.28）式及（2.30）式，可得
∫
∫
== V i Vo oi o
bi bdV b dV
F （2.32） 而由（2.29）式兩側積分得
S ji oj o
i s n dS
F
∫
o= （2.33）
利用Gauss 散度定理，（2.33）等號右側可寫為
V j o
ji
So ji oj o dV
X dS s
n s
∫
o∫
= ∂∂ （2.34）將（2.32）式與（2.34）式代入（2.31）式，可得 0
)
(∂ + =
∫
Vo ∂ j oi oji b dV
X
s （2.35）
由於dVo可為任意範圍，故上式得知
=0
∂ +
∂
j oi
ji b
X
s （2.36）
則上式以變形前座標為基準所推導出之平衡方程式，而對時間微分後 可得
=0
∂ +
∂• •
oi j
ji b
X
s （2.37）
（2.37）式其邊界條件分別為
ji oj
oi s n
f •
• = on Sot （2.38）
i ui
v = on • Sou （2.39） 在上兩式中，Sot表示具有表面力之表面，Sou為具有已知速度之表 面，foi為單位表面積之表面力，而 f•oi為其表面力變化率，vi為速度，
ui為位移。
將上述理論偶合入有限元素法時，必須經由虛功原理推導。假設 在物體之表面So上給予一假想之速度δvi，則表面力 foi所產生的虛功 率δW• 為
S oi i o o
S foi vidS f vdS
W
ot
o
∫
∫
=• =
δ δ
δ （2.40）
對於 Lagrange 應力，考慮邊界條件（2.38）式，可知 Lagrange 應力
率s•ji的虛功原理也是成立的，故可得
o i S ji oj o
S foi vidS s n vdS
o o
δ
δ
∫
∫
• = • （2.41）應用Gauss 散度定理，則
V o
j ji i o
V i j ji V o
j ji i o
i S ji oj
X dV s v
dV X v
s X dV
v dS s
v n s
o o
o o
∫
∫
∫
∫
∂ + ∂
∂
= ∂
∂
= ∂
• •
• •
) ( )
(
δ δ δ δ
（2.42）
結合（2.37）、（2.41）與（2.42）三式，整理得到虛功原理方程式為
o V oi i
o S oi i
V sji LijdVo f v dS b vdV
o o
o
δ δ
δ
∫ ∫
∫
• = • + • （2.43）其中
j
ij X i
L =∂∂δv δ ( )
由於（2.43）式係基於變形前之座標系統下推導而得，故稱為 Total Lagrangian Formulation （TLF）。
2.5 有限變形之 Update Lagrangian Formulation
上節所述之 TLF 理論係以to時刻為變形前之座標Xi為基準所建 立之虛功原理方程式。本節則討論每一變形後時刻t，以變形後之座 標xi為基準而加以考慮之情形。因此前節所述之座標Xj以xj取代，Vo 及So分別以V 及S置換，而 f•oi與b•oi則以 f•i與b•i取代，即可得 Updated Lagrangian Formulation（ULF）虛功原理方程式如下：
