長岡工業高等専門学校研究紀要

論第50巻（2014）

文

仮想領域法を用いた有限要素解析によるフライス加工時における工具切れ刃先端近傍の間接的温度推定法に関する提案－ エンドミルの径と工具切れ刃先端近傍温度の関係に関する考察－

倉橋 貴彦1・羽根川 英2・井伊 優生3・杉澤 元次郎4・井山 徹郎5・近藤 俊美6
1 機械工学科（Department of Mechanical Engineering, Nagaoka National College of Technology）
2 機械工学科（Department of Mechanical Engineering, Nagaoka National College of Technology）
3 機械工学科（Department of Mechanical Engineering, Nagaoka National College of Technology）
4 機械工学科（Department of Mechanical Engineering, Nagaoka National College of Technology）
5 機械工学科（Department of Mechanical Engineering, Nagaoka National College of Technology）
6 機械工学科（Department of Mechanical Engineering, Nagaoka National College of Technology）（2014年度より：1, 2長岡技術科学大学

機械系）(現：3株式会社クラレ)

Proposal of Indirect Estimation Method with Respect to Temperaturenear Tip of Cutting Edge of Tool in Milling
Based on Fictitious Domain and Finite Element Methods- Remarks for Relationship between Diameter of Endmill and Temperaturenear Tip of Cutting Edge of Tool Takahiko KURAHASHI1, Suguru HANEGAWA2, Masaki II3, Motojiro SUGISAWA4,Tetsuro IYAMA5 and Toshimi KONDO6Abstract
In this study, an analysis of temperature distribution on material surface in milling is carried out based onfictitious domain and finite element methods. To know characteristics of tool wear, it is important toknow temperature at the tip of cutting edge in milling. Therefore, temperature near tip of cutting edge inmilling is indirectly estimated such that the observed temperature on the material surface is close to thetemperature obtained by numerical simulation. In case of numerical simulation considering the rotationand the movement of tool, we have to solve the moving boundary problem, and, in general, the adaptive
remeshing method is employed. On the other hand, if the fictitious domain method is employed, it is notnecessary to apply remeshing technique. In this study, the temperature near tip of cutting edge isestimated based on the results of finite element analysis using the fictitious domain method.
Key Words : fictitious domain method, finite element method, moving boundary problem, milling−65−

倉橋貴彦 ・ 羽根川

英 ・ 井伊優生 ・ 杉澤元次郎 ・ 井山徹郎 ・ 近藤俊美

１． はじめにフライス加工時における工具温度に着目した研究は数多く行われている 1)-3)．切削時の刃先温度は工具寿命の低下にも繋がるため，高速切削下や連続加工時間が増加する場合を対象とした刃先温度管理システムの構築 4)や，低熱伝導性の材料加工のための工具温度管理システムに関する研究が行われている5)．実際の切削時に，工具先端に熱電対を付け温度計測を行うことが可能であれば実測も行うことができるが，工具が回転体であるため工具先端近傍の温(a) 全体領域図度を計測することは困難である．そこで，本論文では，数値シミュレーションの結果から間接的に工具先端近傍温度の推定を行う方法について提案する．
フライス加工時における工具や材料の温度をシミュレーションにより求める場合，回転物体と移動物体を考慮した解析を行う必要がある．このような工具の解析を行う場合は，あらかじめ用意したメッシュの節点上を熱源の点が移動するものとし，節点上において温度を直接与える場合がほとんどである．しかし，工具先端の描く軌跡上に必ずしも節点が存在しているとは限らないため，要素内に工具の先端の(b) 拡大図点が存在する場合においても，切削対象内の温度の図-1 仮想領域法における有限要素メッシュの例

拡がり方を適切に再現できる様に計算が行える方が望ましいと考えられる．そのため，切削対象となる素法に基づく解析結果と切削時における材料表面温工作物と工具の領域を分け，２つの領域を使用した度の計測値の比較を行う．また，工具の摩耗特性等，熱伝導解析を考え，仮想領域法 6),7),8)を使用した有限工具を管理する上で，切削時における工具の先端近要素法に基づく熱伝導解析について検討を行う．

傍の温度を知る事は非常に重要なことである．そこ

