行政院國家科學委員會專題研究計畫 成果報告
子計畫一：金屬板材突緣引伸成形製程之分析
計畫類別： 整合型計畫
計畫編號： NSC91-2212-E-011-025-
執行期間： 91 年 08 月 01 日至 92 年 07 月 31 日 執行單位： 國立臺灣科技大學機械工程系
計畫主持人： 向四海 共同主持人： 黃佑民
計畫參與人員： 郭哲良、陳聰嘉、周村憲
報告類型： 完整報告
處理方式： 本計畫可公開查詢
中 華 民 國 92 年 10 月 23 日
行政院國家科學委員會補助專題研究計畫 成 果 報 告
□期中進度報告
智慧型金屬板材成形 CAE 軟體之開發與整合研究(3/3) ---子計畫一 金屬板材突緣引伸成形製程之分析
計畫類別：□ 個別型計畫 整合型計畫 計畫編號：NSC 89－2212－E－011－038－
NSC 90－2212－E－011－014－
NSC 91－2212－E－011－025－
執行期間： 89 年 08 月 01 日 至 92 年 07 月 31 日
計畫主持人：向 四 海 共同主持人：黃 佑 民
計畫參與人員：郭哲良、陳聰嘉、周村憲
成果報告類型(依經費核定清單規定繳交)：□精簡報告 完整報告
本成果報告包括以下應繳交之附件：
□赴國外出差或研習心得報告一份
□赴大陸地區出差或研習心得報告一份
□出席國際學術會議心得報告及發表之論文各一份
□國際合作研究計畫國外研究報告書一份
處理方式：除產學合作研究計畫、提升產業技術及人才培育研究計畫、
列管計畫及下列情形者外，得立即公開查詢
□涉及專利或其他智慧財產權，□一年□二年後可公開查詢
執行單位：國立台灣科技大學 機械工程學系
中 華 民 國 九十二 年 十 月 三十一 日
摘 要
本計畫為整合型計畫之子計畫一，以開發金屬板材成形三維解析 CAE 軟體有關前處理(pre-processor)的部份為主，使用 I-DEAS 套裝軟 體建立工具(包括沖頭、模具與壓料板)與加工件等有限元素之三維模 式，並將其節點座標值、元素號碼及板材邊界條件等數據檔輸出並轉 換成金屬板材成形三維解析 CAE 軟體能接受之輸入格式，並撰寫 CAE 軟體有關邊界條件與材料機械性質輸入部份之程式等以建立完 整的 CAE 軟體之前處理部份，並與其它子計畫整合成完整的金屬板 材成形三維解析 CAE 軟體，且將其應用於金屬板材突緣引伸成形之 分析。
關鍵詞：金屬成形、CAE 軟體前處理、突緣引伸
Abstract
The main purpose of this sub proposal is to develop part of preprocessor of CAE software of metal forming processes. Considering to use IDEAS package to build up 3D tools (including punch, die and blank holder) as well as sheet metal. The CAE computer code will be established to analyze the 3D flanging drawing process of sheet metal.
Keyword: metal forming, pre-processor, flanging drawing.
目 錄
摘 要... I ABSTRACT... II
一、前 言...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 邊界條件 ...37
4.2 材料參數 ...38
4.3 數值模擬與分析 ...38
五、參考文獻 ...48
六、計劃成果自評 ...50
圖表索引
圖 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 帽型引伸加工之工件與邊界條件 ...41
圖 4-2 帽型引伸加工之沖頭 ...41
圖 4-3 帽型引伸加工之模具 ...42
圖 4-4 帽型引伸加工之壓料板 ...42
表 4-1 帽型引伸成形分析，模具與料片之網格分割數據...43
圖 4-5(A) 帽型引伸衝程 20MM厚度圖...43
圖 4-5(B) 帽型引伸衝程 40MM厚度圖...44
圖 4-6(A) 帽型引伸衝程 20MM應力圖...44
圖 4-6(B) 帽型引伸衝程 40MM應力圖...45
圖 4-7(A) 帽型引伸衝程 20MM應變圖...45
圖 4-7(B) 帽型引伸衝程 40MM應變圖...46
圖 4-8(A) 帽型引伸衝程 20MM反力圖...46 圖 4-8(B) 帽型引伸衝程 40MM反力圖...47
一、前 言
國內汽車工業、家電工業等產業，其產品製造的過程中，金屬板 材成形的加工製程所佔比率相當高，且隨著消費者導向的時代來臨，
新規範的產品開發、設計越來越多，所開發的產品也朝著少量而多樣 化的方向發展，每一新產品之生命週期越來越短。以汽車為例，每隔 幾年各汽車製造廠為了滿足顧客的需求，必須要將汽車改款，使同一 款式的產品週期越來越短，因此，以往靠經驗與感覺來設計→試作→
試驗→修正設計→再試作→再試驗等步驟之傳統設計方法浪費工 時，降低產品的競爭力。藉由 CAE 軟體可將設計→試作→試驗→修 正設計→再試作→再試驗等反覆的步驟，利用電腦支援而縮短開發及 設計時程因而降低成本，亦可提升新產品的品質與設計最佳化。
彎曲成形、拉伸成形及深引伸成形為板金成形之基本加工方法，
而本子計畫所探討突緣(flanging drawing)成形係屬於引伸成形加工之 一種，其成形加工過程是將料片放置於沖模上，而後放置壓料板將料 片周緣夾固於沖模與壓料板之間，藉由沖頭對沖模之相對位移作引伸 沖壓成形，料片空孔會隨著沖頭衝程的增加而逐漸擴大，此外隨著圓 周不斷擴張而造成厚度逐漸變薄，甚至達到成形極限，此時空孔圓周 厚度將急速薄化，並造成頸縮現象，而後發生破裂。突緣成形加工為 一重要板金加工技術，廣泛應用於工業零組件之製造上。因此，若能
預先模擬突緣成形之製程，將可作為實際加工時的參考，並取代一般 傳統模具工廠中所採用的試誤法，以得到精確的產品。
金屬板金成形之數值模擬為一複雜且困難的工作，因成形之過程 為非線性問題，而造成此非線性的因素大致可分為三個部分：第一為 工件變形的大位移、大旋轉與大應變所造成幾何上的非線性行為；第 二為因變形而產生的材料非線性構成關係；第三為成形製程中，料片 與模具接觸界面之幾何外形的函數之非線性關係。上述的各種非線性 問題，皆增加了金屬板金成形製程模擬的複雜性。本研究之方法乃結 合有限變形理論 ULF(updated Lagrangian formulation) 的觀念，採用 Cauchy 應力的 Jaumann rate 關係式，並忽略元素的密度或體積之變化 下，建立一個 incremental updated Lagrangian formulation 之彈塑性大 變形、大應變有限元素解析程式。本程式係以沖頭之增量位移為每一 變形增量步驟之起始增量值，文中採用 Yamada 的 rmin方法來判斷材 料內任一元素的彈塑性狀態變換問題，並將此方法擴展至料片上節點 之接觸、分離的判斷，摩擦方向的轉換及有關最大應變增量、最大旋 轉增量的線性化處理，以決定板金成形時每一增量的位移、應變、應 力與除荷後料片之回彈現象與工件最後形狀。
近年來，仍有許多研究學者將有限元素法應用於分析板金成形問 題上，且電腦的快速發展，使得其運算能力與速度大幅提昇，在處理
複雜且高階非線性問題時可縮短其計算的時間，而這些研究成果可提 供工業界作為板金成形實際加工之參考依據。而與本研究相關文獻有 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)及線 性內插法所描述的一般曲面的剛－黏塑性有限元素分析程式所完 成，並驗證此方法較以往所用的三維有限元素分析法更有效率。
Zhiheng 與 Zheng [10]探討中空圓形料片在引伸成形之趨勢，並採用 數學模式找出中空板材在引伸成形時空孔沒變形的情況。Makinouchi 與 Kawka [11]研究有限元素模擬模鍛(stamping)模具之設計，藉由數 值模擬來預測與估計板金成形過程中的缺陷以替代傳統的試誤法，其 理論背景與數值分析是由有限元素程式完成。Kawka 與 Makinouchi [12]亦提出退化殼元素的構想在靜態顯函有限元素程式中執行，主要 探討各種不同的積分法則比較，分別為完全積分(FI)、減化積分(RI)、
選擇減化積分(SRI)、假設應變場(ASF)與穩定矩陣(SM)法則等。Leu [13]以彈塑性有限元素法探討孔突緣成形。Lin 與 Tseng[14]為了解決 金屬成形在彈塑性區間的接觸問題，而應用一表面彈簧元素並將其元 素之對稱性剛性矩陣疊加於整體剛性矩陣中，組成有限元素之分析方 法，經由實驗結果與其數值分析結果之驗證，證明其方法可準確地模 擬實驗結果。Kumagai、Saiki 與 Meng[15]選用 Al/Cu 雙層板金以研究 具有引薄作用之突緣成形，利用 45°錐角之錐形沖頭作突緣成形實 驗，並將 FEM 模擬結果與實驗結果相比較，其顯示剛塑性有限元素 法可有效地預測雙金屬板金在突緣成形中之成形製程及最終形狀。
Huang 與 Chien [16]在軸對稱條件下，進行孔突緣成形極限之正交異 向性分析，經由數值模擬與其實驗結果相比較後，當初始孔徑小於 19mm 時，其料片於不同之初始孔徑與成形後不發生頸縮現象時之擴 孔直徑成線性關係，而當初始孔徑大於 19mm 時，其擴孔直徑與沖頭 直徑相同，此關係式可提供孔凸緣成形加工用模具設計之參考。Huang 與 Chien [17]亦探討孔突緣成形時沖頭的幾何外形與成形極限之關 係，經由數值模擬與其實驗結果相比較後，得知沖頭圓弧角半徑大小 對孔凸緣之成形極限的影響並不是很明顯，然而其圓弧角半徑大小會 影響其最大沖頭負荷。
本計畫為整合型計畫，其總體目標在於開發一套分析三維金屬板 材成形製程分析之有限元素 CAE 軟體，並於完成後能推廣且應用於 工業界。其各子計畫間之整體相關架構及研究目標，如圖 1-1 所示。
本計畫之目的，乃應用彈塑性增量理論有限元素法發展一套模擬金屬 板材突緣引伸成形製程的三維數值模擬程式，並考慮材料異向性對突 緣引伸成形製程之影響。藉由使用 I-DEAS 套裝軟體建立工具(包括沖 頭、模具與壓料板)與加工件等有限元素之三維模式，轉成金屬板材 成形三維解析 CAE 軟體能接受之輸入格式，並與其他子計畫所發展 的程式進行整合與測試。
智慧型金屬板材成形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 k
ij 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 k
ij 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 dV
dV
j
i = =
∂∂
= det[ ]
0
（2.8）
其中 J 為[J 的行列式值。且由質量不變定律可得 ]
ρ ρ0
=
J （2.9）
其中ρo與ρ分別為變形前與變形後之密度。
2.3 有限變形之應力與應力率
在處理有限變形之彈塑性問題時，應力則須給予明確的定義，一 般常使用之應力有 Lagrange 應力、Kirchhoff 應力和 Cauchy 應力三 種。本節對此三種應力之定義加以說明，並推導在 ULF 下應力率之 間的關係。如圖 2-2 所示，左側為變形前之狀態，右側為變形後之狀 態。
（1）Lagrange 應力(sij)
} { ] [ } {
) ;
(
o T i
o oj ji L
oi s n dS dF dF s dS
dF = = = （2.10）
其中sji為 Lagrange 應力張量，具有非對稱形式，以變形前之面 積dS 為基準，其中o n 為oj dS 之單位法線向量之分量。 o
（2）Kirchhoff 應力(tij)
} { ] [ ] [ } {
) ;
(
o T i T
o oj ji K
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 oj i 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）式對時間的導數，可得到各種應力的變化率 ]
][
[ ] ][
[ ] [
•
•
• = t J + t J
s （2.17）
] [ ] [ ] [ ] [ ] [ ] [ ] [
•
•
•
• = J G T σ + J G T σ + J G T σ
s （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]
[ ]
[ ][ ]] ][
[ ] [ ] ][
[ ] [ ] [
•
•
• + +
= T T σ T T T σ T T T σ T
σo （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
j
k 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）式，可得
∫
∫ =
= Vo oi o
V i
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 o
ji b dV
X
s （2.35）
由於dVo可為任意範圍，故上式得知
=0
∂ +
∂
oi j
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為單位表面積之表面力，而 foi
• 為其表面力變化率，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 應力
率sji
• 的虛功原理也是成立的，故可得
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 i
ij X
L =∂∂δv δ ( )
由於（2.43）式係基於變形前之座標系統下推導而得，故稱為 Total Lagrangian Formulation （TLF）。
2.5 有限變形之 Update Lagrangian Formulation
上節所述之 TLF 理論係以to時刻為變形前之座標Xi為基準所建 立之虛功原理方程式。本節則討論每一變形後時刻t，以變形後之座 標xi為基準而加以考慮之情形。因此前節所述之座標Xj以xj取代，Vo 及So分別以V 及S置換，而 foi
• 與boi
• 則以 fi
• 與bi
• 取代，即可得 Updated
Lagrangian Formulation（ULF）虛功原理方程式如下：
dV v b dS
v f dV
L
s i
V i
S i i
V jiδ ij ∫ δ ∫ δ
∫ • = • + • （2.44）
由於 ULF 係以 t 時刻之變形後狀態為參考狀態，並於微小 dt 之時刻 後，在更新其變形狀態為參考狀態。因此，對於邊界之描述與處理較 為簡易，而不像 TLF 必須時時刻刻參考其初始狀態，故 ULF 為多數 研究人員所採用。本研究亦以（2.44）式為基礎，來建立有限元素法 之統制方程式。在 ULF 時，[t]=[σ]，可知
] ][
[ ] ][
[ ] ][
[ ] ][
[ ] [ ] [
•
•
• − + + −
= t L σ ε σ σ L σ ε to
（2.45）
代入[t•]=[s•]+[L][σ]之關係，得到 L T
t
s] [ ] [ ][ ] [ ][ ] [ ][ ]
[• = o − ε• σ − σ ε• + σ （2.46）
或
jk kj ik
ik ik ij kj
ij t L
s• =o −σ ε• −σ ε• +σ
結合toij =toji，σij =σji，toij =σoij及 ij ji
•
• =ε
ε 的對稱特性，可導出下式
∫V(σoij−2σikε•kj)δ ε•ijdV +∫VσjkLikδLijdV =∫S •fiδvidS +∫Vb•iδvidV （2.47）
其中σo ij =σ•ij−ωikσkj +σikωkj。
2.6 材料之彈塑性構成關係式
在考慮彈性負荷、塑性除荷或必須考慮殘留應力以及除荷工件之 最後形狀時，彈性變形行為之因素必須加以考慮。塑性變形的特性之 一是受變形歷程影響，一般在解析材料之塑性變形時，皆是採用各瞬
