塑性加工の温度解析
東京大学生産技術研究所
柳本
潤
1.はじめに
塑性加工において、温度に関連する問題は数多く存在する。例えば、
■熱間加工変形抵抗(流動応力)
■冷間加工変形抵抗(流動応力)
■冷間・熱間変形加工後のミクロ材質・結晶粒径・転位密度
■金型寿命・摩耗
■潤滑性能
等の問題は、被加工材の塑性変形のみならず、塑性変形に伴って発生する温度分布ない
しは温度の値によって影響を受ける。また塑性加工に関連する各種の技術課題は、塑性
変形・温度と、材料学的因子・化学的因子とが絡み合って存在している。
工具・被加工材界面
被加工材
熱伝導
塑性変形
工具
負荷・面圧
潤滑特性
表面酸化
軟化・回復・再結晶
熱伝導
熱伝導
被加工材の性質
工具材料軟化
図1 塑性加工に関連する技術課題と熱伝導との関係
図1に、塑性加工に関連するいくつかの技術課題と熱伝導との関係を記した。他の技術
課題を含め、およそ全ての課題に熱伝導と、これによりもたらされる温度が関わってい
るといっても過言ではない。課題毎に温度の影響度には差があるので、常にこれを意識
している必要があるわけではないが、過去の塑性加工理論の進展を振り返ってみると、
温度解析の重要性は一貫して高まっていると考えて良い。
ここでは、塑性加工の温度解析について、有限要素法の利用を前提とした事柄をまと
める。
- 176 -
2.熱伝導方程式
本章ではエネルギーの保存則より、熱伝導方程式の導出を行う1)。
2.1 エネルギー保存則
単位質量あたり物質が有する全エネルギーを E とする。 E は、
1
ui ui
□運動エネルギー
2
□内部エネルギー e
□位置エネルギー Ω
の総和である。さて、物質が占める領域の中に、基準体積(基準6面体) dV = dxdydz を
考えると、基準体積に含まれる物質の持つエネルギーの総和は、 ρEdxdydz であり、その
時間変化率
∂
(ρEdxdydz ) = ∂(ρE ) dxdydz はエネルギー保存則より、以下の通りに表され
∂t
∂t
る。
物質の持つエネルギーの総和 E の時間変化率=
単位時間あたりに面を介すことなく供給される熱量 Q&
+単位時間あたりに圧力・粘性力によって外部からされる仕事 W&
-単位時間あたりに対流によって外部に流出するエネルギー E&
-単位時間あたりに熱伝導によって外部に流失するエネルギー q&
2.2 熱伝導
熱伝導による流失は良く知られた Fourier の法則によって表される。Fourier の法則
は、「ある単位面積を介して伝わる単位時間あたりの熱量は温度勾配に比例する」とい
∂T
うことである。例えば x 軸に垂直な単位面積を介して x 方向のみに伝わる熱量は − κ
∂x
と表される。 κ は熱伝導率である。Fourier の法則を基に基準体積について q& を評価す
ると、次式が得られる。
⎧ ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞ ∂ ⎛ ∂T ⎞⎫
∂
⎟⎟ + ⎜ κ
q& = − ⎨ ⎜ κ
⎟ + ⎜⎜ κ
⎟⎬dxdydz = −
∂xi
⎩ ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠⎭
⎛ ∂T ⎞
⎜⎜ κ
⎟⎟
⎝ ∂xi ⎠
(1)
2.2 対流
対流は基準体積への物質の出入りであり、当然エネルギーの輸送を伴う。基準体積を
構成する座標軸より見た物質の速度を u : (u x , u y , u z ) とする。単位時間に輸送されるエネ
- 177 -
ルギーは ρE ⋅ u で表されるが、その変化量を基準体積について評価すると次式が得られ
る。
⎧ ∂ ( ρEu x ) ∂ (ρEu y ) ∂ ( ρEu z ) ⎫
∂ ( ρEu i )
+
+
E& = ⎨
dxdydz
⎬dxdydz =
∂y
∂z ⎭
∂xi
⎩ ∂x
z
(ρE ⋅ u )dxdz
y
(2)
∂ (ρE ⋅ u y ) ⎫
⎧
dy ⎬dxdz
⎨ ρE ⋅ u y +
∂
y
⎩
⎭
y
x
図2 対流(物質移動)により輸送されるエネルギー
図2は、対流(物質移動)により輸送されるエネルギーを、 y 軸方向について記したも
のである。
2.3 応力のなす仕事
応力(圧力・粘性力に起因する)のなす仕事は、以下の通りに導出できる。例えば x −
面について、作用している応力は σ xj 、物質移動は u j によって表されるのでこの面に作
用している応力によってなされる仕事は、 σ xj u j dydz と表される。 x + 面に作用している
∂ (σ xj u j ) ⎫
⎧
応力によってなされる仕事は、 ⎨σ xj u j +
dx ⎬dydz であるから、 x 方向の増減は
∂x
⎩
⎭
∂ (σ xj u j )
∂x
W& =
dxdydz 、同様にして y 方向、 z 方向の増減を考えると次式が得られる。
∂ (σ ij u j )
∂xi
∂u j ⎞
⎛ ∂σ ij
⎟dxdydz
dxdydz = ⎜⎜
u j + σ ij
∂xi ⎟⎠
⎝ ∂xi
(3)
さて平衡条件にある物体については、次式の Cauchy の第1運動法則が成立する。
ρu& i = ρ
Du i ∂σ ji
=
+ ρg i
Dt
∂x j
(4)
- 178 -
式(3)に式(4)を代入すると、以下の式が得られる。
∂u j ⎫
⎧ ⎛ Du i
⎞
− g i ⎟ + σ ij
W& = ⎨ ρu i ⎜
⎬dxdydz
∂xi ⎭
⎠
⎩ ⎝ Dt
(5)
ただし g は物体力である。
2.4 熱伝導方程式の導出
まず、式(2)を展開する。
⎛
∂ ( ρu i ) ⎞
∂E
⎟dxdydz
E& = ⎜⎜ ρu i
+E
∂xi
∂xi ⎟⎠
⎝
(6)
式(6)、(5)、(1)をエネルギー保存則に代入する。
∂u j
∂ (ρu i ) ∂
∂ (ρE )
∂E
∂ρ &
∂E
⎛ Du
⎞
=ρ
+E
= Q + ρu i ⎜ i − g i ⎟ + σ ij
− ρu i
+E
+
∂t
∂t
∂t
∂xi
∂xi
∂xi
∂xi
⎝ Dt
⎠
⎛ ∂T ⎞
⎜⎜ κ
⎟⎟
x
∂
i
⎝
⎠
移項すると、
ρ
∂u j
∂ ( ρu i ) &
∂E
∂E
∂ρ
∂
⎛ Du i
⎞
+ ρu i
+E
+E
= Q + ρu i ⎜
− g i ⎟ + σ ij
+
∂t
∂xi
∂t
∂xi
∂xi ∂xi
⎝ Dt
⎠
連続の式より、
ρ
⎛ ∂T ⎞
⎜⎜ κ
⎟⎟
⎝ ∂xi ⎠
∂ρ ∂ ( ρu i )
+
= 0 なので、
∂t
∂xi
∂u j
∂ ⎛ ∂T ⎞
DE &
⎞
⎛ Du i
⎜κ
⎟
+
= Q + ρu i ⎜
− g i ⎟ + σ ij
∂xi ∂xi ⎜⎝ ∂xi ⎟⎠
Dt
⎠
⎝ Dt
(7)
が得られる。
一方、 E =
1
u i u i + e + Ω であるので、これの物質時間微分を取り、位置エネルギーの
2
物質時間微分と物体力との関係に注意すると、
Du i DΩ De
Du i
DE
De
= ui
+
+
= ui
− ui g i +
Dt
Dt
Dt
Dt
Dt
Dt
が得られる。この式を(7)式に代入して整理すると、次式が得られる。
ρ
⎛ ∂T
DT
∂
De
∂T ⎞ &
⎟⎟ = Q +
+ ui
= ρc
= ρc⎜⎜
Dt
∂xi
Dt
∂xi ⎠
⎝ ∂t
∂u j
⎛ ∂T ⎞
⎜⎜ κ
⎟⎟ + σ ij
∂xi
⎝ ∂xi ⎠
- 179 -
(8)
ただし c は比熱である。(8)式は熱伝導方程式と呼ばれる。右辺第3項は、塑性仕事に伴
う内部発熱を表している。
2.5 観測枠と熱伝導方程式の表示の相違
(8)式は空間に固定された座標系についての熱伝導方程式である。基準体積を物体に固
定した場合、基準点(図2の原点)での速度 u はゼロになる(ただし速度勾配はゼロに
∂T
はならない)。従って(8)式のうち物質輸送項 ρcu i
が省略できるので、熱伝導方程式は
∂xi
以下の式により表される。
ρc
∂
∂T
= Q& +
∂xi
∂t
∂u j
⎛ ∂T ⎞
⎜⎜ κ
⎟⎟ + σ ij
∂xi
⎝ ∂xi ⎠
(9)
2.6 工具と被加工材が接触する部分についての境界条件
工具と被加工材など、異材が接触している部分の熱収支が、それぞれの部分の温度解
析を実行するにあたっての境界条件となる。境界条件を大きくわけると、以下の2通り
がある。
1)直接解析法
接触部に存在する熱抵抗層を直接モデリングする方法で、熱抵抗層の厚さ、熱伝導率、
比熱を条件値・物性値として必要とする方法。この方法では、連続する温度分布が解と
して得られる。
2)等価熱伝達係数を利用する方法
接触部を介して伝達する熱流束 q& ′ を、q& ′ = h(T1 − T2 ) として、物体1,2の温度差 T1 − T2
と等価熱伝達率 h により与える方法。物性値の代わりに等価熱伝達率を利用できるため
簡便ではあるが、解として得られる温度分布が不連続となるため、界面での温度の解析
精度には問題がある。
2.7 その他の境界条件
自由表面で利用される境界条件としては、冷却水への熱伝達、大気への熱伝達(強制
対流、自然対流)、輻射などがあり、条件に応じて使い分けがなされている。
- 180 -
3.有限要素法による熱伝導方程式の解法
本章では、熱伝導方程式(8)を有限要素法にて解析するための基礎式を説明する。
3.1 熱伝導方程式の重み付き残差法による弱形式表示
解くべき領域方程式である熱伝導方程式を、移項して再び記す。
⎛ ∂T
∂
∂T ⎞ &
⎟⎟ − Q −
+ ui
∂xi
∂xi ⎠
⎝ ∂t
ρc⎜⎜
∂u j
⎛ ∂T ⎞
⎜⎜ κ
⎟⎟ − σ ij
=0
∂xi
⎝ ∂xi ⎠
(8’)
スカラー関数によって表される(8’)式を、解析の対象とする領域全てについて満足させる
ことが本来必要である。この様な形式を、「微分方程式の強形式」と呼ぶ。これに対し
て、有限の積分領域の平均値として微分方程式を解こうとする考え方がある。その際の
形式を、「微分方程式の弱形式」と呼ぶ。有限要素法による(8’)式の解析は、まず(8’)式を、
重み付き残差法を利用しつつ弱形式に置き換えることに始まる。ちなみに、変形解析の
有限要素解析で利用される「仮想仕事の原理」は、平衡方程式(4)の弱形式である。
境界条件を満足するスカラー関数を w(xi ) とする。温度場 T (xi ) が領域 V において熱伝
導方程式を満足する時には、次式が成立する。
⎧ ⎛ ∂T
∂T ⎞
∂
w
∫V ⎨ ρc⎜⎜⎝ ∂t + ui ∂xi ⎟⎟⎠ − Q& − ∂xi
⎩
∂u j ⎫
⎛ ∂T ⎞
⎜⎜ κ
⎟⎟ − σ ij
⎬dV = 0
∂xi ⎭
⎝ ∂xi ⎠
(10)
式(10)にガウスの発散定理(式(11))を適用する。
∫ divVdv = ∫ V ⋅ ndA
V
S
(11)
∫ div(wV )dv = ∫ wV ⋅ ndA
V
Vi = κ
S
∂T
と置き、また、 div (wV ) = wdivV + grad (w ) ⋅ V であることを利用すれば、式(10)
∂xi
第4項は以下の通りに変形できる。
∂
∂ ⎛ ∂T ⎞
⎜⎜ κ
⎟⎟dV = ∫
∂xi
∂xi ⎠
i ⎝
V
∫ w ∂x
V
⎛
∂w ∂T
∂T
∂w ∂T
∂T ⎞
⎜⎜ wκ
⎟⎟dV − ∫ κ
dV
ni dS − ∫ κ
dV = ∫ wκ
∂xi ∂xi
∂xi
∂xi ∂xi
∂xi ⎠
⎝
V
S
V
- 181 -
Fourier の法則より、
∂T
∫ κw ∂x
S
i
ni dS = − ∫ wq&dS
S
であるので、これらの式を、式(10)に代入して整理すると、次式が得られる。
∫ w ρc
V
⎛ ∂u j & ⎞
⎛ ∂w ∂T
∂T ⎞
∂T
⎟⎟dV − ∫ w⎜⎜ σ ij
+ Q ⎟⎟dV + ∫ wq&dS = 0
+ wρcu i
dV + ∫ ⎜⎜ κ
∂xi ⎠
∂xi
∂xi ∂xi
∂t
⎝
⎠
V
S
V⎝
(12)
3.2 熱伝導方程式の有限要素表示
式(12)を、有限要素法により離散化する。領域を N 個の要素に分割し、要素 e の節点
温度を θ e と表す。要素内部の任意の点の温度 θ 、ならびに重み関数 w はそれぞれ内挿
{ }
関数を利用して、
θ = ⎣N ⎦{θ e }
(13)
{ }
(14)
w = ⎣W ⎦ w e
と表されるものとする。以後、式(12)に含まれる積分毎に、有限要素表示を行う。なお、
以後は指標記法を利用しないで表示する。
1)第1項
{w }
e T
⎧ ∂θ e ⎫
T
⋅ ρc ∫ ⎣W ⎦ ⎣N ⎦dV ⋅⎨
⎬
⎩ ∂t ⎭
V
{ } { } ⎣W ⎦
ただし、 w = ⎣W ⎦ w e = w e
T
(15-1)
T
であることと、比熱・密度および節点での値は体積積
分の対象とならないことを利用している。
2)第2項
{w } ⋅ κ ∫ ⎛⎜⎜ ∂ ⎣W∂x⎦
⎝
e T
V
T
T
T
∂ ⎣N ⎦ ∂ ⎣W ⎦ ∂ ⎣N ⎦ ∂ ⎣W ⎦ ∂ ⎣N ⎦ ⎞
⎟dV ⋅ θ e
+
+
∂x
∂y
∂y
∂z
∂z ⎟⎠
{ }
- 182 -
∂ N
∂ N
∂ N ⎞
T
T⎛
+ {w e } ⋅ ρc ∫ ⎣W ⎦ ⎜⎜ u x ⎣ ⎦ + u y ⎣ ⎦ + u z ⎣ ⎦ ⎟⎟dV ⋅ {θ e }
∂x
∂y
∂z ⎠
⎝
V
(15-2)
3)第3項
括弧内は全て内部発熱項に相当するのでこれらをまとめて Q& ′ と記す。
{ } ∫ ⎣W ⎦ Q& ′dV
T
T
− we
(15-3)
V
4)第4項
{w } ∫ ⎣W ⎦
e T
T
(15-4)
q& dS
S
この項は、要素境界を通過して流れる熱流束に対応している。
5)要素についての熱伝導方程式
式(15-1)~(15-4)をまとめて項をまとめると、以下の式が得られる。
{w } ⎛⎜⎜ [C ]⎧⎨ ∂∂θt
e T
⎩
⎝
e
{ } { } { }⎞⎟⎟ = 0
⎫
e
e
e
⎬ + [K ] θ − Q& ′ − f
⎭
(16)
⎠
式(16)が境界条件を満足する任意の重み関数について成立するためには、次式が成立し
なければならない。
[C ]⎧⎨ ∂θ
⎫
e
e
⎬ + [K ] θ − Q& ′ = f
⎩ ∂t ⎭
{ } { } { }
e
(17)
e
ただし、
[C ] = ρc ∫ ⎣W ⎦T ⎣N ⎦dV
(18-1)
V
⎛
[K ] = κ ∫ ⎜⎜ ∂ ⎣W ⎦
V
⎝
∂x
T
T
T
∂ ⎣N ⎦ ∂ ⎣W ⎦ ∂ ⎣N ⎦ ∂ ⎣W ⎦ ∂ ⎣N ⎦ ⎞
⎟dV
+
+
∂x
∂y
∂y
∂z
∂z ⎟⎠
∂ N
∂ N
∂ N ⎞
T⎛
+ ρc ∫ ⎣W ⎦ ⎜⎜ u x ⎣ ⎦ + u y ⎣ ⎦ + u z ⎣ ⎦ ⎟⎟dV
∂x
∂y
∂z ⎠
⎝
V
- 183 -
(18-2)
{Q& ′ } = ∫ ⎣W ⎦ Q& ′dV
(18-3)
{f } = − ∫ ⎣W ⎦
(18-4)
T
e
V
e
T
q&dS
S
である。 [C ], [K ] はそれぞれ熱容量マトリックス、熱伝導マトリックスと呼ばれ、(節点
数)×(節点数)のマトリックスとなる。式(18-3)は、内部発熱量を節点での値に換算
したものに相当する。また、式(18-4)では、要素境界を経て要素外部に流れ出す熱量を
節点での値に換算した値となり、熱流束ベクトルと呼ばれる。熱流束ベクトルは、物体
内部の境界面では隣接要素との熱収支より等しい値となることから打ち消しあい、物体
外表面ではこの面を通過して流入(流出)する熱量に等しくなる。このことを利用して、
考えている領域を構成する全ての要素について式(17)の和を取ることにより、領域全体
についての有限要素式が得られる。なお、物体外表面を通過して流入する熱流束は、2.6
項に述べた方法によって評価される。
3.2 時間積分と安定限界
式(17)にて表される熱伝導方程式は、時間刻みについての離散化が必要である。時刻 t
と時刻 t + ∆t における値を利用すれば、第1項については、
[C ]⎧⎨ ∂θ ⎫⎬ ≅
⎩ ∂t ⎭
1
[C ]{{θ (t + ∆t )} − {θ (t )}}
∆t
(19)
と近似できる。なお以後は、領域全体についても成立するので、要素についての添え字
e を省略して表示する。
熱伝導マトリックスに関係する項については、以下の近似を利用する。
[K ]{θ } ≅ α ⋅ [K ]{θ (t + ∆t )} + (1 − α ) ⋅ [K ]{θ (t )}
(20)
ただし α は、 0 ≤ α ≤ 1 の定数である。式(19)、式(20)を式(17)に代入すると、次式が得ら
れる。
{
}
1
1
⎛
⎞
⎛
⎞
⎜α ⋅ [K ] + [C ]⎟{θ (t + ∆t )} = ⎜ − (1 − α )[K ] + [C ]⎟{θ (t )} + Q& ′(t ) + { f (t )}
∆t
∆t
⎝
⎠
⎝
⎠
1) Full Implicit Scheme
α = 1.0 の場合には、式(21)は以下の通りとなる。
- 184 -
(21)
{
}
1
1
⎛
⎞
⎜ [K ] + [C ]⎟{θ (t + ∆t )} = [C ]{θ (t )} + Q& ′(t ) + { f (t )}
∆t
∆t
⎝
⎠
(22-1)
この方法は無条件安定である。
2) Full Explicit Scheme
α = 0.0 の場合には、以下の式が得られる。
{
}
1
[C ]{θ (t + ∆t )} = ⎛⎜ − [K ] + 1 [C ]⎞⎟{θ (t )} + Q& ′(t ) + { f (t )}
∆t
∆t
⎝
⎠
(22-2)
この方法では、左辺の係数に温度勾配(熱伝導)に関係する項が含まれていない。それ
故、時間ステップ ∆t の大きさが、要素のサイズ ∆x と比較して以下の大きさ以下になっ
ていないと、数値解が不安定になる。
∆t ≤
∆x 2 ρc
2 κ
(22-3)
3) Crank-Nicolson Scheme
α = 0.5 の場合には、以下の式が得られる。
{
}
1
1
⎛1
⎞
⎛ 1
⎞
⎜ [K ] + [C ]⎟{θ (t + ∆t )} = ⎜ − [K ] + [C ]⎟{θ (t )} + Q& ′(t ) + { f (t )}
∆t
∆t
⎝2
⎠
⎝ 2
⎠
(22-4)
この方法は無条件安定である。
3.3 重み関数
重み関数についての内挿関数 ⎣W ⎦ としては通常、形状関数 ⎣N ⎦ がそのまま利用される。
この方法を特に、ガラーキン法と呼ぶ。
これに対し、Eular 表示を利用して解析を行う場合、式(18-2)右辺第2項の物質移動
を表す項(移流項)が存在する。この項が大きい場合には、本来無条件安定であるはず
の時間積分法を利用しても、解が不安定になる場合がある。流れ方向速度を u 、流れ方
向の要素長さを l とすれば、ガラーキン法を利用した場合の安定限界は、
l≤
2κ
ρcu
(23)
- 185 -
で表される1)。この要素長さは、圧延に相当する条件では数ミクロンとなり、要素分割
数と計算時間の増大を招く。
これに対して、重み関数の内挿関数 ⎣W ⎦ として、流れ方向に上流の要素の形状関数を
も含めて表し、安定性を増す方法がある。1次元問題について流れ方向上流側の節点を
1、下流側を2とする。下流側に正の方向を取ると、1次の内挿関数を利用する場合の
形状関数 N は、
N1 =
1
(1 − ξ ), N 2 = 1 (1 + ξ )
2
2
(24)
と表される。ここで、重み関数の内挿関数 ⎣W ⎦ として、以下の関数を利用する方法を、
ペトロフ・ガラーキン法と呼ぶ。
W1 = N1 + βϕ (ξ ) =
1
(1 − ξ ) + βϕ (ξ ), W2 = N 2 − βϕ (ξ ) = 1 (1 + ξ ) − βφ (ξ )
2
2
(25)
ただし β は正の定数、φ (ξ ) は両端 ξ = ±1.0 で φ = 0.0 、区間内で正の値を取る関数である。
この場合、式(23)以外に β ≥ 1.0 で解が安定するので、要素分割・増分時間についての制
約が緩和される。この方法は、差分法でしばしば用いられる風上差分法と対応しており、
有限要素法の場合には特に上流有限要素法と呼ばれる。
3.4 その他の注意事項
有限要素法の基本的な構成(スキーム)と安定限界については以上に述べたとおりで
ある。塑性加工の有限要素解析では工具と材料間の接触熱伝導を取り扱う場合が多いが、
この部分での放線方向温度勾配は非常に大きく、そのため、放線方向の要素サイズは十
分に小さく(場合によっては数ミクロン程度)しておかないと、解析精度が損なわれる
場合がある。
境界条件としては、2.6 節・2.7 節に述べたとおりの条件を利用することになる。ただ
し、境界条件に含まれる係数の値を適正としておかないと解析精度が損なわれるので、
その選定には十分な注意が必要であるが、残念ながら信頼すべきデータが少ない。
参考文献
1)若松英士:東京大学博士論文,(2000)
- 186 -
ダウンロード

塑性加工の温度解析(PDF) - 柳本研究室