仮想領域法とは，解析を行う全体領域（本論文でで，本研究では被削材表面において計測した温度履は切削対象の工作物）以外に副領域（本論文では工歴とシミュレーション結果が合うような，工具先端具）を設けて，副領域を全体領域に重ねて解析を行近傍における切削時の温度を間接的に推定する方法

う手法である．図-1の(a)に仮想領域法において用について検討する．解析モデルは板厚の薄いアルミいるメッシュの全体図を，(b)に拡大図を示す．移ニウム合金板を対象とし，厚み方向における熱伝導動物体問題においては要素の再生成を行う必要なくに比べ平面方向の熱伝導が支配的であるとし，平面解析を行うことが可能であり，この点が通常の解析二次元モデルを対象とする．

手法と異なる点であり，利点でもある．また，仮想領域法を使用した研究としては，主に流体中における粒子の落下問題を対象とした研究が行われており，実験時に見られる二つの粒子の落下時の吸引や衝突

２．支配方程式の離散化および背景領域内の各要素に対する副領域の節点探索の現象を良く再現できているとの報告もある 9)．仮想領域法を使用した研究は流れ場を対象のとしたも図-2に示す計算領域の背景領域Ω（二次元平面領のがほとんどであるため，本論文では温度場に対す域）に対して，式(1)に示す熱伝導方程式に対してる仮想領域法の適用を考える．具体的な内容として仮想領域法および有限要素法を適用する．

は，工具先端より温度を与え回転移動する工具を含む熱伝導解析に対して仮想領域法を適用し，有限要−66−仮想領域法を用いた有限要素解析によるフライス加工時における工具切れ刃先端近傍の間接的温度推定法に関する提案

ここに，λはラグランジュの未定乗数を示す．式(5)の各式の両辺に重み関数φ*，λ*を乗じて，背景領域における要素領域Ωeと，副領域における要素領域ωeにおいて積分をすると，重み付き残差方程式は式(6)のように得られる．
 n 1 n 1
 n 1 e   d  tK e   ,ii d  e  d n    d in e  n 1  d  0 in  e  

図-2 計算領域のイメージ図(1)  K,ii  0 in (6)また，式(6)の第1式の左辺第2項にグリーンの定理を適用すると，式(7)のようになる．

ここに，φは温度，Kは熱拡散率を示す．
境界は，第1種境界Γd（温度固定境界）と第2種境界Γn（熱伝達境界）に分けられ，式(2)   d  tK    d    d   d  t   q d in 
 n 1eのように書くことができる．

e n(2)

 n 1e n 1e  ˆ on dh    on n
q 
c


 n 1, i ,i(7)e

ここで，領域Ωeにおける積分に際して，補間関数を式(8)のように定義する．
T 





  N11  N 2 2  N 33  N c
T  N11  N 2 2  N 33  N c

ここに，qは熱流束，hは熱伝達率，ρは領域Ωの材料の密度，cは領域Ωの材料の比熱，φ∞は領域Ωの(8)外部の雰囲気温度を示す．ここで，式(1)に完全陰解法を適用し時間方向に離散化すると，式(3)のよここに，Niは三角形一次の形状関数を示し，添え字うに書くことができる．の1,2,3は各頂点の番号を示す．また，領域ωeにおける積分に際して，補間関数を式(9)のように定義す n1  tK,nii1   nin (3)る．

   x  x  a    x  x  a 

また，副領域（仮想領域）ωにおいては，温度φαが既知であるものとし，式(4)のように与える．
(4)

   a in (9)ここに，δはデルタ関数を示している．式(9)を使用
して，式(6)，(7)における領域ωeの積分を行うと式(10)，(11)のように書くことができる．

ここで，式(3)に副領域ωの温度の条件を考慮するために，ラグランジュの未定乗数法を導入すると，式(5)のように書くことができる．
 n1  tK,nii1  n1   n n1
   0 in in (5)

−67−

e

e

  d     d Nd (10)

i ii 1NdNd

     i 1


i ii 1

ii(11)

倉橋貴彦 ・ 羽根川

英 ・ 井伊優生 ・ 杉澤元次郎 ・ 井山徹郎 ・ 近藤俊美

ここに，Ndは対象とする要素Ωe内に存在する副領域

