(1)接触問題の定式化
75
0
0
全文
(2) 接触問題の定式化. z. z. y. y. x. x (a) Point description. (b) Mesh description. Fig.4.1. Various tool descriptions. る必要があり，また被加工材との連成解析を行わなければいけないため，計算準備の手間や，計算時間 の観点から現状ではまだまだ一般的であるとは言えない．そこで本研究ではより実用的な観点から解析 を行うことを目的として，工具の弾性変形を無視して剛体として扱った解析を行うこととする． 通常の塑性加工 FEM 解析では，CAD 等で用いられるパラメータ表現をそのままのかたちで工具モデ ルとして取り扱うケースは少ない. 3)-5). ．それはこの方法は実際の工具形状に最も近いモデル化が可能で. ある一方で，取り扱いが困難であるためである．したがって通常は工具モデルに関しても何らかの離散 化を行う必要がある．しかし工具面の離散化手法は，後述するように接触探索アルゴリズムや工具面法 線ベクトルの定義と深い関わりがあるため，慎重に決定する必要がある．本研究でベースとしている ITAS3D では，次に示す二つの離散化表現方法を用いることができる．以下にそれぞれの手法の概略， および特徴を示す． (a)点集合による表現 3),6) ある平面(x-y 平面とする)上に切られた格子を考え，各格子点上に工具節点があると考えて，そのとき の z 座標の大きさで工具面を表現する方法である．x,y 方向に規則的に並ぶ点の集合で工具面を表現する ため，このような名前でよばれる．Fig.4.1(a)にその例を示す．この手法では x-y 平面上の格子をのぞけ ば工具が各点の z 座標のみで表現されるため，データが簡潔になる．また工具要素(および節点)が x,y 方 向に規則的であるため，接触探索を非常に簡単かつ高速に行うことができる．しかし x-y 平面上の格子 に基づいて節点を作成するため，局所的に格子間隔を細かくしてもその格子に沿った部分全体が細かく なる．そのため局所的な細分化が難しく，他の手法に比べてデータ量が大きくなる欠点がある．また各 工具について一つの格子点上で一つの節点しか持つことができないという制約から，縦壁やアンダーカ ット形状を有する工具には対応できないという問題がある．Fig.4.1(a)の場合，縦壁部にわずかな勾配を 設けることでいくつかの節点を持たせる工夫をしている． (b)有限要素メッシュによる表現 3),6),7) 工具表面についても材料と同様に有限要素メッシュを用いて離散化する方法であり，現在の FEM シ ―. 4-2. ―.
(3) 接触問題の定式化 ミュレーションソフトウェアで最もよく用いられている手法である．Fig.4.1(b)にその例を示す．通常有 限要素メッシュには，三角形一次要素や四角形一次要素が用いられる．四角形要素は要素内で曲面を持 つことができるため自由度が高く，かつ三角形要素の場合と比べて要素数を削減できるという利点があ る．一方で三角形要素は要素内で曲面を持つことができないという制限から，部位によっては四角形要 素に比べて要素分割を細かくする必要がある．しかし後述するように，四角形要素に比べて接触探索が 非常に簡便になるという利点を持つ． この表現方法は CAD システム上やメッシュジェネレータでモデルを作成することができ，データ作 成等の取り扱いが容易である．また点集合表現とは異なり適用できる形状に制限がなくどんな形状の工 具でも対応でき，また形状に応じて粗密をつけることもできるという利点を持つ．しかしモデル化の精 度は要素分割の仕方に大きく依存し，特にフィレット部などは非常に細かい要素分割が必要となる．ま た点集合による表現のように要素は規則的に配列されていないため，接触探索には多くの工夫を凝らす 必要がある． 本研究で対象とする THF では，第 1 章で示した例のように円形断面から最終的には四角形断面まで変 形させることが多い．そのため縦壁が多く発生する可能性があると考えられる．したがって工具面のモ デリング精度から(b)に示した有限要素メッシュによる表現方法が適していると考えられる．そこで本研 究では有限要素メッシュによる工具面形状モデリングを採用することとする．また，4.7 節で詳述する ように接触探索アルゴリズムの簡便性から，三角形一次要素による離散表現を用いることとする． 本研究では対象としないが，工具の変形を考慮した連成解析を行う場合は，被加工材と接触する面だ けではなく工具全体をモデル化する必要があるため，有限要素メッシュによる表現に限られる．この場 合用いる要素は三角形一次要素等のシェル要素(パッチ)ではなく，六面体要素などのソリッド要素とな る．. 4.2.2. 工具面法線ベクトル. 工具に接触した材料節点は，固着摩擦状態ではない限り工具面上を滑ることはできるが，工具面の垂 直方向(法線方向)へ動くことはできない．したがって工具と接触した材料節点に対しては，工具面の法 線方向へは動くことはできないという変位拘束条件を仮想仕事の原理に導入する(付帯条件として与え る)必要がある．このような変位拘束条件を導入するには，変位を拘束する方向，すなわち接触位置にお ける工具面の法線方向を規定する工具面法線ベクトルが必要となる．そこで，以下に有限要素メッシュ によって離散化された工具表面に関する法線ベクトル算出法を示す． まず基礎となる法線ベクトルの算出を行う．工具面が Fig.4.2(a)に示すように三角形一次要素で離散化 されている場合，各要素は常に平面を構成するので各要素に関する単位法線ベクトルは次式より一義的 に与えられる．. V=. a1 × a2 a1 × a2. (4.1). 一方四角形一次要素の場合は通常平面を構成しない．そこで Fig.4.2(b)に示すように各頂点(工具節点)上 で式(4.1)により単位法線ベクトルを定義し，要素内の任意の位置では次式から内挿を用いて求める．. V = ∑ N αV α. (Number of vertices). (4.2). α =1. ―. 4-3. ―.
(4) 接触問題の定式化. V3 V4 3. 3. V a2 1. 4. V1. V2. 2. a1. 1. (a) Triangular linear element. 2 (b) Quadrilateral bilinear element. Fig.4.2 Definition of tool normal. n V2n V1. Tool normal distribution. Vn. V3n 2 n. 2. θ2. θ1. 1. θ3. 3. (a) Discontinuous normal distribution. 1. 3. (b) Continuous normal distribution. Fig.4.3 Discontinuous and continuous normal distributions. ―. 4-4. ―.
(5) 接触問題の定式化 このとき， V α は工具要素頂点(工具節点)上の単位法線ベクトルであり， N α はこの場合四角形一次要素 の形状関数である． しかし，このように単位法線ベクトルを定義すると，Fig.4.3(a)に示すように一つの節点にその節点が 属する要素の数だけ単位法線ベクトルが定義されることとなる．その結果節点近傍で単位法線ベクトル の分布が不連続となる．この不連続性は変位拘束方向の急激な変化，つまり工具面の急激な勾配変化を もたらし，接触している材料節点の滑る方向を急激に変える．それにより，さまざまな不具合を引き起 こす場合がある．THF においても，第 1 章で示した部品のようにさまざまな方向に大きな曲率を有する 形状が対象となる場合が多いため，その対策は重要であると考える．そこで単位法線ベクトルに連続性 を持たせるため，本研究では Fig.4.3(b)に示すように各節点で定義される単位法線ベクトルの平均をその 節点における唯一の単位法線ベクトルと定義する．例えば図中の節点 n の場合，単位法線ベクトルV n は 次式で定義される． 3. Vn =. ∑V k =1 3. n k. (4.3). ∑V k =1. n k. ここで Vkn は要素 k 内で式(4.1)により定義される節点 n における単位法線ベクトルである．式(4.3)による 定義を用いた場合，実際の工具要素の向きと工具面法線ベクトルの向きは若干異なることになるが，工 具法線ベクトルは工具全体に連続性を持つこととなり，滑らかな工具表面を表すことができる．また， 要素内の任意の位置における法線ベクトルの算出には，式(4.2)を用いればよい．しかし式(4.3)による定 義式では，要素形状によっては定義される法線ベクトルが特定の要素の影響を大きく受ける場合が生じ る．そこで本研究では要素の面積を重みとしてかけることにより，この問題を回避している．すなわち， 3. Vn =. ∑S V k =1 3. k. n k. (4.4). ∑S V k =1. k. n k. で与えられる．このとき Sk は要素 k の面積を表す．このほかにも，Fig.4.3(b)に示すように対象としてい る節点での角度 θ k を重みとして用いる方法 7)も提案されている．. 4.3. 弱形式化 8),9). 本節では，2.5.2 節で示した updated Lagrangian 形式の速度形仮想仕事の原理式に基づいて，接触に関 する仮想仕事を導入する．本研究では塑性加工シミュレーションに特化したプログラムを開発するとい う目的から，次のような仮定を設定する． ・ 工具の変形は無視して剛体とし，材料の変形のみを考慮した定式を行う．したがって材料のみの unilateral contact problem を対象とする． ・ 接触境界を材料節点と工具要素との接触と考える． ・ 接触による拘束条件の導入には取り扱いの容易な Penalty 法を用いる． ・ 液圧以外の表面力は全て接触力に起因するとし，対称面境界条件等による付加的な節点力は定式上 無視する． 本章では，表記上の簡便さから増分形に変換する前の仮想仕事の原理式(2.35)に基づいて定式化を行う． ―. 4-5. ―.
(6) 接触問題の定式化. ê3. ê3. ê2. ê2 ê1 ê1 Material. Punch. Fig.4.4 Contact interface showing local unit vectors. ただし，時間微分記号 ⋅ に代わり増分記号 ∆ を用いることで，容易に本研究で用いる増分形に変換でき る．式(2.35)を再び示すと， T. τo − D ⋅ σ − σ ⋅ D + σ ⋅ LT : δ LdV = t& ⋅ δ vdS + p& ⋅ δ vdS ∫v ( J ) ∫sc ∫sγ . (4.5). 右辺第一項が接触による仮想仕事を表している．以下，本章では簡単のため右辺第一項のみを考える． ただしこの項は，離散化後は接触節点に対してのみ存在することに注意する． 接触節点には，工具内に潜り込まないようにその位置における工具面法線方向への変位拘束条件が与 えられる．本研究では Fig.4.4 のようにそれぞれの接触位置において工具面法線ベクトルを eˆ3 とする直交 デカルト系で表された局所座標系 eˆi を設定する．接触位置における工具面法線ベクトル ê3 は式(4.2)から 求められる．このとき接触力速度ベクトル t& は. (. ) ≡ ( tˆ eˆ. t& = tˆ3eˆ3 + tˆ i eˆi. N 3. + tˆTi eˆi. ). ( i = 1, 2 ). (4.6). と表すことができる．このとき tˆN eˆ3 は法線方向の接触力であり， tˆTi eˆi は摩擦力である．これを式(4.5)に 代入すると，. ∫. Sc. t& ⋅ δ vdS = ∫. Sc. =∫. Sc. ( tˆ eˆ. N 3. (. ). + tˆTi eˆi ⋅ δ vdS. ). tˆ&N eˆ3 + tˆN eˆ&3 + tˆ&Ti eˆi + tˆTi eˆ&i ⋅ δ vdS. (4.7). となる．最右辺第一項は接触境界に関する拘束条件によって生じる， eˆ3 方向の接触力速度に関する項で ある．右辺第二項は工具面法線ベクトル eˆ3 の変化に起因する項である．この項は工具が曲率を有してお り，工具面法線ベクトルが勾配を持つ場合のみ表れる．また第三項，四項は摩擦に起因して発生する項 である．摩擦を無視した解析では，この二つの項は無視できる． 冒頭で示した仮定に基づいて，接触による拘束条件の導入に Penalty 法を用いると， tˆ&N は. tˆ&N = α vˆN. (4.8) ―. 4-6. ―.
(7) 接触問題の定式化 という速度との関係で与えられる．このときαは penalty 数であり，非常に大きな正の値である．本研 究では工具を剛体として扱う，一方向接触問題(unilateral contact problem)を仮定しているので，式(4.8)よ り明らかなように penalty 数に起因する大きな負の値が剛性マトリクスに入ることはない．したがって. αは十分に大きく取ることができるため，本研究ではα=1040 としている．αとしてこの程度の大きさ を用いた場合，Penalty 法によって許容される工具への潜り込みは数値計算上生じない．これを式(4.7)に 代入すると，. ∫. Sc. ( (δ vˆ α vˆ. ). t& ⋅ δ vdS = ∫ α vˆN eˆ3 + tˆN eˆ&3 + tˆ&Ti eˆi + tˆTi eˆ&i ⋅ δ vdS Sc. =∫. Sc. N. N. ). (4.9). + tˆN δ v ⋅ eˆ&3 + δ vˆi tˆ&Ti + tˆTiδ v ⋅ eˆ&i dS. 式(4.9)が本研究における接触に関する仮想仕事である．通常塑性加工解析において式(4.9)を取り扱う場 合，前述のように接触境界を材料節点と工具要素との接触に置き換えて，全ての接触節点についての仮 想仕事の総和をとることにより接触力による仮想仕事を近似的に求める．したがって接触節点 m におけ る集中接触力を m t c とすると，式(4.9)は次式のように近似される．. ∫. Sc. (. t& ⋅ δ vdS ≈ ∑ m t c ⋅ δ m v = ∑ δ m vˆNα m vˆNc + mtˆNc δ m v ⋅ m eˆ&3 + δ m vˆi mtˆ&Ti + mtˆTiδ m v ⋅ m eˆ&i m. m. ). (4.10). このとき左肩添え字 m は，これらの変数が節点 m 上で離散化された量であることを示す．以下，式(4.10) に基づいて詳細な取り扱い方法を示す．右辺第二項の工具曲率に起因する仮想仕事の取り扱いは 4.4 節 に，第三，四項の摩擦に起因する仮想仕事の取り扱いには 4.5 節で示す．また右辺第一項の Penalty 法に よる変位拘束条件の剛性マトリクスへの組み込み方法は 4.6 節にて示す．. 4.4 工具曲率に伴う接触力補正アルゴリズム 塑性加工において，接触した材料節点が工具面上を滑る現象は非常の頻繁に発生する．これはハット 曲げや角筒/円筒絞りなどの深絞り成形を考えれば明らかである．また THF においても，材料は軸押し 込みを受けながら拡管するため工具面上の軸方向および周方向への滑りは頻繁に生じる．例えば工具上 の平坦な部分を滑る場合，接触位置の変化によって工具面法線ベクトルは変化せず，その結果変位が拘 束される方向も変化しない．しかし曲率を有した部分を滑る場合には，接触位置によって工具面法線ベ クトルも変化し，それにより増分の前後で変位が拘束される方向が変化することとなる．したがって接 触節点が工具面上を滑る場合， 4.7 節で示す接触探索アルゴリズムによって接触位置および工具面法線 ベクトルを時々刻々更新する必要がある． 増分ステップ中に工具面法線ベクトルが変化した場合，それに伴って接触力の方向も補正する必要が ある．それは接触力から摩擦力を差し引いた残りの分が，常にその瞬間の位置での ê3 方向を向く必要が あるためである．逆に摩擦を無視した場合は，接触力は常に eˆ3 の方向と一致している必要がある．つま り言い換えれば，接触力は工具面法線ベクトルの変化に伴う共回転ベクトル，いわゆる追従力(follower force)10)として取り扱う必要があるといえる．この接触力の取り扱いが適切でないと，THF どころかハ ット曲げや角筒絞りといった非常に単純な成形過程でさえまともに解析することができない．つまり接 触力の取り扱いは塑性加工解析においては非常に重要であるといえる． ITAS3D ではこの問題に対して，工具面法線ベクトルの変化に伴って生じる接触力の接線方向分(摩擦 がある場合は摩擦力ベクトルとの差)を不釣り合い力として考え，3.6 節で示した一ステップ遅れ補正法 ―. 4-7. ―.
(8) 接触問題の定式化 により補正する方法を用いている．しかし 3.6 節で示したようにこの手法はあくまでも近似的な補正法 であるため，工具面法線ベクトルの変化に起因する不釣り合い力を精度よく補正できる保証はない．ま た，この不釣り合い力を不釣り合い力補正アルゴリズム ALGONEQ のなかで補正する手法が用いられる 場合もある．しかしこの手法に関しては十分な検討が行われていないため，その解析精度に関しては明 らかとなっていない．すなわち現状では，工具面法線ベクトル変化に伴う接触力の変化を不釣り合い力 としてとらえる方法が適切かどうか明らかであるとは言えない． そこで本研究では三角形要素で記述された剛体工具を用いるという仮定の枠内で，不釣り合い力補正 としてではなく増分ステップ中に精度よく補正を行うことのできる新たなアルゴリズムを提案した．本 節ではまず参考のために二次元問題における接触力補正アルゴリズムを示し，続いて従来から用いられ ている不釣り合い力ととらえた補正方法，そして本研究で提案するアルゴリズムを示す．. 4.4.1. 二次元問題における接触力補正アルゴリズム 11),12). Fig.4.5 に工具面上を材料節点が滑る様子を示す．ここでは簡単のため摩擦なしとする．∆urel は材料節 点と工具の相対変位増分ベクトルである．増分前の配置(時刻 t)における工具面法線ベクトル n および接 触力ベクトル t が，増分後(時刻 t ′ = t + ∆t )に n′ ， t ′ に変化するとする．また増分前後それぞれの配置に おいて，n を n 方向，および n′ を n′ 方向とする直交デカルト系表示された局所座標系 l − n ，l ′ − n′ を設定 するとする．このとき接触力増分ベクトル ∆t は次式で与えられる．. t ′ − t=∆t=∆ ( tn n ) = ∆tn n+tn ∆n. (4.11). ここで tn は接触力ベクトルの n 方向成分であり， ∆n は. ∆n=n′ − n. (4.12). t′ t = t + dt n′. l′. tool surface. t. ∆urel. n. ∆θ. t =t. l. Fig.4.5 Sheet node sliding on curved tool surface for 2-D case. ―. 4-8. ―.
(9) 接触問題の定式化 で与えられる．式(4.11)の右辺第二項が工具面法線ベクトル変化に伴う接触力補正項であり，式(4.9)右辺. (. ). 第二項に相当する． n から n′ への回転を表す直交テンソルを R =ei′ ⊗ ei = ( e′j ⋅ ei ) ei ⊗ e j とおくと，式 (4.12)は. ∆n=n′ − n=R ⋅ n − n= ( R − I ) ⋅ n. (4.13). と変形できる．ただし I は単位テンソルを表す．ここで Fig.4.5 に示すように工具面法線ベクトルの回転 角を ∆θ とすると， R は局所座標 l − n 系を参照して. l ′ ⋅ l n′ ⋅ l cos ( ∆θ ) − sin ( ∆θ ) Rij = l ′ ⋅ n n′ ⋅ n = sin ( ∆θ ) cos ( ∆θ ) . (4.14). と与えられる．したがって， l − n 系を参照した ∆n の成分は式(4.13)，(4.14)より. ∆nl cos ( ∆θ ) − 1 − sin ( ∆θ ) nl = cos ( ∆θ ) − 1 nn ∆nn sin ( ∆θ ) となる．ここで一増分ステップでの相対変位増分が十分小さく， ∆θ. sin ( ∆θ ). (4.15). 1 が成り立つとすると，. ∆θ (4.16). cos ( ∆θ ) 1 と近似することができる．また l − n 系では ( nl , nn ) = ( 0, 1) となるので，式(4.15)は. ∆nl 0 = ∆nn ∆θ. −∆θ 0 −∆θ = 0 1 0 . (4.17). となる．一方相対変位増分ベクトル ∆urel は. ∆u rel = ∆ulrel l. (4.18). と表される．工具面が円弧で表すことができると仮定すると，回転角と円弧の関係より式(4.18)を用い て. ∆θ =. ∆ulrel. (4.19). ρ. となる．式(4.17)，(4.19)を，式(4.11)を l − n 系で成分表示した式に代入すると，. ∆ulrel ∆tl 0 ∆nl 0 − = + tn = + tn ρ ∆tn ∆tn ∆nn ∆tn 0 . (4.20). ここで工具変位増分ベクトルを ∆utool と置くと，最終的に次のような式が得られる．. tn ∆tl 0 − ∆ul = + ρ ∆tn ∆tn 0 . tn tool ∆ul + ρ 0 . (4.21). 右辺第二項は材料節点の変位増分の関数であるため，剛性方程式において左辺へ移項し，剛性マトリク ―. 4-9. ―.
(10) 接触問題の定式化 スへ組み込む．すなわち，この節点に関する［2×2］のサブマトリクスに. tn ρ 0. 0 0 . (4.22). を足せばよい．この接触力補正項は対称なマトリクスとなる． 式(4.21)より明らかなように接触力補正項は変位成分 ∆ul のみの関数であり局所座標系では一次元的 な変化しかしない．したがって詳細は省略するが，摩擦を考慮した統一的な接触力補正の定式も容易に 行うことができる．このようにこのアルゴリズムは，接触力補正項を変位増分の関数として表すことで 工具面法線ベクトル変化を考慮した剛性マトリクスを得ることができ，その結果増分後には工具面法線 ベクトルに整合した接触力ベクトルが得られる．. 4.4.2. 三次元問題における接触力補正アルゴリズム. 前節で示した二次元問題における接触力補正項の定式では，工具曲面が円弧近似できると仮定するこ とで，摩擦も含めて非常に簡便な定式化を行うことができた．一方三次元問題では，4.2 節で示したよ うにたいていの場合は工具表面を離散化しているため，二次元問題で用いたような仮定を導入すること はできない．また工具面法線ベクトルで規定される平面内の任意の方向に補正する必要があるので，扱 いが複雑となる．そこで本研究では問題を明確化するため，工具曲率による接触力補正と摩擦は分離し て考え，それぞれ全く異なるアルゴリズムを適用することとする．したがって本節では摩擦を無視して， 三次元問題における工具曲率に関する接触力補正アルゴリズムのみを示す． 本節では，まず ITAS3D において従来から用いられている手法を示し，その問題点を明らかにした上 で，本研究で新たに提案したアルゴリズムを示す．. (1). 不釣り合い力としてとらえた手法 13),14). 本アルゴリズムは，3.6 節で示した不釣り合い力補正アルゴリズムを接触問題まで拡張した手法であ る．つまり，工具面法線ベクトルの変化によって生じる接触力の接線方向成分を不釣り合い力としてと らえ，これを次ステップ以降，もしくは ALGONEQ において補正する手法である．本アルゴリズムの概 念を Fig.4.6 に示す．図中の記号は Fig.4.5 と同様であるが，工具面接平面内の成分が一つ増えるため， ここでは接平面を表す基底ベクトルを l1 , l2 と表現している．またこれら基底ベクトルは，Fig.4.4 におけ る記号とは l1 ↔ eˆ1 , l2 ↔ eˆ2 , n ↔ ê3 のようにそれぞれ対応する．簡単のため，増分前の時刻 t では接触力 ベクトルは厳密に接触位置における工具面法線ベクトル n と同じ方向を向いているとする．増分中に工 具曲率に伴う接触力補正をなにも行わなければ，すなわち式(4.10)の第二項を無視すれば，次式で与え られる A tl ′ が残存する． A. tl′ = A t ′ − ( n′ ⊗ n′ ) ⋅ A t ′= A t ′ − A t n′. (4.23). ここで添え字 A は節点 A に関する接触力であることを示している．本アルゴリズムではこれを. t. A neq cont. = − A tl′. (4.24). neq のように不釣り合い力と考える．一ステップ遅れ補正法では，この不釣り合い力 A t cont に加速度パラメー. タ 1 Rsafe をかけて，次式のように節点剛性方程式の右辺に組み込む． ― 4-10. ―.
(11) 接触問題の定式化. A. t =t. t. A. tn′. A. t ′= At+ Adtn n. n. l2 l1. A. durel. t = t + dt. n′. t. A neq cont. l2′ l1′. A. tl′. Fig.4.6 Definition of non-equilibrated force due to curvature of tool in three dimensional case. A. K ⋅ ∆u= ∆f+ A. A. t Rsafe. A neq cont. (4.25). これを全ての節点について重ね合わせることにより， neq t cont K ⋅ ∆u=∆f+ Rsafe. (4.26). という不釣り合い力補正を考慮した全体剛性方程式を得ることができる．THF 解析の場合は式(3.105)に neq 組み込めばよい．また ALGONEQ15)を接触問題まで拡張して用いる場合には，式(3.120)の右辺に t cont を. 加えて， neq K ⋅ ∆u = f neq + t cont. (4.27). を増分ステップ後に解くこととなる．以下本論文では，式(4.27)に基づく ALGONEQ を. 拡張した. ALGONEQ(extended ALGONEQ) とよぶこととする． これらの手法は，従来から用いられるアルゴリズムを応用することで簡便に接触力補正を行うことが できるという利点を有する．しかしその一方で，これらのアルゴリズムでは増分ごとのきちんとした釣 り合いを考えておらず，定式上，式(4.9)の右辺第二項を無視することと等価である．そのため結果とし neq て t cont には「前ステップで発生した不釣り合い力」だけではなく， 「直前のステップまでに蓄積された不. 釣り合い力」が含まれることとなり，その力学的な意味があいまいとなる．特に加速度パラメータ Rsafe ―. 4-11. ―.
(12) 接触問題の定式化. o A. t′. A. t. A neq. t n′′. A. A. t ′= At+ Adtn n. tl′. Fig.4.7 Relationship between co-rotational contact force vector and the one obtained by one step behind correction method. の設定が補正度合いに大きな影響をおよぼす一ステップ遅れ補正法では，この問題は顕著となる．また これらの手法だけでは増分後に新たな不釣り合い力が発生することを許容する必要がある．さらに，本 来は工具面法線ベクトルの回転に伴う共回転接触力ベクトルとして扱うべきだが，本手法では発生した 接線方向成分を消去するだけである．その結果 Fig.4.7 に示すように本来発生すべき共回転接触力ベクト o. ル A t ′ よりも実際得られる接触力ベクトル A t n′ ′ が小さくなるという問題もある． 以上のような問題点から，接触力変化を不釣り合い力ととらえた補正法だけでは解析精度上問題が残る と考えられる． (2). 整合した接線剛性マトリクス(consistent tangent stiffness matrix). 通常の静的陰解法 FEM では，静的陽解法とは異なり一増分ステップ内であらゆる不釣り合い力がほ ぼ消去されるまで反復計算を行う．それにより，変形後の配置における力の釣り合いを満足することが できる．この際，反復計算の計算効率向上を図ることを目的として，反復計算中の接触位置や工具面法 線ベクトルの変化を考慮した厳密な接線剛性を用いる必要性があることが指摘されている 16),17)．この厳 密な接線剛性とは，反復計算中の接触位置の変化や工具面法線ベクトルの変化などを変位増分の関数と して表して剛性マトリクスに組み込んだ，接触状況の変化に整合した接線剛性マトリクス(consistent tangent stiffness matrix)16)のことを示す．4.4.1 節で示した二次元問題における接触力補正アルゴリズムも この consistent tangent stiffness matrix の枠組みにおける定式であるといえる．そこで本研究ではこの手法 に基づいて，三角形要素で記述された剛体工具との unilateral contact problem に特化した consistent tangent stiffness matrix を新たに提案した． 接触による離散化された仮想仕事を再録すると，次式のようになる．. ∫. Sc. (. t& ⋅ δ vdS ≈ ∑ δ m vˆNα m vˆNc + mtˆNc δ m v ⋅ m eˆ&3 + δ m vˆi mtˆ&Ti + mtˆTiδ m v ⋅ m eˆ&i m. ). (4.10). 工具の曲率による接触力補正項は右辺第二項である．そこで接触節点 m に注目してこの項の具体的な離 散化を行う．まず，工具面法線ベクトルの時間変化 m eˆ& の離散化を考える．以下，左肩の添え字 m はい 3. ずれも節点 m 上で与えられる(定義される)諸量であることを示す．二次元問題の場合と同様に，増分前 ― 4-12. ―.
(13) 接触問題の定式化 の配置(時刻 t)における工具面法線ベクトル m n および集中接触力ベクトル m t c が，増分後(時刻 t ′ = t + ∆t ) に m n′ ， m t c′ に変化するとする．増分前の配置において， m n を m ê とする局所座標系 m ê を設定する．4.2 3. i. 節で示したように，工具要素上の任意の位置における工具面法線ベクトルは各要素節点からの内挿によ って求められる．したがって，三角形一次要素の場合工具面法線ベクトルは要素内で線形に分布した位 置の関数として表すことができる．以上の考察より，工具面法線ベクトルの増分 ∆ m n は微分の連鎖則に より次式のように表すことができる．. ∆ m eˆ3 = ∆ m n =. ∂ m n m ∂ m n m rel ∆ x= ∆ u ∂x ∂x. (4.28). ITAS3D では接触節点に関する自由度は全て局所座標系 eˆi を参照して表示している．そのため以下の定 式中で用いるベクトルは全て m êi 系で表示する必要がある．そこで式(4.28)を. m. êi 系を参照して成分表示. すると，. ∆ m nˆi =. ∂ m nˆi m ∂ m nˆi m rel ∆ xˆ j = ∆ uˆ j ∂xˆ j ∂xˆ j. (4.29). これより，∂ m nˆi ∂xˆ j は変位増分と工具面法線ベクトル増分の関係を表す曲率テンソルにような性質を有 しているといえる．工具面法線ベクトル m n は形状関数を用いた内挿により与えられるので，Fig.4.3 に 示す記号に従うと，. ∆ m nˆi =. 3 ∂ m nˆi m rel ∂N α ∆ uˆ j = ∑ ˆj ∂xˆ j α =1 ∂x. M. Vˆiα ∆ muˆ rel j. (4.30). このとき M Vˆiα は節点 m が接触している工具要素 M の節点α上で定義される工具面法線ベクトル成分で ある．式(4.30)をマトリクス表示すると. ∂N α ˆ ∂x m ∆ nˆ1 1α m ∂N ∆ nˆ2 = ∆ m nˆ ∂xˆ1 3 ∂N α ∂xˆ1. M. Vˆ1α. ∂N α ∂xˆ2. M. Vˆ2α. ∂N α ∂xˆ2. M. Vˆ3α. ∂N α ∂xˆ2. M. Vˆ1α. ∂N α ∂xˆ3. M. Vˆ2α. ∂N α ∂xˆ3. M. Vˆ3α. ∂N α ∂xˆ3. Vˆ1α ∆ m uˆ rel 1 M ˆ α m rel V2 ∆ uˆ2 ∆ m uˆ rel 3 M ˆα V3 M. (4.31). 接触節点 m に対して ∆ m uˆ3rel = 0 という変位拘束条件が与えられるため，式(4.31)は. ∂N α M ˆ α ∂xˆ V1 m ∆ nˆ1 1α m ∂N M ˆ α V2 ∆ nˆ2 = ∆ m nˆ ∂xˆ1 3 0 . ∂N α ∂xˆ2 ∂N α ∂xˆ2. M. Vˆ1α. M. Vˆ2α. 0. 0 ∆ muˆ rel m 1rel 0 ∆ uˆ2 ∆ m uˆ rel 3 0 . (4.32). となる．式(4.32)が工具面法線ベクトルの時間変化の離散化表現であり，変位増分ベクトルの関数とな る．式(4.30)を，式(4.10)を増分形にした関係式の右辺第二項に代入すると， ― 4-13. ―.
(14) 接触問題の定式化. 3 ∂N α tˆ δ m v ⋅ m eˆ&3 = δ m vˆi mtˆNc ∑ α =1 ∂xˆ j . m c N. M. Vˆiα ∆ m û rel j . (4.33). 式(4.33)が本研究で提案した，接触節点 m に対する工具曲率に伴う仮想仕事の離散化表現である．これ を式(4.10)のように全ての接触節点について総和をとることにより，系全体に関する工具曲率に伴う仮 想仕事が得られる．ここで. m c 3 ∂N α ∑m tˆN ∑ α =1 ∂xˆ j . M. Vˆiα ∆ muˆ rel ≡ K C ⋅ ∆urel j . (4.34). とおき，これを全体剛性方程式に加えると，. (. K ⋅ ∆u=∆f + K C ⋅ ∆u − ∆utool ∴. ). (4.35). ( K − K C ) ⋅ ∆u = ∆f − K C ⋅ ∆utool. となる．したがって最終的に式(4.35)を用いることにより，増分ステップ中の工具面法線ベクトル変化 に整合した接触力増分を得ることができる．また THF 解析を行う場合は K C を式(3.106)に組み込めばよ い．このときの K C の存在により，剛性マトリクスは非対称となることに注意する． 式(4.33)で用いる三角形一次要素における形状関数の微分 ∂N α ∂xˆ j は次のように求めることができる． まず微分の連鎖則より. ∂N α ∂N α ∂ξ j = ∂xˆi ∂ξ j ∂xˆi. (sum on j, j = 1, 2 ). (4.36). ここで j = 1, 2 なのは，. N 1 = ξ1 , N 2 = ξ 2 , N 3 = 1 − ξ1 − ξ 2. (4.37). のように独立変数が二つのためである．ここで ξ1 , ξ 2 ( = η ) は三角形要素に関する形状関数(面積座標)であ る．また局所座標系で表示しているため平面内の内挿となり，理論上 x̂3 成分も存在しない．本研究では 工具面法線ベクトルは節点上で平均化しているため，局所座標系 ê1 − ê2 の構成する面と実際に工具要素 が構成する平面が異なる可能性がある．したがって実際は，工具節点座標は x̂3 成分を有している可能性 がある．しかしこの誤差は微小であり，無視できる程度である． 式(4.36)をマトリクス表示すると，. ∂N α ˆ ∂x1 α ∂N ∂xˆ2. ∂ξ1 ˆ ∂x1 = ∂ξ1 ∂xˆ2. ∂ξ 2 ∂N α ∂xˆ1 ∂ξ1 ∂ξ 2 ∂N α ∂xˆ2 ∂ξ 2. . (4.38). ここで Jacobian の計算には，. ― 4-14. ―.
(15) 接触問題の定式化. ∂ξ1 ∂xˆ 1 ∂ξ1 ∂xˆ 2. ∂ξ 2 ∂xˆ1 ∂xˆ1 ∂ξ1 = ∂ξ 2 ∂xˆ1 ∂xˆ2 ∂ξ 2. ∂xˆ 2 ∂ξ1 ∂xˆ 2 ∂ξ 2 . −1. (4.39). で与えられることを利用する．三角形一次要素内の任意の点は 3. α 1 2 3 xˆi = ∑ N α xˆi = ξ1 xˆi + ξ 2 xˆi + (1 − ξ1 − ξ 2 ) xˆi. (4.40). α =1. で表されることから，これを式(4.39)に代入すると，. ∂ξ1 ∂xˆ 1 ∂ξ1 ∂xˆ 2. ∂ξ 2 1 3 ∂xˆ1 xˆ1 − xˆ1 = ∂ξ 2 xˆ12 − xˆ13 ∂xˆ2 . 1 3 xˆ2 − xˆ 2 2 3 xˆ 2 − xˆ2 . −1. (4.41). したがって，式(4.38)は次のようになる．. ∂N α ˆ ∂x1 α ∂N ∂xˆ2. ˆ1 ˆ3 x1 − x1 = 2 3 xˆ1 − xˆ1 . ∂N α 1 3 −1 xˆ2 − xˆ 2 ∂ξ1 α 2 3 xˆ2 − xˆ2 ∂N ∂ξ 2. . (4.42). ところで，α=1,2 の場合. ∂N 1 ∂ξ1 1 1 = for α=1， ∂N 0 ∂ξ 2 . ∂N 2 ∂ξ1 0 2 = for α=2 ∂N 1 ∂ξ 2 . {. と陽的に与えられるので，これらを式(4.42)に代入することで ∂N 1 ∂xˆ1. (4.43). ∂N 1 ∂xˆ2 } ，{∂N 2 ∂xˆ1 ∂N 2 ∂xˆ2 }. を得ることができる．またα=3 の場合は次式より求めることができる．. ∂N 3 ∂ (1 − ξ1 − ξ 2 ) ∂ξ1 ∂ξ 2 ˆ ∂xˆ1 ∂x1 ∂xˆ1 ∂xˆ1 3= = − − ∂N ∂ (1 − ξ1 − ξ 2 ) ∂ξ1 ∂ξ 2 ∂xˆ2 ∂xˆ2 ∂xˆ2 ∂xˆ 2 . (4.44). 以上の計算により，マトリクス K C の成分を求めることができる． 元々consistent tangent stiffness matrix は， 厳密な接線剛性マトリクスを求めることにより陰解法 FEM に おける反復計算の収束性(特に Newton-Raphson 法における二次収束性)を向上させることを目的に提案 されたアルゴリズムである. 16). ．したがって陰解法 FEM においては計算効率上必要であっても，不可欠. なアルゴリズムではない．しかし陽解法 FEM は反復計算を行わずに微小な不釣り合いを許容しながら 進める解法であるため，発生する不釣り合い力を低減させるために厳密な接線剛性マトリクスを求める ことは非常に重要であるといえる．したがって陽解法 FEM では consistent tangent stiffness matrix に求め ― 4-15. ―.
(16) 接触問題の定式化 られる役割は大と考えられる．これまでにも工具が点集合表現されている場合を対象とした定式が Nakamachi らによって提案されている 18),19)． 本来非線形 FEM では増分後の配置における厳密な釣り合いが求められるため，たとえ consistent tangent stiffness matrix を用いたとしても，増分ステップ中の反復計算による不釣り合い消去が必要とな る．しかし本アルゴリズムでは三角形一次要素で記述された剛体工具との unilateral contact problem に特 化した定式を行っているため，増分後も材料節点が同じ工具要素上にとどまっている限り，反復計算を 用いることなく接触力は精度よく補正される． 一方で本アルゴリズムは一工具要素内での内挿を基本とした定式なので，材料節点が工具要素をまた がって滑る場合には補正しきれずに不釣り合い力(接線方向の接触力)が発生する可能性がある．これは 従来の consistent tangent stiffness matrix でも発生する問題でもある 20),21)．また増分ステップ中に予測され た増分後の工具面法線ベクトルと，増分後に接触探索によって得られる工具面法線ベクトルは必ずしも 一致する保証はない．したがってこの不整合性によって不釣り合い力が発生することも考えられる．さ らに，当然陽解法 FEM 特有の二次以上の微小項を無視することによる不釣り合いも発生する．しかし 静的陽解法 FEM では rmin 法によりステップ幅を小さく制限しているため，通常複数の工具要素を一ス テップで移動するような大きな変位増分は発生しない．また工具面法線ベクトルの勾配が大きい部分に ついては通常細かく要素分割されるため，要素をまたがって移動したとしても隣接要素間での工具面法 線ベクトルの勾配は小さい．したがって，上記の理由から発生する不釣り合い力は十分小さいと考えら れる．また実際に補正しきれずに発生する不釣り合い力については，先に示した一ステップ遅れ補正法 によって補正することにより，さらなる高精度化を図ることができると考えられる． 最後に，本論文中で示した三つの手法に関して，その特徴および解くべき方程式についてまとめた表 を Table4.1 に示す．これより最終的に解くべき方程式が異なることから，これら三つの手法を自由に組 み合わせて用いることも可能である．例えば前述のように consistent tangent stiffness matrix と一ステップ neq 遅れ補正法を併用するには，consistent tangent stiffness matrix で用いる剛性方程式の右辺に t cont Rsafe を組. Table 4.1 Comparison of the algorithms for curvature term. Algorithm. One step behind correction method. Treatment of the term. ∫. Sc. tˆ δ m v ⋅ m eˆ&3dS. m c N. Omitted from the principle of the virtual work and resulting non-equilibrated forces are cancelled in the subsequent steps. Stiffness equation to be solved neq t cont K ⋅ ∆u= ∆f+ Rsafe. In regular calculation, Omitted from the principle of the virtual work Extended ALGONEQ. and resulting non-equilibrated forces are In ALGONEQ, cancelled in ALGONEQ15). Consistent tangent stiffness matrix method. K ⋅ ∆u=∆f. neq K ⋅ ∆u = f neq + t cont. Discretized as the function of displacement increment vector ― 4-16. ―. ( K − K C ) ⋅ ∆u = ∆f − K C ⋅ ∆utool.
(17) 接触問題の定式化 み込めばよい．また拡張した ALGONEQ との併用の場合は，consistent tangent stiffness matrix で用いる剛 性方程式を解いた後に発生する不釣り合い力を，拡張した ALGONEQ で補正すればよい．. 4.5 摩擦の取り扱い 塑性加工では，材料と工具との間に働く摩擦力が変形挙動に大きな影響を与えることはよく知られて いる．また材料に発生する応力やひずみの分布が摩擦状態に大きく左右されることから，成形限界や形 状凍結性を予測するうえでも摩擦を考慮することは重要である 6)．これは THF においても当然発生する 問題であり，例えば円管の液圧バルジ型成形において潤滑状態の違いにより変形特性が大きく異なるこ とが淵澤らの研究によって明らかにされている 22)． したがって THF の FEM 解析を行う場合においても， 精度よい解析を行うには摩擦を考慮することが不可欠といえる． FEM における摩擦の取り扱いに関しては，これまでにも多くの定式が提案されている. 23)-26). ．しかし. 実際の摩擦現象は，固着摩擦*1 から滑り摩擦への遷移(またはその逆)や，固着摩擦時の微小突起の変形 による微小変位 23)，面圧変化や潤滑状態の変化に伴う摩擦係数の変化 27),28)など，非常に非線形性の高い 挙動を示す．したがってそのためあらゆる条件に対して統一的に取り扱うことのできるモデルは未だに 提案されていない． そのなかでも，摩擦力は接触する二つの物体間に働く垂直力に比例するとする Coulomb 則は，汎用コー ドも含めて FEM 解析ではよく用いられている. 29). ．これは取り扱いが簡便であるにもかかわらずさまざ. まな条件に適用可能であり，また板成形程度の面圧下では摩擦係数を一定とする仮定がある程度の妥当 性を持って受け入れられているためである 30)．Coulomb 則はその挙動からポテンシャルを用いて弾塑性 構成式と同様の定式を行うことができることも，よく用いられる一因となっている Owen による Penalty 法に基づく定式. 23),31)-34). ．Peric &. 26). はなかでも非常に有名な定式の一つである．しかし摩擦は基本的. に非線形性の高い現象であるため，ほとんどの摩擦モデルにおいて反復計算による不釣り合い補正が不 可欠である． そこで本研究では以上の点を踏まえて，Coulomb 則に基づいた次のような定式を行う．この定式で前提 とする仮定は以下の通りである． ① 固着摩擦時にも微小な変位を許容する，いわゆる疑似固着状態(pseudo-sticking)を仮定する． ② 反復計算を避けるため，滑り摩擦によって発生する摩擦力は一ステップ遅れ補正法により補正する． ①で示した仮定を模式的に図示すると，Fig.4.8 のようになる．ただし図中の縦軸は節点力ベクトルの工 具面法線方向成分と工具面接線方向成分の比，すなわち摩擦係数を示している．したがって固着摩擦か ら滑り摩擦に遷移する瞬間に動摩擦係数 µ d に達するとする．以下に本研究で採用する定式を示す． まず，ある増分ステップ後に接触と判定された材料節点は，次の増分ステップから節点力増分 ∆f と 節点相対変位増分 ∆u に対して次のような固着摩擦状態の条件が与えられる． rel. ∆fˆ1 E t ∆fˆ2 = 0 ˆ 0 ∆f 3 . 0 Et 0. 0 ∆uˆ1rel 0 ∆uˆ2rel En ∆uˆ3rel . (4.45). このとき ∧ は各ベクトルが接触位置における局所座標系 êi を参照した成分表示となっていることを示す． *1. 高面圧の塑性加工において材料が工具に固着されたまま移動する状態を「付着摩擦」あるいは「固着摩擦」とよぶことがあるが，ここ. での定義はそれとは別の意味である．. ― 4-17. ―.
(18) 接触問題の定式化. tt tn. µd. tan −1 ( Et ). u stick. slip. Fig.4.8 Pseudo stick-slip friction model. Et は疑似固着に関する接線方向の penalty 数(剛性)を表している． Et → ∞ とすると理論上の固着状態に 一致する．また En は接触境界に関する変位拘束条件を与える penalty 数であり，式(4.8)におけるαと等 価である．前述のように本研究では En = 10 40 としている． 固着摩擦状態から滑り摩擦状態への遷移に関する判定は次式によって行う． 2 2 2 2 F = fˆ1 + fˆ2 − µ d fˆ3. (4.46). If(F<0)then Keep sticking Else Go to sliding Endif 式(4.46)の F は弾塑性理論における塑性ポテンシャルとの類似性から滑りポテンシャルとよばれること もある 29)．滑り摩擦から固着摩擦への遷移も同様の手順によりに判定することができる． 滑り摩擦状態にある節点には，ある増分ステップ後には次式のような摩擦力 f. fric. が発生しているはず. である．. f. fric. ∆u rel = − µ d fˆ3 ∆u rel. (4.47). このとき ∆urel はその増分ステップで生じた，その節点の相対変位増分ベクトルであり，また局所座標系. êi は増分後の配置で定義される基底ベクトルである．しかし工具面の接線方向に実際に発生している節 点力ベクトル fˆ ê +fˆ ê は反復計算等により不釣り合い補正を行わない限り，式(4.47)の f fric と一致しな 1 1. 2 2. い．特に滑り方向に大きな変化が発生した場合はこの差は大きくなる．そこで， ― 4-18. ―.
(19) 接触問題の定式化. (. ). neq f fric = fˆ1eˆ1 +fˆ2 eˆ2 − f. fric. (4.48). neq を定義し，これを一ステップ遅れ補正法により消去 のように摩擦力に関する不釣り合い力ベクトル f fric. する 13),14)．したがって 4.4.2 節と同様に次の増分ステップにおける剛性方程式の右辺に加えて，. K ⋅ ∆u=∆f+. neq f fric. (4.49). Rsafe. となる． この手法は 4.4.2 節でも示したようにさまざまな問題を含んでいる．しかし前述のように摩擦は非常 に非線形性が高く，どのような定式を導入したとしても反復計算による不釣り合い補正は不可欠である． 一方で接触を考慮した弾塑性問題では，接触に関する反復計算の収束性が大きなボトルネックとなって おり，陰解法 FEM が用いられない大きな要因となっている．このような観点から，この定式は不釣り 合い力を完全に消去できる保証はない一方で陽的に扱うことができ，またある程度妥当性のある結果が 得られることから，陽解法 FEM においては非常に実用的な手法であると言える．また THF では，軸方 向の圧縮変形と周方向の引張り変形が主な変形形態であることから，摩擦の発生する方向は変形中大き く変化しないと考えられる．したがって式(4.48)で定義される摩擦に起因する不釣り合い力ベクトルは それほど大きくならないと考えられ，これより本研究の範疇では本手法は十分適用可能であると考えら れる． しかしその一方で，このアルゴリズムだけでは式(4.10)の右辺第三，四項を完全に無視して増分ステ ップの計算を行うことと等価であり，力学的な意味が不明瞭である．したがって工具曲率に起因する接 触力補正の場合と同様に，増分ステップ中にある程度でも摩擦に起因する不釣り合い力を補正する定式 を導入するなど，さまざまな問題が今後の課題として残っている．. 4.6 接触境界条件の剛性方程式への導入 接触による拘束条件の導入には，Lagrange 未定乗数法や Penalty 法. 8),9),29). ，拡大 Lagrange 未定乗数法. 35). ,また接触に関する Kuhn-Tucker 条件を関数化して仮想仕事の原理に導入する方法 36)など，いくつかの. 方法が提案されている 37)．Lagrange 未定乗数法は拘束条件を厳密に満足させることができ，高い解析精 度が得られる．しかし拘束条件の数だけ新たな未知数(自由度)が増えるため，それに伴って剛性マトリ クスのサイズが大きくなる．また増える未知数の数はステップごとで異なるため，求解時の取り扱いが 複雑となる．さらに，通常のバンド形とは異なり剛性マトリクスが縁取り形となるため，効率的に解く ためには求解時に工夫を要する．一方 Penalty 法は若干の工具への潜り込みを許容しながらも，剛性マ トリクスのサイズが変わることなく，簡便に扱うことができる手法である．しかし適切な penalty 数を どのように選択するかが大きな問題であり，penalty 数の大きさによって解は異なってくる．拡大 Lagrange 未定乗数法はこの両者を組み合わせた手法である． Penalty 法において解が penalty 数に大きく依存する原因は，剛性マトリクスに非常に大きな負の値が 入る可能性があるためである．これは接触力の法線方向成分を次式のように定義する変形体同士の接触 問題でよく発生する問題であり，その結果剛性マトリクスが特異になりやすい．. (. ). (. tˆ&N = α v A − v B ⋅ n A = α vˆNA − vˆNB. ). ― 4-19. (4.50) ―.
(20) 接触問題の定式化 ただし右肩の添え字は材料の種類を表す．一方本研究で取り扱う剛体工具との unilateral contact を仮定 した接触問題では. tˆ&N = α v A ⋅ n A = α vˆNA. (4.51). と定義するため，この定式に起因して負の penalty 数が入ることはない．その結果剛性マトリクスの対 角項に penalty 数が加えられるだけとなり，”Wiping rows and column method”38)とほぼ一致する形となる． したがって本研究の範疇では penalty 数をどれだけ大きくとっても問題はなく，またその結果許容され る潜り込み量はその逆数程度のオーダーとなり，数値計算上はほぼ厳密に接触の拘束条件を満たすこと ができる．以上のような考察から，本研究では Penalty 法を採用する． 本節では，式(4.10)右辺第一項に示す Penalty 法による接触の拘束条件の，剛性方程式への具体的な導 入方法を示す．それに先立ち，ITAS3D では工具と接触した材料節点に関する自由度は接触位置で設定 される局所座標系を参照して記述されることから，剛性マトリクスへの局所座標系の導入方法を示す． 今，ある材料節点 n が工具と接触しているとする．すでに接触点における局所座標系 n êi は求められ，. 全体座標系 ei との間に次のような座標変換マトリクス [T ] が得られているとする．. ∆u1n e1 ⋅ n eˆ1 e1 ⋅ n eˆ2 n n n ∆u2 = e2 ⋅ eˆ1 e2 ⋅ eˆ2 ∆u n e ⋅ n eˆ e ⋅ n eˆ 3 3 1 3 2. ∆uˆ1n e1 ⋅ n eˆ3 ∆uˆ1n e2 ⋅ n eˆ3 ∆uˆ2n ≡ [T ] ∆uˆ2n n e3 ⋅ n eˆ3 ∆uˆ3n ∆uˆ3 . (4.52). ただし， ∆ui は全体座標系 ei を参照した成分， ∆uˆi は局所座標系 n êi を参照した成分である．同様に接触 力と等価な節点力増分ベクトルにも適用できて，. ∆fˆ1n ∆f1n n n ∆f 2 = [T ] ∆fˆ2 ∆f n ˆn 3 ∆f 3 . (4.53). 一方 {∆ui } と {∆f i } はサブ剛性マトリクスを用いて次式のように関係づけられる．. K ijnn {∆uin } = {∆f jn }. (4.54). 式(4.52)，(4.53)を式(4.54)に代入すると，. [T ]. T. { }. K ijnn [T ]{∆uˆin } = ∆fˆjn. →. ([T ]. T. ). { }. K ijnn [T ] {∆uˆin } = ∆fˆjn. (4.55). この変換により，式(4.54)全体を局所座標系で表記することができる．同様に異なる節点間，例えば節 点 n の変位増分と節点 m の等価節点力増分の関係は，. ( K. nm ij. ). [T ] {∆uˆin } = {∆f jm }. (4.56). また逆に，. ([T ]. T. ). { }. K ijmn {∆uim } = ∆fˆjn. (4.57). 以上の変換を行うことにより，節点 n に関する全ての成分を局所座標系で記述することができる．この ような座標変換を全ての接触節点について行うことにより，接触節点に関する成分全てをそれぞれの局 ― 4-20. ―.
(21) 接触問題の定式化 所座標系で表すことができる． 以上のような座標変換は必ずしも必要ではなく，全体座標系のままでも拘束条件を導入することはで きる．しかし上記の変換を行うことにより，変位拘束を局所座標系の第三成分に与えるだけとなり，拘 束条件の導入が容易となる． 続いて，座標変換された剛性方程式への penalty 法の導入方法. 38). を示す．全節点の変位増分ベクトル. 成分を一列に並べてベクトル表示したときに，i 番目の自由度 ∆ui に. ∆ui = ∆ui. (4.58). という変位拘束を行う場合を考える．以下にその手順を示す． (a) 剛性方程式において既知変位増分の成分を右辺へ移項する．また同時に無視できるほど十分小さい 未知数 X を変位増分ベクトル内で ∆ui の位置に導入する．このとき X は許容される工具内への潜り込み 量を意味し，. X =−. ∆fi. (4.59). α. を満たす．ここでαは剛性マトリクス成分に比べて十分大きい penalty 数である．以上の操作を行うこ とにより，剛性方程式は次式のように変形できる．. K11 M Ki1 M K n1. L K1i L K1n ∆u1 ∆f1 K1i ∆ui K1i X ∆f1 K1i ∆ui O M O M M M M M M M L Kii L Kin X = ∆f i − Kii ∆ui + Kii X ≈ ∆f i − Kii ∆ui O M O M M M M M M M L K ni L K nn ∆un ∆f n K ni ∆ui K ni X ∆f n K ni ∆ui . (4.60). (b) 続いて i 成分対角項に penalty 数αを加える． このとき，Kii + α ≈ α という近似が成り立つことから， 式(4.60)は. K11 M Ki1 M K n1. L K1i L K1n ∆u1 ∆f1 K1i ∆ui 0 ∆f1 K1i ∆ui O M O M M M M M M M L α L Kin X = ∆f i − Kii ∆ui + α X = 0 − Kii ∆ui O M O M M M M M M M L K ni L K nn ∆un ∆f n K ni ∆ui 0 ∆f n K ni ∆ui . (4.61). ただし式(4.59)の関係を用いている．式(4.61)を解くことにより拘束条件式(4.58)を満足した変位増分ベク トルが得られる．またこのときの接触力増分の工具面法線方向成分は式(4.59)から容易に. ∆f i = −α X. (4.62). で与えられる． 以上のように，結果的には既知変数を右辺に移動するとともに，剛性マトリクスの対角項に penalty 数を加えるだけで，拘束条件を導入することができる．また式(4.59)の関係を導入して剛性マトリクス を変形することにより，節点力増分の成分も式(4.62)から容易に求めることができ，求解後の処理も非 常に簡単となる．. ― 4-21. ―.
(22) 接触問題の定式化. 4.7 接触探索アルゴリズム 4.7.1. 接触探索アルゴリズムとは. 工具の移動および材料の変形に伴い，工具と材料の接触・離脱の状況は時々刻々と変化する．4.4 節 では特に接触節点が工具面上を滑る場合に必要となる接触力補正アルゴリズムについて，また 4.5 節で は接触した節点に対する摩擦条件の導入方法について，そして 4.6 節では接触に起因する変位境界条件 の剛性マトリクスへの導入方法についてそれぞれ示した．これらの定式の中で前提となっているのは， 「接触・非接触の判定が正確に行われている」， 「接触節点の場合，接触位置における工具面法線ベクト ルの定義が精度よく行われている」といった点である．このような接触状況の変化の様子をシミュレー ション中で再現していくためには，どの材料節点がどの工具と接触しているのか，またどの材料節点が 工具から離脱しようとしているのかを正確に，かつ効率よく判定していくことが重要となる．また，接 触している材料節点が工具面上を滑ることに伴う接触状態の変化にしっかりと追従していくことも，非 常に重要である．そしてこのような接触状態の変化に関する探索，特に自由節点が接触状態へ移る場合 および接触節 点が工具面上を滑る場合 の探索を行うのが，接 触探索アルゴリズム (contact search algorithm)の主目的である．つまり，例えば THF の FEM 解析を考えた場合には，液圧や軸押し込みの負 荷に伴ってどの材料節点が工具のどの部分と接触したのか，また工具面上を滑ることで接触による拘束 方向どのように変わるのか，といった接触状況の変化を追跡する上で，接触探索アルゴリズムの探索精 度が解析結果に大きく影響することは明らかである．また形状が複雑化するほどこの問題は顕著となる． そこで以上を踏まえて本章では，本研究で採用する接触探索アルゴリズムについてその概要を示す． 4.2 節で簡単に触れたように接触探索アルゴリズムは工具面形状の表現方法に大きく依存するが，ここでは 本研究で採用している工具面が有限要素メッシュで表現されている場合を対象とする．なお，工具から の離脱の判定方法については 4.8 節で示すこととする．. 4.7.2. 接触探索アルゴリズムの概要. 工具と接触した材料節点は，滑り摩擦状態にある限り工具面上を滑る(より正確には，工具面の接線方 向へ動く)ことはできるが，接触位置における工具面の法線方向には動くことはできない．接触探索アル ゴリズムの役割は，①材料節点が実際に工具に接触しているかどうかを判定し，②接触状態にある(にな った)材料節点の変位拘束方向を特定する(工具面法線ベクトルを算出する)ことにある．接触による拘束 条件の具体的な剛性方程式へ導入方法は 4.6 節で示した通りである． 剛体工具を仮定した塑性加工シミュレーションでは，通常材料側の節点が工具側の要素と接触してい るかどうかを探索する．このような手法は master-slave algorithm とよばれ，材料節点を slave node，工具 要素を master segment とよぶ．材料モデルとしてソリッド要素を用いる場合，材料表面を構成する節点 が slave node となる．一方シェル要素を用いる場合，板厚を考慮した接触探索用の節点 X nc を次式のよう に新たに定義し，これを slave node とする 39)-41)．. X nc = X n ±. hn Vn 2. (4.63). ここで X n は節点本来の位置ベクトル， Vn はその節点におけるファイバーベクトル， hn はその節点にお けるファイバーベクトルの長さ(板厚)である．＋の場合がシェル上面，すなわち上型との接触探索用， −の場合がシェル下面，すなわち下型との接触探索用節点である． 通常厳密に工具要素上に材料節点が接する瞬間を求めるのではなく，材料節点が工具要素の法線方向 ― 4-22. ―.
(23) 接触問題の定式化. Contact pair. Punch. Sheet node Fig.4.9 Contact pair of a sheet node and a tool element. に設けられた許容値内に含まれた時点で接触したと見なす場合が多い．接触探索の流れは，以下のよう である． (1). Fig.4.9 に示すように，材料節点とペアになる工具要素を見つける．この材料節点-工具要素の組を 接触ペアとよぶ．. (2). 材料節点が，接触ペアとなった工具要素に設けられた許容値の範囲に含まれるかチェックをし，含 まれている場合はその材料節点はその工具要素と接触すると判定する．. (3). 工具要素上の接触位置の位置ベクトル，およびその点での工具法線ベクトルを算出する．. 接触ペアとは，実際に接触する，しないに関わらず後述するローカルサーチアルゴリズムにしたがっ て決められる材料節点と工具要素の対である．詳細は 4.7.4 節にて示す．接触ペアの探索を行うにあた り，全ての材料節点に関して，接触ペアが見つかるまで全ての工具要素を対象として一つずつ探索を行 うことも可能である．しかしこの手法は非常に計算コストがかかり 42)，実部品レベルの大規模な問題の 解析を行う場合，連立方程式を解く時間よりも要する場合もある．そこで通常接触探索は，グローバル サーチ(Global search)とローカルサーチ(Local search)の二つのステップに分けて行われる． グローバルサーチの役割は，各材料節点に対して接触ペアとなりそうな候補工具要素を限定すること であり，計算時間の短縮に主眼をおいている．したがってこれは，先に示した接触探索の流れの(1)の前 に行う作業である．そしてグローバルサーチで選び出された候補のみに対して(1)の作業をローカルサー チで行う．そしてローカルサーチの判定結果に基づき，(1)で選ばれた要素のみに対して(2),(3)の作業を ローカルサーチの後処理として行う．このようにローカルサーチでは接触しているかどうかの最終的な 判定を行うため，厳密な探索が必要となり計算コストが比較的かかる．つまり接触探索に必要な計算時 間はローカルサーチにかかる時間・回数で決まるため，グローバルサーチにおいて事前絞り込みを行う ことで効率のよい(探索時間が短い)，かつ高精度な接触探索を実現しようとするのが狙いである． これまでさまざまな接触探索アルゴリズムが提案されている．特にグローバルサーチについては効率 のよいアルゴリズムが多数提案されている 42)-48)．そこで本研究でもこれまでに提案された既存のグロー バルサーチアルゴリズムを採用している．一方ローカルサーチアルゴリズムについてはグローバルサー チに比べると研究事例は少なく 49)-54)，またそのどれも一長一短があるのが現状である．そこで本研究で ― 4-23. ―.
関連したドキュメント
In Proceedings Fourth International Conference on Inverse Problems in Engineering (Rio de Janeiro, 2002), H. Orlande, Ed., vol. An explicit finite difference method and a new
Chu, “H ∞ filtering for singular systems with time-varying delay,” International Journal of Robust and Nonlinear Control, vol. Gan, “H ∞ filtering for continuous-time
VARIATIONAL AND NUMERICAL ANALYSIS OF THE SIGNORINI’S CONTACT PROBLEM IN VISCOPLASTICITY WITH DAMAGE
Sofonea, Variational and numerical analysis of a quasistatic viscoelastic problem with normal compliance, friction and damage,
These results are motivated by the bounds for real subspaces recently found by Bachoc, Bannai, Coulangeon and Nebe, and the bounds generalize those of Delsarte, Goethals and Seidel
In this section we apply approximate solutions to obtain existence results for weak solutions of the initial-boundary value problem for Navier-Stokes- type
Let F be a simple smooth closed curve and denote its exterior by Aco.. From here our plan is to approximate the solution of the problem P using the finite element method. The
Let F be a simple smooth closed curve and denote its exterior by Aco.. From here our plan is to approximate the solution of the problem P using the finite element method. The
For the time step Δt 0.05 seconds, the decays of the total potential energy and the angular momentum, shown in Figures 11a and 11b, respectively, are practically the same for