x3, y3 の節点の数を示す．ここで，φ*，φは式(8)に示す補間関数により表される．今，図-3のように，対象とする要素Ωe内に副領域の節点が1点ある状態を想定し，有限要素方程式を導く．ここに，副領域の節点座標はa=(xa,ya)とする．結果として，有限要素方

x1, y1 e

程式は式(12)のように書くことができる．本検討で

eは，解析領域全体において重ね合わせた有限要素方程式の解法に，共役勾配法を適用する．

M 0 c   0 0 r n1

M 0 c  
 0 0 r  Δt H B T  c  n 1 0  r  Bn  Δt t   f n a  xa , ya x2 , y2 (12)図-3 背景領域と副領域の重なった状態の例

Start

ここに，式(12)における各係数行列は式(13)のよう

Input of computational conditionsに書くことができる．
M  NNT deH  K N, x NT, x d  N, y NT, y d  eB   N
 1 x  xa , y  y a N 2 x  xa , y  y a N 3 x  xa , y  y a 




T
c  1 2 3r   x  xa , y  yat  qNdef   x  xa , y  yaComputation of gradient of shape function {N,x},{N,y} for each elements in domain ΩComputation of coefficient matrix [M]+Δt [H]in finite element equation for domain Ω

Time loop ：istep = 1～maximum number of time stepSearch of nodes of tool in domain Ω, and computation of coefficient matrix [B]Computation of finite element equation by conjugate gradient method
Update of time step：{c}(n) = {c}(n+1)
Update of location of tool by movement and rotation

(13)

Yes

istep ＜ maximum number of time stepNoEnd

図-4 仮想領域法に基づく有限要素解析のフローチャート

D1 

また，仮想領域法を用いた解析では背景領域の各要素における副領域の各節点を探索する必要がある．

D3 

ここで，図-3のように背景領域のある要素の三点の座標を(x1,y1)，(x2,y2)，(x3,y3)とし，その要素内に副x1  x ay1  y ax2  xa

y2  yax3  x ay3  y ax1  x ay1  y a,D2 x2  xa
x3  x ay2  ya
y3  y a,(14),領域の節点（節点座標(xa,ya)）が存在するものとする．この場合，式(14)に示す行列式を計算すると
D1>0，D2>0，D3>0となり，式(14)に示す行列式が全３．仮想領域法および有限要素法によるフライス加工時の非定常熱伝導解析て正の値となった場合に背景領域内の要素において副領域の節点があることを判定できる．ここに，フ3.1 エンドミル径が32mmの場合に対する検討

ライス加工時の工作物表面の温度分布を算出する計算の流れを図-4に示す．

本節では，フライス加工時の非定常熱伝導解析の例を示す．加工対象の材料をアルミニウム合金(A5052)（100mm×100mm×8mm）とし，本検討では，切込み深さ2mm，送り速度を2mm/s， 直径D=32mmの工具（図-5（(a)は横断面，(b)は縦断面））を使−68−

仮想領域法を用いた有限要素解析によるフライス加工時における工具切れ刃先端近傍の間接的温度推定法に関する提案用意する．矩形領域を領域Ωとし，図-7の左側に示
す領域は工具の断面形状を表したメッシュであり，副領域ωとする．また本検討では，矩形領域の周辺境界は熱伝達境界とし，雰囲気温度φ∞ は21.6℃，熱伝達率hは4W/(m2・℃)，アルミニウム合金の密度
ρと比熱cは2700kg/m3, 899J/(kg・℃)とする 10)．その他の計算条件は，表-1に整理する．
(a) 水平断面

また，同条件下において，工作物表面において熱画像をサーモグラフィーにより撮影した画像を図-8に示す．工具先端近傍の点における温度を測れる解像度は有していないため，画像より判別することは困難であるが，この熱画像より切削前半円において温度が高くなっていることは確認できる．この結果を踏まえ，図-9に示すように，切削している点のみに温度を与える条件とするため，図-9に示す角度θ

(b) 垂直断面が-90°<θ<90°の範囲の際に温度φαを切削先端の点に

図-5 エンドミルの写真（径：32mm）おいて与え，それ以外の場合はφαの条件を与えないものとする．また切削厚さの変化による工具刃先先端における温度の変化はないものとし，切削時においては一定温度とする．

Rotation direction
Movement direction
Heat transfer boundary

図-6 計測点の位置用し，回転数を1240min-1 と設定する．実際の加工では，工具は回転のみで加工対象となる工作物を送ることで切削が行われるが，解析においては，工作物（背景領域）を固定し，工具（副領域）を移動させるものとする．工具により切削する際の温度の計

Heat source points(a)全体図

測位置を図-6に示す．本検討では，切削加工時におけるアルミニウム合金表面の各計測点における計測温度と解析温度が良好に一致するシミュレーションを行うことを目的とする．また解析条件として，切れ刃の軸方向に対する熱源位置の変化，また切削時における溝の発生は考慮せず，上平面から見た板表面における温度分布の解析を行う．本解析においては，工具先端における温度は既知量として取り扱う
ため，工具側の熱伝導は考慮せず，アルミニウム合金側のみの熱伝導解析を行う．
仮想領域法を用いるため，図-7（(a)は全体図，(b)は拡大図）のように2種類の有限要素メッシュを−69−(b) 拡大図
図-7 有限要素メッシュと計算条件

倉橋貴彦 ・ 羽根川

英 ・ 井伊優生 ・ 杉澤元次郎 ・ 井山徹郎 ・ 近藤俊美

ここに，Nobsは5s置きに計測した計測値の数を示し

表-1 計算条件

Time increment Δt,sNumber of time step
Number of nodes in Ω
Number of elements in Ω
Number of nodes in ω
Number of elements in ω
Initial temperature in Ω, ℃
Thermal diffusivity K, mm2/s
Rotation speed of endmill, min-1
Feed speed, mm/s

0.001
40000
10201
20000
933
1768
21.5
99.75
1240
2.0ており，φcomp.A(i) ，φcomp.B(i) は解析により求まった温度，φobs.A(i) ，φobs.B(i) は計測により求まった温度を示している．計測温度は，計測点Point A, B（図6）において，熱電対（株式会社チノー製 熱電対(MC3000)，温度分解能：0.1℃）により収録周期を5sとし工作物表面において計測したデータである．
本検討では，計測点で測定した熱電対による温度の分解能が0.1℃であり，有効桁数も考え0.5℃刻みで検討を行う．工具切れ刃先端における温度φαに対して，残差二乗和を整理した結果を図-10に示す．

Endmill

結果として，φαの温度を38.2℃とした場合に，残差Tip of endmill二乗和の計測結果の値が最小となる結果となった．
φαの温度を38.2℃とした場合の計測点Point A,Bにおける解析結果と計測温度の比較を図-11に，また時刻T=10sごとにおける被削材の表面温度分布およびエンドミルの回転移動位置を図-12に整理する．切削時に工具切れ刃先端における温度φαを与える様30℃にしているため，温度分布は工具の先端が描く軌跡

Aluminum alloy plate

拡散し被削材全体に拡がっている様子を確認できる．

Boundary of aluminum alloy plate

また，工具は中心線上を回転しながら移動している

19℃

上が最も大きな値を示し，工具切れ刃先端の温度が

図-8 サーモグラフィによる T=20sec における材料表面ため，上下において対称な温度分布とはならないことも確認できる．の熱画像

3.2 エンドミル径が16mm,25mmの場合に対する検討
θ

前節の検討に基づき，加工対象の材料をアルミニウム合金(A5052)（100mm×100mm×8mm）とし，切θ=90°

45°

0°

－45°

込み深さ2mm，送り速度を2mm/sとした場合におい－90°：Temperature “φ a” is not given.て，エンドミルの径とエンドミル先端近傍温度の関：Temperature “φ a” is given.

係に関する考察を行う．径が16mm，25mmのエンTemperature “φ a” is given during cut of material,i.e., temperature “φ a” is given in -90° < θ ≦ 90°.

ドミルに対して，前節と同様の検討を行い，エンドミル先端近傍温度を推定し，各エンドミルの径とエ

図-9 副領域 ω における温度 g の与え方に対するンドミル先端近傍温度の関係から，エンドミル先端イメージ図近傍温度の推定式を誘導する．解析条件は表-1に従

数値実験として，工具刃先先端における温度φαうことにし，各径の場合に対する板の初期温度は，切削時の測定温度に基づき，径が 16mmの場合はを30.0℃～34.0℃の間において0.5℃刻みで変え，式(15)により計測値と解析結果の残差二乗和により各17.89℃，径が25mmの場合は26.85℃とする．

計測点における計測値と解析値の一致率について評

ミル先端近傍温度を推定した結果，エンドミルの径価する．

を 16mm とした場合は，エンドミル先端近傍温度

前節の例に従い，エンドミルの径を変えてエンドφα は 24.6℃となり（図-13），またエンドミルの径を 25mm とした場合は，エンドミル先端近傍温度

ETotal  E point A  E point B                  Nobsi 1

Aicomp.

Ai 2obs.

Nobs
i 1Bicomp.

Bi 2obs.

(15)

φα は 33.7℃となった（図-14）．エンドミル径16mm と 25mm の場合における各計測点 A,B での温度履歴を図-15,16 に整理する．実線が仮想領域法に−70−仮想領域法を用いた有限要素解析によるフライス加工時における工具切れ刃先端近傍の間接的温度推定法に関する提案(a-1) T=10sec.
図-10 エンドミル先端温度 φα と残差二乗和 E の関係(D=32mm)(a-2) T=10sec.

(a)

A点における温度履歴(b-1) T=20sec.

(b)

B点における温度履歴

図-11 φα=38.2℃の場合における個々の計測点での温度変化(D=32mm)(b-2) T=20sec.

−71−

倉橋

貴彦 ・ 羽根川

英 ・ 井伊優生 ・ 杉澤元次郎 ・ 井山徹郎 ・ 近藤俊美

基づく有限要素解析の結果であり，誤差バーにより示している結果は温度計測結果を示している．どちらの径の場合も，計測値と若干の差異はあるものの，各計測点における計測温度の平均値の推移と傾向は似ているものと考えられる．
これらの結果をもとに，エンドミルの径と本論文において提案する計算の流れにおいて推定されたエンドミル先端近傍温度の関係を整理すると図-17のようになる．結果の傾向より両変数の間には線形の関係があると考えられたため，線形の方程式により(c-1) T=30sec.

最小二乗近似を行ったところ，エンドミル先端近傍
温 度 φα は エ ン ド ミ ル の 径 を D と し た 場 合「φα=0.8575×D+11.3005」という関係で表せることが わ か っ た ． 本 論 文 に お け る 検 証 で は D=16mm,25mm,32mmのケースで検討をしているため，ここに示したエンドミル先端近傍温度φαとエンドミルの径Dの関係式は，エンドミルの径Dが16mm<D<32mmの間においての関係式であることに注意する必要がある．

(c-2) T=30sec.

図-13 エンドミル先端温度 φα と残差二乗和 E の関係(d-1) T=40sec.

(D=16mm)(d-2) T=40sec.
図-12 φα=38.2℃におけるエンドミルの位置と温度分布(D=32mm)

図-14 エンドミル先端温度 φα と残差二乗和 E の関係(D=25mm)−72−

仮想領域法を用いた有限要素解析によるフライス加工時における工具切れ刃先端近傍の間接的温度推定法に関する提案

(b)
(a)

A点における温度履歴

図-16B点における温度履歴

φα=33.7℃の場合における個々の計測点での温度変化(D=25mm)(b)

B点における温度履歴

図-15 φα=24.6℃の場合における個々の計測点での温度

図-17

エンドミル径とエンドミル先端温度の関係

変化(D=16mm)

づく非定常熱伝導解析プログラムの構築を行った．
背景領域Ωの離散化には三角形一次要素によるガラーキン法，副領域 ω の離散化にはディラックのデルタ関数に基づく方法を適用した．フライス加工時における材料表面温度の計測を行い，シミュレーション結果と計測温度の残差二乗和により評価した結果，エンドミル径が 16mm,25mm,32mmの場合に対して，本論文で取り扱った加工条件下
では，工具先端近傍の温度は 24.6℃，33.7℃，38.2℃と推定することができた．また，エンドミ
ル先端近傍温度 φα とエンドミルの径 D に対する関係を最小二乗法によるフィッティング計算をし(a)たところ，エンドミルの径 D によるエンドミル先A点における温度履歴端近傍温度 φα の推定式を誘導することができた．
加工対象によっては，被削材の温度をある一定温度以上にすることができない場合もある 11)ため，４．おわりに本論文において紹介した方法の様に，事前に加工
本研究では，フライス加工時における材料表面温度時における工具先端近傍の温度を推定しておくこを予測するため，仮想領域法および有限要素法に基とで，熱伝導解析において工作物内部や表面の温−73−

倉橋貴彦 ・ 羽根川

英 ・ 井伊優生 ・ 杉澤元次郎 ・ 井山徹郎 ・ 近藤俊美

テムの構築”，精密工学会誌，Vol.70，No.1, pp.123-

度分布を熱伝導シミュレーションにより予測する

127, 2004．

際等にも利用できるものと考えている．また，将5)

来的な課題としては，エンドミル先端近傍におけ奥田孝一，長谷川素由，布引雅之，”超精密切削加工における工具温度上昇と加工誤差”，日本機械学会論

る温度を実測する方法を考案することや，また実文集(C 編)，Vol.67，No.656, pp.1170-1175, 2001.

測結果を使うことで，数値計算により推定された

6)

エンドミル先端近傍温度と計算条件（メッシュサ

R. Glowinski, T.W. Pan, T.I. Hesla, D.D. Joseph, and J.
Periaux, “A distributed Lagrange multiplier/fictitious

イズや時間増分量）との関係についても考察を行

domain method for flows around moving rigid bodies:

う必要があると考えられる．

Application to particulate flow”, Int. J. Numer. Methods inFluids, Vol.30, No. 8, pp.1043-1066, 1999.
7)

謝辞：本研究を行うにあたりユニオンツール株式

E. Becache, J. Rodriguez, and C. Tsogka, “A Fictitious
Domain Method with Mixed Finite Elements for

会社「H25 年度 科学技術研究費助成金制度」の援

Elastodynamics”, SIAM J. SCI. COMPUT., No.29，Vol.3，助を受けた．ここに謝意を表す．本論文の解析結

pp.1244-1267, 2007.
8)

果は，九州大学情報基盤研究開発センターの高性

R. Glowinski and Q. He, “A Least-Squares/Fictitious

能アプリケーションサーバ SR16000 を使用し計算

Domain Method for Linear Elliptic Problems with Robinを行ったものである．計算機の使用についても，

Boundary Conditions”, Commun. Comput. Phys., Vol.9，No.3, pp.587-606, 2011.

センター関係者の方へ謝意を表す．
9)

J.Periaux, “A Fictitious Domain Approach to the Direct

参考文献
1)

Numerical Simulation of Incompressible Viscous Flow past

加藤隆雄，”蒸着薄膜による工具温度分布の測定”，

Moving Rigid Bodies: Application to Particulate Flow”

日本機械学会論文集(C 編)，Vol.56，No.522, pp.512-

Journal of Computational Physics, Vol.169 (2), pp.363-426,518, 1990.
2)

西一文，笹原弘之，”高速エンドミル加工時の工具温度の解析”，日本機械学会関東支部第 7 期総会講演会
講演論文集, pp.61-62, 2001.

3)

周智鵬，細川晃，山田啓司，上田隆司，”高速小径ボ

2001.
10) 森 北 肇 ， ”伝 熱 工 学 ”， 第 2 版 ， 森 北 出 版 株 式 会社,2005．
11) N.Sugita, T.Osa, M.Mitsuishi, “Analysis and estimation ofールエンドミル加工における工具温度分布の測定”，
第 4 回生産加工・工作機械部門講演会講演論文集，225, pp.65-66, 2002.
4)

R. Glowinski, T.W.Pan, T.I.Hesla, D.D.Joseph andcutting temperature distribution during end milling inrelation to orthopedic surgery”, Medical Engineering andPhysics, Vol.31 (1), pp.101-107, 2009.

新谷一博，高井一輝，角谷均，”バインダレス cBN 工
具を用いた高速切削環境における刃先温度管理シス−74−(2014. 3. 31 受付)